CLAHE的实现和研究

CLAHE算法对于医学图像，特别是医学红外图像的增强效果非常明显。

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

中文方面非常好的资料限制对比度自适应直方图均衡化算法原理、实现及效果

在OpenCV中已经实现了CLAHE,但是它在使用过程中，存在参数选择的问题。为了从根本上搞明白，我参考了网络上的一些代码

主要是来源 http://blog.csdn.net/abcdg/article/details/

实现了基于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