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