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   const double PI = 3.141592653589793238463;
 121   double radius = std::pow(3*volume/(4*PI),1./3);
 122
 123   typename ImageType::Pointer kernel = ComputeMeanFilterKernel<ImageType>(input->GetSpacing(), radius);
 124
 125   // Perform the convolution
 126   typedef itk::ConvolutionImageFilter<ImageType> FilterType;
 127   typename FilterType::Pointer filter = FilterType::New();
 128   filter->SetInput(input);
 129   filter->SetKernelImage(kernel);
 130   filter->Update();
 131   typename ImageType::Pointer output = filter->GetOutput();
 132
 133
 134   typedef itk::ImageRegionConstIteratorWithIndex<ImageType> IteratorType;
 135   typedef itk::ImageRegionConstIteratorWithIndex<MaskImageType> MIteratorType;
 136   IteratorType iters(output, output->GetLargestPossibleRegion());
 137   MIteratorType iterm(mask, mask->GetLargestPossibleRegion());
 138   iters.GoToBegin();
 139   iterm.GoToBegin();
 140   double max = 0.0;
 141   typename ImageType::IndexType index;
 142   while (!iters.IsAtEnd()) {
 143     if (iterm.Get() == 1) { // inside the mask
 144       if (iters.Get() > max) {
 145         max = iters.Get();
 146         index = iters.GetIndex();
 147       }
 148     }
 149     ++iters;
 150     ++iterm;
 151   }
 152   typename ImageType::PointType p;
 153   output->TransformIndexToPhysicalPoint(index, p);
 154   std::cout<<"SUV Peak found in "<< p << " mm with the value " << max << std::endl;
 155 }
 156 //--------------------------------------------------------------------
 157
 158
 159 //--------------------------------------------------------------------
 160 template<class args_info_type>
 161 template<class ImageType>
 162 typename ImageType::Pointer
 163 SUVPeakGenericFilter<args_info_type>::ComputeMeanFilterKernel(const typename ImageType::SpacingType & spacing, double radius)
 164 {
 165   // Some kind of cache to speed up a bit
 166   static std::map<double, typename ImageType::Pointer> cache;
 167   if (cache.find(radius) != cache.end()) {
 168     return cache.find(radius)->second;
 169   }
 170
 171   // Compute a kernel that corresponds to a sphere with 1 inside, 0
 172   // outside and in between proportional to the intersection between
 173   // the pixel and the sphere. Computed by Monte-Carlo because I don't
 174   // know an equation that compute the intersection volume between a
 175   // box and a sphere ...
 176   //auto kernel = ImageType::New();
 177   typename ImageType::Pointer kernel = ImageType::New();
 178
 179   // Size of the kernel in pixel (minimum 3 pixels)
 180   typename ImageType::SizeType size;
 181   size[0] = std::max((int)ceil(radius*2/spacing[0]), 3);
 182   size[1] = std::max((int)ceil(radius*2/spacing[1]), 3);
 183   size[2] = std::max((int)ceil(radius*2/spacing[2]), 3);
 184
 185   // Compute the region, such as the origin is at the center
 186   typename ImageType::IndexType start;
 187   start.Fill(0);
 188   typename ImageType::RegionType region;
 189   region.SetSize(size);
 190   region.SetIndex(start);
 191   kernel->SetRegions(region);
 192   kernel->SetSpacing(spacing);
 193   typename ImageType::PointType origin;
 194   origin[0] = -(double)size[0]/2.0*spacing[0]+spacing[0]/2.0;
 195   origin[1] = -(double)size[1]/2.0*spacing[1]+spacing[1]/2.0;
 196   origin[2] = -(double)size[2]/2.0*spacing[2]+spacing[2]/2.0;
 197   kernel->SetOrigin(origin);
 198   kernel->Allocate();
 199
 200   // Fill the kernel
 201   itk::ImageRegionIteratorWithIndex<ImageType> iter(kernel, region);
 202   typename ImageType::PointType center;
 203   center.Fill(0.0);
 204   typename ImageType::PointType hh; // half a voxel
 205   hh[0] = spacing[0]/2.0;
 206   hh[1] = spacing[1]/2.0;
 207   hh[2] = spacing[2]/2.0;
 208   double h = hh.EuclideanDistanceTo(center); // distance of half a pixel to its center.
 209   std::srand(time(NULL));
 210   double sum = 0.0;
 211   while (!iter.IsAtEnd()) {
 212     typename ImageType::IndexType index = iter.GetIndex();
 213     typename ImageType::PointType p;
 214     kernel->TransformIndexToPhysicalPoint(index, p);
 215     double d = p.EuclideanDistanceTo(center) + h;
 216     if (d<radius) { // inside the sphere
 217       iter.Set(1.0);
 218       sum += 1.0;
 219     }
 220     else { // the box intersect the sphere. We randomly pick point in
 221            // the box and compute the probability to be in/out the
 222            // sphere
 223       int n = 500; // number of samples
 224       double w = 0.0;
 225       //for(auto i=0; i<n; i++) {
 226       for(int i=0; i<n; i++) {
 227         // random position inside the current pixel
 228         typename ImageType::PointType pos;
 229         pos[0] = p[0]+(((double)std::rand()/(double)RAND_MAX)-0.5)*spacing[0];
 230         pos[1] = p[1]+(((double)std::rand()/(double)RAND_MAX)-0.5)*spacing[1];
 231         pos[2] = p[2]+(((double)std::rand()/(double)RAND_MAX)-0.5)*spacing[2];
 232         // distance to center
 233         double distance = pos.EuclideanDistanceTo(center);
 234         // lower/greater than radius
 235         if (distance < radius) w += 1.0;
 236       }
 237       w = w/(double)n;
 238       iter.Set(w);
 239       sum += w;
 240     }
 241     ++iter;
 242   }
 243
 244   // Normalize
 245   iter.GoToBegin();
 246   while (!iter.IsAtEnd()) {
 247     iter.Set(iter.Get()/sum);
 248     ++iter;
 249   }
 250
 251   // Put in cache
 252   cache[radius] = kernel;
 253
 254   return kernel;
 255 }
 256 //--------------------------------------------------------------------
 257
 258 } // end namespace
 259
 260 #endif  //#define CLITKSUVPEAKGENERICFILTER_TXX