diff --git a/hippunfold/config/snakebids.yml b/hippunfold/config/snakebids.yml index 0a2a97f6..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: @@ -358,7 +363,13 @@ cifti_metric_types: - gyrification.dscalar - curvature.dscalar - +# params for detecting and removing outlier vertices +outlierSmoothDist: + unfoldiso: 15 + 0p5mm: 7 + 1mm: 3 + 2mm: 1 +vertexOutlierThreshold: 3 @@ -601,7 +612,7 @@ unfold_vol_ref: unfold_crop_epsilon_fractions: - 0 - 0 - - 0.0625 #1/16 + - 1 # space for uniform unfolded grid: # currently only used for interpolating hipp subfields on surface 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/gifti.smk b/hippunfold/workflow/rules/gifti.smk index 2ae35706..32db0bd1 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=lambda wildcards: config["outlierSmoothDist"][wildcards.density], + 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=lambda wildcards: config["outlierSmoothDist"][wildcards.density], + 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..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 cc923fb1..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,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,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 @@ -219,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( @@ -250,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/hippunfold/workflow/scripts/fillbadvertices.py b/hippunfold/workflow/scripts/fillbadvertices.py new file mode 100644 index 00000000..ee9284c6 --- /dev/null +++ b/hippunfold/workflow/scripts/fillbadvertices.py @@ -0,0 +1,59 @@ +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) 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]