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)