1 简介
The L1 norm regularized least squares method is often used for finding sparse approximate solutions and is widely used in 1-D signal restoration. Basis pursuit denoising (BPD) performs noise reduction in this way. However, the shortcoming of using L1 norm regularization is the underestimation of the true solution. Recently, a class of non-convex penalties have been proposed to improve this situation. This kind of penalty function is non-convex itself, but preserves the convexity property of the whole cost function. This approach has been confirmed to offer good performance in 1-D signal denoising. This paper demonstrates the aforementioned method to 2-D signals (images) and applies it to multisensor image fusion. The problem is posed as an inverse one and a corresponding cost function is judiciously designed to include two data attachment terms. The whole cost function is proved to be convex upon suitably choosing the non-convex penalty, so that the cost function minimization can be tackled by convex optimization approaches, which comprise simple computations. The performance of the proposed method is benchmarked against a number of state-of-the-art image fusion techniques and superior performance is demonstrated both visually and in terms of various assessment measures.
2 部分代码
%Code of image fusion with GMC regularization.
getd = @(p)path(p,path);
getd('Metrics/');
addpath('./utilities');
resultDir = '';
seq = 'clocks';
%original image data
z2colour = [];
if strcmp(seq, 'clocks')
z1 = im2double(imread('clock_A.bmp'));
z2 = im2double(imread('clock_B.bmp'));
elseif strcmp(seq,'Bild2')
commondir = '/Users/eexna/Work/Ultrasound imaging/Fusion/data/bild/';
z1 = im2double(imread([commondir,'Bild2_ir_wh.bmp']));
z2 = im2double(imread([commondir,'Bild2_tv.bmp']));
z2 = rgb2ycbcr(z2);
z2colour = z2;
elseif strcmp(seq, 'quad')
commondir = '/Users/eexna/Work/Ultrasound imaging/Fusion/data/junction/';
z1 = im2double(imread([commondir,'quad_TV_bw_user7_im3.jpg']));
z2 = im2double(imread([commondir,'quad_IR_user4_im5.jpg']));
elseif strcmp(seq, 'Octec')
commondir = '/Users/eexna/Work/Ultrasound imaging/Fusion/data/Octec/';
z1 = im2double(imread([commondir,'Octec_IR_22.jpg']));
z2 = im2double(imread([commondir,'Octec_TV_22.jpg']));
z2 = rgb2ycbcr(z2);
z2colour = z2;
elseif strcmp(seq, 'slika')
z1 = im2double(imread('/Users/eexna/Downloads/TestingImageDataset/testna_slika2a.bmp'));
z2 = im2double(imread('/Users/eexna/Downloads/TestingImageDataset/testna_slika2b.bmp'));
z1 = imresize(z1, [256 256]*2);
z2 = imresize(z2, [256 256]*2);
end
z2 = z2(:,:,1);
z1 = z1(:,:,1);
A = uint8(z1*255);
B = uint8(z2*255);
curreseultDir = [resultDir, seq,'/'];
mkdir(curreseultDir);
% z1 = imresize(z1, [512 512]);
% z2 = imresize(z2, [512 512]);
% z1 = im2double(imread('./original image data/ct.jpg'));
% z1=[z1,zeros(159,1)];
% z1=[z1;zeros(1,160)];
% z2 = im2double(imread('./original image data/mri.jpg'));
% z2=[z2,zeros(159,1)];
% z2=[z2;zeros(1,160)];
% z1 = imread('./original image data/011_CT_1.bmp');
% z2 = imread('./original image data/011_T2.bmp');
% z1= rgb2gray(z1);
% z2= rgb2gray(z2);
% z1=im2double(z1);
% z2=im2double(z2);
% z1 = im2double(imread('./original image data/1_1_f.tif'));
% z1=z1(1:600,1:600);
% z2 = imread('./original image data/1_1_o.tif');
% z2= rgb2gray(z2);
% z2=im2double(z2);
% z2=z2(1:600,1:600);
% Two sensor gain estimation methods
% [recovbete1 recovbete2] = correctsensor(z1, z2);
[recovbete1, recovbete2] = SensorGain(z1, z2);
% Wavelet and inverse wavelet transform
levels = 4;
[a,L] = wavedec2(z1,levels,'Haar');
Wt = @(x) wavedec2(x,levels,'Haar');
Wtt = @(x) waverec2(x,L,'Haar');
%Parameter setting
for lam=0.005 % [ 0.001 0.01 0.05]
tic;
rho=1;
gam=0.8;
miu=1.9 / ( rho * max( 1, gam / (1-gam) ) );
alp=lam*miu;
thk=zeros(1,size(a,2));
v=zeros(1,size(a,2));
%main GMC regularization algorithm
itrs=25;
for j=1:itrs
w=thk-miu*(Wt(recovbete2.*(recovbete2.*Wtt(thk)-z1))+Wt(recovbete1.*(recovbete1.*Wtt(thk)-z2))+gam*Wt(Wtt(v-thk)));
u=v-miu*gam*Wt(Wtt(v-thk));
thk=(max(abs(w)-alp,0)).*sign(w);
% normalise
xtemp = Wtt(thk);
xtemp = (xtemp-min(xtemp(:)))/range(xtemp(:));
thk = Wt(xtemp);
v=(max(abs(u)-alp,0)).*sign(u);
end
final=Wtt(thk);
final = final/max(final(:));%*(0.5*(max(z1(:))+max(z2(:))));
toc
final_old = (final-min(final(:)))/range(final(:));
peval = imqmet(cat(3, A, B), uint8(255*final_old));
cvval = Cvejic_metric(cat(3, A, B),uint8(255*final_old));
pvtrovic = petmetric(cat(3, A, B), uint8(255*final_old));
disp(num2str([peval cvval pvtrovic]))
peval = imqmet(cat(3, A, B), uint8(255*final));
cvval = Cvejic_metric(cat(3, A, B),uint8(255*final));
pvtrovic = petmetric(cat(3, A, B), uint8(255*final));
disp(num2str([peval cvval pvtrovic]))
if ~isempty(z2colour)
final_old(:,:,2:3) = z2colour(:,:,2:3);
final_old = ycbcr2rgb(final_old);
end
imwrite(final_old,[curreseultDir, seq, 'lam',num2str(lam),'.png']);
final_old = final;
end
figure
subplot(131)
imshow(z1);title('图1')
subplot(132)
imshow(z2);;title('图2')
subplot(133)
imshow(final_old);;title('融合图')
3 仿真结果
4 参考文献
[1]梅金金. 基于正则化方法的图像复原与融合研究[D]. 电子科技大学, 2017.
[2] Anantrasirichai N , Zheng R , Selesnick I , et al. Image Fusion via Sparse Regularization with Non-Convex Penalties[J]. 2019.
博主简介:擅长智能优化算法、神经网络预测、信号处理、元胞自动机、图像处理、路径规划、无人机等多种领域的Matlab仿真,相关matlab代码问题可私信交流。
部分理论引用网络文献,若有侵权联系博主删除。