diff --git a/hippunfold/workflow/grad_test/APgradienttest_norm.nii.gz b/hippunfold/workflow/grad_test/APgradienttest_norm.nii.gz new file mode 100644 index 00000000..349b688f Binary files /dev/null and b/hippunfold/workflow/grad_test/APgradienttest_norm.nii.gz differ diff --git a/hippunfold/workflow/grad_test/apply_graddev.py b/hippunfold/workflow/grad_test/apply_graddev.py new file mode 100644 index 00000000..613315ff --- /dev/null +++ b/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]) + + diff --git a/hippunfold/workflow/grad_test/c3dsplit0000.nii.gz b/hippunfold/workflow/grad_test/c3dsplit0000.nii.gz new file mode 100644 index 00000000..ae5411c6 Binary files /dev/null and b/hippunfold/workflow/grad_test/c3dsplit0000.nii.gz differ diff --git a/hippunfold/workflow/grad_test/c3dsplit0001.nii.gz b/hippunfold/workflow/grad_test/c3dsplit0001.nii.gz new file mode 100644 index 00000000..2b9dfd74 Binary files /dev/null and b/hippunfold/workflow/grad_test/c3dsplit0001.nii.gz differ diff --git a/hippunfold/workflow/grad_test/c3dsplit0002.nii.gz b/hippunfold/workflow/grad_test/c3dsplit0002.nii.gz new file mode 100644 index 00000000..497b567c Binary files /dev/null and b/hippunfold/workflow/grad_test/c3dsplit0002.nii.gz differ diff --git a/hippunfold/workflow/grad_test/coords-AP.nii.gz b/hippunfold/workflow/grad_test/coords-AP.nii.gz new file mode 100644 index 00000000..b1e23f5c Binary files /dev/null and b/hippunfold/workflow/grad_test/coords-AP.nii.gz differ diff --git a/hippunfold/workflow/grad_test/grad_4Dto5D_IC.vec0.nii.gz b/hippunfold/workflow/grad_test/grad_4Dto5D_IC.vec0.nii.gz new file mode 100644 index 00000000..88890073 Binary files /dev/null and b/hippunfold/workflow/grad_test/grad_4Dto5D_IC.vec0.nii.gz differ diff --git a/hippunfold/workflow/grad_test/grad_IC_5D_crop.vec0.nii.gz b/hippunfold/workflow/grad_test/grad_IC_5D_crop.vec0.nii.gz new file mode 100644 index 00000000..235bcada Binary files /dev/null and b/hippunfold/workflow/grad_test/grad_IC_5D_crop.vec0.nii.gz differ diff --git a/hippunfold/workflow/grad_test/grad_IC_COtrans_cropT1space.nii.gz b/hippunfold/workflow/grad_test/grad_IC_COtrans_cropT1space.nii.gz new file mode 100644 index 00000000..8737ed23 Binary files /dev/null and b/hippunfold/workflow/grad_test/grad_IC_COtrans_cropT1space.nii.gz differ diff --git a/hippunfold/workflow/grad_test/grad_IC_COtrans_unfspace.nii.gz b/hippunfold/workflow/grad_test/grad_IC_COtrans_unfspace.nii.gz new file mode 100644 index 00000000..9b7f1d1b Binary files /dev/null and b/hippunfold/workflow/grad_test/grad_IC_COtrans_unfspace.nii.gz differ diff --git a/hippunfold/workflow/grad_test/grad_IC_COtrans_unfspace_0000.nii.gz b/hippunfold/workflow/grad_test/grad_IC_COtrans_unfspace_0000.nii.gz new file mode 100644 index 00000000..6f126733 Binary files /dev/null and b/hippunfold/workflow/grad_test/grad_IC_COtrans_unfspace_0000.nii.gz differ diff --git a/hippunfold/workflow/grad_test/grad_IC_COtrans_unfspace_0001.nii.gz b/hippunfold/workflow/grad_test/grad_IC_COtrans_unfspace_0001.nii.gz new file mode 100644 index 00000000..f763d15a Binary files /dev/null and b/hippunfold/workflow/grad_test/grad_IC_COtrans_unfspace_0001.nii.gz differ diff --git a/hippunfold/workflow/grad_test/grad_IC_COtrans_unfspace_0002.nii.gz b/hippunfold/workflow/grad_test/grad_IC_COtrans_unfspace_0002.nii.gz new file mode 100644 index 00000000..087e1524 Binary files /dev/null and b/hippunfold/workflow/grad_test/grad_IC_COtrans_unfspace_0002.nii.gz differ diff --git a/hippunfold/workflow/grad_test/grad_IC_corobl.vec0.nii.gz b/hippunfold/workflow/grad_test/grad_IC_corobl.vec0.nii.gz new file mode 100644 index 00000000..c011e253 Binary files /dev/null and b/hippunfold/workflow/grad_test/grad_IC_corobl.vec0.nii.gz differ diff --git a/hippunfold/workflow/grad_test/grad_IC_unfolded.vec0.nii.gz b/hippunfold/workflow/grad_test/grad_IC_unfolded.vec0.nii.gz new file mode 100644 index 00000000..ff8f801f Binary files /dev/null and b/hippunfold/workflow/grad_test/grad_IC_unfolded.vec0.nii.gz differ diff --git a/hippunfold/workflow/grad_test/rotate_vec.py b/hippunfold/workflow/grad_test/rotate_vec.py new file mode 100644 index 00000000..ca21b1e0 --- /dev/null +++ b/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) diff --git a/hippunfold/workflow/grad_test/snakefile b/hippunfold/workflow/grad_test/snakefile new file mode 100644 index 00000000..d3999fa4 --- /dev/null +++ b/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' + + + diff --git a/hippunfold/workflow/grad_test/split_0000.nii.gz b/hippunfold/workflow/grad_test/split_0000.nii.gz new file mode 100644 index 00000000..41f9ff35 Binary files /dev/null and b/hippunfold/workflow/grad_test/split_0000.nii.gz differ diff --git a/hippunfold/workflow/grad_test/split_0001.nii.gz b/hippunfold/workflow/grad_test/split_0001.nii.gz new file mode 100644 index 00000000..ec438a92 Binary files /dev/null and b/hippunfold/workflow/grad_test/split_0001.nii.gz differ diff --git a/hippunfold/workflow/grad_test/split_0002.nii.gz b/hippunfold/workflow/grad_test/split_0002.nii.gz new file mode 100644 index 00000000..b405e1d2 Binary files /dev/null and b/hippunfold/workflow/grad_test/split_0002.nii.gz differ diff --git a/hippunfold/workflow/grad_test/sub-111312_desc-affine_from-T1w_to-corobl_type-itk_xfm.txt b/hippunfold/workflow/grad_test/sub-111312_desc-affine_from-T1w_to-corobl_type-itk_xfm.txt new file mode 100644 index 00000000..7f03a6dd --- /dev/null +++ b/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 diff --git a/hippunfold/workflow/grad_test/sub-111312_hemi-L_from-corobl_to-unfold_mode-image_graddev.nii.gz b/hippunfold/workflow/grad_test/sub-111312_hemi-L_from-corobl_to-unfold_mode-image_graddev.nii.gz new file mode 100644 index 00000000..483b8ba7 Binary files /dev/null and b/hippunfold/workflow/grad_test/sub-111312_hemi-L_from-corobl_to-unfold_mode-image_graddev.nii.gz differ diff --git a/hippunfold/workflow/grad_test/sub-111312_hemi-L_from-corobl_to-unfold_mode-image_xfm.nii.gz b/hippunfold/workflow/grad_test/sub-111312_hemi-L_from-corobl_to-unfold_mode-image_xfm.nii.gz new file mode 100644 index 00000000..679e2c3b Binary files /dev/null and b/hippunfold/workflow/grad_test/sub-111312_hemi-L_from-corobl_to-unfold_mode-image_xfm.nii.gz differ diff --git a/hippunfold/workflow/grad_test/sub-111312_hemi-L_space-unfold_desc-subfields_dseg.nii.gz b/hippunfold/workflow/grad_test/sub-111312_hemi-L_space-unfold_desc-subfields_dseg.nii.gz new file mode 100644 index 00000000..7f92d656 Binary files /dev/null and b/hippunfold/workflow/grad_test/sub-111312_hemi-L_space-unfold_desc-subfields_dseg.nii.gz differ diff --git a/hippunfold/workflow/scripts/compute_graddev.py b/hippunfold/workflow/scripts/compute_graddev.py index 7af496ca..3acc2276 100644 --- a/hippunfold/workflow/scripts/compute_graddev.py +++ b/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)