查找相交的线和三次样条之间的点

Question

我需要一种方法以编程找到的线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/

Answer 1

代码

if (x >= 0) 
{ // One real root and two imaginaries. 
    x = Math.Sqrt(x); 
    x = CubicRoot(q + x) + CubicRoot(q - x) + p; 
} 
else 
{ // Three real roots. 
    x = Math.Sqrt(Math.Abs(x)); 


    x_1 = Math.Sign(q)*2*(q*q + x*x)^(1/6)*Math.Cos(Math.Atan(x/q)/3) + p; 
    x_2 = Math.Sign(q)*2*(q*q + x*x)^(1/6)*Math.Cos(Math.Atan(x/q)/3 + Math.PI*2/3) + p; 
    x_3 = Math.Sign(q)*2*(q*q + x*x)^(1/6)*Math.Cos(Math.Atan(x/q)/3 + Math.PI*4/3) + p; 
} 

你可以计算()^(1/6)Math.Pow( , 1/6)或您的Math.Sqrt(CubicRoot())或Math.Sqrt(Cbrt())。请参阅following thread on Microsoft forum。

请注意q = 0。 (Atan(x/0) = pi/2 radians. Cos(Atan(x/0)/3) = Sqrt(3)/2)

Answer 2

对于二次方程，存在至少1个单根实根。 使用此方法寻找根http://en.wikipedia.org/wiki/Cubic_function#General_formula_of_roots