src/bezier.cpp

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