Back to main page.

How to do math in the Ubitrack library

Working with boost::uBLAS

Useful links:

• uBLAS operations overview - definitely worth reading to get an overview and to look up function names.

uBLAS lives in the namespace boost::numeric::ublas, which can be rather inconvenient. We suggest to use a shortcut in the cpp file (not in a header! never!) like this:

    namespace ublas = boost::numeric::ublas;
    ublas::vector< double > v = ublas::prod( m, p );

ublas provides classes matrix and vector, but if you know the size of your matrix/vector in advance we recommend to use Ubitrack::Math::Matrix< M, N, T=double > and Ubitrack::Math::Vector< N, T=double > instead. Rationale: Clearer interfaces, slightly higher performance and full Lapack compatibility (see below).

Working with Lapack

To solve exciting numerical problems in Ubitrack, we use the Lapack library in combination with the boost numerical bindings to remove some complexity of working with lapack.

Documentation:

• Lapack User's Guide - to find out which function you need.
• /3rd/boost-bindings/ - to find out if the function is provided by the boost numerical bindings. You will have to implement the wrapper yourself otherwise...
• Fortran source code - look up the comments in the Fortran sources if you are unsure what the arguments mean.

A shortcut to the lapack bindings namespace is provided by:

    namespace lapack = boost::numeric::bindings::lapack;

As Lapack is written in Fortran, all matrices have to be column major, Fortran's default storage format. Math::Matrix is already column major, but if you use ublas::matrix you will have to write something like this:

    ublas::matrix< double, ublas::column_major > A( m, n );

Solutions to common numerical problems

Finding the right singular vector corresponding to the smallest singular value

E.g. in DLT estimation of homographies or projection matrices.

    ublas::matrix< double, ublas::column_major > biqEqSystem( 200, 9 ) = ...;
    Math::Matrix< 9, 9 > Vt;
    Math::Vector< 9 > s;
    ublas::matrix< double, ublas::column_major > U; // dummy, not used
    lapack::gesvd( 'N', 'A', bigEqSystem, s, U, Vt );

The desired vector is found in the last row of the matrix Vt: ublas::row( Vt, 8 )

Note that gesvd will destroy the contents of the matrix bigEqSystem!

Solving symmetric equation systems

E.g. in optimizers where the matrix looks like B^T B.

    Math::Matrix< 9, 9 > A( ublas::prod( ublas::trans( B ), B ) );
    Math::Vector< 9 > b = ...;
    lapack::posv( 'U', A, b );

The vector b is overwritten with the result.

Working with Blas

You rarely need to use blas directly, as all blas functionality is also provided by boost::ublas. The only advantage is the probably higher performance for large matrices. Blas bindings are also provided by boost, similar to the lapack bindings:

    namespace blas = boost::numeric::bindings::blas;

Various tricks

Functions that work on floats OR doubles

When implementing tracking or calibration algorithms, one frequently implements algorithms where it makes sense to have both float (fast) and double (precise) versions. In this case, provide two overloaded versions in the header file and use one template implementation in the cpp file. Example:

    /* my.h */
    UBITRACK_EXPORT Math::Matrix< 3, 3, float > myFunction( const std::vector< float >& );
    UBITRACK_EXPORT Math::Matrix< 3, 3, double > myFunction( const std::vector< double >& );

    /* my.cpp */
    template< typename T >
    Math::Matrix< 3, 3, T > myFunctionImpl( const std::vector< T >& )
    { ... }

    Math::Matrix< 3, 3, float > myFunction( const std::vector< float >& p )
    { return myFunctionImpl( p ); }
    Math::Matrix< 3, 3, double > myFunction( const std::vector< double >& p )
    { return myFunctionImpl( p ); }

Rationale: keeps compilation times of client applications down and reduces dependencies of clients to other libraries (e.g. Lapack).

Functions that return sub-matrices

If you implement a function that creates a matrix that is meant to be used as a sub-matrix, e.g. a Jacobian or part of a larger equation system, you can write the function as a template like this:

    template< class MA >
    void createMyJacobian( ublas::matrix_expression< MA >&  J, ... )
    { 
        typename MA::value_type x = ...;
        J( 0, 0 ) = x; 
    }

To assign, e.g. the first row of a matrix, the function can be called using

    ublas::matrix< float > myMatrix( 30, 52 );
    createMyJacobian( ublas::subrange( myMatrix, 0, 1, 0, 52 ), some_param );