| OLD | NEW |
|---|
| (Empty) |
| 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 "CurveIntersection.h" | |
| 8 #include "CubicUtilities.h" | |
| 9 | |
| 10 /* from http://tom.cs.byu.edu/~tom/papers/cvgip84.pdf 4.1 | |
| 11 * | |
| 12 * This paper proves that Syvester's method can compute the implicit form of | |
| 13 * the quadratic from the parameterzied form. | |
| 14 * | |
| 15 * Given x = a*t*t*t + b*t*t + c*t + d (the parameterized form) | |
| 16 * y = e*t*t*t + f*t*t + g*t + h | |
| 17 * | |
| 18 * we want to find an equation of the implicit form: | |
| 19 * | |
| 20 * A*x^3 + B*x*x*y + C*x*y*y + D*y^3 + E*x*x + F*x*y + G*y*y + H*x + I*y + J = 0 | |
| 21 * | |
| 22 * The implicit form can be expressed as a 6x6 determinant, as shown. | |
| 23 * | |
| 24 * The resultant obtained by Syvester's method is | |
| 25 * | |
| 26 * | a b c (d - x) 0 0 | | |
| 27 * | 0 a b c (d - x) 0 | | |
| 28 * | 0 0 a b c (d - x) | | |
| 29 * | e f g (h - y) 0 0 | | |
| 30 * | 0 e f g (h - y) 0 | | |
| 31 * | 0 0 e f g (h - y) | | |
| 32 * | |
| 33 * which, according to Mathematica, expands as shown below. | |
| 34 * | |
| 35 * Resultant[a*t^3 + b*t^2 + c*t + d - x, e*t^3 + f*t^2 + g*t + h - y, t] | |
| 36 * | |
| 37 * -d^3 e^3 + c d^2 e^2 f - b d^2 e f^2 + a d^2 f^3 - c^2 d e^2 g + | |
| 38 * 2 b d^2 e^2 g + b c d e f g - 3 a d^2 e f g - a c d f^2 g - | |
| 39 * b^2 d e g^2 + 2 a c d e g^2 + a b d f g^2 - a^2 d g^3 + c^3 e^2 h - | |
| 40 * 3 b c d e^2 h + 3 a d^2 e^2 h - b c^2 e f h + 2 b^2 d e f h + | |
| 41 * a c d e f h + a c^2 f^2 h - 2 a b d f^2 h + b^2 c e g h - | |
| 42 * 2 a c^2 e g h - a b d e g h - a b c f g h + 3 a^2 d f g h + | |
| 43 * a^2 c g^2 h - b^3 e h^2 + 3 a b c e h^2 - 3 a^2 d e h^2 + | |
| 44 * a b^2 f h^2 - 2 a^2 c f h^2 - a^2 b g h^2 + a^3 h^3 + 3 d^2 e^3 x - | |
| 45 * 2 c d e^2 f x + 2 b d e f^2 x - 2 a d f^3 x + c^2 e^2 g x - | |
| 46 * 4 b d e^2 g x - b c e f g x + 6 a d e f g x + a c f^2 g x + | |
| 47 * b^2 e g^2 x - 2 a c e g^2 x - a b f g^2 x + a^2 g^3 x + | |
| 48 * 3 b c e^2 h x - 6 a d e^2 h x - 2 b^2 e f h x - a c e f h x + | |
| 49 * 2 a b f^2 h x + a b e g h x - 3 a^2 f g h x + 3 a^2 e h^2 x - | |
| 50 * 3 d e^3 x^2 + c e^2 f x^2 - b e f^2 x^2 + a f^3 x^2 + | |
| 51 * 2 b e^2 g x^2 - 3 a e f g x^2 + 3 a e^2 h x^2 + e^3 x^3 - | |
| 52 * c^3 e^2 y + 3 b c d e^2 y - 3 a d^2 e^2 y + b c^2 e f y - | |
| 53 * 2 b^2 d e f y - a c d e f y - a c^2 f^2 y + 2 a b d f^2 y - | |
| 54 * b^2 c e g y + 2 a c^2 e g y + a b d e g y + a b c f g y - | |
| 55 * 3 a^2 d f g y - a^2 c g^2 y + 2 b^3 e h y - 6 a b c e h y + | |
| 56 * 6 a^2 d e h y - 2 a b^2 f h y + 4 a^2 c f h y + 2 a^2 b g h y - | |
| 57 * 3 a^3 h^2 y - 3 b c e^2 x y + 6 a d e^2 x y + 2 b^2 e f x y + | |
| 58 * a c e f x y - 2 a b f^2 x y - a b e g x y + 3 a^2 f g x y - | |
| 59 * 6 a^2 e h x y - 3 a e^2 x^2 y - b^3 e y^2 + 3 a b c e y^2 - | |
| 60 * 3 a^2 d e y^2 + a b^2 f y^2 - 2 a^2 c f y^2 - a^2 b g y^2 + | |
| 61 * 3 a^3 h y^2 + 3 a^2 e x y^2 - a^3 y^3 | |
| 62 */ | |
| 63 | |
| 64 enum { | |
| 65 xxx_coeff, // A | |
| 66 xxy_coeff, // B | |
| 67 xyy_coeff, // C | |
| 68 yyy_coeff, // D | |
| 69 xx_coeff, | |
| 70 xy_coeff, | |
| 71 yy_coeff, | |
| 72 x_coeff, | |
| 73 y_coeff, | |
| 74 c_coeff, | |
| 75 coeff_count | |
| 76 }; | |
| 77 | |
| 78 #define USE_SYVESTER 0 // if 0, use control-point base parametric form | |
| 79 #if USE_SYVESTER | |
| 80 | |
| 81 // FIXME: factoring version unwritten | |
| 82 // static bool straight_forward = true; | |
| 83 | |
| 84 /* from CubicParameterizationCode.cpp output: | |
| 85 * double A = e * e * e; | |
| 86 * double B = -3 * a * e * e; | |
| 87 * double C = 3 * a * a * e; | |
| 88 * double D = -a * a * a; | |
| 89 */ | |
| 90 static void calc_ABCD(double a, double e, double p[coeff_count]) { | |
| 91 double ee = e * e; | |
| 92 p[xxx_coeff] = e * ee; | |
| 93 p[xxy_coeff] = -3 * a * ee; | |
| 94 double aa = a * a; | |
| 95 p[xyy_coeff] = 3 * aa * e; | |
| 96 p[yyy_coeff] = -aa * a; | |
| 97 } | |
| 98 | |
| 99 /* CubicParameterizationCode.cpp turns Mathematica output into C. | |
| 100 * Rather than edit the lines below, please edit the code there instead. | |
| 101 */ | |
| 102 // start of generated code | |
| 103 static double calc_xx(double a, double b, double c, double d, | |
| 104 double e, double f, double g, double h) { | |
| 105 return | |
| 106 -3 * d * e * e * e | |
| 107 + c * e * e * f | |
| 108 - b * e * f * f | |
| 109 + a * f * f * f | |
| 110 + 2 * b * e * e * g | |
| 111 - 3 * a * e * f * g | |
| 112 + 3 * a * e * e * h; | |
| 113 } | |
| 114 | |
| 115 static double calc_xy(double a, double b, double c, double d, | |
| 116 double e, double f, double g, double h) { | |
| 117 return | |
| 118 -3 * b * c * e * e | |
| 119 + 6 * a * d * e * e | |
| 120 + 2 * b * b * e * f | |
| 121 + a * c * e * f | |
| 122 - 2 * a * b * f * f | |
| 123 - a * b * e * g | |
| 124 + 3 * a * a * f * g | |
| 125 - 6 * a * a * e * h; | |
| 126 } | |
| 127 | |
| 128 static double calc_yy(double a, double b, double c, double d, | |
| 129 double e, double f, double g, double h) { | |
| 130 return | |
| 131 -b * b * b * e | |
| 132 + 3 * a * b * c * e | |
| 133 - 3 * a * a * d * e | |
| 134 + a * b * b * f | |
| 135 - 2 * a * a * c * f | |
| 136 - a * a * b * g | |
| 137 + 3 * a * a * a * h; | |
| 138 } | |
| 139 | |
| 140 static double calc_x(double a, double b, double c, double d, | |
| 141 double e, double f, double g, double h) { | |
| 142 return | |
| 143 3 * d * d * e * e * e | |
| 144 - 2 * c * d * e * e * f | |
| 145 + 2 * b * d * e * f * f | |
| 146 - 2 * a * d * f * f * f | |
| 147 + c * c * e * e * g | |
| 148 - 4 * b * d * e * e * g | |
| 149 - b * c * e * f * g | |
| 150 + 6 * a * d * e * f * g | |
| 151 + a * c * f * f * g | |
| 152 + b * b * e * g * g | |
| 153 - 2 * a * c * e * g * g | |
| 154 - a * b * f * g * g | |
| 155 + a * a * g * g * g | |
| 156 + 3 * b * c * e * e * h | |
| 157 - 6 * a * d * e * e * h | |
| 158 - 2 * b * b * e * f * h | |
| 159 - a * c * e * f * h | |
| 160 + 2 * a * b * f * f * h | |
| 161 + a * b * e * g * h | |
| 162 - 3 * a * a * f * g * h | |
| 163 + 3 * a * a * e * h * h; | |
| 164 } | |
| 165 | |
| 166 static double calc_y(double a, double b, double c, double d, | |
| 167 double e, double f, double g, double h) { | |
| 168 return | |
| 169 -c * c * c * e * e | |
| 170 + 3 * b * c * d * e * e | |
| 171 - 3 * a * d * d * e * e | |
| 172 + b * c * c * e * f | |
| 173 - 2 * b * b * d * e * f | |
| 174 - a * c * d * e * f | |
| 175 - a * c * c * f * f | |
| 176 + 2 * a * b * d * f * f | |
| 177 - b * b * c * e * g | |
| 178 + 2 * a * c * c * e * g | |
| 179 + a * b * d * e * g | |
| 180 + a * b * c * f * g | |
| 181 - 3 * a * a * d * f * g | |
| 182 - a * a * c * g * g | |
| 183 + 2 * b * b * b * e * h | |
| 184 - 6 * a * b * c * e * h | |
| 185 + 6 * a * a * d * e * h | |
| 186 - 2 * a * b * b * f * h | |
| 187 + 4 * a * a * c * f * h | |
| 188 + 2 * a * a * b * g * h | |
| 189 - 3 * a * a * a * h * h; | |
| 190 } | |
| 191 | |
| 192 static double calc_c(double a, double b, double c, double d, | |
| 193 double e, double f, double g, double h) { | |
| 194 return | |
| 195 -d * d * d * e * e * e | |
| 196 + c * d * d * e * e * f | |
| 197 - b * d * d * e * f * f | |
| 198 + a * d * d * f * f * f | |
| 199 - c * c * d * e * e * g | |
| 200 + 2 * b * d * d * e * e * g | |
| 201 + b * c * d * e * f * g | |
| 202 - 3 * a * d * d * e * f * g | |
| 203 - a * c * d * f * f * g | |
| 204 - b * b * d * e * g * g | |
| 205 + 2 * a * c * d * e * g * g | |
| 206 + a * b * d * f * g * g | |
| 207 - a * a * d * g * g * g | |
| 208 + c * c * c * e * e * h | |
| 209 - 3 * b * c * d * e * e * h | |
| 210 + 3 * a * d * d * e * e * h | |
| 211 - b * c * c * e * f * h | |
| 212 + 2 * b * b * d * e * f * h | |
| 213 + a * c * d * e * f * h | |
| 214 + a * c * c * f * f * h | |
| 215 - 2 * a * b * d * f * f * h | |
| 216 + b * b * c * e * g * h | |
| 217 - 2 * a * c * c * e * g * h | |
| 218 - a * b * d * e * g * h | |
| 219 - a * b * c * f * g * h | |
| 220 + 3 * a * a * d * f * g * h | |
| 221 + a * a * c * g * g * h | |
| 222 - b * b * b * e * h * h | |
| 223 + 3 * a * b * c * e * h * h | |
| 224 - 3 * a * a * d * e * h * h | |
| 225 + a * b * b * f * h * h | |
| 226 - 2 * a * a * c * f * h * h | |
| 227 - a * a * b * g * h * h | |
| 228 + a * a * a * h * h * h; | |
| 229 } | |
| 230 // end of generated code | |
| 231 | |
| 232 #else | |
| 233 | |
| 234 /* more Mathematica generated code. This takes a different tack, starting with | |
| 235 the control-point based parametric formulas. The C code is unoptimized -- | |
| 236 in this form, this is a proof of concept (since the other code didn't work) | |
| 237 */ | |
| 238 static double calc_c(double a, double b, double c, double d, | |
| 239 double e, double f, double g, double h) { | |
| 240 return | |
| 241 d*d*d*e*e*e - 3*d*d*(3*c*e*e*f + 3*b*e*(-3*f*f + 2*e*g) + a*(9*f*f*f - 9*e*f*g +
e*e*h)) - | |
| 242 h*(27*c*c*c*e*e - 27*c*c*(3*b*e*f - 3*a*f*f + 2*a*e*g) + | |
| 243 h*(-27*b*b*b*e + 27*a*b*b*f - 9*a*a*b*g + a*a*a*h) + | |
| 244 9*c*(9*b*b*e*g + a*b*(-9*f*g + 3*e*h) + a*a*(3*g*g - 2*f*h))) + | |
| 245 3*d*(9*c*c*e*e*g + 9*b*b*e*(3*g*g - 2*f*h) + 3*a*b*(-9*f*g*g + 6*f*f*h + e*g*
h) + | |
| 246 a*a*(9*g*g*g - 9*f*g*h + e*h*h) + 3*c*(3*b*e*(-3*f*g + e*h) + a*(9*f*f*g -
6*e*g*g - e*f*h))) | |
| 247 ; | |
| 248 } | |
| 249 | |
| 250 // - Power(e - 3*f + 3*g - h,3)*Power(x,3) | |
| 251 static double calc_xxx(double e3f3gh) { | |
| 252 return -e3f3gh * e3f3gh * e3f3gh; | |
| 253 } | |
| 254 | |
| 255 static double calc_y(double a, double b, double c, double d, | |
| 256 double e, double f, double g, double h) { | |
| 257 return | |
| 258 + 3*(6*b*d*d*e*e - d*d*d*e*e + 18*b*b*d*e*f - 18*b*d*d*e*f - | |
| 259 9*b*d*d*f*f - 54*b*b*d*e*g + 12*b*d*d*e*g - 27*b*b*d*g*g - 18*b*b*b*e*h +
18*b*b*d*e*h + | |
| 260 18*b*b*d*f*h + a*a*a*h*h - 9*b*b*b*h*h + 9*c*c*c*e*(e + 2*h) + | |
| 261 a*a*(-3*b*h*(2*g + h) + d*(-27*g*g + 9*g*h - h*(2*e + h) + 9*f*(g + h))) + | |
| 262 a*(9*b*b*h*(2*f + h) - 3*b*d*(6*f*f - 6*f*(3*g - 2*h) + g*(-9*g + h) + e*(
g + h)) + | |
| 263 d*d*(e*e + 9*f*(3*f - g) + e*(-9*f - 9*g + 2*h))) - | |
| 264 9*c*c*(d*e*(e + 2*g) + 3*b*(f*h + e*(f + h)) + a*(-3*f*f - 6*f*h + 2*(g*h
+ e*(g + h)))) + | |
| 265 3*c*(d*d*e*(e + 2*f) + a*a*(3*g*g + 6*g*h - 2*h*(2*f + h)) + 9*b*b*(g*h +
e*(g + h)) + | |
| 266 a*d*(-9*f*f - 18*f*g + 6*g*g + f*h + e*(f + 12*g + h)) + | |
| 267 b*(d*(-3*e*e + 9*f*g + e*(9*f + 9*g - 6*h)) + 3*a*(h*(2*e - 3*g + h) -
3*f*(g + h))))) // *y | |
| 268 ; | |
| 269 } | |
| 270 | |
| 271 static double calc_yy(double a, double b, double c, double d, | |
| 272 double e, double f, double g, double h) { | |
| 273 return | |
| 274 - 3*(18*c*c*c*e - 18*c*c*d*e + 6*c*d*d*e - d*d*d*e + 3*c*d*d*f - 9*c*c*d*g + a*a
*a*h + 9*c*c*c*h - | |
| 275 9*b*b*b*(e + 2*h) - a*a*(d*(e - 9*f + 18*g - 7*h) + 3*c*(2*f - 6*g + h)) + | |
| 276 a*(-9*c*c*(2*e - 6*f + 2*g - h) + d*d*(-7*e + 18*f - 9*g + h) + 3*c*d*(7*e
- 17*f + 3*g + h)) + | |
| 277 9*b*b*(3*c*(e + g + h) + a*(f + 2*h) - d*(e - 2*(f - 3*g + h))) - | |
| 278 3*b*(-(d*d*(e - 6*f + 2*g)) - 3*c*d*(e + 3*f + 3*g - h) + 9*c*c*(e + f + h
) + a*a*(g + 2*h) + | |
| 279 a*(c*(-3*e + 9*f + 9*g + 3*h) + d*(e + 3*f - 17*g + 7*h)))) // *Power(y
,2) | |
| 280 ; | |
| 281 } | |
| 282 | |
| 283 // + Power(a - 3*b + 3*c - d,3)*Power(y,3) | |
| 284 static double calc_yyy(double a3b3cd) { | |
| 285 return a3b3cd * a3b3cd * a3b3cd; | |
| 286 } | |
| 287 | |
| 288 static double calc_xx(double a, double b, double c, double d, | |
| 289 double e, double f, double g, double h) { | |
| 290 return | |
| 291 // + Power(x,2)* | |
| 292 (-3*(-9*b*e*f*f + 9*a*f*f*f + 6*b*e*e*g - 9*a*e*f*g + 27*b*e*f*g - 27*a*f*f*g +
18*a*e*g*g - 54*b*e*g*g + | |
| 293 27*a*f*g*g + 27*b*f*g*g - 18*a*g*g*g + a*e*e*h - 9*b*e*e*h + 3*a*e*f*h
+ 9*b*e*f*h + 9*a*f*f*h - | |
| 294 18*b*f*f*h - 21*a*e*g*h + 51*b*e*g*h - 9*a*f*g*h - 27*b*f*g*h + 18*a*g*
g*h + 7*a*e*h*h - 18*b*e*h*h - 3*a*f*h*h + | |
| 295 18*b*f*h*h - 6*a*g*h*h - 3*b*g*h*h + a*h*h*h + | |
| 296 3*c*(-9*f*f*(g - 2*h) + 3*g*g*h - f*h*(9*g + 2*h) + e*e*(f - 6*g + 6*h)
+ | |
| 297 e*(9*f*g + 6*g*g - 17*f*h - 3*g*h + 3*h*h)) - | |
| 298 d*(e*e*e + e*e*(-6*f - 3*g + 7*h) - 9*(2*f - g)*(f*f + g*g - f*(g + h))
+ | |
| 299 e*(18*f*f + 9*g*g + 3*g*h + h*h - 3*f*(3*g + 7*h)))) ) | |
| 300 ; | |
| 301 } | |
| 302 | |
| 303 // + Power(x,2)*(3*(a - 3*b + 3*c - d)*Power(e - 3*f + 3*g - h,2)*y) | |
| 304 static double calc_xxy(double a3b3cd, double e3f3gh) { | |
| 305 return 3 * a3b3cd * e3f3gh * e3f3gh; | |
| 306 } | |
| 307 | |
| 308 static double calc_x(double a, double b, double c, double d, | |
| 309 double e, double f, double g, double h) { | |
| 310 return | |
| 311 // + x* | |
| 312 (-3*(27*b*b*e*g*g - 27*a*b*f*g*g + 9*a*a*g*g*g - 18*b*b*e*f*h + 18*a*b*f*f*h + 3
*a*b*e*g*h - | |
| 313 27*b*b*e*g*h - 9*a*a*f*g*h + 27*a*b*f*g*h - 9*a*a*g*g*h + a*a*e*h*h - 9
*a*b*e*h*h + | |
| 314 27*b*b*e*h*h + 6*a*a*f*h*h - 18*a*b*f*h*h - 9*b*b*f*h*h + 3*a*a*g*h*h + | |
| 315 6*a*b*g*h*h - a*a*h*h*h + 9*c*c*(e*e*(g - 3*h) - 3*f*f*h + e*(3*f + 2*g
)*h) + | |
| 316 d*d*(e*e*e - 9*f*f*f + 9*e*f*(f + g) - e*e*(3*f + 6*g + h)) + | |
| 317 d*(-3*c*(-9*f*f*g + e*e*(2*f - 6*g - 3*h) + e*(9*f*g + 6*g*g + f*h)) + | |
| 318 a*(-18*f*f*f - 18*e*g*g + 18*g*g*g - 2*e*e*h + 3*e*g*h + 2*e*h*h + 9
*f*f*(3*g + 2*h) + | |
| 319 3*f*(6*e*g - 9*g*g - e*h - 6*g*h)) - 3*b*(9*f*g*g + e*e*(4*g - 3*
h) - 6*f*f*h - | |
| 320 e*(6*f*f + g*(18*g + h) - 3*f*(3*g + 4*h)))) + | |
| 321 3*c*(3*b*(e*e*h + 3*f*g*h - e*(3*f*g - 6*f*h + 6*g*h + h*h)) + | |
| 322 a*(9*f*f*(g - 2*h) + f*h*(-e + 9*g + 4*h) - 3*(2*g*g*h + e*(2*g*g -
4*g*h + h*h))))) ) | |
| 323 ; | |
| 324 } | |
| 325 | |
| 326 static double calc_xy(double a, double b, double c, double d, | |
| 327 double e, double f, double g, double h) { | |
| 328 return | |
| 329 // + x*3* | |
| 330 (-2*a*d*e*e - 7*d*d*e*e + 15*a*d*e*f + 21*d*d*e*f - 9*a*d*f*f - 18*d*d*f*f - 15*
a*d*e*g - | |
| 331 3*d*d*e*g - 9*a*a*f*g + 9*d*d*f*g + 18*a*a*g*g + 9*a*d*g*g + 2*a*a*e*h
- 2*d*d*e*h + | |
| 332 3*a*a*f*h + 15*a*d*f*h - 21*a*a*g*h - 15*a*d*g*h + 7*a*a*h*h + 2*a*d*h*
h - | |
| 333 9*c*c*(2*e*e + 3*f*f + 3*f*h - 2*g*h + e*(-3*f - 4*g + h)) + | |
| 334 9*b*b*(3*g*g - 3*g*h + 2*h*(-2*f + h) + e*(-2*f + 3*g + h)) + | |
| 335 3*b*(3*c*(e*e + 3*e*(f - 3*g) + (9*f - 3*g - h)*h) + a*(6*f*f + e*g - 9
*f*g - 9*g*g - 5*e*h + 9*f*h + 14*g*h - 7*h*h) + | |
| 336 d*(-e*e + 12*f*f - 27*f*g + e*(-9*f + 20*g - 5*h) + g*(9*g + h))) + | |
| 337 3*c*(a*(-(e*f) - 9*f*f + 27*f*g - 12*g*g + 5*e*h - 20*f*h + 9*g*h + h*h
) + | |
| 338 d*(7*e*e + 9*f*f + 9*f*g - 6*g*g - f*h + e*(-14*f - 9*g + 5*h)))) //
*y | |
| 339 ; | |
| 340 } | |
| 341 | |
| 342 // - x*3*Power(a - 3*b + 3*c - d,2)*(e - 3*f + 3*g - h)*Power(y,2) | |
| 343 static double calc_xyy(double a3b3cd, double e3f3gh) { | |
| 344 return -3 * a3b3cd * a3b3cd * e3f3gh; | |
| 345 } | |
| 346 | |
| 347 #endif | |
| 348 | |
| 349 static double (*calc_proc[])(double a, double b, double c, double d, | |
| 350 double e, double f, double g, double h) = { | |
| 351 calc_xx, calc_xy, calc_yy, calc_x, calc_y, calc_c | |
| 352 }; | |
| 353 | |
| 354 #if USE_SYVESTER | |
| 355 /* Control points to parametric coefficients | |
| 356 s = 1 - t | |
| 357 Attt + 3Btts + 3Ctss + Dsss == | |
| 358 Attt + 3B(1 - t)tt + 3C(1 - t)(t - tt) + D(1 - t)(1 - 2t + tt) == | |
| 359 Attt + 3B(tt - ttt) + 3C(t - tt - tt + ttt) + D(1-2t+tt-t+2tt-ttt) == | |
| 360 Attt + 3Btt - 3Bttt + 3Ct - 6Ctt + 3Cttt + D - 3Dt + 3Dtt - Dttt == | |
| 361 D + (3C - 3D)t + (3B - 6C + 3D)tt + (A - 3B + 3C - D)ttt | |
| 362 a = A - 3*B + 3*C - D | |
| 363 b = 3*B - 6*C + 3*D | |
| 364 c = 3*C - 3*D | |
| 365 d = D | |
| 366 */ | |
| 367 | |
| 368 /* http://www.algorithmist.net/bezier3.html | |
| 369 p = 3 * A | |
| 370 q = 3 * B | |
| 371 r = 3 * C | |
| 372 a = A | |
| 373 b = q - p | |
| 374 c = p - 2 * q + r | |
| 375 d = D - A + q - r | |
| 376 | |
| 377 B(t) = a + t * (b + t * (c + t * d)) | |
| 378 | |
| 379 so | |
| 380 | |
| 381 B(t) = a + t*b + t*t*(c + t*d) | |
| 382 = a + t*b + t*t*c + t*t*t*d | |
| 383 */ | |
| 384 static void set_abcd(const double* cubic, double& a, double& b, double& c, | |
| 385 double& d) { | |
| 386 a = cubic[0]; // a = A | |
| 387 b = 3 * cubic[2]; // b = 3*B (compute rest of b lazily) | |
| 388 c = 3 * cubic[4]; // c = 3*C (compute rest of c lazily) | |
| 389 d = cubic[6]; // d = D | |
| 390 a += -b + c - d; // a = A - 3*B + 3*C - D | |
| 391 } | |
| 392 | |
| 393 static void calc_bc(const double d, double& b, double& c) { | |
| 394 b -= 3 * c; // b = 3*B - 3*C | |
| 395 c -= 3 * d; // c = 3*C - 3*D | |
| 396 b -= c; // b = 3*B - 6*C + 3*D | |
| 397 } | |
| 398 | |
| 399 static void alt_set_abcd(const double* cubic, double& a, double& b, double& c, | |
| 400 double& d) { | |
| 401 a = cubic[0]; | |
| 402 double p = 3 * a; | |
| 403 double q = 3 * cubic[2]; | |
| 404 double r = 3 * cubic[4]; | |
| 405 b = q - p; | |
| 406 c = p - 2 * q + r; | |
| 407 d = cubic[6] - a + q - r; | |
| 408 } | |
| 409 | |
| 410 const bool try_alt = true; | |
| 411 | |
| 412 #else | |
| 413 | |
| 414 static void calc_ABCD(double a, double b, double c, double d, | |
| 415 double e, double f, double g, double h, | |
| 416 double p[coeff_count]) { | |
| 417 double a3b3cd = a - 3 * (b - c) - d; | |
| 418 double e3f3gh = e - 3 * (f - g) - h; | |
| 419 p[xxx_coeff] = calc_xxx(e3f3gh); | |
| 420 p[xxy_coeff] = calc_xxy(a3b3cd, e3f3gh); | |
| 421 p[xyy_coeff] = calc_xyy(a3b3cd, e3f3gh); | |
| 422 p[yyy_coeff] = calc_yyy(a3b3cd); | |
| 423 } | |
| 424 #endif | |
| 425 | |
| 426 bool implicit_matches(const Cubic& one, const Cubic& two) { | |
| 427 double p1[coeff_count]; // a'xxx , b'xxy , c'xyy , d'xx , e'xy , f'yy, etc. | |
| 428 double p2[coeff_count]; | |
| 429 #if USE_SYVESTER | |
| 430 double a1, b1, c1, d1; | |
| 431 if (try_alt) | |
| 432 alt_set_abcd(&one[0].x, a1, b1, c1, d1); | |
| 433 else | |
| 434 set_abcd(&one[0].x, a1, b1, c1, d1); | |
| 435 double e1, f1, g1, h1; | |
| 436 if (try_alt) | |
| 437 alt_set_abcd(&one[0].y, e1, f1, g1, h1); | |
| 438 else | |
| 439 set_abcd(&one[0].y, e1, f1, g1, h1); | |
| 440 calc_ABCD(a1, e1, p1); | |
| 441 double a2, b2, c2, d2; | |
| 442 if (try_alt) | |
| 443 alt_set_abcd(&two[0].x, a2, b2, c2, d2); | |
| 444 else | |
| 445 set_abcd(&two[0].x, a2, b2, c2, d2); | |
| 446 double e2, f2, g2, h2; | |
| 447 if (try_alt) | |
| 448 alt_set_abcd(&two[0].y, e2, f2, g2, h2); | |
| 449 else | |
| 450 set_abcd(&two[0].y, e2, f2, g2, h2); | |
| 451 calc_ABCD(a2, e2, p2); | |
| 452 #else | |
| 453 double a1 = one[0].x; | |
| 454 double b1 = one[1].x; | |
| 455 double c1 = one[2].x; | |
| 456 double d1 = one[3].x; | |
| 457 double e1 = one[0].y; | |
| 458 double f1 = one[1].y; | |
| 459 double g1 = one[2].y; | |
| 460 double h1 = one[3].y; | |
| 461 calc_ABCD(a1, b1, c1, d1, e1, f1, g1, h1, p1); | |
| 462 double a2 = two[0].x; | |
| 463 double b2 = two[1].x; | |
| 464 double c2 = two[2].x; | |
| 465 double d2 = two[3].x; | |
| 466 double e2 = two[0].y; | |
| 467 double f2 = two[1].y; | |
| 468 double g2 = two[2].y; | |
| 469 double h2 = two[3].y; | |
| 470 calc_ABCD(a2, b2, c2, d2, e2, f2, g2, h2, p2); | |
| 471 #endif | |
| 472 int first = 0; | |
| 473 for (int index = 0; index < coeff_count; ++index) { | |
| 474 #if USE_SYVESTER | |
| 475 if (!try_alt && index == xx_coeff) { | |
| 476 calc_bc(d1, b1, c1); | |
| 477 calc_bc(h1, f1, g1); | |
| 478 calc_bc(d2, b2, c2); | |
| 479 calc_bc(h2, f2, g2); | |
| 480 } | |
| 481 #endif | |
| 482 if (index >= xx_coeff) { | |
| 483 int procIndex = index - xx_coeff; | |
| 484 p1[index] = (*calc_proc[procIndex])(a1, b1, c1, d1, e1, f1, g1, h1); | |
| 485 p2[index] = (*calc_proc[procIndex])(a2, b2, c2, d2, e2, f2, g2, h2); | |
| 486 } | |
| 487 if (approximately_zero(p1[index]) || approximately_zero(p2[index])) { | |
| 488 first += first == index; | |
| 489 continue; | |
| 490 } | |
| 491 if (first == index) { | |
| 492 continue; | |
| 493 } | |
| 494 if (!AlmostEqualUlps(p1[index] * p2[first], p1[first] * p2[index])) { | |
| 495 return false; | |
| 496 } | |
| 497 } | |
| 498 return true; | |
| 499 } | |
| 500 | |
| 501 static double tangent(const double* cubic, double t) { | |
| 502 double a, b, c, d; | |
| 503 #if USE_SYVESTER | |
| 504 set_abcd(cubic, a, b, c, d); | |
| 505 calc_bc(d, b, c); | |
| 506 #else | |
| 507 coefficients(cubic, a, b, c, d); | |
| 508 #endif | |
| 509 return 3 * a * t * t + 2 * b * t + c; | |
| 510 } | |
| 511 | |
| 512 void tangent(const Cubic& cubic, double t, _Point& result) { | |
| 513 result.x = tangent(&cubic[0].x, t); | |
| 514 result.y = tangent(&cubic[0].y, t); | |
| 515 } | |
| 516 | |
| 517 // unit test to return and validate parametric coefficients | |
| 518 #include "CubicParameterization_TestUtility.cpp" | |
| OLD | NEW |
|---|
Chromium Code Reviews
Issue 867213004:
remove prototype pathops code (Closed)
Base URL: https://skia.googlesource.com/skia.git@master