# Multifocus Image Fusion Using Biogeography-Based Optimization

Authors: Ping Zhang; Chun Fei; Zhenming Peng; Jianping Li; Hongyi Fan
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

For multifocus image fusion in spatial domain, sharper blocks from different source images are selected to fuse a new image. Block size significantly affects the fusion results and a fixed block size is not applicable in various multifocus images. In this paper, a novel multifocus image fusion algorithm using biogeography-based optimization is proposed to obtain the optimal block size. The sharper blocks of each source image are first selected by sum modified Laplacian and morphological filter to contain an initial fused image. Then, the proposed algorithm uses the migration and mutation operation of biogeography-based optimization to search the optimal block size according to the fitness function in respect of spatial frequency. The chaotic search is adopted during iteration to improve optimization precision. The final fused image is constructed based on the optimal block size. Experimental results demonstrate that the proposed algorithm has good quantitative and visual evaluations.


---

## Body

## 1. Introduction

Optical lenses with long focal lengths often suffer from the problem of limited depth of field. It is impossible to get an image that contains all relevant objects in focus. The objects only on the focus plane are sharpness and other objects in front of or behind the focus plane are blurred [1]. Multifocus fusion method which synthesizes multiple images of the same view point under different focal settings can be used to extend depth of field and obtain an all-in focus image. The fused image has more useful information of the view point and is more suitable for many applications than any individual images. Multifocus image fusion technology has played important roles in many fields such as target recognition, remote sensing, medical diagnosis, and military application [2].Many multifocus fusion algorithms have been proposed in recent years. Basically, these fusion algorithms can be categorized into two groups: spatial domain fusion and transform domain fusion [3]. For spatial domain fusion, a new image is fused by directly selecting different regions from source images. Firstly, source images are divided into nonoverlapping blocks. Then, the sharpness values of blocks are calculated based on different sharpness measure methods. Finally, the sharper blocks from different source images are selected to fuse a new image. The common algorithms in spatial domain include average, variance, energy of image gradient (EOG), sum modified Laplacian (SML), and spatial frequency (SF) [4]. Recently, many new algorithms in spatial domain have been proposed to improve efficiency, such as artificial neural network (ANN) [5], pulse coupled neural network (PCNN) [6], independent component analysis (ICA) [7], robust principal component analysis (RPCA) [8], and neighbor distance [9]. For transform domain fusion, a new image is generated with certain frequency transforms. Firstly, source images are converted into a transform domain to obtain the corresponding transform coefficients. Then, the transform coefficients are integrated together based on different fusion rules. Finally, the fused image is constructed by applying the inverse transform. The common algorithms in transform domain are based on pyramid transform, such as Laplacian pyramid (LAP), gradient pyramid (GRP), ratio of low-pass pyramid (RAP), and morphological pyramid (ROP) [10]. Some fusion algorithms based on wavelet transform are proposed which are generally superior to the fusion algorithms based on pyramid transform, such as discrete wavelet transform (DWT) [11], stationary wavelet transform (SWT) [12], and dual-tree complex wavelet transform (DTCWT) [13]. Recently, many new multiscale multiresolution transform algorithms have been widely proposed in image fusion, such as curvelet transform [14], contourlet transform [15], nonsubsampled contourlet transform (NSCT) [16], and nonsubsampled shearlet transform (NSST) [17].The transform domain algorithms have no block artifact, but they usually are complicated and time-consuming to implement. And the multiscale multiresolution transform algorithms are generally shift-variant and sensitive to noise. The spatial domain algorithms are simple to implement and have low complexity [18]. However, block artifacts inevitably exist because the block size is fixed in these algorithms. If the block size is too large, there are both in-focus part and out-of-focus part in the same block. The block’s sharpness may be incorrect and the out-of-focus part may be selected as the part of the fused image when considering the integrity of the segmented part. If the block size is too small, the fused result is sensitive to noise and the computational complexity is too high. Obviously, a fixed block size may not be applicable in various multifocus images.To solve the fixed block size problem, some image fusion algorithms based on optimization are proposed, such as genetic algorithm (GA) [19, 20], particle swarm optimization (PSO) [21], and differential evolution (DE) [22]. These fusion algorithms utilize the global optimization characteristics to obtain the optimal block size, which effectively suppress the block artifact and increase the performance of the fused images. However, the efficiency of these fusion algorithms depends largely on the performance of the optimization algorithms. GA, PSO, and DE are not so satisfactory because of low convergence rate and low optimization accuracy.A novel multifocus image fusion algorithm using biogeography-based optimization is proposed in this paper. Biogeography-based optimization (BBO) [23] is a new swarm intelligence optimization, which is proposed by Dr. Simon in 2008. BBO obtains the global optimum through its migration mechanisms and mutation operation. BBO has fast convergence rate and high search precision compared with GA, PSO, and DE [23]. These advantages make BBO solve more effectively complex optimization problem. It has been applied to image and video processing, such as image classification [24], image matching [25], image segmentation [26], image enhancement [27], and motion estimation for video coding [28]. In this paper, the proposed algorithm utilizes the biogeography-based optimization technique to find the best block size. Moreover, chaotic search is embedded to improve optimization accuracy. Experiments on various multifocus images demonstrate that the proposed algorithm has good performance in terms of quantitative and visual evaluations.The rest of this paper is organized as follows. In Section2, the biogeography-based optimization is briefly reviewed. Section 3 describes the proposed image fusion algorithm in detail. The experiment results and discussions are presented in Section 4. The conclusions are given in Section 5.


## 2. Biogeography-Based Optimization

Dr. Simon proposed biogeography-based optimization (BBO) in 2008 [23]. It is based on the mathematics of biogeography, which describes how species migrate from one island to another, how new species arise, and how species become extinct. BBO, which is similar to GA, PSO, and DE, is a stochastic algorithm to solve optimization problem. It has two important operations, one is migration and the other one is mutation. Species among neighboring islands share their information through migration. Individual species improve their diversity through mutation. Global optimum can be effectively obtained with migration and mutation. BBO has good performance because of its fast convergence rate and high search precision [23].


### 2.1. Migration

In BBO, each solution of optimization problem is regarded as an island. The feature of each solution is regarded as a suitability index variable (SIV) and the fitness value of each solution is regarded as its habitat suitability index (HSI). The higher the HIS of an island, the better the performance of the optimization problem. High HSI and low HSI use the emigration and immigration rates of each solution to probabilistically share information between islands, as in Figure1.Figure 1
Island migration rates versus HSI.The immigration rateλ and the emigration rate μ are the functions of the number of species in the island
(1)
λ
s
=
I
1
-
S
S
max
⁡
,
μ
s
=
E
S
S
max
⁡
,
where E is the maximum emigration rate, I is the maximum immigration rate, S is the number of species of the Sth individual, and S
max
⁡ is the maximum number of species. For simplicity, here, assume E
=
I
=
1; that is, λ
k
+
μ
k
=
E.The basic step of migration is as follows. Firstly, calculate the HSI values of all islands and sort them in descending order. Second, select one island that is needed to immigrate based on immigration rate and choose its adjacent islands based on emigration rate. Then, randomly select SIV value from the adjacent islands to replace the SIV value of that island. Recalculate the HSI values of all islands. The island with highest HSI value is the optimal solution.


### 2.2. Mutation

An island’s HSI may change suddenly due to apparently random events such as disease and natural catastrophes. In BBO, this phenomenon is called mutation. The mutation ratesm
s is determined as follows:
(2)
m
s
=
m
max
⁡
1
-
P
s
P
max
⁡
,
where m
max
⁡ is a mutation parameter and P
max
⁡
=
arg max
P
i, i
=
1
,
…
,
N
,
P
s is probability of species. The relationship between P
s and migration rate is shown as follows:
(3)
P
s
=
-
(
λ
s
+
μ
s
)
P
s
+
μ
s
+
1
P
s
+
1
S
=
0
-
(
λ
s
+
μ
s
)
P
s
+
λ
s
-
1
P
s
-
1
+
μ
s
+
1
P
s
+
1
1
≤
S
≤
S
max
⁡
-
1
-
(
λ
s
+
μ
s
)
P
s
+
λ
s
-
1
P
s
-
1
S
=
S
max
⁡
.This mutation approach makes low HSI solutions likely to mutate, which gives them a chance of improving. It also makes high HSI solutions likely to mutate, which gives them a chance of improving more than they already had. Thus, the BBO mutation strategy tends to increase diversity of islands.


### 2.3. Elitism Strategy

BBO incorporate elitism strategy in order to save the features of the habitat that has the best solution in the iterative process as with GA, PSO, and DE. This prevents the best solutions from being corrupted by immigration or being ruined by mutation. So even if the habitat with highest HSI is destroyed, BBO has saved it and can revert back to it if needed. The proposed algorithm in this paper retains two higher islands as elites. These elites are kept from this generation to the next generation.


## 2.1. Migration

In BBO, each solution of optimization problem is regarded as an island. The feature of each solution is regarded as a suitability index variable (SIV) and the fitness value of each solution is regarded as its habitat suitability index (HSI). The higher the HIS of an island, the better the performance of the optimization problem. High HSI and low HSI use the emigration and immigration rates of each solution to probabilistically share information between islands, as in Figure1.Figure 1
Island migration rates versus HSI.The immigration rateλ and the emigration rate μ are the functions of the number of species in the island
(1)
λ
s
=
I
1
-
S
S
max
⁡
,
μ
s
=
E
S
S
max
⁡
,
where E is the maximum emigration rate, I is the maximum immigration rate, S is the number of species of the Sth individual, and S
max
⁡ is the maximum number of species. For simplicity, here, assume E
=
I
=
1; that is, λ
k
+
μ
k
=
E.The basic step of migration is as follows. Firstly, calculate the HSI values of all islands and sort them in descending order. Second, select one island that is needed to immigrate based on immigration rate and choose its adjacent islands based on emigration rate. Then, randomly select SIV value from the adjacent islands to replace the SIV value of that island. Recalculate the HSI values of all islands. The island with highest HSI value is the optimal solution.


## 2.2. Mutation

An island’s HSI may change suddenly due to apparently random events such as disease and natural catastrophes. In BBO, this phenomenon is called mutation. The mutation ratesm
s is determined as follows:
(2)
m
s
=
m
max
⁡
1
-
P
s
P
max
⁡
,
where m
max
⁡ is a mutation parameter and P
max
⁡
=
arg max
P
i, i
=
1
,
…
,
N
,
P
s is probability of species. The relationship between P
s and migration rate is shown as follows:
(3)
P
s
=
-
(
λ
s
+
μ
s
)
P
s
+
μ
s
+
1
P
s
+
1
S
=
0
-
(
λ
s
+
μ
s
)
P
s
+
λ
s
-
1
P
s
-
1
+
μ
s
+
1
P
s
+
1
1
≤
S
≤
S
max
⁡
-
1
-
(
λ
s
+
μ
s
)
P
s
+
λ
s
-
1
P
s
-
1
S
=
S
max
⁡
.This mutation approach makes low HSI solutions likely to mutate, which gives them a chance of improving. It also makes high HSI solutions likely to mutate, which gives them a chance of improving more than they already had. Thus, the BBO mutation strategy tends to increase diversity of islands.


## 2.3. Elitism Strategy

BBO incorporate elitism strategy in order to save the features of the habitat that has the best solution in the iterative process as with GA, PSO, and DE. This prevents the best solutions from being corrupted by immigration or being ruined by mutation. So even if the habitat with highest HSI is destroyed, BBO has saved it and can revert back to it if needed. The proposed algorithm in this paper retains two higher islands as elites. These elites are kept from this generation to the next generation.


## 3. Multifocus Image Fusion Using Biogeography-Based Optimization

In this paper, a novel multifocus image fusion algorithm using biogeography-based optimization is proposed. Considering the characteristic of image fusion, the island, SIV, and HSI of BBO is regarded as block size, width and height of block size, and image quality assessment, respectively. So, the dimension of optimization problem is only two, the width and height of block size. Here are the key steps of the proposed algorithm. Firstly, it selects the random block size as initial islands for utilizing the global random characteristics of BBO. Then, it calculates the HSI value of each block size. During iteration process, the width and height of the block size are updated through the migration and mutation operations. The iteration stops if the termination condition is met. Finally, the block size with the highest HSI is the optimal block size. The details of the multifocus image fusion using BBO are as follows.


### 3.1. HSI Function Selection

Spatial frequency [4] can measure the overall activity level in images and reflect the ability to express small details in images. The larger the value of spatial frequency in fused image, the better the performance of fusion. In the image fusion algorithms based on GA, PSO, and DE, spatial frequency is chosen as the fitness function [20–22]. It is a good fitness function to assess the fused image quality. In this paper, we also choose the spatial frequency as the HSI function of BBO
(4)
HSI
x
,
y
=
SF
x
,
y
.Spatial frequencySF
(
x
,
y
) of M
×
N size image is defined as follows:
(5)
SF
x
,
y
=
F
R
2
x
,
y
+
F
C
2
x
,
y
,
where
(6)
F
R
x
,
y
=
1
M
×
N
∑
x
=
1
M
∑
y
=
1
N
I
x
,
y
+
1
-
I
x
,
y
2
,
F
C
(
x
,
y
)
=
1
M
×
N
∑
x
=
1
M
∑
y
=
1
N
I
x
+
1
,
y
-
I
x
,
y
2
are the row and column gradients, respectively. I
(
x
,
y
) is the pixel of the image.


### 3.2. Sharpness Measurement

For multifocus image, the sharpness areas of different source images are selected to construct a fused image. How to define the sharpness areas is very important. There are many typical sharpness measurements, such as EOG, SF, and SML [18]. They measure the variation of pixels of blocks in source images to determine the sharpness. The blocks with greater values are considered as in-focus blocks. Experimental results of Huang and Jing [18] show that SML method provides better performance than EOG and SF. It differs from the usual Laplacian operator in that the absolute values of the partial second derivatives are summed instead of their actual values. The SML is defined as follows:
(7)
SML
I
x
,
y
=
∑
w
=
-
W
W
∑
h
=
-
H
H
L
I
x
+
w
,
y
+
h
2
,
where I
(
x
,
y
) is the pixel of the image. W and H is the window size. L
(
I
(
x
,
y
)
) is the modified Laplacian given as follows. Here, step is a variable that is used to set the distance between the central pixel and the pixels are used to compute the second order derivative. In this paper, a value of step
=
1 is found to produce the best results
(8)
L
I
x
,
y
=
2
I
(
x
,
y
)
-
I
(
x
-
step
,
y
)
-
I
(
x
+
step
,
y
)
+
2
I
x
,
y
-
I
x
,
y
-
step
-
I
x
,
y
+
step
.Source imagesA and B are divided into nonoverlapped blocks with the size of m
×
n. Refer to the ith block of source images A and B images by I
A
i and I
B
i, respectively. The focus region R
sm
l
i
A
(
x
,
y
) according to SML is defined as follows:
(9)
R
sm
l
i
A
(
x
,
y
)
=
1
SM
L
i
A
(
x
,
y
)
>
SM
L
i
B
(
x
,
y
)
0
otherwise
.However, there are always some small holes in focus region in practice. Here, the morphological filter is used to smooth blocks and remove the holes and gulfs defects, which is shown as follows:(10)
O
C
I
x
,
y
=
max
⁡
O
I
x
,
y
,
C
I
x
,
y
,
O
I
x
,
y
=
I
x
,
y
-
I
x
,
y
∘
D
,
C
I
x
,
y
=
I
x
,
y
·
D
-
I
x
,
y
.The structure elementD is 5
×
5 “diamond” matrix. O
(
I
(
x
,
y
)
) and C
(
I
(
x
,
y
)
) is the opening and closing operation. The focus region R
o
c
i
A
(
x
,
y
) according to morphological operation is defined as follows:
(11)
R
o
c
i
A
(
x
,
y
)
=
1
O
C
i
A
(
x
,
y
)
>
O
C
i
B
(
x
,
y
)
0
otherwise
.So, the final focus regionR
i
A
(
x
,
y
) is defined as
(12)
R
i
A
(
x
,
y
)
=
1
R
sm
l
i
A
(
x
,
y
)
=
1
,
R
o
c
i
A
(
x
,
y
)
=
1
0
otherwise
.The focus regionR
i
B
(
x
,
y
) is defined same as R
i
A
(
x
,
y
). The blocks of fusion image I
F
i are combined with I
A
i and I
B
i in focused regions
(13)
I
F
i
(
x
,
y
)
=
I
A
i
(
x
,
y
)
if
R
i
A
(
x
,
y
)
=
1
I
B
i
(
x
,
y
)
if
R
i
B
(
x
,
y
)
=
1
.


