How to efficiently enumerate all points of sphere in n-dimensional grid

Question

Say, we have an N-dimensional grid and some point X in it with coordinates (x1, x2, ..., xN). For simplicity we can assume that the grid is unbounded.

Let there be a radius R and a sphere of this radius with center in X, that is the set of all points in grid such that their manhattan distance from X is equal to R.

I suspect that their will be 2*N*R such points.

My question is: how do I enumerate them in efficient and simple way? By "enumerate" I mean the algorithm, which, given N, X and R will produce the list of points which form this sphere (where point is the list of it's coordinates).

UPDATE: Initially I called the metric I used "Hamming distance" by mistake. My apologies to all who answered the question. Thanks to Steve Jessop for pointing this out.

Answer 1

Consider the minimal axis-aligned hypercube that bounds the hypersphere and write a procedure to enumerate the grid points inside the hypercube.

Then you only need a simple filter function that allows you to discard the points that are on the cube but not in the hypersphere.

This is a simple and efficient solution for small dimensions. For instance, for 2D, 20% of the points enumerated for the bounding square are discarded; for 6D, almost 90% of the hypercube points are discarded.

For higher dimensions, you will have to use a more complex approach: loop over every dimension (you may need to use a recursive function if the number of dimensions is variable). For every loop you will have to adjust the minimal and maximal values depending on the values of the already calculated grid components. Well, try doing it for 2D, enumerating the points of a circle and once you understand it, generalizing the procedure to higher dimensions would be pretty simple.

update: errh, wait a minute, you want to use the Manhattan distance. Calling the cross polytope "sphere" may be correct but I found it quite confusing! In any case you can use the same approach.

If you only want to enumerate the points on the hyper-surface of the cross polytope, well, the solution is also very similar, you have to loop over every dimension with appropriate limits. For instance:

for (i = 0; i <= n; i++)
  for (j = 0; j + i <= n; j++)
    ...
       for (l = 0; l + ...+ j + i <= n; l++) {
         m = n - l - ... - j - i;
         printf(pat, i, j, ..., l, m);
       }

For every point generated that way, then you will have to consider all the variations resulting of negating any of the components to cover all the faces and then displace them with the vector X.

update

Perl implementation for the case where X = 0:

#!/usr/bin/perl

use strict;
use warnings;

sub enumerate {
    my ($d, $r) = @_;

    if ($d == 1) {
        return ($r ? ([-$r], [$r]) : [0])
    }
    else {
        my @r;
        for my $i (0..$r) {
            for my $s (enumerate($d - 1, $r - $i)) {
                for my $j ($i ? (-$i, $i) : 0) {
                    push @r, [@$s, $j]
                }
            }
        }
        return @r;
    }
}

@ARGV == 2 or die "Usage:\n  $0 dimension radius\n\n";
my ($d, $r) = @ARGV;
my @r = enumerate($d, $r);
print "[", join(',', @$_), "]\n" for @r;

Answer 2

Input: radius R, dimension D

• Generate all integer partitions of R with cardinality ≤ D
• For each partition, permute it without repetition
• For each permutation, twiddle all the signs

For example, code in python:

from itertools import *

# we have to write this function ourselves because python doesn't have it...
def partitions(n, maxSize):
    if n==0:
        yield []
    else:
        for p in partitions(n-1, maxSize):
            if len(p)<maxSize:
                yield [1] + p
            if p and (len(p)<2 or p[1]>p[0]):
                yield [ p[0]+1 ] + p[1:]

# MAIN CODE    
def points(R, D):
    for part in partitions(R,D):             # e.g. 4->[3,1]
        part = part + [0]*(D-len(part))      # e.g. [3,1]->[3,1,0]    (padding)
        for perm in set(permutations(part)): # e.g. [1,3,0], [1,0,3], ...
            for point in product(*[          # e.g. [1,3,0], [-1,3,0], [1,-3,0], [-...
                  ([-x,x] if x!=0 else [0]) for x in perm
                ]):
                yield point

Demo for radius=4, dimension=3:

>>> result = list( points(4,3) )

>>> result
[(-1, -2, -1), (-1, -2, 1), (-1, 2, -1), (-1, 2, 1), (1, -2, -1), (1, -2, 1), (1, 2, -1), (1, 2, 1), (-2, -1, -1), (-2, -1, 1), (-2, 1, -1), (-2, 1, 1), (2, -1, -1), (2, -1, 1), (2, 1, -1), (2, 1, 1), (-1, -1, -2), (-1, -1, 2), (-1, 1, -2), (-1, 1, 2), (1, -1, -2), (1, -1, 2), (1, 1, -2), (1, 1, 2), (0, -2, -2), (0, -2, 2), (0, 2, -2), (0, 2, 2), (-2, 0, -2), (-2, 0, 2), (2, 0, -2), (2, 0, 2), (-2, -2, 0), (-2, 2, 0), (2, -2, 0), (2, 2, 0), (-1, 0, -3), (-1, 0, 3), (1, 0, -3), (1, 0, 3), (-3, -1, 0), (-3, 1, 0), (3, -1, 0), (3, 1, 0), (0, -1, -3), (0, -1, 3), (0, 1, -3), (0, 1, 3), (-1, -3, 0), (-1, 3, 0), (1, -3, 0), (1, 3, 0), (-3, 0, -1), (-3, 0, 1), (3, 0, -1), (3, 0, 1), (0, -3, -1), (0, -3, 1), (0, 3, -1), (0, 3, 1), (0, -4, 0), (0, 4, 0), (0, 0, -4), (0, 0, 4), (-4, 0, 0), (4, 0, 0)]

>>> len(result)
66

(Above I used set(permutations(...)) to get permutations without repetition, which is not efficient in general, but it might not matter here due to the nature of the points. And if efficiency mattered, you could write your own recursive function in your language of choice.)

This method is efficient because it does not scale with the hypervolume, but just scales with the hypersurface, which is what you're trying to enumerate (might not matter much except for very large radii: e.g. will save you roughly a factor of 100x speed if your radius is 100).