added soft versions of the dilation and erosion morphological filters…

…. (You can access them in filter_mrc using the "-dilation-binary-soft" and "-erosion-binary-soft" arguments.)

bin/filter_mrc/filter_mrc.cpp

@@ -29,7 +29,7 @@ using namespace std;
 
 
 string g_program_name("filter_mrc");
-string g_version_string("0.29.4");
+string g_version_string("0.29.5");
 string g_date_string("2021-7-17");
 
 
@@ -269,7 +269,9 @@ int main(int argc, char **argv) {
     // calculations considerably, so I just assume cube-shaped voxels:
     //assert((voxel_width[0] == voxel_width[1]) &&
     //       (voxel_width[1] == voxel_width[2]));
-    settings.thickness_morphology /= voxel_width[0];
+
+    settings.morphology_r /= voxel_width[0];
+    settings.morphology_rmax /= voxel_width[0];
 
 
     // Mote: max_distance_to_feature is specified by the user in units of voxels

bin/filter_mrc/handlers.cpp

@@ -42,12 +42,13 @@ HandleDilation(const Settings &settings,
                MrcSimple &mask,
                float voxel_width[3])
 {
-  DilateSphere(settings.thickness_morphology,
+  DilateSphere(settings.morphology_r,
                tomo_in.header.nvoxels,
                tomo_in.aaafI,
                tomo_out.aaafI,
                mask.aaafI,
-               settings.thickness_morpholoby_soft_penalty,
+               settings.morphology_rmax,
+               settings.morphology_bmax,
                &cerr);
 }
 
@@ -59,12 +60,13 @@ HandleErosion(const Settings &settings,
               MrcSimple &mask,
               float voxel_width[3])
 {
-  ErodeSphere(settings.thickness_morphology,
+  ErodeSphere(settings.morphology_r,
               tomo_in.header.nvoxels,
               tomo_in.aaafI,
               tomo_out.aaafI,
               mask.aaafI,
-              settings.thickness_morpholoby_soft_penalty,
+              settings.morphology_rmax,
+              settings.morphology_bmax,
               &cerr);
 }
 
@@ -76,12 +78,13 @@ HandleOpening(const Settings &settings,
               MrcSimple &mask,
               float voxel_width[3])
 {
-  OpenSphere(settings.thickness_morphology,
+  OpenSphere(settings.morphology_r,
              tomo_in.header.nvoxels,
              tomo_in.aaafI,
              tomo_out.aaafI,
              mask.aaafI,
-             settings.thickness_morpholoby_soft_penalty,
+             settings.morphology_rmax,
+             settings.morphology_bmax,
              &cerr);
 }
 
@@ -93,13 +96,14 @@ HandleClosing(const Settings &settings,
               MrcSimple &mask,
               float voxel_width[3])
 {
-  CloseSphere(settings.thickness_morphology,
-             tomo_in.header.nvoxels,
-             tomo_in.aaafI,
-             tomo_out.aaafI,
-             mask.aaafI,
-             settings.thickness_morpholoby_soft_penalty,
-             &cerr);
+  CloseSphere(settings.morphology_r,
+              tomo_in.header.nvoxels,
+              tomo_in.aaafI,
+              tomo_out.aaafI,
+              mask.aaafI,
+              settings.morphology_rmax,
+              settings.morphology_bmax,
+              &cerr);
 }
 
 
@@ -110,12 +114,13 @@ HandleTopHatWhite(const Settings &settings,
                   MrcSimple &mask,
                   float voxel_width[3])
 {
-  WhiteTopHatSphere(settings.thickness_morphology,
+  WhiteTopHatSphere(settings.morphology_r,
                     tomo_in.header.nvoxels,
                     tomo_in.aaafI,
                     tomo_out.aaafI,
                     mask.aaafI,
-                    settings.thickness_morpholoby_soft_penalty,
+                    settings.morphology_rmax,
+                    settings.morphology_bmax,
                     &cerr);
 }
 
@@ -127,12 +132,13 @@ HandleTopHatBlack(const Settings &settings,
                   MrcSimple &mask,
                   float voxel_width[3])
 {
-  BlackTopHatSphere(settings.thickness_morphology,
+  BlackTopHatSphere(settings.morphology_r,
                     tomo_in.header.nvoxels,
                     tomo_in.aaafI,
                     tomo_out.aaafI,
                     mask.aaafI,
-                    settings.thickness_morpholoby_soft_penalty,
+                    settings.morphology_rmax,
+                    settings.morphology_bmax,
                     &cerr);
 }
 

bin/filter_mrc/settings.cpp

@@ -44,8 +44,9 @@ Settings::Settings() {
   load_intermediate_fname_base = "";
 
   median_radius = 0.0;
-  thickness_morphology = 0.0;
-  thickness_morpholoby_soft_penalty = std::numeric_limits<float>::infinity();
+  morphology_r = 0.0;
+  morphology_rmax = 0.0;
+  morphology_bmax = 0.0;
 
   // The code is a little bit confusing because I tried to cram as much
   // functionality into the smallest number of parameters possible:
@@ -646,7 +647,7 @@ Settings::ParseArgs(vector<string>& vArgs)
         if ((i+1 >= vArgs.size()) ||
             (vArgs[i+1] == "") || (vArgs[i+1][0] == '-'))
           throw invalid_argument("");
-        thickness_morphology = stof(vArgs[i+1]);
+        morphology_r = stof(vArgs[i+1]);
         filter_type = DILATION;
       }
       catch (invalid_argument& exc) {
@@ -665,7 +666,7 @@ Settings::ParseArgs(vector<string>& vArgs)
         if ((i+1 >= vArgs.size()) ||
             (vArgs[i+1] == "") || (vArgs[i+1][0] == '-'))
           throw invalid_argument("");
-        thickness_morphology = stof(vArgs[i+1]);
+        morphology_r = stof(vArgs[i+1]);
         filter_type = EROSION;
       }
       catch (invalid_argument& exc) {
@@ -683,18 +684,19 @@ Settings::ParseArgs(vector<string>& vArgs)
       try {
         if ((i+2 >= vArgs.size()) ||
             (vArgs[i+1] == "") || (vArgs[i+1][0] == '-') ||
-            (vArgs[i+2] == "") || (vArgs[i+2][0] == '-'))
+            (vArgs[i+2] == "") || (vArgs[i+2][0] == '-') ||
+            (vArgs[i+3] == "") || (vArgs[i+3][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]);
+        morphology_r = stof(vArgs[i+1]);
+        morphology_rmax = stof(vArgs[i+2]);
+        morphology_bmax = stof(vArgs[i+3]);
         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;
+      num_arguments_deleted = 4;
     } //(vArgs[i] == "-dilation-binary-soft")
 
 
@@ -703,20 +705,21 @@ Settings::ParseArgs(vector<string>& vArgs)
              (vArgs[i] == "-erode-binary-soft"))
     {
       try {
-        if ((i+2 >= vArgs.size()) ||
+        if ((i+3 >= vArgs.size()) ||
             (vArgs[i+1] == "") || (vArgs[i+1][0] == '-') ||
-            (vArgs[i+2] == "") || (vArgs[i+2][0] == '-'))
+            (vArgs[i+2] == "") || (vArgs[i+2][0] == '-') ||
+            (vArgs[i+3] == "") || (vArgs[i+3][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]);
+        morphology_r = stof(vArgs[i+1]);
+        morphology_rmax = stof(vArgs[i+2]);
+        morphology_bmax = stof(vArgs[i+3]);
         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;
+      num_arguments_deleted = 4;
     } //(vArgs[i] == "-erosion-binary-soft")
 
 
@@ -763,7 +766,7 @@ Settings::ParseArgs(vector<string>& vArgs)
         if ((i+1 >= vArgs.size()) ||
             (vArgs[i+1] == "") || (vArgs[i+1][0] == '-'))
           throw invalid_argument("");
-        thickness_morphology = stof(vArgs[i+1]);
+        morphology_r = stof(vArgs[i+1]);
         filter_type = OPENING;
       }
       catch (invalid_argument& exc) {
@@ -782,7 +785,7 @@ Settings::ParseArgs(vector<string>& vArgs)
         if ((i+1 >= vArgs.size()) ||
             (vArgs[i+1] == "") || (vArgs[i+1][0] == '-'))
           throw invalid_argument("");
-        thickness_morphology = stof(vArgs[i+1]);
+        morphology_r = stof(vArgs[i+1]);
         filter_type = CLOSING;
       }
       catch (invalid_argument& exc) {
@@ -800,7 +803,7 @@ Settings::ParseArgs(vector<string>& vArgs)
         if ((i+1 >= vArgs.size()) ||
             (vArgs[i+1] == "") || (vArgs[i+1][0] == '-'))
           throw invalid_argument("");
-        thickness_morphology = stof(vArgs[i+1]);
+        morphology_r = stof(vArgs[i+1]);
         filter_type = TOP_HAT_WHITE;
       }
       catch (invalid_argument& exc) {
@@ -818,7 +821,7 @@ Settings::ParseArgs(vector<string>& vArgs)
         if ((i+1 >= vArgs.size()) ||
             (vArgs[i+1] == "") || (vArgs[i+1][0] == '-'))
           throw invalid_argument("");
-        thickness_morphology = stof(vArgs[i+1]);
+        morphology_r = stof(vArgs[i+1]);
         filter_type = TOP_HAT_BLACK;
       }
       catch (invalid_argument& exc) {

bin/filter_mrc/settings.hpp

@@ -102,10 +102,10 @@ class Settings {
   string save_intermediate_fname_base = ""; //save progress of slow calculations
   string load_intermediate_fname_base = ""; //load progress from a previous run?
 
-  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;
+  float median_radius;       // radius used for the median filter
+  float morphology_r;        // distance parameter for morphological operations
+  float morphology_rmax;     // distance parameter for morphological operations
+  float morphology_bmax;     // maximum b value used in soft structure factors
 
   // ---- parameters for "difference of generalized gaussians" filters ----
   //

doc/doc_filter_mrc.md

@@ -1036,16 +1036,16 @@ the size of the large objects that you want removed from the image
 (in physical distance units, not in voxels).*
 
 
-### -erosion-gauss thickness
 ### -dilation-gauss thickness
+### -erosion-gauss thickness
 
 These filters are an alternative implementation of
 [image dilation](https://en.wikipedia.org/wiki/Dilation_(morphology)) and
 [image erosion](https://en.wikipedia.org/wiki/Erosion_(morphology))
 that uses fast Gaussian filters to blur an image followed by thresholding.
 
-*Note: These arguments use filters that are only intended for use with
-binary images (images whose voxels have brightnesses of either 0 or 1).*
+*Note: These arguments are only intended for use with binary images
+(images whose voxels have brightnesses of either 0 or 1).*
 For binary images containing large smooth regions of bright voxels,
 the **-dilation-gauss** and **-erosion-gauss** arguments
 will enlarge or shrink the size of these bright regions
@@ -1065,10 +1065,6 @@ for some morphological tasks *(such as
 [top-hat transform](https://en.wikipedia.org/wiki/Top-hat_transform))*.
 Use with caution.
 
-
-
-
-
 #### Details
 This is implemented by blurring the image, and then applying a threshold
 so that all voxels whose brightness are below a threshold are discarded.
@@ -1085,6 +1081,28 @@ expands or shrinks by a distance of *thickness*.
 (But, technically, this is only true near a flat boundary.)
 
 
+### -dilation-binary-soft r1 r2 bmax
+### -erosion-binary-soft r1 r2 bmax
+
+These arguments will perform a morphological
+[dilation](https://en.wikipedia.org/wiki/Dilation_(morphology)) and
+[erosion](https://en.wikipedia.org/wiki/Erosion_(morphology)).
+Unlike traditional dilation and erosion operations which have a hard distance
+cutoff, this version has a cutoff that varies gradually from *r1* to *r2*.
+(Unless the "-w 1" argument is used, both *r1* and *r2* are expressed in
+units of physical distance, not voxels.  *Special case: If r2=0*,
+a soft edge is used which is only 1-voxel wide.)
+
+***Note:*** These "binary-soft" arguments are only intended for use with
+*binary* images or grayscale images whose voxels have brightness
+in the range from *0* to *bmax*.  (Typically *bmax=1*.)
+
+*(Technical details: This version uses a soft spherical structure factor,
+*b(r)* that varies from *0*, when *r=r1*, up to to *-bmax*, when
+*r=r2*, and -∞, when *r>r2*.  Setting *r2=0* results in a
+spherical structure factor whose brightness varies from *0* to *-bmax*
+on a thin 1-voxel wide shell at the sphere's boundary.)*
+
 
 ### -find-minima and -find-maxima