caryclark@google.com | 235f56a | 2012-09-14 14:19:30 +0000 | [diff] [blame] | 1 | // Another approach is to start with the implicit form of one curve and solve |
| 2 | // (seek implicit coefficients in QuadraticParameter.cpp |
| 3 | // by substituting in the parametric form of the other. |
| 4 | // The downside of this approach is that early rejects are difficult to come by. |
| 5 | // http://planetmath.org/encyclopedia/GaloisTheoreticDerivationOfTheQuarticFormula.html#step |
| 6 | |
| 7 | |
| 8 | #include "CurveIntersection.h" |
| 9 | #include "Intersections.h" |
| 10 | #include "QuadraticParameterization.h" |
| 11 | #include "QuarticRoot.h" |
| 12 | #include "QuadraticUtilities.h" |
| 13 | |
| 14 | /* given the implicit form 0 = Ax^2 + Bxy + Cy^2 + Dx + Ey + F |
| 15 | * and given x = at^2 + bt + c (the parameterized form) |
| 16 | * y = dt^2 + et + f |
skia.committer@gmail.com | 055c7c2 | 2012-09-15 02:01:41 +0000 | [diff] [blame] | 17 | * then |
caryclark@google.com | 235f56a | 2012-09-14 14:19:30 +0000 | [diff] [blame] | 18 | * 0 = A(at^2+bt+c)(at^2+bt+c)+B(at^2+bt+c)(dt^2+et+f)+C(dt^2+et+f)(dt^2+et+f)+D(at^2+bt+c)+E(dt^2+et+f)+F |
| 19 | */ |
| 20 | |
| 21 | static int findRoots(const QuadImplicitForm& i, const Quadratic& q2, double roots[4]) { |
| 22 | double a, b, c; |
| 23 | set_abc(&q2[0].x, a, b, c); |
| 24 | double d, e, f; |
| 25 | set_abc(&q2[0].y, d, e, f); |
| 26 | const double t4 = i.x2() * a * a |
| 27 | + i.xy() * a * d |
| 28 | + i.y2() * d * d; |
| 29 | const double t3 = 2 * i.x2() * a * b |
| 30 | + i.xy() * (a * e + b * d) |
| 31 | + 2 * i.y2() * d * e; |
| 32 | const double t2 = i.x2() * (b * b + 2 * a * c) |
| 33 | + i.xy() * (c * d + b * e + a * f) |
| 34 | + i.y2() * (e * e + 2 * d * f) |
| 35 | + i.x() * a |
| 36 | + i.y() * d; |
| 37 | const double t1 = 2 * i.x2() * b * c |
| 38 | + i.xy() * (c * e + b * f) |
| 39 | + 2 * i.y2() * e * f |
| 40 | + i.x() * b |
| 41 | + i.y() * e; |
| 42 | const double t0 = i.x2() * c * c |
| 43 | + i.xy() * c * f |
| 44 | + i.y2() * f * f |
| 45 | + i.x() * c |
| 46 | + i.y() * f |
| 47 | + i.c(); |
| 48 | return quarticRoots(t4, t3, t2, t1, t0, roots); |
| 49 | } |
| 50 | |
| 51 | static void addValidRoots(const double roots[4], const int count, const int side, Intersections& i) { |
| 52 | int index; |
| 53 | for (index = 0; index < count; ++index) { |
| 54 | if (!approximately_zero_or_more(roots[index]) || !approximately_one_or_less(roots[index])) { |
| 55 | continue; |
| 56 | } |
| 57 | double t = 1 - roots[index]; |
| 58 | if (approximately_less_than_zero(t)) { |
| 59 | t = 0; |
| 60 | } else if (approximately_greater_than_one(t)) { |
| 61 | t = 1; |
| 62 | } |
| 63 | i.insertOne(t, side); |
| 64 | } |
| 65 | } |
| 66 | |
caryclark@google.com | 6aea33f | 2012-10-09 14:11:58 +0000 | [diff] [blame^] | 67 | static bool onlyEndPtsInCommon(const Quadratic& q1, const Quadratic& q2, Intersections& i) { |
| 68 | // the idea here is to see at minimum do a quick reject by rotating all points |
| 69 | // to either side of the line formed by connecting the endpoints |
| 70 | // if the opposite curves points are on the line or on the other side, the |
| 71 | // curves at most intersect at the endpoints |
| 72 | for (int oddMan = 0; oddMan < 3; ++oddMan) { |
| 73 | const _Point* endPt[2]; |
| 74 | for (int opp = 1; opp < 3; ++opp) { |
| 75 | int end = oddMan ^ opp; |
| 76 | if (end == 3) { |
| 77 | end = opp; |
| 78 | } |
| 79 | endPt[opp - 1] = &q1[end]; |
| 80 | } |
| 81 | double origX = endPt[0]->x; |
| 82 | double origY = endPt[0]->y; |
| 83 | double adj = endPt[1]->x - origX; |
| 84 | double opp = endPt[1]->y - origY; |
| 85 | double sign = (q1[oddMan].y - origY) * adj - (q1[oddMan].x - origX) * opp; |
| 86 | assert(!approximately_zero(sign)); |
| 87 | for (int n = 0; n < 3; ++n) { |
| 88 | double test = (q2[n].y - origY) * adj - (q2[n].x - origX) * opp; |
| 89 | if (test * sign > 0) { |
| 90 | goto tryNextHalfPlane; |
| 91 | } |
| 92 | } |
| 93 | for (int i1 = 0; i1 < 3; i1 += 2) { |
| 94 | for (int i2 = 0; i2 < 3; i2 += 2) { |
| 95 | if (q1[i1] == q2[i2]) { |
| 96 | i.insertOne(i1 >> 1, 0); |
| 97 | i.insertOne(i2 >> 1, 1); |
| 98 | } |
| 99 | } |
| 100 | } |
| 101 | assert(i.fUsed < 3); |
| 102 | return true; |
| 103 | tryNextHalfPlane: |
| 104 | ; |
| 105 | } |
| 106 | return false; |
| 107 | } |
| 108 | |
caryclark@google.com | 235f56a | 2012-09-14 14:19:30 +0000 | [diff] [blame] | 109 | bool intersect2(const Quadratic& q1, const Quadratic& q2, Intersections& i) { |
caryclark@google.com | 6aea33f | 2012-10-09 14:11:58 +0000 | [diff] [blame^] | 110 | // if the quads share an end point, check to see if they overlap |
| 111 | |
| 112 | if (onlyEndPtsInCommon(q1, q2, i)) { |
| 113 | assert(i.insertBalanced()); |
| 114 | return i.intersected(); |
| 115 | } |
caryclark@google.com | 235f56a | 2012-09-14 14:19:30 +0000 | [diff] [blame] | 116 | QuadImplicitForm i1(q1); |
| 117 | QuadImplicitForm i2(q2); |
| 118 | if (i1.implicit_match(i2)) { |
| 119 | // FIXME: compute T values |
| 120 | // compute the intersections of the ends to find the coincident span |
| 121 | bool useVertical = fabs(q1[0].x - q1[2].x) < fabs(q1[0].y - q1[2].y); |
| 122 | double t; |
| 123 | if ((t = axialIntersect(q1, q2[0], useVertical)) >= 0) { |
| 124 | i.addCoincident(t, 0); |
| 125 | } |
| 126 | if ((t = axialIntersect(q1, q2[2], useVertical)) >= 0) { |
| 127 | i.addCoincident(t, 1); |
| 128 | } |
| 129 | useVertical = fabs(q2[0].x - q2[2].x) < fabs(q2[0].y - q2[2].y); |
| 130 | if ((t = axialIntersect(q2, q1[0], useVertical)) >= 0) { |
| 131 | i.addCoincident(0, t); |
| 132 | } |
| 133 | if ((t = axialIntersect(q2, q1[2], useVertical)) >= 0) { |
| 134 | i.addCoincident(1, t); |
| 135 | } |
| 136 | assert(i.fCoincidentUsed <= 2); |
| 137 | return i.fCoincidentUsed > 0; |
| 138 | } |
| 139 | double roots1[4], roots2[4]; |
| 140 | int rootCount = findRoots(i2, q1, roots1); |
| 141 | // OPTIMIZATION: could short circuit here if all roots are < 0 or > 1 |
caryclark@google.com | d168874 | 2012-09-18 20:08:37 +0000 | [diff] [blame] | 142 | #ifndef NDEBUG |
skia.committer@gmail.com | c1ad022 | 2012-09-19 02:01:47 +0000 | [diff] [blame] | 143 | int rootCount2 = |
caryclark@google.com | d168874 | 2012-09-18 20:08:37 +0000 | [diff] [blame] | 144 | #endif |
| 145 | findRoots(i1, q2, roots2); |
caryclark@google.com | 235f56a | 2012-09-14 14:19:30 +0000 | [diff] [blame] | 146 | assert(rootCount == rootCount2); |
| 147 | addValidRoots(roots1, rootCount, 0, i); |
| 148 | addValidRoots(roots2, rootCount, 1, i); |
| 149 | _Point pts[4]; |
| 150 | bool matches[4]; |
caryclark@google.com | 6aea33f | 2012-10-09 14:11:58 +0000 | [diff] [blame^] | 151 | int flipCheck[4]; |
caryclark@google.com | d168874 | 2012-09-18 20:08:37 +0000 | [diff] [blame] | 152 | int index, ndex2; |
caryclark@google.com | 6aea33f | 2012-10-09 14:11:58 +0000 | [diff] [blame^] | 153 | int flipIndex = 0; |
caryclark@google.com | d168874 | 2012-09-18 20:08:37 +0000 | [diff] [blame] | 154 | for (ndex2 = 0; ndex2 < i.fUsed2; ++ndex2) { |
| 155 | xy_at_t(q2, i.fT[1][ndex2], pts[ndex2].x, pts[ndex2].y); |
| 156 | matches[ndex2] = false; |
caryclark@google.com | 235f56a | 2012-09-14 14:19:30 +0000 | [diff] [blame] | 157 | } |
| 158 | for (index = 0; index < i.fUsed; ) { |
| 159 | _Point xy; |
| 160 | xy_at_t(q1, i.fT[0][index], xy.x, xy.y); |
caryclark@google.com | d168874 | 2012-09-18 20:08:37 +0000 | [diff] [blame] | 161 | for (ndex2 = 0; ndex2 < i.fUsed2; ++ndex2) { |
| 162 | if (approximately_equal(pts[ndex2].x, xy.x) && approximately_equal(pts[ndex2].y, xy.y)) { |
caryclark@google.com | 6aea33f | 2012-10-09 14:11:58 +0000 | [diff] [blame^] | 163 | assert(flipIndex < 4); |
| 164 | flipCheck[flipIndex++] = ndex2; |
caryclark@google.com | d168874 | 2012-09-18 20:08:37 +0000 | [diff] [blame] | 165 | matches[ndex2] = true; |
caryclark@google.com | 235f56a | 2012-09-14 14:19:30 +0000 | [diff] [blame] | 166 | goto next; |
| 167 | } |
| 168 | } |
| 169 | if (--i.fUsed > index) { |
| 170 | memmove(&i.fT[0][index], &i.fT[0][index + 1], (i.fUsed - index) * sizeof(i.fT[0][0])); |
| 171 | continue; |
| 172 | } |
| 173 | next: |
| 174 | ++index; |
| 175 | } |
caryclark@google.com | d168874 | 2012-09-18 20:08:37 +0000 | [diff] [blame] | 176 | for (ndex2 = 0; ndex2 < i.fUsed2; ) { |
| 177 | if (!matches[ndex2]) { |
| 178 | if (--i.fUsed2 > ndex2) { |
| 179 | memmove(&i.fT[1][ndex2], &i.fT[1][ndex2 + 1], (i.fUsed2 - ndex2) * sizeof(i.fT[1][0])); |
| 180 | memmove(&matches[ndex2], &matches[ndex2 + 1], (i.fUsed2 - ndex2) * sizeof(matches[0])); |
caryclark@google.com | 235f56a | 2012-09-14 14:19:30 +0000 | [diff] [blame] | 181 | continue; |
| 182 | } |
| 183 | } |
caryclark@google.com | d168874 | 2012-09-18 20:08:37 +0000 | [diff] [blame] | 184 | ++ndex2; |
caryclark@google.com | 235f56a | 2012-09-14 14:19:30 +0000 | [diff] [blame] | 185 | } |
caryclark@google.com | 6aea33f | 2012-10-09 14:11:58 +0000 | [diff] [blame^] | 186 | i.fFlip = i.fUsed >= 2 && flipCheck[0] > flipCheck[1]; |
caryclark@google.com | 235f56a | 2012-09-14 14:19:30 +0000 | [diff] [blame] | 187 | assert(i.insertBalanced()); |
| 188 | return i.intersected(); |
| 189 | } |