How many unit squares of a square grid overlap a circle of given radius centered on the origin?

1.6k Views Asked by At

Consider, in the plane, the unit squares with corners having integral rectangular coordinates. Let $N_r$ be the number of these unit squares whose interior is intersected by a circle of radius $r$ centered on the origin. Counting these, I find the sequence $(N_r)_{r\in\mathbb{N}}=(0,4,12,20,28,28,44,52,60,68,68,84,92\ldots)$, or $({1\over 4}N_r)_{r\in\mathbb{N}}=(0,1,3,5,7,7,11,13,15,17,17,21,23,\ldots),$ neither of which appears in the OEIS, nor has searching turned up anything online.

Is there anything published about this sequence? Is it somehow obtainable from the known formulas for the solution of Gauss's circle problem or concerning circle lattice points (i.e., counting lattice points inside or on a circle of radius $r$)?

Here are some examples showing only the first quadrant: enter image description here

Apparently, $\lim_{r\to\infty}({1\over r}N_r)=8$ (but how to prove it?):

plots of N_r and N_r/r

(This is related to an older question, where a comment refers to algorithms for rasterizing a circle, but --although that turned out to be what the asker was looking for-- those algorithms don't seem relevant to the present question, as they generally seem to produce fewer than $N_r$ grid points.)

2

There are 2 best solutions below

6
On BEST ANSWER

Consider the function $$d^2(x, y) = x^2 + y^2 - r^2$$ which describes the signed squared distance between point $(x, y)$ and the circle of radius $r$ centered at origin. It is positive when point $(x, y)$ is outside the circle, negative when inside, and zero when point $(x, y)$ is on the circle.

Let's label lattice cells so that cell $(x, y)$ has vertices $(x, y)$, $(x+1, y)$, $(x+1, y+1)$, and $(x, y+1)$; i.e. lattice cells are labeled by the vertex with smallest coordinates.

The circle of radius $r$ intersects the interior of cell $(x, y)$ if and only if the $$\min\bigr( d(x,y), d(x+1, y), d(x, y+1), d(x+1, y+1) \bigr) \lt 0 \lt \max\bigr( d(x,y), d(x+1, y), d(x, y+1), d(x+1, y+1) \bigr)$$

(For intersecting the cell, i.e. including intersections with the vertices, use $\lt 0 \le$ above. Note that edges $x$ and $y$ belong to cell $(x, y)$, and edges $x+1$ and $y+1$ to cells $(x+1, y)$ and $(x, y+1)$, respectively.)

This applies to all circles of radius $r \in \mathbb{R}$, $r \ge 0$.

Here is a Python program that calculates the number of intersections (with cells if radius is positive, with cell interiors if radius is negative):

# SPDX-License-identifier: CC0-1.0

from math import floor, ceil, sqrt

def circle_intersects(radius, edges=True):

    if edges:
        rmax = floor(abs(radius))
    else:
        rmax = ceil(abs(radius)) - 1

    count = 0

    for y in range(0, rmax+1):  # 0 to rmax, inclusive

        dbase = radius*radius - y*y

        # Start at x outside the circle
        x = ceil(sqrt(dbase)) + 1
        d00 = dbase - x*x

        # Loop while cells intersect the circle
        while x >= 0:
            d01 = d00 - (2*x+1)
            d10 = d00 - (2*y+1)
            d11 = d10 - (2*x+1)

            dmin = min(d00, d01, d10, d11)
            dmax = max(d00, d01, d10, d11)

            if dmin < 0 and dmax > 0:
                count = count + 1
            elif edges and dmax > 0 and dmin == 0:
                count = count + 1
            elif dmin > 0:
                # completely inside the circle; next row
                break

            d00 += 2*x - 1
            x = x - 1            

    # Four identical quadrants
    return 4 * count

if __name__ == '__main__':
    from sys import argv, stderr, exit

    if len(argv) < 2 or argv[1] in ('-h', '--help', '/?'):
        stderr.write("\n")
        stderr.write("Usage: %s [ -h | --help | /? ]\n" % argv[0])
        stderr.write("       %s RADIUS [ RADIUS ... ]\n" % argv[0])
        stderr.write("\n")
        stderr.write("This program calculates the number of regular\n")
        stderr.write("rectangular integer lattice cells a circle or\n")
        stderr.write("radius RADIUS centered at origin intersects.\n")
        stderr.write("If the radius is negative, only intersections\n")
        stderr.write("with the interior of the cells are counted.\n")
        stderr.write("\n")
        exit(1)

    for arg in argv[1:]:
        radius = float(arg)

        count = circle_intersects(abs(radius), '-' not in arg)
        n = 8*ceil(abs(radius))-4 - count
        try:
            print("%s %d   %d %d  %f" % (arg.lstrip(" -"), count, 1*('-' not in arg), n, count/abs(radius)))
        except ZeroDivisionError:
            print("%s %d   %d %d" % (arg.lstrip(" -"), count, 1*('-' not in arg), n))

In the output, first column is the radius, second column the number of cells, third column is 1 if entire cells are counted and 0 if only cell interior, fourth column is the number of cells less than $8 r - 4$, and fifth column is the ratio between the number of cells and the radius (that should approach $8$).

The sequence $N_c(r)$ for nonnegative integer $r$ is 0, 4, 12, 20, 28, 36, 44, 52, 60, 68, 76, 84, 92, 100, 108, 116, 124, 132, 140, 148, 156, 164, 172, 180, 188, 196, 204, 212, 220, 228, 236, 244, 252, 260, 268, 276, 284, 292, 300, 308, 316, 324, 332, 340, 348, 356, 364, 372, 380, 388, 396, 404, 412, 420, 428, 436, 444, 452, 460, 468, 476, 484, 492, 500, 508, 516, 524, 532, 540, 548, 556, 564, 572, 580, 588, 596, 604, 612, 620, 628, 636, 644, 652, 660, 668, 676, 684, 692, 700, 708, 716, 724, 732, 740, 748, 756, 764, 772, 780, 788, 796, 804, 812, 820, 828, 836, 844, 852, 860, 868, 876, 884, 892, 900, 908, 916, 924, 932, 940, 948, 956, 964, 972, 980, 988, 996, and so on; in fact, $$N_c(r) = \begin{cases} 0, & r = 0 \\ 8 r - 4 & 1 \le r \in \mathbb{N} \\ \end{cases}$$

The sequence $N_i(r)$ is 0, 4, 12, 20, 28, 28, 44, 52, 60, 68, 68, 84, 92, 92, 108, 108, 124, 124, 140, 148, 148, 164, 172, 180, 188, 180, 196, 212, 220, 220, 228, 244, 252, 260, 260, 268, 284, 284, 300, 300, 308, 316, 332, 340, 348, 348, 364, 372, 380, 388, 380, 396, 404, 412, 428, 428, 444, 452, 452, 468, 468, 476, 492, 500, 508, 484, 524, 532, 532, 548, 548, 564, 572, 572, 580, 580, 604, 612, 612, 628, 628, 644, 644, 660, 668, 644, 684, 684, 700, 700, 708, 716, 732, 740, 748, 748, 764, 764, 780, 788, 780, 796, 804, 820, 820, 828, 836, 852, 860, 860, 868, 876, 892, 892, 908, 908, 916, 924, 940, 940, 948, 964, 964, 972, 988, 972, and so on.

Interestingly, $N_c(r) = N_i(r)$ for positive integer $r$ ($1 \le r \in \mathbb{N}$), if and only if $r$ is a nonhypotenuse number, i.e. when $r$ is not in OEIS A009003, for at least $r \le 10,000$.

In other words, when $r$ is a hypotenuse number (OEIS A009003), then $N_i(r) \lt N_c(r)$; otherwise $N_i(r) = N_c(r)$.

Apparently, $\lim_{r \to \infty}\left(\frac{N_i(r)}{r}\right) = 8$ (but how to prove it?)

Would it suffice to note that $N_i(r) \le N_c(r)$, and that $N_c(r) = 8 r - 4$?

Anyway, if you examine the cells in octant $0 \le y \lt x$, you'll see that there are exactly $r - 1$ cells intersecting with the circle of radius $r \in \mathbb{N}$, and up to $r - 1$ cells whose interiors intersect with that circle. (There is always exactly one cell in the diagonal $x = y \ge 0$.)

This is equivalent to $N_c(r) = 8(r - 1) + 4 = 8 r - 4$.

Because any row in this octant can have at most two cells intersecting with the circle (because of the tangent of the circle in this octant); and this happens on all columns expect when the circle intersects the integer $x$ coordinate between the two cells), and there are $\left\lceil r \sqrt{\frac{1}{2}} \right\rceil - 1$ rows, we know that $$8 \left\lceil r \sqrt{\frac{1}{2}} \right\rceil - 4 \le N_i(r) \le 8 \left\lceil r \sqrt{\frac{1}{2}} \right\rceil + 8\left\lceil r \left(1 - \sqrt{\frac{1}{2}}\right) \right\rceil - 4$$

