CLAHE的实现和研究
CLAHE算法对于医学图像，特别是医学红外图像的增强效果非常明显。
中文方面非常好的资料 限制对比度自适应直方图均衡化算法原理、实现及效果
在OpenCV中已经实现了CLAHE,但是它在使用过程中，存在参数选择的问题。为了从根本上搞明白，我参考了网络上的一些代码
实现了基于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
