From: tbaudier Date: Tue, 13 Nov 2018 07:53:42 +0000 (+0100) Subject: Merge branch 'master' of git://git.creatis.insa-lyon.fr/clitk X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=commitdiff_plain;h=a9d632132fa5c0da47d7b3079ccca9d673381a1e;hp=88759fd450cc61a1b6818f80548c8632b2a3ffea;p=clitk.git Merge branch 'master' of git://git.creatis.insa-lyon.fr/clitk --- diff --git a/tools/CMakeLists.txt b/tools/CMakeLists.txt index 4259e19..ca366af 100644 --- a/tools/CMakeLists.txt +++ b/tools/CMakeLists.txt @@ -270,6 +270,16 @@ if(CLITK_BUILD_TOOLS) 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 ) diff --git a/tools/clitkRTStructStatistics.cxx b/tools/clitkRTStructStatistics.cxx new file mode 100644 index 0000000..504c779 --- /dev/null +++ b/tools/clitkRTStructStatistics.cxx @@ -0,0 +1,46 @@ +/*=========================================================================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 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 diff --git a/tools/clitkRTStructStatistics.ggo b/tools/clitkRTStructStatistics.ggo new file mode 100644 index 0000000..ebd275b --- /dev/null +++ b/tools/clitkRTStructStatistics.ggo @@ -0,0 +1,9 @@ +#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 diff --git a/tools/clitkRTStructStatisticsGenericFilter.cxx b/tools/clitkRTStructStatisticsGenericFilter.cxx new file mode 100644 index 0000000..bff6956 --- /dev/null +++ b/tools/clitkRTStructStatisticsGenericFilter.cxx @@ -0,0 +1,31 @@ +/*========================================================================= + 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; + + +} + +#endif //define CLITKRTSTRUCTSTATISTICSGENERICFILTER_CXX diff --git a/tools/clitkRTStructStatisticsGenericFilter.h b/tools/clitkRTStructStatisticsGenericFilter.h new file mode 100644 index 0000000..4339c69 --- /dev/null +++ b/tools/clitkRTStructStatisticsGenericFilter.h @@ -0,0 +1,84 @@ +/*========================================================================= + 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 + * @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 ITK_EXPORT RTStructStatisticsGenericFilter: + public clitk::ImageToImageGenericFilter > { + + public: + + // Constructor + RTStructStatisticsGenericFilter (); + + // Types + typedef RTStructStatisticsGenericFilter Self; + typedef ImageToImageGenericFilterBase Superclass; + typedef itk::SmartPointer Pointer; + typedef itk::SmartPointer ConstPointer; + + // New + itkNewMacro(Self); + + + + //-------------------------------------------------------------------- + void SetArgsInfo(const args_info_type & a); + + //-------------------------------------------------------------------- + // Main function called each time the filter is updated + template + void UpdateWithInputImageType(); + + protected: + template void InitializeImageType(); + args_info_type mArgsInfo; + + }; // end class RTStructStatisticsGenericFilter + +} // end namespace +//-------------------------------------------------------------------- + + +#ifndef ITK_MANUAL_INSTANTIATION +#include "clitkRTStructStatisticsGenericFilter.txx" +#endif + +#endif //#define CLITKRTSTRUCTSTATISTICSGENERICFILTER_H + diff --git a/tools/clitkRTStructStatisticsGenericFilter.txx b/tools/clitkRTStructStatisticsGenericFilter.txx new file mode 100644 index 0000000..9418ac5 --- /dev/null +++ b/tools/clitkRTStructStatisticsGenericFilter.txx @@ -0,0 +1,124 @@ +/*========================================================================= + 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 +RTStructStatisticsGenericFilter::RTStructStatisticsGenericFilter() + :ImageToImageGenericFilter("RTStructStatisticsGenericFilter") +{ + InitializeImageType<2>(); + InitializeImageType<3>(); + InitializeImageType<4>(); +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +template +void RTStructStatisticsGenericFilter::InitializeImageType() +{ + ADD_DEFAULT_IMAGE_TYPES(Dim); +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void RTStructStatisticsGenericFilter::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 +template +void RTStructStatisticsGenericFilter::UpdateWithInputImageType() +{ + // Read mask input + typedef itk::Image MaskInputImageType; + typename MaskInputImageType::Pointer mask = NULL; + mask = this->template GetInput(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 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 diff --git a/tools/clitkSUVPeak.cxx b/tools/clitkSUVPeak.cxx new file mode 100644 index 0000000..0ad4000 --- /dev/null +++ b/tools/clitkSUVPeak.cxx @@ -0,0 +1,46 @@ +/*========================================================================= + 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 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 diff --git a/tools/clitkSUVPeak.ggo b/tools/clitkSUVPeak.ggo new file mode 100644 index 0000000..6e9206a --- /dev/null +++ b/tools/clitkSUVPeak.ggo @@ -0,0 +1,11 @@ +#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 diff --git a/tools/clitkSUVPeakGenericFilter.cxx b/tools/clitkSUVPeakGenericFilter.cxx new file mode 100644 index 0000000..49576dc --- /dev/null +++ b/tools/clitkSUVPeakGenericFilter.cxx @@ -0,0 +1,31 @@ +/*========================================================================= + 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; + + +} + +#endif //define CLITKSUVPEAKGENERICFILTER_CXX diff --git a/tools/clitkSUVPeakGenericFilter.h b/tools/clitkSUVPeakGenericFilter.h new file mode 100644 index 0000000..57c4d05 --- /dev/null +++ b/tools/clitkSUVPeakGenericFilter.h @@ -0,0 +1,89 @@ +/*========================================================================= + 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 + * @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 ITK_EXPORT SUVPeakGenericFilter: + public clitk::ImageToImageGenericFilter > { + + public: + + // Constructor + SUVPeakGenericFilter (); + + // Types + typedef SUVPeakGenericFilter Self; + typedef ImageToImageGenericFilterBase Superclass; + typedef itk::SmartPointer Pointer; + typedef itk::SmartPointer ConstPointer; + + // New + itkNewMacro(Self); + + + + //-------------------------------------------------------------------- + void SetArgsInfo(const args_info_type & a); + + //SUVPeak + template typename ImageType::Pointer ComputeMeanFilterKernel(const typename ImageType::SpacingType & spacing, double radius); + + //-------------------------------------------------------------------- + // Main function called each time the filter is updated + template + void UpdateWithInputImageType(); + + protected: + template void InitializeImageType(); + args_info_type mArgsInfo; + + }; // end class SUVPeakGenericFilter + +} // end namespace +//-------------------------------------------------------------------- + + +#ifndef ITK_MANUAL_INSTANTIATION +#include "clitkSUVPeakGenericFilter.txx" +#endif + +#endif //#define CLITKSUVPEAKGENERICFILTER_H + diff --git a/tools/clitkSUVPeakGenericFilter.txx b/tools/clitkSUVPeakGenericFilter.txx new file mode 100644 index 0000000..02ca48d --- /dev/null +++ b/tools/clitkSUVPeakGenericFilter.txx @@ -0,0 +1,229 @@ +/*========================================================================= + 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 + +namespace clitk +{ + +//-------------------------------------------------------------------- +template +SUVPeakGenericFilter::SUVPeakGenericFilter() + :ImageToImageGenericFilter("SUVPeakGenericFilter") +{ + InitializeImageType<2>(); + InitializeImageType<3>(); + InitializeImageType<4>(); +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +template +void SUVPeakGenericFilter::InitializeImageType() +{ + ADD_DEFAULT_IMAGE_TYPES(Dim); +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void SUVPeakGenericFilter::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 +template +void SUVPeakGenericFilter::UpdateWithInputImageType() +{ + // Read input + typename ImageType::Pointer input = this->template GetInput(0); + + //Read mask + typedef itk::Image MaskImageType; + typename MaskImageType::Pointer mask = NULL; + if(mArgsInfo.mask_given) { + mask = this->template GetInput(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(input->GetSpacing(), radius); + + // Perform the convolution + typedef itk::ConvolutionImageFilter FilterType; + typename FilterType::Pointer filter = FilterType::New(); + filter->SetInput(input); + filter->SetKernelImage(kernel); + filter->Update(); + typename ImageType::Pointer output = filter->GetOutput(); + + + typedef itk::ImageRegionConstIteratorWithIndex IteratorType; + typedef itk::ImageRegionConstIteratorWithIndex 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 +template +typename ImageType::Pointer +SUVPeakGenericFilter::ComputeMeanFilterKernel(const typename ImageType::SpacingType & spacing, double radius) +{ + // Some kind of cache to speed up a bit + static std::map 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 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