44 #include <visp3/core/vpHistogram.h> 45 #include <visp3/core/vpImageTools.h> 46 #include <visp3/imgproc/vpImgproc.h> 50 bool isBimodal(
const std::vector<float> &hist_float)
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]) {
79 for (first = 0; first < (size_t)hist.
getSize() && hist[(
unsigned char)first] == 0; first++) {
83 for (last = (
size_t)hist.
getSize() - 1; last > first && hist[(
unsigned char)last] == 0; last--) {
92 std::vector<float> S(last + 1);
93 std::vector<float> W(last + 1);
95 S[0] = (float)hist[0];
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];
104 float C = (float)(last - first);
105 std::vector<float> Smu(last + 1 - first);
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);
113 int bestThreshold = 0;
114 float bestEntropy = std::numeric_limits<float>::max();
116 for (
size_t threshold = first; threshold <= last; threshold++) {
119 for (
size_t i = first; i <= threshold; i++) {
120 entropy += Smu[(size_t)std::abs((
int)i - mu)] * hist[(
unsigned char)i];
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];
128 if (bestEntropy > entropy) {
129 bestEntropy = entropy;
130 bestThreshold = (int)threshold;
134 return bestThreshold;
137 int computeThresholdIntermodes(
const vpHistogram &hist)
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];
163 while (!isBimodal(hist_float)) {
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;
170 hist_float[0] = (hist_float[0] + hist_float[1]) / 2.0f;
173 hist_float[hist_float.size() - 1] = (hist_float.size() - 2 + hist_float.size() - 1) / 2.0f;
178 std::cerr <<
"Intermodes Threshold not found after 10000 iterations!" << std::endl;
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]) {
192 return (
int)std::floor(tt / 2.0);
195 int computeThresholdIsoData(
const vpHistogram &hist,
const unsigned int imageSize)
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];
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)]);
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)]);
230 int computeThresholdMean(
const vpHistogram &hist,
const unsigned int imageSize)
236 for (
unsigned int cpt = 0; cpt < hist.
getSize(); cpt++) {
237 sum_ip += cpt * (float)hist[cpt];
240 return (
int)std::floor(sum_ip / imageSize);
243 int computeThresholdOtsu(
const vpHistogram &hist,
const unsigned int imageSize)
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;
257 float w_B = 0.0f, w_F = 0.0f;
259 float max_sigma_b = 0.0f;
262 for (
int cpt = 0; cpt < 256; cpt++) {
264 if (
vpMath::nul(w_B, std::numeric_limits<float>::epsilon())) {
268 w_F = (int)imageSize - w_B;
269 if (
vpMath::nul(w_F, std::numeric_limits<float>::epsilon())) {
274 float mu_B = sum_ip_all[cpt] / (float)w_B;
275 float mu_F = (mu_T - sum_ip_all[cpt]) / (
float)w_F;
278 float sigma_b_sqr = w_B * w_F * (mu_B - mu_F) * (mu_B - mu_F);
280 if (sigma_b_sqr >= max_sigma_b) {
282 max_sigma_b = sigma_b_sqr;
297 int left_bound = -1, right_bound = -1, max_idx = -1, max_value = 0;
299 for (
int cpt = 0; cpt < (int)hist.
getSize(); cpt++) {
300 if (left_bound == -1 && hist[cpt] > 0) {
301 left_bound = (int)cpt;
304 if (right_bound == -1 && hist[(
int)hist.
getSize() - 1 - cpt] > 0) {
305 right_bound = (int)hist.
getSize() - 1 - cpt;
308 if ((
int)hist[cpt] > max_value) {
309 max_value = (int)hist[cpt];
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;
320 if (max_idx - left_bound < right_bound - max_idx) {
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);
331 left_bound = (int)hist.
getSize() - 1 - right_bound;
332 max_idx = (int)hist.
getSize() - 1 - max_idx;
342 float a = (float)max_value;
343 float b = (float)(left_bound - max_idx);
344 float max_dist = 0.0f;
346 for (
int cpt = left_bound + 1; cpt <= max_idx; cpt++) {
347 float dist = a * cpt + b * hist[cpt];
349 if (dist > max_dist) {
357 threshold = (int)hist.
getSize() - 1 - threshold;
375 const unsigned char backgroundValue,
const unsigned char foregroundValue)
387 threshold = computeThresholdHuang(histogram);
391 threshold = computeThresholdIntermodes(histogram);
395 threshold = computeThresholdIsoData(histogram, I.
getSize());
399 threshold = computeThresholdMean(histogram, I.
getSize());
403 threshold = computeThresholdOtsu(histogram, I.
getSize());
407 threshold = computeThresholdTriangle(histogram);
414 if (threshold != -1) {
Class to compute a gray level image histogram.
Provides simple mathematics computation tools that are not available in the C mathematics library (ma...
static int round(const double x)
static bool nul(double x, double s=0.001)
unsigned int getSize() const
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)