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); // half of the pixel
 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
 387     image->TransformIndexToPhysicalPoint(image->GetLargestPossibleRegion().GetIndex()+
 388                                          image->GetLargestPossibleRegion().GetSize(), p);
 389
 390     return CropImageAlongOneAxis<ImageType>(image, dim, max, p[dim], autoCrop, BG);
 391   }
 392   //--------------------------------------------------------------------
 393
 394
 395   //--------------------------------------------------------------------
 396   template<class ImageType>
 397   typename ImageType::Pointer
 398   CropImageAlongOneAxis(const ImageType * image,
 399                         int dim, double min, double max,
 400                         bool autoCrop, typename ImageType::PixelType BG)
 401   {
 402     // Compute region size
 403     typename ImageType::RegionType region;
 404     typename ImageType::SizeType size = image->GetLargestPossibleRegion().GetSize();
 405
 406     // Starting index
 407     typename ImageType::PointType p = image->GetOrigin(); // not at pixel center !
 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
 412     // Size of the region
 413     // -1 because last point is size -1
 414     double m = image->GetOrigin()[dim] + (size[dim]-1)*image->GetSpacing()[dim];
 415     if (max > m) p[dim] = m; // Check if not outside the image
 416     else p[dim] = max;
 417
 418     typename ImageType::IndexType end;
 419     image->TransformPhysicalPointToIndex(p, end);
 420     size[dim] = abs(end[dim]-start[dim])+1;// +1 because we want to include the point.
 421
 422     // Set region
 423     region.SetIndex(start);
 424     region.SetSize(size);
 425
 426     // Perform Crop
 427     typedef itk::RegionOfInterestImageFilter<ImageType, ImageType> CropFilterType;
 428     typename CropFilterType::Pointer cropFilter = CropFilterType::New();
 429     cropFilter->SetInput(image);
 430     cropFilter->SetRegionOfInterest(region);
 431     cropFilter->Update();
 432     typename ImageType::Pointer result = cropFilter->GetOutput();
 433
 434     // Auto Crop
 435     if (autoCrop) {
 436       result = AutoCrop<ImageType>(result, BG);
 437     }
 438     return result;
 439   }
 440   //--------------------------------------------------------------------
 441
 442
 443   //--------------------------------------------------------------------
 444   template<class ImageType>
 445   void
 446   ComputeCentroids(const ImageType * image,
 447                    typename ImageType::PixelType BG,
 448                    std::vector<typename ImageType::PointType> & centroids)
 449   {
 450     typedef long LabelType;
 451     static const unsigned int Dim = ImageType::ImageDimension;
 452     typedef itk::ShapeLabelObject< LabelType, Dim > LabelObjectType;
 453     typedef itk::LabelMap< LabelObjectType > LabelMapType;
 454     typedef itk::LabelImageToLabelMapFilter<ImageType, LabelMapType> ImageToMapFilterType;
 455     typename ImageToMapFilterType::Pointer imageToLabelFilter = ImageToMapFilterType::New();
 456     typedef itk::ShapeLabelMapFilter<LabelMapType, ImageType> ShapeFilterType;
 457     typename ShapeFilterType::Pointer statFilter = ShapeFilterType::New();
 458     imageToLabelFilter->SetBackgroundValue(BG);
 459     imageToLabelFilter->SetInput(image);
 460     statFilter->SetInput(imageToLabelFilter->GetOutput());
 461     statFilter->Update();
 462     typename LabelMapType::Pointer labelMap = statFilter->GetOutput();
 463
 464     centroids.clear();
 465     typename ImageType::PointType dummy;
 466     centroids.push_back(dummy); // label 0 -> no centroid, use dummy point for BG
 467     //DS FIXME (not useful ! to change ..)
 468     for(uint i=0; i<labelMap->GetNumberOfLabelObjects(); i++) {
 469       int label = labelMap->GetLabels()[i];
 470       centroids.push_back(labelMap->GetLabelObject(label)->GetCentroid());
 471     }
 472   }
 473   //--------------------------------------------------------------------
 474
 475
 476   //--------------------------------------------------------------------
 477   template<class ImageType, class LabelType>
 478   typename itk::LabelMap< itk::ShapeLabelObject<LabelType, ImageType::ImageDimension> >::Pointer
 479   ComputeLabelMap(const ImageType * image,
 480                   typename ImageType::PixelType BG,
 481                   bool computePerimeterFlag)
 482   {
 483     static const unsigned int Dim = ImageType::ImageDimension;
 484     typedef itk::ShapeLabelObject< LabelType, Dim > LabelObjectType;
 485     typedef itk::LabelMap< LabelObjectType > LabelMapType;
 486     typedef itk::LabelImageToLabelMapFilter<ImageType, LabelMapType> ImageToMapFilterType;
 487     typename ImageToMapFilterType::Pointer imageToLabelFilter = ImageToMapFilterType::New();
 488     typedef itk::ShapeLabelMapFilter<LabelMapType, ImageType> ShapeFilterType;
 489     typename ShapeFilterType::Pointer statFilter = ShapeFilterType::New();
 490     imageToLabelFilter->SetBackgroundValue(BG);
 491     imageToLabelFilter->SetInput(image);
 492     statFilter->SetInput(imageToLabelFilter->GetOutput());
 493     statFilter->SetComputePerimeter(computePerimeterFlag);
 494     statFilter->Update();
 495     return statFilter->GetOutput();
 496   }
 497   //--------------------------------------------------------------------
 498
 499
 500   //--------------------------------------------------------------------
 501   template<class ImageType>
 502   void
 503   ComputeCentroids2(const ImageType * image,
 504                    typename ImageType::PixelType BG,
 505                    std::vector<typename ImageType::PointType> & centroids)
 506   {
 507     typedef long LabelType;
 508     static const unsigned int Dim = ImageType::ImageDimension;
 509     typedef itk::ShapeLabelObject< LabelType, Dim > LabelObjectType;
 510     typedef itk::LabelMap< LabelObjectType > LabelMapType;
 511     typedef itk::LabelImageToLabelMapFilter<ImageType, LabelMapType> ImageToMapFilterType;
 512     typename ImageToMapFilterType::Pointer imageToLabelFilter = ImageToMapFilterType::New();
 513     typedef itk::ShapeLabelMapFilter<LabelMapType, ImageType> ShapeFilterType;
 514     typename ShapeFilterType::Pointer statFilter = ShapeFilterType::New();
 515     imageToLabelFilter->SetBackgroundValue(BG);
 516     imageToLabelFilter->SetInput(image);
 517     statFilter->SetInput(imageToLabelFilter->GetOutput());
 518     statFilter->Update();
 519     typename LabelMapType::Pointer labelMap = statFilter->GetOutput();
 520
 521     centroids.clear();
 522     typename ImageType::PointType dummy;
 523     centroids.push_back(dummy); // label 0 -> no centroid, use dummy point
 524     for(uint i=1; i<labelMap->GetNumberOfLabelObjects()+1; i++) {
 525       centroids.push_back(labelMap->GetLabelObject(i)->GetCentroid());
 526     }
 527
 528     for(uint i=1; i<labelMap->GetNumberOfLabelObjects()+1; i++) {
 529       DD(labelMap->GetLabelObject(i)->GetBinaryPrincipalAxes());
 530       DD(labelMap->GetLabelObject(i)->GetBinaryFlatness());
 531       DD(labelMap->GetLabelObject(i)->GetRoundness ());
 532
 533       // search for the point on the boundary alog PA
 534
 535     }
 536
 537   }
 538   //--------------------------------------------------------------------
 539
 540
 541   //--------------------------------------------------------------------
 542   template<class ImageType>
 543   void
 544   PointsUtils<ImageType>::Convert2DTo3D(const PointType2D & p2D,
 545                                         const ImageType * image,
 546                                         const int slice,
 547                                         PointType3D & p3D)
 548   {
 549     IndexType3D index3D;
 550     index3D[0] = index3D[1] = 0;
 551     index3D[2] = image->GetLargestPossibleRegion().GetIndex()[2]+slice;
 552     image->TransformIndexToPhysicalPoint(index3D, p3D);
 553     p3D[0] = p2D[0];
 554     p3D[1] = p2D[1];
 555     //  p3D[2] = p[2];//(image->GetLargestPossibleRegion().GetIndex()[2]+slice)*image->GetSpacing()[2]
 556     //    + image->GetOrigin()[2];
 557   }
 558   //--------------------------------------------------------------------
 559
 560
 561   //--------------------------------------------------------------------
 562   template<class ImageType>
 563   void
 564   PointsUtils<ImageType>::Convert2DMapTo3DList(const MapPoint2DType & map,
 565                                             const ImageType * image,
 566                                             VectorPoint3DType & list)
 567   {
 568     typename MapPoint2DType::const_iterator iter = map.begin();
 569     while (iter != map.end()) {
 570       PointType3D p;
 571       Convert2DTo3D(iter->second, image, iter->first, p);
 572       list.push_back(p);
 573       ++iter;
 574     }
 575   }
 576   //--------------------------------------------------------------------
 577
 578
 579   //--------------------------------------------------------------------
 580   template<class ImageType>
 581   void
 582   PointsUtils<ImageType>::Convert2DListTo3DList(const VectorPoint2DType & p2D,
 583                                                 int slice,
 584                                                 const ImageType * image,
 585                                                 VectorPoint3DType & list)
 586   {
 587     for(uint i=0; i<p2D.size(); i++) {
 588       PointType3D p;
 589       Convert2DTo3D(p2D[i], image, slice, p);
 590       list.push_back(p);
 591     }
 592   }
 593   //--------------------------------------------------------------------
 594
 595
 596   //--------------------------------------------------------------------
 597   template<class ImageType>
 598   void
 599   WriteListOfLandmarks(std::vector<typename ImageType::PointType> points,
 600                        std::string filename)
 601   {
 602     std::ofstream os;
 603     openFileForWriting(os, filename);
 604     os << "LANDMARKS1" << std::endl;
 605     for(uint i=0; i<points.size(); i++) {
 606       const typename ImageType::PointType & p = points[i];
 607       // Write it in the file
 608       os << i << " " << p[0] << " " << p[1] << " " << p[2] << " 0 0 " << std::endl;
 609     }
 610     os.close();
 611   }
 612   //--------------------------------------------------------------------
 613
 614
 615   //--------------------------------------------------------------------
 616   template<class ImageType>
 617   typename ImageType::Pointer
 618   Dilate(const ImageType * image, double radiusInMM,
 619          typename ImageType::PixelType BG,
 620          typename ImageType::PixelType FG,
 621          bool extendSupport)
 622   {
 623     typename ImageType::SizeType r;
 624     for(uint i=0; i<ImageType::ImageDimension; i++)
 625       r[i] = (uint)lrint(radiusInMM/image->GetSpacing()[i]);
 626     return Dilate<ImageType>(image, r, BG, FG, extendSupport);
 627   }
 628   //--------------------------------------------------------------------
 629
 630
 631   //--------------------------------------------------------------------
 632   template<class ImageType>
 633   typename ImageType::Pointer
 634   Dilate(const ImageType * image, typename ImageType::PointType radiusInMM,
 635          typename ImageType::PixelType BG,
 636          typename ImageType::PixelType FG,
 637          bool extendSupport)
 638   {
 639     typename ImageType::SizeType r;
 640     for(uint i=0; i<ImageType::ImageDimension; i++)
 641       r[i] = (uint)lrint(radiusInMM[i]/image->GetSpacing()[i]);
 642     return Dilate<ImageType>(image, r, BG, FG, extendSupport);
 643   }
 644   //--------------------------------------------------------------------
 645
 646
 647   //--------------------------------------------------------------------
 648   template<class ImageType>
 649   typename ImageType::Pointer
 650   Dilate(const ImageType * image, typename ImageType::SizeType radius,
 651          typename ImageType::PixelType BG,
 652          typename ImageType::PixelType FG,
 653          bool extendSupport)
 654   {
 655     // Create kernel for dilatation
 656     typedef itk::BinaryBallStructuringElement<typename ImageType::PixelType,
 657                                               ImageType::ImageDimension> KernelType;
 658     KernelType structuringElement;
 659     structuringElement.SetRadius(radius);
 660     structuringElement.CreateStructuringElement();
 661
 662     typename ImageType::Pointer output;
 663     if (extendSupport) {
 664       typedef itk::ConstantPadImageFilter<ImageType, ImageType> PadFilterType;
 665       typename PadFilterType::Pointer padFilter = PadFilterType::New();
 666       padFilter->SetInput(image);
 667       typename ImageType::SizeType lower;
 668       typename ImageType::SizeType upper;
 669       for(uint i=0; i<3; i++) {
 670         lower[i] = upper[i] = 2*(radius[i]+1);
 671       }
 672       padFilter->SetPadLowerBound(lower);
 673       padFilter->SetPadUpperBound(upper);
 674       padFilter->Update();
 675       output = padFilter->GetOutput();
 676     }
 677
 678     // Dilate  filter
 679     typedef itk::BinaryDilateImageFilter<ImageType, ImageType , KernelType> DilateFilterType;
 680     typename DilateFilterType::Pointer dilateFilter = DilateFilterType::New();
 681     dilateFilter->SetBackgroundValue(BG);
 682     dilateFilter->SetForegroundValue(FG);
 683     dilateFilter->SetBoundaryToForeground(false);
 684     dilateFilter->SetKernel(structuringElement);
 685     if (extendSupport) dilateFilter->SetInput(output);
 686     else dilateFilter->SetInput(image);
 687     dilateFilter->Update();
 688     return dilateFilter->GetOutput();
 689   }
 690   //--------------------------------------------------------------------
 691
 692
 693   //--------------------------------------------------------------------
 694   template<class ImageType>
 695   typename ImageType::Pointer
 696   Opening(const ImageType * image, typename ImageType::SizeType radius,
 697          typename ImageType::PixelType BG,
 698          typename ImageType::PixelType FG)
 699   {
 700     // Kernel
 701     typedef itk::BinaryBallStructuringElement<typename ImageType::PixelType,
 702                                               ImageType::ImageDimension> KernelType;
 703     KernelType structuringElement;
 704     structuringElement.SetRadius(radius);
 705     structuringElement.CreateStructuringElement();
 706
 707     // Filter
 708     typedef itk::BinaryMorphologicalOpeningImageFilter<ImageType, ImageType , KernelType> OpeningFilterType;
 709     typename OpeningFilterType::Pointer open = OpeningFilterType::New();
 710     open->SetInput(image);
 711     open->SetBackgroundValue(BG);
 712     open->SetForegroundValue(FG);
 713     open->SetKernel(structuringElement);
 714     open->Update();
 715     return open->GetOutput();
 716   }
 717   //--------------------------------------------------------------------
 718
 719
 720
 721   //--------------------------------------------------------------------
 722   template<class ValueType, class VectorType>
 723   void ConvertOption(std::string optionName, uint given,
 724                      ValueType * values, VectorType & p,
 725                      uint dim, bool required)
 726   {
 727     if (required && (given == 0)) {
 728       clitkExceptionMacro("The option --" << optionName << " must be set and have 1 or "
 729                           << dim << " values.");
 730     }
 731     if (given == 1) {
 732       for(uint i=0; i<dim; i++) p[i] = values[0];
 733       return;
 734     }
 735     if (given == dim) {
 736       for(uint i=0; i<dim; i++) p[i] = values[i];
 737       return;
 738     }
 739     if (given == 0) return;
 740     clitkExceptionMacro("The option --" << optionName << " must have 1 or "
 741                         << dim << " values.");
 742   }
 743   //--------------------------------------------------------------------
 744
 745
 746   //--------------------------------------------------------------------
 747   /*
 748     http://www.gamedev.net/community/forums/topic.asp?topic_id=542870
 749     Assuming the points are (Ax,Ay) (Bx,By) and (Cx,Cy), you need to compute:
 750     (Bx - Ax) * (Cy - Ay) - (By - Ay) * (Cx - Ax)
 751     This will equal zero if the point C is on the line formed by
 752     points A and B, and will have a different sign depending on the
 753     side. Which side this is depends on the orientation of your (x,y)
 754     coordinates, but you can plug test values for A,B and C into this
 755     formula to determine whether negative values are to the left or to
 756     the right.
 757     => to accelerate, start with formula, when change sign -> stop and fill
 758
 759     offsetToKeep = is used to determine which side of the line we
 760     keep. The point along the mainDirection but 'offsetToKeep' mm away
 761     is kept.
 762
 763   */
 764   template<class ImageType>
 765   void
 766   SliceBySliceSetBackgroundFromLineSeparation(ImageType * input,
 767                                               std::vector<typename ImageType::PointType> & lA,
 768                                               std::vector<typename ImageType::PointType> & lB,
 769                                               typename ImageType::PixelType BG,
 770                                               int mainDirection,
 771                                               double offsetToKeep,
 772                                               bool keepIfEqual)
 773   {
 774     assert((mainDirection==0) || (mainDirection==1));
 775     typename ImageType::PointType offset;
 776     offset[0] = offset[1] = offset[2] = 0.0;
 777     offset[mainDirection] = offsetToKeep;
 778     SliceBySliceSetBackgroundFromLineSeparation<ImageType>(input, lA, lB, BG, offset, keepIfEqual);
 779   }
 780   template<class ImageType>
 781   void
 782   SliceBySliceSetBackgroundFromLineSeparation(ImageType * input,
 783                                               std::vector<typename ImageType::PointType> & lA,
 784                                               std::vector<typename ImageType::PointType> & lB,
 785                                               typename ImageType::PixelType BG,
 786                                               typename ImageType::PointType offsetToKeep,
 787                                               bool keepIfEqual)
 788   {
 789     typedef itk::ImageSliceIteratorWithIndex<ImageType> SliceIteratorType;
 790     SliceIteratorType siter = SliceIteratorType(input, input->GetLargestPossibleRegion());
 791     siter.SetFirstDirection(0);
 792     siter.SetSecondDirection(1);
 793     siter.GoToBegin();
 794     uint i=0;
 795     typename ImageType::PointType A;
 796     typename ImageType::PointType B;
 797     typename ImageType::PointType C;
 798     assert(lA.size() == lB.size());
 799     while ((i<lA.size()) && (!siter.IsAtEnd())) {
 800       // Check that the current slice correspond to the current point
 801       input->TransformIndexToPhysicalPoint(siter.GetIndex(), C);
 802       if ((fabs(C[2] - lA[i][2]))>0.01) { // is !equal with a tolerance of 0.01 mm
 803       }
 804       else {
 805         // Define A,B,C points
 806         A = lA[i];
 807         B = lB[i];
 808         C = A;
 809         // Check that the line is not a point (A=B)
 810         bool p = (A[0] == B[0]) && (A[1] == B[1]);
 811
 812         if (!p) {
 813           //C[mainDirection] += offsetToKeep; // I know I must keep this point
 814           C[0] += offsetToKeep[0];
 815           C[1] += offsetToKeep[1];
 816           //C[2] += offsetToKeep[2];
 817           double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
 818           bool isPositive = s<0;
 819           while (!siter.IsAtEndOfSlice()) {
 820             while (!siter.IsAtEndOfLine()) {
 821               // Very slow, I know ... but image should be very small
 822               input->TransformIndexToPhysicalPoint(siter.GetIndex(), C);
 823               double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
 824               if (s == 0) {
 825                 if (!keepIfEqual) siter.Set(BG); // on the line, we decide to remove
 826               }
 827               if (isPositive) {
 828                 if (s > 0) siter.Set(BG);
 829               }
 830               else {
 831                 if (s < 0) siter.Set(BG);
 832               }
 833               ++siter;
 834             }
 835             siter.NextLine();
 836           } // end loop slice
 837         }
 838
 839         ++i;
 840       } // End of current slice
 841       siter.NextSlice();
 842     }
 843   }
 844   //--------------------------------------------------------------------
 845
 846
 847   //--------------------------------------------------------------------
 848   template<class ImageType>
 849   void
 850   AndNot(ImageType * input,
 851          const ImageType * object,
 852          typename ImageType::PixelType BG)
 853   {
 854     typename ImageType::Pointer o;
 855     bool resized=false;
 856     if (!clitk::HaveSameSizeAndSpacing<ImageType, ImageType>(input, object)) {
 857       o = clitk::ResizeImageLike<ImageType>(object, input, BG);
 858       resized = true;
 859     }
 860
 861     typedef clitk::BooleanOperatorLabelImageFilter<ImageType> BoolFilterType;
 862     typename BoolFilterType::Pointer boolFilter = BoolFilterType::New();
 863     boolFilter->InPlaceOn();
 864     boolFilter->SetInput1(input);
 865     if (resized) boolFilter->SetInput2(o);
 866     else boolFilter->SetInput2(object);
 867     boolFilter->SetBackgroundValue1(BG);
 868     boolFilter->SetBackgroundValue2(BG);
 869     boolFilter->SetOperationType(BoolFilterType::AndNot);
 870     boolFilter->Update();
 871   }
 872   //--------------------------------------------------------------------
 873
 874
 875   //--------------------------------------------------------------------
 876   template<class ImageType>
 877   void
 878   And(ImageType * input,
 879       const ImageType * object,
 880       typename ImageType::PixelType BG)
 881   {
 882     typename ImageType::Pointer o;
 883     bool resized=false;
 884     if (!clitk::HaveSameSizeAndSpacing<ImageType, ImageType>(input, object)) {
 885       o = clitk::ResizeImageLike<ImageType>(object, input, BG);
 886       resized = true;
 887     }
 888
 889     typedef clitk::BooleanOperatorLabelImageFilter<ImageType> BoolFilterType;
 890     typename BoolFilterType::Pointer boolFilter = BoolFilterType::New();
 891     boolFilter->InPlaceOn();
 892     boolFilter->SetInput1(input);
 893     if (resized) boolFilter->SetInput2(o);
 894     else boolFilter->SetInput2(object);
 895     boolFilter->SetBackgroundValue1(BG);
 896     boolFilter->SetBackgroundValue2(BG);
 897     boolFilter->SetOperationType(BoolFilterType::And);
 898     boolFilter->Update();
 899   }
 900   //--------------------------------------------------------------------
 901
 902
 903   //--------------------------------------------------------------------
 904   template<class ImageType>
 905   void
 906   Or(ImageType * input,
 907      const ImageType * object,
 908      typename ImageType::PixelType BG)
 909   {
 910     typename ImageType::Pointer o;
 911     bool resized=false;
 912     if (!clitk::HaveSameSizeAndSpacing<ImageType, ImageType>(input, object)) {
 913       o = clitk::ResizeImageLike<ImageType>(object, input, BG);
 914       resized = true;
 915     }
 916
 917     typedef clitk::BooleanOperatorLabelImageFilter<ImageType> BoolFilterType;
 918     typename BoolFilterType::Pointer boolFilter = BoolFilterType::New();
 919     boolFilter->InPlaceOn();
 920     boolFilter->SetInput1(input);
 921     if (resized) boolFilter->SetInput2(o);
 922     else boolFilter->SetInput2(object);
 923     boolFilter->SetBackgroundValue1(BG);
 924     boolFilter->SetBackgroundValue2(BG);
 925     boolFilter->SetOperationType(BoolFilterType::Or);
 926     boolFilter->Update();
 927   }
 928   //--------------------------------------------------------------------
 929
 930
 931   //--------------------------------------------------------------------
 932   template<class ImageType>
 933   typename ImageType::Pointer
 934   Binarize(const ImageType * input,
 935            typename ImageType::PixelType lower,
 936            typename ImageType::PixelType upper,
 937            typename ImageType::PixelType BG,
 938            typename ImageType::PixelType FG)
 939   {
 940     typedef itk::BinaryThresholdImageFilter<ImageType, ImageType> BinaryThresholdFilterType;
 941     typename BinaryThresholdFilterType::Pointer binarizeFilter = BinaryThresholdFilterType::New();
 942     binarizeFilter->SetInput(input);
 943     binarizeFilter->InPlaceOff();
 944     binarizeFilter->SetLowerThreshold(lower);
 945     binarizeFilter->SetUpperThreshold(upper);
 946     binarizeFilter->SetInsideValue(FG);
 947     binarizeFilter->SetOutsideValue(BG);
 948     binarizeFilter->Update();
 949     return binarizeFilter->GetOutput();
 950   }
 951   //--------------------------------------------------------------------
 952
 953
 954   //--------------------------------------------------------------------
 955   template<class ImageType>
 956   void
 957   GetMinMaxPointPosition(const ImageType * input,
 958                          typename ImageType::PointType & min,
 959                          typename ImageType::PointType & max)
 960   {
 961     typename ImageType::IndexType index = input->GetLargestPossibleRegion().GetIndex();
 962     input->TransformIndexToPhysicalPoint(index, min);
 963     index = index+input->GetLargestPossibleRegion().GetSize();
 964     input->TransformIndexToPhysicalPoint(index, max);
 965   }
 966   //--------------------------------------------------------------------
 967
 968
 969   //--------------------------------------------------------------------
 970   template<class ImageType>
 971   typename ImageType::PointType
 972   FindExtremaPointInAGivenLine(const ImageType * input,
 973                                int dimension,
 974                                bool inverse,
 975                                typename ImageType::PointType p,
 976                                typename ImageType::PixelType BG,
 977                                double distanceMax)
 978   {
 979     // Which direction ?  Increasing or decreasing.
 980     int d=1;
 981     if (inverse) d=-1;
 982
 983     // Transform to pixel index
 984     typename ImageType::IndexType index;
 985     input->TransformPhysicalPointToIndex(p, index);
 986
 987     // Loop while inside the mask;
 988     while (input->GetPixel(index) != BG) {
 989       index[dimension] += d;
 990     }
 991
 992     // Transform back to Physical Units
 993     typename ImageType::PointType result;
 994     input->TransformIndexToPhysicalPoint(index, result);
 995
 996     // Check that is is not too far away
 997     double distance = p.EuclideanDistanceTo(result);
 998     if (distance > distanceMax) {
 999       result = p; // Get back to initial value
1000     }
1001
1002     return result;
1003   }
1004   //--------------------------------------------------------------------
1005
1006
1007   //--------------------------------------------------------------------
1008   template<class PointType>
1009   bool
1010   IsOnTheSameLineSide(PointType C, PointType A, PointType B, PointType like)
1011   {
1012     // Look at the position of point 'like' according to the AB line
1013     double s = (B[0] - A[0]) * (like[1] - A[1]) - (B[1] - A[1]) * (like[0] - A[0]);
1014     bool negative = s<0;
1015
1016     // Look the C position
1017     s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
1018
1019     if (negative && (s<=0)) return true;
1020     if (!negative && (s>=0)) return true;
1021     return false;
1022   }
1023   //--------------------------------------------------------------------
1024
1025
1026   //--------------------------------------------------------------------
1027   /* Consider an input object, for each slice, find the extrema
1028      position according to a given direction and build a line segment
1029      passing throught this point in a given direction.  Output is a
1030      vector of line (from point A to B), for each slice;
1031    */
1032   template<class ImageType>
1033   void
1034   SliceBySliceBuildLineSegmentAccordingToExtremaPosition(const ImageType * input,
1035                                                          typename ImageType::PixelType BG,
1036                                                          int sliceDimension,
1037                                                          int extremaDirection,
1038                                                          bool extremaOppositeFlag,
1039                                                          int lineDirection,
1040                                                          double margin,
1041                                                          std::vector<typename ImageType::PointType> & A,
1042                                                          std::vector<typename ImageType::PointType> & B)
1043   {
1044     // Type of a slice
1045     typedef typename itk::Image<typename ImageType::PixelType, ImageType::ImageDimension-1> SliceType;
1046
1047     // Build the list of slices
1048     std::vector<typename SliceType::Pointer> slices;
1049     clitk::ExtractSlices<ImageType>(input, sliceDimension, slices);
1050
1051     // Build the list of 2D points
1052     std::map<int, typename SliceType::PointType> position2D;
1053     for(uint i=0; i<slices.size(); i++) {
1054       typename SliceType::PointType p;
1055       bool found =
1056         clitk::FindExtremaPointInAGivenDirection<SliceType>(slices[i], BG,
1057                                                             extremaDirection, extremaOppositeFlag, p);
1058       if (found) {
1059         position2D[i] = p;
1060       }
1061     }
1062
1063     // Convert 2D points in slice into 3D points
1064     clitk::PointsUtils<ImageType>::Convert2DMapTo3DList(position2D, input, A);
1065
1066     // Create additional point just right to the previous ones, on the
1067     // given lineDirection, in order to create a horizontal/vertical line.
1068     for(uint i=0; i<A.size(); i++) {
1069       typename ImageType::PointType p = A[i];
1070       p[lineDirection] += 10;
1071       B.push_back(p);
1072       // Margins ?
1073       A[i][extremaDirection] += margin;
1074       B[i][extremaDirection] += margin;
1075     }
1076
1077   }
1078   //--------------------------------------------------------------------
1079
1080
1081   //--------------------------------------------------------------------
1082   template<class ImageType>
1083   typename ImageType::Pointer
1084   SliceBySliceKeepMainCCL(const ImageType * input,
1085                           typename ImageType::PixelType BG,
1086                           typename ImageType::PixelType FG)  {
1087
1088     // Extract slices
1089     const int d = ImageType::ImageDimension-1;
1090     typedef typename itk::Image<typename ImageType::PixelType, d> SliceType;
1091     std::vector<typename SliceType::Pointer> slices;
1092     clitk::ExtractSlices<ImageType>(input, d, slices);
1093
1094     // Labelize and keep the main one
1095     std::vector<typename SliceType::Pointer> o;
1096     for(uint i=0; i<slices.size(); i++) {
1097       o.push_back(clitk::Labelize<SliceType>(slices[i], BG, false, 1));
1098       o[i] = clitk::KeepLabels<SliceType>(o[i], BG, FG, 1, 1, true);
1099     }
1100
1101     // Join slices
1102     typename ImageType::Pointer output;
1103     output = clitk::JoinSlices<ImageType>(o, input, d);
1104     return output;
1105   }
1106   //--------------------------------------------------------------------
1107
1108
1109   //--------------------------------------------------------------------
1110   template<class ImageType>
1111   typename ImageType::Pointer
1112   Clone(const ImageType * input) {
1113     typedef itk::ImageDuplicator<ImageType> DuplicatorType;
1114     typename DuplicatorType::Pointer duplicator = DuplicatorType::New();
1115     duplicator->SetInputImage(input);
1116     duplicator->Update();
1117     return duplicator->GetOutput();
1118   }
1119   //--------------------------------------------------------------------
1120
1121
1122   //--------------------------------------------------------------------
1123   /* Consider an input object, start at A, for each slice (dim1):
1124      - compute the intersection between the AB line and the current slice
1125      - remove what is at lower or greater according to dim2 of this point
1126      - stop at B
1127   */
1128   template<class ImageType>
1129   typename ImageType::Pointer
1130   SliceBySliceSetBackgroundFromSingleLine(const ImageType * input,
1131                                           typename ImageType::PixelType BG,
1132                                           typename ImageType::PointType & A,
1133                                           typename ImageType::PointType & B,
1134                                           int dim1, int dim2, bool removeLowerPartFlag)
1135
1136   {
1137     // Extract slices
1138     typedef typename itk::Image<typename ImageType::PixelType, ImageType::ImageDimension-1> SliceType;
1139     typedef typename SliceType::Pointer SlicePointer;
1140     std::vector<SlicePointer> slices;
1141     clitk::ExtractSlices<ImageType>(input, dim1, slices);
1142
1143     // Start at slice that contains A, and stop at B
1144     typename ImageType::IndexType Ap;
1145     typename ImageType::IndexType Bp;
1146     input->TransformPhysicalPointToIndex(A, Ap);
1147     input->TransformPhysicalPointToIndex(B, Bp);
1148
1149     // Determine slice largest region
1150     typename SliceType::RegionType region = slices[0]->GetLargestPossibleRegion();
1151     typename SliceType::SizeType size = region.GetSize();
1152     typename SliceType::IndexType index = region.GetIndex();
1153
1154     // Line slope
1155     double a = (Bp[dim2]-Ap[dim2])/(Bp[dim1]-Ap[dim1]);
1156     double b = Ap[dim2];
1157
1158     // Loop from slice A to slice B
1159     for(uint i=0; i<(Bp[dim1]-Ap[dim1]); i++) {
1160       // Compute intersection between line AB and current slice for the dim2
1161       double p = a*i+b;
1162       // Change region (lower than dim2)
1163       if (removeLowerPartFlag) {
1164         size[dim2] = p-Ap[dim2];
1165       }
1166       else {
1167         size[dim2] = slices[0]->GetLargestPossibleRegion().GetSize()[dim2]-p;
1168         index[dim2] = p;
1169       }
1170       region.SetSize(size);
1171       region.SetIndex(index);
1172       // Fill region with BG (simple region iterator)
1173       FillRegionWithValue<SliceType>(slices[i+Ap[dim1]], BG, region);
1174       /*
1175       typedef itk::ImageRegionIterator<SliceType> IteratorType;
1176       IteratorType iter(slices[i+Ap[dim1]], region);
1177       iter.GoToBegin();
1178       while (!iter.IsAtEnd()) {
1179         iter.Set(BG);
1180         ++iter;
1181       }
1182       */
1183       // Loop
1184     }
1185
1186     // Merge slices
1187     typename ImageType::Pointer output;
1188     output = clitk::JoinSlices<ImageType>(slices, input, dim1);
1189     return output;
1190   }
1191   //--------------------------------------------------------------------
1192
1193   //--------------------------------------------------------------------
1194   /* Consider an input object, slice by slice, use the point A and set
1195      pixel to BG according to their position relatively to A
1196   */
1197   template<class ImageType>
1198   typename ImageType::Pointer
1199   SliceBySliceSetBackgroundFromPoints(const ImageType * input,
1200                                       typename ImageType::PixelType BG,
1201                                       int sliceDim,
1202                                       std::vector<typename ImageType::PointType> & A,
1203                                       bool removeGreaterThanXFlag,
1204                                       bool removeGreaterThanYFlag)
1205
1206   {
1207     // Extract slices
1208     typedef typename itk::Image<typename ImageType::PixelType, ImageType::ImageDimension-1> SliceType;
1209     typedef typename SliceType::Pointer SlicePointer;
1210     std::vector<SlicePointer> slices;
1211     clitk::ExtractSlices<ImageType>(input, sliceDim, slices);
1212
1213     // Start at slice that contains A
1214     typename ImageType::IndexType Ap;
1215
1216     // Determine slice largest region
1217     typename SliceType::RegionType region = slices[0]->GetLargestPossibleRegion();
1218     typename SliceType::SizeType size = region.GetSize();
1219     typename SliceType::IndexType index = region.GetIndex();
1220
1221     // Loop from slice A to slice B
1222     for(uint i=0; i<A.size(); i++) {
1223       input->TransformPhysicalPointToIndex(A[i], Ap);
1224       uint sliceIndex = Ap[2] - input->GetLargestPossibleRegion().GetIndex()[2];
1225       if ((sliceIndex < 0) || (sliceIndex >= slices.size())) {
1226         continue; // do not consider this slice
1227       }
1228
1229       // Compute region for BG
1230       if (removeGreaterThanXFlag) {
1231         index[0] = Ap[0];
1232         size[0] = region.GetSize()[0]-(index[0]-region.GetIndex()[0]);
1233       }
1234       else {
1235         index[0] = region.GetIndex()[0];
1236         size[0] = Ap[0] - index[0];
1237       }
1238
1239       if (removeGreaterThanYFlag) {
1240         index[1] = Ap[1];
1241         size[1] = region.GetSize()[1]-(index[1]-region.GetIndex()[1]);
1242       }
1243       else {
1244         index[1] = region.GetIndex()[1];
1245         size[1] = Ap[1] - index[1];
1246       }
1247
1248       // Set region
1249       region.SetSize(size);
1250       region.SetIndex(index);
1251
1252       // Fill region with BG (simple region iterator)
1253       FillRegionWithValue<SliceType>(slices[sliceIndex], BG, region);
1254       // Loop
1255     }
1256
1257     // Merge slices
1258     typename ImageType::Pointer output;
1259     output = clitk::JoinSlices<ImageType>(slices, input, sliceDim);
1260     return output;
1261   }
1262   //--------------------------------------------------------------------
1263
1264
1265   //--------------------------------------------------------------------
1266   template<class ImageType>
1267   void
1268   FillRegionWithValue(ImageType * input, typename ImageType::PixelType value, typename ImageType::RegionType & region)
1269   {
1270     typedef itk::ImageRegionIterator<ImageType> IteratorType;
1271     IteratorType iter(input, region);
1272     iter.GoToBegin();
1273     while (!iter.IsAtEnd()) {
1274       iter.Set(value);
1275       ++iter;
1276     }
1277   }
1278   //--------------------------------------------------------------------
1279
1280
1281   //--------------------------------------------------------------------
1282   template<class ImageType>
1283   void
1284   GetMinMaxBoundary(ImageType * input, typename ImageType::PointType & min,
1285                     typename ImageType::PointType & max)
1286   {
1287     typedef typename ImageType::PointType PointType;
1288     typedef typename ImageType::IndexType IndexType;
1289     IndexType min_i, max_i;
1290     min_i = input->GetLargestPossibleRegion().GetIndex();
1291     for(uint i=0; i<ImageType::ImageDimension; i++)
1292       max_i[i] = input->GetLargestPossibleRegion().GetSize()[i] + min_i[i];
1293     input->TransformIndexToPhysicalPoint(min_i, min);
1294     input->TransformIndexToPhysicalPoint(max_i, max);
1295   }
1296   //--------------------------------------------------------------------
1297
1298
1299   //--------------------------------------------------------------------
1300   template<class ImageType>
1301   typename itk::Image<float, ImageType::ImageDimension>::Pointer
1302   DistanceMap(const ImageType * input, typename ImageType::PixelType BG)//,
1303   //              typename itk::Image<float, ImageType::ImageDimension>::Pointer dmap)
1304   {
1305     typedef itk::Image<float,ImageType::ImageDimension> FloatImageType;
1306     typedef itk::SignedMaurerDistanceMapImageFilter<ImageType, FloatImageType> DistanceMapFilterType;
1307     typename DistanceMapFilterType::Pointer filter = DistanceMapFilterType::New();
1308     filter->SetInput(input);
1309     filter->SetUseImageSpacing(true);
1310     filter->SquaredDistanceOff();
1311     filter->SetBackgroundValue(BG);
1312     filter->Update();
1313     return filter->GetOutput();
1314   }
1315   //--------------------------------------------------------------------
1316
1317
1318   //--------------------------------------------------------------------
1319   template<class ImageType>
1320   void
1321   SliceBySliceBuildLineSegmentAccordingToMinimalDistanceBetweenStructures(const ImageType * S1,
1322                                                                           const ImageType * S2,
1323                                                                           typename ImageType::PixelType BG,
1324                                                                           int sliceDimension,
1325                                                                           std::vector<typename ImageType::PointType> & A,
1326                                                                           std::vector<typename ImageType::PointType> & B)
1327   {
1328     // Extract slices
1329     typedef typename itk::Image<typename ImageType::PixelType, 2> SliceType;
1330     typedef typename SliceType::Pointer SlicePointer;
1331     std::vector<SlicePointer> slices_s1;
1332     std::vector<SlicePointer> slices_s2;
1333     clitk::ExtractSlices<ImageType>(S1, sliceDimension, slices_s1);
1334     clitk::ExtractSlices<ImageType>(S2, sliceDimension, slices_s2);
1335
1336     assert(slices_s1.size() == slices_s2.size());
1337
1338     // Prepare dmap
1339     typedef itk::Image<float,2> FloatImageType;
1340     typedef itk::SignedMaurerDistanceMapImageFilter<SliceType, FloatImageType> DistanceMapFilterType;
1341     std::vector<typename FloatImageType::Pointer> dmaps1;
1342     std::vector<typename FloatImageType::Pointer> dmaps2;
1343     typename FloatImageType::Pointer dmap;
1344
1345     // loop on slices
1346     for(uint i=0; i<slices_s1.size(); i++) {
1347       // Compute dmap for S1 *TO PUT IN FONCTION*
1348       dmap = clitk::DistanceMap<SliceType>(slices_s1[i], BG);
1349       dmaps1.push_back(dmap);
1350       //writeImage<FloatImageType>(dmap, "dmap1.mha");
1351       // Compute dmap for S2
1352       dmap = clitk::DistanceMap<SliceType>(slices_s2[i], BG);
1353       dmaps2.push_back(dmap);
1354       //writeImage<FloatImageType>(dmap, "dmap2.mha");
1355
1356       // Look in S2 for the point the closest to S1
1357       typename SliceType::PointType p = ComputeClosestPoint<SliceType>(slices_s1[i], dmaps2[i], BG);
1358       typename ImageType::PointType p3D;
1359       clitk::PointsUtils<ImageType>::Convert2DTo3D(p, S1, i, p3D);
1360       A.push_back(p3D);
1361
1362       // Look in S2 for the point the closest to S1
1363       p = ComputeClosestPoint<SliceType>(slices_s2[i], dmaps1[i], BG);
1364       clitk::PointsUtils<ImageType>::Convert2DTo3D(p, S2, i, p3D);
1365       B.push_back(p3D);
1366
1367     }
1368
1369     /*
1370     // Debug dmap
1371     typedef itk::Image<float,3> FT;
1372     FT::Pointer f = FT::New();
1373     typename FT::Pointer d1 = clitk::JoinSlices<FT>(dmaps1, S1, 2);
1374     typename FT::Pointer d2 = clitk::JoinSlices<FT>(dmaps2, S2, 2);
1375     writeImage<FT>(d1, "d1.mha");
1376     writeImage<FT>(d2, "d2.mha");
1377     */
1378   }
1379   //--------------------------------------------------------------------
1380
1381
1382   //--------------------------------------------------------------------
1383   template<class ImageType>
1384   typename ImageType::PointType
1385   ComputeClosestPoint(const ImageType * input,
1386                       const itk::Image<float, ImageType::ImageDimension> * dmap,
1387                       typename ImageType::PixelType & BG)
1388   {
1389     // Loop dmap + S2, if FG, get min
1390     typedef itk::Image<float,ImageType::ImageDimension> FloatImageType;
1391     typedef itk::ImageRegionConstIteratorWithIndex<ImageType> ImageIteratorType;
1392     typedef itk::ImageRegionConstIterator<FloatImageType> DMapIteratorType;
1393     ImageIteratorType iter1(input, input->GetLargestPossibleRegion());
1394     DMapIteratorType iter2(dmap, dmap->GetLargestPossibleRegion());
1395
1396     iter1.GoToBegin();
1397     iter2.GoToBegin();
1398     double dmin = 100000.0;
1399     typename ImageType::IndexType indexmin;
1400     indexmin.Fill(0);
1401     while (!iter1.IsAtEnd()) {
1402       if (iter1.Get() != BG) {
1403         double d = iter2.Get();
1404         if (d<dmin) {
1405           indexmin = iter1.GetIndex();
1406           dmin = d;
1407         }
1408       }
1409       ++iter1;
1410       ++iter2;
1411     }
1412
1413     // Convert in Point
1414     typename ImageType::PointType p;
1415     input->TransformIndexToPhysicalPoint(indexmin, p);
1416     return p;
1417   }
1418   //--------------------------------------------------------------------
1419
1420
1421   //--------------------------------------------------------------------
1422   template<class ImageType>
1423   typename ImageType::Pointer
1424   RemoveNegativeIndexFromRegion(ImageType * input) {
1425     typedef itk::ChangeInformationImageFilter< ImageType > InfoFilterType;
1426     typename InfoFilterType::Pointer indexChangeFilter = InfoFilterType::New();
1427     indexChangeFilter->ChangeRegionOn();
1428     // The next line is commented because not exist in itk 3
1429     // typename InfoFilterType::OutputImageOffsetValueType indexShift[3];
1430     long indexShift[3];
1431     typename ImageType::IndexType index = input->GetLargestPossibleRegion().GetIndex();
1432     for(uint i=0;i<ImageType::ImageDimension; i++)
1433       indexShift[i] = (index[i]<0 ? -index[i]:0);
1434     typename ImageType::PointType origin;
1435     for(uint i=0;i<ImageType::ImageDimension; i++)
1436     origin[i] = input->GetOrigin()[i] - indexShift[i]*input->GetSpacing()[i];
1437     indexChangeFilter->SetOutputOffset( indexShift );
1438     indexChangeFilter->SetInput(input);
1439     indexChangeFilter->SetOutputOrigin(origin);
1440     indexChangeFilter->ChangeOriginOn();
1441     indexChangeFilter->Update();
1442     return indexChangeFilter->GetOutput();
1443   }
1444   //--------------------------------------------------------------------
1445
1446
1447 } // end of namespace
1448