2. 程式人生 > >最大類間方差法（大津法OTSU）

最大類間方差法（大津法OTSU）

阿新 • • 發佈：2019-01-04

演算法介紹

最大類間方差法是1979年由日本學者大津提出的，是一種自適應閾值確定的方法，又叫大津法，簡稱OTSU，是一種基於全域性的二值化演算法，它是根據影象的灰度特性,將影象分為前景和背景兩個部分。當取最佳閾值時，兩部分之間的差別應該是最大的，在OTSU演算法中所採用的衡量差別的標準就是較為常見的最大類間方差。前景和背景之間的類間方差如果越大，就說明構成影象的兩個部分之間的差別越大，當部分目標被錯分為背景或部分背景被錯分為目標，都會導致兩部分差別變小，當所取閾值的分割使類間方差最大時就意味著錯分概率最小[1]。

記T為前景與背景的分割閾值，前景點數佔影象比例為w 0  ，平均灰度為u 0  ；背景點數佔影象比例為w

 1  ，平均灰度為u 1  ，影象的總平均灰度為u ，前景和背景圖象的方差，則有：

u=w 0 ×u 0 +w 1 ×u 1  
g=w 0 ×(u 0 −u) 2 +w 1 ×(u 1 −u) 2  
聯立上面兩式可得：
g=w 0 ×w 1 ×(u 0 −u 1 ) 2  
或
g=w 0 1−w 0 ×(u 0 −u) 2  
當方差g 最大時，可以認為此時前景和背景差異最大，此時的灰度T是最佳閾值。類間方差法對噪聲以及目標大小十分敏感，它僅對類間方差為單峰的影象產生較好的分割效果。當目標與背景的大小比例懸殊時（例如受光照不均、反光或背景複雜等因素影響），類間方差準則函式可能呈現雙峰或多峰，此時效果不好。直接用OTSU演算法處理自然場景銘牌圖片的部分結果例項如下：

1


圖1.最大類間方差法二值化結果 

程式碼實現

（C語言版，VS2012+opencv249）

#include "stdio.h"
#include "cv.h"
#include "highgui.h"
#include "Math.h"

int Otsu(IplImage* src);

int main()
{
    IplImage* img = cvLoadImage("lena.jpg",0); //獲取灰度影象img
    IplImage* dst = cvCreateImage(cvGetSize(img), 8, 1);
    int threshold = Otsu(img); //呼叫大津法求出最佳閾值
 

    printf("otsu threshold = %d\n", threshold);
    cvThreshold(img, dst, threshold, 255, CV_THRESH_BINARY); //用otsu的閾值二值化

    cvNamedWindow( "img", 1 );
    cvNamedWindow( "dst", 1 );
    cvShowImage("img", img);
    cvShowImage("dst", dst);


    cvWaitKey(-1);

    cvReleaseImage(&img);
    cvReleaseImage(&dst);

    cvDestroyWindow( "img" );
    cvDestroyWindow( "dst" );

    return 0;
}

