lib/maracasVisuLib/include/matrix.h

   1 /*# ---------------------------------------------------------------------
   2 #
   3 # Copyright (c) CREATIS (Centre de Recherche en Acquisition et Traitement de l'Image
   4 #                        pour la Sant�)
   5 # Authors : Eduardo Davila, Frederic Cervenansky, Claire Mouton
   6 # Previous Authors : Laurent Guigues, Jean-Pierre Roux
   7 # CreaTools website : www.creatis.insa-lyon.fr/site/fr/creatools_accueil
   8 #
   9 #  This software is governed by the CeCILL-B license under French law and
  10 #  abiding by the rules of distribution of free software. You can  use,
  11 #  modify and/ or redistribute the software under the terms of the CeCILL-B
  12 #  license as circulated by CEA, CNRS and INRIA at the following URL
  13 #  http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
  14 #  or in the file LICENSE.txt.
  15 #
  16 #  As a counterpart to the access to the source code and  rights to copy,
  17 #  modify and redistribute granted by the license, users are provided only
  18 #  with a limited warranty  and the software's author,  the holder of the
  19 #  economic rights,  and the successive licensors  have only  limited
  20 #  liability.
  21 #
  22 #  The fact that you are presently reading this means that you have had
  23 #  knowledge of the CeCILL-B license and that you accept its terms.
  24 # ------------------------------------------------------------------------ */
  25
  26 ////////////////////////////////////////////////////////////////////////////////
  27 // matrix.h
  28 // Creation : 19/03/2000
  29 // Author   : Leonardo FLOREZ VALENCIA
  30 //               l-florez@uniandes.edu.co
  31 //               lflorez@creatis.insa-lyon.fr
  32 // Copyright (C) 2000-2002 Leonardo FLOREZ VALENCIA
  33 //
  34 // This program is free software; you can redistribute it and/or
  35 // modify it under the terms of the GNU General Public License
  36 // as published by the Free Software Foundation; either version 2
  37 // of the License, or (at your option) any later version.
  38 //
  39 // This program is distributed in the hope that it will be useful,
  40 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  41 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  42 // GNU General Public License for more details.
  43 //
  44 // You should have received a copy of the GNU General Public License
  45 // along with this program; if not, write to the Free Software
  46 // Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
  47 ////////////////////////////////////////////////////////////////////////////////
  48
  49 #ifndef GTMLIB__MATH__MATRIX__HXX
  50 #define GTMLIB__MATH__MATRIX__HXX
  51
  52 #include "mathdefs.h"
  53 #include "vmfuncs.h"
  54 #include "vector.h"
  55
  56 namespace gtm
  57 {
  58     /** TMatrix class.
  59      *
  60      *  This class defines a C++ template to use mathematical matrices of NxM.
  61      *  @see TVector
  62      */
  63     template< class T >
  64     class TMatrix
  65     {
  66     public:
  67
  68         /** Constructors.
  69          *
  70          *  @param N Columns.
  71          *  @param M Rows.
  72          *  @param data Initial data.
  73          *  @param r Copy object (matrix or vector).
  74          *  @param block Memory block.
  75          */
  76         //@{
  77         /// Default constructor.
  78         TMatrix( uint32_t N = 3, uint32_t M = 3, T data = ( T )0 );
  79         /// Copy constructor.
  80         TMatrix( const TMatrix< T >& r );
  81         /// ANSI casting constructor.
  82         TMatrix( T** block, uint32_t N, uint32_t M );
  83         //@}
  84
  85         /// Destructor.
  86         ~TMatrix( ) {
  87             MatrixFreeMemory< T >( _matrix, _N );
  88         };
  89
  90         /** Assignation operators.
  91          *
  92          *  @param r Right object (matrix, vector or scalar).
  93          */
  94         //@{
  95         /// Natural assignation.
  96         TMatrix< T >& operator=( const TMatrix< T >& r );
  97         /// Vector assignation.
  98         TMatrix< T >& operator=( TVector< T >& r );
  99         /// Scalar assignation.
 100         TMatrix< T >& operator=( T r );
 101         //@}
 102
 103         /** Comparation operators.
 104          *
 105          *  @param Right matrix.
 106          */
 107         //@{
 108         /// Equality.
 109         bool operator==( const TMatrix< T >& r );
 110         /// Inequality.
 111         bool operator!=( const TMatrix< T >& r );
 112         //@}
 113
 114         /// Reference operator.
 115         T& operator()( uint32_t i, uint32_t j ) {
 116             return( _matrix[ i ][ j ] );
 117         };
 118         /// Columns
 119         uint32_t GetN( ) {
 120             return( _N );
 121         };
 122         /// Rows
 123         uint32_t GetM( ) {
 124             return( _M );
 125         };
 126         /// Returns the ANSI (C/C++) reference.
 127         T** GetAnsiRef( ) {
 128             return( _matrix );
 129         };
 130         /// Determinant.
 131         T Det( ) {
 132             return( MatrixDeterminant< T >( _matrix, _N ) );
 133         };
 134         /// Loads identity.
 135         void LoadIdentity( ) {
 136             MatrixLoadIdentity< T >( _matrix, GTM_MIN( _N, _M ) );
 137         };
 138
 139         /** Binary operators.
 140          *
 141          *  @param r Right objet (matrix, vector or scalar).
 142          */
 143         //@{
 144         /// Addition.
 145         TMatrix< T > operator+( const TMatrix< T >& r );
 146         /// Substraction.
 147         TMatrix< T > operator-( const TMatrix< T >& r );
 148         /// Product.
 149         TMatrix< T > operator*( const TMatrix< T >& r );
 150         /// Product (vector).
 151         TMatrix< T > operator*( TVector< T >& r );
 152         /// Scalar product
 153         TMatrix< T > operator*( T r );
 154         //@}
 155
 156         /** Self-assignation binary operators.
 157          *
 158          *  @param r Right object (matrix, vector or scalar).
 159          */
 160         //@{
 161         /// Addition.
 162         TMatrix< T >& operator+=( const TMatrix< T >& r ) {
 163             *this = *this + r;
 164             return( *this );
 165         };
 166         /// Substraction.
 167         TMatrix< T >& operator-=( const TMatrix< T >& r ) {
 168             *this = *this - r;
 169             return( *this );
 170         };
 171         /// Product.
 172         TMatrix< T >& operator*=( const TMatrix< T >& r ) {
 173             *this = *this * r;
 174             return( *this );
 175         };
 176         /// Product (vector).
 177         TMatrix< T >& operator*=( TVector< T >& r ) {
 178             *this = *this * r;
 179             return( *this );
 180         };
 181         /// Scalar product.
 182         TMatrix< T >& operator*=( T r ) {
 183             *this = *this * r;
 184             return( *this );
 185         };
 186         //@}
 187
 188         /** Unary operators.
 189          */
 190         //@{
 191         /// Additive inverse.
 192         TMatrix< T > operator-( );
 193         /// Matrix inverse.
 194         TMatrix< T > operator!( );
 195         /// Matrix transpose.
 196         TMatrix< T > operator~( );
 197         //@}
 198
 199     private:
 200
 201         /// Matrix' internal state.
 202         //@{
 203         /// Memory block.
 204         T** _matrix;
 205         /// Columns.
 206         uint32_t _N;
 207         /// Rows.
 208         uint32_t _M;
 209         //@}
 210
 211     };
 212
 213 // -----------------------------------------------------------------------------
 214     template< class T >
 215     TVector< T >& TVector< T >::operator=( TMatrix< T >& r )
 216     {
 217         uint32_t i, j, k, min;
 218
 219         // This min calc. avoids to reserve temporary memory, so, be careful.
 220         min = GTM_MIN( r.GetN( ) * r.GetM( ), _N );
 221         _type = ( r.GetN( ) == 1 )? COL_VECTOR: ROW_VECTOR;
 222         for( i = 0, k = 0; i < r.GetN( ) && k < min; i++ )
 223             for( j = 0; j < r.GetM( ) && k < min; _vector[ k++ ] = r( i, j ), j++ );
 224         return( *this );
 225
 226     }
 227
 228 // -----------------------------------------------------------------------------
 229     template< class T >
 230     TMatrix< T > TVector< T >::operator*( TMatrix< T >& r )
 231     {
 232         TMatrix< T > m = *this;
 233         return( m * r );
 234
 235     }
 236
 237 // -----------------------------------------------------------------------------
 238     template< class T >
 239     TMatrix< T >::TMatrix( uint32_t N, uint32_t M, T data )
 240     {
 241         _N = N;
 242         _M = M;
 243         _matrix = MatrixAllocateMemory< T >( _N, _M );
 244         MatrixAssignScalar< T >( _matrix, data, _N, _M );
 245
 246     }
 247
 248 // -----------------------------------------------------------------------------
 249     template< class T >
 250     TMatrix< T >::TMatrix( const TMatrix< T >& r )
 251     {
 252         _N = r._N;
 253         _M = r._M;
 254         _matrix = MatrixCopyMemory< T >( r._matrix, _N, _M );
 255
 256     }
 257
 258 // -----------------------------------------------------------------------------
 259     template< class T >
 260     TMatrix< T >::TMatrix( T** block, uint32_t N, uint32_t M )
 261     {
 262         _N = N;
 263         _M = M;
 264         _matrix = MatrixCopyMemory< T >( block, N, M );
 265
 266     }
 267
 268 // -----------------------------------------------------------------------------
 269     template< class T >
 270     TMatrix< T >& TMatrix< T >::operator=( const TMatrix< T >& r )
 271     {
 272         if( _N != r._N || _M != r._M ) {
 273
 274             MatrixFreeMemory< T >( _matrix, _N );
 275             _N = r._N;
 276             _M = r._M;
 277             _matrix = MatrixCopyMemory< T >( r._matrix, _N, _M );
 278
 279         } else MatrixAssignMemory< T >( _matrix, r._matrix, _N, _M );
 280         return( *this );
 281
 282     }
 283
 284 // -----------------------------------------------------------------------------
 285     template< class T >
 286     TMatrix< T >& TMatrix< T >::operator=( TVector< T >& r )
 287     {
 288         uint32_t i;
 289         uint32_t n = r.GetN( );
 290         bool column = ( r.GetType( ) == COL_VECTOR );
 291
 292         MatrixFreeMemory< T >( _matrix, _N );
 293         _N = ( column )? 1: n;
 294         _M = ( column )? n: 1;
 295         _matrix = MatrixAllocateMemory< T >( _N, _M );
 296         for( i = 0; i < n; _matrix[ ( column )? 0: i ][ ( column )? i: 0 ] = r( i ), i++ );
 297         return( *this );
 298
 299     }
 300
 301 // -----------------------------------------------------------------------------
 302     template< class T >
 303     TMatrix< T >& TMatrix< T >::operator=( T r )
 304     {
 305         MatrixAssignScalar< T >( _matrix, r, _N, _M );
 306         return( *this );
 307
 308     }
 309
 310 // -----------------------------------------------------------------------------
 311     template< class T >
 312     bool TMatrix< T >::operator==( const TMatrix< T >& r )
 313     {
 314         uint32_t i, j;
 315         bool ret;
 316
 317         for(
 318             i = 0, ret = ( _N == r._N && _M == r._M );
 319             i < _N && ret;
 320             i++
 321             ) for(
 322                 j = 0;
 323                 j < _M && ret;
 324                 ret &= ( _matrix[ i ][ j ] == r._matrix[ i ][ j ] ), j++
 325                 );
 326         return( ret );
 327
 328     }
 329
 330 // -----------------------------------------------------------------------------
 331     template< class T >
 332     bool TMatrix< T >::operator!=( const TMatrix< T >& r )
 333     {
 334         uint32_t i, j;
 335         bool ret;
 336
 337         for(
 338             i = 0, ret = ( _N != r._N || _M != r._M );
 339             i < _N && !ret;
 340             i++
 341             ) for(
 342                 j = 0;
 343                 j < _M && !ret;
 344                 ret |= ( _matrix[ i ][ j ] != r._matrix[ i ][ j ] ), j++
 345                 );
 346         return( ret );
 347
 348     }
 349
 350 // -----------------------------------------------------------------------------
 351     template< class T >
 352     TMatrix< T > TMatrix< T >::operator+( const TMatrix< T >& r )
 353     {
 354         TMatrix< T > ret( _N, _M );
 355
 356         MatrixAdd< T >(
 357             ret._matrix,
 358             _matrix,
 359             r._matrix,
 360             GTM_MIN( _N, r._N ),
 361             GTM_MIN( _M, r._M )
 362             );
 363         return( ret );
 364
 365     }
 366
 367 // -----------------------------------------------------------------------------
 368     template< class T >
 369     TMatrix< T > TMatrix< T >::operator-( const TMatrix< T >& r )
 370     {
 371         TMatrix< T > ret( _N, _M );
 372
 373         MatrixSubtract< T >(
 374             ret._matrix,
 375             _matrix,
 376             r._matrix,
 377             GTM_MIN( _N, r._N ),
 378             GTM_MIN( _M, r._M )
 379             );
 380         return( ret );
 381
 382     }
 383
 384 // -----------------------------------------------------------------------------
 385     template< class T >
 386     TMatrix< T > TMatrix< T >::operator*( const TMatrix< T >& r )
 387     {
 388         TMatrix< T > ret( r._N, _M );
 389
 390         MatrixProduct< T >( ret._matrix, _matrix, r._matrix, _N, _M, r._N );
 391         return( ret );
 392
 393     }
 394
 395 // -----------------------------------------------------------------------------
 396     template< class T >
 397     TMatrix< T > TMatrix< T >::operator*( T r )
 398     {
 399         TMatrix< T > ret( _N, _M );
 400
 401         MatrixScalarProduct< T >( ret._matrix, _matrix, r, _N, _M );
 402         return( ret );
 403
 404     }
 405
 406 // -----------------------------------------------------------------------------
 407     template< class T >
 408     TMatrix< T > TMatrix< T >::operator*( TVector< T >& r )
 409     {
 410         TMatrix< T > m;
 411         m = r;
 412         return( *this * m );
 413
 414     }
 415
 416 // -----------------------------------------------------------------------------
 417     template< class T >
 418     TMatrix< T > TMatrix< T >::operator-( )
 419     {
 420         TMatrix< T > ret( _N, _M );
 421         uint32_t i, j;
 422
 423         for( i = 0; i < _N; i++ )
 424             for( j = 0; j < _M; ret._matrix[ i ][ j ] = ( T )0 - _matrix[ i ][ j ], j++ );
 425         return( ret );
 426
 427     }
 428
 429 // -----------------------------------------------------------------------------
 430     template< class T >
 431     TMatrix< T > TMatrix< T >::operator!( )
 432     {
 433         TMatrix< T > ret( _N, _N );
 434
 435         if( _N <= 4 ) MatrixInverseAdjoint< T >( ret._matrix, _matrix, _N );
 436         else          MatrixInverseGauss< T >( ret._matrix, _matrix, _N );
 437         return( ret );
 438
 439     }
 440
 441 // -----------------------------------------------------------------------------
 442     template< class T >
 443     TMatrix< T > TMatrix< T >::operator~( )
 444     {
 445         TMatrix< T > ret( _M, _N );
 446
 447         MatrixTranspose< T >( ret._matrix, _matrix, _N, _M );
 448         return( ret );
 449
 450     }
 451
 452 } // namespace
 453
 454 #endif // GTMLIB__MATH__MATRIX__HXX
 455
 456 // EOF - matrix.h