BB template & vertex outliers (#333)

* vertex cleaning, especially for DG * effectively resotres parameters to previous defaults * lint * note * restoreing one more parameter * addresses configurable DG surface * lint * surface configuration now more editable * bigbrain filled; not on OSF * renaming & lower resolution * grouped outlier opts * lint * lint * removed redudant hemi conditional; replaced with hippVSdentate * bugfixes --------- Co-authored-by: Jordan DeKraker <[email protected]>
khanlab · Dec 17, 2024 · 48ebd46 · 48ebd46
1 parent f642963
commit 48ebd46
Show file tree

Hide file tree

Showing 6 changed files with 134 additions and 104 deletions.
diff --git 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
@@ -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
@@ -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
@@ -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
@@ -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)