Skip to content
Newer
Older
100644 192 lines (151 sloc) 3.96 KB
17b8b79 @rygorous RG2 license stuff
rygorous authored Apr 12, 2012
1 // This code is in the public domain. See LICENSE for details.
2
c46c5bb @rygorous RG2 Initial checkin
rygorous authored Apr 10, 2012
3 // 2d bezier curves implementation
4
5 #include "stdafx.h"
6 #include "bezier.h"
7 #include "types.h"
8 #include "math3d_2.h"
9 #include "tool.h"
10
11 namespace fr
12 {
13 // ---- cubicBezierSegment2D
14
15 void cubicBezierSegment2D::eval(vector2& out, sF32 t) const
16 {
17 sF32 tt = t * t, ttt = tt * t;
18
19 out.scale(A, 1-3*t+3*tt-ttt);
20 out.addscale(B, 3*t-6*tt+3*ttt);
21 out.addscale(C, 3*tt-3*ttt);
22 out.addscale(D, ttt);
23 }
24
25 vector2 cubicBezierSegment2D::nearestPointOnCurve(const vector2& P) const
26 {
27 sInt nSolutions;
28 vector2 w[6], p;
29 sF32 tCandidate[5];
30
31 convertToBezierForm(w, P);
32 nSolutions=findRoots(w, tCandidate);
33
34 sF32 dist=(P-A).lenSq(), t=0;
35
36 for (sInt i=0; i<nSolutions; i++)
37 {
38 eval(p, tCandidate[i]);
39 sF32 newDist=(P-p).lenSq();
40
41 if (newDist<dist)
42 {
43 dist=newDist;
44 t=tCandidate[i];
45 }
46 }
47
48 sF32 newDist=(P-D).lenSq();
49 if (newDist<dist)
50 t=1.0f;
51
52 return eval(t);
53 }
54
55 // ---- tools
56
57 void cubicBezierSegment2D::convertToBezierForm(vector2 *w, const vector2& P) const
58 {
59 vector2 c[4], d[3];
60 sF64 cdTable[3][4];
61 sInt i;
62
63 c[0].sub(A, P); c[1].sub(B, P); c[2].sub(C, P); c[3].sub(D, P);
64 d[0].sub(B, A); d[1].sub(C, B); d[2].sub(D, C);
65 d[0].scale(3); d[1].scale(3); d[2].scale(3);
66
67 for (sInt row=0; row<=2; row++)
68 {
69 for (sInt col=0; col<=3; col++)
70 cdTable[row][col]=d[row].dot(c[col]);
71 }
72
73 for (i=0; i<=5; i++)
74 {
75 w[i].y=0.0f;
76 w[i].x=i/5.0f;
77 }
78
79 static const sF32 z[3][4]={
80 { 1.0f, 0.6f, 0.3f, 0.1f},
81 { 0.4f, 0.6f, 0.6f, 0.4f},
82 { 0.1f, 0.3f, 0.6f, 1.0f}
83 };
84
85 sInt n=3, m=2, k;
86 for (k=0; k<=n+m; k++)
87 {
88 sInt lb=maximum(0, k-m);
89 sInt ub=minimum(k, n);
90
91 for (i=lb; i<=ub; i++)
92 {
93 sInt j=k-i;
94 w[i+j].y+=cdTable[j][i]*z[j][i];
95 }
96 }
97 }
98
99 sInt cubicBezierSegment2D::findRoots(const vector2 *w, sF32 *t, sInt depth) const
100 {
101 sInt cc=crossingCount(w, 5);
102
103 switch (cc)
104 {
105 case 0:
106 return 0;
107
108 case 1:
109 if (depth>=64)
110 {
111 t[0]=(w[0].x+w[5].x)*0.5;
112 return 1;
113 }
114
115 if (controlPolygonFlatEnough(w, 5))
116 {
117 t[0]=computeXIntercept(w, 5);
118 return 1;
119 }
120 }
121
122 vector2 left[6], right[6], vTemp[6][6];
123 sInt i;
124
125 for (i=0; i<6; i++)
126 vTemp[0][i]=w[i];
127
128 for (i=1; i<6; i++)
129 {
130 for (sInt j=0; j<6-i; j++)
131 vTemp[i][j]=0.5f*(vTemp[i-1][j]+vTemp[i-1][j+1]);
132 }
133
134 for (i=0; i<6; i++)
135 {
136 left[i]=vTemp[i][0];
137 right[5-i]=vTemp[i][5-i];
138 }
139
140 sF32 lt[6], rt[6];
141 sInt lCount=findRoots(left, lt, depth+1);
142 sInt rCount=findRoots(right, rt, depth+1);
143
144 for (i=0; i<lCount; i++)
145 t[i]=lt[i];
146
147 for (i=0; i<rCount; i++)
148 t[i+lCount]=rt[i];
149
150 return lCount+rCount;
151 }
152
153 sInt cubicBezierSegment2D::crossingCount(const vector2 *V, sInt degree) const
154 {
155 sInt i, nCrossings=0;
156
157 for (i=0; i<degree; i++)
158 {
159 if (V[i].y*V[i+1].y<0)
160 nCrossings++;
161 }
162
163 return nCrossings;
164 }
165
166 sBool cubicBezierSegment2D::controlPolygonFlatEnough(const vector2 *V, sInt degree) const
167 {
168 sInt i;
169 sF32 a=V[0].y-V[degree].y, b=V[degree].x-V[0].x, c=V[0].x*V[degree].y-V[degree].x*V[0].y;
170 sF32 invABSquared=1.0f / (a*a+b*b);
171 sF32 maxDistAbove=0, maxDistBelow=0;
172
173 for (i=1; i<degree; i++)
174 {
175 sF32 d=a*V[i].x+b*V[i].y+c, dv=d*d*invABSquared;
176
177 if (d<0)
178 maxDistBelow=minimum(maxDistBelow, -dv);
179 else
180 maxDistAbove=maximum(maxDistAbove, dv);
181 }
182
183 return fabs(maxDistAbove-maxDistBelow)<(-2e-8f * a);
184 }
185
186 sF32 cubicBezierSegment2D::computeXIntercept(const vector2 *V, sInt degree) const
187 {
188 sF32 XNM=V[degree].x-V[0].x, YNM=V[degree].y-V[0].y, XMK=V[0].x, YMK=V[0].y;
189 return -(XNM*YMK-YNM*XMK)/YNM;
190 }
17b8b79 @rygorous RG2 license stuff
rygorous authored Apr 13, 2012
191 }
Something went wrong with that request. Please try again.