CASToR  2.0
Tomographic Reconstruction (PET/SPECT/CT)
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Macros Groups Pages
oMatrix.cc
Go to the documentation of this file.
1 /*
2 This file is part of CASToR.
3 
4  CASToR is free software: you can redistribute it and/or modify it under the
5  terms of the GNU General Public License as published by the Free Software
6  Foundation, either version 3 of the License, or (at your option) any later
7  version.
8 
9  CASToR is distributed in the hope that it will be useful, but WITHOUT ANY
10  WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
11  FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
12  details.
13 
14  You should have received a copy of the GNU General Public License along with
15  CASToR (in file GNU_GPL.TXT). If not, see <http://www.gnu.org/licenses/>.
16 
17 Copyright 2017-2018 all CASToR contributors listed below:
18 
19  --> current contributors: Thibaut MERLIN, Simon STUTE, Didier BENOIT, Claude COMTAT, Marina FILIPOVIC, Mael MILLARDET
20  --> past contributors: Valentin VIELZEUF
21 
22 This is CASToR version 2.0.
23 */
24 
31 #include "oMatrix.hh"
32 
33 // =====================================================================
34 // ---------------------------------------------------------------------
35 // ---------------------------------------------------------------------
36 // =====================================================================
37 /*
38  \brief oMatrix constructor.
39  Initialize the member variables to their default values.
40 */
42 {
43  m_lin = 0;
44  m_col = 0;
45 
46  m2p_Mat = NULL;
47 }
48 
49 
50 
51 // =====================================================================
52 // ---------------------------------------------------------------------
53 // ---------------------------------------------------------------------
54 // =====================================================================
55 /*
56  \param nl : a number of lines
57  \param nc : a number of colons
58  \brief oMatrix constructor.
59  Instanciate a Matrix structure with the number of lines and colons provided in parameters
60 */
61 oMatrix::oMatrix(uint16_t nl, uint16_t nc)
62 {
63  m_lin = nl;
64  m_col = nc;
65  m2p_Mat = new FLTNB *[nl];
66 
67  for(uint16_t l=0 ; l<nl ; l++)
68  m2p_Mat[l] = new FLTNB[nc];
69 }
70 
71 
72 
73 // =====================================================================
74 // ---------------------------------------------------------------------
75 // ---------------------------------------------------------------------
76 // =====================================================================
77 /*
78  \brief oMatrix destructor.
79  Free memory of the oMatrix object.
80 */
82 {
83  for(uint16_t l=0 ; l<m_lin ; l++)
84  if(m2p_Mat[l]) delete[] m2p_Mat[l];
85 
86  if(m2p_Mat) delete[] m2p_Mat;
87 }
88 
89 
90 
91 // =====================================================================
92 // ---------------------------------------------------------------------
93 // ---------------------------------------------------------------------
94 // =====================================================================
95 /*
96  \fn Allocate
97  \param nl : a number of lines
98  \param nc : a number of colons
99  \brief Instanciate a Matrix structure with the number of lines and colons provided in parameters
100 */
101 void oMatrix::Allocate(uint16_t nl, uint16_t nc)
102 {
103  // No verbosity in oMatrix structure. Show this only if CASTOR_DEBUG is enabled (for bug tracking)
104  #ifdef CASTOR_DEBUG
105  Cout("oMatrix::Allocate() ...");
106  #endif
107 
108  // Free memory in case the matrix had already been allocated
109  if(m2p_Mat != NULL)
110  {
111  for(uint16_t l=0 ; l<m_lin ; l++)
112  if(m2p_Mat[l]) delete[] m2p_Mat[l];
113 
114  delete[] m2p_Mat;
115  }
116 
117  m_lin = nl;
118  m_col = nc;
119  m2p_Mat = new FLTNB *[nl];
120 
121  for(uint16_t l=0 ; l<nl ; l++)
122  m2p_Mat[l] = new FLTNB[nc];
123 }
124 
125 
126 
127 // =====================================================================
128 // ---------------------------------------------------------------------
129 // ---------------------------------------------------------------------
130 // =====================================================================
131 /*
132  \fn SetMatriceElt
133  \param l : a line index
134  \param c : a colon index
135  \param a_val : a value to initialize the matrix element with
136  \brief set the matrix element corresponding to the argument indices with the provided value.
137  \return 0 if success, positive value otherwise
138 */
139 int oMatrix::SetMatriceElt(uint16_t l, uint16_t c, FLTNB a_val)
140 {
141  if(l>=m_lin || c>=m_col)
142  {
143  Cerr("***** oMatrix::SetMatriceElt()-> Nb of (lin,col) ("<<l+1<<","<<c+1<<") in parameters ");
144  Cerr("> to the number of (lin,col) of this matrix ("<<m_lin<<","<<m_col<<") !" << endl);
145  return 1;
146  }
147 
148  m2p_Mat[l][c] = a_val;
149 
150  return 0;
151 }
152 
153 
154 
155 // =====================================================================
156 // ---------------------------------------------------------------------
157 // ---------------------------------------------------------------------
158 // =====================================================================
166 FLTNB oMatrix::GetMatriceElt(uint16_t l,uint16_t c)
167 {
168  return m2p_Mat[l][c];
169 }
170 
171 
172 
173 // =====================================================================
174 // ---------------------------------------------------------------------
175 // ---------------------------------------------------------------------
176 // =====================================================================
177 /*
178  \fn Multiplication
179  \param ap_Mtx : input matrix to multiply to the matrix from which the function is called
180  \param ap_MtxResult : output matrix containing the result of the multiplication
181  \brief Multiply the member matrix with the matrix provided in 1st parameter
182  Return the result in the matric provided in 2nd parameter
183  \return 0 if success, positive value otherwise
184 */
185 int oMatrix::Multiplication(oMatrix *ap_Mtx, oMatrix *ap_MtxResult)
186 {
187  HPFLTNB sum=0;
188 
189  if ( m_col != ap_Mtx->m_lin )
190  {
191  Cerr("***** oMatrix::Multiplication()-> Not matching number of colons and lines of the two matrices !");
192  return 1;
193  }
194  else if (ap_MtxResult == NULL)
195  {
196  Cerr("***** oMatrix::Multiplication()-> The resulting matrix has not been allocated !");
197  return 1;
198  }
199  else
200  {
201  for ( uint16_t tl = 0; tl < m_lin; tl++ )
202  for ( uint16_t c = 0; c < ap_Mtx->m_col; c++ )
203  {
204  for ( uint16_t l = 0; l < ap_Mtx->m_lin; l++ )
205  {
206  sum += GetMatriceElt(tl,l) * ap_Mtx->GetMatriceElt(l,c);
207  }
208  ap_MtxResult->SetMatriceElt(tl,c,sum);
209  sum = 0;
210  }
211  }
212 
213  return 0;
214 }
215 
216 
217 
218 
219 // =====================================================================
220 // ---------------------------------------------------------------------
221 // ---------------------------------------------------------------------
222 // =====================================================================
223 /*
224  \fn int oMatrix::Transpose()
225  \param ap_MtxResult : output matrix containing the result of the transposition
226  \brief Transpose the elements of the matrix
227  \return 0 if success, positive value otherwise
228 */
230 {
231  if ( this->m_col != ap_Mtx->m_lin
232  || this->m_lin != ap_Mtx->m_col )
233  {
234  Cerr("***** oMatrix::Transpose()-> Not matching number of colons <-> lines between the input and transpose matrices !");
235  return 1;
236  }
237  if (ap_Mtx == NULL)
238  {
239  Cerr("***** oMatrix::Transpose()-> The resulting matrix has not been allocated !");
240  return 1;
241  }
242 
243  uint16_t lsize = m_lin>m_col ? m_lin : m_col;
244 
245  for ( uint16_t l = 0; l < lsize; l++ )
246  for ( uint16_t c = l; c < lsize; c++ )
247  {
248  HPFLTNB tmp = 0.;
249 
250  if( l < m_lin
251  && c < m_col )
252  tmp = m2p_Mat[l][c];
253 
254  if(l < ap_Mtx->m_lin
255  && c < ap_Mtx->m_col
256  && l < m_col
257  && c < m_lin )
258  ap_Mtx->SetMatriceElt( l , c , m2p_Mat[c][l] );
259 
260  if(l < ap_Mtx->m_col
261  && c < ap_Mtx->m_lin
262  && l < m_lin
263  && c < m_col )
264  ap_Mtx->SetMatriceElt( c , l , tmp );
265  }
266 
267  return 0;
268 }
269 
270 
271 
272 
273 // =====================================================================
274 // ---------------------------------------------------------------------
275 // ---------------------------------------------------------------------
276 // =====================================================================
277 /*
278  \fn int oMatrix::Inverse()
279  \param ap_MtxResult : output matrix containing the inverse matrix
280  \brief Inverse the matrix on which the function is called
281  An error is returned if the matrix is not square,
282  or if there is a nil component on the diagonal
283  The resulting matrix object must be different than the input matrix
284  \return 0 if success, positive value otherwise
285 */
287 {
288  // Check matching nb colons/lines
289  if ( (this->m_col != this->m_lin)
290  || (ap_Mtx->m_col != ap_Mtx->m_lin) )
291  {
292  Cerr("***** oMatrix::Inverse()-> Not matching number of colons <-> lines between the input and transpose matrices !");
293  return 1;
294  }
295  // Check if resulting matrix has been allocated
296  if (ap_Mtx == NULL)
297  {
298  Cerr("***** oMatrix::Inverse()-> The resulting matrix has not been allocated !");
299  return 1;
300  }
301  // Check if resulting matrix is the same as this matrix
302  // Throw error as the resulting matrix will overwrite the components of the input matrix in the process
303  if (this == ap_Mtx)
304  {
305  Cerr("***** oMatrix::Inverse()-> The resulting matrix has not been allocated !");
306  return 1;
307  }
308 
309  for ( uint16_t l = 0; l < m_lin; l++ )
310  for ( uint16_t c = 0; c < m_col; c++ )
311  {
312  if (m2p_Mat[l][c] == 0)
313  {
314  Cerr("***** oMatrix::Inverse()-> One element on the diagonal is nil !");
315  return 1;
316  }
317  }
318 
319  oMatrix* p_tMtx = new oMatrix(m_lin, m_col);
320 
321 
322  // Init
323  for( uint16_t l=0 ; l<m_lin ; l++ )
324  for ( uint16_t c = 0; c < m_col; c++ )
325  {
326  p_tMtx->SetMatriceElt( l, c, 0. );
327  ap_Mtx->SetMatriceElt( l, c, m2p_Mat[l][c] );
328  }
329 
330  FLTNB** pMtx = ap_Mtx->GetMtx();
331 
332  for( uint16_t l=0 ; l<m_lin ; l++ )
333  {
334  for ( uint16_t c = 0; c < m_col; c++ )
335  {
336  p_tMtx->SetMatriceElt(l, l, 1./pMtx[l][l]);
337 
338  if(c != l)
339  p_tMtx->SetMatriceElt(l, c, -pMtx[l][c] / pMtx[l][l]);
340 
341  for( uint16_t r=0 ; r<m_lin ; r++ )
342  {
343  if(r!=l)
344  p_tMtx->SetMatriceElt(r, l, pMtx[r][l]/pMtx[l][l]);
345  if(c!=l && r!=l)
346  p_tMtx->SetMatriceElt(r, c, pMtx[r][c] - pMtx[l][c]*pMtx[r][l]/pMtx[l][l]);
347  }
348 
349  }
350 
351  for(int ll=0 ; ll<m_lin ; ll++)
352  for(int cc=0 ; cc<m_col ; cc++)
353  pMtx[ll][cc] = p_tMtx->GetMatriceElt(ll,cc) ;
354 
355  }
356 
357  delete p_tMtx;
358 
359  return 0;
360 }
361 
362 
363 // =====================================================================
364 // ---------------------------------------------------------------------
365 // ---------------------------------------------------------------------
366 // =====================================================================
367 /*
368  \fn int oMatrix::SetXRotMtx()
369  \param ang : angle in degree
370  \brief Set a (3,3) X-axis rotation matrix using the provided angle
371  \return 0 if sucess, positive value otherwise
372 */
374 {
375  if(m_lin !=3 || m_col != 3)
376  {
377  Cerr("***** oMatrix::SetXRotMtx()-> The matrix must be (3,3) in order to use this function !");
378  return 1;
379  }
380 
381  SetMatriceElt(0,0, 1 );
382  SetMatriceElt(0,1, 0 );
383  SetMatriceElt(0,2, 0 );
384  SetMatriceElt(1,0, 0 );
385  SetMatriceElt(1,1, cos(ang));
386  SetMatriceElt(1,2, -sin(ang));
387  SetMatriceElt(2,0, 0 );
388  SetMatriceElt(2,1, sin(ang));
389  SetMatriceElt(2,2, cos(ang)) ;
390 
391  return 0;
392 }
393 
394 
395 
396 
397 // =====================================================================
398 // ---------------------------------------------------------------------
399 // ---------------------------------------------------------------------
400 // =====================================================================
401 /*
402  \fn int oMatrix::SetYRotMtx()
403  \param ang : angle in degree
404  \brief Set a (3,3) Y-axis rotation matrix using the provided angle
405  \return 0 if sucess, positive value otherwise
406 */
408 {
409  if(m_lin !=3 || m_col != 3)
410  {
411  Cerr("***** oMatrix::SetYRotMtx()-> The matrix must be (3,3) in order to use this function !");
412  return 1;
413  }
414 
415  SetMatriceElt(0,0, cos(ang));
416  SetMatriceElt(0,1, 0 );
417  SetMatriceElt(0,2, sin(ang));
418  SetMatriceElt(1,0, 0 );
419  SetMatriceElt(1,1, 1 );
420  SetMatriceElt(1,2, 0 );
421  SetMatriceElt(2,0, -sin(ang));
422  SetMatriceElt(2,1, 0 );
423  SetMatriceElt(2,2, cos(ang));
424 
425  return 0;
426 }
427 
428 
429 
430 
431 // =====================================================================
432 // ---------------------------------------------------------------------
433 // ---------------------------------------------------------------------
434 // =====================================================================
435 /*
436  \fn int oMatrix::SetZRotMtx()
437  \param ang : angle in degree
438  \brief Set a (3,3) Z-axis rotation matrix using the provided angle
439  \return 0 if sucess, positive value otherwise
440 */
442 {
443  if(m_lin !=3 || m_col != 3)
444  {
445  Cerr("***** oMatrix::SetZRotMtx()-> The matrix must be (3,3) in order to use this function !");
446  return 1;
447  }
448 
449  SetMatriceElt(0,0, cos(ang) );
450  SetMatriceElt(0,1, -sin(ang) );
451  SetMatriceElt(0,2, 0);
452  SetMatriceElt(1,0, sin(ang));
453  SetMatriceElt(1,1, cos(ang) );
454  SetMatriceElt(1,2, 0);
455  SetMatriceElt(2,0, 0);
456  SetMatriceElt(2,1, 0);
457  SetMatriceElt(2,2, 1);
458 
459  return 0;
460 }
461 
462 
463 
464 
465 // =====================================================================
466 // ---------------------------------------------------------------------
467 // ---------------------------------------------------------------------
468 // =====================================================================
469 /*
470  \fn int oMatrix::Describe()
471  \brief Display the element of the matrix
472 */
474 {
475  for ( uint16_t l = 0; l < m_lin; l++ )
476  {
477  for ( uint16_t c = 0; c < m_col; c++ )
478  cout << m2p_Mat[l][c] << " ; ";
479 
480  cout << endl;
481  }
482 }
int Inverse(oMatrix *ap_MtxResult)
Inverse the matrix on which the function is called An error is returned if the matrix is not square...
Definition: oMatrix.cc:286
#define FLTNB
Definition: gVariables.hh:81
#define HPFLTNB
Definition: gVariables.hh:83
FLTNB GetMatriceElt(uint16_t l, uint16_t c)
Definition: oMatrix.cc:166
~oMatrix()
oMatrix destructor. Free memory of the oMatrix object.
Definition: oMatrix.cc:81
FLTNB ** GetMtx()
Definition: oMatrix.hh:96
oMatrix()
oMatrix constructor. Initialize the member variables to their default values.
Definition: oMatrix.cc:41
#define Cerr(MESSAGE)
int SetMatriceElt(uint16_t l, uint16_t c, FLTNB a_val)
Set the matrix element corresponding to the argument indices with the provided value.
Definition: oMatrix.cc:139
uint16_t m_lin
Definition: oMatrix.hh:157
Declaration of class oMatrix.
void Describe()
Display the element of the matrix.
Definition: oMatrix.cc:473
uint16_t m_col
Definition: oMatrix.hh:158
int SetZRotMtx(FLTNB ang)
Set a (3,3) Z-axis rotation matrix using the provided angle.
Definition: oMatrix.cc:441
int SetXRotMtx(FLTNB ang)
Set a (3,3) X-axis rotation matrix using the provided angle.
Definition: oMatrix.cc:373
int SetYRotMtx(FLTNB ang)
Set a (3,3) Y-axis rotation matrix using the provided angle.
Definition: oMatrix.cc:407
Structure designed for basic matrices operations.
Definition: oMatrix.hh:42
#define Cout(MESSAGE)
int Multiplication(oMatrix *ap_Mtx, oMatrix *ap_MtxResult)
Multiply the member matrix with the matrix provided in 1st parameter Return the result in the matric ...
Definition: oMatrix.cc:185
FLTNB ** m2p_Mat
Definition: oMatrix.hh:160
int Transpose(oMatrix *a_MtxResult)
Transpose the elements of the matrix.
Definition: oMatrix.cc:229
void Allocate(uint16_t nl, uint16_t nc)
Instanciate a Matrix structure with the number of lines and colons provided in parameters.
Definition: oMatrix.cc:101