I am looking for an algorithm (best case written in Python) that I can use to ‘grow’ clusters/regions of interest in grey matter space. Is anyone aware of such an algorithm? The idea would be to give the algorithm a set of seed coordinates and a set of region sizes (i.e. number of voxels).
1.) The algorithm should then grow clusters from the seed locations until it reaches the defined size.
2.) It also should respect the space defined by the grey matter mask (that is, don’t let the clusters grow inside white matter space and don’t let them grow outside the whole brain space).
3.) Finally, clusters should be not allowed to overlap
I think you can at least hack the searchlight implementation that is in Nilearn if you’re open to having your region size based on mm distance from seed.
Basically your algorithm should calculate the affinity matrix from your grey matter mask, and for each seed return the voxels in this affinity matrix that are at a distance lower than your threshold.
Hi Thomas,
yeah, my first idea was also to use nilearn.input_data.nifti_spheres_masker._apply_mask_and_get_affinity for that, e.g.:
# one sphere of radius 10mm will be displayed
seeds = [[44,0,26]]
radius = 30
# apply sphere masking on the mask
_, X = _apply_mask_and_get_affinity(seeds=seeds,
niimg=None,
radius=radius,
allow_overlap=True,
mask_img=mask_img
)
# inverse-transform X that contains masked images with spheres
masker = NiftiMasker(mask_img=mask_img)
masker.fit()
img = masker.inverse_transform(np.sum(X,axis=0).A.ravel())
# show results
plot_roi(img)
I guess the radius that has to be set given a number of voxels could be inferred from using:
def _get_sphere_radius_from_volume(V):
return (3 * V /(4 * math.pi)) ** (1. / 3)
However, using such an approach does not guarantee that you end up with a sphere that contains your defined number of voxels. If the seed is placed near white matter or near the brain surface (and the sphere is large enough) it will ‘stop there’ (see image above).
Or did you mean something different?
I guess you could change slightly _apply_mask_and_get_affinity() to use number of nearest neighbor instead of radius. l110 from current nilearn source you could define clf = neighbors.NearestNeighbors( n_neighbors=1000) instead of the radius.
(see sklearn.neighbors.NearestNeighbors — scikit-learn 0.24.1 documentation)
For the overlapping constraint just set allow overlap=False. Then either you’re problem has a solution or the function will raise an error and you should adjust you’re desired number of voxels in each region.
Hi Thomas,
thanks for your help! Based on your suggestions I came up with the following code:
## Get In-Mask coordinates ####################################################
def _get_mask_coords(mask_img):
'''Get in-mask voxel coordinates'''
mask, affine = masking._load_mask_img(mask_img)
# get data indices for all '1' voxels inside mask data array
mask_data_indices = np.asarray(np.nonzero(mask)).T.tolist()
# return coordinates for those '1' voxels
mask_coords = np.asarray(list(zip(*mask_data_indices)))
mask_coords = coord_transform(mask_coords[0],
mask_coords[1],
mask_coords[2],
affine)
mask_coords = np.asarray(mask_coords).T
return mask_coords
mask_coords = _get_mask_coords(mask_img)
## Draw random seed coordinates ###############################################
def _draw_n_random_seeds(coord_array,n_seeds,rng):
'''Draw a defined number of random coordinates from an array from voxel
coordinates'''
coord_array = mask_coords
random_row_indices = rng.randint(coord_array.shape[0],size=n_seeds)
random_seeds = coord_array[random_row_indices].tolist()
return random_seeds
for idx in range(0,10):
overlap = True
while overlap:
# draw random seeds
random_seeds = _draw_n_random_seeds(coord_array=mask_coords,
n_seeds=n_seeds,
rng=rng)
cluster_arrays = []
for random_seed,cluster_size in zip(random_seeds,cluster_sizes):
clf = neighbors.NearestNeighbors(n_neighbors=cluster_size)
nearest_neighbors = clf.fit(mask_coords).kneighbors_graph([random_seed])
cluster_arrays.append(nearest_neighbors)
cluster_arrays_summed = sum(cluster_arrays).A
overlap = np.any((cluster_arrays_summed >= 2))
# inverse-transform cluster arrays
masker = NiftiMasker(mask_img=mask_img)
masker.fit()
random_clusters_mask_img = masker.inverse_transform(cluster_arrays_summed.ravel())
# show results
plot_roi(random_clusters_mask_img)
print(np.count_nonzero(random_clusters_mask_img.get_fdata()))
print(np.unique(random_clusters_mask_img.get_fdata()))
So this creates ten binary images with random clusters of fixed sizes. The ‘overlapping’ problem is handled via the while-loop (the algorithm places random seeds somewhere in grey matter space, grows spheres of fixed sizes using n_neighbours, and checks if they are overlapping. If they do, it tries again). Using n_neighbours now guarantees that I get exactly the number of voxels that I specified for each sphere.
The only problem that is still not solved is that the NearestNeighbour-algorithm does not guarantee that each cluster is a connected component, so it can happen that the next nearest grey-matter voxel is a voxel that is not directly connected to the last voxel (creating ‘disconnected’ clusters).
I guess if I also set this constraint I would have to write my own algorithm where this rule is somehow implemented (“the next nearest voxel must always be connected to the last voxel”). I will leave this thread open for now, maybe I will come up with a solution that I can post here 
Thanks again for your help!
Greetings,
Johannes
Great that you managed !
If you add this constraint (each cluster must be a connected component) then I think you need to check nilearn/regions/rena_clustering.py especially the recursive_neighbor_agglomeration() function.
The matrix “connectivity” is what you need to enforce this local connectedness. Then the aggregation scheme has to be tailored with your constraints.
Good luck
Ha! Interesting! I’ll have a look at that thx!
Hi Thomas,
I have been trying to get an intuition on how I could use ReNA clustering for my purpose, but I am struggling with the code…I can’t find a way to define seeds here? I also don’t really understand the purpose of the connectivity matrix…
Greetings,
Johannes