转自：https://blog.csdn.net/dcrmg/article/details/52216622
大津法是一种图像灰度自适应的阈值分割算法，是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时退出
-
{
-
continue;
-
}
-
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;
-
}
原图像：
该幅图像计算出来的大津阈值是104；
用这个阈值分割的图像：
跟Opencv threshold方法中使用CV_THRESH_OTSU参数计算出来的分割图像一致：
直方图直观理解
大津算法可以从图像直方图上有一个更为直观的理解：大津阈值大致上是直方图两个峰值之间低谷的值。
对上述代码稍加修改，增加画出直方图部分：
-
#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};
//灰度直方图，下标是灰度值，保存内容是灰度值对应的像素点总数
-
int Histogram1[
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]]++;
-
Histogram1[data[i*image.step+j]]++;
-
}
-
}
-
-
//***********画出图像直方图********************************
-
Mat image1(255,255,CV_8UC3);
-
for(
int i=
0;i<
255;i++)
-
{
-
Histogram1[i]=Histogram1[i]%
200;
-
line(image1,Point(i,
235),Point(i,
235-Histogram1[i]),Scalar(
255,
0,
0),
1,
8,
0);
-
if(i%
50==
0)
-
{
-
char ch[
255];
-
sprintf(ch,
"%d",i);
-
string str=ch;
-
putText(image1,str,Point(i,
250),
1,
1,Scalar(
0,
0,
255));
-
}
-
}
-
//***********画出图像直方图********************************
-
-
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;
-
}
-
}
-
//画出以T为阈值的分割线
-
line(image1,Point(T,
235),Point(T,
0),Scalar(
0,
0,
255),
2,
8);
-
imshow(
"直方图",image1);
-
return T;
-
}
为显示清晰，本次使用一幅对比明显的灰度图：
OTSU分割效果：
对应阈值和直方图：
以上图像黑白对比度非常明显，从直方图上也可以看到只有两个波峰，求得的OTSU阈值为102。
上图中红色的竖线标识出了OTSu阈值分割线，显见，阈值大致位于两个波峰的低谷之间。
5742
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
