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, c2.kPointCount, isLinear);
199}
200
201bool SkDCubic::hullIntersects(const SkDQuad& quad, bool* isLinear) const {
202 return hullIntersects(quad.fPts, quad.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 // 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 }
317 default:
318 ;
319 }
320 return 0;
321}
322
323bool SkDCubic::monotonicInX() const {
324 return precisely_between(fPts[0].fX, fPts[1].fX, fPts[3].fX)
325 && precisely_between(fPts[0].fX, fPts[2].fX, fPts[3].fX);
326}
327
328bool SkDCubic::monotonicInY() const {
329 return precisely_between(fPts[0].fY, fPts[1].fY, fPts[3].fY)
330 && precisely_between(fPts[0].fY, fPts[2].fY, fPts[3].fY);
331}
332
333void SkDCubic::otherPts(int index, const SkDPoint* o1Pts[kPointCount - 1]) const {
334 int offset = (int) !SkToBool(index);
335 o1Pts[0] = &fPts[offset];
336 o1Pts[1] = &fPts[++offset];
337 o1Pts[2] = &fPts[++offset];
338}
339
340int SkDCubic::searchRoots(double extremeTs[6], int extrema, double axisIntercept,
341 SearchAxis xAxis, double* validRoots) const {
342 extrema += findInflections(&extremeTs[extrema]);
343 extremeTs[extrema++] = 0;
344 extremeTs[extrema] = 1;
345 SkASSERT(extrema < 6);
346 SkTQSort(extremeTs, extremeTs + extrema);
347 int validCount = 0;
348 for (int index = 0; index < extrema; ) {
349 double min = extremeTs[index];
350 double max = extremeTs[++index];
351 if (min == max) {
352 continue;
353 }
354 double newT = binarySearch(min, max, axisIntercept, xAxis);
355 if (newT >= 0) {
356 if (validCount >= 3) {
357 return 0;
358 }
359 validRoots[validCount++] = newT;
360 }
361 }
362 return validCount;
363}
364
365// cubic roots
366
367static const double PI = 3.141592653589793;
368
369// from SkGeometry.cpp (and Numeric Solutions, 5.6)
370int SkDCubic::RootsValidT(double A, double B, double C, double D, double t[3]) {
371 double s[3];
372 int realRoots = RootsReal(A, B, C, D, s);
373 int foundRoots = SkDQuad::AddValidTs(s, realRoots, t);
374 for (int index = 0; index < realRoots; ++index) {
375 double tValue = s[index];
376 if (!approximately_one_or_less(tValue) && between(1, tValue, 1.00005)) {
377 for (int idx2 = 0; idx2 < foundRoots; ++idx2) {
378 if (approximately_equal(t[idx2], 1)) {
379 goto nextRoot;
380 }
381 }
382 SkASSERT(foundRoots < 3);
383 t[foundRoots++] = 1;
384 } else if (!approximately_zero_or_more(tValue) && between(-0.00005, tValue, 0)) {
385 for (int idx2 = 0; idx2 < foundRoots; ++idx2) {
386 if (approximately_equal(t[idx2], 0)) {
387 goto nextRoot;
388 }
389 }
390 SkASSERT(foundRoots < 3);
391 t[foundRoots++] = 0;
392 }
393nextRoot:
394 ;
395 }
396 return foundRoots;
397}
398
399int SkDCubic::RootsReal(double A, double B, double C, double D, double s[3]) {
400#ifdef SK_DEBUG
401 // create a string mathematica understands
402 // GDB set print repe 15 # if repeated digits is a bother
403 // set print elements 400 # if line doesn't fit
404 char str[1024];
405 sk_bzero(str, sizeof(str));
406 SK_SNPRINTF(str, sizeof(str), "Solve[%1.19g x^3 + %1.19g x^2 + %1.19g x + %1.19g == 0, x]",
407 A, B, C, D);
408 SkPathOpsDebug::MathematicaIze(str, sizeof(str));
409#if ONE_OFF_DEBUG && ONE_OFF_DEBUG_MATHEMATICA
410 SkDebugf("%s\n", str);
411#endif
412#endif
413 if (approximately_zero(A)
414 && approximately_zero_when_compared_to(A, B)
415 && approximately_zero_when_compared_to(A, C)
416 && approximately_zero_when_compared_to(A, D)) { // we're just a quadratic
417 return SkDQuad::RootsReal(B, C, D, s);
418 }
419 if (approximately_zero_when_compared_to(D, A)
420 && approximately_zero_when_compared_to(D, B)
421 && approximately_zero_when_compared_to(D, C)) { // 0 is one root
422 int num = SkDQuad::RootsReal(A, B, C, s);
423 for (int i = 0; i < num; ++i) {
424 if (approximately_zero(s[i])) {
425 return num;
426 }
427 }
428 s[num++] = 0;
429 return num;
430 }
431 if (approximately_zero(A + B + C + D)) { // 1 is one root
432 int num = SkDQuad::RootsReal(A, A + B, -D, s);
433 for (int i = 0; i < num; ++i) {
434 if (AlmostDequalUlps(s[i], 1)) {
435 return num;
436 }
437 }
438 s[num++] = 1;
439 return num;
440 }
441 double a, b, c;
442 {
443 double invA = 1 / A;
444 a = B * invA;
445 b = C * invA;
446 c = D * invA;
447 }
448 double a2 = a * a;
449 double Q = (a2 - b * 3) / 9;
450 double R = (2 * a2 * a - 9 * a * b + 27 * c) / 54;
451 double R2 = R * R;
452 double Q3 = Q * Q * Q;
453 double R2MinusQ3 = R2 - Q3;
454 double adiv3 = a / 3;
455 double r;
456 double* roots = s;
457 if (R2MinusQ3 < 0) { // we have 3 real roots
458 // the divide/root can, due to finite precisions, be slightly outside of -1...1
459 double theta = acos(SkTPin(R / sqrt(Q3), -1., 1.));
460 double neg2RootQ = -2 * sqrt(Q);
461
462 r = neg2RootQ * cos(theta / 3) - adiv3;
463 *roots++ = r;
464
465 r = neg2RootQ * cos((theta + 2 * PI) / 3) - adiv3;
466 if (!AlmostDequalUlps(s[0], r)) {
467 *roots++ = r;
468 }
469 r = neg2RootQ * cos((theta - 2 * PI) / 3) - adiv3;
470 if (!AlmostDequalUlps(s[0], r) && (roots - s == 1 || !AlmostDequalUlps(s[1], r))) {
471 *roots++ = r;
472 }
473 } else { // we have 1 real root
474 double sqrtR2MinusQ3 = sqrt(R2MinusQ3);
475 double A = fabs(R) + sqrtR2MinusQ3;
476 A = SkDCubeRoot(A);
477 if (R > 0) {
478 A = -A;
479 }
480 if (A != 0) {
481 A += Q / A;
482 }
483 r = A - adiv3;
484 *roots++ = r;
485 if (AlmostDequalUlps((double) R2, (double) Q3)) {
486 r = -A / 2 - adiv3;
487 if (!AlmostDequalUlps(s[0], r)) {
488 *roots++ = r;
489 }
490 }
491 }
492 return static_cast<int>(roots - s);
493}
494
495// OPTIMIZE? compute t^2, t(1-t), and (1-t)^2 and pass them to another version of derivative at t?
496SkDVector SkDCubic::dxdyAtT(double t) const {
497 SkDVector result = { derivative_at_t(&fPts[0].fX, t), derivative_at_t(&fPts[0].fY, t) };
498 if (result.fX == 0 && result.fY == 0) {
499 if (t == 0) {
500 result = fPts[2] - fPts[0];
501 } else if (t == 1) {
502 result = fPts[3] - fPts[1];
503 } else {
504 // incomplete
505 SkDebugf("!c");
506 }
507 if (result.fX == 0 && result.fY == 0 && zero_or_one(t)) {
508 result = fPts[3] - fPts[0];
509 }
510 }
511 return result;
512}
513
514// OPTIMIZE? share code with formulate_F1DotF2
515int SkDCubic::findInflections(double tValues[]) const {
516 double Ax = fPts[1].fX - fPts[0].fX;
517 double Ay = fPts[1].fY - fPts[0].fY;
518 double Bx = fPts[2].fX - 2 * fPts[1].fX + fPts[0].fX;
519 double By = fPts[2].fY - 2 * fPts[1].fY + fPts[0].fY;
520 double Cx = fPts[3].fX + 3 * (fPts[1].fX - fPts[2].fX) - fPts[0].fX;
521 double Cy = fPts[3].fY + 3 * (fPts[1].fY - fPts[2].fY) - fPts[0].fY;
522 return SkDQuad::RootsValidT(Bx * Cy - By * Cx, Ax * Cy - Ay * Cx, Ax * By - Ay * Bx, tValues);
523}
524
525static void formulate_F1DotF2(const double src[], double coeff[4]) {
526 double a = src[2] - src[0];
527 double b = src[4] - 2 * src[2] + src[0];
528 double c = src[6] + 3 * (src[2] - src[4]) - src[0];
529 coeff[0] = c * c;
530 coeff[1] = 3 * b * c;
531 coeff[2] = 2 * b * b + c * a;
532 coeff[3] = a * b;
533}
534
535/** SkDCubic'(t) = At^2 + Bt + C, where
536 A = 3(-a + 3(b - c) + d)
537 B = 6(a - 2b + c)
538 C = 3(b - a)
539 Solve for t, keeping only those that fit between 0 < t < 1
540*/
541int SkDCubic::FindExtrema(const double src[], double tValues[2]) {
542 // we divide A,B,C by 3 to simplify
543 double a = src[0];
544 double b = src[2];
545 double c = src[4];
546 double d = src[6];
547 double A = d - a + 3 * (b - c);
548 double B = 2 * (a - b - b + c);
549 double C = b - a;
550
551 return SkDQuad::RootsValidT(A, B, C, tValues);
552}
553
554/* from SkGeometry.cpp
555 Looking for F' dot F'' == 0
556
557 A = b - a
558 B = c - 2b + a
559 C = d - 3c + 3b - a
560
561 F' = 3Ct^2 + 6Bt + 3A
562 F'' = 6Ct + 6B
563
564 F' dot F'' -> CCt^3 + 3BCt^2 + (2BB + CA)t + AB
565*/
566int SkDCubic::findMaxCurvature(double tValues[]) const {
567 double coeffX[4], coeffY[4];
568 int i;
569 formulate_F1DotF2(&fPts[0].fX, coeffX);
570 formulate_F1DotF2(&fPts[0].fY, coeffY);
571 for (i = 0; i < 4; i++) {
572 coeffX[i] = coeffX[i] + coeffY[i];
573 }
574 return RootsValidT(coeffX[0], coeffX[1], coeffX[2], coeffX[3], tValues);
575}
576
577SkDPoint SkDCubic::ptAtT(double t) const {
578 if (0 == t) {
579 return fPts[0];
580 }
581 if (1 == t) {
582 return fPts[3];
583 }
584 double one_t = 1 - t;
585 double one_t2 = one_t * one_t;
586 double a = one_t2 * one_t;
587 double b = 3 * one_t2 * t;
588 double t2 = t * t;
589 double c = 3 * one_t * t2;
590 double d = t2 * t;
591 SkDPoint result = {a * fPts[0].fX + b * fPts[1].fX + c * fPts[2].fX + d * fPts[3].fX,
592 a * fPts[0].fY + b * fPts[1].fY + c * fPts[2].fY + d * fPts[3].fY};
593 return result;
594}
595
596/*
597 Given a cubic c, t1, and t2, find a small cubic segment.
598
599 The new cubic is defined as points A, B, C, and D, where
600 s1 = 1 - t1
601 s2 = 1 - t2
602 A = c[0]*s1*s1*s1 + 3*c[1]*s1*s1*t1 + 3*c[2]*s1*t1*t1 + c[3]*t1*t1*t1
603 D = c[0]*s2*s2*s2 + 3*c[1]*s2*s2*t2 + 3*c[2]*s2*t2*t2 + c[3]*t2*t2*t2
604
605 We don't have B or C. So We define two equations to isolate them.
606 First, compute two reference T values 1/3 and 2/3 from t1 to t2:
607
608 c(at (2*t1 + t2)/3) == E
609 c(at (t1 + 2*t2)/3) == F
610
611 Next, compute where those values must be if we know the values of B and C:
612
613 _12 = A*2/3 + B*1/3
614 12_ = A*1/3 + B*2/3
615 _23 = B*2/3 + C*1/3
616 23_ = B*1/3 + C*2/3
617 _34 = C*2/3 + D*1/3
618 34_ = C*1/3 + D*2/3
619 _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
620 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
621 _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
622 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
623 _1234 = (A*4/9 + B*4/9 + C*1/9)*2/3 + (B*4/9 + C*4/9 + D*1/9)*1/3
624 = A*8/27 + B*12/27 + C*6/27 + D*1/27
625 = E
626 1234_ = (A*1/9 + B*4/9 + C*4/9)*1/3 + (B*1/9 + C*4/9 + D*4/9)*2/3
627 = A*1/27 + B*6/27 + C*12/27 + D*8/27
628 = F
629 E*27 = A*8 + B*12 + C*6 + D
630 F*27 = A + B*6 + C*12 + D*8
631
632Group the known values on one side:
633
634 M = E*27 - A*8 - D = B*12 + C* 6
635 N = F*27 - A - D*8 = B* 6 + C*12
636 M*2 - N = B*18
637 N*2 - M = C*18
638 B = (M*2 - N)/18
639 C = (N*2 - M)/18
640 */
641
642static double interp_cubic_coords(const double* src, double t) {
643 double ab = SkDInterp(src[0], src[2], t);
644 double bc = SkDInterp(src[2], src[4], t);
645 double cd = SkDInterp(src[4], src[6], t);
646 double abc = SkDInterp(ab, bc, t);
647 double bcd = SkDInterp(bc, cd, t);
648 double abcd = SkDInterp(abc, bcd, t);
649 return abcd;
650}
651
652SkDCubic SkDCubic::subDivide(double t1, double t2) const {
653 if (t1 == 0 || t2 == 1) {
654 if (t1 == 0 && t2 == 1) {
655 return *this;
656 }
657 SkDCubicPair pair = chopAt(t1 == 0 ? t2 : t1);
658 SkDCubic dst = t1 == 0 ? pair.first() : pair.second();
659 return dst;
660 }
661 SkDCubic dst;
662 double ax = dst[0].fX = interp_cubic_coords(&fPts[0].fX, t1);
663 double ay = dst[0].fY = interp_cubic_coords(&fPts[0].fY, t1);
664 double ex = interp_cubic_coords(&fPts[0].fX, (t1*2+t2)/3);
665 double ey = interp_cubic_coords(&fPts[0].fY, (t1*2+t2)/3);
666 double fx = interp_cubic_coords(&fPts[0].fX, (t1+t2*2)/3);
667 double fy = interp_cubic_coords(&fPts[0].fY, (t1+t2*2)/3);
668 double dx = dst[3].fX = interp_cubic_coords(&fPts[0].fX, t2);
669 double dy = dst[3].fY = interp_cubic_coords(&fPts[0].fY, t2);
670 double mx = ex * 27 - ax * 8 - dx;
671 double my = ey * 27 - ay * 8 - dy;
672 double nx = fx * 27 - ax - dx * 8;
673 double ny = fy * 27 - ay - dy * 8;
674 /* bx = */ dst[1].fX = (mx * 2 - nx) / 18;
675 /* by = */ dst[1].fY = (my * 2 - ny) / 18;
676 /* cx = */ dst[2].fX = (nx * 2 - mx) / 18;
677 /* cy = */ dst[2].fY = (ny * 2 - my) / 18;
678 // FIXME: call align() ?
679 return dst;
680}
681
682void SkDCubic::subDivide(const SkDPoint& a, const SkDPoint& d,
683 double t1, double t2, SkDPoint dst[2]) const {
684 SkASSERT(t1 != t2);
685 // this approach assumes that the control points computed directly are accurate enough
686 SkDCubic sub = subDivide(t1, t2);
687 dst[0] = sub[1] + (a - sub[0]);
688 dst[1] = sub[2] + (d - sub[3]);
689 if (t1 == 0 || t2 == 0) {
690 align(0, 1, t1 == 0 ? &dst[0] : &dst[1]);
691 }
692 if (t1 == 1 || t2 == 1) {
693 align(3, 2, t1 == 1 ? &dst[0] : &dst[1]);
694 }
695 if (AlmostBequalUlps(dst[0].fX, a.fX)) {
696 dst[0].fX = a.fX;
697 }
698 if (AlmostBequalUlps(dst[0].fY, a.fY)) {
699 dst[0].fY = a.fY;
700 }
701 if (AlmostBequalUlps(dst[1].fX, d.fX)) {
702 dst[1].fX = d.fX;
703 }
704 if (AlmostBequalUlps(dst[1].fY, d.fY)) {
705 dst[1].fY = d.fY;
706 }
707}
708
709bool SkDCubic::toFloatPoints(SkPoint* pts) const {
710 const double* dCubic = &fPts[0].fX;
711 SkScalar* cubic = &pts[0].fX;
712 for (int index = 0; index < kPointCount * 2; ++index) {
713 cubic[index] = SkDoubleToScalar(dCubic[index]);
714 if (SkScalarAbs(cubic[index]) < FLT_EPSILON_ORDERABLE_ERR) {
715 cubic[index] = 0;
716 }
717 }
718 return SkScalarsAreFinite(&pts->fX, kPointCount * 2);
719}
720
721double SkDCubic::top(const SkDCubic& dCurve, double startT, double endT, SkDPoint*topPt) const {
722 double extremeTs[2];
723 double topT = -1;
724 int roots = SkDCubic::FindExtrema(&fPts[0].fY, extremeTs);
725 for (int index = 0; index < roots; ++index) {
726 double t = startT + (endT - startT) * extremeTs[index];
727 SkDPoint mid = dCurve.ptAtT(t);
728 if (topPt->fY > mid.fY || (topPt->fY == mid.fY && topPt->fX > mid.fX)) {
729 topT = t;
730 *topPt = mid;
731 }
732 }
733 return topT;
734}
735
736int SkTCubic::intersectRay(SkIntersections* i, const SkDLine& line) const {
737 return i->intersectRay(fCubic, line);
738}
739
740bool SkTCubic::hullIntersects(const SkDQuad& quad, bool* isLinear) const {
741 return quad.hullIntersects(fCubic, isLinear);
742}
743
744bool SkTCubic::hullIntersects(const SkDConic& conic, bool* isLinear) const {
745 return conic.hullIntersects(fCubic, isLinear);
746}
747
748void SkTCubic::setBounds(SkDRect* rect) const {
749 rect->setBounds(fCubic);
750}
751