/*# ---------------------------------------------------------------------
#
# Copyright (c) CREATIS (Centre de Recherche en Acquisition et Traitement de l'Image
#                        pour la Sant�)
# Authors : Eduardo Davila, Frederic Cervenansky, Claire Mouton
# Previous Authors : Laurent Guigues, Jean-Pierre Roux
# CreaTools website : www.creatis.insa-lyon.fr/site/fr/creatools_accueil
#
#  This software is governed by the CeCILL-B license under French law and
#  abiding by the rules of distribution of free software. You can  use,
#  modify and/ or redistribute the software under the terms of the CeCILL-B
#  license as circulated by CEA, CNRS and INRIA at the following URL
#  http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
#  or in the file LICENSE.txt.
#
#  As a counterpart to the access to the source code and  rights to copy,
#  modify and redistribute granted by the license, users are provided only
#  with a limited warranty  and the software's author,  the holder of the
#  economic rights,  and the successive licensors  have only  limited
#  liability.
#
#  The fact that you are presently reading this means that you have had
#  knowledge of the CeCILL-B license and that you accept its terms.
# ------------------------------------------------------------------------ */

////////////////////////////////////////////////////////////////////////////////
// matrix.h
// Creation : 19/03/2000
// Author   : Leonardo FLOREZ VALENCIA
//               l-florez@uniandes.edu.co
//               lflorez@creatis.insa-lyon.fr
// Copyright (C) 2000-2002 Leonardo FLOREZ VALENCIA
//
// This program is free software; you can redistribute it and/or
// modify it under the terms of the GNU General Public License
// as published by the Free Software Foundation; either version 2
// of the License, or (at your option) any later version.
//
// This program 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 General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
////////////////////////////////////////////////////////////////////////////////

#ifndef GTMLIB__MATH__MATRIX__HXX
#define GTMLIB__MATH__MATRIX__HXX

#include "mathdefs.h"
#include "vmfuncs.h"
#include "vector.h"

namespace gtm
{
    /** TMatrix class.
     *
     *  This class defines a C++ template to use mathematical matrices of NxM.
     *  @see TVector
     */
    template< class T >
    class TMatrix
    {
    public:
	
	/** Constructors.
	 *
	 *  @param N Columns.
	 *  @param M Rows.
	 *  @param data Initial data.
	 *  @param r Copy object (matrix or vector).
	 *  @param block Memory block.
	 */
	//@{
	/// Default constructor.
	TMatrix( uint32_t N = 3, uint32_t M = 3, T data = ( T )0 );
	/// Copy constructor.
	TMatrix( const TMatrix< T >& r );
	/// ANSI casting constructor.
	TMatrix( T** block, uint32_t N, uint32_t M );
	//@}

	/// Destructor.
	~TMatrix( ) {
	    MatrixFreeMemory< T >( _matrix, _N );
	};

	/** Assignation operators.
	 *
	 *  @param r Right object (matrix, vector or scalar).
	 */
	//@{
	/// Natural assignation.
	TMatrix< T >& operator=( const TMatrix< T >& r );
	/// Vector assignation.
	TMatrix< T >& operator=( TVector< T >& r );
	/// Scalar assignation.
	TMatrix< T >& operator=( T r );
	//@}
	
	/** Comparation operators.
	 *
	 *  @param Right matrix.
	 */
	//@{
	/// Equality.
	bool operator==( const TMatrix< T >& r );
	/// Inequality.
	bool operator!=( const TMatrix< T >& r );
	//@}

	/// Reference operator.
	T& operator()( uint32_t i, uint32_t j ) {
	    return( _matrix[ i ][ j ] );
	};
	/// Columns
	uint32_t GetN( ) {
	    return( _N );
	};
	/// Rows
	uint32_t GetM( ) {
	    return( _M );
	};
	/// Returns the ANSI (C/C++) reference.
	T** GetAnsiRef( ) {
	    return( _matrix );
	};
	/// Determinant.
	T Det( ) {
	    return( MatrixDeterminant< T >( _matrix, _N ) );
	};
	/// Loads identity.
	void LoadIdentity( ) {
	    MatrixLoadIdentity< T >( _matrix, GTM_MIN( _N, _M ) );
	};

