The bilateral solver for quality estimation based multi-focus image fusion

TL;DR

This paper introduces BS-QEBIF, a fast bilateral solver-based multi-focus image fusion method that produces high-quality all-in-focus images by smoothing focus-level maps while preserving edges, outperforming existing methods.

Contribution

The paper presents a novel fast bilateral solver approach for smoothing focus-level maps in multi-focus image fusion, improving quality and boundary handling.

Findings

Outperforms 13 other fusion methods objectively and subjectively.

Produces all-in-focus images with preserved details and no residual errors.

Effectively handles focus-level-changed boundary regions without artifacts.

Abstract

In this work, a fast Bilateral Solver for Quality Estimation Based multi-focus Image Fusion method (BS-QEBIF) is proposed. The all-in-focus image is generated by pixel-wise summing up the multi-focus source images with their focus-levels maps as weights. Since the visual quality of an image patch is highly correlated with its focus level, the focus-level maps are preliminarily obtained based on visual quality scores, as pre-estimations. However, the pre-estimations are not ideal. Thus the fast bilateral solver is then adopted to smooth the pre-estimations, and edges in the multi-focus source images can be preserved simultaneously. The edge-preserving smoothed results are utilized as final focus-level maps. Moreover, this work provides a confidence-map solution for the unstable fusion in the focus-level-changed boundary regions. Experiments were conducted on $25$ pairs of source images. The proposed BS-QEBIF outperforms the other $13$ fusion methods objectively and subjectively. The all-in-focus image produced by the proposed method can well maintain the details in the multi-focus source images and does not suffer from any residual errors. Experimental results show that BS-QEBIF can handle the focus-level-changed boundary regions without any blocking, ringing and blurring artifacts.

Tables

Table 1. TABLE I: Notations used in this work

IQA	Image Quality Assessment.
F	Generated all-in-focus image
$ℐ$	Multi-focus source images.
$𝒲$	Focus-level maps utilized as weights to do the fusion.
$ℳ$	Preliminary focus-level estimation.
$𝒮$	Estimated visual quality score maps of $ℐ$.
$𝒞$	Confidence map to measure the reliability of different local regions.
Conv	Convolutional neural layer.
DMOS	Difference mean opinion score: subjective visual quality scores provided in the IQA datasets.
FC	Fully connected layer.
${\text{I}}_i$	$i^{th}$ multi-focus source image in $ℐ$.
$N$	Number of multi-focus images.
$L_{\text{QNN}}$	Loss function for training QNN.
$n_{\text{QNN}}$	Mini-batch size when training QNN.
P	Image patch utilized to estimate the visual quality of its centering pixel.
Pool	Pooling layer.
QNN	Quality deep neural network to estimate visual quality.
${𝒲}^{\prime }$	Edge-preserving smoothing results using the fast bilateral solver.
$(x,y)$	Spatial coordinates.
${\theta }_j$	Parameters in QNN.
${\phi }_1$	Learning based visual quality estimation process.
${\phi }_2$	Edge preserving smoothing process.

Table 2. TABLE II: Performance of different fusion methods on 25 25 25 pairs of multi-focus images.

