enh: Add -fill-wm-holes option to subdivide-brain-image

Applications/src/subdivide-brain-image.cc

@@ -2,7 +2,7 @@
  * Medical Image Registration ToolKit (MIRTK)
  *
  * Copyright 2016 Imperial College London
- * Copyright 2016 Andreas Schuh
+ * Copyright 2016-2020 Andreas Schuh
  *
  * Licensed under the Apache License, Version 2.0 (the "License");
  * you may not use this file except in compliance with the License.
@@ -125,6 +125,8 @@ void PrintHelp(const char *name)
   cout << "      No. of iterations used to close holes in cerebellum segmentation. (default: 0)\n";
   cout << "  -brainstem-and-cerebellum, -cerebellum-and-brainstem, -bscb, -cbbs [on|off]\n";
   cout << "      Whether to merge brainstem and cerebellum. (default: off)\n";
+  cout << "  -fill-wm-holes [on|off]\n";
+  cout << "      Whether to fill interhemisphere holes in WM that are closed in xy plane. (default: off)\n";
   PrintStandardOptions(cout);
   cout << endl;
 }
@@ -676,6 +678,39 @@ Plane MedialCuttingPlane(const ByteImage &regions)
   return plane;
 }
 
+// -----------------------------------------------------------------------------
+bool IsInterhemisphereWhiteMatterHole(const ByteImage &regions, int k, const GreyImage &cc, GreyPixel hid)
+{
+  NeighborhoodOffsets offsets(&cc, CONNECTIVITY_4);
+  bool touches_lh = false;
+  bool touches_rh = false;
+  bool touches_non_wm = false;
+  const int nvox = cc.NumberOfVoxels();
+  const int k_nvox = k * nvox;
+  for (int vox = 0; vox < nvox; ++vox) {
+    if (cc(vox) == hid) {
+      for (int n = 0; n < offsets.Size(); ++n) {
+        const int nbr = vox + offsets(n);
+        if (cc.IsInside(nbr)) {
+          if (cc(nbr) == 0) {
+            const GreyPixel label = regions(k_nvox + nbr);
+            if (label == LH) {
+              touches_lh = true;
+            } else if (label == RH) {
+              touches_rh = true;
+            } else if (label != BG) {
+              touches_non_wm = true;
+            }
+          }
+        } else {
+          return false;
+        }
+      }
+    }
+  }
+  return touches_lh && touches_rh && !touches_non_wm;
+}
+
 // -----------------------------------------------------------------------------
 bool TouchesOtherHemisphere(const GreyImage &cc, const GreyImage &other, int vox)
 {
@@ -965,6 +1000,7 @@ int main(int argc, char *argv[])
   int  ymargin     = 0;     // Margin in y direction after resampling in no. of voxels
   int  zmargin     = 0;     // Margin in z direction after resampling in no. of voxels
   bool merge_bs_cb = false; // Whether to merge brainstem and cerebellum segments
+  bool fill_wm_holes = false; // Whether to fill holes combined LH/RH WM region
 
   if (NUM_POSARGS == 1) {
     output_name = POSARG(1);
@@ -1051,6 +1087,7 @@ int main(int argc, char *argv[])
     else if (OPTION("-cerebellum-closing")) {
       PARSE_ARGUMENT(cb_closing);
     }
+    else HANDLE_BOOLEAN_OPTION("fill-wm-holes", fill_wm_holes);
     else HANDLE_COMMON_OR_UNKNOWN_OPTION();
   }
 
@@ -1192,6 +1229,13 @@ int main(int argc, char *argv[])
 
   if (debug) {
     regions.Write("debug_regions.nii.gz");
+    if (debug > 1) {
+      rhmask.Write("debug_regions_rhmask.nii.gz");
+      lhmask.Write("debug_regions_lhmask.nii.gz");
+      bsmask.Write("debug_regions_bsmask.nii.gz");
+      cbmask.Write("debug_regions_cbmask.nii.gz");
+      sbmask.Write("debug_regions_sbmask.nii.gz");
+    }
   }
 
   if (sb_closing > 0 || bs_closing > 0 || cb_closing > 0) {
@@ -1420,6 +1464,50 @@ int main(int argc, char *argv[])
     bounds[0] = max(0,               bounds[0] - 20);
     bounds[1] = min(regions.X() - 1, bounds[1] + 20);
 
+    // ---------------------------------------------------------------------------
+    // Fill closed holes in xy slices of union of left and right WM
+    if (fill_wm_holes) {
+      GreyImage holes(attr._x, attr._y);
+      for (int k = bounds[4]; k <= bounds[5]; ++k) {
+        // Create binary BG mask
+        for (int j = 0; j < attr._y; ++j)
+        for (int i = 0; i < attr._x; ++i) {
+          if (regions(i, j, k) == BG) {
+            holes(i, j) = 1;
+          } else {
+            holes(i, j) = 0;
+          }
+        }
+        // Label connected BG components
+        ConnectedComponents<GreyPixel> cc(CC_LargestFirst, CONNECTIVITY_4);
+        cc.Input(&holes);
+        cc.Output(&holes);
+        cc.Run();
+        // Discard BG components which touch non-WM at their boundary
+        for (GreyPixel hid = 1; hid <= cc.NumberOfComponents(); ++hid) {
+          if (!IsInterhemisphereWhiteMatterHole(regions, k, holes, hid)) {
+            cc.DeleteComponent(hid);
+          }
+        }
+        // Fill remaining BG components with RH/LH labels
+        for (int j = 0; j < attr._y; ++j)
+        for (int i = 0; i < attr._x; ++i) {
+          if (holes(i, j)) {
+            p = Point(i, j, k);
+            regions.ImageToWorld(p);
+            if (rl_plane.SignedDistance(p) < 0.) {
+              regions(i, j, k) = LH;
+            } else {
+              regions(i, j, k) = RH;
+            }
+          }
+        }
+      }
+      if (debug) {
+        regions.Write("debug_output+wm_closed.nii.gz");
+      }
+    }
+
     // -------------------------------------------------------------------------
     // 1. Fill small holes in xz slices nearby interhemispheric bounding box