Generating integral points within a spherical shell

Question

Given a spherical shell volume (i.e. one sphere subtracted by a smaller sphere of same center) I need to calculate all integral positions within this volume.

I can do this easily by calculating all integral points within the outer radius (by a cartesian product of all integers up to the outer diameter) and subsequently removing the points within the inner radius.

However, in some situations the radii can be rather large, such that generating all the points within the outer radius can be more time consuming than I would like. So I would like to eliminate the processing time wasted calculating the points within the inner radius (or at least minimize the number necessary), if possible.

Here is some current code, however the radii can be larger.

import numpy as n
import numpy.linalg as la

def cartesian_product(*arrays):
    la = len(arrays)
    dtype = n.result_type(*arrays)
    arr = n.empty([len(a) for a in arrays] + [la], dtype=dtype)
    for i, a in enumerate(n.ix_(*arrays)):
        arr[..., i] = a
    return arr.reshape(-1, la)

def points_in_shell(inner_radius, outer_radius):
    cutoff = int(n.ceil(outer_radius))
    int_range = n.arange(-cutoff, cutoff+1)
    pos_initial = cartesian_product(int_range, int_range, int_range)
    pos_norm = la.norm(pos_initial, axis=1)

    valid_idx = n.where((pos_norm >= inner_radius) & (pos_norm <= outer_radius))
    pos_final = pos_initial[valid_idx]
    return pos_final

vals = points_in_shell(89.634, 95.254)

Any thoughts?

A simple list comprehension might be better, if your shell is thin. — llllllllll

Paul Panzer Paul Panzer · Accepted Answer · 2018-03-05T00:58:20

Update: Adapted to updated Question and also made solution fully sparse.

def f_sparser_oct(IR, OR):
    OR2, IR2 = OR*OR, IR*IR
    out = np.ones((int(OR)+1, 3), dtype=int)
    out[:, 0] = i = np.arange(int(OR)+1)
    i2 = i*i
    jh = np.sqrt(OR2-i2).astype(int)
    out = np.repeat(out, jh+1, axis=0)
    I, J, K = out.T
    J[0] = 0
    J[np.cumsum(jh[:-1]+1)] = -jh[:-1]
    J[:] = np.cumsum(J)
    kh = np.sqrt(OR2-i2[I]-i2[J]).astype(int)
    kl = np.ceil(np.sqrt(np.maximum(0, IR2-i2[I]-i2[J]))).astype(int)
    out = np.repeat(out, 1+kh-kl, axis=0)
    I, J, K = out.T
    K[0] = kl[0]
    K[np.cumsum(1 + kh[:-1] - kl[:-1])] = kl[1:] - kh[:-1]
    K[:] = np.cumsum(K)
    return out

def f_sparser_full(IR, OR):
    OR2, IR2 = OR*OR, IR*IR
    out = np.ones((2*int(OR)+1, 3), dtype=int)
    out[:, 0] = i = np.arange(-int(OR), int(OR)+1)
    i2 = i*i
    i2i = np.r_[i2[int(OR):], i2[:int(OR)]]
    jh = np.sqrt(OR2-i2).astype(int)
    out = np.repeat(out, 2*jh+1, axis=0)
    I, J, K = out.T
    J[0] = -jh[0]
    J[np.cumsum(2*jh[:-1]+1)] = -jh[1:] - jh[:-1]
    J[:] = np.cumsum(J)
    kh = np.sqrt(OR2-i2i[I]-i2i[J]).astype(int)
    kl = np.ceil(np.sqrt(np.maximum(0, IR2-i2i[I]-i2i[J]))).astype(int)
    szs = 2*(1 + kh - kl)-(kl==0)
    out = np.repeat(out, szs, axis=0)
    I, J, K = out.T
    K[0] = -kh[0]
    SZS = np.cumsum(szs)
    K[SZS[:-1]] = -kh[1:] - kh[:-1]
    K[SZS[kl!=0] - 1 - kh[kl!=0] + kl[kl!=0]] = 2 * kl[kl!=0]
    K[:] = np.cumsum(K)
    return out

Demo:

>>> np.all(f_sparser_full(89.634, 95.254) == vals)
True
>>> t0 = time.perf_counter(); f_sparser_full(998, 1000).shape; t1 = time.perf_counter()
(25092914, 3)
>>> t1-t0
0.969647804973647

Update ends.

Here are two quick numpy suggestions: The first is very short but uses a lot of memory. The second has a loop but uses less memory:

import numpy as np

def f_where(IR, OR):
    i, j, k = np.ogrid[:OR+1, :OR+1, :OR+1]
    R2 = i*i + j*j + k*k
    return np.argwhere((R2>=IR*IR) & (R2<=OR*OR))

def f_sparse(IR, OR):
    out = []
    for I in range(int(OR)+1):
        O = np.sqrt(OR*OR - I*I)
        i, j, k = np.ogrid[I:I+1, :O+1, :O+1]
        R2 = i*i + j*j + k*k
        out.append(np.argwhere((R2>=IR*IR) & (R2<=OR*OR))+(I, 0, 0))
    return np.concatenate(out, axis=0)

Demo:

>>> np.all(f_sparse(89.634, 95.254) == vals)
True
>>> 
>>> import time
>>> 
>>> t0 = time.perf_counter(); f_where(298.5, 300.0).shape; t1 = time.perf_counter()
(211366, 3)
>>> t1-t0
0.20139808906242251
>>> t0 = time.perf_counter(); f_sparse(298.5, 300.0).shape; t1 = time.perf_counter()
(211366, 3)
>>> t1-t0
0.1110449800034985
>>> t0 = time.perf_counter(); f_sparse(888.513, 900.099).shape; t1 = time.perf_counter()
(14577694, 3)
>>> t1-t0
3.6043728749500588

Generating integral points within a spherical shell

2 Answers