Overlap of disk and grid cell

circleseuclidean-geometrygeometry

I have a disk of radius $r$ centered at the origin. Space is divided into square grid cells of size $1$, with the origin $(0, 0)$ being at the center of one such grid cell. For each grid cell in the plane, I need to know its overlap $f$ with the disk, meaning the fraction (or area) of the cell that is enclosed by the circle of radius $r$.

The figure below illustrates the setup, with non-zero $f$ values written within each such cell.

Question

What is the general formula for $f$, given the radius $r$ and some cell specified by its integer coordinates $(i, j)$?

What I already know

I have the following for finding out whether a cell $(i, j)$ is completely outside or contained completely within the disk:
$$
\begin{align}
\biggl(|i| – \frac{1}{2}\biggr)^2 + \biggl(|j| – \frac{1}{2}\biggr)^2 &\ge r^2 \Rightarrow f(i, j) = 0 \quad \text{(outside)} \\
\biggl(|i| + \frac{1}{2}\biggr)^2 + \biggl(|j| + \frac{1}{2}\biggr)^2 &\le r^2 \Rightarrow f(i, j) = 1 \quad \text{(inside)} \\
\end{align}
$$
When none of these are true, we have $0 < f < 1$, for which I'm in need of a formula.

I've also realized that the symmetry generally allows me to consider only the positive quadrant, i.e. $i\rightarrow |i|$, $j\rightarrow |j|$. The problem still seems hard as the circle cuts the cells in different ways. I think a piecewise function is required.

The figure above is generated with the below Python script. Here $f$ is calculated (estimated) using the naive_count() function, which places a large number of points within the cell and checks whether each of them is within the disk or not.

import numpy as np
import matplotlib.pyplot as plt

# Configuration
gridsize = 10
r = 2.8

# Plot grid
nyquist = gridsize//2
for i in range(-nyquist + 1, nyquist + 1):
    plt.plot([i - 0.5, i - 0.5], [-nyquist + 0.5, nyquist - 0.5], 'k-')
for j in range(-nyquist + 1, nyquist + 1):
    plt.plot([-nyquist + 0.5, nyquist - 0.5], [j - 0.5, j - 0.5], 'k-')
plt.plot(0, 0, 'k.')

# Plot circle
theta = np.linspace(0, 2*np.pi, 1000)
plt.plot(r*np.cos(theta), r*np.sin(theta), 'C0')

def naive_count(r, i, j):
    N = 100
    count = 0
    for     x in np.linspace(i - 0.5, i + 0.5, N):
        for y in np.linspace(j - 0.5, j + 0.5, N):
            if x**2 + y**2 < r**2:
                count += 1
    return count/N**2

# Fill in each grid cell with a value
for     i in range(-int(np.ceil(r - 0.5)), int(np.ceil(r - 0.5)) + 1):
    for j in range(-int(np.ceil(r - 0.5)), int(np.ceil(r - 0.5)) + 1):
        I = abs(i)
        J = abs(j)
        if (I - 0.5)**2 + (J - 0.5)**2 >= r**2:
            # Completely outside
            continue
        if (I + 0.5)**2 + (J + 0.5)**2 <= r**2:
            # Completely within
            frac = 1
        else:
            # Partly within
            frac = naive_count(r, i, j)
        plt.text(i, j + 0.2, f'{frac:.2f}', ha='center', fontsize=9)
        plt.fill(
            [i + 0.5, i - 0.5, i - 0.5, i + 0.5],
            [j + 0.5, j + 0.5, j - 0.5, j - 0.5],
            color=[1 - frac, 1, 1 - frac], zorder=-1,
        )

# Finalize plot
plt.xlabel(r'$i$')
plt.ylabel(r'$j$')
plt.axis('image')
plt.xlim(-nyquist + 0.5, nyquist - 0.5)
plt.ylim(-nyquist + 0.5, nyquist - 0.5)
plt.tight_layout()
plt.savefig('grid.png')

Best Answer

Let us first assume that the cell intersects with the disc like you do, and let us also restrict ourselves to the (+,+) quadrant. The integral we need to solve is then $$ I = \int_{x_0}^{x_0 + 1} dx \int _{y_0}^{\sqrt{r^2-x^2}} dy \, 1 = \int_{x_0}^{x_0 + 1} dx \left( \sqrt{r^2-x^2} - y_0\right) \,, $$ where $(x_0,y_0)$ denote the lower left corner of the grid cell. However, this is not yet correct, since we need to (sometimes) restrict the integration region such that $y_0 < \sqrt{r^2-x^2} < y_0 + 1$ as shown in this sketch:

sketch of integration regions

The new integration limits are denoted by $\tilde x_0$ and $\tilde x_1$ and can be computed. We then just need the indefinite integral which is given by

$$ \int dx \left( \sqrt{r^2-x^2} - y_0\right) = \frac{1}{2} x \sqrt{r^2-x^2}-\frac{1}{2} r^2 \tan ^{-1}\left(\frac{x \sqrt{r^2-x^2}}{x^2-r^2}\right)-x y_0 \,. $$

If $\tilde x_0 > x_0$, we also need to add the missing rectangle to the integral, which is simply $\tilde x_0 - x_0$.

Finally, we put it together in a Python function to test it in your script:

def grid_disc_intersection(r, i, j):
    # Treat special case
    if i==j==0:
        if r >= 1/np.sqrt(2.):
            return 1.0
        elif r <= 0.5:
            return np.pi*r**2
        else:
            return np.sqrt(-1 + 4*r**2) + 2*r**2*(np.arccsc(2*r) - np.arctan(np.sqrt(-1 + 4*r**2)))
    minabsij = min(abs(i), abs(j))
    maxabsij = max(abs(i), abs(j))
    
    x0 = minabsij - 0.5
    y0 = maxabsij - 0.5
    
    # Compute integration region which is not neccessarily [x0, x0 + 1]
    x0t = x0
    if (y0 + 1) < r:
        xcross = np.sqrt(r**2 - (y0 + 1)**2)
        if x0t < xcross < x0t + 1:
            x0t = xcross
    
    x1t = x0 + 1
    if y0 < r:
        xcross = np.sqrt(r**2 - y0**2)
        if x0t < xcross < x0 + 1:
            x1t = xcross
    
    # Define indefinite integral
    def F(x):
        return (x*np.sqrt(r**2 - x**2))/2. - x*y0 - (r**2*np.arctan((x*np.sqrt(r**2 - x**2))/(-r**2 + x**2)))/2.
    
    # Apply fundamental theorem of calculus
    return F(x1t) - F(x0t) + (x0t - x0)

Related Solutions

Smallest hexagonal lattice with no more than one point per unit Cartesian cell, as a function of rotation angle

There are three regular tilings of the Euclidean plane, and they correspond to regular triangular lattice, regular square lattice, and regular hexagonal lattice.

In order for a lattice to be able to describe any point configuration with at most one point per cell, with pairwise distances greater than $r_\min$, the radius of the lattice cell circumcircle must be at most $2 r_\min$. In other words, the maximum distance between lattice cell vertices, or "cell diameter", must be $r_\min$ or less. Otherwise, a single lattice cell could contain two points (near vertices of the cell) with their distance greater than $r_\min$.

The properties of the three lattices with the same "cell diameter" (cell circumcircle diameter $d = r_\min$) are $$\begin{array}{l|ll|ll} \text{Lattice} & \text{Side length} & ~ & \text{Cell area} & ~ \\ \hline \text{Triangular} & a_△ = r_\min & \approx 1.0000 \, r_\min & A_△ = \frac{\sqrt{3}}{4} r_\min^2 & \approx 0.433 \, r_\min^2 \\ \text{Square} & a_□ = \frac{r_\min}{\sqrt{2}} & \approx 0.7071 \, r_\min & A_□ = \frac{1}{2} r_\min^2 & \approx 0.500 \, r_\min^2 \\ \text{Hexagonal} & a_⬡ = \frac{r_\min}{2} & \approx 0.5000 \, r_\min & A_⬡ = \frac{3 \sqrt{3}}{8} r_\min^2 & \approx 0.650 \, r_\min^2 \\ \end{array}$$ Maximum density of points is reached with a hexagonal lattice with $a = 2/\sqrt{3} \; r_\min$, where each cell has area $A_\min = \sqrt{3}/2 \, r_\min^2 \approx 0.866 \, r_\min^2$. From this we can estimate that for triangular lattices, at most half the cells will be populated with points; for square lattices, $1/\sqrt{3} \approx 0.5/0.866 = 58\%$ of cells will be populated; and for hexagonal lattices, $3/4 = 75\%$ of cells will be populated.

The following illustration shows these lattices, with shaded areas showing the region where existing points can limit where/whether the center white cell can contain a new point or not.
Image illustrating the three minimum/optimum lattices
Without subdividing the cell further, one needs to examine up to 24 triangular cells, 20 square cells, or 18 hexagonal cells.

If a larger side length $a$ is used for a lattice, some valid point configurations cannot be described using only one point per lattice cell. (That is, with larger $a$, not all valid point sets can be generated when limiting to at most a single point per cell.)

The problem at hand is to find $a_x$ for one lattice, given a larger than minimum other lattice $a_y$ and the relative rotation $\theta$ between the two. Because of the above, I believe this is the wrong way to try to fix the underlying problem.

What I would suggest, is to use the desired lattices (with arbitrary spacings, as long as they are not less than the ones stated above wrt. the minimum pairwise distance $r_\min$), but additionally verify and add each point to a minimal square lattice (which is computationally the simplest to implement, and likely most efficient on current hardware).

That is, an arbitrary lattices can be used to generate the points, but only when the minimal square lattice also allows each new point. While the minimal square lattice has smaller cells than the lattices used to generate points, it allows all possible point configurations (with pairwise distances greater than $r_\min$), and is therefore well suited for verification (that a point to be added does not break the minimum pairwise distance rule).

Calculating the Area of Polygons in the Intersection of Two Grids

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

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

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;