Intersect GM mask with a sphere with specific radius

Hi all.

I have a GM mask in MNI space that includes 150k GM voxels. Link to the mask: https://gofile.io/d/oFKrkL

I want to intersect a sphere with radius r with the GM mask in order to estimate neighbourhood information in the following way. For a specific GM voxel, the neighborhood is defined as the voxels within a radius r based on the voxel coordinates in mm.

I just need to store this information. For example, voxel #0 has neighbors the voxels #10,#20, and #30 etc…

I have used a KDTree to do that but it not really accurate so I am looking for an alternative.

import numpy as np, nibabel as nib

# Load GM mask
gm_mask_header = nib.load("mask.nii") # https://gofile.io/d/oFKrkL
gm_mask = gm_mask_header.get_fdata()                        
gm_mask_bool = np.array(gm_mask, dtype=bool)                
brain_mask = gm_mask_bool.copy()

# Get scanner space coordinates of the voxels in mm
array_indices = np.indices(gm_mask.shape)
vox_indices = array_indices.transpose(1, 2, 3, 0) 
mm = nib.affines.apply_affine(gm_mask_header.affine, vox_indices)
mm = mm.transpose(3, 0, 1, 2)

# Find indices where the GM tissue is located
i_idx,j_idx,k_idx = np.where(brain_mask)

# GM only scanner space coordinates of the voxels in mm
mm_prime = mm[:, i_idx,j_idx,k_idx]

print(mm_prime.shape)
(3, 159459)

Hello! I didn’t have to time to look into details but Nilearn provide a searchlight implementation in python that may give you some ideas :

Thanks for the answer. I will definitely check that