caryclark@google.com | 07393ca | 2013-04-08 11:47:37 +0000 | [diff] [blame^] | 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 "SkIntersections.h" |
| 8 | #include "SkPathOpsLine.h" |
| 9 | |
| 10 | /* Determine the intersection point of two lines. This assumes the lines are not parallel, |
| 11 | and that that the lines are infinite. |
| 12 | From http://en.wikipedia.org/wiki/Line-line_intersection |
| 13 | */ |
| 14 | SkDPoint SkIntersections::Line(const SkDLine& a, const SkDLine& b) { |
| 15 | double axLen = a[1].fX - a[0].fX; |
| 16 | double ayLen = a[1].fY - a[0].fY; |
| 17 | double bxLen = b[1].fX - b[0].fX; |
| 18 | double byLen = b[1].fY - b[0].fY; |
| 19 | double denom = byLen * axLen - ayLen * bxLen; |
| 20 | SkASSERT(denom); |
| 21 | double term1 = a[1].fX * a[0].fY - a[1].fY * a[0].fX; |
| 22 | double term2 = b[1].fX * b[0].fY - b[1].fY * b[0].fX; |
| 23 | SkDPoint p; |
| 24 | p.fX = (term1 * bxLen - axLen * term2) / denom; |
| 25 | p.fY = (term1 * byLen - ayLen * term2) / denom; |
| 26 | return p; |
| 27 | } |
| 28 | |
| 29 | int SkIntersections::computePoints(const SkDLine& line, int used) { |
| 30 | fPt[0] = line.xyAtT(fT[0][0]); |
| 31 | if ((fUsed = used) == 2) { |
| 32 | fPt[1] = line.xyAtT(fT[0][1]); |
| 33 | } |
| 34 | return 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 SkIntersections::intersect(const SkDLine& a, const SkDLine& b) { |
| 44 | double axLen = a[1].fX - a[0].fX; |
| 45 | double ayLen = a[1].fY - a[0].fY; |
| 46 | double bxLen = b[1].fX - b[0].fX; |
| 47 | double byLen = b[1].fY - b[0].fY; |
| 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].fY - b[0].fY; |
| 56 | double ab0x = a[0].fX - b[0].fX; |
| 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 fUsed = 0; |
| 68 | } |
| 69 | fT[0][0] = numerA; |
| 70 | fT[1][0] = numerB; |
| 71 | fPt[0] = a.xyAtT(numerA); |
| 72 | return computePoints(a, 1); |
| 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].fY - ayLen * a[0].fX, |
| 81 | axLen * b[0].fY - ayLen * b[0].fX)) { |
| 82 | return fUsed = 0; |
| 83 | } |
| 84 | const double* aPtr; |
| 85 | const double* bPtr; |
| 86 | if (fabs(axLen) > fabs(ayLen) || fabs(bxLen) > fabs(byLen)) { |
| 87 | aPtr = &a[0].fX; |
| 88 | bPtr = &b[0].fX; |
| 89 | } else { |
| 90 | aPtr = &a[0].fY; |
| 91 | bPtr = &b[0].fY; |
| 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 fUsed = 0; |
| 104 | } |
| 105 | fT[0][0] = 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 fUsed = 0; |
| 111 | } |
| 112 | fT[0][0] = SkTMax<double>(SkTMin<double>(at0, 1.0), 0.0); |
| 113 | fT[0][1] = SkTMax<double>(SkTMin<double>(at1, 1.0), 0.0); |
| 114 | } |
| 115 | double bDenom = b0 - b1; |
| 116 | if (approximately_zero(bDenom)) { |
| 117 | fT[1][0] = fT[1][1] = 0; |
| 118 | } else { |
| 119 | int bIn = aDenom * bDenom < 0; |
| 120 | fT[1][bIn] = SkTMax<double>(SkTMin<double>((b0 - a0) / bDenom, 1.0), 0.0); |
| 121 | fT[1][!bIn] = SkTMax<double>(SkTMin<double>((b0 - a1) / bDenom, 1.0), 0.0); |
| 122 | } |
| 123 | bool second = fabs(fT[0][0] - fT[0][1]) > FLT_EPSILON; |
| 124 | SkASSERT((fabs(fT[1][0] - fT[1][1]) <= FLT_EPSILON) ^ second); |
| 125 | return computePoints(a, 1 + second); |
| 126 | } |
| 127 | |
| 128 | int SkIntersections::horizontal(const SkDLine& line, double y) { |
| 129 | double min = line[0].fY; |
| 130 | double max = line[1].fY; |
| 131 | if (min > max) { |
| 132 | SkTSwap(min, max); |
| 133 | } |
| 134 | if (min > y || max < y) { |
| 135 | return fUsed = 0; |
| 136 | } |
| 137 | if (AlmostEqualUlps(min, max)) { |
| 138 | fT[0][0] = 0; |
| 139 | fT[0][1] = 1; |
| 140 | return fUsed = 2; |
| 141 | } |
| 142 | fT[0][0] = (y - line[0].fY) / (line[1].fY - line[0].fY); |
| 143 | return fUsed = 1; |
| 144 | } |
| 145 | |
| 146 | // OPTIMIZATION Given: dy = line[1].fY - line[0].fY |
| 147 | // and: xIntercept / (y - line[0].fY) == (line[1].fX - line[0].fX) / dy |
| 148 | // then: xIntercept * dy == (line[1].fX - line[0].fX) * (y - line[0].fY) |
| 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].fX - line[0].fX) * (y - line[0].fY) <= 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 SkIntersections::horizontal(const SkDLine& line, double left, double right, double y) { |
| 155 | int result = horizontal(line, y); |
| 156 | if (result != 1) { |
| 157 | SkASSERT(result == 0); // FIXME: this is incorrect if result == 2 |
| 158 | return result; |
| 159 | } |
| 160 | double xIntercept = line[0].fX + fT[0][0] * (line[1].fX - line[0].fX); |
| 161 | if (xIntercept > right || xIntercept < left) { |
| 162 | return fUsed = 0; |
| 163 | } |
| 164 | return result; |
| 165 | } |
| 166 | |
| 167 | int SkIntersections::horizontal(const SkDLine& line, double left, double right, |
| 168 | double y, bool flipped) { |
| 169 | int result = horizontal(line, y); |
| 170 | switch (result) { |
| 171 | case 0: |
| 172 | break; |
| 173 | case 1: { |
| 174 | double xIntercept = line[0].fX + fT[0][0] * (line[1].fX - line[0].fX); |
| 175 | if (xIntercept > right || xIntercept < left) { |
| 176 | return fUsed = 0; |
| 177 | } |
| 178 | fT[1][0] = (xIntercept - left) / (right - left); |
| 179 | break; |
| 180 | } |
| 181 | case 2: |
| 182 | double a0 = line[0].fX; |
| 183 | double a1 = line[1].fX; |
| 184 | double b0 = flipped ? right : left; |
| 185 | double b1 = flipped ? left : right; |
| 186 | // FIXME: share common code below |
| 187 | double at0 = (a0 - b0) / (a0 - a1); |
| 188 | double at1 = (a0 - b1) / (a0 - a1); |
| 189 | if ((at0 < 0 && at1 < 0) || (at0 > 1 && at1 > 1)) { |
| 190 | return fUsed = 0; |
| 191 | } |
| 192 | fT[0][0] = SkTMax<double>(SkTMin<double>(at0, 1.0), 0.0); |
| 193 | fT[0][1] = SkTMax<double>(SkTMin<double>(at1, 1.0), 0.0); |
| 194 | int bIn = (a0 - a1) * (b0 - b1) < 0; |
| 195 | fT[1][bIn] = SkTMax<double>(SkTMin<double>((b0 - a0) / (b0 - b1), 1.0), 0.0); |
| 196 | fT[1][!bIn] = SkTMax<double>(SkTMin<double>((b0 - a1) / (b0 - b1), 1.0), 0.0); |
| 197 | bool second = fabs(fT[0][0] - fT[0][1]) > FLT_EPSILON; |
| 198 | SkASSERT((fabs(fT[1][0] - fT[1][1]) <= FLT_EPSILON) ^ second); |
| 199 | return computePoints(line, 1 + second); |
| 200 | } |
| 201 | if (flipped) { |
| 202 | // OPTIMIZATION: instead of swapping, pass original line, use [1].fX - [0].fX |
| 203 | for (int index = 0; index < result; ++index) { |
| 204 | fT[1][index] = 1 - fT[1][index]; |
| 205 | } |
| 206 | } |
| 207 | return computePoints(line, result); |
| 208 | } |
| 209 | |
| 210 | int SkIntersections::vertical(const SkDLine& line, double x) { |
| 211 | double min = line[0].fX; |
| 212 | double max = line[1].fX; |
| 213 | if (min > max) { |
| 214 | SkTSwap(min, max); |
| 215 | } |
| 216 | if (min > x || max < x) { |
| 217 | return fUsed = 0; |
| 218 | } |
| 219 | if (AlmostEqualUlps(min, max)) { |
| 220 | fT[0][0] = 0; |
| 221 | fT[0][1] = 1; |
| 222 | return fUsed = 2; |
| 223 | } |
| 224 | fT[0][0] = (x - line[0].fX) / (line[1].fX - line[0].fX); |
| 225 | return fUsed = 1; |
| 226 | } |
| 227 | |
| 228 | int SkIntersections::vertical(const SkDLine& line, double top, double bottom, |
| 229 | double x, bool flipped) { |
| 230 | int result = vertical(line, x); |
| 231 | switch (result) { |
| 232 | case 0: |
| 233 | break; |
| 234 | case 1: { |
| 235 | double yIntercept = line[0].fY + fT[0][0] * (line[1].fY - line[0].fY); |
| 236 | if (yIntercept > bottom || yIntercept < top) { |
| 237 | return fUsed = 0; |
| 238 | } |
| 239 | fT[1][0] = (yIntercept - top) / (bottom - top); |
| 240 | break; |
| 241 | } |
| 242 | case 2: |
| 243 | double a0 = line[0].fY; |
| 244 | double a1 = line[1].fY; |
| 245 | double b0 = flipped ? bottom : top; |
| 246 | double b1 = flipped ? top : bottom; |
| 247 | // FIXME: share common code above |
| 248 | double at0 = (a0 - b0) / (a0 - a1); |
| 249 | double at1 = (a0 - b1) / (a0 - a1); |
| 250 | if ((at0 < 0 && at1 < 0) || (at0 > 1 && at1 > 1)) { |
| 251 | return fUsed = 0; |
| 252 | } |
| 253 | fT[0][0] = SkTMax<double>(SkTMin<double>(at0, 1.0), 0.0); |
| 254 | fT[0][1] = SkTMax<double>(SkTMin<double>(at1, 1.0), 0.0); |
| 255 | int bIn = (a0 - a1) * (b0 - b1) < 0; |
| 256 | fT[1][bIn] = SkTMax<double>(SkTMin<double>((b0 - a0) / (b0 - b1), 1.0), 0.0); |
| 257 | fT[1][!bIn] = SkTMax<double>(SkTMin<double>((b0 - a1) / (b0 - b1), 1.0), 0.0); |
| 258 | bool second = fabs(fT[0][0] - fT[0][1]) > FLT_EPSILON; |
| 259 | SkASSERT((fabs(fT[1][0] - fT[1][1]) <= FLT_EPSILON) ^ second); |
| 260 | return computePoints(line, 1 + second); |
| 261 | break; |
| 262 | } |
| 263 | if (flipped) { |
| 264 | // OPTIMIZATION: instead of swapping, pass original line, use [1].fY - [0].fY |
| 265 | for (int index = 0; index < result; ++index) { |
| 266 | fT[1][index] = 1 - fT[1][index]; |
| 267 | } |
| 268 | } |
| 269 | return computePoints(line, result); |
| 270 | } |
| 271 | |
| 272 | // from http://www.bryceboe.com/wordpress/wp-content/uploads/2006/10/intersect.py |
| 273 | // 4 subs, 2 muls, 1 cmp |
| 274 | static bool ccw(const SkDPoint& A, const SkDPoint& B, const SkDPoint& C) { |
| 275 | return (C.fY - A.fY) * (B.fX - A.fX) > (B.fY - A.fY) * (C.fX - A.fX); |
| 276 | } |
| 277 | |
| 278 | // 16 subs, 8 muls, 6 cmps |
| 279 | bool SkIntersections::Test(const SkDLine& a, const SkDLine& b) { |
| 280 | return ccw(a[0], b[0], b[1]) != ccw(a[1], b[0], b[1]) |
| 281 | && ccw(a[0], a[1], b[0]) != ccw(a[0], a[1], b[1]); |
| 282 | } |