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