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 "clitkCropLikeImageFilter.h"
23 #include "clitkResampleImageWithOptionsFilter.h"
25 #include <itkConvolutionImageFilter.h>
30 //--------------------------------------------------------------------
31 template<class args_info_type>
32 SUVPeakGenericFilter<args_info_type>::SUVPeakGenericFilter()
33 :ImageToImageGenericFilter<Self>("SUVPeakGenericFilter")
35 InitializeImageType<2>();
36 InitializeImageType<3>();
37 InitializeImageType<4>();
39 //--------------------------------------------------------------------
42 //--------------------------------------------------------------------
43 template<class args_info_type>
44 template<unsigned int Dim>
45 void SUVPeakGenericFilter<args_info_type>::InitializeImageType()
47 ADD_DEFAULT_IMAGE_TYPES(Dim);
49 //--------------------------------------------------------------------
52 //--------------------------------------------------------------------
53 template<class args_info_type>
54 void SUVPeakGenericFilter<args_info_type>::SetArgsInfo(const args_info_type & a)
59 this->SetIOVerbose(mArgsInfo.verbose_flag);
61 if (mArgsInfo.input_given) this->AddInputFilename(mArgsInfo.input_arg);
63 if (mArgsInfo.mask_given) this->AddInputFilename(mArgsInfo.mask_arg);
65 //--------------------------------------------------------------------
68 //--------------------------------------------------------------------
69 template<class args_info_type>
70 template<class ImageType>
71 void SUVPeakGenericFilter<args_info_type>::UpdateWithInputImageType()
74 typename ImageType::Pointer input = this->template GetInput<ImageType>(0);
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;
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);
94 mask = resampler->GetOutput();
96 typedef clitk::CropLikeImageFilter<MaskImageType> FilterType;
97 typename FilterType::Pointer crop = FilterType::New();
99 crop->SetCropLikeImage(input);
101 mask = crop->GetOutput();
105 std::cerr << "Mask image has a different size/spacing than input. Abort. (Use option --allow_resize)" << std::endl;
111 mask = MaskImageType::New();
112 mask->SetRegions(input->GetLargestPossibleRegion());
113 mask->SetOrigin(input->GetOrigin());
114 mask->SetSpacing(input->GetSpacing());
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);
125 typename ImageType::Pointer kernel = ComputeMeanFilterKernel<ImageType>(input->GetSpacing(), radius);
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);
133 typename ImageType::Pointer output = filter->GetOutput();
136 typedef itk::ImageRegionConstIteratorWithIndex<ImageType> IteratorType;
137 typedef itk::ImageRegionConstIteratorWithIndex<MaskImageType> MIteratorType;
138 IteratorType iters(output, output->GetLargestPossibleRegion());
139 MIteratorType iterm(mask, mask->GetLargestPossibleRegion());
143 typename ImageType::IndexType index;
144 while (!iters.IsAtEnd()) {
145 if (iterm.Get() == 1) { // inside the mask
146 if (iters.Get() > max) {
148 index = iters.GetIndex();
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;
158 //--------------------------------------------------------------------
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)
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;
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();
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);
187 // Compute the region, such as the origin is at the center
188 typename ImageType::IndexType start;
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);
203 itk::ImageRegionIteratorWithIndex<ImageType> iter(kernel, region);
204 typename ImageType::PointType center;
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));
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
222 else { // the box intersect the sphere. We randomly pick point in
223 // the box and compute the probability to be in/out the
225 int n = 500; // number of samples
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;
248 while (!iter.IsAtEnd()) {
249 iter.Set(iter.Get()/sum);
254 cache[radius] = kernel;
258 //--------------------------------------------------------------------
262 #endif //#define CLITKSUVPEAKGENERICFILTER_TXX