int Otsu(IplImage* src)  
{  
    int height=src->height;  
    int width=src->width;      

    //histogram  
    float histogram[256] = {0};  
    for(int i=0; i < height; i++)
    {  
        unsigned char* p=(unsigned char*)src->imageData + src->widthStep * i;  
        for(int j = 0; j < width; j++) 
        {  
            histogram[*p++]++;  
        }  
    }  

    //normalize histogram & average pixel value 
    int size = height * width;  
    float u =0;
    for(int i = 0; i < 256; i++)
    {  
        histogram[i] = histogram[i] / size;  
        u += i * histogram[i];  //整幅影象的平均灰度
    }  

    int threshold;    
    float maxVariance=0;  
    float w0 = 0, avgValue  = 0;
    for(int i = 0; i < 256; i++) 
    {  
        w0 += histogram[i];  //假設當前灰度i為閾值, 0~i 灰度畫素所佔整幅影象的比例即前景比例
        avgValue  += i * histogram[i]; //avgValue/w0 = u0

        float t = avgValue/w0 - u;  //t=u0-u
        float variance = t * t * w0 /(1 - w0);  
        if(variance > maxVariance) 
        {  
            maxVariance = variance;  
            threshold = i;  
        }  
    }  

    return threshold;  
}
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26
  • 27
  • 28
  • 29
  • 30
  • 31
  • 32
  • 33
  • 34
  • 35
  • 36
  • 37
  • 38
  • 39
  • 40
  • 41
  • 42
  • 43
  • 44
  • 45
  • 46
  • 47
  • 48
  • 49
  • 50
  • 51
  • 52
  • 53
  • 54
  • 55
  • 56
  • 57
  • 58
  • 59
  • 60
  • 61
  • 62
  • 63
  • 64
  • 65
  • 66
  • 67
  • 68
  • 69
  • 70
  • 71
  • 72
  • 73
  • 74
  • 75
  • 76

程式碼執行結果：
2

圖2.otsu程序運行結果 

參考資料：
[1]Otsu N. A threshold selection method from gray-level histograms[J]. Automatica, 1975, 11(285-296): 23-27. 

OTSU演算法：就是計算出灰度圖最佳閾值的演算法

1.先對灰度圖進行直方圖計算並歸一化處理，得到0-255之間每個畫素在灰度圖中出現的概率，即表示為某個畫素在灰度圖中出現了n個，灰度圖總的畫素點為N個，則這個畫素的出現概率為Pi=n/N

2.每個灰度圖可以由閾值k將灰度圖分為A，B兩大類，很容易得到A，B類在灰度圖中的出現概率以及灰度均值

3.計算灰度圖A,B類得類間方差，在最佳閾值K處，求得的類間方差最大，也就是類間方差最大的那個時刻的閾值就為灰度圖的最佳閾值

程式碼實現：

以下三種演算法時間複雜度依次降低，前兩種演算法自己所寫，最後一種演算法從網上看到大神寫的，膜拜！當然如果你想偷懶，可是更直接點用cvThreshold一步實現cvThreshold(srcCopy,thresholdImage,0,255,CV_THRESH_OTSU)，花費時間與演算法三相當。

演算法一：

處理一張640*480的圖片花費時間大概在8000毫秒左右，非常低效，但是是一種非常直接的，顯而易見的實現方法。

[cpp] view plain copy print?
  1. int otsu(IplImage* src)  
  2. {  
  3.     CvScalar cs = cvSum(src);                             //影象總灰度值  
  4.     int Width = src->width;  
  5.     int Height = src->height;  
  6.     double U_t = 1.0*cs.val[0]/(Width*Height);                  //影象的平均灰度  
  7.     int threshold = 0;  
  8.     double delta = 0;  
  9.     for(int k=0; k<256; k++)  
  10.     {  
  11.         double sumOfGray = 0;  
  12.         int num = 0;  
  13.         for(int i=0; i < Width; i++)  
  14.         {  
  15.             for(int j=0; j < Height;j++)  
  16.             {  
  17.                 int t = cvGet2D(src, j,i).val[0];  
  18.                 if(t < k)  
  19.                 {  
  20.                     num++;             //統計灰度小於k的畫素點的個數  
  21.                     sumOfGray += t; //灰度值小於k的畫素點的灰度和  
  22.                 }  
  23.             }  
  24.         }  
  25.         double w0 = 1.0*num/(Width*Height);                     //灰度值小於閾值k的畫素的概率  
  26.         double w1 = 1.0 - w0;                               //灰度值大於k畫素的概率  
  27.         double u0 = 1.0*sumOfGray/num;                     //灰度值小於閾值k的畫素的平均灰度值  
  28.         double u1 = 1.0*(cs.val[0]-sumOfGray)/(Width*Height-num); //灰度值大於閾值k的畫素的平均灰度值  
  29.         double delta_tmp = w0*pow(u0-U_t,2) + w1*pow(u1-U_t,2);    //類間方差  
  30.         if(delta_tmp > delta)  
  31.         {  
  32.             delta = delta_tmp;  
  33.             threshold = k;  
  34.         }  
  35.     }  
  36.     cout<<cs.val[0]<<endl;  
  37.     cout<<U_t<<endl;  
  38.     return threshold;  
  39. }  
