Home | History | Annotate | Download | only in pathops
      1 /*
      2 http://stackoverflow.com/questions/2009160/how-do-i-convert-the-2-control-points-of-a-cubic-curve-to-the-single-control-poi
      3 */
      4 
      5 /*
      6 Let's call the control points of the cubic Q0..Q3 and the control points of the quadratic P0..P2.
      7 Then for degree elevation, the equations are:
      8 
      9 Q0 = P0
     10 Q1 = 1/3 P0 + 2/3 P1
     11 Q2 = 2/3 P1 + 1/3 P2
     12 Q3 = P2
     13 In your case you have Q0..Q3 and you're solving for P0..P2. There are two ways to compute P1 from
     14  the equations above:
     15 
     16 P1 = 3/2 Q1 - 1/2 Q0
     17 P1 = 3/2 Q2 - 1/2 Q3
     18 If this is a degree-elevated cubic, then both equations will give the same answer for P1. Since
     19  it's likely not, your best bet is to average them. So,
     20 
     21 P1 = -1/4 Q0 + 3/4 Q1 + 3/4 Q2 - 1/4 Q3
     22 
     23 
     24 SkDCubic defined by: P1/2 - anchor points, C1/C2 control points
     25 |x| is the euclidean norm of x
     26 mid-point approx of cubic: a quad that shares the same anchors with the cubic and has the
     27  control point at C = (3C2 - P2 + 3C1 - P1)/4
     28 
     29 Algorithm
     30 
     31 pick an absolute precision (prec)
     32 Compute the Tdiv as the root of (cubic) equation
     33 sqrt(3)/18  |P2 - 3C2 + 3C1 - P1|/2  Tdiv ^ 3 = prec
     34 if Tdiv < 0.5 divide the cubic at Tdiv. First segment [0..Tdiv] can be approximated with by a
     35  quadratic, with a defect less than prec, by the mid-point approximation.
     36  Repeat from step 2 with the second resulted segment (corresponding to 1-Tdiv)
     37 0.5<=Tdiv<1 - simply divide the cubic in two. The two halves can be approximated by the mid-point
     38  approximation
     39 Tdiv>=1 - the entire cubic can be approximated by the mid-point approximation
     40 
     41 confirmed by (maybe stolen from)
     42 http://www.caffeineowl.com/graphics/2d/vectorial/cubic2quad01.html
     43 // maybe in turn derived from  http://www.cccg.ca/proceedings/2004/36.pdf
     44 // also stored at http://www.cis.usouthal.edu/~hain/general/Publications/Bezier/bezier%20cccg04%20paper.pdf
     45 
     46 */
     47 
     48 #include "SkPathOpsCubic.h"
     49 #include "SkPathOpsLine.h"
     50 #include "SkPathOpsQuad.h"
     51 #include "SkReduceOrder.h"
     52 #include "SkTArray.h"
     53 #include "SkTSort.h"
     54 
     55 #define USE_CUBIC_END_POINTS 1
     56 
     57 static double calc_t_div(const SkDCubic& cubic, double precision, double start) {
     58     const double adjust = sqrt(3.) / 36;
     59     SkDCubic sub;
     60     const SkDCubic* cPtr;
     61     if (start == 0) {
     62         cPtr = &cubic;
     63     } else {
     64         // OPTIMIZE: special-case half-split ?
     65         sub = cubic.subDivide(start, 1);
     66         cPtr = &sub;
     67     }
     68     const SkDCubic& c = *cPtr;
     69     double dx = c[3].fX - 3 * (c[2].fX - c[1].fX) - c[0].fX;
     70     double dy = c[3].fY - 3 * (c[2].fY - c[1].fY) - c[0].fY;
     71     double dist = sqrt(dx * dx + dy * dy);
     72     double tDiv3 = precision / (adjust * dist);
     73     double t = SkDCubeRoot(tDiv3);
     74     if (start > 0) {
     75         t = start + (1 - start) * t;
     76     }
     77     return t;
     78 }
     79 
     80 SkDQuad SkDCubic::toQuad() const {
     81     SkDQuad quad;
     82     quad[0] = fPts[0];
     83     const SkDPoint fromC1 = {(3 * fPts[1].fX - fPts[0].fX) / 2, (3 * fPts[1].fY - fPts[0].fY) / 2};
     84     const SkDPoint fromC2 = {(3 * fPts[2].fX - fPts[3].fX) / 2, (3 * fPts[2].fY - fPts[3].fY) / 2};
     85     quad[1].fX = (fromC1.fX + fromC2.fX) / 2;
     86     quad[1].fY = (fromC1.fY + fromC2.fY) / 2;
     87     quad[2] = fPts[3];
     88     return quad;
     89 }
     90 
     91 static bool add_simple_ts(const SkDCubic& cubic, double precision, SkTArray<double, true>* ts) {
     92     double tDiv = calc_t_div(cubic, precision, 0);
     93     if (tDiv >= 1) {
     94         return true;
     95     }
     96     if (tDiv >= 0.5) {
     97         ts->push_back(0.5);
     98         return true;
     99     }
    100     return false;
    101 }
    102 
    103 static void addTs(const SkDCubic& cubic, double precision, double start, double end,
    104         SkTArray<double, true>* ts) {
    105     double tDiv = calc_t_div(cubic, precision, 0);
    106     double parts = ceil(1.0 / tDiv);
    107     for (double index = 0; index < parts; ++index) {
    108         double newT = start + (index / parts) * (end - start);
    109         if (newT > 0 && newT < 1) {
    110             ts->push_back(newT);
    111         }
    112     }
    113 }
    114 
    115 // flavor that returns T values only, deferring computing the quads until they are needed
    116 // FIXME: when called from recursive intersect 2, this could take the original cubic
    117 // and do a more precise job when calling chop at and sub divide by computing the fractional ts.
    118 // it would still take the prechopped cubic for reduce order and find cubic inflections
    119 void SkDCubic::toQuadraticTs(double precision, SkTArray<double, true>* ts) const {
    120     SkReduceOrder reducer;
    121     int order = reducer.reduce(*this, SkReduceOrder::kAllow_Quadratics, SkReduceOrder::kFill_Style);
    122     if (order < 3) {
    123         return;
    124     }
    125     double inflectT[5];
    126     int inflections = findInflections(inflectT);
    127     SkASSERT(inflections <= 2);
    128     if (!endsAreExtremaInXOrY()) {
    129         inflections += findMaxCurvature(&inflectT[inflections]);
    130         SkASSERT(inflections <= 5);
    131     }
    132     SkTQSort<double>(inflectT, &inflectT[inflections - 1]);
    133     // OPTIMIZATION: is this filtering common enough that it needs to be pulled out into its
    134     // own subroutine?
    135     while (inflections && approximately_less_than_zero(inflectT[0])) {
    136         memmove(inflectT, &inflectT[1], sizeof(inflectT[0]) * --inflections);
    137     }
    138     int start = 0;
    139     do {
    140         int next = start + 1;
    141         if (next >= inflections) {
    142             break;
    143         }
    144         if (!approximately_equal(inflectT[start], inflectT[next])) {
    145             ++start;
    146             continue;
    147         }
    148         memmove(&inflectT[start], &inflectT[next], sizeof(inflectT[0]) * (--inflections - start));
    149     } while (true);
    150     while (inflections && approximately_greater_than_one(inflectT[inflections - 1])) {
    151         --inflections;
    152     }
    153     SkDCubicPair pair;
    154     if (inflections == 1) {
    155         pair = chopAt(inflectT[0]);
    156         int orderP1 = reducer.reduce(pair.first(), SkReduceOrder::kNo_Quadratics,
    157                 SkReduceOrder::kFill_Style);
    158         if (orderP1 < 2) {
    159             --inflections;
    160         } else {
    161             int orderP2 = reducer.reduce(pair.second(), SkReduceOrder::kNo_Quadratics,
    162                     SkReduceOrder::kFill_Style);
    163             if (orderP2 < 2) {
    164                 --inflections;
    165             }
    166         }
    167     }
    168     if (inflections == 0 && add_simple_ts(*this, precision, ts)) {
    169         return;
    170     }
    171     if (inflections == 1) {
    172         pair = chopAt(inflectT[0]);
    173         addTs(pair.first(), precision, 0, inflectT[0], ts);
    174         addTs(pair.second(), precision, inflectT[0], 1, ts);
    175         return;
    176     }
    177     if (inflections > 1) {
    178         SkDCubic part = subDivide(0, inflectT[0]);
    179         addTs(part, precision, 0, inflectT[0], ts);
    180         int last = inflections - 1;
    181         for (int idx = 0; idx < last; ++idx) {
    182             part = subDivide(inflectT[idx], inflectT[idx + 1]);
    183             addTs(part, precision, inflectT[idx], inflectT[idx + 1], ts);
    184         }
    185         part = subDivide(inflectT[last], 1);
    186         addTs(part, precision, inflectT[last], 1, ts);
    187         return;
    188     }
    189     addTs(*this, precision, 0, 1, ts);
    190 }
    191