CLAHE算法对于医学图像，特别是医学红外图像的增强效果非常明显。
CLAHE https://en.wikipedia.org/wiki/Adaptive_histogram_equalization
中文方面非常好的资料 限制对比度自适应直方图均衡化算法原理、实现及效果
在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
CLAHE算法对于医学图像，特别是医学红外图像的增强效果非常明显。
CLAHE https://en.wikipedia.org/wiki/Adaptive_histogram_equalization
中文方面非常好的资料 限制对比度自适应直方图均衡化算法原理、实现及效果
在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