Skip to content

Commit f750ef5

Browse files
joeljonssonjoeljonsson
committed
move autoparameter methods to separate file
1 parent e34fa89 commit f750ef5

File tree

2 files changed

+162
-153
lines changed

2 files changed

+162
-153
lines changed

‎src/algorithm/APRConverter.hpp‎

Lines changed: 2 additions & 153 deletions
Original file line numberDiff line numberDiff line change
@@ -11,6 +11,7 @@
1111

1212
#include <list>
1313

14+
#include "AutoParameters.hpp"
1415
#include "data_structures/APR/APR.hpp"
1516
#include "data_structures/APR/particles/ParticleData.hpp"
1617
#include "data_structures/Mesh/PixelData.hpp"
@@ -124,9 +125,6 @@ class APRConverter {
124125

125126
void generateDatastructures(APR& aAPR);
126127

127-
template<typename T,typename S>
128-
void autoParametersLiEntropy(const PixelData<T> &image, const PixelData<T> &localIntensityScale, const PixelData<S> &grad);
129-
130128
template<typename T>
131129
bool check_input_dimensions(PixelData<T> &input_image);
132130

@@ -579,7 +577,7 @@ inline bool APRConverter<ImageType>::get_apr_cpu(APR &aAPR, PixelData<T> &input_
579577

580578
if (par.auto_parameters) {
581579
method_timer.start_timer("autoParameters");
582-
autoParametersLiEntropy(local_scale_temp2, local_scale_temp, grad_temp);
580+
autoParametersLiEntropy(par, local_scale_temp2, local_scale_temp, grad_temp, bspline_offset, verbose);
583581
aAPR.parameters = par;
584582
method_timer.stop_timer();
585583
}
@@ -624,155 +622,6 @@ inline bool APRConverter<ImageType>::get_apr(APR &aAPR, PixelData<T> &input_imag
624622
}
625623

626624

627-
template<typename T>
628-
void compute_means(const std::vector<T>& data, float threshold, float& mean_back, float& mean_fore) {
629-
float sum_fore=0.f, sum_back=0.f;
630-
size_t count_fore=0, count_back=0;
631-
632-
#ifdef HAVE_OPENMP
633-
#pragma omp parallel for default(none) shared(data, threshold) reduction(+:sum_fore, sum_back, count_fore, count_back)
634-
#endif
635-
for(size_t idx = 0; idx < data.size(); ++idx) {
636-
if(data[idx] > threshold) {
637-
sum_fore += data[idx];
638-
count_fore++;
639-
} else {
640-
sum_back += data[idx];
641-
count_back++;
642-
}
643-
}
644-
mean_fore = sum_fore / count_fore;
645-
mean_back = sum_back / count_back;
646-
}
647-
648-
649-
/**
650-
* Compute threshold value by Li's iterative Minimum Cross Entropy method [1]
651-
*
652-
* Note: it is assumed that the input elements are non-negative (as is the case for the gradient and local intensity
653-
* scale). To apply the method to data that may be negative, subtract the minimum value from the input, and then add
654-
* that value to the computed threshold.
655-
*
656-
* [1] Li, C. H., & Tam, P. K. S. (1998). "An iterative algorithm for minimum cross entropy thresholding."
657-
* Pattern recognition letters, 19(8), 771-776.
658-
*/
659-
template<typename T>
660-
float threshold_li(const std::vector<T> &input) {
661-
662-
if(input.empty()) { return 0.f; } // return 0 if input is empty
663-
664-
const T image_min = *std::min_element(input.begin(), input.end());
665-
const T image_max = *std::max_element(input.begin(), input.end());
666-
667-
if(image_min == image_max) { return image_min; } // if all inputs are equal, return that value
668-
669-
float tolerance = 0.5f; // tolerance 0.5 should suffice for integer inputs
670-
671-
// For floating point inputs we set the tolerance to the lesser of 0.01 and a fraction of the data range
672-
// This could be improved, by instead taking e.g. half the smallest difference between any two non-equal elements
673-
if(std::is_floating_point<T>::value) {
674-
float range = image_max - image_min;
675-
tolerance = std::min(0.01f, range*1e-4f);
676-
}
677-
678-
// Take the mean of the input as initial value
679-
float t_next = std::accumulate(input.begin(), input.end(), 0.f) / (float) input.size();
680-
float t_curr = -2.f * tolerance; //this ensures that the first iteration is performed, since t_next is non-negative
681-
682-
// For integer inputs, we have to ensure a large enough initial value, such that mean_back > 0
683-
if(!std::is_floating_point<T>::value) {
684-
// if initial value is <1, try to increase it to 1.5 unless the range is too narrow
685-
if(t_next < 1.f) {
686-
t_next = std::min(1.5f, (image_min+image_max)/2.f);
687-
}
688-
}
689-
690-
// Perform Li iterations until convergence
691-
while(std::abs(t_next - t_curr) > tolerance) {
692-
t_curr = t_next;
693-
694-
// Compute averages of values above and below the current threshold
695-
float mean_back, mean_fore;
696-
compute_means(input, t_curr, mean_back, mean_fore);
697-
698-
// Handle the edge case where all values < t_curr are 0
699-
if(mean_back == 0) {
700-
std::wcout << "log(0) encountered in APRConverter::threshold_li, returning current threshold" << std::endl;
701-
return t_curr;
702-
}
703-
704-
// Update the threshold (one-point iteration)
705-
t_next = (mean_fore - mean_back) / (std::log(mean_fore) - std::log(mean_back));
706-
}
707-
return t_next;
708-
}
709-
710-
711-
template<typename ImageType>
712-
template<typename T,typename S>
713-
void APRConverter<ImageType>::autoParametersLiEntropy(const PixelData<T> &image,
714-
const PixelData<T> &localIntensityScale,
715-
const PixelData<S> &grad) {
716-
717-
fine_grained_timer.start_timer("autoparameters: subsample buffers");
718-
719-
std::vector<S> grad_subsampled;
720-
std::vector<T> lis_subsampled;
721-
722-
{ // intentional scope
723-
/// First we extract the gradient and local intensity scale values at all locations where the image
724-
/// intensity exceeds the intensity threshold. This allows better adaptation in certain cases, e.g.
725-
/// when there is significant background AND signal noise/autofluorescence (provided that par.Ip_th
726-
/// is set appropriately).
727-
std::vector<S> grad_foreground(grad.mesh.size());
728-
std::vector<T> lis_foreground(localIntensityScale.mesh.size());
729-
730-
const auto threshold = par.Ip_th + bspline_offset;
731-
size_t counter = 0;
732-
for(size_t idx = 0; idx < grad.mesh.size(); ++idx) {
733-
if(image.mesh[idx] > threshold) {
734-
grad_foreground[counter] = grad.mesh[idx];
735-
lis_foreground[counter] = localIntensityScale.mesh[idx];
736-
counter++;
737-
}
738-
}
739-
740-
const size_t num_foreground_pixels = counter;
741-
742-
grad_foreground.resize(num_foreground_pixels); //setting size to non-zero elements.
743-
lis_foreground.resize(num_foreground_pixels);
744-
745-
/// Then we uniformly subsample these signals, as we typically don't need all elements to compute the thresholds
746-
const size_t num_elements = std::min((size_t)32*512*512, num_foreground_pixels); //arbitrary number.
747-
const size_t delta = num_foreground_pixels / num_elements;
748-
grad_subsampled.resize(num_elements);
749-
lis_subsampled.resize(num_elements);
750-
751-
#ifdef HAVE_OPENMP
752-
#pragma omp parallel for default(shared)
753-
#endif
754-
for(size_t idx = 0; idx < num_elements; ++idx) {
755-
grad_subsampled[idx] = grad_foreground[idx*delta];
756-
lis_subsampled[idx] = lis_foreground[idx*delta];
757-
}
758-
}
759-
fine_grained_timer.stop_timer();
760-
761-
fine_grained_timer.start_timer("autoparameters: compute gradient threshold");
762-
par.grad_th = threshold_li(grad_subsampled);
763-
fine_grained_timer.stop_timer();
764-
765-
fine_grained_timer.start_timer("autoparameters: compute sigma threshold");
766-
par.sigma_th = threshold_li(lis_subsampled);
767-
fine_grained_timer.stop_timer();
768-
769-
if(verbose) {
770-
std::cout << "Automatic parameter tuning found sigma_th = " << par.sigma_th <<
771-
" and grad_th = " << par.grad_th << std::endl;
772-
}
773-
}
774-
775-
776625
/**
777626
* Checks if the memory dimension (y) is filled
778627
*/

‎src/algorithm/AutoParameters.hpp‎

Lines changed: 160 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,160 @@
1+
//
2+
// Created by joel on 26.04.22.
3+
//
4+
5+
#ifndef APR_AUTOPARAMETERS_HPP
6+
#define APR_AUTOPARAMETERS_HPP
7+
8+
#include "APRParameters.hpp"
9+
#include "data_structures/Mesh/PixelData.hpp"
10+
11+
#include <stdio.h>
12+
#include <vector>
13+
#include <algorithm>
14+
#include <numeric>
15+
#include <cmath>
16+
17+
template<typename T>
18+
void compute_means(const std::vector<T>& data, float threshold, float& mean_back, float& mean_fore) {
19+
float sum_fore=0.f, sum_back=0.f;
20+
size_t count_fore=0, count_back=0;
21+
22+
#ifdef HAVE_OPENMP
23+
#pragma omp parallel for default(none) shared(data, threshold) reduction(+:sum_fore, sum_back, count_fore, count_back)
24+
#endif
25+
for(size_t idx = 0; idx < data.size(); ++idx) {
26+
if(data[idx] > threshold) {
27+
sum_fore += data[idx];
28+
count_fore++;
29+
} else {
30+
sum_back += data[idx];
31+
count_back++;
32+
}
33+
}
34+
mean_fore = sum_fore / count_fore;
35+
mean_back = sum_back / count_back;
36+
}
37+
38+
39+
/**
40+
* Compute threshold value by Li's iterative Minimum Cross Entropy method [1]
41+
*
42+
* Note: it is assumed that the input elements are non-negative (as is the case for the gradient and local intensity
43+
* scale). To apply the method to data that may be negative, subtract the minimum value from the input, and then add
44+
* that value to the computed threshold.
45+
*
46+
* [1] Li, C. H., & Tam, P. K. S. (1998). "An iterative algorithm for minimum cross entropy thresholding."
47+
* Pattern recognition letters, 19(8), 771-776.
48+
*/
49+
template<typename T>
50+
float threshold_li(const std::vector<T> &input) {
51+
52+
if(input.empty()) { return 0.f; } // return 0 if input is empty
53+
54+
const T image_min = *std::min_element(input.begin(), input.end());
55+
const T image_max = *std::max_element(input.begin(), input.end());
56+
57+
if(image_min == image_max) { return image_min; } // if all inputs are equal, return that value
58+
59+
float tolerance = 0.5f; // tolerance 0.5 should suffice for integer inputs
60+
61+
// For floating point inputs we set the tolerance to the lesser of 0.01 and a fraction of the data range
62+
// This could be improved, by instead taking e.g. half the smallest difference between any two non-equal elements
63+
if(std::is_floating_point<T>::value) {
64+
float range = image_max - image_min;
65+
tolerance = std::min(0.01f, range*1e-4f);
66+
}
67+
68+
// Take the mean of the input as initial value
69+
float t_next = std::accumulate(input.begin(), input.end(), 0.f) / (float) input.size();
70+
float t_curr = -2.f * tolerance; //this ensures that the first iteration is performed, since t_next is non-negative
71+
72+
// For integer inputs, we have to ensure a large enough initial value, such that mean_back > 0
73+
if(!std::is_floating_point<T>::value) {
74+
// if initial value is <1, try to increase it to 1.5 unless the range is too narrow
75+
if(t_next < 1.f) {
76+
t_next = std::min(1.5f, (image_min+image_max)/2.f);
77+
}
78+
}
79+
80+
// Perform Li iterations until convergence
81+
while(std::abs(t_next - t_curr) > tolerance) {
82+
t_curr = t_next;
83+
84+
// Compute averages of values above and below the current threshold
85+
float mean_back, mean_fore;
86+
compute_means(input, t_curr, mean_back, mean_fore);
87+
88+
// Handle the edge case where all values < t_curr are 0
89+
if(mean_back == 0) {
90+
std::wcout << "log(0) encountered in APRConverter::threshold_li, returning current threshold" << std::endl;
91+
return t_curr;
92+
}
93+
94+
// Update the threshold (one-point iteration)
95+
t_next = (mean_fore - mean_back) / (std::log(mean_fore) - std::log(mean_back));
96+
}
97+
return t_next;
98+
}
99+
100+
101+
template<typename S, typename T, typename U>
102+
void autoParametersLiEntropy(APRParameters& par,
103+
const PixelData<S>& image,
104+
const PixelData<T>& localIntensityScale,
105+
const PixelData<U>& grad,
106+
const float bspline_offset,
107+
const bool verbose=false) {
108+
109+
std::vector<U> grad_subsampled;
110+
std::vector<T> lis_subsampled;
111+
112+
{ // intentional scope
113+
/// First we extract the gradient and local intensity scale values at all locations where the image
114+
/// intensity exceeds the intensity threshold. This allows better adaptation in certain cases, e.g.
115+
/// when there is significant background AND signal noise/autofluorescence (provided that par.Ip_th
116+
/// is set appropriately).
117+
std::vector<U> grad_foreground(grad.mesh.size());
118+
std::vector<T> lis_foreground(localIntensityScale.mesh.size());
119+
120+
const auto threshold = par.Ip_th + bspline_offset;
121+
size_t counter = 0;
122+
for(size_t idx = 0; idx < grad.mesh.size(); ++idx) {
123+
if(image.mesh[idx] > threshold) {
124+
grad_foreground[counter] = grad.mesh[idx];
125+
lis_foreground[counter] = localIntensityScale.mesh[idx];
126+
counter++;
127+
}
128+
}
129+
130+
const size_t num_foreground_pixels = counter;
131+
132+
grad_foreground.resize(num_foreground_pixels); //setting size to non-zero elements.
133+
lis_foreground.resize(num_foreground_pixels);
134+
135+
/// Then we uniformly subsample these signals, as we typically don't need all elements to compute the thresholds
136+
const size_t num_elements = std::min((size_t)32*512*512, num_foreground_pixels); //arbitrary number.
137+
const size_t delta = num_foreground_pixels / num_elements;
138+
grad_subsampled.resize(num_elements);
139+
lis_subsampled.resize(num_elements);
140+
141+
#ifdef HAVE_OPENMP
142+
#pragma omp parallel for default(shared)
143+
#endif
144+
for(size_t idx = 0; idx < num_elements; ++idx) {
145+
grad_subsampled[idx] = grad_foreground[idx*delta];
146+
lis_subsampled[idx] = lis_foreground[idx*delta];
147+
}
148+
}
149+
150+
/// Compute thresholds using Li's iterative minimum cross entropy algorithm
151+
par.grad_th = threshold_li(grad_subsampled);
152+
par.sigma_th = threshold_li(lis_subsampled);
153+
154+
if(verbose) {
155+
std::cout << "Automatic parameter tuning found sigma_th = " << par.sigma_th <<
156+
" and grad_th = " << par.grad_th << std::endl;
157+
}
158+
}
159+
160+
#endif //APR_AUTOPARAMETERS_HPP

0 commit comments

Comments
 (0)