tools/clitkSUVPeakGenericFilter.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 #ifndef CLITKSUVPEAKGENERICFILTER_TXX
  19 #define CLITKSUVPEAKGENERICFILTER_TXX
  20
  21 #include "clitkImageCommon.h"
  22 #include "clitkCropLikeImageFilter.h"
  23 #include "clitkResampleImageWithOptionsFilter.h"
  24
  25 #include <itkConvolutionImageFilter.h>
  26
  27 namespace clitk
  28 {
  29
  30 //--------------------------------------------------------------------
  31 template<class args_info_type>
  32 SUVPeakGenericFilter<args_info_type>::SUVPeakGenericFilter()
  33   :ImageToImageGenericFilter<Self>("SUVPeakGenericFilter")
  34 {
  35   InitializeImageType<2>();
  36   InitializeImageType<3>();
  37   InitializeImageType<4>();
  38 }
  39 //--------------------------------------------------------------------
  40
  41
  42 //--------------------------------------------------------------------
  43 template<class args_info_type>
  44 template<unsigned int Dim>
  45 void SUVPeakGenericFilter<args_info_type>::InitializeImageType()
  46 {
  47   ADD_DEFAULT_IMAGE_TYPES(Dim);
  48 }
  49 //--------------------------------------------------------------------
  50
  51
  52 //--------------------------------------------------------------------
  53 template<class args_info_type>
  54 void SUVPeakGenericFilter<args_info_type>::SetArgsInfo(const args_info_type & a)
  55 {
  56   mArgsInfo=a;
  57
  58   // Set value
  59   this->SetIOVerbose(mArgsInfo.verbose_flag);
  60
  61   if (mArgsInfo.input_given) this->AddInputFilename(mArgsInfo.input_arg);
  62
  63   if (mArgsInfo.mask_given)  this->AddInputFilename(mArgsInfo.mask_arg);
  64   }
  65 //--------------------------------------------------------------------
  66
  67
  68 //--------------------------------------------------------------------
  69 template<class args_info_type>
  70 template<class ImageType>
  71 void SUVPeakGenericFilter<args_info_type>::UpdateWithInputImageType()
  72 {
  73   // Read input
  74   typename ImageType::Pointer input = this->template GetInput<ImageType>(0);
  75
  76   //Read mask
  77   typedef itk::Image<unsigned char, ImageType::ImageDimension> MaskImageType;
  78   typename MaskImageType::Pointer mask;
  79   if(mArgsInfo.mask_given) {
  80       mask = this->template GetInput<MaskImageType>(1);
  81       // Check mask sampling/size
  82       if (!HaveSameSizeAndSpacing<MaskImageType, ImageType>(mask, input)) {
  83         if (mArgsInfo.allow_resize_flag) {
  84           if (mArgsInfo.verbose_flag) {
  85             std::cout << "Resize mask image like input" << std::endl;
  86           }
  87           typedef clitk::ResampleImageWithOptionsFilter<MaskImageType> ResamplerType;
  88           typename ResamplerType::Pointer resampler = ResamplerType::New();
  89           resampler->SetInput(mask); //By default the interpolation in NN, Ok for mask
  90           resampler->SetOutputSpacing(input->GetSpacing());
  91           resampler->SetOutputOrigin(mask->GetOrigin());
  92           resampler->SetGaussianFilteringEnabled(false);
  93           resampler->Update();
  94           mask = resampler->GetOutput();
  95
  96           typedef clitk::CropLikeImageFilter<MaskImageType> FilterType;
  97           typename FilterType::Pointer crop = FilterType::New();
  98           crop->SetInput(mask);
  99           crop->SetCropLikeImage(input);
 100           crop->Update();
 101           mask = crop->GetOutput();
 102
 103         }
 104         else {
 105           std::cerr << "Mask image has a different size/spacing than input. Abort. (Use option --allow_resize)" << std::endl;
 106           exit(-1);
 107         }
 108       }
 109   }
 110   else {
 111       mask = MaskImageType::New();
 112       mask->SetRegions(input->GetLargestPossibleRegion());
 113       mask->SetOrigin(input->GetOrigin());
 114       mask->SetSpacing(input->GetSpacing());
 115       mask->Allocate();
 116       mask->FillBuffer(1);
 117   }
 118
 119   double volume = 1000; //1 cc into mc
 120   if (mArgsInfo.volume_given)
 121     volume *= mArgsInfo.volume_arg;
 122   const double PI = 3.141592653589793238463;
 123   double radius = std::pow(3*volume/(4*PI),1./3);
 124
 125   typename ImageType::Pointer kernel = ComputeMeanFilterKernel<ImageType>(input->GetSpacing(), radius);
 126
 127   // Perform the convolution
 128   typedef itk::ConvolutionImageFilter<ImageType> FilterType;
 129   typename FilterType::Pointer filter = FilterType::New();
 130   filter->SetInput(input);
 131   filter->SetKernelImage(kernel);
 132   filter->Update();
 133   typename ImageType::Pointer output = filter->GetOutput();
 134
 135
 136   typedef itk::ImageRegionConstIteratorWithIndex<ImageType> IteratorType;
 137   typedef itk::ImageRegionConstIteratorWithIndex<MaskImageType> MIteratorType;
 138   IteratorType iters(output, output->GetLargestPossibleRegion());
 139   MIteratorType iterm(mask, mask->GetLargestPossibleRegion());
 140   iters.GoToBegin();
 141   iterm.GoToBegin();
 142   double max = 0.0;
 143   typename ImageType::IndexType index;
 144   while (!iters.IsAtEnd()) {
 145     if (iterm.Get() == 1) { // inside the mask
 146       if (iters.Get() > max) {
 147         max = iters.Get();
 148         index = iters.GetIndex();
 149       }
 150     }
 151     ++iters;
 152     ++iterm;
 153   }
 154   typename ImageType::PointType p;
 155   output->TransformIndexToPhysicalPoint(index, p);
 156   std::cout<<"SUV Peak found in "<< p << " mm with the value " << max << std::endl;
 157 }
 158 //--------------------------------------------------------------------
 159
 160
 161 //--------------------------------------------------------------------
 162 template<class args_info_type>
 163 template<class ImageType>
 164 typename ImageType::Pointer
 165 SUVPeakGenericFilter<args_info_type>::ComputeMeanFilterKernel(const typename ImageType::SpacingType & spacing, double radius)
 166 {
 167   // Some kind of cache to speed up a bit
 168   static std::map<double, typename ImageType::Pointer> cache;
 169   if (cache.find(radius) != cache.end()) {
 170     return cache.find(radius)->second;
 171   }
 172
 173   // Compute a kernel that corresponds to a sphere with 1 inside, 0
 174   // outside and in between proportional to the intersection between
 175   // the pixel and the sphere. Computed by Monte-Carlo because I don't
 176   // know an equation that compute the intersection volume between a
 177   // box and a sphere ...
 178   //auto kernel = ImageType::New();
 179   typename ImageType::Pointer kernel = ImageType::New();
 180
 181   // Size of the kernel in pixel (minimum 3 pixels)
 182   typename ImageType::SizeType size;
 183   size[0] = std::max((int)ceil(radius*2/spacing[0]), 3);
 184   size[1] = std::max((int)ceil(radius*2/spacing[1]), 3);
 185   size[2] = std::max((int)ceil(radius*2/spacing[2]), 3);
 186
 187   // Compute the region, such as the origin is at the center
 188   typename ImageType::IndexType start;
 189   start.Fill(0);
 190   typename ImageType::RegionType region;
 191   region.SetSize(size);
 192   region.SetIndex(start);
 193   kernel->SetRegions(region);
 194   kernel->SetSpacing(spacing);
 195   typename ImageType::PointType origin;
 196   origin[0] = -(double)size[0]/2.0*spacing[0]+spacing[0]/2.0;
 197   origin[1] = -(double)size[1]/2.0*spacing[1]+spacing[1]/2.0;
 198   origin[2] = -(double)size[2]/2.0*spacing[2]+spacing[2]/2.0;
 199   kernel->SetOrigin(origin);
 200   kernel->Allocate();
 201
 202   // Fill the kernel
 203   itk::ImageRegionIteratorWithIndex<ImageType> iter(kernel, region);
 204   typename ImageType::PointType center;
 205   center.Fill(0.0);
 206   typename ImageType::PointType hh; // half a voxel
 207   hh[0] = spacing[0]/2.0;
 208   hh[1] = spacing[1]/2.0;
 209   hh[2] = spacing[2]/2.0;
 210   double h = hh.EuclideanDistanceTo(center); // distance of half a pixel to its center.
 211   std::srand(time(NULL));
 212   double sum = 0.0;
 213   while (!iter.IsAtEnd()) {
 214     typename ImageType::IndexType index = iter.GetIndex();
 215     typename ImageType::PointType p;
 216     kernel->TransformIndexToPhysicalPoint(index, p);
 217     double d = p.EuclideanDistanceTo(center) + h;
 218     if (d<radius) { // inside the sphere
 219       iter.Set(1.0);
 220       sum += 1.0;
 221     }
 222     else { // the box intersect the sphere. We randomly pick point in
 223            // the box and compute the probability to be in/out the
 224            // sphere
 225       int n = 500; // number of samples
 226       double w = 0.0;
 227       //for(auto i=0; i<n; i++) {
 228       for(int i=0; i<n; i++) {
 229         // random position inside the current pixel
 230         typename ImageType::PointType pos;
 231         pos[0] = p[0]+(((double)std::rand()/(double)RAND_MAX)-0.5)*spacing[0];
 232         pos[1] = p[1]+(((double)std::rand()/(double)RAND_MAX)-0.5)*spacing[1];
 233         pos[2] = p[2]+(((double)std::rand()/(double)RAND_MAX)-0.5)*spacing[2];
 234         // distance to center
 235         double distance = pos.EuclideanDistanceTo(center);
 236         // lower/greater than radius
 237         if (distance < radius) w += 1.0;
 238       }
 239       w = w/(double)n;
 240       iter.Set(w);
 241       sum += w;
 242     }
 243     ++iter;
 244   }
 245
 246   // Normalize
 247   iter.GoToBegin();
 248   while (!iter.IsAtEnd()) {
 249     iter.Set(iter.Get()/sum);
 250     ++iter;
 251   }
 252
 253   // Put in cache
 254   cache[radius] = kernel;
 255
 256   return kernel;
 257 }
 258 //--------------------------------------------------------------------
 259
 260 } // end namespace
 261
 262 #endif  //#define CLITKSUVPEAKGENERICFILTER_TXX