Merge pull request #440 from lepsa/closing-loops

Adding checks when looking for loops to handle NaNs.

CHANGELOG.md

@@ -7,6 +7,7 @@
     * `torus(r1, r2)`
     * `ellipsoid(a, b, c)`
   * Adding vector-matrix, matrix-vector, and matrix-matrix multiplication support to `*` [#414](https://github.com/Haskell-Things/ImplicitCAD/issues/414)
+  * Several mathematical functions that have undefined values for some inputs or tend to infinity have been given finite but extremely large bounds.
 
 * Haskell interface changes
   * Added matching primitives for `cone`, `torus`, and `ellipsoid`
@@ -16,6 +17,7 @@
   * Migrating StateC and StateE to a ReaderT/WriterT/StateT transformer stack, rather than being just StateT. [#432](https://github.com/Haskell-Things/ImplicitCAD/pull/432)
   * Fixing an off by one error in variable stack lookups. [#431](https://github.com/Haskell-Things/ImplicitCAD/issues/431)
   * Fixing exponent operator precedence. [#428](https://github.com/Haskell-Things/ImplicitCAD/issues/428)
+  * Fixing some crashes relating to unclosed loops when generating surfaces. [#373](https://github.com/Haskell-Things/ImplicitCAD/issues/373)
 
 # Version [0.4.0.0](https://github.com/Haskell-Things/ImplicitCAD/compare/v0.3.0.0...v0.4.0.0) (2022-06-06)
 

Graphics/Implicit/Definitions.hs

@@ -10,6 +10,9 @@
 {-# LANGUAGE UndecidableInstances #-}
 {-# LANGUAGE FlexibleContexts #-}
 {-# LANGUAGE DeriveGeneric #-}
+{-# LANGUAGE TemplateHaskell #-}
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
 
 -- Definitions of the types used when modeling, and a few operators.
 module Graphics.Implicit.Definitions (
@@ -18,6 +21,7 @@ module Graphics.Implicit.Definitions (
     ℝ,
     ℝ2,
     ℝ3,
+    ℝ3' (ℝ3'),
     minℝ,
     ComponentWiseMultable,
     (⋯*),
@@ -71,7 +75,7 @@ where
 
 import GHC.Generics (Generic)
 
-import Prelude (Ord, Eq, atan2, asin, pi, (>=), signum, abs, (+), (-), RealFloat, (==), ($), flip, Semigroup((<>)), Monoid (mempty), Double, Either(Left, Right), Bool(True, False), (*), (/), fromIntegral, Float, realToFrac)
+import Prelude (Ord, Eq, atan2, asin, pi, (>=), signum, abs, (+), (-), RealFloat, (==), ($), flip, Semigroup((<>)), Monoid (mempty), Double, Either(Left, Right), Bool(True, False), (*), (/), fromIntegral, Float, realToFrac, (&&), isNaN, (||))
 
 import Graphics.Implicit.FastIntUtil as F (Fastℕ(Fastℕ), fromFastℕ, toFastℕ)
 
@@ -87,6 +91,7 @@ import Control.Applicative (Applicative(liftA2))
 
 import Text.Show.Combinators
     ( Show(showsPrec, show), (@|), showApp, showCon, PrecShowS)
+import Control.Lens (makeWrapped)
 
 -- | A type synonym for 'Double'. When used in the context of positions or
 -- sizes, measured in units of millimeters. When used as in the context of
@@ -102,6 +107,19 @@ type ℝ2 = V2 ℝ
 -- measured in radians.
 type ℝ3 = V3 ℝ
 
+-- ℝ3 except that we also check if values are NaN because those aren't
+-- equal under the normal floating point equivalence.
+newtype ℝ3' = ℝ3' (V3 ℝ)
+$(makeWrapped ''ℝ3')
+instance Eq ℝ3' where
+  ℝ3' a == ℝ3' b = eqNaNs a b
+
+eqNaNs :: ℝ3 -> ℝ3 -> Bool
+eqNaNs (V3 a b c) (V3 a' b' c') =
+    eqNaN a a' && eqNaN b b' && eqNaN c c'
+eqNaN :: RealFloat a => a -> a -> Bool
+eqNaN a b = (isNaN a && isNaN b) || (a == b)
+
 -- | A give up point for dividing ℝs, and for the maximum difference between abs(n) and abs(-n).
 minℝ :: ℝ
 -- for Doubles.
@@ -221,7 +239,7 @@ data SharedObj obj f a
   | Mirror (f a) obj -- ^ Mirror across the line whose normal is defined by the vector
   | Outset ℝ obj
   | Shell ℝ obj
-  | EmbedBoxedObj ((f a) -> a, ((f a), (f a)))
+  | EmbedBoxedObj (f a -> a, (f a, f a))
   | WithRounding ℝ obj
   deriving (Generic)
 

Graphics/Implicit/Export/Render.hs

@@ -9,9 +9,9 @@
 -- export getContour and getMesh, which returns the edge of a 2D object, or the surface of a 3D object, respectively.
 module Graphics.Implicit.Export.Render (getMesh, getContour) where
 
-import Prelude(error, (-), ceiling, ($), (+), (*), max, div, tail, fmap, reverse, (.), foldMap, min, Int, (<>), (<$>))
+import Prelude(error, (-), ceiling, ($), (+), (*), max, div, tail, fmap, reverse, (.), foldMap, min, Int, (<>), (<$>), traverse)
 
-import Graphics.Implicit.Definitions (ℝ, ℕ, Fastℕ, ℝ2, ℝ3, TriangleMesh, Obj2, SymbolicObj2, Obj3, SymbolicObj3, Polyline(getSegments), (⋯/), fromℕtoℝ, fromℕ)
+import Graphics.Implicit.Definitions (ℝ, ℕ, Fastℕ, ℝ2, ℝ3, TriangleMesh, Obj2, SymbolicObj2, Obj3, SymbolicObj3, Polyline(getSegments), (⋯/), fromℕtoℝ, fromℕ, ℝ3' (ℝ3'))
 
 import Graphics.Implicit.Export.Symbolic.Rebound2 (rebound2)
 
@@ -71,6 +71,7 @@ import Control.Parallel.Strategies (using, rdeepseq, parBuffer)
 import Graphics.Implicit.Export.Render.HandlePolylines (cleanLoopsFromSegs)
 import Data.Maybe (fromMaybe)
 import Graphics.Implicit.Primitives (getImplicit)
+import Control.Lens (_Wrapped, view, over, _Just)
 
 -- Set the default types for the numbers in this file.
 default (ℕ, Fastℕ, ℝ)
@@ -109,15 +110,15 @@ getMesh res@(V3 xres yres zres) symObj =
         par3DList lenx leny lenz =
             [[[ sample mx my mz
             | mx <- [0..lenx] ] | my <- [0..leny] ] | mz <- [0..lenz] ]
-              `using` parBuffer (max 1 $ div (fromℕ $ (lenx+leny+lenz)) forcesteps) rdeepseq
+              `using` parBuffer (max 1 $ div (fromℕ (lenx+leny+lenz)) forcesteps) rdeepseq
 
         -- sample our object(s) at the given point.
         sample :: ℕ -> ℕ -> ℕ -> ℝ
         sample mx my mz = obj $
               V3
-                (x1 + rx*(fromℕtoℝ mx))
-                (y1 + ry*(fromℕtoℝ my))
-                (z1 + rz*(fromℕtoℝ mz))
+                (x1 + rx*fromℕtoℝ mx)
+                (y1 + ry*fromℕtoℝ my)
+                (z1 + rz*fromℕtoℝ mz)
 
         -- (1) Calculate mid points on X, Y, and Z axis in 3D space.
         midsZ = [[[
@@ -168,7 +169,13 @@ getMesh res@(V3 xres yres zres) symObj =
         minres = xres `min` yres `min` zres
         sqTris = [[[
             foldMap (tesselateLoop minres obj) $
-              fromMaybe (error "unclosed loop in paths given") $ getLoops $
+              fromMaybe (error "unclosed loop in paths given") $
+              -- Shove the ℝ3s into ℝ3's to get the NaN checks, then
+              -- unwrap everything. This should mostly compile away
+              -- given that it is lensy and passing a newtype instance
+              -- around. `getLoops` is the function actually doing the
+              -- work we care about
+              over (_Just . traverse . traverse . traverse) (view _Wrapped) . getLoops . over (traverse . traverse) ℝ3' $
                         segX''' <>
                    mapR segX''T <>
                    mapR segY''' <>
@@ -263,19 +270,19 @@ getContour res@(V2 xres yres) symObj =
 -- utility functions
 
 injX :: ℝ -> Polyline -> [ℝ3]
-injX val polyline = (prepend val) <$> getSegments polyline
+injX val polyline = prepend val <$> getSegments polyline
   where
     prepend :: ℝ -> ℝ2 -> ℝ3
     prepend a (V2 b c) = V3 a b c
 
 injY :: ℝ -> Polyline -> [ℝ3]
-injY val polyline = (insert val)  <$> getSegments polyline
+injY val polyline = insert val  <$> getSegments polyline
   where
     insert :: ℝ -> ℝ2 -> ℝ3
     insert b (V2 a c) = V3 a b c
 
 injZ :: ℝ -> Polyline -> [ℝ3]
-injZ val polyline = (postfix val) <$> getSegments polyline
+injZ val polyline = postfix val <$> getSegments polyline
   where
     postfix :: ℝ -> ℝ2 -> ℝ3
     postfix c (V2 a b) = V3 a b c

Graphics/Implicit/ExtOpenScad/Default.hs

@@ -13,7 +13,7 @@
 module Graphics.Implicit.ExtOpenScad.Default (defaultObjects) where
 
 -- be explicit about where we pull things in from.
-import Prelude (Bool(True, False), Maybe(Just, Nothing), ($), (<>), (<$>), fmap, pi, sin, cos, tan, asin, acos, atan, sinh, cosh, tanh, abs, signum, fromInteger, (.), floor, ceiling, round, exp, log, sqrt, max, min, atan2, (**), flip, (<), (>), (<=), (>=), (==), (/=), (&&), (||), not, show, foldl, (*), (/), mod, (+), zipWith, (-), otherwise, id, foldMap, fromIntegral, IO, pure, Int)
+import Prelude (Bool(True, False), Maybe(Just, Nothing), ($), (<>), (<$>), fmap, pi, sin, cos, tan, asin, acos, atan, sinh, cosh, tanh, abs, signum, fromInteger, (.), floor, ceiling, round, exp, log, sqrt, max, min, atan2, (**), flip, (<), (>), (<=), (>=), (==), (/=), (&&), (||), not, show, foldl, (*), (/), mod, (+), zipWith, (-), otherwise, id, foldMap, fromIntegral, IO, pure, Int, isNaN, negate, RealFloat, Ord)
 import qualified Prelude as P (length)
 
 import Graphics.Implicit.Definitions (ℝ, ℕ)
@@ -42,6 +42,10 @@ import System.Random (randomRIO)
 import Data.Maybe (maybe)
 import Data.Tuple (snd)
 import Linear.Matrix ((!*!), (*!), (!*))
+import Graphics.Implicit.MathUtil (infty)
+
+clamp :: Ord a => (a, a) -> a -> a
+clamp (lower, upper) a = min upper (max lower a)
 
 defaultObjects :: Bool -> VarLookup
 defaultObjects withCSG = VarLookup $ fromList $
@@ -58,14 +62,27 @@ defaultConstants = (\(a,b) -> (a, toOObj (b :: ℝ))) <$>
     [(Symbol "pi", pi),
      (Symbol "PI", pi)]
 
+-- Values and functions for dealing with NaNs and Infinities.
+minimumValue :: ℝ
+minimumValue = -1e100
+maximumValue :: ℝ
+maximumValue = 1e100
+nanNegInf :: RealFloat a => a -> a
+nanNegInf x = if isNaN x then -infty else x
+signedNaNInf :: RealFloat a => a -> a -> a
+signedNaNInf x y = if isNaN y then signum x * infty else y
+
 defaultFunctions :: [(Symbol, OVal)]
 defaultFunctions = (\(a,b) -> (a, toOObj ( b :: ℝ -> ℝ))) <$>
     [
         (Symbol "sin",   sin),
         (Symbol "cos",   cos),
         (Symbol "tan",   tan),
-        (Symbol "asin",  asin),
-        (Symbol "acos",  acos),
+        -- If the value is NaN, set it to the signed infinity of the input
+        -- and then clamp the values so that infinity doesn't propagate.
+        (Symbol "asin",  \x -> clamp (minimumValue, maximumValue) . signedNaNInf x $ asin x),
+        -- same as asin, but we need to invert the input sign when clamping
+        (Symbol "acos",  \x -> clamp (minimumValue, maximumValue) . signedNaNInf (negate x) $ acos x),
         (Symbol "atan",  atan),
         (Symbol "sinh",  sinh),
         (Symbol "cosh",  cosh),
@@ -76,10 +93,14 @@ defaultFunctions = (\(a,b) -> (a, toOObj ( b :: ℝ -> ℝ))) <$>
         (Symbol "ceil",  fromInteger . ceiling ),
         (Symbol "round", fromInteger . round ),
         (Symbol "exp",   exp),
-        (Symbol "ln",    log),
-        (Symbol "log",   log),
+        -- Log is undefined for negative values, so we are taking those NaNs
+        -- and -Infinity values and clamping them to a very negative, but
+        -- finite, value.
+        (Symbol "ln",    clamp (minimumValue, infty) . nanNegInf . log),
+        (Symbol "log",   clamp (minimumValue, infty) . nanNegInf . log),
         (Symbol "sign",  signum),
-        (Symbol "sqrt",  sqrt)
+        -- same as log, but clamping to 0 rather than a very large negative value
+        (Symbol "sqrt",  clamp (0, infty) . nanNegInf . sqrt)
     ]
 
 defaultFunctions2 :: [(Symbol, OVal)]
@@ -325,12 +346,12 @@ defaultPolymorphicFunctions =
 
         divide = OFunc $ \case
             (ONum a) -> OFunc $ \case
-                (ONum b) -> ONum (a/b)
+                (ONum b) -> ONum (clamp (minimumValue, maximumValue) $ a/b)
                 b        -> errorAsAppropriate "divide" (ONum a) b
             a -> OFunc $ \case
                 b -> div' a b
 
-        div' (ONum a)  (ONum b) = ONum  (a/b)
+        div' (ONum a)  (ONum b) = ONum  (clamp (minimumValue, maximumValue) $  a/b)
         div' (OList a) (ONum b) = OList (fmap (\x -> div' x (ONum b)) a)
         div' a         b        = errorAsAppropriate "divide" a b
 

tests/GoldenSpec/Spec.hs

@@ -170,3 +170,27 @@ spec = describe "golden tests" $ do
       , ellipsoid 10 15 20
       , translate (V3 0 0 25) $ cone 20 20
       ]
+  golden "closing-paths-1" 0.5 $
+    extrudeM
+      (Left 0)
+      (C1 1)
+      (Left 0)
+      (circle 1)
+      -- limit height or this races off to infinity
+      -- and makes the tests take forever, eating all
+      -- your RAM while it tries to calculate the surface
+      $ Right $ \(V2 x y) -> min 100 $ 1 / sqrt (x ^ 2 + y ^ 2)
+  golden "closing-paths-2" 1 $
+    extrudeM
+      -- Note, this have a gap from the base plane and the extruded
+      -- object, and the base of the extruded object is going to be
+      -- missing. This is because we are directly using haskell
+      -- functions rather than the clamped / Infinity / NaN checked
+      -- versions that the parser is set up to provide.
+      -- However, this is still useful in that it doesn't crash on
+      -- unclosed loops.
+      (pure $ \h -> 35 * log (h*2*pi/30))
+      (C1 1)
+      (Left 0)
+      (union [circle 10])
+      $ Left 40