#ifndef MATRIX3D3_HH
#define MATRIX3D3_HH
/////////////////////////////////////////////////////////////////////////
//! file="amma/Geometry/AnGeo3/Matrix3d.hh"
//! lib=Mag3
//! userlevel=Normal
//! author="Radek Marik"
//! date="26.12.1993"
//! docentry="Geometry.3-D;Basic Types.Numerical.Specialised"
//! rcsid="$Id: Matrix3d.hh,v 1.10 1999/12/02 15:47:51 ees1rr Exp $"

#include "amma/Boolean.hh"
#include "amma/StdType.hh" //RealT
#include "amma/Index.hh"
#include "amma/SBfAcc.hh"

class Point3dC;
class Vector3dC;
class MatrixC;
class ostream;
class istream;

// -------------------------------------------------------------------
// ******** Matrix3d3C ***********************************************
// -------------------------------------------------------------------

//: Matrix 3 x 3
// The class Matrix3d3C represents a matrix 3x3 of real numbers. The indexes
// start from 0.

class Matrix3d3C
{
public:
  // Type definitions.
  // =================

  enum InputDataT {ROWS, COLUMNS};
  // The kind of input data.
  
public:

    // Constructors, assigment, and destructor.
    // ========================================

    inline Matrix3d3C();
    // Creates the zero matrix.

    inline Matrix3d3C(RealT b00, RealT b01, RealT b02,
                        RealT b10, RealT b11, RealT b12,
                        RealT b20, RealT b21, RealT b22
                       );
    // Creates a new matrix and assign the values b[ij].

    Matrix3d3C(const Vector3dC & v0,
               const Vector3dC & v1,
               const Vector3dC & v2,
               const InputDataT  dType = ROWS);
    // Creates a new matrix and assign the values according to the vectors.
    // If 'dType' is equal to ROWS, the vectors are treated as rows of
    // the matrix. If 'dType' is equal to COLUMNS, the vector are treated as
    // columns of the matrix.

    Matrix3d3C(const Point3dC & p0,
               const Point3dC & p1,
               const Point3dC & p2,
               const InputDataT dType = ROWS);
    // Creates a new matrix and assign the values according to the points
    // coordinates. If 'dType' is equal to ROWS, the points are treated
    // as rows of the matrix. If 'dType' is equal to COLUMNS, the points
    // are treated as columns of the matrix.


    inline Matrix3d3C(const Matrix3d3C & mat);
    // Copy constructor.

    Matrix3d3C(const MatrixC & mat);
    // Copy constructor from the N-dimensional matrix.

    Matrix3d3C(istream & inS);
    // Creates a matrix form the input stream.

    inline const Matrix3d3C & operator=(const Matrix3d3C & mat);
    // Assigment.

    // Access to the matrix and its elements.
    // ======================================

    inline SizeT RDim() const;
    // Returns the number of rows.

    inline SizeT CDim() const;
    // Returns the number of columns.

    inline const Matrix3d3C & Matrix() const;
    // Access to the constant matrix.

    inline Matrix3d3C & Matrix();                     
    // Access to the matrix.

    inline const SizeBufferAccessC<RealT>  operator[](IndexT r) const;
    // Access to the row of the constant matrix.

    inline SizeBufferAccessC<RealT>  operator[](IndexT r);
    // Access to the row of the matrix.

    inline RealT operator()(IndexT i, IndexT j) const; 
    // Returns the value of the element a[ij].

    inline RealT & operator()(IndexT i, IndexT j); 
    // Access to the value of the element a[ij].

    inline RealT A00() const;
    // Returns the value of the element a[00].

    inline RealT A01() const;
    // Returns the value of the element a[01].

    inline RealT A02() const;
    // Returns the value of the element a[02].

    inline RealT A10() const;
    // Returns the value of the element a[10].

    inline RealT A11() const;
    // Returns the value of the element a[11].

    inline RealT A12() const;
    // Returns the value of the element a[12].

    inline RealT A20() const;
    // Returns the value of the element a[20].

    inline RealT A21() const;
    // Returns the value of the element a[21].

    inline RealT A22() const;
    // Returns the value of the element a[22].

    inline RealT & A00();
    // Access to the value of the element a[00].

    inline RealT & A01();
    // Access to the value of the element a[01].

    inline RealT & A02();
    // Access to the value of the element a[02].

    inline RealT & A10();
    // Access to the value of the element a[10].

    inline RealT & A11();
    // Access to the value of the element a[11].

    inline RealT & A12();
    // Access to the value of the element a[12].

    inline RealT & A20();
    // Access to the value of the element a[20].

    inline RealT & A21();
    // Access to the value of the element a[21].

    inline RealT & A22();
    // Access to the value of the element a[22].

    Vector3dC Column(IndexT i) const;
    // Returns the column 'i'.
    
    Vector3dC Row(IndexT i) const;
    // Return the row 'i'.

    Vector3dC Diagonal() const;
    // Returns the vector of main diagonal item.

    // Setting of matrix elements.
    // ===========================

    inline void Set(const Matrix3d3C & mat);
    // Sets all elements of the matrix according to the matrix 'mat'.

    inline Matrix3d3C & SetZero();
    // Sets all elements to be 0.0.

    Matrix3d3C & SetRow(IndexT i, const Point3dC & p);
    // Sets the elements of the row 'i' according to the elements of 
    // the point 'p'.

    Matrix3d3C & SetColumn(IndexT i, const Point3dC & p);
    // Sets the elements of the column 'i' according to the elements
    // of the point 'p'.

    Matrix3d3C & SwapRows(const IndexT r1, const IndexT r2);
    // Swaps the rows 'r1' and 'r2'.

    Matrix3d3C & SwapCols(const IndexT c1, const IndexT c2);
    // Swaps the columns 'c1' and 'c2'.

    // Logical operators.
    // ==================

    inline BooleanT operator==(const Matrix3d3C & mat) const;
    // Returns TRUE if all elements of the matrix are equal 
    // to the elements of the matrix 'mat'.

    // Arithmetical operations.
    // ========================

    inline const Matrix3d3C & operator+=(const Matrix3d3C & mat);
    // Adds the matrix 'mat' to this matrix.

    inline const Matrix3d3C & operator-=(const Matrix3d3C & mat);
    // Subtracts the matrix 'mat' from this matrix.

    inline const Matrix3d3C & operator*=(RealT lambda);
    // Multiplies all elements by the scalar 'lambda'.

    inline const Matrix3d3C & operator/=(RealT lambda);
    // Divides all elements by the scalar 'lambda'.

    inline Matrix3d3C operator+(const Matrix3d3C & mat) const;
    // Sums both matrixes and returns the result.

    inline Matrix3d3C operator-(const Matrix3d3C & mat) const;
    // Subtracts the matrix 'mat' from this matrix and returns the result.

    inline Matrix3d3C operator*(RealT lambda) const;
    // Multiplies all elements of this matrix by the scalar 'lambda' and
    // returns the result.

    inline const Matrix3d3C & AddToDiagonal(RealT a);
    // The value 'a' is added to all diagonal elements.

    inline const Matrix3d3C & SubtractFromDiagonal(RealT a);
    // The value 'a' is subtracted from all diagonal elements.

    // Matrix operators.
    // =================

    inline Matrix3d3C T() const;
    // Returns the matrix, which is the transposition of this matrix.

    inline RealT Trace() const;
    // Returns the trace of the matrix.

    inline RealT EuclideanNorm() const;
    // Returns the root of the summ of squares of the elements of the matrix.

    Matrix3d3C I() const;
    // Returns the inversion of this matrix.

    inline RealT Det() const;
    // Returns the determinant of the matrix.

    Point3dC operator*(const Point3dC & point) const;
    // Multiplies this matrix by the point 'p'. The returned point
    // is equal to (*this) * p.

    Vector3dC operator*(const Vector3dC & vector) const;
    // Multiplies the matrix by the vector 'v'. The returned vector
    // is equal to (*this) * v.

    inline Matrix3d3C operator*(const Matrix3d3C & mat) const;
    // Multiplies this matrix by the matrix 'mat'. The result matrix
    // is returned.

    Point3dC LinearSolver(const Point3dC & rightSide) const;
    // Returns the solution of 3 linear equations determined 
    // by the matrix and the right side vector 'rightSide'

    static Matrix3d3C ProductT(const Vector3dC & a);
    // Returns the matrix that is result of a*a.T().

    // Error checking.
    // ===============

    inline void ErrOutOfRange(IndexT i, const char * functionName) const;
    // Triggers the error handling if the index 'i' is smaller than 0
    // or bigger than 2. The name of the calling function 'functionName'
    // is passed to the error hangler.

  private:

    // Object representation.
    // ======================
    enum {dim = 3}; // the dimension of matrix

    RealT a[dim][dim]; // The matrix elements.
};    

ostream & operator<<(ostream & outS, const Matrix3d3C & m);
// Saves the matrix 'm' into the output stream 'outS'.

istream & operator>>(istream & inS, Matrix3d3C & matrix);
// Sets the matrix 'm' elements according to the values read from
// input stream 'inS'.

// -----------------------------------------------------------------
// ******** Matrix3dC **********************************************
// -----------------------------------------------------------------

#include "amma/Error.hh"

inline 
void 
Matrix3d3C::ErrOutOfRange(IndexT i,
                          const char * functionName) const
//========================================================
{
  if ((i < 0) || (i >= ((IntT) dim))) // Extra bracketing for Visual C++.
  {
    errAMMA << "the index out of range <0," << (dim-1) << ">: " << i;
    errAMMA.Function(functionName).Exit();
  }
}

inline 
SizeT 
Matrix3d3C::RDim() const
//======================
{ return dim; }

inline 
SizeT 
Matrix3d3C::CDim() const
//======================
{ return dim; }

inline
const SizeBufferAccessC<RealT> 
Matrix3d3C::operator[](IndexT r) const
//========================================
{
#ifdef AMMA_CHECK
    ErrOutOfRange(r, "Matrix3d3C::operator[](IndexT i) const");
#endif
  return SizeBufferAccessC<RealT>((RealT*)(a[r.V()]),CDim());
}

inline
SizeBufferAccessC<RealT> 
Matrix3d3C::operator[](IndexT r)
//================================
{
#ifdef AMMA_CHECK
    ErrOutOfRange(r, "Matrix3d3C::operator[](IndexT i) const");
#endif
  return SizeBufferAccessC<RealT>((RealT*)(a[r.V()]),CDim());
}

inline 
RealT
Matrix3d3C::operator()(IndexT i,IndexT j) const
//======================================================================
{
  #ifdef AMMA_CHECK
    ErrOutOfRange(i, "Matrix3d3C::operator()(IndexT i, IndexT) const");
    ErrOutOfRange(j, "Matrix3d3C::operator()(IndexT, IndexT j) const");
  #endif
  return a[i.V()][j.V()];
}

inline 
RealT &
Matrix3d3C::operator()(IndexT i, IndexT j)
//================================================================
{
  #ifdef AMMA_CHECK
    ErrOutOfRange(i, "Matrix3d3C::operator()(IndexT i, IndexT) const");
    ErrOutOfRange(j, "Matrix3d3C::operator()(IndexT, IndexT j) const");
  #endif
  return a[i.V()][j.V()];
}

inline 
RealT 
Matrix3d3C::A00() const
//=====================
{ return a[0][0]; }

inline 
RealT 
Matrix3d3C::A01() const
//=====================
{ return a[0][1]; }

inline 
RealT 
Matrix3d3C::A02() const
//=====================
{ return a[0][2]; }

inline 
RealT 
Matrix3d3C::A10() const
//=====================
{ return a[1][0]; }

inline 
RealT 
Matrix3d3C::A11() const
//=====================
{ return a[1][1]; }

inline 
RealT 
Matrix3d3C::A12() const
//=====================
{ return a[1][2]; }

inline 
RealT 
Matrix3d3C::A20() const
//====================
{ return a[2][0]; }

inline 
RealT 
Matrix3d3C::A21() const
//=====================
{ return a[2][1]; }

inline 
RealT 
Matrix3d3C::A22() const
//=====================
{ return a[2][2]; }

inline 
RealT & 
Matrix3d3C::A00()
//===============
{ return a[0][0]; }

inline 
RealT & 
Matrix3d3C::A01()
//===============
{ return a[0][1]; }

