CLAHE的實現和研究

CLAHE算法對於醫學圖像，特別是醫學紅外圖像的增強效果非常明顯。

CLAHE  https://en.wikipedia.org/wiki/Adaptive_histogram_equalization

中文方面非常好的資料 限制對比度自適應直方圖均衡化算法原理、實現及效果

在OpenCV中已經實現了CLAHE,但是它在使用過程中，存在參數選擇的問題。為了從根本上搞明白，我參考了網絡上的一些代碼

主要是來源 http://blog.csdn.net/abcd1992719g/article/details/25483395

實現了基於OpenCV的CLAHE實現和研究。從最基本的開始做，分別實現HE算法，AHE算法，CLHE算法和CLAHE算法。素材分別采用了手部和手臂的紅外圖片，同時調用OpenCV生成代碼和自己編寫代碼進行比對。

調用代碼和實現效果:

int _tmain( int argc, _TCHAR * argv[])
{
     //讀入灰度的手部圖像
    Mat src  = imread( "arm.jpg", 0);
    Mat dst  = src.clone();
    Mat HT_OpenCV;
    Mat HT_GO;
    Mat AHE_GO;
    Mat CLHE_GO;
    Mat CLAHE_Without_Interpolation;
    Mat CLAHE_OpenCV;
    Mat CLAHE_GO;
    Mat matInter;
     ////OpenCV HT 方法
    cv : :equalizeHist(src,HT_OpenCV);
     ////GO HT方法
    HT_GO  = eaualizeHist_GO(src);
     ////GO AHE方法
    AHE_GO  = aheGO(src);
     ////GO CLHE方法
    CLHE_GO  = clheGO(src);
     ////clahe不計算差值
    CLAHE_Without_Interpolation  = claheGoWithoutInterpolation(src);
     ////OpenCV CLAHE 方法
    Ptr <cv : :CLAHE > clahe  = createCLAHE(); //默認參數
    clahe - >apply(src, CLAHE_OpenCV);
     ////GO CLAHE方法
    CLAHE_GO  = claheGO(src);
 
     ////結果顯示
    imshow( "原始圖像",src);
    imshow( "OpencvHT",HT_OpenCV);
    imshow( "GOHT",HT_GO);
    imshow( "GOAHE",AHE_GO);
    imshow( "GOCLHE",CLHE_GO);
    imshow( "GOCLAHE",CLAHE_GO);
    imshow( "CLAHE_Without_Interpolation",CLAHE_Without_Interpolation);
    imshow( "OpencvCLAHE",CLAHE_OpenCV);
    waitKey();
     return  0;
}

原始圖像

GOCLAHE效果

OpenCV CLAHE效果

HE算法： Mat eaualizeHist_GO(Mat src)

{
     int width  = src.cols;
     int height = src.rows;
    Mat HT_GO  = src.clone();
     int tmp[ 256]  ={ 0};
     float C[ 256]  = { 0. 0};
     int total  = width *height;  
     for ( int i = 0 ;i <src.rows;i ++)
    {
         for ( int j = 0;j <src.cols;j ++)
        {
             int index  = src.at <uchar >(i,j);
            tmp[index]  ++;
        }
    }
     //計算累積函數  
     for( int i  =  0;i  <  256 ; i ++){  
         if(i  ==  0)  
            C[i]  =  1.0f  * tmp[i]  / total;  
         else  
            C[i]  = C[i - 1]  +  1.0f  * tmp[i]  / total;  
    }  
     //這里的累積函數分配的方法非常直觀高效
     for( int i  =  0;i  < src.rows;i ++){  
         for( int j  =  0;j  < src.cols;j ++){      
             int index  = src.at <uchar >(i,j);
            HT_GO.at <uchar >(i,j)  = C[index]  *  255  ;
        }  
    }  
     return HT_GO;
}

 

AHE算法:

