JIGSAW : THE ISIS 3 BUNDLE ADJUSTMENT FOR EXTRATERRESTRIAL PHOTOGRAMMETRY

The Integrated Software for Imagers and Spectrometers (ISIS) package was developed by the Astrogeology Science Center of the U.S. Geological Survey in the late 1980s for the cartographic and scientific processing of planetary image data. Initial support was implemented for the Galileo NIMS instrument. Constantly evolving, ISIS has since added support for numerous missions, including most recently the Lunar Reconnaissance Orbiter, MESSENGER, and Dawn missions, plus support for the Metric Cameras flown onboard Apollo 15, 16, and 17. To address the challenges posed by extraterrestrial photogrammetry, the ISIS3 bundle adjustment module, jigsaw, is evolving as well. Here, we report on the current state of jigsaw including improvements such as the implementation of sparse matrix methods, parameter weighting, error propagation, and automated robust outlier detection. Details from the recent processing of Apollo Metric Camera images and from recent missions such as LRO and MESSENGER are given. Finally, we outline future plans for jigsaw, including the implementation of sequential estimation; free network adjustment; augmentation of the functional model of the bundle adjustment to solve for camera interior orientation parameters and target body parameters of shape, spin axis position, and spin rate; the modeling of jitter in line scan sensors; and the combined adjustment of images from a variety of platforms. * Corresponding author.


INTRODUCTION
The Integrated Software for Imagers and Spectrometers (ISIS) package (http://isis.astrogeology.usgs.gov/)was developed by the Astrogeology Science Center (ASC) of the U.S. Geological Survey (USGS) in the late 1980s for the cartographic and scientific processing of planetary image data (Gaddis et al., 1997;Anderson et al., 2004).Initially supporting the Galileo NIMS instrument, ISIS has since added support for numerous missions including the Viking orbiters; Mariner 9/10; Clementine; Voyager I/II; Galileo SSI; Pathfinder IMP; Mars Global Surveyor MOC, TES, and MOLA; 2001 Mars Odyssey THEMIS; Lunar Orbiter; Cassini ISS, VIMS, and RADAR; Mars Express HRSC; Mars Exploration Rover (MER) Pancam, Navcam, and Micro-Imager; and Mars Reconnaissance Orbiter (MRO) HiRISE and CTX.ISIS was originally developed in FORTRAN for the VAX/VMS operating system, requiring hardware display buffers for interactive graphics.Programming of ISIS2, beginning in 1992, was primarily in C and FORTRAN using the X-Window Application Programming Interface for the UNIX operating system.This was modernized with the move in late 2001 to ISIS3, using C++ as the primary programming language.Though ISIS2 is still used to support some missions and datasets, the work discussed in this paper relates specifically to ISIS3.Support has been recently added for the Lunar Reconnaissance Orbiter (LRO), MESSENGER, and Dawn missions, and for the Metric Cameras flown onboard Apollo 15, 16, and 17.
The digital image mosaic and digital terrain model (DTM) are two key cartographic products generated from planetary image data.Errors in spacecraft navigation and attitude data lead to misregistration between images in the mosaic and to errors in the DTM.To minimize these errors, images are controlled using the least-squares bundle adjustment, well known as a key element in the photogrammetric workflow (Brown, 1958).Given a set of two-dimensional point feature measurements from multiple images, this process, in its most general form, simultaneously solves for parameters of image position and orientation in addition to triangulated ground point coordinates.The ISIS3 bundle adjustment module, jigsaw, has undergone recent significant development and continues to evolve to meet the challenges posed by extraterrestrial photogrammetry.In the following, we begin with a discussion of some of these challenges and then report on the current state of jigsaw.Examples from the ongoing processing of Apollo Metric Camera images and images from recent missions such as LRO and MESSENGER are given.Finally we close with a summary of future plans for jigsaw.

CHALLENGES IN EXTRATERRESTRIAL PHOTOGRAMMETRY
Extraterrestrial photogrammetry is often characterized in terms of the challenges it presents.Kirk, Howington-Kraus, and Rosiek (2000) divide these into three primary categories; multiplicity of datasets, novel sensors, and unusual geometric and coverage properties.
Planetary imaging systems are typically developed for a single mission, and historically have not been designed for the purpose of rigorous geodetic mapping (notable exceptions being the

THE PHOTOGRAMMETRIC PROCESS IN ISIS3
After ingestion into ISIS3, images are typically calibrated radiometrically and may undergo additional processing for noise removal.Rigorous sensor models define the relationship between ground and image and include system and camera calibration parameters.Initial estimates of image EO parameters are from spacecraft navigation and attitude data which is available in ISIS3 in the form of SPICE (Spacecraft, Planet, Instrument, C-matrix (pointing), and Events) kernels consistent with the standards of NASA's Navigation and Ancillary Information Facility (NAIF) (Acton, 1996).In the past (prior to the mid 1990's), operators measured tie points between overlapping images manually.This is now accomplished primarily through automated image matching techniques.Operators can still measure tie points manually, if necessary, and manually correct errors occasionally produced by automated methods.Prior to bundle adjustment, initial ground coordinates of tie points are generated by intersection with an available shape model, such as a sphere, ellipsoid, or global, regional, or local DTM.
At the time of writing, the correction of jitter in line scan sensors has not been fully addressed within ISIS.The University of Arizona (UA) has developed a method to correct for jitter in HiRISE images of Mars.This approach utilizes a jitter modeling tool created at UA together with dedicated ISIS tools (Mattson et al., 2009).The HiRISE camera consists of 14 Charge Coupled Devices (CCDs) arranged in a staggered array on a fixed focal plane.The cross-track overlap and along-track time separation of the CCDs are used to collect information on spacecraft motion.Pixel offsets are measured between feature positions in the overlap area of CCD pairs.From these offsets, jitter is estimated in the x and y directions and interpreted in terms of camera pointing angles.These are combined with existing pointing data to create an improved NAIF camera pointing kernel.For applications such as map projection, bundle adjustment, and stereo analysis, images resampled with the improved kernel are ideally jitter-free.The smooth orientation data can be bundle adjusted as usual to bring the image into agreement with ground control.This technique may also be applicable to LRO narrow-angle camera (NAC) images.It depends crucially on the fact that stereo parallax can be neglected when comparing features seen in multiple detectors during a single observation.To model jitter in line scan cameras with significant stereo convergence between detectors (e.g.HRSC) it will be necessary to solve simultaneously for ground point coordinates and jitter motion.

JIGSAW
Jigsaw is currently capable of adjusting image data from frame and line scan cameras and the Chandrayaan-1/LRO Mini-RF radar sensor.Although Magellan SAR and Cassini RADAR images can be processed in ISIS3, they are produced in mapprojected form and are handled as maps rather than images.Currently, only Mini-RF radar images can be bundle adjusted (Kirk and Howington-Kraus, 2008).
Using jigsaw, one may solve for camera pointing alone, camera position alone, or pointing and position simultaneously.Three dimensional coordinates (latitude, longitude, and radius) of all ground points are also determined in the adjustment.Through rigorous weighting, parameters may be held fixed, allowed to adjust freely, or constrained with a priori precision information.
The functional model for the bundle adjustment, defining the relationship between image and object space coordinate systems, is the standard collinearity condition or its radar equivalent, in which a pixel corresponds to a known circle in space rather than to a line.This model is modified for timedependent line scan and radar sensors to allow for changes in EO parameters between neighboring lines (e.g.Ebner, Kornus, and Ohlhof, 1992;Fraser, 2000).In this case, EO parameters are modeled as functions of time.In jigsaw, individual 2 nd or higher order polynomials are utilized as interpolation functions for each EO parameter.In the adjustment of line scan images, it is well-known that high correlations typically occur amongst polynomial coefficients, notably between position and attitude parameters.To control these correlations, the selection and weighting of these parameters must be carefully considered (Fraser, 2000).Alternatives to polynomials such as spline functions are being considered for more accurate representation of EO parameters over long scans.
Several notable improvements to jigsaw have recently been implemented.As mentioned above, the rigorous weighting of parameters allows the a priori precision of spacecraft position, attitude data, and ground point coordinates to be included in the adjustment if available.Similarly, the weight of an image point measurement in the bundle adjustment is computed from its associated a priori precision.In ISIS3, a single precision value can be assigned to all image measurements with the option of specifying any individually.In this way, image points of different precision (e.g., manual vs. automated measurement) can be weighted appropriately.Automated outlier rejection of image measurements via robust M-estimation techniques reduces the need for tedious manual residual analysis (Koch, 1999).The variance-covariance matrix from rigorous error propagation yields estimates of uncertainties for the estimated parameters and any correlations between them.Such estimates are critical to the proper comparison of datasets from different instruments and missions.
With the ever-increasing amount of image data produced by planetary imagers, perhaps the most significant developments in jigsaw are those that have improved computational efficiency and reduced memory requirements.This in turn has dramatically increased productivity.Adjustments that may have previously required many hours or even days are completed significantly faster, for some cases in just minutes.These improvements take advantage of the sparse structure of the normal equations matrix common to aerial and satellite-based photogrammetric networks.Prior to describing these sparse matrix methods, we briefly review least-squares estimation in the Gauss-Markov model and the solution to the normal equations.

Least-squares Estimation in the Gauss-Markov Model
The subsequent development follows closely Uotila (1986).
Given an m x 1 vector l of observables and an n x 1 vector x of unknown parameters, the linear (or linearized) functional model relating the two is Where ε is an m x 1 vector of true errors and A is the m x n design matrix (m ≥ n).Generally, A is of full rank.Given that the expectation of the error vector in (1) is zero, the Gauss-Markov model with full rank is: ( ) ( ) where D is the dispersion operator, l ∑ is the positive-definite covariance matrix of the observations, 2 0 σ is the a priori variance of unit weight, and P is the observational weight matrix.The residual vector is given by The "cap" symbol ( ˆ) indicates an estimated value.The leastsquares method yields the estimates x ˆ and 2 0 σ such that the Euclidean length of the residual vector is minimized.This leads directly to the normal equations: The estimates x ˆ and 2 0 σ are where n and u are the numbers of observation equations and unknown parameters, and n-u is the degrees of freedom.
N is the normal equations matrix.The variance-covariance matrix of the estimated parameters is given by ( )

Solving the Normal Equations
The basic form of the normal equations organized by parameter type facilitates a significant reduction in size of the matrix to be solved through parameter elimination (Brown, 1958).The standard approach is to eliminate point parameters.Following Brown (1976), the normal equation system can be written as Figure 1 shows the normal equations matrix structure for an example network of five images and four points.( ) The solution vector for the image parameters Finally, through back-substitution, the solution vector for the object point parameters is computed from In the typical aerial or satellite photogrammetric network the N matrix and the reduced normal equations matrix are sparse.
In the following section, we discuss how this sparsity may be further exploited.

Sparse Matrix Methods
Jigsaw uses the CHOLMOD package (Chen et al., 2008;Davis, 2006) for solving the sparse, reduced normal equations.CHOLMOD (for CHOLesky MODification), is a highperformance package for the factorization of sparse, symmetric positive definite matrices and for the modification (update/downdate) of an existing sparse Cholesky factorization.
Forming the sparse normal equations matrix efficiently in the bundle adjustment is challenging.The approach adopted here, and the subsequent discussion, follow closely Konolige (2010).
CHOLMOD uses the compressed column storage (CCS) matrix format which is efficient in memory storage but inefficient to fill in an arbitrary manner because the insertion of every nonzero entry requires that all succeeding entries be shifted.To build the reduced normal equations matrix, an interim sparse matrix structure is required that can be efficiently populated in a random fashion and can be traversed by column in row order to subsequently construct the CCS matrix required by CHOLMOD.We use a type of Block Compressed Column Storage (BCCS) which consists of an array of map containers, each representing a column of square matrix blocks in the reduced normal equations (Figure 2).The key into each column map is the block's row index.The value at each key is a square dense matrix with a dimension equivalent to the number of exterior orientation parameters used for the image.Zero blocks are not stored.The BCCS matrix is created only in the first iteration of the bundle adjustment; in subsequent iterations it need only be repopulated.As the normal equations matrix is symmetric, only the upper triangular portion need be stored.

BUNDLE ADJUSTMENT EXAMPLES
In this section we provide examples illustrating the significant time savings resulting from the implementation of sparse matrix methods discussed above.  1 and details of each follow below.

Mercury: MESSENGER
The Mercury Dual Imaging System (MDIS) on the MErcury Surface, Space ENvironment, GEochemistry, and Ranging (MESSENGER) spacecraft (Solomon et al., 2001) consists of the monochrome NAC and multispectral WAC frame cameras which are co-aligned and mounted together on a pivot (Hawkins et al., 2007).MDIS acquired images through three flybys of Mercury and has continued to do so since MESSENGER entered orbit around the planet in March 2011 (Solomon et al., 2008).Using these images, the Astrogeology Science Center (ASC) is constructing a global, controlled monochrome base map and DTM of Mercury with ISIS3 (Becker et al., 2012).
The most recent bundle adjustment of MDIS data consisted of 21,498 NAC and WAC G-filter images acquired from a few days after orbit insertion to 15 November, 2011.The network, one of the largest the ASC has processed, contains 263,762 ground points and 1,857,190 image point measurements (3,714,380 observations -two for each measurement).The 855,780 unknown parameters in the adjustment included the coordinates for all ground points (latitude, longitude, and radius) and three attitude angles for all images.The ratio of observations to unknowns is ~4:1.Attitude angles and radius point coordinates were constrained with a priori precisions of 0.2 degrees and 10,000 m respectively.Latitude and longitude point coordinates were free to adjust.The dimension of the reduced normal equations matrix was 64,494.The network exhibits a high degree of sparsity, indicated by the percentage of non-zero entries (0.12%).The total time per iteration was 104 seconds.Of that, 10 seconds were required by CHOLMOD for matrix analysis, factorization, and solution.

The Moon's South Pole: LRO NAC
With images from the LRO left and right NAC line scan cameras (Robinson et al., 2010), the ASC has produced controlled polar mosaics of the Moon in support of NASA's Lunar Mapping and Modeling Project (LMMP) (Nall et al., 2010).These are the largest geodetically controlled lunar or planetary mosaics ever produced by the ASC (Lee et al., 2012), and we believe may be the largest such mosaics (in numbers of pixels) ever produced.
The South Pole network contains 3827 images, 527,756 ground points, and 3,363,623 image measurements (6,727,246 observations).The number of unknown parameters was 1,606,230 and the ratio of observations to unknowns ~4:1.While all position and attitude exterior orientation (EO) parameters were modeled by 2 nd degree polynomials, only the constant terms for both were included in this adjustment.
Position and attitude parameters were constrained with a priori precisions of 1000m and 3.0 degrees respectively.Thirty-seven ground points were constrained in all three coordinates to a Lunar Orbiter Laser Altimeter (Smith, 2010) derived DTM.For the remainder, only the radius coordinate was constrained to the DTM.The dimension of the reduced normal equations matrix is 22,962 and 3% of the matrix entries are non-zero.Time per iteration was 724 seconds and CHOLMOD required 121 seconds for matrix analysis, factorization, and solution.This network contains roughly a sixth of the images in the MESSENGER network, but twice as many ground points and image measurements.A definitive explanation as to why the solution time here seems out of proportion to that of the MESSENGER example and the Apollo Metric example to follow will require more study.It seems likely however that the explanation must in part be related to the relatively higher number of non-zero elements in this matrix.

The Moon: Apollo Metric Cameras
In the 1970s, an integrated photogrammetric mapping system was among the experiments in the scientific instrument module (SIM) flown on Apollo Missions 15, 16 and 17.This included a frame Metric Camera (MC) for mapping, a high-resolution Panoramic Camera, and a star camera and laser altimeter to provide supporting data.All SIM cameras were film-based.Nearly 25% of the lunar surface was photographed in stereo at resolutions of one to a few meters (Livingston et al., 1980).Later, significant effort was put forth by the Defense Mapping Agency, the National Oceanic and Atmospheric Administration, USGS, and NASA to geodetically control these photographs (Doyle et al., 1976).In 2007, NASA's Johnson Space Center and Arizona State University (ASU) scanned the original MC negatives at film-grain resolution (Robinson et al., 2008).This work has facilitated an ongoing collaboration between the ASC, the Intelligent Robotics Group of the NASA Ames Research Center (ARC), and ASU, to photogrammetrically and geodetically control that subset of these images with nadirpointed geometry and in the future, to include oblique images.
The processing of the nadir-pointed imagery has been accomplished with the highly automated Vision Workbench and Ames Stereo Pipeline software (developed by ARC) (Moratto et al., 2010)  As a more definitive illustration of the time savings achievable from the application of sparse matrix methods, we cite one further Apollo Metric example for which we have available processing times from a previously implemented dense matrix Cholesky approach.This control network consists of 1462 Apollo 15 images and 168,182 ground points.The time to factorize and solve the reduced normal equations in this case was 31 minutes using the dense Cholesky method.Using CHOLMOD this time was reduced to 0.59 seconds.

FUTURE WORK
Here we summarize some capabilities planned for future implementation in jigsaw.Efficient processing of the ever larger data sets generated by planetary missions today requires automating many of the tasks traditionally performed by a human operator.For the bundle adjustment this can be accomplished in part through the implementation of techniques such as robust estimation and sequential estimation.As mentioned in Section 4, robust estimation via M-estimators has recently been implemented in jigsaw.Sequential estimation facilitates the update of an existing data set with new measurements as they become available and the removal of rejected measurements.Free network adjustment (Blaha, 1982) is advantageous for the analysis of measurement residuals free of contamination from errors in constrained parameters.
Variance component estimation (Koch, 1986) is important as mixing observations from different sources with unknown precisions is commonplace.We will implement the ability to solve for camera interior orientation parameters and target body parameters of shape, spin axis position, and spin rate.We also intend to pursue the integration and adjustment of data from sensors on multiple platforms including orbiters, descent stages, landers, and rovers.This involves not only multiple coordinate frames but a wide variety of accurately known constraints and less accurately known transformations between them.As mentioned previously, alternatives to general polynomials such as spline functions are being investigated for the representation of EO parameters in images from time-dependent sensors.Finally, based on promising preliminary tests on short images using polynomial fitting, we plan to implement jitter modeling for line scan sensors within the bundle adjustment via spline fitting.

ACKNOWLEDGEMENTS
We gratefully acknowledge the MESSENGER mission; the LRO mission, the LROC Science Team, and the Mini-RF Science Team; Mark Robinson from the School of Earth & Space Exploration at Arizona State University; and the NASA Ames Intelligent Robotics Group, in particular Terry Fong, Zachary Moratto, Ara Nefian, and Michael Broxton for their invaluable support and collaboration in the processing of Apollo Metric Camera images.This work has been funded under several NASA programs, including the Planetary Figure 1 shows the normal equations matrix structure for an example network of five images and four points.N & and c & are contributions to the normal equations with respect to image parameters and N & & and c& & with respect to point parameters.The solution vectors of image and point parameters are x & and x& & .N & and N & & are block diagonal.N & & consists of 3x3 submatrices along the diagonal, one for each ground point.N & consists of submatrices on the diagonal for each image, the size of each determined by the number of EO parameters being solved for the image.The elimination of x& & from (9) yields the following reduced normal equation system, the size of which depends only upon the number of image parameters.
Figure 1: Structure of full normal equations matrix for an example network with five images and four points.Gray blocks are non-zero, white are zero.

Figure 2 :
Figure 2: Five-image network illustrating the relationship between an upper triangular matrix in dense storage (left) and the same matrix in Block Compressed Column Storage (BCCS) (right).Zero-blocks (in white on the left) are not stored in BCCS format.Numbers inside blocks on the right are row indices.
variety of sensors separated by long time periods is often required.Beyond the relatively simple frame camera, these include line scan, push frame, and radar sensors for which the time variation of orientation parameters must be modeled.For line scan and push frame sensors, the accurate determination of image position and orientation may also be complicated by unmeasured, high frequency spacecraft motions, or jitter.Sensors may be mounted on orbiting spacecraft or at the surface on stationary platforms or roving vehicles.Images scanned from film cameras on missions long past may be processed together with recent digital images.
Metric Camera system flown onboard Apollo missions 15-17 and the High Resolution Stereo Camera (HRSC) on the Mars Express spacecraft).Processing combinations of images acquired from a Bundle adjustment results are given for MESSENGER narrow-angle camera (NAC) and wide-angle camera (WAC) images of the planet Mercury, LRO NAC images of the South Pole of the Moon, and a combined network of images from the Metric Cameras flown onboard Apollo missions 15 and 16.We focus here only on the times required per iteration and for the solution of the reduced normal equations.All examples were run on an Intel Xeon X5650 with 16 GB of RAM, 12 MB of cache, and a clock speed of 2.66 GHz.The examples are summarized in Table

Table 1 :
Jigsaw bundle adjustment times for iteration and sparse matrix solution.
used in conjunction with ISIS3.Here we report bundle adjustment results for a combined network of Apollo 15 and 16 images.Work on adjusting images from all three missions using jigsaw is ongoing.The Apollo 15/16 network contains 2664 images, 306,457 ground points, and 922,573 image measurements (1,845,146 observations).Both position and attitude EO parameters were included in the adjustment.There are 935,355 unknown parameters and the ratio of observations to unknowns is ~2:1.The dimension of the reduced normal equations matrix is 15,984.2% of the matrix entries are non-zero.Time per iteration was 65 seconds.Matrix analysis, factorization, and solution were completed in 2 seconds by CHOLMOD.