QuadraticImplicit.cpp revision d1688744d537d928699b6069f99c4470a0f6e772
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 17 * then 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 21static 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 51static 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 67bool intersect2(const Quadratic& q1, const Quadratic& q2, Intersections& i) { 68 QuadImplicitForm i1(q1); 69 QuadImplicitForm i2(q2); 70 if (i1.implicit_match(i2)) { 71 // FIXME: compute T values 72 // compute the intersections of the ends to find the coincident span 73 bool useVertical = fabs(q1[0].x - q1[2].x) < fabs(q1[0].y - q1[2].y); 74 double t; 75 if ((t = axialIntersect(q1, q2[0], useVertical)) >= 0) { 76 i.addCoincident(t, 0); 77 } 78 if ((t = axialIntersect(q1, q2[2], useVertical)) >= 0) { 79 i.addCoincident(t, 1); 80 } 81 useVertical = fabs(q2[0].x - q2[2].x) < fabs(q2[0].y - q2[2].y); 82 if ((t = axialIntersect(q2, q1[0], useVertical)) >= 0) { 83 i.addCoincident(0, t); 84 } 85 if ((t = axialIntersect(q2, q1[2], useVertical)) >= 0) { 86 i.addCoincident(1, t); 87 } 88 assert(i.fCoincidentUsed <= 2); 89 return i.fCoincidentUsed > 0; 90 } 91 double roots1[4], roots2[4]; 92 int rootCount = findRoots(i2, q1, roots1); 93 // OPTIMIZATION: could short circuit here if all roots are < 0 or > 1 94#ifndef NDEBUG 95 int rootCount2 = 96#endif 97 findRoots(i1, q2, roots2); 98 assert(rootCount == rootCount2); 99 addValidRoots(roots1, rootCount, 0, i); 100 addValidRoots(roots2, rootCount, 1, i); 101 _Point pts[4]; 102 bool matches[4]; 103 int index, ndex2; 104 for (ndex2 = 0; ndex2 < i.fUsed2; ++ndex2) { 105 xy_at_t(q2, i.fT[1][ndex2], pts[ndex2].x, pts[ndex2].y); 106 matches[ndex2] = false; 107 } 108 for (index = 0; index < i.fUsed; ) { 109 _Point xy; 110 xy_at_t(q1, i.fT[0][index], xy.x, xy.y); 111 for (ndex2 = 0; ndex2 < i.fUsed2; ++ndex2) { 112 if (approximately_equal(pts[ndex2].x, xy.x) && approximately_equal(pts[ndex2].y, xy.y)) { 113 matches[ndex2] = true; 114 goto next; 115 } 116 } 117 if (--i.fUsed > index) { 118 memmove(&i.fT[0][index], &i.fT[0][index + 1], (i.fUsed - index) * sizeof(i.fT[0][0])); 119 continue; 120 } 121 next: 122 ++index; 123 } 124 for (ndex2 = 0; ndex2 < i.fUsed2; ) { 125 if (!matches[ndex2]) { 126 if (--i.fUsed2 > ndex2) { 127 memmove(&i.fT[1][ndex2], &i.fT[1][ndex2 + 1], (i.fUsed2 - ndex2) * sizeof(i.fT[1][0])); 128 memmove(&matches[ndex2], &matches[ndex2 + 1], (i.fUsed2 - ndex2) * sizeof(matches[0])); 129 continue; 130 } 131 } 132 ++ndex2; 133 } 134 assert(i.insertBalanced()); 135 return i.intersected(); 136} 137