Intersecting two ellipses

gistfile1.ts

interface IEllipse2D {
    o: IVec2D;
    rx: number;
    ry: number;
    rotation: number
}


/**
 * Intersect two arbitrarily rotated ellipses. In the general case this boils down to solving a
 * quartic equation. This can have complex results, that are ignored.
 *
 * Special cases like intersecting two circles or intersecting two congruent ellipses are considered.
 * The latter can be reduced to intersecting a line with one of the ellipses.
 *
 * This could numerically be improved by not carrying so many intermediate results, I guess. Also the
 * rotation to avoid problems with ellipses that are not rotated introduces numeric error as well.
 *
 * References:
 * http://elliotnoma.wordpress.com/2013/04/10/a-closed-form-solution-for-the-intersections-of-two-ellipses/
 * http://en.wikipedia.org/wiki/Quartic_function#General_formula_for_roots
 *
 * ... and some Wolfram alpha.
 *
 * @param ellipse1
 * @param ellipse2
 * @returns {null}
 */
export function intersectEE(ellipse1:IEllipse2D, ellipse2:IEllipse2D):IVec2D[] {


    /**
     * Create a general quadratic function for the ellipse a x^2 + b x y + c y^2 + d x + e y + c = 0
     * @param origin
     * @param rx
     * @param ry
     * @param rotation
     * @returns {{a: number, b: number, c: number, d: number, e: number, f: number}}
     */
    function quadratic(origin:IVec2D, rx:number, ry:number, rotation:number):any {
        var o = origin;

        var a = o.x;
        var b = rx * rx;
        var c = o.y;
        var d = ry * ry;
        var A = Math.cos(-rotation);
        var B = Math.sin(-rotation);

        return {
            a: A * A / b + B * B / d,                       // x^2
            b: 2 * A * B / d - 2 * A * B / b,               // x * y
            c: A * A / d + B * B / b,                       // y^2
            d: (2 * A * B * c - 2 * a * A * A) / b + (-2 * a * B * B - 2 * A * B * c) / d,    // x
            e: (2 * a * A * B - 2 * B * B * c) / b  + (-2 * a * A * B - 2 * A * A * c) / d,  // y
            f: (a * a * A * A - 2 * a * A * B * c + B * B * c * c) / b + (a * a * B * B + 2 * a * A * B * c + A * A * c * c) / d - 1  // Const
        }
    }

    /**
     * Calculate the coefficient of the quartic.
     * @param el1
     * @param el2
     * @returns {{z0: number, z1: number, z2: number, z3: number, z4: number}}
     */
    function quartics(el1, el2) {
        var a1 = el1.a;
        var b1 = el1.b;
        var c1 = el1.c;
        var d1 = el1.d;
        var e1 = el1.e;
        var f1 = el1.f;

        var a2 = el2.a;
        var b2 = el2.b;
        var c2 = el2.c;
        var d2 = el2.d;
        var e2 = el2.e;
        var f2 = el2.f;


        return {
            z0: f1 * a1 * d2 * d2 + a1 * a1 * f2 * f2 - d1 * a1 * d2 * f2 + a2 * a2 * f1 * f1 - 2 * a1 * f2 * a2 * f1 - d1 * d2 * a2 * f1 + a2 * d1 * d1 * f2,

            z1: e2 * d1 * d1 * a2 - f2 * d2 * a1 * b1 - 2 * a1 * f2 * a2 * e1 - f1 * a2 * b2 * d1 + 2 * d2 * b2 * a1 * f1 + 2 * e2 * f2 * a1 * a1
                + d2 * d2 * a1 * e1 - e2 * d2 * a1 * d1 - 2 * a1 * e2 * a2 * f1 - f1 * a2 * d2 * b1 + 2 * f1 * e1 * a2 * a2 - f2 * b2 * a1 * d1
                - e1 * a2 * d2 * d1 + 2 * f2 * b1 * a2 * d1,

            z2: e2 * e2 * a1 * a1 + 2 * c2 * f2 * a1 * a1 - e1 * a2 * d2 * b1 + f2 * a2 * b1 * b1 - e1 * a2 * b2 * d1 - f2 * b2 * a1 * b1 - 2 * a1 * e2 * a2 * e1
                + 2 * d2 * b2 * a1 * e1 - c2 * d2 * a1 * d1 - 2 * a1 * c2 * a2 * f1 + b2 * b2 * a1 * f1 + 2 * e2 * b1 * a2 * d1 + e1 * e1 * a2 * a2 - c1 * a2 * d2 * d1
                - e2 * b2 * a1 * d1 + 2 * f1 * c1 * a2 * a2 - f1 * a2 * b2 * b1 + c2 * d1 * d1 * a2 + d2 * d2 * a1 * c1 - e2 * d2 * a1 * b1 - 2 * a1 * f2 * a2 * c1,

            z3: -2 * a1 * a2 * c1 * e2 + e2 * a2 * b1 * b1 + 2 * c2 * b1 * a2 * d1 - c1 * a2 * b2 * d1 + b2 * b2 * a1 * e1 -e2 * b2 * a1 * b1 - 2 * a1 * c2 * a2 * e1
                - e1 * a2 * b2 * b1 - c2 * b2 * a1 * d1 + 2 * e2 * c2 * a1 * a1 + 2 * e1 * c1 * a2 * a2 - c1 * a2 * d2 * b1 + 2 * d2 * b2 * a1 * c1 - c2 * d2 * a1 * b1,

            z4: a1 * a1 * c2 * c2 - 2 * a1 * c2 * a2 * c1 + a2 * a2 * c1 * c1 - b1 * a1 * b2 * c2 - b1 * b2 * a2 * c1 + b1 * b1 * a2 * c2 + c1 * a1 * b2 * b2
        }
    }


    /**
     * This basically calculates the rational roots of the quartic
     * @param quartics
     * @returns {Array}
     */
    function getY(quartics) {
        var a = quartics.z4;
        var b = quartics.z3;
        var c = quartics.z2;
        var d = quartics.z1;
        var e = quartics.z0;


        var delta = 256 * a * a * a * e * e * e - 192 * a * a * b * d * e * e - 128 * a * a * c * c * e * e
                  + 144 * a * a * c * d * d * e - 27 * a * a * d * d * d * d
                  + 144 * a * b * b * c * e * e - 6 * a * b * b * d * d * e - 80 * a * b * c * c * d * e
                  + 18 * a * b * c * d * d * d + 16 * a * c * c * c * c * e
                  - 4 * a * c * c * c * d * d - 27 * b * b * b * b * e * e + 18 * b * b * b * c * d * e
                     - 4 * b * b * b * d * d * d - 4 * b * b * c * c * c * e + b * b * c * c * d * d;

        var P = 8 * a * c - 3 * b * b;
        var D = 64 * a * a * a * e - 16 * a * a * c * c + 16 * a * b * b * c - 16 * a * a * b * d - 3 * b * b * b * b;

        var d0 = c * c - 3 * b * d + 12 * a * e;
        var d1 = 2 * c * c * c - 9 * b * c * d + 27 * b * b * e + 27 * a * d * d - 72 * a * c * e;

        var p = (8 * a * c - 3 * b * b) / (8 * a * a);
        var q = (b * b * b - 4 * a * b * c + 8 * a * a * d) / (8 * a * a * a);

        var Q, S;
        var phi = Math.acos(d1 / (2 * Math.sqrt(d0 * d0 * d0)));

        if ((isNaN(phi) && is0(d1))) { // if (delta < 0) I guess the new test is ok because we're only interested in real roots
            Q = d1 + Math.sqrt(d1 * d1 - 4 * d0 * d0 * d0);
            Q = Q / 2;
            Q = Math.pow(Q, 1/3);
            S = 0.5 * Math.sqrt(- 2 / 3 * p  + (1 / (3 * a)) * (Q + d0 / Q));
        } else {
            S = 0.5 * Math.sqrt(- 2 / 3 * p  + 2 / (3 * a) * Math.sqrt(d0) * Math.cos(phi / 3));
        }

        var y = [];

        if (!is0(S)) {
            var R = - 4 * S * S - 2 * p + q / S;
            if (is0(R)) {
                R = 0;
            }

            if (R > 0) {
                R = 0.5 * Math.sqrt(R);
                y.push(- b / (4 * a) - S + R);
                y.push(- b / (4 * a) - S - R);
            } else if (R === 0) {
                y.push(- b / (4 * a) - S);
            }

            R = - 4 * S * S - 2 * p - q / S;
            if (is0(R)) {
                R = 0;
            }
            if (R > 0) {
                R = 0.5 * Math.sqrt(R);
                y.push(- b / (4 * a) + S + R);
                y.push(- b / (4 * a) + S - R);
            } else if (R === 0) {
                y.push(- b / (4 * a) + S);
            }
        }


        return y;
    }

    /**
     * Calculate the x coordinates for the given y coordinates
     * @param y
     * @param e1        Quartics of ellipse1
     * @param e2        Quartics of ellipse2
     * @returns {IVec2D[]}
     */
    function calculatePoints(y, e1, e2):IVec2D[] {
        var a = e1.a;
        var b = e1.b;
        var c = e1.c;
        var d = e1.d;
        var e = e1.e;
        var f = e1.f;

        var a1 = e2.a;
        var b1 = e2.b;
        var c1 = e2.c;
        var d1 = e2.d;
        var e1 = e2.e;
        var fq = e2.f;


        var r:IVec2D[] = [];
        for (var i = 0; i < y.length; i++) {

            var x = -(a * fq + a * c1 * y[i] * y[i] - a1 * c * y[i] * y[i] + a * e1 * y[i] - a1 * e * y[i] - a1 * f)
                     / (a * b1 * y[i] + a * d1 - a1 * b * y[i] - a1 *d);

            r.push(createVec2D(x, y[i]))
        }

        return r;
    }

    /**
     * Calculates the line that runs through the intersection points of two
     * congruent ellipses with the same rotation.
     * @param rotation
     * @param rx
     * @param ry
     * @param o1
     * @param o2
     * @returns {ILine2D}
     */
    function getLine(rotation:number, rx:number, ry:number, o1:IVec2D, o2:IVec2D):ILine2D {
        var A = Math.cos(rotation);
        var B = Math.sin(rotation);
        var b = rx * rx;
        var d = ry * ry;
        var a = o1.x;
        var c = o1.y;
        var o = o2.x;
        var p = o2.y;

        var AA = A * A / b + B * B / d;
        var BB = -2 * A * B / b + 2 * A * B / d;
        var CC = B * B / b + A * A / d;

        var U = -2 * AA * a + BB * c;
        var V = AA * a * a + BB * a * c + CC * c * c;
        var W = BB * a + 2 * CC * c;

        var X = -2 * AA * o + BB * p;
        var Y = BB * o + 2 * CC * p;
        var Z = AA * o * o + BB * o * p + CC * p * p;

        return createLine(createVec2D(U - X, Y - W), Z - V);
    }



    var e1, e2;
    if (eq(ellipse1.rx, ellipse2.ry) && eq(ellipse2.rx, ellipse2.ry)) {
        //Special case: Two circles
        return intersectCC({
            o: ellipse1.o,
            r: ellipse1.rx
        }, {
            o: ellipse2.o,
            r: ellipse2.rx
        })
    } else
    if ((eq(ellipse1.rotation, ellipse2.rotation) || (eq(Math.abs(ellipse1.rotation - ellipse2.rotation), Math.PI)))
        && eq(ellipse1.rx, ellipse2.rx) && eq(ellipse1.ry, ellipse2.ry)) {
        // Special cases congruent ellipses incl. rotation
        // There are at max two intersection points: We can construct a line that runs through these points
        var l = getLine(ellipse1.rotation, ellipse1.rx, ellipse1.ry, ellipse1.o, ellipse2.o);
        return intersectLE(l, ellipse1);
    } else if ((eq(Math.abs(ellipse1.rotation - ellipse2.rotation), Math.PI / 2) || (eq(Math.abs(ellipse1.rotation - ellipse2.rotation), Math.PI * 3 / 2)))&&
        // Special cases congruent ellipses incl. rotation but one is 90 rotated and the radius sizes are swapped
        // There are at max two intersection points: We can construct a line that runs through these points
        eq(ellipse1.rx, ellipse2.ry) && eq(ellipse1.ry, ellipse2.rx)) {
        var l = getLine(ellipse1.rotation, ellipse1.rx, ellipse1.ry, ellipse1.o, ellipse2.o);
        return intersectLE(l, ellipse1);
    } else {
        // General case
        /*
         Test for one situation:
         One of the ellipses axis is parallel to x- or y-axis.
         To avoid special cases testing in getY we simply rotate everything around ellipse1 origin by something and later
         rotate the results back.
         */
        var mPI1 = ellipse1.rotation % Math.PI;
        var mPI2 = ellipse2.rotation % Math.PI;
        var corr = 0;

        if (is0(mPI1) || is0(mPI2)) {
            corr = 0.05;
        }

        if (corr) {
            e1 = quadratic(ellipse1.o, ellipse1.rx, ellipse1.ry, ellipse1.rotation + corr);
            e2 = quadratic(rotateVector(ellipse2.o, corr, ellipse1.o), ellipse2.rx, ellipse2.ry, ellipse2.rotation + corr);
        } else {
            e1 = quadratic(ellipse1.o, ellipse1.rx, ellipse1.ry, ellipse1.rotation);
            e2 = quadratic(ellipse2.o, ellipse2.rx, ellipse2.ry, ellipse2.rotation);
        }

        var q = quartics(e1, e2);
        var y = getY(q);

        var v = calculatePoints(y, e1, e2);

        if (corr) {
            for (var i = 0; i < v.length; i++) {
                v[i] = rotateVector(v[i], -corr, ellipse1.o);
            }
        }

    }
    return v;
}