ANIMA  4.0
animaApproximateMCMSmoothingCostFunction.cxx
Go to the documentation of this file.
2 
3 namespace anima
4 {
5 
6 void
8 ::SetReferenceModels(const std::vector <MCMPointer> &refModels, const std::vector <GradientType> &gradients,
9  const double &smallDelta, const double &bigDelta, const std::vector <double> &gradientStrengths)
10 {
11  unsigned int numRefModels = refModels.size();
12  m_ReferenceModelSignalValues.resize(numRefModels);
13  unsigned int numSamples = gradients.size();
14  std::vector <double> modelSignalValues(numSamples,0);
15 
16  for (unsigned int i = 0;i < numRefModels;++i)
17  {
18  for (unsigned int j = 0;j < numSamples;++j)
19  modelSignalValues[j] = refModels[i]->GetPredictedSignal(smallDelta, bigDelta,
20  gradientStrengths[j], gradients[j]);
21 
22  m_ReferenceModelSignalValues[i] = modelSignalValues;
23  }
24 
25  m_UpdatedData = true;
26 }
27 
28 void
30 ::SetMovingModels(const std::vector<MCMPointer> &movingModels, const std::vector <GradientType> &gradients,
31  const double &smallDelta, const double &bigDelta, const std::vector <double> &gradientStrengths)
32 {
33  unsigned int numMovingModels = movingModels.size();
34  m_MovingModelSignalValues.resize(numMovingModels);
35  unsigned int numSamples = gradients.size();
36  std::vector <double> modelSignalValues(numSamples,0);
37 
38  for (unsigned int i = 0;i < numMovingModels;++i)
39  {
40  for (unsigned int j = 0;j < numSamples;++j)
41  modelSignalValues[j] = movingModels[i]->GetPredictedSignal(smallDelta, bigDelta,
42  gradientStrengths[j], gradients[j]);
43 
44  m_MovingModelSignalValues[i] = modelSignalValues;
45  }
46 
47  m_UpdatedData = true;
48 }
49 
50 void
52 ::SetSmallDelta(double val)
53 {
54  m_SmallDelta = val;
55  m_UpdatedData = true;
56 }
57 
58 void
60 ::SetBigDelta(double val)
61 {
62  m_BigDelta = val;
63  m_UpdatedData = true;
64 }
65 
66 void
68 ::SetGradientStrengths(const std::vector <double> &val)
69 {
70  m_GradientStrengths = val;
71  m_UpdatedData = true;
72 }
73 
74 void
76 ::SetGradientDirections(const std::vector <GradientType> &val)
77 {
78  m_GradientDirections = val;
79  m_UpdatedData = true;
80 }
81 
82 void
84 ::SetBValueWeightIndexes(const std::vector <unsigned int> &val)
85 {
86  m_BValueWeightIndexes = val;
87  m_UpdatedData = true;
88 }
89 
90 void
92 ::SetSphereWeights(const std::vector <double> &val)
93 {
94  m_SphereWeights = val;
95  m_UpdatedData = true;
96 }
97 
100 ::GetValue(const ParametersType &parameters) const
101 {
102  unsigned int numDataPoints = m_ReferenceModelSignalValues.size();
103  if (numDataPoints == 0)
104  return 0.0;
105 
106  double gaussianSigma = parameters[0] * m_ParameterScale;
107  MeasureType outputValue = 0;
108 
109  if (m_UpdatedData)
110  {
111  m_ConstantTerm = 0;
112  for (unsigned int j = 0;j < m_GradientDirections.size();++j)
113  {
114  double constantAddonValue = 0;
115  for (unsigned int l = 0;l < numDataPoints;++l)
116  constantAddonValue += m_MovingModelSignalValues[l][j] * m_MovingModelSignalValues[l][j];
117 
118  m_ConstantTerm += m_SphereWeights[m_BValueWeightIndexes[j]] * constantAddonValue;
119  }
120  }
121 
122  outputValue = m_ConstantTerm;
123 
124  for (unsigned int j = 0;j < m_GradientDirections.size();++j)
125  {
126  double bValue = anima::GetBValueFromAcquisitionParameters(m_SmallDelta, m_BigDelta, m_GradientStrengths[j]);
127  double gaussianValue = 0;
128  for (unsigned int i = 0;i < m_GradientDirections[j].size();++i)
129  gaussianValue += m_GradientDirections[j][i] * m_GradientDirections[j][i];
130 
131  gaussianValue = std::exp (- bValue * gaussianSigma * gaussianValue);
132 
133  double addonValue = 0;
134  for (unsigned int l = 0;l < numDataPoints;++l)
135  {
136  addonValue -= 2.0 * m_MovingModelSignalValues[l][j] * m_ReferenceModelSignalValues[l][j];
137  addonValue += gaussianValue * m_ReferenceModelSignalValues[l][j] * m_ReferenceModelSignalValues[l][j];
138  }
139 
140  outputValue += m_SphereWeights[m_BValueWeightIndexes[j]] * gaussianValue * addonValue;
141  }
142 
143  if (outputValue <= 0)
144  outputValue = 0;
145 
146  return outputValue / numDataPoints;
147 }
148 
149 void
151 ::GetDerivative(const ParametersType &parameters, DerivativeType &derivative) const
152 {
153  unsigned int numDataPoints = m_ReferenceModelSignalValues.size();
154  derivative.set_size(this->GetNumberOfParameters());
155  derivative.fill(0);
156 
157  if (numDataPoints == 0)
158  return;
159 
160  double gaussianSigma = parameters[0] * m_ParameterScale;
161 
162  for (unsigned int j = 0;j < m_GradientDirections.size();++j)
163  {
164  double bValue = anima::GetBValueFromAcquisitionParameters(m_SmallDelta, m_BigDelta, m_GradientStrengths[j]);
165  double addonValue = 0;
166  double gaussianDotProduct = 0;
167  for (unsigned int i = 0;i < m_GradientDirections[j].size();++i)
168  gaussianDotProduct += m_GradientDirections[j][i] * m_GradientDirections[j][i];
169  gaussianDotProduct *= bValue;
170 
171  double gaussianDerivativeValue = - gaussianDotProduct * std::exp (- gaussianSigma * gaussianDotProduct);
172  double gaussianSqDerivativeValue = - 2.0 * gaussianDotProduct * std::exp (- 2.0 * gaussianSigma * gaussianDotProduct);
173 
174  for (unsigned int l = 0;l < numDataPoints;++l)
175  {
176  double firstTerm = gaussianSqDerivativeValue * m_ReferenceModelSignalValues[l][j] * m_ReferenceModelSignalValues[l][j];
177  double secondTerm = - 2.0 * gaussianDerivativeValue * m_MovingModelSignalValues[l][j] - m_ReferenceModelSignalValues[l][j];
178 
179  addonValue += firstTerm + secondTerm;
180  }
181 
182  derivative[0] += m_SphereWeights[m_BValueWeightIndexes[j]] * addonValue;
183  }
184 
185  derivative[0] /= numDataPoints;
186 }
187 
188 } // end namespace anima
virtual MeasureType GetValue(const ParametersType &parameters) const ITK_OVERRIDE
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...
void SetMovingModels(const std::vector< MCMPointer > &movingModels, const std::vector< GradientType > &gradients, const double &smallDelta, const double &bigDelta, const std::vector< double > &gradientStrengths)
void SetBValueWeightIndexes(const std::vector< unsigned int > &val)
void SetGradientDirections(const std::vector< GradientType > &val)
virtual void GetDerivative(const ParametersType &parameters, DerivativeType &derivative) const ITK_OVERRIDE
void SetReferenceModels(const std::vector< MCMPointer > &refModels, const std::vector< GradientType > &gradients, const double &smallDelta, const double &bigDelta, const std::vector< double > &gradientStrengths)