From: delmon <delmon>
Date: Wed, 23 Mar 2011 13:38:26 +0000 (+0000)
Subject: pgueth:
X-Git-Tag: v1.2.0~132
X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=commitdiff_plain;h=18fd13ca69bbe5a6d91e070fef22ca9c64c92063;hp=5918ba0a0f29bb3a3ed51386be8e6bddaddefba9;p=clitk.git

pgueth:
  + Add clitkGammaIndex computation.
---

diff --git a/tools/CMakeLists.txt b/tools/CMakeLists.txt
index ac74780..b9d2441 100644
--- a/tools/CMakeLists.txt
+++ b/tools/CMakeLists.txt
@@ -96,6 +96,10 @@ IF (CLITK_BUILD_TOOLS)
     ADD_EXECUTABLE(clitkGuerreroVentilation clitkGuerreroVentilation.cxx clitkGuerreroVentilationGenericFilter.cxx ${clitkGuerreroVentilation_GGO_C})
     TARGET_LINK_LIBRARIES(clitkGuerreroVentilation clitkCommon ITKIO)
 
+    WRAP_GGO(clitkGammaIndex_GGO_C clitkGammaIndex.ggo)
+    ADD_EXECUTABLE(clitkGammaIndex clitkGammaIndex.cxx ${clitkGammaIndex_GGO_C})
+    TARGET_LINK_LIBRARIES(clitkGammaIndex ITKBasicFilters ITKCommon ITKIO)
+
     ADD_EXECUTABLE(clitkImageArithm clitkImageArithm.cxx ${clitkImageArithm_GGO_C})
     TARGET_LINK_LIBRARIES(clitkImageArithm clitkImageArithmImageLib clitkCommon ITKIO )
 
diff --git a/tools/clitkGammaIndex.cxx b/tools/clitkGammaIndex.cxx
new file mode 100644
index 0000000..d6dfaa7
--- /dev/null
+++ b/tools/clitkGammaIndex.cxx
@@ -0,0 +1,308 @@
+#include <iostream>
+using std::cout;
+using std::endl;
+
+#include <itkImageFileReader.h>
+#include <itkImageFileWriter.h>
+#include <itkMinimumMaximumImageCalculator.h>
+#include <itkShiftScaleImageFilter.h>
+#include <itkImageConstIteratorWithIndex.h>
+#include <itkChangeInformationImageFilter.h>
+#include <itkDanielssonDistanceMapImageFilter.h>
+
+#include "clitkGammaIndex_ggo.h"
+
+const unsigned int image_dim = 2;
+
+typedef itk::Image<float,image_dim> Image;
+typedef itk::ImageRegionIteratorWithIndex<Image> ImageIterator;
+typedef itk::ImageRegionConstIteratorWithIndex<Image> ImageConstIterator;
+typedef itk::ImageFileReader<Image> Reader;
+typedef itk::MinimumMaximumImageCalculator<Image> MinMaxer;
+typedef itk::ShiftScaleImageFilter<Image,Image> Normalizer;
+typedef itk::ChangeInformationImageFilter<Image> Scaler;
+
+typedef itk::Image<unsigned char,image_dim+1> ImageBin;
+typedef itk::Image<float,image_dim+1> ImageMap;
+typedef itk::DanielssonDistanceMapImageFilter<ImageBin,ImageMap> Mapper;
+
+template <typename ImageType>
+void SaveImage(const ImageType *image, const std::string &filename) {
+  typedef typename itk::ImageFileWriter<ImageType> Writer;
+  typename Writer::Pointer writer = Writer::New();
+  writer->SetFileName(filename);
+  writer->SetInput(image);
+  writer->Update();
+}
+
+ImageBin::Pointer AllocateImageBin(const Image *reference, const Image *target, unsigned int dose_size, float dose_max) {
+  const Image::RegionType reference_region = reference->GetLargestPossibleRegion();
+  const Image::RegionType target_region = target->GetLargestPossibleRegion();
+  const Image::SpacingType reference_spacing = reference->GetSpacing();
+  const Image::SpacingType target_spacing = target->GetSpacing();
+  const Image::PointType reference_origin = reference->GetOrigin();
+  const Image::PointType target_origin = target->GetOrigin();
+
+  assert(reference_region == target_region);
+  assert(reference_spacing == target_spacing);
+  assert(reference_origin == target_origin);
+
+  ImageBin::RegionType region;
+  ImageBin::SpacingType spacing;
+  ImageBin::PointType origin;
+  for (unsigned long kk=0; kk<Image::GetImageDimension(); kk++) {
+    region.SetSize(kk,reference_region.GetSize(kk));
+    region.SetIndex(kk,reference_region.GetIndex(kk));
+    spacing[kk] = reference_spacing[kk];
+    origin[kk] = reference_origin[kk];
+  }
+  const unsigned long dose_dimension = Image::GetImageDimension();
+  region.SetSize(dose_dimension,dose_size);
+  region.SetIndex(dose_dimension,0);
+  spacing[dose_dimension] = dose_max/(dose_size-1);
+  origin[dose_dimension] = 0;
+
+  ImageBin::Pointer map = ImageBin::New();
+  map->SetRegions(region);
+  map->SetSpacing(spacing);
+  map->SetOrigin(origin);
+  map->Allocate();
+  map->FillBuffer(0);
+
+  return map;
+}
+
+Image::Pointer AllocateImageGamma(const Image *target) {
+  const Image::RegionType target_region = target->GetLargestPossibleRegion();
+  const Image::SpacingType target_spacing = target->GetSpacing();
+  const Image::PointType target_origin = target->GetOrigin();
+
+  Image::Pointer gamma = Image::New();
+  gamma->SetRegions(target_region);
+  gamma->SetSpacing(target_spacing);
+  gamma->SetOrigin(target_origin);
+  gamma->Allocate();
+  gamma->FillBuffer(0);
+
+  return gamma;
+}
+
+void FillImageBin(ImageBin *image_bin, const Image *reference) {
+  ImageConstIterator iterator(reference,reference->GetLargestPossibleRegion());
+  iterator.GoToBegin();
+  while (!iterator.IsAtEnd()) {
+    Image::PixelType value = iterator.Get();
+    Image::PointType point;
+    reference->TransformIndexToPhysicalPoint(iterator.GetIndex(),point);
+
+    ImageBin::PointType point_bin;
+    for (unsigned long kk=0; kk<Image::GetImageDimension(); kk++) {
+      point_bin[kk] = point[kk];
+    }
+    point_bin[Image::GetImageDimension()] = value;
+
+    ImageBin::IndexType index_bin;
+#ifdef NDEBUG
+    image_bin->TransformPhysicalPointToIndex(point_bin,index_bin);
+#else
+    bool found = image_bin->TransformPhysicalPointToIndex(point_bin,index_bin);
+    assert(found);
+#endif
+
+    image_bin->SetPixel(index_bin,255);
+    
+    ++iterator;
+  }
+}
+
+void FillImageGamma(Image *gamma, const Image *target, const ImageMap *distance) {
+  ImageIterator gamma_iterator(gamma,gamma->GetLargestPossibleRegion());
+  ImageConstIterator target_iterator(target,target->GetLargestPossibleRegion());
+  gamma_iterator.GoToBegin();
+  target_iterator.GoToBegin();
+  while (!target_iterator.IsAtEnd()) {
+    assert(!gamma_iterator.IsAtEnd());
+
+    Image::PixelType value = target_iterator.Get();
+    Image::PointType point;
+    target->TransformIndexToPhysicalPoint(target_iterator.GetIndex(),point);
+
+    ImageMap::PointType point_map;
+    point_map.Fill(0);
+    for (unsigned long kk=0; kk<Image::GetImageDimension(); kk++) {
+      point_map[kk] = point[kk];
+    }
+    point_map[Image::GetImageDimension()] = value;
+
+    ImageMap::IndexType index_map;
+#ifdef NDEBUG
+    distance->TransformPhysicalPointToIndex(point_map,index_map);
+#else
+    bool found = distance->TransformPhysicalPointToIndex(point_map,index_map);
+    assert(found);
+#endif
+
+    gamma_iterator.Set(distance->GetPixel(index_map));
+    
+    ++gamma_iterator;
+    ++target_iterator;
+  }
+}
+
+void TuneScaler(Scaler *scaler,float space_margin) {
+  Scaler::InputImageType::PointType origin = scaler->GetInput()->GetOrigin();
+  Scaler::InputImageType::SpacingType spacing = scaler->GetInput()->GetSpacing();
+  for (unsigned int kk=0; kk<Scaler::InputImageType::GetImageDimension(); kk++) {
+    origin[kk] /= space_margin;
+    spacing[kk] /= space_margin;
+  }
+
+  scaler->SetOutputSpacing(spacing);
+  scaler->SetOutputOrigin(origin);
+  scaler->ChangeSpacingOn();
+  scaler->ChangeOriginOn();
+}
+
+Image::PixelType GetImageMaximum(const Image *image) {
+    MinMaxer::Pointer minmaxer = MinMaxer::New();
+    minmaxer->SetImage(image);
+    minmaxer->ComputeMaximum();
+    return minmaxer->GetMaximum();
+}
+
+void ComputeGammaRatio(const Image *image) {
+  ImageConstIterator iterator(image,image->GetLargestPossibleRegion());
+  iterator.GoToBegin();
+  unsigned long total = 0;
+  unsigned long over_one = 0;
+  while (!iterator.IsAtEnd()) {
+    Image::PixelType value = iterator.Get();
+
+    if (value>1) over_one++;
+    total++;
+
+    ++iterator;
+  }
+
+  cout << "total=" << total << endl;
+  cout << "over_one=" << over_one << endl;
+  cout << "ratio=" << static_cast<float>(over_one)/total << endl;
+}
+
+int main(int argc,char * argv[])
+{
+  args_info_clitkGammaIndex args_info;
+
+  if (cmdline_parser_clitkGammaIndex(argc, argv, &args_info) != 0)
+    exit(1);
+
+  if (!args_info.absolute_dose_margin_given && !args_info.relative_dose_margin_given) {
+    std::cerr << "Specify either relative or absolute dose margin" << endl;
+    exit(1);
+  }
+
+  if (args_info.isodose_number_arg <= 0) {
+    std::cerr << "Specify a valid isodose number (>0)" << endl;
+    exit(1);
+  }
+
+  bool verbose = args_info.verbose_flag;
+
+  std::string reference_filename(args_info.reference_arg);
+  std::string target_filename(args_info.target_arg);
+  std::string gamma_filename(args_info.output_arg);
+  float space_margin = args_info.spatial_margin_arg;
+  float dose_rel_margin = args_info.relative_dose_margin_arg;
+  float dose_margin = args_info.absolute_dose_margin_arg;
+  bool use_dose_margin = args_info.absolute_dose_margin_given;
+  unsigned int dose_size = args_info.isodose_number_arg;
+
+  if (verbose) {
+    cout << "reference_filename=" << reference_filename << endl;
+    cout << "target_filename=" << target_filename << endl;
+    cout << "gamma_filename=" << gamma_filename << endl;
+    cout << "dose_size=" << dose_size << endl;
+    cout << "space_margin=" << space_margin << endl;
+    if (use_dose_margin) cout << "dose_margin=" << dose_margin << endl;
+    else cout << "dose_rel_margin=" << dose_rel_margin << endl;
+  }
+
+  // load images
+  Reader::Pointer reference_reader = Reader::New();
+  Reader::Pointer target_reader = Reader::New();
+  {
+    reference_reader->SetFileName(reference_filename);
+    target_reader->SetFileName(target_filename);
+    reference_reader->Update();
+    target_reader->Update();
+  }
+
+  // intensity normalisation
+  if (!use_dose_margin) {
+    MinMaxer::PixelType reference_max = GetImageMaximum(reference_reader->GetOutput());
+    //MinMaxer::PixelType target_max = GetImageMaximum(target_reader->GetOutput());
+
+    dose_margin = reference_max*dose_rel_margin;
+
+    if (verbose) cout << "dose_margin=" << dose_margin << endl;
+  }
+
+  // scale intensity
+  Normalizer::Pointer reference_normalizer = Normalizer::New();
+  Normalizer::Pointer target_normalizer = Normalizer::New();
+  {
+    reference_normalizer->SetShift(0);
+    reference_normalizer->SetScale(1/dose_margin);
+    reference_normalizer->SetInput(reference_reader->GetOutput());
+    reference_normalizer->Update();
+
+    target_normalizer->SetShift(0);
+    target_normalizer->SetScale(1/dose_margin);
+    target_normalizer->SetInput(target_reader->GetOutput());
+    target_normalizer->Update();
+
+    //cout << "scale=" << reference_normalizer->GetScale() << "/" << target_normalizer->GetScale() << endl;
+    //cout << "shift=" << reference_normalizer->GetShift() << "/" << target_normalizer->GetShift() << endl;
+  }
+
+  // normalize space coordinates
+  Scaler::Pointer reference_scaler = Scaler::New();
+  Scaler::Pointer target_scaler = Scaler::New();
+  {
+    reference_scaler->SetInput(reference_normalizer->GetOutput());
+    TuneScaler(reference_scaler,space_margin);
+    reference_scaler->Update();
+
+    target_scaler->SetInput(target_normalizer->GetOutput());
+    TuneScaler(target_scaler,space_margin);
+    target_scaler->Update();
+
+    //SaveImage(reference_scaler->GetOutput(),"norm_reference.mhd");
+    //SaveImage(reference_scaler->GetOutput(),"norm_target.mhd");
+  }
+
+  // compute hyper surface plane
+  float reference_dose_scaled_max = GetImageMaximum(reference_scaler->GetOutput());
+  float target_dose_scaled_max = GetImageMaximum(target_scaler->GetOutput());
+  float dose_scaled_max = reference_dose_scaled_max > target_dose_scaled_max ? reference_dose_scaled_max : target_dose_scaled_max;
+  ImageBin::Pointer image_bin = AllocateImageBin(reference_scaler->GetOutput(),target_scaler->GetOutput(),dose_size,dose_scaled_max);
+  FillImageBin(image_bin,reference_scaler->GetOutput());
+  //SaveImage(image_bin.GetPointer(),"surface.mhd");
+
+  // compute distance map
+  Mapper::Pointer mapper = Mapper::New();
+  mapper->UseImageSpacingOn();
+  mapper->SetInput(image_bin);
+  mapper->Update();
+  //SaveImage(mapper->GetOutput(),"distance.mhd");
+
+  // extract gamma index from distance map
+  Image::Pointer image_gamma = AllocateImageGamma(target_normalizer->GetOutput());
+  FillImageGamma(image_gamma,target_scaler->GetOutput(),mapper->GetOutput());
+  SaveImage(image_gamma.GetPointer(),gamma_filename);
+
+  if (verbose) ComputeGammaRatio(image_gamma);
+
+  return 0;
+}
+
diff --git a/tools/clitkGammaIndex.ggo b/tools/clitkGammaIndex.ggo
new file mode 100644
index 0000000..1e6870d
--- /dev/null
+++ b/tools/clitkGammaIndex.ggo
@@ -0,0 +1,16 @@
+#File  clitkGammaIndex.ggo
+package "clitkGammaIndex"
+version "1.0"
+purpose "Compute the gamma index map for target dose map vs reference dose map. If verbose is given, the ratio of gamma>1 pixels of the total pixel in the image is computed. Absolute dose margin has priority over relative dose margin. The relative dose margin is relative to the maximum dose in the reference image. Reference and target images must share the same origin, the same spacing and the same size."
+
+option "verbose"     v 	"Verbose"		          flag	   off
+
+option "reference" i 	"Reference dose map"	  string   yes
+option "target"	   j 	"Target dose map"	  string   yes
+option "output"    o 	"Output image filename"	  string   yes
+
+option "spatial-margin"   s 	"Spatial margin [mm]"	  double   yes
+option "relative-dose-margin"   r 	"Dose margin relative to max dose in reference [%]"	  double   no
+option "absolute-dose-margin"   d 	"Absolute dose margin [Gray]"	  double   no
+option "isodose-number"   n 	"Number of isodose computed in distance map"	 int default="128"  no
+