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