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