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   SetDistanceMaxToAnteriorPartOfTheSpine(10);
  68   SetOutputMediastinumFilename("mediastinum.mhd");
  69
  70   UseBonesOff();
  71 }
  72 //--------------------------------------------------------------------
  73
  74
  75 //--------------------------------------------------------------------
  76 template <class ImageType>
  77 void
  78 clitk::ExtractMediastinumFilter<ImageType>::
  79 SetInputPatientLabelImage(const MaskImageType * image, MaskImagePixelType bg)
  80 {
  81   this->SetNthInput(0, const_cast<MaskImageType *>(image));
  82   m_BackgroundValuePatient = bg;
  83 }
  84 //--------------------------------------------------------------------
  85
  86
  87 //--------------------------------------------------------------------
  88 template <class ImageType>
  89 void
  90 clitk::ExtractMediastinumFilter<ImageType>::
  91 SetInputLungLabelImage(const MaskImageType * image, MaskImagePixelType bg,
  92                        MaskImagePixelType fgLeft, MaskImagePixelType fgRight)
  93 {
  94   this->SetNthInput(1, const_cast<MaskImageType *>(image));
  95   m_BackgroundValueLung = bg;
  96   m_ForegroundValueLeftLung = fgLeft;
  97   m_ForegroundValueRightLung = fgRight;
  98 }
  99 //--------------------------------------------------------------------
 100
 101
 102 //--------------------------------------------------------------------
 103 template <class ImageType>
 104 void
 105 clitk::ExtractMediastinumFilter<ImageType>::
 106 SetInputBonesLabelImage(const MaskImageType * image, MaskImagePixelType bg)
 107 {
 108   this->SetNthInput(2, const_cast<MaskImageType *>(image));
 109   m_BackgroundValueBones = bg;
 110 }
 111 //--------------------------------------------------------------------
 112
 113
 114 //--------------------------------------------------------------------
 115 template <class ImageType>
 116 void
 117 clitk::ExtractMediastinumFilter<ImageType>::
 118 SetInputTracheaLabelImage(const MaskImageType * image, MaskImagePixelType bg)
 119 {
 120   this->SetNthInput(3, const_cast<MaskImageType *>(image));
 121   m_BackgroundValueTrachea = bg;
 122 }
 123 //--------------------------------------------------------------------
 124
 125
 126 //--------------------------------------------------------------------
 127 template <class ImageType>
 128 template<class ArgsInfoType>
 129 void
 130 clitk::ExtractMediastinumFilter<ImageType>::
 131 SetArgsInfo(ArgsInfoType mArgsInfo)
 132 {
 133   SetVerboseOption_GGO(mArgsInfo);
 134   SetVerboseStep_GGO(mArgsInfo);
 135   SetWriteStep_GGO(mArgsInfo);
 136   SetVerboseWarningOff_GGO(mArgsInfo);
 137
 138   SetIntermediateSpacing_GGO(mArgsInfo);
 139   SetFuzzyThreshold1_GGO(mArgsInfo);
 140   SetFuzzyThreshold2_GGO(mArgsInfo);
 141   SetFuzzyThreshold3_GGO(mArgsInfo);
 142
 143   SetAFDBFilename_GGO(mArgsInfo);
 144   SetDistanceMaxToAnteriorPartOfTheSpine_GGO(mArgsInfo);
 145   SetUseBones_GGO(mArgsInfo);
 146
 147   SetLowerThreshold_GGO(mArgsInfo);
 148   SetUpperThreshold_GGO(mArgsInfo);
 149 }
 150 //--------------------------------------------------------------------
 151
 152
 153 //--------------------------------------------------------------------
 154 template <class ImageType>
 155 void
 156 clitk::ExtractMediastinumFilter<ImageType>::
 157 GenerateInputRequestedRegion()
 158 {
 159   //DD("GenerateInputRequestedRegion");
 160   // Do not call default
 161   //  Superclass::GenerateInputRequestedRegion();
 162   //  DD("End GenerateInputRequestedRegion");
 163 }
 164
 165 //--------------------------------------------------------------------
 166
 167
 168 //--------------------------------------------------------------------
 169 template <class ImageType>
 170 void
 171 clitk::ExtractMediastinumFilter<ImageType>::
 172 SetInput(const ImageType * image)
 173 {
 174   this->SetNthInput(0, const_cast<ImageType *>(image));
 175 }
 176 //--------------------------------------------------------------------
 177
 178
 179 //--------------------------------------------------------------------
 180 template <class ImageType>
 181 void
 182 clitk::ExtractMediastinumFilter<ImageType>::
 183 GenerateOutputInformation() {
 184   //  DD("GenerateOutputInformation");
 185   // Do not call default
 186   //  Superclass::GenerateOutputInformation();
 187
 188   //--------------------------------------------------------------------
 189   // Get input pointers
 190   LoadAFDB();
 191   ImageConstPointer input = dynamic_cast<const ImageType*>(itk::ProcessObject::GetInput(0));
 192   MaskImagePointer patient = GetAFDB()->template GetImage <MaskImageType>("patient");
 193   MaskImagePointer lung = GetAFDB()->template GetImage <MaskImageType>("lungs");
 194   MaskImagePointer bones = GetAFDB()->template GetImage <MaskImageType>("bones");
 195   MaskImagePointer trachea = GetAFDB()->template GetImage <MaskImageType>("trachea");
 196
 197   //--------------------------------------------------------------------
 198   // Step 1: Crop support (patient) to lung extend in RL
 199   StartNewStep("Crop support like lungs along LR");
 200   typedef clitk::CropLikeImageFilter<MaskImageType> CropFilterType;
 201   typename CropFilterType::Pointer cropFilter = CropFilterType::New();
 202   cropFilter->SetInput(patient);
 203   cropFilter->SetCropLikeImage(lung, 0);// Indicate that we only crop in X (Left-Right) axe
 204   cropFilter->Update();
 205   output = cropFilter->GetOutput();
 206   this->template StopCurrentStep<MaskImageType>(output);
 207
 208   //--------------------------------------------------------------------
 209   // Step 2: Crop support (previous) to bones extend in AP
 210   if (GetUseBones()) {
 211     StartNewStep("Crop support like bones along AP");
 212     cropFilter = CropFilterType::New();
 213     cropFilter->SetInput(output);
 214     cropFilter->SetCropLikeImage(bones, 1);// Indicate that we only crop in Y (Ant-Post) axe
 215     cropFilter->Update();
 216     output = cropFilter->GetOutput();
 217     this->template StopCurrentStep<MaskImageType>(output);
 218   }
 219
 220   //--------------------------------------------------------------------
 221   // Step 3: patient minus lungs, minus bones, minus trachea
 222   StartNewStep("Patient contours minus lungs, bones, trachea");
 223   typedef clitk::BooleanOperatorLabelImageFilter<MaskImageType> BoolFilterType;
 224   typename BoolFilterType::Pointer boolFilter = BoolFilterType::New();
 225   boolFilter->InPlaceOn();
 226   boolFilter->SetInput1(output);
 227   boolFilter->SetInput2(lung);
 228   boolFilter->SetOperationType(BoolFilterType::AndNot);
 229   boolFilter->Update();
 230   if (GetUseBones()) {
 231     boolFilter->SetInput1(boolFilter->GetOutput());
 232     boolFilter->SetInput2(bones);
 233     boolFilter->SetOperationType(BoolFilterType::AndNot);
 234     boolFilter->Update();
 235   }
 236   boolFilter->SetInput1(boolFilter->GetOutput());
 237   boolFilter->SetInput2(trachea);
 238   boolFilter->SetOperationType(BoolFilterType::AndNot);
 239   boolFilter->Update();
 240   output = boolFilter->GetOutput();
 241
 242   // Auto crop to gain support area
 243   output = clitk::AutoCrop<MaskImageType>(output, GetBackgroundValue());
 244   this->template StopCurrentStep<MaskImageType>(output);
 245
 246   //--------------------------------------------------------------------
 247   // Step 4: LR limits from lung (need separate lung ?)
 248   // Get separate lung images to get only the right and left lung
 249   // (because RelativePositionPropImageFilter only consider fg=1);
 250   // (label must be '1' because right is greater than left).  (WE DO
 251   // NOT NEED TO SEPARATE ? )
 252   StartNewStep("Left/Right limits with lungs");
 253   /*
 254     ImagePointer right_lung = clitk::SetBackground<MaskImageType, MaskImageType>(lung, lung, 2, 0);
 255     ImagePointer left_lung = clitk::SetBackground<MaskImageType, MaskImageType>(lung, lung, 1, 0);
 256     writeImage<MaskImageType>(right_lung, "right.mhd");
 257     writeImage<MaskImageType>(left_lung, "left.mhd");
 258   */
 259   typedef clitk::AddRelativePositionConstraintToLabelImageFilter<MaskImageType> RelPosFilterType;
 260   typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New();
 261   relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId());
 262   relPosFilter->VerboseStepOff();
 263   relPosFilter->WriteStepOff();
 264   relPosFilter->SetInput(output);
 265   //relPosFilter->SetInputObject(left_lung);
 266   relPosFilter->SetInputObject(lung);
 267   relPosFilter->SetOrientationType(RelPosFilterType::LeftTo); // warning left lung is at right ;)
 268   relPosFilter->SetIntermediateSpacing(GetIntermediateSpacing());
 269   relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold1());
 270   relPosFilter->Update();
 271   output = relPosFilter->GetOutput();
 272   //writeImage<MaskImageType>(right_lung, "step4-left.mhd");
 273
 274   relPosFilter->SetInput(output);
 275   relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId());
 276   relPosFilter->VerboseStepOff();
 277   relPosFilter->WriteStepOff();
 278   //relPosFilter->SetInputObject(right_lung);
 279   relPosFilter->SetInputObject(lung);
 280   relPosFilter->SetOrientationType(RelPosFilterType::RightTo);
 281   relPosFilter->SetIntermediateSpacing(GetIntermediateSpacing());
 282   relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold1());
 283   relPosFilter->Update();
 284   output = relPosFilter->GetOutput();
 285   this->template StopCurrentStep<MaskImageType>(output);
 286
 287   //--------------------------------------------------------------------
 288   // Step 5: AP limits from bones
 289   // Separate the bones in the ant-post middle
 290   MaskImageConstPointer bones_ant;
 291   MaskImageConstPointer bones_post;
 292   MaskImagePointType middle_AntPost__position;
 293   if (GetUseBones()) {
 294     StartNewStep("Ant/Post limits with bones");
 295     middle_AntPost__position[0] = middle_AntPost__position[2] = 0;
 296     middle_AntPost__position[1] = bones->GetOrigin()[1]+(bones->GetLargestPossibleRegion().GetSize()[1]*bones->GetSpacing()[1])/2.0;
 297     MaskImageIndexType index_bones_middle;
 298     bones->TransformPhysicalPointToIndex(middle_AntPost__position, index_bones_middle);
 299
 300     // Split bone image first into two parts (ant and post), and crop
 301     // lateraly to get vertebral
 302     typedef itk::RegionOfInterestImageFilter<MaskImageType, MaskImageType> ROIFilterType;
 303     //  typedef itk::ExtractImageFilter<MaskImageType, MaskImageType> ROIFilterType;
 304     typename ROIFilterType::Pointer roiFilter = ROIFilterType::New();
 305     MaskImageRegionType region = bones->GetLargestPossibleRegion();
 306     MaskImageSizeType size = region.GetSize();
 307     MaskImageIndexType index = region.GetIndex();
 308     // ANT part
 309     // crop LR to keep 1/4 center part
 310     index[0] = size[0]/4+size[0]/8;
 311     size[0] = size[0]/4;
 312     // crop AP to keep first (ant) part
 313     size[1] =  index_bones_middle[1]; //size[1]/2.0;
 314     region.SetSize(size);
 315     region.SetIndex(index);
 316     roiFilter->SetInput(bones);
 317     roiFilter->SetRegionOfInterest(region);
 318     roiFilter->ReleaseDataFlagOff();
 319     roiFilter->Update();
 320     bones_ant = roiFilter->GetOutput();
 321     writeImage<MaskImageType>(bones_ant, "b_ant.mhd");
 322     // POST part
 323     roiFilter = ROIFilterType::New();
 324     index[1] = bones->GetLargestPossibleRegion().GetIndex()[1] + size[1]-1;
 325     size[1] =  bones->GetLargestPossibleRegion().GetSize()[1] - size[1];
 326     region.SetIndex(index);
 327     region.SetSize(size);
 328     roiFilter->SetInput(bones);
 329     roiFilter->SetRegionOfInterest(region);
 330     roiFilter->ReleaseDataFlagOff();
 331     roiFilter->Update();
 332     bones_post = roiFilter->GetOutput();
 333     writeImage<MaskImageType>(bones_post, "b_post.mhd");
 334
 335     // Go !
 336     relPosFilter->SetCurrentStepNumber(0);
 337     relPosFilter->ResetPipeline();// = RelPosFilterType::New();
 338     relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId());
 339     relPosFilter->VerboseStepOff();
 340     relPosFilter->WriteStepOff();
 341     relPosFilter->SetInput(output);
 342     relPosFilter->SetInputObject(bones_post);
 343     relPosFilter->SetOrientationType(RelPosFilterType::AntTo);
 344     relPosFilter->SetIntermediateSpacing(GetIntermediateSpacing());
 345     relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold2());
 346     relPosFilter->Update();
 347     output = relPosFilter->GetOutput();
 348     writeImage<MaskImageType>(output, "post.mhd");
 349
 350     relPosFilter->SetInput(relPosFilter->GetOutput());
 351     relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId());
 352     relPosFilter->VerboseStepOff();
 353     relPosFilter->WriteStepOff();
 354     relPosFilter->SetInput(output);
 355     relPosFilter->SetInputObject(bones_ant);
 356     relPosFilter->SetOrientationType(RelPosFilterType::PostTo);
 357     relPosFilter->SetIntermediateSpacing(GetIntermediateSpacing());
 358     relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold2());
 359     relPosFilter->Update();
 360     output = relPosFilter->GetOutput();
 361     this->template StopCurrentStep<MaskImageType>(output);
 362   }
 363
 364   //--------------------------------------------------------------------
 365   // Step 6: Get CCL
 366   StartNewStep("Keep main connected component");
 367   output = clitk::Labelize<MaskImageType>(output, GetBackgroundValue(), false, 500);
 368   // output = RemoveLabels<MaskImageType>(output, BG, param->GetLabelsToRemove());
 369   output = clitk::KeepLabels<MaskImageType>(output, GetBackgroundValue(),
 370                                             GetForegroundValue(), 1, 1, 0);
 371   this->template StopCurrentStep<MaskImageType>(output);
 372
 373   //--------------------------------------------------------------------
 374   // Step 7 : Slice by Slice to optimize posterior part
 375   // Warning slice does not necesseraly correspond between 'output' and 'bones'
 376   typedef clitk::ExtractSliceFilter<MaskImageType> ExtractSliceFilterType;
 377   typename ExtractSliceFilterType::Pointer extractSliceFilter = ExtractSliceFilterType::New();
 378   typedef typename ExtractSliceFilterType::SliceType SliceType;
 379   std::vector<typename SliceType::Pointer> mSlices;
 380   if (GetUseBones()) {
 381     StartNewStep("Rafine posterior part according to vertebral body");
 382     extractSliceFilter->SetInput(bones_post);
 383     extractSliceFilter->SetDirection(2);
 384     extractSliceFilter->Update();
 385     std::vector<double> mVertebralAntPositionBySlice;
 386     extractSliceFilter->GetOutputSlices(mSlices);
 387     for(unsigned int i=0; i<mSlices.size(); i++) {
 388       mSlices[i] = Labelize<SliceType>(mSlices[i], 0, true, 10);
 389       mSlices[i] = KeepLabels<SliceType>(mSlices[i],
 390                                          GetBackgroundValue(),
 391                                          GetForegroundValue(), 1, 2, true); // keep two first
 392       // Find most anterior point (start of the vertebral)
 393       typename itk::ImageRegionIteratorWithIndex<SliceType>
 394         iter(mSlices[i], mSlices[i]->GetLargestPossibleRegion());
 395       iter.GoToBegin();
 396       bool stop = false;
 397       while (!stop) {
 398         if (iter.Get() != GetBackgroundValue())
 399           stop = true; // not foreground because we keep two main label
 400         ++iter;
 401         if (iter.IsAtEnd()) stop = true;
 402       }
 403       if (!iter.IsAtEnd()) {
 404         typename SliceType::PointType p;
 405         mSlices[i]->TransformIndexToPhysicalPoint(iter.GetIndex(),p);
 406         mVertebralAntPositionBySlice.push_back(p[1]);
 407       }
 408       else {
 409         mVertebralAntPositionBySlice.push_back(bones_post->GetOrigin()[1]+(bones->GetLargestPossibleRegion().GetSize()[1]*bones->GetSpacing()[1]));
 410         DD(mVertebralAntPositionBySlice.back());
 411         DD("ERROR ?? NO FG in bones here ?");
 412       }
 413     }
 414
 415     // Cut Post position slice by slice
 416     {
 417       MaskImageRegionType region;
 418       MaskImageSizeType size;
 419       MaskImageIndexType start;
 420       size[2] = 1;
 421       start[0] = output->GetLargestPossibleRegion().GetIndex()[0];
 422       for(unsigned int i=0; i<mSlices.size(); i++) {
 423         // Compute index
 424         MaskImagePointType point;
 425         point[0] = 0;
 426
 427         //TODO 10 mm OPTION
 428
 429         point[1] = mVertebralAntPositionBySlice[i]+GetDistanceMaxToAnteriorPartOfTheSpine();// ADD ONE CM
 430         point[2] = bones_post->GetOrigin()[2]+(bones_post->GetLargestPossibleRegion().GetIndex()[2]+i)*bones_post->GetSpacing()[2];
 431         MaskImageIndexType index;
 432         output->TransformPhysicalPointToIndex(point, index);
 433         // Compute region
 434         start[2] = index[2];
 435         start[1] = output->GetLargestPossibleRegion().GetIndex()[1]+index[1];
 436         size[0] = output->GetLargestPossibleRegion().GetSize()[0];
 437         size[1] = output->GetLargestPossibleRegion().GetSize()[1]-start[1];
 438         region.SetSize(size);
 439         region.SetIndex(start);
 440         // Fill Region
 441         if (output->GetLargestPossibleRegion().IsInside(start))  {
 442           itk::ImageRegionIteratorWithIndex<MaskImageType> it(output, region);
 443           it.GoToBegin();
 444           while (!it.IsAtEnd()) {
 445             it.Set(GetBackgroundValue());
 446             ++it;
 447           }
 448         }
 449       }
 450     }
 451     this->template StopCurrentStep<MaskImageType>(output);
 452   }
 453
 454   //--------------------------------------------------------------------
 455   // Step 8: Trial segmentation KMeans
 456   if (0) {
 457     StartNewStep("K means");
 458     // Take input, crop like current mask
 459     typedef CropLikeImageFilter<ImageType> CropLikeFilterType;
 460     typename CropLikeFilterType::Pointer cropLikeFilter = CropLikeFilterType::New();
 461     cropLikeFilter->SetInput(input);
 462     cropLikeFilter->SetCropLikeImage(output);
 463     cropLikeFilter->Update();
 464     ImagePointer working_input = cropLikeFilter->GetOutput();
 465     writeImage<ImageType>(working_input, "crop-input.mhd");
 466     // Set bG at -1000
 467     working_input = clitk::SetBackground<ImageType, MaskImageType>(working_input, output, GetBackgroundValue(), -1000);
 468     writeImage<ImageType>(working_input, "crop-input2.mhd");
 469     // Kmeans
 470     typedef itk::ScalarImageKmeansImageFilter<ImageType> KMeansFilterType;
 471     typename KMeansFilterType::Pointer kmeansFilter = KMeansFilterType::New();
 472     kmeansFilter->SetInput(working_input);
 473     //  const unsigned int numberOfInitialClasses = 3;
 474     // const unsigned int useNonContiguousLabels = 0;
 475     kmeansFilter->AddClassWithInitialMean(-1000);
 476     kmeansFilter->AddClassWithInitialMean(30);
 477     kmeansFilter->AddClassWithInitialMean(-40);  // ==> I want this one
 478     DD("Go!");
 479     kmeansFilter->Update();
 480     DD("End");
 481     typename KMeansFilterType::ParametersType estimatedMeans = kmeansFilter->GetFinalMeans();
 482     const unsigned int numberOfClasses = estimatedMeans.Size();
 483     for ( unsigned int i = 0 ; i < numberOfClasses ; ++i ) {
 484       std::cout << "cluster[" << i << "] ";
 485       std::cout << "    estimated mean : " << estimatedMeans[i] << std::endl;
 486     }
 487     MaskImageType::Pointer kmeans = kmeansFilter->GetOutput();
 488     kmeans = clitk::SetBackground<MaskImageType, MaskImageType>(kmeans, kmeans,
 489                                                                 1, GetBackgroundValue());
 490     writeImage<MaskImageType>(kmeans, "kmeans.mhd");
 491     // Get final results, and remove from current mask
 492     boolFilter = BoolFilterType::New();
 493     boolFilter->InPlaceOn();
 494     boolFilter->SetInput1(output);
 495     boolFilter->SetInput2(kmeans);
 496     boolFilter->SetOperationType(BoolFilterType::And);
 497     boolFilter->Update();
 498     output = boolFilter->GetOutput();
 499     writeImage<MaskImageType>(output, "out-kmean.mhd");
 500     this->template StopCurrentStep<MaskImageType>(output);
 501
 502     // TODO -> FillMASK ?
 503     // comment speed ? mask ? 2 class ?
 504
 505
 506     //TODO
 507     // Confidence connected ?
 508
 509   }
 510
 511   //--------------------------------------------------------------------
 512   // Step 8: Lower limits from lung (need separate lung ?)
 513   if (1) {
 514     // StartNewStep("Trial : minus segmented struct");
 515     // MaskImagePointer heart = GetAFDB()->template GetImage <MaskImageType>("heart");
 516     // boolFilter = BoolFilterType::New();
 517     // boolFilter->InPlaceOn();
 518     // boolFilter->SetInput1(output);
 519     // boolFilter->SetInput2(heart);
 520     // boolFilter->SetOperationType(BoolFilterType::AndNot);
 521     // boolFilter->Update();
 522     //  output = boolFilter->GetOutput(); // not needed because InPlace
 523
 524     // Not below the heart
 525     // relPosFilter = RelPosFilterType::New();
 526     // relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId());
 527     // relPosFilter->VerboseStepOff();
 528     // relPosFilter->WriteStepOff();
 529     // relPosFilter->SetInput(output);
 530     // relPosFilter->SetInputObject(heart);
 531     // relPosFilter->SetOrientationType(RelPosFilterType::SupTo);
 532     // relPosFilter->SetIntermediateSpacing(GetIntermediateSpacing());
 533     // relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold3());
 534     // relPosFilter->Update();
 535     // output = relPosFilter->GetOutput();
 536   }
 537
 538   //--------------------------------------------------------------------
 539   // Step 8: Lower limits from lung (need separate lung ?)
 540   if (0) {
 541     StartNewStep("Lower limits with lungs");
 542     // TODO BOFFF ????
 543     relPosFilter = RelPosFilterType::New();
 544     relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId());
 545     relPosFilter->VerboseStepOff();
 546     relPosFilter->WriteStepOff();
 547     relPosFilter->SetInput(output);
 548     //  relPosFilter->SetInputObject(left_lung);
 549     relPosFilter->SetInputObject(lung);
 550     relPosFilter->SetOrientationType(RelPosFilterType::SupTo);
 551     relPosFilter->SetIntermediateSpacing(GetIntermediateSpacing());
 552     relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold3());
 553     relPosFilter->Update();
 554     output = relPosFilter->GetOutput();
 555     this->template StopCurrentStep<MaskImageType>(output);
 556   }
 557
 558   //--------------------------------------------------------------------
 559   // Step 10: Slice by Slice CCL
 560   StartNewStep("Slice by Slice keep only one component");
 561   typedef clitk::ExtractSliceFilter<MaskImageType> ExtractSliceFilterType;
 562   //  typename ExtractSliceFilterType::Pointer
 563   extractSliceFilter = ExtractSliceFilterType::New();
 564   extractSliceFilter->SetInput(output);
 565   extractSliceFilter->SetDirection(2);
 566   extractSliceFilter->Update();
 567   typedef typename ExtractSliceFilterType::SliceType SliceType;
 568   //  std::vector<typename SliceType::Pointer>
 569   mSlices.clear();
 570   extractSliceFilter->GetOutputSlices(mSlices);
 571   for(unsigned int i=0; i<mSlices.size(); i++) {
 572     mSlices[i] = Labelize<SliceType>(mSlices[i], 0, true, 100);
 573     mSlices[i] = KeepLabels<SliceType>(mSlices[i], 0, 1, 1, 1, true);
 574   }
 575   typedef itk::JoinSeriesImageFilter<SliceType, MaskImageType> JoinSeriesFilterType;
 576   typename JoinSeriesFilterType::Pointer joinFilter = JoinSeriesFilterType::New();
 577   joinFilter->SetOrigin(output->GetOrigin()[2]);
 578   joinFilter->SetSpacing(output->GetSpacing()[2]);
 579   for(unsigned int i=0; i<mSlices.size(); i++) {
 580     joinFilter->PushBackInput(mSlices[i]);
 581   }
 582   joinFilter->Update();
 583   output = joinFilter->GetOutput();
 584   this->template StopCurrentStep<MaskImageType>(output);
 585
 586   //--------------------------------------------------------------------
 587   // Step 9: Binarize to remove too high HU
 588   // --> warning CCL slice by slice must be done before
 589   StartNewStep("Remove hypersignal (bones and injected part");
 590   // Crop initial ct like current support
 591   typedef CropLikeImageFilter<ImageType> CropLikeFilterType;
 592   typename CropLikeFilterType::Pointer cropLikeFilter = CropLikeFilterType::New();
 593   cropLikeFilter->SetInput(input);
 594   cropLikeFilter->SetCropLikeImage(output);
 595   cropLikeFilter->Update();
 596   ImagePointer working_input = cropLikeFilter->GetOutput();
 597   writeImage<ImageType>(working_input, "crop-ct.mhd");
 598   // Binarize
 599   typedef itk::BinaryThresholdImageFilter<ImageType, MaskImageType> InputBinarizeFilterType;
 600   typename InputBinarizeFilterType::Pointer binarizeFilter=InputBinarizeFilterType::New();
 601   binarizeFilter->SetInput(working_input);
 602   binarizeFilter->SetLowerThreshold(GetLowerThreshold());
 603   binarizeFilter->SetUpperThreshold(GetUpperThreshold());
 604   binarizeFilter->SetInsideValue(this->GetBackgroundValue());  // opposite
 605   binarizeFilter->SetOutsideValue(this->GetForegroundValue()); // opposite
 606   binarizeFilter->Update();
 607   MaskImagePointer working_bin = binarizeFilter->GetOutput();
 608   writeImage<MaskImageType>(working_bin, "bin.mhd");
 609   // Remove from support
 610   boolFilter = BoolFilterType::New();
 611   boolFilter->InPlaceOn();
 612   boolFilter->SetInput1(output);
 613   boolFilter->SetInput2(working_bin);
 614   boolFilter->SetOperationType(BoolFilterType::AndNot);
 615   boolFilter->Update();
 616   output = boolFilter->GetOutput();
 617   StopCurrentStep<MaskImageType>(output);
 618
 619   //--------------------------------------------------------------------
 620   // Step 10 : AutoCrop
 621   StartNewStep("AutoCrop");
 622   output = clitk::AutoCrop<MaskImageType>(output, GetBackgroundValue());
 623   this->template StopCurrentStep<MaskImageType>(output);
 624
 625   // Bones ? pb with RAM ? FillHoles ?
 626
 627   // how to do with post part ? spine /lung ?
 628   // POST the spine (should be separated from the rest)
 629   /// DO THAT ---->>
 630   // histo Y on the whole bones_post (3D) -> result is the Y center on the spine (?)
 631   // by slice on bones_post
 632   //       find the most ant point in the center
 633   //       from this point go to post until out of bones.
 634   //
 635
 636
 637   // End, set the real size
 638   this->GetOutput(0)->SetRegions(output->GetLargestPossibleRegion());
 639   this->GetOutput(0)->SetLargestPossibleRegion(output->GetLargestPossibleRegion());
 640   this->GetOutput(0)->SetRequestedRegion(output->GetLargestPossibleRegion());
 641   this->GetOutput(0)->SetBufferedRegion(output->GetLargestPossibleRegion());
 642 }
 643 //--------------------------------------------------------------------
 644
 645
 646 //--------------------------------------------------------------------
 647 template <class ImageType>
 648 void
 649 clitk::ExtractMediastinumFilter<ImageType>::
 650 GenerateData()
 651 {
 652   DD("GenerateData");
 653   this->GraftOutput(output);
 654   // Store image filenames into AFDB
 655   GetAFDB()->SetImageFilename("mediastinum", this->GetOutputMediastinumFilename());
 656   WriteAFDB();
 657 }
 658 //--------------------------------------------------------------------
 659
 660
 661 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX