ANIMA  4.0
animaTensorCompartment.cxx
Go to the documentation of this file.
2 
3 #include <itkSymmetricEigenAnalysis.h>
4 #include <animaMCMConstants.h>
5 
6 namespace anima
7 {
8 
9 double TensorCompartment::GetFourierTransformedDiffusionProfile(double smallDelta, double bigDelta, double gradientStrength, const Vector3DType &gradient)
10 {
11  this->UpdateDiffusionTensor();
12 
13  double quadForm = 0;
14 
15  for (unsigned int i = 0;i < m_SpaceDimension;++i)
16  {
17  quadForm += m_DiffusionTensor(i,i) * gradient[i] * gradient[i];
18  for (unsigned int j = i+1;j < m_SpaceDimension;++j)
19  quadForm += 2 * m_DiffusionTensor(i,j) * gradient[i] * gradient[j];
20  }
21 
22  double bValue = anima::GetBValueFromAcquisitionParameters(smallDelta, bigDelta, gradientStrength);
23 
24  return std::exp(- bValue * quadForm);
25 }
26 
27 TensorCompartment::ListType &TensorCompartment::GetSignalAttenuationJacobian(double smallDelta, double bigDelta, double gradientStrength, const Vector3DType &gradient)
28 {
29  this->UpdateDiffusionTensor();
30 
32 
33  double signalAttenuation = this->GetFourierTransformedDiffusionProfile(smallDelta, bigDelta, gradientStrength, gradient);
34  double innerProd1 = anima::ComputeScalarProduct(gradient, m_EigenVector1);
35  double innerProd2 = anima::ComputeScalarProduct(gradient, m_EigenVector2);
36 
37  double DgTe1DTheta = m_CosTheta * (gradient[0] * m_CosPhi + gradient[1] * m_SinPhi) - gradient[2] * m_SinTheta;
38  double DgTe1DPhi = m_SinTheta * (gradient[1] * m_CosPhi - gradient[0] * m_SinPhi);
39 
40  double DgTe2DTheta = m_SinAlpha * innerProd1;
41  double DgTe2DPhi = gradient[0] * (m_CosTheta * m_SinPhi * m_SinAlpha - m_CosPhi * m_CosAlpha) - gradient[1] * (m_SinPhi * m_CosAlpha + m_CosTheta * m_CosPhi * m_SinAlpha);
42  double DgTe2DAlpha = gradient[0] * (m_SinPhi * m_SinAlpha - m_CosTheta * m_CosPhi * m_CosAlpha) - gradient[1] * (m_CosPhi * m_SinAlpha + m_CosTheta * m_SinPhi * m_CosAlpha) + gradient[2] * m_SinTheta * m_CosAlpha;
43 
44  double diffAxialRadial2 = this->GetAxialDiffusivity() - this->GetRadialDiffusivity2();
45  double diffRadialDiffusivities = this->GetRadialDiffusivity1() - this->GetRadialDiffusivity2();
46 
47  // Derivative w.r.t. theta
48  double bValue = anima::GetBValueFromAcquisitionParameters(smallDelta, bigDelta, gradientStrength);
49  m_JacobianVector[0] = -2.0 * bValue * (diffAxialRadial2 * innerProd1 * DgTe1DTheta + diffRadialDiffusivities * innerProd2 * DgTe2DTheta) * signalAttenuation;
50 
51  // Derivative w.r.t. phi
52  m_JacobianVector[1] = -2.0 * bValue * (diffAxialRadial2 * innerProd1 * DgTe1DPhi + diffRadialDiffusivities * innerProd2 * DgTe2DPhi) * signalAttenuation;
53 
54  // Derivative w.r.t. alpha
55  m_JacobianVector[2] = -2.0 * bValue * diffRadialDiffusivities * innerProd2 * DgTe2DAlpha * signalAttenuation;
56 
57  if (m_EstimateDiffusivities)
58  {
59  // Derivative w.r.t. to d1
60  m_JacobianVector[3] = - bValue * innerProd1 * innerProd1 * signalAttenuation;
61 
62  // Derivative w.r.t. to d2
63  m_JacobianVector[4] = - bValue * (innerProd1 * innerProd1 + innerProd2 * innerProd2) * signalAttenuation;
64 
65  // Derivative w.r.t. to d3
66  m_JacobianVector[5] = - bValue * signalAttenuation;
67  }
68 
69  return m_JacobianVector;
70 }
71 
73 {
75 
76  double resVal = - 1.5 * std::log(2.0 * M_PI) - 0.5 * std::log(m_TensorDeterminant);
77 
78  double quadForm = 0;
79 
80  for (unsigned int i = 0;i < m_SpaceDimension;++i)
81  {
82  quadForm += m_InverseDiffusionTensor(i,i) * sample[i] * sample[i];
83  for (unsigned int j = i+1;j < m_SpaceDimension;++j)
84  quadForm += 2 * m_InverseDiffusionTensor(i,j) * sample[i] * sample[j];
85  }
86 
87  resVal -= quadForm / 2.0;
88 
89  return resVal;
90 }
91 
93 {
94  if (num != this->GetOrientationTheta())
95  {
96  m_ModifiedTensor = true;
97  m_ModifiedAngles = true;
99  }
100 }
101 
103 {
106 }
107 
109 {
110  if (num != this->GetOrientationPhi())
111  {
112  m_ModifiedTensor = true;
113  m_ModifiedAngles = true;
115  }
116 }
117 
119 {
122 }
123 
125 {
126  if (num != this->GetPerpendicularAngle())
127  {
128  m_ModifiedTensor = true;
129  m_ModifiedAngles = true;
131  }
132 }
133 
135 {
138 }
139 
141 {
142  if (num != this->GetAxialDiffusivity())
143  {
144  m_ModifiedTensor = true;
146  }
147 }
148 
150 {
153 }
154 
156 {
157  if (num != this->GetRadialDiffusivity1())
158  {
159  m_ModifiedTensor = true;
161  }
162 }
163 
165 {
168 }
169 
171 {
172  if (num != this->GetRadialDiffusivity2())
173  {
174  m_ModifiedTensor = true;
176  }
177 }
178 
180 {
183 }
184 
186 {
187  if (params.size() != this->GetNumberOfParameters())
188  return;
189 
190  this->SetOrientationTheta(params[0]);
191  this->SetOrientationPhi(params[1]);
192  this->SetPerpendicularAngle(params[2]);
193 
194  if (m_EstimateDiffusivities)
195  {
196  this->SetAxialDiffusivity(params[3] + params[4] + params[5]);
197  this->SetRadialDiffusivity1(params[4] + params[5]);
198  this->SetRadialDiffusivity2(params[5]);
199  }
200 }
201 
203 {
205 
209 
210  if (m_EstimateDiffusivities)
211  {
215  }
216 
217  return m_ParametersVector;
218 }
219 
221 {
224 
225  if (m_EstimateDiffusivities)
226  {
229  }
230 
232 }
233 
235 {
237 
241 
242  if (m_EstimateDiffusivities)
243  {
247  }
248 
250 }
251 
253 {
254  if (m_EstimateDiffusivities == arg)
255  return;
256 
257  m_EstimateDiffusivities = arg;
258  m_ChangedConstraints = true;
259 }
260 
262 {
263  if (compartmentVector.GetSize() != this->GetCompartmentSize())
264  itkExceptionMacro("The input vector size does not match the size of the compartment");
265 
266  anima::GetTensorFromVectorRepresentation(compartmentVector,m_WorkVnlMatrix1,m_SpaceDimension);
267  m_DiffusionTensor = m_WorkVnlMatrix1;
268 
269  m_TensorDeterminant = vnl_determinant <double> (m_WorkVnlMatrix1);
270 
271  if (m_TensorDeterminant <= 0)
272  m_TensorDeterminant = 0;
273 
274  m_ModifiedTensor = false;
275  m_ModifiedAngles = false;
276  m_UpdateInverseTensor = true;
277  m_UpdatedCompartment = true;
278 }
279 
281 {
282  if (!m_UpdatedCompartment)
283  return;
284 
285  vnl_matrix <double> eVecs(m_SpaceDimension,m_SpaceDimension,0);
286  vnl_diag_matrix <double> eVals(m_SpaceDimension);
287 
288  itk::SymmetricEigenAnalysis < Matrix3DType,vnl_diag_matrix <double>,vnl_matrix <double> > eigSys(m_SpaceDimension);
289 
290  eigSys.ComputeEigenValuesAndVectors(m_DiffusionTensor,eVals,eVecs);
291 
292  Vector3DType e1, sphDir;
293 
294  for (unsigned int i = 0;i < m_SpaceDimension;++i)
295  e1[i] = eVecs(2,i);
296 
298 
299  double theta = sphDir[0];
300  double phi = sphDir[1];
301 
302  double alpha = 0;
303  if ((theta != 0)&&(theta != M_PI))
304  alpha = std::atan2(eVecs(1,2),- eVecs(0,2));
305  else
306  {
307  // Use two first components of 2nd eigenvector
308  // Infinite number of solutions for phi from main eigenvector, set it as alpha - pi/2
309  alpha = std::atan2(1 - eVecs(1,0), eVecs(1,1));
310  phi = alpha - M_PI / 2.0;
311  }
312 
313  if (alpha < 0)
314  alpha += 2.0 * M_PI;
315 
316  if (phi < 0)
317  phi += 2.0 * M_PI;
318 
319  m_UpdatedCompartment = false;
320 
321  this->SetOrientationTheta(theta);
322  this->SetOrientationPhi(phi);
323  this->SetPerpendicularAngle(alpha);
324 
325  double axialDiff = std::max(anima::MCMDiffusivityLowerBound,std::min(eVals[2],anima::MCMDiffusivityUpperBound));
326  double radialDiff1 = std::max(anima::MCMDiffusivityLowerBound,std::min(eVals[1],anima::MCMDiffusivityUpperBound));
327  double radialDiff2 = std::max(anima::MCMDiffusivityLowerBound,std::min(eVals[0],anima::MCMDiffusivityUpperBound));
328 
329  this->SetAxialDiffusivity(axialDiff);
330  this->SetRadialDiffusivity1(radialDiff1);
331  this->SetRadialDiffusivity2(radialDiff2);
332 }
333 
335 {
336  return 6;
337 }
338 
340 {
341  if (!m_ChangedConstraints)
342  return m_NumberOfParameters;
343 
344  m_NumberOfParameters = this->GetCompartmentSize();
345 
346  if (!m_EstimateDiffusivities)
347  m_NumberOfParameters -= 3;
348 
349  m_ChangedConstraints = false;
350  return m_NumberOfParameters;
351 }
352 
354 {
355  if (m_CompartmentVector.GetSize() != this->GetCompartmentSize())
356  m_CompartmentVector.SetSize(this->GetCompartmentSize());
357 
358  this->UpdateDiffusionTensor();
359 
360  anima::GetVectorRepresentation(m_DiffusionTensor.GetVnlMatrix().as_matrix(),m_CompartmentVector,this->GetCompartmentSize());
361  return m_CompartmentVector;
362 }
363 
364 void TensorCompartment::Reorient(vnl_matrix <double> &orientationMatrix, bool affineTransform)
365 {
366  this->UpdateDiffusionTensor();
367 
368  m_WorkVnlMatrix1 = m_DiffusionTensor.GetVnlMatrix().as_matrix();
369 
370  if (!affineTransform)
371  anima::RotateSymmetricMatrix(m_WorkVnlMatrix1,orientationMatrix,m_WorkVnlMatrix2);
372  else
373  {
374  // PPD re-orientation
375  vnl_matrix <double> ppdOrientationMatrix(3,3);
376  typedef itk::Matrix <double, 3, 3> EigVecMatrixType;
377  typedef vnl_vector_fixed <double,3> EigValVectorType;
378  itk::SymmetricEigenAnalysis < EigVecMatrixType, EigValVectorType, EigVecMatrixType> eigenComputer(3);
379  EigVecMatrixType eigVecs;
380  EigValVectorType eigVals;
381 
382  eigenComputer.ComputeEigenValuesAndVectors(m_WorkVnlMatrix1,eigVals,eigVecs);
383  anima::ExtractPPDRotationFromJacobianMatrix(orientationMatrix,ppdOrientationMatrix,eigVecs);
384  anima::RotateSymmetricMatrix(m_WorkVnlMatrix1,ppdOrientationMatrix,m_WorkVnlMatrix2);
385  }
386 
388  anima::GetVectorRepresentation(m_WorkVnlMatrix2,tmpVec);
389 
390  this->SetCompartmentVector(tmpVec);
391 }
392 
394 {
395  this->UpdateDiffusionTensor();
396  return m_DiffusionTensor;
397 }
398 
400 {
401  if (!m_ModifiedTensor)
402  return;
403 
404  this->UpdateAngleConfiguration();
405  double axialDiff = this->GetAxialDiffusivity();
406  double radialDiff1 = this->GetRadialDiffusivity1();
407  double radialDiff2 = this->GetRadialDiffusivity2();
408 
409  m_DiffusionTensor.Fill(0.0);
410 
411  for (unsigned int i = 0;i < m_SpaceDimension;++i)
412  m_DiffusionTensor(i,i) = radialDiff2;
413 
414  for (unsigned int i = 0;i < m_SpaceDimension;++i)
415  for (unsigned int j = i;j < m_SpaceDimension;++j)
416  {
417  m_DiffusionTensor(i,j) += m_EigenVector1[i] * m_EigenVector1[j] * (axialDiff - radialDiff2);
418  m_DiffusionTensor(i,j) += m_EigenVector2[i] * m_EigenVector2[j] * (radialDiff1 - radialDiff2);
419 
420  if (i != j)
421  m_DiffusionTensor(j,i) = m_DiffusionTensor(i,j);
422  }
423 
424  m_TensorDeterminant = axialDiff * radialDiff1 * radialDiff2;
425 
426  m_UpdateInverseTensor = true;
427  m_ModifiedTensor = false;
428 }
429 
431 {
432  this->UpdateDiffusionTensor();
433 
434  if (!m_UpdateInverseTensor)
435  return;
436 
437  if (m_TensorDeterminant == 0.0)
438  {
439  m_InverseDiffusionTensor.Fill(0.0);
440  m_UpdateInverseTensor = false;
441  return;
442  }
443 
444  m_WorkVnlMatrix1.set_size(m_SpaceDimension,m_SpaceDimension);
445  m_WorkVnlMatrix1 = m_DiffusionTensor.GetVnlMatrix().as_matrix();
446  m_leCalculator->GetTensorPower(m_WorkVnlMatrix1, m_WorkVnlMatrix2, -1.0);
447  m_InverseDiffusionTensor = m_WorkVnlMatrix2;
448 
449  m_UpdateInverseTensor = false;
450 }
451 
453 {
454  if (!m_ModifiedAngles)
455  return;
456 
457  m_CosTheta = std::cos(this->GetOrientationTheta());
458  m_SinTheta = std::sin(this->GetOrientationTheta());
459  m_CosPhi = std::cos(this->GetOrientationPhi());
460  m_SinPhi = std::sin(this->GetOrientationPhi());
461  m_CosAlpha = std::cos(this->GetPerpendicularAngle());
462  m_SinAlpha = std::sin(this->GetPerpendicularAngle());
463 
464  m_EigenVector1[0] = m_SinTheta * m_CosPhi;
465  m_EigenVector1[1] = m_SinTheta * m_SinPhi;
466  m_EigenVector1[2] = m_CosTheta;
467 
468  m_EigenVector2[0] = -1.0 * (m_CosTheta * m_CosPhi * m_SinAlpha + m_SinPhi * m_CosAlpha);
469  m_EigenVector2[1] = m_CosPhi * m_CosAlpha - m_CosTheta * m_SinPhi * m_SinAlpha;
470  m_EigenVector2[2] = m_SinTheta * m_SinAlpha;
471 
472  m_ModifiedAngles = false;
473 }
474 
476 {
477  double l1 = this->GetAxialDiffusivity();
478  double l2 = this->GetRadialDiffusivity1();
479  double l3 = this->GetRadialDiffusivity2();
480 
481  double numFA = std::sqrt ((l1 - l2) * (l1 - l2) + (l2 - l3) * (l2 - l3) + (l3 - l1) * (l3 - l1));
482  double denomFA = std::sqrt (l1 * l1 + l2 * l2 + l3 * l3);
483 
484  double fa = 0;
485  if (denomFA != 0)
486  fa = std::sqrt(0.5) * (numFA / denomFA);
487 
488  return fa;
489 }
490 
492 {
493  double l1 = this->GetAxialDiffusivity();
494  double l2 = this->GetRadialDiffusivity1();
495  double l3 = this->GetRadialDiffusivity2();
496 
497  return (l1 + l2 + l3) / 3.0;
498 }
499 
501 {
502  return this->GetAxialDiffusivity();
503 }
504 
506 {
507  return (this->GetRadialDiffusivity1() + this->GetRadialDiffusivity2()) / 2.0;
508 }
509 
510 } //end namespace anima
void SetAxialDiffusivity(double num) ITK_OVERRIDE
const double MCMZeroLowerBound
Default zero lower bound (in case we want something else than zero one day)
std::vector< double > ListType
virtual void SetAxialDiffusivity(double num)
double GetPerpendicularAngle() ITK_OVERRIDE
ListType m_ParametersLowerBoundsVector
Vector holding current parameters lower bounds.
virtual double GetFourierTransformedDiffusionProfile(double smallDelta, double bigDelta, double gradientStrength, const Vector3DType &gradient) ITK_OVERRIDE
Superclass::ModelOutputVectorType ModelOutputVectorType
Superclass::Vector3DType Vector3DType
const double MCMAzimuthAngleUpperBound
Azimuth angle upper bound.
virtual ListType & GetSignalAttenuationJacobian(double smallDelta, double bigDelta, double gradientStrength, const Vector3DType &gradient) ITK_OVERRIDE
virtual void SetRadialDiffusivity2(double num)
void SetOrientationTheta(double num) ITK_OVERRIDE
virtual void SetParametersFromVector(const ListType &params) ITK_OVERRIDE
Various methods for optimization parameters setting and getting.
virtual double GetOrientationPhi()
virtual ListType & GetParameterUpperBounds() ITK_OVERRIDE
void UpdateDiffusionTensor()
Update diffusion tensor value from parameters.
virtual double GetLogDiffusionProfile(const Vector3DType &sample) ITK_OVERRIDE
virtual void SetRadialDiffusivity1(double num)
void SetRadialDiffusivity1(double num) ITK_OVERRIDE
double GetApparentParallelDiffusivity() ITK_OVERRIDE
Get approximation to parallel diffusivity of the compartment (may be different from axial diffusivity...
void Reorient(vnl_matrix< double > &orientationMatrix, bool affineTransform) ITK_OVERRIDE
Reorient the fascicle compartment using a matrix, the flag specifies if the transform is affine or ri...
void SetOrientationPhi(double num) ITK_OVERRIDE
virtual void SetPerpendicularAngle(double num)
virtual double GetAxialDiffusivity()
const double MCMRadialDiffusivityUpperBound
Radial diffusivity upper bound.
ListType m_ParametersVector
Vector holding current parameters vector.
virtual double GetPerpendicularAngle()
double GetBValueFromAcquisitionParameters(double smallDelta, double bigDelta, double gradientStrength)
Given gyromagnetic ratio in rad/s/T, gradient strength in T/mm and deltas in s, computes b-value in s...
virtual double GetOrientationTheta()
const double MCMPolarAngleUpperBound
Polar angle upper bound (used in tensor for now)
void GetVectorRepresentation(const vnl_matrix< T1 > &tensor, itk::VariableLengthVector< T2 > &vector, unsigned int vecDim=0, bool scale=false)
void SetPerpendicularAngle(double num) ITK_OVERRIDE
unsigned int GetCompartmentSize() ITK_OVERRIDE
Number of parameters to describe the model, these parameters will be self-contained, i.e. include fixed parameters for example.
void RotateSymmetricMatrix(T1 &tensor, T2 &rotationMatrix, T3 &rotated_tensor, unsigned int tensorDim)
Rotates a symmetric matrix by performing R^T D R where R is a rotation matrix and D the symmetric mat...
const Matrix3DType & GetDiffusionTensor() ITK_OVERRIDE
Get compartment as a 3D tensor (default behavior: throw exception if not supported by the compartment...
virtual void SetOrientationTheta(double num)
double GetOrientationTheta() ITK_OVERRIDE
void GetTensorFromVectorRepresentation(const itk::VariableLengthVector< T1 > &vector, vnl_matrix< T2 > &tensor, unsigned int tensDim=0, bool scale=false)
unsigned int GetNumberOfParameters() ITK_OVERRIDE
Number of optimized parameters: smaller than total number of parameters.
ListType m_ParametersUpperBoundsVector
Vector holding current parameters upper bounds.
ModelOutputVectorType & GetCompartmentVector() ITK_OVERRIDE
Get compartment overall description vector, mainly for writing, should be self-contained.
void UpdateInverseDiffusionTensor()
Update inverse diffusion tensor value from parameters.
double GetOrientationPhi() ITK_OVERRIDE
Superclass::Matrix3DType Matrix3DType
double GetApparentPerpendicularDiffusivity() ITK_OVERRIDE
Get approximation to perpendicular diffusivity of the compartment (may be different from radial diffu...
virtual double GetRadialDiffusivity2()
double GetApparentFractionalAnisotropy() ITK_OVERRIDE
ListType m_JacobianVector
Vector holding current jacobian value.
double GetRadialDiffusivity1() ITK_OVERRIDE
void SetRadialDiffusivity2(double num) ITK_OVERRIDE
const double MCMDiffusivityUpperBound
Diffusivity upper bound.
virtual double GetRadialDiffusivity1()
double GetAxialDiffusivity() ITK_OVERRIDE
double ComputeScalarProduct(const VectorType &v1, const VectorType &v2, const unsigned int NDimension)
virtual ListType & GetParametersAsVector() ITK_OVERRIDE
virtual void SetOrientationPhi(double num)
const double MCMDiffusivityLowerBound
Diffusivity lower bound for estimation.
void SetCompartmentVector(ModelOutputVectorType &compartmentVector) ITK_OVERRIDE
Set compartment overall description vector, for setting automatically the individual parameters when ...
void TransformCartesianToSphericalCoordinates(const VectorType &v, VectorType &resVec)
void UpdateAngleConfiguration()
Update angles&#39; sine and cosine values from parameters.
ModelOutputVectorType m_CompartmentVector
Vector to hold working value of compartment vector.
static const unsigned int m_SpaceDimension
void ExtractPPDRotationFromJacobianMatrix(vnl_matrix< RealType > &jacobianMatrix, vnl_matrix< RealType > &rotationMatrix, MatrixType &eigenVectors)
double GetApparentMeanDiffusivity() ITK_OVERRIDE
const double MCMAxialDiffusivityAddonLowerBound
Axial diffusivity add on to lower bound (used to ensure a minimal anisotropy to the anisotropic compa...
double GetRadialDiffusivity2() ITK_OVERRIDE
virtual ListType & GetParameterLowerBounds() ITK_OVERRIDE