IMPLEMENTATION OF A SELF-CONSISTENT STEREO PROCESSING CHAIN FOR 3 D STEREO RECONSTRUCTION OF THE LUNAR LANDING SITES

The department for Planetary Geodesy at Technical University Berlin is developing routines for photogrammetric processing of planetary image data to derive 3D representations of planetary surfaces. The Integrated Software for Imagers and Spectrometers (ISIS) software (Anderson et al., 2004), developed by USGS, Flagstaff, is readily available, open source, and very well documented. Hence, ISIS was chosen as a prime processing platform and tool kit. However, ISIS does not provide a full photogrammetric stereo processing chain. Several components like image matching, bundle block adjustment (until recently) or digital terrain model (DTM) interpolation from 3D object points are missing. Our group aims to complete this photogrammetric stereo processing chain by implementing the missing components, taking advantage of already existing ISIS classes and functionality. We report here on the current status of the development of our stereo processing chain and its first application on the Lunar Apollo landing sites.


INTRODUCTION
One of the main goals in planetary geodesy is to retrieve a threedimensional (3D) description of planetary surfaces and bodies.Many space exploration missions set out their journey having visual sensors (cameras) on board to observe the foreign worlds.From these sensors, many images of planetary bodies have been captured and are available to use.These planetary images can be used to produce cartographic products including digital terrain models (DTM), map-projected imagery and in some cases even 3D models of the corresponding bodies.These data products are suitable for further science analysis and mission planning.
Our group at the department for Planetary Geodesy at Technical University Berlin is developing routines and softwares.Our aim is to photogrammetrically process the available planetary image data, derive 3D representations and provide the results to the scientific community.The Integrated Software for Imagers and Spectrometers (ISIS) software (Anderson et al., 2004), developed by USGS, Flagstaff, was chosen as a prime processing platform and tool kit.It is a specialized image processing package which provides basic image processing operations as well as mission specific data operations.However, ISIS does not provide a full photogrammetric stereo processing chain.Several components like image matching, bundle block adjustment (until recently) or digital terrain model (DTM) interpolation from 3D object points are missing.Our group aims to complete this photogrammetric stereo processing chain by implementing the missing components, taking advantage of already existing ISIS classes and functionality.
In this research we present our stereo processing chain and its first application on the Lunar Apollo landing sites.The processing chain contains several different steps namely dense image matching, object point estimation and DTM interpolation, resulting in local surface representations.Absolute global position of this local DTM is refined by applying the Lunar Orbiter Laser Altimeter (LOLA) (Smith et al., 2011) data set.Currently, LOLA provides the accepted global reference for the Earth moon.

