OTSU大津法图像分割自适应阈值c++实现（亲测可用）

内容参考原文《A Threshold Selection Method from Gray-Level Histograms》
引用自公众号： 写bug的程旭源

最大类间方差是由日本学者大津(Nobuyuki Otsu)于1979年提出，是一种自适应的阈值确定方法。算法假设图像像素能够根据阈值，被分成背景[background]和目标[objects]两部分。然后，计算该最佳阈值来区分这两类像素，使得两类像素区分度最大【用方差表达，具体公式见后】。OTSU的扩展算法，可进行多级阈值处理，称为“Multi Otsu method”【题外话】

设原始灰度级为M，灰度级为i的像素点个数为ni，对灰度直方图进行归一化：




opencv实现代码：

// An highlighted block
// m_otsu.cpp : 定义控制台应用程序的入口点。
//
#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;
	const char* Image1Name = "D:\\C++demo\\imgSeg\\img2.jpg";
	image = imread(Image1Name);
	cvtColor(image, image, CV_RGB2GRAY);
	//CV_WINDOW_NORMAL参数：表示用户可以改变窗口大小
	namedWindow("原图窗口", CV_WINDOW_NORMAL);
	imshow("原图窗口", image);
	
	//设置图像兴趣域，如果不需要可以删除这部分
	//设置ROI区域 从点（0，875）为左上角，1900宽，1800长
	Rect rect(0, 930, 1900,1800);
	//在原图上画出ROI
	rectangle(image, rect, Scalar(0, 255, 255), 2, 8);
	Mat imageROI = image(rect);
	namedWindow("imageROI", CV_WINDOW_NORMAL);
	imshow("imageROI", imageROI);

	Mat imageOutput;
	Mat imageOtsu;
	int thresholdValue = OtsuAlgThreshold(imageROI);
	cout << "类间方差为： " << thresholdValue << endl;
	//用thresholdValue 处理图片，阈值可以自己再调节，此处-4效果更好了
	threshold(image, imageOutput, thresholdValue - 4 , 255, CV_THRESH_BINARY);
	Opencv Otsu算法 作为对比
	threshold(image, imageOtsu, 0, 255, CV_THRESH_OTSU); 
	namedWindow("Output Image", CV_WINDOW_NORMAL);
	imshow("Output Image", imageOutput);
	namedWindow("Opencv Otsu", CV_WINDOW_NORMAL);
	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数组下标是灰度值，保存内容是灰度值对应像素点数
	//遍历 rows和cols
	for (int i = 0; i < image.rows; i++) 
	{
		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;
}