Merge branch 'master' of /home/dsarrut/clitk3.server

diff --git a/CMakeLists.txt b/CMakeLists.txt
index a8b6c0f6a837dbf1a4c13b10905c0f9b164398a2..e552d978311bade44b56c8c6f96d04d012fb6e41 100644 (file)
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -113,6 +113,7 @@ ENDIF(BUILD_TESTING)
 #=========================================================
 # Install scripts when running make install
 SET(SCRIPTS 
+  scripts/calculate_motion_amplitude.sh
   scripts/midp_common.sh
   scripts/registration.sh
   scripts/create_midP.sh

diff --git a/scripts/calculate_motion_amplitude.sh b/scripts/calculate_motion_amplitude.sh
new file mode 100755 (executable)
index 0000000..9fa0aa9
--- /dev/null
+++ b/scripts/calculate_motion_amplitude.sh
@@ -0,0 +1,52 @@
+#! /bin/bash +x 
+
+############# main 
+
+if echo $* | grep "\-h"; then
+  echo Usage: calculate_motion_amplitude.sh VECTOR_FILE RTSTRUCT_FILE RTSTRUCT_ROI_NUMBER RTSTRUCT_REF_IMAGE
+  echo "  VECTOR_FILE: mhd of the vector field from where to extract motion information (may also be 4D)"
+  echo "  RTSTRUCT_FILE: dicom with contours"
+  echo "  RTSTRUCT_ROI: number of the contour whose motion to analyze"
+  echo "  RTSTRUCT_REF_IMAGE: mhd of the reference image used in the contour delineation"
+
+  exit 0
+fi
+
+if [ $# != 4 ]; then
+  echo Invalid arguments. Type \'`basename $0` -h\' for help
+  exit -1
+fi
+
+vector_file=$1
+rtstruct_file=$2
+rtstruct_roi=$3
+rtstruct_ref_image=$4
+
+if ! cat $vector_file | grep -q "Channels = 3"; then
+  echo Vector file must have 3 channels.
+  exit -2
+fi
+
+# create ROI mask from rtstruct
+roi_mask=roi_${rtstruct_roi}.mhd
+clitkDicomRTStruct2BinaryImage -i ${rtstruct_file} -o ${roi_mask} -j ${rtstruct_ref_image} -r ${rtstruct_roi} 2> /tmp/err.txt
+if cat /tmp/err.txt | grep -q "No ROI"; then
+  echo Invalid contour number.
+  exit -3
+fi
+
+# guarantees the same sampling for roi mask and vector image
+roi_mask2=resampled_${roi_mask}
+clitkResampleImage -i ${roi_mask} -o ${roi_mask2} --like ${vector_file}
+
+# calculate motion amplitudes along the 3 image axes
+dir=( sagittal coronal axial )
+for i in 0 1 2; do
+  echo Motion along ${dir[$i]} direction
+  clitkImageStatistics -i ${vector_file} -m ${roi_mask2} -c $i --verbose 2> /dev/null | tail -n 8 | head -n 6
+  echo
+done
+
+rm `basename $roi_mask .mhd`.{mhd,raw}
+rm `basename $roi_mask2 .mhd`.{mhd,raw}
+rm /tmp/err.txt

diff --git a/tools/CMakeLists.txt b/tools/CMakeLists.txt
index d785a49f9cd5c211f235c5c2b7cd4d8d361a61bc..8bf0cbc26cf5ebf2a0a443040fac459a03bc1f1e 100644 (file)
--- a/tools/CMakeLists.txt
+++ b/tools/CMakeLists.txt
@@ -35,6 +35,11 @@ IF (CLITK_BUILD_TOOLS)
   TARGET_LINK_LIBRARIES(clitkDicom2Image clitkCommon ${ITK_LIBRARIES})
   SET(TOOLS_INSTALL ${TOOLS_INSTALL} clitkDicom2Image)
 
+  WRAP_GGO(clitkDicomWave2Text_GGO_C clitkDicomWave2Text.ggo)
+  ADD_EXECUTABLE(clitkDicomWave2Text clitkDicomWave2Text.cxx ${clitkDicomWave2Text_GGO_C})
+  TARGET_LINK_LIBRARIES(clitkDicomWave2Text clitkCommon ${ITK_LIBRARIES})
+  SET(TOOLS_INSTALL ${TOOLS_INSTALL} clitkDicomWave2Text)
+
   WRAP_GGO(clitkImageInfo_GGO_C clitkImageInfo.ggo)
   ADD_EXECUTABLE(clitkImageInfo clitkImageInfo.cxx ${clitkImageInfo_GGO_C})
   TARGET_LINK_LIBRARIES(clitkImageInfo clitkCommon ${ITK_LIBRARIES})

diff --git a/tools/clitkDicomWave2Text.cxx b/tools/clitkDicomWave2Text.cxx
new file mode 100644 (file)
index 0000000..2ceef2d
--- /dev/null
+++ b/tools/clitkDicomWave2Text.cxx
@@ -0,0 +1,104 @@
+#include "clitkDicomWave2Text.h"
+#include "clitkDicomWave2Text_ggo.h"
+
+//gdcm include
+#include "gdcmUtil.h"
+#include "gdcmFile.h"
+#include "gdcmBinEntry.h"
+#include "gdcmValEntry.h"
+#include "gdcmSeqEntry.h"
+#include "gdcmSQItem.h"
+#include "gdcmSerieHelper.h"
+
+#include <iostream>
+#include <fstream>
+
+//==========================================================================================================================
+
+int main(int argc, char * argv[]) {
+
+//-----------------------------------------------------------------------------
+// init command line
+GGO(clitkDicomWave2Text, args_info);
+//-----------------------------------------------------------------------------
+
+//-----------------------------------------------------------------------
+// opening dicom input file
+gdcm::File * mDCMFile = new gdcm::File();
+mDCMFile->SetFileName(args_info.InputFile_arg);
+mDCMFile->AddForceLoadElement(0x01e1,0x1018); //Load wave data no matter its size
+
+if (!mDCMFile->OpenFile ()) {
+       std::cerr << "Sorry, the file does not exist or does not appear to be a DICOM file. Abort." << std::endl;
+    exit(0);
+}
+mDCMFile->Load();
+std::cout << "File:   "<< args_info.InputFile_arg << "   loaded !"<< std::endl;
+
+
+//-----------------------------------------------------------------------
+// read data
+gdcm::DocEntrySet* item = mDCMFile;
+#define T short
+
+gdcm::BinEntry* entry = item->GetBinEntry(0x01e1,0x1018);
+DD(entry);
+T* buffer = reinterpret_cast<T*>(entry->GetBinArea());
+int length=item->GetEntryLength(0x01e1,0x1018)/sizeof(T);
+std::ofstream text_file(args_info.OutputFile_arg, std::ios::out | std::ios::trunc);
+if(text_file)
+{
+    std::string stat="";
+    text_file << "COMPLETE_WAVE" << '\t' << "MASK"       << '\t' << "AQUISITION_PROFIL" << '\t' << "END-INHALE" << '\t' << "END-EXHALE" << '\t' << "AQUISITION_WAVE" << '\t' << "WAVE_STATISTICS" << '\t' << "MASK"        << std::endl;
+    for (int i=0;i<length-76;i+=2)
+    {
+           if ( i < 74 )
+           {
+                   switch(i)
+                   {
+                     case 68 :
+                       stat="Total points";
+                       break;
+                     case 70 :
+                       stat="Sampling rate (Hz)";
+                       break;
+                     default :
+                       stat="";
+                       break;
+                   }
+                     
+                   if (buffer[i+75] == 0)
+                       text_file << buffer[i+74]    << '\t' << buffer[i+75] << '\t' << 0                   << '\t' << "  "         << '\t' << "  "         << '\t' << "  "              << '\t' << buffer[i]         << '\t' << buffer[i+1] << '\t' << stat << std::endl;
+                   if (buffer[i+75] == 16384)
+                       text_file << buffer[i+74]    << '\t' << buffer[i+75] << '\t' << 0                   << '\t' << buffer[i+74] << '\t' << "  "         << '\t' << "  "              << '\t' << buffer[i]         << '\t' << buffer[i+1] << '\t' << stat << std::endl;
+                   if (buffer[i+75] == 256)
+                       text_file << buffer[i+74]    << '\t' << buffer[i+75] << '\t' << 0                   << '\t' << "  "         << '\t' << buffer[i+74] << '\t' << "  "              << '\t' << buffer[i]         << '\t' << buffer[i+1] << '\t' << stat << std::endl;
+                   if (buffer[i+75] == -32768)
+                       text_file << buffer[i+74]    << '\t' << buffer[i+75] << '\t' << 1                   << '\t' << "  "         << '\t' << "  "         << '\t' << buffer[i+74]      << '\t' << buffer[i]         << '\t' << buffer[i+1] << '\t' << stat << std::endl;
+                   if (buffer[i+75] == -16384)
+                       text_file << buffer[i+74]    << '\t' << buffer[i+75] << '\t' << 1                   << '\t' << buffer[i+74] << '\t' << "  "         << '\t' << buffer[i+74]      << '\t' << buffer[i]         << '\t' << buffer[i+1] << '\t' << stat << std::endl;
+                   if (buffer[i+75] == -32512)
+                       text_file << buffer[i+74]    << '\t' << buffer[i+75] << '\t' << 1                   << '\t' << "  "         << '\t' << buffer[i+74] << '\t' << buffer[i+74]      << '\t' << buffer[i]         << '\t' << buffer[i+1] << '\t' << stat << std::endl;
+           }
+           else
+           {
+                   if (buffer[i+75] == 0)
+                       text_file << buffer[i+74]    << '\t' << buffer[i+75] << '\t' << 0                   << '\t' << "  "         << '\t' << "  "         << '\t' << "  "              << '\t' << "  "              << '\t' << "  "        << std::endl;
+                   if (buffer[i+75] == 16384)
+                       text_file << buffer[i+74]    << '\t' << buffer[i+75] << '\t' << 0                   << '\t' << buffer[i+74] << '\t' << "  "         << '\t' << "  "              << '\t' << "  "              << '\t' << "  "        << std::endl;
+                   if (buffer[i+75] == 256)
+                       text_file << buffer[i+74]    << '\t' << buffer[i+75] << '\t' << 0                   << '\t' << "  "         << '\t' << buffer[i+74] << '\t' << "  "              << '\t' << "  "              << '\t' << "  "        << std::endl;
+                   if (buffer[i+75] == -32768)
+                       text_file << buffer[i+74]    << '\t' << buffer[i+75] << '\t' << 1                   << '\t' << "  "         << '\t' << "  "         << '\t' << buffer[i+74]      << '\t' << "  "              << '\t' << "  "        << std::endl;
+                   if (buffer[i+75] == -16384)
+                       text_file << buffer[i+74]    << '\t' << buffer[i+75] << '\t' << 1                   << '\t' << buffer[i+74] << '\t' << "  "         << '\t' << buffer[i+74]      << '\t' << "  "              << '\t' << "  "        << std::endl;
+                   if (buffer[i+75] == -32512)
+                       text_file << buffer[i+74]    << '\t' << buffer[i+75] << '\t' << 1                   << '\t' << "  "         << '\t' << buffer[i+74] << '\t' << buffer[i+74]      << '\t' << "  "              << '\t' << "  "        << std::endl;
+           }
+    }
+    text_file.close();
+}
+else
+    std::cerr << "Error openning "<< args_info.OutputFile_arg << std::endl;
+       
+}

diff --git a/tools/clitkDicomWave2Text.ggo b/tools/clitkDicomWave2Text.ggo
new file mode 100644 (file)
index 0000000..26542a0
--- /dev/null
+++ b/tools/clitkDicomWave2Text.ggo
@@ -0,0 +1,8 @@
+# file DicomWave2Text
+package "DicomWave2Text"
+version "1.0"
+purpose "Extract data from a Dicom wave file to a text file"
+option "config"         -    "config file"                     string  no
+option "InputFile"      i    "Dicom inputfile name"                       string  no
+option "OutputFile"     o   "Text outputfile name"                       string  no
+option "Verbose"        v    "verbose"                         int  no

diff --git a/tools/clitkDicomWave2Text.h b/tools/clitkDicomWave2Text.h
new file mode 100644 (file)
index 0000000..fa49267
--- /dev/null
+++ b/tools/clitkDicomWave2Text.h
@@ -0,0 +1,30 @@
+#ifndef CLITKDICOMWAVE2TEXT_H
+#define CLITKDICOMWAVE2TEXT_H
+
+#include "clitkCommon.h"
+//gdcm include
+#include "gdcmUtil.h"
+#include "gdcmFile.h"
+#include "gdcmValEntry.h"
+#include "gdcmSeqEntry.h"
+#include "gdcmSQItem.h"
+#include "gdcmSerieHelper.h"
+
+
+namespace clitk {
+
+  //---------------------------------------------------------------------
+  class DicomWave2Text
+  {
+
+      public:
+      //constructor;
+      DicomWave2Text();
+      //destructor;
+      ~DicomWave2Text();
+
+  };
+  
+} // end namespace
+
+#endif

diff --git a/tools/clitkImageStatistics.ggo b/tools/clitkImageStatistics.ggo
index 863dfbd2319364036fa1e30fe9e337f441eab8bf..bbb127b6fbdd30be1543c55a2b320ebf6c2bd37e 100644 (file)
--- a/tools/clitkImageStatistics.ggo
+++ b/tools/clitkImageStatistics.ggo
@@ -1,12 +1,14 @@
 #File clitkImageStatistics.ggo
 package "clitkImageStatistics"
-version "1.0"
-purpose "Compute statistics on an image, or on part of an image specified by a mask and label(s)"
+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 image filename"            string        yes
+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        no
 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

diff --git a/tools/clitkImageStatisticsGenericFilter.cxx b/tools/clitkImageStatisticsGenericFilter.cxx
index 2982b200f04a702b2a65bd3146fc1bed4ff41051..ddfe5e0293661220055a8aabd95259273d759403 100644 (file)
--- a/tools/clitkImageStatisticsGenericFilter.cxx
+++ b/tools/clitkImageStatisticsGenericFilter.cxx
@@ -50,20 +50,79 @@ namespace clitk
   void ImageStatisticsGenericFilter::Update()
   {
     // Read the Dimension and PixelType
-    int Dimension;
+    int Dimension, Components;
     std::string PixelType;
-    ReadImageDimensionAndPixelType(m_InputFileName, Dimension, 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) UpdateWithDim<2>(PixelType);
-    else if(Dimension==3) UpdateWithDim<3>(PixelType);
-    // else if (Dimension==4)UpdateWithDim<4>(PixelType); 
-    else 
-      {
-       std::cout<<"Error, Only for 2 or 3  Dimensions!!!"<<std::endl ;
-       return;
+    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;
+    }
   }
 
 

diff --git a/tools/clitkImageStatisticsGenericFilter.h b/tools/clitkImageStatisticsGenericFilter.h
index b61beb89006e941fed751778cb81b8e72ed7f659..6511c0e0a0bcc30114f5606a5cc3a016e5a31dc3 100644 (file)
--- a/tools/clitkImageStatisticsGenericFilter.h
+++ b/tools/clitkImageStatisticsGenericFilter.h
@@ -72,7 +72,7 @@ namespace clitk
     {
       m_ArgsInfo=a;
       m_Verbose=m_ArgsInfo.verbose_flag;
-      m_InputFileName=m_ArgsInfo.input_arg;
+      m_InputFileName=m_ArgsInfo.input_arg[0];
     }
     
     
@@ -93,8 +93,8 @@ namespace clitk
     //----------------------------------------  
     // Templated members
     //----------------------------------------  
-    template <unsigned int Dimension>  void UpdateWithDim(std::string PixelType);
-    template <unsigned int Dimension, class PixelType>  void UpdateWithDimAndPixelType();
+    template <unsigned int Dimension, unsigned int Components>  void UpdateWithDim(std::string PixelType);
+    template <unsigned int Dimension, class PixelType, unsigned int Components>  void UpdateWithDimAndPixelType();
 
 
     //----------------------------------------  

diff --git a/tools/clitkImageStatisticsGenericFilter.txx b/tools/clitkImageStatisticsGenericFilter.txx
index 0f52d866973a739af9dbf0cce69da6044168f89d..d1b0d4c4abf845ce98e50054e8d259093523e37a 100644 (file)
--- a/tools/clitkImageStatisticsGenericFilter.txx
+++ b/tools/clitkImageStatisticsGenericFilter.txx
@@ -18,6 +18,9 @@
 #ifndef clitkImageStatisticsGenericFilter_txx
 #define clitkImageStatisticsGenericFilter_txx
 
+#include "itkNthElementImageAdaptor.h"
+#include "itkJoinSeriesImageFilter.h"
+
 /* =================================================
  * @file   clitkImageStatisticsGenericFilter.txx
  * @author 
@@ -26,6 +29,7 @@
  * @brief 
  * 
  ===================================================*/
+#include "clitkImageStatisticsGenericFilter.h"
 
 
 namespace clitk
@@ -34,15 +38,15 @@ namespace clitk
   //-------------------------------------------------------------------
   // Update with the number of dimensions
   //-------------------------------------------------------------------
-  template<unsigned int Dimension>
+  template<unsigned int Dimension, unsigned int Components>
   void 
   ImageStatisticsGenericFilter::UpdateWithDim(std::string PixelType)
   {
-    if (m_Verbose) std::cout << "Image was detected to be "<<Dimension<<"D and "<< PixelType<<"..."<<std::endl;
+    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>(); 
+      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;
@@ -51,7 +55,7 @@ namespace clitk
     
     else if (PixelType == "unsigned_char"){ 
       if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and unsigned_char..." << std::endl;
-      UpdateWithDimAndPixelType<Dimension, unsigned char>();
+      UpdateWithDimAndPixelType<Dimension, unsigned char, Components>();
     }
     
     //     else if (PixelType == "char"){ 
@@ -60,7 +64,7 @@ namespace clitk
     //     }
     else {
       if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and float..." << std::endl;
-      UpdateWithDimAndPixelType<Dimension, float>();
+      UpdateWithDimAndPixelType<Dimension, float, Components>();
     }
   }
 
@@ -68,15 +72,15 @@ namespace clitk
   //-------------------------------------------------------------------
   // Update with the number of dimensions and the pixeltype
   //-------------------------------------------------------------------
-  template <unsigned int Dimension, class  PixelType> 
+  template <unsigned int Dimension, class  PixelType, unsigned int Components> 
   void 
   ImageStatisticsGenericFilter::UpdateWithDimAndPixelType()
   {
 
     // ImageTypes
-    typedef itk::Image<PixelType, Dimension> InputImageType;
-    typedef itk::Image<unsigned int, Dimension> LabelImageType;
-    typedef itk::Image<PixelType, Dimension> OutputImageType;
+    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;
@@ -84,31 +88,64 @@ namespace clitk
     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<InputImageType, LabelImageType> StatisticsImageFilterType;
+    typedef itk::LabelStatisticsImageFilter<InputImageAdaptorType, LabelImageType> StatisticsImageFilterType;
     typename StatisticsImageFilterType::Pointer statisticsFilter=StatisticsImageFilterType::New();
-    statisticsFilter->SetInput(input);
+    statisticsFilter->SetInput(input_adaptor);
 
     // Label image
     typename LabelImageType::Pointer labelImage;
-    if (m_ArgsInfo.mask_given)
-      {
-       typedef itk::ImageFileReader<LabelImageType> LabelImageReaderType;
-       typename LabelImageReaderType::Pointer labelImageReader=LabelImageReaderType::New();
-       labelImageReader->SetFileName(m_ArgsInfo.mask_arg);
-       labelImageReader->Update();
-       labelImage= labelImageReader->GetOutput();
+    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
-      { 
-       labelImage=LabelImageType::New();
-       labelImage->SetRegions(input->GetLargestPossibleRegion());
-       labelImage->SetOrigin(input->GetOrigin());
-       labelImage->SetSpacing(input->GetSpacing());
-       labelImage->Allocate();
-       labelImage->FillBuffer(m_ArgsInfo.label_arg[0]);
+      else {
+        typedef itk::ImageFileReader<LabelImageType> LabelImageReaderType;
+        typename LabelImageReaderType::Pointer labelImageReader=LabelImageReaderType::New();
+        labelImageReader->SetFileName(m_ArgsInfo.mask_arg);
+        labelImageReader->Update();
+        labelImage= labelImageReader->GetOutput();
       }
+
+    }
+    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);
 
     // For each Label
@@ -119,77 +156,91 @@ namespace clitk
     else
       numberOfLabels=1;
 
-    for (unsigned int k=0; k< numberOfLabels; k++)
-      {
-       label=m_ArgsInfo.label_arg[k];
-
-       std::cout<<std::endl;
-       if (m_Verbose) std::cout<<"-------------"<<std::endl;
-       if (m_Verbose) std::cout<<"| Label: "<<label<<"  |"<<std::endl;
-       if (m_Verbose) std::cout<<"-------------"<<std::endl;
-       
-       // Histograms
-       if (m_ArgsInfo.histogram_given)
-         {
-           statisticsFilter->SetUseHistograms(true);
-           statisticsFilter->SetHistogramParameters(m_ArgsInfo.bins_arg, m_ArgsInfo.lower_arg, m_ArgsInfo.upper_arg);
-         }
-       statisticsFilter->Update();
-
-
-       // Output
-       if (m_Verbose) std::cout<<"N° of pixels: ";
-       std::cout<<statisticsFilter->GetCount(label)<<std::endl;
-
-       if (m_Verbose) std::cout<<"Mean: ";
-       std::cout<<statisticsFilter->GetMean(label)<<std::endl;
+    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;
+    }
     
-       if (m_Verbose) std::cout<<"SD: ";
-       std::cout<<statisticsFilter->GetSigma(label)<<std::endl;
-
-       if (m_Verbose) std::cout<<"Variance: ";
-       std::cout<<statisticsFilter->GetVariance(label)<<std::endl;
-
-       if (m_Verbose) std::cout<<"Min: ";
-       std::cout<<statisticsFilter->GetMinimum(label)<<std::endl;
-
-       if (m_Verbose) std::cout<<"Max: ";
-       std::cout<<statisticsFilter->GetMaximum(label)<<std::endl;
-
-       if (m_Verbose) std::cout<<"Sum: ";
-       std::cout<<statisticsFilter->GetSum(label)<<std::endl;
-
-       if (m_Verbose) std::cout<<"Bounding box: ";
-       for(unsigned int i =0; i <statisticsFilter->GetBoundingBox(label).size(); i++)
-         std::cout<<statisticsFilter->GetBoundingBox(label)[i]<<" ";
-       std::cout<<std::endl;
-
-
-       // 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;
-         }
+    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];
+
+        std::cout<<std::endl;
+        if (m_Verbose) std::cout<<"-------------"<<std::endl;
+        if (m_Verbose) std::cout<<"| Label: "<<label<<"  |"<<std::endl;
+        if (m_Verbose) std::cout<<"-------------"<<std::endl;
+
+        // Histograms
+        if (m_ArgsInfo.histogram_given) {
+          statisticsFilter->SetUseHistograms(true);
+          statisticsFilter->SetHistogramParameters(m_ArgsInfo.bins_arg, m_ArgsInfo.lower_arg, m_ArgsInfo.upper_arg);
+        }
+        statisticsFilter->Update();
+
+        // Output
+        if (m_Verbose) std::cout<<"N° of pixels: ";
+          std::cout<<statisticsFilter->GetCount(label)<<std::endl;
+
+        if (m_Verbose) std::cout<<"Mean: ";
+          std::cout<<statisticsFilter->GetMean(label)<<std::endl;
+
+        if (m_Verbose) std::cout<<"SD: ";
+          std::cout<<statisticsFilter->GetSigma(label)<<std::endl;
+
+        if (m_Verbose) std::cout<<"Variance: ";
+          std::cout<<statisticsFilter->GetVariance(label)<<std::endl;
+
+        if (m_Verbose) std::cout<<"Min: ";
+          std::cout<<statisticsFilter->GetMinimum(label)<<std::endl;
+
+        if (m_Verbose) std::cout<<"Max: ";
+          std::cout<<statisticsFilter->GetMaximum(label)<<std::endl;
+
+        if (m_Verbose) std::cout<<"Sum: ";
+          std::cout<<statisticsFilter->GetSum(label)<<std::endl;
+
+        if (m_Verbose) std::cout<<"Bounding box: ";
+        
+        for(unsigned int i =0; i <statisticsFilter->GetBoundingBox(label).size(); i++)
+          std::cout<<statisticsFilter->GetBoundingBox(label)[i]<<" ";
+        std::cout<<std::endl;
+
+        // 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;
+        }
       }
-    
+    }
+
     return;
-    
+
   }