1/*
2 * Copyright 2013 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 "PathOpsTestCommon.h"
8#include "SkIntersections.h"
9#include "SkOpSegment.h"
10#include "SkPathOpsTriangle.h"
11#include "SkRandom.h"
12#include "SkTArray.h"
13#include "SkTSort.h"
14#include "Test.h"
15
16static bool gPathOpsAngleIdeasVerbose = false;
17static bool gPathOpsAngleIdeasEnableBruteCheck = false;
18
19class PathOpsAngleTester {
20public:
21    static int ConvexHullOverlaps(const SkOpAngle& lh, const SkOpAngle& rh) {
22        return lh.convexHullOverlaps(rh);
23    }
24
25    static int EndsIntersect(const SkOpAngle& lh, const SkOpAngle& rh) {
26        return lh.endsIntersect(rh);
27    }
28};
29
30struct TRange {
31    double tMin1;
32    double tMin2;
33    double t1;
34    double t2;
35    double tMin;
36    double a1;
37    double a2;
38    bool ccw;
39};
40
41static double testArc(skiatest::Reporter* reporter, const SkDQuad& quad, const SkDQuad& arcRef,
42        int octant) {
43    SkDQuad arc = arcRef;
44    SkDVector offset = {quad[0].fX, quad[0].fY};
45    arc[0] += offset;
46    arc[1] += offset;
47    arc[2] += offset;
48    SkIntersections i;
49    i.intersect(arc, quad);
50    if (i.used() == 0) {
51        return -1;
52    }
53    int smallest = -1;
54    double t = 2;
55    for (int idx = 0; idx < i.used(); ++idx) {
56        if (i[0][idx] > 1 || i[0][idx] < 0) {
57            i.reset();
58            i.intersect(arc, quad);
59        }
60        if (i[1][idx] > 1 || i[1][idx] < 0) {
61            i.reset();
62            i.intersect(arc, quad);
63        }
64        if (t > i[1][idx]) {
65            smallest = idx;
66            t = i[1][idx];
67        }
68    }
69    REPORTER_ASSERT(reporter, smallest >= 0);
70    REPORTER_ASSERT(reporter, t >= 0 && t <= 1);
71    return i[1][smallest];
72}
73
74static void orderQuads(skiatest::Reporter* reporter, const SkDQuad& quad, double radius,
75        SkTArray<double, false>* tArray) {
76    double r = radius;
77    double s = r * SK_ScalarTanPIOver8;
78    double m = r * SK_ScalarRoot2Over2;
79    // construct circle from quads
80    const SkDQuad circle[8] = {{{{ r,  0}, { r, -s}, { m, -m}}},
81                                {{{ m, -m}, { s, -r}, { 0, -r}}},
82                                {{{ 0, -r}, {-s, -r}, {-m, -m}}},
83                                {{{-m, -m}, {-r, -s}, {-r,  0}}},
84                                {{{-r,  0}, {-r,  s}, {-m,  m}}},
85                                {{{-m,  m}, {-s,  r}, { 0,  r}}},
86                                {{{ 0,  r}, { s,  r}, { m,  m}}},
87                                {{{ m,  m}, { r,  s}, { r,  0}}}};
88    for (int octant = 0; octant < 8; ++octant) {
89        double t = testArc(reporter, quad, circle[octant], octant);
90        if (t < 0) {
91            continue;
92        }
93        for (int index = 0; index < tArray->count(); ++index) {
94            double matchT = (*tArray)[index];
95            if (approximately_equal(t, matchT)) {
96                goto next;
97            }
98        }
99        tArray->push_back(t);
100next:   ;
101    }
102}
103
104static double quadAngle(skiatest::Reporter* reporter, const SkDQuad& quad, double t) {
105    const SkDVector& pt = quad.ptAtT(t) - quad[0];
106    double angle = (atan2(pt.fY, pt.fX) + SK_ScalarPI) * 8 / (SK_ScalarPI * 2);
107    REPORTER_ASSERT(reporter, angle >= 0 && angle <= 8);
108    return angle;
109}
110
111static bool angleDirection(double a1, double a2) {
112    double delta = a1 - a2;
113    return (delta < 4 && delta > 0) || delta < -4;
114}
115
116static void setQuadHullSweep(const SkDQuad& quad, SkDVector sweep[2]) {
117    sweep[0] = quad[1] - quad[0];
118    sweep[1] = quad[2] - quad[0];
119}
120
121static double distEndRatio(double dist, const SkDQuad& quad) {
122    SkDVector v[] = {quad[2] - quad[0], quad[1] - quad[0], quad[2] - quad[1]};
123    double longest = SkTMax(v[0].length(), SkTMax(v[1].length(), v[2].length()));
124    return longest / dist;
125}
126
127static bool checkParallel(skiatest::Reporter* reporter, const SkDQuad& quad1, const SkDQuad& quad2) {
128    SkDVector sweep[2], tweep[2];
129    setQuadHullSweep(quad1, sweep);
130    setQuadHullSweep(quad2, tweep);
131    // if the ctrl tangents are not nearly parallel, use them
132    // solve for opposite direction displacement scale factor == m
133    // initial dir = v1.cross(v2) == v2.x * v1.y - v2.y * v1.x
134    // displacement of q1[1] : dq1 = { -m * v1.y, m * v1.x } + q1[1]
135    // straight angle when : v2.x * (dq1.y - q1[0].y) == v2.y * (dq1.x - q1[0].x)
136    //                       v2.x * (m * v1.x + v1.y) == v2.y * (-m * v1.y + v1.x)
137    // - m * (v2.x * v1.x + v2.y * v1.y) == v2.x * v1.y - v2.y * v1.x
138    // m = (v2.y * v1.x - v2.x * v1.y) / (v2.x * v1.x + v2.y * v1.y)
139    // m = v1.cross(v2) / v1.dot(v2)
140    double s0dt0 = sweep[0].dot(tweep[0]);
141    REPORTER_ASSERT(reporter, s0dt0 != 0);
142    double s0xt0 = sweep[0].crossCheck(tweep[0]);
143    double m = s0xt0 / s0dt0;
144    double sDist = sweep[0].length() * m;
145    double tDist = tweep[0].length() * m;
146    bool useS = fabs(sDist) < fabs(tDist);
147    double mFactor = fabs(useS ? distEndRatio(sDist, quad1) : distEndRatio(tDist, quad2));
148    if (mFactor < 5000) {  // empirically found limit
149        return s0xt0 < 0;
150    }
151    SkDVector m0 = quad1.ptAtT(0.5) - quad1[0];
152    SkDVector m1 = quad2.ptAtT(0.5) - quad2[0];
153    return m0.crossCheck(m1) < 0;
154}
155
156/* returns
157   -1 if overlaps
158    0 if no overlap cw
159    1 if no overlap ccw
160*/
161static int quadHullsOverlap(skiatest::Reporter* reporter, const SkDQuad& quad1,
162        const SkDQuad& quad2) {
163    SkDVector sweep[2], tweep[2];
164    setQuadHullSweep(quad1, sweep);
165    setQuadHullSweep(quad2, tweep);
166    double s0xs1 = sweep[0].crossCheck(sweep[1]);
167    double s0xt0 = sweep[0].crossCheck(tweep[0]);
168    double s1xt0 = sweep[1].crossCheck(tweep[0]);
169    bool tBetweenS = s0xs1 > 0 ? s0xt0 > 0 && s1xt0 < 0 : s0xt0 < 0 && s1xt0 > 0;
170    double s0xt1 = sweep[0].crossCheck(tweep[1]);
171    double s1xt1 = sweep[1].crossCheck(tweep[1]);
172    tBetweenS |= s0xs1 > 0 ? s0xt1 > 0 && s1xt1 < 0 : s0xt1 < 0 && s1xt1 > 0;
173    double t0xt1 = tweep[0].crossCheck(tweep[1]);
174    if (tBetweenS) {
175        return -1;
176    }
177    if ((s0xt0 == 0 && s1xt1 == 0) || (s1xt0 == 0 && s0xt1 == 0)) {  // s0 to s1 equals t0 to t1
178        return -1;
179    }
180    bool sBetweenT = t0xt1 > 0 ? s0xt0 < 0 && s0xt1 > 0 : s0xt0 > 0 && s0xt1 < 0;
181    sBetweenT |= t0xt1 > 0 ? s1xt0 < 0 && s1xt1 > 0 : s1xt0 > 0 && s1xt1 < 0;
182    if (sBetweenT) {
183        return -1;
184    }
185    // if all of the sweeps are in the same half plane, then the order of any pair is enough
186    if (s0xt0 >= 0 && s0xt1 >= 0 && s1xt0 >= 0 && s1xt1 >= 0) {
187        return 0;
188    }
189    if (s0xt0 <= 0 && s0xt1 <= 0 && s1xt0 <= 0 && s1xt1 <= 0) {
190        return 1;
191    }
192    // if the outside sweeps are greater than 180 degress:
193        // first assume the inital tangents are the ordering
194        // if the midpoint direction matches the inital order, that is enough
195    SkDVector m0 = quad1.ptAtT(0.5) - quad1[0];
196    SkDVector m1 = quad2.ptAtT(0.5) - quad2[0];
197    double m0xm1 = m0.crossCheck(m1);
198    if (s0xt0 > 0 && m0xm1 > 0) {
199        return 0;
200    }
201    if (s0xt0 < 0 && m0xm1 < 0) {
202        return 1;
203    }
204    REPORTER_ASSERT(reporter, s0xt0 != 0);
205    return checkParallel(reporter, quad1, quad2);
206}
207
208static double radianSweep(double start, double end) {
209    double sweep = end - start;
210    if (sweep > SK_ScalarPI) {
211        sweep -= 2 * SK_ScalarPI;
212    } else if (sweep < -SK_ScalarPI) {
213        sweep += 2 * SK_ScalarPI;
214    }
215    return sweep;
216}
217
218static bool radianBetween(double start, double test, double end) {
219    double startToEnd = radianSweep(start, end);
220    double startToTest = radianSweep(start, test);
221    double testToEnd = radianSweep(test, end);
222    return (startToTest <= 0 && testToEnd <= 0 && startToTest >= startToEnd) ||
223        (startToTest >= 0 && testToEnd >= 0 && startToTest <= startToEnd);
224}
225
226static bool orderTRange(skiatest::Reporter* reporter, const SkDQuad& quad1, const SkDQuad& quad2,
227        double r, TRange* result) {
228    SkTArray<double, false> t1Array, t2Array;
229    orderQuads(reporter, quad1, r, &t1Array);
230    orderQuads(reporter,quad2, r, &t2Array);
231    if (!t1Array.count() || !t2Array.count()) {
232        return false;
233    }
234    SkTQSort<double>(t1Array.begin(), t1Array.end() - 1);
235    SkTQSort<double>(t2Array.begin(), t2Array.end() - 1);
236    double t1 = result->tMin1 = t1Array[0];
237    double t2 = result->tMin2 = t2Array[0];
238    double a1 = quadAngle(reporter,quad1, t1);
239    double a2 = quadAngle(reporter,quad2, t2);
240    if (approximately_equal(a1, a2)) {
241        return false;
242    }
243    bool refCCW = angleDirection(a1, a2);
244    result->t1 = t1;
245    result->t2 = t2;
246    result->tMin = SkTMin(t1, t2);
247    result->a1 = a1;
248    result->a2 = a2;
249    result->ccw = refCCW;
250    return true;
251}
252
253static bool equalPoints(const SkDPoint& pt1, const SkDPoint& pt2, double max) {
254    return approximately_zero_when_compared_to(pt1.fX - pt2.fX, max)
255            && approximately_zero_when_compared_to(pt1.fY - pt2.fY, max);
256}
257
258static double maxDist(const SkDQuad& quad) {
259    SkDRect bounds;
260    bounds.setBounds(quad);
261    SkDVector corner[4] = {
262        { bounds.fLeft - quad[0].fX, bounds.fTop - quad[0].fY },
263        { bounds.fRight - quad[0].fX, bounds.fTop - quad[0].fY },
264        { bounds.fLeft - quad[0].fX, bounds.fBottom - quad[0].fY },
265        { bounds.fRight - quad[0].fX, bounds.fBottom - quad[0].fY }
266    };
267    double max = 0;
268    for (unsigned index = 0; index < SK_ARRAY_COUNT(corner); ++index) {
269        max = SkTMax(max, corner[index].length());
270    }
271    return max;
272}
273
274static double maxQuad(const SkDQuad& quad) {
275    double max = 0;
276    for (int index = 0; index < 2; ++index) {
277        max = SkTMax(max, fabs(quad[index].fX));
278        max = SkTMax(max, fabs(quad[index].fY));
279    }
280    return max;
281}
282
283static bool bruteMinT(skiatest::Reporter* reporter, const SkDQuad& quad1, const SkDQuad& quad2,
284        TRange* lowerRange, TRange* upperRange) {
285    double maxRadius = SkTMin(maxDist(quad1), maxDist(quad2));
286    double maxQuads = SkTMax(maxQuad(quad1), maxQuad(quad2));
287    double r = maxRadius / 2;
288    double rStep = r / 2;
289    SkDPoint best1 = {SK_ScalarInfinity, SK_ScalarInfinity};
290    SkDPoint best2 = {SK_ScalarInfinity, SK_ScalarInfinity};
291    int bestCCW = -1;
292    double bestR = maxRadius;
293    upperRange->tMin = 0;
294    lowerRange->tMin = 1;
295    do {
296        do {  // find upper bounds of single result
297            TRange tRange;
298            bool stepUp = orderTRange(reporter, quad1, quad2, r, &tRange);
299            if (stepUp) {
300                SkDPoint pt1 = quad1.ptAtT(tRange.t1);
301                if (equalPoints(pt1, best1, maxQuads)) {
302                    break;
303                }
304                best1 = pt1;
305                SkDPoint pt2 = quad2.ptAtT(tRange.t2);
306                if (equalPoints(pt2, best2, maxQuads)) {
307                    break;
308                }
309                best2 = pt2;
310                if (gPathOpsAngleIdeasVerbose) {
311                    SkDebugf("u bestCCW=%d ccw=%d bestMin=%1.9g:%1.9g r=%1.9g tMin=%1.9g\n",
312                            bestCCW, tRange.ccw, lowerRange->tMin, upperRange->tMin, r,
313                            tRange.tMin);
314                }
315                if (bestCCW >= 0 && bestCCW != (int) tRange.ccw) {
316                    if (tRange.tMin < upperRange->tMin) {
317                        upperRange->tMin = 0;
318                    } else {
319                        stepUp = false;
320                    }
321                }
322                if (upperRange->tMin < tRange.tMin) {
323                    bestCCW = tRange.ccw;
324                    bestR = r;
325                    *upperRange = tRange;
326                }
327                if (lowerRange->tMin > tRange.tMin) {
328                    *lowerRange = tRange;
329                }
330            }
331            r += stepUp ? rStep : -rStep;
332            rStep /= 2;
333        } while (rStep > FLT_EPSILON);
334        if (bestCCW < 0) {
335            REPORTER_ASSERT(reporter, bestR < maxRadius);
336            return false;
337        }
338        double lastHighR = bestR;
339        r = bestR / 2;
340        rStep = r / 2;
341        do {  // find lower bounds of single result
342            TRange tRange;
343            bool success = orderTRange(reporter, quad1, quad2, r, &tRange);
344            if (success) {
345                if (gPathOpsAngleIdeasVerbose) {
346                    SkDebugf("l bestCCW=%d ccw=%d bestMin=%1.9g:%1.9g r=%1.9g tMin=%1.9g\n",
347                            bestCCW, tRange.ccw, lowerRange->tMin, upperRange->tMin, r,
348                            tRange.tMin);
349                }
350                if (bestCCW != (int) tRange.ccw || upperRange->tMin < tRange.tMin) {
351                    bestCCW = tRange.ccw;
352                    *upperRange = tRange;
353                    bestR = lastHighR;
354                    break;  // need to establish a new upper bounds
355                }
356                SkDPoint pt1 = quad1.ptAtT(tRange.t1);
357                SkDPoint pt2 = quad2.ptAtT(tRange.t2);
358                if (equalPoints(pt1, best1, maxQuads)) {
359                    goto breakOut;
360                }
361                best1 = pt1;
362                if (equalPoints(pt2, best2, maxQuads)) {
363                    goto breakOut;
364                }
365                best2 = pt2;
366                if (equalPoints(pt1, pt2, maxQuads)) {
367                    success = false;
368                } else {
369                    if (upperRange->tMin < tRange.tMin) {
370                        *upperRange = tRange;
371                    }
372                    if (lowerRange->tMin > tRange.tMin) {
373                        *lowerRange = tRange;
374                    }
375                }
376                lastHighR = SkTMin(r, lastHighR);
377            }
378            r += success ? -rStep : rStep;
379            rStep /= 2;
380        } while (rStep > FLT_EPSILON);
381    } while (rStep > FLT_EPSILON);
382breakOut:
383    if (gPathOpsAngleIdeasVerbose) {
384        SkDebugf("l a2-a1==%1.9g\n", lowerRange->a2 - lowerRange->a1);
385    }
386    return true;
387}
388
389static void bruteForce(skiatest::Reporter* reporter, const SkDQuad& quad1, const SkDQuad& quad2,
390        bool ccw) {
391    if (!gPathOpsAngleIdeasEnableBruteCheck) {
392        return;
393    }
394    TRange lowerRange, upperRange;
395    bool result = bruteMinT(reporter, quad1, quad2, &lowerRange, &upperRange);
396    REPORTER_ASSERT(reporter, result);
397    double angle = fabs(lowerRange.a2 - lowerRange.a1);
398    REPORTER_ASSERT(reporter, angle > 3.998 || ccw == upperRange.ccw);
399}
400
401static bool bruteForceCheck(skiatest::Reporter* reporter, const SkDQuad& quad1,
402        const SkDQuad& quad2, bool ccw) {
403    TRange lowerRange, upperRange;
404    bool result = bruteMinT(reporter, quad1, quad2, &lowerRange, &upperRange);
405    REPORTER_ASSERT(reporter, result);
406    return ccw == upperRange.ccw;
407}
408
409class PathOpsSegmentTester {
410public:
411    static void ConstructQuad(SkOpSegment* segment, SkPoint shortQuad[3]) {
412        segment->debugConstructQuad(shortQuad);
413    }
414};
415
416static void makeSegment(const SkDQuad& quad, SkPoint shortQuad[3], SkOpSegment* result) {
417    shortQuad[0] = quad[0].asSkPoint();
418    shortQuad[1] = quad[1].asSkPoint();
419    shortQuad[2] = quad[2].asSkPoint();
420    PathOpsSegmentTester::ConstructQuad(result, shortQuad);
421}
422
423static void testQuadAngles(skiatest::Reporter* reporter, const SkDQuad& quad1, const SkDQuad& quad2,
424        int testNo) {
425    SkPoint shortQuads[2][3];
426    SkOpSegment seg[2];
427    makeSegment(quad1, shortQuads[0], &seg[0]);
428    makeSegment(quad2, shortQuads[1], &seg[1]);
429    int realOverlap = PathOpsAngleTester::ConvexHullOverlaps(*seg[0].debugLastAngle(),
430            *seg[1].debugLastAngle());
431    const SkDPoint& origin = quad1[0];
432    REPORTER_ASSERT(reporter, origin == quad2[0]);
433    double a1s = atan2(origin.fY - quad1[1].fY, quad1[1].fX - origin.fX);
434    double a1e = atan2(origin.fY - quad1[2].fY, quad1[2].fX - origin.fX);
435    double a2s = atan2(origin.fY - quad2[1].fY, quad2[1].fX - origin.fX);
436    double a2e = atan2(origin.fY - quad2[2].fY, quad2[2].fX - origin.fX);
437    bool oldSchoolOverlap = radianBetween(a1s, a2s, a1e)
438        || radianBetween(a1s, a2e, a1e) || radianBetween(a2s, a1s, a2e)
439        || radianBetween(a2s, a1e, a2e);
440    int overlap = quadHullsOverlap(reporter, quad1, quad2);
441    bool realMatchesOverlap = realOverlap == overlap || SK_ScalarPI - fabs(a2s - a1s) < 0.002;
442    if (realOverlap != overlap) {
443        SkDebugf("\nSK_ScalarPI - fabs(a2s - a1s) = %1.9g\n", SK_ScalarPI - fabs(a2s - a1s));
444    }
445    if (!realMatchesOverlap) {
446        DumpQ(quad1, quad2, testNo);
447    }
448    REPORTER_ASSERT(reporter, realMatchesOverlap);
449    if (oldSchoolOverlap != (overlap < 0)) {
450        overlap = quadHullsOverlap(reporter, quad1, quad2);  // set a breakpoint and debug if assert fires
451        REPORTER_ASSERT(reporter, oldSchoolOverlap == (overlap < 0));
452    }
453    SkDVector v1s = quad1[1] - quad1[0];
454    SkDVector v1e = quad1[2] - quad1[0];
455    SkDVector v2s = quad2[1] - quad2[0];
456    SkDVector v2e = quad2[2] - quad2[0];
457    double vDir[2] = { v1s.cross(v1e), v2s.cross(v2e) };
458    bool ray1In2 = v1s.cross(v2s) * vDir[1] <= 0 && v1s.cross(v2e) * vDir[1] >= 0;
459    bool ray2In1 = v2s.cross(v1s) * vDir[0] <= 0 && v2s.cross(v1e) * vDir[0] >= 0;
460    if (overlap >= 0) {
461        // verify that hulls really don't overlap
462        REPORTER_ASSERT(reporter, !ray1In2);
463        REPORTER_ASSERT(reporter, !ray2In1);
464        bool ctrl1In2 = v1e.cross(v2s) * vDir[1] <= 0 && v1e.cross(v2e) * vDir[1] >= 0;
465        REPORTER_ASSERT(reporter, !ctrl1In2);
466        bool ctrl2In1 = v2e.cross(v1s) * vDir[0] <= 0 && v2e.cross(v1e) * vDir[0] >= 0;
467        REPORTER_ASSERT(reporter, !ctrl2In1);
468        // check answer against reference
469        bruteForce(reporter, quad1, quad2, overlap > 0);
470    }
471    // continue end point rays and see if they intersect the opposite curve
472    SkDLine rays[] = {{{origin, quad2[2]}}, {{origin, quad1[2]}}};
473    const SkDQuad* quads[] = {&quad1, &quad2};
474    SkDVector midSpokes[2];
475    SkIntersections intersect[2];
476    double minX, minY, maxX, maxY;
477    minX = minY = SK_ScalarInfinity;
478    maxX = maxY = -SK_ScalarInfinity;
479    double maxWidth = 0;
480    bool useIntersect = false;
481    double smallestTs[] = {1, 1};
482    for (unsigned index = 0; index < SK_ARRAY_COUNT(quads); ++index) {
483        const SkDQuad& q = *quads[index];
484        midSpokes[index] = q.ptAtT(0.5) - origin;
485        minX = SkTMin(SkTMin(SkTMin(minX, origin.fX), q[1].fX), q[2].fX);
486        minY = SkTMin(SkTMin(SkTMin(minY, origin.fY), q[1].fY), q[2].fY);
487        maxX = SkTMax(SkTMax(SkTMax(maxX, origin.fX), q[1].fX), q[2].fX);
488        maxY = SkTMax(SkTMax(SkTMax(maxY, origin.fY), q[1].fY), q[2].fY);
489        maxWidth = SkTMax(maxWidth, SkTMax(maxX - minX, maxY - minY));
490        intersect[index].intersectRay(q, rays[index]);
491        const SkIntersections& i = intersect[index];
492        REPORTER_ASSERT(reporter, i.used() >= 1);
493        bool foundZero = false;
494        double smallT = 1;
495        for (int idx2 = 0; idx2 < i.used(); ++idx2) {
496            double t = i[0][idx2];
497            if (t == 0) {
498                foundZero = true;
499                continue;
500            }
501            if (smallT > t) {
502                smallT = t;
503            }
504        }
505        REPORTER_ASSERT(reporter, foundZero == true);
506        if (smallT == 1) {
507            continue;
508        }
509        SkDVector ray = q.ptAtT(smallT) - origin;
510        SkDVector end = rays[index][1] - origin;
511        if (ray.fX * end.fX < 0 || ray.fY * end.fY < 0) {
512            continue;
513        }
514        double rayDist = ray.length();
515        double endDist = end.length();
516        double delta = fabs(rayDist - endDist) / maxWidth;
517        if (delta > 1e-4) {
518            useIntersect ^= true;
519        }
520        smallestTs[index] = smallT;
521    }
522    bool firstInside;
523    if (useIntersect) {
524        int sIndex = (int) (smallestTs[1] < 1);
525        REPORTER_ASSERT(reporter, smallestTs[sIndex ^ 1] == 1);
526        double t = smallestTs[sIndex];
527        const SkDQuad& q = *quads[sIndex];
528        SkDVector ray = q.ptAtT(t) - origin;
529        SkDVector end = rays[sIndex][1] - origin;
530        double rayDist = ray.length();
531        double endDist = end.length();
532        SkDVector mid = q.ptAtT(t / 2) - origin;
533        double midXray = mid.crossCheck(ray);
534        if (gPathOpsAngleIdeasVerbose) {
535            SkDebugf("rayDist>endDist:%d sIndex==0:%d vDir[sIndex]<0:%d midXray<0:%d\n",
536                    rayDist > endDist, sIndex == 0, vDir[sIndex] < 0, midXray < 0);
537        }
538        SkASSERT(SkScalarSignAsInt(SkDoubleToScalar(midXray))
539            == SkScalarSignAsInt(SkDoubleToScalar(vDir[sIndex])));
540        firstInside = (rayDist > endDist) ^ (sIndex == 0) ^ (vDir[sIndex] < 0);
541    } else if (overlap >= 0) {
542        return;  // answer has already been determined
543    } else {
544        firstInside = checkParallel(reporter, quad1, quad2);
545    }
546    if (overlap < 0) {
547        SkDEBUGCODE(int realEnds =)
548                PathOpsAngleTester::EndsIntersect(*seg[0].debugLastAngle(),
549                *seg[1].debugLastAngle());
550        SkASSERT(realEnds == (firstInside ? 1 : 0));
551    }
552    bruteForce(reporter, quad1, quad2, firstInside);
553}
554
555DEF_TEST(PathOpsAngleOverlapHullsOne, reporter) {
556//    gPathOpsAngleIdeasVerbose = true;
557    const SkDQuad quads[] = {
558{{{939.4808349609375, 914.355224609375}, {-357.7921142578125, 590.842529296875}, {736.8936767578125, -350.717529296875}}},
559{{{939.4808349609375, 914.355224609375}, {-182.85418701171875, 634.4552001953125}, {-509.62615966796875, 576.1182861328125}}}
560    };
561    for (int index = 0; index < (int) SK_ARRAY_COUNT(quads); index += 2) {
562        testQuadAngles(reporter, quads[index], quads[index + 1], 0);
563    }
564}
565
566DEF_TEST(PathOpsAngleOverlapHulls, reporter) {
567    if (!gPathOpsAngleIdeasVerbose) {  // takes a while to run -- so exclude it by default
568        return;
569    }
570    SkRandom ran;
571    for (int index = 0; index < 100000; ++index) {
572        if (index % 1000 == 999) SkDebugf(".");
573        SkDPoint origin = {ran.nextRangeF(-1000, 1000), ran.nextRangeF(-1000, 1000)};
574        SkDQuad quad1 = {{origin, {ran.nextRangeF(-1000, 1000), ran.nextRangeF(-1000, 1000)},
575            {ran.nextRangeF(-1000, 1000), ran.nextRangeF(-1000, 1000)}}};
576        if (quad1[0] == quad1[2]) {
577            continue;
578        }
579        SkDQuad quad2 = {{origin, {ran.nextRangeF(-1000, 1000), ran.nextRangeF(-1000, 1000)},
580            {ran.nextRangeF(-1000, 1000), ran.nextRangeF(-1000, 1000)}}};
581        if (quad2[0] == quad2[2]) {
582            continue;
583        }
584        SkIntersections i;
585        i.intersect(quad1, quad2);
586        REPORTER_ASSERT(reporter, i.used() >= 1);
587        if (i.used() > 1) {
588            continue;
589        }
590        testQuadAngles(reporter, quad1, quad2, index);
591    }
592}
593
594DEF_TEST(PathOpsAngleBruteT, reporter) {
595    if (!gPathOpsAngleIdeasVerbose) {  // takes a while to run -- so exclude it by default
596        return;
597    }
598    SkRandom ran;
599    double smaller = SK_Scalar1;
600    SkDQuad small[2];
601    SkDEBUGCODE(int smallIndex);
602    for (int index = 0; index < 100000; ++index) {
603        SkDPoint origin = {ran.nextRangeF(-1000, 1000), ran.nextRangeF(-1000, 1000)};
604        SkDQuad quad1 = {{origin, {ran.nextRangeF(-1000, 1000), ran.nextRangeF(-1000, 1000)},
605            {ran.nextRangeF(-1000, 1000), ran.nextRangeF(-1000, 1000)}}};
606        if (quad1[0] == quad1[2]) {
607            continue;
608        }
609        SkDQuad quad2 = {{origin, {ran.nextRangeF(-1000, 1000), ran.nextRangeF(-1000, 1000)},
610            {ran.nextRangeF(-1000, 1000), ran.nextRangeF(-1000, 1000)}}};
611        if (quad2[0] == quad2[2]) {
612            continue;
613        }
614        SkIntersections i;
615        i.intersect(quad1, quad2);
616        REPORTER_ASSERT(reporter, i.used() >= 1);
617        if (i.used() > 1) {
618            continue;
619        }
620        TRange lowerRange, upperRange;
621        bool result = bruteMinT(reporter, quad1, quad2, &lowerRange, &upperRange);
622        REPORTER_ASSERT(reporter, result);
623        double min = SkTMin(upperRange.t1, upperRange.t2);
624        if (smaller > min) {
625            small[0] = quad1;
626            small[1] = quad2;
627            SkDEBUGCODE(smallIndex = index);
628            smaller = min;
629        }
630    }
631#ifdef SK_DEBUG
632    DumpQ(small[0], small[1], smallIndex);
633#endif
634}
635
636DEF_TEST(PathOpsAngleBruteTOne, reporter) {
637//    gPathOpsAngleIdeasVerbose = true;
638    const SkDQuad quads[] = {
639{{{-770.8492431640625, 948.2369384765625}, {-853.37066650390625, 972.0301513671875}, {-200.62042236328125, -26.7174072265625}}},
640{{{-770.8492431640625, 948.2369384765625}, {513.602783203125, 578.8681640625}, {960.641357421875, -813.69757080078125}}},
641{{{563.8267822265625, -107.4566650390625}, {-44.67724609375, -136.57452392578125}, {492.3856201171875, -268.79644775390625}}},
642{{{563.8267822265625, -107.4566650390625}, {708.049072265625, -100.77789306640625}, {-48.88226318359375, 967.9022216796875}}},
643{{{598.857421875, 846.345458984375}, {-644.095703125, -316.12921142578125}, {-97.64599609375, 20.6158447265625}}},
644{{{598.857421875, 846.345458984375}, {715.7142333984375, 955.3599853515625}, {-919.9478759765625, 691.611328125}}},
645    };
646    TRange lowerRange, upperRange;
647    bruteMinT(reporter, quads[0], quads[1], &lowerRange, &upperRange);
648    bruteMinT(reporter, quads[2], quads[3], &lowerRange, &upperRange);
649    bruteMinT(reporter, quads[4], quads[5], &lowerRange, &upperRange);
650}
651
652/*
653The sorting problem happens when the inital tangents are not a true indicator of the curve direction
654Nearly always, the initial tangents do give the right answer,
655so the trick is to figure out when the initial tangent cannot be trusted.
656If the convex hulls of both curves are in the same half plane, and not overlapping, sorting the
657hulls is enough.
658If the hulls overlap, and have the same general direction, then intersect the shorter end point ray
659with the opposing curve, and see on which side of the shorter curve the opposing intersection lies.
660Otherwise, if the control vector is extremely short, likely the point on curve must be computed
661If moving the control point slightly can change the sign of the cross product, either answer could
662be "right".
663We need to determine how short is extremely short. Move the control point a set percentage of
664the largest length to determine how stable the curve is vis-a-vis the initial tangent.
665*/
666
667static const SkDQuad extremeTests[][2] = {
668    {
669        {{{-708.0077926931004,-154.61669472244046},
670            {-707.9234268635319,-154.30459999551294},
671            {505.58447265625,-504.9130859375}}},
672        {{{-708.0077926931004,-154.61669472244046},
673            {-711.127526325141,-163.9446090624656},
674            {-32.39227294921875,-906.3277587890625}}},
675    }, {
676        {{{-708.0077926931004,-154.61669472244046},
677            {-708.2875337527566,-154.36676458635623},
678            {505.58447265625,-504.9130859375}}},
679        {{{-708.0077926931004,-154.61669472244046},
680            {-708.4111557216864,-154.5366642875255},
681            {-32.39227294921875,-906.3277587890625}}},
682    }, {
683        {{{-609.0230951752058,-267.5435593490574},
684            {-594.1120809906336,-136.08492475411555},
685            {505.58447265625,-504.9130859375}}},
686        {{{-609.0230951752058,-267.5435593490574},
687            {-693.7467719138988,-341.3259237831895},
688            {-32.39227294921875,-906.3277587890625}}}
689    }, {
690        {{{-708.0077926931004,-154.61669472244046},
691            {-707.9994640658723,-154.58588461064852},
692            {505.58447265625,-504.9130859375}}},
693        {{{-708.0077926931004,-154.61669472244046},
694            {-708.0239418990758,-154.6403553507124},
695            {-32.39227294921875,-906.3277587890625}}}
696    }, {
697        {{{-708.0077926931004,-154.61669472244046},
698            {-707.9993222215099,-154.55999389855003},
699            {68.88981098017803,296.9273945411635}}},
700        {{{-708.0077926931004,-154.61669472244046},
701            {-708.0509091919608,-154.64675214697067},
702            {-777.4871194247767,-995.1470120113145}}}
703    }, {
704        {{{-708.0077926931004,-154.61669472244046},
705            {-708.0060491116379,-154.60889321524968},
706            {229.97088707895057,-430.0569357467175}}},
707        {{{-708.0077926931004,-154.61669472244046},
708            {-708.013911296257,-154.6219143988058},
709            {138.13162892614037,-573.3689311737394}}}
710    }, {
711        {{{-543.2570545751013,-237.29243831075053},
712            {-452.4119186056987,-143.47223056267802},
713            {229.97088707895057,-430.0569357467175}}},
714        {{{-543.2570545751013,-237.29243831075053},
715            {-660.5330371214436,-362.0016148388},
716            {138.13162892614037,-573.3689311737394}}},
717    },
718};
719
720static double endCtrlRatio(const SkDQuad quad) {
721    SkDVector longEdge = quad[2] - quad[0];
722    double longLen = longEdge.length();
723    SkDVector shortEdge = quad[1] - quad[0];
724    double shortLen = shortEdge.length();
725    return longLen / shortLen;
726}
727
728static void computeMV(const SkDQuad& quad, const SkDVector& v, double m, SkDVector mV[2]) {
729        SkDPoint mPta = {quad[1].fX - m * v.fY, quad[1].fY + m * v.fX};
730        SkDPoint mPtb = {quad[1].fX + m * v.fY, quad[1].fY - m * v.fX};
731        mV[0] = mPta - quad[0];
732        mV[1] = mPtb - quad[0];
733}
734
735static double mDistance(skiatest::Reporter* reporter, bool agrees, const SkDQuad& q1,
736        const SkDQuad& q2) {
737    if (1 && agrees) {
738        return SK_ScalarMax;
739    }
740    // how close is the angle from inflecting in the opposite direction?
741    SkDVector v1 = q1[1] - q1[0];
742    SkDVector v2 = q2[1] - q2[0];
743    double dir = v1.crossCheck(v2);
744    REPORTER_ASSERT(reporter, dir != 0);
745    // solve for opposite direction displacement scale factor == m
746    // initial dir = v1.cross(v2) == v2.x * v1.y - v2.y * v1.x
747    // displacement of q1[1] : dq1 = { -m * v1.y, m * v1.x } + q1[1]
748    // straight angle when : v2.x * (dq1.y - q1[0].y) == v2.y * (dq1.x - q1[0].x)
749    //                       v2.x * (m * v1.x + v1.y) == v2.y * (-m * v1.y + v1.x)
750    // - m * (v2.x * v1.x + v2.y * v1.y) == v2.x * v1.y - v2.y * v1.x
751    // m = (v2.y * v1.x - v2.x * v1.y) / (v2.x * v1.x + v2.y * v1.y)
752    // m = v1.cross(v2) / v1.dot(v2)
753    double div = v1.dot(v2);
754    REPORTER_ASSERT(reporter, div != 0);
755    double m = dir / div;
756    SkDVector mV1[2], mV2[2];
757    computeMV(q1, v1, m, mV1);
758    computeMV(q2, v2, m, mV2);
759    double dist1 = v1.length() * m;
760    double dist2 = v2.length() * m;
761    if (gPathOpsAngleIdeasVerbose) {
762        SkDebugf("%c r1=%1.9g r2=%1.9g m=%1.9g dist1=%1.9g dist2=%1.9g "
763                " dir%c 1a=%1.9g 1b=%1.9g 2a=%1.9g 2b=%1.9g\n", agrees ? 'T' : 'F',
764                endCtrlRatio(q1), endCtrlRatio(q2), m, dist1, dist2, dir > 0 ? '+' : '-',
765                mV1[0].crossCheck(v2), mV1[1].crossCheck(v2),
766                mV2[0].crossCheck(v1), mV2[1].crossCheck(v1));
767    }
768    if (1) {
769        bool use1 = fabs(dist1) < fabs(dist2);
770        if (gPathOpsAngleIdeasVerbose) {
771            SkDebugf("%c dist=%1.9g r=%1.9g\n", agrees ? 'T' : 'F', use1 ? dist1 : dist2,
772                use1 ? distEndRatio(dist1, q1) : distEndRatio(dist2, q2));
773        }
774        return fabs(use1 ? distEndRatio(dist1, q1) : distEndRatio(dist2, q2));
775    }
776    return SK_ScalarMax;
777}
778
779static void midPointAgrees(skiatest::Reporter* reporter, const SkDQuad& q1, const SkDQuad& q2,
780        bool ccw) {
781    SkDPoint mid1 = q1.ptAtT(0.5);
782    SkDVector m1 = mid1 - q1[0];
783    SkDPoint mid2 = q2.ptAtT(0.5);
784    SkDVector m2 = mid2 - q2[0];
785    REPORTER_ASSERT(reporter, ccw ? m1.crossCheck(m2) < 0 : m1.crossCheck(m2) > 0);
786}
787
788DEF_TEST(PathOpsAngleExtreme, reporter) {
789    if (!gPathOpsAngleIdeasVerbose) {  // takes a while to run -- so exclude it by default
790        return;
791    }
792    double maxR = SK_ScalarMax;
793    for (int index = 0; index < (int) SK_ARRAY_COUNT(extremeTests); ++index) {
794        const SkDQuad& quad1 = extremeTests[index][0];
795        const SkDQuad& quad2 = extremeTests[index][1];
796        if (gPathOpsAngleIdeasVerbose) {
797            SkDebugf("%s %d\n", __FUNCTION__, index);
798        }
799        REPORTER_ASSERT(reporter, quad1[0] == quad2[0]);
800        SkIntersections i;
801        i.intersect(quad1, quad2);
802        REPORTER_ASSERT(reporter, i.used() == 1);
803        REPORTER_ASSERT(reporter, i.pt(0) == quad1[0]);
804        int overlap = quadHullsOverlap(reporter, quad1, quad2);
805        REPORTER_ASSERT(reporter, overlap >= 0);
806        SkDVector sweep[2], tweep[2];
807        setQuadHullSweep(quad1, sweep);
808        setQuadHullSweep(quad2, tweep);
809        double s0xt0 = sweep[0].crossCheck(tweep[0]);
810        REPORTER_ASSERT(reporter, s0xt0 != 0);
811        bool ccw = s0xt0 < 0;
812        bool agrees = bruteForceCheck(reporter, quad1, quad2, ccw);
813        maxR = SkTMin(maxR, mDistance(reporter, agrees, quad1, quad2));
814        if (agrees) {
815            continue;
816        }
817        midPointAgrees(reporter, quad1, quad2, !ccw);
818        SkDQuad q1 = quad1;
819        SkDQuad q2 = quad2;
820        double loFail = 1;
821        double hiPass = 2;
822        // double vectors until t passes
823        do {
824            q1[1].fX = quad1[0].fX * (1 - hiPass) + quad1[1].fX * hiPass;
825            q1[1].fY = quad1[0].fY * (1 - hiPass) + quad1[1].fY * hiPass;
826            q2[1].fX = quad2[0].fX * (1 - hiPass) + quad2[1].fX * hiPass;
827            q2[1].fY = quad2[0].fY * (1 - hiPass) + quad2[1].fY * hiPass;
828            agrees = bruteForceCheck(reporter, q1, q2, ccw);
829            maxR = SkTMin(maxR, mDistance(reporter, agrees, q1, q2));
830            if (agrees) {
831                break;
832            }
833            midPointAgrees(reporter, quad1, quad2, !ccw);
834            loFail = hiPass;
835            hiPass *= 2;
836        } while (true);
837        // binary search to find minimum pass
838        double midTest = (loFail + hiPass) / 2;
839        double step = (hiPass - loFail) / 4;
840        while (step > FLT_EPSILON) {
841            q1[1].fX = quad1[0].fX * (1 - midTest) + quad1[1].fX * midTest;
842            q1[1].fY = quad1[0].fY * (1 - midTest) + quad1[1].fY * midTest;
843            q2[1].fX = quad2[0].fX * (1 - midTest) + quad2[1].fX * midTest;
844            q2[1].fY = quad2[0].fY * (1 - midTest) + quad2[1].fY * midTest;
845            agrees = bruteForceCheck(reporter, q1, q2, ccw);
846            maxR = SkTMin(maxR, mDistance(reporter, agrees, q1, q2));
847            if (!agrees) {
848                midPointAgrees(reporter, quad1, quad2, !ccw);
849            }
850            midTest += agrees ? -step : step;
851            step /= 2;
852        }
853#ifdef SK_DEBUG
854//        DumpQ(q1, q2, 999);
855#endif
856    }
857    if (gPathOpsAngleIdeasVerbose) {
858        SkDebugf("maxR=%1.9g\n", maxR);
859    }
860}
861