6 template <
typename TInputImage>
9 this->SetNthInput(m_NbInputs, const_cast<TInputImage*>(mask));
10 m_IndexMask = m_NbInputs;
14 template <
typename TInputImage>
17 this->SetNthInput(m_NbInputs, const_cast<TInputImage*>(image));
18 m_IndexImageT1 = m_NbInputs;
23 template <
typename TInputImage>
26 this->SetNthInput(m_NbInputs, const_cast<TInputImage*>(image));
27 m_IndexImageT2 = m_NbInputs;
32 template <
typename TInputImage>
35 this->SetNthInput(m_NbInputs, const_cast<TInputImage*>(image));
36 m_IndexImageDP = m_NbInputs;
41 template <
typename TInputImage>
44 this->SetNthInput(m_NbInputs, const_cast<TInputImage*>(image));
45 m_IndexImageFLAIR = m_NbInputs;
50 template <
typename TInputImage>
53 this->SetNthInput(m_NbInputs, const_cast<TInputImage*>(image));
54 m_IndexImageT1Gd = m_NbInputs;
59 template <
typename TInputImage>
62 this->SetNthInput(m_NbInputs, const_cast<ImageTypeD*>(atlas));
63 m_IndexAtlasCSF = m_NbInputs;
67 template <
typename TInputImage>
70 this->SetNthInput(m_NbInputs, const_cast<ImageTypeD*>(atlas));
71 m_IndexAtlasGM = m_NbInputs;
75 template <
typename TInputImage>
78 this->SetNthInput(m_NbInputs, const_cast<ImageTypeD*>(atlas));
79 m_IndexAtlasWM = m_NbInputs;
83 template <
typename TInputImage>
86 this->SetNthInput(m_NbInputs, const_cast<ImageTypeD*>(image));
87 m_IndexLesionPrior = m_NbInputs;
91 template <
typename TInputImage>
94 this->SetNthInput(m_NbInputs, const_cast<ImageTypeUC*>(image));
95 m_IndexSourcesMask = m_NbInputs;
99 template <
typename TInputImage>
102 this->SetNthInput(m_NbInputs, const_cast<ImageTypeUC*>(image));
103 m_IndexSinksMask = m_NbInputs;
109 template <
typename TInputImage>
112 return static_cast< const TInputImage *
> 113 ( this->itk::ProcessObject::GetInput(m_IndexMask) );
116 template <
typename TInputImage>
119 return static_cast< const TInputImage *
> 120 ( this->itk::ProcessObject::GetInput(m_IndexImageT1) );
123 template <
typename TInputImage>
126 return static_cast< const TInputImage *
> 127 ( this->itk::ProcessObject::GetInput(m_IndexImageT2) );
130 template <
typename TInputImage>
133 return static_cast< const TInputImage *
> 134 ( this->itk::ProcessObject::GetInput(m_IndexImageDP) );
137 template <
typename TInputImage>
140 return static_cast< const TInputImage *
> 141 ( this->itk::ProcessObject::GetInput(m_IndexImageFLAIR) );
143 template <
typename TInputImage>
146 return static_cast< const TInputImage *
> 147 ( this->itk::ProcessObject::GetInput(m_IndexImageT1Gd) );
151 template <
typename TInputImage>
155 ( this->itk::ProcessObject::GetInput(m_IndexAtlasCSF) );
158 template <
typename TInputImage>
162 ( this->itk::ProcessObject::GetInput(m_IndexAtlasGM) );
165 template <
typename TInputImage>
169 ( this->itk::ProcessObject::GetInput(m_IndexAtlasWM) );
172 template <
typename TInputImage>
175 return static_cast< const ImageTypeD *
> (this->itk::ProcessObject::GetInput(m_IndexLesionPrior));
178 template <
typename TInputImage>
182 ( this->itk::ProcessObject::GetInput(m_IndexSourcesMask) );
184 template <
typename TInputImage>
188 ( this->itk::ProcessObject::GetInput(m_IndexSinksMask) );
192 template <
typename TInputImage>
195 itk::DataObject::Pointer output;
197 if( (idx < 11) || (idx == 17) )
199 output = ( TOutputImage::New() ).GetPointer();
201 else if(( 11 <= idx ) && ( idx < 17 ))
203 output = ( ImageTypeD::New() ).GetPointer();
207 std::cerr <<
"No output " << idx << std::endl;
210 return output.GetPointer();
213 template <
typename TInputImage>
216 return dynamic_cast< TOutputImage*
>( this->itk::ProcessObject::GetOutput(0) );
218 template <
typename TInputImage>
221 return dynamic_cast< TOutputImage*
>( this->itk::ProcessObject::GetOutput(1) );
223 template <
typename TInputImage>
226 return dynamic_cast< TOutputImage*
>( this->itk::ProcessObject::GetOutput(2) );
228 template <
typename TInputImage>
231 return dynamic_cast< TOutputImage*
>( this->itk::ProcessObject::GetOutput(3) );
233 template <
typename TInputImage>
236 return dynamic_cast< TOutputImage*
>( this->itk::ProcessObject::GetOutput(4) );
238 template <
typename TInputImage>
241 return dynamic_cast< TOutputImage*
>( this->itk::ProcessObject::GetOutput(5) );
243 template <
typename TInputImage>
246 return dynamic_cast< TOutputImage*
>( this->itk::ProcessObject::GetOutput(6) );
248 template <
typename TInputImage>
251 return dynamic_cast< TOutputImage*
>( this->itk::ProcessObject::GetOutput(7) );
253 template <
typename TInputImage>
256 return dynamic_cast< TOutputImage*
>( this->itk::ProcessObject::GetOutput(8) );
258 template <
typename TInputImage>
261 return dynamic_cast< TOutputImage*
>( this->itk::ProcessObject::GetOutput(9) );
263 template <
typename TInputImage>
266 return dynamic_cast< TOutputImage*
>( this->itk::ProcessObject::GetOutput(10) );
268 template <
typename TInputImage>
271 return dynamic_cast< ImageTypeD*
>( this->itk::ProcessObject::GetOutput(11) );
273 template <
typename TInputImage>
276 return dynamic_cast< ImageTypeD*
>( this->itk::ProcessObject::GetOutput(12) );
278 template <
typename TInputImage>
281 return dynamic_cast< ImageTypeD*
>( this->itk::ProcessObject::GetOutput(13) );
283 template <
typename TInputImage>
286 return dynamic_cast< ImageTypeD*
>( this->itk::ProcessObject::GetOutput(14) );
288 template <
typename TInputImage>
291 return dynamic_cast< ImageTypeD*
>( this->itk::ProcessObject::GetOutput(15) );
293 template <
typename TInputImage>
296 return dynamic_cast< ImageTypeD*
>( this->itk::ProcessObject::GetOutput(16) );
298 template <
typename TInputImage>
301 return dynamic_cast< TOutputImage*
>( this->itk::ProcessObject::GetOutput(17) );
305 template <
typename TInputImage>
309 if( m_OutputLesionFilename !=
"" )
311 std::cout <<
"Writing Lesion output image to: " << m_OutputLesionFilename << std::endl;
312 anima::writeImage<TOutputImage>(m_OutputLesionFilename, this->GetOutputLesions());
315 if( m_LesionSegmentationType!=
manualGC )
318 if( m_OutputIntensityImage1Filename !=
"" )
320 std::cout <<
"Writing intensity output image 1 to: " << m_OutputIntensityImage1Filename << std::endl;
321 anima::writeImage<TOutputImage>(m_OutputIntensityImage1Filename, this->GetOutputIntensityImage1());
323 if( m_OutputIntensityImage2Filename !=
"" )
325 std::cout <<
"Writing intensity output image 2 to: " << m_OutputIntensityImage2Filename << std::endl;
326 anima::writeImage<TOutputImage>(m_OutputIntensityImage2Filename, this->GetOutputIntensityImage2());
328 if( m_OutputCSFFilename !=
"" )
330 std::cout <<
"Writing CSF output image to: " << m_OutputCSFFilename << std::endl;
331 anima::writeImage<TOutputImage>(m_OutputCSFFilename, this->GetOutputCSF());
333 if( m_OutputGMFilename !=
"" )
335 std::cout <<
"Writing GM output image to: " << m_OutputGMFilename << std::endl;
336 anima::writeImage<TOutputImage>(m_OutputGMFilename, this->GetOutputGM());
338 if( m_OutputWMFilename !=
"" )
340 std::cout <<
"Writing WM output image to: " << m_OutputWMFilename << std::endl;
341 anima::writeImage<TOutputImage>(m_OutputWMFilename, this->GetOutputWM());
343 if( m_OutputWholeFilename !=
"" )
345 std::cout <<
"Writing segmentation output image to: " << m_OutputWholeFilename << std::endl;
346 anima::writeImage<TOutputImage>(m_OutputWholeFilename, this->GetOutputWholeSeg());
352 if( m_OutputFuzzyObjectFilename !=
"" )
354 std::cout <<
"Writing fuzzy object output image to: " << m_OutputFuzzyObjectFilename << std::endl;
355 anima::writeImage<ImageTypeD>(m_OutputFuzzyObjectFilename, this->GetOutputFuzzyObjectImage());
359 if( m_LesionSegmentationType==
strem )
361 if( m_OutputStremFilename !=
"" )
363 std::cout <<
"Writing strem output image to: " << m_OutputStremFilename << std::endl;
364 anima::writeImage<TOutputImage>(m_OutputStremFilename, this->GetOutputStrem());
366 if( m_OutputStremCSFFilename !=
"" )
368 std::cout <<
"Writing CSF strem output image to: " << m_OutputStremCSFFilename << std::endl;
369 anima::writeImage<TOutputImage>(m_OutputStremCSFFilename, this->GetOutputStremCSF());
371 if( m_OutputStremGMFilename !=
"" )
373 std::cout <<
"Writing GM strem output image to: " << m_OutputStremGMFilename << std::endl;
374 anima::writeImage<TOutputImage>(m_OutputStremGMFilename, this->GetOutputStremGM());
376 if( m_OutputStremWMFilename !=
"" )
378 std::cout <<
"Writing WM strem output image to: " << m_OutputStremWMFilename << std::endl;
379 anima::writeImage<TOutputImage>(m_OutputStremWMFilename, this->GetOutputStremWM());
384 m_GraphCutFilter->WriteOutputs();
389 template <
typename TInputImage>
393 std::cout <<
"Check input entries..." << std::endl;
394 if(m_LesionSegmentationType!=
manualGC)
397 itkExceptionMacro(
"-- Error in anima::Gc_Strem_MS_Lesions_Segmentation: automatic segmentation requires 3 modalities, exiting...");
399 if(this->GetInputImageT1().IsNull())
400 itkExceptionMacro(
"-- Error in anima::Gc_Strem_MS_Lesions_Segmentation: automatic segmentation requires T1 image, exiting...");
402 if(m_Modalities == 3)
408 if(this->GetInputImageT2().IsNotNull()){m_UseT2=
true;}
409 if(this->GetInputImageDP().IsNotNull()){m_UseDP=
true;}
410 if(this->GetInputImageFLAIR().IsNotNull()){m_UseFLAIR=
true;}
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...");
421 if(this->GetSourcesMask().IsNull() || this->GetSinksMask().IsNull())
422 itkExceptionMacro(
"-- Error in anima::Gc_Strem_MS_Lesions_Segmentation: seed mask missing, exiting...");
425 template <
typename TInputImage>
429 std::cout <<
"Rescale input images..." << std::endl;
431 ImageTypeUC::Pointer InputImage_T2_UC = ImageTypeUC::New();
432 ImageTypeUC::Pointer InputImage_DP_UC = ImageTypeUC::New();
433 ImageTypeUC::Pointer InputImage_FLAIR_UC = ImageTypeUC::New();
435 double desiredMinimum=0,desiredMaximum=255;
437 if(this->GetInputImageT1().IsNotNull())
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();
449 if(this->GetInputImageT2().IsNotNull())
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();
461 if(this->GetInputImageFLAIR().IsNotNull())
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();
473 if(this->GetInputImageDP().IsNotNull())
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();
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);
496 while (!itMaskUC.IsAtEnd())
498 if(itInputImage_Mask.Get()!=0)
506 if(m_LesionSegmentationType!=
manualGC)
508 if(m_UseT2 && m_UseDP)
510 m_InputImage_1_UC = InputImage_T2_UC;
511 m_InputImage_2_UC = InputImage_DP_UC;
513 if(m_UseT2 && m_UseFLAIR)
515 m_InputImage_1_UC = InputImage_T2_UC;
516 m_InputImage_2_UC = InputImage_FLAIR_UC;
518 if(m_UseDP && m_UseFLAIR)
520 m_InputImage_1_UC = InputImage_DP_UC;
521 m_InputImage_2_UC = InputImage_FLAIR_UC;
523 if(m_InitMethodType==1)
525 m_InputImage_1_UC = InputImage_FLAIR_UC;
526 m_InputImage_2_UC = InputImage_DP_UC;
532 template <
typename TInputImage>
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 );
543 ComputeSolutionProcess->SetInputCSFAtlas( this->GetInputCSFAtlas() );
544 ComputeSolutionProcess->SetInputGMAtlas( this->GetInputGMAtlas() );
545 ComputeSolutionProcess->SetInputWMAtlas( this->GetInputWMAtlas() );
547 ComputeSolutionProcess->SetSolutionWriteFilename( m_SolutionWriteFilename );
548 ComputeSolutionProcess->SetSolutionReadFilename( m_SolutionReadFilename );
549 ComputeSolutionProcess->SetGaussianModel( m_GaussianModel );
550 ComputeSolutionProcess->SetAlphas( m_Alphas );
552 ComputeSolutionProcess->SetUseT2( m_UseT2 );
553 ComputeSolutionProcess->SetUseDP( m_UseDP );
554 ComputeSolutionProcess->SetUseFLAIR( m_UseFLAIR );
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 );
564 ComputeSolutionProcess->SetVerbose( m_Verbose );
568 ComputeSolutionProcess->Update();
570 catch (itk::ExceptionObject &e)
572 std::cerr << e << std::endl;
576 m_GaussianModel = ComputeSolutionProcess->GetGaussianModel();
578 std::cout <<
"Computing mahalanobis images..." << std::endl;
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() );
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 );
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() );
608 template <
typename TInputImage>
612 std::cout <<
"Compute strem threshold images..." << std::endl;
614 const unsigned char insideValue = 1;
615 const unsigned char outsideValue = 0;
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 );
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() );
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 );
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() );
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 );
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() );
674 MinimumFilterTypeUC::Pointer filtermin1 = MinimumFilterTypeUC::New();
675 MinimumFilterTypeUC::Pointer filtermin2 = MinimumFilterTypeUC::New();
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);
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() );
693 template <
typename TInputImage>
697 std::cout <<
"Computing fuzzy rules..." << std::endl;
699 MinimumFilterTypeF::Pointer filtermin1 = MinimumFilterTypeF::New();
700 MinimumFilterTypeF::Pointer filtermin2 = MinimumFilterTypeF::New();
702 const double minimumOutputValue = 0;
703 const double maximumOutputValue = 1;
705 const unsigned int indexWM = 2;
706 GaussianFunctionType::MeanVectorType mean = (m_GaussianModel[indexWM])->GetMean();
707 GaussianFunctionType::CovarianceMatrixType covar = (m_GaussianModel[indexWM])->GetCovariance();
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]) );
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 );
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]) );
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 );
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 );
743 if (this->GetInputLesionPrior().IsNotNull())
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());
750 while (!lesionPriorIterator.IsAtEnd())
752 double mahaValue = mahaMinimumIterator.Get();
753 double lesionPriorValue = lesionPriorIterator.Get();
755 mahaMinimumIterator.Set(mahaValue * (1.0 - m_LesionPriorProportion) + lesionPriorValue * m_LesionPriorProportion);
757 ++mahaMinimumIterator;
758 ++lesionPriorIterator;
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 );
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() );
778 m_FuzzyObject = ImageTypeD::New();
779 m_FuzzyObject->Graft(maskFilter->GetOutput());
782 template <
typename TInputImage>
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 );
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;}
804 m_GraphCutFilter->SetMask( m_MaskUC );
806 if( m_LesionSegmentationType==
manualGC )
809 m_GraphCutFilter->SetInputSeedSourcesMask( this->GetSourcesMask() );
810 m_GraphCutFilter->SetInputSeedSinksMask( this->GetSinksMask() );
814 this->CreateFuzzyRuleImage();
817 m_TLinksFilter->SetAlpha( m_Alpha );
818 m_TLinksFilter->SetMultiVarSources( m_MultiVarSources );
819 m_TLinksFilter->SetMultiVarSinks( m_MultiVarSinks );
821 m_TLinksFilter->SetNumberOfWorkUnits( this->GetNumberOfWorkUnits() );
822 m_TLinksFilter->SetTol( m_Tol );
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;}
831 m_TLinksFilter->SetInputSeedSourcesMask( this->GetSourcesMask() );
832 m_TLinksFilter->SetInputSeedSinksMask( this->GetSinksMask() );
834 if(this->GetSourcesMask().IsNotNull() && m_LesionSegmentationType==
gcemAndManualGC)
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 );
842 if(this->GetSinksMask().IsNotNull() && m_LesionSegmentationType==
gcemAndManualGC)
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 );
852 if(this->GetSourcesMask().IsNotNull() && m_LesionSegmentationType==
gcemAndManualGC)
854 m_GraphCutFilter->SetInputSeedSourcesProba( m_FilterMaxSources->GetOutput() );
858 m_GraphCutFilter->SetInputSeedSourcesProba( m_FuzzyObject );
861 if(this->GetSinksMask().IsNotNull() && m_LesionSegmentationType==
gcemAndManualGC)
863 m_GraphCutFilter->SetInputSeedSinksProba( m_FilterMaxSinks->GetOutput() );
865 else if (this->GetInputLesionPrior().IsNotNull())
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());
872 while (!lesionPriorIterator.IsAtEnd())
874 double mahaValue = sinkProbasIterator.Get();
875 double lesionPriorValue = lesionPriorIterator.Get();
877 sinkProbasIterator.Set(mahaValue * (1.0 - m_LesionPriorProportion) + (1.0 - lesionPriorValue) * m_LesionPriorProportion);
879 ++sinkProbasIterator;
880 ++lesionPriorIterator;
883 m_GraphCutFilter->SetInputSeedSinksProba(m_MahalanobisFilter->GetOutputMahaMaximum());
887 m_GraphCutFilter->SetInputSeedSinksProba(m_MahalanobisFilter->GetOutputMahaMaximum());
890 m_GraphCutFilter->GraftOutput( this->GetOutputGraphCut() );
891 m_GraphCutFilter->Update();
892 this->GraftNthOutput( 17 , m_GraphCutFilter->GetOutput() );
895 template <
typename TInputImage>
899 std::cout <<
"Computing heuristic rules..." << std::endl;
901 ImageTypeUC::Pointer OutliersIntense = ImageTypeUC::New();
902 unsigned int IndexT2ImageInModel, IndexDPImageInModel, IndexFLAIRImageInModel;
903 const unsigned char insideValue = 1;
904 const unsigned char outsideValue = 0;
906 if( m_LesionSegmentationType !=
manualGC )
908 std::cout <<
"Computing hyper-intensity images..." << std::endl;
921 GaussianFunctionType::CovarianceMatrixType covarWM = (m_GaussianModel[m_IndexWMinModel])->GetCovariance();
922 GaussianFunctionType::MeanVectorType meanWM = (m_GaussianModel[m_IndexWMinModel])->GetMean();
924 if(m_UseT2 && m_UseDP)
926 IndexT2ImageInModel = 1;
927 double wmMeanT2 = meanWM [IndexT2ImageInModel];
928 double wmCovarT2 = covarWM( IndexT2ImageInModel, IndexT2ImageInModel );
929 m_HyperIntensityThreshold1 = wmMeanT2 + m_IntensityT2Factor * sqrt( wmCovarT2 );
931 IndexDPImageInModel = 2;
932 double wmMeanDP = meanWM[IndexDPImageInModel];
933 double wmCovarDP = covarWM( IndexDPImageInModel, IndexDPImageInModel );
934 m_HyperIntensityThreshold2 = wmMeanDP + m_IntensityDPFactor * sqrt( wmCovarDP );
936 if(m_UseT2 && m_UseFLAIR)
938 IndexT2ImageInModel = 1;
939 double wmMeanT2 = meanWM[IndexT2ImageInModel];
940 double wmCovarT2 = covarWM( IndexT2ImageInModel, IndexT2ImageInModel );
941 m_HyperIntensityThreshold1 = wmMeanT2 + m_IntensityT2Factor * sqrt( wmCovarT2 );
943 IndexFLAIRImageInModel = 2;
944 double wmMeanFLAIR = meanWM[IndexFLAIRImageInModel];
945 double wmCovarFLAIR = covarWM( IndexFLAIRImageInModel, IndexFLAIRImageInModel );
946 m_HyperIntensityThreshold2 = wmMeanFLAIR + m_IntensityFLAIRFactor * sqrt( wmCovarFLAIR );
948 if(m_UseDP && m_UseFLAIR)
950 IndexDPImageInModel = 1;
951 double wmMeanDP = meanWM[IndexDPImageInModel];
952 double wmCovarDP = covarWM( IndexDPImageInModel, IndexDPImageInModel );
953 m_HyperIntensityThreshold1 = wmMeanDP + m_IntensityDPFactor * sqrt( wmCovarDP );
955 IndexFLAIRImageInModel = 2;
956 double wmMeanFLAIR = meanWM[IndexFLAIRImageInModel];
957 double wmCovarFLAIR = covarWM( IndexFLAIRImageInModel, IndexFLAIRImageInModel );
958 m_HyperIntensityThreshold2 = wmMeanFLAIR + m_IntensityFLAIRFactor * sqrt( wmCovarFLAIR );
962 m_HyperIntensityThreshold1 = wmMeanFLAIR + m_IntensityFLAIRFactor * sqrt( wmCovarFLAIR );
963 m_HyperIntensityThreshold2 = wmMeanDP + m_IntensityFLAIRFactor * sqrt( wmCovarDP );
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 );
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() );
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 );
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() );
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 );
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();
1020 OutliersIntense = filtermin2->GetOutput();
1024 OutliersIntense = m_LesionsDetectionImage;
1029 if( m_RemoveBorder )
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();
1040 ImageTypeUC::Pointer bigLesionsImage = ImageTypeUC::New();
1042 if( m_LesionSegmentationType !=
manualGC )
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 );
1055 CheckStructureNeighborFilterFilter->SetInputMap( thresholdFilterMapWM->GetOutput() );
1056 CheckStructureNeighborFilterFilter->SetNumberOfWorkUnits( this->GetNumberOfWorkUnits() );
1057 CheckStructureNeighborFilterFilter->SetTol( m_Tol );
1058 if( m_RemoveBorder )
1060 CheckStructureNeighborFilterFilter->SetInputClassification( BorderMaskFilter->GetOutputNonTouchingBorder() );
1064 CheckStructureNeighborFilterFilter->SetInputClassification( OutliersIntense );
1066 CheckStructureNeighborFilterFilter->SetRatio( m_RatioContourWM );
1067 CheckStructureNeighborFilterFilter->Update();
1068 bigLesionsImage = CheckStructureNeighborFilterFilter->GetOutput();
1072 if( m_RemoveBorder )
1074 bigLesionsImage = BorderMaskFilter->GetOutputNonTouchingBorder();
1078 bigLesionsImage = OutliersIntense ;
1087 ImageTypeUC::SpacingType spacing = m_MaskUC->GetSpacing();
1088 ImageTypeUC::SpacingValueType spacingTot = spacing[0];
1089 for (
unsigned int i = 1; i < 3;++i)
1091 spacingTot *= spacing[i];
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 )
1102 std::cout <<
"Removing lesions smaller than " << minSizeInVoxel <<
"mm3 ..." << std::endl;
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 );
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();
1121 m_LabeledLesions = ImageTypeInt::New();
1122 m_LabeledLesions->Graft(relabelFilterSize->GetOutput());
1125 template <
typename TInputImage>
1129 std::cout <<
"Computing normal appearing brain tissues..." << std::endl;
1131 ImageCSF->SetRegions( m_MaskUC->GetLargestPossibleRegion() );
1132 ImageCSF->CopyInformation( m_MaskUC );
1133 ImageCSF->Allocate();
1134 ImageCSF->FillBuffer(0);
1137 ImageGM->SetRegions( m_MaskUC->GetLargestPossibleRegion() );
1138 ImageGM->CopyInformation( m_MaskUC );
1139 ImageGM->Allocate();
1140 ImageGM->FillBuffer(0);
1143 ImageWM->SetRegions( m_MaskUC->GetLargestPossibleRegion() );
1144 ImageWM->CopyInformation( m_MaskUC );
1145 ImageWM->Allocate();
1146 ImageWM->FillBuffer(0);
1149 ImageWholeSeg->SetRegions( m_MaskUC->GetLargestPossibleRegion() );
1150 ImageWholeSeg->CopyInformation( m_MaskUC );
1151 ImageWholeSeg->Allocate();
1152 ImageWholeSeg->FillBuffer(0);
1155 ImageLesions->SetRegions( m_MaskUC->GetLargestPossibleRegion() );
1156 ImageLesions->CopyInformation( m_MaskUC );
1157 ImageLesions->Allocate();
1158 ImageLesions->FillBuffer(0);
1160 OutputIteratorType itOutputImageCSF( ImageCSF, ImageCSF->GetLargestPossibleRegion() );
1163 OutputIteratorType itOutputImageWholeSeg( ImageWholeSeg, ImageWholeSeg->GetLargestPossibleRegion() );
1164 OutputIteratorType itOutputImageLesions( ImageLesions, ImageLesions->GetLargestPossibleRegion() );
1165 ImageIteratorTypeInt itImageLesions( m_LabeledLesions, m_LabeledLesions->GetLargestPossibleRegion() );
1168 while (!itMask.IsAtEnd())
1170 if((itMask.Get()!=0) && (itImageLesions.Get() != 0))
1172 itOutputImageLesions.Set(1);
1173 itOutputImageWholeSeg.Set(static_cast<OutputPixelType>(m_LabelLesions));
1175 ++itOutputImageLesions;
1176 ++itOutputImageWholeSeg;
1181 itOutputImageWholeSeg.GoToBegin();
1182 if(m_LesionSegmentationType!=
manualGC)
1184 double minimumThWM = 0.3;
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());
1191 itImageLesions.GoToBegin();
1193 while (!itMask.IsAtEnd())
1195 if((itMask.Get()!=0) && (itImageLesions.Get() == 0))
1197 if( (itmahaCSF.Get()>itmahaGM.Get()) && (itmahaCSF.Get()>itmahaWM.Get()) )
1199 itOutputImageCSF.Set(1);
1200 itOutputImageWholeSeg.Set(static_cast<OutputPixelType>(m_LabelCSF));
1202 else if( (itmahaGM.Get()>itmahaCSF.Get()) && (itmahaGM.Get()>itmahaWM.Get()) )
1204 itOutputImageGM.Set(1);
1205 itOutputImageWholeSeg.Set(static_cast<OutputPixelType>(m_LabelGM));
1207 else if( (itmahaWM.Get()>itmahaGM.Get()) && (itmahaWM.Get()>itmahaCSF.Get()) && (itmahaWM.Get()>minimumThWM) )
1209 itOutputImageWM.Set(1);
1210 itOutputImageWholeSeg.Set(static_cast<OutputPixelType>(m_LabelWM));
1216 ++itOutputImageWholeSeg;
1217 ++itOutputImageLesions;
1228 template <
typename TInputImage>
1232 this->CheckInputImages();
1234 this->RescaleImages();
1237 if (m_LesionSegmentationType !=
manualGC)
1239 this->ComputeAutomaticInitialization();
1243 if( m_LesionSegmentationType ==
strem )
1245 this->StremThreshold();
1246 m_LesionsDetectionImage = this->GetOutputStrem();
1251 m_LesionsDetectionImage = this->GetOutputGraphCut();
1255 this->ApplyHeuristicRules();
1258 this->ComputeNABT();
itk::SmartPointer< Self > Pointer
OutputImageType::Pointer OutputImagePointer
itk::Image< PixelTypeD, 3 > ImageTypeD
itk::ImageRegionIterator< OutputImageType > OutputIteratorType
itk::Image< PixelTypeUC, 3 > ImageTypeUC
itk::SmartPointer< Self > Pointer
itk::ImageRegionIterator< ImageTypeUC > ImageIteratorTypeUC
itk::ImageRegionIterator< ImageTypeInt > ImageIteratorTypeInt
itk::ImageRegionConstIterator< InputImageType > ImageIteratorTypeConstInput
itk::SmartPointer< Self > Pointer
Class performing lesion segmentation.
virtual void SetAlpha(double _arg)
itk::ImageRegionIterator< ImageTypeD > ImageIteratorTypeD
itk::Image< unsigned char, 3 > TOutputImage