ANIMA  4.0
animaMultiT2RelaxometryEstimationImageFilter.hxx
Go to the documentation of this file.
1 #pragma once
3 
4 #include <itkImageRegionIterator.h>
5 #include <itkImageRegionConstIterator.h>
6 
8 #include <animaNLOPTOptimizers.h>
11 
12 namespace anima
13 {
14 
15 template <class TPixelScalarType>
16 void
17 MultiT2RelaxometryEstimationImageFilter <TPixelScalarType>
18 ::BeforeThreadedGenerateData()
19 {
20  if (m_NLEstimation && (!m_InitialT2Map || !m_InitialM0Map))
21  itkExceptionMacro("Missing inputs for non-local estimation");
22 
23  Superclass::BeforeThreadedGenerateData();
24 
25  m_T2CompartmentValues.resize(m_NumberOfT2Compartments);
26  double logStart = std::log(m_LowerT2Bound);
27  double logEnd = std::log(m_UpperT2Bound);
28  double step = (logEnd - logStart) / (m_NumberOfT2Compartments - 1.0);
29 
30  double logValue = logStart;
31  m_T2CompartmentValues[0] = std::exp(logValue);
32  for (unsigned int i = 1;i < m_NumberOfT2Compartments;++i)
33  {
34  logValue += step;
35  m_T2CompartmentValues[i] = std::exp(logValue);
36  }
37 
38  this->GetB1OutputImage()->FillBuffer(1.0);
39 
40  m_T2OutputImage = VectorOutputImageType::New();
41  m_T2OutputImage->Initialize();
42  m_T2OutputImage->SetRegions(this->GetInput(0)->GetLargestPossibleRegion());
43  m_T2OutputImage->SetSpacing (this->GetInput(0)->GetSpacing());
44  m_T2OutputImage->SetOrigin (this->GetInput(0)->GetOrigin());
45  m_T2OutputImage->SetDirection (this->GetInput(0)->GetDirection());
46  m_T2OutputImage->SetVectorLength(m_NumberOfT2Compartments);
47  m_T2OutputImage->Allocate();
48 
49  OutputVectorType zero(m_NumberOfT2Compartments);
50  zero.Fill(0.0);
51  m_T2OutputImage->FillBuffer(zero);
52 
53  if (m_NLEstimation)
54  {
55  this->PrepareNLPatchSearchers();
56 
57  typedef itk::ImageRegionConstIterator <InputImageType> InIteratorType;
58  typedef itk::ImageRegionIterator <VectorOutputImageType> VectorIteratorType;
59 
60  InIteratorType m0Iterator(m_InitialM0Map, m_InitialM0Map->GetLargestPossibleRegion());
61  VectorIteratorType t2Iterator(m_InitialT2Map, m_InitialT2Map->GetLargestPossibleRegion());
62  OutputVectorType tmpData;
63 
64  while (!m0Iterator.IsAtEnd())
65  {
66  tmpData = t2Iterator.Get();
67  double m0Value = m0Iterator.Get();
68 
69  for (unsigned int i = 0;i < m_NumberOfT2Compartments;++i)
70  tmpData[i] *= m0Value;
71 
72  t2Iterator.Set(tmpData);
73  ++m0Iterator;
74  ++t2Iterator;
75  }
76  }
77 }
78 
79 template <class TPixelScalarType>
80 void
82 ::CheckComputationMask()
83 {
84  if (this->GetComputationMask())
85  return;
86 
87  typedef itk::ImageRegionConstIterator <InputImageType> IteratorType;
88  typedef itk::ImageRegionIterator <MaskImageType> MaskIteratorType;
89 
90  std::vector <IteratorType> inItrs(this->GetNumberOfIndexedInputs());
91  for (unsigned int i = 0;i < this->GetNumberOfIndexedInputs();++i)
92  inItrs[i] = IteratorType(this->GetInput(i),this->GetOutput(0)->GetLargestPossibleRegion());
93 
94  typename MaskImageType::Pointer maskImage = MaskImageType::New();
95  maskImage->Initialize();
96  maskImage->SetRegions(this->GetInput(0)->GetLargestPossibleRegion());
97  maskImage->SetSpacing (this->GetInput(0)->GetSpacing());
98  maskImage->SetOrigin (this->GetInput(0)->GetOrigin());
99  maskImage->SetDirection (this->GetInput(0)->GetDirection());
100  maskImage->Allocate();
101 
102  MaskIteratorType maskItr (maskImage,this->GetOutput(0)->GetLargestPossibleRegion());
103  while (!maskItr.IsAtEnd())
104  {
105  double averageVal = 0;
106  for (unsigned int i = 0;i < this->GetNumberOfIndexedInputs();++i)
107  averageVal += inItrs[i].Get();
108 
109  averageVal /= this->GetNumberOfIndexedInputs();
110 
111  bool maskPoint = (averageVal <= m_AverageSignalThreshold);
112 
113  if (maskPoint)
114  maskItr.Set(0);
115  else
116  maskItr.Set(1);
117 
118  for (unsigned int i = 0;i < this->GetNumberOfIndexedInputs();++i)
119  ++inItrs[i];
120 
121  ++maskItr;
122  }
123 
124  this->SetComputationMask(maskImage);
125 }
126 
127 template <class TPixelScalarType>
128 void
130 ::DynamicThreadedGenerateData(const OutputImageRegionType &outputRegionForThread)
131 {
132  typedef itk::ImageRegionConstIterator <InputImageType> ImageConstIteratorType;
133  typedef itk::ImageRegionIterator <InputImageType> ImageIteratorType;
134  typedef itk::ImageRegionIterator <VectorOutputImageType> VectorImageIteratorType;
135 
136  VectorImageIteratorType outT2Iterator(this->GetT2OutputImage(),outputRegionForThread);
137  ImageIteratorType outM0Iterator(this->GetM0OutputImage(),outputRegionForThread);
138  ImageIteratorType outMWFIterator(this->GetMWFOutputImage(),outputRegionForThread);
139  ImageIteratorType outB1Iterator(this->GetB1OutputImage(),outputRegionForThread);
140  ImageIteratorType outCostIterator(this->GetCostOutputImage(),outputRegionForThread);
141 
142  VectorImageIteratorType initT2Iterator;
143  ImageIteratorType initM0Iterator;
144  ImageIteratorType initB1Iterator;
145 
146  if (m_NLEstimation)
147  {
148  initT2Iterator = VectorImageIteratorType(m_InitialT2Map, outputRegionForThread);
149  initM0Iterator = ImageIteratorType(m_InitialM0Map, outputRegionForThread);
150  initB1Iterator = ImageIteratorType(m_InitialB1Map, outputRegionForThread);
151  }
152 
153  unsigned int numInputs = this->GetNumberOfIndexedInputs();
154  std::vector <ImageConstIteratorType> inIterators;
155  for (unsigned int i = 0;i < numInputs;++i)
156  inIterators.push_back(ImageConstIteratorType(this->GetInput(i),outputRegionForThread));
157 
158  typedef itk::ImageRegionIterator <MaskImageType> MaskIteratorType;
159  MaskIteratorType maskItr(this->GetComputationMask(),outputRegionForThread);
160 
161  ImageIteratorType t1MapItr;
162  if (m_T1Map)
163  t1MapItr = ImageIteratorType(m_T1Map,outputRegionForThread);
164 
165  T2VectorType signalValues(numInputs);
166  T2VectorType signalValuesExtended(numInputs + m_NumberOfT2Compartments);
167  T2VectorType t2OptimizedWeights(m_NumberOfT2Compartments);
168  T2VectorType priorDistribution(m_NumberOfT2Compartments);
169  priorDistribution.fill(0);
170 
171  OutputVectorType outputT2Weights(m_NumberOfT2Compartments);
172 
173  DataMatrixType AMatrix(numInputs,m_NumberOfT2Compartments,0);
174  DataMatrixType AMatrixExtended(numInputs + m_NumberOfT2Compartments,m_NumberOfT2Compartments,0);
175 
176  anima::EPGSignalSimulator epgSimulator;
177  epgSimulator.SetEchoSpacing(m_EchoSpacing);
178  epgSimulator.SetNumberOfEchoes(numInputs);
179  epgSimulator.SetExcitationFlipAngle(m_T2ExcitationFlipAngle);
180 
181  anima::EPGSignalSimulator::RealVectorType epgSignalValues(numInputs);
182 
183  NNLSOptimizerPointer nnlsOpt = NNLSOptimizerType::New();
184 
185  typedef anima::NLOPTOptimizers B1OptimizerType;
186  typedef anima::MultiT2EPGRelaxometryCostFunction B1CostFunctionType;
187 
188  typename B1CostFunctionType::Pointer cost = B1CostFunctionType::New();
189  cost->SetEchoSpacing(m_EchoSpacing);
190  cost->SetExcitationFlipAngle(m_T2ExcitationFlipAngle);
191  cost->SetT2FlipAngles(m_T2FlipAngles);
192  cost->SetM0Value(1.0);
193 
194  unsigned int dimension = cost->GetNumberOfParameters();
195  itk::Array<double> lowerBounds(dimension);
196  itk::Array<double> upperBounds(dimension);
197  B1OptimizerType::ParametersType p(dimension);
198  lowerBounds[0] = 0;
199  upperBounds[0] = 1;
200 
201  // NL specific variables
202  std::vector <double> workDataWeights;
203  std::vector <OutputVectorType> workDataSamples;
204 
205  unsigned int threadId = this->GetSafeThreadId();
206 
207  while (!maskItr.IsAtEnd())
208  {
209  outputT2Weights.Fill(0);
210 
211  if (maskItr.Get() == 0)
212  {
213  outT2Iterator.Set(outputT2Weights);
214  outM0Iterator.Set(0);
215  outMWFIterator.Set(0.0);
216  outB1Iterator.Set(0.0);
217  outCostIterator.Set(0.0);
218 
219  ++maskItr;
220  ++outT2Iterator;
221  ++outM0Iterator;
222  ++outMWFIterator;
223  ++outB1Iterator;
224  ++outCostIterator;
225 
226  for (unsigned int i = 0;i < numInputs;++i)
227  ++inIterators[i];
228 
229  if (m_T1Map)
230  ++t1MapItr;
231 
232  if (m_NLEstimation)
233  {
234  ++initB1Iterator;
235  ++initT2Iterator;
236  ++initM0Iterator;
237  }
238 
239  continue;
240  }
241 
242  double previousB1Value = -1.0;
243  double t1Value = 1000;
244  double m0Value = 0.0;
245 
246  for (unsigned int i = 0;i < numInputs;++i)
247  {
248  signalValues[i] = inIterators[i].Get();
249  signalValuesExtended[i] = signalValues[i];
250  }
251 
252  if (m_T1Map)
253  {
254  t1Value = t1MapItr.Get();
255  if (t1Value <= 0.0)
256  t1Value = 1000;
257  }
258 
259  cost->SetT1Value(t1Value);
260  cost->SetT2Values(m_T2CompartmentValues);
261  cost->SetT2RelaxometrySignals(signalValues);
262 
263  double b1Value;
264 
265  if (m_InitialB1Map)
266  b1Value = initB1Iterator.Get();
267  else
268  b1Value = outB1Iterator.Get();
269 
270  if (m_NLEstimation)
271  {
272  outputT2Weights = initT2Iterator.Get();
273  this->ComputeTikhonovPrior(maskItr.GetIndex(),outputT2Weights,m_NLPatchSearchers[threadId],priorDistribution,
274  workDataWeights, workDataSamples);
275  }
276 
277  unsigned int numGlobalIterations = 0;
278  double residual = 0;
279 
280  while ((std::abs(b1Value - previousB1Value) > m_B1Tolerance)&&(numGlobalIterations < 100))
281  {
282  ++numGlobalIterations;
283  previousB1Value = b1Value;
284  // T2 weights estimation
285  AMatrix.fill(0);
286  AMatrixExtended.fill(0);
287  for (unsigned int i = 0;i < m_NumberOfT2Compartments;++i)
288  {
289  epgSignalValues = epgSimulator.GetValue(t1Value,m_T2CompartmentValues[i],b1Value,1.0);
290  for (unsigned int j = 0;j < numInputs;++j)
291  {
292  AMatrix(j,i) = epgSignalValues[j];
293  AMatrixExtended(j,i) = epgSignalValues[j];
294  }
295  }
296 
297  // Regular NNLS optimization
298  nnlsOpt->SetDataMatrix(AMatrix);
299  nnlsOpt->SetPoints(signalValues);
300 
301  nnlsOpt->StartOptimization();
302  t2OptimizedWeights = nnlsOpt->GetCurrentPosition();
303  residual = nnlsOpt->GetCurrentResidual();
304 
305  // Regularized NNLS with or without prior
306  double normT2Weights = 0;
307  for (unsigned int i = 0;i < m_NumberOfT2Compartments;++i)
308  normT2Weights += (t2OptimizedWeights[i] - priorDistribution[i]) * (t2OptimizedWeights[i] - priorDistribution[i]);
309 
310  if (m_RegularizationIntensity > 1.0)
311  {
312  double lambdaSq = (m_RegularizationIntensity - 1.0) * residual / normT2Weights;
313  residual = this->ComputeTikhonovRegularizedSolution(nnlsOpt,AMatrixExtended,signalValuesExtended,
314  lambdaSq,priorDistribution,t2OptimizedWeights);
315  }
316 
317  outCostIterator.Set(residual);
318 
319  m0Value = 0.0;
320  for (unsigned int i = 0;i < m_NumberOfT2Compartments;++i)
321  m0Value += t2OptimizedWeights[i];
322 
323  for (unsigned int i = 0;i < m_NumberOfT2Compartments;++i)
324  outputT2Weights[i] = t2OptimizedWeights[i] / m0Value;
325 
326  // B1 value estimation
327  cost->SetT2Weights(t2OptimizedWeights);
328 
329  B1OptimizerType::Pointer b1Optimizer = B1OptimizerType::New();
330  b1Optimizer->SetAlgorithm(NLOPT_LN_BOBYQA);
331  b1Optimizer->SetXTolRel(1.0e-4);
332  b1Optimizer->SetFTolRel(1.0e-6);
333  b1Optimizer->SetMaxEval(500);
334  b1Optimizer->SetVectorStorageSize(2000);
335 
336  b1Optimizer->SetLowerBoundParameters(lowerBounds);
337  b1Optimizer->SetUpperBoundParameters(upperBounds);
338  p[0] = b1Value;
339  b1Optimizer->SetInitialPosition(p);
340  b1Optimizer->SetMaximize(false);
341  b1Optimizer->SetCostFunction(cost);
342 
343  b1Optimizer->StartOptimization();
344  p = b1Optimizer->GetCurrentPosition();
345  b1Value = p[0];
346  }
347 
348  outM0Iterator.Set(m0Value);
349  outT2Iterator.Set(outputT2Weights);
350  outCostIterator.Set(residual);
351  double mwfValue = 0;
352  for (unsigned int i = 0;i < m_NumberOfT2Compartments;++i)
353  {
354  if (m_T2CompartmentValues[i] > m_MyelinThreshold)
355  break;
356 
357  mwfValue += outputT2Weights[i];
358  }
359 
360  outMWFIterator.Set(mwfValue);
361  outB1Iterator.Set(b1Value);
362 
363  this->IncrementNumberOfProcessedPoints();
364  ++maskItr;
365  ++outT2Iterator;
366  ++outM0Iterator;
367  ++outMWFIterator;
368  ++outB1Iterator;
369  ++outCostIterator;
370 
371  for (unsigned int i = 0;i < numInputs;++i)
372  ++inIterators[i];
373 
374  if (m_T1Map)
375  ++t1MapItr;
376 
377  if (m_NLEstimation)
378  {
379  ++initB1Iterator;
380  ++initT2Iterator;
381  ++initM0Iterator;
382  }
383  }
384 
385  this->SafeReleaseThreadId(threadId);
386 }
387 
388 template <class TPixelScalarType>
389 void
391 ::PrepareNLPatchSearchers()
392 {
393  unsigned int numThreads = this->GetNumberOfWorkUnits();
394  m_NLPatchSearchers.resize(numThreads);
395  unsigned int maxAbsDisp = floor((double)(m_SearchNeighborhood / m_SearchStepSize)) * m_SearchStepSize;
396 
397  // Prepare mean and variance images
398  unsigned int numInputs = this->GetNumberOfIndexedInputs();
399  typedef anima::MeanAndVarianceImagesFilter<InputImageType, OutputImageType> MeanAndVarianceImagesFilterType;
400 
401  typename InputImageType::SizeType radius;
402  for (unsigned int j = 0;j < InputImageType::ImageDimension;++j)
403  radius[j] = m_PatchHalfSize;
404 
405  for (unsigned int i = 0;i < numInputs;++i)
406  {
407  typename MeanAndVarianceImagesFilterType::Pointer filter = MeanAndVarianceImagesFilterType::New();
408  filter->SetInput(this->GetInput(i));
409 
410  filter->SetRadius(radius);
411  filter->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits());
412  filter->Update();
413 
414  OutputImagePointer meanImage = filter->GetMeanImage();
415  meanImage->DisconnectPipeline();
416  OutputImagePointer varImage = filter->GetVarImage();
417  varImage->DisconnectPipeline();
418 
419  for (unsigned int j = 0;j < numThreads;++j)
420  {
421  m_NLPatchSearchers[j].AddMeanImage(meanImage);
422  m_NLPatchSearchers[j].AddVarImage(varImage);
423  m_NLPatchSearchers[j].AddPatchTestImage(this->GetInput(i));
424  }
425  }
426 
427  // Prepare noise covariance values
428  typename InputImageRegionType::IndexType baseIndex;
429  typename InputImageRegionType::IndexType valueIndex;
430  std::vector <double> noiseCovariances(numInputs,0);
431 
432  for (unsigned int i = 0;i < numInputs;++i)
433  {
434  typedef itk::ImageRegionConstIteratorWithIndex< InputImageType > InIteratorType;
435  typename InputImageType::RegionType largestRegion = this->GetInput(i)->GetLargestPossibleRegion();
436  InIteratorType dataIterator (this->GetInput(i), largestRegion);
437 
438  double averageLocalSignal, diffSignal;
439  double baseSignal;
440 
441  double averageCovariance = 0;
442  unsigned int numEstimations = 0;
443  unsigned int numLocalPixels = 2 * InputImageType::ImageDimension;
444 
445  while (!dataIterator.IsAtEnd())
446  {
447  baseSignal = static_cast<double>(dataIterator.Get());
448  baseIndex = dataIterator.GetIndex();
449  averageLocalSignal = 0;
450 
451  for (unsigned int d = 0; d < InputImageType::ImageDimension; ++d)
452  {
453  valueIndex = baseIndex;
454  int tmpIndex = baseIndex[d] - m_LocalNeighborhood;
455  valueIndex[d] = std::max(tmpIndex,0);
456  averageLocalSignal += static_cast<double> (this->GetInput(i)->GetPixel(valueIndex));
457 
458  valueIndex = baseIndex;
459  tmpIndex = baseIndex[d] + m_LocalNeighborhood;
460  int maxIndex = largestRegion.GetSize()[d] - 1;
461  valueIndex[d] = std::min(tmpIndex, maxIndex);
462  averageLocalSignal += static_cast<double> (this->GetInput(i)->GetPixel(valueIndex));
463  }
464 
465  averageLocalSignal /= numLocalPixels;
466  diffSignal = sqrt(numLocalPixels / (numLocalPixels + 1.0)) * (baseSignal - averageLocalSignal);
467 
468  averageCovariance += diffSignal * diffSignal;
469 
470  ++numEstimations;
471  ++dataIterator;
472  }
473 
474  // Now divide by number of estimations and compute average variance
475  noiseCovariances[i] = averageCovariance / numEstimations;
476  }
477 
478  //Set searcher parameters
479  for (unsigned int i = 0;i < numThreads;++i)
480  {
481  m_NLPatchSearchers[i].SetPatchHalfSize(m_PatchHalfSize);
482  m_NLPatchSearchers[i].SetSearchStepSize(m_SearchStepSize);
483  m_NLPatchSearchers[i].SetMaxAbsDisp(maxAbsDisp);
484  m_NLPatchSearchers[i].SetInputImage(m_InitialT2Map);
485  m_NLPatchSearchers[i].SetBetaParameter(m_BetaParameter);
486  m_NLPatchSearchers[i].SetWeightThreshold(m_WeightThreshold);
487  m_NLPatchSearchers[i].SetMeanMinThreshold(m_MeanMinThreshold);
488  m_NLPatchSearchers[i].SetVarMinThreshold(m_VarMinThreshold);
489  m_NLPatchSearchers[i].SetNoiseCovariances(noiseCovariances);
490  }
491 }
492 
493 template <class TPixelScalarType>
494 void
496 ::ComputeTikhonovPrior(const IndexType &refIndex, OutputVectorType &refDistribution,
497  PatchSearcherType &nlPatchSearcher, T2VectorType &priorDistribution,
498  std::vector <double> &workDataWeights, std::vector <OutputVectorType> &workDataSamples)
499 {
500  nlPatchSearcher.UpdateAtPosition(refIndex);
501 
502  workDataWeights = nlPatchSearcher.GetDatabaseWeights();
503  workDataSamples = nlPatchSearcher.GetDatabaseSamples();
504 
505  priorDistribution.fill(0);
506 
507  double sumWeights = 0;
508  double maxWeight = -1;
509  for (unsigned int i = 0;i < workDataSamples.size();++i)
510  {
511  bool isZero = true;
512  for (unsigned int j = 0;j < m_NumberOfT2Compartments;++j)
513  {
514  if (workDataSamples[i][j] != 0)
515  {
516  isZero = false;
517  break;
518  }
519  }
520 
521  if (isZero)
522  continue;
523 
524  for (unsigned int j = 0;j < m_NumberOfT2Compartments;++j)
525  priorDistribution[j] += workDataWeights[i] * workDataSamples[i][j];
526 
527  if (maxWeight < workDataWeights[i])
528  maxWeight = workDataWeights[i];
529 
530  sumWeights += workDataWeights[i];
531  }
532 
533  if (maxWeight < 0)
534  maxWeight = 1.0;
535 
536  for (unsigned int j = 0;j < m_NumberOfT2Compartments;++j)
537  priorDistribution[j] += maxWeight * refDistribution[j];
538 
539  sumWeights += maxWeight;
540  priorDistribution /= sumWeights;
541 }
542 
543 template <class TPixelScalarType>
544 double
546 ::ComputeTikhonovRegularizedSolution(anima::NNLSOptimizer *nnlsOpt, DataMatrixType &AMatrix,
547  T2VectorType &signalValues, double lambdaSq,
548  T2VectorType &priorDistribution, T2VectorType &t2OptimizedWeights)
549 {
550  unsigned int rowSize = AMatrix.rows() - AMatrix.cols();
551  unsigned int colSize = AMatrix.cols();
552 
553  double lambda = std::sqrt(lambdaSq);
554 
555  for (unsigned int i = 0;i < colSize;++i)
556  {
557  signalValues[rowSize + i] = lambda * priorDistribution[i];
558  for (unsigned int j = 0;j < colSize;++j)
559  {
560  if (i != j)
561  AMatrix(rowSize + i,j) = 0;
562  else
563  AMatrix(rowSize + i,j) = lambda;
564  }
565  }
566 
567  nnlsOpt->SetDataMatrix(AMatrix);
568  nnlsOpt->SetPoints(signalValues);
569 
570  nnlsOpt->StartOptimization();
571 
572  t2OptimizedWeights = nnlsOpt->GetCurrentPosition();
573 
574  return nnlsOpt->GetCurrentResidual();
575 }
576 
577 } // end namespace anima
Applies an variance filter to an image.
virtual const std::vector< double > & GetDatabaseWeights() const
Implements an ITK wrapper for the NLOPT library.
void UpdateAtPosition(const IndexType &dataIndex)
void SetPoints(const ParametersType &points)
Cost function for estimating B1 from T2 relaxometry acquisition, following a multi-T2 EPG decay model...
void SetDataMatrix(const MatrixType &data)
std::vector< double > RealVectorType
Implements multi-peak T2 relaxometry estimation (with or without regularization)
void StartOptimization() ITK_OVERRIDE
Non negative least squares optimizer. Implements Lawson et al method, of squared problem is activated...
virtual const std::vector< PixelType > & GetDatabaseSamples() const
RealVectorType & GetValue(double t1Value, double t2Value, double flipAngle, double m0Value)
Get EPG values at given point.
void SetNumberOfEchoes(unsigned int val)