OSTU原理及MATLAB和C++算法实现

对于一幅图像，设当前景与背景的分割阈值为 t 时，前景点占图像比例为 w0，均值为 u0，背景点占图像比例为 w1，均值为 u1。则整个图像的均值为 u = w0*u0+w1*u1。建立目标函数g(t)=w0*(u0-u)^2+w1*(u1-u)^2，g(t)就是当分割阈值为t时的类间方差表达式。OTSU算法使得g(t)取得全局最大值，当g(t)为最大时所对应的t称为最佳阈值。OTSU算法又称为最大类间方差法。

最初一直不懂什么叫类间方差，其实用通俗的语言来说就是，把1到255这些像素分成两类，如1到125为一类，125到255为一类，因为分类可以分成很多种，所以可以找到一种最大区分开两类的一种方法就是OSTU算法了，而怎么找到最大的分类方法呢？就是按照 u = w0*u0+w1*u1，以及g(t)=w0*(u0-u)^2+w1*(u1-u)^2来分类了。

#ifndef _OSTU_
#define _OSTU_
#include <cv.h>
#include <highgui.h>
using namespace cv;
//输入图像，输出阈值
namespace OSTU{
int cvOSTU(IplImage *image)
{
assert(NULL != image);
int width = image->width;
int height = image->height;
int x = 0, y = 0;
int pixelCount[256];
float pixelPro[256];
int i, j, pixelSum = width * height, threshold = 0;
uchar* data = (uchar*)image->imageData;
//初始化
for (i = 0; i < 256; i++)
{
pixelCount[i] = 0;
pixelPro[i] = 0;
}
//统计灰度级中每个像素在整幅图像中的个数
for (i = y; i < height; i++)
{
for (j = x; j < width; j++)
{
pixelCount[data[i * image->widthStep + j]]++;
}
}
//计算每个像素在整幅图像中的比例
for (i = 0; i < 256; i++)
{
pixelPro[i] = (float)(pixelCount[i]) / (float)(pixelSum);
}
//经典ostu算法,得到前景和背景的分割
//遍历灰度级[0,255],计算出方差最大的灰度值,为最佳阈值
float w0, w1, u0tmp, u1tmp, u0, u1, u, deltaTmp, deltaMax = 0;
for (i = 0; i < 256; i++)
{
w0 = w1 = u0tmp = u1tmp = u0 = u1 = u = deltaTmp = 0;
for (j = 0; j < 256; j++)
{
if (j <= i) //背景部分
{
//以i为阈值分类，第一类总的概率
w0 += pixelPro[j];
u0tmp += j * pixelPro[j];
}
else       //前景部分
{
//以i为阈值分类，第二类总的概率
w1 += pixelPro[j];
u1tmp += j * pixelPro[j];
}
}
u0 = u0tmp / w0;        //第一类的平均灰度
u1 = u1tmp / w1;        //第二类的平均灰度
u = u0tmp + u1tmp;      //整幅图像的平均灰度
//计算类间方差
deltaTmp = w0 * (u0 - u)*(u0 - u) + w1 * (u1 - u)*(u1 - u);
//找出最大类间方差以及对应的阈值
if (deltaTmp > deltaMax)
{
deltaMax = deltaTmp;
threshold = i;
}
}
//返回最佳阈值;
return threshold;
}

int OSTU(Mat image){
return cvOSTU(cvCloneImage(&(IplImage)image));
}
};
#endif // !_OSTU_

close all;
clear all;
clc;
G = imread('D:\test.bmp');
figure;
imshow(G);
title('原图');
%I = rgb2gray(G);
I=G;
[m,n] = size(I);
Hist = zeros(256);%直方图
dHist = zeros(256);%各像素值所占比例
variance = zeros(256);%方差
PXD = 0;%阈值初始化
for i = 1:m
for j = 1:n
Hist(uint8(I(i,j))+1) = Hist(uint8(I(i,j))+1) + 1;%由于matlab矩阵下标最小为1，所以这里稍微变换一下取下标为uint8(I(i,j))+1代表灰度值uint8(I(i,j))对应的直方图下标
end
end
for i = 1:256
dHist(i) = Hist(i)/(m*n);
end
for PXD = 1:256
w0 = 0;
w1 = 0;
g0 = 0;
g1 = 0;
for i = 1:PXD
g0 = g0 + (i-1)*dHist(i);
w0 = w0 + dHist(i);
end
for i = PXD+1 : 256
g1 = g1 + (i-1)*dHist(i);
w1 = w1 + dHist(i);
end
variance(PXD) = w0*w1*(g0 - g1)*(g0 - g1);%考虑到速度问题，此处用了一个等价算法，其效果等同于最大化类间方差
end
PXD = 1;
for i = 1:256
if variance(PXD) < variance(i)
PXD = i;
end
end
for i = 1:m
for j = 1:n
if I(i,j) > PXD
I(i,j) = 255;
else
I(i,j) = 0;
end
end
end
imageBW = I;
figure;
imshow(imageBW);
title('Ostu二值化');

此外，matlab中还有现成的Ostu阈值分割函数​，用法如下：
I = imread('D:\test.bmp');
level = graythresh(I);    %也就是原理中循环查找使得类间方差最大化的阈值步骤
BW = im2bw(I,level);        %找到阈值二值化即可
figure, imshow(BW)