2 #include <vnl/algo/vnl_determinant.h> 11 unsigned int numRefModels = refModels.size();
12 m_ReferenceModels.resize(numRefModels);
13 m_ReferenceModelWeights.resize(numRefModels);
14 m_ReferenceNumberOfIsotropicCompartments.resize(numRefModels);
16 for (
unsigned int i = 0;i < numRefModels;++i)
18 unsigned int numCompartments = refModels[i]->GetNumberOfCompartments();
19 std::vector <TensorType> compartmentTensors(numCompartments);
20 std::vector <double> compartmentWeights(numCompartments);
22 unsigned int numIsoCompartments = refModels[i]->GetNumberOfIsotropicCompartments();
23 for (
unsigned int j = 0;j < numCompartments;++j)
25 double currentWeight = refModels[i]->GetCompartmentWeight(j);
26 if (currentWeight > 0)
28 compartmentTensors[pos] = refModels[i]->GetCompartment(j)->GetDiffusionTensor();
29 compartmentWeights[pos] = currentWeight;
32 else if (j < refModels[i]->GetNumberOfIsotropicCompartments())
36 compartmentTensors.resize(pos);
37 compartmentWeights.resize(pos);
39 m_ReferenceModels[i] = compartmentTensors;
40 m_ReferenceModelWeights[i] = compartmentWeights;
41 m_ReferenceNumberOfIsotropicCompartments[i] = 0;
44 m_UpdatedReferenceData =
true;
45 m_RecomputeConstantTerm =
true;
52 unsigned int numMovingModels = movingModels.size();
53 m_MovingModels.resize(numMovingModels);
54 m_MovingModelWeights.resize(numMovingModels);
56 for (
unsigned int i = 0;i < numMovingModels;++i)
58 unsigned int numCompartments = movingModels[i]->GetNumberOfCompartments();
59 std::vector <TensorType> compartmentTensors(numCompartments);
60 std::vector <double> compartmentWeights(numCompartments);
62 for (
unsigned int j = 0;j < numCompartments;++j)
64 double currentWeight = movingModels[i]->GetCompartmentWeight(j);
65 if (currentWeight > 0)
67 compartmentTensors[pos] = movingModels[i]->GetCompartment(j)->GetDiffusionTensor();
68 compartmentWeights[pos] = currentWeight;
73 compartmentTensors.resize(pos);
74 compartmentWeights.resize(pos);
76 m_MovingModels[i] = compartmentTensors;
77 m_MovingModelWeights[i] = compartmentWeights;
80 m_UpdatedMovingData =
true;
81 m_RecomputeConstantTerm =
true;
88 unsigned int numDataPoints = m_ReferenceModels.size();
89 if (numDataPoints == 0)
92 this->UpdateDeterminants();
94 double gaussianSigma = parameters[0] / m_TensorsScale;
97 double pi3half = std::pow(M_PI,1.5);
98 double pi23half = std::pow(2.0 * M_PI,1.5);
100 if (m_RecomputeConstantTerm)
103 for (
unsigned int l = 0;l < numDataPoints;++l)
105 unsigned int pos = 0;
106 unsigned int posRefMov = 0;
107 unsigned int refModelSize = m_ReferenceModelWeights[l].size();
108 unsigned int movingModelSize = m_MovingModelWeights[l].size();
109 unsigned int numRefIsoCompartments = m_ReferenceNumberOfIsotropicCompartments[l];
111 for (
unsigned int i = 0;i < numRefIsoCompartments;++i)
113 m_ConstantTerm += m_ReferenceModelWeights[l][i] * m_ReferenceModelWeights[l][i] * pi3half / std::sqrt(m_ReferenceReferenceDeterminants[l][pos]);
116 for (
unsigned int j = i+1;j < numRefIsoCompartments;++j)
118 m_ConstantTerm += 2.0 * m_ReferenceModelWeights[l][i] * m_ReferenceModelWeights[l][j] * pi23half / std::sqrt(m_ReferenceReferenceDeterminants[l][pos]);
122 pos += refModelSize - numRefIsoCompartments;
124 for (
unsigned int j = 0;j < movingModelSize;++j)
126 m_ConstantTerm -= 2.0 * m_ReferenceModelWeights[l][i] * m_MovingModelWeights[l][j] * pi23half / std::sqrt(m_ReferenceMovingDeterminants[l][posRefMov]);
132 for (
unsigned int i = 0;i < movingModelSize;++i)
134 m_ConstantTerm += m_MovingModelWeights[l][i] * m_MovingModelWeights[l][i] * pi3half / std::sqrt(m_MovingMovingDeterminants[l][pos]);
137 for (
unsigned int j = i+1;j < movingModelSize;++j)
139 m_ConstantTerm += 2.0 * m_MovingModelWeights[l][i] * m_MovingModelWeights[l][j] * pi23half / std::sqrt(m_MovingMovingDeterminants[l][pos]);
145 m_RecomputeConstantTerm =
false;
148 outputValue = m_ConstantTerm;
150 for (
unsigned int l = 0;l < numDataPoints;++l)
152 unsigned int refModelSize = m_ReferenceModelWeights[l].size();
153 unsigned int movingModelSize = m_MovingModelWeights[l].size();
154 unsigned int numRefIsoCompartments = m_ReferenceNumberOfIsotropicCompartments[l];
156 unsigned int posRefRef = 0;
157 unsigned int posRefMov = 0;
159 for (
unsigned int i = 0;i < refModelSize;++i)
161 if (i >= numRefIsoCompartments)
163 double refSigmaDeterminant = m_ReferenceReferenceDeterminants[l][posRefRef] + gaussianSigma * (m_ReferenceReferenceTraceDifferences[l][posRefRef] + gaussianSigma * (m_ReferenceReferenceTraces[l][posRefRef] + gaussianSigma));
164 outputValue += m_ReferenceModelWeights[l][i] * m_ReferenceModelWeights[l][i] * pi3half / std::sqrt(refSigmaDeterminant);
166 for (
unsigned int j = 0;j < movingModelSize;++j)
168 double refMovSigmaDeterminant = m_ReferenceMovingDeterminants[l][posRefMov] + gaussianSigma * (m_ReferenceMovingTraceDifferences[l][posRefMov] + gaussianSigma * (m_ReferenceMovingTraces[l][posRefMov] + gaussianSigma));
170 outputValue -= 2.0 * m_ReferenceModelWeights[l][i] * m_MovingModelWeights[l][j] * pi23half / std::sqrt(refMovSigmaDeterminant);
175 posRefMov += movingModelSize;
179 for (
unsigned int j = i+1;j < refModelSize;++j)
181 double factor = (i >= numRefIsoCompartments) + (j >= numRefIsoCompartments);
188 double alphaValue = factor * gaussianSigma;
189 double refSigmaDeterminant = m_ReferenceReferenceDeterminants[l][posRefRef] + alphaValue * (m_ReferenceReferenceTraceDifferences[l][posRefRef] + alphaValue * (m_ReferenceReferenceTraces[l][posRefRef] + alphaValue));
191 outputValue += 2.0 * m_ReferenceModelWeights[l][i] * m_ReferenceModelWeights[l][j] * pi23half / std::sqrt(refSigmaDeterminant);
197 if (outputValue <= 0)
200 return outputValue / numDataPoints;
208 unsigned int numDataPoints = m_ReferenceModels.size();
209 double squaredMatrixTrace;
211 if (m_UpdatedReferenceData)
213 m_ReferenceReferenceDeterminants.resize(numDataPoints);
214 m_ReferenceReferenceTraces.resize(numDataPoints);
215 m_ReferenceReferenceTraceDifferences.resize(numDataPoints);
217 for (
unsigned int l = 0;l < numDataPoints;++l)
219 unsigned int refModelSize = m_ReferenceModelWeights[l].size();
220 unsigned int numElements = refModelSize * (refModelSize + 1) / 2;
221 m_ReferenceReferenceDeterminants[l].resize(numElements);
222 m_ReferenceReferenceTraces[l].resize(numElements);
223 m_ReferenceReferenceTraceDifferences[l].resize(numElements);
225 unsigned int pos = 0;
226 for (
unsigned int i = 0;i < refModelSize;++i)
228 workMatrix = m_ReferenceModels[l][i];
229 for (
unsigned int k = 0;k < 3;++k)
230 workMatrix(k,k) += 1.0 / (2.0 * m_LowPassGaussianSigma);
231 m_ReferenceReferenceDeterminants[l][pos] = vnl_determinant <double> (workMatrix.GetVnlMatrix().as_matrix());
232 this->ComputeMatrixTraces(m_ReferenceModels[l][i],m_ReferenceReferenceTraces[l][pos],squaredMatrixTrace);
233 m_ReferenceReferenceTraceDifferences[l][pos] = 0.5 * (m_ReferenceReferenceTraces[l][pos] * m_ReferenceReferenceTraces[l][pos] - squaredMatrixTrace);
237 for (
unsigned int j = i+1;j < refModelSize;++j)
239 workMatrix = m_ReferenceModels[l][i] + m_ReferenceModels[l][j];
240 for (
unsigned int k = 0;k < 3;++k)
241 workMatrix(k,k) += 1.0 / m_LowPassGaussianSigma;
242 m_ReferenceReferenceDeterminants[l][pos] = vnl_determinant <double> (workMatrix.GetVnlMatrix().as_matrix());
243 this->ComputeMatrixTraces(workMatrix,m_ReferenceReferenceTraces[l][pos],squaredMatrixTrace);
244 m_ReferenceReferenceTraceDifferences[l][pos] = 0.5 * (m_ReferenceReferenceTraces[l][pos] * m_ReferenceReferenceTraces[l][pos] - squaredMatrixTrace);
252 if (m_UpdatedMovingData)
254 m_MovingMovingDeterminants.resize(numDataPoints);
256 for (
unsigned int l = 0;l < numDataPoints;++l)
258 unsigned int movingModelSize = m_MovingModelWeights[l].size();
259 unsigned int numElements = movingModelSize * (movingModelSize + 1) / 2;
260 m_MovingMovingDeterminants[l].resize(numElements);
262 unsigned int pos = 0;
263 for (
unsigned int i = 0;i < movingModelSize;++i)
265 workMatrix = m_MovingModels[l][i];
266 for (
unsigned int k = 0;k < 3;++k)
267 workMatrix(k,k) += 1.0 / (2.0 * m_LowPassGaussianSigma);
268 m_MovingMovingDeterminants[l][pos] = vnl_determinant <double> (workMatrix.GetVnlMatrix().as_matrix());
271 for (
unsigned int j = i+1;j < movingModelSize;++j)
273 workMatrix = m_MovingModels[l][i] + m_MovingModels[l][j];
274 for (
unsigned int k = 0;k < 3;++k)
275 workMatrix(k,k) += 1.0 / m_LowPassGaussianSigma;
276 m_MovingMovingDeterminants[l][pos] = vnl_determinant <double> (workMatrix.GetVnlMatrix().as_matrix());
283 if (m_UpdatedMovingData || m_UpdatedReferenceData)
285 m_ReferenceMovingDeterminants.resize(numDataPoints);
286 m_ReferenceMovingTraces.resize(numDataPoints);
287 m_ReferenceMovingTraceDifferences.resize(numDataPoints);
289 for (
unsigned int l = 0;l < numDataPoints;++l)
291 unsigned int refModelSize = m_ReferenceModelWeights[l].size();
292 unsigned int movingModelSize = m_MovingModelWeights[l].size();
293 unsigned int numElements = refModelSize * movingModelSize;
294 m_ReferenceMovingDeterminants[l].resize(numElements);
295 m_ReferenceMovingTraces[l].resize(numElements);
296 m_ReferenceMovingTraceDifferences[l].resize(numElements);
298 unsigned int pos = 0;
299 for (
unsigned int i = 0;i < refModelSize;++i)
301 for (
unsigned int j = 0;j < movingModelSize;++j)
303 workMatrix = m_ReferenceModels[l][i] + m_MovingModels[l][j];
304 for (
unsigned int k = 0;k < 3;++k)
305 workMatrix(k,k) += 1.0 / m_LowPassGaussianSigma;
306 m_ReferenceMovingDeterminants[l][pos] = vnl_determinant <double> (workMatrix.GetVnlMatrix().as_matrix());
307 this->ComputeMatrixTraces(workMatrix,m_ReferenceMovingTraces[l][pos],squaredMatrixTrace);
308 m_ReferenceMovingTraceDifferences[l][pos] = 0.5 * (m_ReferenceMovingTraces[l][pos] * m_ReferenceMovingTraces[l][pos] - squaredMatrixTrace);
316 m_UpdatedReferenceData =
false;
317 m_UpdatedMovingData =
false;
325 squaredMatrixTrace = 0;
327 for (
unsigned int j = 0;j < TensorType::RowDimensions;++j)
329 matrixTrace += matrix(j,j);
330 for (
unsigned int k = j;k < TensorType::ColumnDimensions;++k)
336 squaredMatrixTrace += factor * matrix(j,k) * matrix(j,k);
345 unsigned int numDataPoints = m_ReferenceModels.size();
346 derivative.set_size(this->GetNumberOfParameters());
349 if (numDataPoints == 0)
352 this->UpdateDeterminants();
354 double gaussianSigma = parameters[0] / m_TensorsScale;
355 double sq2pi3half = std::sqrt(2.0) * std::pow(M_PI,1.5);
356 double spi3half = std::pow(M_PI,1.5) / 16;
358 for (
unsigned int l = 0;l < numDataPoints;++l)
360 unsigned int refModelSize = m_ReferenceModelWeights[l].size();
361 unsigned int movingModelSize = m_MovingModelWeights[l].size();
362 unsigned int numRefIsoCompartments = m_ReferenceNumberOfIsotropicCompartments[l];
364 unsigned int posRefRef = 0;
365 unsigned int posRefMov = 0;
367 for (
unsigned int i = 0;i < refModelSize;++i)
369 if (i >= numRefIsoCompartments)
371 double refSigmaDeterminant = m_ReferenceReferenceDeterminants[l][posRefRef] + gaussianSigma * (m_ReferenceReferenceTraceDifferences[l][posRefRef] + gaussianSigma * (m_ReferenceReferenceTraces[l][posRefRef] + gaussianSigma));
373 double traceDiff = 2.0 * m_ReferenceReferenceTraceDifferences[l][posRefRef] + 4.0 * gaussianSigma * m_ReferenceReferenceTraces[l][posRefRef] + 6.0 * gaussianSigma * gaussianSigma;
375 double matrixDetCube = refSigmaDeterminant * refSigmaDeterminant * refSigmaDeterminant;
377 derivative[0] -= m_ReferenceModelWeights[l][i] * m_ReferenceModelWeights[l][i] * spi3half * traceDiff / std::sqrt(matrixDetCube);
379 for (
unsigned int j = 0;j < movingModelSize;++j)
381 double refMovSigmaDeterminant = m_ReferenceMovingDeterminants[l][posRefMov] + gaussianSigma * (m_ReferenceMovingTraceDifferences[l][posRefMov] + gaussianSigma * (m_ReferenceMovingTraces[l][posRefMov] + gaussianSigma));
383 traceDiff = 2.0 * m_ReferenceMovingTraceDifferences[l][posRefMov] + 4.0 * gaussianSigma * m_ReferenceMovingTraces[l][posRefMov] + 6.0 * gaussianSigma * gaussianSigma;
384 matrixDetCube = refMovSigmaDeterminant * refMovSigmaDeterminant * refMovSigmaDeterminant;
386 derivative[0] += m_ReferenceModelWeights[l][i] * m_MovingModelWeights[l][j] * sq2pi3half * traceDiff / std::sqrt(matrixDetCube);
391 posRefMov += movingModelSize;
395 for (
unsigned int j = i+1;j < refModelSize;++j)
397 double factor = (i >= numRefIsoCompartments) + (j >= numRefIsoCompartments);
404 double alphaValue = factor * gaussianSigma;
405 double refSigmaDeterminant = m_ReferenceReferenceDeterminants[l][posRefRef] + alphaValue * (m_ReferenceReferenceTraceDifferences[l][posRefRef] + alphaValue * (m_ReferenceReferenceTraces[l][posRefRef] + alphaValue));
407 double traceDiff = 2.0 * m_ReferenceReferenceTraceDifferences[l][posRefRef] + 4.0 * alphaValue * m_ReferenceReferenceTraces[l][posRefRef] + 6.0 * alphaValue * alphaValue;
408 double matrixDetCube = refSigmaDeterminant * refSigmaDeterminant * refSigmaDeterminant;
410 derivative[0] -= factor * m_ReferenceModelWeights[l][i] * m_ReferenceModelWeights[l][j] * sq2pi3half * traceDiff / std::sqrt(matrixDetCube);
416 derivative[0] /= numDataPoints;
void SetMovingModels(const std::vector< MCMPointer > &movingModels)
Superclass::ParametersType ParametersType
void SetReferenceModels(const std::vector< MCMPointer > &refModels)
Superclass::MeasureType MeasureType
Superclass::DerivativeType DerivativeType
virtual void GetDerivative(const ParametersType ¶meters, DerivativeType &derivative) const ITK_OVERRIDE
anima::BaseCompartment::Matrix3DType TensorType
virtual MeasureType GetValue(const ParametersType ¶meters) const ITK_OVERRIDE
void ComputeMatrixTraces(const TensorType &matrix, double &matrixTrace, double &squaredMatrixTrace) const
void UpdateDeterminants() const