This is a problem that comes up once in a while. I'm curious if anyone knows 
a good solution.

You have two objects and circles that bound them. For each circle, you know 
its position and its radius squared. (You don't have the radius unless you 
take the square root, which is undesirable because it's slow, and caching it 
isn't a good option.) You want to know if these circles intersect.

Obviously, it's easy to get the distance between the circle's centers 
squared. What you want to test is whether dist < radius1 + radius2, but all 
you have is distSq, radius1Sq, and radius2Sq.

The best I can do is to square both sides of the test:

distSq < (radius1 + radius2)^2

And then use the fact that (a+b)^2 = a^2 + 2*a*b + b^2 to get:

distSq < radius1Sq + radius2Sq + 2 * radius1 * radius2

And then, since the radii are positive, the right hand side can be changed 
to the following, but the test will now succeed more often than it should:

distSq < radius1Sq + radius2Sq + 2 * max( radius1Sq, radius2Sq )

So now we have a test that doesn't require taking any square roots, but it 
isn't as "tight" as it could be - it will return true sometimes when the 
circles aren't intersecting.

Can this be improved?

 - Slime
 [ http://www.slimeland.com/ ]

Post a reply to this message

If you're willing to tolerate an intersection function which has the 
possibility of error, it seems like just using bounding squares instead 
of bounding circles would give better error bounds if you expect your 
circles to often be of very different sizes.

With regard to the actual question you asked, I can't think of a good 
way without square roots off of the top of my head, but if you're 
willing to tolerate a bit of approximation error then computing them 
needn't be slow:

http://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Approximations_that_depend_on_IEEE_representation

Post a reply to this message

> 
> distSq < radius1Sq + radius2Sq + 2 * max( radius1Sq, radius2Sq )
> 
> So now we have a test that doesn't require taking any square roots, but it 
> isn't as "tight" as it could be - it will return true sometimes when the 
> circles aren't intersecting.
> 
> Can this be improved?
> 
>  - Slime
>  [ http://www.slimeland.com/ ] 
> 
> 

If the loose test indicates possible intersection then perform the more 
expensive exact test using Sqrt().  If you can eliminate the majority of 
cases cheaply then the higher cost is amortised.

In a similar case but where I have radius1, radius2 I do the following 
to reduce the number of cases requiring more arithmetic operations:

float rSum = radius1 + radius2;

float xDiff = x1 - x2;
if (xDiff > rSum) return false;

float yDiff = y1 - y2;
if (yDiff > rSum) return false;

if (xDiff + yDiff) < rSum return true;

return (xDiff * xDiff + yDiff * yDiff < rSum * rsum)

Sqrt() isn't needed at all but of course this is a different scenario to 
yours.

As always, you would need to measure the performance of different 
schemes on your actual data and processor.  Too many condition tests and 
intermediate operations can end up being more expensive than Sqrt() 
given modern FP capabilities.

Post a reply to this message

Slime wrote:

> The best I can do is to square both sides of the test:
> 
> distSq < (radius1 + radius2)^2
> 
> And then use the fact that (a+b)^2 = a^2 + 2*a*b + b^2 to get:
> 
> distSq < radius1Sq + radius2Sq + 2 * radius1 * radius2
> 
> And then, since the radii are positive, the right hand side can be changed 
> to the following, but the test will now succeed more often than it should:
> 
> distSq < radius1Sq + radius2Sq + 2 * max( radius1Sq, radius2Sq )
> 
> So now we have a test that doesn't require taking any square roots, but it 
> isn't as "tight" as it could be - it will return true sometimes when the 
> circles aren't intersecting.
> 
> Can this be improved?

If you're that hell-bent on avoiding the square root, it can be 
approximated as follows:

   Suppose that the most significant non-zero digit of x is in the 2^n 
place. Let y = x / 2^(n/2). Now Sqrt(x) is approximately x / y.

Usually this estimated value is greater than the true square root (but 
always by less than 45%). However, exact powers of 2 seem to come out 
slightly below the true value. (By about 10^-16.) I don't know if that's 
a glitch in my test program or what... It shouldn't be insummountable 
though.

Post a reply to this message

> distSq < radius1Sq + radius2Sq + 2 * max( radius1Sq, radius2Sq )
>
> So now we have a test that doesn't require taking any square roots, but it 
> isn't as "tight" as it could be - it will return true sometimes when the 
> circles aren't intersecting.
>
> Can this be improved?

How about more simply:

distSq < 4 * max( radius1Sq, radius2Sq )

Followed by the expensive square root operation to check if they really 
intersect?

Obviously the benefits will depend on many factors, like how similarly sized 
your circles are, and what % of them intersect.

Is there any history to the comparisons that you could make use of?  Eg a 
large number of circles that are being repeatedly tested for intersection 
after they move a small amount?  If so you can get much bigger speed-ups 
using other methods (eg keeping the list of circles ordered by x 
coordinate).

Post a reply to this message

scott wrote:

> Obviously the benefits will depend on many factors, like how similarly 
> sized your circles are, and what % of them intersect.

I would imagine it's quite possible that the extra time required to 
perform a square root operation is more than made up for by the decrease 
in unecessary work due to false intersections. Square root isn't *that* 
slow... But you'd have to measure it to find out for sure.

Post a reply to this message

On 12/10/2009 04:34 AM, Slime wrote:

> And then, since the radii are positive, the right hand side can be changed
> to the following, but the test will now succeed more often than it should:
>
> distSq<  radius1Sq + radius2Sq + 2 * max( radius1Sq, radius2Sq )
>
> So now we have a test that doesn't require taking any square roots, but it
> isn't as "tight" as it could be - it will return true sometimes when the
> circles aren't intersecting.
>
> Can this be improved?
>

Have you tried this condition?

((distSq < (radius1Sq+radius2Sq))
or (4*radius1Sq*radius2Sq > (distSq - (radius1Sq+radius2Sq))^2))

Post a reply to this message

I wrote a small bit of code to test what the speed of the various 
approximation methods was.  The code was complied with g++ and the times 
given for each method are for 100000000 calls to the function.  First, 
the functions I tested:

// Naive method
inline bool isect1(double aSq, double bSq, double dSq) {
     return sqrt(dSq) < sqrt(aSq) + sqrt(bSq);
}

// Slime's approximate method
inline bool isect2(double aSq, double bSq, double dSq) {
     return dSq < aSq + bSq + 2*std::max(aSq,bSq);
}

// Suggested by bart, but with a bugfix.
// Avoids a sqrt, but uses a few multiplies
inline bool isect3(double aSq, double bSq, double dSq) {
     double t1 = aSq+bSq;
     double t2 = dSq - t1
     return (dSq < t1) || (4*aSq*bSq > t2*t2);
}

// Naive method using a fast approximate sqrt
inline bool isect4(double aSq, double bSq, double dSq) {
     return fastSqrt(dSq) < fastSqrt(aSq) + fastSqrt(bSq);
}

and now the times under various optimization levels:


---------- O0
time for isect1 = 9.66935
time for isect2 = 2.03143
time for isect3 = 2.18224
time for isect4 = 7.03525

---------- O1
time for isect1 = 5.13777
time for isect2 = 0.413757
time for isect3 = 0.413756
time for isect4 = 0.0376559

---------- O2
time for isect1 = 0.564358
time for isect2 = 0.276049
time for isect3 = 0.305847
time for isect4 = 0.0381372

---------- O3
time for isect1 = 0.564529
time for isect2 = 0.159936
time for isect3 = 0.206905
time for isect4 = 0.0376384


It looks like my initial suggestion implemented in isect4 is the fastest 
when you have optimization turned on, by a margin which rather surprises 
me.  It takes about as long as simply computing the incorrect formula 
'dSq < aSq + bSq', which seems almost impossibly fast.  Nevertheless 
I've checked the results and it seems to work pretty well despite the 
approximation in the sqrt call.  If you don't want to live with the 
approximation, however, bart's method seems to be the best choice.

Post a reply to this message

Kevin Wampler wrote:
> // Suggested by bart, but with a bugfix.

On second glance, the bug was in how I read his math, and his equation 
was entirely correct.

Post a reply to this message

Kevin Wampler wrote:

> I wrote a small bit of code to test what the speed of the various
> approximation methods was.  The code was complied with g++ and the times
> given for each method are for 100000000 calls to the function.  First,
> the functions I tested:
> 
> // Naive method
> inline bool isect1(double aSq, double bSq, double dSq) {
>      return sqrt(dSq) < sqrt(aSq) + sqrt(bSq);
> }
> 
> // Slime's approximate method
> inline bool isect2(double aSq, double bSq, double dSq) {
>      return dSq < aSq + bSq + 2*std::max(aSq,bSq);
> }
> 
> // Suggested by bart, but with a bugfix.
> // Avoids a sqrt, but uses a few multiplies
> inline bool isect3(double aSq, double bSq, double dSq) {
>      double t1 = aSq+bSq;
>      double t2 = dSq - t1
>      return (dSq < t1) || (4*aSq*bSq > t2*t2);
> }
> 
> // Naive method using a fast approximate sqrt
> inline bool isect4(double aSq, double bSq, double dSq) {
>      return fastSqrt(dSq) < fastSqrt(aSq) + fastSqrt(bSq);
> }
> 
> and now the times under various optimization levels:
> 
> 
> ---------- O0
> time for isect1 = 9.66935
> time for isect2 = 2.03143
> time for isect3 = 2.18224
> time for isect4 = 7.03525
> 
> ---------- O1
> time for isect1 = 5.13777
> time for isect2 = 0.413757
> time for isect3 = 0.413756
> time for isect4 = 0.0376559
> 
> ---------- O2
> time for isect1 = 0.564358
> time for isect2 = 0.276049
> time for isect3 = 0.305847
> time for isect4 = 0.0381372
> 
> ---------- O3
> time for isect1 = 0.564529
> time for isect2 = 0.159936
> time for isect3 = 0.206905
> time for isect4 = 0.0376384
> 
> 
> It looks like my initial suggestion implemented in isect4 is the fastest
> when you have optimization turned on, by a margin which rather surprises
> me.  It takes about as long as simply computing the incorrect formula
> 'dSq < aSq + bSq', which seems almost impossibly fast.  Nevertheless
> I've checked the results and it seems to work pretty well despite the
> approximation in the sqrt call.  If you don't want to live with the
> approximation, however, bart's method seems to be the best choice.

How are you doing the "run 100000000 times" part?

The compiler may be optimizing knowing you're runing it on the same input 
values each time. For example, only calling sqrt() 3 times and not 
300000000.

Post a reply to this message

> Have you tried this condition?
>
> ((distSq < (radius1Sq+radius2Sq))
> or (4*radius1Sq*radius2Sq > (distSq - (radius1Sq+radius2Sq))^2))

Smart. I feel stupid now! =) I see the equivalence, but I'm curious how long 
you took to come up with that. Was it obvious?

 - Slime
 [ http://www.slimeland.com/ ]

Post a reply to this message

> How about more simply:
>
> distSq < 4 * max( radius1Sq, radius2Sq )
>
> Followed by the expensive square root operation to check if they really 
> intersect?

Well, I'd have to weigh the cost of doing the square root against the cost 
of just going forward and doing the more expensive test. Although I posed 
this as a stand-alone question, the actual work I am doing is checking if a 
line segment intersects a triangle.

There's probably a bunch of ways to make my code faster, but I became 
curious about this particular problem.

 - Slime
 [ http://www.slimeland.com/ ]

Post a reply to this message

>   Suppose that the most significant non-zero digit of x is in the 2^n 
> place. Let y = x / 2^(n/2). Now Sqrt(x) is approximately x / y.

Is there a way to quickly get n for a given number?

 - Slime
 [ http://www.slimeland.com/ ]

Post a reply to this message

> float xDiff = x1 - x2;
> if (xDiff > rSum) return false;
>
> float yDiff = y1 - y2;
> if (yDiff > rSum) return false;
>
> if (xDiff + yDiff) < rSum return true;
>
> return (xDiff * xDiff + yDiff * yDiff < rSum * rsum)

I would want to make sure the additional branches don't make this slower. 
You also have to take the absolute value of xDiff and yDiff for the first 
two checks, which may add slowness. However, it may be worth it if you know 
it's extremely uncommon for the circles to overlap in the X or Y direction 
(which I suppose is a reasonable assumption).

 - Slime
 [ http://www.slimeland.com/ ]

Post a reply to this message

> If you're willing to tolerate an intersection function which has the 
> possibility of error, it seems like just using bounding squares instead of 
> bounding circles would give better error bounds if you expect your circles 
> to often be of very different sizes.

I generally prefer bounding circle/sphere checks to bounding box checks 
because box intersection tests (at least naive ones) require multiple 
branches. Also, in my particular application I'm dealing with a lot of line 
segments, for which I have to create the bounding boxes, which I think 
requires at least two more branches (to check if x1 < x2 and if y1 < y2). Of 
course I'd have to try it to be sure.

 - Slime
 [ http://www.slimeland.com/ ]

Post a reply to this message

Thanks for doing all that work! =) Definitely looks like Bart's method is 
the best. I would be scared of using an approximation in case it failed to 
determine an important intersection. What is your implementation of 
fastSqrt?

 - Slime
 [ http://www.slimeland.com/ ]

Post a reply to this message

http://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Approximations_that_depend_on_IEEE_representation

Craziness. I remember finding that invSqrt function a few years ago and 
looking it up. I'm always a little worried about portability of this sort of 
thing, though.

 - Slime
 [ http://www.slimeland.com/ ]

Post a reply to this message

Nicolas Alvarez wrote:
> The compiler may be optimizing knowing you're runing it on the same input 
> values each time. For example, only calling sqrt() 3 times and not 
> 300000000.


Ahh, good point.  After fixing this it looks like isect4 is faster than 
isect3, but by a pretty slim margin, so bart's method would definitely 
be the way to go.

Post a reply to this message

Slime wrote:
> Thanks for doing all that work! =) Definitely looks like Bart's method is 
> the best. I would be scared of using an approximation in case it failed to 
> determine an important intersection. What is your implementation of 
> fastSqrt?

I just used the one given on Wikipedia.

Post a reply to this message

> Although I posed this as a stand-alone question, the actual work I am 
> doing is checking if a line segment intersects a triangle.

So let me understand this, you are constructing two circles, one around the 
line and one around the triangle, then seeing the circles intersect, and if 
they do proceeding to check if the line actually intersects the triangle?

Wouldn't it be faster to just directly check the line against the triangle?

Post a reply to this message

>>   Suppose that the most significant non-zero digit of x is in the 2^n 
>> place. Let y = x / 2^(n/2). Now Sqrt(x) is approximately x / y.
> 
> Is there a way to quickly get n for a given number?

If you don't mind depending on a specific floating-point 
representation... sure. Just find out which bits of the number contain 
the exponent. (For example, for a double-precision IEEE float, the 
bitmask is 0x7FFF000000000000, and you then shift out those zeros and 
subtract 0x3FFF to get the real value.)

That's what makes this so efficient. You just need the exponent field 
from the float value, so no need to calculate anything. (Indeed, as an 
even rougher approximation, just halve the exponent...)

Post a reply to this message

scott wrote:

> So let me understand this, you are constructing two circles, one around 
> the line and one around the triangle, then seeing the circles intersect, 
> and if they do proceeding to check if the line actually intersects the 
> triangle?
> 
> Wouldn't it be faster to just directly check the line against the triangle?

It would.

Now, if the circle encloses a few *thousand* triangles, it's probably 
way, way faster to check the circle first...

Post a reply to this message

> So let me understand this, you are constructing two circles, one around 
> the line and one around the triangle, then seeing the circles intersect, 
> and if they do proceeding to check if the line actually intersects the 
> triangle?

Well, the one around the triangle is precomputed. Otherwise, yes.

> Wouldn't it be faster to just directly check the line against the 
> triangle?

Honestly, I don't know. I'm not sure my code for the triangle-line check is 
very good. =)

As far as I can tell, to check for an intersection between a triangle and a 
line segment, you need to check that one of the triangle's vertices is on an 
opposite side of the line from the other two, and then check that both of 
the line segment's endpoints are on opposite sides of at least one of the 
triangle's edges, or that both of them are inside the triangle.

I did most of these checks with cross products, and then found I had some 
precision problems in some cases, and the code got more complicated when I 
tried to avoid them. Since intersections are rare, I figured a bounding 
check would be the easiest way to improve performance.

 - Slime
 [ http://www.slimeland.com/ ]

Post a reply to this message

>> Wouldn't it be faster to just directly check the line against the 
>> triangle?
> 
> Honestly, I don't know. I'm not sure my code for the triangle-line check is 
> very good. =)
> 
> As far as I can tell, to check for an intersection between a triangle and a 
> line segment, you need to check that one of the triangle's vertices is on an 
> opposite side of the line from the other two, and then check that both of 
> the line segment's endpoints are on opposite sides of at least one of the 
> triangle's edges, or that both of them are inside the triangle.
> 
> I did most of these checks with cross products, and then found I had some 
> precision problems in some cases, and the code got more complicated when I 
> tried to avoid them. Since intersections are rare, I figured a bounding 
> check would be the easiest way to improve performance.

Er... yeah, that's not terribly good. ;-)

Do a quick google for "ray triangle intersection". There are several 
methods for doing this kind of thing, some of them highly optimised.

Post a reply to this message

> Er... yeah, that's not terribly good. ;-)
>
> Do a quick google for "ray triangle intersection". There are several 
> methods for doing this kind of thing, some of them highly optimised.

That tends to give results for a 3D problem, whereas my problem is 2D.

 - Slime
 [ http://www.slimeland.com/ ]

Post a reply to this message

Slime wrote:
>> Er... yeah, that's not terribly good. ;-)
>>
>> Do a quick google for "ray triangle intersection". There are several 
>> methods for doing this kind of thing, some of them highly optimised.
> 
> That tends to give results for a 3D problem, whereas my problem is 2D.

Many of the methods work for 2D as well though.

But now I'm curious... the cross product is only defined for 3D vectors, 
but you said you used a lot of cross products?

Post a reply to this message

> But now I'm curious... the cross product is only defined for 3D vectors, 
> but you said you used a lot of cross products?

Treat your 2D working area as a plane in 3D space (with Z coordinate zero if 
you like), then you can do cross products, and checking whether the result 
points "in to" the plane or "out of" the plane can be pretty useful (eg to 
find out which side of a line a point is on).

Post a reply to this message

scott wrote:
>> But now I'm curious... the cross product is only defined for 3D 
>> vectors, but you said you used a lot of cross products?
> 
> Treat your 2D working area as a plane in 3D space (with Z coordinate 
> zero if you like), then you can do cross products, and checking whether 
> the result points "in to" the plane or "out of" the plane can be pretty 
> useful (eg to find out which side of a line a point is on).

Oh, right.

...or you could just take the dot product instead, which requires only 2 
multiplies and one addition as opposed to four multiplies and two 
subtractions...

Post a reply to this message

> ...or you could just take the dot product instead, which requires only 2 
> multiplies and one addition as opposed to four multiplies and two 
> subtractions...

How is that useful?

Post a reply to this message

scott wrote:
>> ...or you could just take the dot product instead, which requires only 
>> 2 multiplies and one addition as opposed to four multiplies and two 
>> subtractions...
> 
> How is that useful?

Because you can use it to determine which side of the line a point is on?

Post a reply to this message

> Because you can use it to determine which side of the line a point is on?

How?

Example: Line through (0,0) and (0,1), and points either side of the line at 
(1,0) and (-1,0)

Post a reply to this message

scott wrote:
>> Because you can use it to determine which side of the line a point is on?
> 
> How?
> 
> Example: Line through (0,0) and (0,1), and points either side of the 
> line at (1,0) and (-1,0)

Suppose you have a line from point A to point B. Compute the vector AB, 
make it unital, and rotate through 90°; call the result V. Now compute 
the dot product V . A and call it k.

Which side of the line is point X on? Well, take the dot product V . X, 
and subtract k. The result is negative on one side, positive on the 
other, and zero if X is on the line itself (or at least, parallel to it).

...so basically, it's a ray/plane intersection test, except the "plane" 
is a 2D slide - a line.

A = (0,0)
B = (0,1)
AB = (0,1) - (0,0) = (0,1) [already unital]
V = (0,1) * {(0,-1), (1,0)} = (-1,0)
k = (0,0) . (-1,0) = 0

X = (1,0)
X . V - k = (1,0) . (-1,0) - 0 = -1

Y = (-1,0)
Y . V - k = (-1,0) . (-1,0) - 0 = +1

QED.

Post a reply to this message

> Suppose you have a line from point A to point B. Compute the vector AB, 
> make it unital, and rotate through 90°; call the result V. Now compute the 
> dot product V . A and call it k.
>
> Which side of the line is point X on? Well, take the dot product V . X, 
> and subtract k. The result is negative on one side, positive on the other, 
> and zero if X is on the line itself (or at least, parallel to it).
>
> ...so basically, it's a ray/plane intersection test, except the "plane" is 
> a 2D slide - a line.

Sorry, but that seems much more complicated and using up more instructions 
compared to the cross product:

> A = (0,0)
> B = (0,1)
> AB = (0,1) - (0,0) = (0,1) [already unital]

No need to make AB unital for the cross product method (saves a square root 
and divide).

> V = (0,1) * {(0,-1), (1,0)} = (-1,0)
> k = (0,0) . (-1,0) = 0

No need to calculate V or k for the cross product method (saves some 
multiplies and adds).

> X = (1,0)
> X . V - k = (1,0) . (-1,0) - 0 = -1

X cross AB = (1,0) x (0,1) = (1)(1) - (0)(0) = 1
(saves one subtraction compared to your method)

> Y = (-1,0)
> Y . V - k = (-1,0) . (-1,0) - 0 = +1

Y cross AB = (-1,0) x (0,1) = (-1)(1) - (0)(0) = -1
(saves one subtraction compared to your method)

Post a reply to this message

> No need to make AB unital for the cross product method (saves a square 
> root and divide).

Actually, thinking about it, you wouldn't need to for the dot product 
either. I should have seen that...

>> X = (1,0)
>> X . V - k = (1,0) . (-1,0) - 0 = -1
> 
> X cross AB = (1,0) x (0,1) = (1)(1) - (0)(0) = 1
> (saves one subtraction compared to your method)

Oh, I see. You're only computing the Z ordinate of the cross product? 
Well, I guess that could be quicker, yes...

Post a reply to this message

> Oh, I see. You're only computing the Z ordinate of the cross product? 

Well the other 2 ordinates are always zero, so there seems little point :-)

Post a reply to this message

scott wrote:
>> Oh, I see. You're only computing the Z ordinate of the cross product? 
> 
> Well the other 2 ordinates are always zero, so there seems little point :-)

...that's true, of course. I had it in my head that you would only 
compute the X and Y ordinates (which would take *more* effort than a 
simple dot product). But, of course, that would be pointless.

I should go sit in the corner now. :'{

Post a reply to this message

On 12/11/2009 02:43 AM, Slime wrote:
>> Have you tried this condition?
>>
>> ((distSq<  (radius1Sq+radius2Sq))
>> or (4*radius1Sq*radius2Sq>  (distSq - (radius1Sq+radius2Sq))^2))
>
> Smart. I feel stupid now! =) I see the equivalence, but I'm curious how long
> you took to come up with that. Was it obvious?
>
>   - Slime
>   [ http://www.slimeland.com/ ]
>
>
Well, it comes naturally from
distSq < radius1Sq + radius2Sq + 2 * radius1 * radius2
=>
distSq - (radius1Sq + radius2Sq ) <  2 * radius1 * radius2
=>
(distSq - (radius1Sq + radius2Sq ))^2 <  4 * radius1Sq * radius2Sq,

but the downside is that you have an overflow risk here to check

Post a reply to this message

bart <bart@home> wrote:
> Have you tried this condition?

> ((distSq < (radius1Sq+radius2Sq))
> or (4*radius1Sq*radius2Sq > (distSq - (radius1Sq+radius2Sq))^2))

  I was curious to see how much of a difference sqrt() really makes compared
to avoiding it, so I made this C++ file:

//-------------------------------------------------------------
#include <cmath>
#include <algorithm>

bool circlesCollide1(double cx1, double cy1, double rs1,
                     double cx2, double cy2, double rs2)
{
    const double vx = cx2 - cx1, vy = cy2 - cy1;
    return std::sqrt(vx*vx + vy*vy) < std::sqrt(rs1) + std::sqrt(rs2);
}

bool circlesCollide2(double cx1, double cy1, double rs1,
                     double cx2, double cy2, double rs2)
{
    const double vx = cx2 - cx1, vy = cy2 - cy1;
    const double dist2 = vx*vx + vy*vy;
    const double tmp = dist2 - (rs1 + rs2);
    return dist2 < rs1 + rs2 || 4.0*rs1*rs1 > tmp*tmp;
}
//-------------------------------------------------------------

and then a main program to test the speed:

//-------------------------------------------------------------
#include <iostream>
#include <ctime>

bool circlesCollide1(double cx1, double cy1, double rs1,
                     double cx2, double cy2, double rx2);
bool circlesCollide2(double cx1, double cy1, double rs1,
                     double cx2, double cy2, double rx2);

int main()
{
    double cx1 = -10.0, cy1 = -5.0, rs1_1 = 2.0*2.0, rs1_2 = 16.0*16.0;
    double cx2 = 5.0, cy2 = 10.0, rs2_1 = 2.5*2.5, rs2_2 = 18.0*18.0;
    const int iterations = 200000000;

    std::clock_t iClock = std::clock();
    for(int i = 0; i < iterations; ++i)
    {
        circlesCollide1(cx1, cy1, rs1_1, cx2, cy2, rs2_1);
        circlesCollide1(cx1, cy1, rs1_2, cx2, cy2, rs2_2);
    }
    std::clock_t eClock = std::clock();
    std::cout << "circlesCollide1: " << iterations << " iterations took "
              << double(eClock-iClock) / CLOCKS_PER_SEC << " seconds."
              << std::endl;

    iClock = std::clock();
    for(int i = 0; i < iterations; ++i)
    {
        circlesCollide2(cx1, cy1, rs1_1, cx2, cy2, rs2_1);
        circlesCollide2(cx1, cy1, rs1_2, cx2, cy2, rs2_2);
    }
    eClock = std::clock();
    std::cout << "circlesCollide2: " << iterations << " iterations took "
              << double(eClock-iClock) / CLOCKS_PER_SEC << " seconds."
              << std::endl;
}
//-------------------------------------------------------------

  I compiled with "-O3 -march=native" and got this result:

circlesCollide1: 200000000 iterations took 16.23 seconds.
circlesCollide2: 200000000 iterations took 5.42 seconds.

  (I didn't test that the functions give the correct result, though.)

-- 
                                                          - Warp

Post a reply to this message

On 12/11/2009 02:43 AM, Slime wrote:
>> Have you tried this condition?
>>
>> ((distSq<  (radius1Sq+radius2Sq))
>> or (4*radius1Sq*radius2Sq>  (distSq - (radius1Sq+radius2Sq))^2))
>
> Smart. I feel stupid now! =) I see the equivalence, but I'm curious how long
> you took to come up with that. Was it obvious?
>

hmm, it looks that this sqrt-free test
will equally work in 3D for sphere-sphere intersection test.

Post a reply to this message

> Suppose you have a line from point A to point B. Compute the vector AB, 
> make it unital, and rotate through 90°; call the result V. Now compute the 
> dot product V . A and call it k.

If you do this in 2D, minus the unnecessary normalization, the math is 
equivalent to doing the 3D cross product while ignoring the components you 
know to be 0.

inline vec_t Vec2Cross( const vec2_t v1, const vec2_t v2 )
{
 return v1[0] * v2[1] - v1[1] * v2[0];
}

If the result is positive, v2 points to the "left" of v1, otherwise it 
points to the "right".

 - Slime
 [ http://www.slimeland.com/ ]

Post a reply to this message

Warp wrote:

>   I was curious to see how much of a difference sqrt() really makes compared
> to avoiding it, so I made this C++ file:

>   I compiled with "-O3 -march=native" and got this result:
> 
> circlesCollide1: 200000000 iterations took 16.23 seconds.
> circlesCollide2: 200000000 iterations took 5.42 seconds.

I would imagine this sort of thing varies by processor, but yeah... I'm 
often curious to know how "expensive" various math operations are in 
relation to each other. (Of course, real programs are also affected by 
cache issues and which combinations of math ops you perform together and 
branch prediction issues and so forth... But typically you can't do much 
about that.)

Post a reply to this message

Invisible wrote:

> If you're that hell-bent on avoiding the square root, it can be 
> approximated as follows:
> 
>   Suppose that the most significant non-zero digit of x is in the 2^n 
> place. Let y = x / 2^(n/2). Now Sqrt(x) is approximately x / y.
> 
> Usually this estimated value is greater than the true square root (but 
> always by less than 45%). However, exact powers of 2 seem to come out 
> slightly below the true value. (By about 10^-16.) I don't know if that's 
> a glitch in my test program or what... It shouldn't be insurmountable 
> though.

One would have to be very hell-bent indeed to employ, in lieu of sqrt(), 
functions that are probably just as computationally expensive, like ln() 
and exp().  Unless you have a quick way of getting n and 2^(n/2) without 
ln() and exp(), you'll need them.

Now if you know the internal representation of your floats (or doubles), 
you can slap together a union like

union MajorKludge {
	double dvalue;
	char carray[sizeof double];
};

and get both n and 2^(n/2) using some very simple bitwise logical 
operations (as well as code whose portability will be highly suspect).

Regards,
John

Post a reply to this message

John VanSickle wrote:
> Invisible wrote:
> 
>> If you're that hell-bent on avoiding the square root, it can be 
>> approximated as follows:
>>
>>   Suppose that the most significant non-zero digit of x is in the 2^n 
>> place. Let y = x / 2^(n/2). Now Sqrt(x) is approximately x / y.
>>
>> Usually this estimated value is greater than the true square root (but 
>> always by less than 45%). However, exact powers of 2 seem to come out 
>> slightly below the true value. (By about 10^-16.) I don't know if 
>> that's a glitch in my test program or what... It shouldn't be 
>> insurmountable though.
> 
> One would have to be very hell-bent indeed to employ, in lieu of sqrt(), 
> functions that are probably just as computationally expensive, like ln() 
> and exp().  Unless you have a quick way of getting n and 2^(n/2) without 
> ln() and exp(), you'll need them.
> 
> Now if you know the internal representation of your floats (or doubles), 

Indeed, the *only* reason you would ever use a construction like this is 
because it's very fast and efficient to get at these numbers using the 
binary representation of a float.

-- 
http://blog.orphi.me.uk/
http://www.zazzle.com/MathematicalOrchid*

Post a reply to this message