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