itk/clitkSegmentationUtils.txx

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