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 flipped) { 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 flipped) { 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 } 297