Home | History | Annotate | Download | only in pathops
      1 /*
      2  * Copyright 2012 Google Inc.
      3  *
      4  * Use of this source code is governed by a BSD-style license that can be
      5  * found in the LICENSE file.
      6  */
      7 #include "SkIntersections.h"
      8 #include "SkLineParameters.h"
      9 #include "SkPathOpsCubic.h"
     10 #include "SkPathOpsQuad.h"
     11 #include "SkPathOpsTriangle.h"
     12 
     13 // from http://blog.gludion.com/2009/08/distance-to-quadratic-bezier-curve.html
     14 // (currently only used by testing)
     15 double SkDQuad::nearestT(const SkDPoint& pt) const {
     16     SkDVector pos = fPts[0] - pt;
     17     // search points P of bezier curve with PM.(dP / dt) = 0
     18     // a calculus leads to a 3d degree equation :
     19     SkDVector A = fPts[1] - fPts[0];
     20     SkDVector B = fPts[2] - fPts[1];
     21     B -= A;
     22     double a = B.dot(B);
     23     double b = 3 * A.dot(B);
     24     double c = 2 * A.dot(A) + pos.dot(B);
     25     double d = pos.dot(A);
     26     double ts[3];
     27     int roots = SkDCubic::RootsValidT(a, b, c, d, ts);
     28     double d0 = pt.distanceSquared(fPts[0]);
     29     double d2 = pt.distanceSquared(fPts[2]);
     30     double distMin = SkTMin(d0, d2);
     31     int bestIndex = -1;
     32     for (int index = 0; index < roots; ++index) {
     33         SkDPoint onQuad = ptAtT(ts[index]);
     34         double dist = pt.distanceSquared(onQuad);
     35         if (distMin > dist) {
     36             distMin = dist;
     37             bestIndex = index;
     38         }
     39     }
     40     if (bestIndex >= 0) {
     41         return ts[bestIndex];
     42     }
     43     return d0 < d2 ? 0 : 1;
     44 }
     45 
     46 bool SkDQuad::pointInHull(const SkDPoint& pt) const {
     47     return ((const SkDTriangle&) fPts).contains(pt);
     48 }
     49 
     50 SkDPoint SkDQuad::top(double startT, double endT) const {
     51     SkDQuad sub = subDivide(startT, endT);
     52     SkDPoint topPt = sub[0];
     53     if (topPt.fY > sub[2].fY || (topPt.fY == sub[2].fY && topPt.fX > sub[2].fX)) {
     54         topPt = sub[2];
     55     }
     56     if (!between(sub[0].fY, sub[1].fY, sub[2].fY)) {
     57         double extremeT;
     58         if (FindExtrema(sub[0].fY, sub[1].fY, sub[2].fY, &extremeT)) {
     59             extremeT = startT + (endT - startT) * extremeT;
     60             SkDPoint test = ptAtT(extremeT);
     61             if (topPt.fY > test.fY || (topPt.fY == test.fY && topPt.fX > test.fX)) {
     62                 topPt = test;
     63             }
     64         }
     65     }
     66     return topPt;
     67 }
     68 
     69 int SkDQuad::AddValidTs(double s[], int realRoots, double* t) {
     70     int foundRoots = 0;
     71     for (int index = 0; index < realRoots; ++index) {
     72         double tValue = s[index];
     73         if (approximately_zero_or_more(tValue) && approximately_one_or_less(tValue)) {
     74             if (approximately_less_than_zero(tValue)) {
     75                 tValue = 0;
     76             } else if (approximately_greater_than_one(tValue)) {
     77                 tValue = 1;
     78             }
     79             for (int idx2 = 0; idx2 < foundRoots; ++idx2) {
     80                 if (approximately_equal(t[idx2], tValue)) {
     81                     goto nextRoot;
     82                 }
     83             }
     84             t[foundRoots++] = tValue;
     85         }
     86 nextRoot:
     87         {}
     88     }
     89     return foundRoots;
     90 }
     91 
     92 // note: caller expects multiple results to be sorted smaller first
     93 // note: http://en.wikipedia.org/wiki/Loss_of_significance has an interesting
     94 //  analysis of the quadratic equation, suggesting why the following looks at
     95 //  the sign of B -- and further suggesting that the greatest loss of precision
     96 //  is in b squared less two a c
     97 int SkDQuad::RootsValidT(double A, double B, double C, double t[2]) {
     98     double s[2];
     99     int realRoots = RootsReal(A, B, C, s);
    100     int foundRoots = AddValidTs(s, realRoots, t);
    101     return foundRoots;
    102 }
    103 
    104 /*
    105 Numeric Solutions (5.6) suggests to solve the quadratic by computing
    106        Q = -1/2(B + sgn(B)Sqrt(B^2 - 4 A C))
    107 and using the roots
    108       t1 = Q / A
    109       t2 = C / Q
    110 */
    111 // this does not discard real roots <= 0 or >= 1
    112 int SkDQuad::RootsReal(const double A, const double B, const double C, double s[2]) {
    113     const double p = B / (2 * A);
    114     const double q = C / A;
    115     if (approximately_zero(A) && (approximately_zero_inverse(p) || approximately_zero_inverse(q))) {
    116         if (approximately_zero(B)) {
    117             s[0] = 0;
    118             return C == 0;
    119         }
    120         s[0] = -C / B;
    121         return 1;
    122     }
    123     /* normal form: x^2 + px + q = 0 */
    124     const double p2 = p * p;
    125     if (!AlmostEqualUlps(p2, q) && p2 < q) {
    126         return 0;
    127     }
    128     double sqrt_D = 0;
    129     if (p2 > q) {
    130         sqrt_D = sqrt(p2 - q);
    131     }
    132     s[0] = sqrt_D - p;
    133     s[1] = -sqrt_D - p;
    134     return 1 + !AlmostEqualUlps(s[0], s[1]);
    135 }
    136 
    137 bool SkDQuad::isLinear(int startIndex, int endIndex) const {
    138     SkLineParameters lineParameters;
    139     lineParameters.quadEndPoints(*this, startIndex, endIndex);
    140     // FIXME: maybe it's possible to avoid this and compare non-normalized
    141     lineParameters.normalize();
    142     double distance = lineParameters.controlPtDistance(*this);
    143     return approximately_zero(distance);
    144 }
    145 
    146 SkDCubic SkDQuad::toCubic() const {
    147     SkDCubic cubic;
    148     cubic[0] = fPts[0];
    149     cubic[2] = fPts[1];
    150     cubic[3] = fPts[2];
    151     cubic[1].fX = (cubic[0].fX + cubic[2].fX * 2) / 3;
    152     cubic[1].fY = (cubic[0].fY + cubic[2].fY * 2) / 3;
    153     cubic[2].fX = (cubic[3].fX + cubic[2].fX * 2) / 3;
    154     cubic[2].fY = (cubic[3].fY + cubic[2].fY * 2) / 3;
    155     return cubic;
    156 }
    157 
    158 SkDVector SkDQuad::dxdyAtT(double t) const {
    159     double a = t - 1;
    160     double b = 1 - 2 * t;
    161     double c = t;
    162     SkDVector result = { a * fPts[0].fX + b * fPts[1].fX + c * fPts[2].fX,
    163             a * fPts[0].fY + b * fPts[1].fY + c * fPts[2].fY };
    164     return result;
    165 }
    166 
    167 // OPTIMIZE: assert if caller passes in t == 0 / t == 1 ?
    168 SkDPoint SkDQuad::ptAtT(double t) const {
    169     if (0 == t) {
    170         return fPts[0];
    171     }
    172     if (1 == t) {
    173         return fPts[2];
    174     }
    175     double one_t = 1 - t;
    176     double a = one_t * one_t;
    177     double b = 2 * one_t * t;
    178     double c = t * t;
    179     SkDPoint result = { a * fPts[0].fX + b * fPts[1].fX + c * fPts[2].fX,
    180             a * fPts[0].fY + b * fPts[1].fY + c * fPts[2].fY };
    181     return result;
    182 }
    183 
    184 /*
    185 Given a quadratic q, t1, and t2, find a small quadratic segment.
    186 
    187 The new quadratic is defined by A, B, and C, where
    188  A = c[0]*(1 - t1)*(1 - t1) + 2*c[1]*t1*(1 - t1) + c[2]*t1*t1
    189  C = c[3]*(1 - t1)*(1 - t1) + 2*c[2]*t1*(1 - t1) + c[1]*t1*t1
    190 
    191 To find B, compute the point halfway between t1 and t2:
    192 
    193 q(at (t1 + t2)/2) == D
    194 
    195 Next, compute where D must be if we know the value of B:
    196 
    197 _12 = A/2 + B/2
    198 12_ = B/2 + C/2
    199 123 = A/4 + B/2 + C/4
    200     = D
    201 
    202 Group the known values on one side:
    203 
    204 B   = D*2 - A/2 - C/2
    205 */
    206 
    207 static double interp_quad_coords(const double* src, double t) {
    208     double ab = SkDInterp(src[0], src[2], t);
    209     double bc = SkDInterp(src[2], src[4], t);
    210     double abc = SkDInterp(ab, bc, t);
    211     return abc;
    212 }
    213 
    214 bool SkDQuad::monotonicInY() const {
    215     return between(fPts[0].fY, fPts[1].fY, fPts[2].fY);
    216 }
    217 
    218 SkDQuad SkDQuad::subDivide(double t1, double t2) const {
    219     SkDQuad dst;
    220     double ax = dst[0].fX = interp_quad_coords(&fPts[0].fX, t1);
    221     double ay = dst[0].fY = interp_quad_coords(&fPts[0].fY, t1);
    222     double dx = interp_quad_coords(&fPts[0].fX, (t1 + t2) / 2);
    223     double dy = interp_quad_coords(&fPts[0].fY, (t1 + t2) / 2);
    224     double cx = dst[2].fX = interp_quad_coords(&fPts[0].fX, t2);
    225     double cy = dst[2].fY = interp_quad_coords(&fPts[0].fY, t2);
    226     /* bx = */ dst[1].fX = 2*dx - (ax + cx)/2;
    227     /* by = */ dst[1].fY = 2*dy - (ay + cy)/2;
    228     return dst;
    229 }
    230 
    231 void SkDQuad::align(int endIndex, SkDPoint* dstPt) const {
    232     if (fPts[endIndex].fX == fPts[1].fX) {
    233         dstPt->fX = fPts[endIndex].fX;
    234     }
    235     if (fPts[endIndex].fY == fPts[1].fY) {
    236         dstPt->fY = fPts[endIndex].fY;
    237     }
    238 }
    239 
    240 SkDPoint SkDQuad::subDivide(const SkDPoint& a, const SkDPoint& c, double t1, double t2) const {
    241     SkASSERT(t1 != t2);
    242     SkDPoint b;
    243 #if 0
    244     // this approach assumes that the control point computed directly is accurate enough
    245     double dx = interp_quad_coords(&fPts[0].fX, (t1 + t2) / 2);
    246     double dy = interp_quad_coords(&fPts[0].fY, (t1 + t2) / 2);
    247     b.fX = 2 * dx - (a.fX + c.fX) / 2;
    248     b.fY = 2 * dy - (a.fY + c.fY) / 2;
    249 #else
    250     SkDQuad sub = subDivide(t1, t2);
    251     SkDLine b0 = {{a, sub[1] + (a - sub[0])}};
    252     SkDLine b1 = {{c, sub[1] + (c - sub[2])}};
    253     SkIntersections i;
    254     i.intersectRay(b0, b1);
    255     if (i.used() == 1) {
    256         b = i.pt(0);
    257     } else {
    258         SkASSERT(i.used() == 2 || i.used() == 0);
    259         b = SkDPoint::Mid(b0[1], b1[1]);
    260     }
    261 #endif
    262     if (t1 == 0 || t2 == 0) {
    263         align(0, &b);
    264     }
    265     if (t1 == 1 || t2 == 1) {
    266         align(2, &b);
    267     }
    268     if (precisely_subdivide_equal(b.fX, a.fX)) {
    269         b.fX = a.fX;
    270     } else if (precisely_subdivide_equal(b.fX, c.fX)) {
    271         b.fX = c.fX;
    272     }
    273     if (precisely_subdivide_equal(b.fY, a.fY)) {
    274         b.fY = a.fY;
    275     } else if (precisely_subdivide_equal(b.fY, c.fY)) {
    276         b.fY = c.fY;
    277     }
    278     return b;
    279 }
    280 
    281 /* classic one t subdivision */
    282 static void interp_quad_coords(const double* src, double* dst, double t) {
    283     double ab = SkDInterp(src[0], src[2], t);
    284     double bc = SkDInterp(src[2], src[4], t);
    285     dst[0] = src[0];
    286     dst[2] = ab;
    287     dst[4] = SkDInterp(ab, bc, t);
    288     dst[6] = bc;
    289     dst[8] = src[4];
    290 }
    291 
    292 SkDQuadPair SkDQuad::chopAt(double t) const
    293 {
    294     SkDQuadPair dst;
    295     interp_quad_coords(&fPts[0].fX, &dst.pts[0].fX, t);
    296     interp_quad_coords(&fPts[0].fY, &dst.pts[0].fY, t);
    297     return dst;
    298 }
    299 
    300 static int valid_unit_divide(double numer, double denom, double* ratio)
    301 {
    302     if (numer < 0) {
    303         numer = -numer;
    304         denom = -denom;
    305     }
    306     if (denom == 0 || numer == 0 || numer >= denom) {
    307         return 0;
    308     }
    309     double r = numer / denom;
    310     if (r == 0) {  // catch underflow if numer <<<< denom
    311         return 0;
    312     }
    313     *ratio = r;
    314     return 1;
    315 }
    316 
    317 /** Quad'(t) = At + B, where
    318     A = 2(a - 2b + c)
    319     B = 2(b - a)
    320     Solve for t, only if it fits between 0 < t < 1
    321 */
    322 int SkDQuad::FindExtrema(double a, double b, double c, double tValue[1]) {
    323     /*  At + B == 0
    324         t = -B / A
    325     */
    326     return valid_unit_divide(a - b, a - b - b + c, tValue);
    327 }
    328 
    329 /* Parameterization form, given A*t*t + 2*B*t*(1-t) + C*(1-t)*(1-t)
    330  *
    331  * a = A - 2*B +   C
    332  * b =     2*B - 2*C
    333  * c =             C
    334  */
    335 void SkDQuad::SetABC(const double* quad, double* a, double* b, double* c) {
    336     *a = quad[0];      // a = A
    337     *b = 2 * quad[2];  // b =     2*B
    338     *c = quad[4];      // c =             C
    339     *b -= *c;          // b =     2*B -   C
    340     *a -= *b;          // a = A - 2*B +   C
    341     *b -= *c;          // b =     2*B - 2*C
    342 }
    343