No more overlapping points in griddata interpolation in create_warps.py

hippunfold/workflow/rules/autotop.smk

@@ -88,6 +88,7 @@ rule create_warps:
         warpitk_native2unfold= bids(root='work',**config['subj_wildcards'],suffix='autotop/WarpITK_native2unfold.nii',desc='cropped',space='corobl',hemi='{hemi,Lflip|R}',modality='{modality}'),
         #warp_unfold2native_extrap = bids(root='work',**config['subj_wildcards'],suffix='autotop/Warp_unfold2native_extrapolateNearest.nii',desc='cropped',space='corobl',hemi='{hemi,Lflip|R}',modality='{modality}'),
     group: 'subj'
+    log: bids(root='logs',**config['subj_wildcards'],hemi='{hemi,Lflip|R}',modality='seg-{modality}',suffix='create_warps.txt')
     script: '../scripts/create_warps.py'
 
 # TODO: add this extrapolation of the warp file, extrapolate_warp_unfold2native.m
@@ -104,6 +105,7 @@ rule map_to_full_grid:
         warpitk_native2unfold= bids(root='work',**config['subj_wildcards'],suffix='autotop/WarpITK_native2unfold.nii',desc='cropped',space='corobl',hemi='{hemi}',modality='{modality}'),
         unfold_ref = bids(root='work',space='unfold',suffix='refvol.nii.gz',**config['subj_wildcards'])
     params:
+        dims = config['unfold_vol_ref']['dims'],
         script = os.path.join(config['snakemake_dir'],'workflow','scripts','mapUnfoldToFullGrid.sh'),
         warp_dir = bids(root='work',**config['subj_wildcards'],suffix='autotop',desc='cropped',space='corobl',hemi='{hemi,Lflip|R}',modality='{modality}')
     output:
@@ -119,7 +121,7 @@ rule map_to_full_grid:
     log: bids(root='logs',**config['subj_wildcards'],space='corobl',hemi='{hemi,Lflip|R}',modality='seg{modality}',suffix='mapUnfoldToFullGrid.txt')
     shell:
         'SINGULARITYENV_ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS={threads} '
-        '{params.script} {input.coords_ap} {input.coords_pd} {input.coords_io} {input.unfold_ref} {params.warp_dir} &> {log}'
+        '{params.script} {input.coords_ap} {input.coords_pd} {input.coords_io} {input.unfold_ref} {params.warp_dir} {params.dims} &> {log}'
 
 
 rule compose_warps_corobl2unfold_rhemi:

hippunfold/workflow/scripts/create_warps.py

@@ -2,19 +2,19 @@
 import numpy as np
 from  scipy.interpolate import griddata
 
+logfile = open(snakemake.log[0], 'w')
+print(f'Start', file=logfile, flush=True)
+
 def convert_warp_to_itk(warp):
     """ Convert warp to ITK convention by negating X and Y"""
     warp_itk = warp.copy()
     warp_itk[:,:,:,0,0] = -warp_itk[:,:,:,0,0]
     warp_itk[:,:,:,0,1] = -warp_itk[:,:,:,0,1]
     return warp_itk
 
-
-
-
 def summary(name, array):
     """simple function to print stats on an array"""
-    print(f'{name}: shape={array.shape} mean={array.mean()} max={array.max()} min={array.min()}')
+    print(f'{name}: shape={array.shape} mean={np.nanmean(array)} max={np.nanmax(array)} min={np.nanmin(array)}, numNaNs={np.count_nonzero(np.isnan(array))}, type={array.dtype.name}', file=logfile, flush=True)
     return
 
 #params:
@@ -34,15 +34,13 @@ def summary(name, array):
 coord_io = coord_io_nib.get_fdata()
 
 #get mask of coords  (note: this leaves out coord=0)
-mask = (coord_ap > 0) & (coord_pd > 0) & (coord_io > 0) # matlab: mask = (coord_ap>0 & coord_pd>0 & coord_io>0);
+mask = (coord_ap > 0) | (coord_pd > 0) # some points were lost, especially IO
 num_mask_voxels = np.sum(mask>0)
-print(f'num_mask_voxels {num_mask_voxels}')
+print(f'num_mask_voxels {num_mask_voxels}', file=logfile, flush=True)
 
 #get indices of mask voxels
 idxgm = np.flatnonzero(mask) #matlab: idxgm = find(mask ==1);
 summary('idxgm',idxgm)
-
-print(f'idxgm shape: {idxgm.shape}')
 sz = mask.shape
 
 # Part 1: unfold2native warps
@@ -62,30 +60,25 @@ def summary(name, array):
 
 #unravel indices of mask voxels into subscripts...
 (i_L,j_L,k_L) = np.unravel_index(idxgm,sz)  # matlab: [i_L,j_L,k_L]=ind2sub(sz,idxgm);
-
 summary('i_L',i_L)
 
 #... and stack into vectors ...
 native_coords_mat = np.vstack((i_L,j_L,k_L,np.ones(i_L.shape))) #matlab: native_coords_mat = [i_L-1, j_L-1, k_L-1,ones(size(i_L))]';
-
 summary('native_coords_mat',native_coords_mat)
 
 
 #... then,apply native image affine to get world coords ...
