vertex cleaning, especially for DG

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

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

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:

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",
 )
 

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)