Visual Servoing Platform  version 3.6.1 under development (2024-05-21)
vpThreshold.cpp
1 /*
2  * ViSP, open source Visual Servoing Platform software.
3  * Copyright (C) 2005 - 2023 by Inria. All rights reserved.
4  *
5  * This software is free software; you can redistribute it and/or modify
6  * it under the terms of the GNU General Public License as published by
7  * the Free Software Foundation; either version 2 of the License, or
8  * (at your option) any later version.
9  * See the file LICENSE.txt at the root directory of this source
10  * distribution for additional information about the GNU GPL.
11  *
12  * For using ViSP with software that can not be combined with the GNU
13  * GPL, please contact Inria about acquiring a ViSP Professional
14  * Edition License.
15  *
16  * See https://visp.inria.fr for more information.
17  *
18  * This software was developed at:
19  * Inria Rennes - Bretagne Atlantique
20  * Campus Universitaire de Beaulieu
21  * 35042 Rennes Cedex
22  * France
23  *
24  * If you have questions regarding the use of this file, please contact
25  * Inria at visp@inria.fr
26  *
27  * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
28  * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
29  *
30  * Description:
31  * Automatic thresholding functions.
32  */
33 
39 #include <visp3/core/vpHistogram.h>
40 #include <visp3/core/vpImageTools.h>
41 #include <visp3/imgproc/vpImgproc.h>
42 
43 namespace
44 {
45 bool isBimodal(const std::vector<float> &hist_float)
46 {
47  int modes = 0;
48 
49  size_t hist_float_size_m_1 = hist_float.size() - 1;
50  for (size_t cpt = 1; cpt < hist_float_size_m_1; ++cpt) {
51  if ((hist_float[cpt - 1] < hist_float[cpt]) && (hist_float[cpt] > hist_float[cpt + 1])) {
52  ++modes;
53  }
54 
55  if (modes > 2) {
56  return false;
57  }
58  }
59 
60  return (modes == 2);
61 }
62 
63 int computeThresholdHuang(const vpHistogram &hist)
64 {
65  // Code ported from the AutoThreshold ImageJ plugin:
66  // Implements Huang's fuzzy thresholding method
67  // Uses Shannon's entropy function (one can also use Yager's entropy
68  // function) Huang L.-K. and Wang M.-J.J. (1995) "Image Thresholding by
69  // Minimizing the Measures of Fuzziness" Pattern Recognition, 28(1): 41-51
70  // Reimplemented (to handle 16-bit efficiently) by Johannes Schindelin Jan
71  // 31, 2011
72 
73  // Find first and last non-empty bin
74  size_t first, last;
75  size_t hist_size = hist.getSize();
76  for (first = 0; (first < hist_size) && (hist[static_cast<unsigned char>(first)] == 0); ++first) {
77  // do nothing
78  }
79 
80  for (last = (hist_size - 1); (last > first) && (hist[static_cast<unsigned char>(last)] == 0); --last) {
81  // do nothing
82  }
83 
84  if (first == last) {
85  return 0;
86  }
87 
88  // Calculate the cumulative density and the weighted cumulative density
89  std::vector<float> S(last + 1);
90  std::vector<float> W(last + 1);
91 
92  S[0] = static_cast<float>(hist[0]);
93  W[0] = 0.0f;
94  for (size_t i = std::max<size_t>(static_cast<size_t>(1), first); i <= last; ++i) {
95  S[i] = S[i - 1] + hist[static_cast<unsigned char>(i)];
96  W[i] = W[i - 1] + (i * static_cast<float>(hist[static_cast<unsigned char>(i)]));
97  }
98 
99  // Precalculate the summands of the entropy given the absolute difference x
100  // - mu (integral)
101  float C = static_cast<float>(last - first);
102  std::vector<float> Smu((last + 1) - first);
103 
104  size_t smu_size = Smu.size();
105  for (size_t i = 1; i < smu_size; ++i) {
106  float mu = 1 / (1 + (i / C));
107  Smu[i] = (-mu * std::log(mu)) - ((1 - mu) * std::log(1 - mu));
108  }
109 
110  // Calculate the threshold
111  int bestThreshold = 0;
112  float bestEntropy = std::numeric_limits<float>::max();
113 
114  for (size_t threshold = first; threshold <= last; ++threshold) {
115  float entropy = 0;
116  int mu = vpMath::round(W[threshold] / S[threshold]);
117  for (size_t i = first; i <= threshold; ++i) {
118  entropy += Smu[static_cast<size_t>(std::abs(static_cast<int>(i) - mu))] * hist[static_cast<unsigned char>(i)];
119  }
120 
121  mu = vpMath::round((W[last] - W[threshold]) / (S[last] - S[threshold]));
122  for (size_t i = threshold + 1; i <= last; ++i) {
123  entropy += Smu[static_cast<size_t>(std::abs(static_cast<int>(i) - mu))] * hist[static_cast<unsigned char>(i)];
124  }
125 
126  if (bestEntropy > entropy) {
127  bestEntropy = entropy;
128  bestThreshold = static_cast<int>(threshold);
129  }
130  }
131 
132  return bestThreshold;
133 }
134 
135 int computeThresholdIntermodes(const vpHistogram &hist)
136 {
137  if (hist.getSize() < 3) {
138  return -1;
139  }
140 
141  // Code based on the AutoThreshold ImageJ plugin:
142  // J. M. S. Prewitt and M. L. Mendelsohn, "The analysis of cell images," in
143  // Annals of the New York Academy of Sciences, vol. 128, pp. 1035-1053,
144  // 1966. ported to ImageJ plugin by G.Landini from Antti Niemisto's Matlab
145  // code (GPL) Original Matlab code Copyright (C) 2004 Antti Niemisto See
146  // http://www.cs.tut.fi/~ant/histthresh/ for an excellent slide presentation
147  // and the original Matlab code.
148  //
149  // Assumes a bimodal histogram. The histogram needs is smoothed (using a
150  // running average of size 3, iteratively) until there are only two local
151  // maxima. j and k Threshold t is (j+k)/2. Images with histograms having
152  // extremely unequal peaks or a broad and flat valley are unsuitable for this
153  // method.
154 
155  std::vector<float> hist_float(hist.getSize());
156  unsigned int hist_size = hist.getSize();
157  for (unsigned int cpt = 0; cpt < hist_size; ++cpt) {
158  hist_float[cpt] = static_cast<float>(hist[cpt]);
159  }
160 
161  int iter = 0;
162  while (!isBimodal(hist_float)) {
163  // Smooth with a 3 point running mean filter
164  size_t hist_float_size_m_1 = hist_float.size() - 1;
165  for (size_t cpt = 1; cpt < hist_float_size_m_1; ++cpt) {
166  hist_float[cpt] = (hist_float[cpt - 1] + hist_float[cpt] + hist_float[cpt + 1]) / 3;
167  }
168 
169  // First value
170  hist_float[0] = (hist_float[0] + hist_float[1]) / 2.0f;
171 
172  // Last value
173  hist_float[hist_float.size() - 1] = (((hist_float.size() - 2) + hist_float.size()) - 1) / 2.0f;
174 
175  ++iter;
176 
177  if (iter > 10000) {
178  std::cerr << "Intermodes Threshold not found after 10000 iterations!" << std::endl;
179  return -1;
180  }
181  }
182 
183  // The threshold is the mean between the two peaks.
184  int tt = 0;
185  size_t hist_float_size_m_1 = hist_float.size() - 1;
186  for (size_t cpt = 1; cpt < hist_float_size_m_1; ++cpt) {
187  if ((hist_float[cpt - 1] < hist_float[cpt]) && (hist_float[cpt] > hist_float[cpt + 1])) {
188  // Mode
189  tt += static_cast<int>(cpt);
190  }
191  }
192 
193  return static_cast<int>(std::floor(tt / 2.0)); // vpMath round of tt div 2.0
194 }
195 
196 int computeThresholdIsoData(const vpHistogram &hist, unsigned int imageSize)
197 {
198  int threshold = 0;
199 
200  // Code based on BSD Matlab isodata implementation by zephyr
201  // STEP 1: Compute mean intensity of image from histogram, set T=mean(I)
202  std::vector<float> cumsum(hist.getSize(), 0.0f);
203  std::vector<float> sum_ip(hist.getSize(), 0.0f);
204  cumsum[0] = static_cast<float>(hist[0]);
205 
206  unsigned int hist_size = hist.getSize();
207  for (unsigned int cpt = 1; cpt < hist_size; ++cpt) {
208  sum_ip[cpt] = (cpt * static_cast<float>(hist[cpt])) + sum_ip[cpt - 1];
209  cumsum[cpt] = static_cast<float>(hist[cpt]) + cumsum[cpt - 1];
210  }
211 
212  int T = vpMath::round(sum_ip[255] / imageSize);
213 
214  // STEP 2: compute Mean above T (MAT) and Mean below T (MBT) using T from
215  float MBT = sum_ip[static_cast<size_t>(T - 2)] / cumsum[static_cast<size_t>(T - 2)];
216  float MAT = (sum_ip.back() - sum_ip[static_cast<size_t>(T - 1)]) / (cumsum.back() - cumsum[static_cast<size_t>(T - 1)]);
217 
218  int T2 = vpMath::round((MAT + MBT) / 2.0f);
219 
220  //% STEP 3 to n: repeat step 2 if T(i)~=T(i-1)
221  while (std::abs(T2 - T) >= 1) {
222  MBT = sum_ip[static_cast<size_t>(T2 - 2)] / cumsum[static_cast<size_t>(T2 - 2)];
223  MAT = (sum_ip.back() - sum_ip[static_cast<size_t>(T2 - 1)]) / (cumsum.back() - cumsum[static_cast<size_t>(T2 - 1)]);
224 
225  T = T2;
226  T2 = vpMath::round((MAT + MBT) / 2.0f);
227  threshold = T2;
228  }
229 
230  return threshold;
231 }
232 
233 int computeThresholdMean(const vpHistogram &hist, unsigned int imageSize)
234 {
235  // C. A. Glasbey, "An analysis of histogram-based thresholding algorithms,"
236  // CVGIP: Graphical Models and Image Processing, vol. 55, pp. 532-537, 1993.
237  // The threshold is the mean of the greyscale data
238  float sum_ip = 0.0f;
239  unsigned int hist_size = hist.getSize();
240  for (unsigned int cpt = 0; cpt < hist_size; ++cpt) {
241  sum_ip += cpt * static_cast<float>(hist[cpt]);
242  }
243 
244  return static_cast<int>(std::floor(sum_ip / imageSize));
245 }
246 
247 int computeThresholdOtsu(const vpHistogram &hist, unsigned int imageSize)
248 {
249  // Otsu, N (1979), "A threshold selection method from gray-level
250  // histograms", IEEE Trans. Sys., Man., Cyber. 9: 62-66,
251  // doi:10.1109/TSMC.1979.4310076
252 
253  float mu_T = 0.0f;
254  float sum_ip_all[256];
255  int hist_size = static_cast<int>(hist.getSize());
256  for (int cpt = 0; cpt < hist_size; ++cpt) {
257  mu_T += cpt * static_cast<float>(hist[cpt]);
258  sum_ip_all[cpt] = mu_T;
259  }
260 
261  // Weight Background / Foreground
262  float w_B = 0.0f, w_F = 0.0f;
263 
264  float max_sigma_b = 0.0f;
265  int threshold = 0;
266 
267  bool w_f_eq_nul = false;
268  for (int cpt = 0; (cpt < 256) && (w_f_eq_nul == false); ++cpt) {
269  w_B += hist[cpt];
270  bool w_b_eq_nul = vpMath::nul(w_B, std::numeric_limits<float>::epsilon());
271  if (w_b_eq_nul == false) {
272 
273  w_F = static_cast<int>(imageSize) - w_B;
274  w_f_eq_nul = vpMath::nul(w_F, std::numeric_limits<float>::epsilon());
275  if (w_f_eq_nul == false) {
276 
277  // Mean Background / Foreground
278  float mu_B = sum_ip_all[cpt] / static_cast<float>(w_B);
279  float mu_F = (mu_T - sum_ip_all[cpt]) / static_cast<float>(w_F);
280  // If there is a case where (w_B * w_F) exceed FLT_MAX, normalize
281  // histogram
282  float sigma_b_sqr = w_B * w_F * (mu_B - mu_F) * (mu_B - mu_F);
283 
284  if (sigma_b_sqr >= max_sigma_b) {
285  threshold = cpt;
286  max_sigma_b = sigma_b_sqr;
287  }
288  }
289  // else exit the loop
290  }
291  }
292 
293  return threshold;
294 }
295 
296 int computeThresholdTriangle(vpHistogram &hist)
297 {
298  int threshold = 0;
299 
300  // Zack, G. W., Rogers, W. E. and Latt, S. A., 1977,
301  // Automatic Measurement of Sister Chromatid Exchange Frequency,
302  // Journal of Histochemistry and Cytochemistry 25 (7), pp. 741-753
303 
304  int left_bound = -1, right_bound = -1, max_idx = -1, max_value = 0;
305  // Find max value index and left / right most index
306  int hist_size = static_cast<int>(hist.getSize());
307  for (int cpt = 0; cpt < hist_size; ++cpt) {
308  if ((left_bound == -1) && (hist[cpt] > 0)) {
309  left_bound = static_cast<int>(cpt);
310  }
311 
312  if ((right_bound == -1) && (hist[static_cast<int>(hist.getSize()) - 1 - cpt] > 0)) {
313  right_bound = static_cast<int>(hist.getSize()) - 1 - cpt;
314  }
315 
316  if (static_cast<int>(hist[cpt]) > max_value) {
317  max_value = static_cast<int>(hist[cpt]);
318  max_idx = cpt;
319  }
320  }
321 
322  // First / last index when hist(cpt) == 0
323  left_bound = left_bound > 0 ? (left_bound - 1) : left_bound;
324  right_bound = right_bound < (static_cast<int>(hist.getSize()) - 1) ? (right_bound + 1) : right_bound;
325 
326  // Use the largest bound
327  bool flip = false;
328  if ((max_idx - left_bound) < (right_bound - max_idx)) {
329  // Flip histogram to get the largest bound to the left
330  flip = true;
331 
332  int cpt_left = 0;
333  int cpt_right = static_cast<int>(hist.getSize()) - 1;
334  for (; cpt_left < cpt_right; ++cpt_left, --cpt_right) {
335  unsigned int temp = hist[cpt_left];
336  hist.set(cpt_left, hist[cpt_right]);
337  hist.set(cpt_right, temp);
338  }
339 
340  left_bound = static_cast<int>(hist.getSize()) - 1 - right_bound;
341  max_idx = static_cast<int>(hist.getSize()) - 1 - max_idx;
342  }
343 
344  // Distance from a point to a line defined by two points:
345  //\textbf{distance} \left ( P_1, P_2, \left ( x_0,y_0 \right ) \right )
346  // = \frac{ \left | \left ( y_2-y_1 \right ) x_0 - \left ( x_2-x_1 \right )
347  // y_0 + x_2 y_1 - y_2 x_1 \right | }
348  // { \sqrt{ \left ( y_2 - y_1 \right )^{2} + \left ( x_2 - x_1 \right
349  // )^{2} } }
350  // Constants are ignored
351  float a = static_cast<float>(max_value); // y_2 - y_1
352  float b = static_cast<float>(left_bound - max_idx); //-(x_2 - x_1)
353  float max_dist = 0.0f;
354 
355  for (int cpt = left_bound + 1; cpt <= max_idx; ++cpt) {
356  float dist = (a * cpt) + (b * hist[cpt]);
357 
358  if (dist > max_dist) {
359  max_dist = dist;
360  threshold = cpt;
361  }
362  }
363  --threshold;
364 
365  if (flip) {
366  threshold = static_cast<int>(hist.getSize()) - 1 - threshold;
367  }
368 
369  return threshold;
370 }
371 } // namespace
372 
373 namespace vp
374 {
376  const unsigned char backgroundValue, const unsigned char foregroundValue)
377 {
378  if (I.getSize() == 0) {
379  return 0;
380  }
381 
382  // Compute image histogram
383  vpHistogram histogram(I);
384  int threshold = -1;
385 
386  switch (method) {
388  threshold = computeThresholdHuang(histogram);
389  break;
390 
392  threshold = computeThresholdIntermodes(histogram);
393  break;
394 
396  threshold = computeThresholdIsoData(histogram, I.getSize());
397  break;
398 
399  case AUTO_THRESHOLD_MEAN:
400  threshold = computeThresholdMean(histogram, I.getSize());
401  break;
402 
403  case AUTO_THRESHOLD_OTSU:
404  threshold = computeThresholdOtsu(histogram, I.getSize());
405  break;
406 
408  threshold = computeThresholdTriangle(histogram);
409  break;
410 
411  default:
412  break;
413  }
414 
415  if (threshold != -1) {
416  // Threshold
417  vpImageTools::binarise(I, static_cast<unsigned char>(threshold), static_cast<unsigned char>(255), backgroundValue, foregroundValue,
418  foregroundValue);
419  }
420 
421  return threshold;
422 }
423 };
Class to compute a gray level image histogram.
Definition: vpHistogram.h:108
unsigned getSize() const
Definition: vpHistogram.h:280
void set(const unsigned char level, unsigned int value)
Definition: vpHistogram.h:226
static void binarise(vpImage< Type > &I, Type threshold1, Type threshold2, Type value1, Type value2, Type value3, bool useLUT=true)
Definition: vpImageTools.h:470
unsigned int getSize() const
Definition: vpImage.h:224
static bool nul(double x, double threshold=0.001)
Definition: vpMath.h:445
static int round(double x)
Definition: vpMath.h:405
VISP_EXPORT unsigned char autoThreshold(vpImage< unsigned char > &I, const vp::vpAutoThresholdMethod &method, const unsigned char backgroundValue=0, const unsigned char foregroundValue=255)
vpAutoThresholdMethod
Definition: vpImgproc.h:67
@ AUTO_THRESHOLD_ISODATA
Definition: vpImgproc.h:76
@ AUTO_THRESHOLD_HUANG
Definition: vpImgproc.h:68
@ AUTO_THRESHOLD_INTERMODES
Definition: vpImgproc.h:72
@ AUTO_THRESHOLD_TRIANGLE
Definition: vpImgproc.h:88
@ AUTO_THRESHOLD_MEAN
Definition: vpImgproc.h:80
@ AUTO_THRESHOLD_OTSU
Definition: vpImgproc.h:84