Нахождение точки пересечения прямой и кубического сплайна

Мне нужен способ программно найти точку пересечения между прямой f (x), которая берет начало от начала координат, и кубическим сплайном, определенным 4 точками. Линия гарантированно пересекает центральный сегмент сплайна, между X1 и X2.

Я попробовал несколько подходов, но, похоже, не могу получить ожидаемого результата. Я подозреваю, что моя проблема кроется где-то в моем обращении с комплексными числами.

Может ли кто-нибудь найти проблему с моим кодом, или предложить другой подход?

    private Vector2 CubicInterpolatedIntersection(float y0, float y1, 

                            float y2, float y3, float lineSlope, lineYOffset)

    {

        // f(x) = lineSlope * x + lineYOffset

        // g(x) = (a0 * x^3) + (a1 * x^2) + (a2 * x) + a3



        // Calculate Catmull-Rom coefficients for g(x) equation as found 

        // in reference (1). These 

        double a0, a1, a2, a3;

        a0 = -0.5 * y0 + 1.5 * y1 - 1.5 * y2 + 0.5 * y3;

        a1 = y0 - 2.5 * y1 + 2 * y2 - 0.5 * y3;

        a2 = -0.5 * y0 + 0.5 * y2;

        a3 = y1;



        // To find POI, let 'g(x) - f(x) = 0'. Simplified equation is:

        // (a0 * x^3) + (a1 * x^2) + ((a2 - lineSlope) * x) 

        //                         + (a3 - lineYOffset) = 0



        // Put in standard form 'ax^3 + bx^2 + cx + d = 0'

        double a, b, c, d;

        a = a0;

        b = a1;

        c = a2 - lineSlope;

        d = a3 - lineYOffset;





        // Solve for roots using cubic equation as found in reference (2).

        //  x = {q + [q^2 + (r-p^2)^3]^(1/2)}^(1/3) 

        //         + {q - [q^2 + (r-p^2)^3]^(1/2)}^(1/3) + p

        //  Where...

        double p, q, r;

        p = -b / (3 * a);

        q = p * p * p + (b * c - 3 * a * d) / (6 * a * a);

        r = c / (3 * a);



        //Solve the equation starting from the center.

        double x, x2;

        x = r - p * p;

        x = x * x * x + q * q;

        // Here's where I'm not sure. The cubic equation contains a square 

        // root. So if x is negative at this point, then we need to proceed

        // with complex numbers.

        if (x >= 0)

        {

            x = Math.Sqrt(x);

            x = CubicRoot(q + x) + CubicRoot(q - x) + p;

        }

        else

        {

            x = Math.Sqrt(Math.Abs(x)); 

            // x now represents the imaginary component of 

            //                       a complex number 'a + b*i'

            // We need to take the cubic root of 'q + x' and 'q - x'

            // Since x is imaginary, we have two complex numbers in 

            // standard form. Therefore, we take the cubic root of 

            // the magnitude of these complex numbers

            x = CubicRoot(Math.Sqrt(q * q + x * x)) + 

                 CubicRoot(Math.Sqrt(q * q + -x * -x)) + p;

        }



        // At this point, x should hold the x-value of 

        //      the point of intersection.

        // Now use it to solve g(x).



        x2 = x * x;

        return new Vector2((float)Math.Abs(x), 

           (float)Math.Abs(a0 * x * x2 + a1 * x2 + a2 * x + a3));

    }

Ссылки:

1. 
   http://paulbourke.net/miscellaneous/interpolation/


2. 
   http://www.math.vanderbilt.edu/~schectex / courses / cubic /