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