### 3.3. Chaotic Search

To further enhance the optimization accuracy of BBO, the modified chaotic search algorithm is adopted during each late period of iteration. Chaotic search algorithm makes full use of the ergodicity, randomness, and regularity of chaotic motion. It obtains the optimal solution by traversing every state in a small area without repetition [29]. Therefore, further chaotic search can generate several neighborhood islands based on the island with the highest HSI value to improve the performance of the proposed algorithm.The basic process is that the chaotic variables are firstly generated based on the current optimal solution. Then, the fitness value of each variable is calculated and the optimal solution of these chaotic variables is compared with the previous global optimal solution. Better solution is chosen as the current global optimal solution. This paper adopts typical logistic mapping, as shown in (14), where t
=
1,2
,
…
,
T, and  T is the maximum of iteration, which is set as 5. η is the control parameter. The system of function in (14) with η
=
4, y
1
∈
(
0,1
), and y
1
≠
{
0.25,0.5,0.75
} is a chaotic system
(14)
y
t
+
1
=
η
y
t
1
-
y
t
.
The chaotic search is as follows.Step 1.
The current optimal solutionP
b
e
s
t (the highest HSI point of current kth iteration) is mapped to domain [
0,1
] of logistic formula (15). The search area is confined in [
L
min
⁡
,
L
max
⁡
]
=
[
-
4,4
]
(15)
y
1
=
P
b
e
s
t
-
L
min
⁡
L
max
⁡
-
L
min
⁡
.Step 2.
The chaotic variabley
t is generated by logistic function in (14) and then returned by inverse mapping in term of the following equations:
(16)
p
b
e
s
t
t
=
L
min
⁡
+
L
max
⁡
-
L
min
⁡
*
y
t
.Each fitness value of each solution ofp
b
e
s
t
t
=
(
p
b
e
s
t
1
,
p
b
e
s
t
2
,
…
,
p
b
e
s
t
T
) is calculated and the fitness values of the optimal solution of them with P
b
e
s
t are compared. If the former is better than P
b
e
s
t, it replaces P
b
e
s
t as global optimal solution; otherwise, return
(17)
P
b
e
s
t
=
p
b
e
s
t
t
HSI
(
p
b
e
s
t
t
)
<
HSI
(
P
b
e
s
t
)
P
b
e
s
t
otherwise
.


### 3.4. Fusion Algorithm Description

The multifocus image fusion algorithm using biogeography-based optimization is summarized as follows.Step 1.
Define some initial parameters, such as the maximum number of speciesP
N and maximum iteration number T
max
⁡.Step 2.
Randomly generate initial islands,Island
=
{
Islan
d
1
,
Islan
d
2
,
…
,
Islan
d
P
N
}. Each island contains two parameters: Islan
d
i
=
{
m
i
,
n
i
}. m
i and n
i are width and height of ith block size, respectively. Divide source images into nonoverlapping blocks according these initial islands.Step 3.
Calculate the sharpness value of blocks to define focus regions and use these focus regions to construct a fused image. Compute the HSI functions setHS
I
i
(
x
,
y
)
=
{
HS
I
1
(
x
,
y
)
,
HS
I
2
(
x
,
y
)
,
…
,
HS
I
P
N
(
x
,
y
)
} of the fused image set I
F
i
(
x
,
y
)
=
{
I
F
1
(
x
,
y
)
,
I
F
2
(
x
,
y
)
,
…
,
I
F
P
N
(
x
,
y
)
}.Step 4.
BBO and chaotic searches are executed to update the state of islands as follows.(1)
Calculate the immigration rateλ, the emigration rate μ, and mutation rates m
s of each island.
(2)
Migration operation is as follows: the parameters{
m
i
,
n
i
} of Islan
d
i will be changed to generate new parameters based on immigration rate and emigration rate.
(3)
Mutation operation is as follows: if the mutation rate is not zero, mutate the island with the highest HSI value.
(4)
Chaotic search is as follows: generate chaotic variablesy
t based on current optimal solution and obtain new solutions of p
b
e
s
t
t
=
(
p
b
e
s
t
1
,
p
b
e
s
t
2
,
…
,
p
b
e
s
t
T
). Update the global optimal solution.Step 5.
Repeat Steps  3 and 4 until maximum iteration number is reached. Finally, the optimal block size is obtained based on the highest HSI to construct the optimal fused image.


## 3.1. HSI Function Selection

Spatial frequency [4] can measure the overall activity level in images and reflect the ability to express small details in images. The larger the value of spatial frequency in fused image, the better the performance of fusion. In the image fusion algorithms based on GA, PSO, and DE, spatial frequency is chosen as the fitness function [20–22]. It is a good fitness function to assess the fused image quality. In this paper, we also choose the spatial frequency as the HSI function of BBO
(4)
HSI
x
,
y
=
SF
x
,
y
.Spatial frequencySF
(
x
,
y
) of M
×
N size image is defined as follows:
(5)
SF
x
,
y
=
F
R
2
x
,
y
+
F
C
2
x
,
y
,
where
(6)
F
R
x
,
y
=
1
M
×
N
∑
x
=
1
M
∑
y
=
1
N
I
x
,
y
+
1
-
I
x
,
y
2
,
F
C
(
x
,
y
)
=
1
M
×
N
∑
x
=
1
M
∑
y
=
1
N
I
x
+
1
,
y
-
I
x
,
y
2
are the row and column gradients, respectively. I
(
x
,
y
) is the pixel of the image.


## 3.2. Sharpness Measurement

For multifocus image, the sharpness areas of different source images are selected to construct a fused image. How to define the sharpness areas is very important. There are many typical sharpness measurements, such as EOG, SF, and SML [18]. They measure the variation of pixels of blocks in source images to determine the sharpness. The blocks with greater values are considered as in-focus blocks. Experimental results of Huang and Jing [18] show that SML method provides better performance than EOG and SF. It differs from the usual Laplacian operator in that the absolute values of the partial second derivatives are summed instead of their actual values. The SML is defined as follows:
(7)
SML
I
x
,
y
=
∑
w
=
-
W
W
∑
h
=
-
H
H
L
I
x
+
w
,
y
+
h
2
,
where I
(
x
,
y
) is the pixel of the image. W and H is the window size. L
(
I
(
x
,
y
)
) is the modified Laplacian given as follows. Here, step is a variable that is used to set the distance between the central pixel and the pixels are used to compute the second order derivative. In this paper, a value of step
=
1 is found to produce the best results
(8)
L
I
x
,
y
=
2
I
(
x
,
y
)
-
I
(
x
-
step
,
y
)
-
I
(
x
+
step
,
y
)
+
2
I
x
,
y
-
I
x
,
y
-
step
-
I
x
,
y
+
step
.Source imagesA and B are divided into nonoverlapped blocks with the size of m
×
n. Refer to the ith block of source images A and B images by I
A
i and I
B
i, respectively. The focus region R
sm
l
i
A
(
x
,
y
) according to SML is defined as follows:
(9)
R
sm
l
i
A
(
x
,
y
)
=
1
SM
L
i
A
(
x
,
y
)
>
SM
L
i
B
(
x
,
y
)
0
otherwise
.However, there are always some small holes in focus region in practice. Here, the morphological filter is used to smooth blocks and remove the holes and gulfs defects, which is shown as follows:(10)
O
C
I
x
,
y
=
max
⁡
O
I
x
,
y
,
C
I
x
,
y
,
O
I
x
,
y
=
I
x
,
y
-
I
x
,
y
∘
D
,
C
I
x
,
y
=
I
x
,
y
·
D
-
I
x
,
y
.The structure elementD is 5
×
5 “diamond” matrix. O
(
I
(
x
,
y
)
) and C
(
I
(
x
,
y
)
) is the opening and closing operation. The focus region R
o
c
i
A
(
x
,
y
) according to morphological operation is defined as follows:
(11)
R
o
c
i
A
(
x
,
y
)
=
1
O
C
i
A
(
x
,
y
)
>
O
C
i
B
(
x
,
y
)
0
otherwise
.So, the final focus regionR
i
A
(
x
,
y
) is defined as
(12)
R
i
A
(
x
,
y
)
=
1
R
sm
l
i
A
(
x
,
y
)
=
1
,
R
o
c
i
A
(
x
,
y
)
=
1
0
otherwise
.The focus regionR
i
B
(
x
,
y
) is defined same as R
i
A
(
x
,
y
). The blocks of fusion image I
F
i are combined with I
A
i and I
B
i in focused regions
(13)
I
F
i
(
x
,
y
)
=
I
A
i
(
x
,
y
)
if
R
i
A
(
x
,
y
)
=
1
I
B
i
(
x
,
y
)
if
R
i
B
(
x
,
y
)
=
1
.


## 3.3. Chaotic Search

To further enhance the optimization accuracy of BBO, the modified chaotic search algorithm is adopted during each late period of iteration. Chaotic search algorithm makes full use of the ergodicity, randomness, and regularity of chaotic motion. It obtains the optimal solution by traversing every state in a small area without repetition [29]. Therefore, further chaotic search can generate several neighborhood islands based on the island with the highest HSI value to improve the performance of the proposed algorithm.The basic process is that the chaotic variables are firstly generated based on the current optimal solution. Then, the fitness value of each variable is calculated and the optimal solution of these chaotic variables is compared with the previous global optimal solution. Better solution is chosen as the current global optimal solution. This paper adopts typical logistic mapping, as shown in (14), where t
=
1,2
,
…
,
T, and  T is the maximum of iteration, which is set as 5. η is the control parameter. The system of function in (14) with η
=
4, y
1
∈
(
0,1
), and y
1
≠
{
0.25,0.5,0.75
} is a chaotic system
(14)
y
t
+
1
=
η
y
t
1
-
y
t
.
The chaotic search is as follows.Step 1.
The current optimal solutionP
b
e
s
t (the highest HSI point of current kth iteration) is mapped to domain [
0,1
] of logistic formula (15). The search area is confined in [
L
min
⁡
,
L
max
⁡
]
=
[
-
4,4
]
(15)
y
1
=
P
b
e
s
t
-
L
min
⁡
L
max
⁡
-
L
min
⁡
.Step 2.
The chaotic variabley
t is generated by logistic function in (14) and then returned by inverse mapping in term of the following equations:
(16)
p
b
e
s
t
t
=
L
min
⁡
+
L
max
⁡
-
L
min
⁡
*
y
t
.Each fitness value of each solution ofp
b
e
s
t
t
=
(
p
b
e
s
t
1
,
p
b
e
s
t
2
,
…
,
p
b
e
s
t
T
) is calculated and the fitness values of the optimal solution of them with P
b
e
s
t are compared. If the former is better than P
b
e
s
t, it replaces P
b
e
s
t as global optimal solution; otherwise, return
(17)
P
b
e
s
t
=
p
b
e
s
t
t
HSI
(
p
b
e
s
t
t
)
<
HSI
(
P
b
e
s
t
)
P
b
e
s
t
otherwise
.


## 3.4. Fusion Algorithm Description

The multifocus image fusion algorithm using biogeography-based optimization is summarized as follows.Step 1.
Define some initial parameters, such as the maximum number of speciesP
N and maximum iteration number T
max
⁡.Step 2.
Randomly generate initial islands,Island
=
{
Islan
d
1
,
Islan
d
2
,
…
,
Islan
d
P
N
}. Each island contains two parameters: Islan
d
i
=
{
m
i
,
n
i
}. m
i and n
i are width and height of ith block size, respectively. Divide source images into nonoverlapping blocks according these initial islands.Step 3.
Calculate the sharpness value of blocks to define focus regions and use these focus regions to construct a fused image. Compute the HSI functions setHS
I
i
(
x
,
y
)
=
{
HS
I
1
(
x
,
y
)
,
HS
I
2
(
x
,
y
)
,
…
,
HS
I
P
N
(
x
,
y
)
} of the fused image set I
F
i
(
x
,
y
)
=
{
I
F
1
(
x
,
y
)
,
I
F
2
(
x
,
y
)
,
…
,
I
F
P
N
(
x
,
y
)
}.Step 4.
BBO and chaotic searches are executed to update the state of islands as follows.(1)
Calculate the immigration rateλ, the emigration rate μ, and mutation rates m
s of each island.
(2)
Migration operation is as follows: the parameters{
m
i
,
n
i
} of Islan
d
i will be changed to generate new parameters based on immigration rate and emigration rate.
(3)
Mutation operation is as follows: if the mutation rate is not zero, mutate the island with the highest HSI value.
(4)
Chaotic search is as follows: generate chaotic variablesy
t based on current optimal solution and obtain new solutions of p
b
e
s
t
t
=
(
p
b
e
s
t
1
,
p
b
e
s
t
2
,
…
,
p
b
e
s
t
T
). Update the global optimal solution.Step 5.
Repeat Steps  3 and 4 until maximum iteration number is reached. Finally, the optimal block size is obtained based on the highest HSI to construct the optimal fused image.


## 4. Experiments

The proposed algorithm is called CSBBO. In order to verify the performance of CSBBO algorithm, such classical algorithms include the pulse coupled neural network (PCNN) [6] in spatial domain, nonsubsampled contourlet transform (NSCT) [16], and nonsubsampled shearlet transform (NSST) [17] in frequency domain, genetic algorithm (GA) [20], particle swarm optimization (PSO) [21], and differential evolution (DE) [22] based on intelligence optimization.Two kinds of experiment results, with a reference image and without a reference image, are shown in this paper. Different objective evaluation is used in these two kinds of experiments. A reference image is regarded as a ground truth image and various algorithms are assessed using the difference or structure similarity between the fused image and reference image. Without reference image, various algorithms are assessed using the information which transfers from the source images to the fused image. In the implementation, it is assumed that two multifocus source images are registered before the image fusion process. All multifocus source images and reference images are downloaded fromhttp://www.imagefusion.org/.


### 4.1. Fusion with Reference Image

The parameters of various algorithms are set as follows. For PCNN,α
L
=
0.06931, α
θ
=
0.2, and β
=
0.2. The maximal iterative number is η
=
200, V
L
=
1.0, and V
θ
=
20 [6]. For NSTC and NSST, the decomposition level is four and the directions from coarse scale to finer scale are (
4,8
,
8,16
) [15, 17]. The low-pass subband coefficients and the band-pass subband coefficients are simply merged by “averaging” scheme and the “absolute maximum choosing” scheme, respectively. For GA, CrossoverRate
=
0.01 and MutationRate
=
0.01 [20]. For PSO, c
1
=
c
2
=
2 and weight
=
0.9 [21]. For DE, the version is rand/1/bin,
CrossoverRate
=
0.6, and ScalingFactor
=
0.9 [22].With the reference image, root mean squared error (RMSE), peak signal-to-noise ratio (PSNR), and structural similarity metric (SSIM) [30] are used as quantitative assessment metrics to compare various fusion algorithms. RMSE and PSNR are commonly used to evaluate the difference between the fused image and the reference image. SSIM commonly evaluates the structure similarity between the fused image and reference image. The higher PSNR and SSIM values and the less RMSE value indicate better fusion results. They are defined as
(18)
RMSE
=
1
M
×
N
∑
i
=
1
M
∑
j
=
1
N
R
i
,
j
-
F
i
,
j
2
,
PSNR
=
10
lo
g
10
M
×
N
2
∑
i
=
1
M
∑
j
=
1
N
R
i
,
j
-
F
i
,
j
2
,
SSIM
R
,
F
=
2
μ
r
μ
f
+
c
1
μ
r
2
+
μ
f
2
+
c
1
α
·
2
σ
r
σ
f
+
c
2
σ
r
2
+
σ
f
2
+
c
2
β
·
σ
r
f
+
c
3
σ
r
σ
f
+
c
3
γ
,
where R
(
i
,
j
) and F
(
i
,
j
) are the pixel values of the reference image and the fused image, respectively. The image size is M
×
N. μ
r, μ
f, σ
r
2, σ
f
2, σ
r
f is mean, variance, and covariance, respectively. c
1, c
2, c
3 are constants closing to zero.It is impossible to guarantee that an optimal solution could be reached after a number of iterations because of the stochastic behavior of GA, PSO, DE, and BBO [20–22]. Moreover, the optimal solution almost could not be improved in late of iterations because of the premature behavior of GA, PSO, DE, and BBO. Too much iteration will result in very high computational complexity. The experimental result in [22] has shown that the performances of image fusion based on GA and DE cannot be significantly improved after 20 iterations. There is the same trend of the image fusion based on PSO, BBO, and CSBBO in our experiments. Therefore, the maximum iteration number of all intelligence optimization algorithms is set to 20 in order to reduce complexity.In order to analyze the parameters of CSBBO, population number, and maximum mutation parameter, the experiment of “Pepper” image fusion is given in Table1. Three images of “Pepper” are shown as in Figure 2. Figure 2(a) is the reference image in focus everywhere. Figures 2(b) and 2(c) focus on center and surround, respectively. All images have 512 × 512 pixels with 256 level gray scales. Considering the randomness of four intelligence optimization algorithms, experiments are repeated 30 times and the average results are saved. From Table 1, we can see that the difference between RMSE values obtained for different population number (PN) is not very high. Then, PN is set to 10 taking into account the efficiency and complexity. And the most suitable mutation parameter m
max
⁡ that gives the best RMSE value is 0.2 when PN is 10.Table 1
Average RMSE results obtained by CSBBO for Pepper images fusion.
m
max
⁡
=
0.1
m
max
⁡
=
0.2
m
max
⁡
=
0.3
m
max
⁡
=
0.4
PN = 5
1.3576
1.3440
1.3763
1.3863
PN = 10
1.3019
1.2922
1.3528
1.3577
PN = 15
1.2854
1.2718
1.3403
1.3325Pepper images: (a) reference image; (b) focus on center; (c) focus on surround.
(a)
(b)
(c)The fused results of various algorithms are shown in Figure3. The value inside the parentheses represents the optimal block size. Each algorithm can obtain a good fused image, which is almost the same as the reference image. For clearer comparison, the difference results between the fused image and the reference image are shown in Figure 4. Ideally, the difference should be zero. So, less residual features in the difference results means better performance of the fusion algorithm. It can be observed that the difference images obtained by PCNN, NSCT, and NSST have more residual features than GA, PSO, DE, BBO, and CSBBO. This indicates that these fusion algorithms based on block optimization of intelligent techniques can improve the performance of image fusion. Moreover, we can see that the difference image obtained by CSBBO has the least residue among them.Fusion results for Pepper images: (a) PCNN; (b) NSCT; (c) NSST; (d) GA; (e) PSO; (f) DE; (g) BBO; (h) CSBBO.
(a)
PCNN
(b)
NSCT
(c)
NSST
(d)
GA (13 × 8)
(e)
PSO (10 × 15)
(f)
DE (9 × 14)
(g)
BBO(8 × 9)
(h)
CSBBO (10 × 10)Difference results for Pepper images: (a)~(h) is the difference between Figure2(a) and Figures 3(a)~3(h), respectively.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Table2 shows the quantitative evaluation results of various image fusion algorithms with the reference image. As we can see, GA, PSO, DE, BBO, and CSBBO based on intelligent optimization are able to obtain higher PSNR, SSIM values and less RMSE value. Their performances are superior to PCNN, NSCT, and NSST. Among them, the performance of BBO is slightly better than GA, PSO, and DE. And CSBBO has the highest PSNR and SSIM values and the least RMSE value, which can obtain more useful features from the source images and has better performance than other fusion algorithms.Table 2
Evaluation of various image fusion algorithms with the reference image.
PCNN
NSCT
NSST
GA
PSO
DE
BBO
CSBBO
RMSE
9.9285
1.3422
1.3396
1.3013
1.2942
1.3010
1.2953
1.2922
PSNR
38.195
46.886
46.894
47.020
47.036
47.021
47.044
47.051
SSIM
0.9442
0.9986
0.9986
0.9989
0.9990
0.9989
0.9990
0.9991


### 4.2. Fusion without Reference Image

