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