RELATIVE ORIENTATION AND MODIFIED PIECEWISE EPIPOLAR RESAMPLING FOR HIGH RESOLUTION SATELLITE IMAGES

High resolution, optical satellite sensors are boosted to a new era in the last few years, because satellite stereo images at half meter or even 30cm resolution are available. Nowadays, high resolution satellite image data have been commonly used for Digital Surface Model (DSM) generation and 3D reconstruction. It is common that the Rational Polynomial Coefficients (RPCs) provided by the vendors have rough precision and there is no ground control information available to refine the RPCs. Therefore, we present two relative orientation methods by using corresponding image points only: the first method will use quasi ground control information, which is generated from the corresponding points and rough RPCs, for the bias-compensation model; the second method will estimate the relative pointing errors on the matching image and remove this error by an affine model. Both methods do not need ground control information and are applied for the entire image. To get very dense point clouds, the Semi-Global Matching (SGM) method is an efficient tool. However, before accomplishing the matching process the epipolar constraints are required. In most conditions, satellite images have very large dimensions, contrary to the epipolar geometry generation and image resampling, which is usually carried out in small tiles. This paper also presents a modified piecewise epipolar resampling method for the entire image without tiling. The quality of the proposed relative orientation and epipolar resampling method are evaluated, and finally sub-pixel accuracy has been achieved in our work.


INTRODUCTION
Accompanied by the emergence of dense image matching method like Semi-Global Matching (SGM) method (Hirschmuller, 2008), more and more researchers shift their interests to use photogrammetric method for Digital Surface Model (DSM) generation and 3D reconstruction.Since the Earth observing satellite products have grown to a new era, many commercial satellite vendors provide very high resolution stereo images at sub-meter resolution.For instance, the latest DigitalGlobe's satellite, Worldview-4 is able to provide panchromatic imagery of 31cm.The satellite stereo images are showing the attractive prospects in dense image matching area.Thus since recent years, a lot of researchers have done experiments on dense image matching, DSM generation and 3D reconstruction with high resolution satellite imagery (d'Angelo and Reinartz, 2011;Wohlfeil et al., 2012;Franchis et al., 2014;Gong and Fritsch, 2016;Ghuffar, 2016;Shean et al., 2016;Rita et al., 2017).Among various dense image matching methods, SGM method is a popular and efficient way.In order to apply SGM method to satellite stereo imagery, image orientation and rectification are essential pre-procedures.
Different from standard photogrametry, commercial satellite imagery vendors usually deliver Rational Polynomial Coefficients (RPCs) instead of ex-or interior elements to the customers.RPCs represent a generalized and direct relationship between image and object coordinates.This relationship is pure mathematic and without any physical meanings.Former researches have verified, that the RPC model can maintain the accuracy as the same as the rigorous sensor model (Grodecki and Dial, 2001;Hanley and Fraser, 2001).However, RPCs from vendors are not always reliable, so that it needs to be refined by some orientation methods.Grodecki and Dial (2003) proposed the bias-compensated RPCs bundle block adjustment.Through this method, an additional affine model is utilized to compensate the bias in image space for RPC refinement.It is an absolute orientation procedure, so that a proper number of ground control points (GCPs) is required.Sometimes, the GCPs may not available, so Franchis et al. (2014) proposed a relative orientation method without any GCPs.They point out, that if the RPCs are erroneous, for a point x, its corresponding epipolar line will not pass the corresponding point x'.This shift between the corresponding point and its epipolar line is called the "relative pointing error".When the image is small (e.g.1000*1000 pixels), several corresponding point pairs are given and their relative pointing errors are measured as simple translations.The relative pointing error of the small image is removed by the median of the translations.But for large regions, the local pointing correction is not valid.As mentioned in their work, big range images have to be divided into several small tiles.And for each tile, the local pointing corrections are calculated.Then all those local translations are used to estimate an affine transformation that corrects the global relative point error.They only employed the local pointing error correction because they conduct the image rectification in small tiles.Ghuffar (2016) also verified this local relative orientation method and its related tile-wise epipolar resampling.
Epipolar resampling procedure is important because the corresponding points locate on the same row in the generated epipolar images, which reduces the search range of matching form 2D to 1D space.This character improves the efficiency of dense image matching significantly.Unlike the traditional frame images, high resolution satellite images are hard to generate the epipolar geometry because of the changing perspective center and the attitude.As Kim (2000) explained in his work, the epipolar curves of satellite pushbroom sensor are more like hyperbola curves than straight lines, and the epipolar pairs only exist locally.It is well-known, that the RPCs provide a direct relationship between object and image space.Thus, the projection-trajectory method supported by the RPCs is commonly used to find the corresponding epipolar pairs (Wang, et al., 2010).
Several different solutions are researched to solve the problem of establishing satellite image's epipolar geometry and the image resampling.For example, Wang et al. (2011) defined a Project Reference Plane (PRP) in local vertical coordinate system to build the epipolar geometry.In this method, selected point in master image is projected to two height levels upper and lower to the PRP by the RPCs' projection-trajectory method.Then the back projective points of these two height levels in slave image are calculated, and they are projected to the height level, where PRP locates.These two projected points on PRP are used to define the approximate direction of the epipolar curve.An affine model is estimated by the epipolar curves, and it is applied for transformation from the original images to the resampled epipolar images on PRP.Based on this method, Koh et al. (2016) proposed a piecewise epipolar resampling method.Their method sets several control points on the PRP first, and approximates the piecewise epipolar curve between two control points.Then they use a fifth order polynomial function to implement the epipolar transformation.Different from those two methods introduced before, Oh (2011) conduct the epipolar resampling in image space instead of object space.According to his work, the center points of the images are applied to obtain the orthogonal line to the track.Then for each image, the start points on this orthogonal line are selected with proper interval (1000 pixels).Epipolar curves are expanded from the start point.The curves are aprroximated as segments, and the length is defined by proper height range.This height range is equal to or larger than the actual terrain elevation range.When all the epipolar segments on left and right images are derived, each epipolar curve pair are assigned to a constant row for y-parallaxes removal.The epipolar image resampling is done by high order polynomial transformation or interpolation methods.
Based on the image space resampling method, this paper proposes a modified piecewise epipolar geometry and resampling strategy for an entire image without tiling.Since the RPCs' accuracy has effects on the epipolar pair's quality, we also present two global relative orientation methods for large range satellite images.The methodology of the two relative orientation methods is described in the next section.In section 3 we introduce the method of the modified piecewise epipolar resampling.The experiments are carried on QuickBird and WorldView-2 satellite imagery and their results are presented in section 4. Finally, some conclusions are drawn in section 5.

GLOBAL RELATIVE ORIENTATION
The RPCs provided by the data vendors may cause a bias between the calculated coordinates and the true location.This bias presents the systematic error, and it varies from sub-pixel to tens of pixels for different data.The qualities of the epipolar image generation and subsequent procedures are all depended on the accuracy of the RPCs.Therefore, for most satellite stereo imagery, the orientation procedure is essential, and it is almost equal to the refinement of the RPCs.The additional bias-compensation bundle block adjustment is a convenient and accurate method to refine the RPCs and remove the bias (Grodecki and Dial, 2003;Fraser and Hanley, 2005).As we know, the limitation of this method is the requirement of GCPs.In this section, two different relative orientation methods are presented.Both methods are free of GCPs and are applied to an entire image without tiling.

Relative Bias-compensated Model
It is well-known, that the 80-parameter RPCs build the bridge between object and image coordinates as a ratio of two thirdorder polynomial functions.Grodecki and Dial (2003) suggest using an additional affine model to do the bundle block adjustment when the image size is large.The additional biascompensated model is defined in a comprehensive form as follow: The bias-compensated bundle block adjustment is an absolute orientation procedure and, practically, the minimum number of GCPs should be 4 to 6 (Fraser and Hanley, 2005).
For the situation, that the ground control information is not available, we propose a relative additional bias-compensated model.Given a set of correspondences in both images, we first calculate the corresponding object coordinates by forward intersection.It should be noticed that these calculated object coordinates have shifts to the true location in the object space, because the RPCs have not been refined yet.The coordinate system, which contains the calculated object coordinates, is named as the shifted object-coordinate system (SOS).In the next step, these generated object points are used as the "quasi ground control points" to do the bias-compensated bundle block adjustment.The SOS is used as a reference for the orientation, so the adjustment removes the relative but not absolute bias affecting the RPCs.As done in the absolute bias-compensated bundle block adjustment, the relative model also applied the additional affine model to compensate the relative bias for each image, and at least 4 to 6 correspondences are demanded.

Global Relative Pointing Error Correction
According to the projection-trajectory method (Wang, et al., 2010), the projection relationship between image and object space is built via RPCs.Figure 1 presents the sketch of the method.Point P L in left image is projected to two different height level H1, H2 in the object space, and X1 and X2 are the intersections.The projection from image to object space is called the forward projection (FP).Then the object points are projected to the right image.The projection from object to image space is called the backward projection (BP).The projected image points P R ', P R '' are applied to approximate the epipolar curve EP.
Point P R , the corresponding point of P L , is supposed to locate on the epipolar curve EP.When the RPCs are not accurate enough, the distance between the epipolar curve EP and corresponding points P R is the relative pointing error (Franchis et al., 2014).An example of the relative pointing error is exhibited in Figure 2. Several authors (Franchis et al., 2014;Ghuffar, 2016) have reported, that the relative pointing error is modeled as a simple translation in small tiles.Since we will resample the epipolar stereo images without tiling in the later procedure, the global relative pointing error correction is applied.
We select the affine transformation to model the relative pointing error of the whole image scene.The corresponding points distributed in the whole scene of the stereo image are given.First, we choose one of the stereo images as the master image, and generate the local epipolar curves in the slave image according to the project-trajectory method.Then we measure the distances between the corresponding points and the epipolar curves.The global relative point error is formed as: Where, ds, dl are the distances from the corresponding points to the epipolar curves in sample and line directions.s and l are the image coordinates of the corresponding points.A 0 , A 1 , A 2 , A 3 , A 4 , A 5 are the affine model parameters.At last, ε L , ε S are random errors.The affine parameters are estimated by the method of least-square.

MODIFIED PIECEWISE EPIPOLAR RESAMPLING
The projection-trajectory method, which is used to find the corresponding epipolar curves, has been explained in the previous section.By choosing one point on the generated epipolar curve, and redo the same procedure, we can obtain the epipolar curve pair.According to this method, the generated epipolar curve pairs only exist in local area, and they are determined by the height levels for projections.In order to implement the epipolar resampling over the entire image, our work proposes a modified piecewise epipolar resampling strategy based on Oh's (2011) method.Figure 3  1L .According to the character of the satellites imagery, the epipolar curves can be modelled as straight line in small area.So the local epipolar curve is approximated to segment.The epipolar segment grows to a proper length and stop, so that the epipolar curve approximation can be preserved in good quality.The length of the segment is defined by given height range.We select the end point of the segment P 1L ' as the new start point to generate the next segment until the epipolar segments are outside of the original image.In the case of the right image or slave image, we also select the point along y-axis as the initial point S 1R , and we derive the initial epipolar curve's direction EP iniR .When the left image find the start points of epipolar segment (P 1L , P 1L '…), a conjugate points on the right image are set as the start points (P 2L , P 2L '…).The conjugate points are calculated by the projection-trajectory method with a given height level.This height level only affects the x-parallax.
The last step is resampling the epipolar stereo images, and it is shown in Figure 3, step 3.The former steps have built the geometry of satellite's epipolar curve pairs.The interval between the initial points Δy L in left image is set as one pixel.
The interval of the right image's initial points is Δy R , which is set as the length equal to the ground sample distance (GSD) of the left image.As Figure 3, step 2 has shown, the epipolar segments are not located exactly along the direction of x-axis, and for each approximate epipolar curve pair there exists a yparallax.In order to solve this problem, we align each epipolar curve pair to the same row in epipolar image coordinate system.Since the epipolar segments are relocated as a straight line along the x-axis, the stereo images are resampled along the epipolar segments.In left image, the resample distance along the epipolar curve is equal to one pixel like Δy L .For the right image, the resample distance is equal to Δy R .A bi-cubic interpolation is applied for the image resampling.

Data Description
Two different data sets are applied in our work.The first test dataset is a sample QuickBird stereo imagery, which covers the Melbourne, Australia area.The GSD of QuickBird data is 70cm and the size is ca.5500*3500 pixels.The other dataset is a sample of WorldView-2 stereo imagery.The image pair comprises the city area of Munich, Germany, and the GSD is 0.5m.The size of the WorldView-2 stereo image pair is about 9000*24000 pixels.

Experimental Results of Global Relative Orientation
Thirteen even distributed corresponding points are generated by the Envi software for the QuickBird stereo imagery.To apply the relative bias-compensated model, five of these corresponding points are selected to calculate the object coordinates in the SOS.We use five relative control points and the remained eight tie points to accomplish the relative biascompensated RPC bundle block adjustment.We obtain two affine models for the left and right images.For the global relative pointing error correction, we select the left image as the reference.The thirteen corresponding points are applied to calculate the relative pointing errors.First we estimate the affine model for the correction, and then remove the relative corresponding errors in the right image.Another twenty-seven corresponding points are generated in the same way, and they are utilized as check points for the evaluation.
In the case of the WorldView-2 satellite stereo imagery, 28 corresponding points are created by Envi.First we apply the relative bias-compensated model.Twelve of the corresponding points are the relative control points, and the object coordinates in SOS are calculated.Together with the other 16 tie points, we do the relative bias-compensated bundle block adjustment.As to the global relative pointing error correction, we set the left image as the reference again.We then calculate the relative bias between 28 corresponding points and their associated epipolar curves.An affine model is estimated by the relative bias to remove the relative pointing errors.Fifty corresponding points are produced as check points in the WorldView-2 stereo imagery.
The relative pointing errors of the check points are calculated before and after the global relative orientations.The relative pointing error vectors of the check points in the QuickBird stereo imagery are depicted in Figure 4, and The visualization of the relative pointing error vectors are depicted in Figure 5. Considering better visualization, every error vector has been multiplied by 500.In the top-right corner, the sketch gives the scale of one pixel.The red points are the corresponding points and the blue lines indicate the bias to the corresponding epipolar lines.
According to Figure 4(a), the relative pointing error is significant before the relative orientation.After we apply the relative bias-compensated model (Figure 4 For quantitative analysis, the elevation results of both QuickBird (QB) and WorldView-2 (WV2) data are shown in Table 1.In Table 1, the column "Method 1" presents the result of the relative bias-compensated model, and the column "Method 2" shows the result of the global pointing error correction.The column "without RO" gives the result without any relative orientation methods.We find that the RPCs quality of the QuickBird test data is very poor, because it causes about 45 pixels relative pointing errors before the relative orientation.For WorldView-2 data, before the orientation, the maximum bias is one pixel, the minimum bias is 0.007 pixels and the RMSE is 0.45 pixels.It indicates that the RPCs of WorldView-2 data set is accurate.We also find that when we apply two proposed relative orientation methods, the relative pointing error of QuickBird data is decreased significantly and the relative pointing error of WorldView-2 is also reduced slightly.Moreover, the global relative pointing error correction refines the RPCs with the same quality as the relative bias-compensated model.The accuracy of these two methods for two test dataset is very close: the maximum error is smaller than 0.9 pixels and the minimum error is close to zero; the root-mean-square errors (RMSE) of both methods are 0.37 pixels.

Experimental Results of Piecewise Epipolar Resampling
As explained in section 3, we first define the epipolar image coordinate system and find the start points.Then we generate the epipolar curves piecewise as multiple segments.The corresponding epipolar segments are adjusted to the same row in the epipolar image space.The epipolar images are resampled along every epipolar segments.The interval between each epipolar curve and the sample distance along the epipolar segment are both equal to the GSD of the left image.
Following this procedure, we obtain the epipolar stereo image pairs of the QuickBird and WorldView-2 dataset.The two epipolar images can be overlaid to generate the stereo anaglyph.
The generated anaglyph image of QuickBird dataset is displayed in Figure 6.We also generate the anaglyph image of WorldView-2 dataset as shown in Figure 8.The y-parallax can be checked as the difference of the row coordinates between the red (left view) and cyan (right view) colour.
In order to investigate the y-parallax over the entire image, four areas locate at the corners of the epipolar stereo image pair are extracted for both QuickBird and WorldView-2 epipolar images.
The extracted areas are shown in Figure 6 and Figure 8 as yellow numbered rectangles.We generate the analygph subimages for two test datasets.The sub-images of QuickBird data are shown in Figure 7, and the analygph sub-images of WorldView-2 data are displayed in Figure 9.The number at the left-top of each sub-image indicates the related area on the entire image.The yellow lines in the sub-images are assisting to observe whether the corresponding points in the left scene and the right scene are located on the same row.According to Figure 7, the y-parallax of four corner areas is close to zero in the generated QuickBird epipolar images.Figure 9 also shows that all four corner areas achieve very small y-parallax, which proves that the modified piecewise epipolar resampling method works well for the WorldView-2 data.The experimental results show that the modified piecewise epipolar resampling method can generate the epipolar images for both test datasets without tiling.The y-parallax of the epipolar image pair is at subpixel level.

CONCLUSION
To remove the relative pointing error caused by inaccurate RPCs, our work proposes the relative bias-compensated model and the global relative pointing error correction method.Both methods need no GCPs and are applied for the entire image without tiling.Experiments are taken on QuickBird and WorldView-2 data.It has been verified in this paper, that two proposed relative orientation methods can refine the quality of the RPCs for the entire image.Both methods can decrease the relative pointing error to sub-pixel level.
The proposed modified piecewise epipolar resampling method can successfully generate the epipolar curve pairs and resample the epipolar stereo images over the entire image without tiling.
The experiments show that the y-parallax of image's corner area is close to zero.For the entire image, the RMSE of the yparallax is at sub-pixel level.
This paper has built the pipeline to refine the RPCs without ground control information and generate the epipolar images for entire images without tiling.No matter the image is large or it is small, the result can achieve sub-pixel accuracy.In the future, we will take more experiments on the higher resolution and large size satellite stereo imagery to test the proposed methods.We will also verify our methods on the ISPRS benchmark.

Figure 1 .
Figure 1.Projection-trajectory method depicts the steps of the proposed resampling method.The left column shows the procedures for the left image (master image), and the right column exhibits procedures for the right image (slave image).The first step is to define the coordinate system of the epipolar image.According to Figure 3, the centre point of the left image C L is used to calculate the direction of left image's initial epipolar curve EP iniL and its orthogonal line OR iniL .We find four lines EP iniL ', EP iniL '', OR iniL ' and OR iniL '', which are parallel to EP iniL and OR iniL and pass the four corners of the original image.These four lines can define the scope of the left epipolar images.Line OR iniL ' is selected as the y-axis, and line EP iniL ' is selected as the x-axis of epipolar image coordinate system.The intersected point of x-axis X epiL and y-axis Y epiL is the origin point O epiL of the left epipolar image coordinate system.Symmetrically, center point of right image C R derives the initial epipolar curve EP iniR and its orthogonal line OR iniR .Then we define the boundary of the right image and set the origin point O epiR , x-axis X epiR and y-axis Y epiR .The next step is to generate the epipolar curve pair.In Oh's method, he derives the line perpendicular to the epipolar curve of the left image's centre point first.Then the initial points are set along the perpendicular line with a predefined interval (e.g.1000 pixels).The epipolar curve pairs will be generated piecewise from the initial point.Instead of this way, we establish the initial points along the boundary of the left image as shown in Figure 3, step 2. For the master or left image, we select a point S 1L along the y-axis Y epiL as the initial point, extent the line along the direction parallel to the initial epipolar curve EP iniL .When the extented lines reach the boundary of the image, we mark the intersections as a start point P 1L .Then the new local epipolar curves are generated by projection-trajectory method from P

Figure 4 .Figure 5 .
Figure 4. Relative pointing error of QuickBird data: (a) before orientation (b) relative bias-compensated method (c) global relative pointing error correction

Table 1 .
Relative pointing errors of QuickBird dataThe experimental results have proven that, both the relative bias-compensated model and the global relative pointing error correction are accurate and effective methods for relative orientation.No matter the RPCs of the test data is bad (QuickBird) or good (WorldView-2), the proposed relative orientation methods improve the RPCs model and reduce the relative pointing error to the sub-pixel level.

Table 2 .
According to Table2, the RMSE of the yparallax of 40 QuickBird checkpoints is 0.67 pixels.The yparallax's RMSE of 260 WorldView-2 corresponding points are 0.5 pixels.