vpMatrix.h

/****************************************************************************
 *
 * $Id: vpMatrix.h 4574 2014-01-09 08:48:51Z fspindle $
 *
 * This file is part of the ViSP software.
 * Copyright (C) 2005 - 2014 by INRIA. All rights reserved.
 * 
 * This software is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public License
 * ("GPL") version 2 as published by the Free Software Foundation.
 * See the file LICENSE.txt at the root directory of this source
 * distribution for additional information about the GNU GPL.
 *
 * For using ViSP with software that can not be combined with the GNU
 * GPL, please contact INRIA about acquiring a ViSP Professional 
 * Edition License.
 *
 * See http://www.irisa.fr/lagadic/visp/visp.html for more information.
 * 
 * This software was developed at:
 * INRIA Rennes - Bretagne Atlantique
 * Campus Universitaire de Beaulieu
 * 35042 Rennes Cedex
 * France
 * http://www.irisa.fr/lagadic
 *
 * If you have questions regarding the use of this file, please contact
 * INRIA at visp@inria.fr
 * 
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 *
 *
 * Description:
 * Matrix manipulation.
 *
 * Authors:
 * Eric Marchand
 *
 *****************************************************************************/



#ifndef vpMatrix_H
#define vpMatrix_H

#include <visp/vpConfig.h>
#include <visp/vpException.h>
#include <visp/vpTime.h>

#ifdef VISP_HAVE_GSL
#  include <gsl/gsl_math.h>
#  include <gsl/gsl_eigen.h>
#endif

#include <iostream>
#include <math.h>

class vpRowVector;
class vpColVector;
class vpTranslationVector;


class vpColVector;
class vpTranslationVector;
class vpRowVector;
class vpHomography;



class VISP_EXPORT vpMatrix
{
 public:
  typedef enum {
    LU_DECOMPOSITION     
  } vpDetMethod;


protected:
  unsigned int rowNum;
  unsigned int colNum;

 public:
  double *data;
 protected:
  double **rowPtrs;

  unsigned int dsize;
  unsigned int trsize;

 public:
  vpMatrix() ;
  vpMatrix(unsigned int r, unsigned int c) ;

  vpMatrix(const vpMatrix &m, unsigned int r, unsigned int c, 
       unsigned int nrows, unsigned int ncols) ;

  vpMatrix(const vpHomography& H);

  virtual ~vpMatrix();

  void  init() ;

  void kill() ;

  // Initialize an identity matrix n-by-n
  void eye(unsigned int n) ;
  // Initialize an identity matrix m-by-n
  void eye(unsigned int m, unsigned int n) ;
  void setIdentity(const double & val=1.0) ;

  //---------------------------------
  // Set/get Matrix size
  //---------------------------------
  inline unsigned int getRows() const { return rowNum ;}
  inline unsigned int getCols() const { return colNum; }

  // Set the size of the matrix A, initialization with a zero matrix
  void resize(const unsigned int nrows, const unsigned int ncols, 
          const bool nullify = true);
  
  double getMinValue() const;
  
  double getMaxValue() const;

  //---------------------------------
  // Printing
  //---------------------------------

  friend VISP_EXPORT std::ostream &operator << (std::ostream &s,const vpMatrix &m);

  int print(std::ostream& s, unsigned int length, char const* intro=0);
  std::ostream & matlabPrint(std::ostream & os) const;
  std::ostream & maplePrint(std::ostream & os) const;
  std::ostream & csvPrint(std::ostream & os) const;
  std::ostream & cppPrint(std::ostream & os, const char * matrixName = NULL, bool octet = false) const;

  void printSize() { std::cout << getRows() <<" x " << getCols() <<"  " ; }
  

  //---------------------------------
  // Copy / assignment
  //---------------------------------
  vpMatrix (const vpMatrix& m);

  // Assignment from an array
  vpMatrix &operator<<(double*);

  vpMatrix &operator=(const vpMatrix &B);
  vpMatrix &operator=(const vpHomography &H);
  vpMatrix &operator=(const double x);
  void diag(const vpColVector &A);

  //---------------------------------
  // Access/modification operators
  //---------------------------------
  inline double *operator[](unsigned int i) { return rowPtrs[i]; }
  inline double *operator[](unsigned int i) const {return rowPtrs[i];}

  //---------------------------------
  // Matrix operations (Static).
  //---------------------------------  

  static void add2Matrices(const vpMatrix &A, const vpMatrix &B, vpMatrix &C);
  static void add2WeightedMatrices(const vpMatrix &A, const double &wA, const vpMatrix &B,const double &wB, vpMatrix &C);
  static vpMatrix computeCovarianceMatrix(const vpMatrix &A, const vpColVector &x, const vpColVector &b);
  static vpMatrix computeCovarianceMatrix(const vpMatrix &A, const vpColVector &x, const vpColVector &b, const vpMatrix &w);
  static void computeHLM(const vpMatrix &H, const double &alpha, vpMatrix &HLM);

  // Create a diagonal matrix with the element of a vector DAii = Ai
  static void createDiagonalMatrix(const vpColVector &A, vpMatrix &DA)  ;
  // Insert matrix A in the current matrix at the given position (r, c).
  void insert(const vpMatrix&A, const unsigned int r, const unsigned int c);
  // Insert matrix B in matrix A at the given position (r, c).
  static vpMatrix insert(const vpMatrix &A,const  vpMatrix &B, const unsigned int r, const unsigned int c) ;
  // Insert matrix B in matrix A (not modified) at the given position (r, c), the result is given in matrix C.
  static void insert(const vpMatrix &A, const vpMatrix &B, vpMatrix &C, const unsigned int r, const unsigned int c) ;

  // Juxtapose to matrices C = [ A B ]
  static vpMatrix juxtaposeMatrices(const vpMatrix &A,const  vpMatrix &B) ;
  // Juxtapose to matrices C = [ A B ]
  static void juxtaposeMatrices(const vpMatrix &A,const  vpMatrix &B, vpMatrix &C) ;
  // Compute Kronecker product matrix
  static void kron(const vpMatrix  &m1, const vpMatrix  &m2 , vpMatrix  &out);

  // Compute Kronecker product matrix
  static vpMatrix kron(const vpMatrix  &m1, const vpMatrix  &m2 );

  static inline bool loadMatrix(std::string filename, vpMatrix &M,
                                const bool binary = false, char *Header = NULL)
  {
    return vpMatrix::loadMatrix(filename.c_str(), M, binary, Header);
  }
  static bool loadMatrix(const char *filename, vpMatrix &M, const bool binary = false, char *Header = NULL);
  static void mult2Matrices(const vpMatrix &A, const vpMatrix &B, vpMatrix &C);
  static void multMatrixVector(const vpMatrix &A, const vpColVector &b, vpColVector &c);
  static void negateMatrix(const vpMatrix &A, vpMatrix &C);
  static inline bool saveMatrix(std::string filename, const vpMatrix &M,
        const bool binary = false,
        const char *Header = "")
  {
    return vpMatrix::saveMatrix(filename.c_str(), M, binary, Header);
  }
  static bool saveMatrix(const char *filename, const vpMatrix &M, const bool binary = false, const char *Header = "");

  static inline bool saveMatrixYAML(std::string filename, const vpMatrix &M, const char *header = "")
  {
    return vpMatrix::saveMatrixYAML(filename.c_str(), M, header);
  }
  static bool saveMatrixYAML(const char *filename, const vpMatrix &M, const char *header = "");
  static inline bool loadMatrixYAML(std::string filename, vpMatrix &M, char *header = NULL)
  {
    return vpMatrix::loadMatrixYAML(filename.c_str(), M, header);
  }
  static bool loadMatrixYAML(const char *filename, vpMatrix &M, char *header = NULL);


  // Stack the matrix A below the current one, copy if not initialized this = [ this A ]^T
  void stackMatrices(const vpMatrix &A);
  static vpMatrix stackMatrices(const vpMatrix &A,const  vpMatrix &B) ;
  static void stackMatrices(const vpMatrix &A,const  vpMatrix &B, vpMatrix &C) ;
  static void sub2Matrices(const vpMatrix &A, const vpMatrix &B, vpMatrix &C);
  

  //---------------------------------
  // Matrix operations.
  //---------------------------------
  // operation A = A + B
  vpMatrix &operator+=(const vpMatrix &B);
  // operation A = A - B
  vpMatrix &operator-=(const vpMatrix &B);

  vpMatrix operator*(const vpMatrix &B) const;
  vpMatrix operator*(const vpHomography &H) const;
  vpMatrix operator+(const vpMatrix &B) const;
  vpMatrix operator-(const vpMatrix &B) const;
  vpMatrix operator-() const;

  //---------------------------------
  // Matrix/vector operations.
  //---------------------------------

  vpColVector operator*(const vpColVector &b) const;
  // operation c = A * b (A is unchanged, c and b are translation vectors)
  vpTranslationVector operator*(const vpTranslationVector  &b) const;
  //---------------------------------
  // Matrix/real operations.
  //---------------------------------

  vpMatrix &operator+=(const double x);
  vpMatrix &operator-=(const double x);
  vpMatrix &operator*=(const double x);
  vpMatrix &operator/=(double x);

  // Cij = Aij * x (A is unchanged)
  vpMatrix operator*(const double x) const;
  // Cij = Aij / x (A is unchanged)
  vpMatrix operator/(const double x) const;

  double sumSquare() const;

  // return the determinant of the matrix.
  double det(vpDetMethod method = LU_DECOMPOSITION) const;
  
  //Compute the exponential matrix of a square matrix
  vpMatrix expm();

  //-------------------------------------------------
  // Columns, Rows extraction, SubMatrix
  //-------------------------------------------------
  vpRowVector row(const unsigned int i);
  vpColVector column(const unsigned int j);
  void init(const vpMatrix &m, unsigned int r, unsigned int c, 
        unsigned int nrows, unsigned int ncols);

  //-------------------------------------------------
  // transpose
  //-------------------------------------------------
  // Compute the transpose C = A^T
  vpMatrix t() const;

  // Compute the transpose C = A^T
  vpMatrix transpose()const;
  void  transpose(vpMatrix & C )const;
    
  vpMatrix AAt() const;
  void AAt(vpMatrix &B) const;
   
  vpMatrix AtA() const;
  void AtA(vpMatrix &B) const;


  //-------------------------------------------------
  // Kronecker product
  //-------------------------------------------------
  
  // Stacks columns of a matrix in a vector
  void stackColumns(vpColVector  &out );

  // Stacks columns of a matrix in a vector
  vpColVector stackColumns();

  // Stacks columns of a matrix in a vector
  void stackRows(vpRowVector  &out );

  // Stacks columns of a matrix in a vector
  vpRowVector stackRows();
  
  // Compute Kronecker product matrix 
  void kron(const vpMatrix  &m1 , vpMatrix  &out);
  
  // Compute Kronecker product matrix 
  vpMatrix kron(const vpMatrix  &m1);
  
  //-------------------------------------------------
  // Matrix inversion
  //-------------------------------------------------
#ifndef DOXYGEN_SHOULD_SKIP_THIS
  void LUDcmp(unsigned int* perm, int& d);
  void LUBksb(unsigned int* perm, vpColVector& b);

#endif // doxygen should skip this
  // inverse matrix A using the LU decomposition 
  vpMatrix inverseByLU() const;
#if defined(VISP_HAVE_LAPACK)
  // inverse matrix A using the Cholesky decomposition (only for real symmetric matrices)
  vpMatrix inverseByCholesky() const;
  //lapack implementation of inverse by Cholesky
  vpMatrix inverseByCholeskyLapack() const;
  // inverse matrix A using the QR decomposition
  vpMatrix inverseByQR() const;
  //lapack implementation of inverse by QR
  vpMatrix inverseByQRLapack() const;
#endif
  vpMatrix pseudoInverse(double svThreshold=1e-6)  const;
  unsigned int pseudoInverse(vpMatrix &Ap, double svThreshold=1e-6)  const;
  unsigned int pseudoInverse(vpMatrix &Ap, vpColVector &sv, double svThreshold=1e-6) const ;
  unsigned int pseudoInverse(vpMatrix &Ap,
                             vpColVector &sv, double svThreshold,
                             vpMatrix &ImA,
                             vpMatrix &ImAt) const ;
  unsigned int pseudoInverse(vpMatrix &Ap,
                             vpColVector &sv, double svThreshold,
                             vpMatrix &ImA,
                             vpMatrix &ImAt,
                             vpMatrix &kerA) const ;

  //-------------------------------------------------
  // SVD decomposition
  //-------------------------------------------------

#ifndef DOXYGEN_SHOULD_SKIP_THIS
  void svdFlake(vpColVector& w, vpMatrix& v);
  void svdNr(vpColVector& w, vpMatrix& v);
#ifdef VISP_HAVE_GSL
  void svdGsl(vpColVector& w, vpMatrix& v);
#endif
#if (VISP_HAVE_OPENCV_VERSION >= 0x020101) // Require opencv >= 2.1.1
  void svdOpenCV(vpColVector& w, vpMatrix& v);
#endif
#ifdef VISP_HAVE_LAPACK
  void svdLapack(vpColVector& W, vpMatrix& V);
#endif
  void SVBksb(const vpColVector& w,
          const vpMatrix& v,
          const vpColVector& b, vpColVector& x);
#endif

  // singular value decomposition SVD

  void svd(vpColVector& w, vpMatrix& v);

  // solve Ax=B using the SVD decomposition (usage A = solveBySVD(B,x) )
  void solveBySVD(const vpColVector &B, vpColVector &x) const ;
  // solve Ax=B using the SVD decomposition (usage  x=A.solveBySVD(B))
  vpColVector solveBySVD(const vpColVector &B) const ;

  unsigned int kernel(vpMatrix &KerA, double svThreshold=1e-6);
  double cond();

  //-------------------------------------------------
  // Eigen values and vectors
  //-------------------------------------------------

  // compute the eigen values using the Gnu Scientific library
  vpColVector eigenValues();
  void eigenValues(vpColVector &evalue, vpMatrix &evector);

  // -------------------------
  // Norms
  // -------------------------
  // Euclidean norm ||x||=sqrt(sum(x_i^2))
  double euclideanNorm () const;
  // Infinity norm ||x||=max(sum(fabs(x_i)))
  double infinityNorm () const;

 private:
  double detByLU() const;

};






VISP_EXPORT vpMatrix operator*(const double &x, const vpMatrix &A) ;

VISP_EXPORT vpColVector operator*(const double &x, const vpColVector &A) ;



#endif


/*
 * Local variables:
 * c-basic-offset: 2
 * End:
 */