1 简介为了改善图像融合效率,针对当前图像融合方法存在的局限性,提出一种基于小波变换的多聚焦图像融合方法.首先确定最优的小波分解层数,采用小波变换对图像进行多次分解;然后考虑低频子带与高频子带各自的特点,选择局部平均梯度准则作为高频子带的融合规则,低频子带采用3个系数的平均值作为融合规则;最后通过仿真实验对图像融合的有效性进行测试.实验结果表明,该方法获得了更理想的图像融合结果,提高了融合后的图像质量,且融合效果明显优于对比图像融合方法.对二维图像进行 J 层的小波分解可得 3J1 个不同的频带其中包含 3J 个高频带和一个低频带。并且能量主要集中于低频部分而高频部分代表了图像的细节信息。下面以两幅图像的融合为例说明基于小波变换的图像融合原理。如图 1 所示对源图像 1 和源图像 2 进行小波分解即用低通滤波器 L 和高通滤波器 H 分别对两源图像的水平方向和垂直方向进行滤波使源图像分解为水平低频垂直低频(LL)、水平低频垂直高频(LH)、水平高频垂直低频(HL)、水平高频垂直高频(HH)的 4 个子图像再根据需要利用低通滤波器 L 和高通滤波器 H 对LL 子图像重复上面的过程这样就建立各图像的小波塔形分解。接着对分解后的低频子图像和高频子图像根据需要进行融合处理得到融合后的小波金字塔。最后对融合后的小波金字塔进行小波逆变换即可得到融合结果。​2 部分代码function J imedgefuse( para, varargin ) %IMEDGEFUSE composite images taken with different depth-of-fields % as a result a pan focus image is generated. % % Usage: % % J IMEDGEFUSE( para, fname, fname, ..., fname ) % J IMEDGEFUSE( para, Image, Image, ..., Image ) % % input arguments : % % para - parameters (the detail is described later) % image - gray scale images (R^MxN) and rgb color images (R^MxNx3) % fname - file names of input images % % output arguments : % % J - a fused image % % parameters : % % para.nLv 6 - the total number of decompositions. % para.propagate 1 - weight propagation among subbands from parent to children. % para.denoise 1 - smoothing at each subbands for noise reduction. % para.debug_mode 1 - show debug message. % % % % % % % parameters % global g_para; set_parameters( para ); g_para.nImg length( varargin ); % number of images nImg g_para.nImg; nLv g_para.nLv; % total number of decomposition % % Wavelet decomposition for the Y of YCbCr % DEBUG_MSG( Reading images and applying wavelet to them (Y of YCbCr) ); for k 1:nImg % RGB color images are converted into the YCrCb color. I my_imread( varargin{ k } ); if k 1 % initialization of the coefficient array [ey, ex, nLv] output_size( size(I,1), size(I,2), nLv ); C zeros(ey,ex,nImg); end [C(:,:,k), S] wavedec97( I(:,:,1), nLv ); DEBUG_MSG( . ); end DEBUG_MSG( OK\n ); W gen3DWeightMatrix( C, S ); F genNewSubbands( C, W ); J(:,:,1) waverec97( F, S ); if g_para.bColor 0 return end % % Wavelet decomposition for the Cb of YCbCr % DEBUG_MSG( Reading pictures and applying wavelet to them (Cb of YCbCr) ); for k 1:nImg I my_imread( varargin{ k } ); [C(:,:,k), S] wavedec97( I(:,:,2), nLv ); DEBUG_MSG( . ); end DEBUG_MSG( OK\n ); F genNewSubbands( C, W ); J(:,:,2) waverec97( F, S ); % % Wavelet decomposition for the Cr of YCbCr % DEBUG_MSG( Reading pictures and applying wavelet to them (Cr of YCbCr) ); for k 1:nImg I my_imread( varargin{ k } ); [C(:,:,k), S] wavedec97( I(:,:,3 ), nLv ); DEBUG_MSG( . ); end DEBUG_MSG( OK\n ); F genNewSubbands( C, W ); J(:,:,3) waverec97( F, S ); % % Final conversion % J ycbcr2rgb( J ); end % %-------------------------------------------------------------------------- % function set_parameters( para ) global g_para; g_para para; % initialization % total number of decompositions if ~isfield( g_para, nLv ) g_para.nLv 5; end % weight propagation among subbands from parent to children if ~isfield( g_para, propagate ) g_para.propagate 1; end % weight smoothing at each subband if ~isfield( g_para, denoise ) g_para.denoise 1; end % hide debug message if ~isfield( g_para, debug_mode ) g_para.debug_mode false(1); end g_para.bColor true(1); % color image end % %-------------------------------------------------------------------------- % function I my_imread( I ) % % read image from a file. % RGB color images are converted into the YCrCb color. % if ischar( I ) I imread( char( I ) ); end if size( I, 3 ) ~ 3 g_para.bColor false; end if ~isa( I, double ) I im2double( I ); end if size( I, 3 ) 3 I rgb2ycbcr( I ); end end % %-------------------------------------------------------------------------- % function W gen3DWeightMatrix( C, S ) % % C : wavelet coefficient image % S : wavelet sub-bands size % global g_para; DEBUG_MSG( Creating a weight matrix...\n ); nLv size(S, 1) - 1; % Weight matrix W zeros( size(C) ); % G is low pass filter for denoising G fspecial( gaussian, 3, 3 ); % almost box filter DEBUG_MSG( Converting a map of usage index to a ratio map ); for iLv nLv:-1:1 sy S( iLv, 1 ); sx S( iLv, 2 ); W( 1:sy , sx1:sxsx, :) detailCoef2Weight( C( 1:sy , sx1:sxsx, : ), G); W(sy1:sysy, 1:sx , :) detailCoef2Weight( C(sy1:sysy, 1:sx , : ), G); W(sy1:sysy, sx1:sxsx, :) detailCoef2Weight( C(sy1:sysy, sx1:sxsx, : ), G); DEBUG_MSG( . ); end DEBUG_MSG( OK\n ); % % Wavelet parents-to-child relationship % if ( g_para.propagate 1 ) DEBUG_MSG( Reweighting Dc sub-bands from children to parent ); for iLv nLv-1:-1:1 sy_p S( iLv, 1 ); sx_p S( iLv, 2 ); sy_c S( iLv1, 1 ); sx_c S( iLv1, 2 ); W( 1:sy_p, sx_p1:sx_psx_p, : ) ... reweightDetailSubband( ... W( 1:sy_p, sx_p1:sx_psx_p, : ), ... W( 1:sy_c, sx_c1:sx_csx_c, : ) ); W( sy_p1:sy_psy_p, 1:sx_p, : ) ... reweightDetailSubband( ... W( sy_p1:sy_psy_p, 1:sx_p, : ), ... W( sy_c1:sy_csy_c, 1:sx_c, : ) ); W( sy_p1:sy_psy_p, sx_p1:sx_psx_p, : ) ... reweightDetailSubband( ... W( sy_p1:sy_psy_p, sx_p1:sx_psx_p, : ), ... W( sy_c1:sy_csy_c, sx_c1:sx_csx_c, : ) ); DEBUG_MSG( . ); end DEBUG_MSG( OK\n ); end % % Create the weight map of the approximation coefficient sub-band % DEBUG_MSG( Creating the weight map of the approximation sub-band ); sy S(1,1); sx S(1,2); W(1:sy, 1:sx, :) reweightApproxSubband( ... W( 1:sy , sx1:sxsx, :), ... W(sy1:sysy, 1:sx , :), ... W(sy1:sysy, sx1:sxsx, :) ); DEBUG_MSG( OK\n ); end % %-------------------------------------------------------------------------- % function W detailCoef2Weight( C, G ) % % G is low pass filter % global g_para; [~, K] max(abs(C), [], 3); % convert index to usage ratio W zeros( size(C) ); for k 1:size( C, 3 ) W(:,:,k) (K k); end if ( g_para.denoise 1 ) W imfilter( W, G, symmetric ); end end % %-------------------------------------------------------------------------- % function Wp reweightDetailSubband( Wp, Wc ) % % Deblur detail sub-bands by passing up and conversing the weight vectors % global g_para; nImg g_para.nImg; Wc Wc(1:2:end, :, :) Wc(2:2:end, :, :); % ex 1 2 3 4 3 7 Wc Wc(:, 1:2:end, :) Wc(:, 1:2:end, :); Wc 0.25 * Wc; Wp Wp .* Wc; W_sum sum( Wp, 3 ); Mask W_sum ~ 0; W_sum W_sum ~Mask; % interpolate 1 to avoid showing divide by 0 W_sum repmat( W_sum, [1,1,nImg] ); Mask repmat( Mask, [1,1,nImg] ); Wp Mask .* ( Wp ./ W_sum ) ~Mask .* Wc; end % %-------------------------------------------------------------------------- % function Wa reweightApproxSubband( Wh, Wv, Wd ) % % Deblur the approximation sub-band by passing up % and conversing the weight vectors % Wa (1/3) * (Wh Wv Wd); end % %-------------------------------------------------------------------------- % function F genNewSubbands( C, W ) % % Make new sub-bands filled up with focused coefficient % F sum(C.*W, 3); end % %-------------------------------------------------------------------------- % % Funcionts for wavelet transform % %-------------------------------------------------------------------------- % function [ey, ex, n] output_size( sy, sx, n ) % This function gives the size of output image obtained from % the specified decomopsition level. n_max floor( log2( min(sy,sx) ) ) - 1; if (n n_max) n n_max; DEBUG_MSG( \nWarning! : Specified Lv is higher than the decomposable Lv.\n ); DEBUG_MSG( The max Lv %d is used instead.\n, n ); end n2 2^n; % Extend image symmetrically ey n2 * floor( ( sy-1 )/n2 ) n2; ex n2 * floor( ( sx-1 )/n2 ) n2; end % %-------------------------------------------------------------------------- % function [C, S] wavedec97(A, n) %WAVEDEC97 performs the multi stage 2D wavelet decomposition % using lifting scheme. % % [C, S] WAVEDEC97(A, N) % % input argument % - A : Approximation coefficients ( Image ) % - N : Number of decomposition level % % output argument % - C : Coefficient Image % - S : Size of the Image ( margin top and left of each subbands in the whole image ) % % ex. S [y1, x1, y2, x2] % margin-top:y1; margin-left:x1; height:y2-y11; width:x2-x11; % % also see WAVEREC97 % % size check [sy,sx] size( A ); [ey,ex,n] output_size( sy, sx, n ); % image padding dy ey - sy; dx ex - sx; pady 1:sy; padx 1:sx; if ( dy 0 ) pady [pady, fliplr( pady( end-dy:end-1 ) )]; end if ( dx 0 ) padx [padx, fliplr( padx( end-dx:end-1 ) )]; end if ( dy 0 || dx 0 ) A A( pady, padx ); end % Initialize output arguments C zeros(ey,ex); S [sy,sx]; % Aterative decomposition for iL 1:n [A, H, V, D] dwt97( A ); hy ey * 0.5; hx ex * 0.5; C(1:hy, hx1:ex) H; C(hy1:ey, 1:hx ) V; C(hy1:ey, hx1:ex) D; S [hy,hx; S]; ey hy; ex hx; end C( 1:hy, 1:hx ) A; end % %-------------------------------------------------------------------------- % function [A, H, V, D] dwt97( A ) l [-1.58613432, -0.05298011854, 0.8829110762, 0.4435068522; -1.58613432, -0.05298011854, 0.8829110762, 0.4435068522]; s 1.149604398; % scale factor [L, H] filterdownl( A, l, s ); [V, D] filterdownl( H, l, s ); D D; V V; [A, H] filterdownl( L, l, s ); H H; A A; end % %-------------------------------------------------------------------------- % function [A, D] filterdownl( X, l, s ) [sy2, sx] size( X ); sy 0.5*sy2; A X(1:2:end,:); D X(2:2:end,:); pado [1:sy, sy]; pade [1, 1:sy]; D D filter2( l(:,1), A(pado,:), valid ); A A filter2( l(:,2), D(pade,:), valid ); D D filter2( l(:,3), A(pado,:), valid ); A A filter2( l(:,4), D(pade,:), valid ); A A * s; D D *(1/s); end % %-------------------------------------------------------------------------- % function A waverec97( C, S ) %WAVEREC97 performs the multi stage 2 dimensional wavelet reconstruction % using lifting scheme. % A waverec97( C, S ) % % input argument % - C : Coefficient Image % - S : Size of the Image (margin-top and left of each subbands in the whole image) % % ex. S [y1, x1, y2, x2 ] % margin-top:y1; margin-left:x1; height:y2-y11; width:x2-x11; % % output argument % - A : Recounstruction Image % % also see WAVEDEC97 N length( S )-1; A C( 1:S(1,1), 1:S(1,2) ); for i 1:N sy S(i,1); sx S(i,2); sy2 sy sy; sx2 sx sx; H C( 1:sy , sx1:sx2 ); V C( sy1:sy2, 1:sx ); D C( sy1:sy2, sx1:sx2 ); A idwt97( A, H, V, D ); end A A( 1:S( N1, 1 ), 1:S( N1, 2 ) ); end % %-------------------------------------------------------------------------- % function A idwt97( A, H, V, D ) l [-1.58613432, -0.05298011854, 0.8829110762, 0.4435068522; -1.58613432, -0.05298011854, 0.8829110762, 0.4435068522]; s 1.149604398; L filterupl( A, H, l, s ); H filterupl( V, D, l, s ); A filterupl( L, H, l, s ); end % %-------------------------------------------------------------------------- % function X filterupl( A, D, l, s ) [sy,sx] size( A ); sy2 sy * 2; pado [1:sy, sy]; pade [1, 1:sy]; D s * D; A (1/s) * A; A A - filter2( l(:,4), D(pade,:), valid ); D D - filter2( l(:,3), A(pado,:), valid ); A A - filter2( l(:,2), D(pade,:), valid ); D D - filter2( l(:,1), A(pado,:), valid ); X( [1:2:sy2, 2:2:sy2], : ) [A; D]; end % %-------------------------------------------------------------------------- % function DEBUG_MSG( format, varargin ) % % show message when debug mode % global g_para; if g_para.debug_mode fprintf( format, varargin{:} ); end end3 仿真结果4 参考文献[1]管飚. 基于小波变换的多聚焦图像融合方法. 吉林大学学报理学版 55.4(2017):6.博主简介擅长智能优化算法、神经网络预测、信号处理、元胞自动机、图像处理、路径规划、无人机等多种领域的Matlab仿真相关matlab代码问题可私信交流。部分理论引用网络文献若有侵权联系博主删除。