边缘保持滤波 matlab,基于L0边缘保持滤波的图像融合 matlab代码

本文程序是根据"Zhao J, Feng H, Xu Z, et al. Detail enhanced multi-source fusion using visual weight map extraction based on multi scale edge preserving decomposition[J]. Optics Communications, 2013, 287: 45-52."编写的，最大的问题是运行时间太长了，有五六分钟，还有就是最后的融合效果并不是很理想。

具体matlab程序为：

clear all;

%path_A = '..\IR and Vis\Camp_Vis.jpg'; path_B = '..\IR and Vis\Camp_IR.jpg';

%path_A = '..\IR and Vis\Dune_Vis.jpg'; path_B = '..\IR and Vis\Dune_IR.jpg';

path_A = '..\Medical Images\source1_1.tif'; path_B = '..\Medical Images\source1_2.tif';

img1 = double(imread(path_A))/255.0;

img2 = double(imread(path_B))/255.0;

if(size(img1,3)~=1)

img1=rgb2gray(img1);

end

if(size(img2,3)~=1)

img2=rgb2gray(img2);

end

%% ---------- L0 Decomposition --------------

nLevel = 3;

lambda=[0.0005 0.003 0.018];

f1 = cell(1, nLevel+1);

d1 = cell(1, nLevel+1);

f1{1}=img1;

for L = 2:nLevel+1,

f1{L} = L0Smoothing(f1{L-1},lambda(L-1));

d1{L} = f1{L-1} - f1{L};

end

% figure;imshow([mat2gray(f1{2}),mat2gray(f1{3}),mat2gray(f1{4});...

% mat2gray(d1{2}),mat2gray(d1{3}),mat2gray(d1{4})]);

f2 = cell(1, nLevel+1);

d2 = cell(1,nLevel+1);

f2{1}=img2;

for L = 2:nLevel+1,

f2{L} = L0Smoothing(f2{L-1},lambda(L-1));

d2{L} = f2{L-1} - f2{L};

end

% figure;imshow([mat2gray(f2{2}),mat2gray(f2{3}),mat2gray(f2{4});...

% mat2gray(d2{2}),mat2gray(d2{3}),mat2gray(d2{4})]);

%% ---------- Visual Weight Map(VWM) --------------

V_d1 = cell(1,nLevel+1);

V_d2 = cell(1,nLevel+1);

V_f1=zeros(size(img1));

V_f2=zeros(size(img2));

%VWM of img1 detail layer

for L = 2:nLevel+1,

V_d1{L}=zeros(size(img1));

img1_detail=uint8(mat2gray(d1{L})*255);

for j=0:255,

for i=0:255,

V_d1{L}(img1_detail==j)=V_d1{L}(img1_detail==j)+length(find(img1_detail==i))*abs(j-i);

end

end

V_d1{L}=mat2gray(V_d1{L});

end

%VWM of img2 detail layer

for L = 2:nLevel+1,

V_d2{L}=zeros(size(img2));

img2_detail=uint8(mat2gray(d2{L})*255);

for j=0:255,

for i=0:255,

V_d2{L}(img2_detail==j)=V_d2{L}(img2_detail==j)+length(find(img2_detail==i))*abs(j-i);

end

end

V_d2{L}=mat2gray(V_d2{L});

end

%VWM of img1 base layer

img1_base=uint8(255*f1{nLevel+1});

for j=0:255,

for i=0:255,

V_f1(img1_base==j)=V_f1(img1_base==j)+length(find(img1_base==i))*abs(j-i);

end

end

V_f1=mat2gray(V_f1);

%VWM of img2 base layer

img2_base=uint8(255*f2{nLevel+1});

for j=0:255,

for i=0:255,

V_f2(img2_base==j)=V_f2(img2_base==j)+length(find(img2_base==i))*abs(j-i);

end

end

V_f2=mat2gray(V_f2);

% figure;imshow(V_f2);

