自适应阈值分割—大津法（OTSU算法）C++实现

大津法是一种图像灰度自适应的阈值分割算法，是1979年由日本学者大津提出，并由他的名字命名的。大津法按照图像上灰度值的分布，将图像分成背景和前景两部分看待，前景就是我们要按照阈值分割出来的部分。背景和前景的分界值就是我们要求出的阈值。遍历不同的阈值，计算不同阈值下对应的背景和前景之间的类内方差，当类内方差取得极大值时，此时对应的阈值就是大津法（OTSU算法）所求的阈值。

何为类间方差？

对于图像I(x,y)，前景(即目标)和背景的分割阈值记作T，属于前景的像素点数占整幅图像的比例记为ω0，其平均灰度μ0；背景像素点数占整幅图像的比例为ω1，其平均灰度为μ1。图像的总平均灰度记为μ，类间方差记为g。

假设图像的背景较暗，并且图像的大小为M×N，图像中像素的灰度值小于阈值T的像素个数记作N0，像素灰度大于阈值T的像素个数记作N1，则有：

　　　　　　ω0=N0/ M×N (1)

　　　　　　ω1=N1/ M×N (2)

　　　　　　N0+N1=M×N (3)

　　　　　　ω0+ω1=1　　　 (4)

　　　　　　μ=ω0*μ0+ω1*μ1 (5)

　　　　　　g=ω0(μ0-μ)^2+ω1(μ1-μ)^2 (6)

将式(5)代入式(6),得到等价公式:

　　　　　　g=ω0ω1(μ0-μ1)^2 (7)　这个就是类间方差的公式表述

采用遍历的方法得到使类间方差g最大的阈值T,即为所求。

Otsu实现思路

1. 计算0~255各灰阶对应的像素个数，保存至一个数组中，该数组下标是灰度值，保存内容是当前灰度值对应像素数

2. 计算背景图像的平均灰度、背景图像像素数所占比例

3. 计算前景图像的平均灰度、前景图像像素数所占比例

4. 遍历0~255各灰阶，计算并寻找类间方差极大值

C++代码实现：

#include <opencv2/highgui/highgui.hpp>

#include <opencv2/imgproc/imgproc.hpp>

#include <opencv2/core/core.hpp>

#include <iostream>

using namespace cv;

using namespace std;

//***************Otsu算法通过求类间方差极大值求自适应阈值******************

int OtsuAlgThreshold(const Mat image);

int main(int argc,char *argv[])

{

	Mat image=imread(argv[1]);

	imshow("SoureImage",image);

	cvtColor(image,image,CV_RGB2GRAY);

	Mat imageOutput;

	Mat imageOtsu;

	int thresholdValue=OtsuAlgThreshold(image);

	cout<<"类间方差为： "<<thresholdValue<<endl;

	threshold(image,imageOutput,thresholdValue,255,CV_THRESH_BINARY);

	threshold(image,imageOtsu,0,255,CV_THRESH_OTSU); //Opencv Otsu算法

	//imshow("SoureImage",image);

	imshow("Output Image",imageOutput);

	imshow("Opencv Otsu",imageOtsu);

	waitKey();

	return 0;

}

int OtsuAlgThreshold(const Mat image)

{

	if(image.channels()!=1)

	{

		cout<<"Please input Gray-image!"<<endl;

		return 0;

	}

	int T=0; //Otsu算法阈值

	double varValue=0; //类间方差中间值保存

	double w0=0; //前景像素点数所占比例

	double w1=0; //背景像素点数所占比例

	double u0=0; //前景平均灰度

	double u1=0; //背景平均灰度

	double Histogram[256]={0}; //灰度直方图，下标是灰度值，保存内容是灰度值对应的像素点总数

	uchar *data=image.data;

	double totalNum=image.rows*image.cols; //像素总数

	//计算灰度直方图分布，Histogram数组下标是灰度值，保存内容是灰度值对应像素点数

	for(int i=0;i<image.rows;i++)   //为表述清晰，并没有把rows和cols单独提出来

	{

		for(int j=0;j<image.cols;j++)

		{

			Histogram[data[i*image.step+j]]++;

		}

	}

	for(int i=0;i<255;i++)

	{

		//每次遍历之前初始化各变量

		w1=0;		u1=0;		w0=0;		u0=0;

		//***********背景各分量值计算**************************

		for(int j=0;j<=i;j++) //背景部分各值计算

		{

			w1+=Histogram[j];  //背景部分像素点总数

			u1+=j*Histogram[j]; //背景部分像素总灰度和

		}

		if(w1==0) //背景部分像素点数为0时退出

		{

			break;

		}

		u1=u1/w1; //背景像素平均灰度

		w1=w1/totalNum; // 背景部分像素点数所占比例

		//***********背景各分量值计算**************************

		//***********前景各分量值计算**************************

		for(int k=i+1;k<255;k++)

		{

			w0+=Histogram[k];  //前景部分像素点总数

			u0+=k*Histogram[k]; //前景部分像素总灰度和

		}

		if(w0==0) //前景部分像素点数为0时退出

		{

			break;

		}

		u0=u0/w0; //前景像素平均灰度

		w0=w0/totalNum; // 前景部分像素点数所占比例

		//***********前景各分量值计算**************************

		//***********类间方差计算******************************

		double varValueI=w0*w1*(u1-u0)*(u1-u0); //当前类间方差计算

		if(varValue<varValueI)

		{

			varValue=varValueI;

			T=i;

		}

	}

	return T;

}

原图像：