-
-print(f'affine: {coord_ap_nib.affine}, affine shape: {coord_ap_nib.affine.shape}')
-
-native_coords_phys = coord_ap_nib.affine @ native_coords_mat
+aff = coord_ap_nib.affine
+print(f'affine: {aff}, affine shape: {aff.shape}', file=logfile, flush=True)
+native_coords_phys = aff @ native_coords_mat
 native_coords_phys = np.transpose(native_coords_phys[:3,:])
-
 summary('native_coords_phys',native_coords_phys)
 
 # get unfolded grid from file:
 unfold_grid_phys = unfold_phys_coords_nib.get_fdata() 
 
 #unfolded space dims
 unfold_dims = unfold_grid_phys.shape[:3]
-
 summary('unfold_grid_phys',unfold_grid_phys)
 
 
@@ -96,7 +89,11 @@ def summary(name, array):
 # we have points defined by coord_flat_{ap,pd,io}, and corresponding value as native_coords_phys[:,i]
 # and we want to interpolate on a grid in the unfolded space
 
-points = (coord_flat_ap,coord_flat_pd,coord_flat_io)
+# add some noise to avoid perfectly overlapping datapoints!
+points = (coord_flat_ap + (np.random.rand(coord_flat_ap.shape[0])-0.5)*1e-6,
+    coord_flat_pd + (np.random.rand(coord_flat_ap.shape[0])-0.5)*1e-6,
+    coord_flat_io + (np.random.rand(coord_flat_ap.shape[0])-0.5)*1e-6)
+
 # get unfolded grid (from 0 to 1, not world coords), using meshgrid:
 #  note: indexing='ij' to swap the ordering of x and y 
 (unfold_gx, unfold_gy, unfold_gz) = np.meshgrid(np.linspace(0,1,unfold_dims[0]),
@@ -116,8 +113,6 @@ def summary(name, array):
                         values=native_coords_phys[:,0],
                         xi=unfold_xi,
                         method=interp_method)
-
-
 summary('interp_ap',interp_ap)
 interp_pd = griddata(points,
                         values=native_coords_phys[:,1],
@@ -128,7 +123,6 @@ def summary(name, array):
                         values=native_coords_phys[:,2],
                         xi=unfold_xi,
                         method=interp_method)
-
 summary('interp_io',interp_ap)
 
 
@@ -139,24 +133,27 @@ def summary(name, array):
 mapToNative[:,:,:,0,0] = interp_ap
 mapToNative[:,:,:,0,1] = interp_pd
 mapToNative[:,:,:,0,2] = interp_io
-
+# TODO: interpolate nans more better
+mapToNative[np.isnan(mapToNative)] = 0
 summary('mapToNative',mapToNative)
 
 # mapToNative has the absolute coordinates, but we want them relative to the 
 # unfolded grid, so we subtract it out:
 displacementToNative = mapToNative - unfold_grid_phys
-
 summary('dispacementToNative',displacementToNative)
 
 
 # write to file
-warp_unfold2native_nib = nib.Nifti1Image(displacementToNative,
+dt = unfold_phys_coords_nib.get_fdata()
+dt = dt.dtype.name
+warp_unfold2native_nib = nib.Nifti1Image(displacementToNative.astype(dt),
                                         unfold_phys_coords_nib.affine,
                                         unfold_phys_coords_nib.header)
 warp_unfold2native_nib.to_filename(snakemake.output.warp_unfold2native)
 
 # write itk transform to file
-warpitk_native2unfold_nib = nib.Nifti1Image(convert_warp_to_itk(displacementToNative),
+f = convert_warp_to_itk(displacementToNative)
+warpitk_native2unfold_nib = nib.Nifti1Image(f.astype(dt),
                                         unfold_phys_coords_nib.affine,
                                         unfold_phys_coords_nib.header)
 
@@ -224,13 +221,16 @@ def summary(name, array):
 summary('native_to_unfold',native_to_unfold)
 
 #now, can write it to file
-warp_native2unfold_nib = nib.Nifti1Image(native_to_unfold,
+dt = coord_ap_nib.get_fdata()
+dt = dt.dtype.name
+warp_native2unfold_nib = nib.Nifti1Image(native_to_unfold.astype(dt),
                                         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),
+f = convert_warp_to_itk(native_to_unfold)
+warpitk_unfold2native_nib = nib.Nifti1Image(f.astype(dt),
                                         coord_ap_nib.affine,
                                         coord_ap_nib.header)
 warpitk_unfold2native_nib.to_filename(snakemake.output.warpitk_unfold2native)

hippunfold/workflow/scripts/mapUnfoldToFullGrid.sh

@@ -11,10 +11,13 @@ coords_PD=$2
 coords_IO=$3
 in_unfold_ref=$4
 warps_dir=$5
+N_AP=$6
+N_PD=$7
+N_IO=$8
 
-if [ "$#" -lt 5 ]
+if [ "$#" -lt 8 ]
 then
-    echo "Usage: $0 coords_ap coords_pd coords_io unfold_ref out_warps_dir"
+    echo "Usage: $0 coords_ap coords_pd coords_io unfold_ref out_warps_dir N_AP N_PD N_IO"
     exit 1
 fi
 
@@ -24,11 +27,6 @@ fill_val=-1
 pad_cmd="-pad 32x32x32vox 32x32x32vox $fill_val"
 
 
-N_AP=256
-N_PD=128
-N_IO=16
-
-
 in_warpitk_native2unfold=${warps_dir}/WarpITK_native2unfold.nii
 out_norm_coords=${warps_dir}/unfold_norm_coords.nii