FOREWORD A NOVEL SUB-PIXEL MATCHING ALGORITHM BASED ON PHASE CORRELATION USING PEAK CALCULATION Junfeng Xie a , Fan Mo a b , Chao Yang a , Pin Li a d , Shiqiang Tian a c , a Satellite Surveying and Mapping Application Center of China, NO.1 Baisheng Village, Beijing - xiejf@sasmac.cn,yangc@sasmac.cn b Information Engineering University, No.62 Kexue Road, Zhengzhou - surveymofan@163.com c Chang'an University, Middle-section of Nan'er Huan Road, Xi'an - 835301221@qq.com d Liaoning Project Technology University, People Street, Fuxin - 1076760488@qq.com Commission I, WG I/3 KEY WORDS: Image Matching, Phase Correlation, Peak Calculation, Window Constraint, Correlation Coefficient ABSTRACT: 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. 1. 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 cross- correlation 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. 2. METHOD 2.1 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 :    , ,g x y f x x y y   （1） Convert two sides of formula (1) through Fourier transform, and consider the property of Fourier transform, the resolve can be expressed as:      , , i u x v y G u v F u v e      （2） 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).             , , , , , i u x v y F u v G u v Q u v e F u v G u v         （3） The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLI-B1, 2016 XXIII ISPRS Congress, 12–19 July 2016, Prague, Czech Republic This contribution has been peer-reviewed. doi:10.5194/isprsarchives-XLI-B1-253-2016 253 mailto:835301221@qq.com 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  ,x y   can be obtained in place  ,x y  :        , , i u x v yx y IFT Q u v IFT e        （4） 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. 2.2 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:   1 x 2x 3x X Y 3 y 2 y 1 y x P 1 l 2 l A B C 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  3 3,C x y with 1 3 y y . Make a straight line 1 l from C to B and make a vertical line across B , and we can attain the angle between the two line is  , according to the principle of energy accord with symmetrical distribution in matrix of pulse function , there must be a line 2 l symmetrical with line 1 l about vertical line that across the peek point P . Therefore, make a straight line 2l from A with an angle 90  , and make 2l and 1l meet at point P , the difference 2 x x between abscissa x of P and abscissa 2 x of B is sub-pixel offsets. There are two geometric relationships according to the figure above and the formula can be expressed as follows: 3 1 3 1 3 2 3 3 3 2 x x x x y y y y y y y y x x x x             （5） Formula (5) can be deduced and simplified as:       2 3 2 2 1 3 1 2 3 2 10 9 3 12 9 0 y y x y y y x y y y         （6） Formula above can be simplified as follows: 2 0ax bx c   （7） In the above formula，  3 22a y y  , 2 1 310 9b y y y   , 1 2 3 3 12 9c y y y   . The solution of the formula above can be attained as follows: 2 4 2 b b ac x a     （8） The solution belong to the interval of 1 2 x x x  is requested. Similarly, when points  1 1,A x y and  3 3,C x y include 1 3 y y , the formula can be expressed as follows:       2 1 2 2 1 3 1 2 3 2 6 7 5 4 0 y y x y y y x y y y         （9） The solution x of the formula (9) can be attained according to formula (7) and (8), which belong to the interval of 2 3 x x x  are requested. In addition, when points  1 1,A x y and  3 3,C x y include 1 3 y y , the images to be matched are regarded as no offset, and 0x  . 2.3 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: 1 x 2 x 3 x Figure 2. Window constraint sketch When the peak of 3x is far less than the peak of 1x , 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 2x 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. The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLI-B1, 2016 XXIII ISPRS Congress, 12–19 July 2016, Prague, Czech Republic This contribution has been peer-reviewed. doi:10.5194/isprsarchives-XLI-B1-253-2016 254 http://dict.youdao.com/w/abscissa/ http://dict.youdao.com/w/abscissa/ Sub-pixel matching algorithm based on phase correlation is implemented by peak value calculating, and its implementation process is shown as the following diagram: Image Image Image block Image block frequency frequency cross-power spectrum impulse function Peak integer value the integer part Peak value calculation The result of sub- pixel matching L R f g F G Q  Fouriertransform fourier inverse transformation 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. 3. 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. 3.1 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. The “A” image with 1000 1000m m pixels is cropped. Similarly, The “B” image with 1000 1000m m pixels cropped from point  1, 1X Y  in the same raw image ; A and B images are separately down sampled to 1000×1000 pixels, then we get image a and image b. Theoretically, the deviant between image a and image b is 1 1 , m m       , and m is the sampling rate. We choose four frame simulated images with different ground features, as are shown in figure 3: m a image b image 3 5 10 20 Figure 3. Simulated images The images of figure 3 are tested with sub-pixel phase correlation matching based on curved surface fitting and sub- pixel 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. Table 1 the result of curved surface fitting phase correlation matching algorithm m Matching error (%) (0, 0.05) pixel (0.05, 01) pixel (0.1, 1) pixel 3 43.2 43.5 13.3 5 20.7 56.6 22.7 10 59.2 36.4 4.4 20 79.7 19.6 0.7 Table2 the algorithm of phase correlation matching adapted by peak calculation The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLI-B1, 2016 XXIII ISPRS Congress, 12–19 July 2016, Prague, Czech Republic This contribution has been peer-reviewed. doi:10.5194/isprsarchives-XLI-B1-253-2016 255 file:///D:/软件/Dict/6.3.67.7016/resultui/frame/javascript:void(0); file:///D:/软件/Dict/6.3.67.7016/resultui/frame/javascript:void(0); javascript:void(0); file:///D:/软件/Dict/6.3.67.7016/resultui/frame/javascript:void(0); file:///D:/软件/Dict/6.3.67.7016/resultui/frame/javascript:void(0); javascript:void(0); javascript:void(0); m Matching error (%) (0, 0.05) pixel (0.05, 01) pixel (0.1, 1) pixel 3 79.2 16.3 4.5 5 66.9 23.6 9.5 10 72.8 21.8 5.4 20 75.6 20.9 3.5 3.2 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 block 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.3GHz processor, single thread operation. Algorithm is written by Visual Studio 2010 platform using MFC. 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. Figure 4. Time difference 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. 3.3 Experiment III The multi-spectral sensor of ZY-3 consist of four parallel sensors about four bands that include blue, green, red and near- infrared, and each sensor has three staggered CCD arrays, as shown in figure 4. B1 B4 B3 B2    152 pixel 128 pixel 128 pixel Figure 4. Instalment relationship The 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). 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 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 The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLI-B1, 2016 XXIII ISPRS Congress, 12–19 July 2016, Prague, Czech Republic This contribution has been peer-reviewed. doi:10.5194/isprsarchives-XLI-B1-253-2016 256 javascript:void(0); javascript:void(0); javascript:void(0); javascript:void(0); 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 two- dimensional 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 multi- spectral 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 the 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. ACKNOWLEDGEMENTS Acknowledgements of support for the public welfare special surveying and mapping industry (NO.201412001, NO.201512012) and natural science research fund (NO.41301525, 41571440), and the scientific research plan for academic and technical youth leaders of State Bureau of Surveying and Mapping (NO. 201607), youth science and technology of surveying and mapping project (NO. 1461501900202), and the Major Projects of High Resolution Earth Observation System. REFERENCES ARMIN GRUEN. Development and Status of Image Matching in Photogrammetry[J]. The Photogrammetry Record, 2012, 27(137): 36-57. Dazhao Fan, Erhua Shen, Lu Li et al. Small Baseline Stereo Matching Method Based on Phase Correlation[J]. Journal of Geomatics Science and Technology, 2013, 30(2):154-157. Harold S Stone, Michael T, Ee-Chien Chang, et al. A Fast Direct Fourier-Based Algorithm for Subpixel Registration of Image[J]. IEEE Trans. GEOSCIENCE AND REMOTE SENSING, 2001, 39(10):2235-2243. K. Ito, H. Nakajima, K. Kobayashi et al. A Fingerprint Matching Algorithm Using Phase-Only Correlation[J]. IEICE Trans. Fundam. Electron. Commun. Comput. Sci., 2004, 87(3):682-691. Kenji TAKITA, Tatsuo HIGUCHI. High-Accuracy Subpixel Image Registration Based on Phase-Only Correlayion[J]. IEICE TRANS. FUNDAMENTALS, 2003, E86-A(8):1925-1934. Kuglin, C.D., Hines, D.C. The Phase Correlation Image Alignment Method[J]. Proc. IEEE. Int’ Conf. on Cybernetics Society, 1975, 163-165. S. Leprince, S. Barbot, F. Ayoub, et al. Automatic and Precise Orthorectification, Coregistration, and Subpixel Correlation of Satellite Image, Application to Ground Deformation Measurements[J]. IEEE Trans. Geosci. Remote Sens., 2007, 45(6):1529-1558. T. Heid, A. Kääb. Evaluation of Existing Image Matching Methods for Deriving Glacier Surface Displacement Globally from Optical Satellite Imagery[J]. Remote Sens. Environ., 2012, 118:339-355. W. S. Hoge. A Subpixel Identification Extension to the Phase Correlation Method [MRI Application][J]. IEEE Trans. Med. Image., 2003, 22(2):277-280. Xiaohua Tong, Zhen Ye, Yusheng Xu, et al. A Novel Subpixel Phase Correlation Method Using Singular Value Decomposition and Unified Random Sample Consensus[J]. IEEE TRANSCATION ON GEOSCIENCE AND REMOTE SENSING, 2015, 53(8):4143-4156. Tong X, Ye Z, Xu Y, et al. Framework of Jitter Detection and Compensation for High Resolution Satellites[J]. Remote Sensing, 2014, 6(5):3944-3964. The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLI-B1, 2016 XXIII ISPRS Congress, 12–19 July 2016, Prague, Czech Republic This contribution has been peer-reviewed. doi:10.5194/isprsarchives-XLI-B1-253-2016 257