featureAssociation.h

Go to the documentation of this file.

#ifndef DLR_COMPUTERVISION_FEATUREASSOCIATION_H
#define DLR_COMPUTERVISION_FEATUREASSOCIATION_H

#include <dlrNumeric/array2D.h>

namespace dlr {

  namespace computerVision {

    template<class Iterator0, class Iterator1, class Functor>
    std::vector< std::pair<size_t, size_t> >
    associateFeaturesScott91(Iterator0 sequence0Begin, Iterator0 sequence0End,
                             Iterator1 sequence1Begin, Iterator1 sequence1End,
                             Functor similarityFunctor);
    
  } // namespace computerVision
    
} // namespace dlr


/* ============ Definitions of inline & template functions ============ */


#include <cmath>
#include <dlrLinearAlgebra/linearAlgebra.h>
#include <dlrNumeric/maxRecorder.h>
#include <dlrNumeric/utilities.h>


namespace dlr {

  namespace computerVision {


    // This function template implements the feature association
    // algorithm of Guy Scott and H. Christopher Longuet-Higgins.
    template<class Iterator0, class Iterator1, class Functor>
    std::vector< std::pair<size_t, size_t> >
    associateFeaturesScott91(Iterator0 sequence0Begin, Iterator0 sequence0End,
                             Iterator1 sequence1Begin, Iterator1 sequence1End,
                             Functor similarityFunctor)
    {
      // Count the number of features in each sequence.  Even if this
      // operation is O(n), it will be dominated by the more expensive
      // steps below.
      size_t sequence0Length = sequence0End - sequence0Begin;
      size_t sequence1Length = sequence1End - sequence1Begin;

      // Compute the similarity matrix.
      dlr::numeric::Array2D<double> GMatrix(sequence0Length, sequence1Length);
      Iterator0 begin0 = sequence0Begin;
      for(size_t rr = 0; rr < sequence0Length; ++rr) {
        Iterator1 begin1 = sequence1Begin;
        for(size_t cc = 0; cc < sequence1Length; ++cc) {
          GMatrix(rr, cc) = similarityFunctor(*begin0, *begin1);
          ++begin1;
        }
        ++begin0;
      }

      // Compute a new similarity matrix, comprising only rotations
      // and reflections, that is "similar" to G (see the paper for
      // more details).
      dlr::numeric::Array2D<double> uMatrix;
      dlr::numeric::Array1D<double> sigmaArray;
      dlr::numeric::Array2D<double> vTransposeMatrix;
      dlr::linearAlgebra::singularValueDecomposition(
        GMatrix, uMatrix, sigmaArray, vTransposeMatrix);
      dlr::numeric::Array2D<double> PMatrix =
        dlr::numeric::matrixMultiply(uMatrix, vTransposeMatrix);

      // Find the max elements for each row and column of P.
      std::vector< dlr::numeric::MaxRecorder<double, size_t> > rowMaxes(
        sequence0Length);
      std::vector< dlr::numeric::MaxRecorder<double, size_t> > columnMaxes(
        sequence1Length);
      for(size_t rr = 0; rr < sequence0Length; ++rr) {
        for(size_t cc = 0; cc < sequence1Length; ++cc) {
          double similarity = PMatrix(rr, cc);
          rowMaxes[rr].test(similarity, cc);
          columnMaxes[cc].test(similarity, rr);
        }
      }

      // Valid correspondences are those for which the similarity is
      // max of both row and column.
      std::vector< std::pair<size_t, size_t> > result;
      for(size_t rr = 0; rr < sequence0Length; ++rr) {
        size_t bestColumn = rowMaxes[rr].getPayload();
        if(columnMaxes[bestColumn].getPayload() == rr) {
          result.push_back(std::make_pair(rr, bestColumn));
        }
      }

      return result;
    }

  } // namespace computerVision
    
} // namespace dlr

#endif /* #ifndef DLR_COMPUTERVISION_FEATUREASSOCIATION_H */