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

Tried implementing Lukas Kalbertodt's solution (which is, by far the best exact analytic solution I've found online), but found a few critical bugs. Created a testing ground (linked below, also extended to handle n-sided polygons), but the solution (and commented fixes) is here:

float signedEdgeDist(vec2 a, vec2 b, vec2 center, vec3 p) {
  // original answer didn't specify how signed distance was calculated,
  // had inconsistent commentary "positive if inside, positive if outside"
  // and left the definition of "inside" and "outside" ambiguous and confusing
  vec2 dir = normalize(b - a);
  vec2 n = a + dir * dot(p.xy - a, dir);
  float s = sign(dot(p.xy - n, n - center));
  // max -1 handles the case of the circle being entirely within the polygon
  // attempts to clamp to < 1 cause false overlaps when circle entirely
  // outside the polygon for reasons I don't understand 
  return max(-1., distance(p.xy, n) * s / p.z);
}

float angle(vec2 a, vec2 b) {
  return acos(clamp(dot(normalize(a), normalize(b)), -1.0, 1.0));
}

float segmentArea(float ed) {
  // bug in original answer, saturate(ed + 1.0) was used 
  float h = clamp(ed + 1.0, 0.0, 2.0); 
  return acos(1.0 - h) - (1.0 - h) * sqrt(2.0 * h - h * h);
}

float excessArea(float d0, float d1, float angle) {
  // bug in original answer, angle was interior angle, needs to be exterior
  angle = PI - angle; 
  vec2 i0 = vec2(-cos(asin(d0)), -d0);
  vec2 i1 = vec2(cos(angle + asin(-d1)), sin(angle + asin(-d1)));
  vec2 p = vec2(-cos(angle) / sin(angle) * d0 + d1 / sin(angle), -d0);

  if (dot(p, p) > 1.0) {
    // bug in original answer, was p.x < 0.0
    bool b0 = p.x < i0.x; 
    bool b1 = p.y > i1.y;
    if (b0 == b1)
      return 0.0;
    if (b0)
      return segmentArea(d1);
    return segmentArea(d0);
  }

  float triangle_area = 0.5 * distance(p, i0) * distance(p, i1) * sin(PI - angle);
  float segment_angle = asin(d0) - angle - asin(-d1) + PI;
  float segment_area = 0.5 * (segment_angle - sin(segment_angle));
  return triangle_area + segment_area;
}

float overlap(Polygon p, vec3 s) {
  float edges[p.sides];
  for (int i = 0; i < p.sides; i++) {
    edges[i] = signedEdgeDist(p.corners[i], p.corners[(i + 1) % p.sides], p.center, s);
  }
  float inv_area = 0.0;
  for (int j = 0; j < p.sides; j++) {
    inv_area += segmentArea(edges[j]);
    inv_area -= excessArea(edges[j], edges[(j + 1) % p.sides], p.angles[j]);
  }
  return clamp((PI - inv_area) / PI, 0.0, 1.0);
}

Functional implementation here (click and drag to move the circle and the degree of overlap determines the background and overlapping area shade): Circle-polygon analytic overlap

Answer 2

The following is an exact solution which (in my opinion) is a lot simpler than the two topmost answers. Full code can be found at the bottom of this answer.

The idea is to calculate the area of the circle that's not overlapped by the triangle:

• Take the sum of the circle area outside of each side of the triangle. These are always circular segments (or 0). Green areas in the image below.
• Subtract the areas which we "overcounted", i.e. where the circular segments overlap. Red areas in the image below.

In the end, wie just subtract this value from the total area of the circle to arrive at the final solution.

In detail

(The code is written in WGSL, but translating to other languages should be straight forward.)

Inputs and assumptions

Without loss of generality, this algorithm assumes:

• The circle is a unit circle centered at the origin.
• The triangle is defined by:
  • ed0, ed1, ed2: For each edge, the signed distance of the origin from that edge. (It's positive if the origin is "outside" that edge, and positive if inside. The edge is a tangent of the circle if the signed distance is 1 or -1.)
  • angle01, angle12, angle20: The internal angles of the triangle. Of course, they sum to 180°, so angle20 = 180° - angle01 - angle12. angleXY is the angle between edge X and Y.

In this image, the signed distances are represented by the colored lines. If the line is dashed, the distance is positive, negative otherwise. (Left: blue ≈ -0.15, orange ≈ -0.15, cyan ≈ -0.55; right: blue ≈ 0.4, orange ≈ 0.4, cyan ≈ -1.1).

These input parameters might seem strange, but:

• They can be easily calculated from the 2D vertices, if you have those.
• They make the calculation easier.
• Having this abstract of a representation does not require you to have 2D points in the first place. You can calculate these inputs from 3D points in a plane as well, for example.

Step 1: Calculate green areas

The area of a circular segment with radius 1 and height h is (see Wikipedia):

acos(1 - h) - (1 - h) * sqrt(1 - (1 - h)^2)

Our edge distances ed are not quite the same as the height, but we can simply convert it. This results in:

fn area_outside_edge(ed: f32) -> f32 {
    let h = clamp(ed + 1.0, 0.0, 2.0);
    return acos(1.0 - h) - (1.0 - h) * sqrt(1.0 - pow(1.0 - h, 2.0));
}

Step 2: Subtracting overcount

This is unfortunately a bit tricker. The overcounted area can take two forms (or be 0):

The second form is a special case that only occurs if the vertex is outside the circle. Also, the left case is unfortunately not a circular sector. It's area will be more tricky to calculate.

We will create a function that takes two edge distances and the internal angle of those two edge. We will then calculate the three relevant intersection points:

• p: intersection of both edges
• i0: intersection of given edge 0 with the circle
• i1: intersection of given edge 1 with the circle

The function assumes that edge 0 is completely horizontal. This is fine as the whole problem has a circular symmetry.

Calculating those intersection points involves some trigonometry and is tricky as we need to get the correct one of the two line-circle intersection points). I will not explain this here further as it would make the answer even longer and honestly, I kind of forgot the derivation. Sorry :(

(Geogebra playground for this problem)

fn overcounted_area(d0: f32, d1: f32, angle: f32) -> f32 {
    let i0 = vec2(-cos(asin(d0)), -d0);
    let i1 = vec2(cos(angle + asin(-d1)), sin(angle + asin(-d1)));
    let p  = vec2(-cos(angle) / sin(angle) * d0 + d1 / sin(angle), -d0);

Now let's deal with the general case: the potential deformed pizza slice. It's area is the sum of its triangle part and its circular segment.

    let triangle_area = 0.5 * abs(i0.x - p.x) * distance(p, i1) * sin(PI - angle);
    let circ_segment_angle = asin(d0) - angle - asin(-d1) + PI;
    let circ_segment_area = 0.5 * (circ_segment_angle - sin(circ_segment_angle));
    return triangle_area + circ_segment_area;
}

Yay! Now we only have to deal with the case that the vertex (p) lies outside the circle:

    // dot(p, p) is length(p) squared, which is the same if we only compare to 1.
    if dot(p, p) > 1.0 { 
        // Dealing with four different cases of where the vertex could lie.
        let b0 = p.x < 0.0;
        let b1 = p.y > i1.y;
        if b0 == b1 {
            return 0.0;
        }
        if b0 {
            return area_outside_edge(d1);
        } else {
            return area_outside_edge(d0);
        }
    }

(Again here, there are some complexities in this code that I don't quite remember anymore. It has to do with in what quadrant (spanned by the two triangle edges) the vertex lies.)

Final code

Here is the full code:

fn triangle_circle_intersection(
    ed0: f32, 
    ed1: f32, 
    ed2: f32,
    angle01: f32,
    angle12: f32,
) -> f32 {
    let angle20 = PI - angle01 - angle12;

    return PI - (
          area_outside_edge(ed0) 
        + area_outside_edge(ed1) 
        + area_outside_edge(ed2)
        - overcounted_area(ed0, ed1, angle01)
        - overcounted_area(ed1, ed2, angle12)
        - overcounted_area(ed2, ed0, angle20)
    );
}

fn area_outside_edge(ed: f32) -> f32 {
    let h = saturate(ed + 1.0);
    return acos(1.0 - h) - (1.0 - h) * sqrt(1.0 - pow(1.0 - h, 2.0));
}

fn overcounted_area(d0: f32, d1: f32, angle: f32) -> f32 {
    let i0 = vec2(-cos(asin(d0)), -d0);
    let i1 = vec2(cos(angle + asin(-d1)), sin(angle + asin(-d1)));
    let p  = vec2(-cos(angle) / sin(angle) * d0 + d1 / sin(angle), -d0);

    if dot(p, p) > 1.0 { 
        let b0 = p.x < 0.0;
        let b1 = p.y > i1.y;
        if b0 == b1 {
            return 0.0;
        }
        if b0 {
            return area_outside_edge(d1);
        } else {
            return area_outside_edge(d0);
        }
    }

    let triangle_area = 0.5 * abs(i0.x - p.x) * distance(p, i1) * sin(PI - angle);
    let circ_segment_angle = asin(d0) - angle - asin(-d1) + PI;
    let circ_segment_area = 0.5 * (circ_segment_angle - sin(circ_segment_angle));
    return triangle_area + circ_segment_area;
}

(If you are using this, please add comments! The code above does not contain any as it has this answer as context. But moving it into your code base, you should add comments!)

There are certainly still some areas for optimization:

• If you already have the three vertices, you don't need to calculate p in overcounted_area.
• If two area_outside_edge functions return 0, you don't need to calculate any overcount.
• You only need to calculate overcount for edge-pairs where both edges returned a non-0 value from area_outside_edge.
• If the circle is completely inside the triangle (ed0 > 1.0 && ed1 > 1.0 && ed2 > 1.0), just return PI.
• If the triangle is completely inside the circle, just return the area of the triangle.

---

Please note: I came up with this answer myself and it seemed to fully work. I have not validated it against the other solutions, though. I just wanted to post this solution here as I spend a long time figuring this out and I hope I can save other people some time. The code in this answer may contain some typos as I haven't tested it again. If you can confirm that this code works, please post a comment!

Answer 3

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);
    segmentType = boundaryPathIterator.currentSegment(newVertex);
    while (segmentType != PathIterator.SEG_CLOSE) {
        processSegment(oldVertex, newVertex);
        boundaryPathIterator.next();
        System.arraycopy(newVertex, 0, oldVertex, 0, 2);
        segmentType = boundaryPathIterator.currentSegment(newVertex);
    }
    processSegment(newVertex, firstVertex);
}