	/** Binary operators.
	 *
	 *  @param r Right objet (matrix, vector or scalar).
	 */
	//@{
	/// Addition.
	TMatrix< T > operator+( const TMatrix< T >& r );
	/// Substraction.
	TMatrix< T > operator-( const TMatrix< T >& r );
	/// Product.
	TMatrix< T > operator*( const TMatrix< T >& r );
	/// Product (vector).
	TMatrix< T > operator*( TVector< T >& r );
	/// Scalar product
	TMatrix< T > operator*( T r );
	//@}

	/** Self-assignation binary operators.
	 *
	 *  @param r Right object (matrix, vector or scalar).
	 */
	//@{
	/// Addition.
	TMatrix< T >& operator+=( const TMatrix< T >& r ) {
	    *this = *this + r;
	    return( *this );
	};
	/// Substraction.
	TMatrix< T >& operator-=( const TMatrix< T >& r ) {
	    *this = *this - r;
	    return( *this );
	};
	/// Product.
	TMatrix< T >& operator*=( const TMatrix< T >& r ) {
	    *this = *this * r;
	    return( *this );
	};
	/// Product (vector).
	TMatrix< T >& operator*=( TVector< T >& r ) {
	    *this = *this * r;
	    return( *this );
	};
	/// Scalar product.
	TMatrix< T >& operator*=( T r ) {
	    *this = *this * r;
	    return( *this );
	};
	//@}

	/** Unary operators.
	 */
	//@{
	/// Additive inverse.
	TMatrix< T > operator-( );
	/// Matrix inverse.
	TMatrix< T > operator!( );
	/// Matrix transpose.
	TMatrix< T > operator~( );
	//@}

    private:

