Consider the subset $S$ of $\mathbb{R}^3$ consisting of points whose coordinates are integers (compare Gaussian integers, Euclid's orchard). The view of $S$ from a perspective camera within the space has interesting structure; it has visible density variations that indicate planes intersecting the viewpoint with rational slopes, and the image can be interpreted as many different three-dimensional structures.
This image is actually of a finite subset of $S$, cubical and centered around the camera. (Increasing the number of points would narrow the width of the black empty regions.)
The question: I would like a practical-to-compute function for rendering the infinite-grid version of this image; that is, a function from a direction (and perhaps a small solid angle for the pixel size) to a brightness for the pixel. I'm not sure whether that result should be:
- the distance to the nearest point in $S$ in that direction, or
- the density of points in that direction.
As described in the Wikipedia article for Euclid's orchard, the two-dimensional zero-angle version of the function I'm looking for is Thomae's function; but that function is (a) giving a projection of the two-dimensional analogue of $S$ rather than $S$, (b) giving the view from an infinitely thin ray rather than over a small angle, and (c) is not practical to compute using floating-point arithmetic on a GPU.
Best Answer
A one-dimensional analogue of this question is to compute the maximum of Thomae's function in any given interval $[a,b]$. You can do this in time inversely proportional to the length of the interval, simply by iterating over $q = 1, 2, \dots$ and checking if $[qa, qb]$ contains an integer. This amounts to grouping 2D lattice points into parallel planes $x + y = q$ and checking them in order of increasing $q$.
You can generalize this pretty easily to 3 dimensions. An interval now corresponds to a subset $S$ of the image plane—the support of a pixel, say. As before, group the 3D lattice points into parallel planes at increasing distance from the origin. (What I originally had in mind was a direct translation of the 2D case: $x + y + z = q$ for the positive octant and its reflections for the other octants. But Kevin Reid's approach of $x = q$ for $\lvert x \rvert > \max(\lvert y \rvert, \lvert z \lvert)$ and so on for $y$ and $z$ also works, and is simpler to implement.) For $q = 1, 2, \ldots$, project $S$ onto the plane corresponding to $q$, and check if it contains a lattice point. If so, you have the nearest point in the given direction; if not, try $q+1$.
(It's worth mentioning that in all this, we're treating each pixel as if it were a discrete, square-shaped subset of the image plane, which is, strictly speaking, an incorrect and harmful model of image formation. However, in this case, I'll prefer Alvy Ray Smith's ire to the problem of doing correct sampling and reconstruction of an infinite lattice of point singularities.)