target_link_libraries(clitkImageStatistics clitkCommon)
set(TOOLS_INSTALL ${TOOLS_INSTALL} clitkImageStatistics)
+ WRAP_GGO(clitkSUVPeak_GGO_C clitkSUVPeak.ggo)
+ add_executable(clitkSUVPeak clitkSUVPeak.cxx clitkSUVPeakGenericFilter.cxx ${clitkSUVPeak_GGO_C})
+ target_link_libraries(clitkSUVPeak clitkCommon)
+ set(TOOLS_INSTALL ${TOOLS_INSTALL} clitkSUVPeak)
+
+ WRAP_GGO(clitkRTStructStatistics_GGO_C clitkRTStructStatistics.ggo)
+ add_executable(clitkRTStructStatistics clitkRTStructStatistics.cxx clitkRTStructStatisticsGenericFilter.cxx ${clitkRTStructStatistics_GGO_C})
+ target_link_libraries(clitkRTStructStatistics clitkCommon)
+ set(TOOLS_INSTALL ${TOOLS_INSTALL} clitkRTStructStatistics)
+
WRAP_GGO(clitkVFConvert_GGO_C clitkVFConvert.ggo)
add_executable(clitkVFConvert clitkVFConvert.cxx clitkVFConvertGenericFilter.cxx ${clitkVFConvert_GGO_C})
target_link_libraries(clitkVFConvert clitkCommon )
--- /dev/null
+/*=========================================================================SUVP
+ Program: vv http://www.creatis.insa-lyon.fr/rio/vv
+
+ Authors belong to:
+ - University of LYON http://www.universite-lyon.fr/
+ - Léon Bérard cancer center http://www.centreleonberard.fr
+ - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr
+
+ This software is distributed WITHOUT ANY WARRANTY; without even
+ the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+ PURPOSE. See the copyright notices for more information.
+
+ It is distributed under dual licence
+
+ - BSD See included LICENSE.txt file
+ - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
+===========================================================================**/
+#ifndef CLITKRTSTRUCTSTATISTICS_CXX
+#define CLITKRTSTRUCTSTATISTICS_CXX
+
+// clitk include
+#include "clitkRTStructStatistics_ggo.h"
+#include "clitkRTStructStatisticsGenericFilter.h"
+#include "clitkIO.h"
+
+//--------------------------------------------------------------------
+int main(int argc, char * argv[])
+{
+
+ // Init command line
+ GGO(clitkRTStructStatistics, args_info);
+ CLITK_INIT;
+
+ // Creation of a generic filter
+ typedef clitk::RTStructStatisticsGenericFilter<args_info_clitkRTStructStatistics> FilterType;
+ FilterType::Pointer filter = FilterType::New();
+
+ // Go !
+ filter->SetArgsInfo(args_info);
+ CLITK_TRY_CATCH_EXIT(filter->Update());
+
+ // this is the end my friend
+ return EXIT_SUCCESS;
+} // end main
+
+#endif //define CLITKRTSTRUCTSTATISTICS_CXX
--- /dev/null
+#File clitkRTStructStatistics.ggo
+package "clitkRTStructStatistics"
+version "2.0"
+purpose "Find the centroid/roundness of a binarized image."
+
+option "config" - "Config file" string no
+option "verbose" v "Verbose" flag off
+
+option "input" i "Input image filename (mask)" string yes
--- /dev/null
+/*=========================================================================
+ Program: vv http://www.creatis.insa-lyon.fr/rio/vv
+
+ Authors belong to:
+ - University of LYON http://www.universite-lyon.fr/
+ - Léon Bérard cancer center http://www.centreleonberard.fr
+ - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr
+
+ This software is distributed WITHOUT ANY WARRANTY; without even
+ the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+ PURPOSE. See the copyright notices for more information.
+
+ It is distributed under dual licence
+
+ - BSD See included LICENSE.txt file
+ - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
+===========================================================================**/
+#ifndef CLITKRTSTRUCTSTATISTICSGENERICFILTER_CXX
+#define CLITKRTSTRUCTSTATISTICSGENERICFILTER_CXX
+
+#include "clitkRTStructStatisticsGenericFilter.h"
+
+namespace clitk {
+ // Specialisation
+// template<>
+// class RTStructStatisticsGenericFilter<args_info_clitkC>;
+
+
+}
+
+#endif //define CLITKRTSTRUCTSTATISTICSGENERICFILTER_CXX
--- /dev/null
+/*=========================================================================
+ Program: vv http://www.creatis.insa-lyon.fr/rio/vv
+
+ Authors belong to:
+ - University of LYON http://www.universite-lyon.fr/
+ - Léon Bérard cancer center http://www.centreleonberard.fr
+ - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr
+
+ This software is distributed WITHOUT ANY WARRANTY; without even
+ the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+ PURPOSE. See the copyright notices for more information.
+
+ It is distributed under dual licence
+
+ - BSD See included LICENSE.txt file
+ - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
+===========================================================================**/
+#ifndef CLITKRTSTRUCTSTATISTICSGENERICFILTER_H
+#define CLITKRTSTRUCTSTATISTICSGENERICFILTER_H
+/**
+ -------------------------------------------------------------------
+ * @file clitkRTStructStatisticsGenericFilter.h
+ * @author Thomas Baudier <thomas.baudier@creatis.insa-lyon.fr>
+ * @date 11 Jul 2016 08:37:53
+
+ * @brief
+ -------------------------------------------------------------------*/
+
+// clitk include
+#include "clitkCommon.h"
+#include "clitkImageToImageGenericFilter.h"
+#include "clitkRTStructStatistics_ggo.h"
+
+// itk include
+#include "itkImage.h"
+#include "itkImageIOBase.h"
+
+//--------------------------------------------------------------------
+namespace clitk {
+
+ template<class args_info_type>
+ class ITK_EXPORT RTStructStatisticsGenericFilter:
+ public clitk::ImageToImageGenericFilter<RTStructStatisticsGenericFilter<args_info_type> > {
+
+ public:
+
+ // Constructor
+ RTStructStatisticsGenericFilter ();
+
+ // Types
+ typedef RTStructStatisticsGenericFilter Self;
+ typedef ImageToImageGenericFilterBase Superclass;
+ typedef itk::SmartPointer<Self> Pointer;
+ typedef itk::SmartPointer<const Self> ConstPointer;
+
+ // New
+ itkNewMacro(Self);
+
+
+
+ //--------------------------------------------------------------------
+ void SetArgsInfo(const args_info_type & a);
+
+ //--------------------------------------------------------------------
+ // Main function called each time the filter is updated
+ template<class InputImageType>
+ void UpdateWithInputImageType();
+
+ protected:
+ template<unsigned int Dim> void InitializeImageType();
+ args_info_type mArgsInfo;
+
+ }; // end class RTStructStatisticsGenericFilter
+
+} // end namespace
+//--------------------------------------------------------------------
+
+
+#ifndef ITK_MANUAL_INSTANTIATION
+#include "clitkRTStructStatisticsGenericFilter.txx"
+#endif
+
+#endif //#define CLITKRTSTRUCTSTATISTICSGENERICFILTER_H
+
--- /dev/null
+/*=========================================================================
+ Program: vv http://www.creatis.insa-lyon.fr/rio/vv
+
+ Authors belong to:
+ - University of LYON http://www.universite-lyon.fr/
+ - Léon Bérard cancer center http://www.centreleonberard.fr
+ - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr
+
+ This software is distributed WITHOUT ANY WARRANTY; without even
+ the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+ PURPOSE. See the copyright notices for more information.
+
+ It is distributed under dual licence
+
+ - BSD See included LICENSE.txt file
+ - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
+===========================================================================**/
+#ifndef CLITKRTSTRUCTSTATISTICSGENERICFILTER_TXX
+#define CLITKRTSTRUCTSTATISTICSGENERICFILTER_TXX
+
+#include "clitkImageCommon.h"
+
+#include "itkConnectedComponentImageFilter.h"
+#include "itkLabelImageToShapeLabelMapFilter.h"
+
+namespace clitk
+{
+
+//--------------------------------------------------------------------
+template<class args_info_type>
+RTStructStatisticsGenericFilter<args_info_type>::RTStructStatisticsGenericFilter()
+ :ImageToImageGenericFilter<Self>("RTStructStatisticsGenericFilter")
+{
+ InitializeImageType<2>();
+ InitializeImageType<3>();
+ InitializeImageType<4>();
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template<class args_info_type>
+template<unsigned int Dim>
+void RTStructStatisticsGenericFilter<args_info_type>::InitializeImageType()
+{
+ ADD_DEFAULT_IMAGE_TYPES(Dim);
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template<class args_info_type>
+void RTStructStatisticsGenericFilter<args_info_type>::SetArgsInfo(const args_info_type & a)
+{
+ mArgsInfo=a;
+
+ // Set value
+ this->SetIOVerbose(mArgsInfo.verbose_flag);
+
+ if (mArgsInfo.input_given) this->AddInputFilename(mArgsInfo.input_arg);
+
+ }
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template<class args_info_type>
+template<class ImageType>
+void RTStructStatisticsGenericFilter<args_info_type>::UpdateWithInputImageType()
+{
+ // Read mask input
+ typedef itk::Image<unsigned char, ImageType::ImageDimension> MaskInputImageType;
+ typename MaskInputImageType::Pointer mask = NULL;
+ mask = this->template GetInput<MaskInputImageType>(0);
+
+ //Create the Shape Label Map from the mask
+ typedef itk::Image< unsigned char, ImageType::ImageDimension > OutputImageType;
+ typedef itk::ShapeLabelObject< unsigned char, ImageType::ImageDimension > ShapeLabelObjectType;
+ typedef itk::LabelMap< ShapeLabelObjectType > LabelMapType;
+ typedef itk::ConnectedComponentImageFilter <MaskInputImageType, OutputImageType > ConnectedComponentImageFilterType;
+ typedef itk::LabelImageToShapeLabelMapFilter< OutputImageType, LabelMapType> I2LType;
+
+ typename ConnectedComponentImageFilterType::Pointer connected = ConnectedComponentImageFilterType::New ();
+ connected->SetInput(mask);
+ connected->FullyConnectedOn();
+ connected->Update();
+
+ //Create a map to contain all connectedComponent (even a little pixel)
+ typename I2LType::Pointer i2l = I2LType::New();
+ i2l->SetInput( connected->GetOutput() );
+ i2l->SetComputePerimeter(true);
+ i2l->Update();
+
+ // Retrieve the biggest component
+ LabelMapType *labelMap = i2l->GetOutput();
+ int largestComponent(0);
+ int nbPixel(0);
+ for (unsigned int n = 0; n < labelMap->GetNumberOfLabelObjects(); ++n)
+ {
+ ShapeLabelObjectType *labelObject = labelMap->GetNthLabelObject(n);
+ if (labelObject->GetNumberOfPixels() > nbPixel)
+ {
+ nbPixel = labelObject->GetNumberOfPixels();
+ largestComponent = n;
+ }
+ }
+
+ //Write statitistics on the largest component
+ ShapeLabelObjectType *labelObject = labelMap->GetNthLabelObject(largestComponent);
+ std::cout << " Centroid: " << std::endl;
+ std::cout << labelObject->GetCentroid()[0] << std::endl;
+ std::cout << labelObject->GetCentroid()[1] << std::endl;
+ std::cout << labelObject->GetCentroid()[2] << std::endl;
+ std::cout << " Roundness: " << std::endl;
+ std::cout << labelObject->GetRoundness() << std::endl;
+
+}
+//--------------------------------------------------------------------
+
+
+
+} // end namespace
+
+#endif //#define CLITKRTSTRUCTSTATISTICSGENERICFILTER_TXX
--- /dev/null
+/*=========================================================================
+ Program: vv http://www.creatis.insa-lyon.fr/rio/vv
+
+ Authors belong to:
+ - University of LYON http://www.universite-lyon.fr/
+ - Léon Bérard cancer center http://www.centreleonberard.fr
+ - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr
+
+ This software is distributed WITHOUT ANY WARRANTY; without even
+ the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+ PURPOSE. See the copyright notices for more information.
+
+ It is distributed under dual licence
+
+ - BSD See included LICENSE.txt file
+ - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
+===========================================================================**/
+#ifndef CLITKSUVPEAK_CXX
+#define CLITKSUVPEAK_CXX
+
+// clitk include
+#include "clitkSUVPeak_ggo.h"
+#include "clitkSUVPeakGenericFilter.h"
+#include "clitkIO.h"
+
+//--------------------------------------------------------------------
+int main(int argc, char * argv[])
+{
+
+ // Init command line
+ GGO(clitkSUVPeak, args_info);
+ CLITK_INIT;
+
+ // Creation of a generic filter
+ typedef clitk::SUVPeakGenericFilter<args_info_clitkSUVPeak> FilterType;
+ FilterType::Pointer filter = FilterType::New();
+
+ // Go !
+ filter->SetArgsInfo(args_info);
+ CLITK_TRY_CATCH_EXIT(filter->Update());
+
+ // this is the end my friend
+ return EXIT_SUCCESS;
+} // end main
+
+#endif //define CLITKSUVPEAK_CXX
--- /dev/null
+#File clitkSUVPeak.ggo
+package "clitkSUVPeak"
+version "2.0"
+#This tool supports multiple images on the input, or even 4D, but all images must be of the same type and dimensions.
+purpose "Compute statistics on an image, or on part of an image specified by a mask and label(s). The tool also supports multichannel images, which is useful, e.g., for vector fields. All channels are processed (separately) by default, but only one channel may be chosen."
+
+option "config" - "Config file" string no
+option "verbose" v "Verbose" flag off
+
+option "input" i "Input first image filename" string yes
+option "mask" m "Mask image filename (uchar)" string no
--- /dev/null
+/*=========================================================================
+ Program: vv http://www.creatis.insa-lyon.fr/rio/vv
+
+ Authors belong to:
+ - University of LYON http://www.universite-lyon.fr/
+ - Léon Bérard cancer center http://www.centreleonberard.fr
+ - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr
+
+ This software is distributed WITHOUT ANY WARRANTY; without even
+ the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+ PURPOSE. See the copyright notices for more information.
+
+ It is distributed under dual licence
+
+ - BSD See included LICENSE.txt file
+ - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
+===========================================================================**/
+#ifndef CLITKSUVPEAKGENERICFILTER_CXX
+#define CLITKSUVPEAKGENERICFILTER_CXX
+
+#include "clitkSUVPeakGenericFilter.h"
+
+namespace clitk {
+ // Specialisation
+// template<>
+// class SUVPeakGenericFilter<args_info_clitkSUVPeak>;
+
+
+}
+
+#endif //define CLITKSUVPEAKGENERICFILTER_CXX
--- /dev/null
+/*=========================================================================
+ Program: vv http://www.creatis.insa-lyon.fr/rio/vv
+
+ Authors belong to:
+ - University of LYON http://www.universite-lyon.fr/
+ - Léon Bérard cancer center http://www.centreleonberard.fr
+ - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr
+
+ This software is distributed WITHOUT ANY WARRANTY; without even
+ the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+ PURPOSE. See the copyright notices for more information.
+
+ It is distributed under dual licence
+
+ - BSD See included LICENSE.txt file
+ - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
+===========================================================================**/
+#ifndef CLITKSUVPEAKGENERICFILTER_H
+#define CLITKSUVPEAKGENERICFILTER_H
+/**
+ -------------------------------------------------------------------
+ * @file clitkSUVPeakGenericFilter.h
+ * @author David Sarrut <David.Sarrut@creatis.insa-lyon.fr>
+ * @date 23 Feb 2008 08:37:53
+
+ * @brief
+ -------------------------------------------------------------------*/
+
+// clitk include
+#include "clitkCommon.h"
+#include "clitkImageToImageGenericFilter.h"
+#include "clitkSUVPeak_ggo.h"
+
+// itk include
+#include "itkImage.h"
+#include "itkImageIOBase.h"
+#include "itkImageRegionIterator.h"
+#include "itkImageRegionConstIterator.h"
+
+//--------------------------------------------------------------------
+namespace clitk {
+
+ template<class args_info_type>
+ class ITK_EXPORT SUVPeakGenericFilter:
+ public clitk::ImageToImageGenericFilter<SUVPeakGenericFilter<args_info_type> > {
+
+ public:
+
+ // Constructor
+ SUVPeakGenericFilter ();
+
+ // Types
+ typedef SUVPeakGenericFilter Self;
+ typedef ImageToImageGenericFilterBase Superclass;
+ typedef itk::SmartPointer<Self> Pointer;
+ typedef itk::SmartPointer<const Self> ConstPointer;
+
+ // New
+ itkNewMacro(Self);
+
+
+
+ //--------------------------------------------------------------------
+ void SetArgsInfo(const args_info_type & a);
+
+ //SUVPeak
+ template<class ImageType> typename ImageType::Pointer ComputeMeanFilterKernel(const typename ImageType::SpacingType & spacing, double radius);
+
+ //--------------------------------------------------------------------
+ // Main function called each time the filter is updated
+ template<class InputImageType>
+ void UpdateWithInputImageType();
+
+ protected:
+ template<unsigned int Dim> void InitializeImageType();
+ args_info_type mArgsInfo;
+
+ }; // end class SUVPeakGenericFilter
+
+} // end namespace
+//--------------------------------------------------------------------
+
+
+#ifndef ITK_MANUAL_INSTANTIATION
+#include "clitkSUVPeakGenericFilter.txx"
+#endif
+
+#endif //#define CLITKSUVPEAKGENERICFILTER_H
+
--- /dev/null
+/*=========================================================================
+ Program: vv http://www.creatis.insa-lyon.fr/rio/vv
+
+ Authors belong to:
+ - University of LYON http://www.universite-lyon.fr/
+ - Léon Bérard cancer center http://www.centreleonberard.fr
+ - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr
+
+ This software is distributed WITHOUT ANY WARRANTY; without even
+ the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+ PURPOSE. See the copyright notices for more information.
+
+ It is distributed under dual licence
+
+ - BSD See included LICENSE.txt file
+ - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
+===========================================================================**/
+#ifndef CLITKSUVPEAKGENERICFILTER_TXX
+#define CLITKSUVPEAKGENERICFILTER_TXX
+
+#include "clitkImageCommon.h"
+#include <itkConvolutionImageFilter.h>
+
+namespace clitk
+{
+
+//--------------------------------------------------------------------
+template<class args_info_type>
+SUVPeakGenericFilter<args_info_type>::SUVPeakGenericFilter()
+ :ImageToImageGenericFilter<Self>("SUVPeakGenericFilter")
+{
+ InitializeImageType<2>();
+ InitializeImageType<3>();
+ InitializeImageType<4>();
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template<class args_info_type>
+template<unsigned int Dim>
+void SUVPeakGenericFilter<args_info_type>::InitializeImageType()
+{
+ ADD_DEFAULT_IMAGE_TYPES(Dim);
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template<class args_info_type>
+void SUVPeakGenericFilter<args_info_type>::SetArgsInfo(const args_info_type & a)
+{
+ mArgsInfo=a;
+
+ // Set value
+ this->SetIOVerbose(mArgsInfo.verbose_flag);
+
+ if (mArgsInfo.input_given) this->AddInputFilename(mArgsInfo.input_arg);
+
+ if (mArgsInfo.mask_given) this->AddInputFilename(mArgsInfo.mask_arg);
+ }
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template<class args_info_type>
+template<class ImageType>
+void SUVPeakGenericFilter<args_info_type>::UpdateWithInputImageType()
+{
+ // Read input
+ typename ImageType::Pointer input = this->template GetInput<ImageType>(0);
+
+ //Read mask
+ typedef itk::Image<unsigned char, ImageType::ImageDimension> MaskImageType;
+ typename MaskImageType::Pointer mask = NULL;
+ if(mArgsInfo.mask_given) {
+ mask = this->template GetInput<MaskImageType>(1);
+ }
+ else {
+ mask = MaskImageType::New();
+ mask->SetRegions(input->GetLargestPossibleRegion());
+ mask->SetOrigin(input->GetOrigin());
+ mask->SetSpacing(input->GetSpacing());
+ mask->Allocate();
+ mask->FillBuffer(1);
+ }
+
+ double volume = 1000; //1 cc into mc
+ const double PI = 3.141592653589793238463;
+ double radius = std::pow(3*volume/(4*PI),1./3);
+
+ typename ImageType::Pointer kernel = ComputeMeanFilterKernel<ImageType>(input->GetSpacing(), radius);
+
+ // Perform the convolution
+ typedef itk::ConvolutionImageFilter<ImageType> FilterType;
+ typename FilterType::Pointer filter = FilterType::New();
+ filter->SetInput(input);
+ filter->SetKernelImage(kernel);
+ filter->Update();
+ typename ImageType::Pointer output = filter->GetOutput();
+
+
+ typedef itk::ImageRegionConstIteratorWithIndex<ImageType> IteratorType;
+ typedef itk::ImageRegionConstIteratorWithIndex<MaskImageType> MIteratorType;
+ IteratorType iters(output, output->GetLargestPossibleRegion());
+ MIteratorType iterm(mask, mask->GetLargestPossibleRegion());
+ iters.GoToBegin();
+ iterm.GoToBegin();
+ double max = 0.0;
+ typename ImageType::IndexType index;
+ while (!iters.IsAtEnd()) {
+ if (iterm.Get() == 1) { // inside the mask
+ if (iters.Get() > max) {
+ max = iters.Get();
+ index = iters.GetIndex();
+ }
+ }
+ ++iters;
+ ++iterm;
+ }
+ typename ImageType::PointType p;
+ output->TransformIndexToPhysicalPoint(index, p);
+ std::cout<<"SUV Peak found in "<< p << " with the value " << max << std::endl;
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template<class args_info_type>
+template<class ImageType>
+typename ImageType::Pointer
+SUVPeakGenericFilter<args_info_type>::ComputeMeanFilterKernel(const typename ImageType::SpacingType & spacing, double radius)
+{
+ // Some kind of cache to speed up a bit
+ static std::map<double, typename ImageType::Pointer> cache;
+ if (cache.find(radius) != cache.end()) {
+ return cache.find(radius)->second;
+ }
+
+ // Compute a kernel that corresponds to a sphere with 1 inside, 0
+ // outside and in between proportional to the intersection between
+ // the pixel and the sphere. Computed by Monte-Carlo because I don't
+ // know an equation that compute the intersection volume between a
+ // box and a sphere ...
+ //auto kernel = ImageType::New();
+ typename ImageType::Pointer kernel = ImageType::New();
+
+ // Size of the kernel in pixel (minimum 3 pixels)
+ typename ImageType::SizeType size;
+ size[0] = std::max((int)ceil(radius*2/spacing[0]), 3);
+ size[1] = std::max((int)ceil(radius*2/spacing[1]), 3);
+ size[2] = std::max((int)ceil(radius*2/spacing[2]), 3);
+
+ // Compute the region, such as the origin is at the center
+ typename ImageType::IndexType start;
+ start.Fill(0);
+ typename ImageType::RegionType region;
+ region.SetSize(size);
+ region.SetIndex(start);
+ kernel->SetRegions(region);
+ kernel->SetSpacing(spacing);
+ typename ImageType::PointType origin;
+ origin[0] = -(double)size[0]/2.0*spacing[0]+spacing[0]/2.0;
+ origin[1] = -(double)size[1]/2.0*spacing[1]+spacing[1]/2.0;
+ origin[2] = -(double)size[2]/2.0*spacing[2]+spacing[2]/2.0;
+ kernel->SetOrigin(origin);
+ kernel->Allocate();
+
+ // Fill the kernel
+ itk::ImageRegionIteratorWithIndex<ImageType> iter(kernel, region);
+ typename ImageType::PointType center;
+ center.Fill(0.0);
+ typename ImageType::PointType hh; // half a voxel
+ hh[0] = spacing[0]/2.0;
+ hh[1] = spacing[1]/2.0;
+ hh[2] = spacing[2]/2.0;
+ double h = hh.EuclideanDistanceTo(center); // distance of half a pixel to its center.
+ std::srand(time(NULL));
+ double sum = 0.0;
+ while (!iter.IsAtEnd()) {
+ typename ImageType::IndexType index = iter.GetIndex();
+ typename ImageType::PointType p;
+ kernel->TransformIndexToPhysicalPoint(index, p);
+ double d = p.EuclideanDistanceTo(center) + h;
+ if (d<radius) { // inside the sphere
+ iter.Set(1.0);
+ sum += 1.0;
+ }
+ else { // the box intersect the sphere. We randomly pick point in
+ // the box and compute the probability to be in/out the
+ // sphere
+ int n = 500; // number of samples
+ double w = 0.0;
+ //for(auto i=0; i<n; i++) {
+ for(int i=0; i<n; i++) {
+ // random position inside the current pixel
+ typename ImageType::PointType pos;
+ pos[0] = p[0]+(((double)std::rand()/(double)RAND_MAX)-0.5)*spacing[0];
+ pos[1] = p[1]+(((double)std::rand()/(double)RAND_MAX)-0.5)*spacing[1];
+ pos[2] = p[2]+(((double)std::rand()/(double)RAND_MAX)-0.5)*spacing[2];
+ // distance to center
+ double distance = pos.EuclideanDistanceTo(center);
+ // lower/greater than radius
+ if (distance < radius) w += 1.0;
+ }
+ w = w/(double)n;
+ iter.Set(w);
+ sum += w;
+ }
+ ++iter;
+ }
+
+ // Normalize
+ iter.GoToBegin();
+ while (!iter.IsAtEnd()) {
+ iter.Set(iter.Get()/sum);
+ ++iter;
+ }
+
+ // Put in cache
+ cache[radius] = kernel;
+
+ return kernel;
+}
+//--------------------------------------------------------------------
+
+} // end namespace
+
+#endif //#define CLITKSUVPEAKGENERICFILTER_TXX
git://git.creatis.insa-lyon.fr / clitk.git / commitdiff