diff --git a/hippunfold/workflow/rules/warps.smk b/hippunfold/workflow/rules/warps.smk index cbdd11d8..0e4a40ee 100644 --- a/hippunfold/workflow/rules/warps.smk +++ b/hippunfold/workflow/rules/warps.smk @@ -226,8 +226,6 @@ rule create_warps_hipp: hemi="{hemi}", suffix="create_warps-hipp.txt" ), - container: - config["singularity"]["autotop"] script: "../scripts/create_warps.py" @@ -294,7 +292,6 @@ rule create_warps_dentate: ), labelmap=get_labels_for_laplace, params: - interp_method="linear", gm_labels=lambda wildcards: config["laplace_labels"]["PD"]["sink"], epsilon=lambda wildcards: config["unfold_crop_epsilon_fractions"], resources: @@ -353,8 +350,6 @@ rule create_warps_dentate: hemi="{hemi}", suffix="create_warps-dentate.txt" ), - container: - config["singularity"]["autotop"] script: "../scripts/create_warps.py" diff --git a/hippunfold/workflow/scripts/create_warps.py b/hippunfold/workflow/scripts/create_warps.py index 6562d753..3afa245f 100644 --- a/hippunfold/workflow/scripts/create_warps.py +++ b/hippunfold/workflow/scripts/create_warps.py @@ -1,7 +1,7 @@ import nibabel as nib import numpy as np -from scipy.interpolate import griddata -from scipy.interpolate import NearestNDInterpolator +import naturalneighbor +from scipy.stats import zscore logfile = open(snakemake.log[0], "w") print(f"Start", file=logfile, flush=True) @@ -25,9 +25,6 @@ def summary(name, array): return -# params: -interp_method = snakemake.params.interp_method - # load unfolded coordinate map # unfold_ref_nib = nib.load(snakemake.input.unfold_ref_nii) unfold_phys_coords_nib = nib.load(snakemake.input.unfold_phys_coords_nii) @@ -66,7 +63,7 @@ def summary(name, array): # the unfolded space, using scipy's griddata (equivalent to matlab scatteredInterpolant). -coord_flat_ap = coord_ap[mask == True] # matlab: Laplace_AP = coord_ap(mask==1); +coord_flat_ap = coord_ap[mask == True] coord_flat_pd = coord_pd[mask == True] coord_flat_io = coord_io[mask == True] @@ -75,15 +72,11 @@ def summary(name, array): summary("coord_flat_io", coord_flat_io) # 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); +(i_L, j_L, k_L) = np.unravel_index(idxgm, sz) 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))]'; +native_coords_mat = np.vstack((i_L, j_L, k_L, np.ones(i_L.shape))) summary("native_coords_mat", native_coords_mat) @@ -104,20 +97,13 @@ def summary(name, array): # scattered interpolation / griddata: -# matlab: interp_X = scatteredInterpolant(Laplace_AP,Laplace_PD,Laplace_IO,native_coords_phys(:,1),interp,extrap); - # 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 -# add some noise to avoid perfectly overlapping datapoints! -points = ( - coord_flat_ap * unfold_dims[0] - + (np.random.rand(coord_flat_ap.shape[0]) - 0.5) * 1e-6, - coord_flat_pd * unfold_dims[1] - + (np.random.rand(coord_flat_ap.shape[0]) - 0.5) * 1e-6, - coord_flat_io * unfold_dims[2] - + (np.random.rand(coord_flat_ap.shape[0]) - 0.5) * 1e-6, -) +# unnormalize and add a bit of noise so points don't ever perfectly overlap +coord_flat_ap_unnorm = coord_flat_ap * unfold_dims[0] +coord_flat_pd_unnorm = coord_flat_pd * unfold_dims[1] +coord_flat_io_unnorm = coord_flat_io * unfold_dims[2] # get unfolded grid (from 0 to 1, not world coords), using meshgrid: # note: indexing='ij' to swap the ordering of x and y @@ -134,53 +120,40 @@ def summary(name, array): ), indexing="ij", ) - -# tuple for use in griddata: -unfold_xi = (unfold_gx, unfold_gy, unfold_gz) summary("unfold_gx", unfold_gx) summary("unfold_gy", unfold_gy) summary("unfold_gz", unfold_gz) -# perform the interpolation, filling in outside values as 0 -# TODO: linear vs cubic? we were using "natural" interpolation in matlab -# so far, linear seems close enough.. -interp_ap = griddata( - points, 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], xi=unfold_xi, method=interp_method -) -summary("interp_pd", interp_pd) -interp_io = griddata( - points, values=native_coords_phys[:, 2], xi=unfold_xi, method=interp_method +# perform the interpolation + +points = np.stack( + [ + coord_flat_ap * unfold_dims[0] + + (np.random.rand(coord_flat_ap.shape[0]) - 0.5) * 1e-6, + coord_flat_pd * unfold_dims[1] + + (np.random.rand(coord_flat_ap.shape[0]) - 0.5) * 1e-6, + coord_flat_io * unfold_dims[2] + + (np.random.rand(coord_flat_ap.shape[0]) - 0.5) * 1e-6, + ], + axis=1, ) -summary("interp_io", interp_ap) -# fill NaNs (NN interpolater allows for extrapolation!) -[x, y, z] = np.where(np.invert(np.isnan(interp_ap))) -interp = NearestNDInterpolator( - np.c_[x, y, z], interp_ap[np.invert(np.isnan(interp_ap))] +interp_ap = naturalneighbor.griddata( + points, + native_coords_phys[:, 0], + [[0, unfold_dims[0], 1], [0, unfold_dims[1], 1], [0, unfold_dims[2], 1]], ) -[qx, qy, qz] = np.where(np.isnan(interp_ap)) -fill = interp(qx, qy, qz) -interp_ap[np.isnan(interp_ap)] = fill - -[x, y, z] = np.where(np.invert(np.isnan(interp_pd))) -interp = NearestNDInterpolator( - np.c_[x, y, z], interp_pd[np.invert(np.isnan(interp_pd))] +interp_pd = naturalneighbor.griddata( + points, + native_coords_phys[:, 1], + [[0, unfold_dims[0], 1], [0, unfold_dims[1], 1], [0, unfold_dims[2], 1]], ) -[qx, qy, qz] = np.where(np.isnan(interp_pd)) -fill = interp(qx, qy, qz) -interp_pd[np.isnan(interp_pd)] = fill - -[x, y, z] = np.where(np.invert(np.isnan(interp_io))) -interp = NearestNDInterpolator( - np.c_[x, y, z], interp_io[np.invert(np.isnan(interp_io))] +interp_io = naturalneighbor.griddata( + points, + native_coords_phys[:, 2], + [[0, unfold_dims[0], 1], [0, unfold_dims[1], 1], [0, unfold_dims[2], 1]], ) -[qx, qy, qz] = np.where(np.isnan(interp_io)) -fill = interp(qx, qy, qz) -interp_io[np.isnan(interp_io)] = fill + # prepare maps for writing as warp file: @@ -189,10 +162,15 @@ 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) +z = zscore(mapToNative, axis=None) +bad = np.logical_or(z < -5, z > 5) +mapToNative[bad] = np.nan +summary("mapToNative cleaned", mapToNative) +mapToNative[np.isnan(mapToNative)] = 0 + + # mapToNative has the absolute coordinates, but we want them relative to the # unfolded grid, so we subtract it out: displacementToNative = mapToNative - unfold_grid_phys @@ -228,13 +206,6 @@ def summary(name, array): # The image affine from the unfolded grid takes points from 0 to N to world coords, so # just need to un-normalize, then multiply by affine -# unnormalize -coord_flat_ap_unnorm = coord_flat_ap * unfold_dims[0] -coord_flat_pd_unnorm = coord_flat_pd * unfold_dims[1] -coord_flat_io_unnorm = coord_flat_io * unfold_dims[2] - -summary("coord_flat_ap_unnorm", coord_flat_ap_unnorm) - # reshape for multiplication (affine * vec) coord_flat_unnorm_vec = np.stack( @@ -259,7 +230,7 @@ def summary(name, array): summary("uvw_phys", uvw_phys) # now we have the absolute unfold world coords for each native grid point -# but we need the displacement from the native grid point world coord +# but we need the displacement from the native grid point world coord # the native world coords are in native_coords_phys # so we subtract it diff --git a/pyproject.toml b/pyproject.toml index a4c5ae4c..5bbcd235 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -12,6 +12,7 @@ appdirs = "^1.4.4" Jinja2 = "^3.0.3" pygraphviz = "1.7" Pygments = "^2.10.0" +naturalneighbor = "0.2.1" [tool.poetry.dev-dependencies]