eightPointAlgorithm.cpp

Go to the documentation of this file.

#include <dlrComputerVision/eightPointAlgorithm.h>

namespace num = dlr::numeric;
namespace linAlg = dlr::linearAlgebra;

namespace dlr {

  namespace computerVision {


    // This function is used internally by eightPointAlgorithm() to
    // translate and scale input points so that their mean lies at the
    // origin and they have isotropic unit variance.
    void
    normalizePointSequence(dlr::numeric::Array2D<double> const& inputPoints,
                           dlr::numeric::Array2D<double>& outputPoints,
                           dlr::numeric::Array2D<double>& transform)
    {
      // Following Hartley (and quoting almost directly from the
      // paper), let u_i = [u_i,x , u_i,y , 1.0]^T, and form the
      // matrix
      //
      //   E = sum(u_i * u_i^T).
      //
      // Since this matrix is symmetric and positive definite, we may
      // take its Choleski factorization
      //
      //   E = K * K^T
      // 
      // where K is upper triangular. It follows that
      //
      //   sum(inv(K / sqrt(N)) * u_i * u_i^T * inv(K / sqrt(N))^T) = NI,
      //
      // where I is the identity matrix, [and N is the number of
      // elements in the sum]. [...] Note that K is upper triangular,
      // and so it represents an affine transformation.
      //
      // In other words transforming by K/sqrt(N) puts the centroid of
      // our points at the origin, and normalizes their variance.
      // Just what we need.
      dlr::numeric::Array2D<double> EE =
        num::matrixMultiply(inputPoints.transpose(), inputPoints);
      double sqrtN = std::sqrt(inputPoints.rows());

      num::Array2D<double> KK;
      linAlg::choleskyFactorization(EE, KK, false);
      KK /= sqrtN;
      num::Array2D<double> KKInv = linAlg::inverse(KK);

      // OK, now we've found our normalizing transform.  Apply it
      // to the input points.
      outputPoints = num::matrixMultiply(inputPoints, KKInv.transpose());
      transform = KKInv;
    }

  } // namespace computerVision
    
} // namespace dlr