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 #include "itkConstantPadImageFilter.h"
  41 #include <itkBinaryBallStructuringElement.h>
  42 #include "itkStatisticsLabelObject.h"
  43 #include "itkLabelMap.h"
  44 #include "itkLabelImageToShapeLabelMapFilter.h"
  45 #include "itkLabelMapToLabelImageFilter.h"
  46 #include "itkExtractImageFilter.h"
  47 #include "itkOrientImageFilter.h"
  48 #include "itkSpatialOrientationAdapter.h"
  49 #include "itkImageDuplicator.h"
  50 #include <fcntl.h>
  51
  52 //--------------------------------------------------------------------
  53 template <class ImageType>
  54 clitk::ExtractLungFilter<ImageType>::
  55 ExtractLungFilter():
  56   clitk::FilterBase(),
  57   clitk::FilterWithAnatomicalFeatureDatabaseManagement(),
  58   itk::ImageToImageFilter<ImageType, MaskImageType>()
  59 {
  60   SetNumberOfSteps(10);
  61   m_MaxSeedNumber = 500;
  62
  63   // Default global options
  64   this->SetNumberOfRequiredInputs(1);
  65   SetPatientMaskBackgroundValue(0);
  66   SetBackgroundValue(0); // Must be zero
  67   SetForegroundValue(1);
  68   SetMinimalComponentSize(100);
  69   VerboseRegionGrowingFlagOff();
  70
  71   // Step 1 default values
  72   SetUpperThreshold(-300);
  73   SetLowerThreshold(-1000);
  74   UseLowerThresholdOff();
  75   LabelParamType * p1 = new LabelParamType;
  76   p1->SetFirstKeep(1);
  77   p1->UseLastKeepOff();
  78   p1->AddLabelToRemove(2);
  79   SetLabelizeParameters1(p1);
  80
  81   // Step 2 default values
  82   SetTracheaSeedAlgorithm(0);
  83   SetUpperThresholdForTrachea(-700);
  84   SetMultiplierForTrachea(5);
  85   SetThresholdStepSizeForTrachea(64);
  86   SetNumberOfSlicesToSkipBeforeSearchingSeed(0);
  87   TracheaVolumeMustBeCheckedFlagOn();
  88   SetNumSlices(50);
  89   SetMaxElongation(0.5);
  90   SetSeedPreProcessingThreshold(-400);
  91
  92
  93   // Step 3 default values
  94   SetNumberOfHistogramBins(500);
  95   LabelParamType * p2 = new LabelParamType;
  96   p2->SetFirstKeep(1);
  97   p2->UseLastKeepOff();
  98   // p->AddLabelToRemove(2); // No label to remove by default
  99   SetLabelizeParameters2(p2);
 100
 101   // Step 4 default values
 102   SetRadiusForTrachea(1);
 103   LabelParamType * p3 = new LabelParamType;
 104   p3->SetFirstKeep(1);
 105   p3->SetLastKeep(2);
 106   p3->UseLastKeepOff();
 107   SetLabelizeParameters3(p3);
 108
 109   // Step 5
 110   OpenCloseFlagOff();
 111   SetOpenCloseRadius(1);
 112
 113   // Step 6
 114   FillHolesFlagOn();
 115   AutoCropOn();
 116 }
 117 //--------------------------------------------------------------------
 118
 119
 120 //--------------------------------------------------------------------
 121 template <class ImageType>
 122 void
 123 clitk::ExtractLungFilter<ImageType>::
 124 SetInput(const ImageType * image)
 125 {
 126   this->SetNthInput(0, const_cast<ImageType *>(image));
 127 }
 128 //--------------------------------------------------------------------
 129
 130
 131 //--------------------------------------------------------------------
 132 template <class ImageType>
 133 void
 134 clitk::ExtractLungFilter<ImageType>::
 135 AddSeed(InternalIndexType s)
 136 {
 137   m_Seeds.push_back(s);
 138 }
 139 //--------------------------------------------------------------------
 140
 141
 142 //--------------------------------------------------------------------
 143 template <class ImageType>
 144 void
 145 clitk::ExtractLungFilter<ImageType>::
 146 GenerateOutputInformation()
 147 {
 148   clitk::PrintMemory(GetVerboseMemoryFlag(), "Initial memory"); // OK
 149   Superclass::GenerateOutputInformation();
 150
 151   // Read DB
 152   LoadAFDB();
 153
 154   // Get input pointers
 155   input   = dynamic_cast<const ImageType*>(itk::ProcessObject::GetInput(0));
 156   patient = GetAFDB()->template GetImage <MaskImageType>("Patient");
 157   PrintMemory(GetVerboseMemoryFlag(), "After reading patient"); // OK
 158
 159   //--------------------------------------------------------------------
 160   //--------------------------------------------------------------------
 161   // Crop input like patient image (must have the same spacing)
 162   // It copy the input if the same are the same
 163   StartNewStep("Copy and crop input image to 'patient' extends");
 164   typedef clitk::CropLikeImageFilter<ImageType> CropImageFilter;
 165   typename CropImageFilter::Pointer cropFilter = CropImageFilter::New();
 166   // cropFilter->ReleaseDataFlagOn(); // NO=seg fault !!
 167   cropFilter->SetInput(input);
 168   cropFilter->SetCropLikeImage(patient);
 169   cropFilter->Update();
 170   working_input = cropFilter->GetOutput();
 171   //  cropFilter->Delete(); // NO !!!! if yes, sg fault, Cropfilter is buggy !??
 172   StopCurrentStep<ImageType>(working_input);
 173   PrintMemory(GetVerboseMemoryFlag(), "After crop"); // OK, slightly more than a copy of input
 174
 175   //--------------------------------------------------------------------
 176   //--------------------------------------------------------------------
 177   StartNewStep("Set background to initial image");
 178   working_input = SetBackground<ImageType, MaskImageType>
 179     (working_input, patient, GetPatientMaskBackgroundValue(), -1000, true);
 180
 181   // Pad images with air to prevent patient touching the image border
 182   static const unsigned int Dim = ImageType::ImageDimension;
 183   typedef itk::ConstantPadImageFilter<ImageType, ImageType> PadFilterType;
 184   typename PadFilterType::Pointer padFilter = PadFilterType::New();
 185   padFilter->SetInput(working_input);
 186   padFilter->SetConstant(-1000);
 187   typename ImageType::SizeType bounds;
 188   for (unsigned i = 0; i < Dim - 1; ++i)
 189     bounds[i] = 1;
 190   bounds[Dim - 1] = 0;
 191   padFilter->SetPadLowerBound(bounds);
 192   padFilter->SetPadUpperBound(bounds);
 193   padFilter->Update();
 194   working_input = padFilter->GetOutput();
 195
 196   StopCurrentStep<ImageType>(working_input);
 197   PrintMemory(GetVerboseMemoryFlag(), "After set bg"); // OK, additional mem = 0
 198
 199   /*
 200   // We do not need patient mask anymore, release memory
 201   GetAFDB()->template ReleaseImage<MaskImageType>("Patient");
 202   DD(patient->GetReferenceCount());
 203   PrintMemory(GetVerboseMemoryFlag(), "After delete patient"); // OK, additional mem = 0
 204   patient->Delete();
 205   PrintMemory(GetVerboseMemoryFlag(), "After delete patient"); // OK, additional mem = 0
 206   */
 207
 208   //--------------------------------------------------------------------
 209   //--------------------------------------------------------------------
 210   StartNewStep("Remove Air");
 211   // Check threshold
 212   if (m_UseLowerThreshold) {
 213     if (m_LowerThreshold > m_UpperThreshold) {
 214       clitkExceptionMacro("lower threshold cannot be greater than upper threshold.");
 215     }
 216   }
 217   // Threshold to get air
 218   typedef itk::BinaryThresholdImageFilter<ImageType, MaskImageType> InputBinarizeFilterType;
 219   typename InputBinarizeFilterType::Pointer binarizeFilter = InputBinarizeFilterType::New();
 220   binarizeFilter->SetInput(working_input);
 221   if (m_UseLowerThreshold) binarizeFilter->SetLowerThreshold(m_LowerThreshold);
 222   // binarizeFilter->CanRunInPlace() is false
 223   binarizeFilter->SetUpperThreshold(m_UpperThreshold);
 224   binarizeFilter->SetInsideValue(this->GetForegroundValue());
 225   binarizeFilter->SetOutsideValue(this->GetBackgroundValue());
 226   binarizeFilter->Update();
 227   working_mask = binarizeFilter->GetOutput();
 228   PrintMemory(GetVerboseMemoryFlag(), "After Binarizefilter"); // OK, additional mem is one mask image
 229
 230   // Labelize and keep right labels
 231   working_mask =
 232     Labelize<MaskImageType>
 233     (working_mask, GetBackgroundValue(), true, GetMinimalComponentSize());
 234   PrintMemory(GetVerboseMemoryFlag(), "After Labelize"); // BUG ? additional mem around 1 time the input ?
 235
 236   working_mask = RemoveLabels<MaskImageType>
 237     (working_mask, GetBackgroundValue(), GetLabelizeParameters1()->GetLabelsToRemove());
 238   PrintMemory(GetVerboseMemoryFlag(), "After RemoveLabels"); // OK additional mem = 0
 239
 240   working_mask = KeepLabels<MaskImageType>
 241     (working_mask,
 242      GetBackgroundValue(),
 243      GetForegroundValue(),
 244      GetLabelizeParameters1()->GetFirstKeep(),
 245      GetLabelizeParameters1()->GetLastKeep(),
 246      GetLabelizeParameters1()->GetUseLastKeep());
 247   PrintMemory(GetVerboseMemoryFlag(), "After KeepLabels to create the 'air'");
 248
 249   // Set Air to BG
 250   working_input = SetBackground<ImageType, MaskImageType>
 251     (working_input, working_mask, this->GetForegroundValue(), this->GetBackgroundValue(), true);
 252   PrintMemory(GetVerboseMemoryFlag(), "After SetBackground");
 253
 254   // End
 255   StopCurrentStep<ImageType>(working_input);
 256
 257   //--------------------------------------------------------------------
 258   //--------------------------------------------------------------------
 259   StartNewStep("Search for the trachea");
 260   SearchForTrachea();
 261   PrintMemory(GetVerboseMemoryFlag(), "After SearchForTrachea");
 262   if (m_Seeds.empty()) {
 263     clitkExceptionMacro("No seeds for trachea... Aborting.");
 264   }
 265
 266   //--------------------------------------------------------------------
 267   //--------------------------------------------------------------------
 268   StartNewStep("Extract the lung with Otsu filter");
 269   // Automated Otsu thresholding and relabeling
 270   typedef itk::OtsuThresholdImageFilter<ImageType,MaskImageType> OtsuThresholdImageFilterType;
 271   typename OtsuThresholdImageFilterType::Pointer otsuFilter=OtsuThresholdImageFilterType::New();
 272   otsuFilter->SetInput(working_input);
 273   otsuFilter->SetNumberOfHistogramBins(GetNumberOfHistogramBins());
 274   otsuFilter->SetInsideValue(this->GetForegroundValue());
 275   otsuFilter->SetOutsideValue(this->GetBackgroundValue());
 276   otsuFilter->Update();
 277   working_mask =  otsuFilter->GetOutput();
 278
 279   // Set output
 280   StopCurrentStep<MaskImageType>(working_mask);
 281   PrintMemory(GetVerboseMemoryFlag(), "After Otsufilter");
 282
 283   //--------------------------------------------------------------------
 284   //--------------------------------------------------------------------
 285   StartNewStep("Select labels");
 286   // Keep right labels
 287   working_mask = LabelizeAndSelectLabels<MaskImageType>
 288     (working_mask,
 289      GetBackgroundValue(),
 290      GetForegroundValue(),
 291      false,
 292      GetMinimalComponentSize(),
 293      GetLabelizeParameters2());
 294
 295   // Set output
 296   StopCurrentStep<MaskImageType>(working_mask);
 297   PrintMemory(GetVerboseMemoryFlag(), "After LabelizeAndSelectLabels");
 298
 299   //--------------------------------------------------------------------
 300   //--------------------------------------------------------------------
 301   if (m_Seeds.size() != 0) { // if ==0 ->no trachea found
 302     StartNewStep("Remove the trachea");
 303     // Set the trachea
 304     working_mask = SetBackground<MaskImageType, MaskImageType>
 305       (working_mask, trachea, 1, -1, true);
 306     PrintMemory(GetVerboseMemoryFlag(), "After SetBackground");
 307
 308     // Dilate the trachea
 309     static const unsigned int Dim = ImageType::ImageDimension;
 310     typedef itk::BinaryBallStructuringElement<InternalPixelType, Dim> KernelType;
 311     KernelType structuringElement;
 312     structuringElement.SetRadius(GetRadiusForTrachea());
 313     structuringElement.CreateStructuringElement();
 314     typedef clitk::ConditionalBinaryDilateImageFilter<MaskImageType, MaskImageType, KernelType> ConditionalBinaryDilateImageFilterType;
 315     typename ConditionalBinaryDilateImageFilterType::Pointer dilateFilter = ConditionalBinaryDilateImageFilterType::New();
 316     dilateFilter->SetBoundaryToForeground(false);
 317     dilateFilter->SetKernel(structuringElement);
 318     dilateFilter->SetBackgroundValue (1);
 319     dilateFilter->SetForegroundValue (-1);
 320     dilateFilter->SetInput (working_mask);
 321     dilateFilter->Update();
 322     working_mask = dilateFilter->GetOutput();
 323     PrintMemory(GetVerboseMemoryFlag(), "After dilate");
 324
 325     // Set trachea with dilatation
 326     trachea = SetBackground<MaskImageType, MaskImageType>
 327       (trachea, working_mask, -1, this->GetForegroundValue(), true);
 328
 329     // Remove the trachea
 330     working_mask = SetBackground<MaskImageType, MaskImageType>
 331       (working_mask, working_mask, -1, this->GetBackgroundValue(), true);
 332
 333     // Label
 334     working_mask = LabelizeAndSelectLabels<MaskImageType>
 335       (working_mask,
 336        GetBackgroundValue(),
 337        GetForegroundValue(),
 338        false,
 339        GetMinimalComponentSize(),
 340        GetLabelizeParameters3());
 341
 342     // Set output
 343     StopCurrentStep<MaskImageType>(working_mask);
 344   }
 345
 346   //--------------------------------------------------------------------
 347   //--------------------------------------------------------------------
 348   PrintMemory(GetVerboseMemoryFlag(), "before autocropfilter");
 349   if (m_Seeds.size() != 0) { // if ==0 ->no trachea found
 350     if (GetAutoCrop())
 351       trachea = clitk::AutoCrop<MaskImageType>(trachea, GetBackgroundValue());
 352     else
 353     {
 354       // Remove Padding region
 355       typedef itk::CropImageFilter<MaskImageType, MaskImageType> CropFilterType;
 356       typename CropFilterType::Pointer cropFilter = CropFilterType::New();
 357       cropFilter->SetInput(trachea);
 358       cropFilter->SetLowerBoundaryCropSize(bounds);
 359       cropFilter->SetUpperBoundaryCropSize(bounds);
 360       cropFilter->Update();
 361       trachea = cropFilter->GetOutput();
 362     }
 363     StopCurrentStep<MaskImageType>(trachea);
 364     PrintMemory(GetVerboseMemoryFlag(), "after delete trachea");
 365   }
 366   PrintMemory(GetVerboseMemoryFlag(), "after delete trachea");
 367
 368   //--------------------------------------------------------------------
 369   //--------------------------------------------------------------------
 370   StartNewStep("Cropping lung");
 371   PrintMemory(GetVerboseMemoryFlag(), "Before Autocropfilter");
 372   if (GetAutoCrop())
 373     working_mask = clitk::AutoCrop<MaskImageType>(working_mask, GetBackgroundValue());
 374   else
 375   {
 376     // Remove Padding region
 377     typedef itk::CropImageFilter<MaskImageType, MaskImageType> CropFilterType;
 378     typename CropFilterType::Pointer cropFilter = CropFilterType::New();
 379     cropFilter->SetInput(working_mask);
 380     cropFilter->SetLowerBoundaryCropSize(bounds);
 381     cropFilter->SetUpperBoundaryCropSize(bounds);
 382     cropFilter->Update();
 383     working_mask = cropFilter->GetOutput();
 384   }
 385   StopCurrentStep<MaskImageType>(working_mask);
 386
 387   //--------------------------------------------------------------------
 388   //--------------------------------------------------------------------
 389   // Final OpenClose
 390   if (GetOpenCloseFlag()) {
 391     StartNewStep("Open/Close");
 392     PrintMemory(GetVerboseMemoryFlag(), "Before OpenClose");
 393
 394     // Structuring element
 395     typedef itk::BinaryBallStructuringElement<InternalPixelType, ImageDimension> KernelType;
 396     KernelType structuringElement;
 397     structuringElement.SetRadius(GetOpenCloseRadius());
 398     structuringElement.CreateStructuringElement();
 399
 400     // Open
 401     typedef itk::BinaryMorphologicalOpeningImageFilter<MaskImageType, InternalImageType, KernelType> OpenFilterType;
 402     typename OpenFilterType::Pointer openFilter = OpenFilterType::New();
 403     openFilter->SetInput(working_mask);
 404     openFilter->SetBackgroundValue(GetBackgroundValue());
 405     openFilter->SetForegroundValue(GetForegroundValue());
 406     openFilter->SetKernel(structuringElement);
 407
 408     // Close
 409     typedef itk::BinaryMorphologicalClosingImageFilter<MaskImageType, MaskImageType, KernelType> CloseFilterType;
 410     typename CloseFilterType::Pointer closeFilter = CloseFilterType::New();
 411     closeFilter->SetInput(openFilter->GetOutput());
 412     closeFilter->SetSafeBorder(true);
 413     closeFilter->SetForegroundValue(GetForegroundValue());
 414     closeFilter->SetKernel(structuringElement);
 415     closeFilter->Update();
 416     working_mask = closeFilter->GetOutput();
 417   }
 418
 419   //--------------------------------------------------------------------
 420   //--------------------------------------------------------------------
 421   // Fill Lungs
 422   if (GetFillHolesFlag()) {
 423     StartNewStep("Fill Holes");
 424     PrintMemory(GetVerboseMemoryFlag(), "Before Fill Holes");
 425     typedef clitk::FillMaskFilter<MaskImageType> FillMaskFilterType;
 426     typename FillMaskFilterType::Pointer fillMaskFilter = FillMaskFilterType::New();
 427     fillMaskFilter->SetInput(working_mask);
 428     fillMaskFilter->AddDirection(2);
 429     //fillMaskFilter->AddDirection(1);
 430     fillMaskFilter->Update();
 431     working_mask = fillMaskFilter->GetOutput();
 432     StopCurrentStep<MaskImageType>(working_mask);
 433   }
 434
 435   if (GetSeparateLungsFlag()) {
 436     //--------------------------------------------------------------------
 437     //--------------------------------------------------------------------
 438     StartNewStep("Separate Left/Right lungs");
 439       PrintMemory(GetVerboseMemoryFlag(), "Before Separate");
 440     // Initial label
 441     working_mask = Labelize<MaskImageType>(working_mask,
 442                                           GetBackgroundValue(),
 443                                           false,
 444                                           GetMinimalComponentSize());
 445
 446     PrintMemory(GetVerboseMemoryFlag(), "After Labelize");
 447
 448     // Count the labels
 449     typedef itk::StatisticsImageFilter<MaskImageType> StatisticsImageFilterType;
 450     typename StatisticsImageFilterType::Pointer statisticsImageFilter=StatisticsImageFilterType::New();
 451     statisticsImageFilter->SetInput(working_mask);
 452     statisticsImageFilter->Update();
 453     unsigned int initialNumberOfLabels = statisticsImageFilter->GetMaximum();
 454     working_mask = statisticsImageFilter->GetOutput();
 455
 456     PrintMemory(GetVerboseMemoryFlag(), "After count label");
 457
 458     // Decompose the first label
 459     if (initialNumberOfLabels<2) {
 460       // Structuring element radius
 461       typename ImageType::SizeType radius;
 462       for (unsigned int i=0;i<Dim;i++) radius[i]=1;
 463       typedef clitk::DecomposeAndReconstructImageFilter<MaskImageType,MaskImageType> DecomposeAndReconstructFilterType;
 464       typename DecomposeAndReconstructFilterType::Pointer decomposeAndReconstructFilter=DecomposeAndReconstructFilterType::New();
 465       decomposeAndReconstructFilter->SetInput(working_mask);
 466       decomposeAndReconstructFilter->SetVerbose(false);
 467       decomposeAndReconstructFilter->SetRadius(radius);
 468       decomposeAndReconstructFilter->SetMaximumNumberOfLabels(2);
 469       decomposeAndReconstructFilter->SetMinimumObjectSize(this->GetMinimalComponentSize());
 470       decomposeAndReconstructFilter->SetMinimumNumberOfIterations(1);
 471       decomposeAndReconstructFilter->SetBackgroundValue(this->GetBackgroundValue());
 472       decomposeAndReconstructFilter->SetForegroundValue(this->GetForegroundValue());
 473       decomposeAndReconstructFilter->SetFullyConnected(true);
 474       decomposeAndReconstructFilter->SetNumberOfNewLabels(1);
 475       decomposeAndReconstructFilter->Update();
 476       working_mask = decomposeAndReconstructFilter->GetOutput();
 477     }
 478     PrintMemory(GetVerboseMemoryFlag(), "After decomposeAndReconstructFilter");
 479   }
 480
 481   // Retain labels ('1' is largset lung, so right. '2' is left)
 482   typedef itk::ThresholdImageFilter<MaskImageType> ThresholdImageFilterType;
 483   typename ThresholdImageFilterType::Pointer thresholdFilter = ThresholdImageFilterType::New();
 484   thresholdFilter->SetInput(working_mask);
 485   thresholdFilter->ThresholdAbove(2);
 486   thresholdFilter->SetOutsideValue(this->GetBackgroundValue());
 487   thresholdFilter->Update();
 488   working_mask = thresholdFilter->GetOutput();
 489   StopCurrentStep<MaskImageType> (working_mask);
 490   PrintMemory(GetVerboseMemoryFlag(), "After Thresholdfilter");
 491
 492   // Update output info
 493   //  output = working_mask;
 494   //this->GetOutput(0)->SetRegions(output->GetLargestPossibleRegion());
 495
 496   //  this->GetOutput(0)->SetRegions(working_mask->GetLargestPossibleRegion());
 497
 498 }
 499 //--------------------------------------------------------------------
 500
 501
 502 //--------------------------------------------------------------------
 503 template <class TImageType>
 504 void
 505 clitk::ExtractLungFilter<TImageType>::
 506 GenerateInputRequestedRegion() {
 507   //  DD("GenerateInputRequestedRegion (nothing?)");
 508 }
 509 //--------------------------------------------------------------------
 510
 511
 512 //--------------------------------------------------------------------
 513 template <class ImageType>
 514 void
 515 clitk::ExtractLungFilter<ImageType>::
 516 GenerateData()
 517 {
 518   // Set the output
 519   //  this->GraftOutput(output); // not SetNthOutput
 520   this->GraftOutput(working_mask); // not SetNthOutput
 521   // Store image filenames into AFDB
 522   GetAFDB()->SetImageFilename("Lungs", this->GetOutputLungFilename());
 523   GetAFDB()->SetImageFilename("Trachea", this->GetOutputTracheaFilename());
 524   WriteAFDB();
 525 }
 526 //--------------------------------------------------------------------
 527
 528
 529 //--------------------------------------------------------------------
 530 template <class ImageType>
 531 bool
 532 clitk::ExtractLungFilter<ImageType>::
 533 SearchForTracheaSeed(int skip)
 534 {
 535   if (m_Seeds.size() == 0) { // try to find seed (if not zero, it is given by user)
 536     // Restart the search until a seed is found, skipping more and more slices
 537     bool stop = false;
 538     while (!stop) {
 539       // Search seed (parameters = UpperThresholdForTrachea)
 540       static const unsigned int Dim = ImageType::ImageDimension;
 541       typename InternalImageType::RegionType sliceRegion = working_input->GetLargestPossibleRegion();
 542       typename InternalImageType::SizeType sliceRegionSize = sliceRegion.GetSize();
 543       typename InternalImageType::IndexType sliceRegionIndex = sliceRegion.GetIndex();
 544       sliceRegionIndex[Dim-1]=sliceRegionSize[Dim-1]-skip-5;
 545       sliceRegionSize[Dim-1]=5;
 546       sliceRegion.SetSize(sliceRegionSize);
 547       sliceRegion.SetIndex(sliceRegionIndex);
 548       typedef  itk::ImageRegionConstIterator<ImageType> IteratorType;
 549       IteratorType it(working_input, sliceRegion);
 550       it.GoToBegin();
 551       while (!it.IsAtEnd()) {
 552         if(it.Get() < GetUpperThresholdForTrachea() ) {
 553           AddSeed(it.GetIndex());
 554           // DD(it.GetIndex());
 555         }
 556         ++it;
 557       }
 558
 559       // if we do not found : restart
 560       stop = (m_Seeds.size() != 0);
 561       if (!stop) {
 562         if (GetVerboseStepFlag()) {
 563           std::cout << "\t No seed found this time. I skip some slices and restart." << std::endl;
 564         }
 565         if (skip > 0.5 * working_input->GetLargestPossibleRegion().GetSize()[2]) {
 566           // we want to skip more than a half of the image, it is probably a bug
 567           std::cerr << "Number of slices to skip to find trachea to high = " << skip << std::endl;
 568           stop = true;
 569         }
 570         skip += 5;
 571       }
 572       else {
 573         // DD(m_Seeds[0]);
 574         // DD(m_Seeds.size());
 575       }
 576     }
 577   }
 578   return (m_Seeds.size() != 0);
 579 }
 580 //--------------------------------------------------------------------
 581
 582
 583 bool is_orientation_superior(itk::SpatialOrientation::ValidCoordinateOrientationFlags orientation)
 584 {
 585   itk::SpatialOrientation::CoordinateTerms sup = itk::SpatialOrientation::ITK_COORDINATE_Superior;
 586   bool primary = (orientation & 0x0000ff) == sup;
 587   bool secondary = ((orientation & 0x00ff00) >> 8) == sup;
 588   bool tertiary = ((orientation & 0xff0000) >> 16) == sup;
 589   return primary || secondary || tertiary;
 590 }
 591
 592 //--------------------------------------------------------------------
 593 template <class ImageType>
 594 bool
 595 clitk::ExtractLungFilter<ImageType>::
 596 SearchForTracheaSeed2(int numberOfSlices)
 597 {
 598   if (m_Seeds.size() == 0) { // try to find seed (if not zero, it is given by user)
 599     if (GetVerboseRegionGrowingFlag())
 600       std::cout << "SearchForTracheaSeed2(" << numberOfSlices << ", " << GetMaxElongation() << ")" << std::endl;
 601
 602     typedef unsigned char MaskPixelType;
 603     typedef itk::Image<MaskPixelType, ImageType::ImageDimension> MaskImageType;
 604     typedef itk::Image<typename MaskImageType::PixelType, 2> MaskImageType2D;
 605     typedef itk::BinaryThresholdImageFilter<ImageType, MaskImageType> ThresholdFilterType;
 606     typedef itk::BinaryBallStructuringElement<MaskPixelType, MaskImageType::ImageDimension> KernelType;
 607     typedef itk::BinaryMorphologicalClosingImageFilter<MaskImageType, MaskImageType, KernelType> ClosingFilterType;
 608     typedef itk::BinaryMorphologicalOpeningImageFilter<MaskImageType, MaskImageType, KernelType> OpeningFilterType;
 609     typedef itk::ExtractImageFilter<MaskImageType, MaskImageType2D> SlicerFilterType;
 610     typedef itk::ConnectedComponentImageFilter<MaskImageType2D, MaskImageType2D> LabelFilterType;
 611     typedef itk::ShapeLabelObject<MaskPixelType, MaskImageType2D::ImageDimension> ShapeLabelType;
 612     typedef itk::LabelMap<ShapeLabelType> LabelImageType;
 613     typedef itk::LabelImageToShapeLabelMapFilter<MaskImageType2D, LabelImageType> ImageToLabelMapFilterType;
 614     typedef itk::LabelMapToLabelImageFilter<LabelImageType, MaskImageType2D> LabelMapToImageFilterType;
 615     typedef itk::ImageFileWriter<MaskImageType2D> FileWriterType;
 616
 617     // threshold to isolate airawys and lungs
 618     typename ThresholdFilterType::Pointer threshold = ThresholdFilterType::New();
 619     threshold->SetLowerThreshold(-2000);
 620     threshold->SetUpperThreshold(GetSeedPreProcessingThreshold());
 621     threshold->SetInput(working_input);
 622     threshold->Update();
 623
 624     KernelType kernel_closing, kernel_opening;
 625
 626     // remove small noise
 627     typename OpeningFilterType::Pointer opening = OpeningFilterType::New();
 628     kernel_opening.SetRadius(1);
 629     opening->SetKernel(kernel_opening);
 630     opening->SetInput(threshold->GetOutput());
 631     opening->Update();
 632
 633     typename SlicerFilterType::Pointer slicer = SlicerFilterType::New();
 634     slicer->SetInput(opening->GetOutput());
 635
 636     // label result
 637     typename LabelFilterType::Pointer label_filter = LabelFilterType::New();
 638     label_filter->SetInput(slicer->GetOutput());
 639
 640     // extract shape information from labels
 641     typename ImageToLabelMapFilterType::Pointer label_to_map_filter = ImageToLabelMapFilterType::New();
 642     label_to_map_filter->SetInput(label_filter->GetOutput());
 643
 644     typename LabelMapToImageFilterType::Pointer map_to_label_filter = LabelMapToImageFilterType::New();
 645     typename FileWriterType::Pointer writer = FileWriterType::New();
 646
 647     typename ImageType::IndexType index;
 648     typename ImageType::RegionType region = working_input->GetLargestPossibleRegion();
 649     typename ImageType::SizeType size = region.GetSize();
 650     typename ImageType::SpacingType spacing = working_input->GetSpacing();
 651     typename ImageType::PointType origin = working_input->GetOrigin();
 652
 653     int nslices = min(numberOfSlices, size[2]);
 654     int start = 0, increment = 1;
 655     itk::SpatialOrientationAdapter orientation;
 656     typename ImageType::DirectionType dir = working_input->GetDirection();
 657     if (!is_orientation_superior(orientation.FromDirectionCosines(dir))) {
 658       start = size[2]-1;
 659       increment = -1;
 660     }
 661
 662     typename MaskImageType::PointType image_centre;
 663     image_centre[0] = size[0]/2;
 664     image_centre[1] = size[1]/2;
 665     image_centre[2] = 0;
 666
 667     typedef InternalIndexType SeedType;
 668     SeedType trachea_centre, shape_centre, max_e_centre, prev_e_centre;
 669     typedef std::list<SeedType> PointListType;
 670     typedef std::list<PointListType> SequenceListType;
 671     PointListType* current_sequence = NULL;
 672     SequenceListType sequence_list;
 673
 674     prev_e_centre.Fill(0);
 675     std::ostringstream file_name;
 676     index[0] = index[1] = 0;
 677     size[0] = size[1] = 512;
 678     size[2] = 0;
 679     while (nslices--) {
 680       index[2] = start;
 681       start += increment;
 682
 683       region.SetIndex(index);
 684       region.SetSize(size);
 685       slicer->SetExtractionRegion(region);
 686       slicer->Update();
 687       label_filter->SetInput(slicer->GetOutput());
 688       label_filter->Update();
 689
 690       label_to_map_filter->SetInput(label_filter->GetOutput());
 691       label_to_map_filter->Update();
 692       typename LabelImageType::Pointer label_map = label_to_map_filter->GetOutput();
 693
 694       if (GetWriteStepFlag()) {
 695         map_to_label_filter->SetInput(label_map);
 696         writer->SetInput(map_to_label_filter->GetOutput());
 697         file_name.str("");
 698         file_name << "labels_";
 699         file_name.width(3);
 700         file_name.fill('0');
 701         file_name << index[2] << ".mhd";
 702         writer->SetFileName(file_name.str().c_str());
 703         writer->Update();
 704       }
 705
 706       typename ShapeLabelType::Pointer shape, max_e_shape;
 707       double max_elongation = GetMaxElongation();
 708       double max_size = size[0]*size[1]/128;
 709       double max_e = 0;
 710       int nshapes = 0;
 711       unsigned int nlables = label_map->GetNumberOfLabelObjects();
 712       for (unsigned int j = 0; j < nlables; j++) {
 713         shape = label_map->GetNthLabelObject(j);
 714         if (shape->Size() > 150 && shape->Size() <= max_size) {
 715 #if ITK_VERSION_MAJOR < 4
 716           double e = 1 - 1/shape->GetBinaryElongation();
 717 #else
 718           double e = 1 - 1/shape->GetElongation();
 719 #endif
 720           //double area = 1 - r->Area() ;
 721           if (e < max_elongation) {
 722             nshapes++;
 723             shape_centre[0] = (shape->GetCentroid()[0] - origin[0])/spacing[0];
 724             shape_centre[1] = (shape->GetCentroid()[1] - origin[1])/spacing[1];
 725             shape_centre[2] = index[2];
 726             //double d = 1 - (shape_centre - image_centre).Magnitude()/max_dist;
 727             double dx = shape_centre[0] - image_centre[0];
 728             double d = 1 - dx*2/size[0];
 729             e = e + d;
 730             if (e > max_e)
 731             {
 732               max_e = e;
 733               max_e_shape = shape;
 734               max_e_centre = shape_centre;
 735             }
 736           }
 737         }
 738       }
 739
 740       if (nshapes > 0)
 741       {
 742         itk::Point<typename SeedType::IndexValueType, ImageType::ImageDimension> p1, p2;
 743         p1[0] = max_e_centre[0];
 744         p1[1] = max_e_centre[1];
 745         p1[2] = max_e_centre[2];
 746
 747         p2[0] = prev_e_centre[0];
 748         p2[1] = prev_e_centre[1];
 749         p2[2] = prev_e_centre[2];
 750
 751         double mag = (p2 - p1).GetNorm();
 752         if (GetVerboseRegionGrowingFlag()) {
 753           cout.precision(3);
 754           cout << index[2] << ": ";
 755           cout << "region(" << max_e_centre[0] << ", " << max_e_centre[1] << ", " << max_e_centre[2] << "); ";
 756           cout << "prev_region(" << prev_e_centre[0] << ", " << prev_e_centre[1] << ", " << prev_e_centre[2] << "); ";
 757           cout << "mag(" << mag << "); " << endl;
 758         }
 759
 760         if (mag > 5)
 761         {
 762           PointListType point_list;
 763           point_list.push_back(max_e_centre);
 764           sequence_list.push_back(point_list);
 765           current_sequence = &sequence_list.back();
 766         }
 767         else if (current_sequence)
 768           current_sequence->push_back(max_e_centre);
 769
 770         prev_e_centre= max_e_centre;
 771       }
 772       else {
 773         if (GetVerboseRegionGrowingFlag()) {
 774           cout << "No shapes found at slice " << index[2] << std::endl;
 775         }
 776       }
 777     }
 778
 779     size_t longest = 0;
 780     for (typename SequenceListType::iterator s = sequence_list.begin(); s != sequence_list.end(); s++)
 781     {
 782       if (s->size() > longest)
 783       {
 784         longest = s->size();
 785         trachea_centre = s->front();
 786       }
 787     }
 788
 789     if (longest > 0) {
 790       if (GetVerboseRegionGrowingFlag())
 791         std::cout << "seed at: " << trachea_centre << std::endl;
 792       m_Seeds.push_back(trachea_centre);
 793     }
 794   }
 795
 796   return (m_Seeds.size() != 0);
 797 }
 798 //--------------------------------------------------------------------
 799
 800 //--------------------------------------------------------------------
 801 template <class ImageType>
 802 void
 803 clitk::ExtractLungFilter<ImageType>::
 804 TracheaRegionGrowing()
 805 {
 806   // Explosion controlled region growing
 807   PrintMemory(GetVerboseMemoryFlag(), "Before ExplosionControlledThresholdConnectedImageFilter");
 808   typedef clitk::ExplosionControlledThresholdConnectedImageFilter<ImageType, MaskImageType> ImageFilterType;
 809   typename ImageFilterType::Pointer f= ImageFilterType::New();
 810   f->SetInput(working_input);
 811   f->SetLower(-2000);
 812   f->SetUpper(GetUpperThresholdForTrachea());
 813   f->SetMinimumLowerThreshold(-2000);
 814   //  f->SetMaximumUpperThreshold(0); // MAYBE TO CHANGE ???
 815   f->SetMaximumUpperThreshold(-700); // MAYBE TO CHANGE ???
 816   f->SetAdaptLowerBorder(false);
 817   f->SetAdaptUpperBorder(true);
 818   f->SetMinimumSize(5000);
 819   f->SetReplaceValue(1);
 820   f->SetMultiplier(GetMultiplierForTrachea());
 821   f->SetThresholdStepSize(GetThresholdStepSizeForTrachea());
 822   f->SetMinimumThresholdStepSize(1);
 823   f->SetVerbose(GetVerboseRegionGrowingFlag());
 824   for(unsigned int i=0; i<m_Seeds.size();i++) {
 825     f->AddSeed(m_Seeds[i]);
 826     // DD(m_Seeds[i]);
 827   }
 828   f->Update();
 829   PrintMemory(GetVerboseMemoryFlag(), "After RG update");
 830
 831   // take first (main) connected component
 832   trachea = Labelize<MaskImageType>(f->GetOutput(),
 833                                     GetBackgroundValue(),
 834                                     true,
 835                                     1000);//GetMinimalComponentSize());
 836   PrintMemory(GetVerboseMemoryFlag(), "After Labelize");
 837   trachea = KeepLabels<MaskImageType>(trachea,
 838                                               GetBackgroundValue(),
 839                                               GetForegroundValue(),
 840                                               1, 1, false);
 841   PrintMemory(GetVerboseMemoryFlag(), "After KeepLabels");
 842 }
 843 //--------------------------------------------------------------------
 844
 845
 846 //--------------------------------------------------------------------
 847 template <class ImageType>
 848 double
 849 clitk::ExtractLungFilter<ImageType>::
 850 ComputeTracheaVolume()
 851 {
 852   typedef itk::ImageRegionConstIterator<InternalImageType> IteratorType;
 853   IteratorType iter(trachea, trachea->GetLargestPossibleRegion());
 854   iter.GoToBegin();
 855   double volume = 0.0;
 856   while (!iter.IsAtEnd()) {
 857     if (iter.Get() == this->GetForegroundValue()) volume++;
 858     ++iter;
 859   }
 860
 861   double voxelsize = trachea->GetSpacing()[0]*trachea->GetSpacing()[1]*trachea->GetSpacing()[2];
 862   return volume*voxelsize;
 863 }
 864 //--------------------------------------------------------------------
 865
 866
 867 //--------------------------------------------------------------------
 868 template <class ImageType>
 869 void
 870 clitk::ExtractLungFilter<ImageType>::
 871 SearchForTrachea()
 872 {
 873   // Search for seed among n slices, skip some slices before starting
 874   // if not found -> skip more and restart
 875
 876   // when seed found : perform region growing
 877   // compute trachea volume
 878   // if volume not plausible  -> skip more slices and restart
 879
 880   bool has_seed;
 881   bool stop = false;
 882   double volume = 0.0;
 883   int skip = GetNumberOfSlicesToSkipBeforeSearchingSeed();
 884   while (!stop) {
 885     stop = true;
 886     if (GetTracheaSeedAlgorithm() == 0)
 887       has_seed = SearchForTracheaSeed(skip);
 888     else
 889       has_seed = SearchForTracheaSeed2(GetNumSlices());
 890
 891     if (has_seed) {
 892       TracheaRegionGrowing();
 893       volume = ComputeTracheaVolume()/1000; // assume mm3, so divide by 1000 to get cc
 894       if (GetWriteStepFlag()) {
 895         writeImage<MaskImageType>(trachea, "step-trachea-"+toString(skip)+".mhd");
 896       }
 897       if (GetTracheaVolumeMustBeCheckedFlag()) {
 898         if ((volume > 10) && (volume < 65 )) { // depend on image size ...
 899           // Typical volume 22.59 cm 3 (± 7.69 cm 3 ) [ Leader 2004 ]
 900           if (GetVerboseStepFlag())
 901           {
 902             std::cout << "\t Found trachea with volume " << volume << " cc." << std::endl;
 903           }
 904         }
 905         else
 906           if (GetTracheaSeedAlgorithm() == 0) {
 907             if (GetVerboseStepFlag()) {
 908               std::cout << "\t The volume of the trachea (" << volume
 909                         << " cc) seems not correct. I skip some slices (" << skip << ") and restart to find seeds."
 910                         << std::endl;
 911             }
 912             skip += 5;
 913             stop = false;
 914             // empty the list of seed
 915             m_Seeds.clear();
 916         }
 917         if (skip > 0.5 * working_input->GetLargestPossibleRegion().GetSize()[2]) {
 918           // we want to skip more than a half of the image, it is probably a bug
 919           std::cerr << "2 : Number of slices to skip to find trachea too high = " << skip << std::endl;
 920           stop = true;
 921         }
 922       }
 923       else {
 924         stop = true;
 925       }
 926     }
 927   }
 928
 929   if (volume != 0.0) {
 930     // Set output
 931     StopCurrentStep<MaskImageType>(trachea);
 932   }
 933   else { // Trachea not found
 934     this->SetWarning("* WARNING * No seed found for trachea.");
 935     StopCurrentStep();
 936   }
 937 }
 938 //--------------------------------------------------------------------
 939
 940 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX