CLAHE的实现和研究

CLAHE算法对于医学图像，特别是医学红外图像的增强效果非常明显。
CLAHE  https://en.wikipedia.org/wiki/Adaptive_histogram_equalization
中文方面非常好的资料 限制对比度自适应直方图均衡化算法原理、实现及效果
在OpenCV中已经实现了CLAHE,但是它在使用过程中，存在参数选择的问题。为了从根本上搞明白，我参考了网络上的一些代码
主要是来源http://blog.csdn.net/abcd1992719g/article/details/25483395
实现了基于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
 
目前方向：图像拼接融合、图像识别 联系方式：jsxyhelu@foxmail.com

CLAHE的实现和研究相关推荐

  1. 循序渐进之（五）空间域图像增强之自适应直方图均衡化（AHE）

    循序渐进之(五)空间域图像增强之自适应直方图均衡化(AHE) 文字摘自:对比度受限的自适应直方图均衡化(CLAHE) 直方图均衡化(HE)是一种很常用的直方图类方法,基本思想是通过图像的灰度分布直方图 ...

  2. 结合实例与代码谈数字图像处理都研究什么？

    图像处理(以及机器视觉)在学校里是一个很大的研究方向,很多研究生.博士生都在导师的带领下从事着这方面的研究.另外,就工作而言,也确实有很多这方面的岗位和机会虚位以待.而且这种情势也越来越凸显.那么图像 ...

  3. 基于直方图的图像增强算法（HE、CLAHE、Retinex）

    直方图是图像色彩统计特征的抽象表述.基于直方图可以实现很多有趣的算法.例如,图像增强中利用直方图来调整图像的对比度.有人利用直方图来进行大规模无损数据隐藏.还有人利用梯度直方图HOG来构建图像特征进而 ...

  4. 综述：视频和图像去雾算法以及相关的图像恢复和增强研究

    综述:视频和图像去雾算法以及相关的图像恢复和增强研究 翻译自IEEE的一篇文章<Review of Video and Image Defogging Algorithms and Relate ...

  5. matlab视网膜血管分割,视网膜血管增强与分割算法研究

    摘要: 视网膜血管图像的研究在图像处理领域占据重要的地位,在临床医学领域更具有重大的实际意义.它不仅能及时地反映出高血压.糖尿病和血液病等全身性疾病,还可以很好地帮助医生掌握病人的情况,因此对眼底视网 ...

  6. 关于 CLAHE 的理解及实现

    CLAHE CLAHE 是一种非常有效的直方图均衡算法, 目前网上已经有很多文章进行了说明, 这里说一下自己的理解. CLAHE是怎么来的 直方图均衡是一种简单快速的图像增强方法, 其原理和实现过程以 ...

  7. OpenCV-Python自适应直方图均衡类CLAHE及方法详解

    一.引言 对比度受限的自适应直方图均衡在OpenCV中是通过类CLAHE来提供实现的,老猿没研究过C++中的应用,但OpenCV-Python中应用时与普通的Python类构建对象的机制有所不同,老猿 ...

  8. OpenCV自适应直方图均衡CLAHE图像和分块大小不能整除的处理

    一.引言 最近一个月来都在研究OpenCV 中CLAHE算法的一些问题,如: 图像横向和纵向分块大小与图像的宽和高不能整除怎么处理? CLIP的剪裁是怎么实施的? 解决棋盘效应的具体插值处理过程怎样? ...

  9. 限制对比度自适应直方图均衡（Contrast Limited Adaptive histgram equalization/CLAHE)

    转自:http://www.cnblogs.com/Imageshop/archive/2013/04/07/3006334.html 一.自适应直方图均衡化(Adaptive histgram eq ...

最新文章

  1. mysql创建表选择字段的时候下尽量小
  2. MySQL修改字段名、字段类型
  3. html鼠标点击切换图片,js鼠标点击图片切换效果代码分享
  4. 微信小程序想要最短服务路径
  5. regardless what you do
  6. Mysql错误：服务名无效。 请键入 NET HELPMSG 2185 以获得更多的帮助。
  7. DataTable to byte[]、DataTable to XML(string)
  8. Xshell分屏显示
  9. 玩转NumPy——NumPy数组的切片和索引
  10. Linux|麒麟操作系统实现多路RTMP|RTSP播放
  11. Spring的数据库编程浅入浅出——不吹牛逼不装逼
  12. 快手视频如何一键批量下载
  13. DSP入门：GPIO
  14. Solidworks常用插件介绍
  15. 北海市卫生学校计算机教室,北海市卫生学校：借力智慧校园，开启德育管理新篇章...
  16. python爬取固定酒店评论_爬取携程上酒店评论数据
  17. windows注册表恢复方法
  18. 4.蒙特卡洛（Monte-Carlo, MC）+时序差分（Temporal Difference, TD）
  19. C语言中的细节知识点（五）
  20. 常见的数据结构和数据库的设计方法

热门文章

  1. NUS CS5477 assignment1
  2. SQL Server学习笔记3——WHERE子句
  3. java-php-python-ssm面相高校学生的图书共享平台计算机毕业设计
  4. tfidf+java+权重,使用scikit-learn tfidf计算词语权重
  5. 2020车载凯立德懒人包下载_华为HarmonyOS App开发工具DevEco Studio下载安装及第一个HarmonyOS App实战教程...
  6. python创建学生字典_用python创建简单字典
  7. RxBus的简单使用(易懂)
  8. 关于xp64位系统须注意的四个问题
  9. 音频设备的接线图了解
  10. Servlet过滤器与SpringMVC拦截器