Address #2371

Add SignedMaurerDistanceMapImageFilter as default distance transform method

This change makes SignedMaurerDistanceMapImageFilter the default method for computing distance transforms, replacing ReinitializeLevelSetImageFilter. The new default is faster and more robust, working directly on segmentations without requiring anti-aliasing.

  Changes:
  - Add DistanceTransformType enum to Image class with SignedMaurer and
    ReinitializeLevelSet options
  - Update Image::computeDT() to accept method parameter with SignedMaurer
    as default
  - Expose DistanceTransformType enum in Python bindings
  - Add --method option to CLI compute-dt command (choices: "maurer",
    "levelset")
  - Update existing tests to explicitly use ReinitializeLevelSet method
    to match existing ground truth files
  - Add new tests to verify SignedMaurer method and default behavior
  - Update GroomTests baseline to use new default method
  - Maintain backward compatibility by keeping legacy method available

  The legacy ReinitializeLevelSetImageFilter remains available for cases
  where the old behavior is needed.

Author: akenmorris

Applications/shapeworks/ImageCommands.cpp

@@ -548,6 +548,8 @@ void ComputeDT::buildParser()
   parser.prog(prog).description(desc);
 
   parser.add_option("--isovalue").action("store").type("double").set_default(0.0).help("Level set value that defines the interface between foreground and background [default: %default].");
+  std::list<std::string> methods{"maurer", "levelset"};
+  parser.add_option("--method").action("store").type("choice").choices(methods.begin(), methods.end()).set_default("maurer").help("Distance transform method: 'maurer' (SignedMaurerDistanceMapImageFilter - faster, more robust, default) or 'levelset' (ReinitializeLevelSetImageFilter - legacy) [default: %default].");
 
   Command::buildParser();
 }
@@ -561,8 +563,10 @@ bool ComputeDT::execute(const optparse::Values &options, SharedCommandData &shar
   }
 
   double isoValue = static_cast<double>(options.get("isovalue"));
+  std::string method_str(static_cast<std::string>(options.get("method")));
+  Image::DistanceTransformType method = (method_str == "levelset") ? Image::ReinitializeLevelSet : Image::SignedMaurer;
 
-  sharedData.image.computeDT(isoValue);
+  sharedData.image.computeDT(isoValue, method);
   return true;
 }
 

Libs/Image/Image.cpp

@@ -24,6 +24,7 @@
 #include <itkOrientImageFilter.h>
 #include <itkRegionOfInterestImageFilter.h>
 #include <itkReinitializeLevelSetImageFilter.h>
+#include <itkSignedMaurerDistanceMapImageFilter.h>
 #include <itkRelabelComponentImageFilter.h>
 #include <itkResampleImageFilter.h>
 #include <itkScalableAffineTransform.h>
@@ -628,15 +629,28 @@ Image& Image::binarize(PixelType minVal, PixelType maxVal, PixelType innerVal, P
   return *this;
 }
 
-Image& Image::computeDT(PixelType isoValue) {
-  using FilterType = itk::ReinitializeLevelSetImageFilter<ImageType>;
-  FilterType::Pointer filter = FilterType::New();
-
-  filter->SetInput(this->itk_image_);
-  filter->NarrowBandingOff();
-  filter->SetLevelSetValue(isoValue);
-  filter->Update();
-  this->itk_image_ = filter->GetOutput();
+Image& Image::computeDT(PixelType isoValue, DistanceTransformType method) {
+  if (method == SignedMaurer) {
+    using FilterType = itk::SignedMaurerDistanceMapImageFilter<ImageType, ImageType>;
+    FilterType::Pointer filter = FilterType::New();
+
+    filter->SetInput(this->itk_image_);
+    filter->SetBackgroundValue(isoValue);
+    filter->SetInsideIsPositive(false);
+    filter->SetSquaredDistance(false);
+    filter->SetUseImageSpacing(true);
+    filter->Update();
+    this->itk_image_ = filter->GetOutput();
+  } else {  // ReinitializeLevelSet
+    using FilterType = itk::ReinitializeLevelSetImageFilter<ImageType>;
+    FilterType::Pointer filter = FilterType::New();
+
+    filter->SetInput(this->itk_image_);
+    filter->NarrowBandingOff();
+    filter->SetLevelSetValue(isoValue);
+    filter->Update();
+    this->itk_image_ = filter->GetOutput();
+  }
 
   return *this;
 }

Libs/Image/Image.h

@@ -28,13 +28,14 @@ class Mesh;
 class Image {
  public:
   enum InterpolationType { Linear, NearestNeighbor };
+  enum DistanceTransformType { SignedMaurer, ReinitializeLevelSet };
 
   using PixelType = float;
   using ImageType = itk::Image<PixelType, 3>;
   using StatsPtr = itk::StatisticsImageFilter<ImageType>::Pointer;
   using ImageIterator = itk::ImageRegionIterator<ImageType>;
   using InterpolatorType = itk::LinearInterpolateImageFunction<ImageType>;
-
+  
   // constructors and assignment operators //
   explicit Image(const Dims dims);
   explicit Image(const std::string& pathname) : itk_image_(read(pathname)) {}
@@ -161,7 +162,7 @@ class Image {
                   PixelType innerVal = 1.0, PixelType outerVal = 0.0);
 
   /// computes distance transform volume from a (preferably antialiased) binary image using the specified isovalue
-  Image& computeDT(PixelType isoValue = 0.0);
+  Image& computeDT(PixelType isoValue = 0.0, DistanceTransformType method = SignedMaurer);
 
   /// denoises an image using curvature driven flow using curvature flow image filter
   Image& applyCurvatureFilter(unsigned iterations = 10);

Libs/Python/ShapeworksPython.cpp

@@ -94,6 +94,12 @@ PYBIND11_MODULE(shapeworks_py, m) {
       .value("NearestNeighbor", Image::InterpolationType::NearestNeighbor)
       .export_values();
 
+  // Image::DistanceTransformType
+  py::enum_<Image::DistanceTransformType>(m, "DistanceTransformType")
+      .value("SignedMaurer", Image::DistanceTransformType::SignedMaurer)
+      .value("ReinitializeLevelSet", Image::DistanceTransformType::ReinitializeLevelSet)
+      .export_values();
+
   m.def(
       "mean", [](py::array& field) { return mean(pyToArr(field, false /*take_ownership*/)); },
       "incrementally compute (single-component) mean of field");
@@ -323,7 +329,8 @@ PYBIND11_MODULE(shapeworks_py, m) {
            "outerVal"_a = 0.0)
 
       .def("computeDT", &Image::computeDT,
-           "computes signed distance transform volume from an image at the specified isovalue", "isovalue"_a = 0.0)
+           "computes signed distance transform volume from an image at the specified isovalue using the specified method (SignedMaurer is default, faster and more robust)",
+           "isovalue"_a = 0.0, "method"_a = Image::DistanceTransformType::SignedMaurer)
 
       .def("applyCurvatureFilter", &Image::applyCurvatureFilter,
            "denoises an image using curvature driven flow using curvature flow image filter", "iterations"_a = 10)

Testing/ImageTests/ImageTests.cpp

@@ -280,20 +280,42 @@ TEST(ImageTests, binarizeTest2) {
 
 TEST(ImageTests, computedtTest1) {
   Image image(std::string(TEST_DATA_DIR) + "/1x2x2.nrrd");
-  image.computeDT();
+  // Use legacy method to match existing ground truth
+  image.computeDT(0.0, Image::ReinitializeLevelSet);
   Image ground_truth(std::string(TEST_DATA_DIR) + "/computedt1.nrrd");
 
   ASSERT_TRUE(image == ground_truth);
 }
 
 TEST(ImageTests, computedtTest2) {
   Image image(std::string(TEST_DATA_DIR) + "/1x2x2.nrrd");
-  image.computeDT(1.0);
+  // Use legacy method to match existing ground truth
+  image.computeDT(1.0, Image::ReinitializeLevelSet);
   Image ground_truth(std::string(TEST_DATA_DIR) + "/computedt2.nrrd");
 
   ASSERT_TRUE(image == ground_truth);
 }
 
+// Test the new default SignedMaurer method
+TEST(ImageTests, computedtTest3) {
+  Image image(std::string(TEST_DATA_DIR) + "/1x2x2.nrrd");
+  // Use new default SignedMaurer method - should work without anti-aliasing
+  image.computeDT(0.0, Image::SignedMaurer);
+
+  // Should produce a distance transform (negative inside, positive outside)
+  ASSERT_TRUE(image.isDistanceTransform());
+}
+
+TEST(ImageTests, computedtTest4) {
+  Image image(std::string(TEST_DATA_DIR) + "/1x2x2.nrrd");
+  // Test that default is SignedMaurer
+  Image image_default(std::string(TEST_DATA_DIR) + "/1x2x2.nrrd");
+  image.computeDT(0.0, Image::SignedMaurer);
+  image_default.computeDT(0.0);  // Should use SignedMaurer by default
+
+  ASSERT_TRUE(image == image_default);
+}
+
 TEST(ImageTests, curvatureTest) {
   Image image(std::string(TEST_DATA_DIR) + "/1x2x2.nrrd");
   image.applyCurvatureFilter();

Testing/PythonTests/computedt.py

@@ -6,7 +6,8 @@
 
 def computedtTest1():
   img = Image(os.environ["DATA"] + "/1x2x2.nrrd")
-  img.computeDT()
+  # Use legacy method to match existing ground truth
+  img.computeDT(0.0, DistanceTransformType.ReinitializeLevelSet)
 
   compareImg = Image(os.environ["DATA"] + "/computedt1.nrrd")
 
@@ -16,7 +17,8 @@ def computedtTest1():
 
 def computedtTest2():
   img = Image(os.environ["DATA"] + "/1x2x2.nrrd")
-  img.computeDT(1.0)
+  # Use legacy method to match existing ground truth
+  img.computeDT(1.0, DistanceTransformType.ReinitializeLevelSet)
 
   compareImg = Image(os.environ["DATA"] + "/computedt2.nrrd")
 
@@ -26,7 +28,8 @@ def computedtTest2():
 
 def computedtTest3():
   img = Image(os.environ["DATA"] + "/1x2x2.nrrd")
-  img.computeDT(-1.0)
+  # Use legacy method to match existing ground truth
+  img.computeDT(-1.0, DistanceTransformType.ReinitializeLevelSet)
 
   compareImg = Image(os.environ["DATA"] + "/computedt3.nrrd")
 

Testing/data/optimize/shared/spheres/sphere10_DT.nrrd

@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:849b13d964d22c6835946981fcbe1e9bd1e26625f7356075ba8f19660c73f65f
-size 62838
+oid sha256:1f6f60eca442bc9341779f907019897fa6eecd4880f71a5ed06c57d9cf915d52
+size 220399

Testing/data/optimize/shared/spheres/sphere10_DT_baseline.nrrd

@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:adddca3e41326bdb30cdaba8f0a094d46c544013047b98e9083430ce2db67465
-size 169723
+oid sha256:b7a3c296a7c4bc40c75501a80147b78a962188ad7cce21a194523a7b0b787ab3
+size 165912