clock				disk
Method	$Q_G$[57]	$Q_{NMI}$ [58]	NCIE [59]	Method	$Q_G$ [57]	$Q_{NMI}$ [58]	NCIE [59]
AVE	0.6163	1.0058	0.8283	AVE	0.5018	0.8139	0.8194
LAP [20]	0.7129	1.0030	0.8301	LAP [20]	0.6769	0.8649	0.8242
FSD [21]	0.6808	0.8480	0.8235	FSD [21]	0.6511	0.7411	0.8190
GRP [22]	0.6828	0.8493	0.8236	GRP [22]	0.6577	0.7441	0.8191
RAP [23]	0.5999	0.9975	0.8295	RAP [23]	0.5974	0.8875	0.8247
NSCT-SR [8]	0.7302	1.0623	0.8327	NSCT-SR [8]	0.7061	0.9149	0.8263
NSCT [24]	0.7282	1.0135	0.8306	NSCT [24]	0.7063	0.9172	0.8264
SR [25]	0.7310	1.1098	0.8351	SR [25]	0.7159	1.0527	0.8338
GF [6]	0.7434	1.1123	0.8357	GF [6]	0.7224	1.0539	0.8338
MWGF [38] (m)	0.7396	1.1405	0.8364	MWGF [38] (m)	0.7135	1.0298	0.8327
MWGF [38] (r)	0.7430	1.1982	0.8396	MWGF [38] (r)	0.7146	1.0681	0.8346
IM [39]	0.7461	1.2119	0.8431	IM [39]	0.7208	1.1254	0.8382
CBF [46]	0.7385	1.0969	0.8344	CBF [46]	0.7104	0.9447	0.8277
QEBIF [48]	0.7460	1.1675	0.8381	QEBIF [48]	0.7190	1.0748	0.8353
Ours	0.7518	1.2407	0.8425	Ours	0.7223	1.1834	0.8423
lab				paper
Method	$Q_G$[57]	$Q_{NMI}$ [58]	NCIE [59]	Method	$Q_G$ [57]	$Q_{NMI}$ [58]	NCIE [59]
AVE	0.5484	0.8854	0.8214	AVE	0.1424	0.2882	0.8033
LAP [20]	0.7309	1.0469	0.8329	LAP [20]	0.4768	0.2995	0.8045
FSD [21]	0.7123	0.8436	0.8243	FSD [21]	0.4514	0.2934	0.8044
GRP [22]	0.7166	0.8468	0.8244	GRP [22]	0.4595	0.2932	0.8044
RAP [23]	0.6652	1.0508	0.8329	RAP [23]	0.2053	0.2970	0.8043
NSCT-SR [8]	0.7501	1.0947	0.8352	NSCT-SR [8]	0.5749	0.3109	0.8047
NSCT [24]	0.7500	1.1055	0.8355	NSCT [24]	0.5778	0.3131	0.8047
SR [25]	0.7590	1.1973	0.8405	SR [25]	0.6457	0.7813	0.8182
GF [6]	0.7639	1.1929	0.8400	GF [6]	0.6481	0.5859	0.8109
MWGF [38] (m)	0.7539	1.1913	0.8400	MWGF [38] (m)	0.6337	0.6559	0.8130
MWGF [38] (r)	0.7565	1.2263	0.8419	MWGF [38] (r)	0.6567	0.7061	0.8147
IM [39]	0.7627	1.2764	0.8448	IM [39]	0.6555	0.8130	0.8191
CBF [46]	0.7598	1.1195	0.8363	CBF [46]	0.5301	0.3524	0.8053
QEBIF [48]	0.7610	1.2264	0.8420	QEBIF [48]	0.6422	0.7113	0.8150
Ours	0.7660	1.2895	0.8457	Ours	0.6601	0.8565	0.8212
pepsi				Averaging result of 20 pairs in Lytro dataset [47]
Method	$Q_G$[57]	$Q_{NMI}$ [58]	NCIE [59]	Method	$Q_G$ [57]	$Q_{NMI}$ [58]	NCIE [59]
AVE	0.4833	0.9082	0.8220	AVE	0.3594	0.7172	0.8166
LAP [20]	0.7313	1.0126	0.8310	LAP [20]	0.7336	0.9012	0.8280
FSD [21]	0.6858	0.8417	0.8240	FSD [21]	0.6949	0.7360	0.8205
GRP [22]	0.6909	0.8436	0.8240	GRP [22]	0.6989	0.7379	0.8206
RAP [23]	0.5835	1.0228	0.8314	RAP [23]	0.5613	0.8911	0.8272
NSCT-SR [8]	0.7441	1.0774	0.8340	NSCT-SR [8]	0.7511	0.9874	0.8326
NSCT [24]	0.7449	1.0703	0.8337	NSCT [24]	0.7502	0.9592	0.8311
SR [25]	0.7557	1.1251	0.8367	SR [25]	0.7581	1.1096	0.8402
GF [6]	0.7672	1.1762	0.8392	GF [6]	0.7622	1.0993	0.8396
MWGF [38] (m)	0.7581	1.1829	0.8396	MWGF [38] (m)	0.7479	1.1020	0.8402
MWGF [38] (r)	0.7585	1.2213	0.8419	MWGF [38] (r)	0.7497	1.1264	0.8418
IM [39]	0.7712	1.2834	0.8457	IM [39]	0.7575	1.1465	0.8428
CBF [46]	0.7616	1.1147	0.8358	CBF [46]	0.7570	1.0312	0.8352
QEBIF [48]	0.7695	1.1955	0.8407	QEBIF [48]	0.7550	1.1025	0.8396
Ours	0.7748	1.2955	0.8476	Ours	0.7578	1.1762	0.8446

Table 3. TABLE III: Component analysis

Method	$Q_G$[57]	$Q_{NMI}$ [58]	NCIE [59]
Baseline	0.6992	0.9074	0.8284
Baseline+CM	0.7118	0.9650	0.8314
Baseline+Norm	0.7487	1.1713	0.8435
Ours	0.7532	1.1756	0.8436

Equations

$$\text{F}=\sum _{i=1}^N{\text{W}}_i\circ {\text{I}}_i,i\in \{1,2,\dots ,N\}.$$

$${\text{S}}_i={\phi }_1({\text{I}}_i).$$

$$\mathcal{W}={\phi }_2(\mathcal{I},\mathcal{C},\mathcal{S}).$$

$$ReLU:$$

$$Sigmoid:$$

$$L_{\text{QNN}}=\frac{1}{n_{\text{QNN}}}\sum _{l=1}^{n_{\text{QNN}}}|\text{QNN}({\text{P}}_l)-DMOS|,$$

$$△{\theta }_j=\eta \frac{\partial L_{\text{QNN}}}{\partial {\theta }_j},{\theta }_j^{new}={\theta }_j^{old}-△{\theta }_j$$

