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 = (p1- p0 - 2*(p2-p1) + p3-p2);
  63         Real b = (p1-p0 - (p2-p1))*(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                 return pair<Real, Real>(0,1);
  78         }
  79         pair<Real, Real> tsols = SolveQuadratic(a, b, c);
  80         if (tsols.first > 1) tsols.first = 1;
  81         if (tsols.first < 0) tsols.first = 0;
  82         if (tsols.second > 1) tsols.second = 1;
  83         if (tsols.second < 0) tsols.second = 0;
  84         return tsols;
  85 }
  86
  87 inline bool CompRealByPtr(const Real * a, const Real * b)
  88 {
  89         return (*a) < (*b);
  90 }
  91
  92 /**
  93  * Get Bounds Rectangle of Bezier
  94  */
  95 Rect Bezier::SolveBounds() const
  96 {
  97         Rect result;
  98         pair<Real, Real> tsols = BezierTurningPoints(x0, x1, x2, x3);
  99
 100         Real tp0; Real tp1; Real o;
 101         Evaluate(tp0, o, tsols.first);
 102         Evaluate(tp1, o, tsols.second);
 103
 104         Debug("x: tp0 is %f tp1 is %f", Float(tp0), Float(tp1));
 105
 106         vector<const Real*> v(4);
 107         v[0] = &x0;
 108         v[1] = &x3;
 109         v[2] = &tp0;
 110         v[3] = &tp1;
 111
 112         // Not using a lambda to keep this compiling on cabellera
 113         sort(v.begin(), v.end(), CompRealByPtr);
 114
 115         result.x = *(v[0]);
 116         result.w = *(v[3]) - result.x;
 117
 118         // Do the same thing for y component (wow this is a mess)
 119         tsols = BezierTurningPoints(y0, y1, y2, y3);
 120         Evaluate(o, tp0, tsols.first);
 121         Evaluate(o, tp1, tsols.second);
 122
 123
 124         Debug("y: tp0 is %f tp1 is %f", Float(tp0), Float(tp1));
 125
 126         v[0] = &y0;
 127         v[1] = &y3;
 128         v[2] = &tp0;
 129         v[3] = &tp1;
 130         sort(v.begin(), v.end(), CompRealByPtr);
 131
 132         result.y = *(v[0]);
 133         result.h = *(v[3]) - result.y;
 134
 135         Debug("Solved Bezier %s bounds as %s", Str().c_str(), result.Str().c_str());
 136         return result;
 137 }
 138
 139 } // end namespace