diff options
Diffstat (limited to 'Bezier.hs')
-rw-r--r-- | Bezier.hs | 144 |
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 @@ | |||
1 | {-# LANGUAGE BangPatterns #-} | ||
2 | module Bezier where | ||
3 | |||
4 | import Data.List | ||
5 | import Data.Ord | ||
6 | |||
7 | import qualified Data.Vector.Generic as G | ||
8 | import Numeric.LinearAlgebra | ||
9 | |||
10 | type Plane = Vector Float | ||
11 | |||
12 | ĵ :: Vector Float | ||
13 | ĵ = fromList [0,1,0] | ||
14 | |||
15 | type BufferID = Int | ||
16 | |||
17 | data Polygonization = Polygonization | ||
18 | { curveBufferID :: BufferID | ||
19 | , curveStartIndex :: Int | ||
20 | , curveSegmentCount :: Int | ||
21 | } | ||
22 | |||
23 | type AllocatedRange = Polygonization | ||
24 | |||
25 | data Curve = Curve | ||
26 | { curvePlane :: Plane -- ^ Curve is embedded in this plane. | ||
27 | , curveA :: Vector Float -- ^ 3d start point. | ||
28 | , curveB :: Vector Float -- ^ 2nd 3d control point. | ||
29 | , curveC :: Vector Float -- ^ 3rd 3d control point. | ||
30 | , curveD :: Vector Float -- ^ 3d stop point. | ||
31 | , curveSegments :: Maybe Polygonization | ||
32 | } | ||
33 | |||
34 | type StorePoint m = BufferID -> Int -> Vector Float {- 3d -} -> m () | ||
35 | |||
36 | xz :: Vector Float -> (Float,Float) | ||
37 | xz v = (v!0, v!2) | ||
38 | |||
39 | subdivideCurve :: Monad m => Float -> Curve -> AllocatedRange -> StorePoint m -> m Polygonization | ||
40 | subdivideCurve δ curve range store = do | ||
41 | let ts = sort $ take (curveSegmentCount range - 2) $ sampleCurve δ curve | ||
42 | n <- (\f -> foldr f (return 0) (zip [0..] (0 : ts))) $ \(i,t) _ -> do | ||
43 | let v = bezierEval curve t | ||
44 | store (curveBufferID range) (curveStartIndex range + i) v | ||
45 | return i | ||
46 | let v = bezierEval curve 1.0 | ||
47 | store (curveBufferID range) (curveStartIndex range + n + 1) v | ||
48 | return range { curveSegmentCount = n + 2 } | ||
49 | |||
50 | -- | Return an unsorted list of curve parameter values to use as polygonization | ||
51 | -- vertices. The first argument is a desired maximum error distance to the | ||
52 | -- true curve. | ||
53 | -- | ||
54 | -- Note: Caller will likely want to sort and to augment the result with 0 and 1 | ||
55 | -- to sample the curve's end points. | ||
56 | -- | ||
57 | -- The technique used here is from "Sampling Planar Curves Using | ||
58 | -- Curvature-Based Shape Analysis" by Tatiana and Vitaly Surazhsky. In a | ||
59 | -- nutshell: we divide the half circle into n equal parts and sample the curve | ||
60 | -- whenever the tangent line is at one of these angles to a fixed axis. | ||
61 | sampleCurve :: Real a => a -> Curve -> [Float] | ||
62 | sampleCurve δ curve | G.init (curvePlane curve) /= ĵ = error "TODO: arbitrary plane normal in sampleCurve" | ||
63 | sampleCurve δ curve@Curve{curveA = a, curveB = b, curveC = c, curveD = d} = do | ||
64 | let (x0,y0) = xz a -- We assume the curve is embedded in the xz plane so that | ||
65 | (x1,y1) = xz b -- we can use a technique for 2 dimensional curves. | ||
66 | (x2,y2) = xz c -- TODO: Use the correct projection here for curves | ||
67 | (x3,y3) = xz d -- embedded in other planes. | ||
68 | -- This computation of n was not detailed in the paper, but I figure it | ||
69 | -- ought to be alright when δ is small relative to the curve's bounding | ||
70 | -- box. It was derived using the aproximation tan θ ≃ θ for small θ. | ||
71 | n = ceiling (pi/sqrt(2*realToFrac δ/dist)) | ||
72 | dist = maximum [ norm_2 (curveA curve - x) | x <- [b,c,d]] | ||
73 | tx = x3 - x0 | ||
74 | ty = y3 - y0 | ||
75 | θ₀ = atan2 ty tx -- For our fixed axis, we use the line connecting the | ||
76 | -- curve's end points. This has the nice property | ||
77 | -- that the first two returned sample points will be | ||
78 | -- at the maximum distance the curve makes with its | ||
79 | -- straight line aproximation. | ||
80 | -- | ||
81 | -- Idea: It occurs to me that we can extend the | ||
82 | -- Surazhsky method to 3 dimensions by obtaining a | ||
83 | -- parameterization using polar coordinates about this | ||
84 | -- axis. We might ignore curvature about this axis as | ||
85 | -- less important. | ||
86 | θ <- [ θ₀ + pi/fromIntegral n*fromIntegral i | i <- [ 0 .. n - 1 ] ] | ||
87 | tangents (cos θ) (sin θ) x0 y0 x1 y1 x2 y2 x3 y3 | ||
88 | |||
89 | tangents :: (Ord a, Floating a) => | ||
90 | a -> a -> a -> a -> a -> a -> a -> a -> a -> a -> [a] | ||
91 | tangents tx ty x0 y0 x1 y1 x2 y2 x3 y3 = filter relevent $ quadraticRoot a b c | ||
92 | where | ||
93 | relevent x | x < 1e-12 = False | ||
94 | | 1 - x < 1e-12 = False | ||
95 | | otherwise = True | ||
96 | a = tx*((y3 - y0) + 3*(y1 - y2)) - ty*((x3 - x0) + 3*(x1 - x2)) | ||
97 | b = 2*(tx*((y2 + y0) - 2*y1) - ty*((x2 + x0) - 2*x1)) | ||
98 | c = tx*(y1 - y0) - ty*(x1 - x0) | ||
99 | |||
100 | -- | @quadraticRoot a b c@ find the real roots of the quadratic equation | ||
101 | -- @a x^2 + b x + c = 0@. It will return one, two or zero roots. | ||
102 | quadraticRoot :: (Ord a, Floating a) => a -> a -> a -> [a] | ||
103 | quadraticRoot a b c | ||
104 | | a == 0 && b == 0 = [] | ||
105 | | a == 0 = [-c/b] | ||
106 | | otherwise = let d = b*b - 4*a*c | ||
107 | q | b < 0 = - (b - sqrt d) / 2 | ||
108 | | otherwise = - (b + sqrt d) / 2 | ||
109 | x1 = q/a | ||
110 | x2 = c/q | ||
111 | in case compare d 0 of | ||
112 | LT -> [] | ||
113 | EQ -> [x1] | ||
114 | GT -> [x1,x2] | ||
115 | |||
116 | |||
117 | bezierEval :: Curve -> Float -> Vector Float | ||
118 | bezierEval Curve{curveA=a,curveB=b,curveC=c,curveD=d} t = | ||
119 | fromList [ bernsteinEval x t, bernsteinEval y t, bernsteinEval z t] | ||
120 | where | ||
121 | poly i = fromList [a!i,b!i,c!i,d!i] | ||
122 | x = poly 0 | ||
123 | y = poly 1 | ||
124 | z = poly 2 | ||
125 | |||
126 | -- | Evaluate the bernstein polynomial using the horner rule adapted | ||
127 | -- for bernstein polynomials. | ||
128 | bernsteinEval :: (G.Vector v a, Fractional a) => v a -> a -> a | ||
129 | bernsteinEval v _ | ||
130 | | G.length v == 0 = 0 | ||
131 | bernsteinEval v _ | ||
132 | | G.length v == 1 = G.unsafeHead v | ||
133 | bernsteinEval v t = | ||
134 | go t (fromIntegral n) (G.unsafeIndex v 0 * u) 1 | ||
135 | where u = 1-t | ||
136 | n = fromIntegral $ G.length v - 1 | ||
137 | go !tn !bc !tmp !i | ||
138 | | i == n = tmp + tn*G.unsafeIndex v n | ||
139 | | otherwise = | ||
140 | go (tn*t) -- tn | ||
141 | (bc*fromIntegral (n-i)/(fromIntegral i + 1)) -- bc | ||
142 | ((tmp + tn*bc*G.unsafeIndex v i)*u) -- tmp | ||
143 | (i+1) -- i | ||
144 | |||