{- ORMOLU_DISABLE -} -- Implicit CAD. Copyright (C) 2011, Christopher Olah (chris@colah.ca) -- Copyright (C) 2016, Julia Longtin (julial@turinglace.com) -- Released under the GNU AGPLV3+, see LICENSE module Graphics.Implicit.Export.Render.Interpolate (interpolate) where import Prelude((*), (>), (<), (/=), (+), (-), (/), (==), (&&), abs) import Graphics.Implicit.Definitions (ℝ, Fastℕ, ℝ2) import Linear (V2(V2)) default (Fastℕ, ℝ) -- Consider a function f(x): {- | \ f(x) | - \ |_______\________ x | \ -} -- The purpose of interpolate is to find the value of x where f(x) crosses zero. -- This should be accomplished cheaply and accurately. -- We are given the constraint that x will be between a and b. -- We are also given the values of f at a and b: aval and bval. -- Additionaly, we get f (continuous and differentiable almost everywhere), -- and the resolution of the object (so that we can make decisions about -- how precise we need to be). -- While the output will never be used, interpolate will be called -- in cases where f(x) doesn't cross zero (ie. aval and bval are both -- positive or negative. -- Clarification: If f(x) crosses zero, but doesn't necessarily have -- to do so by intermediate value theorem, it is beyond the scope of this -- function. -- If it doesn't cross zero, we don't actually care what answer we give, -- just that it's cheap. -- FIXME: accept resolution on multiple axises. interpolate :: ℝ2 -> ℝ2 -> (ℝ -> ℝ) -> ℝ -> ℝ interpolate (V2 a aval) (V2 _ bval) _ _ | aval*bval > 0 = a -- The obvious: interpolate (V2 a 0) _ _ _ = a interpolate _ (V2 b 0) _ _ = b -- It may seem, at first, that our task is trivial. -- Just use linear interpolation! -- Unfortunately, there's a nasty failure case {- / / ________#________/____ ________________/ -} -- This is really common for us, for example in cubes, -- where another variable dominates. -- It may even be the case that, because we are so close -- to the side, it looks like we are really close to an -- answer... And we just give it back. -- So we need to detect this. And get free goodies while we're -- at it (shrink domain to guess within fromm (a,b) to (a',b')) -- :) interpolate (V2 a aval) (V2 b bval) f _ = -- Make sure aval > bval, then pass to interpolateLin if aval > bval then interpolateLin 0 (V2 a aval) (V2 b bval) f else interpolateLin 0 (V2 b bval) (V2 a aval) f -- Yay, linear interpolation! -- Try the answer linear interpolation gives us... -- (n is to cut us off if recursion goes too deep) interpolateLin :: Fastℕ -> ℝ2 -> ℝ2 -> (ℝ -> ℝ) -> ℝ interpolateLin n (V2 a aval) (V2 b bval) obj | aval /= bval= let -- Interpolate and evaluate mid :: ℝ mid = a + (b-a)*aval/(aval-bval) midval = obj mid -- Are we done? in if midval == 0 then mid -- else let (a', a'val, b', b'val, improveRatio) = if midval > 0 then (mid, midval, b, bval, midval/aval) else (a, aval, mid, midval, midval/bval) -- some times linear interpolate doesn't work, -- because one side is very close to zero and flat -- we catch it because the interval won't shrink when -- this is the case. To test this, we look at whether -- the replaced point evaluates to substantially closer -- to zero than the previous one. in if improveRatio < 0.3 && n < 4 -- And we continue on. then interpolateLin (n+1) (V2 a' a'val) (V2 b' b'val) obj -- But if not, we switch to binary interpolate, which is -- immune to this problem else interpolateBin (n+1) (V2 a' a'val) (V2 b' b'val) obj -- And a fallback: interpolateLin _ (V2 a _) _ _ = a -- Now for binary searching! interpolateBin :: Fastℕ -> ℝ2 -> ℝ2 -> (ℝ -> ℝ) -> ℝ -- The termination case: interpolateBin 5 (V2 a aval) (V2 b bval) _ = if abs aval < abs bval then a else b -- Otherwise, have fun with mid! interpolateBin n (V2 a aval) (V2 b bval) f = let mid :: ℝ mid = (a+b)/2 midval = f mid in if midval > 0 then interpolateBin (n+1) (V2 mid midval) (V2 b bval) f else interpolateBin (n+1) (V2 a aval) (V2 mid midval) f