Sign in to follow this  
technobot2

[SOLVED+UPDATE] Help debugging a ray-torus intersection (was "So I've got mysel [...]

Recommended Posts

(This mostly boils down to a lot of math, so I guess it belongs here) I just implemented ray-torus intersections for my ray tracer, taking the quartic equation derivation from this paper, and adopting the equation solver algorithm from this thread. Well, so I've got myself a torus... only it's.. uhm.. a bit large [grin]: 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?
//------------------------------------------------------------------------------
// Solves the given quadratic equation Ax^2 + Bx + C = 0, puts the two roots
//  in Root1 and Root2, and retruns True if it found any valid roots:
//
function TTorus.SolveQuadratic(A, B, C: Single; out Root1,
  Root2: Single): Boolean;
// Code partially adopted and ported from:
//  http://www.gamedev.net/community/forums/topic.asp?topic_id=451048
var
  D: Single;
begin
  // Assume we'll find a root:
  Result := True;

  // If a is negligible, we have a linear equation bx + c = 0:

  if ((A = 0) or (Abs(A/B) < SmallNum)) then
  begin
    // If c=0, x must be 0:
    if (C = 0) then
    begin
      Root1 := 0;
      Root2 := 0;
      Exit;
    end;

    // If b=0 (and c<>0), there is no solution:
    if (Abs(B) < SmallNum) then
    begin
      Result := False;
      Exit;
    end;

    // Otherwise, the solution is as usual:
    Root1 := -C / B;
    Root2 := Root1;
    Exit;
  end;

  // Calculate the delta d = (b^2 - 4*a*c):
  D := B * B - 4 * A * C;

  // If the delta is negative, there are no real roots:
  if (D < 0) then
  begin
    Result := False;
    Exit;
  end;

  // There is at least one root:

  A := -0.5 / A;
  if (D > 0) then
  begin
    D := Sqrt(D);
    Root1 := (B + D) * A;
    Root2 := (B - D) * A;
  end
  else // d=0
  begin
    Root1 := B * A;
    Root2 := Root1;
  end;
end; { TTorus.SolveQuadratic }




//------------------------------------------------------------------------------
// Solves the given cubic equation Ax^3 + Bx^2 + Cx + D = 0, puts the
//  smallest non-negative root in Root (if possible), and retruns True if it
//  found any valid roots:
//
function TTorus.SolveCubic(A, B, C, D: Single; out Root: Single): Boolean;
// Code adopted (and adapted) and ported from:
//  http://www.gamedev.net/community/forums/topic.asp?topic_id=451048
const
  One3rd: Single = 1.0 / 3.0;
  MinusOne3rd: Single = -1.0 / 3.0;
  One9th: Single = 1.0 / 9.0;
  One54th: Single = 1.0 / 54.0;
var
  R, RR, Q, QQQ, Theta: Single;
begin
  // If A is negligible, the equation reduces to a quadratic one:
  if ((A = 0) or (Abs(A/B) < SmallNum)) then
  begin
    Result := SolveQuadratic(B, C, D, Root, A);
    if (Result) then
      if ((A >= 0) and (A < Root)) then
        Root := A;
    Exit;
  end;

  // Normalize the equation - convert b, c, d to b/a, c/a, d/a, respectively:
  if (not Equal(A, 1)) then
  begin
    A := 1 / A;
    B := B * A;
    C := C * A;
    D := D * A;
  end;

  // Compute the two factors Q^3 and R^2:
  RR := B*B;
  Q := (RR - C * 3.0) * One9th;                  // (B^2 - 3C)/9
  QQQ := Q*Q*Q;
  R := (2.0*RR*B - 9.0*B*C + 27.0*D) * One54th;  // (2B^3 - 9BC + 27D) / 54
  RR := R*R;

  // If R^2 < Q^3 , we have 3 real roots, else we have 1:

  if (RR < QQQ) then
  begin
    // Find Theta. Note that the sqrt and division are safe, since RR >= 0 and
    //  QQQ > RR, so QQQ > 0. The acos is also safe, since RR/QQQ < 1, and thus
    //  R/sqrt(QQQ) = sqrt(RR/QQQ) < 1:
    Theta := ArcCos(R / Sqrt(QQQ));

    // Prepare for finding the 3 roots. Note that the following sqrt is safe,
    //  since QQQ >= 0, and thus Q >= 0:
    A := -2.0 * Sqrt(Q);
    B := B * MinusOne3rd;

    // Find the 3 roots::
    Root := A * Cos(Theta * One3rd) + B;
    R := A * Cos((Theta + c2PI) * One3rd) + B;
    RR := A * Cos((Theta - c2PI) * One3rd) + B;

    // Select the smallest non-negative root (or just leave the 1st one, if they're
    //  all negative):
    if ((R >= 0) and (R < Root)) then
      Root := R;
    if ((RR >= 0) and (RR < Root)) then
      Root := RR;
  end
  else // RR >= QQQ - 1 real root
  begin

    /////////////////////////////////////////////////
    //                                             // 
    //  Most rays go here, but is this correct???  //
    //                                             // 
    /////////////////////////////////////////////////

    A := -Power(Abs(R) + Sqrt(RR - QQQ), One3rd);
    if (A <> 0) then
    begin
      if (R < 0) then
        A := -A;
      Root := A + (Q / A) + (B * MinusOne3rd);
    end
    else
      Root := B * MinusOne3rd;
  end;

  // Confirm that we found a valid root:
  Result := True;
end; { TTorus.SolveCubic }




//------------------------------------------------------------------------------
// Solves the given quartic equation x^4 + Bx^3 + Cx^2 + Dx + E = 0 (A=1),
//  puts the smallest non-negative root in Root, and retruns True if it
//  found any valid roots:
//
function TTorus.SolveQuartic(B, C, D, E: Single; out Root: Single): Boolean;
// Code adopted (and adapted) and ported from:
//  http://www.gamedev.net/community/forums/topic.asp?topic_id=451048
const
  OneOverSmallNum = 1.0 / SmallNum;
  One16th: Single = 1.0 / 16.0;
  One256th: Single = 1.0 / 256.0;
var
  I, J, K, L: Single;
  Found: Boolean;
begin
  // If a is negligible compared to the others, treat as a cubic equation. Note
  //  that this function assumes that a=1:
  if ((B > OneOverSmallNum) or (C > OneOverSmallNum) or
    (D > OneOverSmallNum)) then
  begin
    Result := SolveCubic(B, C, D, E, Root);
    Exit;
  end;

  // Here we should normalize from b, c, d, e to b/a, c/a, d/a, and e/a, but
  //  this function assumes we already receive the parameters that way. So we
  //  continue from there. First, find the depressed quartic coefficients, I, J,
  //  and K:

  K := B*B;
  L := -3.0 * K;      // -3.0*B*B
  
  I := L*0.125 + C;
  J := K*B*0.125 - B*C*0.5 + D;
  K := K*L*One256th + K*C*One16th - B*D*0.25 + E;

  // Solve a cubic equation z^3 + 2I*z^2 + (I^2 - 4K)*z - J^2 = 0 to find an
  //  intermediate root:
  Result := SolveCubic(1.0, I+I, I*I + (-4.0 * K), -(J*J), Root);

  // If found, process it further to find the actual 4 roots of the quartic eq.:
  if (Result) then
  begin
    // Prepare some parameters:
    B := 0.25 * B;
    K := Sqrt(Root);
    J := J/K;
    L := (I + Root - J) * 0.5;

    // Find the four roots, and put them into Root, C, D, E:
    Result := SolveQuadratic(1.0, K, L, Root, C);
    Found := SolveQuadratic(1.0, -K, L+J, D, E);

    // Select the root that will become the smallest non-negative one after the
    //  final adjustment:

    if (Result) then
      if ((C >= B) and (C < Root)) then
        Root := C;
    if (Found) then
    begin
      if ((D >= B) and (D < Root)) then
        Root := D;
      if ((E >= B) and (E < Root)) then
        Root := E;
    end;

    // Adjust the selected root:
    Root := Root - B;
    // Check if we found a non-negative root:
    Result := Result or Found;
  end;
end; { TTorus.SolveQuartic }




//------------------------------------------------------------------------------
procedure TTorus.Intersect(const Ray: TRay; out Intersection: TIntersection);
  {override}
var
  LocalRay: TRay;
  Beta, Gamma, B, C, D, E: Single;
  F4MajorRSqr, F8MajorRSqr, FMinusMinorRSqr, FRSubtrahend: Single;
begin
  // 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 quadric 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"
    //
    //  A = a^2 = 1
    //  B = 2ab = 2b
    //  C = b^2 + 2ac + 4*R^2*Rdz^2
    //    = b^2 + 2c + R'*Rdz^2
    //  D = 2bc + 8*R^2*Roz*Rdz
    //    = Bc + 2*R'*Roz*Rdz
    //  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.

    // These will be precalculated later:
    F4MajorRSqr := FMajorR * FMajorR;
    FMinusMinorRSqr := -(FMinorR * FMinorR);
    FRSubtrahend := FMinusMinorRSqr - F4MajorRSqr;
    F4MajorRSqr := 4 * F4MajorRSqr;
    F8MajorRSqr := F4MajorRSqr + F4MajorRSqr;

    // Start by finding Beta and Gamma (we know Alpha is 1, since the ray
    //  direction is normalized):
    Beta := 2 * VectorDotProduct(LocalRay.Origin.Vec, LocalRay.Direction.Vec);
    Gamma := VectorNorm(LocalRay.Origin.Vec) + FRSubtrahend;

    // Now find coefficients B thourgh E (we know A is 1):
    B := Beta + Beta;
    C := Beta * Beta + Gamma + Gamma + F4MajorRSqr * LocalRay.Direction.Z *
      LocalRay.Direction.Z;
    D := B * Gamma + F8MajorRSqr * LocalRay.Origin.Z * LocalRay.Direction.Z;
    E := Gamma * Gamma + F4MajorRSqr * (LocalRay.Origin.Z * LocalRay.Origin.Z +
      FMinusMinorRSqr);

    // And finally solve the quartic equation to find the distance to the
    //  intersection. If a root is not found, make sure the distance is set to a
    //  negative value (i.e. no intersection):
    if (not SolveQuartic(B, C, D, E, Distance)) then
      Distance := -1;
  end;
end; { TTorus.Intersect }



[Edited by - technobot2 on July 22, 2007 2:35:30 PM]

Share this post


Link to post
Share on other sites
UPDATE: I finished going over the code. Didn't find any mistakes relative to the stated sources. Furthermore, also confirmed the quartic equation derivation from additional sources. It seems whatever goes wrong is probably in the root calculations, I suspect most likely in the SolveCubic method. But I'm still open to other suggestions.

I've written an alternative SolveCubic method based on the explanations here - the code is listed below. However this also doesn't give me correct results, perhaps because it's not enough for it to return just any root... The result with this new method is that I don't see the torus at all, except for very specific viewpoints, and I get very strange shadows.

The new SolveCubic method:

//------------------------------------------------------------------------------
// Solves the given cubic equation Ax^3 + Bx^2 + Cx + D = 0, puts an
// arbitrary real root in Root (if there is one), and retruns True if it
// found any valid roots:
//
function TTorus.SolveCubic(A, B, C, D: Single; out Root: Single): Boolean;
// Code based on solution method from
// http://www.sosmath.com/algebra/factor/fac11/fac11.html
const
One3rd: Single = 1.0 / 3.0;
MinusOne3rd: Single = -1.0 / 3.0;
var
I, J: Single;
begin
// If A is negligible, the equation reduces to a quadratic one:
if ((A = 0) or (Abs(A/B) < SmallNum)) then
begin
Result := SolveQuadratic(B, C, D, Root, A);
if (Result) then
if ((A >= 0) and (A < Root)) then
Root := A;
Exit;
end;

// Normalize the equation - convert b, c, d to b'=b/a, c'=c/a, d'=d/a,
// respectively:
if (not Equal(A, 1)) then
begin
A := 1 / A;
B := B * A;
C := C * A;
D := D * A;
end;

// 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:

A := B * MinusOne3rd; // -b'/3
I := C + B*A; // c' - b'^2/3
J := D + (C - 2.0*A*A)*A; // 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(1.0, -J, I*I*I, Root, J);

// If a root was found, use it to find a root for the original equation:
if (Result) then
begin
J := Power(Root, One3rd); // t
Root := I/J - J + A; // s-t-b'/3 (reminder: A currently holds -b'/3)
end;
end; { TTorus.SolveCubic }

Share this post


Link to post
Share on other sites
After implementing that new cubic equation solver, and fixing a problem with taking the 3rd root of negative numbers, I seem to be getting somewhere.. Now instead of a huge flat torus, I'm getting what seems to be more of a large peanut - at least it's not flat anymore. I'm guessing somehow the torus surface kinda "bends the other way". So the question now is - what should be the error in the equation(s) to generate this effect?

This is roughly what seems to be happening atm:

____
/ \ <- peanut - what I get
___/ \___
/ \ / \
| | | | <- torus - what I want
\___/ \___/
\ /
\____/

Share this post


Link to post
Share on other sites
Thanks, Paul. I've already been there - that's one of the "additional sources" I mentioned in my 2nd post. As far I can tell, I'm doing the same thing. It's more likely that the polynom solvers are somehow distorting the result. And since the end result seems to be peanut-like, I'm thinking that if someone can tell me the mathematical equation for a peanut (that sounds funny [lol]), maybe it will help me figure out where the problem spot is... In the mean time, I've found another implementation of a 4th degree polynom solver, in python. I'll try to port it a bit later, and see if that makes any difference.

Share this post


Link to post
Share on other sites
The porting effort has paid off! I am finally getting a proper torus that responds correctly to all parameters! [cool] Albeit, there are some major numerical stability issues, which I still need to iron out... Here is the result so far (without lighting, since my normals calculation is still somewhat off):



And here is the current solver code, still very crude at this point:


// Adopoted and ported from:
// http://pyx.svn.sourceforge.net/viewvc/pyx/trunk/pyx/pyx/mathutils.py?view=log


//------------------------------------------------------------------------------
// 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:
function SolveQuadratic(const a1, a0: Single;
const Roots: PSingleArray): Integer;
var
D: Single;
begin
if (not Assigned(Roots)) then
raise Exception.Create('Roots array pointer not assigned');

D := a1*a1 - 4*a0;
if (D < 0) then
begin
Result := 0;
Exit;
end;
D := Sqrt(D);
Roots[0] := 0.5 * (-a1 + D);
Roots[1] := 0.5 * (-a1 - D);
Result := 2;
end; { SolveQuadratic }


//------------------------------------------------------------------------------
// 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:
function SolveCubic(const a2, a1, a0: Single;
const Roots: PSingleArray): Integer;
var
Q, R, S, T, D: Single;
begin
if (not Assigned(Roots)) then
raise Exception.Create('Roots array pointer not assigned');

// see http://mathworld.wolfram.com/CubicFormula.html for details

Q := (3*a1 - a2*a2) / 9.0;
R := (9*a2*a1 - 27*a0 - 2*a2*a2*a2) / 54.0;
D := Q*Q*Q + R*R;

if (D > 0) then // one real and two complex roots
begin
D := Sqrt(D);
if (R + D >= 0) then
S := Power(R + D, 1/3.0)
else
S := -Power(-R - D, 1/3.0);
if (R - D >= 0) then
T := Power(R - D, 1/3.0)
else
T := -Power(D - R, 1/3.0);
Roots[0] := S + T - a2/3.0;
Result := 1;
end
else if (D = 0) then
begin
if (Q = 0) then // one real root (R=0)
begin
Roots[0] := -a2/3.0;
Result := 1;
end
else // two real roots (R>0, Q<0)
begin
S := -Sqrt(-Q);
Roots[0] := 2*S - a2/3.0;
Roots[1] := -S - a2/3.0;
Result := 2;
end;
end
else // three real roots (Q<0)
begin
Q := Sqrt(-Q);
D := R / (Q*Q*Q);
if (D >= 1) then
T := 0 // Theta
else if (D <= -1) then
T := Pi
else
T := ArcCos(D);
Roots[0] := 2 * Q * Cos(T/3.0) - a2/3.0;
Roots[0] := 2 * Q * Cos((T + 4*Pi)/3.0) - a2/3.0;
Roots[0] := 2 * Q * Cos((T + 8*Pi)/3.0) - a2/3.0;
Result := 3;
end;
end; { SolveCubic }


//------------------------------------------------------------------------------
// 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:
function SolveQuartic(const a3, a2, a1, a0: Single;
const Roots: PSingleArray): Integer;
var
Ys: array [0..2] of Single;
I, MinY: Integer;
begin
if (not Assigned(Roots)) then
raise Exception.Create('Roots array pointer not assigned');

// see http://mathworld.wolfram.com/QuarticEquation.html for details

// Get up to three intermediate roots:
Result := SolveCubic(-a2, a1*a3 - 4*a0, 4*a0*a2 - a1*a1 - a0*a3*a3,
PSingleArray(@Ys[0]));
// Eliminate all invalid roots, and find the smalest of the remaining ones:
I := 0;
MinY := 0;
while (I < Result) do
begin
if ((a3*a3 - 4*a2 + 4*Ys[I] >= 0) and (Ys[I]*Ys[I] - 4*a0 >= 0)) then
begin
if (Ys[I] < Ys[MinY]) then
MinY := I;
Inc(I);
end
else
begin
Ys[I] := Ys[Result-1];
Dec(Result);
end;
end;

// Contiue as in original:
if (Result = 0) then
Exit;
if (a3*Ys[MinY] - 2*a1 < 0) then
begin
Result := SolveQuadratic(0.5 * (a3 + Sqrt(a3*a3 - 4*a2 + 4*Ys[MinY])),
0.5 * (Ys[MinY] - Sqrt(Ys[MinY]*Ys[MinY] - 4*a0)),
PSingleArray(@Roots[0]));
Inc(Result, SolveQuadratic(0.5 * (a3 - Sqrt(a3*a3 - 4*a2 + 4*Ys[MinY])),
0.5 * (Ys[MinY] + Sqrt(Ys[MinY]*Ys[MinY] - 4*a0)),
PSingleArray(@Roots[Result])));
end
else
begin
Result := SolveQuadratic(0.5 * (a3 + Sqrt(a3*a3 - 4*a2 + 4*Ys[MinY])),
0.5 * (Ys[MinY] + Sqrt(Ys[MinY]*Ys[MinY] - 4*a0)),
PSingleArray(@Roots[0]));
Inc(Result, SolveQuadratic(0.5 * (a3 - Sqrt(a3*a3 - 4*a2 + 4*Ys[MinY])),
0.5 * (Ys[MinY] - Sqrt(Ys[MinY]*Ys[MinY] - 4*a0)),
PSingleArray(@Roots[Result])));
end;
end; { SolveQuartic }

Share this post


Link to post
Share on other sites
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 . }
{ }
{******************************************************************************}

interface


const
// 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;
end
end; { 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[I] >= 0) and ((Distance < 0) or (Roots[I] < Distance))) then
Distance := Roots[I];
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 }

Share this post


Link to post
Share on other sites
Shame no screen shot of the peanut. Also i would i like to see the mathmatical method for generating the perfect peanut.

Good thread. Enjoyed reading it (even if i didnt understand a lot of it!)

Share this post


Link to post
Share on other sites
Pulled out some of the archived code, so here's that peanut :-) - it's without lighting, coz there were too many shadowing artifacts with the lighting on:



But these results were due to a bug in the polynom solvers, so I don't have a peanut equation for you. If you're interested, maybe it deserves a separate thread? :-)

Share this post


Link to post
Share on other sites

Create an account or sign in to comment

You need to be a member in order to leave a comment

Create an account

Sign up for a new account in our community. It's easy!

Register a new account

Sign in

Already have an account? Sign in here.

Sign In Now

Sign in to follow this