| 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 #include "Intersections.h" | |
| 10 #include "LineUtilities.h" | |
| 11 | |
| 12 /* | |
| 13 Find the interection of a line and cubic by solving for valid t values. | |
| 14 | |
| 15 Analogous to line-quadratic intersection, solve line-cubic intersection by | |
| 16 representing the cubic as: | |
| 17 x = a(1-t)^3 + 2b(1-t)^2t + c(1-t)t^2 + dt^3 | |
| 18 y = e(1-t)^3 + 2f(1-t)^2t + g(1-t)t^2 + ht^3 | |
| 19 and the line as: | |
| 20 y = i*x + j (if the line is more horizontal) | |
| 21 or: | |
| 22 x = i*y + j (if the line is more vertical) | |
| 23 | |
| 24 Then using Mathematica, solve for the values of t where the cubic intersects the | |
| 25 line: | |
| 26 | |
| 27 (in) Resultant[ | |
| 28 a*(1 - t)^3 + 3*b*(1 - t)^2*t + 3*c*(1 - t)*t^2 + d*t^3 - x, | |
| 29 e*(1 - t)^3 + 3*f*(1 - t)^2*t + 3*g*(1 - t)*t^2 + h*t^3 - i*x - j, x] | |
| 30 (out) -e + j + | |
| 31 3 e t - 3 f t - | |
| 32 3 e t^2 + 6 f t^2 - 3 g t^2 + | |
| 33 e t^3 - 3 f t^3 + 3 g t^3 - h t^3 + | |
| 34 i ( a - | |
| 35 3 a t + 3 b t + | |
| 36 3 a t^2 - 6 b t^2 + 3 c t^2 - | |
| 37 a t^3 + 3 b t^3 - 3 c t^3 + d t^3 ) | |
| 38 | |
| 39 if i goes to infinity, we can rewrite the line in terms of x. Mathematica: | |
| 40 | |
| 41 (in) Resultant[ | |
| 42 a*(1 - t)^3 + 3*b*(1 - t)^2*t + 3*c*(1 - t)*t^2 + d*t^3 - i*y - j, | |
| 43 e*(1 - t)^3 + 3*f*(1 - t)^2*t + 3*g*(1 - t)*t^2 + h*t^3 - y, y] | |
| 44 (out) a - j - | |
| 45 3 a t + 3 b t + | |
| 46 3 a t^2 - 6 b t^2 + 3 c t^2 - | |
| 47 a t^3 + 3 b t^3 - 3 c t^3 + d t^3 - | |
| 48 i ( e - | |
| 49 3 e t + 3 f t + | |
| 50 3 e t^2 - 6 f t^2 + 3 g t^2 - | |
| 51 e t^3 + 3 f t^3 - 3 g t^3 + h t^3 ) | |
| 52 | |
| 53 Solving this with Mathematica produces an expression with hundreds of terms; | |
| 54 instead, use Numeric Solutions recipe to solve the cubic. | |
| 55 | |
| 56 The near-horizontal case, in terms of: Ax^3 + Bx^2 + Cx + D == 0 | |
| 57 A = (-(-e + 3*f - 3*g + h) + i*(-a + 3*b - 3*c + d) ) | |
| 58 B = 3*(-( e - 2*f + g ) + i*( a - 2*b + c ) ) | |
| 59 C = 3*(-(-e + f ) + i*(-a + b ) ) | |
| 60 D = (-( e ) + i*( a ) + j ) | |
| 61 | |
| 62 The near-vertical case, in terms of: Ax^3 + Bx^2 + Cx + D == 0 | |
| 63 A = ( (-a + 3*b - 3*c + d) - i*(-e + 3*f - 3*g + h) ) | |
| 64 B = 3*( ( a - 2*b + c ) - i*( e - 2*f + g ) ) | |
| 65 C = 3*( (-a + b ) - i*(-e + f ) ) | |
| 66 D = ( ( a ) - i*( e ) - j ) | |
| 67 | |
| 68 For horizontal lines: | |
| 69 (in) Resultant[ | |
| 70 a*(1 - t)^3 + 3*b*(1 - t)^2*t + 3*c*(1 - t)*t^2 + d*t^3 - j, | |
| 71 e*(1 - t)^3 + 3*f*(1 - t)^2*t + 3*g*(1 - t)*t^2 + h*t^3 - y, y] | |
| 72 (out) e - j - | |
| 73 3 e t + 3 f t + | |
| 74 3 e t^2 - 6 f t^2 + 3 g t^2 - | |
| 75 e t^3 + 3 f t^3 - 3 g t^3 + h t^3 | |
| 76 So the cubic coefficients are: | |
| 77 | |
| 78 */ | |
| 79 | |
| 80 class LineCubicIntersections { | |
| 81 public: | |
| 82 | |
| 83 LineCubicIntersections(const Cubic& c, const _Line& l, Intersections& i) | |
| 84 : cubic(c) | |
| 85 , line(l) | |
| 86 , intersections(i) { | |
| 87 } | |
| 88 | |
| 89 // see parallel routine in line quadratic intersections | |
| 90 int intersectRay(double roots[3]) { | |
| 91 double adj = line[1].x - line[0].x; | |
| 92 double opp = line[1].y - line[0].y; | |
| 93 Cubic r; | |
| 94 for (int n = 0; n < 4; ++n) { | |
| 95 r[n].x = (cubic[n].y - line[0].y) * adj - (cubic[n].x - line[0].x) * opp
; | |
| 96 } | |
| 97 double A, B, C, D; | |
| 98 coefficients(&r[0].x, A, B, C, D); | |
| 99 return cubicRootsValidT(A, B, C, D, roots); | |
| 100 } | |
| 101 | |
| 102 int intersect() { | |
| 103 addEndPoints(); | |
| 104 double rootVals[3]; | |
| 105 int roots = intersectRay(rootVals); | |
| 106 for (int index = 0; index < roots; ++index) { | |
| 107 double cubicT = rootVals[index]; | |
| 108 double lineT = findLineT(cubicT); | |
| 109 if (pinTs(cubicT, lineT)) { | |
| 110 _Point pt; | |
| 111 xy_at_t(line, lineT, pt.x, pt.y); | |
| 112 intersections.insert(cubicT, lineT, pt); | |
| 113 } | |
| 114 } | |
| 115 return intersections.fUsed; | |
| 116 } | |
| 117 | |
| 118 int horizontalIntersect(double axisIntercept, double roots[3]) { | |
| 119 double A, B, C, D; | |
| 120 coefficients(&cubic[0].y, A, B, C, D); | |
| 121 D -= axisIntercept; | |
| 122 return cubicRootsValidT(A, B, C, D, roots); | |
| 123 } | |
| 124 | |
| 125 int horizontalIntersect(double axisIntercept, double left, double right, bool fl
ipped) { | |
| 126 addHorizontalEndPoints(left, right, axisIntercept); | |
| 127 double rootVals[3]; | |
| 128 int roots = horizontalIntersect(axisIntercept, rootVals); | |
| 129 for (int index = 0; index < roots; ++index) { | |
| 130 _Point pt; | |
| 131 double cubicT = rootVals[index]; | |
| 132 xy_at_t(cubic, cubicT, pt.x, pt.y); | |
| 133 double lineT = (pt.x - left) / (right - left); | |
| 134 if (pinTs(cubicT, lineT)) { | |
| 135 intersections.insert(cubicT, lineT, pt); | |
| 136 } | |
| 137 } | |
| 138 if (flipped) { | |
| 139 flip(); | |
| 140 } | |
| 141 return intersections.fUsed; | |
| 142 } | |
| 143 | |
| 144 int verticalIntersect(double axisIntercept, double roots[3]) { | |
| 145 double A, B, C, D; | |
| 146 coefficients(&cubic[0].x, A, B, C, D); | |
| 147 D -= axisIntercept; | |
| 148 return cubicRootsValidT(A, B, C, D, roots); | |
| 149 } | |
| 150 | |
| 151 int verticalIntersect(double axisIntercept, double top, double bottom, bool flip
ped) { | |
| 152 addVerticalEndPoints(top, bottom, axisIntercept); | |
| 153 double rootVals[3]; | |
| 154 int roots = verticalIntersect(axisIntercept, rootVals); | |
| 155 for (int index = 0; index < roots; ++index) { | |
| 156 _Point pt; | |
| 157 double cubicT = rootVals[index]; | |
| 158 xy_at_t(cubic, cubicT, pt.x, pt.y); | |
| 159 double lineT = (pt.y - top) / (bottom - top); | |
| 160 if (pinTs(cubicT, lineT)) { | |
| 161 intersections.insert(cubicT, lineT, pt); | |
| 162 } | |
| 163 } | |
| 164 if (flipped) { | |
| 165 flip(); | |
| 166 } | |
| 167 return intersections.fUsed; | |
| 168 } | |
| 169 | |
| 170 protected: | |
| 171 | |
| 172 void addEndPoints() | |
| 173 { | |
| 174 for (int cIndex = 0; cIndex < 4; cIndex += 3) { | |
| 175 for (int lIndex = 0; lIndex < 2; lIndex++) { | |
| 176 if (cubic[cIndex] == line[lIndex]) { | |
| 177 intersections.insert(cIndex >> 1, lIndex, line[lIndex]); | |
| 178 } | |
| 179 } | |
| 180 } | |
| 181 } | |
| 182 | |
| 183 void addHorizontalEndPoints(double left, double right, double y) | |
| 184 { | |
| 185 for (int cIndex = 0; cIndex < 4; cIndex += 3) { | |
| 186 if (cubic[cIndex].y != y) { | |
| 187 continue; | |
| 188 } | |
| 189 if (cubic[cIndex].x == left) { | |
| 190 intersections.insert(cIndex >> 1, 0, cubic[cIndex]); | |
| 191 } | |
| 192 if (cubic[cIndex].x == right) { | |
| 193 intersections.insert(cIndex >> 1, 1, cubic[cIndex]); | |
| 194 } | |
| 195 } | |
| 196 } | |
| 197 | |
| 198 void addVerticalEndPoints(double top, double bottom, double x) | |
| 199 { | |
| 200 for (int cIndex = 0; cIndex < 4; cIndex += 3) { | |
| 201 if (cubic[cIndex].x != x) { | |
| 202 continue; | |
| 203 } | |
| 204 if (cubic[cIndex].y == top) { | |
| 205 intersections.insert(cIndex >> 1, 0, cubic[cIndex]); | |
| 206 } | |
| 207 if (cubic[cIndex].y == bottom) { | |
| 208 intersections.insert(cIndex >> 1, 1, cubic[cIndex]); | |
| 209 } | |
| 210 } | |
| 211 } | |
| 212 | |
| 213 double findLineT(double t) { | |
| 214 double x, y; | |
| 215 xy_at_t(cubic, t, x, y); | |
| 216 double dx = line[1].x - line[0].x; | |
| 217 double dy = line[1].y - line[0].y; | |
| 218 if (fabs(dx) > fabs(dy)) { | |
| 219 return (x - line[0].x) / dx; | |
| 220 } | |
| 221 return (y - line[0].y) / dy; | |
| 222 } | |
| 223 | |
| 224 void flip() { | |
| 225 // OPTIMIZATION: instead of swapping, pass original line, use [1].y - [0].y | |
| 226 int roots = intersections.fUsed; | |
| 227 for (int index = 0; index < roots; ++index) { | |
| 228 intersections.fT[1][index] = 1 - intersections.fT[1][index]; | |
| 229 } | |
| 230 } | |
| 231 | |
| 232 static bool pinTs(double& cubicT, double& lineT) { | |
| 233 if (!approximately_one_or_less(lineT)) { | |
| 234 return false; | |
| 235 } | |
| 236 if (!approximately_zero_or_more(lineT)) { | |
| 237 return false; | |
| 238 } | |
| 239 if (precisely_less_than_zero(cubicT)) { | |
| 240 cubicT = 0; | |
| 241 } else if (precisely_greater_than_one(cubicT)) { | |
| 242 cubicT = 1; | |
| 243 } | |
| 244 if (precisely_less_than_zero(lineT)) { | |
| 245 lineT = 0; | |
| 246 } else if (precisely_greater_than_one(lineT)) { | |
| 247 lineT = 1; | |
| 248 } | |
| 249 return true; | |
| 250 } | |
| 251 | |
| 252 private: | |
| 253 | |
| 254 const Cubic& cubic; | |
| 255 const _Line& line; | |
| 256 Intersections& intersections; | |
| 257 }; | |
| 258 | |
| 259 int horizontalIntersect(const Cubic& cubic, double left, double right, double y, | |
| 260 double tRange[3]) { | |
| 261 LineCubicIntersections c(cubic, *((_Line*) 0), *((Intersections*) 0)); | |
| 262 double rootVals[3]; | |
| 263 int result = c.horizontalIntersect(y, rootVals); | |
| 264 int tCount = 0; | |
| 265 for (int index = 0; index < result; ++index) { | |
| 266 double x, y; | |
| 267 xy_at_t(cubic, rootVals[index], x, y); | |
| 268 if (x < left || x > right) { | |
| 269 continue; | |
| 270 } | |
| 271 tRange[tCount++] = rootVals[index]; | |
| 272 } | |
| 273 return result; | |
| 274 } | |
| 275 | |
| 276 int horizontalIntersect(const Cubic& cubic, double left, double right, double y, | |
| 277 bool flipped, Intersections& intersections) { | |
| 278 LineCubicIntersections c(cubic, *((_Line*) 0), intersections); | |
| 279 return c.horizontalIntersect(y, left, right, flipped); | |
| 280 } | |
| 281 | |
| 282 int verticalIntersect(const Cubic& cubic, double top, double bottom, double x, | |
| 283 bool flipped, Intersections& intersections) { | |
| 284 LineCubicIntersections c(cubic, *((_Line*) 0), intersections); | |
| 285 return c.verticalIntersect(x, top, bottom, flipped); | |
| 286 } | |
| 287 | |
| 288 int intersect(const Cubic& cubic, const _Line& line, Intersections& i) { | |
| 289 LineCubicIntersections c(cubic, line, i); | |
| 290 return c.intersect(); | |
| 291 } | |
| 292 | |
| 293 int intersectRay(const Cubic& cubic, const _Line& line, Intersections& i) { | |
| 294 LineCubicIntersections c(cubic, line, i); | |
| 295 return c.intersectRay(i.fT[0]); | |
| 296 } | |
| OLD | NEW |