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 /**
  13  * Factorial
  14  * Use dynamic programming / recursion
  15  */
  16 int Factorial(int n)
  17 {
  18         static unordered_map<int, int> dp;
  19         static bool init = false;
  20         if (!init)
  21         {
  22                 init = true;
  23                 dp[0] = 1;
  24         }
  25         auto it = dp.find(n);
  26         if (it != dp.end())
  27                 return it->second;
  28         int result = n*Factorial(n-1);
  29         dp[n] = result;
  30         return result;
  31 }
  32
  33 /**
  34  * Binomial coefficients
  35  */
  36 int BinomialCoeff(int n, int k)
  37 {
  38         return Factorial(n) / (Factorial(k) * Factorial(n-k));
  39 }
  40
  41 /**
  42  * Bernstein Basis Polynomial
  43  */
  44 Real Bernstein(int k, int n, const Real & u)
  45 {
  46         return Real(BinomialCoeff(n, k)) * Power(u, k) * Power(Real(1.0) - u, n-k);
  47 }
  48
  49
  50 /**
  51  * Returns the parametric parameter at the turning point(s)
  52  * In one coordinate direction
  53  */
  54
  55 static pair<Real, Real> BezierTurningPoints(const Real & p0, const Real & p1, const Real & p2, const Real & p3)
  56 {
  57         // straight line
  58         if (p1 == p2 && p2 == p3)
  59         {
  60                 return pair<Real,Real>(0, 1);
  61         }
  62         Real a = (3*(p1-p2) + p3 - p0);
  63         Real b = 2*(p2 - 2*p1 + p0);
  64         Real c = (p1-p0);
  65         if (a == 0)
  66         {
  67                 if (b == 0)
  68                         return pair<Real, Real>(0,1);
  69                 Real t = -c/b;
  70                 if (t > 1) t = 1;
  71                 if (t < 0) t = 0;
  72                 return pair<Real, Real>(t, t);
  73         }
  74         //Debug("a, b, c are %f, %f, %f", Float(a), Float(b), Float(c));
  75         if (b*b - 4*a*c < 0)
  76         {
  77                 //Debug("No real roots");
  78                 return pair<Real, Real>(0,1);
  79         }
  80         pair<Real, Real> tsols = SolveQuadratic(a, b, c);
  81         if (tsols.first > 1) tsols.first = 1;
  82         if (tsols.first < 0) tsols.first = 0;
  83         if (tsols.second > 1) tsols.second = 1;
  84         if (tsols.second < 0) tsols.second = 0;
  85         return tsols;
  86 }
  87
  88 inline bool CompRealByPtr(const Real * a, const Real * b)
  89 {
  90         return (*a) < (*b);
  91 }
  92
  93 /**
  94  * Get Bounds Rectangle of Bezier
  95  */
  96 Rect Bezier::SolveBounds() const
  97 {
  98         Rect result;
  99         pair<Real, Real> tsols = BezierTurningPoints(x0, x1, x2, x3);
 100
 101         Real tp0; Real tp1; Real o;
 102         Evaluate(tp0, o, tsols.first);
 103         Evaluate(tp1, o, tsols.second);
 104
 105         //Debug("x: tp0 is %f tp1 is %f", Float(tp0), Float(tp1));
 106
 107         vector<const Real*> v(4);
 108         v[0] = &x0;
 109         v[1] = &x3;
 110         v[2] = &tp0;
 111         v[3] = &tp1;
 112
 113         // Not using a lambda to keep this compiling on cabellera
 114         sort(v.begin(), v.end(), CompRealByPtr);
 115
 116         result.x = *(v[0]);
 117         result.w = *(v[3]) - result.x;
 118
 119         // Do the same thing for y component (wow this is a mess)
 120         tsols = BezierTurningPoints(y0, y1, y2, y3);
 121         Evaluate(o, tp0, tsols.first);
 122         Evaluate(o, tp1, tsols.second);
 123
 124
 125         //Debug("y: tp0 is %f tp1 is %f", Float(tp0), Float(tp1));
 126
 127         v[0] = &y0;
 128         v[1] = &y3;
 129         v[2] = &tp0;
 130         v[3] = &tp1;
 131         sort(v.begin(), v.end(), CompRealByPtr);
 132
 133         result.y = *(v[0]);
 134         result.h = *(v[3]) - result.y;
 135
 136         //Debug("Solved Bezier %s bounds as %s", Str().c_str(), result.Str().c_str());
 137         return result;
 138 }
 139
 140 } // end namespace