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   CropImageAbove(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   CropImageBelow(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
 512     for(uint i=1; i<labelMap->GetNumberOfLabelObjects()+1; i++) {
 513       centroids.push_back(labelMap->GetLabelObject(i)->GetCentroid());
 514     }
 515   }
 516   //--------------------------------------------------------------------
 517
 518
 519   //--------------------------------------------------------------------
 520   template<class ImageType>
 521   void
 522   ComputeCentroids2(const ImageType * image,
 523                    typename ImageType::PixelType BG,
 524                    std::vector<typename ImageType::PointType> & centroids)
 525   {
 526     typedef long LabelType;
 527     static const unsigned int Dim = ImageType::ImageDimension;
 528     typedef itk::ShapeLabelObject< LabelType, Dim > LabelObjectType;
 529     typedef itk::LabelMap< LabelObjectType > LabelMapType;
 530     typedef itk::LabelImageToLabelMapFilter<ImageType, LabelMapType> ImageToMapFilterType;
 531     typename ImageToMapFilterType::Pointer imageToLabelFilter = ImageToMapFilterType::New();
 532     typedef itk::ShapeLabelMapFilter<LabelMapType, ImageType> ShapeFilterType;
 533     typename ShapeFilterType::Pointer statFilter = ShapeFilterType::New();
 534     imageToLabelFilter->SetBackgroundValue(BG);
 535     imageToLabelFilter->SetInput(image);
 536     statFilter->SetInput(imageToLabelFilter->GetOutput());
 537     statFilter->Update();
 538     typename LabelMapType::Pointer labelMap = statFilter->GetOutput();
 539
 540     centroids.clear();
 541     typename ImageType::PointType dummy;
 542     centroids.push_back(dummy); // label 0 -> no centroid, use dummy point
 543     for(uint i=1; i<labelMap->GetNumberOfLabelObjects()+1; i++) {
 544       centroids.push_back(labelMap->GetLabelObject(i)->GetCentroid());
 545     }
 546
 547     for(uint i=1; i<labelMap->GetNumberOfLabelObjects()+1; i++) {
 548       DD(labelMap->GetLabelObject(i)->GetBinaryPrincipalAxes());
 549       DD(labelMap->GetLabelObject(i)->GetBinaryFlatness());
 550       DD(labelMap->GetLabelObject(i)->GetRoundness ());
 551
 552       // search for the point on the boundary alog PA
 553
 554     }
 555
 556   }
 557   //--------------------------------------------------------------------
 558
 559
 560   //--------------------------------------------------------------------
 561   template<class ImageType>
 562   void
 563   ExtractSlices(const ImageType * image, int direction,
 564                 std::vector<typename itk::Image<typename ImageType::PixelType,
 565                                                 ImageType::ImageDimension-1>::Pointer > & slices)
 566   {
 567     typedef ExtractSliceFilter<ImageType> ExtractSliceFilterType;
 568     typedef typename ExtractSliceFilterType::SliceType SliceType;
 569     typename ExtractSliceFilterType::Pointer
 570       extractSliceFilter = ExtractSliceFilterType::New();
 571     extractSliceFilter->SetInput(image);
 572     extractSliceFilter->SetDirection(direction);
 573     extractSliceFilter->Update();
 574     extractSliceFilter->GetOutputSlices(slices);
 575   }
 576   //--------------------------------------------------------------------
 577
 578
 579   //--------------------------------------------------------------------
 580   template<class ImageType>
 581   void
 582   PointsUtils<ImageType>::Convert2DTo3D(const PointType2D & p2D,
 583                                         const ImageType * image,
 584                                         const int slice,
 585                                         PointType3D & p3D)
 586   {
 587     IndexType3D index3D;
 588     index3D[0] = index3D[1] = 0;
 589     index3D[2] = image->GetLargestPossibleRegion().GetIndex()[2]+slice;
 590     image->TransformIndexToPhysicalPoint(index3D, p3D);
 591     p3D[0] = p2D[0];
 592     p3D[1] = p2D[1];
 593     //  p3D[2] = p[2];//(image->GetLargestPossibleRegion().GetIndex()[2]+slice)*image->GetSpacing()[2]
 594     //    + image->GetOrigin()[2];
 595   }
 596   //--------------------------------------------------------------------
 597
 598
 599   //--------------------------------------------------------------------
 600   template<class ImageType>
 601   void
 602   PointsUtils<ImageType>::Convert2DMapTo3DList(const MapPoint2DType & map,
 603                                             const ImageType * image,
 604                                             VectorPoint3DType & list)
 605   {
 606     typename MapPoint2DType::const_iterator iter = map.begin();
 607     while (iter != map.end()) {
 608       PointType3D p;
 609       Convert2DTo3D(iter->second, image, iter->first, p);
 610       list.push_back(p);
 611       ++iter;
 612     }
 613   }
 614   //--------------------------------------------------------------------
 615
 616
 617   //--------------------------------------------------------------------
 618   template<class ImageType>
 619   void
 620   PointsUtils<ImageType>::Convert2DListTo3DList(const VectorPoint2DType & p2D,
 621                                                 int slice,
 622                                                 const ImageType * image,
 623                                                 VectorPoint3DType & list)
 624   {
 625     for(uint i=0; i<p2D.size(); i++) {
 626       PointType3D p;
 627       Convert2DTo3D(p2D[i], image, slice, p);
 628       list.push_back(p);
 629     }
 630   }
 631   //--------------------------------------------------------------------
 632
 633
 634   //--------------------------------------------------------------------
 635   template<class ImageType>
 636   void
 637   WriteListOfLandmarks(std::vector<typename ImageType::PointType> points,
 638                        std::string filename)
 639   {
 640     std::ofstream os;
 641     openFileForWriting(os, filename);
 642     os << "LANDMARKS1" << std::endl;
 643     for(uint i=0; i<points.size(); i++) {
 644       const typename ImageType::PointType & p = points[i];
 645       // Write it in the file
 646       os << i << " " << p[0] << " " << p[1] << " " << p[2] << " 0 0 " << std::endl;
 647     }
 648     os.close();
 649   }
 650   //--------------------------------------------------------------------
 651
 652
 653   //--------------------------------------------------------------------
 654   template<class ImageType>
 655   typename ImageType::Pointer
 656   Dilate(const ImageType * image, double radiusInMM,
 657          typename ImageType::PixelType BG,
 658          typename ImageType::PixelType FG,
 659          bool extendSupport)
 660   {
 661     typename ImageType::SizeType r;
 662     for(uint i=0; i<ImageType::ImageDimension; i++)
 663       r[i] = (uint)lrint(radiusInMM/image->GetSpacing()[i]);
 664     return Dilate<ImageType>(image, r, BG, FG, extendSupport);
 665   }
 666   //--------------------------------------------------------------------
 667
 668
 669   //--------------------------------------------------------------------
 670   template<class ImageType>
 671   typename ImageType::Pointer
 672   Dilate(const ImageType * image, typename ImageType::PointType radiusInMM,
 673          typename ImageType::PixelType BG,
 674          typename ImageType::PixelType FG,
 675          bool extendSupport)
 676   {
 677     typename ImageType::SizeType r;
 678     for(uint i=0; i<ImageType::ImageDimension; i++)
 679       r[i] = (uint)lrint(radiusInMM[i]/image->GetSpacing()[i]);
 680     return Dilate<ImageType>(image, r, BG, FG, extendSupport);
 681   }
 682   //--------------------------------------------------------------------
 683
 684
 685   //--------------------------------------------------------------------
 686   template<class ImageType>
 687   typename ImageType::Pointer
 688   Dilate(const ImageType * image, typename ImageType::SizeType radius,
 689          typename ImageType::PixelType BG,
 690          typename ImageType::PixelType FG,
 691          bool extendSupport)
 692   {
 693     // Create kernel for dilatation
 694     typedef itk::BinaryBallStructuringElement<typename ImageType::PixelType,
 695                                               ImageType::ImageDimension> KernelType;
 696     KernelType structuringElement;
 697     structuringElement.SetRadius(radius);
 698     structuringElement.CreateStructuringElement();
 699
 700     typename ImageType::Pointer output;
 701     if (extendSupport) {
 702       typedef itk::ConstantPadImageFilter<ImageType, ImageType> PadFilterType;
 703       typename PadFilterType::Pointer padFilter = PadFilterType::New();
 704       padFilter->SetInput(image);
 705       typename ImageType::SizeType lower;
 706       typename ImageType::SizeType upper;
 707       for(uint i=0; i<3; i++) {
 708         lower[i] = upper[i] = 2*(radius[i]+1);
 709       }
 710       padFilter->SetPadLowerBound(lower);
 711       padFilter->SetPadUpperBound(upper);
 712       padFilter->Update();
 713       output = padFilter->GetOutput();
 714     }
 715
 716     // Dilate  filter
 717     typedef itk::BinaryDilateImageFilter<ImageType, ImageType , KernelType> DilateFilterType;
 718     typename DilateFilterType::Pointer dilateFilter = DilateFilterType::New();
 719     dilateFilter->SetBackgroundValue(BG);
 720     dilateFilter->SetForegroundValue(FG);
 721     dilateFilter->SetBoundaryToForeground(false);
 722     dilateFilter->SetKernel(structuringElement);
 723     dilateFilter->SetInput(output);
 724     dilateFilter->Update();
 725     return dilateFilter->GetOutput();
 726   }
 727   //--------------------------------------------------------------------
 728
 729
 730   //--------------------------------------------------------------------
 731   template<class ValueType, class VectorType>
 732   void ConvertOption(std::string optionName, uint given,
 733                      ValueType * values, VectorType & p,
 734                      uint dim, bool required)
 735   {
 736     if (required && (given == 0)) {
 737       clitkExceptionMacro("The option --" << optionName << " must be set and have 1 or "
 738                           << dim << " values.");
 739     }
 740     if (given == 1) {
 741       for(uint i=0; i<dim; i++) p[i] = values[0];
 742       return;
 743     }
 744     if (given == dim) {
 745       for(uint i=0; i<dim; i++) p[i] = values[i];
 746       return;
 747     }
 748     if (given == 0) return;
 749     clitkExceptionMacro("The option --" << optionName << " must have 1 or "
 750                         << dim << " values.");
 751   }
 752   //--------------------------------------------------------------------
 753
 754
 755   //--------------------------------------------------------------------
 756   /*
 757     http://www.gamedev.net/community/forums/topic.asp?topic_id=542870
 758     Assuming the points are (Ax,Ay) (Bx,By) and (Cx,Cy), you need to compute:
 759     (Bx - Ax) * (Cy - Ay) - (By - Ay) * (Cx - Ax)
 760     This will equal zero if the point C is on the line formed by
 761     points A and B, and will have a different sign depending on the
 762     side. Which side this is depends on the orientation of your (x,y)
 763     coordinates, but you can plug test values for A,B and C into this
 764     formula to determine whether negative values are to the left or to
 765     the right.
 766     => to accelerate, start with formula, when change sign -> stop and fill
 767
 768     offsetToKeep = is used to determine which side of the line we
 769     keep. The point along the mainDirection but 'offsetToKeep' mm away
 770     is kept.
 771
 772   */
 773   template<class ImageType>
 774   void
 775   SliceBySliceSetBackgroundFromLineSeparation(ImageType * input,
 776                                               std::vector<typename ImageType::PointType> & lA,
 777                                               std::vector<typename ImageType::PointType> & lB,
 778                                               typename ImageType::PixelType BG,
 779                                               int mainDirection,
 780                                               double offsetToKeep)
 781   {
 782     typedef itk::ImageSliceIteratorWithIndex<ImageType> SliceIteratorType;
 783     SliceIteratorType siter = SliceIteratorType(input,
 784                                                 input->GetLargestPossibleRegion());
 785     siter.SetFirstDirection(0);
 786     siter.SetSecondDirection(1);
 787     siter.GoToBegin();
 788     uint i=0;
 789     typename ImageType::PointType A;
 790     typename ImageType::PointType B;
 791     typename ImageType::PointType C;
 792     assert(lA.size() == lB.size());
 793     while ((i<lA.size()) && (!siter.IsAtEnd())) {
 794       // Check that the current slice correspond to the current point
 795       input->TransformIndexToPhysicalPoint(siter.GetIndex(), C);
 796       // DD(C);
 797       // DD(i);
 798       // DD(lA[i]);
 799       if ((fabs(C[2] - lA[i][2]))>0.01) { // is !equal with a tolerance of 0.01 mm
 800       }
 801       else {
 802         // Define A,B,C points
 803         A = lA[i];
 804         B = lB[i];
 805         C = A;
 806         // DD(A);
 807         // DD(B);
 808         // DD(C);
 809
 810         // Check that the line is not a point (A=B)
 811         bool p = (A[0] == B[0]) && (A[1] == B[1]);
 812
 813         if (!p) {
 814           C[mainDirection] += offsetToKeep; // I know I must keep this point
 815           double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
 816           bool isPositive = s<0;
 817           while (!siter.IsAtEndOfSlice()) {
 818             while (!siter.IsAtEndOfLine()) {
 819               // Very slow, I know ... but image should be very small
 820               input->TransformIndexToPhysicalPoint(siter.GetIndex(), C);
 821               double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
 822               if (s == 0) siter.Set(BG); // on the line, we decide to remove
 823               if (isPositive) {
 824                 if (s > 0) siter.Set(BG);
 825               }
 826               else {
 827                 if (s < 0) siter.Set(BG);
 828               }
 829               ++siter;
 830             }
 831             siter.NextLine();
 832           } // end loop slice
 833         }
 834
 835         ++i;
 836       } // End of current slice
 837       siter.NextSlice();
 838     }
 839   }
 840   //--------------------------------------------------------------------
 841
 842   //--------------------------------------------------------------------
 843   template<class ImageType>
 844   void
 845   AndNot(ImageType * input,
 846          const ImageType * object,
 847          typename ImageType::PixelType BG)
 848   {
 849     typename ImageType::Pointer o;
 850     bool resized=false;
 851     if (!clitk::HaveSameSizeAndSpacing<ImageType, ImageType>(input, object)) {
 852       o = clitk::ResizeImageLike<ImageType>(object, input, BG);
 853       resized = true;
 854     }
 855
 856     typedef clitk::BooleanOperatorLabelImageFilter<ImageType> BoolFilterType;
 857     typename BoolFilterType::Pointer boolFilter = BoolFilterType::New();
 858     boolFilter->InPlaceOn();
 859     boolFilter->SetInput1(input);
 860     if (resized) boolFilter->SetInput2(o);
 861     else boolFilter->SetInput2(object);
 862     boolFilter->SetBackgroundValue1(BG);
 863     boolFilter->SetBackgroundValue2(BG);
 864     boolFilter->SetOperationType(BoolFilterType::AndNot);
 865     boolFilter->Update();
 866   }
 867   //--------------------------------------------------------------------
 868
 869
 870   //--------------------------------------------------------------------
 871   template<class ImageType>
 872   typename ImageType::Pointer
 873   Binarize(const ImageType * input,
 874            typename ImageType::PixelType lower,
 875            typename ImageType::PixelType upper,
 876            typename ImageType::PixelType BG,
 877            typename ImageType::PixelType FG)
 878   {
 879     typedef itk::BinaryThresholdImageFilter<ImageType, ImageType> BinaryThresholdFilterType;
 880     typename BinaryThresholdFilterType::Pointer binarizeFilter = BinaryThresholdFilterType::New();
 881     binarizeFilter->SetInput(input);
 882     binarizeFilter->SetLowerThreshold(lower);
 883     binarizeFilter->SetUpperThreshold(upper);
 884     binarizeFilter->SetInsideValue(FG);
 885     binarizeFilter->SetOutsideValue(BG);
 886     binarizeFilter->Update();
 887     return binarizeFilter->GetOutput();
 888   }
 889   //--------------------------------------------------------------------
 890
 891
 892   //--------------------------------------------------------------------
 893   template<class ImageType>
 894   void
 895   GetMinMaxPointPosition(const ImageType * input,
 896                          typename ImageType::PointType & min,
 897                          typename ImageType::PointType & max)
 898   {
 899     typename ImageType::IndexType index = input->GetLargestPossibleRegion().GetIndex();
 900     input->TransformIndexToPhysicalPoint(index, min);
 901     index = index+input->GetLargestPossibleRegion().GetSize();
 902     input->TransformIndexToPhysicalPoint(index, max);
 903   }
 904   //--------------------------------------------------------------------
 905
 906
 907   //--------------------------------------------------------------------
 908   template<class ImageType>
 909   typename ImageType::PointType
 910   FindExtremaPointInAGivenLine(const ImageType * input,
 911                                int dimension,
 912                                bool inverse,
 913                                typename ImageType::PointType p,
 914                                typename ImageType::PixelType BG,
 915                                double distanceMax)
 916   {
 917     // Which direction ?  Increasing or decreasing.
 918     int d=1;
 919     if (inverse) d=-1;
 920
 921     // Transform to pixel index
 922     typename ImageType::IndexType index;
 923     input->TransformPhysicalPointToIndex(p, index);
 924
 925     // Loop while inside the mask;
 926     while (input->GetPixel(index) != BG) {
 927       index[dimension] += d;
 928     }
 929
 930     // Transform back to Physical Units
 931     typename ImageType::PointType result;
 932     input->TransformIndexToPhysicalPoint(index, result);
 933
 934     // Check that is is not too far away
 935     double distance = p.EuclideanDistanceTo(result);
 936     if (distance > distanceMax) {
 937       result = p; // Get back to initial value
 938     }
 939
 940     return result;
 941   }
 942   //--------------------------------------------------------------------
 943
 944
 945   //--------------------------------------------------------------------
 946   template<class PointType>
 947   bool
 948   IsOnTheSameLineSide(PointType C, PointType A, PointType B, PointType like)
 949   {
 950     // Look at the position of point 'like' according to the AB line
 951     double s = (B[0] - A[0]) * (like[1] - A[1]) - (B[1] - A[1]) * (like[0] - A[0]);
 952     bool negative = s<0;
 953
 954     // Look the C position
 955     s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
 956
 957     if (negative && (s<=0)) return true;
 958     if (!negative && (s>=0)) return true;
 959     return false;
 960   }
 961   //--------------------------------------------------------------------
 962
 963
 964   //--------------------------------------------------------------------
 965   /* Consider an input object, for each slice, find the extrema
 966      position according to a given direction and build a line segment
 967      passing throught this point in a given direction.  Output is a
 968      vector of line (from point A to B), for each slice;
 969    */
 970   template<class ImageType>
 971   void
 972   SliceBySliceBuildLineSegmentAccordingToExtremaPosition(const ImageType * input,
 973                                                          typename ImageType::PixelType BG,
 974                                                          int sliceDimension,
 975                                                          int extremaDirection,
 976                                                          bool extremaOppositeFlag,
 977                                                          int lineDirection,
 978                                                          double margin,
 979                                                          std::vector<typename ImageType::PointType> & A,
 980                                                          std::vector<typename ImageType::PointType> & B)
 981   {
 982     // Type of a slice
 983     typedef typename itk::Image<typename ImageType::PixelType, ImageType::ImageDimension-1> SliceType;
 984
 985     // Build the list of slices
 986     std::vector<typename SliceType::Pointer> slices;
 987     clitk::ExtractSlices<ImageType>(input, sliceDimension, slices);
 988
 989     // Build the list of 2D points
 990     std::map<int, typename SliceType::PointType> position2D;
 991     for(uint i=0; i<slices.size(); i++) {
 992       typename SliceType::PointType p;
 993       bool found =
 994         clitk::FindExtremaPointInAGivenDirection<SliceType>(slices[i], BG,
 995                                                             extremaDirection, extremaOppositeFlag, p);
 996       if (found) {
 997         position2D[i] = p;
 998       }
 999     }
1000
1001     // Convert 2D points in slice into 3D points
1002     clitk::PointsUtils<ImageType>::Convert2DMapTo3DList(position2D, input, A);
1003
1004     // Create additional point just right to the previous ones, on the
1005     // given lineDirection, in order to create a horizontal/vertical line.
1006     for(uint i=0; i<A.size(); i++) {
1007       typename ImageType::PointType p = A[i];
1008       p[lineDirection] += 10;
1009       B.push_back(p);
1010       // Margins ?
1011       A[i][1] += margin;
1012       B[i][1] += margin;
1013     }
1014
1015   }
1016   //--------------------------------------------------------------------
1017
1018
1019 } // end of namespace
1020