From 517b9f7c4400b3a0385c55e9fbc7b73a8ebfcf31 Mon Sep 17 00:00:00 2001 From: Jordan DeKraker Date: Wed, 26 Jun 2024 09:59:56 -0400 Subject: [PATCH 1/6] fresh branch for clean changes --- hippunfold/config/snakebids.yml | 7 +- hippunfold/workflow/rules/common.smk | 43 +++++-- hippunfold/workflow/rules/warps.smk | 6 - hippunfold/workflow/scripts/create_warps.py | 118 ++++++++------------ pyproject.toml | 1 + 5 files changed, 90 insertions(+), 85 deletions(-) diff --git a/hippunfold/config/snakebids.yml b/hippunfold/config/snakebids.yml index 66d61bff..d3a3f1bc 100644 --- a/hippunfold/config/snakebids.yml +++ b/hippunfold/config/snakebids.yml @@ -159,6 +159,7 @@ parse_args: --inject_template: choices: - 'upenn' + - 'dHCP' - 'MBMv2' - 'MBMv3' - 'CIVM' @@ -327,11 +328,15 @@ parse_args: - synthseg_v0.2 - neonateT1w_v2 + + + +# --- surface specific configuration -- + autotop_labels: - 'hipp' - 'dentate' - surf_types: hipp: - midthickness diff --git a/hippunfold/workflow/rules/common.smk b/hippunfold/workflow/rules/common.smk index b0e5db03..4c568f95 100644 --- a/hippunfold/workflow/rules/common.smk +++ b/hippunfold/workflow/rules/common.smk @@ -24,41 +24,69 @@ def get_modality_suffix(modality): def get_final_spec(): - if len(config["hemi"]) == 2: - specs = expand( + specs = expand( + bids( + root=root, + datatype="surf", + den="{density}", + space="{space}", + hemi="{hemi}", + label="{autotop}", + suffix="surfaces.spec", + **config["subj_wildcards"], + ), + density=config["output_density"], + space=ref_spaces, + hemi=config["hemi"], + autotop=config["autotop_labels"], + allow_missing=True, + ) + return specs + + +def get_final_surf(): + gii = [] + gii.extend( + expand( bids( root=root, datatype="surf", den="{density}", + suffix="{surfname}.surf.gii", space="{space}", + hemi="{hemi}", label="{autotop}", - suffix="surfaces.spec", **config["subj_wildcards"], ), density=config["output_density"], space=ref_spaces, + hemi=config["hemi"], autotop=config["autotop_labels"], + surfname=config["surf_types"]["hipp"], allow_missing=True, ) - else: - specs = expand( + ) + gii.extend( + expand( bids( root=root, datatype="surf", den="{density}", + suffix="{surfname}.surf.gii", space="{space}", hemi="{hemi}", label="{autotop}", - suffix="surfaces.spec", **config["subj_wildcards"], ), density=config["output_density"], space=ref_spaces, hemi=config["hemi"], autotop=config["autotop_labels"], + surfname=config["surf_types"]["dentate"], allow_missing=True, ) - return specs + ) + return gii def get_final_subfields(): @@ -313,6 +341,7 @@ def get_final_qc(): def get_final_subj_output(): subj_output = [] subj_output.extend(get_final_spec()) + subj_output.extend(get_final_surf()) subj_output.extend(get_final_subfields()) subj_output.extend(get_final_coords()) subj_output.extend(get_final_transforms()) diff --git a/hippunfold/workflow/rules/warps.smk b/hippunfold/workflow/rules/warps.smk index cbdd11d8..62998078 100644 --- a/hippunfold/workflow/rules/warps.smk +++ b/hippunfold/workflow/rules/warps.smk @@ -167,7 +167,6 @@ rule create_warps_hipp: ), labelmap=get_labels_for_laplace, params: - interp_method="linear", gm_labels=lambda wildcards: config["laplace_labels"]["AP"]["gm"], epsilon=lambda wildcards: config["unfold_crop_epsilon_fractions"], resources: @@ -226,8 +225,6 @@ rule create_warps_hipp: hemi="{hemi}", suffix="create_warps-hipp.txt" ), - container: - config["singularity"]["autotop"] script: "../scripts/create_warps.py" @@ -294,7 +291,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 +349,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 cc923fb1..57e47246 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,74 +97,63 @@ 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 + (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, -) +# 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 epsilon = snakemake.params.epsilon (unfold_gx, unfold_gy, unfold_gz) = np.meshgrid( - np.linspace(0 + float(epsilon[0]), 1 - float(epsilon[0]), unfold_dims[0]), - np.linspace(0 + float(epsilon[1]), 1 - float(epsilon[1]), unfold_dims[1]), - np.linspace(0 + float(epsilon[2]), 1 - float(epsilon[2]), unfold_dims[2]), + np.linspace( + 0 + float(epsilon[0]), unfold_dims[0] - float(epsilon[0]), unfold_dims[0] + ), + np.linspace( + 0 + float(epsilon[1]), unfold_dims[1] - float(epsilon[1]), unfold_dims[1] + ), + np.linspace( + 0 + float(epsilon[2]), unfold_dims[2] - float(epsilon[2]), unfold_dims[2] + ), 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: @@ -180,10 +162,11 @@ 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[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 @@ -219,13 +202,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( @@ -250,7 +226,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 b5821b1c..c40958d4 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] From d85d5a65294654981669c525d468db0742b10651 Mon Sep 17 00:00:00 2001 From: Jordan DeKraker Date: Wed, 26 Jun 2024 10:25:30 -0400 Subject: [PATCH 2/6] naturalneighbor install to hippunfold_deps --- hippunfold/workflow/rules/warps.smk | 4 ++++ pyproject.toml | 1 - 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/hippunfold/workflow/rules/warps.smk b/hippunfold/workflow/rules/warps.smk index 62998078..47501f12 100644 --- a/hippunfold/workflow/rules/warps.smk +++ b/hippunfold/workflow/rules/warps.smk @@ -225,6 +225,8 @@ rule create_warps_hipp: hemi="{hemi}", suffix="create_warps-hipp.txt" ), + container: + config["singularity"]["autotop"] script: "../scripts/create_warps.py" @@ -349,6 +351,8 @@ rule create_warps_dentate: hemi="{hemi}", suffix="create_warps-dentate.txt" ), + container: + config["singularity"]["autotop"] script: "../scripts/create_warps.py" diff --git a/pyproject.toml b/pyproject.toml index c40958d4..b5821b1c 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -12,7 +12,6 @@ appdirs = "^1.4.4" Jinja2 = "^3.0.3" pygraphviz = "1.7" Pygments = "^2.10.0" -naturalneighbor = "0.2.1" [tool.poetry.dev-dependencies] From 24bc9e3bbb7bfadb333d47570acf290dc9b4e34a Mon Sep 17 00:00:00 2001 From: Jordan DeKraker Date: Wed, 26 Jun 2024 12:17:00 -0400 Subject: [PATCH 3/6] Update common.smk --- hippunfold/workflow/rules/common.smk | 51 ++++++++++++++++++---------- 1 file changed, 34 insertions(+), 17 deletions(-) diff --git a/hippunfold/workflow/rules/common.smk b/hippunfold/workflow/rules/common.smk index 4c568f95..6dc2fc60 100644 --- a/hippunfold/workflow/rules/common.smk +++ b/hippunfold/workflow/rules/common.smk @@ -24,23 +24,40 @@ def get_modality_suffix(modality): def get_final_spec(): - specs = expand( - bids( - root=root, - datatype="surf", - den="{density}", - space="{space}", - hemi="{hemi}", - label="{autotop}", - suffix="surfaces.spec", - **config["subj_wildcards"], - ), - density=config["output_density"], - space=ref_spaces, - hemi=config["hemi"], - autotop=config["autotop_labels"], - allow_missing=True, - ) + if len(config["hemi"]) == 2: + specs = expand( + bids( + root=root, + datatype="surf", + den="{density}", + space="{space}", + label="{autotop}", + suffix="surfaces.spec", + **config["subj_wildcards"], + ), + density=config["output_density"], + space=ref_spaces, + autotop=config["autotop_labels"], + allow_missing=True, + ) + else: + specs = expand( + bids( + root=root, + datatype="surf", + den="{density}", + space="{space}", + hemi="{hemi}", + label="{autotop}", + suffix="surfaces.spec", + **config["subj_wildcards"], + ), + density=config["output_density"], + space=ref_spaces, + hemi=config["hemi"], + autotop=config["autotop_labels"], + allow_missing=True, + ) return specs From 35806ce2c57bc071d705d1284b94a7027eeba33d Mon Sep 17 00:00:00 2001 From: Jordan DeKraker Date: Wed, 26 Jun 2024 12:24:55 -0400 Subject: [PATCH 4/6] removed epsilon and unused variables --- hippunfold/config/snakebids.yml | 4 ---- hippunfold/workflow/rules/warps.smk | 2 -- hippunfold/workflow/scripts/create_warps.py | 25 +++++---------------- 3 files changed, 5 insertions(+), 26 deletions(-) diff --git a/hippunfold/config/snakebids.yml b/hippunfold/config/snakebids.yml index d3a3f1bc..c374ab19 100644 --- a/hippunfold/config/snakebids.yml +++ b/hippunfold/config/snakebids.yml @@ -608,10 +608,6 @@ unfold_vol_ref: - '2.5' orient: RPI -unfold_crop_epsilon_fractions: - - 0 - - 0 - - 0.0625 #1/16 # space for uniform unfolded grid: # currently only used for interpolating hipp subfields on surface diff --git a/hippunfold/workflow/rules/warps.smk b/hippunfold/workflow/rules/warps.smk index 47501f12..1dd4a1d4 100644 --- a/hippunfold/workflow/rules/warps.smk +++ b/hippunfold/workflow/rules/warps.smk @@ -168,7 +168,6 @@ rule create_warps_hipp: labelmap=get_labels_for_laplace, params: gm_labels=lambda wildcards: config["laplace_labels"]["AP"]["gm"], - epsilon=lambda wildcards: config["unfold_crop_epsilon_fractions"], resources: mem_mb=16000, output: @@ -294,7 +293,6 @@ rule create_warps_dentate: labelmap=get_labels_for_laplace, params: gm_labels=lambda wildcards: config["laplace_labels"]["PD"]["sink"], - epsilon=lambda wildcards: config["unfold_crop_epsilon_fractions"], resources: mem_mb=16000, output: diff --git a/hippunfold/workflow/scripts/create_warps.py b/hippunfold/workflow/scripts/create_warps.py index 57e47246..7c6df89e 100644 --- a/hippunfold/workflow/scripts/create_warps.py +++ b/hippunfold/workflow/scripts/create_warps.py @@ -105,26 +105,7 @@ def summary(name, array): 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 -epsilon = snakemake.params.epsilon -(unfold_gx, unfold_gy, unfold_gz) = np.meshgrid( - np.linspace( - 0 + float(epsilon[0]), unfold_dims[0] - float(epsilon[0]), unfold_dims[0] - ), - np.linspace( - 0 + float(epsilon[1]), unfold_dims[1] - float(epsilon[1]), unfold_dims[1] - ), - np.linspace( - 0 + float(epsilon[2]), unfold_dims[2] - float(epsilon[2]), unfold_dims[2] - ), - indexing="ij", -) -summary("unfold_gx", unfold_gx) -summary("unfold_gy", unfold_gy) -summary("unfold_gz", unfold_gz) - -# perform the interpolation +# get unfolded grid points = np.stack( [ @@ -137,6 +118,10 @@ def summary(name, array): ], axis=1, ) +summary("points", points) + + +# perform the interpolation interp_ap = naturalneighbor.griddata( points, From c9fa333133173f46066cff3ebd246343903859a4 Mon Sep 17 00:00:00 2001 From: Jordan DeKraker Date: Wed, 26 Jun 2024 12:29:09 -0400 Subject: [PATCH 5/6] lint --- hippunfold/workflow/scripts/create_warps.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/hippunfold/workflow/scripts/create_warps.py b/hippunfold/workflow/scripts/create_warps.py index 7c6df89e..7ee0733a 100644 --- a/hippunfold/workflow/scripts/create_warps.py +++ b/hippunfold/workflow/scripts/create_warps.py @@ -105,7 +105,7 @@ def summary(name, array): coord_flat_pd_unnorm = coord_flat_pd * unfold_dims[1] coord_flat_io_unnorm = coord_flat_io * unfold_dims[2] -# get unfolded grid +# get unfolded grid points = np.stack( [ From 6149b182449b1fba0dbb945079a3887314430c4e Mon Sep 17 00:00:00 2001 From: Jordan DeKraker Date: Wed, 26 Jun 2024 14:02:28 -0400 Subject: [PATCH 6/6] hippunfold_deps:v0.5.1 --- hippunfold/config/snakebids.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/hippunfold/config/snakebids.yml b/hippunfold/config/snakebids.yml index c374ab19..b4768f2a 100644 --- a/hippunfold/config/snakebids.yml +++ b/hippunfold/config/snakebids.yml @@ -370,7 +370,7 @@ cifti_metric_types: #--- workflow specific configuration -- singularity: - autotop: 'docker://khanlab/hippunfold_deps:v0.5.0' + autotop: 'docker://khanlab/hippunfold_deps:v0.5.1' xfm_identity: resources/etc/identity_xfm.txt xfm_identity_itk: resources/etc/identity_xfm_itk.txt