Let's take a second to look at initializeForNewCirclePrivate quickly. This method just sets the instance fields and allows the derived class to store any property of the circle. Its definition is

private void initializeForNewCirclePrivate(Circle circle) {
    currentCircle = circle;
    currentSquareRadius = currentCircle.getRadius() * currentCircle.getRadius();
    initializeForNewCircle(circle);
}

initializeForNewCircle is abstract and one implementation would be for it to store the circles radius to avoid having to do square roots. Anyway back to processCircleShape. After calling initializeForNewCirclePrivate, we check if the polygon is null (which I am interpreting as an empty polygon), and we return if it is null. In this case, our computed area would be zero. If the polygon is not null then we get the PathIterator of the polygon. The argument to the getPathIterator method I call is an affine transformation that can be applied to the path. I don't want to apply one though, so I just pass null.

Next I declare the double[]s that will keep track of the vertices. I must remember the first vertex because the PathIterator only gives me each vertex once, so I have to go back after it has given me the last vertex, and form an edge with this last vertex and the first vertex.

The currentSegment method on the next line puts the next vertex in its argument. It returns a code that tells you when it is out of vertices. This is why the control expression for my while loop is what it is.

Most of the rest of the code of this method is uninteresting logic related to iterating through the vertices. The important thing is that once per iteration of the while loop I call processSegment and then I call processSegment again at the end of the method to process the edge that connects the last vertex to the first vertex.

Let's look at the code for processSegment:

private void processSegment(double[] initialVertex, double[] finalVertex) {
    double[] segmentDisplacement = displacment2D(initialVertex, finalVertex);
    if (segmentDisplacement[0] == 0 && segmentDisplacement[1] == 0) {
        return;
    }
    double segmentLength = Math.sqrt(dotProduct2D(segmentDisplacement, segmentDisplacement));
    double[] centerToInitialDisplacement = new double[]{initialVertex[0] - getCurrentCircle().getCenterX(), initialVertex[1] - getCurrentCircle().getCenterY()};
    final double leftX = dotProduct2D(centerToInitialDisplacement, segmentDisplacement) / segmentLength;
    final double rightX = leftX + segmentLength;
    final double y = wedgeProduct2D(segmentDisplacement, centerToInitialDisplacement) / segmentLength;
    processSegmentStandardGeometry(leftX, rightX, y);
}

In this method I implement the steps to transform an edge into the standard geometry as described above. First I calculate segmentDisplacement, the displacement from the initial vertex to the final vertex. This defines the x axis of the standard geometry. I do an early return if this displacement is zero.

Next I calculate the length of the displacement, because this is necessary to get the x unit vector. Once I have this information, I calculate the displacement from the center of the circle to the initial vertex. The dot product of this with segmentDisplacement gives me leftX which I had been calling xi. Then rightX, which I had been calling xf, is just leftX + segmentLength. Finally I do the wedge product to get y as described above.

Now that I have transformed the problem into the standard geometry, it will be easy to deal with. That is what the processSegmentStandardGeometry method does. Let's look at the code

private void processSegmentStandardGeometry(double leftX, double rightX, double y) {
    if (y * y > getCurrentSquareRadius()) {
        processNonIntersectingRegion(leftX, rightX, y);
    } else {
        final double intersectionX = Math.sqrt(getCurrentSquareRadius() - y * y);
        if (leftX < -intersectionX) {
            final double leftRegionRightEndpoint = Math.min(-intersectionX, rightX);
            processNonIntersectingRegion(leftX, leftRegionRightEndpoint, y);
        }
        if (intersectionX < rightX) {
            final double rightRegionLeftEndpoint = Math.max(intersectionX, leftX);
            processNonIntersectingRegion(rightRegionLeftEndpoint, rightX, y);
        }
        final double middleRegionLeftEndpoint = Math.max(-intersectionX, leftX);
        final double middleRegionRightEndpoint = Math.min(intersectionX, rightX);
        final double middleRegionLength = Math.max(middleRegionRightEndpoint - middleRegionLeftEndpoint, 0);
        processIntersectingRegion(middleRegionLength, y);
    }
}

The first if distinguishes the cases where y is small enough that the edge may intersect the circle. If y is big and there is no possibility of intersection, then I call the method to handle that case. Otherwise I handle the case where intersection is possible.

If intersection is possible, I calculate the x coordinate of intersection, intersectionX, and I divide the edge up into three portions, which correspond to regions 1, 2, and 3 of the standard geometry figure above. First I handle region 1.

To handle region 1, I check if leftX is indeed less than -intersectionX for otherwise there would be no region 1. If there is a region 1, then I need to know when it ends. It ends at the minimum of rightX and -intersectionX. After I have found these x-coordinates, I deal with this non-intersection region.

I do a similar thing to handle region 3.

For region 2, I have to do some logic to check that leftX and rightX do actually bracket some region in between -intersectionX and intersectionX. After finding the region, I only need the length of the region and y, so I pass these two numbers on to an abstract method which handles the region 2.

Now let's look at the code for processNonIntersectingRegion

private void processNonIntersectingRegion(double leftX, double rightX, double y) {
    final double initialTheta = Math.atan2(y, leftX);
    final double finalTheta = Math.atan2(y, rightX);
    double deltaTheta = finalTheta - initialTheta;
    if (deltaTheta < -Math.PI) {
        deltaTheta += 2 * Math.PI;
    } else if (deltaTheta > Math.PI) {
        deltaTheta -= 2 * Math.PI;
    }
    processNonIntersectingRegion(deltaTheta);
}

I simply use atan2 to calculate the difference in angle between leftX and rightX. Then I add code to deal with the discontinuity in atan2, but this is probably unnecessary, because the discontinuity occurs either at 180 degrees or 0 degrees. Then I pass the difference in angle onto an abstract method. Lastly we just have abstract methods and getters:

protected abstract void initialize();

    protected abstract void initializeForNewCircle(Circle circle);

    protected abstract void processNonIntersectingRegion(double deltaTheta);

    protected abstract void processIntersectingRegion(double length, double y);

    protected abstract T getValue();

    protected final Circle getCurrentCircle() {
        return currentCircle;
    }

    protected final double getCurrentSquareRadius() {
        return currentSquareRadius;
    }

}

Now let's look at the extending class, CircleAreaFinder

public class CircleAreaFinder extends CircleShapeIntersectionFinder<Double> {

public static double findAreaOfCircle(Circle circle, Shape shape) {
    CircleAreaFinder circleAreaFinder = new CircleAreaFinder();
    return circleAreaFinder.computeValue(circle, shape);
}

double area;

@Override
protected void initialize() {
    area = 0;
}

@Override
protected void processNonIntersectingRegion(double deltaTheta) {
    area += getCurrentSquareRadius() * deltaTheta / 2;
}

@Override
protected void processIntersectingRegion(double length, double y) {
    area -= length * y / 2;
}

@Override
protected Double getValue() {
    return area;
}

@Override
protected void initializeForNewCircle(Circle circle) {

}

}

It has a field area to keep track of the area. initialize sets area to zero, as expected. When we process a non intersecting edge, we increment the area by R^2 Δθ/2 as we concluded we should above. For an intersecting edge, we decrement the area by y*length/2. This was so that negative values for y correspond to positive areas, as we decided they should.

Now the neat thing is if we want to keep track of the perimeter we don't have to do that much more work. I defined an AreaPerimeter class:

public class AreaPerimeter {

    final double area;
    final double perimeter;

    public AreaPerimeter(double area, double perimeter) {
        this.area = area;
        this.perimeter = perimeter;
    }

    public double getArea() {
        return area;
    }

    public double getPerimeter() {
        return perimeter;
    }

}

and now we just need to extend our abstract class again using AreaPerimeter as the type.

public class CircleAreaPerimeterFinder extends CircleShapeIntersectionFinder<AreaPerimeter> {

    public static AreaPerimeter findAreaPerimeterOfCircle(Circle circle, Shape shape) {
        CircleAreaPerimeterFinder circleAreaPerimeterFinder = new CircleAreaPerimeterFinder();
        return circleAreaPerimeterFinder.computeValue(circle, shape);
    }

    double perimeter;
    double radius;
    CircleAreaFinder circleAreaFinder;

    @Override
    protected void initialize() {
        perimeter = 0;
        circleAreaFinder = new CircleAreaFinder();
    }

    @Override
    protected void initializeForNewCircle(Circle circle) {
        radius = Math.sqrt(getCurrentSquareRadius());
    }

    @Override
    protected void processNonIntersectingRegion(double deltaTheta) {
        perimeter += deltaTheta * radius;
        circleAreaFinder.processNonIntersectingRegion(deltaTheta);
    }

    @Override
    protected void processIntersectingRegion(double length, double y) {
        perimeter += Math.abs(length);
        circleAreaFinder.processIntersectingRegion(length, y);
    }

    @Override
    protected AreaPerimeter getValue() {
        return new AreaPerimeter(circleAreaFinder.getValue(), perimeter);
    }

}

We have a variable perimeter to keep track of the perimeter, we remember the value of the radius to avoid have to call Math.sqrt a lot, and we delegate the calculation of the area to our CircleAreaFinder. We can see that the formulas for the perimeter are easy.

