# A Subpixel Matching Method for Stereovision of Narrow Baseline Remotely Sensed Imagery

Authors: Ning Ma; Peng-fei Sun; Yu-bo Men; Chao-guang Men; Xiang Li
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

In this paper, an accurate and efficient image matching method based on phase correlation is proposed to estimate disparity with subpixel precision, which is used for the stereovision of narrow baseline remotely sensed images. The multistep strategy is adopted in our technical frame; thus the disparity estimation is divided into two steps: integer-pixel prematching and subpixel matching. Firstly, integer-pixel disparity is estimated by employing a cross-based local matching method. Then the relationship of corresponding points is established under the guidance of integer-pixel disparity. The subimages are extracted through selecting the corresponding points as the center. Finally, the subpixel disparity is obtained by matching the subimages utilizing a novel variant of phase correlation approach. The experiment results show that the proposed method can match different kinds of large-sized narrow baseline remotely sensed images and estimate disparity with subpixel precision automatically.


---

## Body

## 1. Introduction

Stereovision is an advanced task in remote sensing and photogrammetry [1]. The aim of stereovision is to estimate the disparity through matching two or more images of same scene in different views and extract digital elevation model (DEM) through the disparity [2]. Intuitively, the disparity represents the displacement vectors between corresponding pixels that horizontally shift from the left image to the right image [3]. The stereovision of narrow baseline remotely sensed imagery is a new research hotspot for stereovision in recent years [4]. Automatic subpixel image matching is the essential technique for narrow baseline stereovision [5]. We will brief review the stereovision with different baseline and subpixel image matching method in the remaining part of this section, respectively.


### 1.1. Principle of Stereovision

The process of DEM extracting can be represented by a stereovision model [6], as shown in Figure 1.B andH express the baseline and the altitude of a satellite.L andR are the left and right images captured by the satellite from different views;h represents the height of ground object to be measured.A andB are the 3D points of the real world scene, which are projected to the three 2D locationsa, a′, and b′ inL andR. The 2D points ofA andB inL are coincident; thus the disparityd of the two corresponding pointsa and a′ equals the distance between a′ and b′ inR.D represents the mapping ofd in the 3D scene; according to the similar triangles of geometry principle, the equation D/B=h/(H-h) is established. Assume that the ground sample distance (GSD) of the sensed image isG (meters/pixel); the heighth of ground object can be generated: h=dG/(B/(H-h)). Because the altitude of satellite is much larger than the height of ground object, that is, H≫h, the heighth can be approximated as h≈dG/(B/H). It can be known that the precision ofh is proportional to the precision of disparity and GSD and is inversely proportional to theB/H ratio. In most real sensors,H andG are limited and fixed; therefore, only two aspects can be used to enhance the precision ofh: (1) Improve the precision of disparity to subpixel level [7]; (2) increase the length of baseline, that is, increasing theB/H ratio. In traditional stereovision of remotely sensed imagery,B/H ratio is usually set to about 1, called stereovision of wide baseline remotely sensed imagery. Even if the disparity precision is only up to the integer-pixel level, the DEM can be extracted accurately.Figure 1
Stereovision model.


### 1.2. Comparison of Wide Baseline and Narrow Baseline Stereo

The stereovision of wide baseline remotely sensed imagery can not perform well in urban areas with dense buildings [8]. The captured process of the wide baseline remotely sensed image pair is shown in Figure 2(a). Due to the long baseline and large viewing angle, there is a large occlusion area and large geometric distortion between the images. In addition, the interval time of two imaging processes is long; the illumination condition has changed greatly. The above interference factors will lead to a high false matching rate in stereovision of wide baseline remotely sensed imagery. Therefore, it requires a lot of manual and postcorrection to extract DEM [9]. In order to solve this problem of DEM extraction in urban areas, stereovision of narrow baseline remotely sensed imagery technique has emerged and gradually become a new research focus in stereovision and remote sensing [10, 11]. The captured process of the narrow baseline remotely sensed image pair is shown in Figure 2(b). Due to the short baseline and small viewing angle, there is a small occlusion area and little geometric distortion between images. When B/H→0 (i.e., a very narrow baseline), then occlusions will be minimized, essentially because each image view becomes more geometrically similar. In addition, the interval time of two imaging processes is short; it can be considered that the illumination condition is similar to the same time. Due to above advantages, the stereovision of narrow baseline remotely sensed imagery is suitable for extracting DEM in the urban areas [8]. Two sets of simulated remotely sensed image pairs provided by Beijing Institute of Space Mechanics and Electricity are shown in Figure 3. Figure 3(a) is simulated wide baseline remotely sensed image pair with aB/H ratio of 1, and Figure 3(b) is simulated narrow baseline remotely sensed image pair with aB/H ratio of 0.05. From the comparison of Figures 3(a) and 3(b), we can see that the occlusion range and the geometric distortion of narrow baseline image pair are less, and there is a higher degree of similarity between left and right image. It is beneficial to improve the accuracy of image matching; thus fully automatic DEM extraction can be achieved through stereovision of narrow baseline remotely sensed imagery.Figure 2
The captured process of the remotely sensed image pair. (a) Traditional wide baseline stereo configuration and (b) narrow baseline stereo configuration.
(a) (b)Figure 3
Two sets of simulated remotely sensed image pairs provided by Beijing Institute of Space Mechanics and Electricity. (a) Simulated wide baseline remotely sensed image pair with aB/H ratio of 1 and (b) simulated narrow baseline remotely sensed image pair with aB/H ratio of 0.05.
(a) (b)However the reduction ofB/H ratio will inevitably lead to the reduction of precision for DEM extracting. If the stereovision of narrow baseline is going to extract the DEM as precise as the stereovision of wide baseline, the disparity needs to achieve 1/T pixel precision, where T is the ratio of wide baseline and narrow baseline [6]. For example, whenB/H is 0.05, the precision of disparity needs to achieve 1/20 pixels to satisfy the DEM extraction. Therefore, the key technique for the stereovision of narrow baseline remotely sensed imagery is the high precision subpixel image matching [5, 7].


### 1.3. A Survey of Subpixel Image Matching Method

Automatic subpixel image matching is one of the most essential techniques in stereovision. It can be divided into three major categories: interpolation based method, fitting based method, and phase correlation method [12]. Interpolation based method includes the original image interpolation method and the matching cost interpolation method. The matching cost interpolation method is the representative algorithm, which combines the advantages of high efficiency and accuracy [13, 14]. The initial cost volume is interpolated by various interpolation function, and the subpixel disparity is estimated by searching the extremum of interpolated cost volume. At present, this method is usually used for the auxiliary system of driverless vehicles, robot navigation, and unmanned aerial vehicle. The original image interpolation method significantly increases the computing load, limits the possible precision to the chosen upsampling rate, and also may introduce interpolation artifacts [5]. The fitting based method generally postprocess the cost volume or disparity plane by fitting method to estimate the subpixel disparity. The cost volume fitting method fits the peak neighborhood of cost volume into a parabola; then the subpixel disparity is estimated by searching the extreme of parabola [15, 16]. The disparity plane fitting method models the disparity field by segmentation constraints, and subpixel disparity is obtained by fitting disparity plane [17–20]. The fitting based method is efficient, but the disparity precision is low. Phase correlation method provides high efficiency and accuracy via fast Fourier transform and other supplementary approaches under ideal conditions. In general, the main peak location of the inverse Fourier transform of the normalized cross-power spectrum was interpolated with the parabolic, Gaussian, and sinc functions to get the subpixel disparity [21–24].In this paper, a subpixel image matching method based on phase correlation is proposed to estimate the disparity with subpixel precision for stereovision of narrow baseline remotely sensed imagery. The rest of this paper is organized as follows. The principle of phase correlation is introduced briefly in Section2. A novel improved phase correlation method with subpixel precision is described in detail in Section 3. The complete algorithm is described in Section 4. Experimental evaluation is presented in Section 5. Finally, the paper is concluded in Section 6.


## 1.1. Principle of Stereovision

The process of DEM extracting can be represented by a stereovision model [6], as shown in Figure 1.B andH express the baseline and the altitude of a satellite.L andR are the left and right images captured by the satellite from different views;h represents the height of ground object to be measured.A andB are the 3D points of the real world scene, which are projected to the three 2D locationsa, a′, and b′ inL andR. The 2D points ofA andB inL are coincident; thus the disparityd of the two corresponding pointsa and a′ equals the distance between a′ and b′ inR.D represents the mapping ofd in the 3D scene; according to the similar triangles of geometry principle, the equation D/B=h/(H-h) is established. Assume that the ground sample distance (GSD) of the sensed image isG (meters/pixel); the heighth of ground object can be generated: h=dG/(B/(H-h)). Because the altitude of satellite is much larger than the height of ground object, that is, H≫h, the heighth can be approximated as h≈dG/(B/H). It can be known that the precision ofh is proportional to the precision of disparity and GSD and is inversely proportional to theB/H ratio. In most real sensors,H andG are limited and fixed; therefore, only two aspects can be used to enhance the precision ofh: (1) Improve the precision of disparity to subpixel level [7]; (2) increase the length of baseline, that is, increasing theB/H ratio. In traditional stereovision of remotely sensed imagery,B/H ratio is usually set to about 1, called stereovision of wide baseline remotely sensed imagery. Even if the disparity precision is only up to the integer-pixel level, the DEM can be extracted accurately.Figure 1
Stereovision model.


## 1.2. Comparison of Wide Baseline and Narrow Baseline Stereo

The stereovision of wide baseline remotely sensed imagery can not perform well in urban areas with dense buildings [8]. The captured process of the wide baseline remotely sensed image pair is shown in Figure 2(a). Due to the long baseline and large viewing angle, there is a large occlusion area and large geometric distortion between the images. In addition, the interval time of two imaging processes is long; the illumination condition has changed greatly. The above interference factors will lead to a high false matching rate in stereovision of wide baseline remotely sensed imagery. Therefore, it requires a lot of manual and postcorrection to extract DEM [9]. In order to solve this problem of DEM extraction in urban areas, stereovision of narrow baseline remotely sensed imagery technique has emerged and gradually become a new research focus in stereovision and remote sensing [10, 11]. The captured process of the narrow baseline remotely sensed image pair is shown in Figure 2(b). Due to the short baseline and small viewing angle, there is a small occlusion area and little geometric distortion between images. When B/H→0 (i.e., a very narrow baseline), then occlusions will be minimized, essentially because each image view becomes more geometrically similar. In addition, the interval time of two imaging processes is short; it can be considered that the illumination condition is similar to the same time. Due to above advantages, the stereovision of narrow baseline remotely sensed imagery is suitable for extracting DEM in the urban areas [8]. Two sets of simulated remotely sensed image pairs provided by Beijing Institute of Space Mechanics and Electricity are shown in Figure 3. Figure 3(a) is simulated wide baseline remotely sensed image pair with aB/H ratio of 1, and Figure 3(b) is simulated narrow baseline remotely sensed image pair with aB/H ratio of 0.05. From the comparison of Figures 3(a) and 3(b), we can see that the occlusion range and the geometric distortion of narrow baseline image pair are less, and there is a higher degree of similarity between left and right image. It is beneficial to improve the accuracy of image matching; thus fully automatic DEM extraction can be achieved through stereovision of narrow baseline remotely sensed imagery.Figure 2
The captured process of the remotely sensed image pair. (a) Traditional wide baseline stereo configuration and (b) narrow baseline stereo configuration.
(a) (b)Figure 3
Two sets of simulated remotely sensed image pairs provided by Beijing Institute of Space Mechanics and Electricity. (a) Simulated wide baseline remotely sensed image pair with aB/H ratio of 1 and (b) simulated narrow baseline remotely sensed image pair with aB/H ratio of 0.05.
(a) (b)However the reduction ofB/H ratio will inevitably lead to the reduction of precision for DEM extracting. If the stereovision of narrow baseline is going to extract the DEM as precise as the stereovision of wide baseline, the disparity needs to achieve 1/T pixel precision, where T is the ratio of wide baseline and narrow baseline [6]. For example, whenB/H is 0.05, the precision of disparity needs to achieve 1/20 pixels to satisfy the DEM extraction. Therefore, the key technique for the stereovision of narrow baseline remotely sensed imagery is the high precision subpixel image matching [5, 7].


## 1.3. A Survey of Subpixel Image Matching Method

Automatic subpixel image matching is one of the most essential techniques in stereovision. It can be divided into three major categories: interpolation based method, fitting based method, and phase correlation method [12]. Interpolation based method includes the original image interpolation method and the matching cost interpolation method. The matching cost interpolation method is the representative algorithm, which combines the advantages of high efficiency and accuracy [13, 14]. The initial cost volume is interpolated by various interpolation function, and the subpixel disparity is estimated by searching the extremum of interpolated cost volume. At present, this method is usually used for the auxiliary system of driverless vehicles, robot navigation, and unmanned aerial vehicle. The original image interpolation method significantly increases the computing load, limits the possible precision to the chosen upsampling rate, and also may introduce interpolation artifacts [5]. The fitting based method generally postprocess the cost volume or disparity plane by fitting method to estimate the subpixel disparity. The cost volume fitting method fits the peak neighborhood of cost volume into a parabola; then the subpixel disparity is estimated by searching the extreme of parabola [15, 16]. The disparity plane fitting method models the disparity field by segmentation constraints, and subpixel disparity is obtained by fitting disparity plane [17–20]. The fitting based method is efficient, but the disparity precision is low. Phase correlation method provides high efficiency and accuracy via fast Fourier transform and other supplementary approaches under ideal conditions. In general, the main peak location of the inverse Fourier transform of the normalized cross-power spectrum was interpolated with the parabolic, Gaussian, and sinc functions to get the subpixel disparity [21–24].In this paper, a subpixel image matching method based on phase correlation is proposed to estimate the disparity with subpixel precision for stereovision of narrow baseline remotely sensed imagery. The rest of this paper is organized as follows. The principle of phase correlation is introduced briefly in Section2. A novel improved phase correlation method with subpixel precision is described in detail in Section 3. The complete algorithm is described in Section 4. Experimental evaluation is presented in Section 5. Finally, the paper is concluded in Section 6.


## 2. Phase Correlation

Phase correlation method is based on the well-known Fourier-domain shift property. This states that a translation between two images in the spatial domain will be expressed in the frequency domain as a linear phase difference between their Fourier transform [5]. Let f1 and f2 be two images with translation relationship, the size of image is M×N, f1=f(x,y), f2=f(x+Δx,y+Δy), and (Δx,Δy) are horizontal and vertical translations, respectively. Let F1(u,v) and F2(u,v) denote the 2D discrete Fourier transform (DFT) of two images:(1)F1u,v=1MN∑x=0M-1∑y=0N-1fx,yexp⁡-j2πuxM+vyN,F2u,v=F1u,vexp⁡-j2πuΔxM+vΔyN,where x and y are spatial domain variables, x=0,1,2,…,M-1, y=0,1,2,…,N-1. u and v are frequency domain variables, u=0,1,2,…,M-1, v=0,1,2,…,N-1. The normalized cross-power spectrum Q(u,v) of two images is calculated as(2)Qu,v=F1u,vF2∗u,vF1u,vF2∗u,v=exp⁡-j2πuΔxM+vΔyN,where F2∗(u,v) is complex conjugate of F2(u,v). The phase correlation PC(x,y) of two images is calculated by 2D inverse discrete Fourier transforms (IDFT) of Q(u,v):(3)PCx,y=F-1Qu,v=δx+Δx,y+Δy,where F-1(·) is IDFT function. δ(x+Δx,y+Δy) is a 2D pulse function which exists as peak at position (Δx,Δy). The peak location (Δx,Δy) corresponds to the translation between images.Phase correlation is a common method to estimate the translation between images. Only phase information of images is utilized; thus it is almost not affected by noise and radiation difference and more robustness than spatial domain method. But the digital image is a discrete function whose Fourier transform is a discrete transform, so only the integer-pixel translation(int(Δx),int(Δy)) can be obtained by searching the peak location.


## 3. Improved Phase Correlation with Subpixel Precision

To estimate the subpixel translation, various algorithms have been proposed. In [23], Foroosh et al. achieved the subpixel estimation by approximating the Dirichlet function derived from the normalized cross-power spectrum. Nagashima et al. propose a fitting method to enhance the matching precision from integer-pixel level to subpixel level, which fits a curve surface by using the neighborhood data around phase correlation peak [21]. However, these methods only take into account the peak and its neighborhood; the distribution of phase correlation function is neglected; thus, the precision of translation estimation is low in stereovision of remotely sensed imagery. An improved phase correlation method used for subpixel estimation is designed in this section. The flow chart is shown in Figure 4. Due to the periodicity of DFT, an image can be considered to “wrap around” at an edge; therefore, discontinuities, which are not supposed to exist in real world, occur at every border in 2D DFT computation [21]. A 2D Hanning window function is applied to reduce the effect of discontinuities at image borders. The 2D Hanning window function is defined as (4)wx,y=1+cos⁡πx/M21+cos⁡πy/N2.Figure 4
Flow chart of improved phase correlation image matching method.The subpixel translation between two images can be obtained by downsampling high resolution images with integer-pixel translation. Therefore, the phase correlation of two images can be approximated by 2D sinc function:(5)PCx,y≈sin⁡πx+Δxπx+Δxsin⁡πy+Δyπy+Δy.In theory, sinc function can approximate the distribution of phase correlation in ideal condition. However the phase correlation peak is decreased in remotely sensed images matching due to the interference of random noise, image aliasing, edge effects, and other influences. Therefore, the sinc function can not describe the distribution of phase correlation accurately. In our method, the sinc function is improved by introducing phase correlation coefficient:(6)PCx,y≈αsin⁡πx+Δxπx+Δxsin⁡πy+Δyπy+Δy,where α is the phase correlation coefficient, α≤1. When the translation is integer-pixel and there is no interference factor, α=1; when the translation is subpixel or there are interference factors, α<1. Due to this improved phase correlation method is applied to the stereovision of narrow baseline remotely sensed imagery, and the image pair is treated by epipolar rectification before matching. Thus, the phase correlation PC(x,y) can be separated in spatial domain, and a 1D function expression is given which only exists as translation in horizontal direction:(7)PCxx≈αsin⁡πx+Δxπx+Δx.The peak position(Δx,0) is the translation which is expected to be estimated. To locate the high precision subpixel peak position, a peak evaluation method based on uniformly spaced sampling is proposed. Assume that p is the real peak position of sampled-data, as shown at the bottom left corner of Figure 4, and p-k and p+k are a pair of points located k pixels away from p. The phase correlation values PCx(p-k), PCx(p), and PCx(p+k) of three sample points used for peak evaluation are satisfying the following relationship:(8a)αsinπp-k+Δxπp-k+Δx=PCxp-k(8b)αsinπp+Δxπp+Δx=PCxp(8c)αsinπp+k+Δxπp+k+Δx=PCxp+k.Add (8c) to (8a); the above equation group can be rewritten as(9a)απsin⁡πp+Δx=PCxpp+Δx(9b)απsin⁡πp-k+Δx+sin⁡πp-k+Δx=PCxp-kp-k+Δx+PCxp+kp+k+Δx.Equation (9b) can be rewritten by using the difference product formula: (10a)απsin⁡πp+Δx=PCxpp+Δx(10b)απ2sin⁡πp+Δxcos⁡πk=PCxp-kp-k+Δx+PCxp+kp+k+Δx.Add (10a) to (10b); the equation group ((10a) and (10b)) can be rewritten as(11)PCxp-kp-k+Δx+PCxp+kp+k+Δx=2·PCxpp+Δxcosπk.The above equation can be expressed as a linear expression ofΔx:(12)PCxp-k+PCxp+k-2·cos⁡πkPCxp·Δx=2·p·cos⁡πkPCxp-p-kPCxp-k+p+kPCxp+k.Through (12), the horizontal translation Δx with subpixel precision can be estimated:(13)Δx=2·p·cos⁡πkPCxp-p-kPCxp-k+p+kPCxp+kPCxp-k+PCxp+k-2·cos⁡πkPCxp.To reduce the influence of noise, the peak evaluation is performed by selecting multiple sets of uniformly spaced sampled-data which are all the data on thex axis, and the least square is applied to solve the optimal translation. Assume that S sets of uniformly spaced sampled-data (pi-ki,pi,pi-ki)i∈S are selected for peak evaluation; pi, ki are determined as follows: p1=p2=⋯=pS=argmaxx⁡PCx(x), k1=1, k2=k1+Δk,…,kS=kS-1+Δk, (Δk=1). Equation (12) can be rewritten as(14)ap1,k1ap2,k2⋮apS,kS·ΔxΔx⋮Δx=bp1,k1bp2,k2⋮bpS,kS,where a(pi,ki)i∈S=PCx(pi-ki)+PCx(pi+ki)-2·cos⁡(πki)PCx(pi), and b(pi,ki)i∈S=2·pi·cos(πki)PCx(pi)-(pi-ki)PCx(pi-ki)+pi+kiPCx(pi+ki). The objective of (14) is min⁡(∑i=1Sb(pi,ki)-a(pi,ki)·Δx2). Assume that A=a(p1,k1),a(p2,k2),…,a(pS,kS)T and B=b(p1,k1),b(p2,k2),…,b(pS,kS)T, then the horizontal translation Δx with subpixel precision can be obtained by solving the least square linear equation:(15)Δx=ATA-1ATB.


## 4. Complete Subpixel Image Matching Algorithm

The input of our algorithm is rectified narrow baseline remotely sensed image pair, and output is the disparity with subpixel precision. Because the stability of phase correlation method is poor when the disparity search range is large [9], therefore, a multistep strategy is adopted in our technical frame, and the disparity estimation is divided into two steps: integer-pixel prematching and subpixel matching. The complete framework of our algorithm is shown in Figure 5. Firstly, integer-pixel disparity is estimated by employing an effective cross-based matching algorithm [25]. Then relationship of corresponding points is established under the guidance of integer-pixel disparity. The subimages are extracted through selecting the corresponding points as center. Finally, subpixel disparity is obtained by matching the subimages utilizing the improved phase correlation proposed in Section 3.Figure 5
The complete framework of our subpixel image matching algorithm.Assume that a narrow baseline remotely sensed image pair consists of left imageL and right image R, which can be expressed as a classic stereo model L(x,y)=R(x+d(x,y),y). The disparity d(x,y) is used to describe the geometric position variation between L and R, which consists of two parts, integer-pixel disparity dz(x,y) and subpixel disparity ds(x,y); that is, d(x,y)=dz(x,y)+ds(x,y). The details of our subpixel image matching algorithm are described as follows from step (a) to step (g).(a) Assume(x0,y0) is the pixel which is to be matched. In the first step, the cross-based support window of (x0,y0) is defined. Firstly, an upright cross is established for (x0,y0). It consists of two orthogonal line segments: horizontal segment φH(x0,y0) and vertical segment φV(x0,y0). A quadruple {h-,h+,v-,v+} denotes the left, right, up, and bottom arm length, respectively, and the length of each arm is determined by searching for the extreme point in that direction based on the color similarity. The upright cross skeleton is defined as(16)φHx0,y0=x,y∣x∈x0-h-,x0+h+,y=y0,φVx0,y0=x,y∣x=x0,y∈y0-v-,y0+v+.Then the cross-based support windowφ(x0,y0) is defined based on the upright cross skeleton. φ(x0,y0) is a combination of each horizontal segment φH(x,y) where (x,y) traverses over the vertical segment φV(x0,y0). The formula of φ(x0,y0) is defined as(17)φx0,y0=⋃x,y∈φVx0,y0φHx,y.(b) The matching cost of(x0,y0) is calculated. Firstly, the initial matching cost Cx,y,dz of each pixel in the support window is calculated:(18)Cx,y,dz=1-β·min⁡Lx,y-Rx+dz,y,τ1+β·min⁡∇xLx,y-∇xRx+dz,y,τ2,where ∇x is horizontal gradient, β balances the color and gradient terms, and τ1, τ2 are truncation values. Then the matching cost Cx0,y0,dzA of (x0,y0) is aggregated:(19)Cx0,y0,dzA=∑x,y∈φx0,y0Cx,y,dz.(c) The integer-pixel disparitydz(x0,y0) of (x0,y0) is estimated by winner-takes-all strategy:(20)dzx0,y0=argmindz∈SdCx0,y0,dzA,where Sd is the searching range of disparity.(d) The relationship of corresponding points is established under the guidance of integer-pixel disparity. Assume that the integer-pixel disparity of pixel(x0,y0) in L is dz(x0,y0); the corresponding point of (x0,y0) in R is (x0+dz(x0,y0),y0).(e) The subimagesLs(x0,y0) and Rs(x0,y0) with same size are extracted through selecting the corresponding points as the center. Because the selection of the two subimages is based on the guidance of the integer-pixel disparity, thus there is only a small range of horizontal translation between Ls(x0,y0) and Rs(x0,y0).(f) The subpixel horizontal translationΔx between Ls(x0,y0) and Rs(x0,y0) is estimated by using the improved phase correlation method proposed in Section 3. The Δx is equivalent to the subpixel disparity ds(x0,y0) of the center pixel (x0,y0).(g) Finally, the disparity result is obtained by combining the integer-pixel disparity and the subpixel disparity; that is,d(x0,y0)=dz(x0,y0)+ds(x0,y0).


## 5. Experiments and Analysis

We test the proposed subpixel matching method for stereovision of narrow baseline remotely sensed imagery on a standard personal computer with Intel(R) Core(TM) i7 CPU, and the algorithm is implemented utilizing VS2010+OpenCV. The results analysis includes two aspects: precision analysis and complexity analysis. Constant parameter settings are used for all experiments. The parameters of this method are set asβ=0.11, τ1=7/255, and τ2=2/255.


### 5.1. Test on Simulated Narrow Baseline Remotely sensed Images Generated by Translation

To evaluate the matching precision of proposed method, six urban scenes of narrow baseline remotely sensed image pairs, including Pleiades, SPOT-5, SPOT-6, WorldView-3, and GF-2, are used to perform image matching. For each image pair, the original image is utilized as the reference image (often called the left image), and the original image is moved 8.738 pixels along thex-axis as the target image (called the right image). Thus the true disparity is known to be 8.738 pixels for all the matching points of each image pair. Table 1 shows the general information of the test image pairs used in these experiments.Table 1
General information of the test image pairs.
Test number Data source Image size GSD (m) Location True disparity (pixels) 1 Pleiades 7590 × 6510 0.5 USA: Washington 8.738 2 SPOT-5 2050 × 2100 2.5 China: Beijing 8.738 3 SPOT-6 2015 × 1950 1.5 China: Beijing 8.738 4 WorldView-3 2000 × 2000 0.3 Spain: Madrid 8.738 5 GeoEye-1 3000 × 3000 0.4 Holland: Amsterdam 8.738 6 GF-2 4500 × 4750 0.8 China: Karamay 8.738Figure6 shows the distribution of matched points for the test narrow baseline remotely sensed image pairs. Because the image size is different for each image pair, the step size of grid is different too. The step sizes of grid in Tests 1–6 are set to 540 × 540 pixels, 160 × 160 pixels, 180 × 180 pixels, 170 × 170 pixels, 300 × 300 pixels, and 380 × 380 pixels.Figure 6
The distribution of matched points for the test narrow baseline remotely sensed image pairs. (a–f) are the results of Pleiades, SPOT-5, SPOT-6, WorldView-3, GeoEye-1, and GF-2 image pairs, respectively.
(a) (b) (c) (d) (e) (f)In order to evaluate the matching precision of our method quantitatively, interpolation based method [14], fitting based method [15], and traditional phase correlation method [23] are used to compare with our method. Root mean square error (RMSE) and mean error of matched points are used as the assessment indices to evaluate these methods. Figure 7 shows quantitative evaluation results for the four subpixel image matching methods. Here the traditional phase correlation is called traditional PC for short. Figure 8 shows statistical results of the quantitative evaluation. Figure 8(a) shows statistical results where the RMSE is less than or equal to 0.05 pixels. Figure 8(b) shows statistical results where the RMSE is greater than 0.05 pixels and less than 0.1 pixels. Figure 8(c) shows statistical results where the RMSE is greater than or equal to 0.1 pixels. From the experimental results, one can see that the two PC methods have a high matching precision, and the precision is better than that of interpolation method and fitting method. In addition, the matching precision is related to the spatial resolution of the images. High spatial resolution corresponds to a high matching precision. Because the peak evaluation method based on uniformly spaced sampling is designed to fit the phase correlation function, our method outperforms the traditional PC method, about 66% matched points of our method with RMSE less than or equal to 0.05 pixels.Figure 7
The quantitative evaluation results for the four subpixel image matching methods. (a) RMSE and (b) mean error.
(a) (b)Figure 8
The statistical results of the quantitative evaluation. (a) Statistical results of RMSE less than or equal to 0.05 pixels; (b) statistical results of RMSE greater than 0.05 pixels and less than 0.1 pixels; and (c) statistical results of RMSE greater than or equal to 0.1 pixels.
(a) (b) (c)


### 5.2. Test on Simulated Narrow Baseline Remotely sensed Images

