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