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