39 #include <visp3/core/vpHistogram.h>
40 #include <visp3/core/vpImageTools.h>
41 #include <visp3/imgproc/vpImgproc.h>
45 bool isBimodal(
const std::vector<float> &hist_float)
49 for (
size_t cpt = 1; cpt < hist_float.size() - 1; cpt++) {
50 if (hist_float[cpt - 1] < hist_float[cpt] && hist_float[cpt] > hist_float[cpt + 1]) {
74 for (first = 0; first < (size_t)hist.
getSize() && hist[(
unsigned char)first] == 0; first++) {
78 for (last = (
size_t)hist.
getSize() - 1; last > first && hist[(
unsigned char)last] == 0; last--) {
87 std::vector<float> S(last + 1);
88 std::vector<float> W(last + 1);
90 S[0] = (float)hist[0];
92 for (
size_t i = std::max((
size_t)1, first); i <= last; i++) {
93 S[i] = S[i - 1] + hist[(
unsigned char)i];
94 W[i] = W[i - 1] + i * (float)hist[(
unsigned char)i];
99 float C = (float)(last - first);
100 std::vector<float> Smu(last + 1 - first);
102 for (
size_t i = 1; i < Smu.size(); i++) {
103 float mu = 1 / (1 + i / C);
104 Smu[i] = -mu * std::log(mu) - (1 - mu) * std::log(1 - mu);
108 int bestThreshold = 0;
109 float bestEntropy = std::numeric_limits<float>::max();
111 for (
size_t threshold = first; threshold <= last; threshold++) {
114 for (
size_t i = first; i <= threshold; i++) {
115 entropy += Smu[(size_t)std::abs((
int)i - mu)] * hist[(
unsigned char)i];
118 mu =
vpMath::round((W[last] - W[threshold]) / (S[last] - S[threshold]));
119 for (
size_t i = threshold + 1; i <= last; i++) {
120 entropy += Smu[(size_t)std::abs((
int)i - mu)] * hist[(
unsigned char)i];
123 if (bestEntropy > entropy) {
124 bestEntropy = entropy;
125 bestThreshold = (int)threshold;
129 return bestThreshold;
132 int computeThresholdIntermodes(
const vpHistogram &hist)
152 std::vector<float> hist_float(hist.
getSize());
153 for (
unsigned int cpt = 0; cpt < hist.
getSize(); cpt++) {
154 hist_float[cpt] = (float)hist[cpt];
158 while (!isBimodal(hist_float)) {
160 for (
size_t cpt = 1; cpt < hist_float.size() - 1; cpt++) {
161 hist_float[cpt] = (hist_float[cpt - 1] + hist_float[cpt] + hist_float[cpt + 1]) / 3;
165 hist_float[0] = (hist_float[0] + hist_float[1]) / 2.0f;
168 hist_float[hist_float.size() - 1] = (hist_float.size() - 2 + hist_float.size() - 1) / 2.0f;
173 std::cerr <<
"Intermodes Threshold not found after 10000 iterations!" << std::endl;
180 for (
size_t cpt = 1; cpt < hist_float.size() - 1; cpt++) {
181 if (hist_float[cpt - 1] < hist_float[cpt] && hist_float[cpt] > hist_float[cpt + 1]) {
187 return (
int)std::floor(tt / 2.0);
190 int computeThresholdIsoData(
const vpHistogram &hist,
unsigned int imageSize)
196 std::vector<float> cumsum(hist.
getSize(), 0.0f);
197 std::vector<float> sum_ip(hist.
getSize(), 0.0f);
198 cumsum[0] = (float)hist[0];
199 for (
unsigned int cpt = 1; cpt < hist.
getSize(); cpt++) {
200 sum_ip[cpt] = cpt * (float)hist[cpt] + sum_ip[cpt - 1];
201 cumsum[cpt] = (float)hist[cpt] + cumsum[cpt - 1];
207 float MBT = sum_ip[(size_t)(T - 2)] / cumsum[(size_t)(T - 2)];
208 float MAT = (sum_ip.back() - sum_ip[(size_t)(T - 1)]) / (cumsum.back() - cumsum[(size_t)(T - 1)]);
213 while (std::abs(T2 - T) >= 1) {
214 MBT = sum_ip[(size_t)(T2 - 2)] / cumsum[(size_t)(T2 - 2)];
215 MAT = (sum_ip.back() - sum_ip[(size_t)(T2 - 1)]) / (cumsum.back() - cumsum[(size_t)(T2 - 1)]);
225 int computeThresholdMean(
const vpHistogram &hist,
unsigned int imageSize)
231 for (
unsigned int cpt = 0; cpt < hist.
getSize(); cpt++) {
232 sum_ip += cpt * (float)hist[cpt];
235 return (
int)std::floor(sum_ip / imageSize);
238 int computeThresholdOtsu(
const vpHistogram &hist,
unsigned int imageSize)
245 float sum_ip_all[256];
246 for (
int cpt = 0; cpt < (int)hist.
getSize(); cpt++) {
247 mu_T += cpt * (float)hist[cpt];
248 sum_ip_all[cpt] = mu_T;
252 float w_B = 0.0f, w_F = 0.0f;
254 float max_sigma_b = 0.0f;
257 for (
int cpt = 0; cpt < 256; cpt++) {
259 if (
vpMath::nul(w_B, std::numeric_limits<float>::epsilon())) {
263 w_F = (int)imageSize - w_B;
264 if (
vpMath::nul(w_F, std::numeric_limits<float>::epsilon())) {
269 float mu_B = sum_ip_all[cpt] / (float)w_B;
270 float mu_F = (mu_T - sum_ip_all[cpt]) / (
float)w_F;
273 float sigma_b_sqr = w_B * w_F * (mu_B - mu_F) * (mu_B - mu_F);
275 if (sigma_b_sqr >= max_sigma_b) {
277 max_sigma_b = sigma_b_sqr;
292 int left_bound = -1, right_bound = -1, max_idx = -1, max_value = 0;
294 for (
int cpt = 0; cpt < (int)hist.
getSize(); cpt++) {
295 if (left_bound == -1 && hist[cpt] > 0) {
296 left_bound = (int)cpt;
299 if (right_bound == -1 && hist[(
int)hist.
getSize() - 1 - cpt] > 0) {
300 right_bound = (int)hist.
getSize() - 1 - cpt;
303 if ((
int)hist[cpt] > max_value) {
304 max_value = (int)hist[cpt];
310 left_bound = left_bound > 0 ? left_bound - 1 : left_bound;
311 right_bound = right_bound < (int)hist.
getSize() - 1 ? right_bound + 1 : right_bound;
315 if (max_idx - left_bound < right_bound - max_idx) {
319 int cpt_left = 0, cpt_right = (int)hist.
getSize() - 1;
320 for (; cpt_left < cpt_right; cpt_left++, cpt_right--) {
321 unsigned int temp = hist[cpt_left];
322 hist.
set(cpt_left, hist[cpt_right]);
323 hist.
set(cpt_right, temp);
326 left_bound = (int)hist.
getSize() - 1 - right_bound;
327 max_idx = (int)hist.
getSize() - 1 - max_idx;
337 float a = (float)max_value;
338 float b = (float)(left_bound - max_idx);
339 float max_dist = 0.0f;
341 for (
int cpt = left_bound + 1; cpt <= max_idx; cpt++) {
342 float dist = a * cpt + b * hist[cpt];
344 if (dist > max_dist) {
352 threshold = (int)hist.
getSize() - 1 - threshold;
362 const unsigned char backgroundValue,
const unsigned char foregroundValue)
374 threshold = computeThresholdHuang(histogram);
378 threshold = computeThresholdIntermodes(histogram);
382 threshold = computeThresholdIsoData(histogram, I.
getSize());
386 threshold = computeThresholdMean(histogram, I.
getSize());
390 threshold = computeThresholdOtsu(histogram, I.
getSize());
394 threshold = computeThresholdTriangle(histogram);
401 if (threshold != -1) {
Class to compute a gray level image histogram.
void set(const unsigned char level, unsigned int value)
unsigned int getSize() const
static bool nul(double x, double threshold=0.001)
static int round(double x)
VISP_EXPORT unsigned char autoThreshold(vpImage< unsigned char > &I, const vp::vpAutoThresholdMethod &method, const unsigned char backgroundValue=0, const unsigned char foregroundValue=255)
@ AUTO_THRESHOLD_INTERMODES
@ AUTO_THRESHOLD_TRIANGLE