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