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