For reference here is the full code of CircleShapeIntersectionFinder

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];
    }

    private Circle currentCircle;
    private double currentSquareRadius;

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

    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);
        segmentType = boundaryPathIterator.currentSegment(newVertex);
        while (segmentType != PathIterator.SEG_CLOSE) {
            processSegment(oldVertex, newVertex);
            boundaryPathIterator.next();
            System.arraycopy(newVertex, 0, oldVertex, 0, 2);
            segmentType = boundaryPathIterator.currentSegment(newVertex);
        }
        processSegment(newVertex, firstVertex);
    }

    private void initializeForNewCirclePrivate(Circle circle) {
        currentCircle = circle;
        currentSquareRadius = currentCircle.getRadius() * currentCircle.getRadius();
        initializeForNewCircle(circle);
    }

    private void processSegment(double[] initialVertex, double[] finalVertex) {
        double[] segmentDisplacement = displacment2D(initialVertex, finalVertex);
        if (segmentDisplacement[0] == 0 && segmentDisplacement[1] == 0) {
            return;
        }
        double segmentLength = Math.sqrt(dotProduct2D(segmentDisplacement, segmentDisplacement));
        double[] centerToInitialDisplacement = new double[]{initialVertex[0] - getCurrentCircle().getCenterX(), initialVertex[1] - getCurrentCircle().getCenterY()};
        final double leftX = dotProduct2D(centerToInitialDisplacement, segmentDisplacement) / segmentLength;
        final double rightX = leftX + segmentLength;
        final double y = wedgeProduct2D(segmentDisplacement, centerToInitialDisplacement) / segmentLength;
        processSegmentStandardGeometry(leftX, rightX, y);
    }

    private void processSegmentStandardGeometry(double leftX, double rightX, double y) {
        if (y * y > getCurrentSquareRadius()) {
            processNonIntersectingRegion(leftX, rightX, y);
        } else {
            final double intersectionX = Math.sqrt(getCurrentSquareRadius() - y * y);
            if (leftX < -intersectionX) {
                final double leftRegionRightEndpoint = Math.min(-intersectionX, rightX);
                processNonIntersectingRegion(leftX, leftRegionRightEndpoint, y);
            }
            if (intersectionX < rightX) {
                final double rightRegionLeftEndpoint = Math.max(intersectionX, leftX);
                processNonIntersectingRegion(rightRegionLeftEndpoint, rightX, y);
            }
            final double middleRegionLeftEndpoint = Math.max(-intersectionX, leftX);
            final double middleRegionRightEndpoint = Math.min(intersectionX, rightX);
            final double middleRegionLength = Math.max(middleRegionRightEndpoint - middleRegionLeftEndpoint, 0);
            processIntersectingRegion(middleRegionLength, y);
        }
    }

    private void processNonIntersectingRegion(double leftX, double rightX, double y) {
        final double initialTheta = Math.atan2(y, leftX);
        final double finalTheta = Math.atan2(y, rightX);
        double deltaTheta = finalTheta - initialTheta;
        if (deltaTheta < -Math.PI) {
            deltaTheta += 2 * Math.PI;
        } else if (deltaTheta > Math.PI) {
            deltaTheta -= 2 * Math.PI;
        }
        processNonIntersectingRegion(deltaTheta);
    }

    protected abstract void initialize();

    protected abstract void initializeForNewCircle(Circle circle);

    protected abstract void processNonIntersectingRegion(double deltaTheta);

    protected abstract void processIntersectingRegion(double length, double y);

    protected abstract T getValue();

    protected final Circle getCurrentCircle() {
        return currentCircle;
    }

    protected final double getCurrentSquareRadius() {
        return currentSquareRadius;
    }

Anyway, that is my description of the algorithm. I think it is nice because it is exact and there aren't really that many cases to check.

Answer 4

For future reference, here is some Matlab code I've written to solve this problem:

https://uk.mathworks.com/matlabcentral/fileexchange/126645-intersection-of-polygon-and-circle

It more or less implements Brian Moths' answer.

Here is the code:

function [A, pointList, isArc, AP] = polygonCircleIntersection(R, P)
% polygonCircleIntersection - calculates area of intersection of circle
% of radius R centered in origin with a 2D polygon with vertices in P, and
% provides shape of intersection
% 
% Inputs:
%   R : radius of circle [scalar]
%   P : 2 X nV matrix of polygon vertices. Each column is a vertex
% Vertices must be in counter-clockwise order or things will go wrong
%
% Outputs:
%   A : area of intersection
%   pointList : list of points on the perimiter of the intersection
%   isArc : array of logicals. isArc(i) is true if the segment between 
%   i-1 and i is an arc.
%   AP : the area of the polygon
%
% Used for FOV analysis
%
% See also: intersect, polyshape, triangulation
%
% Author: Simão da Graça Marto
% e-mail: simao.marto@gmail.com
% Date: 22/03/2023

nV = size(P, 2);
AP = 0;
for i = 1:nV
    AP = AP + shoelace(P(:,i), P(:,mod(i,nV)+1));
end
% SANITY CHECK
% order of points
if(AP<0) 
    error("Polygon must be in counter-clockwise order");
    %If this just means the polygon is not visible, uncomment the
    %following and remove the error:
%     A = 0;
%     pointList = [];
%     isArc = [];
%     return
end
nT2 = sum(P.^2);
isOutside = nT2>R^2;
%1st edge case: all points inside, so polygon fully inside:
if(~any(isOutside))
    A = 0;
    for i = 1:nV
        A = A + shoelace(P(:,i), P(:,mod(i,nV)+1));
    end
    pointList = P;
    isArc = false(1,nV);
    return;
end
% We must start from an outside vertex, so cycle vertices to be in a correct order
shiftK = 1-find(isOutside,1,"first");
P = circshift(P,shiftK, 2);
isOutside = circshift(isOutside,shiftK, 2);
% compute intersection
pointList = []; %list of points forming shape
isArc = [true]; % indicates if i-1 to i is an arc. if not, it's a line segment
outsideCircle = true;
iP = 1;
for i = 1:nV
    x0 = P(:,i);
    d = P(:,mod(i,nV)+1)-x0;
    %segment circle intersection
    [xI, onCircle] = segmentCircleIntersection(x0, d, R);
    pointList = [pointList xI];

    if(~isempty(xI))
        %if point iP is outside circle, then segment from previous point to
        %this one must be an arc
        isArc(iP) = outsideCircle;
        %point iP to iP+1 is always a line segment, by construction
        isArc(iP+1) = false;
    
        outsideCircle = onCircle(2); %if segment goes out of circle, we are now outside circle, otherwise we are inside
        iP = iP + 2; %update point index
    end
end
%edge cases: triangle perimeter fully outside, but is the circle
%inside or outside triangle?
if(isempty(pointList) && all(isOutside))
    polygonContainsOrigin = inpolygon(0,0,P(1,:), P(2,:));
    
    if(polygonContainsOrigin) %2nd edge case: circle fully inside triangle
        A = pi*R^2;
        pointList = [1;0]; %to allow drawing
    else %3rd edge case: triangle fully outside circle
        A = 0;
        pointList = [];
        isArc = [];
    end
    return;
end
%SANITY CHECKS
if(isempty(pointList) && any(isOutside) && ~all(isOutside))
    error("there are no intersections, but some points are inside and others outside")
end
if(any(isArc & circshift(isArc, 1)))
    error("there can't be two arcs in a row")
end
nI = size(pointList, 2);
%compute area as a sum of triangles (shoelace) and arcs
A = 0;
for i = 1:nI
    iPrev = mod(i-2,nI)+1;
    xPrev = pointList(:, iPrev);
    xi = pointList(:, i);
    if(isArc(i))
        thPrev = atan2(xPrev(2),xPrev(1));
        thi = atan2(xi(2),xi(1));
        dth = wrapTo2Pi(thi-thPrev);
        if(dth == 2*pi), dth = 0; end
        A = A + dth*R^2/2;
    else
        A = A + shoelace(xPrev, xi); %can be negative. That is correct.
    end
end
% SANITY CHECK
if(A<0 || A>pi*R^2 || A > AP)
    error("invalid area")
end
end
%Always returns 2 or 0 points. If tangent, none are returned because it
%would have no effect on the area. If segment starts or ends inside circle,
%start or end point are returned. By circle understand it to include its
%interior, so intersection is the segment between the two points returned
function [xI, onCircle] = segmentCircleIntersection(x0, d, R)
    a = sum(d.^2);
    b = 2*x0'*d;
    c = sum(x0.^2) - R^2; %could use precomputed nT2 instead
    D = b^2 - 4*a*c;
    if(D<=0) % line never goes in the circle
        xI = [];
        onCircle = [];
        return
    end
    %intersection points (along segment coordinate)
    ll = (-b + [-1 1]*sqrt(D) )/(2*a);
    % SANITY CHECK
    if( ll(1) >= ll(2) )
        error("A mathematical impossibility. Did you give this thing complex numbers or something?")
    end
    if( ll(2)<0 || ll(1)>1) %intersection fully outside segment
        xI = [];
        onCircle = [];
        return
    end
    %compute intersection
    ll = max(0, ll);
    ll = min(1, ll);
    
    xI = x0 + ll.*d;
    onCircle = [ll(1)>0 ll(2)<1];