STEREO PROCESSING CHAIN
To compute accurate 3D models and subsequently ortho-image maps several methods and software were developed.At the current stage this includes a software for automated matching, an object point calculation tool and a 3D interpolation tool.The general workflow and components needed to derive DTMs and ortho-images is displayed in Figure 1.After the identification of image sets that cover the area of interest, these image pairs are matched and dense depth maps which are defined for all pixels in the overlapping areas are produced.The resulting depth maps are later used as input for the object point calculation tool and 3D objects points coordinates are estimated by applying the collinearity equations.The large numbers of 3D object point coordinates are fed into the DTM interpolation tool.Here the point coordinates are map-projected applying existing ISIS functionality.The single pixel value for the final DTM can be determined by choosing between several interpolation methods like mean, median, inverse distance weighting (IDW), nearest neighbor (NN), intersection accuracy weighting (IAW).The next step is to project the images onto the derived DTMs to produce map-projected ortho-images.After this step the DTMs and ortho-images of area of interest are achieved.However due to the uncertainties in the spacecraft attitude, position, or instrument mounting the derived surface models are likely displaced with respect to the global reference shape models.To detect and minimize these misalignments, the resulting DTMs are co-registered to Lunar Orbiter Laser Altimeter (LOLA) (Smith et al., 2011) tracks that intersect the study area.LOLA is a precise and accepted global geodetic grid of the Moon (Smith et al., 2011) and provides a good reference shape model.The shift values from the co-registration result are applied to the DTMs and ortho-images to achieve the global positions of them.Finally the shifted DTMs and ortho-images are merged and final DTM and ortho-image mosaics are produced.The matching software is optimized for both orbital and closeranged planetary images and compatible with ISIS formats as well as other common formats like Vicar, TIFF, PNG or JPEG.It supports multithreading in order to increase the performance and to handle large images, such as Lunar Reconnaissance Orbiter Camera (LROC) data, efficiently.The aim of this software is to provide accurate and high density conjugate point clouds that will serve as input for the following steps.Since the accuracy of the final 3D reconstruction highly depends on the precise identification of conjugate points in the stereo images, special care was taken during the development of the image matching software.
The matcher is designed to handle uncompressed, radiometrically corrected and non-rectified images.Thus a pre-rectification and an existing DTM of the study area are not necessary.Figure 2 shows the internal workflow of the matcher.It integrates different matching algorithms like a feature based (FB) and an areabased (AB) matching algorithm (Rodehorst, 2004).In order to decrease the high geometrical differences between the nonrectified images, a special approach which grids the input images into defined equal regions is implemented.Each grid undergoes a pre-processing stage and the geometrical differences of these individual grids are minimized by using the FB algorithm.It detects Speeded Up Robust Features (SURF) (Bay et al., 2008) as tie-points and applies a transformation between the grids from these pixel coordinates.SURF are used since they are invariant to scale, rotation and up to some degree to radiometric deformations (Bay et al., 2008).These characteristics make this type of features ideal for planetary images where scale, rotation and radiometric differences are common due to the different spacecraft orientations and illumination conditions.With this process the geometrical relations of the input images are determined and the differences are minimized.An example of this gridding and transformation process can be seen in figure 3. The input images which are shown in figure 3(a) contain a high difference in terms of geometry and disparity.These differences are minimized as shown in figure 3(c) by gridding and applying a transformation with the detected features (fig.3(b)).Dense area based matching can be applied in the resulting grids.This pre-processing process increases the performance and the accuracy of the dense matching and makes it possible to match non-rectified images.Each pixel with its neighborhood is compared with the other patches within a search range from the other pair with Equation 1.The correspondings are estimated with a winner-take-all approach (highest cross-correlation coefficient).This delivers matching results with a 1 pixel accuracy.The result is then refined to sub-pixel accuracy with LSM which is considered to be the most accurate image matching technique with an accuracy of up to 0.01 pixels (Ackermann, 1984).It is an adjustment approach to minimize the pixel value difference between two patches by applying a geometrical and radiometrical model.The patches are treated like discreet functions.The goal is to minimize the differences between the functions of two patches.This defines the functional model and the whole system is solved with respect to the parameters of the chosen geometrical and radiometrical model.As the geometrical model Affine transformation (Eq.2) is used and the transformation parameters, that minimize the difference, are calculated by applying the Gauss Markoff model of Least Squares Adjustment (LSA) (Eq. 3) (Plackett, 1950).The transformation parameters represent the sub-pixel shift between the patches.
x = a1x + a2y + a3 (2a) where a1, a2, a3, b1, b2, b3 = six unknown parameters x , y = transformed coordinates x, y = initial coordinates where = the observations x = the initial parameters A = the design matrix v = the error matrix added to the system Equation 2 shows the general Affine transformation with unknown parameters.These parameters are the unknowns to be estimated in the adjustment and they are updated for each iteration.The observations ( ) in equation 3 are the pixel-wise differences of the patches and its size is equal to the number of pixels of the selected patch size.Thus the size of the design matrix (A) depends on the chosen patch size and the number of unknowns which is six in our affine transformation implementation.The new parameters are calculated for each iteration by equation 4. When the change in the parameters in two consecutive iteration drops below a threshold, the process stops and the current parameters are considered to be the final results.
where xnew = the new transformation parameters x = the initial parameters A = the design matrix = the pixel value difference of the patches After the matching of transformed grids, the result is back transferred with respect to the corresponding transformation model.This pre-processing and matching stages are performed for each grid in the template image and the results from each grid are merged to achieve the final matching result.

