Java Curve Intersect intersectQuadAndCubic(double qx1, double qy1, double qx2, double qy2, double qx3, double qy3, double cx1, double cy1, double cx2, double cy2, double cx3, double cy3, double cx4, double cy4, double[] params)

Description

Checks whether the quad (x1, y1) - (x2, y2) - (x3, y3) and the cubic (cx1, cy1) - (cx2, cy2) - (cx3, cy3) - (cx4, cy4) curves intersect.

License

Apache License

Return

the number of intersection points that lie in the interval.

Declaration

public static int intersectQuadAndCubic(double qx1, double qy1, double qx2, double qy2, double qx3, double qy3,
        double cx1, double cy1, double cx2, double cy2, double cx3, double cy3, double cx4, double cy4,
        double[] params) 

Method Source Code

//package com.java2s;
//License from project: Apache License 

public class Main {
    public static final double EPSILON = Math.pow(10, -14);

    /**// w  w  w  . j  a v  a2s. c om
     * Checks whether the quad (x1, y1) - (x2, y2) - (x3, y3) and the cubic (cx1, cy1) - (cx2, cy2)
     * - (cx3, cy3) - (cx4, cy4) curves intersect. The points of the intersection are saved to
     * {@code params}. Thus {@code params} must be of length at least 6.
     *
     * @return the number of intersection points that lie in the interval.
     */
    public static int intersectQuadAndCubic(double qx1, double qy1, double qx2, double qy2, double qx3, double qy3,
            double cx1, double cy1, double cx2, double cy2, double cx3, double cy3, double cx4, double cy4,
            double[] params) {
        int quantity = 0;
        double[] initParams = new double[3];
        double[] xCoefs1 = new double[3];
        double[] yCoefs1 = new double[3];
        double[] xCoefs2 = new double[4];
        double[] yCoefs2 = new double[4];
        xCoefs1[0] = qx1 - 2 * qx2 + qx3;
        xCoefs1[1] = 2 * qx2 - 2 * qx1;
        xCoefs1[2] = qx1;

        yCoefs1[0] = qy1 - 2 * qy2 + qy3;
        yCoefs1[1] = 2 * qy2 - 2 * qy1;
        yCoefs1[2] = qy1;

        xCoefs2[0] = -cx1 + 3 * cx2 - 3 * cx3 + cx4;
        xCoefs2[1] = 3 * cx1 - 6 * cx2 + 3 * cx3;
        xCoefs2[2] = -3 * cx1 + 3 * cx2;
        xCoefs2[3] = cx1;

        yCoefs2[0] = -cy1 + 3 * cy2 - 3 * cy3 + cy4;
        yCoefs2[1] = 3 * cy1 - 6 * cy2 + 3 * cy3;
        yCoefs2[2] = -3 * cy1 + 3 * cy2;
        yCoefs2[3] = cy1;

        // initialize params[0] and params[1]
        params[0] = params[1] = 0.25f;
        quadAndCubicNewton(xCoefs1, yCoefs1, xCoefs2, yCoefs2, initParams);
        if (initParams[0] <= 1 && initParams[0] >= 0 && initParams[1] >= 0 && initParams[1] <= 1) {
            params[2 * quantity] = initParams[0];
            params[2 * quantity + 1] = initParams[1];
            ++quantity;
        }
        // initialize params
        params[0] = params[1] = 0.5f;
        quadAndCubicNewton(xCoefs1, yCoefs1, xCoefs2, yCoefs2, params);
        if (initParams[0] <= 1 && initParams[0] >= 0 && initParams[1] >= 0 && initParams[1] <= 1) {
            params[2 * quantity] = initParams[0];
            params[2 * quantity + 1] = initParams[1];
            ++quantity;
        }

        params[0] = params[1] = 0.75f;
        quadAndCubicNewton(xCoefs1, yCoefs1, xCoefs2, yCoefs2, params);
        if (initParams[0] <= 1 && initParams[0] >= 0 && initParams[1] >= 0 && initParams[1] <= 1) {
            params[2 * quantity] = initParams[0];
            params[2 * quantity + 1] = initParams[1];
            ++quantity;
        }
        return quantity;
    }

    private static void quadAndCubicNewton(double xCoefs1[], double yCoefs1[], double xCoefs2[], double yCoefs2[],
            double[] params) {
        double t = 0f, s = 0f;
        double t1 = params[0];
        double s1 = params[1];
        double d, dt, ds;

        while (Math.sqrt((t - t1) * (t - t1) + (s - s1) * (s - s1)) > EPSILON) {
            d = -(2 * t * xCoefs1[0] + xCoefs1[1]) * (3 * s * s * yCoefs2[0] + 2 * s * yCoefs2[1] + yCoefs2[2])
                    + (2 * t * yCoefs1[0] + yCoefs1[1])
                            * (3 * s * s * xCoefs2[0] + 2 * s * xCoefs2[1] + xCoefs2[2]);

            dt = (t * t * xCoefs1[0] + t * xCoefs1[1] + xCoefs1[2] + -s * s * s * xCoefs2[0] - s * s * xCoefs2[1]
                    - s * xCoefs2[2] - xCoefs2[3]) * (-3 * s * s * yCoefs2[0] - 2 * s * yCoefs2[1] - yCoefs2[2])
                    + (t * t * yCoefs1[0] + t * yCoefs1[1] + yCoefs1[2] - s * s * s * yCoefs2[0]
                            - s * s * yCoefs2[1] - s * yCoefs2[2] - yCoefs2[3])
                            * (3 * s * s * xCoefs2[0] + 2 * s * xCoefs2[1] + xCoefs2[2]);

            ds = (2 * t * xCoefs1[0] + xCoefs1[1])
                    * (t * t * yCoefs1[0] + t * yCoefs1[1] + yCoefs1[2] - s * s * s * yCoefs2[0]
                            - s * s * yCoefs2[1] - s * yCoefs2[2] - yCoefs2[3])
                    - (2 * t * yCoefs1[0] + yCoefs1[1]) * (t * t * xCoefs1[0] + t * xCoefs1[1] + xCoefs1[2]
                            - s * s * s * xCoefs2[0] - s * s * xCoefs2[1] - s * xCoefs2[2] - xCoefs2[3]);

            t1 = t - dt / d;
            s1 = s - ds / d;
        }
        params[0] = t1;
        params[1] = s1;
    }
}

Related

1. intersectQuads(double x1, double y1, double x2, double y2, double x3, double y3, double qx1, double qy1, double qx2, double qy2, double qx3, double qy3, double[] params)