MeVisLab/Standard/Sources/ML/MLVesselGraph/mlMatrixTemplate.h

Go to the documentation of this file.

/*
 * Copyright 2002, MeVis gGmbH
 * ALL RIGHTS RESERVED
 *
 * THE CONTENT OF THIS WORK CONTAINS CONFIDENTIAL AND PROPRIETARY
 * INFORMATION OF MEVIS GGMBH. ANY DUPLICATION, MODIFICATION,
 * DISTRIBUTION, OR DISCLOSURE IN ANY FORM, IN WHOLE, OR IN PART, IS STRICTLY
 * PROHIBITED WITHOUT THE PRIOR EXPRESS WRITTEN PERMISSION OF MEVIS GGMBH
 *
 */

//======================================================

//======================================================
#ifndef __MatrixTemplate_h__
#define __MatrixTemplate_h__

#include "mlInitSystemML.h"
#ifdef UNIX
#include "mlDebug.h"
#endif

#include "mlSystemWarningsDisable.h"
#include <iostream>
#include <valarray>
#include <numeric>
#include <functional>
#include <algorithm>
#include <iterator>
#include "mlSystemWarningsRestore.h"

//======================================================================
// namespaces
//======================================================================
ML_START_NAMESPACE
//======================================================================
// declaration of class *Slice_iter*
//======================================================================
template<class T>
class Slice_iter {
  std::valarray<T>* v;                            
  std::slice s;                                   
  size_t curr;                                    

  T& ref(size_t i) const { return (*v)[s.start()+i*s.stride()]; }
public:
  Slice_iter(std::valarray<T>* vv, std::slice ss) :v(vv), s(ss), curr(0) { }

  Slice_iter<T> end() const
  {
    Slice_iter<T> t = *this;
    t.curr = s.size();                            // index of one plus last element
    return t;
  }


  Slice_iter<T>& operator++() { curr++; return *this; }                         
  Slice_iter<T> operator++(int) { Slice_iter<T> t = *this; curr++; return t; }  



  T& operator[](size_t i) { return ref(i); }      
  T& operator()(size_t i) { return ref(i); }      
  T& operator*() { return ref(curr); }            



  // caution: friend function has to be implemented here inplace to support
  // explicit template specification according to ansi standart in gcc, i.e to
  // let the compiler know, that the friend functions are realy templated overloads
  // (note: the gcc conform syntax (friend bool operator== <> (...)) is not supported
  // on other compiler)
  friend bool operator== (const Slice_iter<T>& p, const Slice_iter<T>& q)   
  {
    return p.curr==q.curr
      && p.s.stride()==q.s.stride()
      && p.s.start()==q.s.start();
  }

  friend bool operator!= (const Slice_iter<T>& p, const Slice_iter<T>& q)   
  {
    return !(p==q);
  }

  friend bool operator< (const Slice_iter<T>& p, const Slice_iter<T>& q)   
  {
    return p.curr<q.curr
      && p.s.stride()==q.s.stride()
      && p.s.start()==q.s.start();
  }
};


//======================================================================
// declaration of class *Cslice_iter*
//======================================================================
template<class T>
class Cslice_iter
{
  std::valarray<T>* v;                            
  std::slice s;                                   
  size_t curr;                                    

  const T& ref(size_t i) const { return (*v)[s.start()+i*s.stride()]; }
public:
  Cslice_iter(std::valarray<T>* vv, std::slice ss): v(vv), s(ss), curr(0){}

  Cslice_iter<T> end() const
  {
    Cslice_iter<T> t = *this;
    t.curr = s.size();                            // index of one plus last element
    return t;
  }


  Cslice_iter<T>& operator++() { curr++; return *this; }                          
  Cslice_iter<T> operator++(int) { Cslice_iter<T> t = *this; curr++; return t; }  



  const T& operator[](size_t i) const { return ref(i); }  
  const T& operator()(size_t i) const { return ref(i); }  
  const T& operator*() const { return ref(curr); }        



  // caution: friend function has to be implemented here inplace to support
  // explicit template specification according to ansi standart in gcc, i.e to
  // let the compiler know, that the friend functions are realy templated overloads
  // (note: the gcc conform syntax (friend bool operator== <> (...)) is not supported
  // on other compiler)
  friend bool operator==(const Cslice_iter<T>& p, const Cslice_iter<T>& q) 
  {
    return p.curr==q.curr
      && p.s.stride()==q.s.stride()
      && p.s.start()==q.s.start();
  }

  friend bool operator!=(const Cslice_iter<T>& p, const Cslice_iter<T>& q) 
  {
    return !(p==q);
  }

  friend bool operator< (const Cslice_iter<T>& p, const Cslice_iter<T>& q) 
  {
    return p.curr<q.curr
      && p.s.stride()==q.s.stride()
      && p.s.start()==q.s.start();
  }

};

//======================================================================
// declaration of class *MatrixTemplate*
//======================================================================
template <class T>
class MatrixTemplate {
  std::valarray<T>* v;                            
  size_t size_x;                                  
  size_t size_y;                                  
public:
  MatrixTemplate(size_t x, size_t y);             
  MatrixTemplate(const MatrixTemplate<T>&);       
  MatrixTemplate<T>& operator=(const MatrixTemplate<T>&);   
  ~MatrixTemplate();                              

  void freeMatrix();                                   
  void resizeMatrix(size_t x, size_t y, const T& c = T());     
  void overrideMatrix(size_t x, size_t y, const T& c = T());   

  size_t size() const { return size_x*size_y; }   
  size_t sizeX() const { return size_x; }         
  size_t sizeY() const { return size_y; }         

  Slice_iter<T> row(size_t i);                    
  Cslice_iter<T> row(size_t i) const;             

  Slice_iter<T> column(size_t i);                 
  Cslice_iter<T> column(size_t i) const;          


  T& operator()(size_t x, size_t y);              
  T operator()(size_t x, size_t y) const;         

  Slice_iter<T> operator()(size_t i) { return column(i); }          
  Cslice_iter<T> operator()(size_t i) const { return column(i); }   



  Slice_iter<T> operator[](size_t i) { return column(i); }          
  Cslice_iter<T> operator[](size_t i) const { return column(i); }   


  MatrixTemplate<T>& operator*=(T);               

  std::valarray<T>& array() { return *v; }        

  operator T * () { return (T*) (&v); }           
};


//======================================================================
// declaration of class *MatrixSizedTemplate*
//======================================================================
template <class T, size_t SIZE_X, size_t SIZE_Y>
class MatrixSizedTemplate : public MatrixTemplate<T>
{
  public:
    MatrixSizedTemplate () : MatrixTemplate<T>(SIZE_X, SIZE_Y) {;}
};






//=======================================================================
// Implementation of MatrixTemplate - functions
//=======================================================================

template <class T>
inline Slice_iter<T> MatrixTemplate<T>::row(size_t i)
{
  return Slice_iter<T>(v,std::slice(i,size_x,size_y));
}

template <class T>
inline Cslice_iter<T> MatrixTemplate<T>::row(size_t i) const
{
  return Cslice_iter<T>(v,std::slice(i,size_x,size_y));
}

template <class T>
inline Slice_iter<T> MatrixTemplate<T>::column(size_t i)
{
  return Slice_iter<T>(v,std::slice(i*size_y,size_y,1));
}

template <class T>
inline Cslice_iter<T> MatrixTemplate<T>::column(size_t i) const
{
  return Cslice_iter<T>(v,std::slice(i*size_y,size_y,1));
}

template <class T>
MatrixTemplate<T>::MatrixTemplate(size_t x, size_t y)
{
  size_x = x;
  size_y = y;
  v = new std::valarray<T>(size_x*size_y);
}

template <class T>
MatrixTemplate<T>::MatrixTemplate(const MatrixTemplate<T>&mt)
{
  // set size according to transfered object
  size_x = mt.size_x;
  size_y = mt.size_y;
  // assign new valarray of appropriate size
  v = new std::valarray<T>(size_x*size_y);

  // copy memory content
  *v = *(mt.v);
}

template <class T>
MatrixTemplate<T>& MatrixTemplate<T>::operator=(const MatrixTemplate<T>&mt)
{
  if(this != &mt){
    // resize matrixTemplate
    this->resizeMatrix(mt.size_x, mt.size_y);
    // copy memory content
    *v = *(mt.v);
  }
  return *this;
}

template <class T>
MatrixTemplate<T>::~MatrixTemplate()
{
  delete v;
}

template <class T>
void MatrixTemplate<T>::resizeMatrix(size_t x, size_t y, const T& c)
{
  size_x = x;
  size_y = y;
  v->resize(size_x*size_y, c);
}

template <class T>
void MatrixTemplate<T>::overrideMatrix(size_t x, size_t y, const T& c)
{
  size_x = x;
  size_y = y;
  *v = c;
  v->resize(size_x*size_y, c);
}


template <class T>
void MatrixTemplate<T>::freeMatrix()
{
  size_x = 0;
  size_y = 0;
  v->free();
}

template <class T>
T& MatrixTemplate<T>::operator()(size_t x, size_t y)
{
  return column(x)[y];
}



//---------------------------------------------------------------------------------
//---------------------------------------------------------------------------------
template <class T>
T mul(const Cslice_iter<T>& v1, const std::valarray<T>& v2)
{
  T res = 0;
  for (size_t i = 0; i<v2.size(); i++) res+= v1[i]*v2[i];
  return res;
}


template <class T>
std::valarray<T> operator*(const MatrixTemplate<T>& m, const std::valarray<T>& v)
{
  if (m.sizeX()!=v.size()){
    ML_PRINT_ERROR("mlMatrixTemplate.h, operator*() # 1",
                   ML_BAD_PARAMETER,
                   "Wrong number of elements in m*v\n, ignoring problem");
  }

  std::valarray<T> res(m.sizeY());
  for (size_t i = 0; i<m.sizeY(); i++) res[i] = mul(m.row(i),v);
  return res;
}


//std::valarray<T> operator*(const MatrixTemplate<T>& m, std::valarray<T>& v)
template <class T>
std::valarray<T> mul_mv(const MatrixTemplate<T>& m, std::valarray<T>& v)
{
  if (m.sizeX()!=v.size()){
    ML_PRINT_ERROR("mlMatrixTemplate.h, mul_mv()",
                   ML_BAD_PARAMETER,
                   "Wrong number of elements in m*v\n, ignoring problem");
  }

  std::valarray<T> res(m.sizeY());

  for (size_t i = 0; i<m.sizeY(); i++) {
    const Cslice_iter<T>& ri = m.row(i);
    res[i] = inner_product(ri,ri.end(),&v[0],T(0));
  }
  return res;
}



template <class T>
std::valarray<T> operator*(std::valarray<T>& v, const MatrixTemplate<T>& m)
{
  if (v.size()!=m.sizeY()){
    ML_PRINT_ERROR("mlMatrixTemplate.h, operator*() # 2",
                   ML_BAD_PARAMETER,
                   "Wrong number of elements in v*m\n, ignoring problem");
  }

  std::valarray<T> res(m.sizeX());

  for (size_t i = 0; i<m.sizeX(); i++) {
    const Cslice_iter<T>& ci = m.column(i);
    res[i] = inner_product(ci,ci.end(),&v[0],T(0));
  }
  return res;
}

template <class T>
MatrixTemplate<T>& MatrixTemplate<T>::operator*=(T d)
{
  (*v) *= d;
  return *this;
}

template <class T>
std::ostream& operator<<(std::ostream& os, MatrixTemplate<T>& m)
{
  for(int y=0; y<m.sizeY(); y++)
  {
    for(int x=0; x<m.sizeX(); x++)
      os<<m[x][y]<<"\t";
    os << "\n";
  }
  return os;
}





} // namespace ml

#endif //__MatrixTemplate_h__