inline 
RealT & 
Matrix3d3C::A02()
//===============
{ return a[0][2]; }

inline 
RealT & 
Matrix3d3C::A10()
//===============
{ return a[1][0]; }

inline 
RealT & 
Matrix3d3C::A11()
//===============
{ return a[1][1]; }

inline 
RealT & 
Matrix3d3C::A12()
//===============
{ return a[1][2]; }

inline 
RealT & 
Matrix3d3C::A20()
//===============
{ return a[2][0]; }

inline 
RealT & 
Matrix3d3C::A21()
//===============
{ return a[2][1]; }

inline 
RealT & 
Matrix3d3C::A22()
//===============
{ return a[2][2]; }


inline
Matrix3d3C::Matrix3d3C(RealT b00,
                       RealT b01,
                       RealT b02,

                       RealT b10,
                       RealT b11,
                       RealT b12,

                       RealT b20,
                       RealT b21,
                       RealT b22
                      )
//============================================
{
  A00() = b00;   A01() = b01;   A02() = b02;
  A10() = b10;   A11() = b11;   A12() = b12;
  A20() = b20;   A21() = b21;   A22() = b22;
}

inline
Matrix3d3C::Matrix3d3C(const Matrix3d3C & mat)
//============================================
{
  A00() = mat.A00();  A01() = mat.A01();  A02() = mat.A02();
  A10() = mat.A10();  A11() = mat.A11();  A12() = mat.A12();
  A20() = mat.A20();  A21() = mat.A21();  A22() = mat.A22();
}

inline
const Matrix3d3C & 
Matrix3d3C::operator=(const Matrix3d3C & mat)
//===========================================
{
  A00() = mat.A00();
  A01() = mat.A01();
  A02() = mat.A02();

  A10() = mat.A10();
  A11() = mat.A11();
  A12() = mat.A12();

  A20() = mat.A20();
  A21() = mat.A21();
  A22() = mat.A22();

  return *this;
}

inline
void 
Matrix3d3C::Set(const Matrix3d3C & mat)
//=====================================
{ *this = mat; }

inline
Matrix3d3C &
Matrix3d3C::SetZero()
//===================
{
  A00() = 0.0;
  A01() = 0.0;
  A02() = 0.0;

  A10() = 0.0;
  A11() = 0.0;
  A12() = 0.0;

  A20() = 0.0;
  A21() = 0.0;
  A22() = 0.0;

  return *this;
}

inline
Matrix3d3C::Matrix3d3C()
//======================
{ SetZero(); }

inline 
const Matrix3d3C &
Matrix3d3C::Matrix() const
//========================
{ return *this; }

inline 
Matrix3d3C & 
Matrix3d3C::Matrix()                     
//==================
{ return *this; }

inline 
BooleanT 
Matrix3d3C::operator==(const Matrix3d3C & mat) const
//==================================================
{
  return
      (A00() == mat.A00())&&(A01() == mat.A01())&&(A02() == mat.A02()) 
    &&(A10() == mat.A10())&&(A11() == mat.A11())&&(A12() == mat.A12())
    &&(A20() == mat.A20())&&(A21() == mat.A21())&&(A22() == mat.A22());
}
  
inline 
const Matrix3d3C & 
Matrix3d3C::operator+=(const Matrix3d3C & mat)
//============================================
{
  A00() += mat.A00();  A01() += mat.A01();  A02() += mat.A02();
  A10() += mat.A10();  A11() += mat.A11();  A12() += mat.A12(); 
  A20() += mat.A20();  A21() += mat.A21();  A22() += mat.A22();
  return *this;
}

inline 
const Matrix3d3C & 
Matrix3d3C::operator-=(const Matrix3d3C & mat)
//============================================
{
  A00() -= mat.A00();  A01() -= mat.A01();  A02() -= mat.A02();
  A10() -= mat.A10();  A11() -= mat.A11();  A12() -= mat.A12(); 
  A20() -= mat.A20();  A21() -= mat.A21();  A22() -= mat.A22();
  return *this;
}

