From 56bc1e9a7dad269e65f97f07930e98969ed1b756 Mon Sep 17 00:00:00 2001 From: Andrew Jewett Date: Tue, 27 Jul 2021 02:42:55 -0700 Subject: [PATCH] renamed "WatershedDirectional()" to "LabelConnected()", and moved this function into a file named connect.hpp. Updated the related doxygen comments. --- bin/filter_mrc/filter_mrc.cpp | 6 +- bin/filter_mrc/handlers.cpp | 121 +-- bin/filter_mrc/handlers.hpp | 10 +- bin/filter_mrc/settings.cpp | 2 +- bin/filter_mrc/settings.hpp | 2 +- lib/visfd/clustering.hpp | 1 - lib/visfd/connect.hpp | 1401 +++++++++++++++++++++++++++++++++ lib/visfd/segmentation.hpp | 1346 ------------------------------- lib/visfd/visfd.hpp | 1 + 9 files changed, 1473 insertions(+), 1417 deletions(-) delete mode 100644 lib/visfd/clustering.hpp create mode 100644 lib/visfd/connect.hpp diff --git a/bin/filter_mrc/filter_mrc.cpp b/bin/filter_mrc/filter_mrc.cpp index 5b39270..e36a492 100644 --- a/bin/filter_mrc/filter_mrc.cpp +++ b/bin/filter_mrc/filter_mrc.cpp @@ -30,7 +30,7 @@ using namespace std; string g_program_name("filter_mrc"); string g_version_string("0.29.10"); -string g_date_string("2021-7-25"); +string g_date_string("2021-7-26"); @@ -540,11 +540,11 @@ int main(int argc, char **argv) { - else if (settings.filter_type == Settings::WATERSHED_DIRECTIONAL) { + else if (settings.filter_type == Settings::LABEL_CONNECTED) { // Cluster adjacent nearby voxels with high "saliency" into "islands" // neighboring voxels (this is similar to watershed segmentation). - HandleWatershedDirectional(settings, tomo_in, tomo_out, mask, voxel_width); + HandleLabelConnected(settings, tomo_in, tomo_out, mask, voxel_width); } diff --git a/bin/filter_mrc/handlers.cpp b/bin/filter_mrc/handlers.cpp index e631832..772bae0 100644 --- a/bin/filter_mrc/handlers.cpp +++ b/bin/filter_mrc/handlers.cpp @@ -17,13 +17,14 @@ using namespace std; #include using namespace visfd; -#include #include #include #include #include #include +#include #include +#include #include "err.hpp" #include "settings.hpp" #include "file_io.hpp" @@ -1296,11 +1297,11 @@ HandleWatershed(const Settings &settings, void -HandleWatershedDirectional(const Settings &settings, - MrcSimple &tomo_in, - MrcSimple &tomo_out, - MrcSimple &mask, - float voxel_width[3]) +HandleLabelConnected(const Settings &settings, + MrcSimple &tomo_in, + MrcSimple &tomo_out, + MrcSimple &mask, + float voxel_width[3]) { const vector > > *pMustLinkConstraints = nullptr; @@ -1330,32 +1331,32 @@ HandleWatershedDirectional(const Settings &settings, undefined_voxel_brightness = -1; size_t num_clusters = - WatershedDirectional(tomo_in.header.nvoxels, //image size - tomo_in.aaafI, //<-saliency - aaaiClusterId, //<-which cluster does each voxel belong to? (results will be stored here) - mask.aaafI, - settings.connect_threshold_saliency, - static_cast ***>(nullptr), - -std::numeric_limits::infinity(), - -std::numeric_limits::infinity(), - false, //normal default value for this (ignored) parameter - static_cast(nullptr), - -std::numeric_limits::infinity(), - -std::numeric_limits::infinity(), - true, //normal default value for this (ignored) parameter - 1, - static_cast(-1), // an impossible value - &cluster_centers, - &cluster_sizes, - &cluster_saliencies, - RegionSortCriteria::SORT_BY_SIZE, - static_cast(nullptr), - #ifndef DISABLE_STANDARDIZE_VECTOR_DIRECTION - static_cast ***>(nullptr), - #endif - pMustLinkConstraints, - true, //(clusters begin at regions of high saliency) - &cerr); //!< print progress to the user + LabelConnected(tomo_in.header.nvoxels, //image size + tomo_in.aaafI, //<-saliency + aaaiClusterId, //<-which cluster does each voxel belong to? (results will be stored here) + mask.aaafI, + settings.connect_threshold_saliency, + static_cast ***>(nullptr), + -std::numeric_limits::infinity(), + -std::numeric_limits::infinity(), + false, //normal default value for this (ignored) parameter + static_cast(nullptr), + -std::numeric_limits::infinity(), + -std::numeric_limits::infinity(), + true, //normal default value for this (ignored) parameter + 1, + static_cast(-1), // an impossible value + &cluster_centers, + &cluster_sizes, + &cluster_saliencies, + RegionSortCriteria::SORT_BY_SIZE, + static_cast(nullptr), + #ifndef DISABLE_STANDARDIZE_VECTOR_DIRECTION + static_cast ***>(nullptr), + #endif + pMustLinkConstraints, + true, //(clusters begin at regions of high saliency) + &cerr); //!< print progress to the user ptrdiff_t max_label = aaaiClusterId[0][0][0]; for (int iz = 0; iz < image_size[2]; ++iz) @@ -1386,7 +1387,7 @@ HandleWatershedDirectional(const Settings &settings, Dealloc3D(aaaiClusterId); -} //HandleWatershedDirectional() +} //HandleLabelConnected() @@ -1858,32 +1859,32 @@ HandleTV(const Settings &settings, size_t num_clusters = - WatershedDirectional(image_size, //image size - aaafSaliency, - aaaiClusterId, //<-which cluster does each voxel belong to? (results will be stored here) - mask.aaafI, - settings.connect_threshold_saliency, - aaaafDirection, - settings.connect_threshold_vector_saliency, - settings.connect_threshold_vector_neighbor, - false, //eigenvector signs are arbitrary so ignore them - aaaafVoteTensor, - settings.connect_threshold_tensor_saliency, - settings.connect_threshold_tensor_neighbor, - true, //the tensor should be positive definite near the target - 1, - static_cast(-1), // an impossible value - &cluster_centers, - &cluster_sizes, - &cluster_saliencies, - RegionSortCriteria::SORT_BY_SIZE, - static_cast(nullptr), - #ifndef DISABLE_STANDARDIZE_VECTOR_DIRECTION - aaaafDirection, - #endif - pMustLinkConstraints, - true, //(clusters begin at regions of high saliency) - &cerr); //!< print progress to the user + LabelConnected(image_size, //image size + aaafSaliency, + aaaiClusterId, //<-which cluster does each voxel belong to? (results will be stored here) + mask.aaafI, + settings.connect_threshold_saliency, + aaaafDirection, + settings.connect_threshold_vector_saliency, + settings.connect_threshold_vector_neighbor, + false, //eigenvector signs are arbitrary so ignore them + aaaafVoteTensor, + settings.connect_threshold_tensor_saliency, + settings.connect_threshold_tensor_neighbor, + true, //the tensor should be positive definite near the target + 1, + static_cast(-1), // an impossible value + &cluster_centers, + &cluster_sizes, + &cluster_saliencies, + RegionSortCriteria::SORT_BY_SIZE, + static_cast(nullptr), + #ifndef DISABLE_STANDARDIZE_VECTOR_DIRECTION + aaaafDirection, + #endif + pMustLinkConstraints, + true, //(clusters begin at regions of high saliency) + &cerr); //!< print progress to the user ptrdiff_t max_label = aaaiClusterId[0][0][0]; for (int iz = 0; iz < image_size[2]; ++iz) @@ -1926,7 +1927,7 @@ HandleTV(const Settings &settings, // where is the lookup table to indicate which cluster a voxel belongs to? float ***aaafVoxel2Cluster = nullptr; if (settings.cluster_connected_voxels) - aaafVoxel2Cluster = tomo_out.aaafI; //tomo_out was filled by WatershedDirectional() + aaafVoxel2Cluster = tomo_out.aaafI; //tomo_out was filled by LabelConnected() for (int iz = 0; iz < image_size[2]; ++iz) { for (int iy = 0; iy < image_size[1]; ++iy) { diff --git a/bin/filter_mrc/handlers.hpp b/bin/filter_mrc/handlers.hpp index ef0b4a8..ba607b8 100644 --- a/bin/filter_mrc/handlers.hpp +++ b/bin/filter_mrc/handlers.hpp @@ -168,11 +168,11 @@ HandleWatershed(const Settings &settings, void -HandleWatershedDirectional(const Settings &settings, - MrcSimple &tomo_in, - MrcSimple &tomo_out, - MrcSimple &mask, - float voxel_width[3]); +HandleLabelConnected(const Settings &settings, + MrcSimple &tomo_in, + MrcSimple &tomo_out, + MrcSimple &mask, + float voxel_width[3]); void diff --git a/bin/filter_mrc/settings.cpp b/bin/filter_mrc/settings.cpp index 12a57ee..a2437bd 100644 --- a/bin/filter_mrc/settings.cpp +++ b/bin/filter_mrc/settings.cpp @@ -3370,7 +3370,7 @@ Settings::ParseArgs(vector& vArgs) (filter_type != CURVE)) { assert(connect_threshold_saliency!=std::numeric_limits::infinity()); - filter_type = WATERSHED_DIRECTIONAL; + filter_type = LABEL_CONNECTED; } diff --git a/bin/filter_mrc/settings.hpp b/bin/filter_mrc/settings.hpp index d846b95..d122a6e 100644 --- a/bin/filter_mrc/settings.hpp +++ b/bin/filter_mrc/settings.hpp @@ -50,7 +50,7 @@ class Settings { LOG_DOG, // Approximation to Laplacian-of-Guassian using DOG (fast) LOCAL_FLUCTUATIONS, // Report the fluctuation of nearby voxel intensities WATERSHED, // Watershed segmentation - WATERSHED_DIRECTIONAL, // Watershed applied to voxels with tensor attributes + LABEL_CONNECTED, // Connected component labeling (count islands in a sea) BLOB, // Scale-free blob detection CURVE, // Detect 1D curve-like ridges (a filament detector) SURFACE_EDGE, // Detect the edge of a light-dark boundary (gradient) diff --git a/lib/visfd/clustering.hpp b/lib/visfd/clustering.hpp deleted file mode 100644 index 8b13789..0000000 --- a/lib/visfd/clustering.hpp +++ /dev/null @@ -1 +0,0 @@ - diff --git a/lib/visfd/connect.hpp b/lib/visfd/connect.hpp new file mode 100644 index 0000000..dcbcf86 --- /dev/null +++ b/lib/visfd/connect.hpp @@ -0,0 +1,1401 @@ +/// @file connect.hpp +/// @brief A collection of functions for connected component labeling. +/// (https://en.wikipedia.org/wiki/Connected-component_labeling) +/// @author Andrew Jewett +/// @date 2021-7-27 + +#ifndef _CONNECT_HPP +#define _CONNECT_HPP + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +using namespace std; +#include // defines the "VisfdErr" exception type +#include // defines matrix diagonalizer (DiagonalizeSym3()) +#include // defines "Normalize3()" +#include // defines invert_permutation(), AveArray(), ... +#include // defines Alloc2D() and Dealloc2D() +#include // defines Alloc3D() and Dealloc3D() +#include // defines "Filter1D" (used in ApplySeparable()) +#include // defines common 3D image filters +#include // defines FindExtrema() + + +namespace visfd { + + +/// @brief The following enumerated constants are used as arguments to the +/// LabelConnected() function. This function +/// divides an image into multiple regions. Each region is associated +/// with a unique integer (a "label", which is similar to an ID number) +/// We need a way to specify how these labels are assigned to +/// these regions. +/// +/// SORT_BY_SIZE will sort regions by their size (number of voxels). +/// +/// SORT_BY_VALUE will sort regions by brightness values inside +/// the region (typically by their darkest or brightest voxel). + +typedef enum eRegionSortCriteria { + SORT_BY_VALUE, + SORT_BY_SIZE +} RegionSortCriteria; + + + +/// @brief This function is used for connected-component labeling. +/// https://en.wikipedia.org/wiki/Connected-component_labeling +/// It works on binary and grayscale images, as well as images with +/// directional (and tensor) voxel features. +/// +/// This function is used to identify different "islands" of bright +/// voxels separated by a sea of other voxels of low brightness +/// (...or incompatible directions, as explained below). +/// These different islands may correspond to different objects +/// in the source image. However this function can also consider +/// voxel's vector and tensor attributes (if present). This function +/// can decline to merge adjacent voxels into the same "island" +/// unless their directions (and tensors) are sufficiently similar. +/// +/// @warning EXPERIMENTAL CODE. THIS FUNCTION IS COMPLICATED AND THIS MAKES IT +/// DIFFICULT TO USE. THIS FUNCTION'S ARGUMENT LIST MAY CHANGE +/// SIGNIFICANTLY IN THE FUTURE. +/// +/// @details Based on the watershed algorithm, this function first +/// finds the locations of all of the local maxima in brightness (a.k.a. +/// "saliency") in the image. These form the peaks of the islands. +/// Then, starting from these locations of high brightness, it searches +/// outwards for adjacent (neighboring) voxels which are also bright and +/// adds them to the "island". (I sometimes call islands "clusters" or +/// "watershed basins".) The islands continue to grow until they reach +/// voxel brightnesses (or saliencies) which fall below some threshold +/// specified by the caller. This is the "shoreline" of the island. +/// If during this outward-search two islands collide with each other +/// they may be merged into one island. However there is one +/// exception: If the voxels in the image are imbued with directional +/// attributes (ie. if "aaaafDirection" or "aaaafSymmetricTensor" +/// are non-NULL), then the directions of these neighboring voxels +/// must be similar. If the neighboring voxels from either island +/// are not pointing in sufficiently similar directions (as specified +/// by directional thresholds supplied by the caller), then the islands +/// will not be merged into one island, even if both voxels are bright. +/// +/// In other words, this function performs the agglomerative clustering +/// of voxels, considering voxel direction, proximity, and brightness. +/// +/// @return The function returns the number of clusters ("islands") found. +/// After the function is finished, the aaaiDest[][][] array will store +/// the cluster-ID (a.k.a. "island-ID") associated with each voxel. +/// The clusters are sorted by their size (by default), and this is +/// relfected in their cluster-ID numbers ("labels"). (The largest +/// cluster has an ID of 1. Cluster-ID numbers begin at 1, not 0.) +/// (The "pv_cluster_maxima" argument, if != nullptr, will store +/// the location of the saliency maxima (or minima) for each cluster.) +/// +/// @note Voxels below the saliency threshold are ignored, and will +/// assigned to the value of UNDEFINED, which (by default) is -1. +/// (This value can be overridden using the "label_undefined" argument.) +/// +/// @note If a aaafMask array is supplied by the caller, then voxels located +/// in regions where aaafMask[iz][iy][ix]=0 will be ignored. +/// +/// The remaining notes below describe the behavior of the optional +/// aaaafVector, aaaafSymmetricTensor, aaaafVectorStandardized arguments, +/// which are not generally useful in all situations +/// @note If aaaafVector != nullptr, then +/// neighboring voxels whose direction (vector) changes discontinuously, +/// will not be added to an existing cluster. +/// @note If aaaafSymmetricTensor != nullptr, then +/// neighboring voxels whose direction (tensor) changes discontinuously, +/// will not be added to an existing cluster. +/// @note In addition, voxels whose Hessian (2nd derivative matrix) of the +/// aaafSaliency array is INCOMPATIBLE with the corresponding entry from +/// aaaafSymmetricTensor (if!=nullptr) will not be added to any cluster. +/// (Tensors are compared using the normalized Trace-product.) +/// @note In addition, voxels whose principle eigenvector from the Hessian +/// (2nd derivative matrix) of the aaafSaliency array is POINTING IN A +/// SIGNIFICANTLY DIFFERENT DIRECTION FROM the corresponding entry from +/// the aaaafVector array (if != nullptr) will not be added to any cluster. +/// @note If (consider_dot_product_sign == false), the difference in vector +/// directions will be calculated by considering only the magnitude +/// of the dot product between them (not the sign). +/// (This is useful because some image processing algorithms generate +/// vectors which lack polarity. For example, planar surface "ridge" +/// detectors calculate vectors normal to a surface, which are equally +/// likely to point inward or outward.) +/// @note If (consider_dot_product_sign == false) AND +/// if aaaafVector and aaaafVectorStandardized are both non-null, THEN +/// aaaafVectorStandardized array will store a version of aaaafVector +/// array whose signs have been flipped to preserve consistency +/// of directionality between adjacent voxels as much as possible. +/// Why? This is necessary because some methods for generating voxel +/// directions are unable to determine the sign of each voxel's +/// direction. (The eigenvectors of the Hessian are a good example.) +/// By supplying an "aaaafVectorStandardized" array, you can ask this +/// function to choose the sign of these vectors. This will ensure +/// that all of the vectors in the same island are pointing +/// in directions that are locally compatible with each other. +/// Robustness: If any clusters contain any any closed loops +/// that have directions which cannot be chosen consistently, +/// they will be cut at a location where their saliency is weakest. +/// (This will avoid Möbius-strip-like defects with an odd number of +/// twists. Since most noisy images have defects with most one twist, +/// if any, this will eliminate the vast majority of random defects.) + +template + +size_t +LabelConnected(const int image_size[3], //!< #voxels in xyz + Scalar const *const *const *aaafSaliency, //!< brightness of each voxel (a.k.a. "saliency") + Label ***aaaiDest, //!< watershed segmentation results go here + Scalar const *const *const *aaafMask, //!< optional: Ignore voxels whose mask value is 0 + Scalar threshold_saliency=-std::numeric_limits::infinity(), //!< don't consider voxels with saliencies below this value + VectorContainer const *const *const *aaaafVector=nullptr, + Scalar threshold_vector_saliency=-std::numeric_limits::infinity(), //!< voxels with incompatible saliency and vector are ignored (numbers below -1 disable this feature) + Scalar threshold_vector_neighbor=-std::numeric_limits::infinity(), //!< neighboring voxels with incompatible vectors are ignored (numbers below -1 disable this feature) + bool consider_dot_product_sign = true, //!< does the sign of the dot product matter? If not, compare abs(DotProduct()) with threshold_vector variables + TensorContainer const *const *const *aaaafSymmetricTensor=nullptr, + Scalar threshold_tensor_saliency=-std::numeric_limits::infinity(), //!< voxels with incompatible saliency and tensor are ignored (numbers below -1 disable this feature) + Scalar threshold_tensor_neighbor=-std::numeric_limits::infinity(), //!< neighboring voxels with incompatible tensors are ignored (numbers below -1 disable this feature) + bool tensor_is_positive_definite_near_target=true, //!< what is the sign of the principal tensor eigenvalue(s) near a ridge we care about? + int connectivity=1, //!< OPTIONAL: square root of the search radius around each voxel (1=nearest_neighbors, 2=2D_diagonal, 3=3D_diagonal) + Label label_undefined=-1, //!< OPTIONAL: voxels storing this value do not belong to any clusters + vector > *pv_cluster_maxima=nullptr, //!< OPTIONAL: the location of saliency minima or maxima which seeded each cluster + vector *pv_cluster_sizes=nullptr, //!< OPTIONAL: what was the size of each cluster? (either the number of voxels, or the sum of voxel weights) + vector *pv_cluster_saliencies=nullptr, //!< OPTIONAL: what was the saliency (brightness) of each cluster's brightest voxel? + RegionSortCriteria sort_criteria = RegionSortCriteria::SORT_BY_SIZE, //!< which clusters get reported first? (by default, the biggest ones) + Scalar const *const *const *aaafVoxelWeights=nullptr, //!< OPTIONAL: weights of each voxel used when sort_criteria==SORT_BY_SIZE + #ifndef DISABLE_STANDARDIZE_VECTOR_DIRECTION + VectorContainer ***aaaafVectorStandardized=nullptr, //!< OPTIONAL: place to store "standardized" vector directions (with flipped sign to agree with neighboring voxels). + #endif + const vector > > *pMustLinkConstraints=nullptr, //!< OPTIONAL: a list of sets of voxel locations. This insures that voxels in each set will belong to the same cluster. + bool start_from_saliency_maxima=true, //!< start from local maxima? (if false, minima will be used) WARNING: As of 2021-7-26, this function has not yet been tested with the non-default value (false) + ostream *pReportProgress=nullptr) //!< print progress to the user? +{ + selfadjoint_eigen3::EigenOrderType eival_order; + if (start_from_saliency_maxima) + eival_order = selfadjoint_eigen3::DECREASING_EIVALS; //<--first eigenvalue will be the largest eigenvalue + else + eival_order = selfadjoint_eigen3::INCREASING_EIVALS; //<--first eigenvalue will be the smallest (most negative) eigenvalue + + assert(image_size); + assert(aaafSaliency); + assert(aaaiDest); + + if (! consider_dot_product_sign) { + // Weird detail: + // A negative threshold_vector should mean that the caller wants us to + // ignore the vector argument and not attempt to compare vectors + // associated with each voxel. + // (At least that's how I invoke it. -AJ 2019-2-11) + // However if (consider_dot_product_sign==false), then the absolute + // value of the threshold is used, which will be a positive number. + // This will cause voxels whose dot product is less than the absolute + // value of this threshold to be ignored. + // This is not what we want. + // In that case we want to choose a threshold which is as low as possible + // for the absolute value of a dot product, which is 0.0; + if (threshold_vector_saliency < 0) + threshold_vector_saliency = 0.0; + if (threshold_vector_neighbor < 0) + threshold_vector_neighbor = 0.0; + } + + + // 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 + + + Scalar SIGN_FACTOR = 1.0; + if (start_from_saliency_maxima) { + SIGN_FACTOR = -1.0; + } + + vector > extrema_locations; //where is the minima/maxima? + vector extrema_scores; //how bright is this minima or maxima? + vector extrema_nvoxels; //(needed to pacify syntax of FindExtrema3d()) + + // Find all the local minima (or maxima?) in the image. + + _FindExtrema(image_size, + aaafSaliency, + aaafMask, + extrema_locations, + extrema_scores, + extrema_nvoxels, + (! start_from_saliency_maxima), //<-- minima or maxima? + threshold_saliency, + connectivity, + true, // maxima are allowed to be located on the image border + static_cast(nullptr), + pReportProgress); + + ptrdiff_t UNDEFINED = extrema_locations.size() + 1; //an impossible value + ptrdiff_t QUEUED = extrema_locations.size() + 2; //an impossible value + + //initialize aaaiDest[][][] + for (int iz=0; iz // location of the voxel + > + > q; + + + if (pReportProgress) + *pReportProgress << + "-- Clustering voxels belonging to different objects --\n" + "starting from " << extrema_locations.size() << " different local " + << (start_from_saliency_maxima ? "maxima" : "minima") + << endl; + + // Initialize the queue with the voxels at these minima locations + + for (size_t i=0; i < extrema_locations.size(); i++) { + // Create an entry in q for each of the local minima + + // Assign a different integer to each of these minima, starting at 1 + ptrdiff_t which_basin = i; + + int ix = extrema_locations[i][0]; + int iy = extrema_locations[i][1]; + int iz = extrema_locations[i][2]; + + // These entries in the priority queue will be sorted by "score" + // which is the intensity of the image at this voxel location. + Scalar score = extrema_scores[i]; + assert(score == aaafSaliency[iz][iy][ix]); + score *= SIGN_FACTOR; //(enable search for either local minima OR maxima) + + // Note:FindExtrema() should avoid minima above threshold_saliency, + // or maxima below threshold_saliency. We check for that with an assert: + assert(score <= threshold_saliency * SIGN_FACTOR); + + // copy the ix,iy,iz coordinates into an array + array icrds; + // (It seems like there should be a way to do this in 1 line, but I'm + // unfamiliar with the new C++ array initializer syntax and I'm on a plane + // without internet at the moment. So I'll just do it the obvious way.) + icrds[0] = ix; + icrds[1] = iy; + icrds[2] = iz; + + q.push(make_tuple(-score, // <-- entries sorted lexicographically by -score + which_basin, + icrds)); + + + assert(aaaiDest[iz][iy][ix] == UNDEFINED); + aaaiDest[iz][iy][ix] = QUEUED; + + } // for (size_t i=0; i < extrema_locations.size(); i++) + + + // Each cluster contains one or more "basins" + // A "basin" is a group of connected voxels which are all part of the + // same local minima (or maxima). + // Initially, each basin is it's own cluster. + // (Later we may merge multiple basins into a single cluster.) + vector basin2cluster(extrema_locations.size()); + for (size_t i=0; i < basin2cluster.size(); i++) + basin2cluster[i] = i; + + // Inverse lookup table + vector > cluster2basins(extrema_locations.size()); + for (size_t i=0; i < basin2cluster.size(); i++) { + cluster2basins[i] = set(); + cluster2basins[i].insert(i); + } + + + // Count the number of voxels in the image we will need to consider. + // (This will be used for printing out progress estimates later.) + size_t n_voxels_processed = 0; + //size_t n_voxels_image = image_size[0] * image_size[1] * image_size[2]; + size_t n_voxels_image = 0; + for (int iz=0; iz + SIGN_FACTOR * threshold_saliency) + continue; + n_voxels_image++; + } + } + } + + #ifndef DISABLE_STANDARDIZE_VECTOR_DIRECTION + // initialize aaaafVectorStandardized[][][] + if (aaaafVector && aaaafVectorStandardized) { + for (int iz=0; iz basin2polarity(extrema_locations.size(), 1); + bool voxels_discarded_due_to_polarity = false; + Scalar voxel_discarded_due_to_polarity_saliency = -1.0; + int voxel_discarded_due_to_polarity_ix = -1; // an impossible value + int voxel_discarded_due_to_polarity_iy = -1; // an impossible value + int voxel_discarded_due_to_polarity_iz = -1; // an impossible value + #endif // #ifndef DISABLE_STANDARDIZE_VECTOR_DIRECTION + + // Loop over all the voxels on the periphery + // of the voxels with lower intensity (brightness). + // This is a priority-queue of voxels which lie adjacent to + // the voxels which have been already considered. + // As we consider new voxels, add their neighbors to this queue as well + // until we process all voxels in the image. + + //#ifndef NDEBUG + //set > visited; + //#endif + + while (! q.empty()) + { + tuple > p = q.top(); + q.pop(); + + // Figure out the properties of that voxel (position, intensity/score,basin) + Scalar i_score = -std::get<0>(p); //abs(i_score) = voxel intensity = aaafSaliency[iz][iy][ix] + Scalar i_which_basin = std::get<1>(p); // the basin to which that voxel belongs (tentatively) + int ix = std::get<2>(p)[0]; // voxel location + int iy = std::get<2>(p)[1]; // " " + int iz = std::get<2>(p)[2]; // " " + + // Should we ignore this voxel? + + if (i_score > threshold_saliency * SIGN_FACTOR) { + // stop when the voxel brightness(*SIGN_FACTOR) exceeds the threshold_saliency + aaaiDest[iz][iy][ix] = UNDEFINED; + continue; + } + + if (aaafMask && aaafMask[iz][iy][ix] == 0.0) { + // ignore voxel if the user specified a mask and the voxel does not belong + aaaiDest[iz][iy][ix] = UNDEFINED; + continue; + } + + + { // Use inconsistencies in aaafSaliency to discard voxel ix,iy,iz? + + // ----------------------------------------------- + // First, condsider discarding voxel ix,iy,iz due to + // inconsistencies between aaafSaliency and aaaafSymTensor here + // ----------------------------------------------- + + Scalar saliency_hessian3x3[3][3]; + + CalcHessianFiniteDifferences(aaafSaliency, //!< source image + ix, iy, iz, + saliency_hessian3x3, + image_size); + + // Confusing sign compatibility issue: + // We assume that the tensor passed + // by the caller will be positive-definite near a ridge. + // In other words, it's largest eigenvalue will be positive there, + // and its corresponding eigenvector is assumed to point in a direction + // where the 2nd-derivative is largest. + // In that case we must invert the sign of the hessian matrix before we + // compare it with the tensor, because we also assume the saliency is + // bright on a dark background. Hence it's second derivative (ie. its + // Hessian) along that direction will be negative instead of positive + // at that location. Hence if both of these assumptions are true, + // we must multiply the entries in saliency_hessian by -1 before + // comparing them with the tensor. + + if (tensor_is_positive_definite_near_target == + start_from_saliency_maxima) { + for (int di = 0; di < 3; di++) + for (int dj = di; dj < 3; dj++) + saliency_hessian3x3[di][dj] *= -1.0; + } + + // To reduce memory consumption, + // save the resulting 3x3 matrix in a smaller 1-D array (with 6 entries) + // whose index is given by MapIndices_3x3_to_linear[][] + + Scalar saliency_hessian[6]; + + for (int di = 0; di < 3; di++) + for (int dj = di; dj < 3; dj++) + saliency_hessian[ MapIndices_3x3_to_linear[di][dj] ] + = saliency_hessian3x3[di][dj]; + + bool discard_this_voxel = false; + + if (aaaafSymmetricTensor) { + + Scalar tp = TraceProductSym3(saliency_hessian, + aaaafSymmetricTensor[iz][iy][ix]); + Scalar fs = FrobeniusNormSym3(saliency_hessian); + Scalar ft = FrobeniusNormSym3(aaaafSymmetricTensor[iz][iy][ix]); + + if (tp < threshold_tensor_saliency * fs * ft) { + discard_this_voxel = true; + } + } + + // ----------------------------------------------- + // Then, condsider discarding voxel ix,iy,iz due to + // inconsistencies between aaafSaliency and aaaafVector at this location + // ----------------------------------------------- + + Scalar s_eivals[3]; + Scalar s_eivecs[3][3]; + + // the eigenvector of the saliency_hessian that we care about is the + // one with the largest eigenvalue (which is assumed to be positive). + + ConvertFlatSym2Evects3(saliency_hessian, + s_eivals, + s_eivecs, + eival_order); + + if (aaaafVector) { + bool vect_threshold_exceeded = false; + + if (consider_dot_product_sign) { + if (DotProduct3(s_eivecs[0], //principal (first) eigenvector + aaaafVector[iz][iy][ix]) + < + (threshold_vector_saliency * + Length3(s_eivecs[0]) * + Length3(aaaafVector[iz][iy][ix]))) + vect_threshold_exceeded = true; + } + else { + if (SQR(DotProduct3(s_eivecs[0], //principal (first) eigenvector + aaaafVector[iz][iy][ix])) + < + (SQR(threshold_vector_saliency) * + SquaredNorm3(s_eivecs[0]) * + SquaredNorm3(aaaafVector[iz][iy][ix]))) + vect_threshold_exceeded = true; + } + if (vect_threshold_exceeded) { + discard_this_voxel = true; + } + } + + if (discard_this_voxel) { + aaaiDest[iz][iy][ix] = UNDEFINED; + // Are we deleting a voxel which is also a basin local minima/maxima? + if ((ix == extrema_locations[i_which_basin][0]) && + (iy == extrema_locations[i_which_basin][1]) && + (iz == extrema_locations[i_which_basin][2])) + // If so, then this entire basin should be deleted from consieration + // Next line effectively deletes i_which_basin from basin2cluster[] + basin2cluster[i_which_basin] = -1; + continue; + } + + } // Use inconsistencies in aaafSaliency to discard voxel ix,iy,iz? + + + assert(aaaiDest[iz][iy][ix] == QUEUED); + // Now we assign this voxel to the basin + aaaiDest[iz][iy][ix] = i_which_basin; + // (Note: This will prevent the voxel from being visited again.) + + + if (pReportProgress) { + // Every time the amount of progress increases by 1%, inform the user: + n_voxels_processed++; + size_t percentage = (n_voxels_processed*100) / n_voxels_image; + size_t percentage_previous = ((n_voxels_processed-1)*100)/n_voxels_image; + if (percentage != percentage_previous) + *pReportProgress << " percent complete: " << percentage << endl; + } + + + //#ifndef NDEBUG + //array ixiyiz; + //ixiyiz[0] = ix; + //ixiyiz[1] = iy; + //ixiyiz[2] = iz; + //assert(visited.find(ixiyiz) == visited.end()); + //visited.insert(ixiyiz); + //#endif //#ifndef NDEBUG + + + // ---------- Meyer (and Beucher's?) inter-pixel flood algorith: --------- + // Check the voxels that surround voxel (ix,iy,iz) + + // ---- loop over neighbors ---- + // Let (jx,jy,jz) denote the index offset for the neighbor voxels + // (located at ix+jx,iy+jy,iz+jz) + 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)) + continue; + if ((iy_jy < 0) || (image_size[1] <= iy_jy)) + continue; + if ((ix_jx < 0) || (image_size[0] <= ix_jx)) + continue; + + if (aaafMask && (aaafMask[iz_jz][iy_jy][ix_jx] == 0.0)) + continue; + + + { // Difference between voxel ix,iy,iz and ix_jx,iy_jy,iz_jz too large? + + // ----------------------------------------------- + // Then, condsider discarding voxel ix_jx,iy_jy,iz_jz due to + // inconsistencies between aaaafSymTensor[iz][iy][ix] + // and aaaafSymTensor[iz_jz][iy_jy][ix_jx] + // ----------------------------------------------- + if (aaaafSymmetricTensor) { + if (TraceProductSym3(aaaafSymmetricTensor[iz][iy][ix], + aaaafSymmetricTensor[iz_jz][iy_jy][ix_jx]) + < + (threshold_tensor_neighbor * + FrobeniusNormSym3(aaaafSymmetricTensor[iz][iy][ix]) * + FrobeniusNormSym3(aaaafSymmetricTensor[iz_jz][iy_jy][ix_jx]))) + { + continue; + } + } + + // ----------------------------------------------- + // Then, condsider discarding voxel ix_jx,iy_jy,iz_jz due to + // inconsistencies between aaaafVector[iz][iy][ix] + // and aaaafVector[iz_jz][iy_jy][ix_jx] + // ----------------------------------------------- + if (aaaafSymmetricTensor) { + if (consider_dot_product_sign) { + if (DotProduct3(aaaafVector[iz][iy][ix], + aaaafVector[iz_jz][iy_jy][ix_jx]) + < + (threshold_tensor_neighbor * + Length3(aaaafVector[iz][iy][ix])* + Length3(aaaafVector[iz_jz][iy_jy][ix_jx]))) + { + continue; + } + } + else { + if (SQR(DotProduct3(aaaafVector[iz][iy][ix], + aaaafVector[iz_jz][iy_jy][ix_jx])) + < + (SQR(threshold_vector_neighbor) * + SquaredNorm3(aaaafVector[iz][iy][ix])* + SquaredNorm3(aaaafVector[iz_jz][iy_jy][ix_jx]))) + { + continue; + } + } + } + } // Difference between voxel ix,iy,iz and ix_jx,iy_jy,iz_jz too large? + + + if (aaaiDest[iz_jz][iy_jy][ix_jx] == QUEUED) { + continue; + } + else if (aaaiDest[iz_jz][iy_jy][ix_jx] == UNDEFINED) + { + + aaaiDest[iz_jz][iy_jy][ix_jx] = QUEUED; + + // and push this neighboring voxels onto the queue. + // (...if they have not been assigned to a basin yet. This + // insures that the same voxel is never pushed more than once.) + array neighbor_crds; + neighbor_crds[0] = ix_jx; + neighbor_crds[1] = iy_jy; + neighbor_crds[2] = iz_jz; + Scalar neighbor_score = aaafSaliency[iz_jz][iy_jy][ix_jx]*SIGN_FACTOR; + q.push(make_tuple(-neighbor_score, + i_which_basin, + neighbor_crds)); + + + #ifndef DISABLE_STANDARDIZE_VECTOR_DIRECTION + if (aaaafVector && aaaafVectorStandardized && + (! consider_dot_product_sign)) + { + if (DotProduct3(aaaafVectorStandardized[iz][iy][ix], + aaaafVectorStandardized[iz_jz][iy_jy][ix_jx]) + < + 0.0) + { + // Voxels added as the basin is growing should always have + // directions which are compatible with each other. + // (IE the dot product between them should not be negative.) + // (Why? See "Discussion" below) + // If this is not the case, then invert the sign of the newcommers. + // (Note: The voxels within different basins are made compatible + // by flipping the sign of their entry in "basin2polarity[]") + aaaafVectorStandardized[iz_jz][iy_jy][ix_jx][0] *= -1.0; + aaaafVectorStandardized[iz_jz][iy_jy][ix_jx][1] *= -1.0; + aaaafVectorStandardized[iz_jz][iy_jy][ix_jx][2] *= -1.0; + // Discussion: + // It is possible to assume this because the topology of the + // basins that we are growing at this point cannot have loops. + // This is because each neighboring voxel we add in this step + // is of type UNDEFINED. Later we may encounter loops when we + // run into neighboring voxels which have already beem processed. + // We will worry about Möbius loops then. + } + } // if (aaaafVectorStandardized && (! consider_dot_product_sign)) + #endif + + + } // else if (aaaiDest[iz_jz][iy_jy][ix_jx] == UNDEFINED) + else { + + assert(aaaiDest[iz][iy][ix] == i_which_basin); + + ptrdiff_t basin_i = aaaiDest[iz][iy][ix]; + ptrdiff_t basin_j = aaaiDest[iz_jz][iy_jy][ix_jx]; + + ptrdiff_t cluster_i = basin2cluster[basin_i]; + ptrdiff_t cluster_j = basin2cluster[basin_j]; + + #ifndef DISABLE_STANDARDIZE_VECTOR_DIRECTION + bool polarity_match = true; + if (aaaafVectorStandardized && (! consider_dot_product_sign)) + { + if (DotProduct3(aaaafVectorStandardized[iz][iy][ix], + aaaafVectorStandardized[iz_jz][iy_jy][ix_jx]) + * + basin2polarity[basin_i] + * + basin2polarity[basin_j] + < + 0.0) + polarity_match = false; + } + #endif + + if (cluster_i == cluster_j) { + #ifndef DISABLE_STANDARDIZE_VECTOR_DIRECTION + if (! polarity_match) { + voxels_discarded_due_to_polarity = true; + voxel_discarded_due_to_polarity_ix = ix_jx; + voxel_discarded_due_to_polarity_iy = iy_jy; + voxel_discarded_due_to_polarity_iz = iz_jz; + voxel_discarded_due_to_polarity_saliency = aaafSaliency[iz_jz][iy_jy][ix_jx]; + + // In that case throw away this voxel. + // Hopefully this will prevent linking of voxels containing + // vector directors which cannot be reconciled. Example: + // This should prevent surfaces which resemble a Möbius strip + // from forming a complete closed loop. Hopefully it will + // cut the loop at the voxels with the most tenuous connections. + continue; + } + #endif + } + else //if (cluster_i != cluster_j) + { + // -- merge the two clusters to which basin_i and basin_j belong -- + + // (arbitrarily) choose the cluster with the smaller ID number + // to absorb the other cluster, and delete the other cluster. + ptrdiff_t merged_cluster_id = std::min(cluster_i, cluster_j); + ptrdiff_t deleted_cluster_id = std::max(cluster_i, cluster_j); + assert(cluster_i != cluster_j); + + // copy the basins from the deleted cluster into the merged cluster + for (auto p = cluster2basins[deleted_cluster_id].begin(); + p != cluster2basins[deleted_cluster_id].end(); + p++) + { + ptrdiff_t basin_id = *p; + cluster2basins[merged_cluster_id].insert(basin_id); + basin2cluster[basin_id] = merged_cluster_id; + + #ifndef DISABLE_STANDARDIZE_VECTOR_DIRECTION + if (aaaafVector && aaaafVectorStandardized && + (! consider_dot_product_sign)) + { + if (! polarity_match) + basin2polarity[basin_id] *= -1.0; + } + #endif + } + + // -- delete all the basins from the deleted cluster -- + cluster2basins[deleted_cluster_id].clear(); + + } // if (aaaiDest[iz_jz][iy_jy][ix_jx] != aaaiDest[iz][iy][ix]) + } // else clause for "if (aaaiDest[iz_jz][iy_jy][ix_jx] == UNDEFINED)" + } // for (int j = 0; j < num_neighbors; j++) + } //while (q.size() != 0) + + + #ifndef NDEBUG + // DEBUGGING + // All of the voxels should either be assigned to a basin or to "UNDEFINED" + // (but not "QUEUED"). Check for that below: + for (int iz=0; izsize(); i_group++) + { + + // some variables we will need + + ptrdiff_t FIRST_ITER = -1; // an impossible value used below + + //r_i_init = coordinates of most recent voxel selected by the user + array r_i_init={-1,-1,-1}; + //r_i = coordinates of the voxel in aaaiDest[][][] nearest to r_i_init + array r_i; + //basin_i = the ID number for the basin to which voxel r_i belongs + ptrdiff_t basin_i = FIRST_ITER; + + //r_j_init, r_j, basin_j correspond to the PREVIOUSLY processed voxel + ptrdiff_t basin_j = FIRST_ITER; + array r_j_init={-1,-1,-1}; //an impossible value + array r_j={-1,-1,-1}; //an impossible value + + + // Loop over the voxels in each group + // and force them to belong to the same cluster + + for (auto pLocation = (*pMustLinkConstraints)[i_group].begin(); + pLocation != (*pMustLinkConstraints)[i_group].end(); + pLocation++) + { + r_i_init[0] = int(floor((*pLocation)[0]+0.5)); + r_i_init[1] = int(floor((*pLocation)[1]+0.5)); + r_i_init[2] = int(floor((*pLocation)[2]+0.5)); + if (*pReportProgress) { + *pReportProgress << " finding voxel nearest to ("; + for (int d = 0; d < 3; d++) { + *pReportProgress << r_i_init[d]; + if (d+1 != 3) *pReportProgress << ","; + } + *pReportProgress << ") -> ("; + } + + set