From a799140aa74ccd672430fdcaab4fc43ef3271645 Mon Sep 17 00:00:00 2001 From: Agata Krason Date: Tue, 27 May 2014 19:27:30 +0200 Subject: [PATCH] First version of new tool clitkDVH - possibility to compute dose volume histogram %/cc/voxels. --- tools/CMakeLists.txt | 19 +- tools/clitkDVH.cxx | 47 +++++ tools/clitkDVH.ggo | 22 +++ tools/clitkDVHGenericFilter.cxx | 120 ++++++++++++ tools/clitkDVHGenericFilter.h | 136 +++++++++++++ tools/clitkDVHGenericFilter.txx | 338 ++++++++++++++++++++++++++++++++ 6 files changed, 668 insertions(+), 14 deletions(-) create mode 100644 tools/clitkDVH.cxx create mode 100644 tools/clitkDVH.ggo create mode 100644 tools/clitkDVHGenericFilter.cxx create mode 100644 tools/clitkDVHGenericFilter.h create mode 100644 tools/clitkDVHGenericFilter.txx diff --git a/tools/CMakeLists.txt b/tools/CMakeLists.txt index 429b751..0b12205 100644 --- a/tools/CMakeLists.txt +++ b/tools/CMakeLists.txt @@ -36,11 +36,6 @@ IF (CLITK_BUILD_TOOLS) # ) SET(TOOLS_INSTALL clitkDicomInfo) - WRAP_GGO(clitkPointRigidRegistration_GGO_C clitkPointRigidRegistration.ggo) - ADD_EXECUTABLE(clitkPointRigidRegistration clitkPointRigidRegistration.cxx ${clitkPointRigidRegistration_GGO_C}) - TARGET_LINK_LIBRARIES(clitkPointRigidRegistration clitkCommon) - SET(TOOLS_INSTALL ${TOOLS_INSTALL} clitkPointRigidRegistration) - WRAP_GGO(clitkDicom2Image_GGO_C clitkDicom2Image.ggo) ADD_EXECUTABLE(clitkDicom2Image clitkDicom2Image.cxx ${clitkDicom2Image_GGO_C}) TARGET_LINK_LIBRARIES(clitkDicom2Image clitkCommon) @@ -237,6 +232,11 @@ IF (CLITK_BUILD_TOOLS) TARGET_LINK_LIBRARIES(clitkImageStatistics clitkCommon) SET(TOOLS_INSTALL ${TOOLS_INSTALL} clitkImageStatistics) + WRAP_GGO(clitkDVH_GGO_C clitkDVH.ggo) + ADD_EXECUTABLE(clitkDVH clitkDVH.cxx clitkDVHGenericFilter.cxx ${clitkDVH_GGO_C}) + TARGET_LINK_LIBRARIES(clitkDVH clitkVectorArithmLib clitkCommon) + SET(TOOLS_INSTALL ${TOOLS_INSTALL} clitkDVH) + WRAP_GGO(clitkVFConvert_GGO_C clitkVFConvert.ggo) ADD_EXECUTABLE(clitkVFConvert clitkVFConvert.cxx clitkVFConvertGenericFilter.cxx ${clitkVFConvert_GGO_C}) TARGET_LINK_LIBRARIES(clitkVFConvert clitkCommon ) @@ -305,11 +305,6 @@ IF (CLITK_BUILD_TOOLS) TARGET_LINK_LIBRARIES(clitkImageUncertainty clitkCommon ${ITK_LIBRARIES}) SET(TOOLS_INSTALL ${TOOLS_INSTALL} clitkImageUncertainty) - WRAP_GGO(clitkNormalizeImageFilter_GGO_C clitkNormalizeImageFilter.ggo) - ADD_EXECUTABLE(clitkNormalizeImageFilter clitkNormalizeImageFilter.cxx ${clitkNormalizeImageFilter_GGO_C}) - TARGET_LINK_LIBRARIES(clitkNormalizeImageFilter clitkCommon ) - SET(TOOLS_INSTALL ${TOOLS_INSTALL} clitkNormalizeImageFilter) - WRAP_GGO(clitkImageGradientMagnitude_GGO_C clitkImageGradientMagnitude.ggo) ADD_EXECUTABLE(clitkImageGradientMagnitude clitkImageGradientMagnitude.cxx ${clitkImageGradientMagnitude_GGO_C}) TARGET_LINK_LIBRARIES(clitkImageGradientMagnitude clitkCommon ) @@ -325,10 +320,6 @@ IF (CLITK_BUILD_TOOLS) TARGET_LINK_LIBRARIES(clitkImageIntensityWindowing clitkCommon ) SET(TOOLS_INSTALL ${TOOLS_INSTALL} clitkImageIntensityWindowing) - WRAP_GGO(clitkBlurImage_GGO_C clitkBlurImage.ggo) - ADD_EXECUTABLE(clitkBlurImage clitkBlurImage.cxx ${clitkBlurImage_GGO_C}) - TARGET_LINK_LIBRARIES(clitkBlurImage clitkCommon ) - SET(TOOLS_INSTALL ${TOOLS_INSTALL} clitkBlurImage) #========================================================= option(CLITK_USE_ROOT "Build experimental tools using root" OFF) diff --git a/tools/clitkDVH.cxx b/tools/clitkDVH.cxx new file mode 100644 index 0000000..bf069f4 --- /dev/null +++ b/tools/clitkDVH.cxx @@ -0,0 +1,47 @@ +/*========================================================================= + 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 +===========================================================================**/ + +/* ================================================= + * @file clitkHGenericFilter.txx + * @author Agata Krason + * @date 19 November 2013 + * + * @brief Dose Histogram + * + ===================================================*/ + +// clitk +#include "clitkDVH_ggo.h" +#include "clitkIO.h" +#include "clitkDVHGenericFilter.h" + +//-------------------------------------------------------------------- +int main(int argc, char * argv[]) +{ + + // Init command line + GGO(clitkDVH, args_info); + CLITK_INIT; + + // Filter + clitk::DVHGenericFilter::Pointer genericFilter=clitk::DVHGenericFilter::New(); + genericFilter->SetArgsInfo(args_info); + CLITK_TRY_CATCH_EXIT(genericFilter->Update()); + + return EXIT_SUCCESS; +}// end main diff --git a/tools/clitkDVH.ggo b/tools/clitkDVH.ggo new file mode 100644 index 0000000..cc439fc --- /dev/null +++ b/tools/clitkDVH.ggo @@ -0,0 +1,22 @@ +#File clitkDVH.ggo +package "clitkDVH" +version "2.0" +purpose "Dose volume histogram" + + +option "config" - "Config file" string no +option "verbose" v "Verbose" flag off + +option "input" i "Input image filename" string yes multiple +option "channel" c "Image channel to be used in statistics (-1 to process all channels)" int no default="-1" +option "mask" m "Mask image filename (uchar)" string yes +option "label" l "Label(s) in the mask image to consider" int no multiple default="1" +option "histogram" - "Compute histogram, allows median calculation" string no +option "dvhistogram" - "Compute dose volume histogram " string no +option "bins" - "Number of histogram bins" int no default="80" +option "lower" - "Lower histogram bound" double no default="0" +option "upper" - "Upper histogram bound" double no default="80" +option "allow_resize" r "Resize mask if different from input" flag off + + + diff --git a/tools/clitkDVHGenericFilter.cxx b/tools/clitkDVHGenericFilter.cxx new file mode 100644 index 0000000..f781c49 --- /dev/null +++ b/tools/clitkDVHGenericFilter.cxx @@ -0,0 +1,120 @@ +/*========================================================================= + 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 clitkDVHGenericFilter_cxx +#define clitkDVHGenericFilter_cxx + +#include "clitkDVHGenericFilter.h" + +namespace clitk +{ + + //----------------------------------------------------------- + // Constructor + //----------------------------------------------------------- + DVHGenericFilter::DVHGenericFilter() + { + m_Verbose=false; + m_InputFileName=""; + } + //----------------------------------------------------------- + + //----------------------------------------------------------- + // Update + //----------------------------------------------------------- + void DVHGenericFilter::Update() + { + // Read the Dimension and PixelType + int Dimension, Components; + std::string PixelType; + ReadImageDimensionAndPixelType(m_InputFileName, Dimension, PixelType, Components); + + if (m_ArgsInfo.channel_arg < -1 || m_ArgsInfo.channel_arg >= Components) { + std::cout << "Invalid image channel" << std::endl; + return; + } + + if (m_ArgsInfo.mask_given) { + int maskDimension, maskComponents; + std::string maskPixelType; + ReadImageDimensionAndPixelType(m_ArgsInfo.mask_arg, maskDimension, maskPixelType, maskComponents); + if (!(maskDimension == Dimension || maskDimension == (Dimension - 1))) { + std::cout << "Dimension of label mask must be equal to the (d)imension of the input image or d-1." << std::endl; + return; + } + } + + + // Call UpdateWithDim + if (Dimension==2) { + switch (Components) { + case 1: + UpdateWithDim<2,1>(PixelType); + break; + case 2: + UpdateWithDim<2,2>(PixelType); + break; + case 3: + UpdateWithDim<2,3>(PixelType); + break; + default: + std::cout << "Unsupported number of channels" << std::endl; + break; + } + } + else if (Dimension==3) { + switch (Components) { + case 1: + UpdateWithDim<3,1>(PixelType); + break; + case 2: + UpdateWithDim<3,2>(PixelType); + break; + case 3: + UpdateWithDim<3,3>(PixelType); + break; + default: + std::cout << "Unsupported number of channels" << std::endl; + break; + } + } + else if (Dimension==4) { + switch (Components) { + case 1: + UpdateWithDim<4,1>(PixelType); + break; + case 2: + UpdateWithDim<4,2>(PixelType); + break; + case 3: + UpdateWithDim<4,3>(PixelType); + break; + default: + std::cout << "Unsupported number of channels" << std::endl; + break; + } + } + else { + std::cout<<"Error, Only for 2 or 3 Dimensions!!!"< +class ITK_EXPORT DVHGenericFilter: public itk::LightObject +{ + +public: + + + //-------------------------------------------------------------------- + typedef DVHGenericFilter Self; + typedef itk::LightObject Superclass; + typedef itk::SmartPointer Pointer; + typedef itk::SmartPointer ConstPointer; + + //-------------------------------------------------------------------- + // Method for creation through the object factory + // and Run-time type information (and related methods) + itkNewMacro(Self); + itkTypeMacro(DVHGenericFilter, LightObject); + + //-------------------------------------------------------------------- + void SetArgsInfo(const args_info_clitkDVH & a) + { + m_ArgsInfo=a; + m_Verbose=m_ArgsInfo.verbose_flag; + + if(m_ArgsInfo.input_given) + m_InputFileName=m_ArgsInfo.input_arg[0]; + else if(m_ArgsInfo.inputs_num>0) + m_InputFileName=m_ArgsInfo.inputs[0]; + else { + std::cerr << "You must give an input file name" << std::endl; + exit(1); + } + } + //-------------------------------------------------------------------- + // Main function called each time the filter is updated + //template + //void UpdateWithInputImageType(); + void Update(); + + // Set methods + void SetDefaultPixelValue (double value) { mDefaultPixelValue = value ;} + void SetTypeOfOperation (int value) { mTypeOfOperation = value ;} + void SetScalar (double value) { mScalar = value ;} + void EnableOverwriteInputImage(bool b); + + // Get methods + double GetDefaultPixelValue () { return mDefaultPixelValue ;} + int GetTypeOfOperation () { return mTypeOfOperation ;} + double GetScalar () { return mScalar ;} + +protected: + + //-------------------------------------------------------------------- + // Constructor & Destructor + DVHGenericFilter(); + ~DVHGenericFilter() {}; + // template void InitializeImageType(); + // args_info_type mArgsInfo; + + //-------------------------------------------------------------------- + // Tempated members + //-------------------------------------------------------------------- + template void UpdateWithDim(std::string PixelType); + template void UpdateWithDimAndPixelType(); + + //-------------------------------------------------------------------- + // Data members + //-------------------------------------------------------------------- + args_info_clitkDVH m_ArgsInfo; + bool m_Verbose; + std::string m_InputFileName; + + template void InitializeImageType(); + bool mIsOperationUseASecondImage; + double mScalar; + double mDefaultPixelValue; + int mTypeOfOperation; + bool mOverwriteInputImage; + bool mOutputIsFloat; + + template void ComputeImage(Iter1 it, Iter2 ito); + + template void ComputeImage(Iter1 it1, Iter2 it2, Iter3 ito); + + +}; // end class +//-------------------------------------------------------------------- + +} // end namespace clitk + +#ifndef ITK_MANUAL_INSTANTIATION +#include "clitkDVHGenericFilter.txx" +#endif + +#endif // #define clitkDVHGenericFilter_h diff --git a/tools/clitkDVHGenericFilter.txx b/tools/clitkDVHGenericFilter.txx new file mode 100644 index 0000000..bd17a94 --- /dev/null +++ b/tools/clitkDVHGenericFilter.txx @@ -0,0 +1,338 @@ +/*========================================================================= + 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 clitkDVHGenericFilter_txx +#define clitkDVHGenericFilter_txx + +/* ================================================= + * @file clitkDVHGenericFilter.txx + * @author Agata Krason + * @date 20 November 2013 + * + * @brief Dose volume and image histogram + * + ===================================================*/ + +// itk include +#include "itkBinaryThresholdImageFilter.h" +#include "itkMaskImageFilter.h" +#include "itkMaskNegatedImageFilter.h" +#include "itkNthElementImageAdaptor.h" +#include "itkJoinSeriesImageFilter.h" +#include "itkMinimumMaximumImageCalculator.h" + +// clitk include +#include +#include "clitkImageCommon.h" +#include "clitkDVHGenericFilter.h" +#include "clitkCropLikeImageFilter.h" +#include "clitkResampleImageWithOptionsFilter.h" + +//------------------------------------------------------------------- +// +namespace clitk +{ + +//------------------------------------------------------------------- +template +void +DVHGenericFilter::UpdateWithDim(std::string PixelType) +{ + if (m_Verbose) std::cout << "Image was detected to be "<(); + } + else if(PixelType == "unsigned_short"){ + if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and unsigned_short..." << std::endl; + UpdateWithDimAndPixelType(); + } + + else if (PixelType == "unsigned_char"){ + if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and unsigned_char..." << std::endl; + UpdateWithDimAndPixelType(); + } + + else if(PixelType == "double"){ + if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and double..." << std::endl; + UpdateWithDimAndPixelType(); + } + else { + if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and float..." << std::endl; + UpdateWithDimAndPixelType(); + } +} + + +//------------------------------------------------------------------- +// Update with the number of dimensions and the pixeltype +//------------------------------------------------------------------- +template +void +DVHGenericFilter::UpdateWithDimAndPixelType() +{ + + // ImageTypes + typedef unsigned char LabelPixelType; + typedef itk::Image, Dimension> InputImageType; + typedef itk::Image LabelImageType; + + // Read the input + typedef itk::ImageFileReader InputReaderType; + typename InputReaderType::Pointer reader = InputReaderType::New(); + reader->SetFileName( m_InputFileName); + reader->Update(); + typename InputImageType::Pointer input= reader->GetOutput(); + + typedef itk::NthElementImageAdaptor InputImageAdaptorType; + typedef itk::Image OutputImageType; + + typename InputImageAdaptorType::Pointer input_adaptor = InputImageAdaptorType::New(); + input_adaptor->SetImage(input); + + // Filter + typedef itk::LabelStatisticsImageFilter StatisticsImageFilterType; + typename StatisticsImageFilterType::Pointer statisticsFilter=StatisticsImageFilterType::New(); + statisticsFilter->SetInput(input_adaptor); + + // Label image + typename LabelImageType::Pointer labelImage; + if (m_ArgsInfo.mask_given) { + int maskDimension, maskComponents; + std::string maskPixelType; + ReadImageDimensionAndPixelType(m_ArgsInfo.mask_arg, maskDimension, maskPixelType, maskComponents); + + if (maskDimension == Dimension - 1) { + // Due to a limitation of filter itk::LabelStatisticsImageFilter, InputImageType and LabelImageType + // must have the same image dimension. However, we want to support label images with Dl = Di - 1, + // so we need to replicate the label image as many times as the size along dimension Di. + if (m_Verbose) + std::cout << "Replicating label image to match input image's dimension... " << std::endl; + + typedef itk::Image ReducedLabelImageType; + typedef itk::ImageFileReader LabelImageReaderType; + typedef itk::JoinSeriesImageFilter JoinImageFilterType; + + + typename LabelImageReaderType::Pointer labelImageReader=LabelImageReaderType::New(); + labelImageReader->SetFileName(m_ArgsInfo.mask_arg); + labelImageReader->Update(); + + typename JoinImageFilterType::Pointer joinFilter = JoinImageFilterType::New(); + typename InputImageType::SizeType size = input->GetLargestPossibleRegion().GetSize(); + for (unsigned int i = 0; i < size[Dimension - 1]; i++) + joinFilter->PushBackInput(labelImageReader->GetOutput()); + + joinFilter->Update(); + labelImage = joinFilter->GetOutput(); + } + else { + typedef itk::ImageFileReader LabelImageReaderType; + typename LabelImageReaderType::Pointer labelImageReader=LabelImageReaderType::New(); + labelImageReader->SetFileName(m_ArgsInfo.mask_arg); + labelImageReader->Update(); + labelImage= labelImageReader->GetOutput(); + + // Check mask sampling/size + if (!HaveSameSizeAndSpacing(labelImage, input)) { + if (m_ArgsInfo.allow_resize_flag) { + if (m_ArgsInfo.verbose_flag) { + std::cout << "Resize mask image like input" << std::endl; + } + + typedef clitk::ResampleImageWithOptionsFilter ResamplerType; + typename ResamplerType::Pointer resampler = ResamplerType::New(); + resampler->SetInput(labelImage); + resampler->SetOutputSpacing(input->GetSpacing()); + resampler->SetGaussianFilteringEnabled(false); + resampler->Update(); + labelImage = resampler->GetOutput(); + labelImage->GetSpacing(); + typedef clitk::CropLikeImageFilter FilterType; + typename FilterType::Pointer crop = FilterType::New(); + crop->SetInput(labelImage); + crop->SetCropLikeImage(input); + crop->Update(); + labelImage = crop->GetOutput(); + //writeImage(labelImage, "test2.mha"); + + } + else { + std::cerr << "Mask image has a different size/spacing than input. Abort" << std::endl; + exit(-1); + } + } + + } + + } + else { + labelImage=LabelImageType::New(); + labelImage->SetRegions(input->GetLargestPossibleRegion()); + labelImage->SetOrigin(input->GetOrigin()); + labelImage->SetSpacing(input->GetSpacing()); + labelImage->Allocate(); + labelImage->FillBuffer(m_ArgsInfo.label_arg[0]); + } + statisticsFilter->SetLabelInput(labelImage); + + // Check/compute spacing + const typename LabelImageType::SpacingType& spacing = input->GetSpacing(); + double spacing_cc = (spacing[0]*spacing[1]*spacing[2])/1000; + // std::cout<<"Spacing x : "<SelectNthElement(c); + input_adaptor->Update(); + + for (unsigned int k=0; k< numberOfLabels; k++) { + label=m_ArgsInfo.label_arg[k]; + // Histogram + if (m_ArgsInfo.histogram_given) + { + std::cout<<"--------------"<SetUseHistograms(true); + statisticsFilter->SetHistogramParameters(m_ArgsInfo.bins_arg, m_ArgsInfo.lower_arg, m_ArgsInfo.upper_arg); + } + + // DVHistogram + if(m_ArgsInfo.dvhistogram_given) + { + statisticsFilter->SetUseHistograms(true); + statisticsFilter->SetHistogramParameters(m_ArgsInfo.bins_arg, m_ArgsInfo.lower_arg, m_ArgsInfo.upper_arg); + } + + statisticsFilter->Update(); + + // Histogram + if (m_ArgsInfo.histogram_given) + { + if (m_Verbose) std::cout<<"Median: "; + std::cout<GetMedian(label)<GetHistogram(label); + + // Screen + if (m_Verbose) std::cout<<"Histogram: "<GetBinMin(0,i)<<"\t"<GetMeasurement(i,0)<<"\t"<GetBinMax(0,i)<<"\t"<GetFrequency(i)<GetBinMin(0,i)<<"\t"<GetMeasurement(i,0)<<"\t"<GetBinMax(0,i)<<"\t"<GetFrequency(i)<GetHistogram(label); + + // Screen + std::cout<<"Total volume : "; + std::cout<GetCount(label)<<" [No. of voxels]"<GetCount(label))*spacing_cc)<<" [cc]"<GetMean(label)<<" [Gy]"<GetMinimum(label)<<" [Gy]"<GetMaximum(label)<<" [Gy]"<GetFrequency(i))*100)/(statisticsFilter->GetCount(label)); + double popCumulativeVolume = 0; + for(int j=0; jGetFrequency(j)); + } + double cumulativeVolume = popCumulativeVolume + (dvhistogram->GetFrequency(i)); + double percentCumulativeVolume =(cumulativeVolume*100)/(statisticsFilter->GetCount(label)) ; + double ccCumulativeVolume = (popCumulativeVolume + (dvhistogram->GetFrequency(i)))*spacing_cc; + std::cout<GetBinMax(0,i)<<"\t "<GetFrequency(i)<<"\t "<GetFrequency(i))*spacing_cc)<<"\t "<GetCount(label)<<" [No. of voxels]"<GetCount(label))*spacing_cc)<<" [cc]"<GetMean(label)<<" [Gy]"<GetMinimum(label)<<" [Gy]"<GetMaximum(label)<<" [Gy]"<GetFrequency(i))*100)/(statisticsFilter->GetCount(label)); + double popCumulativeVolume = 0; + for(int j=0; jGetFrequency(j)); + } + double cumulativeVolume = popCumulativeVolume + (dvhistogram->GetFrequency(i)); + double percentCumulativeVolume =(cumulativeVolume*100)/(statisticsFilter->GetCount(label)) ; + double ccCumulativeVolume = (popCumulativeVolume + (dvhistogram->GetFrequency(i)))*spacing_cc; + dvhistogramFile<GetBinMax(0,i)<<"\t "<GetFrequency(i)<<"\t "<GetFrequency(i))*spacing_cc)<<"\t "<