Natural digital camera images have no reference images, which contain multiple objects in focus and defocus parts because of their location at different distance from the camera. We choose two gray images, Plane (160 × 160), Clock (256 × 256), and two color images, Flower (480 × 480), Book (512 × 512). They all have two source images with different focus parts, as shown in Figure5, respectively.Source images: (a)Plane left focus; (b)Plane right focus; (c)Clock left focus; (d)Clock right focus. (e)Flower left focus; (f)Flower right focus; (g)Book left focus; (h)Book right focus.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)For these four image pairs, the fused results of various algorithms are presented in Figures6, 8, 10, and 12, respectively. To make better comparisons, the difference results of various algorithms are presented in Figures 7, 9, 11, and 13, respectively.Fusion results forPlane images: (a) PCNN; (b) NSCT; (c) NSST; (d) GA; (e) PSO; (f) DE; (g) BBO; (h) CSBBO.
(a)
PCNN
(b)
NSCT
(c)
NSST
(d)
GA (13 × 3)
(e)
PSO (9 × 6)
(f)
DE (14 × 4)
(g)
BBO (11 × 5)
(h)
CSBBO (9 × 9)Difference results forPlane images: (a)~(h) is the difference between Figure 5(b) and Figures 6(a)~6(h), respectively.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Fusion results forClock images: (a) PCNN; (b) NSCT; (c) NSST; (d) GA; (e) PSO; (f) DE; (g) BBO; (h) CSBBO.
(a)
PCNN
(b)
NSCT
(c)
NSST
(d)
GA (4 × 16)
(e)
PSO (13 × 4)
(f)
DE (16 × 6)
(g)
BBO (13 × 11)
(h)
CSBBO (10 × 8)Difference results forClock images: (a)~(h) is the difference between Figure 5(d) and Figures 8(a)~8(h), respectively.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Fusion results forFlower images: (a) PCNN; (b) NSCT; (c) NSST; (d) GA; (e) PSO; (f) DE; (g) BBO; (h) CSBBO.
(a)
PCNN
(b)
NSCT
(c)
NSST
(d)
GA (31 × 9)
(e)
PSO (7 × 32)
(f)
DE (11 × 29)
(g)
BBO (10 × 26)
(h)
CSBBO (9 × 30)Difference results forFlower images: (a)~(h) is the difference between Figure 5(f) and Figures 10(a)~10(h), respectively.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Fusion results forBook images: (a) PCNN; (b) NSCT; (c) NSST; (d) GA; (e) PSO; (f) DE; (g) BBO; (h) CSBBO.
(a)
PCNN
(b)
NSCT
(c)
NSST
(d)
GA (12 × 27)
(e)
PSO (13 × 19)
(f)
DE (8 × 17)
(g)
BBO (13 × 11)
(h)
CSBBO (15 × 11)Difference results forBook images: (a)~(h) is the difference between Figure 5(h) and Figures 12(a)~12(h), respectively.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)From Figures6 and 8 for gray imagePlane andClock, we can see that each algorithm can get a clear fused image except for PCNN. The fused image obtained by PCNN has no good visual effect because there are some blurred blocks in focus area, such as the front plane and the clock surface. Moreover, the fused image obtained by PCNN does not maintain continuous edge features, such as the book edge and the clock edge in Figure 8. The fused images obtained by NSCT, NSST, GA, PSO, DE, BBO, and CSBBO have almost no difference for visual effects.However, we can distinguish the fusion performance of these algorithms from the difference results from Figures7 and 9. It can be observed that the difference images obtained by NSCT and NSST have more residual features than GA, PSO, DE, BBO, and CSBBO. There is large residual information of front plane obtained by NSCT and NSST in Figures 7(b) and 7(c). The numbers of clock surface are obviously remained in Figures 9(b) and 9(c). From Figures 7(d)~7(g) and Figures 9(d)~9(g), it can be found that the difference results obtained by GA, PSO, DE, BBO, and CSBBO are better than other algorithms. Among them, it is obvious that CSBBO has the least residual features. The difference images obtained by GA, PSO, DE, and BBO have more residual information especially in the edge parts of the focus area, such as Figures 9(d)~9(f), and they also have more block artifacts in Figures 7(d)~7(f). Figures 6(g), 8(g), 7(g), and 9(g) show that CSBBO can detect the focus region accurately and transfer more useful features from source images to the fused image.From Figures10 and 12 for color imageFlower andBook, we can see the performance of each algorithm in vision. The fused image obtained by PCNN has distortion brightness and blurred edges, such as the switch and bricks in Figure 10. The fused images obtained by NSCT and NSST preserve a good brightness and edge features, but they still have some obvious blurred area, such as the words at the top right corner in Figure 12. The fused images obtained by GA, PSO, DE, BBO, and CSBBO improve the visual effects because these algorithms depress efficiently the blurred vision, such as the words at the top right corner in these images, which are clearer in Figure 12. Moreover, we can distinguish the fusion performance of these algorithms from the difference results from Figures 11 and 13. The difference image obtained by PCNN has the most residual features in Figures 11(a) and 13(a). And the difference images obtained by NSCT and NSST have more residual features than GA, PSO, DE, BBO, and CSBBO in Figures 11(b), 11(c), 13(b), and 13(c). Compared with GA, PSO, DE, and BBO, CSBBO has the least residual features. From the above-mentioned results, it can be seen that CSBBO can obtain clearer and more nature fused image than other algorithms.For objective assessment without the reference image, natural multifocus images should be evaluated by nonreference fusion metrics. In this paper, the spatial frequency (SF) [4], feature mutual information (FMI) [31], the edge information based metric Q
F
A
B [32], and similarity based on block metric Q
b
(
A
,
B
,
F
) [33] are used to evaluate the fusion performance of various fusion algorithms. SF evaluates the ability to express small details in the fused images, defined as formula (5). FMI calculates the amount of information conducted from the source images to the fused image, defined as
(19)
FM
I
F
A
B
=
∑
f
,
a
p
F
A
x
,
y
,
z
,
w
log
⁡
p
F
A
x
,
y
,
z
,
w
p
F
x
,
y
·
p
A
z
,
w
+
∑
f
,
b
p
F
B
x
,
y
,
z
,
w
log
⁡
p
F
B
x
,
y
,
z
,
w
p
F
x
,
y
·
p
B
z
,
w
,
where A and B are the two source images. F is the fused image. p
F
A
(
x
,
y
,
z
,
w
) and p
F
B
(
x
,
y
,
z
,
w
) are joint probability distribution. p
A
(
z
,
w
) and p
B
(
z
,
w
) are probability distribution.Q
F
A
B measures the relative amount of edge information that is transferred from the source images into the fused image. It is defined as
(20)
Q
F
A
B
=
∑
m
=
1
M
∑
n
=
1
N
Q
A
F
(
m
,
n
)
w
A
F
(
m
,
n
)
+
∑
m
=
1
M
∑
n
=
1
N
Q
A
F
(
m
,
n
)
w
A
F
(
m
,
n
)
Q
B
F
(
m
,
n
)
w
B
F
(
m
,
n
)
·
∑
m
=
1
M
∑
n
=
1
N
w
A
F
m
,
n
+
w
B
F
m
,
n
-
1
,
where Q
A
F
(
m
,
n
) and Q
B
F
(
m
,
n
) are the edge preservation values for two M
×
N source images. Q
i
F
(
m
,
n
)
=
Q
g
i
F
(
m
,
n
)
Q
o
i
F
(
m
,
n
), (
i
=
A
,
B
), Q
g
i
F
(
m
,
n
), and Q
o
i
F
(
m
,
n
) are the Sobel edge strength and orientation preservation value at location (
m
,
n
), respectively. w
A
F
(
m
,
n
) and w
B
F
(
m
,
n
) are weighted coefficients.Q
b
(
A
,
B
,
F
) is based on the universal image quality index and uses the similarity between blocks of pixels in source images and the fused image as the weighting factors for the metrics. It is defined as
(21)
Q
b
A
,
B
,
F
=
∑
w
∈
W
sim
A
,
B
,
F
∣
w
·
Q
A
,
F
∣
w
-
Q
B
,
F
∣
w
+
Q
B
,
F
∣
w
,
where w is the analysis window and W is the family of all windows. Q
(
i
,
F
∣
w
)  (
i
=
A
,
B
) is the image quality index. sim
(
A
,
B
,
F
∣
w
) is the similarity in spatial domain between the input image and the fused image.For all the nonreference fusion metrics, the larger the value is, the better the fusion performance is. Table3 shows the objective assessments of the fused images obtained by various fusion algorithms without the reference image. From this table, it can be observed that CSBBO outperforms other algorithms in terms of SF, FMI, and Q
F
A
B for all images. For Q
b
(
A
,
B
,
F
), PCNN performs a little bit better than CSBBO forPlane image, but there are obvious defects in the fused image obtained by PCNN in Figure 6(a). Based on the analysis above, it can be concluded that CSBBO has good performance in terms of quantitative and visual evaluations.Table 3
Evaluation of various image fusion algorithms without the reference image.
PCNN
NSCT
NSST
GA
PSO
DE
BBO
CSBBO
Plane
SF
10.7549
12.4720
12.5931
12.6338
12.6364
12.5951
12.6426
12.6454
FMI
0.9258
0.9273
0.9310
0.9362
0.9375
0.9370
0.9372
0.9378
Q
F
A
B
0.7195
0.7298
0.7321
0.7432
0.7524
0.7535
0.7530
0.7546
Q
b
(
A
,
B
,
F
)
0.8944
0.8645
0.8634
0.8566
0.8668
0.8769
0.8762
0.8812
Clock
SF
21.5565
22.7596
23.0138
22.9630
23.0049
23.0087
22.9725
23.0105
FMI
0.8828
0.8805
0.8829
0.8834
0.8840
0.8842
0.8842
0.8843
Q
F
A
B
0.6869
0.6830
0.6889
0.7165
0.7180
0.7186
0.7213
0.7219
Q
b
(
A
,
B
,
F
)
0.8739
0.8633
0.8624
0.8751
0.8745
0.8785
0.8795
0.8803
Flower
SF
14.7468
16.3808
16.6530
16.6400
16.6500
16.6129
16.6396
16.6539
FMI
0.8479
0.8755
0.8774
0.8808
0.8812
0.8813
0.8813
0.8815
Q
F
A
B
0.6506
0.6980
0.6994
0.7161
0.7163
0.7166
0.7172
0.7173
Q
b
(
A
,
B
,
F
)
0.8660
0.8974
0.8959
0.8972
0.8968
0.8978
0.8990
0.8995
Book
SF
27.9575
28.6606
28.9587
28.9555
28.9917
28.9727
28.9969
28.9991
FMI
0.8678
0.8631
0.8653
0.8672
0.8671
0.8670
0.8670
0.8673
Q
F
A
B
0.6961
0.6975
0.7014
0.7206
0.7211
0.7199
0.7205
0.7214
Q
b
(
A
,
B
,
F
)
0.8848
0.8546
0.8563
0.8930
0.8930
0.8939
0.8933
0.8950


## 4.1. Fusion with Reference Image

The parameters of various algorithms are set as follows. For PCNN,α
L
=
0.06931, α
θ
=
0.2, and β
=
0.2. The maximal iterative number is η
=
200, V
L
=
1.0, and V
θ
=
20 [6]. For NSTC and NSST, the decomposition level is four and the directions from coarse scale to finer scale are (
4,8
,
8,16
) [15, 17]. The low-pass subband coefficients and the band-pass subband coefficients are simply merged by “averaging” scheme and the “absolute maximum choosing” scheme, respectively. For GA, CrossoverRate
=
0.01 and MutationRate
=
0.01 [20]. For PSO, c
1
=
c
2
=
2 and weight
=
0.9 [21]. For DE, the version is rand/1/bin,
CrossoverRate
=
0.6, and ScalingFactor
=
0.9 [22].With the reference image, root mean squared error (RMSE), peak signal-to-noise ratio (PSNR), and structural similarity metric (SSIM) [30] are used as quantitative assessment metrics to compare various fusion algorithms. RMSE and PSNR are commonly used to evaluate the difference between the fused image and the reference image. SSIM commonly evaluates the structure similarity between the fused image and reference image. The higher PSNR and SSIM values and the less RMSE value indicate better fusion results. They are defined as
(18)
RMSE
=
1
M
×
N
∑
i
=
1
M
∑
j
=
1
N
R
i
,
j
-
F
i
,
j
2
,
PSNR
=
10
lo
g
10
M
×
N
2
∑
i
=
1
M
∑
j
=
1
N
R
i
,
j
-
F
i
,
j
2
,
SSIM
R
,
F
=
2
μ
r
μ
f
+
c
1
μ
r
2
+
μ
f
2
+
c
1
α
·
2
σ
r
σ
f
+
c
2
σ
r
2
+
σ
f
2
+
c
2
β
·
σ
r
f
+
c
3
σ
r
σ
f
+
c
3
γ
,
where R
(
i
,
j
) and F
(
i
,
j
) are the pixel values of the reference image and the fused image, respectively. The image size is M
×
N. μ
r, μ
f, σ
r
2, σ
f
2, σ
r
f is mean, variance, and covariance, respectively. c
1, c
2, c
3 are constants closing to zero.It is impossible to guarantee that an optimal solution could be reached after a number of iterations because of the stochastic behavior of GA, PSO, DE, and BBO [20–22]. Moreover, the optimal solution almost could not be improved in late of iterations because of the premature behavior of GA, PSO, DE, and BBO. Too much iteration will result in very high computational complexity. The experimental result in [22] has shown that the performances of image fusion based on GA and DE cannot be significantly improved after 20 iterations. There is the same trend of the image fusion based on PSO, BBO, and CSBBO in our experiments. Therefore, the maximum iteration number of all intelligence optimization algorithms is set to 20 in order to reduce complexity.In order to analyze the parameters of CSBBO, population number, and maximum mutation parameter, the experiment of “Pepper” image fusion is given in Table1. Three images of “Pepper” are shown as in Figure 2. Figure 2(a) is the reference image in focus everywhere. Figures 2(b) and 2(c) focus on center and surround, respectively. All images have 512 × 512 pixels with 256 level gray scales. Considering the randomness of four intelligence optimization algorithms, experiments are repeated 30 times and the average results are saved. From Table 1, we can see that the difference between RMSE values obtained for different population number (PN) is not very high. Then, PN is set to 10 taking into account the efficiency and complexity. And the most suitable mutation parameter m
max
⁡ that gives the best RMSE value is 0.2 when PN is 10.Table 1
Average RMSE results obtained by CSBBO for Pepper images fusion.
m
max
⁡
=
0.1
m
max
⁡
=
0.2
m
max
⁡
=
0.3
m
max
⁡
=
0.4
PN = 5
1.3576
1.3440
1.3763
1.3863
PN = 10
1.3019
1.2922
1.3528
1.3577
PN = 15
1.2854
1.2718
1.3403
1.3325Pepper images: (a) reference image; (b) focus on center; (c) focus on surround.
(a)
(b)
(c)The fused results of various algorithms are shown in Figure3. The value inside the parentheses represents the optimal block size. Each algorithm can obtain a good fused image, which is almost the same as the reference image. For clearer comparison, the difference results between the fused image and the reference image are shown in Figure 4. Ideally, the difference should be zero. So, less residual features in the difference results means better performance of the fusion algorithm. It can be observed that the difference images obtained by PCNN, NSCT, and NSST have more residual features than GA, PSO, DE, BBO, and CSBBO. This indicates that these fusion algorithms based on block optimization of intelligent techniques can improve the performance of image fusion. Moreover, we can see that the difference image obtained by CSBBO has the least residue among them.Fusion results for Pepper images: (a) PCNN; (b) NSCT; (c) NSST; (d) GA; (e) PSO; (f) DE; (g) BBO; (h) CSBBO.
(a)
PCNN
(b)
NSCT
(c)
NSST
(d)
GA (13 × 8)
(e)
PSO (10 × 15)
(f)
DE (9 × 14)
(g)
BBO(8 × 9)
(h)
CSBBO (10 × 10)Difference results for Pepper images: (a)~(h) is the difference between Figure2(a) and Figures 3(a)~3(h), respectively.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Table2 shows the quantitative evaluation results of various image fusion algorithms with the reference image. As we can see, GA, PSO, DE, BBO, and CSBBO based on intelligent optimization are able to obtain higher PSNR, SSIM values and less RMSE value. Their performances are superior to PCNN, NSCT, and NSST. Among them, the performance of BBO is slightly better than GA, PSO, and DE. And CSBBO has the highest PSNR and SSIM values and the least RMSE value, which can obtain more useful features from the source images and has better performance than other fusion algorithms.Table 2
Evaluation of various image fusion algorithms with the reference image.
PCNN
NSCT
NSST
GA
PSO
DE
BBO
CSBBO
RMSE
9.9285
1.3422
1.3396
1.3013
1.2942
1.3010
1.2953
1.2922
PSNR
38.195
46.886
46.894
47.020
47.036
47.021
47.044
47.051
SSIM
0.9442
0.9986
0.9986
0.9989
0.9990
0.9989
0.9990
0.9991


## 4.2. Fusion without Reference Image

