From 175e4a0ee82e4db1fade4fbd8b5e55e88c21a826 Mon Sep 17 00:00:00 2001
From: Joe Crayne <joe@jerkface.net>
Date: Sat, 27 Jul 2019 19:31:44 -0400
Subject: Generalized curve-sampling to 3 dimensions.

---
 Bezier.hs | 67 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++----
 1 file changed, 63 insertions(+), 4 deletions(-)

diff --git a/Bezier.hs b/Bezier.hs
index 941ad5b..e526608 100644
--- a/Bezier.hs
+++ b/Bezier.hs
@@ -24,7 +24,7 @@ data Polygonization = Polygonization
 type AllocatedRange = Polygonization
 
 data Curve = Curve
-    { curvePlane    :: Plane        -- ^ Curve is embedded in this plane.
+    { curvePlane    :: Maybe Plane  -- ^ Curve is embedded in this plane.
     , curveA        :: Vector Float -- ^ 3d start point.
     , curveB        :: Vector Float -- ^ 2nd 3d control point.
     , curveC        :: Vector Float -- ^ 3rd 3d control point.
@@ -61,8 +61,24 @@ subdivideCurve δ curve range store = do
 -- nutshell: we divide the half circle into n equal parts and sample the curve
 -- whenever the tangent line is at one of these angles to a fixed axis.
 sampleCurve :: Real a => a -> Curve -> [Float]
-sampleCurve δ curve  | G.init (curvePlane curve) /= ĵ = error "TODO: arbitrary plane normal in sampleCurve"
-sampleCurve δ curve@Curve{curveA = a, curveB = b, curveC = c, curveD = d}  = do
+sampleCurve δ curve
+    | Just plane <- curvePlane curve
+    , G.init plane == ĵ  = sampleCurveXZ δ curve
+    | otherwise          = sampleCurveTwist δ curve
+
+-- | Bounding box diagonal length for a cubic curve's 4 control points.
+curveBBDiag :: Fractional b => Curve -> b
+curveBBDiag Curve{curveA = a, curveB = b, curveC = c, curveD = d} =
+    realToFrac $ norm_2 ( mx - mn )
+ where
+    mx = maximum [a,b,c,d]
+    mn = minimum [a,b,c,d]
+
+
+-- | The Surazhsky curve sampling where the Y coordinates are ignored in order
+-- that we can treat the curve in only two dimensions.
+sampleCurveXZ :: Real a => a -> Curve -> [Float]
+sampleCurveXZ δ curve@Curve{curveA = a, curveB = b, curveC = c, curveD = d}  = do
     let (x0,y0) = xz a -- We assume the curve is embedded in the xz plane so that
         (x1,y1) = xz b -- we can use a technique for 2 dimensional curves.
         (x2,y2) = xz c -- TODO: Use the correct projection here for curves
@@ -71,7 +87,7 @@ sampleCurve δ curve@Curve{curveA = a, curveB = b, curveC = c, curveD = d}  = do
         -- ought to be alright when δ is small relative to the curve's bounding
         -- box.  It was derived using the aproximation tan θ ≃ θ for small θ.
         n = ceiling (pi/sqrt(2*realToFrac δ/dist))
-        dist = maximum [ norm_2 (curveA curve - x) | x <- [b,c,d]]
+            where dist = curveBBDiag curve
         tx = x3 - x0
         ty = y3 - y0
         θ₀ = atan2 ty tx -- For our fixed axis, we use the line connecting the
@@ -88,6 +104,49 @@ sampleCurve δ curve@Curve{curveA = a, curveB = b, curveC = c, curveD = d}  = do
     θ <- [ θ₀ + pi/fromIntegral n*fromIntegral i | i <- [ 0 .. n - 1 ] ]
     tangents (cos θ) (sin θ) x0 y0 x1 y1 x2 y2 x3 y3
 
+twistAbout :: Vector Float -> Vector Float -> Vector Float -> (Float,Float)
+twistAbout ĥ k̂ v = let { x = dot ĥ v
+                       ; y = v - scale x ĥ
+                       }
+                 in (x, (if dot y k̂ <0 then negate else id) $ realToFrac $ norm_2 y)
+
+
+-- | The Surazhsky curve sampling where the middle two control points are
+-- rotated about the line connecting the curve's end points in order that they
+-- become coplanar so that we can treat the curve in only two dimensions.
+sampleCurveTwist :: Real a => a -> Curve -> [Float]
+sampleCurveTwist δ curve@Curve{curveA = a, curveB = b, curveC = c, curveD = d}  = do
+    let (x0,y0) = twistAbout ĥ k̂ a
+        (x1,y1) = twistAbout ĥ k̂ b
+        (x2,y2) = twistAbout ĥ k̂ c
+        (x3,y3) = twistAbout ĥ k̂ d
+        -- ĥ = normalize $ vector $ map realToFrac $ toList (d - a)
+        ĥ = let v = d - a in scale (1 / realToFrac (norm_2 v)) v
+        k̂ = let { k = minimumBy (comparing $ \x -> abs (dot x ĥ)) $ map (subtract a) [b,c]
+                ; v = k - scale (dot k ĥ) ĥ }
+            in scale (1 / realToFrac (norm_2 v)) v
+        -- This computation of n was not detailed in the paper, but I figure it
+        -- ought to be alright when δ is small relative to the curve's bounding
+        -- box.  It was derived using the aproximation tan θ ≃ θ for small θ.
+        n = 1 + ceiling (pi/sqrt(realToFrac δ/dist))
+            where dist = curveBBDiag curve
+        tx = x3 - x0
+        ty = y3 - y0
+        θ₀ = atan2 ty tx -- For our fixed axis, we use the line connecting the
+                         -- curve's end points.  This has the nice property
+                         -- that the first two returned sample points will be
+                         -- at the maximum distance the curve makes with its
+                         -- straight line aproximation.
+                         --
+                         -- Idea: It occurs to me that we can extend the
+                         -- Surazhsky method to 3 dimensions by obtaining a
+                         -- parameterization using polar coordinates about this
+                         -- axis.  We might ignore curvature about this axis as
+                         -- less important.
+    θ <- [ θ₀ + fromIntegral i * pi/fromIntegral n | i <- [ 0 .. n - 1 ] ]
+    tangents (cos θ) (sin θ) x0 y0 x1 y1 x2 y2 x3 y3
+
+
 tangents :: (Ord a, Floating a) =>
                   a -> a -> a -> a -> a -> a -> a -> a -> a -> a -> [a]
 tangents tx ty x0 y0 x1 y1 x2 y2 x3 y3 = filter relevent $ quadraticRoot a b c
-- 
cgit v1.2.3