Thursday, July 5, 2012

Intersection area of a circle and a rectangle


Sometimes there are problems in the world of programming, which are not always easy to solve. One of them is the calculation of the intersection area of a circle and a rectangle. There are many known solutions for this problem, some of them are described here.
some possible cases of circle/rectangle intersection
While the approach of searching for corners of the rectangle which lie within the circle, and then using area calculation based on triangles and circular segments may be the fastest and most accurate approach, it is very difficult to implement. Splitting up the rectangle into two triangles makes this implementation much easier. Many good and reliable examples for this case are out there in the Internet.
Nevertheless, it is not possible to adopt the algorithm to other shapes easily. This is the reason why I decided to play around with computer integration and the calculation of shape areas a little bit. 

Solving the problem with computer integration

Computer integration is only an approximation, but may be adequate for many cases. It is basically a form of numeric integration, where the area below the function to integrate is split up into many narrow rectangles. With this approach, there is no need to do symbolic integration of the function to integrate. Computer integration is also much slower than other intersection methods, but this won't really matter, since nowadays computers are really fast.
The big advantage of solving this problem with computer integration is that someone can calculate the intersection area of virtually any shape, assumed its lower and upper bounds can be represented by a function y = f(x).
approach of computer integration

Defining the functions

The first step is to define our rectangle and our circle as functions. We need a function for the upper bound and a function for the lower bound.
For the circle, the circle equation can be used, where r is the radius and xc and yc define the center position of the circle:  
the circle equation
Keep in mind that this circle equation has in fact two results for each x, since the root function theoretically returns a positive and a negative value, where one can be used for the upper function and the other for the lower function. We will have to care about this later in our source code.
The functions for the rectangle are simply defined by the upper and lower bounds of the rectangle, as our rectangle lies parallel to the axis.
Keep in mind, that the origin of the coordinate system is in the top left corner when programming.
Here is the code of the functions I used for defining circle and rectangle:
  1. // rect is a System.Drawing.RectangleF structure which defines our rectangle  
  2. // x is the x coordinate to calculate f(x) for  
  3. static double UpperRectangleFunction(RectangleF rect, decimal x)  
  4. {  
  5.     return rect.Top;  
  6. }  
  7.   
  8. static double LowerRectangleFunction(RectangleF rect, decimal x)  
  9. {  
  10.     return rect.Bottom;  
  11. }  
  12.   
  13. // m is a System.Drawing.PointF which represents the center of our circle  
  14. // r is the radius of our circle  
  15. // x is the x coordinate to calculate f(x) for  
  16. static double UpperCircleFunction(PointF m, float r, decimal x)  
  17. {  
  18.     return m.Y - Math.Sqrt((r * r) - Math.Pow(((double)x - m.X), 2));  
  19. }  
  20.   
  21. static double LowerCircleFunction(PointF m, float r, decimal x)  
  22. {  
  23.     return m.Y + Math.Sqrt((r * r) - Math.Pow(((double)x - m.X), 2));  
  24. }  

Determining integration boundaries

While this step is not always essential, it can improve performance in almost all cases, but also, this step is a little complex to accomplish. Take into consideration, that this step also protects our functions from values which are out of function range (for instance a X which is greater than the circle radius).  
Below is a plot illustrating this problem. The green lines indicate the optimal horizontal boundaries for our integration.
integration boundaries
In the case of circle and rectangle, someone has to get the horizontal edge of the rectangle (top or bottom) which is nearer at the center of the circle and then get the two x coordinates (lets call them xcl for the left one and xcr for the right one) on the circle for this y coordinate, which can again be done by using a modified circle function.
The left bound of our intersection area is the larger one of the left circle x coordinate (xcl) and the left bound of our rectangle.
The right bound of our intersection area is the smaller one of the right circle x coordinate (xcr) and the right bound of our rectangle.
Below is the source code to determine the left and right boundaries of the integration, including a code snippet to return 0 directly for the case the rectangle lies completely outside the circle.
  1. //Check whether the rectangle lies completely outside of the circle.   
  2. //Note: It is easier to check if a rectangle is outside another rectangle or  
  3. //circle than to check whether it is inside.  
  4. if ((rect.Bottom < m.Y && rect.Right < m.X)  
  5.      && (GetDistance(new PointF(rect.Bottom, rect.Right), m) > r) ||  
  6.    (rect.Top > m.Y && rect.Right < m.X)  
  7.      && (GetDistance(new PointF(rect.Top, rect.Right), m) > r) ||  
  8.    (rect.Bottom < m.Y && rect.Left > m.X)  
  9.      && (GetDistance(new PointF(rect.Bottom, rect.Left), m) > r) ||  
  10.   (rect.Top > m.Y && rect.Left > m.X)  
  11.      && (GetDistance(new PointF(rect.Top, rect.Left), m) > r))  
  12. {  
  13.     return 0; //Terminate fast  
  14. }  
  15.   
  16. //A variable storing the nearest horizontal edge of the rectangle.   
  17. double nearestRectangleEdge = 0;  
  18.   
  19. //Determine what is nearer to the circle center - the rectangle top edge or the rectangle bottom edge  
  20. if (Math.Abs(m.X - rect.Top) > Math.Abs(m.X - rect.Bottom))  
  21. {  
  22.     nearestRectangleEdge = rect.Bottom;  
  23. }  
  24. else  
  25. {  
  26.     nearestRectangleEdge = rect.Top;   
  27. }  
  28.   
  29. //The bounds of our integration  
  30. decimal leftBound = 0;  
  31. decimal rightBound = 0;  
  32.   
  33. if (m.Y >= rect.Top && m.Y <= rect.Bottom)  
  34. {  
  35.     //Take care if the circle's center lies within the rectangle.   
  36.     leftBound = RoundToDecimal(Math.Max(-r + m.X, rect.Left));  
  37.     rightBound = RoundToDecimal(Math.Min(r + m.X, rect.Right));  
  38. }  
  39. else if (r >= Math.Abs(nearestRectangleEdge - m.Y))  
  40. {  
  41.     //If the circle's center lies outside of the rectangle, we can choose optimal bounds.  
  42.     leftBound = RoundToDecimal(Math.Max(-Math.Sqrt(r * r - Math.Abs(Math.Pow(nearestRectangleEdge - m.Y, 2))) + m.X, rect.Left));  
  43.     rightBound = RoundToDecimal(Math.Min(Math.Sqrt(r * r - Math.Abs(Math.Pow(nearestRectangleEdge - m.Y, 2))) + m.X, rect.Right));  
  44. }  
