ANIMA  4.0
animaHierarchicalInitializer.hxx
Go to the documentation of this file.
2 
3 namespace anima
4 {
5 
6 template <typename TInputImage, typename TMaskImage>
8 {
9  this->SetNthInput(0, const_cast<TMaskImage*>(mask));
10 }
11 
12 template <typename TInputImage, typename TMaskImage>
14 {
15  this->SetNthInput(m_NbInputs, const_cast<TInputImage*>(image));
16  m_IndexImage1 = m_NbInputs;
17  m_NbInputs++;
18 }
19 
20 template <typename TInputImage, typename TMaskImage>
22 {
23  this->SetNthInput(m_NbInputs, const_cast<TInputImage*>(image));
24  m_IndexImage2 = m_NbInputs;
25  m_NbInputs++;
26 }
27 
28 template <typename TInputImage, typename TMaskImage>
30 {
31  this->SetNthInput(m_NbInputs, const_cast<TInputImage*>(image));
32  m_IndexImage3 = m_NbInputs;
33  m_NbInputs++;
34 }
35 template <typename TInputImage, typename TMaskImage>
37 {
38  this->SetNthInput(m_NbInputs, const_cast<TInputImage*>(image));
39  m_IndexImage4 = m_NbInputs;
40  m_NbInputs++;
41 }
42 
43 template <typename TInputImage, typename TMaskImage>
45 {
46  this->SetNthInput(m_NbInputs, const_cast<TInputImage*>(image));
47  m_IndexImage5 = m_NbInputs;
48  m_NbInputs++;
49 }
50 
51 
52 template <typename TInputImage, typename TMaskImage>
54 {
55  return static_cast< const TMaskImage * >
56  ( this->itk::ProcessObject::GetInput(0) );
57 }
58 
59 template <typename TInputImage, typename TMaskImage>
61 {
62  return static_cast< const TInputImage * >
63  ( this->itk::ProcessObject::GetInput(m_IndexImage1) );
64 }
65 
66 template <typename TInputImage, typename TMaskImage>
68 {
69  return static_cast< const TInputImage * >
70  ( this->itk::ProcessObject::GetInput(m_IndexImage2) );
71 }
72 
73 template <typename TInputImage, typename TMaskImage>
75 {
76  return static_cast< const TInputImage * >
77  ( this->itk::ProcessObject::GetInput(m_IndexImage3) );
78 }
79 
80 template <typename TInputImage, typename TMaskImage>
82 {
83  return static_cast< const TInputImage * >
84  ( this->itk::ProcessObject::GetInput(m_IndexImage4) );
85 }
86 
87 template <typename TInputImage, typename TMaskImage>
89 {
90  return static_cast< const TInputImage * >
91  ( this->itk::ProcessObject::GetInput(m_IndexImage5) );
92 }
93 
94 template <typename TInputImage, typename TMaskImage>
97 {
98  // Compute EM for T1-w
99  std::vector<double> min1D( 1, 255.0 ); //min values for random Initialization
100  std::vector<double> max1D( 1, 0.0 ); //max values for random Initialization
101 
102  // Decide initializer
104  initiaRandom->SetMinValues( min1D );
105  initiaRandom->SetMaxValues( max1D );
106  initiaRandom->SetDimensionGaussian(1);
107  initiaRandom->SetNbGaussian(3);
108 
109  typename GaussianREMEstimatorType::Pointer estimator= GaussianREMEstimatorType::New();
110  estimator ->SetMaxIterations( 100 );
111  estimator ->SetModelMinDistance( 1e-3 );
112  estimator ->SetMaxIterationsConc( 10 );
113  estimator ->SetStremMode( false );
114  estimator ->SetRejectionRatio( m_Robust );
115  estimator ->SetMask( this->GetMask() );
116  estimator ->SetInputImage1( this->GetInputImage1() );
117  estimator ->SetVerbose( false );
118 
119  std::vector<unsigned int> emSteps( 1, 60 );
120  std::vector<unsigned int> iterSteps( 1, 60 );
121 
122  typename ClassificationStrategyType::Pointer strategy = ClassificationStrategyType::New();
123  strategy->SetEstimator( estimator );
124  strategy->SetInitializer( initiaRandom );
125  strategy->SetStrategy( emSteps, iterSteps );
126  itk::CStyleCommand::Pointer callback = itk::CStyleCommand::New();
127  callback ->SetCallback(eventCallback);
128  strategy ->AddObserver(itk::ProgressEvent(), callback );
129  strategy->Update();
130  std::cout << std::endl;
131 
132  ModelMap solutionClassifier;
133  ModelMap2 solutionClassifierAlpha;
134 
135  if ( !strategy->GetSolutionMap( solutionClassifier, solutionClassifierAlpha ) )
136  {
137  std::cerr << "-- Error: No solution found in first modality" << std::endl;
138  return;
139  }
140 
141  m_Solution1D = solutionClassifier.begin() ->second;
142  m_Solution1DAlphas = solutionClassifierAlpha.begin() ->second;
143 }
144 
145 template <typename TInputImage, typename TMaskImage>
148 {
149  // Classify T1-w
150  typename ImageClassifierType::Pointer classifier = ImageClassifierType::New();
151  classifier->SetInputImage1(this->GetInputImage1());
152  classifier->SetMask(this->GetMask());
153  classifier->SetGaussianModel(m_Solution1D);
154  classifier->SetAlphas(m_Solution1DAlphas);
155  classifier->SetTol(m_Tol);
156  classifier->Update();
157 
158  // creating one image per class (CSF, GM, WM)
159  std::vector<InputIteratorType> classesIt;
160  for(int i = 0; i < 3; i++)
161  {
162  InputImagePointer im = InputImageType::New();
163  im->SetRegions(this->GetMask()->GetLargestPossibleRegion());
164  im->CopyInformation(this->GetMask());
165  im->Allocate();
166  im->FillBuffer(0);
167  m_ImagesClasses.push_back(im);
168 
169  InputIteratorType imIt(im, im->GetLargestPossibleRegion() );
170  classesIt.push_back(imIt);
171  }
172 
173  InputConstIteratorType classifIt(classifier->GetOutput(), classifier->GetOutput()->GetLargestPossibleRegion() );
174  while(!classifIt.IsAtEnd())
175  {
176  if(!(classifIt.Get() == 0 || classifIt.Get() > 3 ))
177  {
178  classesIt[ classifIt.Get() - 1 ].Set(1);
179  }
180  ++classifIt;
181  for(unsigned int i = 0; i < classesIt.size(); i++)
182  {
183  ++classesIt[i];
184  }
185  }
186 }
187 
188 template <typename TInputImage, typename TMaskImage>
191 {
192  const double resolution = 1;
193  const double bandwidth = 5;
194 
195  // selection rules [sequence][class]
196  // +1 brighter max
197  // 0 absolute max
198  // -1 darker max
199  std::vector< std::vector<int> > selectionRules( 3, std::vector<int>( 3, 0 ) );
200  selectionRules[ 1 ][ 0 ] = 1; // T2-w & CSF: brighter max
201  if(m_ThirdIsFLAIR)
202  selectionRules[ 2 ][ 0 ] = -1; // FLAIR & CSF: darker max
203  else
204  selectionRules[ 2 ][ 0 ] = 1; // PD-w & CSF: brighter max
205 
206  m_SelectedMax.resize( 3, std::vector<double>( 3, 0.0 ) );
207  m_Stds.resize(3, std::vector<double>( 3, 0.0 ));
208 
209  typedef double MeasurementType;
210  const unsigned int MeasurementVectorLength = 1;
211  typedef itk::Vector< MeasurementType , MeasurementVectorLength > MeasurementVectorType;
212  typedef itk::Statistics::ListSample< MeasurementVectorType > ListSampleType;
213  typedef itk::Statistics::DenseFrequencyContainer2 FrequencyContainerType;
214  typedef itk::Statistics::Histogram< MeasurementType,FrequencyContainerType > HistogramType;
215  typedef itk::Statistics::SampleToHistogramFilter<ListSampleType,HistogramType> FilterType;
216 
217  InputImagePixelType valMaxInput = std::numeric_limits<InputImagePixelType>::max();
218  InputImagePixelType valMinInput = std::numeric_limits<InputImagePixelType>::min();
219 
220  // Compute mean values
221  for ( unsigned int im = 1; im < m_NumberOfModalities; im++ )
222  {
223  for ( unsigned int c = 0; c < m_ImagesClasses.size(); c++ )
224  {
225  //create the histogram with the selected mask
226  typename MaskFilterType::Pointer maskFilter = MaskFilterType::New();
227  typename MinMaxCalculatorType::Pointer calculator = MinMaxCalculatorType::New();
228  maskFilter->SetDirectionTolerance(m_Tol);
229  maskFilter->SetCoordinateTolerance(m_Tol);
230  maskFilter->SetInput( m_ImagesVector[ im ] );
231  maskFilter->SetMaskImage( m_ImagesClasses[ c ] );
232 
233  maskFilter->SetOutsideValue( valMaxInput );
234  maskFilter->Update();
235  calculator->SetImage( maskFilter->GetOutput() );
236  calculator->ComputeMinimum();
237  InputImagePixelType minimumResult = calculator->GetMinimum();
238 
239  maskFilter->SetOutsideValue( valMinInput );
240  maskFilter->Update();
241  calculator->SetImage( maskFilter->GetOutput() );
242  calculator->ComputeMaximum();
243  InputImagePixelType maximumResult = calculator->GetMaximum();
244 
245  double minV = static_cast<double>(minimumResult);
246  double maxV = static_cast<double>(maximumResult);
247 
248  unsigned int bins = static_cast<unsigned int>(std::floor( maxV - minV + 1.0 ));
249 
250  ListSampleType::Pointer listSample = ListSampleType::New();
251  listSample->SetMeasurementVectorSize( MeasurementVectorLength );
252 
253  InputConstIteratorType tmpIt (m_ImagesVector[im], m_ImagesVector[im]->GetLargestPossibleRegion() );
254  InputIteratorType classesCIt (m_ImagesClasses[ c ], m_ImagesClasses[ c ]->GetLargestPossibleRegion() );
255 
256  MeasurementVectorType mv;
257  while (!classesCIt.IsAtEnd())
258  {
259  if(classesCIt.Get()!=0)
260  {
261  mv[0] = static_cast <MeasurementType> (tmpIt.Get());
262  listSample->PushBack(mv);
263  }
264  ++classesCIt;
265  ++tmpIt;
266  }
267 
268  // creation histo
269  const unsigned int numberOfComponents = 1;
270  HistogramType::SizeType size( numberOfComponents );
271  size.Fill( bins );
272  HistogramType::MeasurementVectorType min( numberOfComponents );
273  HistogramType::MeasurementVectorType max( numberOfComponents );
274  min.Fill( minV );
275  max.Fill( maxV );
276 
277  FilterType::Pointer filter = FilterType::New();
278  filter->SetInput( listSample );
279  filter->SetHistogramSize( size );
280  filter->SetHistogramBinMinimum( min );
281  filter->SetHistogramBinMaximum( max );
282  filter->Update();
283 
284  HistogramType::ConstPointer histogram = filter->GetOutput();
285 
286  double bins2 = static_cast<double>(bins) / resolution;
287  unsigned int smoothedBins = static_cast<unsigned int>(bins2 + 1);
288 
289  HistogramType::SizeType size2(MeasurementVectorLength);
290  size2.Fill(smoothedBins);
291 
292  HistogramType::MeasurementVectorType lowerBound;
293  lowerBound.SetSize(smoothedBins);
294  lowerBound.Fill(minV);
295 
296  HistogramType::MeasurementVectorType upperBound;
297  upperBound.SetSize(smoothedBins);
298  upperBound.Fill(maxV);
299 
300  HistogramType::Pointer smoothed = HistogramType::New();
301  smoothed->SetMeasurementVectorSize(MeasurementVectorLength);
302  smoothed->Initialize(size2, lowerBound, upperBound );
303  smoothed->SetFrequency(0);
304 
305  std::vector<double> smoothedFrequency(smoothedBins,0);
306 
307  for (unsigned int i = 0; i < bins; ++i)
308  {
309  for (unsigned int j = 0; j < smoothedBins; ++j)
310  {
311  // Calculate distance in intensity values
312  double valueI = static_cast<double>(i) / static_cast<double>(bins -1) * (maxV - minV) + minV;
313  double valueJ = static_cast<double>(j) / static_cast<double>(smoothedBins -1) * (maxV - minV) + minV;
314  double kvalue = (valueI - valueJ) / bandwidth;
315  double norm = std::exp( -0.5 * kvalue * kvalue ) / std::sqrt( 2.0 * M_PI );
316  double value = norm * static_cast<double>(histogram->GetFrequency(i));
317  smoothedFrequency[j] += value;
318  }
319  }
320 
321  for (unsigned int i = 0; i < smoothedBins; i++ )
322  {
323  smoothedFrequency[i]/=(static_cast<double>( histogram->GetTotalFrequency()) * bandwidth);
324  }
325 
326  // get maxs
327  std::vector<unsigned int> maxs;
328  const int localMax = 3;
329  for (unsigned int i = 0; i < smoothedBins; i++ )
330  {
331  int testVal = static_cast <int> (i) - localMax;
332  unsigned int jmin = testVal < 0 ? 0 : static_cast <unsigned int> (testVal);
333 
334  unsigned int jmax = ( i + localMax ) > smoothedBins ? smoothedBins : ( i + localMax );
335 
336  bool ismax = true;
337  for (unsigned int j = jmin; j < jmax; j++ )
338  {
339  if ( (i != j) && (smoothedFrequency[i] <= smoothedFrequency[j]) )
340  {
341  ismax = false;
342  continue;
343  }
344  }
345 
346  if ( ismax && ( smoothedFrequency[i] > 0.001 ) )
347  {
348  maxs.push_back( i );
349  }
350  }
351 
352  if(maxs.size()==0)
353  {
354  std::cerr << "--Error in Hierarchical Initializer: no maximum value" << std::endl;
355  exit(-1);
356  }
357  else
358  {
359  // Keeping only the good max according to the rules
360  double value = 0;
361  switch ( selectionRules[ im ][ c ] )
362  {
363  case - 1: // darker max
364  {
365  value = ( static_cast<double>(maxs[0]) / static_cast<double>(bins-1) ) * (maxV - minV) + minV;
366  m_SelectedMax[im][c] = value;
367  break;
368  }
369  case 1: // brigther max
370  {
371  value = ( static_cast<double>(maxs[ maxs.size() - 1 ]) / static_cast<double>(bins-1) ) * (maxV - minV) + minV;
372  m_SelectedMax[ im ][ c ] = value;
373  break;
374  }
375  case 0: // absolute max
376  default:
377  {
378  unsigned int absMaxIndex = 0;
379  unsigned int absMaxValue = 0;
380  for ( unsigned int i = 0; i < maxs.size(); i++ )
381  {
382  if ( maxs[ i ] > absMaxValue )
383  {
384  absMaxValue = maxs[ i ];
385  absMaxIndex = i;
386  }
387  }
388  value = ( static_cast<double>(maxs[ absMaxIndex ]) / static_cast<double>(bins-1) ) * (maxV - minV) + minV;
389  m_SelectedMax[ im ][ c ] = value;
390  break;
391  }
392  }
393  }
394  }
395  }
396 }
397 
398 template <typename TInputImage, typename TMaskImage>
401 {
402  ImageTypeD::Pointer diff = ImageTypeD::New();
403  diff->SetRegions(this->GetMask()->GetLargestPossibleRegion());
404  diff->CopyInformation(this->GetMask());
405  diff->Allocate();
406  diff->FillBuffer(0);
407  ImageIteratorTypeD diffIt(diff, diff->GetLargestPossibleRegion() );
408 
409  std::vector<InputIteratorType> classesItVec;
410  for(unsigned int i = 0; i < m_ImagesClasses.size(); i++)
411  {
412  InputIteratorType classesIt (m_ImagesClasses[i],m_ImagesClasses[i]->GetLargestPossibleRegion() );
413  classesItVec.push_back(classesIt);
414  }
415 
416  std::vector<InputConstIteratorType> imagesVectorItVec;
417  for(unsigned int i = 0; i < m_NumberOfModalities; i++)
418  {
419  InputConstIteratorType imagesVectorIt (m_ImagesVector[i], m_ImagesVector[i]->GetLargestPossibleRegion() );
420  imagesVectorItVec.push_back(imagesVectorIt);
421  }
422 
423  const double estimatorFactor = 1.4918; // factor to compute standard deviation using a robust variance estimator
424 
425  for ( unsigned int im = 1; im < m_NumberOfModalities; im++ )
426  {
427  for ( unsigned int c = 0; c < m_ImagesClasses.size(); c++ )
428  {
429 
430  for(unsigned int i = 0; i < m_ImagesClasses.size(); i++)
431  {
432  classesItVec[i].GoToBegin();
433  }
434  for(unsigned int i = 0; i < m_NumberOfModalities; i++)
435  {
436  imagesVectorItVec[i].GoToBegin();
437  }
438  diffIt.GoToBegin();
439 
440  while (!diffIt.IsAtEnd())
441  {
442 
443  if(classesItVec[c].Get() > 0)
444  {
445  diffIt.Set(std::fabs(imagesVectorItVec[im].Get()-m_SelectedMax[im][c]));
446  }
447 
448  for(unsigned int i = 0; i < m_ImagesClasses.size(); i++)
449  {
450  ++classesItVec[i];
451  }
452  for(unsigned int i = 0; i < m_NumberOfModalities; i++)
453  {
454  ++imagesVectorItVec[i];
455  }
456  ++diffIt;
457  }
458  m_Stds[im][c] = estimatorFactor * regionMedianValue(diff,m_ImagesClasses[c]);
459  }
460  }
461 }
462 
463 template <typename TInputImage, typename TMaskImage>
466 {
467  m_GaussianModel.clear();
468  m_Alphas.clear();
469 
470  for ( unsigned int i = 0; i < m_NbClasses; i++ )
471  {
472  GaussianFunctionType::Pointer densityFunction1 = GaussianFunctionType::New();
473  densityFunction1->SetMeasurementVectorSize( m_NbClasses );
474  GaussianFunctionType::MeanVectorType mean1( m_NbClasses );
475  GaussianFunctionType::CovarianceMatrixType cov1;
476  cov1.SetSize( m_NbClasses, m_NbClasses );
477  cov1.Fill(0);
478 
479  // Take information from T1-w
480  GaussianFunctionType::MeanVectorType t1mean = (m_Solution1D[i])->GetMean();
481  GaussianFunctionType::CovarianceMatrixType t1var = (m_Solution1D[i])->GetCovariance();
482 
483  mean1[0] = t1mean[0];
484  mean1[1] = m_SelectedMax[1][i];
485  mean1[2] = m_SelectedMax[2][i];
486 
487  cov1[0][0] = t1var[0][0];
488  cov1[1][1] = m_Stds[1][i]*m_Stds[1][i];
489  cov1[2][2] = m_Stds[2][i]*m_Stds[2][i];
490 
491  densityFunction1->SetMean( mean1 );
492  densityFunction1->SetCovariance( cov1 );
493 
494  m_GaussianModel.push_back(densityFunction1);
495  m_Alphas.push_back(m_Solution1DAlphas[i]);
496  }
497 }
498 
499 template <typename TInputImage, typename TMaskImage>
501 ::regionMedianValue(itk::Image< double, 3 >::Pointer image, typename TInputImage::Pointer mask )
502 {
503  ImageIteratorTypeD imageIt(image, image->GetLargestPossibleRegion() );
504  InputIteratorType maskIt(mask, mask->GetLargestPossibleRegion() );
505 
506  std::vector<double> regionValueVector;
507 
508  while(!imageIt.IsAtEnd())
509  {
510  if ( maskIt.Get() != 0 )
511  {
512  regionValueVector.push_back(imageIt.Get());
513  }
514  ++imageIt;
515  ++maskIt;
516  }
517  std::sort (regionValueVector.begin(), regionValueVector.end());
518  return regionValueVector[regionValueVector.size()/2];
519 }
520 
521 template <typename TInputImage, typename TMaskImage>
523 {
524  m_ImagesVector.clear();
525 
526  if(m_IndexImage1<m_NbMaxImages){m_ImagesVector.push_back(this->GetInputImage1());}
527  if(m_IndexImage2<m_NbMaxImages){m_ImagesVector.push_back(this->GetInputImage2());}
528  if(m_IndexImage3<m_NbMaxImages){m_ImagesVector.push_back(this->GetInputImage3());}
529  if(m_IndexImage4<m_NbMaxImages){m_ImagesVector.push_back(this->GetInputImage4());}
530  if(m_IndexImage5<m_NbMaxImages){m_ImagesVector.push_back(this->GetInputImage5());}
531 
532  m_NumberOfModalities = m_ImagesVector.size();
533 
534  this->ComputeSolution1D();
535  this->ComputeInitialT1Classification();
536  this->ComputeHistogram();
537  this->ComputeVariances();
538  this->FillInitialGaussianModel();
539 }
540 
541 
542 } // end of namespace
TInputImage::ConstPointer GetInputImage4()
std::map< double, std::vector< double > > ModelMap2
void SetInputImage5(const TInputImage *image)
TInputImage::ConstPointer GetInputImage2()
itk::ImageRegionIterator< ImageTypeD > ImageIteratorTypeD
TInputImage::ConstPointer GetInputImage5()
void SetInputImage2(const TInputImage *image)
TInputImage::ConstPointer GetInputImage3()
void SetInputImage1(const TInputImage *image)
itk::ImageRegionIterator< InputImageType > InputIteratorType
double regionMedianValue(itk::Image< double, 3 >::Pointer image, typename TInputImage::Pointer mask)
InputImageType::PixelType InputImagePixelType
TInputImage::ConstPointer GetInputImage1()
void eventCallback(itk::Object *caller, const itk::EventObject &event, void *clientData)
itk::VariableLengthVector< double > MeasurementVectorType
void SetInputImage4(const TInputImage *image)
void SetInputImage3(const TInputImage *image)
static Pointer New()
itk::SmartPointer< Self > Pointer
itk::ImageRegionConstIterator< InputImageType > InputConstIteratorType
itk::SmartPointer< Self > Pointer
std::map< double, std::vector< GaussianFunctionType::Pointer > > ModelMap