1 /*
2 * Copyright 2018 Google Inc.
3 *
4 * Use of this source code is governed by a BSD-style license that can be
5 * found in the LICENSE file.
6 */
7
8 #include "include/core/SkCubicMap.h"
9
10 #include "include/private/base/SkTPin.h"
11 #include "src/base/SkVx.h"
12
13 #include <algorithm>
14 #include <cmath>
15
eval_poly(float t,float b)16 static float eval_poly(float t, float b) { return b; }
17
18 template <typename... Rest>
eval_poly(float t,float m,float b,Rest...rest)19 static float eval_poly(float t, float m, float b, Rest... rest) {
20 return eval_poly(t, std::fma(m, t, b), rest...);
21 }
22
cubic_solver(float A,float B,float C,float D)23 static float cubic_solver(float A, float B, float C, float D) {
24 #ifdef SK_DEBUG
25 auto valid = [](float t) { return t >= 0 && t <= 1; };
26 #endif
27
28 auto guess_nice_cubic_root = [](float a, float b, float c, float d) { return -d; };
29 float t = guess_nice_cubic_root(A, B, C, D);
30
31 int iters = 0;
32 const int MAX_ITERS = 8;
33 for (; iters < MAX_ITERS; ++iters) {
34 SkASSERT(valid(t));
35 float f = eval_poly(t, A, B, C, D); // f = At^3 + Bt^2 + Ct + D
36 if (std::fabs(f) <= 0.00005f) {
37 break;
38 }
39 float fp = eval_poly(t, 3*A, 2*B, C); // f' = 3At^2 + 2Bt + C
40 float fpp = eval_poly(t, 3*A + 3*A, 2*B); // f'' = 6At + 2B
41
42 float numer = 2 * fp * f;
43 float denom = std::fma(2 * fp, fp, -(f * fpp));
44
45 t -= numer / denom;
46 }
47
48 SkASSERT(valid(t));
49 return t;
50 }
51
nearly_zero(SkScalar x)52 static inline bool nearly_zero(SkScalar x) {
53 SkASSERT(x >= 0);
54 return x <= 0.0000000001f;
55 }
56
compute_t_from_x(float A,float B,float C,float x)57 static float compute_t_from_x(float A, float B, float C, float x) {
58 return cubic_solver(A, B, C, -x);
59 }
60
computeYFromX(float x) const61 float SkCubicMap::computeYFromX(float x) const {
62 x = SkTPin(x, 0.0f, 1.0f);
63
64 if (nearly_zero(x) || nearly_zero(1 - x)) {
65 return x;
66 }
67 if (fType == kLine_Type) {
68 return x;
69 }
70 float t;
71 if (fType == kCubeRoot_Type) {
72 t = std::pow(x / fCoeff[0].fX, 1.0f / 3);
73 } else {
74 t = compute_t_from_x(fCoeff[0].fX, fCoeff[1].fX, fCoeff[2].fX, x);
75 }
76 float a = fCoeff[0].fY;
77 float b = fCoeff[1].fY;
78 float c = fCoeff[2].fY;
79 float y = ((a * t + b) * t + c) * t;
80
81 return y;
82 }
83
coeff_nearly_zero(float delta)84 static inline bool coeff_nearly_zero(float delta) {
85 return std::fabs(delta) <= 0.0000001f;
86 }
87
SkCubicMap(SkPoint p1,SkPoint p2)88 SkCubicMap::SkCubicMap(SkPoint p1, SkPoint p2) {
89 // Clamp X values only (we allow Ys outside [0..1]).
90 p1.fX = std::min(std::max(p1.fX, 0.0f), 1.0f);
91 p2.fX = std::min(std::max(p2.fX, 0.0f), 1.0f);
92
93 auto s1 = skvx::float2::Load(&p1) * 3;
94 auto s2 = skvx::float2::Load(&p2) * 3;
95
96 (1 + s1 - s2).store(&fCoeff[0]);
97 (s2 - s1 - s1).store(&fCoeff[1]);
98 s1.store(&fCoeff[2]);
99
100 fType = kSolver_Type;
101 if (SkScalarNearlyEqual(p1.fX, p1.fY) && SkScalarNearlyEqual(p2.fX, p2.fY)) {
102 fType = kLine_Type;
103 } else if (coeff_nearly_zero(fCoeff[1].fX) && coeff_nearly_zero(fCoeff[2].fX)) {
104 fType = kCubeRoot_Type;
105 }
106 }
107
computeFromT(float t) const108 SkPoint SkCubicMap::computeFromT(float t) const {
109 auto a = skvx::float2::Load(&fCoeff[0]);
110 auto b = skvx::float2::Load(&fCoeff[1]);
111 auto c = skvx::float2::Load(&fCoeff[2]);
112
113 SkPoint result;
114 (((a * t + b) * t + c) * t).store(&result);
115 return result;
116 }
117