]> Creatis software - creaMaracasVisu.git/blob - lib/maracasVisuLib/include/matrix.h
#3418 creaMaracasVisu Feature New Normal - ManualPaint_model with openmp
[creaMaracasVisu.git] / 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