基于ImageJ和OpenCv的16种自动二值化图像代码(C++)

/** @file * @brief 自动计算单通道图像二值化阈值 * @author ImageJ fh * 由于VC编译器Max, Min函数可能会出现问题, 所以有关求最值函数全部自定义在cpp文件中 */ #include "auto_thresholder.h" /*****************************************常用函数*********************************************/ static int CalcMaxValue(int a, int b) { return (a > b) ? a : b; } static double CalcMaxValue(double a, double b) { return (a > b) ? a : b; } static int CalcMinValue(int a, int b) { return (a < b) ? a : b; } static double CalcMinValue(double a, double b) { return (a < b) ? a : b; } void AutoThresholder::CalcHist(const cv::Mat& src, std::vector<int>& data) { CV_Assert(CV_8UC1 == src.type()); data.clear(); // 初始化数组 const int GRAY_LEVEL = 256; for (int i = 0; i < GRAY_LEVEL; ++i) { data.push_back(0); } // 统计每个灰度级及其对应的像素个数 for (int y = 0; y < src.rows; ++y) { for (int x = 0; x < src.cols; ++x) { int value = src.at<uchar>(y, x); data.at(value) += 1; } } } int AutoThresholder::DefaultIsoData(const std::vector<int>& data) { // This is the modified IsoData method used by the "Threshold" widget in "Default" mode int n = data.size(); std::vector<int> data2(n, 0); int mode = 0, maxCount = 0; for (int i = 0; i < n; i++) { int count = data[i]; data2[i] = data[i]; if (data2[i] > maxCount) { maxCount = data2[i]; mode = i; } } int maxCount2 = 0; for (int i = 0; i < n; i++) { if ((data2[i] > maxCount2) && (i != mode)) maxCount2 = data2[i]; } int hmax = maxCount; if ((hmax > (maxCount2 * 2)) && (maxCount2 != 0)) { hmax = (int)(maxCount2 * 1.5); data2[mode] = hmax; } return IJIsoData(data2); } int AutoThresholder::IJIsoData(std::vector<int>& data) { // This is the original ImageJ IsoData implementation, here for backward compatibility. int level; int maxValue = data.size() - 1; double result, sum1, sum2, sum3, sum4; int count0 = data[0]; data[0] = 0; //set to zero so erased areas aren't included int countMax = data[maxValue]; data[maxValue] = 0; int min = 0; while ((data[min] == 0) && (min < maxValue)) min++; int max = maxValue; while ((data[max] == 0) && (max > 0)) max--; if (min >= max) { data[0] = count0; data[maxValue] = countMax; level = data.size() / 2; return level; } int movingIndex = min; int inc = CalcMaxValue(max / 40, 1); do { sum1 = sum2 = sum3 = sum4 = 0.0; for (int i = min; i <= movingIndex; i++) { sum1 += (double)i*data[i]; sum2 += data[i]; } for (int i = (movingIndex + 1); i <= max; i++) { sum3 += (double)i*data[i]; sum4 += data[i]; } result = (sum1 / sum2 + sum3 / sum4) / 2.0; movingIndex++; } while ((movingIndex + 1) <= result && movingIndex < max - 1); data[0] = count0; data[maxValue] = countMax; level = (int)std::round(result); return level; } int AutoThresholder::Huang(const std::vector<int>& data) { const int GRAY_LEVELS = 256; int threshold = -1; int ih = 0; int it = 0; int first_bin = 0; int last_bin = 0; double sum_pix = 0.0; double num_pix = 0.0; double term = 0.0; double ent = 0.0; // entropy double min_ent = 0.0; // min entropy double mu_x = 0.0; /* Determine the first non-zero bin */ first_bin = 0; for (ih = 0; ih < GRAY_LEVELS; ih++) { if (data[ih] != 0) { first_bin = ih; break; } } /* Determine the last non-zero bin */ last_bin = 255; for (ih = 255; ih >= first_bin; ih--) { if (data[ih] != 0) { last_bin = ih; break; } } term = 1.0 / (double)(last_bin - first_bin); std::vector<double> mu_0(GRAY_LEVELS, 0.0); sum_pix = num_pix = 0; for (ih = first_bin; ih < GRAY_LEVELS; ih++) { sum_pix += (double)ih * data[ih]; num_pix += data[ih]; /* NUM_PIX cannot be zero ! */ mu_0[ih] = sum_pix / num_pix; } //double[] mu_1 = new double[256]; std::vector<double> mu_1(GRAY_LEVELS, 0.0); sum_pix = num_pix = 0; for (ih = last_bin; ih > 0; ih--) { sum_pix += (double)ih * data[ih]; num_pix += data[ih]; /* NUM_PIX cannot be zero ! */ mu_1[ih - 1] = sum_pix / (double)num_pix; } /* Determine the threshold that minimizes the fuzzy entropy */ threshold = -1; min_ent = DBL_MAX; for (it = 0; it < GRAY_LEVELS; it++) { ent = 0.0; for (ih = 0; ih <= it; ih++) { /* Equation (4) in Ref. 1 */ mu_x = 1.0 / (1.0 + term * std::abs(ih - mu_0[it])); if (!((mu_x < 1e-06) || (mu_x > 0.999999))) { /* Equation (6) & (8) in Ref. 1 */ ent += data[ih] * (-mu_x * std::log(mu_x) - (1.0 - mu_x) * std::log(1.0 - mu_x)); } } for (ih = it + 1; ih < GRAY_LEVELS; ih++) { /* Equation (4) in Ref. 1 */ mu_x = 1.0 / (1.0 + term * std::abs(ih - mu_1[it])); if (!((mu_x < 1e-06) || (mu_x > 0.999999))) { /* Equation (6) & (8) in Ref. 1 */ ent += data[ih] * (-mu_x * std::log(mu_x) - (1.0 - mu_x) * std::log(1.0 - mu_x)); } } /* No need to divide by NUM_ROWS * NUM_COLS * LOG(2) ! */ if (ent < min_ent) { min_ent = ent; threshold = it; } } return threshold; } static bool BimodalTest(std::vector<double>& y) { int len = y.size(); bool b = false; int modes = 0; for (int k = 1; k < len - 1; k++) { if (y[k - 1] < y[k] && y[k + 1] < y[k]) { modes++; if (modes > 2) return false; } } if (modes == 2) b = true; return b; } int AutoThresholder::Intermodes(const std::vector<int>& data) { int maxbin = -1; for (int i = 0; i < data.size(); i++) { if (data[i] > 0) maxbin = i; } int minbin = -1; for (int i = data.size() - 1; i >= 0; i--) { if (data[i] > 0) minbin = i; } int length = (maxbin - minbin) + 1; std::vector<double> hist(data.size(), 0); for (int i = minbin; i <= maxbin; i++) { hist[i - minbin] = data[i]; } int iter = 0; int threshold = -1; while (!BimodalTest(hist)) { //smooth with a 3 point running mean filter double previous = 0; double current = 0; double next = hist[0]; for (int i = 0; i < length - 1; i++) { previous = current; current = next; next = hist[i + 1]; hist[i] = (previous + current + next) / 3; } hist[length - 1] = (current + next) / 3; iter++; if (iter > 10000) { threshold = -1; return threshold; } } // The threshold is the mean between the two peaks. int tt = 0; for (int i = 1; i < length - 1; i++) { if (hist[i - 1] < hist[i] && hist[i + 1] < hist[i]) { tt += i; } } threshold = (int)std::floor(tt / 2.0); return threshold + minbin; } int AutoThresholder::IsoData(const std::vector<int>& data) { int i = 0; int l = 0; int totl = 0; int g = 0; double toth = 0.0; double h = 0.0; const int GRAY_LEVELS = 256; for (i = 1; i < GRAY_LEVELS; i++) { if (data[i] > 0) { g = i + 1; break; } } while (true) { l = 0; totl = 0; for (i = 0; i < g; i++) { totl = totl + data[i]; l = l + (data[i] * i); } h = 0; toth = 0; for (i = g + 1; i < GRAY_LEVELS; i++) { toth += data[i]; h += ((double)data[i] * i); } if (totl > 0 && toth > 0) { l /= totl; h /= toth; if (g == (int)std::round((l + h) / 2.0)) break; } g++; if (g > 254) return -1; } return g; } int AutoThresholder::Li(const std::vector<int>& data) { int threshold = 0; double num_pixels = 0; double sum_back = 0; /* sum of the background pixels at a given threshold */ double sum_obj = 0; /* sum of the object pixels at a given threshold */ double num_back = 0; /* number of background pixels at a given threshold */ double num_obj = 0; /* number of object pixels at a given threshold */ double old_thresh = 0; double new_thresh = 0; double mean_back = 0; /* mean of the background pixels at a given thre