Natural digital camera images have no reference images, which contain multiple objects in focus and defocus parts because of their location at different distance from the camera. We choose two gray images, Plane (160 × 160), Clock (256 × 256), and two color images, Flower (480 × 480), Book (512 × 512). They all have two source images with different focus parts, as shown in Figure5, respectively.Source images: (a)Plane left focus; (b)Plane right focus; (c)Clock left focus; (d)Clock right focus. (e)Flower left focus; (f)Flower right focus; (g)Book left focus; (h)Book right focus.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)For these four image pairs, the fused results of various algorithms are presented in Figures6, 8, 10, and 12, respectively. To make better comparisons, the difference results of various algorithms are presented in Figures 7, 9, 11, and 13, respectively.Fusion results forPlane images: (a) PCNN; (b) NSCT; (c) NSST; (d) GA; (e) PSO; (f) DE; (g) BBO; (h) CSBBO.
(a)
PCNN
(b)
NSCT
(c)
NSST
(d)
GA (13 × 3)
(e)
PSO (9 × 6)
(f)
DE (14 × 4)
(g)
BBO (11 × 5)
(h)
CSBBO (9 × 9)Difference results forPlane images: (a)~(h) is the difference between Figure 5(b) and Figures 6(a)~6(h), respectively.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Fusion results forClock images: (a) PCNN; (b) NSCT; (c) NSST; (d) GA; (e) PSO; (f) DE; (g) BBO; (h) CSBBO.
(a)
PCNN
(b)
NSCT
(c)
NSST
(d)
GA (4 × 16)
(e)
PSO (13 × 4)
(f)
DE (16 × 6)
(g)
BBO (13 × 11)
(h)
CSBBO (10 × 8)Difference results forClock images: (a)~(h) is the difference between Figure 5(d) and Figures 8(a)~8(h), respectively.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Fusion results forFlower images: (a) PCNN; (b) NSCT; (c) NSST; (d) GA; (e) PSO; (f) DE; (g) BBO; (h) CSBBO.
(a)
PCNN
(b)
NSCT
(c)
NSST
(d)
GA (31 × 9)
(e)
PSO (7 × 32)
(f)
DE (11 × 29)
(g)
BBO (10 × 26)
(h)
CSBBO (9 × 30)Difference results forFlower images: (a)~(h) is the difference between Figure 5(f) and Figures 10(a)~10(h), respectively.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Fusion results forBook images: (a) PCNN; (b) NSCT; (c) NSST; (d) GA; (e) PSO; (f) DE; (g) BBO; (h) CSBBO.
(a)
PCNN
(b)
NSCT
(c)
NSST
(d)
GA (12 × 27)
(e)
PSO (13 × 19)
(f)
DE (8 × 17)
(g)
BBO (13 × 11)
(h)
CSBBO (15 × 11)Difference results forBook images: (a)~(h) is the difference between Figure 5(h) and Figures 12(a)~12(h), respectively.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)From Figures6 and 8 for gray imagePlane andClock, we can see that each algorithm can get a clear fused image except for PCNN. The fused image obtained by PCNN has no good visual effect because there are some blurred blocks in focus area, such as the front plane and the clock surface. Moreover, the fused image obtained by PCNN does not maintain continuous edge features, such as the book edge and the clock edge in Figure 8. The fused images obtained by NSCT, NSST, GA, PSO, DE, BBO, and CSBBO have almost no difference for visual effects.However, we can distinguish the fusion performance of these algorithms from the difference results from Figures7 and 9. It can be observed that the difference images obtained by NSCT and NSST have more residual features than GA, PSO, DE, BBO, and CSBBO. There is large residual information of front plane obtained by NSCT and NSST in Figures 7(b) and 7(c). The numbers of clock surface are obviously remained in Figures 9(b) and 9(c). From Figures 7(d)~7(g) and Figures 9(d)~9(g), it can be found that the difference results obtained by GA, PSO, DE, BBO, and CSBBO are better than other algorithms. Among them, it is obvious that CSBBO has the least residual features. The difference images obtained by GA, PSO, DE, and BBO have more residual information especially in the edge parts of the focus area, such as Figures 9(d)~9(f), and they also have more block artifacts in Figures 7(d)~7(f). Figures 6(g), 8(g), 7(g), and 9(g) show that CSBBO can detect the focus region accurately and transfer more useful features from source images to the fused image.From Figures10 and 12 for color imageFlower andBook, we can see the performance of each algorithm in vision. The fused image obtained by PCNN has distortion brightness and blurred edges, such as the switch and bricks in Figure 10. The fused images obtained by NSCT and NSST preserve a good brightness and edge features, but they still have some obvious blurred area, such as the words at the top right corner in Figure 12. The fused images obtained by GA, PSO, DE, BBO, and CSBBO improve the visual effects because these algorithms depress efficiently the blurred vision, such as the words at the top right corner in these images, which are clearer in Figure 12. Moreover, we can distinguish the fusion performance of these algorithms from the difference results from Figures 11 and 13. The difference image obtained by PCNN has the most residual features in Figures 11(a) and 13(a). And the difference images obtained by NSCT and NSST have more residual features than GA, PSO, DE, BBO, and CSBBO in Figures 11(b), 11(c), 13(b), and 13(c). Compared with GA, PSO, DE, and BBO, CSBBO has the least residual features. From the above-mentioned results, it can be seen that CSBBO can obtain clearer and more nature fused image than other algorithms.For objective assessment without the reference image, natural multifocus images should be evaluated by nonreference fusion metrics. In this paper, the spatial frequency (SF) [4], feature mutual information (FMI) [31], the edge information based metric Q
F
A
B [32], and similarity based on block metric Q
b
(
A
,
B
,
F
) [33] are used to evaluate the fusion performance of various fusion algorithms. SF evaluates the ability to express small details in the fused images, defined as formula (5). FMI calculates the amount of information conducted from the source images to the fused image, defined as
(19)
FM
I
F
A
B
=
∑
f
,
a
p
F
A
x
,
y
,
z
,
w
log
⁡
p
F
A
x
,
y
,
z
,
w
p
F
x
,
y
·
p
A
z
,
w
+
∑
f
,
b
p
F
B
x
,
y
,
z
,
w
log
⁡
p
F
B
x
,
y
,
z
,
w
p
F
x
,
y
·
p
B
z
,
w
,
where A and B are the two source images. F is the fused image. p
F
A
(
x
,
y
,
z
,
w
) and p
F
B
(
x
,
y
,
z
,
w
) are joint probability distribution. p
A
(
z
,
w
) and p
B
(
z
,
w
) are probability distribution.Q
F
A
B measures the relative amount of edge information that is transferred from the source images into the fused image. It is defined as
(20)
Q
F
A
B
=
∑
m
=
1
M
∑
n
=
1
N
Q
A
F
(
m
,
n
)
w
A
F
(
m
,
n
)
+
∑
m
=
1
M
∑
n
=
1
N
Q
A
F
(
m
,
n
)
w
A
F
(
m
,
n
)
Q
B
F
(
m
,
n
)
w
B
F
(
m
,
n
)
·
∑
m
=
1
M
∑
n
=
1
N
w
A
F
m
,
n
+
w
B
F
m
,
n
-
1
,
where Q
A
F
(
m
,
n
) and Q
B
F
(
m
,
n
) are the edge preservation values for two M
×
N source images. Q
i
F
(
m
,
n
)
=
Q
g
i
F
(
m
,
n
)
Q
o
i
F
(
m
,
n
), (
i
=
A
,
B
), Q
g
i
F
(
m
,
n
), and Q
o
i
F
(
m
,
n
) are the Sobel edge strength and orientation preservation value at location (
m
,
n
), respectively. w
A
F
(
m
,
n
) and w
B
F
(
m
,
n
) are weighted coefficients.Q
b
(
A
,
B
,
F
) is based on the universal image quality index and uses the similarity between blocks of pixels in source images and the fused image as the weighting factors for the metrics. It is defined as
(21)
Q
b
A
,
B
,
F
=
∑
w
∈
W
sim
A
,
B
,
F
∣
w
·
Q
A
,
F
∣
w
-
Q
B
,
F
∣
w
+
Q
B
,
F
∣
w
,
where w is the analysis window and W is the family of all windows. Q
(
i
,
F
∣
w
)  (
i
=
A
,
B
) is the image quality index. sim
(
A
,
B
,
F
∣
w
) is the similarity in spatial domain between the input image and the fused image.For all the nonreference fusion metrics, the larger the value is, the better the fusion performance is. Table3 shows the objective assessments of the fused images obtained by various fusion algorithms without the reference image. From this table, it can be observed that CSBBO outperforms other algorithms in terms of SF, FMI, and Q
F
A
B for all images. For Q
b
(
A
,
B
,
F
), PCNN performs a little bit better than CSBBO forPlane image, but there are obvious defects in the fused image obtained by PCNN in Figure 6(a). Based on the analysis above, it can be concluded that CSBBO has good performance in terms of quantitative and visual evaluations.Table 3
Evaluation of various image fusion algorithms without the reference image.
PCNN
NSCT
NSST
GA
PSO
DE
BBO
CSBBO
Plane
SF
10.7549
12.4720
12.5931
12.6338
12.6364
12.5951
12.6426
12.6454
FMI
0.9258
0.9273
0.9310
0.9362
0.9375
0.9370
0.9372
0.9378
Q
F
A
B
0.7195
0.7298
0.7321
0.7432
0.7524
0.7535
0.7530
0.7546
Q
b
(
A
,
B
,
F
)
0.8944
0.8645
0.8634
0.8566
0.8668
0.8769
0.8762
0.8812
Clock
SF
21.5565
22.7596
23.0138
22.9630
23.0049
23.0087
22.9725
23.0105
FMI
0.8828
0.8805
0.8829
0.8834
0.8840
0.8842
0.8842
0.8843
Q
F
A
B
0.6869
0.6830
0.6889
0.7165
0.7180
0.7186
0.7213
0.7219
Q
b
(
A
,
B
,
F
)
0.8739
0.8633
0.8624
0.8751
0.8745
0.8785
0.8795
0.8803
Flower
SF
14.7468
16.3808
16.6530
16.6400
16.6500
16.6129
16.6396
16.6539
FMI
0.8479
0.8755
0.8774
0.8808
0.8812
0.8813
0.8813
0.8815
Q
F
A
B
0.6506
0.6980
0.6994
0.7161
0.7163
0.7166
0.7172
0.7173
Q
b
(
A
,
B
,
F
)
0.8660
0.8974
0.8959
0.8972
0.8968
0.8978
0.8990
0.8995
Book
SF
27.9575
28.6606
28.9587
28.9555
28.9917
28.9727
28.9969
28.9991
FMI
0.8678
0.8631
0.8653
0.8672
0.8671
0.8670
0.8670
0.8673
Q
F
A
B
0.6961
0.6975
0.7014
0.7206
0.7211
0.7199
0.7205
0.7214
Q
b
(
A
,
B
,
F
)
0.8848
0.8546
0.8563
0.8930
0.8930
0.8939
0.8933
0.8950


## 5. Conclusion

For spatial domain fusion, the fixed block size of source images will result in block artifacts in multifocus image fusion. To solve this problem, a novel multifocus image fusion algorithm using biogeography-based optimization is proposed in this paper. An optimized block size of each source image is obtained by BBO and chaotic search. The final fused image consisted of the optimal block size, not a fixed block size. The experiment results with or without a reference image show that CSBBO outperforms the traditional PCNN, NSCT, and NSST algorithms. Furthermore, among image fusion algorithms based on intelligence optimization, CSBBO is superior to GA, PSO, DE, and BBO in terms of visual analysis and quantitative evaluation. In the future, it is worthy to further investigate other types of fitness functions that may affect fusion performance. And it is possible to apply the proposed algorithm to other image fusions of different types of images.



---

*Source: 340675-2015-02-28.xml*