algorithm - Compute the area of intersection between a circle and a triangle?

Question

How does one compute the area of intersection between a triangle (specified as three (X,Y) pairs) and a circle (X,Y,R)? I've done some searching to no avail. This is for work, not school. :)

It would look something like this in C#:

struct { PointF vert[3]; } Triangle;
struct { PointF center; float radius; } Circle;

// returns the area of intersection, e.g.:
// if the circle contains the triangle, return area of triangle
// if the triangle contains the circle, return area of circle
// if partial intersection, figure that out
// if no intersection, return 0
double AreaOfIntersection(Triangle t, Circle c)
{
 ...
}

Answer 1

First I will remind us how to find the area of a polygon. Once we have done this, the algorithm to find the intersection between a polygon and a circle should be easy to understand.

How to find the area of a polygon

Let's look at the case of a triangle, because all the essential logic appears there. Let's assume we have a triangle with vertices (x1,y1), (x2,y2), and (x3,y3) as you go around the triangle counter-clockwise, as shown in the following figure:

Then you can compute the area by the formula

A=(x1 y2 + x2 y3 + x3 y1 - x2y1- x3 y2 - x1y3)/2.

To see why this formula works, let's rearrange it so it is in the form

A=(x1 y2 - x2 y1)/2 + (x2 y3 - x3 y2)/2 + (x3 y1 - x1y3 )/2.

Now the first term is the following area, which is positive in our case:

If it isn't clear that the area of the green region is indeed (x1 y2 - x2 y1)/2, then read this.

The second term is this area, which is again positive:

And the third area is shown in the following figure. This time the area is negative

Adding these three up we get the following picture

We see that the green area that was outside the triangle is cancelled by the red area, so that the net area is just the area of the triangle, and this shows why our formula was true in this case.

What I said above was the intuitive explanation as to why the area formula was correct. A more rigorous explanation would be to observe that when we calculate the area from an edge, the area we get is the same area we would get from integration r^2dθ/2, so we are effectively integrating r^2dθ/2 around the boundary of the polygon, and by stokes theorem, this gives the same result as integrating rdrdθ over the region bounded the polygon. Since integrating rdrdθ over the region bounded by the polygon gives the area, we conclude that our procedure must correctly give the area.

Area of the intersection of a circle with a polygon

Now let's discuss how to find the area of the intersection of a circle of radius R with a polygon as show in the following figure:

We are interested in find the area of the green region. We may, just as in the case of the single polygon, break our calculation into finding an area for each side of the polygon, and then add those areas up.

Our first area will look like:

The second area will look like

And the third area will be

Again, the first two areas are positive in our case while the third one will be negative. Hopefully the cancellations will work out so that the net area is indeed the area we are interested in. Let's see.

Indeed the sum of the areas will be area we are interested in.

Again, we can give a more rigorous explanation of why this works. Let I be the region defined by the intersection and let P be the polygon. Then from the previous discussion, we know that we want to computer the integral of r^2dθ/2 around the boundary of I. However, this difficult to do because it requires finding the intersection.

Instead we did an integral over the polygon. We integrated max(r,R)^2 dθ/2 over the boundary of the polygon. To see why this gives the right answer, let's define a function π which takes a point in polar coordinates (r,θ) to the point (max(r,R),θ). It shouldn't be confusing to refer to the the coordinate functions of π(r)=max(r,R) and π(θ)=θ. Then what we did was to integrate π(r)^2 dθ/2 over the boundary of the polygon.

On the other hand since π(θ)=θ, this is the same as integrating π(r)^2 dπ(θ)/2 over the boundary of the polygon.

Now doing a change of variable, we find that we would get the same answer if we integrated r^2 dθ/2 over the boundary of π(P), where π(P) is the image of P under π.

Using Stokes theorem again we know that integrating r^2 dθ/2 over the boundary of π(P) gives us the area of π(P). In other words it gives the same answer as integrating dxdy over π(P).

Using a change of variable again, we know that integrating dxdy over π(P) is the same as integrating Jdxdy over P, where J is the jacobian of π.

Now we can split the integral of Jdxdy into two regions: the part in the circle and the part outside the circle. Now π leaves points in the circle alone so J=1 there, so the contribution from this part of P is the area of the part of P that lies in the circle i.e., the area of the intersection. The second region is the region outside the circle. There J=0 since π collapses this part down to the boundary of the circle.

Thus what we compute is indeed the area of the intersection.

Now that we are relatively sure we know conceptually how to find the area, let's talk more specifically about how to compute the contribution from a single segment. Let's start by looking at a segment in what I will call "standard geometry". It is shown below.

