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

>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.
>

Well, some combinations may require ten or fifteen iterations: if the
circles come close to touching, but at very oblique angles.

Post a reply to this message

Retsam wrote:
> 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.

Thanks for the code, that's pretty much the approach I was looking at, 
but hadn't got around to implementing it yet.

Post a reply to this message

Edward Coffey wrote:
>
>Thanks for the code, that's pretty much the approach I was looking at,
>but hadn't got around to implementing it yet.
>

I said there were probably scenarios that break my code, and I found one.
Take two circles, where the plane containing points C and D and vector M
also contains or very nearly contains vector N.  I.e. N ~= M + a*(D-C),
where a is arbitrary.  The easiest example is if M=N.

Now assume that the distance between C and D is less than (R1+R2), and that
the two circles are in very close parallel planes.  If this is true, then
the first point found by my algorithm, Q1, will lie roughly on a line
between C and D.  But the closest points (two of them) will be where the
circles would intersect if they were in the same plane.

If you moved point Q1 just a degree to the left or right, the algorithm
would eventually get one of the closest points.  But in this example, it
would bounce back and forth between two points on the line between the two
centers.

The only way I can think of to efficiently get around this is to make three
tests.  Start with a point between the centers C and D, and bail out when
the change in the distance is less than 1 part in 10^6 or 10^8 or whatever.

Then try two other starting points, at 90 degree angles to the original
starting point.  The easiest way to find these points would be to find:

Q1 = D-(Nx(Q1-D)) and Q1 = D+(Nx(Q1-D))

in this case, x means cross product.  Since N is a unit vector, the new
vector should still have a length of R2 (should... I haven't tried it yet).

These other two starting points should move much more quickly toward an
answer.

Also, since there can be two closest points, or two almost equally close
points, the original method might find a "local minimum" that is actually
not the minimum, so using this three point approach (or hell, just do the
two 90 degree off-center points, that might be even better) should
guarantee that you get the "most closest" point.

Post a reply to this message

Retsam wrote:
> Edward Coffey wrote:
> 
>>Thanks for the code, that's pretty much the approach I was looking at,
>>but hadn't got around to implementing it yet.
>>
> 
> 
> I said there were probably scenarios that break my code, and I found one.
> Take two circles, where the plane containing points C and D and vector M
> also contains or very nearly contains vector N.  I.e. N ~= M + a*(D-C),
> where a is arbitrary.  The easiest example is if M=N.

Thankfully my circles are very well behaved. They never touch one 
another, never lie in parallel planes and always have exactly one pair 
of closest points.

Post a reply to this message