一、演算法原理

       OSTU演算法也就是最大類間方差法，也即大津法。是一種選取最佳的閾值分割的方法，是閾值選取中最佳的方法。

按照灰度的特性將影象分成前景和背景兩部分。背景和前景之間的類間方差越大，說明構成影象兩部分的差別就

越大，當部分前景錯分為背景活部分背景錯分為前景都會導致兩部分差別變小。因此，是類間方差最大的分割

意味著錯分概率最小。

       設影象I（x，y）的前景和背景的分割閾值為T，屬於前景畫素點佔整幅影象的比例記為w1，其平均灰度為u1，

屬於背景畫素點佔整幅影象的比例記為w2，其平均灰度為u2。影象的總平均灰度記為u，類間方差記為g。假設

影象的大小為M*N，灰度值是0至255。具體步驟為：

1.首先計算每個畫素在整幅影象中的個數h[k]，然後把每個畫素值的個數除以影象的大小（M*N），得到每個畫素在

在整幅影象中的比例 p[i]。

2.求得前k個畫素值所佔的比例w[k]]和其均值u[k]，同時求的整幅影象畫素的均值uT,然後根據方差公式求的最大類間方差。

3.最大類方差公式為：

在這裡我只求了前景的比例和均值。因為根據方差公式可以推匯出最簡單的公式：

t=Max[(uT*w[k] - u[k])*(uT*w[k] - u[k]) / (w[k]*(1-w[k]))];推導過程如下：

二、ostu程式碼如下

#include <stdafx.h>
#include <math.h>   
#include <cv.h>   
#include <highgui.h>   
#include <stdio.h>   
int hei;   
int wid;    
void otsu(IplImage* A, IplImage* B)   
{   
    int i,j,k;
    long N = hei * wid;   
    int h[256];   
    double p[256],u[256],w[256];   
    for(i = 0; i < 256; i++)   
    {   
        h[i] = 0;   
        p[i] = 0;   
        u[i] = 0;  
        w[i] = 0;  

    }   
    for(i = 0; i < hei; i++)   
    {   
        for(j = 0; j < wid; j++)   
        {   
            for(k = 0; k < 256; k++)   
            {   
                if(((uchar*)(A->imageData + A->widthStep*i))[j] == k)   
                {   
                    h[k]++;   //統計灰度級中每個畫素在整幅影象中的個數
                }   
            }   
        }   
    }   
    for(i = 0; i < 256; i++)   
        p[i] = h[i] / double(N);   //計算每個畫素在整幅影象中的比例
    int T = 0; //閾值初始化為0  
    double uT,thegma2fang;   
    double thegma2fang_max = -10000;   
    for(int k = 0; k < 256; k++)   
    {   
        uT = 0;  
        for(i = 0; i <= k; i++)   
        {   
            u[k] += i*p[i];//u[k]表示前k個畫素值的均值 
            w[k] += p[i]; //w[k]表示前k個畫素值所佔的比例   
        } 
        for(i = 0; i < 256; i++)   
            uT += i*p[i];   //uT表示整幅影象的均值
        thegma2fang = (uT*w[k] - u[k])*(uT*w[k] - u[k]) / (w[k]*(1-w[k]));   
        if(thegma2fang > thegma2fang_max)   
        {   
            thegma2fang_max = thegma2fang;   
            T = k;
        }   
                   
    }  
printf("T=%d",T);//T=117
    for(i = 0; i < hei; i++)   
    {      
        for(j = 0; j < wid; j++)   
        {   
            if(((uchar*)(A->imageData + A->widthStep*i))[j] > T)   
                ((uchar*)(B->imageData + B->widthStep*i))[j] = 255;   
            else   
               ((uchar*)(B->imageData + B->widthStep*i))[j] = 0;   
        }   
    }   
           
}      
int main(int argc, char** argv)   
{   
    IplImage *src = cvLoadImage("lena.jpg"); 
    IplImage *InPut = cvCreateImage(cvSize(src->width,src->height),IPL_DEPTH_8U,1);   
    cvCvtColor(src,InPut,CV_BGR2GRAY);
    IplImage *OutPut = cvCreateImage(cvSize(src->width,src->height),IPL_DEPTH_8U,1);   
    hei = src->height;   
    wid = src->width;     
    otsu(InPut,OutPut);  
    cvNamedWindow( "Resource");   
    cvShowImage( "Resource", InPut );   
    cvNamedWindow( "Result");   
    cvShowImage( "Result", OutPut );   
    cvWaitKey(0);   
    return 0;   
}