# Robust Multisensor Image Matching Using Bayesian Estimated Mutual Information

Authors: Lurong Shen; Xinsheng Huang; Yuzhuang Yan; Yongbin Zheng; Wanying Xu
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

Mutual information (MI) has been widely used in multisensor image matching, but it may lead to mismatch among images with messy background. However, additional prior information can be of great help in improving the matching performance. In this paper, a robust Bayesian estimated mutual information, named as BMI, for multisensor image matching is proposed. This method has been implemented by utilizing the gradient prior information, in which the prior is estimated by the kernel density estimate (KDE) method, and the likelihood is modeled according to the distance of orientations. To further improve the robustness, we restrict the matching within the regions where the corresponding pixels of template image are salient enough. Experiments on several groups of multisensor images show that the proposed method outperforms the standard MI in robustness and accuracy and is similar with Pluim’s method. However, our computation is far more cost saving.


---

## Body

## 1. Introduction

The matching of multisensor images is an important preprocessing step for many applications, such as navigation, data-fusion, and visualization tasks. Due to different physical characteristics of different sensors, the relationship between the intensities of the matching pixels is often complex. Visual features present in one sensor image may not appear in other ones, and vice versa; image contrast may differ from each other in the same regions; multiple intensities in one image may be mapped to a single intensity in other images, and vise versa. The uncorrelated intensities make the normal intensity-based matching methods ineffective; therefore, the matching of multisensor image becomes challenging in computer vision.Previous works on multisensor image matching can be generally categorized into two classes [1]. The first class is to find some invariant representations among multisensor images, such as point feature [2], edge feature [3], contour feature [4] and edge orientation maps [5], then the problem can be solved by utilizing the traditional monosensor matching methods. Nevertheless, the feature extraction has been regarded as one of the most difficult problems in image processing, especially in multisensor cases, as the feature present in one image may not appear in other ones. Additionally, these methods usually omit most of the image details, which may lead to a wrong matching.The second class is to use an invariant similarity measure of images. These methods are conducted directly in the intensity domain, and they have been successfully used in a variety of multisensor registrations. Within these methods, the images are treated as combinations of different sets. Each set stands for a separate gray level, assuming that the same component of the same object is consistent in one image, and the pixels of the same component with intensities nearx in one image correspond to the pixels with intensities near y in the other image. In other words, there may be no pixel-to-pixel correspondence in intensity, but there do exist set-to-set correspondences between two images; namely, ignoring the intensity property, the topology of the corresponding sets must be similar between two images. The corresponding relationship can be well shown by the joint intensity scatter plot (JISP) or the joint histogram (JH) [6], and by measuring their dispersion one can quantify the matching degree. One of the most successful measures is the entropy. The lower the entropy is, the more tightly clustered the scatter plot will be and, hence, the closer the matching will be. Mutual information (MI) [7, 8] or normalized mutual information (NMI) [9] has been considered as the state-of-the-art matching frameworks of the entropy-based methods.Nevertheless, the robustness of the MI-based matching methods is questionable. A mismatched result of the NMI between an infrared (IR) and a visible (VS) image is shown in Figure1, in which statistical relationship between the images is decreased by the messy background. In addition, a “Grey stripe” registration also demonstrates the pitfall of NMI in [10].A mismatching example of MI. (a) An infrared template image, (b) a visible reference image and the matching result, and (c) the matching result of our method. There are many background textures around the target (car) of the infrared image, but only a few in the VS image, which decrease the statistical relationship and lead to a mismatching when using the NMI method. The matching result is shown by superimposing the Sobel edges of the template image directly on the reference image.
(a)
(b)Recently, many approaches have been proposed to improve the MI by incorporating additional spatial information. Studholme et al. [11] introduce a prior of labeled connecting regions as an additional information channel to the MI. Guo and Lu [12] adopt the gradient vector flow (GVF) to enhance the original intensity images before using the MI, which improves the success rate. Kim et al. [13] propose a new method based on a 3D joint histogram (JH), which integrates the intensity information and the edge orientation information. The method provides more robustness. Lee et al. [14] extend [13] by using a modified GVF to extend the capture range of gradient-based MI. Wang and Tian [15] encode the images into 32 bins according to the gradient magnitude and orientation and then compute the MI directly on the encoded field, which improves the computation efficiency of the MI. Liu et al. [16] improved this method with an adaptive combination strategy of an intensity-based MI and a coded gradient-based MI. Fan et al. [17] propose a method that combines the spatial information through a feature-based selection mechanism, such as the Harris corner, and the JH is increased to 4D. Loeckx et al. [18] propose conditional mutual information (cMI) as a new similarity for image registration; cMI starts from a 3D joint histogram incorporating, besides the intensity dimensions, also a spatial dimension expressing the location of the joint intensity pair. Pluim et al. [19] modify the MI-based method by multiplying the MI with a term of the gradient. The gradient term not only seeks to align locations of high gradient magnitude but also aims for a similar orientation of the gradients at these locations. This method is referred to as GNMI in this paper.These improvements for the MI-based methods strongly depend on the invariant characteristics mentioned in the first class, for example, the gradient information of magnitude and orientation, as well as the corner features. Actually, since they take full advantages of the two classes, we classify them as the third class and regard this additional information as prior information in this paper.It should be pointed out that some methods which introduce the prior into MI are problematic. For example, in [15], the pixels are encoded according to their gradient magnitude and orientation. The intensities in the standard MI are substituted with gradient codings. It is known that the gradient orientation appears to be parallel to each other in multisensor images, and they must be distributed along the diagonal of JH in the true matching. However, as the codings function similarly as the intensities of the standard MI in JH, and MI measures the dispersion of JH only, this method still gives a high correlation even though the compact bins are located far away from the diagonal of JH. Moreover, the gradient magnitudes participating in the coding are not so reliable, as the magnitudes vary a lot in multisensor images [13, 14]. Some better methods are presented in [13, 16, 19], in which they treat the measure of gradient vectors between multisensor images as a weight of the MI. However, they suffer from expensive computation.In this paper, a new similarity measure that incorporates the prior information into MI by the Bayesian method is proposed. For simplicity, it is named as Bayesian MI (BMI) here. BMI is deduced by substituting classical probabilities in MI with the Bayesian probability, such as the prior, the likelihood, and the posterior. For example, the classical conditional probability can be regarded as a likelihood probability or a posterior probability in the Bayesian framework and through which we can introduce some subjective experiences in the estimation. To avoid the statistic of JH, the BMI is defined directly on each pixel pair, and the images can then be measured by the mean BMI of the valid matching pixel pairs.Then, an implementation of the BMI is proposed based on the prior information of gradient. Since the magnitude information between multisensor images is not so reliable, we only concern whether their magnitudes are salient or not and treat the salient ones as the saliency. Generally, saliency of image contributes more in matching [20, 21], so we restrict the pixel pairs of matching within a mask where the pixels of the template image are salient.In the implementation, the orientation of each pixel is defined firstly. The saliency pixels are defined as the absolute value of their orientations, while the nonsaliency ones are defined as zero. Secondly, the prior distribution is estimated based on the kernel density estimate (KDE) method, and the likelihood between orientations is modeled by the Gaussian function according to their distance. The closer the orientation is, the more similar the pixel pair will be. Finally, the sequential statistic BMI based on the gradient information is modeled as well.At the end of this paper, experiments are conducted on six groups of images which are captured by different sensors, and the results suggest that the proposed method is more robust and more accurate than the standard NMI and is similar to Pluim’s method [19]. Moreover, our computation is far more cost saving.The rest of this paper is outlined as follows. Section2 briefly introduces the standard MI, and the proposed BMI is deduced by Bayesian formula. Based on the gradient prior information, the application of the BMI is shown in Section 3. In Section 4, some experimental results and discussions are presented. Finally, conclusions are drawn in Section 5.


## 2. Proposed Method

### 2.1. Mutual Information

The MI method is of general significance because it makes few assumptions on image intensities. It assumes neither linear correlation nor even functional correlation but only statistical dependence [10].Given two images,A and B, one can define their joint probability (JP) pAB(a,b) by simply normalizing the JH. Let pA(a) and pB(a) denote the corresponding probabilities, which can be obtained directly from pAB(a,b). The MI between A and B based on Shannon entropy is given by
(1)MI(A,B)=H(A)+H(B)-H(A,B)=∑a,bpAB(a,b)logpAB(a,b)pA(a)pB(b).Apart from (1), MI can also be described in many other forms, as summarized by [22]:
(2)MI(A,B)=H(A)-H(A∣B)=∑a,bpAB(a,b)logpA∣B(a∣b)pA(a),
where H(A) and H(B) are the Shannon entropy of images A and B, which indicate the amount of uncertainty about A and B, respectively. H(A∣B) and H(B∣A) denote the conditional entropy, indicating the uncertainty left in A or B when B or A is known. H(A,B) is the joint entropy, which indicates the total uncertainty of A and B. When b is known, the conditional probability of a is noted as  pAB(a,b); it can be easily deduced from their JP by the classical probability theory.In fact,-log(a) can be regarded as the self-information of a, and the entropy H(a) represents its mean. Similarly, the mutual information between a and b can be defined as
(3)MI(a,b)=logpAB(a,b)pA(a)pB(b),
and MI(A,B) can be deemed as the mathematical expectation of MI(a,b).


