Area of intersection of a circle with a rectangle I want to find the area of a given circle that comes under the region of a given rectangle. I searched many posts on stackoverflow but they are not satisfying. I followed this post
http://www.eex-dev.net/index.php?id=100
But it doesn't seem accurate. I'd appreciate if someone provides me a good accurate solution atleast upto 10^-6 precision. Thanks.
Decompose your intersection area into polygons which are completely inside the circle, and circular segments formed by a chord and a part of the arc. Use e.g. the shoelace formula to compute the area of the polygoms, and the segment area formula for the segments. Compute these to whatever accuracy you want.

The decomposition will need to make many case distinctions. But choosing integration bounds and integration formulas isn't really any easier, in my opinion. There just is no simple way around all those case distinctions.
You put the circle at the origin. The idea is to divide the rectangle 
into four rectangles, and then replace all rectangles 
that do not fall on the first quadrant by an equivalent rectangle, but located on the first quadrant. Then you deal with each
case according to the number of vertex inside the circle. At the end all areas are given in terms of a single function F(U,V)
that arises from the evaluation of the intersection area in the second case,
which is the case when only on vertex is inside the ellipse. Then you sum the four areas. The matlab code is:

A=2;        % Circle of radius 2 (or an ellipse of semi axes A,B)
B=2;
L_x=0.5;    % The  width
L_y=2.75;   % The height  

x_1=0;     %bottom left corner (x_1,y_1)
y_1=-1;  

%This function makes the calculation  

int_area(x_1,y_1,L_x,L_y,A,B)  

%The function definition
function [suma]=int_area(x_1,y_1,L_x,L_y,A,B)  

x(1)=x_1;
y(1)=y_1;
% each one of the remaining vertex  

x(2)=x(1);
y(2)=y(1)+L_y;
x(3)=x(1)+L_x;
y(3)=y(1)+L_y;
x(4)=x(1)+L_x;
y(4)=y(1);  

% the center of the rectangle  

x_m=x(1)+L_x/2;
y_m=y(1)+L_y/2;
suma=0;  

% The original rectangle was divided in four rectangles
% with the new vertex coordinate of the closest vertex to the origin given by (a,b)
% according to the article http://www.dtic.mil/dtic/tr/fulltext/u2/410103.pdf
% a (x- coordinate ) b( y coordinate)
% c new width
% d new height  

for i=1:4  

a(i)=max([ 0 (-1)^( 1/2*(i^2-i) )  x_m - L_x/2 ]);
b(i)=max([ 0 (-1)^( 1/2(i^2+i-2) )y_m - L_y/2 ]);
c(i)=max([ 0 (-1)^( 1/2(i^2-i) )  x_m + L_x/2 - a(i) ]);
d(i)=max([ 0 (-1)^( 1/2(i^2+i-2) )*y_m + L_y/2 - b(i) ]);  

%in case the width and the height are different from 0 otherwise it contributes area zero
if (c(i) !=0 &amp;&amp; d(i) !=0 )

%this is the ellipse equation evaluated on each vertex
eq_1=(a(i)/A)^2 + (b(i)/B)^2;
eq_2=(a(i)/A)^2 + ((b(i)+d(i))/B)^2;
eq_3=((a(i)+c(i))/A)^2 + ((b(i)+d(i))/B)^2;
eq_4=((a(i)+c(i))/A)^2 + (b(i)/B)^2;  

% Area of intersection for each case according to the number of vertex inside the circle (or ellipse)
% S intersection area  

if( eq_1 >=1  &amp;&amp; eq_2 >=1 &amp;&amp; eq_3 >=1 &amp;&amp; eq_4 >=1)
S=0;        %case_I:  All vertex outside
end  

if( eq_1 &lt;1  &amp;&amp; eq_2 >=1 &amp;&amp; eq_3 >=1 &amp;&amp; eq_4 >=1)
S=AB/2 F(a(i)/A,b(i)/B);   %case_II: vertex I inside
end  

if( eq_1 &lt;1  &amp;&amp; eq_2 >=1 &amp;&amp; eq_3 >=1 &amp;&amp; eq_4 &lt;1)
S=AB/2 (F(a(i)/A,b(i)/B)-F((a(i)+c(i))/A,b(i)/B));   %case_III: vertex 1 and 4 inside
end  

if( eq_1 &lt;1  &amp;&amp; eq_2 &lt;1 &amp;&amp; eq_3 >=1 &amp;&amp; eq_4 >=1)
S=AB/2 (F(a(i)/A,b(i)/B)-F(a(i)/A,(b(i)+d(i))/B));   %case_IV: vertex 1 and 2 inside
end  

if( eq_1 &lt;1  &amp;&amp; eq_2 &lt;1 &amp;&amp; eq_3 >=1 &amp;&amp; eq_4 &lt;1)
S=AB/2 (F(a(i)/A,b(i)/B)-F((a(i)+c(i))/A,b(i)/B) - F(a(i)/A,(b(i)+d(i))/B));   %case_V: vertex 3 outside
end  

if( eq_1 &lt;1  &amp;&amp; eq_2 &lt;1 &amp;&amp; eq_3 &lt;1 &amp;&amp; eq_4 &lt;1)
S=c(i)*d(i);     %case_VI: all inside
end  

else
S=0;    %in case the width or the height of the new rectangle is zero
end  

suma=suma+S;   %the total area is suma
end  

end  

function [res] =  F (U,V)
  res=asin( sqrt(1-U^2)sqrt(1-V^2) -UV ) -Usqrt(1-U^2)-Vsqrt(1-V^2)+2*U*V;
end

Note: for some reason some '*' (multiplication symbol) doesn't appear on the screen, but matlab uses it. All after the sign '%' is a comment

Best Regards

Ed.