DGtal  1.0.0
SimpleMatrix.ih
1 /**
2  * This program is free software: you can redistribute it and/or modify
3  * it under the terms of the GNU Lesser General Public License as
4  * published by the Free Software Foundation, either version 3 of the
5  * License, or (at your option) any later version.
6  *
7  * This program is distributed in the hope that it will be useful,
8  * but WITHOUT ANY WARRANTY; without even the implied warranty of
9  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
10  * GNU General Public License for more details.
11  *
12  * You should have received a copy of the GNU General Public License
13  * along with this program. If not, see <http://www.gnu.org/licenses/>.
14  *
15  **/
16 
17 /**
18  * @file SimpleMatrix.ih
19  * @author David Coeurjolly (\c david.coeurjolly@liris.cnrs.fr )
20  * Laboratoire d'InfoRmatique en Image et Systèmes d'information - LIRIS (CNRS, UMR 5205), CNRS, France
21  *
22  * @date 2012/07/10
23  *
24  * Implementation of inline methods defined in SimpleMatrix.h
25  *
26  * This file is part of the DGtal library.
27  */
28 
29 
30 //////////////////////////////////////////////////////////////////////////////
31 #include <cstdlib>
32 //////////////////////////////////////////////////////////////////////////////
33 
34 ///////////////////////////////////////////////////////////////////////////////
35 // IMPLEMENTATION of inline methods.
36 ///////////////////////////////////////////////////////////////////////////////
37 
38 ///////////////////////////////////////////////////////////////////////////////
39 // ----------------------- Standard services ------------------------------
40 
41 
42 
43 /**
44  * Destructor.
45  */
46 template <typename T, DGtal::Dimension TM, DGtal::Dimension TN>
47 inline
48 DGtal::SimpleMatrix<T,TM, TN>::~SimpleMatrix()
49 {
50 }
51 
52 
53 /**
54  * Constructor.
55  */
56 template <typename T, DGtal::Dimension TM, DGtal::Dimension TN>
57 inline
58 DGtal::SimpleMatrix<T,TM, TN>::SimpleMatrix()
59 {
60  for ( DGtal::Dimension i = 0; i < TM*TN; ++i )
61  myValues[ i ] = NumberTraits<Component>::ZERO;
62  const Component one = NumberTraits<Component>::ONE;
63  const Component minus_one = -NumberTraits<Component>::ONE;
64  //Cofactor coefs computation
65  for (DGtal::Dimension i=0; i<TM; i++)
66  for (DGtal::Dimension j=0; j<TN; j++)
67  myCofactorCoefs[i*N+j] = ( (i+j) & 1 )
68  ? minus_one : one;
69 }
70 
71 /**
72  * Constructor.
73  */
74 template <typename T, DGtal::Dimension TM, DGtal::Dimension TN>
75 inline
76 DGtal::SimpleMatrix<T,TM, TN>::SimpleMatrix(std::initializer_list<T> values)
77 {
78  DGtal::Dimension i = 0;
79  for ( const T *p = values.begin(); p != values.end() && i < TM*TN; ++p, ++i )
80  myValues[ i ] = *p;
81  for ( ; i < TM*TN; ++i )
82  myValues[ i ] = NumberTraits<Component>::ZERO;
83 
84  const Component one = NumberTraits<Component>::ONE;
85  const Component minus_one = -NumberTraits<Component>::ONE;
86  // Cofactor coefs computation
87  for ( i = 0; i < TM; i++ )
88  for ( DGtal::Dimension j = 0; j < TN; j++ )
89  myCofactorCoefs[ i * N + j ] = ( ( i + j ) & 1 ) ? minus_one : one;
90 }
91 //------------------------------------------------------------------------------
92 
93 template<typename T, DGtal::Dimension TM, DGtal::Dimension TN>
94 inline
95 DGtal::SimpleMatrix<T ,TM, TN>::SimpleMatrix(const Self& other)
96 {
97  for ( DGtal::Dimension i = 0; i < M*N; ++i )
98  {
99  myValues[ i ] = other.myValues[i];
100  myCofactorCoefs[ i ] = other.myCofactorCoefs[i];
101  }
102 }
103 
104 //---------------------------------------------------------------------------
105 template<typename T, DGtal::Dimension TM, DGtal::Dimension TN>
106 inline
107 void
108 DGtal::SimpleMatrix<T, TM, TN>::clear()
109 {
110  for ( DGtal::Dimension i = 0; i < M*N; ++i )
111  myValues[ i ] = NumberTraits<T>::ZERO;
112 }
113 //---------------------------------------------------------------------------
114 template<typename T, DGtal::Dimension TM, DGtal::Dimension TN>
115 inline
116 void
117 DGtal::SimpleMatrix<T, TM, TN>::constant(const T &aSc)
118 {
119  for ( DGtal::Dimension i = 0; i < M*N; ++i )
120  myValues[ i ] = aSc;
121 }
122 //---------------------------------------------------------------------------
123 template<typename T, DGtal::Dimension M, DGtal::Dimension N>
124 inline
125 void
126 DGtal::SimpleMatrix<T, M,N>::identity( )
127 {
128  BOOST_STATIC_ASSERT( M == N);
129 
130  //fast clear
131  this->clear();
132  for (DGtal::Dimension i=0; i<M; i++)
133  this->setComponent( i, i , NumberTraits<T>::ONE);
134 }
135 //---------------------------------------------------------------------------
136 template<typename T, DGtal::Dimension TM, DGtal::Dimension TN>
137 inline
138 typename DGtal::SimpleMatrix<T, TM, TN>::RowVector
139 DGtal::SimpleMatrix<T, TM, TN>::row(const DGtal::Dimension i) const
140 {
141  ASSERT(i<M);
142  RowVector v;
143  for ( DGtal::Dimension j = 0; j < N; ++j )
144  v[ j ] = this->operator()(i,j);
145  return v;
146 }
147 template<typename T, DGtal::Dimension TM, DGtal::Dimension TN>
148 inline
149 typename DGtal::SimpleMatrix<T, TM, TN>::ColumnVector
150 DGtal::SimpleMatrix<T, TM, TN>::column(const DGtal::Dimension j) const
151 {
152  ASSERT(j<N);
153  ColumnVector v;
154  for ( DGtal::Dimension i = 0; i < M; ++i )
155  v[ i ] = this->operator()(i,j);
156  return v;
157 }
158 //---------------------------------------------------------------------------
159 template<typename T, DGtal::Dimension TM, DGtal::Dimension TN>
160 template<typename TC>
161 inline
162 DGtal::SimpleMatrix<T, TM, TN> &
163 DGtal::SimpleMatrix<T, TM, TN>::operator=(const SimpleMatrix<TC,M,N>& other)
164 {
165  for ( DGtal::Dimension i = 0; i < M*N; ++i )
166  myValues[ i ] = static_cast<T>(other.myValues[i]);
167  return *this;
168 }
169 
170 //------------------------------------------------------------------------------
171 template<typename T, DGtal::Dimension TM, DGtal::Dimension TN>
172 inline
173 DGtal::SimpleMatrix<T, TM, TN>
174 DGtal::SimpleMatrix<T, TM, TN>::operator+(const Self& other) const
175 {
176  SimpleMatrix<T,TM,TN> res;
177  for ( DGtal::Dimension i = 0; i < M*N; ++i )
178  res.myValues[ i ] = this->myValues[i] + other.myValues[i];
179  return res;
180 }
181 //------------------------------------------------------------------------------
182 template<typename T, DGtal::Dimension TM, DGtal::Dimension TN>
183 inline
184 DGtal::SimpleMatrix<T, TM, TN> &
185 DGtal::SimpleMatrix<T, TM, TN>::operator+=(const Self& other)
186 {
187  for ( DGtal::Dimension i = 0; i < M*N; ++i )
188  myValues[ i ] += other.myValues[i];
189  return *this;
190 }
191 //------------------------------------------------------------------------------
192 template<typename T, DGtal::Dimension M, DGtal::Dimension N>
193 inline
194 T
195 DGtal::SimpleMatrix<T, M,N>::cofactor(const DGtal::Dimension i,
196  const DGtal::Dimension j ) const
197 {
198  BOOST_STATIC_ASSERT(M == N);
199  return minorDeterminant(i,j)*myCofactorCoefs[i*N+j];
200 }
201 //------------------------------------------------------------------------------
202 template<typename T, DGtal::Dimension M, DGtal::Dimension N>
203 inline
204 DGtal::SimpleMatrix<T, M,N>
205 DGtal::SimpleMatrix<T, M,N>::cofactor( ) const
206 {
207  DGtal::SimpleMatrix<T, M,N> mat;
208  BOOST_STATIC_ASSERT(M == N);
209 
210  for (DGtal::Dimension i=0; i<M; i++)
211  for (DGtal::Dimension j=0; j<M; j++)
212  mat.setComponent( i, j , cofactor(i,j));
213 
214  return mat;
215 }
216 //------------------------------------------------------------------------------
217 template<typename T, DGtal::Dimension M, DGtal::Dimension N>
218 inline
219 T
220 DGtal::SimpleMatrix<T, M,N>::minorDeterminant(const DGtal::Dimension i,
221  const DGtal::Dimension j) const
222 {
223  return SimpleMatrixSpecializations<Self,M,N>::minorDeterminant(*this,i,j);
224 }
225 //------------------------------------------------------------------------------
226 template<typename T, DGtal::Dimension M, DGtal::Dimension N>
227 inline
228 T
229 DGtal::SimpleMatrix<T, M,N>::determinant() const
230 {
231  return SimpleMatrixSpecializations<Self,M,N>::determinant(*this);
232 }
233 
234  //------------------------------------------------------------------------------
235 template<typename T, DGtal::Dimension TM, DGtal::Dimension TN>
236 inline
237 typename DGtal::SimpleMatrix<T, TM, TN>
238 DGtal::SimpleMatrix<T, TM, TN>::inverse() const
239 {
240  BOOST_STATIC_ASSERT(TM == TN);
241 
242  SimpleMatrix<T,TM,TM> r = cofactor().transpose();
243 
244  //determinant
245  T det = determinant();
246  ASSERT(det != 0);
247  return r/det;
248 }
249 
250 //------------------------------------------------------------------------------
251 template<typename T, DGtal::Dimension TM, DGtal::Dimension TN>
252 inline
253 DGtal::SimpleMatrix<T, TM, TN>
254 DGtal::SimpleMatrix<T, TM, TN>::operator-(const Self& other) const
255 {
256  SimpleMatrix<T,TM,TN> res;
257  for ( DGtal::Dimension i = 0; i < M*N; ++i )
258  res.myValues[ i ] = this->myValues[i] - other.myValues[i];
259  return res;
260 }
261 //------------------------------------------------------------------------------
262 template<typename T, DGtal::Dimension TM, DGtal::Dimension TN>
263 inline
264 DGtal::SimpleMatrix<T, TM, TN> &
265 DGtal::SimpleMatrix<T, TM, TN>::operator-=(const Self& other)
266 {
267  for ( DGtal::Dimension i = 0; i < M*N; ++i )
268  myValues[ i ] -= other.myValues[i];
269  return *this;
270 }
271 //------------------------------------------------------------------------------
272 template<typename T, DGtal::Dimension TM, DGtal::Dimension TN>
273 inline
274 bool
275 DGtal::SimpleMatrix<T, TM, TN>::operator==(const Self& other) const
276 {
277  return myValues == other.myValues;
278 }
279 
280 //------------------------------------------------------------------------------
281 template<typename T, DGtal::Dimension TM, DGtal::Dimension TN>
282 inline
283 typename DGtal::SimpleMatrix<T, TN, TM>
284 DGtal::SimpleMatrix<T, TM, TN>::transpose() const
285 {
286  DGtal::SimpleMatrix<T, TN, TM> res;
287  for (DGtal::Dimension i=0; i<M; i++)
288  for (DGtal::Dimension j=0; j<N; j++)
289  res.setComponent(j,i, this->operator()(i,j));
290  return res;
291 }
292 //------------------------------------------------------------------------------
293 template<typename T, DGtal::Dimension M, DGtal::Dimension N>
294 inline
295 typename DGtal::SimpleMatrix<T,M,N>::ColumnVector
296 DGtal::SimpleMatrix<T, M, N>::operator*(const RowVector& other) const
297 {
298  ColumnVector res;
299  for (DGtal::Dimension i=0; i<M; i++)
300  for (DGtal::Dimension k=0; k<N; k++)
301  res[i] += this->operator()(i, k )*other[k];
302 
303  return res;
304 }
305 //------------------------------------------------------------------------------
306 template<typename T, DGtal::Dimension M, DGtal::Dimension N>
307 inline
308 typename DGtal::SimpleMatrix<T,M,M>
309 DGtal::SimpleMatrix<T, M, N>::operator*(const DGtal::SimpleMatrix<T,N,M>& other) const
310 {
311  SimpleMatrix<T,M,M> res;
312  T e = NumberTraits<T>::ZERO;
313  for (DGtal::Dimension i=0; i<M; i++)
314  for (DGtal::Dimension j=0; j<M; j++)
315  {
316  for (DGtal::Dimension k=0; k<N; k++)
317  {
318  e += this->operator()(i, k )*other(k ,j );
319  }
320 
321  res.setComponent(i,j,e);
322 
323  e = NumberTraits<T>::ZERO;
324  }
325  return res;
326 }
327 //------------------------------------------------------------------------------
328 template<typename T, DGtal::Dimension M, DGtal::Dimension N>
329 inline
330 DGtal::SimpleMatrix<T, M,N> &
331 DGtal::SimpleMatrix<T, M, N>::operator/=(const T& other)
332 {
333  for (DGtal::Dimension i=0; i<M*N; i++)
334  this->myValues[i] /= other;
335 
336  return *this;
337 }//------------------------------------------------------------------------------
338 template<typename T, DGtal::Dimension M, DGtal::Dimension N>
339 inline
340 DGtal::SimpleMatrix<T, M,N>
341 DGtal::SimpleMatrix<T, M, N>::operator/(const T& other) const
342 {
343  Self resultat;
344  for (DGtal::Dimension i=0; i<M*N; i++)
345  resultat.myValues[i] = myValues[i]/other;
346 
347  return resultat;
348 }
349 //------------------------------------------------------------------------------
350 template<typename T, DGtal::Dimension M, DGtal::Dimension N>
351 inline
352 DGtal::SimpleMatrix<T, M,N> &
353 DGtal::SimpleMatrix<T, M,N>::operator*=(const T& other)
354 {
355  for (DGtal::Dimension i=0; i<M*N; i++)
356  this->myValues[i] *= other;
357 
358  return *this;
359 }
360 //------------------------------------------------------------------------------
361 template<typename T, DGtal::Dimension M, DGtal::Dimension N>
362 inline
363 DGtal::SimpleMatrix<T, M,N>
364 DGtal::SimpleMatrix<T, M,N>::operator*(const T& other) const
365 {
366  Self resultat;
367  for (DGtal::Dimension i=0; i<M*N; i++)
368  resultat.myValues[i] = other*myValues[i];
369 
370  return resultat;
371 }
372 
373 //------------------------------------------------------------------------------
374 template<typename T, DGtal::Dimension TM, DGtal::Dimension TN>
375 inline
376 void
377 DGtal::SimpleMatrix<T, TM, TN>::setComponent(const DGtal::Dimension i,
378  const DGtal::Dimension j,
379  const T &aValue)
380 {
381  ASSERT(i<M);
382  ASSERT(j<N);
383  myValues[i*N + j] = aValue;
384 }
385 //------------------------------------------------------------------------------
386 template<typename T, DGtal::Dimension TM, DGtal::Dimension TN>
387 inline
388 T
389 DGtal::SimpleMatrix<T, TM, TN>::operator()(const DGtal::Dimension i,
390  const DGtal::Dimension j) const
391 {
392  ASSERT(i<M);
393  ASSERT(j<N);
394  return myValues[i*N + j];
395 }
396 //------------------------------------------------------------------------------
397 template<typename T, DGtal::Dimension TM, DGtal::Dimension TN>
398 inline
399 T&
400 DGtal::SimpleMatrix<T, TM, TN>::operator()(const DGtal::Dimension i,
401  const DGtal::Dimension j)
402 {
403  ASSERT(i<M);
404  ASSERT(j<N);
405  return myValues[i*N + j];
406 }
407 
408 ///////////////////////////////////////////////////////////////////////////////
409 // Interface - public :
410 
411 /**
412  * Writes/Displays the object on an output stream.
413  * @param out the output stream where the object is written.
414  */
415 template <typename T, DGtal::Dimension TM, DGtal::Dimension TN >
416 inline
417 void
418 DGtal::SimpleMatrix<T,TM,TN>::selfDisplay ( std::ostream & out ) const
419 {
420  out << "[SimpleMatrix] "<<M<<"x"<<N<< " [";
421  for(DGtal::Dimension i = 0; i < M; ++i)
422  {
423  out<<"[";
424  for(DGtal::Dimension j = 0; j < N; ++j)
425  out<< this->operator()(i,j)<<" ";
426  out<<"]";
427  }
428  out<<"]";
429 }
430 
431 /**
432  * Checks the validity/consistency of the object.
433  * @return 'true' if the object is valid, 'false' otherwise.
434  */
435 template <typename T,DGtal::Dimension M,DGtal::Dimension N>
436 inline
437 bool
438 DGtal::SimpleMatrix<T,M,N>::isValid() const
439 {
440  return true;
441 }
442 
443 
444 
445 
446 
447 ///////////////////////////////////////////////////////////////////////////////
448 // Implementation of inline functions //
449 
450 template <typename T,DGtal::Dimension M,DGtal::Dimension N>
451 inline
452 std::ostream&
453 DGtal::operator<< ( std::ostream & out,
454  const SimpleMatrix<T,M,N> & object )
455 {
456  object.selfDisplay( out );
457  return out;
458 }
459 
460 template <typename TComponent, DGtal::Dimension TM, DGtal::Dimension TN>
461 inline
462 DGtal::SimpleMatrix<TComponent, TM, TN>
463 operator* ( const TComponent& scalar, const DGtal::SimpleMatrix<TComponent, TM, TN>& matrix)
464 {
465  return matrix * scalar;
466 }