fixed a dilation and erosion bug (and added some obscure and undocume…

…nted command-line arguments)

bin/filter_mrc/handlers.cpp

@@ -47,7 +47,7 @@ HandleDilation(const Settings &settings,
                tomo_in.aaafI,
                tomo_out.aaafI,
                mask.aaafI,
-               false,
+               settings.thickness_morpholoby_soft_penalty,
                &cerr);
 }
 
@@ -64,7 +64,7 @@ HandleErosion(const Settings &settings,
               tomo_in.aaafI,
               tomo_out.aaafI,
               mask.aaafI,
-              false,
+              settings.thickness_morpholoby_soft_penalty,
               &cerr);
 }
 
@@ -81,7 +81,7 @@ HandleOpening(const Settings &settings,
              tomo_in.aaafI,
              tomo_out.aaafI,
              mask.aaafI,
-             false,
+             settings.thickness_morpholoby_soft_penalty,
              &cerr);
 }
 
@@ -98,7 +98,7 @@ HandleClosing(const Settings &settings,
              tomo_in.aaafI,
              tomo_out.aaafI,
              mask.aaafI,
-             false,
+             settings.thickness_morpholoby_soft_penalty,
              &cerr);
 }
 
@@ -115,7 +115,7 @@ HandleTopHatWhite(const Settings &settings,
                     tomo_in.aaafI,
                     tomo_out.aaafI,
                     mask.aaafI,
-                    false,
+                    settings.thickness_morpholoby_soft_penalty,
                     &cerr);
 }
 
@@ -132,7 +132,7 @@ HandleTopHatBlack(const Settings &settings,
                     tomo_in.aaafI,
                     tomo_out.aaafI,
                     mask.aaafI,
-                    false,
+                    settings.thickness_morpholoby_soft_penalty,
                     &cerr);
 }
 
@@ -151,7 +151,6 @@ HandleMedian(const Settings &settings,
                tomo_in.aaafI,
                tomo_out.aaafI,
                mask.aaafI,
-               false,
                &cerr);
 }
 #endif

bin/filter_mrc/settings.cpp

@@ -10,7 +10,6 @@ using namespace std;
 #ifndef DISABLE_OPENMP
 #include <omp.h>       // (OpenMP-specific)
 #endif
-
 #include "err.hpp"  // defines InputErr exception class
 #include "file_io.hpp"
 #include "settings.hpp"
@@ -44,8 +43,9 @@ Settings::Settings() {
   save_intermediate_fname_base = "";
   load_intermediate_fname_base = "";
 
-  float median_radius = 0.0;
-  float thickness_morphology = 0.0;
+  median_radius = 0.0;
+  thickness_morphology = 0.0;
+  thickness_morpholoby_soft_penalty = std::numeric_limits<float>::infinity();
 
   // The code is a little bit confusing because I tried to cram as much
   // functionality into the smallest number of parameters possible:
@@ -677,6 +677,50 @@ Settings::ParseArgs(vector<string>& vArgs)
 
 
 
+    else if ((vArgs[i] == "-dilation-binary-soft") ||
+             (vArgs[i] == "-dilate-binary-soft"))
+    {
+      try {
+        if ((i+2 >= vArgs.size()) ||
+            (vArgs[i+1] == "") || (vArgs[i+1][0] == '-') ||
+            (vArgs[i+2] == "") || (vArgs[i+2][0] == '-'))
+          throw invalid_argument("");
+        thickness_morphology = stof(vArgs[i+1]);
+        //thickness_morphology_soft_max = stof(vArgs[i+2]);
+        thickness_morpholoby_soft_penalty = stof(vArgs[i+2]);
+        filter_type = DILATION;
+      }
+      catch (invalid_argument& exc) {
+        throw InputErr("Error: The " + vArgs[i] +
+                       " argument must be followed by nonnegative numbers\n");
+      }
+      num_arguments_deleted = 3;
+    } //(vArgs[i] == "-dilation-binary-soft")
+
+
+
+    else if ((vArgs[i] == "-erosion-binary-soft") ||
+             (vArgs[i] == "-erode-binary-soft"))
+    {
+      try {
+        if ((i+2 >= vArgs.size()) ||
+            (vArgs[i+1] == "") || (vArgs[i+1][0] == '-') ||
+            (vArgs[i+2] == "") || (vArgs[i+2][0] == '-'))
+          throw invalid_argument("");
+        thickness_morphology = stof(vArgs[i+1]);
+        //thickness_morphology_soft_max = stof(vArgs[i+2]);
+        thickness_morpholoby_soft_penalty = stof(vArgs[i+2]);
+        filter_type = EROSION;
+      }
+      catch (invalid_argument& exc) {
+        throw InputErr("Error: The " + vArgs[i] +
+                       " argument must be followed by nonnegative numbers\n");
+      }
+      num_arguments_deleted = 3;
+    } //(vArgs[i] == "-erosion-binary-soft")
+
+
+
     else if ((vArgs[i] == "-dilation-gauss") ||
              (vArgs[i] == "-dilate-gauss") ||
              (vArgs[i] == "-erosion-gauss") ||
@@ -694,16 +738,18 @@ Settings::ParseArgs(vector<string>& vArgs)
                        " argument must be followed by a nonnegative number\n");
       }
       filter_type = GAUSS;
+      // First, specify the "thickness" of the opening or closing operation.
+      // This is the same as the "blur_distance" used in Gaussian blurring.
       width_a[0] = blur_distance;
       width_a[1] = width_a[0];
       width_a[2] = width_a[0];
+      // Then append a post-processing threshold filter:
+      use_intensity_map = true;
       if (vArgs[i] == "-dilate-gauss")
         in_threshold_01_a = 0.1572992070502851; // ≈ 1-erf(1)
       else if (vArgs[i] == "-erode-gauss")
         in_threshold_01_a = 0.8427007929497149;  // ≈ erf(1)
-
       in_threshold_01_b = in_threshold_01_a;
-      use_intensity_map = true;
 
       num_arguments_deleted = 2;
     } // if ((vArgs[i] == "-erode-gauss") || (vArgs[i] == "-dilate-gauss"))

bin/filter_mrc/settings.hpp

@@ -104,6 +104,8 @@ class Settings {
 
   float median_radius;        // radius used for the median filter
   float thickness_morphology; // distance parameter for morphological operations
+  float thickness_morpholoby_soft_max;
+  float thickness_morpholoby_soft_penalty;
 
   // ---- parameters for "difference of generalized gaussians" filters ----
   //

doc/doc_filter_mrc.md

@@ -366,8 +366,10 @@ Other, simpler filters and operations are also provided
 
 
 The
-["**-dilation**"](#-dilate-thickness) and
-["**-erosion**"](#-erosion-thickness)
+["**-dilation**"](#-dilate-thickness),
+["**-erosion**"](#-erosion-thickness),
+["**-dilation-gauss**"](#-dilate-thickness), and
+["**-erosion-gauss**"](#-erosion-thickness)
 filters are used for
 [enlarging](https://en.wikipedia.org/wiki/Dilation_(morphology)) or
 [shrinking](https://en.wikipedia.org/wiki/Erosion_(morphology))
@@ -971,15 +973,15 @@ The brightness value of each voxel in a mask image is typically either 0 or 1
 *(depending on whether the voxel is supposed to be
 ignored or considered during subsequent calculations)*.
 
-
 ***WARNING:*** Dilation and erosion are slow for large *thickness* parameters.
 These filters have not been optimized for speed.
-If you have a binary image (images whose voxels have brightnesses of
-either 0 or 1), you can try using the fast Gaussian dilation and erosion
-method, using the
+
+If you have a binary image (images whose voxels have brightnesses 
+in the range from 0 to 1), you can try using the
 [-dilation-gauss](#-dilation-gauss-thickness) and
 [-erosion-gauss](#-erosion-gauss-thickness)
 arguments.
+(These fast filters are insensitive to small features in the image.)
 *(For additional morphological filter operations, see
 [binary_dilate](https://scikit-image.org/docs/dev/api/skimage.morphology.html?highlight=skeletonize#skimage.morphology.binary_dilation)
 [binary_erode](https://scikit-image.org/docs/dev/api/skimage.morphology.html?highlight=skeletonize#skimage.morphology.binary_erosion)

lib/visfd/filter3d.hpp

@@ -1639,7 +1639,6 @@ MedianSphere(Scalar radius,              //!< radius of the sphere (in voxels)
              Scalar const *const *const *aaafSource, //!< source image array
              Scalar ***aaafDest,             //!< filter results stored here
              Scalar const *const *const *aaafMask = nullptr,   //!< optional: ignore voxel ix,iy,iz if aaafMask!=nullptr and aaafMask[iz][iy][ix]==0
-             bool smooth_boundary = false, //!< attempt to produce a rounder smoother structure factor that varies between 0 and -1 at its boundary voxels
              ostream *pReportProgress = nullptr //!< report progress to the user
              )
 {

lib/visfd/morphology.hpp

@@ -344,25 +344,23 @@ DilateSphere(Scalar radius,              //!< radius of the sphere (in voxels)
              Scalar const *const *const *aaafSource, //!< source image array
              Scalar ***aaafDest,         //!< filter results stored here
              Scalar const *const *const *aaafMask = nullptr,   //!< optional: ignore voxel ix,iy,iz if aaafMask!=nullptr and aaafMask[iz][iy][ix]==0
-             bool smooth_boundary = false, //!< attempt to produce a rounder smoother structure factor that varies between 0 and -1 at its boundary voxels
+             Scalar smooth_boundary_penalty = std::numeric_limits<Scalar>::infinity(), //!< structure factor b varies between 0 and this value
              ostream *pReportProgress = nullptr //!< report progress to the user
              )
 {
   // structure_factor = a vector of (ix,iy,iz,b) values, one for each voxel
   vector<tuple<int,int,int,Scalar> > structure_factor;
-
   int Ri = ceil(radius);
   for (int iz = -Ri; iz <= Ri; iz++) {
     for (int iy = -Ri; iy <= Ri; iy++) {
       for (int ix = -Ri; ix <= Ri; ix++) {
         bool add_this_voxel = false;
         Scalar b = 0.0;
-        if (! smooth_boundary) {
+        if (smooth_boundary_penalty == std::numeric_limits<Scalar>::infinity())
+        {
           Scalar r = sqrt(SQR(ix) + SQR(iy) + SQR(iz));
-          if (r <= radius) {
+          if (r <= radius)
             add_this_voxel = true;
-            b = 1.0;
-          }
         }
         else {
           // Calculate the distance of all of the corners to the origin.  The
@@ -390,18 +388,15 @@ DilateSphere(Scalar radius,              //!< radius of the sphere (in voxels)
               }
             }
           }
-          if (r_max < radius) {
+          if (r_max < radius)
             add_this_voxel = true;
-            b = 1.0;
-          }
-          else if (r_min > radius) {
+          else if (r_min > radius)
             add_this_voxel = false;
-            b = 0.0;
-          }
           else {
             add_this_voxel = true;
             b = (r_max - radius) / (r_max - r_min);
             b = -b;
+            b *= smooth_boundary_penalty;
           }
         } // if (smooth_boundary)
         if (add_this_voxel)
@@ -434,25 +429,23 @@ ErodeSphere(Scalar radius,              //!< radius of the sphere (in voxels)
             Scalar const *const *const *aaafSource, //!< source image array
             Scalar ***aaafDest,         //!< filter results stored here
             Scalar const *const *const *aaafMask = nullptr,   //!< optional: ignore voxel ix,iy,iz if aaafMask!=nullptr and aaafMask[iz][iy][ix]==0
-            bool smooth_boundary = false, //!< attempt to produce a rounder smoother structure factor that varies between 0 and -1 at its boundary voxels
+            Scalar smooth_boundary_penalty = std::numeric_limits<Scalar>::infinity(), //!< structure factor b varies between 0 and this value
             ostream *pReportProgress = nullptr //!< report progress to the user
             )
 {
   // structure_factor = a vector of (ix,iy,iz,b) values, one for each voxel
   vector<tuple<int,int,int,Scalar> > structure_factor;
-
   int Ri = ceil(radius);
   for (int iz = -Ri; iz <= Ri; iz++) {
     for (int iy = -Ri; iy <= Ri; iy++) {
       for (int ix = -Ri; ix <= Ri; ix++) {
         bool add_this_voxel = false;
         Scalar b = 0.0;
-        if (! smooth_boundary) {
+        if (smooth_boundary_penalty == std::numeric_limits<Scalar>::infinity())
+        {
           Scalar r = sqrt(SQR(ix) + SQR(iy) + SQR(iz));
-          if (r <= radius) {
+          if (r <= radius)
             add_this_voxel = true;
-            b = 1.0;
-          }
         }
         else {
           // Calculate the distance of all of the corners to the origin.  The
@@ -480,17 +473,15 @@ ErodeSphere(Scalar radius,              //!< radius of the sphere (in voxels)
               }
             }
           }
-          if (r_max < radius) {
+          if (r_max < radius)
             add_this_voxel = true;
-            b = 1.0;
-          }
-          else if (r_min > radius) {
+          else if (r_min > radius)
             add_this_voxel = false;
-            b = 0.0;
-          }
           else {
             add_this_voxel = true;
             b = (r_max - radius) / (r_max - r_min);
+            b = -b;
+            b *= smooth_boundary_penalty;
           }
         } // if (smooth_boundary)
         if (add_this_voxel)
@@ -523,7 +514,7 @@ OpenSphere(Scalar radius,             //!< radius of the sphere (in voxels)
            Scalar const *const *const *aaafSource, //!< source image array
            Scalar ***aaafDest,        //!< filter results stored here
            Scalar const *const *const *aaafMask = nullptr,   //!< optional: ignore voxel ix,iy,iz if aaafMask!=nullptr and aaafMask[iz][iy][ix]==0
-           bool smooth_boundary = false, //!< attempt to produce a round and smooth structure factor that varies between 0 and -1 at its boundary voxels
+           Scalar smooth_boundary_penalty = std::numeric_limits<Scalar>::infinity(), //!< structure factor b varies between 0 and this value
            ostream *pReportProgress = nullptr //!< report progress to the user
            )
 {
@@ -536,15 +527,15 @@ OpenSphere(Scalar radius,             //!< radius of the sphere (in voxels)
               aaafSource,
               aaafTmp,  //<--save the resulting image in aaafTmp
               aaafMask,
-              smooth_boundary,
+              smooth_boundary_penalty,
               pReportProgress);
 
   DilateSphere(radius,
                image_size,
                aaafTmp,  //<--use aaafTmp as the source image
                aaafDest, //<--save the resulting image in aaafDest
                aaafMask,
-               smooth_boundary,
+               smooth_boundary_penalty,
                pReportProgress);
 
   Dealloc3D(aaafTmp);
@@ -563,7 +554,7 @@ CloseSphere(Scalar radius,            //!< radius of the sphere (in voxels)
             Scalar const *const *const *aaafSource, //!< source image array
             Scalar ***aaafDest,         //!< filter results stored here
             Scalar const *const *const *aaafMask = nullptr,   //!< optional: ignore voxel ix,iy,iz if aaafMask!=nullptr and aaafMask[iz][iy][ix]==0
-            bool smooth_boundary = false, //!< attempt to produce a round and smooth structure factor that varies between 0 and -1 at its boundary voxels
+            Scalar smooth_boundary_penalty = std::numeric_limits<Scalar>::infinity(), //!< structure factor b varies between 0 and this value
             ostream *pReportProgress = nullptr //!< report progress to the user
             )
 {
@@ -576,15 +567,15 @@ CloseSphere(Scalar radius,            //!< radius of the sphere (in voxels)
                aaafSource,
                aaafTmp,  //<--save the resulting image in aaafTmp
                aaafMask,
-               smooth_boundary,
+               smooth_boundary_penalty,
                pReportProgress);
 
   ErodeSphere(radius,
               image_size,
               aaafTmp,  //<--use aaafTmp as the source image
               aaafDest, //<--save the resulting image in aaafDest
               aaafMask,
-              smooth_boundary,
+              smooth_boundary_penalty,
               pReportProgress);
 
   Dealloc3D(aaafTmp);
@@ -603,7 +594,7 @@ WhiteTopHatSphere(Scalar radius,             //!< radius of the sphere (in voxel
                   Scalar const *const *const *aaafSource, //!<source image array
                   Scalar ***aaafDest,         //!< filter results stored here
                   Scalar const *const *const *aaafMask = nullptr,   //!< optional: ignore voxel ix,iy,iz if aaafMask!=nullptr and aaafMask[iz][iy][ix]==0
-                  bool smooth_boundary = false, //!< attempt to produce a round and smooth structure factor that varies between 0 and -1 at its boundary voxels
+                  Scalar smooth_boundary_penalty = std::numeric_limits<Scalar>::infinity(), //!< structure factor b varies between 0 and this value
                   ostream *pReportProgress = nullptr //!< report progress to the user
                   )
 {
@@ -617,7 +608,7 @@ WhiteTopHatSphere(Scalar radius,             //!< radius of the sphere (in voxel
              aaafSource,
              aaafTmp,  //<--save the resulting image in aaafTmp
              aaafMask,
-             smooth_boundary,
+             smooth_boundary_penalty,
              pReportProgress);
 
   // Subtract it from the original image brightness
@@ -642,7 +633,7 @@ BlackTopHatSphere(Scalar radius,             //!< radius of the sphere (in voxel
                   Scalar const *const *const *aaafSource,//!< source image array
                   Scalar ***aaafDest,         //!< filter results stored here
                   Scalar const *const *const *aaafMask = nullptr,   //!< optional: ignore voxel ix,iy,iz if aaafMask!=nullptr and aaafMask[iz][iy][ix]==0
-                  bool smooth_boundary = false, //!< attempt to produce a round and smooth structure factor that varies between 0 and -1 at its boundary voxels
+                  Scalar smooth_boundary_penalty = std::numeric_limits<Scalar>::infinity(), //!< structure factor b varies between 0 and this value
                   ostream *pReportProgress = nullptr //!< report progress to the user
                   )
 {
@@ -656,7 +647,7 @@ BlackTopHatSphere(Scalar radius,             //!< radius of the sphere (in voxel
               aaafSource,
               aaafTmp,  //<--save the resulting image in aaafTmp
               aaafMask,
-              smooth_boundary,
+              smooth_boundary_penalty,
               pReportProgress);
 
   // Subtract the original image brightness