图像增强去雾之直方图均衡化/同态滤波/Retinex算法

图像增强去雾之直方图均衡化/同态滤波/Retinex算法

  最近撸了一发图像去雾的算法，主要举四个例子，分别用了全局直方图均衡化，局部直方图均衡化，同态滤波，Retinex增强算法。感兴趣的可以一起讨论一下啊

  1.全局直方图均衡化

  直方图均衡化处理的“中心思想”是把原始图像的灰度直方图从比较集中的某个灰度区间变成在全部灰度范围内的均匀分布。直方图均衡化就是对图像进行非线性拉伸，重新分配图像像素值，使一定灰度范围内的像素数量大致相同。直方图均衡化就是把给定图像的直方图分布改变成“均匀”分布直方图分布。举一个适合做直方图均衡化的去雾例子，原图的雾霾比较均匀，可以看到原灰度直方图的灰度区间比较集中，均衡化以后就几乎变为全部灰度范围的均匀分布了。处理后灰度范围变大，对比度变大，清晰度变大，所以能有效增强图像。

  全局直方图均衡化主要是以下几步：

  1）统计直方图每个灰度出现的次数；

  2）累积归一化直方图；

  3）计算新的像素值。


 

Xs=imread('3.png');
subplot(2,1,1);
imshow(Xs);
title('原图像', 'FontWeight', 'Bold');
R=Xs(:,:,1);
G=Xs(:,:,2);
B=Xs(:,:,3);
M=histeq(R);
N=histeq(G);
L=histeq(B);
In=cat(3,M,N,L);
subplot(2,1,2);
imshow(In);
title('全局直方图处理后的图像', 'FontWeight', 'Bold');

  2.局部直方图均衡化

  举例中的图不像第一张图雾霾分布那么平均，而是集中在某一区域，所以比较适合局部直方图均衡化，简单的上步骤上图吧：

  1）确定模板大小 n*n；

  2）将图像进行扩展，因为对边界处的处理会使得图像无法与模板达到一一对应；

  3）从图像的第一个像素开始，与模版点乘，点乘后的局部区域进行直方图均值化，并将局部的中心元素的作为图像的当前值。


 
 

Img=imread('2.png');
imshow(Img);
g1=GetLocalHisteq(Img(:,:,1));
g2=GetLocalHisteq(Img(:,:,2));
g3=GetLocalHisteq(Img(:,:,3));
ImgS=cat(3,g1,g2,g3);
subplot(1,2,1);imshow(Img);subplot(1,2,2);imshow(ImgS);

function g = GetLocalHisteq(I)
x=mat2gray(I);
f=im2double(x);
w=4;
k=0.06;
M=mean2(f);
z=colfilt(f,[w w],'sliding',@std);
m=colfilt(f,[w w],'sliding',@mean);
A=k*M./z;
g=A.*(f-m)+m;
g=im2uint8(mat2gray(g));

  3.Retinex增强

  Retinex理论的基础理论是物体的颜色是由物体对长波（红色）、中波（绿色）、短波（蓝色）光线的反射能力来决定的，而不是由反射光强度的绝对值来决定的，物体的色彩不受光照非均匀性的影响，具有一致性，即retinex是以色感一致性（颜色恒常性）为基础的。不同于传统的线性、非线性的只能增强图像某一类特征的方法，Retinex可以在动态范围压缩、边缘增强和颜色恒常三个方面打到平衡，因此可以对各种不同类型的图像进行自适应的增强。

  Retinex算法默认一幅给定的图像S(x,y)可以分解为两个不同的图像：反射图像R(x,y)和亮度图像（也有人称之为入射图像）L(x,y)的乘积。而我们把入射图像(低频部分)去掉以后，留下的就是我们希望得到的增强边缘以后的结果了。先来看一下算法的处理流程图：


 
  具体步骤：

  1）原图S(x,y)=反射部分R(x,y)*照射部分L(x,y)

  2）去除照射部分（低频），只留下反射部分（高频/边缘）

      


  这里要重点提一下公式L(x,y)=F(x,y)*S(x,y)

  这不是乘法而是一步卷积操作，F(x,y)是中心环绕函数，其实可以把他理解为低通函数。这一步卷积的目的其实就是通过计算图片中像素点与周围区域的加权平均，来估计出图像中的照度部分(也就是低频部分)，从而把低频部分减掉，只留下高频的边缘部分。

  最后，给出Retinex算法的处理结果与直方图显示（原图实在太不清晰，这已经是我能做出最好的结果了。如有更好的算法可以积极讨论啊）


 

I = imread('4.png');
R = I(:, :, 1);
[N1, M1] = size(R);
R0 = double(R);
Rlog = log(R0+1);
Rfft2 = fft2(R0);

sigma = 250;
F = fspecial('gaussian', [N1,M1], sigma);
Efft = fft2(double(F));

DR0 = Rfft2.* Efft;
DR = ifft2(DR0);

DRlog = log(DR +1);
Rr = Rlog - DRlog;
EXPRr = exp(Rr);
MIN = min(min(EXPRr));
MAX = max(max(EXPRr));
EXPRr = (EXPRr - MIN)/(MAX - MIN);
EXPRr = adapthisteq(EXPRr);

G = I(:, :, 2);

G0 = double(G);
Glog = log(G0+1);
Gfft2 = fft2(G0);

DG0 = Gfft2.* Efft;
DG = ifft2(DG0);

DGlog = log(DG +1);
Gg = Glog - DGlog;
EXPGg = exp(Gg);
MIN = min(min(EXPGg));
MAX = max(max(EXPGg));
EXPGg = (EXPGg - MIN)/(MAX - MIN);
EXPGg = adapthisteq(EXPGg);

B = I(:, :, 3);

B0 = double(B);
Blog = log(B0+1);
Bfft2 = fft2(B0);

DB0 = Bfft2.* Efft;
DB = ifft2(DB0);

DBlog = log(DB+1);
Bb = Blog - DBlog;
EXPBb = exp(Bb);
MIN = min(min(EXPBb));
MAX = max(max(EXPBb));
EXPBb = (EXPBb - MIN)/(MAX - MIN);
EXPBb = adapthisteq(EXPBb);

result = cat(3, EXPRr, EXPGg, EXPBb);
subplot(121), imshow(I);
subplot(122), imshow(result);

  4.同态滤波

  同态滤波和Retinex增强算法其实很相近，都是把一张图片看作照射光和反射光组成的，照射部分图像缓慢变化，属于低频；反射部分主要就是急剧变化的交界处/边缘，也就是高频部分。要想消除噪声，提高对比度，实际上就是衰减低频，增高高频，也就达到了图像增强的目的。

  同态滤波算法流程图如下：

  原图S(x,y)->取对数log->傅里叶变化DFT->频域滤波->反傅里叶IDFT->求指数Exp->输出T(x,y)

  1）原图f(x,y)=照射部分i(x,y)*反射部分r(x,y)

  2）取对数 ln f(x,y)=ln i(x,y)+ln r(x,y)

  3）傅里叶变化: F(x,y)=I(x,y)+R(x,y)

  4）频域滤波选定函数H(x,y)：H(x,y)F(x,y)=H(x,y)I(x,y)+H(x,y)R(x,y)

  图像对数傅里叶变换中的低频部分主要对应照度分量，而高频部分主要对应反射分量。我们可以设计1个对傅里叶变换的高频分量和低频分量影响不同的滤波函数H(x,y)。


 
  5）反傅里叶（反变换到空域）: 


  6）取指数：输出图


  最后给出同态滤波处理过的结果图，个人感觉不如Retinex增强算法的效果好。

[name,path] = uigetfile('4.png');
file = strcat(path,name);
[X,map]=imread(file);
X=double(X);
I=rgb2hsv(X);
H=I(:,:,1);
S=I(:,:,2);
V=I(:,:,3);
%if size(X,3)==3
%   X= rgb2gray(X);
%end
% 装载图片
% 显示这个图片
figure,imshow('4.png');
title('原始图像');
% 构造一个高斯滤波器
f_high = 1.0;
f_low = 0.8;
% 得到一个高斯低通滤波器
gauss_low_filter = fspecial('gaussian', [7 7], 1.414);
matsize = size(gauss_low_filter);
% 由于同态滤波是要滤出高频部分,
% 所以得把这个低通滤波器转换成一个高通滤波器.
% f_high 和 f_low 是控制这个高通滤波器形态的参数.
gauss_high_filter = zeros(matsize);
gauss_high_filter(ceil(matsize(1,1)/2) , ceil(matsize(1,2)/2)) = 1.0;
gauss_high_filter = f_high*gauss_high_filter - (f_high-f_low)*gauss_low_filter;
% 利用对数变换将入射光和反射光部分分开
log_img = log(double(V)+eps);

% 将高斯高通滤波器与对数转换后的图象卷积
high_log_part = imfilter(log_img, gauss_high_filter, 'symmetric', 'conv');
% 由于被处理的图象是经过对数变换的,再用幂变换将图象恢复过来
high_part = exp(high_log_part);
minv = min(min(high_part));
maxv = max(max(high_part));
% 得到的结果图象
temp=(high_part-minv)/(maxv-minv);

S=adapthisteq(S)*2.1;
HSV3(:,:,1)=H;       %保留H不变，开始合成
HSV3(:,:,2)=S;
HSV3(:,:,3)=temp;
rgb2=hsv2rgb(HSV3);    %转换回RGB空间
figure;
imshow(rgb2);
title('同态滤波图像');