From c46cb5e4b825e43a6bb50e4ffa1689faba2ece39 Mon Sep 17 00:00:00 2001 From: Jordan DeKraker Date: Wed, 12 Jun 2024 15:24:58 -0400 Subject: [PATCH 1/9] vertex cleaning, especially for DG --- hippunfold/config/snakebids.yml | 13 +++-- hippunfold/workflow/rules/gifti.smk | 26 +++++---- hippunfold/workflow/rules/warps.smk | 4 +- hippunfold/workflow/scripts/create_warps.py | 12 ++-- .../workflow/scripts/fillbadvertices.py | 58 +++++++++++++++++++ .../workflow/scripts/fillnanvertices.py | 24 -------- 6 files changed, 90 insertions(+), 47 deletions(-) create mode 100644 hippunfold/workflow/scripts/fillbadvertices.py delete mode 100644 hippunfold/workflow/scripts/fillnanvertices.py diff --git a/hippunfold/config/snakebids.yml b/hippunfold/config/snakebids.yml index 0a2a97f6..c391e851 100644 --- a/hippunfold/config/snakebids.yml +++ b/hippunfold/config/snakebids.yml @@ -348,6 +348,7 @@ gifti_metric_types: dentate: - gyrification.shape - curvature.shape + - thickness.shape cifti_metric_types: hipp: @@ -357,8 +358,10 @@ cifti_metric_types: dentate: - gyrification.dscalar - curvature.dscalar + - thickness.dscalar - +outlierSmoothDist: 15 +vertexOutlierThreshold: 3 @@ -564,7 +567,7 @@ unfold_vol_ref: dims: - '256' - '128' - - '16' + - '8' voxdims: - '0.15625' - '0.15625' @@ -599,9 +602,9 @@ unfold_vol_ref: orient: RPI unfold_crop_epsilon_fractions: - - 0 - - 0 - - 0.0625 #1/16 + - 1 + - 1 + - 1 # space for uniform unfolded grid: # currently only used for interpolating hipp subfields on surface diff --git a/hippunfold/workflow/rules/gifti.smk b/hippunfold/workflow/rules/gifti.smk index 2ae35706..fcbf0a07 100644 --- a/hippunfold/workflow/rules/gifti.smk +++ b/hippunfold/workflow/rules/gifti.smk @@ -181,7 +181,7 @@ rule warp_gii_unfold2corobl1: datatype="surf", den="{density}", suffix="{surfname}.surf.gii", - desc="nonancorrect", + desc="nobadcorrect", space="corobl", unfoldreg="none", hemi="{hemi}", @@ -198,21 +198,24 @@ rule warp_gii_unfold2corobl1: " -surface-secondary-type {params.secondary_type}" -# previous rule seems to be where nan vertices emerge, so we'll correct them here immediately after -rule correct_nan_vertices1: +# previous rule seems to be where bad vertices emerge, so we'll correct them here immediately after +rule correct_bad_vertices1: input: gii=bids( root=work, datatype="surf", den="{density}", suffix="{surfname}.surf.gii", - desc="nonancorrect", + desc="nobadcorrect", space="corobl", unfoldreg="none", hemi="{hemi}", label="hipp", **config["subj_wildcards"] ), + params: + dist=config["outlierSmoothDist"], + threshold=config["vertexOutlierThreshold"], output: gii=bids( root=work, @@ -230,7 +233,7 @@ rule correct_nan_vertices1: container: config["singularity"]["autotop"] script: - "../scripts/fillnanvertices.py" + "../scripts/fillbadvertices.py" # morphological features, calculated in native space: @@ -730,7 +733,7 @@ rule warp_gii_unfold2corobl2: datatype="surf", den="{density}", suffix="{surfname}.surf.gii", - desc="nonancorrect", + desc="nobadcorrect", space="corobl", hemi="{hemi}", label="{autotop}", @@ -746,20 +749,23 @@ rule warp_gii_unfold2corobl2: " -surface-secondary-type {params.secondary_type}" -# previous rule seems to be where nan vertices emerge, so we'll correct them here immediately after -rule correct_nan_vertices2: +# previous rule seems to be where bad vertices emerge, so we'll correct them here immediately after +rule correct_bad_vertices2: input: gii=bids( root=work, datatype="surf", den="{density}", suffix="{surfname}.surf.gii", - desc="nonancorrect", + desc="nobadcorrect", space="corobl", hemi="{hemi}", label="{autotop}", **config["subj_wildcards"] ), + params: + dist=config["outlierSmoothDist"], + threshold=config["vertexOutlierThreshold"], output: gii=bids( root=work, @@ -776,7 +782,7 @@ rule correct_nan_vertices2: container: config["singularity"]["autotop"] script: - "../scripts/fillnanvertices.py" + "../scripts/fillbadvertices.py" # needed for if native_modality is corobl diff --git a/hippunfold/workflow/rules/warps.smk b/hippunfold/workflow/rules/warps.smk index cbdd11d8..8e328377 100644 --- a/hippunfold/workflow/rules/warps.smk +++ b/hippunfold/workflow/rules/warps.smk @@ -167,7 +167,7 @@ rule create_warps_hipp: ), labelmap=get_labels_for_laplace, params: - interp_method="linear", + interp_method="nearest", gm_labels=lambda wildcards: config["laplace_labels"]["AP"]["gm"], epsilon=lambda wildcards: config["unfold_crop_epsilon_fractions"], resources: @@ -294,7 +294,7 @@ rule create_warps_dentate: ), labelmap=get_labels_for_laplace, params: - interp_method="linear", + interp_method="nearest", gm_labels=lambda wildcards: config["laplace_labels"]["PD"]["sink"], epsilon=lambda wildcards: config["unfold_crop_epsilon_fractions"], resources: diff --git a/hippunfold/workflow/scripts/create_warps.py b/hippunfold/workflow/scripts/create_warps.py index cc923fb1..7db55e14 100644 --- a/hippunfold/workflow/scripts/create_warps.py +++ b/hippunfold/workflow/scripts/create_warps.py @@ -111,18 +111,18 @@ def summary(name, array): # 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, + 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, ) # 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", ) diff --git a/hippunfold/workflow/scripts/fillbadvertices.py b/hippunfold/workflow/scripts/fillbadvertices.py new file mode 100644 index 00000000..32f8f362 --- /dev/null +++ b/hippunfold/workflow/scripts/fillbadvertices.py @@ -0,0 +1,58 @@ +import nibabel as nib +import numpy as np +from scipy.stats import zscore +import copy + +SDthreshold = snakemake.params.threshold +iters = snakemake.params.dist + +gii = nib.load(snakemake.input.gii) +varr = gii.get_arrays_from_intent("NIFTI_INTENT_POINTSET")[0] +V = varr.data +farr = gii.get_arrays_from_intent("NIFTI_INTENT_TRIANGLE")[0] +F = farr.data + + +# find local outliers by smoothing and then substracting from original +# https://github.com/MICA-MNI/hippomaps/blob/master/hippomaps/utils.py +def avg_neighbours(F, cdat, n): + frows = np.where(F == n)[0] + v = np.unique(F[frows, :]) + cdat = np.reshape(cdat, (len(cdat), -1)) + out = np.nanmean(cdat[v,:], 0) + return out + +def surfdat_smooth(F, cdata, iters=1): + sz = cdata.shape + cdata = cdata.reshape(cdata.shape[0], -1) + cdata_smooth = copy.deepcopy(cdata) + for i in range(iters): + for n in range(len(cdata)): + cdata_smooth[n,:] = avg_neighbours(F, cdata, n) + cdata = copy.deepcopy(cdata_smooth) + return cdata_smooth.reshape(sz) + +Vsmooth = surfdat_smooth(F, V, iters=iters) +Vdiffs = V - Vsmooth +Vdists = np.sqrt((Vdiffs[:,0])**2 + (Vdiffs[:,1])**2 + (Vdiffs[:,2])**2) +Vzscored = zscore(Vdists) +outliers = ((Vzscored > SDthreshold) | (Vzscored < -SDthreshold)) +V[outliers,:] = np.nan + + + +# most nans should be just isolated points, but in case there is an island of nans this will erode it, replacing with decent (but not perfect) guesses of where vertices should be +while np.isnan(np.sum(V)): + # index of vertices containing nan + i = np.where(np.isnan(V)) + ii = np.unique(i[0]) + # replace with the nanmean of neighbouring vertices + newV = V + for n in ii: + f = np.where(F == n) + v = F[f[0]] + vv = np.unique(v) + newV[n, :] = np.nanmean(V[vv, :], 0) + V = newV + +nib.save(gii, snakemake.output.gii) diff --git a/hippunfold/workflow/scripts/fillnanvertices.py b/hippunfold/workflow/scripts/fillnanvertices.py deleted file mode 100644 index 22e7b1db..00000000 --- a/hippunfold/workflow/scripts/fillnanvertices.py +++ /dev/null @@ -1,24 +0,0 @@ -import nibabel as nib -import numpy as np - -gii = nib.load(snakemake.input.gii) -varr = gii.get_arrays_from_intent("NIFTI_INTENT_POINTSET")[0] -V = varr.data -farr = gii.get_arrays_from_intent("NIFTI_INTENT_TRIANGLE")[0] -F = farr.data - -# most nans should be just isolated points, but in case there is an island of nans this will erode it, replacing with decent (but not perfect) guesses of where vertices should be -while np.isnan(np.sum(V)): - # index of vertices containing nan - i = np.where(np.isnan(V)) - ii = np.unique(i[0]) - # replace with the nanmean of neighbouring vertices - newV = V - for n in ii: - f = np.where(F == n) - v = F[f[0]] - vv = np.unique(v) - newV[n, :] = np.nanmean(V[vv, :], 0) - V = newV - -nib.save(gii, snakemake.output.gii) From 3f73bf3d59002883d7bca1ca1cabaedde228d19f Mon Sep 17 00:00:00 2001 From: Jordan DeKraker Date: Wed, 12 Jun 2024 15:31:42 -0400 Subject: [PATCH 2/9] effectively resotres parameters to previous defaults --- hippunfold/config/snakebids.yml | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/hippunfold/config/snakebids.yml b/hippunfold/config/snakebids.yml index c391e851..c347a312 100644 --- a/hippunfold/config/snakebids.yml +++ b/hippunfold/config/snakebids.yml @@ -348,7 +348,6 @@ gifti_metric_types: dentate: - gyrification.shape - curvature.shape - - thickness.shape cifti_metric_types: hipp: @@ -358,7 +357,6 @@ cifti_metric_types: dentate: - gyrification.dscalar - curvature.dscalar - - thickness.dscalar outlierSmoothDist: 15 vertexOutlierThreshold: 3 @@ -567,7 +565,7 @@ unfold_vol_ref: dims: - '256' - '128' - - '8' + - '16' voxdims: - '0.15625' - '0.15625' @@ -602,8 +600,8 @@ unfold_vol_ref: orient: RPI unfold_crop_epsilon_fractions: - - 1 - - 1 + - 0 + - 0 - 1 # space for uniform unfolded grid: From 8f1367c28ff2a63a07dc4e06a39230223c4a23c6 Mon Sep 17 00:00:00 2001 From: Jordan DeKraker Date: Wed, 12 Jun 2024 15:35:38 -0400 Subject: [PATCH 3/9] lint --- hippunfold/workflow/scripts/create_warps.py | 21 +++++++++++++------ .../workflow/scripts/fillbadvertices.py | 13 ++++++------ 2 files changed, 22 insertions(+), 12 deletions(-) diff --git a/hippunfold/workflow/scripts/create_warps.py b/hippunfold/workflow/scripts/create_warps.py index 7db55e14..6562d753 100644 --- a/hippunfold/workflow/scripts/create_warps.py +++ b/hippunfold/workflow/scripts/create_warps.py @@ -111,18 +111,27 @@ def summary(name, array): # 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, + 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, ) # 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]), + 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", ) diff --git a/hippunfold/workflow/scripts/fillbadvertices.py b/hippunfold/workflow/scripts/fillbadvertices.py index 32f8f362..ee9284c6 100644 --- a/hippunfold/workflow/scripts/fillbadvertices.py +++ b/hippunfold/workflow/scripts/fillbadvertices.py @@ -19,26 +19,27 @@ def avg_neighbours(F, cdat, n): frows = np.where(F == n)[0] v = np.unique(F[frows, :]) cdat = np.reshape(cdat, (len(cdat), -1)) - out = np.nanmean(cdat[v,:], 0) + out = np.nanmean(cdat[v, :], 0) return out + def surfdat_smooth(F, cdata, iters=1): sz = cdata.shape cdata = cdata.reshape(cdata.shape[0], -1) cdata_smooth = copy.deepcopy(cdata) for i in range(iters): for n in range(len(cdata)): - cdata_smooth[n,:] = avg_neighbours(F, cdata, n) + cdata_smooth[n, :] = avg_neighbours(F, cdata, n) cdata = copy.deepcopy(cdata_smooth) return cdata_smooth.reshape(sz) + Vsmooth = surfdat_smooth(F, V, iters=iters) Vdiffs = V - Vsmooth -Vdists = np.sqrt((Vdiffs[:,0])**2 + (Vdiffs[:,1])**2 + (Vdiffs[:,2])**2) +Vdists = np.sqrt((Vdiffs[:, 0]) ** 2 + (Vdiffs[:, 1]) ** 2 + (Vdiffs[:, 2]) ** 2) Vzscored = zscore(Vdists) -outliers = ((Vzscored > SDthreshold) | (Vzscored < -SDthreshold)) -V[outliers,:] = np.nan - +outliers = (Vzscored > SDthreshold) | (Vzscored < -SDthreshold) +V[outliers, :] = np.nan # most nans should be just isolated points, but in case there is an island of nans this will erode it, replacing with decent (but not perfect) guesses of where vertices should be From 426058aed37622dd785b9429fffb356852cfbc53 Mon Sep 17 00:00:00 2001 From: Jordan DeKraker Date: Wed, 12 Jun 2024 15:36:43 -0400 Subject: [PATCH 4/9] note --- hippunfold/config/snakebids.yml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/hippunfold/config/snakebids.yml b/hippunfold/config/snakebids.yml index c347a312..02726824 100644 --- a/hippunfold/config/snakebids.yml +++ b/hippunfold/config/snakebids.yml @@ -358,6 +358,8 @@ cifti_metric_types: - gyrification.dscalar - curvature.dscalar +# params for detecting and removing outlier vertices +# NOTE this is currently dependent on surface density outlierSmoothDist: 15 vertexOutlierThreshold: 3 From c8190729fd324c54d24478c5207771f3a1fc6861 Mon Sep 17 00:00:00 2001 From: Jordan DeKraker Date: Wed, 12 Jun 2024 15:44:27 -0400 Subject: [PATCH 5/9] restoreing one more parameter --- hippunfold/workflow/rules/warps.smk | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/hippunfold/workflow/rules/warps.smk b/hippunfold/workflow/rules/warps.smk index 8e328377..cbdd11d8 100644 --- a/hippunfold/workflow/rules/warps.smk +++ b/hippunfold/workflow/rules/warps.smk @@ -167,7 +167,7 @@ rule create_warps_hipp: ), labelmap=get_labels_for_laplace, params: - interp_method="nearest", + interp_method="linear", gm_labels=lambda wildcards: config["laplace_labels"]["AP"]["gm"], epsilon=lambda wildcards: config["unfold_crop_epsilon_fractions"], resources: @@ -294,7 +294,7 @@ rule create_warps_dentate: ), labelmap=get_labels_for_laplace, params: - interp_method="nearest", + interp_method="linear", gm_labels=lambda wildcards: config["laplace_labels"]["PD"]["sink"], epsilon=lambda wildcards: config["unfold_crop_epsilon_fractions"], resources: From f863be6dbd487f86fee8a94e4312de73d31a4bec Mon Sep 17 00:00:00 2001 From: Jordan DeKraker Date: Thu, 13 Jun 2024 11:32:58 -0400 Subject: [PATCH 6/9] addresses configurable DG surface --- hippunfold/workflow/rules/common.smk | 93 ++++++++++++++++++---------- 1 file changed, 59 insertions(+), 34 deletions(-) diff --git a/hippunfold/workflow/rules/common.smk b/hippunfold/workflow/rules/common.smk index b0e5db03..41e29196 100644 --- a/hippunfold/workflow/rules/common.smk +++ b/hippunfold/workflow/rules/common.smk @@ -24,43 +24,67 @@ def get_modality_suffix(modality): def get_final_spec(): - 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, - ) + 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}", + **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, + )) + gii.extend(expand( + bids( + root=root, + datatype="surf", + den="{density}", + suffix="{surfname}.surf.gii", + space="{space}", + hemi="{hemi}", + label="{autotop}", + **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 gii + + def get_final_subfields(): return expand( bids( @@ -313,6 +337,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()) From cbf3b0a5f6065dc398facb27680658fc2ec81dbc Mon Sep 17 00:00:00 2001 From: Jordan DeKraker Date: Thu, 13 Jun 2024 11:35:53 -0400 Subject: [PATCH 7/9] lint --- hippunfold/workflow/rules/common.smk | 76 +++++++++++++++------------- 1 file changed, 40 insertions(+), 36 deletions(-) diff --git a/hippunfold/workflow/rules/common.smk b/hippunfold/workflow/rules/common.smk index 41e29196..4c568f95 100644 --- a/hippunfold/workflow/rules/common.smk +++ b/hippunfold/workflow/rules/common.smk @@ -46,42 +46,46 @@ def get_final_spec(): 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}", - **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, - )) - gii.extend(expand( - bids( - root=root, - datatype="surf", - den="{density}", - suffix="{surfname}.surf.gii", - space="{space}", - hemi="{hemi}", - label="{autotop}", - **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, - )) + gii.extend( + expand( + bids( + root=root, + datatype="surf", + den="{density}", + suffix="{surfname}.surf.gii", + space="{space}", + hemi="{hemi}", + label="{autotop}", + **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, + ) + ) + gii.extend( + expand( + bids( + root=root, + datatype="surf", + den="{density}", + suffix="{surfname}.surf.gii", + space="{space}", + hemi="{hemi}", + label="{autotop}", + **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 gii From 92113225b9ea4cbb6353e2651fe0b66af54080f8 Mon Sep 17 00:00:00 2001 From: Jordan DeKraker Date: Thu, 13 Jun 2024 11:36:28 -0400 Subject: [PATCH 8/9] surface configuration now more editable --- hippunfold/config/snakebids.yml | 14 +++++++++++--- hippunfold/workflow/rules/gifti.smk | 4 ++-- 2 files changed, 13 insertions(+), 5 deletions(-) diff --git a/hippunfold/config/snakebids.yml b/hippunfold/config/snakebids.yml index 02726824..dee56824 100644 --- a/hippunfold/config/snakebids.yml +++ b/hippunfold/config/snakebids.yml @@ -327,11 +327,14 @@ parse_args: - synthseg_v0.2 - neonateT1w_v2 + + +# --- surface specific configuration -- + autotop_labels: - 'hipp' - 'dentate' - surf_types: hipp: - midthickness @@ -339,6 +342,8 @@ surf_types: - outer dentate: - midthickness +# - inner +# - outer gifti_metric_types: hipp: @@ -359,8 +364,11 @@ cifti_metric_types: - curvature.dscalar # params for detecting and removing outlier vertices -# NOTE this is currently dependent on surface density -outlierSmoothDist: 15 +outlierSmoothDist: + unfoldiso: 15 + 0p5mm: 7 + 1mm: 3 + 2mm: 1 vertexOutlierThreshold: 3 diff --git a/hippunfold/workflow/rules/gifti.smk b/hippunfold/workflow/rules/gifti.smk index fcbf0a07..32db0bd1 100644 --- a/hippunfold/workflow/rules/gifti.smk +++ b/hippunfold/workflow/rules/gifti.smk @@ -214,7 +214,7 @@ rule correct_bad_vertices1: **config["subj_wildcards"] ), params: - dist=config["outlierSmoothDist"], + dist=lambda wildcards: config["outlierSmoothDist"][wildcards.density], threshold=config["vertexOutlierThreshold"], output: gii=bids( @@ -764,7 +764,7 @@ rule correct_bad_vertices2: **config["subj_wildcards"] ), params: - dist=config["outlierSmoothDist"], + dist=lambda wildcards: config["outlierSmoothDist"][wildcards.density], threshold=config["vertexOutlierThreshold"], output: gii=bids( From 6d9b08a39926dbc5dc52569600ab6862f3bb386b Mon Sep 17 00:00:00 2001 From: Jordan DeKraker Date: Wed, 26 Jun 2024 09:18:02 -0400 Subject: [PATCH 9/9] natural interpolation is fast and may solve issues with misplaced vertices! --- hippunfold/workflow/rules/warps.smk | 5 - hippunfold/workflow/scripts/create_warps.py | 113 ++++++++------------ pyproject.toml | 1 + 3 files changed, 43 insertions(+), 76 deletions(-) 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]