# 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 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.)

-
Thanks! This seems to be doing what I want, though it may still be too slow (due to the need for high iterations) for real-time rendering. I'm working on getting the right mapping from a direction vector onto the inputs to the two-argument function; would you happen to have a canned solution for that? – Kevin Reid Sep 26 '11 at 17:20
@Kevin: I'm sorry, which two-argument function do you mean? – Rahul Sep 26 '11 at 17:50
Thomae's function has one argument and can be seen as a one-dimensional image of two-dimensional space; the generalization to 3 dimensions you refer to has two arguments as it is a two-dimensional image. – Kevin Reid Sep 26 '11 at 18:07
@Kevin: Oh, that would be whatever your projection matrix is, wouldn't it? For each pixel in your image, you determine the view direction vector, and start checking lattice points; when you find one, that determines the colour of that pixel. – Rahul Sep 26 '11 at 18:31
I don't understand what the generalization you propose is, then. As I said, Thomae's function is itself an image in one less dimension; my generalization is a two-dimensional “square” image, and it is an additional step to map it onto the “spherical” space of direction vectors. Mapping it as if it were the texture on the faces of a cube works, but is not elegant. If this remains unclear, I'll post my current algorithm as an answer. – Kevin Reid Sep 30 '11 at 17:25

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;

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.

-
I'm a little embarrassed to be replying to this after months have passed, but here goes anyway... First, I don't think there's any way to get rid of the branching on x, y, and z, given that those are the symmetry axes of the lattice itself that we're trying to render. The solution I had in mind actually branched on the octant that v lay in, but your approach is better. I think your code is almost perfect as is, apart from the need to define an ad-hoc rad tolerance. For that, there are a couple of different things you could do: – Rahul Dec 20 '11 at 14:09
(1) If you want to guarantee that a lattice point appears as a single pixel, you have to test the point not simply against the view vector but against a rectangular cone around it. Imagine each pixel to be a little square, project it out into space, and see if the lattice point lies in that region. Equivalently, but more easily, project the lattice point onto the image plane and see if it lies within the pixel's square. – Rahul Dec 20 '11 at 14:14
(2) Alternatively, if you want the lattice point to have real size (like the "ideal" result you mention), you could check if the view vector intersects some geometry centered around the lattice point. You could do a cube, but a sphere is much easier to test against and would probably look prettier too. – Rahul Dec 20 '11 at 14:16
Thanks for the clarification and extra info. I’ve accepted your answer since it contains the core principle, but could you edit it with a little more of the information that has come out in the comments? – Kevin Reid Dec 20 '11 at 14:23
Done. (more characters) – Rahul Dec 20 '11 at 14:52