AN ITERATIVE PIXEL-LEVEL IMAGE MATCHING METHOD FOR MARS MAPPING USING APPROXIMATE ORTHOPHOTOS

Mars mapping is essential to the scientific research of the red planet. The special terrain characteristics of Martian surface can be used to develop the targeted image matching method. In this paper, in order to generate high resolution Mars DEM, a pixel-level image matching method for Mars orbital pushbroom images is proposed. The main strategies of our method include: (1) image matching on approximate orthophotos; (2) estimating approximate value of conjugate points by using ground point coordinates of orthophotos; (3) hierarchical image matching; (4) generating DEM and approximate orthophotos at each pyramid level; (5) fast transformation from ground points to image points for pushbroom images. The derived DEM at each pyramid level is used as reference data for the generation of approximate orthophotos at the next pyramid level. With iterative processing, the generated DEM becomes more and more accurate and a very small search window is precise enough for the determination of conjugate points. The images acquired by High Resolution Stereo Camera (HRSC) on European Mars Express were used to verify our method’s feasibility. Experiment results demonstrate that accurate DEM data can be derived with an acceptable time cost by pixel-level image matching.


INTRODUCTION
Currently, Mars exploration is the hot spot of deep space exploration.Both scientific research and practical engineering applications such as landing sites selection and rover navigation requires high resolution topographic data (Di, 2008).The Digital Elevation Model (DEM) can be derived through the photogrammetric processing of orbital images.Several data sources can be used to generate Mars topographic data, including Visual Imaging System (VIS) on each of the two Viking orbiters (Rosiek, 2005), Mars Orbiter Camera (MOC) on Mars Global Surveyor (Shan, 2005), High Resolution Stereo Camera (HRSC) on Mars Express (Albertz, 2005), and High Resolution Imaging Science Experiment (HiRISE) on Mars Reconnaissance Orbiter (Kirk, 2008).
Mars global DEM is essential data for planetary science.The Mars Orbiter Laser Altimeter (MOLA) DEM is widely used due to its global coverage and high absolute accuracy.However, the grid spacing of MOLA DEM is about 500m at the equatorial regions, which cannot meet many practical application requirements.When compared to other orbital images, HRSC exhibits advantages in image resolution and multi-stereo capabilities.There are nine CCD line detectors mounted in parallel in the focal plane of HRSC camera, which can obtain five panchromatic stereo images and four colour images simultaneously (Albertz, 2005).The HRSC has been imaging the Martian surface since January 2000, and is still on operational up to now.After more than ten years on-orbit operation, the HRSC images have almost covered the entire Martian surface.The HRSC team and interested researchers have been investigated the generation of DEM from HRSC images (Hirschmüller, 2006;Heipke, 2007;Gwinner, 2009;Sidiropoulos, 2016).The early released standard HRSC Level-4 data are single orbit products.

* Corresponding author.
Recently, new global mapping program using the complete HRSC mission data was performed by the HRSC team, and multi-orbit DEM products were derived (Sidiropoulos, 2016).The HRSC team adopted the MC-30 global mapping scheme and subdivided Mars into 30 quadrangles.Currently, the first HRSC MC-30 half-tile was completed, and the grid spacing of DEM and DOM are 50m and 12.5m respectively.However, considering that the maximum resolution of HRSC images is up to 10m, it is expected that higher resolution Mars global DEM can be generated through pixel-level image matching.Moreover, it is very meaningful to develop targeted image matching algorithm according to the terrain characteristics of Martian surface (Wu, 2016).In this paper, we propose a pixellevel image matching method to generate high resolution Mars DEM.In section 2, the disadvantages and advantages for Martian surface image matching are discussed.The basic principle and detailed method are described in section 3. The experiment results are presented in section 4. The conclusion is drawn in section 5.

Disadvantages
Due to the special terrain characteristics, traditional matching method may fail or deliver poor results for Martian surface image matching.The main disadvantages of Martian surface image matching are summarized as follows.
1. Low contrast.As shown in Figure 1, the histogram of HRSC images concentrates in a small range.The low signal-tonoise ratio of Martian surface images will reduce the similarity between corresponding image patches.
2. Insufficient feature points.Feature-based matching methods such as Scale-Invariant Feature Transform (SIFT) are more robust than area-based matching methods.Unfortunately, there are insufficient feature points on Martian surface.On the other hand, though a small amount of feature points can be extracted and used as tie points for bundle adjustment, the point density is not enough for surface reconstruction.
3. Poor image quality.Many factors such as imaging instrument, atmospheric environment and illumination conditions have influences on the image quality.In terms of image quality, Martian surface images perform worse than earth observation images.
4. Repetitive pattern.Repetitive pattern is very common on Martian surface.It will result in wrong matched points especially when area-based matching method is used.To address this issue, some constraint conditions such as epipolar line geometry need to be introduced.Meanwhile, precise approximate value of conjugate points is required.

