详细及易读懂的 大津法（OTSU）原理 和 比opencv自带更快的算法实现

OTSU算法原理简述：

最大类间方差是由日本学者大津(Nobuyuki Otsu)于1979年提出，是一种自适应的阈值确定方法。算法假设图像像素能够根据阈值，被分成背景[background]和目标[objects]两部分。然后，计算该最佳阈值来区分这两类像素，使得两类像素区分度最大。

公式：  记 M = 256 单通道灰度分级 Sum = 像素总数

1. 背景像素占比   
2. 前景像素占比
3. 背景的平均灰度值
4. 前景的平均灰度值
5. 0~M灰度区间的灰度累计值
6. 类间方差：
7. 将公式3.4.5带入公式6 可得最终简化公式： 

 3、4的意思：每个灰度（0~255）比例 * 当前灰度的和。

代码：

#include <opencv2\opencv.hpp>
#include <iostream>
#include <time.h>
using namespace std;
using namespace cv;

int myOtsu(Mat & src, Mat &dst)
{
	
	int th = 0;
	const int GrayScale = 256;	//单通道图像总灰度256级
	int pixCount[GrayScale] = { 0 };//每个灰度值所占像素个数
	int pixSum = src.cols * src.rows;//图像总像素点
	float pixPro[GrayScale] = { 0 };//每个灰度值所占总像素比例
	float SumpixPro[GrayScale] = { 0 }; // 比例的和
	float WpixPro[GrayScale] = { 0 }; //比例 * 权重
	float SumWpixPro[GrayScale] = { 0 };//比例 * 权重 的 和
	float w0, w1, u0tmp, u1tmp, u0, u1, deltaTmp, deltaMax = 0;
	double start = getTickCount(); //开始时间
	for (int i = 0; i < src.rows; i++)
	{
		for (int j = 0; j < src.cols; j++)
		{
			int index = i * src.cols + j;
			pixCount[src.data[index]]++;//统计每个灰度级中像素的个数  
		}
	}
	double c1 = getTickCount(); 
	//cout << "c1 >> " << (c1 - start) / getTickFrequency() << endl;//输出时间
	for (int i = 0; i < GrayScale; i++)
	{
		pixPro[i] = pixCount[i] * 1.0 / pixSum;//计算每个灰度级的像素数目占整幅图像的比例  
		WpixPro[i] = i * pixPro[i];
		if (i == 0)
		{
			SumWpixPro[i] += WpixPro[i];
			SumpixPro[i] += pixPro[i];
		}
		else
		{
			SumWpixPro[i] = WpixPro[i] + SumWpixPro[i - 1];
			SumpixPro[i] = pixPro[i] + SumpixPro[i - 1];
		}
	}

	for (int i = 0; i < GrayScale; i++)//遍历所有从0到255灰度级的阈值分割条件，测试哪一个的类间方差最大
	{
		w0 = w1 = u0tmp = u1tmp = u0 = u1 = deltaTmp = 0;
		/*for (int j = 0; j < GrayScale; j++)
		{
			if (j <= i)//背景
			{
				w0 += pixPro[j];
				u0tmp += WpixPro[j];
			}
			else//前景
			{
				w1 += pixPro[j];
				u1tmp += WpixPro[j];
			}
		}*/
		w0 = SumpixPro[i];
		w1 = 1 - w0;
		
		if (w0 == 0 || w1 == 0)
			continue;
		u0tmp = SumWpixPro[i];
		u1tmp = SumWpixPro[255] - SumWpixPro[i];
		
		u0 = u0tmp / w0;
		u1 = u1tmp / w1;
		deltaTmp = (float)(w0 *w1* pow((u0 - u1), 2)); //类间方差公式 g = w1 * w2 * (u1 - u2) ^ 2
		if (deltaTmp > deltaMax)
		{
			deltaMax = deltaTmp;
			th = i;
		}
	}
	double c2 = getTickCount();
	//cout << "c2 >> " << (c2 - c1) / getTickFrequency() << endl;//输出时间
	for (int i = 0; i < src.rows; i++)
	{
		for (int j = 0; j < src.cols; j++)
		{
			int index = i * src.cols + j;
			if (src.data[index] > th)
				dst.data[index] = 255;
			else
				dst.data[index] = 0;
		}
	}
	double c3 = getTickCount();
	//cout << "c3 >> " << (c3 - c2) / getTickFrequency() << endl;//输出时间
	return th;
}

int main()
{
	Mat src = imread("scene.jpg", IMREAD_GRAYSCALE);//单通道读取图像
												   /*my_dst: 自己实现的大津法 得到的处理图像
												   otsu_dst：opencv自带的大津法 得到的处理图像
												   sub：两个处理图像相差图
												   */
	Mat my_dst, otsu_dst, sub;
	/*my_th: 自己实现的大津法 得到的最大类件方差 即阈值
	th：opencv自带的大津法 得到的最大类件方差 即阈值
	*/
	int my_th, th;

	/*计算开销时间，对比两个算法效率*/
	int times = 100;
	
	double start = getTickCount(); //开始时间
	my_dst = Mat::zeros(src.size(), CV_8UC1);
	while(times--)
	{
		
		my_th = myOtsu(src, my_dst);
	}
	double end = getTickCount();   //结束时间
	cout << "myOtsu end algorithm >> " << (end - start) / getTickFrequency() << endl; //输出时间
	cout << "myOtsu threshold >> " << my_th << endl;

	start = getTickCount(); //开始时间
	times = 100;
	while(times--)
	{
		th = threshold(src, otsu_dst, 0, 255, THRESH_OTSU);
	}
	end = getTickCount();   //结束时间
	cout << "Otsu end algorithm >> " << (end - start) / getTickFrequency() << endl;//输出时间
	cout << "Otsu threshold >> " << th << endl;


	waitKey();
	system("pause");
	return 0;
}

原图：用的是1920*1080的图像，循环100次得到测试结果和图像。

结论：可以看到自己实现的otsu和opencv自带的自适应阈值算法效率上相差大约3倍左右。

分析：遍历图片没有加速，opencv内部有加速。

分析时间慢的原因：第一次遍历整张图片，统计每个灰度级像素个数，第二次遍历整张图片进行阈值操作。

似曾相识第一次遍历不就是灰度直方图吗，然后第二次遍历不就是LUT图像压缩。

改进后源码：

#if 1
#include <opencv2\opencv.hpp>
#include <iostream>
#include <time.h>
using namespace std;
using namespace cv;

void getHistogram(Mat &src, int *dst)
{
	Mat hist;
	int channels[1] = { 0 };
	int histSize[1] = { 256 };
	float hranges[2] = { 0, 256.0 };
	const float *ranges[1] = { hranges };
	calcHist(&src, 1, channels, Mat(), hist, 1, histSize, ranges);
	//cout << hist ;
	for (int i = 0; i < 256; i++)
	{ 
		float binVal = hist.at<float>(i);
		dst[i] = int(binVal);
	}
}


int myOtsu(Mat & src, Mat &dst)
{
	
	int th = 0;
	const int GrayScale = 256;	//单通道图像总灰度256级
	int pixCount[GrayScale] = { 0 };//每个灰度值所占像素个数
	int pixSum = src.cols * src.rows;//图像总像素点
	float pixPro[GrayScale] = { 0 };//每个灰度值所占总像素比例
	float SumpixPro[GrayScale] = { 0 }; // 比例的和
	float WpixPro[GrayScale] = { 0 }; //比例 * 权重
	float SumWpixPro[GrayScale] = { 0 };//比例 * 权重 的 和
	float w0, w1, u0tmp, u1tmp, u0, u1, deltaTmp, deltaMax = 0;
	double start = getTickCount(); //开始时间
	
	/*for (int i = 0; i < src.rows; i++)
	{
		for (int j = 0; j < src.cols; j++)
		{
			int index = i * src.cols + j;
			pixCount[src.data[index]]++;//统计每个灰度级中像素的个数  
		}
	}*/
	getHistogram(src, pixCount);

	double c1 = getTickCount(); 
	//cout << "c1 >> " << (c1 - start) / getTickFrequency() << endl;//输出时间
	for (int i = 0; i < GrayScale; i++)
	{
		pixPro[i] = pixCount[i] * 1.0 / pixSum;//计算每个灰度级的像素数目占整幅图像的比例  
		WpixPro[i] = i * pixPro[i];
		if (i == 0)
		{
			SumWpixPro[i] += WpixPro[i];
			SumpixPro[i] += pixPro[i];
		}
		else
		{
			SumWpixPro[i] = WpixPro[i] + SumWpixPro[i - 1];
			SumpixPro[i] = pixPro[i] + SumpixPro[i - 1];
		}
	}

	for (int i = 0; i < GrayScale; i++)//遍历所有从0到255灰度级的阈值分割条件，测试哪一个的类间方差最大
	{
		w0 = w1 = u0tmp = u1tmp = u0 = u1 = deltaTmp = 0;
		/*for (int j = 0; j < GrayScale; j++)
		{
			if (j <= i)//背景
			{
				w0 += pixPro[j];
				u0tmp += WpixPro[j];
			}
			else//前景
			{
				w1 += pixPro[j];
				u1tmp += WpixPro[j];
			}
		}*/
		w0 = SumpixPro[i];
		w1 = 1 - w0;
		
		if (w0 == 0 || w1 == 0)
			continue;
		u0tmp = SumWpixPro[i];
		u1tmp = SumWpixPro[255] - SumWpixPro[i];
		
		u0 = u0tmp / w0;
		u1 = u1tmp / w1;
		deltaTmp = (float)(w0 *w1* pow((u0 - u1), 2)); //类间方差公式 g = w1 * w2 * (u1 - u2) ^ 2
		if (deltaTmp > deltaMax)
		{
			deltaMax = deltaTmp;
			th = i;
		}
	}
	double c2 = getTickCount();
	//cout << "c2 >> " << (c2 - c1) / getTickFrequency() << endl;//输出时间
	uchar lutData[256];
	for (int i = 0; i < 256; i++)
	{ 
		if (i > th)
			lutData[i] = 255;
		else
			lutData[i] = 0;
	}
	Mat lut(1, 256, CV_8UC1, lutData);
	LUT(src, lut, dst);
	/*for (int i = 0; i < src.rows; i++)
	{
		for (int j = 0; j < src.cols; j++)
		{
			int index = i * src.cols + j;
			if (src.data[index] > th)
				dst.data[index] = 255;
			else
				dst.data[index] = 0;
		}
	}*/
	double c3 = getTickCount();
	//cout << "c3 >> " << (c3 - c2) / getTickFrequency() << endl;//输出时间
	return th;
}

int main()
{
	Mat src = imread("scene.jpg", IMREAD_GRAYSCALE);//单通道读取图像
												   /*my_dst: 自己实现的大津法 得到的处理图像
												   otsu_dst：opencv自带的大津法 得到的处理图像
												   sub：两个处理图像相差图
												   */
	//Mat src = Mat(Size(3000, 3000), CV_8UC1, Scalar(127));
	Mat my_dst, otsu_dst, sub;
	/*my_th: 自己实现的大津法 得到的最大类件方差 即阈值
	th：opencv自带的大津法 得到的最大类件方差 即阈值
	*/
	int my_th, th;

	/*计算开销时间，对比两个算法效率*/
	int times = 100;
	
	double start = getTickCount(); //开始时间
	my_dst = Mat::zeros(src.size(), CV_8UC1);
	while(times--)
	{
		
		my_th = myOtsu(src, my_dst);
	}
	double end = getTickCount();   //结束时间
	cout << "myOtsu end algorithm >> " << (end - start) / getTickFrequency() << endl; //输出时间
	cout << "myOtsu threshold >> " << my_th << endl;

	start = getTickCount(); //开始时间
	times = 100;
	while(times--)
	{
		th = threshold(src, otsu_dst, 0, 255, THRESH_OTSU);
	}
	end = getTickCount();   //结束时间
	cout << "Otsu end algorithm >> " << (end - start) / getTickFrequency() << endl;//输出时间
	cout << "Otsu threshold >> " << th << endl;


	waitKey();
	system("pause");
	return 0;
}
#endif

发现优化后的效率比opencv自带的还快。

最后：欢迎大家的批评，很高心与大家分享，谢谢大家。