1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
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
|
