2 * Copyright 2012,2013 Google, Inc. All Rights Reserved.
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
8 * http://www.apache.org/licenses/LICENSE-2.0
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.
16 * Google Author(s): Behdad Esfahbod, Maysum Panju
19 #ifndef GLYPHY_GEOMETRY_HH
20 #define GLYPHY_GEOMETRY_HH
22 // Glyphy is written using C style casts
23 #pragma GCC diagnostic push
24 #pragma GCC diagnostic ignored "-Wold-style-cast"
26 #include "glyphy-common.hh"
31 template <typename Type> struct Pair;
40 /* returns tan (2 * atan (d)) */
41 inline double tan2atan (double d) { return 2 * d / (1 - d*d); }
43 /* returns sin (2 * atan (d)) */
44 inline double sin2atan (double d) { return 2 * d / (1 + d*d); }
46 /* returns cos (2 * atan (d)) */
47 inline double cos2atan (double d) { return (1 - d*d) / (1 + d*d); }
49 template <typename Type>
51 typedef Type ElementType;
53 inline Pair (const Type &first_, const Type &second_) : first (first_), second (second_) {}
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; }
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! */
75 inline bool is_finite (void) const;
76 inline const Point lerp (const double &a, const Point &p) const;
80 inline Vector (double dx_, double dy_) : dx (dx_), dy (dy_) {}
81 inline explicit Vector (const Point &p) : dx (p.x), dy (p.y) {}
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;
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;
106 inline const Vector rebase (const Vector &bx, const Vector &by) const;
107 inline const Vector rebase (const Vector &bx) const;
112 struct SignedVector : Vector {
113 inline SignedVector (const Vector &v, bool negative_) : Vector (v), negative (negative_) {}
115 inline bool operator == (const SignedVector &v) const;
116 inline bool operator != (const SignedVector &v) const;
117 inline const SignedVector operator- (void) const;
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)) {}
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 */
132 inline const Line normalized (void) const;
133 inline const Vector normal (void) const;
135 Vector n; /* line normal */
136 double c; /* n.dx*x + n.dy*y = c */
140 inline Segment (const Point &p0_, const Point &p1_) :
141 p0 (p0_), p1 (p1_) {}
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;
157 inline Arc (const Point &p0_, const Point &p1_, const Point &pm, bool complement) :
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 ¢er, 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; }
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 */
174 inline double radius (void) const;
175 inline const Point center (void) const;
176 inline const Pair<Vector> tangents (void) const;
178 inline Bezier approximate_bezier (double *error) const;
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;
185 inline void extents (glyphy_extents_t &extents) const;
188 double d; /* Depth */
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_) {}
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;
205 Point p0, p1, p2, p3;
209 /* Implementations */
214 inline Point::Point (const Vector &v) {
218 inline bool Point::operator == (const Point &p) const {
219 return x == p.x && y == p.y;
221 inline bool Point::operator != (const Point &p) const {
222 return !(*this == p);
224 inline Point& Point::operator+= (const Vector &v) {
229 inline Point& Point::operator-= (const Vector &v) {
234 inline const Point Point::operator+ (const Vector &v) const {
235 return Point (*this) += v;
237 inline const Point Point::operator- (const Vector &v) const {
238 return Point (*this) -= v;
240 inline const Vector Point::operator- (const Point &p) const {
241 return Vector (x - p.x, y - p.y);
244 inline const Point Point::midpoint (const Point &p) const {
245 return *this + (p - *this) / 2;
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));
252 inline double Point::distance_to_point (const Point &p) const {
253 return ((*this) - p).len ();
256 inline double Point::squared_distance_to_point (const Point &p) const {
257 return ((*this) - p).len2 ();
260 inline bool Point::is_finite (void) const {
261 return std::isfinite (x) && std::isfinite (y);
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
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);
275 inline bool Vector::operator == (const Vector &v) const {
276 return dx == v.dx && dy == v.dy;
278 inline bool Vector::operator != (const Vector &v) const {
279 return !(*this == v);
281 inline const Vector Vector::operator+ (void) const {
284 inline const Vector Vector::operator- (void) const {
285 return Vector (-dx, -dy);
287 inline Vector& Vector::operator+= (const Vector &v) {
292 inline Vector& Vector::operator-= (const Vector &v) {
297 inline Vector& Vector::operator*= (const double &s) {
302 inline Vector& Vector::operator/= (const double &s) {
307 inline const Vector Vector::operator+ (const Vector &v) const {
308 return Vector (*this) += v;
310 inline const Vector Vector::operator- (const Vector &v) const {
311 return Vector (*this) -= v;
313 inline const Vector Vector::operator* (const double &s) const {
314 return Vector (*this) *= s;
316 inline const Vector operator* (const double &s, const Vector &v) {
319 inline const Vector Vector::operator/ (const double &s) const {
320 return Vector (*this) /= s;
322 inline double Vector::operator* (const Vector &v) const { /* dot product */
323 return dx * v.dx + dy * v.dy;
325 inline const Point Vector::operator+ (const Point &p) const {
329 inline bool Vector::is_nonzero (void) const {
332 inline double Vector::len (void) const {
333 return hypot (dx, dy);
335 inline double Vector::len2 (void) const {
336 return dx * dx + dy * dy;
338 inline const Vector Vector::normalized (void) const {
340 return d ? *this / d : *this;
342 inline const Vector Vector::ortho (void) const {
343 return Vector (-dy, dx);
345 inline const Vector Vector::normal (void) const {
346 return ortho ().normalized ();
348 inline double Vector::angle (void) const {
349 return atan2 (dy, dx);
352 inline const Vector Vector::rebase (const Vector &bx,
353 const Vector &by) const {
354 return Vector (*this * bx, *this * by);
356 inline const Vector Vector::rebase (const Vector &bx) const {
357 return rebase (bx, bx.ortho ());
363 inline bool SignedVector::operator == (const SignedVector &v) const {
364 return (const Vector &)(*this) == (const Vector &)(v) && negative == v.negative;
366 inline bool SignedVector::operator != (const SignedVector &v) const {
367 return !(*this == v);
369 inline const SignedVector SignedVector::operator- (void) const {
370 return SignedVector (-(const Vector &)(*this), !negative);
376 inline const Point Line::operator+ (const Line &l) const {
377 double det = n.dx * l.n.dy - n.dy * l.n.dx;
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);
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. *************************************/
388 inline const SignedVector operator- (const Point &p, const Line &l) {
392 inline const Line Line::normalized (void) const {
394 return d ? Line (n / d, c / d) : *this;
396 inline const Vector Line::normal (void) const {
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?? ***********************/
407 inline bool Segment::contains_in_span (const Point &p) const {
411 /* shortest vector from point to line */
413 double mag = -(temp.n * Vector (p) - temp.c) / temp.n.len ();
414 Vector y (temp.n.normalized () * mag);
417 // Check if z is between p0 and p1.
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));
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));
429 inline double Segment::distance_to_point (const Point &p) const {
433 // Check if z is between p0 and p1.
435 if (contains_in_span (p))
436 return -(temp.n * Vector (p) - temp.c) / temp.n.len ();
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);
444 inline double Segment::squared_distance_to_point (const Point &p) const {
448 // Check if z is between p0 and 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);
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);
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)) ;
468 inline bool Arc::operator == (const Arc &a) const {
469 return p0 == a.p0 && p1 == a.p1 && d == a.d;
471 inline bool Arc::operator != (const Arc &a) const {
472 return !(*this == a);
476 inline const SignedVector Arc::operator- (const Point &p) const {
478 if (fabs(d) < 1e-5) {
479 Segment arc_segment (p0, p1);
480 return arc_segment - p;
482 if (wedge_contains_point (p)){
483 Vector difference = (center () - p).normalized () * fabs (p.distance_to_point (center ()) - radius ());
485 return SignedVector (difference, ((p - center ()).len () < radius ()) ^ (d < 0));
487 double d0 = p.squared_distance_to_point (p0);
488 double d1 = p.squared_distance_to_point (p1);
490 Arc other_arc (p0, p1, (1.0 + d) / (1.0 - d)); /********************************* NOT Robust. But works? *****************/
491 Vector normal = center () - (d0 < d1 ? p0 : p1) ;
493 if (normal.len() == 0)
494 return SignedVector (Vector (0, 0), true); /************************************ Check sign of this S.D. *************/
496 return SignedVector (Line (normal.dx, normal.dy, normal * Vector ((d0 < d1 ? p0 : p1))) - p, !other_arc.wedge_contains_point(p));
499 inline const SignedVector operator- (const Point &p, const Arc &a) {
505 inline double Arc::radius (void) const
507 return fabs ((p1 - p0).len () / (2 * sin2atan (d)));
510 inline const Point Arc::center (void) const
512 return (p0.midpoint (p1)) + (p1 - p0).ortho () / (2 * tan2atan (d));
515 inline const Pair<Vector> Arc::tangents (void) const
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);
525 inline Bezier Arc::approximate_bezier (double *error) const
528 Vector pp = dp.ortho ();
531 *error = dp.len () * pow (fabs (d), 5) / (54 * (1 + d*d));
533 dp *= ((1 - d*d) / 3);
536 Point p0s = p0 + dp - pp;
537 Point p1s = p1 - dp - pp;
539 return Bezier (p0, p0s, p1s, p1);
543 inline bool Arc::wedge_contains_point (const Point &p) const
545 Pair<Vector> t = tangents ();
547 return (p - p0) * t.first >= 0 && (p - p1) * t.second <= 0;
549 return (p - p0) * t.first >= 0 || (p - p1) * t.second <= 0;
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);
560 SignedVector difference = *this - p;
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);
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);
576 //SignedVector difference = *this - p;
578 if (wedge_contains_point (p) && fabs(d) > 1e-5) {
579 double answer = p.distance_to_point (center ()) - radius ();
580 return answer * answer;
582 double d1 = p.squared_distance_to_point (p0);
583 double d2 = p.squared_distance_to_point (p1);
584 return (d1 < d2 ? d1 : d2);
587 inline double Arc::extended_dist (const Point &p) const {
588 Point m = p0.lerp (.5, p1);
590 Vector pp = dp.ortho ();
591 float d2 = tan2atan (d);
592 if ((p - m) * (p1 - m) < 0)
593 return (p - p0) * (pp + dp * d2).normalized ();
595 return (p - p1) * (pp - dp * d2).normalized ();
598 inline void Arc::extents (glyphy_extents_t &extents) const {
599 glyphy_extents_clear (&extents);
600 glyphy_extents_add (&extents, &p0);
601 glyphy_extents_add (&extents, &p1);
603 double r = radius ();
604 Point p[4] = {c + r * Vector (-1, 0),
605 c + r * Vector (+1, 0),
606 c + r * Vector ( 0, -1),
607 c + r * Vector ( 0, +1)};
608 for (unsigned int i = 0; i < 4; i++)
609 if (wedge_contains_point (p[i]))
610 glyphy_extents_add (&extents, &p[i]);
616 inline const Point Bezier::point (const double &t) const {
617 Point p01 = p0.lerp (t, p1);
618 Point p12 = p1.lerp (t, p2);
619 Point p23 = p2.lerp (t, p3);
620 Point p012 = p01.lerp (t, p12);
621 Point p123 = p12.lerp (t, p23);
622 Point p0123 = p012.lerp (t, p123);
626 inline const Point Bezier::midpoint (void) const
628 Point p01 = p0.midpoint (p1);
629 Point p12 = p1.midpoint (p2);
630 Point p23 = p2.midpoint (p3);
631 Point p012 = p01.midpoint (p12);
632 Point p123 = p12.midpoint (p23);
633 Point p0123 = p012.midpoint (p123);
637 inline const Vector Bezier::tangent (const double &t) const
639 double t_2_0 = t * t;
640 double t_0_2 = (1 - t) * (1 - t);
642 double _1__4t_1_0_3t_2_0 = 1 - 4 * t + 3 * t_2_0;
643 double _2t_1_0_3t_2_0 = 2 * t - 3 * t_2_0;
645 return Vector (-3 * p0.x * t_0_2
646 +3 * p1.x * _1__4t_1_0_3t_2_0
647 +3 * p2.x * _2t_1_0_3t_2_0
650 +3 * p1.y * _1__4t_1_0_3t_2_0
651 +3 * p2.y * _2t_1_0_3t_2_0
655 inline const Vector Bezier::d_tangent (const double &t) const {
656 return Vector (6 * ((-p0.x + 3*p1.x - 3*p2.x + p3.x) * t + (p0.x - 2*p1.x + p2.x)),
657 6 * ((-p0.y + 3*p1.y - 3*p2.y + p3.y) * t + (p0.y - 2*p1.y + p2.y)));
660 inline double Bezier::curvature (const double &t) const {
661 Vector dpp = tangent (t).ortho ();
662 Vector ddp = d_tangent (t);
663 /* normal vector len squared */
664 double len = dpp.len ();
665 double curvature = (dpp * ddp) / (len * len * len);
669 inline const Pair<Bezier > Bezier::split (const double &t) const {
670 Point p01 = p0.lerp (t, p1);
671 Point p12 = p1.lerp (t, p2);
672 Point p23 = p2.lerp (t, p3);
673 Point p012 = p01.lerp (t, p12);
674 Point p123 = p12.lerp (t, p23);
675 Point p0123 = p012.lerp (t, p123);
676 return Pair<Bezier> (Bezier (p0, p01, p012, p0123),
677 Bezier (p0123, p123, p23, p3));
680 inline const Pair<Bezier > Bezier::halve (void) const
682 Point p01 = p0.midpoint (p1);
683 Point p12 = p1.midpoint (p2);
684 Point p23 = p2.midpoint (p3);
685 Point p012 = p01.midpoint (p12);
686 Point p123 = p12.midpoint (p23);
687 Point p0123 = p012.midpoint (p123);
688 return Pair<Bezier> (Bezier (p0, p01, p012, p0123),
689 Bezier (p0123, p123, p23, p3));
692 inline const Bezier Bezier::segment (const double &t0, const double &t1) const
694 Point p01 = p0.lerp (t0, p1);
695 Point p12 = p1.lerp (t0, p2);
696 Point p23 = p2.lerp (t0, p3);
697 Point p012 = p01.lerp (t0, p12);
698 Point p123 = p12.lerp (t0, p23);
699 Point p0123 = p012.lerp (t0, p123);
701 Point q01 = p0.lerp (t1, p1);
702 Point q12 = p1.lerp (t1, p2);
703 Point q23 = p2.lerp (t1, p3);
704 Point q012 = q01.lerp (t1, q12);
705 Point q123 = q12.lerp (t1, q23);
706 Point q0123 = q012.lerp (t1, q123);
708 return Bezier (p0123,
709 p0123 + (p123 - p0123) * ((t1 - t0) / (1 - t0)),
710 q0123 + (q012 - q0123) * ((t1 - t0) / t1),
714 } /* namespace Geometry */
715 } /* namespace GLyphy */
717 #pragma GCC diagnostic pop
719 #endif /* GLYPHY_GEOMETRY_HH */