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://www.centreleonberard.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 #include "clitkMemoryUsage.h"
  30
  31 // itk
  32 #include "itkBinaryThresholdImageFilter.h"
  33 #include "itkConnectedComponentImageFilter.h"
  34 #include "itkRelabelComponentImageFilter.h"
  35 #include "itkOtsuThresholdImageFilter.h"
  36 #include "itkBinaryThinningImageFilter3D.h"
  37 #include "itkImageIteratorWithIndex.h"
  38 #include "itkBinaryMorphologicalOpeningImageFilter.h"
  39 #include "itkBinaryMorphologicalClosingImageFilter.h"
  40 #include "itkConstantPadImageFilter.h"
  41
  42 //--------------------------------------------------------------------
  43 template <class ImageType>
  44 clitk::ExtractLungFilter<ImageType>::
  45 ExtractLungFilter():
  46   clitk::FilterBase(),
  47   clitk::FilterWithAnatomicalFeatureDatabaseManagement(),
  48   itk::ImageToImageFilter<ImageType, MaskImageType>()
  49 {
  50   SetNumberOfSteps(10);
  51   m_MaxSeedNumber = 500;
  52
  53   // Default global options
  54   this->SetNumberOfRequiredInputs(1);
  55   SetPatientMaskBackgroundValue(0);
  56   SetBackgroundValue(0); // Must be zero
  57   SetForegroundValue(1);
  58   SetMinimalComponentSize(100);
  59   VerboseRegionGrowingFlagOff();
  60
  61   // Step 1 default values
  62   SetUpperThreshold(-300);
  63   SetLowerThreshold(-1000);
  64   UseLowerThresholdOff();
  65   LabelParamType * p1 = new LabelParamType;
  66   p1->SetFirstKeep(1);
  67   p1->UseLastKeepOff();
  68   p1->AddLabelToRemove(2);
  69   SetLabelizeParameters1(p1);
  70
  71   // Step 2 default values
  72   SetUpperThresholdForTrachea(-900);
  73   SetMultiplierForTrachea(5);
  74   SetThresholdStepSizeForTrachea(64);
  75   SetNumberOfSlicesToSkipBeforeSearchingSeed(0);
  76   TracheaVolumeMustBeCheckedFlagOn();
  77
  78   // Step 3 default values
  79   SetNumberOfHistogramBins(500);
  80   LabelParamType * p2 = new LabelParamType;
  81   p2->SetFirstKeep(1);
  82   p2->UseLastKeepOff();
  83   // p->AddLabelToRemove(2); // No label to remove by default
  84   SetLabelizeParameters2(p2);
  85
  86   // Step 4 default values
  87   SetRadiusForTrachea(1);
  88   LabelParamType * p3 = new LabelParamType;
  89   p3->SetFirstKeep(1);
  90   p3->SetLastKeep(2);
  91   p3->UseLastKeepOff();
  92   SetLabelizeParameters3(p3);
  93
  94   // Step 5
  95   OpenCloseFlagOff();
  96   SetOpenCloseRadius(1);
  97
  98   // Step 6
  99   FillHolesFlagOn();
 100   AutoCropOn();
 101 }
 102 //--------------------------------------------------------------------
 103
 104
 105 //--------------------------------------------------------------------
 106 template <class ImageType>
 107 void
 108 clitk::ExtractLungFilter<ImageType>::
 109 SetInput(const ImageType * image)
 110 {
 111   this->SetNthInput(0, const_cast<ImageType *>(image));
 112 }
 113 //--------------------------------------------------------------------
 114
 115
 116 //--------------------------------------------------------------------
 117 template <class ImageType>
 118 void
 119 clitk::ExtractLungFilter<ImageType>::
 120 AddSeed(InternalIndexType s)
 121 {
 122   m_Seeds.push_back(s);
 123 }
 124 //--------------------------------------------------------------------
 125
 126
 127 //--------------------------------------------------------------------
 128 template <class ImageType>
 129 void
 130 clitk::ExtractLungFilter<ImageType>::
 131 GenerateOutputInformation()
 132 {
 133   clitk::PrintMemory(GetVerboseMemoryFlag(), "Initial memory"); // OK
 134   Superclass::GenerateOutputInformation();
 135
 136   // Read DB
 137   LoadAFDB();
 138
 139   // Get input pointers
 140   input   = dynamic_cast<const ImageType*>(itk::ProcessObject::GetInput(0));
 141   patient = GetAFDB()->template GetImage <MaskImageType>("Patient");
 142   PrintMemory(GetVerboseMemoryFlag(), "After reading patient"); // OK
 143
 144   //--------------------------------------------------------------------
 145   //--------------------------------------------------------------------
 146   // Crop input like patient image (must have the same spacing)
 147   // It copy the input if the same are the same
 148   StartNewStep("Copy and crop input image to 'patient' extends");
 149   typedef clitk::CropLikeImageFilter<ImageType> CropImageFilter;
 150   typename CropImageFilter::Pointer cropFilter = CropImageFilter::New();
 151   // cropFilter->ReleaseDataFlagOn(); // NO=seg fault !!
 152   cropFilter->SetInput(input);
 153   cropFilter->SetCropLikeImage(patient);
 154   cropFilter->Update();
 155   working_input = cropFilter->GetOutput();
 156   //  cropFilter->Delete(); // NO !!!! if yes, sg fault, Cropfilter is buggy !??
 157   StopCurrentStep<ImageType>(working_input);
 158   PrintMemory(GetVerboseMemoryFlag(), "After crop"); // OK, slightly more than a copy of input
 159
 160   //--------------------------------------------------------------------
 161   //--------------------------------------------------------------------
 162   StartNewStep("Set background to initial image");
 163   working_input = SetBackground<ImageType, MaskImageType>
 164     (working_input, patient, GetPatientMaskBackgroundValue(), -1000, true);
 165
 166   // Pad images with air to prevent patient touching the image border
 167   static const unsigned int Dim = ImageType::ImageDimension;
 168   typedef itk::ConstantPadImageFilter<ImageType, ImageType> PadFilterType;
 169   typename PadFilterType::Pointer padFilter = PadFilterType::New();
 170   padFilter->SetInput(working_input);
 171   padFilter->SetConstant(-1000);
 172   typename ImageType::SizeType bounds;
 173   for (unsigned i = 0; i < Dim - 1; ++i)
 174     bounds[i] = 1;
 175   bounds[Dim - 1] = 0;
 176   padFilter->SetPadLowerBound(bounds);
 177   padFilter->SetPadUpperBound(bounds);
 178   padFilter->Update();
 179   working_input = padFilter->GetOutput();
 180
 181   StopCurrentStep<ImageType>(working_input);
 182   PrintMemory(GetVerboseMemoryFlag(), "After set bg"); // OK, additional mem = 0
 183
 184   /*
 185   // We do not need patient mask anymore, release memory
 186   GetAFDB()->template ReleaseImage<MaskImageType>("Patient");
 187   DD(patient->GetReferenceCount());
 188   PrintMemory(GetVerboseMemoryFlag(), "After delete patient"); // OK, additional mem = 0
 189   patient->Delete();
 190   PrintMemory(GetVerboseMemoryFlag(), "After delete patient"); // OK, additional mem = 0
 191   */
 192
 193   //--------------------------------------------------------------------
 194   //--------------------------------------------------------------------
 195   StartNewStep("Remove Air");
 196   // Check threshold
 197   if (m_UseLowerThreshold) {
 198     if (m_LowerThreshold > m_UpperThreshold) {
 199       clitkExceptionMacro("lower threshold cannot be greater than upper threshold.");
 200     }
 201   }
 202   // Threshold to get air
 203   typedef itk::BinaryThresholdImageFilter<ImageType, MaskImageType> InputBinarizeFilterType;
 204   typename InputBinarizeFilterType::Pointer binarizeFilter = InputBinarizeFilterType::New();
 205   binarizeFilter->SetInput(working_input);
 206   if (m_UseLowerThreshold) binarizeFilter->SetLowerThreshold(m_LowerThreshold);
 207   // binarizeFilter->CanRunInPlace() is false
 208   binarizeFilter->SetUpperThreshold(m_UpperThreshold);
 209   binarizeFilter->SetInsideValue(this->GetForegroundValue());
 210   binarizeFilter->SetOutsideValue(this->GetBackgroundValue());
 211   binarizeFilter->Update();
 212   working_mask = binarizeFilter->GetOutput();
 213   PrintMemory(GetVerboseMemoryFlag(), "After Binarizefilter"); // OK, additional mem is one mask image
 214
 215   // Labelize and keep right labels
 216   working_mask =
 217     Labelize<MaskImageType>
 218     (working_mask, GetBackgroundValue(), true, GetMinimalComponentSize());
 219   PrintMemory(GetVerboseMemoryFlag(), "After Labelize"); // BUG ? additional mem around 1 time the input ?
 220
 221   working_mask = RemoveLabels<MaskImageType>
 222     (working_mask, GetBackgroundValue(), GetLabelizeParameters1()->GetLabelsToRemove());
 223   PrintMemory(GetVerboseMemoryFlag(), "After RemoveLabels"); // OK additional mem = 0
 224
 225   working_mask = KeepLabels<MaskImageType>
 226     (working_mask,
 227      GetBackgroundValue(),
 228      GetForegroundValue(),
 229      GetLabelizeParameters1()->GetFirstKeep(),
 230      GetLabelizeParameters1()->GetLastKeep(),
 231      GetLabelizeParameters1()->GetUseLastKeep());
 232   PrintMemory(GetVerboseMemoryFlag(), "After KeepLabels to create the 'air'");
 233
 234   // Set Air to BG
 235   working_input = SetBackground<ImageType, MaskImageType>
 236     (working_input, working_mask, this->GetForegroundValue(), this->GetBackgroundValue(), true);
 237   PrintMemory(GetVerboseMemoryFlag(), "After SetBackground");
 238
 239   // End
 240   StopCurrentStep<ImageType>(working_input);
 241
 242   //--------------------------------------------------------------------
 243   //--------------------------------------------------------------------
 244   StartNewStep("Search for the trachea");
 245   SearchForTrachea();
 246   PrintMemory(GetVerboseMemoryFlag(), "After SearchForTrachea");
 247
 248   //--------------------------------------------------------------------
 249   //--------------------------------------------------------------------
 250   StartNewStep("Extract the lung with Otsu filter");
 251   // Automated Otsu thresholding and relabeling
 252   typedef itk::OtsuThresholdImageFilter<ImageType,MaskImageType> 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_mask =  otsuFilter->GetOutput();
 260
 261   // Set output
 262   StopCurrentStep<MaskImageType>(working_mask);
 263   PrintMemory(GetVerboseMemoryFlag(), "After Otsufilter");
 264
 265   //--------------------------------------------------------------------
 266   //--------------------------------------------------------------------
 267   StartNewStep("Select labels");
 268   // Keep right labels
 269   working_mask = LabelizeAndSelectLabels<MaskImageType>
 270     (working_mask,
 271      GetBackgroundValue(),
 272      GetForegroundValue(),
 273      false,
 274      GetMinimalComponentSize(),
 275      GetLabelizeParameters2());
 276
 277   // Set output
 278   StopCurrentStep<MaskImageType>(working_mask);
 279   PrintMemory(GetVerboseMemoryFlag(), "After LabelizeAndSelectLabels");
 280
 281   //--------------------------------------------------------------------
 282   //--------------------------------------------------------------------
 283   if (m_Seeds.size() != 0) { // if ==0 ->no trachea found
 284     StartNewStep("Remove the trachea");
 285     // Set the trachea
 286     working_mask = SetBackground<MaskImageType, MaskImageType>
 287       (working_mask, trachea, 1, -1, true);
 288     PrintMemory(GetVerboseMemoryFlag(), "After SetBackground");
 289
 290     // Dilate the trachea
 291     static const unsigned int Dim = ImageType::ImageDimension;
 292     typedef itk::BinaryBallStructuringElement<InternalPixelType, Dim> KernelType;
 293     KernelType structuringElement;
 294     structuringElement.SetRadius(GetRadiusForTrachea());
 295     structuringElement.CreateStructuringElement();
 296     typedef clitk::ConditionalBinaryDilateImageFilter<MaskImageType, MaskImageType, KernelType> ConditionalBinaryDilateImageFilterType;
 297     typename ConditionalBinaryDilateImageFilterType::Pointer dilateFilter = ConditionalBinaryDilateImageFilterType::New();
 298     dilateFilter->SetBoundaryToForeground(false);
 299     dilateFilter->SetKernel(structuringElement);
 300     dilateFilter->SetBackgroundValue (1);
 301     dilateFilter->SetForegroundValue (-1);
 302     dilateFilter->SetInput (working_mask);
 303     dilateFilter->Update();
 304     working_mask = dilateFilter->GetOutput();
 305     PrintMemory(GetVerboseMemoryFlag(), "After dilate");
 306
 307     // Set trachea with dilatation
 308     trachea = SetBackground<MaskImageType, MaskImageType>
 309       (trachea, working_mask, -1, this->GetForegroundValue(), true);
 310
 311     // Remove the trachea
 312     working_mask = SetBackground<MaskImageType, MaskImageType>
 313       (working_mask, working_mask, -1, this->GetBackgroundValue(), true);
 314
 315     // Label
 316     working_mask = LabelizeAndSelectLabels<MaskImageType>
 317       (working_mask,
 318        GetBackgroundValue(),
 319        GetForegroundValue(),
 320        false,
 321        GetMinimalComponentSize(),
 322        GetLabelizeParameters3());
 323
 324     // Set output
 325     StopCurrentStep<MaskImageType>(working_mask);
 326   }
 327
 328   //--------------------------------------------------------------------
 329   //--------------------------------------------------------------------
 330   PrintMemory(GetVerboseMemoryFlag(), "before autocropfilter");
 331   if (m_Seeds.size() != 0) { // if ==0 ->no trachea found
 332     if (GetAutoCrop())
 333       trachea = clitk::AutoCrop<MaskImageType>(trachea, GetBackgroundValue());
 334     else
 335     {
 336       // Remove Padding region
 337       typedef itk::CropImageFilter<MaskImageType, MaskImageType> CropFilterType;
 338       typename CropFilterType::Pointer cropFilter = CropFilterType::New();
 339       cropFilter->SetInput(trachea);
 340       cropFilter->SetLowerBoundaryCropSize(bounds);
 341       cropFilter->SetUpperBoundaryCropSize(bounds);
 342       cropFilter->Update();
 343       trachea = cropFilter->GetOutput();
 344     }
 345     StopCurrentStep<MaskImageType>(trachea);
 346     PrintMemory(GetVerboseMemoryFlag(), "after delete trachea");
 347   }
 348   PrintMemory(GetVerboseMemoryFlag(), "after delete trachea");
 349
 350   //--------------------------------------------------------------------
 351   //--------------------------------------------------------------------
 352   StartNewStep("Cropping lung");
 353   PrintMemory(GetVerboseMemoryFlag(), "Before Autocropfilter");
 354   if (GetAutoCrop())
 355     working_mask = clitk::AutoCrop<MaskImageType>(working_mask, GetBackgroundValue());
 356   else
 357   {
 358     // Remove Padding region
 359     typedef itk::CropImageFilter<MaskImageType, MaskImageType> CropFilterType;
 360     typename CropFilterType::Pointer cropFilter = CropFilterType::New();
 361     cropFilter->SetInput(working_mask);
 362     cropFilter->SetLowerBoundaryCropSize(bounds);
 363     cropFilter->SetUpperBoundaryCropSize(bounds);
 364     cropFilter->Update();
 365     working_mask = cropFilter->GetOutput();
 366   }
 367   StopCurrentStep<MaskImageType>(working_mask);
 368
 369   //--------------------------------------------------------------------
 370   //--------------------------------------------------------------------
 371   // Final OpenClose
 372   if (GetOpenCloseFlag()) {
 373     StartNewStep("Open/Close");
 374     PrintMemory(GetVerboseMemoryFlag(), "Before OpenClose");
 375
 376     // Structuring element
 377     typedef itk::BinaryBallStructuringElement<InternalPixelType, ImageDimension> KernelType;
 378     KernelType structuringElement;
 379     structuringElement.SetRadius(GetOpenCloseRadius());
 380     structuringElement.CreateStructuringElement();
 381
 382     // Open
 383     typedef itk::BinaryMorphologicalOpeningImageFilter<MaskImageType, InternalImageType, KernelType> OpenFilterType;
 384     typename OpenFilterType::Pointer openFilter = OpenFilterType::New();
 385     openFilter->SetInput(working_mask);
 386     openFilter->SetBackgroundValue(GetBackgroundValue());
 387     openFilter->SetForegroundValue(GetForegroundValue());
 388     openFilter->SetKernel(structuringElement);
 389
 390     // Close
 391     typedef itk::BinaryMorphologicalClosingImageFilter<MaskImageType, MaskImageType, KernelType> CloseFilterType;
 392     typename CloseFilterType::Pointer closeFilter = CloseFilterType::New();
 393     closeFilter->SetInput(openFilter->GetOutput());
 394     closeFilter->SetSafeBorder(true);
 395     closeFilter->SetForegroundValue(GetForegroundValue());
 396     closeFilter->SetKernel(structuringElement);
 397     closeFilter->Update();
 398     working_mask = closeFilter->GetOutput();
 399   }
 400
 401   //--------------------------------------------------------------------
 402   //--------------------------------------------------------------------
 403   // Fill Lungs
 404   if (GetFillHolesFlag()) {
 405     StartNewStep("Fill Holes");
 406     PrintMemory(GetVerboseMemoryFlag(), "Before Fill Holes");
 407     typedef clitk::FillMaskFilter<MaskImageType> FillMaskFilterType;
 408     typename FillMaskFilterType::Pointer fillMaskFilter = FillMaskFilterType::New();
 409     fillMaskFilter->SetInput(working_mask);
 410     fillMaskFilter->AddDirection(2);
 411     //fillMaskFilter->AddDirection(1);
 412     fillMaskFilter->Update();
 413     working_mask = fillMaskFilter->GetOutput();
 414     StopCurrentStep<MaskImageType>(working_mask);
 415   }
 416
 417   //--------------------------------------------------------------------
 418   //--------------------------------------------------------------------
 419   StartNewStep("Separate Left/Right lungs");
 420     PrintMemory(GetVerboseMemoryFlag(), "Before Separate");
 421   // Initial label
 422   working_mask = Labelize<MaskImageType>(working_mask,
 423                                          GetBackgroundValue(),
 424                                          false,
 425                                          GetMinimalComponentSize());
 426
 427   PrintMemory(GetVerboseMemoryFlag(), "After Labelize");
 428
 429   // Count the labels
 430   typedef itk::StatisticsImageFilter<MaskImageType> StatisticsImageFilterType;
 431   typename StatisticsImageFilterType::Pointer statisticsImageFilter=StatisticsImageFilterType::New();
 432   statisticsImageFilter->SetInput(working_mask);
 433   statisticsImageFilter->Update();
 434   unsigned int initialNumberOfLabels = statisticsImageFilter->GetMaximum();
 435   working_mask = statisticsImageFilter->GetOutput();
 436
 437   PrintMemory(GetVerboseMemoryFlag(), "After count label");
 438
 439   // Decompose the first label
 440   if (initialNumberOfLabels<2) {
 441     // Structuring element radius
 442     typename ImageType::SizeType radius;
 443     for (unsigned int i=0;i<Dim;i++) radius[i]=1;
 444     typedef clitk::DecomposeAndReconstructImageFilter<MaskImageType,MaskImageType> DecomposeAndReconstructFilterType;
 445     typename DecomposeAndReconstructFilterType::Pointer decomposeAndReconstructFilter=DecomposeAndReconstructFilterType::New();
 446     decomposeAndReconstructFilter->SetInput(working_mask);
 447     decomposeAndReconstructFilter->SetVerbose(false);
 448     decomposeAndReconstructFilter->SetRadius(radius);
 449     decomposeAndReconstructFilter->SetMaximumNumberOfLabels(2);
 450     decomposeAndReconstructFilter->SetMinimumObjectSize(this->GetMinimalComponentSize());
 451     decomposeAndReconstructFilter->SetMinimumNumberOfIterations(1);
 452     decomposeAndReconstructFilter->SetBackgroundValue(this->GetBackgroundValue());
 453     decomposeAndReconstructFilter->SetForegroundValue(this->GetForegroundValue());
 454     decomposeAndReconstructFilter->SetFullyConnected(true);
 455     decomposeAndReconstructFilter->SetNumberOfNewLabels(1);
 456     decomposeAndReconstructFilter->Update();
 457     working_mask = decomposeAndReconstructFilter->GetOutput();
 458   }
 459   PrintMemory(GetVerboseMemoryFlag(), "After decomposeAndReconstructFilter");
 460
 461   // Retain labels ('1' is largset lung, so right. '2' is left)
 462   typedef itk::ThresholdImageFilter<MaskImageType> ThresholdImageFilterType;
 463   typename ThresholdImageFilterType::Pointer thresholdFilter = ThresholdImageFilterType::New();
 464   thresholdFilter->SetInput(working_mask);
 465   thresholdFilter->ThresholdAbove(2);
 466   thresholdFilter->SetOutsideValue(this->GetBackgroundValue());
 467   thresholdFilter->Update();
 468   working_mask = thresholdFilter->GetOutput();
 469   StopCurrentStep<MaskImageType> (working_mask);
 470   PrintMemory(GetVerboseMemoryFlag(), "After Thresholdfilter");
 471
 472   // Update output info
 473   //  output = working_mask;
 474   //this->GetOutput(0)->SetRegions(output->GetLargestPossibleRegion());
 475
 476   //  this->GetOutput(0)->SetRegions(working_mask->GetLargestPossibleRegion());
 477
 478 }
 479 //--------------------------------------------------------------------
 480
 481
 482 //--------------------------------------------------------------------
 483 template <class TImageType>
 484 void
 485 clitk::ExtractLungFilter<TImageType>::
 486 GenerateInputRequestedRegion() {
 487   //  DD("GenerateInputRequestedRegion (nothing?)");
 488 }
 489 //--------------------------------------------------------------------
 490
 491
 492 //--------------------------------------------------------------------
 493 template <class ImageType>
 494 void
 495 clitk::ExtractLungFilter<ImageType>::
 496 GenerateData()
 497 {
 498   // Set the output
 499   //  this->GraftOutput(output); // not SetNthOutput
 500   this->GraftOutput(working_mask); // not SetNthOutput
 501   // Store image filenames into AFDB
 502   GetAFDB()->SetImageFilename("Lungs", this->GetOutputLungFilename());
 503   GetAFDB()->SetImageFilename("Trachea", this->GetOutputTracheaFilename());
 504   WriteAFDB();
 505 }
 506 //--------------------------------------------------------------------
 507
 508
 509 //--------------------------------------------------------------------
 510 template <class ImageType>
 511 bool
 512 clitk::ExtractLungFilter<ImageType>::
 513 SearchForTracheaSeed(int skip)
 514 {
 515   if (m_Seeds.size() == 0) { // try to find seed (if not zero, it is given by user)
 516     // Restart the search until a seed is found, skipping more and more slices
 517     bool stop = false;
 518     while (!stop) {
 519       // Search seed (parameters = UpperThresholdForTrachea)
 520       static const unsigned int Dim = ImageType::ImageDimension;
 521       typename InternalImageType::RegionType sliceRegion = working_input->GetLargestPossibleRegion();
 522       typename InternalImageType::SizeType sliceRegionSize = sliceRegion.GetSize();
 523       typename InternalImageType::IndexType sliceRegionIndex = sliceRegion.GetIndex();
 524       sliceRegionIndex[Dim-1]=sliceRegionSize[Dim-1]-skip-5;
 525       sliceRegionSize[Dim-1]=5;
 526       sliceRegion.SetSize(sliceRegionSize);
 527       sliceRegion.SetIndex(sliceRegionIndex);
 528       typedef  itk::ImageRegionConstIterator<ImageType> IteratorType;
 529       IteratorType it(working_input, sliceRegion);
 530       it.GoToBegin();
 531       while (!it.IsAtEnd()) {
 532         if(it.Get() < GetUpperThresholdForTrachea() ) {
 533           AddSeed(it.GetIndex());
 534           // DD(it.GetIndex());
 535         }
 536         ++it;
 537       }
 538
 539       // if we do not found : restart
 540       stop = (m_Seeds.size() != 0);
 541       if (!stop) {
 542         if (GetVerboseStepFlag()) {
 543           std::cout << "\t No seed found this time. I skip some slices and restart." << std::endl;
 544         }
 545         if (skip > 0.5 * working_input->GetLargestPossibleRegion().GetSize()[2]) {
 546           // we want to skip more than a half of the image, it is probably a bug
 547           std::cerr << "Number of slices to skip to find trachea to high = " << skip << std::endl;
 548           stop = true;
 549         }
 550         skip += 5;
 551       }
 552       else {
 553         // DD(m_Seeds[0]);
 554         // DD(m_Seeds.size());
 555       }
 556     }
 557   }
 558   return (m_Seeds.size() != 0);
 559 }
 560 //--------------------------------------------------------------------
 561
 562
 563 //--------------------------------------------------------------------
 564 template <class ImageType>
 565 void
 566 clitk::ExtractLungFilter<ImageType>::
 567 TracheaRegionGrowing()
 568 {
 569   // Explosion controlled region growing
 570   PrintMemory(GetVerboseMemoryFlag(), "Before ExplosionControlledThresholdConnectedImageFilter");
 571   typedef clitk::ExplosionControlledThresholdConnectedImageFilter<ImageType, MaskImageType> ImageFilterType;
 572   typename ImageFilterType::Pointer f= ImageFilterType::New();
 573   f->SetInput(working_input);
 574   f->SetLower(-2000);
 575   f->SetUpper(GetUpperThresholdForTrachea());
 576   f->SetMinimumLowerThreshold(-2000);
 577   //  f->SetMaximumUpperThreshold(0); // MAYBE TO CHANGE ???
 578   f->SetMaximumUpperThreshold(-800); // MAYBE TO CHANGE ???
 579   f->SetAdaptLowerBorder(false);
 580   f->SetAdaptUpperBorder(true);
 581   f->SetMinimumSize(5000);
 582   f->SetReplaceValue(1);
 583   f->SetMultiplier(GetMultiplierForTrachea());
 584   f->SetThresholdStepSize(GetThresholdStepSizeForTrachea());
 585   f->SetMinimumThresholdStepSize(1);
 586   f->SetVerbose(GetVerboseRegionGrowingFlag());
 587   for(unsigned int i=0; i<m_Seeds.size();i++) {
 588     f->AddSeed(m_Seeds[i]);
 589     // DD(m_Seeds[i]);
 590   }
 591   f->Update();
 592   PrintMemory(GetVerboseMemoryFlag(), "After RG update");
 593
 594   // take first (main) connected component
 595   trachea = Labelize<MaskImageType>(f->GetOutput(),
 596                                     GetBackgroundValue(),
 597                                     true,
 598                                     1000);//GetMinimalComponentSize());
 599   PrintMemory(GetVerboseMemoryFlag(), "After Labelize");
 600   trachea = KeepLabels<MaskImageType>(trachea,
 601                                               GetBackgroundValue(),
 602                                               GetForegroundValue(),
 603                                               1, 1, false);
 604   PrintMemory(GetVerboseMemoryFlag(), "After KeepLabels");
 605 }
 606 //--------------------------------------------------------------------
 607
 608
 609 //--------------------------------------------------------------------
 610 template <class ImageType>
 611 double
 612 clitk::ExtractLungFilter<ImageType>::
 613 ComputeTracheaVolume()
 614 {
 615   typedef itk::ImageRegionConstIterator<InternalImageType> IteratorType;
 616   IteratorType iter(trachea, trachea->GetLargestPossibleRegion());
 617   iter.GoToBegin();
 618   double volume = 0.0;
 619   while (!iter.IsAtEnd()) {
 620     if (iter.Get() == this->GetForegroundValue()) volume++;
 621     ++iter;
 622   }
 623
 624   double voxelsize = trachea->GetSpacing()[0]*trachea->GetSpacing()[1]*trachea->GetSpacing()[2];
 625   return volume*voxelsize;
 626 }
 627 //--------------------------------------------------------------------
 628
 629
 630 //--------------------------------------------------------------------
 631 template <class ImageType>
 632 void
 633 clitk::ExtractLungFilter<ImageType>::
 634 SearchForTrachea()
 635 {
 636   // Search for seed among n slices, skip some slices before starting
 637   // if not found -> skip more and restart
 638
 639   // when seed found : perform region growing
 640   // compute trachea volume
 641   // if volume not plausible  -> skip more slices and restart
 642
 643   bool stop = false;
 644   double volume = 0.0;
 645   int skip = GetNumberOfSlicesToSkipBeforeSearchingSeed();
 646   while (!stop) {
 647     stop = true;
 648     if (SearchForTracheaSeed(skip)) {
 649       TracheaRegionGrowing();
 650       volume = ComputeTracheaVolume()/1000; // assume mm3, so divide by 1000 to get cc
 651       if (GetWriteStepFlag()) {
 652         writeImage<MaskImageType>(trachea, "step-trachea-"+toString(skip)+".mhd");
 653       }
 654       if (GetTracheaVolumeMustBeCheckedFlag()) {
 655         if ((volume > 10) && (volume < 65 )) { // depend on image size ...
 656           // Typical volume 22.59 cm 3 (± 7.69 cm 3 ) [ Leader 2004 ]
 657           if (GetVerboseStepFlag()) {
 658             std::cout << "\t Found trachea with volume " << volume << " cc." << std::endl;
 659           }
 660         }
 661         else {
 662           if (GetVerboseStepFlag()) {
 663             std::cout << "\t The volume of the trachea (" << volume
 664                       << " cc) seems not correct. I skip some slices (" << skip << ") and restart to find seeds."
 665                       << std::endl;
 666           }
 667           skip += 5;
 668           stop = false;
 669           // empty the list of seed
 670           m_Seeds.clear();
 671         }
 672         if (skip > 0.5 * working_input->GetLargestPossibleRegion().GetSize()[2]) {
 673           // we want to skip more than a half of the image, it is probably a bug
 674           std::cerr << "2 : Number of slices to skip to find trachea too high = " << skip << std::endl;
 675           stop = true;
 676         }
 677       }
 678       else {
 679         stop = true;
 680       }
 681     }
 682   }
 683
 684   if (volume != 0.0) {
 685     // Set output
 686     StopCurrentStep<MaskImageType>(trachea);
 687   }
 688   else { // Trachea not found
 689     this->SetWarning("* WARNING * No seed found for trachea.");
 690     StopCurrentStep();
 691   }
 692 }
 693 //--------------------------------------------------------------------
 694
 695 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX