Path: utzoo!utgpu!news-server.csri.toronto.edu!cs.utexas.edu!swrinde!sdd.hp.com!zaphod.mps.ohio-state.edu!van-bc!ubc-cs!alberta!herald.usask.ca!ccu.umanitoba.ca!ccu.umanitoba.ca!lambert From: lambert@silver.cs.umanitoba.ca (Tim Lambert) Newsgroups: comp.graphics Subject: Re: Wanted: algorithm to compute intersection of square and circle Message-ID: Date: 26 Apr 91 01:30:01 GMT References: <1991Apr25.192400.8438@cfa203.harvard.edu> Sender: news@ccu.umanitoba.ca Organization: Department of Computer Science, University of Manitoba Lines: 85 In-Reply-To: dp@cfa.harvard.edu's message of 25 Apr 91 19:24:00 GMT >>>>> On 25 Apr 91 19:24:00 GMT, dp@cfa.harvard.edu (Dave Plummer) said: > I am looking for an algorithm that finds the area of intersection > between a circle and a square of arbitrary position and size. All right, I'll have a go at it. You could modify one of the standard polygon clipping algorithms to work on circles, to give you the perimeter of the resulting shape and compute the area from that, but a square and a circle seem simple enough that something just for this case would be nicer. First, let's rotate, scale and translate everything so that the circle has centre (0,0) and radius 1. Let the rectangle have corners (l,b) (left bottom) and (r,t) (right top). Now cut the rectangle up into four triangular pieces by connecting the four corners to the origin. (If the origin is not in the square, some of these pieces may be negative, but that's OK.) We'll add up the areas (possibly negative) of the intersection between these triangles and the circle. I'll describe how to work out the area of the intersection between the circle and the bottom triangle (corners (0,0) (l,b) (r,b)). The area of this triangle is -0.5*(r-l)*b. (Note that if b is positive this area is _negative_) Now let's look at where the line y=b intersects the circle. (Drawing a picture for each following case will probably help understand what is going on.) Case 1: b>1 or b<-1 -- no intersection between line and circle. Both (l,b) and (r,b) are outside the circle. The triangle-circle intersection is just a sector of the circle between angles atan2(b,l) and atan2(b,r). Area is atan2(b,r) - atan2(b,l) Note1: if b is positive, this is negative, which is correct. Note2: you could also calculate this area as acos(b^2+r*l/((sqrt(b^2+r^2)*sqrt(b^2+l^2)))) (and then fix the sign) Case 2: -1<=b<=1 -- intersection at (A,b) and (-A,b) where A=sqrt(1-b^2) Now -A<=A and l<=r so there are six possible orderings of -A,A,l,r Case 2a: -A<=l<=r<=A The triangle is entirely inside the circle. Area is -0.5*(r-l)*b Case 2b: A<=l<=r (l,b) and (r,b) are outside the circle. Intersection is a sector as in case 1. Case 2c: l<=r<=-A Same as 2b. Case 2d. -A<=l<=A<=r The intersection between the triangle and circle can be cut into two pieces: the triangle (0,0) (l,b) (A,b) and the sector between the angles atan2(b,A) and atan2(b,l). Area is -0.5*(A-l)*b + atan2(b,r) - atan2(b,A) Case 2e. l<=-A<=r<=A Similar to 2d Area is -0.5*(r+A)*b + atan2(b,-A) - atan2(b,l) Case 2f. l<=-A<=A<=r This is the complicated one. We have to cut it up into two sectors and a triangle. Area is atan2(b,-A) - atan2(b,l) -0.5*(A+A)*b + atan2(b,r) - atan2(b,A) Now code that all up into a function tri_circ(r,l,b) that computes the area between the bottom triangle and the circle. The area of the circle rectangle intersection is: tri_circ(r,l,b) + tri_circ(t,b,r) + tri_circ(l,r,t) + tri_circ(b,t,l) (I just rotate the picture through 90 degrees each time) Well, almost. I assumed r>=l in the above analysis so I'd better reflect about the y axis (which will multiply the area by -1) as well in the last two cases. So the true answer is: tri_circ(r,l,b) + tri_circ(t,b,r) - tri_circ(r,l,t) - tri_circ(t,b,l) If the intersection is empty all the positive and negative areas above will cancel out and the area will be zero (or very close, given floating point arithmetic). This is a good test that you got all the signs right. Tim P.S. Alert readers might be concerned that if (c,d) is in the third quadrant and (e,f) in the second, then atan2(d,c) - atan2(f,e) doesn't give the angle correctly (it's 2*pi too small). Fortunately this can't happen in the above algorithm, since we only compute angles between points on the same horizontal line.