Retsam wrote:
>You know, this has been bugging me, because it just seems like there's got
>to be a way to do it with straight vector math, without resorting to
>sin/cos parameterizations and taking partial derivatives (with respect
>to each parameter) to find local minima.
>
>The best I've been able to come up with so far is this...


Okay, I think I may have a fast and much more accurate approach than most of
the others mentioned.  First of all, it will be a numerical approximation
approach, so no derivatives.

Instead of finding parameterizations of the two circles, and going through
all 360 degrees of both circles and comparing the distances (which leads to
360*360 comparisons), we will simply pick an "arbitrary" point P1 on circle
C, and find the closest point on circle D to that point; call it Q1.  Then
we find the closest point on circle C to Q1, and call that P2.  Repeat
until desired accuracy is found.

So how do we pick our arbitrary point P1.  I say just pick the center, C.
There's a complicated reason why point C is a really good point to start
with.  And yes, I realize that technically it's not on the circle.  Don't
worry, it still works.  The only drawback is it throws off your first
distance measurement, if you care to take it.

Anyway, once we have a point on one circle, to find the closest point on the
other circle to that point, we simply do the following.

Given point P1 (C in this case), find the vector U1 = P1-D.  Assuming N is a
unit normal, make U1 perpendicular to N as follows:
U1 = U1 - N*(U1.N)

But we need U1 to have a length of R2 (we're on circle D), so Q1 = U1 * R2 /
sqrt(U1.U1).  I know, a pesky square root.  in POV-Ray, vlength(U1) =
sqrt(U1.U1).  Oh, and of course, we have to make sure that Q1 is centered
at D, so Q1 = Q1 + D.

Now rinse and repeat :)

//arbitrary circles for testing.
#declare C = <-1, 6, 3>;
#declare Rot1 = <25, 15, 50>;
#declare R1 = 1.2;
#declare M = vrotate(y, Rot1);

#declare D = <6, 3, 4>;
#declare Rot2 = <-55, 0, 30>;
#declare R2 = 1.5;
#declare N = vrotate(y, Rot2);

//Here's the formula itself
#declare Iters=5;
#declare P1 = C;
#while (Iters>0)
  #declare U1 = P1-D;
  #declare U1 = U1 - N*vdot(U1, N);
  #declare Q1 = D + U1 * R2 / vlength(U1);
  #debug concat("Distance is: ", str(vlength(Q1-P1), 5, 9),"\n")

  #declare V1 = Q1-C;
  #declare V1 = V1 - M*vdot(V1, M);
  #declare P1 = C + V1 * R1 / vlength(V1);
  #debug concat("Distance is: ", str(vlength(P1-Q1), 5, 9),"\n")


  #declare Iters = Iters-1;
#end

Notice with the debug messages that three or four iterations should be
plenty, but five or six should give you very, very good accuracy.  In the
sample above, on my machine the accuracy maxes out at 10 significant digits
after the first part of the fifth iteration.

Like any of the formulas, there are probably scenarios that break this
formula, causing divide by zero errors or whatever.  But in general, it is
very, very fast.

Post a reply to this message