Area of intersection between circle and rectangle

Question

I'm looking for a fast way to determine the area of intersection between a rectangle and a circle (I need to do millions of these calculations).

A specific property is that in all cases the circle and rectangle always have 2 points of intersection.

Answer 1

Given 2 points of intersection:

0 vertices is inside the circle: The area of a circular segment

    XXXXX              -------------------    X     X               X            X Circular segment   X       X               XX        XX  +-X-------X--+              XXXXXXXX  |  X     X   | |   XXXXX    | 

1 vertex is inside the circle: The sum of the areas of a circular segment and a triangle.

    XXXXX                   XXXXXXXXX    X     X       Triangle ->X     _-X   X       X                 X   _-  X    X    +--X--+              X _-   X <- Circular segment     X   | X   |              X-  XXX      XXXXX    |              XXXX        |     |  

2 vertices are inside the circle: The sum of the area of two triangles and a circular segment

    XXXXX                   +------------X    X     X                  |      _--'/'X    X    +--X---    Triangle->|   _--   / X   X    |  X                 |_--     /XX <- Circular segment    X   +-X----              +-------XX     XXXXX                 Triangle^ 

3 vertices are inside the circle: The area of the rectangle minus the area of a triangle plus the area of a circular segment

    XXXXX    X  +--X+             XXX   X   |   X         -------XXX-----+ <- Triangle outside  X    |   |X        Rect ''.  XXX  |  X    +---+X                ''.  XX|    X         X                   ''. X <- Circular segment inside    X       X                       ^|X     X     X                         | X      XXXXX 

To calculate these areas:

• Most of the points you'll need to use can be found by finding the intersection of a line and a circle

• The areas you need to compute can be found by computing the area of a circular segment and computing the area of a triangle.

• You can determine if a vertex is inside the circle by calculating if its distance from the center is less than the radius.

Answer 2

I realize this was answered a while ago but I'm solving the same problem and I couldn't find an out-of-the box workable solution I could use. Note that my boxes are axis aligned, this is not quite specified by the OP. The solution below is completely general, and will work for any number of intersections (not only two). Note that if your boxes are not axis-aligned (but still boxes with right angles, rather than general quads), you can take advantage of the circles being round, rotate the coordinates of everything so that the box ends up axis-aligned and then use this code.

I want to use integration - that seems like a good idea. Let's start with writing an obvious formula for plotting a circle:

x = center.x + cos(theta) * radius y = center.y + sin(theta) * radius  ^ | |**###        ** | #*  #      *          * x |#  *  #    *           # y |#  *  #    *      +-----------------------> theta      *  #  *  #       *  #  *  #       *  #*  #        ***### 

This is nice, but I'm unable to integrate the area of that circle over x or y; those are different quantities. I can only integrate over the angle theta, yielding areas of pizza slices. Not what I want. Let's try to change the arguments:

(x - center.x) / radius = cos(theta) // the 1st equation theta = acos((x - center.x) / radius) // value of theta from the 1st equation y = center.y + sin(acos((x - center.x) / radius)) * radius // substitute to the 2nd equation 

That's more like it. Now given the range of x, I can integrate over y, to get an area of the upper half of a circle. This only holds for x in [center.x - radius, center.x + radius] (other values will cause imaginary outputs) but we know that the area outside that range is zero, so that is handled easily. Let's assume unit circle for simplicity, we can always plug the center and radius back later on:

y = sin(acos(x)) // x in [-1, 1] y = sqrt(1 - x * x) // the same thing, arguably faster to compute http://www.wolframalpha.com/input/?i=sin%28acos%28x%29%29+                 ^ y                |             ***|***     <- 1         ****   |   ****       **       |       **      *         |         *     *          |          * ----|----------+----------|-----> x    -1                     1 

This function indeed has an integral of pi/2, since it is an upper half of a unit circle (the area of half a circle is pi * r^2 / 2 and we already said unit, which means r = 1). Now we can calculate area of intersection of a half-circle and an infinitely tall box, standing on the x axis (the center of the circle also lies on the the x axis) by integrating over y:

f(x): integral(sqrt(1 - x * x) * dx) = (sqrt(1 - x * x) * x + asin(x)) / 2 + C // http://www.wolframalpha.com/input/?i=sqrt%281+-+x*x%29 area(x0, x1) = f(max(-1, min(1, x1))) - f(max(-1, min(1, x0))) // the integral is not defined outside [-1, 1] but we want it to be zero out there          ~         ~         |      ^  |         |      |  |         |   ***|***     <- 1         ****###|##|****       **|######|##|    **      *  |######|##|      *     *   |######|##|       * ----|---|------+--|-------|-----> x    -1   x0        x1      1 

This may not be very useful, as infinitely tall boxes are not what we want. We need to add one more parameter in order to be able to free the bottom edge of the infinitely tall box:

g(x, h): integral((sqrt(1 - x * x) - h) * dx) = (sqrt(1 - x * x) * x + asin(x) - 2 * h * x) / 2 + C // http://www.wolframalpha.com/input/?i=sqrt%281+-+x*x%29+-+h area(x0, x1, h) = g(min(section(h), max(-section(h), x1))) - g(min(section(h), max(-section(h), x0)))          ~         ~         |      ^  |         |      |  |         |   ***|***     <- 1         ****###|##|****       **|######|##|    **      *  +------+--+      *   <- h     *          |          * ----|---|------+--|-------|-----> x    -1   x0        x1      1 

Where h is the (positive) distance of the bottom edge of our infinite box from the x axis. The section function calculates the (positive) position of intersection of the unit circle with the horizontal line given by y = h and we could define it by solving:

sqrt(1 - x * x) = h // http://www.wolframalpha.com/input/?i=sqrt%281+-+x+*+x%29+%3D+h section(h): (h < 1)? sqrt(1 - h * h) : 0 // if h is 1 or above, then the section is an empty interval and we want the area integral to be zero                 ^ y                |               ***|***     <- 1         ****   |   ****         **       |       ** -----*---------+---------*------- y = h     *          |          * ----||---------+---------||-----> x    -1|                   |1 -section(h)          section(h) 

Now we can get the things going. So how to calculate the area of intersection of a finite box intersecting a unit circle above the x axis:

area(x0, x1, y0, y1) = area(x0, x1, y0) - area(x0, x1, y1) // where x0 <= x1 and y0 <= y1          ~         ~                              ~         ~         |      ^  |                              |      ^  |         |      |  |                              |      |  |         |   ***|***                              |   ***|***          ****###|##|****                          ****---+--+****      <- y1       **|######|##|    **                      **       |       **      *  +------+--+      *   <- y0            *         |         *     *          |          *                  *          |          * ----|---|------+--|-------|-----> x      ----|---|------+--|-------|-----> x         x0        x1                             x0        x1                  ^                |             ***|***         ****---+--+****      <- y1       **|######|##|    **      *  +------+--+      *   <- y0     *          |          * ----|---|------+--|-------|-----> x         x0        x1 

That's nice. So how about a box which is not above the x axis? I'd say not all the boxes are. Three simple cases arise:

• the box is above the x axis (use the above equation)
• the box is below the x axis (flip the sign of y coordinates and use the above equation)
• the box is intersecting the x axis (split the box to upper and lower half, calculate the area of both using the above and sum them)

Easy enough? Let's write some code:

float section(float h, float r = 1) // returns the positive root of intersection of line y = h with circle centered at the origin and radius r {     assert(r >= 0); // assume r is positive, leads to some simplifications in the formula below (can factor out r from the square root)     return (h < r)? sqrt(r * r - h * h) : 0; // http://www.wolframalpha.com/input/?i=r+*+sin%28acos%28x+%2F+r%29%29+%3D+h }  float g(float x, float h, float r = 1) // indefinite integral of circle segment {     return .5f * (sqrt(1 - x * x / (r * r)) * x * r + r * r * asin(x / r) - 2 * h * x); // http://www.wolframalpha.com/input/?i=r+*+sin%28acos%28x+%2F+r%29%29+-+h }  float area(float x0, float x1, float h, float r) // area of intersection of an infinitely tall box with left edge at x0, right edge at x1, bottom edge at h and top edge at infinity, with circle centered at the origin with radius r {     if(x0 > x1)         std::swap(x0, x1); // this must be sorted otherwise we get negative area     float s = section(h, r);     return g(max(-s, min(s, x1)), h, r) - g(max(-s, min(s, x0)), h, r); // integrate the area }  float area(float x0, float x1, float y0, float y1, float r) // area of the intersection of a finite box with a circle centered at the origin with radius r {     if(y0 > y1)         std::swap(y0, y1); // this will simplify the reasoning     if(y0 < 0) {         if(y1 < 0)             return area(x0, x1, -y0, -y1, r); // the box is completely under, just flip it above and try again         else             return area(x0, x1, 0, -y0, r) + area(x0, x1, 0, y1, r); // the box is both above and below, divide it to two boxes and go again     } else {         assert(y1 >= 0); // y0 >= 0, which means that y1 >= 0 also (y1 >= y0) because of the swap at the beginning         return area(x0, x1, y0, r) - area(x0, x1, y1, r); // area of the lower box minus area of the higher box     } }  float area(float x0, float x1, float y0, float y1, float cx, float cy, float r) // area of the intersection of a general box with a general circle {     x0 -= cx; x1 -= cx;     y0 -= cy; y1 -= cy;     // get rid of the circle center      return area(x0, x1, y0, y1, r); } 

And some basic unit tests:

printf("%f\n", area(-10, 10, -10, 10, 0, 0, 1)); // unit circle completely inside a huge box, area of intersection is pi printf("%f\n", area(-10, 0, -10, 10, 0, 0, 1)); // half of unit circle inside a large box, area of intersection is pi/2 printf("%f\n", area(0, 10, -10, 10, 0, 0, 1)); // half of unit circle inside a large box, area of intersection is pi/2 printf("%f\n", area(-10, 10, -10, 0, 0, 0, 1)); // half of unit circle inside a large box, area of intersection is pi/2 printf("%f\n", area(-10, 10, 0, 10, 0, 0, 1)); // half of unit circle inside a large box, area of intersection is pi/2 printf("%f\n", area(0, 1, 0, 1, 0, 0, 1)); // unit box covering one quadrant of the circle, area of intersection is pi/4 printf("%f\n", area(0, -1, 0, 1, 0, 0, 1)); // unit box covering one quadrant of the circle, area of intersection is pi/4 printf("%f\n", area(0, -1, 0, -1, 0, 0, 1)); // unit box covering one quadrant of the circle, area of intersection is pi/4 printf("%f\n", area(0, 1, 0, -1, 0, 0, 1)); // unit box covering one quadrant of the circle, area of intersection is pi/4 printf("%f\n", area(-.5f, .5f, -.5f, .5f, 0, 0, 10)); // unit box completely inside a huge circle, area of intersection is 1 printf("%f\n", area(-20, -10, -10, 10, 0, 0, 1)); // huge box completely outside a circle (left), area of intersection is 0 printf("%f\n", area(10, 20, -10, 10, 0, 0, 1)); // huge box completely outside a circle (right), area of intersection is 0 printf("%f\n", area(-10, 10, -20, -10, 0, 0, 1)); // huge box completely outside a circle (below), area of intersection is 0 printf("%f\n", area(-10, 10, 10, 20, 0, 0, 1)); // huge box completely outside a circle (above), area of intersection is 0 

The output of this is:

3.141593 1.570796 1.570796 1.570796 1.570796 0.785398 0.785398 0.785398 0.785398 1.000000 -0.000000 0.000000 0.000000 0.000000 

Which seems correct to me. You may want to inline the functions if you don't trust your compiler to do it for you.

This uses 6 sqrt, 4 asin, 8 div, 16 mul and 17 adds for boxes that do not intersect the x axis and a double of that (and 1 more add) for boxes that do. Note that the divisions are by radius and could be reduced to two divisions and a bunch of multiplications. If the box in question intersected the x axis but did not intersect the y axis, you could rotate everything by pi/2 and do the calculation in the original cost.

I'm using it as a reference to debug subpixel-precise antialiased circle rasterizer. It is slow as hell :), I need to calculate the area of intersection of each pixel in the bounding box of the circle with the circle to get alpha. You can sort of see that it works and no numerical artifacts seem to appear. The figure below is a plot of a bunch of circles with the radius increasing from 0.3 px to about 6 px, laid out in a spiral.