blob: 268d7d3f624087d5fefe2734a35c6370ad222843 [file] [log] [blame]
// Another approach is to start with the implicit form of one curve and solve
// (seek implicit coefficients in QuadraticParameter.cpp
// by substituting in the parametric form of the other.
// The downside of this approach is that early rejects are difficult to come by.
// http://planetmath.org/encyclopedia/GaloisTheoreticDerivationOfTheQuarticFormula.html#step
#include "CurveIntersection.h"
#include "Intersections.h"
#include "QuadraticParameterization.h"
#include "QuarticRoot.h"
#include "QuadraticUtilities.h"
/* given the implicit form 0 = Ax^2 + Bxy + Cy^2 + Dx + Ey + F
* and given x = at^2 + bt + c (the parameterized form)
* y = dt^2 + et + f
* then
* 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
*/
static int findRoots(const QuadImplicitForm& i, const Quadratic& q2, double roots[4]) {
double a, b, c;
set_abc(&q2[0].x, a, b, c);
double d, e, f;
set_abc(&q2[0].y, d, e, f);
const double t4 = i.x2() * a * a
+ i.xy() * a * d
+ i.y2() * d * d;
const double t3 = 2 * i.x2() * a * b
+ i.xy() * (a * e + b * d)
+ 2 * i.y2() * d * e;
const double t2 = i.x2() * (b * b + 2 * a * c)
+ i.xy() * (c * d + b * e + a * f)
+ i.y2() * (e * e + 2 * d * f)
+ i.x() * a
+ i.y() * d;
const double t1 = 2 * i.x2() * b * c
+ i.xy() * (c * e + b * f)
+ 2 * i.y2() * e * f
+ i.x() * b
+ i.y() * e;
const double t0 = i.x2() * c * c
+ i.xy() * c * f
+ i.y2() * f * f
+ i.x() * c
+ i.y() * f
+ i.c();
return quarticRoots(t4, t3, t2, t1, t0, roots);
}
static void addValidRoots(const double roots[4], const int count, const int side, Intersections& i) {
int index;
for (index = 0; index < count; ++index) {
if (!approximately_zero_or_more(roots[index]) || !approximately_one_or_less(roots[index])) {
continue;
}
double t = 1 - roots[index];
if (approximately_less_than_zero(t)) {
t = 0;
} else if (approximately_greater_than_one(t)) {
t = 1;
}
i.insertOne(t, side);
}
}
static bool onlyEndPtsInCommon(const Quadratic& q1, const Quadratic& q2, Intersections& i) {
// the idea here is to see at minimum do a quick reject by rotating all points
// to either side of the line formed by connecting the endpoints
// if the opposite curves points are on the line or on the other side, the
// curves at most intersect at the endpoints
for (int oddMan = 0; oddMan < 3; ++oddMan) {
const _Point* endPt[2];
for (int opp = 1; opp < 3; ++opp) {
int end = oddMan ^ opp;
if (end == 3) {
end = opp;
}
endPt[opp - 1] = &q1[end];
}
double origX = endPt[0]->x;
double origY = endPt[0]->y;
double adj = endPt[1]->x - origX;
double opp = endPt[1]->y - origY;
double sign = (q1[oddMan].y - origY) * adj - (q1[oddMan].x - origX) * opp;
assert(!approximately_zero(sign));
for (int n = 0; n < 3; ++n) {
double test = (q2[n].y - origY) * adj - (q2[n].x - origX) * opp;
if (test * sign > 0) {
goto tryNextHalfPlane;
}
}
for (int i1 = 0; i1 < 3; i1 += 2) {
for (int i2 = 0; i2 < 3; i2 += 2) {
if (q1[i1] == q2[i2]) {
i.insertOne(i1 >> 1, 0);
i.insertOne(i2 >> 1, 1);
}
}
}
assert(i.fUsed < 3);
return true;
tryNextHalfPlane:
;
}
return false;
}
bool intersect2(const Quadratic& q1, const Quadratic& q2, Intersections& i) {
// if the quads share an end point, check to see if they overlap
if (onlyEndPtsInCommon(q1, q2, i)) {
assert(i.insertBalanced());
return i.intersected();
}
QuadImplicitForm i1(q1);
QuadImplicitForm i2(q2);
if (i1.implicit_match(i2)) {
// FIXME: compute T values
// compute the intersections of the ends to find the coincident span
bool useVertical = fabs(q1[0].x - q1[2].x) < fabs(q1[0].y - q1[2].y);
double t;
if ((t = axialIntersect(q1, q2[0], useVertical)) >= 0) {
i.addCoincident(t, 0);
}
if ((t = axialIntersect(q1, q2[2], useVertical)) >= 0) {
i.addCoincident(t, 1);
}
useVertical = fabs(q2[0].x - q2[2].x) < fabs(q2[0].y - q2[2].y);
if ((t = axialIntersect(q2, q1[0], useVertical)) >= 0) {
i.addCoincident(0, t);
}
if ((t = axialIntersect(q2, q1[2], useVertical)) >= 0) {
i.addCoincident(1, t);
}
assert(i.fCoincidentUsed <= 2);
return i.fCoincidentUsed > 0;
}
double roots1[4], roots2[4];
int rootCount = findRoots(i2, q1, roots1);
// OPTIMIZATION: could short circuit here if all roots are < 0 or > 1
#ifndef NDEBUG
int rootCount2 =
#endif
findRoots(i1, q2, roots2);
assert(rootCount == rootCount2);
addValidRoots(roots1, rootCount, 0, i);
addValidRoots(roots2, rootCount, 1, i);
_Point pts[4];
bool matches[4];
int flipCheck[4];
int index, ndex2;
int flipIndex = 0;
for (ndex2 = 0; ndex2 < i.fUsed2; ++ndex2) {
xy_at_t(q2, i.fT[1][ndex2], pts[ndex2].x, pts[ndex2].y);
matches[ndex2] = false;
}
for (index = 0; index < i.fUsed; ) {
_Point xy;
xy_at_t(q1, i.fT[0][index], xy.x, xy.y);
for (ndex2 = 0; ndex2 < i.fUsed2; ++ndex2) {
if (approximately_equal(pts[ndex2].x, xy.x) && approximately_equal(pts[ndex2].y, xy.y)) {
assert(flipIndex < 4);
flipCheck[flipIndex++] = ndex2;
matches[ndex2] = true;
goto next;
}
}
if (--i.fUsed > index) {
memmove(&i.fT[0][index], &i.fT[0][index + 1], (i.fUsed - index) * sizeof(i.fT[0][0]));
continue;
}
next:
++index;
}
for (ndex2 = 0; ndex2 < i.fUsed2; ) {
if (!matches[ndex2]) {
if (--i.fUsed2 > ndex2) {
memmove(&i.fT[1][ndex2], &i.fT[1][ndex2 + 1], (i.fUsed2 - ndex2) * sizeof(i.fT[1][0]));
memmove(&matches[ndex2], &matches[ndex2 + 1], (i.fUsed2 - ndex2) * sizeof(matches[0]));
continue;
}
}
++ndex2;
}
i.fFlip = i.fUsed >= 2 && flipCheck[0] > flipCheck[1];
assert(i.insertBalanced());
return i.intersected();
}