	/// Matrix' internal state.
	//@{
	/// Memory block.
	T** _matrix;
	/// Columns.
	uint32_t _N;
	/// Rows.
	uint32_t _M;
	//@}

    };

// -----------------------------------------------------------------------------
    template< class T >
    TVector< T >& TVector< T >::operator=( TMatrix< T >& r )
    {
	uint32_t i, j, k, min;

	// This min calc. avoids to reserve temporary memory, so, be careful.
	min = GTM_MIN( r.GetN( ) * r.GetM( ), _N );
	_type = ( r.GetN( ) == 1 )? COL_VECTOR: ROW_VECTOR;
	for( i = 0, k = 0; i < r.GetN( ) && k < min; i++ )
	    for( j = 0; j < r.GetM( ) && k < min; _vector[ k++ ] = r( i, j ), j++ );
	return( *this );

    }

// -----------------------------------------------------------------------------
    template< class T >
    TMatrix< T > TVector< T >::operator*( TMatrix< T >& r )
    {
	TMatrix< T > m = *this;
	return( m * r );

    }

// -----------------------------------------------------------------------------
    template< class T >
    TMatrix< T >::TMatrix( uint32_t N, uint32_t M, T data )
    {
	_N = N;
	_M = M;
	_matrix = MatrixAllocateMemory< T >( _N, _M );
	MatrixAssignScalar< T >( _matrix, data, _N, _M );

    }

// -----------------------------------------------------------------------------
    template< class T >
    TMatrix< T >::TMatrix( const TMatrix< T >& r )
    {
	_N = r._N;
	_M = r._M;
	_matrix = MatrixCopyMemory< T >( r._matrix, _N, _M );

    }

// -----------------------------------------------------------------------------
    template< class T >
    TMatrix< T >::TMatrix( T** block, uint32_t N, uint32_t M )
    {
	_N = N;
	_M = M;
	_matrix = MatrixCopyMemory< T >( block, N, M );

    }

// -----------------------------------------------------------------------------
    template< class T >
    TMatrix< T >& TMatrix< T >::operator=( const TMatrix< T >& r )
    {
	if( _N != r._N || _M != r._M ) {

	    MatrixFreeMemory< T >( _matrix, _N );
	    _N = r._N;
	    _M = r._M;
	    _matrix = MatrixCopyMemory< T >( r._matrix, _N, _M );

	} else MatrixAssignMemory< T >( _matrix, r._matrix, _N, _M );
	return( *this );
  
    }

// -----------------------------------------------------------------------------
    template< class T >
    TMatrix< T >& TMatrix< T >::operator=( TVector< T >& r )
    {
	uint32_t i;
	uint32_t n = r.GetN( );
	bool column = ( r.GetType( ) == COL_VECTOR );

	MatrixFreeMemory< T >( _matrix, _N );
	_N = ( column )? 1: n;
	_M = ( column )? n: 1;
	_matrix = MatrixAllocateMemory< T >( _N, _M );
	for( i = 0; i < n; _matrix[ ( column )? 0: i ][ ( column )? i: 0 ] = r( i ), i++ );
	return( *this );
  
    }

// -----------------------------------------------------------------------------
    template< class T >
    TMatrix< T >& TMatrix< T >::operator=( T r )
    {
	MatrixAssignScalar< T >( _matrix, r, _N, _M );
	return( *this );

    }

// -----------------------------------------------------------------------------
    template< class T >
    bool TMatrix< T >::operator==( const TMatrix< T >& r )
    {
	uint32_t i, j;
	bool ret;

	for(
	    i = 0, ret = ( _N == r._N && _M == r._M );
	    i < _N && ret;
	    i++
	    ) for(
		j = 0;
		j < _M && ret;
		ret &= ( _matrix[ i ][ j ] == r._matrix[ i ][ j ] ), j++
		);
	return( ret );

    }

// -----------------------------------------------------------------------------
    template< class T >
    bool TMatrix< T >::operator!=( const TMatrix< T >& r )
    {
	uint32_t i, j;
	bool ret;

	for(
	    i = 0, ret = ( _N != r._N || _M != r._M );
	    i < _N && !ret;
	    i++
	    ) for(
		j = 0;
		j < _M && !ret;
		ret |= ( _matrix[ i ][ j ] != r._matrix[ i ][ j ] ), j++
		);
	return( ret );

    }

// -----------------------------------------------------------------------------
    template< class T >
    TMatrix< T > TMatrix< T >::operator+( const TMatrix< T >& r )
    {
	TMatrix< T > ret( _N, _M );

	MatrixAdd< T >(
	    ret._matrix,
	    _matrix,
	    r._matrix,
	    GTM_MIN( _N, r._N ),
	    GTM_MIN( _M, r._M )
	    );
	return( ret );

    }

// -----------------------------------------------------------------------------
    template< class T >
    TMatrix< T > TMatrix< T >::operator-( const TMatrix< T >& r )
    {
	TMatrix< T > ret( _N, _M );
  
	MatrixSubtract< T >(
	    ret._matrix,
	    _matrix,
	    r._matrix,
	    GTM_MIN( _N, r._N ),
	    GTM_MIN( _M, r._M )
	    );
	return( ret );

    }

// -----------------------------------------------------------------------------
    template< class T >
    TMatrix< T > TMatrix< T >::operator*( const TMatrix< T >& r )
    {
	TMatrix< T > ret( r._N, _M );

	MatrixProduct< T >( ret._matrix, _matrix, r._matrix, _N, _M, r._N );
	return( ret );

    }

// -----------------------------------------------------------------------------
    template< class T >
    TMatrix< T > TMatrix< T >::operator*( T r )
    {
	TMatrix< T > ret( _N, _M );

	MatrixScalarProduct< T >( ret._matrix, _matrix, r, _N, _M );
	return( ret );

    }

// -----------------------------------------------------------------------------
    template< class T >
    TMatrix< T > TMatrix< T >::operator*( TVector< T >& r )
    {
	TMatrix< T > m;
	m = r;
	return( *this * m );

    }

// -----------------------------------------------------------------------------
    template< class T >
    TMatrix< T > TMatrix< T >::operator-( )
    {
	TMatrix< T > ret( _N, _M );
	uint32_t i, j;
  
	for( i = 0; i < _N; i++ )
	    for( j = 0; j < _M; ret._matrix[ i ][ j ] = ( T )0 - _matrix[ i ][ j ], j++ );
	return( ret );
  
    }

// -----------------------------------------------------------------------------
    template< class T >
    TMatrix< T > TMatrix< T >::operator!( )
    {
	TMatrix< T > ret( _N, _N );

	if( _N <= 4 ) MatrixInverseAdjoint< T >( ret._matrix, _matrix, _N );
	else          MatrixInverseGauss< T >( ret._matrix, _matrix, _N );
	return( ret );

    }

// -----------------------------------------------------------------------------
    template< class T >
    TMatrix< T > TMatrix< T >::operator~( )
    {
	TMatrix< T > ret( _M, _N );

	MatrixTranspose< T >( ret._matrix, _matrix, _N, _M );
	return( ret );

    }

} // namespace

#endif // GTMLIB__MATH__MATRIX__HXX

// EOF - matrix.h