segmentation/clitkExtractLungFilter.txx

   1 /*=========================================================================
   2   Program:   vv                     http://www.creatis.insa-lyon.fr/rio/vv
   3
   4   Authors belong to:
   5   - University of LYON              http://www.universite-lyon.fr/
   6   - Léon Bérard cancer center       http://oncora1.lyon.fnclcc.fr
   7   - CREATIS CNRS laboratory         http://www.creatis.insa-lyon.fr
   8
   9   This software is distributed WITHOUT ANY WARRANTY; without even
  10   the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  11   PURPOSE.  See the copyright notices for more information.
  12
  13   It is distributed under dual licence
  14
  15   - BSD        See included LICENSE.txt file
  16   - CeCILL-B   http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
  17   ======================================================================-====*/
  18
  19 #ifndef CLITKEXTRACTLUNGSFILTER_TXX
  20 #define CLITKEXTRACTLUNGSFILTER_TXX
  21
  22 // clitk
  23 #include "clitkImageCommon.h"
  24 #include "clitkSetBackgroundImageFilter.h"
  25 #include "clitkSegmentationUtils.h"
  26 #include "clitkAutoCropFilter.h"
  27 #include "clitkCropLikeImageFilter.h"
  28 #include "clitkFillMaskFilter.h"
  29
  30 // itk
  31 #include "itkBinaryThresholdImageFilter.h"
  32 #include "itkConnectedComponentImageFilter.h"
  33 #include "itkRelabelComponentImageFilter.h"
  34 #include "itkOtsuThresholdImageFilter.h"
  35 #include "itkBinaryThinningImageFilter3D.h"
  36 #include "itkImageIteratorWithIndex.h"
  37 #include "itkBinaryMorphologicalOpeningImageFilter.h"
  38 #include "itkBinaryMorphologicalClosingImageFilter.h"
  39
  40 //--------------------------------------------------------------------
  41 template <class ImageType>
  42 clitk::ExtractLungFilter<ImageType>::
  43 ExtractLungFilter():
  44   clitk::FilterBase(),
  45   clitk::FilterWithAnatomicalFeatureDatabaseManagement(),
  46   itk::ImageToImageFilter<ImageType, MaskImageType>()
  47 {
  48   SetNumberOfSteps(10);
  49   m_MaxSeedNumber = 500;
  50
  51   // Default global options
  52   this->SetNumberOfRequiredInputs(1);
  53   SetPatientMaskBackgroundValue(0);
  54   SetBackgroundValue(0); // Must be zero
  55   SetForegroundValue(1);
  56   SetMinimalComponentSize(100);
  57
  58   // Step 1 default values
  59   SetUpperThreshold(-300);
  60   SetLowerThreshold(-1000);
  61   UseLowerThresholdOff();
  62   LabelParamType * p1 = new LabelParamType;
  63   p1->SetFirstKeep(1);
  64   p1->UseLastKeepOff();
  65   p1->AddLabelToRemove(2);
  66   SetLabelizeParameters1(p1);
  67
  68   // Step 2 default values
  69   SetUpperThresholdForTrachea(-900);
  70   SetMultiplierForTrachea(5);
  71   SetThresholdStepSizeForTrachea(64);
  72   SetNumberOfSlicesToSkipBeforeSearchingSeed(0);
  73
  74   // Step 3 default values
  75   SetNumberOfHistogramBins(500);
  76   LabelParamType * p2 = new LabelParamType;
  77   p2->SetFirstKeep(1);
  78   p2->UseLastKeepOff();
  79   // p->AddLabelToRemove(2); // No label to remove by default
  80   SetLabelizeParameters2(p2);
  81
  82   // Step 4 default values
  83   SetRadiusForTrachea(1);
  84   LabelParamType * p3 = new LabelParamType;
  85   p3->SetFirstKeep(1);
  86   p3->SetLastKeep(2);
  87   p3->UseLastKeepOff();
  88   SetLabelizeParameters3(p3);
  89
  90   // Step 5
  91   OpenCloseOff();
  92   SetOpenCloseRadius(1);
  93
  94   // Step 6
  95   FillHolesOn();
  96 }
  97 //--------------------------------------------------------------------
  98
  99
 100 //--------------------------------------------------------------------
 101 template <class ImageType>
 102 void
 103 clitk::ExtractLungFilter<ImageType>::
 104 SetInput(const ImageType * image)
 105 {
 106   this->SetNthInput(0, const_cast<ImageType *>(image));
 107 }
 108 //--------------------------------------------------------------------
 109
 110
 111 //--------------------------------------------------------------------
 112 template <class ImageType>
 113 void
 114 clitk::ExtractLungFilter<ImageType>::
 115 AddSeed(InternalIndexType s)
 116 {
 117   m_Seeds.push_back(s);
 118 }
 119 //--------------------------------------------------------------------
 120
 121
 122 //--------------------------------------------------------------------
 123 template <class ImageType>
 124 template<class ArgsInfoType>
 125 void
 126 clitk::ExtractLungFilter<ImageType>::
 127 SetArgsInfo(ArgsInfoType mArgsInfo)
 128 {
 129   SetVerboseOption_GGO(mArgsInfo);
 130   SetVerboseStep_GGO(mArgsInfo);
 131   SetWriteStep_GGO(mArgsInfo);
 132   SetVerboseWarningOff_GGO(mArgsInfo);
 133
 134   SetAFDBFilename_GGO(mArgsInfo);
 135   SetOutputLungFilename_GGO(mArgsInfo);
 136   SetOutputTracheaFilename_GGO(mArgsInfo);
 137
 138   SetUpperThreshold_GGO(mArgsInfo);
 139   SetLowerThreshold_GGO(mArgsInfo);
 140   SetNumberOfSlicesToSkipBeforeSearchingSeed_GGO(mArgsInfo);
 141   SetLabelizeParameters1_GGO(mArgsInfo);
 142   if (!mArgsInfo.remove1_given) {
 143     GetLabelizeParameters1()->AddLabelToRemove(2);
 144     if (GetVerboseOption()) VerboseOption("remove1", 2);
 145   }
 146
 147   SetUpperThresholdForTrachea_GGO(mArgsInfo);
 148   SetMultiplierForTrachea_GGO(mArgsInfo);
 149   SetThresholdStepSizeForTrachea_GGO(mArgsInfo);
 150   AddSeed_GGO(mArgsInfo);
 151
 152   SetMinimalComponentSize_GGO(mArgsInfo);
 153
 154   SetNumberOfHistogramBins_GGO(mArgsInfo);
 155   SetLabelizeParameters2_GGO(mArgsInfo);
 156
 157   SetRadiusForTrachea_GGO(mArgsInfo);
 158   SetLabelizeParameters3_GGO(mArgsInfo);
 159
 160   SetOpenCloseRadius_GGO(mArgsInfo);
 161   SetOpenClose_GGO(mArgsInfo);
 162
 163   SetFillHoles_GGO(mArgsInfo);
 164 }
 165 //--------------------------------------------------------------------
 166
 167
 168 //--------------------------------------------------------------------
 169 template <class ImageType>
 170 void
 171 clitk::ExtractLungFilter<ImageType>::
 172 GenerateOutputInformation()
 173 {
 174   Superclass::GenerateOutputInformation();
 175
 176   // Read DB
 177   LoadAFDB();
 178
 179   // Get input pointers
 180   input   = dynamic_cast<const ImageType*>(itk::ProcessObject::GetInput(0));
 181   patient = GetAFDB()->template GetImage <MaskImageType>("patient");
 182
 183   //--------------------------------------------------------------------
 184   //--------------------------------------------------------------------
 185   // Crop input like patient image (must have the same spacing)
 186   StartNewStep("Crop input image to 'patient' extends");
 187   typedef clitk::CropLikeImageFilter<ImageType> CropImageFilter;
 188   typename CropImageFilter::Pointer cropFilter = CropImageFilter::New();
 189   cropFilter->SetInput(input);
 190   cropFilter->SetCropLikeImage(patient);
 191   cropFilter->Update();
 192   working_input = cropFilter->GetOutput();
 193   StopCurrentStep<ImageType>(working_input);
 194
 195   //--------------------------------------------------------------------
 196   //--------------------------------------------------------------------
 197   StartNewStep("Set background to initial image");
 198   working_input = SetBackground<ImageType, MaskImageType>
 199     (working_input, patient, GetPatientMaskBackgroundValue(), -1000);
 200   StopCurrentStep<ImageType>(working_input);
 201
 202   //--------------------------------------------------------------------
 203   //--------------------------------------------------------------------
 204   StartNewStep("Remove Air");
 205   // Check threshold
 206   if (m_UseLowerThreshold) {
 207     if (m_LowerThreshold > m_UpperThreshold) {
 208       clitkExceptionMacro("lower threshold cannot be greater than upper threshold.");
 209     }
 210   }
 211   // Threshold to get air
 212   typedef itk::BinaryThresholdImageFilter<ImageType, InternalImageType> InputBinarizeFilterType;
 213   typename InputBinarizeFilterType::Pointer binarizeFilter=InputBinarizeFilterType::New();
 214   binarizeFilter->SetInput(working_input);
 215   if (m_UseLowerThreshold) binarizeFilter->SetLowerThreshold(m_LowerThreshold);
 216   binarizeFilter->SetUpperThreshold(m_UpperThreshold);
 217   binarizeFilter ->SetInsideValue(this->GetForegroundValue());
 218   binarizeFilter ->SetOutsideValue(this->GetBackgroundValue());
 219   binarizeFilter->Update();
 220   working_image = binarizeFilter->GetOutput();
 221
 222   // Labelize and keep right labels
 223   working_image = Labelize<InternalImageType>(working_image, GetBackgroundValue(), true, GetMinimalComponentSize());
 224
 225   working_image = RemoveLabels<InternalImageType>
 226     (working_image, GetBackgroundValue(), GetLabelizeParameters1()->GetLabelsToRemove());
 227
 228   typename InternalImageType::Pointer air = KeepLabels<InternalImageType>
 229     (working_image,
 230      GetBackgroundValue(),
 231      GetForegroundValue(),
 232      GetLabelizeParameters1()->GetFirstKeep(),
 233      GetLabelizeParameters1()->GetLastKeep(),
 234      GetLabelizeParameters1()->GetUseLastKeep());
 235
 236   // Set Air to BG
 237   working_input = SetBackground<ImageType, InternalImageType>
 238     (working_input, air, this->GetForegroundValue(), this->GetBackgroundValue());
 239
 240   // End
 241   StopCurrentStep<ImageType>(working_input);
 242
 243   //--------------------------------------------------------------------
 244   //--------------------------------------------------------------------
 245   StartNewStep("Search for the trachea");
 246   SearchForTrachea();
 247
 248   //--------------------------------------------------------------------
 249   //--------------------------------------------------------------------
 250   StartNewStep("Extract the lung with Otsu filter");
 251   // Automated Otsu thresholding and relabeling
 252   typedef itk::OtsuThresholdImageFilter<ImageType,InternalImageType> OtsuThresholdImageFilterType;
 253   typename OtsuThresholdImageFilterType::Pointer otsuFilter=OtsuThresholdImageFilterType::New();
 254   otsuFilter->SetInput(working_input);
 255   otsuFilter->SetNumberOfHistogramBins(GetNumberOfHistogramBins());
 256   otsuFilter->SetInsideValue(this->GetForegroundValue());
 257   otsuFilter->SetOutsideValue(this->GetBackgroundValue());
 258   otsuFilter->Update();
 259   working_image =  otsuFilter->GetOutput();
 260
 261   // Set output
 262   StopCurrentStep<InternalImageType>(working_image);
 263
 264   //--------------------------------------------------------------------
 265   //--------------------------------------------------------------------
 266   StartNewStep("Select labels");
 267   // Keep right labels
 268   working_image = LabelizeAndSelectLabels<InternalImageType>
 269     (working_image,
 270      GetBackgroundValue(),
 271      GetForegroundValue(),
 272      false,
 273      GetMinimalComponentSize(),
 274      GetLabelizeParameters2());
 275
 276   // Set output
 277   StopCurrentStep<InternalImageType>(working_image);
 278
 279   //--------------------------------------------------------------------
 280   //--------------------------------------------------------------------
 281   if (m_Seeds.size() != 0) { // if ==0 ->no trachea found
 282     StartNewStep("Remove the trachea");
 283     // Set the trachea
 284     working_image = SetBackground<InternalImageType, InternalImageType>
 285       (working_image, trachea_tmp, 1, -1);
 286
 287     // Dilate the trachea
 288     static const unsigned int Dim = ImageType::ImageDimension;
 289     typedef itk::BinaryBallStructuringElement<InternalPixelType, Dim> KernelType;
 290     KernelType structuringElement;
 291     structuringElement.SetRadius(GetRadiusForTrachea());
 292     structuringElement.CreateStructuringElement();
 293     typedef clitk::ConditionalBinaryDilateImageFilter<InternalImageType, InternalImageType, KernelType> ConditionalBinaryDilateImageFilterType;
 294     typename ConditionalBinaryDilateImageFilterType::Pointer dilateFilter = ConditionalBinaryDilateImageFilterType::New();
 295     dilateFilter->SetBoundaryToForeground(false);
 296     dilateFilter->SetKernel(structuringElement);
 297     dilateFilter->SetBackgroundValue (1);
 298     dilateFilter->SetForegroundValue (-1);
 299     dilateFilter->SetInput (working_image);
 300     dilateFilter->Update();
 301     working_image = dilateFilter->GetOutput();
 302
 303     // Set trachea with dilatation
 304     trachea_tmp = SetBackground<InternalImageType, InternalImageType>
 305       (trachea_tmp, working_image, -1, this->GetForegroundValue());
 306
 307     // Remove the trachea
 308     working_image = SetBackground<InternalImageType, InternalImageType>
 309       (working_image, working_image, -1, this->GetBackgroundValue());
 310
 311     // Label
 312     working_image = LabelizeAndSelectLabels<InternalImageType>
 313       (working_image,
 314        GetBackgroundValue(),
 315        GetForegroundValue(),
 316        false,
 317        GetMinimalComponentSize(),
 318        GetLabelizeParameters3());
 319
 320     // Set output
 321     StopCurrentStep<InternalImageType>(working_image);
 322   }
 323
 324   //--------------------------------------------------------------------
 325   //--------------------------------------------------------------------
 326   typedef clitk::AutoCropFilter<InternalImageType> AutoCropFilterType;
 327   typename AutoCropFilterType::Pointer autocropFilter = AutoCropFilterType::New();
 328   if (m_Seeds.size() != 0) { // if ==0 ->no trachea found
 329     StartNewStep("Cropping trachea");
 330     autocropFilter->SetInput(trachea_tmp);
 331     autocropFilter->Update(); // Needed
 332     typedef itk::CastImageFilter<InternalImageType, MaskImageType> CastImageFilterType;
 333     typename CastImageFilterType::Pointer caster= CastImageFilterType::New();
 334     caster->SetInput(autocropFilter->GetOutput());
 335     caster->Update();
 336     trachea = caster->GetOutput();
 337     StopCurrentStep<MaskImageType>(trachea);
 338   }
 339
 340   //--------------------------------------------------------------------
 341   //--------------------------------------------------------------------
 342   StartNewStep("Cropping lung");
 343   typename AutoCropFilterType::Pointer autocropFilter2 = AutoCropFilterType::New(); // Needed to reset pipeline
 344   autocropFilter2->SetInput(working_image);
 345   autocropFilter2->Update();
 346   working_image = autocropFilter2->GetOutput();
 347   StopCurrentStep<InternalImageType>(working_image);
 348
 349   //--------------------------------------------------------------------
 350   //--------------------------------------------------------------------
 351   // Final OpenClose
 352   if (GetOpenClose()) {
 353     StartNewStep("Open/Close");
 354
 355     // Structuring element
 356     typedef itk::BinaryBallStructuringElement<InternalPixelType, ImageDimension> KernelType;
 357     KernelType structuringElement;
 358     structuringElement.SetRadius(GetOpenCloseRadius());
 359     structuringElement.CreateStructuringElement();
 360
 361     // Open
 362     typedef itk::BinaryMorphologicalOpeningImageFilter<InternalImageType, InternalImageType, KernelType> OpenFilterType;
 363     typename OpenFilterType::Pointer openFilter = OpenFilterType::New();
 364     openFilter->SetInput(working_image);
 365     openFilter->SetBackgroundValue(GetBackgroundValue());
 366     openFilter->SetForegroundValue(GetForegroundValue());
 367     openFilter->SetKernel(structuringElement);
 368
 369     // Close
 370     typedef itk::BinaryMorphologicalClosingImageFilter<InternalImageType, InternalImageType, KernelType> CloseFilterType;
 371     typename CloseFilterType::Pointer closeFilter = CloseFilterType::New();
 372     closeFilter->SetInput(openFilter->GetOutput());
 373     closeFilter->SetSafeBorder(true);
 374     closeFilter->SetForegroundValue(GetForegroundValue());
 375     closeFilter->SetKernel(structuringElement);
 376     closeFilter->Update();
 377     working_image = closeFilter->GetOutput();
 378   }
 379
 380   //--------------------------------------------------------------------
 381   //--------------------------------------------------------------------
 382   // Fill Lungs
 383   if (GetFillHoles()) {
 384     StartNewStep("Fill Holes");
 385     /*
 386     typename FillMaskFilterType::Pointer fillMaskFilter = FillMaskFilterType::New();
 387     fillMaskFilter(working_image);
 388     fillMaskFilter->Update();
 389     working_image = fillMaskFilter->GetOutput();
 390     StopCurrentStep<InternalImageType>(working_image);
 391     */
 392   }
 393
 394   //--------------------------------------------------------------------
 395   //--------------------------------------------------------------------
 396   StartNewStep("Separate Left/Right lungs");
 397   // Initial label
 398   working_image = Labelize<InternalImageType>(working_image,
 399                                               GetBackgroundValue(),
 400                                               false,
 401                                               GetMinimalComponentSize());
 402
 403   // Count the labels
 404   typedef itk::StatisticsImageFilter<InternalImageType> StatisticsImageFilterType;
 405   typename StatisticsImageFilterType::Pointer statisticsImageFilter=StatisticsImageFilterType::New();
 406   statisticsImageFilter->SetInput(working_image);
 407   statisticsImageFilter->Update();
 408   unsigned int initialNumberOfLabels = statisticsImageFilter->GetMaximum();
 409   working_image = statisticsImageFilter->GetOutput();
 410
 411   // Decompose the first label
 412   static const unsigned int Dim = ImageType::ImageDimension;
 413   if (initialNumberOfLabels<2) {
 414     // Structuring element radius
 415     typename ImageType::SizeType radius;
 416     for (unsigned int i=0;i<Dim;i++) radius[i]=1;
 417     typedef clitk::DecomposeAndReconstructImageFilter<InternalImageType,InternalImageType> DecomposeAndReconstructFilterType;
 418     typename DecomposeAndReconstructFilterType::Pointer decomposeAndReconstructFilter=DecomposeAndReconstructFilterType::New();
 419     decomposeAndReconstructFilter->SetInput(working_image);
 420     decomposeAndReconstructFilter->SetVerbose(false);
 421     decomposeAndReconstructFilter->SetRadius(radius);
 422     decomposeAndReconstructFilter->SetMaximumNumberOfLabels(2);
 423     decomposeAndReconstructFilter->SetMinimumObjectSize(this->GetMinimalComponentSize());
 424     decomposeAndReconstructFilter->SetMinimumNumberOfIterations(1);
 425     decomposeAndReconstructFilter->SetBackgroundValue(this->GetBackgroundValue());
 426     decomposeAndReconstructFilter->SetForegroundValue(this->GetForegroundValue());
 427     decomposeAndReconstructFilter->SetFullyConnected(true);
 428     decomposeAndReconstructFilter->SetNumberOfNewLabels(1);
 429     decomposeAndReconstructFilter->Update();
 430     working_image = decomposeAndReconstructFilter->GetOutput();
 431   }
 432
 433   // Retain labels ('1' is largset lung, so right. '2' is left)
 434   typedef itk::ThresholdImageFilter<InternalImageType> ThresholdImageFilterType;
 435   typename ThresholdImageFilterType::Pointer thresholdFilter = ThresholdImageFilterType::New();
 436   thresholdFilter->SetInput(working_image);
 437   thresholdFilter->ThresholdAbove(2);
 438   thresholdFilter->SetOutsideValue(this->GetBackgroundValue());
 439   thresholdFilter->Update();
 440   working_image = thresholdFilter->GetOutput();
 441   StopCurrentStep<InternalImageType> (working_image);
 442
 443   // Final Cast
 444   StartNewStep("Cast the lung mask");
 445   typedef itk::CastImageFilter<InternalImageType, MaskImageType> CastImageFilterType;
 446   typename CastImageFilterType::Pointer caster= CastImageFilterType::New();
 447   caster->SetInput(working_image);
 448   caster->Update();
 449   output = caster->GetOutput();
 450
 451   // Update output info
 452   this->GetOutput(0)->SetRegions(output->GetLargestPossibleRegion());
 453 }
 454 //--------------------------------------------------------------------
 455
 456
 457 //--------------------------------------------------------------------
 458 template <class ImageType>
 459 void
 460 clitk::ExtractLungFilter<ImageType>::
 461 GenerateData()
 462 {
 463   // Set the output
 464   this->GraftOutput(output); // not SetNthOutput
 465   // Store image filenames into AFDB
 466   GetAFDB()->SetImageFilename("lungs", this->GetOutputLungFilename());
 467   GetAFDB()->SetImageFilename("trachea", this->GetOutputTracheaFilename());
 468   WriteAFDB();
 469 }
 470 //--------------------------------------------------------------------
 471
 472
 473 //--------------------------------------------------------------------
 474 template <class ImageType>
 475 bool
 476 clitk::ExtractLungFilter<ImageType>::
 477 SearchForTracheaSeed(int skip)
 478 {
 479   if (m_Seeds.size() == 0) { // try to find seed (if not zero, it is given by user)
 480     // Restart the search until a seed is found, skipping more and more slices
 481     bool stop = false;
 482     while (!stop) {
 483       // Search seed (parameters = UpperThresholdForTrachea)
 484       static const unsigned int Dim = ImageType::ImageDimension;
 485       typename InternalImageType::RegionType sliceRegion = working_input->GetLargestPossibleRegion();
 486       typename InternalImageType::SizeType sliceRegionSize = sliceRegion.GetSize();
 487       typename InternalImageType::IndexType sliceRegionIndex = sliceRegion.GetIndex();
 488       sliceRegionIndex[Dim-1]=sliceRegionSize[Dim-1]-skip-5;
 489       sliceRegionSize[Dim-1]=5;
 490       sliceRegion.SetSize(sliceRegionSize);
 491       sliceRegion.SetIndex(sliceRegionIndex);
 492       typedef  itk::ImageRegionConstIterator<ImageType> IteratorType;
 493       IteratorType it(working_input, sliceRegion);
 494       it.GoToBegin();
 495       while (!it.IsAtEnd()) {
 496         if(it.Get() < GetUpperThresholdForTrachea() ) {
 497           AddSeed(it.GetIndex());
 498           // DD(it.GetIndex());
 499         }
 500         ++it;
 501       }
 502
 503       // if we do not found : restart
 504       stop = (m_Seeds.size() != 0);
 505       if (!stop) {
 506         if (GetVerboseStep()) {
 507           std::cout << "\t No seed found this time. I skip some slices and restart." << std::endl;
 508         }
 509         if (skip > 0.5 * working_input->GetLargestPossibleRegion().GetSize()[2]) {
 510           // we want to skip more than a half of the image, it is probably a bug
 511           std::cerr << "Number of slices to skip to find trachea to high = " << skip << std::endl;
 512           stop = true;
 513         }
 514         skip += 5;
 515       }
 516     }
 517   }
 518   return (m_Seeds.size() != 0);
 519 }
 520 //--------------------------------------------------------------------
 521
 522
 523 //--------------------------------------------------------------------
 524 template <class ImageType>
 525 void
 526 clitk::ExtractLungFilter<ImageType>::
 527 TracheaRegionGrowing()
 528 {
 529   // Explosion controlled region growing
 530   typedef clitk::ExplosionControlledThresholdConnectedImageFilter<ImageType, InternalImageType> ImageFilterType;
 531   typename ImageFilterType::Pointer f= ImageFilterType::New();
 532   f->SetInput(working_input);
 533   f->SetVerbose(false);
 534   f->SetLower(-2000);
 535   f->SetUpper(GetUpperThresholdForTrachea());
 536   f->SetMinimumLowerThreshold(-2000);
 537   f->SetMaximumUpperThreshold(0);
 538   f->SetAdaptLowerBorder(false);
 539   f->SetAdaptUpperBorder(true);
 540   f->SetMinimumSize(5000);
 541   f->SetReplaceValue(1);
 542   f->SetMultiplier(GetMultiplierForTrachea());
 543   f->SetThresholdStepSize(GetThresholdStepSizeForTrachea());
 544   f->SetMinimumThresholdStepSize(1);
 545   for(unsigned int i=0; i<m_Seeds.size();i++) {
 546     f->AddSeed(m_Seeds[i]);
 547     // DD(m_Seeds[i]);
 548   }
 549   f->Update();
 550
 551   // take first (main) connected component
 552   trachea_tmp = Labelize<InternalImageType>(f->GetOutput(),
 553                                             GetBackgroundValue(),
 554                                             true,
 555                                             GetMinimalComponentSize());
 556   trachea_tmp = KeepLabels<InternalImageType>(trachea_tmp,
 557                                               GetBackgroundValue(),
 558                                               GetForegroundValue(),
 559                                               1, 1, false);
 560 }
 561 //--------------------------------------------------------------------
 562
 563
 564 //--------------------------------------------------------------------
 565 template <class ImageType>
 566 double
 567 clitk::ExtractLungFilter<ImageType>::
 568 ComputeTracheaVolume()
 569 {
 570   typedef itk::ImageRegionConstIterator<InternalImageType> IteratorType;
 571   IteratorType iter(trachea_tmp, trachea_tmp->GetLargestPossibleRegion());
 572   iter.GoToBegin();
 573   double volume = 0.0;
 574   while (!iter.IsAtEnd()) {
 575     if (iter.Get() == this->GetForegroundValue()) volume++;
 576     ++iter;
 577   }
 578
 579   double voxelsize = trachea_tmp->GetSpacing()[0]*trachea_tmp->GetSpacing()[1]*trachea_tmp->GetSpacing()[2];
 580   return volume*voxelsize;
 581 }
 582 //--------------------------------------------------------------------
 583
 584
 585 //--------------------------------------------------------------------
 586 template <class ImageType>
 587 void
 588 clitk::ExtractLungFilter<ImageType>::
 589 SearchForTrachea()
 590 {
 591   // Search for seed among n slices, skip some slices before starting
 592   // if not found -> skip more and restart
 593
 594   // when seed found : perform region growing
 595   // compute trachea volume
 596   // if volume not plausible  -> skip more slices and restart
 597
 598   bool stop = false;
 599   double volume = 0.0;
 600   int skip = GetNumberOfSlicesToSkipBeforeSearchingSeed();
 601   while (!stop) {
 602     stop = SearchForTracheaSeed(skip);
 603     if (stop) {
 604       TracheaRegionGrowing();
 605       volume = ComputeTracheaVolume()/1000; // assume mm3, so divide by 1000 to get cc
 606       if ((volume > 10) && (volume < 55 )) { // it is ok
 607         // Typical volume 22.59 cm 3 (± 7.69 cm 3 ) [ Leader 2004 ]
 608         if (GetVerboseStep()) {
 609           std::cout << "\t Found trachea with volume " << volume << " cc." << std::endl;
 610         }
 611         stop = true;
 612       }
 613       else {
 614         if (GetVerboseStep()) {
 615           std::cout << "\t The volume of the trachea (" << volume
 616                     << " cc) seems not correct. I skip some slices (" << skip << ") and restart to find seeds."
 617                     << std::endl;
 618         }
 619         skip += 5;
 620         stop = false;
 621         // empty the list of seed
 622         m_Seeds.clear();
 623       }
 624     }
 625   }
 626
 627   if (volume != 0.0) {
 628     // Set output
 629     StopCurrentStep<InternalImageType>(trachea_tmp);
 630   }
 631   else { // Trachea not found
 632     this->SetWarning("* WARNING * No seed found for trachea.");
 633     StopCurrentStep();
 634   }
 635 }
 636 //--------------------------------------------------------------------
 637
 638 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX