---

matrix.h

Go to the documentation of this file.

#ifndef _MATRIX_H_
#define _MATRIX_H_

/*
  Choose at least one matrix package; currently MTL and TNT are
  supported.

  Enabling more than one package is legal and can be useful if you
  want to support more than one matrix implementation, e.g. for
  different source files.  Note that the type generator (by design)
  and the MSEQ definitions (due to MTL's broken iterator concept) only
  use a single implementation regardless.
*/
#define USE_TNT_MATRIX
//#define USE_MTL_MATRIX



#include <numeric>


#ifdef USE_TNT_MATRIX

#define num_cols ncols
#define num_rows nrows
#include "tnt/tnt.h"
#include "tnt/cmat.h"
#include "tnt/vec.h"
#include "tnt/region1d.h"
#include "tnt/region2d.h"
#undef  num_cols
#undef  num_rows
#endif


#ifdef USE_MTL_MATRIX
#include "mtl/mtl.h"
#include "mtl/scaled2D.h"
#include "mtl/matrix.h"
#include "mtl/flatten_iterator.h"
#include "mtl/flatten_iterator.h"
#define MSEQ(container)      mtl::flatten_begin(container),mtl::flatten_end(container)

#define CMSEQ(container)     mtl::const_flatten_begin(container),mtl::const_flatten_end(container)

#define MTL_MSEQ(container)  mtl::flatten_begin(container),mtl::flatten_end(container)

#define MTL_CMSEQ(container) mtl:::const_flatten_begin(container),mtl:::const_flatten_end(container)

#define MBEGIN(container)    mtl::flatten_begin(container)

#define CMBEGIN(container)   mtl::const_flatten_begin(container)
#endif


#ifdef USE_TNT_MATRIX

#define SubMatrixType const_Region

#define MatIndex(i) (i+1)

#define UpperBound(i) (i)

#else

#define SubMatrixType submatrix_type

#define MatIndex(i) (i)

#define UpperBound(i) (i)
#endif



#ifndef  MSEQ

#define  MSEQ(container) (container).begin(),(container).end()

#define CMSEQ(container) (container).begin(),(container).end()

#define  MBEGIN(container)  (container).begin()

#define  CMBEGIN(container) (container).begin()
#endif


/******************************************************************************/
/* Type generator                                                             */
/******************************************************************************/
template<class T=double>
struct MatrixType {


  
#ifdef USE_MTL_MATRIX
  /* Use MTL if available */
  typedef typename mtl::matrix< T, mtl::rectangle<>, mtl::dense<>      >::type rectangular;
  typedef typename mtl::matrix< T, mtl::rectangle<>, mtl::dense<>      >::type circular;
  typedef typename mtl::matrix< T, mtl::rectangle<>, mtl::compressed<> >::type sparse;

#else
  /* Default to TNT */
  typedef TNT::Matrix<T> rectangular;
  typedef TNT::Matrix<T> circular;  
  typedef TNT::Matrix<T> sparse;  
#endif

};



/******************************************************************************/
/* Additional routines for namespace mtl                                      */
/******************************************************************************/

namespace mtl {

template <class Matrix> inline
Matrix& operator*=(Matrix &A, const typename Matrix::value_type B)
{  mtl::scale(A,B);  return A;  }

 

template <class Matrix> inline
Matrix& operator+=(Matrix &A, const typename Matrix::value_type B)
{
  typedef typename Matrix::size_type Subscript;

  for (Subscript i=0; i<A.nrows(); i++)
    for (Subscript j=0; j<A.ncols(); j++)
      A[i][j] += B;

  return A;
}


template <class Matrix> inline
Matrix& operator-=(Matrix &A, const typename Matrix::value_type B)
{  return A += -1*B;  }


} /* namespace mtl */



/******************************************************************************/
/* Additional routines for namespace TNT                                      */
/******************************************************************************/

#ifdef USE_TNT_MATRIX

namespace TNT {
  
template <class T>
inline Matrix<T>& operator+=(Matrix<T>  &A, const T B)
{
  for (Subscript i=0; i<A.nrows(); i++)
    for (Subscript j=0; j<A.ncols(); j++)
      A[i][j] += B;

  return A;
}

template <class T>
inline Matrix<T>& operator-=(Matrix<T>  &A, const T B)
{
  for (Subscript i=0; i<A.nrows(); i++)
    for (Subscript j=0; j<A.ncols(); j++)
      A[i][j] -= B;

  return A;
}

template <class T>
inline Matrix<T>& operator*=(Matrix<T>  &A, const T B)
{
  for (Subscript i=0; i<A.nrows(); i++)
    for (Subscript j=0; j<A.ncols(); j++)
      A[i][j] *= B;

  return A;
}

template <class T>
inline Matrix<T> operator*(const Matrix<T>  &A, const T B)
{
  Subscript M = A.nrows();
  Subscript N = A.ncols();
  
  Matrix<T> tmp(M,N);
  
  for (Subscript i=0; i<M; i++)
    for (Subscript j=0; j<N; j++)
      tmp[i][j] = A[i][j]*B; 
  
  return tmp;
}

template <class T> inline T sum(const Matrix<T>& A)
{  return std::accumulate(A.begin(),A.end(),T(0));  }
 
template <class T> inline T sum(const Vector<T>& A)
{  return std::accumulate(A.begin(),A.end(),T(0));  }
 
template <class T> inline T sum(const_Region2D< Matrix<T> >& A)
{
  Subscript start = A.lbound();
  Subscript Mend  = A.lbound() + A.nrows() - 1;
  Subscript Nend  = A.lbound() + A.ncols() - 1;
  T partialsum = T(0);
  
  for (Subscript i=start; i<=Mend; i++)
    for (Subscript j=start; j<=Nend; j++)
        partialsum += A(i,j);
  
  return partialsum;
}

} /* namespace TNT */
#endif




/******************************************************************************/
/* Routines for namespace mat                                                 */
/******************************************************************************/
namespace mat {

  
template <class Matrix, class InnerC, class OuterC> inline
Matrix matrix_from_nested_containers(const OuterC& n, bool column_major=false)
{
  typedef typename Matrix::size_type Subscript;

  /* Find maximum size of a nested container */
  Subscript max_minor = 0;
  for (typename OuterC::const_iterator i=n.begin(); i!=n.end(); i++) {
    const Subscript new_size =  (*i ? ((*i)->size()) : 0);
    if (new_size > max_minor) max_minor=new_size;
  }

  /* Create and populate new Matrix */
  Matrix m((column_major? max_minor  : n.size()  ),
           (column_major? n.size()   : max_minor));
  int major=0;
  for (    typename OuterC::const_iterator i=n.begin();         i!=n.end(); i++,major++) {
    int minor=0;
    if (*i)
      for (typename InnerC::const_iterator j=(*i)->begin(); j!=(*i)->end(); j++,minor++) {
        if (column_major) 
          m[minor][major] = *j;
        else
          m[major][minor] = *j;
      }
  }
    
  return m;
}



template <class Matrix> inline
typename Matrix::value_type
sum(const Matrix& A)
{  return std::accumulate(CMSEQ(A), typename Matrix::value_type(0) );  }



template <class Matrix> inline
typename Matrix::size_type
size(const Matrix& A)
{  return (A.nrows()*A.ncols());  }



template <class Matrix> inline
Matrix& set(Matrix &A, const typename Matrix::value_type B)
{
  typedef typename Matrix::size_type Subscript;

  for (Subscript i=0; i<A.nrows(); i++)
    for (Subscript j=0; j<A.ncols(); j++)
      A[i][j] = B;

  return A;
}



template <class Matrix> inline
typename Matrix::value_type elem(const Matrix &A, typename Matrix::size_type i, typename Matrix::size_type j)
{
  return A(MatIndex(i),MatIndex(j));
}



template <class Matrix> inline
const typename Matrix::SubMatrixType
submatrix(const Matrix &A,
          typename Matrix::size_type i1, typename Matrix::size_type i2,
          typename Matrix::size_type j1, typename Matrix::size_type j2)
{
#ifdef USE_TNT_MATRIX
  TNT::Index1D I(MatIndex(i1),UpperBound(i2));
  TNT::Index1D J(MatIndex(j1),UpperBound(j2));  
  return A(I,J);
#else
  return A.sub_matrix(MatIndex(i1),UpperBound(i2),MatIndex(j1),UpperBound(j2));
#endif
}



#ifdef USE_TNT_MATRIX
template <class T> inline
TNT::Subscript
size(TNT::const_Region2D< TNT::Matrix<T> >& A)
{  return (A.lbound() + A.nrows() - 1)*(A.lbound() + A.ncols() - 1);  }
#endif
 


template <class Matrix>
typename Matrix::value_type
edge_average(const Matrix& mat)
{
  typedef typename Matrix::value_type T;
  typedef typename Matrix::size_type  Subscript;
  
  Subscript X=mat.nrows();
  Subscript Y=mat.ncols();

  /* Return zero in case of empty matrix to avoid divide by zero */
  if (X<1 || Y<1) return T(0);
  
  typename Matrix::submatrix_type left   = mat.sub_matrix(0,   1,   0,   Y-1);
  typename Matrix::submatrix_type top    = mat.sub_matrix(1,   X,   0,   1  );
  typename Matrix::submatrix_type right  = mat.sub_matrix(X-1, X,   1,   Y  );
  typename Matrix::submatrix_type bottom = mat.sub_matrix(0,   X-1, Y-1, Y  );

  const T         total =  mat::sum(left)+  mat::sum(top)+  mat::sum(right)+  mat::sum(bottom);
  const Subscript count = mat::size(left)+ mat::size(top)+ mat::size(right)+ mat::size(bottom);

  return total/count;
}


#ifdef USE_TNT_MATRIX
template <class T>
T edge_average(const TNT::Matrix<T>& mat)
{
  using namespace TNT;
  
  Subscript X=mat.nrows();
  Subscript Y=mat.ncols();

  /* Image must be at least 2x2 to have an edge */
  if (X<2 || Y<2) return T(0);
  
  const_Region2D< Matrix<double> > left(   mat, 1, 1,   1, Y-1);
  const_Region2D< Matrix<double> > top(    mat, 2, X,   1, 1  );
  const_Region2D< Matrix<double> > right(  mat, X, X,   2, Y  );
  const_Region2D< Matrix<double> > bottom( mat, 1, X-1, Y, Y  );

  const T         total =  TNT::sum(left)+  TNT::sum(top)+  TNT::sum(right)+  TNT::sum(bottom);  
  const Subscript count = mat::size(left)+ mat::size(top)+ mat::size(right)+ mat::size(bottom);

  return total/count;
}
#endif


/* Returns the maximum num_rows of any element in the given matrix or region.
   
  There may be a way to code this using a generic mem_fun() call, but
  it wasn't obvious how to specify that call since we don't know the
  actual Matrix::value_type, i.e. the type of an element.

  On second thought, if there are iterators available, it should work;
  we would just then accept a pair of iterators instead of the
  container itself.  However, I'm not sure if one-dimensional
  iterators are defined for regions in TNT and MTL.
*/
template <class Matrix> inline
typename Matrix::size_type max_nrows(const Matrix &A)
{
  typedef typename Matrix::size_type Subscript;
  Subscript val=0;
  
  for (Subscript i=0; i<A.nrows(); i++)
    for (Subscript j=0; j<A.ncols(); j++)
      if (mat::elem(A,i,j).nrows()>val)
        val = mat::elem(A,i,j).nrows();

  return val;
}



/* Returns the maximum num_cols of any element in the given matrix or region.
   
  As for max_nrows, this routine *should* be unnecessary.
*/
template <class Matrix> inline
typename Matrix::size_type max_ncols(const Matrix &A)
{
  typedef typename Matrix::size_type Subscript;
  Subscript val=0;
  
  for (Subscript i=0; i<A.nrows(); i++)
    for (Subscript j=0; j<A.ncols(); j++)
      if (mat::elem(A,i,j).ncols()>val)
        val = mat::elem(A,i,j).ncols();

  return val;
}


/* Returns the maximum num_rows of any pointed-to element in the given
  matrix or region.
   
  There may be a way to code this using a generic mem_fun() call, but
  it wasn't obvious how to specify that call since we don't know the
  actual Matrix::value_type, i.e. the type of an element.
*/
template <class Matrix> inline
typename Matrix::size_type max_pnrows(const Matrix &A)
{
  typedef typename Matrix::size_type Subscript;
  Subscript val=0;
  
  for (Subscript i=0; i<A.nrows(); i++)
    for (Subscript j=0; j<A.ncols(); j++) {
      typename Matrix::value_type p=mat::elem(A,i,j);
      if (p && (p->nrows()>val))
        val = p->nrows();
    }
  
  return val;
}



/* Returns the maximum num_cols of any pointed-to element
  in the given matrix or region.
   
  As for max_pnrows, this routine *should* be unnecessary.
*/
template <class Matrix> inline
typename Matrix::size_type max_pncols(const Matrix &A)
{
  typedef typename Matrix::size_type Subscript;
  Subscript val=0;
  
  for (Subscript i=0; i<A.nrows(); i++)
    for (Subscript j=0; j<A.ncols(); j++) {
      typename Matrix::value_type p=mat::elem(A,i,j);
      if (p && (p->ncols()>val))
        val = p->ncols();
    }

  return val;
}

template<class Matrix>
void delete_contents(Matrix& A)
{
  typedef typename Matrix::size_type Subscript;
  for (Subscript i=0; i<A.nrows(); i++)
    for (Subscript j=0; j<A.ncols(); j++)
      delete A[i][j];
}


template<class Matrix, class OutputStream>
OutputStream& stream_output(const Matrix& A, OutputStream& stream)
{
  typedef typename Matrix::size_type Subscript;
  for (Subscript i=0; i<A.nrows(); i++) {
    for (Subscript j=0; j<A.ncols(); j++)
      stream << mat::elem(A,i,j) << " ";
    stream << std::endl;
  }

  return stream;
}

template<class Matrix, class InputStream>
InputStream& stream_input(Matrix& A, InputStream& stream)
{
  typedef typename Matrix::size_type Subscript;
  for (Subscript i=0; i<A.nrows(); i++) {
    for (Subscript j=0; j<A.ncols(); j++)
      stream >> A[i][j];
  }

  return stream;
}


/* Avoid namespace pollution */
//#undef MatIndex
#undef UpperBound

} /* namespace mat */
#endif

---