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://oncora1.lyon.fnclcc.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
  28 // itk
  29 #include "itkBinaryThresholdImageFilter.h"
  30 #include "itkConnectedComponentImageFilter.h"
  31 #include "itkRelabelComponentImageFilter.h"
  32 #include "itkOtsuThresholdImageFilter.h"
  33
  34 //--------------------------------------------------------------------
  35 template <class TInputImageType, class TMaskImageType>
  36 clitk::ExtractLungFilter<TInputImageType, TMaskImageType>::
  37 ExtractLungFilter():
  38   clitk::FilterBase(),
  39   itk::ImageToImageFilter<TInputImageType, TMaskImageType>()
  40 {
  41   SetNumberOfSteps(10);
  42   // Default global options
  43   this->SetNumberOfRequiredInputs(2);
  44   SetPatientMaskBackgroundValue(0);
  45   SetBackgroundValue(0); // Must be zero
  46   SetForegroundValue(1);
  47
  48   // Step 1 default values
  49   SetUpperThreshold(-300);
  50   SetLowerThreshold(-1000);
  51   UseLowerThresholdOff();
  52   LabelParamType * p1 = new LabelParamType;
  53   p1->SetFirstKeep(1);
  54   p1->UseLastKeepOff();
  55   p1->AddLabelToRemove(2);
  56   SetLabelizeParameters1(p1);
  57
  58   // Step 2 default values
  59   SetUpperThresholdForTrachea(-900);
  60   SetMultiplierForTrachea(5);
  61   SetThresholdStepSizeForTrachea(64);
  62
  63   // Step 3 default values
  64   SetNumberOfHistogramBins(500);
  65   LabelParamType * p2 = new LabelParamType;
  66   p2->SetFirstKeep(1);
  67   p2->UseLastKeepOff();
  68   // p->AddLabelToRemove(2); // No label to remove by default
  69   SetLabelizeParameters2(p2);
  70
  71   // Step 4 default values
  72   SetRadiusForTrachea(1);
  73   LabelParamType * p3 = new LabelParamType;
  74   p3->SetFirstKeep(1);
  75   p3->SetLastKeep(2);
  76   p3->UseLastKeepOff();
  77   SetLabelizeParameters3(p3);
  78 }
  79 //--------------------------------------------------------------------
  80
  81
  82 //--------------------------------------------------------------------
  83 template <class TInputImageType, class TMaskImageType>
  84 void
  85 clitk::ExtractLungFilter<TInputImageType, TMaskImageType>::
  86 SetInput(const TInputImageType * image)
  87 {
  88   this->SetNthInput(0, const_cast<TInputImageType *>(image));
  89 }
  90 //--------------------------------------------------------------------
  91
  92
  93 //--------------------------------------------------------------------
  94 template <class TInputImageType, class TMaskImageType>
  95 void
  96 clitk::ExtractLungFilter<TInputImageType, TMaskImageType>::
  97 SetInputPatientMask(TMaskImageType * image, MaskImagePixelType bg )
  98 {
  99   this->SetNthInput(1, const_cast<TMaskImageType *>(image));
 100   SetPatientMaskBackgroundValue(bg);
 101 }
 102 //--------------------------------------------------------------------
 103
 104
 105 //--------------------------------------------------------------------
 106 template <class TInputImageType, class TMaskImageType>
 107 void
 108 clitk::ExtractLungFilter<TInputImageType, TMaskImageType>::
 109 AddSeed(InternalIndexType s)
 110 {
 111   m_Seeds.push_back(s);
 112 }
 113 //--------------------------------------------------------------------
 114
 115
 116 //--------------------------------------------------------------------
 117 template <class TInputImageType, class TMaskImageType>
 118 template<class ArgsInfoType>
 119 void
 120 clitk::ExtractLungFilter<TInputImageType, TMaskImageType>::
 121 SetArgsInfo(ArgsInfoType mArgsInfo)
 122 {
 123   SetVerboseOption_GGO(mArgsInfo);
 124   SetVerboseStep_GGO(mArgsInfo);
 125   SetWriteStep_GGO(mArgsInfo);
 126   SetVerboseWarningOff_GGO(mArgsInfo);
 127
 128   SetUpperThreshold_GGO(mArgsInfo);
 129   SetLowerThreshold_GGO(mArgsInfo);
 130   SetLabelizeParameters1_GGO(mArgsInfo);
 131   if (!mArgsInfo.remove1_given) {
 132     GetLabelizeParameters1()->AddLabelToRemove(2);
 133     if (GetVerboseOption()) VerboseOption("remove1", 2);
 134   }
 135
 136   SetUpperThresholdForTrachea_GGO(mArgsInfo);
 137   SetMultiplierForTrachea_GGO(mArgsInfo);
 138   SetThresholdStepSizeForTrachea_GGO(mArgsInfo);
 139   AddSeed_GGO(mArgsInfo);
 140
 141   SetNumberOfHistogramBins_GGO(mArgsInfo);
 142   SetLabelizeParameters2_GGO(mArgsInfo);
 143
 144   SetRadiusForTrachea_GGO(mArgsInfo);
 145   SetLabelizeParameters3_GGO(mArgsInfo);
 146 }
 147 //--------------------------------------------------------------------
 148
 149
 150 //--------------------------------------------------------------------
 151 template <class TInputImageType, class TMaskImageType>
 152 void
 153 clitk::ExtractLungFilter<TInputImageType, TMaskImageType>::
 154 GenerateOutputInformation()
 155 {
 156   input   = dynamic_cast<const TInputImageType*>(itk::ProcessObject::GetInput(0));
 157   Superclass::GenerateOutputInformation();
 158
 159   // Get input pointers
 160   input   = dynamic_cast<const TInputImageType*>(itk::ProcessObject::GetInput(0));
 161   patient = dynamic_cast<const TMaskImageType*>(itk::ProcessObject::GetInput(1));
 162
 163   // Check image
 164   if (!HaveSameSizeAndSpacing<TInputImageType, TMaskImageType>(input, patient)) {
 165     this->SetLastError("* ERROR * the images (input and patient mask) must have the same size & spacing");
 166     return;
 167   }
 168
 169   // Set Number of steps
 170   SetNumberOfSteps(9);
 171
 172   //--------------------------------------------------------------------
 173   //--------------------------------------------------------------------
 174   StartNewStep("Set background to initial image");
 175   working_input = SetBackground<TInputImageType, TMaskImageType>
 176     (input, patient, GetPatientMaskBackgroundValue(), -1000);
 177   StopCurrentStep<InputImageType>(working_input);
 178
 179   //--------------------------------------------------------------------
 180   //--------------------------------------------------------------------
 181   StartNewStep("Remove Air");
 182   // Threshold to get air
 183   typedef itk::BinaryThresholdImageFilter<InputImageType, InternalImageType> InputBinarizeFilterType;
 184   typename InputBinarizeFilterType::Pointer binarizeFilter=InputBinarizeFilterType::New();
 185   binarizeFilter->SetInput(working_input);
 186   if (m_UseLowerThreshold) binarizeFilter->SetLowerThreshold(m_LowerThreshold);
 187   binarizeFilter->SetUpperThreshold(m_UpperThreshold);
 188   binarizeFilter ->SetInsideValue(this->GetForegroundValue());
 189   binarizeFilter ->SetOutsideValue(this->GetBackgroundValue());
 190   binarizeFilter->Update();
 191   working_image = binarizeFilter->GetOutput();
 192
 193   // Labelize and keep right labels
 194   working_image = Labelize<InternalImageType>(working_image, GetBackgroundValue(), true, GetMinimalComponentSize());
 195   working_image = RemoveLabels<InternalImageType>
 196     (working_image, GetBackgroundValue(), GetLabelizeParameters1()->GetLabelsToRemove());
 197   typename InternalImageType::Pointer air = KeepLabels<InternalImageType>
 198     (working_image,
 199      GetBackgroundValue(),
 200      GetForegroundValue(),
 201      GetLabelizeParameters1()->GetFirstKeep(),
 202      GetLabelizeParameters1()->GetLastKeep(),
 203      GetLabelizeParameters1()->GetUseLastKeep());
 204
 205   // Set Air to BG
 206   working_input = SetBackground<TInputImageType, InternalImageType>
 207     (working_input, air, this->GetForegroundValue(), this->GetBackgroundValue());
 208
 209   // End
 210   StopCurrentStep<InputImageType>(working_input);
 211
 212   //--------------------------------------------------------------------
 213   //--------------------------------------------------------------------
 214   StartNewStep("Find the trachea");
 215   //DD(m_Seeds.size());
 216   if (m_Seeds.size() == 0) { // try to find seed
 217     // Search seed (parameters = UpperThresholdForTrachea)
 218     static const unsigned int Dim = InputImageType::ImageDimension;
 219     typename InternalImageType::RegionType sliceRegion = working_input->GetLargestPossibleRegion();
 220     typename InternalImageType::SizeType sliceRegionSize = sliceRegion.GetSize();
 221     typename InternalImageType::IndexType sliceRegionIndex = sliceRegion.GetIndex();
 222     sliceRegionIndex[Dim-1]=sliceRegionSize[Dim-1]-5;
 223     sliceRegionSize[Dim-1]=5;
 224     sliceRegion.SetSize(sliceRegionSize);
 225     sliceRegion.SetIndex(sliceRegionIndex);
 226     //DD(GetUpperThresholdForTrachea());
 227     //DD(sliceRegion);
 228     typedef  itk::ImageRegionConstIterator<InputImageType> IteratorType;
 229     IteratorType it(working_input, sliceRegion);
 230     it.GoToBegin();
 231     while (!it.IsAtEnd()) {
 232       if(it.Get() < GetUpperThresholdForTrachea() ) {
 233         AddSeed(it.GetIndex());
 234         //      DD(it.GetIndex());
 235       }
 236       ++it;
 237     }
 238   }
 239
 240   //DD(m_Seeds.size());
 241   if (m_Seeds.size() != 0) {
 242     // Explosion controlled region growing
 243     typedef clitk::ExplosionControlledThresholdConnectedImageFilter<InputImageType, InternalImageType> ImageFilterType;
 244     typename ImageFilterType::Pointer f= ImageFilterType::New();
 245     f->SetInput(working_input);
 246     f->SetVerbose(false);
 247     f->SetLower(-2000);
 248     f->SetUpper(GetUpperThresholdForTrachea());
 249     f->SetMinimumLowerThreshold(-2000);
 250     f->SetMaximumUpperThreshold(0);
 251     f->SetAdaptLowerBorder(false);
 252     f->SetAdaptUpperBorder(true);
 253     f->SetMinimumSize(5000);
 254     f->SetReplaceValue(1);
 255     f->SetMultiplier(GetMultiplierForTrachea());
 256     f->SetThresholdStepSize(GetThresholdStepSizeForTrachea());
 257     f->SetMinimumThresholdStepSize(1);
 258     for(unsigned int i=0; i<m_Seeds.size();i++) {
 259       //    std::cout<<"Adding seed " <<m_Seeds[i]<<"..."<<std::endl;
 260       f->AddSeed(m_Seeds[i]);
 261     }
 262     f->Update();
 263     trachea_tmp = f->GetOutput();
 264     // Set output
 265     StopCurrentStep<InternalImageType>(trachea_tmp);
 266
 267   }
 268   else { // Trachea not found
 269     this->SetWarning("* WARNING * No seed found for trachea.");
 270     // Set output
 271     StopCurrentStep();
 272   }
 273
 274   //--------------------------------------------------------------------
 275   //--------------------------------------------------------------------
 276   StartNewStep("Extract the lung with Otsu filter");
 277   // Automated Otsu thresholding and relabeling
 278   typedef itk::OtsuThresholdImageFilter<InputImageType,InternalImageType> OtsuThresholdImageFilterType;
 279   typename OtsuThresholdImageFilterType::Pointer otsuFilter=OtsuThresholdImageFilterType::New();
 280   otsuFilter->SetInput(working_input);
 281   otsuFilter->SetNumberOfHistogramBins(GetNumberOfHistogramBins());
 282   otsuFilter->SetInsideValue(this->GetForegroundValue());
 283   otsuFilter->SetOutsideValue(this->GetBackgroundValue());
 284   otsuFilter->Update();
 285   working_image =  otsuFilter->GetOutput();
 286
 287   // Set output
 288   StopCurrentStep<InternalImageType>(working_image);
 289
 290   //--------------------------------------------------------------------
 291   //--------------------------------------------------------------------
 292   StartNewStep("Select labels");
 293   // Keep right labels
 294   working_image = LabelizeAndSelectLabels<InternalImageType>
 295     (working_image,
 296      GetBackgroundValue(),
 297      GetForegroundValue(),
 298      false,
 299      GetMinimalComponentSize(),
 300      GetLabelizeParameters2());
 301
 302   // Set output
 303   StopCurrentStep<InternalImageType>(working_image);
 304
 305   //--------------------------------------------------------------------
 306   //--------------------------------------------------------------------
 307   if (m_Seeds.size() != 0) { // if ==0 ->no trachea found
 308     StartNewStep("Remove the trachea");
 309     // Set the trachea
 310     working_image = SetBackground<InternalImageType, InternalImageType>
 311       (working_image, trachea_tmp, 1, -1);
 312
 313    // Dilate the trachea
 314     static const unsigned int Dim = InputImageType::ImageDimension;
 315     typedef itk::BinaryBallStructuringElement<InternalPixelType, Dim> KernelType;
 316     KernelType structuringElement;
 317     structuringElement.SetRadius(GetRadiusForTrachea());
 318     structuringElement.CreateStructuringElement();
 319     typedef clitk::ConditionalBinaryDilateImageFilter<InternalImageType, InternalImageType, KernelType> ConditionalBinaryDilateImageFilterType;
 320     typename ConditionalBinaryDilateImageFilterType::Pointer dilateFilter = ConditionalBinaryDilateImageFilterType::New();
 321     dilateFilter->SetBoundaryToForeground(false);
 322     dilateFilter->SetKernel(structuringElement);
 323     dilateFilter->SetBackgroundValue (1);
 324     dilateFilter->SetForegroundValue (-1);
 325     dilateFilter->SetInput (working_image);
 326     dilateFilter->Update();
 327     working_image = dilateFilter->GetOutput();
 328
 329     // Set trachea with dilatation
 330     trachea_tmp = SetBackground<InternalImageType, InternalImageType>
 331       (trachea_tmp, working_image, -1, this->GetForegroundValue());
 332
 333     // Remove the trachea
 334     working_image = SetBackground<InternalImageType, InternalImageType>
 335       (working_image, working_image, -1, this->GetBackgroundValue());
 336
 337     // Label
 338     working_image = LabelizeAndSelectLabels<InternalImageType>
 339       (working_image,
 340        GetBackgroundValue(),
 341        GetForegroundValue(),
 342        false,
 343        GetMinimalComponentSize(),
 344        GetLabelizeParameters3());
 345
 346     // Set output
 347     StopCurrentStep<InternalImageType>(working_image);
 348   }
 349
 350
 351   //--------------------------------------------------------------------
 352   //--------------------------------------------------------------------
 353   typedef clitk::AutoCropFilter<InternalImageType> CropFilterType;
 354   typename CropFilterType::Pointer cropFilter = CropFilterType::New();
 355   if (m_Seeds.size() != 0) { // if ==0 ->no trachea found
 356     StartNewStep("Croping trachea");
 357     cropFilter->SetInput(trachea_tmp);
 358     cropFilter->Update(); // Needed
 359     typedef itk::CastImageFilter<InternalImageType, MaskImageType> CastImageFilterType;
 360     typename CastImageFilterType::Pointer caster= CastImageFilterType::New();
 361     caster->SetInput(cropFilter->GetOutput());
 362     caster->Update();
 363     trachea = caster->GetOutput();
 364     StopCurrentStep<MaskImageType>(trachea);
 365   }
 366
 367   //--------------------------------------------------------------------
 368   //--------------------------------------------------------------------
 369   StartNewStep("Croping lung");
 370   typename CropFilterType::Pointer cropFilter2 = CropFilterType::New(); // Needed to reset pipeline
 371   cropFilter2->SetInput(working_image);
 372   cropFilter2->Update();
 373   working_image = cropFilter2->GetOutput();
 374   StopCurrentStep<InternalImageType>(working_image);
 375
 376   //--------------------------------------------------------------------
 377   //--------------------------------------------------------------------
 378   StartNewStep("Separate Left/Right lungs");
 379   // Initial label
 380   working_image = Labelize<InternalImageType>(working_image,
 381                                               GetBackgroundValue(),
 382                                               false,
 383                                               GetMinimalComponentSize());
 384
 385   // Count the labels
 386   typedef itk::StatisticsImageFilter<InternalImageType> StatisticsImageFilterType;
 387   typename StatisticsImageFilterType::Pointer statisticsImageFilter=StatisticsImageFilterType::New();
 388   statisticsImageFilter->SetInput(working_image);
 389   statisticsImageFilter->Update();
 390   unsigned int initialNumberOfLabels = statisticsImageFilter->GetMaximum();
 391   working_image = statisticsImageFilter->GetOutput();
 392
 393   // Decompose the first label
 394   static const unsigned int Dim = InputImageType::ImageDimension;
 395   if (initialNumberOfLabels<2) {
 396     // Structuring element radius
 397     typename InputImageType::SizeType radius;
 398     for (unsigned int i=0;i<Dim;i++) radius[i]=1;
 399     typedef clitk::DecomposeAndReconstructImageFilter<InternalImageType,InternalImageType> DecomposeAndReconstructFilterType;
 400     typename DecomposeAndReconstructFilterType::Pointer decomposeAndReconstructFilter=DecomposeAndReconstructFilterType::New();
 401     decomposeAndReconstructFilter->SetInput(working_image);
 402     decomposeAndReconstructFilter->SetVerbose(false);
 403     decomposeAndReconstructFilter->SetRadius(radius);
 404     decomposeAndReconstructFilter->SetMaximumNumberOfLabels(2);
 405     decomposeAndReconstructFilter->SetMinimumObjectSize(this->GetMinimalComponentSize());
 406     decomposeAndReconstructFilter->SetMinimumNumberOfIterations(1);
 407     decomposeAndReconstructFilter->SetBackgroundValue(this->GetBackgroundValue());
 408     decomposeAndReconstructFilter->SetForegroundValue(this->GetForegroundValue());
 409     decomposeAndReconstructFilter->SetFullyConnected(true);
 410     decomposeAndReconstructFilter->SetNumberOfNewLabels(1);
 411     decomposeAndReconstructFilter->Update();
 412     working_image = decomposeAndReconstructFilter->GetOutput();
 413   }
 414
 415   // Retain labels (lungs)
 416   typedef itk::ThresholdImageFilter<InternalImageType> ThresholdImageFilterType;
 417   typename ThresholdImageFilterType::Pointer thresholdFilter = ThresholdImageFilterType::New();
 418   thresholdFilter->SetInput(working_image);
 419   thresholdFilter->ThresholdAbove(2);
 420   thresholdFilter->SetOutsideValue(this->GetBackgroundValue());
 421   thresholdFilter->Update();
 422   working_image = thresholdFilter->GetOutput();
 423   StopCurrentStep<InternalImageType> (working_image);
 424 }
 425 //--------------------------------------------------------------------
 426
 427
 428 //--------------------------------------------------------------------
 429 template <class TInputImageType, class TMaskImageType>
 430 void
 431 clitk::ExtractLungFilter<TInputImageType, TMaskImageType>::
 432 GenerateData() {
 433   // Final Cast
 434   typedef itk::CastImageFilter<InternalImageType, MaskImageType> CastImageFilterType;
 435   typename CastImageFilterType::Pointer caster= CastImageFilterType::New();
 436   caster->SetInput(working_image);
 437   caster->Update();
 438   // Set output
 439   //this->SetNthOutput(0, caster->GetOutput()); // -> no because redo filter otherwise
 440   this->GraftOutput(caster->GetOutput());
 441   return;
 442 }
 443 //--------------------------------------------------------------------
 444
 445
 446 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX