otsuアルゴリズム実装(opecv 2.0バージョンベース)

出典:http://www.cnblogs.com/xy-kidult/p/3397553.html
otsuはnobuyuki otsuが1979年に提案した階調画像の最適しきい値を探すアルゴリズムである.
その論文名は『A Threshold Selection Method from Gray-Level Histograms』で、googleでpdf版を見つけることができ、大まかに論文を見て、プログラミングを実現することを提案し、重点を把握しやすい.私は最初はotsuの紹介を見てプログラミングを実現しましたが、求めたしきい値が理想的な効果を達成できませんでした.仕方なく、とうとう原作者の論文を読むことにした.
アルゴリズムの簡単な手順:
1各階調値の画像全体の数n[i]を算出する.i値を0~255からとる.
2各階調値が画像に現れる確率p[i],p[i]=n[i]/Nを計算し,Nは画像全体の画素数であり,画像の幅*が高い.
3 0~255から、各数を1つの閾値としてとり、そのときの前景と背景確率の総和とそれらの期待を計算し、計算式は下図を参照する.

4最後に、上記のいくつかの値に基づいて1つのクラス間分散とクラス内分散が得られ、式は置かれず、コードを直接見る.
 

#include <stdio.h>
#include <string>
#include "opencv2/highgui/highgui.hpp"
#include "opencv2/opencv.hpp"


#define MAX_GRAY_VALUE 256
#define MIN_GRAY_VALUE 0


int otsu(cv::Mat dst){

    int i,j;
    int tmp;

    double u0,u1,w0,w1,u, uk;
    
    double cov;
    double maxcov=0.0;
    int maxthread=0;

    int hst[MAX_GRAY_VALUE]={0};
    double pro_hst[MAX_GRAY_VALUE]={0.0};

    int height=dst.cols;
    int width=dst.rows;


    

    //         
    for( i =0 ; i<width; i++ ){
        for( j=0; j<height; j++){
            tmp=dst.at<uchar>(i,j);
            hst[tmp]++;
        }
    }

    //              
    for( i=MIN_GRAY_VALUE ; i<MAX_GRAY_VALUE; i++)
        pro_hst[i]=(double)hst[i]/(double)(width*height);
        

    //       
    u=0.0;
    for( i=MIN_GRAY_VALUE; i<MAX_GRAY_VALUE; i++)
        u += i*pro_hst[i];

    double det=0.0;
    for( i= MIN_GRAY_VALUE; i< MAX_GRAY_VALUE; i++)
        det += (i-u)*(i-u)*pro_hst[i];

    //             ，       
    
    for( i=MIN_GRAY_VALUE; i<MAX_GRAY_VALUE; i++){
        
        w0=0.0; w1=0.0; u0=0.0; u1=0.0; uk=0.0;

        for( j=MIN_GRAY_VALUE; j < i; j++){
            
            uk += j*pro_hst[j];
            w0 += pro_hst[j];
        
        }
        u0=uk/w0;

    
        w1=1-w0;
        u1= (u - uk )/(1-w0);
        

        //      
        cov=w0*w1*(u1-u0)*(u1-u0);

    
        

        if ( cov > maxcov )
        {
            maxcov=cov;
            maxthread=i;
        }
    }

    std::cout<<maxthread<<std::endl;
    return maxthread;
    

}


int main(){
    int width,height;
    int i,j;
    cv::Mat obj=cv::imread("Image16.jpg");

    cv::Mat dst;
    cv::cvtColor(obj,dst,CV_RGB2GRAY);

     height=dst.cols;
     width=dst.rows;
    int thd=otsu(dst);
    
    
    cv::imshow("origin img",dst);

    for( i=0; i < width; i++)
                for( j=0; j< height; j++)
                    if( dst.at<uchar>(i,j) > thd)
                        dst.at<uchar>(i,j)=MAX_GRAY_VALUE-1;
                    else
                        dst.at<uchar>(i,j)=MIN_GRAY_VALUE;

    
    
    cv::imshow("img after ostu",dst);

                
    cv::waitKey(0);
    return 0;
}