From bb48bea4b718c4b7ea20e3be03120455f5155448 Mon Sep 17 00:00:00 2001
From: Jennings Zhang <jenni_zh@protonmail.com>
Date: Thu, 3 Mar 2022 23:31:48 -0500
Subject: [PATCH] Add mincblur

---
 requirements.txt |  1 +
 setup.py         |  2 +-
 skimc.py         | 56 ++++++++++++++++++++++++++++++++++++------------
 3 files changed, 44 insertions(+), 15 deletions(-)

diff --git a/requirements.txt b/requirements.txt
index 11cad27..a601aa7 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -3,3 +3,4 @@ pybicpl~=0.3.0
 numpy~=1.22.2
 nibabel~=3.2.2
 loguru~=0.6.0
+pycivet==0.0.3
diff --git a/setup.py b/setup.py
index e8992a6..2830721 100644
--- a/setup.py
+++ b/setup.py
@@ -2,7 +2,7 @@
 
 setup(
     name             = 'ep-skimage-mcubes-mni',
-    version          = '1.1.0',
+    version          = '1.2.0',
     description      = 'Marching-cubes implementation from scikit-image',
     author           = 'Jennings Zhang',
     author_email     = 'Jennings.Zhang@childrens.harvard.edu',
diff --git a/skimc.py b/skimc.py
index b354886..d461a1a 100755
--- a/skimc.py
+++ b/skimc.py
@@ -1,16 +1,21 @@
 #!/usr/bin/env python
 
 import os
-import nibabel as nib
-from nibabel.affines import apply_affine
-from skimage import measure
-from bicpl import PolygonObj
-
-from pathlib import Path
+import subprocess as sp
 from argparse import ArgumentParser, Namespace
 from concurrent.futures import ThreadPoolExecutor
-from loguru import logger
+from contextlib import contextmanager
+from pathlib import Path
+from tempfile import NamedTemporaryFile
+from typing import Sequence
+
+import nibabel as nib
+from bicpl import PolygonObj
 from chris_plugin import chris_plugin, PathMapper
+from civet.minc import Mask
+from loguru import logger
+from nibabel.affines import apply_affine
+from skimage import measure
 
 parser = ArgumentParser(description='cli description')
 parser.add_argument('-l', '--level', default=None, type=float,
@@ -30,19 +35,38 @@
                          ' method will be used.')
 parser.add_argument('-p', '--pattern', default='**/*.mnc',
                     help='pattern for file names to include')
+parser.add_argument('-f', '--fwhm', type=float, default=3.0,
+                    help='Amount of blurring (mincblur -fwhm) to do as preprocessing. '
+                         'A value of 0 skips blurring.')
+
+
+def __quiet_shell(cmd: Sequence[str | os.PathLike]) -> None:
+    sp.run(cmd, check=True, stdout=sp.DEVNULL, stderr=sp.STDOUT)
+
+
+@contextmanager
+def mincblur(mask_path: Path, fwhm: float) -> str:
+    if fwhm == 0:
+        yield str(mask_path)
+        return
+    with NamedTemporaryFile(suffix='.mnc') as tmp:
+        Mask(mask_path).resamplef64().mincblur(fwhm=fwhm).save(tmp.name, shell=__quiet_shell)
+        yield tmp.name
 
 
 def mcubes(mask_path: Path, surface_path: Path,
+           fwhm: float,
            level: float,
            spacing: tuple[float, float, float],
            step_size: int,
            method: str):
-    mask = nib.load(mask_path)
-    data = mask.get_fdata()
+    with mincblur(mask_path, fwhm) as blurred:
+        blurred_mask = nib.load(blurred)
+    data = blurred_mask.get_fdata()
     verts, faces, normals, values = measure.marching_cubes(
         data, level=level, spacing=spacing, step_size=step_size, method=method, allow_degenerate=False
     )
-    transformed_verts = apply_affine(mask.affine, verts)
+    transformed_verts = apply_affine(blurred_mask.affine, verts)
     obj = PolygonObj.from_data(transformed_verts, faces, normals)
     obj.save(surface_path)
     logger.info('Completed: {} => {}', mask_path, surface_path)
@@ -58,17 +82,21 @@ def mcubes(mask_path: Path, surface_path: Path,
 def main(options: Namespace, inputdir: Path, outputdir: Path):
     spacing = tuple(float(n) for n in options.spacing.split(','))
 
-    logger.info('scikit-image marching_cubes: spacing={} step_size={} method={}',
-                spacing, options.step_size, options.method)
+    logger.info('scikit-image marching_cubes: fwhm={} spacing={} step_size={} method={}',
+                options.fwhm, spacing, options.step_size, options.method)
 
     results = []
     with ThreadPoolExecutor(max_workers=len(os.sched_getaffinity(0))) as pool:
         mapper = PathMapper(inputdir, outputdir, glob=options.pattern, suffix='.obj')
         for mnc, obj in mapper:
-            results.append(pool.submit(mcubes, mnc, obj, options.level, spacing, options.step_size, options.method))
+            results.append(pool.submit(mcubes, mnc, obj,
+                                       options.fwhm, options.level, spacing,
+                                       options.step_size, options.method))
 
     for future in results:
-        future.exception()
+        e = future.exception()
+        if e is not None:
+            raise e
 
 
 if __name__ == '__main__':