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://www.centreleonberard.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::StructuresExtractionFilter<ImageType>()
  49 {
  50   this->SetNumberOfRequiredInputs(1);
  51   SetBackgroundValuePatient(0);
  52   SetBackgroundValueLung(0);
  53   SetBackgroundValueBones(0);
  54   SetForegroundValueLeftLung(1);
  55   SetForegroundValueRightLung(2);
  56   SetDistanceMaxToAnteriorPartOfTheVertebralBody(10);
  57   SetOutputMediastinumFilename("mediastinum.mha");
  58   UseBonesOff();
  59 }
  60 //--------------------------------------------------------------------
  61
  62
  63 //--------------------------------------------------------------------
  64 template <class ImageType>
  65 void
  66 clitk::ExtractMediastinumFilter<ImageType>::
  67 SetInputPatientLabelImage(const MaskImageType * image, MaskImagePixelType bg)
  68 {
  69   this->SetNthInput(0, const_cast<MaskImageType *>(image));
  70   m_BackgroundValuePatient = bg;
  71 }
  72 //--------------------------------------------------------------------
  73
  74
  75 //--------------------------------------------------------------------
  76 template <class ImageType>
  77 void
  78 clitk::ExtractMediastinumFilter<ImageType>::
  79 SetInputLungLabelImage(const MaskImageType * image, MaskImagePixelType bg,
  80                        MaskImagePixelType fgLeft, MaskImagePixelType fgRight)
  81 {
  82   this->SetNthInput(1, const_cast<MaskImageType *>(image));
  83   m_BackgroundValueLung = bg;
  84   m_ForegroundValueLeftLung = fgLeft;
  85   m_ForegroundValueRightLung = fgRight;
  86 }
  87 //--------------------------------------------------------------------
  88
  89
  90 //--------------------------------------------------------------------
  91 template <class ImageType>
  92 void
  93 clitk::ExtractMediastinumFilter<ImageType>::
  94 SetInputBonesLabelImage(const MaskImageType * image, MaskImagePixelType bg)
  95 {
  96   this->SetNthInput(2, const_cast<MaskImageType *>(image));
  97   m_BackgroundValueBones = bg;
  98 }
  99 //--------------------------------------------------------------------
 100
 101
 102 //--------------------------------------------------------------------
 103 template <class ImageType>
 104 void
 105 clitk::ExtractMediastinumFilter<ImageType>::
 106 SetInputTracheaLabelImage(const MaskImageType * image, MaskImagePixelType bg)
 107 {
 108   this->SetNthInput(3, const_cast<MaskImageType *>(image));
 109   m_BackgroundValueTrachea = bg;
 110 }
 111 //--------------------------------------------------------------------
 112
 113 //--------------------------------------------------------------------
 114 template <class ImageType>
 115 void
 116 clitk::ExtractMediastinumFilter<ImageType>::
 117 SetFuzzyThreshold(std::string tag, double value)
 118 {
 119   m_FuzzyThreshold[tag] = value;
 120 }
 121 //--------------------------------------------------------------------
 122
 123
 124 //--------------------------------------------------------------------
 125 template <class ImageType>
 126 double
 127 clitk::ExtractMediastinumFilter<ImageType>::
 128 GetFuzzyThreshold(std::string tag)
 129 {
 130   if (m_FuzzyThreshold.find(tag) == m_FuzzyThreshold.end()) {
 131     clitkExceptionMacro("Could not find options "+tag+" in the list of FuzzyThresholds.");
 132     return 0.0;
 133   }
 134
 135   return m_FuzzyThreshold[tag];
 136 }
 137 //--------------------------------------------------------------------
 138
 139
 140 //--------------------------------------------------------------------
 141 template <class ImageType>
 142 void
 143 clitk::ExtractMediastinumFilter<ImageType>::
 144 GenerateInputRequestedRegion()
 145 {
 146   // DD("GenerateInputRequestedRegion");
 147   // Do not call default
 148   // Superclass::GenerateInputRequestedRegion();
 149   // DD("End GenerateInputRequestedRegion");
 150 }
 151 //--------------------------------------------------------------------
 152
 153
 154 //--------------------------------------------------------------------
 155 template <class ImageType>
 156 void
 157 clitk::ExtractMediastinumFilter<ImageType>::
 158 SetInput(const ImageType * image)
 159 {
 160   this->SetNthInput(0, const_cast<ImageType *>(image));
 161 }
 162 //--------------------------------------------------------------------
 163
 164
 165
 166 //--------------------------------------------------------------------
 167 template <class ImageType>
 168 void
 169 clitk::ExtractMediastinumFilter<ImageType>::
 170 GenerateOutputInformation() {
 171   // Do not call default
 172   // Superclass::GenerateOutputInformation();
 173
 174   //--------------------------------------------------------------------
 175   // Get input pointers
 176   this->LoadAFDB();
 177   ImageConstPointer input = dynamic_cast<const ImageType*>(itk::ProcessObject::GetInput(0));
 178   MaskImagePointer patient, lung, bones, trachea;
 179   patient = this->GetAFDB()->template GetImage <MaskImageType>("Patient");
 180   lung = this->GetAFDB()->template GetImage <MaskImageType>("Lungs");
 181   if (GetUseBones()) {
 182     bones = this->GetAFDB()->template GetImage <MaskImageType>("Bones");
 183   }
 184   trachea = this->GetAFDB()->template GetImage <MaskImageType>("Trachea");
 185
 186   //this->VerboseImageSizeFlagOn();
 187
 188   //--------------------------------------------------------------------
 189   // Step : Crop support (patient) to lung extend in RL
 190   this->StartNewStep("[Mediastinum] Crop support like lungs along LR");
 191   typedef clitk::CropLikeImageFilter<MaskImageType> CropFilterType;
 192   typename CropFilterType::Pointer cropFilter = CropFilterType::New();
 193   cropFilter->SetInput(patient);
 194   cropFilter->SetCropLikeImage(lung, 0);// Indicate that we only crop in X (Left-Right) axe
 195   cropFilter->Update();
 196   output = cropFilter->GetOutput();
 197   this->template StopCurrentStep<MaskImageType>(output);
 198
 199   //--------------------------------------------------------------------
 200   // Step : Crop support (previous) to bones extend in AP
 201   if (GetUseBones()) {
 202     this->StartNewStep("[Mediastinum] Crop support like bones along AP");
 203     cropFilter = CropFilterType::New();
 204     cropFilter->SetInput(output);
 205     cropFilter->SetCropLikeImage(bones, 1);// Indicate that we only crop in Y (Ant-Post) axe
 206     cropFilter->Update();
 207     output = cropFilter->GetOutput();
 208     this->template StopCurrentStep<MaskImageType>(output);
 209   }
 210
 211   //--------------------------------------------------------------------
 212   // Step : patient minus lungs, minus bones, minus trachea
 213   this->StartNewStep("[Mediastinum] Patient contours minus lungs, trachea [and bones]");
 214   typedef clitk::BooleanOperatorLabelImageFilter<MaskImageType> BoolFilterType;
 215   typename BoolFilterType::Pointer boolFilter = BoolFilterType::New();
 216   boolFilter->InPlaceOn();
 217   boolFilter->SetInput1(output);
 218   boolFilter->SetInput2(lung);
 219   boolFilter->SetOperationType(BoolFilterType::AndNot);
 220   boolFilter->Update();
 221   if (GetUseBones()) {
 222     boolFilter->SetInput1(boolFilter->GetOutput());
 223     boolFilter->SetInput2(bones);
 224     boolFilter->SetOperationType(BoolFilterType::AndNot);
 225     boolFilter->Update();
 226   }
 227   boolFilter->SetInput1(boolFilter->GetOutput());
 228   boolFilter->SetInput2(trachea);
 229   boolFilter->SetOperationType(BoolFilterType::AndNot);
 230   boolFilter->Update();
 231   output = boolFilter->GetOutput();
 232
 233   // Auto crop to gain support area
 234   output = clitk::AutoCrop<MaskImageType>(output, this->GetBackgroundValue());
 235   this->template StopCurrentStep<MaskImageType>(output);
 236
 237   //--------------------------------------------------------------------
 238   // Step : LR limits from lung (need separate lung ?)
 239   // Get separate lung images to get only the right and left lung
 240   // (because RelativePositionPropImageFilter only consider fg=1);
 241   // (label must be '1' because right is greater than left).  (WE DO
 242   // NOT NEED TO SEPARATE ? )
 243   this->StartNewStep("[Mediastinum] Left/Right limits with lungs");
 244
 245   // The following cannot be "inplace" because mask is the same than input ...
 246   MaskImagePointer right_lung =
 247     clitk::SetBackground<MaskImageType, MaskImageType>(lung, lung, 2, 0, false);
 248   MaskImagePointer left_lung =
 249     clitk::SetBackground<MaskImageType, MaskImageType>(lung, lung, 1, 0, false);
 250   left_lung = clitk::SetBackground<MaskImageType, MaskImageType>(left_lung, left_lung, 2, 1, false);
 251   right_lung = clitk::ResizeImageLike<MaskImageType>(right_lung, output, this->GetBackgroundValue());
 252   left_lung = clitk::ResizeImageLike<MaskImageType>(left_lung, output, this->GetBackgroundValue());
 253   this->GetAFDB()->template SetImage<MaskImageType>("RightLung", "seg/RightLung.mha",
 254                                                     right_lung, true);
 255   this->GetAFDB()->template SetImage<MaskImageType>("LeftLung", "seg/LeftLung.mha",
 256                                                     left_lung, true);
 257   this->GetAFDB()->Write();
 258   this->template StopCurrentStep<MaskImageType>(output);
 259
 260   //--------------------------------------------------------------------
 261   // Step : AP limits from bones
 262   // Separate the bones in the ant-post middle
 263   MaskImagePointer bones_ant;
 264   MaskImagePointer bones_post;
 265   MaskImagePointType middle_AntPost_position;
 266   middle_AntPost_position.Fill(NumericTraits<MaskImagePointType::ValueType>::Zero);
 267   if (GetUseBones()) {
 268     this->StartNewStep("[Mediastinum] Ant/Post limits with bones");
 269
 270     // To define ant and post part of the bones with a single horizontal line
 271     MaskImageIndexType index_bones_middle;
 272
 273     // Method1: cut in the middle (not optimal)
 274     /*
 275     middle_AntPost_position[0] = middle_AntPost_position[2] = 0;
 276     middle_AntPost_position[1] = bones->GetOrigin()[1]+
 277     (bones->GetLargestPossibleRegion().GetSize()[1]*bones->GetSpacing()[1])/2.0;
 278     DD(middle_AntPost_position);
 279     bones->TransformPhysicalPointToIndex(middle_AntPost_position, index_bones_middle);
 280     */
 281
 282     // Method2: Use VertebralBody, take most ant point
 283     MaskImagePointer VertebralBody = this->GetAFDB()->template GetImage<MaskImageType>("VertebralBody");
 284     FindExtremaPointInAGivenDirection<MaskImageType>(VertebralBody, this->GetBackgroundValue(),
 285                                                      1, true, middle_AntPost_position);
 286     bones->TransformPhysicalPointToIndex(middle_AntPost_position, index_bones_middle);
 287
 288     // Split bone image first into two parts (ant and post), and crop
 289     // lateraly to get vertebral
 290     typedef itk::RegionOfInterestImageFilter<MaskImageType, MaskImageType> ROIFilterType;
 291     //  typedef itk::ExtractImageFilter<MaskImageType, MaskImageType> ROIFilterType;
 292     typename ROIFilterType::Pointer roiFilter = ROIFilterType::New();
 293     MaskImageRegionType region = bones->GetLargestPossibleRegion();
 294     MaskImageSizeType size = region.GetSize();
 295     MaskImageIndexType index = region.GetIndex();
 296     // ANT part
 297     // crop LR to keep 1/4 center part
 298     //    index[0] = size[0]/4+size[0]/8;
 299     //size[0] = size[0]/4;
 300     // crop AP to keep first (ant) part
 301     size[1] =  index_bones_middle[1]; //size[1]/2.0;
 302     region.SetSize(size);
 303     region.SetIndex(index);
 304     roiFilter->SetInput(bones);
 305     roiFilter->SetRegionOfInterest(region);
 306     roiFilter->ReleaseDataFlagOff();
 307     roiFilter->Update();
 308     bones_ant = roiFilter->GetOutput();
 309     //      writeImage<MaskImageType>(bones_ant, "b_ant.mhd");
 310     // POST part
 311     roiFilter = ROIFilterType::New();
 312     index[1] = bones->GetLargestPossibleRegion().GetIndex()[1] + size[1]-1;
 313     size[1] =  bones->GetLargestPossibleRegion().GetSize()[1] - size[1];
 314     region.SetIndex(index);
 315     region.SetSize(size);
 316     roiFilter->SetInput(bones);
 317     roiFilter->SetRegionOfInterest(region);
 318     roiFilter->ReleaseDataFlagOff();
 319     roiFilter->Update();
 320     bones_post = roiFilter->GetOutput();
 321     // writeImage<MaskImageType>(bones_post, "b_post.mhd");
 322
 323     // Now, insert this image in the AFDB ==> (needed because used in the RelativePosition config file)
 324     this->GetAFDB()->template SetImage<MaskImageType>("Bones_Post", "seg/Bones_Post.mha",
 325                                                  bones_post, true);
 326     this->GetAFDB()->template SetImage<MaskImageType>("Bones_Ant", "seg/Bones_Ant.mha",
 327                                                  bones_ant, true);
 328     this->GetAFDB()->Write();
 329
 330     this->template StopCurrentStep<MaskImageType>(output);
 331   }
 332
 333   //--------------------------------------------------------------------
 334   // Remove VertebralBody part
 335   this->StartNewStep("[Mediastinum] Remove VertebralBody");
 336   MaskImagePointer VertebralBody = this->GetAFDB()->template GetImage<MaskImageType>("VertebralBody");
 337   clitk::AndNot<MaskImageType>(output, VertebralBody, this->GetBackgroundValue());
 338   this->template StopCurrentStep<MaskImageType>(output);
 339
 340   //--------------------------------------------------------------------
 341   // Generic RelativePosition processes
 342   output = this->ApplyRelativePositionList("Mediastinum", output);
 343
 344
 345   //--------------------------------------------------------------------
 346   // FIXME --> do not put this limits here !
 347   /*
 348   // Step : SI limits It is better to do this limit *AFTER* the
 349   // RelativePosition to avoid some issue due to superior boundaries.
 350   this->StartNewStep("[Mediastinum] Keep inferior to CricoidCartilag");
 351   // load Cricoid, get centroid, cut above (or below), lower bound
 352   MaskImagePointType p;
 353   try {
 354     MaskImagePointer CricoidCartilag = this->GetAFDB()->template GetImage <MaskImageType>("CricoidCartilag");
 355     p[0] = p[1] = p[2] =  0.0; // to avoid warning
 356     clitk::FindExtremaPointInAGivenDirection<MaskImageType>(CricoidCartilag,
 357                                                             this->GetBackgroundValue(), 2, true, p);
 358   } catch (clitk::ExceptionObject e) {
 359     //DD("CricoidCartilag image not found, try CricoidCartilagZ");
 360     this->GetAFDB()->GetPoint3D("CricoidCartilagPoint", p);
 361   }
 362   output = clitk::CropImageRemoveGreaterThan<MaskImageType>(output, 2, p[2], true, this->GetBackgroundValue());
 363   this->template StopCurrentStep<MaskImageType>(output);
 364   */
 365
 366   //--------------------------------------------------------------------
 367   // Step: Get CCL
 368   this->StartNewStep("[Mediastinum] Keep main connected component");
 369   output = clitk::Labelize<MaskImageType>(output, this->GetBackgroundValue(), false, 500);
 370   // output = RemoveLabels<MaskImageType>(output, BG, param->GetLabelsToRemove());
 371   output = clitk::KeepLabels<MaskImageType>(output, this->GetBackgroundValue(),
 372                                             this->GetForegroundValue(), 1, 1, 0);
 373   this->template StopCurrentStep<MaskImageType>(output);
 374
 375   //--------------------------------------------------------------------
 376   // Step: Remove post part from VertebralBody
 377   this->StartNewStep("[Mediastinum] Remove post part according to VertebralBody");
 378   RemovePostPartOfVertebralBody();
 379   this->template StopCurrentStep<MaskImageType>(output);
 380
 381   //--------------------------------------------------------------------
 382   // Step: Slice by Slice CCL
 383   this->StartNewStep("[Mediastinum] Slice by Slice keep only one component");
 384   typedef clitk::ExtractSliceFilter<MaskImageType> ExtractSliceFilterType;
 385   //  typename ExtractSliceFilterType::Pointer
 386   ExtractSliceFilterType::Pointer extractSliceFilter = ExtractSliceFilterType::New();
 387   extractSliceFilter->SetInput(output);
 388   extractSliceFilter->SetDirection(2);
 389   extractSliceFilter->Update();
 390   typedef typename ExtractSliceFilterType::SliceType SliceType;
 391   std::vector<typename SliceType::Pointer> mSlices;
 392   extractSliceFilter->GetOutputSlices(mSlices);
 393   for(unsigned int i=0; i<mSlices.size(); i++) {
 394     mSlices[i] = Labelize<SliceType>(mSlices[i], 0, true, 100);
 395     mSlices[i] = KeepLabels<SliceType>(mSlices[i], 0, 1, 1, 1, true);
 396   }
 397   typedef itk::JoinSeriesImageFilter<SliceType, MaskImageType> JoinSeriesFilterType;
 398   typename JoinSeriesFilterType::Pointer joinFilter = JoinSeriesFilterType::New();
 399   joinFilter->SetOrigin(output->GetOrigin()[2]);
 400   joinFilter->SetSpacing(output->GetSpacing()[2]);
 401   for(unsigned int i=0; i<mSlices.size(); i++) {
 402     joinFilter->PushBackInput(mSlices[i]);
 403   }
 404   joinFilter->Update();
 405   output = joinFilter->GetOutput();
 406   this->template StopCurrentStep<MaskImageType>(output);
 407
 408   //--------------------------------------------------------------------
 409   // Step 10 : AutoCrop
 410   this->StartNewStep("[Mediastinum] AutoCrop");
 411   output = clitk::AutoCrop<MaskImageType>(output, this->GetBackgroundValue());
 412   this->template StopCurrentStep<MaskImageType>(output);
 413
 414   // End, set the real size
 415   this->GetOutput(0)->SetRegions(output->GetLargestPossibleRegion());
 416   this->GetOutput(0)->SetLargestPossibleRegion(output->GetLargestPossibleRegion());
 417   this->GetOutput(0)->SetRequestedRegion(output->GetLargestPossibleRegion());
 418   this->GetOutput(0)->SetBufferedRegion(output->GetLargestPossibleRegion());
 419 }
 420 //--------------------------------------------------------------------
 421
 422
 423 //--------------------------------------------------------------------
 424 template <class ImageType>
 425 void
 426 clitk::ExtractMediastinumFilter<ImageType>::
 427 GenerateData()
 428 {
 429   this->GraftOutput(output);
 430   // Store image filenames into AFDB
 431   this->GetAFDB()->SetImageFilename("Mediastinum", this->GetOutputMediastinumFilename());
 432   this->WriteAFDB();
 433 }
 434 //--------------------------------------------------------------------
 435
 436
 437 //--------------------------------------------------------------------
 438 template <class ImageType>
 439 void
 440 clitk::ExtractMediastinumFilter<ImageType>::
 441 RemovePostPartOfVertebralBody()
 442 {
 443
 444   /*
 445     Posteriorly, Station 8 abuts the descending aorta and the anterior
 446     aspect of the vertebral body until an imaginary horizontal line
 447     running 1 cm posterior to the anterior border of the vertebral
 448     body (Fig. 3C).
 449
 450     => We use this definition for all the mediastinum
 451
 452    Find most Ant point in VertebralBody. Consider the horizontal line
 453    which is 'DistanceMaxToAnteriorPartOfTheVertebralBody' away from
 454    the most ant point.
 455   */
 456
 457   // Get VertebralBody mask image
 458   MaskImagePointer VertebralBody =
 459     this->GetAFDB()->template GetImage <MaskImageType>("VertebralBody");
 460
 461   // Consider vertebral body slice by slice
 462   std::vector<MaskSlicePointer> vertebralSlices;
 463   clitk::ExtractSlices<MaskImageType>(VertebralBody, 2, vertebralSlices);
 464
 465   // For each slice, compute the most anterior point
 466   std::map<int, MaskSlicePointType> vertebralAntPositionBySlice;
 467   for(uint i=0; i<vertebralSlices.size(); i++) {
 468     MaskSlicePointType p;
 469     bool found = clitk::FindExtremaPointInAGivenDirection<MaskSliceType>(vertebralSlices[i],
 470                                                                      this->GetBackgroundValue(),
 471                                                                      1, true, p);
 472     if (found) {
 473       vertebralAntPositionBySlice[i] = p;
 474     }
 475     else {
 476       // It should not happen ! But sometimes, a contour is missing or
 477       // the VertebralBody is not delineated enough inferiorly ... in
 478       // those cases, we consider the first found slice.
 479       //        std::cerr << "No foreground pixels in this VertebralBody slices !?? I try with the previous/next slice" << std::endl;
 480       // [ Possible alternative -> consider previous limit ]
 481     }
 482   }
 483
 484   // Convert 2D points in slice into 3D points
 485   std::vector<MaskImagePointType> vertebralAntPositions;
 486   clitk::PointsUtils<MaskImageType>::Convert2DMapTo3DList(vertebralAntPositionBySlice,
 487                                                           VertebralBody,
 488                                                           vertebralAntPositions);
 489
 490   // DEBUG : write list of points
 491   clitk::WriteListOfLandmarks<MaskImageType>(vertebralAntPositions,
 492                                              "VertebralBodyMostAnteriorPoints.txt");
 493
 494   // Cut support posteriorly 1cm the most anterior point of the
 495   // VertebralBody. Do nothing below and above the VertebralBody. To
 496   // do that compute several region, slice by slice and fill.
 497   MaskImageRegionType region;
 498   MaskImageSizeType size;
 499   MaskImageIndexType start;
 500   size[2] = 1; // a single slice
 501   start[0] = output->GetLargestPossibleRegion().GetIndex()[0];
 502   size[0] = output->GetLargestPossibleRegion().GetSize()[0];
 503   for(uint i=0; i<vertebralAntPositions.size(); i++) {
 504     typedef typename itk::ImageRegionIterator<MaskImageType> IteratorType;
 505     IteratorType iter =
 506       IteratorType(output, output->GetLargestPossibleRegion());
 507     MaskImageIndexType index;
 508     // Consider some cm posterior to most anterior positions (usually
 509     // 1 cm).
 510     vertebralAntPositions[i][1] += GetDistanceMaxToAnteriorPartOfTheVertebralBody();
 511     // Get index of this point
 512     output->TransformPhysicalPointToIndex(vertebralAntPositions[i], index);
 513     // Compute region (a single slice)
 514     start[2] = index[2];
 515     start[1] = output->GetLargestPossibleRegion().GetIndex()[1]+index[1];
 516     size[1] = output->GetLargestPossibleRegion().GetSize()[1]-start[1];
 517     region.SetSize(size);
 518     region.SetIndex(start);
 519     // Fill region
 520     if (output->GetLargestPossibleRegion().IsInside(start))  {
 521       itk::ImageRegionIterator<MaskImageType> it(output, region);
 522       it.GoToBegin();
 523       while (!it.IsAtEnd()) {
 524         it.Set(this->GetBackgroundValue());
 525         ++it;
 526       }
 527     }
 528   }
 529 }
 530 //--------------------------------------------------------------------
 531
 532
 533 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX