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 #include <itkImageDuplicator.h>
  37 #include <itkSignedMaurerDistanceMapImageFilter.h>
  38
  39 namespace clitk {
  40
  41   //--------------------------------------------------------------------
  42   template<class ImageType>
  43   void ComputeBBFromImageRegion(const ImageType * image,
  44                                 typename ImageType::RegionType region,
  45                                 typename itk::BoundingBox<unsigned long,
  46                                 ImageType::ImageDimension>::Pointer bb) {
  47     typedef typename ImageType::IndexType IndexType;
  48     IndexType firstIndex;
  49     IndexType lastIndex;
  50     for(unsigned int i=0; i<image->GetImageDimension(); i++) {
  51       firstIndex[i] = region.GetIndex()[i];
  52       lastIndex[i] = firstIndex[i]+region.GetSize()[i];
  53     }
  54
  55     typedef itk::BoundingBox<unsigned long,
  56                              ImageType::ImageDimension> BBType;
  57     typedef typename BBType::PointType PointType;
  58     PointType lastPoint;
  59     PointType firstPoint;
  60     image->TransformIndexToPhysicalPoint(firstIndex, firstPoint);
  61     image->TransformIndexToPhysicalPoint(lastIndex, lastPoint);
  62
  63     bb->SetMaximum(lastPoint);
  64     bb->SetMinimum(firstPoint);
  65   }
  66   //--------------------------------------------------------------------
  67
  68
  69   //--------------------------------------------------------------------
  70   template<int Dimension>
  71   void ComputeBBIntersection(typename itk::BoundingBox<unsigned long, Dimension>::Pointer bbo,
  72                              typename itk::BoundingBox<unsigned long, Dimension>::Pointer bbi1,
  73                              typename itk::BoundingBox<unsigned long, Dimension>::Pointer bbi2) {
  74
  75     typedef itk::BoundingBox<unsigned long, Dimension> BBType;
  76     typedef typename BBType::PointType PointType;
  77     PointType lastPoint;
  78     PointType firstPoint;
  79
  80     for(unsigned int i=0; i<Dimension; i++) {
  81       firstPoint[i] = std::max(bbi1->GetMinimum()[i],
  82                                bbi2->GetMinimum()[i]);
  83       lastPoint[i] = std::min(bbi1->GetMaximum()[i],
  84                               bbi2->GetMaximum()[i]);
  85     }
  86
  87     bbo->SetMaximum(lastPoint);
  88     bbo->SetMinimum(firstPoint);
  89   }
  90   //--------------------------------------------------------------------
  91
  92
  93   //--------------------------------------------------------------------
  94   template<class ImageType>
  95   void ComputeRegionFromBB(const ImageType * image,
  96                            const typename itk::BoundingBox<unsigned long,
  97                                                            ImageType::ImageDimension>::Pointer bb,
  98                            typename ImageType::RegionType & region) {
  99     // Types
 100     typedef typename ImageType::IndexType  IndexType;
 101     typedef typename ImageType::PointType  PointType;
 102     typedef typename ImageType::RegionType RegionType;
 103     typedef typename ImageType::SizeType   SizeType;
 104
 105     // Region starting point
 106     IndexType regionStart;
 107     PointType start = bb->GetMinimum();
 108     image->TransformPhysicalPointToIndex(start, regionStart);
 109
 110     // Region size
 111     SizeType regionSize;
 112     PointType maxs = bb->GetMaximum();
 113     PointType mins = bb->GetMinimum();
 114     for(unsigned int i=0; i<ImageType::ImageDimension; i++) {
 115       regionSize[i] = lrint((maxs[i] - mins[i])/image->GetSpacing()[i]);
 116     }
 117
 118     // Create region
 119     region.SetIndex(regionStart);
 120     region.SetSize(regionSize);
 121   }
 122   //--------------------------------------------------------------------
 123
 124   //--------------------------------------------------------------------
 125   template<class ImageType, class TMaskImageType>
 126   typename ImageType::Pointer
 127   SetBackground(const ImageType * input,
 128                 const TMaskImageType * mask,
 129                 typename TMaskImageType::PixelType maskBG,
 130                 typename ImageType::PixelType outValue,
 131                 bool inPlace) {
 132     typedef SetBackgroundImageFilter<ImageType, TMaskImageType, ImageType>
 133       SetBackgroundImageFilterType;
 134     typename SetBackgroundImageFilterType::Pointer setBackgroundFilter
 135       = SetBackgroundImageFilterType::New();
 136     //  if (inPlace) setBackgroundFilter->ReleaseDataFlagOn(); // No seg fault
 137     setBackgroundFilter->SetInPlace(inPlace); // This is important to keep memory low
 138     setBackgroundFilter->SetInput(input);
 139     setBackgroundFilter->SetInput2(mask);
 140     setBackgroundFilter->SetMaskValue(maskBG);
 141     setBackgroundFilter->SetOutsideValue(outValue);
 142     setBackgroundFilter->Update();
 143     return setBackgroundFilter->GetOutput();
 144   }
 145   //--------------------------------------------------------------------
 146
 147
 148   //--------------------------------------------------------------------
 149   template<class ImageType>
 150   int GetNumberOfConnectedComponentLabels(const ImageType * input,
 151                                           typename ImageType::PixelType BG,
 152                                           bool isFullyConnected) {
 153     // Connected Component label
 154     typedef itk::ConnectedComponentImageFilter<ImageType, ImageType> ConnectFilterType;
 155     typename ConnectFilterType::Pointer connectFilter = ConnectFilterType::New();
 156     connectFilter->SetInput(input);
 157     connectFilter->SetBackgroundValue(BG);
 158     connectFilter->SetFullyConnected(isFullyConnected);
 159     connectFilter->Update();
 160
 161     // Return result
 162     return connectFilter->GetObjectCount();
 163   }
 164   //--------------------------------------------------------------------
 165
 166   //--------------------------------------------------------------------
 167   /*
 168     Warning : in this cas, we consider outputType like inputType, not
 169     InternalImageType. Be sure it fits.
 170   */
 171   template<class ImageType>
 172   typename ImageType::Pointer
 173   Labelize(const ImageType * input,
 174            typename ImageType::PixelType BG,
 175            bool isFullyConnected,
 176            int minimalComponentSize) {
 177     // InternalImageType for storing large number of component
 178     typedef itk::Image<int, ImageType::ImageDimension> InternalImageType;
 179
 180     // Connected Component label
 181     typedef itk::ConnectedComponentImageFilter<ImageType, InternalImageType> ConnectFilterType;
 182     typename ConnectFilterType::Pointer connectFilter = ConnectFilterType::New();
 183     //  connectFilter->ReleaseDataFlagOn();
 184     connectFilter->SetInput(input);
 185     connectFilter->SetBackgroundValue(BG);
 186     connectFilter->SetFullyConnected(isFullyConnected);
 187
 188     // Sort by size and remove too small area.
 189     typedef itk::RelabelComponentImageFilter<InternalImageType, ImageType> RelabelFilterType;
 190     typename RelabelFilterType::Pointer relabelFilter = RelabelFilterType::New();
 191     //  relabelFilter->ReleaseDataFlagOn(); // if yes, fail when ExplosionControlledThresholdConnectedImageFilter ???
 192     relabelFilter->SetInput(connectFilter->GetOutput());
 193     relabelFilter->SetMinimumObjectSize(minimalComponentSize);
 194     relabelFilter->Update();
 195
 196     // Return result
 197     typename ImageType::Pointer output = relabelFilter->GetOutput();
 198     return output;
 199   }
 200   //--------------------------------------------------------------------
 201
 202
 203   //--------------------------------------------------------------------
 204   /*
 205     Warning : in this cas, we consider outputType like inputType, not
 206     InternalImageType. Be sure it fits.
 207   */
 208   template<class ImageType>
 209   typename ImageType::Pointer
 210   LabelizeAndCountNumberOfObjects(const ImageType * input,
 211                                   typename ImageType::PixelType BG,
 212                                   bool isFullyConnected,
 213                                   int minimalComponentSize,
 214                                   int & nb) {
 215     // InternalImageType for storing large number of component
 216     typedef itk::Image<int, ImageType::ImageDimension> InternalImageType;
 217
 218     // Connected Component label
 219     typedef itk::ConnectedComponentImageFilter<ImageType, InternalImageType> ConnectFilterType;
 220     typename ConnectFilterType::Pointer connectFilter = ConnectFilterType::New();
 221     //  connectFilter->ReleaseDataFlagOn();
 222     connectFilter->SetInput(input);
 223     connectFilter->SetBackgroundValue(BG);
 224     connectFilter->SetFullyConnected(isFullyConnected);
 225
 226     // Sort by size and remove too small area.
 227     typedef itk::RelabelComponentImageFilter<InternalImageType, ImageType> RelabelFilterType;
 228     typename RelabelFilterType::Pointer relabelFilter = RelabelFilterType::New();
 229     //  relabelFilter->ReleaseDataFlagOn(); // if yes, fail when ExplosionControlledThresholdConnectedImageFilter ???
 230     relabelFilter->SetInput(connectFilter->GetOutput());
 231     relabelFilter->SetMinimumObjectSize(minimalComponentSize);
 232     relabelFilter->Update();
 233
 234     nb = relabelFilter->GetNumberOfObjects();
 235     // DD(relabelFilter->GetOriginalNumberOfObjects());
 236     // DD(relabelFilter->GetSizeOfObjectsInPhysicalUnits()[0]);
 237
 238     // Return result
 239     typename ImageType::Pointer output = relabelFilter->GetOutput();
 240     return output;
 241   }
 242   //--------------------------------------------------------------------
 243
 244
 245
 246   //--------------------------------------------------------------------
 247   template<class ImageType>
 248   typename ImageType::Pointer
 249   RemoveLabels(const ImageType * input,
 250                typename ImageType::PixelType BG,
 251                std::vector<typename ImageType::PixelType> & labelsToRemove) {
 252     assert(labelsToRemove.size() != 0);
 253     typename ImageType::Pointer working_image;// = input;
 254     for (unsigned int i=0; i <labelsToRemove.size(); i++) {
 255       typedef SetBackgroundImageFilter<ImageType, ImageType> SetBackgroundImageFilterType;
 256       typename SetBackgroundImageFilterType::Pointer setBackgroundFilter = SetBackgroundImageFilterType::New();
 257       setBackgroundFilter->SetInput(input);
 258       setBackgroundFilter->SetInput2(input);
 259       setBackgroundFilter->SetMaskValue(labelsToRemove[i]);
 260       setBackgroundFilter->SetOutsideValue(BG);
 261       setBackgroundFilter->Update();
 262       working_image = setBackgroundFilter->GetOutput();
 263     }
 264     return working_image;
 265   }
 266   //--------------------------------------------------------------------
 267
 268
 269   //--------------------------------------------------------------------
 270   template<class ImageType>
 271   typename ImageType::Pointer
 272   KeepLabels(const ImageType * input,
 273              typename ImageType::PixelType BG,
 274              typename ImageType::PixelType FG,
 275              typename ImageType::PixelType firstKeep,
 276              typename ImageType::PixelType lastKeep,
 277              bool useLastKeep) {
 278     typedef itk::BinaryThresholdImageFilter<ImageType, ImageType> BinarizeFilterType;
 279     typename BinarizeFilterType::Pointer binarizeFilter = BinarizeFilterType::New();
 280     binarizeFilter->SetInput(input);
 281     binarizeFilter->SetLowerThreshold(firstKeep);
 282     if (useLastKeep) binarizeFilter->SetUpperThreshold(lastKeep);
 283     binarizeFilter->SetInsideValue(FG);
 284     binarizeFilter->SetOutsideValue(BG);
 285     binarizeFilter->Update();
 286     return binarizeFilter->GetOutput();
 287   }
 288   //--------------------------------------------------------------------
 289
 290
 291   //--------------------------------------------------------------------
 292   template<class ImageType>
 293   typename ImageType::Pointer
 294   LabelizeAndSelectLabels(const ImageType * input,
 295                           typename ImageType::PixelType BG,
 296                           typename ImageType::PixelType FG,
 297                           bool isFullyConnected,
 298                           int minimalComponentSize,
 299                           LabelizeParameters<typename ImageType::PixelType> * param)
 300   {
 301     typename ImageType::Pointer working_image;
 302     working_image = Labelize<ImageType>(input, BG, isFullyConnected, minimalComponentSize);
 303     if (param->GetLabelsToRemove().size() != 0)
 304       working_image = RemoveLabels<ImageType>(working_image, BG, param->GetLabelsToRemove());
 305     working_image = KeepLabels<ImageType>(working_image,
 306                                           BG, FG,
 307                                           param->GetFirstKeep(),
 308                                           param->GetLastKeep(),
 309                                           param->GetUseLastKeep());
 310     return working_image;
 311   }
 312   //--------------------------------------------------------------------
 313
 314
 315   //--------------------------------------------------------------------
 316   template<class ImageType>
 317   typename ImageType::Pointer
 318   ResizeImageLike(const ImageType * input,
 319                   const itk::ImageBase<ImageType::ImageDimension> * like,
 320                   typename ImageType::PixelType backgroundValue)
 321   {
 322     typedef CropLikeImageFilter<ImageType> CropFilterType;
 323     typename CropFilterType::Pointer cropFilter = CropFilterType::New();
 324     cropFilter->SetInput(input);
 325     cropFilter->SetCropLikeImage(like);
 326     cropFilter->SetBackgroundValue(backgroundValue);
 327     cropFilter->Update();
 328     return cropFilter->GetOutput();
 329   }
 330   //--------------------------------------------------------------------
 331
 332
 333   //--------------------------------------------------------------------
 334   template<class MaskImageType>
 335   typename MaskImageType::Pointer
 336   SliceBySliceRelativePosition(const MaskImageType * input,
 337                                const MaskImageType * object,
 338                                int direction,
 339                                double threshold,
 340                                std::string orientation,
 341                                bool uniqueConnectedComponent,
 342                                double spacing,
 343                                bool autocropFlag,
 344                                bool singleObjectCCL)
 345   {
 346     typedef SliceBySliceRelativePositionFilter<MaskImageType> SliceRelPosFilterType;
 347     typename SliceRelPosFilterType::Pointer sliceRelPosFilter = SliceRelPosFilterType::New();
 348     sliceRelPosFilter->VerboseStepFlagOff();
 349     sliceRelPosFilter->WriteStepFlagOff();
 350     sliceRelPosFilter->SetInput(input);
 351     sliceRelPosFilter->SetInputObject(object);
 352     sliceRelPosFilter->SetDirection(direction);
 353     sliceRelPosFilter->SetFuzzyThreshold(threshold);
 354     sliceRelPosFilter->AddOrientationTypeString(orientation);
 355     sliceRelPosFilter->SetIntermediateSpacingFlag((spacing != -1));
 356     sliceRelPosFilter->SetIntermediateSpacing(spacing);
 357     sliceRelPosFilter->SetUniqueConnectedComponentBySliceFlag(uniqueConnectedComponent);
 358     sliceRelPosFilter->ObjectCCLSelectionFlagOff();
 359     sliceRelPosFilter->SetUseTheLargestObjectCCLFlag(singleObjectCCL);
 360     //    sliceRelPosFilter->SetInverseOrientationFlag(inverseflag);
 361     sliceRelPosFilter->SetAutoCropFlag(autocropFlag);
 362     sliceRelPosFilter->IgnoreEmptySliceObjectFlagOn();
 363     sliceRelPosFilter->Update();
 364     return sliceRelPosFilter->GetOutput();
 365   }
 366   //--------------------------------------------------------------------
 367
 368
 369   //--------------------------------------------------------------------
 370   template<class ImageType>
 371   bool
 372   FindExtremaPointInAGivenDirection(const ImageType * input,
 373                                     typename ImageType::PixelType bg,
 374                                     int direction, bool opposite,
 375                                     typename ImageType::PointType & point)
 376   {
 377     typename ImageType::PointType dummy;
 378     return FindExtremaPointInAGivenDirection(input, bg, direction,
 379                                              opposite, dummy, 0, point);
 380   }
 381   //--------------------------------------------------------------------
 382
 383
 384   //--------------------------------------------------------------------
 385   template<class ImageType>
 386   bool
 387   FindExtremaPointInAGivenDirection(const ImageType * input,
 388                                     typename ImageType::PixelType bg,
 389                                     int direction, bool opposite,
 390                                     typename ImageType::PointType refpoint,
 391                                     double distanceMax,
 392                                     typename ImageType::PointType & point)
 393   {
 394     /*
 395       loop over input pixels, store the index in the fg that is max
 396       according to the given direction.
 397     */
 398     typedef itk::ImageRegionConstIteratorWithIndex<ImageType> IteratorType;
 399     IteratorType iter(input, input->GetLargestPossibleRegion());
 400     iter.GoToBegin();
 401     typename ImageType::IndexType max = input->GetLargestPossibleRegion().GetIndex();
 402     if (opposite) max = max+input->GetLargestPossibleRegion().GetSize();
 403     bool found=false;
 404     while (!iter.IsAtEnd()) {
 405       if (iter.Get() != bg) {
 406         bool test = iter.GetIndex()[direction] >  max[direction];
 407         if (opposite) test = !test;
 408         if (test) {
 409           typename ImageType::PointType p;
 410           input->TransformIndexToPhysicalPoint(iter.GetIndex(), p);
 411           if ((distanceMax==0) || (p.EuclideanDistanceTo(refpoint) < distanceMax)) {
 412             max = iter.GetIndex();
 413             found = true;
 414           }
 415         }
 416       }
 417       ++iter;
 418     }
 419     if (!found) return false;
 420     input->TransformIndexToPhysicalPoint(max, point);
 421     return true;
 422   }
 423   //--------------------------------------------------------------------
 424
 425
 426   //--------------------------------------------------------------------
 427   template<class ImageType>
 428   typename ImageType::Pointer
 429   CropImageRemoveGreaterThan(const ImageType * image,
 430                  int dim, double min, bool autoCrop,
 431                  typename ImageType::PixelType BG)
 432   {
 433     return CropImageAlongOneAxis<ImageType>(image, dim,
 434                                             image->GetOrigin()[dim],
 435                                             min,
 436                                             autoCrop, BG);
 437   }
 438   //--------------------------------------------------------------------
 439
 440
 441   //--------------------------------------------------------------------
 442   template<class ImageType>
 443   typename ImageType::Pointer
 444   CropImageRemoveLowerThan(const ImageType * image,
 445                  int dim, double max, bool autoCrop,
 446                  typename ImageType::PixelType BG)
 447   {
 448     typename ImageType::PointType p;
 449     image->TransformIndexToPhysicalPoint(image->GetLargestPossibleRegion().GetIndex()+
 450                                          image->GetLargestPossibleRegion().GetSize(), p);
 451     return CropImageAlongOneAxis<ImageType>(image, dim, max, p[dim], autoCrop, BG);
 452   }
 453   //--------------------------------------------------------------------
 454
 455
 456   //--------------------------------------------------------------------
 457   template<class ImageType>
 458   typename ImageType::Pointer
 459   CropImageAlongOneAxis(const ImageType * image,
 460                         int dim, double min, double max,
 461                         bool autoCrop, typename ImageType::PixelType BG)
 462   {
 463     // Compute region size
 464     typename ImageType::RegionType region;
 465     typename ImageType::SizeType size = image->GetLargestPossibleRegion().GetSize();
 466     typename ImageType::PointType p = image->GetOrigin();
 467     if (min > p[dim]) p[dim] = min; // Check if not outside the image
 468     typename ImageType::IndexType start;
 469     image->TransformPhysicalPointToIndex(p, start);
 470     double m = image->GetOrigin()[dim] + size[dim]*image->GetSpacing()[dim];
 471     if (max > m) p[dim] = m; // Check if not outside the image
 472     else p[dim] = max;
 473     typename ImageType::IndexType end;
 474     image->TransformPhysicalPointToIndex(p, end);
 475     size[dim] = abs(end[dim]-start[dim]);
 476     region.SetIndex(start);
 477     region.SetSize(size);
 478
 479     // Perform Crop
 480     typedef itk::RegionOfInterestImageFilter<ImageType, ImageType> CropFilterType;
 481     typename CropFilterType::Pointer cropFilter = CropFilterType::New();
 482     cropFilter->SetInput(image);
 483     cropFilter->SetRegionOfInterest(region);
 484     cropFilter->Update();
 485     typename ImageType::Pointer result = cropFilter->GetOutput();
 486
 487     // Auto Crop
 488     if (autoCrop) {
 489       result = AutoCrop<ImageType>(result, BG);
 490     }
 491     return result;
 492   }
 493   //--------------------------------------------------------------------
 494
 495
 496   //--------------------------------------------------------------------
 497   template<class ImageType>
 498   void
 499   ComputeCentroids(const ImageType * image,
 500                    typename ImageType::PixelType BG,
 501                    std::vector<typename ImageType::PointType> & centroids)
 502   {
 503     typedef long LabelType;
 504     static const unsigned int Dim = ImageType::ImageDimension;
 505     typedef itk::ShapeLabelObject< LabelType, Dim > LabelObjectType;
 506     typedef itk::LabelMap< LabelObjectType > LabelMapType;
 507     typedef itk::LabelImageToLabelMapFilter<ImageType, LabelMapType> ImageToMapFilterType;
 508     typename ImageToMapFilterType::Pointer imageToLabelFilter = ImageToMapFilterType::New();
 509     typedef itk::ShapeLabelMapFilter<LabelMapType, ImageType> ShapeFilterType;
 510     typename ShapeFilterType::Pointer statFilter = ShapeFilterType::New();
 511     imageToLabelFilter->SetBackgroundValue(BG);
 512     imageToLabelFilter->SetInput(image);
 513     statFilter->SetInput(imageToLabelFilter->GetOutput());
 514     statFilter->Update();
 515     typename LabelMapType::Pointer labelMap = statFilter->GetOutput();
 516
 517     centroids.clear();
 518     typename ImageType::PointType dummy;
 519     centroids.push_back(dummy); // label 0 -> no centroid, use dummy point for BG
 520     //DS FIXME (not useful ! to change ..)
 521     for(uint i=0; i<labelMap->GetNumberOfLabelObjects(); i++) {
 522       int label = labelMap->GetLabels()[i];
 523       centroids.push_back(labelMap->GetLabelObject(label)->GetCentroid());
 524     }
 525   }
 526   //--------------------------------------------------------------------
 527
 528
 529   //--------------------------------------------------------------------
 530   template<class ImageType, class LabelType>
 531   typename itk::LabelMap< itk::ShapeLabelObject<LabelType, ImageType::ImageDimension> >::Pointer
 532   ComputeLabelMap(const ImageType * image,
 533                   typename ImageType::PixelType BG,
 534                   bool computePerimeterFlag)
 535   {
 536     static const unsigned int Dim = ImageType::ImageDimension;
 537     typedef itk::ShapeLabelObject< LabelType, Dim > LabelObjectType;
 538     typedef itk::LabelMap< LabelObjectType > LabelMapType;
 539     typedef itk::LabelImageToLabelMapFilter<ImageType, LabelMapType> ImageToMapFilterType;
 540     typename ImageToMapFilterType::Pointer imageToLabelFilter = ImageToMapFilterType::New();
 541     typedef itk::ShapeLabelMapFilter<LabelMapType, ImageType> ShapeFilterType;
 542     typename ShapeFilterType::Pointer statFilter = ShapeFilterType::New();
 543     imageToLabelFilter->SetBackgroundValue(BG);
 544     imageToLabelFilter->SetInput(image);
 545     statFilter->SetInput(imageToLabelFilter->GetOutput());
 546     statFilter->SetComputePerimeter(computePerimeterFlag);
 547     statFilter->Update();
 548     return statFilter->GetOutput();
 549   }
 550   //--------------------------------------------------------------------
 551
 552
 553   //--------------------------------------------------------------------
 554   template<class ImageType>
 555   void
 556   ComputeCentroids2(const ImageType * image,
 557                    typename ImageType::PixelType BG,
 558                    std::vector<typename ImageType::PointType> & centroids)
 559   {
 560     typedef long LabelType;
 561     static const unsigned int Dim = ImageType::ImageDimension;
 562     typedef itk::ShapeLabelObject< LabelType, Dim > LabelObjectType;
 563     typedef itk::LabelMap< LabelObjectType > LabelMapType;
 564     typedef itk::LabelImageToLabelMapFilter<ImageType, LabelMapType> ImageToMapFilterType;
 565     typename ImageToMapFilterType::Pointer imageToLabelFilter = ImageToMapFilterType::New();
 566     typedef itk::ShapeLabelMapFilter<LabelMapType, ImageType> ShapeFilterType;
 567     typename ShapeFilterType::Pointer statFilter = ShapeFilterType::New();
 568     imageToLabelFilter->SetBackgroundValue(BG);
 569     imageToLabelFilter->SetInput(image);
 570     statFilter->SetInput(imageToLabelFilter->GetOutput());
 571     statFilter->Update();
 572     typename LabelMapType::Pointer labelMap = statFilter->GetOutput();
 573
 574     centroids.clear();
 575     typename ImageType::PointType dummy;
 576     centroids.push_back(dummy); // label 0 -> no centroid, use dummy point
 577     for(uint i=1; i<labelMap->GetNumberOfLabelObjects()+1; i++) {
 578       centroids.push_back(labelMap->GetLabelObject(i)->GetCentroid());
 579     }
 580
 581     for(uint i=1; i<labelMap->GetNumberOfLabelObjects()+1; i++) {
 582       DD(labelMap->GetLabelObject(i)->GetBinaryPrincipalAxes());
 583       DD(labelMap->GetLabelObject(i)->GetBinaryFlatness());
 584       DD(labelMap->GetLabelObject(i)->GetRoundness ());
 585
 586       // search for the point on the boundary alog PA
 587
 588     }
 589
 590   }
 591   //--------------------------------------------------------------------
 592
 593
 594   //--------------------------------------------------------------------
 595   template<class ImageType>
 596   void
 597   PointsUtils<ImageType>::Convert2DTo3D(const PointType2D & p2D,
 598                                         const ImageType * image,
 599                                         const int slice,
 600                                         PointType3D & p3D)
 601   {
 602     IndexType3D index3D;
 603     index3D[0] = index3D[1] = 0;
 604     index3D[2] = image->GetLargestPossibleRegion().GetIndex()[2]+slice;
 605     image->TransformIndexToPhysicalPoint(index3D, p3D);
 606     p3D[0] = p2D[0];
 607     p3D[1] = p2D[1];
 608     //  p3D[2] = p[2];//(image->GetLargestPossibleRegion().GetIndex()[2]+slice)*image->GetSpacing()[2]
 609     //    + image->GetOrigin()[2];
 610   }
 611   //--------------------------------------------------------------------
 612
 613
 614   //--------------------------------------------------------------------
 615   template<class ImageType>
 616   void
 617   PointsUtils<ImageType>::Convert2DMapTo3DList(const MapPoint2DType & map,
 618                                             const ImageType * image,
 619                                             VectorPoint3DType & list)
 620   {
 621     typename MapPoint2DType::const_iterator iter = map.begin();
 622     while (iter != map.end()) {
 623       PointType3D p;
 624       Convert2DTo3D(iter->second, image, iter->first, p);
 625       list.push_back(p);
 626       ++iter;
 627     }
 628   }
 629   //--------------------------------------------------------------------
 630
 631
 632   //--------------------------------------------------------------------
 633   template<class ImageType>
 634   void
 635   PointsUtils<ImageType>::Convert2DListTo3DList(const VectorPoint2DType & p2D,
 636                                                 int slice,
 637                                                 const ImageType * image,
 638                                                 VectorPoint3DType & list)
 639   {
 640     for(uint i=0; i<p2D.size(); i++) {
 641       PointType3D p;
 642       Convert2DTo3D(p2D[i], image, slice, p);
 643       list.push_back(p);
 644     }
 645   }
 646   //--------------------------------------------------------------------
 647
 648
 649   //--------------------------------------------------------------------
 650   template<class ImageType>
 651   void
 652   WriteListOfLandmarks(std::vector<typename ImageType::PointType> points,
 653                        std::string filename)
 654   {
 655     std::ofstream os;
 656     openFileForWriting(os, filename);
 657     os << "LANDMARKS1" << std::endl;
 658     for(uint i=0; i<points.size(); i++) {
 659       const typename ImageType::PointType & p = points[i];
 660       // Write it in the file
 661       os << i << " " << p[0] << " " << p[1] << " " << p[2] << " 0 0 " << std::endl;
 662     }
 663     os.close();
 664   }
 665   //--------------------------------------------------------------------
 666
 667
 668   //--------------------------------------------------------------------
 669   template<class ImageType>
 670   typename ImageType::Pointer
 671   Dilate(const ImageType * image, double radiusInMM,
 672          typename ImageType::PixelType BG,
 673          typename ImageType::PixelType FG,
 674          bool extendSupport)
 675   {
 676     typename ImageType::SizeType r;
 677     for(uint i=0; i<ImageType::ImageDimension; i++)
 678       r[i] = (uint)lrint(radiusInMM/image->GetSpacing()[i]);
 679     return Dilate<ImageType>(image, r, BG, FG, extendSupport);
 680   }
 681   //--------------------------------------------------------------------
 682
 683
 684   //--------------------------------------------------------------------
 685   template<class ImageType>
 686   typename ImageType::Pointer
 687   Dilate(const ImageType * image, typename ImageType::PointType radiusInMM,
 688          typename ImageType::PixelType BG,
 689          typename ImageType::PixelType FG,
 690          bool extendSupport)
 691   {
 692     typename ImageType::SizeType r;
 693     for(uint i=0; i<ImageType::ImageDimension; i++)
 694       r[i] = (uint)lrint(radiusInMM[i]/image->GetSpacing()[i]);
 695     return Dilate<ImageType>(image, r, BG, FG, extendSupport);
 696   }
 697   //--------------------------------------------------------------------
 698
 699
 700   //--------------------------------------------------------------------
 701   template<class ImageType>
 702   typename ImageType::Pointer
 703   Dilate(const ImageType * image, typename ImageType::SizeType radius,
 704          typename ImageType::PixelType BG,
 705          typename ImageType::PixelType FG,
 706          bool extendSupport)
 707   {
 708     // Create kernel for dilatation
 709     typedef itk::BinaryBallStructuringElement<typename ImageType::PixelType,
 710                                               ImageType::ImageDimension> KernelType;
 711     KernelType structuringElement;
 712     structuringElement.SetRadius(radius);
 713     structuringElement.CreateStructuringElement();
 714
 715     typename ImageType::Pointer output;
 716     if (extendSupport) {
 717       typedef itk::ConstantPadImageFilter<ImageType, ImageType> PadFilterType;
 718       typename PadFilterType::Pointer padFilter = PadFilterType::New();
 719       padFilter->SetInput(image);
 720       typename ImageType::SizeType lower;
 721       typename ImageType::SizeType upper;
 722       for(uint i=0; i<3; i++) {
 723         lower[i] = upper[i] = 2*(radius[i]+1);
 724       }
 725       padFilter->SetPadLowerBound(lower);
 726       padFilter->SetPadUpperBound(upper);
 727       padFilter->Update();
 728       output = padFilter->GetOutput();
 729     }
 730
 731     // Dilate  filter
 732     typedef itk::BinaryDilateImageFilter<ImageType, ImageType , KernelType> DilateFilterType;
 733     typename DilateFilterType::Pointer dilateFilter = DilateFilterType::New();
 734     dilateFilter->SetBackgroundValue(BG);
 735     dilateFilter->SetForegroundValue(FG);
 736     dilateFilter->SetBoundaryToForeground(false);
 737     dilateFilter->SetKernel(structuringElement);
 738     if (extendSupport) dilateFilter->SetInput(output);
 739     else dilateFilter->SetInput(image);
 740     dilateFilter->Update();
 741     return dilateFilter->GetOutput();
 742   }
 743   //--------------------------------------------------------------------
 744
 745
 746   //--------------------------------------------------------------------
 747   template<class ImageType>
 748   typename ImageType::Pointer
 749   Opening(const ImageType * image, typename ImageType::SizeType radius,
 750          typename ImageType::PixelType BG,
 751          typename ImageType::PixelType FG)
 752   {
 753     // Kernel
 754     typedef itk::BinaryBallStructuringElement<typename ImageType::PixelType,
 755                                               ImageType::ImageDimension> KernelType;
 756     KernelType structuringElement;
 757     structuringElement.SetRadius(radius);
 758     structuringElement.CreateStructuringElement();
 759
 760     // Filter
 761     typedef itk::BinaryMorphologicalOpeningImageFilter<ImageType, ImageType , KernelType> OpeningFilterType;
 762     typename OpeningFilterType::Pointer open = OpeningFilterType::New();
 763     open->SetInput(image);
 764     open->SetBackgroundValue(BG);
 765     open->SetForegroundValue(FG);
 766     open->SetKernel(structuringElement);
 767     open->Update();
 768     return open->GetOutput();
 769   }
 770   //--------------------------------------------------------------------
 771
 772
 773
 774   //--------------------------------------------------------------------
 775   template<class ValueType, class VectorType>
 776   void ConvertOption(std::string optionName, uint given,
 777                      ValueType * values, VectorType & p,
 778                      uint dim, bool required)
 779   {
 780     if (required && (given == 0)) {
 781       clitkExceptionMacro("The option --" << optionName << " must be set and have 1 or "
 782                           << dim << " values.");
 783     }
 784     if (given == 1) {
 785       for(uint i=0; i<dim; i++) p[i] = values[0];
 786       return;
 787     }
 788     if (given == dim) {
 789       for(uint i=0; i<dim; i++) p[i] = values[i];
 790       return;
 791     }
 792     if (given == 0) return;
 793     clitkExceptionMacro("The option --" << optionName << " must have 1 or "
 794                         << dim << " values.");
 795   }
 796   //--------------------------------------------------------------------
 797
 798
 799   //--------------------------------------------------------------------
 800   /*
 801     http://www.gamedev.net/community/forums/topic.asp?topic_id=542870
 802     Assuming the points are (Ax,Ay) (Bx,By) and (Cx,Cy), you need to compute:
 803     (Bx - Ax) * (Cy - Ay) - (By - Ay) * (Cx - Ax)
 804     This will equal zero if the point C is on the line formed by
 805     points A and B, and will have a different sign depending on the
 806     side. Which side this is depends on the orientation of your (x,y)
 807     coordinates, but you can plug test values for A,B and C into this
 808     formula to determine whether negative values are to the left or to
 809     the right.
 810     => to accelerate, start with formula, when change sign -> stop and fill
 811
 812     offsetToKeep = is used to determine which side of the line we
 813     keep. The point along the mainDirection but 'offsetToKeep' mm away
 814     is kept.
 815
 816   */
 817   template<class ImageType>
 818   void
 819   SliceBySliceSetBackgroundFromLineSeparation(ImageType * input,
 820                                               std::vector<typename ImageType::PointType> & lA,
 821                                               std::vector<typename ImageType::PointType> & lB,
 822                                               typename ImageType::PixelType BG,
 823                                               int mainDirection,
 824                                               double offsetToKeep)
 825   {
 826     assert((mainDirection==0) || (mainDirection==1));
 827     typedef itk::ImageSliceIteratorWithIndex<ImageType> SliceIteratorType;
 828     SliceIteratorType siter = SliceIteratorType(input,
 829                                                 input->GetLargestPossibleRegion());
 830     siter.SetFirstDirection(0);
 831     siter.SetSecondDirection(1);
 832     siter.GoToBegin();
 833     uint i=0;
 834     typename ImageType::PointType A;
 835     typename ImageType::PointType B;
 836     typename ImageType::PointType C;
 837     assert(lA.size() == lB.size());
 838     while ((i<lA.size()) && (!siter.IsAtEnd())) {
 839       // Check that the current slice correspond to the current point
 840       input->TransformIndexToPhysicalPoint(siter.GetIndex(), C);
 841       if ((fabs(C[2] - lA[i][2]))>0.01) { // is !equal with a tolerance of 0.01 mm
 842       }
 843       else {
 844         // Define A,B,C points
 845         A = lA[i];
 846         B = lB[i];
 847         C = A;
 848
 849         // Check that the line is not a point (A=B)
 850         bool p = (A[0] == B[0]) && (A[1] == B[1]);
 851
 852         if (!p) {
 853           C[mainDirection] += offsetToKeep; // I know I must keep this point
 854           double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
 855           bool isPositive = s<0;
 856           while (!siter.IsAtEndOfSlice()) {
 857             while (!siter.IsAtEndOfLine()) {
 858               // Very slow, I know ... but image should be very small
 859               input->TransformIndexToPhysicalPoint(siter.GetIndex(), C);
 860               double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
 861               if (s == 0) siter.Set(BG); // on the line, we decide to remove
 862               if (isPositive) {
 863                 if (s > 0) siter.Set(BG);
 864               }
 865               else {
 866                 if (s < 0) siter.Set(BG);
 867               }
 868               ++siter;
 869             }
 870             siter.NextLine();
 871           } // end loop slice
 872         }
 873
 874         ++i;
 875       } // End of current slice
 876       siter.NextSlice();
 877     }
 878   }
 879   //--------------------------------------------------------------------
 880
 881
 882   //--------------------------------------------------------------------
 883   template<class ImageType>
 884   void
 885   AndNot(ImageType * input,
 886          const ImageType * object,
 887          typename ImageType::PixelType BG)
 888   {
 889     typename ImageType::Pointer o;
 890     bool resized=false;
 891     if (!clitk::HaveSameSizeAndSpacing<ImageType, ImageType>(input, object)) {
 892       o = clitk::ResizeImageLike<ImageType>(object, input, BG);
 893       resized = true;
 894     }
 895
 896     typedef clitk::BooleanOperatorLabelImageFilter<ImageType> BoolFilterType;
 897     typename BoolFilterType::Pointer boolFilter = BoolFilterType::New();
 898     boolFilter->InPlaceOn();
 899     boolFilter->SetInput1(input);
 900     if (resized) boolFilter->SetInput2(o);
 901     else boolFilter->SetInput2(object);
 902     boolFilter->SetBackgroundValue1(BG);
 903     boolFilter->SetBackgroundValue2(BG);
 904     boolFilter->SetOperationType(BoolFilterType::AndNot);
 905     boolFilter->Update();
 906   }
 907   //--------------------------------------------------------------------
 908
 909
 910   //--------------------------------------------------------------------
 911   template<class ImageType>
 912   void
 913   And(ImageType * input,
 914       const ImageType * object,
 915       typename ImageType::PixelType BG)
 916   {
 917     typename ImageType::Pointer o;
 918     bool resized=false;
 919     if (!clitk::HaveSameSizeAndSpacing<ImageType, ImageType>(input, object)) {
 920       o = clitk::ResizeImageLike<ImageType>(object, input, BG);
 921       resized = true;
 922     }
 923
 924     typedef clitk::BooleanOperatorLabelImageFilter<ImageType> BoolFilterType;
 925     typename BoolFilterType::Pointer boolFilter = BoolFilterType::New();
 926     boolFilter->InPlaceOn();
 927     boolFilter->SetInput1(input);
 928     if (resized) boolFilter->SetInput2(o);
 929     else boolFilter->SetInput2(object);
 930     boolFilter->SetBackgroundValue1(BG);
 931     boolFilter->SetBackgroundValue2(BG);
 932     boolFilter->SetOperationType(BoolFilterType::And);
 933     boolFilter->Update();
 934   }
 935   //--------------------------------------------------------------------
 936
 937
 938   //--------------------------------------------------------------------
 939   template<class ImageType>
 940   void
 941   Or(ImageType * input,
 942      const ImageType * object,
 943      typename ImageType::PixelType BG)
 944   {
 945     typename ImageType::Pointer o;
 946     bool resized=false;
 947     if (!clitk::HaveSameSizeAndSpacing<ImageType, ImageType>(input, object)) {
 948       o = clitk::ResizeImageLike<ImageType>(object, input, BG);
 949       resized = true;
 950     }
 951
 952     typedef clitk::BooleanOperatorLabelImageFilter<ImageType> BoolFilterType;
 953     typename BoolFilterType::Pointer boolFilter = BoolFilterType::New();
 954     boolFilter->InPlaceOn();
 955     boolFilter->SetInput1(input);
 956     if (resized) boolFilter->SetInput2(o);
 957     else boolFilter->SetInput2(object);
 958     boolFilter->SetBackgroundValue1(BG);
 959     boolFilter->SetBackgroundValue2(BG);
 960     boolFilter->SetOperationType(BoolFilterType::Or);
 961     boolFilter->Update();
 962   }
 963   //--------------------------------------------------------------------
 964
 965
 966   //--------------------------------------------------------------------
 967   template<class ImageType>
 968   typename ImageType::Pointer
 969   Binarize(const ImageType * input,
 970            typename ImageType::PixelType lower,
 971            typename ImageType::PixelType upper,
 972            typename ImageType::PixelType BG,
 973            typename ImageType::PixelType FG)
 974   {
 975     typedef itk::BinaryThresholdImageFilter<ImageType, ImageType> BinaryThresholdFilterType;
 976     typename BinaryThresholdFilterType::Pointer binarizeFilter = BinaryThresholdFilterType::New();
 977     binarizeFilter->SetInput(input);
 978     binarizeFilter->InPlaceOff();
 979     binarizeFilter->SetLowerThreshold(lower);
 980     binarizeFilter->SetUpperThreshold(upper);
 981     binarizeFilter->SetInsideValue(FG);
 982     binarizeFilter->SetOutsideValue(BG);
 983     binarizeFilter->Update();
 984     return binarizeFilter->GetOutput();
 985   }
 986   //--------------------------------------------------------------------
 987
 988
 989   //--------------------------------------------------------------------
 990   template<class ImageType>
 991   void
 992   GetMinMaxPointPosition(const ImageType * input,
 993                          typename ImageType::PointType & min,
 994                          typename ImageType::PointType & max)
 995   {
 996     typename ImageType::IndexType index = input->GetLargestPossibleRegion().GetIndex();
 997     input->TransformIndexToPhysicalPoint(index, min);
 998     index = index+input->GetLargestPossibleRegion().GetSize();
 999     input->TransformIndexToPhysicalPoint(index, max);
1000   }
1001   //--------------------------------------------------------------------
1002
1003
1004   //--------------------------------------------------------------------
1005   template<class ImageType>
1006   typename ImageType::PointType
1007   FindExtremaPointInAGivenLine(const ImageType * input,
1008                                int dimension,
1009                                bool inverse,
1010                                typename ImageType::PointType p,
1011                                typename ImageType::PixelType BG,
1012                                double distanceMax)
1013   {
1014     // Which direction ?  Increasing or decreasing.
1015     int d=1;
1016     if (inverse) d=-1;
1017
1018     // Transform to pixel index
1019     typename ImageType::IndexType index;
1020     input->TransformPhysicalPointToIndex(p, index);
1021
1022     // Loop while inside the mask;
1023     while (input->GetPixel(index) != BG) {
1024       index[dimension] += d;
1025     }
1026
1027     // Transform back to Physical Units
1028     typename ImageType::PointType result;
1029     input->TransformIndexToPhysicalPoint(index, result);
1030
1031     // Check that is is not too far away
1032     double distance = p.EuclideanDistanceTo(result);
1033     if (distance > distanceMax) {
1034       result = p; // Get back to initial value
1035     }
1036
1037     return result;
1038   }
1039   //--------------------------------------------------------------------
1040
1041
1042   //--------------------------------------------------------------------
1043   template<class PointType>
1044   bool
1045   IsOnTheSameLineSide(PointType C, PointType A, PointType B, PointType like)
1046   {
1047     // Look at the position of point 'like' according to the AB line
1048     double s = (B[0] - A[0]) * (like[1] - A[1]) - (B[1] - A[1]) * (like[0] - A[0]);
1049     bool negative = s<0;
1050
1051     // Look the C position
1052     s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
1053
1054     if (negative && (s<=0)) return true;
1055     if (!negative && (s>=0)) return true;
1056     return false;
1057   }
1058   //--------------------------------------------------------------------
1059
1060
1061   //--------------------------------------------------------------------
1062   /* Consider an input object, for each slice, find the extrema
1063      position according to a given direction and build a line segment
1064      passing throught this point in a given direction.  Output is a
1065      vector of line (from point A to B), for each slice;
1066    */
1067   template<class ImageType>
1068   void
1069   SliceBySliceBuildLineSegmentAccordingToExtremaPosition(const ImageType * input,
1070                                                          typename ImageType::PixelType BG,
1071                                                          int sliceDimension,
1072                                                          int extremaDirection,
1073                                                          bool extremaOppositeFlag,
1074                                                          int lineDirection,
1075                                                          double margin,
1076                                                          std::vector<typename ImageType::PointType> & A,
1077                                                          std::vector<typename ImageType::PointType> & B)
1078   {
1079     // Type of a slice
1080     typedef typename itk::Image<typename ImageType::PixelType, ImageType::ImageDimension-1> SliceType;
1081
1082     // Build the list of slices
1083     std::vector<typename SliceType::Pointer> slices;
1084     clitk::ExtractSlices<ImageType>(input, sliceDimension, slices);
1085
1086     // Build the list of 2D points
1087     std::map<int, typename SliceType::PointType> position2D;
1088     for(uint i=0; i<slices.size(); i++) {
1089       typename SliceType::PointType p;
1090       bool found =
1091         clitk::FindExtremaPointInAGivenDirection<SliceType>(slices[i], BG,
1092                                                             extremaDirection, extremaOppositeFlag, p);
1093       if (found) {
1094         position2D[i] = p;
1095       }
1096     }
1097
1098     // Convert 2D points in slice into 3D points
1099     clitk::PointsUtils<ImageType>::Convert2DMapTo3DList(position2D, input, A);
1100
1101     // Create additional point just right to the previous ones, on the
1102     // given lineDirection, in order to create a horizontal/vertical line.
1103     for(uint i=0; i<A.size(); i++) {
1104       typename ImageType::PointType p = A[i];
1105       p[lineDirection] += 10;
1106       B.push_back(p);
1107       // Margins ?
1108       A[i][extremaDirection] += margin;
1109       B[i][extremaDirection] += margin;
1110     }
1111
1112   }
1113   //--------------------------------------------------------------------
1114
1115
1116   //--------------------------------------------------------------------
1117   template<class ImageType>
1118   typename ImageType::Pointer
1119   SliceBySliceKeepMainCCL(const ImageType * input,
1120                           typename ImageType::PixelType BG,
1121                           typename ImageType::PixelType FG)  {
1122
1123     // Extract slices
1124     const int d = ImageType::ImageDimension-1;
1125     typedef typename itk::Image<typename ImageType::PixelType, d> SliceType;
1126     std::vector<typename SliceType::Pointer> slices;
1127     clitk::ExtractSlices<ImageType>(input, d, slices);
1128
1129     // Labelize and keep the main one
1130     std::vector<typename SliceType::Pointer> o;
1131     for(uint i=0; i<slices.size(); i++) {
1132       o.push_back(clitk::Labelize<SliceType>(slices[i], BG, false, 1));
1133       o[i] = clitk::KeepLabels<SliceType>(o[i], BG, FG, 1, 1, true);
1134     }
1135
1136     // Join slices
1137     typename ImageType::Pointer output;
1138     output = clitk::JoinSlices<ImageType>(o, input, d);
1139     return output;
1140   }
1141   //--------------------------------------------------------------------
1142
1143
1144   //--------------------------------------------------------------------
1145   template<class ImageType>
1146   typename ImageType::Pointer
1147   Clone(const ImageType * input) {
1148     typedef itk::ImageDuplicator<ImageType> DuplicatorType;
1149     typename DuplicatorType::Pointer duplicator = DuplicatorType::New();
1150     duplicator->SetInputImage(input);
1151     duplicator->Update();
1152     return duplicator->GetOutput();
1153   }
1154   //--------------------------------------------------------------------
1155
1156
1157   //--------------------------------------------------------------------
1158   /* Consider an input object, start at A, for each slice (dim1):
1159      - compute the intersection between the AB line and the current slice
1160      - remove what is at lower or greater according to dim2 of this point
1161      - stop at B
1162   */
1163   template<class ImageType>
1164   typename ImageType::Pointer
1165   SliceBySliceSetBackgroundFromSingleLine(const ImageType * input,
1166                                           typename ImageType::PixelType BG,
1167                                           typename ImageType::PointType & A,
1168                                           typename ImageType::PointType & B,
1169                                           int dim1, int dim2, bool removeLowerPartFlag)
1170
1171   {
1172     // Extract slices
1173     typedef typename itk::Image<typename ImageType::PixelType, ImageType::ImageDimension-1> SliceType;
1174     typedef typename SliceType::Pointer SlicePointer;
1175     std::vector<SlicePointer> slices;
1176     clitk::ExtractSlices<ImageType>(input, dim1, slices);
1177
1178     // Start at slice that contains A, and stop at B
1179     typename ImageType::IndexType Ap;
1180     typename ImageType::IndexType Bp;
1181     input->TransformPhysicalPointToIndex(A, Ap);
1182     input->TransformPhysicalPointToIndex(B, Bp);
1183
1184     // Determine slice largest region
1185     typename SliceType::RegionType region = slices[0]->GetLargestPossibleRegion();
1186     typename SliceType::SizeType size = region.GetSize();
1187     typename SliceType::IndexType index = region.GetIndex();
1188
1189     // Line slope
1190     double a = (Bp[dim2]-Ap[dim2])/(Bp[dim1]-Ap[dim1]);
1191     double b = Ap[dim2];
1192
1193     // Loop from slice A to slice B
1194     for(uint i=0; i<(Bp[dim1]-Ap[dim1]); i++) {
1195       // Compute intersection between line AB and current slice for the dim2
1196       double p = a*i+b;
1197       // Change region (lower than dim2)
1198       if (removeLowerPartFlag) {
1199         size[dim2] = p-Ap[dim2];
1200       }
1201       else {
1202         size[dim2] = slices[0]->GetLargestPossibleRegion().GetSize()[dim2]-p;
1203         index[dim2] = p;
1204       }
1205       region.SetSize(size);
1206       region.SetIndex(index);
1207       // Fill region with BG (simple region iterator)
1208       FillRegionWithValue<SliceType>(slices[i+Ap[dim1]], BG, region);
1209       /*
1210       typedef itk::ImageRegionIterator<SliceType> IteratorType;
1211       IteratorType iter(slices[i+Ap[dim1]], region);
1212       iter.GoToBegin();
1213       while (!iter.IsAtEnd()) {
1214         iter.Set(BG);
1215         ++iter;
1216       }
1217       */
1218       // Loop
1219     }
1220
1221     // Merge slices
1222     typename ImageType::Pointer output;
1223     output = clitk::JoinSlices<ImageType>(slices, input, dim1);
1224     return output;
1225   }
1226   //--------------------------------------------------------------------
1227
1228   //--------------------------------------------------------------------
1229   /* Consider an input object, slice by slice, use the point A and set
1230      pixel to BG according to their position relatively to A
1231   */
1232   template<class ImageType>
1233   typename ImageType::Pointer
1234   SliceBySliceSetBackgroundFromPoints(const ImageType * input,
1235                                       typename ImageType::PixelType BG,
1236                                       int sliceDim,
1237                                       std::vector<typename ImageType::PointType> & A,
1238                                       bool removeGreaterThanXFlag,
1239                                       bool removeGreaterThanYFlag)
1240
1241   {
1242     // Extract slices
1243     typedef typename itk::Image<typename ImageType::PixelType, ImageType::ImageDimension-1> SliceType;
1244     typedef typename SliceType::Pointer SlicePointer;
1245     std::vector<SlicePointer> slices;
1246     clitk::ExtractSlices<ImageType>(input, sliceDim, slices);
1247
1248     // Start at slice that contains A
1249     typename ImageType::IndexType Ap;
1250
1251     // Determine slice largest region
1252     typename SliceType::RegionType region = slices[0]->GetLargestPossibleRegion();
1253     typename SliceType::SizeType size = region.GetSize();
1254     typename SliceType::IndexType index = region.GetIndex();
1255
1256     // Loop from slice A to slice B
1257     for(uint i=0; i<A.size(); i++) {
1258       input->TransformPhysicalPointToIndex(A[i], Ap);
1259       uint sliceIndex = Ap[2] - input->GetLargestPossibleRegion().GetIndex()[2];
1260       if ((sliceIndex < 0) || (sliceIndex >= slices.size())) {
1261         continue; // do not consider this slice
1262       }
1263
1264       // Compute region for BG
1265       if (removeGreaterThanXFlag) {
1266         index[0] = Ap[0];
1267         size[0] = region.GetSize()[0]-(index[0]-region.GetIndex()[0]);
1268       }
1269       else {
1270         index[0] = region.GetIndex()[0];
1271         size[0] = Ap[0] - index[0];
1272       }
1273
1274       if (removeGreaterThanYFlag) {
1275         index[1] = Ap[1];
1276         size[1] = region.GetSize()[1]-(index[1]-region.GetIndex()[1]);
1277       }
1278       else {
1279         index[1] = region.GetIndex()[1];
1280         size[1] = Ap[1] - index[1];
1281       }
1282
1283       // Set region
1284       region.SetSize(size);
1285       region.SetIndex(index);
1286
1287       // Fill region with BG (simple region iterator)
1288       FillRegionWithValue<SliceType>(slices[sliceIndex], BG, region);
1289       // Loop
1290     }
1291
1292     // Merge slices
1293     typename ImageType::Pointer output;
1294     output = clitk::JoinSlices<ImageType>(slices, input, sliceDim);
1295     return output;
1296   }
1297   //--------------------------------------------------------------------
1298
1299
1300   //--------------------------------------------------------------------
1301   template<class ImageType>
1302   void
1303   FillRegionWithValue(ImageType * input, typename ImageType::PixelType value, typename ImageType::RegionType & region)
1304   {
1305     typedef itk::ImageRegionIterator<ImageType> IteratorType;
1306     IteratorType iter(input, region);
1307     iter.GoToBegin();
1308     while (!iter.IsAtEnd()) {
1309       iter.Set(value);
1310       ++iter;
1311     }
1312   }
1313   //--------------------------------------------------------------------
1314
1315
1316   //--------------------------------------------------------------------
1317   template<class ImageType>
1318   void
1319   GetMinMaxBoundary(ImageType * input, typename ImageType::PointType & min,
1320                     typename ImageType::PointType & max)
1321   {
1322     typedef typename ImageType::PointType PointType;
1323     typedef typename ImageType::IndexType IndexType;
1324     IndexType min_i, max_i;
1325     min_i = input->GetLargestPossibleRegion().GetIndex();
1326     for(uint i=0; i<ImageType::ImageDimension; i++)
1327       max_i[i] = input->GetLargestPossibleRegion().GetSize()[i] + min_i[i];
1328     input->TransformIndexToPhysicalPoint(min_i, min);
1329     input->TransformIndexToPhysicalPoint(max_i, max);
1330   }
1331   //--------------------------------------------------------------------
1332
1333
1334   //--------------------------------------------------------------------
1335   template<class ImageType>
1336   typename itk::Image<float, ImageType::ImageDimension>::Pointer
1337   DistanceMap(const ImageType * input, typename ImageType::PixelType BG)//,
1338   //              typename itk::Image<float, ImageType::ImageDimension>::Pointer dmap)
1339   {
1340     typedef itk::Image<float,ImageType::ImageDimension> FloatImageType;
1341     typedef itk::SignedMaurerDistanceMapImageFilter<ImageType, FloatImageType> DistanceMapFilterType;
1342     typename DistanceMapFilterType::Pointer filter = DistanceMapFilterType::New();
1343     filter->SetInput(input);
1344     filter->SetUseImageSpacing(true);
1345     filter->SquaredDistanceOff();
1346     filter->SetBackgroundValue(BG);
1347     filter->Update();
1348     return filter->GetOutput();
1349   }
1350   //--------------------------------------------------------------------
1351
1352
1353   //--------------------------------------------------------------------
1354   template<class ImageType>
1355   void
1356   SliceBySliceBuildLineSegmentAccordingToMinimalDistanceBetweenStructures(const ImageType * S1,
1357                                                                           const ImageType * S2,
1358                                                                           typename ImageType::PixelType BG,
1359                                                                           int sliceDimension,
1360                                                                           std::vector<typename ImageType::PointType> & A,
1361                                                                           std::vector<typename ImageType::PointType> & B)
1362   {
1363     // Extract slices
1364     typedef typename itk::Image<typename ImageType::PixelType, 2> SliceType;
1365     typedef typename SliceType::Pointer SlicePointer;
1366     std::vector<SlicePointer> slices_s1;
1367     std::vector<SlicePointer> slices_s2;
1368     clitk::ExtractSlices<ImageType>(S1, sliceDimension, slices_s1);
1369     clitk::ExtractSlices<ImageType>(S2, sliceDimension, slices_s2);
1370
1371     assert(slices_s1.size() == slices_s2.size());
1372
1373     // Prepare dmap
1374     typedef itk::Image<float,2> FloatImageType;
1375     typedef itk::SignedMaurerDistanceMapImageFilter<SliceType, FloatImageType> DistanceMapFilterType;
1376     std::vector<typename FloatImageType::Pointer> dmaps1;
1377     std::vector<typename FloatImageType::Pointer> dmaps2;
1378     typename FloatImageType::Pointer dmap;
1379
1380     // loop on slices
1381     for(uint i=0; i<slices_s1.size(); i++) {
1382       // Compute dmap for S1 *TO PUT IN FONCTION*
1383       dmap = clitk::DistanceMap<SliceType>(slices_s1[i], BG);
1384       dmaps1.push_back(dmap);
1385       writeImage<FloatImageType>(dmap, "dmap1.mha");
1386       // Compute dmap for S2
1387       dmap = clitk::DistanceMap<SliceType>(slices_s2[i], BG);
1388       dmaps2.push_back(dmap);
1389       writeImage<FloatImageType>(dmap, "dmap2.mha");
1390
1391       // Look in S2 for the point the closest to S1
1392       typename SliceType::PointType p = ComputeClosestPoint<SliceType>(slices_s1[i], dmaps2[i], BG);
1393       typename ImageType::PointType p3D;
1394       clitk::PointsUtils<ImageType>::Convert2DTo3D(p, S1, i, p3D);
1395       A.push_back(p3D);
1396
1397       // Look in S2 for the point the closest to S1
1398       p = ComputeClosestPoint<SliceType>(slices_s2[i], dmaps1[i], BG);
1399       clitk::PointsUtils<ImageType>::Convert2DTo3D(p, S2, i, p3D);
1400       B.push_back(p3D);
1401
1402     }
1403
1404     // Debug dmap
1405     /*
1406       typedef itk::Image<float,3> FT;
1407       FT::Pointer f = FT::New();
1408       typename FT::Pointer d1 = clitk::JoinSlices<FT>(dmaps1, S1, 2);
1409       typename FT::Pointer d2 = clitk::JoinSlices<FT>(dmaps2, S2, 2);
1410       writeImage<FT>(d1, "d1.mha");
1411       writeImage<FT>(d2, "d2.mha");
1412     */
1413   }
1414   //--------------------------------------------------------------------
1415
1416
1417   //--------------------------------------------------------------------
1418   template<class ImageType>
1419   typename ImageType::PointType
1420   ComputeClosestPoint(const ImageType * input,
1421                       const itk::Image<float, ImageType::ImageDimension> * dmap,
1422                       typename ImageType::PixelType & BG)
1423   {
1424     // Loop dmap + S2, if FG, get min
1425     typedef itk::Image<float,ImageType::ImageDimension> FloatImageType;
1426     typedef itk::ImageRegionConstIteratorWithIndex<ImageType> ImageIteratorType;
1427     typedef itk::ImageRegionConstIterator<FloatImageType> DMapIteratorType;
1428     ImageIteratorType iter1(input, input->GetLargestPossibleRegion());
1429     DMapIteratorType iter2(dmap, dmap->GetLargestPossibleRegion());
1430
1431     iter1.GoToBegin();
1432     iter2.GoToBegin();
1433     double dmin = 100000.0;
1434     typename ImageType::IndexType indexmin;
1435     while (!iter1.IsAtEnd()) {
1436       if (iter1.Get() != BG) {
1437         double d = iter2.Get();
1438         if (d<dmin) {
1439           indexmin = iter1.GetIndex();
1440           dmin = d;
1441         }
1442       }
1443       ++iter1;
1444       ++iter2;
1445     }
1446
1447     // Convert in Point
1448     typename ImageType::PointType p;
1449     input->TransformIndexToPhysicalPoint(indexmin, p);
1450     return p;
1451   }
1452   //--------------------------------------------------------------------
1453
1454
1455
1456
1457 } // end of namespace
1458