fixed bugs in "-thresh", "-blur-expand", "-blur-contract", "-mask-rec…

…t", "-mask-rect-subtract" for the "filter_mrc" program. One or two bugs in the VISFD library (in "draw.hpp") were also corrected.

bin/filter_mrc/filter_mrc.cpp

@@ -29,7 +29,7 @@ using namespace std;
 
 
 string g_program_name("filter_mrc");
-string g_version_string("0.27.1");
+string g_version_string("0.27.2");
 string g_date_string("2021-7-05");
 
 

bin/filter_mrc/settings.cpp

@@ -159,7 +159,7 @@ Settings::Settings() {
   in_threshold_10_a = 0.0;
   in_threshold_10_b = 0.0;
   out_thresh2_use_clipping = false;
-  out_thresh2_use_clipping_sigma = true;
+  out_thresh2_use_clipping_sigma = false;
   out_thresh_a_value = 0.0;
   out_thresh_b_value = 1.0;
   rescale_min_max_in = false;
@@ -543,9 +543,9 @@ Settings::ParseArgs(vector<string>& vArgs)
       width_a[1] = width_a[0];
       width_a[2] = width_a[0];
       if ((vArgs[i] == "-blur-expand") || (vArgs[i] == "-expand"))
-        in_threshold_01_a = 0.07864960352514255; // ≈ (1-erf(1))/2
+        in_threshold_01_a = 0.1572992070502851; // ≈ 1-erf(1)
       else if ((vArgs[i] == "-blur-contract") || (vArgs[i] == "-contract"))
-        in_threshold_01_a = 0.9213503964748575;  // ≈ (1+erf(1))/2
+        in_threshold_01_a = 0.8427007929497149;  // ≈ erf(1)
 
       in_threshold_01_b = in_threshold_01_a;
       use_intensity_map = true;

doc/doc_filter_mrc.md

@@ -202,7 +202,6 @@ of the surface by approximately 40 Angstroms inward, creating a new file
 filter_mrc -i segmented.rec -o segmented_interior.rec -blur-contract 40
 ```
 *(See [-blur-contract](#--blur-contract-distance) for details.)*
-*(WARNING: The -blur-contract feature has not yet been tested. -Andrew 2021-7-05)*
 
 
 **Note:**
@@ -1507,8 +1506,6 @@ of the Gaussian independently in the x,y,z directions:
 
 ### -bin binsize
 
-***WARNING: THIS FEATURE IS EXPERIMENTAL.  PLEASE REPORT BUGS.  2021-6-21***
-
 Reduce the resolution of the image in each direction by a factor of *binsize*.
 (*binsize* must be a positive integer.)
 The new image will be smaller in each direction by a factof of *binsize*.
@@ -2504,7 +2501,6 @@ This provides a way to see where the mask region is located.
 ### -blur-contract distance
 ### -blur-expand distance
 
-***WARNING: This feature has not yet been tested. -Andrew 2021-7-05***
 
 The **-blur-contract** and **-blur-expand** arguments are useful
 for modifying the size of an existing binary image (eg. mask),
@@ -2538,9 +2534,9 @@ Consequently, you cannot combine these arguments with either the *-gauss*
 or *-thresh* arguments, because:
 
 - "-blur-expand *distance*" is equivalent to
-"-gauss *distance* -thresh 0.0786496", *(where 0.0786496 ≈ (1-erf(1))/2)*
+"-gauss *distance* -thresh 0.157299", *(where 0.157299 ≈ 1-erf(1))*
 - "-blur-contract *distance*" is equivalent to
-"-gauss *distance* -thresh 0.9213504", *(where 0.9213504 ≈ (1+erf(1))/2)*
+"-gauss *distance* -thresh 0.842701", *(where 0.842701 ≈ erf(1))*
 
 The motivation for this strategy is outlined below:
 If you start with a binary image (consisting of 1s and 0s),

lib/visfd/draw.hpp

@@ -96,6 +96,44 @@ DrawRegions(int const image_size[3], //!< image size
             bool negative_means_subtract = false  //!< should we interpret negative brightness as set subtraction?
             )
 {
+
+  // Special case:  Should the background voxels be initialized with 1s?
+  // If the first region has a negative value and negative_means_subtract==true,
+  // then it means that the caller wants to remove (subtract) those voxels from
+  // the set of non-zero background voxels, by setting their brightness to zero.
+  // That doesn't make sense if all of the voxels in the image
+  // already have brightness zero.  So I interpret this special case as a
+  // request to subtract these voxels from a binary image which is filled with
+  // non-zero values.  Since a binary image only contains 1s or 0s, I will
+  // initialize the image full of 1s.  Then in the next step, we will remove
+  // the voxels in the first region (vRegions[0]) from this background of 1s.
+  if (negative_means_subtract &&
+      (vRegions.size() > 0) && (vRegions[0].value < 0)) {
+    bool all_zero = true;
+    for (int iz = 0; iz < image_size[2] && all_zero; iz++) {
+      for (int iy = 0; iy < image_size[1] && all_zero; iy++) {
+        for (int ix = 0; ix < image_size[0] && all_zero; ix++) {
+          if (aaafMask && (aaafMask[iz][iy][ix] == 0.0))
+            continue;
+          if (aaafDest[iz][iy][ix] != 0.0)
+            all_zero = false;
+        }
+      }
+    }
+    if (all_zero) {
+      for (int iz = 0; iz < image_size[2] && all_zero; iz++) {
+        for (int iy = 0; iy < image_size[1] && all_zero; iy++) {
+          for (int ix = 0; ix < image_size[0] && all_zero; ix++) {
+            if (aaafMask && (aaafMask[iz][iy][ix] == 0.0))
+              continue;
+            aaafDest[iz][iy][ix] = 1.0;
+          }
+        }
+      }
+    }
+  } // if (negative_means_subtract && (vRegions[0].value < 0))
+
+
   for (int i=0; i < vRegions.size(); i++) {
 
     switch (vRegions[i].type) {
@@ -152,13 +190,13 @@ DrawRegions(int const image_size[3], //!< image size
         Scalar izmax = floor(vRegions[i].data.rect.zmax + 0.5);
 
         for (int iz=std::max<float>(izmin, 0);
-             iz < std::min<float>(izmax, image_size[2]);
+             iz <= std::min<float>(izmax, image_size[2]-1);
              iz++) {
           for (int iy=std::max<float>(iymin, 0);
-               iy < std::min<float>(iymax, image_size[1]);
+               iy <= std::min<float>(iymax, image_size[1]-1);
                iy++) {
             for (int ix=std::max<float>(ixmin, 0);
-                 ix < std::min<float>(ixmax, image_size[0]);
+                 ix <= std::min<float>(ixmax, image_size[0]-1);
                  ix++) {
               if (aaafMask && (aaafMask[iz][iy][ix] == 0.0))
                 continue;