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
  36 namespace clitk {
  37
  38   //--------------------------------------------------------------------
  39   template<class ImageType>
  40   void ComputeBBFromImageRegion(const ImageType * image,
  41                                 typename ImageType::RegionType region,
  42                                 typename itk::BoundingBox<unsigned long,
  43                                 ImageType::ImageDimension>::Pointer bb) {
  44     typedef typename ImageType::IndexType IndexType;
  45     IndexType firstIndex;
  46     IndexType lastIndex;
  47     for(unsigned int i=0; i<image->GetImageDimension(); i++) {
  48       firstIndex[i] = region.GetIndex()[i];
  49       lastIndex[i] = firstIndex[i]+region.GetSize()[i];
  50     }
  51
  52     typedef itk::BoundingBox<unsigned long,
  53                              ImageType::ImageDimension> BBType;
  54     typedef typename BBType::PointType PointType;
  55     PointType lastPoint;
  56     PointType firstPoint;
  57     image->TransformIndexToPhysicalPoint(firstIndex, firstPoint);
  58     image->TransformIndexToPhysicalPoint(lastIndex, lastPoint);
  59
  60     bb->SetMaximum(lastPoint);
  61     bb->SetMinimum(firstPoint);
  62   }
  63   //--------------------------------------------------------------------
  64
  65
  66   //--------------------------------------------------------------------
  67   template<int Dimension>
  68   void ComputeBBIntersection(typename itk::BoundingBox<unsigned long, Dimension>::Pointer bbo,
  69                              typename itk::BoundingBox<unsigned long, Dimension>::Pointer bbi1,
  70                              typename itk::BoundingBox<unsigned long, Dimension>::Pointer bbi2) {
  71
  72     typedef itk::BoundingBox<unsigned long, Dimension> BBType;
  73     typedef typename BBType::PointType PointType;
  74     PointType lastPoint;
  75     PointType firstPoint;
  76
  77     for(unsigned int i=0; i<Dimension; i++) {
  78       firstPoint[i] = std::max(bbi1->GetMinimum()[i],
  79                                bbi2->GetMinimum()[i]);
  80       lastPoint[i] = std::min(bbi1->GetMaximum()[i],
  81                               bbi2->GetMaximum()[i]);
  82     }
  83
  84     bbo->SetMaximum(lastPoint);
  85     bbo->SetMinimum(firstPoint);
  86   }
  87   //--------------------------------------------------------------------
  88
  89
  90   //--------------------------------------------------------------------
  91   template<class ImageType>
  92   void ComputeRegionFromBB(const ImageType * image,
  93                            const typename itk::BoundingBox<unsigned long,
  94                                                            ImageType::ImageDimension>::Pointer bb,
  95                            typename ImageType::RegionType & region) {
  96     // Types
  97     typedef typename ImageType::IndexType  IndexType;
  98     typedef typename ImageType::PointType  PointType;
  99     typedef typename ImageType::RegionType RegionType;
 100     typedef typename ImageType::SizeType   SizeType;
 101
 102     // Region starting point
 103     IndexType regionStart;
 104     PointType start = bb->GetMinimum();
 105     image->TransformPhysicalPointToIndex(start, regionStart);
 106
 107     // Region size
 108     SizeType regionSize;
 109     PointType maxs = bb->GetMaximum();
 110     PointType mins = bb->GetMinimum();
 111     for(unsigned int i=0; i<ImageType::ImageDimension; i++) {
 112       regionSize[i] = lrint((maxs[i] - mins[i])/image->GetSpacing()[i]);
 113     }
 114
 115     // Create region
 116     region.SetIndex(regionStart);
 117     region.SetSize(regionSize);
 118   }
 119   //--------------------------------------------------------------------
 120
 121   //--------------------------------------------------------------------
 122   template<class ImageType, class TMaskImageType>
 123   typename ImageType::Pointer
 124   SetBackground(const ImageType * input,
 125                 const TMaskImageType * mask,
 126                 typename TMaskImageType::PixelType maskBG,
 127                 typename ImageType::PixelType outValue,
 128                 bool inPlace) {
 129     typedef SetBackgroundImageFilter<ImageType, TMaskImageType, ImageType>
 130       SetBackgroundImageFilterType;
 131     typename SetBackgroundImageFilterType::Pointer setBackgroundFilter
 132       = SetBackgroundImageFilterType::New();
 133     //  if (inPlace) setBackgroundFilter->ReleaseDataFlagOn(); // No seg fault
 134     setBackgroundFilter->SetInPlace(inPlace); // This is important to keep memory low
 135     setBackgroundFilter->SetInput(input);
 136     setBackgroundFilter->SetInput2(mask);
 137     setBackgroundFilter->SetMaskValue(maskBG);
 138     setBackgroundFilter->SetOutsideValue(outValue);
 139     setBackgroundFilter->Update();
 140     return setBackgroundFilter->GetOutput();
 141   }
 142   //--------------------------------------------------------------------
 143
 144
 145   //--------------------------------------------------------------------
 146   template<class ImageType>
 147   int GetNumberOfConnectedComponentLabels(const ImageType * input,
 148                                           typename ImageType::PixelType BG,
 149                                           bool isFullyConnected) {
 150     // Connected Component label
 151     typedef itk::ConnectedComponentImageFilter<ImageType, ImageType> ConnectFilterType;
 152     typename ConnectFilterType::Pointer connectFilter = ConnectFilterType::New();
 153     connectFilter->SetInput(input);
 154     connectFilter->SetBackgroundValue(BG);
 155     connectFilter->SetFullyConnected(isFullyConnected);
 156     connectFilter->Update();
 157
 158     // Return result
 159     return connectFilter->GetObjectCount();
 160   }
 161   //--------------------------------------------------------------------
 162
 163   //--------------------------------------------------------------------
 164   /*
 165     Warning : in this cas, we consider outputType like inputType, not
 166     InternalImageType. Be sure it fits.
 167   */
 168   template<class ImageType>
 169   typename ImageType::Pointer
 170   Labelize(const ImageType * input,
 171            typename ImageType::PixelType BG,
 172            bool isFullyConnected,
 173            int minimalComponentSize) {
 174     // InternalImageType for storing large number of component
 175     typedef itk::Image<int, ImageType::ImageDimension> InternalImageType;
 176
 177     // Connected Component label
 178     typedef itk::ConnectedComponentImageFilter<ImageType, InternalImageType> ConnectFilterType;
 179     typename ConnectFilterType::Pointer connectFilter = ConnectFilterType::New();
 180     //  connectFilter->ReleaseDataFlagOn();
 181     connectFilter->SetInput(input);
 182     connectFilter->SetBackgroundValue(BG);
 183     connectFilter->SetFullyConnected(isFullyConnected);
 184
 185     // Sort by size and remove too small area.
 186     typedef itk::RelabelComponentImageFilter<InternalImageType, ImageType> RelabelFilterType;
 187     typename RelabelFilterType::Pointer relabelFilter = RelabelFilterType::New();
 188     //  relabelFilter->ReleaseDataFlagOn(); // if yes, fail when ExplosionControlledThresholdConnectedImageFilter ???
 189     relabelFilter->SetInput(connectFilter->GetOutput());
 190     relabelFilter->SetMinimumObjectSize(minimalComponentSize);
 191     relabelFilter->Update();
 192
 193     // Return result
 194     typename ImageType::Pointer output = relabelFilter->GetOutput();
 195     return output;
 196   }
 197   //--------------------------------------------------------------------
 198
 199
 200   //--------------------------------------------------------------------
 201   /*
 202     Warning : in this cas, we consider outputType like inputType, not
 203     InternalImageType. Be sure it fits.
 204   */
 205   template<class ImageType>
 206   typename ImageType::Pointer
 207   LabelizeAndCountNumberOfObjects(const ImageType * input,
 208                                   typename ImageType::PixelType BG,
 209                                   bool isFullyConnected,
 210                                   int minimalComponentSize,
 211                                   int & nb) {
 212     // InternalImageType for storing large number of component
 213     typedef itk::Image<int, ImageType::ImageDimension> InternalImageType;
 214
 215     // Connected Component label
 216     typedef itk::ConnectedComponentImageFilter<ImageType, InternalImageType> ConnectFilterType;
 217     typename ConnectFilterType::Pointer connectFilter = ConnectFilterType::New();
 218     //  connectFilter->ReleaseDataFlagOn();
 219     connectFilter->SetInput(input);
 220     connectFilter->SetBackgroundValue(BG);
 221     connectFilter->SetFullyConnected(isFullyConnected);
 222
 223     // Sort by size and remove too small area.
 224     typedef itk::RelabelComponentImageFilter<InternalImageType, ImageType> RelabelFilterType;
 225     typename RelabelFilterType::Pointer relabelFilter = RelabelFilterType::New();
 226     //  relabelFilter->ReleaseDataFlagOn(); // if yes, fail when ExplosionControlledThresholdConnectedImageFilter ???
 227     relabelFilter->SetInput(connectFilter->GetOutput());
 228     relabelFilter->SetMinimumObjectSize(minimalComponentSize);
 229     relabelFilter->Update();
 230
 231     nb = relabelFilter->GetNumberOfObjects();
 232     // DD(relabelFilter->GetOriginalNumberOfObjects());
 233     // DD(relabelFilter->GetSizeOfObjectsInPhysicalUnits()[0]);
 234
 235     // Return result
 236     typename ImageType::Pointer output = relabelFilter->GetOutput();
 237     return output;
 238   }
 239   //--------------------------------------------------------------------
 240
 241
 242
 243   //--------------------------------------------------------------------
 244   template<class ImageType>
 245   typename ImageType::Pointer
 246   RemoveLabels(const ImageType * input,
 247                typename ImageType::PixelType BG,
 248                std::vector<typename ImageType::PixelType> & labelsToRemove) {
 249     assert(labelsToRemove.size() != 0);
 250     typename ImageType::Pointer working_image;// = input;
 251     for (unsigned int i=0; i <labelsToRemove.size(); i++) {
 252       typedef SetBackgroundImageFilter<ImageType, ImageType> SetBackgroundImageFilterType;
 253       typename SetBackgroundImageFilterType::Pointer setBackgroundFilter = SetBackgroundImageFilterType::New();
 254       setBackgroundFilter->SetInput(input);
 255       setBackgroundFilter->SetInput2(input);
 256       setBackgroundFilter->SetMaskValue(labelsToRemove[i]);
 257       setBackgroundFilter->SetOutsideValue(BG);
 258       setBackgroundFilter->Update();
 259       working_image = setBackgroundFilter->GetOutput();
 260     }
 261     return working_image;
 262   }
 263   //--------------------------------------------------------------------
 264
 265
 266   //--------------------------------------------------------------------
 267   template<class ImageType>
 268   typename ImageType::Pointer
 269   KeepLabels(const ImageType * input,
 270              typename ImageType::PixelType BG,
 271              typename ImageType::PixelType FG,
 272              typename ImageType::PixelType firstKeep,
 273              typename ImageType::PixelType lastKeep,
 274              bool useLastKeep) {
 275     typedef itk::BinaryThresholdImageFilter<ImageType, ImageType> BinarizeFilterType;
 276     typename BinarizeFilterType::Pointer binarizeFilter = BinarizeFilterType::New();
 277     binarizeFilter->SetInput(input);
 278     binarizeFilter->SetLowerThreshold(firstKeep);
 279     if (useLastKeep) binarizeFilter->SetUpperThreshold(lastKeep);
 280     binarizeFilter->SetInsideValue(FG);
 281     binarizeFilter->SetOutsideValue(BG);
 282     binarizeFilter->Update();
 283     return binarizeFilter->GetOutput();
 284   }
 285   //--------------------------------------------------------------------
 286
 287
 288   //--------------------------------------------------------------------
 289   template<class ImageType>
 290   typename ImageType::Pointer
 291   LabelizeAndSelectLabels(const ImageType * input,
 292                           typename ImageType::PixelType BG,
 293                           typename ImageType::PixelType FG,
 294                           bool isFullyConnected,
 295                           int minimalComponentSize,
 296                           LabelizeParameters<typename ImageType::PixelType> * param)
 297   {
 298     typename ImageType::Pointer working_image;
 299     working_image = Labelize<ImageType>(input, BG, isFullyConnected, minimalComponentSize);
 300     if (param->GetLabelsToRemove().size() != 0)
 301       working_image = RemoveLabels<ImageType>(working_image, BG, param->GetLabelsToRemove());
 302     working_image = KeepLabels<ImageType>(working_image,
 303                                           BG, FG,
 304                                           param->GetFirstKeep(),
 305                                           param->GetLastKeep(),
 306                                           param->GetUseLastKeep());
 307     return working_image;
 308   }
 309   //--------------------------------------------------------------------
 310
 311
 312   //--------------------------------------------------------------------
 313   template<class ImageType>
 314   typename ImageType::Pointer
 315   ResizeImageLike(const ImageType * input,
 316                   const itk::ImageBase<ImageType::ImageDimension> * like,
 317                   typename ImageType::PixelType backgroundValue)
 318   {
 319     typedef CropLikeImageFilter<ImageType> CropFilterType;
 320     typename CropFilterType::Pointer cropFilter = CropFilterType::New();
 321     cropFilter->SetInput(input);
 322     cropFilter->SetCropLikeImage(like);
 323     cropFilter->SetBackgroundValue(backgroundValue);
 324     cropFilter->Update();
 325     return cropFilter->GetOutput();
 326   }
 327   //--------------------------------------------------------------------
 328
 329
 330   //--------------------------------------------------------------------
 331   template<class MaskImageType>
 332   typename MaskImageType::Pointer
 333   SliceBySliceRelativePosition(const MaskImageType * input,
 334                                const MaskImageType * object,
 335                                int direction,
 336                                double threshold,
 337                                std::string orientation,
 338                                bool uniqueConnectedComponent,
 339                                double spacing,
 340                                bool autocropFlag,
 341                                bool singleObjectCCL)
 342   {
 343     typedef SliceBySliceRelativePositionFilter<MaskImageType> SliceRelPosFilterType;
 344     typename SliceRelPosFilterType::Pointer sliceRelPosFilter = SliceRelPosFilterType::New();
 345     sliceRelPosFilter->VerboseStepFlagOff();
 346     sliceRelPosFilter->WriteStepFlagOff();
 347     sliceRelPosFilter->SetInput(input);
 348     sliceRelPosFilter->SetInputObject(object);
 349     sliceRelPosFilter->SetDirection(direction);
 350     sliceRelPosFilter->SetFuzzyThreshold(threshold);
 351     sliceRelPosFilter->AddOrientationTypeString(orientation);
 352     sliceRelPosFilter->SetIntermediateSpacingFlag((spacing != -1));
 353     sliceRelPosFilter->SetIntermediateSpacing(spacing);
 354     sliceRelPosFilter->SetUniqueConnectedComponentBySlice(uniqueConnectedComponent);
 355     sliceRelPosFilter->SetUseASingleObjectConnectedComponentBySliceFlag(singleObjectCCL);
 356     //    sliceRelPosFilter->SetInverseOrientationFlag(inverseflag);
 357     sliceRelPosFilter->SetAutoCropFlag(autocropFlag);
 358     sliceRelPosFilter->IgnoreEmptySliceObjectFlagOn();
 359     sliceRelPosFilter->Update();
 360     return sliceRelPosFilter->GetOutput();
 361   }
 362   //--------------------------------------------------------------------
 363
 364   //--------------------------------------------------------------------
 365   template<class ImageType>
 366   bool
 367   FindExtremaPointInAGivenDirection(const ImageType * input,
 368                                     typename ImageType::PixelType bg,
 369                                     int direction, bool opposite,
 370                                     typename ImageType::PointType & point)
 371   {
 372     typename ImageType::PointType dummy;
 373     return FindExtremaPointInAGivenDirection(input, bg, direction,
 374                                              opposite, dummy, 0, point);
 375   }
 376   //--------------------------------------------------------------------
 377
 378   //--------------------------------------------------------------------
 379   template<class ImageType>
 380   bool
 381   FindExtremaPointInAGivenDirection(const ImageType * input,
 382                                     typename ImageType::PixelType bg,
 383                                     int direction, bool opposite,
 384                                     typename ImageType::PointType refpoint,
 385                                     double distanceMax,
 386                                     typename ImageType::PointType & point)
 387   {
 388     /*
 389       loop over input pixels, store the index in the fg that is max
 390       according to the given direction.
 391     */
 392     typedef itk::ImageRegionConstIteratorWithIndex<ImageType> IteratorType;
 393     IteratorType iter(input, input->GetLargestPossibleRegion());
 394     iter.GoToBegin();
 395     typename ImageType::IndexType max = input->GetLargestPossibleRegion().GetIndex();
 396     if (opposite) max = max+input->GetLargestPossibleRegion().GetSize();
 397     bool found=false;
 398     while (!iter.IsAtEnd()) {
 399       if (iter.Get() != bg) {
 400         bool test = iter.GetIndex()[direction] >  max[direction];
 401         if (opposite) test = !test;
 402         if (test) {
 403           typename ImageType::PointType p;
 404           input->TransformIndexToPhysicalPoint(iter.GetIndex(), p);
 405           if ((distanceMax==0) || (p.EuclideanDistanceTo(refpoint) < distanceMax)) {
 406             max = iter.GetIndex();
 407             found = true;
 408           }
 409         }
 410       }
 411       ++iter;
 412     }
 413     if (!found) return false;
 414     input->TransformIndexToPhysicalPoint(max, point);
 415     return true;
 416   }
 417   //--------------------------------------------------------------------
 418
 419
 420   //--------------------------------------------------------------------
 421   template<class ImageType>
 422   typename ImageType::Pointer
 423   CropImageRemoveGreaterThan(const ImageType * image,
 424                  int dim, double min, bool autoCrop,
 425                  typename ImageType::PixelType BG)
 426   {
 427     return CropImageAlongOneAxis<ImageType>(image, dim,
 428                                             image->GetOrigin()[dim],
 429                                             min,
 430                                             autoCrop, BG);
 431   }
 432   //--------------------------------------------------------------------
 433
 434
 435   //--------------------------------------------------------------------
 436   template<class ImageType>
 437   typename ImageType::Pointer
 438   CropImageRemoveLowerThan(const ImageType * image,
 439                  int dim, double max, bool autoCrop,
 440                  typename ImageType::PixelType BG)
 441   {
 442     typename ImageType::PointType p;
 443     image->TransformIndexToPhysicalPoint(image->GetLargestPossibleRegion().GetIndex()+
 444                                          image->GetLargestPossibleRegion().GetSize(), p);
 445     return CropImageAlongOneAxis<ImageType>(image, dim, max, p[dim], autoCrop, BG);
 446   }
 447   //--------------------------------------------------------------------
 448
 449
 450   //--------------------------------------------------------------------
 451   template<class ImageType>
 452   typename ImageType::Pointer
 453   CropImageAlongOneAxis(const ImageType * image,
 454                         int dim, double min, double max,
 455                         bool autoCrop, typename ImageType::PixelType BG)
 456   {
 457     // Compute region size
 458     typename ImageType::RegionType region;
 459     typename ImageType::SizeType size = image->GetLargestPossibleRegion().GetSize();
 460     typename ImageType::PointType p = image->GetOrigin();
 461     p[dim] = min;
 462     typename ImageType::IndexType start;
 463     image->TransformPhysicalPointToIndex(p, start);
 464     p[dim] = max;
 465     typename ImageType::IndexType end;
 466     image->TransformPhysicalPointToIndex(p, end);
 467     size[dim] = fabs(end[dim]-start[dim]);
 468     region.SetIndex(start);
 469     region.SetSize(size);
 470
 471     // Perform Crop
 472     typedef itk::RegionOfInterestImageFilter<ImageType, ImageType> CropFilterType;
 473     typename CropFilterType::Pointer cropFilter = CropFilterType::New();
 474     cropFilter->SetInput(image);
 475     cropFilter->SetRegionOfInterest(region);
 476     cropFilter->Update();
 477     typename ImageType::Pointer result = cropFilter->GetOutput();
 478
 479     // Auto Crop
 480     if (autoCrop) {
 481       result = AutoCrop<ImageType>(result, BG);
 482     }
 483     return result;
 484   }
 485   //--------------------------------------------------------------------
 486
 487
 488   //--------------------------------------------------------------------
 489   template<class ImageType>
 490   void
 491   ComputeCentroids(const ImageType * image,
 492                    typename ImageType::PixelType BG,
 493                    std::vector<typename ImageType::PointType> & centroids)
 494   {
 495     typedef long LabelType;
 496     static const unsigned int Dim = ImageType::ImageDimension;
 497     typedef itk::ShapeLabelObject< LabelType, Dim > LabelObjectType;
 498     typedef itk::LabelMap< LabelObjectType > LabelMapType;
 499     typedef itk::LabelImageToLabelMapFilter<ImageType, LabelMapType> ImageToMapFilterType;
 500     typename ImageToMapFilterType::Pointer imageToLabelFilter = ImageToMapFilterType::New();
 501     typedef itk::ShapeLabelMapFilter<LabelMapType, ImageType> ShapeFilterType;
 502     typename ShapeFilterType::Pointer statFilter = ShapeFilterType::New();
 503     imageToLabelFilter->SetBackgroundValue(BG);
 504     imageToLabelFilter->SetInput(image);
 505     statFilter->SetInput(imageToLabelFilter->GetOutput());
 506     statFilter->Update();
 507     typename LabelMapType::Pointer labelMap = statFilter->GetOutput();
 508
 509     centroids.clear();
 510     typename ImageType::PointType dummy;
 511     centroids.push_back(dummy); // label 0 -> no centroid, use dummy point for BG
 512     //DS FIXME (not useful ! to change ..)
 513     DD(labelMap->GetNumberOfLabelObjects());
 514     for(uint i=0; i<labelMap->GetNumberOfLabelObjects(); i++) {
 515       int label = labelMap->GetLabels()[i];
 516       centroids.push_back(labelMap->GetLabelObject(label)->GetCentroid());
 517     }
 518   }
 519   //--------------------------------------------------------------------
 520
 521
 522   //--------------------------------------------------------------------
 523   template<class ImageType>
 524   void
 525   ComputeCentroids2(const ImageType * image,
 526                    typename ImageType::PixelType BG,
 527                    std::vector<typename ImageType::PointType> & centroids)
 528   {
 529     typedef long LabelType;
 530     static const unsigned int Dim = ImageType::ImageDimension;
 531     typedef itk::ShapeLabelObject< LabelType, Dim > LabelObjectType;
 532     typedef itk::LabelMap< LabelObjectType > LabelMapType;
 533     typedef itk::LabelImageToLabelMapFilter<ImageType, LabelMapType> ImageToMapFilterType;
 534     typename ImageToMapFilterType::Pointer imageToLabelFilter = ImageToMapFilterType::New();
 535     typedef itk::ShapeLabelMapFilter<LabelMapType, ImageType> ShapeFilterType;
 536     typename ShapeFilterType::Pointer statFilter = ShapeFilterType::New();
 537     imageToLabelFilter->SetBackgroundValue(BG);
 538     imageToLabelFilter->SetInput(image);
 539     statFilter->SetInput(imageToLabelFilter->GetOutput());
 540     statFilter->Update();
 541     typename LabelMapType::Pointer labelMap = statFilter->GetOutput();
 542
 543     centroids.clear();
 544     typename ImageType::PointType dummy;
 545     centroids.push_back(dummy); // label 0 -> no centroid, use dummy point
 546     for(uint i=1; i<labelMap->GetNumberOfLabelObjects()+1; i++) {
 547       centroids.push_back(labelMap->GetLabelObject(i)->GetCentroid());
 548     }
 549
 550     for(uint i=1; i<labelMap->GetNumberOfLabelObjects()+1; i++) {
 551       DD(labelMap->GetLabelObject(i)->GetBinaryPrincipalAxes());
 552       DD(labelMap->GetLabelObject(i)->GetBinaryFlatness());
 553       DD(labelMap->GetLabelObject(i)->GetRoundness ());
 554
 555       // search for the point on the boundary alog PA
 556
 557     }
 558
 559   }
 560   //--------------------------------------------------------------------
 561
 562
 563   //--------------------------------------------------------------------
 564   template<class ImageType>
 565   void
 566   ExtractSlices(const ImageType * image, int direction,
 567                 std::vector<typename itk::Image<typename ImageType::PixelType,
 568                                                 ImageType::ImageDimension-1>::Pointer > & slices)
 569   {
 570     typedef ExtractSliceFilter<ImageType> ExtractSliceFilterType;
 571     typedef typename ExtractSliceFilterType::SliceType SliceType;
 572     typename ExtractSliceFilterType::Pointer
 573       extractSliceFilter = ExtractSliceFilterType::New();
 574     extractSliceFilter->SetInput(image);
 575     extractSliceFilter->SetDirection(direction);
 576     extractSliceFilter->Update();
 577     extractSliceFilter->GetOutputSlices(slices);
 578   }
 579   //--------------------------------------------------------------------
 580
 581
 582   //--------------------------------------------------------------------
 583   template<class ImageType>
 584   void
 585   PointsUtils<ImageType>::Convert2DTo3D(const PointType2D & p2D,
 586                                         const ImageType * image,
 587                                         const int slice,
 588                                         PointType3D & p3D)
 589   {
 590     IndexType3D index3D;
 591     index3D[0] = index3D[1] = 0;
 592     index3D[2] = image->GetLargestPossibleRegion().GetIndex()[2]+slice;
 593     image->TransformIndexToPhysicalPoint(index3D, p3D);
 594     p3D[0] = p2D[0];
 595     p3D[1] = p2D[1];
 596     //  p3D[2] = p[2];//(image->GetLargestPossibleRegion().GetIndex()[2]+slice)*image->GetSpacing()[2]
 597     //    + image->GetOrigin()[2];
 598   }
 599   //--------------------------------------------------------------------
 600
 601
 602   //--------------------------------------------------------------------
 603   template<class ImageType>
 604   void
 605   PointsUtils<ImageType>::Convert2DMapTo3DList(const MapPoint2DType & map,
 606                                             const ImageType * image,
 607                                             VectorPoint3DType & list)
 608   {
 609     typename MapPoint2DType::const_iterator iter = map.begin();
 610     while (iter != map.end()) {
 611       PointType3D p;
 612       Convert2DTo3D(iter->second, image, iter->first, p);
 613       list.push_back(p);
 614       ++iter;
 615     }
 616   }
 617   //--------------------------------------------------------------------
 618
 619
 620   //--------------------------------------------------------------------
 621   template<class ImageType>
 622   void
 623   PointsUtils<ImageType>::Convert2DListTo3DList(const VectorPoint2DType & p2D,
 624                                                 int slice,
 625                                                 const ImageType * image,
 626                                                 VectorPoint3DType & list)
 627   {
 628     for(uint i=0; i<p2D.size(); i++) {
 629       PointType3D p;
 630       Convert2DTo3D(p2D[i], image, slice, p);
 631       list.push_back(p);
 632     }
 633   }
 634   //--------------------------------------------------------------------
 635
 636
 637   //--------------------------------------------------------------------
 638   template<class ImageType>
 639   void
 640   WriteListOfLandmarks(std::vector<typename ImageType::PointType> points,
 641                        std::string filename)
 642   {
 643     std::ofstream os;
 644     openFileForWriting(os, filename);
 645     os << "LANDMARKS1" << std::endl;
 646     for(uint i=0; i<points.size(); i++) {
 647       const typename ImageType::PointType & p = points[i];
 648       // Write it in the file
 649       os << i << " " << p[0] << " " << p[1] << " " << p[2] << " 0 0 " << std::endl;
 650     }
 651     os.close();
 652   }
 653   //--------------------------------------------------------------------
 654
 655
 656   //--------------------------------------------------------------------
 657   template<class ImageType>
 658   typename ImageType::Pointer
 659   Dilate(const ImageType * image, double radiusInMM,
 660          typename ImageType::PixelType BG,
 661          typename ImageType::PixelType FG,
 662          bool extendSupport)
 663   {
 664     typename ImageType::SizeType r;
 665     for(uint i=0; i<ImageType::ImageDimension; i++)
 666       r[i] = (uint)lrint(radiusInMM/image->GetSpacing()[i]);
 667     return Dilate<ImageType>(image, r, BG, FG, extendSupport);
 668   }
 669   //--------------------------------------------------------------------
 670
 671
 672   //--------------------------------------------------------------------
 673   template<class ImageType>
 674   typename ImageType::Pointer
 675   Dilate(const ImageType * image, typename ImageType::PointType radiusInMM,
 676          typename ImageType::PixelType BG,
 677          typename ImageType::PixelType FG,
 678          bool extendSupport)
 679   {
 680     typename ImageType::SizeType r;
 681     for(uint i=0; i<ImageType::ImageDimension; i++)
 682       r[i] = (uint)lrint(radiusInMM[i]/image->GetSpacing()[i]);
 683     return Dilate<ImageType>(image, r, BG, FG, extendSupport);
 684   }
 685   //--------------------------------------------------------------------
 686
 687
 688   //--------------------------------------------------------------------
 689   template<class ImageType>
 690   typename ImageType::Pointer
 691   Dilate(const ImageType * image, typename ImageType::SizeType radius,
 692          typename ImageType::PixelType BG,
 693          typename ImageType::PixelType FG,
 694          bool extendSupport)
 695   {
 696     // Create kernel for dilatation
 697     typedef itk::BinaryBallStructuringElement<typename ImageType::PixelType,
 698                                               ImageType::ImageDimension> KernelType;
 699     KernelType structuringElement;
 700     structuringElement.SetRadius(radius);
 701     structuringElement.CreateStructuringElement();
 702
 703     typename ImageType::Pointer output;
 704     if (extendSupport) {
 705       typedef itk::ConstantPadImageFilter<ImageType, ImageType> PadFilterType;
 706       typename PadFilterType::Pointer padFilter = PadFilterType::New();
 707       padFilter->SetInput(image);
 708       typename ImageType::SizeType lower;
 709       typename ImageType::SizeType upper;
 710       for(uint i=0; i<3; i++) {
 711         lower[i] = upper[i] = 2*(radius[i]+1);
 712       }
 713       padFilter->SetPadLowerBound(lower);
 714       padFilter->SetPadUpperBound(upper);
 715       padFilter->Update();
 716       output = padFilter->GetOutput();
 717     }
 718
 719     // Dilate  filter
 720     typedef itk::BinaryDilateImageFilter<ImageType, ImageType , KernelType> DilateFilterType;
 721     typename DilateFilterType::Pointer dilateFilter = DilateFilterType::New();
 722     dilateFilter->SetBackgroundValue(BG);
 723     dilateFilter->SetForegroundValue(FG);
 724     dilateFilter->SetBoundaryToForeground(false);
 725     dilateFilter->SetKernel(structuringElement);
 726     dilateFilter->SetInput(output);
 727     dilateFilter->Update();
 728     return dilateFilter->GetOutput();
 729   }
 730   //--------------------------------------------------------------------
 731
 732
 733   //--------------------------------------------------------------------
 734   template<class ValueType, class VectorType>
 735   void ConvertOption(std::string optionName, uint given,
 736                      ValueType * values, VectorType & p,
 737                      uint dim, bool required)
 738   {
 739     if (required && (given == 0)) {
 740       clitkExceptionMacro("The option --" << optionName << " must be set and have 1 or "
 741                           << dim << " values.");
 742     }
 743     if (given == 1) {
 744       for(uint i=0; i<dim; i++) p[i] = values[0];
 745       return;
 746     }
 747     if (given == dim) {
 748       for(uint i=0; i<dim; i++) p[i] = values[i];
 749       return;
 750     }
 751     if (given == 0) return;
 752     clitkExceptionMacro("The option --" << optionName << " must have 1 or "
 753                         << dim << " values.");
 754   }
 755   //--------------------------------------------------------------------
 756
 757
 758   //--------------------------------------------------------------------
 759   /*
 760     http://www.gamedev.net/community/forums/topic.asp?topic_id=542870
 761     Assuming the points are (Ax,Ay) (Bx,By) and (Cx,Cy), you need to compute:
 762     (Bx - Ax) * (Cy - Ay) - (By - Ay) * (Cx - Ax)
 763     This will equal zero if the point C is on the line formed by
 764     points A and B, and will have a different sign depending on the
 765     side. Which side this is depends on the orientation of your (x,y)
 766     coordinates, but you can plug test values for A,B and C into this
 767     formula to determine whether negative values are to the left or to
 768     the right.
 769     => to accelerate, start with formula, when change sign -> stop and fill
 770
 771     offsetToKeep = is used to determine which side of the line we
 772     keep. The point along the mainDirection but 'offsetToKeep' mm away
 773     is kept.
 774
 775   */
 776   template<class ImageType>
 777   void
 778   SliceBySliceSetBackgroundFromLineSeparation(ImageType * input,
 779                                               std::vector<typename ImageType::PointType> & lA,
 780                                               std::vector<typename ImageType::PointType> & lB,
 781                                               typename ImageType::PixelType BG,
 782                                               int mainDirection,
 783                                               double offsetToKeep)
 784   {
 785     typedef itk::ImageSliceIteratorWithIndex<ImageType> SliceIteratorType;
 786     SliceIteratorType siter = SliceIteratorType(input,
 787                                                 input->GetLargestPossibleRegion());
 788     siter.SetFirstDirection(0);
 789     siter.SetSecondDirection(1);
 790     siter.GoToBegin();
 791     uint i=0;
 792     typename ImageType::PointType A;
 793     typename ImageType::PointType B;
 794     typename ImageType::PointType C;
 795     assert(lA.size() == lB.size());
 796     while ((i<lA.size()) && (!siter.IsAtEnd())) {
 797       // Check that the current slice correspond to the current point
 798       input->TransformIndexToPhysicalPoint(siter.GetIndex(), C);
 799       // DD(C);
 800       // DD(i);
 801       // DD(lA[i]);
 802       if ((fabs(C[2] - lA[i][2]))>0.01) { // is !equal with a tolerance of 0.01 mm
 803       }
 804       else {
 805         // Define A,B,C points
 806         A = lA[i];
 807         B = lB[i];
 808         C = A;
 809         // DD(A);
 810         // DD(B);
 811         // DD(C);
 812
 813         // Check that the line is not a point (A=B)
 814         bool p = (A[0] == B[0]) && (A[1] == B[1]);
 815
 816         if (!p) {
 817           C[mainDirection] += offsetToKeep; // I know I must keep this point
 818           double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
 819           bool isPositive = s<0;
 820           while (!siter.IsAtEndOfSlice()) {
 821             while (!siter.IsAtEndOfLine()) {
 822               // Very slow, I know ... but image should be very small
 823               input->TransformIndexToPhysicalPoint(siter.GetIndex(), C);
 824               double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
 825               if (s == 0) siter.Set(BG); // on the line, we decide to remove
 826               if (isPositive) {
 827                 if (s > 0) siter.Set(BG);
 828               }
 829               else {
 830                 if (s < 0) siter.Set(BG);
 831               }
 832               ++siter;
 833             }
 834             siter.NextLine();
 835           } // end loop slice
 836         }
 837
 838         ++i;
 839       } // End of current slice
 840       siter.NextSlice();
 841     }
 842   }
 843   //--------------------------------------------------------------------
 844
 845   //--------------------------------------------------------------------
 846   template<class ImageType>
 847   void
 848   AndNot(ImageType * input,
 849          const ImageType * object,
 850          typename ImageType::PixelType BG)
 851   {
 852     typename ImageType::Pointer o;
 853     bool resized=false;
 854     if (!clitk::HaveSameSizeAndSpacing<ImageType, ImageType>(input, object)) {
 855       o = clitk::ResizeImageLike<ImageType>(object, input, BG);
 856       resized = true;
 857     }
 858
 859     typedef clitk::BooleanOperatorLabelImageFilter<ImageType> BoolFilterType;
 860     typename BoolFilterType::Pointer boolFilter = BoolFilterType::New();
 861     boolFilter->InPlaceOn();
 862     boolFilter->SetInput1(input);
 863     if (resized) boolFilter->SetInput2(o);
 864     else boolFilter->SetInput2(object);
 865     boolFilter->SetBackgroundValue1(BG);
 866     boolFilter->SetBackgroundValue2(BG);
 867     boolFilter->SetOperationType(BoolFilterType::AndNot);
 868     boolFilter->Update();
 869   }
 870   //--------------------------------------------------------------------
 871
 872
 873   //--------------------------------------------------------------------
 874   template<class ImageType>
 875   typename ImageType::Pointer
 876   Binarize(const ImageType * input,
 877            typename ImageType::PixelType lower,
 878            typename ImageType::PixelType upper,
 879            typename ImageType::PixelType BG,
 880            typename ImageType::PixelType FG)
 881   {
 882     typedef itk::BinaryThresholdImageFilter<ImageType, ImageType> BinaryThresholdFilterType;
 883     typename BinaryThresholdFilterType::Pointer binarizeFilter = BinaryThresholdFilterType::New();
 884     binarizeFilter->SetInput(input);
 885     binarizeFilter->InPlaceOff();
 886     binarizeFilter->SetLowerThreshold(lower);
 887     binarizeFilter->SetUpperThreshold(upper);
 888     binarizeFilter->SetInsideValue(FG);
 889     binarizeFilter->SetOutsideValue(BG);
 890     binarizeFilter->Update();
 891     return binarizeFilter->GetOutput();
 892   }
 893   //--------------------------------------------------------------------
 894
 895
 896   //--------------------------------------------------------------------
 897   template<class ImageType>
 898   void
 899   GetMinMaxPointPosition(const ImageType * input,
 900                          typename ImageType::PointType & min,
 901                          typename ImageType::PointType & max)
 902   {
 903     typename ImageType::IndexType index = input->GetLargestPossibleRegion().GetIndex();
 904     input->TransformIndexToPhysicalPoint(index, min);
 905     index = index+input->GetLargestPossibleRegion().GetSize();
 906     input->TransformIndexToPhysicalPoint(index, max);
 907   }
 908   //--------------------------------------------------------------------
 909
 910
 911   //--------------------------------------------------------------------
 912   template<class ImageType>
 913   typename ImageType::PointType
 914   FindExtremaPointInAGivenLine(const ImageType * input,
 915                                int dimension,
 916                                bool inverse,
 917                                typename ImageType::PointType p,
 918                                typename ImageType::PixelType BG,
 919                                double distanceMax)
 920   {
 921     // Which direction ?  Increasing or decreasing.
 922     int d=1;
 923     if (inverse) d=-1;
 924
 925     // Transform to pixel index
 926     typename ImageType::IndexType index;
 927     input->TransformPhysicalPointToIndex(p, index);
 928
 929     // Loop while inside the mask;
 930     while (input->GetPixel(index) != BG) {
 931       index[dimension] += d;
 932     }
 933
 934     // Transform back to Physical Units
 935     typename ImageType::PointType result;
 936     input->TransformIndexToPhysicalPoint(index, result);
 937
 938     // Check that is is not too far away
 939     double distance = p.EuclideanDistanceTo(result);
 940     if (distance > distanceMax) {
 941       result = p; // Get back to initial value
 942     }
 943
 944     return result;
 945   }
 946   //--------------------------------------------------------------------
 947
 948
 949   //--------------------------------------------------------------------
 950   template<class PointType>
 951   bool
 952   IsOnTheSameLineSide(PointType C, PointType A, PointType B, PointType like)
 953   {
 954     // Look at the position of point 'like' according to the AB line
 955     double s = (B[0] - A[0]) * (like[1] - A[1]) - (B[1] - A[1]) * (like[0] - A[0]);
 956     bool negative = s<0;
 957
 958     // Look the C position
 959     s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
 960
 961     if (negative && (s<=0)) return true;
 962     if (!negative && (s>=0)) return true;
 963     return false;
 964   }
 965   //--------------------------------------------------------------------
 966
 967
 968   //--------------------------------------------------------------------
 969   /* Consider an input object, for each slice, find the extrema
 970      position according to a given direction and build a line segment
 971      passing throught this point in a given direction.  Output is a
 972      vector of line (from point A to B), for each slice;
 973    */
 974   template<class ImageType>
 975   void
 976   SliceBySliceBuildLineSegmentAccordingToExtremaPosition(const ImageType * input,
 977                                                          typename ImageType::PixelType BG,
 978                                                          int sliceDimension,
 979                                                          int extremaDirection,
 980                                                          bool extremaOppositeFlag,
 981                                                          int lineDirection,
 982                                                          double margin,
 983                                                          std::vector<typename ImageType::PointType> & A,
 984                                                          std::vector<typename ImageType::PointType> & B)
 985   {
 986     // Type of a slice
 987     typedef typename itk::Image<typename ImageType::PixelType, ImageType::ImageDimension-1> SliceType;
 988
 989     // Build the list of slices
 990     std::vector<typename SliceType::Pointer> slices;
 991     clitk::ExtractSlices<ImageType>(input, sliceDimension, slices);
 992
 993     // Build the list of 2D points
 994     std::map<int, typename SliceType::PointType> position2D;
 995     for(uint i=0; i<slices.size(); i++) {
 996       typename SliceType::PointType p;
 997       bool found =
 998         clitk::FindExtremaPointInAGivenDirection<SliceType>(slices[i], BG,
 999                                                             extremaDirection, extremaOppositeFlag, p);
1000       if (found) {
1001         position2D[i] = p;
1002       }
1003     }
1004
1005     // Convert 2D points in slice into 3D points
1006     clitk::PointsUtils<ImageType>::Convert2DMapTo3DList(position2D, input, A);
1007
1008     // Create additional point just right to the previous ones, on the
1009     // given lineDirection, in order to create a horizontal/vertical line.
1010     for(uint i=0; i<A.size(); i++) {
1011       typename ImageType::PointType p = A[i];
1012       p[lineDirection] += 10;
1013       B.push_back(p);
1014       // Margins ?
1015       A[i][1] += margin;
1016       B[i][1] += margin;
1017     }
1018
1019   }
1020   //--------------------------------------------------------------------
1021
1022
1023 } // end of namespace
1024