Advantages
On the other hand, when compared to earth observation images, Mars surface images show some advantages for image matching.
Obviously, there are no trees or rivers on Martian surface.Moving objects such as cars are impossible to appear on Martian surface as well.Additionally, occlusion caused by high building and tall trees on earth observation images can be avoided on Martian surface images.In a word, in terms of terrain continuity, Martian surface images perform better than earth observation images.Thus, our image matching method will make full use of this characteristic.

Basic Principle
A large number of image matching methods have been proposed in the literature.Some practical matching strategies such as hierarchical image matching can be used to perform Martian surface image matching as well.In the recent years, the Semi-Global Matching (SGM) method proposed by Hirschmüller (2008) won great success.SGM was also used to process Mars HRSC images (Hirschmüller, 2006).However, the pixel-level image matching performed with SGM still consumed a lot of processing time, usually several hours for a single orbit.In order to generate high resolution Mars DEM efficiently, targeted image matching method is still required.
Here, we point out that when image matching is performed for DEM generation, the interior orientation (IO) and exterior orientation (EO) data are usually accurately known.Thus, given a stereopair, two orthophotos can be generated through geometric rectification with reference DEM.If the reference DEM is accurate enough, the coordinate displacements of conjugate points at the overlapping area of the two orthophotos will be very small.This interesting and useful information can be used to estimate the approximate value of conjugate points.
In order to introduce the basic principle of our method in detail, image matching on the third pyramid level of the HRSC Level-3 images was performed with SIFT.The results are presented in Figure 2. The HRSC Level-3 images are generated through geometric rectification with a rough DEM, which can be seen as approximate orthophotos.It is observed that most of the pixel coordinates displacements of conjugate points are less than two pixels in X and Y directions.The point 9 (green cross hair) shows the maximum displacements, and results in less than four pixels.
It is noted that point 9 is located at the mountain area.This implies that inaccurate DEM will introduce large coordinate displacements at the overlapping area of two approximate orthophotos.Furthermore, it can be inferred that when image matching is performed at the fourth pyramid level, the pixel coordinates displacements of conjugate points will be less than two pixels.Therefore, at the fourth pyramid level, using the identical pixel coordinates on the left image as approximate value, the conjugate point on the right image can be determined with 5×5 search window.

Data Processing Chain
The main strategies of our pixel-level image matching method include: (1) pixel-level image matching using normalized cross  The detailed processing steps are as follows.
1.The Integrated System for Imagers and Spectrometers (ISIS) developed by USGS is used to perform data pre-processing.The HRSC stereo images with raw PDS format are imported into ISIS using 'hrsc2isis', and the kernels related to the images are determined with 'spiceinit'.
2. Data pre-processing such as contrast enhancement and image pyramids generation is performed.The optional bundle adjustment process is carried out with 'jigsaw'.
3. The interior orientation (IO) data such as pixel size, focal length and initial exterior orientation (EO) data are extracted from SPICE kernels.The scan line exposure time information is acquired with 'tabledump'.Consequently, the rigorous geometric model for HRSC pushbroom images is constructed.
4. The DEM is generated or refined in an iterative processing way.At the lowest pyramid level, the approximate orthophotos are generated using MOLA DEM.The ground points coordinates are used as approximate value of conjugate points.The conjugate points matched on orthophotos are back-projected to original images.Then, through forward intersection, 3D ground points are generated.The forward intersection residuals are used to eliminate wrong matched points.The DEM generated at current pyramid level is used as reference data to generate approximate orthophotos at the next pyramid level.In order to match more points with pixel-level image matching method, a small NCC threshold is used.
5. Through hierarchical image matching, the generated DEM becomes more and more accurate.At last, the grid spacing of the final DEM product can achieve pixel-level.

Image Matching on Approximate Orthophotos:
As illustrated in the top plot of Figure 4, due to the different viewing angles and image scale, it is difficult to perform image matching on the original HRSC Level-2 images.Using IO, EO and rough DEM, approximate orthophotos can be generated.Here, 'approximate' means that the DEM used for geometric rectification is not accurate enough.As shown in the bottom plot of Figure 4, the image distortions are removed through geometric rectification.Moreover, stereo images are resampled with identical Ground Sample Distance (GSD).Therefore, matching on approximate orthophotos is helpful to improve the success rate and accuracy.Furthermore, there is no need to extract feature points due to that at each pyramid level, image matching is performed with pixel-level.

