Module for subdividing bezier curves.

1 files changed, 144 insertions, 0 deletions

diff --git a/Bezier.hs b/Bezier.hs
new file mode 100644
index 0000000..ed464d8
--- /dev/null
+++ b/Bezier.hs
@@ -0,0 +1,144 @@
+{-# LANGUAGE BangPatterns #-}
+module Bezier where
+
+import Data.List
+import Data.Ord
+
+import qualified Data.Vector.Generic as G
+import Numeric.LinearAlgebra
+
+type Plane = Vector Float
+
+ĵ :: Vector Float
+ĵ = fromList [0,1,0]
+
+type BufferID = Int
+
+data Polygonization = Polygonization
+    { curveBufferID     :: BufferID
+    , curveStartIndex   :: Int
+    , curveSegmentCount :: Int
+    }
+
+type AllocatedRange = Polygonization
+
+data Curve = Curve
+    { curvePlane    :: 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.
+    , curveD        :: Vector Float -- ^ 3d stop point.
+    , curveSegments :: Maybe Polygonization
+    }
+
+type StorePoint m = BufferID -> Int -> Vector Float {- 3d -} -> m ()
+
+xz :: Vector Float -> (Float,Float)
+xz v = (v!0, v!2)
+
+subdivideCurve :: Monad m => Float -> Curve -> AllocatedRange -> StorePoint m -> m Polygonization
+subdivideCurve δ curve range store = do
+    let ts = sort $ take (curveSegmentCount range - 2) $ sampleCurve δ curve
+    n <- (\f -> foldr f (return 0) (zip [0..] (0 : ts))) $ \(i,t) _ -> do
+        let v = bezierEval curve t
+        store (curveBufferID range) (curveStartIndex range + i) v
+        return i
+    let v = bezierEval curve 1.0
+    store (curveBufferID range) (curveStartIndex range + n + 1) v
+    return range { curveSegmentCount = n + 2 }
+
+-- | Return an unsorted list of curve parameter values to use as polygonization
+-- vertices.  The first argument is a desired maximum error distance to the
+-- true curve.
+--
+-- Note: Caller will likely want to sort and to augment the result with 0 and 1
+-- to sample the curve's end points.
+--
+-- The technique used here is from "Sampling Planar Curves Using
+-- Curvature-Based Shape Analysis" by Tatiana and Vitaly Surazhsky.  In a
+-- 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
+    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
+        (x3,y3) = xz d -- embedded in other planes.
+        -- 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 = ceiling (pi/sqrt(2*realToFrac δ/dist))
+        dist = maximum [ norm_2 (curveA curve - x) | x <- [b,c,d]]
+        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.
+    θ <- [ θ₀ + pi/fromIntegral n*fromIntegral i | 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
+ where
+    relevent x | x < 1e-12     = False
+               | 1 - x < 1e-12 = False
+               | otherwise     = True
+    a = tx*((y3 - y0) + 3*(y1 - y2)) - ty*((x3 - x0) + 3*(x1 - x2))
+    b = 2*(tx*((y2 + y0) - 2*y1) - ty*((x2 + x0) - 2*x1))
+    c = tx*(y1 - y0) - ty*(x1 - x0)
+
+-- | @quadraticRoot a b c@ find the real roots of the quadratic equation
+-- @a x^2 + b x + c = 0@.  It will return one, two or zero roots.
+quadraticRoot :: (Ord a, Floating a) => a -> a -> a -> [a]
+quadraticRoot a b c
+  | a == 0 && b == 0 = []
+  | a == 0 = [-c/b]
+  | otherwise = let d = b*b - 4*a*c
+                    q | b < 0     = - (b - sqrt d) / 2
+                      | otherwise = - (b + sqrt d) / 2
+                    x1 = q/a
+                    x2 = c/q
+                in case compare d 0 of
+                    LT -> []
+                    EQ -> [x1]
+                    GT -> [x1,x2]
+
+
+bezierEval :: Curve -> Float -> Vector Float
+bezierEval Curve{curveA=a,curveB=b,curveC=c,curveD=d} t =
+    fromList [ bernsteinEval x t, bernsteinEval y t, bernsteinEval z t]
+ where
+    poly i = fromList [a!i,b!i,c!i,d!i]
+    x = poly 0
+    y = poly 1
+    z = poly 2
+
+-- | Evaluate the bernstein polynomial using the horner rule adapted
+-- for bernstein polynomials.
+bernsteinEval :: (G.Vector v a, Fractional a) => v a -> a -> a
+bernsteinEval v _
+  | G.length v == 0 = 0
+bernsteinEval v _
+  | G.length v == 1 = G.unsafeHead v
+bernsteinEval v t =
+  go t (fromIntegral n) (G.unsafeIndex v 0 * u) 1
+  where u = 1-t
+        n = fromIntegral $ G.length v - 1
+        go !tn !bc !tmp !i
+          | i == n = tmp + tn*G.unsafeIndex v n
+          | otherwise =
+            go (tn*t) -- tn
+            (bc*fromIntegral (n-i)/(fromIntegral i + 1)) -- bc
+            ((tmp + tn*bc*G.unsafeIndex v i)*u) -- tmp
+            (i+1) -- i
+