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)