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