CLAHE算法对于医学图像，特别是医学红外图像的增强效果非常明显。
CLAHE https://en.wikipedia.org/wiki/Adaptive_histogram_equalization
中文方面非常好的资料 限制对比度自适应直方图均衡化算法原理、实现及效果（我自己直接看解释没怎么看懂，不如直接看本篇博文下面的代码）
在OpenCV中已经实现了CLAHE,但是它在使用过程中，存在参数选择的问题。为了从根本上搞明白，我参考了网络上的一些代码
主要是来源 http://blog.csdn.net/abcd1992719g/article/details/25483395 （我觉的这篇博客的java代码更容易使人理解）
看完上面的java代码，我的理解：
首先你得了解直方图均衡化（https://blog.csdn.net/u013066730/article/details/82969768）。
限制对比度自适应直方图均衡化方法步骤为：
将整幅图像切分为5*5的块block（但通常代码中为8*8，但本文以5*5举例），如图一右侧的黑色实线框所示。
图一
接着计算每个block的累积分布，在计算累积分布时，限制一下他的对比度，别让直方图出现特别陡峭的情况，直观图像就如图二所示。具体操作就是将直方图中高于某个阈值a的数值全都找出来构成一个集合B，然后B中的每一个元素b都变为a，其中b高于a的部分进行求和得到c，c除以256（该256是指灰度等级，也就是0-255）得到d，将d加在每一个灰度等级下，其实就是图二中蓝色尖状部分变为平坦状。
图二
累积分布求好后，就需要将原始图像中的每个像素值通过累积分布映射为新的像素值。其中使用了双线性插值的方法，而在使用双线性插值时，需要将原图重新分配，也就是图一的粉色，绿色和紫色区域。其中粉色区域的像素值就直接使用对应block的累积分布直接计算即可得到新的像素值（这里如何计算可以参考直方图均衡化）。绿色区域的像素值a，使用绿色区域中邻近的2个block的累计分布计算得到全新的2个像素值b1,b2，然后根据a与这两个block的距离比例得到求和比例c1,1-c1，所以原始像素a对应的最终像素d = b1*c1+b2*(1-c1)。紫色区域的像素值p，使用紫色区域中邻近的4个block的累积分布计算得到全新的4个像素值q1,q2,q3,q4，然后像素p根据与这4个block的距离远近得到对应比例，第一个block的比例为u,v（就是图一中紫色方框叉距离第一个黑色实心框的x，y的比例），第二个block的比例为1-u,v，第三个block的比例为u,1-v，第四个block的比例为1-u,1-v，所以原始像素p对应的最终结果i=u*v*q1+(1-u)*v*q2+u*(1-v)*q3+(1-u)*(1-v)*q4。
哎，我这里说的也有点乱，我感觉不如直接花点时间看下面的代码，真是一看就懂。
/* * CLAHE * 自适应直方图均衡化 */ public int[][] AHE(int[][] oldmat,int pblock) { int block = pblock; //将图像均匀分成等矩形大小，8行8列64个块是常用的选择 int width_block = width/block; int height_block = height/block; //存储各个直方图 int[][] tmp = new int[block*block][256]; //存储累积函数 float[][] C = new float[block*block][256]; //计算累积函数 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; int total = width_block * height_block; for(int ii = start_x ; ii < end_x ; ii++) { for(int jj = start_y ; jj < end_y ; jj++) { int index = oldmat[jj][ii]; tmp[num][index]++; } } //裁剪操作 int average = width_block * height_block / 255; int LIMIT = 4 * average; int steal = 0; for(int k = 0 ; k < 256 ; k++) { if(tmp[num][k] > LIMIT){ steal += tmp[num][k] - LIMIT; tmp[num][k] = LIMIT; } } int bonus = steal/256; //hand out the steals averagely for(int k = 0 ; k < 256 ; k++) { tmp[num][k] += bonus; } //计算累积分布直方图 for(int k = 0 ; k < 256 ; k++) { if( k == 0) C[num][k] = 1.0f * tmp[num][k] / total; else C[num][k] = C[num][k-1] + 1.0f * tmp[num][k] / total; } } } int[][] new_mat = new int[height][width]; //计算变换后的像素值 //根据像素点的位置，选择不同的计算方法 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; new_mat[j][i] = (int)(C[num][oldmat[j][i]] * 255); }else if(i <= width_block/2 && j >= ((block-1)*height_block + height_block/2)){ int num = block*(block-1); new_mat[j][i] = (int)(C[num][oldmat[j][i]] * 255); }else if(i >= ((block-1)*width_block+width_block/2) && j <= height_block/2){ int num = block-1; new_mat[j][i] = (int)(C[num][oldmat[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; new_mat[j][i] = (int)(C[num][oldmat[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; new_mat[j][i] = (int)((q*C[num1][oldmat[j][i]]+ p*C[num2][oldmat[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; new_mat[j][i] = (int)((q*C[num1][oldmat[j][i]]+ p*C[num2][oldmat[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; new_mat[j][i] = (int)((q*C[num1][oldmat[j][i]]+ p*C[num2][oldmat[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; new_mat[j][i] = (int)((q*C[num1][oldmat[j][i]]+ p*C[num2][oldmat[j][i]])* 255); } //inner area 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); new_mat[j][i] = (int)((u*v*C[num4][oldmat[j][i]] + (1-v)*(1-u)*C[num1][oldmat[j][i]] + u*(1-v)*C[num2][oldmat[j][i]] + v*(1-u)*C[num3][oldmat[j][i]]) * 255); } new_mat[j][i] = new_mat[j][i] + (new_mat[j][i] << 8) + (new_mat[j][i] << 16); } } return new_mat; }
