ANIMA  4.0
animaGcStremMsLesionsSegmentationFilter.hxx
Go to the documentation of this file.
1 #include "animaGcStremMsLesionsSegmentationFilter.h"
2 
3 namespace anima
4 {
5 
6 template <typename TInputImage>
7 void GcStremMsLesionsSegmentationFilter <TInputImage>::SetMask(const TInputImage* mask)
8 {
9  this->SetNthInput(m_NbInputs, const_cast<TInputImage*>(mask));
10  m_IndexMask = m_NbInputs;
11  m_NbInputs++;
12 }
13 
14 template <typename TInputImage>
15 void GcStremMsLesionsSegmentationFilter <TInputImage>::SetInputImageT1(const TInputImage* image)
16 {
17  this->SetNthInput(m_NbInputs, const_cast<TInputImage*>(image));
18  m_IndexImageT1 = m_NbInputs;
19  m_NbInputs++;
20  m_Modalities++;
21 }
22 
23 template <typename TInputImage>
24 void GcStremMsLesionsSegmentationFilter <TInputImage>::SetInputImageT2(const TInputImage* image)
25 {
26  this->SetNthInput(m_NbInputs, const_cast<TInputImage*>(image));
27  m_IndexImageT2 = m_NbInputs;
28  m_NbInputs++;
29  m_Modalities++;
30 }
31 
32 template <typename TInputImage>
33 void GcStremMsLesionsSegmentationFilter <TInputImage>::SetInputImageDP(const TInputImage* image)
34 {
35  this->SetNthInput(m_NbInputs, const_cast<TInputImage*>(image));
36  m_IndexImageDP = m_NbInputs;
37  m_NbInputs++;
38  m_Modalities++;
39 }
40 
41 template <typename TInputImage>
42 void GcStremMsLesionsSegmentationFilter <TInputImage>::SetInputImageFLAIR(const TInputImage* image)
43 {
44  this->SetNthInput(m_NbInputs, const_cast<TInputImage*>(image));
45  m_IndexImageFLAIR = m_NbInputs;
46  m_NbInputs++;
47  m_Modalities++;
48 }
49 
50 template <typename TInputImage>
51 void GcStremMsLesionsSegmentationFilter <TInputImage>::SetInputImageT1Gd(const TInputImage* image)
52 {
53  this->SetNthInput(m_NbInputs, const_cast<TInputImage*>(image));
54  m_IndexImageT1Gd = m_NbInputs;
55  m_NbInputs++;
56  m_Modalities++;
57 }
58 
59 template <typename TInputImage>
60 void GcStremMsLesionsSegmentationFilter <TInputImage>::SetInputCSFAtlas(const ImageTypeD* atlas)
61 {
62  this->SetNthInput(m_NbInputs, const_cast<ImageTypeD*>(atlas));
63  m_IndexAtlasCSF = m_NbInputs;
64  m_NbInputs++;
65 }
66 
67 template <typename TInputImage>
68 void GcStremMsLesionsSegmentationFilter <TInputImage>::SetInputGMAtlas(const ImageTypeD* atlas)
69 {
70  this->SetNthInput(m_NbInputs, const_cast<ImageTypeD*>(atlas));
71  m_IndexAtlasGM = m_NbInputs;
72  m_NbInputs++;
73 }
74 
75 template <typename TInputImage>
76 void GcStremMsLesionsSegmentationFilter <TInputImage>::SetInputWMAtlas(const ImageTypeD* atlas)
77 {
78  this->SetNthInput(m_NbInputs, const_cast<ImageTypeD*>(atlas));
79  m_IndexAtlasWM = m_NbInputs;
80  m_NbInputs++;
81 }
82 
83 template <typename TInputImage>
84 void GcStremMsLesionsSegmentationFilter <TInputImage>::SetInputLesionPrior(const ImageTypeD* image)
85 {
86  this->SetNthInput(m_NbInputs, const_cast<ImageTypeD*>(image));
87  m_IndexLesionPrior = m_NbInputs;
88  m_NbInputs++;
89 }
90 
91 template <typename TInputImage>
92 void GcStremMsLesionsSegmentationFilter <TInputImage>::SetSourcesMask(const ImageTypeUC* image)
93 {
94  this->SetNthInput(m_NbInputs, const_cast<ImageTypeUC*>(image));
95  m_IndexSourcesMask = m_NbInputs;
96  m_NbInputs++;
97 }
98 
99 template <typename TInputImage>
100 void GcStremMsLesionsSegmentationFilter <TInputImage>::SetSinksMask(const ImageTypeUC* image)
101 {
102  this->SetNthInput(m_NbInputs, const_cast<ImageTypeUC*>(image));
103  m_IndexSinksMask = m_NbInputs;
104  m_NbInputs++;
105 }
106 
107 
108 
109 template <typename TInputImage>
110 typename TInputImage::ConstPointer GcStremMsLesionsSegmentationFilter <TInputImage>::GetMask()
111 {
112  return static_cast< const TInputImage * >
113  ( this->itk::ProcessObject::GetInput(m_IndexMask) );
114 }
115 
116 template <typename TInputImage>
117 typename TInputImage::ConstPointer GcStremMsLesionsSegmentationFilter <TInputImage>::GetInputImageT1()
118 {
119  return static_cast< const TInputImage * >
120  ( this->itk::ProcessObject::GetInput(m_IndexImageT1) );
121 }
122 
123 template <typename TInputImage>
124 typename TInputImage::ConstPointer GcStremMsLesionsSegmentationFilter <TInputImage>::GetInputImageT2()
125 {
126  return static_cast< const TInputImage * >
127  ( this->itk::ProcessObject::GetInput(m_IndexImageT2) );
128 }
129 
130 template <typename TInputImage>
131 typename TInputImage::ConstPointer GcStremMsLesionsSegmentationFilter <TInputImage>::GetInputImageDP()
132 {
133  return static_cast< const TInputImage * >
134  ( this->itk::ProcessObject::GetInput(m_IndexImageDP) );
135 }
136 
137 template <typename TInputImage>
138 typename TInputImage::ConstPointer GcStremMsLesionsSegmentationFilter <TInputImage>::GetInputImageFLAIR()
139 {
140  return static_cast< const TInputImage * >
141  ( this->itk::ProcessObject::GetInput(m_IndexImageFLAIR) );
142 }
143 template <typename TInputImage>
144 typename TInputImage::ConstPointer GcStremMsLesionsSegmentationFilter <TInputImage>::GetInputImageT1Gd()
145 {
146  return static_cast< const TInputImage * >
147  ( this->itk::ProcessObject::GetInput(m_IndexImageT1Gd) );
148 }
149 
150 
151 template <typename TInputImage>
152 itk::Image <double,3>::ConstPointer GcStremMsLesionsSegmentationFilter <TInputImage>::GetInputCSFAtlas()
153 {
154  return static_cast< const ImageTypeD * >
155  ( this->itk::ProcessObject::GetInput(m_IndexAtlasCSF) );
156 }
157 
158 template <typename TInputImage>
159 itk::Image <double,3>::ConstPointer GcStremMsLesionsSegmentationFilter <TInputImage>::GetInputGMAtlas()
160 {
161  return static_cast< const ImageTypeD * >
162  ( this->itk::ProcessObject::GetInput(m_IndexAtlasGM) );
163 }
164 
165 template <typename TInputImage>
166 itk::Image <double,3>::ConstPointer GcStremMsLesionsSegmentationFilter <TInputImage>::GetInputWMAtlas()
167 {
168  return static_cast< const ImageTypeD * >
169  ( this->itk::ProcessObject::GetInput(m_IndexAtlasWM) );
170 }
171 
172 template <typename TInputImage>
173 itk::Image <double,3>::ConstPointer GcStremMsLesionsSegmentationFilter <TInputImage>::GetInputLesionPrior()
174 {
175  return static_cast< const ImageTypeD * > (this->itk::ProcessObject::GetInput(m_IndexLesionPrior));
176 }
177 
178 template <typename TInputImage>
179 itk::Image <unsigned char,3>::ConstPointer GcStremMsLesionsSegmentationFilter <TInputImage>::GetSourcesMask()
180 {
181  return static_cast< const ImageTypeUC * >
182  ( this->itk::ProcessObject::GetInput(m_IndexSourcesMask) );
183 }
184 template <typename TInputImage>
185 itk::Image <unsigned char,3>::ConstPointer GcStremMsLesionsSegmentationFilter <TInputImage>::GetSinksMask()
186 {
187  return static_cast< const ImageTypeUC * >
188  ( this->itk::ProcessObject::GetInput(m_IndexSinksMask) );
189 }
190 
191 
192 template <typename TInputImage>
193 itk::DataObject::Pointer GcStremMsLesionsSegmentationFilter <TInputImage>::MakeOutput(unsigned int idx)
194 {
195  itk::DataObject::Pointer output;
196 
197  if( (idx < 11) || (idx == 17) )
198  {
199  output = ( TOutputImage::New() ).GetPointer();
200  }
201  else if(( 11 <= idx ) && ( idx < 17 ))
202  {
203  output = ( ImageTypeD::New() ).GetPointer();
204  }
205  else
206  {
207  std::cerr << "No output " << idx << std::endl;
208  output = NULL;
209  }
210  return output.GetPointer();
211 }
212 
213 template <typename TInputImage>
214 itk::Image <unsigned char, 3>* GcStremMsLesionsSegmentationFilter <TInputImage>::GetOutputLesions()
215 {
216  return dynamic_cast< TOutputImage* >( this->itk::ProcessObject::GetOutput(0) );
217 }
218 template <typename TInputImage>
219 itk::Image <unsigned char, 3>* GcStremMsLesionsSegmentationFilter <TInputImage>::GetOutputCSF()
220 {
221  return dynamic_cast< TOutputImage* >( this->itk::ProcessObject::GetOutput(1) );
222 }
223 template <typename TInputImage>
224 itk::Image <unsigned char, 3>* GcStremMsLesionsSegmentationFilter <TInputImage>::GetOutputGM()
225 {
226  return dynamic_cast< TOutputImage* >( this->itk::ProcessObject::GetOutput(2) );
227 }
228 template <typename TInputImage>
229 itk::Image <unsigned char, 3>* GcStremMsLesionsSegmentationFilter <TInputImage>::GetOutputWM()
230 {
231  return dynamic_cast< TOutputImage* >( this->itk::ProcessObject::GetOutput(3) );
232 }
233 template <typename TInputImage>
234 itk::Image <unsigned char, 3>* GcStremMsLesionsSegmentationFilter <TInputImage>::GetOutputWholeSeg()
235 {
236  return dynamic_cast< TOutputImage* >( this->itk::ProcessObject::GetOutput(4) );
237 }
238 template <typename TInputImage>
239 itk::Image <unsigned char, 3>* GcStremMsLesionsSegmentationFilter <TInputImage>::GetOutputStrem()
240 {
241  return dynamic_cast< TOutputImage* >( this->itk::ProcessObject::GetOutput(5) );
242 }
243 template <typename TInputImage>
244 itk::Image <unsigned char, 3>* GcStremMsLesionsSegmentationFilter <TInputImage>::GetOutputStremCSF()
245 {
246  return dynamic_cast< TOutputImage* >( this->itk::ProcessObject::GetOutput(6) );
247 }
248 template <typename TInputImage>
249 itk::Image <unsigned char, 3>* GcStremMsLesionsSegmentationFilter <TInputImage>::GetOutputStremGM()
250 {
251  return dynamic_cast< TOutputImage* >( this->itk::ProcessObject::GetOutput(7) );
252 }
253 template <typename TInputImage>
254 itk::Image <unsigned char, 3>* GcStremMsLesionsSegmentationFilter <TInputImage>::GetOutputStremWM()
255 {
256  return dynamic_cast< TOutputImage* >( this->itk::ProcessObject::GetOutput(8) );
257 }
258 template <typename TInputImage>
259 itk::Image <unsigned char, 3>* GcStremMsLesionsSegmentationFilter <TInputImage>::GetOutputIntensityImage1()
260 {
261  return dynamic_cast< TOutputImage* >( this->itk::ProcessObject::GetOutput(9) );
262 }
263 template <typename TInputImage>
264 itk::Image <unsigned char, 3>* GcStremMsLesionsSegmentationFilter <TInputImage>::GetOutputIntensityImage2()
265 {
266  return dynamic_cast< TOutputImage* >( this->itk::ProcessObject::GetOutput(10) );
267 }
268 template <typename TInputImage>
269 itk::Image <double,3>* GcStremMsLesionsSegmentationFilter <TInputImage>::GetOutputFuzzyObjectImage()
270 {
271  return dynamic_cast< ImageTypeD* >( this->itk::ProcessObject::GetOutput(11) );
272 }
273 template <typename TInputImage>
274 itk::Image <double,3>* GcStremMsLesionsSegmentationFilter <TInputImage>::GetOutputMahaCSFImage()
275 {
276  return dynamic_cast< ImageTypeD* >( this->itk::ProcessObject::GetOutput(12) );
277 }
278 template <typename TInputImage>
279 itk::Image <double,3>* GcStremMsLesionsSegmentationFilter <TInputImage>::GetOutputMahaGMImage()
280 {
281  return dynamic_cast< ImageTypeD* >( this->itk::ProcessObject::GetOutput(13) );
282 }
283 template <typename TInputImage>
284 itk::Image <double,3>* GcStremMsLesionsSegmentationFilter <TInputImage>::GetOutputMahaWMImage()
285 {
286  return dynamic_cast< ImageTypeD* >( this->itk::ProcessObject::GetOutput(14) );
287 }
288 template <typename TInputImage>
289 itk::Image <double,3>* GcStremMsLesionsSegmentationFilter <TInputImage>::GetOutputMahaMinimumImage()
290 {
291  return dynamic_cast< ImageTypeD* >( this->itk::ProcessObject::GetOutput(15) );
292 }
293 template <typename TInputImage>
294 itk::Image <double,3>* GcStremMsLesionsSegmentationFilter <TInputImage>::GetOutputMahaMaximumImage()
295 {
296  return dynamic_cast< ImageTypeD* >( this->itk::ProcessObject::GetOutput(16) );
297 }
298 template <typename TInputImage>
299 itk::Image <unsigned char, 3>* GcStremMsLesionsSegmentationFilter <TInputImage>::GetOutputGraphCut()
300 {
301  return dynamic_cast< TOutputImage* >( this->itk::ProcessObject::GetOutput(17) );
302 }
303 
304 
305 template <typename TInputImage>
306 void
307 GcStremMsLesionsSegmentationFilter <TInputImage>::WriteOutputs()
308 {
309  if( m_OutputLesionFilename != "" )
310  {
311  std::cout << "Writing Lesion output image to: " << m_OutputLesionFilename << std::endl;
312  anima::writeImage<TOutputImage>(m_OutputLesionFilename, this->GetOutputLesions());
313  }
314 
315  if( m_LesionSegmentationType!=manualGC )
316  {
317  m_MahalanobisFilter->WriteOutputs();
318  if( m_OutputIntensityImage1Filename != "" )
319  {
320  std::cout << "Writing intensity output image 1 to: " << m_OutputIntensityImage1Filename << std::endl;
321  anima::writeImage<TOutputImage>(m_OutputIntensityImage1Filename, this->GetOutputIntensityImage1());
322  }
323  if( m_OutputIntensityImage2Filename != "" )
324  {
325  std::cout << "Writing intensity output image 2 to: " << m_OutputIntensityImage2Filename << std::endl;
326  anima::writeImage<TOutputImage>(m_OutputIntensityImage2Filename, this->GetOutputIntensityImage2());
327  }
328  if( m_OutputCSFFilename != "" )
329  {
330  std::cout << "Writing CSF output image to: " << m_OutputCSFFilename << std::endl;
331  anima::writeImage<TOutputImage>(m_OutputCSFFilename, this->GetOutputCSF());
332  }
333  if( m_OutputGMFilename != "" )
334  {
335  std::cout << "Writing GM output image to: " << m_OutputGMFilename << std::endl;
336  anima::writeImage<TOutputImage>(m_OutputGMFilename, this->GetOutputGM());
337  }
338  if( m_OutputWMFilename != "" )
339  {
340  std::cout << "Writing WM output image to: " << m_OutputWMFilename << std::endl;
341  anima::writeImage<TOutputImage>(m_OutputWMFilename, this->GetOutputWM());
342  }
343  if( m_OutputWholeFilename != "" )
344  {
345  std::cout << "Writing segmentation output image to: " << m_OutputWholeFilename << std::endl;
346  anima::writeImage<TOutputImage>(m_OutputWholeFilename, this->GetOutputWholeSeg());
347  }
348  }
349 
350  if( m_LesionSegmentationType==gcem || m_LesionSegmentationType==gcemAndManualGC )
351  {
352  if( m_OutputFuzzyObjectFilename != "" )
353  {
354  std::cout << "Writing fuzzy object output image to: " << m_OutputFuzzyObjectFilename << std::endl;
355  anima::writeImage<ImageTypeD>(m_OutputFuzzyObjectFilename, this->GetOutputFuzzyObjectImage());
356  }
357  }
358 
359  if( m_LesionSegmentationType==strem )
360  {
361  if( m_OutputStremFilename != "" )
362  {
363  std::cout << "Writing strem output image to: " << m_OutputStremFilename << std::endl;
364  anima::writeImage<TOutputImage>(m_OutputStremFilename, this->GetOutputStrem());
365  }
366  if( m_OutputStremCSFFilename != "" )
367  {
368  std::cout << "Writing CSF strem output image to: " << m_OutputStremCSFFilename << std::endl;
369  anima::writeImage<TOutputImage>(m_OutputStremCSFFilename, this->GetOutputStremCSF());
370  }
371  if( m_OutputStremGMFilename != "" )
372  {
373  std::cout << "Writing GM strem output image to: " << m_OutputStremGMFilename << std::endl;
374  anima::writeImage<TOutputImage>(m_OutputStremGMFilename, this->GetOutputStremGM());
375  }
376  if( m_OutputStremWMFilename != "" )
377  {
378  std::cout << "Writing WM strem output image to: " << m_OutputStremWMFilename << std::endl;
379  anima::writeImage<TOutputImage>(m_OutputStremWMFilename, this->GetOutputStremWM());
380  }
381  }
382  else
383  {
384  m_GraphCutFilter->WriteOutputs();
385  }
386 }
387 
388 
389 template <typename TInputImage>
390 void
391 GcStremMsLesionsSegmentationFilter <TInputImage>::CheckInputImages()
392 {
393  std::cout << "Check input entries..." << std::endl;
394  if(m_LesionSegmentationType!=manualGC)
395  {
396  if(m_Modalities < 3)
397  itkExceptionMacro("-- Error in anima::Gc_Strem_MS_Lesions_Segmentation: automatic segmentation requires 3 modalities, exiting...");
398 
399  if(this->GetInputImageT1().IsNull())
400  itkExceptionMacro("-- Error in anima::Gc_Strem_MS_Lesions_Segmentation: automatic segmentation requires T1 image, exiting...");
401 
402  if(m_Modalities == 3)
403  {
404  m_UseT2 = false;
405  m_UseDP = false;
406  m_UseFLAIR = false;
407 
408  if(this->GetInputImageT2().IsNotNull()){m_UseT2=true;}
409  if(this->GetInputImageDP().IsNotNull()){m_UseDP=true;}
410  if(this->GetInputImageFLAIR().IsNotNull()){m_UseFLAIR=true;}
411  }
412 
413  if(m_Modalities>3)
414  {
415  if( !( (m_UseT2 & m_UseDP & !m_UseFLAIR) || (m_UseT2 & !m_UseDP & m_UseFLAIR) || (!m_UseT2 & m_UseDP & m_UseFLAIR) ) )
416  itkExceptionMacro("-- Error in anima::Gc_Strem_MS_Lesions_Segmentation: 2 images among T2, DP and FLAIR must be choosen for automatic segmentation, exiting...");
417  }
418  }
419  else
420  {
421  if(this->GetSourcesMask().IsNull() || this->GetSinksMask().IsNull())
422  itkExceptionMacro("-- Error in anima::Gc_Strem_MS_Lesions_Segmentation: seed mask missing, exiting...");
423  }
424 }
425 template <typename TInputImage>
426 void
427 GcStremMsLesionsSegmentationFilter <TInputImage>::RescaleImages()
428 {
429  std::cout << "Rescale input images..." << std::endl;
430 
431  ImageTypeUC::Pointer InputImage_T2_UC = ImageTypeUC::New();
432  ImageTypeUC::Pointer InputImage_DP_UC = ImageTypeUC::New();
433  ImageTypeUC::Pointer InputImage_FLAIR_UC = ImageTypeUC::New();
434 
435  double desiredMinimum=0,desiredMaximum=255;
436 
437  if(this->GetInputImageT1().IsNotNull())
438  {
439  typename RescaleFilterType::Pointer rescaleFilter1 = RescaleFilterType::New();
440  rescaleFilter1->SetInput( this->GetInputImageT1() );
441  rescaleFilter1->SetOutputMinimum( desiredMinimum );
442  rescaleFilter1->SetOutputMaximum( desiredMaximum );
443  rescaleFilter1->SetNumberOfWorkUnits( this->GetNumberOfWorkUnits() );
444  rescaleFilter1->SetCoordinateTolerance( m_Tol );
445  rescaleFilter1->SetDirectionTolerance( m_Tol );
446  rescaleFilter1->UpdateLargestPossibleRegion();
447  m_InputImage_T1_UC = rescaleFilter1->GetOutput();
448  }
449  if(this->GetInputImageT2().IsNotNull())
450  {
451  typename RescaleFilterType::Pointer rescaleFilter2 = RescaleFilterType::New();
452  rescaleFilter2->SetInput( this->GetInputImageT2() );
453  rescaleFilter2->SetOutputMinimum( desiredMinimum );
454  rescaleFilter2->SetOutputMaximum( desiredMaximum );
455  rescaleFilter2->SetNumberOfWorkUnits( this->GetNumberOfWorkUnits() );
456  rescaleFilter2->SetCoordinateTolerance( m_Tol );
457  rescaleFilter2->SetDirectionTolerance( m_Tol );
458  rescaleFilter2->UpdateLargestPossibleRegion();
459  InputImage_T2_UC = rescaleFilter2->GetOutput();
460  }
461  if(this->GetInputImageFLAIR().IsNotNull())
462  {
463  typename RescaleFilterType::Pointer rescaleFilter3 = RescaleFilterType::New();
464  rescaleFilter3->SetInput( this->GetInputImageFLAIR() );
465  rescaleFilter3->SetOutputMinimum( desiredMinimum );
466  rescaleFilter3->SetOutputMaximum( desiredMaximum );
467  rescaleFilter3->SetNumberOfWorkUnits( this->GetNumberOfWorkUnits() );
468  rescaleFilter3->SetCoordinateTolerance( m_Tol );
469  rescaleFilter3->SetDirectionTolerance( m_Tol );
470  rescaleFilter3->UpdateLargestPossibleRegion();
471  InputImage_FLAIR_UC = rescaleFilter3->GetOutput();
472  }
473  if(this->GetInputImageDP().IsNotNull())
474  {
475  typename RescaleFilterType::Pointer rescaleFilter4 = RescaleFilterType::New();
476  rescaleFilter4->SetInput( this->GetInputImageDP() );
477  rescaleFilter4->SetOutputMinimum( desiredMinimum );
478  rescaleFilter4->SetOutputMaximum( desiredMaximum );
479  rescaleFilter4->SetNumberOfWorkUnits( this->GetNumberOfWorkUnits() );
480  rescaleFilter4->SetCoordinateTolerance( m_Tol );
481  rescaleFilter4->SetDirectionTolerance( m_Tol );
482  rescaleFilter4->UpdateLargestPossibleRegion();
483  InputImage_DP_UC = rescaleFilter4->GetOutput();
484  }
485 
486  // --------------- Create Brain Mask --------------- //
487  m_MaskUC = ImageTypeUC::New();
488  m_MaskUC->SetRegions( this->GetMask()->GetLargestPossibleRegion() );
489  m_MaskUC->CopyInformation( this->GetMask() );
490  m_MaskUC->Allocate();
491  m_MaskUC->FillBuffer(0);
492 
493  ImageIteratorTypeUC itMaskUC( m_MaskUC, m_MaskUC->GetLargestPossibleRegion() );
494  ImageIteratorTypeConstInput itInputImage_Mask( this->GetMask(), this->GetMask()->GetLargestPossibleRegion() );
495 
496  while (!itMaskUC.IsAtEnd())
497  {
498  if(itInputImage_Mask.Get()!=0)
499  {
500  itMaskUC.Set(1);
501  }
502  ++itMaskUC;
503  ++itInputImage_Mask;
504  }
505 
506  if(m_LesionSegmentationType!=manualGC)
507  {
508  if(m_UseT2 && m_UseDP)
509  {
510  m_InputImage_1_UC = InputImage_T2_UC;
511  m_InputImage_2_UC = InputImage_DP_UC;
512  }
513  if(m_UseT2 && m_UseFLAIR)
514  {
515  m_InputImage_1_UC = InputImage_T2_UC;
516  m_InputImage_2_UC = InputImage_FLAIR_UC;
517  }
518  if(m_UseDP && m_UseFLAIR)
519  {
520  m_InputImage_1_UC = InputImage_DP_UC;
521  m_InputImage_2_UC = InputImage_FLAIR_UC;
522 
523  if(m_InitMethodType==1)
524  {
525  m_InputImage_1_UC = InputImage_FLAIR_UC;
526  m_InputImage_2_UC = InputImage_DP_UC;
527  }
528  }
529  }
530 }
531 
532 template <typename TInputImage>
533 void
534 GcStremMsLesionsSegmentationFilter <TInputImage>::ComputeAutomaticInitialization()
535 {
536  ComputeSolutionType::Pointer ComputeSolutionProcess = ComputeSolutionType::New();
537 
538  ComputeSolutionProcess->SetMask( m_MaskUC );
539  ComputeSolutionProcess->SetInputImage1( m_InputImage_T1_UC );
540  ComputeSolutionProcess->SetInputImage2( m_InputImage_1_UC );
541  ComputeSolutionProcess->SetInputImage3( m_InputImage_2_UC );
542 
543  ComputeSolutionProcess->SetInputCSFAtlas( this->GetInputCSFAtlas() );
544  ComputeSolutionProcess->SetInputGMAtlas( this->GetInputGMAtlas() );
545  ComputeSolutionProcess->SetInputWMAtlas( this->GetInputWMAtlas() );
546 
547  ComputeSolutionProcess->SetSolutionWriteFilename( m_SolutionWriteFilename );
548  ComputeSolutionProcess->SetSolutionReadFilename( m_SolutionReadFilename );
549  ComputeSolutionProcess->SetGaussianModel( m_GaussianModel );
550  ComputeSolutionProcess->SetAlphas( m_Alphas );
551 
552  ComputeSolutionProcess->SetUseT2( m_UseT2 );
553  ComputeSolutionProcess->SetUseDP( m_UseDP );
554  ComputeSolutionProcess->SetUseFLAIR( m_UseFLAIR );
555 
556  ComputeSolutionProcess->SetInitMethodType( m_InitMethodType );
557  ComputeSolutionProcess->SetRejRatioHierar( m_RejRatioHierar );
558  ComputeSolutionProcess->SetMinDistance( m_MinDistance );
559  ComputeSolutionProcess->SetEmIter( m_EmIter );
560  ComputeSolutionProcess->SetRejRatio( m_RejRatio );
561  ComputeSolutionProcess->SetEmIter_concentration( m_EmIter_concentration );
562  ComputeSolutionProcess->SetEM_before_concentration( m_EM_before_concentration );
563 
564  ComputeSolutionProcess->SetVerbose( m_Verbose );
565 
566  try
567  {
568  ComputeSolutionProcess->Update();
569  }
570  catch (itk::ExceptionObject &e)
571  {
572  std::cerr << e << std::endl;
573  exit(-1);
574  }
575 
576  m_GaussianModel = ComputeSolutionProcess->GetGaussianModel();
577 
578  std::cout << "Computing mahalanobis images..." << std::endl;
579 
580  m_MahalanobisFilter->SetGaussianModel( ComputeSolutionProcess->GetGaussianModel() );
581  m_MahalanobisFilter->SetMask( m_MaskUC );
582  m_MahalanobisFilter->SetInputImage1( m_InputImage_T1_UC );
583  m_MahalanobisFilter->SetInputImage2( m_InputImage_1_UC );
584  m_MahalanobisFilter->SetInputImage3( m_InputImage_2_UC );
585  m_MahalanobisFilter->SetTol( m_Tol );
586  m_MahalanobisFilter->SetNumberOfWorkUnits( this->GetNumberOfWorkUnits() );
587 
588  m_MahalanobisFilter->SetOutputMahaCSFFilename( m_OutputMahaCSFFilename );
589  m_MahalanobisFilter->SetOutputMahaGMFilename( m_OutputMahaGMFilename );
590  m_MahalanobisFilter->SetOutputMahaWMFilename( m_OutputMahaWMFilename );
591  m_MahalanobisFilter->SetOutputMahaMaximumFilename( m_OutputMahaMaximumFilename );
592  m_MahalanobisFilter->SetOutputMahaMinimumFilename ( m_OutputMahaMinimumFilename );
593 
594  m_MahalanobisFilter->GraftNthOutput( 0, this->GetOutputMahaCSFImage() );
595  m_MahalanobisFilter->GraftNthOutput( 1, this->GetOutputMahaGMImage() );
596  m_MahalanobisFilter->GraftNthOutput( 2, this->GetOutputMahaWMImage() );
597  m_MahalanobisFilter->GraftNthOutput( 3, this->GetOutputMahaMinimumImage() );
598  m_MahalanobisFilter->GraftNthOutput( 4, this->GetOutputMahaMaximumImage() );
599  m_MahalanobisFilter->Update();
600  this->GraftNthOutput( 12 , m_MahalanobisFilter->GetOutputMahaCSF() );
601  this->GraftNthOutput( 13 , m_MahalanobisFilter->GetOutputMahaGM() );
602  this->GraftNthOutput( 14 , m_MahalanobisFilter->GetOutputMahaWM() );
603  this->GraftNthOutput( 15 , m_MahalanobisFilter->GetOutputMahaMinimum() );
604  this->GraftNthOutput( 16 , m_MahalanobisFilter->GetOutputMahaMaximum() );
605 
606 }
607 
608 template <typename TInputImage>
609 void
610 GcStremMsLesionsSegmentationFilter <TInputImage>::StremThreshold()
611 {
612  std::cout << "Compute strem threshold images..." << std::endl;
613 
614  const unsigned char insideValue = 1;
615  const unsigned char outsideValue = 0;
616 
617  BinaryThresholdImageFilterType_F_UC::Pointer thresholdFilterCSF = BinaryThresholdImageFilterType_F_UC::New();
618  thresholdFilterCSF->SetInput( m_MahalanobisFilter->GetOutputMahaCSF() );
619  thresholdFilterCSF->SetUpperThreshold( m_MahalanobisThCSF );
620  thresholdFilterCSF->SetInsideValue( insideValue );
621  thresholdFilterCSF->SetOutsideValue( outsideValue );
622  thresholdFilterCSF->SetNumberOfWorkUnits( this->GetNumberOfWorkUnits() );
623  thresholdFilterCSF->SetCoordinateTolerance( m_Tol );
624  thresholdFilterCSF->SetDirectionTolerance( m_Tol );
625 
626  MaskFilterType_UC_UC::Pointer maskFilterCSF = MaskFilterType_UC_UC::New();
627  maskFilterCSF->SetInput( thresholdFilterCSF->GetOutput()) ;
628  maskFilterCSF->SetMaskImage( m_MaskUC );
629  maskFilterCSF->SetNumberOfWorkUnits( this->GetNumberOfWorkUnits() );
630  maskFilterCSF->SetCoordinateTolerance( m_Tol );
631  maskFilterCSF->SetDirectionTolerance( m_Tol );
632  maskFilterCSF->GraftOutput( this->GetOutputStremCSF() );
633  maskFilterCSF->Update();
634  this->GraftNthOutput( 6 , maskFilterCSF->GetOutput() );
635 
636  BinaryThresholdImageFilterType_F_UC::Pointer thresholdFilterGM = BinaryThresholdImageFilterType_F_UC::New();
637  thresholdFilterGM->SetInput( m_MahalanobisFilter->GetOutputMahaGM() );
638  thresholdFilterGM->SetUpperThreshold( m_MahalanobisThGM );
639  thresholdFilterGM->SetInsideValue( insideValue );
640  thresholdFilterGM->SetOutsideValue( outsideValue );
641  thresholdFilterGM->SetNumberOfWorkUnits( this->GetNumberOfWorkUnits() );
642  thresholdFilterGM->SetCoordinateTolerance( m_Tol );
643  thresholdFilterGM->SetDirectionTolerance( m_Tol );
644 
645  MaskFilterType_UC_UC::Pointer maskFilterGM = MaskFilterType_UC_UC::New();
646  maskFilterGM->SetInput( thresholdFilterGM->GetOutput() );
647  maskFilterGM->SetMaskImage( m_MaskUC );
648  maskFilterGM->SetNumberOfWorkUnits( this->GetNumberOfWorkUnits() );
649  maskFilterGM->SetCoordinateTolerance( m_Tol );
650  maskFilterGM->SetDirectionTolerance( m_Tol );
651  maskFilterGM->GraftOutput( this->GetOutputStremGM() );
652  maskFilterGM->Update();
653  this->GraftNthOutput( 7 , maskFilterGM->GetOutput() );
654 
655  BinaryThresholdImageFilterType_F_UC::Pointer thresholdFilterWM = BinaryThresholdImageFilterType_F_UC::New();
656  thresholdFilterWM->SetInput( m_MahalanobisFilter->GetOutputMahaWM() );
657  thresholdFilterWM->SetUpperThreshold( m_MahalanobisThWM );
658  thresholdFilterWM->SetInsideValue( insideValue );
659  thresholdFilterWM->SetOutsideValue( outsideValue );
660  thresholdFilterWM->SetNumberOfWorkUnits( this->GetNumberOfWorkUnits() );
661  thresholdFilterWM->SetCoordinateTolerance( m_Tol );
662  thresholdFilterWM->SetDirectionTolerance( m_Tol );
663 
664  MaskFilterType_UC_UC::Pointer maskFilterWM = MaskFilterType_UC_UC::New();
665  maskFilterWM->SetInput( thresholdFilterWM->GetOutput() );
666  maskFilterWM->SetMaskImage( m_MaskUC );
667  maskFilterWM->SetNumberOfWorkUnits( this->GetNumberOfWorkUnits() );
668  maskFilterWM->SetCoordinateTolerance( m_Tol );
669  maskFilterWM->SetDirectionTolerance( m_Tol );
670  maskFilterWM->GraftOutput( this->GetOutputStremWM() );
671  maskFilterWM->Update();
672  this->GraftNthOutput( 8 , maskFilterWM->GetOutput() );
673 
674  MinimumFilterTypeUC::Pointer filtermin1 = MinimumFilterTypeUC::New();
675  MinimumFilterTypeUC::Pointer filtermin2 = MinimumFilterTypeUC::New();
676 
677  filtermin1->SetInput1( maskFilterCSF->GetOutput() );
678  filtermin1->SetInput2( maskFilterGM->GetOutput() );
679  filtermin1->SetNumberOfWorkUnits( this->GetNumberOfWorkUnits() );
680  filtermin1->SetCoordinateTolerance( m_Tol);
681  filtermin1->SetDirectionTolerance( m_Tol);
682 
683  filtermin2->SetInput1( filtermin1->GetOutput() );
684  filtermin2->SetInput2( maskFilterWM->GetOutput() );
685  filtermin2->SetNumberOfWorkUnits( this->GetNumberOfWorkUnits() );
686  filtermin2->SetCoordinateTolerance( m_Tol );
687  filtermin2->SetDirectionTolerance( m_Tol );
688  filtermin2->GraftOutput( this->GetOutputStrem() );
689  filtermin2->Update();
690  this->GraftNthOutput( 5 , filtermin2->GetOutput() );
691 }
692 
693 template <typename TInputImage>
694 void
695 GcStremMsLesionsSegmentationFilter <TInputImage>::CreateFuzzyRuleImage()
696 {
697  std::cout << "Computing fuzzy rules..." << std::endl;
698 
699  MinimumFilterTypeF::Pointer filtermin1 = MinimumFilterTypeF::New();
700  MinimumFilterTypeF::Pointer filtermin2 = MinimumFilterTypeF::New();
701 
702  const double minimumOutputValue = 0;
703  const double maximumOutputValue = 1;
704 
705  const unsigned int indexWM = 2;
706  GaussianFunctionType::MeanVectorType mean = (m_GaussianModel[indexWM])->GetMean();
707  GaussianFunctionType::CovarianceMatrixType covar = (m_GaussianModel[indexWM])->GetCovariance();
708 
709  const unsigned int indexImage1 = 1;
710  unsigned char intensityWindowLowerInputValue1 = static_cast<unsigned char>( mean[indexImage1] + m_FuzzyRuleMin * std::sqrt(covar[indexImage1][indexImage1]) );
711  unsigned char intensityWindowUpperInputValue1 = static_cast<unsigned char>( mean[indexImage1] + m_FuzzyRuleMax * std::sqrt(covar[indexImage1][indexImage1]) );
712 
713  IntensityWindowingImageFilterType::Pointer intensityWindowFilter1 = IntensityWindowingImageFilterType::New();
714  intensityWindowFilter1->SetInput( m_InputImage_1_UC );
715  intensityWindowFilter1->SetWindowMinimum( intensityWindowLowerInputValue1 );
716  intensityWindowFilter1->SetWindowMaximum( intensityWindowUpperInputValue1 );
717  intensityWindowFilter1->SetOutputMinimum( minimumOutputValue );
718  intensityWindowFilter1->SetOutputMaximum( maximumOutputValue );
719  intensityWindowFilter1->SetNumberOfWorkUnits( this->GetNumberOfWorkUnits() );
720  intensityWindowFilter1->SetCoordinateTolerance( m_Tol );
721  intensityWindowFilter1->SetDirectionTolerance( m_Tol );
722 
723  const unsigned int indexImage2 = 2;
724  unsigned char intensityWindowLowerInputValue2 = static_cast<unsigned char>( mean[indexImage2] + m_FuzzyRuleMin * std::sqrt(covar[indexImage2][indexImage2]) );
725  unsigned char intensityWindowUpperInputValue2 = static_cast<unsigned char>( mean[indexImage2] + m_FuzzyRuleMax * std::sqrt(covar[indexImage2][indexImage2]) );
726 
727  IntensityWindowingImageFilterType::Pointer intensityWindowFilter2 = IntensityWindowingImageFilterType::New();
728  intensityWindowFilter2->SetInput( m_InputImage_2_UC );
729  intensityWindowFilter2->SetWindowMinimum( intensityWindowLowerInputValue2 );
730  intensityWindowFilter2->SetWindowMaximum( intensityWindowUpperInputValue2 );
731  intensityWindowFilter2->SetOutputMinimum( minimumOutputValue );
732  intensityWindowFilter2->SetOutputMaximum( maximumOutputValue );
733  intensityWindowFilter2->SetNumberOfWorkUnits( this->GetNumberOfWorkUnits() );
734  intensityWindowFilter2->SetCoordinateTolerance( m_Tol );
735  intensityWindowFilter2->SetDirectionTolerance( m_Tol );
736 
737  filtermin1->SetInput1( intensityWindowFilter1->GetOutput() );
738  filtermin1->SetInput2( intensityWindowFilter2->GetOutput() );
739  filtermin1->SetNumberOfWorkUnits( this->GetNumberOfWorkUnits() );
740  filtermin1->SetCoordinateTolerance( m_Tol );
741  filtermin1->SetDirectionTolerance( m_Tol );
742 
743  if (this->GetInputLesionPrior().IsNotNull())
744  {
745  typedef itk::ImageRegionIterator <ImageTypeD> IteratorType;
746  IteratorType mahaMinimumIterator(m_MahalanobisFilter->GetOutputMahaMinimum(),m_MahalanobisFilter->GetOutputMahaMinimum()->GetLargestPossibleRegion());
747  typedef itk::ImageRegionConstIterator <ImageTypeD> ConstIteratorType;
748  ConstIteratorType lesionPriorIterator(this->GetInputLesionPrior(),m_MahalanobisFilter->GetOutputMahaMinimum()->GetLargestPossibleRegion());
749 
750  while (!lesionPriorIterator.IsAtEnd())
751  {
752  double mahaValue = mahaMinimumIterator.Get();
753  double lesionPriorValue = lesionPriorIterator.Get();
754 
755  mahaMinimumIterator.Set(mahaValue * (1.0 - m_LesionPriorProportion) + lesionPriorValue * m_LesionPriorProportion);
756 
757  ++mahaMinimumIterator;
758  ++lesionPriorIterator;
759  }
760  }
761 
762  filtermin2->SetInput1( filtermin1->GetOutput() );
763  filtermin2->SetInput2(m_MahalanobisFilter->GetOutputMahaMinimum());
764  filtermin2->SetNumberOfWorkUnits( this->GetNumberOfWorkUnits() );
765  filtermin2->SetCoordinateTolerance( m_Tol );
766  filtermin2->SetDirectionTolerance( m_Tol );
767 
768  MaskFilterType_F_UC::Pointer maskFilter = MaskFilterType_F_UC::New();
769  maskFilter->SetInput( filtermin2->GetOutput() ) ;
770  maskFilter->SetMaskImage( m_MaskUC );
771  maskFilter->SetNumberOfWorkUnits( this->GetNumberOfWorkUnits() );
772  maskFilter->SetCoordinateTolerance( m_Tol );
773  maskFilter->SetDirectionTolerance( m_Tol );
774  maskFilter->GraftOutput( this->GetOutputFuzzyObjectImage() );
775  maskFilter->Update();
776  this->GraftNthOutput( 11 , maskFilter->GetOutput() );
777 
778  m_FuzzyObject = ImageTypeD::New();
779  m_FuzzyObject->Graft(maskFilter->GetOutput());
780 }
781 
782 template <typename TInputImage>
783 void
784 GcStremMsLesionsSegmentationFilter <TInputImage>::GraphCut()
785 {
786  m_GraphCutFilter->SetUseSpectralGradient( m_UseSpecGrad );
787  m_GraphCutFilter->SetAlpha( m_Alpha );
788  m_GraphCutFilter->SetSigma( m_Sigma );
789  m_GraphCutFilter->SetMultiVarSources( m_MultiVarSources );
790  m_GraphCutFilter->SetMultiVarSinks( m_MultiVarSinks );
791  m_GraphCutFilter->SetMatrixGradFilename( m_MatrixGradFilename );
792  m_GraphCutFilter->SetOutputFilename( m_OutputGCFilename );
793  m_GraphCutFilter->SetVerbose( m_Verbose );
794  m_GraphCutFilter->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits());
795  m_GraphCutFilter->SetTol( m_Tol );
796 
797  unsigned int index = 0;
798  if(this->GetInputImageT1().IsNotNull()){ m_GraphCutFilter->SetInputImage(index, this->GetInputImageT1() ); ++index;}
799  if(this->GetInputImageT2().IsNotNull()){ m_GraphCutFilter->SetInputImage(index, this->GetInputImageT2() ); ++index;}
800  if(this->GetInputImageDP().IsNotNull()){ m_GraphCutFilter->SetInputImage(index, this->GetInputImageDP() ); ++index;}
801  if(this->GetInputImageFLAIR().IsNotNull()){ m_GraphCutFilter->SetInputImage(index, this->GetInputImageFLAIR() ); ++index;}
802  if(this->GetInputImageT1Gd().IsNotNull()){ m_GraphCutFilter->SetInputImage(index, this->GetInputImageT1Gd() ); ++index;}
803 
804  m_GraphCutFilter->SetMask( m_MaskUC );
805 
806  if( m_LesionSegmentationType==manualGC )
807  {
808  m_GraphCutFilter->SetTLinkMode( singleGaussianTLink );
809  m_GraphCutFilter->SetInputSeedSourcesMask( this->GetSourcesMask() );
810  m_GraphCutFilter->SetInputSeedSinksMask( this->GetSinksMask() );
811  }
812  else
813  {
814  this->CreateFuzzyRuleImage();
815 
816  // Compute TLinks if necessary
817  m_TLinksFilter->SetAlpha( m_Alpha );
818  m_TLinksFilter->SetMultiVarSources( m_MultiVarSources );
819  m_TLinksFilter->SetMultiVarSinks( m_MultiVarSinks );
820  m_TLinksFilter->SetTLinkMode( singleGaussianTLink );
821  m_TLinksFilter->SetNumberOfWorkUnits( this->GetNumberOfWorkUnits() );
822  m_TLinksFilter->SetTol( m_Tol );
823 
824  index = 0;
825  if(this->GetInputImageT1().IsNotNull()){m_TLinksFilter->SetInputImage(index, this->GetInputImageT1() ); ++index;}
826  if(this->GetInputImageT2().IsNotNull()){m_TLinksFilter->SetInputImage(index, this->GetInputImageT2() ); ++index;}
827  if(this->GetInputImageDP().IsNotNull()){m_TLinksFilter->SetInputImage(index, this->GetInputImageDP() ); ++index;}
828  if(this->GetInputImageFLAIR().IsNotNull()){m_TLinksFilter->SetInputImage(index, this->GetInputImageFLAIR() ); ++index;}
829  if(this->GetInputImageT1Gd().IsNotNull()){m_TLinksFilter->SetInputImage(index, this->GetInputImageT1Gd() ); ++index;}
830 
831  m_TLinksFilter->SetInputSeedSourcesMask( this->GetSourcesMask() );
832  m_TLinksFilter->SetInputSeedSinksMask( this->GetSinksMask() );
833 
834  if(this->GetSourcesMask().IsNotNull() && m_LesionSegmentationType==gcemAndManualGC)
835  {
836  m_FilterMaxSources->SetInput1( m_FuzzyObject );
837  m_FilterMaxSources->SetInput2( m_TLinksFilter->GetOutputSources() );
838  m_FilterMaxSources->SetNumberOfWorkUnits( this->GetNumberOfWorkUnits() );
839  m_FilterMaxSources->SetCoordinateTolerance( m_Tol );
840  m_FilterMaxSources->SetDirectionTolerance( m_Tol );
841  }
842  if(this->GetSinksMask().IsNotNull() && m_LesionSegmentationType==gcemAndManualGC)
843  {
844  m_FilterMaxSinks->SetInput1( m_MahalanobisFilter->GetOutputMahaMaximum() );
845  m_FilterMaxSinks->SetInput2( m_TLinksFilter->GetOutputSinks() );
846  m_FilterMaxSinks->SetNumberOfWorkUnits( this->GetNumberOfWorkUnits() );
847  m_FilterMaxSinks->SetCoordinateTolerance( m_Tol );
848  m_FilterMaxSinks->SetDirectionTolerance( m_Tol );
849  }
850 
851  m_GraphCutFilter->SetTLinkMode( stremTLink );
852  if(this->GetSourcesMask().IsNotNull() && m_LesionSegmentationType==gcemAndManualGC)
853  {
854  m_GraphCutFilter->SetInputSeedSourcesProba( m_FilterMaxSources->GetOutput() );
855  }
856  else
857  {
858  m_GraphCutFilter->SetInputSeedSourcesProba( m_FuzzyObject );
859  }
860 
861  if(this->GetSinksMask().IsNotNull() && m_LesionSegmentationType==gcemAndManualGC)
862  {
863  m_GraphCutFilter->SetInputSeedSinksProba( m_FilterMaxSinks->GetOutput() );
864  }
865  else if (this->GetInputLesionPrior().IsNotNull())
866  {
867  typedef itk::ImageRegionIterator <ImageTypeD> IteratorType;
868  IteratorType sinkProbasIterator(m_MahalanobisFilter->GetOutputMahaMaximum(),m_MahalanobisFilter->GetOutputMahaMaximum()->GetLargestPossibleRegion());
869  typedef itk::ImageRegionConstIterator <ImageTypeD> ConstIteratorType;
870  ConstIteratorType lesionPriorIterator(this->GetInputLesionPrior(),m_MahalanobisFilter->GetOutputMahaMaximum()->GetLargestPossibleRegion());
871 
872  while (!lesionPriorIterator.IsAtEnd())
873  {
874  double mahaValue = sinkProbasIterator.Get();
875  double lesionPriorValue = lesionPriorIterator.Get();
876 
877  sinkProbasIterator.Set(mahaValue * (1.0 - m_LesionPriorProportion) + (1.0 - lesionPriorValue) * m_LesionPriorProportion);
878 
879  ++sinkProbasIterator;
880  ++lesionPriorIterator;
881  }
882 
883  m_GraphCutFilter->SetInputSeedSinksProba(m_MahalanobisFilter->GetOutputMahaMaximum());
884  }
885  else
886  {
887  m_GraphCutFilter->SetInputSeedSinksProba(m_MahalanobisFilter->GetOutputMahaMaximum());
888  }
889  }
890  m_GraphCutFilter->GraftOutput( this->GetOutputGraphCut() );
891  m_GraphCutFilter->Update();
892  this->GraftNthOutput( 17 , m_GraphCutFilter->GetOutput() );
893 }
894 
895 template <typename TInputImage>
896 void
897 GcStremMsLesionsSegmentationFilter <TInputImage>::ApplyHeuristicRules()
898 {
899  std::cout << "Computing heuristic rules..." << std::endl;
900 
901  ImageTypeUC::Pointer OutliersIntense = ImageTypeUC::New();
902  unsigned int IndexT2ImageInModel, IndexDPImageInModel, IndexFLAIRImageInModel;
903  const unsigned char insideValue = 1;
904  const unsigned char outsideValue = 0;
905 
906  if( m_LesionSegmentationType != manualGC )
907  {
908  std::cout << "Computing hyper-intensity images..." << std::endl;
909 
921  GaussianFunctionType::CovarianceMatrixType covarWM = (m_GaussianModel[m_IndexWMinModel])->GetCovariance();
922  GaussianFunctionType::MeanVectorType meanWM = (m_GaussianModel[m_IndexWMinModel])->GetMean();
923 
924  if(m_UseT2 && m_UseDP)
925  {
926  IndexT2ImageInModel = 1;
927  double wmMeanT2 = meanWM [IndexT2ImageInModel];
928  double wmCovarT2 = covarWM( IndexT2ImageInModel, IndexT2ImageInModel );
929  m_HyperIntensityThreshold1 = wmMeanT2 + m_IntensityT2Factor * sqrt( wmCovarT2 );
930 
931  IndexDPImageInModel = 2;
932  double wmMeanDP = meanWM[IndexDPImageInModel];
933  double wmCovarDP = covarWM( IndexDPImageInModel, IndexDPImageInModel );
934  m_HyperIntensityThreshold2 = wmMeanDP + m_IntensityDPFactor * sqrt( wmCovarDP );
935  }
936  if(m_UseT2 && m_UseFLAIR)
937  {
938  IndexT2ImageInModel = 1;
939  double wmMeanT2 = meanWM[IndexT2ImageInModel];
940  double wmCovarT2 = covarWM( IndexT2ImageInModel, IndexT2ImageInModel );
941  m_HyperIntensityThreshold1 = wmMeanT2 + m_IntensityT2Factor * sqrt( wmCovarT2 );
942 
943  IndexFLAIRImageInModel = 2;
944  double wmMeanFLAIR = meanWM[IndexFLAIRImageInModel];
945  double wmCovarFLAIR = covarWM( IndexFLAIRImageInModel, IndexFLAIRImageInModel );
946  m_HyperIntensityThreshold2 = wmMeanFLAIR + m_IntensityFLAIRFactor * sqrt( wmCovarFLAIR );
947  }
948  if(m_UseDP && m_UseFLAIR)
949  {
950  IndexDPImageInModel = 1;
951  double wmMeanDP = meanWM[IndexDPImageInModel];
952  double wmCovarDP = covarWM( IndexDPImageInModel, IndexDPImageInModel );
953  m_HyperIntensityThreshold1 = wmMeanDP + m_IntensityDPFactor * sqrt( wmCovarDP );
954 
955  IndexFLAIRImageInModel = 2;
956  double wmMeanFLAIR = meanWM[IndexFLAIRImageInModel];
957  double wmCovarFLAIR = covarWM( IndexFLAIRImageInModel, IndexFLAIRImageInModel );
958  m_HyperIntensityThreshold2 = wmMeanFLAIR + m_IntensityFLAIRFactor * sqrt( wmCovarFLAIR );
959 
960  if( m_InitMethodType == HierarchicalDP )
961  {
962  m_HyperIntensityThreshold1 = wmMeanFLAIR + m_IntensityFLAIRFactor * sqrt( wmCovarFLAIR );
963  m_HyperIntensityThreshold2 = wmMeanDP + m_IntensityFLAIRFactor * sqrt( wmCovarDP );
964  }
965  }
966 
967  BinaryThresholdImageFilterType_UC_UC::Pointer IntensityFilter1 = BinaryThresholdImageFilterType_UC_UC::New();
968  IntensityFilter1->SetInput( m_InputImage_1_UC );
969  IntensityFilter1->SetLowerThreshold( m_HyperIntensityThreshold1 );
970  IntensityFilter1->SetInsideValue( insideValue );
971  IntensityFilter1->SetOutsideValue( outsideValue );
972  IntensityFilter1->SetNumberOfWorkUnits( this->GetNumberOfWorkUnits() );
973  IntensityFilter1->SetCoordinateTolerance( m_Tol );
974  IntensityFilter1->SetDirectionTolerance( m_Tol );
975 
976  MaskFilterType_UC_UC::Pointer maskFilterIntensity1 = MaskFilterType_UC_UC::New();
977  maskFilterIntensity1->SetInput( IntensityFilter1->GetOutput() ) ;
978  maskFilterIntensity1->SetMaskImage( m_MaskUC );
979  maskFilterIntensity1->SetNumberOfWorkUnits( this->GetNumberOfWorkUnits() );
980  maskFilterIntensity1->SetCoordinateTolerance( m_Tol );
981  maskFilterIntensity1->SetDirectionTolerance( m_Tol );
982  maskFilterIntensity1->GraftOutput( this->GetOutputIntensityImage1() );
983  maskFilterIntensity1->Update();
984  this->GraftNthOutput( 9 , maskFilterIntensity1->GetOutput() );
985 
986  BinaryThresholdImageFilterType_UC_UC::Pointer IntensityFilter2 = BinaryThresholdImageFilterType_UC_UC::New();
987  IntensityFilter2->SetInput( m_InputImage_2_UC );
988  IntensityFilter2->SetLowerThreshold( m_HyperIntensityThreshold2 );
989  IntensityFilter2->SetInsideValue( insideValue );
990  IntensityFilter2->SetOutsideValue( outsideValue );
991  IntensityFilter2->SetNumberOfWorkUnits( this->GetNumberOfWorkUnits() );
992  IntensityFilter2->SetCoordinateTolerance( m_Tol );
993  IntensityFilter2->SetDirectionTolerance( m_Tol );
994 
995  MaskFilterType_UC_UC::Pointer maskFilterIntensity2 = MaskFilterType_UC_UC::New();
996  maskFilterIntensity2->SetInput( IntensityFilter2->GetOutput() ) ;
997  maskFilterIntensity2->SetMaskImage( m_MaskUC );
998  maskFilterIntensity2->SetNumberOfWorkUnits( this->GetNumberOfWorkUnits() );
999  maskFilterIntensity2->SetCoordinateTolerance( m_Tol );
1000  maskFilterIntensity2->SetDirectionTolerance( m_Tol );
1001  maskFilterIntensity2->GraftOutput( this->GetOutputIntensityImage2() );
1002  maskFilterIntensity2->Update();
1003  this->GraftNthOutput( 10 , maskFilterIntensity2->GetOutput() );
1004 
1005  MinimumFilterTypeUC::Pointer filtermin1 = MinimumFilterTypeUC::New();
1006  filtermin1->SetInput1( maskFilterIntensity1->GetOutput() );
1007  filtermin1->SetInput2( maskFilterIntensity2->GetOutput() );
1008  filtermin1->SetNumberOfWorkUnits( this->GetNumberOfWorkUnits() );
1009  filtermin1->SetCoordinateTolerance( m_Tol );
1010  filtermin1->SetDirectionTolerance( m_Tol );
1011 
1012  MinimumFilterTypeUC::Pointer filtermin2 = MinimumFilterTypeUC::New();
1013  filtermin2->SetInput1( filtermin1->GetOutput() );
1014  filtermin2->SetInput2( m_LesionsDetectionImage );
1015  filtermin2->SetNumberOfWorkUnits( this->GetNumberOfWorkUnits() );
1016  filtermin2->SetCoordinateTolerance( m_Tol );
1017  filtermin2->SetDirectionTolerance( m_Tol );
1018  filtermin2->Update();
1019 
1020  OutliersIntense = filtermin2->GetOutput();
1021  }
1022  else
1023  {
1024  OutliersIntense = m_LesionsDetectionImage;
1025  }
1026 
1027 
1028  typename RemoveTouchingBorderFilterType::Pointer BorderMaskFilter = RemoveTouchingBorderFilterType::New();
1029  if( m_RemoveBorder )
1030  {
1031  std::cout << "Remove lesions touching border..." << std::endl;
1032  BorderMaskFilter->SetMask( m_MaskUC );
1033  BorderMaskFilter->SetInputImageSeg( OutliersIntense );
1034  BorderMaskFilter->SetVerbose( m_Verbose );
1035  BorderMaskFilter->SetTol( m_Tol );
1036  BorderMaskFilter->SetNumberOfWorkUnits( this->GetNumberOfWorkUnits() );
1037  BorderMaskFilter->Update();
1038  }
1039 
1040  ImageTypeUC::Pointer bigLesionsImage = ImageTypeUC::New();
1041 
1042  if( m_LesionSegmentationType != manualGC )
1043  {
1044  std::cout << "Computing white matter ratio..." << std::endl;
1045  BinaryThresholdImageFilterType_F_UC::Pointer thresholdFilterMapWM = BinaryThresholdImageFilterType_F_UC::New();
1046  thresholdFilterMapWM->SetInput( m_MahalanobisFilter->GetOutputMahaWM() );
1047  thresholdFilterMapWM->SetLowerThreshold( m_ThresoldWMmap );
1048  thresholdFilterMapWM->SetInsideValue( insideValue );
1049  thresholdFilterMapWM->SetOutsideValue( outsideValue );
1050  thresholdFilterMapWM->SetNumberOfWorkUnits( this->GetNumberOfWorkUnits() );
1051  thresholdFilterMapWM->SetCoordinateTolerance( m_Tol );
1052  thresholdFilterMapWM->SetDirectionTolerance( m_Tol );
1053 
1054  CheckStructureNeighborFilterFilterType::Pointer CheckStructureNeighborFilterFilter = CheckStructureNeighborFilterFilterType::New();
1055  CheckStructureNeighborFilterFilter->SetInputMap( thresholdFilterMapWM->GetOutput() );
1056  CheckStructureNeighborFilterFilter->SetNumberOfWorkUnits( this->GetNumberOfWorkUnits() );
1057  CheckStructureNeighborFilterFilter->SetTol( m_Tol );
1058  if( m_RemoveBorder )
1059  {
1060  CheckStructureNeighborFilterFilter->SetInputClassification( BorderMaskFilter->GetOutputNonTouchingBorder() );
1061  }
1062  else
1063  {
1064  CheckStructureNeighborFilterFilter->SetInputClassification( OutliersIntense );
1065  }
1066  CheckStructureNeighborFilterFilter->SetRatio( m_RatioContourWM );
1067  CheckStructureNeighborFilterFilter->Update();
1068  bigLesionsImage = CheckStructureNeighborFilterFilter->GetOutput();
1069  }
1070  else
1071  {
1072  if( m_RemoveBorder )
1073  {
1074  bigLesionsImage = BorderMaskFilter->GetOutputNonTouchingBorder();
1075  }
1076  else
1077  {
1078  bigLesionsImage = OutliersIntense ;
1079  }
1080  }
1081 
1082 
1087  ImageTypeUC::SpacingType spacing = m_MaskUC->GetSpacing();
1088  ImageTypeUC::SpacingValueType spacingTot = spacing[0];
1089  for (unsigned int i = 1; i < 3;++i)
1090  {
1091  spacingTot *= spacing[i];
1092  }
1093 
1094  // Compute minimum lesion size in number of pixels
1095  double epsilon = 10e-6;
1096  double minSizeInVoxelD = m_MinLesionsSize / spacingTot;
1097  minSizeInVoxelD -= epsilon;
1098  double minSizeInVoxelD_ceil = std::ceil( minSizeInVoxelD );
1099  unsigned int minSizeInVoxel = static_cast<unsigned int>( minSizeInVoxelD_ceil );
1100  if( minSizeInVoxel > 1 )
1101  {
1102  std::cout << "Removing lesions smaller than " << minSizeInVoxel << "mm3 ..." << std::endl;
1103  }
1104 
1105  bool connectivity = false;
1106  ConnectedComponentType::Pointer ccFilterSize = ConnectedComponentType::New();
1107  ccFilterSize->SetInput( bigLesionsImage );
1108  ccFilterSize->SetFullyConnected( connectivity );
1109  ccFilterSize->SetNumberOfWorkUnits( this->GetNumberOfWorkUnits() );
1110  ccFilterSize->SetCoordinateTolerance( m_Tol );
1111  ccFilterSize->SetDirectionTolerance( m_Tol );
1112 
1113  RelabelComponentType::Pointer relabelFilterSize = RelabelComponentType::New();
1114  relabelFilterSize->SetInput( ccFilterSize->GetOutput() );
1115  relabelFilterSize->SetMinimumObjectSize( minSizeInVoxel );
1116  relabelFilterSize->SetNumberOfWorkUnits( this->GetNumberOfWorkUnits() );
1117  relabelFilterSize->SetCoordinateTolerance( m_Tol );
1118  relabelFilterSize->SetDirectionTolerance( m_Tol );
1119  relabelFilterSize->Update();
1120 
1121  m_LabeledLesions = ImageTypeInt::New();
1122  m_LabeledLesions->Graft(relabelFilterSize->GetOutput());
1123 }
1124 
1125 template <typename TInputImage>
1126 void
1127 GcStremMsLesionsSegmentationFilter <TInputImage>::ComputeNABT()
1128 {
1129  std::cout << "Computing normal appearing brain tissues..." << std::endl;
1130  OutputImagePointer ImageCSF = this->GetOutputCSF();
1131  ImageCSF->SetRegions( m_MaskUC->GetLargestPossibleRegion() );
1132  ImageCSF->CopyInformation( m_MaskUC );
1133  ImageCSF->Allocate();
1134  ImageCSF->FillBuffer(0);
1135 
1136  OutputImagePointer ImageGM = this->GetOutputGM();
1137  ImageGM->SetRegions( m_MaskUC->GetLargestPossibleRegion() );
1138  ImageGM->CopyInformation( m_MaskUC );
1139  ImageGM->Allocate();
1140  ImageGM->FillBuffer(0);
1141 
1142  OutputImagePointer ImageWM = this->GetOutputWM();
1143  ImageWM->SetRegions( m_MaskUC->GetLargestPossibleRegion() );
1144  ImageWM->CopyInformation( m_MaskUC );
1145  ImageWM->Allocate();
1146  ImageWM->FillBuffer(0);
1147 
1148  OutputImagePointer ImageWholeSeg = this->GetOutputWholeSeg();
1149  ImageWholeSeg->SetRegions( m_MaskUC->GetLargestPossibleRegion() );
1150  ImageWholeSeg->CopyInformation( m_MaskUC );
1151  ImageWholeSeg->Allocate();
1152  ImageWholeSeg->FillBuffer(0);
1153 
1154  OutputImagePointer ImageLesions = this->GetOutputLesions();
1155  ImageLesions->SetRegions( m_MaskUC->GetLargestPossibleRegion() );
1156  ImageLesions->CopyInformation( m_MaskUC );
1157  ImageLesions->Allocate();
1158  ImageLesions->FillBuffer(0);
1159 
1160  OutputIteratorType itOutputImageCSF( ImageCSF, ImageCSF->GetLargestPossibleRegion() );
1161  OutputIteratorType itOutputImageGM( ImageGM, ImageGM->GetLargestPossibleRegion() );
1162  OutputIteratorType itOutputImageWM( ImageWM, ImageWM->GetLargestPossibleRegion() );
1163  OutputIteratorType itOutputImageWholeSeg( ImageWholeSeg, ImageWholeSeg->GetLargestPossibleRegion() );
1164  OutputIteratorType itOutputImageLesions( ImageLesions, ImageLesions->GetLargestPossibleRegion() );
1165  ImageIteratorTypeInt itImageLesions( m_LabeledLesions, m_LabeledLesions->GetLargestPossibleRegion() );
1166 
1167  ImageIteratorTypeUC itMask( m_MaskUC, m_MaskUC->GetLargestPossibleRegion() );
1168  while (!itMask.IsAtEnd())
1169  {
1170  if((itMask.Get()!=0) && (itImageLesions.Get() != 0))
1171  {
1172  itOutputImageLesions.Set(1);
1173  itOutputImageWholeSeg.Set(static_cast<OutputPixelType>(m_LabelLesions));
1174  }
1175  ++itOutputImageLesions;
1176  ++itOutputImageWholeSeg;
1177  ++itImageLesions;
1178  ++itMask;
1179  }
1180 
1181  itOutputImageWholeSeg.GoToBegin();
1182  if(m_LesionSegmentationType!=manualGC)
1183  {
1184  double minimumThWM = 0.3; // set a minimum proba value to avoid classifying grey matter as white matter
1185 
1186  ImageIteratorTypeD itmahaWM(m_MahalanobisFilter->GetOutputMahaWM(),m_MahalanobisFilter->GetOutputMahaWM()->GetLargestPossibleRegion());
1187  ImageIteratorTypeD itmahaGM(m_MahalanobisFilter->GetOutputMahaGM(),m_MahalanobisFilter->GetOutputMahaGM()->GetLargestPossibleRegion());
1188  ImageIteratorTypeD itmahaCSF(m_MahalanobisFilter->GetOutputMahaCSF(),m_MahalanobisFilter->GetOutputMahaCSF()->GetLargestPossibleRegion());
1189 
1190  itMask.GoToBegin();
1191  itImageLesions.GoToBegin();
1192 
1193  while (!itMask.IsAtEnd())
1194  {
1195  if((itMask.Get()!=0) && (itImageLesions.Get() == 0))
1196  {
1197  if( (itmahaCSF.Get()>itmahaGM.Get()) && (itmahaCSF.Get()>itmahaWM.Get()) )
1198  {
1199  itOutputImageCSF.Set(1);
1200  itOutputImageWholeSeg.Set(static_cast<OutputPixelType>(m_LabelCSF));
1201  }
1202  else if( (itmahaGM.Get()>itmahaCSF.Get()) && (itmahaGM.Get()>itmahaWM.Get()) )
1203  {
1204  itOutputImageGM.Set(1);
1205  itOutputImageWholeSeg.Set(static_cast<OutputPixelType>(m_LabelGM));
1206  }
1207  else if( (itmahaWM.Get()>itmahaGM.Get()) && (itmahaWM.Get()>itmahaCSF.Get()) && (itmahaWM.Get()>minimumThWM) )
1208  {
1209  itOutputImageWM.Set(1);
1210  itOutputImageWholeSeg.Set(static_cast<OutputPixelType>(m_LabelWM));
1211  }
1212  }
1213  ++itOutputImageCSF;
1214  ++itOutputImageGM;
1215  ++itOutputImageWM;
1216  ++itOutputImageWholeSeg;
1217  ++itOutputImageLesions;
1218  ++itImageLesions;
1219  ++itmahaWM;
1220  ++itmahaGM;
1221  ++itmahaCSF;
1222  ++itMask;
1223  }
1224  }
1225 }
1226 
1227 
1228 template <typename TInputImage>
1229 void
1230 GcStremMsLesionsSegmentationFilter <TInputImage>::GenerateData()
1231 {
1232  this->CheckInputImages();
1233 
1234  this->RescaleImages();
1235 
1236  // Compute NABT model estimation
1237  if (m_LesionSegmentationType != manualGC)
1238  {
1239  this->ComputeAutomaticInitialization();
1240  }
1241 
1242  // Lesions detection by graph cut or thresholding
1243  if( m_LesionSegmentationType == strem )
1244  {
1245  this->StremThreshold();
1246  m_LesionsDetectionImage = this->GetOutputStrem();
1247  }
1248  else
1249  {
1250  this->GraphCut();
1251  m_LesionsDetectionImage = this->GetOutputGraphCut();
1252  }
1253 
1254  // Remove false positives
1255  this->ApplyHeuristicRules();
1256 
1257  // Compute NABT maps
1258  this->ComputeNABT();
1259 }
1260 
1261 }
anima::GcStremMsLesionsSegmentationFilter::OutputIteratorType
itk::ImageRegionIterator< OutputImageType > OutputIteratorType
Definition: animaGcStremMsLesionsSegmentationFilter.h:95
anima::GcStremMsLesionsSegmentationFilter::ImageIteratorTypeInt
itk::ImageRegionIterator< ImageTypeInt > ImageIteratorTypeInt
Definition: animaGcStremMsLesionsSegmentationFilter.h:109
anima::GcStremMsLesionsSegmentationFilter::ImageIteratorTypeConstInput
itk::ImageRegionConstIterator< InputImageType > ImageIteratorTypeConstInput
Definition: animaGcStremMsLesionsSegmentationFilter.h:90
anima::ComputeSolution::Pointer
itk::SmartPointer< Self > Pointer
Definition: animaComputeSolution.h:33
anima::GcStremMsLesionsSegmentationFilter::SetAlpha
virtual void SetAlpha(double _arg)
Generated by   doxygen 1.8.13