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