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 Cubic 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 "CubicUtilities.h" 49 #include "CurveIntersection.h" 50 #include "LineIntersection.h" 51 #include "TSearch.h" 52 53 const bool AVERAGE_END_POINTS = true; // results in better fitting curves 54 55 #define USE_CUBIC_END_POINTS 1 56 57 static double calcTDiv(const Cubic& cubic, double precision, double start) { 58 const double adjust = sqrt(3) / 36; 59 Cubic sub; 60 const Cubic* cPtr; 61 if (start == 0) { 62 cPtr = &cubic; 63 } else { 64 // OPTIMIZE: special-case half-split ? 65 sub_divide(cubic, start, 1, sub); 66 cPtr = ⊂ 67 } 68 const Cubic& c = *cPtr; 69 double dx = c[3].x - 3 * (c[2].x - c[1].x) - c[0].x; 70 double dy = c[3].y - 3 * (c[2].y - c[1].y) - c[0].y; 71 double dist = sqrt(dx * dx + dy * dy); 72 double tDiv3 = precision / (adjust * dist); 73 double t = cube_root(tDiv3); 74 if (start > 0) { 75 t = start + (1 - start) * t; 76 } 77 return t; 78 } 79 80 void demote_cubic_to_quad(const Cubic& cubic, Quadratic& quad) { 81 quad[0] = cubic[0]; 82 if (AVERAGE_END_POINTS) { 83 const _Point fromC1 = { (3 * cubic[1].x - cubic[0].x) / 2, (3 * cubic[1].y - cubic[0].y) / 2 }; 84 const _Point fromC2 = { (3 * cubic[2].x - cubic[3].x) / 2, (3 * cubic[2].y - cubic[3].y) / 2 }; 85 quad[1].x = (fromC1.x + fromC2.x) / 2; 86 quad[1].y = (fromC1.y + fromC2.y) / 2; 87 } else { 88 lineIntersect((const _Line&) cubic[0], (const _Line&) cubic[2], quad[1]); 89 } 90 quad[2] = cubic[3]; 91 } 92 93 int cubic_to_quadratics(const Cubic& cubic, double precision, SkTDArray<Quadratic>& quadratics) { 94 SkTDArray<double> ts; 95 cubic_to_quadratics(cubic, precision, ts); 96 int tsCount = ts.count(); 97 double t1Start = 0; 98 int order = 0; 99 for (int idx = 0; idx <= tsCount; ++idx) { 100 double t1 = idx < tsCount ? ts[idx] : 1; 101 Cubic part; 102 sub_divide(cubic, t1Start, t1, part); 103 Quadratic q1; 104 demote_cubic_to_quad(part, q1); 105 Quadratic s1; 106 int o1 = reduceOrder(q1, s1, kReduceOrder_TreatAsFill); 107 if (order < o1) { 108 order = o1; 109 } 110 memcpy(quadratics.append(), o1 < 2 ? s1 : q1, sizeof(Quadratic)); 111 t1Start = t1; 112 } 113 return order; 114 } 115 116 static bool addSimpleTs(const Cubic& cubic, double precision, SkTDArray<double>& ts) { 117 double tDiv = calcTDiv(cubic, precision, 0); 118 if (tDiv >= 1) { 119 return true; 120 } 121 if (tDiv >= 0.5) { 122 *ts.append() = 0.5; 123 return true; 124 } 125 return false; 126 } 127 128 static void addTs(const Cubic& cubic, double precision, double start, double end, 129 SkTDArray<double>& ts) { 130 double tDiv = calcTDiv(cubic, precision, 0); 131 double parts = ceil(1.0 / tDiv); 132 for (double index = 0; index < parts; ++index) { 133 double newT = start + (index / parts) * (end - start); 134 if (newT > 0 && newT < 1) { 135 *ts.append() = newT; 136 } 137 } 138 } 139 140 // flavor that returns T values only, deferring computing the quads until they are needed 141 // FIXME: when called from recursive intersect 2, this could take the original cubic 142 // and do a more precise job when calling chop at and sub divide by computing the fractional ts. 143 // it would still take the prechopped cubic for reduce order and find cubic inflections 144 void cubic_to_quadratics(const Cubic& cubic, double precision, SkTDArray<double>& ts) { 145 Cubic reduced; 146 int order = reduceOrder(cubic, reduced, kReduceOrder_QuadraticsAllowed, 147 kReduceOrder_TreatAsFill); 148 if (order < 3) { 149 return; 150 } 151 double inflectT[5]; 152 int inflections = find_cubic_inflections(cubic, inflectT); 153 SkASSERT(inflections <= 2); 154 if (!ends_are_extrema_in_x_or_y(cubic)) { 155 inflections += find_cubic_max_curvature(cubic, &inflectT[inflections]); 156 SkASSERT(inflections <= 5); 157 } 158 QSort<double>(inflectT, &inflectT[inflections - 1]); 159 // OPTIMIZATION: is this filtering common enough that it needs to be pulled out into its 160 // own subroutine? 161 while (inflections && approximately_less_than_zero(inflectT[0])) { 162 memcpy(inflectT, &inflectT[1], sizeof(inflectT[0]) * --inflections); 163 } 164 int start = 0; 165 do { 166 int next = start + 1; 167 if (next >= inflections) { 168 break; 169 } 170 if (!approximately_equal(inflectT[start], inflectT[next])) { 171 ++start; 172 continue; 173 } 174 memcpy(&inflectT[start], &inflectT[next], sizeof(inflectT[0]) * (--inflections - start)); 175 } while (true); 176 while (inflections && approximately_greater_than_one(inflectT[inflections - 1])) { 177 --inflections; 178 } 179 CubicPair pair; 180 if (inflections == 1) { 181 chop_at(cubic, pair, inflectT[0]); 182 int orderP1 = reduceOrder(pair.first(), reduced, kReduceOrder_NoQuadraticsAllowed, 183 kReduceOrder_TreatAsFill); 184 if (orderP1 < 2) { 185 --inflections; 186 } else { 187 int orderP2 = reduceOrder(pair.second(), reduced, kReduceOrder_NoQuadraticsAllowed, 188 kReduceOrder_TreatAsFill); 189 if (orderP2 < 2) { 190 --inflections; 191 } 192 } 193 } 194 if (inflections == 0 && addSimpleTs(cubic, precision, ts)) { 195 return; 196 } 197 if (inflections == 1) { 198 chop_at(cubic, pair, inflectT[0]); 199 addTs(pair.first(), precision, 0, inflectT[0], ts); 200 addTs(pair.second(), precision, inflectT[0], 1, ts); 201 return; 202 } 203 if (inflections > 1) { 204 Cubic part; 205 sub_divide(cubic, 0, inflectT[0], part); 206 addTs(part, precision, 0, inflectT[0], ts); 207 int last = inflections - 1; 208 for (int idx = 0; idx < last; ++idx) { 209 sub_divide(cubic, inflectT[idx], inflectT[idx + 1], part); 210 addTs(part, precision, inflectT[idx], inflectT[idx + 1], ts); 211 } 212 sub_divide(cubic, inflectT[last], 1, part); 213 addTs(part, precision, inflectT[last], 1, ts); 214 return; 215 } 216 addTs(cubic, precision, 0, 1, ts); 217 } 218
