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