大津阈值法原理_ostu阈值分割

//返回阈值的大津阈值法
double Otsu_threshold(const cv::Mat& InputImage)
{ 
   	
	cv::Mat SrcImage = InputImage.clone();
	CV_Assert(SrcImage.type() == CV_8UC1);

	int rows = SrcImage.rows;
	int cols = SrcImage.cols;

	const int L = 256;
	int N = rows * cols;					//灰度图大小(元素个数)
	int n_i[L] = { 
    0 };						//灰度直方图
	
	for (int i = 0; i < rows; ++i)
	{ 
   
		uchar* p = SrcImage.ptr<uchar>(i);
		for (int j = 0; j < cols; ++j)
			n_i[p[j]]++;
	}

	double pn_i[L];							//概率灰度直方图
	for (int i = 0; i < L; ++i)
	{ 
   
		pn_i[i] = (double)n_i[i] / N;
	}

	//全局均值和全局方差
	cv::Mat mat_mean, mat_stddev;
	double gray_mean, gray_sigma;

	cv::meanStdDev(SrcImage, mat_mean, mat_stddev);
	gray_mean = mat_mean.at<double>(0, 0);			//m_G
	//全局方差是用来计算类的可分离测度
	gray_sigma = mat_stddev.at<double>(0, 0) * mat_stddev.at<double>(0, 0);

	//遍历所有灰度级，计算类间方差
	std::vector<double>sigma_ks(L);
	for (int k = 0; k < L; ++k)
	{ 
   
		double p1 = 0.0;				//p1类发生概率
		double m_k = 0.0;				//高达k阶累计平均灰度

		for (int i = 0; i <= k; ++i)
		{ 
   
			p1 += pn_i[i];
			m_k += i * pn_i[i];
		}

		if (p1 == 0.0 || (1 - p1) == 0.0)					//分母不能为0
			sigma_ks[k] = 0.0;
		else
			sigma_ks[k] = (gray_mean * p1 - m_k) * (gray_mean * p1 - m_k) / (p1 * (1 - p1));
	}

	double max_Sigma_k = 0.0;
	std::vector<int>maxval_Ts;
	double Threshold_T = 0;									//最终输出的阈值T
	//找类间方差最大值
	for (int i = 0; i < sigma_ks.size(); ++i)
	{ 
   
		if (sigma_ks[i] > max_Sigma_k)
			max_Sigma_k = sigma_ks[i];
	}
	//找极大值对应的所有灰度值
	for (int i = 0; i < sigma_ks.size(); ++i)
	{ 
   
		if (abs(max_Sigma_k - sigma_ks[i]) < 1e-8)
			maxval_Ts.push_back(i);
	}
	//如果极大值点不唯一，那么取对应各个极大值的各个k的平均值来得到最终阈值threshold_T
	for (int i = 0; i < maxval_Ts.size(); ++i)
		Threshold_T += maxval_Ts[i];

	return Threshold_T / maxval_Ts.size();
}

//-----------------test--------------------//
int main()
{ 
   
	std::string path = "F:\\NoteImage\\Lena.jpg";
	cv::Mat src = imread(path, cv::IMREAD_GRAYSCALE);
	if (!src.data) { 
   
		std::cout << "Could not open or find the image" << std::endl;
		return -1;
	}

	cv::Mat dst;
	//对比一下opencv官方计算结果(显然结果是相同的)
	double thres1 = cv::threshold(src, dst, 0, 255, cv::THRESH_OTSU);
	double thres2 = Otsu_threshold(src);

	std::cout << "opencv = " << thres1 << " my = " << thres2;

	cv::waitKey(0);
	return 0;
}