限制对比度自适应直方图均衡化（自我理解）

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;  
          
    }