phat code "The trick is to forget about writing something good and pull cheese out of your butt." - Necros
Main

Projects

Downloads

Articles

Links

Articles

Graphics

search for in   Tutorials / References
 

Circle-Triangle Intersection Method

AuthorThe_Grey_Beast (Gabriel Ivăncescu)
Emailgrey_predator@yahoo.com
WebsiteNone
ReleasedNov 1 2006
PlatformAll
Languagepseudo/C
Summary

This article presents 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.

 

Article Text

Printable Version / Download Article

Circle-Triangle Intersection Method

 

by The_Grey_Beast (Gabriel Ivăncescu)

 

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