From: Romulo Pinho <pinho@lyon.fnclcc.fr>
Date: Tue, 4 Oct 2011 14:40:04 +0000 (+0200)
Subject: clitkImageStatistics supports multi-channel images
X-Git-Tag: v1.3.0~196^2~5
X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=commitdiff_plain;h=b96634936966a5da8c3d6ebb72b8954fb53b800a;p=clitk.git

clitkImageStatistics supports multi-channel images
---

diff --git a/tools/clitkImageStatistics.ggo b/tools/clitkImageStatistics.ggo
index 863dfbd..bbb127b 100644
--- 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 2982b20..23cf4df 100644
--- 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 b61beb8..6511c0e 100644
--- 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 0f52d86..cf05bb6 100644
--- 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;
-    
+
   }