From: Agata Krason <agatakrason@gmail.com>
Date: Tue, 27 May 2014 17:27:30 +0000 (+0200)
Subject: First version of new tool clitkDVH - possibility to compute dose volume histogram... 
X-Git-Tag: v1.4.0~121
X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=commitdiff_plain;h=a799140aa74ccd672430fdcaab4fc43ef3271645;p=clitk.git

First version of new tool clitkDVH - possibility to compute dose volume histogram %/cc/voxels.
---

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 <agata.krason@creatis.insa-lyon.fr>
+ * @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!!!"<<std::endl ;
+      return;
+    }
+  }
+
+
+} //end clitk
+
+#endif  //#define clitkDVHGenericFilter_cxx
diff --git a/tools/clitkDVHGenericFilter.h b/tools/clitkDVHGenericFilter.h
new file mode 100644
index 0000000..07b4d4c
--- /dev/null
+++ b/tools/clitkDVHGenericFilter.h
@@ -0,0 +1,136 @@
+/*=========================================================================
+  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_H
+#define CLITKDVHGENERICFILTER_H
+
+// clitk include
+#include "clitkIO.h"
+#include "clitkCommon.h"
+#include "clitkImageCommon.h" 
+#include "clitkImageToImageGenericFilter.h"
+#include "clitkDVH_ggo.h"
+
+// itk include
+#include "itkImage.h"
+#include "itkImageIOBase.h"
+#include "itkImageRegionIterator.h"
+#include "itkImageRegionConstIterator.h"
+#include "itkLightObject.h"
+#include "itkLabelStatisticsImageFilter.h"
+#include "itkLabelGeometryImageFilter.h"
+
+//--------------------------------------------------------------------
+namespace clitk
+{
+
+//template<class args_info_type>
+class ITK_EXPORT DVHGenericFilter: public itk::LightObject
+{
+
+public:
+
+
+    //--------------------------------------------------------------------
+    typedef DVHGenericFilter         Self;
+    typedef itk::LightObject               Superclass;
+    typedef itk::SmartPointer<Self>            Pointer;
+    typedef itk::SmartPointer<const Self>      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<class InputImageType>
+    //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<unsigned int Dim> void InitializeImageType();
+   // args_info_type mArgsInfo;
+
+   //--------------------------------------------------------------------
+   // Tempated members
+   //--------------------------------------------------------------------
+   template<unsigned int Dimension, unsigned int Components> void UpdateWithDim(std::string PixelType);
+   template <unsigned int Dimension, class PixelType, unsigned int Components> void UpdateWithDimAndPixelType();
+
+   //--------------------------------------------------------------------
+   // Data members
+   //--------------------------------------------------------------------
+   args_info_clitkDVH m_ArgsInfo;
+   bool m_Verbose;
+   std::string m_InputFileName;
+
+   template<unsigned int Dim> void InitializeImageType();
+   bool mIsOperationUseASecondImage;
+   double mScalar;
+   double mDefaultPixelValue;
+   int mTypeOfOperation;
+   bool mOverwriteInputImage;
+   bool mOutputIsFloat;
+
+   template<class Iter1, class Iter2> void ComputeImage(Iter1 it, Iter2 ito);
+
+   template<class Iter1, class Iter2, class Iter3> 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 <agata.krason@creatis.insa-lyon.fr>
+ * @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 <clitkCommon.h>
+#include "clitkImageCommon.h"
+#include "clitkDVHGenericFilter.h"
+#include "clitkCropLikeImageFilter.h"
+#include "clitkResampleImageWithOptionsFilter.h"
+
+//-------------------------------------------------------------------
+//
+namespace clitk
+{
+
+//-------------------------------------------------------------------
+template<unsigned int Dimension, unsigned int Components>
+void
+DVHGenericFilter::UpdateWithDim(std::string PixelType)
+{
+  if (m_Verbose) std::cout << "Image was detected to be "<<Dimension<<"D and "<< PixelType<<" with " << Components << " channel(s)..."<<std::endl;
+
+  if(PixelType == "short"){  
+    if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and signed short..." << std::endl;
+    UpdateWithDimAndPixelType<Dimension, signed short, Components>(); 
+  }
+  else if(PixelType == "unsigned_short"){  
+    if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and unsigned_short..." << std::endl;
+    UpdateWithDimAndPixelType<Dimension, unsigned short, Components>(); 
+  }
+  
+  else if (PixelType == "unsigned_char"){ 
+    if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and unsigned_char..." << std::endl;
+    UpdateWithDimAndPixelType<Dimension, unsigned char, Components>();
+  }
+      
+  else if(PixelType == "double"){  
+    if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and double..." << std::endl;
+    UpdateWithDimAndPixelType<Dimension, double, Components>(); 
+  }
+  else {
+    if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and float..." << std::endl;
+    UpdateWithDimAndPixelType<Dimension, float, Components>();
+  }
+}
+
+
+//-------------------------------------------------------------------
+// Update with the number of dimensions and the pixeltype
+//-------------------------------------------------------------------
+template <unsigned int Dimension, class  PixelType, unsigned int Components> 
+void 
+DVHGenericFilter::UpdateWithDimAndPixelType()
+{
+
+  // ImageTypes
+  typedef unsigned char LabelPixelType;
+  typedef itk::Image<itk::Vector<PixelType, Components>, Dimension> InputImageType;
+  typedef itk::Image<LabelPixelType, Dimension> LabelImageType;
+  
+  // Read the input
+  typedef itk::ImageFileReader<InputImageType> InputReaderType;
+  typename InputReaderType::Pointer reader = InputReaderType::New();
+  reader->SetFileName( m_InputFileName);
+  reader->Update();
+  typename InputImageType::Pointer input= reader->GetOutput();
+  
+  typedef itk::NthElementImageAdaptor<InputImageType, PixelType> InputImageAdaptorType;
+  typedef itk::Image<PixelType, Dimension> OutputImageType;
+
+  typename InputImageAdaptorType::Pointer input_adaptor = InputImageAdaptorType::New();
+  input_adaptor->SetImage(input);
+  
+  // Filter
+  typedef itk::LabelStatisticsImageFilter<InputImageAdaptorType, LabelImageType> 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<LabelPixelType, Dimension - 1> ReducedLabelImageType;
+      typedef itk::ImageFileReader<ReducedLabelImageType> LabelImageReaderType;
+      typedef itk::JoinSeriesImageFilter<ReducedLabelImageType, LabelImageType> 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<LabelImageType> LabelImageReaderType;
+      typename LabelImageReaderType::Pointer labelImageReader=LabelImageReaderType::New();
+      labelImageReader->SetFileName(m_ArgsInfo.mask_arg);
+      labelImageReader->Update();
+      labelImage= labelImageReader->GetOutput();
+
+      // Check mask sampling/size
+      if (!HaveSameSizeAndSpacing<LabelImageType, InputImageType>(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<LabelImageType> 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<LabelImageType> FilterType;
+          typename FilterType::Pointer crop = FilterType::New();
+          crop->SetInput(labelImage);
+          crop->SetCropLikeImage(input);
+          crop->Update();
+          labelImage = crop->GetOutput();                        
+          //writeImage<LabelImageType>(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 : "<<spacing[0]<<std::endl;
+  // std::cout<<"Spacing y :  "<< spacing[1]<<std::endl;
+  // std::cout<<"Spacing z :  "<< spacing[2]<<std::endl;
+  // std::cout <<"spacing_cc : "<< spacing_cc << std::endl;
+
+  // For each Label
+  typename LabelImageType::PixelType label;
+  unsigned int numberOfLabels;
+  if (m_ArgsInfo.label_given)
+    numberOfLabels=m_ArgsInfo.label_given;
+  else
+    numberOfLabels=1;
+
+  unsigned int firstComponent = 0, lastComponent = 0;
+  if (m_ArgsInfo.channel_arg == -1) {
+    firstComponent = 0; 
+    lastComponent = Components - 1;
+  }
+  else {
+    firstComponent = m_ArgsInfo.channel_arg;
+    lastComponent = m_ArgsInfo.channel_arg;
+  }
+  
+  for (unsigned int c=firstComponent; c<=lastComponent; c++) {
+    if (m_Verbose) std::cout << std::endl << "Processing channel " << c << std::endl;
+    
+    input_adaptor->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<<"--------------"<<std::endl;
+          std::cout<<"| Label:   |"<<label<<" |"<<std::endl;
+          std::cout<<"--------------"<<std::endl;
+
+        statisticsFilter->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<<statisticsFilter->GetMedian(label)<<std::endl;
+
+        typename StatisticsImageFilterType::HistogramPointer histogram =statisticsFilter->GetHistogram(label);
+
+        // Screen
+        if (m_Verbose) std::cout<<"Histogram: "<<std::endl;
+          std::cout<<"# MinBin\tMidBin\tMaxBin\tFrequency"<<std::endl;
+        for( int i =0; i <m_ArgsInfo.bins_arg; i++)
+          std::cout<<histogram->GetBinMin(0,i)<<"\t"<<histogram->GetMeasurement(i,0)<<"\t"<<histogram->GetBinMax(0,i)<<"\t"<<histogram->GetFrequency(i)<<std::endl;
+        // Add to the file
+        std::ofstream histogramFile(m_ArgsInfo.histogram_arg);
+        histogramFile<<"#Histogram: "<<std::endl;
+        histogramFile<<"#MinBin\tMidBin\tMaxBin\tFrequency"<<std::endl;
+        for( int i =0; i <m_ArgsInfo.bins_arg; i++)
+          histogramFile<<histogram->GetBinMin(0,i)<<"\t"<<histogram->GetMeasurement(i,0)<<"\t"<<histogram->GetBinMax(0,i)<<"\t"<<histogram->GetFrequency(i)<<std::endl;
+      }
+
+      // DVH
+	  if(m_ArgsInfo.dvhistogram_given)
+	  {
+          typename StatisticsImageFilterType::HistogramPointer dvhistogram = statisticsFilter->GetHistogram(label);
+
+          // Screen
+          std::cout<<"Total volume : ";
+          std::cout<<statisticsFilter->GetCount(label)<<" [No. of voxels]"<<std::endl;
+          std::cout<<"Total volume : ";
+          std::cout<<((statisticsFilter->GetCount(label))*spacing_cc)<<" [cc]"<<std::endl;
+          std::cout<<"Dose mean: ";
+          std::cout<<statisticsFilter->GetMean(label)<<" [Gy]"<<std::endl;
+          std::cout<<"Dose min: ";
+          std::cout<<statisticsFilter->GetMinimum(label)<<" [Gy]"<<std::endl;
+          std::cout<<"Dose max: ";
+          std::cout<<statisticsFilter->GetMaximum(label)<<" [Gy]"<<std::endl;
+          std::cout<<" "<<std::endl;
+          std::cout<<"#Dose[Gy] Volume_diff[No. of voxels]] Volume_cumul[No. of voxels] Volume_diff[%] Volume_cumul[%] Volume_diff[cc] Volume_cumul[cc]"<<std::endl;
+          for( int i =0; i <m_ArgsInfo.bins_arg; i++)
+          {
+            double percentDiffVolume = ((dvhistogram->GetFrequency(i))*100)/(statisticsFilter->GetCount(label));
+            double popCumulativeVolume = 0;
+            for(int j=0; j<i; j++)
+            {
+                popCumulativeVolume+=(dvhistogram->GetFrequency(j));
+            }
+            double cumulativeVolume = popCumulativeVolume + (dvhistogram->GetFrequency(i));
+            double percentCumulativeVolume =(cumulativeVolume*100)/(statisticsFilter->GetCount(label)) ;
+            double ccCumulativeVolume = (popCumulativeVolume + (dvhistogram->GetFrequency(i)))*spacing_cc;
+            std::cout<<dvhistogram->GetBinMax(0,i)<<"\t   "<<dvhistogram->GetFrequency(i)<<"\t  "<<cumulativeVolume<<"\t  "<<percentDiffVolume<<"\t "<<percentCumulativeVolume<<"\t  "<<((dvhistogram->GetFrequency(i))*spacing_cc)<<"\t  "<<ccCumulativeVolume<<std::endl;
+          }
+
+          // Add to the file
+          std::ofstream dvhistogramFile(m_ArgsInfo.dvhistogram_arg);
+          dvhistogramFile<<"Total volume : ";
+          dvhistogramFile<<statisticsFilter->GetCount(label)<<" [No. of voxels]"<<std::endl;
+          dvhistogramFile<<"Total volume : ";
+          dvhistogramFile<<((statisticsFilter->GetCount(label))*spacing_cc)<<" [cc]"<<std::endl;
+          dvhistogramFile<<"Dose mean: ";
+          dvhistogramFile<<statisticsFilter->GetMean(label)<<" [Gy]"<<std::endl;
+          dvhistogramFile<<"Dose min: ";
+          dvhistogramFile<<statisticsFilter->GetMinimum(label)<<" [Gy]"<<std::endl;
+          dvhistogramFile<<"Dose max: ";
+          dvhistogramFile<<statisticsFilter->GetMaximum(label)<<" [Gy]"<<std::endl;
+          dvhistogramFile<<"  "<<std::endl;
+          dvhistogramFile<<"#Dose[Gy] Volume_diff[No. of voxels]] Volume_cumul[No. of voxels] Volume_diff[%] Volume_cumul[%] Volume_diff[cc] Volume_cumul[cc]"<<std::endl;
+          for( int i =0; i <m_ArgsInfo.bins_arg; i++)
+          {
+             double percentDiffVolume = ((dvhistogram->GetFrequency(i))*100)/(statisticsFilter->GetCount(label));
+             double popCumulativeVolume = 0;
+             for(int j=0; j<i; j++)
+             {
+                 popCumulativeVolume+=(dvhistogram->GetFrequency(j));
+             }
+             double cumulativeVolume = popCumulativeVolume + (dvhistogram->GetFrequency(i));
+             double percentCumulativeVolume =(cumulativeVolume*100)/(statisticsFilter->GetCount(label)) ;
+             double ccCumulativeVolume = (popCumulativeVolume + (dvhistogram->GetFrequency(i)))*spacing_cc;
+             dvhistogramFile<<dvhistogram->GetBinMax(0,i)<<"\t "<<dvhistogram->GetFrequency(i)<<"\t "<<cumulativeVolume<<"\t "<<percentDiffVolume<<"\t "<<percentCumulativeVolume<<"\t "<<((dvhistogram->GetFrequency(i))*spacing_cc)<<"\t "<<ccCumulativeVolume<<std::endl;
+         }
+	  }
+    }
+  }
+
+  return;
+
+}
+
+}//end clitk
+
+#endif //#define clitkDVHGenericFilter_txx