otsu算法详解

参考：https://blog.csdn.net/leonardohaig/article/details/120269341

简介

OTSU算法称为最大类间方差法，也称为大津法，由大津再1979年提出，被认为是图像分割中阈值选取的最佳算法，在图像处理中得到广泛应用。它是根据图像的灰度特性，将图像分割成前景区域（object）和背景区域（background），其基本思想是根据分割的阈值计算一个类间方差，将类间方差最大时对应的阈值作为最佳阈值。

类间方差推导

设图像尺寸是MxN，其二值化的的阈值为k，该阈值将图像分成背景和前景两个区域，其中背景区域的像素数量是$N_0$，前景区域的像素数量为$N_1$。背景像素点的占比设为$w_0$，灰度均值为$u_0$，前景像素点的占比设为$w_1$，灰度均值为$u_1$。整幅图像的灰度均值为$u$，则：

具体推导请看参考

理想情况下，对于同一类，其类内方差是很小的，同时对于背景和前景的类间方差是很大的。因此只要遍历找到使得类间方差最大的阈值T即为所求。

算法步骤与实现

步骤：

1. 根据图像的尺寸确定图像总的像素点数量MxN
2. 遍历图像灰度级0-255中的每一个灰度值L，假设L将图像分为背景(Background)和目标(Object)两部分，然后计算类间方差
3. 遍历过程中使得类间方差最大的阈值T即为所求

代码：

#include <iostream>
#include <opencv2/opencv.hpp>
using namespace std;
using namespace cv;
 
/*
	功能：大津法阈值分割
	输入：gray img
	输出：thresh
*/
int thresholdOTSU(const Mat& img)
{
	CV_Assert(img.type() == CV_8UC1);
	int thresh = -1;
	
	int pixelNum = img.rows * img.cols;

	// 图像直方图统计
	int pixelCount[256] = { 0 };
	Size size(img.size());
	if (img.isContinuous()) {
		size.width *= size.height;
		size.height = 1;
	}
	for (int i = 0; i < size.height; i++) {
		const uchar* sdata = img.ptr(i);
		for (int j = 0; j < size.width; j++) {
			pixelCount[sdata[j]] += 1;
		}
	}

	// 穷举0-255的灰度阈值，找到使得类间方差最大的阈值k
	const float Ppixelnum = 1.0f / (img.rows * img.cols);
	float sigma_max = -FLT_MAX;
	for (int k = 0; k < 256; k++) {
		int backNum = 0, objectNum = 0;
		float backSum = 0, objectSum = 0;
		for (int i = 0; i <= k; i++) {
			backNum += pixelCount[i];
			backSum += i * pixelCount[i];
		}
		for (int i = k + 1; i < 256; i++) {
			objectNum += pixelCount[i];
			objectSum += i * pixelCount[i];
		}

		//求类间方差
		const float w0 = backNum * Ppixelnum;
		const float w1 = objectNum * Ppixelnum;
		const float u0 = backSum / backNum;
		const float u1 = objectSum / objectNum;
		float sigma = w0 * w1 * pow(u0 - u1, 2);
		if (sigma > sigma_max) {
			sigma_max = sigma;
			thresh = k;
		}
	}

	return thresh;
}

void ShowHist(const cv::Mat& image, cv::Mat& hist_img, int scale = 2)
{
	int bins = 256;
	int hist_size[] = { bins };
	float range[] = { 0, 255 };
	const float* ranges[] = { range };
	Mat hist;
	int channels[] = { 0 };

	cv::calcHist(&image, 1, channels, cv::Mat(), // do not use mask
		hist, 1, hist_size, ranges,
		true, // the histogram is uniform
		false);

	double max_val;
	cv::minMaxLoc(hist, nullptr, &max_val, nullptr, nullptr);
	const int hist_height = image.rows;
	hist_img = cv::Mat::zeros(hist_height, bins * scale, CV_8UC3);
	for (int i = 0; i < bins; i++)
	{
		float bin_val = hist.at<float>(i);
		int intensity = cvRound(bin_val * hist_height / max_val);  //要绘制的高度
		cv::rectangle(hist_img, cv::Point(i * scale, hist_height - 1),
			cv::Point((i + 1) * scale - 1, hist_height - intensity),
			CV_RGB(255, 255, 255));
	}
}

int main()
{
	Mat img = imread("test_ostu.bmp",IMREAD_GRAYSCALE);
	imshow("img", img);

	//求解直方图
	Mat hist_img;
	int scale = 2;
	ShowHist(img, hist_img, scale);

	//ostu二值化
	Mat ourdst, cvdst;

	auto thresh = thresholdOTSU(img);
	threshold(img, ourdst, thresh, 255, CV_THRESH_BINARY);

	//在直方图上绘制阈值位置
	line(hist_img, Point(thresh * scale, 0), Point(thresh * scale, img.rows - 1), Scalar(255, 0, 0), 2);

	auto thresh_cv = threshold(img, cvdst, 0, 255, THRESH_OTSU);
	cout << thresh << " " << thresh_cv << endl;
	
	imshow("hist", hist_img);
	imshow("ourdst", ourdst);
	imshow("cvdst", cvdst);

	waitKey(0);
	return 0;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118

结果请看参考

在参考的那篇博客中，除了从基本思想上对类间方差进行了推导之外，还从概率论角度说明了最大化类间方差等价于最小化类内方差，同时给出了另一条计算类间方差的公式。此外，还扩展到了多级大津法，也就是运用同样的思想，但可以将图像分割成多个类（大于两个）。

算法评价

• 应用：是求图像全局阈值的最佳方法。
• 优点：计算简单快速，不受图像亮度和对比度的影响。
• 缺点：对图像噪声敏感；当前景和背景大小比例（面积）悬殊时效果不好。
• 解释：当图像中的前景与背景的面积相差很大时，直方图没有明显的双峰，或者两个峰的大小相差很大，分割效果不佳，或者前景与背景的灰度有较大的重叠时也不能准确的将前景与背景分开。导致这种现象出现的原因是该方法忽略了图像的空间信息，同时该方法将图像的灰度分布作为分割图像的依据，因而对噪声也相当敏感。所以，在实际应用中，总是将其与其他方法结合起来使用，比如高斯滤波。