diff --git a/hippunfold/workflow/rules/autotop.smk b/hippunfold/workflow/rules/autotop.smk index bca28d87..b0465298 100644 --- a/hippunfold/workflow/rules/autotop.smk +++ b/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' - - diff --git a/hippunfold/workflow/rules/compose_warps.smk b/hippunfold/workflow/rules/compose_warps.smk index 16c23862..5f0f0e68 100644 --- a/hippunfold/workflow/rules/compose_warps.smk +++ b/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""" diff --git a/hippunfold/workflow/scripts/compute_graddev.py b/hippunfold/workflow/scripts/compute_graddev.py index 3acc2276..e656aa9e 100644 --- a/hippunfold/workflow/scripts/compute_graddev.py +++ b/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) diff --git a/hippunfold/workflow/scripts/compute_jacobian.py b/hippunfold/workflow/scripts/compute_jacobian.py index 9e3737cd..46f65054 100644 --- a/hippunfold/workflow/scripts/compute_jacobian.py +++ b/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() diff --git a/hippunfold/workflow/scripts/compute_jacobian_simple.py b/hippunfold/workflow/scripts/compute_jacobian_simple.py new file mode 100644 index 00000000..3441fd9f --- /dev/null +++ b/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) + + + diff --git a/hippunfold/workflow/scripts/create_warps.py b/hippunfold/workflow/scripts/create_warps.py index 97a301cd..84c35b07 100644 --- a/hippunfold/workflow/scripts/create_warps.py +++ b/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)