1 /*
2  * Copyright (C) 2008 Apple Inc. All Rights Reserved.
3  *
4  * Redistribution and use in source and binary forms, with or without
5  * modification, are permitted provided that the following conditions
6  * are met:
7  * 1. Redistributions of source code must retain the above copyright
8  *    notice, this list of conditions and the following disclaimer.
9  * 2. Redistributions in binary form must reproduce the above copyright
10  *    notice, this list of conditions and the following disclaimer in the
11  *    documentation and/or other materials provided with the distribution.
12  *
13  * THIS SOFTWARE IS PROVIDED BY APPLE INC. ``AS IS'' AND ANY
14  * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
15  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
16  * PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL APPLE INC. OR
17  * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
18  * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
19  * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
20  * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY
21  * OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
22  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
23  * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
24  */
25 
26 #pragma once
27 
28 #include <cmath>
29 #include <tuple>
30 
31 namespace mbgl {
32 namespace util {
33 
34 struct UnitBezier {
35     // Calculate the polynomial coefficients, implicit first and last control points are (0,0) and (1,1).
UnitBeziermbgl::util::UnitBezier36     constexpr UnitBezier(double p1x, double p1y, double p2x, double p2y)
37         : cx(3.0 * p1x)
38         , bx(3.0 * (p2x - p1x) - (3.0 * p1x))
39         , ax(1.0 - (3.0 * p1x) - (3.0 * (p2x - p1x) - (3.0 * p1x)))
40         , cy(3.0 * p1y)
41         , by(3.0 * (p2y - p1y) - (3.0 * p1y))
42         , ay(1.0 - (3.0 * p1y) - (3.0 * (p2y - p1y) - (3.0 * p1y))) {
43     }
44 
getP1mbgl::util::UnitBezier45     std::pair<double, double> getP1() const {
46         return { cx / 3.0, cy / 3.0 };
47     }
48 
getP2mbgl::util::UnitBezier49     std::pair<double, double> getP2() const {
50         return {
51             (bx + (3.0 * cx / 3.0) + cx) / 3.0,
52             (by + (3.0 * cy / 3.0) + cy) / 3.0,
53         };
54     }
55 
sampleCurveXmbgl::util::UnitBezier56     double sampleCurveX(double t) const {
57         // `ax t^3 + bx t^2 + cx t' expanded using Horner's rule.
58         return ((ax * t + bx) * t + cx) * t;
59     }
60 
sampleCurveYmbgl::util::UnitBezier61     double sampleCurveY(double t) const {
62         return ((ay * t + by) * t + cy) * t;
63     }
64 
sampleCurveDerivativeXmbgl::util::UnitBezier65     double sampleCurveDerivativeX(double t) const {
66         return (3.0 * ax * t + 2.0 * bx) * t + cx;
67     }
68 
69     // Given an x value, find a parametric value it came from.
solveCurveXmbgl::util::UnitBezier70     double solveCurveX(double x, double epsilon) const {
71         double t0;
72         double t1;
73         double t2;
74         double x2;
75         double d2;
76         int i;
77 
78         // First try a few iterations of Newton's method -- normally very fast.
79         for (t2 = x, i = 0; i < 8; ++i) {
80             x2 = sampleCurveX(t2) - x;
81             if (fabs (x2) < epsilon)
82                 return t2;
83             d2 = sampleCurveDerivativeX(t2);
84             if (fabs(d2) < 1e-6)
85                 break;
86             t2 = t2 - x2 / d2;
87         }
88 
89         // Fall back to the bisection method for reliability.
90         t0 = 0.0;
91         t1 = 1.0;
92         t2 = x;
93 
94         if (t2 < t0)
95             return t0;
96         if (t2 > t1)
97             return t1;
98 
99         while (t0 < t1) {
100             x2 = sampleCurveX(t2);
101             if (fabs(x2 - x) < epsilon)
102                 return t2;
103             if (x > x2)
104                 t0 = t2;
105             else
106                 t1 = t2;
107             t2 = (t1 - t0) * .5 + t0;
108         }
109 
110         // Failure.
111         return t2;
112     }
113 
solvembgl::util::UnitBezier114     double solve(double x, double epsilon) const {
115         return sampleCurveY(solveCurveX(x, epsilon));
116     }
117 
operator ==mbgl::util::UnitBezier118     bool operator==(const UnitBezier& rhs) const {
119         return std::tie(cx, bx, ax, cy, by, ay) ==
120             std::tie(rhs.cx, rhs.bx, rhs.ax, rhs.cy, rhs.by, rhs.ay);
121     }
122 
123 private:
124     const double cx;
125     const double bx;
126     const double ax;
127 
128     const double cy;
129     const double by;
130     const double ay;
131 };
132 
133 } // namespace util
134 } // namespace mbgl
135