-
Notifications
You must be signed in to change notification settings - Fork 0
Generating parcels and computing summaries for subsequent analysis
Zaki Alasmar edited this page Jun 7, 2022
·
2 revisions
import numpy
import nighres as nr
import nibabel
m2 = nb.load('mask_file.nii.gz')
num_parcels = 1000 #checked 5000, and 10000 as well
distance_mm = 100.0 #if this is too small, you might miss some voxels
num_vox = d2.sum()
print(num_vox)
selected_vox_idxs = np.round(np.linspace(0,num_vox-1,num=num_parcels)).astype(int)
#select the mask locations based on the index you generated
#give each voxel a unique index (1-based)
vox_locs = np.array(np.where(d2))
initial_idx_mask = np.zeros(d2.shape).astype(int)
idx_val = 1
for idx in selected_vox_idxs:
x = vox_locs[0,idx]
y = vox_locs[1,idx]
z = vox_locs[2,idx]
initial_idx_mask[x,y,z] = idx_val
idx_val += 1
intial_idx_mask_img = nb.Nifti1Image(initial_idx_mask,affine=m2.affine,header=m2.header)
intial_idx_mask_img.set_data_dtype(int)
res = nr.intensity.intensity_propagation(intial_idx_mask_img,domain=domain_mask,combine='min',target='zero',distance_mm=distance_mm)
# this is the number of voxels that are missing from the indexed mask that were originally in the binary mask (this can happen b/c things are not connected)
(res['result'].get_fdata()[m2.get_fdata().astype(bool)]==0).sum()import NeuralABC_tools as nabc #The NeuralABC_tools module should be loaded before conda
import time
mask_f = "res['result']" #my parcellated mask
d_mask_parcel = nb.load(mask_f).get_fdata().astype(int) #note that 0 is ALWAYS background, so be careful about the indexing that you do later on!
parcel_idxs = np.unique(d_mask_parcel) #parcel labels
d_parcel = np.zeros(parcel_idxs.shape[0])
df = someData #df is an array of shape nVoxels X nSubjects (voxels as rows and subjects as columns)
d_vals = np.zeros((d_parcel.shape[0]-1, df.shape[1])) #output array
start = time.time() #just to time it, can be removed if not needed
for sub in range(df.shape[1]):
print(sub, end=" ")
df_sub_img = np.zeros(mask_d.shape)
df_sub_img[mask_d] = df[:,sub]
output = nabc.bin_statistics_by_index(d_mask_parcel,df_sub_img, statistic="median",remove_nans=True) #the options for statistic are = {'mean','std','median','count','sum','min','max',function}, if remove_nans is true, then it outputs nan{'mean,median,std ...'}
d_vals[:, sub] = output["stat"]
end = time.time()
print(f"It took {(end - start) / 60} minutes to finish {df.shape[1]} participants")How to map summarized values (or their processing derivatives/statistics) back into the space from which they came
statistic_d = np.hstack([0,statistic_output]) #Add 0 as a background
statistic_image = nabc.map_vals_to_index(d_mask_parcel, statistic_d)
out_statistic_image = nb.Nifti1Image(dxt_img,nb.load(mask_f).affine,nb.load(mask_f).header)
out_statistic_image.to_filename('NameOfMyFile.nii.gz')