Sorry for re-raising an oldish thread, but I wanted to give you guys an update and post the now-working code in case someone may find it useful.
I've been working on this more slowly since my last post due to time constraints, but yesterday I finally got the torus to function properly. Bellow is a screenshot and the current version of my code. The solver code isn't 100% polished etc, but it works.
For those who are interested, here is the details of my torus adventures since the last post. Those who aren't interested, may skip to the screenshot and code.
First, I implemented a hybrid approach: cubic (and quadratic) solver from the python source, and quartic solver from the thread I linked to earlier (with some smalls adjustments to fit the python source's system of returning all roots). After removing the check for too large B, C, or D from the start of the quartic solver, this gave me a correct torus with most of the numerical errors eliminated. Apparently my quardic solver was correct all along, except for that erroneous check. But I was still getting numerical errors when (apparently) the ray was at some constant angle with respect to the torus surface.
Then I went back to my own version of the cubic solver, based on the link I gave earlier. I found an error in my implementation, where I've put a plus instead of a minus. Fixing this gave me correct results with fewer numerical errors than the python solver. I'm still using the python code's approach of returning all the roots in an array, but that's about all that was left from the porting effort I wrote about in my last post (that and the quadratic equation solver, but that's easy to reproduce). The non-python method also seemed to be ~30% faster, though it may be view dependant..
Next I tried fine-tuning the intermediate root I get from the cubic solver using the newton-raphson method. At first glance, this seemed to eliminate all the remaining numerical errors. However, further examinations showed that numerical error was in fact still present when viewing the torus at a large distance. Solving entirely using newton-raphson delayed the appearance of the errors until larger distances, but did not eliminate them, and was much slower than the (semi) analytic method.
The solution was to use higher-precision floating point for the solvers (double instead of single precision). The faster analytic method gave sufficiently accurate results in double precision mode.
Then I found a bug where overlapping sections of the torus retruned "no intersection" (false negative). This seemed to be due to some innacuracy in the quartic solver or the equation parameter calculations. I found that as a result of these innacuracies, the cubic solver (on which the quartic solver relies) did not find any roots (since it used an intermediate quadratic equation that came out unsolvable). Since I know that a cubic equation must have at least one real root, I added code for the cubic solver to fall back to the newton-raphson method if the analytic method fails. This solved the remaining problem, but noticeably slowed down the interesction test.
Finally, adding a preliminary bounding sphere inetersection test returned much of the speed (perhaps even made it faster than the original), since much of the time was being wasted on calculating misses.
Screenshot:
Polynom solvers and related utilities:
unit Polynoms;{******************************************************************************}{ Polynoms - polynom solvers and utilities }{ }{ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -}{ }{ Author: Michael Kuntzman [technobot@users.sourceforge.net] }{ Original file name: 'Polynoms.pas' }{ Created on: Jul 1st, 2007 }{ }{ Parts ported and adapted from }{ http://pyx.svn.sourceforge.net/viewvc/pyx/trunk/pyx/pyx/ }{ mathutils.py?view=log , written by: }{ - Michael Schindler [m-schindler@users.sourceforge.net] }{ - Andr Wobst [wobsta@users.sourceforge.net] }{ }{ Additional parts ported and adapted from }{ http://www.gamedev.net/community/forums/topic.asp?topic_id=451048 }{ }{ Further additional parts based on }{ http://www.sosmath.com/algebra/factor/fac11/fac11.html and }{ http://www.math.ubc.ca/~duncan/newtonmethod.pdf }{ }{ The contents of this file is subject to the GNU General Public License as }{ published by the Free Software Foundation; either version 2 of the License, }{ or (at your option) any later version. A copy of the license is available }{ at http://www.fsf.org/licensing/licenses/gpl.txt . }{ }{******************************************************************************}interfaceconst // The maximum allowed polynom degree: MaxPolDegree = 10;type // Single precision floating point tends to produce significant numerical // error in certain situations. Double precision seems to give adequate // results: TFloat = Double; TPolynomDegree = 0..MaxPolDegree; // An array of up to 255 floats: TFloatArray = array [TPolynomDegree] of TFloat; PFloatArray = ^TFloatArray;var // The maximum number of iterations allowed for the FineTuneRoot function: MaxFTRIterations: Cardinal = 20; //---------------------------------------------------------------------------- // Differentiates the given polynom. A are the original polynom coefficients, // DiffA will hold the coefficients of the differential, and Degree is the // highest power of X in the original polynom. IMPORTANT: This function does // NOT check the float array pointers. It is the user's responsibility to // provide valid arrays of suffcient size!: procedure DifferentiatePol(const A, DiffA: PFloatArray; const Degree: TPolynomDegree); //---------------------------------------------------------------------------- // Divides the given polynom by (x-x1), where x1 is the given root. A are the // polynom coefficients, NewA will hold the coefficients of the resulting // polynom, and Degree is the highest power of X in the original polynom. The // resulting polynom will be of degree Degree-1. IMPORTANT: This function // does NOT check the float array pointers, or the supplied root. It is the // user's responsibility to provide valid arrays of suffcient size, and a // correct root!: procedure DividePol(const A, NewA: PFloatArray; const Root: TFloat; Degree: TPolynomDegree); //---------------------------------------------------------------------------- // Evaluates the given polynom at the given X. A are the polynom coefficients, // and Degree is the highest power of X in the polynom. IMPORTANT: This // function does NOT check the float array pointer. It is the user's // responsibility to provide a valid array of suffcient size!: function EvaluatePol(const A: PFloatArray; const X: TFloat; const Degree: TPolynomDegree): TFloat; //---------------------------------------------------------------------------- // Finds a root of the given polynomial in the given interval. Returns True on // success, or False on failure. For best results, there should be only one // root in the given interval, and the interval should be small. Otherwise // the function may not find any root, or may take longer to find it. For // large intervals, the MaxFTRIterations global variable may need to be // increased. A are the polynom coefficients of the polynom, X1 and X2 define // the interval (X2>X1), MaxError (>0) is the maximum allowed error in X, // Root will hold the root if it is found, and Degree is the highest power of // X in the polynom. IMPORTANT: This function does NOT check the values of // X1, X2, and MaxError, or the float array pointer. It is the user's // responsibility to provide correct values and a valid array of correct // size!: function FindRoot(const A: PFloatArray; const X1, X2, MaxError: TFloat; out Root: TFloat; const Degree: TPolynomDegree): Boolean; //---------------------------------------------------------------------------- // Fine-tunes the proposed root of the given polynomial using the secant // variation of the newton-raphson method. Returns True on success, or False // on failure. If the function fails, Root remains unchnaged. For best // results, the suggested root should be as close as possible to the actual // root. Otherwise the function may not find the root, or may take longer to // find it. For large distances, the MaxFTRIterations global variable may // need to be increased. A are the polynom coefficients, MaxError (>0) is the // maximum allowed error in Root, and Degree is the highest power of X in the // polynom. The code is based on the explanations in // http://www.math.ubc.ca/~duncan/newtonmethod.pdf. IMPORTANT: This function // does NOT check the values of X1, X2, and MaxError, or the float array // pointer. It is the user's responsibility to provide correct values and a // valid array of correct size!: function FineTuneRoot(const A: PFloatArray; const MaxError: TFloat; var Root: TFloat; const Degree: TPolynomDegree): Boolean; //---------------------------------------------------------------------------- // Puts the real roots of x^3 + a2*x^2 + a1*x + a0 = 0 in Roots, and returns // the number of roots found. The roots may not be unique. The code is based // mostly on the explanations from // http://www.sosmath.com/algebra/factor/fac11/fac11.html, with some // principles adapted from http://pyx.svn.sourceforge.net/viewvc/pyx/trunk/ // pyx/pyx/mathutils.py?view=log. IMPORTANT: This function does NOT check the // float array pointer. It is the user's responsibility to provide a valid // array of suffcient size!: function SolveCubic(const a2, a1, a0: TFloat; const Roots: PFloatArray): Integer; //---------------------------------------------------------------------------- // Solves the given polynom within the specified interval, puts its real roots // within the interval in Roots, and returns the number of roots found. The // roots are each unique, and are returned in increasing order (from low to // high). For best results, the interval should be as small as possible. // Otherwise the function may take longer to find the roots, or may not find // them at all if the MaxFTRIterations global variable is set too low. For // polynoms of the 4th or lower degree, the functions SolveQartic, // SolveCubic, and SolveQuadratic are sometimes preferable. They are usually // faster and are more general (no interval required), but may be less // accurate. A combination of those functions with FineTuneRoot may prove to // be a good compromise in some cases. A are the polynom coefficients, X1 and // X2 define the interval (X2>X1), MaxError (>0) is the maximum allowed error // in the roots, and Degree is the highest power of X in the polynom. // IMPORTANT: This function does NOT check the values of X1, X2, and // MaxError, or the float array pointers. It is the user's responsibility to // provide correct values and valid arrays of correct size!: function SolvePolynom(const A, Roots: PFloatArray; const X1, X2, MaxError: TFloat; const Degree: TPolynomDegree): Integer; //---------------------------------------------------------------------------- // Puts the real roots of x^2 + a1*x + a0 = 0 in Roots, and returns the number // of roots found. The roots may not be unique. Code adapted from // http://pyx.svn.sourceforge.net/viewvc/pyx/trunk/pyx/pyx/ // mathutils.py?view=log. IMPORTANT: This function does NOT check the float // array pointer. It is the user's responsibility to provide a valid array of // suffcient size!: function SolveQuadratic(const a1, a0: TFloat; const Roots: PFloatArray): Integer; //---------------------------------------------------------------------------- // Puts the real roots of x^4 + a3*x^3 + a2*x^2 + a1*x + a0 = 0 in Roots, and // returns the number of roots found. The roots may not be unique. Code based // mostly on http://www.gamedev.net/community/forums/ // topic.asp?topic_id=451048, with some principles adapted from // http://pyx.svn.sourceforge.net/viewvc/pyx/trunk/pyx/pyx/ // mathutils.py?view=log. IMPORTANT: This function does NOT check the float // array pointer. It is the user's responsibility to provide a valid array of // suffcient size!: function SolveQuartic(const a3, a2, a1, a0: TFloat; const Roots: PFloatArray): Integer;implementation//------------------------------------------------------------------------------// Internal function - computes Base^Exponent:function Power(const Base, Exponent: TFloat): TFloat;begin if (Exponent = 0.0) then Result := 1.0 // n^0 = 1 else if ((Base = 0.0) and (Exponent > 0.0)) then Result := 0.0 // 0^n = 0, n > 0 else Result := Exp(Exponent * Ln(Base));end;//------------------------------------------------------------------------------// Differentiates the given polynom. A are the original polynom coefficients,// DiffA will hold the coefficients of the differential, and Degree is the// highest power of X in the original polynom. IMPORTANT: This function does// NOT check the float array pointers. It is the user's responsibility to// provide valid arrays of suffcient size!:procedure DifferentiatePol(const A, DiffA: PFloatArray; const Degree: TPolynomDegree);var D: TPolynomDegree; Mul: TFloat;begin case Degree of 0: Exit; 1: DiffA[0] := A[1]; 2: begin DiffA[0] := A[1]; DiffA[1] := A[2] * 2.0; end; 3: begin DiffA[0] := A[1]; DiffA[1] := A[2] * 2.0; DiffA[2] := A[3] * 3.0; end; 4: begin DiffA[0] := A[1]; DiffA[1] := A[2] * 2.0; DiffA[2] := A[3] * 3.0; DiffA[3] := A[4] * 4.0; end; else DiffA[0] := A[1]; DiffA[1] := A[2] * 2.0; DiffA[2] := A[3] * 3.0; DiffA[3] := A[4] * 4.0; Mul := 5.0; for D := 5 to Degree do begin DiffA[D-1] := A[D] * Mul; Mul := Mul + 1.0; end; end;end; { DifferentiatePol }//------------------------------------------------------------------------------// Divides the given polynom by (x-x1), where x1 is the given root. A are the// polynom coefficients, NewA will hold the coefficients of the resulting// polynom, and Degree is the highest power of X in the original polynom. The// resulting polynom will be of degree Degree-1. IMPORTANT: This function// does NOT check the float array pointers, or the supplied root. It is the// user's responsibility to provide valid arrays of suffcient size, and a// correct root!:procedure DividePol(const A, NewA: PFloatArray; const Root: TFloat; Degree: TPolynomDegree);var RemainderAn: TFloat;begin // We actually doing long division here. First, find the "remainder" for the // intitial degree. This is actually the the entire polynom as is: RemainderAn := A[Degree]; // Now calculate the lower degrees: while (Degree > 0) do begin // Move to the next degree: Dec(Degree); // The result for this degree is the remainder from the last one (since // we're dividing by x-x1): NewA[Degree] := RemainderAn; // Calculate the next remainder: RemainderAn := A[Degree] - RemainderAn * Root; end;end; { DividePol }//------------------------------------------------------------------------------// Evaluates the given polynom at the given X. A are the polynom coefficients,// and Degree is the highest power of X in the polynom. IMPORTANT: This// function does NOT check the float array pointer. It is the user's// responsibility to provide a valid array of suffcient size!:function EvaluatePol(const A: PFloatArray; const X: TFloat; const Degree: TPolynomDegree): TFloat;var D: TPolynomDegree; XX: TFloat;begin case Degree of 0: Result := A[0]; 1: Result := A[0] + A[1]*X; 2: Result := A[0] + A[1]*X + A[2]*X*X; 3: Result := A[0] + A[1]*X + A[2]*X*X + A[3]*X*X*X; 4: Result := A[0] + A[1]*X + A[2]*X*X + A[3]*X*X*X + A[4]*X*X*X*X; else XX := X*X*X*X; Result := A[0] + A[1]*X + A[2]*X*X + A[3]*X*X*X + A[4]*XX; for D := 5 to Degree do begin XX := XX * X; Result := Result + A[D]*XX; end; end;end; { EvaluatePol }//------------------------------------------------------------------------------// Finds a root of the given polynomial in the given interval. Returns True on// success, or False on failure. For best results, there should be only one// root in the given interval, and the interval should be small. Otherwise// the function may not find any root, or may take longer to find it. For// large intervals, the MaxFTRIterations global variable may need to be// increased. A are the polynom coefficients of the polynom, X1 and X2 define// the interval (X2>X1), MaxError (>0) is the maximum allowed error in X,// Root will hold the root if it is found, and Degree is the highest power of// X in the polynom. IMPORTANT: This function does NOT check the values of// X1, X2, and MaxError, or the float array pointer. It is the user's// responsibility to provide correct values and a valid array of correct// size!:function FindRoot(const A: PFloatArray; const X1, X2, MaxError: TFloat; out Root: TFloat; const Degree: TPolynomDegree): Boolean;var X: TFloat;begin X := 0.5 * (X1 + X2); Result := FineTuneRoot(A, MaxError, X, Degree); if (Result) then begin if ((X >= X1) and (X <= X2)) then Root := X else Result := False; endend; { FindRoot }//------------------------------------------------------------------------------// Fine-tunes the proposed root of the given polynomial using the secant// variation of the newton-raphson method. Returns True on success, or False// on failure. If the function fails, Root remains unchnaged. For best// results, the suggested root should be as close as possible to the actual// root. Otherwise the function may not find the root, or may take longer to// find it. For large distances, the MaxFTRIterations global variable may// need to be increased. A are the polynom coefficients, MaxError (>0) is the// maximum allowed error in Root, and Degree is the highest power of X in the// polynom. The code is based on the explanations in// http://www.math.ubc.ca/~duncan/newtonmethod.pdf. IMPORTANT: This function// does NOT check the values of X1, X2, and MaxError, or the float array// pointer. It is the user's responsibility to provide correct values and a// valid array of correct size!:function FineTuneRoot(const A: PFloatArray; const MaxError: TFloat; var Root: TFloat; const Degree: TPolynomDegree): Boolean;const MaxCautionLevel = 2;var Xn_1, Xn, X, Fn_1, Fn: TFloat; Iteration: Cardinal; CautionLevel: Byte;begin // Assume failure: Result := False; // Init: Xn := Root; // X_n Fn := EvaluatePol(A, Xn, Degree); // F(X_n) Xn_1 := Root + 0.001; // X_n-1 Fn_1 := EvaluatePol(A, Xn_1, Degree); // F(X_n-1) CautionLevel := 0; // Try to find a better root: Iteration := 0; while (True) do begin // Count the current iteration: Inc(Iteration); // Calculate the difference in the function values: Fn_1 := Fn_1 - Fn; // If that's zero, we're an extremum, and can't continue: if (Fn_1 = 0) then Exit; // Calculate the new root: X := Xn - Fn * (Xn_1 - Xn) / Fn_1; // Fn_1 hold the the difference now // If we didn't get much of an improvement, we're done: if (Abs(X - Xn) < MaxError) then begin // Return the root we found, and confirm success: Root := X; Result := True; Exit; end; // Check if we're doing ok so far - i.e. that we're not moving further and // further away from the root, and we haven't been trying for too long: if (Iteration >= MaxFTRIterations) then Exit; // we've reached our max. better quit now. if (Abs(X - Xn) <= Abs(Xn - Xn_1)) then CautionLevel := 0 // we're getting closer. relax. else begin Inc(CautionLevel); // we're getting further away. something may be wrong. if (CautionLevel >= MaxCautionLevel) then Exit; end; // Move down the series: Xn_1 := Xn; Fn_1 := Fn; Xn := X; Fn := EvaluatePol(A, Xn, Degree); end; // If we got here, then we failed to find the root:end; { FineTuneRoot }//------------------------------------------------------------------------------// Puts the real roots of x^3 + a2*x^2 + a1*x + a0 = 0 in Roots, and returns// the number of roots found. The roots may not be unique. The code is based// mostly on the explanations from// http://www.sosmath.com/algebra/factor/fac11/fac11.html, with some// principles adapted from http://pyx.svn.sourceforge.net/viewvc/pyx/trunk/// pyx/pyx/mathutils.py?view=log. IMPORTANT: This function does NOT check the// float array pointer. It is the user's responsibility to provide a valid// array of suffcient size!:function SolveCubic(const a2, a1, a0: TFloat; const Roots: PFloatArray): Integer;const One3rd: TFloat = 1.0 / 3.0; MinusOne3rd: TFloat = -1.0 / 3.0; MaxError = 0.01; RangeEnd = 1e100; RangeStart = -RangeEnd; RootOffset = 100;var A: array [0..3] of TFloat; M3B, // Minus 3rd B, i.e. -b/3.0 I, J: TFloat;begin // We now have: x^3 + b'*x^2 + c'*x + d' = 0. Find the coefficients of the // corresponding depressed cubic, y^3 + I*y + J = 0: M3B := a2 * MinusOne3rd; // -b'/3a = -b'/3 I := a1 + a2 * M3B; // c' - b'^2/3 J := a0 + (a1 - 2.0 * M3B * M3B) * M3B; // d' + (2*b'^3 / 27) - (c'*b'/3) // Now we must find s and t such that I = 3st and J = t^3 - s^3. Then y=s-t // is a real root of the depressed equation, and x = y - b'/3 is a real root // of the original equation. To find s and t, we make use of the fact that if // u=t^3, then u^2 - J*u - I^3/27 = 0. It is more convinient to work with // I' = I/3 = st, with which the above eq. becomes u^2 - J*u - I'^3 = 0: I := I * One3rd; Result := SolveQuadratic(-J, -(I*I*I), Roots); // If a root was found, use it to find the 1st real root of our cubic // equation. Otherwise we must use other means, since we know a cubic // equation always has at least one root: if (Result > 0) then begin // Using the subtraction (-) root of the quadratic equation seems to give // better numerical stability: if (Result > 1) then Roots[0] := Roots[1]; // Find the root: if (Abs(Roots[0]) > 1e-6) then // t^3 <> 0 begin // First find t: if (Roots[0] >= 0) then J := Power(Roots[0], One3rd) else J := -Power(-Roots[0], One3rd); // Now use t to find the root: Roots[0] := I/J - J + M3B; // s-t-b'/3 end else // t^3 ~ 0 begin // t^3 = 0, so J = -s^3, and y=s, and x = y-b'/3 = s-b'/3: if (J < 0) then Roots[0] := M3B + Power(-J, One3rd) else Roots[0] := M3B - Power(J, One3rd); end; // Confirm that we found a root: Result := 1; end else begin // Set the polynom coefficients: A[3] := 1.0; A[2] := a2; A[1] := a1; A[0] := a0; // Find the extrema of the cubic. To do that, we diferentiate and solve the // resulting quadratic. If we start with x^3 + a2*x^2 + a1*x + a0 = 0, then // after differentiation we get 3*x^2 + 2*a2*x + a1 = 0. We want the // highest coefficient to be 1, so we divide by 3, and get // x^2 + (2*a2/3)*x + a1/3 = 0. Then we solve that: Result := SolveQuadratic(2*One3rd*a2, One3rd*a1, Roots); // Next we need to guess an approximate root, as a starting point for the // newton-raphson method. We know that the cubic curve can have either 0 or // 2 extrema. If there are 0 extrema, any guess would do, so we choose 0. // If there are 2 extrema, it's a bit more complicated: if (Result = 0) then Roots[0] := 0 else begin // Find the values of our cubic curve at the two extrema: I := EvaluatePol(PFloatArray(@A[0]), Roots[0], 3); J := EvaluatePol(PFloatArray(@A[0]), Roots[1], 3); // If one of the extrema is above the x axis, and the other below, the // root must be between them, so we choose the average x. Otherwise we // need to step a bit sideways in x, on the side of the curve that goes // toward the x axis. To decide which side that is, we check which // extremum is higher and whether they are both above or below the x // axis: if (I * J < 0) then Roots[0] := 0.5 * (Roots[0] + Roots[1]) else if (I < 0) then begin // Both of the extrema are below the x axis. If the left extremum is // higher, we need to go to the right (else to the left): if (I > J) then Roots[0] := Roots[1] + RootOffset else Roots[0] := Roots[0] - RootOffset; end else begin // Both of the extrema are above the x axis. If the left extremum is // higher, we need to go to the left (else to the right): if (I > J) then Roots[0] := Roots[0] - RootOffset else Roots[0] := Roots[1] + RootOffset; end; end; // Now try to fine-tune our guess: if (FineTuneRoot(PFloatArray(@A[0]), MaxError, Roots[0], 3)) then Result := 1 else Result := 0; end; // If we couldn't find a root, we must solve the entire equation the long way // (this should rarely happen, if ever). Otherwise we can proceed to finding // the other two roots, which may be either real or complex. To find them we // will divide the polynomial by (x-x1), and solve the resulting quadratic // equation: if (Result > 0) then // found a root (Result should be 1) begin // First find the coefficients of the quadratic equation. Note // that if we work with the normalzied cubic, the 1st coefficient (a) will // again be 1: I := a2 + Roots[0]; J := a1 + I * Roots[0]; // Now solve the quadratic equation putting the two roots in A and B, and if // the roots are real, try to select the smallest positive one of the // three: Result := Result + SolveQuadratic(I, J, PFloatArray(@Roots[Result])); end else // couldn't find a root - solve the hard way: Result := SolvePolynom(PFloatArray(@A[0]), Roots, RangeStart, RangeEnd, MaxError, 3);end; { SolveCubic }//------------------------------------------------------------------------------// Solves the given polynom within the specified interval, puts its real roots// within the interval in Roots, and returns the number of roots found. The// roots are each unique, and are returned in increasing order (from low to// high). For best results, the interval should be as small as possible.// Otherwise the function may take longer to find the roots, or may not find// them at all if the MaxFTRIterations global variable is set too low. For// polynoms of the 4th or lower degree, the functions SolveQartic,// SolveCubic, and SolveQuadratic are sometimes preferable. They are usually// faster and are more general (no interval required), but may be less// accurate. A combination of those functions with FineTuneRoot may prove to// be a good compromise in some cases. A are the polynom coefficients, X1 and// X2 define the interval (X2>X1), MaxError (>0) is the maximum allowed error// in the roots, and Degree is the highest power of X in the polynom.// IMPORTANT: This function does NOT check the values of X1, X2, and// MaxError, or the float array pointers. It is the user's responsibility to// provide correct values and valid arrays of correct size!:function SolvePolynom(const A, Roots: PFloatArray; const X1, X2, MaxError: TFloat; const Degree: TPolynomDegree): Integer;var DiffA, Ex: TFloatArray; D, Extrema: Integer;begin case Degree of 0: Result := 0; 1: begin Roots[0] := -A[0] / A[1]; if ((Roots[0] >= X1) and (Roots[0] <= X2)) then Result := 1 else Result := 0; end; else // Find the extrema of the polynom in the given interval: DifferentiatePol(A, PFloatArray(@DiffA[0]), Degree); Extrema := SolvePolynom(PFloatArray(@DiffA[0]), PFloatArray(@Ex[0]), X1, X2, MaxError, Degree-1); // If we didn't find any extrema, there may be only one root: if (Extrema = 0) then begin Result := 0; if (FindRoot(A, X1,X2, MaxError, Roots[0], Degree)) then Inc(Result); Exit; end; // First root: Result := 0; if (FindRoot(A, X1,Ex[0], MaxError, Roots[0], Degree)) then Inc(Result); // Middle roots: for D := 0 to Extrema-2 do if (FindRoot(A, Ex[D],Ex[D+1], MaxError, Roots[Result], Degree)) then Inc(Result); // Last root: if (FindRoot(A, Ex[Extrema-1],X2, MaxError, Roots[Result], Degree)) then Inc(Result); end;end; { SolvePolynom }//------------------------------------------------------------------------------// Puts the real roots of x^2 + a1*x + a0 = 0 in Roots, and returns the number// of roots found. The roots may not be unique. Code adapted from// http://pyx.svn.sourceforge.net/viewvc/pyx/trunk/pyx/pyx/// mathutils.py?view=log. IMPORTANT: This function does NOT check the float// array pointer. It is the user's responsibility to provide a valid array of// suffcient size!:function SolveQuadratic(const a1, a0: TFloat; const Roots: PFloatArray): Integer;var D: TFloat;begin D := a1*a1 - 4*a0; if (D >= 0) then begin D := Sqrt(D); Roots[0] := 0.5 * (-a1 + D); Roots[1] := 0.5 * (-a1 - D); Result := 2; end else Result := 0;end; { SolveQuadratic }//------------------------------------------------------------------------------// Puts the real roots of x^4 + a3*x^3 + a2*x^2 + a1*x + a0 = 0 in Roots, and// returns the number of roots found. The roots may not be unique. Code based// mostly on http://www.gamedev.net/community/forums/// topic.asp?topic_id=451048, with some principles adapted from// http://pyx.svn.sourceforge.net/viewvc/pyx/trunk/pyx/pyx/// mathutils.py?view=log. IMPORTANT: This function does NOT check the float// array pointer. It is the user's responsibility to provide a valid array of// suffcient size!:function SolveQuartic(const a3, a2, a1, a0: TFloat; const Roots: PFloatArray): Integer;const One16th: TFloat = 1.0 / 16.0; One256th: TFloat = 1.0 / 256.0;var A: array [0..3] of TFloat; I, J, K, L: TFloat; Index: Integer;begin // Here we should normalize from b, c, d, e to b/a, c/a, d/a, and e/a, such // that we get an equation of the form x^4 + b'*x^3 + c'*x^2 + d'*x + e' = 0. // However, this function assumes we already receive the parameters that way, // so we continue from there. First, find the coefficients I, J, and K of the // depressed quartic y^4 + I*y^2 + J*y + K = 0, where x = y - b'/4: K := a3*a3; L := -3.0 * K; // -3.0*b'^2 I := L*0.125 + a2; // -3*b'^2/8 + c' J := K*a3*0.125 - a3*a2*0.5 + a1; // b'^3/8 - b'*c'/2 + d' K := K*L*One256th + K*a2*One16th - a3*a1*0.25 + a0; // -3*b'^4/256 // + b'^2*c'/16 // - b'*d'/4 + e' // Solve a cubic equation z^3 + 2I*z^2 + (I^2 - 4K)*z - J^2 = 0 to find an // intermediate root: A[2] := I+I; A[1] := I*I + (-4.0 * K); A[0] := -(J*J); Result := SolveCubic(A[2], A[1], A[0], Roots); // Discard all the non-positive roots since we don't want to take a square // root of a negative number or divide by zero. Put the smallest one of the // remaining roots in slot 0: Index := 0; while (Index < Result) do begin if (Roots[Index] > 0) then begin if (Roots[Index] < Roots[0]) then Roots[0] := Roots[Index]; Inc(Index); end else begin Roots[Index] := Roots[Result-1]; Dec(Result); end; end; // If found, process it further to find the actual 4 roots of the quartic eq.: if (Result > 0) then begin // We need a more accurate root than what SolveCubic can give us, so we must // fine-tune it a bit: A[3] := 1.0; FineTuneRoot(PFloatArray(@A[0]), 1e-4, Roots[0], 3); // Prepare the parameters for our two final quadratic equations: K := Sqrt(Roots[0]); J := J/K; L := (I + Roots[0] - J) * 0.5; // Find the four roots: Result := SolveQuadratic(K, L, PFloatArray(@Roots[0])); Inc(Result, SolveQuadratic(-K, L+J, PFloatArray(@Roots[Result]))); // Adjust the roots: L := -0.25 * a3; Index := 0; while (Index < Result) do begin Roots[Index] := Roots[Index] + L; Inc(Index); end; end;end; { SolveQuartic }end.
Ray-torus intersection test and torus normals calculation:
//------------------------------------------------------------------------------procedure TTorus.Intersect(const Ray: TRay; out Intersection: TIntersection; const FilterHeuristic: TRayFilterHeuristic = nil); {override}var A: array [0..4] of TFloat; Roots: array [0..3] of TFloat; Beta, Gamma: TFloat; I, Count: Integer; LocalRay: TRay;begin // First test against the bounding sphere - maybe it will give us a quick // "no": FBoundingSphere.Intersect(Ray, Intersection); if (Intersection.Distance < 0) then Exit; // Pre-init, assuming no intersection: inherited; // Check if there is an interesction, and find where it is: with Intersection do begin // Make a copy of the ray, and transform it to the torus' local space: LocalRay.Origin := GlobalToLocal(Ray.Origin); LocalRay.Direction := GlobalToLocal(Ray.Direction, True); // Now we need to cacluate the paremeters for the quartic equation. They are // calculated as follows: // // Ro = local ray origin = (Rox, Roy, Roz); // Rd = local ray direction = (Rdx, Rdy, Rdz); // R = major radius of torus; // r = minor radius of torus; // // R' = 4*R^2 // r' = -r^2 // R" = r' - R^2 // // a = Alpha = Rd.Rd = 1; // b = Beta = 2*Ro.Rd // c = Gamma = Ro.Ro - r^2 - R^2 // = Ro.Ro + R" // // a4 = A = a^2 = 1 // a3 = B = 2ab = 2b // a2 = C = b^2 + 2ac + 4*R^2*Rdz^2 // = b^2 + 2c + R'*Rdz^2 // a1 = D = 2bc + 8*R^2*Roz*Rdz // = Bc + 2*R'*Roz*Rdz // a0 = E = c^2 + 4*R^2*Roz^2 - 4*R^2*r^2 // = c^2 + R'*(Roz^2 + r') // // where the euqation is: Ax^4 + Bx^3 + Bx^2 + Cx + D = 0, with A=1. // // Start by finding Beta and Gamma (we know Alpha is 1, since the ray // direction is normalized). While your at it, also find the square // distance from the torus center to the ray origin: Beta := 2 * VectorDotProduct(LocalRay.Origin.Vec, LocalRay.Direction.Vec); Gamma := VectorNorm(LocalRay.Origin.Vec) + FRSubtrahend; // Now find coefficients B through E (we know A is 1): A[4] := 1.0; A[3] := Beta + Beta; A[2] := Beta * Beta + Gamma + Gamma + F4MajorRSqr * LocalRay.Direction.Z * LocalRay.Direction.Z; A[1] := A[3] * Gamma + F8MajorRSqr * LocalRay.Origin.Z * LocalRay.Direction.Z; A[0] := Gamma * Gamma + F4MajorRSqr * (LocalRay.Origin.Z * LocalRay.Origin.Z + FMinusMinorRSqr); // And finally solve the quartic equation to find the distance to the // intersection. Usually just fine-tuning the distance we got from the // bounding-cylinder intersection will work, and is a bit faster. When that // fails, we resort to the full-blown quartic equation solver, since it is // more robust. Note: for some reason, the quartic solver also fails in // some situations, but luckily FineTuneRoot handles them properly: Count := SolveQuartic(A[3], A[2], A[1], A[0], PFloatArray(@Roots[0])); // left here for comparison: // Count := SolvePolynom(PFloatArray(@A[0]), PFloatArray(@Roots[0]), 0,1000, SmallNum, 4); // If a root was not found, make sure the distance is set to a negative // value. If it was found, select the smallest positive one: if (Count = 0) then Distance := -1 else begin Distance := Roots[0]; I := 1; while (I < Count) do begin if ((Roots >= 0) and ((Distance < 0) or (Roots < Distance))) then Distance := Roots; Inc(I); end; end; end;end; { TTorus.Intersect }//------------------------------------------------------------------------------procedure TTorus.CalcNormal(const Point: T3DVector; out Normal: T3DVector); {override}var LocalPoint: T3DVector; Ratio, Len: Single;begin LocalPoint := GlobalToLocal(Point); Len := LocalPoint.X * LocalPoint.X + LocalPoint.Y * LocalPoint.Y; if (Len = 0) then begin Normal.X := 0; Normal.Y := 0; end else begin Ratio := 1 - FMajorR / Sqrt(Len); Normal.X := LocalPoint.X * Ratio; Normal.Y := LocalPoint.Y * Ratio; end; Normal.Z := LocalPoint.Z; Normal := LocalToGlobal(Normal, True); NormalizeVector(Normal.Vec);end; { TTorus.CalcNormal }
Also, it seems to be flat, it's dimensions seem to depend on the viewpoint rather than the ones I specify, and it pretty much ignores the orientation I give it, using just the global coordinates system (even though I've used the exact same transformation code before and I know it works)...
I haven't yet finished going over all my code, but it seems that at least the solver part is correct - or at least it matches what was done in the aforementioned thread.. Although to be honest I'm not 100% sure about the the cubic equation solver method, especially the case of RR > QQQ - and that is where most if not all of the rays end up. My code follows, if anyone cares to take a look, though be warned - it's in pascal, not c++. And in the mean time, anyone has any ideas what may be causing this behavior?