For $1 \le r \in \mathbb{N}$, that upper limit simplifies to $8 r - 4$.

The logical reason why the number of cells tends to the upper limit is that the difference to the upper limit only occurs when the circle passes through a point with integer coordinates: at that point, there is one cell above and to the left of it, and one cell below and to the right of it, and therefore one less cell in this octant than the upper limit would state. This is also why there is a difference to the upper limit only when the radius is a hypotenuse number: only then are there points $(\chi, \gamma)$ on the circle with $0 \lt \gamma \lt \chi$ with $\chi \in \mathbb{N}$ and $\gamma \in \mathbb{N}$.

I am not good enough in maths to state all of the above in a form that would be acceptable as a proof, sorry. I only know this from rasterising circles, especially with antialiasing...

0
On

This is to supplement the accepted answer by sketching a "geometrical" argument that $N_r=8r-a(r)$, where $a(r)$ is the number of lattice points on a circle of radius $r$ centered on the origin.

Here we suppose horizontal and vertical gridlines connect all the lattice points that define the corners of the unit squares, and let an "overlap square" be any one of these unit squares whose interior is intersected by the circle.

First, by inspection it is clear that the circle touches exactly $8r$ gridlines (i.e. $2r$ gridlines per quadrant), noting that to touch a lattice point is to touch two gridlines simultaneously.

Second, there is exactly one overlap cell per touch point, because one new overlap square is entered upon passing through any touch point (which may be on either one or two gridlines).

Finally, the number of gridlines touched equals the number of touch points plus the number of lattice points touched (again because to touch a lattice point is to touch two gridlines simultaneously). Thus, $8r = N_r + a(r)$, and the required result follows.

Note that $a(r)=S(r^2)$, and both are described with a variety of algorithms in OEIS:

$a(n)$ is the number lattice points on a circle of radius $n$ (OEIS #A046109): $(1, 4, 4, 4, 4, 12, 4, 4, 4, 4, 12, 4, 4, 12, 4, 12, 4, 12, 4, 4, 12, 4, 4, 4, 4, 20,...)$

$S(n)$ is the number lattice points on a circle of radius $\sqrt{n}$ (OEIS #A004018): $(1, 4, 4, 0, 4, 8, 0, 0, 4, 4, 8, 0, 0, 8, 0, 0, 4, 8, 4, 0, 8, 0, 0, 0, 0, 12, 8, 0, 0,... )$