ANIMA  4.0
animaDTIEstimationImageFilter.hxx
Go to the documentation of this file.
1 #pragma once
3 
4 #include <itkImageRegionIterator.h>
5 #include <itkImageRegionIteratorWithIndex.h>
6 #include <itkImageRegionConstIterator.h>
7 
8 #include <vnl/algo/vnl_determinant.h>
9 #include <vnl/algo/vnl_matrix_inverse.h>
10 #include <vnl/algo/vnl_symmetric_eigensystem.h>
11 
12 #include <nlopt.hpp>
13 
14 #include <animaVectorOperations.h>
15 #include <animaBaseTensorTools.h>
16 
17 namespace anima
18 {
19 
20 template <class InputPixelScalarType, class OutputPixelScalarType>
21 void
23 ::AddGradientDirection(unsigned int i, vnl_vector_fixed<double,3> &grad)
24 {
25  if (i == m_GradientDirections.size())
26  m_GradientDirections.push_back(grad);
27  else if (i > m_GradientDirections.size())
28  std::cerr << "Trying to add a direction not contiguous... Add directions contiguously (0,1,2,3,...)..." << std::endl;
29  else
30  m_GradientDirections[i] = grad;
31 }
32 
33 template <class InputPixelScalarType, class OutputPixelScalarType>
34 void
37 {
38  if (this->GetComputationMask())
39  return;
40 
41  typedef itk::ImageRegionConstIterator <InputImageType> B0IteratorType;
42  typedef itk::ImageRegionIterator <MaskImageType> MaskIteratorType;
43 
44  unsigned int firstB0Index = 0;
45  while (m_BValuesList[firstB0Index] > 10)
46  ++firstB0Index;
47 
48  B0IteratorType b0Itr(this->GetInput(firstB0Index),this->GetOutput()->GetLargestPossibleRegion());
49 
50  typename MaskImageType::Pointer maskImage = MaskImageType::New();
51  maskImage->Initialize();
52  maskImage->SetRegions(this->GetInput(0)->GetLargestPossibleRegion());
53  maskImage->SetSpacing (this->GetInput(0)->GetSpacing());
54  maskImage->SetOrigin (this->GetInput(0)->GetOrigin());
55  maskImage->SetDirection (this->GetInput(0)->GetDirection());
56  maskImage->Allocate();
57 
58  MaskIteratorType maskItr (maskImage,this->GetOutput()->GetLargestPossibleRegion());
59  while (!b0Itr.IsAtEnd())
60  {
61  if (b0Itr.Get() > m_B0Threshold)
62  maskItr.Set(1);
63  else
64  maskItr.Set(0);
65 
66  ++b0Itr;
67  ++maskItr;
68  }
69 
70  this->SetComputationMask(maskImage);
71 }
72 
73 template <class InputPixelScalarType, class OutputPixelScalarType>
74 void
77 {
78  // Override the method in itkImageSource, so we can set the vector length of
79  // the output itk::VectorImage
80 
81  this->Superclass::GenerateOutputInformation();
82 
83  OutputImageType *output = this->GetOutput();
84  output->SetVectorLength(m_NumberOfComponents);
85 }
86 
87 template <class InputPixelScalarType, class OutputPixelScalarType>
88 void
91 {
92  if (m_BValuesList.size() != this->GetNumberOfIndexedInputs())
93  {
94  std::string error("There should be the same number of input images and input b-values... ");
95  error += std::to_string (m_BValuesList.size());
96  error += " ";
97  error += std::to_string (this->GetNumberOfIndexedInputs());
98 
99  throw itk::ExceptionObject(__FILE__, __LINE__,error,ITK_LOCATION);
100  }
101 
102  itk::ImageRegionIterator <OutputImageType> fillOut(this->GetOutput(),this->GetOutput()->GetLargestPossibleRegion());
103  typename OutputImageType::PixelType emptyModelVec(m_NumberOfComponents);
104  emptyModelVec.Fill(0);
105 
106  while (!fillOut.IsAtEnd())
107  {
108  fillOut.Set(emptyModelVec);
109  ++fillOut;
110  }
111 
112  m_EstimatedB0Image = OutputB0ImageType::New();
113  m_EstimatedB0Image->Initialize();
114  m_EstimatedB0Image->SetOrigin(this->GetOutput()->GetOrigin());
115  m_EstimatedB0Image->SetSpacing(this->GetOutput()->GetSpacing());
116  m_EstimatedB0Image->SetDirection(this->GetOutput()->GetDirection());
117  m_EstimatedB0Image->SetRegions(this->GetOutput()->GetLargestPossibleRegion());
118 
119  m_EstimatedB0Image->Allocate();
120  m_EstimatedB0Image->FillBuffer(0.0);
121 
122  m_EstimatedVarianceImage = OutputB0ImageType::New();
123  m_EstimatedVarianceImage->Initialize();
124  m_EstimatedVarianceImage->SetOrigin(this->GetOutput()->GetOrigin());
125  m_EstimatedVarianceImage->SetSpacing(this->GetOutput()->GetSpacing());
126  m_EstimatedVarianceImage->SetDirection(this->GetOutput()->GetDirection());
127  m_EstimatedVarianceImage->SetRegions(this->GetOutput()->GetLargestPossibleRegion());
128 
129  m_EstimatedVarianceImage->Allocate();
130  m_EstimatedVarianceImage->FillBuffer(0.0);
131 
132  vnl_matrix <double> initSolverSystem(this->GetNumberOfIndexedInputs(),m_NumberOfComponents + 1);
133  initSolverSystem.fill(0.0);
134  for (unsigned int i = 0;i < this->GetNumberOfIndexedInputs();++i)
135  {
136  initSolverSystem(i,0) = 1.0;
137  unsigned int pos = 1;
138  for (unsigned int j = 0;j < 3;++j)
139  {
140  for (unsigned int k = 0;k < j;++k)
141  {
142  initSolverSystem(i,pos) = - 2.0 * m_BValuesList[i] * m_GradientDirections[i][j] * m_GradientDirections[i][k];
143  ++pos;
144  }
145 
146  initSolverSystem(i,pos) = - m_BValuesList[i] * m_GradientDirections[i][j] * m_GradientDirections[i][j];
147  ++pos;
148  }
149  }
150 
151  vnl_matrix_inverse <double> inverter (initSolverSystem);
152  m_InitialMatrixSolver = inverter.pinverse();
153 
154  Superclass::BeforeThreadedGenerateData();
155 }
156 
157 template <class InputPixelScalarType, class OutputPixelScalarType>
158 void
161 {
162  typedef itk::ImageRegionConstIterator <InputImageType> ImageIteratorType;
163 
164  unsigned int numInputs = this->GetNumberOfIndexedInputs();
165  std::vector <ImageIteratorType> inIterators;
166  for (unsigned int i = 0;i < numInputs;++i)
167  inIterators.push_back(ImageIteratorType(this->GetInput(i),outputRegionForThread));
168 
169  typedef itk::ImageRegionConstIterator <MaskImageType> MaskIteratorType;
170  MaskIteratorType maskIterator(this->GetComputationMask(),outputRegionForThread);
171 
172  typedef itk::ImageRegionIterator <OutputImageType> OutImageIteratorType;
173  OutImageIteratorType outIterator(this->GetOutput(),outputRegionForThread);
174 
175  typedef itk::ImageRegionIterator <OutputB0ImageType> OutB0ImageIteratorType;
176  OutB0ImageIteratorType outB0Iterator(m_EstimatedB0Image,outputRegionForThread);
177  OutB0ImageIteratorType outVarianceIterator(m_EstimatedVarianceImage,outputRegionForThread);
178 
179  typedef typename OutputImageType::PixelType OutputPixelType;
180  std::vector <double> dwi (numInputs,0);
181  std::vector <double> lnDwi (numInputs,0);
182  OutputPixelType resVec(m_NumberOfComponents);
183 
184  std::vector <double> predictedValues(numInputs,0);
186  data.rotationMatrix.set_size(3,3);
187  data.workEigenValues.set_size(3);
188  data.workTensor.set_size(3,3);
189 
190  typedef itk::SymmetricEigenAnalysis < vnl_matrix <double>, vnl_diag_matrix<double>, vnl_matrix <double> > EigenAnalysisType;
191  EigenAnalysisType eigen(3);
192 
193  while (!outIterator.IsAtEnd())
194  {
195  resVec.Fill(0.0);
196 
197  if (maskIterator.Get() == 0)
198  {
199  outIterator.Set(resVec);
200  outB0Iterator.Set(0);
201  outVarianceIterator.Set(0);
202 
203  for (unsigned int i = 0;i < numInputs;++i)
204  ++inIterators[i];
205 
206  ++maskIterator;
207  ++outIterator;
208  ++outB0Iterator;
209  ++outVarianceIterator;
210 
211  continue;
212  }
213 
214  for (unsigned int i = 0;i < numInputs;++i)
215  {
216  dwi[i] = inIterators[i].Get();
217  lnDwi[i] = std::log(std::max(1.0e-6,dwi[i]));
218  }
219 
220  for (unsigned int i = 0;i < m_NumberOfComponents;++i)
221  {
222  for (unsigned int j = 0;j < numInputs;++j)
223  resVec[i] += m_InitialMatrixSolver(i+1,j) * lnDwi[j];
224  }
225 
227 
228  eigen.ComputeEigenValuesAndVectors(data.workTensor,data.workEigenValues,data.rotationMatrix);
229  if (vnl_determinant (data.rotationMatrix) < 0)
230  data.rotationMatrix *= -1;
231 
232  std::vector <double> optimizedValue(m_NumberOfComponents, 0.0);
233 
234  double cThetaControl = 0;
235  for (unsigned int i = 0;i < 3;++i)
236  cThetaControl += data.rotationMatrix(i,i);
237 
238  if (std::abs(cThetaControl + 1.0) > 1.0e-5)
239  anima::Get3DRotationLogarithm(data.rotationMatrix,optimizedValue);
240 
241  double minValue = 1.0e-7;
242  double maxValue = 1.0e-2;
243  for (unsigned int i = 0;i < 3;++i)
244  {
245  optimizedValue[i] += M_PI;
246  int num2Pi = std::floor(optimizedValue[i] / (2.0 * M_PI));
247  optimizedValue[i] -= 2.0 * M_PI * num2Pi + M_PI;
248 
249  optimizedValue[i + 3] = std::min(maxValue, std::max(data.workEigenValues[i], minValue));
250  }
251 
252  // NLOPT optimization
253  nlopt::opt opt(nlopt::LN_BOBYQA, m_NumberOfComponents);
254 
255  std::vector <double> lowerBounds(m_NumberOfComponents, - M_PI);
256  for (unsigned int i = 0;i < 3;++i)
257  lowerBounds[i + 3] = minValue;
258 
259  opt.set_lower_bounds(lowerBounds);
260 
261  std::vector <double> upperBounds(m_NumberOfComponents, M_PI);
262  for (unsigned int i = 0;i < 3;++i)
263  upperBounds[i + 3] = maxValue;
264 
265  opt.set_upper_bounds(upperBounds);
266  opt.set_xtol_rel(1e-4);
267  opt.set_ftol_rel(1e-4);
268  opt.set_maxeval(2500);
269 
270  double minf;
271 
272  data.filter = this;
273  data.dwi = dwi;
274  data.predictedValues = predictedValues;
275 
276  opt.set_min_objective(OptimizationFunction, &data);
277 
278  try
279  {
280  opt.optimize(optimizedValue, minf);
281  }
282  catch(nlopt::roundoff_limited& e)
283  {
284  bool failedOpt = false;
285  for (unsigned int i = 0;i < 6;++i)
286  {
287  if (!std::isfinite(optimizedValue[i]))
288  {
289  failedOpt = true;
290  break;
291  }
292  }
293 
294  if (failedOpt)
295  {
296  resVec.Fill(0.0);
297  outIterator.Set(resVec);
298  outB0Iterator.Set(0);
299  outVarianceIterator.Set(0);
300 
301  for (unsigned int i = 0;i < numInputs;++i)
302  ++inIterators[i];
303 
304  ++maskIterator;
305  ++outIterator;
306  ++outB0Iterator;
307  ++outVarianceIterator;
308  this->IncrementNumberOfProcessedPoints();
309 
310  continue;
311  }
312  }
313 
314  anima::Get3DRotationExponential(optimizedValue,data.rotationMatrix);
315  for (unsigned int i = 0;i < 3;++i)
316  data.workEigenValues[i] = optimizedValue[3 + i];
317 
320 
321  double outVarianceValue;
322  double outB0Value = this->ComputeB0AndVarianceFromTensorVector(data.workTensor,dwi,outVarianceValue);
323 
324  outIterator.Set(resVec);
325  outB0Iterator.Set(outB0Value);
326  outVarianceIterator.Set(outVarianceValue);
327 
328  for (unsigned int i = 0;i < numInputs;++i)
329  ++inIterators[i];
330 
331  this->IncrementNumberOfProcessedPoints();
332  ++maskIterator;
333  ++outIterator;
334  ++outB0Iterator;
335  ++outVarianceIterator;
336  }
337 }
338 
339 template <class InputPixelScalarType, class OutputPixelScalarType>
340 double
342 ::ComputeB0AndVarianceFromTensorVector(const vnl_matrix <double> &tensorValue, const std::vector <double> &dwiSignal, double &outVarianceValue)
343 {
344  double sumSquaredObservedSignals = 0;
345  double sumPredictedPerObservedSignals = 0;
346  double sumSquaredPredictedSignals = 0;
347  unsigned int numInputs = dwiSignal.size();
348 
349  for (unsigned int i = 0;i < numInputs;++i)
350  {
351  double observedValue = dwiSignal[i];
352  double predictedValue = 0;
353  double bValue = m_BValuesList[i];
354 
355  if (bValue == 0)
356  predictedValue = 1;
357  else
358  {
359  for (unsigned int j = 0;j < 3;++j)
360  {
361  predictedValue += tensorValue(j,j) * m_GradientDirections[i][j] * m_GradientDirections[i][j];
362  for (unsigned int k = j + 1;k < 3;++k)
363  predictedValue += 2.0 * tensorValue(j,k) * m_GradientDirections[i][j] * m_GradientDirections[i][k];
364  }
365 
366  predictedValue = std::exp(- bValue * predictedValue);
367  }
368 
369  sumSquaredObservedSignals += observedValue * observedValue;
370  sumSquaredPredictedSignals += predictedValue * predictedValue;
371  sumPredictedPerObservedSignals += predictedValue * observedValue;
372  }
373 
374  double outB0Value = sumPredictedPerObservedSignals / sumSquaredPredictedSignals;
375  outVarianceValue = (sumSquaredObservedSignals - 2.0 * outB0Value * sumPredictedPerObservedSignals + outB0Value * outB0Value * sumSquaredPredictedSignals) / (numInputs - 1.0);
376  return outB0Value;
377 }
378 
379 template <class InputPixelScalarType, class OutputPixelScalarType>
380 double
382 ::OptimizationFunction(const std::vector<double> &x, std::vector<double> &grad, void *func_data)
383 {
384  OptimizationDataStructure *data = static_cast <OptimizationDataStructure *> (func_data);
385 
386  return data->filter->ComputeCostAtPosition(x,data->dwi,data->predictedValues,data->rotationMatrix,data->workTensor,data->workEigenValues);
387 }
388 
389 template <class InputPixelScalarType, class OutputPixelScalarType>
390 double
392 ::ComputeCostAtPosition(const std::vector<double> &x, const std::vector <double> &observedData,
393  std::vector <double> &predictedValues, vnl_matrix <double> &rotationMatrix,
394  vnl_matrix<double> &workTensor, vnl_diag_matrix <double> &workEigenValues)
395 {
396  anima::Get3DRotationExponential(x,rotationMatrix);
397  workEigenValues.set_size(3);
398  for (unsigned int i = 0;i < 3;++i)
399  workEigenValues[i] = x[3 + i];
400 
401  workTensor.set_size(3,3);
402  anima::RecomposeTensor(workEigenValues,rotationMatrix,workTensor);
403 
404  for (unsigned int i = 0;i < this->GetNumberOfIndexedInputs();++i)
405  {
406  predictedValues[i] = 0;
407  double bValue = m_BValuesList[i];
408  if (bValue == 0)
409  {
410  predictedValues[i] = 1;
411  continue;
412  }
413 
414  for (unsigned int j = 0;j < 3;++j)
415  {
416  predictedValues[i] += workTensor(j,j) * m_GradientDirections[i][j] * m_GradientDirections[i][j];
417  for (unsigned int k = j + 1;k < 3;++k)
418  predictedValues[i] += 2.0 * workTensor(j,k) * m_GradientDirections[i][j] * m_GradientDirections[i][k];
419  }
420 
421  predictedValues[i] = std::exp(- bValue * predictedValues[i]);
422  }
423 
424  double b0Val = 0;
425  double normConstant = 0;
426 
427  for (unsigned int i = 0;i < this->GetNumberOfIndexedInputs();++i)
428  {
429  b0Val += predictedValues[i] * observedData[i];
430  normConstant += predictedValues[i] * predictedValues[i];
431  }
432 
433  b0Val /= normConstant;
434 
435  double costValue = 0;
436  for (unsigned int i = 0;i < this->GetNumberOfIndexedInputs();++i)
437  {
438  double predVal = b0Val * predictedValues[i];
439  costValue += (predVal - observedData[i]) * (predVal - observedData[i]);
440  }
441 
442  return costValue;
443 }
444 
445 } // end of namespace anima
void AddGradientDirection(unsigned int i, vnl_vector_fixed< double, 3 > &grad)
double ComputeB0AndVarianceFromTensorVector(const vnl_matrix< double > &tensorValue, const std::vector< double > &dwiSignal, double &outVarianceValue)
void Get3DRotationLogarithm(const vnl_matrix< T > &rotationMatrix, std::vector< T > &outputAngles)
Computation of a 3D rotation matrix logarithm. Rodrigues&#39; explicit formula.
Superclass::OutputImageRegionType OutputImageRegionType
itk::VectorImage< OutputPixelScalarType, 3 > OutputImageType
static double OptimizationFunction(const std::vector< double > &x, std::vector< double > &grad, void *func_data)
void DynamicThreadedGenerateData(const OutputImageRegionType &outputRegionForThread) ITK_OVERRIDE
void GetVectorRepresentation(const vnl_matrix< T1 > &tensor, itk::VariableLengthVector< T2 > &vector, unsigned int vecDim=0, bool scale=false)
void GetTensorFromVectorRepresentation(const itk::VariableLengthVector< T1 > &vector, vnl_matrix< T2 > &tensor, unsigned int tensDim=0, bool scale=false)
void RecomposeTensor(vnl_diag_matrix< T1 > &eigs, vnl_matrix< T1 > &eigVecs, vnl_matrix< T2 > &resMatrix)
Recompose tensor from values extracted using SymmetricEigenAnalysis (vnl_symmetric_eigensystem transp...
void Get3DRotationExponential(const std::vector< T > &angles, vnl_matrix< T > &outputRotation)
Computation of a 3D rotation matrix exponential. Rodrigues&#39; explicit formula.
double ComputeCostAtPosition(const std::vector< double > &x, const std::vector< double > &observedData, std::vector< double > &predictedValues, vnl_matrix< double > &rotationMatrix, vnl_matrix< double > &workTensor, vnl_diag_matrix< double > &workEigenValues)