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
  24 // itk
  25 #include <itkConnectedComponentImageFilter.h>
  26 #include <itkRelabelComponentImageFilter.h>
  27 #include <itkBinaryThresholdImageFilter.h>
  28 #include <itkPasteImageFilter.h>
  29
  30 //--------------------------------------------------------------------
  31 template<class ImageType>
  32 void clitk::ComputeBBFromImageRegion(typename ImageType::Pointer image,
  33                                      typename ImageType::RegionType region,
  34                                      typename itk::BoundingBox<unsigned long,
  35                                                                ImageType::ImageDimension>::Pointer bb) {
  36   typedef typename ImageType::IndexType IndexType;
  37   IndexType firstIndex;
  38   IndexType lastIndex;
  39   for(unsigned int i=0; i<image->GetImageDimension(); i++) {
  40     firstIndex[i] = region.GetIndex()[i];
  41     lastIndex[i] = firstIndex[i]+region.GetSize()[i];
  42   }
  43
  44   typedef itk::BoundingBox<unsigned long,
  45                            ImageType::ImageDimension> BBType;
  46   typedef typename BBType::PointType PointType;
  47   PointType lastPoint;
  48   PointType firstPoint;
  49   image->TransformIndexToPhysicalPoint(firstIndex, firstPoint);
  50   image->TransformIndexToPhysicalPoint(lastIndex, lastPoint);
  51
  52   bb->SetMaximum(lastPoint);
  53   bb->SetMinimum(firstPoint);
  54 }
  55 //--------------------------------------------------------------------
  56
  57
  58 //--------------------------------------------------------------------
  59 template<int Dimension>
  60 void clitk::ComputeBBIntersection(typename itk::BoundingBox<unsigned long, Dimension>::Pointer bbo,
  61                                   typename itk::BoundingBox<unsigned long, Dimension>::Pointer bbi1,
  62                                   typename itk::BoundingBox<unsigned long, Dimension>::Pointer bbi2) {
  63
  64   typedef itk::BoundingBox<unsigned long, Dimension> BBType;
  65   typedef typename BBType::PointType PointType;
  66   PointType lastPoint;
  67   PointType firstPoint;
  68
  69   for(unsigned int i=0; i<Dimension; i++) {
  70     firstPoint[i] = std::max(bbi1->GetMinimum()[i],
  71                              bbi2->GetMinimum()[i]);
  72     lastPoint[i] = std::min(bbi1->GetMaximum()[i],
  73                             bbi2->GetMaximum()[i]);
  74   }
  75
  76   bbo->SetMaximum(lastPoint);
  77   bbo->SetMinimum(firstPoint);
  78 }
  79 //--------------------------------------------------------------------
  80
  81
  82 //--------------------------------------------------------------------
  83 template<class ImageType>
  84 void clitk::ComputeRegionFromBB(typename ImageType::Pointer image,
  85                                 const typename itk::BoundingBox<unsigned long,
  86                                 ImageType::ImageDimension>::Pointer bb,
  87                                 typename ImageType::RegionType & region) {
  88   // Types
  89   typedef typename ImageType::IndexType  IndexType;
  90   typedef typename ImageType::PointType  PointType;
  91   typedef typename ImageType::RegionType RegionType;
  92   typedef typename ImageType::SizeType   SizeType;
  93
  94   // Region starting point
  95   IndexType regionStart;
  96   PointType start = bb->GetMinimum();
  97   image->TransformPhysicalPointToIndex(start, regionStart);
  98
  99   // Region size
 100   SizeType regionSize;
 101   PointType maxs = bb->GetMaximum();
 102   PointType mins = bb->GetMinimum();
 103   for(unsigned int i=0; i<ImageType::ImageDimension; i++) {
 104     // DD(maxs[i]);
 105     // DD(mins[i]);
 106     // DD((maxs[i] - mins[i])/image->GetSpacing()[i]);
 107     regionSize[i] = lrint((maxs[i] - mins[i])/image->GetSpacing()[i]);
 108     // DD(regionSize[i]);
 109   }
 110
 111   // Create region
 112   region.SetIndex(regionStart);
 113   region.SetSize(regionSize);
 114 }
 115 //--------------------------------------------------------------------
 116
 117 //--------------------------------------------------------------------
 118 template<class ImageType, class TMaskImageType>
 119 typename ImageType::Pointer
 120 clitk::SetBackground(//typename ImageType::ConstPointer input,
 121                      const ImageType * input,
 122                      const TMaskImageType * mask,
 123                      typename TMaskImageType::PixelType maskBG,
 124                      typename ImageType::PixelType outValue) {
 125   typedef clitk::SetBackgroundImageFilter<ImageType, TMaskImageType, ImageType> SetBackgroundImageFilterType;
 126   typename SetBackgroundImageFilterType::Pointer setBackgroundFilter = SetBackgroundImageFilterType::New();
 127   setBackgroundFilter->SetInput(input);
 128   setBackgroundFilter->SetInput2(mask);
 129   setBackgroundFilter->SetMaskValue(maskBG);
 130   setBackgroundFilter->SetOutsideValue(outValue);
 131   setBackgroundFilter->Update();
 132   return setBackgroundFilter->GetOutput();
 133 }
 134 //--------------------------------------------------------------------
 135
 136
 137 //--------------------------------------------------------------------
 138 template<class ImageType>
 139 int clitk::GetNumberOfConnectedComponentLabels(typename ImageType::Pointer input,
 140                                                typename ImageType::PixelType BG,
 141                                                bool isFullyConnected) {
 142   // Connected Component label
 143   typedef itk::ConnectedComponentImageFilter<ImageType, ImageType> ConnectFilterType;
 144   typename ConnectFilterType::Pointer connectFilter = ConnectFilterType::New();
 145   connectFilter->SetInput(input);
 146   connectFilter->SetBackgroundValue(BG);
 147   connectFilter->SetFullyConnected(isFullyConnected);
 148   connectFilter->Update();
 149
 150   // Return result
 151   return connectFilter->GetObjectCount();
 152 }
 153 //--------------------------------------------------------------------
 154
 155 //--------------------------------------------------------------------
 156 template<class ImageType>
 157 typename ImageType::Pointer
 158 clitk::Labelize(const ImageType * input,
 159                 typename ImageType::PixelType BG,
 160                 bool isFullyConnected,
 161                 int minimalComponentSize) {
 162   // InternalImageType for storing large number of component
 163   typedef itk::Image<int, ImageType::ImageDimension> InternalImageType;
 164
 165   // Connected Component label
 166   typedef itk::ConnectedComponentImageFilter<ImageType, InternalImageType> ConnectFilterType;
 167   typename ConnectFilterType::Pointer connectFilter = ConnectFilterType::New();
 168   connectFilter->SetInput(input);
 169   connectFilter->SetBackgroundValue(BG);
 170   connectFilter->SetFullyConnected(isFullyConnected);
 171
 172   // Sort by size and remove too small area.
 173   typedef itk::RelabelComponentImageFilter<InternalImageType, ImageType> RelabelFilterType;
 174   typename RelabelFilterType::Pointer relabelFilter = RelabelFilterType::New();
 175   relabelFilter->InPlaceOn();
 176   relabelFilter->SetInput(connectFilter->GetOutput());
 177   relabelFilter->SetMinimumObjectSize(minimalComponentSize);
 178   relabelFilter->Update();
 179
 180   // Return result
 181   return relabelFilter->GetOutput();
 182 }
 183 //--------------------------------------------------------------------
 184
 185
 186 //--------------------------------------------------------------------
 187 template<class ImageType>
 188 typename ImageType::Pointer
 189 clitk::RemoveLabels(typename ImageType::Pointer input,
 190                     typename ImageType::PixelType BG,
 191                     std::vector<typename ImageType::PixelType> & labelsToRemove) {
 192   typename ImageType::Pointer working_image = input;
 193   for (unsigned int i=0; i <labelsToRemove.size(); i++) {
 194     typedef clitk::SetBackgroundImageFilter<ImageType, ImageType> SetBackgroundImageFilterType;
 195     typename SetBackgroundImageFilterType::Pointer setBackgroundFilter = SetBackgroundImageFilterType::New();
 196     setBackgroundFilter->SetInput(input);
 197     setBackgroundFilter->SetInput2(input);
 198     setBackgroundFilter->SetMaskValue(labelsToRemove[i]);
 199     setBackgroundFilter->SetOutsideValue(BG);
 200     setBackgroundFilter->Update();
 201     working_image = setBackgroundFilter->GetOutput();
 202   }
 203   return working_image;
 204 }
 205 //--------------------------------------------------------------------
 206
 207
 208 //--------------------------------------------------------------------
 209 template<class ImageType>
 210 typename ImageType::Pointer
 211 clitk::KeepLabels(const ImageType * input,
 212                   typename ImageType::PixelType BG,
 213                   typename ImageType::PixelType FG,
 214                   typename ImageType::PixelType firstKeep,
 215                   typename ImageType::PixelType lastKeep,
 216                   bool useLastKeep) {
 217   typedef itk::BinaryThresholdImageFilter<ImageType, ImageType> BinarizeFilterType;
 218   typename BinarizeFilterType::Pointer binarizeFilter = BinarizeFilterType::New();
 219   binarizeFilter->SetInput(input);
 220   binarizeFilter->SetLowerThreshold(firstKeep);
 221   if (useLastKeep) binarizeFilter->SetUpperThreshold(lastKeep);
 222   binarizeFilter->SetInsideValue(FG);
 223   binarizeFilter->SetOutsideValue(BG);
 224   binarizeFilter->Update();
 225   return binarizeFilter->GetOutput();
 226 }
 227 //--------------------------------------------------------------------
 228
 229
 230 //--------------------------------------------------------------------
 231 template<class ImageType>
 232 typename ImageType::Pointer
 233 clitk::LabelizeAndSelectLabels(typename ImageType::Pointer input,
 234                                typename ImageType::PixelType BG,
 235                                typename ImageType::PixelType FG,
 236                                bool isFullyConnected,
 237                                int minimalComponentSize,
 238                                LabelizeParameters<typename ImageType::PixelType> * param)
 239 {
 240   typename ImageType::Pointer working_image;
 241   working_image = Labelize<ImageType>(input, BG, isFullyConnected, minimalComponentSize);
 242   working_image = RemoveLabels<ImageType>(working_image, BG, param->GetLabelsToRemove());
 243   working_image = KeepLabels<ImageType>(working_image,
 244                                         BG, FG,
 245                                         param->GetFirstKeep(),
 246                                         param->GetLastKeep(),
 247                                         param->GetUseLastKeep());
 248   return working_image;
 249 }
 250 //--------------------------------------------------------------------
 251
 252
 253 //--------------------------------------------------------------------
 254 template<class ImageType>
 255 typename ImageType::Pointer
 256 clitk::ResizeImageLike(typename ImageType::Pointer input,
 257                        typename ImageType::Pointer like,
 258                        typename ImageType::PixelType backgroundValue)
 259 {
 260   typedef clitk::CropLikeImageFilter<ImageType> CropFilterType;
 261   typename CropFilterType::Pointer cropFilter = CropFilterType::New();
 262   cropFilter->SetInput(input);
 263   cropFilter->SetCropLikeImage(like);
 264   cropFilter->SetBackgroundValue(backgroundValue);
 265   cropFilter->Update();
 266   return cropFilter->GetOutput();
 267 }
 268 //--------------------------------------------------------------------
 269
 270
 271 //--------------------------------------------------------------------
 272 template<class MaskImageType>
 273 typename MaskImageType::Pointer
 274 clitk::SliceBySliceRelativePosition(const MaskImageType * input,
 275                                     const MaskImageType * object,
 276                                     int direction,
 277                                     double threshold,
 278                                     std::string orientation,
 279                                     bool uniqueConnectedComponent,
 280                                     double spacing,
 281                                     bool notflag)
 282 {
 283   typedef clitk::SliceBySliceRelativePositionFilter<MaskImageType> SliceRelPosFilterType;
 284   typename SliceRelPosFilterType::Pointer sliceRelPosFilter = SliceRelPosFilterType::New();
 285   sliceRelPosFilter->VerboseStepOff();
 286   sliceRelPosFilter->WriteStepOff();
 287   sliceRelPosFilter->SetInput(input);
 288   sliceRelPosFilter->SetInputObject(object);
 289   sliceRelPosFilter->SetDirection(direction);
 290   sliceRelPosFilter->SetFuzzyThreshold(threshold);
 291   sliceRelPosFilter->SetOrientationTypeString(orientation);
 292   sliceRelPosFilter->SetResampleBeforeRelativePositionFilter((spacing != -1));
 293   sliceRelPosFilter->SetIntermediateSpacing(spacing);
 294   sliceRelPosFilter->SetUniqueConnectedComponentBySlice(uniqueConnectedComponent);
 295   sliceRelPosFilter->SetNotFlag(notflag);
 296   //  sliceRelPosFilter->SetAutoCropFlag(true); ??
 297   sliceRelPosFilter->Update();
 298   return sliceRelPosFilter->GetOutput();
 299 }
 300 //--------------------------------------------------------------------
 301
 302 //--------------------------------------------------------------------
 303 template<class SliceType>
 304 typename SliceType::PointType
 305 clitk::FindExtremaPointInAGivenDirection(const SliceType * input,
 306                                          typename SliceType::PixelType bg,
 307                                          int direction,
 308                                          bool notFlag,
 309                                          typename SliceType::PointType point,
 310                                          double distanceMax)
 311 {
 312   /*
 313     loop over input pixels, store the index in the fg that is max
 314     according to the given direction.
 315   */
 316   typedef itk::ImageRegionConstIteratorWithIndex<SliceType> IteratorType;
 317   IteratorType iter(input, input->GetLargestPossibleRegion());
 318   iter.GoToBegin();
 319   typename SliceType::IndexType max = input->GetLargestPossibleRegion().GetIndex();
 320   if (notFlag) max = max+input->GetLargestPossibleRegion().GetSize();
 321   while (!iter.IsAtEnd()) {
 322     if (iter.Get() != bg) {
 323       bool test = iter.GetIndex()[direction] >  max[direction];
 324       if (notFlag) test = !test;
 325       if (test) {
 326         typename SliceType::PointType p;
 327         input->TransformIndexToPhysicalPoint(iter.GetIndex(), p);
 328         if ((distanceMax==0) || (p.EuclideanDistanceTo(point) < distanceMax)) {
 329           max = iter.GetIndex();
 330         }
 331       }
 332     }
 333     ++iter;
 334   }
 335   typename SliceType::PointType p;
 336   input->TransformIndexToPhysicalPoint(max, p);
 337   return p;
 338 }
 339 //--------------------------------------------------------------------
 340
 341 //--------------------------------------------------------------------
 342 template<class ImageType>
 343 typename ImageType::Pointer
 344 clitk::CropImageAlongOneAxis(typename ImageType::Pointer image,
 345                              int dim, double min, double max,
 346                              bool autoCrop,
 347                              typename ImageType::PixelType BG)
 348 {
 349   // Compute region size
 350   typename ImageType::RegionType region;
 351   typename ImageType::SizeType size = image->GetLargestPossibleRegion().GetSize();
 352   typename ImageType::PointType p = image->GetOrigin();
 353   p[dim] = min;
 354   typename ImageType::IndexType start;
 355   image->TransformPhysicalPointToIndex(p, start);
 356   p[dim] = max;
 357   typename ImageType::IndexType end;
 358   image->TransformPhysicalPointToIndex(p, end);
 359   size[dim] = fabs(end[dim]-start[dim]);
 360   region.SetIndex(start);
 361   region.SetSize(size);
 362   // Perform Crop
 363   typedef itk::RegionOfInterestImageFilter<ImageType, ImageType> CropFilterType;
 364   typename CropFilterType::Pointer cropFilter = CropFilterType::New();
 365   cropFilter->SetInput(image);
 366   cropFilter->SetRegionOfInterest(region);
 367   cropFilter->Update();
 368   typename ImageType::Pointer result = cropFilter->GetOutput();
 369   // Auto Crop
 370   if (autoCrop) {
 371     result = clitk::AutoCrop<ImageType>(result, BG);
 372   }
 373   return result;
 374 }
 375 //--------------------------------------------------------------------