1

Let there be two identical $n\times n$ grids $A$ and $B$, each consisting of unit cells. Assume grid $B$ is translated by $(x, y)$ relative to grid $A$ and then rotated by an angle $\theta$. How can I calculate (or approximate) the area of all the small polygons created through the intersection of these two grids?

EDIT. The motivation for this problem is that each cell of grid $B$ is assigned a value and I want to use them to assign area-averaged values to cells of $A$ that are fully or partly covered by $B$.

My current solution is to divide each cell of $B$ into $d\times d$ smaller cells of the same value, and then add $\frac{1}{d ^ 2}$ times the value of the smaller cell to the corresponding cell of $A$ where its center-point falls into. Though it provides a good approximation, this can be computationally intensive for larger grids.

Goodarz Mehr
  • 2,309

1 Answers1

1

enter image description here

Fig. 1 : The "big picture" displaying many different cases.

enter image description here

Fig. 2 : An octagonal case (Fig. 2 and Fig. 3 are created by the program below "playing on" translation vector $V$).

enter image description here

Fig. 3 : The case of a quadrilateral.

Your problem can be reduced to a classical issue : the determination of the intersection of 2 "filled" convex polygons $(A)$ and $(B)$ with vertices $A_1,A_2...$ and $B_1,B_2...$.

It is based on 2 principles (implemented as 2 functions in the Matlab program below)

  • ("getting new points by side intersections") Testing whether side $A_pA_{p+1}$ intersects side $B_qB_{q+1}$ (of course, most of them have no intersection). Build a list $(L)$ of all these intersections.

  • ("keeping inside vertices") Add to list $(L)$ all vertices $A_p$ which are inside polygon $(B)$ and, in a symmetrical way, all vertices $B_q$ which are inside polygon $(A)$.

Points of $(L)$ are the looked-for vertices of the intersection.

I have written (see below) a Matlab program in the case where polygons $(A)$ and $(B)$ are squares. $(A)$ has been taken WLOG to be the square $[0,1] \times [0,1]$ ; square $(B)$ is obtained from square $(A)$ by a rotation followed by a translation.

This program can be translated into another language ; you will just have to account for some "idiosyncrasies" of Matlab language I have used. For example, the last lines of the program use the "convhull" function ("convhull" = convex hull) yielding an anticlockwise ordering of points in $(L)$ and a computation of the common area.

Change the translation vector $V$ in order to see different cases.

clear all;close all;hold on;axis equal
%%%%%
function I=inters(A,Ap,B,Bp) ; %  intersection of line segments [AA'] and [BB']
I=[-1;-1]; % case where there is no intersection
AAp=Ap-A;BBp=Bp-B;
AB=B-A;AAp=Ap-A;
d=det([AAp,BBp]);
if d!=0
   L=det([AB,BBp])/d;M=-det([AAp,AB])/d;
   if (0<=L) && (L<=1) && (0<=M) && (M<=1) 
      I=A+L*AAp;
   end
end;
end;
%%%%%
function b=pointinsq(P,C); % boolean b=1 iff point P is inside square C
% criteria : all det(AP_k,AP_{k+1}) must be >0
   b=1;
   for k=1:4
      b=b && det([C(:,k)-P,C(:,k+1)-P])>=0;
   end;  
end; % 
%%%%%
t=pi/8;c=cos(t);s=sin(t);
mA=[0,1,1,0,0
    0,0,1,1,0]; % first square
    Ro=[c -s
        s c]; % rotation matrix
     V=[0.3;-0.2]; % transl. vector
mB=Ro*mA+V*ones(1,5); % second square
L=[];
for k=1:4 : % Adding to list the good vertices
   if pointinsq(mB(:,k),mA);
      L=[L,mB(:,k)]; % concatenation
   end
   if pointinsq(mA(:,k),mB);
      L=[L,mA(:,k)];
   end;
end;
for p=1:4
   for q=1:4
      I=inters(mA(:,p),mA(:,p+1),...
               mB(:,q),mB(:,q+1));
      if pointinsq(I,mA)
         L=[L,I];
      end;
   end;
end;
plot(mA(1,:),mA(2,:),'r');
plot(mB(1,:),mB(2,:),'b');
if size(L,2)>0, % if nonvoid intersection
   plot(L(1,:),L(2,:),'*k');
   [K,Area]=convhull(L(1,:),L(2,:));
   text(mean(L(1,:))-0.3,mean(Li(2,:)),['Area : ',num2str(Area)]);
   L=L(:,K);
   fill(L(1,:),L(2,:),'g');
end;
Jean Marie
  • 81,803