itk/clitkSegmentationUtils.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 // clitk
  20 #include "clitkSetBackgroundImageFilter.h"
  21 #include "clitkSliceBySliceRelativePositionFilter.h"
  22 #include "clitkCropLikeImageFilter.h"
  23 #include "clitkMemoryUsage.h"
  24
  25 // itk
  26 #include <itkConnectedComponentImageFilter.h>
  27 #include <itkRelabelComponentImageFilter.h>
  28 #include <itkBinaryThresholdImageFilter.h>
  29 #include <itkPasteImageFilter.h>
  30 #include <itkStatisticsLabelMapFilter.h>
  31 #include <itkBinaryBallStructuringElement.h>
  32 #include <itkBinaryDilateImageFilter.h>
  33 #include <itkConstantPadImageFilter.h>
  34 #include <itkImageSliceIteratorWithIndex.h>
  35 #include <itkBinaryMorphologicalOpeningImageFilter.h>
  36 #include <itkImageDuplicator.h>
  37
  38 namespace clitk {
  39
  40   //--------------------------------------------------------------------
  41   template<class ImageType>
  42   void ComputeBBFromImageRegion(const ImageType * image,
  43                                 typename ImageType::RegionType region,
  44                                 typename itk::BoundingBox<unsigned long,
  45                                 ImageType::ImageDimension>::Pointer bb) {
  46     typedef typename ImageType::IndexType IndexType;
  47     IndexType firstIndex;
  48     IndexType lastIndex;
  49     for(unsigned int i=0; i<image->GetImageDimension(); i++) {
  50       firstIndex[i] = region.GetIndex()[i];
  51       lastIndex[i] = firstIndex[i]+region.GetSize()[i];
  52     }
  53
  54     typedef itk::BoundingBox<unsigned long,
  55                              ImageType::ImageDimension> BBType;
  56     typedef typename BBType::PointType PointType;
  57     PointType lastPoint;
  58     PointType firstPoint;
  59     image->TransformIndexToPhysicalPoint(firstIndex, firstPoint);
  60     image->TransformIndexToPhysicalPoint(lastIndex, lastPoint);
  61
  62     bb->SetMaximum(lastPoint);
  63     bb->SetMinimum(firstPoint);
  64   }
  65   //--------------------------------------------------------------------
  66
  67
  68   //--------------------------------------------------------------------
  69   template<int Dimension>
  70   void ComputeBBIntersection(typename itk::BoundingBox<unsigned long, Dimension>::Pointer bbo,
  71                              typename itk::BoundingBox<unsigned long, Dimension>::Pointer bbi1,
  72                              typename itk::BoundingBox<unsigned long, Dimension>::Pointer bbi2) {
  73
  74     typedef itk::BoundingBox<unsigned long, Dimension> BBType;
  75     typedef typename BBType::PointType PointType;
  76     PointType lastPoint;
  77     PointType firstPoint;
  78
  79     for(unsigned int i=0; i<Dimension; i++) {
  80       firstPoint[i] = std::max(bbi1->GetMinimum()[i],
  81                                bbi2->GetMinimum()[i]);
  82       lastPoint[i] = std::min(bbi1->GetMaximum()[i],
  83                               bbi2->GetMaximum()[i]);
  84     }
  85
  86     bbo->SetMaximum(lastPoint);
  87     bbo->SetMinimum(firstPoint);
  88   }
  89   //--------------------------------------------------------------------
  90
  91
  92   //--------------------------------------------------------------------
  93   template<class ImageType>
  94   void ComputeRegionFromBB(const ImageType * image,
  95                            const typename itk::BoundingBox<unsigned long,
  96                                                            ImageType::ImageDimension>::Pointer bb,
  97                            typename ImageType::RegionType & region) {
  98     // Types
  99     typedef typename ImageType::IndexType  IndexType;
 100     typedef typename ImageType::PointType  PointType;
 101     typedef typename ImageType::RegionType RegionType;
 102     typedef typename ImageType::SizeType   SizeType;
 103
 104     // Region starting point
 105     IndexType regionStart;
 106     PointType start = bb->GetMinimum();
 107     image->TransformPhysicalPointToIndex(start, regionStart);
 108
 109     // Region size
 110     SizeType regionSize;
 111     PointType maxs = bb->GetMaximum();
 112     PointType mins = bb->GetMinimum();
 113     for(unsigned int i=0; i<ImageType::ImageDimension; i++) {
 114       regionSize[i] = lrint((maxs[i] - mins[i])/image->GetSpacing()[i]);
 115     }
 116
 117     // Create region
 118     region.SetIndex(regionStart);
 119     region.SetSize(regionSize);
 120   }
 121   //--------------------------------------------------------------------
 122
 123   //--------------------------------------------------------------------
 124   template<class ImageType, class TMaskImageType>
 125   typename ImageType::Pointer
 126   SetBackground(const ImageType * input,
 127                 const TMaskImageType * mask,
 128                 typename TMaskImageType::PixelType maskBG,
 129                 typename ImageType::PixelType outValue,
 130                 bool inPlace) {
 131     typedef SetBackgroundImageFilter<ImageType, TMaskImageType, ImageType>
 132       SetBackgroundImageFilterType;
 133     typename SetBackgroundImageFilterType::Pointer setBackgroundFilter
 134       = SetBackgroundImageFilterType::New();
 135     //  if (inPlace) setBackgroundFilter->ReleaseDataFlagOn(); // No seg fault
 136     setBackgroundFilter->SetInPlace(inPlace); // This is important to keep memory low
 137     setBackgroundFilter->SetInput(input);
 138     setBackgroundFilter->SetInput2(mask);
 139     setBackgroundFilter->SetMaskValue(maskBG);
 140     setBackgroundFilter->SetOutsideValue(outValue);
 141     setBackgroundFilter->Update();
 142     return setBackgroundFilter->GetOutput();
 143   }
 144   //--------------------------------------------------------------------
 145
 146
 147   //--------------------------------------------------------------------
 148   template<class ImageType>
 149   int GetNumberOfConnectedComponentLabels(const ImageType * input,
 150                                           typename ImageType::PixelType BG,
 151                                           bool isFullyConnected) {
 152     // Connected Component label
 153     typedef itk::ConnectedComponentImageFilter<ImageType, ImageType> ConnectFilterType;
 154     typename ConnectFilterType::Pointer connectFilter = ConnectFilterType::New();
 155     connectFilter->SetInput(input);
 156     connectFilter->SetBackgroundValue(BG);
 157     connectFilter->SetFullyConnected(isFullyConnected);
 158     connectFilter->Update();
 159
 160     // Return result
 161     return connectFilter->GetObjectCount();
 162   }
 163   //--------------------------------------------------------------------
 164
 165   //--------------------------------------------------------------------
 166   /*
 167     Warning : in this cas, we consider outputType like inputType, not
 168     InternalImageType. Be sure it fits.
 169   */
 170   template<class ImageType>
 171   typename ImageType::Pointer
 172   Labelize(const ImageType * input,
 173            typename ImageType::PixelType BG,
 174            bool isFullyConnected,
 175            int minimalComponentSize) {
 176     // InternalImageType for storing large number of component
 177     typedef itk::Image<int, ImageType::ImageDimension> InternalImageType;
 178
 179     // Connected Component label
 180     typedef itk::ConnectedComponentImageFilter<ImageType, InternalImageType> ConnectFilterType;
 181     typename ConnectFilterType::Pointer connectFilter = ConnectFilterType::New();
 182     //  connectFilter->ReleaseDataFlagOn();
 183     connectFilter->SetInput(input);
 184     connectFilter->SetBackgroundValue(BG);
 185     connectFilter->SetFullyConnected(isFullyConnected);
 186
 187     // Sort by size and remove too small area.
 188     typedef itk::RelabelComponentImageFilter<InternalImageType, ImageType> RelabelFilterType;
 189     typename RelabelFilterType::Pointer relabelFilter = RelabelFilterType::New();
 190     //  relabelFilter->ReleaseDataFlagOn(); // if yes, fail when ExplosionControlledThresholdConnectedImageFilter ???
 191     relabelFilter->SetInput(connectFilter->GetOutput());
 192     relabelFilter->SetMinimumObjectSize(minimalComponentSize);
 193     relabelFilter->Update();
 194
 195     // Return result
 196     typename ImageType::Pointer output = relabelFilter->GetOutput();
 197     return output;
 198   }
 199   //--------------------------------------------------------------------
 200
 201
 202   //--------------------------------------------------------------------
 203   /*
 204     Warning : in this cas, we consider outputType like inputType, not
 205     InternalImageType. Be sure it fits.
 206   */
 207   template<class ImageType>
 208   typename ImageType::Pointer
 209   LabelizeAndCountNumberOfObjects(const ImageType * input,
 210                                   typename ImageType::PixelType BG,
 211                                   bool isFullyConnected,
 212                                   int minimalComponentSize,
 213                                   int & nb) {
 214     // InternalImageType for storing large number of component
 215     typedef itk::Image<int, ImageType::ImageDimension> InternalImageType;
 216
 217     // Connected Component label
 218     typedef itk::ConnectedComponentImageFilter<ImageType, InternalImageType> ConnectFilterType;
 219     typename ConnectFilterType::Pointer connectFilter = ConnectFilterType::New();
 220     //  connectFilter->ReleaseDataFlagOn();
 221     connectFilter->SetInput(input);
 222     connectFilter->SetBackgroundValue(BG);
 223     connectFilter->SetFullyConnected(isFullyConnected);
 224
 225     // Sort by size and remove too small area.
 226     typedef itk::RelabelComponentImageFilter<InternalImageType, ImageType> RelabelFilterType;
 227     typename RelabelFilterType::Pointer relabelFilter = RelabelFilterType::New();
 228     //  relabelFilter->ReleaseDataFlagOn(); // if yes, fail when ExplosionControlledThresholdConnectedImageFilter ???
 229     relabelFilter->SetInput(connectFilter->GetOutput());
 230     relabelFilter->SetMinimumObjectSize(minimalComponentSize);
 231     relabelFilter->Update();
 232
 233     nb = relabelFilter->GetNumberOfObjects();
 234     // DD(relabelFilter->GetOriginalNumberOfObjects());
 235     // DD(relabelFilter->GetSizeOfObjectsInPhysicalUnits()[0]);
 236
 237     // Return result
 238     typename ImageType::Pointer output = relabelFilter->GetOutput();
 239     return output;
 240   }
 241   //--------------------------------------------------------------------
 242
 243
 244
 245   //--------------------------------------------------------------------
 246   template<class ImageType>
 247   typename ImageType::Pointer
 248   RemoveLabels(const ImageType * input,
 249                typename ImageType::PixelType BG,
 250                std::vector<typename ImageType::PixelType> & labelsToRemove) {
 251     assert(labelsToRemove.size() != 0);
 252     typename ImageType::Pointer working_image;// = input;
 253     for (unsigned int i=0; i <labelsToRemove.size(); i++) {
 254       typedef SetBackgroundImageFilter<ImageType, ImageType> SetBackgroundImageFilterType;
 255       typename SetBackgroundImageFilterType::Pointer setBackgroundFilter = SetBackgroundImageFilterType::New();
 256       setBackgroundFilter->SetInput(input);
 257       setBackgroundFilter->SetInput2(input);
 258       setBackgroundFilter->SetMaskValue(labelsToRemove[i]);
 259       setBackgroundFilter->SetOutsideValue(BG);
 260       setBackgroundFilter->Update();
 261       working_image = setBackgroundFilter->GetOutput();
 262     }
 263     return working_image;
 264   }
 265   //--------------------------------------------------------------------
 266
 267
 268   //--------------------------------------------------------------------
 269   template<class ImageType>
 270   typename ImageType::Pointer
 271   KeepLabels(const ImageType * input,
 272              typename ImageType::PixelType BG,
 273              typename ImageType::PixelType FG,
 274              typename ImageType::PixelType firstKeep,
 275              typename ImageType::PixelType lastKeep,
 276              bool useLastKeep) {
 277     typedef itk::BinaryThresholdImageFilter<ImageType, ImageType> BinarizeFilterType;
 278     typename BinarizeFilterType::Pointer binarizeFilter = BinarizeFilterType::New();
 279     binarizeFilter->SetInput(input);
 280     binarizeFilter->SetLowerThreshold(firstKeep);
 281     if (useLastKeep) binarizeFilter->SetUpperThreshold(lastKeep);
 282     binarizeFilter->SetInsideValue(FG);
 283     binarizeFilter->SetOutsideValue(BG);
 284     binarizeFilter->Update();
 285     return binarizeFilter->GetOutput();
 286   }
 287   //--------------------------------------------------------------------
 288
 289
 290   //--------------------------------------------------------------------
 291   template<class ImageType>
 292   typename ImageType::Pointer
 293   LabelizeAndSelectLabels(const ImageType * input,
 294                           typename ImageType::PixelType BG,
 295                           typename ImageType::PixelType FG,
 296                           bool isFullyConnected,
 297                           int minimalComponentSize,
 298                           LabelizeParameters<typename ImageType::PixelType> * param)
 299   {
 300     typename ImageType::Pointer working_image;
 301     working_image = Labelize<ImageType>(input, BG, isFullyConnected, minimalComponentSize);
 302     if (param->GetLabelsToRemove().size() != 0)
 303       working_image = RemoveLabels<ImageType>(working_image, BG, param->GetLabelsToRemove());
 304     working_image = KeepLabels<ImageType>(working_image,
 305                                           BG, FG,
 306                                           param->GetFirstKeep(),
 307                                           param->GetLastKeep(),
 308                                           param->GetUseLastKeep());
 309     return working_image;
 310   }
 311   //--------------------------------------------------------------------
 312
 313
 314   //--------------------------------------------------------------------
 315   template<class ImageType>
 316   typename ImageType::Pointer
 317   ResizeImageLike(const ImageType * input,
 318                   const itk::ImageBase<ImageType::ImageDimension> * like,
 319                   typename ImageType::PixelType backgroundValue)
 320   {
 321     typedef CropLikeImageFilter<ImageType> CropFilterType;
 322     typename CropFilterType::Pointer cropFilter = CropFilterType::New();
 323     cropFilter->SetInput(input);
 324     cropFilter->SetCropLikeImage(like);
 325     cropFilter->SetBackgroundValue(backgroundValue);
 326     cropFilter->Update();
 327     return cropFilter->GetOutput();
 328   }
 329   //--------------------------------------------------------------------
 330
 331
 332   //--------------------------------------------------------------------
 333   template<class MaskImageType>
 334   typename MaskImageType::Pointer
 335   SliceBySliceRelativePosition(const MaskImageType * input,
 336                                const MaskImageType * object,
 337                                int direction,
 338                                double threshold,
 339                                std::string orientation,
 340                                bool uniqueConnectedComponent,
 341                                double spacing,
 342                                bool autocropFlag,
 343                                bool singleObjectCCL)
 344   {
 345     typedef SliceBySliceRelativePositionFilter<MaskImageType> SliceRelPosFilterType;
 346     typename SliceRelPosFilterType::Pointer sliceRelPosFilter = SliceRelPosFilterType::New();
 347     sliceRelPosFilter->VerboseStepFlagOff();
 348     sliceRelPosFilter->WriteStepFlagOff();
 349     sliceRelPosFilter->SetInput(input);
 350     sliceRelPosFilter->SetInputObject(object);
 351     sliceRelPosFilter->SetDirection(direction);
 352     sliceRelPosFilter->SetFuzzyThreshold(threshold);
 353     sliceRelPosFilter->AddOrientationTypeString(orientation);
 354     sliceRelPosFilter->SetIntermediateSpacingFlag((spacing != -1));
 355     sliceRelPosFilter->SetIntermediateSpacing(spacing);
 356     sliceRelPosFilter->SetUniqueConnectedComponentBySliceFlag(uniqueConnectedComponent);
 357     sliceRelPosFilter->ObjectCCLSelectionFlagOff();
 358     sliceRelPosFilter->SetUseTheLargestObjectCCLFlag(singleObjectCCL);
 359     //    sliceRelPosFilter->SetInverseOrientationFlag(inverseflag);
 360     sliceRelPosFilter->SetAutoCropFlag(autocropFlag);
 361     sliceRelPosFilter->IgnoreEmptySliceObjectFlagOn();
 362     sliceRelPosFilter->Update();
 363     return sliceRelPosFilter->GetOutput();
 364   }
 365   //--------------------------------------------------------------------
 366
 367
 368   //--------------------------------------------------------------------
 369   template<class ImageType>
 370   bool
 371   FindExtremaPointInAGivenDirection(const ImageType * input,
 372                                     typename ImageType::PixelType bg,
 373                                     int direction, bool opposite,
 374                                     typename ImageType::PointType & point)
 375   {
 376     typename ImageType::PointType dummy;
 377     return FindExtremaPointInAGivenDirection(input, bg, direction,
 378                                              opposite, dummy, 0, point);
 379   }
 380   //--------------------------------------------------------------------
 381
 382
 383   //--------------------------------------------------------------------
 384   template<class ImageType>
 385   bool
 386   FindExtremaPointInAGivenDirection(const ImageType * input,
 387                                     typename ImageType::PixelType bg,
 388                                     int direction, bool opposite,
 389                                     typename ImageType::PointType refpoint,
 390                                     double distanceMax,
 391                                     typename ImageType::PointType & point)
 392   {
 393     /*
 394       loop over input pixels, store the index in the fg that is max
 395       according to the given direction.
 396     */
 397     typedef itk::ImageRegionConstIteratorWithIndex<ImageType> IteratorType;
 398     IteratorType iter(input, input->GetLargestPossibleRegion());
 399     iter.GoToBegin();
 400     typename ImageType::IndexType max = input->GetLargestPossibleRegion().GetIndex();
 401     if (opposite) max = max+input->GetLargestPossibleRegion().GetSize();
 402     bool found=false;
 403     while (!iter.IsAtEnd()) {
 404       if (iter.Get() != bg) {
 405         bool test = iter.GetIndex()[direction] >  max[direction];
 406         if (opposite) test = !test;
 407         if (test) {
 408           typename ImageType::PointType p;
 409           input->TransformIndexToPhysicalPoint(iter.GetIndex(), p);
 410           if ((distanceMax==0) || (p.EuclideanDistanceTo(refpoint) < distanceMax)) {
 411             max = iter.GetIndex();
 412             found = true;
 413           }
 414         }
 415       }
 416       ++iter;
 417     }
 418     if (!found) return false;
 419     input->TransformIndexToPhysicalPoint(max, point);
 420     return true;
 421   }
 422   //--------------------------------------------------------------------
 423
 424
 425   //--------------------------------------------------------------------
 426   template<class ImageType>
 427   typename ImageType::Pointer
 428   CropImageRemoveGreaterThan(const ImageType * image,
 429                  int dim, double min, bool autoCrop,
 430                  typename ImageType::PixelType BG)
 431   {
 432     return CropImageAlongOneAxis<ImageType>(image, dim,
 433                                             image->GetOrigin()[dim],
 434                                             min,
 435                                             autoCrop, BG);
 436   }
 437   //--------------------------------------------------------------------
 438
 439
 440   //--------------------------------------------------------------------
 441   template<class ImageType>
 442   typename ImageType::Pointer
 443   CropImageRemoveLowerThan(const ImageType * image,
 444                  int dim, double max, bool autoCrop,
 445                  typename ImageType::PixelType BG)
 446   {
 447     typename ImageType::PointType p;
 448     image->TransformIndexToPhysicalPoint(image->GetLargestPossibleRegion().GetIndex()+
 449                                          image->GetLargestPossibleRegion().GetSize(), p);
 450     return CropImageAlongOneAxis<ImageType>(image, dim, max, p[dim], autoCrop, BG);
 451   }
 452   //--------------------------------------------------------------------
 453
 454
 455   //--------------------------------------------------------------------
 456   template<class ImageType>
 457   typename ImageType::Pointer
 458   CropImageAlongOneAxis(const ImageType * image,
 459                         int dim, double min, double max,
 460                         bool autoCrop, typename ImageType::PixelType BG)
 461   {
 462     // Compute region size
 463     typename ImageType::RegionType region;
 464     typename ImageType::SizeType size = image->GetLargestPossibleRegion().GetSize();
 465     typename ImageType::PointType p = image->GetOrigin();
 466     p[dim] = min;
 467     typename ImageType::IndexType start;
 468     image->TransformPhysicalPointToIndex(p, start);
 469     p[dim] = max;
 470     typename ImageType::IndexType end;
 471     image->TransformPhysicalPointToIndex(p, end);
 472     size[dim] = abs(end[dim]-start[dim]);
 473     region.SetIndex(start);
 474     region.SetSize(size);
 475
 476     // Perform Crop
 477     typedef itk::RegionOfInterestImageFilter<ImageType, ImageType> CropFilterType;
 478     typename CropFilterType::Pointer cropFilter = CropFilterType::New();
 479     cropFilter->SetInput(image);
 480     cropFilter->SetRegionOfInterest(region);
 481     cropFilter->Update();
 482     typename ImageType::Pointer result = cropFilter->GetOutput();
 483
 484     // Auto Crop
 485     if (autoCrop) {
 486       result = AutoCrop<ImageType>(result, BG);
 487     }
 488     return result;
 489   }
 490   //--------------------------------------------------------------------
 491
 492
 493   //--------------------------------------------------------------------
 494   template<class ImageType>
 495   void
 496   ComputeCentroids(const ImageType * image,
 497                    typename ImageType::PixelType BG,
 498                    std::vector<typename ImageType::PointType> & centroids)
 499   {
 500     typedef long LabelType;
 501     static const unsigned int Dim = ImageType::ImageDimension;
 502     typedef itk::ShapeLabelObject< LabelType, Dim > LabelObjectType;
 503     typedef itk::LabelMap< LabelObjectType > LabelMapType;
 504     typedef itk::LabelImageToLabelMapFilter<ImageType, LabelMapType> ImageToMapFilterType;
 505     typename ImageToMapFilterType::Pointer imageToLabelFilter = ImageToMapFilterType::New();
 506     typedef itk::ShapeLabelMapFilter<LabelMapType, ImageType> ShapeFilterType;
 507     typename ShapeFilterType::Pointer statFilter = ShapeFilterType::New();
 508     imageToLabelFilter->SetBackgroundValue(BG);
 509     imageToLabelFilter->SetInput(image);
 510     statFilter->SetInput(imageToLabelFilter->GetOutput());
 511     statFilter->Update();
 512     typename LabelMapType::Pointer labelMap = statFilter->GetOutput();
 513
 514     centroids.clear();
 515     typename ImageType::PointType dummy;
 516     centroids.push_back(dummy); // label 0 -> no centroid, use dummy point for BG
 517     //DS FIXME (not useful ! to change ..)
 518     for(uint i=0; i<labelMap->GetNumberOfLabelObjects(); i++) {
 519       int label = labelMap->GetLabels()[i];
 520       centroids.push_back(labelMap->GetLabelObject(label)->GetCentroid());
 521     }
 522   }
 523   //--------------------------------------------------------------------
 524
 525
 526   //--------------------------------------------------------------------
 527   template<class ImageType, class LabelType>
 528   typename itk::LabelMap< itk::ShapeLabelObject<LabelType, ImageType::ImageDimension> >::Pointer
 529   ComputeLabelMap(const ImageType * image,
 530                   typename ImageType::PixelType BG,
 531                   bool computePerimeterFlag)
 532   {
 533     static const unsigned int Dim = ImageType::ImageDimension;
 534     typedef itk::ShapeLabelObject< LabelType, Dim > LabelObjectType;
 535     typedef itk::LabelMap< LabelObjectType > LabelMapType;
 536     typedef itk::LabelImageToLabelMapFilter<ImageType, LabelMapType> ImageToMapFilterType;
 537     typename ImageToMapFilterType::Pointer imageToLabelFilter = ImageToMapFilterType::New();
 538     typedef itk::ShapeLabelMapFilter<LabelMapType, ImageType> ShapeFilterType;
 539     typename ShapeFilterType::Pointer statFilter = ShapeFilterType::New();
 540     imageToLabelFilter->SetBackgroundValue(BG);
 541     imageToLabelFilter->SetInput(image);
 542     statFilter->SetInput(imageToLabelFilter->GetOutput());
 543     statFilter->SetComputePerimeter(computePerimeterFlag);
 544     statFilter->Update();
 545     return statFilter->GetOutput();
 546   }
 547   //--------------------------------------------------------------------
 548
 549
 550   //--------------------------------------------------------------------
 551   template<class ImageType>
 552   void
 553   ComputeCentroids2(const ImageType * image,
 554                    typename ImageType::PixelType BG,
 555                    std::vector<typename ImageType::PointType> & centroids)
 556   {
 557     typedef long LabelType;
 558     static const unsigned int Dim = ImageType::ImageDimension;
 559     typedef itk::ShapeLabelObject< LabelType, Dim > LabelObjectType;
 560     typedef itk::LabelMap< LabelObjectType > LabelMapType;
 561     typedef itk::LabelImageToLabelMapFilter<ImageType, LabelMapType> ImageToMapFilterType;
 562     typename ImageToMapFilterType::Pointer imageToLabelFilter = ImageToMapFilterType::New();
 563     typedef itk::ShapeLabelMapFilter<LabelMapType, ImageType> ShapeFilterType;
 564     typename ShapeFilterType::Pointer statFilter = ShapeFilterType::New();
 565     imageToLabelFilter->SetBackgroundValue(BG);
 566     imageToLabelFilter->SetInput(image);
 567     statFilter->SetInput(imageToLabelFilter->GetOutput());
 568     statFilter->Update();
 569     typename LabelMapType::Pointer labelMap = statFilter->GetOutput();
 570
 571     centroids.clear();
 572     typename ImageType::PointType dummy;
 573     centroids.push_back(dummy); // label 0 -> no centroid, use dummy point
 574     for(uint i=1; i<labelMap->GetNumberOfLabelObjects()+1; i++) {
 575       centroids.push_back(labelMap->GetLabelObject(i)->GetCentroid());
 576     }
 577
 578     for(uint i=1; i<labelMap->GetNumberOfLabelObjects()+1; i++) {
 579       DD(labelMap->GetLabelObject(i)->GetBinaryPrincipalAxes());
 580       DD(labelMap->GetLabelObject(i)->GetBinaryFlatness());
 581       DD(labelMap->GetLabelObject(i)->GetRoundness ());
 582
 583       // search for the point on the boundary alog PA
 584
 585     }
 586
 587   }
 588   //--------------------------------------------------------------------
 589
 590
 591   //--------------------------------------------------------------------
 592   template<class ImageType>
 593   void
 594   PointsUtils<ImageType>::Convert2DTo3D(const PointType2D & p2D,
 595                                         const ImageType * image,
 596                                         const int slice,
 597                                         PointType3D & p3D)
 598   {
 599     IndexType3D index3D;
 600     index3D[0] = index3D[1] = 0;
 601     index3D[2] = image->GetLargestPossibleRegion().GetIndex()[2]+slice;
 602     image->TransformIndexToPhysicalPoint(index3D, p3D);
 603     p3D[0] = p2D[0];
 604     p3D[1] = p2D[1];
 605     //  p3D[2] = p[2];//(image->GetLargestPossibleRegion().GetIndex()[2]+slice)*image->GetSpacing()[2]
 606     //    + image->GetOrigin()[2];
 607   }
 608   //--------------------------------------------------------------------
 609
 610
 611   //--------------------------------------------------------------------
 612   template<class ImageType>
 613   void
 614   PointsUtils<ImageType>::Convert2DMapTo3DList(const MapPoint2DType & map,
 615                                             const ImageType * image,
 616                                             VectorPoint3DType & list)
 617   {
 618     typename MapPoint2DType::const_iterator iter = map.begin();
 619     while (iter != map.end()) {
 620       PointType3D p;
 621       Convert2DTo3D(iter->second, image, iter->first, p);
 622       list.push_back(p);
 623       ++iter;
 624     }
 625   }
 626   //--------------------------------------------------------------------
 627
 628
 629   //--------------------------------------------------------------------
 630   template<class ImageType>
 631   void
 632   PointsUtils<ImageType>::Convert2DListTo3DList(const VectorPoint2DType & p2D,
 633                                                 int slice,
 634                                                 const ImageType * image,
 635                                                 VectorPoint3DType & list)
 636   {
 637     for(uint i=0; i<p2D.size(); i++) {
 638       PointType3D p;
 639       Convert2DTo3D(p2D[i], image, slice, p);
 640       list.push_back(p);
 641     }
 642   }
 643   //--------------------------------------------------------------------
 644
 645
 646   //--------------------------------------------------------------------
 647   template<class ImageType>
 648   void
 649   WriteListOfLandmarks(std::vector<typename ImageType::PointType> points,
 650                        std::string filename)
 651   {
 652     std::ofstream os;
 653     openFileForWriting(os, filename);
 654     os << "LANDMARKS1" << std::endl;
 655     for(uint i=0; i<points.size(); i++) {
 656       const typename ImageType::PointType & p = points[i];
 657       // Write it in the file
 658       os << i << " " << p[0] << " " << p[1] << " " << p[2] << " 0 0 " << std::endl;
 659     }
 660     os.close();
 661   }
 662   //--------------------------------------------------------------------
 663
 664
 665   //--------------------------------------------------------------------
 666   template<class ImageType>
 667   typename ImageType::Pointer
 668   Dilate(const ImageType * image, double radiusInMM,
 669          typename ImageType::PixelType BG,
 670          typename ImageType::PixelType FG,
 671          bool extendSupport)
 672   {
 673     typename ImageType::SizeType r;
 674     for(uint i=0; i<ImageType::ImageDimension; i++)
 675       r[i] = (uint)lrint(radiusInMM/image->GetSpacing()[i]);
 676     return Dilate<ImageType>(image, r, BG, FG, extendSupport);
 677   }
 678   //--------------------------------------------------------------------
 679
 680
 681   //--------------------------------------------------------------------
 682   template<class ImageType>
 683   typename ImageType::Pointer
 684   Dilate(const ImageType * image, typename ImageType::PointType radiusInMM,
 685          typename ImageType::PixelType BG,
 686          typename ImageType::PixelType FG,
 687          bool extendSupport)
 688   {
 689     typename ImageType::SizeType r;
 690     for(uint i=0; i<ImageType::ImageDimension; i++)
 691       r[i] = (uint)lrint(radiusInMM[i]/image->GetSpacing()[i]);
 692     return Dilate<ImageType>(image, r, BG, FG, extendSupport);
 693   }
 694   //--------------------------------------------------------------------
 695
 696
 697   //--------------------------------------------------------------------
 698   template<class ImageType>
 699   typename ImageType::Pointer
 700   Dilate(const ImageType * image, typename ImageType::SizeType radius,
 701          typename ImageType::PixelType BG,
 702          typename ImageType::PixelType FG,
 703          bool extendSupport)
 704   {
 705     // Create kernel for dilatation
 706     typedef itk::BinaryBallStructuringElement<typename ImageType::PixelType,
 707                                               ImageType::ImageDimension> KernelType;
 708     KernelType structuringElement;
 709     structuringElement.SetRadius(radius);
 710     structuringElement.CreateStructuringElement();
 711
 712     typename ImageType::Pointer output;
 713     if (extendSupport) {
 714       typedef itk::ConstantPadImageFilter<ImageType, ImageType> PadFilterType;
 715       typename PadFilterType::Pointer padFilter = PadFilterType::New();
 716       padFilter->SetInput(image);
 717       typename ImageType::SizeType lower;
 718       typename ImageType::SizeType upper;
 719       for(uint i=0; i<3; i++) {
 720         lower[i] = upper[i] = 2*(radius[i]+1);
 721       }
 722       padFilter->SetPadLowerBound(lower);
 723       padFilter->SetPadUpperBound(upper);
 724       padFilter->Update();
 725       output = padFilter->GetOutput();
 726     }
 727
 728     // Dilate  filter
 729     typedef itk::BinaryDilateImageFilter<ImageType, ImageType , KernelType> DilateFilterType;
 730     typename DilateFilterType::Pointer dilateFilter = DilateFilterType::New();
 731     dilateFilter->SetBackgroundValue(BG);
 732     dilateFilter->SetForegroundValue(FG);
 733     dilateFilter->SetBoundaryToForeground(false);
 734     dilateFilter->SetKernel(structuringElement);
 735     if (extendSupport) dilateFilter->SetInput(output);
 736     else dilateFilter->SetInput(image);
 737     dilateFilter->Update();
 738     return dilateFilter->GetOutput();
 739   }
 740   //--------------------------------------------------------------------
 741
 742
 743   //--------------------------------------------------------------------
 744   template<class ImageType>
 745   typename ImageType::Pointer
 746   Opening(const ImageType * image, typename ImageType::SizeType radius,
 747          typename ImageType::PixelType BG,
 748          typename ImageType::PixelType FG)
 749   {
 750     // Kernel
 751     typedef itk::BinaryBallStructuringElement<typename ImageType::PixelType,
 752                                               ImageType::ImageDimension> KernelType;
 753     KernelType structuringElement;
 754     structuringElement.SetRadius(radius);
 755     structuringElement.CreateStructuringElement();
 756
 757     // Filter
 758     typedef itk::BinaryMorphologicalOpeningImageFilter<ImageType, ImageType , KernelType> OpeningFilterType;
 759     typename OpeningFilterType::Pointer open = OpeningFilterType::New();
 760     open->SetInput(image);
 761     open->SetBackgroundValue(BG);
 762     open->SetForegroundValue(FG);
 763     open->SetKernel(structuringElement);
 764     open->Update();
 765     return open->GetOutput();
 766   }
 767   //--------------------------------------------------------------------
 768
 769
 770
 771   //--------------------------------------------------------------------
 772   template<class ValueType, class VectorType>
 773   void ConvertOption(std::string optionName, uint given,
 774                      ValueType * values, VectorType & p,
 775                      uint dim, bool required)
 776   {
 777     if (required && (given == 0)) {
 778       clitkExceptionMacro("The option --" << optionName << " must be set and have 1 or "
 779                           << dim << " values.");
 780     }
 781     if (given == 1) {
 782       for(uint i=0; i<dim; i++) p[i] = values[0];
 783       return;
 784     }
 785     if (given == dim) {
 786       for(uint i=0; i<dim; i++) p[i] = values[i];
 787       return;
 788     }
 789     if (given == 0) return;
 790     clitkExceptionMacro("The option --" << optionName << " must have 1 or "
 791                         << dim << " values.");
 792   }
 793   //--------------------------------------------------------------------
 794
 795
 796   //--------------------------------------------------------------------
 797   /*
 798     http://www.gamedev.net/community/forums/topic.asp?topic_id=542870
 799     Assuming the points are (Ax,Ay) (Bx,By) and (Cx,Cy), you need to compute:
 800     (Bx - Ax) * (Cy - Ay) - (By - Ay) * (Cx - Ax)
 801     This will equal zero if the point C is on the line formed by
 802     points A and B, and will have a different sign depending on the
 803     side. Which side this is depends on the orientation of your (x,y)
 804     coordinates, but you can plug test values for A,B and C into this
 805     formula to determine whether negative values are to the left or to
 806     the right.
 807     => to accelerate, start with formula, when change sign -> stop and fill
 808
 809     offsetToKeep = is used to determine which side of the line we
 810     keep. The point along the mainDirection but 'offsetToKeep' mm away
 811     is kept.
 812
 813   */
 814   template<class ImageType>
 815   void
 816   SliceBySliceSetBackgroundFromLineSeparation(ImageType * input,
 817                                               std::vector<typename ImageType::PointType> & lA,
 818                                               std::vector<typename ImageType::PointType> & lB,
 819                                               typename ImageType::PixelType BG,
 820                                               int mainDirection,
 821                                               double offsetToKeep)
 822   {
 823     assert((mainDirection==0) || (mainDirection==1));
 824     typedef itk::ImageSliceIteratorWithIndex<ImageType> SliceIteratorType;
 825     SliceIteratorType siter = SliceIteratorType(input,
 826                                                 input->GetLargestPossibleRegion());
 827     siter.SetFirstDirection(0);
 828     siter.SetSecondDirection(1);
 829     siter.GoToBegin();
 830     uint i=0;
 831     typename ImageType::PointType A;
 832     typename ImageType::PointType B;
 833     typename ImageType::PointType C;
 834     assert(lA.size() == lB.size());
 835     while ((i<lA.size()) && (!siter.IsAtEnd())) {
 836       // Check that the current slice correspond to the current point
 837       input->TransformIndexToPhysicalPoint(siter.GetIndex(), C);
 838       if ((fabs(C[2] - lA[i][2]))>0.01) { // is !equal with a tolerance of 0.01 mm
 839       }
 840       else {
 841         // Define A,B,C points
 842         A = lA[i];
 843         B = lB[i];
 844         C = A;
 845
 846         // Check that the line is not a point (A=B)
 847         bool p = (A[0] == B[0]) && (A[1] == B[1]);
 848
 849         if (!p) {
 850           C[mainDirection] += offsetToKeep; // I know I must keep this point
 851           double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
 852           bool isPositive = s<0;
 853           while (!siter.IsAtEndOfSlice()) {
 854             while (!siter.IsAtEndOfLine()) {
 855               // Very slow, I know ... but image should be very small
 856               input->TransformIndexToPhysicalPoint(siter.GetIndex(), C);
 857               double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
 858               if (s == 0) siter.Set(BG); // on the line, we decide to remove
 859               if (isPositive) {
 860                 if (s > 0) siter.Set(BG);
 861               }
 862               else {
 863                 if (s < 0) siter.Set(BG);
 864               }
 865               ++siter;
 866             }
 867             siter.NextLine();
 868           } // end loop slice
 869         }
 870
 871         ++i;
 872       } // End of current slice
 873       siter.NextSlice();
 874     }
 875   }
 876   //--------------------------------------------------------------------
 877
 878
 879   //--------------------------------------------------------------------
 880   template<class ImageType>
 881   void
 882   AndNot(ImageType * input,
 883          const ImageType * object,
 884          typename ImageType::PixelType BG)
 885   {
 886     typename ImageType::Pointer o;
 887     bool resized=false;
 888     if (!clitk::HaveSameSizeAndSpacing<ImageType, ImageType>(input, object)) {
 889       o = clitk::ResizeImageLike<ImageType>(object, input, BG);
 890       resized = true;
 891     }
 892
 893     typedef clitk::BooleanOperatorLabelImageFilter<ImageType> BoolFilterType;
 894     typename BoolFilterType::Pointer boolFilter = BoolFilterType::New();
 895     boolFilter->InPlaceOn();
 896     boolFilter->SetInput1(input);
 897     if (resized) boolFilter->SetInput2(o);
 898     else boolFilter->SetInput2(object);
 899     boolFilter->SetBackgroundValue1(BG);
 900     boolFilter->SetBackgroundValue2(BG);
 901     boolFilter->SetOperationType(BoolFilterType::AndNot);
 902     boolFilter->Update();
 903   }
 904   //--------------------------------------------------------------------
 905
 906
 907   //--------------------------------------------------------------------
 908   template<class ImageType>
 909   void
 910   And(ImageType * input,
 911       const ImageType * object,
 912       typename ImageType::PixelType BG)
 913   {
 914     typename ImageType::Pointer o;
 915     bool resized=false;
 916     if (!clitk::HaveSameSizeAndSpacing<ImageType, ImageType>(input, object)) {
 917       o = clitk::ResizeImageLike<ImageType>(object, input, BG);
 918       resized = true;
 919     }
 920
 921     typedef clitk::BooleanOperatorLabelImageFilter<ImageType> BoolFilterType;
 922     typename BoolFilterType::Pointer boolFilter = BoolFilterType::New();
 923     boolFilter->InPlaceOn();
 924     boolFilter->SetInput1(input);
 925     if (resized) boolFilter->SetInput2(o);
 926     else boolFilter->SetInput2(object);
 927     boolFilter->SetBackgroundValue1(BG);
 928     boolFilter->SetBackgroundValue2(BG);
 929     boolFilter->SetOperationType(BoolFilterType::And);
 930     boolFilter->Update();
 931   }
 932   //--------------------------------------------------------------------
 933
 934
 935   //--------------------------------------------------------------------
 936   template<class ImageType>
 937   void
 938   Or(ImageType * input,
 939      const ImageType * object,
 940      typename ImageType::PixelType BG)
 941   {
 942     typename ImageType::Pointer o;
 943     bool resized=false;
 944     if (!clitk::HaveSameSizeAndSpacing<ImageType, ImageType>(input, object)) {
 945       o = clitk::ResizeImageLike<ImageType>(object, input, BG);
 946       resized = true;
 947     }
 948
 949     typedef clitk::BooleanOperatorLabelImageFilter<ImageType> BoolFilterType;
 950     typename BoolFilterType::Pointer boolFilter = BoolFilterType::New();
 951     boolFilter->InPlaceOn();
 952     boolFilter->SetInput1(input);
 953     if (resized) boolFilter->SetInput2(o);
 954     else boolFilter->SetInput2(object);
 955     boolFilter->SetBackgroundValue1(BG);
 956     boolFilter->SetBackgroundValue2(BG);
 957     boolFilter->SetOperationType(BoolFilterType::Or);
 958     boolFilter->Update();
 959   }
 960   //--------------------------------------------------------------------
 961
 962
 963   //--------------------------------------------------------------------
 964   template<class ImageType>
 965   typename ImageType::Pointer
 966   Binarize(const ImageType * input,
 967            typename ImageType::PixelType lower,
 968            typename ImageType::PixelType upper,
 969            typename ImageType::PixelType BG,
 970            typename ImageType::PixelType FG)
 971   {
 972     typedef itk::BinaryThresholdImageFilter<ImageType, ImageType> BinaryThresholdFilterType;
 973     typename BinaryThresholdFilterType::Pointer binarizeFilter = BinaryThresholdFilterType::New();
 974     binarizeFilter->SetInput(input);
 975     binarizeFilter->InPlaceOff();
 976     binarizeFilter->SetLowerThreshold(lower);
 977     binarizeFilter->SetUpperThreshold(upper);
 978     binarizeFilter->SetInsideValue(FG);
 979     binarizeFilter->SetOutsideValue(BG);
 980     binarizeFilter->Update();
 981     return binarizeFilter->GetOutput();
 982   }
 983   //--------------------------------------------------------------------
 984
 985
 986   //--------------------------------------------------------------------
 987   template<class ImageType>
 988   void
 989   GetMinMaxPointPosition(const ImageType * input,
 990                          typename ImageType::PointType & min,
 991                          typename ImageType::PointType & max)
 992   {
 993     typename ImageType::IndexType index = input->GetLargestPossibleRegion().GetIndex();
 994     input->TransformIndexToPhysicalPoint(index, min);
 995     index = index+input->GetLargestPossibleRegion().GetSize();
 996     input->TransformIndexToPhysicalPoint(index, max);
 997   }
 998   //--------------------------------------------------------------------
 999
1000
1001   //--------------------------------------------------------------------
1002   template<class ImageType>
1003   typename ImageType::PointType
1004   FindExtremaPointInAGivenLine(const ImageType * input,
1005                                int dimension,
1006                                bool inverse,
1007                                typename ImageType::PointType p,
1008                                typename ImageType::PixelType BG,
1009                                double distanceMax)
1010   {
1011     // Which direction ?  Increasing or decreasing.
1012     int d=1;
1013     if (inverse) d=-1;
1014
1015     // Transform to pixel index
1016     typename ImageType::IndexType index;
1017     input->TransformPhysicalPointToIndex(p, index);
1018
1019     // Loop while inside the mask;
1020     while (input->GetPixel(index) != BG) {
1021       index[dimension] += d;
1022     }
1023
1024     // Transform back to Physical Units
1025     typename ImageType::PointType result;
1026     input->TransformIndexToPhysicalPoint(index, result);
1027
1028     // Check that is is not too far away
1029     double distance = p.EuclideanDistanceTo(result);
1030     if (distance > distanceMax) {
1031       result = p; // Get back to initial value
1032     }
1033
1034     return result;
1035   }
1036   //--------------------------------------------------------------------
1037
1038
1039   //--------------------------------------------------------------------
1040   template<class PointType>
1041   bool
1042   IsOnTheSameLineSide(PointType C, PointType A, PointType B, PointType like)
1043   {
1044     // Look at the position of point 'like' according to the AB line
1045     double s = (B[0] - A[0]) * (like[1] - A[1]) - (B[1] - A[1]) * (like[0] - A[0]);
1046     bool negative = s<0;
1047
1048     // Look the C position
1049     s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
1050
1051     if (negative && (s<=0)) return true;
1052     if (!negative && (s>=0)) return true;
1053     return false;
1054   }
1055   //--------------------------------------------------------------------
1056
1057
1058   //--------------------------------------------------------------------
1059   /* Consider an input object, for each slice, find the extrema
1060      position according to a given direction and build a line segment
1061      passing throught this point in a given direction.  Output is a
1062      vector of line (from point A to B), for each slice;
1063    */
1064   template<class ImageType>
1065   void
1066   SliceBySliceBuildLineSegmentAccordingToExtremaPosition(const ImageType * input,
1067                                                          typename ImageType::PixelType BG,
1068                                                          int sliceDimension,
1069                                                          int extremaDirection,
1070                                                          bool extremaOppositeFlag,
1071                                                          int lineDirection,
1072                                                          double margin,
1073                                                          std::vector<typename ImageType::PointType> & A,
1074                                                          std::vector<typename ImageType::PointType> & B)
1075   {
1076     // Type of a slice
1077     typedef typename itk::Image<typename ImageType::PixelType, ImageType::ImageDimension-1> SliceType;
1078
1079     // Build the list of slices
1080     std::vector<typename SliceType::Pointer> slices;
1081     clitk::ExtractSlices<ImageType>(input, sliceDimension, slices);
1082
1083     // Build the list of 2D points
1084     std::map<int, typename SliceType::PointType> position2D;
1085     for(uint i=0; i<slices.size(); i++) {
1086       typename SliceType::PointType p;
1087       bool found =
1088         clitk::FindExtremaPointInAGivenDirection<SliceType>(slices[i], BG,
1089                                                             extremaDirection, extremaOppositeFlag, p);
1090       if (found) {
1091         position2D[i] = p;
1092       }
1093     }
1094
1095     // Convert 2D points in slice into 3D points
1096     clitk::PointsUtils<ImageType>::Convert2DMapTo3DList(position2D, input, A);
1097
1098     // Create additional point just right to the previous ones, on the
1099     // given lineDirection, in order to create a horizontal/vertical line.
1100     for(uint i=0; i<A.size(); i++) {
1101       typename ImageType::PointType p = A[i];
1102       p[lineDirection] += 10;
1103       B.push_back(p);
1104       // Margins ?
1105       A[i][extremaDirection] += margin;
1106       B[i][extremaDirection] += margin;
1107     }
1108
1109   }
1110   //--------------------------------------------------------------------
1111
1112
1113   //--------------------------------------------------------------------
1114   template<class ImageType>
1115   typename ImageType::Pointer
1116   SliceBySliceKeepMainCCL(const ImageType * input,
1117                           typename ImageType::PixelType BG,
1118                           typename ImageType::PixelType FG)  {
1119
1120     // Extract slices
1121     const int d = ImageType::ImageDimension-1;
1122     typedef typename itk::Image<typename ImageType::PixelType, d> SliceType;
1123     std::vector<typename SliceType::Pointer> slices;
1124     clitk::ExtractSlices<ImageType>(input, d, slices);
1125
1126     // Labelize and keep the main one
1127     std::vector<typename SliceType::Pointer> o;
1128     for(uint i=0; i<slices.size(); i++) {
1129       o.push_back(clitk::Labelize<SliceType>(slices[i], BG, false, 1));
1130       o[i] = clitk::KeepLabels<SliceType>(o[i], BG, FG, 1, 1, true);
1131     }
1132
1133     // Join slices
1134     typename ImageType::Pointer output;
1135     output = clitk::JoinSlices<ImageType>(o, input, d);
1136     return output;
1137   }
1138   //--------------------------------------------------------------------
1139
1140
1141   //--------------------------------------------------------------------
1142   template<class ImageType>
1143   typename ImageType::Pointer
1144   Clone(const ImageType * input) {
1145     typedef itk::ImageDuplicator<ImageType> DuplicatorType;
1146     typename DuplicatorType::Pointer duplicator = DuplicatorType::New();
1147     duplicator->SetInputImage(input);
1148     duplicator->Update();
1149     return duplicator->GetOutput();
1150   }
1151   //--------------------------------------------------------------------
1152
1153
1154   //--------------------------------------------------------------------
1155   /* Consider an input object, start at A, for each slice (dim1):
1156      - compute the intersection between the AB line and the current slice
1157      - remove what is at lower or greater according to dim2 of this point
1158      - stop at B
1159   */
1160   template<class ImageType>
1161   typename ImageType::Pointer
1162   SliceBySliceSetBackgroundFromSingleLine(const ImageType * input,
1163                                           typename ImageType::PixelType BG,
1164                                           typename ImageType::PointType & A,
1165                                           typename ImageType::PointType & B,
1166                                           int dim1, int dim2, bool removeLowerPartFlag)
1167
1168   {
1169     // Extract slices
1170     typedef typename itk::Image<typename ImageType::PixelType, ImageType::ImageDimension-1> SliceType;
1171     typedef typename SliceType::Pointer SlicePointer;
1172     std::vector<SlicePointer> slices;
1173     clitk::ExtractSlices<ImageType>(input, dim1, slices);
1174
1175     // Start at slice that contains A, and stop at B
1176     typename ImageType::IndexType Ap;
1177     typename ImageType::IndexType Bp;
1178     input->TransformPhysicalPointToIndex(A, Ap);
1179     input->TransformPhysicalPointToIndex(B, Bp);
1180
1181     // Determine slice largest region
1182     typename SliceType::RegionType region = slices[0]->GetLargestPossibleRegion();
1183     typename SliceType::SizeType size = region.GetSize();
1184     typename SliceType::IndexType index = region.GetIndex();
1185
1186     // Line slope
1187     double a = (Bp[dim2]-Ap[dim2])/(Bp[dim1]-Ap[dim1]);
1188     double b = Ap[dim2];
1189
1190     // Loop from slice A to slice B
1191     for(uint i=0; i<(Bp[dim1]-Ap[dim1]); i++) {
1192       // Compute intersection between line AB and current slice for the dim2
1193       double p = a*i+b;
1194       // Change region (lower than dim2)
1195       if (removeLowerPartFlag) {
1196         size[dim2] = p-Ap[dim2];
1197       }
1198       else {
1199         size[dim2] = slices[0]->GetLargestPossibleRegion().GetSize()[dim2]-p;
1200         index[dim2] = p;
1201       }
1202       region.SetSize(size);
1203       region.SetIndex(index);
1204       // Fill region with BG (simple region iterator)
1205       FillRegionWithValue<SliceType>(slices[i+Ap[dim1]], BG, region);
1206       /*
1207       typedef itk::ImageRegionIterator<SliceType> IteratorType;
1208       IteratorType iter(slices[i+Ap[dim1]], region);
1209       iter.GoToBegin();
1210       while (!iter.IsAtEnd()) {
1211         iter.Set(BG);
1212         ++iter;
1213       }
1214       */
1215       // Loop
1216     }
1217
1218     // Merge slices
1219     typename ImageType::Pointer output;
1220     output = clitk::JoinSlices<ImageType>(slices, input, dim1);
1221     return output;
1222   }
1223   //--------------------------------------------------------------------
1224
1225   //--------------------------------------------------------------------
1226   /* Consider an input object, slice by slice, use the point A and set
1227      pixel to BG according to their position relatively to A
1228   */
1229   template<class ImageType>
1230   typename ImageType::Pointer
1231   SliceBySliceSetBackgroundFromPoints(const ImageType * input,
1232                                       typename ImageType::PixelType BG,
1233                                       int sliceDim,
1234                                       std::vector<typename ImageType::PointType> & A,
1235                                       bool removeGreaterThanXFlag,
1236                                       bool removeGreaterThanYFlag)
1237
1238   {
1239     // Extract slices
1240     typedef typename itk::Image<typename ImageType::PixelType, ImageType::ImageDimension-1> SliceType;
1241     typedef typename SliceType::Pointer SlicePointer;
1242     std::vector<SlicePointer> slices;
1243     clitk::ExtractSlices<ImageType>(input, sliceDim, slices);
1244
1245     // Start at slice that contains A
1246     typename ImageType::IndexType Ap;
1247
1248     // Determine slice largest region
1249     typename SliceType::RegionType region = slices[0]->GetLargestPossibleRegion();
1250     typename SliceType::SizeType size = region.GetSize();
1251     typename SliceType::IndexType index = region.GetIndex();
1252
1253     // Loop from slice A to slice B
1254     for(uint i=0; i<A.size(); i++) {
1255       input->TransformPhysicalPointToIndex(A[i], Ap);
1256       uint sliceIndex = Ap[2] - input->GetLargestPossibleRegion().GetIndex()[2];
1257       if ((sliceIndex < 0) || (sliceIndex >= slices.size())) {
1258         continue; // do not consider this slice
1259       }
1260
1261       // Compute region for BG
1262       if (removeGreaterThanXFlag) {
1263         index[0] = Ap[0];
1264         size[0] = region.GetSize()[0]-(index[0]-region.GetIndex()[0]);
1265       }
1266       else {
1267         index[0] = region.GetIndex()[0];
1268         size[0] = Ap[0] - index[0];
1269       }
1270
1271       if (removeGreaterThanYFlag) {
1272         index[1] = Ap[1];
1273         size[1] = region.GetSize()[1]-(index[1]-region.GetIndex()[1]);
1274       }
1275       else {
1276         index[1] = region.GetIndex()[1];
1277         size[1] = Ap[1] - index[1];
1278       }
1279
1280       // Set region
1281       region.SetSize(size);
1282       region.SetIndex(index);
1283
1284       // Fill region with BG (simple region iterator)
1285       FillRegionWithValue<SliceType>(slices[sliceIndex], BG, region);
1286       // Loop
1287     }
1288
1289     // Merge slices
1290     typename ImageType::Pointer output;
1291     output = clitk::JoinSlices<ImageType>(slices, input, sliceDim);
1292     return output;
1293   }
1294   //--------------------------------------------------------------------
1295
1296
1297   //--------------------------------------------------------------------
1298   template<class ImageType>
1299   void
1300   FillRegionWithValue(ImageType * input, typename ImageType::PixelType value, typename ImageType::RegionType & region)
1301   {
1302     typedef itk::ImageRegionIterator<ImageType> IteratorType;
1303     IteratorType iter(input, region);
1304     iter.GoToBegin();
1305     while (!iter.IsAtEnd()) {
1306       iter.Set(value);
1307       ++iter;
1308     }
1309   }
1310   //--------------------------------------------------------------------
1311
1312
1313   //--------------------------------------------------------------------
1314   template<class ImageType>
1315   void
1316   GetMinMaxBoundary(ImageType * input, typename ImageType::PointType & min,
1317                     typename ImageType::PointType & max)
1318   {
1319     typedef typename ImageType::PointType PointType;
1320     typedef typename ImageType::IndexType IndexType;
1321     IndexType min_i, max_i;
1322     min_i = input->GetLargestPossibleRegion().GetIndex();
1323     for(uint i=0; i<ImageType::ImageDimension; i++)
1324       max_i[i] = input->GetLargestPossibleRegion().GetSize()[i] + min_i[i];
1325     input->TransformIndexToPhysicalPoint(min_i, min);
1326     input->TransformIndexToPhysicalPoint(max_i, max);
1327   }
1328   //--------------------------------------------------------------------
1329
1330 } // end of namespace
1331