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)