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 "Intersections.h" 9 #include "LineIntersection.h" 10 #include "LineUtilities.h" 11 12 /* Determine the intersection point of two lines. This assumes the lines are not parallel, 13 and that that the lines are infinite. 14 From http://en.wikipedia.org/wiki/Line-line_intersection 15 */ 16 void lineIntersect(const _Line& a, const _Line& b, _Point& p) { 17 double axLen = a[1].x - a[0].x; 18 double ayLen = a[1].y - a[0].y; 19 double bxLen = b[1].x - b[0].x; 20 double byLen = b[1].y - b[0].y; 21 double denom = byLen * axLen - ayLen * bxLen; 22 SkASSERT(denom); 23 double term1 = a[1].x * a[0].y - a[1].y * a[0].x; 24 double term2 = b[1].x * b[0].y - b[1].y * b[0].x; 25 p.x = (term1 * bxLen - axLen * term2) / denom; 26 p.y = (term1 * byLen - ayLen * term2) / denom; 27 } 28 29 static int computePoints(const _Line& a, int used, Intersections& i) { 30 i.fPt[0] = xy_at_t(a, i.fT[0][0]); 31 if ((i.fUsed = used) == 2) { 32 i.fPt[1] = xy_at_t(a, i.fT[0][1]); 33 } 34 return i.fUsed; 35 } 36 37 /* 38 Determine the intersection point of two line segments 39 Return FALSE if the lines don't intersect 40 from: http://paulbourke.net/geometry/lineline2d/ 41 */ 42 43 int intersect(const _Line& a, const _Line& b, Intersections& i) { 44 double axLen = a[1].x - a[0].x; 45 double ayLen = a[1].y - a[0].y; 46 double bxLen = b[1].x - b[0].x; 47 double byLen = b[1].y - b[0].y; 48 /* Slopes match when denom goes to zero: 49 axLen / ayLen == bxLen / byLen 50 (ayLen * byLen) * axLen / ayLen == (ayLen * byLen) * bxLen / byLen 51 byLen * axLen == ayLen * bxLen 52 byLen * axLen - ayLen * bxLen == 0 ( == denom ) 53 */ 54 double denom = byLen * axLen - ayLen * bxLen; 55 double ab0y = a[0].y - b[0].y; 56 double ab0x = a[0].x - b[0].x; 57 double numerA = ab0y * bxLen - byLen * ab0x; 58 double numerB = ab0y * axLen - ayLen * ab0x; 59 bool mayNotOverlap = (numerA < 0 && denom > numerA) || (numerA > 0 && denom < numerA) 60 || (numerB < 0 && denom > numerB) || (numerB > 0 && denom < numerB); 61 numerA /= denom; 62 numerB /= denom; 63 if ((!approximately_zero(denom) || (!approximately_zero_inverse(numerA) 64 && !approximately_zero_inverse(numerB))) && !sk_double_isnan(numerA) 65 && !sk_double_isnan(numerB)) { 66 if (mayNotOverlap) { 67 return 0; 68 } 69 i.fT[0][0] = numerA; 70 i.fT[1][0] = numerB; 71 i.fPt[0] = xy_at_t(a, numerA); 72 return computePoints(a, 1, i); 73 } 74 /* See if the axis intercepts match: 75 ay - ax * ayLen / axLen == by - bx * ayLen / axLen 76 axLen * (ay - ax * ayLen / axLen) == axLen * (by - bx * ayLen / axLen) 77 axLen * ay - ax * ayLen == axLen * by - bx * ayLen 78 */ 79 // FIXME: need to use AlmostEqualUlps variant instead 80 if (!approximately_equal_squared(axLen * a[0].y - ayLen * a[0].x, 81 axLen * b[0].y - ayLen * b[0].x)) { 82 return 0; 83 } 84 const double* aPtr; 85 const double* bPtr; 86 if (fabs(axLen) > fabs(ayLen) || fabs(bxLen) > fabs(byLen)) { 87 aPtr = &a[0].x; 88 bPtr = &b[0].x; 89 } else { 90 aPtr = &a[0].y; 91 bPtr = &b[0].y; 92 } 93 double a0 = aPtr[0]; 94 double a1 = aPtr[2]; 95 double b0 = bPtr[0]; 96 double b1 = bPtr[2]; 97 // OPTIMIZATION: restructure to reject before the divide 98 // e.g., if ((a0 - b0) * (a0 - a1) < 0 || abs(a0 - b0) > abs(a0 - a1)) 99 // (except efficient) 100 double aDenom = a0 - a1; 101 if (approximately_zero(aDenom)) { 102 if (!between(b0, a0, b1)) { 103 return 0; 104 } 105 i.fT[0][0] = i.fT[0][1] = 0; 106 } else { 107 double at0 = (a0 - b0) / aDenom; 108 double at1 = (a0 - b1) / aDenom; 109 if ((at0 < 0 && at1 < 0) || (at0 > 1 && at1 > 1)) { 110 return 0; 111 } 112 i.fT[0][0] = SkTMax(SkTMin(at0, 1.0), 0.0); 113 i.fT[0][1] = SkTMax(SkTMin(at1, 1.0), 0.0); 114 } 115 double bDenom = b0 - b1; 116 if (approximately_zero(bDenom)) { 117 i.fT[1][0] = i.fT[1][1] = 0; 118 } else { 119 int bIn = aDenom * bDenom < 0; 120 i.fT[1][bIn] = SkTMax(SkTMin((b0 - a0) / bDenom, 1.0), 0.0); 121 i.fT[1][!bIn] = SkTMax(SkTMin((b0 - a1) / bDenom, 1.0), 0.0); 122 } 123 bool second = fabs(i.fT[0][0] - i.fT[0][1]) > FLT_EPSILON; 124 SkASSERT((fabs(i.fT[1][0] - i.fT[1][1]) <= FLT_EPSILON) ^ second); 125 return computePoints(a, 1 + second, i); 126 } 127 128 int horizontalIntersect(const _Line& line, double y, double tRange[2]) { 129 double min = line[0].y; 130 double max = line[1].y; 131 if (min > max) { 132 SkTSwap(min, max); 133 } 134 if (min > y || max < y) { 135 return 0; 136 } 137 if (AlmostEqualUlps(min, max)) { 138 tRange[0] = 0; 139 tRange[1] = 1; 140 return 2; 141 } 142 tRange[0] = (y - line[0].y) / (line[1].y - line[0].y); 143 return 1; 144 } 145 146 // OPTIMIZATION Given: dy = line[1].y - line[0].y 147 // and: xIntercept / (y - line[0].y) == (line[1].x - line[0].x) / dy 148 // then: xIntercept * dy == (line[1].x - line[0].x) * (y - line[0].y) 149 // Assuming that dy is always > 0, the line segment intercepts if: 150 // left * dy <= xIntercept * dy <= right * dy 151 // thus: left * dy <= (line[1].x - line[0].x) * (y - line[0].y) <= right * dy 152 // (clever as this is, it does not give us the t value, so may be useful only 153 // as a quick reject -- and maybe not then; it takes 3 muls, 3 adds, 2 cmps) 154 int horizontalLineIntersect(const _Line& line, double left, double right, 155 double y, double tRange[2]) { 156 int result = horizontalIntersect(line, y, tRange); 157 if (result != 1) { 158 // FIXME: this is incorrect if result == 2 159 return result; 160 } 161 double xIntercept = line[0].x + tRange[0] * (line[1].x - line[0].x); 162 if (xIntercept > right || xIntercept < left) { 163 return 0; 164 } 165 return result; 166 } 167 168 int horizontalIntersect(const _Line& line, double left, double right, 169 double y, bool flipped, Intersections& intersections) { 170 int result = horizontalIntersect(line, y, intersections.fT[0]); 171 switch (result) { 172 case 0: 173 break; 174 case 1: { 175 double xIntercept = line[0].x + intersections.fT[0][0] 176 * (line[1].x - line[0].x); 177 if (xIntercept > right || xIntercept < left) { 178 return 0; 179 } 180 intersections.fT[1][0] = (xIntercept - left) / (right - left); 181 break; 182 } 183 case 2: 184 #if 0 // sorting edges fails to preserve original direction 185 double lineL = line[0].x; 186 double lineR = line[1].x; 187 if (lineL > lineR) { 188 SkTSwap(lineL, lineR); 189 } 190 double overlapL = SkTMax(left, lineL); 191 double overlapR = SkTMin(right, lineR); 192 if (overlapL > overlapR) { 193 return 0; 194 } 195 if (overlapL == overlapR) { 196 result = 1; 197 } 198 intersections.fT[0][0] = (overlapL - line[0].x) / (line[1].x - line[0].x); 199 intersections.fT[1][0] = (overlapL - left) / (right - left); 200 if (result > 1) { 201 intersections.fT[0][1] = (overlapR - line[0].x) / (line[1].x - line[0].x); 202 intersections.fT[1][1] = (overlapR - left) / (right - left); 203 } 204 #else 205 double a0 = line[0].x; 206 double a1 = line[1].x; 207 double b0 = flipped ? right : left; 208 double b1 = flipped ? left : right; 209 // FIXME: share common code below 210 double at0 = (a0 - b0) / (a0 - a1); 211 double at1 = (a0 - b1) / (a0 - a1); 212 if ((at0 < 0 && at1 < 0) || (at0 > 1 && at1 > 1)) { 213 return 0; 214 } 215 intersections.fT[0][0] = SkTMax(SkTMin(at0, 1.0), 0.0); 216 intersections.fT[0][1] = SkTMax(SkTMin(at1, 1.0), 0.0); 217 int bIn = (a0 - a1) * (b0 - b1) < 0; 218 intersections.fT[1][bIn] = SkTMax(SkTMin((b0 - a0) / (b0 - b1), 219 1.0), 0.0); 220 intersections.fT[1][!bIn] = SkTMax(SkTMin((b0 - a1) / (b0 - b1), 221 1.0), 0.0); 222 bool second = fabs(intersections.fT[0][0] - intersections.fT[0][1]) 223 > FLT_EPSILON; 224 SkASSERT((fabs(intersections.fT[1][0] - intersections.fT[1][1]) 225 <= FLT_EPSILON) ^ second); 226 return computePoints(line, 1 + second, intersections); 227 #endif 228 break; 229 } 230 if (flipped) { 231 // OPTIMIZATION: instead of swapping, pass original line, use [1].x - [0].x 232 for (int index = 0; index < result; ++index) { 233 intersections.fT[1][index] = 1 - intersections.fT[1][index]; 234 } 235 } 236 return computePoints(line, result, intersections); 237 } 238 239 static int verticalIntersect(const _Line& line, double x, double tRange[2]) { 240 double min = line[0].x; 241 double max = line[1].x; 242 if (min > max) { 243 SkTSwap(min, max); 244 } 245 if (min > x || max < x) { 246 return 0; 247 } 248 if (AlmostEqualUlps(min, max)) { 249 tRange[0] = 0; 250 tRange[1] = 1; 251 return 2; 252 } 253 tRange[0] = (x - line[0].x) / (line[1].x - line[0].x); 254 return 1; 255 } 256 257 int verticalIntersect(const _Line& line, double top, double bottom, 258 double x, bool flipped, Intersections& intersections) { 259 int result = verticalIntersect(line, x, intersections.fT[0]); 260 switch (result) { 261 case 0: 262 break; 263 case 1: { 264 double yIntercept = line[0].y + intersections.fT[0][0] 265 * (line[1].y - line[0].y); 266 if (yIntercept > bottom || yIntercept < top) { 267 return 0; 268 } 269 intersections.fT[1][0] = (yIntercept - top) / (bottom - top); 270 break; 271 } 272 case 2: 273 #if 0 // sorting edges fails to preserve original direction 274 double lineT = line[0].y; 275 double lineB = line[1].y; 276 if (lineT > lineB) { 277 SkTSwap(lineT, lineB); 278 } 279 double overlapT = SkTMax(top, lineT); 280 double overlapB = SkTMin(bottom, lineB); 281 if (overlapT > overlapB) { 282 return 0; 283 } 284 if (overlapT == overlapB) { 285 result = 1; 286 } 287 intersections.fT[0][0] = (overlapT - line[0].y) / (line[1].y - line[0].y); 288 intersections.fT[1][0] = (overlapT - top) / (bottom - top); 289 if (result > 1) { 290 intersections.fT[0][1] = (overlapB - line[0].y) / (line[1].y - line[0].y); 291 intersections.fT[1][1] = (overlapB - top) / (bottom - top); 292 } 293 #else 294 double a0 = line[0].y; 295 double a1 = line[1].y; 296 double b0 = flipped ? bottom : top; 297 double b1 = flipped ? top : bottom; 298 // FIXME: share common code above 299 double at0 = (a0 - b0) / (a0 - a1); 300 double at1 = (a0 - b1) / (a0 - a1); 301 if ((at0 < 0 && at1 < 0) || (at0 > 1 && at1 > 1)) { 302 return 0; 303 } 304 intersections.fT[0][0] = SkTMax(SkTMin(at0, 1.0), 0.0); 305 intersections.fT[0][1] = SkTMax(SkTMin(at1, 1.0), 0.0); 306 int bIn = (a0 - a1) * (b0 - b1) < 0; 307 intersections.fT[1][bIn] = SkTMax(SkTMin((b0 - a0) / (b0 - b1), 308 1.0), 0.0); 309 intersections.fT[1][!bIn] = SkTMax(SkTMin((b0 - a1) / (b0 - b1), 310 1.0), 0.0); 311 bool second = fabs(intersections.fT[0][0] - intersections.fT[0][1]) 312 > FLT_EPSILON; 313 SkASSERT((fabs(intersections.fT[1][0] - intersections.fT[1][1]) 314 <= FLT_EPSILON) ^ second); 315 return computePoints(line, 1 + second, intersections); 316 #endif 317 break; 318 } 319 if (flipped) { 320 // OPTIMIZATION: instead of swapping, pass original line, use [1].y - [0].y 321 for (int index = 0; index < result; ++index) { 322 intersections.fT[1][index] = 1 - intersections.fT[1][index]; 323 } 324 } 325 return computePoints(line, result, intersections); 326 } 327 328 // from http://www.bryceboe.com/wordpress/wp-content/uploads/2006/10/intersect.py 329 // 4 subs, 2 muls, 1 cmp 330 static bool ccw(const _Point& A, const _Point& B, const _Point& C) { 331 return (C.y - A.y) * (B.x - A.x) > (B.y - A.y) * (C.x - A.x); 332 } 333 334 // 16 subs, 8 muls, 6 cmps 335 bool testIntersect(const _Line& a, const _Line& b) { 336 return ccw(a[0], b[0], b[1]) != ccw(a[1], b[0], b[1]) 337 && ccw(a[0], a[1], b[0]) != ccw(a[0], a[1], b[1]); 338 } 339
