diff --git a/hippunfold/config/snakebids.yml b/hippunfold/config/snakebids.yml index 96e08641..e6f75028 100644 --- a/hippunfold/config/snakebids.yml +++ b/hippunfold/config/snakebids.yml @@ -278,6 +278,7 @@ parse_args: - synthseg_v0.2 - neonateT1w_v2 + --enable-bids-validation: help: | Enable validation of BIDS dataset. BIDS validation would be performed @@ -294,7 +295,6 @@ parse_args: version: "1.5.2-pre.2" - # --- surface specific configuration -- autotop_labels: @@ -308,6 +308,8 @@ surf_types: - outer dentate: - midthickness +# - inner +# - outer gifti_metric_types: hipp: @@ -327,7 +329,14 @@ cifti_metric_types: - gyrification.dscalar - curvature.dscalar - +# params for detecting and removing outlier vertices +outlier_opts: + outlierSmoothDist: + unfoldiso: 15 + 0p5mm: 7 + 1mm: 3 + 2mm: 1 + vertexOutlierThreshold: 3 @@ -373,6 +382,7 @@ template_based_segmentation: - R - L + template_files: CITI168: T1w: T1w_head_700um.nii.gz @@ -427,12 +437,14 @@ template_files: dseg: tpl-ABAv3_hemi-{hemi}_space-corobl_desc-tissuemanual_dseg.nii.gz coords: tpl-ABAv3_dir-{dir}_hemi-{hemi}_space-corobl_label-{autotop}_desc-laplace_coords.nii.gz bigbrain: - T1w: tpl-bbhist_100um_T1w.nii.gz - xfm_corobl: tpl-bbhist_from-native_to-corobl_type-itk_affine.txt - crop_ref: tpl-bbhist_hemi-{hemi}_space-corobl_desc-tissuemanual_40um_dseg.nii.gz - Mask_crop: tpl-bbhist_hemi-{hemi}_space-corobl_desc-tissuemanual_40um_dseg.nii.gz - dseg: tpl-bbhist_hemi-{hemi}_space-corobl_desc-tissuemanual_40um_dseg.nii.gz - coords: tpl-bbhist_dir-{dir}_hemi-{hemi}_space-corobl_label-{autotop}_desc-laplace_coords.nii.gz + T1w: tpl-bigbrain_1mm_merker.nii.gz + T2w: tpl-bigbrain_1mm_merker.nii.gz + xfm_corobl: tpl-bigbrain_from-native_to-corobl_type-itk_affine.txt + crop_ref: tpl-bigbrain_hemi-{hemi}_space-corobl_100um_merker.nii.gz + Mask_crop: tpl-bigbrain_hemi-{hemi}_space-corobl_desc-tissuemanual_100um_dseg.nii.gz + dseg: tpl-bigbrain_hemi-{hemi}_space-corobl_desc-tissuemanual_100um_dseg.nii.gz + coords: tpl-bigbrain_dir-{dir}_hemi-{hemi}_space-corobl_label-{autotop}_desc-laplace_coords.nii.gz + atlas_files: multihist7: @@ -591,6 +603,11 @@ unfold_vol_ref: - '2.5' orient: RPI +unfold_crop_epsilon_fractions: + - 0 + - 0 + - 1 #in unfolded space voxels + # 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 61fe338f..91fd1c23 100644 --- a/hippunfold/workflow/rules/common.smk +++ b/hippunfold/workflow/rules/common.smk @@ -40,55 +40,17 @@ def get_modality_suffix(modality): def get_final_spec(): - if len(config["hemi"]) == 2: - specs = inputs[config["modality"]].expand( - bids( - root=root, - datatype="surf", - den="{density}", - space="{space}", - label="{autotop}", - suffix="surfaces.spec", - **inputs.subj_wildcards, - ), - density=config["output_density"], - space=ref_spaces, - autotop=config["autotop_labels"], - allow_missing=True, - ) - else: - specs = inputs[config["modality"]].expand( - bids( - root=root, - datatype="surf", - den="{density}", - space="{space}", - hemi="{hemi}", - label="{autotop}", - suffix="surfaces.spec", - **inputs.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( + specs = [] + specs.extend( inputs[config["modality"]].expand( bids( root=root, datatype="surf", den="{density}", - suffix="{surfname}.surf.gii", space="{space}", hemi="{hemi}", label="{autotop}", + suffix="surfaces.spec", **inputs.subj_wildcards, ), density=config["output_density"], @@ -99,16 +61,16 @@ def get_final_surf(): allow_missing=True, ) ) - gii.extend( + specs.extend( inputs[config["modality"]].expand( bids( root=root, datatype="surf", den="{density}", - suffix="{surfname}.surf.gii", space="{space}", hemi="{hemi}", label="{autotop}", + suffix="surfaces.spec", **inputs.subj_wildcards, ), density=config["output_density"], @@ -119,7 +81,7 @@ def get_final_surf(): allow_missing=True, ) ) - return gii + return specs def get_final_subfields(): @@ -373,7 +335,6 @@ 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 515b1d16..8f26acca 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,26 @@ 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", **inputs.subj_wildcards ), + params: + dist=lambda wildcards: config["outlier_opts"]["outlierSmoothDist"][ + wildcards.density + ], + threshold=config["outlier_opts"]["vertexOutlierThreshold"], output: gii=bids( root=work, @@ -230,7 +235,7 @@ rule correct_nan_vertices1: container: config["singularity"]["autotop"] script: - "../scripts/fillnanvertices.py" + "../scripts/fillbadvertices.py" # morphological features, calculated in native space: @@ -730,7 +735,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 +751,25 @@ 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}", **inputs.subj_wildcards ), + params: + dist=lambda wildcards: config["outlier_opts"]["outlierSmoothDist"][ + wildcards.density + ], + threshold=config["outlier_opts"]["vertexOutlierThreshold"], output: gii=bids( root=work, @@ -776,7 +786,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/scripts/create_warps.py b/hippunfold/workflow/scripts/create_warps.py index 7ee0733a..6f9105d3 100644 --- a/hippunfold/workflow/scripts/create_warps.py +++ b/hippunfold/workflow/scripts/create_warps.py @@ -100,23 +100,30 @@ def summary(name, array): # 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 -# 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 - -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, +# 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, +) + +# 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("points", points) 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)