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