itk/clitkSegmentationUtils.txx

   1 /*=========================================================================
   2   Program:   vv                     http://www.creatis.insa-lyon.fr/rio/vv
   3
   4   Authors belong to:
   5   - University of LYON              http://www.universite-lyon.fr/
   6   - Léon Bérard cancer center       http://oncora1.lyon.fnclcc.fr
   7   - CREATIS CNRS laboratory         http://www.creatis.insa-lyon.fr
   8
   9   This software is distributed WITHOUT ANY WARRANTY; without even
  10   the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  11   PURPOSE.  See the copyright notices for more information.
  12
  13   It is distributed under dual licence
  14
  15   - BSD        See included LICENSE.txt file
  16   - CeCILL-B   http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
  17   ======================================================================-====*/
  18
  19 // clitk
  20 #include "clitkSetBackgroundImageFilter.h"
  21 #include "clitkSliceBySliceRelativePositionFilter.h"
  22 #include "clitkCropLikeImageFilter.h"
  23 #include "clitkMemoryUsage.h"
  24
  25 // itk
  26 #include <itkConnectedComponentImageFilter.h>
  27 #include <itkRelabelComponentImageFilter.h>
  28 #include <itkBinaryThresholdImageFilter.h>
  29 #include <itkPasteImageFilter.h>
  30 #include <itkStatisticsLabelMapFilter.h>
  31 #include <itkBinaryBallStructuringElement.h>
  32 #include <itkBinaryDilateImageFilter.h>
  33 #include <itkConstantPadImageFilter.h>
  34 #include <itkImageSliceIteratorWithIndex.h>
  35
  36 namespace clitk {
  37
  38   //--------------------------------------------------------------------
  39   template<class ImageType>
  40   void ComputeBBFromImageRegion(const ImageType * image,
  41                                 typename ImageType::RegionType region,
  42                                 typename itk::BoundingBox<unsigned long,
  43                                 ImageType::ImageDimension>::Pointer bb) {
  44     typedef typename ImageType::IndexType IndexType;
  45     IndexType firstIndex;
  46     IndexType lastIndex;
  47     for(unsigned int i=0; i<image->GetImageDimension(); i++) {
  48       firstIndex[i] = region.GetIndex()[i];
  49       lastIndex[i] = firstIndex[i]+region.GetSize()[i];
  50     }
  51
  52     typedef itk::BoundingBox<unsigned long,
  53                              ImageType::ImageDimension> BBType;
  54     typedef typename BBType::PointType PointType;
  55     PointType lastPoint;
  56     PointType firstPoint;
  57     image->TransformIndexToPhysicalPoint(firstIndex, firstPoint);
  58     image->TransformIndexToPhysicalPoint(lastIndex, lastPoint);
  59
  60     bb->SetMaximum(lastPoint);
  61     bb->SetMinimum(firstPoint);
  62   }
  63   //--------------------------------------------------------------------
  64
  65
  66   //--------------------------------------------------------------------
  67   template<int Dimension>
  68   void ComputeBBIntersection(typename itk::BoundingBox<unsigned long, Dimension>::Pointer bbo,
  69                              typename itk::BoundingBox<unsigned long, Dimension>::Pointer bbi1,
  70                              typename itk::BoundingBox<unsigned long, Dimension>::Pointer bbi2) {
  71
  72     typedef itk::BoundingBox<unsigned long, Dimension> BBType;
  73     typedef typename BBType::PointType PointType;
  74     PointType lastPoint;
  75     PointType firstPoint;
  76
  77     for(unsigned int i=0; i<Dimension; i++) {
  78       firstPoint[i] = std::max(bbi1->GetMinimum()[i],
  79                                bbi2->GetMinimum()[i]);
  80       lastPoint[i] = std::min(bbi1->GetMaximum()[i],
  81                               bbi2->GetMaximum()[i]);
  82     }
  83
  84     bbo->SetMaximum(lastPoint);
  85     bbo->SetMinimum(firstPoint);
  86   }
  87   //--------------------------------------------------------------------
  88
  89
  90   //--------------------------------------------------------------------
  91   template<class ImageType>
  92   void ComputeRegionFromBB(const ImageType * image,
  93                            const typename itk::BoundingBox<unsigned long,
  94                                                            ImageType::ImageDimension>::Pointer bb,
  95                            typename ImageType::RegionType & region) {
  96     // Types
  97     typedef typename ImageType::IndexType  IndexType;
  98     typedef typename ImageType::PointType  PointType;
  99     typedef typename ImageType::RegionType RegionType;
 100     typedef typename ImageType::SizeType   SizeType;
 101
 102     // Region starting point
 103     IndexType regionStart;
 104     PointType start = bb->GetMinimum();
 105     image->TransformPhysicalPointToIndex(start, regionStart);
 106
 107     // Region size
 108     SizeType regionSize;
 109     PointType maxs = bb->GetMaximum();
 110     PointType mins = bb->GetMinimum();
 111     for(unsigned int i=0; i<ImageType::ImageDimension; i++) {
 112       regionSize[i] = lrint((maxs[i] - mins[i])/image->GetSpacing()[i]);
 113     }
 114
 115     // Create region
 116     region.SetIndex(regionStart);
 117     region.SetSize(regionSize);
 118   }
 119   //--------------------------------------------------------------------
 120
 121   //--------------------------------------------------------------------
 122   template<class ImageType, class TMaskImageType>
 123   typename ImageType::Pointer
 124   SetBackground(const ImageType * input,
 125                 const TMaskImageType * mask,
 126                 typename TMaskImageType::PixelType maskBG,
 127                 typename ImageType::PixelType outValue,
 128                 bool inPlace) {
 129     typedef SetBackgroundImageFilter<ImageType, TMaskImageType, ImageType>
 130       SetBackgroundImageFilterType;
 131     typename SetBackgroundImageFilterType::Pointer setBackgroundFilter
 132       = SetBackgroundImageFilterType::New();
 133     //  if (inPlace) setBackgroundFilter->ReleaseDataFlagOn(); // No seg fault
 134     setBackgroundFilter->SetInPlace(inPlace); // This is important to keep memory low
 135     setBackgroundFilter->SetInput(input);
 136     setBackgroundFilter->SetInput2(mask);
 137     setBackgroundFilter->SetMaskValue(maskBG);
 138     setBackgroundFilter->SetOutsideValue(outValue);
 139     setBackgroundFilter->Update();
 140     return setBackgroundFilter->GetOutput();
 141   }
 142   //--------------------------------------------------------------------
 143
 144
 145   //--------------------------------------------------------------------
 146   template<class ImageType>
 147   int GetNumberOfConnectedComponentLabels(const ImageType * input,
 148                                           typename ImageType::PixelType BG,
 149                                           bool isFullyConnected) {
 150     // Connected Component label
 151     typedef itk::ConnectedComponentImageFilter<ImageType, ImageType> ConnectFilterType;
 152     typename ConnectFilterType::Pointer connectFilter = ConnectFilterType::New();
 153     connectFilter->SetInput(input);
 154     connectFilter->SetBackgroundValue(BG);
 155     connectFilter->SetFullyConnected(isFullyConnected);
 156     connectFilter->Update();
 157
 158     // Return result
 159     return connectFilter->GetObjectCount();
 160   }
 161   //--------------------------------------------------------------------
 162
 163   //--------------------------------------------------------------------
 164   /*
 165     Warning : in this cas, we consider outputType like inputType, not
 166     InternalImageType. Be sure it fits.
 167   */
 168   template<class ImageType>
 169   typename ImageType::Pointer
 170   Labelize(const ImageType * input,
 171            typename ImageType::PixelType BG,
 172            bool isFullyConnected,
 173            int minimalComponentSize) {
 174     // InternalImageType for storing large number of component
 175     typedef itk::Image<int, ImageType::ImageDimension> InternalImageType;
 176
 177     // Connected Component label
 178     typedef itk::ConnectedComponentImageFilter<ImageType, InternalImageType> ConnectFilterType;
 179     typename ConnectFilterType::Pointer connectFilter = ConnectFilterType::New();
 180     //  connectFilter->ReleaseDataFlagOn();
 181     connectFilter->SetInput(input);
 182     connectFilter->SetBackgroundValue(BG);
 183     connectFilter->SetFullyConnected(isFullyConnected);
 184
 185     // Sort by size and remove too small area.
 186     typedef itk::RelabelComponentImageFilter<InternalImageType, ImageType> RelabelFilterType;
 187     typename RelabelFilterType::Pointer relabelFilter = RelabelFilterType::New();
 188     //  relabelFilter->ReleaseDataFlagOn(); // if yes, fail when ExplosionControlledThresholdConnectedImageFilter ???
 189     relabelFilter->SetInput(connectFilter->GetOutput());
 190     relabelFilter->SetMinimumObjectSize(minimalComponentSize);
 191     relabelFilter->Update();
 192
 193     // Return result
 194     typename ImageType::Pointer output = relabelFilter->GetOutput();
 195     return output;
 196   }
 197   //--------------------------------------------------------------------
 198
 199
 200   //--------------------------------------------------------------------
 201   /*
 202     Warning : in this cas, we consider outputType like inputType, not
 203     InternalImageType. Be sure it fits.
 204   */
 205   template<class ImageType>
 206   typename ImageType::Pointer
 207   LabelizeAndCountNumberOfObjects(const ImageType * input,
 208                                   typename ImageType::PixelType BG,
 209                                   bool isFullyConnected,
 210                                   int minimalComponentSize,
 211                                   int & nb) {
 212     // InternalImageType for storing large number of component
 213     typedef itk::Image<int, ImageType::ImageDimension> InternalImageType;
 214
 215     // Connected Component label
 216     typedef itk::ConnectedComponentImageFilter<ImageType, InternalImageType> ConnectFilterType;
 217     typename ConnectFilterType::Pointer connectFilter = ConnectFilterType::New();
 218     //  connectFilter->ReleaseDataFlagOn();
 219     connectFilter->SetInput(input);
 220     connectFilter->SetBackgroundValue(BG);
 221     connectFilter->SetFullyConnected(isFullyConnected);
 222
 223     // Sort by size and remove too small area.
 224     typedef itk::RelabelComponentImageFilter<InternalImageType, ImageType> RelabelFilterType;
 225     typename RelabelFilterType::Pointer relabelFilter = RelabelFilterType::New();
 226     //  relabelFilter->ReleaseDataFlagOn(); // if yes, fail when ExplosionControlledThresholdConnectedImageFilter ???
 227     relabelFilter->SetInput(connectFilter->GetOutput());
 228     relabelFilter->SetMinimumObjectSize(minimalComponentSize);
 229     relabelFilter->Update();
 230
 231     nb = relabelFilter->GetNumberOfObjects();
 232     // DD(relabelFilter->GetOriginalNumberOfObjects());
 233     // DD(relabelFilter->GetSizeOfObjectsInPhysicalUnits()[0]);
 234
 235     // Return result
 236     typename ImageType::Pointer output = relabelFilter->GetOutput();
 237     return output;
 238   }
 239   //--------------------------------------------------------------------
 240
 241
 242
 243   //--------------------------------------------------------------------
 244   template<class ImageType>
 245   typename ImageType::Pointer
 246   RemoveLabels(const ImageType * input,
 247                typename ImageType::PixelType BG,
 248                std::vector<typename ImageType::PixelType> & labelsToRemove) {
 249     assert(labelsToRemove.size() != 0);
 250     typename ImageType::Pointer working_image;// = input;
 251     for (unsigned int i=0; i <labelsToRemove.size(); i++) {
 252       typedef SetBackgroundImageFilter<ImageType, ImageType> SetBackgroundImageFilterType;
 253       typename SetBackgroundImageFilterType::Pointer setBackgroundFilter = SetBackgroundImageFilterType::New();
 254       setBackgroundFilter->SetInput(input);
 255       setBackgroundFilter->SetInput2(input);
 256       setBackgroundFilter->SetMaskValue(labelsToRemove[i]);
 257       setBackgroundFilter->SetOutsideValue(BG);
 258       setBackgroundFilter->Update();
 259       working_image = setBackgroundFilter->GetOutput();
 260     }
 261     return working_image;
 262   }
 263   //--------------------------------------------------------------------
 264
 265
 266   //--------------------------------------------------------------------
 267   template<class ImageType>
 268   typename ImageType::Pointer
 269   KeepLabels(const ImageType * input,
 270              typename ImageType::PixelType BG,
 271              typename ImageType::PixelType FG,
 272              typename ImageType::PixelType firstKeep,
 273              typename ImageType::PixelType lastKeep,
 274              bool useLastKeep) {
 275     typedef itk::BinaryThresholdImageFilter<ImageType, ImageType> BinarizeFilterType;
 276     typename BinarizeFilterType::Pointer binarizeFilter = BinarizeFilterType::New();
 277     binarizeFilter->SetInput(input);
 278     binarizeFilter->SetLowerThreshold(firstKeep);
 279     if (useLastKeep) binarizeFilter->SetUpperThreshold(lastKeep);
 280     binarizeFilter->SetInsideValue(FG);
 281     binarizeFilter->SetOutsideValue(BG);
 282     binarizeFilter->Update();
 283     return binarizeFilter->GetOutput();
 284   }
 285   //--------------------------------------------------------------------
 286
 287
 288   //--------------------------------------------------------------------
 289   template<class ImageType>
 290   typename ImageType::Pointer
 291   LabelizeAndSelectLabels(const ImageType * input,
 292                           typename ImageType::PixelType BG,
 293                           typename ImageType::PixelType FG,
 294                           bool isFullyConnected,
 295                           int minimalComponentSize,
 296                           LabelizeParameters<typename ImageType::PixelType> * param)
 297   {
 298     typename ImageType::Pointer working_image;
 299     working_image = Labelize<ImageType>(input, BG, isFullyConnected, minimalComponentSize);
 300     working_image = RemoveLabels<ImageType>(working_image, BG, param->GetLabelsToRemove());
 301     working_image = KeepLabels<ImageType>(working_image,
 302                                           BG, FG,
 303                                           param->GetFirstKeep(),
 304                                           param->GetLastKeep(),
 305                                           param->GetUseLastKeep());
 306     return working_image;
 307   }
 308   //--------------------------------------------------------------------
 309
 310
 311   //--------------------------------------------------------------------
 312   template<class ImageType>
 313   typename ImageType::Pointer
 314   ResizeImageLike(const ImageType * input,
 315                   const itk::ImageBase<ImageType::ImageDimension> * like,
 316                   typename ImageType::PixelType backgroundValue)
 317   {
 318     typedef CropLikeImageFilter<ImageType> CropFilterType;
 319     typename CropFilterType::Pointer cropFilter = CropFilterType::New();
 320     cropFilter->SetInput(input);
 321     cropFilter->SetCropLikeImage(like);
 322     cropFilter->SetBackgroundValue(backgroundValue);
 323     cropFilter->Update();
 324     return cropFilter->GetOutput();
 325   }
 326   //--------------------------------------------------------------------
 327
 328
 329   //--------------------------------------------------------------------
 330   template<class MaskImageType>
 331   typename MaskImageType::Pointer
 332   SliceBySliceRelativePosition(const MaskImageType * input,
 333                                const MaskImageType * object,
 334                                int direction,
 335                                double threshold,
 336                                std::string orientation,
 337                                bool uniqueConnectedComponent,
 338                                double spacing,
 339                                bool inverseflag)
 340   {
 341     typedef SliceBySliceRelativePositionFilter<MaskImageType> SliceRelPosFilterType;
 342     typename SliceRelPosFilterType::Pointer sliceRelPosFilter = SliceRelPosFilterType::New();
 343     sliceRelPosFilter->VerboseStepFlagOff();
 344     sliceRelPosFilter->WriteStepFlagOff();
 345     sliceRelPosFilter->SetInput(input);
 346     sliceRelPosFilter->SetInputObject(object);
 347     sliceRelPosFilter->SetDirection(direction);
 348     sliceRelPosFilter->SetFuzzyThreshold(threshold);
 349     sliceRelPosFilter->AddOrientationTypeString(orientation);
 350     sliceRelPosFilter->SetIntermediateSpacingFlag((spacing != -1));
 351     sliceRelPosFilter->SetIntermediateSpacing(spacing);
 352     sliceRelPosFilter->SetUniqueConnectedComponentBySlice(uniqueConnectedComponent);
 353     sliceRelPosFilter->SetInverseOrientationFlag(inverseflag);
 354     //  sliceRelPosFilter->SetAutoCropFlag(true); ??
 355     sliceRelPosFilter->Update();
 356     return sliceRelPosFilter->GetOutput();
 357   }
 358   //--------------------------------------------------------------------
 359
 360   //--------------------------------------------------------------------
 361   template<class ImageType>
 362   bool
 363   FindExtremaPointInAGivenDirection(const ImageType * input,
 364                                     typename ImageType::PixelType bg,
 365                                     int direction, bool opposite,
 366                                     typename ImageType::PointType & point)
 367   {
 368     typename ImageType::PointType dummy;
 369     return FindExtremaPointInAGivenDirection(input, bg, direction,
 370                                              opposite, dummy, 0, point);
 371   }
 372   //--------------------------------------------------------------------
 373
 374   //--------------------------------------------------------------------
 375   template<class ImageType>
 376   bool
 377   FindExtremaPointInAGivenDirection(const ImageType * input,
 378                                     typename ImageType::PixelType bg,
 379                                     int direction, bool opposite,
 380                                     typename ImageType::PointType refpoint,
 381                                     double distanceMax,
 382                                     typename ImageType::PointType & point)
 383   {
 384     /*
 385       loop over input pixels, store the index in the fg that is max
 386       according to the given direction.
 387     */
 388     typedef itk::ImageRegionConstIteratorWithIndex<ImageType> IteratorType;
 389     IteratorType iter(input, input->GetLargestPossibleRegion());
 390     iter.GoToBegin();
 391     typename ImageType::IndexType max = input->GetLargestPossibleRegion().GetIndex();
 392     if (opposite) max = max+input->GetLargestPossibleRegion().GetSize();
 393     bool found=false;
 394     while (!iter.IsAtEnd()) {
 395       if (iter.Get() != bg) {
 396         bool test = iter.GetIndex()[direction] >  max[direction];
 397         if (opposite) test = !test;
 398         if (test) {
 399           typename ImageType::PointType p;
 400           input->TransformIndexToPhysicalPoint(iter.GetIndex(), p);
 401           if ((distanceMax==0) || (p.EuclideanDistanceTo(refpoint) < distanceMax)) {
 402             max = iter.GetIndex();
 403             found = true;
 404           }
 405         }
 406       }
 407       ++iter;
 408     }
 409     if (!found) return false;
 410     input->TransformIndexToPhysicalPoint(max, point);
 411     return true;
 412   }
 413   //--------------------------------------------------------------------
 414
 415
 416   //--------------------------------------------------------------------
 417   template<class ImageType>
 418   typename ImageType::Pointer
 419   CropImageAbove(const ImageType * image,
 420                  int dim, double min, bool autoCrop,
 421                  typename ImageType::PixelType BG)
 422   {
 423     return CropImageAlongOneAxis<ImageType>(image, dim,
 424                                             image->GetOrigin()[dim],
 425                                             min,
 426                                             autoCrop, BG);
 427   }
 428   //--------------------------------------------------------------------
 429
 430
 431   //--------------------------------------------------------------------
 432   template<class ImageType>
 433   typename ImageType::Pointer
 434   CropImageBelow(const ImageType * image,
 435                  int dim, double max, bool autoCrop,
 436                  typename ImageType::PixelType BG)
 437   {
 438     typename ImageType::PointType p;
 439     image->TransformIndexToPhysicalPoint(image->GetLargestPossibleRegion().GetIndex()+
 440                                          image->GetLargestPossibleRegion().GetSize(), p);
 441     return CropImageAlongOneAxis<ImageType>(image, dim, max, p[dim], autoCrop, BG);
 442   }
 443   //--------------------------------------------------------------------
 444
 445
 446   //--------------------------------------------------------------------
 447   template<class ImageType>
 448   typename ImageType::Pointer
 449   CropImageAlongOneAxis(const ImageType * image,
 450                         int dim, double min, double max,
 451                         bool autoCrop, typename ImageType::PixelType BG)
 452   {
 453     // Compute region size
 454     typename ImageType::RegionType region;
 455     typename ImageType::SizeType size = image->GetLargestPossibleRegion().GetSize();
 456     typename ImageType::PointType p = image->GetOrigin();
 457     p[dim] = min;
 458     typename ImageType::IndexType start;
 459     image->TransformPhysicalPointToIndex(p, start);
 460     p[dim] = max;
 461     typename ImageType::IndexType end;
 462     image->TransformPhysicalPointToIndex(p, end);
 463     size[dim] = fabs(end[dim]-start[dim]);
 464     region.SetIndex(start);
 465     region.SetSize(size);
 466
 467     // Perform Crop
 468     typedef itk::RegionOfInterestImageFilter<ImageType, ImageType> CropFilterType;
 469     typename CropFilterType::Pointer cropFilter = CropFilterType::New();
 470     cropFilter->SetInput(image);
 471     cropFilter->SetRegionOfInterest(region);
 472     cropFilter->Update();
 473     typename ImageType::Pointer result = cropFilter->GetOutput();
 474
 475     // Auto Crop
 476     if (autoCrop) {
 477       result = AutoCrop<ImageType>(result, BG);
 478     }
 479     return result;
 480   }
 481   //--------------------------------------------------------------------
 482
 483
 484   //--------------------------------------------------------------------
 485   template<class ImageType>
 486   void
 487   ComputeCentroids(const ImageType * image,
 488                    typename ImageType::PixelType BG,
 489                    std::vector<typename ImageType::PointType> & centroids)
 490   {
 491     typedef long LabelType;
 492     static const unsigned int Dim = ImageType::ImageDimension;
 493     typedef itk::ShapeLabelObject< LabelType, Dim > LabelObjectType;
 494     typedef itk::LabelMap< LabelObjectType > LabelMapType;
 495     typedef itk::LabelImageToLabelMapFilter<ImageType, LabelMapType> ImageToMapFilterType;
 496     typename ImageToMapFilterType::Pointer imageToLabelFilter = ImageToMapFilterType::New();
 497     typedef itk::ShapeLabelMapFilter<LabelMapType, ImageType> ShapeFilterType;
 498     typename ShapeFilterType::Pointer statFilter = ShapeFilterType::New();
 499     imageToLabelFilter->SetBackgroundValue(BG);
 500     imageToLabelFilter->SetInput(image);
 501     statFilter->SetInput(imageToLabelFilter->GetOutput());
 502     statFilter->Update();
 503     typename LabelMapType::Pointer labelMap = statFilter->GetOutput();
 504
 505     centroids.clear();
 506     typename ImageType::PointType dummy;
 507     centroids.push_back(dummy); // label 0 -> no centroid, use dummy point
 508     for(uint i=1; i<labelMap->GetNumberOfLabelObjects()+1; i++) {
 509       centroids.push_back(labelMap->GetLabelObject(i)->GetCentroid());
 510     }
 511   }
 512   //--------------------------------------------------------------------
 513
 514
 515   //--------------------------------------------------------------------
 516   template<class ImageType>
 517   void
 518   ExtractSlices(const ImageType * image, int direction,
 519                 std::vector<typename itk::Image<typename ImageType::PixelType,
 520                                                 ImageType::ImageDimension-1>::Pointer > & slices)
 521   {
 522     typedef ExtractSliceFilter<ImageType> ExtractSliceFilterType;
 523     typedef typename ExtractSliceFilterType::SliceType SliceType;
 524     typename ExtractSliceFilterType::Pointer
 525       extractSliceFilter = ExtractSliceFilterType::New();
 526     extractSliceFilter->SetInput(image);
 527     extractSliceFilter->SetDirection(direction);
 528     extractSliceFilter->Update();
 529     extractSliceFilter->GetOutputSlices(slices);
 530   }
 531   //--------------------------------------------------------------------
 532
 533
 534   //--------------------------------------------------------------------
 535   template<class ImageType>
 536   typename ImageType::Pointer
 537   JoinSlices(std::vector<typename itk::Image<typename ImageType::PixelType,
 538                                              ImageType::ImageDimension-1>::Pointer > & slices,
 539              const ImageType * input,
 540              int direction) {
 541     typedef typename itk::Image<typename ImageType::PixelType, ImageType::ImageDimension-1> SliceType;
 542     typedef itk::JoinSeriesImageFilter<SliceType, ImageType> JoinSeriesFilterType;
 543     typename JoinSeriesFilterType::Pointer joinFilter = JoinSeriesFilterType::New();
 544     joinFilter->SetOrigin(input->GetOrigin()[direction]);
 545     joinFilter->SetSpacing(input->GetSpacing()[direction]);
 546     for(unsigned int i=0; i<slices.size(); i++) {
 547       joinFilter->PushBackInput(slices[i]);
 548     }
 549     joinFilter->Update();
 550     return joinFilter->GetOutput();
 551   }
 552   //--------------------------------------------------------------------
 553
 554
 555   //--------------------------------------------------------------------
 556   template<class ImageType>
 557   void
 558   PointsUtils<ImageType>::Convert2DTo3D(const PointType2D & p2D,
 559                                         const ImageType * image,
 560                                         const int slice,
 561                                         PointType3D & p3D)
 562   {
 563     IndexType3D index3D;
 564     index3D[0] = index3D[1] = 0;
 565     index3D[2] = image->GetLargestPossibleRegion().GetIndex()[2]+slice;
 566     image->TransformIndexToPhysicalPoint(index3D, p3D);
 567     p3D[0] = p2D[0];
 568     p3D[1] = p2D[1];
 569     //  p3D[2] = p[2];//(image->GetLargestPossibleRegion().GetIndex()[2]+slice)*image->GetSpacing()[2]
 570     //    + image->GetOrigin()[2];
 571   }
 572   //--------------------------------------------------------------------
 573
 574
 575   //--------------------------------------------------------------------
 576   template<class ImageType>
 577   void
 578   PointsUtils<ImageType>::Convert2DTo3DList(const MapPoint2DType & map,
 579                                             const ImageType * image,
 580                                             VectorPoint3DType & list)
 581   {
 582     typename MapPoint2DType::const_iterator iter = map.begin();
 583     while (iter != map.end()) {
 584       PointType3D p;
 585       Convert2DTo3D(iter->second, image, iter->first, p);
 586       list.push_back(p);
 587       ++iter;
 588     }
 589   }
 590   //--------------------------------------------------------------------
 591
 592   //--------------------------------------------------------------------
 593   template<class ImageType>
 594   void
 595   WriteListOfLandmarks(std::vector<typename ImageType::PointType> points,
 596                        std::string filename)
 597   {
 598     std::ofstream os;
 599     openFileForWriting(os, filename);
 600     os << "LANDMARKS1" << std::endl;
 601     for(uint i=0; i<points.size(); i++) {
 602       const typename ImageType::PointType & p = points[i];
 603       // Write it in the file
 604       os << i << " " << p[0] << " " << p[1] << " " << p[2] << " 0 0 " << std::endl;
 605     }
 606     os.close();
 607   }
 608   //--------------------------------------------------------------------
 609
 610
 611   //--------------------------------------------------------------------
 612   template<class ImageType>
 613   typename ImageType::Pointer
 614   Dilate(const ImageType * image, double radiusInMM,
 615          typename ImageType::PixelType BG,
 616          typename ImageType::PixelType FG,
 617          bool extendSupport)
 618   {
 619     typename ImageType::SizeType r;
 620     for(uint i=0; i<ImageType::ImageDimension; i++)
 621       r[i] = (uint)lrint(radiusInMM/image->GetSpacing()[i]);
 622     return Dilate<ImageType>(image, r, BG, FG, extendSupport);
 623   }
 624   //--------------------------------------------------------------------
 625
 626
 627   //--------------------------------------------------------------------
 628   template<class ImageType>
 629   typename ImageType::Pointer
 630   Dilate(const ImageType * image, typename ImageType::PointType radiusInMM,
 631          typename ImageType::PixelType BG,
 632          typename ImageType::PixelType FG,
 633          bool extendSupport)
 634   {
 635     typename ImageType::SizeType r;
 636     for(uint i=0; i<ImageType::ImageDimension; i++)
 637       r[i] = (uint)lrint(radiusInMM[i]/image->GetSpacing()[i]);
 638     return Dilate<ImageType>(image, r, BG, FG, extendSupport);
 639   }
 640   //--------------------------------------------------------------------
 641
 642
 643   //--------------------------------------------------------------------
 644   template<class ImageType>
 645   typename ImageType::Pointer
 646   Dilate(const ImageType * image, typename ImageType::SizeType radius,
 647          typename ImageType::PixelType BG,
 648          typename ImageType::PixelType FG,
 649          bool extendSupport)
 650   {
 651     // Create kernel for dilatation
 652     typedef itk::BinaryBallStructuringElement<typename ImageType::PixelType,
 653                                               ImageType::ImageDimension> KernelType;
 654     KernelType structuringElement;
 655     structuringElement.SetRadius(radius);
 656     structuringElement.CreateStructuringElement();
 657
 658     typename ImageType::Pointer output;
 659     if (extendSupport) {
 660       typedef itk::ConstantPadImageFilter<ImageType, ImageType> PadFilterType;
 661       typename PadFilterType::Pointer padFilter = PadFilterType::New();
 662       padFilter->SetInput(image);
 663       typename ImageType::SizeType lower;
 664       typename ImageType::SizeType upper;
 665       for(uint i=0; i<3; i++) {
 666         lower[i] = upper[i] = 2*(radius[i]+1);
 667       }
 668       padFilter->SetPadLowerBound(lower);
 669       padFilter->SetPadUpperBound(upper);
 670       padFilter->Update();
 671       output = padFilter->GetOutput();
 672     }
 673
 674     // Dilate  filter
 675     typedef itk::BinaryDilateImageFilter<ImageType, ImageType , KernelType> DilateFilterType;
 676     typename DilateFilterType::Pointer dilateFilter = DilateFilterType::New();
 677     dilateFilter->SetBackgroundValue(BG);
 678     dilateFilter->SetForegroundValue(FG);
 679     dilateFilter->SetBoundaryToForeground(false);
 680     dilateFilter->SetKernel(structuringElement);
 681     dilateFilter->SetInput(output);
 682     dilateFilter->Update();
 683     return dilateFilter->GetOutput();
 684   }
 685   //--------------------------------------------------------------------
 686
 687
 688   //--------------------------------------------------------------------
 689   template<class ValueType, class VectorType>
 690   void ConvertOption(std::string optionName, uint given,
 691                      ValueType * values, VectorType & p,
 692                      uint dim, bool required)
 693   {
 694     if (required && (given == 0)) {
 695       clitkExceptionMacro("The option --" << optionName << " must be set and have 1 or "
 696                           << dim << " values.");
 697     }
 698     if (given == 1) {
 699       for(uint i=0; i<dim; i++) p[i] = values[0];
 700       return;
 701     }
 702     if (given == dim) {
 703       for(uint i=0; i<dim; i++) p[i] = values[i];
 704       return;
 705     }
 706     if (given == 0) return;
 707     clitkExceptionMacro("The option --" << optionName << " must have 1 or "
 708                         << dim << " values.");
 709   }
 710   //--------------------------------------------------------------------
 711
 712
 713   //--------------------------------------------------------------------
 714   /*
 715     http://www.gamedev.net/community/forums/topic.asp?topic_id=542870
 716     Assuming the points are (Ax,Ay) (Bx,By) and (Cx,Cy), you need to compute:
 717     (Bx - Ax) * (Cy - Ay) - (By - Ay) * (Cx - Ax)
 718     This will equal zero if the point C is on the line formed by
 719     points A and B, and will have a different sign depending on the
 720     side. Which side this is depends on the orientation of your (x,y)
 721     coordinates, but you can plug test values for A,B and C into this
 722     formula to determine whether negative values are to the left or to
 723     the right.
 724     => to accelerate, start with formula, when change sign -> stop and fill
 725
 726     offsetToKeep = is used to determine which side of the line we
 727     keep. The point along the mainDirection but 'offsetToKeep' mm away
 728     is kept.
 729
 730   */
 731   template<class ImageType>
 732   void
 733   SliceBySliceSetBackgroundFromLineSeparation(ImageType * input,
 734                                               std::vector<typename ImageType::PointType> & lA,
 735                                               std::vector<typename ImageType::PointType> & lB,
 736                                               typename ImageType::PixelType BG,
 737                                               int mainDirection,
 738                                               double offsetToKeep)
 739   {
 740     typedef itk::ImageSliceIteratorWithIndex<ImageType> SliceIteratorType;
 741     SliceIteratorType siter = SliceIteratorType(input,
 742                                                 input->GetLargestPossibleRegion());
 743     siter.SetFirstDirection(0);
 744     siter.SetSecondDirection(1);
 745     siter.GoToBegin();
 746     uint i=0;
 747     typename ImageType::PointType A;
 748     typename ImageType::PointType B;
 749     typename ImageType::PointType C;
 750     assert(lA.size() == lB.size());
 751     while ((i<lA.size()) && (!siter.IsAtEnd())) {
 752       // Check that the current slice correspond to the current point
 753       input->TransformIndexToPhysicalPoint(siter.GetIndex(), C);
 754       // DD(C);
 755       // DD(i);
 756       // DD(lA[i]);
 757       if ((fabs(C[2] - lA[i][2]))>0.01) { // is !equal with a tolerance of 0.01 mm
 758       }
 759       else {
 760         // Define A,B,C points
 761         A = lA[i];
 762         B = lB[i];
 763         C = A;
 764         // DD(A);
 765         // DD(B);
 766         // DD(C);
 767
 768         // Check that the line is not a point (A=B)
 769         bool p = (A[0] == B[0]) && (A[1] == B[1]);
 770
 771         if (!p) {
 772           C[mainDirection] += offsetToKeep; // I know I must keep this point
 773           double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
 774           bool isPositive = s<0;
 775           while (!siter.IsAtEndOfSlice()) {
 776             while (!siter.IsAtEndOfLine()) {
 777               // Very slow, I know ... but image should be very small
 778               input->TransformIndexToPhysicalPoint(siter.GetIndex(), C);
 779               double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
 780               if (s == 0) siter.Set(BG); // on the line, we decide to remove
 781               if (isPositive) {
 782                 if (s > 0) siter.Set(BG);
 783               }
 784               else {
 785                 if (s < 0) siter.Set(BG);
 786               }
 787               ++siter;
 788             }
 789             siter.NextLine();
 790           } // end loop slice
 791         }
 792
 793         ++i;
 794       } // End of current slice
 795       siter.NextSlice();
 796     }
 797   }
 798   //--------------------------------------------------------------------
 799
 800   //--------------------------------------------------------------------
 801   template<class ImageType>
 802   void
 803   AndNot(ImageType * input,
 804          const ImageType * object,
 805          typename ImageType::PixelType BG)
 806   {
 807     typename ImageType::Pointer o;
 808     bool resized=false;
 809     if (!clitk::HaveSameSizeAndSpacing<ImageType, ImageType>(input, object)) {
 810       o = clitk::ResizeImageLike<ImageType>(object, input, BG);
 811       resized = true;
 812     }
 813
 814     typedef clitk::BooleanOperatorLabelImageFilter<ImageType> BoolFilterType;
 815     typename BoolFilterType::Pointer boolFilter = BoolFilterType::New();
 816     boolFilter->InPlaceOn();
 817     boolFilter->SetInput1(input);
 818     if (resized) boolFilter->SetInput2(o);
 819     else boolFilter->SetInput2(object);
 820     boolFilter->SetBackgroundValue1(BG);
 821     boolFilter->SetBackgroundValue2(BG);
 822     boolFilter->SetOperationType(BoolFilterType::AndNot);
 823     boolFilter->Update();
 824   }
 825   //--------------------------------------------------------------------
 826
 827
 828   //--------------------------------------------------------------------
 829   template<class ImageType>
 830   typename ImageType::Pointer
 831   Binarize(const ImageType * input,
 832            typename ImageType::PixelType lower,
 833            typename ImageType::PixelType upper,
 834            typename ImageType::PixelType BG=0,
 835            typename ImageType::PixelType FG=1)
 836   {
 837     typedef itk::BinaryThresholdImageFilter<ImageType, ImageType> BinaryThresholdFilterType;
 838     typename BinaryThresholdFilterType::Pointer binarizeFilter = BinaryThresholdFilterType::New();
 839     binarizeFilter->SetInput(input);
 840     binarizeFilter->SetLowerThreshold(lower);
 841     binarizeFilter->SetUpperThreshold(upper);
 842     binarizeFilter->SetInsideValue(FG);
 843     binarizeFilter->SetOutsideValue(BG);
 844     binarizeFilter->Update();
 845     return binarizeFilter->GetOutput();
 846   }
 847   //--------------------------------------------------------------------
 848
 849
 850   //--------------------------------------------------------------------
 851   template<class ImageType>
 852   void
 853   GetMinMaxPointPosition(const ImageType * input,
 854                          typename ImageType::PointType & min,
 855                          typename ImageType::PointType & max)
 856   {
 857     typename ImageType::IndexType index = input->GetLargestPossibleRegion().GetIndex();
 858     input->TransformIndexToPhysicalPoint(index, min);
 859     index = index+input->GetLargestPossibleRegion().GetSize();
 860     input->TransformIndexToPhysicalPoint(index, max);
 861   }
 862   //--------------------------------------------------------------------
 863
 864
 865   //--------------------------------------------------------------------
 866   template<class ImageType>
 867   typename ImageType::PointType
 868   FindExtremaPointInAGivenLine(const ImageType * input,
 869                                int dimension,
 870                                bool inverse,
 871                                typename ImageType::PointType p,
 872                                typename ImageType::PixelType BG,
 873                                double distanceMax)
 874   {
 875     // Which direction ?  Increasing or decreasing.
 876     int d=1;
 877     if (inverse) d=-1;
 878
 879     // Transform to pixel index
 880     typename ImageType::IndexType index;
 881     input->TransformPhysicalPointToIndex(p, index);
 882
 883     // Loop while inside the mask;
 884     while (input->GetPixel(index) != BG) {
 885       index[dimension] += d;
 886     }
 887
 888     // Transform back to Physical Units
 889     typename ImageType::PointType result;
 890     input->TransformIndexToPhysicalPoint(index, result);
 891
 892     // Check that is is not too far away
 893     double distance = p.EuclideanDistanceTo(result);
 894     if (distance > distanceMax) {
 895       result = p; // Get back to initial value
 896     }
 897
 898     return result;
 899   }
 900   //--------------------------------------------------------------------
 901
 902
 903   //--------------------------------------------------------------------
 904   template<class PointType>
 905   bool
 906   IsOnTheSameLineSide(PointType C, PointType A, PointType B, PointType like)
 907   {
 908     // Look at the position of point 'like' according to the AB line
 909     double s = (B[0] - A[0]) * (like[1] - A[1]) - (B[1] - A[1]) * (like[0] - A[0]);
 910     bool negative = s<0;
 911
 912     // Look the C position
 913     s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
 914
 915     if (negative && (s<=0)) return true;
 916     if (!negative && (s>=0)) return true;
 917     return false;
 918   }
 919   //--------------------------------------------------------------------
 920
 921
 922 } // end of namespace
 923