Home | History | Annotate | Download | only in pathops
      1 // Another approach is to start with the implicit form of one curve and solve
      2 // (seek implicit coefficients in QuadraticParameter.cpp
      3 // by substituting in the parametric form of the other.
      4 // The downside of this approach is that early rejects are difficult to come by.
      5 // http://planetmath.org/encyclopedia/GaloisTheoreticDerivationOfTheQuarticFormula.html#step
      6 
      7 #include "SkDQuadImplicit.h"
      8 #include "SkIntersections.h"
      9 #include "SkPathOpsLine.h"
     10 #include "SkQuarticRoot.h"
     11 #include "SkTArray.h"
     12 #include "SkTSort.h"
     13 
     14 /* given the implicit form 0 = Ax^2 + Bxy + Cy^2 + Dx + Ey + F
     15  * and given x = at^2 + bt + c  (the parameterized form)
     16  *           y = dt^2 + et + f
     17  * then
     18  * 0 = A(at^2+bt+c)(at^2+bt+c)+B(at^2+bt+c)(dt^2+et+f)+C(dt^2+et+f)(dt^2+et+f)+D(at^2+bt+c)+E(dt^2+et+f)+F
     19  */
     20 
     21 static int findRoots(const SkDQuadImplicit& i, const SkDQuad& quad, double roots[4],
     22         bool oneHint, bool flip, int firstCubicRoot) {
     23     SkDQuad flipped;
     24     const SkDQuad& q = flip ? (flipped = quad.flip()) : quad;
     25     double a, b, c;
     26     SkDQuad::SetABC(&q[0].fX, &a, &b, &c);
     27     double d, e, f;
     28     SkDQuad::SetABC(&q[0].fY, &d, &e, &f);
     29     const double t4 =     i.x2() *  a * a
     30                     +     i.xy() *  a * d
     31                     +     i.y2() *  d * d;
     32     const double t3 = 2 * i.x2() *  a * b
     33                     +     i.xy() * (a * e +     b * d)
     34                     + 2 * i.y2() *  d * e;
     35     const double t2 =     i.x2() * (b * b + 2 * a * c)
     36                     +     i.xy() * (c * d +     b * e + a * f)
     37                     +     i.y2() * (e * e + 2 * d * f)
     38                     +     i.x()  *  a
     39                     +     i.y()  *  d;
     40     const double t1 = 2 * i.x2() *  b * c
     41                     +     i.xy() * (c * e + b * f)
     42                     + 2 * i.y2() *  e * f
     43                     +     i.x()  *  b
     44                     +     i.y()  *  e;
     45     const double t0 =     i.x2() *  c * c
     46                     +     i.xy() *  c * f
     47                     +     i.y2() *  f * f
     48                     +     i.x()  *  c
     49                     +     i.y()  *  f
     50                     +     i.c();
     51     int rootCount = SkReducedQuarticRoots(t4, t3, t2, t1, t0, oneHint, roots);
     52     if (rootCount < 0) {
     53         rootCount = SkQuarticRootsReal(firstCubicRoot, t4, t3, t2, t1, t0, roots);
     54     }
     55     if (flip) {
     56         for (int index = 0; index < rootCount; ++index) {
     57             roots[index] = 1 - roots[index];
     58         }
     59     }
     60     return rootCount;
     61 }
     62 
     63 static int addValidRoots(const double roots[4], const int count, double valid[4]) {
     64     int result = 0;
     65     int index;
     66     for (index = 0; index < count; ++index) {
     67         if (!approximately_zero_or_more(roots[index]) || !approximately_one_or_less(roots[index])) {
     68             continue;
     69         }
     70         double t = 1 - roots[index];
     71         if (approximately_less_than_zero(t)) {
     72             t = 0;
     73         } else if (approximately_greater_than_one(t)) {
     74             t = 1;
     75         }
     76         valid[result++] = t;
     77     }
     78     return result;
     79 }
     80 
     81 static bool only_end_pts_in_common(const SkDQuad& q1, const SkDQuad& q2) {
     82 // the idea here is to see at minimum do a quick reject by rotating all points
     83 // to either side of the line formed by connecting the endpoints
     84 // if the opposite curves points are on the line or on the other side, the
     85 // curves at most intersect at the endpoints
     86     for (int oddMan = 0; oddMan < 3; ++oddMan) {
     87         const SkDPoint* endPt[2];
     88         for (int opp = 1; opp < 3; ++opp) {
     89             int end = oddMan ^ opp;  // choose a value not equal to oddMan
     90             if (3 == end) {  // and correct so that largest value is 1 or 2
     91                 end = opp;
     92             }
     93             endPt[opp - 1] = &q1[end];
     94         }
     95         double origX = endPt[0]->fX;
     96         double origY = endPt[0]->fY;
     97         double adj = endPt[1]->fX - origX;
     98         double opp = endPt[1]->fY - origY;
     99         double sign = (q1[oddMan].fY - origY) * adj - (q1[oddMan].fX - origX) * opp;
    100         if (approximately_zero(sign)) {
    101             goto tryNextHalfPlane;
    102         }
    103         for (int n = 0; n < 3; ++n) {
    104             double test = (q2[n].fY - origY) * adj - (q2[n].fX - origX) * opp;
    105             if (test * sign > 0 && !precisely_zero(test)) {
    106                 goto tryNextHalfPlane;
    107             }
    108         }
    109         return true;
    110 tryNextHalfPlane:
    111         ;
    112     }
    113     return false;
    114 }
    115 
    116 // returns false if there's more than one intercept or the intercept doesn't match the point
    117 // returns true if the intercept was successfully added or if the
    118 // original quads need to be subdivided
    119 static bool add_intercept(const SkDQuad& q1, const SkDQuad& q2, double tMin, double tMax,
    120                           SkIntersections* i, bool* subDivide) {
    121     double tMid = (tMin + tMax) / 2;
    122     SkDPoint mid = q2.ptAtT(tMid);
    123     SkDLine line;
    124     line[0] = line[1] = mid;
    125     SkDVector dxdy = q2.dxdyAtT(tMid);
    126     line[0] -= dxdy;
    127     line[1] += dxdy;
    128     SkIntersections rootTs;
    129     rootTs.allowNear(false);
    130     int roots = rootTs.intersect(q1, line);
    131     if (roots == 0) {
    132         if (subDivide) {
    133             *subDivide = true;
    134         }
    135         return true;
    136     }
    137     if (roots == 2) {
    138         return false;
    139     }
    140     SkDPoint pt2 = q1.ptAtT(rootTs[0][0]);
    141     if (!pt2.approximatelyEqual(mid)) {
    142         return false;
    143     }
    144     i->insertSwap(rootTs[0][0], tMid, pt2);
    145     return true;
    146 }
    147 
    148 static bool is_linear_inner(const SkDQuad& q1, double t1s, double t1e, const SkDQuad& q2,
    149                             double t2s, double t2e, SkIntersections* i, bool* subDivide) {
    150     SkDQuad hull = q1.subDivide(t1s, t1e);
    151     SkDLine line = {{hull[2], hull[0]}};
    152     const SkDLine* testLines[] = { &line, (const SkDLine*) &hull[0], (const SkDLine*) &hull[1] };
    153     const size_t kTestCount = SK_ARRAY_COUNT(testLines);
    154     SkSTArray<kTestCount * 2, double, true> tsFound;
    155     for (size_t index = 0; index < kTestCount; ++index) {
    156         SkIntersections rootTs;
    157         rootTs.allowNear(false);
    158         int roots = rootTs.intersect(q2, *testLines[index]);
    159         for (int idx2 = 0; idx2 < roots; ++idx2) {
    160             double t = rootTs[0][idx2];
    161 #if 0 // def SK_DEBUG   // FIXME : accurate for error = 16, error of 17.5 seen
    162 // {{{136.08723965397621, 1648.2814535211637}, {593.49031197259478, 1190.8784277439891}, {593.49031197259478, 544.0128173828125}}}
    163 // {{{-968.181396484375, 544.0128173828125}, {592.2825927734375, 870.552490234375}, {593.435302734375, 557.8828125}}}
    164 
    165             SkDPoint qPt = q2.ptAtT(t);
    166             SkDPoint lPt = testLines[index]->ptAtT(rootTs[1][idx2]);
    167             SkASSERT(qPt.approximatelyDEqual(lPt));
    168 #endif
    169             if (approximately_negative(t - t2s) || approximately_positive(t - t2e)) {
    170                 continue;
    171             }
    172             tsFound.push_back(rootTs[0][idx2]);
    173         }
    174     }
    175     int tCount = tsFound.count();
    176     if (tCount <= 0) {
    177         return true;
    178     }
    179     double tMin, tMax;
    180     if (tCount == 1) {
    181         tMin = tMax = tsFound[0];
    182     } else {
    183         SkASSERT(tCount > 1);
    184         SkTQSort<double>(tsFound.begin(), tsFound.end() - 1);
    185         tMin = tsFound[0];
    186         tMax = tsFound[tsFound.count() - 1];
    187     }
    188     SkDPoint end = q2.ptAtT(t2s);
    189     bool startInTriangle = hull.pointInHull(end);
    190     if (startInTriangle) {
    191         tMin = t2s;
    192     }
    193     end = q2.ptAtT(t2e);
    194     bool endInTriangle = hull.pointInHull(end);
    195     if (endInTriangle) {
    196         tMax = t2e;
    197     }
    198     int split = 0;
    199     SkDVector dxy1, dxy2;
    200     if (tMin != tMax || tCount > 2) {
    201         dxy2 = q2.dxdyAtT(tMin);
    202         for (int index = 1; index < tCount; ++index) {
    203             dxy1 = dxy2;
    204             dxy2 = q2.dxdyAtT(tsFound[index]);
    205             double dot = dxy1.dot(dxy2);
    206             if (dot < 0) {
    207                 split = index - 1;
    208                 break;
    209             }
    210         }
    211     }
    212     if (split == 0) {  // there's one point
    213         if (add_intercept(q1, q2, tMin, tMax, i, subDivide)) {
    214             return true;
    215         }
    216         i->swap();
    217         return is_linear_inner(q2, tMin, tMax, q1, t1s, t1e, i, subDivide);
    218     }
    219     // At this point, we have two ranges of t values -- treat each separately at the split
    220     bool result;
    221     if (add_intercept(q1, q2, tMin, tsFound[split - 1], i, subDivide)) {
    222         result = true;
    223     } else {
    224         i->swap();
    225         result = is_linear_inner(q2, tMin, tsFound[split - 1], q1, t1s, t1e, i, subDivide);
    226     }
    227     if (add_intercept(q1, q2, tsFound[split], tMax, i, subDivide)) {
    228         result = true;
    229     } else {
    230         i->swap();
    231         result |= is_linear_inner(q2, tsFound[split], tMax, q1, t1s, t1e, i, subDivide);
    232     }
    233     return result;
    234 }
    235 
    236 static double flat_measure(const SkDQuad& q) {
    237     SkDVector mid = q[1] - q[0];
    238     SkDVector dxy = q[2] - q[0];
    239     double length = dxy.length();  // OPTIMIZE: get rid of sqrt
    240     return fabs(mid.cross(dxy) / length);
    241 }
    242 
    243 // FIXME ? should this measure both and then use the quad that is the flattest as the line?
    244 static bool is_linear(const SkDQuad& q1, const SkDQuad& q2, SkIntersections* i) {
    245     double measure = flat_measure(q1);
    246     // OPTIMIZE: (get rid of sqrt) use approximately_zero
    247     if (!approximately_zero_sqrt(measure)) {
    248         return false;
    249     }
    250     return is_linear_inner(q1, 0, 1, q2, 0, 1, i, NULL);
    251 }
    252 
    253 // FIXME: if flat measure is sufficiently large, then probably the quartic solution failed
    254 // avoid imprecision incurred with chopAt
    255 static void relaxed_is_linear(const SkDQuad* q1, double s1, double e1, const SkDQuad* q2,
    256         double s2, double e2, SkIntersections* i) {
    257     double m1 = flat_measure(*q1);
    258     double m2 = flat_measure(*q2);
    259     i->reset();
    260     const SkDQuad* rounder, *flatter;
    261     double sf, midf, ef, sr, er;
    262     if (m2 < m1) {
    263         rounder = q1;
    264         sr = s1;
    265         er = e1;
    266         flatter = q2;
    267         sf = s2;
    268         midf = (s2 + e2) / 2;
    269         ef = e2;
    270     } else {
    271         rounder = q2;
    272         sr = s2;
    273         er = e2;
    274         flatter = q1;
    275         sf = s1;
    276         midf = (s1 + e1) / 2;
    277         ef = e1;
    278     }
    279     bool subDivide = false;
    280     is_linear_inner(*flatter, sf, ef, *rounder, sr, er, i, &subDivide);
    281     if (subDivide) {
    282         relaxed_is_linear(flatter, sf, midf, rounder, sr, er, i);
    283         relaxed_is_linear(flatter, midf, ef, rounder, sr, er, i);
    284     }
    285     if (m2 < m1) {
    286         i->swapPts();
    287     }
    288 }
    289 
    290 // each time through the loop, this computes values it had from the last loop
    291 // if i == j == 1, the center values are still good
    292 // otherwise, for i != 1 or j != 1, four of the values are still good
    293 // and if i == 1 ^ j == 1, an additional value is good
    294 static bool binary_search(const SkDQuad& quad1, const SkDQuad& quad2, double* t1Seed,
    295                           double* t2Seed, SkDPoint* pt) {
    296     double tStep = ROUGH_EPSILON;
    297     SkDPoint t1[3], t2[3];
    298     int calcMask = ~0;
    299     do {
    300         if (calcMask & (1 << 1)) t1[1] = quad1.ptAtT(*t1Seed);
    301         if (calcMask & (1 << 4)) t2[1] = quad2.ptAtT(*t2Seed);
    302         if (t1[1].approximatelyEqual(t2[1])) {
    303             *pt = t1[1];
    304     #if ONE_OFF_DEBUG
    305             SkDebugf("%s t1=%1.9g t2=%1.9g (%1.9g,%1.9g) == (%1.9g,%1.9g)\n", __FUNCTION__,
    306                     t1Seed, t2Seed, t1[1].fX, t1[1].fY, t2[1].fX, t2[1].fY);
    307     #endif
    308             return true;
    309         }
    310         if (calcMask & (1 << 0)) t1[0] = quad1.ptAtT(SkTMax(0., *t1Seed - tStep));
    311         if (calcMask & (1 << 2)) t1[2] = quad1.ptAtT(SkTMin(1., *t1Seed + tStep));
    312         if (calcMask & (1 << 3)) t2[0] = quad2.ptAtT(SkTMax(0., *t2Seed - tStep));
    313         if (calcMask & (1 << 5)) t2[2] = quad2.ptAtT(SkTMin(1., *t2Seed + tStep));
    314         double dist[3][3];
    315         // OPTIMIZE: using calcMask value permits skipping some distance calcuations
    316         //   if prior loop's results are moved to correct slot for reuse
    317         dist[1][1] = t1[1].distanceSquared(t2[1]);
    318         int best_i = 1, best_j = 1;
    319         for (int i = 0; i < 3; ++i) {
    320             for (int j = 0; j < 3; ++j) {
    321                 if (i == 1 && j == 1) {
    322                     continue;
    323                 }
    324                 dist[i][j] = t1[i].distanceSquared(t2[j]);
    325                 if (dist[best_i][best_j] > dist[i][j]) {
    326                     best_i = i;
    327                     best_j = j;
    328                 }
    329             }
    330         }
    331         if (best_i == 1 && best_j == 1) {
    332             tStep /= 2;
    333             if (tStep < FLT_EPSILON_HALF) {
    334                 break;
    335             }
    336             calcMask = (1 << 0) | (1 << 2) | (1 << 3) | (1 << 5);
    337             continue;
    338         }
    339         if (best_i == 0) {
    340             *t1Seed -= tStep;
    341             t1[2] = t1[1];
    342             t1[1] = t1[0];
    343             calcMask = 1 << 0;
    344         } else if (best_i == 2) {
    345             *t1Seed += tStep;
    346             t1[0] = t1[1];
    347             t1[1] = t1[2];
    348             calcMask = 1 << 2;
    349         } else {
    350             calcMask = 0;
    351         }
    352         if (best_j == 0) {
    353             *t2Seed -= tStep;
    354             t2[2] = t2[1];
    355             t2[1] = t2[0];
    356             calcMask |= 1 << 3;
    357         } else if (best_j == 2) {
    358             *t2Seed += tStep;
    359             t2[0] = t2[1];
    360             t2[1] = t2[2];
    361             calcMask |= 1 << 5;
    362         }
    363     } while (true);
    364 #if ONE_OFF_DEBUG
    365     SkDebugf("%s t1=%1.9g t2=%1.9g (%1.9g,%1.9g) != (%1.9g,%1.9g) %s\n", __FUNCTION__,
    366         t1Seed, t2Seed, t1[1].fX, t1[1].fY, t1[2].fX, t1[2].fY);
    367 #endif
    368     return false;
    369 }
    370 
    371 static void lookNearEnd(const SkDQuad& q1, const SkDQuad& q2, int testT,
    372         const SkIntersections& orig, bool swap, SkIntersections* i) {
    373     if (orig.used() == 1 && orig[!swap][0] == testT) {
    374         return;
    375     }
    376     if (orig.used() == 2 && orig[!swap][1] == testT) {
    377         return;
    378     }
    379     SkDLine tmpLine;
    380     int testTIndex = testT << 1;
    381     tmpLine[0] = tmpLine[1] = q2[testTIndex];
    382     tmpLine[1].fX += q2[1].fY - q2[testTIndex].fY;
    383     tmpLine[1].fY -= q2[1].fX - q2[testTIndex].fX;
    384     SkIntersections impTs;
    385     impTs.intersectRay(q1, tmpLine);
    386     for (int index = 0; index < impTs.used(); ++index) {
    387         SkDPoint realPt = impTs.pt(index);
    388         if (!tmpLine[0].approximatelyPEqual(realPt)) {
    389             continue;
    390         }
    391         if (swap) {
    392             i->insert(testT, impTs[0][index], tmpLine[0]);
    393         } else {
    394             i->insert(impTs[0][index], testT, tmpLine[0]);
    395         }
    396     }
    397 }
    398 
    399 int SkIntersections::intersect(const SkDQuad& q1, const SkDQuad& q2) {
    400     fMax = 4;
    401     // if the quads share an end point, check to see if they overlap
    402     for (int i1 = 0; i1 < 3; i1 += 2) {
    403         for (int i2 = 0; i2 < 3; i2 += 2) {
    404             if (q1[i1].asSkPoint() == q2[i2].asSkPoint()) {
    405                 insert(i1 >> 1, i2 >> 1, q1[i1]);
    406             }
    407         }
    408     }
    409     SkASSERT(fUsed < 3);
    410     if (only_end_pts_in_common(q1, q2)) {
    411         return fUsed;
    412     }
    413     if (only_end_pts_in_common(q2, q1)) {
    414         return fUsed;
    415     }
    416     // see if either quad is really a line
    417     // FIXME: figure out why reduce step didn't find this earlier
    418     if (is_linear(q1, q2, this)) {
    419         return fUsed;
    420     }
    421     SkIntersections swapped;
    422     swapped.setMax(fMax);
    423     if (is_linear(q2, q1, &swapped)) {
    424         swapped.swapPts();
    425         *this = swapped;
    426         return fUsed;
    427     }
    428     SkIntersections copyI(*this);
    429     lookNearEnd(q1, q2, 0, *this, false, &copyI);
    430     lookNearEnd(q1, q2, 1, *this, false, &copyI);
    431     lookNearEnd(q2, q1, 0, *this, true, &copyI);
    432     lookNearEnd(q2, q1, 1, *this, true, &copyI);
    433     int innerEqual = 0;
    434     if (copyI.fUsed >= 2) {
    435         SkASSERT(copyI.fUsed <= 4);
    436         double width = copyI[0][1] - copyI[0][0];
    437         int midEnd = 1;
    438         for (int index = 2; index < copyI.fUsed; ++index) {
    439             double testWidth = copyI[0][index] - copyI[0][index - 1];
    440             if (testWidth <= width) {
    441                 continue;
    442             }
    443             midEnd = index;
    444         }
    445         for (int index = 0; index < 2; ++index) {
    446             double testT = (copyI[0][midEnd] * (index + 1)
    447                     + copyI[0][midEnd - 1] * (2 - index)) / 3;
    448             SkDPoint testPt1 = q1.ptAtT(testT);
    449             testT = (copyI[1][midEnd] * (index + 1) + copyI[1][midEnd - 1] * (2 - index)) / 3;
    450             SkDPoint testPt2 = q2.ptAtT(testT);
    451             innerEqual += testPt1.approximatelyEqual(testPt2);
    452         }
    453     }
    454     bool expectCoincident = copyI.fUsed >= 2 && innerEqual == 2;
    455     if (expectCoincident) {
    456         reset();
    457         insertCoincident(copyI[0][0], copyI[1][0], copyI.fPt[0]);
    458         int last = copyI.fUsed - 1;
    459         insertCoincident(copyI[0][last], copyI[1][last], copyI.fPt[last]);
    460         return fUsed;
    461     }
    462     SkDQuadImplicit i1(q1);
    463     SkDQuadImplicit i2(q2);
    464     int index;
    465     bool flip1 = q1[2] == q2[0];
    466     bool flip2 = q1[0] == q2[2];
    467     bool useCubic = q1[0] == q2[0];
    468     double roots1[4];
    469     int rootCount = findRoots(i2, q1, roots1, useCubic, flip1, 0);
    470     // OPTIMIZATION: could short circuit here if all roots are < 0 or > 1
    471     double roots1Copy[4];
    472     int r1Count = addValidRoots(roots1, rootCount, roots1Copy);
    473     SkDPoint pts1[4];
    474     for (index = 0; index < r1Count; ++index) {
    475         pts1[index] = q1.ptAtT(roots1Copy[index]);
    476     }
    477     double roots2[4];
    478     int rootCount2 = findRoots(i1, q2, roots2, useCubic, flip2, 0);
    479     double roots2Copy[4];
    480     int r2Count = addValidRoots(roots2, rootCount2, roots2Copy);
    481     SkDPoint pts2[4];
    482     for (index = 0; index < r2Count; ++index) {
    483         pts2[index] = q2.ptAtT(roots2Copy[index]);
    484     }
    485     if (r1Count == r2Count && r1Count <= 1) {
    486         if (r1Count == 1 && used() == 0) {
    487             if (pts1[0].approximatelyEqual(pts2[0])) {
    488                 insert(roots1Copy[0], roots2Copy[0], pts1[0]);
    489             } else {
    490                 // find intersection by chasing t
    491                 if (binary_search(q1, q2, roots1Copy, roots2Copy, pts1)) {
    492                     insert(roots1Copy[0], roots2Copy[0], pts1[0]);
    493                 }
    494             }
    495         }
    496         return fUsed;
    497     }
    498     int closest[4];
    499     double dist[4];
    500     bool foundSomething = false;
    501     for (index = 0; index < r1Count; ++index) {
    502         dist[index] = DBL_MAX;
    503         closest[index] = -1;
    504         for (int ndex2 = 0; ndex2 < r2Count; ++ndex2) {
    505             if (!pts2[ndex2].approximatelyEqual(pts1[index])) {
    506                 continue;
    507             }
    508             double dx = pts2[ndex2].fX - pts1[index].fX;
    509             double dy = pts2[ndex2].fY - pts1[index].fY;
    510             double distance = dx * dx + dy * dy;
    511             if (dist[index] <= distance) {
    512                 continue;
    513             }
    514             for (int outer = 0; outer < index; ++outer) {
    515                 if (closest[outer] != ndex2) {
    516                     continue;
    517                 }
    518                 if (dist[outer] < distance) {
    519                     goto next;
    520                 }
    521                 closest[outer] = -1;
    522             }
    523             dist[index] = distance;
    524             closest[index] = ndex2;
    525             foundSomething = true;
    526         next:
    527             ;
    528         }
    529     }
    530     if (r1Count && r2Count && !foundSomething) {
    531         relaxed_is_linear(&q1, 0, 1, &q2, 0, 1, this);
    532         return fUsed;
    533     }
    534     int used = 0;
    535     do {
    536         double lowest = DBL_MAX;
    537         int lowestIndex = -1;
    538         for (index = 0; index < r1Count; ++index) {
    539             if (closest[index] < 0) {
    540                 continue;
    541             }
    542             if (roots1Copy[index] < lowest) {
    543                 lowestIndex = index;
    544                 lowest = roots1Copy[index];
    545             }
    546         }
    547         if (lowestIndex < 0) {
    548             break;
    549         }
    550         insert(roots1Copy[lowestIndex], roots2Copy[closest[lowestIndex]],
    551                 pts1[lowestIndex]);
    552         closest[lowestIndex] = -1;
    553     } while (++used < r1Count);
    554     return fUsed;
    555 }
    556