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