Visual Servoing Platform version 3.6.0
Loading...
Searching...
No Matches
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
43namespace
44{
45bool 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
62int 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
132int 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
190int 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
225int 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
238int 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
284int 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
359namespace 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
386 threshold = computeThresholdMean(histogram, I.getSize());
387 break;
388
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.
unsigned getSize() const
void set(const unsigned char level, unsigned int value)
static void binarise(vpImage< Type > &I, Type threshold1, Type threshold2, Type value1, Type value2, Type value3, bool useLUT=true)
Definition of the vpImage class member functions.
Definition vpImage.h:135
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