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