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