segmentation/clitkExtractMediastinumFilter.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 CLITKEXTRACTMEDIASTINUMFILTER_TXX
  20 #define CLITKEXTRACTMEDIASTINUMFILTER_TXX
  21
  22 // clitk
  23 #include "clitkCommon.h"
  24 #include "clitkExtractMediastinumFilter.h"
  25 #include "clitkAddRelativePositionConstraintToLabelImageFilter.h"
  26 #include "clitkSliceBySliceRelativePositionFilter.h"
  27 #include "clitkSegmentationUtils.h"
  28 #include "clitkExtractAirwaysTreeInfoFilter.h"
  29 #include "clitkCropLikeImageFilter.h"
  30
  31 // std
  32 #include <deque>
  33
  34 // itk
  35 #include "itkStatisticsLabelMapFilter.h"
  36 #include "itkLabelImageToStatisticsLabelMapFilter.h"
  37 #include "itkRegionOfInterestImageFilter.h"
  38 #include "itkBinaryThresholdImageFilter.h"
  39 #include "itkScalarImageKmeansImageFilter.h"
  40
  41 // itk ENST
  42 #include "RelativePositionPropImageFilter.h"
  43
  44 //--------------------------------------------------------------------
  45 template <class ImageType>
  46 clitk::ExtractMediastinumFilter<ImageType>::
  47 ExtractMediastinumFilter():
  48   clitk::FilterBase(),
  49   clitk::FilterWithAnatomicalFeatureDatabaseManagement(),
  50   itk::ImageToImageFilter<ImageType, MaskImageType>()
  51 {
  52   this->SetNumberOfRequiredInputs(1);
  53
  54   SetBackgroundValuePatient(0);
  55   SetBackgroundValueLung(0);
  56   SetBackgroundValueBones(0);
  57   SetForegroundValueLeftLung(1);
  58   SetForegroundValueRightLung(2);
  59   SetBackgroundValue(0);
  60   SetForegroundValue(1);
  61
  62   SetIntermediateSpacing(6);
  63   SetFuzzyThreshold1(0.5);
  64   SetFuzzyThreshold2(0.6);
  65   SetFuzzyThreshold3(0.05);
  66
  67   SetOutputMediastinumFilename("mediastinum.mhd");
  68
  69   UseBonesOff();
  70 }
  71 //--------------------------------------------------------------------
  72
  73
  74 //--------------------------------------------------------------------
  75 template <class ImageType>
  76 void
  77 clitk::ExtractMediastinumFilter<ImageType>::
  78 SetInputPatientLabelImage(const MaskImageType * image, MaskImagePixelType bg)
  79 {
  80   this->SetNthInput(0, const_cast<MaskImageType *>(image));
  81   m_BackgroundValuePatient = bg;
  82 }
  83 //--------------------------------------------------------------------
  84
  85
  86 //--------------------------------------------------------------------
  87 template <class ImageType>
  88 void
  89 clitk::ExtractMediastinumFilter<ImageType>::
  90 SetInputLungLabelImage(const MaskImageType * image, MaskImagePixelType bg,
  91                        MaskImagePixelType fgLeft, MaskImagePixelType fgRight)
  92 {
  93   this->SetNthInput(1, const_cast<MaskImageType *>(image));
  94   m_BackgroundValueLung = bg;
  95   m_ForegroundValueLeftLung = fgLeft;
  96   m_ForegroundValueRightLung = fgRight;
  97 }
  98 //--------------------------------------------------------------------
  99
 100
 101 //--------------------------------------------------------------------
 102 template <class ImageType>
 103 void
 104 clitk::ExtractMediastinumFilter<ImageType>::
 105 SetInputBonesLabelImage(const MaskImageType * image, MaskImagePixelType bg)
 106 {
 107   this->SetNthInput(2, const_cast<MaskImageType *>(image));
 108   m_BackgroundValueBones = bg;
 109 }
 110 //--------------------------------------------------------------------
 111
 112
 113 //--------------------------------------------------------------------
 114 template <class ImageType>
 115 void
 116 clitk::ExtractMediastinumFilter<ImageType>::
 117 SetInputTracheaLabelImage(const MaskImageType * image, MaskImagePixelType bg)
 118 {
 119   this->SetNthInput(3, const_cast<MaskImageType *>(image));
 120   m_BackgroundValueTrachea = bg;
 121 }
 122 //--------------------------------------------------------------------
 123
 124
 125 //--------------------------------------------------------------------
 126 template <class ImageType>
 127 void
 128 clitk::ExtractMediastinumFilter<ImageType>::
 129 GenerateInputRequestedRegion()
 130 {
 131   // DD("GenerateInputRequestedRegion");
 132   // Do not call default
 133   // Superclass::GenerateInputRequestedRegion();
 134   // DD("End GenerateInputRequestedRegion");
 135 }
 136
 137 //--------------------------------------------------------------------
 138
 139
 140 //--------------------------------------------------------------------
 141 template <class ImageType>
 142 void
 143 clitk::ExtractMediastinumFilter<ImageType>::
 144 SetInput(const ImageType * image)
 145 {
 146   this->SetNthInput(0, const_cast<ImageType *>(image));
 147 }
 148 //--------------------------------------------------------------------
 149
 150
 151 //--------------------------------------------------------------------
 152 template <class ImageType>
 153 void
 154 clitk::ExtractMediastinumFilter<ImageType>::
 155 GenerateOutputInformation() {
 156   // Do not call default
 157   // Superclass::GenerateOutputInformation();
 158
 159   //--------------------------------------------------------------------
 160   // Get input pointers
 161   clitk::PrintMemory(GetVerboseMemoryFlag(), "Initial memory"); // OK
 162   LoadAFDB();
 163   ImageConstPointer input = dynamic_cast<const ImageType*>(itk::ProcessObject::GetInput(0));
 164   MaskImagePointer patient = GetAFDB()->template GetImage <MaskImageType>("Patient");
 165   MaskImagePointer lung = GetAFDB()->template GetImage <MaskImageType>("Lungs");
 166   MaskImagePointer bones;
 167   if (GetUseBones()) {
 168     bones = GetAFDB()->template GetImage <MaskImageType>("Bones");
 169   }
 170   MaskImagePointer trachea = GetAFDB()->template GetImage <MaskImageType>("Trachea");
 171
 172   clitk::PrintMemory(GetVerboseMemoryFlag(), "After read patient, lung");
 173
 174   //--------------------------------------------------------------------
 175   // Step 1: Crop support (patient) to lung extend in RL
 176   StartNewStep("Crop support like lungs along LR");
 177   typedef clitk::CropLikeImageFilter<MaskImageType> CropFilterType;
 178   typename CropFilterType::Pointer cropFilter = CropFilterType::New();
 179   cropFilter->SetInput(patient);
 180   cropFilter->SetCropLikeImage(lung, 0);// Indicate that we only crop in X (Left-Right) axe
 181   cropFilter->Update();
 182   output = cropFilter->GetOutput();
 183   this->template StopCurrentStep<MaskImageType>(output);
 184
 185   //--------------------------------------------------------------------
 186   // Step 2: Crop support (previous) to bones extend in AP
 187   if (GetUseBones()) {
 188     StartNewStep("Crop support like bones along AP");
 189     cropFilter = CropFilterType::New();
 190     cropFilter->SetInput(output);
 191     cropFilter->SetCropLikeImage(bones, 1);// Indicate that we only crop in Y (Ant-Post) axe
 192     cropFilter->Update();
 193     output = cropFilter->GetOutput();
 194     this->template StopCurrentStep<MaskImageType>(output);
 195   }
 196
 197   //--------------------------------------------------------------------
 198   // Step 3: patient minus lungs, minus bones, minus trachea
 199   StartNewStep("Patient contours minus lungs, trachea [and bones]");
 200   typedef clitk::BooleanOperatorLabelImageFilter<MaskImageType> BoolFilterType;
 201   typename BoolFilterType::Pointer boolFilter = BoolFilterType::New();
 202   boolFilter->InPlaceOn();
 203   boolFilter->SetInput1(output);
 204   boolFilter->SetInput2(lung);
 205   boolFilter->SetOperationType(BoolFilterType::AndNot);
 206   boolFilter->Update();
 207   if (GetUseBones()) {
 208     boolFilter->SetInput1(boolFilter->GetOutput());
 209     boolFilter->SetInput2(bones);
 210     boolFilter->SetOperationType(BoolFilterType::AndNot);
 211     boolFilter->Update();
 212   }
 213   boolFilter->SetInput1(boolFilter->GetOutput());
 214   boolFilter->SetInput2(trachea);
 215   boolFilter->SetOperationType(BoolFilterType::AndNot);
 216   boolFilter->Update();
 217   output = boolFilter->GetOutput();
 218
 219   // Auto crop to gain support area
 220   output = clitk::AutoCrop<MaskImageType>(output, GetBackgroundValue());
 221   this->template StopCurrentStep<MaskImageType>(output);
 222
 223   //--------------------------------------------------------------------
 224   // Step 4: LR limits from lung (need separate lung ?)
 225   // Get separate lung images to get only the right and left lung
 226   // (because RelativePositionPropImageFilter only consider fg=1);
 227   // (label must be '1' because right is greater than left).  (WE DO
 228   // NOT NEED TO SEPARATE ? )
 229   StartNewStep("Left/Right limits with lungs");
 230
 231   // The following cannot be "inplace" because mask is the same than input ...
 232   MaskImagePointer right_lung =
 233     clitk::SetBackground<MaskImageType, MaskImageType>(lung, lung, 2, 0, false);
 234   MaskImagePointer left_lung =
 235     clitk::SetBackground<MaskImageType, MaskImageType>(lung, lung, 1, 0, false);
 236   right_lung = clitk::ResizeImageLike<MaskImageType>(right_lung, output, GetBackgroundValue());
 237   left_lung = clitk::ResizeImageLike<MaskImageType>(left_lung, output, GetBackgroundValue());
 238   // writeImage<MaskImageType>(right_lung, "right.mhd");
 239   // writeImage<MaskImageType>(left_lung, "left.mhd");
 240
 241   typedef clitk::AddRelativePositionConstraintToLabelImageFilter<MaskImageType> RelPosFilterType;
 242   typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New();
 243   relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId());
 244   relPosFilter->VerboseStepFlagOff();
 245   relPosFilter->WriteStepFlagOff();
 246   relPosFilter->SetInput(output);
 247   relPosFilter->SetInputObject(left_lung);
 248   //  relPosFilter->SetInputObject(lung);
 249   relPosFilter->AddOrientationType(RelPosFilterType::AtRightTo); // warning left lung is at right ;)
 250   relPosFilter->SetIntermediateSpacing(GetIntermediateSpacing());
 251   relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold1());
 252   relPosFilter->Update();
 253   output = relPosFilter->GetOutput();
 254   //writeImage<MaskImageType>(right_lung, "step4-left.mhd");
 255
 256   relPosFilter = RelPosFilterType::New();
 257   relPosFilter->SetInput(output);
 258   relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId());
 259   relPosFilter->VerboseStepFlagOff();
 260   relPosFilter->WriteStepFlagOff();
 261   relPosFilter->SetInput(output);
 262   relPosFilter->SetInputObject(right_lung);
 263   //relPosFilter->SetInputObject(lung);
 264   relPosFilter->AddOrientationType(RelPosFilterType::AtLeftTo);
 265   relPosFilter->SetIntermediateSpacing(GetIntermediateSpacing());
 266   relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold1());
 267   relPosFilter->Update();
 268   output = relPosFilter->GetOutput();
 269   this->template StopCurrentStep<MaskImageType>(output);
 270
 271   //--------------------------------------------------------------------
 272   // Step 5: AP limits from bones
 273   // Separate the bones in the ant-post middle
 274   MaskImageConstPointer bones_ant;
 275   MaskImageConstPointer bones_post;
 276   MaskImagePointType middle_AntPost__position;
 277   if (GetUseBones()) {
 278     StartNewStep("Ant/Post limits with bones");
 279     middle_AntPost__position[0] = middle_AntPost__position[2] = 0;
 280     middle_AntPost__position[1] = bones->GetOrigin()[1]+(bones->GetLargestPossibleRegion().GetSize()[1]*bones->GetSpacing()[1])/2.0;
 281     MaskImageIndexType index_bones_middle;
 282     bones->TransformPhysicalPointToIndex(middle_AntPost__position, index_bones_middle);
 283
 284     // Split bone image first into two parts (ant and post), and crop
 285     // lateraly to get vertebral
 286     typedef itk::RegionOfInterestImageFilter<MaskImageType, MaskImageType> ROIFilterType;
 287     //  typedef itk::ExtractImageFilter<MaskImageType, MaskImageType> ROIFilterType;
 288     typename ROIFilterType::Pointer roiFilter = ROIFilterType::New();
 289     MaskImageRegionType region = bones->GetLargestPossibleRegion();
 290     MaskImageSizeType size = region.GetSize();
 291     MaskImageIndexType index = region.GetIndex();
 292     // ANT part
 293     // crop LR to keep 1/4 center part
 294     index[0] = size[0]/4+size[0]/8;
 295     size[0] = size[0]/4;
 296     // crop AP to keep first (ant) part
 297     size[1] =  index_bones_middle[1]; //size[1]/2.0;
 298     region.SetSize(size);
 299     region.SetIndex(index);
 300     roiFilter->SetInput(bones);
 301     roiFilter->SetRegionOfInterest(region);
 302     roiFilter->ReleaseDataFlagOff();
 303     roiFilter->Update();
 304     bones_ant = roiFilter->GetOutput();
 305     //    writeImage<MaskImageType>(bones_ant, "b_ant.mhd");
 306     // POST part
 307     roiFilter = ROIFilterType::New();
 308     index[1] = bones->GetLargestPossibleRegion().GetIndex()[1] + size[1]-1;
 309     size[1] =  bones->GetLargestPossibleRegion().GetSize()[1] - size[1];
 310     region.SetIndex(index);
 311     region.SetSize(size);
 312     roiFilter->SetInput(bones);
 313     roiFilter->SetRegionOfInterest(region);
 314     roiFilter->ReleaseDataFlagOff();
 315     roiFilter->Update();
 316     bones_post = roiFilter->GetOutput();
 317     //    writeImage<MaskImageType>(bones_post, "b_post.mhd");
 318
 319     // Go !
 320     relPosFilter->SetCurrentStepNumber(0);
 321     relPosFilter->ResetPipeline();// = RelPosFilterType::New();
 322     relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId());
 323     relPosFilter->VerboseStepFlagOff();
 324     relPosFilter->WriteStepFlagOff();
 325     relPosFilter->SetInput(output);
 326     relPosFilter->SetInputObject(bones_post);
 327     relPosFilter->AddOrientationType(RelPosFilterType::AntTo);
 328     relPosFilter->SetIntermediateSpacing(GetIntermediateSpacing());
 329     relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold2());
 330     relPosFilter->Update();
 331     output = relPosFilter->GetOutput();
 332     //    writeImage<MaskImageType>(output, "post.mhd");
 333
 334     relPosFilter->SetInput(relPosFilter->GetOutput());
 335     relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId());
 336     relPosFilter->VerboseStepFlagOff();
 337     relPosFilter->WriteStepFlagOff();
 338     relPosFilter->SetInput(output);
 339     relPosFilter->SetInputObject(bones_ant);
 340     relPosFilter->AddOrientationType(RelPosFilterType::PostTo);
 341     relPosFilter->SetIntermediateSpacing(GetIntermediateSpacing());
 342     relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold2());
 343     relPosFilter->Update();
 344     output = relPosFilter->GetOutput();
 345     this->template StopCurrentStep<MaskImageType>(output);
 346   }
 347
 348   //--------------------------------------------------------------------
 349   // Step 6: Get CCL
 350   StartNewStep("Keep main connected component");
 351   output = clitk::Labelize<MaskImageType>(output, GetBackgroundValue(), false, 500);
 352   // output = RemoveLabels<MaskImageType>(output, BG, param->GetLabelsToRemove());
 353   output = clitk::KeepLabels<MaskImageType>(output, GetBackgroundValue(),
 354                                             GetForegroundValue(), 1, 1, 0);
 355   this->template StopCurrentStep<MaskImageType>(output);
 356
 357
 358   //--------------------------------------------------------------------
 359   // Step 8: Trial segmentation KMeans
 360   if (0) {
 361     StartNewStep("K means");
 362     // Take input, crop like current mask
 363     typedef CropLikeImageFilter<ImageType> CropLikeFilterType;
 364     typename CropLikeFilterType::Pointer cropLikeFilter = CropLikeFilterType::New();
 365     cropLikeFilter->SetInput(input);
 366     cropLikeFilter->SetCropLikeImage(output);
 367     cropLikeFilter->Update();
 368     ImagePointer working_input = cropLikeFilter->GetOutput();
 369     writeImage<ImageType>(working_input, "crop-input.mhd");
 370     // Set bG at -1000
 371     working_input = clitk::SetBackground<ImageType, MaskImageType>(working_input, output, GetBackgroundValue(), -1000, true);
 372     writeImage<ImageType>(working_input, "crop-input2.mhd");
 373     // Kmeans
 374     typedef itk::ScalarImageKmeansImageFilter<ImageType> KMeansFilterType;
 375     typename KMeansFilterType::Pointer kmeansFilter = KMeansFilterType::New();
 376     kmeansFilter->SetInput(working_input);
 377     //  const unsigned int numberOfInitialClasses = 3;
 378     // const unsigned int useNonContiguousLabels = 0;
 379     kmeansFilter->AddClassWithInitialMean(-1000);
 380     kmeansFilter->AddClassWithInitialMean(30);
 381     kmeansFilter->AddClassWithInitialMean(-40);  // ==> I want this one
 382     DD("Go!");
 383     kmeansFilter->Update();
 384     DD("End");
 385     typename KMeansFilterType::ParametersType estimatedMeans = kmeansFilter->GetFinalMeans();
 386     const unsigned int numberOfClasses = estimatedMeans.Size();
 387     for ( unsigned int i = 0 ; i < numberOfClasses ; ++i ) {
 388       std::cout << "cluster[" << i << "] ";
 389       std::cout << "    estimated mean : " << estimatedMeans[i] << std::endl;
 390     }
 391     MaskImageType::Pointer kmeans = kmeansFilter->GetOutput();
 392     kmeans = clitk::SetBackground<MaskImageType, MaskImageType>(kmeans, kmeans,
 393                                                                 1, GetBackgroundValue(), true);
 394     writeImage<MaskImageType>(kmeans, "kmeans.mhd");
 395     // Get final results, and remove from current mask
 396     boolFilter = BoolFilterType::New();
 397     boolFilter->InPlaceOn();
 398     boolFilter->SetInput1(output);
 399     boolFilter->SetInput2(kmeans);
 400     boolFilter->SetOperationType(BoolFilterType::And);
 401     boolFilter->Update();
 402     output = boolFilter->GetOutput();
 403     writeImage<MaskImageType>(output, "out-kmean.mhd");
 404     this->template StopCurrentStep<MaskImageType>(output);
 405
 406     // TODO -> FillMASK ?
 407     // comment speed ? mask ? 2 class ?
 408
 409
 410     //TODO
 411     // Confidence connected ?
 412
 413   }
 414
 415   //--------------------------------------------------------------------
 416   // Step 8: Lower limits from lung (need separate lung ?)
 417   if (0) {
 418     // StartNewStep("Trial : minus segmented struct");
 419     // MaskImagePointer heart = GetAFDB()->template GetImage <MaskImageType>("heart");
 420     // boolFilter = BoolFilterType::New();
 421     // boolFilter->InPlaceOn();
 422     // boolFilter->SetInput1(output);
 423     // boolFilter->SetInput2(heart);
 424     // boolFilter->SetOperationType(BoolFilterType::AndNot);
 425     // boolFilter->Update();
 426     //  output = boolFilter->GetOutput(); // not needed because InPlace
 427
 428     // Not below the heart
 429     // relPosFilter = RelPosFilterType::New();
 430     // relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId());
 431     // relPosFilter->VerboseStepFlagOff();
 432     // relPosFilter->WriteStepFlagOff();
 433     // relPosFilter->SetInput(output);
 434     // relPosFilter->SetInputObject(heart);
 435     // relPosFilter->SetOrientationType(RelPosFilterType::SupTo);
 436     // relPosFilter->SetIntermediateSpacing(GetIntermediateSpacing());
 437     // relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold3());
 438     // relPosFilter->Update();
 439     // output = relPosFilter->GetOutput();
 440   }
 441
 442   //--------------------------------------------------------------------
 443   // Step 8: Lower limits from lung (need separate lung ?)
 444   if (0) {
 445     StartNewStep("Lower limits with lungs");
 446     // TODO BOFFF ????
 447     relPosFilter = RelPosFilterType::New();
 448     relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId());
 449     relPosFilter->VerboseStepFlagOff();
 450     relPosFilter->WriteStepFlagOff();
 451     relPosFilter->SetInput(output);
 452     //  relPosFilter->SetInputObject(left_lung);
 453     relPosFilter->SetInputObject(lung);
 454     relPosFilter->AddOrientationType(RelPosFilterType::SupTo);
 455     relPosFilter->SetIntermediateSpacing(GetIntermediateSpacing());
 456     relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold3());
 457     relPosFilter->Update();
 458     output = relPosFilter->GetOutput();
 459     this->template StopCurrentStep<MaskImageType>(output);
 460   }
 461
 462   //--------------------------------------------------------------------
 463   // Step 10: Slice by Slice CCL
 464   StartNewStep("Slice by Slice keep only one component");
 465   typedef clitk::ExtractSliceFilter<MaskImageType> ExtractSliceFilterType;
 466   //  typename ExtractSliceFilterType::Pointer
 467   ExtractSliceFilterType::Pointer extractSliceFilter = ExtractSliceFilterType::New();
 468   extractSliceFilter->SetInput(output);
 469   extractSliceFilter->SetDirection(2);
 470   extractSliceFilter->Update();
 471   typedef typename ExtractSliceFilterType::SliceType SliceType;
 472   std::vector<typename SliceType::Pointer> mSlices;
 473   extractSliceFilter->GetOutputSlices(mSlices);
 474   for(unsigned int i=0; i<mSlices.size(); i++) {
 475     mSlices[i] = Labelize<SliceType>(mSlices[i], 0, true, 100);
 476     mSlices[i] = KeepLabels<SliceType>(mSlices[i], 0, 1, 1, 1, true);
 477   }
 478   typedef itk::JoinSeriesImageFilter<SliceType, MaskImageType> JoinSeriesFilterType;
 479   typename JoinSeriesFilterType::Pointer joinFilter = JoinSeriesFilterType::New();
 480   joinFilter->SetOrigin(output->GetOrigin()[2]);
 481   joinFilter->SetSpacing(output->GetSpacing()[2]);
 482   for(unsigned int i=0; i<mSlices.size(); i++) {
 483     joinFilter->PushBackInput(mSlices[i]);
 484   }
 485   joinFilter->Update();
 486   output = joinFilter->GetOutput();
 487   this->template StopCurrentStep<MaskImageType>(output);
 488
 489   //--------------------------------------------------------------------
 490   // Step 9: Binarize to remove too high HU
 491   // --> warning CCL slice by slice must be done before
 492   if (0) {
 493     StartNewStep("Remove hypersignal (bones and injected part");
 494     // Crop initial ct like current support
 495     typedef CropLikeImageFilter<ImageType> CropLikeFilterType;
 496     typename CropLikeFilterType::Pointer cropLikeFilter = CropLikeFilterType::New();
 497     cropLikeFilter->SetInput(input);
 498     cropLikeFilter->SetCropLikeImage(output);
 499     cropLikeFilter->Update();
 500     ImagePointer working_input = cropLikeFilter->GetOutput();
 501     //  writeImage<ImageType>(working_input, "crop-ct.mhd");
 502     // Binarize
 503     typedef itk::BinaryThresholdImageFilter<ImageType, MaskImageType> InputBinarizeFilterType;
 504     typename InputBinarizeFilterType::Pointer binarizeFilter=InputBinarizeFilterType::New();
 505     binarizeFilter->SetInput(working_input);
 506     binarizeFilter->SetLowerThreshold(GetLowerThreshold());
 507     binarizeFilter->SetUpperThreshold(GetUpperThreshold());
 508     binarizeFilter->SetInsideValue(this->GetBackgroundValue());  // opposite
 509     binarizeFilter->SetOutsideValue(this->GetForegroundValue()); // opposite
 510     binarizeFilter->Update();
 511     MaskImagePointer working_bin = binarizeFilter->GetOutput();
 512     // writeImage<MaskImageType>(working_bin, "bin.mhd");
 513     // Remove from support
 514     boolFilter = BoolFilterType::New();
 515     boolFilter->InPlaceOn();
 516     boolFilter->SetInput1(output);
 517     boolFilter->SetInput2(working_bin);
 518     boolFilter->SetOperationType(BoolFilterType::AndNot);
 519     boolFilter->Update();
 520     output = boolFilter->GetOutput();
 521     StopCurrentStep<MaskImageType>(output);
 522   }
 523
 524   //--------------------------------------------------------------------
 525   // Step 10 : AutoCrop
 526   StartNewStep("AutoCrop");
 527   output = clitk::AutoCrop<MaskImageType>(output, GetBackgroundValue());
 528   this->template StopCurrentStep<MaskImageType>(output);
 529
 530   // Bones ? pb with RAM ? FillHoles ?
 531
 532   // how to do with post part ? spine /lung ?
 533   // POST the spine (should be separated from the rest)
 534   /// DO THAT ---->>
 535   // histo Y on the whole bones_post (3D) -> result is the Y center on the spine (?)
 536   // by slice on bones_post
 537   //       find the most ant point in the center
 538   //       from this point go to post until out of bones.
 539   //
 540
 541
 542   // End, set the real size
 543   this->GetOutput(0)->SetRegions(output->GetLargestPossibleRegion());
 544   this->GetOutput(0)->SetLargestPossibleRegion(output->GetLargestPossibleRegion());
 545   this->GetOutput(0)->SetRequestedRegion(output->GetLargestPossibleRegion());
 546   this->GetOutput(0)->SetBufferedRegion(output->GetLargestPossibleRegion());
 547 }
 548 //--------------------------------------------------------------------
 549
 550
 551 //--------------------------------------------------------------------
 552 template <class ImageType>
 553 void
 554 clitk::ExtractMediastinumFilter<ImageType>::
 555 GenerateData()
 556 {
 557   this->GraftOutput(output);
 558   // Store image filenames into AFDB
 559   GetAFDB()->SetImageFilename("Mediastinum", this->GetOutputMediastinumFilename());
 560   WriteAFDB();
 561 }
 562 //--------------------------------------------------------------------
 563
 564
 565 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX