tools/clitkHistogramImageGenericFilter.cxx

   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 clitkHistogramImageGenericFilter_cxx
  19 #define clitkHistogramImageGenericFilter_cxx
  20
  21 /* =================================================
  22  * @file   clitkHistogramImageGenericFilter.cxx
  23  * @author Thomas Baudier <thomas.baudier@creatis.insa-lyon.fr>
  24  * @date   22 dec 2015
  25  *
  26  * @brief
  27  *
  28  ===================================================*/
  29
  30 #include "clitkHistogramImageGenericFilter.h"
  31
  32 // itk include
  33 #include <itkImageToHistogramFilter.h>
  34 #include <itkStatisticsImageFilter.h>
  35
  36 #include <clitkCommon.h>
  37
  38
  39
  40 namespace clitk
  41 {
  42
  43 //--------------------------------------------------------------------
  44 HistogramImageGenericFilter::HistogramImageGenericFilter():
  45   ImageToImageGenericFilter<Self>("HistogramImage")
  46 {
  47   InitializeImageType<2>();
  48   InitializeImageType<3>();
  49   InitializeImageType<4>();
  50   mBinSize = 100;
  51 }
  52 //--------------------------------------------------------------------
  53
  54
  55 //--------------------------------------------------------------------
  56 template<unsigned int Dim>
  57 void HistogramImageGenericFilter::InitializeImageType()
  58 {
  59   ADD_DEFAULT_IMAGE_TYPES(Dim);
  60 }
  61 //--------------------------------------------------------------------
  62
  63
  64 //--------------------------------------------------------------------
  65 vtkFloatArray* HistogramImageGenericFilter::GetArrayX()
  66 {
  67   return(mArrayX);
  68 }
  69 //--------------------------------------------------------------------
  70
  71
  72 //--------------------------------------------------------------------
  73 vtkFloatArray* HistogramImageGenericFilter::GetArrayY()
  74 {
  75   return(mArrayY);
  76 }
  77 //--------------------------------------------------------------------
  78
  79
  80 //--------------------------------------------------------------------
  81 void HistogramImageGenericFilter::SetArgsInfo(const args_info_type & a)
  82 {
  83   mArgsInfo=a;
  84   if (mArgsInfo.verbose_given)
  85     SetIOVerbose(mArgsInfo.verbose_flag);
  86   if (mArgsInfo.imagetypes_given && mArgsInfo.imagetypes_flag)
  87     this->PrintAvailableImageTypes();
  88
  89   if (mArgsInfo.input_given) {
  90     SetInputFilename(mArgsInfo.input_arg);
  91   }
  92   if (mArgsInfo.output_given) {
  93     SetOutputFilename(mArgsInfo.output_arg);
  94   }
  95   if (mArgsInfo.size_given) {
  96     SetSizeBin(mArgsInfo.size_arg);
  97   }
  98 }
  99 //--------------------------------------------------------------------
 100
 101
 102 //--------------------------------------------------------------------
 103 void HistogramImageGenericFilter::SetSizeBin(const double size)
 104 {
 105   mBinSize = size;
 106 }
 107 //--------------------------------------------------------------------
 108
 109
 110 //--------------------------------------------------------------------
 111 // Update with the number of dimensions and the pixeltype
 112 //--------------------------------------------------------------------
 113 template<class InputImageType>
 114 void
 115 HistogramImageGenericFilter::UpdateWithInputImageType()
 116 {
 117
 118   // Reading input
 119   typename InputImageType::Pointer input = this->template GetInput<InputImageType>(0);
 120   typedef typename InputImageType::PixelType PixelType;
 121   typedef typename InputImageType::IndexType IndexType;
 122   typedef typename itk::Statistics::ImageToHistogramFilter<InputImageType> HistogramFilterType;
 123   typedef typename itk::StatisticsImageFilter<InputImageType> StatisticsImageFilterType;
 124
 125   //compute range
 126   typename StatisticsImageFilterType::Pointer statisticsImageFilter = StatisticsImageFilterType::New ();
 127   statisticsImageFilter->SetInput(input);
 128   statisticsImageFilter->Update();
 129   double range = statisticsImageFilter->GetMaximum() - statisticsImageFilter->GetMinimum();
 130   //compute bin number
 131   typedef typename HistogramFilterType::HistogramSizeType SizeType;
 132   SizeType binNumber(1);
 133   if (!mArgsInfo.size_given)
 134     mBinSize = std::log(range);
 135   if (mBinSize == 0)
 136     mBinSize = 1;
 137   binNumber[0] = (int)(range/mBinSize);
 138   if (binNumber[0] == 0)
 139     binNumber[0] = 1;
 140
 141   //compute histogram
 142   typename HistogramFilterType::Pointer histogramFilter = HistogramFilterType::New();
 143   histogramFilter->SetHistogramSize(binNumber);
 144   histogramFilter->SetAutoMinimumMaximum(true);
 145   typename HistogramFilterType::HistogramMeasurementVectorType lowerBound(1);
 146   typename HistogramFilterType::HistogramMeasurementVectorType upperBound(1);
 147   lowerBound[0] = statisticsImageFilter->GetMinimum();
 148   upperBound[0] = statisticsImageFilter->GetMaximum();
 149   histogramFilter->SetHistogramBinMinimum(lowerBound);
 150   histogramFilter->SetHistogramBinMaximum(upperBound);
 151   histogramFilter->SetInput(input);
 152   histogramFilter->Update();
 153
 154   mArrayX = vtkSmartPointer<vtkFloatArray>::New();
 155   mArrayY = vtkSmartPointer<vtkFloatArray>::New();
 156
 157   for(unsigned int i = 0; i < histogramFilter->GetOutput()->GetSize()[0]; ++i)
 158   {
 159     mArrayY->InsertNextTuple1(histogramFilter->GetOutput()->GetFrequency(i));
 160     mArrayX->InsertNextTuple1(statisticsImageFilter->GetMinimum() + (i+0.5)*mBinSize);
 161   }
 162 }
 163 //--------------------------------------------------------------------
 164
 165 //------------------------------------------------------------------------------
 166 void HistogramImageGenericFilter::SaveAs()
 167 {
 168   // Output
 169   std::string textFileName = GetOutputFilename();
 170   ofstream fileOpen(textFileName.c_str(), std::ofstream::trunc);
 171
 172   if(!fileOpen) {
 173       cerr << "Error during saving" << endl;
 174       return;
 175   }
 176
 177   int i(0);
 178   fileOpen << "Value represents the number of voxels around the corresponding intensity (by default the windows size around intensity is log(range))" << endl;
 179   fileOpen << "Intensity" << "\t" << "Value" << endl;
 180
 181   while (i<mArrayX->GetNumberOfTuples()) {
 182       fileOpen << mArrayX->GetTuple(i)[0] << "\t" << mArrayY->GetTuple(i)[0] << endl;
 183       ++i;
 184   }
 185
 186   fileOpen.close();
 187 }
 188 //------------------------------------------------------------------------------
 189
 190
 191 }//end clitk
 192
 193 #endif  //#define clitkHistogramImageGenericFilter_cxx