From e3825f34360b64bd26874150d12a56bf7dd1bfe1 Mon Sep 17 00:00:00 2001 From: Andrew Jewett Date: Thu, 8 Jul 2021 02:01:19 -0700 Subject: [PATCH] reduced the influence of noise near mask boundaries, fixed a few small bugs, and updated the documentation --- README.md | 3 +- bin/filter_mrc/feature_variants.hpp | 5 +- bin/filter_mrc/filter3d_variants.hpp | 291 +++++++++++- bin/filter_mrc/filter_mrc.cpp | 6 +- bin/filter_mrc/handlers.cpp | 25 +- bin/filter_mrc/handlers_unsupported.cpp | 8 +- bin/filter_mrc/settings.cpp | 19 + bin/filter_mrc/settings.hpp | 1 + bin/pval_mrc/pval_mrc.cpp | 2 +- doc/doc_filter_mrc.md | 143 ++++-- lib/visfd/feature.hpp | 401 +++++++--------- lib/visfd/feature_implementation.hpp | 608 ++---------------------- lib/visfd/filter2d.hpp | 124 ----- lib/visfd/filter3d.hpp | 187 +------- lib/visfd/morphology.hpp | 226 +++++++++ lib/visfd/morphology_implementation.hpp | 596 +++++++++++++++++++++++ lib/visfd/segmentation.hpp | 2 +- lib/visfd/visfd.hpp | 7 +- lib/visfd/visfd_utils.hpp | 306 +++--------- 19 files changed, 1545 insertions(+), 1415 deletions(-) create mode 100644 lib/visfd/morphology.hpp create mode 100644 lib/visfd/morphology_implementation.hpp diff --git a/README.md b/README.md index 73b1cf6..3b64a1e 100644 --- a/README.md +++ b/README.md @@ -28,8 +28,7 @@ Multiprocessor support is implemented using ## Alternatives to VISFD Some of the features of VISFD are also available from popular python libraries, including: -[scikit-image](https://scikit-image.org), -[opencv](https://opencv.org), +[scikit-image](https://scikit-image.org) and [scipy.ndimage](https://docs.scipy.org/doc/scipy/reference/ndimage.html) *(MRC files can be read in python using the diff --git a/bin/filter_mrc/feature_variants.hpp b/bin/filter_mrc/feature_variants.hpp index 4c552b7..64b7571 100644 --- a/bin/filter_mrc/feature_variants.hpp +++ b/bin/filter_mrc/feature_variants.hpp @@ -24,6 +24,7 @@ #include using namespace std; #include +#include #include "filter3d_variants.hpp" @@ -130,7 +131,7 @@ BlobDogNM(int const image_size[3], //! + /// @brief This is a version of the function of the same name defined in /// "filter2d.h". This version allows the user to control the width of the /// interval considered for the filter in 2 ways: "filter_truncate_ratio", @@ -112,6 +113,125 @@ GenFilterGenGauss3D(Scalar width[3], //!<"σ_x", "σ_y", "σ_z", parameters +/// @brief +/// Create a 2D filter and fill it with a difference of (generalized) Gaussians: +/// This version requires that the caller has already created individual +/// filters for the two gaussians. +/// All this function does is subtract one filter from the other (and rescale). + +template + +static Filter2D +_GenFilterDogg2D(Scalar width_a[2], //!< "a" parameter in formula + Scalar width_b[2], //!< "b" parameter in formula + Filter2D& filterXY_A, //!< filters for the two + Filter2D& filterXY_B, //!< gaussians + Scalar *pA=nullptr, //!< optional:report A,B coeffs to user + Scalar *pB=nullptr, //!< optional:report A,B coeffs to user + ostream *pReportProgress = nullptr //!< optional:report filter details to the user? + ) +{ + + Scalar A, B; + //A, B = height of the central peak + A = filterXY_A.aafH[0][0]; + B = filterXY_B.aafH[0][0]; + + // The "difference of gaussians" filter is the difference between + // these two (generalized) gaussian filters. + int halfwidth[2]; + halfwidth[0] = std::max(filterXY_A.halfwidth[0], filterXY_B.halfwidth[0]); + halfwidth[1] = std::max(filterXY_A.halfwidth[1], filterXY_B.halfwidth[1]); + Filter2D filter(halfwidth); + + if (pReportProgress) + *pReportProgress << "Array of 2D filter entries:" << endl; + + for (int iy=-halfwidth[1]; iy<=halfwidth[1]; iy++) { + for (int ix=-halfwidth[0]; ix<=halfwidth[0]; ix++) { + filter.aafH[iy][ix] = 0.0; + + // The two filters may have different widths, so we have to check + // that ix and iy lie within the domain of these two filters before + // adding or subtracting their values from the final GDOG filter. + if (((-filterXY_A.halfwidth[0]<=ix) && (ix<=filterXY_A.halfwidth[0])) && + ((-filterXY_A.halfwidth[1]<=iy) && (iy<=filterXY_A.halfwidth[1]))) + + filter.aafH[iy][ix] +=filterXY_A.aafH[iy][ix]; // /(A-B); COMMENTING OUT + + // COMMENTING OUT: (The factor of 1/(A-B) insures that the central peak has height 1) + + if (((-filterXY_B.halfwidth[0]<=ix) && (ix<=filterXY_B.halfwidth[0])) && + ((-filterXY_B.halfwidth[1]<=iy) && (iy<=filterXY_B.halfwidth[1]))) + + filter.aafH[iy][ix] -=filterXY_B.aafH[iy][ix]; // /(A-B); COMMENTING OUT + + + + //FOR DEBUGGING REMOVE EVENTUALLY + //if (pReportProgress) + // *pReportProgress << "GenDogg2D: aafH["< +// Create a 2D filter and fill it with a difference of (generalized) Gaussians: +Filter2D +GenFilterDogg2D(Scalar width_a[2], //"a" parameter in formula + Scalar width_b[2], //"b" parameter in formula + Scalar m_exp, //"m" parameter in formula + Scalar n_exp, //"n" parameter in formula + int halfwidth[2], + Scalar *pA = nullptr, //optional:report A,B coeffs to user + Scalar *pB = nullptr, //optional:report A,B coeffs to user + ostream *pReportProgress = nullptr) +{ + Filter2D filterXY_A = + GenFilterGenGauss2D(width_a, //"a_x", "a_y" gaussian width parameters + m_exp, //"n" exponent parameter + halfwidth); + //pReportProgress); + + Filter2D filterXY_B = + GenFilterGenGauss2D(width_b, //"b_x", "b_y" gaussian width parameters + n_exp, //"n" exponent parameter + halfwidth); + //pReportProgress); + + return _GenFilterDogg2D(width_a, + width_b, + filterXY_A, filterXY_B, + pA, + pB, + pReportProgress); +} //GenFilterDogg2D(...halfwidth...) + + + /// @brief This is a version of the function of the same name defined in /// "filter2d.h". This version allows the user to control the width of the /// interval considered for the filter in 2 ways: "filter_truncate_ratio", @@ -155,8 +275,6 @@ GenFilterDogg2D(Scalar width_a[2], //!< "a" parameter in formula return _GenFilterDogg2D(width_a, width_b, - m_exp, - n_exp, filterXY_A, filterXY_B, pA, pB, @@ -166,6 +284,175 @@ GenFilterDogg2D(Scalar width_a[2], //!< "a" parameter in formula +/// @brief Create a 3D filter and fill it with a difference of (generalized) Gaussians +/// +/// This version requires that the caller has already created individual +/// filters for the two gaussians. +/// All this function does is subtract one filter from the other (and rescale). +/// @note: DEPRECIATION WARNING: It's not clear if this filter is useful. +/// I may delete this function in the future. +/// @note: THIS FUNCTION WAS NOT INTENDED FOR PUBLIC USE +template + +static Filter3D +_GenFilterDogg3D(Scalar width_a[3], //!< "a" parameter in formula + Scalar width_b[3], //!< "b" parameter in formula + Scalar m_exp, //!< "m" parameter in formula + Scalar n_exp, //!< "n" parameter in formula + Filter3D& filter_A, //!< filters for the two + Filter3D& filter_B, //!< gaussians + Scalar *pA=nullptr, //!< optional:report A,B coeffs to user + Scalar *pB=nullptr, //!< optional:report A,B coeffs to user + ostream *pReportEquation = nullptr //!< optional: report equation to the user + ) +{ + Scalar A, B; + //A, B = height of the central peak + A = filter_A.aaafH[0][0][0]; + B = filter_B.aaafH[0][0][0]; + + + // The "difference of gaussians" filter is the difference between + // these two (generalized) gaussian filters. + int halfwidth[3]; + halfwidth[0] = std::max(filter_A.halfwidth[0], filter_B.halfwidth[0]); + halfwidth[1] = std::max(filter_A.halfwidth[1], filter_B.halfwidth[1]); + halfwidth[2] = std::max(filter_A.halfwidth[2], filter_B.halfwidth[2]); + Filter3D filter(halfwidth); + + //FOR DEBUGGING REMOVE EVENTUALLY + //if (pReportEquation) + // *pReportEquation << "Array of 3D filter entries:" << endl; + + + for (int iz=-halfwidth[2]; iz<=halfwidth[2]; iz++) { + for (int iy=-halfwidth[1]; iy<=halfwidth[1]; iy++) { + for (int ix=-halfwidth[0]; ix<=halfwidth[0]; ix++) { + + filter.aaafH[iz][iy][ix] = 0.0; + + // The two filters may have different widths, so we have to check + // that ix,iy and iz lie within the domain of these two filters before + // adding or subtracting their values from the final GDoG filter. + if (((-filter_A.halfwidth[0]<=ix) && (ix<=filter_A.halfwidth[0])) && + ((-filter_A.halfwidth[1]<=iy) && (iy<=filter_A.halfwidth[1])) && + ((-filter_A.halfwidth[2]<=iz) && (iz<=filter_A.halfwidth[2]))) + + filter.aaafH[iz][iy][ix] += + filter_A.aaafH[iz][iy][ix]; // /(A-B); COMMENTING OUT + + + + // COMMENTING OUT: (The factor of 1/(A-B) insures that the central peak has height 1) + + + if (((-filter_B.halfwidth[0]<=ix) && (ix<=filter_B.halfwidth[0])) && + ((-filter_B.halfwidth[1]<=iy) && (iy<=filter_B.halfwidth[1])) && + ((-filter_B.halfwidth[2]<=iz) && (iz<=filter_B.halfwidth[2]))) + + filter.aaafH[iz][iy][ix] -= + filter_B.aaafH[iz][iy][ix]; // /(A-B); COMMENTING OUT + + //*pReportEquation << aaafH[iz][iy][ix]; + // + //if (ix == 0) pReportEquation << "\n"; else pReportEquation << " "; + + } // for (int ix=-halfwidth[0]; ix<=halfwidth[0]; ix++) { + } // for (int iy=-halfwidth[1]; iy<=halfwidth[1]; iy++) { + } // for (int iz=-halfwidth[2]; iz<=halfwidth[2]; iz++) { + + + // COMMENTING OUT the factor of 1/(A-B): + //A = A/(A-B); + //B = B/(A-B); + + if (pA && pB) { + *pA = A; // Rescale A and B numbers returned to the caller + *pB = B; // (because we divided the array entries by (A-B) earlier) + } + + if (pReportEquation) { + *pReportEquation << "\n"; + *pReportEquation << " Filter Used:\n" + " h(x,y,z) = h_a(x,y,z) - h_b(x,y,z)\n" + " h_a(x,y,z) = A*exp(-((x/a_x)^2 + (y/a_y)^2 + (z/a_z)^2)^(m/2))\n" + " h_b(x,y,z) = B*exp(-((x/b_x)^2 + (y/b_y)^2 + (z/b_z)^2)^(n/2))\n" + " ... where A = " << A << "\n" + " B = " << B << "\n" + " m = " << m_exp << "\n" + " n = " << n_exp << "\n" + " (a_x, a_y, a_z) = " + << "(" << width_a[0] + << " " << width_a[1] + << " " << width_a[2] << ")\n" + " (b_x, b_y, b_z) = " + << "(" << width_b[0] + << " " << width_b[1] + << " " << width_b[2] << ")\n"; + *pReportEquation << " You can plot a slice of this function\n" + << " in the X direction using:\n" + " draw_filter_1D.py -dogg " << A << " " << B + << " " << width_a[0] << " " << width_b[0] + << " " << m_exp << " " << n_exp << endl; + *pReportEquation << " and in the Y direction using:\n" + " draw_filter_1D.py -dogg " << A << " " << B + << " " << width_a[1] << " " << width_b[1] + << " " << m_exp << " " << n_exp << endl; + *pReportEquation << " and in the Z direction using:\n" + " draw_filter_1D.py -dogg " << A << " " << B + << " " << width_a[2] << " " << width_b[2] + << " " << m_exp << " " << n_exp << endl; + } //if (pReportEquation) + + return filter; +} //_GenFilterDogg3D() + + + + +/// @brief Create a 3D filter and fill it with a difference of (generalized) Gaussians +/// @verbatim h(x,y,z) = A*exp(-(r/a)^m) - B*exp(-(r/b)^n) @endverbatim +/// where @verbatim r = sqrt(x^2 + y^2 + z^2) @endverbatim +/// and "A" and "B" are determined by normalization of each term independently +/// DEPRECIATION WARNING: It's not clear if this type if filter is useful. +/// I may delete this function in the future. + +template + +Filter3D +GenFilterDogg3D(Scalar width_a[3], //!< "a" parameter in formula + Scalar width_b[3], //!< "b" parameter in formula + Scalar m_exp, //!< "m" parameter in formula + Scalar n_exp, //!< "n" parameter in formula + int halfwidth[3], //!< the width of the filter + Scalar *pA=nullptr, //!< optional:report A,B coeffs to user + Scalar *pB=nullptr, //!< optional:report A,B coeffs to user + ostream *pReportEquation = nullptr //!< optional: print params used? + ) +{ + Filter3D filter_A = + GenFilterGenGauss3D(width_a, //"a_x", "a_y" gaussian width parameters + m_exp, //"m" exponent parameter + halfwidth); + + Filter3D filter_B = + GenFilterGenGauss3D(width_b, //"b_x", "b_y" gaussian width parameters + n_exp, //"n" exponent parameter + halfwidth); + + // The next function is defind in "filter3d_implementation.hpp" + return _GenFilterDogg3D(width_a, + width_b, + m_exp, + n_exp, + filter_A, filter_B, + pA, + pB, + pReportEquation); +} //GenFilterDogg3D(...halfwidth...) + + + /// @brief This is a version of the function of the same name defined in /// "filter3d.hpp". This version allows the user to control the width of the /// interval considered for the filter in 2 ways: "filter_truncate_ratio", diff --git a/bin/filter_mrc/filter_mrc.cpp b/bin/filter_mrc/filter_mrc.cpp index e5b5b03..4ba326a 100644 --- a/bin/filter_mrc/filter_mrc.cpp +++ b/bin/filter_mrc/filter_mrc.cpp @@ -29,8 +29,8 @@ using namespace std; string g_program_name("filter_mrc"); -string g_version_string("0.27.3"); -string g_date_string("2021-7-06"); +string g_version_string("0.27.4"); +string g_date_string("2021-7-08"); @@ -232,7 +232,7 @@ int main(int argc, char **argv) { for (int iz=0; iz #include #include +#include #include #include "err.hpp" #include "settings.hpp" @@ -63,7 +64,7 @@ HandleGGauss(Settings settings, tomo_in.aaafI, tomo_out.aaafI, mask.aaafI, - true, + settings.normalize_near_boundaries, &cerr); } //HandleGGauss() @@ -89,7 +90,7 @@ HandleGauss(Settings settings, settings.width_a, settings.filter_truncate_ratio, settings.filter_truncate_threshold, - true, + settings.normalize_near_boundaries, &cerr); cerr << " Filter Used:\n" @@ -369,7 +370,7 @@ HandleBlobsNonmaxSuppression(Settings settings, settings.nonmax_min_radial_separation_ratio, settings.nonmax_max_volume_overlap_large, settings.nonmax_max_volume_overlap_small, - PRIORITIZE_HIGH_MAGNITUDE_SCORES, + SORT_DECREASING_MAGNITUDE, &cerr); cerr << " " << crds.size() << " blobs remaining" << endl; @@ -388,7 +389,7 @@ HandleBlobsNonmaxSuppression(Settings settings, scores, settings.training_pos_crds, settings.training_neg_crds, - PRIORITIZE_HIGH_MAGNITUDE_SCORES, + SORT_DECREASING_MAGNITUDE, &settings.score_lower_bound, &settings.score_upper_bound, &cerr); @@ -455,7 +456,7 @@ HandleBlobScoreSupervisedMulti(Settings settings, settings.multi_training_neg_crds, &threshold_lower_bound, &threshold_upper_bound, - PRIORITIZE_HIGH_MAGNITUDE_SCORES, + SORT_DECREASING_MAGNITUDE, &cerr); } //HandleBlobScoreSupervisedMulti() @@ -949,7 +950,7 @@ HandleExtrema(Settings settings, settings.nonmax_min_radial_separation_ratio, settings.nonmax_max_volume_overlap_large, settings.nonmax_max_volume_overlap_small, - PRIORITIZE_LOW_SCORES, + SORT_INCREASING, &cerr); } @@ -961,7 +962,7 @@ HandleExtrema(Settings settings, settings.nonmax_min_radial_separation_ratio, settings.nonmax_max_volume_overlap_large, settings.nonmax_max_volume_overlap_small, - PRIORITIZE_HIGH_SCORES, + SORT_DECREASING, &cerr); } @@ -1023,7 +1024,7 @@ HandleLocalFluctuations(Settings settings, settings.template_background_exponent, settings.filter_truncate_ratio, settings.filter_truncate_threshold, - true, + settings.normalize_near_boundaries, &cerr); tomo_out.FindMinMaxMean(); @@ -1269,7 +1270,7 @@ HandleTV(Settings settings, mask.aaafI, settings.width_b[0], truncate_halfwidth, - true, + settings.normalize_near_boundaries, &cerr); truncate_halfwidth = floor(settings.width_a[0] * @@ -1281,7 +1282,7 @@ HandleTV(Settings settings, mask.aaafI, settings.width_a[0], truncate_halfwidth, - true, + settings.normalize_near_boundaries, &cerr); } @@ -1523,8 +1524,8 @@ HandleTV(Settings settings, mask.aaafI, (settings.filter_type == Settings::CURVE), //settings.hessian_score_threshold, - true, // (do normalize near rectangular image bounaries) - false, // (diagonalize each tensor afterwards?) + false, // don't boost votes from voxels near mask boundary + false, // don't diagonalize each tensor afterwards &cerr); } // if (settings.load_intermediate_fname_base == "") diff --git a/bin/filter_mrc/handlers_unsupported.cpp b/bin/filter_mrc/handlers_unsupported.cpp index 69890bd..b63feed 100644 --- a/bin/filter_mrc/handlers_unsupported.cpp +++ b/bin/filter_mrc/handlers_unsupported.cpp @@ -94,7 +94,7 @@ HandleDoggXY(Settings settings, afIZorig, afIZnew, afMask, - true); + settings.normalize_near_boundaries); //It would be wasteful to allocate a temporary array to store this //Instead store the result of the convolution in the original array: @@ -526,7 +526,7 @@ HandleBootstrapDogg(Settings settings, //Initialize the random number generator we will need later: RANDOM_INIT(settings.bs_random_seed); - + for (int i_sim = 0; i_sim < settings.bs_ntests; i_sim++) { cerr << "\n" @@ -854,14 +854,14 @@ HandleTemplateGGauss(Settings settings, template_background_sigma, settings.filter_truncate_ratio, settings.filter_truncate_threshold, - true); + settings.normalize_near_boundaries); } else { w.Apply(tomo_in.header.nvoxels, tomo_in.aaafI, P.aaafI, // <-- save result here mask.aaafI, - true, + settings.normalize_near_boundaries, &cerr); } diff --git a/bin/filter_mrc/settings.cpp b/bin/filter_mrc/settings.cpp index 9727bc7..e813a46 100644 --- a/bin/filter_mrc/settings.cpp +++ b/bin/filter_mrc/settings.cpp @@ -31,6 +31,7 @@ Settings::Settings() { out_file_name = ""; out_file_overwrite = false; mask_file_name = ""; + normalize_near_boundaries = true; mask_select = 1; use_mask_select = false; mask_out = 0.0; @@ -310,6 +311,24 @@ Settings::ParseArgs(vector& vArgs) } // if (vArgs[i] == "-mask") + + else if (vArgs[i] == "-normalize-filters") + { + InputErr arg_err("Error: The " + vArgs[i] + + " argument must be followed by \"yes\" or \"no\".\n"); + if ((i+1 >= vArgs.size()) || (vArgs[i+1] == "") || (vArgs[i+1][0] == '-')) + throw arg_err; + if (vArgs[i+1] == "yes") + normalize_near_boundaries = true; + if (vArgs[i+1] == "no") + normalize_near_boundaries = false; + else + throw arg_err; + num_arguments_deleted = 2; + } // if (vArgs[i] == "-mask") + + + else if (vArgs[i] == "-mask-select") { try { diff --git a/bin/filter_mrc/settings.hpp b/bin/filter_mrc/settings.hpp index 0af7acc..8b35729 100644 --- a/bin/filter_mrc/settings.hpp +++ b/bin/filter_mrc/settings.hpp @@ -81,6 +81,7 @@ class Settings { string out_file_name; // name of the image file we want to create bool out_file_overwrite; // allow out_file_name to equal in_file_name? // Mask parameters are used to select (ignore) voxels from the original image. + bool normalize_near_boundaries; // compensate when mask boundary cuts filter? string mask_file_name; // name of an image file used for masking bool use_mask_select; // do we select voxels with a specific value? int mask_select; // select only voxels from the input with this value diff --git a/bin/pval_mrc/pval_mrc.cpp b/bin/pval_mrc/pval_mrc.cpp index 0b50122..e3df4e7 100644 --- a/bin/pval_mrc/pval_mrc.cpp +++ b/bin/pval_mrc/pval_mrc.cpp @@ -17,7 +17,7 @@ using namespace std; //#include #include -#include +#include using namespace visfd; #include "err.hpp" diff --git a/doc/doc_filter_mrc.md b/doc/doc_filter_mrc.md index d02d074..b3382cc 100644 --- a/doc/doc_filter_mrc.md +++ b/doc/doc_filter_mrc.md @@ -45,8 +45,7 @@ whose boundaries are smooth and gradual as opposed to jagged and rectangular, ### Alternatives to filter_mrc Many of the features of **filter_mrc** are also available from popular python libraries, including: -[scikit-image](https://scikit-image.org), -[opencv](https://opencv.org), +[scikit-image](https://scikit-image.org) and [scipy.ndimage](https://docs.scipy.org/doc/scipy/reference/ndimage.html) MRC files can be read in python using the @@ -214,7 +213,7 @@ of the surface by approximately 40 Angstroms inward, creating a new file ``` filter_mrc -i segmented.rec -o segmented_interior.rec -erode-gauss 40 ``` -*(See [-erode-gauss](#-erode-gauss-distance) for details.)* +*(See [-erode-gauss](#-erode-gauss-thickness) for details.)* **Note:** @@ -467,7 +466,54 @@ and improve the signal-to-noise ratio. (enough to store at least 9 copies of the original image in 32bit-float mode).* -### -edge thickness +#### Detection and masking + +You may wish to ignore certain voxels which are not part of the region +that you want to consider. +This can reduce the computation time. +In some cases, it can also improve detection accuracy. +You can do this using the ["**-mask**"](#-mask-MRC_FILE) argument +and [other similar arguments](#Masking). + +You can use an image mask during the initial stages of detection +(for example, using the [-surface](#-surface-thickness), +[-curve](#-surface-thickness), +[-edge](#-surface-thickness), and +[-tv](#-tv-σ_ratio) arguments). +However since detection is very sensitive to noise near the mask boundary, +it may be better to refrain from using a mask until later on, +(for example, when using the +[-connect](#-connect-threshold) and +["**-surface-normals-file**"](#-surface-normals-file-PLY_FILE) +arguments). Otherwise, you may detect many faint, minor objects or spurious +noise near this boundary instead of the surface or curve that you are seeking. +Alternatively, to get around this problem you can dilate (expand) the size +of the mask region to include these voxels near the periphery. +*(One way to do this is using the +[-dilate-gauss](#-dilate-gauss-thickness) +argument with a large thickness parameter. This is discussed below. +Alternatively, if your only goal is to reduce the computation time, you can +use the [-mask-rect(#-mask-rect-xmin-xmax-ymin-ymax-zmin-zmax) and +[-mask-sphere](#-mask-sphere-x0-y0-z0-r) arguments to create a mask +region large enough to enclose the object of interest.)* + +It may also help to vary the *size of the region* in the mask file using +[erosion](https://en.wikipedia.org/wiki/Erosion_(morphology)) or +[dilation](https://en.wikipedia.org/wiki/Dilation_(morphology)). +(If you are using the "-mask" argument to load a mask image file, +you can do this by creating a new mask file using *filter_mrc* again with the +[-erode-gauss](#-erode-gauss-thickness) or +[-dilate-gauss](#-dilate-gauss-thickness) arguments.) +Eroding the mask can elimate traces of any structure at the edge +of the mask boundary that might otherwise confuse or distract the detector. +But, as mentioned earlier, dilating the mask is frequently more effective +at improving the detection accuracy near the mask boundary. + +For noisy or cluttered images, you may have to experiment with all of these +options to maximize the chance of success. + + +### -edge distance ***WARNING: THIS FEATURE HAS NOT BEEN TESTED AND PROBABLY DOES NOT WORK 2021-6-21*** @@ -476,7 +522,7 @@ If the "**-edge**" filter is selected, the program will attempt to detect the sharp boundary surfaces between light regions and a dark volumetric regions in the image. (This filter detects gradients in image brightness.) -Unfortunately real images are noisy. Image noise will create many spurious +Unfortunately real images are noisy. Typical images will have many spurious bumps which must be filtered or smoothed away before edge detection begins. Consequently the user must specify a *thickness* argument which indicates amount of smoothing to use to erase features that you wish to ignore @@ -484,7 +530,8 @@ amount of smoothing to use to erase features that you wish to ignore This means that light or dark patches in the image which are smaller than *thickness* will not be visible to the edge-detector. Using a larger *thickness* will reduce the noise in the image, but may also -smooth over small (high-frequency) features in image you wish to detect. +smooth over small (high-frequency) features in boundary surface +that you want to detect. As with **-surface** argument, the edge detection fidelity can be improved using tensor voting (with the [-tv](#-tv-σ_ratio) argument). @@ -497,6 +544,11 @@ such as Angstroms, not in voxels.) *This filter requires a very large amount of memory (enough to store at least 9 copies of the original image in 32bit-float mode).* +#### edge detection and masking +It is recommended that you refrain from using any of the "-mask" arguments +***during*** edge detection. +See [here](#Detection-and-masking) for alternative strategies. + ## Detecting curves @@ -534,6 +586,11 @@ argument. *This filter requires a very large amount of memory (enough to store at least 9 copies of the original image in 32bit-float mode).* +#### curve detection and masking +It is recommended that you refrain from using any of the "-mask" arguments +***during*** curve detection. +See [here](#Detection-and-masking) for alternative strategies. + ### -tv σ_ratio @@ -961,16 +1018,23 @@ process of running tensor-voting each time. We can do this using the "-load-progress FILENAME" argument. The *FILENAME* argument should match the argument you used with "-save-progress" earlier. -Note that when using "-load-progress", you must respecify all -of the same arguments that you used earlier when you used -"-save-progress". -(Note: The "-out FNAME" argument is an exception. -When you use "-load-progress", you are permitted to change the "-out" argument.) +Be warned that files generated by "-save-progress" will consume +a substantial amount of disk space. You should delete them eventually. + + +Note that when using "-load-progress", you must again specify +all the same arguments that you used earlier when you used +"-save-progress" +(including the "-in" "-surface", "-edge", and "-curve", and "-tv" arguments). +When you do this, you are permitted to make some changes which are listed below. +When you use "-load-progress", you are permitted to change the "-out" argument. +You can also add an image mask (or change the mask, for example +using the "-mask", "-mask-rect", "-mask-sphere" arguments). +Adding a mask is useful if you want to restrict the detection +of surface features to certain regions inside the image, for example. +(It makes sense to start using a mask now because +[masks are not typically used during the initial stages of detection](#Surface-detection-and-masking).) -Be warned that using "-save-progress" will also consume substantial disk space -because each of these 6 files files may be up to 4 times larger than the -size of the original tomogram file. -(You should probably delete these files eventually.) @@ -1165,7 +1229,6 @@ Note: If **type** is set to *all*, then two different files will be generated whose names will end in ".minima.txt" and ".maxima.txt", respectively. - #### Blob detection and masking It is recommended that you refrain from specifying a mask @@ -1173,21 +1236,19 @@ It is recommended that you refrain from specifying a mask (for example, by using the ["**-mask**"](#-mask-MRC_FILE) argument). -Doing so could cause blobs which are near the periphery of the mask -to be ignored. +This is because the detector is more sensitive to noise in the image +near the mask boundary (or image boundary). Instead, it is recommended that you use the -["**-mask**"](#-mask-MRC_FILE) -argument -*later on, after you have finished blob detection*. +"-mask" argument *later on, after you have finished blob detection*. (This was shown in [example 2](#Example-2) above.) -Sometimes you might want to use the -["**-mask**"](#-mask-MRC_FILE) -argument during blob detection because it can -accelerate the search for blobs -(by ignoring voxels outside the mask). -However don't do this if you care -about blobs near the periphery. +Sometimes you might want to use the "-mask" argument during blob detection +because it can accelerate the search for blobs (by ignoring voxels outside +the mask). However don't do this if you care about blobs near the periphery. +Alternatively, to get around this problem you can dilate (expand) the size +of the mask region to include these voxels near the periphery. +*(One way to do this is using the +[-dilate-gauss](#-dilate-gauss-thickness) argument.)* #### Non-max suppression: automatic disposal of blobs @@ -2453,7 +2514,7 @@ This command can be used to *remove* the voxels in a existing mask which belong to a rectangular region bounded by the 6 numeric arguments: *xmin*, *xmax*, *ymin*, *ymax*, *zmin*, and *zmax*. This command can be issued multiple times as part of an ordered chain of -["**-mask**"](#-mask-MRC_FILE), +[-mask](#-mask-MRC_FILE), [-mask-rect](#-mask-rect-xmin-xmax-ymin-ymax-zmin-zmax), [-mask-sphere](#-mask-sphere-x0-y0-z0-r), and [-mask-sphere-subtract](#-mask-sphere-subtract-x0-y0-z0-r) @@ -2497,7 +2558,7 @@ not a number). This command can be used to *remove* the voxels in a existing mask which belong to a spherical region centered at *x0*, *y0*, *z0*, with radius *r*. This command can be issued multiple times as part of an ordered chain of -["**-mask**"](#-mask-MRC_FILE), +[-mask](#-mask-MRC_FILE), [-mask-sphere](#-mask-sphere-x0-y0-z0-r), and [-mask-rect](#-mask-rect-xmin-xmax-ymin-ymax-zmin-zmax), [-mask-rect-subtract](#-mask-rect-subtract-xmin-xmax-ymin-ymax-zmin-zmax), @@ -2517,8 +2578,8 @@ This provides a way to see where the mask region is located. -### -erode-gauss distance -### -dilate-gauss distance +### -erode-gauss thickness +### -dilate-gauss thickness *(Note: This is a crude implementation of @@ -2550,14 +2611,14 @@ ignored or considered during subsequent calculations)*. I will refer to an image (of 1s and 0s) as a "binary image". The **-dilate-gauss** and **-erode-gauss** arguments will enlarge or shrink the size of the region stored in a binary image by -a distance of roughly *distance* (which has units of physical distance +a distance of roughly *thickness* (which has units of physical distance unless the "*-w 1*" argument is used). Note that, due to blurring, features in the image which are smaller or -narrower than *distance* will be lost after this operation is completed. +narrower than *thickness* will be lost after this operation is completed. Consquently, for large distances, it might work better to -repeat this operation multiple time with the same small *distance* argument, -instead of using a single large *distance* argument. +repeat this operation multiple time with the same small *thickness* argument, +instead of using a single large *thickness* argument. You will have to experiment to see what works best. #### Details @@ -2566,10 +2627,10 @@ so that all voxels whose brightness are below a threshold are discarded. Consequently, you cannot combine these arguments with either the *-gauss* or *-thresh* arguments, because: -- "-dilate-gauss *distance*" is equivalent to -"-gauss *distance* -thresh 0.157299", *(where 0.157299 ≈ 1-erf(1))* -- "-erode-gauss *distance*" is equivalent to -"-gauss *distance* -thresh 0.842701", *(where 0.842701 ≈ erf(1))* +- "-dilate-gauss *thickness*" is equivalent to +"-gauss *thickness* -thresh 0.157299", *(where 0.157299 ≈ 1-erf(1))* +- "-erode-gauss *thickness*" is equivalent to +"-gauss *thickness* -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), @@ -2590,7 +2651,7 @@ brightness less than, say, 0.25, then this is a more forgiving criteria. This will select *more* voxels, effectively expanding the size of the selected region from the binary image. The thresholds above were chosen so that the selected region -expands or shrinks by roughly a distance of *distance*. +expands or shrinks by roughly a distance of *thickness*. ## Miscellaneous diff --git a/lib/visfd/feature.hpp b/lib/visfd/feature.hpp index c92f14e..13e2d4f 100644 --- a/lib/visfd/feature.hpp +++ b/lib/visfd/feature.hpp @@ -18,7 +18,7 @@ using namespace std; #include // defines the "VisfdErr" exception type #include // defines matrix diagonalizer (DiagonalizeSym3()) -#include // defines invert_permutation(), SortBlobs(), FindSpheres() ... +#include // defines invert_permutation(), FindSpheres() #include // defines Alloc2D() and Dealloc2D() #include // defines Alloc3D() and Dealloc3D() #include // defines "Filter1D" (used in ApplySeparable()) @@ -31,197 +31,6 @@ namespace visfd { -/// @brief Find all of the local minima and local maxima in an image (aaafI). -/// The locations of minima and maxima are stored in the -/// *pv_minima_crds and *pv_maxima_crds arrays, and sorted in -/// increasing and decreasing order respectively. (IE they are sorted -/// so that the most significant local minima or maxima will appear -/// first in the list.) -/// If either pv_minima_crds or pv_maxima_crds is nullptr, then -/// the minima or maxima will be ignored. -/// The optional pv_minima_scores and pv_maxima_scores store the -/// The caller can automatically discard minima or maxima which -/// are not sufficiently low or high, by supplying thresholds. -/// The optional aaafMask array (if not nullptr) can be used to ignore -/// certain voxels in the image (whose aaafMask entries are zero). -/// @note Local minima or maxima on the boundaries of the image -/// (or near the edge of the mask) -/// are not as trustworthy since some of the neighboring voxels -/// will not be available for comparison. These minima and maxima -/// can be ignored by setting allow_borders=false. The number of -/// neighbors around every voxel which are considered (eg, 6, 18, 26) -/// can be controlled using the "connectivity" argument. - -template - -void -FindExtrema(int const image_size[3], //!< size of the image in x,y,z directions - Scalar const *const *const *aaafI, //!< image array aaafI[iz][iy][ix] - Scalar const *const *const *aaafMask, //!< optional: ignore voxel ix,iy,iz if aaafMask!=nullptr and aaafMask[iz][iy][ix]==0 - vector > *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) - vector *pv_maxima_scores, //!< store corresponding maxima aaafI[iz][iy][ix] values here (if not nullptr) - vector *pv_minima_nvoxels, //!< store number of voxels in each minima (usually 1) - vector *pv_maxima_nvoxels, //!< store number of voxels in each maxima (usually 1) - Scalar minima_threshold=std::numeric_limits::infinity(), //!< ignore minima with intensities greater than this - Scalar maxima_threshold=-std::numeric_limits::infinity(), //!< ignore maxima with intensities lessr 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? - ) -{ - bool find_minima = (pva_minima_crds != nullptr); - bool find_maxima = (pva_maxima_crds != nullptr); - - vector minima_indices; - vector maxima_indices; - vector minima_scores; - vector maxima_scores; - vector *pv_minima_indices = nullptr; - vector *pv_maxima_indices = nullptr; - if (find_minima) { - pv_minima_indices = &minima_indices; - if (! pv_minima_scores) - pv_minima_scores = &minima_scores; - } - if (find_maxima) { - pv_maxima_indices = &maxima_indices; - if (! pv_maxima_scores) - pv_maxima_scores = &maxima_scores; - } - - _FindExtrema(image_size, - aaafI, - aaafMask, - pv_minima_indices, - pv_maxima_indices, - pv_minima_scores, - pv_maxima_scores, - pv_minima_nvoxels, - pv_maxima_nvoxels, - minima_threshold, - maxima_threshold, - connectivity, - allow_borders, - aaaiDest, - pReportProgress); - - // Now convert the indices back to x,y,z coordinates - if (pva_minima_crds) { - size_t N = minima_indices.size(); - assert(pva_minima_crds); - pva_minima_crds->resize(N); - for (size_t n = 0; n < N; n++) { - size_t i = minima_indices[n]; - // convert from a 1D index (i) to 3-D indices (ix, iy, iz) - size_t ix = i % image_size[0]; - i /= image_size[0]; - size_t iy = i % image_size[1]; - i /= image_size[1]; - size_t iz = i; - (*pva_minima_crds)[n][0] = ix; - (*pva_minima_crds)[n][1] = iy; - (*pva_minima_crds)[n][2] = iz; - } - } - if (pva_maxima_crds) { - size_t N = maxima_indices.size(); - assert(pva_maxima_crds); - pva_maxima_crds->resize(N); - for (size_t n = 0; n < N; n++) { - size_t i = maxima_indices[n]; - // convert from a 1D index (i) to 3-D indices (ix, iy, iz) - size_t ix = i % image_size[0]; - i /= image_size[0]; - size_t iy = i % image_size[1]; - i /= image_size[1]; - size_t iz = i; - (*pva_maxima_crds)[n][0] = ix; - (*pva_maxima_crds)[n][1] = iy; - (*pva_maxima_crds)[n][2] = iz; - } - } -} //FindExtrema() - - - - - -/// @brief The following version of this function seeks either minima -/// or maxima, but not both. (If you want both, use the other version. -/// That version is faster than invoking this function twice.) -/// See the description of that version for details. - -template - -void -FindExtrema(int const image_size[3], //!< size of input image array - Scalar const *const *const *aaafI, //!< input image array - Scalar const *const *const *aaafMask, //!< if not nullptr then zero entries indicate which voxels to ignore - vector > &extrema_crds, //!< store the location of each extrema - vector &extrema_scores, //!< store the brightness of each extrema - vector &extrema_nvoxels, //!< store number of voxels in each extrema (usually 1) - bool seek_minima=true, //!< search for minima or maxima? - 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? -{ - // NOTE: - // C++ will not allow us to supply nullptr to a function that expects a pointer - // to a template expression: Template argument deduction/substitution fails. - // We need to re-cast "nullptr" as a pointer with the correct type. - // One way to do that is to define these new versions of nullptr: - vector > *null_vai3 = nullptr; - vector *null_vf = nullptr; - vector *null_vi = nullptr; - - if (seek_minima) { - FindExtrema(image_size, - aaafI, - aaafMask, - &extrema_crds, // store minima locations here - null_vai3, // <-- don't search for maxima - &extrema_scores, // store minima values here - null_vf, // <-- don't search for maxima - &extrema_nvoxels, // store number of voxels in each minima - null_vi, // <-- don't search for maxima - threshold, - -std::numeric_limits::infinity(), - connectivity, - allow_borders, - aaaiDest, - pReportProgress); - } - else { - if (threshold == std::numeric_limits::infinity()) - threshold = -std::numeric_limits::infinity(); - FindExtrema(image_size, - aaafI, - aaafMask, - null_vai3, // <-- don't search for minima - &extrema_crds, // store maxima locations here - null_vf, // <-- don't search for minima - &extrema_scores, // store maxima values here - null_vi, // <-- don't search for minima - &extrema_nvoxels, // store number of voxels in each maxima - std::numeric_limits::infinity(), - threshold, - connectivity, - allow_borders, - aaaiDest, - pReportProgress); - } -} // FindExtrema() - - - - - - /// @brief Find all scale-invariant blobs in the image as a function of Gaussian /// width (ie. the "sigma" parameter), regardless of overlap. (A blob's /// "sigma" parameter represents the optimal width of the Gaussian blur @@ -718,39 +527,177 @@ BlobDogD(int const image_size[3], //! +void +SortBlobs(vector >& blob_crds,//!< x,y,z of each blob's center + vector& blob_diameters, //!< the width of each blob + vector& blob_scores, //!< the score for each blob + bool ascending_order = true, + bool ignore_score_sign = true, + vector *pPermutation = nullptr, //!< optional: return the new sorted order to the caller + ostream *pReportProgress = nullptr //!< optional: report progress to the user? + ) +{ + size_t n_blobs = blob_crds.size(); + assert(n_blobs == blob_diameters.size()); + assert(n_blobs == blob_scores.size()); + vector > score_index(n_blobs); + for (size_t i = 0; i < n_blobs; i++) { + if (ignore_score_sign) + score_index[i] = make_tuple(std::fabs(blob_scores[i]), i); + else + score_index[i] = make_tuple(blob_scores[i], i); + } + + if (n_blobs == 0) + return; + if (pReportProgress) + *pReportProgress << "-- Sorting blobs according to their scores... "; + if (ascending_order) + sort(score_index.begin(), + score_index.end()); + else + sort(score_index.rbegin(), + score_index.rend()); + + vector permutation(n_blobs); + for (size_t i = 0; i < score_index.size(); i++) + permutation[i] = get<1>(score_index[i]); + + if (pPermutation != nullptr) + *pPermutation = permutation; + + score_index.clear(); + apply_permutation(permutation, blob_crds); + apply_permutation(permutation, blob_diameters); + apply_permutation(permutation, blob_scores); + if (pReportProgress) + *pReportProgress << "done --" << endl; -/// @brief Calculate the volume of overlap between two spheres of radius -/// Ri and Rj separated by a distance of rij. +} //SortBlobs() -template -Scalar CalcSphereOverlap(Scalar rij,//! + +void +SortBlobs(vector >& blob_crds,//!< x,y,z of each blob's center + vector& blob_diameters, //!< the width of each blob + vector& blob_scores, //!< the score for each blob + SortCriteria sort_blob_criteria, //!< give priority to high or low scoring blobs? + bool ascending_order = true, + vector *pPermutation = nullptr, //!< optional: return the new sorted order to the caller + ostream *pReportProgress = nullptr //!< optional: report progress to the user? + ) { - // WLOG, assume Ri <= Rj. Insure that below - if (Ri > Rj) { - Scalar tmp = Ri; - Ri = Rj; - Rj = tmp; - } - if (rij <= Ri) { - return (4*M_PI/3) * Ri*Ri*Ri; + if (sort_blob_criteria == SORT_DECREASING) + SortBlobs(blob_crds, + blob_diameters, + blob_scores, + ascending_order, + false, + pPermutation, + pReportProgress); + else if (sort_blob_criteria == SORT_INCREASING) + SortBlobs(blob_crds, + blob_diameters, + blob_scores, + ! ascending_order, + false, + pPermutation, + pReportProgress); + else if (sort_blob_criteria == SORT_DECREASING_MAGNITUDE) + SortBlobs(blob_crds, + blob_diameters, + blob_scores, + ascending_order, + true, + pPermutation, + pReportProgress); + else if (sort_blob_criteria == SORT_INCREASING_MAGNITUDE) + SortBlobs(blob_crds, + blob_diameters, + blob_scores, + ! ascending_order, + true, + pPermutation, + pReportProgress); +} //SortBlobs() + + + + +/// @brief This variant of the function also takes care of discarding +/// training data which is not sufficiently close to any of the blobs. +/// It also throws an exception if the remaining training data is empty. +/// @return This function has no return value. +/// The results are stored in: +/// out_training_crds, +/// out_training_accepted, +/// out_training_scores + +template +static void + +FindBlobScores(const vector >& in_training_crds, //!< locations of blob-like things in training set + const vector& in_training_accepted, //!< classify each blob-like thing as "accepted" (true) or "rejected" (false) + vector >& out_training_crds, //!< same as in_training_crds after discarding entries too far away from any existing blobs + vector& out_training_accepted, //!< same as in_training_accepted after discarding entries too far away from any existing blobs + vector& out_training_scores, //!< store the scores of the remaining blob-like things here + const vector >& blob_crds, //!< location of each blob (in voxels, sorted by score in increasing priority) + const vector& blob_diameters, //!< diameger of each blob (sorted by score in increasing priority) + const vector& blob_scores, //!< priority of each blob (sorted by score in increasing priority) + SortCriteria sort_blob_criteria = SORT_DECREASING_MAGNITUDE, //!< give priority to high or low scoring blobs? + ostream *pReportProgress = nullptr //!< report progress back to the user? + ) +{ + // Figure out the score for each blob: + // in_training_crds[i] == position of ith training set data + // in_training_which_blob[i] == which blob contains this position? + // in_training_score[i] = score of this blob + + vector in_training_which_blob; // which blob (sphere) contains this position? + vector in_training_scores; // what is the score of that blob? + + // The next function is defined in "feather_implementation.hpp" + _FindBlobScores(in_training_crds, + in_training_scores, + in_training_which_blob, + blob_crds, + blob_diameters, + blob_scores, + sort_blob_criteria, + pReportProgress); + + size_t _N = in_training_crds.size(); + assert(_N == in_training_scores.size()); + assert(_N == in_training_accepted.size()); + + + // Discard training data that was not suffiently close to one of the blobs. + // Store the remaining training set data in the following variables: + const size_t UNOCCUPIED = 0; + for (size_t i=0; i < _N; i++) { + if (in_training_which_blob[i] != UNOCCUPIED) { + out_training_crds.push_back(in_training_crds[i]); + out_training_scores.push_back(in_training_scores[i]); + out_training_accepted.push_back(in_training_accepted[i]); + } } - // "xi" and "xj" are the distances from the sphere centers - // to the plane where the two spheres intersect. - Scalar xi = 0.5 * (1.0/rij) * (rij*rij + Ri*Ri - Rj*Rj); - Scalar xj = 0.5 * (1.0/rij) * (rij*rij + Rj*Rj - Ri*Ri); - Scalar volume_overlap = - (M_PI/3)*( Ri*Ri*Ri * (2 - (xi/Ri) * (3 - SQR(xi/Ri))) + - Rj*Rj*Rj * (2 - (xj/Rj) * (3 - SQR(xj/Rj))) ); - return volume_overlap; -} + // N = the number of training data left after discarding datum which + // do not lie within any blobs. (not within any blob radii) + size_t N = out_training_crds.size(); + assert(N == out_training_scores.size()); + assert(N == out_training_accepted.size()); + +} //FindBlobScores() @@ -767,11 +714,11 @@ Scalar CalcSphereOverlap(Scalar rij,//! // defines _FindExtrema() @@ -784,7 +731,7 @@ DiscardOverlappingBlobs(vector >& blob_crds, //!< location of ea Scalar min_radial_separation_ratio, //!< discard blobs if closer than this (ratio of sum of radii) Scalar max_volume_overlap_large=std::numeric_limits::infinity(), //!< discard blobs which overlap too much with the large blob (disabled by default; 1.0 would also do this) Scalar max_volume_overlap_small=std::numeric_limits::infinity(), //!< discard blobs which overlap too much with the small blob (disabled by default; 1.0 would also do this) - SortBlobCriteria sort_blob_criteria=PRIORITIZE_HIGH_MAGNITUDE_SCORES, //!< give priority to high or low scoring blobs? (See explanation above) + SortCriteria sort_blob_criteria=SORT_DECREASING_MAGNITUDE, //!< give priority to high or low scoring blobs? (See explanation above) ostream *pReportProgress = nullptr, //!< report progress back to the user? int scale=6 //! >& blob_crds, //!< locati const vector >& training_set_neg, //!< locations of blob-like things we want to ignore Scalar *pthreshold_lower_bound = nullptr, //!< return threshold to the caller Scalar *pthreshold_upper_bound = nullptr, //!< return threshold to the caller - SortBlobCriteria sort_blob_criteria = PRIORITIZE_HIGH_MAGNITUDE_SCORES, //!< give priority to high or low scoring blobs? + SortCriteria sort_blob_criteria = SORT_DECREASING_MAGNITUDE, //!< give priority to high or low scoring blobs? ostream *pReportProgress = nullptr //!< report progress back to the user? ) { @@ -1140,7 +1087,7 @@ ChooseBlobScoreThresholdsMulti(const vector > >& blob_crd const vector > >& training_set_neg, //!< locations of blob-like things Scalar *pthreshold_lower_bound = nullptr, //!< return threshold to the caller Scalar *pthreshold_upper_bound = nullptr, //!< return threshold to the caller - SortBlobCriteria sort_blob_criteria = PRIORITIZE_HIGH_MAGNITUDE_SCORES, //!< give priority to high or low scoring blobs? + SortCriteria sort_blob_criteria = SORT_DECREASING_MAGNITUDE, //!< give priority to high or low scoring blobs? ostream *pReportProgress = nullptr //!< report progress back to the user? ) { @@ -1202,7 +1149,7 @@ DiscardBlobsByScoreSupervised(vector >& blob_crds, //!< location vector& blob_scores, //!< priority of each blob const vector >& training_set_pos, //!< locations of blob-like things we are looking for const vector >& training_set_neg, //!< locations of blob-like things we want to ignore - SortBlobCriteria sort_blob_criteria = PRIORITIZE_HIGH_MAGNITUDE_SCORES, //!< give priority to high or low scoring blobs? + SortCriteria sort_blob_criteria = SORT_DECREASING_MAGNITUDE, //!< give priority to high or low scoring blobs? Scalar *pthreshold_lower_bound = nullptr, //!< optional: return threshold to the caller Scalar *pthreshold_upper_bound = nullptr, //!< optional: return threshold to the caller ostream *pReportProgress = nullptr //!< report progress back to the user? @@ -2667,7 +2614,7 @@ class TV3D { // instead of 2-dimensional surfaces (membranes, ...) // then the stick direction is assumed to be along the // of the curve, (as opposed to perpendicular to the 2D surface). - // As such + // Consequently: angle_dependence2 = sin2; } @@ -2841,7 +2788,7 @@ class TV3D { // instead of 2-dimensional surfaces (membranes, ...) // then the stick direction is assumed to be along the // of the curve, (as opposed to perpendicular to the 2D surface). - // As such + // Consequently: angle_dependence2 = sin2; } diff --git a/lib/visfd/feature_implementation.hpp b/lib/visfd/feature_implementation.hpp index 85b6d0d..77b91e9 100644 --- a/lib/visfd/feature_implementation.hpp +++ b/lib/visfd/feature_implementation.hpp @@ -31,566 +31,62 @@ namespace visfd { -/// @brief The following function finds EITHER local minima OR local maxima -/// in an image. +/// @brief Figure out the score for blobs located at each position in crds[]. +/// Terminology: Every blob has a position, diameter, and a score. +/// The j'th "blob" is a sphere, centerd at blob_crds[j], with diameter +/// blob_diameters[j]. Associated with every blob is a "score" +/// (a number) stored in blob_scores[j]. +/// If crds[i] lies inside one of the blobs, store the corresponding +/// score for that blob in scores[i]. (If crds[i] lies inside more +/// than one sphere, give priority spheres occuring later in the list.) +/// If crds[i] lies outside any of the spheres, store -infinity there. /// -/// For this version, the location where each minima or maxima occurs (ix,iy,iz) -/// is represented by an integer (called an "index") which equals: -/// @code -/// index = ix + iy*image_size[0] + iz*image_size[0]*image_size[1] -/// @endcode -/// where ix,iy,iz are the coordinates of the corresponding voxel in the image, -/// and image_size[] stores the size of the image in the x,y,z directions. -/// If pv_minima_indices!=nullptr, then *pv_minima_indices will store a list -/// of the indices corresponding to the locations of the local minima. -/// If pv_maxima_indices!=nullptr, then *pv_maxima_indices will store a list -/// of the indices corresponding to the locations of the local maxima. -/// The corresponding voxel intensities (brightness values) will be stored in -/// *pv_minima_scores and *pv_maxima_scores (assuming they are != nullptr). -/// Thresholds can be used to discard minima or maxima whose corresponding -/// voxel intensities are not sufficiently low or high, respectively. -/// If the aaafMask[][][] is not equal to nullptr, then local minima and maxima -/// will be ignored if the corresponding entry in aaafMask[][][] equals 0. +/// @returns void. Results are stored in "scores". /// -/// @note: THIS VERSION OF THE FUNCTION WAS NOT INTENDED FOR PUBLIC USE. +/// @note THIS FUNCTION WAS NOT INTENDED FOR PUBLIC USE -template +template static void -_FindExtrema(int const image_size[3], //!< size of the image in x,y,z directions - Scalar const *const *const *aaafSource, //!< image array aaafSource[iz][iy][ix] - Scalar const *const *const *aaafMask, //!< optional: ignore voxel ix,iy,iz if aaafMask!=nullptr and aaafMask[iz][iy][ix]==0 - vector *pv_minima_indices, //!< a list of integers uniquely identifying the location of each minima - vector *pv_maxima_indices, //!< a list of integers uniquely identifying the location of each maxima - vector *pv_minima_scores, //!< store corresponding minima aaafSource[iz][iy][ix] values here (if not nullptr) - vector *pv_maxima_scores, //!< store corresponding maxima aaafSource[iz][iy][ix] values here (if not nullptr) - vector *pv_minima_nvoxels, //!< store number of voxels in each minima (usually 1) - vector *pv_maxima_nvoxels, //!< store number of voxels in each maxima (usually 1) - Scalar minima_threshold = std::numeric_limits::infinity(), // Ignore minima which are not sufficiently low - Scalar maxima_threshold = -std::numeric_limits::infinity(), // Ignore maxima which are not sufficiently 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? +_FindBlobScores(const vector >& crds, //!< locations of blob-like things we are looking for + vector& scores, //!< stores the score of the blob (sphere) contains that position + vector& sphere_ids, //!< which blob (sphere) contains this position? + const vector >& blob_crds, //!< location of the center of each spherical blob (sorted in order of increasing priority) + const vector& blob_diameters, //!< diameter of each blob + const vector& blob_scores, //!< "score" of each blob (a number) + SortCriteria sort_blob_criteria = SORT_DECREASING_MAGNITUDE, //!< give priority to high or low scoring blobs? + ostream *pReportProgress = nullptr //!< report progress back to the user? + ) { - assert(aaafSource); - - bool find_minima = (pv_minima_indices != nullptr); - bool find_maxima = (pv_maxima_indices != nullptr); - - vector minima_indices; //store minima indices here - vector maxima_indices; //store maxima indices here - vector minima_scores; //store the brightness of the minima here - vector maxima_scores; //store the brightness of the maxima here - vector minima_nvoxels; //store number of voxels in each minima here - vector maxima_nvoxels; //store number of voxels in each maxima here - if (pv_minima_indices == nullptr) - pv_minima_indices = &minima_indices; - if (pv_maxima_indices == nullptr) - pv_maxima_indices = &maxima_indices; - if (pv_minima_scores == nullptr) - pv_minima_scores = &minima_scores; - if (pv_maxima_scores == nullptr) - pv_maxima_scores = &maxima_scores; - if (pv_minima_nvoxels == nullptr) - pv_minima_nvoxels = &minima_nvoxels; - if (pv_maxima_nvoxels == nullptr) - pv_maxima_nvoxels = &maxima_nvoxels; - - - ptrdiff_t ***aaaiExtrema; // 3-D array storing the plateau to which each voxel belongs (if any) - ptrdiff_t *aiExtrema; // the same array, stored contiguously in 1-D - Alloc3D(image_size, - &aiExtrema, - &aaaiExtrema); - - // We will assign voxels in aaaiExtrema[][][] to the following values: - ptrdiff_t NEITHER = 0;//(means this voxel is neither a local minima or maxima) - ptrdiff_t UNDEFINED = image_size[0] * image_size[1] * image_size[2] + 1; - ptrdiff_t QUEUED = image_size[0] * image_size[1] * image_size[2] + 2; - // ... or we will assign the voxel to an integer denoting which - // local minima or maxima it belongs to. - //Note: Both "UNDEFINED" and "QUEUED" are impossible values (in this respect). - // There cannot be that many local maxima or minima in the image. - - - for (int iz = 0; iz < image_size[2]; iz++) - for (int iy = 0; iy < image_size[1]; iy++) - for (int ix = 0; ix < image_size[0]; ix++) - aaaiExtrema[iz][iy][ix] = UNDEFINED; - - - if (pReportProgress) - *pReportProgress - << " -- Attempting to allocate space for one more image. --\n" - << " -- (If this crashes your computer, find a computer with --\n" - << " -- more RAM and use \"ulimit\", OR use a smaller image.) --\n" - << "\n"; - - if (pReportProgress) - *pReportProgress << "---- searching for local minima & maxima ----\n"; - - // Figure out which neighbors to consider when searching neighboring voxels - int (*neighbors)[3] = nullptr; //a pointer to a fixed-length array of 3 ints - int num_neighbors = 0; - { - // How big is the search neighborhood around each minima? - int r_neigh_sqd = connectivity; - int r_neigh = floor(sqrt(r_neigh_sqd)); - - vector > vNeighbors; - for (int jz = -r_neigh; jz <= r_neigh; jz++) { - for (int jy = -r_neigh; jy <= r_neigh; jy++) { - for (int jx = -r_neigh; jx <= r_neigh; jx++) { - array j_xyz; - j_xyz[0] = jx; - j_xyz[1] = jy; - j_xyz[2] = jz; - if ((jx == 0) && (jy == 0) && (jz == 0)) - continue; - else if (jx*jx+jy*jy+jz*jz > r_neigh_sqd) - continue; - else - vNeighbors.push_back(j_xyz); - } - } - } - // Convert the list of neighbors vNeigh to a contiguous 2D C-style array: - num_neighbors = vNeighbors.size(); - neighbors = new int[num_neighbors][3]; - for (int j = 0; j < num_neighbors; j++) - for (int d = 0; d < 3; d++) - neighbors[j][d] = vNeighbors[j][d]; - // We will use neighbors[][] from now on.. - } // ...done figuring out the list of voxel neighbors - - - // Now search the image for minima, maxima. - // Do not assume that the local minima or maxima are point-like. - // They could be "plateaus". - // So at every voxel location ix,iy,iz, search for neighboring voxels - // with identical intensities (brightness). These voxels collectively are - // a "plateau". Find the neighbors of these voxels to determine whether - // this "plateau" is also a local minima or maxima. - for (int iz0 = 0; iz0 < image_size[2]; iz0++) { - - if (pReportProgress) - *pReportProgress << " searching for extrema at z="< > visited; - //#endif - - - // It's possible that a local minima or maxima consists of more than one - // adjacent voxels with identical intensities (brightnesses). - // These voxels belong to a "plateau" of voxels of the same height - // which are above (or below) all of the surrounding voxels. - // (Note: Most of the time this plateau consists of only one voxel.) - // The "q_plateau" variable is a data structure to keep track of which - // voxels belong to the current plateau. - queue > q_plateau; - array i_xyz0; - i_xyz0[0] = ix0;//note:this can be shortened to 1 line with make_array() - i_xyz0[1] = iy0; - i_xyz0[2] = iz0; - q_plateau.push(i_xyz0); - size_t n_plateau = 0; // the number of voxels in this plateau - - // We also need a reverse lookup table to check whether a given voxel - // is in the queue (or has already been previously assigned). - // That's what aaaiExtrema[][][] is. - // If aaaiExtrema[iz][iy][ix] == QUEUED, then we know it's in the queue. - // We need to come up with an impossible value for QUEUED: - aaaiExtrema[iz0][iy0][ix0] = QUEUED; //make sure we don't visit it twice - - while (! q_plateau.empty()) - { - array p = q_plateau.front(); - q_plateau.pop(); - int ix = p[0]; - int iy = p[1]; - int iz = p[2]; - n_plateau++; - assert(aaaiExtrema[iz][iy][ix] == QUEUED); - - - //#ifndef NDEBUG - //assert(visited.find(p) == visited.end()); - //visited.insert(p); - //#endif - - - for (int j = 0; j < num_neighbors; j++) { - int jx = neighbors[j][0]; - int jy = neighbors[j][1]; - int jz = neighbors[j][2]; - int iz_jz = iz + jz; - int iy_jy = iy + jy; - int ix_jx = ix + jx; - if (((iz_jz < 0) || (image_size[2] <= iz_jz)) - || - ((iy_jy < 0) || (image_size[1] <= iy_jy)) - || - ((ix_jx < 0) || (image_size[0] <= ix_jx)) - || - (aaafMask && (aaafMask[iz_jz][iy_jy][ix_jx] == 0.0))) - { - if (! allow_borders) { - is_minima = false; - is_maxima = false; - } - continue; - } - - - - if (aaafSource[iz_jz][iy_jy][ix_jx] == aaafSource[iz][iy][ix]) - { - if (aaaiExtrema[iz_jz][iy_jy][ix_jx] == UNDEFINED)//don't visit twice - { - // then add this voxel to the q_plateau - array ij_xyz; - ij_xyz[0] = ix_jx; - ij_xyz[1] = iy_jy; - ij_xyz[2] = iz_jz; - q_plateau.push(ij_xyz); - aaaiExtrema[iz_jz][iy_jy][ix_jx] = QUEUED; - - //#ifndef NDEBUG - //assert(visited.find(ij_xyz) == visited.end()); - //#endif - } - } - else { - if (aaafSource[iz+jz][iy+jy][ix+jx] < aaafSource[iz][iy][ix]) - is_minima = false; - else if (aaafSource[iz+jz][iy+jy][ix+jx] > aaafSource[iz][iy][ix]) - is_maxima = false; - else - assert(false); - } - } //for (int j = 0; j < num_neighbors; j++) - } // while (! q_plateau.empty()) - - - // If this voxel is either a minima or a maxima, add it to the list. - if (is_minima && find_minima) { - if (((! aaafMask) || (aaafMask[iz0][iy0][ix0] != 0)) && - ((aaafSource[iz0][iy0][ix0] <= minima_threshold) - // || (maxima_threshold < minima_threshold) - )) - { - // convert from a 3D location (ix,iy,iz) to a 1D index ("index") - IntegerIndex index = ix0 + image_size[0]*(iy0 + image_size[1]*iz0); - //minima_crds.push_back(ixiyiz); - pv_minima_indices->push_back(index); - pv_minima_scores->push_back(aaafSource[iz0][iy0][ix0]); - pv_minima_nvoxels->push_back(n_plateau); - } - } - if (is_maxima && find_maxima) { - if (((! aaafMask) || (aaafMask[iz0][iy0][ix0] != 0)) && - ((aaafSource[iz0][iy0][ix0] >= maxima_threshold) - // || (maxima_threshold < minima_threshold) - )) - { - // convert from a 3D location (ix,iy,iz) to a 1D index ("index") - IntegerIndex index = ix0 + image_size[0]*(iy0 + image_size[1]*iz0); - //maxima_crds.push_back(ixiyiz); - pv_maxima_indices->push_back(index); - pv_maxima_scores->push_back(aaafSource[iz0][iy0][ix0]); - pv_maxima_nvoxels->push_back(n_plateau); - } - } - - - - //#ifndef NDEBUG - //visited.clear(); - //#endif - - - - //now loop through those same voxels again to reset the aaaiExtrema array - - // What to store in aaaiExtrema[][][] for voxels in this plateau? - ptrdiff_t assigned_plateau_label; - if (is_maxima) - assigned_plateau_label = pv_maxima_scores->size(); - else if (is_minima) - assigned_plateau_label = -pv_minima_scores->size(); - else - assigned_plateau_label = NEITHER; - aaaiExtrema[iz0][iy0][ix0] = assigned_plateau_label; - - - assert(q_plateau.empty()); - q_plateau.push(i_xyz0); - - - while (! q_plateau.empty()) - { - array p = q_plateau.front(); - q_plateau.pop(); - int ix = p[0]; - int iy = p[1]; - int iz = p[2]; - assert(aaaiExtrema[iz][iy][ix] != QUEUED); - - //#ifndef NDEBUG - //assert(visited.find(p) == visited.end()); - //visited.insert(p); - //#endif - - - for (int j = 0; j < num_neighbors; j++) { - int jx = neighbors[j][0]; - int jy = neighbors[j][1]; - int jz = neighbors[j][2]; - int iz_jz = iz + jz; - int iy_jy = iy + jy; - int ix_jx = ix + jx; - if (((iz_jz < 0) || (image_size[2] <= iz_jz)) - || - ((iy_jy < 0) || (image_size[1] <= iy_jy)) - || - ((ix_jx < 0) || (image_size[0] <= ix_jx)) - || - (aaafMask && (aaafMask[iz_jz][iy_jy][ix_jx] == 0.0))) - continue; - if (aaaiExtrema[iz_jz][iy_jy][ix_jx] == QUEUED) //don't visit twice - { - assert(aaafSource[iz_jz][iy_jy][ix_jx] == aaafSource[iz][iy][ix]); - - // then add this voxel to the q_plateau - array ij_xyz; - ij_xyz[0] = ix_jx; - ij_xyz[1] = iy_jy; - ij_xyz[2] = iz_jz; - q_plateau.push(ij_xyz); - - // Commenting out: - //assert(! (is_minima && is_maxima)); <- wrong - // Actually this is possible. For example in an image where all - // the voxels are identical, all voxels are BOTH minima a maxima. - // This is such a rare, pathelogical case, I don't worry about it. - // (The aaaiExtrema[][][] array is only used to prevent visiting - // the same voxel twice. It doesn't matter what is stored there - // as long as it fulfills this role.) - - aaaiExtrema[iz_jz][iy_jy][ix_jx] = assigned_plateau_label; - } - } //for (int j = 0; j < num_neighbors; j++) - } // while (! q_plateau.empty()) - } // for (int ix0 = 0; ix0 < image_size[0]; ix0++) - } // for (int iy0 = 0; iy0 < image_size[1]; iy0++) - } // for (int iz0 = 0; iz0 < image_size[2]; iz0++) - - - // Sort the minima and maxima in increasing and decreasing order, respectively - IntegerIndex n_minima, n_maxima; - bool descending_order; - - if (pv_minima_indices) { - // Sort the minima in increasing order - n_minima = pv_minima_indices->size(); - if (n_minima > 0) { - descending_order = false; - vector > score_index(n_minima); - for (IntegerIndex i = 0; i < n_minima; i++) - score_index[i] = make_tuple((*pv_minima_scores)[i], i); - if (pReportProgress) - *pReportProgress << "-- Sorting minima according to their scores... "; - if (descending_order) - sort(score_index.rbegin(), - score_index.rend()); - else - sort(score_index.begin(), - score_index.end()); - vector permutation(n_minima); - for (IntegerIndex i = 0; i < score_index.size(); i++) - permutation[i] = get<1>(score_index[i]); - score_index.clear(); - apply_permutation(permutation, *pv_minima_indices); - apply_permutation(permutation, *pv_minima_scores); - apply_permutation(permutation, *pv_minima_nvoxels); - // Optional: The minima in the image are not in sorted order either. Fix? - vector perm_inv; - invert_permutation(permutation, perm_inv); - for (int iz = 0; iz < image_size[2]; iz++) { - for (int iy = 0; iy < image_size[1]; iy++) { - for (int ix = 0; ix < image_size[0]; ix++) { - if (aaafMask && (aaafMask[iz][iy][ix] == 0)) - continue; - if (aaaiExtrema[iz][iy][ix] < 0) - aaaiExtrema[iz][iy][ix]=-perm_inv[(-aaaiExtrema[iz][iy][ix])-1]-1; - } - } - } - if (pReportProgress) - *pReportProgress << "done --" << endl; - } - } - - if (pv_maxima_indices) { - // Sort the maxima in decreasing order - n_maxima = pv_maxima_indices->size(); - if (n_maxima > 0) { - descending_order = true; - vector > score_index(n_maxima); - for (IntegerIndex i = 0; i < n_maxima; i++) - score_index[i] = make_tuple((*pv_maxima_scores)[i], i); - if (pReportProgress) - *pReportProgress << "-- Sorting maxima according to their scores... "; - if (descending_order) - sort(score_index.rbegin(), - score_index.rend()); - else - sort(score_index.begin(), - score_index.end()); - vector permutation(n_maxima); - for (IntegerIndex i = 0; i < score_index.size(); i++) - permutation[i] = get<1>(score_index[i]); - score_index.clear(); - apply_permutation(permutation, *pv_maxima_indices); - apply_permutation(permutation, *pv_maxima_scores); - apply_permutation(permutation, *pv_maxima_nvoxels); - if (pReportProgress) - *pReportProgress << "done --" << endl; - // Optional: The maxima in the image are not in sorted order either. Fix? - vector perm_inv; - invert_permutation(permutation, perm_inv); - for (int iz = 0; iz < image_size[2]; iz++) { - for (int iy = 0; iy < image_size[1]; iy++) { - for (int ix = 0; ix < image_size[0]; ix++) { - if (aaafMask && (aaafMask[iz][iy][ix] == 0)) - continue; - if (aaaiExtrema[iz][iy][ix] > 0) - aaaiExtrema[iz][iy][ix] = perm_inv[aaaiExtrema[iz][iy][ix]-1]+1; - } - } - } - } + // Sort the blobs by score in order of increasing priority + vector > blob_crds_sorted = blob_crds; + vector blob_diameters_sorted = blob_diameters; + vector blob_scores_sorted = blob_scores; + + SortBlobs(blob_crds_sorted, + blob_diameters_sorted, + blob_scores_sorted, + sort_blob_criteria, + true, + nullptr, + pReportProgress); + + // Figure out which sphere is located at each position + FindSpheres(crds, + sphere_ids, + blob_crds_sorted, + blob_diameters_sorted, + pReportProgress); + + assert(sphere_ids.size() == crds.size()); + scores.resize(sphere_ids.size(), + -std::numeric_limits::infinity()); //<-impossible score + const size_t UNOCCUPIED = 0; + for (size_t i=0; i < sphere_ids.size(); i++) { + size_t which_blob = sphere_ids[i]; + if (which_blob != UNOCCUPIED) + scores[i] = blob_scores_sorted[which_blob-1]; } - - if (aaaiDest) { - for (int iz = 0; iz < image_size[2]; iz++) { - for (int iy = 0; iy < image_size[1]; iy++) { - for (int ix = 0; ix < image_size[0]; ix++) { - if (aaafMask && (aaafMask[iz][iy][ix] == 0.0)) - continue; // don't modify voxels outside the mask - aaaiDest[iz][iy][ix] = aaaiExtrema[iz][iy][ix]; - // Until now, minima in the image are represented by negative integers - // and maxima are represented by positive integers. - if (((! find_minima) || (! find_maxima)) && - (aaaiExtrema[iz][iy][ix] < 0)) - // If we only asked for minima OR maxima (not both) - // then just use positive integers only. - aaaiExtrema[iz][iy][ix] = -aaaiExtrema[iz][iy][ix]; - } - } - } - } - - Dealloc3D(image_size, - &aiExtrema, - &aaaiExtrema); - - delete [] neighbors; - -} // _FindExtrema() - - - - -/// @brief The following function finds EITHER local minima OR local maxima -/// in an image. -/// -/// For this version, the location where each minima or maxima occurs (ix,iy,iz) -/// is represented by an integer (called an "index") which equals: -/// @code -/// index = ix + iy*image_size[0] + iz*image_size[0]*image_size[1] -/// @endcode -/// where ix,iy,iz are the coordinates of the corresponding voxel in the image. -/// -/// @note: THIS VERSION OF THE FUNCTION WAS NOT INTENDED FOR PUBLIC USE. - -template -static void -_FindExtrema(int const image_size[3], //!< size of the image in x,y,z directions - Scalar const *const *const *aaafI, //!< image array aaafI[iz][iy][ix] - Scalar const *const *const *aaafMask, //!< optional: ignore voxel ix,iy,iz if aaafMask!=nullptr and aaafMask[iz][iy][ix]==0 - vector &extrema_indices, //!< a list of integers uniquely identifying the location of each minima or maxima - vector &extrema_scores, //!< corresponding voxel brightness at that location - vector &extrema_nvoxels, //!< how many voxels belong to each minima or maxima? - bool seek_minima=true, //!< search for minima or maxima? - 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? -{ - // NOTE: - // C++ will not allow us to supply nullptr to a function that expects a pointer - // to a template expression: Template argument deduction/substitution fails. - // We need to re-cast "nullptr" as a pointer with the correct type. - // One way to do that is to define these new versions of nullptr: - vector *null_vI = nullptr; - vector *null_vi = nullptr; - vector *null_vf = nullptr; - if (seek_minima) { - _FindExtrema(image_size, - aaafI, - aaafMask, - &extrema_indices, // store maxima locations here - null_vI, // <-- don't search for maxima - &extrema_scores, // store minima values here - null_vf, // <-- don't search for maxima - &extrema_nvoxels, // store number of voxels in each minima - null_vi, // <-- don't search for maxima - threshold, - -std::numeric_limits::infinity(), - connectivity, - allow_borders, - aaaiDest, - pReportProgress); - } - else { - if (threshold == std::numeric_limits::infinity()) - threshold = -std::numeric_limits::infinity(); - _FindExtrema(image_size, - aaafI, - aaafMask, - null_vI, // <-- don't search for minima_crds, - &extrema_indices, // store maxima locations here - null_vf, // <-- don't search for minima_scores, - &extrema_scores, // store maxima values here - null_vi, // <-- don't search for minima - &extrema_nvoxels, // store number of voxels in each maxima - std::numeric_limits::infinity(), - threshold, - connectivity, - allow_borders, - aaaiDest, - pReportProgress); - } - -} //_FindExtrema() - +} //_FindBlobScores() @@ -819,7 +315,7 @@ _ChooseBlobScoreThresholds(const vector >& blob_crds, //!< locat const vector& training_accepted, //!< classify each blob as "accepted" (true) or "rejected" (false) Scalar *pthreshold_lower_bound = nullptr, //!< return threshold to the caller Scalar *pthreshold_upper_bound = nullptr, //!< return threshold to the caller - SortBlobCriteria sort_blob_criteria = PRIORITIZE_HIGH_MAGNITUDE_SCORES, //!< give priority to high or low scoring blobs? + SortCriteria sort_blob_criteria = SORT_DECREASING_MAGNITUDE, //!< give priority to high or low scoring blobs? ostream *pReportProgress = nullptr //!< report progress back to the user? ) { @@ -886,7 +382,7 @@ _ChooseBlobScoreThresholdsMulti(const vector > >& blob_cr const vector >& training_accepted, //!< classify each blob as "accepted" (true) or "rejected" (false) Scalar *pthreshold_lower_bound = nullptr, //!< return threshold to the caller Scalar *pthreshold_upper_bound = nullptr, //!< return threshold to the caller - SortBlobCriteria sort_blob_criteria = PRIORITIZE_HIGH_MAGNITUDE_SCORES, //!< give priority to high or low scoring blobs? + SortCriteria sort_blob_criteria = SORT_DECREASING_MAGNITUDE, //!< give priority to high or low scoring blobs? ostream *pReportProgress = nullptr //!< report progress back to the user? ) { diff --git a/lib/visfd/filter2d.hpp b/lib/visfd/filter2d.hpp index 2dbc55a..83d8e00 100644 --- a/lib/visfd/filter2d.hpp +++ b/lib/visfd/filter2d.hpp @@ -427,130 +427,6 @@ GenFilterGenGauss2D(Scalar width[2], //!< "s_x", "s_y" parameters - -/// @brief -/// Create a 2D filter and fill it with a difference of (generalized) Gaussians: -/// This version requires that the caller has already created individual -/// filters for the two gaussians. -/// All this function does is subtract one filter from the other (and rescale). - -template - -static Filter2D -_GenFilterDogg2D(Scalar width_a[2], //!< "a" parameter in formula - Scalar width_b[2], //!< "b" parameter in formula - Scalar m_exp, //!< "m" parameter in formula - Scalar n_exp, //!< "n" parameter in formula - Filter2D& filterXY_A, //!< filters for the two - Filter2D& filterXY_B, //!< gaussians - Scalar *pA=nullptr, //!< optional:report A,B coeffs to user - Scalar *pB=nullptr, //!< optional:report A,B coeffs to user - ostream *pReportProgress = nullptr //!< optional:report filter details to the user? - ) -{ - - Scalar A, B; - //A, B = height of the central peak - A = filterXY_A.aafH[0][0]; - B = filterXY_B.aafH[0][0]; - - // The "difference of gaussians" filter is the difference between - // these two (generalized) gaussian filters. - int halfwidth[2]; - halfwidth[0] = std::max(filterXY_A.halfwidth[0], filterXY_B.halfwidth[0]); - halfwidth[1] = std::max(filterXY_A.halfwidth[1], filterXY_B.halfwidth[1]); - Filter2D filter(halfwidth); - - if (pReportProgress) - *pReportProgress << "Array of 2D filter entries:" << endl; - - for (int iy=-halfwidth[1]; iy<=halfwidth[1]; iy++) { - for (int ix=-halfwidth[0]; ix<=halfwidth[0]; ix++) { - filter.aafH[iy][ix] = 0.0; - - // The two filters may have different widths, so we have to check - // that ix and iy lie within the domain of these two filters before - // adding or subtracting their values from the final GDOG filter. - if (((-filterXY_A.halfwidth[0]<=ix) && (ix<=filterXY_A.halfwidth[0])) && - ((-filterXY_A.halfwidth[1]<=iy) && (iy<=filterXY_A.halfwidth[1]))) - - filter.aafH[iy][ix] +=filterXY_A.aafH[iy][ix]; // /(A-B); COMMENTING OUT - - // COMMENTING OUT: (The factor of 1/(A-B) insures that the central peak has height 1) - - if (((-filterXY_B.halfwidth[0]<=ix) && (ix<=filterXY_B.halfwidth[0])) && - ((-filterXY_B.halfwidth[1]<=iy) && (iy<=filterXY_B.halfwidth[1]))) - - filter.aafH[iy][ix] -=filterXY_B.aafH[iy][ix]; // /(A-B); COMMENTING OUT - - - - //FOR DEBUGGING REMOVE EVENTUALLY - //if (pReportProgress) - // *pReportProgress << "GenDogg2D: aafH["< -// Create a 2D filter and fill it with a difference of (generalized) Gaussians: -Filter2D -GenFilterDogg2D(Scalar width_a[2], //"a" parameter in formula - Scalar width_b[2], //"b" parameter in formula - Scalar m_exp, //"m" parameter in formula - Scalar n_exp, //"n" parameter in formula - int halfwidth[2], - Scalar *pA = nullptr, //optional:report A,B coeffs to user - Scalar *pB = nullptr, //optional:report A,B coeffs to user - ostream *pReportProgress = nullptr) -{ - Filter2D filterXY_A = - GenFilterGenGauss2D(width_a, //"a_x", "a_y" gaussian width parameters - m_exp, //"n" exponent parameter - halfwidth); - //pReportProgress); - - Filter2D filterXY_B = - GenFilterGenGauss2D(width_b, //"b_x", "b_y" gaussian width parameters - n_exp, //"n" exponent parameter - halfwidth); - //pReportProgress); - - return _GenFilterDogg2D(width_a, - width_b, - m_exp, - n_exp, - filterXY_A, filterXY_B, - pA, - pB, - pReportProgress); -} //GenFilterDogg2D(...halfwidth...) - - - } //namespace visfd diff --git a/lib/visfd/filter3d.hpp b/lib/visfd/filter3d.hpp index b11ea5c..54ff4f6 100644 --- a/lib/visfd/filter3d.hpp +++ b/lib/visfd/filter3d.hpp @@ -7,21 +7,14 @@ #define _FILTER3D_HPP #include +#include #include #include #include -#include -#include -#include -#include -#include using namespace std; -#include // defines the "VisfdErr" exception type #include // defines Alloc3D() and Dealloc3D() #include // defines "Filter1D" (used in ApplySeparable()) -#include // defines invert_permutation(), AveArray(), ... -#include // defines matrix diagonalizer (DiagonalizeSym3()) -#include // defines DotProduct3(),CrossProduct(),quaternions... +#include // defines the "VisfdErr" exception type @@ -400,7 +393,7 @@ class Filter3D { /// @param pDenominator=if you want to store the sum of the weights considered, pass a pointer to a number /// (useful if the sum was not complete due to some voxels being masked out, /// or because the filter extends beyond the boundaries of the image) - /// @note: THIS FUNCTION WAS NOT INTENDED FOR PUBLIC USE. + /// @note: This is a private function, not intended for public use. Scalar ApplyToVoxel(Integer ix, Integer iy, @@ -1866,180 +1859,6 @@ LocalFluctuationsByRadius(Integer const image_size[3], //!< number of voxels in - -// --------------- MISCELLANEOUS ODDS AND ENDS --------------- - - - -/// @brief Create a 3D filter and fill it with a difference of (generalized) Gaussians -/// -/// This version requires that the caller has already created individual -/// filters for the two gaussians. -/// All this function does is subtract one filter from the other (and rescale). -/// @note: DEPRECIATION WARNING: It's not clear if this filter is useful. -/// I may delete this function in the future. -/// @note: THIS FUNCTION WAS NOT INTENDED FOR PUBLIC USE - -template - -static -Filter3D -_GenFilterDogg3D(Scalar width_a[3], //!< "a" parameter in formula - Scalar width_b[3], //!< "b" parameter in formula - Scalar m_exp, //!< "m" parameter in formula - Scalar n_exp, //!< "n" parameter in formula - Filter3D& filter_A, //!< filters for the two - Filter3D& filter_B, //!< gaussians - Scalar *pA=nullptr, //!< optional:report A,B coeffs to user - Scalar *pB=nullptr, //!< optional:report A,B coeffs to user - ostream *pReportEquation = nullptr //!< optional: report equation to the user - ) -{ - Scalar A, B; - //A, B = height of the central peak - A = filter_A.aaafH[0][0][0]; - B = filter_B.aaafH[0][0][0]; - - - // The "difference of gaussians" filter is the difference between - // these two (generalized) gaussian filters. - int halfwidth[3]; - halfwidth[0] = std::max(filter_A.halfwidth[0], filter_B.halfwidth[0]); - halfwidth[1] = std::max(filter_A.halfwidth[1], filter_B.halfwidth[1]); - halfwidth[2] = std::max(filter_A.halfwidth[2], filter_B.halfwidth[2]); - Filter3D filter(halfwidth); - - //FOR DEBUGGING REMOVE EVENTUALLY - //if (pReportEquation) - // *pReportEquation << "Array of 3D filter entries:" << endl; - - - for (int iz=-halfwidth[2]; iz<=halfwidth[2]; iz++) { - for (int iy=-halfwidth[1]; iy<=halfwidth[1]; iy++) { - for (int ix=-halfwidth[0]; ix<=halfwidth[0]; ix++) { - - filter.aaafH[iz][iy][ix] = 0.0; - - // The two filters may have different widths, so we have to check - // that ix,iy and iz lie within the domain of these two filters before - // adding or subtracting their values from the final GDoG filter. - if (((-filter_A.halfwidth[0]<=ix) && (ix<=filter_A.halfwidth[0])) && - ((-filter_A.halfwidth[1]<=iy) && (iy<=filter_A.halfwidth[1])) && - ((-filter_A.halfwidth[2]<=iz) && (iz<=filter_A.halfwidth[2]))) - - filter.aaafH[iz][iy][ix] += - filter_A.aaafH[iz][iy][ix]; // /(A-B); COMMENTING OUT - - - - // COMMENTING OUT: (The factor of 1/(A-B) insures that the central peak has height 1) - - - if (((-filter_B.halfwidth[0]<=ix) && (ix<=filter_B.halfwidth[0])) && - ((-filter_B.halfwidth[1]<=iy) && (iy<=filter_B.halfwidth[1])) && - ((-filter_B.halfwidth[2]<=iz) && (iz<=filter_B.halfwidth[2]))) - - filter.aaafH[iz][iy][ix] -= - filter_B.aaafH[iz][iy][ix]; // /(A-B); COMMENTING OUT - - //*pReportEquation << aaafH[iz][iy][ix]; - // - //if (ix == 0) pReportEquation << "\n"; else pReportEquation << " "; - - } // for (int ix=-halfwidth[0]; ix<=halfwidth[0]; ix++) { - } // for (int iy=-halfwidth[1]; iy<=halfwidth[1]; iy++) { - } // for (int iz=-halfwidth[2]; iz<=halfwidth[2]; iz++) { - - - // COMMENTING OUT the factor of 1/(A-B): - //A = A/(A-B); - //B = B/(A-B); - - if (pA && pB) { - *pA = A; // Rescale A and B numbers returned to the caller - *pB = B; // (because we divided the array entries by (A-B) earlier) - } - - if (pReportEquation) { - *pReportEquation << "\n"; - *pReportEquation << " Filter Used:\n" - " h(x,y,z) = h_a(x,y,z) - h_b(x,y,z)\n" - " h_a(x,y,z) = A*exp(-((x/a_x)^2 + (y/a_y)^2 + (z/a_z)^2)^(m/2))\n" - " h_b(x,y,z) = B*exp(-((x/b_x)^2 + (y/b_y)^2 + (z/b_z)^2)^(n/2))\n" - " ... where A = " << A << "\n" - " B = " << B << "\n" - " m = " << m_exp << "\n" - " n = " << n_exp << "\n" - " (a_x, a_y, a_z) = " - << "(" << width_a[0] - << " " << width_a[1] - << " " << width_a[2] << ")\n" - " (b_x, b_y, b_z) = " - << "(" << width_b[0] - << " " << width_b[1] - << " " << width_b[2] << ")\n"; - *pReportEquation << " You can plot a slice of this function\n" - << " in the X direction using:\n" - " draw_filter_1D.py -dogg " << A << " " << B - << " " << width_a[0] << " " << width_b[0] - << " " << m_exp << " " << n_exp << endl; - *pReportEquation << " and in the Y direction using:\n" - " draw_filter_1D.py -dogg " << A << " " << B - << " " << width_a[1] << " " << width_b[1] - << " " << m_exp << " " << n_exp << endl; - *pReportEquation << " and in the Z direction using:\n" - " draw_filter_1D.py -dogg " << A << " " << B - << " " << width_a[2] << " " << width_b[2] - << " " << m_exp << " " << n_exp << endl; - } //if (pReportEquation) - - return filter; -} //_GenFilterDogg3D() - - - -/// @brief Create a 3D filter and fill it with a difference of (generalized) Gaussians -/// @verbatim h(x,y,z) = A*exp(-(r/a)^m) - B*exp(-(r/b)^n) @endverbatim -/// where @verbatim r = sqrt(x^2 + y^2 + z^2) @endverbatim -/// and "A" and "B" are determined by normalization of each term independently -/// DEPRECIATION WARNING: It's not clear if this type if filter is useful. -/// I may delete this function in the future. - -template - -Filter3D -GenFilterDogg3D(Scalar width_a[3], //!< "a" parameter in formula - Scalar width_b[3], //!< "b" parameter in formula - Scalar m_exp, //!< "m" parameter in formula - Scalar n_exp, //!< "n" parameter in formula - int halfwidth[3], //!< the width of the filter - Scalar *pA=nullptr, //!< optional:report A,B coeffs to user - Scalar *pB=nullptr, //!< optional:report A,B coeffs to user - ostream *pReportEquation = nullptr //!< optional: print params used? - ) -{ - Filter3D filter_A = - GenFilterGenGauss3D(width_a, //"a_x", "a_y" gaussian width parameters - m_exp, //"m" exponent parameter - halfwidth); - - Filter3D filter_B = - GenFilterGenGauss3D(width_b, //"b_x", "b_y" gaussian width parameters - n_exp, //"n" exponent parameter - halfwidth); - - return _GenFilterDogg3D(width_a, - width_b, - m_exp, - n_exp, - filter_A, filter_B, - pA, - pB, - pReportEquation); -} //GenFilterDogg3D(...halfwidth...) - - - } //namespace visfd diff --git a/lib/visfd/morphology.hpp b/lib/visfd/morphology.hpp new file mode 100644 index 0000000..7129fa3 --- /dev/null +++ b/lib/visfd/morphology.hpp @@ -0,0 +1,226 @@ +/// @file morphology.hpp +/// @brief Functions relevant to image morphology. +/// @author Andrew Jewett +/// @date 2021-7-07 + +#ifndef _MORPHOLOGY_HPP +#define _MORPHOLOGY_HPP + + +#include +#include +#include +#include +#include +#include +using namespace std; +#include // defines the "VisfdErr" exception type +#include // defines matrix diagonalizer (DiagonalizeSym3()) +#include // defines invert_permutation(), FindSpheres() +#include // defines Alloc2D() and Dealloc2D() +#include // defines Alloc3D() and Dealloc3D() +#include // defines "Filter1D" (used in ApplySeparable()) +#include // defines common 3D image filters +#include + + + + +namespace visfd { + + + +/// @brief Find all of the local minima and local maxima in an image (aaafI). +/// The locations of minima and maxima are stored in the +/// *pv_minima_crds and *pv_maxima_crds arrays, and sorted in +/// increasing and decreasing order respectively. (IE they are sorted +/// so that the most significant local minima or maxima will appear +/// first in the list.) +/// If either pv_minima_crds or pv_maxima_crds is nullptr, then +/// the minima or maxima will be ignored. +/// The optional pv_minima_scores and pv_maxima_scores store the +/// The caller can automatically discard minima or maxima which +/// are not sufficiently low or high, by supplying thresholds. +/// The optional aaafMask array (if not nullptr) can be used to ignore +/// certain voxels in the image (whose aaafMask entries are zero). +/// @note Local minima or maxima on the boundaries of the image +/// (or near the edge of the mask) +/// are not as trustworthy since some of the neighboring voxels +/// will not be available for comparison. These minima and maxima +/// can be ignored by setting allow_borders=false. The number of +/// neighbors around every voxel which are considered (eg, 6, 18, 26) +/// can be controlled using the "connectivity" argument. + +template + +void +FindExtrema(int const image_size[3], //!< size of the image in x,y,z directions + Scalar const *const *const *aaafI, //!< image array aaafI[iz][iy][ix] + Scalar const *const *const *aaafMask, //!< optional: ignore voxel ix,iy,iz if aaafMask!=nullptr and aaafMask[iz][iy][ix]==0 + vector > *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) + vector *pv_maxima_scores, //!< store corresponding maxima aaafI[iz][iy][ix] values here (if not nullptr) + vector *pv_minima_nvoxels, //!< store number of voxels in each minima (usually 1) + vector *pv_maxima_nvoxels, //!< store number of voxels in each maxima (usually 1) + Scalar minima_threshold=std::numeric_limits::infinity(), //!< ignore minima with intensities greater than this + Scalar maxima_threshold=-std::numeric_limits::infinity(), //!< ignore maxima with intensities lessr 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? + ) +{ + bool find_minima = (pva_minima_crds != nullptr); + bool find_maxima = (pva_maxima_crds != nullptr); + + vector minima_indices; + vector maxima_indices; + vector minima_scores; + vector maxima_scores; + vector *pv_minima_indices = nullptr; + vector *pv_maxima_indices = nullptr; + if (find_minima) { + pv_minima_indices = &minima_indices; + if (! pv_minima_scores) + pv_minima_scores = &minima_scores; + } + if (find_maxima) { + pv_maxima_indices = &maxima_indices; + if (! pv_maxima_scores) + pv_maxima_scores = &maxima_scores; + } + + // This function is defined in visfd_utils_implementation.hpp + _FindExtrema(image_size, + aaafI, + aaafMask, + pv_minima_indices, + pv_maxima_indices, + pv_minima_scores, + pv_maxima_scores, + pv_minima_nvoxels, + pv_maxima_nvoxels, + minima_threshold, + maxima_threshold, + connectivity, + allow_borders, + aaaiDest, + pReportProgress); + + // Now convert the indices back to x,y,z coordinates + if (pva_minima_crds) { + size_t N = minima_indices.size(); + assert(pva_minima_crds); + pva_minima_crds->resize(N); + for (size_t n = 0; n < N; n++) { + size_t i = minima_indices[n]; + // convert from a 1D index (i) to 3-D indices (ix, iy, iz) + size_t ix = i % image_size[0]; + i /= image_size[0]; + size_t iy = i % image_size[1]; + i /= image_size[1]; + size_t iz = i; + (*pva_minima_crds)[n][0] = ix; + (*pva_minima_crds)[n][1] = iy; + (*pva_minima_crds)[n][2] = iz; + } + } + if (pva_maxima_crds) { + size_t N = maxima_indices.size(); + assert(pva_maxima_crds); + pva_maxima_crds->resize(N); + for (size_t n = 0; n < N; n++) { + size_t i = maxima_indices[n]; + // convert from a 1D index (i) to 3-D indices (ix, iy, iz) + size_t ix = i % image_size[0]; + i /= image_size[0]; + size_t iy = i % image_size[1]; + i /= image_size[1]; + size_t iz = i; + (*pva_maxima_crds)[n][0] = ix; + (*pva_maxima_crds)[n][1] = iy; + (*pva_maxima_crds)[n][2] = iz; + } + } +} //FindExtrema() + + + + + +/// @brief The following version of this function seeks either minima +/// or maxima, but not both. (If you want both, use the other version. +/// That version is faster than invoking this function twice.) +/// See the description of that version for details. + +template + +void +FindExtrema(int const image_size[3], //!< size of input image array + Scalar const *const *const *aaafI, //!< input image array + Scalar const *const *const *aaafMask, //!< if not nullptr then zero entries indicate which voxels to ignore + vector > &extrema_crds, //!< store the location of each extrema + vector &extrema_scores, //!< store the brightness of each extrema + vector &extrema_nvoxels, //!< store number of voxels in each extrema (usually 1) + bool seek_minima=true, //!< search for minima or maxima? + 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? +{ + // NOTE: + // C++ will not allow us to supply nullptr to a function that expects a pointer + // to a template expression: Template argument deduction/substitution fails. + // We need to re-cast "nullptr" as a pointer with the correct type. + // One way to do that is to define these new versions of nullptr: + vector > *null_vai3 = nullptr; + vector *null_vf = nullptr; + vector *null_vi = nullptr; + + if (seek_minima) { + FindExtrema(image_size, + aaafI, + aaafMask, + &extrema_crds, // store minima locations here + null_vai3, // <-- don't search for maxima + &extrema_scores, // store minima values here + null_vf, // <-- don't search for maxima + &extrema_nvoxels, // store number of voxels in each minima + null_vi, // <-- don't search for maxima + threshold, + -std::numeric_limits::infinity(), + connectivity, + allow_borders, + aaaiDest, + pReportProgress); + } + else { + if (threshold == std::numeric_limits::infinity()) + threshold = -std::numeric_limits::infinity(); + FindExtrema(image_size, + aaafI, + aaafMask, + null_vai3, // <-- don't search for minima + &extrema_crds, // store maxima locations here + null_vf, // <-- don't search for minima + &extrema_scores, // store maxima values here + null_vi, // <-- don't search for minima + &extrema_nvoxels, // store number of voxels in each maxima + std::numeric_limits::infinity(), + threshold, + connectivity, + allow_borders, + aaaiDest, + pReportProgress); + } +} // FindExtrema() + + + +} //namespace visfd + + + +#endif //#ifndef _MORPHOLOGY_HPP diff --git a/lib/visfd/morphology_implementation.hpp b/lib/visfd/morphology_implementation.hpp new file mode 100644 index 0000000..0f2462c --- /dev/null +++ b/lib/visfd/morphology_implementation.hpp @@ -0,0 +1,596 @@ +/// @file morphology_implementation.hpp +/// @brief IMPLEMENTATION DETAILS. THIS CODE IS NOT INTENDED FOR PUBLIC USE +/// (The functions defined here are invoked by "morphology.hpp") +/// @author Andrew Jewett +/// @date 2021-7-07 + +#ifndef _MORPHOLOGY_IMPLEMENTATION_HPP +#define _MORPHOLOGY_IMPLEMENTATION_HPP + +#include +#include +#include +#include +#include +#include +#include +#include +#include +using namespace std; +#include // defines the "VisfdErr" exception type +#include // defines invert_permutation(), FindSpheres() +#include // defines Alloc2D() and Dealloc2D() +#include // defines Alloc3D() and Dealloc3D() +#include // defines "Filter1D" (used in ApplySeparable()) +#include // defines common 3D image filters + + +namespace visfd { + +/// @brief The following function finds EITHER local minima OR local maxima +/// in an image. +/// +/// For this version, the location where each minima or maxima occurs (ix,iy,iz) +/// is represented by an integer (called an "index") which equals: +/// @code +/// index = ix + iy*image_size[0] + iz*image_size[0]*image_size[1] +/// @endcode +/// where ix,iy,iz are the coordinates of the corresponding voxel in the image, +/// and image_size[] stores the size of the image in the x,y,z directions. +/// If pv_minima_indices!=nullptr, then *pv_minima_indices will store a list +/// of the indices corresponding to the locations of the local minima. +/// If pv_maxima_indices!=nullptr, then *pv_maxima_indices will store a list +/// of the indices corresponding to the locations of the local maxima. +/// The corresponding voxel intensities (brightness values) will be stored in +/// *pv_minima_scores and *pv_maxima_scores (assuming they are != nullptr). +/// Thresholds can be used to discard minima or maxima whose corresponding +/// voxel intensities are not sufficiently low or high, respectively. +/// If the aaafMask[][][] is not equal to nullptr, then local minima and maxima +/// will be ignored if the corresponding entry in aaafMask[][][] equals 0. +/// +/// @note: THIS VERSION OF THE FUNCTION WAS NOT INTENDED FOR PUBLIC USE. + +template +static void +_FindExtrema(int const image_size[3], //!< size of the image in x,y,z directions + Scalar const *const *const *aaafSource, //!< image array aaafSource[iz][iy][ix] + Scalar const *const *const *aaafMask, //!< optional: ignore voxel ix,iy,iz if aaafMask!=nullptr and aaafMask[iz][iy][ix]==0 + vector *pv_minima_indices, //!< a list of integers uniquely identifying the location of each minima + vector *pv_maxima_indices, //!< a list of integers uniquely identifying the location of each maxima + vector *pv_minima_scores, //!< store corresponding minima aaafSource[iz][iy][ix] values here (if not nullptr) + vector *pv_maxima_scores, //!< store corresponding maxima aaafSource[iz][iy][ix] values here (if not nullptr) + vector *pv_minima_nvoxels, //!< store number of voxels in each minima (usually 1) + vector *pv_maxima_nvoxels, //!< store number of voxels in each maxima (usually 1) + Scalar minima_threshold = std::numeric_limits::infinity(), // Ignore minima which are not sufficiently low + Scalar maxima_threshold = -std::numeric_limits::infinity(), // Ignore maxima which are not sufficiently 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? +{ + assert(aaafSource); + + bool find_minima = (pv_minima_indices != nullptr); + bool find_maxima = (pv_maxima_indices != nullptr); + + vector minima_indices; //store minima indices here + vector maxima_indices; //store maxima indices here + vector minima_scores; //store the brightness of the minima here + vector maxima_scores; //store the brightness of the maxima here + vector minima_nvoxels; //store number of voxels in each minima here + vector maxima_nvoxels; //store number of voxels in each maxima here + if (pv_minima_indices == nullptr) + pv_minima_indices = &minima_indices; + if (pv_maxima_indices == nullptr) + pv_maxima_indices = &maxima_indices; + if (pv_minima_scores == nullptr) + pv_minima_scores = &minima_scores; + if (pv_maxima_scores == nullptr) + pv_maxima_scores = &maxima_scores; + if (pv_minima_nvoxels == nullptr) + pv_minima_nvoxels = &minima_nvoxels; + if (pv_maxima_nvoxels == nullptr) + pv_maxima_nvoxels = &maxima_nvoxels; + + + ptrdiff_t ***aaaiExtrema; // 3-D array storing the plateau to which each voxel belongs (if any) + ptrdiff_t *aiExtrema; // the same array, stored contiguously in 1-D + Alloc3D(image_size, + &aiExtrema, + &aaaiExtrema); + + // We will assign voxels in aaaiExtrema[][][] to the following values: + ptrdiff_t NEITHER = 0;//(means this voxel is neither a local minima or maxima) + ptrdiff_t UNDEFINED = image_size[0] * image_size[1] * image_size[2] + 1; + ptrdiff_t QUEUED = image_size[0] * image_size[1] * image_size[2] + 2; + // ... or we will assign the voxel to an integer denoting which + // local minima or maxima it belongs to. + //Note: Both "UNDEFINED" and "QUEUED" are impossible values (in this respect). + // There cannot be that many local maxima or minima in the image. + + + for (int iz = 0; iz < image_size[2]; iz++) + for (int iy = 0; iy < image_size[1]; iy++) + for (int ix = 0; ix < image_size[0]; ix++) + aaaiExtrema[iz][iy][ix] = UNDEFINED; + + + if (pReportProgress) + *pReportProgress + << " -- Attempting to allocate space for one more image. --\n" + << " -- (If this crashes your computer, find a computer with --\n" + << " -- more RAM and use \"ulimit\", OR use a smaller image.) --\n" + << "\n"; + + if (pReportProgress) + *pReportProgress << "---- searching for local minima & maxima ----\n"; + + // Figure out which neighbors to consider when searching neighboring voxels + int (*neighbors)[3] = nullptr; //a pointer to a fixed-length array of 3 ints + int num_neighbors = 0; + { + // How big is the search neighborhood around each minima? + int r_neigh_sqd = connectivity; + int r_neigh = floor(sqrt(r_neigh_sqd)); + + vector > vNeighbors; + for (int jz = -r_neigh; jz <= r_neigh; jz++) { + for (int jy = -r_neigh; jy <= r_neigh; jy++) { + for (int jx = -r_neigh; jx <= r_neigh; jx++) { + array j_xyz; + j_xyz[0] = jx; + j_xyz[1] = jy; + j_xyz[2] = jz; + if ((jx == 0) && (jy == 0) && (jz == 0)) + continue; + else if (jx*jx+jy*jy+jz*jz > r_neigh_sqd) + continue; + else + vNeighbors.push_back(j_xyz); + } + } + } + // Convert the list of neighbors vNeigh to a contiguous 2D C-style array: + num_neighbors = vNeighbors.size(); + neighbors = new int[num_neighbors][3]; + for (int j = 0; j < num_neighbors; j++) + for (int d = 0; d < 3; d++) + neighbors[j][d] = vNeighbors[j][d]; + // We will use neighbors[][] from now on.. + } // ...done figuring out the list of voxel neighbors + + + // Now search the image for minima, maxima. + // Do not assume that the local minima or maxima are point-like. + // They could be "plateaus". + // So at every voxel location ix,iy,iz, search for neighboring voxels + // with identical intensities (brightness). These voxels collectively are + // a "plateau". Find the neighbors of these voxels to determine whether + // this "plateau" is also a local minima or maxima. + for (int iz0 = 0; iz0 < image_size[2]; iz0++) { + + if (pReportProgress) + *pReportProgress << " searching for extrema at z="< > visited; + //#endif + + + // It's possible that a local minima or maxima consists of more than one + // adjacent voxels with identical intensities (brightnesses). + // These voxels belong to a "plateau" of voxels of the same height + // which are above (or below) all of the surrounding voxels. + // (Note: Most of the time this plateau consists of only one voxel.) + // The "q_plateau" variable is a data structure to keep track of which + // voxels belong to the current plateau. + queue > q_plateau; + array i_xyz0; + i_xyz0[0] = ix0;//note:this can be shortened to 1 line with make_array() + i_xyz0[1] = iy0; + i_xyz0[2] = iz0; + q_plateau.push(i_xyz0); + size_t n_plateau = 0; // the number of voxels in this plateau + + // We also need a reverse lookup table to check whether a given voxel + // is in the queue (or has already been previously assigned). + // That's what aaaiExtrema[][][] is. + // If aaaiExtrema[iz][iy][ix] == QUEUED, then we know it's in the queue. + // We need to come up with an impossible value for QUEUED: + aaaiExtrema[iz0][iy0][ix0] = QUEUED; //make sure we don't visit it twice + + while (! q_plateau.empty()) + { + array p = q_plateau.front(); + q_plateau.pop(); + int ix = p[0]; + int iy = p[1]; + int iz = p[2]; + n_plateau++; + assert(aaaiExtrema[iz][iy][ix] == QUEUED); + + + //#ifndef NDEBUG + //assert(visited.find(p) == visited.end()); + //visited.insert(p); + //#endif + + + for (int j = 0; j < num_neighbors; j++) { + int jx = neighbors[j][0]; + int jy = neighbors[j][1]; + int jz = neighbors[j][2]; + int iz_jz = iz + jz; + int iy_jy = iy + jy; + int ix_jx = ix + jx; + if (((iz_jz < 0) || (image_size[2] <= iz_jz)) + || + ((iy_jy < 0) || (image_size[1] <= iy_jy)) + || + ((ix_jx < 0) || (image_size[0] <= ix_jx)) + || + (aaafMask && (aaafMask[iz_jz][iy_jy][ix_jx] == 0.0))) + { + if (! allow_borders) { + is_minima = false; + is_maxima = false; + } + continue; + } + + + + if (aaafSource[iz_jz][iy_jy][ix_jx] == aaafSource[iz][iy][ix]) + { + if (aaaiExtrema[iz_jz][iy_jy][ix_jx] == UNDEFINED)//don't visit twice + { + // then add this voxel to the q_plateau + array ij_xyz; + ij_xyz[0] = ix_jx; + ij_xyz[1] = iy_jy; + ij_xyz[2] = iz_jz; + q_plateau.push(ij_xyz); + aaaiExtrema[iz_jz][iy_jy][ix_jx] = QUEUED; + + //#ifndef NDEBUG + //assert(visited.find(ij_xyz) == visited.end()); + //#endif + } + } + else { + if (aaafSource[iz+jz][iy+jy][ix+jx] < aaafSource[iz][iy][ix]) + is_minima = false; + else if (aaafSource[iz+jz][iy+jy][ix+jx] > aaafSource[iz][iy][ix]) + is_maxima = false; + else + assert(false); + } + } //for (int j = 0; j < num_neighbors; j++) + } // while (! q_plateau.empty()) + + + // If this voxel is either a minima or a maxima, add it to the list. + if (is_minima && find_minima) { + if (((! aaafMask) || (aaafMask[iz0][iy0][ix0] != 0)) && + ((aaafSource[iz0][iy0][ix0] <= minima_threshold) + // || (maxima_threshold < minima_threshold) + )) + { + // convert from a 3D location (ix,iy,iz) to a 1D index ("index") + IntegerIndex index = ix0 + image_size[0]*(iy0 + image_size[1]*iz0); + //minima_crds.push_back(ixiyiz); + pv_minima_indices->push_back(index); + pv_minima_scores->push_back(aaafSource[iz0][iy0][ix0]); + pv_minima_nvoxels->push_back(n_plateau); + } + } + if (is_maxima && find_maxima) { + if (((! aaafMask) || (aaafMask[iz0][iy0][ix0] != 0)) && + ((aaafSource[iz0][iy0][ix0] >= maxima_threshold) + // || (maxima_threshold < minima_threshold) + )) + { + // convert from a 3D location (ix,iy,iz) to a 1D index ("index") + IntegerIndex index = ix0 + image_size[0]*(iy0 + image_size[1]*iz0); + //maxima_crds.push_back(ixiyiz); + pv_maxima_indices->push_back(index); + pv_maxima_scores->push_back(aaafSource[iz0][iy0][ix0]); + pv_maxima_nvoxels->push_back(n_plateau); + } + } + + + + //#ifndef NDEBUG + //visited.clear(); + //#endif + + + + //now loop through those same voxels again to reset the aaaiExtrema array + + // What to store in aaaiExtrema[][][] for voxels in this plateau? + ptrdiff_t assigned_plateau_label; + if (is_maxima) + assigned_plateau_label = pv_maxima_scores->size(); + else if (is_minima) + assigned_plateau_label = -pv_minima_scores->size(); + else + assigned_plateau_label = NEITHER; + aaaiExtrema[iz0][iy0][ix0] = assigned_plateau_label; + + + assert(q_plateau.empty()); + q_plateau.push(i_xyz0); + + + while (! q_plateau.empty()) + { + array p = q_plateau.front(); + q_plateau.pop(); + int ix = p[0]; + int iy = p[1]; + int iz = p[2]; + assert(aaaiExtrema[iz][iy][ix] != QUEUED); + + //#ifndef NDEBUG + //assert(visited.find(p) == visited.end()); + //visited.insert(p); + //#endif + + + for (int j = 0; j < num_neighbors; j++) { + int jx = neighbors[j][0]; + int jy = neighbors[j][1]; + int jz = neighbors[j][2]; + int iz_jz = iz + jz; + int iy_jy = iy + jy; + int ix_jx = ix + jx; + if (((iz_jz < 0) || (image_size[2] <= iz_jz)) + || + ((iy_jy < 0) || (image_size[1] <= iy_jy)) + || + ((ix_jx < 0) || (image_size[0] <= ix_jx)) + || + (aaafMask && (aaafMask[iz_jz][iy_jy][ix_jx] == 0.0))) + continue; + if (aaaiExtrema[iz_jz][iy_jy][ix_jx] == QUEUED) //don't visit twice + { + assert(aaafSource[iz_jz][iy_jy][ix_jx] == aaafSource[iz][iy][ix]); + + // then add this voxel to the q_plateau + array ij_xyz; + ij_xyz[0] = ix_jx; + ij_xyz[1] = iy_jy; + ij_xyz[2] = iz_jz; + q_plateau.push(ij_xyz); + + // Commenting out: + //assert(! (is_minima && is_maxima)); <- wrong + // Actually this is possible. For example in an image where all + // the voxels are identical, all voxels are BOTH minima a maxima. + // This is such a rare, pathelogical case, I don't worry about it. + // (The aaaiExtrema[][][] array is only used to prevent visiting + // the same voxel twice. It doesn't matter what is stored there + // as long as it fulfills this role.) + + aaaiExtrema[iz_jz][iy_jy][ix_jx] = assigned_plateau_label; + } + } //for (int j = 0; j < num_neighbors; j++) + } // while (! q_plateau.empty()) + } // for (int ix0 = 0; ix0 < image_size[0]; ix0++) + } // for (int iy0 = 0; iy0 < image_size[1]; iy0++) + } // for (int iz0 = 0; iz0 < image_size[2]; iz0++) + + + // Sort the minima and maxima in increasing and decreasing order, respectively + IntegerIndex n_minima, n_maxima; + bool descending_order; + + if (pv_minima_indices) { + // Sort the minima in increasing order + n_minima = pv_minima_indices->size(); + if (n_minima > 0) { + descending_order = false; + vector > score_index(n_minima); + for (IntegerIndex i = 0; i < n_minima; i++) + score_index[i] = make_tuple((*pv_minima_scores)[i], i); + if (pReportProgress) + *pReportProgress << "-- Sorting minima according to their scores... "; + if (descending_order) + sort(score_index.rbegin(), + score_index.rend()); + else + sort(score_index.begin(), + score_index.end()); + vector permutation(n_minima); + for (IntegerIndex i = 0; i < score_index.size(); i++) + permutation[i] = get<1>(score_index[i]); + score_index.clear(); + apply_permutation(permutation, *pv_minima_indices); + apply_permutation(permutation, *pv_minima_scores); + apply_permutation(permutation, *pv_minima_nvoxels); + // Optional: The minima in the image are not in sorted order either. Fix? + vector perm_inv; + invert_permutation(permutation, perm_inv); + for (int iz = 0; iz < image_size[2]; iz++) { + for (int iy = 0; iy < image_size[1]; iy++) { + for (int ix = 0; ix < image_size[0]; ix++) { + if (aaafMask && (aaafMask[iz][iy][ix] == 0)) + continue; + if (aaaiExtrema[iz][iy][ix] < 0) + aaaiExtrema[iz][iy][ix]=-perm_inv[(-aaaiExtrema[iz][iy][ix])-1]-1; + } + } + } + if (pReportProgress) + *pReportProgress << "done --" << endl; + } + } + + if (pv_maxima_indices) { + // Sort the maxima in decreasing order + n_maxima = pv_maxima_indices->size(); + if (n_maxima > 0) { + descending_order = true; + vector > score_index(n_maxima); + for (IntegerIndex i = 0; i < n_maxima; i++) + score_index[i] = make_tuple((*pv_maxima_scores)[i], i); + if (pReportProgress) + *pReportProgress << "-- Sorting maxima according to their scores... "; + if (descending_order) + sort(score_index.rbegin(), + score_index.rend()); + else + sort(score_index.begin(), + score_index.end()); + vector permutation(n_maxima); + for (IntegerIndex i = 0; i < score_index.size(); i++) + permutation[i] = get<1>(score_index[i]); + score_index.clear(); + apply_permutation(permutation, *pv_maxima_indices); + apply_permutation(permutation, *pv_maxima_scores); + apply_permutation(permutation, *pv_maxima_nvoxels); + if (pReportProgress) + *pReportProgress << "done --" << endl; + // Optional: The maxima in the image are not in sorted order either. Fix? + vector perm_inv; + invert_permutation(permutation, perm_inv); + for (int iz = 0; iz < image_size[2]; iz++) { + for (int iy = 0; iy < image_size[1]; iy++) { + for (int ix = 0; ix < image_size[0]; ix++) { + if (aaafMask && (aaafMask[iz][iy][ix] == 0)) + continue; + if (aaaiExtrema[iz][iy][ix] > 0) + aaaiExtrema[iz][iy][ix] = perm_inv[aaaiExtrema[iz][iy][ix]-1]+1; + } + } + } + } + } + + if (aaaiDest) { + for (int iz = 0; iz < image_size[2]; iz++) { + for (int iy = 0; iy < image_size[1]; iy++) { + for (int ix = 0; ix < image_size[0]; ix++) { + if (aaafMask && (aaafMask[iz][iy][ix] == 0.0)) + continue; // don't modify voxels outside the mask + aaaiDest[iz][iy][ix] = aaaiExtrema[iz][iy][ix]; + // Until now, minima in the image are represented by negative integers + // and maxima are represented by positive integers. + if (((! find_minima) || (! find_maxima)) && + (aaaiExtrema[iz][iy][ix] < 0)) + // If we only asked for minima OR maxima (not both) + // then just use positive integers only. + aaaiExtrema[iz][iy][ix] = -aaaiExtrema[iz][iy][ix]; + } + } + } + } + + Dealloc3D(image_size, + &aiExtrema, + &aaaiExtrema); + + delete [] neighbors; + +} // _FindExtrema() + + + + +/// @brief The following function finds EITHER local minima OR local maxima +/// in an image. +/// +/// For this version, the location where each minima or maxima occurs (ix,iy,iz) +/// is represented by an integer (called an "index") which equals: +/// @code +/// index = ix + iy*image_size[0] + iz*image_size[0]*image_size[1] +/// @endcode +/// where ix,iy,iz are the coordinates of the corresponding voxel in the image. +/// +/// @note: THIS VERSION OF THE FUNCTION WAS NOT INTENDED FOR PUBLIC USE. + +template +static void +_FindExtrema(int const image_size[3], //!< size of the image in x,y,z directions + Scalar const *const *const *aaafI, //!< image array aaafI[iz][iy][ix] + Scalar const *const *const *aaafMask, //!< optional: ignore voxel ix,iy,iz if aaafMask!=nullptr and aaafMask[iz][iy][ix]==0 + vector &extrema_indices, //!< a list of integers uniquely identifying the location of each minima or maxima + vector &extrema_scores, //!< corresponding voxel brightness at that location + vector &extrema_nvoxels, //!< how many voxels belong to each minima or maxima? + bool seek_minima=true, //!< search for minima or maxima? + 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? +{ + // NOTE: + // C++ will not allow us to supply nullptr to a function that expects a pointer + // to a template expression: Template argument deduction/substitution fails. + // We need to re-cast "nullptr" as a pointer with the correct type. + // One way to do that is to define these new versions of nullptr: + vector *null_vI = nullptr; + vector *null_vi = nullptr; + vector *null_vf = nullptr; + if (seek_minima) { + _FindExtrema(image_size, + aaafI, + aaafMask, + &extrema_indices, // store maxima locations here + null_vI, // <-- don't search for maxima + &extrema_scores, // store minima values here + null_vf, // <-- don't search for maxima + &extrema_nvoxels, // store number of voxels in each minima + null_vi, // <-- don't search for maxima + threshold, + -std::numeric_limits::infinity(), + connectivity, + allow_borders, + aaaiDest, + pReportProgress); + } + else { + if (threshold == std::numeric_limits::infinity()) + threshold = -std::numeric_limits::infinity(); + _FindExtrema(image_size, + aaafI, + aaafMask, + null_vI, // <-- don't search for minima_crds, + &extrema_indices, // store maxima locations here + null_vf, // <-- don't search for minima_scores, + &extrema_scores, // store maxima values here + null_vi, // <-- don't search for minima + &extrema_nvoxels, // store number of voxels in each maxima + std::numeric_limits::infinity(), + threshold, + connectivity, + allow_borders, + aaaiDest, + pReportProgress); + } + +} //_FindExtrema() + + + +} //namespace visfd + + + +#endif //#ifndef _MORPHOLOGY_IMPLEMENTATION_HPP diff --git a/lib/visfd/segmentation.hpp b/lib/visfd/segmentation.hpp index 98d2546..0ced42e 100644 --- a/lib/visfd/segmentation.hpp +++ b/lib/visfd/segmentation.hpp @@ -24,7 +24,7 @@ using namespace std; #include // defines Alloc3D() and Dealloc3D() #include // defines "Filter1D" (used in ApplySeparable()) #include // defines common 3D image filters - +#include // defines FindExtrema() namespace visfd { diff --git a/lib/visfd/visfd.hpp b/lib/visfd/visfd.hpp index 59c6df9..529a959 100644 --- a/lib/visfd/visfd.hpp +++ b/lib/visfd/visfd.hpp @@ -9,12 +9,13 @@ // higher level functions and classes are defined in these header files: -#include // simple feature detectors (blob, surface,...) +#include // feature detectors (blob, surface, edge, curve..) #include // functions used for segmentation (eg watershed) +#include // functions relevant to image morphology #include // functions for clustering nearby voxels (useful - // to separate different objects in an image) + // to distinguish different objects in an image) #include // functions for drawing and annotation -#include // common 3D filters +#include // common 3D (linear) filters #include // resize the image // lower level functions and classes (upon which the files above depend): diff --git a/lib/visfd/visfd_utils.hpp b/lib/visfd/visfd_utils.hpp index 988c632..4fb1c15 100644 --- a/lib/visfd/visfd_utils.hpp +++ b/lib/visfd/visfd_utils.hpp @@ -37,16 +37,16 @@ static int SGN(Scalar val) { -/// @brief invert a permutation -template -void -invert_permutation(const vector& p, - vector& p_inv) -{ - p_inv.resize(p.size()); - for (Integer i = 0; i < p.size(); i++) - p_inv[p[i]] = i; -} +/// @brief A variable of type "SortCriteria" is often passed as an +/// argument to any function that needs to sort lists of numbers. +typedef enum eSortCriteria { + DO_NOT_SORT, + SORT_DECREASING, + SORT_INCREASING, + SORT_DECREASING_MAGNITUDE, + SORT_INCREASING_MAGNITUDE +} SortCriteria; + @@ -67,6 +67,49 @@ apply_permutation(const vector& p, } +/// @brief invert a permutation +template +void +invert_permutation(const vector& p, + vector& p_inv) +{ + p_inv.resize(p.size()); + for (Integer i = 0; i < p.size(); i++) + p_inv[p[i]] = i; +} + + + +/// @brief Calculate the volume of overlap between two spheres of radius +/// Ri and Rj separated by a distance of rij. + +template + +Scalar CalcSphereOverlap(Scalar rij,//! Rj) { + Scalar tmp = Ri; + Ri = Rj; + Rj = tmp; + } + if (rij <= Ri) { + return (4*M_PI/3) * Ri*Ri*Ri; + } + + // "xi" and "xj" are the distances from the sphere centers + // to the plane where the two spheres intersect. + Scalar xi = 0.5 * (1.0/rij) * (rij*rij + Ri*Ri - Rj*Rj); + Scalar xj = 0.5 * (1.0/rij) * (rij*rij + Rj*Rj - Ri*Ri); + Scalar volume_overlap = + (M_PI/3)*( Ri*Ri*Ri * (2 - (xi/Ri) * (3 - SQR(xi/Ri))) + + Rj*Rj*Rj * (2 - (xj/Rj) * (3 - SQR(xj/Rj))) ); + return volume_overlap; +} + /// @brief This function was intended to be used as a way to allow end-users @@ -195,124 +238,6 @@ FindNearestVoxel(int const image_size[3], //!< #voxels in xyz -/// @brief sort blobs by their scores - -template - -void -SortBlobs(vector >& blob_crds,//!< x,y,z of each blob's center - vector& blob_diameters, //!< the width of each blob - vector& blob_scores, //!< the score for each blob - bool ascending_order = true, - bool ignore_score_sign = true, - vector *pPermutation = nullptr, //!< optional: return the new sorted order to the caller - ostream *pReportProgress = nullptr //!< optional: report progress to the user? - ) -{ - size_t n_blobs = blob_crds.size(); - assert(n_blobs == blob_diameters.size()); - assert(n_blobs == blob_scores.size()); - vector > score_index(n_blobs); - for (size_t i = 0; i < n_blobs; i++) { - if (ignore_score_sign) - score_index[i] = make_tuple(std::fabs(blob_scores[i]), i); - else - score_index[i] = make_tuple(blob_scores[i], i); - } - - if (n_blobs == 0) - return; - - if (pReportProgress) - *pReportProgress << "-- Sorting blobs according to their scores... "; - if (ascending_order) - sort(score_index.begin(), - score_index.end()); - else - sort(score_index.rbegin(), - score_index.rend()); - - vector permutation(n_blobs); - for (size_t i = 0; i < score_index.size(); i++) - permutation[i] = get<1>(score_index[i]); - - if (pPermutation != nullptr) - *pPermutation = permutation; - - score_index.clear(); - apply_permutation(permutation, blob_crds); - apply_permutation(permutation, blob_diameters); - apply_permutation(permutation, blob_scores); - if (pReportProgress) - *pReportProgress << "done --" << endl; - -} //SortBlobs() - - - -/// @brief these options tell "DiscardOverlappingBlobs" how to give -/// priority to different blobs based on their score (ie, minima or maxima?) -typedef enum eSortBlobCriteria { - DO_NOT_SORT, - PRIORITIZE_HIGH_SCORES, - PRIORITIZE_LOW_SCORES, - PRIORITIZE_HIGH_MAGNITUDE_SCORES, - PRIORITIZE_LOW_MAGNITUDE_SCORES -} SortBlobCriteria; - - -/// @brief sort blobs by their scores according to "sort_blob_criteria" - -template - -void -SortBlobs(vector >& blob_crds,//!< x,y,z of each blob's center - vector& blob_diameters, //!< the width of each blob - vector& blob_scores, //!< the score for each blob - SortBlobCriteria sort_blob_criteria, //!< give priority to high or low scoring blobs? - bool ascending_order = true, - vector *pPermutation = nullptr, //!< optional: return the new sorted order to the caller - ostream *pReportProgress = nullptr //!< optional: report progress to the user? - ) -{ - if (sort_blob_criteria == PRIORITIZE_HIGH_SCORES) - SortBlobs(blob_crds, - blob_diameters, - blob_scores, - ascending_order, - false, - pPermutation, - pReportProgress); - else if (sort_blob_criteria == PRIORITIZE_LOW_SCORES) - SortBlobs(blob_crds, - blob_diameters, - blob_scores, - ! ascending_order, - false, - pPermutation, - pReportProgress); - else if (sort_blob_criteria == PRIORITIZE_HIGH_MAGNITUDE_SCORES) - SortBlobs(blob_crds, - blob_diameters, - blob_scores, - ascending_order, - true, - pPermutation, - pReportProgress); - else if (sort_blob_criteria == PRIORITIZE_LOW_MAGNITUDE_SCORES) - SortBlobs(blob_crds, - blob_diameters, - blob_scores, - ! ascending_order, - true, - pPermutation, - pReportProgress); -} //SortBlobs() - - - - - /// @brief Suppose we have already detected a list of objects ("blobs") in /// an image, and recorded their locations and sizes @@ -430,131 +355,6 @@ FindSpheres(const vector >& crds, //!< locations of blob-like th -/// @brief Figure out the score for blobs located at each position in crds[]. -/// Terminology: Every blob has a position, diameter, and a score. -/// The j'th "blob" is a sphere, centerd at blob_crds[j], with diameter -/// blob_diameters[j]. Associated with every blob is a "score" -/// (a number) stored in blob_scores[j]. -/// If crds[i] lies inside one of the blobs, store the corresponding -/// score for that blob in scores[i]. (If crds[i] lies inside more -/// than one sphere, give priority spheres occuring later in the list.) -/// If crds[i] lies outside any of the spheres, store -infinity there. -/// -/// @returns void. Results are stored in "scores". -/// -/// @note THIS FUNCTION WAS NOT INTENDED FOR PUBLIC USE - -template -static void -_FindBlobScores(const vector >& crds, //!< locations of blob-like things we are looking for - vector& scores, //!< stores the score of the blob (sphere) contains that position - vector& sphere_ids, //!< which blob (sphere) contains this position? - const vector >& blob_crds, //!< location of the center of each spherical blob (sorted in order of increasing priority) - const vector& blob_diameters, //!< diameter of each blob - const vector& blob_scores, //!< "score" of each blob (a number) - SortBlobCriteria sort_blob_criteria = PRIORITIZE_HIGH_MAGNITUDE_SCORES, //!< give priority to high or low scoring blobs? - ostream *pReportProgress = nullptr //!< report progress back to the user? - ) -{ - // Sort the blobs by score in order of increasing priority - vector > blob_crds_sorted = blob_crds; - vector blob_diameters_sorted = blob_diameters; - vector blob_scores_sorted = blob_scores; - - SortBlobs(blob_crds_sorted, - blob_diameters_sorted, - blob_scores_sorted, - sort_blob_criteria, - true, - nullptr, - pReportProgress); - - // Figure out which sphere is located at each position - FindSpheres(crds, - sphere_ids, - blob_crds_sorted, - blob_diameters_sorted, - pReportProgress); - - assert(sphere_ids.size() == crds.size()); - scores.resize(sphere_ids.size(), - -std::numeric_limits::infinity()); //<-impossible score - const size_t UNOCCUPIED = 0; - for (size_t i=0; i < sphere_ids.size(); i++) { - size_t which_blob = sphere_ids[i]; - if (which_blob != UNOCCUPIED) - scores[i] = blob_scores_sorted[which_blob-1]; - } -} //_FindBlobScores() - - - -/// @brief This variant of the function also takes care of discarding -/// training data which is not sufficiently close to any of the blobs. -/// It also throws an exception if the remaining training data is empty. -/// @return This function has no return value. -/// The results are stored in: -/// out_training_crds, -/// out_training_accepted, -/// out_training_scores - -template -static void - -FindBlobScores(const vector >& in_training_crds, //!< locations of blob-like things in training set - const vector& in_training_accepted, //!< classify each blob-like thing as "accepted" (true) or "rejected" (false) - vector >& out_training_crds, //!< same as in_training_crds after discarding entries too far away from any existing blobs - vector& out_training_accepted, //!< same as in_training_accepted after discarding entries too far away from any existing blobs - vector& out_training_scores, //!< store the scores of the remaining blob-like things here - const vector >& blob_crds, //!< location of each blob (in voxels, sorted by score in increasing priority) - const vector& blob_diameters, //!< diameger of each blob (sorted by score in increasing priority) - const vector& blob_scores, //!< priority of each blob (sorted by score in increasing priority) - SortBlobCriteria sort_blob_criteria = PRIORITIZE_HIGH_MAGNITUDE_SCORES, //!< give priority to high or low scoring blobs? - ostream *pReportProgress = nullptr //!< report progress back to the user? - ) -{ - // Figure out the score for each blob: - // in_training_crds[i] == position of ith training set data - // in_training_which_blob[i] == which blob contains this position? - // in_training_score[i] = score of this blob - - vector in_training_which_blob; // which blob (sphere) contains this position? - vector in_training_scores; // what is the score of that blob? - - _FindBlobScores(in_training_crds, - in_training_scores, - in_training_which_blob, - blob_crds, - blob_diameters, - blob_scores, - sort_blob_criteria, - pReportProgress); - - size_t _N = in_training_crds.size(); - assert(_N == in_training_scores.size()); - assert(_N == in_training_accepted.size()); - - - // Discard training data that was not suffiently close to one of the blobs. - // Store the remaining training set data in the following variables: - const size_t UNOCCUPIED = 0; - for (size_t i=0; i < _N; i++) { - if (in_training_which_blob[i] != UNOCCUPIED) { - out_training_crds.push_back(in_training_crds[i]); - out_training_scores.push_back(in_training_scores[i]); - out_training_accepted.push_back(in_training_accepted[i]); - } - } - - // N = the number of training data left after discarding datum which - // do not lie within any blobs. (not within any blob radii) - size_t N = out_training_crds.size(); - assert(N == out_training_scores.size()); - assert(N == out_training_accepted.size()); - -} //FindBlobScores() - - // @brief This function chooses the "score" threshold which maximizes the // number of correctly labelled training data.