%% ---------- Multi-scale Combination --------------

% base layer combination

F_base=0.5*((f1{nLevel+1}.*V_f1+f2{nLevel+1}.*(1-V_f1))+(f1{nLevel+1}.*V_f2+f2{nLevel+1}.*(1-V_f2)));

% detail layer combination

F=cell(1,nLevel+1);

for L = 2:nLevel+1,

F{L}=0.5*((d1{L}.*V_d1{L}+d2{L}.*(1-V_d1{L}))+(d1{L}.*V_d2{L}+d2{L}.*(1-V_d2{L})));

end

%% ---------- Fusion Result --------------

beta=[0.8 0.6 0.6 0.5];

F_out=zeros(size(img1));

for L = 2:nLevel+1,

F_out=F_out+beta(L-1)*F{L};

if L==nLevel+1,

F_out=F_out+beta(L)*F_base;

end

end

figure;imshow([img1,img2]);

figure;imshow(F_out,[]);

其中用到的L0Smoothing()函数如下：

% Distribution code Version 1.0 -- 09/23/2011 by Jiaya Jia Copyright 2011, The Chinese University of Hong Kong.

%

% The Code is created based on the method described in the following paper

% [1] "Image Smoothing via L0 Gradient Minimization", Li Xu, Cewu Lu, Yi Xu, Jiaya Jia, ACM Transactions on Graphics,

% (SIGGRAPH Asia 2011), 2011.

%

% The code and the algorithm are for non-comercial use only.

function S = L0Smoothing(Im, lambda, kappa)

% L0Smooth - Image Smoothing via L0 Gradient Minimization

% S = L0Smooth(Im, lambda, kappa) performs L0 graidient smoothing of input

% image Im, with smoothness weight lambda and rate kappa.

%

% Paras:

% @Im : Input UINT8 image, both grayscale and color images are acceptable.

% @lambda: Smoothing parameter controlling the degree of smooth. (See [1])

% Typically it is within the range [1e-3, 1e-1], 2e-2 by default.

% @kappa : Parameter that controls the rate. (See [1])

% Small kappa results in more iteratioins and with sharper edges.

% We select kappa in (1, 2].

% kappa = 2 is suggested for natural images.

%

% Example

% ==========

% Im = imread('pflower.jpg');

% S = L0Smooth(Im); % Default Parameters (lambda = 2e-2, kappa = 2)

% figure, imshow(Im), figure, imshow(S);

if ~exist('kappa','var')

kappa = 2.0;

end

if ~exist('lambda','var')

lambda = 2e-2;

end

S = im2double(Im);

betamax = 1e5;

fx = [1, -1];

fy = [1; -1];

[N,M,D] = size(Im);

sizeI2D = [N,M];

otfFx = psf2otf(fx,sizeI2D);

otfFy = psf2otf(fy,sizeI2D);

Normin1 = fft2(S);

Denormin2 = abs(otfFx).^2 + abs(otfFy ).^2;

if D>1

Denormin2 = repmat(Denormin2,[1,1,D]);

end

beta = 2*lambda;

while beta < betamax

Denormin = 1 + beta*Denormin2;

% h-v subproblem

h = [diff(S,1,2), S(:,1,:) - S(:,end,:)];

v = [diff(S,1,1); S(1,:,:) - S(end,:,:)];

if D==1

t = (h.^2+v.^2)

else

t = sum((h.^2+v.^2),3)

t = repmat(t,[1,1,D]);

end

h(t)=0; v(t)=0;

% S subproblem

Normin2 = [h(:,end,:) - h(:, 1,:), -diff(h,1,2)];

Normin2 = Normin2 + [v(end,:,:) - v(1, :,:); -diff(v,1,1)];

FS = (Normin1 + beta*fft2(Normin2))./Denormin;

S = real(ifft2(FS));

beta = beta*kappa;

%fprintf('.');

end

%fprintf('\n');

end

下面是效果图：