diff --git a/bin/filter_mrc/handlers.cpp b/bin/filter_mrc/handlers.cpp index 31d1e11..2915e03 100644 --- a/bin/filter_mrc/handlers.cpp +++ b/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 diff --git a/bin/filter_mrc/settings.cpp b/bin/filter_mrc/settings.cpp index 5dab22b..6782210 100644 --- a/bin/filter_mrc/settings.cpp +++ b/bin/filter_mrc/settings.cpp @@ -10,7 +10,6 @@ using namespace std; #ifndef DISABLE_OPENMP #include // (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::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& 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& 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")) diff --git a/bin/filter_mrc/settings.hpp b/bin/filter_mrc/settings.hpp index c43c455..b04c2c3 100644 --- a/bin/filter_mrc/settings.hpp +++ b/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 ---- // diff --git a/doc/doc_filter_mrc.md b/doc/doc_filter_mrc.md index aa0e757..be424d9 100644 --- a/doc/doc_filter_mrc.md +++ b/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) diff --git a/lib/visfd/filter3d.hpp b/lib/visfd/filter3d.hpp index b333fd1..04e855e 100644 --- a/lib/visfd/filter3d.hpp +++ b/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 ) { diff --git a/lib/visfd/morphology.hpp b/lib/visfd/morphology.hpp index 6ca8582..cda42a3 100644 --- a/lib/visfd/morphology.hpp +++ b/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::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 > 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::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::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 > 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::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::infinity(), //!< structure factor b varies between 0 and this value ostream *pReportProgress = nullptr //!< report progress to the user ) { @@ -536,7 +527,7 @@ 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, @@ -544,7 +535,7 @@ OpenSphere(Scalar radius, //!< radius of the sphere (in voxels) 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::infinity(), //!< structure factor b varies between 0 and this value ostream *pReportProgress = nullptr //!< report progress to the user ) { @@ -576,7 +567,7 @@ 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, @@ -584,7 +575,7 @@ CloseSphere(Scalar radius, //!< radius of the sphere (in voxels) 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, //!::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::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 diff --git a/lib/visfd/visfd_utils.hpp b/lib/visfd/visfd_utils.hpp index 39067d8..50fbab3 100644 --- a/lib/visfd/visfd_utils.hpp +++ b/lib/visfd/visfd_utils.hpp @@ -30,6 +30,7 @@ namespace visfd { template constexpr Scalar SQR(Scalar x) { return x*x; } + template constexpr int SGN(Scalar val) { return (static_cast(0) < val) - (val < static_cast(0));