### 2.2. Bayesian Mutual Information.

BMI is deduced directly from the MI by regarding the conditional probabilitypA∣B(a∣b) in (2) as a posterior probability, which can be estimated according to the prior information based on the Bayesian method.LetpA(a) be the prior probability; then, according to the Bayesian formula, we have
(4)p^A∣B(a∣b)=p^B∣A(b∣a)pA(a)p^B(b)=p^B∣A(b∣a)pA(a)∑ap^B∣A(b∣a)pA(a),
where p^B|A(b∣a) and p^B(b) denote the likelihood probability and evidence probability, respectively. p^B(b) can be estimated through the joint of prior and likelihood by the total probability formula. Here, p^ is used for distinguishing the Bayesian probabilities from the classical probability p. Using (2) and (4), the BMI can be defined as
(5)BMI(A,B)=∑a,bpAB(a,b)logp^B∣A(b∣a)∑ap^B∣A(b∣a)pA(a).Similar to the MI, (5) needs the statistic of JP beforehand. For convenience, we refine (5) into a sequential statistic form. As pAB(a,b)=(1/N)∑x,yδ(ax,y-a)·δ(bx,y-b), where N is the total pixel number in the matching, δ(x) is the Kronecker function, ax,y and bx,y are the pixel values of A and B at (x,y), respectively. Thus, (5) can be rewritten as
(6)BMI(A,B)=1N∑x,yBMI(ax,y,bx,y)=1N∑x,ylogp^B∣A(bx,y∣ax,y)p^B(bx,y)=1N∑x,ylogp^B∣A(bx,y∣ax,y)∑ax,yp^B∣A(bx,y∣ax,y)pA(ax,y),
where BMI(ax,y,bx,y) denotes the BMI of the pixel pair (ax,y,bx,y). Obviously, a higher BMI BMI(ax,y,bx,y) can be achieved from higher likelihood and lower evidence. Equation (6) is applicable regardless of whether a and b are discrete or continuous, while (5) is for the discrete form only.BMI can be comprehended favourably using the communication theory when regardingax,y, p^B∣A(bx,y∣ax,y) and bx,y as the message source, the channel, and the sink; respectively, BMI is just the evaluation of the correct rate after transmission [23].It should be noticed thatpA(ax,y), p^B∣A(bx,y∣ax,y), p^B(b), and their logarithmic values can be calculated offline; moreover, as log(xy)=log(x)+log(y), the cost of our computation can be remarkably saved by utilizing the lookup-table (LUT) technique.


## 2.1. Mutual Information

The MI method is of general significance because it makes few assumptions on image intensities. It assumes neither linear correlation nor even functional correlation but only statistical dependence [10].Given two images,A and B, one can define their joint probability (JP) pAB(a,b) by simply normalizing the JH. Let pA(a) and pB(a) denote the corresponding probabilities, which can be obtained directly from pAB(a,b). The MI between A and B based on Shannon entropy is given by
(1)MI(A,B)=H(A)+H(B)-H(A,B)=∑a,bpAB(a,b)logpAB(a,b)pA(a)pB(b).Apart from (1), MI can also be described in many other forms, as summarized by [22]:
(2)MI(A,B)=H(A)-H(A∣B)=∑a,bpAB(a,b)logpA∣B(a∣b)pA(a),
where H(A) and H(B) are the Shannon entropy of images A and B, which indicate the amount of uncertainty about A and B, respectively. H(A∣B) and H(B∣A) denote the conditional entropy, indicating the uncertainty left in A or B when B or A is known. H(A,B) is the joint entropy, which indicates the total uncertainty of A and B. When b is known, the conditional probability of a is noted as  pAB(a,b); it can be easily deduced from their JP by the classical probability theory.In fact,-log(a) can be regarded as the self-information of a, and the entropy H(a) represents its mean. Similarly, the mutual information between a and b can be defined as
(3)MI(a,b)=logpAB(a,b)pA(a)pB(b),
and MI(A,B) can be deemed as the mathematical expectation of MI(a,b).


## 2.2. Bayesian Mutual Information.

