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