In standard geometry, the edge goes horizontally from left to right. It is described by three numbers: xi, the x-coordinate where the edge starts, xf, the x-coordinate where the edge ends, and y, the y coordinate of the edge.

Now we see that if |y| < R, as in the figure, then the edge will intersect the circle at the points (-xint,y) and (xint,y) where xint = (R^2-y^2)^(1/2). Then the area we need to calculate is broken up into three pieces labelled in the figure. To get the areas of regions 1 and 3, we can use arctan to get the angles of the various points and then equate the area to R^2 Δθ/2. So for example we would set θi = atan2(y,xi) and θl = atan2(y,-xint). Then the area of region one is R^2 (θl-θi)/2. We can obtain the area of region 3 similarly.

The area of region 2 is just the area of a triangle. However, we must be careful about sign. We want the area shown to be positive so we will say the area is -(xint - (-xint))y/2.

Another thing to keep in mind is that in general, xi does not have to be less than -xint and xf does not have to be greater than xint.

The other case to consider is |y| > R. This case is simpler, because there is only one piece which is similar to region 1 in the figure.

Now that we know how to compute the area from an edge in standard geometry, the only thing left to do is describe how to transform any edge into standard geometry.

But this just a simple change of coordinates. Given some with initial vertex vi and final vertex vf, the new x unit vector will be the unit vector pointing from vi to vf. Then xi is just the displacement of vi from the center of the circle dotted into x, and xf is just xi plus the distance between vi and vf. Meanwhile y is given by the wedge product of x with the displacement of vi from the center of the circle.

Code

That completes the description of the algorithm, now it is time to write some code. I will use java.

First off, since we are working with circles, we should have a circle class

public class Circle {

    final Point2D center;
    final double radius;

    public Circle(double x, double y, double radius) {
        center = new Point2D.Double(x, y);
        this.radius = radius;
    }

    public Circle(Point2D.Double center, double radius) {
        this(center.getX(), center.getY(), radius);
    }

    public Point2D getCenter() {
        return new Point2D.Double(getCenterX(), getCenterY());
    }

    public double getCenterX() {
        return center.getX();
    }

    public double getCenterY() {
        return center.getY();
    }

    public double getRadius() {
        return radius;
    }

}

For polygons, I will use java's Shape class. Shapes have a PathIterator that I can use to iterate through the edges of the polygon.

Now for the actual work. I will separate the logic of iterating through the edges, putting the edges in standard geometry etc, from the logic of computing the area once this is done. The reason for this is that you may in the future want to compute something else besides or in addition to the area and you want to be able to reuse the code having to deal with iterating through the edges.

So I have a generic class which computes some property of class T about our polygon circle intersection.

public abstract class CircleShapeIntersectionFinder<T> {

It has three static methods that just help compute geometry:

private static double[] displacment2D(final double[] initialPoint, final double[] finalPoint) {
    return new double[]{finalPoint[0] - initialPoint[0], finalPoint[1] - initialPoint[1]};
}

private static double wedgeProduct2D(final double[] firstFactor, final double[] secondFactor) {
    return firstFactor[0] * secondFactor[1] - firstFactor[1] * secondFactor[0];
}

static private double dotProduct2D(final double[] firstFactor, final double[] secondFactor) {
    return firstFactor[0] * secondFactor[0] + firstFactor[1] * secondFactor[1];
}

There are two instance fields, a Circle which just keeps a copy of the circle, and the currentSquareRadius, which keeps a copy of the square radius. This may seem odd, but the class I am using is actually equipped to find the areas of a whole collection of circle-polygon intersections. That is why I am referring to one of the circles as "current".

private Circle currentCircle;
private double currentSquareRadius;

Next comes the method for computing what we want to compute:

public final T computeValue(Circle circle, Shape shape) {
    initialize();
    processCircleShape(circle, shape);
    return getValue();
}

initialize() and getValue() are abstract. initialize() would set the variable that is keeping a total of the area to zero, and getValue() would just return the area. The definition for processCircleShape is

private void processCircleShape(Circle circle, final Shape cellBoundaryPolygon) {
    initializeForNewCirclePrivate(circle);
    if (cellBoundaryPolygon == null) {
        return;
    }
    PathIterator boundaryPathIterator = cellBoundaryPolygon.getPathIterator(null);
    double[] firstVertex = new double[2];
    double[] oldVertex = new double[2];
    double[] newVertex = new double[2];
    int segmentType = boundaryPathIterator.currentSegment(firstVertex);
    if (segmentType != PathIterator.SEG_MOVETO) {
        throw new AssertionError();
    }
    System.arraycopy(firstVertex, 0, newVertex, 0, 2);
    boundaryPathIterator.next();
    System.arraycopy(newVertex, 0, oldVertex, 0, 2);
    segment