tools/clitkScintivolStats.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 CLITKSCINTIVOLSTATS_CXX
  19 #define CLITKSCINTIVOLSTATS_CXX
  20
  21 // clitk include
  22 #include "clitkScintivolStats_ggo.h"
  23 #include "clitkIO.h"
  24 #include "clitkImageCommon.h"
  25 #include "clitkCommon.h"
  26
  27 #include "itkImageFileReader.h"
  28 #include "itkBinaryImageToLabelMapFilter.h"
  29 #include "itkLabelMapToLabelImageFilter.h"
  30 #include "itkLabelStatisticsImageFilter.h"
  31 #include "itkStatisticsImageFilter.h"
  32 #include "itkExtractImageFilter.h"
  33
  34 #include <iostream>
  35
  36 //-------------------------------------------------------------------=
  37 int main(int argc, char * argv[])
  38 {
  39
  40   // init command line
  41   GGO(clitkScintivolStats, args_info);
  42   CLITK_INIT;
  43
  44   typedef itk::Image<double, 3> Input3DType;
  45   typedef itk::Image<double, 2> Input2DType;
  46   typedef itk::Image<char, 3> Input3DMaskType;
  47   typedef itk::Image<char, 2> Input2DMaskType;
  48   typedef itk::ImageFileReader<Input3DType> InputReader3DType;
  49   typedef itk::ImageFileReader<Input3DMaskType> InputReader3DMaskType;
  50   typedef itk::ImageFileReader<Input2DMaskType> InputReader2DMaskType;
  51   typedef itk::LabelStatisticsImageFilter<Input3DType, Input3DMaskType> LabelStatistics3DImageFilterType;
  52   typedef itk::LabelStatisticsImageFilter<Input2DType, Input2DMaskType> LabelStatistics2DImageFilterType;
  53   typedef itk::StatisticsImageFilter<Input2DType> Statistics2DImageFilterType;
  54   typedef itk::ExtractImageFilter<Input3DType, Input2DType> ExtractImageFilter;
  55
  56   //Determine nbAcquisitionDynamic1 and nbAcquisitionDynamic2 (ie. the number of slice in dynamic1 and dynamic2)
  57   //Start to open the dynamic1 and dynamic2 and get the 3rd dimension size
  58   InputReader3DType::Pointer readerDynamic1 = InputReader3DType::New();
  59   readerDynamic1->SetFileName(args_info.dynamic1_arg);
  60   readerDynamic1->Update();
  61   Input3DType::Pointer dynamic1 = readerDynamic1->GetOutput();
  62   int nbAcquisitionDynamic1 = dynamic1->GetLargestPossibleRegion().GetSize()[2];
  63
  64   InputReader3DType::Pointer readerDynamic2 = InputReader3DType::New();
  65   readerDynamic2->SetFileName(args_info.dynamic2_arg);
  66   readerDynamic2->Update();
  67   Input3DType::Pointer dynamic2 = readerDynamic2->GetOutput();
  68   int nbAcquisitionDynamic2 = dynamic2->GetLargestPossibleRegion().GetSize()[2];
  69
  70   //Open the csv file
  71   //If it's scatter or attenuation correction, just append results to the csv file
  72   //If not open it normally and write the number of dynamic acquisition for the first and the second images (given by the ggo and from the dicom tag)
  73   std::ofstream csvFile;
  74   if (args_info.append_flag)
  75     csvFile.open (args_info.output_arg, std::ios::app);
  76   else {
  77     csvFile.open (args_info.output_arg);
  78     csvFile << "1;" << nbAcquisitionDynamic1 << ";" << nbAcquisitionDynamic2 << "\n";
  79   }
  80
  81   //Read Tomo image, total and remnant masks
  82   InputReader3DType::Pointer readerTomo = InputReader3DType::New();
  83   readerTomo->SetFileName(args_info.tomo_arg);
  84   readerTomo->Update();
  85   Input3DType::Pointer tomo = readerTomo->GetOutput();
  86
  87   InputReader3DMaskType::Pointer readerTotalLiver = InputReader3DMaskType::New();
  88   readerTotalLiver->SetFileName(args_info.totalLiverMask_arg);
  89   readerTotalLiver->Update();
  90
  91   InputReader3DMaskType::Pointer readerRemnantLiver = InputReader3DMaskType::New();
  92   readerRemnantLiver->SetFileName(args_info.remnantLiverMask_arg);
  93   readerRemnantLiver->Update();
  94
  95   //Find number of counts in tomo for total Liver and remnant Liver
  96   LabelStatistics3DImageFilterType::Pointer labelStatisticsImageFilterTotalLiver = LabelStatistics3DImageFilterType::New();
  97   labelStatisticsImageFilterTotalLiver->SetLabelInput(readerTotalLiver->GetOutput());
  98   labelStatisticsImageFilterTotalLiver->SetInput(tomo);
  99   labelStatisticsImageFilterTotalLiver->Update();
 100
 101   LabelStatistics3DImageFilterType::Pointer labelStatisticsImageFilterRemnantLiver = LabelStatistics3DImageFilterType::New();
 102   labelStatisticsImageFilterRemnantLiver->SetLabelInput(readerRemnantLiver->GetOutput());
 103   labelStatisticsImageFilterRemnantLiver->SetInput(tomo);
 104   labelStatisticsImageFilterRemnantLiver->Update();
 105
 106   //Write them in the csv file
 107   csvFile << args_info.acquisitionTimeTomo_arg << ";" << labelStatisticsImageFilterTotalLiver->GetSum(1) << ";" << labelStatisticsImageFilterRemnantLiver->GetSum(1) << "\n";
 108
 109   //Read liver and heart masks for dynamic1
 110   InputReader2DMaskType::Pointer readerLiver = InputReader2DMaskType::New();
 111   readerLiver->SetFileName(args_info.liverMask_arg);
 112   readerLiver->Update();
 113
 114   InputReader2DMaskType::Pointer readerHeart = InputReader2DMaskType::New();
 115   readerHeart->SetFileName(args_info.heartMask_arg);
 116   readerHeart->Update();
 117
 118   for (unsigned int i=0; i<nbAcquisitionDynamic1; ++i) {
 119     //Compute time frame in s
 120     double timeFrame = args_info.frameDurationDynamic1_arg * i + args_info.frameDurationDynamic1_arg;
 121
 122     //Extract the Region of interest of the dynamic1
 123     ExtractImageFilter::Pointer extractSlice = ExtractImageFilter::New();
 124     Input3DType::IndexType start;
 125     start[0] = 0;
 126     start[1] = 0;
 127     start[2] = i;
 128     Input3DType::SizeType size;
 129     size[0] = dynamic1->GetLargestPossibleRegion().GetSize()[0];
 130     size[1] = dynamic1->GetLargestPossibleRegion().GetSize()[1];
 131     size[2] = 0;
 132     Input3DType::RegionType desiredRegion;
 133     desiredRegion.SetSize(size);
 134     desiredRegion.SetIndex(start);
 135     extractSlice->SetExtractionRegion(desiredRegion);
 136     extractSlice->SetInput(dynamic1);
 137 #if ITK_VERSION_MAJOR >= 4
 138     extractSlice->SetDirectionCollapseToIdentity();
 139 #endif
 140     extractSlice->Update();
 141
 142     //Find number of counts in dynamic1 slice for Liver and heart
 143     LabelStatistics2DImageFilterType::Pointer labelStatisticsImageFilterLiver = LabelStatistics2DImageFilterType::New();
 144     labelStatisticsImageFilterLiver->SetLabelInput(readerLiver->GetOutput());
 145     labelStatisticsImageFilterLiver->SetInput(extractSlice->GetOutput());
 146     labelStatisticsImageFilterLiver->SetCoordinateTolerance(0.001);
 147     labelStatisticsImageFilterLiver->SetDirectionTolerance(0.001);
 148     labelStatisticsImageFilterLiver->Update();
 149
 150     LabelStatistics2DImageFilterType::Pointer labelStatisticsImageFilterHeart = LabelStatistics2DImageFilterType::New();
 151     labelStatisticsImageFilterHeart->SetLabelInput(readerHeart->GetOutput());
 152     labelStatisticsImageFilterHeart->SetInput(extractSlice->GetOutput());
 153     labelStatisticsImageFilterHeart->SetCoordinateTolerance(0.001);
 154     labelStatisticsImageFilterHeart->SetDirectionTolerance(0.001);
 155     labelStatisticsImageFilterHeart->Update();
 156
 157     //Find number of counts in dynamic1 slice
 158     Statistics2DImageFilterType::Pointer statisticsImageFilterDynamic1 = Statistics2DImageFilterType::New();
 159     statisticsImageFilterDynamic1->SetInput(extractSlice->GetOutput());
 160     statisticsImageFilterDynamic1->Update();
 161
 162     //Write them in the csv file
 163     csvFile << timeFrame << ";" << labelStatisticsImageFilterLiver->GetSum(1) << ";" << labelStatisticsImageFilterHeart->GetSum(1) << ";" << statisticsImageFilterDynamic1->GetSum() << "\n";
 164   }
 165
 166   //Read parenchyma mask for dynamic2
 167   InputReader2DMaskType::Pointer readerParenchyma = InputReader2DMaskType::New();
 168   readerParenchyma->SetFileName(args_info.parenchymaMask_arg);
 169   readerParenchyma->Update();
 170
 171   for (unsigned int i=0; i<nbAcquisitionDynamic2; ++i) {
 172     //Compute time frame in s
 173     double timeFrame = args_info.frameDurationDynamic2_arg * i + args_info.acquisitionTimeDynamic2_arg + args_info.frameDurationDynamic2_arg;
 174
 175     //Extract the Region of interest of the dynamic1
 176     ExtractImageFilter::Pointer extractSlice = ExtractImageFilter::New();
 177     Input3DType::IndexType start;
 178     start[0] = 0;
 179     start[1] = 0;
 180     start[2] = i;
 181     Input3DType::SizeType size;
 182     size[0] = dynamic2->GetLargestPossibleRegion().GetSize()[0];
 183     size[1] = dynamic2->GetLargestPossibleRegion().GetSize()[1];
 184     size[2] = 0;
 185     Input3DType::RegionType desiredRegion;
 186     desiredRegion.SetSize(size);
 187     desiredRegion.SetIndex(start);
 188     extractSlice->SetExtractionRegion(desiredRegion);
 189     extractSlice->SetInput(dynamic2);
 190 #if ITK_VERSION_MAJOR >= 4
 191     extractSlice->SetDirectionCollapseToIdentity();
 192 #endif
 193     extractSlice->Update();
 194
 195     //Find number of counts in dynamic2 slice for parenchyma
 196     LabelStatistics2DImageFilterType::Pointer labelStatisticsImageFilterParenchyma = LabelStatistics2DImageFilterType::New();
 197     labelStatisticsImageFilterParenchyma->SetLabelInput(readerParenchyma->GetOutput());
 198     labelStatisticsImageFilterParenchyma->SetInput(extractSlice->GetOutput());
 199     labelStatisticsImageFilterParenchyma->SetCoordinateTolerance(0.001);
 200     labelStatisticsImageFilterParenchyma->SetDirectionTolerance(0.001);
 201     labelStatisticsImageFilterParenchyma->Update();
 202
 203     //Write them in the csv file
 204     csvFile << timeFrame << ";" << labelStatisticsImageFilterParenchyma->GetSum(1) << "\n";
 205   }
 206
 207   csvFile.close();
 208
 209   // this is the end my friend
 210   return 0;
 211 }
 212 //-------------------------------------------------------------------=
 213
 214 #endif /* end #define CLITKSCINTIVOLSTATS_CXX */