tools/clitkGammaIndex3DGenericFilter.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 CLITKGAMMAINDEX3DGENERICFILTER_TXX
  19 #define CLITKGAMMAINDEX3DGENERICFILTER_TXX
  20
  21 #include "clitkImageCommon.h"
  22
  23 #include "itkApproximateSignedDistanceMapImageFilter.h"
  24 #include "itkAddImageFilter.h"
  25 #include <itkSquareImageFilter.h>
  26 #include <itkDivideImageFilter.h>
  27 #include <itkSqrtImageFilter.h>
  28 #include <itkMinimumMaximumImageCalculator.h>
  29
  30 namespace clitk
  31 {
  32
  33 //--------------------------------------------------------------------
  34 template<class args_info_type>
  35 GammaIndex3DGenericFilter<args_info_type>::GammaIndex3DGenericFilter()
  36   :ImageToImageGenericFilter<Self>("GammaIndex3DGenericFilter"), mDose(3), mDistance(3)
  37 {
  38   InitializeImageType<2>();
  39   InitializeImageType<3>();
  40   //InitializeImageType<4>();
  41 }
  42 //--------------------------------------------------------------------
  43
  44
  45 //--------------------------------------------------------------------
  46 template<class args_info_type>
  47 template<unsigned int Dim>
  48 void GammaIndex3DGenericFilter<args_info_type>::InitializeImageType()
  49 {
  50   ADD_DEFAULT_IMAGE_TYPES(Dim);
  51 }
  52 //--------------------------------------------------------------------
  53
  54
  55 //--------------------------------------------------------------------
  56 template<class args_info_type>
  57 void GammaIndex3DGenericFilter<args_info_type>::SetArgsInfo(const args_info_type & a)
  58 {
  59   mArgsInfo=a;
  60
  61   // Set value
  62   this->SetIOVerbose(mArgsInfo.verbose_flag);
  63   mDistance = mArgsInfo.distance_arg;
  64   mDose = mArgsInfo.dose_arg;
  65
  66   if (mArgsInfo.imagetypes_flag) this->PrintAvailableImageTypes();
  67
  68   if (mArgsInfo.input1_given) this->AddInputFilename(mArgsInfo.input1_arg);
  69   if (mArgsInfo.input2_given) this->AddInputFilename(mArgsInfo.input2_arg);
  70
  71   if (mArgsInfo.output_given) this->SetOutputFilename(mArgsInfo.output_arg);
  72 }
  73 //--------------------------------------------------------------------
  74
  75
  76 //--------------------------------------------------------------------
  77 template<class args_info_type>
  78 template<class ImageType>
  79 void GammaIndex3DGenericFilter<args_info_type>::UpdateWithInputImageType()
  80 {
  81   // Read input1
  82   typename ImageType::Pointer input1 = this->template GetInput<ImageType>(0);
  83
  84   // Set input image iterator
  85   typedef itk::ImageRegionIterator<ImageType> IteratorType;
  86   IteratorType it(input1, input1->GetLargestPossibleRegion());
  87
  88   // typedef input2
  89   typename ImageType::Pointer input2 = this->template GetInput<ImageType>(1);
  90
  91   //Create output
  92   typename ImageType::Pointer output = ImageType::New();
  93   output->SetRegions(input1->GetLargestPossibleRegion());
  94   output->SetOrigin(input1->GetOrigin());
  95   output->SetSpacing(input1->GetSpacing());
  96   output->Allocate();
  97
  98   int sizeTotal = input1->GetLargestPossibleRegion().GetSize(0)*input1->GetLargestPossibleRegion().GetSize(1)*input1->GetLargestPossibleRegion().GetSize(2);
  99   int sizeTemp(0);
 100   double minSpacing(std::numeric_limits<double>::max());
 101   for (unsigned int temp=0; temp < input1->GetSpacing().Size(); ++temp)
 102   {
 103     if  (input1->GetSpacing()[temp] < minSpacing)
 104       minSpacing = input1->GetSpacing()[temp];
 105   }
 106
 107
 108   while (!it.IsAtEnd())
 109   {
 110     double doseRef = it.Get();
 111
 112     //Create a region around it.Get()
 113     typename ImageType::RegionType smallRegion;
 114     typename ImageType::RegionType tempRegion;
 115     smallRegion=input2->GetLargestPossibleRegion();
 116
 117     typename ImageType::SizeType sizeRegion;
 118     sizeRegion.Fill(2*mDistance/minSpacing+1); //5 voxels de cote
 119
 120     typename ImageType::OffsetType offsetRegion;
 121     offsetRegion.Fill((sizeRegion[0]-1)/2);
 122
 123     typename ImageType::IndexType startRegion;
 124     startRegion=it.GetIndex();
 125     startRegion=startRegion-offsetRegion;
 126
 127     tempRegion.SetIndex(startRegion);
 128     tempRegion.SetSize(sizeRegion);
 129     smallRegion.Crop(tempRegion);
 130
 131     //compute the dose difference
 132     typedef itk::AddImageFilter <ImageType, ImageType, ImageType> AddImageFilterType;
 133     typename AddImageFilterType::Pointer addImageFilter = AddImageFilterType::New();
 134     addImageFilter->SetInput(input2);
 135     addImageFilter->SetConstant2(-doseRef);
 136     addImageFilter->GetOutput()->SetRequestedRegion(smallRegion);
 137     addImageFilter->Update();
 138
 139     //compute the distance
 140     typename ImageType::Pointer distanceMap = ImageType::New();
 141     distanceMap->SetRegions(smallRegion);
 142     distanceMap->SetOrigin(input1->GetOrigin());
 143     distanceMap->SetSpacing(input1->GetSpacing());
 144     distanceMap->Allocate();
 145     distanceMap->FillBuffer(0.0);
 146     //distanceMap->SetPixel(it.GetIndex(),1.0);
 147
 148     typedef  itk::ApproximateSignedDistanceMapImageFilter< ImageType, ImageType  > ApproximateSignedDistanceMapImageFilterType;
 149     typename ApproximateSignedDistanceMapImageFilterType::Pointer approximateSignedDistanceMapImageFilter = ApproximateSignedDistanceMapImageFilterType::New();
 150     approximateSignedDistanceMapImageFilter->SetInput(distanceMap);
 151     approximateSignedDistanceMapImageFilter->SetInsideValue(1.0);
 152     approximateSignedDistanceMapImageFilter->SetOutsideValue(0.0);
 153     approximateSignedDistanceMapImageFilter->GetOutput()->SetRequestedRegion(smallRegion);
 154     approximateSignedDistanceMapImageFilter->Update();
 155
 156     //Square and Divide by dose & distance criteria
 157     typedef itk::SquareImageFilter <ImageType, ImageType> SquareImageFilterType;
 158     typename SquareImageFilterType::Pointer squareImageFilter = SquareImageFilterType::New();
 159     squareImageFilter->SetInput(addImageFilter->GetOutput());
 160     squareImageFilter->GetOutput()->SetRequestedRegion(smallRegion);
 161     squareImageFilter->Update();
 162
 163     typedef itk::DivideImageFilter <ImageType, ImageType, ImageType > DivideImageFilterType;
 164     typename DivideImageFilterType::Pointer divideImageFilter = DivideImageFilterType::New ();
 165     divideImageFilter->SetInput1(squareImageFilter->GetOutput());
 166     divideImageFilter->SetConstant2(mDose*mDose);
 167     divideImageFilter->GetOutput()->SetRequestedRegion(smallRegion);
 168     divideImageFilter->Update();
 169
 170     typename SquareImageFilterType::Pointer squareImageFilter2 = SquareImageFilterType::New();
 171     squareImageFilter2->SetInput(approximateSignedDistanceMapImageFilter->GetOutput());
 172     squareImageFilter2->GetOutput()->SetRequestedRegion(smallRegion);
 173     squareImageFilter2->Update();
 174
 175     typename DivideImageFilterType::Pointer divideImageFilter2 = DivideImageFilterType::New ();
 176     divideImageFilter2->SetInput1(squareImageFilter2->GetOutput());
 177     divideImageFilter2->SetConstant2(mDistance*mDistance);
 178     divideImageFilter2->GetOutput()->SetRequestedRegion(smallRegion);
 179     divideImageFilter2->Update();
 180
 181     //Sum the 2 images and take the square root
 182     typename AddImageFilterType::Pointer addImageFilter2 = AddImageFilterType::New();
 183     addImageFilter2->SetInput1(divideImageFilter2->GetOutput());
 184     addImageFilter2->SetInput2(divideImageFilter2->GetOutput());
 185     addImageFilter2->GetOutput()->SetRequestedRegion(smallRegion);
 186     addImageFilter2->Update();
 187
 188     typedef itk::SqrtImageFilter<ImageType, ImageType> SqrtFilterType;
 189     typename SqrtFilterType::Pointer sqrtFilter = SqrtFilterType::New();
 190     sqrtFilter->SetInput(addImageFilter2->GetOutput());
 191     sqrtFilter->GetOutput()->SetRequestedRegion(smallRegion);
 192     sqrtFilter->Update();
 193
 194     //find min
 195     typedef itk::MinimumMaximumImageCalculator <ImageType> ImageCalculatorFilterType;
 196     typename ImageCalculatorFilterType::Pointer imageCalculatorFilter = ImageCalculatorFilterType::New ();
 197     imageCalculatorFilter->SetImage(sqrtFilter->GetOutput());
 198     imageCalculatorFilter->SetRegion(smallRegion);
 199     imageCalculatorFilter->ComputeMinimum();
 200
 201     //Set the value
 202     output->SetPixel(it.GetIndex(),imageCalculatorFilter->GetMinimum());
 203     ++it;
 204
 205     ++sizeTemp;
 206     std::cout << (double)sizeTemp / sizeTotal * 100.0 << std::endl;
 207   }
 208
 209   this->template SetNextOutput<ImageType>(output);
 210
 211 }
 212 //--------------------------------------------------------------------
 213
 214 } // end namespace
 215
 216 #endif  //#define CLITKGAMMAINDEX3DGENERICFILTER_TXX