segmentation/clitkExtractLungFilter.txx

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