To evaluate the precision of DEM extraction by using the disparity estimated by our method, a simulated narrow baseline remotely sensed image pair provided by Beijing Institute of Space Mechanics and Electricity is used to perform image matching. The simulated narrow baseline remotely sensed image pair is generated by SE-WORKBENCH software. The attached parameters file of the image pair not only provides the true disparity, but also provides the true elevation of targets and buildings in the scene. Therefore, the precision of disparity and elevation estimated by our method can be evaluated carefully. Table2 shows the general information of the test image pair used in these experiments. TheB/H ratio of the image pair is 0.05, and spatial resolution is 0.3 meters/pixel. Overlap ratio of left and right images is 60%, and the overlap area includes 22 targets and 40 buildings. Figure 9 shows the simulated narrow baseline remotely sensed image pair and distribution of targets and buildings. Figure 9(a) is the image pair, Figure 9(b) is the distribution of targets, and Figure 9(c) is the distribution of buildings. The targets and buildings are marked with red numbers. Figure 10 shows the matching results. Figure 10(a) is the distribution of matched points for the buildings. Figure 10(b) is the distribution of matched points for the targets.Table 2
The general information of the test image pair.
Target number Elevation (meter) True disparity (pixel) Building number Elevation (meter) True disparity Building number Elevation (meter) True disparity (pixel) 1 73 12.167 1 111.42 18.570 21 97.3 16.217 2 65 10.833 2 131.54 21.923 22 131.33 21.888 3 63 10.500 3 147.13 24.522 23 165.92 27.653 4 65 10.833 4 153.77 25.628 24 134.85 22.475 5 67 11.167 5 123.43 20.572 25 85.14 14.190 6 61 10.167 6 104.77 17.462 26 92.28 15.380 7 67 11.167 7 98.22 16.370 27 71.51 11.918 8 82 13.667 8 88.81 14.802 28 83.5 13.917 9 74 12.333 9 75.07 12.512 29 169.2 28.200 10 72 12.000 10 98.46 16.410 30 131.22 21.870 11 75 12.500 11 79.79 13.298 31 98.42 16.403 12 65 10.833 12 72.94 12.157 32 98.42 16.403 13 60 10.000 13 155.81 25.968 33 154.57 25.762 14 64 10.667 14 111.81 18.635 34 115.5 19.250 15 63 10.500 15 133.81 22.302 35 105.21 17.535 16 61 10.167 16 100.81 16.802 36 133.07 22.178 17 76 12.667 17 92.21 15.368 37 157.59 26.265 18 64 10.667 18 120.38 20.063 38 66.3 11.050 19 62 10.333 19 162.39 27.065 39 72.5 12.083 20 66 11.000 20 145.29 24.215 40 72.44 12.073 21 92 15.333 — — — — — — 22 55 9.167 — — — — — —Figure 9
The simulated narrow baseline remotely sensed image pair and distribution of targets and buildings. (a) The image pair with aB/H ratio of 0.05. (b) The distribution of targets. (c) The distribution of buildings.
(a) (b) (c)Figure 10
The matching results. (a) The distribution of matched points for the buildings. (b) The distribution of matched points for the targets.
(a) (b)To evaluate the matching precision of our method quantitatively, interpolation based method [14], fitting based method [15], and traditional phase correlation method [23] are used to compare with our method. RMSE of disparity and mean elevation error are used as the assessment indices to evaluate the precision of four methods, where the elevation error is the absolute difference between the true and the extracted elevation. The simulated narrow baseline remotely sensed image pair is matched 300 times by each subpixel image matching method. Figures 11 and 12 show the statistical results of four methods for targets and buildings, respectively. The results show that the RMSE and the mean elevation error obtained by our method are lower than the other three methods for most targets and buildings. Table 3 shows the average values of RMSE and mean elevation error. The true average elevation of 22 targets is 67.818 m, and the true average elevation of 40 buildings is 112.353 m. According to the statistical results in Table 3, the mean elevation error of targets estimated by our method is 0.303 m, and the mean elevation error of buildings estimated by our method is 0.406 m. Such precision can satisfy the requirements for stereovision of narrow baseline remotely sensed imagery; therefore, it is proved that the image matching method based on phase correlation proposed in this paper has practical application value.Table 3
The RMSE and mean elevation error of four methods.
Algorithms Targets Buildings RMSE (pixel) Mean elevation error (meter) RMSE (pixel) Mean elevation error (meter) Interpolation method [14] 0.103 0.619 0.114 0.681 Fitting method [15] 0.169 1.012 0.133 0.797 Traditional PC method [23] 0.061 0.368 0.071 0.422 Our method 0.051 0.303 0.068 0.406Figure 11
The statistical results of targets. (a) Comparison of RMSE. (b) Comparison of mean elevation error.
(a) (b)Figure 12
The statistical results of buildings. (a) Comparison of RMSE. (b) Comparison of mean elevation error.
(a) (b)


### 5.3. Computational Complexity

Assume that there areN1 points needed to be matched. In integer-pixel prematching step of our algorithm, the complexity of constructing the cross-based support window for each point is O(1), and the total complexity for all points of this prematching step is ON1·d, where d is the disparity range.In the subpixel matching step, the core algorithm is based on the phase correlation; therefore, the complexity is unrelated to the disparity range but related to the size of subimage. Assume that there areN2 pixels in the extracted subimage. The complexity of fast DFT is ON2·log⁡N2, the complexity of Hanning window is ON2, the complexity of normalized cross-power spectrum is ON2, and the complexity of IDFT of normalized cross-power spectrum is ON2·log⁡N2. The peak evaluation is a least square fitting operation for S sets of uniformly spaced sampled-data essentially. The minimum value of S is 1; that is, the peak position is located at the most left or right of sampled-data. However, the maximum value of S is (N2-1)/2; this means the peak position is located at the center of sampled-data; thus S≤(N2-1)/2. The complexity of subpixel peak position evaluation utilizing the least square method is O(N2-1)/22, which is equivalent to ON2. The complexity of a single point matching in subpixel matching step is ON2·log⁡N2; thus the complexity of total points is ON1·N2·log⁡N2. The complexity of complete image matching algorithm is ON1·d+ON1·N2·log⁡N2. Because the complexity of subpixel matching is higher than that of integer-pixel prematching, the complexity of our algorithm is ON1·N2·log⁡N2.In order to clearly illustrate the complexity of our algorithm, we compare with the classical normalized cross correlation (NCC) algorithm [2]. Assume that the support window of NCC is φ; in the practical image matching, φ≈N2≪N1. When log⁡N2≈d, the complexity of our method is similar to that of NCC. Generally, the size of subimage is 30 × 30~40 × 40 pixels; the number of pixels in subimage is 900~1600; thus the value range of log⁡N2 is 9.814~10.644. The disparity range of narrow baseline stereo image pair is usually within 20 pixels. In our experiment, the disparity range is 17 pixels; log⁡N2≈d is established; therefore, the complexity of our method is similar to that of the classical NCC method.


## 5.1. Test on Simulated Narrow Baseline Remotely sensed Images Generated by Translation

To evaluate the matching precision of proposed method, six urban scenes of narrow baseline remotely sensed image pairs, including Pleiades, SPOT-5, SPOT-6, WorldView-3, and GF-2, are used to perform image matching. For each image pair, the original image is utilized as the reference image (often called the left image), and the original image is moved 8.738 pixels along thex-axis as the target image (called the right image). Thus the true disparity is known to be 8.738 pixels for all the matching points of each image pair. Table 1 shows the general information of the test image pairs used in these experiments.Table 1
General information of the test image pairs.
Test number Data source Image size GSD (m) Location True disparity (pixels) 1 Pleiades 7590 × 6510 0.5 USA: Washington 8.738 2 SPOT-5 2050 × 2100 2.5 China: Beijing 8.738 3 SPOT-6 2015 × 1950 1.5 China: Beijing 8.738 4 WorldView-3 2000 × 2000 0.3 Spain: Madrid 8.738 5 GeoEye-1 3000 × 3000 0.4 Holland: Amsterdam 8.738 6 GF-2 4500 × 4750 0.8 China: Karamay 8.738Figure6 shows the distribution of matched points for the test narrow baseline remotely sensed image pairs. Because the image size is different for each image pair, the step size of grid is different too. The step sizes of grid in Tests 1–6 are set to 540 × 540 pixels, 160 × 160 pixels, 180 × 180 pixels, 170 × 170 pixels, 300 × 300 pixels, and 380 × 380 pixels.Figure 6
The distribution of matched points for the test narrow baseline remotely sensed image pairs. (a–f) are the results of Pleiades, SPOT-5, SPOT-6, WorldView-3, GeoEye-1, and GF-2 image pairs, respectively.
(a) (b) (c) (d) (e) (f)In order to evaluate the matching precision of our method quantitatively, interpolation based method [14], fitting based method [15], and traditional phase correlation method [23] are used to compare with our method. Root mean square error (RMSE) and mean error of matched points are used as the assessment indices to evaluate these methods. Figure 7 shows quantitative evaluation results for the four subpixel image matching methods. Here the traditional phase correlation is called traditional PC for short. Figure 8 shows statistical results of the quantitative evaluation. Figure 8(a) shows statistical results where the RMSE is less than or equal to 0.05 pixels. Figure 8(b) shows statistical results where the RMSE is greater than 0.05 pixels and less than 0.1 pixels. Figure 8(c) shows statistical results where the RMSE is greater than or equal to 0.1 pixels. From the experimental results, one can see that the two PC methods have a high matching precision, and the precision is better than that of interpolation method and fitting method. In addition, the matching precision is related to the spatial resolution of the images. High spatial resolution corresponds to a high matching precision. Because the peak evaluation method based on uniformly spaced sampling is designed to fit the phase correlation function, our method outperforms the traditional PC method, about 66% matched points of our method with RMSE less than or equal to 0.05 pixels.Figure 7
The quantitative evaluation results for the four subpixel image matching methods. (a) RMSE and (b) mean error.
(a) (b)Figure 8
The statistical results of the quantitative evaluation. (a) Statistical results of RMSE less than or equal to 0.05 pixels; (b) statistical results of RMSE greater than 0.05 pixels and less than 0.1 pixels; and (c) statistical results of RMSE greater than or equal to 0.1 pixels.
(a) (b) (c)


## 5.2. Test on Simulated Narrow Baseline Remotely sensed Images

