Added test case to debug transformation to unfolded space

hippunfold/workflow/grad_test/apply_graddev.py

@@ -0,0 +1,32 @@
+import nibabel as nib
+import numpy as np
+
+#Load in vector
+img = nib.load(snakemake.input[0])
+imgdat = img.get_fdata()
+affine = img.affine
+
+#Load in grad_dev
+grad = nib.load(snakemake.input[1])
+graddev = grad.get_fdata()
+
+nodunf = np.zeros((imgdat.shape))
+
+#Reshape grad_dev to 3x3 using column style
+newgrad = graddev.reshape(graddev.shape[0],graddev.shape[1],graddev.shape[2],3,3,order='F')
+
+#Add back identity matrix
+grad_new = newgrad + np.identity(3)
+
+#At each point in unfolded space, multiply the vector by the grad_dev
+for ii in range(imgdat.shape[0]):
+    for jj in range(imgdat.shape[1]):
+        for kk in range(imgdat.shape[2]):
+            nodunf[ii,jj,kk,:] = grad_new[ii,jj,kk,:,:] @ imgdat[ii,jj,kk,:]      
+
+
+
+newer = nib.Nifti1Image(nodunf, affine)
+nib.save(newer,snakemake.output[0])
+
+

hippunfold/workflow/grad_test/rotate_vec.py

@@ -0,0 +1,47 @@
+import pandas as pd
+import numpy as np
+import nibabel as nib
+
+#Designed to take in an ITK affine transform and apply it to a bvector file. Will save resulting transformed bvec to savingpath
+
+transformfile = snakemake.input[1]
+data = pd.read_csv(transformfile)
+#This part is very specific to the way ITK provides transforms
+x = data.iloc[2][0][12:]
+x = x.split()
+arr = np.array(x)
+
+trans = arr[0:9]
+translation = arr[9:]
+
+transform = np.zeros(9)
+
+for i in range(len(trans)):
+    transform[i] = float(trans[i])  
+
+transform = transform.reshape((3,3))
+
+vecfile = snakemake.input[0]
+vec = nib.load(vecfile)
+nodvec = vec.get_fdata()
+affine = vec.affine
+
+#coords = snakemake.input[2]
+#coord = nib.load(coords)
+#affine = coord.affine
+
+
+nodvecCO = np.zeros((nodvec.shape))
+
+for ii in range(nodvec.shape[0]):
+    for jj in range(nodvec.shape[1]):
+        for kk in range(nodvec.shape[2]):
+            nodvecCO[ii,jj,kk,:] = transform @ nodvec[ii,jj,kk,:]      
+
+
+
+
+savingpath = snakemake.output[0]
+
+newer = nib.Nifti1Image(nodvecCO, affine)
+nib.save(newer,savingpath)

hippunfold/workflow/grad_test/snakefile

@@ -0,0 +1,111 @@
+from os.path import join
+
+in_scratchdir = '/scratch/bkarat/hippdev/hippunfold/hippunfold/workflow/grad_test'
+
+
+hemi = ['L','R']
+
+#Plan->Get vec to CO either through ants here or apply manually and then transform those manuals to corobl. Then get those corobl vectors to unfold using XFM, then apply grad dev
+
+
+#Takes in gradient vector and manually applies rotation from T1 to CO orientation
+rule resample_crop_L_to_corobl:
+    input: 
+        crop0 = in_scratchdir + '/APgradienttest_norm.nii.gz',
+	transform = in_scratchdir + '/sub-111312_desc-affine_from-T1w_to-corobl_type-itk_xfm.txt'
+    output:
+        cor = in_scratchdir + '/grad_IC_COtrans_cropT1space.nii.gz'
+    resources:
+    	mem_mb= 1000
+    script: 'rotate_vec.py'
+
+#Splits the T1 space, CO oriented vector
+rule split_nod:
+    input: 
+        nod = in_scratchdir + '/grad_IC_COtrans_cropT1space.nii.gz'
+    params: 
+        out_prefix = in_scratchdir + '/split_'
+    output:
+        nodsplit = expand(in_scratchdir + '/split_{idx:04d}.nii.gz',idx=range(3))
+    resources:
+    	mem_mb= 1000
+    shell: 'fslsplit {input.nod} {params.out_prefix} -t'
+
+#Makes the vector 5D for ants
+rule push_5D_vec:
+    input: 
+        split = expand(in_scratchdir + '/split_{idx:04d}.nii.gz',idx=range(3))
+    output:
+        push = in_scratchdir + '/grad_4Dto5D_IC.vec0.nii.gz'
+    resources:
+    	mem_mb= 800
+    shell: 'c3d {input.split} -omc {output.push}'
+
+#Take CO oriented vectors and now transform the space from T1 to CO
+rule nod_to_corobl_L:
+    input: 
+        nod = in_scratchdir + '/grad_4Dto5D_IC.vec0.nii.gz',
+	ref = in_scratchdir + '/coords-AP.nii.gz',
+	trans = in_scratchdir + '/sub-111312_desc-affine_from-T1w_to-corobl_type-itk_xfm.txt'
+    output:
+        nod = in_scratchdir + '/grad_IC_5D_crop.vec0.nii.gz'
+    resources:
+    	mem_mb= 2000
+    shell: 'antsApplyTransforms -d 3 -e 1 -n NearestNeighbor -i {input.nod} -r {input.ref} -t {input.trans} -o {output.nod}'
+
+#Split again
+rule pull_4D_vec:
+    input: 
+        split = in_scratchdir + '/grad_IC_5D_crop.vec0.nii.gz'
+    output:
+        push = expand(in_scratchdir + '/c3dsplit{idx:04d}.nii.gz',idx=range(3))
+    resources:
+    	mem_mb= 800
+    shell: 'c3d -mcs {input.split} -oo {output.push}'
+
+#Back to 4D vector. Inspecting this output looks correct, the vectors are oriented AP in CO space
+rule merge_L:
+    input: 
+        crop0 = expand(in_scratchdir + '/c3dsplit{idx:04d}.nii.gz',idx=range(3))
+    output:
+        push = in_scratchdir + '/grad_IC_corobl.vec0.nii.gz'
+    resources:
+    	mem_mb= 800
+    shell: 'fslmerge -t {output.push} {input.crop0}'
+
+#Takes in each X,Y,Z component of the CO vector and moves it to unfolded space (with the CO orientation still)
+rule resample_CO_to_unfold_L:
+    input: 
+        nod = in_scratchdir + '/c3dsplit{idx}.nii.gz',
+	ref = in_scratchdir + '/sub-111312_hemi-L_space-unfold_desc-subfields_dseg.nii.gz',
+	trans = in_scratchdir + '/sub-111312_hemi-L_from-corobl_to-unfold_mode-image_xfm.nii.gz'
+    output:
+        nod = in_scratchdir + '/grad_IC_COtrans_unfspace_{idx}.nii.gz'
+    resources:
+    	mem_mb= 1000
+    shell: 'antsApplyTransforms -d 3 -i {input.nod} -r {input.ref} -t {input.trans}  -o {output.nod}'
+
+#Merge CO orientation unfolded space components
+rule merge_unf_L:
+    input: 
+        crop0 = expand(in_scratchdir + '/grad_IC_COtrans_unfspace_{idx:04d}.nii.gz',idx=range(3))
+    output:
+        push = in_scratchdir + '/grad_IC_COtrans_unfspace.nii.gz'
+    resources:
+    	mem_mb= 700
+    shell: 'fslmerge -t {output.push} {input.crop0}'
+
+
+#Finally, manually rotate vectors to unfolded orientation using graddev
+rule rotate_corobl_to_unf_L:
+    input: 
+        crop0 = in_scratchdir + '/grad_IC_COtrans_unfspace.nii.gz',
+	grad = in_scratchdir + '/sub-111312_hemi-L_from-corobl_to-unfold_mode-image_graddev.nii.gz'
+    output:
+        cor = in_scratchdir + '/grad_IC_unfolded.vec0.nii.gz'
+    resources:
+    	mem_mb= 1000
+    script: 'apply_graddev.py'
+
+
+

hippunfold/workflow/grad_test/sub-111312_desc-affine_from-T1w_to-corobl_type-itk_xfm.txt

@@ -0,0 +1,5 @@
+#Insight Transform File V1.0
+#Transform 0
+Transform: MatrixOffsetTransformBase_double_3_3
+Parameters: 0.945101 -0.0218978960145 0.00041464939890000115 0.0141472 0.7493499974345 -0.5679727159843 0.0218898 0.5800489902646999 0.7564422383045 0.2739182978189994 3.7043449122720005 -9.401099845824998
+FixedParameters: 0 0 0

hippunfold/workflow/scripts/compute_graddev.py

@@ -14,8 +14,8 @@
 #Reshape to 3x3 for decomposition
 grads = jac_data.reshape(-1,3,3,order='F')
 
-graddev = grads[:,:,:]/jac_nib.affine[0][0]
-
+#graddev = grads[:,:,:]/jac_nib.affine[0][0]
+graddev = grads
 #graddev = jac_data.reshape(-1,3,3,order='F')
 
 temp_graddev=np.copy(graddev)