Tensor3D3x3x3.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 __Tensor3D3x3x3_hh__
#define __Tensor3D3x3x3_hh__

#include "Tensor3D.hh"
#include "Matrix3x3.hh"
#include "Vector3.hh"

namespace BIAS {
    /**
     * \class Tensor3D3x3x3
     * \brief
     * \author Felix Woelk
     */
    template <class T=double>
    class BIASMathBase_EXPORT  Tensor3D3x3x3 : public Tensor3D<T>
  {
  public:
    Tensor3D3x3x3<T>();    
    virtual ~Tensor3D3x3x3<T>();

    void NewSize(const int i, const int j, const int k,
                 const T& def=(T)0.0);
    void newsize(const int i, const int j, const int k, 
                 const T& def=(T)0.0);

    /** computes the contraction of index 1 from left with contravariant
        vector , i.e. res[k][j] = sum_i v[i] * (*this)[k][j]
        @author woelk 07/2004 */
    inline void LeftContravariantContraction1(Vector3<T>& v, 
                                              Matrix3x3<T>& res);
    inline void LeftContravariantContraction2(Vector3<T>& v, 
                                              Matrix3x3<T>& res);
    inline void LeftContravariantContraction3(Vector3<T>& v, 
                                              Matrix3x3<T>& res);
    inline Matrix3x3<T> LeftContravariantContraction1(Vector3<T>& v)
    {
      Matrix3x3<T> M;
      LeftContravariantContraction1(v,M);
      return M;
    }
    
    inline Matrix3x3<T> LeftContravariantContraction2(Vector3<T>& v)
    {
      Matrix3x3<T> M;
      LeftContravariantContraction2(v,M);
      return M;
    }

    inline Matrix3x3<T> LeftContravariantContraction3(Vector3<T>& v)
    {
      Matrix3x3<T> M;
      LeftContravariantContraction3(v,M);
      return M;
    }

    /** @brief Compute S_{ijk} = T_{pjk} M_{pi}
        @author mdunda */
    void postmultiply1(const BIAS::Matrix3x3<T> &M,
                       Tensor3D3x3x3<T> &result) const;

    /** @brief Compute S_{ijk} = T_{ipk} M_{pj}
        @author mdunda */
    void postmultiply2(const BIAS::Matrix3x3<T> &M,
                       Tensor3D3x3x3<T> &result) const;

    /** @brief Compute S_{ijk} = T_{ijp} M_{pk}
        @author mdunda */
    void postmultiply3(const BIAS::Matrix3x3<T> &M,
                       Tensor3D3x3x3<T> &result) const;

    /** @brief compute S_{ijk} = M_{ip} T_{pjk}
        @author mdunda */
    void premultiply1(const BIAS::Matrix3x3<T> &M,
                      Tensor3D3x3x3<T> &result) const;

    /** @brief compute S_{ijk} = M_{jp} T_{ipk}
        @author mdunda */
    void premultiply2(const BIAS::Matrix3x3<T> &M,
                      Tensor3D3x3x3<T> &result) const;

    /** @brief compute S_{ijk} = M_{kp} T_{ijp}
        @author mdunda */
    void premultiply3(const BIAS::Matrix3x3<T> &M,
                      Tensor3D3x3x3<T> &result) const;

    inline void SetFromMatrices(const Matrix3x3<T>& T1, const Matrix3x3<T>& T2,
                                const Matrix3x3<T>& T3);

    inline void GetMatrices(Matrix3x3<T>& T1, Matrix3x3<T>& T2, 
                            Matrix3x3<T>& T3) const;

    inline Tensor3D3x3x3<T>& operator=(const Tensor3D3x3x3<T>& t)
    { Tensor3D<T>::operator=(t); return *this; }

  protected:
  
  }; // class

  //////////////////////////////////////////////////////////////////////////
  ///            implementation
  //////////////////////////////////////////////////////////////////////////

  template <class T> 
  void Tensor3D3x3x3<T>::SetFromMatrices(const Matrix3x3<T>& T1,
                                         const Matrix3x3<T>& T2,
                                         const Matrix3x3<T>& T3)
  {
    size_t size=9*sizeof(T);
    memcpy((void *)this->_Data, (void *)T1.GetData(), size);
    memcpy((void *)(this->_Data+9), (void *)T2.GetData(), size);
    memcpy((void *)(this->_Data+18), (void *)T3.GetData(), size);
  }

  template <class T> 
  void Tensor3D3x3x3<T>::GetMatrices(Matrix3x3<T>& T1, Matrix3x3<T>& T2, 
                                     Matrix3x3<T>& T3) const
  {
    size_t size=9*sizeof(T);
    memcpy((void *)T1.GetData(), (void *)this->_Data, size);
    memcpy((void *)T2.GetData(), (void *)(this->_Data+9), size);
    memcpy((void *)T3.GetData(), (void *)(this->_Data+18), size);
  }

  template <class T> 
  void Tensor3D3x3x3<T>::LeftContravariantContraction1(Vector3<T>& v, 
                                                       Matrix3x3<T>& res)
  {
//     res.SetZero();
//     for(int i=0; i<3; i++)
//       for(int j=0; j<3; j++)
//         for(int k=0; k<3; k++)
//           res[j][k] += v[i] * _DataArray[i][j][k];
    Matrix3x3<T> T2, T3; 
    GetMatrices(res, T2, T3);
    res*=v[0];
    T2*=v[1];
    T3*=v[2];
    res+=T2;
    res+=T3;
  }

  template <class T> 
  void Tensor3D3x3x3<T>::LeftContravariantContraction2(Vector3<T>& v, 
                                                       Matrix3x3<T>& res)
  {
    res.SetZero();
    for(int i=0; i<3; i++)
      for(int j=0; j<3; j++)
        for(int k=0; k<3; k++)
          res[i][k] += v[j] * this->_DataArray[i][j][k];
  }

  template <class T> 
  void Tensor3D3x3x3<T>::LeftContravariantContraction3(Vector3<T>& v, 
                                                       Matrix3x3<T>& res)
  {
    res.SetZero();
    for(int i=0; i<3; i++)
      for(int j=0; j<3; j++)
        for(int k=0; k<3; k++)
          res[i][j] += v[k] * this->_DataArray[i][j][k];
  }


} // namespace



#endif // __Tensor3D3x3x3_hh__