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