滤波及对图像处理的去雾的研究

本周复习了滤波,对图像处理的去雾作了进一步研究

去雾

暗原色先验快速去雾

大气散射模型

大气散射模型描述了雾化图像的退化过程：
I(x)=J(x)t(x)+A(1-t(x));

I是观测图像的强度，

J是景物光线的强度，

A是无穷远处的大气光，

t称为透射率。

去雾的目标就是从I中复原J。方程中的第一项J(x)t(x)叫做直接衰减项，A(1−t(x))是大气光成分。

暗原色先验

暗原色先验是HEKai-ming等人发现的，是通过对大量户外无雾图像的统计观察得出的：在绝大多数图像的局部区域里，某一些像素总会有至少一个颜色通道具有很低的值。换言之，该区域光强度的最小值是很小的数。对图像J，定义：

代码

clc;
clear all;
img_name='2.jpg';
%原始图像
I=double(imread(img_name))/255;
%图像每个像元每个通道的灰度值均要除以255
%获取图像大小
[h,w,c]=size(I);%h,w代表行和列
win_size=7;%窗口大小
img_size=w*h;
figure,imshow(I);
win_dark=ones(h,w);%定义一个矩阵有h行w列
%计算分块darkchannel
for j=1+win_size:w-win_size
    for i=1+win_size:h-win_size
        m_pos_min=min(I(i,j,:));%每个(i,j)像元三通道的最小值
      for n=j-win_size:j+win_size                 %该像元的15领域
            for m=i-win_size:i+win_size
                if(win_dark(m,n)>m_pos_min)
                    win_dark(m,n)=m_pos_min;
                end
            end
        end
    end
end
figure(2),imshow(win_dark);%显示分块darkchannel
%选定精确dark value坐标
win_t=1-0.95*win_dark;%得到的win_t还是160行，243列
win_b=zeros(img_size,1);%定义一个矩阵,img_size行，1列，里面数值全为0
for ci=1:h
    for cj=1:w
        if(rem(ci-8,15)<1)%rem表示取余数
            if(rem(ci-8,15)<1)     
                win_b(ci*w+cj)=win_dark(ci*w+cj);     
            end
        end
    end
end
neb_size = 9; 
win_size = 1;  
epsilon = 0.0000001; %指定矩阵形状
indsM=reshape([1:img_size],h,w);%是不是win_dark 
%计算矩阵L
tlen = img_size*neb_size^2;
row_inds=zeros(tlen ,1);
col_inds=zeros(tlen,1);
vals=zeros(tlen,1); 
len=0;  
for j=1+win_size:w-win_size     
    for i=win_size+1:h-win_size         
        if(rem(ci-8,15)<1)             
            if(rem(cj-8,15)<1)                
                continue;            
            end
        end
        win_inds=indsM(i-win_size:i+win_size,j-win_size:j+win_size);    
        win_inds=win_inds(:);%列显示      
        winI=I(i-win_size:i+win_size,j-win_size:j+win_size,:);
              winI=reshape(winI,neb_size,c); %三个通道被拉平成为一个二维矩阵 3*9
      win_mu=mean(winI,1)';  %求每一列的均值 如果第二个参数为2 则为求每一行的均值  //矩阵变向量
      win_var=inv(winI'*winI/neb_size-win_mu*win_mu' +epsilon/neb_size*eye(c)); %求方差
      winI=winI-repmat(win_mu',neb_size,1);%求离差
      tvals=(1+winI*win_var*winI')/neb_size;% 求论文所指的矩阵L
      row_inds(1+len:neb_size^2+len)=reshape(repmat(win_inds,1,neb_size),...
                                             neb_size^2,1);
      col_inds(1+len:neb_size^2+len)=reshape(repmat(win_inds',neb_size,1),...
                                             neb_size^2,1);
      vals(1+len:neb_size^2+len)=tvals(:);
      len=len+neb_size^2;
    end
end 
 
 
vals=vals(1:len);
row_inds=row_inds(1:len);
col_inds=col_inds(1:len);
%创建稀疏矩阵
A=sparse(row_inds,col_inds,vals,img_size,img_size);
%求行的总和 sumA为列向量
sumA=sum(A,2);
%spdiags(sumA(:),0,img_size,img_size) 创建img_size大小的稀疏矩阵其元素是sumA中的列元素放在由0指定的对角线位置上。
A=spdiags(sumA(:),0,img_size,img_size)-A;

  %创建稀疏矩阵
  D=spdiags(win_b(:),0,img_size,img_size);
  lambda=1;
  x=(A+lambda*D)(lambda*win_b(:).*win_b(:));
 
   %去掉0-1范围以外的数
  alpha=max(min(reshape(x,h,w),1),0);
 
figure, imshow(alpha);
A=220/255; %大气光没有去计算
%去雾
       
for i=1:c
    for j=1:h
        for l=1:w
            dehaze(j,l,i)=(I(j,l,i)-A)/alpha(j,l)+A;
        end
    end
end
figure, imshow(dehaze);

右边圆柱的细节没有得到保留

接着用retinex实现，相比较

retinex实现颜色有些失真

Retinex实现图像去雾

步骤

读入图像→归一化→设置高斯函数参数及矩阵→高斯函数和输入图像矩阵卷积→取对数→和输入图像矩阵的对数相差→取指数→输出图像

基于SSR算法便可以实现一个基本的图像去雾程序。下面的MATLAB代码是完全按照上面的思路来实现的。只是在最后，对R(x,y)作对比度增强，以得到最终的去雾图像。此外，因为这里处理的是彩色图像，所以我们对R、G、B三个通道分别进行了处理。

实现代码

I = imread('canon.jpg');

R = I(:, :, 1);
[N1, M1] = size(R);
R0 = double(R);
Rlog = log(R0+1);
%图像I做对数变换前，需要转化为double型并做归一化。log(I+1)是为了满足真数大于0，以防计算无意义。特别提醒，一般不建议采用log(eps+I)方式
Rfft2 = fft2(R0); %对图像进行傅立叶变换,转换到频域中 

sigma = 250;
F = fspecial('gaussian', [N1,M1], sigma);
Efft = fft2(double(F));% 对高斯滤波函数进行二维傅里叶变换

DR0 = Rfft2.* Efft;% 对R分量与高斯滤波函数进行卷积运算
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);
% 对增强后的图像R、G、B分量进行融合
result = cat(3, EXPRr, EXPGg, EXPBb);
subplot(121), imshow(I);
subplot(122), imshow(result);

直方图去雾

f = imread('2.jpg');
hsi = rgb2hsi(f);
h = hsi (:, :, 1);
s = hsi (:, :, 2);
i = hsi (:, :, 3);
g = cat(3, h, s, adapthisteq(i));
g = hsi2rgb(g);
subplot(1, 2, 1), imshow(f);
subplot(1, 2, 2), imshow(g);

高斯滤波和双边滤波

matlab自带高斯滤波

clear all;clc;close all
img=imread('1.jpg');
gray=rgb2gray(img);                                     %把彩色图片转化成灰度图片
figure(1),imshow(gray),title('彩色原始图像转灰色图像）'); %显示原始图像
gray=imnoise(gray,'gaussian',0,0.002);                  %加入均值为0，方差为0.001的高斯噪声
figure(2),imshow(gray),title('加入高斯噪声之后的图像');   %显示加入高斯噪声之后的图像

% 用matlab系统函数进行高斯滤波
sigma=0.5;                                                %滤波器的标准值，单位为像素
hsize=[8 8];                                              %模板尺寸
gsseq=fspecial('gaussian',hsize,sigma);                   %生成高斯序列
Y1=filter2(gsseq,gray)/255;                               %用生成的高斯序列进行滤波
figure(3),imshow(Y1),title('用Matlab自带函数进行高斯滤波'); %显示滤波后的图像

傅里叶变换实现高斯滤波

clear all;clc;close all
img=imread('1.jpg');
gray=rgb2gray(img);                                     %把彩色图片转化成灰度图片
figure,imshow(gray),title('彩色原始图像转灰色图像）'); %显示原始图像
gray=imnoise(gray,'gaussian',0,0.002);                  %加入均值为0，方差为0.001的高斯噪声
figure,imshow(gray),title('加入高斯噪声之后的图像');   %显示加入高斯噪声之后的图像
gray=double(gray);
gray=fft2(gray);%二维傅里叶变换
gray=fftshift(gray);%频率居中
[m,n]=size(gray);
d0=90%标准差
m1=fix(m/2);
n1=fix(n/2);
for i=1:m
    for j=1:n
        d=sqrt((i-m1)^2+(j-n1)^2);%计算像素点到图像中心的距离
        h(i,j)=exp(-d^2/2/d0^2);%高斯滤波器
    end
end
g=gray.*h;%将图像进行高斯滤波，频域上表现为为两个函数相乘
g=ifftshift(g);                   %频域圆周移位
g=ifft2(g);                       %二维傅里叶反变换
g=mat2gray(real(g));              %归一化
figure,imshow(g),title('用重新编写的程序进行高斯滤波');%显示滤波后的图像

双边滤波

function B = bfltGray(A,w,sigma_d,sigma_r)
% 计算距离因子权重
[X,Y] = meshgrid(-w:w,-w:w);
%创建核距离矩阵
%e.g.
%[x,y]=meshgrid(-1:1,-1:1)
% 
% x =
% 
%     -1     0     1
%     -1     0     1
%     -1     0     1
% 
% 
% y =
% 
%     -1    -1    -1
%      0     0     0
%      1     1     1

%计算定义域核
G = exp(-(X.^2+Y.^2)/(2*sigma_d^2));

%创建进度条
h = waitbar(0,'Applying bilateral filter...');
set(h,'Name','Bilateral Filter Progress');

% 应用双边滤波
%计算值域核H 并与定义域核G 乘积得到双边权重函数F
dim = size(A);
B = zeros(dim);
for i = 1:dim(1)
   for j = 1:dim(2)
         %边界限制
         iMin = max(i-w,1);
         iMax = min(i+w,dim(1));
         jMin = max(j-w,1);
         jMax = min(j+w,dim(2));
         
         %定义当前核所作用的区域为(iMin:iMax,jMin:jMax)
         I = A(iMin:iMax,jMin:jMax);%提取该区域的源图像值赋给I
      
         % 计算亮度因子权重
         H = exp(-(I-A(i,j)).^2/(2*sigma_r^2));
      
         % 计算双边滤波结果
         F = H.*G((iMin:iMax)-i+w+1,(jMin:jMax)-j+w+1);
         B(i,j) = sum(F(:).*I(:))/sum(F(:));
               
   end
   waitbar(i/dim(1));
end
% 结束进度条
close(h);

%A为给定图像，归一化到[0,1]的矩阵
%w为双边滤波器（核）的边长/2
%定义域方差σd记为SIGMA(1),值域方差σr记为SIGMA(2)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function B = bfilter2(A,w,sigma)
% 检验给定图像是否存在并且有效
if ~exist('A','var') || isempty(A)
   error('Input image A is undefined or invalid.');
end
if ~isfloat(A) || ~sum([1,3] == size(A,3)) || ...
      min(A(:)) < 0 || max(A(:)) > 1
   error(['Input image A must be a double precision ',...
          'matrix of size NxMx1 or NxMx3 on the closed ',...
          'interval [0,1].']);      
end

% 检验双边滤波器的半宽是否符合要求
if ~exist('w','var') || isempty(w) || ...
      numel(w) ~= 1 || w < 1
   w = 5;
end
w = ceil(w);

% 检验sigma参数是否符合要求
if ~exist('sigma','var') || isempty(sigma) || ...
      numel(sigma) ~= 2 || sigma(1) <= 0 || sigma(2) <= 0
   sigma = [3 0.1];
end

%选择彩色模式或灰度模式
if size(A,3) == 1
   B = bfltGray(A,w,sigma(1),sigma(2));
else
   B = bfltColor(A,w,sigma(1),sigma(2));
end

%调用示例
I=imread('noise.tif'); %读入图片
I=double(I)/255;  %转为double型并归一化
w     = 5;        % 双边滤波器半宽，w越大平滑作用越强
sigma = [3 0.1];  % 空间距离方差σd记为SIGMA(1),像素亮度方差σr记为SIGMA(2),即空间邻近度因子和亮度相似度因子的衰减程度
I1=bfilter2(I,w,sigma);                     %双边滤波器滤波
figure(1),imshow(I),title('原始图像');       %作出原始图像
figure(2),imshow(I1),title('双边滤波后的图像')%作出双边滤波后的图像