Recently I had to find the neighbors in a 3 by 3 by 3 grid around a given point in 3D space. J.D. Smith’s Dimensional Juggling Tutorial has a lot of similar use of the REFORM and REBIN routines, but mostly for the purpose of “inflating” a dimension of an array. This problem called for cyclically repeating the same values in different patterns.

I needed the nearest neighbors (i.e. inside a 3 by 3 by 3 window), but other window sizes might be needed.

n = 3
offsets = lindgen(n) - (n - 1) / 2

So offsets is just [-1, 0, 1]. The magic is:

xoffsets = reform(rebin(offsets, n, n^2), n^3)
yoffsets = reform(rebin(offsets, n^2, n), n^3)
zoffsets = rebin(offsets, n^3)

which does the REFORM and REBIN in the opposite order from the Dimensional Juggling Tutorial. You expand into an extra dimension, then get rid of that extra dimension which creates the repeating patterns. For example, try

IDL> print, reform(rebin([-1, 0, 1], 3, 9), 27)
-1 0 1 -1 0 1 -1 0 1
-1 0 1 -1 0 1 -1 0 1
-1 0 1 -1 0 1 -1 0 1
IDL> print, reform(rebin([-1, 0, 1], 9, 3), 27)
-1 -1 -1 0 0 0 1 1 1
-1 -1 -1 0 0 0 1 1 1
-1 -1 -1 0 0 0 1 1 1
IDL> print, rebin([-1, 0, 1], 27)
-1 -1 -1 -1 -1 -1 -1 -1 -1
0 0 0 0 0 0 0 0 0
1 1 1 1 1 1 1 1 1

You can use these like:

neighbors = volume[x0 + xoffsets, y0 + yoffsets, z0 + zoffsets]

There is one particularity of IDL that is important in this problem. What happens when x0, y0, or z0 is 0 and you subtract 1 or more from it? Or if x0, y0, or z0 is the last valid value for a dimension and you add 1 to it? By default, IDL will not throw an error, but will use the nearest valid value. For example,

IDL> arr = findgen(5)
IDL> print, arr[[-1, 0, 3, 4, 5]]
0.00000 0.00000 3.00000 4.00000 4.00000

I never knew why this behavior would be desirable, but here it saves the day. By the way, you can turn it off with

compile_opt strictarrsubs

which then will cause an error on the same statement:

IDL> print, arr[[-1, 0, 3, 4, 5]]
% Array used to subscript array contains out of range subscript: ARR.
% Execution halted at: $MAIN$