The function GetDistance is simply a function which uses trigonometry to get the distance between two points in the coordinate plane:
  1. static double GetDistance(PointF p1, PointF p2)  
  2. {  
  3.     return Math.Sqrt(Math.Pow(p1.X - p2.X, 2) + Math.Pow(p1.Y - p2.Y, 2));  
  4. }  

The struggle with floating point types

Our readers may have already noticed a strange function with the name RoundToDecimal in my source code.
This function is simply a workaround to the inaccuracy of float and double data types regarding decimal places. This inaccuracy is normally not a big problem, but can yield strange results when used for comparisons, which we will need later when integrating.
The function RoundToDecimal simply rounds a float or double number to thousand parts. Here is the code:
  1. static decimal RoundToDecimal(double d)  
  2. {  
  3.     return Math.Round((decimal)d * 1000) / 1000;  
  4. }  

Putting it all together

Now we have got functions which define our shapes and also safe integration boundaries which are guaranteed to lie within both of the shapes, the circle and the rectangle.
All we have to do now is to choose a resolution to use for the computer integration and loop through the integration bounds. The resolution is the width of the tiny rectangles which are used for approximation by the program.
To integrate the area of one such rectangle for a coordinate x, calculate the distance between the y values of the curves which are nearer to the y coordinate of the intersection areas center and multiply it with the width of the rectangle, the resolution.  
In other words, calculate UpperRectangleFunction(x) and UpperCircleFunction(x), take the maximum value and calculate LowerRectangleFunction(x) and LowerCircleFunction and take the minimum value. Then subtract the first value from the second and multiply the result with the resolution.
calculation of a single section within the computer integration
To integrate the whole area of intersection, you just have to loop trough the boundaries we calculated before and sum up all rectangles.
This should be no problem now. Here is the source code:
  1. decimal Resolution = 0.01M;  
  2.   
  3. //Loop trough the intersection area and sum up the area  
  4. for (decimal i = leftBound + Resolution; i <= rightBound; i += Resolution)  
  5. {  
  6.     upperBound = Math.Max(UpperRectangleFunction(rect, i - Resolution / 2), UpperCircleFunction(m, r, i - Resolution / 2));  
  7.     lowerBound = Math.Min(LowerRectangleFunction(rect, i - Resolution / 2), LowerCircleFunction(m, r, i - Resolution / 2));  
  8.   
  9.     a += ((decimal)lowerBound - (decimal)upperBound) * Resolution;  
  10. }  
  11.   
  12. return (float)a;  

Test cases

Here are some simple test cases as proof of concept. As you easily can see, the algorithm gets inaccurate beyond the second decimal place, so it wont suit very accurate applications, but at least it was fun and interesting to implement. Keep in mind that this algorithm is easy to adopt to more complex shapes: you only have to re-define the functions and the boundary calculation.
Rectangle fits in circle,  intersection area has to be as large as the rectangle area:
Rectangle: x=0 y=0 width=10 height=10 Circle: x=5 y=5 r=7,0711 Rectangle Area: 100 Circle Area: 157,081073204053 Intersection Area: 100
Circle fits in rectangle, intersection area has to be equal to the circle area:
Rectangle: x=0 y=0 width=10 height=10 Circle: x=5 y=5 r=5 Rectangle Area: 100 Circle Area: 78,5398163397448 Intersection Area: 78,5405883789063
Circle fits halve into rectangle, intersection area has to be equal to the half circle area:
Rectangle: x=0 y=0 width=10 height=5 Circle: x=5 y=0 r=5 Rectangle Area: 50 Circle Area: 78,5398163397448 Intersection Area: 39,2702941894531
Circle fits one fourth into rectangle, intersection area has to be equal to quarter circle area:
Rectangle: x=0 y=0 width=5 height=5 Circle: x=5 y=0 r=5 Rectangle Area: 25 Circle Area: 78,5398163397448 Intersection Area: 19,6351470947266
Rectangle lies outside of the circle, intersection area has to be zero.
Rectangle: x=7 y=7 width=5 height=5 Circle: x=0 y=0 r=5 Rectangle Area: 25 Circle Area: 78,5398163397448 Intersection Area: 0

Conclusion and download

It can be seen that the approximation of complex problems with simple algorithms is a valuable option when writing computer programs, since someone can harness processor speed. Developing a solution which generates estimates can sometimes be very useful, especially when there is no time for complex approaches, which would be fast and accurate but are difficult to implement.

No comments:

Post a Comment