QUANTITATIVE ASSESSMENT OF THE PROJECTION TRAJECTORY-BASED EPIPOLARITY MODEL AND EPIPOLAR IMAGE RESAMPLING FOR LINEAR-ARRAY SATELLITE IMAGES

Epipolar geometry rectification is one of the critical issues in photogrammetry, which is a strong corresponding searching constraint in dense image matching for 3D reconstruction. In this paper, the properties of the projection trajectory-based epipolarity model are analyzed quantitatively, and the approximate straight line and parallelism property of the epipolar curve are discussed comprehensively using the linear pushbroom satellite images, i.e. IKONOS, GeoEye images. Based on the analysis of the epipolar line properties, a practical method for epipolar resampling developed. In this method, the pixelwise relationship is established between the original and the epipolar images. The experiments on TH-1 images show that quasi rigorous epipolar images can be resampled using our proposed method for both along-track images. After epipolar geometry rectification, the vertical parallaxes at checkpoints can achieve sub-pixel level accuracy, thus demonstrating the correctness and applicability of the proposed method.


INTRODUCTION
Linear-array pushbroom imaging sensors have been the leading payload of current high-resolution optical satellites. As full used sensors, it has fueled applications in a variety of domains, including geolocation, digital surface model (DSM), 3D reconstruction and mapping by delivering satellite images and geolocation models (i.e. a rigorous model or a Rational Polynomial Coefficients (RPC) model as a polynomial approximation from different satellite image vendors) (Y. Zhang et al., 2004).
Since the geometry of the linear pushbroom imaging sensors model is more complicated when comparing with the collinearity equations of the perspective cameras. The methods of the traditional frame perspective imaging model are no longer universally applicable for the satellite images. Linear-array pushbroom images are characterized "line of central projection" dynamic imaging and the exterior orientation elements vary for each scanning line. There is no rigorous definition of the epipolar line as the frame perspective image. Moreover, the coverage area of high-resolution satellite images is much larger than aerial or street-view images, which results in a vast number of pixels (billion-level pixels). This is how the considerable gap between the photogrammetry and computer vision was created over the past few decades, which should be bridged to benefit each other in scientific development. Previous studies analyzed less on the geometric property of the epipolarity model or the central projection models in the aerial images are adopted directly. Research on satellite image epipolarity model and epipolar geometry rectification method has been one of the hot topics in the field of aerospace photogrammetry and remote sensing.
In the early 1980s, Z. Zhang and Zhou (1989) proposed a method for generating an approximate epipolar line based on polynomial * Corresponding author fitting of the corresponding point coordinates for the SPOT across-track stereo images. Kim (2000) analyzed the projection trajectory (PT) method based on a rigid model. It has been clear about the epipolar curve shape; at the same time, it concluded the approximate linear characteristic of the epipolar curve within a certain range, which laid a foundation for the subsequent application. Morgan et al. (2006) and Habib et al. (2005) proposed a method for generating an approximate epipolar line based on a parallel projection model, according to the small imaging field angle characteristics of high-resolution satellite imagery. The sensor model based on parallel projection was adopted for satellite image, then the epipolar line model in the form of a straight line was proposed based on at least four pairs of corresponding points. Wang et al. (2011) obtained the direction of the epipolar line on the reference plane by the projection trajectory method and projected along the epipolar line direction to get an approximate epipolar line.
By adopting this method, the rigorous coordinates corresponding relationship between the original image pixel and the epipolar image pixel can be easily derived from the orientation data of satellite stereo imagery. Based on the Rational Polynomial Functions (RPF), G. Zhang et al. (2010) put forward a practical way of generating an epipolar image by the projection trajectory method and rebuilding its geometric model. In general, some deficiencies still exist in the above methods: (1) most of the methods are a kind of approximate method, and only the final resampled epipolar images are evaluated, lacking quantitative analysis in the whole process, especially on the epipolar lines; (2) present epipolar resampling methods are complicated for the epipolar geometry rectification, Level-2 images and DSM generation.
The remainder of this work is organized as, 1) In Section 2 the principle and foundation of the projection trajectory-based epipolarity model are introduced, and the properties of epipolar curve in the model is analyzed to address the problems mentioned above; 2) In Section 3, a novel epipolar rectification method is proposed based on the model leveraging the approximate conjugate property of epipolar line, followed by accuracy assessment and discussion on TH-1 dataset; 3) Section 4 concludes the work.

The Projection Trajectory-based Epipolarity Model
Given the geometry of the frame perspective images can no longer be dealing with the linear pushbroom satellite images because differences existing in the imaging models, many researchers in aerospace photogrammetry community have made unremitting efforts and attempts to build the epipolarity model of linear pushbroom images. The most rigorous epipolarity model of the linear pushbroom image is the projection trajectory based on RPF.
As shown in figure 1, if all the points of the optical ray from the ground point through the left imaging centre 1 ( 1 , 1 , 1 ) to the image point on the left image, are projected to the right image 2 through the right imaging center 2 ( 2 , 2 , 2 ) , the projection trajectory of these points will form a curve on the right image. This curve is termed an epipolar line of the point . If is the corresponding point of , it is easy to know that always locates on the epipolar curve . This is the geometric definition of the projection trajectory-based epipolar line for the linear pushbroom satellite images. It is easy to know that for the centre perspective camera, the epipolar line based on the projection trajectory method is linear. The results conclude that the epipolarity model based on the projection trajectory method is completely consistent with the traditional model. Furthermore, the epipolarity model based on the projection trajectory method is more general. The existing studies have shown that the epipolar curve based on the projection trajectory method has the following properties: (1) The epipolar lines are hyperbola curves for the linear pushbroom satellite stereo images and can be handled as approximately straight lines on a small scale.
(2) The corresponding points of point in an epipolar line and its adjacent points in a certain range on an image are all located in the epipolar line of the point . This conclusion is established at the local scale. According to this conclusion and conclusion (1), the search process of stereo image pairs matching can be constrained to one-dimension from two-dimension.
(3) The corresponding epipolar line pair for the corresponding points exists. If two points are the corresponding points, then the corresponding epipolar lines are pairwise. Furthermore, all the corresponding points on the two epipolar lines are pairwise. The conclusion can be deduced by the conclusion (2) and the conjugate property of stereo image pairs.

Analysis of the PT-based Epipolarity Model
Previous analysis indicates that, on the local scale, the properties of the linear pushbroom satellite image epipolar line model are similar to the epipolar line of the traditional frame perspective image. So how about the linear approximation property in the image range? It needs further quantitative analysis to accurately grasp the epipolar line application properties and provide a theoretical basis for subsequent epipolar geometry rectification. Here two typical linear pushbroom commercial satellite images of IKONOS and GeoEye are selected as the experimental datasets for approximation linear property analysis of the epipolar curve.

Experimental dataset:
Experimental dataset 1 consists of two IKONOS images covering a part area of Beijing, China, as is shown in Figure 2, which was after radiometric correction on the original image without any other geometry processing. The size of the image is 5057×8063 pixels, and the orientation parameters files are provided in standard RPC format.
Experimental dataset 2 is two GeoEye images, as is shown in Figure 3. The size of the image is 26928×15668 pixels, and the corresponding orientation parameters files are provided in standard RPC format.

Left image
Right image In the analysis part on the PT-based epipolarity model of this work, a number of evenly distributed points are selected from the left image; indeed, left images are divided into 10×10 blocks, the centre point of each image block is selected as the test point, 100 test points are selected totally. The point distribution on the whole image is about ten lines multiply ten columns. According to image point coordinates of each test point on the left image and the RPC parameters of this image, straight space line through the projection centre and the image point is determined by the projection trajectory method. Then using the RPC parameters of the right image, the points on this space straight line are projected to the right image one by one. The projection trajectory of these points in the right image is the corresponding epipolar line. For the epipolar curve, an approximate line which is closest to this curve is determined by the least square method. Then the maximum distance between all points on the curve and the fitted straight line is computed. For the 100 test points, the above steps are executed gradually. With these steps, the curvature of the 100 lines, corresponding to these 100 test points, and the maximum distance between each curve and the fitted straight line are reported and analyzed.

