Add quadtree back to the Makefile
[ipdf/code.git] / src / bezier.cpp
1 #include "bezier.h"
2
3 #include <unordered_map>
4 #include <cmath>
5 #include <algorithm>
6
7 using namespace std;
8
9 namespace IPDF
10 {
11
12 vector<Real> SolveQuadratic(const Real & a, const Real & b, const Real & c, const Real & min, const Real & max)
13 {
14         vector<Real> roots; roots.reserve(2);
15         if (a == 0 && b != 0)
16         {
17                 roots.push_back(-c/b);
18                 return roots;
19         }
20         Real disc(b*b - Real(4)*a*c);
21         if (disc < 0)
22         {
23                 return roots;
24         }
25         else if (disc == 0)
26         {
27                 Real x(-b/Real(2)*a);
28                 if (x >= min && x <= max)
29                         roots.push_back(x);
30                 return roots;
31         }
32         
33         Real x0((-b - Sqrt(b*b - Real(4)*a*c))/(Real(2)*a));
34         Real x1((-b + Sqrt(b*b - Real(4)*a*c))/(Real(2)*a));
35         if (x0 > x1)
36         {
37                 Real tmp(x0);
38                 x0 = x1;
39                 x1 = tmp;
40         }
41         if (x0 >= min && x0 <= max)
42                 roots.push_back(x0);
43         if (x1 >= min && x1 <= max)
44                 roots.push_back(x1);
45         return roots;
46 }
47
48 /**
49  * Finds the root (if it exists) in a monotonicly in(de)creasing segment of a Cubic
50  */
51
52 static void CubicSolveSegment(vector<Real> & roots, const Real & a, const Real & b, const Real & c, const Real & d, Real & tl, Real & tu, const Real & delta)
53 {
54         Real l = a*tl*tl*tl + b*tl*tl + c*tl + d;
55         Real u = a*tu*tu*tu + b*tu*tu + c*tu + d;
56         if ((l < 0 && u < 0) || (l > 0 && u > 0))
57                 return;
58         
59         bool negative = (u < l); // lower point > 0, upper point < 0
60         while (tu - tl > delta)
61         {
62                 Real t(tu+tl);
63                 t /= 2;
64                 Real m = a*t*t*t + b*t*t + c*t + d;
65                 if (m > 0)
66                 {
67                         if (negative)
68                                 tl = t;
69                         else
70                                 tu = t;
71                 }
72                 else if (negative)
73                 {
74                         tu = t;
75                 }
76                 else
77                 {
78                         tl = t;
79                 }
80                 //Debug("Delta is %f (%f - %f -> %f)", tu-tl, tu, tl, t);
81         }
82         roots.push_back(tl);
83 }
84 vector<Real> SolveCubic(const Real & a, const Real & b, const Real & c, const Real & d, const Real & min, const Real & max, const Real & delta)
85 {
86         vector<Real> roots; roots.reserve(3);
87         Real tu(max);
88         Real tl(min);
89         vector<Real> turns(SolveQuadratic(a*3, b*2, c));
90         //Debug("%u turning points", turns.size());
91         for (unsigned i = 1; i < turns.size(); ++i)
92         {
93                 tu = turns[i];
94                 CubicSolveSegment(roots, a, b, c, d, tl, tu,delta);
95                 tl = turns[i];
96         }
97         tu = max;
98         CubicSolveSegment(roots, a, b, c, d, tl, tu,delta);
99         return roots;
100 }
101
102 /**
103  * Factorial
104  * Use dynamic programming / recursion
105  */
106 int Factorial(int n)
107 {
108         static unordered_map<int, int> dp;
109         static bool init = false;
110         if (!init)
111         {
112                 init = true;
113                 dp[0] = 1;
114         }
115         auto it = dp.find(n);
116         if (it != dp.end())
117                 return it->second;
118         int result = n*Factorial(n-1);
119         dp[n] = result;
120         return result;
121 }
122
123 /**
124  * Binomial coefficients
125  */
126 int BinomialCoeff(int n, int k)
127 {
128         return Factorial(n) / (Factorial(k) * Factorial(n-k));
129 }
130
131 /**
132  * Bernstein Basis Polynomial
133  */
134 Real Bernstein(int k, int n, const Real & u)
135 {
136         return Real(BinomialCoeff(n, k)) * Power(u, k) * Power(Real(1.0) - u, n-k);
137 }
138
139
140 /**
141  * Returns the parametric parameter at the turning point(s)
142  * In one coordinate direction
143  */
144
145 pair<Real, Real> BezierTurningPoints(const Real & p0, const Real & p1, const Real & p2, const Real & p3)
146 {
147         // straight line
148         if (p1 == p2 && p2 == p3)
149         {
150                 return pair<Real,Real>(0, 1);
151         }
152         Real a = ((p1-p2)*3 + p3 - p0);
153         Real b = (p2 - p1*2 + p0)*2;
154         Real c = (p1-p0);
155         if (a == 0)
156         {
157                 if (b == 0)
158                         return pair<Real, Real>(0,1);
159                 Real t = -c/b;
160                 if (t > 1) t = 1;
161                 if (t < 0) t = 0;
162                 return pair<Real, Real>(t, t);
163         }
164         //Debug("a, b, c are %f, %f, %f", Float(a), Float(b), Float(c));
165         if (b*b - a*c*4 < 0)
166         {
167                 //Debug("No real roots");
168                 return pair<Real, Real>(0,1);
169         }
170         vector<Real> tsols = SolveQuadratic(a, b, c);
171         if (tsols.size() == 1)
172                 return pair<Real,Real>(tsols[0], tsols[0]);
173         else if (tsols.size() == 0)
174                 return pair<Real, Real>(0,1);
175         
176         return pair<Real,Real>(tsols[0], tsols[1]);
177         
178 }
179
180 inline bool CompRealByPtr(const Real * a, const Real * b) 
181 {
182         return (*a) < (*b);
183 }
184
185 /**
186  * Get top most *point* on Bezier curve
187  */
188 pair<Real,Real> Bezier::GetTop() const
189 {
190         pair<Real, Real> tsols = BezierTurningPoints(y0,y1,y2,y3);
191         Real tx0; Real ty0;
192         Real tx1; Real ty1;
193         Evaluate(tx0, ty0, tsols.first);
194         Evaluate(tx1, ty1, tsols.second);
195         vector<const Real*> v(4);
196         v[0] = &y0;
197         v[1] = &y3;
198         v[2] = &ty0;
199         v[3] = &ty1;
200         sort(v.begin(), v.end(), CompRealByPtr);
201         pair<Real,Real> result;
202         result.second = *v[0];
203         if (v[0] == &y0)
204         {
205                 result.first = x0;
206         }
207         else if (v[0] == &y3)
208         {
209                 result.first = x3;
210         }
211         else if (v[0] == &ty0)
212         {
213                 result.first = tx0;
214         }
215         else if (v[0] == &ty1)
216         {
217                 result.first = tx1;
218         }
219         return result;
220 }
221
222 /**
223  * Get bottom most *point* on Bezier curve
224  */
225 pair<Real,Real> Bezier::GetBottom() const
226 {
227         pair<Real, Real> tsols = BezierTurningPoints(y0,y1,y2,y3);
228         Real tx0; Real ty0;
229         Real tx1; Real ty1;
230         Evaluate(tx0, ty0, tsols.first);
231         Evaluate(tx1, ty1, tsols.second);
232         vector<const Real*> v(4);
233         v[0] = &y0;
234         v[1] = &y3;
235         v[2] = &ty0;
236         v[3] = &ty1;
237         sort(v.begin(), v.end(), CompRealByPtr);
238         pair<Real,Real> result;
239         result.second = *v[3];
240         if (v[3] == &y0)
241         {
242                 result.first = x0;
243         }
244         else if (v[3] == &y3)
245         {
246                 result.first = x3;
247         }
248         else if (v[3] == &ty0)
249         {
250                 result.first = tx0;
251         }
252         else if (v[3] == &ty1)
253         {
254                 result.first = tx1;
255         }
256         return result;
257 }
258
259 /**
260  * Get left most *point* on Bezier curve
261  */
262 pair<Real,Real> Bezier::GetLeft() const
263 {
264         pair<Real, Real> tsols = BezierTurningPoints(x0,x1,x2,x3);
265         Real tx0; Real ty0;
266         Real tx1; Real ty1;
267         Evaluate(tx0, ty0, tsols.first);
268         Evaluate(tx1, ty1, tsols.second);
269         vector<const Real*> v(4);
270         v[0] = &x0;
271         v[1] = &x3;
272         v[2] = &tx0;
273         v[3] = &tx1;
274         sort(v.begin(), v.end(), CompRealByPtr);
275         pair<Real,Real> result;
276         result.first = *v[0];
277         if (v[0] == &x0)
278         {
279                 result.second = y0;
280         }
281         else if (v[0] == &x3)
282         {
283                 result.second = y3;
284         }
285         else if (v[0] == &tx0)
286         {
287                 result.second = ty0;
288         }
289         else if (v[0] == &tx1)
290         {
291                 result.second = ty1;
292         }
293         return result;
294 }
295
296
297 /**
298  * Get left most *point* on Bezier curve
299  */
300 pair<Real,Real> Bezier::GetRight() const
301 {
302         pair<Real, Real> tsols = BezierTurningPoints(x0,x1,x2,x3);
303         Real tx0; Real ty0;
304         Real tx1; Real ty1;
305         Evaluate(tx0, ty0, tsols.first);
306         Evaluate(tx1, ty1, tsols.second);
307         vector<const Real*> v(4);
308         v[0] = &x0;
309         v[1] = &x3;
310         v[2] = &tx0;
311         v[3] = &tx1;
312         sort(v.begin(), v.end(), CompRealByPtr);
313         pair<Real,Real> result;
314         result.first = *v[3];
315         if (v[3] == &x0)
316         {
317                 result.second = y0;
318         }
319         else if (v[3] == &x3)
320         {
321                 result.second = y3;
322         }
323         else if (v[3] == &tx0)
324         {
325                 result.second = ty0;
326         }
327         else if (v[3] == &tx1)
328         {
329                 result.second = ty1;
330         }
331         return result;
332 }
333
334 vector<Real> Bezier::SolveXParam(const Real & x) const
335 {
336         Real d(x0 - x);
337         Real c((x1 - x0)*Real(3));
338         Real b((x2 - x1)*Real(3) - c);
339         Real a(x3 -x0 - c - b);
340         vector<Real> results(SolveCubic(a, b, c, d));
341         for (unsigned i = 0; i < results.size(); ++i)
342         {
343                 Vec2 p;
344                 Evaluate(p.x, p.y, results[i]);
345         }
346         return results;
347 }
348
349
350 vector<Real> Bezier::SolveYParam(const Real & y) const
351 {
352         Real d(y0 - y);
353         Real c((y1 - y0)*Real(3));
354         Real b((y2 - y1)*Real(3) - c);
355         Real a(y3 -y0 - c - b);
356         vector<Real> results(SolveCubic(a, b, c, d));
357         for (unsigned i = 0; i < results.size(); ++i)
358         {
359                 Vec2 p;
360                 Evaluate(p.x, p.y, results[i]);
361         }
362         return results;
363 }
364
365 vector<Vec2> Bezier::Evaluate(const vector<Real> & u) const
366 {
367         vector<Vec2> result(u.size());
368         for (unsigned i = 0; i < u.size(); ++i)
369         {
370                 Evaluate(result[i].x, result[i].y, u[i]);
371         }
372         return result;
373 }
374
375 /**
376  * Get Bounds Rectangle of Bezier
377  */
378 Rect Bezier::SolveBounds() const
379 {
380         Rect result;
381         pair<Real, Real> tsols = BezierTurningPoints(x0, x1, x2, x3);
382         
383         Real tp0; Real tp1; Real o;
384         Evaluate(tp0, o, tsols.first);
385         Evaluate(tp1, o, tsols.second);
386         
387         //Debug("x: tp0 is %f tp1 is %f", Float(tp0), Float(tp1));
388         
389         vector<const Real*> v(4);
390         v[0] = &x0;
391         v[1] = &x3;
392         v[2] = &tp0;
393         v[3] = &tp1;
394         
395         // Not using a lambda to keep this compiling on cabellera
396         sort(v.begin(), v.end(), CompRealByPtr);
397
398         result.x = *(v[0]);
399         result.w = *(v[3]) - result.x;
400
401         // Do the same thing for y component (wow this is a mess)
402         tsols = BezierTurningPoints(y0, y1, y2, y3);
403         Evaluate(o, tp0, tsols.first);
404         Evaluate(o, tp1, tsols.second);
405         
406         
407         //Debug("y: tp0 is %f tp1 is %f", Float(tp0), Float(tp1));
408         
409         v[0] = &y0;
410         v[1] = &y3;
411         v[2] = &tp0;
412         v[3] = &tp1;
413         sort(v.begin(), v.end(), CompRealByPtr);
414         
415         result.y = *(v[0]);
416         result.h = *(v[3]) - result.y;
417         
418         //Debug("Solved Bezier %s bounds as %s", Str().c_str(), result.Str().c_str());
419         return result;
420 }
421
422 } // end namespace

UCC git Repository :: git.ucc.asn.au