int otsu(IplImage* src)
{
	CvScalar cs = cvSum(src);                             //影象總灰度值
	int Width = src->width;
	int Height = src->height;
	double U_t = 1.0*cs.val[0]/(Width*Height);                  //影象的平均灰度
	int threshold = 0;
	double delta = 0;
	for(int k=0; k<256; k++)
	{
		double sumOfGray = 0;
		int num = 0;
		for(int i=0; i < Width; i++)
		{
			for(int j=0; j < Height;j++)
			{
				int t = cvGet2D(src, j,i).val[0];
				if(t < k)
				{
					num++;             //統計灰度小於k的畫素點的個數
					sumOfGray += t; //灰度值小於k的畫素點的灰度和
				}
			}
		}
		double w0 = 1.0*num/(Width*Height);                     //灰度值小於閾值k的畫素的概率
		double w1 = 1.0 - w0;                               //灰度值大於k畫素的概率
		double u0 = 1.0*sumOfGray/num;                     //灰度值小於閾值k的畫素的平均灰度值
		double u1 = 1.0*(cs.val[0]-sumOfGray)/(Width*Height-num); //灰度值大於閾值k的畫素的平均灰度值

		double delta_tmp = w0*pow(u0-U_t,2) + w1*pow(u1-U_t,2);    //類間方差
		
		if(delta_tmp > delta)
		{
			delta = delta_tmp;
			threshold = k;
		}
	}
	cout<<cs.val[0]<<endl;
	cout<<U_t<<endl;
	return threshold;
}

演算法二：

借用opencv中的直方圖函式，大大減少迴圈數，同一張圖片花費時間大概在16毫秒左右。

[cpp] view plain copy print?
  1. int otsu(IplImage* src)    
  2. {     
  3.     int hist_size = 256;    //直方圖尺寸      
  4.     int hist_height = 256;      
  5.     float range[] = {0,255};  //灰度級的範圍      
  6.     float* ranges[]={range};      
  7.     //建立一維直方圖，統計影象在[0 255]畫素的均勻分佈      
  8.     CvHistogram* gray_hist = cvCreateHist(1,&hist_size,CV_HIST_ARRAY,ranges,1);      
  9.     //計算灰度影象的一維直方圖      
  10.     cvCalcHist(&src,gray_hist,0,0);      
  11.     //歸一化直方圖      
  12.     cvNormalizeHist(gray_hist,1.0);      
  13.     int Width = src->width;    
  14.     int Height = src->height;     
  15.     int threshold = 0;    
  16.     double delta = 0;  
  17.     CvScalar sumOfGray = cvSum(src);  
  18.     double U_t = sumOfGray.val[0]/(Width*Height);  
  19.     for(int k=0; k<256; k++)    
  20.     {    
  21.         double u0 = 0,  u1 = 0, w0 = 0, w1 = 0, delta_tmp = 0, u00 = 0, u11 = 0;    
  22.         for(int i=0; i<k; i++)    
  23.         {    
  24.             u00 += cvQueryHistValue_1D(gray_hist,i)*i;   //灰度小於閾值k的畫素的平均灰度值    
  25.             w0 += cvQueryHistValue_1D(gray_hist,i);     //灰度小於閾值k的畫素的概率    
  26.         }  
  27.         u0 = u00/w0;  
  28.         for(int j=k; j<256; j++)    
  29.         {    
  30.             u11 += cvQueryHistValue_1D(gray_hist,j)*j;   //灰度大於閾值k的畫素的平均灰度值    
  31.             w1 += cvQueryHistValue_1D(gray_hist,j);     //灰度小於閾值k的畫素的概率    
  32.         }  
  33.         u1 = u11/w1;  
  34.         delta_tmp = w0*pow(u0-U_t,2) + w1*pow(u1-U_t,2);    //類間方差    
  35.         if(delta_tmp > delta)    
  36.         {    
  37.             delta = delta_tmp;    
  38.             threshold = k;    
  39.         }    
  40.     }   
  41.     cvReleaseHist(&gray_hist);  
  42.     return threshold;     
  43. }  
int otsu(IplImage* src)  
{   
    int hist_size = 256;    //直方圖尺寸    
    int hist_height = 256;    
    float range[] = {0,255};  //灰度級的範圍    
    float* ranges[]={range};    
    //建立一維直方圖，統計影象在[0 255]畫素的均勻分佈    
    CvHistogram* gray_hist = cvCreateHist(1,&hist_size,CV_HIST_ARRAY,ranges,1);    
    //計算灰度影象的一維直方圖    
    cvCalcHist(&src,gray_hist,0,0);    
    //歸一化直方圖    
    cvNormalizeHist(gray_hist,1.0);    
      
    int Width = src->width;  
    int Height = src->height;   
    int threshold = 0;  
    double delta = 0;
	CvScalar sumOfGray = cvSum(src);
	double U_t = sumOfGray.val[0]/(Width*Height);
	
    for(int k=0; k<256; k++)  
    {  
        double u0 = 0,  u1 = 0, w0 = 0, w1 = 0, delta_tmp = 0, u00 = 0, u11 = 0;  
         
        for(int i=0; i<k; i++)  
        {  
            u00 += cvQueryHistValue_1D(gray_hist,i)*i;   //灰度小於閾值k的畫素的平均灰度值  
            w0 += cvQueryHistValue_1D(gray_hist,i);     //灰度小於閾值k的畫素的概率  
        }
		
		u0 = u00/w0;

        for(int j=k; j<256; j++)  
        {  
            u11 += cvQueryHistValue_1D(gray_hist,j)*j;   //灰度大於閾值k的畫素的平均灰度值  
            w1 += cvQueryHistValue_1D(gray_hist,j);     //灰度小於閾值k的畫素的概率  
        }

		u1 = u11/w1;

        delta_tmp = w0*pow(u0-U_t,2) + w1*pow(u1-U_t,2);    //類間方差  
        if(delta_tmp > delta)  
        {  
            delta = delta_tmp;  
            threshold = k;  
        }  
    } 

	cvReleaseHist(&gray_hist);
    return threshold;   
}

演算法三：

一次迴圈實現，其中的公式需要自己用筆推導一下，非常精妙！完成一次運算的時間毫秒已經不能計量了，顯示為0！

[cpp] view plain copy print?
  1. int otsu(IplImage* src)    
  2. {     
  3.     int hist_size = 256;    //直方圖尺寸      
  4.     int hist_height = 256;      
  5.     float range[] = {0,255};  //灰度級的範圍      
  6.     float* ranges[]={range};      
  7.     //建立一維直方圖，統計影象在[0 255]畫素的均勻分佈      
  8.     CvHistogram* gray_hist = cvCreateHist(1,&hist_size,CV_HIST_ARRAY,ranges,1);      
  9.     //計算灰度影象的一維直方圖      
  10.     cvCalcHist(&src,gray_hist,0,0);      
  11.     //歸一化直方圖      
  12.     cvNormalizeHist(gray_hist,1.0);      
  13.     int Width = src->width;    
  14.     int Height = src->height;     
  15.     int threshold = 0;    
  16.     double delta = 0;  
  17.     double U_t = 0;  
  18.     for(int m=0; m<256; m++)  
  19.     {  
  20.         U_t += cvQueryHistValue_1D(gray_hist,m)*m;  
  21.     }  
  22.     double u = 0, w = 0;   
  23.     for(int k=0; k<256; k++)    
  24.     {    
  25.         u += cvQueryHistValue_1D(gray_hist,k)*k;   //    
  26.         w += cvQueryHistValue_1D(gray_hist,k);      //灰度大於閾值k的畫素的概率    
  27.         double t = U_t * w - u;  
  28.         double delta_tmp = t * t / (w * (1 - w) );   
  29.         if(delta_tmp > delta)    
  30.         {    
  31.             delta = delta_tmp;    
  32.             threshold = k;    
  33.         }    
  34.     }   
  35.     cvReleaseHist(&gray_hist);  
  36.     return threshold;     
  37. }   