Estimating Approximate Value of Conjugate Points:
The traditional image matching methods use constraint conditions such as epipolar line or affine transformation to estimate or predict the approximate value of conjugate points.
Using our matching strategies, there is no need to estimate the approximate value of conjugate points by complicated technique.
As presented in section 3.1, the ground points coordinates of orthophotos can be used to estimate the approximate value of conjugate points.Due to that the image matching is performed on orthophotos, approximate value of conjugate points is determined with an accuracy of several pixels, depending on the resolution and accuracy of reference DEM used in the geometric rectification process.
Given an image point  on the left image, the pixel coordinates of point  are (, ) and the 2D ground point coordinates of point  are (, ).Thus, (, ) can be calculated with equation 1. Obviously, the point estimation accuracy is influenced by the reference DEM and EO used in the geometric rectification process, assuming that the IO is accurate enough.

Hierarchical Matching:
Hierarchical matching is widely used in most of the image matching algorithms.Similarly, our method adopts this practical matching strategy as well.At each image pyramid level, after the pixel-level image matching is completed, the DEM is generated and used as reference data for the geometric rectification at the next pyramid level.
Suppose the original image resolution of HRSC Level-2 images is 25m and image pyramids are generated with four levels.Thus, the image resolution for 1~4 pyramid levels is 50m, 100m, 200m and 400m respectively.At the highest pyramid level, pixel-level image matching can generate DEM with grid spacing of 400m, which is still higher than that of MOLA DEM.The DEM derived at the fourth pyramid level is used as reference data to generate orthophotos at the third pyramid level.Consequently, through hierarchical image matching, the generated DEM are refined iteratively, and the point prediction accuracy of conjugate points will be more accurate.As shown in Figure 6, in order to perform back-projection for pushbroom images, the best scan line for ground point  must be determined.This requires that the rigorous geometric model for pushbroom images is constructed.Due to the special imaging principle of pushbroom images, each scan line has six exterior orientation elements.The rigorous geometric model is constructed with the extended collinear equation for pushbroom images, and the detailed expressions are given in equation 3. are rotation matrix elements of EO.To search the best scan line for ground point , the geometric constraints in object space can be used to improve efficiency (Wang, 2009).For HRSC pushbroom images, Geng (2014) presented that the computation efficiency of back-projection can achieve one million points per second on a normal PC, which can meet the image matching requirements well.

Forward Intersection for Pushbroom Images:
Using the matched points on original stereo images, 3D ground points are generated through forward intersection.When compared to the frame images, forward intersection for pushbroom images is more complicated.The detailed processing procedures of forward intersection are illustrated in Figure 7. Assuming that  1 and  2 are a pair of conjugate points, ( 1 ,  1 ) and ( 2 ,  2 ) are pixel coordinates for  1 and  2 , ( 1 ,  1 ) and ( 2 ,  2 ) are focal plane coordinates, (  1 ,  1 , − 1 ) and (  2 ,  2 , − 2 ) are image space coordinates respectively, and (, , ) are the 3D ground point coordinates.Firstly, the pixel coordinates are converted to focal plane coordinates by using the transformation equation obtained from the SPICE kernels.Secondly, the exact exterior orientation elements EO1 and EO2 for sane line  1 and  2 are interpolated by using the EO data.Then, the focal plane coordinates are converted to image space coordinates.At last , the ground point ( , ,  ) is calculated with extended collinear equation.It is noteworthy that the ground points are defined in body-fixed Martian Cartesian coordinates.In order to generate DEM, map projection and DEM interpolation are required.

Image space coordinates
Image space coordinates Figure 7.The processing procedure of forward intersection for pushbroom images

Wrong Matched Points Elimination
It is no doubt that wrong matched points cannot be avoided.Therefore, it is necessary to detect and eliminate the outliers.Here, we use the residuals of forward intersection to eliminate the wrong matched points.Given a pair of conjugate points on stereo images, in the forward intersection process there are three unknowns and four observation equations.Consequently, the residuals of ground point coordinates can be calculated.Assuming that the IO and EO are accurate enough, the residuals are determined by the image matching accuracy.As shown in Figure 8, image point  1 and  2 are a pair of conjugate points, and  is the corresponding ground point.If the image point  2 ′ is wrongly determined as the conjugate point of  1 , the forward intersection will deliver large residuals.Therefore, we give a threshold such as two GSD at each pyramid level to eliminate wrong matched points.

EXPERIMENT RESULTS
The software development for image matching and geometric rectification was carried out using Visual Studio 2013 and Qt 5.4.2 on Windows 7 platform.The rigorous geometric model and forward intersection for pushbroom images are implemented based on the open source photogrammetric software DGAP (Stallman, 2008).The tests were performed on a laptop with Intel Core i5 CPU and 8 GB RAM capacity.

