diff --git a/bin/filter_mrc/filter_mrc.cpp b/bin/filter_mrc/filter_mrc.cpp index a12af62..032db96 100644 --- a/bin/filter_mrc/filter_mrc.cpp +++ b/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 diff --git a/bin/filter_mrc/handlers.cpp b/bin/filter_mrc/handlers.cpp index 2915e03..5d1a28c 100644 --- a/bin/filter_mrc/handlers.cpp +++ b/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); } diff --git a/bin/filter_mrc/settings.cpp b/bin/filter_mrc/settings.cpp index 6782210..5d58066 100644 --- a/bin/filter_mrc/settings.cpp +++ b/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::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& 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& 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& 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& 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& 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& 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& 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& 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) { diff --git a/bin/filter_mrc/settings.hpp b/bin/filter_mrc/settings.hpp index b04c2c3..1919e90 100644 --- a/bin/filter_mrc/settings.hpp +++ b/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 ---- // diff --git a/doc/doc_filter_mrc.md b/doc/doc_filter_mrc.md index be424d9..6845daa 100644 --- a/doc/doc_filter_mrc.md +++ b/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 diff --git a/lib/visfd/morphology.hpp b/lib/visfd/morphology.hpp index cda42a3..29bda9a 100644 --- a/lib/visfd/morphology.hpp +++ b/lib/visfd/morphology.hpp @@ -56,7 +56,7 @@ template > *pva_minima_crds, //!< store minima locations (ix,iy,iz) here (if not nullptr) vector > *pva_maxima_crds, //!< store maxima locations (ix,iy,iz) here (if not nullptr) vector *pv_minima_scores, //!< store corresponding minima aaafI[iz][iy][ix] values here (if not nullptr) @@ -67,8 +67,8 @@ FindExtrema(int const image_size[3], //!< size of the image in x,y,z di Scalar maxima_threshold=-std::numeric_limits::infinity(), //!< ignore maxima with intensities lower than this int connectivity=3, //!< square root of search radius around each voxel (1=nearest_neighbors, 2=diagonal2D, 3=diagonal3D) bool allow_borders=true, //!< if true, plateaus that touch the image border (or mask boundary) are valid extrema - Label ***aaaiDest=nullptr, //!< optional: create an image showing where the extrema are? - ostream *pReportProgress=nullptr //!< optional: print progress to the user? + Label ***aaaiDest=nullptr, //!< OPTIONAL: create an image showing where the extrema are? + ostream *pReportProgress=nullptr //!< OPTIONAL: print progress to the user? ) { bool find_minima = (pva_minima_crds != nullptr); @@ -167,8 +167,8 @@ FindExtrema(int const image_size[3], //!< size of input image array Scalar threshold=std::numeric_limits::infinity(), // Ignore minima or maxima which are not sufficiently low or high int connectivity=3, //!< square root of search radius around each voxel (1=nearest_neighbors, 2=diagonal2D, 3=diagonal3D) bool allow_borders=true, //!< if true, plateaus that touch the image border (or mask boundary) are valid extrema - Label ***aaaiDest=nullptr, //!< optional: create an image showing where the extrema are? - ostream *pReportProgress=nullptr) //!< print progress to the user? + Label ***aaaiDest=nullptr, //!< OPTIONAL: create an image showing where the extrema are? + ostream *pReportProgress=nullptr) //!< OPTIONAL: print progress to the user? { // NOTE: // C++ will not allow us to supply nullptr to a function that expects a pointer @@ -236,9 +236,9 @@ Dilate(vector > structure_factor, // a list of (ix,iy, const int image_size[3], //!< size of the image in x,y,z directions 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 - - ostream *pReportProgress = nullptr) + Scalar const *const *const *aaafMask = nullptr, //!< OPTIONAL: ignore voxel ix,iy,iz if aaafMask!=nullptr and aaafMask[iz][iy][ix]==0 + ostream *pReportProgress = nullptr //!< OPTIONAL: print progress to the user? + ) { for (int iz=0; iz < image_size[2]; iz++) { if (pReportProgress) @@ -293,9 +293,9 @@ Erode(vector > structure_factor, // a list of (ix,iy,i const int image_size[3], //!< size of the image in x,y,z directions 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 - - ostream *pReportProgress = nullptr) + Scalar const *const *const *aaafMask = nullptr, //!< OPTIONAL: ignore voxel ix,iy,iz if aaafMask!=nullptr and aaafMask[iz][iy][ix]==0 + ostream *pReportProgress = nullptr //!< OPTIONAL: print progress to the user? + ) { for (int iz=0; iz < image_size[2]; iz++) { if (pReportProgress) @@ -339,29 +339,39 @@ Erode(vector > structure_factor, // a list of (ix,iy,i /// https://en.wikipedia.org/wiki/Dilation_(morphology)#Grayscale_dilation template void -DilateSphere(Scalar radius, //!< radius of the sphere (in voxels) - int const image_size[3], //!< size of the image in x,y,z directions +DilateSphere(Scalar radius, //!< radius of the sphere (in voxels) + int const image_size[3], //!< size of the image in x,y,z directions 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 - 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 + Scalar const *const *const *aaafMask = nullptr, //!< OPTIONAL: ignore voxel ix,iy,iz if aaafMask!=nullptr and aaafMask[iz][iy][ix]==0 + Scalar radius_max = 0.0, //!< OPTIONAL: smoothly vary between two radii? + Scalar bmax = 0.0, //!< OPTIONAL: structure factor b varies between 0 and this value + ostream *pReportProgress = nullptr //!< OPTIONAL: 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); + int Ri = ceil(std::max(radius, radius_max)); 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_penalty == std::numeric_limits::infinity()) - { + if (bmax == 0.0) { Scalar r = sqrt(SQR(ix) + SQR(iy) + SQR(iz)); if (r <= radius) add_this_voxel = true; } + else if (radius_max > radius) { + Scalar r = sqrt(SQR(ix) + SQR(iy) + SQR(iz)); + if (r <= radius) + add_this_voxel = true; + else if (r <= radius_max) { + add_this_voxel = true; + b = -(r - radius) / (radius_max - radius); + b *= bmax; + } + } else { // Calculate the distance of all of the corners to the origin. The // corners of the cube are located at the corners of voxel ix,iy,iz @@ -394,11 +404,10 @@ DilateSphere(Scalar radius, //!< radius of the sphere (in voxels) add_this_voxel = false; else { add_this_voxel = true; - b = (r_max - radius) / (r_max - r_min); - b = -b; - b *= smooth_boundary_penalty; + b = -(r_max - radius) / (r_max - r_min); + b *= bmax; } - } // if (smooth_boundary) + } // if ((bmax != 0) && (radius_max > radius)) if (add_this_voxel) structure_factor.push_back(tuple(ix,iy,iz,b)); } // for (int ix=0; ix < image_size[0]; ix++) @@ -406,7 +415,8 @@ DilateSphere(Scalar radius, //!< radius of the sphere (in voxels) } // for (int iz=0; iz < image_size[2]; iz++) if (pReportProgress) - *pReportProgress << "DilateSphere(r="< void -ErodeSphere(Scalar radius, //!< radius of the sphere (in voxels) - int const image_size[3], //!< size of the image in x,y,z directions +ErodeSphere(Scalar radius, //!< radius of the sphere (in voxels) + int const image_size[3], //!< size of the image in x,y,z directions 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 - 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 + Scalar const *const *const *aaafMask = nullptr, //!< OPTIONAL: ignore voxel ix,iy,iz if aaafMask!=nullptr and aaafMask[iz][iy][ix]==0 + Scalar radius_max = 0.0, //!< OPTIONAL: smoothly vary between two radii? + Scalar bmax = 0.0, //!< OPTIONAL: structure factor b varies between 0 and this value + ostream *pReportProgress = nullptr //!< OPTIONAL: 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); + int Ri = ceil(std::max(radius, radius_max)); 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_penalty == std::numeric_limits::infinity()) - { + if (bmax == 0.0) { Scalar r = sqrt(SQR(ix) + SQR(iy) + SQR(iz)); if (r <= radius) add_this_voxel = true; } + else if (radius_max > radius) { + Scalar r = sqrt(SQR(ix) + SQR(iy) + SQR(iz)); + if (r <= radius) + add_this_voxel = true; + else if (r <= radius_max) { + add_this_voxel = true; + b = -(r - radius) / (radius_max - radius); + b *= bmax; + } + } else { // Calculate the distance of all of the corners to the origin. The // corners of the cube are located at the corners of voxel ix,iy,iz @@ -479,11 +499,10 @@ ErodeSphere(Scalar radius, //!< radius of the sphere (in voxels) add_this_voxel = false; else { add_this_voxel = true; - b = (r_max - radius) / (r_max - r_min); - b = -b; - b *= smooth_boundary_penalty; + b = -(r_max - radius) / (r_max - r_min); + b *= bmax; } - } // if (smooth_boundary) + } // if ((bmax != 0) && (radius_max > radius)) if (add_this_voxel) structure_factor.push_back(tuple(ix,iy,iz,b)); } // for (int ix=0; ix < image_size[0]; ix++) @@ -491,7 +510,8 @@ ErodeSphere(Scalar radius, //!< radius of the sphere (in voxels) } // for (int iz=0; iz < image_size[2]; iz++) if (pReportProgress) - *pReportProgress << "ErodeSphere(r="<::infinity(), //!< structure factor b varies between 0 and this value - ostream *pReportProgress = nullptr //!< report progress to the user + Scalar const *const *const *aaafMask = nullptr, //!< OPTIONAL: ignore voxel ix,iy,iz if aaafMask!=nullptr and aaafMask[iz][iy][ix]==0 + Scalar radius_max = 0.0, //!< OPTIONAL: smoothly vary between two radii? + Scalar bmax = 0.0, //!< OPTIONAL: structure factor b varies between 0 and this value + ostream *pReportProgress = nullptr //!< OPTIONAL: report progress to the user ) { float ***aaafTmp; //temporary space to store image after erosion @@ -527,7 +548,8 @@ OpenSphere(Scalar radius, //!< radius of the sphere (in voxels) aaafSource, aaafTmp, //<--save the resulting image in aaafTmp aaafMask, - smooth_boundary_penalty, + radius_max, + bmax, pReportProgress); DilateSphere(radius, @@ -535,7 +557,8 @@ 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_penalty, + radius_max, + bmax, pReportProgress); Dealloc3D(aaafTmp); @@ -553,9 +576,10 @@ CloseSphere(Scalar radius, //!< radius of the sphere (in voxels) const int image_size[3], //!< size of the image in x,y,z directions 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 - 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 + Scalar const *const *const *aaafMask = nullptr, //!< OPTIONAL: ignore voxel ix,iy,iz if aaafMask!=nullptr and aaafMask[iz][iy][ix]==0 + Scalar radius_max = 0.0, //!< OPTIONAL: smoothly vary between two radii? + Scalar bmax = 0.0, //!< OPTIONAL: structure factor b varies between 0 and this value + ostream *pReportProgress = nullptr //!< OPTIONAL: report progress to the user ) { float ***aaafTmp; //temporary space to store image after dilation @@ -567,7 +591,8 @@ CloseSphere(Scalar radius, //!< radius of the sphere (in voxels) aaafSource, aaafTmp, //<--save the resulting image in aaafTmp aaafMask, - smooth_boundary_penalty, + radius_max, + bmax, pReportProgress); ErodeSphere(radius, @@ -575,7 +600,8 @@ 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_penalty, + radius_max, + bmax, pReportProgress); Dealloc3D(aaafTmp); @@ -593,9 +619,10 @@ WhiteTopHatSphere(Scalar radius, //!< radius of the sphere (in voxel const int image_size[3], //!< size of the image in x,y,z directions Scalar const *const *const *aaafSource, //!::infinity(), //!< structure factor b varies between 0 and this value - ostream *pReportProgress = nullptr //!< report progress to the user + Scalar const *const *const *aaafMask = nullptr, //!< OPTIONAL: ignore voxel ix,iy,iz if aaafMask!=nullptr and aaafMask[iz][iy][ix]==0 + Scalar radius_max = 0.0, //!< OPTIONAL: smoothly vary between two radii? + Scalar bmax = 0.0, //!< OPTIONAL: structure factor b varies between 0 and this value + ostream *pReportProgress = nullptr //!< OPTIONAL: report progress to the user ) { float ***aaafTmp; //temporary space to store image after erosion @@ -608,7 +635,8 @@ WhiteTopHatSphere(Scalar radius, //!< radius of the sphere (in voxel aaafSource, aaafTmp, //<--save the resulting image in aaafTmp aaafMask, - smooth_boundary_penalty, + radius_max, + bmax, pReportProgress); // Subtract it from the original image brightness @@ -632,9 +660,10 @@ BlackTopHatSphere(Scalar radius, //!< radius of the sphere (in voxel const int image_size[3], //!< size of the image in x,y,z directions 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 - 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 + Scalar const *const *const *aaafMask = nullptr, //!< OPTIONAL: ignore voxel ix,iy,iz if aaafMask!=nullptr and aaafMask[iz][iy][ix]==0 + Scalar radius_max = 0.0, //!< OPTIONAL: smoothly vary between two radii? + Scalar bmax = 0.0, //!< OPTIONAL: structure factor b varies between 0 and this value + ostream *pReportProgress = nullptr //!< OPTIONAL: report progress to the user ) { float ***aaafTmp; //temporary space to store image after erosion @@ -647,7 +676,8 @@ BlackTopHatSphere(Scalar radius, //!< radius of the sphere (in voxel aaafSource, aaafTmp, //<--save the resulting image in aaafTmp aaafMask, - smooth_boundary_penalty, + radius_max, + bmax, pReportProgress); // Subtract the original image brightness