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 "SkIntersections.h" |
| 8 #include "SkPathOpsLine.h" |
| 9 #include "SkPathOpsQuad.h" |
| 10 |
| 11 /* |
| 12 Find the interection of a line and quadratic by solving for valid t values. |
| 13 |
| 14 From http://stackoverflow.com/questions/1853637/how-to-find-the-mathematical-fun
ction-defining-a-bezier-curve |
| 15 |
| 16 "A Bezier curve is a parametric function. A quadratic Bezier curve (i.e. three |
| 17 control points) can be expressed as: F(t) = A(1 - t)^2 + B(1 - t)t + Ct^2 where |
| 18 A, B and C are points and t goes from zero to one. |
| 19 |
| 20 This will give you two equations: |
| 21 |
| 22 x = a(1 - t)^2 + b(1 - t)t + ct^2 |
| 23 y = d(1 - t)^2 + e(1 - t)t + ft^2 |
| 24 |
| 25 If you add for instance the line equation (y = kx + m) to that, you'll end up |
| 26 with three equations and three unknowns (x, y and t)." |
| 27 |
| 28 Similar to above, the quadratic is represented as |
| 29 x = a(1-t)^2 + 2b(1-t)t + ct^2 |
| 30 y = d(1-t)^2 + 2e(1-t)t + ft^2 |
| 31 and the line as |
| 32 y = g*x + h |
| 33 |
| 34 Using Mathematica, solve for the values of t where the quadratic intersects the |
| 35 line: |
| 36 |
| 37 (in) t1 = Resultant[a*(1 - t)^2 + 2*b*(1 - t)*t + c*t^2 - x, |
| 38 d*(1 - t)^2 + 2*e*(1 - t)*t + f*t^2 - g*x - h, x] |
| 39 (out) -d + h + 2 d t - 2 e t - d t^2 + 2 e t^2 - f t^2 + |
| 40 g (a - 2 a t + 2 b t + a t^2 - 2 b t^2 + c t^2) |
| 41 (in) Solve[t1 == 0, t] |
| 42 (out) { |
| 43 {t -> (-2 d + 2 e + 2 a g - 2 b g - |
| 44 Sqrt[(2 d - 2 e - 2 a g + 2 b g)^2 - |
| 45 4 (-d + 2 e - f + a g - 2 b g + c g) (-d + a g + h)]) / |
| 46 (2 (-d + 2 e - f + a g - 2 b g + c g)) |
| 47 }, |
| 48 {t -> (-2 d + 2 e + 2 a g - 2 b g + |
| 49 Sqrt[(2 d - 2 e - 2 a g + 2 b g)^2 - |
| 50 4 (-d + 2 e - f + a g - 2 b g + c g) (-d + a g + h)]) / |
| 51 (2 (-d + 2 e - f + a g - 2 b g + c g)) |
| 52 } |
| 53 } |
| 54 |
| 55 Using the results above (when the line tends towards horizontal) |
| 56 A = (-(d - 2*e + f) + g*(a - 2*b + c) ) |
| 57 B = 2*( (d - e ) - g*(a - b ) ) |
| 58 C = (-(d ) + g*(a ) + h ) |
| 59 |
| 60 If g goes to infinity, we can rewrite the line in terms of x. |
| 61 x = g'*y + h' |
| 62 |
| 63 And solve accordingly in Mathematica: |
| 64 |
| 65 (in) t2 = Resultant[a*(1 - t)^2 + 2*b*(1 - t)*t + c*t^2 - g'*y - h', |
| 66 d*(1 - t)^2 + 2*e*(1 - t)*t + f*t^2 - y, y] |
| 67 (out) a - h' - 2 a t + 2 b t + a t^2 - 2 b t^2 + c t^2 - |
| 68 g' (d - 2 d t + 2 e t + d t^2 - 2 e t^2 + f t^2) |
| 69 (in) Solve[t2 == 0, t] |
| 70 (out) { |
| 71 {t -> (2 a - 2 b - 2 d g' + 2 e g' - |
| 72 Sqrt[(-2 a + 2 b + 2 d g' - 2 e g')^2 - |
| 73 4 (a - 2 b + c - d g' + 2 e g' - f g') (a - d g' - h')]) / |
| 74 (2 (a - 2 b + c - d g' + 2 e g' - f g')) |
| 75 }, |
| 76 {t -> (2 a - 2 b - 2 d g' + 2 e g' + |
| 77 Sqrt[(-2 a + 2 b + 2 d g' - 2 e g')^2 - |
| 78 4 (a - 2 b + c - d g' + 2 e g' - f g') (a - d g' - h')])/ |
| 79 (2 (a - 2 b + c - d g' + 2 e g' - f g')) |
| 80 } |
| 81 } |
| 82 |
| 83 Thus, if the slope of the line tends towards vertical, we use: |
| 84 A = ( (a - 2*b + c) - g'*(d - 2*e + f) ) |
| 85 B = 2*(-(a - b ) + g'*(d - e ) ) |
| 86 C = ( (a ) - g'*(d ) - h' ) |
| 87 */ |
| 88 |
| 89 |
| 90 class LineQuadraticIntersections { |
| 91 public: |
| 92 LineQuadraticIntersections(const SkDQuad& q, const SkDLine& l, SkIntersectio
ns* i) |
| 93 : quad(q) |
| 94 , line(l) |
| 95 , intersections(i) { |
| 96 } |
| 97 |
| 98 int intersectRay(double roots[2]) { |
| 99 /* |
| 100 solve by rotating line+quad so line is horizontal, then finding the root
s |
| 101 set up matrix to rotate quad to x-axis |
| 102 |cos(a) -sin(a)| |
| 103 |sin(a) cos(a)| |
| 104 note that cos(a) = A(djacent) / Hypoteneuse |
| 105 sin(a) = O(pposite) / Hypoteneuse |
| 106 since we are computing Ts, we can ignore hypoteneuse, the scale factor: |
| 107 | A -O | |
| 108 | O A | |
| 109 A = line[1].fX - line[0].fX (adjacent side of the right triangle) |
| 110 O = line[1].fY - line[0].fY (opposite side of the right triangle) |
| 111 for each of the three points (e.g. n = 0 to 2) |
| 112 quad[n].fY' = (quad[n].fY - line[0].fY) * A - (quad[n].fX - line[0].fX)
* O |
| 113 */ |
| 114 double adj = line[1].fX - line[0].fX; |
| 115 double opp = line[1].fY - line[0].fY; |
| 116 double r[3]; |
| 117 for (int n = 0; n < 3; ++n) { |
| 118 r[n] = (quad[n].fY - line[0].fY) * adj - (quad[n].fX - line[0].fX) *
opp; |
| 119 } |
| 120 double A = r[2]; |
| 121 double B = r[1]; |
| 122 double C = r[0]; |
| 123 A += C - 2 * B; // A = a - 2*b + c |
| 124 B -= C; // B = -(b - c) |
| 125 return SkDQuad::RootsValidT(A, 2 * B, C, roots); |
| 126 } |
| 127 |
| 128 int intersect() { |
| 129 addEndPoints(); |
| 130 double rootVals[2]; |
| 131 int roots = intersectRay(rootVals); |
| 132 for (int index = 0; index < roots; ++index) { |
| 133 double quadT = rootVals[index]; |
| 134 double lineT = findLineT(quadT); |
| 135 if (PinTs(&quadT, &lineT)) { |
| 136 SkDPoint pt = line.xyAtT(lineT); |
| 137 intersections->insert(quadT, lineT, pt); |
| 138 } |
| 139 } |
| 140 return intersections->used(); |
| 141 } |
| 142 |
| 143 int horizontalIntersect(double axisIntercept, double roots[2]) { |
| 144 double D = quad[2].fY; // f |
| 145 double E = quad[1].fY; // e |
| 146 double F = quad[0].fY; // d |
| 147 D += F - 2 * E; // D = d - 2*e + f |
| 148 E -= F; // E = -(d - e) |
| 149 F -= axisIntercept; |
| 150 return SkDQuad::RootsValidT(D, 2 * E, F, roots); |
| 151 } |
| 152 |
| 153 int horizontalIntersect(double axisIntercept, double left, double right, boo
l flipped) { |
| 154 addHorizontalEndPoints(left, right, axisIntercept); |
| 155 double rootVals[2]; |
| 156 int roots = horizontalIntersect(axisIntercept, rootVals); |
| 157 for (int index = 0; index < roots; ++index) { |
| 158 double quadT = rootVals[index]; |
| 159 SkDPoint pt = quad.xyAtT(quadT); |
| 160 double lineT = (pt.fX - left) / (right - left); |
| 161 if (PinTs(&quadT, &lineT)) { |
| 162 intersections->insert(quadT, lineT, pt); |
| 163 } |
| 164 } |
| 165 if (flipped) { |
| 166 intersections->flip(); |
| 167 } |
| 168 return intersections->used(); |
| 169 } |
| 170 |
| 171 int verticalIntersect(double axisIntercept, double roots[2]) { |
| 172 double D = quad[2].fX; // f |
| 173 double E = quad[1].fX; // e |
| 174 double F = quad[0].fX; // d |
| 175 D += F - 2 * E; // D = d - 2*e + f |
| 176 E -= F; // E = -(d - e) |
| 177 F -= axisIntercept; |
| 178 return SkDQuad::RootsValidT(D, 2 * E, F, roots); |
| 179 } |
| 180 |
| 181 int verticalIntersect(double axisIntercept, double top, double bottom, bool
flipped) { |
| 182 addVerticalEndPoints(top, bottom, axisIntercept); |
| 183 double rootVals[2]; |
| 184 int roots = verticalIntersect(axisIntercept, rootVals); |
| 185 for (int index = 0; index < roots; ++index) { |
| 186 double quadT = rootVals[index]; |
| 187 SkDPoint pt = quad.xyAtT(quadT); |
| 188 double lineT = (pt.fY - top) / (bottom - top); |
| 189 if (PinTs(&quadT, &lineT)) { |
| 190 intersections->insert(quadT, lineT, pt); |
| 191 } |
| 192 } |
| 193 if (flipped) { |
| 194 intersections->flip(); |
| 195 } |
| 196 return intersections->used(); |
| 197 } |
| 198 |
| 199 protected: |
| 200 // add endpoints first to get zero and one t values exactly |
| 201 void addEndPoints() { |
| 202 for (int qIndex = 0; qIndex < 3; qIndex += 2) { |
| 203 for (int lIndex = 0; lIndex < 2; lIndex++) { |
| 204 if (quad[qIndex] == line[lIndex]) { |
| 205 intersections->insert(qIndex >> 1, lIndex, line[lIndex]); |
| 206 } |
| 207 } |
| 208 } |
| 209 } |
| 210 |
| 211 void addHorizontalEndPoints(double left, double right, double y) { |
| 212 for (int qIndex = 0; qIndex < 3; qIndex += 2) { |
| 213 if (quad[qIndex].fY != y) { |
| 214 continue; |
| 215 } |
| 216 if (quad[qIndex].fX == left) { |
| 217 intersections->insert(qIndex >> 1, 0, quad[qIndex]); |
| 218 } |
| 219 if (quad[qIndex].fX == right) { |
| 220 intersections->insert(qIndex >> 1, 1, quad[qIndex]); |
| 221 } |
| 222 } |
| 223 } |
| 224 |
| 225 void addVerticalEndPoints(double top, double bottom, double x) { |
| 226 for (int qIndex = 0; qIndex < 3; qIndex += 2) { |
| 227 if (quad[qIndex].fX != x) { |
| 228 continue; |
| 229 } |
| 230 if (quad[qIndex].fY == top) { |
| 231 intersections->insert(qIndex >> 1, 0, quad[qIndex]); |
| 232 } |
| 233 if (quad[qIndex].fY == bottom) { |
| 234 intersections->insert(qIndex >> 1, 1, quad[qIndex]); |
| 235 } |
| 236 } |
| 237 } |
| 238 |
| 239 double findLineT(double t) { |
| 240 SkDPoint xy = quad.xyAtT(t); |
| 241 double dx = line[1].fX - line[0].fX; |
| 242 double dy = line[1].fY - line[0].fY; |
| 243 if (fabs(dx) > fabs(dy)) { |
| 244 return (xy.fX - line[0].fX) / dx; |
| 245 } |
| 246 return (xy.fY - line[0].fY) / dy; |
| 247 } |
| 248 |
| 249 static bool PinTs(double* quadT, double* lineT) { |
| 250 if (!approximately_one_or_less(*lineT)) { |
| 251 return false; |
| 252 } |
| 253 if (!approximately_zero_or_more(*lineT)) { |
| 254 return false; |
| 255 } |
| 256 if (precisely_less_than_zero(*quadT)) { |
| 257 *quadT = 0; |
| 258 } else if (precisely_greater_than_one(*quadT)) { |
| 259 *quadT = 1; |
| 260 } |
| 261 if (precisely_less_than_zero(*lineT)) { |
| 262 *lineT = 0; |
| 263 } else if (precisely_greater_than_one(*lineT)) { |
| 264 *lineT = 1; |
| 265 } |
| 266 return true; |
| 267 } |
| 268 |
| 269 private: |
| 270 const SkDQuad& quad; |
| 271 const SkDLine& line; |
| 272 SkIntersections* intersections; |
| 273 }; |
| 274 |
| 275 // utility for pairs of coincident quads |
| 276 static double horizontalIntersect(const SkDQuad& quad, const SkDPoint& pt) { |
| 277 LineQuadraticIntersections q(quad, *(static_cast<SkDLine*>(0)), |
| 278 static_cast<SkIntersections*>(0)); |
| 279 double rootVals[2]; |
| 280 int roots = q.horizontalIntersect(pt.fY, rootVals); |
| 281 for (int index = 0; index < roots; ++index) { |
| 282 double t = rootVals[index]; |
| 283 SkDPoint qPt = quad.xyAtT(t); |
| 284 if (AlmostEqualUlps(qPt.fX, pt.fX)) { |
| 285 return t; |
| 286 } |
| 287 } |
| 288 return -1; |
| 289 } |
| 290 |
| 291 static double verticalIntersect(const SkDQuad& quad, const SkDPoint& pt) { |
| 292 LineQuadraticIntersections q(quad, *(static_cast<SkDLine*>(0)), |
| 293 static_cast<SkIntersections*>(0)); |
| 294 double rootVals[2]; |
| 295 int roots = q.verticalIntersect(pt.fX, rootVals); |
| 296 for (int index = 0; index < roots; ++index) { |
| 297 double t = rootVals[index]; |
| 298 SkDPoint qPt = quad.xyAtT(t); |
| 299 if (AlmostEqualUlps(qPt.fY, pt.fY)) { |
| 300 return t; |
| 301 } |
| 302 } |
| 303 return -1; |
| 304 } |
| 305 |
| 306 double SkIntersections::Axial(const SkDQuad& q1, const SkDPoint& p, bool vertica
l) { |
| 307 if (vertical) { |
| 308 return verticalIntersect(q1, p); |
| 309 } |
| 310 return horizontalIntersect(q1, p); |
| 311 } |
| 312 |
| 313 int SkIntersections::horizontal(const SkDQuad& quad, double left, double right,
double y, |
| 314 bool flipped) { |
| 315 LineQuadraticIntersections q(quad, *(static_cast<SkDLine*>(0)), this); |
| 316 return q.horizontalIntersect(y, left, right, flipped); |
| 317 } |
| 318 |
| 319 int SkIntersections::vertical(const SkDQuad& quad, double top, double bottom, do
uble x, |
| 320 bool flipped) { |
| 321 LineQuadraticIntersections q(quad, *(static_cast<SkDLine*>(0)), this); |
| 322 return q.verticalIntersect(x, top, bottom, flipped); |
| 323 } |
| 324 |
| 325 int SkIntersections::intersect(const SkDQuad& quad, const SkDLine& line) { |
| 326 LineQuadraticIntersections q(quad, line, this); |
| 327 return q.intersect(); |
| 328 } |
| 329 |
| 330 int SkIntersections::intersectRay(const SkDQuad& quad, const SkDLine& line) { |
| 331 LineQuadraticIntersections q(quad, line, this); |
| 332 return q.intersectRay(fT[0]); |
| 333 } |
OLD | NEW |