相關推薦

最大類間方差法、大津法（ostu）

一、演算法原理        OSTU演算法也就是最大類間方差法，也即大津法。是一種選取最佳的閾值分割的方法，是閾值選取中最佳的方法。 按照灰度的特性將影象分成前景和背景兩部分。背景和前景之間的類間方差越大，說明構成影象兩部分的差別就 越大，當部分

Otsu（大津法，最大類間方差法）

雖然很早就看過這本綠皮書，但是當時是剛入門的小菜鳥，根本就不知道這就是大名鼎鼎的大津法，當時只是覺得Otsu好奇怪的英文名字。 現在就來重新看看這個所謂的大津法: nin_ini​表示灰度級為i的畫素數。影象中的畫素總數MN為MN=n0+n1+...+nL−1M

最大類間方差法（大津法OTSU）

演算法介紹 最大類間方差法是1979年由日本學者大津提出的，是一種自適應閾值確定的方法，又叫大津法，簡稱OTSU，是一種基於全域性的二值化演算法，它是根據影象的灰度特性,將影象分為前景和背景兩個部分。當取最佳閾值時，兩部分之間的差別應該是最大的，在OTSU演算法中所採用的衡量差別的標準就是較為常見的最大類間

自適應閾值分割(最大類間方差法、大津法、OTSU)

最大類間方差法是由日本學者大津(Nobuyuki Otsu)於1979年提出的,是一種自適應的閾值確定的方法,又叫大津法,簡稱OTSU。它是按影象的灰度特性,將影象分成背景和目標兩部分，或者說，是尋找一個閾值為K，將影象的顏色分為1,2.....K和K+1.....256

MATLAB實現最大類間方差演算法

Otsu演算法(大律法或最大類間方差法) 最大類間方差法是由日本學者大津(Nobuyuki Otsu)於1979年提出的,是一種自適應的閾值確定的方法,又叫大津法,簡稱OTSU。它是按影象的灰度特性,將影象分成背景和目標2部分。背景和目標之間的類間方差越大,說

最大方差和最小協方差解釋（線性代數看PCA）

數學基礎 限制 兩個 ali 從大到小 其中 我們 特征 向上 轉自： 最大方差和最小協方差解釋（線性代數看PCA） PCA降維 ——最大方差和最小協方差聯合解釋（線性代數看PCA） 註：根據網上資料整理而得，歡迎討論

Ostu最大類間差方法

Ostu方法又名最大類間差方法，通過統計整個影象的直方圖特性來實現全域性閾值T的自動選取，其演算法步驟為： 1) 先計算影象的直方圖，即將影象所有的畫素點按照0~255共256個bin，統計落在每

斯坦福大學公開課機器學習：advice for applying machine learning | learning curves （改進學習算法：高偏差和高方差與學習曲線的關系）

繪制 學習曲線 pos 情況 但我 容量 繼續 並且 inf 繪制學習曲線非常有用，比如你想檢查你的學習算法，運行是否正常。或者你希望改進算法的表現或效果。那麽學習曲線就是一種很好的工具。學習曲線可以判斷某一個學習算法，是偏差、方差問題，或是二者皆有。 為了繪制一條學習曲

