A NOVEL SUB-PIXEL MATCHING ALGORITHM BASED ON PHASE CORRELATION USING PEAK CALCULATION

The matching accuracy of homonymy points of stereo images is a key point in the development of photogrammetry, which influences the geometrical accuracy of the image products. This paper presents a novel sub-pixel matching method phase correlation using peak calculation to improve the matching accuracy. The peak theoretic centre that means to sub-pixel deviation can be acquired by Peak Calculation (PC) according to inherent geometrical relationship, which is generated by inverse normalized cross-power spectrum, and the mismatching points are rejected by two strategies: window constraint, which is designed by matching window and geometric constraint, and correlation coefficient, which is effective for satellite images used for mismatching points removing. After above, a lot of high-precise homonymy points can be left. Lastly, three experiments are taken to verify the accuracy and efficiency of the presented method. Excellent results show that the presented method is better than traditional phase correlation matching methods based on surface fitting in these aspects of accuracy and efficiency, and the accuracy of the proposed phase correlation matching algorithm can reach 0.1 pixel with a higher calculation efficiency.


INTRODUCTION
Image matching, in the field of digital photography, is one of the important research topics (ARMIN GRUEN, 2012), and its accuracy directly restricts the development of photogrammetry to full digital photogrammetry, also influences the geometry accuracy of the subsequent geometric processing.Phase correlation matching converts stereo images to frequency domain through Fourier transform, and acquires tie points through the processing of frequency domain information (Kuglin, C.D, 1975).Compared with the traditional crosscorrelation and other high-precision image matching method, phase correlation matching has better accuracy and reliability (T.Heid, 2012).Except being applied in digital photography, as the advantage of phase correlation matching, and it has been applied to other areas such as medical imaging (W. S. Hoge), computer vision (K.Ito, 2004) and environmental change monitoring (S.Leprince, 2007), etc. Classic phase correlation matching can attain pixel precision, at present, and we can improve the pixel precision to sub-pixel precision through three kinds of optimization strategy, such as the fitting interpolation method (Kenji TAKITA, 2003), the singular value decomposition method (Xiaohua Tong, 2015) and the local upward sampling method.However, fitting function attained by the least squares estimate was affected by side lobe energy easily, which will bring amount of calculation; singular value decomposition of cross-power spectrum will cause phase unwrapping fuzzy more or less, which make it unable to attain the exact offsets for the cumulative system error.Local upward sampling method is limited by sampling ratio.However, this paper presents a high-precise sub-pixel matching method, based on symmetrical distribution of symmetry energy through peak calculation, to enhance the matching accuracy.The method is based on traditional phase correlation matching and improved, and the peak location acquired by calculation according to inherent geometry.Lastly, the mismatching points are rejected by window constraint.The related theory is simple, but it has been confirmed that the algorithm has high matching precision with less calculation through experiments of simulation data.

Classic Phase Correlation Matching
Image matching based on phase correlation method employs the Fourier transform to transform the image to be matched to frequency domain for cross-correlation.It only uses power spectrum phase information of the image blocks in the frequency domain, and reduces the effect of the image content, such as pixel value.So it has a good reliability.The principle of phase correlation algorithm is based on the characteristics of the Fourier transform.Two image blocks will reflect a linear phase angle in the frequency domain, if they only have offsets.When offsets, x  and y  , exist between image block g and image block f : Convert two sides of formula (1) through Fourier transform, and consider the property of Fourier transform, the resolve can be expressed as: Where, G and F is the Fourier transform matrix of image g and f respectively.
The cross-power spectral function can be acquired by formula (2).
Where, " *" is matrix dot product and G is conjugate function of G .
Inverse Fourier transform can be used in cross-power spectrum above and Dirac function   , xy Where, IFT is the function of inverse fast Fourier transform.If two blocks of image show the same area, peek value of the pulse function can be calculated in the place   , x y   , and other value, around the peak value, will be far less than the peek value and close to zero.

Principle of Peak Calculation
Classical phase correlation matching only can attain pixel precision according to acquire the path and row of the maximum value in matrix of impulse function.This paper presents a sub-pixel phase correlation matching method, which calculates power peak value base on the symmetrical distribution around energy peak.In the algorithm, the order of value around the peak determines formula, so there are two different conditions to calculate the peak point.When the peak value is greater than the back value, peak calculation algorithm sketch is shown in Figure 1: There are two geometric relationships according to the figure above and the formula can be expressed as follows: Formula ( 5) can be deduced and simplified as: Formula above can be simplified as follows: In the above formula， The solution of the formula above can be attained as follows: The solution belong to the interval of 12 x x x  is requested.
Similarly, when points , C x y include 13 yy  , the formula can be expressed as follows: The solution x of the formula ( 9) can be attained according to formula ( 7) and ( 8), which belong to the interval of yy  , the images to be matched are regarded as no offset, and 0 x  .

Matching Window Constraint
The window mentioned to constraint the result is designed as 3 ×1 or 1×3(set it according to row or line), after getting the result of sub-pixel matching, and it can use window constraint to control the value of matching.The ketch of window constraint is shown as Figure 2: When the peak of 3 x is far less than the peak of 1 x , and 1 x is close to the limitation of peak 2 x , it is shown as figure 2, according to geometric knowledge, we know the distance of x deviates 2 x must be less than 0.5, so we can use this constraint to control the value of sub-pixel, using it to act as removed strategy of mismatching points of sub-pixel matching fractional part.

Matching window constraint
Figure 3 the process of phase correlation matching using peak calculation Stereo images L and R respectively refer to image data blocks f and g , and obtaining F and G according to the two- dimensional fast Fourier transform(adding harming window).The cross-power spectrum Q is obtained through formula (3), and  is obtained using the inverse transformation of Fast Fourier Transform(FFT) for Q(adding harming window).After above, integer part of the sub-pixel can be obtained through getting maximum value in matrix of impulse function, and error values are removed by using the window constraint; We can get sub-pixel matching results through the above process.

THE EXPERIMENT AND ANALYSIS
Three experiments are designed to verify effectiveness of algorithm presented.Compared with curved surface fitting phase correlation matching algorithm using simulation data in accuracy and speed, the effectiveness of algorithm presented is verified.These simulation data are obtained by down sampling, thus the absolute matching precision of algorithm presented can be measured.Lastly, the multi-spectral Images of ZY-3 satellite, the first civilian stereo mapping satellite in China, is tested to verify the effectiveness according to detection of satellite attitude jitter.

Experiment I
The high resolution remote sensing images of ZY-3 are selected as raw images, and the spatial resolution of panchromatic image is 2.1 meters.The test images are simulated by down sampling raw images with several rates.The point (X,Y) in the raw image is considered as the beginning point of the simulation image.  1 and table2.According to these results, the stability and accuracy of curved surface fitting phase correlation matching algorithm is worse than the phase correlation matching algorithm using peak calculation, and the accuracy of the presented algorithm can reach 0.1 pixel, and the accuracy of it mostly is better than 0.05 pixel.

Experiment II
In order to verify the speed advantage of the algorithm presented, the same matching window (32×32) and calculation window (3×3) are adopted.
The speed difference of two algorithms mainly displays in calculation of peak value, as curved surface fitting needs the least squares fitting, so its calculation is very complicated.Because calculation of matrix could bring bulk single precision floating-point calculation and it will be greater than one hundred times, the times of multiplication is greater than six hundred times, the times of addition is greater than five hundred times; However, the algorithm presented in this paper only needs simple calculation, the times of multiplication is less than thirty times, and the times of addition is less than twenty five times.The presented algorithm needs a very small amount of calculation than curved surface fitting phase correlation matching algorithm in theory.
Experiment applies 100 × 100 pixel image data block and gradually increasing to 900×900 pixel image data block, a total of nine images(increasing 100 pixels every time), contrasting speed of the algorithm presented with curved surface fitting phase correlation.Computing environment is Intel the second generation of core i7-2820QM@2.3GHzprocessor, single thread operation.Algorithm is written by Visual Studio 2010 platform using Because Fourier transform transformation of phase correlation matching needing large amount of calculation, in order to reflect the relative calculation speed, the computation time of the presented algorithm will subtract the other one.The difference of time is shown as figure 4. As shown in figure 4, the horizontal axis is the image size, and the unit is pixel.The vertical axis is the consuming time, and the unit is millisecond.For the proposed algorithm, as the image block size increases, we can find that the advantage of processing speed is more obvious.It is seen from the above two experiments that the algorithm presented in this paper can acquire better matching accuracy, higher stability and smaller amount of calculation time.

Experiment III
The multi-spectral sensor of ZY-3 consist of four parallel sensors about four bands that include blue, green, red and nearinfrared, and each sensor has three staggered CCD arrays, as shown in figure 4. We apply raw multi-spectral images without any geometric rectification to detect the attitude jitter based on the physical infrastructure above.According to the experiment III, the actual feasibility of algorithm presented will be tested.The test images as shown in figure 5  We can find that the periodic change, caused by satellite platform jitter, is obvious in these parallax images.The regularity, which is reflected in the cross-track parallax, is better than the along-track one, and the main reason may be that along-track direction is affected by more factors.We process the parallax image by mean the each line parallax to get a twodimensional curve, and the two parallax curves can be obtained and showed as figure 7. a) cross-track b)along-track Figure 7. Parallax curves We can find that approximately 0.67Hz (considering the imaging line time 0.8 ms) attitude jitter is detected in multispectral stereo images of ZY-3.The accuracy of matching must reach 0.1 pixel, so the platform jitter can be detected by presented matching method.From the experiment III, the algorithm presented is verified that it could reach the matching accuracy for actual image product.

CONCLUSION
This paper presents a novel phase correlation sub-pixel algorithm based on peak calculation of pulse function matrix, and it takes the geometric relationships into account to infer the corresponding mathematical formula.The amount of calculation is smaller than other methods as without the least squares adjustment, and the theory is simple but complete.According to experimental results, we can conclude the accuracy of algorithm presented in this paper can achieve 0.1 pixel, and meet the need of high-precision matching application.

BB
Figure 1.Peak calculation sketch Where,   2 2 , B x y is peak point in matrix of pulse function, and between its surrounding points   1 1 , A x y and   33 , C x y with Sub-pixel matching algorithm based on phase correlation is implemented by peak value calculating, and its implementation process is shown as the following diagram:

,
Figure 3. Simulated images The images of figure 3 are tested with sub-pixel phase correlation matching based on curved surface fitting and subpixel phase correlation matching based on peak calculation, and these matching results are showed in table 1 and table2.According to these results, the stability and accuracy of curved surface fitting phase correlation matching algorithm is worse than the phase correlation matching algorithm using peak calculation, and the accuracy of the presented algorithm can reach 0.1 pixel, and the accuracy of it mostly is better than 0.05 pixel.

Figure
Figure 4. Time difference Figure 4. Instalment relationshipThe tiny physical distance between the parallel CCD arrays, which are in same scanning column in the direction of flight, will bring parallax matched among them if there is a certain attitude jitter in the satellite platform(Tong X, 2014).
are blue and green band image respectively.a) blue band image b)green band image Figure 5. Experimental images We implement the presented algorithm to process the stereo images, and a lot of tie-points can be obtained pixel by pixel.The dense pixel offset between image a and image b in line and column could form parallax maps in cross-track and along-track direction.The parallax maps are showed as follows: a) cross-track b)along-track Figure 6.Parallax images