$${\text{M}}_i(x,y)=\{\begin{array}{rl}1,& \text{if min(S1(x,y),…,SN(x,y))=Si(x,y)},\\ 0,& \text{otherwise},\end{array}$$

$$\begin{array}{r}\underset{{\text{W}}_i^{\prime }}{\mathrm{argmin}}\sum _{\text{r}}{\text{C}}_i(\text{r})({\text{W}}_i^{\prime }(\text{r})-{\text{M}}_i(\text{r}){)}^2+\frac{\lambda }{2}\sum _{\text{r},\text{k}}{\hat{B}}_i(\text{r},\text{k})({\text{W}}_i^{\prime }(\text{r})-{\text{W}}_i^{\prime }(\text{k}){)}^2\end{array}$$

$${\hat{B}}_{\text{r},\text{k}}=exp(-\frac{({\text{I}}_i\text{(r)}-{\text{I}}_i\text{(k)}{)}^2}{2{\sigma }_{in}^2}-\frac{||[{\text{I}}_i^x\text{(r)},{\text{I}}_i^y\text{(r)}]-[{\text{I}}_i^x\text{(k)},{\text{I}}_i^y\text{(k)}]|{|}^2}{2{\sigma }_{xy}^2})$$

$$Q_G=\frac{{\sum }_{x=1}^{N_1}{\sum }_{y=1}^{N_2}(Q^{{\text{I}}_1\text{F}}(x,y){\tau }^{{\text{I}}_1}(x,y)+Q^{{\text{I}}_2\text{F}}(x,y){\tau }^{{\text{I}}_2}(x,y))}{{\sum }_{x=1}^{N_1}{\sum }_{y=1}^{N_2}({\tau }^{{\text{I}}_1}(x,y)+{\tau }^{{\text{I}}_2}(x,y))}$$

$$Q_{NMI}=2[\frac{MI({\text{I}}_1,\text{F})}{H({\text{I}}_1)+H(\text{F})}+\frac{MI({\text{I}}_2,\text{F})}{H({\text{I}}_2)+H(\text{F})}]$$

$$NCIE=1-H_R=1+\sum _{k=1}^K\frac{{\lambda }_k^R}{K}lo{g}_{256}\frac{{\lambda }_k^R}{K}$$

Full text

I Introduction

Various photographs are taken every day. If the light from an object point is not well converged in the focal plane, the object would be out-of-focus which leads to image blurring and loss of information. Take multi-focus source image A and B in Fig. 1 (a) as an example. These images are taken from the same scene. The out-of-focus effect leads to information loss in the left clock in A and the right clock in B. The multi-focus image fusion task [1, 2] aims at generating an all-in-focus image based on the source images, where both clocks can be in-focus in this example.

The accurate estimation of focus-level of source images is crucial in generating the all-in-focus image. In this work, we advocate utilizing image quality assessment (IQA) task to do focus-level estimation. The IQA task targets on accurately estimating image visual quality, where the increase in image blurring level makes its visual quality worse. As an example shown in Fig. 1 (b), the focus levels of the images decrease from D to F. Their visual quality also decrease from left to right. Hence, visual quality can be utilized to coarsely measure image focus levels, and generate preliminary focus-level estimation maps, pre-estimation, by comparing the visual quality scores among all source images. For example in Fig. 1 (c), the pre-estimation of image G is shown as image H.

As shown in image H, most estimation results within the object regions are correct. However, the estimation results in the focus-level-changed boundary regions and background out-of-focus regions are not accurate. Thus pre-estimations are not appropriate to be directly used as focus-level maps. In this work, the state-of-the-art edge-preserving smoothing filter, the fast bilateral solver [3], is first introduced to the image fusion task. With the fast bilateral solver, the focus-level map I is generated by smoothing H and preserving edges in G simultaneously. The focus-level map, image I, is used to generate the all-in-focus image.

To better estimate the focus-levels, we further analyze the source images. These source images are supposed to be focused at different regions, and one region cannot be focused in all multi-focus images. However, it is possible for some regions, especially the background, to be out of focus in all multi-focus source images. For example in Fig. 1 (a), image C illustrates the absolute visual-quality difference between source image A and B, where the brighter color corresponds to the larger difference. In the yellow square patches, the difference between A and B is big representing their big difference in visual quality scores and focus levels. The focus-level estimations in these regions are more likely to be accurate, and the square regions can be regarded as reliable regions. While in the rectangular and circular patches, the differences are small. The rectangular patches in both source images are out-of-focus, and the circular patches locate at the focus-level-changed boundary regions which cover both in-focus and out-of-focus regions. It is more challenging to measure the slightly focus-level difference for the rectangular and circular patches, and they should be regarded as less reliable regions. In this work, a confidence map is proposed to measure the reliability of different regions. With the confidence map, different regions can be treated differently and adaptively to generate more accurate focus-level maps.

To sum up, a novel multi-focus image fusion method, the fast Bilateral Solver for Quality Estimation Based Image Fusion (BS-QEBIF), is proposed. We follow the focus-level weighted summation pipeline where the all-in-focus image is generated by pixel-wise summing up all multi-focus images with their focus-levels as weights. The contributions are mainly three-fold. First, the visual quality is first introduced to image fusion task to help estimate image focus levels. With rich images with subjective quality scores provided in IQA datasets, various supervised training methods can be explored in the fusion procedure. Second, the fast bilateral solver [3] is first adopted to smooth the pre-estimation results and generate the focus-level maps as weights. Instead of directly applying the fast bilateral solver to do the fusion, we explore different settings and design the confidence map accordingly. Benefiting from the bilateral solver, good and robust fusion results can be obtained. Finally, the confidence map is proposed to help improve focus-level estimation. The higher confidence scores correspond to the more reliable regions. The proposed method employs this concept in the focus-level estimation and effectively improves the fusion performance, especially in the focus-level-changed boundary regions.

II Related work

Various multi-focus image fusion methods were proposed in the last decades. Generally speaking, the image fusion methods can be divided into three categories, defocus-modeling methods, transform-based methods, and spatial-frequency methods. The defocus-modelling methods [4, 5] defocus the input images by a designed filter to reserve the blur effect. The transform-based methods [6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25] employed various transformations to do the fusion. The* spatial-frequency* methods concentrate on accurately measuring the focus or sharpness levels of the multi-focus images [26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37]. For example, the simplest method is averaging (AVG) where the multi-focus source images are averaged to generate the all-in-focus image. Besides, MWGF [38] and IM [39] presented multi-scale structure-based and morphological filtering-based focus estimations respectively. The proposed method also concentrates on exploring the focus level estimations and belongs to spatial-frequency methods.

In the last decades, neural network and learning based methods were widely employed in image processing tasks and achieved impressive success. Some image fusion methods also try to utilize machine learning to help focus level estimation. For example, some methods [40, 41, 42, 30] adopt pulse coupled neural networks (PCNN) to do image fusion where PCNN does not require the training procedure. Some other methods [43, 44, 45] train neural networks or dictionaries where focus level estimation was regarded as a classification problem to divide a local region into an in-focus or out-of-focus one. To do the training, [44, 45] collect some all-in-focus images and manually blur them using the Gaussian filters, while [40] manually labels in-focus/out-of-focus regions from source images as training sets.

However, all above methods have not made fully use of machine learning and trained a deep neural network. The main reason is the lack of large amount of training data with labels. In this work, we first introduce a highly-related task, Image Quality Assessment (IQA), to help effectively estimate focus levels. With the large amount of training data with labels provided in the IQA dataset, various deep neural networks could be explored.

Besides, different edge-preserving smoothing filters were explored in earlier works to make full use of spatial consistency during the fusion process. CBF [46] adopts the cross bilateral filter (CBF) to compute the weights to measure the strength of details in multi-focus source images. Li et al. [6], Nejati et al. [47] and Guan et al. [48] adopted the guided filter to do the fusion. In this work, a recent proposed state-of-the-art edge-preserving smoothing filter, the fast bilateral solver, is first adopted in image fusion task. Moreover, the confidence map is proposed to help the focus-level estimation with the fast bilateral solver in BS-QEBIF.

III The proposed BS-QEBIF method

The pipeline of the proposed BS-QEBIF method is summarized in Fig. 2. The all-in-focus image F is generated by pixel-wise summarizing the Hadamard products [49] of multi-focus source images \cal{I}$$=\{\textbf{I}_{1},\dots,\textbf{I}_{N}\} with their corresponding focus-level maps \cal{W}$$=\{\textbf{W}_{1},\dots,\textbf{W}_{N}\} as weights, i.e.

[TABLE]

$N$ is the number of multi-focus source images. All notations used in this work are summarized in TABLE I. The focus-level weights $\cal{W}$ are mainly obtained via steps $\varphi_{1}$ and $\varphi_{2}$.

In $\varphi_{1}$, the visual quality score for each pixel is estimated to generate the visual quality score maps \cal{S}$$=\{\textbf{S}_{1},\dots,\textbf{S}_{N}\}.

[TABLE]

Next, the focus-level maps of all images $\cal{W}$ are estimated based on visual quality score maps $\cal{S}$ and the source images $\cal{I}$. During this process, the confidence maps $\cal{C}$ are employed to improve the estimation, i.e.

[TABLE]

In the following sections, $\varphi_{1}$ (Sec. III-A) and $\varphi_{2}$ (Sec. III-B) would be introduced in details.

III-A Estimation of visual quality score maps $\varphi_{1}$

As discussed earlier, the focus level of a region can be approximately estimated by its visual quality score. Thus a Quality deep Neural Network (QNN) is proposed to estimate visual quality. In QNN, the large amount of training data with labels are from IQA datasets.

Architecture of QNN is illustrated in Fig. 3. To avoid the blocking artifacts in quality score map $\textbf{S}_{i}$, every pixel in $\textbf{I}_{i}$ is evaluated by QNN. As the input of QNN, each pixel in $\textbf{I}_{i}$ is first normalized by removing the local mean and dividing by the standard deviation. The local mean and standard deviation are computed within a neighboring area of size $7\times 7$ which is popularly utilized in state-of-the-art IQA models [50, 51, 52, 53]. After normalization, an image patch $\textbf{P}_{i}(x,y)$ of size $32\times 32$ centering at $(x,y)$ is utilized as the input of QNN to estimate the visual quality of $\textbf{I}_{i}(x,y)$, i.e. $\textbf{S}_{i}(x,y)=\texttt{QNN}(\textbf{P}_{i}(x,y))$. The centering position of a $32\times 32$ patch is set as $(17,17)$ in this work. Details of the three layer types used in QNN, i.e. convolutional neural layer (Conv), pooling layer (Pool) and fully connected layer (FC), are introduced below.

The convolutional neural layer (Conv) convolves the image patch P with the $50$ learned filters and obtain $50$ response maps. The filters are of size $7\times 7$. The stride size is $1$. 2. 2.

The pooling layer (Pool) is a way of sub-sampling and has been widely used in deep learning. In this work, both max-pooling and min-pooling are adopted simultaneously for each response map. Both the kernel size and stride size are $26$. 3. 3.

The fully connected layer (FC) takes the flattened results of previous layer as input and connects it with all neurons in this FC layer. The non-linear activation function is applied to each filter response. In QNN, Rectified Linear Unit (ReLU [54]) and Sigmoid are utilized for $\mathrm{FC_{1}}(100)$ and $\mathrm{FC_{2}}(1)$ respectively. The activation function of ReLU and Sigmoid are illustrated as follows.

[TABLE]

where $a$ and $b$ represent the input and output of the activation function respectively.

Training procedure of QNN: The parameters in QNN are learned through back-propagation. The key to success for such supervised learning task is the use of large amount of labeled training images. Since the principle of out-of-focus blurring is similar to the Gaussian blurring effect, images blurred by Gaussian filter in the IQA dataset, the LIVE dataset [55], are utilized as training samples. The subjective visual quality scores provided in the dataset [55], difference mean opinion scores (DMOS), were utilized as training labels, where the smaller DMOS represents better visual quality. Since these Gaussian images were generated with a circular-symmetric 2D Gaussian kernel of a certain standard deviation sigma $\sigma$, all pixels in one image can be regarded as of the same focus levels and share the same label. In the implementation, these blurry images are non-overlappingly divided into image patches P of size $32\times 32$ to serve as the training input of QNN, while the corresponding DMOS of the whole image are utilized as the training labels. The loss function for training QNN is

[TABLE]

where $n_{\texttt{QNN}}$ is the mini-batch size [56] and $l$ is the index of P. For any parameter $\theta_{j}$ in the QNN, the updated process is illustrated as

[TABLE]

where $\eta$ represents the learning rate, $\frac{\partial L_{\texttt{QNN}}}{\partial\theta_{j}}$ is the corresponding derivative. It should be noted that the current QNN model was trained on gray images. It is easy to extend the model to process color images.

III-B Estimation of focus-level maps $\varphi_{2}$

The estimation of focus-level maps can be divided into two steps. First, the preliminary focus-level estimation results, pre-estimation $\cal{M}$, can be estimated by regarding the image with best visual quality as the in-focused one. However, the pre-estimations contain much noise and are not appropriate to be directly used as the weights in the fusion. Therefore, the fast bilateral solver [3] is employed in the second step to edge-preserved smooth $\cal{M}$ and generate focus-level maps $\cal{W}$. During the second stage, the confidence map $\cal{C}$ and the normalization are employed.

III-B1 Pre-estimation $\cal{M}$

The pre-estimation of focus levels $\cal{M}$ are estimated based on the visual quality score maps $\cal{S}$. The in-focus image is estimated to be the one with best visual quality. i.e.

[TABLE]

where the smaller value in $\textbf{S}_{i}$ indicates better visual quality. Thus $\textbf{M}_{i}(x,y)=1$ indicates that $\textbf{I}_{i}$ is approximated as the in-focused source image at position $(x,y)$.

An example of $\cal{M}$ and the fusion result directly using $\cal{M}$ are shown in Fig. 4. The fusion result shown in image (b) is not satisfactory since some unwanted sudden changes exist. In the object region, $\cal{M}$ is expected to be smooth, such as the right clock in image (a) since they are of the same focus level. While in some focus-level-changed boundary region, the lines in $\cal{I}$ are expected to be well preserved in $\cal{M}$, such as the lines between the two clocks. These unwanted sudden changes in $\cal{M}$ result in the unwanted sudden changes in the fusion results as shown in (b).

Reason of inaccuracy : These unwanted sudden changes in $\cal{M}$ are mainly result from the instability of QNN when the focus-level changes of inputs are small. Fig. 5 illustrates the change of estimated visual quality scores using QNN (y-axis) with different focus levels (x-axis). In this experiment, $100$ image patches distorted by Gaussian blur are utilized for evaluation, where the averaging visual quality estimation scores of these $100$ patches are reported in y-axis. It can be seen that the overall trend is monotonous. However, the estimation results of small focus-level changing images are not stable as shown in the enlarged detail block in the lower right corner. Therefore, when using QNN to evaluate several regions with similar quality scores, the results may not be stable.

III-B2 Edge-preserving smoothing $\cal{M}$ to generate $\cal{W}$

In order to smooth within object regions and preserve boundary edges, an edge-preserving smoothing filter, the fast bilateral solver [3], is utilized to generate the final estimation of focus-level maps $\cal{W}$.

Besides, a lot of information is lost when calculating $\cal{M}$ in $\varphi_{1}$. Because for any same-sign quality difference, $\cal{M}$ would always be given the same decision. To make full use of information in quality score maps, the confidence map is proposed to pixel-wise measure the reliability of $\cal{M}$, and help improve focus level estimation accuracy. In addition, the normalization is employed. In this section, the confidence map, fast bilateral solver and normalization are introduced.

Confidence map: As shown in Fig. 5, the quality estimation results tend to be more accurate in regions with larger focus-level difference than those with smaller difference. Based on such observation, $\cal{C}$ are generated based on the quality difference as illustrated in the upper part of Fig. 6. For the orange patches, the difference in their quality maps is big. There is a big chance that the estimation result in this region is accurate. On the contrary, things become uncertain when the difference is small, which mainly results from two reasons. First, both patches may be out of focus such as the gray ones and result in unreliable results in $\cal{M}$. Second, the patches may locate at the focus-level-changed boundary region such as the blue ones. To be specific, the estimation results in the boundary patches are not reliable as illustrated in the ‘Boundary regions’ block of Fig. 6. In (a), the yellow and red point is out-of-focus and in-focus respectively. However, the patches utilized for estimating the quality scores of these two points cover both the in-focus and out-of-focus regions as shown in (b). Since the two patches are of similar focus levels, the estimated visual quality scores for these two points are similar as shown in (c), which are not reliable. Therefore, the corresponding points in $\cal{M}$ should be regarded as unreliable regions.

Considering the above conditions, the confidence map divides the whole image as reliable and unreliable regions according to their visual-quality-score difference. The computation of $\cal{C}$ is summarized in Algorithm 1. Thr = $0.1$ in the implementation. The confidence maps of all source images are the same.

The fast bilateral solver: In the fast bilateral solver [3], the pre-estimations $\cal{M}$ serves as the filter input which need to be smoothed. The multi-focus images $\cal{I}$ are utilized as the reference images that guide which lines to be preserved. Meanwhile, the confidence maps $\cal{C}$ are adopted to measure the reliability of filter input $\cal{M}$.

In order to preserve the lines in $\cal{I}$ and smooth $\cal{M}$, an optimization problem is constructed [3] which includes two terms, a data-fidelity term and a reference image based smoothness term, i.e.

[TABLE]

where $\textbf{W}^{\prime}_{i}$ is the $i^{th}$ edge-preserving smoothed result. r and k are utilized to represent pixel $(x_{r},y_{r})$ and $(x_{k},y_{k})$ respectively in this subsection for simplify. $\textbf{W}^{\prime}_{i}(\textbf{r})$ is a scalar representing the value at pixel r in map $\textbf{W}^{\prime}_{i}$.

The data-fidelity term $\textbf{C}_{i}(\textbf{r})(\textbf{W}^{\prime}_{i}(\textbf{r})-\textbf{M}_{i}(\textbf{r}))^{2}$ help minimize the squared residual between filter input $\cal{M}$ and output $\cal{W^{\prime}}$ while weighted by confidence map $\cal{C}$. A bigger weight is given to a reliable region of the input $\cal{M}$ compared with that given to an unreliable region.

The reference image based smoothness term $\frac{\lambda}{2}\sum_{\textbf{r},\textbf{k}}{}$ $\hat{B}_{i}({\textbf{r},\textbf{k}})(\textbf{W}^{\prime}_{i}(\textbf{r})-\textbf{W}^{\prime}_{i}(\textbf{k}))^{2}$ tried to minimize the difference between any two pixels, r and k, in $\textbf{W}^{\prime}_{i}$ weighted by $\hat{B}_{i}({\textbf{r},\textbf{k}})$. Specifically, $\hat{B}_{i}({\textbf{r},\textbf{k}})$ is a bistochastized version of a bilateral affinity matrix which measures the similarity between r and k in $\textbf{I}_{i}$ in terms of spatial position and intensity, i.e.

[TABLE]

where $\textbf{I}_{i}(\textbf{r})$ and $[\textbf{I}_{i}^{x}{\textbf{(r)}},\textbf{I}_{i}^{y}\textbf{(r)}]$ represents intensity and the spatial position of pixel r in $\textbf{I}_{i}$. $\sigma_{xy}$ and $\sigma_{in}$ are the parameters control the extent of the spatial and intensity support of the filter respectively.

Specifically, the weight $\hat{B}_{i}({\textbf{r},\textbf{k}})$ is big if r and k have similar intensities or are spatially close in $\textbf{I}_{i}$. Thus the big weights given to the similar intensity help minimize the difference in r and k in $\textbf{W}_{i}$, and results in edge-preserving in lines and smoothness in plane regions. Besides, the closer pixels may have higher correlation compared with not related ones in the minimization. To speed up, the fast bilateral solver [3] treats $\hat{B}$ as a “splat / blur / slice” procedure so that the optimization process can be solved in the bilateral-space and be reduced into a simple least-square problem. It is proved that the fast bilateral solver is fast, robust and has been successfully generalized to many new domains, such as stereo, depth super resolution, colorization and semantic segmentation [3].

The fusion results are highly correlated with different settings of the bilateral solver. First, directly applying the fast bilateral solver may lead to obvious block artifact as shown in the first row of Fig. 7 (a). Specifically, the block size is equal to $\sigma_{xy}$. This is because $\sigma_{xy}$ helps measure the spatial similarity of the filter inputs. In the implementation, the input is non-overlappingly divided into several blocks for measuring the spatial difference. Pixels belonging to the same block and two neighboring block would have a big difference in $\hat{B}$ and result in the block artifact of the filter output.

In order to eliminate the blocking artifact, the block motion process is proposed as illustrated in the ‘Block motion process’ of Fig. 7. In this example, the input image, image (a), is of size $6\times 6$. In order to measure the spatial similarity in (10), the filter input is non-overlappingly divided into blocks of size $\sigma_{xy}\times\sigma_{xy}=2\times 2$ beginning from the left-up corner. Thus the result shown in (a) suffers from the block artifact from left-up corner. To eliminate the blocking artifact using block motion process, the fast bilateral solver is processed again and get another filter output as (b). This bilateral solver process is repeated for $\sigma_{xy}=2$ times according to the block sizes. During each process, the beginning point of dividing blocks shifted one-pixel in each direction from the red point in (a) to the yellow point in (b). The final fusion result is obtained by averaging over all the filter outputs as shown in Fig. 7 (c). The filter output of ‘clock’ with block motion process is shown in (b) of Fig. 7. Benefiting from the block motion process, the block artifacts have been effectively removed.

Besides, the choice of $\sigma_{xy}$ has influence on the filter outputs. As shown in Fig. 7, $\sigma_{xy}$ corresponds to the respective filed. In the image fusion task, the estimation results may be not accurate in some boundary regions if $\sigma_{xy}$ is too small. Besides, the speed of the bilateral solver is highly depended on $\sigma_{xy}$. Thus $\sigma_{xy}$ cannot be too big. To increase the operation speed and maintain the performance, $\sigma_{xy}=8$ in the implementation. Besides, $\lambda=64$.

Normalization: Via the fast bilateral solver [3], the edge-preserving smoothed result $\cal{W^{\prime}}$ are calculated. Next, the normalization is proceed to get the final focus-level maps $\cal{W}$ as weights.

At each position $(x,y)$, all values in $N$ $\textbf{W}_{i}^{\prime}(x,y),i\in\{1,2,\dots,N\}$ are normalized to sum as 1 first. Next, a non-linear transformation, sigmoid, with mean = 0.5, $\sigma$ = 40, is employed. In this way, most values in $\cal{W}$ would concentrated around the maximum value $1$ and minimum value [math], which would help improve the fusion results.

IV Experiment

Several experiments were conducted to evaluate the performance of the proposed method on $25$ pairs of source images subjectively and objectively. First, experimental configurations are introduced in Sec. IV-A, including multi-focus source images and objective evaluation metrics. Moreover, component analysis results are reported in Sec. IV-C.

IV-A Experimental configurations

IV-A1 Multi-focus source images

$25$ pairs of source images were utilized for the evaluation. The first $5$ pairs are commonly used images pairs, ‘clock’, ‘pepsi’, ‘lab’, ‘paper’ and ‘disk’, as shown in Fig. 8. The remaining ones are from the Lytro multi-focus dataset [47]. These image pairs vary in multiple aspects, which provides a good representation of various situations encountered in practice.

IV-A2 Objective evaluation metrics

In this work, three widely used evaluation metrics, $Q_{G}$ [57], $Q_{NMI}$ [58] and NCIE [59], were adopted. Default parameters provided by respective metrics were used and bigger score indicates better performance. These metrics evaluate the fusion results in multiple respects and provide a solid evaluation.

$Q_{G}$ [57] measures the success of edge information transferred from the multi-focus source images to the fusion results. Take \cal{I}$$=\{\textbf{I}_{1},\textbf{I}_{2}\} as an example, $Q_{G}$ can be defined as,

[TABLE]

where $N_{1}$ and $N_{2}$ represents the number of coordinates in spatial domain. $Q^{\textbf{I}_{1}\textbf{F}}$ is the sum of $Q^{\textbf{I}_{1}\textbf{F}}_{g}$ and $Q^{\textbf{I}_{1}\textbf{F}}_{o}$ of $\textbf{I}_{1}$, where $Q^{\textbf{I}_{1}\textbf{F}}_{g}$ and $Q^{\textbf{I}_{1}\textbf{F}}_{o}$ represent element-wise product of edge strength and orientation preservation value respectively. $Q^{\textbf{I}_{2}\textbf{F}}$ is defined similarly for $\textbf{I}_{2}$. $\tau^{\textbf{I}_{1}}$ and $\tau^{\textbf{I}_{2}}$ represents weights for $\textbf{I}_{1}$ and $\textbf{I}_{2}$ respectively. 2. 2.

Normalized mutual information $Q_{NMI}$ [58] measures the success of mutual information transferred from the multi-focus source images to the fusion results based on the information theory. The normalized mutual information is more stable than traditional mutual information metric [60]. $Q_{NMI}$ can be defined as,

[TABLE]

where $H(\textbf{I}_{1}),H(\textbf{I}_{2})$ and $H(\textbf{F})$ represent the marginal entropy of source image $\textbf{I}_{1},\textbf{I}_{2}$ and F respectively. $MI(\textbf{I}_{1},\textbf{F})=H(\textbf{I}_{1})+H(\textbf{F})-H(\textbf{I}_{1},\textbf{F})$ measures the mutual information between $\textbf{I}_{1}$ and F. $H(\textbf{I}_{1},\textbf{F})$ is the joint entropy between $\textbf{I}_{1}$ and F. $MI(\textbf{I}_{2},\textbf{F})$ can be calculated similarly. 3. 3.

Nonlinear Correlation Information Entropy NCIE [59] measures the nonlinear correlation of source images $\cal{I}$ and the fused image F.

[TABLE]

where $H_{R}$ represents the nonlinear joint entropy. $K=N+1$ is the number of concerned variables including the source images $\cal{I}$ and the fusion result F. $R$ is the nonlinear correlation matrix of the concerned variables, including $\cal{I}$ and F. $\lambda_{k}^{R}$ is the $k^{th}$ eigenvalue of $R$.

IV-A3 Image fusion methods for comparison

13 traditional and state-of-the-art multi-focus image fusion methods were utilized for comparison, including AVE, LAP [20], FSD [21], GRP [22], RAP [23], NSCT-SR [8], NSCT [24], SR [25], GF [6], MWGF [38], IM [39], CBF [46] and QEBIF [48]. Code of implementing these methods were obtained online or directly from the authors, and their default parameters are used to make a fair comparison. Especially, MWGF [38] provides two sets of parameters for two cases. In the implementation, we provide the results under both assumptions and adopt (m) and (r) for mis-registered and well registered cases respectively. Since the proposed method targets on fusing gray images, all color images were converted into gray ones before the fusion process in all methods.

IV-B Fusion Results

The proposed method was evaluated on all $25$ pairs of source images, objectively and subjectively.

IV-B1 Objective evaluation

TABLE II summarized $Q_{G}$ [57], $Q_{NMI}$ [58] and NCIE [59] of the proposed BS-QEBIF method and the $13$ other image fusion methods. The results on each source image pair are grouped together with the name of the source images ahead. Results of the first five pairs were summarized individually. The averaging results of all $20$ pairs in Lytro dataset [47] were summarized in the last block in the table. The best performing methods are shown in bold.

The proposed method achieves the best performance in almost all situations. The consistent satisfactory performance on $25$ pairs of source images prove the stability of the proposed method. Moreover, the good performance of the proposed method in terms of $Q_{G}$ [57] and $Q_{NMI}$ [58] indicate its strong ability to maintain sharpness information and mutual information from the source images. Besides, the fusion result of the proposed method is highly correlated with the source images which can be concluded from NCIE [59]. Nevertheless, the proposed method BS-QEBIF with the fast bilateral solver works better than QEBIF. It shows that the fast bilateral solver is more suitable than the guided filter in the proposed pipeline.

IV-B2 Subjective evaluation

To further evaluate the performance of the proposed method, subjective evaluation were conducted by visualizing the fusion results as shown in Fig 8, 9 and 10. Fig. 8 illustrates the fusion result of the proposed method on the first 5 pairs of source images. The fusion results can well preserve the details in the source images.

Fig. 8 and Fig. 9 analyze the results in two major challenges in image fusion task, mis-registration and fusion at focus-level-changed boundary regions. The fusion results of the proposed method are compared with those of 13 other fusion methods. Specifically, the source images and the generated fusion results are shown in (a) and (b-o) respectively. For MWGF [38] in TABLE II, the result of (r) performs better than (m). Thus MWGF [38] (r) is adopted for visualization in image (k). For each method, the fused image F (first line) and its normalized difference with $I_{1}$, $\textbf{D}_{norm}$ (second line), are shown. $\textbf{D}_{norm}$ is calculated by $\textbf{D}_{norm}=\frac{\textbf{D}-\textbf{D}_{min}}{\textbf{D}_{max}-\textbf{D}_{min}}$, where $\textbf{D}=\textbf{F}-\textbf{I}_{1}$ is the difference map. $\textbf{D}_{max}$ and $\textbf{D}_{min}$ denote the maximum and minimum values among all $14$ difference maps, which ensures the fairness of the comparison.

In Fig 9, multi-focus source images ‘lab’ suffer from the mis-registration, especially on the man’s head as shown in (a). Since $\textbf{I}_{1}$ focuses on the man’s head, thus this area of $\textbf{D}_{norm}$ should be zero in the ideal case. It can be seen that AVE, LAP [20], FSD [21], GRP [22], RAP [23], NSCT-SR [8], NSCT [24] and CBF [46] failed to handle the mis-registered well, because big residual errors can be observed in the blocks in their normalized difference maps. SR [25], GF [6] and QEBIF [48] works better, but few residual errors still exist in the white block. The proposed method, MWGF [38] and IM [39] can handle thus mis-registered problems well since no residual errors exist in the block. This example indicates the proposed method can handle source images with mis-registration well and fusion result does not suffer from the blocking, ringing or blurring artifacts, which may occur in other methods.

Similar results can be observed in Fig. 10. The left block concentrates on measuring the performance of mis-registration. In this case, only the proposed method and MWGF [38] work well and do not suffer from any residual errors. The right blocks concentrate on measuring the accuracy in the focus-level-changed boundary region. The left and right parts within the block are focused on $\textbf{I}_{1}$ and $\textbf{I}_{2}$ respectively, Therefore, the left part in $\textbf{D}_{norm}$ should be zero in the ideal case. However, only the proposed method and IM [39] work well and accurately fuse in the focus-level-changed boundary regions. For all other methods, clear horizontal lines can be observed in the left part of $\textbf{D}_{norm}$. Considering these two blocks, the proposed method works the best on ‘pepsi’.

From these examples, it can be seen that the proposed method is effective and works well subjectively. The fusion result of the proposed method does not suffer from any residual errors with source images mis-registered. Besides, it can deal well with the focus-level-changed boundary regions, avoid blocking, ringing and blurring artifacts.

IV-C Component analysis

In this section, the effectiveness of the proposed confidence map (CM) and normalization (Norm) are evaluated objectively and subjectively. Objectively, they were evaluated on all 25 pairs of source image in term of $Q_{G}$ [57], $Q_{NMI}$ [58] and NCIE [59]. The averaging results on these 25 pairs were summarized in TABLE III, where ‘baseline’ corresponds to the fusion results using $\cal{W}^{\prime}$ with values in the confidence map were set to identical; ‘baseline+CM’ corresponds to the fusion results using $\cal{W}^{\prime}$ with the proposed confidence map; ‘baseline+Norm’ corresponds to the fusion result while adding normalization on the ‘baseline’. ‘Ours’ represents the proposed method which utilizes the confidence map and normalization.

The proposed confidence map (CM) and normalization (Norm) are proved to be effective by comparing the results in ‘Baseline+CM’ and ‘Baseline+Norm’ to ‘Baseline’ respectively. Performance gains are achieved in both cases. ‘Ours’ combines them together and achieved the best results indicating ‘CM’ and ‘Norm’ cooperate well with each other. Besides, ‘CM’ and ‘Norm’ are effective to transfer mutual information from source images to the fusion results concluded from $29.6\%$ performance gain in $Q_{NMI}$ [58]. Subjectively, fusion results of ‘clock’ under different settings are shown in Fig. 11. The fusion results of (a) ‘Baseline’, (b) ‘Baseline+CM’, (c) ‘Baseline+Norm’ and (d) ‘Ours’ are shown. Without normalization, image (a) and (b) suffers from the residual errors on the focus-level-changed boundary region as shown in the upper block. Benefiting from the normalization process, image (c) is much smoother and appears better. However, (c) still cannot handle focus-level-changed boundary regions accuratly as we can see on the bottom block. The confidence map effectively improve the fusion results as we can see from the bottom block in (d), where the fusion errors have be effectively eliminated. To sum up, CM and Norm are both effective. Especially, the proposed confidence map achieved the desired goal (discussed in ‘Boundary regions’ of Fig. 6) to improve focus-level-changed boundary regions even it does not achieved an apparent performance improvement during objective evaluation.

V Conclusions

In this work, a quality estimation based multi-focus image fusion method with the fast bilateral solver (BS-QEBIF) is proposed. The visual quality is estimated to help measure the focus levels since the visual quality of an image is highly correlated with its focus level. In addition, the confidence map is proposed to measure the reliability of different local regions. The fusion results of BS-QEBIF can well maintain the details in the multi-focus images and do not suffer from the ringing or blocking artifacts.