clitkImageStatistics supports multi-channel images

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..23cf4df06eb75901590f2f7b9603960b9696ad8f 100644 (file)
--- a/tools/clitkImageStatisticsGenericFilter.cxx
+++ b/tools/clitkImageStatisticsGenericFilter.cxx
@@ -50,20 +50,69 @@ 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;
+    }
 
     
     // 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..cf05bb65d642eb61680a3b49039c8619f4c0b5db 100644 (file)
--- a/tools/clitkImageStatisticsGenericFilter.txx
+++ b/tools/clitkImageStatisticsGenericFilter.txx
@@ -18,6 +18,8 @@
 #ifndef clitkImageStatisticsGenericFilter_txx
 #define clitkImageStatisticsGenericFilter_txx
 
+#include "itkNthElementImageAdaptor.h"
+
 /* =================================================
  * @file   clitkImageStatisticsGenericFilter.txx
  * @author 
@@ -34,15 +36,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 +53,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 +62,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 +70,14 @@ 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<itk::Vector<PixelType, Components>, Dimension> InputImageType;
     typedef itk::Image<unsigned int, Dimension> LabelImageType;
-    typedef itk::Image<PixelType, Dimension> OutputImageType;
     
     // Read the input
     typedef itk::ImageFileReader<InputImageType> InputReaderType;
@@ -84,31 +85,35 @@ 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();
-      }
-    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]);
-      }
+    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();
+    }
+    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 +124,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;
-    
+
   }