1 简介
为了改善图像融合效率,针对当前图像融合方法存在的局限性,提出一种基于小波变换的多聚焦图像融合方法.首先确定最优的小波分解层数,采用小波变换对图像进行多次分解;然后考虑低频子带与高频子带各自的特点,选择局部平均梯度准则作为高频子带的融合规则,低频子带采用3个系数的平均值作为融合规则;最后通过仿真实验对图像融合的有效性进行测试.实验结果表明,该方法获得了更理想的图像融合结果,提高了融合后的图像质量,且融合效果明显优于对比图像融合方法.
对二维图像进行 J 层的小波分解，可得 3J+1 个不同的频带，其中包含 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 , sx+1:sx+sx, :) = detailCoef2Weight( C( 1:sy , sx+1:sx+sx, : ), G);
W(sy+1:sy+sy, 1:sx , :) = detailCoef2Weight( C(sy+1:sy+sy, 1:sx , : ), G);
W(sy+1:sy+sy, sx+1:sx+sx, :) = detailCoef2Weight( C(sy+1:sy+sy, sx+1:sx+sx, : ), 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( iLv+1, 1 ); sx_c = S( iLv+1, 2 );
W( 1:sy_p, sx_p+1:sx_p+sx_p, : ) ...
= reweightDetailSubband( ...
W( 1:sy_p, sx_p+1:sx_p+sx_p, : ), ...
W( 1:sy_c, sx_c+1:sx_c+sx_c, : ) );
W( sy_p+1:sy_p+sy_p, 1:sx_p, : ) ...
= reweightDetailSubband( ...
W( sy_p+1:sy_p+sy_p, 1:sx_p, : ), ...
W( sy_c+1:sy_c+sy_c, 1:sx_c, : ) );
W( sy_p+1:sy_p+sy_p, sx_p+1:sx_p+sx_p, : ) ...
= reweightDetailSubband( ...
W( sy_p+1:sy_p+sy_p, sx_p+1:sx_p+sx_p, : ), ...
W( sy_c+1:sy_c+sy_c, sx_c+1:sx_c+sx_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 , sx+1:sx+sx, :), ...
W(sy+1:sy+sy, 1:sx , :), ...
W(sy+1:sy+sy, sx+1:sx+sx, :) );
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-y1+1; width:x2-x1+1;
%
% 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, hx+1:ex) = H;
C(hy+1:ey, 1:hx ) = V;
C(hy+1:ey, hx+1: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-y1+1; width:x2-x1+1;
%
% 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 , sx+1:sx2 );
V = C( sy+1:sy2, 1:sx );
D = C( sy+1:sy2, sx+1:sx2 );
A = idwt97( A, H, V, D );
end
A = A( 1:S( N+1, 1 ), 1:S( N+1, 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
end
3 仿真结果
4 参考文献
[1]管飚. "基于小波变换的多聚焦图像融合方法." 吉林大学学报：理学版 55.4(2017):6.
博主简介：擅长智能优化算法、神经网络预测、信号处理、元胞自动机、图像处理、路径规划、无人机等多种领域的Matlab仿真，相关matlab代码问题可私信交流。
部分理论引用网络文献，若有侵权联系博主删除。