Object Point Calculation
A software which is capable of calculating the object points from the matching results is being developed.It uses the projective geometric relationship between the images and the captured scene.Figure 4 (Julia, 2011) shows the basic idea of projecting a scene to the image plane.In the figure 4 point O is the camera center with the focal length of c; x, y, z is the camera and X, Y, Z the global coordinate system.P is a point in the scene in the global coordinate which is projected to point p in the image plane.It can be noted that the ray between the camera center and the object point is the difference between the vector of object point and the vector of camera center.This geometric relation can be expressed as in equation 5.The rotation matrix and the camera positions can be retrieved from SPICE (Acton, 1996) kernels.The scale factor (m) is an unknown value which varies for each object point.If only one image is available then only the direction to an object point (P ) can be determined but not its absolute position.The object coordinate of point P can only be computed if this ray intersects with other rays from other images in space.That is why two or more images which provide two or more rays for each object points are needed.
where PX , PY and PZ = the object point coordinates OX , OY and OZ = camera center coordinates R = the rotation matrix between the camera coordinate and the global coordinate system m = the scale factor px, py = the coordinate of image point p in image coordinate system c = the focal length of the camera Figure 5 depicts the situation if we only consider 2 images.In an ideal case the two rays would intersect in space defining the object point at this intersection.However, due to the unavoidable uncertainties from camera position, camera orientation, camera distortion, or correspondence detection, these two rays do not intersect.In order to find the point that is close to both rays, the values of the scale factors (m1 and m2) that correspond to the minimum distance between the rays should be solved by minimizing the normal equation 6.
where X 2 = the norm to be minimized O1X , O1Y and O1Z = the first camera center O2X , O2Y and O2Z = the second camera center R1 and R2 = the rotation matrices between the cameras and the global coordinate system m1 and m2 = the scale factors for each camera p1x, p1y =the image coordinate of object point P in the first image p2x, p2y = the image coordinate of object point P in the second image c1, c2 = the focal lengths of the cameras To solve the ray intersection problem, equation 6 is differentiated with respect to m1 and m2 and the result is set to zero.After solving it for m1 and m2, these values are plugged into the ray equations to obtain the point for each ray that is closest to the other ray.Mid point of the two point locations yields the depth estimate.
The object point calculation tool parses all the matched points from the previous step of the stereo processing chain and applies the explained approach to estimate the object points.It uses already existing camera models within the ISIS frame.SPICE kernels provide the camera interior and exterior orientations.A first implementation is finalized and needs further testing.It is also envisaged to integrate a bundle adjustment module within this tool to minimize the uncertainties of camera orientations and avoid the possible offsets and distortion between the derived surface model and the absolute frame.

Interpolation Tool
The 3D object point coordinates, that are derived from the matching and subsequent object point calculation, are then used for the DTM interpolation.The body-centric coordinates are mapprojected into a pre-defined cube file, which serves as a target container.The pixel values of the target cube defined by the object point height at this position.It is possible that several object points define the value of one pixel of the target position.Different interpolation methods like mean, median, inverse distance weighting (IDW), nearest neighbor (NN), intersection accuracy weighting (IAW) are implemented to determine exactly one value for the resulting pixel.The input data can be provided in non-sequential order and there are no specific requirements in terms of spatial distribution or homogeneity of the distribution of the points.Furthermore the interpolator has a gap filling feature and export capability to different formats like CUB, TIFF, Erdas RAW. Figure 6 schematizes the DTM creation process.The object points are converted from 3D body-centric coordinates to map coordinates and a template DTM is filled with the values from coinciding points by applying different interpolation methods.
Figure 6: Interpolation Tool -The points with different colors represents the object points that fall into one pixel

TESTS
The stereo chain was tested with the Apollo 15 landing site.We used LROC images of 1576 and 1577 with 0.51 meters per pixel ground resolution.The two right images, the two left images and one right/left image pair were matched to provide DTMs covering parts of the Apollo 15 landing site area (fig.7).Table 1 shows the statistical summary of the image matching and triangulation results.As can be seen, the rate of successful matches is very high and thus the coverage of the depth maps are satisfactory.Moreover, the average triangulation error are relatively low.This error states a relative error and should not be confused with the absolute error, that relates to the uncertainty with respect to a global reference.It only suggests that the object points are in a good consistency and says nothing about the deviation of DTMs from the real surface values.

Co-Registration to the LOLA Reference
While the DTMs benefit from the very good internal consistency, small offsets and possible distortion between the derived surface models remain.The models are also not referenced to a common reference frame as uncertainties in the spacecraft attitude, position, or instrument mounting on the spacecraft are present.While a bundle block approach to correct for these small misalignments is currently under development, we have studied a different approach using LOLA data as reference.The resulting DTMs are co-registered to Lunar Orbiter Laser Altimeter (LOLA) tracks that intersect the study area.During this process a grid search is performed looking for a best fit between the LOLA profile and the DTM heights on pixel accuracy level.This is followed by a least squares fit of the track with respect to the DTM to obtain sub-pixel accuracy for the fit.More information about the LOLA co-registration can be found in Gläser et al. (2013).Results of the fit between DTMs and LOLA tracks are shown in table 2. These shift parameters are applied to each DTM and orthoimage to register them to the LOLA tracks.The shifted DTMs and ortho-images are then mosaiced with the existing ISIS tools and final products of the Apollo 15 landing site is achieved.

Figure
Figure 1: Stereo processing flowchart 2.1 Matching Software

Figure 4 :
Figure 4: Imaging of a surface

Figure 5 :
Figure 5: Two images imaging geometry

Table 1 :
Statistical summary of the image matching results.

Table 2 :
The shift values of the DTMs with respect to LOLA tracks