Computing the projection of an infinite 3D grid of points 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.
 A: 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.)
A: While this is an answer to the stated question in that it produces (almost) the desired result, I am posting it mainly to clarify my discussion with Rahul Narain in his comments, to show the the (GLSL) code I have written based on his description.
This is the “plus one dimension” generalization of Thomae's-function-over-an-interval as I interpreted Rahul Narain's description.
Notes: rad is currently arbitrary and should in principle be computed from the screen resolution. ceil(x0*fi) <= floor(x1*fi) answers “is there an integer in $[x_0i, x_1i]$?”. The return value's complications over being simply $1/i$ are just tweaks for the visual result I want.
const float rad = 0.001;
const int iter = int(0.5/rad);

float thomae(float x, float y) {
  float x0 = x - rad;
  float y0 = y - rad;
  float x1 = x + rad;
  float y1 = y + rad;
  for (int i = 1; i < iter; i++) {
    float fi = float(i);
    if (ceil(x0*fi) <= floor(x1*fi) && ceil(y0*fi) <= floor(y1*fi)) {
      return min(1.0, 10.0*pow(fi, -0.8));
    }
  }
  return 0.0;
}

thomae(x,y) computes a two-dimensional image. Therefore, to render it as surrounding the viewpoint, I project as if inside a cube:
    float spot;
    {
      vec3 v = vFixedOrientationPosition;
      vec3 a = abs(v);
      if (a.x > a.y && a.x > a.z) {
        spot = thomae(v.z/v.x, v.y/v.x);
      } else if (a.y > a.x && a.y > a.z) {
        spot = thomae(v.x/v.y, v.z/v.y);
      } else if (a.z > a.x && a.z > a.y) {
        spot = thomae(v.x/v.z, v.y/v.z);
      }
    }

v is the input direction vector; spot is the output brightness/nearness value used in later rendering. This branching code is the part I would like to get rid of by using a more natural “three-dimensional” calculation.
Another way in which this fails to be the desired result is that if rad spans more than one pixel, the visible large spot is a square on the projected cube (you can see this by observing that the condition in thomae() is for a rectangular region), whereas the “ideal” result would be cubical, I assume.
