QuadraticImplicit.cpp revision 055c7c299cb47eebd360b809ad58a0006e2e55f7 
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    int rootCount2 = findRoots(i1, q2, roots2);
95    assert(rootCount == rootCount2);
96    addValidRoots(roots1, rootCount, 0, i);
97    addValidRoots(roots2, rootCount, 1, i);
98    _Point pts[4];
99    bool matches[4];
100    int index;
101    for (index = 0; index < i.fUsed2; ++index) {
102        xy_at_t(q2, i.fT[1][index], pts[index].x, pts[index].y);
103        matches[index] = false;
104    }
105    for (index = 0; index < i.fUsed; ) {
106        _Point xy;
107        xy_at_t(q1, i.fT[0][index], xy.x, xy.y);
108        for (int inner = 0; inner < i.fUsed2; ++inner) {
109             if (approximately_equal(pts[inner].x, xy.x) && approximately_equal(pts[inner].y, xy.y)) {
110                matches[index] = true;
111                goto next;
112             }
113        }
114        if (--i.fUsed > index) {
115            memmove(&i.fT[0][index], &i.fT[0][index + 1], (i.fUsed - index) * sizeof(i.fT[0][0]));
116            continue;
117        }
118    next:
119        ++index;
120    }
121    for (index = 0; index < i.fUsed2; ) {
122        if (!matches[index]) {
123             if (--i.fUsed2 > index) {
124                memmove(&i.fT[1][index], &i.fT[1][index + 1], (i.fUsed2 - index) * sizeof(i.fT[1][0]));
125                continue;
126             }
127        }
128        ++index;
129    }
130    assert(i.insertBalanced());
131    return i.intersected();
132}
133