tools/clitkImageStatisticsGenericFilter.txx

   1 /*=========================================================================
   2   Program:   vv                     http://www.creatis.insa-lyon.fr/rio/vv
   3
   4   Authors belong to:
   5   - University of LYON              http://www.universite-lyon.fr/
   6   - Léon Bérard cancer center       http://www.centreleonberard.fr
   7   - CREATIS CNRS laboratory         http://www.creatis.insa-lyon.fr
   8
   9   This software is distributed WITHOUT ANY WARRANTY; without even
  10   the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  11   PURPOSE.  See the copyright notices for more information.
  12
  13   It is distributed under dual licence
  14
  15   - BSD        See included LICENSE.txt file
  16   - CeCILL-B   http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
  17 ===========================================================================**/
  18 #ifndef clitkImageStatisticsGenericFilter_txx
  19 #define clitkImageStatisticsGenericFilter_txx
  20
  21 #include "itkNthElementImageAdaptor.h"
  22
  23 /* =================================================
  24  * @file   clitkImageStatisticsGenericFilter.txx
  25  * @author
  26  * @date
  27  *
  28  * @brief
  29  *
  30  ===================================================*/
  31
  32
  33 namespace clitk
  34 {
  35
  36   //-------------------------------------------------------------------
  37   // Update with the number of dimensions
  38   //-------------------------------------------------------------------
  39   template<unsigned int Dimension, unsigned int Components>
  40   void
  41   ImageStatisticsGenericFilter::UpdateWithDim(std::string PixelType)
  42   {
  43     if (m_Verbose) std::cout << "Image was detected to be "<<Dimension<<"D and "<< PixelType<<" with " << Components << " channel(s)..."<<std::endl;
  44
  45     if(PixelType == "short"){
  46       if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and signed short..." << std::endl;
  47       UpdateWithDimAndPixelType<Dimension, signed short, Components>();
  48     }
  49     //    else if(PixelType == "unsigned_short"){
  50     //       if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and unsigned_short..." << std::endl;
  51     //       UpdateWithDimAndPixelType<Dimension, unsigned short>();
  52     //     }
  53
  54     else if (PixelType == "unsigned_char"){
  55       if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and unsigned_char..." << std::endl;
  56       UpdateWithDimAndPixelType<Dimension, unsigned char, Components>();
  57     }
  58
  59     //     else if (PixelType == "char"){
  60     //       if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and signed_char..." << std::endl;
  61     //       UpdateWithDimAndPixelType<Dimension, signed char>();
  62     //     }
  63     else {
  64       if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and float..." << std::endl;
  65       UpdateWithDimAndPixelType<Dimension, float, Components>();
  66     }
  67   }
  68
  69
  70   //-------------------------------------------------------------------
  71   // Update with the number of dimensions and the pixeltype
  72   //-------------------------------------------------------------------
  73   template <unsigned int Dimension, class  PixelType, unsigned int Components>
  74   void
  75   ImageStatisticsGenericFilter::UpdateWithDimAndPixelType()
  76   {
  77
  78     // ImageTypes
  79     typedef itk::Image<itk::Vector<PixelType, Components>, Dimension> InputImageType;
  80     typedef itk::Image<unsigned int, Dimension> LabelImageType;
  81
  82     // Read the input
  83     typedef itk::ImageFileReader<InputImageType> InputReaderType;
  84     typename InputReaderType::Pointer reader = InputReaderType::New();
  85     reader->SetFileName( m_InputFileName);
  86     reader->Update();
  87     typename InputImageType::Pointer input= reader->GetOutput();
  88
  89     typedef itk::NthElementImageAdaptor<InputImageType, PixelType> InputImageAdaptorType;
  90     typedef itk::Image<PixelType, Dimension> OutputImageType;
  91
  92     typename InputImageAdaptorType::Pointer input_adaptor = InputImageAdaptorType::New();
  93     input_adaptor->SetImage(input);
  94
  95     // Filter
  96     typedef itk::LabelStatisticsImageFilter<InputImageAdaptorType, LabelImageType> StatisticsImageFilterType;
  97     typename StatisticsImageFilterType::Pointer statisticsFilter=StatisticsImageFilterType::New();
  98     statisticsFilter->SetInput(input_adaptor);
  99
 100     // Label image
 101     typename LabelImageType::Pointer labelImage;
 102     if (m_ArgsInfo.mask_given) {
 103       typedef itk::ImageFileReader<LabelImageType> LabelImageReaderType;
 104       typename LabelImageReaderType::Pointer labelImageReader=LabelImageReaderType::New();
 105       labelImageReader->SetFileName(m_ArgsInfo.mask_arg);
 106       labelImageReader->Update();
 107       labelImage= labelImageReader->GetOutput();
 108     }
 109     else {
 110       labelImage=LabelImageType::New();
 111       labelImage->SetRegions(input->GetLargestPossibleRegion());
 112       labelImage->SetOrigin(input->GetOrigin());
 113       labelImage->SetSpacing(input->GetSpacing());
 114       labelImage->Allocate();
 115       labelImage->FillBuffer(m_ArgsInfo.label_arg[0]);
 116     }
 117     statisticsFilter->SetLabelInput(labelImage);
 118
 119     // For each Label
 120     typename LabelImageType::PixelType label;
 121     unsigned int numberOfLabels;
 122     if (m_ArgsInfo.label_given)
 123       numberOfLabels=m_ArgsInfo.label_given;
 124     else
 125       numberOfLabels=1;
 126
 127     unsigned int firstComponent = 0, lastComponent = 0;
 128     if (m_ArgsInfo.channel_arg == -1) {
 129       firstComponent = 0;
 130       lastComponent = Components - 1;
 131     }
 132     else {
 133       firstComponent = m_ArgsInfo.channel_arg;
 134       lastComponent = m_ArgsInfo.channel_arg;
 135     }
 136
 137     for (unsigned int c=firstComponent; c<=lastComponent; c++) {
 138       if (m_Verbose) std::cout << std::endl << "Processing channel " << c << std::endl;
 139
 140       input_adaptor->SelectNthElement(c);
 141       input_adaptor->Update();
 142
 143       for (unsigned int k=0; k< numberOfLabels; k++) {
 144         label=m_ArgsInfo.label_arg[k];
 145
 146         std::cout<<std::endl;
 147         if (m_Verbose) std::cout<<"-------------"<<std::endl;
 148         if (m_Verbose) std::cout<<"| Label: "<<label<<"  |"<<std::endl;
 149         if (m_Verbose) std::cout<<"-------------"<<std::endl;
 150
 151         // Histograms
 152         if (m_ArgsInfo.histogram_given) {
 153           statisticsFilter->SetUseHistograms(true);
 154           statisticsFilter->SetHistogramParameters(m_ArgsInfo.bins_arg, m_ArgsInfo.lower_arg, m_ArgsInfo.upper_arg);
 155         }
 156         statisticsFilter->Update();
 157
 158         // Output
 159         if (m_Verbose) std::cout<<"N° of pixels: ";
 160           std::cout<<statisticsFilter->GetCount(label)<<std::endl;
 161
 162         if (m_Verbose) std::cout<<"Mean: ";
 163           std::cout<<statisticsFilter->GetMean(label)<<std::endl;
 164
 165         if (m_Verbose) std::cout<<"SD: ";
 166           std::cout<<statisticsFilter->GetSigma(label)<<std::endl;
 167
 168         if (m_Verbose) std::cout<<"Variance: ";
 169           std::cout<<statisticsFilter->GetVariance(label)<<std::endl;
 170
 171         if (m_Verbose) std::cout<<"Min: ";
 172           std::cout<<statisticsFilter->GetMinimum(label)<<std::endl;
 173
 174         if (m_Verbose) std::cout<<"Max: ";
 175           std::cout<<statisticsFilter->GetMaximum(label)<<std::endl;
 176
 177         if (m_Verbose) std::cout<<"Sum: ";
 178           std::cout<<statisticsFilter->GetSum(label)<<std::endl;
 179
 180         if (m_Verbose) std::cout<<"Bounding box: ";
 181
 182         for(unsigned int i =0; i <statisticsFilter->GetBoundingBox(label).size(); i++)
 183           std::cout<<statisticsFilter->GetBoundingBox(label)[i]<<" ";
 184         std::cout<<std::endl;
 185
 186         // Histogram
 187         if (m_ArgsInfo.histogram_given)
 188         {
 189           if (m_Verbose) std::cout<<"Median: ";
 190           std::cout<<statisticsFilter->GetMedian(label)<<std::endl;
 191
 192           typename StatisticsImageFilterType::HistogramPointer histogram =statisticsFilter->GetHistogram(label);
 193
 194           // Screen
 195           if (m_Verbose) std::cout<<"Histogram: "<<std::endl;
 196             std::cout<<"# MinBin\tMidBin\tMaxBin\tFrequency"<<std::endl;
 197           for( int i =0; i <m_ArgsInfo.bins_arg; i++)
 198             std::cout<<histogram->GetBinMin(0,i)<<"\t"<<histogram->GetMeasurement(i,0)<<"\t"<<histogram->GetBinMax(0,i)<<"\t"<<histogram->GetFrequency(i)<<std::endl;
 199
 200           // Add to the file
 201           std::ofstream histogramFile(m_ArgsInfo.histogram_arg);
 202           histogramFile<<"#Histogram: "<<std::endl;
 203           histogramFile<<"#MinBin\tMidBin\tMaxBin\tFrequency"<<std::endl;
 204           for( int i =0; i <m_ArgsInfo.bins_arg; i++)
 205             histogramFile<<histogram->GetBinMin(0,i)<<"\t"<<histogram->GetMeasurement(i,0)<<"\t"<<histogram->GetBinMax(0,i)<<"\t"<<histogram->GetFrequency(i)<<std::endl;
 206         }
 207       }
 208     }
 209
 210     return;
 211
 212   }
 213
 214
 215 }//end clitk
 216
 217 #endif //#define clitkImageStatisticsGenericFilter_txx