[image fusion] sparse regularization based on non convex penalty to realize image fusion with matlab code

1 brief introduction

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 Part of the code

%Code of image fusion with GMC regularization.getd = @(p)path(p,path);getd('Metrics/');addpath('./utilities');​resultDir = '';seq = 'clocks';​%original image dataz2colour = [];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);endz2 = 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 transformlevels = 4;[a,L] = wavedec2(z1,levels,'Haar');​Wt  = @(x) wavedec2(x,levels,'Haar');Wtt = @(x) waverec2(x,L,'Haar');​%Parameter settingfor 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;   endfiguresubplot(131)imshow(z1);title(' chart 1')subplot(132)​imshow(z2);;title(' chart 2')subplot(133)​imshow(final_old);;title(' Fusion map ')

3 Simulation results

4 reference

[1] Meijinjin . Research on image restoration and fusion based on regularization method [D]. University of electronic technology , 2017.

[2] Anantrasirichai N ,  Zheng R ,  Selesnick I , et al. Image Fusion via Sparse Regularization with Non-Convex Penalties[J].  2019.

About bloggers ： Good at intelligent optimization algorithms 、 Neural networks predict 、 signal processing 、 Cellular automata 、 The image processing 、 Path planning 、 UAV and other fields Matlab Simulation , relevant matlab Code problems can be exchanged by private letter .

Some theories cite network literature , If there is infringement, contact the blogger to delete .