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 "SkOpContour.h"
10#include "SkOpSegment.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(SkOpAngle& lh, SkOpAngle& rh) {
22        return lh.convexHullOverlaps(&rh);
23    }
24
25    static int EndsIntersect(SkOpAngle& lh, 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
409static void makeSegment(SkOpContour* contour, const SkDQuad& quad, SkPoint shortQuad[3],
410        SkChunkAlloc* allocator) {
411    shortQuad[0] = quad[0].asSkPoint();
412    shortQuad[1] = quad[1].asSkPoint();
413    shortQuad[2] = quad[2].asSkPoint();
414    contour->addQuad(shortQuad, allocator);
415}
416
417static void testQuadAngles(skiatest::Reporter* reporter, const SkDQuad& quad1, const SkDQuad& quad2,
418        int testNo, SkChunkAlloc* allocator) {
419    SkPoint shortQuads[2][3];
420
421    SkOpContourHead contour;
422    SkOpGlobalState state(nullptr, &contour  SkDEBUGPARAMS(nullptr));
423    contour.init(&state, false, false);
424    makeSegment(&contour, quad1, shortQuads[0], allocator);
425    makeSegment(&contour, quad1, shortQuads[1], allocator);
426    SkOpSegment* seg1 = contour.first();
427    seg1->debugAddAngle(0, 1, allocator);
428    SkOpSegment* seg2 = seg1->next();
429    seg2->debugAddAngle(0, 1, allocator);
430    int realOverlap = PathOpsAngleTester::ConvexHullOverlaps(*seg1->debugLastAngle(),
431            *seg2->debugLastAngle());
432    const SkDPoint& origin = quad1[0];
433    REPORTER_ASSERT(reporter, origin == quad2[0]);
434    double a1s = atan2(origin.fY - quad1[1].fY, quad1[1].fX - origin.fX);
435    double a1e = atan2(origin.fY - quad1[2].fY, quad1[2].fX - origin.fX);
436    double a2s = atan2(origin.fY - quad2[1].fY, quad2[1].fX - origin.fX);
437    double a2e = atan2(origin.fY - quad2[2].fY, quad2[2].fX - origin.fX);
438    bool oldSchoolOverlap = radianBetween(a1s, a2s, a1e)
439        || radianBetween(a1s, a2e, a1e) || radianBetween(a2s, a1s, a2e)
440        || radianBetween(a2s, a1e, a2e);
441    int overlap = quadHullsOverlap(reporter, quad1, quad2);
442    bool realMatchesOverlap = realOverlap == overlap || SK_ScalarPI - fabs(a2s - a1s) < 0.002;
443    if (realOverlap != overlap) {
444        SkDebugf("\nSK_ScalarPI - fabs(a2s - a1s) = %1.9g\n", SK_ScalarPI - fabs(a2s - a1s));
445    }
446    if (!realMatchesOverlap) {
447        DumpQ(quad1, quad2, testNo);
448    }
449    REPORTER_ASSERT(reporter, realMatchesOverlap);
450    if (oldSchoolOverlap != (overlap < 0)) {
451        overlap = quadHullsOverlap(reporter, quad1, quad2);  // set a breakpoint and debug if assert fires
452        REPORTER_ASSERT(reporter, oldSchoolOverlap == (overlap < 0));
453    }
454    SkDVector v1s = quad1[1] - quad1[0];
455    SkDVector v1e = quad1[2] - quad1[0];
456    SkDVector v2s = quad2[1] - quad2[0];
457    SkDVector v2e = quad2[2] - quad2[0];
458    double vDir[2] = { v1s.cross(v1e), v2s.cross(v2e) };
459    bool ray1In2 = v1s.cross(v2s) * vDir[1] <= 0 && v1s.cross(v2e) * vDir[1] >= 0;
460    bool ray2In1 = v2s.cross(v1s) * vDir[0] <= 0 && v2s.cross(v1e) * vDir[0] >= 0;
461    if (overlap >= 0) {
462        // verify that hulls really don't overlap
463        REPORTER_ASSERT(reporter, !ray1In2);
464        REPORTER_ASSERT(reporter, !ray2In1);
465        bool ctrl1In2 = v1e.cross(v2s) * vDir[1] <= 0 && v1e.cross(v2e) * vDir[1] >= 0;
466        REPORTER_ASSERT(reporter, !ctrl1In2);
467        bool ctrl2In1 = v2e.cross(v1s) * vDir[0] <= 0 && v2e.cross(v1e) * vDir[0] >= 0;
468        REPORTER_ASSERT(reporter, !ctrl2In1);
469        // check answer against reference
470        bruteForce(reporter, quad1, quad2, overlap > 0);
471    }
472    // continue end point rays and see if they intersect the opposite curve
473    SkDLine rays[] = {{{origin, quad2[2]}}, {{origin, quad1[2]}}};
474    const SkDQuad* quads[] = {&quad1, &quad2};
475    SkDVector midSpokes[2];
476    SkIntersections intersect[2];
477    double minX, minY, maxX, maxY;
478    minX = minY = SK_ScalarInfinity;
479    maxX = maxY = -SK_ScalarInfinity;
480    double maxWidth = 0;
481    bool useIntersect = false;
482    double smallestTs[] = {1, 1};
483    for (unsigned index = 0; index < SK_ARRAY_COUNT(quads); ++index) {
484        const SkDQuad& q = *quads[index];
485        midSpokes[index] = q.ptAtT(0.5) - origin;
486        minX = SkTMin(SkTMin(SkTMin(minX, origin.fX), q[1].fX), q[2].fX);
487        minY = SkTMin(SkTMin(SkTMin(minY, origin.fY), q[1].fY), q[2].fY);
488        maxX = SkTMax(SkTMax(SkTMax(maxX, origin.fX), q[1].fX), q[2].fX);
489        maxY = SkTMax(SkTMax(SkTMax(maxY, origin.fY), q[1].fY), q[2].fY);
490        maxWidth = SkTMax(maxWidth, SkTMax(maxX - minX, maxY - minY));
491        intersect[index].intersectRay(q, rays[index]);
492        const SkIntersections& i = intersect[index];
493        REPORTER_ASSERT(reporter, i.used() >= 1);
494        bool foundZero = false;
495        double smallT = 1;
496        for (int idx2 = 0; idx2 < i.used(); ++idx2) {
497            double t = i[0][idx2];
498            if (t == 0) {
499                foundZero = true;
500                continue;
501            }
502            if (smallT > t) {
503                smallT = t;
504            }
505        }
506        REPORTER_ASSERT(reporter, foundZero == true);
507        if (smallT == 1) {
508            continue;
509        }
510        SkDVector ray = q.ptAtT(smallT) - origin;
511        SkDVector end = rays[index][1] - origin;
512        if (ray.fX * end.fX < 0 || ray.fY * end.fY < 0) {
513            continue;
514        }
515        double rayDist = ray.length();
516        double endDist = end.length();
517        double delta = fabs(rayDist - endDist) / maxWidth;
518        if (delta > 1e-4) {
519            useIntersect ^= true;
520        }
521        smallestTs[index] = smallT;
522    }
523    bool firstInside;
524    if (useIntersect) {
525        int sIndex = (int) (smallestTs[1] < 1);
526        REPORTER_ASSERT(reporter, smallestTs[sIndex ^ 1] == 1);
527        double t = smallestTs[sIndex];
528        const SkDQuad& q = *quads[sIndex];
529        SkDVector ray = q.ptAtT(t) - origin;
530        SkDVector end = rays[sIndex][1] - origin;
531        double rayDist = ray.length();
532        double endDist = end.length();
533        SkDVector mid = q.ptAtT(t / 2) - origin;
534        double midXray = mid.crossCheck(ray);
535        if (gPathOpsAngleIdeasVerbose) {
536            SkDebugf("rayDist>endDist:%d sIndex==0:%d vDir[sIndex]<0:%d midXray<0:%d\n",
537                    rayDist > endDist, sIndex == 0, vDir[sIndex] < 0, midXray < 0);
538        }
539        SkASSERT(SkScalarSignAsInt(SkDoubleToScalar(midXray))
540            == SkScalarSignAsInt(SkDoubleToScalar(vDir[sIndex])));
541        firstInside = (rayDist > endDist) ^ (sIndex == 0) ^ (vDir[sIndex] < 0);
542    } else if (overlap >= 0) {
543        return;  // answer has already been determined
544    } else {
545        firstInside = checkParallel(reporter, quad1, quad2);
546    }
547    if (overlap < 0) {
548        SkDEBUGCODE(int realEnds =)
549                PathOpsAngleTester::EndsIntersect(*seg1->debugLastAngle(),
550                *seg2->debugLastAngle());
551        SkASSERT(realEnds == (firstInside ? 1 : 0));
552    }
553    bruteForce(reporter, quad1, quad2, firstInside);
554}
555
556DEF_TEST(PathOpsAngleOverlapHullsOne, reporter) {
557    SkChunkAlloc allocator(4096);
558//    gPathOpsAngleIdeasVerbose = true;
559    const SkDQuad quads[] = {
560{{{939.4808349609375, 914.355224609375}, {-357.7921142578125, 590.842529296875}, {736.8936767578125, -350.717529296875}}},
561{{{939.4808349609375, 914.355224609375}, {-182.85418701171875, 634.4552001953125}, {-509.62615966796875, 576.1182861328125}}}
562    };
563    for (int index = 0; index < (int) SK_ARRAY_COUNT(quads); index += 2) {
564        testQuadAngles(reporter, quads[index], quads[index + 1], 0, &allocator);
565    }
566}
567
568DEF_TEST(PathOpsAngleOverlapHulls, reporter) {
569    SkChunkAlloc allocator(4096);
570    if (!gPathOpsAngleIdeasVerbose) {  // takes a while to run -- so exclude it by default
571        return;
572    }
573    SkRandom ran;
574    for (int index = 0; index < 100000; ++index) {
575        if (index % 1000 == 999) SkDebugf(".");
576        SkDPoint origin = {ran.nextRangeF(-1000, 1000), ran.nextRangeF(-1000, 1000)};
577        SkDQuad quad1 = {{origin, {ran.nextRangeF(-1000, 1000), ran.nextRangeF(-1000, 1000)},
578            {ran.nextRangeF(-1000, 1000), ran.nextRangeF(-1000, 1000)}}};
579        if (quad1[0] == quad1[2]) {
580            continue;
581        }
582        SkDQuad quad2 = {{origin, {ran.nextRangeF(-1000, 1000), ran.nextRangeF(-1000, 1000)},
583            {ran.nextRangeF(-1000, 1000), ran.nextRangeF(-1000, 1000)}}};
584        if (quad2[0] == quad2[2]) {
585            continue;
586        }
587        SkIntersections i;
588        i.intersect(quad1, quad2);
589        REPORTER_ASSERT(reporter, i.used() >= 1);
590        if (i.used() > 1) {
591            continue;
592        }
593        testQuadAngles(reporter, quad1, quad2, index, &allocator);
594    }
595}
596
597DEF_TEST(PathOpsAngleBruteT, reporter) {
598    if (!gPathOpsAngleIdeasVerbose) {  // takes a while to run -- so exclude it by default
599        return;
600    }
601    SkRandom ran;
602    double smaller = SK_Scalar1;
603    SkDQuad small[2];
604    SkDEBUGCODE(int smallIndex);
605    for (int index = 0; index < 100000; ++index) {
606        SkDPoint origin = {ran.nextRangeF(-1000, 1000), ran.nextRangeF(-1000, 1000)};
607        SkDQuad quad1 = {{origin, {ran.nextRangeF(-1000, 1000), ran.nextRangeF(-1000, 1000)},
608            {ran.nextRangeF(-1000, 1000), ran.nextRangeF(-1000, 1000)}}};
609        if (quad1[0] == quad1[2]) {
610            continue;
611        }
612        SkDQuad quad2 = {{origin, {ran.nextRangeF(-1000, 1000), ran.nextRangeF(-1000, 1000)},
613            {ran.nextRangeF(-1000, 1000), ran.nextRangeF(-1000, 1000)}}};
614        if (quad2[0] == quad2[2]) {
615            continue;
616        }
617        SkIntersections i;
618        i.intersect(quad1, quad2);
619        REPORTER_ASSERT(reporter, i.used() >= 1);
620        if (i.used() > 1) {
621            continue;
622        }
623        TRange lowerRange, upperRange;
624        bool result = bruteMinT(reporter, quad1, quad2, &lowerRange, &upperRange);
625        REPORTER_ASSERT(reporter, result);
626        double min = SkTMin(upperRange.t1, upperRange.t2);
627        if (smaller > min) {
628            small[0] = quad1;
629            small[1] = quad2;
630            SkDEBUGCODE(smallIndex = index);
631            smaller = min;
632        }
633    }
634#ifdef SK_DEBUG
635    DumpQ(small[0], small[1], smallIndex);
636#endif
637}
638
639DEF_TEST(PathOpsAngleBruteTOne, reporter) {
640//    gPathOpsAngleIdeasVerbose = true;
641    const SkDQuad quads[] = {
642{{{-770.8492431640625, 948.2369384765625}, {-853.37066650390625, 972.0301513671875}, {-200.62042236328125, -26.7174072265625}}},
643{{{-770.8492431640625, 948.2369384765625}, {513.602783203125, 578.8681640625}, {960.641357421875, -813.69757080078125}}},
644{{{563.8267822265625, -107.4566650390625}, {-44.67724609375, -136.57452392578125}, {492.3856201171875, -268.79644775390625}}},
645{{{563.8267822265625, -107.4566650390625}, {708.049072265625, -100.77789306640625}, {-48.88226318359375, 967.9022216796875}}},
646{{{598.857421875, 846.345458984375}, {-644.095703125, -316.12921142578125}, {-97.64599609375, 20.6158447265625}}},
647{{{598.857421875, 846.345458984375}, {715.7142333984375, 955.3599853515625}, {-919.9478759765625, 691.611328125}}},
648    };
649    TRange lowerRange, upperRange;
650    bruteMinT(reporter, quads[0], quads[1], &lowerRange, &upperRange);
651    bruteMinT(reporter, quads[2], quads[3], &lowerRange, &upperRange);
652    bruteMinT(reporter, quads[4], quads[5], &lowerRange, &upperRange);
653}
654
655/*
656The sorting problem happens when the inital tangents are not a true indicator of the curve direction
657Nearly always, the initial tangents do give the right answer,
658so the trick is to figure out when the initial tangent cannot be trusted.
659If the convex hulls of both curves are in the same half plane, and not overlapping, sorting the
660hulls is enough.
661If the hulls overlap, and have the same general direction, then intersect the shorter end point ray
662with the opposing curve, and see on which side of the shorter curve the opposing intersection lies.
663Otherwise, if the control vector is extremely short, likely the point on curve must be computed
664If moving the control point slightly can change the sign of the cross product, either answer could
665be "right".
666We need to determine how short is extremely short. Move the control point a set percentage of
667the largest length to determine how stable the curve is vis-a-vis the initial tangent.
668*/
669
670static const SkDQuad extremeTests[][2] = {
671    {
672        {{{-708.0077926931004,-154.61669472244046},
673            {-707.9234268635319,-154.30459999551294},
674            {505.58447265625,-504.9130859375}}},
675        {{{-708.0077926931004,-154.61669472244046},
676            {-711.127526325141,-163.9446090624656},
677            {-32.39227294921875,-906.3277587890625}}},
678    }, {
679        {{{-708.0077926931004,-154.61669472244046},
680            {-708.2875337527566,-154.36676458635623},
681            {505.58447265625,-504.9130859375}}},
682        {{{-708.0077926931004,-154.61669472244046},
683            {-708.4111557216864,-154.5366642875255},
684            {-32.39227294921875,-906.3277587890625}}},
685    }, {
686        {{{-609.0230951752058,-267.5435593490574},
687            {-594.1120809906336,-136.08492475411555},
688            {505.58447265625,-504.9130859375}}},
689        {{{-609.0230951752058,-267.5435593490574},
690            {-693.7467719138988,-341.3259237831895},
691            {-32.39227294921875,-906.3277587890625}}}
692    }, {
693        {{{-708.0077926931004,-154.61669472244046},
694            {-707.9994640658723,-154.58588461064852},
695            {505.58447265625,-504.9130859375}}},
696        {{{-708.0077926931004,-154.61669472244046},
697            {-708.0239418990758,-154.6403553507124},
698            {-32.39227294921875,-906.3277587890625}}}
699    }, {
700        {{{-708.0077926931004,-154.61669472244046},
701            {-707.9993222215099,-154.55999389855003},
702            {68.88981098017803,296.9273945411635}}},
703        {{{-708.0077926931004,-154.61669472244046},
704            {-708.0509091919608,-154.64675214697067},
705            {-777.4871194247767,-995.1470120113145}}}
706    }, {
707        {{{-708.0077926931004,-154.61669472244046},
708            {-708.0060491116379,-154.60889321524968},
709            {229.97088707895057,-430.0569357467175}}},
710        {{{-708.0077926931004,-154.61669472244046},
711            {-708.013911296257,-154.6219143988058},
712            {138.13162892614037,-573.3689311737394}}}
713    }, {
714        {{{-543.2570545751013,-237.29243831075053},
715            {-452.4119186056987,-143.47223056267802},
716            {229.97088707895057,-430.0569357467175}}},
717        {{{-543.2570545751013,-237.29243831075053},
718            {-660.5330371214436,-362.0016148388},
719            {138.13162892614037,-573.3689311737394}}},
720    },
721};
722
723static double endCtrlRatio(const SkDQuad quad) {
724    SkDVector longEdge = quad[2] - quad[0];
725    double longLen = longEdge.length();
726    SkDVector shortEdge = quad[1] - quad[0];
727    double shortLen = shortEdge.length();
728    return longLen / shortLen;
729}
730
731static void computeMV(const SkDQuad& quad, const SkDVector& v, double m, SkDVector mV[2]) {
732        SkDPoint mPta = {quad[1].fX - m * v.fY, quad[1].fY + m * v.fX};
733        SkDPoint mPtb = {quad[1].fX + m * v.fY, quad[1].fY - m * v.fX};
734        mV[0] = mPta - quad[0];
735        mV[1] = mPtb - quad[0];
736}
737
738static double mDistance(skiatest::Reporter* reporter, bool agrees, const SkDQuad& q1,
739        const SkDQuad& q2) {
740    if (1 && agrees) {
741        return SK_ScalarMax;
742    }
743    // how close is the angle from inflecting in the opposite direction?
744    SkDVector v1 = q1[1] - q1[0];
745    SkDVector v2 = q2[1] - q2[0];
746    double dir = v1.crossCheck(v2);
747    REPORTER_ASSERT(reporter, dir != 0);
748    // solve for opposite direction displacement scale factor == m
749    // initial dir = v1.cross(v2) == v2.x * v1.y - v2.y * v1.x
750    // displacement of q1[1] : dq1 = { -m * v1.y, m * v1.x } + q1[1]
751    // straight angle when : v2.x * (dq1.y - q1[0].y) == v2.y * (dq1.x - q1[0].x)
752    //                       v2.x * (m * v1.x + v1.y) == v2.y * (-m * v1.y + v1.x)
753    // - m * (v2.x * v1.x + v2.y * v1.y) == v2.x * v1.y - v2.y * v1.x
754    // m = (v2.y * v1.x - v2.x * v1.y) / (v2.x * v1.x + v2.y * v1.y)
755    // m = v1.cross(v2) / v1.dot(v2)
756    double div = v1.dot(v2);
757    REPORTER_ASSERT(reporter, div != 0);
758    double m = dir / div;
759    SkDVector mV1[2], mV2[2];
760    computeMV(q1, v1, m, mV1);
761    computeMV(q2, v2, m, mV2);
762    double dist1 = v1.length() * m;
763    double dist2 = v2.length() * m;
764    if (gPathOpsAngleIdeasVerbose) {
765        SkDebugf("%c r1=%1.9g r2=%1.9g m=%1.9g dist1=%1.9g dist2=%1.9g "
766                " dir%c 1a=%1.9g 1b=%1.9g 2a=%1.9g 2b=%1.9g\n", agrees ? 'T' : 'F',
767                endCtrlRatio(q1), endCtrlRatio(q2), m, dist1, dist2, dir > 0 ? '+' : '-',
768                mV1[0].crossCheck(v2), mV1[1].crossCheck(v2),
769                mV2[0].crossCheck(v1), mV2[1].crossCheck(v1));
770    }
771    if (1) {
772        bool use1 = fabs(dist1) < fabs(dist2);
773        if (gPathOpsAngleIdeasVerbose) {
774            SkDebugf("%c dist=%1.9g r=%1.9g\n", agrees ? 'T' : 'F', use1 ? dist1 : dist2,
775                use1 ? distEndRatio(dist1, q1) : distEndRatio(dist2, q2));
776        }
777        return fabs(use1 ? distEndRatio(dist1, q1) : distEndRatio(dist2, q2));
778    }
779    return SK_ScalarMax;
780}
781
782static void midPointAgrees(skiatest::Reporter* reporter, const SkDQuad& q1, const SkDQuad& q2,
783        bool ccw) {
784    SkDPoint mid1 = q1.ptAtT(0.5);
785    SkDVector m1 = mid1 - q1[0];
786    SkDPoint mid2 = q2.ptAtT(0.5);
787    SkDVector m2 = mid2 - q2[0];
788    REPORTER_ASSERT(reporter, ccw ? m1.crossCheck(m2) < 0 : m1.crossCheck(m2) > 0);
789}
790
791DEF_TEST(PathOpsAngleExtreme, reporter) {
792    if (!gPathOpsAngleIdeasVerbose) {  // takes a while to run -- so exclude it by default
793        return;
794    }
795    double maxR = SK_ScalarMax;
796    for (int index = 0; index < (int) SK_ARRAY_COUNT(extremeTests); ++index) {
797        const SkDQuad& quad1 = extremeTests[index][0];
798        const SkDQuad& quad2 = extremeTests[index][1];
799        if (gPathOpsAngleIdeasVerbose) {
800            SkDebugf("%s %d\n", __FUNCTION__, index);
801        }
802        REPORTER_ASSERT(reporter, quad1[0] == quad2[0]);
803        SkIntersections i;
804        i.intersect(quad1, quad2);
805        REPORTER_ASSERT(reporter, i.used() == 1);
806        REPORTER_ASSERT(reporter, i.pt(0) == quad1[0]);
807        int overlap = quadHullsOverlap(reporter, quad1, quad2);
808        REPORTER_ASSERT(reporter, overlap >= 0);
809        SkDVector sweep[2], tweep[2];
810        setQuadHullSweep(quad1, sweep);
811        setQuadHullSweep(quad2, tweep);
812        double s0xt0 = sweep[0].crossCheck(tweep[0]);
813        REPORTER_ASSERT(reporter, s0xt0 != 0);
814        bool ccw = s0xt0 < 0;
815        bool agrees = bruteForceCheck(reporter, quad1, quad2, ccw);
816        maxR = SkTMin(maxR, mDistance(reporter, agrees, quad1, quad2));
817        if (agrees) {
818            continue;
819        }
820        midPointAgrees(reporter, quad1, quad2, !ccw);
821        SkDQuad q1 = quad1;
822        SkDQuad q2 = quad2;
823        double loFail = 1;
824        double hiPass = 2;
825        // double vectors until t passes
826        do {
827            q1[1].fX = quad1[0].fX * (1 - hiPass) + quad1[1].fX * hiPass;
828            q1[1].fY = quad1[0].fY * (1 - hiPass) + quad1[1].fY * hiPass;
829            q2[1].fX = quad2[0].fX * (1 - hiPass) + quad2[1].fX * hiPass;
830            q2[1].fY = quad2[0].fY * (1 - hiPass) + quad2[1].fY * hiPass;
831            agrees = bruteForceCheck(reporter, q1, q2, ccw);
832            maxR = SkTMin(maxR, mDistance(reporter, agrees, q1, q2));
833            if (agrees) {
834                break;
835            }
836            midPointAgrees(reporter, quad1, quad2, !ccw);
837            loFail = hiPass;
838            hiPass *= 2;
839        } while (true);
840        // binary search to find minimum pass
841        double midTest = (loFail + hiPass) / 2;
842        double step = (hiPass - loFail) / 4;
843        while (step > FLT_EPSILON) {
844            q1[1].fX = quad1[0].fX * (1 - midTest) + quad1[1].fX * midTest;
845            q1[1].fY = quad1[0].fY * (1 - midTest) + quad1[1].fY * midTest;
846            q2[1].fX = quad2[0].fX * (1 - midTest) + quad2[1].fX * midTest;
847            q2[1].fY = quad2[0].fY * (1 - midTest) + quad2[1].fY * midTest;
848            agrees = bruteForceCheck(reporter, q1, q2, ccw);
849            maxR = SkTMin(maxR, mDistance(reporter, agrees, q1, q2));
850            if (!agrees) {
851                midPointAgrees(reporter, quad1, quad2, !ccw);
852            }
853            midTest += agrees ? -step : step;
854            step /= 2;
855        }
856#ifdef SK_DEBUG
857//        DumpQ(q1, q2, 999);
858#endif
859    }
860    if (gPathOpsAngleIdeasVerbose) {
861        SkDebugf("maxR=%1.9g\n", maxR);
862    }
863}
864