committing to save progress on a WIP

hippunfold/workflow/rules/autotop.smk

@@ -127,56 +127,4 @@ rule map_to_full_grid:
         '{params.script} {input.coords_ap} {input.coords_pd} {input.coords_io} {input.unfold_ref} {params.warp_dir} &> {log}'
 
 
-<<<<<<< HEAD
-rule compose_warps_corobl2unfold_rhemi:
-    """ Compose corobl to unfold (unfold-template), for right hemi (ie no flip)"""
-    input:
-        native2unfold = bids(root='work',**config['subj_wildcards'],suffix='autotop/WarpITK_native2unfold.nii',desc='cropped',space='corobl',hemi='R',modality='{modality}'),
-        unfold2unfoldtemplate = bids(root='work',**config['subj_wildcards'],suffix='autotop/WarpITK_unfold2unfoldtemplate_0Warp.nii.gz',desc='cropped',space='corobl',hemi='R',modality='{modality}'),
-        ref = bids(root='work',space='unfold',suffix='refvol.nii.gz',**config['subj_wildcards'])
-    output:
-        corobl2unfold = bids(root='work',datatype='seg_{modality}',**config['subj_wildcards'],suffix='xfm.nii.gz',hemi='R',from_='corobl',to='unfold',mode='image')
-    container: config['singularity']['autotop'] 
-    group: 'subj'
-    shell: 'ComposeMultiTransform 3 {output.corobl2unfold} -R {input.ref} {input.unfold2unfoldtemplate} {input.native2unfold}'
-
-
-rule compose_warps_corobl2unfold_lhemi:
-    """ Compose corobl to unfold (unfold-template), for left hemi (ie with flip)"""
-    input:
-        native2unfold = bids(root='work',**config['subj_wildcards'],suffix='autotop/WarpITK_native2unfold.nii',desc='cropped',space='corobl',hemi='Lflip',modality='{modality}'),
-        unfold2unfoldtemplate = bids(root='work',**config['subj_wildcards'],suffix='autotop/WarpITK_unfold2unfoldtemplate_0Warp.nii.gz',desc='cropped',space='corobl',hemi='Lflip',modality='{modality}'),
-        ref = bids(root='work',space='unfold',suffix='refvol.nii.gz',**config['subj_wildcards']),
-        flipLR_xfm = os.path.join(config['snakemake_dir'],'resources','desc-flipLR_type-itk_xfm.txt')
-    output:
-        corobl2unfold = bids(root='work',datatype='seg_{modality}',**config['subj_wildcards'],suffix='xfm.nii.gz',hemi='L',from_='corobl',to='unfold',mode='image')
-    container: config['singularity']['autotop'] 
-    group: 'subj'
-    shell: 'ComposeMultiTransform 3 {output.corobl2unfold} -R {input.ref} {input.unfold2unfoldtemplate} {input.native2unfold} {input.flipLR_xfm}'
-
-
-
-rule compose_warps_t1_to_unfold:
-    """ Compose warps from T1w to unfold """
-    input:
-        corobl2unfold = bids(root='work',datatype='seg_{modality}',**config['subj_wildcards'],suffix='xfm.nii.gz',hemi='{hemi}',from_='corobl',to='unfold',mode='image'),
-        ref = bids(root='work',space='unfold',suffix='refvol.nii.gz',**config['subj_wildcards']),
-        t1w2corobl = bids(root='work',datatype='anat',**config['subj_wildcards'],suffix='xfm.txt',from_='T1w',to='corobl',desc='affine',type_='itk'),
-    output: 
-        bids(root='results',datatype='seg_{modality}',**config['subj_wildcards'],suffix='xfm.nii.gz',hemi='{hemi}',from_='T1w',to='unfold',mode='image')
-    container: config['singularity']['autotop'] 
-    group: 'subj'
-    shell: 'ComposeMultiTransform 3 {output} -R {input.ref} {input.corobl2unfold} {input.t1w2corobl}'
-
-
-
-rule compute_jacobian:
-    input:
-        bids(root='results',datatype='seg_{modality}',**config['subj_wildcards'],suffix='xfm.nii.gz',hemi='{hemi}',from_='T1w',to='unfold',mode='image')
-    output:
-        bids(root='results',datatype='seg_{modality}',**config['subj_wildcards'],suffix='jacobian.nii.gz',hemi='{hemi}',from_='T1w',to='unfold',mode='image')
-    group: 'subj'
-    script: '../scripts/compute_jacobian.py'
-
-
 

hippunfold/workflow/rules/compose_warps.smk

@@ -31,25 +31,25 @@ rule compose_warps_corobl2unfold_lhemi:
 rule compose_warps_unfold2corobl_rhemi:
     """ Compose unfold (unfold-template) to corobl, for right hemi (ie no flip)"""
     input:
-        unfold2native = bids(root='work',**config['subj_wildcards'],suffix='autotop/WarpITK_unfold2native_RPI.nii',desc='cropped',space='corobl',hemi='R',modality='{modality}'),
+        unfold2native = bids(root='work',**config['subj_wildcards'],suffix='autotop/WarpITK_unfold2native.nii',desc='cropped',space='corobl',hemi='R',modality='{modality}'),
         unfoldtemplate2unfold = bids(root='work',**config['subj_wildcards'],suffix='autotop/WarpITK_unfold2unfoldtemplate_0InverseWarp.nii.gz',desc='cropped',space='corobl',hemi='R',modality='{modality}'),
         unfold_ref = bids(root='work',space='unfold',suffix='refvol.nii.gz',**config['subj_wildcards']),
-        ref = bids(root='work',**config['subj_wildcards'],suffix='autotop/coords-AP.nii.gz',desc='cropped',space='corobl',hemi='R',modality='{modality}')
+        ref = bids(root='work',datatype='seg_{modality}',dir='AP',suffix='coords.nii.gz',space='corobl',hemi='R', **config['subj_wildcards']),
     output:
         unfold2corobl = bids(root='results',datatype='seg_{modality}',**config['subj_wildcards'],suffix='xfm.nii.gz',hemi='R',from_='unfold',to='corobl',mode='image')
     container: config['singularity']['ants']
     group: 'subj'
-    shell: 'antsApplyTransforms -o [{output.unfold2corobl},1] -r {input.unfold2native} -t {input.unfold2native} -t {input.unfoldtemplate2unfold} -i {input.unfold_ref} -v'
+    shell: 'antsApplyTransforms -o [{output.unfold2corobl},1] -r {input.ref} -t {input.unfold2native} -t {input.unfoldtemplate2unfold} -i {input.unfold_ref} -v'
 
 
 rule compose_warps_unfold2corobl_lhemi:
     """ Compose unfold (unfold-template) to corobl, for left hemi (ie with flip)"""
     input:
-        unfold2native = bids(root='work',**config['subj_wildcards'],suffix='autotop/WarpITK_unfold2native_RPI.nii',desc='cropped',space='corobl',hemi='Lflip',modality='{modality}'),
+        unfold2native = bids(root='work',**config['subj_wildcards'],suffix='autotop/WarpITK_unfold2native.nii',desc='cropped',space='corobl',hemi='Lflip',modality='{modality}'),
         unfoldtemplate2unfold = bids(root='work',**config['subj_wildcards'],suffix='autotop/WarpITK_unfold2unfoldtemplate_0InverseWarp.nii.gz',desc='cropped',space='corobl',hemi='Lflip',modality='{modality}'),
         flipLR_xfm = os.path.join(config['snakemake_dir'],'resources','desc-flipLR_type-itk_xfm.txt'),
         unfold_ref = bids(root='work',space='unfold',suffix='refvol.nii.gz',**config['subj_wildcards']),
-        ref = bids(root='work',**config['subj_wildcards'],suffix='autotop/coords-AP.nii.gz',desc='cropped',space='corobl',hemi='L',modality='{modality}') #use unflipped ref here
+        ref = bids(root='work',datatype='seg_{modality}',dir='AP',suffix='coords.nii.gz',space='corobl',hemi='L', **config['subj_wildcards']),
     output:
         unfold2corobl = bids(root='results',datatype='seg_{modality}',**config['subj_wildcards'],suffix='xfm.nii.gz',hemi='L',from_='unfold',to='corobl',mode='image')
     container: config['singularity']['ants']
@@ -92,7 +92,7 @@ rule compute_warp_jacobian:
     input: '{prefix}_xfm.nii.gz'
     output: '{prefix}_jacobian.nii.gz'
     group: 'subj'
-    script: '../scripts/compute_jacobian.py'
+    script: '../scripts/compute_jacobian_simple.py'
 
 rule compute_graddev:
     """ Generic rule for computing graddev from jacobian"""

hippunfold/workflow/scripts/compute_graddev.py

@@ -11,23 +11,47 @@
 
 #shape of Jacobian is (256, 128, 16, 9)
 
+print(jac_data.shape)
+
 #Reshape to 3x3 for decomposition
 grads = jac_data.reshape(-1,3,3,order='F')
 
+print(grads.shape)
+
+test_indices = [1000,2000]
+
+
 #graddev = grads[:,:,:]/jac_nib.affine[0][0]
 graddev = grads
 #graddev = jac_data.reshape(-1,3,3,order='F')
 
+
+
 temp_graddev=np.copy(graddev)
 graddev[:]=np.NaN
 
+
+
 #Polar decomposition of jacobian where we only take rotation
-for g in range(0,graddev.shape[0]):
+for g in test_indices:  #range(0,graddev.shape[0]):
+    print(f'index: {g}')    
     m=temp_graddev[g,:,:]
+    print(f'm: {m}')
     x = ~np.isnan(m)
+    print(f'x: {x}')
     if x.all(): 
         rot, deform = polar(m)
+        print(f'rot:')
+        print(f'{rot}')
+        print(f'deform:')
+        print(f'{deform}')
         graddev[g,:,:]= rot
+        print('eye:')
+        print(f'{np.identity(3)}')
+        print('rot - eye:')
+        sub_eye = rot-np.identity(3)
+        print(f'{sub_eye}')
+
 
 grad_hold = graddev - np.identity(3)
 

hippunfold/workflow/scripts/compute_jacobian.py

@@ -35,7 +35,7 @@ def nanChecker(A,B,E=None):
 
 #find which mask to use based on warp
 if 'to-corobl' in in_warp:
-    mask = f'work/{subnum}/{subnum}_hemi-{hemi}_space-corobl_desc-cropped_modality-T2w_autotop/coords-AP.nii.gz'
+    mask = f'work/{subnum}/seg_T2w/{subnum}_dir-AP_hemi-{hemi}_space-corobl_coords.nii.gz'
     masker = nib.load(mask)
     mask = masker.get_fdata()
 

hippunfold/workflow/scripts/compute_jacobian_simple.py

@@ -0,0 +1,59 @@
+import os
+import numpy as np
+import nibabel as nib
+
+in_warp = snakemake.input[0]
+out_jac = snakemake.output[0]
+
+warp_nib = nib.load(in_warp)
+warp_data = warp_nib.get_fdata()
+print(warp_data.shape)
+
+#shape is (256, 128, 16, 1, 3)
+
+#remove the singleton dimension to make easier to work with
+warp_data = np.squeeze(warp_data)
+
+# for u, v, w compute gradient
+gx_u, gx_v, gx_w = np.gradient(warp_data[:,:,:,0])
+gy_u, gy_v, gy_w = np.gradient(warp_data[:,:,:,1])
+gz_u, gz_v, gz_w = np.gradient(warp_data[:,:,:,2])
+
+print(warp_data[62,101,59,0])
+print(warp_data[62,101,59,1])
+print(warp_data[62,101,59,2])
+
+print(gx_u[62,101,59])
+print(gx_v[62,101,59])
+print(gx_w[62,101,59])
+print(gy_u[62,101,59])
+print(gy_v[62,101,59])
+print(gy_w[62,101,59])
+print(gz_u[62,101,59])
+print(gz_v[62,101,59])
+print(gz_w[62,101,59])
+
+
+
+
+#make shape vec for output
+shape_jac = np.array(warp_data.shape)
+shape_jac[3] = 9
+
+jac_data = np.zeros(shape_jac)
+jac_data[:,:,:,0] = gx_u
+jac_data[:,:,:,1] = gy_u
+jac_data[:,:,:,2] = gz_u
+jac_data[:,:,:,3] = gx_v
+jac_data[:,:,:,4] = gy_v
+jac_data[:,:,:,5] = gz_v
+jac_data[:,:,:,6] = gx_w
+jac_data[:,:,:,7] = gy_w
+jac_data[:,:,:,8] = gz_w
+
+jac_nib = nib.Nifti1Image(jac_data, warp_nib.affine)
+
+nib.save(jac_nib,out_jac)
+
+
+

hippunfold/workflow/scripts/create_warps.py

@@ -225,14 +225,14 @@ def summary(name, array):
 
 #now, can write it to file
 warp_native2unfold_nib = nib.Nifti1Image(native_to_unfold,
-                                        coord_ap_nib.affine,
-                                        coord_ap_nib.header)
+                                        coord_ap_nib.affine) #,
+#                                        coord_ap_nib.header)
 warp_native2unfold_nib.to_filename(snakemake.output.warp_native2unfold)
 
 #and save ITK warp too
 warpitk_unfold2native_nib = nib.Nifti1Image(convert_warp_to_itk(native_to_unfold),
-                                        coord_ap_nib.affine,
-                                        coord_ap_nib.header)
+                                        coord_ap_nib.affine) #,
+#                                        coord_ap_nib.header)
 warpitk_unfold2native_nib.to_filename(snakemake.output.warpitk_unfold2native)