Students T Robust Bundle Adjustment Algorithm

Student’ s T Rob ust Bundle Adjustment Algorithm Aleksandr Y . Aravkin Uni versity of British Columbia V ancouver , BC saravkin@eos.ubc.ca Michael Styer Stanford Uni versity Stanford, CA mstyer@cs.stanford.edu Zachary Moratto N ASA Ames Moffett Field zachary.m.moratto@nasa.gov Ara Nefian Carnegie Mellon Uni versity and N ASA Ames Mof fett Field, CA ara.nefian@nasa.gov Michael Broxton Carnegie Mellon Uni versity and N ASA Ames Mof fett Field, CA michael.broxton@nasa.gov Abstract Bundle adjustment (BA) is the pr oblem of r efining a vi- sual r econstruction to pr oduce better structure and vie wing parameter estimates. This problem is often formulated as a nonlinear least squares pr oblem, wher e data arises fr om inter est point matching . Mismatched inter est points cause serious pr oblems in this appr oach, as a single mismatch will affect the entir e reconstruction. In this paper , we pro- pose a novel r obust Student’s t BA algorithm (RST -BA). W e model r epr ojection err ors using the heavy tailed Student’ s t-distribution, and use an implicit trust r e gion method to compute the maximum a posteriori (MAP) estimate of the camera and viewing parameters in this model. The result- ing algorithm exploits the sparse structur e essential for r e- constructing multi-image scenarios, has the same time com- plexity as standar d L 2 bundle adjustment ( L 2 -BA), and can be implemented with minimal chang es to the standar d least squar es framework. W e show that the RST -B A is mor e ac- curate than either L 2 -BA or L 2 -BA with a σ -edit rule for outlier r emoval for a range of simulated err or gener ation scenarios. The new method has also been used to r econ- struct lunar topogr aphy using data fr om the NASA Apollo 15 orbiter , and we present visual and quantitative compar- isons of RST -B A and L 2 -BA methods on this application. In particular , using the RST -B A algorithm we were able to r e- construct a DEM fr om unpr ocessed data with many outliers and no gr ound contr ol points, which was not possible with the L 2 -BA method. 1. Introduction Bundle adjustment is a large sparse geometric parame- ter estimation problem, in which parameters are 3D feature coordinates and camera poses. Classically , b undle adjust- ment is formulated in the nonlinear least squares frame work (see [ 28 ] and sources cited within). The goal of bundle adjustment is to refine a visual reconstruction by identify- ing a sparse cloud of tie-point features in multiple images, matching tie-points common to several images, and adjust- ing camera and world point parameters simultaneously to improv e the reconstruction. Prior information about cam- era parameters and ground control points are incorporated using the Bayesian modeling framework, and camera pa- rameters and 3D world coordinates of the tie-point features are estimated using the matched tie-point features as data. The moti vating problem for our work in robust bundle adjustment is to produce a robust reconstruction in the pres- ence of outliers generated due to the misidentification of tie-points across images. It is virtually impossible to en- sure that automated algorithms for finding and matching tie- point features always generate correct correspondences. A single wrong identification can lead to a large phantom er- ror which dominates other ‘good’ data in the least squares framew ork. T o deri ve our approach, we be gin with a statistical model for the reprojection errors in pix el space as well as ini- tial camera parameter errors and ground control point er - rors, and then find a maximum a posteriori estimate for this model using an implicit trust region optimization method. W e model both reprojection errors and prior uncertainty on cameras and ground control points as distributed according 1 − 4 − 3 − 2 − 1 0 1 2 3 4 v N (0 , 1) L (0 , 1) T ( ν = 1) − 4 − 3 − 2 − 1 0 1 2 3 4 v 0 . 5 v 2 k √ 2 | v k | log(1 + v 2 k ) − 4 − 3 − 2 − 1 0 1 2 3 4 v v k √ 2 v k / | v k | 2 v k / (1 + v 2 k ) Figure 1. Gaussian, Laplace, and Student’ s t Densities, Corresponding Negati ve Log Likelihoods, and Influence Functions (for scalar v k ) to the Student’ s t-distribution. There are many approaches that try to incorporate the idea of ‘rob ust cost’ or ‘fat tails’ into the bundle adjust- ment framew ork (see [ 28 ], [ 11 ]). Most of these schemes require some iterative re-weighting of the least squares ob- jectiv e. Howe ver , these approaches lack a progress metric or any optimality guarantees, which limits possibilities for algorithm design and con vergence criteria. W e begin in- stead with distributional assumptions on reprojection errors, deriv e a corresponding maximum a posteriori (MAP) opti- mization problem, and dev elop an algorithm to solve it. Our approach takes advantage of the sparse structure described in [ 16 ] and [ 12 ], and determines the best camera parame- ters by solving a maximum likelihood (ML) problem. using an algorithm that still follo ws the main idea in [ 18 ] while exploiting the sparse structure in [ 11 ]. The resulting robust Student’ s t B A (RST -BA) preserv es the sparse structure and efficiency of the method used in [ 16 ], allo wing the robust implementation to run in a com- parable time to the standard B A method. Compared with automated outlier remov al schemes and ad-hoc re weighting using ‘robust-cost’ functions, RST -B A is faster and more accurate both on synthetic data and in real applications. Moreov er, it is straightforward to implement with minimal changes to the standard least-squares framew ork. Mismatches in tie-point features is a well known prob- lem, and there are a variety of approaches in the literature for robustification and dealing with outliers. Threshold- based outlier removal is a common approach with many variants ([ 27 ], [ 22 ], [ 20 ], [ 21 ], [ 4 ], [ 9 ], [ 15 ], [ 13 ]). After one or more iterations of bundle adjustment, observations with residuals exceeding some predetermined threshold are remov ed from the dataset and the b undle adjustment opti- mization is run again on the remaining observations. The choice of threshold depends on the particular problem and dataset. [ 27 ], for instance, uses a threshold of 1.5 stan- dard deviations from the mean residual error , while [ 21 ] discards residuals with size to v ariance ratio greater than 3. Generally , some σ -edit rule is applied, with a typical range between 1 to 2 standard deviations. Howe ver , there are problems with using σ -edit rules even in the linear re- gression conte xt [ 19 , Chapter 1], and it is particularly limit- ing in nonlinear problems because outliers affect the initial fit, which is then used for classification of outliers and in- liers. Another common approach is iterative reweighting ([ 8 ], [ 7 ], [ 20 ], [ 23 ]). In this method, the observations are weighted at each iteration in in verse proportion to their residual errors. Observations may also be remov ed accord- ing to a σ -edit rule with a higher threshold [ 20 ]. Closely related to this approach is that of ‘robust cost functions’, which is discussed in [ 11 ] and implemented in [ 3 ] and [ 4 ]. The main problems with ad hoc re weighting according to such cost functions is the lack of a conv ergence theory that guarantees the algorithm will stop at a stationary point of the objecti ve, and in many cases the lack of an explicit ob- jectiv e function. Model-based bundle adjustment is yet another rob ust technique, in which prior information about the surface un- der reconstruction is used to aid in the adjustment [ 7 ]. [ 29 ] proposes a bundle adjustment method which does not in- volv e solving for the camera parameters, in order to reduce numerical instability and decrease sensitivity to noise. W e do not consider these approaches, since we hav e no prior information on the surf ace we are reconstructing, and re- finement of the camera parameters is essential to our appli- cation. The paper proceeds as follows. In Section 2 we formu- late standard b undle adjustment as a nonlinear least squares problem that can be solved with L 2 -B A. In Section 3 we de- scribe the multiv ariate Student’ s t-distrib ution. In Section 4 , we formulate robust bundle adjustment as the maximum a posteriori likelihood problem and dev elop the RST -B A al- gorithm to solve it. In Section 5 , we use simulated data to compare RST -BA with the standard L 2 -B A algorithm, and with L 2 -B A combined with a 2 -sigma edit rule for remov- ing outliers. In Section 6 , we show the results of applying the RST -B A algorithm to the problem of reconstructing lu- nar topography from N ASA Apollo 15 orbital imagery , and compare the results with results obtained from L 2 -B A for both processed and unprocessed data. 2. Mathematical Model for B A Suppose that we hav e m images taken by a moving cam- era, or equi v alently by m cameras with different poses (lo- cations and attitudes). Suppose that multiple tie-point fea- tures are identified in each image, and that a single feature may appear in se veral images. Denote by S the control net- work of indices describing the image tie-point matching, so that ( i, j ) ∈ S if feature i was seen in image j . Let x j de- note the pose of the j th camera, and let y i denote the 3D world coordiantes of feature i . Let z ij be the pixel coor- dinates of the projection h ( x j , y i ) of feature i onto image j ,  ij be the reprojection error z ij − h ( x j , y i ) , and Σ ij be the co variance matrix associated to  ij . Define x 0 j and y 0 i to be prior estimates of camera parameters and ground control points, with covariances Ω j and Φ i , respectively . All the data e xcept for { x j } and { y i } are known and gi ven. The standard statistical model for the B A problem is as follo ws: z ij = h ( x j , y i ) +  ij  ij ∼ (0 , Σ ij ) x j ∼ ( x 0 j , Ω j ) y i ∼ ( y 0 i , Φ i ) (1) The standard L 2 approach to bundle adjustment is to assume that  ij , x j , and y i are all normally distributed, and then the maximum a posteriori likelihood solution is equiv alent to minimizing the objectiv e min 1 2 X ( i,j ) ∈ S [ z ij − h ( x j , y i )] T Σ − 1 ij [ z ij − h ( x j , y i )] + X j [ x 0 j − x j ] T Ω − 1 j [ x 0 j − x j ] (2) + X i [ y 0 i − y i ] T Φ − 1 i [ y 0 i − y i ] w .r .t. { x j } , { y i } . The prior terms for camera parameters and ground con- trol points are added to deal with gauge freedom. Ground control points can be used when 3D coordinates of certain tie-point features are well known, which is the case for Lu- nar topography data. Note that if there is no prior informa- tion on a particular x j or y i , we simply set the correspond- ing Ω − 1 j or Φ − 1 i to 0 in ( 2 ). The standard approach to bundle adjustment is to min- imize the objecti ve ( 2 ) using implicit trust region meth- ods, and in particular variants of the Levenber g-Marquardt method are very popular (see [ 18 ], [ 11 ], [ 26 ], [ 25 ], [ 6 ] for more details on these methods). F or our implementation of L 2 -B A we use a particular v ariant of the Levenber g- Marquard detailed in [ 18 , Algorithm 3.16], which is also used in the SBA implementation [ 16 ]. The method of choosing a cloud of points that ‘links’ the images together giv es rise to a sparse structure, and we exploit this structure as described in [ 11 , Algorithm A6.4]. 3. Student’ s t A pproach W e introduce the follo wing notation: for a vector u ∈ R n and any positive definite matrix M ∈ R n × n , let k u k M := √ u T M u . W e use the following generalization of the Student’ s t-distribution: p (  | µ ) = Γ( s + m 2 ) Γ( s 2 ) det[ π sR ] 1 / 2  1 + k  − µ k 2 R − 1 s  − ( s + m ) 2 (3) where µ is the mean parameter, s is the de grees of freedom, m is the dimesion of the vector  , R is a positi ve definite matrix, and √ R or R 1 / 2 denotes a Choleski factor; i.e., √ R √ R T = R 1 / 2 R T / 2 = R . A comparison of this dis- tribution with the Gaussian distribution assumed in ( 2 ) and the Laplace distribution appears in Figure 1 . Note that the Student’ s t-distribution has much thicker tails than the oth- ers, and that its influence function is redescending (see [ 19 ] for a discussion of influence functions). The main idea of the RST -B A algorithm is to assume that reprojection errors  ij come from the Student’ s t- distribution ( 3 ). W e also assume the Student’ s t-distribution prior for the initial camera parameters { x 0 j } and ground control points { y 0 i } . The intuition behind this approach is that extreme observations are much more likely in the Stu- dent’ s t model than in the Gaussian model. Therefore, a large residual will affect the ov erall fit less if fitting is done in model ( 3 ). See [ 1 ] for more details. 4. Maximum Likelihood Formulation Maximizing the likelihood for our model ( 1 ) is equi va- lent to minimizing the associated negati ve log likelihood − log p ( {  ij } ) − log p ( { x 0 j − x j } ) − log p ( { y 0 i − y i } ) Dropping the terms that do not depend on { x j } or { y i } our objectiv e is 1 2 P ( i,j ) ∈ S ( s ij + 2) log  1 + 1 s ij k z ij − h ( x j , y i ) k 2 Σ − 1 ij  + 1 2 P j ( r j + 6) log  1 + 1 r j k x 0 j − x j k 2 Ω − 1 j  + 1 2 P i ( q i + 3) log h 1 + 1 q i k y 0 i − y i k 2 Φ − 1 i i (4) T able 1. Relativ e mean µ and standard deviation σ of MSE calculated over 1000 runs for L 2 -B A, L 2 -B A with the 2 σ -edit rule (2 σ -BA), and RST -B A, presented as: µ ( σ ). Error values for world points and camera XYZ parameters respectiv ely are presented relativ e to the error incurred by L 2 -B A in the nominal case, shown in bold, i.e. where reprojection errors added were distributed as N (0 , 1) . Noise T ype W orld Points Camera XYZ L 2 -B A 2 σ -B A RST -B A L 2 -B A 2 σ -B A RST -B A N (0 , 1) 1.0 (1.3) 1.0 (1.2) 1.0 (1.3) 1.0 (0.9) 0.8 (0.7) 0.7 (1.4) . 95 N (0 , 1) + . 05 N (0 , 4) 1.3 (1.6) 1.2 (1.5) 1.1 (1.4) 6.3 (7.6) 2.7 (4.1) 3.5 (13.3) . 9 N (0 , 1) + . 1 N (0 , 4) 1.5 (1.7) 1.5 (1.9) 1.4 (1.7) 11.5 (12.5) 5.6 (7.6) 5.9 (18.1) . 95 N (0 , 1) + . 05 N (0 , 10) 2.7 (3.4) 1.8 (2.0) 1.2 (1.4) 69 (62) 23 (27) 7.3 (23) . 9 N (0 , 1) + . 1 N (0 , 10) 3.6 (4.6) 2.7 (3.0) 1.4 (1.6) 101 (76) 49 (42) 16.5 (34) . 95 N (0 , 1) + . 05 N (0 , 50) 39 (45) 21 (30) 1.9 (1.7) 580 (380) 306 (242) 12 (23) . 9 N (0 , 1) + . 1 N (0 , 50) 60 (63) 44 (47) 2.5 (2.1) 740 (510) 470 (300) 20 (36) Student’ s t, df = 4 12.3 (13.5) 12.2 (15.1) 8.9 (10.2) 240 (150) 190 (130) 38 (60) where s ij , r j , and q i are known de grees of freedom param- eters in model ( 3 ) associated to observ ations z ij , prior cam- era parameters x 0 j , and ground control points y 0 i , respec- tiv ely . The constants 2 , 6 and 3 that appear in ( 4 ) are the di- mensions of the pixel coordinates, camera poses, and world points, respectiv ely . Minimizing objectiv e ( 4 ) provides maximum a posteri- ori (MAP) likelihood estimates for parameter vectors { x j } and { y i } in the Student’ s t model ( 3 ). Now we describe an implicit trust region algorithm for minimizing ( 4 ). Giv en a sequence of column vectors { u k } and matrices { T k } we use the notation v ec( { u k } ) =      u 1 u 2 . . . u N      , diag( { T k } ) =       T 1 0 · · · 0 0 T 2 . . . . . . . . . . . . . . . 0 0 · · · 0 T N       W e define c = v ec( { x j } , { y i } ) . W e will now refer to objective ( 4 ) as F ( c ) . In order to min- imize ( 4 ), we implement an iterativ e method of the form ˆ c = c k −  H k  − 1 ∇ F ( c ) c k +1 = ( ˆ c if F ( ˆ c ) < F ( c k ) c k otherwise. (5) where k inde xes the iterations, and H k is a particular pos- itiv e definite matrix described below , which one may think of as a Hessian approximation to ∇ 2 F ( c k ) . Let J k =  A k B k  , where A k ij = ∂ x j h ( x k j , y k i ) and B k ij = ∂ y i h ( x k j , y k i ) . Define weights g k i = q i +3 q i + k y 0 i − y k i k Φ − 1 j , % k j = r j +6 r j + k x 0 j − x k j k Ω − 1 j ρ k ij = r s ij +2 s ij + k  k ij k Σ − 1 ij (6) Let ˜ A k ij = ρ k ij A k ij and ˜ B k ij = ρ k ij B k ij . Let ˜ A k and ˜ B k be matrices with block components ˜ A k ij and ˜ B k ij , respecti vely , and define ˜ J k =  ˜ A k ˜ B k  . Define H k =  ˜ J k  T Σ − 1  ˜ J k  + diag  { % k j Ω j }  + diag  { g k i Φ i }  + λ k I (7) where λ k is a regularization parameter similar to the Lev enberg-Marquardt method, and is updated according to the rule defined in Algorithm 3.16 of [ 18 ]. Specifically , λ increases quickly when the iteration ( 5 ) f ails to improv e the objectiv e function ( 4 ), and otherwise is adjusted according to the rule λ k +1 = λ k max( 1 3 , 1 − (2 φ k − 1) 3 ) (8) where φ k is the ratio of improvement predicted by the quadratic model with Hessian H k to the actual improv e- ment F k − F k +1 . From this presentation, it is clear that the RST -B A algo- rithm can be implemented by a simple reweighting of the data structures already present in L 2 -B A, and so RST -BA takes about the same time per iteration as L 2 -B A. The algo- rithm terminates when all the components of ∇ F are belo w 10 − 6 . In practice, a hard iteration limit is set, since the problems are large and it is rarely necessary to solve them exactly . W e followed this approach in testing and simula- tion. 5. Numerical Simulations The RST -B A code used for the simulated and real tests is currently implemented as part of Nasa V isionW ork- bench [ 10 ]. Since our target application is the reconstruc- tion of the lunar surface from Apollo orbital imager data, our synthetic data was modeled in a similar context. Camera positions were generated in a sequence incremented along the camera x -axis, with the z -axis of the camera coordi- nate system defined to point toward the lunar surface. The x -increment was calculated to yield the desired overlap be- tween camera fields of view , to guarantee that each point on the surface was seen by at least two cameras. Gi ven spec- ifications for the camera ele vation and location, a synthetic surface region was calculated, bounded by the camera field of view in the y -direction, the combined fields of view of the second through the penultimate cameras in the x -direction, and an estimate of minimum and maximum lunar surface height in the z -direction. 3 -dimensional world points were then randomly generated within the volume bounded by this surface. Finally , each generated 3 -dimensional point was projected into the image plane of each camera in which it was visible, gi ving a set of image coordinates for each point and camera pair . After generation of the synthetic 3D world points, we added sev eral kinds of noise to the “observ ations” made in the simulated system. In the lunar surface reconstruc- tion context, observ ations include the image coordinates of each visible point, and extrinsic camera parameters that are known up with some precision from the Apollo mission telemetry . Our data generator perturbs image coordinates, camera positions, and camera pose according to nominal Gaussian distrib utions with specified v ariance in order to simulate measurement noise and camera uncertainty . T o test the rob ustness of our algorithms against mistakes in the data, we also introduced outliers in the simulated errors ac- cording to error schemes we describe below . 1. Nominal conditions: The reprojection errors were gen- erated using the normal distribution  ij ∼ N (0 , 1) . 2. Contaminated normal: The reprojection errors were generated using a mixture of two normals, i.e.,  ij ∼ (1 − p ) N (0 , 0 . 25) + p N (0 , φ ) for v alues of p ∈ { 0 . 05 , 0 . 1 } and values of φ ∈ { 4 , 10 , 50 } . 3. Student’ s t-distrib ution: The reprojection errors were generated using Student’ s t-distribution with d f = 4 . For each run of each experiment, we ran the L 2 -B A algo- rithm as the baseline, along with L 2 -B A combined with a 2 σ -edit rule (removing ‘outliers’ that were two standards of deviation away from the mean and refitting), and the RST - B A algorithm. All degrees of freedom parameters for RST - B A were set at 4 for all of the experiments. The results for our simulated fitting are presented in T a- ble 1 . Each experiment w as performed 1000 times, and we provide the relati ve median Mean Squared Error (MSE) value and standard de viation for the difference between ‘ground truth’ and the final estimates of the algorithms, for the 3D world points data and for the camera location ( x, y , z ) data. W e left the camera pose parameters fixed at their true values during the experiment, by placing a very strong prior on them. The relativ e MSE for the world points is defined by 1 N N X k =1 k X k − ˆ X k 2 2 ! / ( MSE 0 ) (9) where N is the total number of 3D world points, X k is the k -th ‘true’ world point, ˆ X k is the estimate, and MSE 0 is the MSE of the baseline BA method in nominal conditions. The relative MSE measure for camera coordinates is simi- larly defined. The L 2 -B A method with the 2 σ -edit rule works as well or better than L 2 -B A alone. When the variance of the out- liers is very large, the 2 σ -edit rule cuts the relative error nearly in half, for both w orld points and camera parameters. The RST -BA algorithm works about as well as the 2 σ -edit rule for cases with small outliers, but as the variance of the outliers gro ws, RST -B A cuts the relati ve error by a f actor of 30 – an order of magnitude improvement o ver the 2 σ - edit rule. When the errors are actually generated from the Student’ s t-distrib ution, RST -B A cuts the camera error by a factor of 6 relativ e to the standard, and achiev es a small improv ement for the world points. 6. Bundle Adjustment in Orbital Imagery T o check the performance of the RST -BA algorithm on real data, we used imagery captured by the Apollo Met- ric Camera (AMC) on board the NASA Apollo 15 orbiter . The AMC is a frame camera with a 74 degree field of view that captured snapshots of the Moon’ s surface at re gular in- tervals. This resulted in overlap between images of 80%. The Apollo-era satellite tracking network was highly in- accurate by today’ s standards, with errors estimated to be 2.04-km for satellite station positions and 0.002 degrees for pose estimates. This creates a need for refinement via bun- dle adjustment in order to create consistent 3D models be- tween stereo pairs (see Figure 2 ). The specific frames pro- cessed were AS15-M-1089 through AS15-M-1159, which were part of Apollo 15’ s 33rd orbital re volution[ 14 ]. In order to test the ef fecti veness of L 2 -B A against RST - B A, two datasets of image measurements were created: 1. Pr ocessed Apollo tie-point data was created with ex- tensiv e processing and cleaning. First, tie points were automatically detected with the SURF [ 2 ] algorithm. Then outliers were removed using the RANSA C algo- rithm [ 5 ]. Finally the tie points were thinned do wn to 500 matches between pairs by removing the weak- est matches while ensuring that the tie points remained ev enly distributed across each image. (a) (b) (c) (d) Figure 2. Surface reconstruction from Orbit 33 images. From top to bottom: (a) L 2 -B A, processed data set; (b) RST -B A, processed data set; (c) L 2 -B A, unprocessed data set; (d) RST -BA, unpro- cessed data set. Red indicates high elev ation, blue indicates lo w elev ation. Black indicates ele vations out of the range of the color map. Ground control points were not used in these experiments. 2. Unpr ocessed Apollo tie-point data was created using the Interest Point detection algorithm based on SIFT [ 17 ]. No outlier rejection was done in this case, yield- ing a data set with up to 50% outliers. L 2 -B A and RST -BA were tested with both processed and unprocessed data sets. Results of these tests are shown in T able 2 . Here, triangulation error is a measure of the av er- age distance between the closest point of intersection of tw o forward projected rays for a set of tie-points. A decrease in triangulation error indicates a substantial improvement in the self-consistency of the DEMs in the data set. After bundle adjustment was complete, we processed the imagery using stereo reconstruction tools [ 24 ] to produce dense topography of the lunar surface. This 3D reconstruc- tion used the improved camera extrinsic parameters from bundle adjustment to produce a more consistent, seamless mosaic of 3D topographic models. Figure 2 sho ws these results with v arious bundle adjusment tests. T opography is represented by a color-map with red indicating high ele va- tion and blue low ele vation. The original (unadjusted) camera parameters show clear discontinuities between adjacent models that are due to the uncertainties in the original Apollo tracking data. While Pr ocessed data yields reasonable results when ground con- trol points are used, we did not use these points to empha- size that RST -B A can be used in their absence while L 2 - B A cannot. W ithout ground control points, L 2 -B A found a ‘kink’ in the DEM, which is responsible for the black sec- tions visible in Figure 2 while RST -B A was able to recon- struct the DEM. The results from the unpr ocessed data set, which con- tained nearly 50% outliers, sho w a stark dif ference between the two approaches. The L 2 -B A algorithm failed to cre- ate any improvement. Instead, the outliers caused sev ere and unpredictable distorition of the camera parameters. The RST -B A algorithm, on the other hand, was nearly unaf- fected, and produced results remarkably similar to those produced from the pr ocessed data set. T able 2 shows that median triangulation error were only slightly higher for the unpr ocessed data than they were for the pr ocessed data. This data suggests that RST -B A is significantly more robust to outliers than the standard bundle adjustment technique. 7. Conclusion W e have proposed RST -B A, a robust bundle adjustment algorithm, based on the Student’ s t distribution, for per- forming bundle adjustment in the presence of outliers in tie-point matching. RST -BA preserves the sparse structure, and hence the speed, of L 2 -B A, and can be implemented by simple modifications to the L 2 -B A algorithm. Our test re- sults on both synthetic and real data show that when the data hav e been preprocessed to remove outliers in the tie-point matches, RST -B A outperforms L 2 -B A by a small margin, and on unpreprocessed data with numerous outliers, RST - B A outperforms L 2 -B A by a significant margin. RST -BA demonstrates significant adv antages in both speed and ac- curacy of results ov er both L 2 -B A and L 2 -B A with a 2 σ - edit rule, and can be used to reconstruct DEMs without data preprocessing and without ground control points. In future work, we will perform an extended comparison of RST -B A with “robust” methods such as Cauchy re-weighting in ad- dition to the 2 σ -edit rule. W e will also work on the problem of estimating the degrees of freedom parameters, which are currently assumed to be known by the RST -B A algorithm. 8. Acknowledgements W e would like to thank our colleagues at the Arizona State University for supplying high resolution scans of the Apollo Metric Camera images. This work was funded by the N ASA Lunar Adv anced Science and Exploration Research (LASER) program grant #07-LASER07-0148, N ASA Advanced Information Systems Research (AISR) program grant #06-AISRP06-0142, and by the NASA ESMD Lunar Mapping and Modeling Program (LMMP). W e would also like to thank James Burke and Bradley Bell at the Univ ersity of W ashington for guidance in dev eloping the optimization algorithm. References [1] A. Aravkin. Robust Methods with Applications to Kalman Smoothing and Bundle Adjustment . PhD thesis, Univ ersity of W ashington, Seattle, W A, June 2010. 3 [2] H. Bay , A. Ess, T . T uyetelaars, and L. V angool. Speeded-up robust features (surf). Computer V ision and Image Under- standing , 110(3):346, 2008. 5 [3] S. S. Brandt and U. Ziese. Robust alignment of transmission electron microscope tilt series. P attern Recognition, Interna- tional Confer ence on , 4:683–686, 2006. 2 [4] A. Cumani and A. Guiducci. V isual odometry for robust rov er na vigation by binocular stereo. In Pr oc. WSEAS IC- SSIP , 2006. 2 [5] M. A. Fischler and R. C. Bolles. Random sample consen- sus: a paradigm for model fitting with applications to image analysis and automated cartography . Communications of the A CM , 24(6):381, 1981. 5 [6] Fletcher . Practical Methods of Optimization . Joh Wile y and Sons, second edition, 1987. 3 [7] P . Fua. Regularized bundle-adjustment to model heads from image sequences without calibration data. Int. J. Comput. V ision , 38(2):153–171, 2000. 2 [8] P . Fua, R. Plaenkers, and D. Thalmann. From synthesis to analysis: Fitting human animation models to image data. Computer Graphics International Confer ence , 0:4, 1999. 2 [9] N. Fukaya, T . Anai, H. Sato, N. Kochi, M. Y amada, and H. Otani. Application of rob ust regression for e xterior orien- tation of video images. In ISPRS08 , 2008. 2 [10] M. Hancher, M. Broxton, T . Fong, and Contributors. N ASA V isionW orkbench . https://github .com/visionworkbench/visionworkbench . 4 [11] R. I. Hartley and A. Zisserman. Multiple V ie w Geometry in Computer V ision . Cambridge Uni versity Press, ISBN: 0521540518, second edition, 2004. 2 , 3 [12] R. I. Hartley and A. Zisserman. Multiple V ie w Geometry in Computer V ision . Cambridge Uni versity Press, ISBN: 0521540518, second edition, 2004. 2 [13] Q. K e and T . Kanade. Quasiconv ex optimization for robust geometric reconstruction. In IEEE International Confer ence on Computer V ision (ICCV 2005) , volume 2, pages 986–993, October 2005. 2 [14] S. Lawrence, M. Robinson, M. Broxton, J. Stopar , W . Close, J. Grunsfeld, R. Ingram, L. Jefferson, S. Locke, R. Mitchell, T . Scarsella, M. White, M. Hager , T . W atters, E. Bo wman- Cisneros, J. Danton, and J. Garvin. The apollo digital image archiv e: New research and data products. In Pr oceedings of the NLSI Lunar Science Confer ence , page 2066, 2008. 5 [15] T . Lee. Robust 3d street-view reconstruction using sky motion estimation. In Pr oceedings of the IEEE Interna- tional W orkshop on 3-D Digital Imaging and Modeling (3DIM2009) in conjunction with ICCV , October 2009. 2 [16] M. Lourakis and A. Ar gyros. Sba: A software package for sparse bundle adjustment. ACM T ransactions on Mathemat- ical Softwar e , 36(1):30, 2009. 2 , 3 [17] D. G. Lowe. Distincti ve image features from scale-inv ariant keypoints. International Journal of Computer V ision , 60(2):91, 2004. 6 [18] K. Madsen, H. Nielsen, and O. Tinglef f. Methods for Non- Linear Least Squar es Pr oblems . T echnical University of Denmark. Springer , 2nd edition, 2004. 2 , 3 , 4 [19] R. A. Maronna, D. Martin, and Y ohai. Robust Statistics . W i- ley Series in Probability and Statistics. John W iley and Sons, 2006. 2 , 3 [20] H. Mayer . Robust least-squares adjustment-based orientation and auto-calibration of wide-baseline image sequences. In O. Hellwich, I. Niini, C. Ressl, V . Rodehorst, D. Scharstein, and P . Sturm, editors, T owards Benchmarking Automated Calibration, Orientation and Surface Reconstruction fr om Images , Beijing, China, 2005. ICCV 2005. 2 [21] H. Mayer . 3d reconstruction and visualization of urban scenes from uncalibrated wide-baseline image sequences (mayer) 2006. In IAPRS V olume XXXVI, P art 5 , Dresden, 2006. 2 [22] H. Mayer , M. Mosch, and J. Peipe. Comparison of pho- togrammetric and computer vision techniques - 3d recon- struction and visualization of wartbur g castle. In CIP A 2003 19th International Symposium , pages 103–108, An- talya, 2003. 2 [23] E. Mouragnon, M. Lhuillier , M. Dhome, F . Deke yser , and P . Sayd. Generic and real-time structure from motion us- ing local bundle adjustment. Image and V ision Computing , 27(8):1178–1193, 2009. 2 [24] A. Nefian, K. Husmann, M. Broxton, V . T o, , M. Lundy , and M. Hancher . A bayesian formulation of sub-pixel refinement in stereo orbital imagery . In Int’l Conf on Image Pr ocessing , 2009. 6 [25] J. Nocedal and S. J. Wright. Numerical Optimization . Springer Series in Operations Research. Springer , 1999. 3 Dataset Algorithm Start T riangulation Error End T riangulation Error Min Median Max Min Median Max Processed L 2 -B A 0.0616 714.18 134234 0.0008 36.309 178636 RST -B A 0.0616 714.18 134234 0.0002 14.329 132966 Unprocessed L 2 -B A 0 729.87 239452 0 4409.9 242179 RST -B A 0 729.87 239452 0 14.404 240647 T able 2. Triangulation errors for the Apollo DEM reconstruction, using pr ocessed and unpr ocessed Apollo tie-point data. RST -BA performs much better than L 2 -B A for unprocessed data, and better than L 2 -B A for processed data. [26] G. Seber and C. W ild. Nonlinear Re gr ession . John W iley and Sons, 1988. 3 [27] N. S ¨ underhauf, K. K onolige, S. Lacroix, and P . Protzel. V i- sual odometry using sparse b undle adjustment on an au- tonomous outdoor vehicle. In P . Levi, M. Schanz, R. Lafrenz, and V . A vrutin, editors, AMS , Informatik Aktuell, pages 157– 163. Springer , 2005. 2 [28] B. Triggs, P . F . McLauchlan, R. I. Hartle y , and A. W . Fitzgib- bon. Bundle adjustment - a modern synthesis. In ICCV ’99: Pr oceedings of the International W orkshop on V ision Algo- rithms , pages 298–372, London, UK, 2000. Springer -V erlag. 1 , 2 [29] J. Zhang, M. Boutin, and D. G. Aliaga. Robust b undle ad- justment for structure from motion. In Imag e Pr ocessing, 2006 IEEE International Conference on , pages 2185–2188, Oct. 2006. 2