From 17e2f7825dec84c34d5cb891dfef31c2241918b8 Mon Sep 17 00:00:00 2001
From: Julia Longtin <julia.longtin@gmail.com>
Date: Mon, 23 Feb 2026 16:17:38 +0000
Subject: [PATCH 01/28] add polyhedron support.

---
 Graphics/Implicit/Canon.hs                   |   3 +
 Graphics/Implicit/Definitions.hs             |   3 +
 Graphics/Implicit/Export/SymbolicFormats.hs  |  22 ++-
 Graphics/Implicit/Export/TextBuilderUtils.hs |   5 +
 Graphics/Implicit/ExtOpenScad/Primitives.hs  |  31 +++-
 Graphics/Implicit/ObjectUtil/GetBox3.hs      |  22 ++-
 Graphics/Implicit/ObjectUtil/GetImplicit3.hs | 162 +++++++++++++++++--
 Graphics/Implicit/Primitives.hs              |  15 +-
 Makefile                                     |   2 +-
 tests/Graphics/Implicit/Test/Instances.hs    |   9 +-
 10 files changed, 239 insertions(+), 35 deletions(-)

diff --git a/Graphics/Implicit/Canon.hs b/Graphics/Implicit/Canon.hs
index 0cb72d62..1f234145 100644
--- a/Graphics/Implicit/Canon.hs
+++ b/Graphics/Implicit/Canon.hs
@@ -83,6 +83,7 @@ import Graphics.Implicit.Definitions
   , SymbolicObj3
       ( Cube
       , Cylinder
+      , Polyhedron
       , Extrude
       , ExtrudeM
       , ExtrudeOnEdgeOf
@@ -176,6 +177,7 @@ fmapObj3
 fmapObj3 f _ _ (Cube v) = f $ Cube v
 fmapObj3 f _ _ (Sphere r) = f $ Sphere r
 fmapObj3 f _ _ (Cylinder r1 r2 h) = f $ Cylinder r1 r2 h
+fmapObj3 f _ _ (Polyhedron points faces) = f $ Polyhedron points faces
 fmapObj3 f _ _ (Torus r1 r2) = f $ Torus r1 r2
 fmapObj3 f _ _ (Ellipsoid a b c) = f $ Ellipsoid a b c
 fmapObj3 f _ _ (BoxFrame b e) = f $ BoxFrame b e
@@ -239,6 +241,7 @@ instance EqObj SymbolicObj3 where
   Ellipsoid a1 b1 c1 =^= Ellipsoid a2 b2 c2 = a1 == a2 && b1 == b2 && c1 == c2
   Cylinder r1a r2a ha =^= Cylinder r1b r2b hb = r1a == r1b && r2a == r2b && ha == hb
   BoxFrame b1 e1 =^= BoxFrame b2 e2 = b1 == b2 && e1 == e2
+  Polyhedron p1 f1 =^= Polyhedron p2 f2 = p1 == p2 && f1 == f2
   Link a1 b1 c1 =^= Link a2 b2 c2 = a1 == a2 && b1 == b2 && c1 == c2
   Rotate3 x a =^= Rotate3 y b = x == y && a =^= b
   Transform3 x a =^= Transform3 y b = x == y && a =^= b
diff --git a/Graphics/Implicit/Definitions.hs b/Graphics/Implicit/Definitions.hs
index 5153e0dc..4b2282d0 100644
--- a/Graphics/Implicit/Definitions.hs
+++ b/Graphics/Implicit/Definitions.hs
@@ -56,6 +56,7 @@ module Graphics.Implicit.Definitions (
         Cube,
         Sphere,
         Cylinder,
+        Polyhedron,
         Rotate3,
         Transform3,
         Torus,
@@ -331,6 +332,7 @@ data SymbolicObj3 =
       Cube ℝ3 -- rounding, size.
     | Sphere ℝ -- radius
     | Cylinder ℝ ℝ ℝ --
+    | Polyhedron [ℝ3] [(ℕ,ℕ,ℕ)] -- virtexes, triangles-by-index
     -- Simple transforms
     | Rotate3 (Quaternion ℝ) SymbolicObj3
     | Transform3 (M44 ℝ) SymbolicObj3
@@ -364,6 +366,7 @@ instance Show SymbolicObj3 where
     -- centered.
     Cube sz -> showCon "cube" @| False @| sz
     Sphere d -> showCon "sphere" @| d
+    Polyhedron points tris -> showCon "polyhedron" @| points @| tris
     BoxFrame b e -> showCon "boxFrame" @| b @| e
     Link le r1 r2 -> showCon "link" @| le @| r1 @| r2
     -- NB: The arguments to 'Cylinder' are backwards compared to 'cylinder' and
diff --git a/Graphics/Implicit/Export/SymbolicFormats.hs b/Graphics/Implicit/Export/SymbolicFormats.hs
index 6c5018e9..7d09e00c 100644
--- a/Graphics/Implicit/Export/SymbolicFormats.hs
+++ b/Graphics/Implicit/Export/SymbolicFormats.hs
@@ -11,8 +11,9 @@ module Graphics.Implicit.Export.SymbolicFormats (scad2, scad3) where
 
 import Prelude((.), fmap, Either(Left, Right), ($), (*), (-), (/), pi, error, (+), (==), take, floor, (&&), const, (<>), (<$>))
 
-import Graphics.Implicit.Definitions(ℝ, SymbolicObj2(Shared2, Square, Circle, Polygon, Rotate2, Transform2, Slice), SymbolicObj3(Shared3, Cube, Sphere, Cylinder, BoxFrame, Rotate3, Transform3, Extrude, ExtrudeM, RotateExtrude, ExtrudeOnEdgeOf, Torus, Ellipsoid, Link), isScaleID, SharedObj(Empty, Full, Complement, UnionR, IntersectR, DifferenceR, Translate, Scale, Mirror, Outset, Shell, EmbedBoxedObj, WithRounding), quaternionToEuler)
-import Graphics.Implicit.Export.TextBuilderUtils(Text, bf)
+import Graphics.Implicit.Definitions(ℝ, SymbolicObj2(Shared2, Square, Circle, Polygon, Rotate2, Transform2, Slice), SymbolicObj3(Shared3, Cube, Sphere, Cylinder, Polyhedron, BoxFrame, Rotate3, Transform3, Extrude, ExtrudeM, RotateExtrude, ExtrudeOnEdgeOf, Torus, Ellipsoid, Link), isScaleID, SharedObj(Empty, Full, Complement, UnionR, IntersectR, DifferenceR, Translate, Scale, Mirror, Outset, Shell, EmbedBoxedObj, WithRounding), quaternionToEuler)
+
+import Graphics.Implicit.Export.TextBuilderUtils(Text, bf, bℕ)
 
 -- For constructing vectors of ℝs.
 import Linear (V2(V2), V3(V3), V4(V4))
@@ -131,11 +132,16 @@ buildS3 _ (Cube (V3 w d h)) = call "cube" [pretty $ bf w, pretty $ bf d, pretty
 
 buildS3 _ (Sphere r) = callNaked "sphere" ["r = " <> pretty (bf r)] []
 
-buildS3 _ (Cylinder h r1 r2) = callNaked "cylinder" [
-                              "r1 = " <> pretty (bf r1)
-                             ,"r2 = " <> pretty (bf r2)
-                             , pretty $ bf h
-                             ] []
+buildS3 _ (Cylinder h r1 r2) = callNaked "cylinder" ["r1 = " <> pretty (bf r1)
+                                                    ,"r2 = " <> pretty (bf r2)
+                                                    , pretty $ bf h
+                                                    ] []
+
+buildS3 _ (Polyhedron points tris) = callNaked "polyhedron" ["points = [" <> (fold $ intersperse "," $ renderPoint <$> points) <> "] faces = [" <> (fold $ intersperse "," $ renderTri <$> tris) <> "]" ] []
+  where
+    renderPoint (V3 v1 v2 v3) = "[" <> pretty (bf v1) <> "," <> pretty (bf v2) <> "," <> pretty (bf v3) <> "]"
+    renderTri (n1,n2,n3) = "[" <> pretty (bℕ n1) <> "," <> pretty (bℕ n2) <> "," <> pretty (bℕ n3) <> "]"
+
 buildS3 res (Rotate3 q obj) =
   let (V3 x y z) = quaternionToEuler q
    in call "rotate" [pretty $ bf (rad2deg x), pretty $ bf (rad2deg y), pretty $ bf (rad2deg z)] [buildS3 res obj]
@@ -180,7 +186,7 @@ buildS2 res (Shared2 obj) = buildShared res obj
 
 buildS2 _ (Circle r) = call "circle" [pretty $ bf r] []
 
-buildS2 _ (Polygon points) = call "polygon" (fmap bvect points) []
+buildS2 _ (Polygon points) = call "polygon" (bvect <$> points) []
 
 buildS2 res (Rotate2 r obj)     = call "rotate" [pretty $ bf (rad2deg r)] [buildS2 res obj]
 
diff --git a/Graphics/Implicit/Export/TextBuilderUtils.hs b/Graphics/Implicit/Export/TextBuilderUtils.hs
index efc4ac4f..f9daae56 100644
--- a/Graphics/Implicit/Export/TextBuilderUtils.hs
+++ b/Graphics/Implicit/Export/TextBuilderUtils.hs
@@ -13,6 +13,7 @@ module Graphics.Implicit.Export.TextBuilderUtils (
      toLazyText,
      -- some special case Builders.
      bf,
+     bℕ,
      buildTruncFloat,
      buildℕ,
      buildInt
@@ -37,6 +38,10 @@ toLazyText = toLazyTextWith defaultChunkSize
 bf :: ℝ -> Text
 bf value = toLazyText . formatRealFloat Exponent Nothing $ fromℝtoFloat value
 
+-- | Serialize a value as an Integer.
+bℕ :: ℕ -> Text
+bℕ = toLazyText . decimal
+
 -- | Serialize a float with four decimal places
 buildTruncFloat :: ℝ -> Builder
 buildTruncFloat = formatRealFloat Fixed $ Just 4
diff --git a/Graphics/Implicit/ExtOpenScad/Primitives.hs b/Graphics/Implicit/ExtOpenScad/Primitives.hs
index cf0ce1c5..f551a84f 100644
--- a/Graphics/Implicit/ExtOpenScad/Primitives.hs
+++ b/Graphics/Implicit/ExtOpenScad/Primitives.hs
@@ -15,7 +15,7 @@
 -- Export one set containing all of the primitive modules.
 module Graphics.Implicit.ExtOpenScad.Primitives (primitiveModules) where
 
-import Prelude((.), Either(Left, Right), Bool(True, False), Maybe(Just, Nothing), ($), pure, either, id, (-), (==), (&&), (<), (*), cos, sin, pi, (/), (>), const, uncurry, (/=), (||), not, null, fmap, (<>), otherwise, error, (<*>), (<$>))
+import Prelude((.), Either(Left, Right), Bool(True, False), Maybe(Just, Nothing), ($), pure, show, either, id, (-), (==), (&&), (<), (*), cos, sin, pi, (/), (>), const, uncurry, (/=), (||), not, null, fmap, (<>), otherwise, error, (<*>), (<$>))
 
 import Graphics.Implicit.Definitions (ℝ, ℝ2, ℝ3, ℕ, SymbolicObj2, SymbolicObj3, ExtrudeMScale(C1), fromℕtoℝ, isScaleID)
 
@@ -27,10 +27,12 @@ import qualified Graphics.Implicit.ExtOpenScad.Util.ArgParser as GIEUA (argument
 
 import Graphics.Implicit.ExtOpenScad.Util.OVal (OTypeMirror, caseOType, divideObjs, (<||>))
 
-import Graphics.Implicit.ExtOpenScad.Util.StateC (errorC)
+import Graphics.Implicit.ExtOpenScad.Util.StateC (errorC, warnC)
 
 -- Note the use of a qualified import, so we don't have the functions in this file conflict with what we're importing.
-import qualified Graphics.Implicit.Primitives as Prim (withRounding, sphere, rect3, rect, translate, circle, polygon, extrude, cylinder2, union, unionR, intersect, intersectR, difference, differenceR, rotate, slice, transform, rotate3V, rotate3, transform3, scale, extrudeM, rotateExtrude, shell, mirror, pack3, pack2, torus, ellipsoid, cone)
+import qualified Graphics.Implicit.Primitives as Prim (withRounding, sphere, rect3, rect, translate, circle, polygon, polyhedron, extrude, cylinder2, union, unionR, intersect, intersectR, difference, differenceR, rotate, slice, transform, rotate3V, rotate3, transform3, scale, extrudeM, rotateExtrude, shell, mirror, pack3, pack2, torus, ellipsoid, cone)
+
+import Data.List (concatMap)
 
 import Control.Monad (when, mplus)
 
@@ -60,6 +62,7 @@ primitiveModules =
   , onModIze torus [([("r1", noDefault), ("r2", hasDefault)], noSuite)]
   , onModIze ellipsoid [([("a", noDefault), ("b", hasDefault), ("c", hasDefault)], noSuite)]
   , onModIze polygon [([("points", noDefault)], noSuite)]
+  , onModIze polyhedron [([("points", noDefault), ("faces", noDefault)], noSuite)]
   , onModIze union [([("r", hasDefault)], requiredSuite)]
   , onModIze intersect [([("r", hasDefault)], requiredSuite)]
   , onModIze difference [([("r", hasDefault)], requiredSuite)]
@@ -265,6 +268,25 @@ cylinder = moduleWithoutSuite "cylinder" $ \_ _ -> do
         in shift obj3
         else shift $ Prim.cylinder2 r1 r2 dh
 
+polyhedron :: (Symbol, SourcePosition -> [OVal] -> ArgParser (StateC [OVal]))
+polyhedron = moduleWithoutSuite "polyhedron" $ \_ _ -> do
+    example "polyhedron(points=[[0,0,0], [2,0,0], [2,2,0], [0,2,0], [1, 1, 2]], faces=[[0,1,2,3], [0,5,1], [1,5,2], [2,5,3], [3,5,4], [4,5,0]]);"
+    -- arguments
+    points :: [ℝ3] <- argument "points" `defaultTo` [] `doc` "list of points to construct faces from"
+    faces :: [[ℕ]] <- argument "faces" `defaultTo` [] `doc` "list of sets of indices into points, used to create faces on the polyhedron."
+    -- A tri is constructed of three indexes into the points.
+    -- This decomposes our faces into tris.
+    let
+      tris = concatMap trianglesFromFace faces
+      in
+      addObj3 $ Prim.polyhedron points tris
+      where
+        -- FIXME: use warnC here, instead of error.
+        trianglesFromFace :: [ℕ] -> [(ℕ,ℕ,ℕ)]
+        trianglesFromFace [p1,p2,p3] = [(p1,p2,p3)]
+        trianglesFromFace (p1:p2:p3:xs) = ((p1,p2,p3):trianglesFromFace (p2:p3:xs))
+        trianglesFromFace ys = error $ "too few points:" <> (show ys) <> "\n"
+
 cone :: (Symbol, SourcePosition -> [OVal] -> ArgParser (StateC [OVal]))
 cone = moduleWithoutSuite "cone" $ \_ _ -> do
     example "cone(r=10, h=30, center=true);"
@@ -351,8 +373,7 @@ circle = moduleWithoutSuite "circle" $ \_ _ -> do
         else Prim.polygon
             [V2 (r*cos θ) (r*sin θ) | θ <- [2*pi*fromℕtoℝ n/fromℕtoℝ sides | n <- [0 .. sides - 1]]]
 
--- | FIXME: 3D Polygons?
---   FIXME: handle rectangles that are not grid alligned.
+-- | FIXME: handle rectangles that are not grid alligned.
 --   FIXME: allow for rounding of polygon corners, specification of vertex ordering.
 --   FIXME: polygons have to have more than two points, or do not generate geometry, and generate an error.
 polygon :: (Symbol, SourcePosition -> [OVal] -> ArgParser (StateC [OVal]))
diff --git a/Graphics/Implicit/ObjectUtil/GetBox3.hs b/Graphics/Implicit/ObjectUtil/GetBox3.hs
index d2604dad..472902eb 100644
--- a/Graphics/Implicit/ObjectUtil/GetBox3.hs
+++ b/Graphics/Implicit/ObjectUtil/GetBox3.hs
@@ -5,13 +5,19 @@
 
 module Graphics.Implicit.ObjectUtil.GetBox3 (getBox3) where
 
-import Prelude(uncurry, pure, Bool(False), Either (Left, Right), (==), max, (/), (-), (+), fmap, unzip, ($), (<$>), (.), minimum, maximum, min, (>), (*), (<), abs, either, const, otherwise, take, fst, snd)
+import Prelude(foldl, uncurry, pure, Bool(False), Either (Left, Right), (==), max, (/), (-), (+), fmap, unzip, ($), (<$>), (.), minimum, maximum, min, (>), (*), (<), abs, either, const, otherwise, take, fst, snd)
+
+-- For Maybe types.
+import Data.Maybe (fromMaybe, Maybe(Just, Nothing))
+
+import Linear (V2(V2), V3(V3))
+import qualified Linear (rotate, point, normalizePoint, (!*))
 
 import Graphics.Implicit.Definitions
     ( Fastℕ,
       fromFastℕ,
       ExtrudeMScale(C2, C1),
-      SymbolicObj3(Shared3, Cube, Sphere, Cylinder, Rotate3, Transform3, Extrude, ExtrudeOnEdgeOf, ExtrudeM, RotateExtrude, Torus, Ellipsoid, BoxFrame, Link),
+      SymbolicObj3(Shared3, Cube, Sphere, Cylinder, Polyhedron, Rotate3, Transform3, Extrude, ExtrudeOnEdgeOf, ExtrudeM, RotateExtrude, Torus, Ellipsoid, BoxFrame, Link),
       Box3,
       ℝ,
       fromFastℕtoℝ,
@@ -21,9 +27,6 @@ import Graphics.Implicit.ObjectUtil.GetBox2 (getBox2, getBox2R)
 
 import Graphics.Implicit.ObjectUtil.GetBoxShared (corners, pointsBox, getBoxShared)
 
-import Linear (V2(V2), V3(V3))
-import qualified Linear
-
 -- FIXME: many variables are being ignored here. no rounding for intersect, or difference.. etc.
 
 -- Get a Box3 around the given object.
@@ -33,6 +36,15 @@ getBox3 (Shared3 obj) = getBoxShared obj
 getBox3 (Cube size) = (pure 0, size)
 getBox3 (Sphere r) = (pure (-r), pure r)
 getBox3 (Cylinder h r1 r2) = (V3 (-r) (-r) 0, V3 r r h ) where r = max r1 r2
+getBox3 (Polyhedron points _) = (minimum_point, maximum_point)
+  where
+    (minimum_point, maximum_point) = fromMaybe (V3 0 0 0, V3 0 0 0) maybeVs
+    maybeVs :: (Maybe (V3 ℝ,V3 ℝ))
+    maybeVs = foldl findMinMax Nothing points
+      where
+        findMinMax :: (Maybe (V3 ℝ,V3 ℝ)) -> V3 ℝ -> (Maybe (V3 ℝ,V3 ℝ))
+        findMinMax Nothing newV3 = Just (newV3, newV3)
+        findMinMax (Just (V3 minx miny minz,V3 maxx maxy maxz)) (V3 newx newy newz) = Just (V3 (min minx newx) (min miny newy) (min minz newz), V3 (max maxx newx) (max maxy newy) (max maxz newz))
 getBox3 (Torus r1 r2) =
   let r = r1 + r2
   in (V3 (-r) (-r) (-r2), V3 r r r2)
diff --git a/Graphics/Implicit/ObjectUtil/GetImplicit3.hs b/Graphics/Implicit/ObjectUtil/GetImplicit3.hs
index db239303..a69aac42 100644
--- a/Graphics/Implicit/ObjectUtil/GetImplicit3.hs
+++ b/Graphics/Implicit/ObjectUtil/GetImplicit3.hs
@@ -5,38 +5,52 @@
 
 module Graphics.Implicit.ObjectUtil.GetImplicit3 (getImplicit3) where
 
-import Prelude (id, (||), (/=), either, round, fromInteger, Either(Left, Right), abs, (-), (/), (*), sqrt, (+), atan2, max, cos, minimum, ($), sin, pi, (.), Bool(True, False), ceiling, floor, pure, (==), otherwise, (**), min, Num, Applicative)
+-- Import only what we need from the prelude.
+import Prelude (abs, atan2, cos, ceiling, either, error, floor, fromInteger, id, max, min, minimum, negate, odd, otherwise, pi, pure, round, show, sin, sqrt, (||), (/=), Either(Left, Right), (<), (<=), (<>), (>), (>=), (&&), (-), (/), (*), (+), ($), (.), Bool(True, False), (==), (**), Num, Applicative, (<$>))
+
+import Control.Lens ((^.))
+
+import qualified Data.Either as Either (either)
+
+import Data.List(genericIndex, length)
+
+import Linear (V2(V2), V3(V3), _xy, _z, distance, dot)
+import qualified Linear (conjugate, inv44, normalizePoint, normalize, point, rotate, Metric)
+
+import Linear.V3 (cross)
+
+-- Matrix times column vector.
+import Linear.Matrix ((!*))
+
+-- Matrix times scalar.
+import Linear.Vector ((*^))
 
 import Graphics.Implicit.Definitions
     ( objectRounding,
       ObjectContext,
       ℕ,
-      SymbolicObj3(Cube, Sphere, Cylinder, Rotate3, Transform3, Extrude,
+      SymbolicObj3(Cube, Sphere, Cylinder, Polyhedron, Rotate3, Transform3, Extrude,
                    ExtrudeM, ExtrudeOnEdgeOf, RotateExtrude, Shared3, Torus, Ellipsoid, BoxFrame, Link),
       Obj3,
       ℝ2,
       ℝ,
       fromℕtoℝ,
+      fromℕ,
       toScaleFn,
       ℝ3 )
 
-import Graphics.Implicit.MathUtil ( rmax, rmaximum )
+-- For handling extrusion of 2D shapes to 3D.
+import {-# SOURCE #-} Graphics.Implicit.Primitives (getImplicit)
 
-import qualified Data.Either as Either (either)
+import Graphics.Implicit.MathUtil (rmax, rmaximum)
 
--- Use getImplicit for handling extrusion of 2D shapes to 3D.
 import Graphics.Implicit.ObjectUtil.GetImplicitShared (getImplicitShared)
-import Linear (V2(V2), V3(V3), _xy, _z)
-import qualified Linear
-
-import {-# SOURCE #-} Graphics.Implicit.Primitives (getImplicit)
-import Control.Lens ((^.))
 
 default (ℝ)
 
 -- Length similar to the opengl version, needed for some of the shape definitions
 openglLength :: (Linear.Metric f, Num (f ℝ), Applicative f) => f ℝ -> ℝ
-openglLength v = Linear.distance (abs v) $ pure 0
+openglLength v = distance (abs v) $ pure 0
 
 -- Component wise maximum. This is what the opengl language is doing, so we need
 -- it for the function as defined by the blog above.
@@ -59,6 +73,28 @@ getImplicit3 _ (Cylinder h r1 r2) = \(V3 x y z) ->
         θ = atan2 (r2-r1) h
     in
         max (d * cos θ) (abs (z-h/2) - (h/2))
+-- FIXME: Make Polyhedron correct by construction.
+getImplicit3 _ (Polyhedron [] _) = error "Asked to find distance to an empty polygon. No points."
+getImplicit3 _ (Polyhedron _ []) = error "Asked to find distance to an empty polygon. No tris."
+getImplicit3 _ (Polyhedron points tris) = \(point) ->
+  if insideMesh triangles point
+  then          minimum $ distancePointToTriangle point <$> triangles
+  else negate $ minimum $ distancePointToTriangle point <$> triangles
+  where
+    -- decompose our tris into triangles.
+    triangles = findTriangle points <$> tris
+    -- FIXME: make these indices correct by construction?
+    findTriangle :: [V3 ℝ] -> (ℕ,ℕ,ℕ) -> (V3 ℝ, V3 ℝ, V3 ℝ)
+    findTriangle virtices (i1,i2,i3)
+      | outOfRange i1 = error $ "bad virtex index(out of range): " <> show i1 <> "\n"
+      | outOfRange i2 = error $ "bad virtex index(out of range): " <> show i2 <> "\n"
+      | outOfRange i3 = error $ "bad virtex index(out of range): " <> show i3 <> "\n"
+      -- FIXME: there are many more degenerate forms of polyhedron possible here. move polyhedron to only holding a mesh?
+      | otherwise = (genericIndex virtices i1, genericIndex virtices i2, genericIndex virtices i3)
+        where
+          -- FIXME: >=BASE-4.21: replace this with compareLength once debian stable ships 4.21.
+          outOfRange :: ℕ -> Bool
+          outOfRange v = v < 0 || length virtices < fromℕ v
 getImplicit3 _ (BoxFrame b e) = \p' ->
     let p@(V3 px py pz) = abs p' - b
         V3 qx qy qz = abs (p + pure e) - pure e
@@ -77,7 +113,7 @@ getImplicit3 _ (Link le r1 r2) = \(V3 px py pz) ->
 getImplicit3 ctx (Rotate3 q symbObj) =
     getImplicit3 ctx symbObj . Linear.rotate (Linear.conjugate q)
 getImplicit3 ctx (Transform3 m symbObj) =
-    getImplicit3 ctx symbObj . Linear.normalizePoint . (Linear.inv44 m Linear.!*) . Linear.point
+    getImplicit3 ctx symbObj . Linear.normalizePoint . (Linear.inv44 m !*) . Linear.point
 -- 2D Based
 getImplicit3 ctx (Extrude h symbObj) =
     let
@@ -190,3 +226,105 @@ getImplicit3 ctx (RotateExtrude totalRotation translate rotate symbObj) =
                     (obj rz_pos)
               else obj rz_pos
 getImplicit3 ctx (Shared3 obj) = getImplicitShared ctx obj
+
+-- | Check to see if a point is inside or outside of a mesh.
+insideMesh :: [(V3 ℝ, V3 ℝ, V3 ℝ)] -> V3 ℝ -> Bool
+insideMesh triangles point = odd $ length $ foundIntersection point unsafeSafeRay <$> triangles
+  where
+    -- our assumed safe ray. it's not really safe, but.. gotta start somewhere.
+    unsafeSafeRay :: V3 ℝ
+    unsafeSafeRay = Linear.normalize $ V3 (sqrt 2) (sqrt 3) (sqrt 5)
+    foundIntersection :: V3 ℝ -> V3 ℝ -> (V3 ℝ, V3 ℝ, V3 ℝ) -> Bool
+    foundIntersection source direction (v1, v2, v3)
+      -- shouldn't happen, triangle is on the same plane as unsafeSafeRay
+      | abs determinate <= 0 = False
+      -- the intersection is not on the triangle (distance along vec12 out of range)
+      | u <= 0 || u > 1      = False
+      | v <= 0 || v > 1      = False
+      | otherwise            = True
+      where
+        -- Our edge vectors. We have picked v1 to address the space by, for convenience.
+        vec12 = v2 - v1
+        vec13 = v3 - v1
+        -- The normal of the plane formed by vec13 and our 'safe' vector. at 90 degree angles to both.
+        norm13d = direction `cross` vec13
+        -- The determinate. A volume. AKA, how close to parallel are the triangle, and the plane with normal norm13d.
+        determinate = vec12 `dot` norm13d
+        -- An edge vector in our 'safe' direction, from v1.
+        vec1s = source - v1
+        -- The distance along vec12 to find the intersection of a ray from source in direction direction, and the plane our triangle is on.
+        u = (vec1s `dot` norm13d) / determinate
+        -- The normal of the plane fromed by vec12 and vec1s.
+        norm121s = vec1s `cross` vec12
+        -- The distance along vec13 to find the intersertion of a ray from source in direction direction, and the plane our triangle is on.
+        v = (direction `dot` norm121s) / determinate
+
+-- With inspiration from: https://github.com/RenderKit/embree/blob/master/tutorials/common/math/closest_point.h
+-- FIXME: Sensitive to really skinny triangles; detect these, and select a different 'v1'?.
+distancePointToTriangle :: V3 ℝ -> (V3 ℝ,V3 ℝ,V3 ℝ) -> ℝ
+distancePointToTriangle point (virtex1, virtex2, virtex3) = distance point closestPointToTriangleCenteredSorted
+  where
+    -- FIXME: Here are two precision error mitigation wrappers. Find out if they're needed.
+    -- Reorder triangles such that we use the corner across from the longest side to address the space.
+    closestPointToTriangleCenteredSorted :: V3 ℝ
+    closestPointToTriangleCenteredSorted = closestPointToTriangleCentered adjustedTriangle
+      where
+        adjustedTriangle
+          | abLength >= bcLength && abLength >= caLength = (virtex3, virtex1, virtex2)
+          | abLength >= caLength                         = (virtex1, virtex2, virtex3)
+          | otherwise                                    = (virtex2, virtex3, virtex1)
+          where
+            -- Really, using length-squared. don't have to abs it, don't have to sqrt it.
+            abLength = (virtex2-virtex1) `dot` (virtex2-virtex1)
+            bcLength = (virtex3-virtex2) `dot` (virtex3-virtex2)
+            caLength = (virtex1-virtex3) `dot` (virtex1-virtex3)
+    -- Force closestPointToTriangle to work at the origin
+    closestPointToTriangleCentered :: (V3 ℝ,V3 ℝ,V3 ℝ) -> V3 ℝ
+    closestPointToTriangleCentered (vir1, vir2, vir3) = originDistance + closestPointToTriangle adjustedTriangle adjustedPoint
+      where
+        originDistance = 1/3 *^ (vir1 + vir2 + vir3)
+        adjustedTriangle = (vir1 - originDistance, vir2 - originDistance, vir3 - originDistance)
+        adjustedPoint = point - originDistance
+    closestPointToTriangle :: (V3 ℝ,V3 ℝ,V3 ℝ) -> V3 ℝ -> V3 ℝ
+    closestPointToTriangle (v1, v2, v3) p
+      -- Closest to the virtices
+      | d1 <= 0 && d2 <=  0 = v1
+      | d3 <= 0 && d4 <= d3 = v2
+      | d5 <= 0 && d6 <= d5 = v3
+      -- Closest to the edges
+      | va <= 0 && d1 >=  0 && d3 <= 0 = v1 + (d1 / (d1 - d3)) *^ vec12
+      | vb <= 0 && d2 >=  0 && d6 <= 0 = v1 + (d2 / (d2 - d6)) *^ vec13
+      | vc <= 0 && dx >=  0 && dy >= 0 = v2 + (dx / (dx + dy)) *^ vec23
+      -- On the triangle's surface
+      | otherwise = v1 + (vb * denom *^ vec12) + (vc * denom *^ vec13)
+      where
+        -- The distance along edge12 and edge23, for segment V1 -> P when translated onto the triangle's plane.
+        -- (P when translaned? Read: a line is drawn down to the plane the triangle is on, from p, to a point that is at a right angle with said line.)
+        d1 = vec12 `dot` vec1p
+        d2 = vec13 `dot` vec1p
+        -- Our edge vectors. We have picked v1 to address the space by, for convenience.
+        vec12 = v2 - v1
+        vec13 = v3 - v1
+        -- A segment between our point, and chosen virtex.
+        vec1p =  p - v1
+        -- Distance along edge12 and edge23, for segment V2 -> P when translated onto the triangle's plane.
+        d3 = vec12 `dot` vec2p
+        d4 = vec13 `dot` vec2p
+        -- A segment between our point, and the second virtex.
+        vec2p =  p - v2
+        -- Distance along edge12 and edge23, for segment V3 -> P when translated onto the triangle's plane.
+        d5 = vec12 `dot` vec3p
+        d6 = vec13 `dot` vec3p
+        -- A segment between our point, and the third virtex.
+        vec3p =  p - v3
+        -- An edge vector, along edge23.
+        vec23 = v3 - v2
+        -- The fractional denenominators.
+        va = d1 * d4 - d3 * d2
+        vb = d2 * d5 - d6 * d1
+        vc = d3 * d6 - d5 * d4
+        -- Two convienience values, to make the spacing on the formulas above work.
+        dx = d4 - d3
+        dy = d5 - d6
+        -- The denominator.
+        denom = 1 / va + vb + vc
diff --git a/Graphics/Implicit/Primitives.hs b/Graphics/Implicit/Primitives.hs
index 9d73f161..9ddb6447 100644
--- a/Graphics/Implicit/Primitives.hs
+++ b/Graphics/Implicit/Primitives.hs
@@ -36,6 +36,7 @@ module Graphics.Implicit.Primitives (
                                      ellipsoid,
                                      square, rect,
                                      polygon,
+                                     polyhedron,
                                      rotateExtrude,
                                      rotate3,
                                      rotateQ,
@@ -58,7 +59,7 @@ module Graphics.Implicit.Primitives (
 import Prelude(Applicative, Eq, Foldable, Num, abs, (<), otherwise, Num, (-), (*), (/), (.), negate, Bool(True, False), Maybe(Just, Nothing), Either, fmap, ($), (<=), (&&), max, Ord)
 
 import Graphics.Implicit.Canon (canonicalize2, canonicalize3)
-import Graphics.Implicit.Definitions (ObjectContext, ℝ, ℝ2, ℝ3, Box2,
+import Graphics.Implicit.Definitions (ObjectContext, ℝ, ℝ2, ℝ3, ℕ, Box2,
                                       SharedObj(Empty,
                                                 Full,
                                                 Translate,
@@ -87,6 +88,7 @@ import Graphics.Implicit.Definitions (ObjectContext, ℝ, ℝ2, ℝ3, Box2,
                                                    Cube,
                                                    Sphere,
                                                    Cylinder,
+                                                   Polyhedron,
                                                    Torus,
                                                    BoxFrame,
                                                    Rotate3,
@@ -131,9 +133,16 @@ cube
 cube False size = Cube size
 cube True  size = translate (fmap (negate . (/ 2)) size) $ Cube size
 
+-- | A polyhedron
+polyhedron
+    :: [ℝ3]           -- ^ Points
+    -> [(ℕ,ℕ,ℕ)]      -- ^ triangles, resolved through indexing Points
+    -> SymbolicObj3   -- ^ Resulting polyhedron
+polyhedron points tris = Polyhedron points tris
+
 -- | A conical frustum --- ie. a cylinder with different radii at either end.
-cylinder2 ::
-    ℝ                   -- ^ Radius of the cylinder
+cylinder2
+    :: ℝ                -- ^ Radius of the cylinder
     -> ℝ                -- ^ Second radius of the cylinder
     -> ℝ                -- ^ Height of the cylinder
     -> SymbolicObj3     -- ^ Resulting cylinder
diff --git a/Makefile b/Makefile
index 3e1b8864..e10ac901 100644
--- a/Makefile
+++ b/Makefile
@@ -152,7 +152,7 @@ dist: $(TARGETS)
 
 # Generate examples.
 examples: $(EXTOPENSCADBIN)
-	cd Examples && for each in `find ./ -name '*scad' -type f | sort`; do { echo $$each ; ../$(EXTOPENSCADBIN) $(SCADOPTS) $$each $(RTSOPTS); } done
+	cd Examples && for each in `find ./ -name '*scad' -type f | sort`; do { echo ../$(EXTOPENSCADBIN) $(SCADOPTS) $$each $(RTSOPTS); ../$(EXTOPENSCADBIN) $(SCADOPTS) $$each $(RTSOPTS); } done
 	# NOTE: on debian, if this fails to find the linear package, run: 'apt install libghc-linear-dev libghc-show-combinators-dev libghc-blaze-svg-dev libghc-data-default-dev libghc-juicypixels-dev'
 	cd Examples && for each in `find ./ -name '*.hs' -type f | sort`; do { filename=$(basename "$$each"); filename="$${filename%.*}"; cd ..; $(GHC) $(EXAMPLEOPTS) Examples/$$filename.hs -o Examples/$$filename; cd Examples; echo $$filename; $$filename +RTS -t ; } done
 
diff --git a/tests/Graphics/Implicit/Test/Instances.hs b/tests/Graphics/Implicit/Test/Instances.hs
index 87ece632..d55628e5 100644
--- a/tests/Graphics/Implicit/Test/Instances.hs
+++ b/tests/Graphics/Implicit/Test/Instances.hs
@@ -20,7 +20,7 @@
 
 module Graphics.Implicit.Test.Instances (Observe, observe, (=~=), arbitraryNonZeroV) where
 
-import Prelude (Applicative, (.), not, abs, fmap, Bool(False, True), Bounded, Double, Integer, fromIntegral, (*), (/), (^), round, Enum, Show(show), unlines, Ord, compare, Eq, (==), pure, RealFloat(isNaN), Int, Double, ($), (<), div, (<*>), (<$>), (+), (<>), (<=))
+import Prelude (Applicative, (.), not, return, abs, fmap, Bool(False, True), Bounded, Double, Integer, fromIntegral, (*), (/), (^), round, Enum, Show(show), unlines, Ord, compare, Eq, (==), pure, RealFloat(isNaN), Int, Double, ($), (<), div, (<*>), (<$>), (+), (<>), (<=))
 #if MIN_VERSION_base(4,17,0)
 import Prelude (type(~))
 #endif
@@ -50,6 +50,7 @@ import Graphics.Implicit.Definitions
       ℝ,
       ℝ2,
       ℝ3,
+      ℕ,
       SharedObj(Outset, Translate, Scale, UnionR, IntersectR,
                 DifferenceR, Shell, WithRounding) )
 
@@ -215,6 +216,12 @@ instance Arbitrary (Quaternion ℝ) where
       then discard
       else pure $ axisAngle v q
 
+instance Arbitrary ℕ where
+--  shrink = genericShrink
+  arbitrary = do
+    n <- getPositive <$> arbitrary
+    return n
+
 ------------------------------------------------------------------------------
 -- Minimum of quickspec(s) Observe class and instances required for implicit testsuite
 -- BSD3 Copyright: 2009-2019 Nick Smallbone

From 5dbedd30fcfd35593fbc766563643cd44c54da19 Mon Sep 17 00:00:00 2001
From: Julia Longtin <julia.longtin@gmail.com>
Date: Mon, 23 Feb 2026 16:51:11 +0000
Subject: [PATCH 02/28] first bug fixes.

---
 Graphics/Implicit/ExtOpenScad/Primitives.hs  | 2 +-
 Graphics/Implicit/ObjectUtil/GetImplicit3.hs | 7 ++++---
 2 files changed, 5 insertions(+), 4 deletions(-)

diff --git a/Graphics/Implicit/ExtOpenScad/Primitives.hs b/Graphics/Implicit/ExtOpenScad/Primitives.hs
index f551a84f..8334957e 100644
--- a/Graphics/Implicit/ExtOpenScad/Primitives.hs
+++ b/Graphics/Implicit/ExtOpenScad/Primitives.hs
@@ -284,7 +284,7 @@ polyhedron = moduleWithoutSuite "polyhedron" $ \_ _ -> do
         -- FIXME: use warnC here, instead of error.
         trianglesFromFace :: [ℕ] -> [(ℕ,ℕ,ℕ)]
         trianglesFromFace [p1,p2,p3] = [(p1,p2,p3)]
-        trianglesFromFace (p1:p2:p3:xs) = ((p1,p2,p3):trianglesFromFace (p2:p3:xs))
+        trianglesFromFace (p1:p2:p3:xs) = ((p1,p2,p3):trianglesFromFace (p1:p3:xs))
         trianglesFromFace ys = error $ "too few points:" <> (show ys) <> "\n"
 
 cone :: (Symbol, SourcePosition -> [OVal] -> ArgParser (StateC [OVal]))
diff --git a/Graphics/Implicit/ObjectUtil/GetImplicit3.hs b/Graphics/Implicit/ObjectUtil/GetImplicit3.hs
index a69aac42..855bfc73 100644
--- a/Graphics/Implicit/ObjectUtil/GetImplicit3.hs
+++ b/Graphics/Implicit/ObjectUtil/GetImplicit3.hs
@@ -78,8 +78,8 @@ getImplicit3 _ (Polyhedron [] _) = error "Asked to find distance to an empty pol
 getImplicit3 _ (Polyhedron _ []) = error "Asked to find distance to an empty polygon. No tris."
 getImplicit3 _ (Polyhedron points tris) = \(point) ->
   if insideMesh triangles point
-  then          minimum $ distancePointToTriangle point <$> triangles
-  else negate $ minimum $ distancePointToTriangle point <$> triangles
+  then negate $ minimum $ distancePointToTriangle point <$> triangles
+  else          minimum $ distancePointToTriangle point <$> triangles
   where
     -- decompose our tris into triangles.
     triangles = findTriangle points <$> tris
@@ -241,6 +241,7 @@ insideMesh triangles point = odd $ length $ foundIntersection point unsafeSafeRa
       -- the intersection is not on the triangle (distance along vec12 out of range)
       | u <= 0 || u > 1      = False
       | v <= 0 || v > 1      = False
+      | u + v > 1            = False
       | otherwise            = True
       where
         -- Our edge vectors. We have picked v1 to address the space by, for convenience.
@@ -327,4 +328,4 @@ distancePointToTriangle point (virtex1, virtex2, virtex3) = distance point close
         dx = d4 - d3
         dy = d5 - d6
         -- The denominator.
-        denom = 1 / va + vb + vc
+        denom = 1 / (va + vb + vc)

From 6c5134c64714bbea610d0c5fffccaa4a4680c462 Mon Sep 17 00:00:00 2001
From: Julia Longtin <julia.longtin@gmail.com>
Date: Mon, 23 Feb 2026 17:14:57 +0000
Subject: [PATCH 03/28] avoid dividing 1 by a small number, divide v[a-c]
 instead.

---
 Graphics/Implicit/ObjectUtil/GetImplicit3.hs | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/Graphics/Implicit/ObjectUtil/GetImplicit3.hs b/Graphics/Implicit/ObjectUtil/GetImplicit3.hs
index 855bfc73..ee46909b 100644
--- a/Graphics/Implicit/ObjectUtil/GetImplicit3.hs
+++ b/Graphics/Implicit/ObjectUtil/GetImplicit3.hs
@@ -297,7 +297,7 @@ distancePointToTriangle point (virtex1, virtex2, virtex3) = distance point close
       | vb <= 0 && d2 >=  0 && d6 <= 0 = v1 + (d2 / (d2 - d6)) *^ vec13
       | vc <= 0 && dx >=  0 && dy >= 0 = v2 + (dx / (dx + dy)) *^ vec23
       -- On the triangle's surface
-      | otherwise = v1 + (vb * denom *^ vec12) + (vc * denom *^ vec13)
+      | otherwise = (va / denom *^ v1) + (vb / denom *^ v2) + (vc / denom *^ v3)
       where
         -- The distance along edge12 and edge23, for segment V1 -> P when translated onto the triangle's plane.
         -- (P when translaned? Read: a line is drawn down to the plane the triangle is on, from p, to a point that is at a right angle with said line.)
@@ -328,4 +328,4 @@ distancePointToTriangle point (virtex1, virtex2, virtex3) = distance point close
         dx = d4 - d3
         dy = d5 - d6
         -- The denominator.
-        denom = 1 / (va + vb + vc)
+        denom = va + vb + vc

From 47621500e468747203286228dd950cd43f8d5465 Mon Sep 17 00:00:00 2001
From: Julia Longtin <julia.longtin@gmail.com>
Date: Mon, 23 Feb 2026 21:25:43 +0000
Subject: [PATCH 04/28] cover more edge cases.

---
 Graphics/Implicit/ObjectUtil/GetImplicit3.hs | 41 +++++++++++++-------
 1 file changed, 27 insertions(+), 14 deletions(-)

diff --git a/Graphics/Implicit/ObjectUtil/GetImplicit3.hs b/Graphics/Implicit/ObjectUtil/GetImplicit3.hs
index ee46909b..6d0ca000 100644
--- a/Graphics/Implicit/ObjectUtil/GetImplicit3.hs
+++ b/Graphics/Implicit/ObjectUtil/GetImplicit3.hs
@@ -12,7 +12,9 @@ import Control.Lens ((^.))
 
 import qualified Data.Either as Either (either)
 
-import Data.List(genericIndex, length)
+import Data.List(filter, genericIndex, length)
+
+import Data.Ord (clamp)
 
 import Linear (V2(V2), V3(V3), _xy, _z, distance, dot)
 import qualified Linear (conjugate, inv44, normalizePoint, normalize, point, rotate, Metric)
@@ -229,7 +231,7 @@ getImplicit3 ctx (Shared3 obj) = getImplicitShared ctx obj
 
 -- | Check to see if a point is inside or outside of a mesh.
 insideMesh :: [(V3 ℝ, V3 ℝ, V3 ℝ)] -> V3 ℝ -> Bool
-insideMesh triangles point = odd $ length $ foundIntersection point unsafeSafeRay <$> triangles
+insideMesh triangles point = odd . length . filter (==True) $ foundIntersection point unsafeSafeRay <$> triangles
   where
     -- our assumed safe ray. it's not really safe, but.. gotta start somewhere.
     unsafeSafeRay :: V3 ℝ
@@ -237,13 +239,16 @@ insideMesh triangles point = odd $ length $ foundIntersection point unsafeSafeRa
     foundIntersection :: V3 ℝ -> V3 ℝ -> (V3 ℝ, V3 ℝ, V3 ℝ) -> Bool
     foundIntersection source direction (v1, v2, v3)
       -- shouldn't happen, triangle is on the same plane as unsafeSafeRay
-      | abs determinate <= 0 = False
+      | abs determinate <= eps = False
       -- the intersection is not on the triangle (distance along vec12 out of range)
-      | u <= 0 || u > 1      = False
-      | v <= 0 || v > 1      = False
-      | u + v > 1            = False
-      | otherwise            = True
+      | u <= eps || u > 1+eps  = False
+      | v <= eps || v > 1+eps  = False
+      | u + v > 1+eps          = False
+      | otherwise              = t > 0
       where
+        -- fudge factor.
+        eps :: ℝ
+        eps = 1e-12
         -- Our edge vectors. We have picked v1 to address the space by, for convenience.
         vec12 = v2 - v1
         vec13 = v3 - v1
@@ -259,6 +264,8 @@ insideMesh triangles point = odd $ length $ foundIntersection point unsafeSafeRa
         norm121s = vec1s `cross` vec12
         -- The distance along vec13 to find the intersertion of a ray from source in direction direction, and the plane our triangle is on.
         v = (direction `dot` norm121s) / determinate
+        -- The distance along the ray to find the intersertion of a ray from source in direction direction, and the plane our triangle is on.
+        t = (vec13 `dot` norm121s) / determinate
 
 -- With inspiration from: https://github.com/RenderKit/embree/blob/master/tutorials/common/math/closest_point.h
 -- FIXME: Sensitive to really skinny triangles; detect these, and select a different 'v1'?.
@@ -289,16 +296,19 @@ distancePointToTriangle point (virtex1, virtex2, virtex3) = distance point close
     closestPointToTriangle :: (V3 ℝ,V3 ℝ,V3 ℝ) -> V3 ℝ -> V3 ℝ
     closestPointToTriangle (v1, v2, v3) p
       -- Closest to the virtices
-      | d1 <= 0 && d2 <=  0 = v1
-      | d3 <= 0 && d4 <= d3 = v2
-      | d5 <= 0 && d6 <= d5 = v3
+      | d1 <=  eps && d2 <= eps = v1
+      | d3 >= -eps && d4 <=  d3 = v2
+      | d5 >= -eps && d6 <=  d5 = v3
       -- Closest to the edges
-      | va <= 0 && d1 >=  0 && d3 <= 0 = v1 + (d1 / (d1 - d3)) *^ vec12
-      | vb <= 0 && d2 >=  0 && d6 <= 0 = v1 + (d2 / (d2 - d6)) *^ vec13
-      | vc <= 0 && dx >=  0 && dy >= 0 = v2 + (dx / (dx + dy)) *^ vec23
+      | vc <= eps && d1 >= -eps && d3 <=  eps = v1 + (d1 / (d1 - d3)) *^ vec12
+      | vb <= eps && d2 >= -eps && d6 <=  eps = v1 + (d2 / (d2 - d6)) *^ vec13
+      | va <= eps && dx >= -eps && dy >= -eps = v2 + (dx / (dx + dy)) *^ vec23
       -- On the triangle's surface
-      | otherwise = (va / denom *^ v1) + (vb / denom *^ v2) + (vc / denom *^ v3)
+      | otherwise = v1 + v *^ vec12 + w *^ vec13
       where
+        -- fudge factor.
+        eps :: ℝ
+        eps = 1e-12
         -- The distance along edge12 and edge23, for segment V1 -> P when translated onto the triangle's plane.
         -- (P when translaned? Read: a line is drawn down to the plane the triangle is on, from p, to a point that is at a right angle with said line.)
         d1 = vec12 `dot` vec1p
@@ -329,3 +339,6 @@ distancePointToTriangle point (virtex1, virtex2, virtex3) = distance point close
         dy = d5 - d6
         -- The denominator.
         denom = va + vb + vc
+        -- barycentric results, where we actually intersect.
+        v = clamp (0,1) $ vb / denom
+        w = clamp (0,1) $ vc / denom

From ff39421121eb3b829e66bc8423de53b38a12c931 Mon Sep 17 00:00:00 2001
From: Julia Longtin <julia.longtin@gmail.com>
Date: Mon, 23 Feb 2026 23:56:47 +0000
Subject: [PATCH 05/28] add a note for next time I am debugging.

---
 implicit.cabal | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/implicit.cabal b/implicit.cabal
index 49f403fd..9671dbd6 100644
--- a/implicit.cabal
+++ b/implicit.cabal
@@ -53,6 +53,8 @@ Common binstuff
                 -threaded
                 -rtsopts "-with-rtsopts -N -qg -t"
                 -optc-O3
+                -- Call the heap checker, even when memory is not being allocated. needed for control-c to work right when debugging.
+                -- -fno-omit-yields
                 -- see GHC manual 8.2.1 section 6.5.1.
                 -feager-blackholing
                 -Wall
@@ -73,6 +75,8 @@ Common libstuff
     Default-extensions: NoImplicitPrelude
     Ghc-options:
                 -optc-O3
+                -- Call the heap checker, even when memory is not being allocated. needed for control-c to work right when debugging.
+                -- -fno-omit-yields
                 -- see GHC manual 8.2.1 section 6.5.1.
                 -feager-blackholing
                 -Wall

From 904ea9cc9f122892d250e0f1896d1e0b7b8741dd Mon Sep 17 00:00:00 2001
From: Julia Longtin <julia.longtin@gmail.com>
Date: Mon, 23 Feb 2026 23:57:10 +0000
Subject: [PATCH 06/28] rotated math is less spikey?

---
 Graphics/Implicit/ObjectUtil/GetImplicit3.hs | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/Graphics/Implicit/ObjectUtil/GetImplicit3.hs b/Graphics/Implicit/ObjectUtil/GetImplicit3.hs
index 6d0ca000..40027d53 100644
--- a/Graphics/Implicit/ObjectUtil/GetImplicit3.hs
+++ b/Graphics/Implicit/ObjectUtil/GetImplicit3.hs
@@ -273,14 +273,14 @@ distancePointToTriangle :: V3 ℝ -> (V3 ℝ,V3 ℝ,V3 ℝ) -> ℝ
 distancePointToTriangle point (virtex1, virtex2, virtex3) = distance point closestPointToTriangleCenteredSorted
   where
     -- FIXME: Here are two precision error mitigation wrappers. Find out if they're needed.
-    -- Reorder triangles such that we use the corner across from the longest side to address the space.
+    -- Reorder triangles such that we use one of the corners of the longest side to address the space in barycentric coordinates. Yes, this is wrong, but it's much more numerically stable?
     closestPointToTriangleCenteredSorted :: V3 ℝ
     closestPointToTriangleCenteredSorted = closestPointToTriangleCentered adjustedTriangle
       where
         adjustedTriangle
-          | abLength >= bcLength && abLength >= caLength = (virtex3, virtex1, virtex2)
-          | abLength >= caLength                         = (virtex1, virtex2, virtex3)
-          | otherwise                                    = (virtex2, virtex3, virtex1)
+          | abLength >= bcLength && abLength >= caLength = (virtex2, virtex3, virtex1)
+          | abLength >= caLength                         = (virtex3, virtex1, virtex2)
+          | otherwise                                    = (virtex1, virtex2, virtex3)
           where
             -- Really, using length-squared. don't have to abs it, don't have to sqrt it.
             abLength = (virtex2-virtex1) `dot` (virtex2-virtex1)

From 148cbb703dfdcd451c0a84fe5ea6877f66cdebce Mon Sep 17 00:00:00 2001
From: Julia Longtin <julia.longtin@gmail.com>
Date: Tue, 24 Feb 2026 18:55:32 +0000
Subject: [PATCH 07/28] adjust EPS values and handling, remove clamping.

---
 Graphics/Implicit/ObjectUtil/GetImplicit3.hs | 49 ++++++++++----------
 1 file changed, 24 insertions(+), 25 deletions(-)

diff --git a/Graphics/Implicit/ObjectUtil/GetImplicit3.hs b/Graphics/Implicit/ObjectUtil/GetImplicit3.hs
index 40027d53..a9b3c1b5 100644
--- a/Graphics/Implicit/ObjectUtil/GetImplicit3.hs
+++ b/Graphics/Implicit/ObjectUtil/GetImplicit3.hs
@@ -14,8 +14,6 @@ import qualified Data.Either as Either (either)
 
 import Data.List(filter, genericIndex, length)
 
-import Data.Ord (clamp)
-
 import Linear (V2(V2), V3(V3), _xy, _z, distance, dot)
 import qualified Linear (conjugate, inv44, normalizePoint, normalize, point, rotate, Metric)
 
@@ -230,6 +228,7 @@ getImplicit3 ctx (RotateExtrude totalRotation translate rotate symbObj) =
 getImplicit3 ctx (Shared3 obj) = getImplicitShared ctx obj
 
 -- | Check to see if a point is inside or outside of a mesh.
+-- FIXME: Replace this with a winding number based approach. See: https://github.com/marmakoide/inside-3d-mesh/blob/master/README.md
 insideMesh :: [(V3 ℝ, V3 ℝ, V3 ℝ)] -> V3 ℝ -> Bool
 insideMesh triangles point = odd . length . filter (==True) $ foundIntersection point unsafeSafeRay <$> triangles
   where
@@ -239,16 +238,17 @@ insideMesh triangles point = odd . length . filter (==True) $ foundIntersection
     foundIntersection :: V3 ℝ -> V3 ℝ -> (V3 ℝ, V3 ℝ, V3 ℝ) -> Bool
     foundIntersection source direction (v1, v2, v3)
       -- shouldn't happen, triangle is on the same plane as unsafeSafeRay
-      | abs determinate <= eps = False
+      | abs determinate < eps = False
       -- the intersection is not on the triangle (distance along vec12 out of range)
-      | u <= eps || u > 1+eps  = False
-      | v <= eps || v > 1+eps  = False
-      | u + v > 1+eps          = False
-      | otherwise              = t > 0
+      | u < eps || u > 1-eps  = False
+      | v < eps || v > 1-eps  = False
+      | u + v > 1-eps-eps     = False
+      | otherwise             = t > 0
       where
         -- fudge factor.
+        -- adjusted from: 6,7,9,12,15,17,18==(43,171), 16 ==(38,171)
         eps :: ℝ
-        eps = 1e-12
+        eps = 1e-16
         -- Our edge vectors. We have picked v1 to address the space by, for convenience.
         vec12 = v2 - v1
         vec13 = v3 - v1
@@ -262,25 +262,23 @@ insideMesh triangles point = odd . length . filter (==True) $ foundIntersection
         u = (vec1s `dot` norm13d) / determinate
         -- The normal of the plane fromed by vec12 and vec1s.
         norm121s = vec1s `cross` vec12
-        -- The distance along vec13 to find the intersertion of a ray from source in direction direction, and the plane our triangle is on.
+        -- The distance along vec13 to find the intersection of a ray from source in direction direction, and the plane our triangle is on.
         v = (direction `dot` norm121s) / determinate
-        -- The distance along the ray to find the intersertion of a ray from source in direction direction, and the plane our triangle is on.
+        -- The distance along the ray to find the intersection of a ray from source in direction direction, and the plane our triangle is on.
         t = (vec13 `dot` norm121s) / determinate
 
 -- With inspiration from: https://github.com/RenderKit/embree/blob/master/tutorials/common/math/closest_point.h
--- FIXME: Sensitive to really skinny triangles; detect these, and select a different 'v1'?.
 distancePointToTriangle :: V3 ℝ -> (V3 ℝ,V3 ℝ,V3 ℝ) -> ℝ
 distancePointToTriangle point (virtex1, virtex2, virtex3) = distance point closestPointToTriangleCenteredSorted
   where
-    -- FIXME: Here are two precision error mitigation wrappers. Find out if they're needed.
-    -- Reorder triangles such that we use one of the corners of the longest side to address the space in barycentric coordinates. Yes, this is wrong, but it's much more numerically stable?
+    -- Reorder triangles such that we use the corner away from the longest side to address the space in barycentric coordinates.
     closestPointToTriangleCenteredSorted :: V3 ℝ
     closestPointToTriangleCenteredSorted = closestPointToTriangleCentered adjustedTriangle
       where
         adjustedTriangle
-          | abLength >= bcLength && abLength >= caLength = (virtex2, virtex3, virtex1)
-          | abLength >= caLength                         = (virtex3, virtex1, virtex2)
-          | otherwise                                    = (virtex1, virtex2, virtex3)
+          | abLength >= bcLength && abLength >= caLength = (virtex3, virtex1, virtex2)
+          | abLength >= caLength                         = (virtex1, virtex2, virtex3)
+          | otherwise                                    = (virtex2, virtex3, virtex1)
           where
             -- Really, using length-squared. don't have to abs it, don't have to sqrt it.
             abLength = (virtex2-virtex1) `dot` (virtex2-virtex1)
@@ -297,18 +295,19 @@ distancePointToTriangle point (virtex1, virtex2, virtex3) = distance point close
     closestPointToTriangle (v1, v2, v3) p
       -- Closest to the virtices
       | d1 <=  eps && d2 <= eps = v1
-      | d3 >= -eps && d4 <=  d3 = v2
-      | d5 >= -eps && d6 <=  d5 = v3
+      | d3  > -eps && d4 <=  d3 = v2
+      | d6  > -eps && d5 <=  d6 = v3
       -- Closest to the edges
-      | vc <= eps && d1 >= -eps && d3 <=  eps = v1 + (d1 / (d1 - d3)) *^ vec12
-      | vb <= eps && d2 >= -eps && d6 <=  eps = v1 + (d2 / (d2 - d6)) *^ vec13
-      | va <= eps && dx >= -eps && dy >= -eps = v2 + (dx / (dx + dy)) *^ vec23
+      | vc <= eps && d1 > -eps && d3 <=  eps = v1 + (d1 / (d1 - d3)) *^ vec12
+      | vb <= eps && d2 > -eps && d6 <=  eps = v1 + (d2 / (d2 - d6)) *^ vec13
+      | va <= eps && dx > -eps && dy  > -eps = v2 + (dx / (dx + dy)) *^ vec23
       -- On the triangle's surface
       | otherwise = v1 + v *^ vec12 + w *^ vec13
       where
-        -- fudge factor.
+        -- fudge factor. chosen by tuning with our unit pyramid.
+        -- 3,11 = 52 , 12 == 50, 13,14,15,18 == 45
         eps :: ℝ
-        eps = 1e-12
+        eps = 1e-13
         -- The distance along edge12 and edge23, for segment V1 -> P when translated onto the triangle's plane.
         -- (P when translaned? Read: a line is drawn down to the plane the triangle is on, from p, to a point that is at a right angle with said line.)
         d1 = vec12 `dot` vec1p
@@ -340,5 +339,5 @@ distancePointToTriangle point (virtex1, virtex2, virtex3) = distance point close
         -- The denominator.
         denom = va + vb + vc
         -- barycentric results, where we actually intersect.
-        v = clamp (0,1) $ vb / denom
-        w = clamp (0,1) $ vc / denom
+        v = vb / denom
+        w = vc / denom

From 40be7ced4aadbda7c6c431f2d328b3f5b38a6693 Mon Sep 17 00:00:00 2001
From: Julia Longtin <julia.longtin@gmail.com>
Date: Wed, 25 Feb 2026 17:04:18 +0000
Subject: [PATCH 08/28] add example for Polyhedron.

---
 Examples/example26-Polyhedron.scad | 24 ++++++++++++++++++++++++
 1 file changed, 24 insertions(+)
 create mode 100644 Examples/example26-Polyhedron.scad

diff --git a/Examples/example26-Polyhedron.scad b/Examples/example26-Polyhedron.scad
new file mode 100644
index 00000000..2471f125
--- /dev/null
+++ b/Examples/example26-Polyhedron.scad
@@ -0,0 +1,24 @@
+module pyramid(base, height) {
+  half = base / 2;
+
+  polyhedron(
+    points = [
+      [-half, -half, 0],
+      [ half, -half, 0],
+      [ half,  half, 0],
+      [-half,  half, 0],
+      [0, 0, height]
+    ],
+    faces = [
+      [0, 1, 2, 3],
+      [0, 1, 4],
+      [1, 2, 4],
+      [2, 3, 4],
+      [3, 0, 4]
+    ]
+  );
+}
+
+pyramid_base = 12;
+pyramid_height = 14;
+pyramid(pyramid_base, pyramid_height);

From 656eec02972f39df0a166f6a99209b327552aa12 Mon Sep 17 00:00:00 2001
From: Julia Longtin <julia.longtin@gmail.com>
Date: Sat, 28 Feb 2026 02:19:40 +0000
Subject: [PATCH 09/28] remove warning, and unnecessary import.

---
 Graphics/Implicit/ExtOpenScad/Primitives.hs | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Graphics/Implicit/ExtOpenScad/Primitives.hs b/Graphics/Implicit/ExtOpenScad/Primitives.hs
index 8334957e..d429a6c6 100644
--- a/Graphics/Implicit/ExtOpenScad/Primitives.hs
+++ b/Graphics/Implicit/ExtOpenScad/Primitives.hs
@@ -27,7 +27,7 @@ import qualified Graphics.Implicit.ExtOpenScad.Util.ArgParser as GIEUA (argument
 
 import Graphics.Implicit.ExtOpenScad.Util.OVal (OTypeMirror, caseOType, divideObjs, (<||>))
 
-import Graphics.Implicit.ExtOpenScad.Util.StateC (errorC, warnC)
+import Graphics.Implicit.ExtOpenScad.Util.StateC (errorC)
 
 -- Note the use of a qualified import, so we don't have the functions in this file conflict with what we're importing.
 import qualified Graphics.Implicit.Primitives as Prim (withRounding, sphere, rect3, rect, translate, circle, polygon, polyhedron, extrude, cylinder2, union, unionR, intersect, intersectR, difference, differenceR, rotate, slice, transform, rotate3V, rotate3, transform3, scale, extrudeM, rotateExtrude, shell, mirror, pack3, pack2, torus, ellipsoid, cone)

From c7a53dc0edc735eb300ceea32f8e66b7295550b8 Mon Sep 17 00:00:00 2001
From: Julia Longtin <julia.longtin@gmail.com>
Date: Sat, 7 Mar 2026 15:09:41 +0000
Subject: [PATCH 10/28] split ONModule into with and without suite variants,
 simplify Primitives.hs, and use errorC/warnC from Primitives.hs

---
 Graphics/Implicit/ExtOpenScad/Definitions.hs  |  34 ++--
 .../Implicit/ExtOpenScad/Eval/Statement.hs    |  38 ++--
 Graphics/Implicit/ExtOpenScad/Primitives.hs   | 181 ++++++++++--------
 Graphics/Implicit/ExtOpenScad/Util/OVal.hs    |   5 +-
 Graphics/Implicit/ObjectUtil/GetImplicit3.hs  |   2 +-
 programs/docgen.hs                            |   7 +-
 6 files changed, 148 insertions(+), 119 deletions(-)

diff --git a/Graphics/Implicit/ExtOpenScad/Definitions.hs b/Graphics/Implicit/ExtOpenScad/Definitions.hs
index 73e4b20f..352943d4 100644
--- a/Graphics/Implicit/ExtOpenScad/Definitions.hs
+++ b/Graphics/Implicit/ExtOpenScad/Definitions.hs
@@ -20,7 +20,7 @@ module Graphics.Implicit.ExtOpenScad.Definitions (ArgParser(AP, APTest, APBranch
                                                   Expr(LitE, Var, ListE, LamE, (:$)),
                                                   StatementI(StatementI),
                                                   Statement(DoNothing, NewModule, Include, If, ModuleCall, (:=)),
-                                                  OVal(OIO, ONum, OBool, OString, OList, OFunc, OUndefined, OUModule, ONModule, OVargsModule, OError, OObj2, OObj3),
+                                                  OVal(OIO, ONum, OBool, OString, OList, OFunc, OUndefined, OUModule, ONModule, ONModuleWithSuite, OVargsModule, OError, OObj2, OObj3),
                                                   TestInvariant(EulerCharacteristic),
                                                   SourcePosition(SourcePosition),
                                                   StateC,
@@ -37,7 +37,7 @@ module Graphics.Implicit.ExtOpenScad.Definitions (ArgParser(AP, APTest, APBranch
                                                   CanCompState'
                                                   ) where
 
-import Prelude(Eq, Show, Ord, Maybe(Just), Bool(True, False), IO, FilePath, (==), show, ($), (<>), and, zipWith, Int, (<$>))
+import Prelude(Eq, Show, Ord, Maybe, Bool(True, False), IO, FilePath, (==), show, ($), (<>), and, zipWith, Int, (<$>))
 
 -- Resolution of the world, Integer type, and symbolic languages for 2D and 3D objects.
 import Graphics.Implicit.Definitions (ℝ, ℕ, Fastℕ, SymbolicObj2, SymbolicObj3, fromFastℕ)
@@ -119,6 +119,7 @@ data ArgParser a
                  -- ^ For returns: @APTerminator (return value)@
                  | APFail Text
                  -- ^ For failure: @APFail (error message)@
+                 -- NOTE: we don't use APFail to fail parsing in Primitives.hs, we have errorC.
                  | APExample Text (ArgParser a)
                  -- ^ An example, then next
                  | APTest Text [TestInvariant] (ArgParser a)
@@ -197,8 +198,10 @@ data OVal = OUndefined
          | OIO (IO OVal)
          -- Name, arguments, argument parsers.
          | OUModule Symbol (Maybe [(Symbol, Bool)]) (VarLookup -> ArgParser (StateC [OVal]))
-         -- Name, implementation, arguments, whether the module accepts/requires a suite.
-         | ONModule Symbol (SourcePosition -> [OVal] -> ArgParser (StateC [OVal])) [([(Symbol, Bool)],  Maybe Bool)]
+         -- Name, implementation, arguments.
+         | ONModule Symbol (SourcePosition -> ArgParser (StateC [OVal])) [[(Symbol, Bool)]]
+         -- Name, implementation, arguments.
+         | ONModuleWithSuite Symbol (SourcePosition -> [OVal] -> ArgParser (StateC [OVal])) [[(Symbol, Bool)]]
          | OVargsModule Symbol (Symbol -> SourcePosition -> [(Maybe Symbol, OVal)] -> [StatementI] -> ([StatementI] -> StateC ()) -> StateC ())
          | OObj3 SymbolicObj3
          | OObj2 SymbolicObj2
@@ -230,18 +233,23 @@ instance Show OVal where
         showArg (Symbol a, hasDefault) = if hasDefault
                                          then a
                                          else a <> "=..."
-        showInstances :: [([(Symbol, Bool)], Maybe Bool)] -> Text
+        showInstances :: [[(Symbol, Bool)]] -> Text
         showInstances [] = ""
         showInstances [oneInstance] = "module " <> name <> showInstance oneInstance
         showInstances multipleInstances = "Module " <> name <> "[ " <> intercalate ", " (showInstance <$> multipleInstances) <> " ]"
-        showInstance :: ([(Symbol, Bool)], Maybe Bool) -> Text
-        showInstance (arguments, suiteInfo) = " (" <> intercalate ", " (showArg <$> arguments) <> ") {}" <> showSuiteInfo suiteInfo
-        showSuiteInfo :: Maybe Bool -> Text
-        showSuiteInfo suiteInfo = case suiteInfo of
-                          Just requiresSuite -> if requiresSuite
-                                                then " requiring suite {}"
-                                                else " accepting suite {}"
-                          _ -> ""
+        showInstance :: [(Symbol, Bool)] -> Text
+        showInstance arguments = " (" <> intercalate ", " (showArg <$> arguments) <> ") "
+    show (ONModuleWithSuite (Symbol name) _ instances) = unpack $ showInstances instances
+      where
+        showArg (Symbol a, hasDefault) = if hasDefault
+                                         then a
+                                         else a <> "=..."
+        showInstances :: [[(Symbol, Bool)]] -> Text
+        showInstances [] = ""
+        showInstances [oneInstance] = "module " <> name <> showInstance oneInstance
+        showInstances multipleInstances = "Module " <> name <> "[ " <> intercalate ", " (showInstance <$> multipleInstances) <> " ]"
+        showInstance :: [(Symbol, Bool)] -> Text
+        showInstance arguments = " (" <> intercalate ", " (showArg <$> arguments) <> ") {} requiring suite {}"
     show (OVargsModule (Symbol name) _) = "varargs module " <> unpack name
     show (OError msg) = unpack $ "Execution Error:\n" <> msg
     show (OObj2 obj) = "<obj2: " <> show obj <> ">"
diff --git a/Graphics/Implicit/ExtOpenScad/Eval/Statement.hs b/Graphics/Implicit/ExtOpenScad/Eval/Statement.hs
index d2701e09..4955c71d 100644
--- a/Graphics/Implicit/ExtOpenScad/Eval/Statement.hs
+++ b/Graphics/Implicit/ExtOpenScad/Eval/Statement.hs
@@ -13,7 +13,7 @@ import Graphics.Implicit.ExtOpenScad.Definitions (
                                                   Statement(Include, (:=), If, NewModule, ModuleCall, DoNothing),
                                                   Pattern(Name),
                                                   Expr(LitE),
-                                                  OVal(OBool, OUModule, ONModule, OVargsModule),
+                                                  OVal(OBool, OUModule, ONModule, ONModuleWithSuite, OVargsModule),
                                                   VarLookup(VarLookup),
                                                   StatementI(StatementI),
                                                   Symbol(Symbol),
@@ -117,14 +117,23 @@ runStatementI (StatementI sourcePos (ModuleCall (Symbol name) argsExpr suite)) =
             Just (ONModule _ implementation forms) -> do
               possibleInstances <- selectInstances forms
               let
-                suiteInfo = case possibleInstances of
-                  [(_, suiteInfoFound)]   -> suiteInfoFound
-                  []                      -> Nothing
-                  ((_, suiteInfoFound):_) -> suiteInfoFound
               when (null possibleInstances) (do
                                                 errorC sourcePos $ "no instance of " <> name <> " found to match given parameters.\nInstances available:\n" <> pack (show (ONModule (Symbol name) implementation forms))
-                                                traverse_ ((`checkOptions` True) . Just . fst) forms
+                                                traverse_ ((`checkOptions` True) . Just) forms
                                             )
+              -- Evaluate all of the arguments.
+              evaluatedArgs <- evalArgs argsExpr
+              -- Run the module.
+              let
+                argsMapped = argMap evaluatedArgs $ implementation sourcePos
+              for_ (pack <$> snd argsMapped) $ errorC sourcePos
+              fromMaybe (pure []) (fst argsMapped)
+            Just (ONModuleWithSuite _ implementation forms) -> do
+              possibleInstances <- selectInstances forms
+              let
+              when (null possibleInstances) $ do
+                                              errorC sourcePos $ "no instance of " <> name <> " found to match given parameters.\nInstances available:\n" <> pack (show (ONModuleWithSuite (Symbol name) implementation forms))
+                                              traverse_ ((`checkOptions` True) . Just) forms
               -- Ignore this for now, because all instances we define have the same suite requirements.
               {-
               when (length possibleInstances > 1) (do
@@ -135,14 +144,9 @@ runStatementI (StatementI sourcePos (ModuleCall (Symbol name) argsExpr suite)) =
               evaluatedArgs <- evalArgs argsExpr
               -- Evaluate the suite.
               vals <- runSuiteCapture varlookup suite
-              suiteResults <- case suiteInfo of
-                              Just True -> do
-                                when (null vals) (errorC sourcePos "Suite required, but none provided.")
-                                pure vals
-                              Just False -> pure vals
-                              _ -> do
-                                when (suite /= []) (errorC sourcePos $ "Suite provided, but module " <> name <> " does not accept one. Perhaps a missing semicolon?")
-                                pure []
+              suiteResults <- do
+                              when (null vals) (errorC sourcePos "Suite required, but none provided.")
+                              pure vals
               -- Run the module.
               let
                 argsMapped = argMap evaluatedArgs $ implementation sourcePos suiteResults
@@ -164,12 +168,12 @@ runStatementI (StatementI sourcePos (ModuleCall (Symbol name) argsExpr suite)) =
                 pure []
         pushVals newVals
           where
-            selectInstances :: [([(Symbol, Bool)], Maybe Bool)] -> StateC [([(Symbol, Bool)], Maybe Bool)]
+            selectInstances :: [[(Symbol, Bool)]] -> StateC [[(Symbol, Bool)]]
             selectInstances instances = do
               validInstances <- for instances
-                    ( \(args, suiteInfo) -> do
+                    ( \args -> do
                         res <- checkOptions (Just args) False
-                        pure $ if res then Just (args, suiteInfo) else Nothing
+                        pure $ if res then Just args else Nothing
                     )
               pure $ catMaybes validInstances
             checkOptions :: Maybe [(Symbol, Bool)] -> Bool -> StateC Bool
diff --git a/Graphics/Implicit/ExtOpenScad/Primitives.hs b/Graphics/Implicit/ExtOpenScad/Primitives.hs
index d429a6c6..4956b46f 100644
--- a/Graphics/Implicit/ExtOpenScad/Primitives.hs
+++ b/Graphics/Implicit/ExtOpenScad/Primitives.hs
@@ -15,11 +15,11 @@
 -- Export one set containing all of the primitive modules.
 module Graphics.Implicit.ExtOpenScad.Primitives (primitiveModules) where
 
-import Prelude((.), Either(Left, Right), Bool(True, False), Maybe(Just, Nothing), ($), pure, show, either, id, (-), (==), (&&), (<), (*), cos, sin, pi, (/), (>), const, uncurry, (/=), (||), not, null, fmap, (<>), otherwise, error, (<*>), (<$>))
+import Prelude(concat, mapM, (.), Either(Left, Right), Bool(True, False), Maybe(Just, Nothing), ($), pure, show, either, id, (-), (==), (&&), (<), (*), cos, sin, pi, (/), (>), const, uncurry, (/=), (||), not, null, fmap, (<>), otherwise, error, (<*>), (<$>))
 
 import Graphics.Implicit.Definitions (ℝ, ℝ2, ℝ3, ℕ, SymbolicObj2, SymbolicObj3, ExtrudeMScale(C1), fromℕtoℝ, isScaleID)
 
-import Graphics.Implicit.ExtOpenScad.Definitions (OVal (OObj2, OObj3, ONModule), ArgParser(APFail), Symbol(Symbol), StateC, SourcePosition)
+import Graphics.Implicit.ExtOpenScad.Definitions (OVal (OObj2, OObj3, ONModule, ONModuleWithSuite), ArgParser, Symbol(Symbol), StateC, SourcePosition)
 
 import Graphics.Implicit.ExtOpenScad.Util.ArgParser (doc, defaultTo, example, test, eulerCharacteristic)
 
@@ -27,16 +27,15 @@ import qualified Graphics.Implicit.ExtOpenScad.Util.ArgParser as GIEUA (argument
 
 import Graphics.Implicit.ExtOpenScad.Util.OVal (OTypeMirror, caseOType, divideObjs, (<||>))
 
-import Graphics.Implicit.ExtOpenScad.Util.StateC (errorC)
+import Graphics.Implicit.ExtOpenScad.Util.StateC (errorC, warnC)
 
 -- Note the use of a qualified import, so we don't have the functions in this file conflict with what we're importing.
 import qualified Graphics.Implicit.Primitives as Prim (withRounding, sphere, rect3, rect, translate, circle, polygon, polyhedron, extrude, cylinder2, union, unionR, intersect, intersectR, difference, differenceR, rotate, slice, transform, rotate3V, rotate3, transform3, scale, extrudeM, rotateExtrude, shell, mirror, pack3, pack2, torus, ellipsoid, cone)
 
-import Data.List (concatMap)
-
-import Control.Monad (when, mplus)
+import Control.Monad (mplus)
 
 import Data.Text.Lazy (Text)
+import qualified Data.Text.Lazy as DTL (pack)
 
 import Control.Lens ((^.))
 import Linear (_m33, M34, M44, V2(V2), V3(V3), V4(V4))
@@ -52,44 +51,52 @@ argument a = GIEUA.argument (Symbol a)
 primitiveModules :: [(Symbol, OVal)]
 primitiveModules =
   [
-    onModIze sphere [([("r", noDefault)], noSuite), ([("d", noDefault)], noSuite)]
-  , onModIze cube [([("x", noDefault), ("y", noDefault), ("z", noDefault), ("center", hasDefault), ("r", hasDefault)], noSuite),([("size", noDefault), ("center", hasDefault), ("r", hasDefault)], noSuite)]
-  , onModIze square [([("x", noDefault), ("y", noDefault), ("center", hasDefault), ("r", hasDefault)], noSuite), ([("size", noDefault), ("center", hasDefault), ("r", hasDefault)], noSuite)]
-  , onModIze cylinder [([("r", hasDefault), ("h", hasDefault), ("r1", hasDefault), ("r2", hasDefault), ("$fn", hasDefault), ("center", hasDefault)], noSuite),
-                       ([("d", hasDefault), ("h", hasDefault), ("d1", hasDefault), ("d2", hasDefault), ("$fn", hasDefault), ("center", hasDefault)], noSuite)]
-  , onModIze circle [([("r", noDefault), ("$fn", hasDefault)], noSuite), ([("d", noDefault), ("$fn", hasDefault)], noSuite)]
-  , onModIze cone [([("r", noDefault), ("h", hasDefault), ("center", hasDefault)], noSuite), ([("d", noDefault), ("h", hasDefault), ("center", hasDefault)], noSuite)]
-  , onModIze torus [([("r1", noDefault), ("r2", hasDefault)], noSuite)]
-  , onModIze ellipsoid [([("a", noDefault), ("b", hasDefault), ("c", hasDefault)], noSuite)]
-  , onModIze polygon [([("points", noDefault)], noSuite)]
-  , onModIze polyhedron [([("points", noDefault), ("faces", noDefault)], noSuite)]
-  , onModIze union [([("r", hasDefault)], requiredSuite)]
-  , onModIze intersect [([("r", hasDefault)], requiredSuite)]
-  , onModIze difference [([("r", hasDefault)], requiredSuite)]
-  , onModIze translate [([("x", noDefault), ("y", noDefault), ("z", noDefault)], requiredSuite), ([("v", noDefault)], requiredSuite)]
-  , onModIze rotate [([("a", noDefault), ("v", hasDefault)], requiredSuite)]
-  , onModIze scale [([("v", noDefault)], requiredSuite)]
-  , onModIze extrude [([("height", hasDefault), ("center", hasDefault), ("twist", hasDefault), ("scale", hasDefault), ("translate", hasDefault), ("r", hasDefault)], requiredSuite)]
-  , onModIze rotateExtrude [([("angle", hasDefault), ("r", hasDefault), ("translate", hasDefault), ("rotate", hasDefault)], requiredSuite)]
-  , onModIze shell [([("w", noDefault)], requiredSuite)]
-  , onModIze projection [([("cut", hasDefault)], requiredSuite)]
-  , onModIze pack [([("size", noDefault), ("sep", noDefault)], requiredSuite)]
-  , onModIze unit [([("unit", noDefault)], requiredSuite)]
-  , onModIze mirror [([("x", noDefault), ("y", noDefault), ("z", noDefault)], requiredSuite), ([("v", noDefault)], requiredSuite)]
-  , onModIze multmatrix [([("m", noDefault)], requiredSuite)]
+    consModule sphere [[("r", noDefault)], [("d", noDefault)]]
+  , consModule cube [[("x", noDefault), ("y", noDefault), ("z", noDefault), ("center", hasDefault), ("r", hasDefault)],[("size", noDefault), ("center", hasDefault), ("r", hasDefault)]]
+  , consModule square [[("x", noDefault), ("y", noDefault), ("center", hasDefault), ("r", hasDefault)], [("size", noDefault), ("center", hasDefault), ("r", hasDefault)]]
+  , consModule cylinder [[("r", hasDefault), ("h", hasDefault), ("r1", hasDefault), ("r2", hasDefault), ("$fn", hasDefault), ("center", hasDefault)],
+                                   [("d", hasDefault), ("h", hasDefault), ("d1", hasDefault), ("d2", hasDefault), ("$fn", hasDefault), ("center", hasDefault)]]
+  , consModule circle [[("r", noDefault), ("$fn", hasDefault)], [("d", noDefault), ("$fn", hasDefault)]]
+  , consModule cone [[("r", noDefault), ("h", hasDefault), ("center", hasDefault)], [("d", noDefault), ("h", hasDefault), ("center", hasDefault)]]
+  , consModule torus [[("r1", noDefault), ("r2", hasDefault)]]
+  , consModule ellipsoid [[("a", noDefault), ("b", hasDefault), ("c", hasDefault)]]
+  , consModule polygon [[("points", noDefault)]]
+  , consModule polyhedron [[("points", noDefault), ("faces", noDefault)]]
+  , consModuleWithSuite union [[("r", hasDefault)]]
+  , consModuleWithSuite intersect [[("r", hasDefault)]]
+  , consModuleWithSuite difference [[("r", hasDefault)]]
+  , consModuleWithSuite translate [[("x", noDefault), ("y", noDefault), ("z", noDefault)], [("v", noDefault)]]
+  , consModuleWithSuite rotate [[("a", noDefault), ("v", hasDefault)]]
+  , consModuleWithSuite scale [[("v", noDefault)]]
+  , consModuleWithSuite extrude [[("height", hasDefault), ("center", hasDefault), ("twist", hasDefault), ("scale", hasDefault), ("translate", hasDefault), ("r", hasDefault)]]
+  , consModuleWithSuite rotateExtrude [[("angle", hasDefault), ("r", hasDefault), ("translate", hasDefault), ("rotate", hasDefault)]]
+  , consModuleWithSuite shell [[("w", noDefault)]]
+  , consModuleWithSuite projection [[("cut", hasDefault)]]
+  , consModuleWithSuite pack [[("size", noDefault), ("sep", noDefault)]]
+  , consModuleWithSuite unit [[("unit", noDefault)]]
+  , consModuleWithSuite mirror [[("x", noDefault), ("y", noDefault), ("z", noDefault)], [("v", noDefault)]]
+  , consModuleWithSuite multmatrix [[("m", noDefault)]]
   ]
   where
     hasDefault = True
     noDefault = False
-    noSuite :: Maybe Bool
-    noSuite = Nothing
-    requiredSuite = Just True
-    onModIze func rawInstances = (name, ONModule name implementation instances)
+    consModuleWithSuite :: (Symbol, SourcePosition -> [OVal] -> ArgParser (StateC [OVal])) -> [[(Text, Bool)]] -> (Symbol, OVal)
+    consModuleWithSuite func rawInstances = (name, ONModuleWithSuite name implementation instances)
       where
         (name, implementation) = func
         instances = fmap fixup rawInstances
-        fixup :: ([(Text, Bool)], Maybe Bool) -> ([(Symbol, Bool)], Maybe Bool)
-        fixup (args, suiteInfo) = (fmap fixupArgs args, suiteInfo)
+        fixup :: [(Text, Bool)] -> [(Symbol, Bool)]
+        fixup args = fmap fixupArgs args
+          where
+            fixupArgs :: (Text, Bool) -> (Symbol, Bool)
+            fixupArgs (symbol, maybeDefault) = (Symbol symbol, maybeDefault)
+    consModule :: (Symbol, SourcePosition -> ArgParser (StateC [OVal])) -> [[(Text, Bool)]] -> (Symbol, OVal)
+    consModule func rawInstances = (name, ONModule name implementation instances)
+      where
+        (name, implementation) = func
+        instances = fmap fixup rawInstances
+        fixup :: [(Text, Bool)] -> [(Symbol, Bool)]
+        fixup (args) = fmap fixupArgs args
           where
             fixupArgs :: (Text, Bool) -> (Symbol, Bool)
             fixupArgs (symbol, maybeDefault) = (Symbol symbol, maybeDefault)
@@ -97,8 +104,8 @@ primitiveModules =
 -- | sphere is a module without a suite.
 --   this means that the parser will look for this like
 --   sphere(args...);
-sphere :: (Symbol, SourcePosition -> [OVal] -> ArgParser (StateC [OVal]))
-sphere = moduleWithoutSuite "sphere" $ \_ _ -> do
+sphere :: (Symbol, SourcePosition -> ArgParser (StateC [OVal]))
+sphere = moduleWithoutSuite "sphere" $ \_ -> do
     example "sphere(3);"
     example "sphere(r=5);"
     -- arguments:
@@ -120,8 +127,8 @@ sphere = moduleWithoutSuite "sphere" $ \_ _ -> do
 
 -- | FIXME: square1, square2 like cylinder has?
 --   FIXME: translate for square2?
-cube :: (Symbol, SourcePosition -> [OVal] -> ArgParser (StateC [OVal]))
-cube = moduleWithoutSuite "cube" $ \_ _ -> do
+cube :: (Symbol, SourcePosition -> ArgParser (StateC [OVal]))
+cube = moduleWithoutSuite "cube" $ \_ -> do
     -- examples
     example "cube(size = [2,3,4], center = true, r = 0.5);"
     example "cube(4);"
@@ -165,8 +172,8 @@ cube = moduleWithoutSuite "cube" $ \_ _ -> do
     -- Implementation
     addObj3 $ Prim.withRounding r $ Prim.rect3 (V3 x1 y1 z1) (V3 x2 y2 z2)
 
-square :: (Symbol, SourcePosition -> [OVal] -> ArgParser (StateC [OVal]))
-square = moduleWithoutSuite "square" $ \_ _ -> do
+square :: (Symbol, SourcePosition -> ArgParser (StateC [OVal]))
+square = moduleWithoutSuite "square" $ \_ -> do
     -- examples
     example "square(x=[-2,2], y=[-1,5]);"
     example "square(size = [3,4], center = true, r = 0.5);"
@@ -206,8 +213,8 @@ square = moduleWithoutSuite "square" $ \_ _ -> do
     -- Implementation
     addObj2 $ Prim.withRounding r $ Prim.rect (V2 x1 y1) (V2 x2 y2)
 
-cylinder :: (Symbol, SourcePosition -> [OVal] -> ArgParser (StateC [OVal]))
-cylinder = moduleWithoutSuite "cylinder" $ \_ _ -> do
+cylinder :: (Symbol, SourcePosition -> ArgParser (StateC [OVal]))
+cylinder = moduleWithoutSuite "cylinder" $ \_ -> do
     example "cylinder(r=10, h=30, center=true);"
     example "cylinder(r1=4, r2=6, h=10);"
     example "cylinder(r=5, h=10, $fn = 6);"
@@ -268,27 +275,31 @@ cylinder = moduleWithoutSuite "cylinder" $ \_ _ -> do
         in shift obj3
         else shift $ Prim.cylinder2 r1 r2 dh
 
-polyhedron :: (Symbol, SourcePosition -> [OVal] -> ArgParser (StateC [OVal]))
-polyhedron = moduleWithoutSuite "polyhedron" $ \_ _ -> do
+polyhedron :: (Symbol, SourcePosition -> ArgParser (StateC [OVal]))
+polyhedron = moduleWithoutSuite "polyhedron" $ \sourcePos -> do
     example "polyhedron(points=[[0,0,0], [2,0,0], [2,2,0], [0,2,0], [1, 1, 2]], faces=[[0,1,2,3], [0,5,1], [1,5,2], [2,5,3], [3,5,4], [4,5,0]]);"
     -- arguments
     points :: [ℝ3] <- argument "points" `defaultTo` [] `doc` "list of points to construct faces from"
     faces :: [[ℕ]] <- argument "faces" `defaultTo` [] `doc` "list of sets of indices into points, used to create faces on the polyhedron."
-    -- A tri is constructed of three indexes into the points.
-    -- This decomposes our faces into tris.
-    let
-      tris = concatMap trianglesFromFace faces
-      in
-      addObj3 $ Prim.polyhedron points tris
+    pure $ do
+      -- A tri is constructed of three indexes into the points.
+      -- This decomposes our faces into tris.
+      tris <- fmap concat $ mapM (trianglesFromFace sourcePos) faces
+      pure [OObj3 $ Prim.polyhedron points tris]
       where
-        -- FIXME: use warnC here, instead of error.
-        trianglesFromFace :: [ℕ] -> [(ℕ,ℕ,ℕ)]
-        trianglesFromFace [p1,p2,p3] = [(p1,p2,p3)]
-        trianglesFromFace (p1:p2:p3:xs) = ((p1,p2,p3):trianglesFromFace (p1:p3:xs))
-        trianglesFromFace ys = error $ "too few points:" <> (show ys) <> "\n"
-
-cone :: (Symbol, SourcePosition -> [OVal] -> ArgParser (StateC [OVal]))
-cone = moduleWithoutSuite "cone" $ \_ _ -> do
+        trianglesFromFace :: SourcePosition -> [ℕ] -> StateC [(ℕ,ℕ,ℕ)]
+        trianglesFromFace _         []         = pure []
+        trianglesFromFace sourcePos [p1]       = do
+                                                 warnC sourcePos $ "only one point found: " <> (DTL.pack $ show p1) <> "\n"
+                                                 pure []
+        trianglesFromFace sourcePos [p1,p2]    = do
+                                                 warnC sourcePos $ "only two points found: " <> (DTL.pack $ show p1) <> "\n" <> (DTL.pack $ show p2) <> "\n"
+                                                 pure []
+        trianglesFromFace _         [p1,p2,p3] = pure [(p1,p2,p3)]
+        trianglesFromFace sourcePos (p1:p2:p3:xs) = ((p1,p2,p3):) <$> trianglesFromFace sourcePos (p1:p3:xs)
+
+cone :: (Symbol, SourcePosition -> ArgParser (StateC [OVal]))
+cone = moduleWithoutSuite "cone" $ \_ -> do
     example "cone(r=10, h=30, center=true);"
     -- arguments
     r <- do
@@ -318,8 +329,8 @@ cone = moduleWithoutSuite "cone" $ \_ _ -> do
             else Prim.translate (V3 0 0 h1)
     addObj3 . shift $ Prim.cone r dh
 
-torus :: (Symbol, SourcePosition -> [OVal] -> ArgParser (StateC [OVal]))
-torus = moduleWithoutSuite "torus" $ \_ _ -> do
+torus :: (Symbol, SourcePosition -> ArgParser (StateC [OVal]))
+torus = moduleWithoutSuite "torus" $ \_ -> do
     example "torus(r1=10, r2=5);"
     -- arguments
     (r1, r2) <- (,)
@@ -332,8 +343,8 @@ torus = moduleWithoutSuite "torus" $ \_ _ -> do
     -- based on the args.
     addObj3 $ Prim.torus r1 r2
 
-ellipsoid :: (Symbol, SourcePosition -> [OVal] -> ArgParser (StateC [OVal]))
-ellipsoid = moduleWithoutSuite "ellipsoid" $ \_ _ -> do
+ellipsoid :: (Symbol, SourcePosition -> ArgParser (StateC [OVal]))
+ellipsoid = moduleWithoutSuite "ellipsoid" $ \_ -> do
     example "ellipsoid(a=1, b=2, c=3);"
     -- arguments
     (a, b, c) <- (,,)
@@ -347,8 +358,8 @@ ellipsoid = moduleWithoutSuite "ellipsoid" $ \_ _ -> do
     -- based on the args.
     addObj3 $ Prim.ellipsoid a b c
 
-circle :: (Symbol, SourcePosition -> [OVal] -> ArgParser (StateC [OVal]))
-circle = moduleWithoutSuite "circle" $ \_ _ -> do
+circle :: (Symbol, SourcePosition -> ArgParser (StateC [OVal]))
+circle = moduleWithoutSuite "circle" $ \_ -> do
     example "circle(r=10); // circle"
     example "circle(r=5, $fn=6); //hexagon"
     -- Arguments
@@ -376,8 +387,8 @@ circle = moduleWithoutSuite "circle" $ \_ _ -> do
 -- | FIXME: handle rectangles that are not grid alligned.
 --   FIXME: allow for rounding of polygon corners, specification of vertex ordering.
 --   FIXME: polygons have to have more than two points, or do not generate geometry, and generate an error.
-polygon :: (Symbol, SourcePosition -> [OVal] -> ArgParser (StateC [OVal]))
-polygon = moduleWithoutSuite "polygon" $ \_ _ -> do
+polygon :: (Symbol, SourcePosition -> ArgParser (StateC [OVal]))
+polygon = moduleWithoutSuite "polygon" $ \_ -> do
     example "polygon ([(0,0), (0,10), (10,0)]);"
     points :: [ℝ2]  <- argument "points"
                         `doc` "vertices of the polygon"
@@ -433,22 +444,26 @@ intersect = moduleWithSuite "intersection" $ \_ children -> do
         else objReduce  Prim.intersect      Prim.intersect     children
 
 difference :: (Symbol, SourcePosition -> [OVal] -> ArgParser (StateC [OVal]))
-difference = moduleWithSuite "difference" $ \_ children -> do
-    when (null children) $ APFail "Call to 'difference' requires at least one child"
+difference = moduleWithSuite "difference" $ \sourcePos children -> do
     r :: ℝ <- argument "r"
         `defaultTo` 0
         `doc` "Radius of rounding for the difference interface"
-    pure $ pure $ if r > 0
-        then objReduce (unsafeUncurry (Prim.differenceR r)) (unsafeUncurry (Prim.differenceR r)) children
-        else objReduce (unsafeUncurry  Prim.difference)     (unsafeUncurry  Prim.difference)     children
-  where
-    unsafeUncurry :: (a -> [a] -> c) -> [a] -> c
-    unsafeUncurry f = uncurry f . unsafeUncons
-
-    unsafeUncons :: [a] -> (a, [a])
-    unsafeUncons (a : as) = (a, as)
-    -- NOTE: This error is guarded against during the @null children@ check in the function body.
-    unsafeUncons _ = error "difference requires at least one element; zero given"
+    pure $ do
+      if (null children)
+        then do
+             errorC sourcePos "difference requires at least one element; none given."
+             pure []
+        else pure $ if r > 0
+                    then objReduce (unsafeUncurry (Prim.differenceR r)) (unsafeUncurry (Prim.differenceR r)) children
+                    else objReduce (unsafeUncurry  Prim.difference)     (unsafeUncurry  Prim.difference)     children
+        where
+          unsafeUncurry :: (a -> [a] -> c) -> [a] -> c
+          unsafeUncurry f = uncurry f . unsafeUncons
+
+          unsafeUncons :: [a] -> (a, [a])
+          unsafeUncons (a : as) = (a, as)
+          -- NOTE: This error is guarded against during the @null children@ check in the function body.
+          unsafeUncons _ = error "difference requires at least one element; zero given"
 
 translate :: (Symbol, SourcePosition -> [OVal] -> ArgParser (StateC [OVal]))
 translate = moduleWithSuite "translate" $ \_ children -> do
@@ -711,7 +726,7 @@ multmatrix = moduleWithSuite "multmatrix" $ \_ children -> do
 moduleWithSuite :: Text -> (SourcePosition -> [OVal] -> ArgParser (StateC [OVal])) -> (Symbol, SourcePosition -> [OVal] -> ArgParser (StateC [OVal]))
 moduleWithSuite name modArgMapper = (Symbol name, modArgMapper)
 
-moduleWithoutSuite :: Text -> (SourcePosition -> [OVal] -> ArgParser (StateC [OVal])) -> (Symbol, SourcePosition -> [OVal] -> ArgParser (StateC [OVal]))
+moduleWithoutSuite :: Text -> (SourcePosition -> ArgParser (StateC [OVal])) -> (Symbol, SourcePosition -> ArgParser (StateC [OVal]))
 moduleWithoutSuite name modArgMapper = (Symbol name, modArgMapper)
 
 addObj2 :: SymbolicObj2 -> ArgParser (StateC [OVal])
diff --git a/Graphics/Implicit/ExtOpenScad/Util/OVal.hs b/Graphics/Implicit/ExtOpenScad/Util/OVal.hs
index 657670d5..f7ca0e8f 100644
--- a/Graphics/Implicit/ExtOpenScad/Util/OVal.hs
+++ b/Graphics/Implicit/ExtOpenScad/Util/OVal.hs
@@ -18,7 +18,7 @@ import Prelude(Maybe(Just, Nothing), Bool(True, False), Either(Left,Right), (==)
 
 import Graphics.Implicit.Definitions(V2, ℝ, ℝ2, ℕ, SymbolicObj2, SymbolicObj3, ExtrudeMScale(C1, C2, Fn), fromℕtoℝ)
 
-import Graphics.Implicit.ExtOpenScad.Definitions (OVal(ONum, OBool, OString, OList, OFunc, OUndefined, OUModule, ONModule, OVargsModule, OError, OObj2, OObj3, OIO))
+import Graphics.Implicit.ExtOpenScad.Definitions (OVal(ONum, OBool, OString, OList, OFunc, OUndefined, OUModule, ONModule, ONModuleWithSuite, OVargsModule, OError, OObj2, OObj3, OIO))
 
 import Control.Monad (msum)
 
@@ -168,7 +168,8 @@ oTypeStr (OString        _ ) = "String"
 oTypeStr (OFunc          _ ) = "Function"
 oTypeStr (OIO            _ ) = "IO"
 oTypeStr (OUModule {}      ) = "User Defined Module"
-oTypeStr (ONModule {}      ) = "Built-in Module"
+oTypeStr (ONModuleWithSuite {}) = "Built-in Module with suite"
+oTypeStr (ONModule {}) = "Built-in Module"
 oTypeStr (OVargsModule _ _ ) = "VargsModule"
 oTypeStr (OError         _ ) = "Error"
 oTypeStr (OObj2          _ ) = "2D Object"
diff --git a/Graphics/Implicit/ObjectUtil/GetImplicit3.hs b/Graphics/Implicit/ObjectUtil/GetImplicit3.hs
index a9b3c1b5..91e204fa 100644
--- a/Graphics/Implicit/ObjectUtil/GetImplicit3.hs
+++ b/Graphics/Implicit/ObjectUtil/GetImplicit3.hs
@@ -245,8 +245,8 @@ insideMesh triangles point = odd . length . filter (==True) $ foundIntersection
       | u + v > 1-eps-eps     = False
       | otherwise             = t > 0
       where
-        -- fudge factor.
         -- adjusted from: 6,7,9,12,15,17,18==(43,171), 16 ==(38,171)
+        -- fudge factor.
         eps :: ℝ
         eps = 1e-16
         -- Our edge vectors. We have picked v1 to address the space by, for convenience.
diff --git a/programs/docgen.hs b/programs/docgen.hs
index 4f244422..8eefca75 100644
--- a/programs/docgen.hs
+++ b/programs/docgen.hs
@@ -6,7 +6,7 @@
 
 import Prelude(IO, Show, String, Int, Maybe(Just,Nothing), Eq, return, ($), show, fmap, (<>), putStrLn, filter, zip, null, undefined, const, Bool(True,False), fst, (.), head, length, (/=), (+), error, print, drop)
 import Graphics.Implicit.ExtOpenScad.Primitives (primitiveModules)
-import Graphics.Implicit.ExtOpenScad.Definitions (ArgParser(AP,APFail,APExample,APTest,APTerminator,APBranch), Symbol(Symbol), OVal(ONModule), SourcePosition(SourcePosition), StateC)
+import Graphics.Implicit.ExtOpenScad.Definitions (ArgParser(AP,APFail,APExample,APTest,APTerminator,APBranch), Symbol(Symbol), OVal(ONModule, ONModuleWithSuite), SourcePosition(SourcePosition), StateC)
 
 import qualified Control.Exception as Ex (catch, SomeException)
 import Control.Monad (forM_)
@@ -114,9 +114,10 @@ main = do
           dumpPrimitive moduleName moduleDocList 0
           where
             getArgParserFrom :: (Symbol, OVal) -> ArgParser(StateC [OVal])
-            getArgParserFrom (_, ONModule _ implementation _) = implementation sourcePosition []
-              where sourcePosition = SourcePosition 0 0 "docgen"
+            getArgParserFrom (_, ONModule _ implementation _) = implementation sourcePosition
+            getArgParserFrom (_, ONModuleWithSuite _ implementation _) = implementation sourcePosition []
             getArgParserFrom (_, _) = error "bad value in primitive array."
+            sourcePosition = SourcePosition 0 0 "docgen"
 
 -- | the format we extract documentation into
 data Doc = Doc String [DocPart]

From 60ff96b4226fedbc3b01b5ded652a6a9098a61d2 Mon Sep 17 00:00:00 2001
From: Julia Longtin <julia.longtin@gmail.com>
Date: Sat, 7 Mar 2026 15:10:17 +0000
Subject: [PATCH 11/28] reformat for clarity.

---
 Graphics/Implicit/ExtOpenScad/Util/OVal.hs | 26 +++++++++++-----------
 1 file changed, 13 insertions(+), 13 deletions(-)

diff --git a/Graphics/Implicit/ExtOpenScad/Util/OVal.hs b/Graphics/Implicit/ExtOpenScad/Util/OVal.hs
index f7ca0e8f..2f42dadb 100644
--- a/Graphics/Implicit/ExtOpenScad/Util/OVal.hs
+++ b/Graphics/Implicit/ExtOpenScad/Util/OVal.hs
@@ -160,20 +160,20 @@ instance OTypeMirror ExtrudeMScale where
 
 -- A string representing each type.
 oTypeStr :: OVal -> Text
-oTypeStr OUndefined         = "Undefined"
-oTypeStr (OBool          _ ) = "Bool"
-oTypeStr (ONum           _ ) = "Number"
-oTypeStr (OList          _ ) = "List"
-oTypeStr (OString        _ ) = "String"
-oTypeStr (OFunc          _ ) = "Function"
-oTypeStr (OIO            _ ) = "IO"
-oTypeStr (OUModule {}      ) = "User Defined Module"
+oTypeStr OUndefined             = "Undefined"
+oTypeStr (OBool             _ ) = "Bool"
+oTypeStr (ONum              _ ) = "Number"
+oTypeStr (OList             _ ) = "List"
+oTypeStr (OString           _ ) = "String"
+oTypeStr (OFunc             _ ) = "Function"
+oTypeStr (OIO               _ ) = "IO"
+oTypeStr (OUModule {}         ) = "User Defined Module"
 oTypeStr (ONModuleWithSuite {}) = "Built-in Module with suite"
-oTypeStr (ONModule {}) = "Built-in Module"
-oTypeStr (OVargsModule _ _ ) = "VargsModule"
-oTypeStr (OError         _ ) = "Error"
-oTypeStr (OObj2          _ ) = "2D Object"
-oTypeStr (OObj3          _ ) = "3D Object"
+oTypeStr (ONModule {}         ) = "Built-in Module"
+oTypeStr (OVargsModule    _ _ ) = "VargsModule"
+oTypeStr (OError            _ ) = "Error"
+oTypeStr (OObj2             _ ) = "2D Object"
+oTypeStr (OObj3             _ ) = "3D Object"
 
 getErrors :: OVal -> Maybe Text
 getErrors (OError er) = Just er

From c2109a9cee16acd31c8ba33aeadc37684260bd83 Mon Sep 17 00:00:00 2001
From: Julia Longtin <julia.longtin@gmail.com>
Date: Tue, 10 Mar 2026 01:02:25 +0000
Subject: [PATCH 12/28] much closer. some shadowing to deal with.

---
 Graphics/Implicit/ExtOpenScad/Primitives.hs  |  96 +++++++++-
 Graphics/Implicit/ObjectUtil/GetImplicit3.hs | 178 ++++++++-----------
 2 files changed, 168 insertions(+), 106 deletions(-)

diff --git a/Graphics/Implicit/ExtOpenScad/Primitives.hs b/Graphics/Implicit/ExtOpenScad/Primitives.hs
index 4956b46f..e153a6a0 100644
--- a/Graphics/Implicit/ExtOpenScad/Primitives.hs
+++ b/Graphics/Implicit/ExtOpenScad/Primitives.hs
@@ -15,10 +15,12 @@
 -- Export one set containing all of the primitive modules.
 module Graphics.Implicit.ExtOpenScad.Primitives (primitiveModules) where
 
-import Prelude(concat, mapM, (.), Either(Left, Right), Bool(True, False), Maybe(Just, Nothing), ($), pure, show, either, id, (-), (==), (&&), (<), (*), cos, sin, pi, (/), (>), const, uncurry, (/=), (||), not, null, fmap, (<>), otherwise, error, (<*>), (<$>))
+import Prelude(any, concat, elem, foldr, head, mapM, (.), Either(Left, Right), Bool(True, False), Maybe(Just, Nothing), ($), pure, show, either, id, (-), (==), (&&), (<), (*), cos, sin, pi, (/), (>), const, uncurry, (/=), (||), not, null, fmap, (<>), otherwise, error, (<*>), (<$>))
 
 import Graphics.Implicit.Definitions (ℝ, ℝ2, ℝ3, ℕ, SymbolicObj2, SymbolicObj3, ExtrudeMScale(C1), fromℕtoℝ, isScaleID)
 
+import Graphics.Implicit.Export.Util (centroid)
+
 import Graphics.Implicit.ExtOpenScad.Definitions (OVal (OObj2, OObj3, ONModule, ONModuleWithSuite), ArgParser, Symbol(Symbol), StateC, SourcePosition)
 
 import Graphics.Implicit.ExtOpenScad.Util.ArgParser (doc, defaultTo, example, test, eulerCharacteristic)
@@ -32,13 +34,24 @@ import Graphics.Implicit.ExtOpenScad.Util.StateC (errorC, warnC)
 -- Note the use of a qualified import, so we don't have the functions in this file conflict with what we're importing.
 import qualified Graphics.Implicit.Primitives as Prim (withRounding, sphere, rect3, rect, translate, circle, polygon, polyhedron, extrude, cylinder2, union, unionR, intersect, intersectR, difference, differenceR, rotate, slice, transform, rotate3V, rotate3, transform3, scale, extrudeM, rotateExtrude, shell, mirror, pack3, pack2, torus, ellipsoid, cone)
 
-import Control.Monad (mplus)
+import Control.Monad (foldM, mplus)
+
+import Data.Foldable (toList)
+
+import Data.List (genericIndex)
+
+import Data.Sequence (Seq, deleteAt, filter, fromList)
+import qualified Data.Sequence as DS (null)
 
 import Data.Text.Lazy (Text)
 import qualified Data.Text.Lazy as DTL (pack)
 
+import Data.Maybe (isJust)
+
 import Control.Lens ((^.))
-import Linear (_m33, M34, M44, V2(V2), V3(V3), V4(V4))
+
+import Linear (_m33, cross, dot, M34, M44, V2(V2), V3(V3), V4(V4))
+
 import Linear.Affine (qdA)
 
 default (ℝ)
@@ -275,6 +288,9 @@ cylinder = moduleWithoutSuite "cylinder" $ \_ -> do
         in shift obj3
         else shift $ Prim.cylinder2 r1 r2 dh
 
+-- A triangle, with the points in an external index.
+type Tri = (ℕ,ℕ,ℕ)
+
 polyhedron :: (Symbol, SourcePosition -> ArgParser (StateC [OVal]))
 polyhedron = moduleWithoutSuite "polyhedron" $ \sourcePos -> do
     example "polyhedron(points=[[0,0,0], [2,0,0], [2,2,0], [0,2,0], [1, 1, 2]], faces=[[0,1,2,3], [0,5,1], [1,5,2], [2,5,3], [3,5,4], [4,5,0]]);"
@@ -283,12 +299,15 @@ polyhedron = moduleWithoutSuite "polyhedron" $ \sourcePos -> do
     faces :: [[ℕ]] <- argument "faces" `defaultTo` [] `doc` "list of sets of indices into points, used to create faces on the polyhedron."
     pure $ do
       -- A tri is constructed of three indexes into the points.
-      -- This decomposes our faces into tris.
       tris <- fmap concat $ mapM (trianglesFromFace sourcePos) faces
-      pure [OObj3 $ Prim.polyhedron points tris]
+      woundTris <- reWindTriangles sourcePos points tris
+      pure [OObj3 $ Prim.polyhedron points woundTris]
       where
-        trianglesFromFace :: SourcePosition -> [ℕ] -> StateC [(ℕ,ℕ,ℕ)]
-        trianglesFromFace _         []         = pure []
+        -- decompose our faces into tris.
+        trianglesFromFace :: SourcePosition -> [ℕ] -> StateC [Tri]
+        trianglesFromFace sourcePos []         = do
+                                                 errorC sourcePos "no point found when trying to generate triangles from a face.\n"
+                                                 pure []
         trianglesFromFace sourcePos [p1]       = do
                                                  warnC sourcePos $ "only one point found: " <> (DTL.pack $ show p1) <> "\n"
                                                  pure []
@@ -297,6 +316,69 @@ polyhedron = moduleWithoutSuite "polyhedron" $ \sourcePos -> do
                                                  pure []
         trianglesFromFace _         [p1,p2,p3] = pure [(p1,p2,p3)]
         trianglesFromFace sourcePos (p1:p2:p3:xs) = ((p1,p2,p3):) <$> trianglesFromFace sourcePos (p1:p3:xs)
+        -- | Ensure our triangles are wound in the same direction
+        reWindTriangles :: SourcePosition -> [ℝ3] -> [Tri] -> StateC [Tri]
+        reWindTriangles _         _ [] = pure []
+        -- Really, forces them to have the same winding as the first triangle, from us putting [tri] here.
+        reWindTriangles sourcePos points (tri:moreTris) = windTriangles [safeTri] (fromList moreTris)
+          where
+            polyCentroid = centroid points
+            -- The first triangle, flipped based on comparing two centroids.
+            safeTri
+              | (triCentroid - polyCentroid) `dot` triNorm < 0 = flipTri tri
+              | otherwise                                      = tri
+              where
+                (p1,p2,p3) = tri
+                (v1,v2,v3) = (genericIndex points p1,genericIndex points p2,genericIndex points p3)
+                -- The norm of the safe triangle.
+                triNorm    = (v2-v1) `cross` (v3-v1)
+                triCentroid = centroid [v1,v2,v3]
+                -- type conversion.
+            flipTri (p1,p2,p3) = (p1,p3,p2)
+            -- | Wind the triangles.
+            windTriangles :: [Tri] -> Seq Tri -> StateC [Tri]
+            windTriangles visited unvisited
+              -- We're done.
+              | DS.null unvisited = pure visited
+              | otherwise = do
+                  (newVisited, newUnvisited) <- foldM (classifyTri visited) ([], unvisited) (toList unvisited)
+                  if null newVisited
+                    then do
+                         warnC sourcePos $ "Had to pick a new root"
+                         windTriangles (visited <> [head $ toList newUnvisited]) (deleteAt 0 newUnvisited)
+                    else windTriangles (visited <> newVisited)                   newUnvisited
+            -- | Compare one unvisited triangle against all visited triangles.
+            -- If it's a neighbor of a visited triangle, correct it's winding and throw a warning if we had to flip it, then move it from unvisited to found (after flipping).
+            classifyTri :: [Tri] -> ([Tri], Seq Tri) -> Tri -> StateC ([Tri], Seq Tri)
+            classifyTri visited (found, remaining) triUnderTest =
+              case res of
+                Just triFound -> do -- See if we flipped our tri, and if so, throw a warning.
+                                 if triFound /= triUnderTest
+                                   then warnC sourcePos $ "Flipped face detected with vertices " <> (DTL.pack $ show triUnderTest)
+                                   else pure ()
+                                 pure (found <> [triFound], filter (/= triUnderTest) remaining)
+                Nothing ->       pure (found, remaining)
+              where
+                res = foldr (\tri state -> firstNeighborFilter tri triUnderTest state) Nothing visited 
+                -- | A short-circuiting filter we fold over visited, and grab the first neighboring tri.
+                firstNeighborFilter :: Tri -> Tri -> Maybe Tri -> Maybe Tri
+                firstNeighborFilter src triUnderTest res
+                  | isJust res = res
+                  | otherwise  = maybeWindNeighbor src triUnderTest
+                -- | Checks whether a triangle under test is a neighbor of the given triangle, and if it is, returns if after ensuring it is wound in the proper direction.
+                maybeWindNeighbor :: Tri -> Tri -> Maybe Tri
+                maybeWindNeighbor src triUnderTest
+                -- A correctly wound neighbor will have the opposite direction, when it is referring to a given edge.
+                  | oppositeNeighbor = Just           triUnderTest
+                -- A neighbor that needs flipped will refer to an edge in the same direction as our given triangle.
+                  | sameNeighbor     = Just $ flipTri triUnderTest
+                  | otherwise = Nothing
+                  where
+                    oppositeNeighbor = any (\edge -> flip edge `elem` edgesOfTri src) $ edgesOfTri triUnderTest
+                      where
+                        flip (a,b) = (b,a)
+                    sameNeighbor     = any (\edge ->      edge `elem` edgesOfTri src) $ edgesOfTri triUnderTest
+                    edgesOfTri (p1,p2,p3) = [(p1,p2), (p2,p3), (p3,p1)]
 
 cone :: (Symbol, SourcePosition -> ArgParser (StateC [OVal]))
 cone = moduleWithoutSuite "cone" $ \_ -> do
diff --git a/Graphics/Implicit/ObjectUtil/GetImplicit3.hs b/Graphics/Implicit/ObjectUtil/GetImplicit3.hs
index 91e204fa..4f2b4954 100644
--- a/Graphics/Implicit/ObjectUtil/GetImplicit3.hs
+++ b/Graphics/Implicit/ObjectUtil/GetImplicit3.hs
@@ -3,16 +3,18 @@
 -- Implicit CAD. Copyright (C) 2011, Christopher Olah (chris@colah.ca)
 -- Released under the GNU AGPLV3+, see LICENSE
 
-module Graphics.Implicit.ObjectUtil.GetImplicit3 (getImplicit3) where
+module Graphics.Implicit.ObjectUtil.GetImplicit3 (getImplicit3, distancePointToTriangle, closestPointToTriangle) where
 
 -- Import only what we need from the prelude.
-import Prelude (abs, atan2, cos, ceiling, either, error, floor, fromInteger, id, max, min, minimum, negate, odd, otherwise, pi, pure, round, show, sin, sqrt, (||), (/=), Either(Left, Right), (<), (<=), (<>), (>), (>=), (&&), (-), (/), (*), (+), ($), (.), Bool(True, False), (==), (**), Num, Applicative, (<$>))
+import Prelude (abs, atan2, cos, ceiling, either, error, floor, fromInteger, id, max, min, minimum, negate, otherwise, pi, pure, round, show, sin, sqrt, (||), (/=), Either(Left, Right), (<), (<=), (<>), (>), (>=), (&&), (-), (/), (*), (+), ($), (.), Bool(True, False), (==), (**), Num, Applicative, (<$>))
 
 import Control.Lens ((^.))
 
 import qualified Data.Either as Either (either)
 
-import Data.List(filter, genericIndex, length)
+import Data.List (genericIndex, length, minimumBy)
+
+import Data.Ord (compare)
 
 import Linear (V2(V2), V3(V3), _xy, _z, distance, dot)
 import qualified Linear (conjugate, inv44, normalizePoint, normalize, point, rotate, Metric)
@@ -77,10 +79,24 @@ getImplicit3 _ (Cylinder h r1 r2) = \(V3 x y z) ->
 getImplicit3 _ (Polyhedron [] _) = error "Asked to find distance to an empty polygon. No points."
 getImplicit3 _ (Polyhedron _ []) = error "Asked to find distance to an empty polygon. No tris."
 getImplicit3 _ (Polyhedron points tris) = \(point) ->
-  if insideMesh triangles point
-  then negate $ minimum $ distancePointToTriangle point <$> triangles
-  else          minimum $ distancePointToTriangle point <$> triangles
+  let
+    (res, closestTri) = unsignedDistanceAndTriangleClosestTo point
+  in
+    if pointOnOutside point closestTri
+    then          res
+    else negate $ res
   where
+    unsignedDistanceAndTriangleClosestTo point = minimumBy (\(a,_) (b,_) -> a `compare` b) $ distTris point
+    distTris point = (\a -> (distancePointToTriangle point a, a)) <$> triangles
+    -- NOTE: May get slightly different values, depending on the V selected.
+    normOfTri :: (V3 ℝ,V3 ℝ,V3 ℝ) -> V3 ℝ
+    normOfTri (v1,v2,v3) = Linear.normalize $ (v2-v1) `cross` (v3-v1)
+    firstPointOfTri (v1,_,_) = v1
+    pointOnOutside point closestTri = (point - firstPointOfTri closestTri) `dot` normOfTri closestTri >= -eps
+      where
+        -- fudge factor.
+        eps :: ℝ
+        eps = 1e-13
     -- decompose our tris into triangles.
     triangles = findTriangle points <$> tris
     -- FIXME: make these indices correct by construction?
@@ -94,7 +110,7 @@ getImplicit3 _ (Polyhedron points tris) = \(point) ->
         where
           -- FIXME: >=BASE-4.21: replace this with compareLength once debian stable ships 4.21.
           outOfRange :: ℕ -> Bool
-          outOfRange v = v < 0 || length virtices < fromℕ v
+          outOfRange v = v < 0 || length virtices <= fromℕ v
 getImplicit3 _ (BoxFrame b e) = \p' ->
     let p@(V3 px py pz) = abs p' - b
         V3 qx qy qz = abs (p + pure e) - pure e
@@ -227,53 +243,14 @@ getImplicit3 ctx (RotateExtrude totalRotation translate rotate symbObj) =
               else obj rz_pos
 getImplicit3 ctx (Shared3 obj) = getImplicitShared ctx obj
 
--- | Check to see if a point is inside or outside of a mesh.
--- FIXME: Replace this with a winding number based approach. See: https://github.com/marmakoide/inside-3d-mesh/blob/master/README.md
-insideMesh :: [(V3 ℝ, V3 ℝ, V3 ℝ)] -> V3 ℝ -> Bool
-insideMesh triangles point = odd . length . filter (==True) $ foundIntersection point unsafeSafeRay <$> triangles
-  where
-    -- our assumed safe ray. it's not really safe, but.. gotta start somewhere.
-    unsafeSafeRay :: V3 ℝ
-    unsafeSafeRay = Linear.normalize $ V3 (sqrt 2) (sqrt 3) (sqrt 5)
-    foundIntersection :: V3 ℝ -> V3 ℝ -> (V3 ℝ, V3 ℝ, V3 ℝ) -> Bool
-    foundIntersection source direction (v1, v2, v3)
-      -- shouldn't happen, triangle is on the same plane as unsafeSafeRay
-      | abs determinate < eps = False
-      -- the intersection is not on the triangle (distance along vec12 out of range)
-      | u < eps || u > 1-eps  = False
-      | v < eps || v > 1-eps  = False
-      | u + v > 1-eps-eps     = False
-      | otherwise             = t > 0
-      where
-        -- adjusted from: 6,7,9,12,15,17,18==(43,171), 16 ==(38,171)
-        -- fudge factor.
-        eps :: ℝ
-        eps = 1e-16
-        -- Our edge vectors. We have picked v1 to address the space by, for convenience.
-        vec12 = v2 - v1
-        vec13 = v3 - v1
-        -- The normal of the plane formed by vec13 and our 'safe' vector. at 90 degree angles to both.
-        norm13d = direction `cross` vec13
-        -- The determinate. A volume. AKA, how close to parallel are the triangle, and the plane with normal norm13d.
-        determinate = vec12 `dot` norm13d
-        -- An edge vector in our 'safe' direction, from v1.
-        vec1s = source - v1
-        -- The distance along vec12 to find the intersection of a ray from source in direction direction, and the plane our triangle is on.
-        u = (vec1s `dot` norm13d) / determinate
-        -- The normal of the plane fromed by vec12 and vec1s.
-        norm121s = vec1s `cross` vec12
-        -- The distance along vec13 to find the intersection of a ray from source in direction direction, and the plane our triangle is on.
-        v = (direction `dot` norm121s) / determinate
-        -- The distance along the ray to find the intersection of a ray from source in direction direction, and the plane our triangle is on.
-        t = (vec13 `dot` norm121s) / determinate
-
+-- | You see, what I thought I'd do is put a raytracer inside of a raytracer... what could go wrong...
 -- With inspiration from: https://github.com/RenderKit/embree/blob/master/tutorials/common/math/closest_point.h
 distancePointToTriangle :: V3 ℝ -> (V3 ℝ,V3 ℝ,V3 ℝ) -> ℝ
-distancePointToTriangle point (virtex1, virtex2, virtex3) = distance point closestPointToTriangleCenteredSorted
+distancePointToTriangle point triangle@(virtex1, virtex2, virtex3) = distance point closestPointToTriangleCenteredSorted
   where
     -- Reorder triangles such that we use the corner away from the longest side to address the space in barycentric coordinates.
     closestPointToTriangleCenteredSorted :: V3 ℝ
-    closestPointToTriangleCenteredSorted = closestPointToTriangleCentered adjustedTriangle
+    closestPointToTriangleCenteredSorted = closestPointToTriangle triangle point -- closestPointToTriangleCentered adjustedTriangle
       where
         adjustedTriangle
           | abLength >= bcLength && abLength >= caLength = (virtex3, virtex1, virtex2)
@@ -291,53 +268,56 @@ distancePointToTriangle point (virtex1, virtex2, virtex3) = distance point close
         originDistance = 1/3 *^ (vir1 + vir2 + vir3)
         adjustedTriangle = (vir1 - originDistance, vir2 - originDistance, vir3 - originDistance)
         adjustedPoint = point - originDistance
-    closestPointToTriangle :: (V3 ℝ,V3 ℝ,V3 ℝ) -> V3 ℝ -> V3 ℝ
-    closestPointToTriangle (v1, v2, v3) p
-      -- Closest to the virtices
-      | d1 <=  eps && d2 <= eps = v1
-      | d3  > -eps && d4 <=  d3 = v2
-      | d6  > -eps && d5 <=  d6 = v3
-      -- Closest to the edges
-      | vc <= eps && d1 > -eps && d3 <=  eps = v1 + (d1 / (d1 - d3)) *^ vec12
-      | vb <= eps && d2 > -eps && d6 <=  eps = v1 + (d2 / (d2 - d6)) *^ vec13
-      | va <= eps && dx > -eps && dy  > -eps = v2 + (dx / (dx + dy)) *^ vec23
-      -- On the triangle's surface
-      | otherwise = v1 + v *^ vec12 + w *^ vec13
-      where
-        -- fudge factor. chosen by tuning with our unit pyramid.
-        -- 3,11 = 52 , 12 == 50, 13,14,15,18 == 45
-        eps :: ℝ
-        eps = 1e-13
-        -- The distance along edge12 and edge23, for segment V1 -> P when translated onto the triangle's plane.
-        -- (P when translaned? Read: a line is drawn down to the plane the triangle is on, from p, to a point that is at a right angle with said line.)
-        d1 = vec12 `dot` vec1p
-        d2 = vec13 `dot` vec1p
-        -- Our edge vectors. We have picked v1 to address the space by, for convenience.
-        vec12 = v2 - v1
-        vec13 = v3 - v1
-        -- A segment between our point, and chosen virtex.
-        vec1p =  p - v1
-        -- Distance along edge12 and edge23, for segment V2 -> P when translated onto the triangle's plane.
-        d3 = vec12 `dot` vec2p
-        d4 = vec13 `dot` vec2p
-        -- A segment between our point, and the second virtex.
-        vec2p =  p - v2
-        -- Distance along edge12 and edge23, for segment V3 -> P when translated onto the triangle's plane.
-        d5 = vec12 `dot` vec3p
-        d6 = vec13 `dot` vec3p
-        -- A segment between our point, and the third virtex.
-        vec3p =  p - v3
-        -- An edge vector, along edge23.
-        vec23 = v3 - v2
-        -- The fractional denenominators.
-        va = d1 * d4 - d3 * d2
-        vb = d2 * d5 - d6 * d1
-        vc = d3 * d6 - d5 * d4
-        -- Two convienience values, to make the spacing on the formulas above work.
-        dx = d4 - d3
-        dy = d5 - d6
-        -- The denominator.
-        denom = va + vb + vc
-        -- barycentric results, where we actually intersect.
-        v = vb / denom
-        w = vc / denom
+
+closestPointToTriangle :: (V3 ℝ,V3 ℝ,V3 ℝ) -> V3 ℝ -> V3 ℝ
+closestPointToTriangle (v1, v2, v3) p
+  -- Closest to the virtices
+  | d1 <= 0 && d2 <=  0 = v1
+  | d3 >= 0 && d4 <= d3 = v2
+  | d6 >= 0 && d5 <= d6 = v3
+  -- Nearest to the edges
+  | va <= 0 && d1 > 0 && d3 <= 0 = v1 + (d1 / (d1 - d3)) *^ vec12
+  | vb <= 0 && d2 > 0 && d6 <= 0 = v1 + (d2 / (d2 - d6)) *^ vec13
+  | vc <= 0 && dx > 0 && dy  > 0 = v2 + (dx / (dx + dy)) *^ vec23
+  -- Exactly on an edge, don't bother dividing by zero, please.
+  | denom == 0 = p
+  -- On the triangle's surface
+  | otherwise = v1 + v *^ vec12 + w *^ vec13
+  where
+    -- fudge factor. chosen by tuning with our unit pyramid.
+    -- 3,11 = 52 , 12 == 50, 13,14,15,18 == 45
+    eps :: ℝ
+    eps = 1e-13
+    -- The distance along edge12 and edge23, for segment V1 -> P when translated onto the triangle's plane.
+    -- (P when translaned? Read: a line is drawn down to the plane the triangle is on, from p, to a point that is at a right angle with said line.)
+    d1 = vec12 `dot` vec1p
+    d2 = vec13 `dot` vec1p
+    -- Our edge vectors. We have picked v1 to address the space by, for convenience.
+    vec12 = v2 - v1
+    vec13 = v3 - v1
+    -- A segment between our point, and chosen virtex.
+    vec1p =  p - v1
+    -- Distance along edge12 and edge23, for segment V2 -> P when translated onto the triangle's plane.
+    d3 = vec12 `dot` vec2p
+    d4 = vec13 `dot` vec2p
+    -- A segment between our point, and the second virtex.
+    vec2p =  p - v2
+    -- Distance along edge12 and edge23, for segment V3 -> P when translated onto the triangle's plane.
+    d5 = vec12 `dot` vec3p
+    d6 = vec13 `dot` vec3p
+    -- A segment between our point, and the third virtex.
+    vec3p =  p - v3
+    -- An edge vector, along edge23.
+    vec23 = v3 - v2
+    -- The fractional denenominators.
+    va = d1 * d4 - d3 * d2
+    vb = d2 * d5 - d6 * d1
+    vc = d3 * d6 - d5 * d4
+    -- Two convienience values, to make the spacing on the formulas above work.
+    dx = d4 - d3
+    dy = d5 - d6
+    -- The denominator.
+    denom = va + vb + vc
+    -- barycentric results, where we actually intersect.
+    v = vb / denom
+    w = va / denom

From 60889245d464a0a6f3220384c60725e89cee3b8a Mon Sep 17 00:00:00 2001
From: Julia Longtin <julia.longtin@gmail.com>
Date: Wed, 11 Mar 2026 16:56:09 +0000
Subject: [PATCH 13/28] convex polyhedrons now work.

---
 Graphics/Implicit/ObjectUtil/GetImplicit3.hs | 167 +++++++++++++------
 1 file changed, 113 insertions(+), 54 deletions(-)

diff --git a/Graphics/Implicit/ObjectUtil/GetImplicit3.hs b/Graphics/Implicit/ObjectUtil/GetImplicit3.hs
index 4f2b4954..6d5f8398 100644
--- a/Graphics/Implicit/ObjectUtil/GetImplicit3.hs
+++ b/Graphics/Implicit/ObjectUtil/GetImplicit3.hs
@@ -6,16 +6,24 @@
 module Graphics.Implicit.ObjectUtil.GetImplicit3 (getImplicit3, distancePointToTriangle, closestPointToTriangle) where
 
 -- Import only what we need from the prelude.
-import Prelude (abs, atan2, cos, ceiling, either, error, floor, fromInteger, id, max, min, minimum, negate, otherwise, pi, pure, round, show, sin, sqrt, (||), (/=), Either(Left, Right), (<), (<=), (<>), (>), (>=), (&&), (-), (/), (*), (+), ($), (.), Bool(True, False), (==), (**), Num, Applicative, (<$>))
+import Prelude (abs, atan2, acos, cos, ceiling, either, error, floor, fromInteger, fromIntegral, id, max, min, minimum, negate, otherwise, pi, pure, round, show, snd, sin, sqrt, sum, (||), (/=), Either(Left, Right), (<), (<=), (<>), (>), (>=), (&&), (-), (/), (*), (+), (++), ($), (.), Bool(True, False), (==), (**), Num, Applicative, Int, (<$>), Eq)
 
 import Control.Lens ((^.))
 
 import qualified Data.Either as Either (either)
 
-import Data.List (genericIndex, length, minimumBy)
+import Data.Foldable (toList)
+
+import Data.List (concatMap, genericIndex, length, minimumBy)
+
+import Data.Map (fromListWith, lookup, Map)
+
+import Data.Maybe (fromMaybe)
 
 import Data.Ord (compare)
 
+import Data.Sequence (fromList, mapWithIndex)
+
 import Linear (V2(V2), V3(V3), _xy, _z, distance, dot)
 import qualified Linear (conjugate, inv44, normalizePoint, normalize, point, rotate, Metric)
 
@@ -80,37 +88,58 @@ getImplicit3 _ (Polyhedron [] _) = error "Asked to find distance to an empty pol
 getImplicit3 _ (Polyhedron _ []) = error "Asked to find distance to an empty polygon. No tris."
 getImplicit3 _ (Polyhedron points tris) = \(point) ->
   let
-    (res, closestTri) = unsignedDistanceAndTriangleClosestTo point
+    ((feature,res), closestTri) = unsignedDistanceAndTriangleClosestTo point
   in
-    if pointOnOutside point closestTri
+    if pointOnOutside point (findTriangle points closestTri) closestTri feature
     then          res
     else negate $ res
   where
-    unsignedDistanceAndTriangleClosestTo point = minimumBy (\(a,_) (b,_) -> a `compare` b) $ distTris point
-    distTris point = (\a -> (distancePointToTriangle point a, a)) <$> triangles
-    -- NOTE: May get slightly different values, depending on the V selected.
-    normOfTri :: (V3 ℝ,V3 ℝ,V3 ℝ) -> V3 ℝ
-    normOfTri (v1,v2,v3) = Linear.normalize $ (v2-v1) `cross` (v3-v1)
+    unsignedDistanceAndTriangleClosestTo point = minimumBy (\((_,a),_) ((_,b),_) -> a `compare` b) $ featDistTriangles point
+    featDistTriangles point = (\a -> (distancePointToTriangle point (findTriangle points a), a)) <$> tris
     firstPointOfTri (v1,_,_) = v1
-    pointOnOutside point closestTri = (point - firstPointOfTri closestTri) `dot` normOfTri closestTri >= -eps
+    pointOnOutside :: ℝ3 -> (ℝ3,ℝ3,ℝ3) -> (ℕ,ℕ,ℕ) -> ClosestFeature -> Bool
+    pointOnOutside point closestTriangle closestTri feature = (point - firstPointOfTri closestTriangle) `dot` (weighedNormish points closestTri feature) >= -eps
       where
         -- fudge factor.
         eps :: ℝ
         eps = 1e-13
-    -- decompose our tris into triangles.
-    triangles = findTriangle points <$> tris
-    -- FIXME: make these indices correct by construction?
-    findTriangle :: [V3 ℝ] -> (ℕ,ℕ,ℕ) -> (V3 ℝ, V3 ℝ, V3 ℝ)
-    findTriangle virtices (i1,i2,i3)
-      | outOfRange i1 = error $ "bad virtex index(out of range): " <> show i1 <> "\n"
-      | outOfRange i2 = error $ "bad virtex index(out of range): " <> show i2 <> "\n"
-      | outOfRange i3 = error $ "bad virtex index(out of range): " <> show i3 <> "\n"
-      -- FIXME: there are many more degenerate forms of polyhedron possible here. move polyhedron to only holding a mesh?
-      | otherwise = (genericIndex virtices i1, genericIndex virtices i2, genericIndex virtices i3)
-        where
-          -- FIXME: >=BASE-4.21: replace this with compareLength once debian stable ships 4.21.
-          outOfRange :: ℕ -> Bool
-          outOfRange v = v < 0 || length virtices <= fromℕ v
+        triSeq = fromList tris
+        -- For each edge, the tri indexes that share that edge: 
+        triByEdge :: Map (ℕ,ℕ) [Int]
+        triByEdge = fromListWith (++) edgeTris
+          where
+            edgeTris = concatMap edgesOfTri $ toList $ mapWithIndex (,) triSeq
+            edgesOfTri :: (Int,(ℕ,ℕ,ℕ)) -> [((ℕ,ℕ),[Int])]
+            edgesOfTri (i,(p1,p2,p3)) = [(sortEdge p1 p2,[i]),(sortEdge p2 p3,[i]),(sortEdge p3 p1,[i])]
+        sortEdge a b = (min a b, max a b)
+        -- For each vertex, the tri indexes that contain that vertex: 
+        triByVertex :: Map ℕ [Int]
+        triByVertex = fromListWith (++) vertexTris
+          where
+            vertexTris = concatMap vertexesOfTri $ toList $ mapWithIndex (,) triSeq
+            vertexesOfTri :: (Int,(ℕ,ℕ,ℕ)) -> [(ℕ,[Int])]
+            vertexesOfTri (i,(p1,p2,p3)) = [(p1,[i]),(p2,[i]),(p3,[i])]  
+        -- Get the normalized average of a set of triangles, referred to by index.
+        averageNorm triangles triIndexes = Linear.normalize $ sum $ normOfTriangle . genericIndex triangles <$> triIndexes
+        weighedNormish :: [ℝ3] -> (ℕ,ℕ,ℕ) -> ClosestFeature -> ℝ3
+        weighedNormish points tri@(p1,p2,p3) closest
+          | closest == FeatFace = normOfTriangle closestTriangle
+          | closest == FeatEdge12 = averageNorm triangles $ fromMaybe [] $ lookup (sortEdge p1 p2) triByEdge 
+          | closest == FeatEdge13 = averageNorm triangles $ fromMaybe [] $ lookup (sortEdge p1 p3) triByEdge 
+          | closest == FeatEdge23 = averageNorm triangles $ fromMaybe [] $ lookup (sortEdge p2 p3) triByEdge
+          | closest == FeatVertex1 = Linear.normalize $ sum $ angleWeighed (genericIndex points p1) <$> fromMaybe [] (lookup p1 triByVertex)
+          | closest == FeatVertex2 = Linear.normalize $ sum $ angleWeighed (genericIndex points p2) <$> fromMaybe [] (lookup p2 triByVertex)
+          | closest == FeatVertex3 = Linear.normalize $ sum $ angleWeighed (genericIndex points p3) <$> fromMaybe [] (lookup p3 triByVertex)
+          | otherwise = normOfTriangle closestTriangle
+          where
+            closestTriangle = findTriangle points tri
+        angleWeighed :: ℝ3 -> Int -> ℝ3
+        angleWeighed vertex triNo = angleAt vertex triangle *^ normOfTriangle triangle
+          where
+            triangle = findTriangle points $ genericIndex tris triNo
+        -- decompose our tris into triangles.
+        triangles = findTriangle points <$> tris
+
 getImplicit3 _ (BoxFrame b e) = \p' ->
     let p@(V3 px py pz) = abs p' - b
         V3 qx qy qz = abs (p + pure e) - pure e
@@ -225,9 +254,9 @@ getImplicit3 ctx (RotateExtrude totalRotation translate rotate symbObj) =
                         [0 .. floor $ (totalRotation - θ) / tau]
             n <- ns
             let
-                θvirt = fromℕtoℝ n * tau + θ
-                (V2 rshift zshift) = translate' θvirt
-                twist = rotate' θvirt
+                θvert = fromℕtoℝ n * tau + θ
+                (V2 rshift zshift) = translate' θvert
+                twist = rotate' θvert
                 rz_pos = if twists
                         then let
                             (c,s) = (cos twist, sin twist)
@@ -238,51 +267,81 @@ getImplicit3 ctx (RotateExtrude totalRotation translate rotate symbObj) =
             pure $
               if capped
               then rmax round'
-                    (abs (θvirt - (totalRotation / 2)) - (totalRotation / 2))
+                    (abs (θvert - (totalRotation / 2)) - (totalRotation / 2))
                     (obj rz_pos)
               else obj rz_pos
 getImplicit3 ctx (Shared3 obj) = getImplicitShared ctx obj
 
+-- The closest part of a triangle to a point.
+data ClosestFeature = FeatVertex1 | FeatVertex2 | FeatVertex3 | FeatEdge12 | FeatEdge13 | FeatEdge23 | FeatFace
+  deriving Eq
+
+-- FIXME: make these indices correct by construction?
+findTriangle :: [V3 ℝ] -> (ℕ,ℕ,ℕ) -> (V3 ℝ, V3 ℝ, V3 ℝ)
+findTriangle vertices (i1,i2,i3)
+  | outOfRange i1 = error $ "bad vertex index(out of range): " <> show i1 <> "\n"
+  | outOfRange i2 = error $ "bad vertex index(out of range): " <> show i2 <> "\n"
+  | outOfRange i3 = error $ "bad vertex index(out of range): " <> show i3 <> "\n"
+  -- FIXME: there are many more degenerate forms of polyhedron possible here. move polyhedron to only holding a mesh?
+  | otherwise = (genericIndex vertices i1, genericIndex vertices i2, genericIndex vertices i3)
+  where
+    -- FIXME: >=BASE-4.21: replace this with compareLength once debian stable ships 4.21.
+    outOfRange :: ℕ -> Bool
+    outOfRange v = v < 0 || length vertices <= fromℕ v
+
+normOfTriangle :: (ℝ3,ℝ3,ℝ3) -> ℝ3
+normOfTriangle (v1,v2,v3) = Linear.normalize $ (v2-v1) `cross` (v3-v1)
+
+angleAt :: ℝ3 -> (ℝ3,ℝ3,ℝ3) -> ℝ
+angleAt vertex (v1,v2,v3)
+  | vertex == v1 = acos $ clamp $ Linear.normalize (v2-v1) `dot` Linear.normalize (v3-v1)
+  | vertex == v2 = acos $ clamp $ Linear.normalize (v1-v2) `dot` Linear.normalize (v3-v2)
+  | vertex == v3 = acos $ clamp $ Linear.normalize (v1-v3) `dot` Linear.normalize (v2-v3)
+  where
+    clamp = max (-1) . min 1
+
 -- | You see, what I thought I'd do is put a raytracer inside of a raytracer... what could go wrong...
 -- With inspiration from: https://github.com/RenderKit/embree/blob/master/tutorials/common/math/closest_point.h
-distancePointToTriangle :: V3 ℝ -> (V3 ℝ,V3 ℝ,V3 ℝ) -> ℝ
-distancePointToTriangle point triangle@(virtex1, virtex2, virtex3) = distance point closestPointToTriangleCenteredSorted
+distancePointToTriangle :: V3 ℝ -> (V3 ℝ,V3 ℝ,V3 ℝ) -> (ClosestFeature, ℝ)
+distancePointToTriangle point triangle@(vertex1, vertex2, vertex3) = (resFeature, distance point res)
   where
+    (resFeature, res) = closestPointToTriangleCenteredSorted
     -- Reorder triangles such that we use the corner away from the longest side to address the space in barycentric coordinates.
-    closestPointToTriangleCenteredSorted :: V3 ℝ
+    closestPointToTriangleCenteredSorted :: (ClosestFeature, ℝ3)
     closestPointToTriangleCenteredSorted = closestPointToTriangle triangle point -- closestPointToTriangleCentered adjustedTriangle
       where
         adjustedTriangle
-          | abLength >= bcLength && abLength >= caLength = (virtex3, virtex1, virtex2)
-          | abLength >= caLength                         = (virtex1, virtex2, virtex3)
-          | otherwise                                    = (virtex2, virtex3, virtex1)
+          | abLength >= bcLength && abLength >= caLength = (vertex3, vertex1, vertex2)
+          | abLength >= caLength                         = (vertex1, vertex2, vertex3)
+          | otherwise                                    = (vertex2, vertex3, vertex1)
           where
             -- Really, using length-squared. don't have to abs it, don't have to sqrt it.
-            abLength = (virtex2-virtex1) `dot` (virtex2-virtex1)
-            bcLength = (virtex3-virtex2) `dot` (virtex3-virtex2)
-            caLength = (virtex1-virtex3) `dot` (virtex1-virtex3)
+            abLength = (vertex2-vertex1) `dot` (vertex2-vertex1)
+            bcLength = (vertex3-vertex2) `dot` (vertex3-vertex2)
+            caLength = (vertex1-vertex3) `dot` (vertex1-vertex3)
     -- Force closestPointToTriangle to work at the origin
-    closestPointToTriangleCentered :: (V3 ℝ,V3 ℝ,V3 ℝ) -> V3 ℝ
-    closestPointToTriangleCentered (vir1, vir2, vir3) = originDistance + closestPointToTriangle adjustedTriangle adjustedPoint
+    closestPointToTriangleCentered :: (V3 ℝ,V3 ℝ,V3 ℝ) -> (ClosestFeature, ℝ3)
+    closestPointToTriangleCentered (ver1, ver2, ver3) = (resFeature, originDistance + res)
       where
-        originDistance = 1/3 *^ (vir1 + vir2 + vir3)
-        adjustedTriangle = (vir1 - originDistance, vir2 - originDistance, vir3 - originDistance)
+        (resFeature, res) = closestPointToTriangle adjustedTriangle adjustedPoint
+        originDistance = 1/3 *^ (ver1 + ver2 + ver3)
+        adjustedTriangle = (ver1 - originDistance, ver2 - originDistance, ver3 - originDistance)
         adjustedPoint = point - originDistance
 
-closestPointToTriangle :: (V3 ℝ,V3 ℝ,V3 ℝ) -> V3 ℝ -> V3 ℝ
+closestPointToTriangle :: (V3 ℝ,V3 ℝ,V3 ℝ) -> V3 ℝ -> (ClosestFeature, V3 ℝ)
 closestPointToTriangle (v1, v2, v3) p
-  -- Closest to the virtices
-  | d1 <= 0 && d2 <=  0 = v1
-  | d3 >= 0 && d4 <= d3 = v2
-  | d6 >= 0 && d5 <= d6 = v3
+  -- Closest to the vertices
+  | d1 <= 0 && d2 <=  0 = (FeatVertex1, v1)
+  | d3 >= 0 && d4 <= d3 = (FeatVertex2, v2)
+  | d6 >= 0 && d5 <= d6 = (FeatVertex3, v3)
   -- Nearest to the edges
-  | va <= 0 && d1 > 0 && d3 <= 0 = v1 + (d1 / (d1 - d3)) *^ vec12
-  | vb <= 0 && d2 > 0 && d6 <= 0 = v1 + (d2 / (d2 - d6)) *^ vec13
-  | vc <= 0 && dx > 0 && dy  > 0 = v2 + (dx / (dx + dy)) *^ vec23
+  | va <= 0 && d1 > 0 && d3 <= 0 = (FeatEdge12, v1 + (d1 / (d1 - d3)) *^ vec12)
+  | vb <= 0 && d2 > 0 && d6 <= 0 = (FeatEdge13, v1 + (d2 / (d2 - d6)) *^ vec13)
+  | vc <= 0 && dx > 0 && dy  > 0 = (FeatEdge23, v2 + (dx / (dx + dy)) *^ vec23)
   -- Exactly on an edge, don't bother dividing by zero, please.
-  | denom == 0 = p
+  | denom == 0 = (FeatFace, p)
   -- On the triangle's surface
-  | otherwise = v1 + v *^ vec12 + w *^ vec13
+  | otherwise = (FeatFace, v1 + v *^ vec12 + w *^ vec13)
   where
     -- fudge factor. chosen by tuning with our unit pyramid.
     -- 3,11 = 52 , 12 == 50, 13,14,15,18 == 45
@@ -295,17 +354,17 @@ closestPointToTriangle (v1, v2, v3) p
     -- Our edge vectors. We have picked v1 to address the space by, for convenience.
     vec12 = v2 - v1
     vec13 = v3 - v1
-    -- A segment between our point, and chosen virtex.
+    -- A segment between our point, and chosen vertex.
     vec1p =  p - v1
     -- Distance along edge12 and edge23, for segment V2 -> P when translated onto the triangle's plane.
     d3 = vec12 `dot` vec2p
     d4 = vec13 `dot` vec2p
-    -- A segment between our point, and the second virtex.
+    -- A segment between our point, and the second vertex.
     vec2p =  p - v2
     -- Distance along edge12 and edge23, for segment V3 -> P when translated onto the triangle's plane.
     d5 = vec12 `dot` vec3p
     d6 = vec13 `dot` vec3p
-    -- A segment between our point, and the third virtex.
+    -- A segment between our point, and the third vertex.
     vec3p =  p - v3
     -- An edge vector, along edge23.
     vec23 = v3 - v2

From b3c0c7754afd8c7548ecf2360a524c214987e64c Mon Sep 17 00:00:00 2001
From: Julia Longtin <julia.longtin@gmail.com>
Date: Wed, 11 Mar 2026 22:35:42 +0000
Subject: [PATCH 14/28] add utility file for containing triangle logic.

---
 Graphics/Implicit/TriUtil.hs | 147 +++++++++++++++++++++++++++++++++++
 1 file changed, 147 insertions(+)
 create mode 100644 Graphics/Implicit/TriUtil.hs

diff --git a/Graphics/Implicit/TriUtil.hs b/Graphics/Implicit/TriUtil.hs
new file mode 100644
index 00000000..8cd2d7ab
--- /dev/null
+++ b/Graphics/Implicit/TriUtil.hs
@@ -0,0 +1,147 @@
+-- Implicit CAD. Copyright (C) 2011, Christopher Olah (chris@colah.ca)
+-- Copyright 2015 2016, Mike MacHenry (mike.machenry@gmail.com)
+-- Copyright 2014-2026, Julia Longtin (julia.longtin@gmail.com)
+-- Released under the GNU AGPLV3+, see LICENSE
+
+module Graphics.Implicit.TriUtil (
+  angleAt,
+  closestPointToTriangle,
+  distancePointToTriangle,
+  findTriangle,
+  normOfTriangle,
+  ClosestFeature(FeatFace,
+                 FeatVertex1, FeatVertex2, FeatVertex3,
+                 FeatEdge12,  FeatEdge13,  FeatEdge23),
+  Tri
+  ) where
+
+import Prelude (acos, error, length, max, min, otherwise, show, (>), (<), (&&), (<=), (>=),($), (.), (/), (+), (-), (*), (==), (||), (<>), Bool, Eq)
+
+import Graphics.Implicit.Definitions (
+  fromℕ,
+  ℝ,
+  ℝ3,
+  ℕ)
+
+import Data.List (genericIndex)
+
+import Linear (dot, distance)
+import qualified Linear (normalize)
+
+-- The cross product.
+import Linear.V3 (cross)
+
+-- Matrix times scalar.
+import Linear.Vector ((*^))
+
+type Tri = (ℕ,ℕ,ℕ)
+
+type Triangle = (ℝ3, ℝ3, ℝ3)
+
+-- The closest part of a triangle to a point.
+data ClosestFeature = FeatVertex1 | FeatVertex2 | FeatVertex3 | FeatEdge12 | FeatEdge13 | FeatEdge23 | FeatFace
+  deriving Eq
+
+-- FIXME: make these indices correct by construction?
+findTriangle :: [ℝ3] -> Tri -> Triangle
+findTriangle vertices (i1,i2,i3)
+  | outOfRange i1 = error $ "bad vertex index(out of range): " <> show i1 <> "\n"
+  | outOfRange i2 = error $ "bad vertex index(out of range): " <> show i2 <> "\n"
+  | outOfRange i3 = error $ "bad vertex index(out of range): " <> show i3 <> "\n"
+  -- FIXME: there are many more degenerate forms of polyhedron possible here. move polyhedron to only holding a mesh?
+  | otherwise = (genericIndex vertices i1, genericIndex vertices i2, genericIndex vertices i3)
+  where
+    -- FIXME: >=BASE-4.21: replace this with compareLength once debian stable ships 4.21.
+    outOfRange :: ℕ -> Bool
+    outOfRange v = v < 0 || length vertices <= fromℕ v
+
+normOfTriangle :: Triangle -> ℝ3
+normOfTriangle (v1,v2,v3) = Linear.normalize $ (v2-v1) `cross` (v3-v1)
+
+angleAt :: ℝ3 -> Triangle -> ℝ
+angleAt vertex (v1,v2,v3)
+  | vertex == v1 = acos $ clamp $ Linear.normalize (v2-v1) `dot` Linear.normalize (v3-v1)
+  | vertex == v2 = acos $ clamp $ Linear.normalize (v1-v2) `dot` Linear.normalize (v3-v2)
+  | vertex == v3 = acos $ clamp $ Linear.normalize (v1-v3) `dot` Linear.normalize (v2-v3)
+  | otherwise = error $ "tried to get angleAt with a point not one of the vertexes: " <> show vertex <> "\n"
+  where
+    clamp = max (-1) . min 1
+
+-- | You see, what I thought I'd do is put a raytracer inside of a raytracer... what could go wrong...
+-- With inspiration from: https://github.com/RenderKit/embree/blob/master/tutorials/common/math/closest_point.h
+distancePointToTriangle :: ℝ3 -> Triangle -> (ClosestFeature, ℝ)
+distancePointToTriangle point triangle@(vertex1, vertex2, vertex3) = (resFeature, distance point res)
+  where
+    (resFeature, res) = closestPointToTriangleCenteredSorted
+    -- Reorder triangles such that we use the corner away from the longest side to address the space in barycentric coordinates.
+    closestPointToTriangleCenteredSorted :: (ClosestFeature, ℝ3)
+-- FIXME: test the following, rather than what's here.
+--    closestPointToTriangleCenteredSorted = closestPointToTriangleCentered adjustedTriangle
+    closestPointToTriangleCenteredSorted = closestPointToTriangle triangle point
+      where
+        adjustedTriangle
+          | abLength >= bcLength && abLength >= caLength = (vertex3, vertex1, vertex2)
+          | abLength >= caLength                         = (vertex1, vertex2, vertex3)
+          | otherwise                                    = (vertex2, vertex3, vertex1)
+          where
+            -- Really, using length-squared. don't have to abs it, don't have to sqrt it.
+            abLength = (vertex2-vertex1) `dot` (vertex2-vertex1)
+            bcLength = (vertex3-vertex2) `dot` (vertex3-vertex2)
+            caLength = (vertex1-vertex3) `dot` (vertex1-vertex3)
+    -- Force closestPointToTriangle to work at the origin
+    closestPointToTriangleCentered :: Triangle -> (ClosestFeature, ℝ3)
+    closestPointToTriangleCentered (ver1, ver2, ver3) = (resFeature, originDistance + res)
+      where
+        (resFeature, res) = closestPointToTriangle adjustedTriangle adjustedPoint
+        originDistance = 1/3 *^ (ver1 + ver2 + ver3)
+        adjustedTriangle = (ver1 - originDistance, ver2 - originDistance, ver3 - originDistance)
+        adjustedPoint = point - originDistance
+
+closestPointToTriangle :: Triangle -> ℝ3 -> (ClosestFeature, ℝ3)
+closestPointToTriangle (v1, v2, v3) p
+  -- Closest to the vertices
+  | d1 <= 0 && d2 <=  0 = (FeatVertex1, v1)
+  | d3 >= 0 && d4 <= d3 = (FeatVertex2, v2)
+  | d6 >= 0 && d5 <= d6 = (FeatVertex3, v3)
+  -- Nearest to the edges
+  | va <= 0 && d1 > 0 && d3 <= 0 = (FeatEdge12, v1 + (d1 / (d1 - d3)) *^ vec12)
+  | vb <= 0 && d2 > 0 && d6 <= 0 = (FeatEdge13, v1 + (d2 / (d2 - d6)) *^ vec13)
+  | vc <= 0 && dx > 0 && dy  > 0 = (FeatEdge23, v2 + (dx / (dx + dy)) *^ vec23)
+  -- Exactly on an edge, don't bother dividing by zero, please.
+  | denom == 0 = (FeatFace, p)
+  -- On the triangle's surface
+  | otherwise = (FeatFace, v1 + v *^ vec12 + w *^ vec13)
+  where
+    -- The distance along edge12 and edge23, for segment V1 -> P when translated onto the triangle's plane.
+    -- (P when translaned? Read: a line is drawn down to the plane the triangle is on, from p, to a point that is at a right angle with said line.)
+    d1 = vec12 `dot` vec1p
+    d2 = vec13 `dot` vec1p
+    -- Our edge vectors. We have picked v1 to address the space by, for convenience.
+    vec12 = v2 - v1
+    vec13 = v3 - v1
+    -- A segment between our point, and chosen vertex.
+    vec1p =  p - v1
+    -- Distance along edge12 and edge23, for segment V2 -> P when translated onto the triangle's plane.
+    d3 = vec12 `dot` vec2p
+    d4 = vec13 `dot` vec2p
+    -- A segment between our point, and the second vertex.
+    vec2p =  p - v2
+    -- Distance along edge12 and edge23, for segment V3 -> P when translated onto the triangle's plane.
+    d5 = vec12 `dot` vec3p
+    d6 = vec13 `dot` vec3p
+    -- A segment between our point, and the third vertex.
+    vec3p =  p - v3
+    -- An edge vector, along edge23.
+    vec23 = v3 - v2
+    -- The fractional denenominators.
+    va = d1 * d4 - d3 * d2
+    vb = d2 * d5 - d6 * d1
+    vc = d3 * d6 - d5 * d4
+    -- Two convienience values, to make the spacing on the formulas above work.
+    dx = d4 - d3
+    dy = d5 - d6
+    -- The denominator.
+    denom = va + vb + vc
+    -- barycentric results, where we actually intersect.
+    v = vb / denom
+    w = va / denom

From 883092c4e5de7842fd133a4ad59a2d78dbf43f43 Mon Sep 17 00:00:00 2001
From: Julia Longtin <julia.longtin@gmail.com>
Date: Wed, 11 Mar 2026 22:36:05 +0000
Subject: [PATCH 15/28] use triangle utility file.

---
 Graphics/Implicit/ExtOpenScad/Primitives.hs  |   5 +-
 Graphics/Implicit/ObjectUtil/GetImplicit3.hs | 117 +------------------
 2 files changed, 8 insertions(+), 114 deletions(-)

diff --git a/Graphics/Implicit/ExtOpenScad/Primitives.hs b/Graphics/Implicit/ExtOpenScad/Primitives.hs
index e153a6a0..247740cf 100644
--- a/Graphics/Implicit/ExtOpenScad/Primitives.hs
+++ b/Graphics/Implicit/ExtOpenScad/Primitives.hs
@@ -31,6 +31,8 @@ import Graphics.Implicit.ExtOpenScad.Util.OVal (OTypeMirror, caseOType, divideOb
 
 import Graphics.Implicit.ExtOpenScad.Util.StateC (errorC, warnC)
 
+import Graphics.Implicit.TriUtil (Tri)
+
 -- Note the use of a qualified import, so we don't have the functions in this file conflict with what we're importing.
 import qualified Graphics.Implicit.Primitives as Prim (withRounding, sphere, rect3, rect, translate, circle, polygon, polyhedron, extrude, cylinder2, union, unionR, intersect, intersectR, difference, differenceR, rotate, slice, transform, rotate3V, rotate3, transform3, scale, extrudeM, rotateExtrude, shell, mirror, pack3, pack2, torus, ellipsoid, cone)
 
@@ -288,9 +290,6 @@ cylinder = moduleWithoutSuite "cylinder" $ \_ -> do
         in shift obj3
         else shift $ Prim.cylinder2 r1 r2 dh
 
--- A triangle, with the points in an external index.
-type Tri = (ℕ,ℕ,ℕ)
-
 polyhedron :: (Symbol, SourcePosition -> ArgParser (StateC [OVal]))
 polyhedron = moduleWithoutSuite "polyhedron" $ \sourcePos -> do
     example "polyhedron(points=[[0,0,0], [2,0,0], [2,2,0], [0,2,0], [1, 1, 2]], faces=[[0,1,2,3], [0,5,1], [1,5,2], [2,5,3], [3,5,4], [4,5,0]]);"
diff --git a/Graphics/Implicit/ObjectUtil/GetImplicit3.hs b/Graphics/Implicit/ObjectUtil/GetImplicit3.hs
index 6d5f8398..75159af5 100644
--- a/Graphics/Implicit/ObjectUtil/GetImplicit3.hs
+++ b/Graphics/Implicit/ObjectUtil/GetImplicit3.hs
@@ -1,9 +1,9 @@
--- Copyright 2014 2015 2016, Julia Longtin (julial@turinglace.com)
--- Copyright 2015 2016, Mike MacHenry (mike.machenry@gmail.com)
 -- Implicit CAD. Copyright (C) 2011, Christopher Olah (chris@colah.ca)
+-- Copyright 2015 2016, Mike MacHenry (mike.machenry@gmail.com)
+-- Copyright 2014 2015 2016, Julia Longtin (julia.longtin@gmail.com)
 -- Released under the GNU AGPLV3+, see LICENSE
 
-module Graphics.Implicit.ObjectUtil.GetImplicit3 (getImplicit3, distancePointToTriangle, closestPointToTriangle) where
+module Graphics.Implicit.ObjectUtil.GetImplicit3 (getImplicit3) where
 
 -- Import only what we need from the prelude.
 import Prelude (abs, atan2, acos, cos, ceiling, either, error, floor, fromInteger, fromIntegral, id, max, min, minimum, negate, otherwise, pi, pure, round, show, snd, sin, sqrt, sum, (||), (/=), Either(Left, Right), (<), (<=), (<>), (>), (>=), (&&), (-), (/), (*), (+), (++), ($), (.), Bool(True, False), (==), (**), Num, Applicative, Int, (<$>), Eq)
@@ -27,6 +27,7 @@ import Data.Sequence (fromList, mapWithIndex)
 import Linear (V2(V2), V3(V3), _xy, _z, distance, dot)
 import qualified Linear (conjugate, inv44, normalizePoint, normalize, point, rotate, Metric)
 
+-- The cross product.
 import Linear.V3 (cross)
 
 -- Matrix times column vector.
@@ -52,6 +53,8 @@ import Graphics.Implicit.Definitions
 -- For handling extrusion of 2D shapes to 3D.
 import {-# SOURCE #-} Graphics.Implicit.Primitives (getImplicit)
 
+import Graphics.Implicit.TriUtil (angleAt, distancePointToTriangle, findTriangle, normOfTriangle, ClosestFeature(..), Tri)
+
 import Graphics.Implicit.MathUtil (rmax, rmaximum)
 
 import Graphics.Implicit.ObjectUtil.GetImplicitShared (getImplicitShared)
@@ -272,111 +275,3 @@ getImplicit3 ctx (RotateExtrude totalRotation translate rotate symbObj) =
               else obj rz_pos
 getImplicit3 ctx (Shared3 obj) = getImplicitShared ctx obj
 
--- The closest part of a triangle to a point.
-data ClosestFeature = FeatVertex1 | FeatVertex2 | FeatVertex3 | FeatEdge12 | FeatEdge13 | FeatEdge23 | FeatFace
-  deriving Eq
-
--- FIXME: make these indices correct by construction?
-findTriangle :: [V3 ℝ] -> (ℕ,ℕ,ℕ) -> (V3 ℝ, V3 ℝ, V3 ℝ)
-findTriangle vertices (i1,i2,i3)
-  | outOfRange i1 = error $ "bad vertex index(out of range): " <> show i1 <> "\n"
-  | outOfRange i2 = error $ "bad vertex index(out of range): " <> show i2 <> "\n"
-  | outOfRange i3 = error $ "bad vertex index(out of range): " <> show i3 <> "\n"
-  -- FIXME: there are many more degenerate forms of polyhedron possible here. move polyhedron to only holding a mesh?
-  | otherwise = (genericIndex vertices i1, genericIndex vertices i2, genericIndex vertices i3)
-  where
-    -- FIXME: >=BASE-4.21: replace this with compareLength once debian stable ships 4.21.
-    outOfRange :: ℕ -> Bool
-    outOfRange v = v < 0 || length vertices <= fromℕ v
-
-normOfTriangle :: (ℝ3,ℝ3,ℝ3) -> ℝ3
-normOfTriangle (v1,v2,v3) = Linear.normalize $ (v2-v1) `cross` (v3-v1)
-
-angleAt :: ℝ3 -> (ℝ3,ℝ3,ℝ3) -> ℝ
-angleAt vertex (v1,v2,v3)
-  | vertex == v1 = acos $ clamp $ Linear.normalize (v2-v1) `dot` Linear.normalize (v3-v1)
-  | vertex == v2 = acos $ clamp $ Linear.normalize (v1-v2) `dot` Linear.normalize (v3-v2)
-  | vertex == v3 = acos $ clamp $ Linear.normalize (v1-v3) `dot` Linear.normalize (v2-v3)
-  where
-    clamp = max (-1) . min 1
-
--- | You see, what I thought I'd do is put a raytracer inside of a raytracer... what could go wrong...
--- With inspiration from: https://github.com/RenderKit/embree/blob/master/tutorials/common/math/closest_point.h
-distancePointToTriangle :: V3 ℝ -> (V3 ℝ,V3 ℝ,V3 ℝ) -> (ClosestFeature, ℝ)
-distancePointToTriangle point triangle@(vertex1, vertex2, vertex3) = (resFeature, distance point res)
-  where
-    (resFeature, res) = closestPointToTriangleCenteredSorted
-    -- Reorder triangles such that we use the corner away from the longest side to address the space in barycentric coordinates.
-    closestPointToTriangleCenteredSorted :: (ClosestFeature, ℝ3)
-    closestPointToTriangleCenteredSorted = closestPointToTriangle triangle point -- closestPointToTriangleCentered adjustedTriangle
-      where
-        adjustedTriangle
-          | abLength >= bcLength && abLength >= caLength = (vertex3, vertex1, vertex2)
-          | abLength >= caLength                         = (vertex1, vertex2, vertex3)
-          | otherwise                                    = (vertex2, vertex3, vertex1)
-          where
-            -- Really, using length-squared. don't have to abs it, don't have to sqrt it.
-            abLength = (vertex2-vertex1) `dot` (vertex2-vertex1)
-            bcLength = (vertex3-vertex2) `dot` (vertex3-vertex2)
-            caLength = (vertex1-vertex3) `dot` (vertex1-vertex3)
-    -- Force closestPointToTriangle to work at the origin
-    closestPointToTriangleCentered :: (V3 ℝ,V3 ℝ,V3 ℝ) -> (ClosestFeature, ℝ3)
-    closestPointToTriangleCentered (ver1, ver2, ver3) = (resFeature, originDistance + res)
-      where
-        (resFeature, res) = closestPointToTriangle adjustedTriangle adjustedPoint
-        originDistance = 1/3 *^ (ver1 + ver2 + ver3)
-        adjustedTriangle = (ver1 - originDistance, ver2 - originDistance, ver3 - originDistance)
-        adjustedPoint = point - originDistance
-
-closestPointToTriangle :: (V3 ℝ,V3 ℝ,V3 ℝ) -> V3 ℝ -> (ClosestFeature, V3 ℝ)
-closestPointToTriangle (v1, v2, v3) p
-  -- Closest to the vertices
-  | d1 <= 0 && d2 <=  0 = (FeatVertex1, v1)
-  | d3 >= 0 && d4 <= d3 = (FeatVertex2, v2)
-  | d6 >= 0 && d5 <= d6 = (FeatVertex3, v3)
-  -- Nearest to the edges
-  | va <= 0 && d1 > 0 && d3 <= 0 = (FeatEdge12, v1 + (d1 / (d1 - d3)) *^ vec12)
-  | vb <= 0 && d2 > 0 && d6 <= 0 = (FeatEdge13, v1 + (d2 / (d2 - d6)) *^ vec13)
-  | vc <= 0 && dx > 0 && dy  > 0 = (FeatEdge23, v2 + (dx / (dx + dy)) *^ vec23)
-  -- Exactly on an edge, don't bother dividing by zero, please.
-  | denom == 0 = (FeatFace, p)
-  -- On the triangle's surface
-  | otherwise = (FeatFace, v1 + v *^ vec12 + w *^ vec13)
-  where
-    -- fudge factor. chosen by tuning with our unit pyramid.
-    -- 3,11 = 52 , 12 == 50, 13,14,15,18 == 45
-    eps :: ℝ
-    eps = 1e-13
-    -- The distance along edge12 and edge23, for segment V1 -> P when translated onto the triangle's plane.
-    -- (P when translaned? Read: a line is drawn down to the plane the triangle is on, from p, to a point that is at a right angle with said line.)
-    d1 = vec12 `dot` vec1p
-    d2 = vec13 `dot` vec1p
-    -- Our edge vectors. We have picked v1 to address the space by, for convenience.
-    vec12 = v2 - v1
-    vec13 = v3 - v1
-    -- A segment between our point, and chosen vertex.
-    vec1p =  p - v1
-    -- Distance along edge12 and edge23, for segment V2 -> P when translated onto the triangle's plane.
-    d3 = vec12 `dot` vec2p
-    d4 = vec13 `dot` vec2p
-    -- A segment between our point, and the second vertex.
-    vec2p =  p - v2
-    -- Distance along edge12 and edge23, for segment V3 -> P when translated onto the triangle's plane.
-    d5 = vec12 `dot` vec3p
-    d6 = vec13 `dot` vec3p
-    -- A segment between our point, and the third vertex.
-    vec3p =  p - v3
-    -- An edge vector, along edge23.
-    vec23 = v3 - v2
-    -- The fractional denenominators.
-    va = d1 * d4 - d3 * d2
-    vb = d2 * d5 - d6 * d1
-    vc = d3 * d6 - d5 * d4
-    -- Two convienience values, to make the spacing on the formulas above work.
-    dx = d4 - d3
-    dy = d5 - d6
-    -- The denominator.
-    denom = va + vb + vc
-    -- barycentric results, where we actually intersect.
-    v = vb / denom
-    w = va / denom

From ea99224022e1702e4d1af9493beb493cde91c42e Mon Sep 17 00:00:00 2001
From: Julia Longtin <julia.longtin@gmail.com>
Date: Wed, 11 Mar 2026 23:11:56 +0000
Subject: [PATCH 16/28] cleanup warnings, build problems.

---
 Graphics/Implicit/ExtOpenScad/Primitives.hs | 39 ++++++++++-----------
 implicit.cabal                              |  1 +
 2 files changed, 20 insertions(+), 20 deletions(-)

diff --git a/Graphics/Implicit/ExtOpenScad/Primitives.hs b/Graphics/Implicit/ExtOpenScad/Primitives.hs
index 247740cf..e01fb849 100644
--- a/Graphics/Implicit/ExtOpenScad/Primitives.hs
+++ b/Graphics/Implicit/ExtOpenScad/Primitives.hs
@@ -318,21 +318,20 @@ polyhedron = moduleWithoutSuite "polyhedron" $ \sourcePos -> do
         -- | Ensure our triangles are wound in the same direction
         reWindTriangles :: SourcePosition -> [ℝ3] -> [Tri] -> StateC [Tri]
         reWindTriangles _         _ [] = pure []
-        -- Really, forces them to have the same winding as the first triangle, from us putting [tri] here.
-        reWindTriangles sourcePos points (tri:moreTris) = windTriangles [safeTri] (fromList moreTris)
+        -- Really, forces them to have the same winding as the first triangle, from us putting [safeTri] here.
+        reWindTriangles sourcePos points (firstTri:moreTris) = windTriangles [safeTri] (fromList moreTris)
           where
-            polyCentroid = centroid points
             -- The first triangle, flipped based on comparing two centroids.
             safeTri
-              | (triCentroid - polyCentroid) `dot` triNorm < 0 = flipTri tri
-              | otherwise                                      = tri
+              | (triCentroid - polyCentroid) `dot` triNorm < 0 = flipTri firstTri
+              | otherwise                                      = firstTri
               where
-                (p1,p2,p3) = tri
+                (p1,p2,p3) = firstTri
                 (v1,v2,v3) = (genericIndex points p1,genericIndex points p2,genericIndex points p3)
                 -- The norm of the safe triangle.
                 triNorm    = (v2-v1) `cross` (v3-v1)
                 triCentroid = centroid [v1,v2,v3]
-                -- type conversion.
+            polyCentroid = centroid points
             flipTri (p1,p2,p3) = (p1,p3,p2)
             -- | Wind the triangles.
             windTriangles :: [Tri] -> Seq Tri -> StateC [Tri]
@@ -343,7 +342,7 @@ polyhedron = moduleWithoutSuite "polyhedron" $ \sourcePos -> do
                   (newVisited, newUnvisited) <- foldM (classifyTri visited) ([], unvisited) (toList unvisited)
                   if null newVisited
                     then do
-                         warnC sourcePos $ "Had to pick a new root"
+                         warnC sourcePos $ "Had to pick a new root, incomplete polyhedron?"
                          windTriangles (visited <> [head $ toList newUnvisited]) (deleteAt 0 newUnvisited)
                     else windTriangles (visited <> newVisited)                   newUnvisited
             -- | Compare one unvisited triangle against all visited triangles.
@@ -352,31 +351,31 @@ polyhedron = moduleWithoutSuite "polyhedron" $ \sourcePos -> do
             classifyTri visited (found, remaining) triUnderTest =
               case res of
                 Just triFound -> do -- See if we flipped our tri, and if so, throw a warning.
-                                 if triFound /= triUnderTest
-                                   then warnC sourcePos $ "Flipped face detected with vertices " <> (DTL.pack $ show triUnderTest)
-                                   else pure ()
-                                 pure (found <> [triFound], filter (/= triUnderTest) remaining)
+                  if triFound /= triUnderTest
+                    then warnC sourcePos $ "Flipped face detected with vertices " <> (DTL.pack $ show triUnderTest)
+                    else pure ()
+                  pure (found <> [triFound], filter (/= triUnderTest) remaining)
                 Nothing ->       pure (found, remaining)
               where
                 res = foldr (\tri state -> firstNeighborFilter tri triUnderTest state) Nothing visited 
                 -- | A short-circuiting filter we fold over visited, and grab the first neighboring tri.
                 firstNeighborFilter :: Tri -> Tri -> Maybe Tri -> Maybe Tri
-                firstNeighborFilter src triUnderTest res
-                  | isJust res = res
-                  | otherwise  = maybeWindNeighbor src triUnderTest
+                firstNeighborFilter src testTri maybeRes
+                  | isJust maybeRes = maybeRes
+                  | otherwise  = maybeWindNeighbor src testTri
                 -- | Checks whether a triangle under test is a neighbor of the given triangle, and if it is, returns if after ensuring it is wound in the proper direction.
                 maybeWindNeighbor :: Tri -> Tri -> Maybe Tri
-                maybeWindNeighbor src triUnderTest
+                maybeWindNeighbor src testTri
                 -- A correctly wound neighbor will have the opposite direction, when it is referring to a given edge.
-                  | oppositeNeighbor = Just           triUnderTest
+                  | oppositeNeighbor = Just           testTri
                 -- A neighbor that needs flipped will refer to an edge in the same direction as our given triangle.
-                  | sameNeighbor     = Just $ flipTri triUnderTest
+                  | sameNeighbor     = Just $ flipTri testTri
                   | otherwise = Nothing
                   where
-                    oppositeNeighbor = any (\edge -> flip edge `elem` edgesOfTri src) $ edgesOfTri triUnderTest
+                    oppositeNeighbor = any (\edge -> flip edge `elem` edgesOfTri src) $ edgesOfTri testTri
                       where
                         flip (a,b) = (b,a)
-                    sameNeighbor     = any (\edge ->      edge `elem` edgesOfTri src) $ edgesOfTri triUnderTest
+                    sameNeighbor     = any (\edge ->      edge `elem` edgesOfTri src) $ edgesOfTri testTri
                     edgesOfTri (p1,p2,p3) = [(p1,p2), (p2,p3), (p3,p1)]
 
 cone :: (Symbol, SourcePosition -> ArgParser (StateC [OVal]))
diff --git a/implicit.cabal b/implicit.cabal
index 9671dbd6..069fa03b 100644
--- a/implicit.cabal
+++ b/implicit.cabal
@@ -149,6 +149,7 @@ Library
                     Graphics.Implicit.Export.Resolution
     Other-modules:
                   Graphics.Implicit.FastIntUtil
+                  Graphics.Implicit.TriUtil
                   Graphics.Implicit.IntegralUtil
                   Graphics.Implicit.ObjectUtil.GetBox2
                   Graphics.Implicit.ObjectUtil.GetBox3

From 7878b3068b62f67128f31cd95bbaf4e8130d35c9 Mon Sep 17 00:00:00 2001
From: Julia Longtin <julia.longtin@gmail.com>
Date: Thu, 12 Mar 2026 01:31:41 +0000
Subject: [PATCH 17/28] fix the precision optimizing transforms in
 distancePointToTriangle.

---
 Graphics/Implicit/TriUtil.hs | 70 ++++++++++++++++++++++--------------
 1 file changed, 44 insertions(+), 26 deletions(-)

diff --git a/Graphics/Implicit/TriUtil.hs b/Graphics/Implicit/TriUtil.hs
index 8cd2d7ab..627e28d5 100644
--- a/Graphics/Implicit/TriUtil.hs
+++ b/Graphics/Implicit/TriUtil.hs
@@ -3,9 +3,12 @@
 -- Copyright 2014-2026, Julia Longtin (julia.longtin@gmail.com)
 -- Released under the GNU AGPLV3+, see LICENSE
 
+-- You see, what I thought I'd do is put a raytracer inside of a raytracer... what could go wrong...
+-- With inspiration from: https://github.com/RenderKit/embree/blob/master/tutorials/common/math/closest_point.h
+
 module Graphics.Implicit.TriUtil (
   angleAt,
-  closestPointToTriangle,
+  closestFeatureToTriangle,
   distancePointToTriangle,
   findTriangle,
   normOfTriangle,
@@ -65,40 +68,55 @@ angleAt vertex (v1,v2,v3)
   | vertex == v3 = acos $ clamp $ Linear.normalize (v1-v3) `dot` Linear.normalize (v2-v3)
   | otherwise = error $ "tried to get angleAt with a point not one of the vertexes: " <> show vertex <> "\n"
   where
+    clamp :: ℝ -> ℝ
     clamp = max (-1) . min 1
 
--- | You see, what I thought I'd do is put a raytracer inside of a raytracer... what could go wrong...
--- With inspiration from: https://github.com/RenderKit/embree/blob/master/tutorials/common/math/closest_point.h
 distancePointToTriangle :: ℝ3 -> Triangle -> (ClosestFeature, ℝ)
-distancePointToTriangle point triangle@(vertex1, vertex2, vertex3) = (resFeature, distance point res)
+distancePointToTriangle point (vertex1, vertex2, vertex3) = (adjustedFeature, distance point pointOnFeature)
   where
-    (resFeature, res) = closestPointToTriangleCenteredSorted
-    -- Reorder triangles such that we use the corner away from the longest side to address the space in barycentric coordinates.
-    closestPointToTriangleCenteredSorted :: (ClosestFeature, ℝ3)
--- FIXME: test the following, rather than what's here.
---    closestPointToTriangleCenteredSorted = closestPointToTriangleCentered adjustedTriangle
-    closestPointToTriangleCenteredSorted = closestPointToTriangle triangle point
+    (resFeature, pointOnFeature) = closestFeatureToTriangleCentered adjustedTriangle point
+    -- First math precision transform: change which adressing system we use for the triangle, ensuring the far side is 'away' from the virtex we address from.
+    adjustedFeature
+      | adjustedTriangle == (vertex3, vertex1, vertex2) =
+        case resFeature of
+          FeatVertex1 -> FeatVertex3
+          FeatVertex2 -> FeatVertex1
+          FeatVertex3 -> FeatVertex2
+          FeatEdge12  -> FeatEdge13
+          FeatEdge13  -> FeatEdge23
+          FeatEdge23  -> FeatEdge12
+          FeatFace -> FeatFace
+      | adjustedTriangle == (vertex2, vertex3, vertex1) =
+        case resFeature of
+          FeatVertex1 -> FeatVertex2
+          FeatVertex2 -> FeatVertex3
+          FeatVertex3 -> FeatVertex1
+          FeatEdge12  -> FeatEdge23
+          FeatEdge13  -> FeatEdge12
+          FeatEdge23  -> FeatEdge13
+          FeatFace -> FeatFace
+      | otherwise = resFeature
+    adjustedTriangle
+      | abLength >= bcLength && abLength >= caLength = (vertex3, vertex1, vertex2)
+      | abLength >= caLength                         = (vertex1, vertex2, vertex3)
+      | otherwise                                    = (vertex2, vertex3, vertex1)
       where
-        adjustedTriangle
-          | abLength >= bcLength && abLength >= caLength = (vertex3, vertex1, vertex2)
-          | abLength >= caLength                         = (vertex1, vertex2, vertex3)
-          | otherwise                                    = (vertex2, vertex3, vertex1)
-          where
-            -- Really, using length-squared. don't have to abs it, don't have to sqrt it.
-            abLength = (vertex2-vertex1) `dot` (vertex2-vertex1)
-            bcLength = (vertex3-vertex2) `dot` (vertex3-vertex2)
-            caLength = (vertex1-vertex3) `dot` (vertex1-vertex3)
-    -- Force closestPointToTriangle to work at the origin
-    closestPointToTriangleCentered :: Triangle -> (ClosestFeature, ℝ3)
-    closestPointToTriangleCentered (ver1, ver2, ver3) = (resFeature, originDistance + res)
+        -- Really, using length-squared. don't have to abs it, don't have to sqrt it.
+        abLength = (vertex2-vertex1) `dot` (vertex2-vertex1)
+        bcLength = (vertex3-vertex2) `dot` (vertex3-vertex2)
+        caLength = (vertex1-vertex3) `dot` (vertex1-vertex3)
+    -- Second math precision transform: Force closestFeatureToTriangle to work near the origin by translating our query, and then translating the response.
+    closestFeatureToTriangleCentered :: Triangle -> ℝ3 -> (ClosestFeature, ℝ3)
+    closestFeatureToTriangleCentered (ver1, ver2, ver3) inpoint = (resFeature, originDistance + res)
       where
-        (resFeature, res) = closestPointToTriangle adjustedTriangle adjustedPoint
+        (resFeature, res) = closestFeatureToTriangle adjustedTriangle adjustedPoint
         originDistance = 1/3 *^ (ver1 + ver2 + ver3)
         adjustedTriangle = (ver1 - originDistance, ver2 - originDistance, ver3 - originDistance)
-        adjustedPoint = point - originDistance
+        adjustedPoint = inpoint - originDistance
 
-closestPointToTriangle :: Triangle -> ℝ3 -> (ClosestFeature, ℝ3)
-closestPointToTriangle (v1, v2, v3) p
+-- | Find the closest part of a triangle (edge, center, vertex) to a given point , along with the point on the closest part that is closest to the given point.
+closestFeatureToTriangle :: Triangle -> ℝ3 -> (ClosestFeature, ℝ3)
+closestFeatureToTriangle (v1, v2, v3) p
   -- Closest to the vertices
   | d1 <= 0 && d2 <=  0 = (FeatVertex1, v1)
   | d3 >= 0 && d4 <= d3 = (FeatVertex2, v2)

From b8eedc138a389a2757cef13d4f200529c2903199 Mon Sep 17 00:00:00 2001
From: Julia Longtin <julia.longtin@gmail.com>
Date: Thu, 12 Mar 2026 01:47:50 +0000
Subject: [PATCH 18/28] remove warnings.

---
 Graphics/Implicit/TriUtil.hs | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/Graphics/Implicit/TriUtil.hs b/Graphics/Implicit/TriUtil.hs
index 627e28d5..10c6de5b 100644
--- a/Graphics/Implicit/TriUtil.hs
+++ b/Graphics/Implicit/TriUtil.hs
@@ -45,7 +45,7 @@ type Triangle = (ℝ3, ℝ3, ℝ3)
 data ClosestFeature = FeatVertex1 | FeatVertex2 | FeatVertex3 | FeatEdge12 | FeatEdge13 | FeatEdge23 | FeatFace
   deriving Eq
 
--- FIXME: make these indices correct by construction?
+-- FIXME: Make these indices correct by construction?
 findTriangle :: [ℝ3] -> Tri -> Triangle
 findTriangle vertices (i1,i2,i3)
   | outOfRange i1 = error $ "bad vertex index(out of range): " <> show i1 <> "\n"
@@ -107,12 +107,12 @@ distancePointToTriangle point (vertex1, vertex2, vertex3) = (adjustedFeature, di
         caLength = (vertex1-vertex3) `dot` (vertex1-vertex3)
     -- Second math precision transform: Force closestFeatureToTriangle to work near the origin by translating our query, and then translating the response.
     closestFeatureToTriangleCentered :: Triangle -> ℝ3 -> (ClosestFeature, ℝ3)
-    closestFeatureToTriangleCentered (ver1, ver2, ver3) inpoint = (resFeature, originDistance + res)
+    closestFeatureToTriangleCentered (ver1, ver2, ver3) inpoint = (feature, originDistance + res)
       where
-        (resFeature, res) = closestFeatureToTriangle adjustedTriangle adjustedPoint
+        (feature, res) = closestFeatureToTriangle translatedTriangle translatedPoint
+        translatedTriangle = (ver1 - originDistance, ver2 - originDistance, ver3 - originDistance)
+        translatedPoint = inpoint - originDistance
         originDistance = 1/3 *^ (ver1 + ver2 + ver3)
-        adjustedTriangle = (ver1 - originDistance, ver2 - originDistance, ver3 - originDistance)
-        adjustedPoint = inpoint - originDistance
 
 -- | Find the closest part of a triangle (edge, center, vertex) to a given point , along with the point on the closest part that is closest to the given point.
 closestFeatureToTriangle :: Triangle -> ℝ3 -> (ClosestFeature, ℝ3)

From 54622e353d2d5ec95af5e29e338cf0439c9072b3 Mon Sep 17 00:00:00 2001
From: Julia Longtin <julia.longtin@gmail.com>
Date: Thu, 12 Mar 2026 02:31:57 +0000
Subject: [PATCH 19/28] remove warnings.

---
 Graphics/Implicit/ObjectUtil/GetImplicit3.hs | 33 ++++++++------------
 Graphics/Implicit/TriUtil.hs                 |  6 +++-
 2 files changed, 18 insertions(+), 21 deletions(-)

diff --git a/Graphics/Implicit/ObjectUtil/GetImplicit3.hs b/Graphics/Implicit/ObjectUtil/GetImplicit3.hs
index 75159af5..1a5de593 100644
--- a/Graphics/Implicit/ObjectUtil/GetImplicit3.hs
+++ b/Graphics/Implicit/ObjectUtil/GetImplicit3.hs
@@ -6,7 +6,7 @@
 module Graphics.Implicit.ObjectUtil.GetImplicit3 (getImplicit3) where
 
 -- Import only what we need from the prelude.
-import Prelude (abs, atan2, acos, cos, ceiling, either, error, floor, fromInteger, fromIntegral, id, max, min, minimum, negate, otherwise, pi, pure, round, show, snd, sin, sqrt, sum, (||), (/=), Either(Left, Right), (<), (<=), (<>), (>), (>=), (&&), (-), (/), (*), (+), (++), ($), (.), Bool(True, False), (==), (**), Num, Applicative, Int, (<$>), Eq)
+import Prelude (abs, atan2, cos, ceiling, either, error, floor, fromInteger, id, max, min, minimum, negate, otherwise, pi, pure, round, sin, sqrt, sum, (||), (/=), Either(Left, Right), (>=), (-), (/), (*), (+), (++), ($), (.), Bool(True, False), (==), (**), Num, Applicative, Int, (<$>))
 
 import Control.Lens ((^.))
 
@@ -14,7 +14,7 @@ import qualified Data.Either as Either (either)
 
 import Data.Foldable (toList)
 
-import Data.List (concatMap, genericIndex, length, minimumBy)
+import Data.List (concatMap, genericIndex, minimumBy)
 
 import Data.Map (fromListWith, lookup, Map)
 
@@ -27,9 +27,6 @@ import Data.Sequence (fromList, mapWithIndex)
 import Linear (V2(V2), V3(V3), _xy, _z, distance, dot)
 import qualified Linear (conjugate, inv44, normalizePoint, normalize, point, rotate, Metric)
 
--- The cross product.
-import Linear.V3 (cross)
-
 -- Matrix times column vector.
 import Linear.Matrix ((!*))
 
@@ -46,14 +43,13 @@ import Graphics.Implicit.Definitions
       ℝ2,
       ℝ,
       fromℕtoℝ,
-      fromℕ,
       toScaleFn,
       ℝ3 )
 
 -- For handling extrusion of 2D shapes to 3D.
 import {-# SOURCE #-} Graphics.Implicit.Primitives (getImplicit)
 
-import Graphics.Implicit.TriUtil (angleAt, distancePointToTriangle, findTriangle, normOfTriangle, ClosestFeature(..), Tri)
+import Graphics.Implicit.TriUtil (angleAt, distancePointToTriangle, findTriangle, normOfTriangle, ClosestFeature(FeatVertex1, FeatVertex2, FeatVertex3, FeatEdge12, FeatEdge13, FeatEdge23, FeatFace), Tri, Triangle)
 
 import Graphics.Implicit.MathUtil (rmax, rmaximum)
 
@@ -100,8 +96,8 @@ getImplicit3 _ (Polyhedron points tris) = \(point) ->
     unsignedDistanceAndTriangleClosestTo point = minimumBy (\((_,a),_) ((_,b),_) -> a `compare` b) $ featDistTriangles point
     featDistTriangles point = (\a -> (distancePointToTriangle point (findTriangle points a), a)) <$> tris
     firstPointOfTri (v1,_,_) = v1
-    pointOnOutside :: ℝ3 -> (ℝ3,ℝ3,ℝ3) -> (ℕ,ℕ,ℕ) -> ClosestFeature -> Bool
-    pointOnOutside point closestTriangle closestTri feature = (point - firstPointOfTri closestTriangle) `dot` (weighedNormish points closestTri feature) >= -eps
+    pointOnOutside :: ℝ3 -> Triangle -> Tri -> ClosestFeature -> Bool
+    pointOnOutside point closestTriangle closestTri feature = (point - firstPointOfTri closestTriangle) `dot` (weighedNormish closestTri feature) >= -eps
       where
         -- fudge factor.
         eps :: ℝ
@@ -112,7 +108,7 @@ getImplicit3 _ (Polyhedron points tris) = \(point) ->
         triByEdge = fromListWith (++) edgeTris
           where
             edgeTris = concatMap edgesOfTri $ toList $ mapWithIndex (,) triSeq
-            edgesOfTri :: (Int,(ℕ,ℕ,ℕ)) -> [((ℕ,ℕ),[Int])]
+            edgesOfTri :: (Int,Tri) -> [((ℕ,ℕ),[Int])]
             edgesOfTri (i,(p1,p2,p3)) = [(sortEdge p1 p2,[i]),(sortEdge p2 p3,[i]),(sortEdge p3 p1,[i])]
         sortEdge a b = (min a b, max a b)
         -- For each vertex, the tri indexes that contain that vertex: 
@@ -120,29 +116,26 @@ getImplicit3 _ (Polyhedron points tris) = \(point) ->
         triByVertex = fromListWith (++) vertexTris
           where
             vertexTris = concatMap vertexesOfTri $ toList $ mapWithIndex (,) triSeq
-            vertexesOfTri :: (Int,(ℕ,ℕ,ℕ)) -> [(ℕ,[Int])]
+            vertexesOfTri :: (Int,Tri) -> [(ℕ,[Int])]
             vertexesOfTri (i,(p1,p2,p3)) = [(p1,[i]),(p2,[i]),(p3,[i])]  
         -- Get the normalized average of a set of triangles, referred to by index.
-        averageNorm triangles triIndexes = Linear.normalize $ sum $ normOfTriangle . genericIndex triangles <$> triIndexes
-        weighedNormish :: [ℝ3] -> (ℕ,ℕ,ℕ) -> ClosestFeature -> ℝ3
-        weighedNormish points tri@(p1,p2,p3) closest
+        averageNorm triIndexes = Linear.normalize $ sum $ normOfTriangle . genericIndex triangles <$> triIndexes
+        weighedNormish :: Tri -> ClosestFeature -> ℝ3
+        weighedNormish (p1,p2,p3) closest
           | closest == FeatFace = normOfTriangle closestTriangle
-          | closest == FeatEdge12 = averageNorm triangles $ fromMaybe [] $ lookup (sortEdge p1 p2) triByEdge 
-          | closest == FeatEdge13 = averageNorm triangles $ fromMaybe [] $ lookup (sortEdge p1 p3) triByEdge 
-          | closest == FeatEdge23 = averageNorm triangles $ fromMaybe [] $ lookup (sortEdge p2 p3) triByEdge
+          | closest == FeatEdge12 = averageNorm $ fromMaybe [] $ lookup (sortEdge p1 p2) triByEdge 
+          | closest == FeatEdge13 = averageNorm $ fromMaybe [] $ lookup (sortEdge p1 p3) triByEdge 
+          | closest == FeatEdge23 = averageNorm $ fromMaybe [] $ lookup (sortEdge p2 p3) triByEdge
           | closest == FeatVertex1 = Linear.normalize $ sum $ angleWeighed (genericIndex points p1) <$> fromMaybe [] (lookup p1 triByVertex)
           | closest == FeatVertex2 = Linear.normalize $ sum $ angleWeighed (genericIndex points p2) <$> fromMaybe [] (lookup p2 triByVertex)
           | closest == FeatVertex3 = Linear.normalize $ sum $ angleWeighed (genericIndex points p3) <$> fromMaybe [] (lookup p3 triByVertex)
           | otherwise = normOfTriangle closestTriangle
-          where
-            closestTriangle = findTriangle points tri
         angleWeighed :: ℝ3 -> Int -> ℝ3
         angleWeighed vertex triNo = angleAt vertex triangle *^ normOfTriangle triangle
           where
             triangle = findTriangle points $ genericIndex tris triNo
         -- decompose our tris into triangles.
         triangles = findTriangle points <$> tris
-
 getImplicit3 _ (BoxFrame b e) = \p' ->
     let p@(V3 px py pz) = abs p' - b
         V3 qx qy qz = abs (p + pure e) - pure e
diff --git a/Graphics/Implicit/TriUtil.hs b/Graphics/Implicit/TriUtil.hs
index 10c6de5b..84e00648 100644
--- a/Graphics/Implicit/TriUtil.hs
+++ b/Graphics/Implicit/TriUtil.hs
@@ -15,7 +15,8 @@ module Graphics.Implicit.TriUtil (
   ClosestFeature(FeatFace,
                  FeatVertex1, FeatVertex2, FeatVertex3,
                  FeatEdge12,  FeatEdge13,  FeatEdge23),
-  Tri
+  Tri,
+  Triangle
   ) where
 
 import Prelude (acos, error, length, max, min, otherwise, show, (>), (<), (&&), (<=), (>=),($), (.), (/), (+), (-), (*), (==), (||), (<>), Bool, Eq)
@@ -46,6 +47,7 @@ data ClosestFeature = FeatVertex1 | FeatVertex2 | FeatVertex3 | FeatEdge12 | Fea
   deriving Eq
 
 -- FIXME: Make these indices correct by construction?
+-- | Reconstitune a Triangle from a Tri, and our points array.
 findTriangle :: [ℝ3] -> Tri -> Triangle
 findTriangle vertices (i1,i2,i3)
   | outOfRange i1 = error $ "bad vertex index(out of range): " <> show i1 <> "\n"
@@ -58,9 +60,11 @@ findTriangle vertices (i1,i2,i3)
     outOfRange :: ℕ -> Bool
     outOfRange v = v < 0 || length vertices <= fromℕ v
 
+-- | Find the normal of a given Triangle
 normOfTriangle :: Triangle -> ℝ3
 normOfTriangle (v1,v2,v3) = Linear.normalize $ (v2-v1) `cross` (v3-v1)
 
+-- find the angle of the corner of the triangle containing a given vertex.
 angleAt :: ℝ3 -> Triangle -> ℝ
 angleAt vertex (v1,v2,v3)
   | vertex == v1 = acos $ clamp $ Linear.normalize (v2-v1) `dot` Linear.normalize (v3-v1)

From fb2f678ed6ce3f0d58bf76642d54bfa9d4b6e69a Mon Sep 17 00:00:00 2001
From: Julia Longtin <julia.longtin@gmail.com>
Date: Sat, 14 Mar 2026 20:04:06 +0000
Subject: [PATCH 20/28] tiny fixes.

---
 CHANGELOG.md                                  |  1 +
 Graphics/Implicit/Export/SymbolicFormats.hs   |  2 +-
 .../Implicit/ExtOpenScad/Eval/Statement.hs    |  5 ++-
 Graphics/Implicit/ExtOpenScad/Primitives.hs   | 35 ++++++++++---------
 Graphics/Implicit/ObjectUtil/GetBox3.hs       | 16 +++++----
 Graphics/Implicit/ObjectUtil/GetImplicit3.hs  |  8 ++---
 Graphics/Implicit/Primitives.hs               |  2 +-
 7 files changed, 37 insertions(+), 32 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 37948a12..5e97cee0 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,6 +1,7 @@
 # Version [next](https://github.com/Haskell-Things/ImplicitCAD/compare/v0.4.1.0...master) (202Y-MM-DD)
 
 * ExtOpenScad interface changes
+  * Added convex `polyhedron()` support [#497](https://github.com/Haskell-Things/ImplicitCAD/Pull/497)
   * Added `projection(cut=true)` support [#448](https://github.com/Haskell-Things/ImplicitCAD/pull/448)
 
 * Haskell interface changes
diff --git a/Graphics/Implicit/Export/SymbolicFormats.hs b/Graphics/Implicit/Export/SymbolicFormats.hs
index 7d09e00c..1869389a 100644
--- a/Graphics/Implicit/Export/SymbolicFormats.hs
+++ b/Graphics/Implicit/Export/SymbolicFormats.hs
@@ -137,7 +137,7 @@ buildS3 _ (Cylinder h r1 r2) = callNaked "cylinder" ["r1 = " <> pretty (bf r1)
                                                     , pretty $ bf h
                                                     ] []
 
-buildS3 _ (Polyhedron points tris) = callNaked "polyhedron" ["points = [" <> (fold $ intersperse "," $ renderPoint <$> points) <> "] faces = [" <> (fold $ intersperse "," $ renderTri <$> tris) <> "]" ] []
+buildS3 _ (Polyhedron points tris) = callNaked "polyhedron" ["points = [" <> (fold $ intersperse "," $ renderPoint <$> points) <> "], faces = [" <> (fold $ intersperse "," $ renderTri <$> tris) <> "]" ] []
   where
     renderPoint (V3 v1 v2 v3) = "[" <> pretty (bf v1) <> "," <> pretty (bf v2) <> "," <> pretty (bf v3) <> "]"
     renderTri (n1,n2,n3) = "[" <> pretty (bℕ n1) <> "," <> pretty (bℕ n2) <> "," <> pretty (bℕ n3) <> "]"
diff --git a/Graphics/Implicit/ExtOpenScad/Eval/Statement.hs b/Graphics/Implicit/ExtOpenScad/Eval/Statement.hs
index 4955c71d..6c6fb95d 100644
--- a/Graphics/Implicit/ExtOpenScad/Eval/Statement.hs
+++ b/Graphics/Implicit/ExtOpenScad/Eval/Statement.hs
@@ -116,21 +116,20 @@ runStatementI (StatementI sourcePos (ModuleCall (Symbol name) argsExpr suite)) =
               fromMaybe (pure []) (fst argsMapped)
             Just (ONModule _ implementation forms) -> do
               possibleInstances <- selectInstances forms
-              let
               when (null possibleInstances) (do
                                                 errorC sourcePos $ "no instance of " <> name <> " found to match given parameters.\nInstances available:\n" <> pack (show (ONModule (Symbol name) implementation forms))
                                                 traverse_ ((`checkOptions` True) . Just) forms
                                             )
               -- Evaluate all of the arguments.
               evaluatedArgs <- evalArgs argsExpr
+              when (suite /= []) (errorC sourcePos $ "Suite provided, but module " <> name <> " does not accept one. Perhaps a missing semicolon?")
               -- Run the module.
               let
                 argsMapped = argMap evaluatedArgs $ implementation sourcePos
               for_ (pack <$> snd argsMapped) $ errorC sourcePos
-              fromMaybe (pure []) (fst argsMapped)
+              fromMaybe (pure []) $ fst argsMapped
             Just (ONModuleWithSuite _ implementation forms) -> do
               possibleInstances <- selectInstances forms
-              let
               when (null possibleInstances) $ do
                                               errorC sourcePos $ "no instance of " <> name <> " found to match given parameters.\nInstances available:\n" <> pack (show (ONModuleWithSuite (Symbol name) implementation forms))
                                               traverse_ ((`checkOptions` True) . Just) forms
diff --git a/Graphics/Implicit/ExtOpenScad/Primitives.hs b/Graphics/Implicit/ExtOpenScad/Primitives.hs
index e01fb849..76a15054 100644
--- a/Graphics/Implicit/ExtOpenScad/Primitives.hs
+++ b/Graphics/Implicit/ExtOpenScad/Primitives.hs
@@ -15,7 +15,7 @@
 -- Export one set containing all of the primitive modules.
 module Graphics.Implicit.ExtOpenScad.Primitives (primitiveModules) where
 
-import Prelude(any, concat, elem, foldr, head, mapM, (.), Either(Left, Right), Bool(True, False), Maybe(Just, Nothing), ($), pure, show, either, id, (-), (==), (&&), (<), (*), cos, sin, pi, (/), (>), const, uncurry, (/=), (||), not, null, fmap, (<>), otherwise, error, (<*>), (<$>))
+import Prelude(any, concat, elem, error, foldr, head, mapM, (.), Either(Left, Right), Bool(True, False), Maybe(Just, Nothing), ($), pure, show, either, id, (-), (==), (&&), (<), (*), cos, sin, pi, (/), (>), const, uncurry, (/=), (||), not, null, fmap, (<>), otherwise, (<*>), (<$>))
 
 import Graphics.Implicit.Definitions (ℝ, ℝ2, ℝ3, ℕ, SymbolicObj2, SymbolicObj3, ExtrudeMScale(C1), fromℕtoℝ, isScaleID)
 
@@ -42,14 +42,14 @@ import Data.Foldable (toList)
 
 import Data.List (genericIndex)
 
+import Data.Maybe (fromMaybe, isJust)
+
 import Data.Sequence (Seq, deleteAt, filter, fromList)
 import qualified Data.Sequence as DS (null)
 
 import Data.Text.Lazy (Text)
 import qualified Data.Text.Lazy as DTL (pack)
 
-import Data.Maybe (isJust)
-
 import Control.Lens ((^.))
 
 import Linear (_m33, cross, dot, M34, M44, V2(V2), V3(V3), V4(V4))
@@ -58,6 +58,8 @@ import Linear.Affine (qdA)
 
 default (ℝ)
 
+-- FIXME: `defaultTo` is used inconsistently. The line between defaults and examples is a bit blurry.
+
 -- | Use the old syntax when defining arguments.
 argument :: OTypeMirror desiredType => Text -> ArgParser desiredType
 argument a = GIEUA.argument (Symbol a)
@@ -292,10 +294,10 @@ cylinder = moduleWithoutSuite "cylinder" $ \_ -> do
 
 polyhedron :: (Symbol, SourcePosition -> ArgParser (StateC [OVal]))
 polyhedron = moduleWithoutSuite "polyhedron" $ \sourcePos -> do
-    example "polyhedron(points=[[0,0,0], [2,0,0], [2,2,0], [0,2,0], [1, 1, 2]], faces=[[0,1,2,3], [0,5,1], [1,5,2], [2,5,3], [3,5,4], [4,5,0]]);"
+    example "polyhedron(points=[[0,0,0], [2,0,0], [2,2,0], [0,2,0], [1, 1, 2]], faces=[[0,1,2,3], [0,4,1], [1,4,2], [2,4,3], [3,4,0]]);"
     -- arguments
-    points :: [ℝ3] <- argument "points" `defaultTo` [] `doc` "list of points to construct faces from"
-    faces :: [[ℕ]] <- argument "faces" `defaultTo` [] `doc` "list of sets of indices into points, used to create faces on the polyhedron."
+    points :: [ℝ3] <- argument "points" `doc` "list of points to construct faces from"
+    faces :: [[ℕ]] <- argument "faces" `doc` "list of sets of indices into points, used to create faces on the polyhedron."
     pure $ do
       -- A tri is constructed of three indexes into the points.
       tris <- fmap concat $ mapM (trianglesFromFace sourcePos) faces
@@ -305,13 +307,13 @@ polyhedron = moduleWithoutSuite "polyhedron" $ \sourcePos -> do
         -- decompose our faces into tris.
         trianglesFromFace :: SourcePosition -> [ℕ] -> StateC [Tri]
         trianglesFromFace sourcePos []         = do
-                                                 errorC sourcePos "no point found when trying to generate triangles from a face.\n"
+                                                 warnC sourcePos "no point found when trying to generate triangles from a face.\n"
                                                  pure []
         trianglesFromFace sourcePos [p1]       = do
-                                                 warnC sourcePos $ "only one point found: " <> (DTL.pack $ show p1) <> "\n"
+                                                 errorC sourcePos $ "only one point found: " <> (DTL.pack $ show p1) <> "\n"
                                                  pure []
         trianglesFromFace sourcePos [p1,p2]    = do
-                                                 warnC sourcePos $ "only two points found: " <> (DTL.pack $ show p1) <> "\n" <> (DTL.pack $ show p2) <> "\n"
+                                                 errorC sourcePos $ "only two points found: " <> (DTL.pack $ show p1) <> "\n" <> (DTL.pack $ show p2) <> "\n"
                                                  pure []
         trianglesFromFace _         [p1,p2,p3] = pure [(p1,p2,p3)]
         trianglesFromFace sourcePos (p1:p2:p3:xs) = ((p1,p2,p3):) <$> trianglesFromFace sourcePos (p1:p3:xs)
@@ -357,7 +359,7 @@ polyhedron = moduleWithoutSuite "polyhedron" $ \sourcePos -> do
                   pure (found <> [triFound], filter (/= triUnderTest) remaining)
                 Nothing ->       pure (found, remaining)
               where
-                res = foldr (\tri state -> firstNeighborFilter tri triUnderTest state) Nothing visited 
+                res = foldr (\tri state -> firstNeighborFilter tri triUnderTest state) Nothing visited
                 -- | A short-circuiting filter we fold over visited, and grab the first neighboring tri.
                 firstNeighborFilter :: Tri -> Tri -> Maybe Tri -> Maybe Tri
                 firstNeighborFilter src testTri maybeRes
@@ -537,14 +539,13 @@ difference = moduleWithSuite "difference" $ \sourcePos children -> do
                     then objReduce (unsafeUncurry (Prim.differenceR r)) (unsafeUncurry (Prim.differenceR r)) children
                     else objReduce (unsafeUncurry  Prim.difference)     (unsafeUncurry  Prim.difference)     children
         where
-          unsafeUncurry :: (a -> [a] -> c) -> [a] -> c
-          unsafeUncurry f = uncurry f . unsafeUncons
-
-          unsafeUncons :: [a] -> (a, [a])
-          unsafeUncons (a : as) = (a, as)
+          unsafeUncons :: [a] -> Maybe (a, [a])
+          unsafeUncons (a : as) = Just (a, as)
           -- NOTE: This error is guarded against during the @null children@ check in the function body.
-          unsafeUncons _ = error "difference requires at least one element; zero given"
-
+          unsafeUncons _ = Nothing
+          -- This error should not be reachable.
+          unsafeUncurry :: (a -> [a] -> c) -> [a] -> c
+          unsafeUncurry f = uncurry f . fromMaybe (error "Impossible error: difference requires at least one element; zero given") .  unsafeUncons
 translate :: (Symbol, SourcePosition -> [OVal] -> ArgParser (StateC [OVal]))
 translate = moduleWithSuite "translate" $ \_ children -> do
     example "translate ([2,3]) circle (4);"
diff --git a/Graphics/Implicit/ObjectUtil/GetBox3.hs b/Graphics/Implicit/ObjectUtil/GetBox3.hs
index 472902eb..92a7a06b 100644
--- a/Graphics/Implicit/ObjectUtil/GetBox3.hs
+++ b/Graphics/Implicit/ObjectUtil/GetBox3.hs
@@ -5,11 +5,13 @@
 
 module Graphics.Implicit.ObjectUtil.GetBox3 (getBox3) where
 
-import Prelude(foldl, uncurry, pure, Bool(False), Either (Left, Right), (==), max, (/), (-), (+), fmap, unzip, ($), (<$>), (.), minimum, maximum, min, (>), (*), (<), abs, either, const, otherwise, take, fst, snd)
+import Prelude(uncurry, pure, Bool(False), Either (Left, Right), (==), max, (/), (-), (+), fmap, unzip, ($), (<$>), (.), minimum, maximum, min, (>), (*), (<), abs, either, const, otherwise, take, fst, snd)
 
 -- For Maybe types.
 import Data.Maybe (fromMaybe, Maybe(Just, Nothing))
 
+import Data.Foldable (foldl')
+
 import Linear (V2(V2), V3(V3))
 import qualified Linear (rotate, point, normalizePoint, (!*))
 
@@ -20,6 +22,8 @@ import Graphics.Implicit.Definitions
       SymbolicObj3(Shared3, Cube, Sphere, Cylinder, Polyhedron, Rotate3, Transform3, Extrude, ExtrudeOnEdgeOf, ExtrudeM, RotateExtrude, Torus, Ellipsoid, BoxFrame, Link),
       Box3,
       ℝ,
+      ℝ2,
+      ℝ3,
       fromFastℕtoℝ,
       toScaleFn )
 
@@ -39,10 +43,10 @@ getBox3 (Cylinder h r1 r2) = (V3 (-r) (-r) 0, V3 r r h ) where r = max r1 r2
 getBox3 (Polyhedron points _) = (minimum_point, maximum_point)
   where
     (minimum_point, maximum_point) = fromMaybe (V3 0 0 0, V3 0 0 0) maybeVs
-    maybeVs :: (Maybe (V3 ℝ,V3 ℝ))
-    maybeVs = foldl findMinMax Nothing points
+    maybeVs :: (Maybe (ℝ3,ℝ3))
+    maybeVs = foldl' findMinMax Nothing points
       where
-        findMinMax :: (Maybe (V3 ℝ,V3 ℝ)) -> V3 ℝ -> (Maybe (V3 ℝ,V3 ℝ))
+        findMinMax :: (Maybe (ℝ3,ℝ3)) -> ℝ3 -> (Maybe (ℝ3,ℝ3))
         findMinMax Nothing newV3 = Just (newV3, newV3)
         findMinMax (Just (V3 minx miny minz,V3 maxx maxy maxz)) (V3 newx newy newz) = Just (V3 (min minx newx) (min miny newy) (min minz newz), V3 (max maxx newx) (max maxy newy) (max maxz newz))
 getBox3 (Torus r1 r2) =
@@ -117,9 +121,9 @@ getBox3 (ExtrudeM twist scale translate symbObj height) =
 
         (tminx, tmaxx, tminy, tmaxy) =
           let
-            tvalsx :: (ℝ -> V2 ℝ) -> [ℝ]
+            tvalsx :: (ℝ -> ℝ2) -> [ℝ]
             tvalsx tfun = fmap (fst . unpack . tfun) hrange
-            tvalsy :: (ℝ -> V2 ℝ) -> [ℝ]
+            tvalsy :: (ℝ -> ℝ2) -> [ℝ]
             tvalsy tfun = fmap (snd . unpack . tfun) hrange
           in case translate of
             Left  (V2 tvalx tvaly) -> (tvalx, tvalx, tvaly, tvaly)
diff --git a/Graphics/Implicit/ObjectUtil/GetImplicit3.hs b/Graphics/Implicit/ObjectUtil/GetImplicit3.hs
index 1a5de593..84a90137 100644
--- a/Graphics/Implicit/ObjectUtil/GetImplicit3.hs
+++ b/Graphics/Implicit/ObjectUtil/GetImplicit3.hs
@@ -103,7 +103,7 @@ getImplicit3 _ (Polyhedron points tris) = \(point) ->
         eps :: ℝ
         eps = 1e-13
         triSeq = fromList tris
-        -- For each edge, the tri indexes that share that edge: 
+        -- For each edge, the tri indexes that share that edge:
         triByEdge :: Map (ℕ,ℕ) [Int]
         triByEdge = fromListWith (++) edgeTris
           where
@@ -111,7 +111,7 @@ getImplicit3 _ (Polyhedron points tris) = \(point) ->
             edgesOfTri :: (Int,Tri) -> [((ℕ,ℕ),[Int])]
             edgesOfTri (i,(p1,p2,p3)) = [(sortEdge p1 p2,[i]),(sortEdge p2 p3,[i]),(sortEdge p3 p1,[i])]
         sortEdge a b = (min a b, max a b)
-        -- For each vertex, the tri indexes that contain that vertex: 
+        -- For each vertex, the tri indexes that contain that vertex:
         triByVertex :: Map ℕ [Int]
         triByVertex = fromListWith (++) vertexTris
           where
@@ -123,8 +123,8 @@ getImplicit3 _ (Polyhedron points tris) = \(point) ->
         weighedNormish :: Tri -> ClosestFeature -> ℝ3
         weighedNormish (p1,p2,p3) closest
           | closest == FeatFace = normOfTriangle closestTriangle
-          | closest == FeatEdge12 = averageNorm $ fromMaybe [] $ lookup (sortEdge p1 p2) triByEdge 
-          | closest == FeatEdge13 = averageNorm $ fromMaybe [] $ lookup (sortEdge p1 p3) triByEdge 
+          | closest == FeatEdge12 = averageNorm $ fromMaybe [] $ lookup (sortEdge p1 p2) triByEdge
+          | closest == FeatEdge13 = averageNorm $ fromMaybe [] $ lookup (sortEdge p1 p3) triByEdge
           | closest == FeatEdge23 = averageNorm $ fromMaybe [] $ lookup (sortEdge p2 p3) triByEdge
           | closest == FeatVertex1 = Linear.normalize $ sum $ angleWeighed (genericIndex points p1) <$> fromMaybe [] (lookup p1 triByVertex)
           | closest == FeatVertex2 = Linear.normalize $ sum $ angleWeighed (genericIndex points p2) <$> fromMaybe [] (lookup p2 triByVertex)
diff --git a/Graphics/Implicit/Primitives.hs b/Graphics/Implicit/Primitives.hs
index 9ddb6447..8bdaaf28 100644
--- a/Graphics/Implicit/Primitives.hs
+++ b/Graphics/Implicit/Primitives.hs
@@ -138,7 +138,7 @@ polyhedron
     :: [ℝ3]           -- ^ Points
     -> [(ℕ,ℕ,ℕ)]      -- ^ triangles, resolved through indexing Points
     -> SymbolicObj3   -- ^ Resulting polyhedron
-polyhedron points tris = Polyhedron points tris
+polyhedron = Polyhedron
 
 -- | A conical frustum --- ie. a cylinder with different radii at either end.
 cylinder2

From e53aef2706b9851dde5187d6e16532eeab564437 Mon Sep 17 00:00:00 2001
From: Julia Longtin <julia.longtin@gmail.com>
Date: Sun, 15 Mar 2026 14:10:11 +0000
Subject: [PATCH 21/28] move items out of pointOnOutside, so they aren't
 recomputed during the lambda.

---
 Graphics/Implicit/ExtOpenScad/Primitives.hs  |  2 +-
 Graphics/Implicit/ObjectUtil/GetImplicit3.hs | 39 ++++++++++----------
 2 files changed, 21 insertions(+), 20 deletions(-)

diff --git a/Graphics/Implicit/ExtOpenScad/Primitives.hs b/Graphics/Implicit/ExtOpenScad/Primitives.hs
index 76a15054..c17101d5 100644
--- a/Graphics/Implicit/ExtOpenScad/Primitives.hs
+++ b/Graphics/Implicit/ExtOpenScad/Primitives.hs
@@ -335,7 +335,7 @@ polyhedron = moduleWithoutSuite "polyhedron" $ \sourcePos -> do
                 triCentroid = centroid [v1,v2,v3]
             polyCentroid = centroid points
             flipTri (p1,p2,p3) = (p1,p3,p2)
-            -- | Wind the triangles.
+            -- | Wind the triangles. For two triangles sharing an edge, said edge MUST be expressed by the two triangles in the opposite order.
             windTriangles :: [Tri] -> Seq Tri -> StateC [Tri]
             windTriangles visited unvisited
               -- We're done.
diff --git a/Graphics/Implicit/ObjectUtil/GetImplicit3.hs b/Graphics/Implicit/ObjectUtil/GetImplicit3.hs
index 84a90137..78eeeddb 100644
--- a/Graphics/Implicit/ObjectUtil/GetImplicit3.hs
+++ b/Graphics/Implicit/ObjectUtil/GetImplicit3.hs
@@ -96,28 +96,31 @@ getImplicit3 _ (Polyhedron points tris) = \(point) ->
     unsignedDistanceAndTriangleClosestTo point = minimumBy (\((_,a),_) ((_,b),_) -> a `compare` b) $ featDistTriangles point
     featDistTriangles point = (\a -> (distancePointToTriangle point (findTriangle points a), a)) <$> tris
     firstPointOfTri (v1,_,_) = v1
+    -- For each edge, the tri indexes that share that edge:
+    triByEdge :: Map (ℕ,ℕ) [Int]
+    triByEdge = fromListWith (++) edgeTris
+      where
+        edgeTris = concatMap edgesOfTri $ toList $ mapWithIndex (,) triSeq
+        edgesOfTri :: (Int,Tri) -> [((ℕ,ℕ),[Int])]
+        edgesOfTri (i,(p1,p2,p3)) = [(sortEdge p1 p2,[i]),(sortEdge p2 p3,[i]),(sortEdge p3 p1,[i])]
+    -- For each vertex, the tri indexes that contain that vertex:
+    triByVertex :: Map ℕ [Int]
+    triByVertex = fromListWith (++) vertexTris
+      where
+        vertexTris = concatMap vertexesOfTri $ toList $ mapWithIndex (,) triSeq
+        vertexesOfTri :: (Int,Tri) -> [(ℕ,[Int])]
+        vertexesOfTri (i,(p1,p2,p3)) = [(p1,[i]),(p2,[i]),(p3,[i])]
+    -- Sequence our Tris.
+    triSeq = fromList tris
+    -- Decompose our tris into triangles.
+    triangles = findTriangle points <$> tris
+    sortEdge a b = (min a b, max a b)
     pointOnOutside :: ℝ3 -> Triangle -> Tri -> ClosestFeature -> Bool
     pointOnOutside point closestTriangle closestTri feature = (point - firstPointOfTri closestTriangle) `dot` (weighedNormish closestTri feature) >= -eps
       where
         -- fudge factor.
         eps :: ℝ
         eps = 1e-13
-        triSeq = fromList tris
-        -- For each edge, the tri indexes that share that edge:
-        triByEdge :: Map (ℕ,ℕ) [Int]
-        triByEdge = fromListWith (++) edgeTris
-          where
-            edgeTris = concatMap edgesOfTri $ toList $ mapWithIndex (,) triSeq
-            edgesOfTri :: (Int,Tri) -> [((ℕ,ℕ),[Int])]
-            edgesOfTri (i,(p1,p2,p3)) = [(sortEdge p1 p2,[i]),(sortEdge p2 p3,[i]),(sortEdge p3 p1,[i])]
-        sortEdge a b = (min a b, max a b)
-        -- For each vertex, the tri indexes that contain that vertex:
-        triByVertex :: Map ℕ [Int]
-        triByVertex = fromListWith (++) vertexTris
-          where
-            vertexTris = concatMap vertexesOfTri $ toList $ mapWithIndex (,) triSeq
-            vertexesOfTri :: (Int,Tri) -> [(ℕ,[Int])]
-            vertexesOfTri (i,(p1,p2,p3)) = [(p1,[i]),(p2,[i]),(p3,[i])]  
         -- Get the normalized average of a set of triangles, referred to by index.
         averageNorm triIndexes = Linear.normalize $ sum $ normOfTriangle . genericIndex triangles <$> triIndexes
         weighedNormish :: Tri -> ClosestFeature -> ℝ3
@@ -133,9 +136,7 @@ getImplicit3 _ (Polyhedron points tris) = \(point) ->
         angleWeighed :: ℝ3 -> Int -> ℝ3
         angleWeighed vertex triNo = angleAt vertex triangle *^ normOfTriangle triangle
           where
-            triangle = findTriangle points $ genericIndex tris triNo
-        -- decompose our tris into triangles.
-        triangles = findTriangle points <$> tris
+            triangle = genericIndex triangles triNo
 getImplicit3 _ (BoxFrame b e) = \p' ->
     let p@(V3 px py pz) = abs p' - b
         V3 qx qy qz = abs (p + pure e) - pure e

From f840ed1f5cb2b33d921adba12dc2bf45fa94bff4 Mon Sep 17 00:00:00 2001
From: Julia Longtin <julia.longtin@gmail.com>
Date: Sun, 15 Mar 2026 21:57:22 +0000
Subject: [PATCH 22/28] provide function for finding insideness based on
 winding numbers.

---
 Graphics/Implicit/TriUtil.hs | 35 +++++++++++++++++++++++++++++++++--
 1 file changed, 33 insertions(+), 2 deletions(-)

diff --git a/Graphics/Implicit/TriUtil.hs b/Graphics/Implicit/TriUtil.hs
index 84e00648..9abc6761 100644
--- a/Graphics/Implicit/TriUtil.hs
+++ b/Graphics/Implicit/TriUtil.hs
@@ -12,6 +12,7 @@ module Graphics.Implicit.TriUtil (
   distancePointToTriangle,
   findTriangle,
   normOfTriangle,
+  pointOnOutsideByWinding,
   ClosestFeature(FeatFace,
                  FeatVertex1, FeatVertex2, FeatVertex3,
                  FeatEdge12,  FeatEdge13,  FeatEdge23),
@@ -19,7 +20,7 @@ module Graphics.Implicit.TriUtil (
   Triangle
   ) where
 
-import Prelude (acos, error, length, max, min, otherwise, show, (>), (<), (&&), (<=), (>=),($), (.), (/), (+), (-), (*), (==), (||), (<>), Bool, Eq)
+import Prelude (abs, acos, atan2, error, length, max, min, otherwise, pi, show, sum, (<$>), (>), (<), (&&), (<=), (>=),($), (.), (/), (+), (-), (*), (==), (||), (<>), Bool, Eq)
 
 import Graphics.Implicit.Definitions (
   fromℕ,
@@ -29,7 +30,7 @@ import Graphics.Implicit.Definitions (
 
 import Data.List (genericIndex)
 
-import Linear (dot, distance)
+import Linear (dot, distance, norm)
 import qualified Linear (normalize)
 
 -- The cross product.
@@ -64,6 +65,36 @@ findTriangle vertices (i1,i2,i3)
 normOfTriangle :: Triangle -> ℝ3
 normOfTriangle (v1,v2,v3) = Linear.normalize $ (v2-v1) `cross` (v3-v1)
 
+-- | Determine if the point is on the outside of the object.
+-- Note: We cannot terminate early.. because all of the triangles must be counted.
+pointOnOutsideByWinding :: ℝ3 -> [Triangle] -> Bool
+pointOnOutsideByWinding point triangles = abs res < pi
+  where
+    res = sum $ halfSolidAngleOfTriangle point <$> triangles
+
+-- Find half of the signed solid angle of a given Triangle as seen from point.
+-- A solid angle measures how much of the field of view this triangle covers. think: arcs in astronomy.
+halfSolidAngleOfTriangle :: ℝ3 -> Triangle -> ℝ
+halfSolidAngleOfTriangle point (v1,v2,v3)
+  | normv1p <= eps || normv2p <= eps || normv3p <= eps = 0
+  | otherwise                                          = atan2 numerator denominator
+    where
+      eps :: ℝ
+      eps = 1e-14
+      v1p = v1 - point
+      v2p = v2 - point
+      v3p = v3 - point
+      -- length of v1p, v2p, v3p
+      normv1p = norm v1p
+      normv2p = norm v2p
+      normv3p = norm v3p
+      -- The Van Oosterom-Strackee Algorithm.
+      numerator = v1p `dot` (v2p `cross` v3p)
+      denominator = normv1p * normv2p * normv3p
+                  + (v1p `dot` v2p) * normv3p
+                  + (v2p `dot` v3p) * normv1p
+                  + (v3p `dot` v1p) * normv2p
+
 -- find the angle of the corner of the triangle containing a given vertex.
 angleAt :: ℝ3 -> Triangle -> ℝ
 angleAt vertex (v1,v2,v3)

From 92612de61ba908b50f9266701f02a91a51a54360 Mon Sep 17 00:00:00 2001
From: Julia Longtin <julia.longtin@gmail.com>
Date: Sun, 15 Mar 2026 22:49:44 +0000
Subject: [PATCH 23/28] use pointOnOutsideByWinding.

---
 Graphics/Implicit/ObjectUtil/GetImplicit3.hs | 66 +++-----------------
 1 file changed, 8 insertions(+), 58 deletions(-)

diff --git a/Graphics/Implicit/ObjectUtil/GetImplicit3.hs b/Graphics/Implicit/ObjectUtil/GetImplicit3.hs
index 78eeeddb..5c7075f0 100644
--- a/Graphics/Implicit/ObjectUtil/GetImplicit3.hs
+++ b/Graphics/Implicit/ObjectUtil/GetImplicit3.hs
@@ -6,33 +6,22 @@
 module Graphics.Implicit.ObjectUtil.GetImplicit3 (getImplicit3) where
 
 -- Import only what we need from the prelude.
-import Prelude (abs, atan2, cos, ceiling, either, error, floor, fromInteger, id, max, min, minimum, negate, otherwise, pi, pure, round, sin, sqrt, sum, (||), (/=), Either(Left, Right), (>=), (-), (/), (*), (+), (++), ($), (.), Bool(True, False), (==), (**), Num, Applicative, Int, (<$>))
+import Prelude (abs, atan2, cos, ceiling, either, error, floor, fromInteger, id, max, min, minimum, negate, otherwise, pi, pure, round, sin, sqrt, (||), (/=), Either(Left, Right), (-), (/), (*), (+), ($), (.), Bool(True, False), (==), (**), Num, Applicative, (<$>))
 
 import Control.Lens ((^.))
 
 import qualified Data.Either as Either (either)
 
-import Data.Foldable (toList)
-
-import Data.List (concatMap, genericIndex, minimumBy)
-
-import Data.Map (fromListWith, lookup, Map)
-
-import Data.Maybe (fromMaybe)
+import Data.List (minimumBy)
 
 import Data.Ord (compare)
 
-import Data.Sequence (fromList, mapWithIndex)
-
-import Linear (V2(V2), V3(V3), _xy, _z, distance, dot)
-import qualified Linear (conjugate, inv44, normalizePoint, normalize, point, rotate, Metric)
+import Linear (V2(V2), V3(V3), _xy, _z, distance)
+import qualified Linear (conjugate, inv44, normalizePoint, point, rotate, Metric)
 
 -- Matrix times column vector.
 import Linear.Matrix ((!*))
 
--- Matrix times scalar.
-import Linear.Vector ((*^))
-
 import Graphics.Implicit.Definitions
     ( objectRounding,
       ObjectContext,
@@ -49,7 +38,7 @@ import Graphics.Implicit.Definitions
 -- For handling extrusion of 2D shapes to 3D.
 import {-# SOURCE #-} Graphics.Implicit.Primitives (getImplicit)
 
-import Graphics.Implicit.TriUtil (angleAt, distancePointToTriangle, findTriangle, normOfTriangle, ClosestFeature(FeatVertex1, FeatVertex2, FeatVertex3, FeatEdge12, FeatEdge13, FeatEdge23, FeatFace), Tri, Triangle)
+import Graphics.Implicit.TriUtil (distancePointToTriangle, findTriangle, pointOnOutsideByWinding)
 
 import Graphics.Implicit.MathUtil (rmax, rmaximum)
 
@@ -87,56 +76,17 @@ getImplicit3 _ (Polyhedron [] _) = error "Asked to find distance to an empty pol
 getImplicit3 _ (Polyhedron _ []) = error "Asked to find distance to an empty polygon. No tris."
 getImplicit3 _ (Polyhedron points tris) = \(point) ->
   let
-    ((feature,res), closestTri) = unsignedDistanceAndTriangleClosestTo point
+    ((_,res), _) = unsignedDistanceAndTriangleClosestTo point
   in
-    if pointOnOutside point (findTriangle points closestTri) closestTri feature
+--    if pointOnOutside point (findTriangle points closestTri) closestTri feature
+    if pointOnOutsideByWinding point triangles
     then          res
     else negate $ res
   where
     unsignedDistanceAndTriangleClosestTo point = minimumBy (\((_,a),_) ((_,b),_) -> a `compare` b) $ featDistTriangles point
     featDistTriangles point = (\a -> (distancePointToTriangle point (findTriangle points a), a)) <$> tris
-    firstPointOfTri (v1,_,_) = v1
-    -- For each edge, the tri indexes that share that edge:
-    triByEdge :: Map (ℕ,ℕ) [Int]
-    triByEdge = fromListWith (++) edgeTris
-      where
-        edgeTris = concatMap edgesOfTri $ toList $ mapWithIndex (,) triSeq
-        edgesOfTri :: (Int,Tri) -> [((ℕ,ℕ),[Int])]
-        edgesOfTri (i,(p1,p2,p3)) = [(sortEdge p1 p2,[i]),(sortEdge p2 p3,[i]),(sortEdge p3 p1,[i])]
-    -- For each vertex, the tri indexes that contain that vertex:
-    triByVertex :: Map ℕ [Int]
-    triByVertex = fromListWith (++) vertexTris
-      where
-        vertexTris = concatMap vertexesOfTri $ toList $ mapWithIndex (,) triSeq
-        vertexesOfTri :: (Int,Tri) -> [(ℕ,[Int])]
-        vertexesOfTri (i,(p1,p2,p3)) = [(p1,[i]),(p2,[i]),(p3,[i])]
-    -- Sequence our Tris.
-    triSeq = fromList tris
     -- Decompose our tris into triangles.
     triangles = findTriangle points <$> tris
-    sortEdge a b = (min a b, max a b)
-    pointOnOutside :: ℝ3 -> Triangle -> Tri -> ClosestFeature -> Bool
-    pointOnOutside point closestTriangle closestTri feature = (point - firstPointOfTri closestTriangle) `dot` (weighedNormish closestTri feature) >= -eps
-      where
-        -- fudge factor.
-        eps :: ℝ
-        eps = 1e-13
-        -- Get the normalized average of a set of triangles, referred to by index.
-        averageNorm triIndexes = Linear.normalize $ sum $ normOfTriangle . genericIndex triangles <$> triIndexes
-        weighedNormish :: Tri -> ClosestFeature -> ℝ3
-        weighedNormish (p1,p2,p3) closest
-          | closest == FeatFace = normOfTriangle closestTriangle
-          | closest == FeatEdge12 = averageNorm $ fromMaybe [] $ lookup (sortEdge p1 p2) triByEdge
-          | closest == FeatEdge13 = averageNorm $ fromMaybe [] $ lookup (sortEdge p1 p3) triByEdge
-          | closest == FeatEdge23 = averageNorm $ fromMaybe [] $ lookup (sortEdge p2 p3) triByEdge
-          | closest == FeatVertex1 = Linear.normalize $ sum $ angleWeighed (genericIndex points p1) <$> fromMaybe [] (lookup p1 triByVertex)
-          | closest == FeatVertex2 = Linear.normalize $ sum $ angleWeighed (genericIndex points p2) <$> fromMaybe [] (lookup p2 triByVertex)
-          | closest == FeatVertex3 = Linear.normalize $ sum $ angleWeighed (genericIndex points p3) <$> fromMaybe [] (lookup p3 triByVertex)
-          | otherwise = normOfTriangle closestTriangle
-        angleWeighed :: ℝ3 -> Int -> ℝ3
-        angleWeighed vertex triNo = angleAt vertex triangle *^ normOfTriangle triangle
-          where
-            triangle = genericIndex triangles triNo
 getImplicit3 _ (BoxFrame b e) = \p' ->
     let p@(V3 px py pz) = abs p' - b
         V3 qx qy qz = abs (p + pure e) - pure e

From 90d819159d12473194a9110eba418b410c27f78c Mon Sep 17 00:00:00 2001
From: Julia Longtin <julia.longtin@gmail.com>
Date: Sun, 15 Mar 2026 23:00:12 +0000
Subject: [PATCH 24/28] should work on concave or convex

---
 CHANGELOG.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 5e97cee0..41371535 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,7 +1,7 @@
 # Version [next](https://github.com/Haskell-Things/ImplicitCAD/compare/v0.4.1.0...master) (202Y-MM-DD)
 
 * ExtOpenScad interface changes
-  * Added convex `polyhedron()` support [#497](https://github.com/Haskell-Things/ImplicitCAD/Pull/497)
+  * Added `polyhedron()` support [#497](https://github.com/Haskell-Things/ImplicitCAD/Pull/497)
   * Added `projection(cut=true)` support [#448](https://github.com/Haskell-Things/ImplicitCAD/pull/448)
 
 * Haskell interface changes

From a93ac47e27944c37b61980915a331561979d4af2 Mon Sep 17 00:00:00 2001
From: Julia Longtin <julia.longtin@gmail.com>
Date: Mon, 16 Mar 2026 15:55:58 +0000
Subject: [PATCH 25/28] tiny changes.

---
 CHANGELOG.md                                 | 2 +-
 Graphics/Implicit/Definitions.hs             | 2 +-
 Graphics/Implicit/ObjectUtil/GetImplicit3.hs | 4 ++--
 Graphics/Implicit/TriUtil.hs                 | 2 +-
 4 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 41371535..426d319d 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,7 +1,7 @@
 # Version [next](https://github.com/Haskell-Things/ImplicitCAD/compare/v0.4.1.0...master) (202Y-MM-DD)
 
 * ExtOpenScad interface changes
-  * Added `polyhedron()` support [#497](https://github.com/Haskell-Things/ImplicitCAD/Pull/497)
+  * Added `polyhedron()` support [#497](https://github.com/Haskell-Things/ImplicitCAD/pull/497)
   * Added `projection(cut=true)` support [#448](https://github.com/Haskell-Things/ImplicitCAD/pull/448)
 
 * Haskell interface changes
diff --git a/Graphics/Implicit/Definitions.hs b/Graphics/Implicit/Definitions.hs
index 4b2282d0..73f87aa6 100644
--- a/Graphics/Implicit/Definitions.hs
+++ b/Graphics/Implicit/Definitions.hs
@@ -332,7 +332,7 @@ data SymbolicObj3 =
       Cube ℝ3 -- rounding, size.
     | Sphere ℝ -- radius
     | Cylinder ℝ ℝ ℝ --
-    | Polyhedron [ℝ3] [(ℕ,ℕ,ℕ)] -- virtexes, triangles-by-index
+    | Polyhedron [ℝ3] [(ℕ,ℕ,ℕ)] -- vertexes, triangles-by-index
     -- Simple transforms
     | Rotate3 (Quaternion ℝ) SymbolicObj3
     | Transform3 (M44 ℝ) SymbolicObj3
diff --git a/Graphics/Implicit/ObjectUtil/GetImplicit3.hs b/Graphics/Implicit/ObjectUtil/GetImplicit3.hs
index 5c7075f0..3d89dae7 100644
--- a/Graphics/Implicit/ObjectUtil/GetImplicit3.hs
+++ b/Graphics/Implicit/ObjectUtil/GetImplicit3.hs
@@ -72,8 +72,8 @@ getImplicit3 _ (Cylinder h r1 r2) = \(V3 x y z) ->
     in
         max (d * cos θ) (abs (z-h/2) - (h/2))
 -- FIXME: Make Polyhedron correct by construction.
-getImplicit3 _ (Polyhedron [] _) = error "Asked to find distance to an empty polygon. No points."
-getImplicit3 _ (Polyhedron _ []) = error "Asked to find distance to an empty polygon. No tris."
+getImplicit3 _ (Polyhedron [] _) = error "Asked to find distance to an empty polyhedron. No points given."
+getImplicit3 _ (Polyhedron _ []) = error "Asked to find distance to an empty polyhedron. No faces given."
 getImplicit3 _ (Polyhedron points tris) = \(point) ->
   let
     ((_,res), _) = unsignedDistanceAndTriangleClosestTo point
diff --git a/Graphics/Implicit/TriUtil.hs b/Graphics/Implicit/TriUtil.hs
index 9abc6761..c2ec3a2f 100644
--- a/Graphics/Implicit/TriUtil.hs
+++ b/Graphics/Implicit/TriUtil.hs
@@ -48,7 +48,7 @@ data ClosestFeature = FeatVertex1 | FeatVertex2 | FeatVertex3 | FeatEdge12 | Fea
   deriving Eq
 
 -- FIXME: Make these indices correct by construction?
--- | Reconstitune a Triangle from a Tri, and our points array.
+-- | Reconstitute a Triangle from a Tri, and our points array.
 findTriangle :: [ℝ3] -> Tri -> Triangle
 findTriangle vertices (i1,i2,i3)
   | outOfRange i1 = error $ "bad vertex index(out of range): " <> show i1 <> "\n"

From fb80e97fcbdf589ad031dc7c4ffeb3515a1f0f2c Mon Sep 17 00:00:00 2001
From: Julia Longtin <julia.longtin@gmail.com>
Date: Tue, 17 Mar 2026 21:52:26 +0000
Subject: [PATCH 26/28] use the triangles from triangles, rather than tris.

---
 Graphics/Implicit/ObjectUtil/GetImplicit3.hs | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Graphics/Implicit/ObjectUtil/GetImplicit3.hs b/Graphics/Implicit/ObjectUtil/GetImplicit3.hs
index 3d89dae7..5a1707bf 100644
--- a/Graphics/Implicit/ObjectUtil/GetImplicit3.hs
+++ b/Graphics/Implicit/ObjectUtil/GetImplicit3.hs
@@ -84,7 +84,7 @@ getImplicit3 _ (Polyhedron points tris) = \(point) ->
     else negate $ res
   where
     unsignedDistanceAndTriangleClosestTo point = minimumBy (\((_,a),_) ((_,b),_) -> a `compare` b) $ featDistTriangles point
-    featDistTriangles point = (\a -> (distancePointToTriangle point (findTriangle points a), a)) <$> tris
+    featDistTriangles point = (\triangle -> (distancePointToTriangle point triangle, triangle)) <$> triangles
     -- Decompose our tris into triangles.
     triangles = findTriangle points <$> tris
 getImplicit3 _ (BoxFrame b e) = \p' ->

From 30f2ccf69e5cd24c945e8147929df43f5af4c4cd Mon Sep 17 00:00:00 2001
From: Julia Longtin <julia.longtin@gmail.com>
Date: Thu, 19 Mar 2026 18:52:13 +0000
Subject: [PATCH 27/28] remove unnecessary import.

---
 Graphics/Implicit/ObjectUtil/GetImplicit3.hs | 6 ++----
 1 file changed, 2 insertions(+), 4 deletions(-)

diff --git a/Graphics/Implicit/ObjectUtil/GetImplicit3.hs b/Graphics/Implicit/ObjectUtil/GetImplicit3.hs
index 5a1707bf..114575a9 100644
--- a/Graphics/Implicit/ObjectUtil/GetImplicit3.hs
+++ b/Graphics/Implicit/ObjectUtil/GetImplicit3.hs
@@ -10,8 +10,6 @@ import Prelude (abs, atan2, cos, ceiling, either, error, floor, fromInteger, id,
 
 import Control.Lens ((^.))
 
-import qualified Data.Either as Either (either)
-
 import Data.List (minimumBy)
 
 import Data.Ord (compare)
@@ -175,12 +173,12 @@ getImplicit3 ctx (RotateExtrude totalRotation translate rotate symbObj) =
             || either is360m (\f -> is360m (f 0 - f totalRotation)) rotate
         round' = objectRounding ctx
         translate' :: ℝ -> ℝ2
-        translate' = Either.either
+        translate' = either
                 (\(V2 a b) θ -> V2 (a*θ/totalRotation) (b*θ/totalRotation))
                 id
                 translate
         rotate' :: ℝ -> ℝ
-        rotate' = Either.either
+        rotate' = either
                 (\t θ -> t*θ/totalRotation )
                 id
                 rotate

From 3e227a472c0a32fff02c4001a16182954eb41b93 Mon Sep 17 00:00:00 2001
From: Julia Longtin <julia.longtin@gmail.com>
Date: Thu, 19 Mar 2026 18:57:44 +0000
Subject: [PATCH 28/28] add a note, and make sure we can handle polygons with
 more than 4 billion faces. ;)

---
 Graphics/Implicit/ExtOpenScad/Primitives.hs | 3 ++-
 Graphics/Implicit/TriUtil.hs                | 4 ++--
 2 files changed, 4 insertions(+), 3 deletions(-)

diff --git a/Graphics/Implicit/ExtOpenScad/Primitives.hs b/Graphics/Implicit/ExtOpenScad/Primitives.hs
index c17101d5..bd824892 100644
--- a/Graphics/Implicit/ExtOpenScad/Primitives.hs
+++ b/Graphics/Implicit/ExtOpenScad/Primitives.hs
@@ -295,7 +295,8 @@ cylinder = moduleWithoutSuite "cylinder" $ \_ -> do
 polyhedron :: (Symbol, SourcePosition -> ArgParser (StateC [OVal]))
 polyhedron = moduleWithoutSuite "polyhedron" $ \sourcePos -> do
     example "polyhedron(points=[[0,0,0], [2,0,0], [2,2,0], [0,2,0], [1, 1, 2]], faces=[[0,1,2,3], [0,4,1], [1,4,2], [2,4,3], [3,4,0]]);"
-    -- arguments
+    -- Arguments
+    -- FIXME: find a way to mark an arguement as non-empty!
     points :: [ℝ3] <- argument "points" `doc` "list of points to construct faces from"
     faces :: [[ℕ]] <- argument "faces" `doc` "list of sets of indices into points, used to create faces on the polyhedron."
     pure $ do
diff --git a/Graphics/Implicit/TriUtil.hs b/Graphics/Implicit/TriUtil.hs
index c2ec3a2f..cd80fc0d 100644
--- a/Graphics/Implicit/TriUtil.hs
+++ b/Graphics/Implicit/TriUtil.hs
@@ -20,7 +20,7 @@ module Graphics.Implicit.TriUtil (
   Triangle
   ) where
 
-import Prelude (abs, acos, atan2, error, length, max, min, otherwise, pi, show, sum, (<$>), (>), (<), (&&), (<=), (>=),($), (.), (/), (+), (-), (*), (==), (||), (<>), Bool, Eq)
+import Prelude (abs, acos, atan2, error, length, max, min, otherwise, pi, show, sum, toInteger, (<$>), (>), (<), (&&), (<=), (>=),($), (.), (/), (+), (-), (*), (==), (||), (<>), Bool, Eq)
 
 import Graphics.Implicit.Definitions (
   fromℕ,
@@ -59,7 +59,7 @@ findTriangle vertices (i1,i2,i3)
   where
     -- FIXME: >=BASE-4.21: replace this with compareLength once debian stable ships 4.21.
     outOfRange :: ℕ -> Bool
-    outOfRange v = v < 0 || length vertices <= fromℕ v
+    outOfRange v = v < 0 || toInteger (length vertices) <= fromℕ v
 
 -- | Find the normal of a given Triangle
 normOfTriangle :: Triangle -> ℝ3