tools/clitkGammaIndex.cxx

   1 #include <iostream>
   2 using std::cout;
   3 using std::endl;
   4
   5 #include <itkImageFileReader.h>
   6 #include <itkImageFileWriter.h>
   7 #include <itkMinimumMaximumImageCalculator.h>
   8 #include <itkShiftScaleImageFilter.h>
   9 #include <itkImageConstIteratorWithIndex.h>
  10 #include <itkChangeInformationImageFilter.h>
  11 #include <itkSignedMaurerDistanceMapImageFilter.h>
  12
  13 #include "clitkGammaIndex_ggo.h"
  14
  15 const unsigned int image_dim = 2;
  16
  17 typedef itk::Image<float,image_dim> Image;
  18 typedef itk::ImageRegionIteratorWithIndex<Image> ImageIterator;
  19 typedef itk::ImageRegionConstIteratorWithIndex<Image> ImageConstIterator;
  20 typedef itk::ImageFileReader<Image> Reader;
  21 typedef itk::MinimumMaximumImageCalculator<Image> MinMaxer;
  22 typedef itk::ShiftScaleImageFilter<Image,Image> Normalizer;
  23 typedef itk::ChangeInformationImageFilter<Image> Scaler;
  24
  25 typedef itk::Image<unsigned char,image_dim+1> ImageBin;
  26 typedef itk::Image<float,image_dim+1> ImageMap;
  27 typedef itk::SignedMaurerDistanceMapImageFilter<ImageBin,ImageMap> Mapper;
  28
  29 template <typename ImageType>
  30 void SaveImage(const ImageType *image, const std::string &filename) {
  31   typedef typename itk::ImageFileWriter<ImageType> Writer;
  32   typename Writer::Pointer writer = Writer::New();
  33   writer->SetFileName(filename);
  34   writer->SetInput(image);
  35   writer->Update();
  36 }
  37
  38 ImageBin::Pointer AllocateImageBin(const Image *reference, const Image *target, unsigned int dose_size, float dose_max) {
  39   const Image::RegionType reference_region = reference->GetLargestPossibleRegion();
  40   const Image::RegionType target_region = target->GetLargestPossibleRegion();
  41   const Image::SpacingType reference_spacing = reference->GetSpacing();
  42   const Image::SpacingType target_spacing = target->GetSpacing();
  43   const Image::PointType reference_origin = reference->GetOrigin();
  44   const Image::PointType target_origin = target->GetOrigin();
  45
  46   assert(reference_region == target_region);
  47   assert(reference_spacing == target_spacing);
  48   assert(reference_origin == target_origin);
  49
  50   ImageBin::RegionType region;
  51   ImageBin::SpacingType spacing;
  52   ImageBin::PointType origin;
  53   for (unsigned long kk=0; kk<Image::GetImageDimension(); kk++) {
  54     region.SetSize(kk,reference_region.GetSize(kk));
  55     region.SetIndex(kk,reference_region.GetIndex(kk));
  56     spacing[kk] = reference_spacing[kk];
  57     origin[kk] = reference_origin[kk];
  58   }
  59   const unsigned long dose_dimension = Image::GetImageDimension();
  60   region.SetSize(dose_dimension,dose_size);
  61   region.SetIndex(dose_dimension,0);
  62   spacing[dose_dimension] = dose_max/(dose_size-1);
  63   origin[dose_dimension] = 0;
  64
  65   ImageBin::Pointer map = ImageBin::New();
  66   map->SetRegions(region);
  67   map->SetSpacing(spacing);
  68   map->SetOrigin(origin);
  69   map->Allocate();
  70   map->FillBuffer(0);
  71
  72   return map;
  73 }
  74
  75 Image::Pointer AllocateImageGamma(const Image *target) {
  76   const Image::RegionType target_region = target->GetLargestPossibleRegion();
  77   const Image::SpacingType target_spacing = target->GetSpacing();
  78   const Image::PointType target_origin = target->GetOrigin();
  79
  80   Image::Pointer gamma = Image::New();
  81   gamma->SetRegions(target_region);
  82   gamma->SetSpacing(target_spacing);
  83   gamma->SetOrigin(target_origin);
  84   gamma->Allocate();
  85   gamma->FillBuffer(0);
  86
  87   return gamma;
  88 }
  89
  90 void FillImageBin(ImageBin *image_bin, const Image *reference) {
  91   ImageConstIterator iterator(reference,reference->GetLargestPossibleRegion());
  92   iterator.GoToBegin();
  93   while (!iterator.IsAtEnd()) {
  94     Image::PixelType value = iterator.Get();
  95     Image::PointType point;
  96     reference->TransformIndexToPhysicalPoint(iterator.GetIndex(),point);
  97
  98     ImageBin::PointType point_bin;
  99     for (unsigned long kk=0; kk<Image::GetImageDimension(); kk++) {
 100       point_bin[kk] = point[kk];
 101     }
 102     point_bin[Image::GetImageDimension()] = value;
 103
 104     ImageBin::IndexType index_bin;
 105 #ifdef NDEBUG
 106     image_bin->TransformPhysicalPointToIndex(point_bin,index_bin);
 107 #else
 108     bool found = image_bin->TransformPhysicalPointToIndex(point_bin,index_bin);
 109     assert(found);
 110 #endif
 111
 112     while (index_bin[Image::GetImageDimension()] >=0 ) {
 113       image_bin->SetPixel(index_bin,255);
 114       index_bin[Image::GetImageDimension()]--;
 115     }
 116
 117     ++iterator;
 118   }
 119 }
 120
 121 void FillImageGamma(Image *gamma, const Image *target, const ImageMap *distance) {
 122   ImageIterator gamma_iterator(gamma,gamma->GetLargestPossibleRegion());
 123   ImageConstIterator target_iterator(target,target->GetLargestPossibleRegion());
 124   gamma_iterator.GoToBegin();
 125   target_iterator.GoToBegin();
 126   while (!target_iterator.IsAtEnd()) {
 127     assert(!gamma_iterator.IsAtEnd());
 128
 129     Image::PixelType value = target_iterator.Get();
 130     Image::PointType point;
 131     target->TransformIndexToPhysicalPoint(target_iterator.GetIndex(),point);
 132
 133     ImageMap::PointType point_map;
 134     point_map.Fill(0);
 135     for (unsigned long kk=0; kk<Image::GetImageDimension(); kk++) {
 136       point_map[kk] = point[kk];
 137     }
 138     point_map[Image::GetImageDimension()] = value;
 139
 140     ImageMap::IndexType index_map;
 141 #ifdef NDEBUG
 142     distance->TransformPhysicalPointToIndex(point_map,index_map);
 143 #else
 144     bool found = distance->TransformPhysicalPointToIndex(point_map,index_map);
 145     assert(found);
 146 #endif
 147
 148     gamma_iterator.Set(fabsf(distance->GetPixel(index_map)));
 149
 150     ++gamma_iterator;
 151     ++target_iterator;
 152   }
 153 }
 154
 155 void TuneScaler(Scaler *scaler,float space_margin) {
 156   Scaler::InputImageType::PointType origin = scaler->GetInput()->GetOrigin();
 157   Scaler::InputImageType::SpacingType spacing = scaler->GetInput()->GetSpacing();
 158   for (unsigned int kk=0; kk<Scaler::InputImageType::GetImageDimension(); kk++) {
 159     origin[kk] /= space_margin;
 160     spacing[kk] /= space_margin;
 161   }
 162
 163   scaler->SetOutputSpacing(spacing);
 164   scaler->SetOutputOrigin(origin);
 165   scaler->ChangeSpacingOn();
 166   scaler->ChangeOriginOn();
 167 }
 168
 169 Image::PixelType GetImageMaximum(const Image *image) {
 170     MinMaxer::Pointer minmaxer = MinMaxer::New();
 171     minmaxer->SetImage(image);
 172     minmaxer->ComputeMaximum();
 173     return minmaxer->GetMaximum();
 174 }
 175
 176 void ComputeGammaRatio(const Image *image) {
 177   ImageConstIterator iterator(image,image->GetLargestPossibleRegion());
 178   iterator.GoToBegin();
 179   unsigned long total = 0;
 180   unsigned long over_one = 0;
 181   while (!iterator.IsAtEnd()) {
 182     Image::PixelType value = iterator.Get();
 183
 184     if (value>1) over_one++;
 185     total++;
 186
 187     ++iterator;
 188   }
 189
 190   cout << "total=" << total << endl;
 191   cout << "over_one=" << over_one << endl;
 192   cout << "ratio=" << static_cast<float>(over_one)/total << endl;
 193 }
 194
 195 int main(int argc,char * argv[])
 196 {
 197   args_info_clitkGammaIndex args_info;
 198
 199   if (cmdline_parser_clitkGammaIndex(argc, argv, &args_info) != 0)
 200     exit(1);
 201
 202   if (!args_info.absolute_dose_margin_given && !args_info.relative_dose_margin_given) {
 203     std::cerr << "Specify either relative or absolute dose margin" << endl;
 204     exit(1);
 205   }
 206
 207   if (args_info.isodose_number_arg <= 0) {
 208     std::cerr << "Specify a valid isodose number (>0)" << endl;
 209     exit(1);
 210   }
 211
 212   bool verbose = args_info.verbose_flag;
 213
 214   std::string reference_filename(args_info.reference_arg);
 215   std::string target_filename(args_info.target_arg);
 216   std::string gamma_filename(args_info.output_arg);
 217   float space_margin = args_info.spatial_margin_arg;
 218   float dose_rel_margin = args_info.relative_dose_margin_arg;
 219   float dose_margin = args_info.absolute_dose_margin_arg;
 220   bool use_dose_margin = args_info.absolute_dose_margin_given;
 221   unsigned int dose_size = args_info.isodose_number_arg;
 222
 223   if (verbose) {
 224     cout << "reference_filename=" << reference_filename << endl;
 225     cout << "target_filename=" << target_filename << endl;
 226     cout << "gamma_filename=" << gamma_filename << endl;
 227     cout << "dose_size=" << dose_size << endl;
 228     cout << "space_margin=" << space_margin << endl;
 229     if (use_dose_margin) cout << "dose_margin=" << dose_margin << endl;
 230     else cout << "dose_rel_margin=" << dose_rel_margin << endl;
 231   }
 232
 233   // load images
 234   Reader::Pointer reference_reader = Reader::New();
 235   Reader::Pointer target_reader = Reader::New();
 236   {
 237     reference_reader->SetFileName(reference_filename);
 238     target_reader->SetFileName(target_filename);
 239     reference_reader->Update();
 240     target_reader->Update();
 241   }
 242
 243   // intensity normalisation
 244   if (!use_dose_margin) {
 245     MinMaxer::PixelType reference_max = GetImageMaximum(reference_reader->GetOutput());
 246     //MinMaxer::PixelType target_max = GetImageMaximum(target_reader->GetOutput());
 247
 248     dose_margin = reference_max*dose_rel_margin;
 249
 250     if (verbose) cout << "dose_margin=" << dose_margin << endl;
 251   }
 252
 253   // scale intensity
 254   Normalizer::Pointer reference_normalizer = Normalizer::New();
 255   Normalizer::Pointer target_normalizer = Normalizer::New();
 256   {
 257     reference_normalizer->SetShift(0);
 258     reference_normalizer->SetScale(1/dose_margin);
 259     reference_normalizer->SetInput(reference_reader->GetOutput());
 260     reference_normalizer->Update();
 261
 262     target_normalizer->SetShift(0);
 263     target_normalizer->SetScale(1/dose_margin);
 264     target_normalizer->SetInput(target_reader->GetOutput());
 265     target_normalizer->Update();
 266
 267     //cout << "scale=" << reference_normalizer->GetScale() << "/" << target_normalizer->GetScale() << endl;
 268     //cout << "shift=" << reference_normalizer->GetShift() << "/" << target_normalizer->GetShift() << endl;
 269   }
 270
 271   // normalize space coordinates
 272   Scaler::Pointer reference_scaler = Scaler::New();
 273   Scaler::Pointer target_scaler = Scaler::New();
 274   {
 275     reference_scaler->SetInput(reference_normalizer->GetOutput());
 276     TuneScaler(reference_scaler,space_margin);
 277     reference_scaler->Update();
 278
 279     target_scaler->SetInput(target_normalizer->GetOutput());
 280     TuneScaler(target_scaler,space_margin);
 281     target_scaler->Update();
 282
 283     //SaveImage(reference_scaler->GetOutput(),"norm_reference.mhd");
 284     //SaveImage(reference_scaler->GetOutput(),"norm_target.mhd");
 285   }
 286
 287   // compute hyper surface plane
 288   float reference_dose_scaled_max = GetImageMaximum(reference_scaler->GetOutput());
 289   float target_dose_scaled_max = GetImageMaximum(target_scaler->GetOutput());
 290   float dose_scaled_max = reference_dose_scaled_max > target_dose_scaled_max ? reference_dose_scaled_max : target_dose_scaled_max;
 291   ImageBin::Pointer image_bin = AllocateImageBin(reference_scaler->GetOutput(),target_scaler->GetOutput(),dose_size,dose_scaled_max);
 292   FillImageBin(image_bin,reference_scaler->GetOutput());
 293   //SaveImage(image_bin.GetPointer(),"surface.mhd");
 294
 295   // compute distance map
 296   Mapper::Pointer mapper = Mapper::New();
 297   mapper->InsideIsPositiveOn();
 298   mapper->SquaredDistanceOff();
 299   mapper->UseImageSpacingOn();
 300   mapper->SetInput(image_bin);
 301   mapper->Update();
 302   //SaveImage(mapper->GetOutput(),"distance.mhd");
 303
 304   // extract gamma index from distance map
 305   Image::Pointer image_gamma = AllocateImageGamma(target_normalizer->GetOutput());
 306   FillImageGamma(image_gamma,target_scaler->GetOutput(),mapper->GetOutput());
 307   SaveImage(image_gamma.GetPointer(),gamma_filename);
 308
 309   if (verbose) ComputeGammaRatio(image_gamma);
 310
 311   return 0;
 312 }
 313