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
  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   PrintMemory(GetVerboseMemoryFlag(), "After count label");
 396
 397   // Decompose the first label
 398   static const unsigned int Dim = ImageType::ImageDimension;
 399   if (initialNumberOfLabels<2) {
 400     // Structuring element radius
 401     typename ImageType::SizeType radius;
 402     for (unsigned int i=0;i<Dim;i++) radius[i]=1;
 403     typedef clitk::DecomposeAndReconstructImageFilter<MaskImageType,MaskImageType> DecomposeAndReconstructFilterType;
 404     typename DecomposeAndReconstructFilterType::Pointer decomposeAndReconstructFilter=DecomposeAndReconstructFilterType::New();
 405     decomposeAndReconstructFilter->SetInput(working_mask);
 406     decomposeAndReconstructFilter->SetVerbose(false);
 407     decomposeAndReconstructFilter->SetRadius(radius);
 408     decomposeAndReconstructFilter->SetMaximumNumberOfLabels(2);
 409     decomposeAndReconstructFilter->SetMinimumObjectSize(this->GetMinimalComponentSize());
 410     decomposeAndReconstructFilter->SetMinimumNumberOfIterations(1);
 411     decomposeAndReconstructFilter->SetBackgroundValue(this->GetBackgroundValue());
 412     decomposeAndReconstructFilter->SetForegroundValue(this->GetForegroundValue());
 413     decomposeAndReconstructFilter->SetFullyConnected(true);
 414     decomposeAndReconstructFilter->SetNumberOfNewLabels(1);
 415     decomposeAndReconstructFilter->Update();
 416     working_mask = decomposeAndReconstructFilter->GetOutput();
 417   }
 418   PrintMemory(GetVerboseMemoryFlag(), "After decomposeAndReconstructFilter");
 419
 420   // Retain labels ('1' is largset lung, so right. '2' is left)
 421   typedef itk::ThresholdImageFilter<MaskImageType> ThresholdImageFilterType;
 422   typename ThresholdImageFilterType::Pointer thresholdFilter = ThresholdImageFilterType::New();
 423   thresholdFilter->SetInput(working_mask);
 424   thresholdFilter->ThresholdAbove(2);
 425   thresholdFilter->SetOutsideValue(this->GetBackgroundValue());
 426   thresholdFilter->Update();
 427   working_mask = thresholdFilter->GetOutput();
 428   StopCurrentStep<MaskImageType> (working_mask);
 429   PrintMemory(GetVerboseMemoryFlag(), "After Thresholdfilter");
 430
 431   // Update output info
 432   //  output = working_mask;
 433   //this->GetOutput(0)->SetRegions(output->GetLargestPossibleRegion());
 434
 435   //  this->GetOutput(0)->SetRegions(working_mask->GetLargestPossibleRegion());
 436
 437 }
 438 //--------------------------------------------------------------------
 439
 440
 441 //--------------------------------------------------------------------
 442 template <class TImageType>
 443 void
 444 clitk::ExtractLungFilter<TImageType>::
 445 GenerateInputRequestedRegion() {
 446   //  DD("GenerateInputRequestedRegion (nothing?)");
 447 }
 448 //--------------------------------------------------------------------
 449
 450
 451 //--------------------------------------------------------------------
 452 template <class ImageType>
 453 void
 454 clitk::ExtractLungFilter<ImageType>::
 455 GenerateData()
 456 {
 457   // Set the output
 458   //  this->GraftOutput(output); // not SetNthOutput
 459   this->GraftOutput(working_mask); // not SetNthOutput
 460   // Store image filenames into AFDB
 461   GetAFDB()->SetImageFilename("Lungs", this->GetOutputLungFilename());
 462   GetAFDB()->SetImageFilename("Trachea", this->GetOutputTracheaFilename());
 463   WriteAFDB();
 464 }
 465 //--------------------------------------------------------------------
 466
 467
 468 //--------------------------------------------------------------------
 469 template <class ImageType>
 470 bool
 471 clitk::ExtractLungFilter<ImageType>::
 472 SearchForTracheaSeed(int skip)
 473 {
 474   if (m_Seeds.size() == 0) { // try to find seed (if not zero, it is given by user)
 475     // Restart the search until a seed is found, skipping more and more slices
 476     bool stop = false;
 477     while (!stop) {
 478       // Search seed (parameters = UpperThresholdForTrachea)
 479       static const unsigned int Dim = ImageType::ImageDimension;
 480       typename InternalImageType::RegionType sliceRegion = working_input->GetLargestPossibleRegion();
 481       typename InternalImageType::SizeType sliceRegionSize = sliceRegion.GetSize();
 482       typename InternalImageType::IndexType sliceRegionIndex = sliceRegion.GetIndex();
 483       sliceRegionIndex[Dim-1]=sliceRegionSize[Dim-1]-skip-5;
 484       sliceRegionSize[Dim-1]=5;
 485       sliceRegion.SetSize(sliceRegionSize);
 486       sliceRegion.SetIndex(sliceRegionIndex);
 487       typedef  itk::ImageRegionConstIterator<ImageType> IteratorType;
 488       IteratorType it(working_input, sliceRegion);
 489       it.GoToBegin();
 490       while (!it.IsAtEnd()) {
 491         if(it.Get() < GetUpperThresholdForTrachea() ) {
 492           AddSeed(it.GetIndex());
 493           // DD(it.GetIndex());
 494         }
 495         ++it;
 496       }
 497
 498       // if we do not found : restart
 499       stop = (m_Seeds.size() != 0);
 500       if (!stop) {
 501         if (GetVerboseStepFlag()) {
 502           std::cout << "\t No seed found this time. I skip some slices and restart." << std::endl;
 503         }
 504         if (skip > 0.5 * working_input->GetLargestPossibleRegion().GetSize()[2]) {
 505           // we want to skip more than a half of the image, it is probably a bug
 506           std::cerr << "Number of slices to skip to find trachea to high = " << skip << std::endl;
 507           stop = true;
 508         }
 509         skip += 5;
 510       }
 511       else {
 512         // DD(m_Seeds[0]);
 513         // DD(m_Seeds.size());
 514       }
 515     }
 516   }
 517   return (m_Seeds.size() != 0);
 518 }
 519 //--------------------------------------------------------------------
 520
 521
 522 //--------------------------------------------------------------------
 523 template <class ImageType>
 524 void
 525 clitk::ExtractLungFilter<ImageType>::
 526 TracheaRegionGrowing()
 527 {
 528   // Explosion controlled region growing
 529   PrintMemory(GetVerboseMemoryFlag(), "Before ExplosionControlledThresholdConnectedImageFilter");
 530   typedef clitk::ExplosionControlledThresholdConnectedImageFilter<ImageType, MaskImageType> ImageFilterType;
 531   typename ImageFilterType::Pointer f= ImageFilterType::New();
 532   f->SetInput(working_input);
 533   f->SetLower(-2000);
 534   f->SetUpper(GetUpperThresholdForTrachea());
 535   f->SetMinimumLowerThreshold(-2000);
 536   //  f->SetMaximumUpperThreshold(0); // MAYBE TO CHANGE ???
 537   f->SetMaximumUpperThreshold(-800); // MAYBE TO CHANGE ???
 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   f->SetVerbose(GetVerboseRegionGrowingFlag());
 546   for(unsigned int i=0; i<m_Seeds.size();i++) {
 547     f->AddSeed(m_Seeds[i]);
 548     // DD(m_Seeds[i]);
 549   }
 550   f->Update();
 551   PrintMemory(GetVerboseMemoryFlag(), "After RG update");
 552
 553   // take first (main) connected component
 554   trachea = Labelize<MaskImageType>(f->GetOutput(),
 555                                     GetBackgroundValue(),
 556                                     true,
 557                                     1000);//GetMinimalComponentSize());
 558   PrintMemory(GetVerboseMemoryFlag(), "After Labelize");
 559   trachea = KeepLabels<MaskImageType>(trachea,
 560                                               GetBackgroundValue(),
 561                                               GetForegroundValue(),
 562                                               1, 1, false);
 563   PrintMemory(GetVerboseMemoryFlag(), "After KeepLabels");
 564 }
 565 //--------------------------------------------------------------------
 566
 567
 568 //--------------------------------------------------------------------
 569 template <class ImageType>
 570 double
 571 clitk::ExtractLungFilter<ImageType>::
 572 ComputeTracheaVolume()
 573 {
 574   typedef itk::ImageRegionConstIterator<InternalImageType> IteratorType;
 575   IteratorType iter(trachea, trachea->GetLargestPossibleRegion());
 576   iter.GoToBegin();
 577   double volume = 0.0;
 578   while (!iter.IsAtEnd()) {
 579     if (iter.Get() == this->GetForegroundValue()) volume++;
 580     ++iter;
 581   }
 582
 583   double voxelsize = trachea->GetSpacing()[0]*trachea->GetSpacing()[1]*trachea->GetSpacing()[2];
 584   return volume*voxelsize;
 585 }
 586 //--------------------------------------------------------------------
 587
 588
 589 //--------------------------------------------------------------------
 590 template <class ImageType>
 591 void
 592 clitk::ExtractLungFilter<ImageType>::
 593 SearchForTrachea()
 594 {
 595   // Search for seed among n slices, skip some slices before starting
 596   // if not found -> skip more and restart
 597
 598   // when seed found : perform region growing
 599   // compute trachea volume
 600   // if volume not plausible  -> skip more slices and restart
 601
 602   bool stop = false;
 603   double volume = 0.0;
 604   int skip = GetNumberOfSlicesToSkipBeforeSearchingSeed();
 605   while (!stop) {
 606     stop = SearchForTracheaSeed(skip);
 607     if (stop) {
 608       TracheaRegionGrowing();
 609       volume = ComputeTracheaVolume()/1000; // assume mm3, so divide by 1000 to get cc
 610       if (GetWriteStepFlag()) {
 611         writeImage<MaskImageType>(trachea, "step-trachea-"+toString(skip)+".mhd");
 612       }
 613       if (GetTracheaVolumeMustBeCheckedFlag()) {
 614         if ((volume > 10) && (volume < 65 )) { // depend on image size ...
 615           // Typical volume 22.59 cm 3 (± 7.69 cm 3 ) [ Leader 2004 ]
 616           if (GetVerboseStepFlag()) {
 617             std::cout << "\t Found trachea with volume " << volume << " cc." << std::endl;
 618           }
 619           stop = true;
 620         }
 621         else {
 622           if (GetVerboseStepFlag()) {
 623             std::cout << "\t The volume of the trachea (" << volume
 624                       << " cc) seems not correct. I skip some slices (" << skip << ") and restart to find seeds."
 625                       << std::endl;
 626           }
 627           skip += 5;
 628           stop = false;
 629           // empty the list of seed
 630           m_Seeds.clear();
 631         }
 632       }
 633       else {
 634         stop = true;
 635       }
 636     }
 637   }
 638
 639   if (volume != 0.0) {
 640     // Set output
 641     StopCurrentStep<MaskImageType>(trachea);
 642   }
 643   else { // Trachea not found
 644     this->SetWarning("* WARNING * No seed found for trachea.");
 645     StopCurrentStep();
 646   }
 647 }
 648 //--------------------------------------------------------------------
 649
 650 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX