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)