KWStyle - itkDiffusionTensor3D.txx

Matrix View

Description

---

1		/*=========================================================================
2
3		Program:   Insight Segmentation & Registration Toolkit
4		Module:    $RCSfile: itkDiffusionTensor3D.txx.html,v $
5		Language:  C++
6		Date:      $Date: 2006/01/17 19:15:34 $
7		Version:   $Revision: 1.4 $
8
9		Copyright (c) Insight Software Consortium. All rights reserved.
10		See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
11
12		This software is distributed WITHOUT ANY WARRANTY; without even
13		the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
14		PURPOSE.  See the above copyright notices for more information.
15
16		=========================================================================*/
17	DEF	#ifndef _itkDiffusionTensor3D_txx
18	DEF	#define _itkDiffusionTensor3D_txx
19
20		#include "itkDiffusionTensor3D.h"
21		#include "itkNumericTraits.h"
22
23		namespace itk
24		{
25
26
27		/*
28	IND	** Default Constructor
29	IND	**/
30		template<class T>
31		DiffusionTensor3D<T>
32		::DiffusionTensor3D()
33		{
34		}
35
36
37		/*
38	IND	** Constructor with initialization
39	IND	**/
40		template<class T>
41		DiffusionTensor3D<T>
42		::DiffusionTensor3D( const Self & r ):SymmetricSecondRankTensor<T,3>(r)
43		{
44		}
45
46
47	EML
48		/*
49	IND	** Constructor with initialization
50	IND	**/
51		template<class T>
52		DiffusionTensor3D<T>
53		::DiffusionTensor3D( const Superclass & r ):SymmetricSecondRankTensor<T,3>(r)
54		{
55		}
56
57
58	EML
59		/*
60	IND	** Constructor with initialization
61	IND	**/
62		template<class T>
63		DiffusionTensor3D<T>
64		::DiffusionTensor3D( const ComponentType & r ):SymmetricSecondRankTensor<T,3>(r)
65		{
66		}
67
68
69	EML
70	EML
71		/*
72	IND	** Constructor with initialization
73	IND	**/
74		template<class T>
75		DiffusionTensor3D<T>
76	LEN	::DiffusionTensor3D( const ComponentArrayType r ):SymmetricSecondRankTensor<T,3>(r)
77		{
78		}
79
80
81		/*
82	IND	** Assignment Operator
83	IND	**/
84		template<class T>
85		DiffusionTensor3D<T>&
86		DiffusionTensor3D<T>
87		::operator= (const Self& r)
88		{
89		Superclass::operator=(r);
90		return *this;
91		}
92
93
94	EML
95		/*
96	IND	** Assignment Operator
97	IND	**/
98		template<class T>
99		DiffusionTensor3D<T>&
100		DiffusionTensor3D<T>
101		::operator= (const ComponentType & r)
102		{
103		Superclass::operator=(r);
104		return *this;
105		}
106
107
108	EML
109		/*
110	IND	** Assignment Operator
111	IND	**/
112		template<class T>
113		DiffusionTensor3D<T>&
114		DiffusionTensor3D<T>
115		::operator= (const ComponentArrayType r)
116		{
117		Superclass::operator=(r);
118		return *this;
119		}
120
121
122	EML
123		/*
124	IND	** Assignment Operator
125	IND	**/
126		template<class T>
127		DiffusionTensor3D<T>&
128		DiffusionTensor3D<T>
129		::operator= (const Superclass & r)
130		{
131		Superclass::operator=(r);
132		return *this;
133		}
134
135
136	EML
137		/*
138	IND	** Get the Trace, specialized version for 3D.
139	IND	**
140	IND	** Note that the indices are related to the fact
141	IND	** that we store only the upper-right triangle of
142	IND	** the matrix. Like
143	IND	**
144	IND	**       | 0  1  2  |
145	IND	**       | X  3  4  |
146	IND	**       | X  X  5  |
147	IND	**
148	IND	** The trace is therefore the sum of the components
149	IND	** M[0], M[3] and M[5].
150	IND	**
151	IND	**/
152		template<class T>
153		typename DiffusionTensor3D<T>::AccumulateValueType
154		DiffusionTensor3D<T>
155		::GetTrace() const
156		{
157		AccumulateValueType trace = (*this)[0];
158		trace += (*this)[3];
159		trace += (*this)[5];
160		return trace;
161		}
162
163
164		/**
165		*  Compute the value of fractional anisotropy
166		*/
167		template<class T>
168		typename DiffusionTensor3D<T>::RealValueType
169		DiffusionTensor3D<T>
170		::GetFractionalAnisotropy() const
171		{
172		// Computed as
173		// FA = sqrt(1.5*sum(sum(N.*N))/sum((sum(D.*D))))
174		// where N = D - ((1/3)*trace(D)*eye(3,3))
175	LEN	// equation (28) in http://lmi.bwh.harvard.edu/papers/pdfs/2002/westinMEDIA02.pdf
176		const RealValueType isp   = this->GetInnerScalarProduct();
177
178		if( isp > 0.0 )
179		{
180		const RealValueType trace = this->GetTrace();
181		const RealValueType anisotropy = 3.0 * isp - trace * trace;
182		const RealValueType fractionalAnisotropy =
183	IND	********static_cast< RealValueType >( sqrt( anisotropy / ( 2.0 * isp ) ) );
184		return fractionalAnisotropy;
185		}
186
187	SEM,IND	***return 0.0 ;
188
189	IND	***/*
190	IND	***// Computed as
191	LEN,IND	***// FA = sqrt(1.5 * ( \sum_i ( lambda_i - lambda_mean )^2 ) / \sum_i ( lambda_i^2 ) )
192	LEN,IND	***// as in http://splweb.bwh.harvard.edu:8000/pages/papers/martha/DTI_Tech354.pdf
193	IND	***// [lambda = eig(A)].
194	IND	***EigenValuesArrayType eigenValues;
195	IND	***ComputeEigenValues( eigenValues );
196	IND	***eigenValues[0] = vnl_math_abs(eigenValues[0]);
197	IND	***eigenValues[1] = vnl_math_abs(eigenValues[1]);
198	IND	***eigenValues[2] = vnl_math_abs(eigenValues[2]);
199	IND	***const RealValueType norm_E = vnl_math_sqr(eigenValues[0])
200	IND	******************************+ vnl_math_sqr(eigenValues[1])
201	IND	******************************+ vnl_math_sqr(eigenValues[2]);
202
203	IND	***if( norm_E > 0.0 )
204	IND	*****{
205	IND	*****const RealValueType anisotropy =
206	IND	**********************vnl_math_sqr(eigenValues[0] - eigenValues[1]) +
207	IND	**********************vnl_math_sqr(eigenValues[1] - eigenValues[2]) +
208	IND	**********************vnl_math_sqr(eigenValues[2] - eigenValues[0]);
209	LEN,IND	*****const RealValueType fractionalAnisotropy = vcl_sqrt( 0.5 * anisotropy/norm_E);
210	IND	*****return fractionalAnisotropy;
211	IND	*****}
212
213	IND	***return 0.0;
214	IND	****/
215
216		}
217
218
219		/**
220		*  Compute the value of relative anisotropy
221		*/
222		template<class T>
223		typename DiffusionTensor3D<T>::RealValueType
224		DiffusionTensor3D<T>
225		::GetRelativeAnisotropy() const
226		{
227		const RealValueType trace = this->GetTrace();
228		const RealValueType isp   = this->GetInnerScalarProduct();
229
230		// Avoid negative trace and traces small enough to look like a division by
231		// zero.
232		if( trace < NumericTraits< RealValueType >::min() )
233		{
234		return NumericTraits<RealValueType>::Zero;
235		}
236
237		const RealValueType anisotropy = 3.0 * isp - trace * trace;
238
239		if( anisotropy  < NumericTraits< RealValueType >::Zero )
240		{
241		return NumericTraits<RealValueType>::Zero;
242		}
243
244		const RealValueType relativeAnisotropySquared =
245	IND	********static_cast< RealValueType >( anisotropy / ( sqrt( 3.0 ) * trace ) );
246
247		const RealValueType relativeAnisotropy =
248	IND	********static_cast< RealValueType >( sqrt( relativeAnisotropySquared ) );
249
250		return relativeAnisotropy;
251		}
252
253
254	EML
255		/**
256		*  Compute the inner scalar product
257		*/
258		template<class T>
259		typename DiffusionTensor3D<T>::RealValueType
260		DiffusionTensor3D<T>
261		::GetInnerScalarProduct() const
262		{
263
264		const RealValueType xx = (*this)[0];
265		const RealValueType xy = (*this)[1];
266		const RealValueType xz = (*this)[2];
267		const RealValueType yy = (*this)[3];
268		const RealValueType yz = (*this)[4];
269		const RealValueType zz = (*this)[5];
270
271		return ( xx*xx + yy*yy + zz*zz + 2.0*(xy*xy + xz*xz + yz*yz) );
272
273		}
274
275
276	EML
277		} // end namespace itk
278
279		#endif
280

---

Generated by KWStyle 1.0b on Tuesday January,17 at 02:14:06PM

© Kitware Inc.