inline 
const Matrix3d3C & 
Matrix3d3C::operator*=(RealT lambda)
//==================================
{
  A00() *= lambda;  A01() *= lambda;  A02() *= lambda;
  A10() *= lambda;  A11() *= lambda;  A12() *= lambda; 
  A20() *= lambda;  A21() *= lambda;  A22() *= lambda;
  return *this;
}

inline 
const Matrix3d3C & 
Matrix3d3C::operator/=(RealT lambda)
//==================================
{
  A00() /= lambda;  A01() /= lambda;  A02() /= lambda;
  A10() /= lambda;  A11() /= lambda;  A12() /= lambda; 
  A20() /= lambda;  A21() /= lambda;  A22() /= lambda;
  return *this;
}

inline 
Matrix3d3C 
Matrix3d3C::operator+(const Matrix3d3C & mat) const
//=================================================
{
  return Matrix3d3C(A00()+mat.A00(), A01()+mat.A01(), A02()+mat.A02(),
                    A10()+mat.A10(), A11()+mat.A11(), A12()+mat.A12(),
                    A20()+mat.A20(), A21()+mat.A21(), A22()+mat.A22()
                   );
}

inline 
Matrix3d3C 
Matrix3d3C::operator-(const Matrix3d3C & mat) const
//=================================================
{
  return Matrix3d3C(A00()-mat.A00(), A01()-mat.A01(), A02()-mat.A02(),
                    A10()-mat.A10(), A11()-mat.A11(), A12()-mat.A12(),
                    A20()-mat.A20(), A21()-mat.A21(), A22()-mat.A22()
                   );
}

inline 
Matrix3d3C 
Matrix3d3C::operator*(RealT lambda) const
//=======================================
{
  return Matrix3d3C(A00()*lambda, A01()*lambda, A02()*lambda,
                    A10()*lambda, A11()*lambda, A12()*lambda,
                    A20()*lambda, A21()*lambda, A22()*lambda
                   );
}

inline 
Matrix3d3C 
Matrix3d3C::operator*(const Matrix3d3C & mat) const
//=================================================
{
  return Matrix3d3C(A00()*mat.A00()+A01()*mat.A10()+A02()*mat.A20(),
                    A00()*mat.A01()+A01()*mat.A11()+A02()*mat.A21(),
                    A00()*mat.A02()+A01()*mat.A12()+A02()*mat.A22(),
 
                    A10()*mat.A00()+A11()*mat.A10()+A12()*mat.A20(),
                    A10()*mat.A01()+A11()*mat.A11()+A12()*mat.A21(),
                    A10()*mat.A02()+A11()*mat.A12()+A12()*mat.A22(),
 
                    A20()*mat.A00()+A21()*mat.A10()+A22()*mat.A20(),
                    A20()*mat.A01()+A21()*mat.A11()+A22()*mat.A21(),
                    A20()*mat.A02()+A21()*mat.A12()+A22()*mat.A22()
                   );
}

inline 
const Matrix3d3C & 
Matrix3d3C::AddToDiagonal(RealT a)
//================================
{
  A00() += a;
  A11() += a;
  A22() += a;
  return *this;
}

inline 
const Matrix3d3C &
Matrix3d3C::SubtractFromDiagonal(RealT a)
//=======================================
{
  A00() -= a;
  A11() -= a;
  A22() -= a;
  return *this;
}

inline 
Matrix3d3C 
Matrix3d3C::T() const
//===================
{
  return Matrix3d3C(A00(), A10(), A20(),
                    A01(), A11(), A21(),
                    A02(), A12(), A22()
                   );
}

inline 
RealT 
Matrix3d3C::Trace() const
//=======================
{ return A00()+A11()+A22(); }

inline 
RealT 
Matrix3d3C::EuclideanNorm() const
//=======================
{  return sqrt(A00()*A00() + A10()*A10() + A20()*A20() +
                    A01()*A01() + A11()*A11() + A21()*A21() +
                    A02()*A02() + A12()*A12() + A22()*A22());
}

inline 
RealT
Matrix3d3C::Det() const
//=====================
{
  return    A00() * (A11()*A22() - A12()*A21())
          + A10() * (A21()*A02() - A22()*A01())
          + A20() * (A01()*A12() - A02()*A11());
}

#endif

// IAPS - Image analysis program system. 
// End of include file Matrix3d3.hh