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