HMatrix.hh

/* 
 This file is part of the BIAS library (Basic ImageAlgorithmS).

 Copyright (C) 2003-2009    (see file CONTACT for details)
 Multimediale Systeme der Informationsverarbeitung
 Institut fuer Informatik
 Christian-Albrechts-Universitaet Kiel


 BIAS is free software; you can redistribute it and/or modify
 it under the terms of the GNU Lesser General Public License as published by
 the Free Software Foundation; either version 2.1 of the License, or
 (at your option) any later version.

 BIAS is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 GNU Lesser General Public License for more details.

 You should have received a copy of the GNU Lesser General Public License
 along with BIAS; if not, write to the Free Software
 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 */

#ifndef __HMatrix_hh__
#define __HMatrix_hh__
#include "bias_config.h" 

#include <Base/Geometry/HMatrixBase.hh>
#include <Base/Image/Image.hh>
#include <Geometry/PMatrix.hh>

namespace BIAS {

/** @class HMatrix
 @ingroup g_geometry
 @brief a 3x3 Matrix describing projective transformations between planes
 */
class BIASGeometry_EXPORT HMatrix : public HMatrixBase
{
public:
  HMatrix();

  explicit HMatrix(const MatrixInitType& i)
  : HMatrixBase(i) {};

  HMatrix(Matrix3x3<double> Mat);

  /** @brief homography induced by a plane (e.g. infinite homography) 
   @author koeser 06/2008 */
  inline HMatrix(const PMatrix& P1, const PMatrix& P2,
      const HomgPlane3D& scenePlane=HomgPlane3D(0,0,0,1)) {
    MapAcrossPlane(P1,P2,scenePlane);
  }

  ~HMatrix();

  HMatrix& operator=(const HMatrix& mat);

  /** maps the rectangle (0,0),(OrigWidth,OrigHeight) with this H
   and computes mapped range as corner + size. JW 07/2003 */
  void GetMappedRange(const unsigned int OrigWidth,
      const unsigned int OrigHeight,
      BIAS::Vector2<double> & corner,
      double & sizeWidth, double & sizeHeight) const;

  /** @return the squared distance between the mapped point from
   the first picture and the measured point in the second picture
   @author mdunad 02 2004 */
  double GetMappedError(const BIAS::HomgPoint2D &PointPicOne,
      const BIAS::HomgPoint2D &PointPicTwo) const;

  /** @brief returns jacobian of H: R^2 --> R^2
   
   H as a mapping of 2D to 2D is not a linear mapping, because of
   the nonlinear homogenization. This function returns the 2x2 Jacobian
   matrix at the point x. You have to pass homogenized points and a
   correctly sized jacobian

   @param x homogenized point (real point in R^2) for the jacobian
   @param Jac (output) 2x2 jacobian at x
   @author koeser 04/2005 */
  void GetJacobian(const HomgPoint2D& x, Matrix<double>& Jac) const;

  /** @brief returns 1st order Taylor expansion of homography at x 
   @author koeser 06/2008 */
  HMatrix GetLinearized(const HomgPoint2D& x) const;


  /** @brief compute distance between y mapped with homography or first order
      taylor approximation at development point x0 in pixel */
  double ComputeLinearizationError(const HomgPoint2D& x0,
                                   const HomgPoint2D& y) const;

  /** @brief Set (*this) to map coordinates from camera p1 to p2 using the
   given scene plane. Infinite homography by default. 
   @author koeser 06/2008 */
  void MapAcrossPlane(const PMatrix& P1, const PMatrix& P2,
      const HomgPlane3D& scenePlane=HomgPlane3D(0,0,0,1));

  /** @brief decompose affine matrix into rotations and nonisotropic scalings
   *
   *  if result=0 then  A = theta * phi^-1 * d * phi
   *  where theta and phi are rotation matrices
   *
   * @author koeser 05/2007 */
  static int FactorizeAffineMatrixRLeft(const Matrix2x2<double>& A,
      Matrix2x2<double>& phi,
      Matrix2x2<double>& theta,
      double& d1, double& d2);

  /** @brief decompose affine matrix into rotations and nonisotropic scalings
   *
   *  if result=0 then  A = phi^-1 * d * phi * theta
   *  where theta and phi are rotation matrices
   *
   * @author koeser 05/2007 */
  static int FactorizeAffineMatrixRRight(const Matrix2x2<double>& A,
      Matrix2x2<double>& phi,
      Matrix2x2<double>& theta,
      double& d1, double& d2);

  /** @brief compose affine matrix from rotations and nonisotropic scalings
   *
   *  if result=0 then  A = theta * phi^-1 * d * phi
   *  where theta and phi must be rotation matrices
   *
   * @author koeser 10/2007 */
  static int ComposeAffineMatrixRLeft(Matrix2x2<double>& A,
      const Matrix2x2<double>& phi,
      const Matrix2x2<double>& theta,
      const double& d1, const double& d2);

  /** @brief compose affine matrix from rotations and nonisotropic scalings
   *
   *  if result=0 then  A = phi^-1 * d * phi * theta
   *  where theta and phi must be rotation matrices
   *
   * @author koeser 10/2007 */
  static int ComposeAffineMatrixRRight(Matrix2x2<double>& A,
      const Matrix2x2<double>& phi,
      const Matrix2x2<double>& theta,
      const double& d1, const double& d2);

  /** @brief compute reprojection error for both vectors of HomgPoint2D
   *
   * Computes sum of squared distances of projected points to real ones in both directions
   * 
   * @author sedlazeck 05/2008 */
  double GetReprojectionError(std::vector<BIAS::HomgPoint2D>& p1,
      std::vector<BIAS::HomgPoint2D>& p2);

};

}
// namespace

#endif // __HMatrix_hh__