[ACCEPTED]-Area of intersection between circle and rectangle-area

Accepted answer
Score: 64

Given 2 points of intersection:

0 vertices is inside 13 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 12 circle: The sum of the areas of a circular 11 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 10 sum of the area of two triangles and a circular 9 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 8 of the rectangle minus the area of a triangle 7 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 6 these areas:

Score: 14

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

I want to use integration 64 - that seems like a good idea. Let's start 63 with writing an obvious formula for plotting 62 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 61 integrate the area of that circle over x or 60 y; those are different quantities. I can 59 only integrate over the angle theta, yielding 58 areas of pizza slices. Not what I want. Let's 57 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 56 like it. Now given the range of x, I can 55 integrate over y, to get an area of the upper 54 half of a circle. This only holds for x in 53 [center.x - radius, center.x + radius] (other values will cause imaginary outputs) but 52 we know that the area outside that range 51 is zero, so that is handled easily. Let's 50 assume unit circle for simplicity, we can 49 always plug the center and radius back later 48 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 47 of pi/2, since it is an upper half of a unit 46 circle (the area of half a circle is pi * r^2 / 2 and 45 we already said unit, which means r = 1). Now we 44 can calculate area of intersection of a 43 half-circle and an infinitely tall box, standing 42 on the x axis (the center of the circle 41 also lies on the the x axis) by integrating 40 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 39 tall boxes are not what we want. We need 38 to add one more parameter in order to be 37 able to free the bottom edge of the infinitely 36 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 35 of the bottom edge of our infinite box from 34 the x axis. The section function calculates the 33 (positive) position of intersection of the 32 unit circle with the horizontal line given 31 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 30 we can get the things going. So how to calculate 29 the area of intersection of a finite box 28 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 27 nice. So how about a box which is not above 26 the x axis? I'd say not all the boxes are. Three 25 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 24 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 23 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 22 me. You may want to inline the functions 21 if you don't trust your compiler to do it 20 for you.

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

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

circles

Score: 5

I hope its not bad form to post an answer 14 to such an old question. I looked over the 13 above solutions and worked out an algorithm 12 which is similar to Daniels first answer, but 11 a good bit tighter.

In short, assume the 10 full area is in the rectangle, subtract 9 off the four segments in the external half 8 planes, then add any the areas in the four 7 external quadrants, discarding trivial cases 6 along the way.

pseudocde (my actual code 5 is only ~12 lines..)

find the signed (negative out) normalized distance from the circle center
to each of the infinitely extended rectangle edge lines,
ie.
d_1=(xcenter-xleft)/r
d_2=(ycenter-ybottom)/r
etc

for convenience order 1,2,3,4 around the edge. If the rectangle is not
aligned with the cartesian coordinates this step is more complicated but
the remainder of the algorithm is the same

If ANY d_i <=- 1 return 0

if ALL d_i >=  1 return Pi r^2

this leave only one remaining fully outside case: circle center in
an external quadrant, and distance to corner greater than circle radius:

for each adjacent i,j (ie. i,j=1,2;2,3;3,4;4,1)
     if d_i<=0 and d_j <= 0 and d_i^2+d_j^2 > 1 return 0

now begin with full circle area  and subtract any areas in the
four external half planes

Area= Pi r^2
for each  d_i>-1
     a_i=arcsin( d_i )  #save a_i for next step
     Area -= r^2/2 (Pi - 2 a_i - sin(2 a_i)) 

At this point note we have double counted areas in the four external
quadrants, so add back in:

for each adjacent i,j
   if  d_i < 1 and   d_j < 1  and d_i^2+d_j^2 < 1
       Area += r^2/4 (Pi- 2 a_i - 2 a_j -sin(2 a_i) -sin(2 a_j) + 4 sin(a_i) sin(a_j))

return Area

Incidentally, that last 4 formula for the area of a circle contained 3 in a plane quadrant is readily derived as 2 the sum of a circular segment, two right 1 triangles and a rectangle.

Enjoy.

Score: 3

The following is how to calculate the overlapping 19 area between circle and rectangle where 18 the center of circle lies outside the rectangle. Other 17 cases can be reduced to this problem.

The 16 area can be calculate by integrating the 15 circle equation y = sqrt[a^2 - (x-h)^2] + k where a is radius, (h,k) is 14 circle center, to find the area under curve. You 13 may use computer integration where the area 12 is divided into many small rectangle and 11 calculating the sum of them, or just use 10 closed form here.

alt text

And here is a C# source 9 implementing the concept above. Note that 8 there is a special case where the specified 7 x lies outside the boundaries of the circle. I 6 just use a simple workaround here (which 5 is not producing the correct answers in 4 all cases)

public static void RunSnippet()
{
    // test code
    double a,h,k,x1,x2;
    a = 10;
    h = 4;
    k = 0;
    x1 = -100;
    x2 = 100;

    double r1 = Integrate(x1, a, h, k);
    double r2 = Integrate(x2, a, h, k);

    Console.WriteLine(r2 - r1);

}

private static double Integrate(double x, double a,double h, double k)
{
    double a0 = a*a - (h-x)*(h-x);

    if(a0 <= 0.0){
        if(k == 0.0)
            return Math.PI * a * a / 4.0 * Math.Sign(x);
        else
            throw new Exception("outside boundaries");
    }

    double a1 = Math.Sqrt(a*a - (h-x)*(h-x)) * (h-x);
    double area = 0.5 * Math.Atan(a1 / ((h-x)*(h-x) - a*a))*a*a - 0.5 * a1 + k * x;
    return area;
}

Note: This problem is very similar 3 to one in Google Code Jam 2008 Qualification round problem: Fly Swatter. You can click on the 2 score links to download the source code 1 of the solution too.

Score: 2

Thanks for the answers,

I forgot to mention 9 that area estimatations were enough. That; s 8 why in the end, after looking at all the 7 options, I went with monte-carlo estimation 6 where I generate random points in the circle 5 and test if they're in the box.

In my case 4 this is likely more performant. (I have 3 a grid placed over the circle and I have 2 to measure the ratio of the circle belonging 1 to each of the grid-cells. )

Thanks

Score: 1

Perhaps you can use the answer to this question, where 9 the area of intersection between a circle 8 and a triangle is asked. Split your rectangle 7 into two triangles and use the algorithms 6 described there.

Another way is to draw a 5 line between the two intersection points, this 4 splits your area into a polygon (3 or 4 3 edges) and a circular segment, for both you should be able 2 to find libraries easier or do the calculations 1 yourself.

Score: 1

Here is another solution for the problem:

public static bool IsIntersected(PointF circle, float radius, RectangleF rectangle)
{

        var rectangleCenter = new PointF((rectangle.X +  rectangle.Width / 2),
                                         (rectangle.Y + rectangle.Height / 2));

        var w = rectangle.Width  / 2;
        var h = rectangle.Height / 2;

        var dx = Math.Abs(circle.X - rectangleCenter.X);
        var dy = Math.Abs(circle.Y - rectangleCenter.Y);

        if (dx > (radius + w) || dy > (radius + h)) return false;


        var circleDistance = new PointF
                                 {
                                     X = Math.Abs(circle.X - rectangle.X - w),
                                     Y = Math.Abs(circle.Y - rectangle.Y - h)
                                 };


        if (circleDistance.X <= (w))
        {
            return true;
        }

        if (circleDistance.Y <= (h))
        {
            return true;
        }

        var cornerDistanceSq = Math.Pow(circleDistance.X - w, 2) + 
                    Math.Pow(circleDistance.Y - h, 2);

        return (cornerDistanceSq <= (Math.Pow(radius, 2)));
}

0

More Related questions