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 #include "itkBinaryThinningImageFilter3D.h"
  34 #include "itkImageIteratorWithIndex.h"
  35
  36
  37 //--------------------------------------------------------------------
  38 template <class ImageType, class MaskImageType>
  39 clitk::ExtractLungFilter<ImageType, MaskImageType>::
  40 ExtractLungFilter():
  41   clitk::FilterBase(),
  42   itk::ImageToImageFilter<ImageType, MaskImageType>()
  43 {
  44   SetNumberOfSteps(10);
  45   m_MaxSeedNumber = 500;
  46
  47   // Default global options
  48   this->SetNumberOfRequiredInputs(2);
  49   SetPatientMaskBackgroundValue(0);
  50   SetBackgroundValue(0); // Must be zero
  51   SetForegroundValue(1);
  52   SetMinimalComponentSize(100);
  53
  54   // Step 1 default values
  55   SetUpperThreshold(-300);
  56   SetLowerThreshold(-1000);
  57   UseLowerThresholdOff();
  58   LabelParamType * p1 = new LabelParamType;
  59   p1->SetFirstKeep(1);
  60   p1->UseLastKeepOff();
  61   p1->AddLabelToRemove(2);
  62   SetLabelizeParameters1(p1);
  63
  64   // Step 2 default values
  65   SetUpperThresholdForTrachea(-900);
  66   SetMultiplierForTrachea(5);
  67   SetThresholdStepSizeForTrachea(64);
  68   SetNumberOfSlicesToSkipBeforeSearchingSeed(0);
  69
  70   // Step 3 default values
  71   SetNumberOfHistogramBins(500);
  72   LabelParamType * p2 = new LabelParamType;
  73   p2->SetFirstKeep(1);
  74   p2->UseLastKeepOff();
  75   // p->AddLabelToRemove(2); // No label to remove by default
  76   SetLabelizeParameters2(p2);
  77
  78   // Step 4 default values
  79   SetRadiusForTrachea(1);
  80   LabelParamType * p3 = new LabelParamType;
  81   p3->SetFirstKeep(1);
  82   p3->SetLastKeep(2);
  83   p3->UseLastKeepOff();
  84   SetLabelizeParameters3(p3);
  85
  86   // Step 5 : find bronchial bifurcations
  87   FindBronchialBifurcationsOn();
  88 }
  89 //--------------------------------------------------------------------
  90
  91
  92 //--------------------------------------------------------------------
  93 template <class ImageType, class MaskImageType>
  94 void
  95 clitk::ExtractLungFilter<ImageType, MaskImageType>::
  96 SetInput(const ImageType * image)
  97 {
  98   this->SetNthInput(0, const_cast<ImageType *>(image));
  99 }
 100 //--------------------------------------------------------------------
 101
 102
 103 //--------------------------------------------------------------------
 104 template <class ImageType, class MaskImageType>
 105 void
 106 clitk::ExtractLungFilter<ImageType, MaskImageType>::
 107 SetInputPatientMask(MaskImageType * image, MaskImagePixelType bg )
 108 {
 109   this->SetNthInput(1, const_cast<MaskImageType *>(image));
 110   SetPatientMaskBackgroundValue(bg);
 111 }
 112 //--------------------------------------------------------------------
 113
 114
 115 //--------------------------------------------------------------------
 116 template <class ImageType, class MaskImageType>
 117 void
 118 clitk::ExtractLungFilter<ImageType, MaskImageType>::
 119 AddSeed(InternalIndexType s)
 120 {
 121   m_Seeds.push_back(s);
 122 }
 123 //--------------------------------------------------------------------
 124
 125
 126 //--------------------------------------------------------------------
 127 template <class ImageType, class MaskImageType>
 128 template<class ArgsInfoType>
 129 void
 130 clitk::ExtractLungFilter<ImageType, MaskImageType>::
 131 SetArgsInfo(ArgsInfoType mArgsInfo)
 132 {
 133   SetVerboseOption_GGO(mArgsInfo);
 134   SetVerboseStep_GGO(mArgsInfo);
 135   SetWriteStep_GGO(mArgsInfo);
 136   SetVerboseWarningOff_GGO(mArgsInfo);
 137
 138   SetUpperThreshold_GGO(mArgsInfo);
 139   SetLowerThreshold_GGO(mArgsInfo);
 140   SetNumberOfSlicesToSkipBeforeSearchingSeed_GGO(mArgsInfo);
 141   SetLabelizeParameters1_GGO(mArgsInfo);
 142   if (!mArgsInfo.remove1_given) {
 143     GetLabelizeParameters1()->AddLabelToRemove(2);
 144     if (GetVerboseOption()) VerboseOption("remove1", 2);
 145   }
 146
 147   SetUpperThresholdForTrachea_GGO(mArgsInfo);
 148   SetMultiplierForTrachea_GGO(mArgsInfo);
 149   SetThresholdStepSizeForTrachea_GGO(mArgsInfo);
 150   AddSeed_GGO(mArgsInfo);
 151
 152   SetMinimalComponentSize_GGO(mArgsInfo);
 153
 154   SetNumberOfHistogramBins_GGO(mArgsInfo);
 155   SetLabelizeParameters2_GGO(mArgsInfo);
 156
 157   SetRadiusForTrachea_GGO(mArgsInfo);
 158   SetLabelizeParameters3_GGO(mArgsInfo);
 159 }
 160 //--------------------------------------------------------------------
 161
 162
 163 //--------------------------------------------------------------------
 164 template <class ImageType, class MaskImageType>
 165 void
 166 clitk::ExtractLungFilter<ImageType, MaskImageType>::
 167 GenerateOutputInformation()
 168 {
 169   Superclass::GenerateOutputInformation();
 170   //this->GetOutput(0)->SetRequestedRegion(this->GetOutput(0)->GetLargestPossibleRegion());
 171
 172   // Get input pointers
 173   patient = dynamic_cast<const MaskImageType*>(itk::ProcessObject::GetInput(1));
 174   input   = dynamic_cast<const ImageType*>(itk::ProcessObject::GetInput(0));
 175
 176   // Check image
 177   if (!HaveSameSizeAndSpacing<ImageType, MaskImageType>(input, patient)) {
 178     this->SetLastError("* ERROR * the images (input and patient mask) must have the same size & spacing");
 179     return;
 180   }
 181
 182   //--------------------------------------------------------------------
 183   //--------------------------------------------------------------------
 184   StartNewStepOrStop("Set background to initial image");
 185   working_input = SetBackground<ImageType, MaskImageType>
 186     (input, patient, GetPatientMaskBackgroundValue(), -1000);
 187   StopCurrentStep<ImageType>(working_input);
 188
 189   //--------------------------------------------------------------------
 190   //--------------------------------------------------------------------
 191   StartNewStepOrStop("Remove Air");
 192   // Check threshold
 193   if (m_UseLowerThreshold) {
 194     if (m_LowerThreshold > m_UpperThreshold) {
 195       this->SetLastError("ERROR: lower threshold cannot be greater than upper threshold.");
 196       return;
 197     }
 198   }
 199   // Threshold to get air
 200   typedef itk::BinaryThresholdImageFilter<ImageType, InternalImageType> InputBinarizeFilterType;
 201   typename InputBinarizeFilterType::Pointer binarizeFilter=InputBinarizeFilterType::New();
 202   binarizeFilter->SetInput(working_input);
 203   if (m_UseLowerThreshold) binarizeFilter->SetLowerThreshold(m_LowerThreshold);
 204   binarizeFilter->SetUpperThreshold(m_UpperThreshold);
 205   binarizeFilter ->SetInsideValue(this->GetForegroundValue());
 206   binarizeFilter ->SetOutsideValue(this->GetBackgroundValue());
 207   binarizeFilter->Update();
 208   working_image = binarizeFilter->GetOutput();
 209
 210   // Labelize and keep right labels
 211   working_image = Labelize<InternalImageType>(working_image, GetBackgroundValue(), true, GetMinimalComponentSize());
 212
 213   working_image = RemoveLabels<InternalImageType>
 214     (working_image, GetBackgroundValue(), GetLabelizeParameters1()->GetLabelsToRemove());
 215
 216   typename InternalImageType::Pointer air = KeepLabels<InternalImageType>
 217     (working_image,
 218      GetBackgroundValue(),
 219      GetForegroundValue(),
 220      GetLabelizeParameters1()->GetFirstKeep(),
 221      GetLabelizeParameters1()->GetLastKeep(),
 222      GetLabelizeParameters1()->GetUseLastKeep());
 223
 224   // Set Air to BG
 225   working_input = SetBackground<ImageType, InternalImageType>
 226     (working_input, air, this->GetForegroundValue(), this->GetBackgroundValue());
 227
 228   // End
 229   StopCurrentStep<ImageType>(working_input);
 230
 231   //--------------------------------------------------------------------
 232   //--------------------------------------------------------------------
 233   StartNewStepOrStop("Search for the trachea");
 234   SearchForTrachea();
 235
 236   //--------------------------------------------------------------------
 237   //--------------------------------------------------------------------
 238   StartNewStepOrStop("Extract the lung with Otsu filter");
 239   // Automated Otsu thresholding and relabeling
 240   typedef itk::OtsuThresholdImageFilter<ImageType,InternalImageType> OtsuThresholdImageFilterType;
 241   typename OtsuThresholdImageFilterType::Pointer otsuFilter=OtsuThresholdImageFilterType::New();
 242   otsuFilter->SetInput(working_input);
 243   otsuFilter->SetNumberOfHistogramBins(GetNumberOfHistogramBins());
 244   otsuFilter->SetInsideValue(this->GetForegroundValue());
 245   otsuFilter->SetOutsideValue(this->GetBackgroundValue());
 246   otsuFilter->Update();
 247   working_image =  otsuFilter->GetOutput();
 248
 249   // Set output
 250   StopCurrentStep<InternalImageType>(working_image);
 251
 252   //--------------------------------------------------------------------
 253   //--------------------------------------------------------------------
 254   StartNewStepOrStop("Select labels");
 255   // Keep right labels
 256   working_image = LabelizeAndSelectLabels<InternalImageType>
 257     (working_image,
 258      GetBackgroundValue(),
 259      GetForegroundValue(),
 260      false,
 261      GetMinimalComponentSize(),
 262      GetLabelizeParameters2());
 263
 264   // Set output
 265   StopCurrentStep<InternalImageType>(working_image);
 266
 267   //--------------------------------------------------------------------
 268   //--------------------------------------------------------------------
 269   if (m_Seeds.size() != 0) { // if ==0 ->no trachea found
 270     StartNewStepOrStop("Remove the trachea");
 271     // Set the trachea
 272     working_image = SetBackground<InternalImageType, InternalImageType>
 273       (working_image, trachea_tmp, 1, -1);
 274
 275     // Dilate the trachea
 276     static const unsigned int Dim = ImageType::ImageDimension;
 277     typedef itk::BinaryBallStructuringElement<InternalPixelType, Dim> KernelType;
 278     KernelType structuringElement;
 279     structuringElement.SetRadius(GetRadiusForTrachea());
 280     structuringElement.CreateStructuringElement();
 281     typedef clitk::ConditionalBinaryDilateImageFilter<InternalImageType, InternalImageType, KernelType> ConditionalBinaryDilateImageFilterType;
 282     typename ConditionalBinaryDilateImageFilterType::Pointer dilateFilter = ConditionalBinaryDilateImageFilterType::New();
 283     dilateFilter->SetBoundaryToForeground(false);
 284     dilateFilter->SetKernel(structuringElement);
 285     dilateFilter->SetBackgroundValue (1);
 286     dilateFilter->SetForegroundValue (-1);
 287     dilateFilter->SetInput (working_image);
 288     dilateFilter->Update();
 289     working_image = dilateFilter->GetOutput();
 290
 291     // Set trachea with dilatation
 292     trachea_tmp = SetBackground<InternalImageType, InternalImageType>
 293       (trachea_tmp, working_image, -1, this->GetForegroundValue());
 294
 295     // Remove the trachea
 296     working_image = SetBackground<InternalImageType, InternalImageType>
 297       (working_image, working_image, -1, this->GetBackgroundValue());
 298
 299     // Label
 300     working_image = LabelizeAndSelectLabels<InternalImageType>
 301       (working_image,
 302        GetBackgroundValue(),
 303        GetForegroundValue(),
 304        false,
 305        GetMinimalComponentSize(),
 306        GetLabelizeParameters3());
 307
 308     // Set output
 309     StopCurrentStep<InternalImageType>(working_image);
 310   }
 311
 312   //--------------------------------------------------------------------
 313   //--------------------------------------------------------------------
 314   typedef clitk::AutoCropFilter<InternalImageType> CropFilterType;
 315   typename CropFilterType::Pointer cropFilter = CropFilterType::New();
 316   if (m_Seeds.size() != 0) { // if ==0 ->no trachea found
 317     StartNewStepOrStop("Croping trachea");
 318     cropFilter->SetInput(trachea_tmp);
 319     cropFilter->Update(); // Needed
 320     typedef itk::CastImageFilter<InternalImageType, MaskImageType> CastImageFilterType;
 321     typename CastImageFilterType::Pointer caster= CastImageFilterType::New();
 322     caster->SetInput(cropFilter->GetOutput());
 323     caster->Update();
 324     trachea = caster->GetOutput();
 325     StopCurrentStep<MaskImageType>(trachea);
 326   }
 327
 328   //--------------------------------------------------------------------
 329   //--------------------------------------------------------------------
 330   StartNewStepOrStop("Croping lung");
 331   typename CropFilterType::Pointer cropFilter2 = CropFilterType::New(); // Needed to reset pipeline
 332   cropFilter2->SetInput(working_image);
 333   cropFilter2->Update();
 334   working_image = cropFilter2->GetOutput();
 335   StopCurrentStep<InternalImageType>(working_image);
 336
 337   //--------------------------------------------------------------------
 338   //--------------------------------------------------------------------
 339   StartNewStepOrStop("Separate Left/Right lungs");
 340   // Initial label
 341   working_image = Labelize<InternalImageType>(working_image,
 342                                               GetBackgroundValue(),
 343                                               false,
 344                                               GetMinimalComponentSize());
 345
 346   // Count the labels
 347   typedef itk::StatisticsImageFilter<InternalImageType> StatisticsImageFilterType;
 348   typename StatisticsImageFilterType::Pointer statisticsImageFilter=StatisticsImageFilterType::New();
 349   statisticsImageFilter->SetInput(working_image);
 350   statisticsImageFilter->Update();
 351   unsigned int initialNumberOfLabels = statisticsImageFilter->GetMaximum();
 352   working_image = statisticsImageFilter->GetOutput();
 353
 354   // Decompose the first label
 355   static const unsigned int Dim = ImageType::ImageDimension;
 356   if (initialNumberOfLabels<2) {
 357     // Structuring element radius
 358     typename ImageType::SizeType radius;
 359     for (unsigned int i=0;i<Dim;i++) radius[i]=1;
 360     typedef clitk::DecomposeAndReconstructImageFilter<InternalImageType,InternalImageType> DecomposeAndReconstructFilterType;
 361     typename DecomposeAndReconstructFilterType::Pointer decomposeAndReconstructFilter=DecomposeAndReconstructFilterType::New();
 362     decomposeAndReconstructFilter->SetInput(working_image);
 363     decomposeAndReconstructFilter->SetVerbose(false);
 364     decomposeAndReconstructFilter->SetRadius(radius);
 365     decomposeAndReconstructFilter->SetMaximumNumberOfLabels(2);
 366     decomposeAndReconstructFilter->SetMinimumObjectSize(this->GetMinimalComponentSize());
 367     decomposeAndReconstructFilter->SetMinimumNumberOfIterations(1);
 368     decomposeAndReconstructFilter->SetBackgroundValue(this->GetBackgroundValue());
 369     decomposeAndReconstructFilter->SetForegroundValue(this->GetForegroundValue());
 370     decomposeAndReconstructFilter->SetFullyConnected(true);
 371     decomposeAndReconstructFilter->SetNumberOfNewLabels(1);
 372     decomposeAndReconstructFilter->Update();
 373     working_image = decomposeAndReconstructFilter->GetOutput();
 374   }
 375
 376   // Retain labels ('1' is largset lung, so right. '2' is left)
 377   typedef itk::ThresholdImageFilter<InternalImageType> ThresholdImageFilterType;
 378   typename ThresholdImageFilterType::Pointer thresholdFilter = ThresholdImageFilterType::New();
 379   thresholdFilter->SetInput(working_image);
 380   thresholdFilter->ThresholdAbove(2);
 381   thresholdFilter->SetOutsideValue(this->GetBackgroundValue());
 382   thresholdFilter->Update();
 383   working_image = thresholdFilter->GetOutput();
 384   StopCurrentStep<InternalImageType> (working_image);
 385
 386   // Final Cast
 387   StartNewStepOrStop("Final cast");
 388   typedef itk::CastImageFilter<InternalImageType, MaskImageType> CastImageFilterType;
 389   typename CastImageFilterType::Pointer caster= CastImageFilterType::New();
 390   caster->SetInput(working_image);
 391   caster->Update();
 392   output = caster->GetOutput();
 393
 394   // Update output info
 395   this->GetOutput(0)->SetRegions(output->GetLargestPossibleRegion());
 396
 397   // Try to extract bifurcation in the trachea (bronchi)
 398   // STILL EXPERIMENTAL
 399   if (m_Seeds.size() != 0) { // if ==0 ->no trachea found
 400
 401     if (GetFindBronchialBifurcations()) {
 402       StartNewStepOrStop("Find bronchial bifurcations");
 403       // Step 1 : extract skeleton
 404       // Define the thinning filter
 405       typedef itk::BinaryThinningImageFilter3D<MaskImageType, MaskImageType> ThinningFilterType;
 406       typename ThinningFilterType::Pointer thinningFilter = ThinningFilterType::New();
 407       thinningFilter->SetInput(trachea);
 408       thinningFilter->Update();
 409       typename MaskImageType::Pointer skeleton = thinningFilter->GetOutput();
 410       writeImage<MaskImageType>(skeleton, "skeleton.mhd");
 411
 412       // Step 2 : tracking
 413       DD("tracking");
 414
 415       // Step 2.1 : find first point for tracking
 416       typedef itk::ImageRegionConstIteratorWithIndex<MaskImageType> IteratorType;
 417       IteratorType it(skeleton, skeleton->GetLargestPossibleRegion());
 418       it.GoToReverseBegin();
 419       while ((!it.IsAtEnd()) && (it.Get() == GetBackgroundValue())) {
 420         --it;
 421       }
 422       if (it.IsAtEnd()) {
 423         this->SetLastError("ERROR: first point in the skeleton not found ! Abort");
 424         return;
 425       }
 426       DD(skeleton->GetLargestPossibleRegion().GetIndex());
 427       typename MaskImageType::IndexType index = it.GetIndex();
 428       DD(index);
 429
 430       // Step 2.2 : initialize neighborhooditerator
 431       typedef itk::NeighborhoodIterator<MaskImageType> NeighborhoodIteratorType;
 432       typename NeighborhoodIteratorType::SizeType radius;
 433       radius.Fill(1);
 434       NeighborhoodIteratorType nit(radius, skeleton, skeleton->GetLargestPossibleRegion());
 435       DD(nit.GetSize());
 436       DD(nit.Size());
 437
 438       // Find first label number (must be different from BG and FG)
 439       typename MaskImageType::PixelType label = GetForegroundValue()+1;
 440       while ((label == GetBackgroundValue()) || (label == GetForegroundValue())) { label++; }
 441       DD(label);
 442
 443       // Track from the first point
 444       std::vector<BifurcationType> listOfBifurcations;
 445       TrackFromThisIndex(listOfBifurcations, skeleton, index, label);
 446       DD("end track");
 447       DD(listOfBifurcations.size());
 448       writeImage<MaskImageType>(skeleton, "skeleton2.mhd");
 449
 450       for(unsigned int i=0; i<listOfBifurcations.size(); i++) {
 451         typename MaskImageType::PointType p;
 452         skeleton->TransformIndexToPhysicalPoint(listOfBifurcations[i].index, p);
 453         DD(p);
 454       }
 455
 456     }
 457   }
 458 }
 459 //--------------------------------------------------------------------
 460
 461
 462 //--------------------------------------------------------------------
 463 template <class ImageType, class MaskImageType>
 464 void
 465 clitk::ExtractLungFilter<ImageType, MaskImageType>::
 466 GenerateData()
 467 {
 468   // Do not put some "startnewstep" here, because the object if
 469   // modified and the filter's pipeline it do two times. But it is
 470   // required to quit if MustStop was set before.
 471   if (GetMustStop()) return;
 472
 473   // If everything goes well, set the output
 474   this->GraftOutput(output); // not SetNthOutput
 475 }
 476 //--------------------------------------------------------------------
 477
 478
 479 //--------------------------------------------------------------------
 480 template <class ImageType, class MaskImageType>
 481 void
 482 clitk::ExtractLungFilter<ImageType, MaskImageType>::
 483 TrackFromThisIndex(std::vector<BifurcationType> & listOfBifurcations,
 484                    MaskImagePointer skeleton,
 485                    MaskImageIndexType index,
 486                    MaskImagePixelType label)
 487 {
 488   DD("TrackFromThisIndex");
 489   DD(index);
 490   DD((int)label);
 491   // Create NeighborhoodIterator
 492   typedef itk::NeighborhoodIterator<MaskImageType> NeighborhoodIteratorType;
 493   typename NeighborhoodIteratorType::SizeType radius;
 494   radius.Fill(1);
 495   NeighborhoodIteratorType nit(radius, skeleton, skeleton->GetLargestPossibleRegion());
 496
 497   // Track
 498   std::vector<typename NeighborhoodIteratorType::IndexType> listOfTrackedPoint;
 499   bool stop = false;
 500   while (!stop) {
 501     nit.SetLocation(index);
 502     // DD((int)nit.GetCenterPixel());
 503     nit.SetCenterPixel(label);
 504     listOfTrackedPoint.clear();
 505     for(unsigned int i=0; i<nit.Size(); i++) {
 506       if (i != nit.GetCenterNeighborhoodIndex ()) { // Do not observe the current point
 507         //          DD(nit.GetIndex(i));
 508         if (nit.GetPixel(i) == GetForegroundValue()) { // if this is foreground, we continue the tracking
 509           // DD(nit.GetIndex(i));
 510           listOfTrackedPoint.push_back(nit.GetIndex(i));
 511         }
 512       }
 513     }
 514     // DD(listOfTrackedPoint.size());
 515     if (listOfTrackedPoint.size() == 1) {
 516       index = listOfTrackedPoint[0];
 517     }
 518     else {
 519       if (listOfTrackedPoint.size() == 2) {
 520         BifurcationType bif(index, label, label+1, label+2);
 521         listOfBifurcations.push_back(bif);
 522         TrackFromThisIndex(listOfBifurcations, skeleton, listOfTrackedPoint[0], label+1);
 523         TrackFromThisIndex(listOfBifurcations, skeleton, listOfTrackedPoint[1], label+2);
 524       }
 525       else {
 526         if (listOfTrackedPoint.size() > 2) {
 527           std::cerr << "too much bifurcation points ... ?" << std::endl;
 528           exit(0);
 529         }
 530         // Else this it the end of the tracking
 531       }
 532       stop = true;
 533     }
 534   }
 535 }
 536 //--------------------------------------------------------------------
 537
 538
 539 //--------------------------------------------------------------------
 540 template <class ImageType, class MaskImageType>
 541 bool
 542 clitk::ExtractLungFilter<ImageType, MaskImageType>::
 543 SearchForTracheaSeed(int skip)
 544 {
 545   if (m_Seeds.size() == 0) { // try to find seed (if not zero, it is given by user)
 546     // Restart the search until a seed is found, skipping more and more slices
 547     bool stop = false;
 548     while (!stop) {
 549       // Search seed (parameters = UpperThresholdForTrachea)
 550       static const unsigned int Dim = ImageType::ImageDimension;
 551       typename InternalImageType::RegionType sliceRegion = working_input->GetLargestPossibleRegion();
 552       typename InternalImageType::SizeType sliceRegionSize = sliceRegion.GetSize();
 553       typename InternalImageType::IndexType sliceRegionIndex = sliceRegion.GetIndex();
 554       sliceRegionIndex[Dim-1]=sliceRegionSize[Dim-1]-skip-5;
 555       sliceRegionSize[Dim-1]=5;
 556       sliceRegion.SetSize(sliceRegionSize);
 557       sliceRegion.SetIndex(sliceRegionIndex);
 558       typedef  itk::ImageRegionConstIterator<ImageType> IteratorType;
 559       IteratorType it(working_input, sliceRegion);
 560       it.GoToBegin();
 561       while (!it.IsAtEnd()) {
 562         if(it.Get() < GetUpperThresholdForTrachea() ) {
 563           AddSeed(it.GetIndex());
 564           // DD(it.GetIndex());
 565         }
 566         ++it;
 567       }
 568
 569       // if we do not found : restart
 570       stop = (m_Seeds.size() != 0);
 571       if (!stop) {
 572         if (GetVerboseStep()) {
 573           std::cout << "\t No seed found this time. I skip some slices and restart." << std::endl;
 574         }
 575         if (skip > 0.5 * working_input->GetLargestPossibleRegion().GetSize()[2]) {
 576           // we want to skip more than a half of the image, it is probably a bug
 577           std::cerr << "Number of slices to skip to find trachea to high = " << skip << std::endl;
 578           stop = true;
 579         }
 580         skip += 5;
 581       }
 582     }
 583   }
 584   return (m_Seeds.size() != 0);
 585 }
 586 //--------------------------------------------------------------------
 587
 588
 589 //--------------------------------------------------------------------
 590 template <class ImageType, class MaskImageType>
 591 void
 592 clitk::ExtractLungFilter<ImageType, MaskImageType>::
 593 TracheaRegionGrowing()
 594 {
 595   // Explosion controlled region growing
 596   typedef clitk::ExplosionControlledThresholdConnectedImageFilter<ImageType, InternalImageType> ImageFilterType;
 597   typename ImageFilterType::Pointer f= ImageFilterType::New();
 598   f->SetInput(working_input);
 599   f->SetVerbose(false);
 600   f->SetLower(-2000);
 601   f->SetUpper(GetUpperThresholdForTrachea());
 602   f->SetMinimumLowerThreshold(-2000);
 603   f->SetMaximumUpperThreshold(0);
 604   f->SetAdaptLowerBorder(false);
 605   f->SetAdaptUpperBorder(true);
 606   f->SetMinimumSize(5000);
 607   f->SetReplaceValue(1);
 608   f->SetMultiplier(GetMultiplierForTrachea());
 609   f->SetThresholdStepSize(GetThresholdStepSizeForTrachea());
 610   f->SetMinimumThresholdStepSize(1);
 611   for(unsigned int i=0; i<m_Seeds.size();i++) {
 612     f->AddSeed(m_Seeds[i]);
 613     // DD(m_Seeds[i]);
 614   }
 615   f->Update();
 616
 617   // take first (main) connected component
 618   writeImage<InternalImageType>(f->GetOutput(), "t1.mhd");
 619   trachea_tmp = Labelize<InternalImageType>(f->GetOutput(),
 620                                             GetBackgroundValue(),
 621                                             true,
 622                                             GetMinimalComponentSize());
 623   writeImage<InternalImageType>(trachea_tmp, "t2.mhd");
 624   trachea_tmp = KeepLabels<InternalImageType>(trachea_tmp,
 625                                               GetBackgroundValue(),
 626                                               GetForegroundValue(),
 627                                               1, 1, false);
 628   writeImage<InternalImageType>(trachea_tmp, "t3.mhd");
 629 }
 630 //--------------------------------------------------------------------
 631
 632
 633 //--------------------------------------------------------------------
 634 template <class ImageType, class MaskImageType>
 635 double
 636 clitk::ExtractLungFilter<ImageType, MaskImageType>::
 637 ComputeTracheaVolume()
 638 {
 639   typedef itk::ImageRegionConstIterator<InternalImageType> IteratorType;
 640   IteratorType iter(trachea_tmp, trachea_tmp->GetLargestPossibleRegion());
 641   iter.GoToBegin();
 642   double volume = 0.0;
 643   while (!iter.IsAtEnd()) {
 644     if (iter.Get() == this->GetForegroundValue()) volume++;
 645     ++iter;
 646   }
 647
 648   double voxelsize = trachea_tmp->GetSpacing()[0]*trachea_tmp->GetSpacing()[1]*trachea_tmp->GetSpacing()[2];
 649   return volume*voxelsize;
 650 }
 651 //--------------------------------------------------------------------
 652
 653
 654 //--------------------------------------------------------------------
 655 template <class ImageType, class MaskImageType>
 656 void
 657 clitk::ExtractLungFilter<ImageType, MaskImageType>::
 658 SearchForTrachea()
 659 {
 660   // Search for seed among n slices, skip some slices before starting
 661   // if not found -> skip more and restart
 662
 663   // when seed found : perform region growing
 664   // compute trachea volume
 665   // if volume not plausible  -> skip more slices and restart
 666
 667   bool stop = false;
 668   double volume = 0.0;
 669   int skip = GetNumberOfSlicesToSkipBeforeSearchingSeed();
 670   while (!stop) {
 671     stop = SearchForTracheaSeed(skip);
 672     if (stop) {
 673       TracheaRegionGrowing();
 674       volume = ComputeTracheaVolume(); // assume mm3
 675       if ((volume > 10000) && (volume < 55000 )) { // it is ok
 676         // Typical volume 22.59 cm 3 (± 7.69 cm 3 ) [ Leader 2004 ]
 677         if (GetVerboseStep()) {
 678           std::cout << "\t Found trachea with volume " << volume << " cc." << std::endl;
 679         }
 680         stop = true;
 681       }
 682       else {
 683         if (GetVerboseStep()) {
 684           std::cout << "\t The volume of the trachea (" << volume
 685                     << " cc) seems not correct. I skip some slices (" << skip << ") and restart to find seeds."
 686                     << std::endl;
 687         }
 688         skip += 5;
 689         stop = false;
 690         // empty the list of seed
 691         m_Seeds.clear();
 692       }
 693     }
 694   }
 695
 696   if (volume != 0.0) {
 697     // Set output
 698     StopCurrentStep<InternalImageType>(trachea_tmp);
 699   }
 700   else { // Trachea not found
 701     this->SetWarning("* WARNING * No seed found for trachea.");
 702     StopCurrentStep();
 703   }
 704 }
 705 //--------------------------------------------------------------------
 706
 707 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX