Overlap of disk and grid cell 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')

 A: 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)