To evaluate the precision of DEM extraction by using the disparity estimated by our method, a simulated narrow baseline remotely sensed image pair provided by Beijing Institute of Space Mechanics and Electricity is used to perform image matching. The simulated narrow baseline remotely sensed image pair is generated by SE-WORKBENCH software. The attached parameters file of the image pair not only provides the true disparity, but also provides the true elevation of targets and buildings in the scene. Therefore, the precision of disparity and elevation estimated by our method can be evaluated carefully. Table2 shows the general information of the test image pair used in these experiments. TheB/H ratio of the image pair is 0.05, and spatial resolution is 0.3 meters/pixel. Overlap ratio of left and right images is 60%, and the overlap area includes 22 targets and 40 buildings. Figure 9 shows the simulated narrow baseline remotely sensed image pair and distribution of targets and buildings. Figure 9(a) is the image pair, Figure 9(b) is the distribution of targets, and Figure 9(c) is the distribution of buildings. The targets and buildings are marked with red numbers. Figure 10 shows the matching results. Figure 10(a) is the distribution of matched points for the buildings. Figure 10(b) is the distribution of matched points for the targets.Table 2
The general information of the test image pair.
Target number Elevation (meter) True disparity (pixel) Building number Elevation (meter) True disparity Building number Elevation (meter) True disparity (pixel) 1 73 12.167 1 111.42 18.570 21 97.3 16.217 2 65 10.833 2 131.54 21.923 22 131.33 21.888 3 63 10.500 3 147.13 24.522 23 165.92 27.653 4 65 10.833 4 153.77 25.628 24 134.85 22.475 5 67 11.167 5 123.43 20.572 25 85.14 14.190 6 61 10.167 6 104.77 17.462 26 92.28 15.380 7 67 11.167 7 98.22 16.370 27 71.51 11.918 8 82 13.667 8 88.81 14.802 28 83.5 13.917 9 74 12.333 9 75.07 12.512 29 169.2 28.200 10 72 12.000 10 98.46 16.410 30 131.22 21.870 11 75 12.500 11 79.79 13.298 31 98.42 16.403 12 65 10.833 12 72.94 12.157 32 98.42 16.403 13 60 10.000 13 155.81 25.968 33 154.57 25.762 14 64 10.667 14 111.81 18.635 34 115.5 19.250 15 63 10.500 15 133.81 22.302 35 105.21 17.535 16 61 10.167 16 100.81 16.802 36 133.07 22.178 17 76 12.667 17 92.21 15.368 37 157.59 26.265 18 64 10.667 18 120.38 20.063 38 66.3 11.050 19 62 10.333 19 162.39 27.065 39 72.5 12.083 20 66 11.000 20 145.29 24.215 40 72.44 12.073 21 92 15.333 — — — — — — 22 55 9.167 — — — — — —Figure 9
The simulated narrow baseline remotely sensed image pair and distribution of targets and buildings. (a) The image pair with aB/H ratio of 0.05. (b) The distribution of targets. (c) The distribution of buildings.
(a) (b) (c)Figure 10
The matching results. (a) The distribution of matched points for the buildings. (b) The distribution of matched points for the targets.
(a) (b)To evaluate the matching precision of our method quantitatively, interpolation based method [14], fitting based method [15], and traditional phase correlation method [23] are used to compare with our method. RMSE of disparity and mean elevation error are used as the assessment indices to evaluate the precision of four methods, where the elevation error is the absolute difference between the true and the extracted elevation. The simulated narrow baseline remotely sensed image pair is matched 300 times by each subpixel image matching method. Figures 11 and 12 show the statistical results of four methods for targets and buildings, respectively. The results show that the RMSE and the mean elevation error obtained by our method are lower than the other three methods for most targets and buildings. Table 3 shows the average values of RMSE and mean elevation error. The true average elevation of 22 targets is 67.818 m, and the true average elevation of 40 buildings is 112.353 m. According to the statistical results in Table 3, the mean elevation error of targets estimated by our method is 0.303 m, and the mean elevation error of buildings estimated by our method is 0.406 m. Such precision can satisfy the requirements for stereovision of narrow baseline remotely sensed imagery; therefore, it is proved that the image matching method based on phase correlation proposed in this paper has practical application value.Table 3
The RMSE and mean elevation error of four methods.
Algorithms Targets Buildings RMSE (pixel) Mean elevation error (meter) RMSE (pixel) Mean elevation error (meter) Interpolation method [14] 0.103 0.619 0.114 0.681 Fitting method [15] 0.169 1.012 0.133 0.797 Traditional PC method [23] 0.061 0.368 0.071 0.422 Our method 0.051 0.303 0.068 0.406Figure 11
The statistical results of targets. (a) Comparison of RMSE. (b) Comparison of mean elevation error.
(a) (b)Figure 12
The statistical results of buildings. (a) Comparison of RMSE. (b) Comparison of mean elevation error.
(a) (b)


## 5.3. Computational Complexity

Assume that there areN1 points needed to be matched. In integer-pixel prematching step of our algorithm, the complexity of constructing the cross-based support window for each point is O(1), and the total complexity for all points of this prematching step is ON1·d, where d is the disparity range.In the subpixel matching step, the core algorithm is based on the phase correlation; therefore, the complexity is unrelated to the disparity range but related to the size of subimage. Assume that there areN2 pixels in the extracted subimage. The complexity of fast DFT is ON2·log⁡N2, the complexity of Hanning window is ON2, the complexity of normalized cross-power spectrum is ON2, and the complexity of IDFT of normalized cross-power spectrum is ON2·log⁡N2. The peak evaluation is a least square fitting operation for S sets of uniformly spaced sampled-data essentially. The minimum value of S is 1; that is, the peak position is located at the most left or right of sampled-data. However, the maximum value of S is (N2-1)/2; this means the peak position is located at the center of sampled-data; thus S≤(N2-1)/2. The complexity of subpixel peak position evaluation utilizing the least square method is O(N2-1)/22, which is equivalent to ON2. The complexity of a single point matching in subpixel matching step is ON2·log⁡N2; thus the complexity of total points is ON1·N2·log⁡N2. The complexity of complete image matching algorithm is ON1·d+ON1·N2·log⁡N2. Because the complexity of subpixel matching is higher than that of integer-pixel prematching, the complexity of our algorithm is ON1·N2·log⁡N2.In order to clearly illustrate the complexity of our algorithm, we compare with the classical normalized cross correlation (NCC) algorithm [2]. Assume that the support window of NCC is φ; in the practical image matching, φ≈N2≪N1. When log⁡N2≈d, the complexity of our method is similar to that of NCC. Generally, the size of subimage is 30 × 30~40 × 40 pixels; the number of pixels in subimage is 900~1600; thus the value range of log⁡N2 is 9.814~10.644. The disparity range of narrow baseline stereo image pair is usually within 20 pixels. In our experiment, the disparity range is 17 pixels; log⁡N2≈d is established; therefore, the complexity of our method is similar to that of the classical NCC method.


## 6. Conclusion

Through the analysis of narrow baseline remotely sensed imagery stereovision, a subpixel image matching approach based on improved phase correlation is proposed. The cross-based local matching method is employed for prematching the image pair, and the obtained integer-pixel disparity reduces the search range of subpixel matching. A high precision subpixel matching step is implemented under the guidance of integer-pixel disparity. A peak evaluation method based on uniformly spaced sampling is designed to improve the precision of disparity estimation. The experimental results show that our method has superior performance on precision and efficiency, and it is beneficial to the stereovision of narrow baseline remotely sensed imagery.



---

*Source: 7901692-2017-02-23.xml*