BMI is deduced directly from the MI by regarding the conditional probabilitypA∣B(a∣b) in (2) as a posterior probability, which can be estimated according to the prior information based on the Bayesian method.LetpA(a) be the prior probability; then, according to the Bayesian formula, we have
(4)p^A∣B(a∣b)=p^B∣A(b∣a)pA(a)p^B(b)=p^B∣A(b∣a)pA(a)∑ap^B∣A(b∣a)pA(a),
where p^B|A(b∣a) and p^B(b) denote the likelihood probability and evidence probability, respectively. p^B(b) can be estimated through the joint of prior and likelihood by the total probability formula. Here, p^ is used for distinguishing the Bayesian probabilities from the classical probability p. Using (2) and (4), the BMI can be defined as
(5)BMI(A,B)=∑a,bpAB(a,b)logp^B∣A(b∣a)∑ap^B∣A(b∣a)pA(a).Similar to the MI, (5) needs the statistic of JP beforehand. For convenience, we refine (5) into a sequential statistic form. As pAB(a,b)=(1/N)∑x,yδ(ax,y-a)·δ(bx,y-b), where N is the total pixel number in the matching, δ(x) is the Kronecker function, ax,y and bx,y are the pixel values of A and B at (x,y), respectively. Thus, (5) can be rewritten as
(6)BMI(A,B)=1N∑x,yBMI(ax,y,bx,y)=1N∑x,ylogp^B∣A(bx,y∣ax,y)p^B(bx,y)=1N∑x,ylogp^B∣A(bx,y∣ax,y)∑ax,yp^B∣A(bx,y∣ax,y)pA(ax,y),
where BMI(ax,y,bx,y) denotes the BMI of the pixel pair (ax,y,bx,y). Obviously, a higher BMI BMI(ax,y,bx,y) can be achieved from higher likelihood and lower evidence. Equation (6) is applicable regardless of whether a and b are discrete or continuous, while (5) is for the discrete form only.BMI can be comprehended favourably using the communication theory when regardingax,y, p^B∣A(bx,y∣ax,y) and bx,y as the message source, the channel, and the sink; respectively, BMI is just the evaluation of the correct rate after transmission [23].It should be noticed thatpA(ax,y), p^B∣A(bx,y∣ax,y), p^B(b), and their logarithmic values can be calculated offline; moreover, as log(xy)=log(x)+log(y), the cost of our computation can be remarkably saved by utilizing the lookup-table (LUT) technique.


## 3.  Implementation

In this section, we are about to implement the proposed BMI based on the prior information of the gradient. Actually, BMI is universal as long as the prior and the likelihood can be modeled, such as the intensity and the local statistic between monosensor matching, the gradient orientation and magnitude, and the edges but not the intensity of multisensor images.As mentioned in the introduction, the saliency includes most of the information [24]; thus, it would contribute more in matching. In [21], Qin et al. extract the salient vector for each pixel using the principal axis analysis (PAA). The more similar the corresponding vectors are, the higher weight the pixel pairs will be in JP. In [20], Kim et al. extract the salient regions from expanding the edges; then the MI is calculated on those regions in which the images are both salient. This method enhances the statistical correlation between the IR and VS images.In this study, we simply define the saliency as those pixels whose gradient magnitude are higher than a threshold, which is correlated to the mean of the full image. Different from Kim’s scheme, the proposed BMI is computed only on those pixel pairs in which the template image is salient, regardless of the reference image.LetR and T denote the reference image and the template image, respectively, and Ruv denotes the subimage cropped from the top-left corner (u,v) of R, which has the same size as T. Both of the two images are filtered by a 5×5 Gaussian window to suppress the noise.The realization of the BMI includes the gradient orientation mapping, the prior probability estimation, and the likelihood probability modelingIt should be noted that, in the following implementation, it is assumed that there is only translation between two images, so are the experiments. This assumption does not affect manifesting the effect of our method.


### 3.1. Gradient Orientation Mapping.

Letθx,yT and θx,yRuv be the gradient orientation maps of T and Ruv, and (x,y) is the coordinate of pixel. Taking θx,yT for an example, it is defined as
(7)θx,yT={angle(Vxy),if(∇Ty>0∩Gxy>t),p,if(∇Ty=0∩Gxy>t),0,ifGxy<t,Vxy=sgn(∇Ty)[∇Tx∇Ty],Gxy=|Vxy|,
where angle(x) is a function which returns the phase of vector x, Gxy is its magnitude, and sgn(·) is the sign function, which reverses the vectors in the 3rd and 4th quadrant and limits the θx,yT to [0,π]. τ is the threshold which separates the gradient of salient structures from noises, and
(8)τ=aN∑x,yGxy,
where a is an empirical percentage of the mean magnitude, N is the total pixel number of I. Practically, a∈[0.8,1.2]. Here, the fact that gradient orientations are parallel to each other between multisensor images is utilized. In other words, their orientations are either the same or opposite. When the gradient magnitude is small, θx,yT is set to zero, and the pixels with horizontal orientation are defined as π.LetΔθ be the resolution angle; then we can get the discrete form of (7)
(9)cx,yT=round(θx,yT+0.5Δθ).In this paper, the pixels ofθx,yT=0 are named as zero-code pixels, and those of θx,yT≠0 are named as nonzero-code pixels, and those nonzero-code pixels are defined as saliency. The saliency of T forms the mask in which the matching is conducted, and we denote the set of the corresponding coordinations (x,y) as V.


### 3.2. Prior Probability Density Estimation

Normally, the prior probability can be defined by the histogram ofCx,yT. However, when the size of the template image is small (small sample size problem), the histogram will be incapable to character the probability precisely by the classical probability theory, which requires as many samples as possible. For this problem, kernel density estimate (KDE), or Parzen window method [25, 26], is a good solution. KDE has been widely used for nonparametric distribution estimation [27]. Therefore, the probability density of  θx,yT>0 can be written as
(10)fT(θ)=1Ch1∑(x,y)∈VK1(∥θ-θx,yTh1∥2),
where K1(·) is the kernel function and h1 corresponds to the bandwidth. The larger h1 is, the smoother fT(θ) will be. Constant Ch1 makes ∫fT(θ)dθ=1.Actually, the KDE method distributes eachθx,yT to its neighbors by a bounded K1(·) with h1 and makes fT(θ) more accurate and more robust than the discrete one. Frequently used kernel functions include Uniform, Triangle, Epanechnikov, and Normal [28]. In this study, the Normal Gaussian kernel function is employed as kernel function.


### 3.3. Likelihood Probability Density Model

Likelihood is a subjective probability in our work. It is a connection between the prior and the evidence. In this section, we sort the valid pixel pairs under the mask into two classes: nonzero-code to nonzero-code matching and nonzero-code to zero-code matching, as shown in Figure2. Their likelihoods will be modeled, respectively.Figure 2
Valid pixel pairs in matching are sorted into two classes: nonzero-code to nonzero-code and nonzero-code to zero-code, which are marked with “1” and “2,” respectively. The saliency (nonzero-code) is marked by diagonal lines.As for the first class, its likelihood can be modeled according to their distance betweenθx,yT and θx,yRuv. It is clear that the closer the orientation is, the higher the likelihood will be. Taking two random orientations, θr and θt, as an example, their likelihood probability density can be modeled as
(11)f^Ruv∣T(θr∣θt)=1Ch2K2(∠(θr,θt)h2),θr,θt>0,
where K2(·) is a hat-like function which is high in the middle but low on both sides, for example, the Gaussian function. h2 defines its bandwidth. ∠(θr,θt) is the intersection angle between the two orientations, and 0≤∠(θr,θt)≤π/2, which is calculated by
(12)∠(θr,θt)=min{|θr-θt-π|,|θr-θt|,|θr-θt+π|}.Ch2 is a constant which ensures f^Ruv|T(θr|θt) satisfy the boundary condition of ∫f^Ruv∣T(θr∣θt)dθr=1.Figure3 shows an example of f^Ruv∣T(θr∣θt) when h2=12/π. It indicates that the bigger the intersection angle is, the lower the likelihood will be, and it approaches to 0 as ∠(θr,θt) tend to π/2. Combining (10) and (11), the evidence probability density can be obtained:
(13)f^Ruv(θ)=∫fT(θt)·f^Ruv∣T(θ∣θt)dθt.Figure 3
The likelihood probability densityf^Ruv∣T(θr∣θt) versus the intersection angle ∠(θr,θt) based on the Gaussian function h2=12/π.The continuous function,f^Ruv(θ), is just a smoothed form of fT(θ) by using the function K2, and they will be over-lapped when h2≤π/180, as shown in Figure 4.Figure 4
Curves offT(θ) and f^Ruv(θ) under different h2. Apparently, f^Ruv(θ) is a smoothed result of fT(θ). The bigger h2 is, the smoother f^Ruv(θ) will be.With regard to the nonzero-code to zero-code matching, it can be treated as mismatching, and its likelihood can be defined as 0. Finally, we have(14)f^Ruv∣T(θr∣θt)={1Ch2K2(∠(θr,θt)h2),if(θr,θt>0),0,if(θr=0,θt>0).


### 3.4. Practical Issues

Now, the BMI ofRuv and T can be derived by substituting ax,y, bx,y, A,  and B in (6) with θx,yT, θx,yRuv, T, and  Ruv, respectively. First, some issues should be noticed.(1)
BMI will be defined as 0 if the likelihood equals 0, such as the second class.(2)
We can see that there do exist some pixel pairs that  satisfyf^Ruv∣T(θx,yRuv∣θx,yT)<f^Ruv(θx,yRuv), which results in a negative BMI. Sometimes, the mean BMI may be negative, too, especially in those mismatched locations. But it never occurs in the standard MI [22]. In this study, the pixel pairs with negative BMI are regarded as mismatched, and their BMIs are limited to zero, so that they are consistent with the second class.Finally, the BMI-based matching formula betweenRuv and T can be derived as
(15)BMI(T,Ruv)=1N×∑(x,y)∈V{logf^Ruv∣T(θx,yRuv∣θx,yT)f^Ruv(θx,yRuv),if(f^Ruv∣Tf^Ruv>1,θx,yRuv≠0),0,otherwise.The true matching of the template imageT can be obtained by maximizing the BMI(T,Ruv) through out the interval (u,v) of the reference image R.


## 3.1. Gradient Orientation Mapping.

Letθx,yT and θx,yRuv be the gradient orientation maps of T and Ruv, and (x,y) is the coordinate of pixel. Taking θx,yT for an example, it is defined as
(7)θx,yT={angle(Vxy),if(∇Ty>0∩Gxy>t),p,if(∇Ty=0∩Gxy>t),0,ifGxy<t,Vxy=sgn(∇Ty)[∇Tx∇Ty],Gxy=|Vxy|,
where angle(x) is a function which returns the phase of vector x, Gxy is its magnitude, and sgn(·) is the sign function, which reverses the vectors in the 3rd and 4th quadrant and limits the θx,yT to [0,π]. τ is the threshold which separates the gradient of salient structures from noises, and
(8)τ=aN∑x,yGxy,
where a is an empirical percentage of the mean magnitude, N is the total pixel number of I. Practically, a∈[0.8,1.2]. Here, the fact that gradient orientations are parallel to each other between multisensor images is utilized. In other words, their orientations are either the same or opposite. When the gradient magnitude is small, θx,yT is set to zero, and the pixels with horizontal orientation are defined as π.LetΔθ be the resolution angle; then we can get the discrete form of (7)
(9)cx,yT=round(θx,yT+0.5Δθ).In this paper, the pixels ofθx,yT=0 are named as zero-code pixels, and those of θx,yT≠0 are named as nonzero-code pixels, and those nonzero-code pixels are defined as saliency. The saliency of T forms the mask in which the matching is conducted, and we denote the set of the corresponding coordinations (x,y) as V.


## 3.2. Prior Probability Density Estimation

Normally, the prior probability can be defined by the histogram ofCx,yT. However, when the size of the template image is small (small sample size problem), the histogram will be incapable to character the probability precisely by the classical probability theory, which requires as many samples as possible. For this problem, kernel density estimate (KDE), or Parzen window method [25, 26], is a good solution. KDE has been widely used for nonparametric distribution estimation [27]. Therefore, the probability density of  θx,yT>0 can be written as
(10)fT(θ)=1Ch1∑(x,y)∈VK1(∥θ-θx,yTh1∥2),
where K1(·) is the kernel function and h1 corresponds to the bandwidth. The larger h1 is, the smoother fT(θ) will be. Constant Ch1 makes ∫fT(θ)dθ=1.Actually, the KDE method distributes eachθx,yT to its neighbors by a bounded K1(·) with h1 and makes fT(θ) more accurate and more robust than the discrete one. Frequently used kernel functions include Uniform, Triangle, Epanechnikov, and Normal [28]. In this study, the Normal Gaussian kernel function is employed as kernel function.


## 3.3. Likelihood Probability Density Model

Likelihood is a subjective probability in our work. It is a connection between the prior and the evidence. In this section, we sort the valid pixel pairs under the mask into two classes: nonzero-code to nonzero-code matching and nonzero-code to zero-code matching, as shown in Figure2. Their likelihoods will be modeled, respectively.Figure 2
Valid pixel pairs in matching are sorted into two classes: nonzero-code to nonzero-code and nonzero-code to zero-code, which are marked with “1” and “2,” respectively. The saliency (nonzero-code) is marked by diagonal lines.As for the first class, its likelihood can be modeled according to their distance betweenθx,yT and θx,yRuv. It is clear that the closer the orientation is, the higher the likelihood will be. Taking two random orientations, θr and θt, as an example, their likelihood probability density can be modeled as
(11)f^Ruv∣T(θr∣θt)=1Ch2K2(∠(θr,θt)h2),θr,θt>0,
where K2(·) is a hat-like function which is high in the middle but low on both sides, for example, the Gaussian function. h2 defines its bandwidth. ∠(θr,θt) is the intersection angle between the two orientations, and 0≤∠(θr,θt)≤π/2, which is calculated by
(12)∠(θr,θt)=min{|θr-θt-π|,|θr-θt|,|θr-θt+π|}.Ch2 is a constant which ensures f^Ruv|T(θr|θt) satisfy the boundary condition of ∫f^Ruv∣T(θr∣θt)dθr=1.Figure3 shows an example of f^Ruv∣T(θr∣θt) when h2=12/π. It indicates that the bigger the intersection angle is, the lower the likelihood will be, and it approaches to 0 as ∠(θr,θt) tend to π/2. Combining (10) and (11), the evidence probability density can be obtained:
(13)f^Ruv(θ)=∫fT(θt)·f^Ruv∣T(θ∣θt)dθt.Figure 3
The likelihood probability densityf^Ruv∣T(θr∣θt) versus the intersection angle ∠(θr,θt) based on the Gaussian function h2=12/π.The continuous function,f^Ruv(θ), is just a smoothed form of fT(θ) by using the function K2, and they will be over-lapped when h2≤π/180, as shown in Figure 4.Figure 4
Curves offT(θ) and f^Ruv(θ) under different h2. Apparently, f^Ruv(θ) is a smoothed result of fT(θ). The bigger h2 is, the smoother f^Ruv(θ) will be.With regard to the nonzero-code to zero-code matching, it can be treated as mismatching, and its likelihood can be defined as 0. Finally, we have(14)f^Ruv∣T(θr∣θt)={1Ch2K2(∠(θr,θt)h2),if(θr,θt>0),0,if(θr=0,θt>0).


## 3.4. Practical Issues

Now, the BMI ofRuv and T can be derived by substituting ax,y, bx,y, A,  and B in (6) with θx,yT, θx,yRuv, T, and  Ruv, respectively. First, some issues should be noticed.(1)
BMI will be defined as 0 if the likelihood equals 0, such as the second class.(2)
We can see that there do exist some pixel pairs that  satisfyf^Ruv∣T(θx,yRuv∣θx,yT)<f^Ruv(θx,yRuv), which results in a negative BMI. Sometimes, the mean BMI may be negative, too, especially in those mismatched locations. But it never occurs in the standard MI [22]. In this study, the pixel pairs with negative BMI are regarded as mismatched, and their BMIs are limited to zero, so that they are consistent with the second class.Finally, the BMI-based matching formula betweenRuv and T can be derived as
(15)BMI(T,Ruv)=1N×∑(x,y)∈V{logf^Ruv∣T(θx,yRuv∣θx,yT)f^Ruv(θx,yRuv),if(f^Ruv∣Tf^Ruv>1,θx,yRuv≠0),0,otherwise.The true matching of the template imageT can be obtained by maximizing the BMI(T,Ruv) through out the interval (u,v) of the reference image R.


## 4. Experimental Validation

To evaluate the performance of the proposed method, firstly, ten pairs of IR and VS images are used in the test as shown in Figures5 and 6, respectively. The IR images are captured by a long-wave IR camera, whose sizes are 108×168, and they are taken as the template image; the VS images are captured by a CCD camera, whose sizes are 256×256, and they are taken as the reference images. All the images are captured on a sunny afternoon, with only some random translation transform between them. As the ground truth is unknown for the sets of real images, it is estimated by manually matching in Photoshop. We also test these four methods on two pairs of image scropped from [29], as shown in Figures 7 and 8.Figure 5
Ten groups of IR reference images.Figure 6
Ten groups of VS reference images.(Group 11) A pair of downtown images captured at night and in daytime.
(a)
(b)(Group 12) A pair of images of Mauna Loa.
(a)
(b)The matching results obtained from the proposed method are compared with the NMI, GNMI [19], and GCMI methods [15]. These methods are chosen because they have been widely used in multisensor image matching, and they provide relatively good matching performance. In the NMI and GNMI methods, the densities are estimated by JH with a fixed number of 64 * 64 bins for both images. The intensity values are linearly mapped into these bins. The parameters used in these methods are  h1=10π/180/3, h2=90π/180/3, and a=0.8, the matching results of our method are shown in Figure 10.The comparison results of the NMI, GNMI, GCMI, and our method are shown in Figure9. A global searching strategy was applied to get a full view of the correlation curve. The figures indicate the following.(1)
The BMI outperforms the NMI significantly. It can be seen that (1) the difference between the highest and the 2nd higher peak, which measures the interference suppressing ability, is more salient in BMI than in NMI through all the tests. In particular, the NMI mismatches in the 8th, 9th, and 11st test, while the BMI is still applicable. (2) The highest peak of BMI, which measures the position accuracy, is much sharper than that of NMI. The matching results are shown in Table 1, from which we can see that all the methods have matching errors, but our method is more accurate.(2)
Because GNMI, GCMI, and BMI adopt gradient information as measure, the performances of our method are similar to the GNMI and GCMI; they all have similar matching accuracy and the salient highest peak. However, as discussed previously, the computation of our method is remarkably cost saving when the LUT technique is used. In the MATLAB environment, the cost ratio of the four methods is about 1 : 2 : 3 : 4 (BMI : GCMI : NMI : GNMI), and the proportion can be further improved when using other program tools, such as the C and C++ language.Table 1
Comparison of the position accuracy.
Group
Ground truth
NMI/error
GNMI/error
GCMI/error
BMI/error
1
(72, 44)
(72, 37)/(0, 7)
(73, 44)/(1, 1)
(73, 45)/(1, 1)
(73, 44)/(1, 0)
2
(66, 16)
(65, 15)/(1, 1)
(65, 15)/(1, 1)
(65, 15)/(1, 1)
(65, 15)/(1, 1)
3
(38, 57)
(39, 57)/(1, 0)
(39, 57)/(1, 0)
(39, 57)/(1, 0)
(39, 57)/(1, 0)
4
(104, 43)
(105, 43)/(1, 0)
(105, 41)/(1, 2)
(105, 42)/(1, 1)
(105, 42)/(1, 1)
5
(49, 45)
(49, 44)/(0, 1)
(49, 44)/(0, 1)
(49, 44)/(0, 1)
(49, 44)/(0, 1)
6
(99, 14)
(103, 13)/(4, 1)
(99, 13)/(0, 1)
(99, 13)/(0, 1)
(99, 14)/(0, 0)
7
(69, 28)
(70, 26)/(1, 2)
(70, 26)/(1, 2)
(70, 25)/(1, 3)
(70, 26)/(1, 2)
8
(107, 55)
(142, 54)/(35, 1)
(107, 55)/(0, 0)
(107, 55)/(0, 0)
(107, 55)/(0, 0)
9
(118, 38)
(143, 78)/(35, 40)
(119, 38)/(1, 0)
(118, 37)/(0, 1)
(119, 37)/(1, 1)
10
(97, 41)
(97, 41)/(0, 0)
(96, 41)/(1, 0)
(96, 41)/(1, 0)
(96, 41)/(1, 0)
11
(192, 396)
(33, 398)/(159, 2)
(192, 396)/(0, 0)
(191, 396)/(1, 0)
(192, 396)/(0, 0)
12
(55, 306)
(55, 306)/(0, 0)
(55, 306)/(0, 0)
(55, 306)/(0, 0)
(55, 306)/(0, 0)
Mean error
(—/—)
(19.7, 5.4)
(0.6, 0.7)
(0.7, 0.8)
(0.6, 0.5)
The bold data means mismatching pairs.Matching curves of the BMI (the 4th row) compared with the NMI (the 1st row), GNMI (the 2nd row), and GCMI (the 3rd row).
(a)
(b)Figure 10
Matching results of the BMI.Additional experiments have also been done to test the robustness according to the parameter settings, includinga, h1, and h2, which refer to (10), (11), and (14), respectively. The results are listed below.(1)
Whena<0.8, the noise increases and it makes the matching curve of BMI not so smooth, but the main peak of the curve remains unchanged. When a>0.8, the salient regions shrink, and the position accuracy decreases.(2)
The change ofh1 has no substantial influence on f^Ruv(θ), as long as h1<h2.(3)
h2 determines the cover range of the likelihood, and the higher h2 is, the smoother the curve will be and, hence, the more robust the BMI will be.


## 5. Conclusions

In this paper, a robust Bayesian estimated mutual information method for multisensor images matching is proposed. Our method is universal as long as the prior and the likelihood can be modeled, and we have implemented it by using the prior information of gradient, in which the prior probability density is estimated by the KDE method, and the likelihoods are modeled, respectively, according to the class of the matching pixel pair. Our method is compared with the standard NMI and GNMI by experimenting on some multisensor images, and the results show that (1) our method provides higher robustness and position accuracy than NMI; (2) as the keys of BMI can be precalculated offline, our method is remarkably cost saving; (3) the parameters of the BMI can be set easily, which facilitate its application.Further studies include improving the matching accuracy by suppressing the speckle noise and refining the saliency mask. Moreover, to enhance its robustness, additional information can be introduced according to the practical conditions.



---

*Source: 414785-2013-10-31.xml*