lib/maracasVisuLib/include/vmfuncs.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 // vmfuncs.h
  28 // Creation : 02/02/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__VMFUNCS__HXX
  50 #define GTMLIB__MATH__VMFUNCS__HXX
  51
  52 #include <iostream>
  53 #include <math.h>
  54 #include <memory.h>
  55 #include "mathdefs.h"
  56
  57 #include <creaSystem.h>
  58
  59 using namespace std;
  60
  61 namespace gtm
  62 {
  63 // -----------------------------------------------------------------------------
  64 //
  65 /// Functions for vectors.
  66 //@{
  67 // -----------------------------------------------------------------------------
  68     /** Stream out function for vectors.
  69      *
  70      *  @param o Output stream.
  71      *  @param v ANSI array.
  72      *  @param n Cardinality.
  73      *  @param col Should I print it in column format?
  74      */
  75     template< class T >
  76     inline
  77     void VectorPrint( std::ostream& o, T* v, uint32_t n, bool col )
  78     {
  79         uint32_t i;
  80
  81         o << n << ": [";
  82         if( col ) o << endl;
  83         else      o << " ";
  84         for( uint32_t i = 0; i < n; i++ ) {
  85
  86             o << v[ i ];
  87             if( col ) o << endl;
  88             else      o << " ";
  89
  90         } // rof
  91         o << "]" << endl;
  92
  93     }
  94
  95 // -----------------------------------------------------------------------------
  96     /** Allocation function.
  97      *
  98      *  @param n Cardinality.
  99      */
 100     template< class T >
 101     inline
 102     T* VectorAllocateMemory( uint32_t n )
 103     {
 104         T* v = new T[ n ];
 105
 106         memset( v, 0, sizeof( T ) * n );
 107         return( v );
 108
 109     }
 110
 111 // -----------------------------------------------------------------------------
 112     /** Frees vector's memory.
 113      *
 114      *  @param v Vector.
 115      */
 116     template< class T >
 117     inline
 118     void VectorFreeMemory( T* v )
 119     {
 120         if( v ) delete [] v;
 121
 122     }
 123
 124 // -----------------------------------------------------------------------------
 125     /** Copy.
 126      *
 127      *  @param v Target vector.
 128      *  @param v_o Source vector.
 129      *  @param n Cardinality.
 130      */
 131     template< class T >
 132     inline
 133     void VectorAssignMemory( T* v, T* v_o, uint32_t n )
 134     {
 135         memcpy( v, v_o, sizeof( T ) * n );
 136
 137     }
 138
 139 // -----------------------------------------------------------------------------
 140     /** Allocation & copy.
 141      *
 142      *  @param v Vector to be copied.
 143      *  @param n Cardinality.
 144      */
 145     template< class T >
 146     inline
 147     T* VectorCopyMemory( T* v, uint32_t n )
 148     {
 149         T* n_v = VectorAllocateMemory< T >( n );
 150         VectorAssignMemory< T >( n_v, v, n );
 151         return( n_v );
 152
 153     }
 154
 155 // -----------------------------------------------------------------------------
 156     /** Scalar assignation.
 157      *
 158      *  @param v Vector.
 159      *  @param s Scalar value.
 160      *  @param n Cardinality.
 161      */
 162     template< class T >
 163     inline
 164     void VectorAssignScalar( T* v, T s, uint32_t n )
 165     {
 166         uint32_t i;
 167
 168         for( i = 0; i < n; v[ i ] = s, i++ );
 169
 170     }
 171
 172 // -----------------------------------------------------------------------------
 173     /** Convert a vector in a matrix.
 174      *
 175      *  @param ma Matrix.
 176      *  @param v Vector.
 177      *  @param n Columns.
 178      *  @param m Rows.
 179      *  @param col Is it a column vector?
 180      */
 181     template< class T >
 182     inline
 183     void VectorMatrixCast( T** ma, T* v, uint32_t n, uint32_t m, bool col )
 184     {
 185         uint32_t i;
 186
 187         for( i = 0; i < ( ( col )? m: n ); i++ )
 188             ma[ ( col )? n: i ][ ( col )? i: m ] = v[ i ];
 189
 190     }
 191
 192 // -----------------------------------------------------------------------------
 193     /** Performs a vectorial addition.
 194      *
 195      *  @param v Result.
 196      *  @param v_l Left operand.
 197      *  @param v_r Right operand.
 198      *  @param n Cardinality.
 199      */
 200     template< class T >
 201     inline
 202     void VectorAdd( T* v, T* v_l, T* v_r, uint32_t n )
 203     {
 204         uint32_t i;
 205
 206         for( i = 0; i < n; v[ i ] = v_l[ i ] + v_r[ i ], i++ );
 207
 208     }
 209
 210 // -----------------------------------------------------------------------------
 211     /** Performs a vectorial substraction.
 212      *
 213      *  @param v Result.
 214      *  @param v_l Left operand.
 215      *  @param v_r Right operand.
 216      *  @param n Cardinality.
 217      */
 218     template< class T >
 219     inline
 220     void VectorSubtract( T* v, T* v_l, T* v_r, uint32_t n )
 221     {
 222         uint32_t i;
 223
 224         for( i = 0; i < n; v[ i ] = v_l[ i ] - v_r[ i ], i++ );
 225
 226     }
 227
 228 // -----------------------------------------------------------------------------
 229     /** Performs a vectorial cross product.
 230      *
 231      *  @param v Result.
 232      *  @param v_l Left operand.
 233      *  @param v_r Right operand.
 234      */
 235     template< class T >
 236     inline
 237     void VectorCrossProduct( T* v, T* v_l, T* v_r )
 238     {
 239         v[ 0 ] = ( v_l[ 1 ] * v_r[ 2 ] ) - ( v_l[ 2 ] * v_r[ 1 ] );
 240         v[ 1 ] = ( v_l[ 2 ] * v_r[ 0 ] ) - ( v_l[ 0 ] * v_r[ 2 ] );
 241         v[ 2 ] = ( v_l[ 0 ] * v_r[ 1 ] ) - ( v_l[ 1 ] * v_r[ 0 ] );
 242
 243     }
 244
 245 // -----------------------------------------------------------------------------
 246     /** Performs a vectorial dot product.
 247      *
 248      *  @param v_l Left operand.
 249      *  @param v_r Right operand.
 250      *  @param n Cardinality.
 251      *  @return Dot product.
 252      */
 253     template< class T >
 254     inline
 255     T VectorDotProduct( T* v_l, T* v_r, uint32_t n )
 256     {
 257         T ret;
 258         uint32_t i;
 259
 260         for( i = 0, ret = ( T )0; i < n; ret = ret + ( v_l[ i ] * v_r[ i ] ), i++ );
 261         return( ret );
 262
 263     }
 264
 265 // -----------------------------------------------------------------------------
 266     /** Performs a vector-scalar product.
 267      *
 268      *  @param v Result.
 269      *  @param v_l Left operand.
 270      *  @param s Scalar.
 271      *  @param n Cardinality.
 272      */
 273     template< class T >
 274     inline
 275     void VectorScalarProduct( T* v, T* v_l, T s, uint32_t n )
 276     {
 277         uint32_t i;
 278
 279         for( i = 0; i < n; v[ i ] = v_l[ i ] * s, i++ );
 280
 281     }
 282
 283 // -----------------------------------------------------------------------------
 284     /** Calculates vectorial L2 norm.
 285      *
 286      *  @param v Vector.
 287      *  @param n Cardinality.
 288      */
 289     template< class T >
 290     inline
 291     T VectorNorm( T* v, uint32_t n )
 292     {
 293         return( ( T )( sqrt( ( double )( VectorDotProduct< T >( v, v, n ) ) ) ) );
 294
 295     }
 296
 297 // -----------------------------------------------------------------------------
 298     /** Calculates a normal vector, usign the L2 norm.
 299      *
 300      *  @param v Result.
 301      *  @param v_o Original vector.
 302      *  @param n Cardinality.
 303      */
 304     template< class T >
 305     inline
 306     void VectorNormalize( T* v, T* v_o, uint32_t n )
 307     {
 308         uint32_t i;
 309         T norm = VectorNorm< T >( v_o, n );
 310
 311         norm = ( norm == ( T )0 )? ( T )1: norm;
 312         for( i = 0; i < n; v[ i ] = v_o[ i ] / norm, i++ );
 313
 314     }
 315
 316 // -----------------------------------------------------------------------------
 317 //@}
 318 /// Functions for matrices.
 319 //@{
 320 // -----------------------------------------------------------------------------
 321     /** Stream out function for matrices.
 322      *
 323      *  @param o Output stream.
 324      *  @param ma Matrix.
 325      *  @param n Columns.
 326      *  @param m Rows.
 327      */
 328     template< class T >
 329     inline
 330     void MatrixPrint( std::ostream& o, T** ma, uint32_t n, uint32_t m )
 331     {
 332         uint32_t i, j;
 333
 334         o << n << " x " << m << endl;
 335         for( j = 0; j < m; j++ ) {
 336
 337             for( i = 0; i < n; o << ma[ i ][ j ] << " ", i++ );
 338             o << endl;
 339
 340         } // rof
 341
 342     }
 343
 344 // -----------------------------------------------------------------------------
 345     /** Allocates memory for a matrix.
 346      *
 347      *  @param n Columns.
 348      *  @param m Rows.
 349      */
 350     template< class T >
 351     inline
 352     T** MatrixAllocateMemory( uint32_t n, uint32_t m )
 353     {
 354         uint32_t i;
 355         T** ma = new T*[ n ];
 356
 357         for( i = 0; i < n; i++ ) {
 358
 359             ma[ i ] = new T[ m ];
 360             memset( ma[ i ], 0, sizeof( T ) * m );
 361
 362         } // rof
 363         return( ma );
 364
 365     }
 366
 367 // -----------------------------------------------------------------------------
 368     /** Copies.
 369      *
 370      *  @param ma Target matrix.
 371      *  @param ma_o Source matrix.
 372      *  @param n Columns.
 373      *  @param m Rows.
 374      */
 375     template< class T >
 376     inline
 377     void MatrixAssignMemory( T** ma, T** ma_o, uint32_t n, uint32_t m )
 378     {
 379         uint32_t i;
 380
 381         for( i = 0; i < n; i++ )
 382             memcpy( ma[ i ], ma_o[ i ], sizeof( T ) * m );
 383
 384     }
 385
 386 // -----------------------------------------------------------------------------
 387     /** Allocates & copies.
 388      *
 389      *  @param ma Matrix.
 390      *  @param n Columns.
 391      *  @param m Rows.
 392      */
 393     template< class T >
 394     inline
 395     T** MatrixCopyMemory( T** ma, uint32_t n, uint32_t m )
 396     {
 397         T** n_ma = MatrixAllocateMemory< T >( n, m );
 398         MatrixAssignMemory< T >( n_ma, ma, n, m );
 399         return( n_ma );
 400
 401     }
 402
 403 // -----------------------------------------------------------------------------
 404     /** Scalar assignation.
 405      *
 406      *  @param ma Matrix.
 407      *  @param s Scalar.
 408      *  @param n Columns.
 409      *  @param m Rows.
 410      */
 411     template< class T >
 412     inline
 413     void MatrixAssignScalar( T** ma, T s, uint32_t n, uint32_t m )
 414     {
 415         uint32_t i, j;
 416
 417         for( i = 0; i < n; i++ )
 418             for( j = 0; j < m; j++ )
 419                 ma[ i ][ j ] = s;
 420
 421     }
 422
 423 // -----------------------------------------------------------------------------
 424     /** Converts a matrix in a vector.
 425      *
 426      *  @param v Vector.
 427      *  @param ma Matrix.
 428      *  @param n Columns.
 429      *  @param m Rows.
 430      *  @param col Convert to a column vector?
 431      */
 432     template< class T >
 433     inline
 434     void MatrixVectorCast( T* v, T** ma, uint32_t n, uint32_t m, bool col )
 435     {
 436         uint32_t i;
 437
 438         for( i = 0; i < ( ( col )? m: n ); i++ )
 439             v[ i ] = ma[ ( col )? n: i ][ ( col )? i: m ];
 440
 441     }
 442
 443 // -----------------------------------------------------------------------------
 444     /** Frees matrix memory.
 445      *
 446      *  @param ma Matrix.
 447      *  @param n Columns.
 448      */
 449     template< class T >
 450     inline
 451     void MatrixFreeMemory( T** ma, uint32_t n )
 452     {
 453         uint32_t i;
 454
 455         if( ma ) {
 456
 457             for( i = 0; i < n; i++ )
 458                 if( ma[ i ] ) delete [] ma[ i ];
 459             delete [] ma;
 460
 461         } // fi
 462
 463     }
 464
 465 // -----------------------------------------------------------------------------
 466     /** Performs a matricial addition.
 467      *
 468      *  @param ma Result.
 469      *  @param ma_l Left operand.
 470      *  @param ma_r Right operand.
 471      *  @param n Columns.
 472      *  @param m Rows.
 473      */
 474     template< class T >
 475     inline
 476     void MatrixAdd( T** ma, T** ma_l, T** ma_r, uint32_t n, uint32_t m )
 477     {
 478         uint32_t i, j;
 479
 480         for( i = 0; i < n; i++ )
 481             for( j = 0; j < m; ma[ i ][ j ] = ma_l[ i ][ j ] + ma_r[ i ][ j ], j++ );
 482
 483     }
 484
 485 // -----------------------------------------------------------------------------
 486     /** Performs a matricial substraction.
 487      *
 488      *  @param ma Result.
 489      *  @param ma_l Left operand.
 490      *  @param ma_r Right operand.
 491      *  @param n Columns.
 492      *  @param m Rows.
 493      */
 494     template< class T >
 495     inline
 496     void MatrixSubtract( T** ma, T** ma_l, T** ma_r, uint32_t n, uint32_t m )
 497     {
 498         uint32_t i, j;
 499
 500         for( i = 0; i < n; i++ )
 501             for( j = 0; j < m; ma[ i ][ j ] = ma_l[ i ][ j ] - ma_r[ i ][ j ], j++ );
 502
 503     }
 504
 505 // -----------------------------------------------------------------------------
 506     /** Performs a matricial product.
 507      *
 508      *  @param ma Result.
 509      *  @param ma_l Left operand.
 510      *  @param ma_r Right operand.
 511      *  @param n Left columns.
 512      *  @param m Left rows/right columns.
 513      *  @param nr Right columns.
 514      */
 515     template< class T >
 516     inline
 517     void MatrixProduct( T** ma, T** ma_l, T** ma_r, uint32_t n, uint32_t m, uint32_t nr )
 518     {
 519         unsigned i, j, k;
 520
 521         for( i = 0; i < nr; i++ )
 522             for( j = 0; j < m; j++ )
 523                 for(
 524                     k = 0, ma[ i ][ j ] = ( T )0;
 525                     k < n;
 526                     ma[ i ][ j ] = ma[ i ][ j ] + ( ma_l[ k ][ j ] * ma_r[ i ][ k ] ), k++
 527                     );
 528
 529     }
 530
 531 // -----------------------------------------------------------------------------
 532     /** Performs a scalar-matrix product.
 533      *
 534      *  @param ma Result.
 535      *  @param ma_l Left operand.
 536      *  @param s Scalar.
 537      *  @param n Columns.
 538      *  @param m Rows.
 539      */
 540     template< class T >
 541     inline
 542     void MatrixScalarProduct( T** ma, T** ma_l, T s, uint32_t n, uint32_t m )
 543     {
 544         uint32_t i, j;
 545
 546         for( i = 0; i < n; i++ )
 547             for( j = 0; j < m; ma[ i ][ j ] = ma_l[ i ][ j ] * s, j++ );
 548
 549     }
 550
 551 // -----------------------------------------------------------------------------
 552     /** Gets the matrix cofactor.
 553      *
 554      *  @param ma Result.
 555      *  @param ma_o Matrix.
 556      *  @param i Column index.
 557      *  @param j Row index.
 558      *  @param n Columns.
 559      *  @param m Rows.
 560      */
 561     template< class T >
 562     inline
 563     void MatrixCofactor( T** ma, T** ma_o, uint32_t i, uint32_t j, uint32_t n, uint32_t m )
 564     {
 565         uint32_t k, l;
 566
 567         for( k = 0; k < i; k++ ) {
 568
 569             for( l = 0;     l < j; l++ ) ma[ k ][ l ]     = ma_o[ k ][ l ];
 570             for( l = j + 1; l < m; l++ ) ma[ k ][ l - 1 ] = ma_o[ k ][ l ];
 571
 572         } // rof
 573
 574         for( k = i + 1; k < n; k++ ) {
 575
 576             for( l = 0;     l < j; l++ ) ma[ k - 1 ][ l ]     = ma_o[ k ][ l ];
 577             for( l = j + 1; l < m; l++ ) ma[ k - 1 ][ l - 1 ] = ma_o[ k ][ l ];
 578
 579         } // rof
 580
 581     }
 582
 583 // -----------------------------------------------------------------------------
 584     /** Gets the matrix determinant (square matrices only).
 585      *
 586      *  @param ma Matrix.
 587      *  @param n Dimension.
 588      */
 589     template< class T >
 590     inline
 591     T MatrixDeterminant( T** ma, uint32_t n )
 592     {
 593         uint32_t k;
 594         T** tmp = NULL;
 595         T ret = ( T )0, c;
 596
 597         //  All these cases are for speed-up calculations.
 598         if( n == 1 ) {
 599
 600             ret = ma[ 0 ][ 0 ];
 601
 602         } else if( n == 2 ) {
 603
 604             ret =
 605                 ( ma[ 0 ][ 0 ] * ma[ 1 ][ 1 ] ) -
 606                 ( ma[ 1 ][ 0 ] * ma[ 0 ][ 1 ] );
 607
 608         } else if( n == 3 ) {
 609
 610             ret =
 611                 ( ma[ 0 ][ 0 ] * ma[ 1 ][ 1 ] * ma[ 2 ][ 2 ] ) -
 612                 ( ma[ 0 ][ 0 ] * ma[ 2 ][ 1 ] * ma[ 1 ][ 2 ] ) -
 613                 ( ma[ 1 ][ 0 ] * ma[ 0 ][ 1 ] * ma[ 2 ][ 2 ] ) +
 614                 ( ma[ 1 ][ 0 ] * ma[ 2 ][ 1 ] * ma[ 0 ][ 2 ] ) +
 615                 ( ma[ 2 ][ 0 ] * ma[ 0 ][ 1 ] * ma[ 1 ][ 2 ] ) -
 616                 ( ma[ 2 ][ 0 ] * ma[ 1 ][ 1 ] * ma[ 0 ][ 2 ] );
 617
 618         } else {
 619
 620             tmp = MatrixAllocateMemory< T >( n - 1, n - 1 );
 621             for( k = 0, c = ( T )1, ret = ( T )0; k < n; k++ ) {
 622
 623                 MatrixCofactor< T >( tmp, ma, k, 0, n, n );
 624                 ret = ret + ( c * ( ma[ k ][ 0 ] * MatrixDeterminant< T >( tmp, n - 1 ) ) );
 625                 c = c * ( T )( 0 - 1 );
 626
 627             } // rof
 628             MatrixFreeMemory< T >( tmp, n - 1 );
 629
 630         } // fi
 631         return( ret );
 632
 633     }
 634
 635 // -----------------------------------------------------------------------------
 636     /** Calculates the inverse using the adjoint method (only square matrices).
 637      *
 638      *  @param ma Result.
 639      *  @param ma_o Matrix.
 640      *  @param n Dimension.
 641      */
 642     template< class T >
 643     inline
 644     void MatrixInverseAdjoint( T** ma, T** ma_o, uint32_t n )
 645     {
 646         uint32_t i, j;
 647         T** tmp;
 648         T c;
 649
 650         tmp = MatrixAllocateMemory< T >( n - 1, n - 1 );
 651         for( i = 0; i < n; i++ ) {
 652
 653             for( j = 0; j < n; j++ ) {
 654
 655                 c = ( ( i + j ) % 2 == 0 )? ( T )1: ( T )( 0 - 1 );
 656                 MatrixCofactor< T >( tmp, ma_o, i, j, n, n );
 657                 ma[ j ][ i ] = c * MatrixDeterminant< T >( tmp, n - 1 );
 658
 659             } // rof
 660
 661         } // rof
 662         MatrixFreeMemory< T >( tmp, n - 1 );
 663
 664         c = ( T )1 / MatrixDeterminant< T >( ma_o, n );
 665         MatrixScalarProduct< T >( ma, ma, c, n, n );
 666
 667     }
 668
 669 // -----------------------------------------------------------------------------
 670     /** Calculates the inverse using the gauss method (only square matrices).
 671      *
 672      *  @param ma Result.
 673      *  @param ma_o Matrix.
 674      *  @param n Dimension.
 675      *  @warning Adapted from a java-implementation: http://www.nauticom.net/www/jdtaft/JavaMatrix.htm
 676      */
 677     template< class T >
 678     inline
 679     void MatrixInverseGauss( T** ma, T** ma_o, uint32_t n )
 680     {
 681         uint32_t i, j, k, n2 = 2 * n;
 682         T** ma_a = MatrixAllocateMemory< T >( n2, n );
 683         T a, b;
 684
 685         //  Augmented matrix initialization
 686         for( i = 0; i < n; i++ )
 687             for( j = 0; j < n; ma_a[ i ][ j ] = ma_o[ i ][ j ], j++ );
 688         for( i = n; i < n2; i++ )
 689             for( j = 0; j < n; ma_a[ i ][ j ] = ( i - n == j )? ( T )1: ( T )0, j++ );
 690
 691         //  Reduction
 692         for( j = 0; j < n; j++ ) {
 693
 694             a = ma_a[ j ][ j ];
 695             if( a != 0 ) for( i = 0; i < n2; ma_a[ i ][ j ] = ma_a[ i ][ j ] / a, i++ );
 696             for( k = 0; k < n; k++ ) {
 697
 698                 if( ( k - j ) != 0 ) {
 699
 700                     b = ma_a[ j ][ k ];
 701                     for(
 702                         i = 0;
 703                         i < n2;
 704                         ma_a[ i ][ k ] = ma_a[ i ][ k ] - ( b * ma_a[ i ][ j ] ), i++
 705                         );
 706
 707                 } // fi
 708
 709             } // rof
 710
 711         } // rof
 712
 713         //  Result assignation
 714         MatrixAssignMemory< T >( ma, ma_a, n, n );
 715         MatrixFreeMemory< T >( ma_a, n2 );
 716
 717     }
 718
 719 // -----------------------------------------------------------------------------
 720     /** Gets the transpose matrix.
 721      *
 722      *  @param ma Matrix result.
 723      *  @param ma_o Matrix.
 724      *  @param n Columns.
 725      *  @param m Rows.
 726      */
 727     template< class T >
 728     inline
 729     void MatrixTranspose( T** ma, T** ma_o, uint32_t n, uint32_t m )
 730     {
 731         uint32_t i, j;
 732
 733         for( i = 0; i < m; i++ )
 734             for( j = 0; j < n; ma[ i ][ j ] = ma_o[ j ][ i ], j++ );
 735
 736     }
 737
 738 // -----------------------------------------------------------------------------
 739     /** Load a square matrix with the identity.
 740      *
 741      *  @param ma Matrix.
 742      *  @param n Dimension.
 743      */
 744     template< class T >
 745     inline
 746     void MatrixLoadIdentity( T** ma, uint32_t n )
 747     {
 748         uint32_t i, j;
 749
 750         for( i = 0; i < n; i++ )
 751             for( j = 0; j < n; ma[ i ][ j ] = ( i == j )? ( T )1: ( T )0, j++ );
 752
 753     }
 754
 755 }
 756 //@}
 757
 758 #endif // GTMLIB__MATH__VMFUNCS__HXX
 759
 760 // EOF - vmfuncs.h