end
%Shoelace formula for area of a straight segment with edges connected to
%origin, forming a triangle.
%can be negative if area is meant to be subtracted
function A = shoelace(p1, p2)
    A = (p1(1)*p2(2) - p1(2)*p2(1))/2;
end

Answer 5

If just one of the triangle's line segments intersects the circle, the pure math solution isn't too hard. Once you know when the two points of intersection are, you can use the distance formula to find the chord length.

According to these equations:

ϑ = 2 sin⁻¹(0.5 c / r)
A = 0.5 r² (ϑ - sin(ϑ))

where c is the chord length, r is the radius, ϑ becomes the angle through the center, and A is the area. Note that this solution breaks if more than half the circle is cut off.

It's probably not worth the effort if you just need an approximation, since it makes several assumptions about what the actual intersection looks like.

Answer 6

If you want an exact solution (or at least as exact as you can get using floating-point arithmetic) then this is going to involve a lot of legwork, because there are so many cases to consider.

I count nine different cases (categorized in the figure below by the number of vertices of the triangle inside the circle, and the number of edges of the triangle that intersect or are contained in the circle):

(However, this kind of enumeration of geometric cases is well known to be tricky, and it wouldn't surprise me at all if I missed one or two!)

So the approach is:

1. Determine for each vertex of the triangle if it's inside the circle. I'm going to assume you know how to do that.

2. Determine for each edge of the triangle if it intersects the circle. (I wrote up one method here, or see any computational geometry book.) You'll need to compute the point or points of intersection (if any) for use in step 4.

3. Determine which of the nine cases you have.

4. Compute the area of the intersection. Cases 1, 2, and 9 are easy. In the remaining six cases I've drawn dashed lines to show how to partition the area of intersection into triangles and circular segments based on the original vertices of the triangle, and on the points of intersection you computed in step 2.

This algorithm is going to be rather delicate and prone to errors that affect only one of the cases, so make sure you have test cases that cover all nine cases (and I suggest permuting the vertices of the test triangles too). Pay particular attention to cases in which one of the vertices of the triangle is on the edge of the circle.

If you don't need an exact solution, then rasterizing the figures and counting the pixels in the intersection (as suggested by a couple of other respondents) seems like a much easier approach to code, and correspondingly less prone to errors.

Answer 7

I'm almost a year and a half late, but I thought maybe people will be interested in code here that I wrote which I think does this correctly. Look in function IntersectionArea near the bottom. The general approach is to pick off the convex polygon circumscribed by the circle, and then deal with the little circular caps.

Answer 8

Since your shapes are convex, you can use Monte Carlo area estimation.

Draw a box around the circle and triangle.

Choose random points in the box and keep a count of how many fall in the circle, and how many fall in both the circle and triangle.

Area of Intersection ≅ Area of circle * # points in circle and triangle / # points in circle

Stop choosing points when the estimated area doesn't change by more than a certain amount over a certain number of rounds, or just choose a fixed number of points based on the area of the box. The area estimate should converge pretty fast unless one of your shapes has very little area.

Note: Here's how you determine if a point is in a triangle: Barycentric coordinates

Answer 9

I think you shouldn't approximate circle as some set of triangles, instead of that you can approximate it's shape with a polygon. The naive algorithm can look like:

1. Convert you circle to polygon with some desired number of vertices.
2. Calculate the intersection of two polygons (converted circle and a triangle).
3. Calculate square of that intersection.

You can optimize this algorithm by combining step 2 and step 3 into single function.

Read this links:
Area of convex polygon
Intersection of convex polygons

Answer 10

Assuming you're talking integer pixels, not real, the naive implementation would be to loop through every pixel of the triangle and check the distance from the circle's center against its radius.

It's not a cute formula, or particularly fast, but it does get the job done.