Visual Servoing Platform  version 3.6.0 under development (2023-09-25)
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  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]) {
51  modes++;
52  }
53 
54  if (modes > 2) {
55  return false;
56  }
57  }
58 
59  return (modes == 2);
60 }
61 
62 int computeThresholdHuang(const vpHistogram &hist)
63 {
64  // Code ported from the AutoThreshold ImageJ plugin:
65  // Implements Huang's fuzzy thresholding method
66  // Uses Shannon's entropy function (one can also use Yager's entropy
67  // function) Huang L.-K. and Wang M.-J.J. (1995) "Image Thresholding by
68  // Minimizing the Measures of Fuzziness" Pattern Recognition, 28(1): 41-51
69  // Reimplemented (to handle 16-bit efficiently) by Johannes Schindelin Jan
70  // 31, 2011
71 
72  // Find first and last non-empty bin
73  size_t first, last;
74  for (first = 0; first < (size_t)hist.getSize() && hist[(unsigned char)first] == 0; first++) {
75  // do nothing
76  }
77 
78  for (last = (size_t)hist.getSize() - 1; last > first && hist[(unsigned char)last] == 0; last--) {
79  // do nothing
80  }
81 
82  if (first == last) {
83  return 0;
84  }
85 
86  // Calculate the cumulative density and the weighted cumulative density
87  std::vector<float> S(last + 1);
88  std::vector<float> W(last + 1);
89 
90  S[0] = (float)hist[0];
91  W[0] = 0.0f;
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];
95  }
96 
97  // Precalculate the summands of the entropy given the absolute difference x
98  // - mu (integral)
99  float C = (float)(last - first);
100  std::vector<float> Smu(last + 1 - first);
101 
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);
105  }
106 
107  // Calculate the threshold
108  int bestThreshold = 0;
109  float bestEntropy = std::numeric_limits<float>::max();
110 
111  for (size_t threshold = first; threshold <= last; threshold++) {
112  float entropy = 0;
113  int mu = vpMath::round(W[threshold] / S[threshold]);
114  for (size_t i = first; i <= threshold; i++) {
115  entropy += Smu[(size_t)std::abs((int)i - mu)] * hist[(unsigned char)i];
116  }
117 
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];
121  }
122 
123  if (bestEntropy > entropy) {
124  bestEntropy = entropy;
125  bestThreshold = (int)threshold;
126  }
127  }
128 
129  return bestThreshold;
130 }
131 
132 int computeThresholdIntermodes(const vpHistogram &hist)
133 {
134  if (hist.getSize() < 3) {
135  return -1;
136  }
137 
138  // Code based on the AutoThreshold ImageJ plugin:
139  // J. M. S. Prewitt and M. L. Mendelsohn, "The analysis of cell images," in
140  // Annals of the New York Academy of Sciences, vol. 128, pp. 1035-1053,
141  // 1966. ported to ImageJ plugin by G.Landini from Antti Niemisto's Matlab
142  // code (GPL) Original Matlab code Copyright (C) 2004 Antti Niemisto See
143  // http://www.cs.tut.fi/~ant/histthresh/ for an excellent slide presentation
144  // and the original Matlab code.
145  //
146  // Assumes a bimodal histogram. The histogram needs is smoothed (using a
147  // running average of size 3, iteratively) until there are only two local
148  // maxima. j and k Threshold t is (j+k)/2. Images with histograms having
149  // extremely unequal peaks or a broad and flat valley are unsuitable for this
150  // method.
151 
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];
155  }
156 
157  int iter = 0;
158  while (!isBimodal(hist_float)) {
159  // Smooth with a 3 point running mean filter
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;
162  }
163 
164  // First value
165  hist_float[0] = (hist_float[0] + hist_float[1]) / 2.0f;
166 
167  // Last value
168  hist_float[hist_float.size() - 1] = (hist_float.size() - 2 + hist_float.size() - 1) / 2.0f;
169 
170  iter++;
171 
172  if (iter > 10000) {
173  std::cerr << "Intermodes Threshold not found after 10000 iterations!" << std::endl;
174  return -1;
175  }
176  }
177 
178  // The threshold is the mean between the two peaks.
179  int tt = 0;
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]) {
182  // Mode
183  tt += (int)cpt;
184  }
185  }
186 
187  return (int)std::floor(tt / 2.0); // vpMath::round(tt / 2.0);
188 }
189 
190 int computeThresholdIsoData(const vpHistogram &hist, unsigned int imageSize)
191 {
192  int threshold = 0;
193 
194  // Code based on BSD Matlab isodata implementation by zephyr
195  // STEP 1: Compute mean intensity of image from histogram, set T=mean(I)
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];
202  }
203 
204  int T = vpMath::round(sum_ip[255] / imageSize);
205 
206  // STEP 2: compute Mean above T (MAT) and Mean below T (MBT) using T from
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)]);
209 
210  int T2 = vpMath::round((MAT + MBT) / 2.0f);
211 
212  //% STEP 3 to n: repeat step 2 if T(i)~=T(i-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)]);
216 
217  T = T2;
218  T2 = vpMath::round((MAT + MBT) / 2.0f);
219  threshold = T2;
220  }
221 
222  return threshold;
223 }
224 
225 int computeThresholdMean(const vpHistogram &hist, unsigned int imageSize)
226 {
227  // C. A. Glasbey, "An analysis of histogram-based thresholding algorithms,"
228  // CVGIP: Graphical Models and Image Processing, vol. 55, pp. 532-537, 1993.
229  // The threshold is the mean of the greyscale data
230  float sum_ip = 0.0f;
231  for (unsigned int cpt = 0; cpt < hist.getSize(); cpt++) {
232  sum_ip += cpt * (float)hist[cpt];
233  }
234 
235  return (int)std::floor(sum_ip / imageSize);
236 }
237 
238 int computeThresholdOtsu(const vpHistogram &hist, unsigned int imageSize)
239 {
240  // Otsu, N (1979), "A threshold selection method from gray-level
241  // histograms", IEEE Trans. Sys., Man., Cyber. 9: 62-66,
242  // doi:10.1109/TSMC.1979.4310076
243 
244  float mu_T = 0.0f;
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;
249  }
250 
251  // Weight Background / Foreground
252  float w_B = 0.0f, w_F = 0.0f;
253 
254  float max_sigma_b = 0.0f;
255  int threshold = 0;
256 
257  for (int cpt = 0; cpt < 256; cpt++) {
258  w_B += hist[cpt];
259  if (vpMath::nul(w_B, std::numeric_limits<float>::epsilon())) {
260  continue;
261  }
262 
263  w_F = (int)imageSize - w_B;
264  if (vpMath::nul(w_F, std::numeric_limits<float>::epsilon())) {
265  break;
266  }
267 
268  // Mean Background / Foreground
269  float mu_B = sum_ip_all[cpt] / (float)w_B;
270  float mu_F = (mu_T - sum_ip_all[cpt]) / (float)w_F;
271  // If there is a case where (w_B * w_F) exceed FLT_MAX, normalize
272  // histogram
273  float sigma_b_sqr = w_B * w_F * (mu_B - mu_F) * (mu_B - mu_F);
274 
275  if (sigma_b_sqr >= max_sigma_b) {
276  threshold = cpt;
277  max_sigma_b = sigma_b_sqr;
278  }
279  }
280 
281  return threshold;
282 }
283 
284 int computeThresholdTriangle(vpHistogram &hist)
285 {
286  int threshold = 0;
287 
288  // Zack, G. W., Rogers, W. E. and Latt, S. A., 1977,
289  // Automatic Measurement of Sister Chromatid Exchange Frequency,
290  // Journal of Histochemistry and Cytochemistry 25 (7), pp. 741-753
291 
292  int left_bound = -1, right_bound = -1, max_idx = -1, max_value = 0;
293  // Find max value index and left / right most index
294  for (int cpt = 0; cpt < (int)hist.getSize(); cpt++) {
295  if (left_bound == -1 && hist[cpt] > 0) {
296  left_bound = (int)cpt;
297  }
298 
299  if (right_bound == -1 && hist[(int)hist.getSize() - 1 - cpt] > 0) {
300  right_bound = (int)hist.getSize() - 1 - cpt;
301  }
302 
303  if ((int)hist[cpt] > max_value) {
304  max_value = (int)hist[cpt];
305  max_idx = cpt;
306  }
307  }
308 
309  // First / last index when hist(cpt) == 0
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;
312 
313  // Use the largest bound
314  bool flip = false;
315  if (max_idx - left_bound < right_bound - max_idx) {
316  // Flip histogram to get the largest bound to the left
317  flip = true;
318 
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);
324  }
325 
326  left_bound = (int)hist.getSize() - 1 - right_bound;
327  max_idx = (int)hist.getSize() - 1 - max_idx;
328  }
329 
330  // Distance from a point to a line defined by two points:
331  //\textbf{distance} \left ( P_1, P_2, \left ( x_0,y_0 \right ) \right )
332  // = \frac{ \left | \left ( y_2-y_1 \right ) x_0 - \left ( x_2-x_1 \right )
333  // y_0 + x_2 y_1 - y_2 x_1 \right | }
334  // { \sqrt{ \left ( y_2 - y_1 \right )^{2} + \left ( x_2 - x_1 \right
335  // )^{2} } }
336  // Constants are ignored
337  float a = (float)max_value; // y_2 - y_1
338  float b = (float)(left_bound - max_idx); //-(x_2 - x_1)
339  float max_dist = 0.0f;
340 
341  for (int cpt = left_bound + 1; cpt <= max_idx; cpt++) {
342  float dist = a * cpt + b * hist[cpt];
343 
344  if (dist > max_dist) {
345  max_dist = dist;
346  threshold = cpt;
347  }
348  }
349  threshold--;
350 
351  if (flip) {
352  threshold = (int)hist.getSize() - 1 - threshold;
353  }
354 
355  return threshold;
356 }
357 } // namespace
358 
359 namespace vp
360 {
362  const unsigned char backgroundValue, const unsigned char foregroundValue)
363 {
364  if (I.getSize() == 0) {
365  return 0;
366  }
367 
368  // Compute image histogram
369  vpHistogram histogram(I);
370  int threshold = -1;
371 
372  switch (method) {
374  threshold = computeThresholdHuang(histogram);
375  break;
376 
378  threshold = computeThresholdIntermodes(histogram);
379  break;
380 
382  threshold = computeThresholdIsoData(histogram, I.getSize());
383  break;
384 
385  case AUTO_THRESHOLD_MEAN:
386  threshold = computeThresholdMean(histogram, I.getSize());
387  break;
388 
389  case AUTO_THRESHOLD_OTSU:
390  threshold = computeThresholdOtsu(histogram, I.getSize());
391  break;
392 
394  threshold = computeThresholdTriangle(histogram);
395  break;
396 
397  default:
398  break;
399  }
400 
401  if (threshold != -1) {
402  // Threshold
403  vpImageTools::binarise(I, (unsigned char)threshold, (unsigned char)255, backgroundValue, foregroundValue,
404  foregroundValue);
405  }
406 
407  return threshold;
408 }
409 };
Class to compute a gray level image histogram.
Definition: vpHistogram.h:108
unsigned getSize() const
Definition: vpHistogram.h:262
void set(const unsigned char level, unsigned int value)
Definition: vpHistogram.h:225
static void binarise(vpImage< Type > &I, Type threshold1, Type threshold2, Type value1, Type value2, Type value3, bool useLUT=true)
Definition: vpImageTools.h:455
unsigned int getSize() const
Definition: vpImage.h:223
static bool nul(double x, double threshold=0.001)
Definition: vpMath.h:360
static int round(double x)
Definition: vpMath.h:323
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:66
@ AUTO_THRESHOLD_ISODATA
Definition: vpImgproc.h:75
@ AUTO_THRESHOLD_HUANG
Definition: vpImgproc.h:67
@ AUTO_THRESHOLD_INTERMODES
Definition: vpImgproc.h:71
@ AUTO_THRESHOLD_TRIANGLE
Definition: vpImgproc.h:87
@ AUTO_THRESHOLD_MEAN
Definition: vpImgproc.h:79
@ AUTO_THRESHOLD_OTSU
Definition: vpImgproc.h:83