How to efficiently list all points of a sphere in an 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 not limited.

Let there be a radius R and a ball of this radius centered at X, i.e. the set of all grid points such that their Manhattan distance from X is R.

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

My question is: how can I list them in an efficient and easy way? By "enumerate" I mean an algorithm that, given N, X and R, will create a list of points that form this sphere (where a point is a list of its coordinates).

UPDATE: I called the metric first, mistakenly using the “Hamming distance”. I apologize to everyone who answered this question. Thanks to Steve Jessop for this.

Answer 1

Consider the smallest axis-aligned hypercube that limits the hypersphere and records the procedure for listing the grid points within the hypercube.

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

It is a simple and effective solution for small sizes. For example, for 2D, 20% of the points listed for the bounding box are discarded; for 6D, almost 90% of the hypercube points are discarded.

For higher dimensions, you will have to use a more complex approach: a cycle for each dimension (you may need to use a recursive function if the number of dimensions is variable). For each cycle, you will need to adjust the minimum and maximum values ​​depending on the values ​​of the already calculated grid components. Well, try doing this for 2D by listing the points of the circle, and once you understand it, generalizing the procedure to higher sizes will be pretty simple.

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

If you only want to list the points on the hypersurfaces of the cross-polytope, well, the solution is also very similar, you need to iterate over all the dimensions with the 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 each point generated in this way, you will have to consider all the variations leading to the negation of any of the components in order to cover all the faces and then extrude them with the vector X.

Update

Perl implementation for the case when 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

You can work recursively from the center, counting the zero distance once and working on symmetry. This Python implementation works on the bottom size of the stem vector and implements one one-dimensional slice at a time. You can also do the opposite, but this will mean iteration on partial hyperspheres. Although mathematically the same, the effectiveness of both approaches is highly language dependent.

If you know the power of the target space in advance, I would recommend writing an iterative implementation.

Below are the points on the hyper-LEGO R = 16 block in six dimensions of approximately 200 ms on my laptop. Of course, productivity decreases rapidly with size or larger areas.

 def lapp(lst, el): lst2 = list(lst) lst2.append(el) return lst2 def hypersphere(n, r, stem = [ ]): mystem = lapp(stem, 0) if 1 == n: ret = [ mystem ] for d in range(1, r+1): ret.append(lapp(stem, d)) ret.append(lapp(stem, -d)) else: ret = hypersphere(n-1, r, mystem) for d in range(1, r+1): mystem[-1] = d ret.extend(hypersphere(n-1, rd, mystem)) mystem[-1] = -d ret.extend(hypersphere(n-1, rd, mystem)) return ret 

(This implementation assumes that the hypersphere is centered at the origin. It would be easier to translate all the points after the transfer along the center coordinates).

Answer 3

Input: radius R, size D

• Generate all integer sections from R with cardinality ≤ D
• Rearrange it for each section without repeating
• For each permutation, twist all the signs

For example, the 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): # eg 4->[3,1] part = part + [0]*(D-len(part)) # eg [3,1]->[3,1,0] (padding) for perm in set(permutations(part)): # eg [1,3,0], [1,0,3], ... for point in product(*[ # eg [1,3,0], [-1,3,0], [1,-3,0], [-... ([-x,x] if x!=0 else [0]) for x in perm ]): yield point 

Demo radius = 4, size = 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 repeating, which is ineffective at all, but here it may not matter 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 effective because it doesn’t scale with a hyperwave, but simply scales with the hypersurface you are trying to list (maybe it doesn’t matter much except for very large radii: for example, you will save about one factor 100x speed if your radius equal to 100).