tools/clitkDVHGenericFilter.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 clitkDVHGenericFilter_txx
  19 #define clitkDVHGenericFilter_txx
  20
  21 /* =================================================
  22  * @file   clitkDVHGenericFilter.txx
  23  * @author Agata Krason <agata.krason@creatis.insa-lyon.fr>
  24  * @date   20 November 2013
  25  *
  26  * @brief Dose volume and image histogram
  27  *
  28  ===================================================*/
  29
  30 // itk include
  31 #include "itkBinaryThresholdImageFilter.h"
  32 #include "itkMaskImageFilter.h"
  33 #include "itkMaskNegatedImageFilter.h"
  34 #include "itkNthElementImageAdaptor.h"
  35 #include "itkJoinSeriesImageFilter.h"
  36 #include "itkMinimumMaximumImageCalculator.h"
  37
  38 // clitk include
  39 #include <clitkCommon.h>
  40 #include "clitkImageCommon.h"
  41 #include "clitkDVHGenericFilter.h"
  42 #include "clitkCropLikeImageFilter.h"
  43 #include "clitkResampleImageWithOptionsFilter.h"
  44
  45 //-------------------------------------------------------------------
  46 //
  47 namespace clitk
  48 {
  49
  50 //-------------------------------------------------------------------
  51 template<unsigned int Dimension, unsigned int Components>
  52 void
  53 DVHGenericFilter::UpdateWithDim(std::string PixelType)
  54 {
  55   if (m_Verbose) std::cout << "Image was detected to be "<<Dimension<<"D and "<< PixelType<<" with " << Components << " channel(s)..."<<std::endl;
  56
  57   if(PixelType == "short"){
  58     if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and signed short..." << std::endl;
  59     UpdateWithDimAndPixelType<Dimension, signed short, Components>();
  60   }
  61   else if(PixelType == "unsigned_short"){
  62     if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and unsigned_short..." << std::endl;
  63     UpdateWithDimAndPixelType<Dimension, unsigned short, Components>();
  64   }
  65
  66   else if (PixelType == "unsigned_char"){
  67     if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and unsigned_char..." << std::endl;
  68     UpdateWithDimAndPixelType<Dimension, unsigned char, Components>();
  69   }
  70
  71   else if(PixelType == "double"){
  72     if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and double..." << std::endl;
  73     UpdateWithDimAndPixelType<Dimension, double, Components>();
  74   }
  75   else {
  76     if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and float..." << std::endl;
  77     UpdateWithDimAndPixelType<Dimension, float, Components>();
  78   }
  79 }
  80
  81
  82 //-------------------------------------------------------------------
  83 // Update with the number of dimensions and the pixeltype
  84 //-------------------------------------------------------------------
  85 template <unsigned int Dimension, class  PixelType, unsigned int Components>
  86 void
  87 DVHGenericFilter::UpdateWithDimAndPixelType()
  88 {
  89
  90   // ImageTypes
  91   typedef unsigned char LabelPixelType;
  92   typedef itk::Image<itk::Vector<PixelType, Components>, Dimension> InputImageType;
  93   typedef itk::Image<LabelPixelType, Dimension> LabelImageType;
  94
  95   // Read the input
  96   typedef itk::ImageFileReader<InputImageType> InputReaderType;
  97   typename InputReaderType::Pointer reader = InputReaderType::New();
  98   reader->SetFileName( m_InputFileName);
  99   reader->Update();
 100   typename InputImageType::Pointer input= reader->GetOutput();
 101
 102   typedef itk::NthElementImageAdaptor<InputImageType, PixelType> InputImageAdaptorType;
 103   typedef itk::Image<PixelType, Dimension> OutputImageType;
 104
 105   typename InputImageAdaptorType::Pointer input_adaptor = InputImageAdaptorType::New();
 106   input_adaptor->SetImage(input);
 107
 108   // Filter
 109   typedef itk::LabelStatisticsImageFilter<InputImageAdaptorType, LabelImageType> StatisticsImageFilterType;
 110   typename StatisticsImageFilterType::Pointer statisticsFilter=StatisticsImageFilterType::New();
 111   statisticsFilter->SetInput(input_adaptor);
 112
 113   // Label image
 114   typename LabelImageType::Pointer labelImage;
 115   if (m_ArgsInfo.mask_given) {
 116     int maskDimension, maskComponents;
 117     std::string maskPixelType;
 118     ReadImageDimensionAndPixelType(m_ArgsInfo.mask_arg, maskDimension, maskPixelType, maskComponents);
 119
 120     if (maskDimension == Dimension - 1) {
 121       // Due to a limitation of filter itk::LabelStatisticsImageFilter, InputImageType and LabelImageType
 122       // must have the same image dimension. However, we want to support label images with Dl = Di - 1,
 123       // so we need to replicate the label image as many times as the size along dimension Di.
 124       if (m_Verbose)
 125         std::cout << "Replicating label image to match input image's dimension... " << std::endl;
 126
 127       typedef itk::Image<LabelPixelType, Dimension - 1> ReducedLabelImageType;
 128       typedef itk::ImageFileReader<ReducedLabelImageType> LabelImageReaderType;
 129       typedef itk::JoinSeriesImageFilter<ReducedLabelImageType, LabelImageType> JoinImageFilterType;
 130
 131
 132       typename LabelImageReaderType::Pointer labelImageReader=LabelImageReaderType::New();
 133       labelImageReader->SetFileName(m_ArgsInfo.mask_arg);
 134       labelImageReader->Update();
 135
 136       typename JoinImageFilterType::Pointer joinFilter = JoinImageFilterType::New();
 137       typename InputImageType::SizeType size = input->GetLargestPossibleRegion().GetSize();
 138       for (unsigned int i = 0; i < size[Dimension - 1]; i++)
 139         joinFilter->PushBackInput(labelImageReader->GetOutput());
 140
 141       joinFilter->Update();
 142       labelImage = joinFilter->GetOutput();
 143     }
 144     else {
 145       typedef itk::ImageFileReader<LabelImageType> LabelImageReaderType;
 146       typename LabelImageReaderType::Pointer labelImageReader=LabelImageReaderType::New();
 147       labelImageReader->SetFileName(m_ArgsInfo.mask_arg);
 148       labelImageReader->Update();
 149       labelImage= labelImageReader->GetOutput();
 150
 151       // Check mask sampling/size
 152       if (!HaveSameSizeAndSpacing<LabelImageType, InputImageType>(labelImage, input)) {
 153         if (m_ArgsInfo.allow_resize_flag) {
 154           if (m_ArgsInfo.verbose_flag) {
 155             std::cout << "Resize mask image like input" << std::endl;
 156           }
 157
 158           typedef clitk::ResampleImageWithOptionsFilter<LabelImageType> ResamplerType;
 159           typename ResamplerType::Pointer resampler = ResamplerType::New();
 160           resampler->SetInput(labelImage);
 161           resampler->SetOutputSpacing(input->GetSpacing());
 162           resampler->SetGaussianFilteringEnabled(false);
 163           resampler->Update();
 164           labelImage = resampler->GetOutput();
 165           labelImage->GetSpacing();
 166           typedef clitk::CropLikeImageFilter<LabelImageType> FilterType;
 167           typename FilterType::Pointer crop = FilterType::New();
 168           crop->SetInput(labelImage);
 169           crop->SetCropLikeImage(input);
 170           crop->Update();
 171           labelImage = crop->GetOutput();
 172           //writeImage<LabelImageType>(labelImage, "test2.mha");
 173
 174         }
 175         else {
 176           std::cerr << "Mask image has a different size/spacing than input. Abort" << std::endl;
 177           exit(-1);
 178         }
 179       }
 180
 181     }
 182
 183   }
 184   else {
 185     labelImage=LabelImageType::New();
 186     labelImage->SetRegions(input->GetLargestPossibleRegion());
 187     labelImage->SetOrigin(input->GetOrigin());
 188     labelImage->SetSpacing(input->GetSpacing());
 189     labelImage->Allocate();
 190     labelImage->FillBuffer(m_ArgsInfo.label_arg[0]);
 191   }
 192   statisticsFilter->SetLabelInput(labelImage);
 193
 194   // Check/compute spacing
 195   const typename LabelImageType::SpacingType& spacing = input->GetSpacing();
 196   double spacing_cc = (spacing[0]*spacing[1]*spacing[2])/1000;
 197   // std::cout<<"Spacing x : "<<spacing[0]<<std::endl;
 198   // std::cout<<"Spacing y :  "<< spacing[1]<<std::endl;
 199   // std::cout<<"Spacing z :  "<< spacing[2]<<std::endl;
 200   // std::cout <<"spacing_cc : "<< spacing_cc << std::endl;
 201
 202   // For each Label
 203   typename LabelImageType::PixelType label;
 204   unsigned int numberOfLabels;
 205   if (m_ArgsInfo.label_given)
 206     numberOfLabels=m_ArgsInfo.label_given;
 207   else
 208     numberOfLabels=1;
 209
 210   unsigned int firstComponent = 0, lastComponent = 0;
 211   if (m_ArgsInfo.channel_arg == -1) {
 212     firstComponent = 0;
 213     lastComponent = Components - 1;
 214   }
 215   else {
 216     firstComponent = m_ArgsInfo.channel_arg;
 217     lastComponent = m_ArgsInfo.channel_arg;
 218   }
 219
 220   for (unsigned int c=firstComponent; c<=lastComponent; c++) {
 221     if (m_Verbose) std::cout << std::endl << "Processing channel " << c << std::endl;
 222
 223     input_adaptor->SelectNthElement(c);
 224     input_adaptor->Update();
 225
 226     for (unsigned int k=0; k< numberOfLabels; k++) {
 227       label=m_ArgsInfo.label_arg[k];
 228       // Histogram
 229       if (m_ArgsInfo.histogram_given)
 230       {
 231           std::cout<<"--------------"<<std::endl;
 232           std::cout<<"| Label:   |"<<label<<" |"<<std::endl;
 233           std::cout<<"--------------"<<std::endl;
 234
 235         statisticsFilter->SetUseHistograms(true);
 236         statisticsFilter->SetHistogramParameters(m_ArgsInfo.bins_arg, m_ArgsInfo.lower_arg, m_ArgsInfo.upper_arg);
 237       }
 238
 239       // DVHistogram
 240           if(m_ArgsInfo.dvhistogram_given)
 241           {
 242         statisticsFilter->SetUseHistograms(true);
 243         statisticsFilter->SetHistogramParameters(m_ArgsInfo.bins_arg, m_ArgsInfo.lower_arg, m_ArgsInfo.upper_arg);
 244       }
 245
 246       statisticsFilter->Update();
 247
 248       // Histogram
 249       if (m_ArgsInfo.histogram_given)
 250       {
 251           if (m_Verbose) std::cout<<"Median: ";
 252           std::cout<<statisticsFilter->GetMedian(label)<<std::endl;
 253
 254         typename StatisticsImageFilterType::HistogramPointer histogram =statisticsFilter->GetHistogram(label);
 255
 256         // Screen
 257         if (m_Verbose) std::cout<<"Histogram: "<<std::endl;
 258           std::cout<<"# MinBin\tMidBin\tMaxBin\tFrequency"<<std::endl;
 259         for( int i =0; i <m_ArgsInfo.bins_arg; i++)
 260           std::cout<<histogram->GetBinMin(0,i)<<"\t"<<histogram->GetMeasurement(i,0)<<"\t"<<histogram->GetBinMax(0,i)<<"\t"<<histogram->GetFrequency(i)<<std::endl;
 261         // Add to the file
 262         std::ofstream histogramFile(m_ArgsInfo.histogram_arg);
 263         histogramFile<<"#Histogram: "<<std::endl;
 264         histogramFile<<"#MinBin\tMidBin\tMaxBin\tFrequency"<<std::endl;
 265         for( int i =0; i <m_ArgsInfo.bins_arg; i++)
 266           histogramFile<<histogram->GetBinMin(0,i)<<"\t"<<histogram->GetMeasurement(i,0)<<"\t"<<histogram->GetBinMax(0,i)<<"\t"<<histogram->GetFrequency(i)<<std::endl;
 267       }
 268
 269       // DVH
 270           if(m_ArgsInfo.dvhistogram_given)
 271           {
 272           typename StatisticsImageFilterType::HistogramPointer dvhistogram = statisticsFilter->GetHistogram(label);
 273
 274           // Screen
 275           std::cout<<"Total volume : ";
 276           std::cout<<statisticsFilter->GetCount(label)<<" [No. of voxels]"<<std::endl;
 277           std::cout<<"Total volume : ";
 278           std::cout<<((statisticsFilter->GetCount(label))*spacing_cc)<<" [cc]"<<std::endl;
 279           std::cout<<"Dose mean: ";
 280           std::cout<<statisticsFilter->GetMean(label)<<" [Gy]"<<std::endl;
 281           std::cout<<"Dose min: ";
 282           std::cout<<statisticsFilter->GetMinimum(label)<<" [Gy]"<<std::endl;
 283           std::cout<<"Dose max: ";
 284           std::cout<<statisticsFilter->GetMaximum(label)<<" [Gy]"<<std::endl;
 285           std::cout<<" "<<std::endl;
 286           std::cout<<"#Dose[Gy] Volume_diff[No. of voxels]] Volume_cumul[No. of voxels] Volume_diff[%] Volume_cumul[%] Volume_diff[cc] Volume_cumul[cc]"<<std::endl;
 287           for( int i =0; i <m_ArgsInfo.bins_arg; i++)
 288           {
 289             double percentDiffVolume = ((dvhistogram->GetFrequency(i))*100)/(statisticsFilter->GetCount(label));
 290             double popCumulativeVolume = 0;
 291             for(int j=0; j<i; j++)
 292             {
 293                 popCumulativeVolume+=(dvhistogram->GetFrequency(j));
 294             }
 295             double cumulativeVolume = popCumulativeVolume + (dvhistogram->GetFrequency(i));
 296             double percentCumulativeVolume =(cumulativeVolume*100)/(statisticsFilter->GetCount(label)) ;
 297             double ccCumulativeVolume = (popCumulativeVolume + (dvhistogram->GetFrequency(i)))*spacing_cc;
 298             std::cout<<dvhistogram->GetBinMax(0,i)<<"\t   "<<dvhistogram->GetFrequency(i)<<"\t  "<<cumulativeVolume<<"\t  "<<percentDiffVolume<<"\t "<<percentCumulativeVolume<<"\t  "<<((dvhistogram->GetFrequency(i))*spacing_cc)<<"\t  "<<ccCumulativeVolume<<std::endl;
 299           }
 300
 301           // Add to the file
 302           std::ofstream dvhistogramFile(m_ArgsInfo.dvhistogram_arg);
 303           dvhistogramFile<<"Total volume : ";
 304           dvhistogramFile<<statisticsFilter->GetCount(label)<<" [No. of voxels]"<<std::endl;
 305           dvhistogramFile<<"Total volume : ";
 306           dvhistogramFile<<((statisticsFilter->GetCount(label))*spacing_cc)<<" [cc]"<<std::endl;
 307           dvhistogramFile<<"Dose mean: ";
 308           dvhistogramFile<<statisticsFilter->GetMean(label)<<" [Gy]"<<std::endl;
 309           dvhistogramFile<<"Dose min: ";
 310           dvhistogramFile<<statisticsFilter->GetMinimum(label)<<" [Gy]"<<std::endl;
 311           dvhistogramFile<<"Dose max: ";
 312           dvhistogramFile<<statisticsFilter->GetMaximum(label)<<" [Gy]"<<std::endl;
 313           dvhistogramFile<<"  "<<std::endl;
 314           dvhistogramFile<<"#Dose[Gy] Volume_diff[No. of voxels]] Volume_cumul[No. of voxels] Volume_diff[%] Volume_cumul[%] Volume_diff[cc] Volume_cumul[cc]"<<std::endl;
 315           for( int i =0; i <m_ArgsInfo.bins_arg; i++)
 316           {
 317              double percentDiffVolume = ((dvhistogram->GetFrequency(i))*100)/(statisticsFilter->GetCount(label));
 318              double popCumulativeVolume = 0;
 319              for(int j=0; j<i; j++)
 320              {
 321                  popCumulativeVolume+=(dvhistogram->GetFrequency(j));
 322              }
 323              double cumulativeVolume = popCumulativeVolume + (dvhistogram->GetFrequency(i));
 324              double percentCumulativeVolume =(cumulativeVolume*100)/(statisticsFilter->GetCount(label)) ;
 325              double ccCumulativeVolume = (popCumulativeVolume + (dvhistogram->GetFrequency(i)))*spacing_cc;
 326              dvhistogramFile<<dvhistogram->GetBinMax(0,i)<<"\t "<<dvhistogram->GetFrequency(i)<<"\t "<<cumulativeVolume<<"\t "<<percentDiffVolume<<"\t "<<percentCumulativeVolume<<"\t "<<((dvhistogram->GetFrequency(i))*spacing_cc)<<"\t "<<ccCumulativeVolume<<std::endl;
 327          }
 328           }
 329     }
 330   }
 331
 332   return;
 333
 334 }
 335
 336 }//end clitk
 337
 338 #endif //#define clitkDVHGenericFilter_txx