1 /*=========================================================================
2 Program: vv http://www.creatis.insa-lyon.fr/rio/vv
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
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.
13 It is distributed under dual licence
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
21 #include "clitkImageCommon.h"
22 #include <itkConvolutionImageFilter.h>
27 //--------------------------------------------------------------------
28 template<class args_info_type>
29 SUVPeakGenericFilter<args_info_type>::SUVPeakGenericFilter()
30 :ImageToImageGenericFilter<Self>("SUVPeakGenericFilter")
32 InitializeImageType<2>();
33 InitializeImageType<3>();
34 InitializeImageType<4>();
36 //--------------------------------------------------------------------
39 //--------------------------------------------------------------------
40 template<class args_info_type>
41 template<unsigned int Dim>
42 void SUVPeakGenericFilter<args_info_type>::InitializeImageType()
44 ADD_DEFAULT_IMAGE_TYPES(Dim);
46 //--------------------------------------------------------------------
49 //--------------------------------------------------------------------
50 template<class args_info_type>
51 void SUVPeakGenericFilter<args_info_type>::SetArgsInfo(const args_info_type & a)
56 this->SetIOVerbose(mArgsInfo.verbose_flag);
58 if (mArgsInfo.input_given) this->AddInputFilename(mArgsInfo.input_arg);
60 if (mArgsInfo.mask_given) this->AddInputFilename(mArgsInfo.mask_arg);
62 //--------------------------------------------------------------------
65 //--------------------------------------------------------------------
66 template<class args_info_type>
67 template<class ImageType>
68 void SUVPeakGenericFilter<args_info_type>::UpdateWithInputImageType()
71 typename ImageType::Pointer input = this->template GetInput<ImageType>(0);
74 typedef itk::Image<unsigned char, ImageType::ImageDimension> MaskImageType;
75 typename MaskImageType::Pointer mask = NULL;
76 if(mArgsInfo.mask_given) {
77 mask = this->template GetInput<MaskImageType>(1);
80 mask = MaskImageType::New();
81 mask->SetRegions(input->GetLargestPossibleRegion());
82 mask->SetOrigin(input->GetOrigin());
83 mask->SetSpacing(input->GetSpacing());
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);
92 typename ImageType::Pointer kernel = ComputeMeanFilterKernel<ImageType>(input->GetSpacing(), radius);
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);
100 typename ImageType::Pointer output = filter->GetOutput();
103 typedef itk::ImageRegionConstIteratorWithIndex<ImageType> IteratorType;
104 typedef itk::ImageRegionConstIteratorWithIndex<MaskImageType> MIteratorType;
105 IteratorType iters(output, output->GetLargestPossibleRegion());
106 MIteratorType iterm(mask, mask->GetLargestPossibleRegion());
110 typename ImageType::IndexType index;
111 while (!iters.IsAtEnd()) {
112 if (iterm.Get() == 1) { // inside the mask
113 if (iters.Get() > max) {
115 index = iters.GetIndex();
121 typename ImageType::PointType p;
122 output->TransformIndexToPhysicalPoint(index, p);
123 std::cout<<"SUV Peak found in "<< p << " with the value " << max << std::endl;
125 //--------------------------------------------------------------------
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)
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;
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();
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);
154 // Compute the region, such as the origin is at the center
155 typename ImageType::IndexType start;
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);
170 itk::ImageRegionIteratorWithIndex<ImageType> iter(kernel, region);
171 typename ImageType::PointType center;
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));
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
189 else { // the box intersect the sphere. We randomly pick point in
190 // the box and compute the probability to be in/out the
192 int n = 500; // number of samples
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;
215 while (!iter.IsAtEnd()) {
216 iter.Set(iter.Get()/sum);
221 cache[radius] = kernel;
225 //--------------------------------------------------------------------
229 #endif //#define CLITKSUVPEAKGENERICFILTER_TXX