A LEAST SQUARE ALGORITHM FOR GEOMETRIC MATCHING OF REMOTE SENSED IMAGES

The aim of geometric matching is to extract the geometric transformation parameters between the corresponding images. That is useful for photogrammetric mapping, deformation detection, and flying platform's posture analyses, etc. It is different compare with ordinary feature based image matching succeed by selecting feature points correctly, the proposed method takes all the pixels within the corresponding images to participate the matching procedure for calculating the geometric parameters by least square criterion. The principle of the algorithm, such as the gray corresponding equation, the information quantity inequation and procedure of least square solution are introduced in detail. Particularly, the wavelet analyses for gray signal and calculating the information quantity by signal to noise ratio. Finally, a series of sequential images obtained by a low-altitude helicopter equipped with a video camera was used to test and verify the validity and reliability of the theory and algorithm in this paper. Two typical results are got according to the relative orientation elements model and parallax grid model. The conclusion is got in comparing APM with ordinary feature point method by the information quantity inequation. 1. INSTRUCTIONS The geometric matching is happened in the case, that the available corresponding images are confirmed taking from a same object surface. The aim of image matching is not to estimate the degree of similarity, but for solving the geometric relationship to get the transformation parameters. It is useful for photogrammetric mapping, deformation detection, flying platform’s posture analyses, registration and mosaic of remote sensing images, etc. Photogrammetric image matching is a typical geometric matching, the research works of this field begin in 1959(Zhizhuo Wang, 1979), and the integrated theory and experience are rich extremely(Konecny and Pape, 1981). One of the best outstanding and widely applied methods may be the matching method based on feature points. However, most pixels of the image are not participated in forming the results, only the selected interesting points and the feature from around pixels are participated the matching, therefore, the mistake matching due to the poor texture of the local image areas are avoided mostly. Many famous results have published, such as Moravec Operator(Moravec, 1977), Harris Operator(Harris and Stephens, 1988), Förstner Operator(Förstner, 1986), Scale Invariant Feature Transform (SIFT) algorithm(Lowe, 1999, Lowe, 2004), High Precision Lease Square Matching(Ackermann, 1984), etc. Although some imperfections still exist, those methods are applied successful, effective and widely spread. It is different compare with above mentioned feature point matching methods; the following proposed method takes all the pixels within the corresponding images to participate the matching procedure for calculation parameters of the geometric models by least square criterion, abbreviation as APM (All Pixels-participated Matching). * Corresponding author 2. PRINCIPLE OF THE ALGORITHM 2.1 Gray Corresponding Equation Assuming the gray signals of the corresponding images to be Gr(x, y) and G[F(x, y)]. Here Gr(x, y) is the gray value at the position (x, y) in the reference image, F(x, y) is the corresponding position of (x, y) in another image. In normal case, noise nr(x, y) and n[F(x, y)] are included in the gray values. The gray corresponding equation can be described as Equation (1): Gr(x, y) − nr(x, y) = G[F(x, y)] − n[F(x, y)] (1) Combining nr and n to be N(x, y), then the Equation (1) can be described as Equation (2): Gr(x, y) − N(x, y) = G[F(x, y)] (2) For least square solution at later procedure, F(x, y) should be decomposed in differential form of the parameters as described in Equation (3): { Fx(x, y) = x0 + ∂x ∂u1 ∆u1+ ∂x ∂u2 ∆u2 + ∂x ∂u3 ∆u3+⋯ Fy(x, y) = y0 + ∂y ∂u1 ∆u1 + ∂y ∂u2 ∆u2 + ∂y ∂u3 ∆u3+⋯ (3) where ui = u0i + ∆ui Therefore, the G[F(x, y)] can be described as Equation (4): G[F(x, y)] = G(x0, y0) + ∂G ∂x ( ∂x ∂u1 ∆u1+ ∂x ∂u2 ∆u2 + ∂x ∂u3 ∆u3+ ⋯)+ ∂G ∂y ( ∂y ∂u1 ∆u1 + ∂y ∂u2 ∆u2 + ∂y ∂u3 ∆u3+⋯) (4) ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume V-2-2020, 2020 XXIV ISPRS Congress (2020 edition) This contribution has been peer-reviewed. The double-blind peer-review was conducted on the basis of the full paper. https://doi.org/10.5194/isprs-annals-V-2-2020-121-2020 | © Authors 2020. CC BY 4.0 License. 121 two kinds of typical examples are as follow: { Fx(x, y) = x0 + ∆bx + x f ∆bz + (f + x2 f )∆φ + xy f ∆ω − y∆κ Fy(x, y) = y0 + ∆by + y f ∆bz + (f + y2 f )∆ω + xy f ∆φ + x∆κ (5) this is the relative orientation equation of photogrammetry. Another, the parallax grid as described in Equation (5): { Fx (x, y) = xi+1−x b yi+1−y b pij + x−xi b yi+1−y b pi+1,j + xi+1−x b y−yi b pi,j+1 + x−xi b y−yi b pi+1,j+1 Fy(x, y) = xi+1−x b yi+1−y b qij + x−xi b yi+1−y b qi+1,j + xi+1−x b y−yi b qi,j+1 + x−xi b y−yi b qi+1,j+1 (6) Here (pij , qij) are the parallaxes in x and y direction in every grid position, as the unknown parameters, b = xi+1 − xi = yi+1 − yi is the interval of the grid. The bilinear interpolation method is applied; of course other suited method is allowed. 2.2 Information Quantity Inequation Many methods can be used for solving the gray corresponding equation (2). From the view of information theory(Chunzhan Tao, 2005, Zongjian Lin, 2006, Zhuyun Fu, 2006), a necessary condition need to be satisfied, that the information quantity (HG) providing from given corresponding images, must larger than the uncertain (HF ) of the unknown parameters of the geometric model. HG = tx ∙ ty ∙ log ( σG 2 σN 2) ≥ ∑ pi s i=1 log 1 pi = HF (7) Here, tx ∙ ty = T is the total amount of pixels participated into image matching, σG 2 and σN 2 are the mean square variance of the signal and noise, σG 2 σN 2 ⁄ be called signal to noise ratio, and log(σG 2 σN 2 ⁄ ) is the mean quantity of information of the pixel. That can be calculated from measurements of sample images, and the geometric model has been solved as the known in this case. e.g Signal of G1(x1, y1) = Signal of G2(x2, y2) (8) Therefore, σG1 2 , σG2 2 and σ(G1−G2) 2 can be calculated by G1(x, y), G2(x, y) and [ G1(x1, y1) − G2(x2, y2) ] respectively. The relationship can be described as Equation (9): SNR = σG 2 σN 2 = σG1 2 +σG2 2

The geometric matching is happened in the case, that the available corresponding images are confirmed taking from a same object surface. The aim of image matching is not to estimate the degree of similarity, but for solving the geometric relationship to get the transformation parameters. It is useful for photogrammetric mapping, deformation detection, flying platform's posture analyses, registration and mosaic of remote sensing images, etc. Photogrammetric image matching is a typical geometric matching, the research works of this field begin in 1959 (Zhizhuo Wang, 1979), and the integrated theory and experience are rich extremely (Konecny and Pape, 1981). One of the best outstanding and widely applied methods may be the matching method based on feature points. However, most pixels of the image are not participated in forming the results, only the selected interesting points and the feature from around pixels are participated the matching, therefore, the mistake matching due to the poor texture of the local image areas are avoided mostly. Many famous results have published, such as Moravec Operator (Moravec, 1977), Harris Operator(Harris and Stephens, 1988), Förstner Operator(Förstner, 1986, Scale Invariant Feature Transform (SIFT) algorithm (Lowe, 1999, Lowe, 2004, High Precision Lease Square Matching (Ackermann, 1984), etc. Although some imperfections still exist, those methods are applied successful, effective and widely spread. It is different compare with above mentioned feature point matching methods; the following proposed method takes all the pixels within the corresponding images to participate the matching procedure for calculation parameters of the geometric models by least square criterion, abbreviation as APM (All Pixels-participated Matching).

Gray Corresponding Equation
Assuming the gray signals of the corresponding images to be ( , ) and [ ( , )]. Here ( , ) is the gray value at the position ( , ) in the reference image, ( , ) is the corresponding position of ( , ) in another image. In normal case, noise ( , ) and [ ( , )] are included in the gray values. The gray corresponding equation can be described as Equation (1): Combining and to be ( , ), then the Equation (1) can be described as Equation (2): For least square solution at later procedure, ( , ) should be decomposed in differential form of the parameters as described in Equation (3): where = 0 + ∆ Therefore, the [ ( , )] can be described as Equation (4) this is the relative orientation equation of photogrammetry.
Another, the parallax grid as described in Equation (5): Here ( , ) are the parallaxes in x and y direction in every grid position, as the unknown parameters, = +1 − = +1 − is the interval of the grid. The bilinear interpolation method is applied; of course other suited method is allowed.

Information Quantity Inequation
Many methods can be used for solving the gray corresponding equation (2). From the view of information theory (Chunzhan Tao, 2005, Zongjian Lin, 2006, Zhuyun Fu, 2006, a necessary condition need to be satisfied, that the information quantity ( ) providing from given corresponding images, must larger than the uncertain ( ) of the unknown parameters of the geometric model.
Here, • = is the total amount of pixels participated into image matching, 2 and 2 are the mean square variance of the signal and noise, 2 2 ⁄ be called signal to noise ratio, and log( 2 2 ⁄ ) is the mean quantity of information of the pixel. That can be calculated from measurements of sample images, and the geometric model has been solved as the known in this case. e.g 1 ( 1 , 1 ) = 2 ( 2 , 2 ) Therefore, 1 2 , 2 2 and ( 1 − 2 ) 2 can be calculated by 1 ( , ), 2 ( , ) and [ 1 ( 1 , 1 ) − 2 ( 2 , 2 ) ] respectively. The relationship can be described as Equation (9): The uncertainty is calculated by the entropy formula. For example, when need only to solve a point within eight neighborhood of the approximate position ( , ), it means = 9, = 0.111, then = 3.176 . Figure. 1 shows a pair of gray signal profile in direction cutting from a corresponding image pair. By processing using wavelet analyses, a series of waveforms in different scales are shown, and the 2 , 2 are calculated listed in Table. 1. It is clear that the SNR at higher resolution levels are lower than at lower resolution levels. Particularly, the detail components 1 , 2 , 3 show negative value of log ( 2 2 ⁄ ), that means, such components have negative quantity of information due to the noise power over the signal power. It may lead mistake matching to produce much trouble for later procedure. In this example, the scale of 4 or 5 may be selected, 5 is better. The multi resolution analysis is based on binary scale wavelet with base function as Equation (10): Contrast with Figure 1 and Table 1, the relationship of multi resolution signals can be described as Equation (11): 0 = 1 + 1 = 2 + 2 + 1 = 3 + 3 + 2 + 1 = 4 + 4 + 3 + 2 + 1 = 5 + 5 + 4 + 3 + 2 + 1 (11) Figure 2. Parallax grid Based on the scale factor = 2 , the interval of parallax grid in equation (6) can be determined, that the interval should larger than . e.g. > , shown as Fig.2. Considering the resolution and scales, the information quantity inequation, equation (7), should be completed as equation (12):

Least Square Solution
After multi resolution analyses, scale factor and are determined, the suitable components of gray signal g ( , ) and g ( 0 , 0 ) is produced instead of ( , ) and ( , ). For every pixel, the gray corresponding equation can be listed as: Here, g = ġ, g = ġ are the gradient of and direction, which can be calculated from gray values. Reform (13) to be: There are • total amount of pixels to form • equations. Using the matrix to described the equations, all the equations can be writing as: Solving by least square criterion as Equation (16):

Pre-processing of Images
For removing the system difference between the corresponding images, pre-processing may need, such as histogram equalization, etc. Otherwise need to add radiation difference parameters into the Equation (2).

Getting the Approximates of Geometric Parameters
Least square matching need approximates of the geometric parameters. That can be got from many way, such as the pyramid structure processing, Radon transform based one-dimensional correlation, etc.

Wavelet Analyses and Low-pass Filtering
It is convenient to use the wavelet software box to make multiresolution analyses. Then based on the signal to noise ratio to select the suitable resolution level of signal, and to determine the interval or detailed degrees of geometric model. In practice, lowpass filtering of gray signals may not be restricted in wavelet; other low-pass filters such as Gauss Low-pass Filter may suitable also.

Geometric Model Design
Beside the above mentioned two kinds of typical geometric parameters, there are many ways to design the geometric models according to requirements of application flexibly. The Fourier transform or Wavelet transform may be very interesting. Take easy for example; one-dimensional Fourier transform can be described as follow equation: Here, , are the unknown parameters.

Experimental Dataset
A series of sequential images taking from a video camera mounted on a helicopter are selected for experiment. Two adjacent images are show at Figure.3.
The parallax grid interval is determined = 100 > , the searching area within 8 × 8 pixels neighbourhood.
The redundant information quantity is useful for enhancing the reliability and for supporting the least square solution to improve the precision of results. For example, when need to refine the results into 0.1 pixel precision, then the should be recalculated as: According to the parallax grid model, 10 * 18 grid parallax values are resulted shown in Figure 4. According to the relative orientation model, 6 parameters are resulted shown in Table.2. Three adjacent images are selected for compare the experimental results. Table 3 shows the errors of 6 parameters.   Table 3. Relative orientation errors (pixel/degree) Using the resulted parameters to make rectification of three images, to make mosaic by every 1/3 part to form a picture shown at Figure 5, that has no any gaps can be found.

Making Combined Wide Angle and Long focal length Aerial Camera
Imaging angle and focal length are the important parameters for aerial camera. When the flying altitude is giving, limited by the weather condition, the long focal length provide high resolution imaging and the wide angle provide high efficiency of photographic production. Therefore, a combined camera with wide angle 84° and long focal length 200mm, named CKAC300 have been developed, the structure diagram and flying platform of CKAC300 shown in Figure 6. The APM algorithm was used to transform the 36 central projective images, from Canon 5DS, to be a single central projective image. Figure 6. Structure diagram and flying platform of CKAC300

Developing the Video Monitor for Geological Deformation of the Slopes
Based on the high precision of image matching results from APM, a video monitor for geological deformation have been developed. A series of video images was taking from a slope surface by a fixed digital camera, and the design structure and field installation status shown in Figure 7. Then, every adjacent image pair to be make matching by APM, to get the deformation value with 0.1~1.0mm precision at a distance of 100 meters

CONCLUSIONS
It is complex to compare the feature points method and APM method. But from information theory view, that the similarities and differences of the two methods can be described by information quantity inequation, the formula (7) or (12). The first method selected the points with larger SNR from the whole image; the APM selected signal components of resolutionfrequency with larger SNR from all pixels of the whole image. The first method decomposes the geometric model to be every isolated point for reducing the uncertainty of geometric model; APM selected a finite number of parameters representation the geometric model through multi-resolution analyses for reducing the uncertainty. The common point of the two method is the necessary condition of information quantity inequation, that the information quantity providing from corresponding images must larger than the uncertainty of the unknown parameters of the geometric model. Only the principle of APM method are expound above.
Combining with the practical applications, there are many detailed technical problems need to study in the future research.