最大类间方差（大津法）详解

参考： https://www.cnblogs.com/xiaomanon/p/4110006.html

一、算法介绍

• 最大类间方差法是1979年由日本学者大津提出的，是一种自适应阈值确定的方法，又叫大津法，简称OTSU，是一种基于全局的二值化算法。因为按照大津法求得的阈值进行图像二值化分割后，前景与背景图像的类间方差最大。

• 它被认为是图像分割中阈值选取的最佳算法，计算简单，不受图像亮度和对比度的影响，因此在数字图像处理上得到了广泛的应用。它是按图像的灰度特性，将图像分成背景和前景两部分。在OTSU算法中所采用的衡量差别的标准就是较为常见的最大类间方差。因方差是灰度分布均匀性的一种度量，背景和前景之间的类间方差越大，说明构成图像的两部分的差别越大。当部分前景错分为背景或部分背景错分为前景都会导致两部分差别变小。因此，使类间方差最大的分割意味着错分概率最小。

1. 应用：是求图像全局阈值的最佳方法，应用不言而喻，适用于大部分需要求图像全局阈值的场合。

2. 优点：计算简单快速，不受图像亮度和对比度的影响。

3. 缺点：对图像噪声敏感；只能针对单一目标分割；当目标和背景大小比例悬殊、类间方差函数可能呈现双峰或者多峰，这个时候效果不好。

二、算法思想

• 假设一幅图像有L个灰度级[1,2,…,L]。灰度级为i的像素点的个数为ni，那么总的像素点个数就应该为N=n1+n2+…+nL。为了讨论方便，我们使用归一化的灰度级直方图并且视为这幅图像的概率分布：
  
• 现在假设我们通过一个灰度级为k的门限将这些像素点划分为两类：C0和C1(背景和目标，或者反之亦然)；C0表示灰度级为[1,…,k]的像素点，C1表示灰度级为[k+1,…,L]的像素点。那么，每一类出现的概率以及各类的平均灰度级分别由下面的式子给出：
  
  
  其中：
  
  分别为灰度级从1到k的累积出现概率和平均灰度级（一阶累积矩）；
  而：
  
  是整幅图像的平均灰度级。
• 我们可以很容易验证，对于任意选定的k，都有：
  
  这两类的类内方差由下面的公式给出：
  
  这需要二阶累积矩(second-order cumulative moment，统计学概念)。
• 为了评价(灰度级k)这个门限“好”的程度，我们需要引入判别式分析中使用的判别式标准来测量（类的分离性测量）：
  
  其中：
  
  又根据式(9)，可以得出：
  
  这三个式子分别是类内方差、类间方差和灰度级的总方差。
• 然后，我们的问题就简化为一个优化问题，即寻找一个门限k使(12)式中给出的一个目标函数取最大值。
• 这个观点是出于这样一个猜想，一个好的阈值将会把灰度级分为两类，那么反过来说，就是如果一个门限能够在灰度级上将图像分割为最好的两类的话，那么这个门限就是最好的门限。
• 上面给出的判别式标准是分别求取λ、κ和η的最大值。然而，对于κ而言，它又等于另外一个，比如κ=λ+1；而对于λ而言，又有η=λ/(λ+1)，因为始终存在下面的基本关系：
  
  
  

三、效果展示

四、基于OpenCV实现

#include <iostream>
#include <opencv2/core.hpp>
#include <opencv2/highgui.hpp>
#include <opencv2/imgproc.hpp>
 
int Otsu(cv::Mat& src, cv::Mat& dst, int thresh){
	const int Grayscale = 256;
	int graynum[Grayscale] = { 0 };
	int r = src.rows;
	int c = src.cols;
	for (int i = 0; i < r; ++i){
		const uchar* ptr = src.ptr<uchar>(i);
		for (int j = 0; j < c; ++j){        //直方图统计
			graynum[ptr[j]]++;
		}
	}
 
    double P[Grayscale] = { 0 };   
	double PK[Grayscale] = { 0 };
	double MK[Grayscale] = { 0 };
	double srcpixnum = r*c, sumtmpPK = 0, sumtmpMK = 0;
	for (int i = 0; i < Grayscale; ++i){
		P[i] = graynum[i] / srcpixnum;   //每个灰度级出现的概率
		PK[i] = sumtmpPK + P[i];         //概率累计和 
		sumtmpPK = PK[i];
		MK[i] = sumtmpMK + i*P[i];       //灰度级的累加均值                                                                                                                                                                                                                                                                                                                                                                                                        
		sumtmpMK = MK[i];
	}
	
	//计算类间方差
	double Var=0;
	for (int k = 0; k < Grayscale; ++k){
		if ((MK[Grayscale-1] * PK[k] - MK[k])*(MK[Grayscale-1] * PK[k] - MK[k]) / (PK[k] * (1 - PK[k])) > Var){
			Var = (MK[Grayscale-1] * PK[k] - MK[k])*(MK[Grayscale-1] * PK[k] - MK[k]) / (PK[k] * (1 - PK[k]));
			thresh = k;
		}
	}
 
	//阈值处理
	src.copyTo(dst);
	for (int i = 0; i < r; ++i){
	    uchar* ptr = dst.ptr<uchar>(i);
		for (int j = 0; j < c; ++j){
			if (ptr[j]> thresh)
				ptr[j] = 255;
			else
				ptr[j] = 0;
		}
	}
	return thresh;
}
 
 
int main(){
	cv::Mat src = cv::imread("I:\\Learning-and-Practice\\2019Change\\Image process algorithm\\Img\\Fig1039(a)(polymersomes).tif");
	if (src.empty()){
		return -1;
	}
	if (src.channels() > 1)
		cv::cvtColor(src, src, CV_RGB2GRAY);
 
	cv::Mat dst,dst2;
	int thresh=0;
	double t2 = (double)cv::getTickCount();
	thresh=Otsu(src , dst, thresh); //Otsu
	std::cout << "Mythresh=" << thresh << std::endl;
	t2 = (double)cv::getTickCount() - t2;
	double time2 = (t2 *1000.) / ((double)cv::getTickFrequency());
	std::cout << "my_process=" << time2 << " ms. " << std::endl << std::endl;
    double  Otsu = 0;
 
	Otsu=cv::threshold(src, dst2, Otsu, 255, CV_THRESH_OTSU + CV_THRESH_BINARY);
	std::cout << "OpenCVthresh=" << Otsu << std::endl;
 
	cv::namedWindow("src", CV_WINDOW_NORMAL);
	cv::imshow("src", src);
	cv::namedWindow("dst", CV_WINDOW_NORMAL);
	cv::imshow("dst", dst);
	cv::namedWindow("dst2", CV_WINDOW_NORMAL);
	cv::imshow("dst2", dst2);
	//cv::imwrite("I:\\Learning-and-Practice\\2019Change\\Image process algorithm\\Image Filtering\\MeanFilter\\TXT.jpg",dst);
	cv::waitKey(0);
}

五、算法评价

• 就最大类间方差算法而言，优点是算法简单，当目标与背景的面积相差不大时，能够有效地对图像进行分割。但是，当图像中的目标与背景的面积相差很大时，表现为直方图没有明显的双峰，或者两个峰的大小相差很大，分割效果不佳，或者目标与背景的灰度有较大的重叠时也不能准确的将目标与背景分开。导致这种现象出现的原因是该方法忽略了图像的空间信息，同时该方法将图像的灰度分布作为分割图像的依据，因而对噪声也相当敏感。所以，在实际应用中，总是将其与其他方法结合起来使用。