中值阈值和otsu（大津法）的源码实现

中值阈值和大津算法的原理及手动实现。

OTSU

介绍

大津法（OTSU）又名最大类间差法，由日本学者大津于1979年提出。被认为是图像分割中阈值选取的最佳算法，计算简单，不受图像亮度和对比度的影响。

大津法是按图像的灰度特征，把图像分成前景和背景两部分。因方差是灰度分布均匀性的一种度量,背景和前景之间的类间方差越大,说明构成图像的两部分的差别越大,当部分前景错分为背景或部分背景错分为前景都会导致两部分差别变小。因此使用类间方差最大的分割意味着错分概率最小。

原理

对于图像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,即为所求。

代码

//使用大津法Mat的阈值
int CvUtils::GetMatOTSU(Mat& img)
{
  //判断如果不是单通道直接返回128
  if (img.channels() > 1) return 128;
  int rows = img.rows;
  int cols = img.cols;
  //定义数组
  float mathists[256] = { 0 };
  //遍历计算0-255的个数
  for (int row = 0; row < rows; ++row) {
    for (int col = 0; col < cols; ++col) {
      int val = img.at<uchar>(row, col);
      mathists[val]++;
    }
  }
 
  //定义灰度级像素在整个图像中的比例
  float grayPro[256] = { 0 };
  int matSize = rows * cols;
  for (int i = 0; i < 256; ++i) {
    grayPro[i] = (float)mathists[i] / (float)matSize;
  }
 
  //大津法OTSU，前景与背景分割，计算出方差最大的灰度值
  int calcval;
  int calcMax = 0;
  for (int i = 0; i < 256; ++i) {
    float w0 = 0, w1 = 0, u0tmp = 0, u1tmp = 0, u0 = 0, u1 = 0, u = 0, calctmp = 0;
      
    for (int k = 0; k < 256; ++k) {
      float curGray = grayPro[k];
      //计算背景部分
      if (k <= i) {
        //以i为阈值分类，第一类总的概率
        w0 += curGray;
        u0tmp += curGray * k;
      }
      //计算前景部分
      else {
        //以i为阈值分类，第一类总的概率
        w1 += curGray;
        u1tmp += curGray * k;
      }
    }
 
    //求出第一类和第二类的平均灰度
    u0 = u0tmp / w0;
    u1 = u1tmp / w1;
    //求出整幅图像的平均灰度
    u = u0tmp + u1tmp;
 
    //计算类间方差
    calctmp = w0 * pow((u0 - u), 2) + w1 * pow((u1 - u), 2);
 
    //更新最大类间方差，并设置阈值
    if (calctmp > calcMax) {
      calcMax = calctmp;
      calcval = i;
    }
  }
 
  return calcval;
}

中值阈值

要实现自动阈值，方法就是求出图像的灰度直方图，直方图中找出中位数，然后根据中位数值设定一个标准差值，用中位数的值加上标准差来求出高低阈值。

实现思路：

• 图像转灰度图
• 求出灰度直方图，并找到中位数
• 根据中位数和设定的sigma值求出高低阈值
• 使用Canny边缘检测

代码

//求Mat的中位数
int CvUtils::GetMatMidVal(Mat& img)
{
  //判断如果不是单通道直接返回128
  if (img.channels() > 1) return 128;
  int rows = img.rows;
  int cols = img.cols;
  //定义数组
  float mathists[256] = { 0 };
  //遍历计算0-255的个数
  for (int row = 0; row < rows; ++row) {
    for (int col = 0; col < cols; ++col) {
      int val = img.at<uchar>(row, col);
      mathists[val]++;
    }
  }

  int calcval = rows * cols / 2;
  int tmpsum = 0;
  for (int i = 0; i < 255; ++i) {
    tmpsum += mathists[i];
    if (tmpsum > calcval) {
      return i;
    }
  }
  return 0;
}

高低阈值

代码

通过设定默认sigma的值，求出高低阈值进行canny。

void GetMinMaxthreshold(Mat& img, int& minval, int& maxval,int type, float sigma)
{
	int midval;
	if (type == 0)
	{
		midval = getmatmidval(img);
	}
	else {
		midval = GetMatOTSU(img);
	}
	// 计算低阈值
	minval = saturate_cast<uchar>((1.0 - sigma) * midval);
	//计算高阈值
	maxval = saturate_cast<uchar>((1.0 + sigma) * midval);
}

测试代码

测试了中值阈值和大津法对图像进行canny边缘提取，还可以加入
threshold(gray, dst5, 0, 255, THRESH_BINARY+THRESH_OTSU);
测试调用opencv库实现的大津法和手动实现大津法的结果进行比对。

#include<opencv.hpp>
#include<iostream>

using namespace std;
using namespace cv;

int getmatmidval(Mat& img)
{
	if (img.channels() > 1) return 128;
	int rows = img.rows;
	int cols = img.cols;
	//定义数组
	float mathists[256] = { 0 };
	//遍历计算0-255的个数
	for (int row = 0; row < rows; ++row) {
		uchar* ptr = img.ptr<uchar>(row);
		for (int col = 0; col < cols; ++col) {
			int val = *ptr;
			mathists[val]++;
			ptr++;
		}
	}

	int calval = rows * cols / 2;
	int tempsum = 0;
	for (int i = 0; i < 255; ++i)
	{
		tempsum += mathists[i];
		if (tempsum > calval) return i;
	}
	return 0;
}
int GetMatOTSU(Mat& img)
{

	if (img.channels() > 1) return 128;
	int rows = img.rows;
	int cols = img.cols;
	//定义数组
	float mathists[256] = { 0 };
	//遍历计算0-255的个数
	for (int row = 0; row < rows; ++row) {
		for (int col = 0; col < cols; ++col) {
			int val = img.at<uchar>(row, col);
			mathists[val]++;
		}
	}
	//定义灰度级像素在整个图像中的比例
	float grayPro[256] = { 0 };
	int matSize = rows * cols;
	for (int i = 0; i < 256; ++i) {
		grayPro[i] = (float)mathists[i] / (float)matSize;
	}

	//大津法OTSU，前景与背景分割，计算出方差最大的灰度值
	int calcval;
	int calcMax = 0;
	for (int i = 0; i < 256; ++i) {
		float w0 = 0, w1 = 0, u0tmp = 0, u1tmp = 0, u0 = 0, u1 = 0, u = 0, calctmp = 0;

		for (int k = 0; k < 256; ++k) {
			float curGray = grayPro[k];
			//计算背景部分
			if (k <= i) {
				//以i为阈值分类，第一类总的概率
				w0 += curGray;
				u0tmp += curGray * k;
			}
			//计算前景部分
			else {
				//以i为阈值分类，第一类总的概率
				w1 += curGray;
				u1tmp += curGray * k;
			}
		}

		//求出第一类和第二类的平均灰度
		u0 = u0tmp / w0;
		u1 = u1tmp / w1;
		//求出整幅图像的平均灰度
		u = u0tmp + u1tmp;

		//计算类间方差
		calctmp = w0 * pow((u0 - u), 2) + w1 * pow((u1 - u), 2);

		//更新最大类间方差，并设置阈值
		if (calctmp > calcMax) {
			calcMax = calctmp;
			calcval = i;
		}
	}


	return calcval;
}
void GetMinMaxthreshold(Mat& img, int& minval, int& maxval,int type, float sigma)
{
	int midval;
	if (type == 0)
	{
		midval = getmatmidval(img);
	}
	else {
		midval = GetMatOTSU(img);
	}
	// 计算低阈值
	minval = saturate_cast<uchar>((1.0 - sigma) * midval);
	//计算高阈值
	maxval = saturate_cast<uchar>((1.0 + sigma) * midval);
}


int main() {
	Mat src = imread("1.png", 1);
	Mat gray;
	cvtColor(src, gray, COLOR_BGR2GRAY);
	GaussianBlur(gray, gray, Size(3, 3), 0.5, 0.5);
	int minthreshold, maxthreshold;
	GetMinMaxthreshold(gray, minthreshold, maxthreshold,0, 0.3);
	int minthreshold1, maxthreshold1;
	GetMinMaxthreshold(gray, minthreshold1, maxthreshold1, 1, 0.3);
	Mat dst1, dst2;
	Canny(gray, dst1, minthreshold, maxthreshold);
	Canny(gray, dst2, minthreshold1, maxthreshold1);
	imshow("1", dst1);
	imshow("2", dst2);
	waitKey();
	return 0;

}