MATLAB implementation example of wavelet image fusion (add image demo 080428)

A few days ago, we discussed a Wavelet Image Fusion Algorithm Based on Low Frequency Fusion Strategy (http://blog.csdn.net/chenyusiyuan/archive/2007/11/14/1883987.aspx). Today we will talk about several details of Algorithm Programming implementation.

(1) high-frequency band convergence
First, high-frequency coefficient fusion adopts the rule of getting a large value based on the absolute value of the pixel. Note that the object to be compared is the absolute value of the pixel of the two original images, the value of the pixel in the merged image is the value of the point with the larger absolute value of the two original images, that is, the actual value of the point with the larger absolute value. The corresponding MATLAB program is:

% I1, I2 is two original images
Y1 = mywavedec2 (I1, dim );
Y2 = mywavedec2 (I2, dim );

[R, C] = size (Y1); % Image Fusion Based on the Low-Frequency Fusion Algorithm
For I = 1: R % first, take the value of the absolute wavelet decomposition coefficient corresponding to the two source images as the decomposition coefficient of the fusion images.
For j = 1: c
If (ABS (Y1 (I, j)> = ABS (Y2 (I, j )))
Y3 (I, j) = Y1 (I, j );
Elseif (ABS (Y1 (I, j) <ABS (Y2 (I, j )))
Y3 (I, j) = Y2 (I, j );
End
% Y3 (I, j) = max (ABS (Y1 (I, j), ABS (Y2 (I, j); % this line of code is incorrect, take the absolute value rather than the actual value
End
End

(2) objects of Low Frequency Fusion
The object of Low-Frequency Fusion is the LL-N sub-image of the N-level decomposed image, the corresponding program code is:

LLA = Y1 (1: R/(2 ^ dim), 1: C/(2 ^ dim); % call the lowfrefus function to combine the wavelet factorization coefficients of the low-frequency part.
LLB = Y2 (1: R/(2 ^ dim), 1: C/(2 ^ dim ));
Y3 (1: R/(2 ^ dim), 1: C/(2 ^ dim) = lowfrefus (LLA, LLB );

(3) program code of the Low-Frequency Fusion Policy
The implementation of this part of the program code is not very difficult, according to the formula listed in the algorithm principle can quickly write the corresponding code. However, we have to deal with some details, for example, determining the region size, the relationship between the region boundary and the image boundary, determining the weights of each point in the region center and the region, and the specific location of the region center in the original image.

Function Y = lowfrefus (A, B );
% Function Y = lowfrefus (A, B) on the input two wavelet decomposition coefficient matrices, according to the low frequency fusion algorithm, the low frequency wavelet decomposition coefficient of the fusion image is obtained.

[Row, Col] = size (a); % calculate the number of rows and columns in the factorization Coefficient Matrix
Alpha = 0.5; % Alpha is the threshold for variance matching comparison
For I = 1: row % based on the low-frequency fusion algorithm, first obtain the regional variance and variance matching degree centered on point P in matrix A and matrix B.
For j = 1: Col % and then determine the wavelet decomposition coefficient of the fusion image based on the comparison of variance matching degree and threshold value.
[M2p (I, j), GA (I, j), GB (I, j)] = area_var_match (a, B, [I, j]);
Wmin = 0.5-0.5 * (1-m2p (I, j)/(1-alpha ));
Wmax = 1-wmin;
If m2p (I, j) <Alpha % m2p indicates variance matching degree
If GA (I, j)> = GB (I, j) % if the matching degree is smaller than the threshold value, the decomposition coefficient of the corresponding point with a large regional variance is used as the decomposition coefficient of the fusion image.
Y (I, j) = a (I, j );
Else
Y (I, j) = B (I, j );
End
Else % if the matching degree is greater than the threshold value, the corresponding decomposition coefficient is obtained using the weighted average method.
If GA (I, j)> = GB (I, j)
Y (I, j) = wmax * a (I, j) + wmin * B (I, j );
Else
Y (I, j) = wmin * a (I, j) + wmax * B (I, j );
End
End
End
End



Function W = weivec (x, P );
% Function W = weivec (X, p) calculates the corresponding weight of each point in the matrix when the vertex P is centered on the input R * C matrix.
% The closer the distance from point P, the larger the weight. Weights are obtained by adding the Gaussian distribution of rows and columns in a weighted manner.

[R, C] = size (X );
P1 = P (1); P2 = P (2 );
Sig = 1;
For I = 1: R
For j = 1: c
W (I, j) = 0.5 * (gaussmf (I, [sig P1]) + gaussmf (J, [sig P2]);
End
End



Function [m2p, Ga, GB] = area_var_match (A, B, P );
The % function area_var_match calculates the regional variance and the regional variance matching degree of the two input matrices centered on point P.

Level = 1; % set the area size
[Suba, MPa, NPA] = submat (A, P, level); % The submat function takes the vertex P as the center of the input matrix and the order number (2 * level + 1) as a child Matrix
[Subb, MPB, NPB] = submat (B, P, level );

[R, C] = size (Suba );
W = weivec (Suba, [MPa NPA]); % obtain the weight distribution of the submatrix

Avera = sum (Suba)/(R * C); % calculate the average value of the sub-Matrix
Averb = sum (subb)/(R * C );

GA = sum (W. * (Suba-avera). ^ 2); % calculate the regional variance of the submatrix
GB = sum (W. * (subb-averb). ^ 2 ));
If (GA = 0) & (GB = 0) % calculate the regional variance matching degree of the two submatrices
M2p = 0;
Else
M2p = 2 * sum (W. * ABS (Suba-avera). * ABS (subb-averb)/(GA + GB );
End



Function [SMAT, MP, NP] = submat (x, P, level );
% Function submat: Take the input matrix with the dot P as the center and the square matrix with the order (2 * level + 1) as the output sub-matrix.

[Row, Col] = size (X );
M = P (1); n = P (2 );

If (M> row) | (n> col)
Error ('point P is out of matrix X! ');
Return;
End
If (2 * level + 1)> row) | (2 * level + 1)> col)
Error ('too large sample area level! ');
Return;
End

% Set the submatrix Boundary Value
Up = m-level; down = m + level;
Left = N-level; Right = N + level;
% If a Submatrix boundary value exceeds the corresponding boundary of the input matrix, the boundary is processed,
% Refers to translation in the opposite direction after the boundary is exceeded, so that it exactly overlaps with the boundary
If left <1
Right = right + 1-left;
Left = 1;
End
If right> col
Left = left + col-right;
Right = Col;
End
If up <1
Down = down + 1-up;
Up = 1;
End
If down> row
Up = up + row-down;
Down = row;
End

% Obtain the child matrix as the output, and calculate the position of vertex P in the Child matrix of the output.
SMAT = x (up: down, left: right );
MP = m-up + 1; Np = N-left + 1;

The following is an example of a specific image. Figure 1 is an image with a clear vision and a blurred vision. Figure 2 is an image with a clear vision, and figure 3 is a fusion.