Mat aheGO(Mat src, int _step  =  8)
{
    Mat AHE_GO  = src.clone();
     int block  = _step;
     int width  = src.cols;
     int height  = src.rows;
     int width_block  = width /block;  //每個小格子的長和寬
     int height_block  = height /block;
     //存儲各個直方圖  
     int tmp2[ 8 * 8][ 256]  ={ 0};
     float C2[ 8 * 8][ 256]  = { 0. 0};
     //分塊
     int total  = width_block  * height_block; 
     for ( int i = 0;i <block;i ++)
    {
         for ( int j = 0;j <block;j ++)
        {
             int start_x  = i *width_block;
             int end_x  = start_x  + width_block;
             int start_y  = j *height_block;
             int end_y  = start_y  + height_block;
             int num  = i +block *j;  
             //遍歷小塊,計算直方圖
             for( int ii  = start_x ; ii  < end_x ; ii ++)  
            {  
                 for( int jj  = start_y ; jj  < end_y ; jj ++)  
                {  
                     int index  =src.at <uchar >(jj,ii);
                    tmp2[num][index] ++;  
                }  
            } 
             //計算累積分布直方圖  
             for( int k  =  0 ; k  <  256 ; k ++)  
            {  
                 if( k  ==  0)  
                    C2[num][k]  =  1.0f  * tmp2[num][k]  / total;  
                 else  
                    C2[num][k]  = C2[num][k - 1]  +  1.0f  * tmp2[num][k]  / total;  
            }  
        }
    }
     //將統計結果寫入
     for ( int i = 0;i <block;i ++)
    {
         for ( int j = 0;j <block;j ++)
        {
             int start_x  = i *width_block;
             int end_x  = start_x  + width_block;
             int start_y  = j *height_block;
             int end_y  = start_y  + height_block;
             int num  = i +block *j;  
             //遍歷小塊,計算直方圖
             for( int ii  = start_x ; ii  < end_x ; ii ++)  
            {  
                 for( int jj  = start_y ; jj  < end_y ; jj ++)  
                {  
                     int index  =src.at <uchar >(jj,ii);
                     //結果直接寫入AHE_GO中去
                    AHE_GO.at <uchar >(jj,ii)  = C2[num][index]  *  255  ;
                }  
            } 
        }
    }
     return AHE_GO;
}

CLHE算法：

//這里是在全局直方圖加入“限制對比度”方法
Mat clheGO(Mat src, int _step  =  8)
{
     int width  = src.cols;
     int height = src.rows;
    Mat CLHE_GO  = src.clone();
     int tmp[ 256]  ={ 0};
     float C[ 256]  = { 0. 0};
     int total  = width *height;  
     for ( int i = 0 ;i <src.rows;i ++)
    {
         for ( int j = 0;j <src.cols;j ++)
        {
             int index  = src.at <uchar >(i,j);
            tmp[index]  ++;
        }
    }
     /////////////////////////限制對比度計算部分,注意這個地方average的計算不一定科學
     int average  = width  * height  /  255 / 64;  
     int LIMIT  =  4  * average;  
     int steal  =  0;  
     for( int k  =  0 ; k  <  256 ; k ++)  
    {  
         if(tmp[k]  > LIMIT){  
            steal  += tmp[k]  - LIMIT;  
            tmp[k]  = LIMIT;  
        }  
    }  
     int bonus  = steal / 256;  
     //hand out the steals averagely  
     for( int k  =  0 ; k  <  256 ; k ++)  
    {  
        tmp[k]  += bonus;  
    }  
     ///////////////////////////////////////////
     //計算累積函數  
     for( int i  =  0;i  <  256 ; i ++){  
         if(i  ==  0)  
            C[i]  =  1.0f  * tmp[i]  / total;  
         else  
            C[i]  = C[i - 1]  +  1.0f  * tmp[i]  / total;  
    }  
     //這里的累積函數分配的方法非常直觀高效
     for( int i  =  0;i  < src.rows;i ++){  
         for( int j  =  0;j  < src.cols;j ++){      
             int index  = src.at <uchar >(i,j);
            CLHE_GO.at <uchar >(i,j)  = C[index]  *  255  ;
        }  
    }  
     return CLHE_GO;
}

CLAHE不包括插值算法：

Mat claheGoWithoutInterpolation(Mat src,  int _step  =  8)
{
    Mat CLAHE_GO  = src.clone();
     int block  = _step; //pblock
     int width  = src.cols;
     int height = src.rows;
     int width_block  = width /block;  //每個小格子的長和寬
     int height_block  = height /block;
     //存儲各個直方圖  
     int tmp2[ 8 * 8][ 256]  ={ 0};
     float C2[ 8 * 8][ 256]  = { 0. 0};
     //分塊
     int total  = width_block  * height_block; 
     for ( int i = 0;i <block;i ++)
    {
         for ( int j = 0;j <block;j ++)
        {
             int start_x  = i *width_block;
             int end_x  = start_x  + width_block;
             int start_y  = j *height_block;
             int end_y  = start_y  + height_block;
             int num  = i +block *j;  
             //遍歷小塊,計算直方圖
             for( int ii  = start_x ; ii  < end_x ; ii ++)  
            {  
                 for( int jj  = start_y ; jj  < end_y ; jj ++)  
                {  
                     int index  =src.at <uchar >(jj,ii);
                    tmp2[num][index] ++;  
                }  
            } 
             //裁剪和增加操作，也就是clahe中的cl部分
             //這里的參數 對應《Gem》上面 fCliplimit  = 4  , uiNrBins  = 255
             int average  = width_block  * height_block  /  255;  
             int LIMIT  =  4  * average;  
             int steal  =  0;  
             for( int k  =  0 ; k  <  256 ; k ++)  
            {  
                 if(tmp2[num][k]  > LIMIT){  
                    steal  += tmp2[num][k]  - LIMIT;  
                    tmp2[num][k]  = LIMIT;  
                }  
            }  
             int bonus  = steal / 256;  
             //hand out the steals averagely  
             for( int k  =  0 ; k  <  256 ; k ++)  
            {  
                tmp2[num][k]  += bonus;  
            }  
             //計算累積分布直方圖  
             for( int k  =  0 ; k  <  256 ; k ++)  
            {  
                 if( k  ==  0)  
                    C2[num][k]  =  1.0f  * tmp2[num][k]  / total;  
                 else  
                    C2[num][k]  = C2[num][k - 1]  +  1.0f  * tmp2[num][k]  / total;  
            }  
        }
    }
     //計算變換后的像素值  
     //將統計結果寫入
     for ( int i = 0;i <block;i ++)
    {
         for ( int j = 0;j <block;j ++)
        {
             int start_x  = i *width_block;
             int end_x  = start_x  + width_block;
             int start_y  = j *height_block;
             int end_y  = start_y  + height_block;
             int num  = i +block *j;  
             //遍歷小塊,計算直方圖
             for( int ii  = start_x ; ii  < end_x ; ii ++)  
            {  
                 for( int jj  = start_y ; jj  < end_y ; jj ++)  
                {  
                     int index  =src.at <uchar >(jj,ii);
                     //結果直接寫入AHE_GO中去
                    CLAHE_GO.at <uchar >(jj,ii)  = C2[num][index]  *  255  ;
                }  
            } 
        }
    
     }  
     return CLAHE_GO;
}

CLAHE算法：

Mat claheGO(Mat src, int _step  =  8)
{
    Mat CLAHE_GO  = src.clone();
     int block  = _step; //pblock
     int width  = src.cols;
     int height = src.rows;
     int width_block  = width /block;  //每個小格子的長和寬
     int height_block  = height /block;
     //存儲各個直方圖  
     int tmp2[ 8 * 8][ 256]  ={ 0};
     float C2[ 8 * 8][ 256]  = { 0. 0};
     //分塊
     int total  = width_block  * height_block; 
     for ( int i = 0;i <block;i ++)
    {
         for ( int j = 0;j <block;j ++)
        {
             int start_x  = i *width_block;
             int end_x  = start_x  + width_block;
             int start_y  = j *height_block;
             int end_y  = start_y  + height_block;
             int num  = i +block *j;  
             //遍歷小塊,計算直方圖
             for( int ii  = start_x ; ii  < end_x ; ii ++)  
            {  
                 for( int jj  = start_y ; jj  < end_y ; jj ++)  
                {  
                     int index  =src.at <uchar >(jj,ii);
                    tmp2[num][index] ++;  
                }  
            } 
             //裁剪和增加操作，也就是clahe中的cl部分
             //這里的參數 對應《Gem》上面 fCliplimit  = 4  , uiNrBins  = 255
             int average  = width_block  * height_block  /  255;  
             //關於參數如何選擇，需要進行討論。不同的結果進行討論
             //關於全局的時候，這里的這個cl如何算，需要進行討論 
             int LIMIT  =  40  * average;  
             int steal  =  0;  
             for( int k  =  0 ; k  <  256 ; k ++)  
            {  
                 if(tmp2[num][k]  > LIMIT){  
                    steal  += tmp2[num][k]  - LIMIT;  
                    tmp2[num][k]  = LIMIT;  
                }  
            }  
             int bonus  = steal / 256;  
             //hand out the steals averagely  
             for( int k  =  0 ; k  <  256 ; k ++)  
            {  
                tmp2[num][k]  += bonus;  
            }  
             //計算累積分布直方圖  
             for( int k  =  0 ; k  <  256 ; k ++)  
            {  
                 if( k  ==  0)  
                    C2[num][k]  =  1.0f  * tmp2[num][k]  / total;  
                 else  
                    C2[num][k]  = C2[num][k - 1]  +  1.0f  * tmp2[num][k]  / total;  
            }  
        }
    }
     //計算變換后的像素值  
     //根據像素點的位置，選擇不同的計算方法  
     for( int  i  =  0 ; i  < width; i ++)  
    {  
         for( int j  =  0 ; j  < height; j ++)  
        {  
             //four coners  
             if(i  < = width_block / 2  && j  < = height_block / 2)  
            {  
                 int num  =  0;  
                CLAHE_GO.at <uchar >(j,i)  = ( int)(C2[num][CLAHE_GO.at <uchar >(j,i)]  *  255);  
            } else  if(i  < = width_block / 2  && j  > = ((block - 1) *height_block  + height_block / 2)){  
                 int num  = block *(block - 1);  
                CLAHE_GO.at <uchar >(j,i)  = ( int)(C2[num][CLAHE_GO.at <uchar >(j,i)]  *  255);  
            } else  if(i  > = ((block - 1) *width_block +width_block / 2)  && j  < = height_block / 2){  
                 int num  = block - 1;  
                CLAHE_GO.at <uchar >(j,i)  = ( int)(C2[num][CLAHE_GO.at <uchar >(j,i)]  *  255);  
            } else  if(i  > = ((block - 1) *width_block +width_block / 2)  && j  > = ((block - 1) *height_block  + height_block / 2)){  
                 int num  = block *block - 1;  
                CLAHE_GO.at <uchar >(j,i)  = ( int)(C2[num][CLAHE_GO.at <uchar >(j,i)]  *  255);  
            }  
             //four edges except coners  
             else  if( i  < = width_block / 2 )  
            {  
                 //線性插值  
                 int num_i  =  0;  
                 int num_j  = (j  - height_block / 2) /height_block;  
                 int num1  = num_j *block  + num_i;  
                 int num2  = num1  + block;  
                 float p  =  (j  - (num_j *height_block +height_block / 2)) /( 1.0f *height_block);  
                 float q  =  1 -p;  
                CLAHE_GO.at <uchar >(j,i)  = ( int)((q *C2[num1][CLAHE_GO.at <uchar >(j,i)] + p *C2[num2][CLAHE_GO.at <uchar >(j,i)]) *  255);  
            } else  if( i  > = ((block - 1) *width_block +width_block / 2)){  
                 //線性插值  
                 int num_i  = block - 1;  
                 int num_j  = (j  - height_block / 2) /height_block;  
                 int num1  = num_j *block  + num_i;  
                 int num2  = num1  + block;  
                 float p  =  (j  - (num_j *height_block +height_block / 2)) /( 1.0f *height_block);  
                 float q  =  1 -p;  
                CLAHE_GO.at <uchar >(j,i)  = ( int)((q *C2[num1][CLAHE_GO.at <uchar >(j,i)] + p *C2[num2][CLAHE_GO.at <uchar >(j,i)]) *  255);  
            } else  if( j  < = height_block / 2 ){  
                 //線性插值  
                 int num_i  = (i  - width_block / 2) /width_block;  
                 int num_j  =  0;  
                 int num1  = num_j *block  + num_i;  
                 int num2  = num1  +  1;  
                 float p  =  (i  - (num_i *width_block +width_block / 2)) /( 1.0f *width_block);  
                 float q  =  1 -p;  
                CLAHE_GO.at <uchar >(j,i)  = ( int)((q *C2[num1][CLAHE_GO.at <uchar >(j,i)] + p *C2[num2][CLAHE_GO.at <uchar >(j,i)]) *  255);  
            } else  if( j  > = ((block - 1) *height_block  + height_block / 2) ){  
                 //線性插值  
                 int num_i  = (i  - width_block / 2) /width_block;  
                 int num_j  = block - 1;  
                 int num1  = num_j *block  + num_i;  
                 int num2  = num1  +  1;  
                 float p  =  (i  - (num_i *width_block +width_block / 2)) /( 1.0f *width_block);  
                 float q  =  1 -p;  
                CLAHE_GO.at <uchar >(j,i)  = ( int)((q *C2[num1][CLAHE_GO.at <uchar >(j,i)] + p *C2[num2][CLAHE_GO.at <uchar >(j,i)]) *  255);  
            }  
             //雙線性插值
             else{  
                 int num_i  = (i  - width_block / 2) /width_block;  
                 int num_j  = (j  - height_block / 2) /height_block;  
                 int num1  = num_j *block  + num_i;  
                 int num2  = num1  +  1;  
                 int num3  = num1  + block;  
                 int num4  = num2  + block;  
                 float u  = (i  - (num_i *width_block +width_block / 2)) /( 1.0f *width_block);  
                 float v  = (j  - (num_j *height_block +height_block / 2)) /( 1.0f *height_block);  
                CLAHE_GO.at <uchar >(j,i)  = ( int)((u *v *C2[num4][CLAHE_GO.at <uchar >(j,i)]  +   
                    ( 1 -v) *( 1 -u) *C2[num1][CLAHE_GO.at <uchar >(j,i)]  +  
                    u *( 1 -v) *C2[num2][CLAHE_GO.at <uchar >(j,i)]  +  
                    v *( 1 -u) *C2[num3][CLAHE_GO.at <uchar >(j,i)])  *  255);  
            }  
             //最后這步，類似高斯平滑
            CLAHE_GO.at <uchar >(j,i)  = CLAHE_GO.at <uchar >(j,i)  + (CLAHE_GO.at <uchar >(j,i)  <<  8)  + (CLAHE_GO.at <uchar >(j,i)  <<  16);         
        }  
    }  
   return CLAHE_GO;
}

原始圖像

GOCLAHE效果

OpenCV CLAHE效果

從結果上來看，GOCLAHE方法和OpenCV提供的CLAHE方法是一樣的。

再放一組圖片

代碼實現之后，留下兩個問題：

集中在這段代碼

            //這里的參數 對應《Gem》上面 fCliplimit  = 4  , uiNrBins  = 255
             int average  = width_block  * height_block  /  255;  
             //關於參數如何選擇，需要進行討論。不同的結果進行討論
             //關於全局的時候，這里的這個cl如何算，需要進行討論 
             int LIMIT  =  40  * average;  
             int steal  =  0;  

1、在進行CLAHE中CL的計算，也就是限制對比度的計算的時候，參數的選擇缺乏依據。在原始的《GEMS》中提供的參數中，  fCliplimit  = 4  , uiNrBins  = 255. 但是在OpenCV的默認參數中，這里是40.就本例而言，如果從結果上反推，我看10比較好。這里參數的選擇缺乏依據；

2、CLHE是可以用來進行全局直方圖增強的，那么這個時候，這個 average 如何計算，肯定不是width * height/255,這樣就太大了，算出來的LIMIT根本沒有辦法獲得。

但是就實現血管增強的效果而言，這些結果是遠遠不夠的。一般來說，對於CLAHE計算出來的結果，進行Frangi增強或者使用超分辨率增強？結果就是要把血管區域強化出來。

p.s:

arm.jpg 和 hand.jpg