/src/skia/src/pathops/SkPathOpsCubic.cpp
Line | Count | Source (jump to first uncovered line) |
1 | | /* |
2 | | * Copyright 2012 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 | | #include "include/private/SkTPin.h" |
8 | | #include "src/core/SkGeometry.h" |
9 | | #include "src/core/SkTSort.h" |
10 | | #include "src/pathops/SkLineParameters.h" |
11 | | #include "src/pathops/SkPathOpsConic.h" |
12 | | #include "src/pathops/SkPathOpsCubic.h" |
13 | | #include "src/pathops/SkPathOpsCurve.h" |
14 | | #include "src/pathops/SkPathOpsLine.h" |
15 | | #include "src/pathops/SkPathOpsQuad.h" |
16 | | #include "src/pathops/SkPathOpsRect.h" |
17 | | |
18 | | const int SkDCubic::gPrecisionUnit = 256; // FIXME: test different values in test framework |
19 | | |
20 | 2.57M | void SkDCubic::align(int endIndex, int ctrlIndex, SkDPoint* dstPt) const { |
21 | 2.57M | if (fPts[endIndex].fX == fPts[ctrlIndex].fX) { |
22 | 571k | dstPt->fX = fPts[endIndex].fX; |
23 | 571k | } |
24 | 2.57M | if (fPts[endIndex].fY == fPts[ctrlIndex].fY) { |
25 | 365k | dstPt->fY = fPts[endIndex].fY; |
26 | 365k | } |
27 | 2.57M | } |
28 | | |
29 | | // give up when changing t no longer moves point |
30 | | // also, copy point rather than recompute it when it does change |
31 | | double SkDCubic::binarySearch(double min, double max, double axisIntercept, |
32 | 20.8M | SearchAxis xAxis) const { |
33 | 20.8M | double t = (min + max) / 2; |
34 | 20.8M | double step = (t - min) / 2; |
35 | 20.8M | SkDPoint cubicAtT = ptAtT(t); |
36 | 20.8M | double calcPos = (&cubicAtT.fX)[xAxis]; |
37 | 20.8M | double calcDist = calcPos - axisIntercept; |
38 | 993M | do { |
39 | 993M | double priorT = std::max(min, t - step); |
40 | 993M | SkDPoint lessPt = ptAtT(priorT); |
41 | 993M | if (approximately_equal_half(lessPt.fX, cubicAtT.fX) |
42 | 34.1M | && approximately_equal_half(lessPt.fY, cubicAtT.fY)) { |
43 | 14.4M | return -1; // binary search found no point at this axis intercept |
44 | 14.4M | } |
45 | 978M | double lessDist = (&lessPt.fX)[xAxis] - axisIntercept; |
46 | | #if DEBUG_CUBIC_BINARY_SEARCH |
47 | | SkDebugf("t=%1.9g calc=%1.9g dist=%1.9g step=%1.9g less=%1.9g\n", t, calcPos, calcDist, |
48 | | step, lessDist); |
49 | | #endif |
50 | 978M | double lastStep = step; |
51 | 978M | step /= 2; |
52 | 978M | if (calcDist > 0 ? calcDist > lessDist : calcDist < lessDist) { |
53 | 479M | t = priorT; |
54 | 499M | } else { |
55 | 499M | double nextT = t + lastStep; |
56 | 499M | if (nextT > max) { |
57 | 1.14k | return -1; |
58 | 1.14k | } |
59 | 499M | SkDPoint morePt = ptAtT(nextT); |
60 | 499M | if (approximately_equal_half(morePt.fX, cubicAtT.fX) |
61 | 12.1M | && approximately_equal_half(morePt.fY, cubicAtT.fY)) { |
62 | 1.25M | return -1; // binary search found no point at this axis intercept |
63 | 1.25M | } |
64 | 498M | double moreDist = (&morePt.fX)[xAxis] - axisIntercept; |
65 | 498M | if (calcDist > 0 ? calcDist <= moreDist : calcDist >= moreDist) { |
66 | 14.1M | continue; |
67 | 14.1M | } |
68 | 484M | t = nextT; |
69 | 484M | } |
70 | 963M | SkDPoint testAtT = ptAtT(t); |
71 | 963M | cubicAtT = testAtT; |
72 | 963M | calcPos = (&cubicAtT.fX)[xAxis]; |
73 | 963M | calcDist = calcPos - axisIntercept; |
74 | 977M | } while (!approximately_equal(calcPos, axisIntercept)); |
75 | 5.14M | return t; |
76 | 20.8M | } |
77 | | |
78 | | // get the rough scale of the cubic; used to determine if curvature is extreme |
79 | 985k | double SkDCubic::calcPrecision() const { |
80 | 985k | return ((fPts[1] - fPts[0]).length() |
81 | 985k | + (fPts[2] - fPts[1]).length() |
82 | 985k | + (fPts[3] - fPts[2]).length()) / gPrecisionUnit; |
83 | 985k | } |
84 | | |
85 | | /* classic one t subdivision */ |
86 | 52.6M | static void interp_cubic_coords(const double* src, double* dst, double t) { |
87 | 52.6M | double ab = SkDInterp(src[0], src[2], t); |
88 | 52.6M | double bc = SkDInterp(src[2], src[4], t); |
89 | 52.6M | double cd = SkDInterp(src[4], src[6], t); |
90 | 52.6M | double abc = SkDInterp(ab, bc, t); |
91 | 52.6M | double bcd = SkDInterp(bc, cd, t); |
92 | 52.6M | double abcd = SkDInterp(abc, bcd, t); |
93 | | |
94 | 52.6M | dst[0] = src[0]; |
95 | 52.6M | dst[2] = ab; |
96 | 52.6M | dst[4] = abc; |
97 | 52.6M | dst[6] = abcd; |
98 | 52.6M | dst[8] = bcd; |
99 | 52.6M | dst[10] = cd; |
100 | 52.6M | dst[12] = src[6]; |
101 | 52.6M | } |
102 | | |
103 | 35.3M | SkDCubicPair SkDCubic::chopAt(double t) const { |
104 | 35.3M | SkDCubicPair dst; |
105 | 35.3M | if (t == 0.5) { |
106 | 8.98M | dst.pts[0] = fPts[0]; |
107 | 8.98M | dst.pts[1].fX = (fPts[0].fX + fPts[1].fX) / 2; |
108 | 8.98M | dst.pts[1].fY = (fPts[0].fY + fPts[1].fY) / 2; |
109 | 8.98M | dst.pts[2].fX = (fPts[0].fX + 2 * fPts[1].fX + fPts[2].fX) / 4; |
110 | 8.98M | dst.pts[2].fY = (fPts[0].fY + 2 * fPts[1].fY + fPts[2].fY) / 4; |
111 | 8.98M | dst.pts[3].fX = (fPts[0].fX + 3 * (fPts[1].fX + fPts[2].fX) + fPts[3].fX) / 8; |
112 | 8.98M | dst.pts[3].fY = (fPts[0].fY + 3 * (fPts[1].fY + fPts[2].fY) + fPts[3].fY) / 8; |
113 | 8.98M | dst.pts[4].fX = (fPts[1].fX + 2 * fPts[2].fX + fPts[3].fX) / 4; |
114 | 8.98M | dst.pts[4].fY = (fPts[1].fY + 2 * fPts[2].fY + fPts[3].fY) / 4; |
115 | 8.98M | dst.pts[5].fX = (fPts[2].fX + fPts[3].fX) / 2; |
116 | 8.98M | dst.pts[5].fY = (fPts[2].fY + fPts[3].fY) / 2; |
117 | 8.98M | dst.pts[6] = fPts[3]; |
118 | 8.98M | return dst; |
119 | 8.98M | } |
120 | 26.3M | interp_cubic_coords(&fPts[0].fX, &dst.pts[0].fX, t); |
121 | 26.3M | interp_cubic_coords(&fPts[0].fY, &dst.pts[0].fY, t); |
122 | 26.3M | return dst; |
123 | 26.3M | } |
124 | | |
125 | 72.7M | void SkDCubic::Coefficients(const double* src, double* A, double* B, double* C, double* D) { |
126 | 72.7M | *A = src[6]; // d |
127 | 72.7M | *B = src[4] * 3; // 3*c |
128 | 72.7M | *C = src[2] * 3; // 3*b |
129 | 72.7M | *D = src[0]; // a |
130 | 72.7M | *A -= *D - *C + *B; // A = -a + 3*b - 3*c + d |
131 | 72.7M | *B += 3 * *D - 2 * *C; // B = 3*a - 6*b + 3*c |
132 | 72.7M | *C -= 3 * *D; // C = -3*a + 3*b |
133 | 72.7M | } |
134 | | |
135 | 0 | bool SkDCubic::endsAreExtremaInXOrY() const { |
136 | 0 | return (between(fPts[0].fX, fPts[1].fX, fPts[3].fX) |
137 | 0 | && between(fPts[0].fX, fPts[2].fX, fPts[3].fX)) |
138 | 0 | || (between(fPts[0].fY, fPts[1].fY, fPts[3].fY) |
139 | 0 | && between(fPts[0].fY, fPts[2].fY, fPts[3].fY)); |
140 | 0 | } |
141 | | |
142 | | // Do a quick reject by rotating all points relative to a line formed by |
143 | | // a pair of one cubic's points. If the 2nd cubic's points |
144 | | // are on the line or on the opposite side from the 1st cubic's 'odd man', the |
145 | | // curves at most intersect at the endpoints. |
146 | | /* if returning true, check contains true if cubic's hull collapsed, making the cubic linear |
147 | | if returning false, check contains true if the the cubic pair have only the end point in common |
148 | | */ |
149 | 96.0M | bool SkDCubic::hullIntersects(const SkDPoint* pts, int ptCount, bool* isLinear) const { |
150 | 96.0M | bool linear = true; |
151 | 96.0M | char hullOrder[4]; |
152 | 96.0M | int hullCount = convexHull(hullOrder); |
153 | 96.0M | int end1 = hullOrder[0]; |
154 | 96.0M | int hullIndex = 0; |
155 | 96.0M | const SkDPoint* endPt[2]; |
156 | 96.0M | endPt[0] = &fPts[end1]; |
157 | 341M | do { |
158 | 341M | hullIndex = (hullIndex + 1) % hullCount; |
159 | 341M | int end2 = hullOrder[hullIndex]; |
160 | 341M | endPt[1] = &fPts[end2]; |
161 | 341M | double origX = endPt[0]->fX; |
162 | 341M | double origY = endPt[0]->fY; |
163 | 341M | double adj = endPt[1]->fX - origX; |
164 | 341M | double opp = endPt[1]->fY - origY; |
165 | 341M | int oddManMask = other_two(end1, end2); |
166 | 341M | int oddMan = end1 ^ oddManMask; |
167 | 341M | double sign = (fPts[oddMan].fY - origY) * adj - (fPts[oddMan].fX - origX) * opp; |
168 | 341M | int oddMan2 = end2 ^ oddManMask; |
169 | 341M | double sign2 = (fPts[oddMan2].fY - origY) * adj - (fPts[oddMan2].fX - origX) * opp; |
170 | 341M | if (sign * sign2 < 0) { |
171 | 9.20M | continue; |
172 | 9.20M | } |
173 | 332M | if (approximately_zero(sign)) { |
174 | 23.3M | sign = sign2; |
175 | 23.3M | if (approximately_zero(sign)) { |
176 | 20.0M | continue; |
177 | 20.0M | } |
178 | 312M | } |
179 | 312M | linear = false; |
180 | 312M | bool foundOutlier = false; |
181 | 535M | for (int n = 0; n < ptCount; ++n) { |
182 | 520M | double test = (pts[n].fY - origY) * adj - (pts[n].fX - origX) * opp; |
183 | 520M | if (test * sign > 0 && !precisely_zero(test)) { |
184 | 297M | foundOutlier = true; |
185 | 297M | break; |
186 | 297M | } |
187 | 520M | } |
188 | 312M | if (!foundOutlier) { |
189 | 14.8M | return false; |
190 | 14.8M | } |
191 | 297M | endPt[0] = endPt[1]; |
192 | 297M | end1 = end2; |
193 | 327M | } while (hullIndex); |
194 | 81.2M | *isLinear = linear; |
195 | 81.2M | return true; |
196 | 96.0M | } |
197 | | |
198 | 86.8M | bool SkDCubic::hullIntersects(const SkDCubic& c2, bool* isLinear) const { |
199 | 86.8M | return hullIntersects(c2.fPts, SkDCubic::kPointCount, isLinear); |
200 | 86.8M | } |
201 | | |
202 | 9.18M | bool SkDCubic::hullIntersects(const SkDQuad& quad, bool* isLinear) const { |
203 | 9.18M | return hullIntersects(quad.fPts, SkDQuad::kPointCount, isLinear); |
204 | 9.18M | } |
205 | | |
206 | 4.31M | bool SkDCubic::hullIntersects(const SkDConic& conic, bool* isLinear) const { |
207 | | |
208 | 4.31M | return hullIntersects(conic.fPts, isLinear); |
209 | 4.31M | } |
210 | | |
211 | 3.31M | bool SkDCubic::isLinear(int startIndex, int endIndex) const { |
212 | 3.31M | if (fPts[0].approximatelyDEqual(fPts[3])) { |
213 | 5.03k | return ((const SkDQuad *) this)->isLinear(0, 2); |
214 | 5.03k | } |
215 | 3.31M | SkLineParameters lineParameters; |
216 | 3.31M | lineParameters.cubicEndPoints(*this, startIndex, endIndex); |
217 | | // FIXME: maybe it's possible to avoid this and compare non-normalized |
218 | 3.31M | lineParameters.normalize(); |
219 | 3.31M | double tiniest = std::min(std::min(std::min(std::min(std::min(std::min(std::min(fPts[0].fX, fPts[0].fY), |
220 | 3.31M | fPts[1].fX), fPts[1].fY), fPts[2].fX), fPts[2].fY), fPts[3].fX), fPts[3].fY); |
221 | 3.31M | double largest = std::max(std::max(std::max(std::max(std::max(std::max(std::max(fPts[0].fX, fPts[0].fY), |
222 | 3.31M | fPts[1].fX), fPts[1].fY), fPts[2].fX), fPts[2].fY), fPts[3].fX), fPts[3].fY); |
223 | 3.31M | largest = std::max(largest, -tiniest); |
224 | 3.31M | double distance = lineParameters.controlPtDistance(*this, 1); |
225 | 3.31M | if (!approximately_zero_when_compared_to(distance, largest)) { |
226 | 2.92M | return false; |
227 | 2.92M | } |
228 | 386k | distance = lineParameters.controlPtDistance(*this, 2); |
229 | 386k | return approximately_zero_when_compared_to(distance, largest); |
230 | 386k | } |
231 | | |
232 | | // from http://www.cs.sunysb.edu/~qin/courses/geometry/4.pdf |
233 | | // c(t) = a(1-t)^3 + 3bt(1-t)^2 + 3c(1-t)t^2 + dt^3 |
234 | | // c'(t) = -3a(1-t)^2 + 3b((1-t)^2 - 2t(1-t)) + 3c(2t(1-t) - t^2) + 3dt^2 |
235 | | // = 3(b-a)(1-t)^2 + 6(c-b)t(1-t) + 3(d-c)t^2 |
236 | 70.7M | static double derivative_at_t(const double* src, double t) { |
237 | 70.7M | double one_t = 1 - t; |
238 | 70.7M | double a = src[0]; |
239 | 70.7M | double b = src[2]; |
240 | 70.7M | double c = src[4]; |
241 | 70.7M | double d = src[6]; |
242 | 70.7M | return 3 * ((b - a) * one_t * one_t + 2 * (c - b) * t * one_t + (d - c) * t * t); |
243 | 70.7M | } |
244 | | |
245 | 2.73M | int SkDCubic::ComplexBreak(const SkPoint pointsPtr[4], SkScalar* t) { |
246 | 2.73M | SkDCubic cubic; |
247 | 2.73M | cubic.set(pointsPtr); |
248 | 2.73M | if (cubic.monotonicInX() && cubic.monotonicInY()) { |
249 | 1.73M | return 0; |
250 | 1.73M | } |
251 | 1.00M | double tt[2], ss[2]; |
252 | 1.00M | SkCubicType cubicType = SkClassifyCubic(pointsPtr, tt, ss); |
253 | 1.00M | switch (cubicType) { |
254 | 279k | case SkCubicType::kLoop: { |
255 | 279k | const double &td = tt[0], &te = tt[1], &sd = ss[0], &se = ss[1]; |
256 | 279k | if (roughly_between(0, td, sd) && roughly_between(0, te, se)) { |
257 | 12.0k | t[0] = static_cast<SkScalar>((td * se + te * sd) / (2 * sd * se)); |
258 | 12.0k | return (int) (t[0] > 0 && t[0] < 1); |
259 | 12.0k | } |
260 | 266k | } |
261 | 266k | [[fallthrough]]; // fall through if no t value found |
262 | 867k | case SkCubicType::kSerpentine: |
263 | 992k | case SkCubicType::kLocalCusp: |
264 | 995k | case SkCubicType::kCuspAtInfinity: { |
265 | 995k | double inflectionTs[2]; |
266 | 995k | int infTCount = cubic.findInflections(inflectionTs); |
267 | 995k | double maxCurvature[3]; |
268 | 995k | int roots = cubic.findMaxCurvature(maxCurvature); |
269 | | #if DEBUG_CUBIC_SPLIT |
270 | | SkDebugf("%s\n", __FUNCTION__); |
271 | | cubic.dump(); |
272 | | for (int index = 0; index < infTCount; ++index) { |
273 | | SkDebugf("inflectionsTs[%d]=%1.9g ", index, inflectionTs[index]); |
274 | | SkDPoint pt = cubic.ptAtT(inflectionTs[index]); |
275 | | SkDVector dPt = cubic.dxdyAtT(inflectionTs[index]); |
276 | | SkDLine perp = {{pt - dPt, pt + dPt}}; |
277 | | perp.dump(); |
278 | | } |
279 | | for (int index = 0; index < roots; ++index) { |
280 | | SkDebugf("maxCurvature[%d]=%1.9g ", index, maxCurvature[index]); |
281 | | SkDPoint pt = cubic.ptAtT(maxCurvature[index]); |
282 | | SkDVector dPt = cubic.dxdyAtT(maxCurvature[index]); |
283 | | SkDLine perp = {{pt - dPt, pt + dPt}}; |
284 | | perp.dump(); |
285 | | } |
286 | | #endif |
287 | 995k | if (infTCount == 2) { |
288 | 11.3k | for (int index = 0; index < roots; ++index) { |
289 | 10.9k | if (between(inflectionTs[0], maxCurvature[index], inflectionTs[1])) { |
290 | 9.85k | t[0] = maxCurvature[index]; |
291 | 9.85k | return (int) (t[0] > 0 && t[0] < 1); |
292 | 9.85k | } |
293 | 10.9k | } |
294 | 985k | } else { |
295 | 985k | int resultCount = 0; |
296 | | // FIXME: constant found through experimentation -- maybe there's a better way.... |
297 | 985k | double precision = cubic.calcPrecision() * 2; |
298 | 2.17M | for (int index = 0; index < roots; ++index) { |
299 | 1.18M | double testT = maxCurvature[index]; |
300 | 1.18M | if (0 >= testT || testT >= 1) { |
301 | 402k | continue; |
302 | 402k | } |
303 | | // don't call dxdyAtT since we want (0,0) results |
304 | 783k | SkDVector dPt = { derivative_at_t(&cubic.fPts[0].fX, testT), |
305 | 783k | derivative_at_t(&cubic.fPts[0].fY, testT) }; |
306 | 783k | double dPtLen = dPt.length(); |
307 | 783k | if (dPtLen < precision) { |
308 | 141k | t[resultCount++] = testT; |
309 | 141k | } |
310 | 783k | } |
311 | 985k | if (!resultCount && infTCount == 1) { |
312 | 227k | t[0] = inflectionTs[0]; |
313 | 227k | resultCount = (int) (t[0] > 0 && t[0] < 1); |
314 | 227k | } |
315 | 985k | return resultCount; |
316 | 985k | } |
317 | 416 | break; |
318 | 416 | } |
319 | 95 | default: |
320 | 95 | break; |
321 | 511 | } |
322 | 511 | return 0; |
323 | 511 | } |
324 | | |
325 | 127M | bool SkDCubic::monotonicInX() const { |
326 | 127M | return precisely_between(fPts[0].fX, fPts[1].fX, fPts[3].fX) |
327 | 123M | && precisely_between(fPts[0].fX, fPts[2].fX, fPts[3].fX); |
328 | 127M | } |
329 | | |
330 | 127M | bool SkDCubic::monotonicInY() const { |
331 | 127M | return precisely_between(fPts[0].fY, fPts[1].fY, fPts[3].fY) |
332 | 115M | && precisely_between(fPts[0].fY, fPts[2].fY, fPts[3].fY); |
333 | 127M | } |
334 | | |
335 | 30.6M | void SkDCubic::otherPts(int index, const SkDPoint* o1Pts[kPointCount - 1]) const { |
336 | 30.6M | int offset = (int) !SkToBool(index); |
337 | 30.6M | o1Pts[0] = &fPts[offset]; |
338 | 30.6M | o1Pts[1] = &fPts[++offset]; |
339 | 30.6M | o1Pts[2] = &fPts[++offset]; |
340 | 30.6M | } |
341 | | |
342 | | int SkDCubic::searchRoots(double extremeTs[6], int extrema, double axisIntercept, |
343 | 13.4M | SearchAxis xAxis, double* validRoots) const { |
344 | 13.4M | extrema += findInflections(&extremeTs[extrema]); |
345 | 13.4M | extremeTs[extrema++] = 0; |
346 | 13.4M | extremeTs[extrema] = 1; |
347 | 13.4M | SkASSERT(extrema < 6); |
348 | 13.4M | SkTQSort(extremeTs, extremeTs + extrema + 1); |
349 | 13.4M | int validCount = 0; |
350 | 42.3M | for (int index = 0; index < extrema; ) { |
351 | 28.8M | double min = extremeTs[index]; |
352 | 28.8M | double max = extremeTs[++index]; |
353 | 28.8M | if (min == max) { |
354 | 8.00M | continue; |
355 | 8.00M | } |
356 | 20.8M | double newT = binarySearch(min, max, axisIntercept, xAxis); |
357 | 20.8M | if (newT >= 0) { |
358 | 5.14M | if (validCount >= 3) { |
359 | 484 | return 0; |
360 | 484 | } |
361 | 5.14M | validRoots[validCount++] = newT; |
362 | 5.14M | } |
363 | 20.8M | } |
364 | 13.4M | return validCount; |
365 | 13.4M | } |
366 | | |
367 | | // cubic roots |
368 | | |
369 | | static const double PI = 3.141592653589793; |
370 | | |
371 | | // from SkGeometry.cpp (and Numeric Solutions, 5.6) |
372 | 73.7M | int SkDCubic::RootsValidT(double A, double B, double C, double D, double t[3]) { |
373 | 73.7M | double s[3]; |
374 | 73.7M | int realRoots = RootsReal(A, B, C, D, s); |
375 | 73.7M | int foundRoots = SkDQuad::AddValidTs(s, realRoots, t); |
376 | 271M | for (int index = 0; index < realRoots; ++index) { |
377 | 197M | double tValue = s[index]; |
378 | 197M | if (!approximately_one_or_less(tValue) && between(1, tValue, 1.00005)) { |
379 | 501k | for (int idx2 = 0; idx2 < foundRoots; ++idx2) { |
380 | 140k | if (approximately_equal(t[idx2], 1)) { |
381 | 30.0k | goto nextRoot; |
382 | 30.0k | } |
383 | 140k | } |
384 | 361k | SkASSERT(foundRoots < 3); |
385 | 361k | t[foundRoots++] = 1; |
386 | 197M | } else if (!approximately_zero_or_more(tValue) && between(-0.00005, tValue, 0)) { |
387 | 1.72M | for (int idx2 = 0; idx2 < foundRoots; ++idx2) { |
388 | 814k | if (approximately_equal(t[idx2], 0)) { |
389 | 71.5k | goto nextRoot; |
390 | 71.5k | } |
391 | 814k | } |
392 | 914k | SkASSERT(foundRoots < 3); |
393 | 914k | t[foundRoots++] = 0; |
394 | 914k | } |
395 | 197M | nextRoot: |
396 | 197M | ; |
397 | 197M | } |
398 | 73.7M | return foundRoots; |
399 | 73.7M | } |
400 | | |
401 | 73.7M | int SkDCubic::RootsReal(double A, double B, double C, double D, double s[3]) { |
402 | | #ifdef SK_DEBUG |
403 | | #if ONE_OFF_DEBUG && ONE_OFF_DEBUG_MATHEMATICA |
404 | | // create a string mathematica understands |
405 | | // GDB set print repe 15 # if repeated digits is a bother |
406 | | // set print elements 400 # if line doesn't fit |
407 | | char str[1024]; |
408 | | sk_bzero(str, sizeof(str)); |
409 | | SK_SNPRINTF(str, sizeof(str), "Solve[%1.19g x^3 + %1.19g x^2 + %1.19g x + %1.19g == 0, x]", |
410 | | A, B, C, D); |
411 | | SkPathOpsDebug::MathematicaIze(str, sizeof(str)); |
412 | | SkDebugf("%s\n", str); |
413 | | #endif |
414 | | #endif |
415 | 73.7M | if (approximately_zero(A) |
416 | 7.30M | && approximately_zero_when_compared_to(A, B) |
417 | 268k | && approximately_zero_when_compared_to(A, C) |
418 | 264k | && approximately_zero_when_compared_to(A, D)) { // we're just a quadratic |
419 | 259k | return SkDQuad::RootsReal(B, C, D, s); |
420 | 259k | } |
421 | 73.4M | if (approximately_zero_when_compared_to(D, A) |
422 | 9.30M | && approximately_zero_when_compared_to(D, B) |
423 | 9.23M | && approximately_zero_when_compared_to(D, C)) { // 0 is one root |
424 | 6.57M | int num = SkDQuad::RootsReal(A, B, C, s); |
425 | 18.2M | for (int i = 0; i < num; ++i) { |
426 | 12.1M | if (approximately_zero(s[i])) { |
427 | 432k | return num; |
428 | 432k | } |
429 | 12.1M | } |
430 | 6.14M | s[num++] = 0; |
431 | 6.14M | return num; |
432 | 66.9M | } |
433 | 66.9M | if (approximately_zero(A + B + C + D)) { // 1 is one root |
434 | 10.6M | int num = SkDQuad::RootsReal(A, A + B, -D, s); |
435 | 30.3M | for (int i = 0; i < num; ++i) { |
436 | 20.0M | if (AlmostDequalUlps(s[i], 1)) { |
437 | 322k | return num; |
438 | 322k | } |
439 | 20.0M | } |
440 | 10.3M | s[num++] = 1; |
441 | 10.3M | return num; |
442 | 56.2M | } |
443 | 56.2M | double a, b, c; |
444 | 56.2M | { |
445 | 56.2M | double invA = 1 / A; |
446 | 56.2M | a = B * invA; |
447 | 56.2M | b = C * invA; |
448 | 56.2M | c = D * invA; |
449 | 56.2M | } |
450 | 56.2M | double a2 = a * a; |
451 | 56.2M | double Q = (a2 - b * 3) / 9; |
452 | 56.2M | double R = (2 * a2 * a - 9 * a * b + 27 * c) / 54; |
453 | 56.2M | double R2 = R * R; |
454 | 56.2M | double Q3 = Q * Q * Q; |
455 | 56.2M | double R2MinusQ3 = R2 - Q3; |
456 | 56.2M | double adiv3 = a / 3; |
457 | 56.2M | double r; |
458 | 56.2M | double* roots = s; |
459 | 56.2M | if (R2MinusQ3 < 0) { // we have 3 real roots |
460 | | // the divide/root can, due to finite precisions, be slightly outside of -1...1 |
461 | 45.5M | double theta = acos(SkTPin(R / sqrt(Q3), -1., 1.)); |
462 | 45.5M | double neg2RootQ = -2 * sqrt(Q); |
463 | | |
464 | 45.5M | r = neg2RootQ * cos(theta / 3) - adiv3; |
465 | 45.5M | *roots++ = r; |
466 | | |
467 | 45.5M | r = neg2RootQ * cos((theta + 2 * PI) / 3) - adiv3; |
468 | 45.5M | if (!AlmostDequalUlps(s[0], r)) { |
469 | 45.5M | *roots++ = r; |
470 | 45.5M | } |
471 | 45.5M | r = neg2RootQ * cos((theta - 2 * PI) / 3) - adiv3; |
472 | 45.5M | if (!AlmostDequalUlps(s[0], r) && (roots - s == 1 || !AlmostDequalUlps(s[1], r))) { |
473 | 45.4M | *roots++ = r; |
474 | 45.4M | } |
475 | 10.7M | } else { // we have 1 real root |
476 | 10.7M | double sqrtR2MinusQ3 = sqrt(R2MinusQ3); |
477 | 10.7M | A = fabs(R) + sqrtR2MinusQ3; |
478 | 10.7M | A = SkDCubeRoot(A); |
479 | 10.7M | if (R > 0) { |
480 | 5.21M | A = -A; |
481 | 5.21M | } |
482 | 10.7M | if (A != 0) { |
483 | 10.7M | A += Q / A; |
484 | 10.7M | } |
485 | 10.7M | r = A - adiv3; |
486 | 10.7M | *roots++ = r; |
487 | 10.7M | if (AlmostDequalUlps((double) R2, (double) Q3)) { |
488 | 1.34M | r = -A / 2 - adiv3; |
489 | 1.34M | if (!AlmostDequalUlps(s[0], r)) { |
490 | 1.33M | *roots++ = r; |
491 | 1.33M | } |
492 | 1.34M | } |
493 | 10.7M | } |
494 | 56.2M | return static_cast<int>(roots - s); |
495 | 56.2M | } |
496 | | |
497 | | // OPTIMIZE? compute t^2, t(1-t), and (1-t)^2 and pass them to another version of derivative at t? |
498 | 34.5M | SkDVector SkDCubic::dxdyAtT(double t) const { |
499 | 34.5M | SkDVector result = { derivative_at_t(&fPts[0].fX, t), derivative_at_t(&fPts[0].fY, t) }; |
500 | 34.5M | if (result.fX == 0 && result.fY == 0) { |
501 | 398k | if (t == 0) { |
502 | 151k | result = fPts[2] - fPts[0]; |
503 | 247k | } else if (t == 1) { |
504 | 242k | result = fPts[3] - fPts[1]; |
505 | 5.35k | } else { |
506 | | // incomplete |
507 | 5.35k | SkDebugf("!c"); |
508 | 5.35k | } |
509 | 398k | if (result.fX == 0 && result.fY == 0 && zero_or_one(t)) { |
510 | 2.41k | result = fPts[3] - fPts[0]; |
511 | 2.41k | } |
512 | 398k | } |
513 | 34.5M | return result; |
514 | 34.5M | } |
515 | | |
516 | | // OPTIMIZE? share code with formulate_F1DotF2 |
517 | 19.3M | int SkDCubic::findInflections(double tValues[]) const { |
518 | 19.3M | double Ax = fPts[1].fX - fPts[0].fX; |
519 | 19.3M | double Ay = fPts[1].fY - fPts[0].fY; |
520 | 19.3M | double Bx = fPts[2].fX - 2 * fPts[1].fX + fPts[0].fX; |
521 | 19.3M | double By = fPts[2].fY - 2 * fPts[1].fY + fPts[0].fY; |
522 | 19.3M | double Cx = fPts[3].fX + 3 * (fPts[1].fX - fPts[2].fX) - fPts[0].fX; |
523 | 19.3M | double Cy = fPts[3].fY + 3 * (fPts[1].fY - fPts[2].fY) - fPts[0].fY; |
524 | 19.3M | return SkDQuad::RootsValidT(Bx * Cy - By * Cx, Ax * Cy - Ay * Cx, Ax * By - Ay * Bx, tValues); |
525 | 19.3M | } |
526 | | |
527 | 1.99M | static void formulate_F1DotF2(const double src[], double coeff[4]) { |
528 | 1.99M | double a = src[2] - src[0]; |
529 | 1.99M | double b = src[4] - 2 * src[2] + src[0]; |
530 | 1.99M | double c = src[6] + 3 * (src[2] - src[4]) - src[0]; |
531 | 1.99M | coeff[0] = c * c; |
532 | 1.99M | coeff[1] = 3 * b * c; |
533 | 1.99M | coeff[2] = 2 * b * b + c * a; |
534 | 1.99M | coeff[3] = a * b; |
535 | 1.99M | } |
536 | | |
537 | | /** SkDCubic'(t) = At^2 + Bt + C, where |
538 | | A = 3(-a + 3(b - c) + d) |
539 | | B = 6(a - 2b + c) |
540 | | C = 3(b - a) |
541 | | Solve for t, keeping only those that fit between 0 < t < 1 |
542 | | */ |
543 | 33.7M | int SkDCubic::FindExtrema(const double src[], double tValues[2]) { |
544 | | // we divide A,B,C by 3 to simplify |
545 | 33.7M | double a = src[0]; |
546 | 33.7M | double b = src[2]; |
547 | 33.7M | double c = src[4]; |
548 | 33.7M | double d = src[6]; |
549 | 33.7M | double A = d - a + 3 * (b - c); |
550 | 33.7M | double B = 2 * (a - b - b + c); |
551 | 33.7M | double C = b - a; |
552 | | |
553 | 33.7M | return SkDQuad::RootsValidT(A, B, C, tValues); |
554 | 33.7M | } |
555 | | |
556 | | /* from SkGeometry.cpp |
557 | | Looking for F' dot F'' == 0 |
558 | | |
559 | | A = b - a |
560 | | B = c - 2b + a |
561 | | C = d - 3c + 3b - a |
562 | | |
563 | | F' = 3Ct^2 + 6Bt + 3A |
564 | | F'' = 6Ct + 6B |
565 | | |
566 | | F' dot F'' -> CCt^3 + 3BCt^2 + (2BB + CA)t + AB |
567 | | */ |
568 | 995k | int SkDCubic::findMaxCurvature(double tValues[]) const { |
569 | 995k | double coeffX[4], coeffY[4]; |
570 | 995k | int i; |
571 | 995k | formulate_F1DotF2(&fPts[0].fX, coeffX); |
572 | 995k | formulate_F1DotF2(&fPts[0].fY, coeffY); |
573 | 4.97M | for (i = 0; i < 4; i++) { |
574 | 3.98M | coeffX[i] = coeffX[i] + coeffY[i]; |
575 | 3.98M | } |
576 | 995k | return RootsValidT(coeffX[0], coeffX[1], coeffX[2], coeffX[3], tValues); |
577 | 995k | } |
578 | | |
579 | 2.71G | SkDPoint SkDCubic::ptAtT(double t) const { |
580 | 2.71G | if (0 == t) { |
581 | 24.8M | return fPts[0]; |
582 | 24.8M | } |
583 | 2.68G | if (1 == t) { |
584 | 33.6M | return fPts[3]; |
585 | 33.6M | } |
586 | 2.65G | double one_t = 1 - t; |
587 | 2.65G | double one_t2 = one_t * one_t; |
588 | 2.65G | double a = one_t2 * one_t; |
589 | 2.65G | double b = 3 * one_t2 * t; |
590 | 2.65G | double t2 = t * t; |
591 | 2.65G | double c = 3 * one_t * t2; |
592 | 2.65G | double d = t2 * t; |
593 | 2.65G | SkDPoint result = {a * fPts[0].fX + b * fPts[1].fX + c * fPts[2].fX + d * fPts[3].fX, |
594 | 2.65G | a * fPts[0].fY + b * fPts[1].fY + c * fPts[2].fY + d * fPts[3].fY}; |
595 | 2.65G | return result; |
596 | 2.65G | } |
597 | | |
598 | | /* |
599 | | Given a cubic c, t1, and t2, find a small cubic segment. |
600 | | |
601 | | The new cubic is defined as points A, B, C, and D, where |
602 | | s1 = 1 - t1 |
603 | | s2 = 1 - t2 |
604 | | A = c[0]*s1*s1*s1 + 3*c[1]*s1*s1*t1 + 3*c[2]*s1*t1*t1 + c[3]*t1*t1*t1 |
605 | | D = c[0]*s2*s2*s2 + 3*c[1]*s2*s2*t2 + 3*c[2]*s2*t2*t2 + c[3]*t2*t2*t2 |
606 | | |
607 | | We don't have B or C. So We define two equations to isolate them. |
608 | | First, compute two reference T values 1/3 and 2/3 from t1 to t2: |
609 | | |
610 | | c(at (2*t1 + t2)/3) == E |
611 | | c(at (t1 + 2*t2)/3) == F |
612 | | |
613 | | Next, compute where those values must be if we know the values of B and C: |
614 | | |
615 | | _12 = A*2/3 + B*1/3 |
616 | | 12_ = A*1/3 + B*2/3 |
617 | | _23 = B*2/3 + C*1/3 |
618 | | 23_ = B*1/3 + C*2/3 |
619 | | _34 = C*2/3 + D*1/3 |
620 | | 34_ = C*1/3 + D*2/3 |
621 | | _123 = (A*2/3 + B*1/3)*2/3 + (B*2/3 + C*1/3)*1/3 = A*4/9 + B*4/9 + C*1/9 |
622 | | 123_ = (A*1/3 + B*2/3)*1/3 + (B*1/3 + C*2/3)*2/3 = A*1/9 + B*4/9 + C*4/9 |
623 | | _234 = (B*2/3 + C*1/3)*2/3 + (C*2/3 + D*1/3)*1/3 = B*4/9 + C*4/9 + D*1/9 |
624 | | 234_ = (B*1/3 + C*2/3)*1/3 + (C*1/3 + D*2/3)*2/3 = B*1/9 + C*4/9 + D*4/9 |
625 | | _1234 = (A*4/9 + B*4/9 + C*1/9)*2/3 + (B*4/9 + C*4/9 + D*1/9)*1/3 |
626 | | = A*8/27 + B*12/27 + C*6/27 + D*1/27 |
627 | | = E |
628 | | 1234_ = (A*1/9 + B*4/9 + C*4/9)*1/3 + (B*1/9 + C*4/9 + D*4/9)*2/3 |
629 | | = A*1/27 + B*6/27 + C*12/27 + D*8/27 |
630 | | = F |
631 | | E*27 = A*8 + B*12 + C*6 + D |
632 | | F*27 = A + B*6 + C*12 + D*8 |
633 | | |
634 | | Group the known values on one side: |
635 | | |
636 | | M = E*27 - A*8 - D = B*12 + C* 6 |
637 | | N = F*27 - A - D*8 = B* 6 + C*12 |
638 | | M*2 - N = B*18 |
639 | | N*2 - M = C*18 |
640 | | B = (M*2 - N)/18 |
641 | | C = (N*2 - M)/18 |
642 | | */ |
643 | | |
644 | 639M | static double interp_cubic_coords(const double* src, double t) { |
645 | 639M | double ab = SkDInterp(src[0], src[2], t); |
646 | 639M | double bc = SkDInterp(src[2], src[4], t); |
647 | 639M | double cd = SkDInterp(src[4], src[6], t); |
648 | 639M | double abc = SkDInterp(ab, bc, t); |
649 | 639M | double bcd = SkDInterp(bc, cd, t); |
650 | 639M | double abcd = SkDInterp(abc, bcd, t); |
651 | 639M | return abcd; |
652 | 639M | } |
653 | | |
654 | 128M | SkDCubic SkDCubic::subDivide(double t1, double t2) const { |
655 | 128M | if (t1 == 0 || t2 == 1) { |
656 | 48.8M | if (t1 == 0 && t2 == 1) { |
657 | 14.0M | return *this; |
658 | 14.0M | } |
659 | 34.8M | SkDCubicPair pair = chopAt(t1 == 0 ? t2 : t1); |
660 | 18.2M | SkDCubic dst = t1 == 0 ? pair.first() : pair.second(); |
661 | 34.8M | return dst; |
662 | 34.8M | } |
663 | 79.9M | SkDCubic dst; |
664 | 79.9M | double ax = dst[0].fX = interp_cubic_coords(&fPts[0].fX, t1); |
665 | 79.9M | double ay = dst[0].fY = interp_cubic_coords(&fPts[0].fY, t1); |
666 | 79.9M | double ex = interp_cubic_coords(&fPts[0].fX, (t1*2+t2)/3); |
667 | 79.9M | double ey = interp_cubic_coords(&fPts[0].fY, (t1*2+t2)/3); |
668 | 79.9M | double fx = interp_cubic_coords(&fPts[0].fX, (t1+t2*2)/3); |
669 | 79.9M | double fy = interp_cubic_coords(&fPts[0].fY, (t1+t2*2)/3); |
670 | 79.9M | double dx = dst[3].fX = interp_cubic_coords(&fPts[0].fX, t2); |
671 | 79.9M | double dy = dst[3].fY = interp_cubic_coords(&fPts[0].fY, t2); |
672 | 79.9M | double mx = ex * 27 - ax * 8 - dx; |
673 | 79.9M | double my = ey * 27 - ay * 8 - dy; |
674 | 79.9M | double nx = fx * 27 - ax - dx * 8; |
675 | 79.9M | double ny = fy * 27 - ay - dy * 8; |
676 | 79.9M | /* bx = */ dst[1].fX = (mx * 2 - nx) / 18; |
677 | 79.9M | /* by = */ dst[1].fY = (my * 2 - ny) / 18; |
678 | 79.9M | /* cx = */ dst[2].fX = (nx * 2 - mx) / 18; |
679 | 79.9M | /* cy = */ dst[2].fY = (ny * 2 - my) / 18; |
680 | | // FIXME: call align() ? |
681 | 79.9M | return dst; |
682 | 79.9M | } |
683 | | |
684 | | void SkDCubic::subDivide(const SkDPoint& a, const SkDPoint& d, |
685 | 6.15M | double t1, double t2, SkDPoint dst[2]) const { |
686 | 6.15M | SkASSERT(t1 != t2); |
687 | | // this approach assumes that the control points computed directly are accurate enough |
688 | 6.15M | SkDCubic sub = subDivide(t1, t2); |
689 | 6.15M | dst[0] = sub[1] + (a - sub[0]); |
690 | 6.15M | dst[1] = sub[2] + (d - sub[3]); |
691 | 6.15M | if (t1 == 0 || t2 == 0) { |
692 | 760k | align(0, 1, t1 == 0 ? &dst[0] : &dst[1]); |
693 | 1.25M | } |
694 | 6.15M | if (t1 == 1 || t2 == 1) { |
695 | 952k | align(3, 2, t1 == 1 ? &dst[0] : &dst[1]); |
696 | 1.31M | } |
697 | 6.15M | if (AlmostBequalUlps(dst[0].fX, a.fX)) { |
698 | 247k | dst[0].fX = a.fX; |
699 | 247k | } |
700 | 6.15M | if (AlmostBequalUlps(dst[0].fY, a.fY)) { |
701 | 483k | dst[0].fY = a.fY; |
702 | 483k | } |
703 | 6.15M | if (AlmostBequalUlps(dst[1].fX, d.fX)) { |
704 | 421k | dst[1].fX = d.fX; |
705 | 421k | } |
706 | 6.15M | if (AlmostBequalUlps(dst[1].fY, d.fY)) { |
707 | 599k | dst[1].fY = d.fY; |
708 | 599k | } |
709 | 6.15M | } |
710 | | |
711 | 477k | bool SkDCubic::toFloatPoints(SkPoint* pts) const { |
712 | 477k | const double* dCubic = &fPts[0].fX; |
713 | 477k | SkScalar* cubic = &pts[0].fX; |
714 | 4.29M | for (int index = 0; index < kPointCount * 2; ++index) { |
715 | 3.81M | cubic[index] = SkDoubleToScalar(dCubic[index]); |
716 | 3.81M | if (SkScalarAbs(cubic[index]) < FLT_EPSILON_ORDERABLE_ERR) { |
717 | 48.3k | cubic[index] = 0; |
718 | 48.3k | } |
719 | 3.81M | } |
720 | 477k | return SkScalarsAreFinite(&pts->fX, kPointCount * 2); |
721 | 477k | } |
722 | | |
723 | 0 | double SkDCubic::top(const SkDCubic& dCurve, double startT, double endT, SkDPoint*topPt) const { |
724 | 0 | double extremeTs[2]; |
725 | 0 | double topT = -1; |
726 | 0 | int roots = SkDCubic::FindExtrema(&fPts[0].fY, extremeTs); |
727 | 0 | for (int index = 0; index < roots; ++index) { |
728 | 0 | double t = startT + (endT - startT) * extremeTs[index]; |
729 | 0 | SkDPoint mid = dCurve.ptAtT(t); |
730 | 0 | if (topPt->fY > mid.fY || (topPt->fY == mid.fY && topPt->fX > mid.fX)) { |
731 | 0 | topT = t; |
732 | 0 | *topPt = mid; |
733 | 0 | } |
734 | 0 | } |
735 | 0 | return topT; |
736 | 0 | } |
737 | | |
738 | 44.6M | int SkTCubic::intersectRay(SkIntersections* i, const SkDLine& line) const { |
739 | 44.6M | return i->intersectRay(fCubic, line); |
740 | 44.6M | } |
741 | | |
742 | 2.20M | bool SkTCubic::hullIntersects(const SkDQuad& quad, bool* isLinear) const { |
743 | 2.20M | return quad.hullIntersects(fCubic, isLinear); |
744 | 2.20M | } |
745 | | |
746 | 2.65M | bool SkTCubic::hullIntersects(const SkDConic& conic, bool* isLinear) const { |
747 | 2.65M | return conic.hullIntersects(fCubic, isLinear); |
748 | 2.65M | } |
749 | | |
750 | 122M | void SkTCubic::setBounds(SkDRect* rect) const { |
751 | 122M | rect->setBounds(fCubic); |
752 | 122M | } |