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 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 top most *point* on Bezier curve
  95  */
  96 pair<Real,Real> Bezier::GetTop() const
  97 {
  98         pair<Real, Real> tsols = BezierTurningPoints(y0,y1,y2,y3);
  99         Real tx0; Real ty0;
 100         Real tx1; Real ty1;
 101         Evaluate(tx0, ty0, tsols.first);
 102         Evaluate(tx1, ty1, tsols.second);
 103         vector<const Real*> v(4);
 104         v[0] = &y0;
 105         v[1] = &y3;
 106         v[2] = &ty0;
 107         v[3] = &ty1;
 108         sort(v.begin(), v.end(), CompRealByPtr);
 109         pair<Real,Real> result;
 110         result.second = *v[0];
 111         if (v[0] == &y0)
 112         {
 113                 result.first = x0;
 114         }
 115         else if (v[0] == &y3)
 116         {
 117                 result.first = x3;
 118         }
 119         else if (v[0] == &ty0)
 120         {
 121                 result.first = tx0;
 122         }
 123         else if (v[0] == &ty1)
 124         {
 125                 result.first = tx1;
 126         }
 127         return result;
 128 }
 129
 130 /**
 131  * Get bottom most *point* on Bezier curve
 132  */
 133 pair<Real,Real> Bezier::GetBottom() const
 134 {
 135         pair<Real, Real> tsols = BezierTurningPoints(y0,y1,y2,y3);
 136         Real tx0; Real ty0;
 137         Real tx1; Real ty1;
 138         Evaluate(tx0, ty0, tsols.first);
 139         Evaluate(tx1, ty1, tsols.second);
 140         vector<const Real*> v(4);
 141         v[0] = &y0;
 142         v[1] = &y3;
 143         v[2] = &ty0;
 144         v[3] = &ty1;
 145         sort(v.begin(), v.end(), CompRealByPtr);
 146         pair<Real,Real> result;
 147         result.second = *v[3];
 148         if (v[3] == &y0)
 149         {
 150                 result.first = x0;
 151         }
 152         else if (v[3] == &y3)
 153         {
 154                 result.first = x3;
 155         }
 156         else if (v[3] == &ty0)
 157         {
 158                 result.first = tx0;
 159         }
 160         else if (v[3] == &ty1)
 161         {
 162                 result.first = tx1;
 163         }
 164         return result;
 165 }
 166
 167 /**
 168  * Get left most *point* on Bezier curve
 169  */
 170 pair<Real,Real> Bezier::GetLeft() const
 171 {
 172         pair<Real, Real> tsols = BezierTurningPoints(x0,x1,x2,x3);
 173         Real tx0; Real ty0;
 174         Real tx1; Real ty1;
 175         Evaluate(tx0, ty0, tsols.first);
 176         Evaluate(tx1, ty1, tsols.second);
 177         vector<const Real*> v(4);
 178         v[0] = &x0;
 179         v[1] = &x3;
 180         v[2] = &tx0;
 181         v[3] = &tx1;
 182         sort(v.begin(), v.end(), CompRealByPtr);
 183         pair<Real,Real> result;
 184         result.first = *v[0];
 185         if (v[0] == &x0)
 186         {
 187                 result.second = y0;
 188         }
 189         else if (v[0] == &x3)
 190         {
 191                 result.second = y3;
 192         }
 193         else if (v[0] == &tx0)
 194         {
 195                 result.second = ty0;
 196         }
 197         else if (v[0] == &tx1)
 198         {
 199                 result.second = ty1;
 200         }
 201         return result;
 202 }
 203
 204
 205 /**
 206  * Get left most *point* on Bezier curve
 207  */
 208 pair<Real,Real> Bezier::GetRight() const
 209 {
 210         pair<Real, Real> tsols = BezierTurningPoints(x0,x1,x2,x3);
 211         Real tx0; Real ty0;
 212         Real tx1; Real ty1;
 213         Evaluate(tx0, ty0, tsols.first);
 214         Evaluate(tx1, ty1, tsols.second);
 215         vector<const Real*> v(4);
 216         v[0] = &x0;
 217         v[1] = &x3;
 218         v[2] = &tx0;
 219         v[3] = &tx1;
 220         sort(v.begin(), v.end(), CompRealByPtr);
 221         pair<Real,Real> result;
 222         result.first = *v[3];
 223         if (v[3] == &x0)
 224         {
 225                 result.second = y0;
 226         }
 227         else if (v[3] == &x3)
 228         {
 229                 result.second = y3;
 230         }
 231         else if (v[3] == &tx0)
 232         {
 233                 result.second = ty0;
 234         }
 235         else if (v[3] == &tx1)
 236         {
 237                 result.second = ty1;
 238         }
 239         return result;
 240 }
 241
 242 /**
 243  * Get Bounds Rectangle of Bezier
 244  */
 245 Rect Bezier::SolveBounds() const
 246 {
 247         Rect result;
 248         pair<Real, Real> tsols = BezierTurningPoints(x0, x1, x2, x3);
 249
 250         Real tp0; Real tp1; Real o;
 251         Evaluate(tp0, o, tsols.first);
 252         Evaluate(tp1, o, tsols.second);
 253
 254         //Debug("x: tp0 is %f tp1 is %f", Float(tp0), Float(tp1));
 255
 256         vector<const Real*> v(4);
 257         v[0] = &x0;
 258         v[1] = &x3;
 259         v[2] = &tp0;
 260         v[3] = &tp1;
 261
 262         // Not using a lambda to keep this compiling on cabellera
 263         sort(v.begin(), v.end(), CompRealByPtr);
 264
 265         result.x = *(v[0]);
 266         result.w = *(v[3]) - result.x;
 267
 268         // Do the same thing for y component (wow this is a mess)
 269         tsols = BezierTurningPoints(y0, y1, y2, y3);
 270         Evaluate(o, tp0, tsols.first);
 271         Evaluate(o, tp1, tsols.second);
 272
 273
 274         //Debug("y: tp0 is %f tp1 is %f", Float(tp0), Float(tp1));
 275
 276         v[0] = &y0;
 277         v[1] = &y3;
 278         v[2] = &tp0;
 279         v[3] = &tp1;
 280         sort(v.begin(), v.end(), CompRealByPtr);
 281
 282         result.y = *(v[0]);
 283         result.h = *(v[3]) - result.y;
 284
 285         //Debug("Solved Bezier %s bounds as %s", Str().c_str(), result.Str().c_str());
 286         return result;
 287 }
 288
 289 } // end namespace