ANIMA  4.0
animaB1GMMRelaxometryCostFunction.cxx
Go to the documentation of this file.
3 #include <animaBaseTensorTools.h>
5 
6 #include <boost/math/special_functions/erf.hpp>
7 
8 namespace anima
9 {
10 
13 {
14  m_TestedParameters.SetSize(this->GetNumberOfParameters());
15 
16  for (unsigned int i = 0;i < this->GetNumberOfParameters();++i)
17  m_TestedParameters[i] = parameters[i];
18 
19  this->PrepareDataForLLS();
20 
22 
23  return m_SigmaSquare;
24 }
25 
26 void
28  DerivativeType &derivative) const
29 {
30  itkExceptionMacro("Derivative not handled here. Too much work worth nothing.")
31 }
32 
33 void
35 {
36  unsigned int numT2Signals = m_T2RelaxometrySignals.size();
37  unsigned int numDistributions = m_GaussianMeans.size();
38 
39  m_T2SignalSimulator.SetNumberOfEchoes(numT2Signals);
40  m_T2SignalSimulator.SetEchoSpacing(m_EchoSpacing);
41  m_T2SignalSimulator.SetExcitationFlipAngle(m_ExcitationFlipAngle);
42 
43  // Contains int_{space of ith distribution} EPG(t2, b1, jth echo) G(t2, mu_i, sigma_i) d t2
44  m_PredictedSignalAttenuations.set_size(numT2Signals,numDistributions);
45 
47  B1GMMDistributionIntegrand t2Integrand(m_T2SignalSimulator,epgVectors);
48  t2Integrand.SetT1Value(m_T1Value);
49  t2Integrand.SetFlipAngle(m_TestedParameters[0]);
50 
52 
53  for (unsigned int i = 0;i < numDistributions;++i)
54  {
55  t2Integrand.SetGaussianMean(m_GaussianMeans[i]);
56  t2Integrand.SetGaussianVariance(m_GaussianVariances[i]);
57  epgVectors.clear();
58 
59  double minInterestZoneValue = std::max(0.0, m_GaussianMeans[i] - 5.0 * std::sqrt(m_GaussianVariances[i]));
60  double theoreticalMinValue = m_GaussianMeans[i] - i * std::sqrt(m_GaussianVariances[i]);
61  double maxInterestZoneValue = m_GaussianMeans[i] + 5.0 * std::sqrt(m_GaussianVariances[i]);
62  if (theoreticalMinValue < 0.0)
63  maxInterestZoneValue -= theoreticalMinValue;
64 
65  glQuad.SetInterestZone(minInterestZoneValue, maxInterestZoneValue);
66 
67  double truncationAbscissa = 2.0 * m_GaussianMeans[i];
68  double truncatedGaussianIntegral = 0.5 * (1.0 + boost::math::erf((truncationAbscissa - m_GaussianMeans[i]) / (std::sqrt(2.0 * m_GaussianVariances[i]))));
69 
70  for (unsigned int j = 0;j < numT2Signals;++j)
71  {
72  t2Integrand.SetEchoNumber(j);
73  m_PredictedSignalAttenuations(j,i) = glQuad.GetIntegralValue(t2Integrand) / truncatedGaussianIntegral;
74  }
75  }
76 
77  m_CholeskyMatrix.set_size(numDistributions,numDistributions);
78  m_FSignals.SetSize(numDistributions);
79 
80  for (unsigned int i = 0;i < numDistributions;++i)
81  {
82  m_FSignals[i] = 0.0;
83  for (unsigned int j = 0;j < numT2Signals;++j)
84  m_FSignals[i] += m_PredictedSignalAttenuations(j,i) * m_T2RelaxometrySignals[j];
85 
86  for (unsigned int j = i;j < numDistributions;++j)
87  {
88  m_CholeskyMatrix(i,j) = 0;
89  for (unsigned int k = 0;k < numT2Signals;++k)
90  m_CholeskyMatrix(i,j) += m_PredictedSignalAttenuations(k,i) * m_PredictedSignalAttenuations(k,j);
91 
92  if (i != j)
93  m_CholeskyMatrix(j,i) = m_CholeskyMatrix(i,j);
94  }
95  }
96 
97  m_CholeskySolver.SetInputMatrix(m_CholeskyMatrix);
98  m_CholeskySolver.PerformDecomposition();
99  m_OptimalT2Weights = m_CholeskySolver.SolveLinearSystem(m_FSignals);
100 
101  bool nnlsNeeded = false;
102 
103  for (unsigned int i = 0;i < numDistributions;++i)
104  {
105  if (m_OptimalT2Weights[i] <= 0.0)
106  {
107  nnlsNeeded = true;
108  break;
109  }
110  }
111 
112  if (nnlsNeeded)
113  {
114  m_NNLSBordersOptimizer->SetDataMatrix(m_CholeskyMatrix);
115  m_NNLSBordersOptimizer->SetPoints(m_FSignals);
116  m_NNLSBordersOptimizer->SetSquaredProblem(true);
117  m_NNLSBordersOptimizer->StartOptimization();
118 
119  m_OptimalT2Weights = m_NNLSBordersOptimizer->GetCurrentPosition();
120  }
121 
122  m_Residuals.set_size(numT2Signals);
123 }
124 
125 void
127 {
128  unsigned int numT2Signals = m_T2RelaxometrySignals.size();
129  unsigned int numDistributions = m_GaussianMeans.size();
130 
131  m_SigmaSquare = 0.0;
132 
133  for (unsigned int i = 0;i < numT2Signals;++i)
134  {
135  m_Residuals[i] = m_T2RelaxometrySignals[i];
136  for (unsigned int j = 0;j < numDistributions;++j)
137  m_Residuals[i] -= m_PredictedSignalAttenuations(i,j) * m_OptimalT2Weights[j];
138 
139  m_SigmaSquare += m_Residuals[i] * m_Residuals[i];
140  }
141 
142  m_SigmaSquare /= numT2Signals;
143 }
144 
145 } // end namespace anima
double GetIntegralValue(FunctionType integrand)
Computes the Gauss Laguerre quadrature of a function defined from R^+ into R. Recenters the function ...
void SolveLinearLeastSquares() const
Computes maximum likelihood estimates of weights.
unsigned int GetNumberOfParameters() const ITK_OVERRIDE
void PerformDecomposition()
Performs LDL decomposition from input matrix.
void SetInputMatrix(const MatrixType &matrix)
Set input matrix to decompose.
virtual MeasureType GetValue(const ParametersType &parameters) const ITK_OVERRIDE
virtual void GetDerivative(const ParametersType &parameters, DerivativeType &derivative) const ITK_OVERRIDE
std::map< double, anima::EPGSignalSimulator::RealVectorType > EPGVectorsMapType
VectorType & SolveLinearSystem(const VectorType &b)
Solves linear system Ax=b and outputs result in new variable.
Integrand to compute the internal integral per distribution in B1GMMRelaxometryCostFunction.
void SetInterestZone(double minVal, double maxVal)
Specifies region on which the main part of the function is to be seen. If not specified, R^+ is the region.
void SetNumberOfEchoes(unsigned int val)