图像处理增强算法代码

2017-04-20 15:58:45 charlene_bo 阅读数 22862

关于图像增强必须清楚的基本概念

1.图像增强的目的:

1)改善图像的视觉效果,
2)转换为更适合于人或机器分析处理的形式
3)突出对人或机器分析有意义的信息
4)抑制无用信息,提高图像的使用价值
5)增强后的图像并不一定保真


2,图像增强的方法分类:

1)从处理对象分类:灰度图像,(伪)彩色图像
2)从处理策略分类:全局处理,局部处理(ROI ROI,Region of Interest Interest)
3)从处理方法分类:空间域(点域运算,即灰度变换;邻域方法,即空域滤波),频域方法
4)从处理目的分类:图像锐化,平滑去噪,灰度调整(对比度增强)


3,图像增强的方法之对比度增强

1)灰度变换法

线性变换(已实现)
对数变换(已实现)
指数变换(已实现)

2)直方图调整法
直方图均衡化(已实现)
直方图匹配(未实现)



一,直方图均衡化 

直方图均衡化的英文名称是Histogram Equalization. 

  图像对比度增强的方法可以分成两类:一类是直接对比度增强方法;另一类是间接对比度增强方法。直方图拉伸和直方图均衡化是两种最常见的间接对比度增强方法。直方图拉伸是通过对比度拉伸对直方图进行调整,从而“扩大”前景和背景灰度的差别,以达到增强对比度的目的,这种方法可以利用线性或非线性的方法来实现;直方图均衡化则通过使用累积函数对灰度值进行“调整”以实现对比度的增强。
  直方图均衡化处理的“中心思想”是把原始图像的灰度直方图从比较集中的某个灰度区间变成在全部灰度范围内的均匀分布。直方图均衡化就是对图像进行非线性拉伸,重新分配图像像素值,使一定灰度范围内的像素数量大致相同。直方图均衡化就是把给定图像的直方图分布改变成“均匀”分布直方图分布。
  缺点: 
  1)变换后图像的灰度级减少,某些细节消失; 
  2)某些图像,如直方图有高峰,经处理后对比度不自然的过分增强。 
  直方图均衡化是图像处理领域中利用图像直方图对对比度进行调整的方法。 
  这种方法通常用来增加许多图像的局部对比度,尤其是当图像的有用数据的对比度相当接近的时候。通过这种方法,亮度可以更好地在直方图上分布。这样就可以用于增强局部的对比度而不影响整体的对比度,直方图均衡化通过有效地扩展常用的亮度来实现这种功能。
  这种方法对于背景和前景都太亮或者太暗的图像非常有用,这种方法尤其是可以带来X光图像中更好的骨骼结构显示以及曝光过度或者曝光不足照片中更好的细节。这种方法的一个主要优势是它是一个相当直观的技术并且是可逆操作,如果已知均衡化函数,那么就可以恢复原始的直方图,并且计算量也不大。这种方法的一个缺点是它对处理的数据不加选择,它可能会增加背景杂讯的对比度并且降低有用信号的对比度。
  

关于编程实现,同样是不调用matlab库函数,自己编程实现。这样可以更深刻地理解直方图均衡化技术,提高编程能力。

实现代码(matlab):


  1. clc;  
  2. close all;  
  3. clear all;  
  4.    
  5. src_img = imread('flyman_gray.bmp');    
  6.   
  7. figure (1)   
  8. subplot(321),imshow(src_img),title('原图像');%显示原始图像    
  9. subplot(322),imhist(src_img),title('原图像直方图');%显示原始图像直方图    
  10.   
  11. matlab_eq=histeq(src_img);         %利用matlab的函数直方图均衡化  
  12. subplot(323),imshow(matlab_eq),title('matlab直方图均衡化原图像');%显示原始图像    
  13. subplot(324),imhist(matlab_eq),title('matlab均衡化后的直方图');%显示原始图像直方图   
  14.   
  15. dst_img=myHE(src_img);             %利用自己写的函数直方图均衡化  
  16. subplot(325),imshow(dst_img),title('手写均衡化效果');%显示原始图像    
  17. subplot(326),imhist(dst_img),title('手写均衡化直方图');%显示原始图像直方图   

直方图均衡化函数的实现:


  1. function dst_img=myHE(src_img)    
  2.   
  3. [height,width] = size(src_img);  
  4. dst_img=uint8(zeros(height,width));  
  5. %进行像素灰度统计;      
  6. NumPixel = zeros(1,256);%统计各灰度数目,共256个灰度级      
  7. for i = 1:height      
  8.     for j = 1: width      
  9.         NumPixel(src_img(i,j) + 1) = NumPixel(src_img(i,j) + 1) + 1;%对应灰度值像素点数量增加一      
  10.     end      
  11. end      
  12. %计算灰度分布密度      
  13. ProbPixel = zeros(1,256);      
  14. for i = 1:256      
  15.     ProbPixel(i) = NumPixel(i) / (height * width * 1.0);      
  16. end      
  17. %计算累计直方图分布      
  18. CumuPixel = zeros(1,256);      
  19. for i = 1:256      
  20.     if i == 1      
  21.         CumuPixel(i) = ProbPixel(i);      
  22.     else      
  23.         CumuPixel(i) = CumuPixel(i - 1) + ProbPixel(i);      
  24.     end      
  25. end      
  26.     
  27. % 指定范围进行均衡化    
  28. % pixel_max=max(max(I));    
  29. % pixel_min=min(min(I));    
  30. pixel_max=255;    
  31. pixel_min=0;    
  32. %对灰度值进行映射(均衡化)      
  33. for i = 1:height      
  34.     for j = 1: width      
  35.         dst_img(i,j) = CumuPixel(src_img(i,j)+1)*(pixel_max-pixel_min)+pixel_min;      
  36.     end      
  37. end      
  38. return;  



为什们和matlab的直方图不一样呢???



二,指数变换

指数变换(Power-Law )的公式:S=c*R^r,通过合理的选择c和r可以压缩灰度范围,算法以c=1.0/255.0, r=2实现。
要做该图像增强变换需要先做归一化,再指数变换,最后反归一化
增强效果展示:可以看见,改增强算法并不能很好的将像素尽可能的碾平。
指数增强参考程序为:
  1. clc;  
  2. close all;  
  3. clear all;   
  4.      
  5. % -------------Gamma Transformations-----------------    
  6. %f = imread('Fig0316(4)(bottom_left).tif');     
  7. f = imread('seed.tif');     
  8. Gamma = 0.4;    
  9. g2 = myExpEnhance(f,Gamma);    
  10.   
  11. figure();    
  12. subplot(221);  imshow(f);  xlabel('a).Original Image');    
  13. subplot(222),imhist(f),title('原图像直方图');%显示原始图像直方图    
  14. subplot(223);  imshow(g2);  xlabel('b).Gamma Transformations \gamma = 0.4');    
  15. subplot(224),imhist(g2),title('增强图像直方图');%显示原始图像直方图   
指数增强核心函数为:
  1. function dst_img=myExpEnhance(src_img,Gamma)    
  2. src_img = mat2gray(src_img,[0 255]);%将图像矩阵A中介于amin和amax的数据归一化处理, 其余小于amin的元素都变为0, 大于amax的元素都变为1。    
  3. C = 1;    
  4. g2 = C*(src_img.^Gamma);   
  5. %反归一化  
  6. max=255;  
  7. min=0;  
  8. dst_img=uint8(g2*(max-min)+min);  




三,对数变换

       对数变换主要用于将图像的低灰度值部分扩展,将其高灰度值部分压缩,以达到强调图像低灰度部分的目的。变换方法由下式给出。

这里的对数变换,底数为(v+1),实际计算的时候,需要用换底公式。其输入范围为归一化的【0-1】,其输出也为【0-1】。对于不同的底数,其对应的变换曲线如下图所示。

底数越大,对低灰度部分的强调就越强,对高灰度部分的压缩也就越强。相反的,如果想强调高灰度部分,则用反对数函数就可以了。看下面的实验就可以很直观的理解,下图是某图像的二维傅里叶变换图像,其为了使其灰度部分较为明显,一般都会使用灰度变换处理一下。

效果图:


参考代码:
  1. clc;  
  2. close all;  
  3. clear all;   
  4.   
  5. %-------------Log Transformations-----------------  
  6. f = imread('seed.tif');  
  7.   
  8. g_1 = myLogEnhance(f,10);  
  9. g_2 = myLogEnhance(f,100);  
  10. g_3 = myLogEnhance(f,200);  
  11.   
  12. figure();  
  13. subplot(2,2,1);  
  14. imshow(f);xlabel('a).Original Image');  
  15.   
  16. subplot(2,2,2);  
  17. imshow(g_1);xlabel('b).Log Transformations v=10');  
  18.   
  19. subplot(2,2,3);  
  20. imshow(g_2);xlabel('c).Log Transformations v=100');  
  21.   
  22. subplot(2,2,4);  
  23. imshow(g_3);  
  24. xlabel('d).Log Transformations v=200');  

对数变换核心函数
  1. function dst_img=myLogEnhance(src_img,v)   
  2. c=1.0;  
  3. src_img = mat2gray(src_img,[0 255]);  
  4. g =c*log2(1 + v*src_img)/log2(v+1);  
  5. %反归一化  
  6. max=255;  
  7. min=0;  
  8. dst_img=uint8(g*(max-min)+min);  





四,灰度拉伸

灰度拉伸也用于强调图像的某个部分,与伽马变换与对数变换不同的是,灰度拉升可以改善图像的动态范围。可以将原来低对比度的图像拉伸为高对比度图像。实现灰度拉升的方法很多,其中最简单的一种就是线性拉伸。而这里介绍的方法稍微复杂一些。灰度拉伸所用数学式如下所示。

同样的,其输入r为【0-1】,其输出s也为【0-1】。这个式子再熟悉不过了,跟巴特沃斯高通滤波器像极了,其输入输出关系也大致能猜到是个什么形状的。但是,这里就出现一个问题了,输入为0时候,式子无意义了。所以,在用Matlab计算的时候,将其变为如下形式。

这里的eps,就是Matlab里面,一个很小数。如此做的话,式子变得有意义了。但是,其输入范围为【0-1】的时候,其输出范围变为了。输出范围大致为【0-1】,为了精确起见,使用mat2gray函数将其归一化到精确的[0-1]。调用格式如下。



五,线性拉伸

为了突出感兴趣的目标或者灰度区间,相对抑制那些不感兴趣的灰度区域,可采用分段线性法,常用的是三段线性变换




参考程序:

  1. clc;  
  2. close all;  
  3. clear all;   
  4.   
  5. I=imread('seed.tif');   
  6. [m,n,k]=size(I);  
  7. figure (1)  
  8. imshow('seed.tif');title(' 原图像');   
  9. mid=mean(mean(I));  
  10. %横轴  
  11. fa=20; fb=80;  
  12. %纵轴  
  13. ga=50; gb=230;  
  14.   
  15. J=myLinearEnhance(I,fa,fb,ga,gb);  
  16. figure (2)  
  17. imshow(J);title(' 线性拉伸图像');   
  18.   
  19. pixel_f=1:256;  
  20. pixel_g=zeros(1,256);  
  21.   
  22. %三段斜率,小于1表示该段将会被收缩  
  23. k1=double(ga/fa);   
  24. k2=(gb- ga)/(fb- fa);  
  25. k3=(256- gb)/(256- fb);  
  26. for i=1:256  
  27.     if i <= fa  
  28.         pixel_g(i)= k1*i;  
  29.     elseif fa < i && i <= fb  
  30.         pixel_g(i)= k2*( i- fa)+ ga;  
  31.     else  
  32.         pixel_g(i)= k3*( i - fb)+ gb;  
  33.     end  
  34. end  
  35. figure (3)  
  36. plot(pixel_f,pixel_g);  


核心函数:

  1. function dst_img=myLinearEnhance(src_img,fa,fb,ga,gb)    
  2.   
  3. [height,width] = size(src_img);  
  4. dst_img=uint8(zeros(height,width));  
  5.   
  6. src_img=double(src_img);  
  7.   
  8. %三段斜率  
  9. k1=ga/fa;   
  10. k2=(gb- ga)/(fb- fa);  
  11. k3=(255- gb)/(255- fb);  
  12. for i=1:height  
  13.     for j=1:width  
  14.             if src_img(i,j) <= fa  
  15.                 dst_img(i,j)= k1*src_img(i,j);  
  16.             elseif fa < src_img(i,j) && src_img(i,j) <= fb  
  17.                 dst_img(i,j)= k2*( src_img(i,j)- fa)+ ga;  
  18.             else  
  19.                 dst_img(i,j)= k3*( src_img(i,j)- fb)+ gb;  
  20.             end  
  21.     end  
  22. end  
  23. dst_img=uint8(dst_img);   




附录:

附录网上的另一份讲解:
直方图均衡化算法分为三个步骤,第一步是统计直方图每个灰度级出现的次数,第二步是累计归一化的直方图,第三步是计算新的像素值。
第一步:
for(i=0;i<height;i++)
for(j=0;j<width;j++)
n[s[i][j]]++;

for(i=0;i<L;i++)
p[i]=n[i]/(width*height);

这里,n[i]表示的是灰度级为i的像素的个数,L表示的是最大灰度级,width和height分别表示的是原始图像的宽度和高度,所以,p[i]表示的就是灰度级为i的像素在整幅图像中出现的概率(其实就是p[]这个数组存储的就是这幅图像的归一化之后的直方图)。
第二步:
for(i=0;i<=L;i++)
for(j=0;j<=i;j++)
c[i]+=p[j];

c[]这个数组存储的就是累计的归一化直方图。
第三步:
max=min=s[0][0];
for(i=0;i<height;i++)
for(j=0;j<width;j++)
if(max<s[i][j]){
max=s[i][j];
}else if(min>s[i][j]){
min=s[i][j];
}

找出像素的最大值和最小值。
for(i=0;i<height;i++)
for(j=0;j<width;j++)
t[i][j]=c[s[i][j]]*(max-min)+min;

t[][]就是最终直方图均衡化之后的结果。


收录优秀代码:

这份代码写得不错,学习了,原博客地址见参考资源【3】!

  1. #include <stdio.h>  
  2. #include <iostream>  
  3. #include "fftw3.h"  
  4. #include "string"  
  5. #include "vector"  
  6. #include <windows.h>  
  7. #include <opencv2/legacy/legacy.hpp>  
  8. #include <opencv2/nonfree/nonfree.hpp>//opencv_nonfree模块:包含一些拥有专利的算法,如SIFT、SURF函数源码。   
  9. #include "opencv2/core/core.hpp"  
  10. #include "opencv2/features2d/features2d.hpp"  
  11. #include "opencv2/highgui/highgui.hpp"  
  12. #include <opencv2/nonfree/features2d.hpp>  
  13.   
  14. using namespace cv;  
  15. using namespace std;  
  16.   
  17. class hisEqt  
  18. {  
  19. public:  
  20.     hisEqt::hisEqt();  
  21.     hisEqt::~hisEqt();  
  22. public:  
  23.     int w;  
  24.     int h;  
  25.     int nlen;  
  26.   
  27.     int *pHis;  
  28.     float *pdf;  
  29.   
  30.     //=====求像素分布概率密度====    
  31.     void  getPdf();  
  32.   
  33.     //======统计像素个数=======    
  34.     void getHis(unsigned char*imgdata);  
  35.   
  36.     //==========画统计分布直方图===============    
  37.     void drawHistogram(const float*pdf,Mat &hist1);    
  38.   
  39.     //===========直方图均衡化==========    
  40.     void hisBal();  
  41.   
  42.     //====直方图均衡化后的图像===    
  43.     void imgBal(unsigned char* img);  
  44. };  
  45.   
  46.   
  47. hisEqt::hisEqt() :nlen(0){  
  48.     pHis = new int[256 * sizeof(int)];  
  49.     memset(pHis, 0, 256 * sizeof(int));  
  50.     pdf = new float[255 * sizeof(float)];  
  51.     memset(pdf, 0, 255 * sizeof(float));  
  52. }  
  53.   
  54. hisEqt::~hisEqt(){  
  55.     delete[]pHis;  
  56.     delete[]pdf;  
  57. }  
  58.   
  59.   
  60. //======统计像素个数=======    
  61. void hisEqt::getHis(unsigned char*imgdata){  
  62.     for (int i = 0; i<nlen; i++)  
  63.     {  
  64.         pHis[imgdata[i]]++;  
  65.     }  
  66. }  
  67.   
  68.   
  69. //=====求像素分布概率密度====    
  70. void hisEqt::getPdf(){  
  71.     for (int k = 0; k<256; k++)  
  72.     {  
  73.         pdf[k] = pHis[k] / float(nlen);  
  74.     }  
  75. }  
  76.   
  77. //===========直方图均衡化==========    
  78. void hisEqt::hisBal(){  
  79.     for (int k = 1; k<256; k++)  
  80.     {  
  81.         pdf[k] += pdf[k - 1];  
  82.     }  
  83.     for (int k = 0; k<256; k++)  
  84.     {  
  85.         pHis[k] = 255 * pdf[k];  
  86.     }  
  87. }  
  88.   
  89. //====直方图均衡化    
  90. void hisEqt::imgBal(unsigned char* img){  
  91.     for (int i = 0; i<nlen; i++)  
  92.     {  
  93.         img[i] = pHis[img[i]];  
  94.     }  
  95. }  
  96.   
  97.   
  98. void hisEqt::drawHistogram(const float *pdf, Mat& hist1){  
  99.     for (int k = 0; k<256; k++)  
  100.     {  
  101.         if (k % 2 == 0)  
  102.         {  
  103.             Point a(k, 255), b(k, 255 - pdf[k] * 2550);  
  104.             line(hist1,  
  105.                 a,  
  106.                 b,  
  107.                 Scalar(0, 0, 255),  
  108.                 1);  
  109.         }  
  110.         else  
  111.         {  
  112.             Point a(k, 255), b(k, 255 - pdf[k] * 2550);  
  113.             line(hist1,  
  114.                 a,  
  115.                 b,  
  116.                 Scalar(0, 255, 0),  
  117.                 1);  
  118.         }  
  119.     }  
  120. }  
  121.   
  122.   
  123. int main()  
  124. {  
  125.     Mat image = imread("Fig0651(a)(flower_no_compression).tif");  
  126.     if (!image.data)  
  127.         return -1;  
  128.   
  129.     Mat hist2(256, 256, CV_8UC3, Scalar(0, 0, 0));  
  130.     Mat hist1(256, 256, CV_8UC3, Scalar(0, 0, 0));  
  131.   
  132.     Mat imgOut = Mat(image.rows, image.cols, CV_8UC3, Scalar(0, 0, 0));  
  133.     vector<Mat> planes;  
  134.     int chn = image.channels();  
  135.     if (chn == 3)  
  136.     {  
  137.         split(image, planes);  
  138.     }  
  139.     while (chn)  
  140.     {  
  141.         chn--;  
  142.         unsigned char* imageData = new unsigned char[sizeof(unsigned char)*(image.cols*image.rows)];  
  143.         memcpy(imageData, planes[chn].data, planes[chn].cols*planes[chn].rows);  
  144.         hisEqt his;//自定义的类  
  145.         his.nlen = image.rows*image.cols;  
  146.         his.getHis(imageData);  
  147.         his.getPdf();  
  148.   
  149.         //  //======画原图直方图并保存============    
  150.         his.drawHistogram(his.pdf, hist1);  
  151.         string pic_name = "hisline";  
  152.         pic_name = pic_name + to_string(chn);  
  153.         pic_name=pic_name+  ".jpg";  
  154.         imwrite(pic_name, hist1);  
  155.   
  156.         his.hisBal();  
  157.         his.getPdf();  
  158.         //  //======画均衡化后直方图并保存============    
  159.         his.drawHistogram(his.pdf, hist2);  
  160.         string pic_name0 = "his_balanceline";  
  161.         pic_name0 = pic_name0 + to_string(chn);  
  162.         pic_name0 = pic_name0 + ".jpg";  
  163.         imwrite(pic_name0, hist2);  
  164.   
  165.         //  //=====图像均衡化===    
  166.         his.imgBal(imageData);  
  167.         memcpy(planes[chn].data, imageData, planes[chn].cols*planes[chn].rows);  
  168.         delete[] imageData;  
  169.         imageData = NULL;  
  170.     }  
  171.     merge(planes, imgOut);//单通道合并  
  172.     imwrite("result.jpg", imgOut);  
  173.     return 0;  
  174. }  
2018-07-21 17:12:52 guanzhen3657 阅读数 34690

常见的8种图像增强算法及其opencv实现

1.直方图均衡化

       直方图均衡化是图像处理领域中利用图像直方图对对比度进行调整的方法。 
  这种方法通常用来增加许多图像的局部对比度。这种方法对于背景和前景都太亮或者太暗的图像非常有用,这种方法尤其是可以带来X光图像中更好的骨骼结构显示以及曝光过度或者曝光不足照片中更好的细节。这种方法的一个主要优势是它是一个相当直观的技术并且是可逆操作。

       参考来源 openCV直方图均衡化https://blog.csdn.net/zhangfuliang123/article/details/74170894

      首先openCV没有直方图直接显示的函数,所以我们需要创建直方图来自定义绘图,函数如下:

/*
     width:直方图宽度     height:直方图高度     scale:
*/
IplImage* showImageHistogram(IplImage** image, int width, int height ,int scale){
 
    int dims = 1;
    int histSize = 256;
    float frange[] = { 0, 255 };
    float* ranges[] = { frange };
    //创建一个直方图     CV_HIST_ARRAY多维密集数组
    CvHistogram*  hist = cvCreateHist(dims, &histSize, CV_HIST_ARRAY, ranges);
    //根据输入图像计算直方图
    cvCalcHist(image, hist);
 
 
    //绘制直方图区域
    IplImage* histImage = cvCreateImage(cvSize(width*scale, height), IPL_DEPTH_8U, 1);
    //直方图背景区域置位白色
    cvRectangle(histImage, cvPoint(0, 0), cvPoint(histImage->width, histImage->height), CV_RGB(255,255,255), CV_FILLED);
    //获取最大值
    float maxHistValue = 0;
    cvGetMinMaxHistValue(hist, NULL, &maxHistValue, NULL,NULL);
    //绘制各个灰度级的直方图
    for (int i = 0; i < width; i++){
        float value = cvQueryHistValue_1D(hist, i);
        int drawHeight = cvRound((value / maxHistValue) * height);
        cvRectangle(histImage, cvPoint(i*scale, height-1), cvPoint((i+1)*scale -1, height-drawHeight), cvScalar(i, 0, 0, 0), CV_FILLED);
    }
    return histImage;
 
}


直接处理灰度图:

