Merge branch 'devel/master' into tizen
[platform/core/uifw/dali-adaptor.git] / 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 */