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
 263   //--------------------------------------------------------------------
 264   //--------------------------------------------------------------------
 265   StartNewStep("Extract the lung with Otsu filter");
 266   // Automated Otsu thresholding and relabeling
 267   typedef itk::OtsuThresholdImageFilter<ImageType,MaskImageType> OtsuThresholdImageFilterType;
 268   typename OtsuThresholdImageFilterType::Pointer otsuFilter=OtsuThresholdImageFilterType::New();
 269   otsuFilter->SetInput(working_input);
 270   otsuFilter->SetNumberOfHistogramBins(GetNumberOfHistogramBins());
 271   otsuFilter->SetInsideValue(this->GetForegroundValue());
 272   otsuFilter->SetOutsideValue(this->GetBackgroundValue());
 273   otsuFilter->Update();
 274   working_mask =  otsuFilter->GetOutput();
 275
 276   // Set output
 277   StopCurrentStep<MaskImageType>(working_mask);
 278   PrintMemory(GetVerboseMemoryFlag(), "After Otsufilter");
 279
 280   //--------------------------------------------------------------------
 281   //--------------------------------------------------------------------
 282   StartNewStep("Select labels");
 283   // Keep right labels
 284   working_mask = LabelizeAndSelectLabels<MaskImageType>
 285     (working_mask,
 286      GetBackgroundValue(),
 287      GetForegroundValue(),
 288      false,
 289      GetMinimalComponentSize(),
 290      GetLabelizeParameters2());
 291
 292   // Set output
 293   StopCurrentStep<MaskImageType>(working_mask);
 294   PrintMemory(GetVerboseMemoryFlag(), "After LabelizeAndSelectLabels");
 295
 296   //--------------------------------------------------------------------
 297   //--------------------------------------------------------------------
 298   if (m_Seeds.size() != 0) { // if ==0 ->no trachea found
 299     StartNewStep("Remove the trachea");
 300     // Set the trachea
 301     working_mask = SetBackground<MaskImageType, MaskImageType>
 302       (working_mask, trachea, 1, -1, true);
 303     PrintMemory(GetVerboseMemoryFlag(), "After SetBackground");
 304
 305     // Dilate the trachea
 306     static const unsigned int Dim = ImageType::ImageDimension;
 307     typedef itk::BinaryBallStructuringElement<InternalPixelType, Dim> KernelType;
 308     KernelType structuringElement;
 309     structuringElement.SetRadius(GetRadiusForTrachea());
 310     structuringElement.CreateStructuringElement();
 311     typedef clitk::ConditionalBinaryDilateImageFilter<MaskImageType, MaskImageType, KernelType> ConditionalBinaryDilateImageFilterType;
 312     typename ConditionalBinaryDilateImageFilterType::Pointer dilateFilter = ConditionalBinaryDilateImageFilterType::New();
 313     dilateFilter->SetBoundaryToForeground(false);
 314     dilateFilter->SetKernel(structuringElement);
 315     dilateFilter->SetBackgroundValue (1);
 316     dilateFilter->SetForegroundValue (-1);
 317     dilateFilter->SetInput (working_mask);
 318     dilateFilter->Update();
 319     working_mask = dilateFilter->GetOutput();
 320     PrintMemory(GetVerboseMemoryFlag(), "After dilate");
 321
 322     // Set trachea with dilatation
 323     trachea = SetBackground<MaskImageType, MaskImageType>
 324       (trachea, working_mask, -1, this->GetForegroundValue(), true);
 325
 326     // Remove the trachea
 327     working_mask = SetBackground<MaskImageType, MaskImageType>
 328       (working_mask, working_mask, -1, this->GetBackgroundValue(), true);
 329
 330     // Label
 331     working_mask = LabelizeAndSelectLabels<MaskImageType>
 332       (working_mask,
 333        GetBackgroundValue(),
 334        GetForegroundValue(),
 335        false,
 336        GetMinimalComponentSize(),
 337        GetLabelizeParameters3());
 338
 339     // Set output
 340     StopCurrentStep<MaskImageType>(working_mask);
 341   }
 342
 343   //--------------------------------------------------------------------
 344   //--------------------------------------------------------------------
 345   PrintMemory(GetVerboseMemoryFlag(), "before autocropfilter");
 346   if (m_Seeds.size() != 0) { // if ==0 ->no trachea found
 347     if (GetAutoCrop())
 348       trachea = clitk::AutoCrop<MaskImageType>(trachea, GetBackgroundValue());
 349     else
 350     {
 351       // Remove Padding region
 352       typedef itk::CropImageFilter<MaskImageType, MaskImageType> CropFilterType;
 353       typename CropFilterType::Pointer cropFilter = CropFilterType::New();
 354       cropFilter->SetInput(trachea);
 355       cropFilter->SetLowerBoundaryCropSize(bounds);
 356       cropFilter->SetUpperBoundaryCropSize(bounds);
 357       cropFilter->Update();
 358       trachea = cropFilter->GetOutput();
 359     }
 360     StopCurrentStep<MaskImageType>(trachea);
 361     PrintMemory(GetVerboseMemoryFlag(), "after delete trachea");
 362   }
 363   PrintMemory(GetVerboseMemoryFlag(), "after delete trachea");
 364
 365   //--------------------------------------------------------------------
 366   //--------------------------------------------------------------------
 367   StartNewStep("Cropping lung");
 368   PrintMemory(GetVerboseMemoryFlag(), "Before Autocropfilter");
 369   if (GetAutoCrop())
 370     working_mask = clitk::AutoCrop<MaskImageType>(working_mask, GetBackgroundValue());
 371   else
 372   {
 373     // Remove Padding region
 374     typedef itk::CropImageFilter<MaskImageType, MaskImageType> CropFilterType;
 375     typename CropFilterType::Pointer cropFilter = CropFilterType::New();
 376     cropFilter->SetInput(working_mask);
 377     cropFilter->SetLowerBoundaryCropSize(bounds);
 378     cropFilter->SetUpperBoundaryCropSize(bounds);
 379     cropFilter->Update();
 380     working_mask = cropFilter->GetOutput();
 381   }
 382   StopCurrentStep<MaskImageType>(working_mask);
 383
 384   //--------------------------------------------------------------------
 385   //--------------------------------------------------------------------
 386   // Final OpenClose
 387   if (GetOpenCloseFlag()) {
 388     StartNewStep("Open/Close");
 389     PrintMemory(GetVerboseMemoryFlag(), "Before OpenClose");
 390
 391     // Structuring element
 392     typedef itk::BinaryBallStructuringElement<InternalPixelType, ImageDimension> KernelType;
 393     KernelType structuringElement;
 394     structuringElement.SetRadius(GetOpenCloseRadius());
 395     structuringElement.CreateStructuringElement();
 396
 397     // Open
 398     typedef itk::BinaryMorphologicalOpeningImageFilter<MaskImageType, InternalImageType, KernelType> OpenFilterType;
 399     typename OpenFilterType::Pointer openFilter = OpenFilterType::New();
 400     openFilter->SetInput(working_mask);
 401     openFilter->SetBackgroundValue(GetBackgroundValue());
 402     openFilter->SetForegroundValue(GetForegroundValue());
 403     openFilter->SetKernel(structuringElement);
 404
 405     // Close
 406     typedef itk::BinaryMorphologicalClosingImageFilter<MaskImageType, MaskImageType, KernelType> CloseFilterType;
 407     typename CloseFilterType::Pointer closeFilter = CloseFilterType::New();
 408     closeFilter->SetInput(openFilter->GetOutput());
 409     closeFilter->SetSafeBorder(true);
 410     closeFilter->SetForegroundValue(GetForegroundValue());
 411     closeFilter->SetKernel(structuringElement);
 412     closeFilter->Update();
 413     working_mask = closeFilter->GetOutput();
 414   }
 415
 416   //--------------------------------------------------------------------
 417   //--------------------------------------------------------------------
 418   // Fill Lungs
 419   if (GetFillHolesFlag()) {
 420     StartNewStep("Fill Holes");
 421     PrintMemory(GetVerboseMemoryFlag(), "Before Fill Holes");
 422     typedef clitk::FillMaskFilter<MaskImageType> FillMaskFilterType;
 423     typename FillMaskFilterType::Pointer fillMaskFilter = FillMaskFilterType::New();
 424     fillMaskFilter->SetInput(working_mask);
 425     fillMaskFilter->AddDirection(2);
 426     //fillMaskFilter->AddDirection(1);
 427     fillMaskFilter->Update();
 428     working_mask = fillMaskFilter->GetOutput();
 429     StopCurrentStep<MaskImageType>(working_mask);
 430   }
 431
 432   if (GetSeparateLungsFlag()) {
 433     //--------------------------------------------------------------------
 434     //--------------------------------------------------------------------
 435     StartNewStep("Separate Left/Right lungs");
 436       PrintMemory(GetVerboseMemoryFlag(), "Before Separate");
 437     // Initial label
 438     working_mask = Labelize<MaskImageType>(working_mask,
 439                                           GetBackgroundValue(),
 440                                           false,
 441                                           GetMinimalComponentSize());
 442
 443     PrintMemory(GetVerboseMemoryFlag(), "After Labelize");
 444
 445     // Count the labels
 446     typedef itk::StatisticsImageFilter<MaskImageType> StatisticsImageFilterType;
 447     typename StatisticsImageFilterType::Pointer statisticsImageFilter=StatisticsImageFilterType::New();
 448     statisticsImageFilter->SetInput(working_mask);
 449     statisticsImageFilter->Update();
 450     unsigned int initialNumberOfLabels = statisticsImageFilter->GetMaximum();
 451     working_mask = statisticsImageFilter->GetOutput();
 452
 453     PrintMemory(GetVerboseMemoryFlag(), "After count label");
 454
 455     // Decompose the first label
 456     if (initialNumberOfLabels<2) {
 457       // Structuring element radius
 458       typename ImageType::SizeType radius;
 459       for (unsigned int i=0;i<Dim;i++) radius[i]=1;
 460       typedef clitk::DecomposeAndReconstructImageFilter<MaskImageType,MaskImageType> DecomposeAndReconstructFilterType;
 461       typename DecomposeAndReconstructFilterType::Pointer decomposeAndReconstructFilter=DecomposeAndReconstructFilterType::New();
 462       decomposeAndReconstructFilter->SetInput(working_mask);
 463       decomposeAndReconstructFilter->SetVerbose(false);
 464       decomposeAndReconstructFilter->SetRadius(radius);
 465       decomposeAndReconstructFilter->SetMaximumNumberOfLabels(2);
 466       decomposeAndReconstructFilter->SetMinimumObjectSize(this->GetMinimalComponentSize());
 467       decomposeAndReconstructFilter->SetMinimumNumberOfIterations(1);
 468       decomposeAndReconstructFilter->SetBackgroundValue(this->GetBackgroundValue());
 469       decomposeAndReconstructFilter->SetForegroundValue(this->GetForegroundValue());
 470       decomposeAndReconstructFilter->SetFullyConnected(true);
 471       decomposeAndReconstructFilter->SetNumberOfNewLabels(1);
 472       decomposeAndReconstructFilter->Update();
 473       working_mask = decomposeAndReconstructFilter->GetOutput();
 474     }
 475     PrintMemory(GetVerboseMemoryFlag(), "After decomposeAndReconstructFilter");
 476   }
 477
 478   // Retain labels ('1' is largset lung, so right. '2' is left)
 479   typedef itk::ThresholdImageFilter<MaskImageType> ThresholdImageFilterType;
 480   typename ThresholdImageFilterType::Pointer thresholdFilter = ThresholdImageFilterType::New();
 481   thresholdFilter->SetInput(working_mask);
 482   thresholdFilter->ThresholdAbove(2);
 483   thresholdFilter->SetOutsideValue(this->GetBackgroundValue());
 484   thresholdFilter->Update();
 485   working_mask = thresholdFilter->GetOutput();
 486   StopCurrentStep<MaskImageType> (working_mask);
 487   PrintMemory(GetVerboseMemoryFlag(), "After Thresholdfilter");
 488
 489   // Update output info
 490   //  output = working_mask;
 491   //this->GetOutput(0)->SetRegions(output->GetLargestPossibleRegion());
 492
 493   //  this->GetOutput(0)->SetRegions(working_mask->GetLargestPossibleRegion());
 494
 495 }
 496 //--------------------------------------------------------------------
 497
 498
 499 //--------------------------------------------------------------------
 500 template <class TImageType>
 501 void
 502 clitk::ExtractLungFilter<TImageType>::
 503 GenerateInputRequestedRegion() {
 504   //  DD("GenerateInputRequestedRegion (nothing?)");
 505 }
 506 //--------------------------------------------------------------------
 507
 508
 509 //--------------------------------------------------------------------
 510 template <class ImageType>
 511 void
 512 clitk::ExtractLungFilter<ImageType>::
 513 GenerateData()
 514 {
 515   // Set the output
 516   //  this->GraftOutput(output); // not SetNthOutput
 517   this->GraftOutput(working_mask); // not SetNthOutput
 518   // Store image filenames into AFDB
 519   GetAFDB()->SetImageFilename("Lungs", this->GetOutputLungFilename());
 520   GetAFDB()->SetImageFilename("Trachea", this->GetOutputTracheaFilename());
 521   WriteAFDB();
 522 }
 523 //--------------------------------------------------------------------
 524
 525
 526 //--------------------------------------------------------------------
 527 template <class ImageType>
 528 bool
 529 clitk::ExtractLungFilter<ImageType>::
 530 SearchForTracheaSeed(int skip)
 531 {
 532   if (m_Seeds.size() == 0) { // try to find seed (if not zero, it is given by user)
 533     // Restart the search until a seed is found, skipping more and more slices
 534     bool stop = false;
 535     while (!stop) {
 536       // Search seed (parameters = UpperThresholdForTrachea)
 537       static const unsigned int Dim = ImageType::ImageDimension;
 538       typename InternalImageType::RegionType sliceRegion = working_input->GetLargestPossibleRegion();
 539       typename InternalImageType::SizeType sliceRegionSize = sliceRegion.GetSize();
 540       typename InternalImageType::IndexType sliceRegionIndex = sliceRegion.GetIndex();
 541       sliceRegionIndex[Dim-1]=sliceRegionSize[Dim-1]-skip-5;
 542       sliceRegionSize[Dim-1]=5;
 543       sliceRegion.SetSize(sliceRegionSize);
 544       sliceRegion.SetIndex(sliceRegionIndex);
 545       typedef  itk::ImageRegionConstIterator<ImageType> IteratorType;
 546       IteratorType it(working_input, sliceRegion);
 547       it.GoToBegin();
 548       while (!it.IsAtEnd()) {
 549         if(it.Get() < GetUpperThresholdForTrachea() ) {
 550           AddSeed(it.GetIndex());
 551           // DD(it.GetIndex());
 552         }
 553         ++it;
 554       }
 555
 556       // if we do not found : restart
 557       stop = (m_Seeds.size() != 0);
 558       if (!stop) {
 559         if (GetVerboseStepFlag()) {
 560           std::cout << "\t No seed found this time. I skip some slices and restart." << std::endl;
 561         }
 562         if (skip > 0.5 * working_input->GetLargestPossibleRegion().GetSize()[2]) {
 563           // we want to skip more than a half of the image, it is probably a bug
 564           std::cerr << "Number of slices to skip to find trachea to high = " << skip << std::endl;
 565           stop = true;
 566         }
 567         skip += 5;
 568       }
 569       else {
 570         // DD(m_Seeds[0]);
 571         // DD(m_Seeds.size());
 572       }
 573     }
 574   }
 575   return (m_Seeds.size() != 0);
 576 }
 577 //--------------------------------------------------------------------
 578
 579
 580 bool is_orientation_superior(itk::SpatialOrientation::ValidCoordinateOrientationFlags orientation)
 581 {
 582   itk::SpatialOrientation::CoordinateTerms sup = itk::SpatialOrientation::ITK_COORDINATE_Superior;
 583   bool primary = (orientation & 0x0000ff) == sup;
 584   bool secondary = ((orientation & 0x00ff00) >> 8) == sup;
 585   bool tertiary = ((orientation & 0xff0000) >> 16) == sup;
 586   return primary || secondary || tertiary;
 587 }
 588
 589 //--------------------------------------------------------------------
 590 template <class ImageType>
 591 bool
 592 clitk::ExtractLungFilter<ImageType>::
 593 SearchForTracheaSeed2(int numberOfSlices)
 594 {
 595   if (m_Seeds.size() == 0) { // try to find seed (if not zero, it is given by user)
 596     if (GetVerboseRegionGrowingFlag())
 597       std::cout << "SearchForTracheaSeed2(" << numberOfSlices << ", " << GetMaxElongation() << ")" << std::endl;
 598
 599     typedef unsigned char MaskPixelType;
 600     typedef itk::Image<MaskPixelType, ImageType::ImageDimension> MaskImageType;
 601     typedef itk::Image<typename MaskImageType::PixelType, 2> MaskImageType2D;
 602     typedef itk::BinaryThresholdImageFilter<ImageType, MaskImageType> ThresholdFilterType;
 603     typedef itk::BinaryBallStructuringElement<MaskPixelType, MaskImageType::ImageDimension> KernelType;
 604     typedef itk::BinaryMorphologicalClosingImageFilter<MaskImageType, MaskImageType, KernelType> ClosingFilterType;
 605     typedef itk::BinaryMorphologicalOpeningImageFilter<MaskImageType, MaskImageType, KernelType> OpeningFilterType;
 606     typedef itk::ExtractImageFilter<MaskImageType, MaskImageType2D> SlicerFilterType;
 607     typedef itk::ConnectedComponentImageFilter<MaskImageType2D, MaskImageType2D> LabelFilterType;
 608     typedef itk::ShapeLabelObject<MaskPixelType, MaskImageType2D::ImageDimension> ShapeLabelType;
 609     typedef itk::LabelMap<ShapeLabelType> LabelImageType;
 610     typedef itk::LabelImageToShapeLabelMapFilter<MaskImageType2D, LabelImageType> ImageToLabelMapFilterType;
 611     typedef itk::LabelMapToLabelImageFilter<LabelImageType, MaskImageType2D> LabelMapToImageFilterType;
 612     typedef itk::ImageFileWriter<MaskImageType2D> FileWriterType;
 613
 614     // threshold to isolate airawys and lungs
 615     typename ThresholdFilterType::Pointer threshold = ThresholdFilterType::New();
 616     threshold->SetLowerThreshold(-2000);
 617     threshold->SetUpperThreshold(GetSeedPreProcessingThreshold());
 618     threshold->SetInput(working_input);
 619     threshold->Update();
 620
 621     KernelType kernel_closing, kernel_opening;
 622
 623     // remove small noise
 624     typename OpeningFilterType::Pointer opening = OpeningFilterType::New();
 625     kernel_opening.SetRadius(1);
 626     opening->SetKernel(kernel_opening);
 627     opening->SetInput(threshold->GetOutput());
 628     opening->Update();
 629
 630     typename SlicerFilterType::Pointer slicer = SlicerFilterType::New();
 631     slicer->SetInput(opening->GetOutput());
 632
 633     // label result
 634     typename LabelFilterType::Pointer label_filter = LabelFilterType::New();
 635     label_filter->SetInput(slicer->GetOutput());
 636
 637     // extract shape information from labels
 638     typename ImageToLabelMapFilterType::Pointer label_to_map_filter = ImageToLabelMapFilterType::New();
 639     label_to_map_filter->SetInput(label_filter->GetOutput());
 640
 641     typename LabelMapToImageFilterType::Pointer map_to_label_filter = LabelMapToImageFilterType::New();
 642     typename FileWriterType::Pointer writer = FileWriterType::New();
 643
 644     typename ImageType::IndexType index;
 645     typename ImageType::RegionType region = working_input->GetLargestPossibleRegion();
 646     typename ImageType::SizeType size = region.GetSize();
 647     typename ImageType::SpacingType spacing = working_input->GetSpacing();
 648     typename ImageType::PointType origin = working_input->GetOrigin();
 649
 650     int nslices = min(numberOfSlices, size[2]);
 651     int start = 0, increment = 1;
 652     itk::SpatialOrientationAdapter orientation;
 653     typename ImageType::DirectionType dir = working_input->GetDirection();
 654     if (!is_orientation_superior(orientation.FromDirectionCosines(dir))) {
 655       start = size[2]-1;
 656       increment = -1;
 657     }
 658
 659     typename MaskImageType::PointType image_centre;
 660     image_centre[0] = size[0]/2;
 661     image_centre[1] = size[1]/2;
 662     image_centre[2] = 0;
 663
 664     typedef InternalIndexType SeedType;
 665     SeedType trachea_centre, shape_centre, max_e_centre, prev_e_centre;
 666     typedef std::list<SeedType> PointListType;
 667     typedef std::list<PointListType> SequenceListType;
 668     PointListType* current_sequence = NULL;
 669     SequenceListType sequence_list;
 670
 671     prev_e_centre.Fill(0);
 672     std::ostringstream file_name;
 673     index[0] = index[1] = 0;
 674     size[0] = size[1] = 512;
 675     size[2] = 0;
 676     while (nslices--) {
 677       index[2] = start;
 678       start += increment;
 679
 680       region.SetIndex(index);
 681       region.SetSize(size);
 682       slicer->SetExtractionRegion(region);
 683       slicer->Update();
 684       label_filter->SetInput(slicer->GetOutput());
 685       label_filter->Update();
 686
 687       label_to_map_filter->SetInput(label_filter->GetOutput());
 688       label_to_map_filter->Update();
 689       typename LabelImageType::Pointer label_map = label_to_map_filter->GetOutput();
 690
 691       if (GetWriteStepFlag()) {
 692         map_to_label_filter->SetInput(label_map);
 693         writer->SetInput(map_to_label_filter->GetOutput());
 694         file_name.str("");
 695         file_name << "labels_";
 696         file_name.width(3);
 697         file_name.fill('0');
 698         file_name << index[2] << ".mhd";
 699         writer->SetFileName(file_name.str().c_str());
 700         writer->Update();
 701       }
 702
 703       typename ShapeLabelType::Pointer shape, max_e_shape;
 704       double max_elongation = GetMaxElongation();
 705       double max_size = size[0]*size[1]/128;
 706       double max_e = 0;
 707       int nshapes = 0;
 708       unsigned int nlables = label_map->GetNumberOfLabelObjects();
 709       for (unsigned int j = 0; j < nlables; j++) {
 710         shape = label_map->GetNthLabelObject(j);
 711         if (shape->Size() > 150 && shape->Size() <= max_size) {
 712 #if ITK_VERSION_MAJOR < 4
 713           double e = 1 - 1/shape->GetBinaryElongation();
 714 #else
 715           double e = 1 - 1/shape->GetElongation();
 716 #endif
 717           //double area = 1 - r->Area() ;
 718           if (e < max_elongation) {
 719             nshapes++;
 720             shape_centre[0] = (shape->GetCentroid()[0] - origin[0])/spacing[0];
 721             shape_centre[1] = (shape->GetCentroid()[1] - origin[1])/spacing[1];
 722             shape_centre[2] = index[2];
 723             //double d = 1 - (shape_centre - image_centre).Magnitude()/max_dist;
 724             double dx = shape_centre[0] - image_centre[0];
 725             double d = 1 - dx*2/size[0];
 726             e = e + d;
 727             if (e > max_e)
 728             {
 729               max_e = e;
 730               max_e_shape = shape;
 731               max_e_centre = shape_centre;
 732             }
 733           }
 734         }
 735       }
 736
 737       if (nshapes > 0)
 738       {
 739         itk::Point<typename SeedType::IndexValueType, ImageType::ImageDimension> p1, p2;
 740         p1[0] = max_e_centre[0];
 741         p1[1] = max_e_centre[1];
 742         p1[2] = max_e_centre[2];
 743
 744         p2[0] = prev_e_centre[0];
 745         p2[1] = prev_e_centre[1];
 746         p2[2] = prev_e_centre[2];
 747
 748         double mag = (p2 - p1).GetNorm();
 749         if (GetVerboseRegionGrowingFlag()) {
 750           cout.precision(3);
 751           cout << index[2] << ": ";
 752           cout << "region(" << max_e_centre[0] << ", " << max_e_centre[1] << ", " << max_e_centre[2] << "); ";
 753           cout << "prev_region(" << prev_e_centre[0] << ", " << prev_e_centre[1] << ", " << prev_e_centre[2] << "); ";
 754           cout << "mag(" << mag << "); " << endl;
 755         }
 756
 757         if (mag > 5)
 758         {
 759           PointListType point_list;
 760           point_list.push_back(max_e_centre);
 761           sequence_list.push_back(point_list);
 762           current_sequence = &sequence_list.back();
 763         }
 764         else if (current_sequence)
 765           current_sequence->push_back(max_e_centre);
 766
 767         prev_e_centre= max_e_centre;
 768       }
 769     }
 770
 771     size_t longest = 0;
 772     for (typename SequenceListType::iterator s = sequence_list.begin(); s != sequence_list.end(); s++)
 773     {
 774       if (s->size() > longest)
 775       {
 776         longest = s->size();
 777         trachea_centre = s->front();
 778       }
 779     }
 780
 781     if (GetVerboseRegionGrowingFlag())
 782       std::cout << "seed at: " << trachea_centre << std::endl;
 783     m_Seeds.push_back(trachea_centre);
 784   }
 785
 786   return (m_Seeds.size() != 0);
 787 }
 788 //--------------------------------------------------------------------
 789
 790 //--------------------------------------------------------------------
 791 template <class ImageType>
 792 void
 793 clitk::ExtractLungFilter<ImageType>::
 794 TracheaRegionGrowing()
 795 {
 796   // Explosion controlled region growing
 797   PrintMemory(GetVerboseMemoryFlag(), "Before ExplosionControlledThresholdConnectedImageFilter");
 798   typedef clitk::ExplosionControlledThresholdConnectedImageFilter<ImageType, MaskImageType> ImageFilterType;
 799   typename ImageFilterType::Pointer f= ImageFilterType::New();
 800   f->SetInput(working_input);
 801   f->SetLower(-2000);
 802   f->SetUpper(GetUpperThresholdForTrachea());
 803   f->SetMinimumLowerThreshold(-2000);
 804   //  f->SetMaximumUpperThreshold(0); // MAYBE TO CHANGE ???
 805   f->SetMaximumUpperThreshold(-700); // MAYBE TO CHANGE ???
 806   f->SetAdaptLowerBorder(false);
 807   f->SetAdaptUpperBorder(true);
 808   f->SetMinimumSize(5000);
 809   f->SetReplaceValue(1);
 810   f->SetMultiplier(GetMultiplierForTrachea());
 811   f->SetThresholdStepSize(GetThresholdStepSizeForTrachea());
 812   f->SetMinimumThresholdStepSize(1);
 813   f->SetVerbose(GetVerboseRegionGrowingFlag());
 814   for(unsigned int i=0; i<m_Seeds.size();i++) {
 815     f->AddSeed(m_Seeds[i]);
 816     // DD(m_Seeds[i]);
 817   }
 818   f->Update();
 819   PrintMemory(GetVerboseMemoryFlag(), "After RG update");
 820
 821   // take first (main) connected component
 822   trachea = Labelize<MaskImageType>(f->GetOutput(),
 823                                     GetBackgroundValue(),
 824                                     true,
 825                                     1000);//GetMinimalComponentSize());
 826   PrintMemory(GetVerboseMemoryFlag(), "After Labelize");
 827   trachea = KeepLabels<MaskImageType>(trachea,
 828                                               GetBackgroundValue(),
 829                                               GetForegroundValue(),
 830                                               1, 1, false);
 831   PrintMemory(GetVerboseMemoryFlag(), "After KeepLabels");
 832 }
 833 //--------------------------------------------------------------------
 834
 835
 836 //--------------------------------------------------------------------
 837 template <class ImageType>
 838 double
 839 clitk::ExtractLungFilter<ImageType>::
 840 ComputeTracheaVolume()
 841 {
 842   typedef itk::ImageRegionConstIterator<InternalImageType> IteratorType;
 843   IteratorType iter(trachea, trachea->GetLargestPossibleRegion());
 844   iter.GoToBegin();
 845   double volume = 0.0;
 846   while (!iter.IsAtEnd()) {
 847     if (iter.Get() == this->GetForegroundValue()) volume++;
 848     ++iter;
 849   }
 850
 851   double voxelsize = trachea->GetSpacing()[0]*trachea->GetSpacing()[1]*trachea->GetSpacing()[2];
 852   return volume*voxelsize;
 853 }
 854 //--------------------------------------------------------------------
 855
 856
 857 //--------------------------------------------------------------------
 858 template <class ImageType>
 859 void
 860 clitk::ExtractLungFilter<ImageType>::
 861 SearchForTrachea()
 862 {
 863   // Search for seed among n slices, skip some slices before starting
 864   // if not found -> skip more and restart
 865
 866   // when seed found : perform region growing
 867   // compute trachea volume
 868   // if volume not plausible  -> skip more slices and restart
 869
 870   bool has_seed;
 871   bool stop = false;
 872   double volume = 0.0;
 873   int skip = GetNumberOfSlicesToSkipBeforeSearchingSeed();
 874   while (!stop) {
 875     stop = true;
 876     if (GetTracheaSeedAlgorithm() == 0)
 877       has_seed = SearchForTracheaSeed(skip);
 878     else
 879       has_seed = SearchForTracheaSeed2(GetNumSlices());
 880
 881     if (has_seed) {
 882       TracheaRegionGrowing();
 883       volume = ComputeTracheaVolume()/1000; // assume mm3, so divide by 1000 to get cc
 884       if (GetWriteStepFlag()) {
 885         writeImage<MaskImageType>(trachea, "step-trachea-"+toString(skip)+".mhd");
 886       }
 887       if (GetTracheaVolumeMustBeCheckedFlag()) {
 888         if ((volume > 10) && (volume < 65 )) { // depend on image size ...
 889           // Typical volume 22.59 cm 3 (± 7.69 cm 3 ) [ Leader 2004 ]
 890           if (GetVerboseStepFlag())
 891           {
 892             std::cout << "\t Found trachea with volume " << volume << " cc." << std::endl;
 893           }
 894         }
 895         else {
 896           if (GetVerboseStepFlag()) {
 897             std::cout << "\t The volume of the trachea (" << volume
 898                       << " cc) seems not correct. I skip some slices (" << skip << ") and restart to find seeds."
 899                       << std::endl;
 900           }
 901           skip += 5;
 902           stop = false;
 903           // empty the list of seed
 904           m_Seeds.clear();
 905         }
 906         if (skip > 0.5 * working_input->GetLargestPossibleRegion().GetSize()[2]) {
 907           // we want to skip more than a half of the image, it is probably a bug
 908           std::cerr << "2 : Number of slices to skip to find trachea too high = " << skip << std::endl;
 909           stop = true;
 910         }
 911       }
 912       else {
 913         stop = true;
 914       }
 915     }
 916   }
 917
 918   if (volume != 0.0) {
 919     // Set output
 920     StopCurrentStep<MaskImageType>(trachea);
 921   }
 922   else { // Trachea not found
 923     this->SetWarning("* WARNING * No seed found for trachea.");
 924     StopCurrentStep();
 925   }
 926 }
 927 //--------------------------------------------------------------------
 928
 929 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX