www.pudn.com > Dark-Channel-Prior-Codes.rar > hescodes_1.m, change:2014-01-12,size:6495b


%何凯明的 "Single Image Haze Removal Using Dark Channel Prior" 论文的代码实现 
 
 
clc;clear;warning off;close all; %#ok<WNOFF> 
win_size = 15;%局部窗口大小为15*15 
Omege=0.95;%The nice property of this modi?cation is that we adaptively keep more haze for the distant objects.\ 
Beta=0.1;%I want to use this and caculate the depth map 
 
img_name='02.png'; 
IM=double(imread(img_name))/255;%影像规范化到0-1 
%figure(1), imshow(IM);title('原图'); 
[M,N,P]=size(IM); 
img_size=N*M;%图像大小 
dehaze=zeros(M,N,P); 
 
%得到黑色通道影像 
%------------------------------------------------------------------------------------------------------- 
if ~exist(sprintf('(%s)_TxyEstimated.mat',img_name),'file') 
    tempFilter=zeros(M,N,P); 
    for i=1:P 
        tempFilter(:,:,i)=ordfilt2(IM(:,:,i),1,true(win_size),'symmetric');%symmetric 边缘处取对称的有效数据 
    end 
    DarkChannel=min(tempFilter,[],3); 
    TxyEstimated=1-Omege*DarkChannel;%Omege 对于较较远目标,调整保留一定量的雾霾 
    save(sprintf('(%s)_TxyEstimated.mat',img_name),'TxyEstimated','DarkChannel') 
else 
    load(sprintf('(%s)_TxyEstimated.mat',img_name)); 
end 
 
figure(2), imshow(TxyEstimated);title('Estimate the transmission(t)'); 
 
%利用黑色通道来提高大气光的估计精度 
DarkChannel_vector=reshape(DarkChannel,M*N,1); 
[~,Index]=sort(DarkChannel_vector,1,'descend'); 
Location=Index(1:round(M*N*0.001));%获得黑色通道前千分之一的最亮像素的图像位置 
temp=reshape(IM,M*N,P); 
temp=temp(Location,:);Length=size(temp,1); 
[AtomsphericLight,Index]=max(temp);%天空亮度 
fprintf(1,sprintf('AtomsphericLight is %f\n',AtomsphericLight)); 
temp=rem(Index,Length);%取余.得到所在行数 
LocalUse=Location(temp); 
 
figure(3),imshow(IM,[]);hold on;title('AtomsphericLight','FontSize',10); 
plot(ceil(Location/M),rem(Location,M),'r.','MarkerSize',3); 
plot(ceil(LocalUse/M),rem(LocalUse,M),'rs','MarkerSize',30,'MarkerEdgeColor',[1,0,0],'LineWidth',3); 
plot(ceil(LocalUse/M),rem(LocalUse,M),'r+','MarkerSize',10,'MarkerEdgeColor',[0,1,0]); 
text(0,M+20,sprintf('AtomsphericLight:%f,%f,%f',AtomsphericLight(1),AtomsphericLight(2),AtomsphericLight(3)),'FontSize',15,'Color',[.9,.2,.4]); 
hold off; 
pause(1); 
 
 
 
%选定精确dark_value坐标 
%--------------------------------------------------------------------------------------------- 
if ~exist(sprintf('(%s)_win_b.mat',img_name),'file') 
    win_b = zeros(img_size,1); 
    for ci=1:M 
        for cj=1:N 
             
            if(rem(ci-8,15)<1) 
                if(rem(cj-8,15)<1)%rem取余 
                    %这么做的目的是-在全矩阵中得到稀疏点--均匀分布在矩阵中 
                    win_b(ci*N+cj)=TxyEstimated(ci*N+cj); 
                end 
            end 
             
        end 
    end 
    save(sprintf('(%s)_win_b.mat',img_name),'win_b'); 
else 
    load(sprintf('(%s)_win_b.mat',img_name)); 
end 
 
 
%下面计算Soft Matting 
if ~exist(sprintf('(%s)_uint.mat',img_name),'file') 
    neb_size = 9;%he number of pixels in the window |Omega_k| 
    win_size = 1; 
    epsilon = 0.0000001;% a regularizing parameter 
    %指定矩阵形状 
    indsM=reshape([1:img_size],M,N); %#ok<NBRAK>    %MRF 标号场 
    %计算矩阵L 
    tlen = img_size*neb_size^2;%M*N*9^2 
    row_inds=zeros(tlen ,1); 
    col_inds=zeros(tlen,1); 
    vals=zeros(tlen,1); 
    len=0; 
    temp=(M-2)*(N-2);temp1=0; 
     
    for j=1+win_size:N-win_size 
        for i=win_size+1:M-win_size 
             
            temp1=temp1+1; 
            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);%indsM矩阵的一个3*3小块 
            win_inds=win_inds(:);%转换成为一列,9*1 
            winI=IM(i-win_size:i+win_size,j-win_size:j+win_size,:);%对应彩色影像的3*3*3小立方体影像 
            winI=reshape(winI,neb_size,P);%排成9*3 
            win_mu=mean(winI,1)';%按照winI的列求均值,然后转置,3*1 
            % win_var=inv(winI'*winI/neb_size-win_mu*win_mu' +epsilon/neb_size*eye(P));%△ 
            win_var=winI'*winI/neb_size-win_mu*win_mu' +epsilon/neb_size*eye(P);%△ 
             
            winI=winI-repmat(win_mu',neb_size,1);   % I_i-Mu_k 
            tvals=(1+winI/win_var*winI')/neb_size; 
             
            row_inds(1+len:neb_size^2+len)=reshape(repmat(win_inds,1,neb_size),neb_size^2,1);%81*9 
            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;%写入位置进行更改,每次写入81个数 
            fprintf(1,'进度:%f%% (%6d/%6d).\n',temp1/temp*100,temp1,temp); 
        end 
    end 
    vals=vals(1:len); 
    row_inds=row_inds(1:len); 
    col_inds=col_inds(1:len); 
    save(sprintf('(%s)_uint.mat',img_name),'vals','row_inds','col_inds'); 
else 
    load(sprintf('(%s)_uint.mat',img_name)); 
end 
 
 
%使用向量row_inds,col_inds建立一个img_size*img_size的稀疏矩阵 A(row_inds(k),col_inds(k)) = vals(k) 
if ~exist(sprintf('(%s)_A.mat',img_name),'file') 
    A=sparse(row_inds,col_inds,vals,img_size,img_size); 
    save(sprintf('(%s)_A.mat',img_name),'A'); 
else 
    load(sprintf('(%s)_A.mat',img_name)); 
end 
sumA=sum(A,2);%求稀疏矩阵的每一行的总和 
%spdiags(sumA(:),0,img_size,img_size);从矩阵sum(A)的每一列按照0偏移创建一个img_size*img_size大小稀疏存储的矩阵。 
%也就是对角线元素为sumA(:),一个系数存储的矩阵 
L=spdiags(sumA(:),0,img_size,img_size)-A; 
U=spdiags(win_b(:),0,img_size,img_size);% U is an identity matrix of the same size as L 
% Lambda=0.0001; 
Lambda=0.001; 
% The aurthor set a small value on λ (10?4 in our experiments) so that t is softly constrained by t?. 
if ~exist(sprintf('(%s)_x.mat',img_name),'file') 
    x=(L+Lambda*U)\(Lambda*win_b(:).*win_b(:));% ☆求解 
    save(sprintf('(%s)_x.mat',img_name),'x'); 
else 
    load(sprintf('(%s)_x.mat',img_name)); 
end 
 
%去掉0-1范围以外的数 
%--------------------------------------------------------------------------------------------- 
alpha=max(min( reshape(x,M,N),1 ),0.1);% △ alpha is the "refined transmission map"! 
figure(4), imshow(alpha,[]);title('Refined transmission map ( Soft Matting Technology Used)'); 
%A=200/255; 
A=AtomsphericLight; 
 
 
for i=1:P, 
    dehaze(:,:,i)=(IM(:,:,i)-A(i))./alpha+A(i); 
end 
imwrite(dehaze,'temp.png'); 
figure(5), imshow(dehaze,[]);title('修复的结果'); 
 
% the depth map 
DepthMap=log(alpha)./Beta*(-1); 
figure(6), imshow(DepthMap,[]);title('DepthMap');colormap('hot');