Circle-Triangle Intersection Method
Foreword
In this article I present a method for checking whether a circle intersects a triangle (i.e true/false). It does not
need precalculated normals for each edge of the triangle, nor does it need to "normalize" the normals -- in other words,
it doesn't require square roots nor normal precalculation. The method is the result of my own thoughts -- thus I have
never read an article about circle-triangle intersection. However, as this algorithm may be common, it would be the case
that this will not be "original" at all -- but rest assured as I designed it myself without ever reading these kind of
algorithms, and that's a good thing since I can explain exactly what was in my mind and each steps which otherwise
could've been impossible to explain, since it wouldn't be me thinking about the method from scratch. :)
I am the kind of guy that thinks tutorials should be easy to follow, with examples and diagrams, not short and straight
to the point -- this is because, in my opinion, some algorithms are so blurred inside heavy mathematical formulas that I
cannot simply comprehend them at all. Some of the articles I read about other things were based on complex mathematical
terms without even describing the "bare bones" and how it got there. This article will focus on simplicity; the heavy
stuff (i.e stuff that I didn't figure out very quickly) will have more detailed explanations, and since I know how I
designed this, a step-by-step approach is taken in consideration.
I like tutorials to be understandable, enjoyable; not to make the readers just go into the code section and copy &
paste that code. It just bugs me -- the pseudo-code I use is very much C-ish, except for a few things, but is very much
natural and should be easy to grasp the idea. The code is not intended to be copied, rather to give you an example of how
it might look. I encourage everyone to write their own routine, specific for their environment or platform, and even apply
further specific-optimizations to it. All the low-level optimizations are best done in assembly anyway; don't bash
optimization or assembly language, as you will see the design phase is the one difficult; then you just have to focus on
the code for optimizations. However, choosing assembly at the first stage (i.e when designing) is not a good idea since
you will change design often, and you'll have to start all over again. It's best if we first decide on our design and,
when we choose a method and start to code it, then implement and optimize it. But be better than others and write it
yourself for your specific needs, that's just my advice :)
Anyway without further ado let's start, shall we?
Analyzing the Goals
Now what do we wish to accomplish? Oh yes, if a circle and a triangle intersect... let's analyze this "goal". We notice
the triangle uses 3 vectors for it's vertices, and the circle just a centre and a radius. Now, for a solid
circle (i.e not just it's outer-most ring) to intersect a triangle, the centre can be anywhere near or inside the triangle
so the triangle is within the radius.
As we know, any point on the same circle is of equal distance from the centre. By distance we mean the length of the
vector with origin located at the circle centre.
What do we mean with a triangle? 3 points connected through segments (edges). Of course, we use a solid triangle and a
solid circle (a solid circle is also called a "ball", but that's another story), so anything within these edges is within
the triangle, because it's solid (i.e not wireframe).
The first test that might arise in our heads is the one checking whether each of the triangle's vertices is within the
circle's radius -- this means that it intersects the circle, right? Take a look at the following fig:
At the pic above it is easy to observe that the triangle is intersecting the circle, because a vertex of the triangle
is inside the circle. This test is very simple, but of course doesn't handle the full intersection. The circle can be
totally inside the triangle, thus not intersecting vertices and neither edges.
Now, are there any other cases where the circle might intersect the triangle? Actually yes, the circle may intersect an
edge, no vertex, and it's centre could be outside the triangle. To simplify the vizualization, look at the following fig:
Ok, now that we observed three cases that can occur, it's natural to ask ourselves: is it necessary to test all these
cases? Of course not, as you have unlimited designs you can choose, some more efficient than others, but nevertheless you
don't need to do these three tests if you work out a completely different algorithm. However it is always a good idea to
analyze the situations that can occur before we decide on our algorithms because if we later discover that something else
can happen and our algorithm can't handle it, it might require a rewrite.
My method also uses three tests, but that is not the reason why I have shown the above "cases" which can occur. Now that
we analyzed the goals, it's time to move on to the intersection stages. My method presented below will use three tests to
check for the three cases discussed above. Again this is not the only way... actually I thought about another method (which
scales the triangle by an offset amount representing the radius, however it required that I either normalize or precalculate
normals which was not my goal).
Test 1: Vertex Within Circle
This test checks if a triangle's vertex is inside the circle or not -- it is actually simpler than you think. A
circle's radius represents a size (scalar). Now, if you create a vector straight up (parallel to the Y-axis) and make
it's size equal the radius, then rotating it around the centre will actually make a "circle" shape. Look at the following
fig:
In the above fig, the radius is initially vertical, fading as it is rotated in clockwise order (just to give you an
impression of continuation for the rotation). Sorry if the image doesn't make sense... It is easy to see that, as long as
a point relative to the circle's centre (which means it's origin is at the circle's centre) is smaller than the radius, it
falls inside the circle, right?
First let's define a few variables: c1, c2 and c3. These represent each vertex of the triangle, but
with the origin at the circle's centre. To do that we need to subtract the vectors. Simple enough:
c1 = v1 - centre
c2 = v2 - centre
c3 = v3 - centre
v1, v2 and v3 are the triangle's vertices, and centre is the circle's centre. Note however
that the centre and the vertices must have the same origin. Usually, when you do operation with vectors, they should be
adjusted to have the same origin before doing the calculation. It is also true in this case. Before attempting to subtract,
make sure that the triangle's vertices have the same origin. Manipulating the origin is easy: use the addition and
subtraction of vectors.. but that's not the scope of this article :P
Now that we computed the vertices relative to the circle's centre, we can easily see if their size is smaller than the
circle's radius. Here's a pic:
In the above pic the vertex we test is outside the circle because c1's length is larger than the radius
(r). I'm sure nobody had a problem seeing this, it's simple! We can calculate the length with Pythagorean
theorem (sqrt is the square root):
sqrt(x² + y²)
and then see if it is smaller than the radius.
if sqrt(c1x*c1x + c1y*c1y) <= radius
return true
c1x is the x component of c1, and c1y is the y component, of course. We should test
these for the other vertices as well. I know some of you might be thinking to evade that square root, but don't worry:
that's the part of optimizations -- but for now let's focus on the design. :)
Test 2: Circle Centre Within Triangle
Now that we know the circle doesn't intersect any of the triangle's vertices, we can go and test whether the circle's
centre is inside the triangle or not. This isn't actually very difficult -- it's enough to see if the centre, which
is a point, is on the positive side of all lines that make up the triangle. An edge is a segment of a line, but it doesn't
matter because we check if it is in front of ALL the three lines. This condition is true only when the point tested
(centre) is within the triangle. It doesn't matter that a line is infinite, because we test all of three of them and the
point has to be in front of ALL to actually be within the triangle. To make things clearer, look at the following diagram:
It is easy to see that all the line equations agree only within the triangle. Imagine that each normal points inward.
Then the test is positive if the point is on the side with the normal, negative otherwise. Now, examine all the three lines
and test all of them at once. You'll see that the result will only be positive within the triangle, since the point we test
will only be in front of each line at that area.
NOTE: The normals point "inward" because I assume we have our vertices in clockwise order. If you have your
vertices set up in counter-clockwise order, please reverse the comparison or operators specific to the sign of the result
(i.e instead of > 0, use < 0 because in counter-clockwise all the normals point outward, hence the test is negative for all
three lines only when we're inside the triangle this time; but you can also flip the equation without changing the
comparison at all). This line equation stuff is pretty basic, that's why I won't cover it in detail here. If you don't
know what a line equation is, imagine it's the same as a plane (in 3D), only that a line is a 2D plane. That's basically it.
Ok, now how do we implement all this stuff? Well, let's take a look at the line equation:
Ax + By + C = 0
A and B are the x and y components of the normal of the line, x and y are the
coordinates of any point on the line, and C is the distance from the origin to the line. This origin is the origin
of the triangle's vertices... just think of it as the absolute position (you know :)
What does this equation mean, exactly? It simply says that any point on the line makes up the line, which is more than
logical (i.e you don't create a line with points outside it, do you?). But in fact, the result (i.e 0 in the
equation) is the "signed" distance from the tested point to the line. Why have I said "signed"? Because if it's on the
positive side of the normal to the line, then the distance is positive, otherwise negative. That's also why the equation
says it's 0 because any point on the line is actually on the line, so it's distance to the line is 0, right?
Let's see an example:
N is the normal of the line. a is obviously further than c. b is on the line. Since a
is on the positive side of the normal, it has a positive distance. Since b is on the line, it's distance is 0.
c is on the negative side so it'll have a negative distance, but also a smaller distance than a. A simple
example for distances with arbitrary values would be: a:2, b:0, c:-0.5 I hope you get it, because
it's quite easy :)
Ok, now it's up to the math... Some of you might be screaming: "Hey, I thought your method does not involve normals?!?".
Yes, well in fact it doesn't need to have any normals precalculated, because it avoids the "normalization" of the normals
which is slow because of the square root involved. Well then how are we going to do this line equation if we don't have
normals, huh?
The trick is to ask ourselves: "Do we REALLY need the distance from the point to the line?". Answer: "Nope, we don't.
We only need the sign (positive/negative)."
But how do we get the correct sign without normals? Actually we DO calculate the normal, but we do NOT normalize it.
Calculating this "fake" normal is not very difficult. This is actually common stuff, so I won't go into details, but
anyway the normal of a line defined by 2 points (x1,y1) and (x2,y2) is:
A = (y2 - y1)
B = (x1 - x2)
Alternatively, just for the sake of it, B can be rewritten as the following (note that this does not change the
output at all):
B = -(x2 - x1)
An interesting sidenote is that you can flip the sign of these equations and you'll get the flipped normal. Thus, for
example if you have your vertices in counter-clockwise order, and you use the equations above you'll get all the
normals pointing outward the triangle. To solve this, you can flip the sign of the equations, i.e:
A = -(y2 - y1) = (y1 - y2)
B = -(x1 - x2) = (x2 - x1)
Now you'll get the normal pointing inward if you had your vertices in counter-clockwise order :)
As it can be seen from above, you can do all sorts of manipulations, some do not modify the output but may nevertheless
be useful for specific optimization purposes, other manipulations such as the last explained one can lead to a different
output. All is basic math.
Back to the point. We also have to calculate C, right? It turns out to not be necessary to do that, because I
chose a different trick. Let's examine the line equation. "Ax + By" actually means the dot product between
the normal and the point, doesn't it? Then we can also write the equation as:
N·p + C = 0
Where N is the normal and p is the point we test (circle centre in our case). But the dot product
is also:
N·p = |N|*|p|*cos(theta)
Where |N| is the length of N, |p| is the length of p, and θ is the angle
between N and p. To get rid of C we need to make sure N⋅p is 0 because the equation says the
result must be 0. When will N⋅p be 0? When either one of the lengths is 0 (which isn't happening in this case), or
when cos(θ) is 0. If the angle between N and p is 90 degrees, then it will be 0 since
cos(90°)=0. We have to dot the normal and a vector perpendicular to it. The normal is already perpendicular
to the line! This means the other "point" must be a vector parallel to the line. And this is the mighty segment between
two points on the line:
In the diagram above, V is the edge between v1 and v2. We can also say that V is the
vector v2 with the origin at v1. This "edge" vector is perpendicular to the normal, of course. So if we
dot these two vectors, we will arrive at the line equation. Thus we can rewrite the equation as follows:
N·(v2 - v1) = 0
Ok, so what? What good is this? It turns out that, by replacing v2 with the point we wish to test (circle
centre), then we'll get a positive result if the point is on the positive side of the normal, negative otherwise. This
is exactly the same "sign" as the line equation. And we don't even need to normalize the normal, because it's length is
always positive, so it won't affect the sign. And we're only interested in the sign. Thus for example:
N·(centre - v1) = sign
If we recall how the normal to a line between two points (x1,y1) and (x2,y2) is calculated, and using
the "two points" that define the line as the two vertices of the triangle that define the edge, of course, then we will
arrive at:
Nx edgex Ny edgey
(v2y - v1y)(centrex - v1x) + (v1x - v2x)(centrey - v1y) = sign
Or if you want to rearrange Ny and have the same output, then you can flip the + operator and the components
inside the parantheses of Ny too:
Nx edgex Ny edgey
(v2y - v1y)(centrex - v1x) - (v2x - v1x)(centrey - v1y) = sign
Mathematically, it is the same thing :)
Until here we only tested against the first edge (v2 - v1). Here's our first algo that we use to see if the
centre is on the positive side of ALL the edges:
if ((v2y - v1y)*(centrex - v1x) - (v2x - v1x)*(centrey - v1y)) >= 0 AND
((v3y - v2y)*(centrex - v2x) - (v3x - v2x)*(centrey - v2y)) >= 0 AND
((v1y - v3y)*(centrex - v3x) - (v1x - v3x)*(centrex - v3x)) >= 0
return true
And that's it! That's our "circle centre within triangle" test! Remember that you can play around with flipping
signs to get the opposite normal and/or result. You can also flip the comparison, to make it work if you need to. Also
remember that the arrangement of your vertices is important (i.e clockwise, counter-clockwise). It WILL affect the
result. If it doesn't work in a way, try flipping the comparison or operation (i.e multiply it by -1 and see what
result you get). But don't do both. Reverse the comparison OR flip the sign, but not both because you'll arrive at the
same result. It's all basic math :)
NOTE: This algo is made to work if your vertices are in clockwise order! If you have them in
counter-clockwise order, then multiply the left side by -1 or simply reverse the comparison
(i.e <= 0). Again you shouldn't do both of them. One fact why you should choose this is because you can take
advantage of specific optimizations from this manipulation. Here's a sample about multiplying the left side by -1
and rearranging the terms a bit to be more optimal (without the unary '-'):
if ((v2x - v1x)*(centrey - v1y) - (v2y - v1y)*(centrex - v1x)) >= 0 AND
((v3x - v2x)*(centrey - v2y) - (v3y - v2y)*(centrex - v2x)) >= 0 AND
((v1x - v3x)*(centrex - v3x) - (v1y - v3y)*(centrex - v3x)) >= 0
return true
Test 3: Circle Intersects Edge
We've come to the third test, and the worst of all :)
Now really, it isn't that hard.. It's just a matter of basic trigonometry. But first things first...
What we're trying to do is see if the circle intersects an edge of the triangle. We know at this point that the
circle's centre is outside. How are we going to do this? A logical situation would be to check "the shortest distance
from the centre to the edge", and see whether it is within the circle's radius. This is what we'll base our algorithm
on, but of course doesn't handle everything. If you don't know, the big secret is that the perpendicular from a point
to a line IS the shortest distance. So that means we just have to drop a perpendicular from the centre to the edge (line):
Uh huh, but how do we calculate this "perpendicular"? Actually, since we do not have normals, it is not as trivial as
you might think. But let's look at what we need: We do not need the perpendicular vector, but instead only it's
length. Ok, but how does this help? Let's just think a bit here...
Aha, what if we draw an imaginary triangle with points C (centre), v1 and O (perpendicular
point)? This imaginary triangle is also a right-triangle, because the perpendicular forms a 90 degree angle:
In this triangle we know the hypotenuse (opposite side of the 90 degree angle) which is p. What do I
mean about "we know"? Well, we know it because we can calculate it trivially since it's actually the length of the
vector between the centre and v1 (and even if this might require a square root at first glance, it can be
optimized). On the other hand, we do not know the point O, and calculating it is not at all efficient; but we do not
even have to, since we're not interested in the point location at all, only in the length.
If we sum this up, we realize we only know the hypotenuse, so it's not enough information to solve the
triangle (i.e we do not know d nor k).. but look at the angle formed by (centre, v1 and
O). This angle is also the angle between the edge vector (e1) and the vector between the centre
and v1 (c1).
What kind of operation can we use in this angle-situation? The first one that should probably ring a bell is the
dot product. Why? Recall that the dot product between c1 and e1 is:
c1·e1 = c1x*e1x + c1y*e1y
But it is also:
c1·e1 = |c1|*|e1|*cos(theta)
Where θ is the angle between them. Wait a minute. Did I say angle between them??? Yup, and it also
happens that it uses the cos function. Ok, let's think about what we can do with this info, hmm...
Well, we know that our triangle is a right-triangle, correct? So, the trigonometric functions work there, this
means we can find a use for that cos... Let's see, cos is Adjacent Side/Hypotenuse if we recall this
from basic trigonometric math, and this means for our angle k/p. Yup, that's what the cosine for that angle
is.. If we replace it in our dot product equation:
c1·e1 = |c1|*|e1|*(k/p)
As we have noted our variables, |c1| is p, so:
p*|e1|*k
c1·e1 = --------
p
Therefore we can simplify this equation without p. Now we need to find k. Why? We need d, not
k. Yes and no.. We'll see later that we DO need k. And besides, after we find out k, it's just a
matter of Pythagorean theorem in our triangle to solve d, since we know 2 sides and wish to find the third. But
for now let's find k:
c1·e1
k = -----
|e1|
Now that we have k we can calculate d easily with Pythagorean theorem. This says:
Side² + Side² = Hypotenuse²
And since we know a side (k) and the hypotenuse which is (p) we can replace 'em in the equation
and solve for the other side. Simple enough:
d² = p² - k²
Now all we have to do is check if d is within the radius... The above formula used the squared d
but we can get the root with a square root:
if sqrt(p*p - k*k) <= radius
return true
Ok, before you start screamin' at me, just remember that this implementation is not optimized, so that square root is
still there :P Don't worry, we'll optimize this greatly and achieve no square roots at all; I wouldn't want to break my
word :)
Alas we can't leave the design yet -- otherwise there would be bugs.. What?!? Well it turns out that we aren't done
yet.. This algo will actually report true more than we'd like. Why? It so happens that our algo handles the
intersection of the circle with a line, rather than an edge. An edge is not infinite, so we'll get
intersections far beyond the edge (on the same line).. Look at these pics to see the two sides where the circle can
intersect the line, but not the edge:
Perpendicular dropped at negative side of edge
Perpendicular dropped beyond the edge (positive side)
Handling each of these is done separately. Let's start with the negative side.
This check is not done on the perpendicular point as you might think, because we do not even calculate that point,
since it's pretty expensive. Well then let's focus on the other information available, namely the circle's centre.
Ok, what can we see about the centre when it is on the negative side? Well, first let's think about the vector c1
which is actually the vector from v1 to the centre.
Now for a bit of thinking and analyzing. Let's imagine two "circles", one that has it's perpendicular ON the edge
(which can intersect the triangle, if it is close enough), and the other circle having the perpendicular outside the
edge (negative side). Here's a pic to simplify this:
In the diagram above we see two circles, one which has it's perpendicular on the edge (the one at the bottom), and the
other one which is obviously not intersecting the edge, but rather the line (the circle at the top). The dotted
perpendicular line on the diagram represents the perpendicular to the edge, going through the point v1. The angles
between the two vectors going from v1 to each circle centres and the edge are also displayed. Go ahead and think a
bit about it :)
The answer is: "If the angle between the edge vector and the vector from v1 to the centre is more than
90 degrees, then the perpendicular from the centre will fall on the negative side of the edge". Just look at the diagram.
It is obvious that the top circle's centre (green vector) forms a big angle with the edge vector, bigger than 90° anyway.
If it's exactly 90° then the perpendicular from that point is exactly on the origin (which is v1). Otherwise given
a smaller angle, it falls on the positive side of the edge. If you don't understand, look again at the diagram and analyze
it, with my words in your head, to see that analyzing graphically is a very good way of discovering things :)
Ok, now that we know that the angle is the one we need to check, the question of how this can be accomplished arises.
How do we check the angle? We already did that with the dot product above. Let's see:
c1·e1 = |c1|*|e1|*cos(theta)
Where θ is the angle between them, hence the one we wish to test, right? Okay, but how can we test it
because it's inside that damn cos. Fortunately, cos is very helpful here :)
Remember the trigonometric unit circle (radius=1). cos is the X coordinate on that circle. The
angle is counted as counter-clockwise -- that is you go counter-clockwise on the circle and, the more you go, the bigger
the angle. I'm not going to explain this in detail since it's very basic trigonometric stuff, and it's probably explained
much better elsewhere. And I assume you know this. Also remember we're only interested in the first 180° because the
other half of the circle is exactly symmetric for the cos function and even the dot product has no clue :)
Let's see what cos is for angles bigger than 90°. Well, on the unit circle it's obvious that the
X coord for those angles (remember we're NOT interested in angles > 180°) is negative, hence we only have to see
if the result of the dot product is negative. And since we already calculated the dot product because we had to
find k (explained above), then we can simply put a condition after we calculate the dot product to see
whether it is negative or not. And that's it for this case:
if c1·e1 > 0
Then it falls on the positive side of the edge.. Phew :)
NOTE: If you're wondering why we can simply see if the dot product is negative since there are other
variables involved besides cos (i.e |c1| and |e1|) then rest assured because |c1| and
|e1| are only lengths and they are ALWAYS positive, so they will NOT affect the sign. Only cos
decides the sign of the result.
However we also have to test if the perpendicular is dropped beyond the edge. This is very easy. Recall that we
calculated k which is the length of the segment formed by the perpendicular and v1. If you're not sure,
look at our "imaginary" triangle above (at the beginning of this section).
The only test we have to do is to check if k is larger than the edge's length, correct? Because if k
is bigger, then the perpendicular certainly falls beyond the edge.. Simple enough:
if k < len
After these two tests can we finally go ahead with the "perpendicular length is within radius" test :)
NOTE: Some of you might be wondering: "What happens at the corners?". Our edge-intersection algo doesn't
handle the corners well. That's because the "shortest" distance is not necessarily the "only" distance from the centre
to the edge. The corners are problematic because the edge is not infinite, so the perpendicular may as well be outside
while the circle might still intersect the edge. Look at this diagram to vizualize it easier:
Fortunately we are already done. It so happens that our first test (vertex within circle) handles this. At the
corners, the circle intersects the vertices, so we don't have to worry one bit -- it won't even reach our
edge-intersection algo because it will immediately see that the circle is intersecting the vertex. Phew :) Finally
we're done with the design. We should seriously start to optimize this, but anyway here's the code so far...
First Working Implementation
;
; TEST 1: Vertex within circle
;
c1x = v1x - centrex
c1y = v1y - centrey
if sqrt(c1x*c1x + c1y*c1y) <= radius
return true
c2x = v2x - centrex
c2y = v2y - centrey
if sqrt(c2x*c2x + c2y*c2y) <= radius
return true
c3x = v3x - centrex
c3y = v3y - centrey
if sqrt(c3x*c3x + c3y*c3y) <= radius
return true
;
; TEST 2: Circle centre within triangle
;
; NOTE: This works for clockwise ordered vertices!
;
if ((v2y - v1y)*(centrex - v1x) - (v2x - v1x)*(centrey - v1y)) >= 0 AND
((v3y - v2y)*(centrex - v2x) - (v3x - v2x)*(centrey - v2y)) >= 0 AND
((v1y - v3y)*(centrex - v3x) - (v1x - v3x)*(centrex - v3x)) >= 0
return true
;
; TEST 3: Circle intersects edge
;
; Get the dot product...
;
c1x = centrex - v1x
c1y = centrey - v1y
e1x = v2x - v1x
e1y = v2y - v1y
k = c1x*e1x + c1y*e1y
if k > 0
{
len = sqrt(e1x*e1x + e1y*e1y)
k = k/len
if k < len
{
if sqrt(c1x*c1x + c1y*c1y - k*k) <= radius
return true
}
}
; Second edge
c2x = centrex - v2x
c2y = centrey - v2y
e2x = v3x - v2x
e2y = v3y - v2y
k = c2x*e2x + c2y*e2y
if k > 0
{
len = sqrt(e2x*e2x + e2y*e2y)
k = k/len
if k < len
{
if sqrt(c2x*c2x + c2y*c2y - k*k) <= radius
return true
}
}
; Third edge
c3x = centrex - v3x
c3y = centrey - v3y
e3x = v1x - v3x
e3y = v1y - v3y
k = c3x*e3x + c3y*e3y
if k > 0
{
len = sqrt(e3x*e3x + e3y*e3y)
k = k/len
if k < len
{
if sqrt(c3x*c3x + c3y*c3y - k*k) <= radius
return true
}
}
; We're done, no intersection
return false
Optimizations
This piece of code is way too slow. Especially the edge test -- I believe it is even slower than calculating the
normals on the fly. However this can be optimized greatly. This is what we'll now focus on: code. We'll not think
anymore about math, etc... let's analyze the code.
The straightforward optimizations that do not require low-level attention should be attended first, at least in my
opinion. We should attempt to reuse as many variables as possible, and never calculate anything twice, mainly if it
requires a lot of operations or expensive operations. In C you can use "references" to other local variables such that
you reuse them. You can also write it in assembly -- there you will have the greatest control over your variables. But
it is a good idea to reuse them as much as possible.
We can see that c1, c2 and c3 are required both at tests 1 and 3. In fact, they are also
required for the second test! Just look at the second test, where you will find expressions of the form
(centrex - v1x). That's exactly the c1 we need. However, if you take a look at the first test we'll see
that c1 is defined as (v1x - centrex) which isn't the same thing...
But in fact what do we need at the first test? We need c1's length, not it's sense. We do not care if it is
facing away from the centre or to the centre, the length is the same. So, we can rewrite the first test as:
c1x = centrex - v1x
c1y = centrey - v1y
if sqrt(c1x*c1x + c1y*c1y) <= radius
return true
c2x = centrex - v2x
c2y = centrey - v2y
if sqrt(c2x*c2x + c2y*c2y) <= radius
return true
c3x = centrex - v3x
c3y = centrey - v3y
if sqrt(c3x*c3x + c3y*c3y) <= radius
return true
Now that we precalculated c1, c2 and c3, we can go and replace 'em in the second test
accordingly and also remove their calculation from the third test. This is a good demonstration of "don't calculate
anything twice" optimization. And also a bit of "reuse as many variables as possible". Ok let's see the second
test in action:
if ((v2y - v1y)*c1x - (v2x - v1x)*c1y) >= 0 AND
((v3y - v2y)*c2x - (v3x - v2x)*c2y) >= 0 AND
((v1y - v3y)*c3x - (v1x - v3x)*c3y) >= 0
return true
Ok, it starts to be more compact.. However we can also precalculate (v2y - v1y) and the like.. why? That's
exactly the edges we calculate in the third test. By calculating 'em in the second test, we can use them in the second
test and also in the third test. Again, try to not calculate anything twice :)
Still don't get it? Look at the third test where we define our edges (e1, e2 and e3). And now
look at the second test. There are the same expressions found there, so why calculate them twice? Let's put this in the
second test. We just have to remove their declaration from the third test. Here's the layout of the code:
c1x = centrex - v1x
c1y = centrey - v1y
if sqrt(c1x*c1x + c1y*c1y) <= radius
return true
c2x = centrex - v2x
c2y = centrey - v2y
if sqrt(c2x*c2x + c2y*c2y) <= radius
return true
c3x = centrex - v3x
c3y = centrey - v3y
if sqrt(c3x*c3x + c3y*c3y) <= radius
return true
e1x = v2x - v1x
e1y = v2y - v1y
e2x = v3x - v2x
e2y = v3y - v2y
e3x = v1x - v3x
e3y = v1y - v3y
if e1y*c1x - e1x*c1y >= 0 AND
e2y*c2x - e2x*c2y >= 0 AND
e3y*c3x - e3x*c3y >= 0
return true
k = c1x*e1x + c1y*e1y
if k > 0
{
len = sqrt(e1x*e1x + e1y*e1y)
k = k/len
if k < len
{
if sqrt(c1x*c1x + c1y*c1y - k*k) <= radius
return true
}
}
k = c2x*e2x + c2y*e2y
if k > 0
{
len = sqrt(e2x*e2x + e2y*e2y)
k = k/len
if k < len
{
if sqrt(c2x*c2x + c2y*c2y - k*k) <= radius
return true
}
}
k = c3x*e3x + c3y*e3y
if k > 0
{
len = sqrt(e3x*e3x + e3y*e3y)
k = k/len
if k < len
{
if sqrt(c3x*c3x + c3y*c3y - k*k) <= radius
return true
}
}
return false
Ok, these straightforward optimizations won't make it much faster, but it's always a good idea to apply optimizations
wherever possible. Now that we organized this up a bit, it's time to focus on the square roots. Don't worry we'll get
back to these kind of optimizations but first let's finish off the square roots :)
By far the easiest approach to whack the square roots from the first test is to just test everything against the
squared result. For example, in math:
x <= r
Can be rewritten as:
x² <= r²
Which is the same thing, right? So it turns out we just have to precalculate the squared radius and use it whenever
the need arises. Thus in the first test:
c1x = centrex - v1x
c1y = centrey - v1y
; Precalculate squared radius
radiusSqr = radius*radius
if c1x*c1x + c1y*c1y <= radiusSqr
return true
c2x = centrex - v2x
c2y = centrey - v2y
if c2x*c2x + c2y*c2y <= radiusSqr
return true
c3x = centrex - v3x
c3y = centrey - v3y
if c3x*c3x + c3y*c3y <= radiusSqr
return true
A similar approach is done for the third test too, where the comparison against the radius is done the same way. Also,
as a tip, if you have a constant-size circle and don't need the "standard" radius at all (this algo doesn't need it),
then you can replace it with the squared radius and precalculate it ONLY once for the entire application. Instead of
keeping the radius you keep the squared radius. Or you can keep them both if you really need the standard radius.
Anyway I'll best leave this to you :)
The second test has no square root so let's pass it. The third test, however, is bloated with square roots. Given
our previous method of eliminating the square roots, here's how the third test looks for the first edge:
k = c1x*e1x + c1y*e1y
if k > 0
{
len = sqrt(e1x*e1x + e1y*e1y)
k = k/len
if k < len
{
if c1x*c1x + c1y*c1y - k*k <= radiusSqr
return true
}
}
Well we still have a square root lying somewhere... but what if we do the same trick? I mean, instead of calculating
k we'll calculate the squared k. I chose the same name instead of kSqr, but you can do however you
like. Instead of calculating len we'll calculate the squared len. When we divide k by this squared
len, we'll not do it mathematically correct. Look at this:
k
---
len
Is not the same as:
k
----
len²
If we square everything, it should be:
k²
----
len²
So it means we have to use k = k*k / len to get the squared k. Simple enough:
k = c1x*e1x + c1y*e1y
if k > 0
{
len = e1x*e1x + e1y*e1y ; squared len
k = k*k/len ; squared k
if k < len ; squared test; same
{
if c1x*c1x + c1y*c1y - k <= radiusSqr
return true
}
}
In case you're wondering why, at the last condition, we used k instead of k*k, rest assured as
k is already squared. In my opinion, this method fits perfectly: we calculate the squared k and squared
len to get rid of the square roots, and the Pythagorean theorem also requires the squared k. Like I said,
it fits perfectly :)
Wahoo, we just got rid of the square roots. Quite an impressive achievement. We have absolutely no square root at all,
nor need any precalculation of normals.. in short, it's quite fast. :)
There are also other general-purpose optimizations possible. Told ya' we'll get back to them. Let's see the first test
where we calculate c1's length. We check if this length (actually squared length) is within the (squared)
radius. Ok, it's fine. But now look at test 3. We calculate c1's length there too to see if the (squared) shortest
distance is within the (squared) radius. So we calculate it twice. Let's precalculate it. The names used are c1sqr,
c2sqr and c3sqr. Ok now here's the code so far:
;
; TEST 1: Vertex within circle
;
c1x = centrex - v1x
c1y = centrey - v1y
radiusSqr = radius*radius
c1sqr = c1x*c1x + c1y*c1y
if c1sqr <= radiusSqr
return true
c2x = centrex - v2x
c2y = centrey - v2y
c2sqr = c2x*c2x + c2y*c2y
if c2sqr <= radiusSqr
return true
c3x = centrex - v3x
c3y = centrey - v3y
c3sqr = c3x*c3x + c3y*c3y
if c3sqr <= radiusSqr
return true
;
; TEST 2: Circle centre within triangle
;
; NOTE: This works for clockwise ordered vertices!
;
;
; Calculate edges
;
e1x = v2x - v1x
e1y = v2y - v1y
e2x = v3x - v2x
e2y = v3y - v2y
e3x = v1x - v3x
e3y = v1y - v3y
if e1y*c1x - e1x*c1y >= 0 AND
e2y*c2x - e2x*c2y >= 0 AND
e3y*c3x - e3x*c3y >= 0
return true
;
; TEST 3: Circle intersects edge
;
; Get the dot product...
;
k = c1x*e1x + c1y*e1y
if k > 0
{
len = e1x*e1x + e1y*e1y ; squared len
k = k*k/len ; squared k
if k < len ; squared test
{
if c1sqr - k <= radiusSqr
return true
}
}
; Second edge
k = c2x*e2x + c2y*e2y
if k > 0
{
len = e2x*e2x + e2y*e2y
k = k*k/len
if k < len
{
if c2sqr - k <= radiusSqr
return true
}
}
; Third edge
k = c3x*e3x + c3y*e3y
if k > 0
{
len = e3x*e3x + e3y*e3y
k = k*k/len
if k < len
{
if c3sqr - k <= radiusSqr
return true
}
}
; We're done, no intersection
return false
Further specific-wise optimizations should be handled now, as they are a bit more middle-level. This code is really
just pseudo-code. It was intended to be understandable, not copied. I encourage you to write your own routine, and even
apply specific optimizations for your needs; i.e my code was general-purpose. Now if you are blurred and don't know where
to start from, then I'll give some short examples here :)
For instance, we use a comparison with c1 against the squared radius two times -- one for the first
test and one for the third test. We could optimize this for specific purpose; that is we could modify c1sqr such
that the comparison is tested against 0. In math:
x <= y
We can move y around in the equation, flipping it's sign when we move it to the other side. Just like a normal
equation:
x - y <= 0
Why 0? There's nothing on the right after we move y on the left side, that's why.. this is very basic math,
but just to remind you of it...
Now we can subtract from c1sqr the squared radius, so that we can test c1sqr against 0. How does this
help? First let's see what happens:
c1x = centrex - v1x
c1y = centrey - v1y
radiusSqr = radius*radius
c1sqr = c1x*c1x + c1y*c1y - radiusSqr
if c1sqr <= 0
return true
Ok this code doesn't look more optimal, but recall the third test that uses this condition:
if c1sqr - k <= radiusSqr
return true
With our previous modification, c1sqr already contains the subtracted radiusSqr, so the comparison
becomes:
if c1sqr - k <= 0
return true
Which, by moving k to the right, becomes:
if c1sqr <= k
return true
Which requires only a comparison instead of a subtraction and a comparison. Now we do not even need to store the
squared radius after we compute c1sqr, c2sqr and c3sqr, so we can reuse the variable for
storing other information, for example c3sqr. Look at this sample:
radiusSqr = radius*radius
c1sqr = c1x*c1x + c1y*c1y - radiusSqr
c2sqr = c2x*c2x + c2y*c2y - radiusSqr
c3sqr = c3x*c3x + c3y*c3y - radiusSqr ; at this point we do not need radiusSqr anymore
So we can reuse radiusSqr (the variable in memory, I mean...). With C this can be done with
references (i.e references are actually nothing -- they just provide an alternative access to a variable, or memory
location, through a different name). Why do we use references instead of plain and simple radiusSqr? Because
c3sqr is a much more convenient name for this purpose :)
radiusSqr = radius*radius ; here is the variable declaration
c1sqr = c1x*c1x + c1y*c1y - radiusSqr ; another variable
c2sqr = c2x*c2x + c2y*c2y - radiusSqr ; another variable
&c3sqr = radiusSqr ; declare reference to radiusSqr
c3sqr = c3x*c3x + c3y*c3y - radiusSqr
The last assignment will be stored in radiusSqr's location in memory. Thus it is the same as if you wrote:
radiusSqr = c3x*c3x + c3y*c3y - radiusSqr
References are very powerful optimizing tools that you can use. Even more powerful as this example because you can use
different types (i.e access the lower 2 bytes of unsigned long variable a through a unsigned short
reference b). Unlike what some people say, references do NOT eat any memory or do not generate extra code at ALL.
They are nothing -- just name aliases. It's same as if you used #define directives, though some compilers will
optimize better when you use references :)
In assembly this whole reference thing is very obvious and simple, since I suppose you know a thing or two about
low-level programming (asm). You can write this algo in assembly and have even more control over your code.
We can extend this further with a sneaky trick to get rid of the division in the third test. How? Well, the third
test looks like this (note that this takes in account the previous optimization):
k = k*k/len
if k < len
{
if c1sqr <= k
return true
}
What happens if we don't divide? Look at how the first assignment looks mathematically:
k²
--- < len
len
We can take out the division if we multiply the right side with the denominator:
k² < len * len
Thus now we have:
k = k*k
if k < len*len
{
if c1sqr <= k/len
return true
}
In case you're wondering why there is a k/len in the second comparison, well it's because we didn't
divide k by len as we were supposed to, and we didn't fix the second comparison yet. How about we fix it,
eh? Ok, we can use the same thing really:
k
c1sqr <= ---
len
Which, by using the same method, gives us the following code:
k = k*k
if k < len*len
{
if c1sqr*len <= k
return true
}
We can also optimize this even more. Look at the first assignment: it calculates the squared k. And then look
at the condition that follows: it checks if the squared k is less than the squared len. So:
k² < len²
But, by doing the inverse of our trick that we used in eliminating the square roots, we can write this as:
k < len
Same thing, right? Instead of squaring everything, we took the square root of everything, which is the same thing
for a comparison, right? Ok, so we have this piece of code:
if k < len ; non-squared k checked against len
{
if c1sqr*len <= k*k ; we need to calculate squared k
return true
}
Note however that you can have more overflow problems than with the divide method, because you work with much larger
numbers. This is especially an issue if you use fixed point formats. Be sure to keep overflowing under control. This
overflow handling, and more assembly-level or specific optimizations are left to you as an exercise.
Now I will briefly discuss optimizations that you can do to the second test too. If we look at it, it's easy to
realize that we can use the same "comparison" trick here too. Let's take a look at it:
if e1y*c1x - e1x*c1y >= 0 AND
e2y*c2x - e2x*c2y >= 0 AND
e3y*c3x - e3x*c3y >= 0
return true
We can move the second statements to the right side, thus having only a comparison instead of (on some platforms) a
subtraction and a comparison:
if e1y*c1x >= e1x*c1y AND
e2y*c2x >= e2x*c2y AND
e3y*c3x >= e3x*c3y
return true
Or if you want to evade three conditions, you can optimize this better. However I'll give you only one example as
this tends to be specific-purpose.
Ok, assume we have our numbers in a simple format that supports bitwise logical operators (and, or, xor, not),
preferably fixed point format. And assume that the most significant bit (i.e left-most bit in the number) is
1 when the number is negative, and 0 when positive (examples include two's complement,
one's complement, signed magnitude and other representations). Ok, my example is actually based on x86
architecture with fixed point types (integers). Now, what can we do in this situation? Simple.
We can use a bitwise operation to indicate that ALL the numbers had the most-significant bit 0 (positive). What
kind of bitwise operation? Answer: OR (written as '|'). Why? Because if at least one of the line equations
is negative (1 in the most significant bit), then the result MUST be negative, so we can be sure it certainly
isn't inside the triangle. If ALL three of them were positive (0 in most significant bit) then the result's most
significant bit should also be 0 and in this case we know it is inside the triangle:
k = (e1y*c1x - e1x*c1y) | (e2y*c2x - e2x*c2y) | (e3y*c3x - e3x*c3y)
if (signed)k > 0
return true
Why have I used the term "signed"? Well because, if you used that representation, the signed value of k
should be positive if the most significant bit is 0, thus >= 0. If the most significant bit is 1, then the
number is < 0. If you're still confused, you can also use:
if (k & 0x80000000) == 0
return true
This is assuming k is 32-bits wide. 0x80000000 means a hexadecimal value with a 1 only at the
most significant bit, and & is the bitwise AND operator. You can do all sorts of these low-level
optimizations. But anyway much of this stuff is better done in assembly. Be creative and inspire yourself with sneaky
algos :)
I believe I have covered enough optimization examples. As a brief overall view of our optimizations:
;
; TEST 1: Vertex within circle
;
c1x = centrex - v1x
c1y = centrey - v1y
radiusSqr = radius*radius
c1sqr = c1x*c1x + c1y*c1y - radiusSqr
if c1sqr <= 0
return true
c2x = centrex - v2x
c2y = centrey - v2y
c2sqr = c2x*c2x + c2y*c2y - radiusSqr
if c2sqr <= 0
return true
c3x = centrex - v3x
c3y = centrey - v3y
&c3sqr = radiusSqr ; reference to radiusSqr
c3sqr = c3x*c3x + c3y*c3y - radiusSqr
if c3sqr <= 0
return true
;
; TEST 2: Circle centre within triangle
;
;
; Calculate edges
;
e1x = v2x - v1x
e1y = v2y - v1y
e2x = v3x - v2x
e2y = v3y - v2y
e3x = v1x - v3x
e3y = v1y - v3y
if signed((e1y*c1x - e1x*c1y) | (e2y*c2x - e2x*c2y) | (e3y*c3x - e3x*c3y)) >= 0
return true
;
; TEST 3: Circle intersects edge
;
k = c1x*e1x + c1y*e1y
if k > 0
{
len = e1x*e1x + e1y*e1y ; squared len
if k < len
{
if c1sqr * len <= k*k
return true
}
}
; Second edge
k = c2x*e2x + c2y*e2y
if k > 0
{
len = e2x*e2x + e2y*e2y
if k < len
{
if c2sqr * len <= k*k
return true
}
}
; Third edge
k = c3x*e3x + c3y*e3y
if k > 0
{
len = e3x*e3x + e3y*e3y
if k < len
{
if c3sqr * len <= k*k
return true
}
}
; We're done, no intersection
return false
This was only an example to give you an overall impression of our optimizations. Note that I included the "specific"
optimizations which may not work for your specific needs -- this was only an example. Go ahead and adjust it to work for
you, or optimize it specifically. :)
Well, what's to say? We have completed a method of checking whether a circle intersects a triangle, analyzed it, and
optimized it... the rest is up to you to do what you want with it.
Summary
This article focused on a step-by-step approach -- that is, how you should (in my opinion) create a method. First you
have to design it, and think in high-level terms. This is because you don't need to be confused at this important stage.
After you have a workable version, go and optimize it like hell. And remember, the design is the difficult part, not
writing the code, nor the low-level optimizations -- you should really take the time to optimize or implement it in
assembly. Even if it isn't a really noticeable result, it's always worth optimizing. You'll see that assembly or other
low-level languages do not take very much extra time. It's the design that's difficult.
If you find any errors, or want to comment about articles, methods and optimizations, or simply need to contact me,
drop me a mail at
grey_predator@yahoo.com.
However as I'm not very frequent on the net, I cannot guarantee immediate
response, so please do not get frustrated if I do not reply very soon. I will do my best.
I can also be found on forums. My nicks are usually The_Grey_Beast or Bidirectional. However some people
just make idiotic jokes and use other people's nicks and then claim they're the respective guy. This is not a very
mature thing and others should judge a guy by his personality, not nick. However I usually use those two nicks.
What else to say? I apologize if I forgot something, it's so hard to keep it tight and yet explain the details.
Please understand that you are not allowed to modify this document without my permission. The reason is to keep
things under control (i.e not to have 500 versions of this doc around the net). I hope you understand. If you need to
modify or do other things with it, please contact me first.
Best wishes and happy coding :) |