void histGrayChange(){
    const char* picName = "test.tif";//test.tif lenaRGB.tif
    //采用IplImage *img = cvLoadImage(picName)默认值是CV_LOAD_IMAGE_COLOR  读取无论原始图像的通道数是多少,都将被转换为3个通道读入。
    //IplImage *img = cvLoadImage(picName);
 
    //******以灰度图像读入,强制转换为单通道*****
    IplImage *img = cvLoadImage(picName,CV_LOAD_IMAGE_GRAYSCALE);
    if (img == NULL){
        cout << "Load File Failed." << endl;
    }
    cout << "ChannelL:" << img->nChannels;
 
 
    IplImage* imgDst = cvCreateImage(cvGetSize(img), IPL_DEPTH_8U, 1);
    //直方图均衡化
    cvEqualizeHist(img, imgDst);
    
    cvNamedWindow("Origin", CV_WINDOW_AUTOSIZE);
    cvShowImage("Origin", img);
 
    cvNamedWindow("Result", CV_WINDOW_AUTOSIZE);
    cvShowImage("Result", imgDst);
 
    //
    int histImageWidth = 255;
    int histImageHeight = 150;
    int histImageScale = 2;
    IplImage *histImage1 = showImageHistogram(&img, histImageWidth, histImageHeight, histImageScale);
    cvNamedWindow("Hist1", CV_WINDOW_AUTOSIZE);
    cvShowImage("Hist1", histImage1);
 
    IplImage *histImage2 = showImageHistogram(&imgDst, histImageWidth, histImageHeight, histImageScale);
    cvNamedWindow("Hist2", CV_WINDOW_AUTOSIZE);
    cvShowImage("Hist2", histImage2);
 
 
    cvWaitKey();
 
    cvDestroyWindow("Origin"); cvReleaseImage(&img);
    cvDestroyWindow("Result"); cvReleaseImage(&imgDst);
 
}
 

也可以实现对彩色图片的直方图均衡化,代码就在前面的链接里,不贴了。


 

2、对数图像增强算法

      对数图像增强是图像增强的一种常见方法,其公式为: S = c log(r+1),其中c是常数(以下算法c=255/(log(256)),这样可以实现整个画面的亮度增大。

void LogEnhance(IplImage* img, IplImage* dst)
{
    // 由于oldPixel:[1,256],则可以先保存一个查找表
    uchar lut[256] ={0};
 
    double temp = 255/log(256);
 
    for ( int i =0; i<255; i++)
    {
        lut[i] = (uchar)(temp* log(i+1)+0.5);
    }
 
    for( int row =0; row <img->height; row++)
    {
        uchar *data = (uchar*)img->imageData+ row* img->widthStep;
        uchar *dstData = (uchar*)dst->imageData+ row* dst->widthStep;
 
        for ( int col = 0; col<img->width; col++)
        {
            for( int k=0; k<img->nChannels; k++)
            {
                uchar t1 = data[col*img->nChannels+k];                
                dstData[col*img->nChannels+k] = lut[t1];
            }
        }        
    }    
}

3、指数图像增强算法

 

      指数图像增强的表达为:S = cR^r,通过合理的选择c和r可以压缩灰度范围,算法以c=1.0/255.0, r=2实现。

void ExpEnhance(IplImage* img, IplImage* dst)
{
    // 由于oldPixel:[1,256],则可以先保存一个查找表
    uchar lut[256] ={0};
 
    double temp = 1.0/255.0;
 
    for ( int i =0; i<255; i++)
    {
        lut[i] = (uchar)(temp*i*i+0.5);
    }
 
    for( int row =0; row <img->height; row++)
    {
        uchar *data = (uchar*)img->imageData+ row* img->widthStep;
        uchar *dstData = (uchar*)dst->imageData+ row* dst->widthStep;
 
        for ( int col = 0; col<img->width; col++)
        {
            for( int k=0; k<img->nChannels; k++)
            {
                uchar t1 = data[col*img->nChannels+k];                
                dstData[col*img->nChannels+k] = lut[t1];
            }
        }        
    }    
}
 

4、加Masaic算法(马赛克)

 

        在日常中有时候保密或其他需要将图像马赛克,下面的算法实现图像马赛克功能(原理:用中心像素来表示邻域像素)。

 

uchar getPixel( IplImage* img, int row, int col, int k)
{
    return ((uchar*)img->imageData + row* img->widthStep)[col*img->nChannels +k];
}
 
void setPixel( IplImage* img, int row, int col, int k, uchar val)
{
    ((uchar*)img->imageData + row* img->widthStep)[col*img->nChannels +k] = val;
}

// nSize:为尺寸大小,奇数
// 将邻域的值用中心像素的值替换
void Masic(IplImage* img, IplImage* dst, int nSize)
{
    int offset = (nSize-1)/2;
    for ( int row = offset; row <img->height - offset; row= row+offset)
    {
        for( int col= offset; col<img->width - offset; col = col+offset)
        {
            int val0 = getPixel(img, row, col, 0);
            int val1 = getPixel(img, row, col, 1);
            int val2 = getPixel(img, row, col, 2);
            for ( int m= -offset; m<offset; m++)
            {
                for ( int n=-offset; n<offset; n++)
                {
                    setPixel(dst, row+m, col+n, 0, val0);
                    setPixel(dst, row+m, col+n, 1, val1);
                    setPixel(dst, row+m, col+n, 2, val2);
                }
            }
        }
    }
}
 

5、曝光过度问题处理

 

      对于曝光过度问题,可以通过计算当前图像的反相(255-image),然后取当前图像和反相图像的较小者为当前像素位置的值。

// 过度曝光原理:图像翻转,然后求原图与反图的最小值
void ExporeOver(IplImage* img, IplImage* dst)
{
    for( int row =0; row <img->height; row++)
    {
        uchar *data = (uchar*)img->imageData+ row* img->widthStep;
        uchar *dstData = (uchar*)dst->imageData+ row* dst->widthStep;
 
        for ( int col = 0; col<img->width; col++)
        {
            for( int k=0; k<img->nChannels; k++)
            {
                uchar t1 = data[col*img->nChannels+k];
                uchar t2 = 255 - t1;
                dstData[col*img->nChannels+k] = min(t1,t2);
            }
        }        
    }
}
 

6、高反差保留

      高反差保留主要是将图像中颜色、明暗反差较大两部分的交界处保留下来,比如图像中有一个人和一块石头,那么石头的轮廓线和人的轮廓线以及面部、服装等有明显线条的地方会变被保留,儿其他大面积无明显明暗变化的地方则生成中灰色。其表达形式为:dst = r*(img - Blur(img))。

 

Mat HighPass(Mat img)
{
    Mat temp;
    GaussianBlur(img, temp,Size(7,7),1.6,1.6);
 
    int r=3;    
    Mat diff = img + r*(img-temp); //高反差保留算法
    return diff;
}

测试代码:

int main(int argc, char* argv[])
{
    const char* Path = "02.bmp";
    IplImage *img = cvLoadImage(Path,CV_LOAD_IMAGE_ANYCOLOR);
    IplImage *dst = cvCreateImage(cvGetSize(img), img->depth, img->nChannels);
    cout<<"输入你要选择的操作:"<<endl;
    cout<<"1、曝光过度"<<endl;
    cout<<"2、加马赛克"<<endl;
    cout<<"3、对数增强"<<endl;
    cout<<"4、指数增强"<<endl;
    cout<<"请输入你的选择:";
    int choice = 1;
    cin>>choice;
    switch (choice)
    {
    case 1: 
        ExporeOver(img, dst);   //这四个算法中总觉得某个算法有问题
        break;
    case 2: 
        Masic(img, dst, 21);
        break;
    case 3: 
        LogEnhance(img, dst);
        break;
    case 4:
        ExpEnhance(img, dst);
        break;
    default:
        cout<<"输入错误"<<endl;
        break;              
    }
    cvSaveImage("dst.jpg",dst);
    cvNamedWindow("SRC",1);
    cvNamedWindow("DST", 1);
    cvShowImage("SRC", img);
    cvShowImage("DST", dst);
    cvWaitKey();
    return 0;
}

7.拉普拉斯算子图像增强

使用中心为5的8邻域拉普拉斯算子与图像卷积可以达到锐化增强图像的目的。

拉普拉斯算子增强局部的图像对比度的opencv代码:

#include <opencv2/highgui/highgui.hpp>
#include <opencv2/imgproc/imgproc.hpp>
#include <iostream>
 
using namespace cv;
 
int main(int argc, char *argv[])
{
    Mat image = imread("/Users/shandiangou/Downloads/lena.png");
    if (image.empty())
    {
        std::cout << "打开图片失败,请检查" << std::endl;
        return -1;
    }
    imshow("原图像", image);
    waitKey();
    Mat imageEnhance;
    Mat kernel = (Mat_<float>(3, 3) << 0, -1, 0, 0, 5, 0, 0, -1, 0);
    filter2D(image, imageEnhance, CV_8UC3, kernel);
    imshow("拉普拉斯算子图像增强效果", imageEnhance);
    waitKey();
    return 0;
}
 

参考博客:直方图均衡化、拉普拉斯算子图像增强、Gamma校正https://blog.csdn.net/sinat_28296297/article/details/77972023

                   OpenCV_基于Laplacian算子的图像边缘增强https://blog.csdn.net/icvpr/article/details/8502949

8.Gamma校正

       伽马变换主要用于图像的校正,将灰度过高或者灰度过低的图片进行修正,增强对比度。伽马变换对图像的修正作用其实就是通过增强低灰度或高灰度的细节实现的。

实现代码

#include <opencv2/core/core.hpp>
#include <opencv2/imgproc/imgproc.hpp>
#include <opencv2/highgui/highgui.hpp>
 
#include <iostream>
 
using namespace cv;
using namespace std;
 
// Normalizes a given image into a value range between 0 and 255.
Mat norm_0_255(const Mat& src) {
    // Create and return normalized image:
    Mat dst;
    switch(src.channels()) {
        case 1:
            cv::normalize(src, dst, 0, 255, NORM_MINMAX, CV_8UC1);
            break;
        case 3:
            cv::normalize(src, dst, 0, 255, NORM_MINMAX, CV_8UC3);
            break;
        default:
            src.copyTo(dst);
            break;
    }
    return dst;
}
 
int main(int argc, const char *argv[]) {
    // Get filename to the source image:
    // Load image & get skin proportions:
    Mat image = imread("/Users/shandiangou/Downloads/guobao.jpeg");
    // Convert to floating point:
    Mat X;
    image.convertTo(X, CV_32FC1);
    //image.convertTo(X, CV_32F);
    // Start preprocessing:
    Mat I;
    float gamma = 3;
    pow(X, gamma, I);
    
    
    // Draw it on screen:
    imshow("Original Image", image);
    waitKey();
    imshow("Gamma correction image", norm_0_255(I));
    // Show the images:
    waitKey();
    // Success!
    return 0;  
}

暂时总结的就这8类,其他的之后再添

2016-01-07 17:20:43 EbowTang 阅读数 19513

      小波变换下的图像对比度增强技术实质上是通过小波变换把图像信号分解成不同子带,针对不同子带应用不同的算法来增强不同频率范围内的图像分量,突出不同尺度下的近似和细节,从而达到增强图像层次感的目的。


       根据小波的多分辨率分析原理将图像进行多级二维离散小波变换,可以将图像分解成图像近似信号的低频子带和图像细节信号的高频子带。其中,图像中大部分的噪声和一些边缘细节都属于高频子带,而低频子带主要表征图像的近似信号。为了能够在增强图像的同时减少噪声的影响,可以对低频子带进行非线性图像增强,用以增强目标的对比度,抑制背景;而对高频部分进行小波去噪处理,减少噪声对图像的影响。最后小波重构得到增强的红外图像。基于小波变换的红外图像增强模型如下图所示




1,基于小波变换的简单线性变换对比度增强算法:

load woman;  
subplot(121);
image(X);colormap(map);title('原始图像');%画出原图像

[c,s]=wavedec2(X,3,'sym4');  %进行二层小波分解

len=length(c);
justdet = prod(s(1,:));%截取细节系数起始位置(不处理近似系数)  
%处理低频分解系数,突出轮廓
for I =1:justdet
    if(c( I )>250)
      c( I )=1.5*c( I );
    end
end
%处理高频分解系数,弱化细节
for I =justdet:len
    if(c( I ) < 150)
      c( I )=0.75*c( I );
    end
end

nx=waverec2(c,s,'sym4');%分解系数重构

subplot(122);
image(nx);title('增强图像')%画出增强图像


利用上面代码进行两层分解得到的增强效果图:



利用上面代码进行三层分解得到的增强效果图:



通过对小波系数进行简单的线性变来增强图像对比度,虽然能够取得一定的效果,根据对低频和高频系数增加不同的倍数,可以得到不同的视觉效果,但倍数的选择更多受到不同人的主观影响,不能做到最优的选取。

另外,简单的使用该方法还是会在增强图像的同时加大噪声,影响对比度的增强效果。



2,小波变换和直方图均衡化结合的图像增强算法

通俗一点说,直方图均衡化就是对图像进行非线性拉伸,重新分配图像像素值,使一定灰度范围内的像素数量大致相同。直方图均衡化就是把给定图像的直方图分布改变成“均匀”分布直方图分布。一般而言,经过直方图均衡化处理之后,灰度值出现概率较小的像素会被合并,从而导致该图像的部分灰度值被压缩,而灰度值出现概率较大的像素则会被拉伸。因此灰度值的合并必然引起来图像细节的丢失,而灰度值拉伸过程中必然会使得图像中的噪声被局部的增强。尤其是对于红外图像中背景和目标的灰度值非常接近的情况下,必然使得红外图像中的噪声被放大,细节丢失。

图像中经过小波变换后得到小波系数中的低频信息实际上表征了图像的低频信息,也即图像的整体轮廓。

而图像中的高频小波系数则包含了图像的细节信息,包括噪声、细节和边缘信息等。

因此通过只对低频小波系数进行直方图均衡化的方法就会很好的改善图像的噪声信息被放大的特性,因为图像的噪声信息已经大部分被包含在高频小波系数中。而对于高频小波系数而言,可以通过类似硬阈值去噪的方法,去除图像中的噪声,并对非噪声的高频小波系数进行增强,从而强化图像中的细节部分。这样就用小波变换的方法很好的弥补了直方图均衡化算法中增强噪声和丢失细节的问题,从而更好的完成了图像的增强。


I=imread('oct.bmp');
figure 
imshow(I);title('yuan')


[ca,ch,cv,cd]=dwt2(I,'db6');
X=imadjust(uint8(ca),[],[],1.3);
xx=idwt2(uint8(X),ch,cv,cd,'db6');
figure 
imshow(uint8(xx));



参考资源:

【1】一种基于小波变换的非线性红外图像增强算法刘兴淼,王仕成,赵 静(陕西省第二炮兵工程学院,陕西 西安 710025)


2018-11-28 22:50:27 entoon 阅读数 1712

在实际应用当中,有时候需要进行图像增强来改善图像的视觉效果。在此问题处理当中,按照颜色可以分为灰度图像增强和彩色图像增强。按照作用域分类,可以分为空域处理和频域处理。

图像空域处理方法通常有灰度变换,直方图均衡,图像平滑和锐化。频域处理有DFT变换,采用滤波的方法进行图像增强。现有的方法自适应的效果都比较差,这里提出一种模糊自适应的方法,利用遗传算法完成图像的增强。

1,选着合适的评价函数,得到图像质量的适应度函数值

2,利用遗传算法,选择,交叉,变异 优化处理

3,得到优化结果,完成图像的自适应增强。
本案例采用beta函数来实现图像灰度变换曲线的自动拟合。相关代码如下
function Incmp_Beta_Result=IncmpBeta(a,b,X)
[m n]=size(X);
aaa=Gammln(a+b)-Gammln(a)-Gammln(b);
for i=1:m
for j=1:n
if X(i,j)<0 | X(i,j)>1
helpdlg(‘变量X的取值范围不在0和1之间’,‘错误!’);
return
end
if X(i,j)==0 | X(i,j)==1
bt(i,j)=0;
else
bt(i,j)=exp(aaa+alog(X(i,j))+blog(1-X(i,j)));
end
end
end

for i=1:m
for j=1:n
if X(i,j)<((a+1)/(a+b+2))
Incmp_Beta_Result(i,j)=bt(i,j)*Betacf(a,b,X(i,j))/a;
else
Incmp_Beta_Result(i,j)=1-bt(i,j)*Betacf(b,a,1-X(i,j))/b;
end
end
end

%//////////////Betacf过程//////////////////////////////
function Betacf_Result=Betacf(a,b,x)
itmax=100;
eps=0.0000003;
am=1;
bm=1;
az=1;
qab=a+b;
qap=a+1;
qam=a-1;
bz=1-qabx/qap;
for m=1:itmax
em=m;
tem=em+em;
d=em
(b-m)x/((qam+tem)(a+tem));
ap=az+dam;
bp=bz+d
bm;
d=-(a+em)(qab+em)x/((a+tem)(qap+tem));
aap=ap+d
az;
bpp=bp+dbz;
aold=az;
am=ap/bpp;
bm=bp/bpp;
az=aap/bpp;
bz=1;
if abs(az-aold)<eps
abs(az)
Betacf_Result=az;
break;
end
end
Betacf_Result=az;

%/////////////计算gama函数程序////////////////////////
function Gammln_result=Gammln(xx)
cof(1)=76.18009173;
cof(2)=-86.50532033;
cof(3)=24.01409822;
cof(4)=-1.231739516;
cof(5)=0.00120858003;
cof(6)=-0.00000536382;
stp=2.50662827465;
half=0.5;
one=1.0;
fpf=5.5;
x=xx-one;
tmp=x+fpf;
tmp=(x+half)log(tmp)-tmp;
ser=one;
for j=1:6
x=x+one;
ser=ser+cof(j)./x;
end
Gammln_result=tmp+log(stp
ser);
%----------programs end here---------------

适应度函数
function [sol eval]=fitnesszq_gray(sol,options)
%实现灰度图像增强的适应度函数
%%遗传算法的适应度函数。公式3-23的编程实现。
%////////////////////////////////////////////////////////////////////
a=sol(1);
b=sol(2);
P=imread(‘my.png’);
[M,N]=size§;
PP=double§;
n=M*N;
%求beta函数
syms t x
B=double(int((t.(a-1)).*((1-t).(b-1)),t,0,1));%分母
lmin=double(min(min§));
lmax=double(max(max§));
%原图像 P 归一化处理得到 G
G=(PP-lmin)/(lmax-lmin);
GP=IncmpBeta(a,b,G);
FP=round((lmax-lmin).GP+lmin);%原图像的变换
count=imhist(uint8(FP));
th=5;
NN=sum(count>th);%求像素个数大于一给定阈值的灰度级的数量
H=0;
co=count/(M
N);
%cc=co>0;
for i=1:256
if co(i)>0
H=H+co(i).*log2(co(i));
end
end
%H=-sum(co.*log2(co));
HH=-H;%求熵
delta=sum(sum(FP.2))/n-(sum(sum(FP))/n).2;%求方差
U=FP./(lmax-1);%模糊变换函数
pmr1=0;pmr2=0;
for i=1:M
for j=1:N-1
pmr1=pmr1+abs(U(i,j)-U(i,j+1));
end
end
for i=1:M-1
for j=1:N
pmr2=pmr2+abs(U(i,j)-U(i+1,j));
end
end
pmr=pmr1+pmr2;
area=sum(sum(U));%模糊几何量
Comp=area./(pmr.^2);
%本节设计的适应度函数
Fitness=log10(NN.*HH.*delta/Comp) %没有加;是为了在运行窗口可以看到运行中的结果

eval=Fitness;
%----------all programs end here---------------
%%%%来自私人定制程序文案工作室%%%%
https://weidian.com/?userid=1808020072&wfr=wx&sfr=app&source=shop&from=singlemessage

2016-02-18 17:35:36 sky200543012 阅读数 25864

转自:http://www.cnblogs.com/sleepwalker/p/3676600.html?utm_source=tuicool&utm_medium=referral

Retinex图像增强算法

前一段时间研究了一下图像增强算法,发现Retinex理论在彩色图像增强、图像去雾、彩色图像恢复方面拥有很好的效果,下面介绍一下我对该算法的理解。

Retinex理论

Retinex理论始于Land和McCann于20世纪60年代作出的一系列贡献,其基本思想是人感知到某点的颜色和亮度并不仅仅取决于该点进入人眼的绝对光线,还和其周围的颜色和亮度有关。Retinex这个词是由视网膜(Retina)和大脑皮层(Cortex)两个词组合构成的.Land之所以设计这个词,是为了表明他不清楚视觉系统的特性究竟取决于此两个生理结构中的哪一个,抑或是与两者都有关系。

Land的Retinex模型是建立在以下的基础之上的:

一、真实世界是无颜色的,我们所感知的颜色是光与物质的相互作用的结果。我们见到的水是无色的,但是水膜—肥皂膜却是显现五彩缤纷,那是薄膜表面光干涉的结果;

二、每一颜色区域由给定波长的红、绿、蓝三原色构成的;

三、三原色决定了每个单位区域的颜色。

Retinex 理论的基本内容是物体的颜色是由物体对长波(红)、中波(绿)和短波(蓝)光线的反射能力决定的,而不是由反射光强度的绝对值决定的;物体的色彩不受光照非均性的影响,具有一致性,即Retinex理论是以色感一致性(颜色恒常性)为基础的。如下图所示,观察者所看到的物体的图像S是由物体表面对入射光L反射得到的,反射率R由物体本身决定,不受入射光L变化。

clip_image002[4]

图1.Retinex理论中图像的构成

Retinex理论的基本假设是原始图像S是光照图像L和反射率图像R的乘积,即可表示为下式的形式:

clip_image004[4]

基于Retinex的图像增强的目的就是从原始图像S中估计出光照L,从而分解出R,消除光照不均的影响,以改善图像的视觉效果,正如人类视觉系统那样。在处理中,通常将图像转至对数域,即clip_image006[4],从而将乘积关系转换为和的关系:

clip_image008[4]

Retinex方法的核心就是估测照度L,从图像S中估测L分量,并去除L分量,得到原始反射分量R,即:

clip_image010[4]

函数clip_image012[4]实现对照度L的估计(可以去这么理解,实际很多都是直接估计r分量)。

Retinex理论的理解

如果大家看论文,那么在接下去的篇幅当中,肯定会介绍两个经典的Retinex算法:基于路径的Retinex以及基于中心/环绕Retinex。在介绍两个经典的Retinex算法之前,我先来讲一点个人的理解,以便第一次接触该理论的朋友能够更快速地理解。当然,如果我的理解有问题,也请大家帮忙指出。

Retinex理论就我理解,与降噪类似,该理论的关键就是合理地假设了图像的构成。如果将观察者看到的图像看成是一幅带有乘性噪声的图像,那么入射光的分量就是一种乘性的,相对均匀,且变换缓慢的噪声。Retinex算法所做的就是合理地估计图像中各个位置的噪声,并除去它。

在极端情况下,我们大可以认为整幅图像中的分量都是均匀的,那么最简单的估计照度L的方式就是在将图像变换到对数域后对整幅图像求均值。因此,我设计了以下算法来验证自己的猜想,流程如下:

(1) 将图像变换到对数域clip_image014[4] ;

(2) 归一化去除加性分量clip_image016[4] ;

(3) 对步骤3得到的结果求指数,反变换到实数域clip_image018[4]

这里为了简化描述,省略了对图像本身格式的变换,算法用Matlab实现:

复制代码
% ImOriginal:原始图像
% type:'add'表示分量是加性的,如雾天图像;'mult'表示分量是乘性的,如对照度的估计
[m,n,z] = size(ImOriginal);
ImOut = uint8(zeros(m,n,z));
for i = 1:z
    if strcmp(type,'add')
        ImChannel = double(ImOriginal(:,:,i))+eps;
    elseif strcmp(type,'mult')
        ImChannel = log(double(ImOriginal(:,:,i))+eps);
    else
        error('type must be ''add'' or ''mult''');
    end
    ImOut(:,:,i) = EnhanceOneChannel(ImChannel);
end
ImOut = max(min(ImOut,255), 0);
end

function ImOut = EnhanceOneChannel(ImChannel)
% 计算计算单个通道的反射分量
% 1.对全图进行照射分量估计
% 2.减去照射分量
% 3.灰度拉伸
ImChannel = ImChannel./max(ImChannel(:));
ImRetinex = round(exp(ImChannel.*5.54));
ImOut = uint8(ImRetinex);
end
复制代码

为了验证算法的有效性,这里使用经典的Retinex算法与我所用的算法进行对比试验,效果如下:

clip_image002[6]

图2.测试原图

clip_image004[6]

图3.经典Retinex算法结果

clip_image006[6]

图4.上述方法结果

从对比中可以看到,对于去除照度,还原图像本身来讲,效果还可以,并且不会在边缘位置产生光晕现象。缺点就是在去除照度分量L过程中,保留的反射分量R我在上述算法中使用归一化后直接进行反变换。这一步的作用可以近似看成去除一个均匀的直流分量,即均匀的照度分量。由于操作都是全局的,这里默认假设了所有位置的照射分量都是相同的,因此在灰度拉伸的时候没有照顾到局部的特性,图像整体亮度偏暗。当然,全局的照度估计对于图像的增强肯定有相当的局限性,其增强效果在色彩的还原和亮度处理等方面还是有一定缺陷的。

个人认为,Retinex算法的关键还是正确的分析了噪声的性质。相信很多人都看到利用基于Retinex的图像水下增强、基于Retinex的图像去雾等等,我也好奇,那就试试吧。大雾图片嘛谁没有,前一阵子大雾天,没事拍了几张照片,终于用上了,请看对比图:

clip_image008[6]

图5.有雾原图

clip_image010[6]

图6.经典Retinex去雾效果

clip_image012[6]

图7.上述方法去雾效果

还是老规矩,这时候对比试验还是最能说明效果的,为此选了一幅干扰很大的图像。基本上各位要是显示器比较差一点,从原图当中是很难看出大雾后面的东西。从去雾效果来看,上述方法的效果并不比经典算法要差,至少在去雾的效果上,本实验结果从主观上给人的感觉还是不错的。

在上述案例中,最重要的就是正确分析有雾图像的结构,与Retinex理论一开始的核心思想有区别的是,在针对这种加性的干扰时,经典的Retinex算法在处理过程中,其实仅仅是利用其估计加性的干扰分量;当然,抛开Retinex理论对照度、反射率对最终图像形成的核心思想(如图1),后续最重要的就是对这个加性的干扰的估计了。

对于有雾的图像,我们大可以看作透过一块磨砂玻璃去看一幅清晰的图像,这样大家就能很好理解为什么认为在这个案例中,将雾的干扰认为是一个加性的了。诸如后面两个经典的算法,所有的这类算法归根结底就是更好地利用原图像中的像素点去估计原始照度。从上面例程上可以看出,使用一个全局估计对局部的增强是比较差的,如果存在照度不均匀(雾的浓度不均匀),或者背景颜色亮度很高等情况时,处理结果会趋向恶劣,效果比较差。

当然,经典也不是完美,从图3中可以看到,经典的算法容易出现光晕效果(蓝色书本文字周围一圈白色),因此后续对于照度估计和去除光晕等问题又有很多基于Retinex理论的变体算法,这里暂不进行介绍。下面开始介绍两个经典的算法,查看Matlab代码下载点击:

http://www.cs.sfu.ca/~colour/publications/IST-2000/

McCann Retinex算法

McCann Retinex 算法是McCann和Frankle一起提出的一种Retinex算法,该算法是一种基于多重迭代策略的Retinex算法,单个点的像素值取决于一条特定路径的环绕的结果,经过多次迭代逼近理想值。通过比较螺旋式路径上的各像素点的灰度值来估计和去除图像的照度分量。

clip_image002[8]

图8.McCann Retinex算法路径选择示意图

这个图是参照人家论文的,不过我准备像论文那样讲,因为太复杂了,不容易懂。从图中我们可以看到,算法沿着这个螺旋式的路径选取用于估计的点,并且越靠近预测的中心点选取的点数越多。这个从实际物理上也解释的通,靠的近的像素点与中心像素点的相关性肯定要要比远处的点要高。

该算法估测图像经过以下几个步骤:

1. 将原图像变换到对数域clip_image004[8],彩色图像对各通道都进行对数变换;

2. 初始化常数图像矩阵clip_image006[8],该矩阵作为进行迭代运算的初始值;

3. 计算路径,如上图8所示,这里令clip_image008[8]为路径上的点,从远到近排列;

4. 对路径上的像素点按照如下公式运算:

clip_image010[8]

clip_image012[8]

公式所表示的大致意思为:从远到近,中心点像素值减去路径上的像素值得到的差值的一半与前一时刻的估计值之间的和。最终,中心像素点的像素大致的形式为

clip_image014[6]

其中clip_image016[6]表示中心位置最终的反射率估计,clip_image018[6]为常数值为转换后的图像中的最大值,在步骤2中被确定。从这里将clip_image020按照Retinex理论进行分解,最终公式中可以看到,最终照度分量被去除了,而中心位置的反射率由路径上各点的反射率之间的差值进行估计,并且从轨迹上可以看到,靠的越近的在最终估计的时候所占比重越大。

就我个人理解,这类估计算法,实质是将中心位置和周围亮度之间的差异作为最终中心位置的反射率的估计。如果中心位置亮度本身高,那么最终的结果还是高亮度的;如果中心位置亮度低,那么中心与其它点的差肯定是负值,最终的结果clip_image016[7]就比较小,亮度就比较低。这就要求路径上的点需要能够很好地代表整幅图像的特性,如果图像本身很规则,那么可能中央与周围计算的结果无法估测该点像素本身的灰度特性,结果就与预期的可能不一样了。如下图9可以看到,算法实质是估计中心和周边的差值,因此图中原本黑色区域的图像由于与周边差别很小,因此呈现高亮度;而在小型的矩形周围随着距离的增大,差别渐渐变小,因此亮度逐渐升高,无法体现原本黑色像素点处原本的亮度特性。

clip_image022

图9.算法的测试图(来自Retinex in Matlab)

从原始算法中可以看到,还有一个重要的步骤,就是迭代。迭代的作用是尽可能保留中央周边差的分量,原本每次只保留中央周边差的一半(见步骤4中最后的除2的处理),迭代次数越多,保留的分量就越多,迭代n次保留的分量就是clip_image024,这样局部的信息就更多,相当于降低动态范围的压缩。这样的操作可以使图像更加自然,但是会增加运算量,下图可以更加明显地看出迭代次数的影响:

clip_image026

图10.迭代次数的影响

为了方便各位自己研究,下面我给出该算法源码供大家参考:

复制代码
function Retinex = retinex_frankle_mccann(L, nIterations)

% RETINEX_FRANKLE_McCANN: 
%         Computes the raw Retinex output from an intensity image, based on the
%         original model described in:
%         Frankle, J. and McCann, J., "Method and Apparatus for Lightness Imaging"
%         US Patent #4,384,336, May 17, 1983
%
% INPUT:  L           - logarithmic single-channel intensity image to be processed
%         nIterations - number of Retinex iterations
%
% OUTPUT: Retinex     - raw Retinex output
%
% NOTES:  - The input image is assumed to be logarithmic and in the range [0..1]
%         - To obtain the retinex "sensation" prediction, a look-up-table needs to
%         be applied to the raw retinex output
%         - For colour images, apply the algorithm individually for each channel
%
% AUTHORS: Florian Ciurea, Brian Funt and John McCann. 
%          Code developed at Simon Fraser University.
%
% For information about the code see: Brian Funt, Florian Ciurea, and John McCann
% "Retinex in Matlab," by Proceedings of the IS&T/SID Eighth Color Imaging 
% Conference: Color Science, Systems and Applications, 2000, pp 112-121.
%
% paper available online at http://www.cs.sfu.ca/~colour/publications/IST-2000/
%
% Copyright 2000. Permission granted to use and copy the code for research and 
% educational purposes only.  Sale of the code is not permitted. The code may be 
% redistributed so long as the original source and authors are cited.

global RR IP OP NP Maximum
RR = L;
Maximum = max(L(:));                                 % maximum color value in the image
[nrows, ncols] = size(L);

shift = 2^(fix(log2(min(nrows, ncols)))-1);          % initial shift
OP = Maximum*ones(nrows, ncols);                     % initialize Old Product

while (abs(shift) >= 1)
   for i = 1:nIterations
      CompareWith(0, shift);                         % horizontal step
      CompareWith(shift, 0);                         % vertical step
   end
   shift = -shift/2;                                 % update the shift
end
Retinex = NP;

function CompareWith(s_row, s_col)
global RR IP OP NP Maximum
IP = OP;
if (s_row + s_col > 0)
   IP((s_row+1):end, (s_col+1):end) = OP(1:(end-s_row), 1:(end-s_col)) + ...
   RR((s_row+1):end, (s_col+1):end) - RR(1:(end-s_row), 1:(end-s_col));
else
   IP(1:(end+s_row), 1:(end+s_col)) = OP((1-s_row):end, (1-s_col):end) + ...
   RR(1:(end+s_row),1:(end+s_col)) - RR((1-s_row):end, (1-s_col):end);
end
IP(IP > Maximum) = Maximum;                          % The Reset operation
NP = (IP + OP)/2;                                    % average with the previous Old Product
OP = NP;                                             % get ready for the next comparison
复制代码

McCann99 Retinex算法

McCann99 Retinex算法本质上与McCann Retinex算法没有区别,两者不同的就是在McCann99算法中,不再是使用螺旋式的路径来选取估计像素点,而是使用图像金字塔的方式逐层选取像素。算法同样经取点、比较、平均这几步进行迭代运算。

图像金字塔的概念很简单,就是对图像进行下采样,以多分辨率的形式表示图像。最顶层的图像分辨率最低,底层的最高(一般为原图)。

clip_image001

图11.图像金字塔示意图(来自网络)

如上图所示,McCann99算法就是从顶层开始让每一个像素点与其8个相邻的像素进行比较,估计反射率分量clip_image003;前一层计算结束之后对估计的反射率分类进行插值运算,使上一层估计的结果clip_image003[1]的图像与金字塔下一层图像尺寸相同,并再一次进行相同的比较运算;最终,对原始图像进行8邻域的比较结束后就能得到最终的结果,即增强后的图像了。其中,比较操作与McCann算法相同,都是相减后求平均,见上一节步骤4中描述。

因此,McCann99算法,此处简化描述为以下几步:

1. 将原图像变换到对数域clip_image005,彩色图像对各通道都进行对数变换;

2. 初始化(计算图像金字塔层数;初始化常数图像矩阵clip_image007,该矩阵作为进行迭代运算的初始值);

3. 从顶层开始,到最后一层进行8邻域比较运算,运算规则与MccCann Retinex算法相同,见上一节步骤4;

4. 第n层运算结束后对第n层的运算结果clip_image009进行插值,变成原来的两倍,与n+1层大小相同(此处默认n越大越靠近底层);

5. 当最底层计算完毕得到的clip_image011即最终增强后的图像。

为了方便各位自己研究,下面我给出该算法源码供大家参考:

复制代码

function Retinex = retinex_mccann99(L, nIterations)


global OPE RRE Maximum
[nrows ncols] = size(L)
;                             % get size of the input image
nLayers = ComputeLayers(nrows, ncols);               % compute the number of pyramid layers
nrows = nrows/(2^nLayers);                           % size of image to process for layer 0
ncols = ncols/(2^nLayers);
if (nrows*ncols > 25)                                % not processing images of area > 25
  error(
'invalid image size.')                       % at first layer
end
Maximum = max(L(:))
;                                 % maximum color value in the image
OP = Maximum*ones([nrows ncols]);                    % initialize Old Product
for layer = 0:nLayers
   RR = ImageDownResolution(L, 
2^(nLayers-layer));   % reduce input to required layer size
   
   OPE = [zeros(nrows,
1) OP zeros(nrows,1)];         % pad OP with additional columns
   OPE = [zeros(1,ncols+2); OPE; zeros(1,ncols+2)];  % and rows
   RRE = [RR(:,1) RR RR(:,end)];                     % pad RR with additional columns
   RRE = [RRE(1,:); RRE; RRE(end,:)];                % and rows
   
   for iter = 
1:nIterations
     CompareWithNeighbor(-
10);                     % North
     CompareWithNeighbor(-11);                     % North-East
     CompareWithNeighbor(01);                      % East
     CompareWithNeighbor(11);                      % South-East
     CompareWithNeighbor(10);                      % South
     CompareWithNeighbor(1, -1);                     % South-West
     CompareWithNeighbor(0, -1);                     % West
     CompareWithNeighbor(-1, -1);                    % North-West
   end
   NP = OPE(
2:(end-1), 2:(end-1));
   OP = NP(:, [fix(1:0.5:ncols) ncols]);             %%% these two lines are equivalent with 
   OP = OP([fix(1:0.5:nrows) nrows], :);             %%% OP = imresize(NP, 2) if using Image
   nrows = 2*nrows; ncols = 2*ncols;                 % Processing Toolbox in MATLAB
end
Retinex = NP
;

function CompareWithNeighbor(dif_row, dif_col)
global OPE RRE Maximum

% Ratio-Product operation
IP = OPE(
2+dif_row:(end-1+dif_row), 2+dif_col:(end-1+dif_col)) + ...
     RRE(
2:(end-1),2:(end-1)) - RRE(2+dif_row:(end-1+dif_row), 2+dif_col:(end-1+dif_col));
     
IP(IP > Maximum) = Maximum
;                          % The Reset step

% ignore the results obtained 
in the rows or columns for which the neighbors are undefined
if (dif_col == -
1) IP(:,1) = OPE(2:(end-1),2); end
if (dif_col == +1) IP(:,end) = OPE(2:(end-1),end-1); end
if (dif_row == -1) IP(1,:) = OPE(22:(end-1)); end
if (dif_row == +1) IP(end,:) = OPE(end-12:(end-1)); end
NP = (OPE(2:(end-1),2:(end-1)) + IP)/2;              % The Averaging operation
OPE(2:(end-1), 2:(end-1)) = NP;

function Layers = ComputeLayers(nrows, ncols)
power = 
2^fix(log2(gcd(nrows, ncols)));              % start from the Greatest Common Divisor
while(power > 1 & ((rem(nrows, power) ~= 0) | (rem(ncols, power) ~= 0)))
   power = power/
2;                                  % and find the greatest common divisor
end                                                  % that is a power of 2
Layers = log2(power)
;

function Result = ImageDownResolution(A, blocksize)
[rows, cols] = size(A)
;                              % the input matrix A is viewed as
result_rows = rows/blocksize;                        % a series of square blocks
result_cols = cols/blocksize;                        % of size = blocksize
Result = zeros([result_rows result_cols]);
for crt_row = 1:result_rows                          % then each pixel is computed as
   for crt_col = 
1:result_cols                       % the average of each such block
      Result(crt_row, crt_col) = mean2(A(
1+(crt_row-1)*blocksize:crt_row*blocksize, ...
                                       
1+(crt_col-1)*blocksize:crt_col*blocksize));
   end
en
复制代码

总结

      在上述两个经典的估计算法之后经过反对数变换就可以恢复成原格式的图像。通过试验可以发现,这些算法都还是有一定缺陷的,对于照度的估计后续还有很多算法如:单尺度Retinex (SSR)、多尺度Retinex (MSR)、可变框架的Retinex等等;还有针对增强后的光晕现象先用进行照度分割,然后在再计算等等方法,本质上都大同小异。这些改进算法这里不再进行介绍,有兴趣的朋友可以去下载些论文看看。有不明白的或者本文有错误的地方也希望各位能够指出,以免误导后面的读者。