直方图均衡去雾
%自带函数
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')