经典去雾算法

直方图均衡去雾

%自带函数
clear; clc;close all
tic
I=imread('test6.jpg'); %return RGB

I1(:,:,1)=histeq(I(:,:,1));%灰度级数目64
I1(:,:,2)=histeq(I(:,:,2));
I1(:,:,3)=histeq(I(:,:,3));
%M=rgb2gray(I);%将RGB图转为灰度图
%N=rgb2gray(I1);
t=toc
figure;
%subplot(321);
imshow(I),title('原图');
%subplot(322);
figure;
imshow(I1),title('全局均衡化后的图像');
%subplot(323);imhist(I),title('原图直方图');%
%subplot(324);imhist(I1),title('全局均衡化后的图像的直方图');%直方图描述的是图像中具有各灰度级的出现的概率(像素的个数),其横坐标为灰度级,纵坐标为图像中具有该灰度级的像素个数。由于灰度级的大小为0-255,故横坐标的数值范围为0-255.
%subplot(325);imshow(M);title('原图的灰度图');
%subplot(326);imshow(N);title('全局均衡化后的图像的灰度图');
%imwrite(I1,'test6whole.jpg')

 

%手动代码
clear all;clc;
tic
image=imread('test7.jpg');                %读入BMP彩色图像文件
%imwrite(rgb2gray(PS),'PicSampleGray.bmp'); %将彩色图片灰度化并保存
R=image(:,:,1);                          %灰度化后的数据存入数组
[m,n]=size(R);  %测量图像尺寸参数
GPR=zeros(1,256);                            %预创建存放灰度出现概率的向量
for k=0:255
     GPR(k+1)=length(find(R==k))/(m*n);      %计算每级灰度出现的概率,将其存入GP中相应位置
end
%直方图均衡化
S1R=zeros(1,256);
for i=1:256
     for j=1:i
          S1R(i)=GPR(j)+S1R(i);                 %计算SR
     end
end
S2R=round((S1R*256)+0.5);                          %将Sk归到相近级的灰度

PAR=R;
for i=0:255
     PAR(find(R==i))=S2R(i+1);                %将各个像素归一化后的灰度值赋给这个像素
end

G=image(:,:,2);                          %灰度化后的数据存入数组
[m,n]=size(R);                             %测量图像尺寸参数
GPG=zeros(1,256);                            %预创建存放灰度出现概率的向量
for k=0:255
     GPG(k+1)=length(find(G==k))/(m*n);      %计算每级灰度出现的概率,将其存入GP中相应位置
end
%直方图均衡化
S1G=zeros(1,256);
for i=1:256
     for j=1:i
          S1G(i)=GPG(j)+S1G(i);                 %计算Sk
     end
end
S2G=round((S1G*256)+0.5);                          %将Sk归到相近级的灰度

PAG=G;
for i=0:255
     PAG(find(G==i))=S2G(i+1);                %将各个像素归一化后的灰度值赋给这个像素
end
%figure,imshow(PAG)                           %显示均衡化后的图像
%title('均衡化后图像')

B=image(:,:,3);                          %灰度化后的数据存入数组

%绘制直方图
[m,n]=size(B);                             %测量图像尺寸参数
GPB=zeros(1,256);                            %预创建存放灰度出现概率的向量
for k=0:255
     GPB(k+1)=length(find(B==k))/(m*n);      %计算每级灰度出现的概率,将其存入GP中相应位置
end

S1B=zeros(1,256);
for i=1:256
     for j=1:i
          S1B(i)=GPB(j)+S1B(i);                 %计算Sk
     end
end
S2B=round((S1B*256)+0.5);                          %将Sk归到相近级的灰度

PAB=B;
for i=0:255
     PAB(find(B==i))=S2B(i+1);                %将各个像素归一化后的灰度值赋给这个像素
end
out=cat(3,PAR,PAG,PAB);
t=toc
figure;
imshowpair(image, out, 'montage');
%imwrite(out,'test6whole2.jpg');


暗通道去雾

%导向滤波
clear;clc;close all
w0=0.95;
t0=0.001;
path0= 'D:\learnapp\matlab\bin\code\defog\persontask\fog1\';
I=imread([path0,'1258_0.85_0.2.jpg']);



[h,w,~]=size(I);
kenlRatio = 0.01;
minAtomsLight = 240;
min_I= zeros(h,w);

for y=1:h

    for x=1:w
        min_I(y,x) = min(I(y,x,:));
         end

end
%figure,imhist(uint8(min_I)), title('Min(R,G,B)');%每个像素RGB分量中的最小值,形成灰度图

krnlsz = floor(max([3, w*kenlRatio, h*kenlRatio]));%最小值滤波半径
dark_I = ordfilt2(min_I,1,ones(krnlsz,krnlsz),'symmetric'); 

dark_channel=double(dark_I);
%%

A= min([minAtomsLight, max(max(dark_I))]); 

%%

t=1-w0*(dark_channel/A);  %t 

t1=max(t,t0); %清除透射率图中可能出现的零值
%figure;imshow(uint8(t*255));title('(c)透射率图');



J = zeros(h,w,3);
I1 = double(I);
J(:,:,1) = (I1(:,:,1) - (1-t1)*A)./t1;
J(:,:,2) = (I1(:,:,2) - (1-t1)*A)./t1;
J(:,:,3) = (I1(:,:,3) - (1-t1)*A)./t1;
out=uint8(J);
%J(:,:,1) = A + (I1(:,:,1)-A)./t1; J(:,:,2) = A + (I1(:,:,2)-A)./t1; J(:,:,3) = A + (I1(:,:,3)-A)./t1;
%figure;imshow(out); title('(d)去雾图像');


r = krnlsz*4;%当r比较小的时候,在透射率图中基本看不到什么细节信息,因此恢复处的图像边缘处不明显。在前面进行最小值后暗通道的图像成一块一块的,为了使透射率图更加精细,建议这个r的取值不小于进行最小值滤波的半径的4倍
eps = 10^-6;%为了防止计算中除以0的错误以及为了使得某些计算结果不至于过大,一般建议取值0.001或者更小。

% filtered = guidedfilter_color(double(img)/255, t_d, r, eps);%如果使用的彩色RGB图做导向图,计算时间上会增加不少,所的到的透射率图的边缘会比灰度图所处理的保留了更多的细节,效果上略微比灰度图好
filtered = guidedfilter(double(rgb2gray(I))/255, t1, r, eps);%导向图,透射率图 
t_d = filtered;
%figure;imshow(t_d),title('(e)导向滤波的透射率图');

J(:,:,1) = (I1(:,:,1) - (1-t_d)*A)./t_d;
J(:,:,2) = (I1(:,:,2) - (1-t_d)*A)./t_d;
J(:,:,3) = (I1(:,:,3) - (1-t_d)*A)./t_d;
out=uint8(J);
figure;
imshowpair(I,out,'montage');










function q = guidedfilter(I, p, r, eps)
%   GUIDEDFILTER   O(1) time implementation of guided filter.
%
%   - guidance image: I (should be a gray-scale/single channel image)
%   - filtering input image: p (should be a gray-scale/single channel image)
%   - local window radius: r
%   - regularization parameter: eps

[hei, wid] = size(I);
N = boxfilter(ones(hei, wid), r); % the size of each local patch; N=(2r+1)^2 except for boundary pixels.


mean_I = boxfilter(I, r) ./ N;
mean_p = boxfilter(p, r) ./ N;
mean_Ip = boxfilter(I.*p, r) ./ N;
cov_Ip = mean_Ip - mean_I .* mean_p; % this is the covariance of (I, p) in each local patch.

mean_II = boxfilter(I.*I, r) ./ N;
var_I = mean_II - mean_I .* mean_I;%方差variance,平方的期望-期望的平方

a = cov_Ip ./ (var_I + eps); % Eqn. (5) in the paper;
b = mean_p - a .* mean_I; % Eqn. (6) in the paper;

mean_a = boxfilter(a, r) ./ N;
mean_b = boxfilter(b, r) ./ N;

q = mean_a .* I + mean_b; % Eqn. (8) in the paper;
end









function imDst = boxfilter(imSrc, r)

%   BOXFILTER   O(1) time box filtering using cumulative sum
%
%   - Definition imDst(x, y)=sum(sum(imSrc(x-r:x+r,y-r:y+r)));
%   - Running time independent of r; 
%   - Equivalent to the function: colfilt(imSrc, [2*r+1, 2*r+1], 'sliding', @sum);
%   - But much faster.

[hei, wid] = size(imSrc);
imDst = zeros(size(imSrc));

%cumulative sum over Y axis
imCum = cumsum(imSrc, 1);%在列方向求积累和
%difference over Y axis
imDst(1:r+1, :) = imCum(1+r:2*r+1, :);
imDst(r+2:hei-r, :) = imCum(2*r+2:hei, :) - imCum(1:hei-2*r-1, :);
imDst(hei-r+1:hei, :) = repmat(imCum(hei, :), [r, 1]) - imCum(hei-2*r:hei-r-1, :);

%cumulative sum over X axis
imCum = cumsum(imDst, 2);
%difference over Y axis
imDst(:, 1:r+1) = imCum(:, 1+r:2*r+1);
imDst(:, r+2:wid-r) = imCum(:, 2*r+2:wid) - imCum(:, 1:wid-2*r-1);
imDst(:, wid-r+1:wid) = repmat(imCum(:, wid), [1, r]) - imCum(:, wid-2*r:wid-r-1);
end

Retinex去雾

%MSR
clear all;
tic
I = imread('11.jpg');
R = I(:, :, 1);              %R通道的二维图像
G = I(:, :, 2);              %G通道的二维图像
B = I(:, :, 3);              %B通道的二维图像
R0 = im2double(R);        %将图像的数据类型转换为double类型
G0 = im2double(G);
B0 = im2double(B);
[m,n] = size(R);
Rlog = log(R0+1);
Rfft2 = fft2(R0);

sigma1 = 30;          %尺度参数之一
F1 = fspecial('gaussian', [m,n], sigma1);       %建立高斯低通滤波算子
Efft1 = fft2(double(F1));
DR0 = Rfft2.* Efft1; %滤波算子和原始R通道图像进行卷积,得到空间平滑图像
DR = ifft2(DR0);
DRlog = log(DR +1);
Rr1 = Rlog - DRlog;  %R通道图像减去空间平滑图像,得到R通道的物体反射特性
 

sigma2 = 90;        %尺度参数之一,同sigma1
F2 = fspecial('gaussian', [m,n], sigma2);
Efft2 = fft2(double(F2));
DR0 = Rfft2.* Efft2;
DR = ifft2(DR0);
DRlog = log(DR +1);
Rr2 = Rlog - DRlog;

sigma3 = 256;
F3 = fspecial('gaussian', [m,n], sigma3);
Efft3 = fft2(double(F3));
DR0 = Rfft2.* Efft3;
DR = ifft2(DR0);
DRlog = log(DR +1);
Rr3 = Rlog - DRlog;

Rr = (Rr1 + Rr2 +Rr3)/3;  %为三个尺度取权重因子
       
EXPRr = exp(Rr);
MIN = min(min(EXPRr));
MAX = max(max(EXPRr));
EXPRr = (EXPRr - MIN)/(MAX - MIN)*255;   %将输出量化
EXPRr=uint8(EXPRr);
EXPRr = adapthisteq(EXPRr);     %进行直方图均衡


%G通道与R通道实现思路相同
Glog = log(G0+1);
Gfft2 = fft2(G0);

DG0 = Gfft2.* Efft1;
DG = ifft2(DG0);
DGlog = log(DG +1);
Gg1 = Glog - DGlog;

DG0 = Gfft2.* Efft2;
DG = ifft2(DG0);
DGlog = log(DG +1);
Gg2 = Glog - DGlog;

DG0 = Gfft2.* Efft3;
DG = ifft2(DG0);
DGlog = log(DG +1);
Gg3 = Glog - DGlog;

Gg = (Gg1 + Gg2 +Gg3)/3;
EXPGg = exp(Gg);
MIN = min(min(EXPGg));
MAX = max(max(EXPGg));
EXPGg = (EXPGg - MIN)/(MAX - MIN)*255;
EXPGg=uint8(EXPGg);
EXPGg = adapthisteq(EXPGg);

%B通道与R通道实现思路相同
Blog = log(B0+1);
Bfft2 = fft2(B0);

DB0 = Bfft2.* Efft1;
DB = ifft2(DB0);
DBlog = log(DB +1);
Bb1 = Blog - DBlog;

DB0 = Bfft2.* Efft2;
DB = ifft2(DB0);
DBlog = log(DB +1);
Bb2 = Blog - DBlog;

DB0 = Bfft2.* Efft3;
DB = ifft2(DB0);
DBlog = log(DB +1);
Bb3 = Blog - DBlog;
Bb = (Bb1 + Bb2 +Bb3)/3;

EXPBb = exp(Bb);
MIN = min(min(EXPBb));
MAX = max(max(EXPBb));
EXPBb = (EXPBb - MIN)/(MAX - MIN)*255;
EXPBb=uint8(EXPBb);
EXPBb = adapthisteq(EXPBb);

out = cat(3, EXPGg, EXPGg, EXPBb);  %将R、G、B三个图像的矩阵连结在一起
t=toc
%figure(1);imshow(I);title('原图');
%figure(2);imshow(out);title('msr');
%imwrite(result,'test3msr.jpg')
imshowpair(I, out, 'montage');
%MSRCR
clear;close all;

I = imread('11.jpg');         %读取图像
R = I(:, :, 1);              %R通道的二维图像
G = I(:, :, 2);              %G通道的二维图像
B = I(:, :, 3);              %B通道的二维图像
R0 = im2double(R);        %将图像的数据类型转换为double类型
G0 = im2double(G);
B0 = im2double(B);
[m,n] = size(R);

a = 50;                %MSRCR算法中的调整参数
II = imadd(R0, G0);             
II = imadd(II, B0);       %将R、G、B三个通道的图叠加在一起
Ir = immultiply(R0, a);    %R通道的图像乘以调整参数
Ig = immultiply(G0, a);    %G通道的图像乘以调整参数
Ib = immultiply(B0, a);    %B通道的图像乘以调整参数

Rlog = log(R0+1);
Rfft2 = fft2(R0);

sigma1 = 64;          %尺度参数之一
F1 = fspecial('gaussian', [m,n], sigma1);       %建立高斯低通滤波算子
Efft1 = fft2(double(F1));
DR0 = Rfft2.* Efft1; %滤波算子和原始R通道图像进行卷积,得到空间平滑图像
DR = ifft2(DR0);
DRlog = log(DR +1);
Rr1 = Rlog - DRlog;  %R通道图像减去空间平滑图像,得到R通道的物体反射特性 
 

sigma2 = 128;        %尺度参数之一,同sigma1
F2 = fspecial('gaussian', [m,n], sigma2);
Efft2 = fft2(double(F2));
DR0 = Rfft2.* Efft2;
DR = ifft2(DR0);
DRlog = log(DR +1);
Rr2 = Rlog - DRlog;

sigma3 = 256;
F3 = fspecial('gaussian', [m,n], sigma3);
Efft3 = fft2(double(F3));
DR0 = Rfft2.* Efft3;
DR = ifft2(DR0);
DRlog = log(DR +1);
Rr3 = Rlog - DRlog;

Rr = (Rr1 + Rr2 +Rr3)/3;  %为三个尺度取权重因子
       
alpha = imdivide(Ir, II);   % alpha为彩色恢复因子
alpha = log(alpha+1);    

Rr = immultiply(alpha, Rr);  %将彩色恢复因子与得到的R通道物体反射特性相乘
EXPRr = exp(Rr);
MIN = min(min(EXPRr));
MAX = max(max(EXPRr));
EXPRr = (EXPRr - MIN)/(MAX - MIN)*255;   %将输出量化
EXPRr=uint8(EXPRr);
EXPRr = adapthisteq(EXPRr);     %进行直方图均衡

%G通道与R通道实现思路相同
Glog = log(G0+1);
Gfft2 = fft2(G0);

DG0 = Gfft2.* Efft1;
DG = ifft2(DG0);
DGlog = log(DG +1);
Gg1 = Glog - DGlog;

DG0 = Gfft2.* Efft2;
DG = ifft2(DG0);
DGlog = log(DG +1);
Gg2 = Glog - DGlog;

DG0 = Gfft2.* Efft3;
DG = ifft2(DG0);
DGlog = log(DG +1);
Gg3 = Glog - DGlog;

Gg = (Gg1 + Gg2 +Gg3)/3;

alpha = imdivide(Ig, II);
alpha = log(alpha+1);

Gg = immultiply(alpha, Gg);
EXPGg = exp(Gg);
MIN = min(min(EXPGg));
MAX = max(max(EXPGg));
EXPGg = (EXPGg - MIN)/(MAX - MIN)*255;
EXPGg=uint8(EXPGg);
EXPGg = adapthisteq(EXPGg);

%B通道与R通道实现思路相同
Blog = log(B0+1);
Bfft2 = fft2(B0);

DB0 = Bfft2.* Efft1;
DB = ifft2(DB0);
DBlog = log(DB +1);
Bb1 = Blog - DBlog;

DB0 = Bfft2.* Efft2;
DB = ifft2(DB0);
DBlog = log(DB +1);
Bb2 = Blog - DBlog;

DB0 = Bfft2.* Efft3;
DB = ifft2(DB0);
DBlog = log(DB +1);
Bb3 = Blog - DBlog;
Bb = (Bb1 + Bb2 +Bb3)/3;

alpha = imdivide(Ib, II);
alpha = log(alpha+1);

Bb = immultiply(alpha, Bb);
EXPBb = exp(Bb);
MIN = min(min(EXPBb));
MAX = max(max(EXPBb));
EXPBb = (EXPBb - MIN)/(MAX - MIN)*255;
EXPBb=uint8(EXPBb);
EXPBb = adapthisteq(EXPBb);

out = cat(3, EXPRr, EXPGg, EXPBb);  %将R、G、B三个图像的矩阵连结在一起

%figure;
%subplot(2,2,1);imshow(I);title('原图')   %输出原图像
%subplot(222);
imshowpair(I,out,'montage');title('msrcr')
%imwrite(out,'test5msrcr.jpg')

  • 4
    点赞
  • 22
    收藏
    觉得还不错? 一键收藏
  • 4
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 4
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值