From 231522feb65774b26310720ed476989c1f349b86 Mon Sep 17 00:00:00 2001 From: tbaudier Date: Mon, 21 Nov 2016 14:46:08 +0100 Subject: [PATCH] Add gamma index tool for 3D images Quite time-consuming even if it uses requested region instead of largest region in ITK filters --- tools/CMakeLists.txt | 5 + tools/clitkGammaIndex3D.cxx | 46 +++++ tools/clitkGammaIndex3D.ggo | 16 ++ tools/clitkGammaIndex3DGenericFilter.cxx | 31 ++++ tools/clitkGammaIndex3DGenericFilter.h | 96 ++++++++++ tools/clitkGammaIndex3DGenericFilter.txx | 216 +++++++++++++++++++++++ 6 files changed, 410 insertions(+) create mode 100644 tools/clitkGammaIndex3D.cxx create mode 100644 tools/clitkGammaIndex3D.ggo create mode 100644 tools/clitkGammaIndex3DGenericFilter.cxx create mode 100644 tools/clitkGammaIndex3DGenericFilter.h create mode 100644 tools/clitkGammaIndex3DGenericFilter.txx diff --git a/tools/CMakeLists.txt b/tools/CMakeLists.txt index 422d8f9..5d67f4b 100644 --- a/tools/CMakeLists.txt +++ b/tools/CMakeLists.txt @@ -90,6 +90,11 @@ if(CLITK_BUILD_TOOLS) add_executable(clitkSplitImage clitkSplitImage.cxx clitkSplitImageGenericFilter.cxx ${clitkSplitImage_GGO_C}) target_link_libraries(clitkSplitImage clitkCommon ) set(TOOLS_INSTALL ${TOOLS_INSTALL} clitkSplitImage) + + WRAP_GGO(clitkGammaIndex3D_GGO_C clitkGammaIndex3D.ggo) + add_executable(clitkGammaIndex3D clitkGammaIndex3D.cxx clitkGammaIndex3DGenericFilter.cxx ${clitkGammaIndex3D_GGO_C}) + target_link_libraries(clitkGammaIndex3D clitkCommon ) + set(TOOLS_INSTALL ${TOOLS_INSTALL} clitkGammaIndex3D) WRAP_GGO(clitkWriteDicomSeries_GGO_C clitkWriteDicomSeries.ggo) add_executable(clitkWriteDicomSeries clitkWriteDicomSeries.cxx ${clitkWriteDicomSeries_GGO_C}) diff --git a/tools/clitkGammaIndex3D.cxx b/tools/clitkGammaIndex3D.cxx new file mode 100644 index 0000000..560c1cd --- /dev/null +++ b/tools/clitkGammaIndex3D.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 CLITKGAMMAINDEX3D_CXX +#define CLITKGAMMAINDEX3D_CXX + +// clitk include +#include "clitkGammaIndex3D_ggo.h" +#include "clitkGammaIndex3DGenericFilter.h" +#include "clitkIO.h" + +//-------------------------------------------------------------------- +int main(int argc, char * argv[]) +{ + + // Init command line + GGO(clitkGammaIndex3D, args_info); + CLITK_INIT; + + // Creation of a generic filter + typedef clitk::GammaIndex3DGenericFilter 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 CLITKGAMMAINDEX3D_CXX diff --git a/tools/clitkGammaIndex3D.ggo b/tools/clitkGammaIndex3D.ggo new file mode 100644 index 0000000..0f52d5a --- /dev/null +++ b/tools/clitkGammaIndex3D.ggo @@ -0,0 +1,16 @@ +#File clitkGammaIndex3D.ggo +package "clitkGammaIndex3D" +version "1.0" +purpose "Compute the gamma index for 2D and 3D images" + + +option "config" - "Config file" string no +option "verbose" v "Verbose" flag off +option "imagetypes" - "Display allowed image types" flag off + +option "input1" i "Input first dose image filename" string yes +option "input2" j "Input second dose image filename" string yes +option "output" o "Output image filename" string yes + +option "distance" d "Distance criteria (mm)" double yes +option "dose" p "Dose criteria (%)" double yes diff --git a/tools/clitkGammaIndex3DGenericFilter.cxx b/tools/clitkGammaIndex3DGenericFilter.cxx new file mode 100644 index 0000000..da944fb --- /dev/null +++ b/tools/clitkGammaIndex3DGenericFilter.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 CLITKGAMMAINDEX3DGENERICFILTER_CXX +#define CLITKGAMMAINDEX3DGENERICFILTER_CXX + +#include "clitkGammaIndex3DGenericFilter.h" + +namespace clitk { + // Specialisation +// template<> +// class GammaIndex3DGenericFilter; + + +} + +#endif //define CLITKGAMMAINDEX3DGENERICFILTER_CXX diff --git a/tools/clitkGammaIndex3DGenericFilter.h b/tools/clitkGammaIndex3DGenericFilter.h new file mode 100644 index 0000000..79c1f45 --- /dev/null +++ b/tools/clitkGammaIndex3DGenericFilter.h @@ -0,0 +1,96 @@ +/*========================================================================= + 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 CLITKGAMMAINDEX3DGENERICFILTER_H +#define CLITKGAMMAINDEX3DGENERICFILTER_H +/** + ------------------------------------------------------------------- + * @file clitkGammaIndex3DGenericFilter.h + * @author David Sarrut + * @date 23 Feb 2008 08:37:53 + + * @brief + -------------------------------------------------------------------*/ + +// clitk include +#include "clitkCommon.h" +#include "clitkImageToImageGenericFilter.h" +#include "clitkGammaIndex3D_ggo.h" + +// itk include +#include "itkImage.h" +#include "itkImageIOBase.h" +#include "itkImageRegionIterator.h" +#include "itkImageRegionConstIterator.h" + +//-------------------------------------------------------------------- +namespace clitk { + + template + class ITK_EXPORT GammaIndex3DGenericFilter: + public clitk::ImageToImageGenericFilter > { + + public: + + // Constructor + GammaIndex3DGenericFilter (); + + // Types + typedef GammaIndex3DGenericFilter Self; + typedef ImageToImageGenericFilterBase Superclass; + typedef itk::SmartPointer Pointer; + typedef itk::SmartPointer ConstPointer; + + // New + itkNewMacro(Self); + + //-------------------------------------------------------------------- + void SetArgsInfo(const args_info_type & a); + + // Set methods + void SetDistance (double distance) { mDistance = distance ;} + void SetDose (double dose) { mDose = dose ;} + + // Get methods + double GetDistance () { return mDistance ;} + double GetDose () { return mDose ;} + + //-------------------------------------------------------------------- + // Main function called each time the filter is updated + template + void UpdateWithInputImageType(); + + protected: + template void InitializeImageType(); + double mDistance; + double mDose; + args_info_type mArgsInfo; + + //-------------------------------------------------------------------- + + }; // end class GammaIndex3DGenericFilter + +} // end namespace +//-------------------------------------------------------------------- + + +#ifndef ITK_MANUAL_INSTANTIATION +#include "clitkGammaIndex3DGenericFilter.txx" +#endif + +#endif //#define CLITKGAMMAINDEX3DMGENERICFILTER_H + diff --git a/tools/clitkGammaIndex3DGenericFilter.txx b/tools/clitkGammaIndex3DGenericFilter.txx new file mode 100644 index 0000000..9137c9e --- /dev/null +++ b/tools/clitkGammaIndex3DGenericFilter.txx @@ -0,0 +1,216 @@ +/*========================================================================= + 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 CLITKGAMMAINDEX3DGENERICFILTER_TXX +#define CLITKGAMMAINDEX3DGENERICFILTER_TXX + +#include "clitkImageCommon.h" + +#include "itkApproximateSignedDistanceMapImageFilter.h" +#include "itkAddImageFilter.h" +#include +#include +#include +#include + +namespace clitk +{ + +//-------------------------------------------------------------------- +template +GammaIndex3DGenericFilter::GammaIndex3DGenericFilter() + :ImageToImageGenericFilter("GammaIndex3DGenericFilter"), mDose(3), mDistance(3) +{ + InitializeImageType<2>(); + InitializeImageType<3>(); + //InitializeImageType<4>(); +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +template +void GammaIndex3DGenericFilter::InitializeImageType() +{ + ADD_DEFAULT_IMAGE_TYPES(Dim); +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void GammaIndex3DGenericFilter::SetArgsInfo(const args_info_type & a) +{ + mArgsInfo=a; + + // Set value + this->SetIOVerbose(mArgsInfo.verbose_flag); + mDistance = mArgsInfo.distance_arg; + mDose = mArgsInfo.dose_arg; + + if (mArgsInfo.imagetypes_flag) this->PrintAvailableImageTypes(); + + if (mArgsInfo.input1_given) this->AddInputFilename(mArgsInfo.input1_arg); + if (mArgsInfo.input2_given) this->AddInputFilename(mArgsInfo.input2_arg); + + if (mArgsInfo.output_given) this->SetOutputFilename(mArgsInfo.output_arg); +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +template +void GammaIndex3DGenericFilter::UpdateWithInputImageType() +{ + // Read input1 + typename ImageType::Pointer input1 = this->template GetInput(0); + + // Set input image iterator + typedef itk::ImageRegionIterator IteratorType; + IteratorType it(input1, input1->GetLargestPossibleRegion()); + + // typedef input2 + typename ImageType::Pointer input2 = this->template GetInput(1); + + //Create output + typename ImageType::Pointer output = ImageType::New(); + output->SetRegions(input1->GetLargestPossibleRegion()); + output->SetOrigin(input1->GetOrigin()); + output->SetSpacing(input1->GetSpacing()); + output->Allocate(); + + int sizeTotal = input1->GetLargestPossibleRegion().GetSize(0)*input1->GetLargestPossibleRegion().GetSize(1)*input1->GetLargestPossibleRegion().GetSize(2); + int sizeTemp(0); + double minSpacing(std::numeric_limits::max()); + for (unsigned int temp=0; temp < input1->GetSpacing().Size(); ++temp) + { + if (input1->GetSpacing()[temp] < minSpacing) + minSpacing = input1->GetSpacing()[temp]; + } + + + while (!it.IsAtEnd()) + { + double doseRef = it.Get(); + + //Create a region around it.Get() + typename ImageType::RegionType smallRegion; + typename ImageType::RegionType tempRegion; + smallRegion=input2->GetLargestPossibleRegion(); + + typename ImageType::SizeType sizeRegion; + sizeRegion.Fill(2*mDistance/minSpacing+1); //5 voxels de cote + + typename ImageType::OffsetType offsetRegion; + offsetRegion.Fill((sizeRegion[0]-1)/2); + + typename ImageType::IndexType startRegion; + startRegion=it.GetIndex(); + startRegion=startRegion-offsetRegion; + + tempRegion.SetIndex(startRegion); + tempRegion.SetSize(sizeRegion); + smallRegion.Crop(tempRegion); + + //compute the dose difference + typedef itk::AddImageFilter AddImageFilterType; + typename AddImageFilterType::Pointer addImageFilter = AddImageFilterType::New(); + addImageFilter->SetInput(input2); + addImageFilter->SetConstant2(-doseRef); + addImageFilter->GetOutput()->SetRequestedRegion(smallRegion); + addImageFilter->Update(); + + //compute the distance + typename ImageType::Pointer distanceMap = ImageType::New(); + distanceMap->SetRegions(smallRegion); + distanceMap->SetOrigin(input1->GetOrigin()); + distanceMap->SetSpacing(input1->GetSpacing()); + distanceMap->Allocate(); + distanceMap->FillBuffer(0.0); + //distanceMap->SetPixel(it.GetIndex(),1.0); + + typedef itk::ApproximateSignedDistanceMapImageFilter< ImageType, ImageType > ApproximateSignedDistanceMapImageFilterType; + typename ApproximateSignedDistanceMapImageFilterType::Pointer approximateSignedDistanceMapImageFilter = ApproximateSignedDistanceMapImageFilterType::New(); + approximateSignedDistanceMapImageFilter->SetInput(distanceMap); + approximateSignedDistanceMapImageFilter->SetInsideValue(1.0); + approximateSignedDistanceMapImageFilter->SetOutsideValue(0.0); + approximateSignedDistanceMapImageFilter->GetOutput()->SetRequestedRegion(smallRegion); + approximateSignedDistanceMapImageFilter->Update(); + + //Square and Divide by dose & distance criteria + typedef itk::SquareImageFilter SquareImageFilterType; + typename SquareImageFilterType::Pointer squareImageFilter = SquareImageFilterType::New(); + squareImageFilter->SetInput(addImageFilter->GetOutput()); + squareImageFilter->GetOutput()->SetRequestedRegion(smallRegion); + squareImageFilter->Update(); + + typedef itk::DivideImageFilter DivideImageFilterType; + typename DivideImageFilterType::Pointer divideImageFilter = DivideImageFilterType::New (); + divideImageFilter->SetInput1(squareImageFilter->GetOutput()); + divideImageFilter->SetConstant2(mDose*mDose); + divideImageFilter->GetOutput()->SetRequestedRegion(smallRegion); + divideImageFilter->Update(); + + typename SquareImageFilterType::Pointer squareImageFilter2 = SquareImageFilterType::New(); + squareImageFilter2->SetInput(approximateSignedDistanceMapImageFilter->GetOutput()); + squareImageFilter2->GetOutput()->SetRequestedRegion(smallRegion); + squareImageFilter2->Update(); + + typename DivideImageFilterType::Pointer divideImageFilter2 = DivideImageFilterType::New (); + divideImageFilter2->SetInput1(squareImageFilter2->GetOutput()); + divideImageFilter2->SetConstant2(mDistance*mDistance); + divideImageFilter2->GetOutput()->SetRequestedRegion(smallRegion); + divideImageFilter2->Update(); + + //Sum the 2 images and take the square root + typename AddImageFilterType::Pointer addImageFilter2 = AddImageFilterType::New(); + addImageFilter2->SetInput1(divideImageFilter2->GetOutput()); + addImageFilter2->SetInput2(divideImageFilter2->GetOutput()); + addImageFilter2->GetOutput()->SetRequestedRegion(smallRegion); + addImageFilter2->Update(); + + typedef itk::SqrtImageFilter SqrtFilterType; + typename SqrtFilterType::Pointer sqrtFilter = SqrtFilterType::New(); + sqrtFilter->SetInput(addImageFilter2->GetOutput()); + sqrtFilter->GetOutput()->SetRequestedRegion(smallRegion); + sqrtFilter->Update(); + + //find min + typedef itk::MinimumMaximumImageCalculator ImageCalculatorFilterType; + typename ImageCalculatorFilterType::Pointer imageCalculatorFilter = ImageCalculatorFilterType::New (); + imageCalculatorFilter->SetImage(sqrtFilter->GetOutput()); + //imageCalculatorFilter->GetOutput()->SetRequestedRegion(smallRegion); + imageCalculatorFilter->ComputeMinimum(); + + //Set the value + output->SetPixel(it.GetIndex(),imageCalculatorFilter->GetMinimum()); + ++it; + + ++sizeTemp; + std::cout << (double)sizeTemp / sizeTotal * 100.0 << std::endl; + } + + this->template SetNextOutput(output); + +} +//-------------------------------------------------------------------- + +} // end namespace + +#endif //#define CLITKGAMMAINDEX3DGENERICFILTER_TXX -- 2.45.1