1 /*=========================================================================
2 *
3 * Copyright Insight Software Consortium
4 *
5 * Licensed under the Apache License, Version 2.0 (the "License");
6 * you may not use this file except in compliance with the License.
7 * You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0.txt
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 *
17 *=========================================================================*/
18 #ifndef itkSymmetricSecondRankTensor_hxx
19 #define itkSymmetricSecondRankTensor_hxx
20
21 #include "itkNumericTraitsTensorPixel.h"
22
23 namespace itk
24 {
25 /**
26 * Assignment Operator from a scalar constant
27 */
28 template< typename T, unsigned int NDimension >
29 SymmetricSecondRankTensor< T, NDimension > &
30 SymmetricSecondRankTensor< T, NDimension >
operator =(const ComponentType & r)31 ::operator=(const ComponentType & r)
32 {
33 BaseArray::operator=(&r);
34 return *this;
35 }
36
37 /**
38 * Assigment from a plain array
39 */
40 template< typename T, unsigned int NDimension >
41 SymmetricSecondRankTensor< T, NDimension > &
42 SymmetricSecondRankTensor< T, NDimension >
operator =(const ComponentArrayType r)43 ::operator=(const ComponentArrayType r)
44 {
45 BaseArray::operator=(r);
46 return *this;
47 }
48
49 /**
50 * Returns a temporary copy of a vector
51 */
52 template< typename T, unsigned int NDimension >
53 SymmetricSecondRankTensor< T, NDimension >
54 SymmetricSecondRankTensor< T, NDimension >
operator +(const Self & r) const55 ::operator+(const Self & r) const
56 {
57 Self result;
58
59 for ( unsigned int i = 0; i < InternalDimension; i++ )
60 {
61 result[i] = ( *this )[i] + r[i];
62 }
63 return result;
64 }
65
66 /**
67 * Returns a temporary copy of a vector
68 */
69 template< typename T, unsigned int NDimension >
70 SymmetricSecondRankTensor< T, NDimension >
71 SymmetricSecondRankTensor< T, NDimension >
operator -(const Self & r) const72 ::operator-(const Self & r) const
73 {
74 Self result;
75
76 for ( unsigned int i = 0; i < InternalDimension; i++ )
77 {
78 result[i] = ( *this )[i] - r[i];
79 }
80 return result;
81 }
82
83 /**
84 * Performs addition in place
85 */
86 template< typename T, unsigned int NDimension >
87 const SymmetricSecondRankTensor< T, NDimension > &
88 SymmetricSecondRankTensor< T, NDimension >
operator +=(const Self & r)89 ::operator+=(const Self & r)
90 {
91 for ( unsigned int i = 0; i < InternalDimension; i++ )
92 {
93 ( *this )[i] += r[i];
94 }
95 return *this;
96 }
97
98 /**
99 * Performs subtraction in place
100 */
101 template< typename T, unsigned int NDimension >
102 const SymmetricSecondRankTensor< T, NDimension > &
103 SymmetricSecondRankTensor< T, NDimension >
operator -=(const Self & r)104 ::operator-=(const Self & r)
105 {
106 for ( unsigned int i = 0; i < InternalDimension; i++ )
107 {
108 ( *this )[i] -= r[i];
109 }
110 return *this;
111 }
112
113 /**
114 * Performs multiplication by a scalar, in place
115 */
116 template< typename T, unsigned int NDimension >
117 const SymmetricSecondRankTensor< T, NDimension > &
118 SymmetricSecondRankTensor< T, NDimension >
operator *=(const RealValueType & r)119 ::operator*=(const RealValueType & r)
120 {
121 for ( unsigned int i = 0; i < InternalDimension; i++ )
122 {
123 ( *this )[i] *= r;
124 }
125 return *this;
126 }
127
128 /**
129 * Performs division by a scalar, in place
130 */
131 template< typename T, unsigned int NDimension >
132 const SymmetricSecondRankTensor< T, NDimension > &
133 SymmetricSecondRankTensor< T, NDimension >
operator /=(const RealValueType & r)134 ::operator/=(const RealValueType & r)
135 {
136 for ( unsigned int i = 0; i < InternalDimension; i++ )
137 {
138 ( *this )[i] /= r;
139 }
140 return *this;
141 }
142
143 /**
144 * Performs multiplication with a scalar
145 */
146 template< typename T, unsigned int NDimension >
147 SymmetricSecondRankTensor< T, NDimension >
148 SymmetricSecondRankTensor< T, NDimension >
operator *(const RealValueType & r) const149 ::operator*(const RealValueType & r) const
150 {
151 Self result;
152
153 for ( unsigned int i = 0; i < InternalDimension; i++ )
154 {
155 result[i] = ( *this )[i] * r;
156 }
157 return result;
158 }
159
160 /**
161 * Performs division by a scalar
162 */
163 template< typename T, unsigned int NDimension >
164 SymmetricSecondRankTensor< T, NDimension >
165 SymmetricSecondRankTensor< T, NDimension >
operator /(const RealValueType & r) const166 ::operator/(const RealValueType & r) const
167 {
168 Self result;
169
170 for ( unsigned int i = 0; i < InternalDimension; i++ )
171 {
172 result[i] = ( *this )[i] / r;
173 }
174 return result;
175 }
176
177 /**
178 * Matrix notation access to elements
179 */
180 template< typename T, unsigned int NDimension >
181 const typename SymmetricSecondRankTensor< T, NDimension >::ValueType &
182 SymmetricSecondRankTensor< T, NDimension >
operator ()(unsigned int row,unsigned int col) const183 ::operator()(unsigned int row, unsigned int col) const
184 {
185 unsigned int k;
186
187 if ( row < col )
188 {
189 k = row * Dimension + col - row * ( row + 1 ) / 2;
190 }
191 else
192 {
193 k = col * Dimension + row - col * ( col + 1 ) / 2;
194 }
195
196 if ( k >= InternalDimension )
197 {
198 k = 0;
199 }
200
201 return ( *this )[k];
202 }
203
204 /**
205 * Matrix notation access to elements
206 */
207 template< typename T, unsigned int NDimension >
208 typename SymmetricSecondRankTensor< T, NDimension >::ValueType &
209 SymmetricSecondRankTensor< T, NDimension >
operator ()(unsigned int row,unsigned int col)210 ::operator()(unsigned int row, unsigned int col)
211 {
212 unsigned int k;
213
214 if ( row < col )
215 {
216 k = row * Dimension + col - row * ( row + 1 ) / 2;
217 }
218 else
219 {
220 k = col * Dimension + row - col * ( col + 1 ) / 2;
221 }
222
223 if ( k >= InternalDimension )
224 {
225 k = 0;
226 }
227
228 return ( *this )[k];
229 }
230
231 /**
232 * Set the Tensor to an Identity.
233 * Set ones in the diagonal and zeroes every where else.
234 */
235 template< typename T, unsigned int NDimension >
236 void
237 SymmetricSecondRankTensor< T, NDimension >
SetIdentity()238 ::SetIdentity()
239 {
240 this->Fill(NumericTraits< T >::ZeroValue());
241 for ( unsigned int i = 0; i < Dimension; i++ )
242 {
243 ( *this )( i, i ) = NumericTraits< T >::OneValue();
244 }
245 }
246
247 /**
248 * Get the Trace
249 */
250 template< typename T, unsigned int NDimension >
251 typename SymmetricSecondRankTensor< T, NDimension >::AccumulateValueType
252 SymmetricSecondRankTensor< T, NDimension >
GetTrace() const253 ::GetTrace() const
254 {
255 AccumulateValueType trace = NumericTraits< AccumulateValueType >::ZeroValue();
256 unsigned int k = 0;
257
258 for ( unsigned int i = 0; i < Dimension; i++ )
259 {
260 trace += ( *this )[k];
261 k += ( Dimension - i );
262 }
263 return trace;
264 }
265
266 /**
267 * Compute Eigen Values
268 */
269 template< typename T, unsigned int NDimension >
270 void
271 SymmetricSecondRankTensor< T, NDimension >
ComputeEigenValues(EigenValuesArrayType & eigenValues) const272 ::ComputeEigenValues(EigenValuesArrayType & eigenValues) const
273 {
274 SymmetricEigenAnalysisType symmetricEigenSystem;
275
276 MatrixType tensorMatrix;
277
278 for ( unsigned int row = 0; row < Dimension; row++ )
279 {
280 for ( unsigned int col = 0; col < Dimension; col++ )
281 {
282 tensorMatrix[row][col] = ( *this )( row, col );
283 }
284 }
285
286 symmetricEigenSystem.ComputeEigenValues(tensorMatrix, eigenValues);
287 }
288
289 /**
290 * Compute Eigen analysis, it returns an array with eigen values
291 * and a Matrix with eigen vectors
292 */
293 template< typename T, unsigned int NDimension >
294 void
295 SymmetricSecondRankTensor< T, NDimension >
ComputeEigenAnalysis(EigenValuesArrayType & eigenValues,EigenVectorsMatrixType & eigenVectors) const296 ::ComputeEigenAnalysis(EigenValuesArrayType & eigenValues,
297 EigenVectorsMatrixType & eigenVectors) const
298 {
299 SymmetricEigenAnalysisType symmetricEigenSystem;
300
301 MatrixType tensorMatrix;
302
303 for ( unsigned int row = 0; row < Dimension; row++ )
304 {
305 for ( unsigned int col = 0; col < Dimension; col++ )
306 {
307 tensorMatrix[row][col] = ( *this )( row, col );
308 }
309 }
310
311 symmetricEigenSystem.ComputeEigenValuesAndVectors(
312 tensorMatrix, eigenValues, eigenVectors);
313 }
314
315 /**
316 * Set the Tensor to a Rotated version of the current tensor.
317 * matrix * self * Transpose(matrix)
318 *
319 */
320 template<typename T,unsigned int NDimension>
321 template <typename TMatrixValueType>
322 SymmetricSecondRankTensor<T,NDimension>
323 SymmetricSecondRankTensor<T,NDimension>
Rotate(const Matrix<TMatrixValueType,NDimension,NDimension> & m) const324 ::Rotate( const Matrix<TMatrixValueType, NDimension, NDimension> & m ) const
325 {
326 Self result;
327 using RotationMatrixType = Matrix<double, NDimension, NDimension>;
328 RotationMatrixType SCT; //self * Transpose(m)
329 for(unsigned int r=0; r<NDimension; r++)
330 {
331 for(unsigned int c=0; c<NDimension; c++)
332 {
333 double sum = 0.0;
334 for(unsigned int t=0; t<NDimension; t++)
335 {
336 sum += (*this)(r,t) * m(c,t);
337 }
338 SCT(r,c) = sum;
339 }
340 }
341 //self = m * sct;
342 for(unsigned int r=0; r<NDimension; r++)
343 {
344 for(unsigned int c=0; c<NDimension; c++)
345 {
346 double sum = 0.0;
347 for(unsigned int t=0; t<NDimension; t++)
348 {
349 sum += m(r,t) * SCT(t,c);
350 }
351 (result)(r,c) = static_cast<T>( sum );
352 }
353 }
354 return result;
355 }
356
357 /**
358 * Pre-multiply the Tensor by a Matrix
359 */
360 template< typename T, unsigned int NDimension >
361 typename SymmetricSecondRankTensor< T, NDimension >::MatrixType
362 SymmetricSecondRankTensor< T, NDimension >
PreMultiply(const MatrixType & m) const363 ::PreMultiply(const MatrixType & m) const
364 {
365 MatrixType result;
366
367 using AccumulateType = typename NumericTraits< T >::AccumulateType;
368 for ( unsigned int r = 0; r < NDimension; r++ )
369 {
370 for ( unsigned int c = 0; c < NDimension; c++ )
371 {
372 AccumulateType sum = NumericTraits< AccumulateType >::ZeroValue();
373 for ( unsigned int t = 0; t < NDimension; t++ )
374 {
375 sum += m(r, t) * ( *this )( t, c );
376 }
377 result(r, c) = static_cast< T >( sum );
378 }
379 }
380 return result;
381 }
382
383 /**
384 * Post-multiply the Tensor by a Matrix
385 */
386 template< typename T, unsigned int NDimension >
387 typename SymmetricSecondRankTensor< T, NDimension >::MatrixType
388 SymmetricSecondRankTensor< T, NDimension >
PostMultiply(const MatrixType & m) const389 ::PostMultiply(const MatrixType & m) const
390 {
391 MatrixType result;
392
393 using AccumulateType = typename NumericTraits< T >::AccumulateType;
394 for ( unsigned int r = 0; r < NDimension; r++ )
395 {
396 for ( unsigned int c = 0; c < NDimension; c++ )
397 {
398 AccumulateType sum = NumericTraits< AccumulateType >::ZeroValue();
399 for ( unsigned int t = 0; t < NDimension; t++ )
400 {
401 sum += ( *this )( r, t ) * m(t, c);
402 }
403 result(r, c) = static_cast< T >( sum );
404 }
405 }
406 return result;
407 }
408
409 /**
410 * Print content to an ostream
411 */
412 template< typename T, unsigned int NDimension >
413 std::ostream &
operator <<(std::ostream & os,const SymmetricSecondRankTensor<T,NDimension> & c)414 operator<<(std::ostream & os, const SymmetricSecondRankTensor< T, NDimension > & c)
415 {
416 for ( unsigned int i = 0; i < c.GetNumberOfComponents(); i++ )
417 {
418 os << static_cast< typename NumericTraits< T >::PrintType >( c[i] ) << " ";
419 }
420 return os;
421 }
422
423 /**
424 * Read content from an istream
425 */
426 template< typename T, unsigned int NDimension >
427 std::istream &
operator >>(std::istream & is,SymmetricSecondRankTensor<T,NDimension> & dt)428 operator>>(std::istream & is, SymmetricSecondRankTensor< T, NDimension > & dt)
429 {
430 for ( unsigned int i = 0; i < dt.GetNumberOfComponents(); i++ )
431 {
432 is >> dt[i];
433 }
434 return is;
435 }
436 } // end namespace itk
437
438 #endif
439