DMRITool  v0.1.1-139-g860d86b4
Diffusion MRI Tool
itkDiffusionTensor.hxx
Go to the documentation of this file.
1 
18 #ifndef __itkDiffusionTensor_hxx
19 #define __itkDiffusionTensor_hxx
20 
21 #include "itkDiffusionTensor.h"
22 #include "utl.h"
23 
24 #include "vnl/vnl_matrix.h"
25 #include <vnl/vnl_math.h>
26 #include "vnl/vnl_vector.h"
27 #include "vnl/vnl_trace.h"
28 #include "vnl/algo/vnl_determinant.h"
29 
30 
31 namespace itk
32 {
33 
37 template< class TPrecision >
38 const typename DiffusionTensor< TPrecision >::ValueType &
39 DiffusionTensor< TPrecision >
40 ::operator()(unsigned int row, unsigned int col) const
41 {
42  unsigned int k;
43 
44  if ( row < col )
45  {
46  k = row * NDimension + col - row * ( row + 1 ) / 2;
47  }
48  else
49  {
50  k = col * NDimension + row - col * ( col + 1 ) / 2;
51  }
52 
53  utlAssert(k < NDegreesOfFreedom, "wrong dimension");
54 
55  return ( *this )[k];
56 }
57 
61 template< class TPrecision >
62 typename DiffusionTensor< TPrecision >::ValueType &
63 DiffusionTensor< TPrecision >
64 ::operator()(unsigned int row, unsigned int col)
65 {
66  unsigned int k;
67 
68  if ( row < col )
69  {
70  k = row * NDimension + col - row * ( row + 1 ) / 2;
71  }
72  else
73  {
74  k = col * NDimension + row - col * ( col + 1 ) / 2;
75  }
76  utlAssert(k < NDegreesOfFreedom, "wrong dimension");
77 
78  return ( *this )[k];
79 }
80 
81 template<class TPrecision >
82 void
83 DiffusionTensor<TPrecision>
84 ::GetEigenValuesVectors(vnl_diag_matrix<TPrecision> & eigenValues, vnl_matrix<TPrecision> & eigenVectors) const
85 {
86  EigenValuesArrayType values;
87  EigenVectorsMatrixType vectors;
88  this->ComputeEigenAnalysis(values, vectors);
89 
90  eigenValues.set_size(NDimension);
91  for ( int i = 0; i < NDimension; i += 1 )
92  eigenValues[i] = values[i];
93  eigenVectors = vectors.GetVnlMatrix().transpose();
94 }
95 
96 template<class TPrecision>
97 template<class TArrayType, class TMatrixType >
98 void
99 DiffusionTensor<TPrecision>
100 ::GetEigenValuesVectorsAnalytic(TArrayType& eigenValues, TMatrixType& eigenVectors) const
101 {
102  if (IsDiagonal(1e-10))
103  {
104  for ( int i = 0; i < 3; ++i )
105  {
106  eigenValues[i] = (*this)(i,i);
107  eigenVectors(i,i) = 1.0;
108  for ( int j = 0; j < 3; ++j )
109  {
110  if (i!=j)
111  eigenVectors(i,j) = 0.0;
112  }
113  }
114 
115  // sort three values in ascending order
116  if (eigenValues[0]>eigenValues[1])
117  {
118  std::swap(eigenValues[0], eigenValues[1]);
119  for ( int i = 0; i < 3; ++i )
120  std::swap(eigenVectors(i,0), eigenVectors(i,1));
121  }
122  if (eigenValues[1]>eigenValues[2])
123  {
124  std::swap(eigenValues[1], eigenValues[2]);
125  for ( int i = 0; i < 3; ++i )
126  std::swap(eigenVectors(i,1), eigenVectors(i,2));
127  }
128  if (eigenValues[0]>eigenValues[1])
129  {
130  std::swap(eigenValues[0], eigenValues[1]);
131  for ( int i = 0; i < 3; ++i )
132  std::swap(eigenVectors(i,0), eigenVectors(i,1));
133  }
134  return;
135  }
136 
137  std::vector<double> detCoef(4);
138  std::vector<std::complex<double> > lambdaVec;
139 
140  const TPrecision* p = this->GetDataPointer();
141  detCoef[0] = p[2]*p[2]*p[3] - 2.0*p[1]*p[2]*p[4] + p[0]*p[4]*p[4] + p[1]*p[1]*p[5] - p[0]*p[3]*p[5];
142  detCoef[1] = -p[1]*p[1] - p[2]*p[2] + p[0]*p[3] - p[4]*p[4] + p[0]*p[5] + p[3]*p[5];
143  detCoef[2] = -p[0] - p[3] - p[5];
144  detCoef[3] = 1;
145  lambdaVec = utl::PolynomialRoot(detCoef);
146 
147  for ( int i = 0; i < 3; ++i )
148  eigenValues[i] = std::real(lambdaVec[i]);
149 
150  // sort three values in ascending order
151  if (eigenValues[0]>eigenValues[1])
152  std::swap(eigenValues[0], eigenValues[1]);
153  if (eigenValues[1]>eigenValues[2])
154  std::swap(eigenValues[1], eigenValues[2]);
155  if (eigenValues[0]>eigenValues[1])
156  std::swap(eigenValues[0], eigenValues[1]);
157 
158  bool numericalZero =false;
159  if (std::fabs(p[2])<1e-5)
160  numericalZero = true;
161 
162  if (!numericalZero)
163  {
164  double norm, v1,v2,v3;
165  // double p2=p[2]!=0?p[2]:1e-10, p1=p[1]!=0?p[1]:1e-10;
166  double p2=p[2], p1=p[1];
167 
168  double cebf = -p2*p[4]+p1*p[5];
169  double cdbe = -p2*p[3]+p1*p[4];
170  for ( int i = 0; i < 3; ++i )
171  {
172  v1 = -(p[5]-eigenValues[i])*(cdbe+p2*eigenValues[i]) + p[4]*(cebf-p1*eigenValues[i]);
173  v2 = -p2*(cebf-p1*eigenValues[i]);
174  v3 = p2*(cdbe+p2*eigenValues[i]);
175  norm = std::sqrt(v1*v1 + v2*v2 + v3*v3);
176  if (std::fabs(p2)<1e-5 || std::fabs(cdbe+p2*eigenValues[i])<1e-5 )
177  {
178  // NOTE: when p[2] or cdbe+p2*eigenValues[i] equal 0, we have divide by zero problem.
179  numericalZero=true;
180  break;
181  }
182  eigenVectors(0,i) = v1/norm;
183  eigenVectors(1,i) = v2/norm;
184  eigenVectors(2,i) = v3/norm;
185  }
186  }
187 
188  if (numericalZero)
189  {
190  vnl_diag_matrix<TPrecision> val;
191  vnl_matrix<double> vec;
192  GetEigenValuesVectors(val, vec);
193  for ( int i = 0; i < 3; ++i )
194  for ( int j = 0; j < 3; ++j )
195  eigenVectors(i,j) = vec(i,j);
196  }
197 }
198 
199 template<class TPrecision>
200 template<class TArrayType >
201 void
202 DiffusionTensor<TPrecision>
203 ::SetEigenValues(const TArrayType& array)
204 {
205  if (this->IsDiagonal())
206  {
207  for ( int i = 0; i < NDimension; i += 1 )
208  (*this)(i,i) = array[i];
209  }
210  else
211  {
212  vnl_diag_matrix<TPrecision> eigenValues(NDimension);
213  vnl_matrix<TPrecision> eigenVectors(NDimension, NDimension);
214  this->GetEigenValuesVectorsAnalytic(eigenValues, eigenVectors);
215  for ( int i = 0; i < NDimension; i += 1 )
216  eigenValues[i] = array[i];
217  *this = eigenVectors * eigenValues.asMatrix() * eigenVectors.transpose();
218  }
219 }
220 
221 template<class TPrecision>
222 template<class TMatrixType >
223 void
224 DiffusionTensor<TPrecision>
225 ::SetEigenVectors(const TMatrixType& vectors)
226 {
227  vnl_matrix<TPrecision> eigenVectors(NDimension, NDimension);
228  for ( int i = 0; i < NDimension; i += 1 )
229  for ( int j = 0; j < NDimension; j += 1 )
230  eigenVectors(i,j) = (TPrecision)vectors(i,j);
231  if (this->IsDiagonal())
232  *this = eigenVectors * (*this) * eigenVectors.transpose();
233  else
234  {
235  vnl_diag_matrix<TPrecision> eigenValues(NDimension);
236  this->GetEigenValuesVectorsAnalytic(eigenValues, eigenVectors);
237  *this = eigenVectors * eigenValues.asMatrix() * eigenVectors.transpose();
238  }
239 }
240 
241 
242 template<class TPrecision >
243 vnl_matrix<TPrecision>
244 DiffusionTensor<TPrecision>
245 ::GetVnlMatrix(void) const
246 {
247  TPrecision* block = new TPrecision[9];
248  for(unsigned int i=0;i<NDimension;i++)
249  for(unsigned int j=0;j<NDimension;j++)
250  block[NDimension*i+j] = (*this)(i,j);
251 
252  vnl_matrix< TPrecision > result = vnl_matrix< TPrecision >(block, NDimension, NDimension);
253  delete [] block;
254  return result;
255 }
256 
257 template<class TPrecision>
258 void
259 DiffusionTensor<TPrecision>
260 ::SetVnlMatrix( const vnl_matrix<TPrecision> & mat)
261 {
262  for(unsigned int i=0; i<NDimension; i++)
263  {
264  for(unsigned int j=i; j<NDimension; j++)
265  (*this)[NDimension*i-(i+1)*i/2+j] = mat(i,j);
266  }
267 }
268 
269 template<class TPrecision>
270 DiffusionTensor<TPrecision>&
271 DiffusionTensor<TPrecision>
272 ::Rotate (const typename Self::MatrixType& R)
273 {
274  this->SetVnlMatrix ( this->GetRotate (R.GetVnlMatrix() ) );
275  return *this;
276 }
277 
278 template<class TPrecision>
279 DiffusionTensor<TPrecision>
280 DiffusionTensor<TPrecision>
281 ::GetRotate (const typename Self::MatrixType& R) const
282 {
283  Self res;
284  res.SetVnlMatrix ( this->GetRotate (R.GetVnlMatrix() ) );
285  return res;
286 }
287 
288 
289 template<class TPrecision>
290 vnl_matrix<TPrecision>
291 DiffusionTensor<TPrecision>
292 ::GetRotate (const vnl_matrix<TPrecision> & R) const
293 {
294  vnl_matrix<TPrecision> res;
295  res = ( R*(this->GetVnlMatrix())*R.transpose() );
296  return res;
297 }
298 
299 template < class TPrecision>
300 double
301 DiffusionTensor<TPrecision>
302 ::GetQuadraticForm ( const double x, const double y, const double z) const
303 {
304  return (*this)[0]*x*x + 2*(*this)[1]*x*y + 2*(*this)[2]*x*z +
305  (*this)[3]*y*y + 2*(*this)[4]*y*z + (*this)[5]*z*z;
306 }
307 
308 template < class TPrecision>
309 template < class TVectorType >
310 void
311 DiffusionTensor<TPrecision>
312 ::GetSphericalSamples ( TVectorType& samples, const utl::NDArray<TPrecision,2>& gradients) const
313 {
314  for ( int i = 0; i < gradients.Rows(); i += 1 )
315  {
316  TPrecision x=gradients(i,0), y=gradients(i,1), z=gradients(i,2);
317  double quadratic = GetQuadraticForm(x,y,z);
318  samples[i] = GetQuadraticForm(x,y,z);
319  }
320 }
321 
322 template < class TPrecision>
323 template < class TVectorType >
324 void
325 DiffusionTensor<TPrecision>
326 ::GetDWISamples ( TVectorType& dwisignal, const utl::NDArray<TPrecision,2>& gradients, const std::vector<TPrecision>& bValues) const
327 {
328  utlAssert(gradients.Rows()==bValues.size(), "wrong size! gradients.Rows()="<< gradients.Rows() << ", bValues.size()="<<bValues.size());
329 
330  for ( int i = 0; i < gradients.Rows(); i += 1 )
331  {
332  TPrecision x=gradients(i,0), y=gradients(i,1), z=gradients(i,2);
333  double quadratic = GetQuadraticForm(x,y,z);
334  dwisignal[i] = std::exp(-1.0*bValues[i]*quadratic);
335  }
336 } // ----- end of method itkDiffusionTensor<TPrecision>::GetDWISignal -----
337 
338 template < class TPrecision >
339 template < class TVectorType >
340 void
341 DiffusionTensor<TPrecision>
342 ::GetODFSamples ( TVectorType& odf, const utl::NDArray<TPrecision,2>& gradients, const int& odfOrder, const bool& isNormalize) const
343 {
344  DiffusionTensor<TPrecision> D_inv = this->Inv();
345  TPrecision D_det = this->GetDeterminant();
346  for ( int i = 0; i < gradients.Rows(); i += 1 )
347  {
348  TPrecision x=gradients(i,0), y=gradients(i,1), z=gradients(i,2);
349  TPrecision quadratic_inv = D_inv(0,0)*x*x + 2*D_inv(0,1)*x*y + 2*D_inv(0,2)*x*z +
350  D_inv(1,1)*y*y + 2*D_inv(1,2)*y*z + D_inv(2,2)*z*z;
351 
352  odf[i] = 1.0 / ( std::pow(D_det, 0.5) * std::pow(quadratic_inv,0.5*(odfOrder+1)) ) ;
353  if (odfOrder==2)
354  odf[i] *= 1.0/(4*M_PI);
355  }
356 
357  // when odfOrder==2, it is naturally normalized.
358  if (isNormalize && odfOrder!=2)
359  {
360  TPrecision normFactor = 4*M_PI*utl::GetSumOfVector<TVectorType>(odf,gradients.Rows()) / gradients.Rows();
361  if (normFactor!=0)
362  {
363  for ( int i = 0; i < gradients.Rows(); i += 1 )
364  odf[i] /= normFactor;
365  }
366  }
367 } // ----- end of method itkDiffusionTensor<TPrecision>::GetDWISignal -----
368 
369 template < class TPrecision >
370 template < class TVectorType >
371 void
372 DiffusionTensor<TPrecision>
373 ::GetEAPSamples ( TVectorType& eap, const utl::NDArray<TPrecision,2>& gradients, const std::vector<TPrecision>& rValues, const TPrecision& tau) const
374 {
375  utlAssert(gradients.Rows()==rValues.size(), "wrong size! gradients.Rows()="<< gradients.Rows() << ", rValues.size()="<<rValues.size());
376 
377  DiffusionTensor<TPrecision> D_inv = this->Inv();
378  TPrecision D_det = this->GetDeterminant();
379  for ( int i = 0; i < gradients.Rows(); i += 1 )
380  {
381  TPrecision x=gradients(i,0), y=gradients(i,1), z=gradients(i,2);
382  TPrecision quadratic_inv = D_inv(0,0)*x*x + 2*D_inv(0,1)*x*y + 2*D_inv(0,2)*x*z +
383  D_inv(1,1)*y*y + 2*D_inv(1,2)*y*z + D_inv(2,2)*z*z;
384  eap[i] = 1.0 / (std::pow(4.0*M_PI*tau,(double)1.5)*std::sqrt(D_det) ) * std::exp(-1.0/(4.0*tau)*rValues[i]*rValues[i]*quadratic_inv);
385  }
386 } // ----- end of method itkDiffusionTensor<TPrecision>::GetDWISignal -----
387 
388 template < class TPrecision >
389 typename DiffusionTensor<TPrecision>::RealValueType
390 DiffusionTensor<TPrecision>
391 ::GetReturnToOrigin (const TPrecision tau) const
392 {
393  RealValueType D_det = this->GetDeterminant();
394  return 1.0 / (std::pow(4*M_PI*tau,(double)1.5)*std::sqrt(D_det) );
395 } // ----- end of method itkDiffusionTensor<TPrecision>::GetDWISignal -----
396 
397 template < class TPrecision >
398 typename DiffusionTensor<TPrecision>::RealValueType
399 DiffusionTensor<TPrecision>
400 ::GetMeanSquaredDisplacement (const TPrecision tau) const
401 {
402  return this->GetTrace()*2.0*tau;
403 } // ----- end of method itkDiffusionTensor<TPrecision>::GetDWISignal -----
404 
405 template < class TPrecision >
406 typename DiffusionTensor<TPrecision>::RealValueType
407 DiffusionTensor<TPrecision>
408 ::GetGA () const
409 {
410  vnl_diag_matrix<TPrecision> eigenValues(NDimension);
411  vnl_matrix<TPrecision> eigenVectors(NDimension, NDimension);
412  this->GetEigenValuesVectorsAnalytic(eigenValues, eigenVectors);
413 
414  vnl_diag_matrix<TPrecision> logeigenValues(NDimension);
415  for ( int i = 0; i < NDimension; i += 1 )
416  logeigenValues[i] = std::log(eigenValues[i]);
417 
418  TPrecision meanEigenValues=0.0, ga=0.0;
419  for ( int i = 0; i < NDimension; i += 1 )
420  meanEigenValues += logeigenValues[i];
421  meanEigenValues /= NDimension;
422  for ( int i = 0; i < NDimension; i += 1 )
423  {
424  double tmp = logeigenValues[i]-meanEigenValues;
425  ga += tmp*tmp;
426  }
427  return std::sqrt(ga);
428 }
429 
430 template < class TPrecision >
431 typename DiffusionTensor<TPrecision>::RealValueType
432 DiffusionTensor<TPrecision>
433 ::GetFA () const
434 {
435  return this->GetFractionalAnisotropy();
436 }
437 
438 template < class TPrecision >
439 typename DiffusionTensor<TPrecision>::RealValueType
440 DiffusionTensor<TPrecision>
441 ::GetADC () const
442 {
443  return this->GetTrace()/RealValueType(NDimension);
444 }
445 
446 template < class TPrecision >
447 typename DiffusionTensor<TPrecision>::RealValueType
448 DiffusionTensor<TPrecision>
449 ::GetRA () const
450 {
451  return this->GetRelativeAnisotropy();
452 }
453 
454 template < class TPrecision >
455 typename DiffusionTensor<TPrecision>::RealValueType
456 DiffusionTensor<TPrecision>
457 ::GetMODE () const
458 {
459  double md = GetMD();
460 
461  Self dTensor(*this);
462  dTensor[0] -= md;
463  dTensor[3] -= md;
464  dTensor[5] -= md;
465 
466  double dtNorm = dTensor.GetNorm();
467  dTensor /= dtNorm;
468  return 3.0*std::sqrt(6.0)*dTensor.GetDeterminant();
469 }
470 
471 template<class TPrecision>
472 bool
473 DiffusionTensor<TPrecision>
474 ::IsPositive () const
475 {
476  utlAssert (this->IsFinite(), "Tensor is not finite");
477  EigenValuesArrayType eigenValues;
478  this->ComputeEigenValues(eigenValues);
479  if( eigenValues[0] <= NumericTraits<double>::Zero )
480  return false;
481  return true;
482 }
483 
484 template<class TPrecision>
485 void
486 DiffusionTensor<TPrecision>
487 ::Positivize () const
488 {
489  if (IsPositive())
490  return;
491 
492  vnl_diag_matrix<TPrecision> eigenValues(NDimension);
493  vnl_matrix<TPrecision> eigenVectors(NDimension, NDimension);
494  this->GetEigenValuesVectorsAnalytic(eigenValues, eigenVectors);
495 
496  for ( int i = 0; i < NDimension; i += 1 )
497  {
498  if (eigenValues[i]<0)
499  eigenValues[i]=0;
500  }
501  *this = eigenVectors* eigenValues.asMatrix()* eigenVectors.transpose();
502 }
503 
504 template<class TPrecision>
505 bool
506 DiffusionTensor<TPrecision>
507 ::IsZero (const double eps) const
508 {
509  for ( int i = 0; i < NDegreesOfFreedom; i += 1 )
510  {
511  if (std::fabs((*this)[i])>eps)
512  return false;
513  }
514  return true;
515 }
516 
517 template<class TPrecision>
518 bool
519 DiffusionTensor<TPrecision>
520 ::IsDiagonal (const double eps) const
521 {
522  for ( int i = 0; i < NDimension; i += 1 )
523  for ( int j = i+1; j < NDimension; j += 1 )
524  {
525  if (std::fabs((*this)(i,j))>eps)
526  return false;
527  }
528  return true;
529 }
530 
531 template<class TPrecision>
532 bool
533 DiffusionTensor<TPrecision>
534 ::IsFinite () const
535 {
536  return this->GetVnlMatrix().is_finite();
537 }
538 
539 template<class TPrecision>
540 bool
541 DiffusionTensor<TPrecision>
542 ::HasNans () const
543 {
544  return this->GetVnlMatrix().has_nans();
545 }
546 
547 template<class TPrecision>
548 double
549 DiffusionTensor<TPrecision>
550 ::GetDeterminant () const
551 {
552  const TPrecision* p = this->GetDataPointer();
553  return -p[2]*p[2]*p[3] + 2.0*p[1]*p[2]*p[4] - p[0]*p[4]*p[4] - p[1]*p[1]*p[5] + p[0]*p[3]*p[5];
554 }
555 
556 template<class TPrecision>
557 double
558 DiffusionTensor<TPrecision>
559 ::GetNorm () const
560 {
561  const TPrecision* p = this->GetDataPointer();
562  return std::sqrt(p[0]*p[0]+2.0*p[1]*p[1]+2.0*p[2]*p[2] + p[3]*p[3]+ 2.0*p[4]*p[4] + p[5]*p[5]);
563 }
564 
565 template<class TPrecision>
566 DiffusionTensor<TPrecision>
567 DiffusionTensor<TPrecision>
568 ::Log () const
569 {
570  utlException (!this->IsFinite(), "Tensor is not finite");
571  Self result;
572  vnl_diag_matrix<TPrecision> eigenValues(NDimension);
573  vnl_matrix<TPrecision> eigenVectors(NDimension, NDimension);
574  this->GetEigenValuesVectorsAnalytic(eigenValues, eigenVectors);
575  for(unsigned int i=0;i<NDimension;i++)
576  {
577  utlException (eigenValues[i] < 0.0, "Negative eigenvalue encountered.");
578  eigenValues[i] = vcl_log (eigenValues[i]);
579  }
580 
581  vnl_matrix<TPrecision> vnl_result = eigenVectors * eigenValues.asMatrix() * eigenVectors.transpose();
582  result.SetVnlMatrix ( vnl_result );
583  return result;
584 }
585 
586 template<class TPrecision>
587 DiffusionTensor<TPrecision>
588 DiffusionTensor<TPrecision>
589 ::Exp () const
590 {
591  utlAssert(this->IsFinite(),"Tensor is not finite");
592  Self result;
593  vnl_diag_matrix<TPrecision> eigenValues(NDimension);
594  vnl_matrix<TPrecision> eigenVectors(NDimension, NDimension);
595  this->GetEigenValuesVectorsAnalytic(eigenValues, eigenVectors);
596  for(unsigned int i=0;i<NDimension;i++)
597  eigenValues[i] = vcl_exp (eigenValues[i]);
598 
599  vnl_matrix<TPrecision> vnl_result = eigenVectors * eigenValues.asMatrix() * eigenVectors.transpose();
600  result.SetVnlMatrix ( vnl_result );
601  return result;
602 }
603 
604 template<class TPrecision>
605 DiffusionTensor<TPrecision>
606 DiffusionTensor<TPrecision>
607 ::Pow (const double& n) const
608 {
609  utlAssert (this->IsFinite(), "Tensor is not finite");
610  Self result;
611 
612  vnl_diag_matrix<TPrecision> eigenValues(NDimension);
613  vnl_matrix<TPrecision> eigenVectors(NDimension, NDimension);
614  this->GetEigenValuesVectorsAnalytic(eigenValues, eigenVectors);
615  // std::cout << *this << std::endl << std::flush;
616  // std::cout << "eigenValues = " << eigenValues << std::endl << std::flush;
617  // std::cout << "eigenVectors = " << eigenVectors << std::endl << std::flush;
618 
619  // typedef vnl_symmetric_eigensystem< TPrecision > SymEigenSystemType;
620  // SymEigenSystemType eig (this->GetVnlMatrix());
621  // std::cout << "eig.D = " << eig.D << std::endl << std::flush;
622  // std::cout << "eig.V = " << eig.V << std::endl << std::flush;
623 
624  for(unsigned int i=0;i<NDimension;i++)
625  {
626  utlException(n<0 && std::abs(eigenValues[i])<=1e-10, "pow="<<n<<", eigenValues["<<i<<"]="<<eigenValues[i] << " close to zero");
627  eigenValues[i] = static_cast<TPrecision> (vcl_pow (static_cast<double>(eigenValues[i]), n));
628  }
629 
630  vnl_matrix<TPrecision> vnl_result = eigenVectors * eigenValues.asMatrix() * eigenVectors.transpose();
631  result.SetVnlMatrix ( vnl_result );
632  return result;
633 }
634 
635 template<class TPrecision>
636 DiffusionTensor<TPrecision>
637 DiffusionTensor<TPrecision>
638 ::Inv (void) const
639 {
640  double det = GetDeterminant();
641  utlGlobalException(det==0, "det=0, cannot be inverted");
642  double detInv=1.0/det;
643  const TPrecision* p = this->GetDataPointer();
644  DiffusionTensor<TPrecision> result;
645  result[0] = (-p[4]*p[4] + p[3]*p[5]) *detInv;
646  result[1] = ( p[2]*p[4] - p[1]*p[5]) *detInv;
647  result[2] = (-p[2]*p[3] + p[1]*p[4]) *detInv;
648  result[3] = (-p[2]*p[2] + p[0]*p[5]) *detInv;
649  result[4] = ( p[1]*p[2] - p[0]*p[4]) *detInv;
650  result[5] = (-p[1]*p[1] + p[0]*p[3]) *detInv;
651  return result;
652 }
653 
654 
655 template<class TPrecision>
656 DiffusionTensor<TPrecision>
657 DiffusionTensor<TPrecision>
658 ::Sqrt (void) const
659 {
660  return this->Pow (0.5);
661 }
662 
663 
664 template<class TPrecision>
665 DiffusionTensor<TPrecision>
666 DiffusionTensor<TPrecision>
667 ::InvSqrt () const
668 {
669  return this->Pow (-0.5);
670 }
671 
672 template<typename TPrecision>
673 TPrecision
674 DiffusionTensor<TPrecision>
675 ::EuclideanDistance(const DiffusionTensor<TPrecision>& X) const
676 {
677  return (*this - X).GetInnerScalarProduct();
678 }
679 
680 template<typename TPrecision>
681 TPrecision
682 DiffusionTensor<TPrecision>
683 ::KLDistance(const DiffusionTensor<TPrecision>& X) const
684 {
685  const DiffusionTensor<TPrecision> ithis = this->Inv();
686  const DiffusionTensor<TPrecision> iX = X.Inv();
687  const TPrecision tr = (ithis*X + iX*(*this)).GetTrace();
688  return 0.5*std::sqrt(tr - 2*NDimension);
689 }
690 
691 template<typename TPrecision>
692 TPrecision
693 DiffusionTensor<TPrecision>
694 ::GeodesicDistance(const DiffusionTensor<TPrecision>& X) const
695 {
696  const DiffusionTensor<TPrecision> isqrtT = this->InvSqrt();
697  const DiffusionTensor<TPrecision> A = (isqrtT * X * isqrtT).Log();
698  return std::sqrt(0.5*(A*A).GetTrace());
699 }
700 
701 template<typename TPrecision>
702 TPrecision
703 DiffusionTensor<TPrecision>
704 ::LogEucDistance(const DiffusionTensor<TPrecision>& X) const
705 {
706  DiffusionTensor<TPrecision> logThis = this->Log();
707  DiffusionTensor<TPrecision> logX = X.Log();
708  return (logThis-logX).GetInnerScalarProduct();
709 }
710 
711 
712 }
713 
714 #endif
itk::DiffusionTensor::RealValueType
Superclass::RealValueType RealValueType
Definition: itkDiffusionTensor.h:64
utl::NDArray
NDArray is a N-Dimensional array class (row-major, c version)
Definition: utlFunctors.h:131
itk::DiffusionTensor::GetEAPSamples
void GetEAPSamples(TVectorType &vec, const utl::NDArray< TPrecision, 2 > &gradients, const std::vector< TPrecision > &rValues, const TPrecision &tau) const
Definition: itkDiffusionTensor.hxx:373
itk::DiffusionTensor::GetMODE
RealValueType GetMODE() const
Definition: itkDiffusionTensor.hxx:457
utl.h
helper functions specifically used in dmritool
itk::DiffusionTensor::GeodesicDistance
TPrecision GeodesicDistance(const DiffusionTensor< TPrecision > &X) const
Definition: itkDiffusionTensor.hxx:694
itk::DiffusionTensor::LogEucDistance
TPrecision LogEucDistance(const DiffusionTensor< TPrecision > &X) const
Definition: itkDiffusionTensor.hxx:704
utl::Pow
auto Pow(const TLeft &lhs, const TRight &rhs) -> decltype(utl::F< utl::Functor::Pow< Expr2ValueType< TLeft, TRight >> >(lhs, rhs))
Definition: utlFunctors.h:385
itk::DiffusionTensor::EuclideanDistance
TPrecision EuclideanDistance(const DiffusionTensor< TPrecision > &X) const
Definition: itkDiffusionTensor.hxx:675
itk::DiffusionTensor::GetRA
RealValueType GetRA() const
Definition: itkDiffusionTensor.hxx:449
itk::DiffusionTensor::GetDWISamples
void GetDWISamples(TVectorType &vec, const utl::NDArray< TPrecision, 2 > &gradients, const std::vector< TPrecision > &bValues) const
Definition: itkDiffusionTensor.hxx:326
utl::Log
auto Log(const ExprT &expr) -> decltype(utl::F< utl::Functor::Log< typename ExprT::ValueType > >(expr))
Definition: utlFunctors.h:359
utlException
#define utlException(cond, expout)
Definition: utlCoreMacro.h:548
itk::DiffusionTensor::GetReturnToOrigin
RealValueType GetReturnToOrigin(const TPrecision tau) const
Definition: itkDiffusionTensor.hxx:391
itk::DiffusionTensor::Pow
Self Pow(const double &n) const
Definition: itkDiffusionTensor.hxx:607
itk::DiffusionTensor::GetEigenValuesVectors
void GetEigenValuesVectors(vnl_diag_matrix< TPrecision > &eigenValues, vnl_matrix< TPrecision > &eigenVectors) const
Return an array containing EigenValues in ascending order, and a matrix containing the corresponding ...
Definition: itkDiffusionTensor.hxx:84
itk::DiffusionTensor::GetEigenValuesVectorsAnalytic
void GetEigenValuesVectorsAnalytic(TArrayType &eigenValues, TMatrixType &eigenVectors) const
analytic way to calculate eigenValues (in ascending order) and eigenVectors. In mathematica, Eigenvectors[( { {a, b, c}, {b, d, e}, {c, e, f} } )]
Definition: itkDiffusionTensor.hxx:100
itk::DiffusionTensor::GetDeterminant
double GetDeterminant() const
Definition: itkDiffusionTensor.hxx:550
utl::PolynomialRoot
std::vector< std::complex< double > > PolynomialRoot(const std::vector< double > &coef)
Definition: utlMath.h:1100
itk::DiffusionTensor::GetFA
RealValueType GetFA() const
Definition: itkDiffusionTensor.hxx:433
itk::DiffusionTensor::Rotate
Self & Rotate(const typename Self::MatrixType &mat)
Definition: itkDiffusionTensor.hxx:272
itk::DiffusionTensor::operator()
ValueType & operator()(unsigned int row, unsigned int col)
Definition: itkDiffusionTensor.hxx:64
itk::DiffusionTensor::GetADC
RealValueType GetADC() const
Definition: itkDiffusionTensor.hxx:441
spams::abs
T abs(const T x)
template version of the fabs function
M_PI
#define M_PI
Definition: utlCoreMacro.h:57
utlGlobalException
#define utlGlobalException(cond, expout)
Definition: utlCoreMacro.h:372
itk::DiffusionTensor::EigenValuesArrayType
Superclass::EigenValuesArrayType EigenValuesArrayType
Definition: itkDiffusionTensor.h:60
itk::DiffusionTensor::SetEigenValues
void SetEigenValues(const TArrayType &array)
Definition: itkDiffusionTensor.hxx:203
itk::DiffusionTensor::GetVnlMatrix
vnl_matrix< TPrecision > GetVnlMatrix() const
Definition: itkDiffusionTensor.hxx:245
itk::DiffusionTensor::EigenVectorsMatrixType
Superclass::EigenVectorsMatrixType EigenVectorsMatrixType
Definition: itkDiffusionTensor.h:61
utlAssert
#define utlAssert(cond, expout)
Definition: utlCoreMacro.h:550
itk::DiffusionTensor::SetVnlMatrix
void SetVnlMatrix(const vnl_matrix< TPrecision > &mat)
Definition: itkDiffusionTensor.hxx:260
itk::DiffusionTensor::IsZero
bool IsZero(const double eps=1e-10) const
Definition: itkDiffusionTensor.hxx:507
itk::DiffusionTensor::IsDiagonal
bool IsDiagonal(const double eps=1e-10) const
Definition: itkDiffusionTensor.hxx:520
itk::DiffusionTensor::GetGA
RealValueType GetGA() const
Definition: itkDiffusionTensor.hxx:408
itk::DiffusionTensor::GetSphericalSamples
void GetSphericalSamples(TVectorType &vec, const utl::NDArray< TPrecision, 2 > &gradients) const
Definition: itkDiffusionTensor.hxx:312
itk::DiffusionTensor::GetMeanSquaredDisplacement
RealValueType GetMeanSquaredDisplacement(const TPrecision tau) const
Definition: itkDiffusionTensor.hxx:400
itk::DiffusionTensor::GetODFSamples
void GetODFSamples(TVectorType &vec, const utl::NDArray< TPrecision, 2 > &gradients, const int &odfOrder=2, const bool &isNormalize=false) const
Definition: itkDiffusionTensor.hxx:342
itk::DiffusionTensor::GetRotate
Self GetRotate(const typename Self::MatrixType &mat) const
Definition: itkDiffusionTensor.hxx:281
itk::DiffusionTensor::KLDistance
TPrecision KLDistance(const DiffusionTensor< TPrecision > &X) const
Definition: itkDiffusionTensor.hxx:683
itk::DiffusionTensor::GetQuadraticForm
double GetQuadraticForm(const double x, const double y, const double z) const
Definition: itkDiffusionTensor.hxx:302
itk::DiffusionTensor
tensor with some useful functions
Definition: itkDiffusionTensor.h:46
itk::DiffusionTensor::ValueType
Superclass::ValueType ValueType
Definition: itkDiffusionTensor.h:57
itk::DiffusionTensor::SetEigenVectors
void SetEigenVectors(const TMatrixType &vectors)
Definition: itkDiffusionTensor.hxx:225
Generated by   doxygen 1.8.11