sba FAQ

sba is a C/C++ package for generic sparse bundle adjustment that is distributed under the GNU General Public License. sba is generic in the sense that it provides increased flexibility regarding the choice of parameters and functional relations defining the image projections involved in bundle adjustment. For more details, refer to the related paper published in ACM TOMS.

Assume that we are given initial estimates of the 3D coordinates corresponding to a set of points seen in a series of images, as well as initial estimates for the viewing parameters pertaining to each image. Bundle Adjustment (BA) is a large optimization problem involving the simultaneous refinement of the 3D structure and viewing parameters (i.e. camera pose and possibly intrinsic calibration and radial distortion), aiming to obtain a reconstruction which is optimal under certain assumptions regarding the noise pertaining to the observed image features: If the image error is zero-mean Gaussian, then BA is the Maximum Likelihood Estimator. Its name refers to the “bundles” of light rays originating from each 3D feature and converging on each camera's optical center, which are adjusted optimally with respect to both the structure and viewing parameters. sba uses a custom implementation of the Levenberg-Marquardt nonlinear least squares algorithm to solve the sparse, large-scale optimization problem associated with BA.

Because of the lack of interaction among parameters for different 3D points and cameras, the linear systems (i.e. normal equations) that have to be solved during BA contain lots of zeroes and are characterized by a sparse block structure. This structure can be exploited by avoiding storing and operating on zero elements, thus gaining considerable computational benefits. Click here to see a typical example of this sparse block structure. By the way, these data have originated from Oxford's corridor sequence.

BA involves the solution of a large-scale minimization problem, typically involving several thousand variables. Optimization algorithms work by iteratively linearizing the function to be minimized in the neighborhood of the current estimate and obtaining linear systems whose solution determines an increment to the current estimate. Most optimization codes (e.g. minpack, nl2sol/n2g, etc) assume that those linear systems are dense, i.e. contain very few zeroes. The computational complexity associated with solving such systems is O(N**3) in the number of unknowns. Recalling that BA involves the solution of a large-scale minimization problem (typically involving several thousand variables), it is clear that most general-purpose optimization codes would incur prohibitive processing and memory storage overheads when applied to BA.

An excellent overview of the application of BA to vision-based reconstruction is given by Triggs et al. in Bundle Adjustment: A Modern Synthesis, Vision Algorithms: Theory and Practice, LNCS Vol. 1883, pp. 298-372, Springer-Verlag, 2000. For a much shorter exposition, you can refer to my bundle adjustment article in Wikipedia.

Consult the lecture notes entitled Methods for Non-Linear Least Squares Problems, by K. Madsen, H.B. Nielsen and O. Tingleff, Technical University of Denmark, 2004.

sba gives to its user full control over the definition of the parameters describing cameras and 3D structure. Therefore, it can support a wide range of instances of the multiple view reconstruction problem such as projective reconstruction, Euclidean reconstruction, intrinsic camera parameters refinement, etc. sba also provides routines for dealing with intersectioning/resectioning problems, in which the camera pose and scene structure are respectively held constant.

sba can be invaluable to researchers and practitioners in the fields of computational vision, robotics, image-based graphics, photogrammetry, surveying, cartography, etc. It is being used in several labs around the globe and is currently the only such software worldwide to be offered in source under the GNU GPL license with the option of commercial licensing. If you want to know more about how is sba being used, click here for a sample list of published papers employing it.

To be able to use all sba functions, you should have installed LAPACK or an equivalent vendor library. Check http://www.netlib.org/clapack for an f2c'ed free version. It has been reported that the precompiled MSWin libraries at that site are broken and should be rebuilt using the included project files. You might also try these precompiled MSWin LAPACK/BLAS libraries.

First, you have to make sure that LAPACK is installed at your site. If you have to install it, please follow the installation notes included with your LAPACK distribution. The second step is to compile sba itself. The tarfile contains makefiles for Unix/Linux using gcc (Makefile) and MSWin using Visual Studio (Makefile.vc). Please read the comments within those files for more information. Based on the supplied Makefiles, it should be straightforward to compile sba using any ANSI-compliant C compiler.

When building eucsbademo, certain LAPACK/BLAS libraries cause the linker to fail with an error like undefined reference to `MAIN__', or similar. This is due to these libraries calling a Fortran main() routine for initialization. There are several workarounds to this problem:

• Under gcc, pass the option -u MAIN__ to the linker.
• Define a dummy MAIN__ function as int MAIN__(int argc, char **argv) { abort(); return 1; }
• Under Windows, use these precompiled libraries.

See also the discussion here.

It seems that you have an old version of LAPACK. dgesdd was introduced in LAPACK 3.0. If you do not wish to upgrade, an alternative is to use the somewhat slower dgesvd. sba_lapack.c has the equivalent calls to dgesvd within comments.

A detailed description of the various routines included in sba v.1.5 can be found in M.I.A. Lourakis and A.A. Argyros, 2009. SBA: A software package for generic sparse bundle adjustment. ACM TOMS 36(1), pp. 1-30. Be warned that versions after 1.5 might differ slightly. There's also an older technical report, please refer to ICS/FORTH TR-340: The Design and Implementation of a Generic Sparse Bundle Adjustment Software Package Based on the Levenberg-Marquardt Algorithm, by M.I.A. Lourakis and A.A. Argyros, 2004.
Detailed comments explaining the role of each argument are also included in the source code of each function.

As mentioned above, sba is quite flexible in the choice of parameters and functional relations describing cameras, 3D structure and image projections. Therefore, it can support several flavors of BA involving arbitrary projective or affine cameras, partially or fully intrinsically calibrated cameras, exterior orientation (i.e. pose) estimation from fixed 3D points, 3D reconstruction from extrinsically calibrated images, refinement of intrinsic calibration parameters, etc. To adapt sba to a particular problem, its user has to choose an appropriate parameterization for cameras and 3D structure and then supply code implementing the corresponding projection function and, optionally, its Jacobian. The demo program included with sba illustrates this in more detail dealing with the problem of Euclidean BA.

Since BA involves the solution of a minimization problem that is solved iteratively, initial estimates of the various parameters should be supplied to sba. Such estimates define an initial 3D reconstruction and can be computed with any structure and motion estimation vision algorithm such as the ones described in Hartley and Zisserman's textbook.

Short answer: No.
Longer answer: Assuming that the localization error of image point features is zero-mean Gaussian, BA is the Maximum Likelihood Estimator. It is exactly this property that makes BA so powerful for simultaneously refining the 3D structure and viewing parameters. Note, however, that the previous assumption does not include cases where gross outliers (i.e. mismatched features) are present. Such outlying features should be detected and eliminated using geometry-based techniques before BA. More details on such techniques can, for example, be found here.

By passing NULL as the value of the argument corresponding to the Jacobian, the former is computed numerically with the aid of forward finite differences. Note, however, that this choice is not recommended: For achieving maximum efficiency, you are strongly advised to supply a function evaluating the Jacobian analytically.

Assume that you want to pass two arrays, one of doubles and one of integers. A quick but not very elegant solution is to use global variables. The simplest way to avoid using globals is to declare a structure as

struct mydata{
   double dar[XXX];
   int iar[YYY];
};

where XXX and YYY denote the appropriate array sizes. Then, define a structure variable with

struct mydata data;

and fill it:

data.dar[0]=7.0;
data.iar[0]=-17;
// etc

Following this, call the appropriate sba routine passing it the address of data as the adata argument, e.g.

ret=sba_motstr_levmar(..., proj, projac, (void *)&data, itmax, verbose, opts, info);

Your proj and projac routines should interpret the supplied data using type casting:

struct mydata *dptr;

dptr=(struct mydata *)adata; // adata is passed as void *

// supplied data can now be accessed as dptr->dar[0], etc
 

File sba_chkjac.c includes routines for checking the correctness (i.e. consistency with the user-supplied projection functions) of the Jacobians passed to the expert BA functions. These routines can be invoked indirectly by specifying 0 as the value of itmax in all sba_XXX_levmar() and sba_XXX_levmar_x() routines, which return immediately with a value of 0 in this case. Verification routines print suspicious gradients (i.e. Jacobian rows) to stderr. It should be noted, however, that these routines are not 100% foolproof since depending on the point of evaluation, they might complain even for correct Jacobians. A solution to this is to try to check the correctness of the Jacobians evaluated at some other point(s). Refer to the comments in sba_chkjac.c for more details. Also, keep in mind that these routines are meant to verify analytic Jacobians rather than finite difference numerical approximations.

The full message looks like this:

LAPACK error: the leading minor of order XX is not positive definite,
the factorization could not be completed for dpotf2/dpotrf in sba_Axb_Chol()

This means that the Cholesky-based solver for the augmented normal equations fails at some iteration. In theory, the matrix of the system being solved should be positive definite, thus Cholesky should always work. Due to numerical error, however, sometimes this is not true and the Cholesky factorization cannot be computed. Nevertheless, no harm is done since when an iteration fails, the damping of the normal equations is increased so in a future iteration their matrix will be positive definite. In other words, this error does not affect convergence. To prevent sba from complaining, comment the calls to sba_Axb_Chol() in sba_levmar.c and uncomment the ones to sba_Axb_LU(). These changes will result in using the LU decomposition in place of Cholesky for solving the normal equations and will get rid of the warning. Note, however, that LU decomposition takes longer to compute, so the modifications will render sba a bit slower.

This error return occurs when sba's routines detect one or more invalid (i.e., NaN or Inf) values in the predicted projections generated by the user-supplied projection function. If you are using your own projection functions, make sure that they have been coded correctly. For the projection functions included with sba's demo program, this is typically due to invalid initial estimates for certain 3D points; however, it can also be attributed to erroneous initial camera parameters. To have the offending projections printed on the screen, re-run eucsbademo with a verbosity level 2 or higer. This can be achieved by setting the verbose variable in eucsbademo.c to 2 and recompiling.

Yes. In some applications, the parameters of certain 3D points are precisely known and should not be refined during BA. Such points, for example, are the so-called ground control points used in satellite imagery. Starting with version 1.6, sba supports the use of fixed points through argument ncon.

Yes. sba supports the use of cameras whose parameters are not changed during BA through argument mcon.

As explained in Q7, sba is very flexible with respect to the parametrization chosen for bundle adjustment. To demonstrate its use for Euclidean bundle adjustment, the demo directory contains the eucsbademo program. Up to version 1.3, eucsbademo supported only the case of camera intrinsics that are indentical for all images and remain fixed throughout bundle adjustment. Starting in version 1.3, eucsbademo also supports the case of camera intrinsics that are different among images and are allowed to vary during bundle adjustment. In addition to the above, support for varying distortion parameters was added in eucsbademo version 1.6. Note that in order to have the distortion parameters refined via BA, the intrinsics should vary and be refined as well.

Yes. The eucsbademo program supports BA with varying camera intrinsics, different for each camera (see also Q24). eucsbademo supports different numbers of intrinsic parameters that are to be kept fixed during BA, as follows: fixed skew, fixed aspect ratio and skew, fixed aspect ratio, skew and principal point. More details can be found in the first lines of eucsbademo.c.

Yes. Similarly to the number of fixed intrinsics (see Q25), eucsbademo supports different numbers of distortion parameters that are to remain unchanged during BA. More details can be found near the top of eucsbademo.c.

Normaly the demo program produces no output; only some statistics regarding the minimization are printed. To have the motion and/or structure parameters printed to stdout after the BA, edit function sba_driver() in eucsbademo.c and uncomment one of the assignments setting the prnt variable to one of BA_MOTSTRUCT, BA_MOT or BA_STRUCT.

The demo program assumes a right-handed camera coordinate system with the Z axis aligned with the camera optical axis and pointing towards the scene. The positive X axis points right and the positive Y axis points down. If you are familiar with OpenGL, the Y and Z axes point in opposite directions compared to OpenGL. Another right-handed system is assumed for the world coordinate system. Please examine this, this and the following schematic for graphical illustrations.

   Z /
    /
   /     X
  +------->
  |
  |
Y |
  v

The covariance matrix of an image point can be only approximately estimated; the most common practice is to employ the spatial derivatives of image intensity. More details are, for example, given by Brooks et al. in What value covariance information in estimating vision parameters?, ICCV01, vol. I, pp. 302-308, Vancouver, IEEE Press, 2001. In the case of multiscale features, the work of Zeisl et al. Estimation of Location Uncertainty for Scale Invariant Feature Points, BMVC 2009 [source code here] can be employed.

If your application uses 3x3 rotation matrices to represent rotations, they can be converted to quaternions (and vice-versa) suitable for sba using the Perl utility functions included in quats.pl in the utils directory.

The eucsbademo demo program included with sba represents rotation using quaternions. More details regarding quaternions and rotation can be found in Berthold Horn's notes here.

Yes. Julien Seinturier has written JSBA, which is an sba wrapper for Java. This allows sba to be used within pure Java code without the need to provide native libraries. Please refer to JSBA's web page for more details and code.

Yes. Starting with version 1.4, sba includes a MEX-file that can be used to interface it with matlab. See sba.c in the matlab directory of sba's source distribution. Documentation for using sba with matlab is provided in file sba.m. Be warned that using sba from matlab is much slower compared to directly calling it from your C or C++ code! Assuming that the MEX-file has been compiled (see also Q34), file eucsbademo.m contains examples of its use with matlab. Try it with the command matlab < eucsbademo.m. The MEX-file might also work with GNU octave, although this has not been tested; if you succeed with it, please let us know.

The answer depends on your operating system/compiler. The Makefile in the matlab directory works with GCC under Linux, whereas Makefile.w32 works with MSVC under Windows. You might also need to run "mex -setup" from the matlab or OS command prompt to configure your system to build MEX-files. In any case, the provided makefiles can be used as a starting point for compiling on your own system.

The answer depends on the choice made for implementing the user-supplied functions for image projections and their Jacobians. If the user chooses to code them in matlab, then the matlab MEX interface has to make repeated calls to matlab to evaluate them through the mexCallMATLAB() function. This is inherently inefficient, resulting in slow overall execution. Performance can be improved significantly by coding the image projection and Jacobian functions in a language such as C and packaging them in a shared (dynamic in Windows terminology) library. The MEX interface can then be instructed to load such functions from the library using a funcname@libname format, please refer to sba.m and the example in projac.c for more details.

sba is distributed under GPL with the intention of allowing its free use for research and personal purposes. GPL allows you to modify the original code into a new program and then sell copies of the modified program commercially, but only under the terms of the GPL. Thus, you must make the source code available to the users of the program as described in the GPL, and they must be allowed to redistribute and modify it as described in the GPL. Clearly, such requirements are in conflict with using sba as part of a proprietary application. However, it is possible to obtain a paid license for proprietary commercial use under terms different than those of GPL. Please contact the author for further information on such licensing. More answers to questions related to the GPL license can be found here.