third-party/glyphy/glyphy-geometry.hh

   1 /*
   2  * Copyright 2012,2013 Google, Inc. All Rights Reserved.
   3  *
   4  * Licensed under the Apache License, Version 2.0 (the "License");
   5  * you may not use this file except in compliance with the License.
   6  * You may obtain a copy of the License at
   7  *
   8  *     http://www.apache.org/licenses/LICENSE-2.0
   9  *
  10  * Unless required by applicable law or agreed to in writing, software
  11  * distributed under the License is distributed on an "AS IS" BASIS,
  12  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  13  * See the License for the specific language governing permissions and
  14  * limitations under the License.
  15  *
  16  * Google Author(s): Behdad Esfahbod, Maysum Panju
  17  */
  18
  19 #ifndef GLYPHY_GEOMETRY_HH
  20 #define GLYPHY_GEOMETRY_HH
  21
  22 // Glyphy is written using C style casts
  23 #pragma GCC diagnostic push
  24 #pragma GCC diagnostic ignored "-Wold-style-cast"
  25
  26 #include "glyphy-common.hh"
  27
  28 namespace GLyphy {
  29 namespace Geometry {
  30
  31 template <typename Type> struct Pair;
  32 struct Vector;
  33 struct SignedVector;
  34 struct Point;
  35 struct Line;
  36 struct Segment;
  37 struct Arc;
  38 struct Bezier;
  39
  40 /* returns tan (2 * atan (d)) */
  41 inline double tan2atan (double d) { return 2 * d / (1 - d*d); }
  42
  43 /* returns sin (2 * atan (d)) */
  44 inline double sin2atan (double d) { return 2 * d / (1 + d*d); }
  45
  46 /* returns cos (2 * atan (d)) */
  47 inline double cos2atan (double d) { return (1 - d*d) / (1 + d*d); }
  48
  49 template <typename Type>
  50 struct Pair {
  51   typedef Type ElementType;
  52
  53   inline Pair (const Type &first_, const Type &second_) : first (first_), second (second_) {}
  54
  55   Type first, second;
  56 };
  57
  58 struct Point : glyphy_point_t {
  59   inline Point (double x_, double y_) { x = x_; y = y_; }
  60   inline explicit Point (const Vector &v);
  61   inline Point (const glyphy_point_t &p) { *(glyphy_point_t *)this = p; }
  62
  63   inline bool operator == (const Point &p) const;
  64   inline bool operator != (const Point &p) const;
  65   inline Point& operator+= (const Vector &v);
  66   inline Point& operator-= (const Vector &v);
  67   inline const Point operator+ (const Vector &v) const;
  68   inline const Point operator- (const Vector &v) const;
  69   inline const Vector operator- (const Point &p) const;
  70   inline const Point midpoint (const Point &p) const;
  71   inline const Line bisector (const Point &p) const;
  72   inline double distance_to_point (const Point &p) const; /* distance to point! */
  73   inline double squared_distance_to_point (const Point &p) const; /* square of distance to point! */
  74
  75   inline bool is_finite (void) const;
  76   inline const Point lerp (const double &a, const Point &p) const;
  77 };
  78
  79 struct Vector {
  80   inline Vector (double dx_, double dy_) : dx (dx_), dy (dy_) {}
  81   inline explicit Vector (const Point &p) : dx (p.x), dy (p.y) {}
  82
  83   inline bool operator == (const Vector &v) const;
  84   inline bool operator != (const Vector &v) const;
  85   inline const Vector operator+ (void) const;
  86   inline const Vector operator- (void) const;
  87   inline Vector& operator+= (const Vector &v);
  88   inline Vector& operator-= (const Vector &v);
  89   inline Vector& operator*= (const double &s);
  90   inline Vector& operator/= (const double &s);
  91   inline const Vector operator+ (const Vector &v) const;
  92   inline const Vector operator- (const Vector &v) const;
  93   inline const Vector operator* (const double &s) const;
  94   inline const Vector operator/ (const double &s) const;
  95   inline double operator* (const Vector &v) const; /* dot product */
  96   inline const Point operator+ (const Point &p) const;
  97
  98   inline bool is_nonzero (void) const;
  99   inline double len (void) const;
 100   inline double len2 (void) const;
 101   inline const Vector normalized (void) const;
 102   inline const Vector ortho (void) const;
 103   inline const Vector normal (void) const; /* ortho().normalized() */
 104   inline double angle (void) const;
 105
 106   inline const Vector rebase (const Vector &bx, const Vector &by) const;
 107   inline const Vector rebase (const Vector &bx) const;
 108
 109   double dx, dy;
 110 };
 111
 112 struct SignedVector : Vector {
 113   inline SignedVector (const Vector &v, bool negative_) : Vector (v), negative (negative_) {}
 114
 115   inline bool operator == (const SignedVector &v) const;
 116   inline bool operator != (const SignedVector &v) const;
 117   inline const SignedVector operator- (void) const;
 118
 119   bool negative;
 120 };
 121
 122 struct Line {
 123   inline Line (double a_, double b_, double c_) : n (a_, b_), c (c_) {}
 124   inline Line (Vector n_, double c_) : n (n_), c (c_) {}
 125   inline Line (const Point &p0, const Point &p1) :
 126                n ((p1 - p0).ortho ()), c (n * Vector (p0)) {}
 127
 128   inline const Point operator+ (const Line &l) const; /* line intersection! */
 129   inline const SignedVector operator- (const Point &p) const; /* shortest vector from point to line */
 130
 131
 132   inline const Line normalized (void) const;
 133   inline const Vector normal (void) const;
 134
 135   Vector n; /* line normal */
 136   double c; /* n.dx*x + n.dy*y = c */
 137 };
 138
 139 struct Segment {
 140   inline Segment (const Point &p0_, const Point &p1_) :
 141                   p0 (p0_), p1 (p1_) {}
 142
 143   inline const SignedVector operator- (const Point &p) const; /* shortest vector from point to ***line*** */
 144   inline double distance_to_point (const Point &p) const; /* shortest distance from point to segment */
 145   inline double squared_distance_to_point (const Point &p) const; /* shortest distance squared from point to segment */
 146   inline bool contains_in_span (const Point &p) const; /* is p in the stripe formed by sliding this segment? */
 147   inline double max_distance_to_arc (const Arc &a) const;
 148
 149
 150   Point p0;
 151   Point p1;
 152 };
 153
 154
 155
 156 struct Arc {
 157   inline Arc (const Point &p0_, const Point &p1_, const Point &pm, bool complement) :
 158               p0 (p0_), p1 (p1_),
 159               d (p0_ == pm || p1_ == pm ? 0 :
 160                  tan (((p1_-pm).angle () - (p0_-pm).angle ()) / 2 - (complement ? 0 : M_PI_2))) {}
 161   inline Arc (const Point &p0_, const Point &p1_, const double &d_) :
 162               p0 (p0_), p1 (p1_), d (d_) {}
 163   inline Arc (const Point &center, double radius, const double &a0, const double &a1, bool complement) :
 164               p0 (center + Vector (cos(a0),sin(a0)) * radius),
 165               p1 (center + Vector (cos(a1),sin(a1)) * radius),
 166               d (tan ((a1 - a0) / 4 - (complement ? 0 : M_PI_2))) {}
 167   inline Arc (const glyphy_arc_t &a) : p0 (a.p0), p1 (a.p1), d (a.d) {}
 168   inline operator glyphy_arc_t (void) const { glyphy_arc_t a = {p0, p1, d}; return a; }
 169
 170   inline bool operator == (const Arc &a) const;
 171   inline bool operator != (const Arc &a) const;
 172   inline const SignedVector operator- (const Point &p) const; /* shortest vector from point to arc */
 173
 174   inline double radius (void) const;
 175   inline const Point center (void) const;
 176   inline const Pair<Vector> tangents (void) const;
 177
 178   inline Bezier approximate_bezier (double *error) const;
 179
 180   inline bool wedge_contains_point (const Point &p) const;
 181   inline double distance_to_point (const Point &p) const;
 182   inline double squared_distance_to_point (const Point &p) const;
 183   inline double extended_dist (const Point &p) const;
 184
 185   inline void extents (glyphy_extents_t &extents) const;
 186
 187   Point p0, p1;
 188   double d; /* Depth */
 189 };
 190
 191 struct Bezier {
 192   inline Bezier (const Point &p0_, const Point &p1_,
 193                  const Point &p2_, const Point &p3_) :
 194                  p0 (p0_), p1 (p1_), p2 (p2_), p3 (p3_) {}
 195
 196   inline const Point point (const double &t) const;
 197   inline const Point midpoint (void) const;
 198   inline const Vector tangent (const double &t) const;
 199   inline const Vector d_tangent (const double &t) const;
 200   inline double curvature (const double &t) const;
 201   inline const Pair<Bezier> split (const double &t) const;
 202   inline const Pair<Bezier> halve (void) const;
 203   inline const Bezier segment (const double &t0, const double &t1) const;
 204
 205   Point p0, p1, p2, p3;
 206 };
 207
 208
 209 /* Implementations */
 210
 211
 212 /* Point */
 213
 214 inline Point::Point (const Vector &v) {
 215   x = v.dx;
 216   y = v.dy;
 217 }
 218 inline bool Point::operator == (const Point &p) const {
 219   return x == p.x && y == p.y;
 220 }
 221 inline bool Point::operator != (const Point &p) const {
 222   return !(*this == p);
 223 }
 224 inline Point& Point::operator+= (const Vector &v) {
 225   x += v.dx;
 226   y += v.dy;
 227   return *this;
 228 }
 229 inline Point& Point::operator-= (const Vector &v) {
 230   x -= v.dx;
 231   y -= v.dy;
 232   return *this;
 233 }
 234 inline const Point Point::operator+ (const Vector &v) const {
 235   return Point (*this) += v;
 236 }
 237 inline const Point Point::operator- (const Vector &v) const {
 238   return Point (*this) -= v;
 239 }
 240 inline const Vector Point::operator- (const Point &p) const {
 241   return Vector (x - p.x, y - p.y);
 242 }
 243
 244 inline const Point Point::midpoint (const Point &p) const {
 245   return *this + (p - *this) / 2;
 246 }
 247 inline const Line Point::bisector (const Point &p) const {
 248   Vector d = p - *this;
 249   return Line (d.dx * 2, d.dy * 2, d * Vector (p) + d * Vector (*this));
 250 }
 251
 252 inline double Point::distance_to_point (const Point &p) const {
 253   return ((*this) - p).len ();
 254 }
 255
 256 inline double Point::squared_distance_to_point (const Point &p) const {
 257   return ((*this) - p).len2 ();
 258 }
 259
 260 inline bool Point::is_finite (void) const {
 261   return std::isfinite (x) && std::isfinite (y);
 262 }
 263 inline const Point Point::lerp (const double &a, const Point &p) const {
 264   /* The following two cases are special-cased to get better floating
 265    * point stability.  We require that points that are the same be
 266    * bit-equal. */
 267   if (a == 0)   return *this;
 268   if (a == 1.0) return p;
 269   return Point ((1-a) * x + a * p.x, (1-a) * y + a * p.y);
 270 }
 271
 272
 273 /* Vector */
 274
 275 inline bool Vector::operator == (const Vector &v) const {
 276   return dx == v.dx && dy == v.dy;
 277 }
 278 inline bool Vector::operator != (const Vector &v) const {
 279   return !(*this == v);
 280 }
 281 inline const Vector Vector::operator+ (void) const {
 282   return *this;
 283 }
 284 inline const Vector Vector::operator- (void) const {
 285   return Vector (-dx, -dy);
 286 }
 287 inline Vector& Vector::operator+= (const Vector &v) {
 288   dx += v.dx;
 289   dy += v.dy;
 290   return *this;
 291 }
 292 inline Vector& Vector::operator-= (const Vector &v) {
 293   dx -= v.dx;
 294   dy -= v.dy;
 295   return *this;
 296 }
 297 inline Vector& Vector::operator*= (const double &s) {
 298   dx *= s;
 299   dy *= s;
 300   return *this;
 301 }
 302 inline Vector& Vector::operator/= (const double &s) {
 303   dx /= s;
 304   dy /= s;
 305   return *this;
 306 }
 307 inline const Vector Vector::operator+ (const Vector &v) const {
 308   return Vector (*this) += v;
 309 }
 310 inline const Vector Vector::operator- (const Vector &v) const {
 311   return Vector (*this) -= v;
 312 }
 313 inline const Vector Vector::operator* (const double &s) const {
 314   return Vector (*this) *= s;
 315 }
 316 inline const Vector operator* (const double &s, const Vector &v) {
 317   return v * s;
 318 }
 319 inline const Vector Vector::operator/ (const double &s) const {
 320   return Vector (*this) /= s;
 321 }
 322 inline double Vector::operator* (const Vector &v) const { /* dot product */
 323   return dx * v.dx + dy * v.dy;
 324 }
 325 inline const Point Vector::operator+ (const Point &p) const {
 326   return p + *this;
 327 }
 328
 329 inline bool Vector::is_nonzero (void) const {
 330   return dx || dy;
 331 }
 332 inline double Vector::len (void) const {
 333   return hypot (dx, dy);
 334 }
 335 inline double Vector::len2 (void) const {
 336   return dx * dx + dy * dy;
 337 }
 338 inline const Vector Vector::normalized (void) const {
 339   double d = len ();
 340   return d ? *this / d : *this;
 341 }
 342 inline const Vector Vector::ortho (void) const {
 343   return Vector (-dy, dx);
 344 }
 345 inline const Vector Vector::normal (void) const {
 346   return ortho ().normalized ();
 347 }
 348 inline double Vector::angle (void) const {
 349   return atan2 (dy, dx);
 350 }
 351
 352 inline const Vector Vector::rebase (const Vector &bx,
 353                                     const Vector &by) const {
 354   return Vector (*this * bx, *this * by);
 355 }
 356 inline const Vector Vector::rebase (const Vector &bx) const {
 357   return rebase (bx, bx.ortho ());
 358 }
 359
 360
 361 /* SignedVector */
 362
 363 inline bool SignedVector::operator == (const SignedVector &v) const {
 364   return (const Vector &)(*this) == (const Vector &)(v) && negative == v.negative;
 365 }
 366 inline bool SignedVector::operator != (const SignedVector &v) const {
 367   return !(*this == v);
 368 }
 369 inline const SignedVector SignedVector::operator- (void) const {
 370   return SignedVector (-(const Vector &)(*this), !negative);
 371 }
 372
 373
 374 /* Line */
 375
 376 inline const Point Line::operator+ (const Line &l) const {
 377   double det = n.dx * l.n.dy - n.dy * l.n.dx;
 378   if (!det)
 379     return Point (GLYPHY_INFINITY, GLYPHY_INFINITY);
 380   return Point ((c * l.n.dy - n.dy * l.c) / det,
 381                        (n.dx * l.c - c * l.n.dx) / det);
 382 }
 383 inline const SignedVector Line::operator- (const Point &p) const {
 384   double mag = -(n * Vector (p) - c) / n.len ();
 385   return SignedVector (n.normalized () * mag, mag < 0); /******************************************************************************************* FIX. *************************************/
 386 }
 387
 388 inline const SignedVector operator- (const Point &p, const Line &l) {
 389   return -(l - p);
 390 }
 391
 392 inline const Line Line::normalized (void) const {
 393   double d = n.len ();
 394   return d ? Line (n / d, c / d) : *this;
 395 }
 396 inline const Vector Line::normal (void) const {
 397   return n;
 398 }
 399
 400 /* Segment */
 401 inline const SignedVector Segment::operator- (const Point &p) const {
 402   /* shortest vector from point to line */
 403   return p - Line (p1, p0); /************************************************************************************************** Should the order (p1, p0) depend on d?? ***********************/
 404 }
 405
 406 /* Segment */
 407 inline bool Segment::contains_in_span (const Point &p) const {
 408   if (p0 == p1)
 409     return false;
 410
 411   /* shortest vector from point to line */
 412   Line temp (p0, p1);
 413   double mag = -(temp.n * Vector (p) - temp.c) / temp.n.len ();
 414   Vector y (temp.n.normalized () * mag);
 415   Point z = y + p;
 416
 417   // Check if z is between p0 and p1.
 418
 419   if (fabs (p1.y - p0.y) > fabs (p1.x - p0.x)) {
 420     return ((z.y - p0.y > 0 && p1.y - p0.y > z.y - p0.y) ||
 421             (z.y - p0.y < 0 && p1.y - p0.y < z.y - p0.y));
 422   }
 423   else {
 424     return ((0 < z.x - p0.x && z.x - p0.x < p1.x - p0.x) ||
 425             (0 > z.x - p0.x && z.x - p0.x > p1.x - p0.x));
 426   }
 427 }
 428
 429 inline double Segment::distance_to_point (const Point &p) const {
 430   if (p0 == p1)
 431     return 0;
 432
 433   // Check if z is between p0 and p1.
 434   Line temp (p0, p1);
 435   if (contains_in_span (p))
 436     return -(temp.n * Vector (p) - temp.c) / temp.n.len ();
 437
 438   double dist_p_p0 = p.distance_to_point (p0);
 439   double dist_p_p1 = p.distance_to_point (p1);
 440   return (dist_p_p0 < dist_p_p1 ? dist_p_p0 : dist_p_p1) * (-(temp.n * Vector (p) - temp.c) < 0 ? -1 : 1);
 441 }
 442
 443
 444 inline double Segment::squared_distance_to_point (const Point &p) const {
 445   if (p0 == p1)
 446     return 0;
 447
 448   // Check if z is between p0 and p1.
 449   Line temp (p0, p1);
 450   if (contains_in_span (p))
 451     return (temp.n * Vector (p) - temp.c) * (temp.n * Vector (p) - temp.c) / (temp.n * temp.n);
 452
 453   double dist_p_p0 = p.squared_distance_to_point (p0);
 454   double dist_p_p1 = p.squared_distance_to_point (p1);
 455   return (dist_p_p0 < dist_p_p1 ? dist_p_p0 : dist_p_p1);
 456 }
 457
 458
 459 inline double Segment::max_distance_to_arc (const Arc &a) const {
 460   double max_distance = fabs(a.distance_to_point(p0)) ;
 461   return  max_distance >  fabs(a.distance_to_point(p1)) ? max_distance : fabs(a.distance_to_point(p1)) ;
 462 }
 463
 464
 465
 466 /* Arc */
 467
 468 inline bool Arc::operator == (const Arc &a) const {
 469   return p0 == a.p0 && p1 == a.p1 && d == a.d;
 470 }
 471 inline bool Arc::operator != (const Arc &a) const {
 472   return !(*this == a);
 473 }
 474
 475
 476 inline const SignedVector Arc::operator- (const Point &p) const {
 477
 478   if (fabs(d) < 1e-5) {
 479     Segment arc_segment (p0, p1);
 480     return arc_segment - p;
 481   }
 482   if (wedge_contains_point (p)){
 483     Vector difference = (center () - p).normalized () * fabs (p.distance_to_point (center ()) - radius ());
 484
 485     return SignedVector  (difference, ((p - center ()).len () < radius ()) ^ (d < 0));
 486   }
 487   double d0 = p.squared_distance_to_point (p0);
 488   double d1 = p.squared_distance_to_point (p1);
 489
 490   Arc other_arc (p0, p1, (1.0 + d) / (1.0 - d));  /********************************* NOT Robust. But works? *****************/
 491   Vector normal = center () - (d0 < d1 ? p0 : p1) ;
 492
 493   if (normal.len() == 0)
 494     return SignedVector (Vector (0, 0), true);    /************************************ Check sign of this S.D. *************/
 495
 496   return SignedVector (Line (normal.dx, normal.dy, normal * Vector ((d0 < d1 ? p0 : p1))) - p, !other_arc.wedge_contains_point(p));
 497 }
 498
 499 inline const SignedVector operator- (const Point &p, const Arc &a) {
 500   return -(a - p);
 501 }
 502
 503
 504
 505 inline double Arc::radius (void) const
 506 {
 507   return fabs ((p1 - p0).len () / (2 * sin2atan (d)));
 508 }
 509
 510 inline const Point Arc::center (void) const
 511 {
 512   return (p0.midpoint (p1)) + (p1 - p0).ortho () / (2 * tan2atan (d));
 513 }
 514
 515 inline const Pair<Vector> Arc::tangents (void) const
 516 {
 517   Vector dp = (p1 - p0) * .5;
 518   Vector pp = dp.ortho () * -sin2atan (d);
 519   dp = dp * cos2atan (d);
 520   return Pair<Vector> (dp + pp, dp - pp);
 521 }
 522
 523
 524
 525 inline Bezier Arc::approximate_bezier (double *error) const
 526 {
 527   Vector dp = p1 - p0;
 528   Vector pp = dp.ortho ();
 529
 530   if (error)
 531     *error = dp.len () * pow (fabs (d), 5) / (54 * (1 + d*d));
 532
 533   dp *= ((1 - d*d) / 3);
 534   pp *= (2 * d / 3);
 535
 536   Point p0s = p0 + dp - pp;
 537   Point p1s = p1 - dp - pp;
 538
 539   return Bezier (p0, p0s, p1s, p1);
 540 }
 541
 542
 543 inline bool Arc::wedge_contains_point (const Point &p) const
 544 {
 545   Pair<Vector> t = tangents ();
 546   if (fabs (d) <= 1)
 547     return (p - p0) * t.first  >= 0 && (p - p1) * t.second <= 0;
 548   else
 549     return (p - p0) * t.first  >= 0 || (p - p1) * t.second <= 0;
 550 }
 551
 552
 553 /* Distance may not always be positive, but will be to an endpoint whenever necessary. */
 554 inline double Arc::distance_to_point (const Point &p) const {
 555   if (fabs(d) < 1e-5) {
 556     Segment arc_segment (p0, p1);
 557     return arc_segment.distance_to_point (p);
 558   }
 559
 560   SignedVector difference = *this - p;
 561
 562   if (wedge_contains_point (p) && fabs(d) > 1e-5)
 563     return fabs (p.distance_to_point (center ()) - radius ()) * (difference.negative ? -1 : 1);
 564   double d1 = p.squared_distance_to_point (p0);
 565   double d2 = p.squared_distance_to_point (p1);
 566   return (d1 < d2 ? sqrt(d1) : sqrt(d2)) * (difference.negative ? -1 : 1);
 567 }
 568
 569 /* Distance will be to an endpoint whenever necessary. */
 570 inline double Arc::squared_distance_to_point (const Point &p) const {
 571   if (fabs(d) < 1e-5) {
 572     Segment arc_segment (p0, p1);
 573     return arc_segment.squared_distance_to_point (p);
 574   }
 575
 576   //SignedVector difference = *this - p;
 577
 578   if (wedge_contains_point (p) && fabs(d) > 1e-5) {
 579     double answer = p.distance_to_point (center ()) - radius ();
 580     return answer * answer;
 581   }
 582   double d1 = p.squared_distance_to_point (p0);
 583   double d2 = p.squared_distance_to_point (p1);
 584   return (d1 < d2 ? d1 : d2);
 585 }
 586
 587 inline double Arc::extended_dist (const Point &p) const {
 588   Point m = p0.lerp (.5, p1);
 589   Vector dp = p1 - p0;
 590   Vector pp = dp.ortho ();
 591   /*
 592    * In the original file, there was no explicit cast to float.
 593    */
 594   float d2 = static_cast<float>( tan2atan (d) );
 595   if ((p - m) * (p1 - m) < 0)
 596     return (p - p0) * (pp + dp * d2).normalized ();
 597   else
 598     return (p - p1) * (pp - dp * d2).normalized ();
 599 }
 600
 601 inline void Arc::extents (glyphy_extents_t &extents) const {
 602   glyphy_extents_clear (&extents);
 603   glyphy_extents_add (&extents, &p0);
 604   glyphy_extents_add (&extents, &p1);
 605   Point c = center ();
 606   double r = radius ();
 607   Point p[4] = {c + r * Vector (-1,  0),
 608                 c + r * Vector (+1,  0),
 609                 c + r * Vector ( 0, -1),
 610                 c + r * Vector ( 0, +1)};
 611   for (unsigned int i = 0; i < 4; i++)
 612     if (wedge_contains_point (p[i]))
 613       glyphy_extents_add (&extents, &p[i]);
 614 }
 615
 616
 617 /* Bezier */
 618
 619 inline const Point Bezier::point (const double &t) const {
 620   Point p01 = p0.lerp (t, p1);
 621   Point p12 = p1.lerp (t, p2);
 622   Point p23 = p2.lerp (t, p3);
 623   Point p012 = p01.lerp (t, p12);
 624   Point p123 = p12.lerp (t, p23);
 625   Point p0123 = p012.lerp (t, p123);
 626   return p0123;
 627 }
 628
 629 inline const Point Bezier::midpoint (void) const
 630 {
 631   Point p01 = p0.midpoint (p1);
 632   Point p12 = p1.midpoint (p2);
 633   Point p23 = p2.midpoint (p3);
 634   Point p012 = p01.midpoint (p12);
 635   Point p123 = p12.midpoint (p23);
 636   Point p0123 = p012.midpoint (p123);
 637   return p0123;
 638 }
 639
 640 inline const Vector Bezier::tangent (const double &t) const
 641 {
 642   double t_2_0 = t * t;
 643   double t_0_2 = (1 - t) * (1 - t);
 644
 645   double _1__4t_1_0_3t_2_0 = 1 - 4 * t + 3 * t_2_0;
 646   double _2t_1_0_3t_2_0    =     2 * t - 3 * t_2_0;
 647
 648   return Vector (-3 * p0.x * t_0_2
 649                         +3 * p1.x * _1__4t_1_0_3t_2_0
 650                         +3 * p2.x * _2t_1_0_3t_2_0
 651                         +3 * p3.x * t_2_0,
 652                         -3 * p0.y * t_0_2
 653                         +3 * p1.y * _1__4t_1_0_3t_2_0
 654                         +3 * p2.y * _2t_1_0_3t_2_0
 655                         +3 * p3.y * t_2_0);
 656 }
 657
 658 inline const Vector Bezier::d_tangent (const double &t) const {
 659   return Vector (6 * ((-p0.x + 3*p1.x - 3*p2.x + p3.x) * t + (p0.x - 2*p1.x + p2.x)),
 660                         6 * ((-p0.y + 3*p1.y - 3*p2.y + p3.y) * t + (p0.y - 2*p1.y + p2.y)));
 661 }
 662
 663 inline double Bezier::curvature (const double &t) const {
 664   Vector dpp = tangent (t).ortho ();
 665   Vector ddp = d_tangent (t);
 666   /* normal vector len squared */
 667   double len = dpp.len ();
 668   double curvature = (dpp * ddp) / (len * len * len);
 669   return curvature;
 670 }
 671
 672 inline const Pair<Bezier > Bezier::split (const double &t) const {
 673   Point p01 = p0.lerp (t, p1);
 674   Point p12 = p1.lerp (t, p2);
 675   Point p23 = p2.lerp (t, p3);
 676   Point p012 = p01.lerp (t, p12);
 677   Point p123 = p12.lerp (t, p23);
 678   Point p0123 = p012.lerp (t, p123);
 679   return Pair<Bezier> (Bezier (p0, p01, p012, p0123),
 680                        Bezier (p0123, p123, p23, p3));
 681 }
 682
 683 inline const Pair<Bezier > Bezier::halve (void) const
 684 {
 685   Point p01 = p0.midpoint (p1);
 686   Point p12 = p1.midpoint (p2);
 687   Point p23 = p2.midpoint (p3);
 688   Point p012 = p01.midpoint (p12);
 689   Point p123 = p12.midpoint (p23);
 690   Point p0123 = p012.midpoint (p123);
 691   return Pair<Bezier> (Bezier (p0, p01, p012, p0123),
 692                        Bezier (p0123, p123, p23, p3));
 693 }
 694
 695 inline const Bezier Bezier::segment (const double &t0, const double &t1) const
 696 {
 697   Point p01 = p0.lerp (t0, p1);
 698   Point p12 = p1.lerp (t0, p2);
 699   Point p23 = p2.lerp (t0, p3);
 700   Point p012 = p01.lerp (t0, p12);
 701   Point p123 = p12.lerp (t0, p23);
 702   Point p0123 = p012.lerp (t0, p123);
 703
 704   Point q01 = p0.lerp (t1, p1);
 705   Point q12 = p1.lerp (t1, p2);
 706   Point q23 = p2.lerp (t1, p3);
 707   Point q012 = q01.lerp (t1, q12);
 708   Point q123 = q12.lerp (t1, q23);
 709   Point q0123 = q012.lerp (t1, q123);
 710
 711   return Bezier (p0123,
 712                  p0123 + (p123 - p0123) * ((t1 - t0) / (1 - t0)),
 713                  q0123 + (q012 - q0123) * ((t1 - t0) / t1),
 714                  q0123);
 715 }
 716
 717 } /* namespace Geometry */
 718 } /* namespace GLyphy */
 719
 720 #pragma GCC diagnostic pop
 721
 722 #endif /* GLYPHY_GEOMETRY_HH */