Merge pull request #151 from AdaptiveParticles/autoparameters

Add automatic parameter tuning based on minimum cross-entropy thresholding

Author: joeljonsson

src/algorithm/APRConverter.hpp

@@ -113,6 +113,9 @@ class APRConverter {
     template<typename T,typename S>
     void autoParameters(const PixelData<T> &localIntensityScale,const PixelData<S> &grad);
 
+    template<typename T,typename S>
+    void autoParametersLiEntropy(const PixelData<T> &image, const PixelData<T> &localIntensityScale, const PixelData<S> &grad);
+
     template<typename T>
     bool check_input_dimensions(PixelData<T> &input_image);
 
@@ -459,7 +462,11 @@ inline bool APRConverter<ImageType>::get_apr(APR &aAPR, PixelData<T>& input_imag
     computation_timer.start_timer("apply_parameters");
 
     if( par.auto_parameters ) {
-        autoParameters(local_scale_temp,grad_temp);
+        method_timer.start_timer("autoParameters");
+//        autoParameters(local_scale_temp,grad_temp);
+        autoParametersLiEntropy(local_scale_temp2, local_scale_temp, grad_temp);
+        aAPR.parameters = par;
+        method_timer.stop_timer();
     }
 
     applyParameters(aAPR,par);
@@ -587,92 +594,152 @@ inline bool APRConverter<ImageType>::get_apr(APR &aAPR, PixelData<T>& input_imag
     return true;
 }
 
-template<typename ImageType>
-template<typename T,typename S>
-void APRConverter<ImageType>::autoParameters(const PixelData<T> &localIntensityScale,const PixelData<S> &grad){
-    /*
-     *  Assumes a dark background. Please use the Python interactive parameter selection for more detailed approaches.
-     *
-     *  Finds the flatest (1% of the image) and estimates the noise level in the gradient and sets the grad threshold from that.
-     */
+template<typename T>
+void compute_means(const std::vector<T>& data, float threshold, float& mean_back, float& mean_fore) {
+    float sum_fore=0.f, sum_back=0.f;
+    size_t count_fore=0, count_back=0;
 
+#ifdef HAVE_OPENMP
+#pragma omp parallel for default(none) shared(data, threshold) reduction(+:sum_fore, sum_back, count_fore, count_back)
+#endif
+    for(size_t idx = 0; idx < data.size(); ++idx) {
+        if(data[idx] > threshold) {
+            sum_fore += data[idx];
+            count_fore++;
+        } else {
+            sum_back += data[idx];
+            count_back++;
+        }
+    }
+    mean_fore = sum_fore / count_fore;
+    mean_back = sum_back / count_back;
+}
 
-    //need to select some pixels. we need some buffer room in the image.
-    const size_t total_required_pixels = std::min((size_t) 10*512*512,localIntensityScale.mesh.size()/2);
-    const size_t delta = localIntensityScale.mesh.size()/total_required_pixels - 1;
 
-    std::vector<T> lis_buffer(total_required_pixels);
-    std::vector<S> grad_buffer(total_required_pixels);
+/**
+ * Compute threshold value by Li's iterative Minimum Cross Entropy method [1]
+ *
+ * Note: it is assumed that the input elements are non-negative (as is the case for the gradient and local intensity
+ * scale). To apply the method to data that may be negative, subtract the minimum value from the input, and then add
+ * that value to the computed threshold.
+ *
+ * [1] Li, C. H., & Tam, P. K. S. (1998). "An iterative algorithm for minimum cross entropy thresholding."
+ *     Pattern recognition letters, 19(8), 771-776.
+ */
+template<typename T>
+float threshold_li(const std::vector<T> &input) {
 
-    uint64_t counter = 0;
-    uint64_t counter_sampled = 0;
+    if(input.empty()) { return 0.f; }     // return 0 if input is empty
 
-    while((counter < localIntensityScale.mesh.size()) && (counter_sampled < total_required_pixels)){
+    const T image_min = *std::min_element(input.begin(), input.end());
+    const T image_max = *std::max_element(input.begin(), input.end());
 
-        lis_buffer[counter_sampled] = localIntensityScale.mesh[counter];
-        grad_buffer[counter_sampled] = grad.mesh[counter];
+    if(image_min == image_max) { return image_min; }  // if all inputs are equal, return that value
 
-        counter+=delta;
-        counter_sampled++;
+    float tolerance = 0.5f;   // tolerance 0.5 should suffice for integer inputs
 
+    // For floating point inputs we set the tolerance to the lesser of 0.01 and a fraction of the data range
+    // This could be improved, by instead taking e.g. half the smallest difference between any two non-equal elements
+    if(std::is_floating_point<T>::value) {
+        float range = image_max - image_min;
+        tolerance = std::min(0.01f, range*1e-4f);
     }
 
+    // Take the mean of the input as initial value
+    float t_next = std::accumulate(input.begin(), input.end(), 0.f) / (float) input.size();
+    float t_curr = -2.f * tolerance; //this ensures that the first iteration is performed, since t_next is non-negative
 
-    //float min_lis = *std::min_element(lis_buffer.begin(),lis_buffer.end());
-    float max_lis = *std::max_element(lis_buffer.begin(),lis_buffer.end());
-
-    std::vector<uint64_t> hist_lis;
-    hist_lis.resize(std::ceil(max_lis)+1,0);
-
-    for (uint64_t i = 0; i < total_required_pixels; ++i) {
-        auto lis_val = std::floor(lis_buffer[i]);
-        hist_lis[lis_val]++;
+    // For integer inputs, we have to ensure a large enough initial value, such that mean_back > 0
+    if(!std::is_floating_point<T>::value) {
+        // if initial value is <1, try to increase it to 1.5 unless the range is too narrow
+        if(t_next < 1.f) {
+            t_next = std::min(1.5f, (image_min+image_max)/2.f);
+        }
     }
 
-    //Then find 5% therhold, and take the grad values from that.
-    uint64_t prop_values = 0.05*total_required_pixels;
+    // Perform Li iterations until convergence
+    while(std::abs(t_next - t_curr) > tolerance) {
+        t_curr = t_next;
 
-    uint64_t cumsum = 0;
-    uint64_t freq_val=0;
-    while (cumsum < prop_values){
-        cumsum+= hist_lis[freq_val];
-        freq_val++;
+        // Compute averages of values above and below the current threshold
+        float mean_back, mean_fore;
+        compute_means(input, t_curr, mean_back, mean_fore);
+
+        // Handle the edge case where all values < t_curr are 0
+        if(mean_back == 0) {
+            std::wcout << "log(0) encountered in APRConverter::threshold_li, returning current threshold" << std::endl;
+            return t_curr;
+        }
 
+        // Update the threshold (one-point iteration)
+        t_next = (mean_fore - mean_back) / (std::log(mean_fore) - std::log(mean_back));
     }
+    return t_next;
+}
 
-    std::vector<S> grad_hist;
-    float grad_max = *std::max_element(grad_buffer.begin(),grad_buffer.end());
-    //float grad_min = *std::max_element(grad_buffer.begin(),grad_buffer.end());
 
-    grad_hist.resize(std::ceil(grad_max));
-    uint64_t grad_counter = 0;
+template<typename ImageType>
+template<typename T,typename S>
+void APRConverter<ImageType>::autoParametersLiEntropy(const PixelData<T> &image,
+                                                      const PixelData<T> &localIntensityScale,
+                                                      const PixelData<S> &grad) {
+
+    fine_grained_timer.start_timer("autoparameters: subsample buffers");
+
+    std::vector<S> grad_subsampled;
+    std::vector<T> lis_subsampled;
+
+    {   // intentional scope
+        /// First we extract the gradient and local intensity scale values at all locations where the image
+        /// intensity exceeds the intensity threshold. This allows better adaptation in certain cases, e.g.
+        /// when there is significant background AND signal noise/autofluorescence (provided that par.Ip_th
+        /// is set appropriately).
+        std::vector<S> grad_foreground(grad.mesh.size());
+        std::vector<T> lis_foreground(localIntensityScale.mesh.size());
+
+        const auto threshold = par.Ip_th + bspline_offset;
+        size_t counter = 0;
+        for(size_t idx = 0; idx < grad.mesh.size(); ++idx) {
+            if(image.mesh[idx] > threshold) {
+                grad_foreground[counter] = grad.mesh[idx];
+                lis_foreground[counter] = localIntensityScale.mesh[idx];
+                counter++;
+            }
+        }
 
-    for (uint64_t i = 0; i < total_required_pixels; ++i) {
-        auto val = lis_buffer[i];
+        const size_t num_foreground_pixels = counter;
 
-        if(val <= freq_val){
-            auto grad_val = grad_buffer[i];
-            grad_hist[std::floor(grad_val)]++;
-            grad_counter++;
-        }
-    }
+        grad_foreground.resize(num_foreground_pixels); //setting size to non-zero elements.
+        lis_foreground.resize(num_foreground_pixels);
 
-    auto max_it = std::max_element(grad_hist.begin(),grad_hist.end());
-    uint64_t mode = std::distance(grad_hist.begin(),max_it);
+        /// Then we uniformly subsample these signals, as we typically don't need all elements to compute the thresholds
+        const size_t num_elements = std::min((size_t)32*512*512, num_foreground_pixels); //arbitrary number.
+        const size_t delta = num_foreground_pixels / num_elements;
+        grad_subsampled.resize(num_elements);
+        lis_subsampled.resize(num_elements);
 
-    float grad_th = std::round(4*mode); //magic numbers.
+#ifdef HAVE_OPENMP
+#pragma omp parallel for default(shared)
+#endif
+        for(size_t idx = 0; idx < num_elements; ++idx) {
+            grad_subsampled[idx] = grad_foreground[idx*delta];
+            lis_subsampled[idx] = lis_foreground[idx*delta];
+        }
+    }
+    fine_grained_timer.stop_timer();
 
-    par.grad_th = grad_th;
-    par.sigma_th = freq_val;
-    par.sigma_th_max = 1;
+    fine_grained_timer.start_timer("autoparameters: compute gradient threshold");
+    par.grad_th = threshold_li(grad_subsampled);
+    fine_grained_timer.stop_timer();
 
-    std::cout << "Used parameters: " << std::endl;
-    std::cout << "I_th: " << par.Ip_th << std::endl;
-    std::cout << "sigma_th: " << par.sigma_th << std::endl;
-    std::cout << "grad_th: " << par.grad_th << std::endl;
-    std::cout << "relative error (E): " << par.rel_error << std::endl;
-    std::cout << "lambda: " << par.lambda << std::endl;
+    fine_grained_timer.start_timer("autoparameters: compute sigma threshold");
+    par.sigma_th = threshold_li(lis_subsampled);
+    fine_grained_timer.stop_timer();
 
+    if(verbose) {
+        std::cout << "Automatic parameter tuning found sigma_th = " << par.sigma_th <<
+                     " and grad_th = " << par.grad_th << std::endl;
+    }
 }
 
 

test/APRTest.cpp

@@ -287,10 +287,18 @@ bool test_auto_parameters(TestData& test_data){
 
     aprConverter.get_apr(apr, test_data.img_original);
 
-    //This test doesn't do anything except check if it runs.
+    auto par = apr.get_apr_parameters();
 
-    return success;
+    //checks if run, and the values are positive
+    if(par.sigma_th <= 0){
+        success = false;
+    }
 
+    if(par.grad_th <= 0){
+        success = false;
+    }
+
+    return success;
 
 }