本周复习了滤波,对图像处理的去雾作了进一步研究
去雾
暗原色先验快速去雾
大气散射模型
大气散射模型描述了雾化图像的退化过程：
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('双边滤波后的图像')%作出双边滤波后的图像