Test Datasets
Two orbits HRSC Level-2 images are used to test our method.SPICE kernels and images are ingested into ISIS firstly.Then, image pyramids are generated with four levels by using bi-cubic interpolation.The Level-2 images are orthorectified using equirectangular projection.The Martian reference datum is defined using a sphere body with axis of 3396.19 km.The test datasets information is listed in table 1.The image width of orbit 5273 and orbit 5124 are 5176 and 2584 respectively, and each of two orbits has more than 16000 scanlines.The orbit 5273 images are located at Gale crater, which covers the landing area of Curiosity.The search window size is 3×5 and the match window size is 9×9.Stereo matching is performed between S1 and S2 images.The processing time of the image matching process are measured.It is noteworthy that the forward intersection and DEM interpolation with large-scale point clouds requires a lot of processing time as well.

Experiment Results
Point Prediction Accuracy: One main advantage of our method is using ground point coordinates of orthophotos to predict the approximate value of conjugate points.The pixel coordinates displacements between the predicted value and the matched value are calculated, and the results are illustrated in Figure 9.It is observed that at the original image resolution, the pixel coordinates displacements were less than two pixels.This indicates that our method provides very precise approximate value of conjugate points, which is helpful to improve the image matching efficiency and accuracy.

Discussions
Our image matching method has the advantage that the a priori knowledge at each pyramid level is utilized to the most extent.In addition, the matching strategies take account of the terrain continuity characteristic, and are especially suitable for Martian surface image matching.However, the generated DEM exhibits some deficiencies.At low contrast area, area-based matching may fail and results in pointless region.

CONCLUSION
The objective of this paper is to generate higher resolution Mars DEM through stereo photogrammetry.A targeted pixel-level image matching method for Mars mapping is proposed.The experiment was carried out using only S1 and S2 images.Indeed, more stereo viewing angles images such as nadir, P1 and P2 can be introduced to improve geometric accuracy.Obviously, more images require more data processing time.In order to make our algorithm more practical, the method shall be optimized with parallel processing abilities in the future.Additionally, it is necessary to introduce other image matching strategies into our method to address the topography discontinuity issue at low contrast area.
Figure 2. Image matching results on the third pyramid level of HRSC Level-3 images using SIFT method.

Figure 3 .
Figure 3.The pixel-level image matching and DEM generation procedure.
1) where ( 0 ,  0 ) are the left bottom corner point coordinates of the orthophoto,  and  are the pixel resolution in  and  directions respectively.The pixel coordinates (′, ′) of the conjugate point i' on the right image can be estimated with equation 2. ′ = ( −  0 )/  ′ = ( −  0 )/(2) Conjugate Points: Figure 5 illustrates the matching results on orthophotos and original images.The coordinates transformation from orthophoto to original image is indeed a back-projection process.Given an image point  on orthophoto, the corresponding 3D ground point coordinates are calculated with the DEM used in the geometric rectification process.Through back-projection, the conjugate point ′ on the original image is determined.(a) conjugate points on orthophotos.(b) conjugate points on original Level-2 images.

Figure 5 .
Figure 5. back-projection conjugate points from orthophotos to original images.

Figure 6 .
Figure 6.Best scan line searching for pushbroom image.

Figure 8 .
Figure 8.Using forward intersection residuals to eliminate wrong matched points

Figure 9 .Figure 11 .Figure 13 .
Figure 9. Point prediction accuracy 4.2.1 Image Matching Results: The top plot of Figure 10 shows the image matching results on approximate orthophotos.The matched points on approximate orthophotos are backprojected to original HRSC Level-2 images, and the results are shown in the bottom plot of Figure 10.It is noteworthy that due to low contrast of Martian surface imagery, the success rate of image matching cannot achieve 100% and pointless regions are observed.
correlation coefficient (NCC) on approximate orthophotos; (2) estimating approximate value of conjugate points by using ground point coordinates of orthophotos; (3) hierarchical image matching and DEM generation at each pyramid level; (4) the DEM generated at current pyramid level is used as reference data for the orthophoto generation at the next pyramid level; and (5) fast coordinates transformation from 3D ground points to 2D image points.Due to that the HRSC images are suitable for the generation of Mars global DEM, the pixel-level image matching method is designed for the HRSC images in this paper.Indeed, the method can be used to process other orbital images such as MOC and HiRISE.The procedure of pixel-level image matching method for Mars mapping is illustrated in Figure3.
Pixel-level image matching on orthophotos at pyramid level iBack-projection conjugate points from orthophotos to original images 3D ground points generation by forward intersectionError elimination Error eliminationGenerate DEM at pyramid level i

Table 1 .
Information of test datasets.