Results and analysis on the IKONOS dataset：Table
1 describe the slopes of the straight lines, fitted from the epipolar curves, corresponding to the 100 evenly distributed points in the IKONOS dataset. Table 2 describes the maximum distance between all the 100 points on the curves and the straight lines fitted from these curves, which describes the linear fitting characteristic of the epipolar curve. It should be noted that the columns stand for the order of the 100 points in coordinate direction, while the rows stand for the order of the 100 points in y coordinate direction. It is of importance to analyze the conjugate property of the PT-based epipolarity model. Thus, in order to verify the conjugate characteristic of the epipolar line, the same experimental method is adopted from the right image to the left image as shown in Table 3 and Table 4.  Table 2 reports that considering 100 points distributed evenly on the image, the fitting error is very little, and the maximum value is 0.15623 pixels.     experimental method is adopted from the right image to the left image. The results are shown in Table 7 and Table 8. Table 5 reports that for GeoEye satellite image, the maximum orientation value of the epipolar line on the right image is 75.16487 degrees, and the minimum value is 74.14043 degrees, with a variation range being about 0.015 degrees. If the average slope (75.157 degrees) is regarded as the arranged orientation of the epipolar image, the maximal error exits at the edge of the image, which is about 1.96 pixels. Table 5. Slopes of the fitting straight lines on the right image of the GeoEye dataset (degree) Table 6. Maximum distances of the fitting straight lines on the right image of the GeoEye dataset (pixel) Table 6 reports that considering 100 points distributed evenly on the right image of the GeoEye dataset, the fitting error is very little, and the maximum value is 0.039 pixel.    Comparing with the results on the IKONOS dataset, the epipolar line slope variation range of GeoEye image is larger. That is to say, the parallelism property of the epipolar line is weak. However, the approximate linear property of all epipolar lines is similar to the IKONOS dataset. For the IKONOS image, when the precision range is in the sub-pixel level, epipolar lines on the image can be considered parallel to each other. However, for the GeoEye dataset, epipolar lines are not parallel. This is partially due to the size of the experimental dataset changes greatly, while the size of the IKONOS image used in the experiment is only 5057×8063, and GeoEye is 26928×15668. Experiment results on the IKONOS dataset report that in the whole image range, the PT-based epipolarity model has good approximate linear property and approximate parallel relationship between each epipolar line. Based on the approximate parallel property, the rectification of the epipolar geometry becomes intuitive. Rotated directly by the epipolar line slope angle, the epipolar image can be resampled with very little vertical parallaxes.
It can be concluded from the above analysis: the projection trajectory-based epipolarity model has distinct different parallel properties on linear pushbroom satellite images with different sizes. Nevertheless, in a local range in the image, the approximate linear property is satisfied for the sub-pixel level requirement. Consequently, this projection trajectory-based epipolar curve can be handled as a straight line in the local-scale image.

EPIPOLAR IMAGE RESAMPLING
It is of importance to generate the epipolar images, which is one of the critical constraints for stereo epipolar images matching in the practical work of surveying and mapping. Traditional methods of generating epipolar images can be basically divided into two categories: one is based on geometric correction of digital images; The second is based on the coplanar condition.
The method based on digital image geometric correction is essentially a digital correction, projecting the epipolar line of the oblique image to the ortho-ready image plane to obtain corresponding epipolar lines of ortho-ready images. The coplanar condition-based method can directly determine the direction of the epipolar line on the oblique image, then determine the location of the epipolar line on the oblique image according to the coordinates of one epipolar line's start point and the direction of epipolar line.
Referring to the epipolar image generation method based on coplane condition, we proposed a method to rectify the epipolar geometry utilizing the properties of the projection trajectorybased epipolarity model as discussed in the aforementioned sections. Referring to the method of traditional frame perspective remote sensing images, which samples the epipolar image from the oblique image directly, the epipolar image of the whole area is acquired through processing line by line. According to the epipolar model assessed above, when the point in the left image is known, the optical ray determined by left image point and the projection centre S 1 of the left image is projected on the right image and formed a curve l; then, a straight line is fitted to approximate the curve, and the straight line is decided as an epipolar line on the right image. In the same way, the optical ray determined by right image point q and the projection centre S 2 of the right image is projected on the left image to form a curve ′ . A straight line is fitted to approximate the curve, and the line is epipolar line ′ corresponding to l . Finally, the equation parameters of l and ′ are used to resample epipolar lines in the original image directly.
In addition, for the generated epipolar images, the epipolar image coordinates can be converted from the original image coordinates by the parameters of the two straight lines. It is critical information that should be recorded and thus can be applied for computing digital surface models on the translation stage from the original image pairs.

The Experimental Dataset
A Tianhui-1 (TH-1) satellite dataset is selected in this section, as shown in Figure 4. This dataset contains three scenes of images and the corresponding RPC files. TH-1 is China's first transmissive optical stereo mapping satellite. The satellite was launched on August 24, 2010, which carries on three kinds of loads including three linear array stereo mapping cameras, a high-resolution camera and four bands multi-spectral camera. The resolution of three linear-array panchromatic images is better than 5 meters (including forward, nadir and back). The resolution of the multi-spectral image is better than 10 meters. The resolution of the high-resolution image is better than 2 meters.
In this experiment, one scene of TH-1 three linear array image in Baoji, Shaanxi region in July 2011 is chosen as is shown in Table 9 below. The cover area of this dataset is 3600 square kilometres (60 km×60 km). There is no cloud coverage in the images. The elevation difference is moderate, which has both mountain and urban areas.
a. Forward b. Nadir c.Backward Figure 4. Overview of the TH-1 Satellite Images

Methodology for Epipolar Image Resampling
Taking an example of the TH-1 satellite stereo image pair generation with forward and nadir, the detailed steps are below: (1) Starting from the first column in the forward image, the middle image point (0, 2 ⁄ ) in the column direction is selected. Based on the PT-based epipolarity model, the line determined by point and the projection centre S of the left image is projected on the right image. Considering the computation, in this paper only two object space points 1 and 2 of this line on elevation ℎ and ℎ are projected on the right image to determine the two image points 1 and 2 .
(2) Determine whether the point 1 and 2 are in the range of the right image. If not, return to (1) and continue to start from the next column. According to these two points, a straight line ′ is determined and the line ′ + ′ + = 0 is recorded. (3) Determine the midpoint q between 1 and 2 on the right image. According to the PT-based epipolarity model, the optical ray determined by the point q and the projection centre ′ of the right image is projected on the left image.
Considering the computation simpleness, in this work only two object space points 3 and 4 of this optical ray on elevation ℎ and ℎ are projected on the left image to determine the two image points 3 and 4 . (4) Determine whether the point 3 and 4 are in the range of the left image or not. If not, return to (3) and continue to start from the next column. According to these two points, a straight line is determined and recorded: Using Where and y are the pixel coordinates on the original image, which are the number of columns and rows respectively; And ′ and ′ are the pixel coordinates on the epipolar image, which are the number of columns and rows respectively; height represents the rows number of the original image; is the number of epipolar equations; is the number of straight lines; ， ， indicate the parameters of the line equations.

Result Assessment
The epipolar images for TH-1satellite images are generated by the method proposed in Section 3.1, including the forward image, the nadir image, and the backward image. As shown in Figure 5.
a. Forward b. Nadir c.Backward Figure 5. Computed Epipolar Images In order to verify the vertical parallax of the epipolar image, firstly the stereo mapping module of VirtuoZo was adopted to select 31 pairs of evenly distributed corresponding in the epipolar stereo image pair. Then the coordinates and its vertical parallax of these corresponding image points in the epipolar stereo image pair are computed, and the results are shown in Table 9.  Table 9. The Vertical Parallax evaluation 1 is the column value of the corresponding point on the forward epipolar image, and 1 is the row value; 2 is the column value of the corresponding point on the back epipolar image, while 2 is the row value. The vertical parallax (abs ( 2 -1 ) ) of the epipolar image is computed to describe the precise of epipolar images, especially the left and right edge points. Table 1 indicates that the vertical parallax of corresponding image points involved in the inspection can achieve accuracy at sub-pixel and root mean square error can achieve accuracy up to 0.590 pixels. It shows that this method can generate a quasi rigorous satellite epipolar image.

CONCLUSION
The contribution of this work is two folds. On the one hand, we explored the application of the rational function model in the projection trajectory-based epipolarity model for linear pushbroom satellite images. The projection trajectory-based epipolarity model is analyzed quantitatively using IKONOS and GeoEye datasets. The approximate linear and parallelism property of the epipolar curve in the image is analyzed.
On the other hand, for the good approximate linear property of epipolar images, referring to the epipolar line sampling method of traditional frame perspective images, a method of epipolar image generation is proposed, which is similar to the method based on coplanar conditions. Based on the projection trajectorybased epipolarity model, treating the epipolar line as the approximate straight line, this method established a corresponding rigid relationship of coordinates between the original image and epipolar image by the orientation parameters of stereo images, and then carried out the epipolar image resample directly on the oblique image, and finally generated quasi rigorous epipolar image pairs with the vertical parallaxes in sub-pixel level.
The proposed method is practical and straightforward, which is performed on the TH-1 satellites images in this work. At present, the proposed method has been successfully applied to the followup stereo mapping and automatic DSM generation software for the TH-1 satellite images.