【51NOD-0】1089 最長回文子串 V2（Manacher算法）

lose 最長回文子串 gif () none print struct hide pac 【算法】回文樹 #include<cstdio> #include<algorithm> #include<cstring> using na

51nod 1089 最長回文子串 V2（Manacher算法）

clu 記得 file 越界 str tool algorithm scanf ++i 1089 最長回文子串 V2（Manacher算法） 基準時間限制：1 秒 空間限制：131072 KB 分值: 0 難度：基礎題 收藏 關註

求一般圖的最小頂點覆蓋集問題的混合貪婪算法（近似算法）

混合 論文 之前 blog tor col mes 偏差 cnblogs 之前準備做hiho一下的時候，網上查關於無向圖的最大獨立集； 看到了一篇論文，說是能“求一般圖的最小頂點覆蓋集問題”的混合貪婪算法； 我一看覺得挺牛逼的啊，跑去研究了大半天的這篇論文，發現實際還是一個

最簡單的排序算法（C和C++實現）

最簡單的排序算法（C和C++實現）1、算法思想如下圖：把待排序的數都存在對應的數組的下標中，如果待排序的數有重復的，就用對應的數組加一，最後把數組的下標打印出來即可。2、源碼（C）如下：#include <stdio.h>int main (void){ int a[100], i, j,

淺談協方差矩陣（轉）

但是 資料 sample 可能 tlab 原來 data- 一個 處理 原文地址：http://pinkyjie.com/2010/08/31/covariance/ 今天看論文的時候又看到了協方差矩陣這個破東西，以前看模式分類的時候就特困擾，沒想到現在還是搞不清楚，索性開

求兩個數的最大公約數，輾轉相除法與更相減損法（遞歸叠代）

叠代 div 余數 公約數 穩定 log test 算法 復雜度 問題：給出兩個數a和b，求出他們的最大公約數(greatest common divisor)。 解法一：輾轉相除法，又叫歐幾裏得算法。兩個正整數a和b（a>b），他們的最大公約數等於a除以b的余數和b

統計學中的協方差矩陣（陣列訊號基礎）

在處理陣列訊號的時候，為了獲得空間訊號維度的相關性，以估計目標的資訊。故使用協方差矩陣能夠獲得這些，因為協方差矩陣是每一維度下（也就是陣元）訊號的相關性。當兩個維度相關時，訊號的協方差也是最大的。   一、統計學的基本概念   統計學裡最基本的概念就是樣本的均值、

懲罰函式法（內點法、外點法）求解約束優化問題最優值 matlab

1、 用外點法求下列問題的最優解 方法一：外點牛頓法： clc m=zeros(1,50);a=zeros(1,50);b=zeros(1,50);f0=zeros(1,50);%a b為最優點座標，f0為最優點函式值，f1 f2最優點梯度。 syms x1 x2 e;  &n

多階矩在影象中的含義（方差，偏度，峰度）

設影象共有N點（影象塊則為w*h個畫素點），第i點的灰度值為Xi，其均值為X，則這些特徵的含義如下： 方差（Variance）：是一種衡量樣本分佈均勻性的尺度標準。 計算公式為： 偏度（Skew

協方差實現（python）

import numpy as np if __name__ == '__main__': a = np.array([[1,2,3],[3,1,1]]) c_ = np.mean(a, axis =0) #按列計算 x = a - c_ # 廣播機制 cov =

字串最小表示法（LeetCode 899. Orderly Queue ）

A string S of lowercase letters is given.  Then, we may make any number of moves. In each move, we choose one of the first K letters (sta

輸出從大到小排序的成績（氣泡排序法）

使用氣泡排序法： #include <iostream> usingnamespace std; int main()   //輸出從大到小排序的成績 {     int n