From: dsarrut <dsarrut>
Date: Fri, 25 Mar 2011 09:30:51 +0000 (+0000)
Subject: clitkGammaIndex uses vtkCellLocator
X-Git-Tag: v1.2.0~120
X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=commitdiff_plain;h=82d0db04054ac02ddbd921427ddc2edf73d475a0;p=clitk.git

clitkGammaIndex uses vtkCellLocator
---

diff --git a/tools/CMakeLists.txt b/tools/CMakeLists.txt
index b9d2441..789241b 100644
--- a/tools/CMakeLists.txt
+++ b/tools/CMakeLists.txt
@@ -98,7 +98,7 @@ IF (CLITK_BUILD_TOOLS)
 
     WRAP_GGO(clitkGammaIndex_GGO_C clitkGammaIndex.ggo)
     ADD_EXECUTABLE(clitkGammaIndex clitkGammaIndex.cxx ${clitkGammaIndex_GGO_C})
-    TARGET_LINK_LIBRARIES(clitkGammaIndex ITKBasicFilters ITKCommon ITKIO)
+    TARGET_LINK_LIBRARIES(clitkGammaIndex vtkCommon vtkFiltering vtkGraphics vtkIO vtkImaging)
 
     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
index 39cbf6f..c784122 100644
--- a/tools/clitkGammaIndex.cxx
+++ b/tools/clitkGammaIndex.cxx
@@ -1,195 +1,148 @@
+#include <vtkPoints.h>
+#include <vtkCellArray.h>
+#include <vtkPointData.h>
+#include <vtkImageData.h>
+#include <vtkMetaImageReader.h>
+#include <vtkMetaImageWriter.h>
+#include <vtkPNGReader.h>
+#include <vtkPolyData.h>
+#include <vtkCellLocator.h>
 #include <iostream>
-using std::cout;
+#include <cmath>
+#include <cassert>
 using std::endl;
-
-#include <itkImageFileReader.h>
-#include <itkImageFileWriter.h>
-#include <itkMinimumMaximumImageCalculator.h>
-#include <itkShiftScaleImageFilter.h>
-#include <itkImageConstIteratorWithIndex.h>
-#include <itkChangeInformationImageFilter.h>
-#include <itkSignedMaurerDistanceMapImageFilter.h>
+using std::cout;
 
 #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::SignedMaurerDistanceMapImageFilter<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);
+vtkImageData *loadImage(const std::string &filename) {
+  vtkImageReader2 *reader = vtkMetaImageReader::New();
+  //vtkImageReader2 *reader = vtkPNGReader::New();
+  reader->SetFileName(filename.c_str());
+  reader->Update();
+
+  vtkImageData *image = reader->GetOutput();
+  image->Register(NULL);
+  reader->Delete();
+
+  return image;
+}
+
+void saveImage(vtkImageData *image,const std::string &filename) {
+  cout << "saving " << filename << endl;
+  vtkImageWriter *writer = vtkMetaImageWriter::New();
+  writer->SetFileName(filename.c_str());
   writer->SetInput(image);
-  writer->Update();
+  writer->Write();
+  writer->Delete();
 }
 
-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];
+void insertTriangles(vtkCellArray *cells, vtkPoints *points, const vtkIdType ids[4]) {
+  double p0[3]; points->GetPoint(ids[0],p0);
+  double p1[3]; points->GetPoint(ids[1],p1);
+  double p2[3]; points->GetPoint(ids[2],p2);
+  double p3[3]; points->GetPoint(ids[3],p3);
+  //cout << "----------------------------------" << endl;
+  //cout << "p0=[" << p0[0] << "," << p0[1] << "," << p0[2] << "]" << endl;
+  //cout << "p1=[" << p1[0] << "," << p1[1] << "," << p1[2] << "]" << endl;
+  //cout << "p2=[" << p2[0] << "," << p2[1] << "," << p2[2] << "]" << endl;
+  //cout << "p3=[" << p3[0] << "," << p3[1] << "," << p3[2] << "]" << endl;
+
+  double center[] = {0,0,0};
+  for (int kk=0; kk<3; kk++) {
+    center[kk] += p0[kk];
+    center[kk] += p1[kk];
+    center[kk] += p2[kk];
+    center[kk] += p3[kk];
+    center[kk] /= 4.;
   }
-  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();
+  vtkIdType center_id = points->InsertNextPoint(center);
+  //cout << "center=[" << center[0] << "," << center[1] << "," << center[2] << "]" << endl;
+
+  cells->InsertNextCell(3);
+  cells->InsertCellPoint(ids[0]);
+  cells->InsertCellPoint(ids[1]);
+  cells->InsertCellPoint(center_id);
 
-  Image::Pointer gamma = Image::New();
-  gamma->SetRegions(target_region);
-  gamma->SetSpacing(target_spacing);
-  gamma->SetOrigin(target_origin);
-  gamma->Allocate();
-  gamma->FillBuffer(0);
+  cells->InsertNextCell(3);
+  cells->InsertCellPoint(ids[1]);
+  cells->InsertCellPoint(ids[3]);
+  cells->InsertCellPoint(center_id);
 
-  return gamma;
+  cells->InsertNextCell(3);
+  cells->InsertCellPoint(ids[3]);
+  cells->InsertCellPoint(ids[2]);
+  cells->InsertCellPoint(center_id);
+
+  cells->InsertNextCell(3);
+  cells->InsertCellPoint(ids[2]);
+  cells->InsertCellPoint(ids[0]);
+  cells->InsertCellPoint(center_id);
 }
 
-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;
+double getMaximum(vtkImageData *image) {
+  bool first = true;
+  double maximum = 0;
 
-    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
+  for (int kk=0; kk<image->GetNumberOfPoints(); kk++) {
+    double value = image->GetPointData()->GetScalars()->GetTuple1(kk);
 
-    while (index_bin[Image::GetImageDimension()] >=0 ) {
-      image_bin->SetPixel(index_bin,255);
-      index_bin[Image::GetImageDimension()]--;
+    if (first) {
+      maximum = value;
+      first = false;
+      continue;
     }
-    
-    ++iterator;
+
+    if (maximum<value) maximum = value;
   }
-}
 
-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;
+  return maximum;
+}
 
-    ImageMap::IndexType index_map;
-#ifdef NDEBUG
-    distance->TransformPhysicalPointToIndex(point_map,index_map);
+vtkPolyData *buildPlane(vtkImageData *image,double spatial_margin,double dose_margin) {
+  vtkPoints *points = vtkPoints::New();
+  for (int kk=0; kk<image->GetNumberOfPoints(); kk++) {
+    double *point = image->GetPoint(kk);
+    double value = image->GetPointData()->GetScalars()->GetTuple1(kk);
+    assert(value>=0);
+    assert(point[2]==0);
+    point[2] = value;
+
+    point[0] /= spatial_margin;
+    point[1] /= spatial_margin;
+    point[2] /= dose_margin;
+
+#ifndef NDEBUG
+    vtkIdType point_id = points->InsertNextPoint(point);
+    assert(kk==point_id);
 #else
-    bool found = distance->TransformPhysicalPointToIndex(point_map,index_map);
-    assert(found);
+    points->InsertNextPoint(point);
 #endif
-
-    gamma_iterator.Set(fabsf(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;
+  vtkCellArray *cells = vtkCellArray::New();
+  for (int kk=0; kk<image->GetNumberOfCells(); kk++) {
+    vtkCell *cell = image->GetCell(kk);
+    assert(cell->GetNumberOfPoints()==4);
+    insertTriangles(cells,points,cell->GetPointIds()->GetPointer(0));
   }
 
-  scaler->SetOutputSpacing(spacing);
-  scaler->SetOutputOrigin(origin);
-  scaler->ChangeSpacingOn();
-  scaler->ChangeOriginOn();
-}
+  vtkPolyData *data = vtkPolyData::New();
+  data->SetPoints(points);
+  data->SetPolys(cells);
+  points->Delete();
+  cells->Delete();
 
-Image::PixelType GetImageMaximum(const Image *image) {
-    MinMaxer::Pointer minmaxer = MinMaxer::New();
-    minmaxer->SetImage(image);
-    minmaxer->ComputeMaximum();
-    return minmaxer->GetMaximum();
+  return data;
 }
 
-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;
+void assert2D(vtkImageData *image) {
+#ifndef NDEBUG
+  int *extent = image->GetWholeExtent();
+  assert(extent[4]==0);
+  assert(extent[5]==0);
+#endif
 }
 
 int main(int argc,char * argv[])
@@ -204,109 +157,98 @@ int main(int argc,char * argv[])
     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;
+  double space_margin = args_info.spatial_margin_arg;
+  double dose_rel_margin = args_info.relative_dose_margin_arg;
+  double 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();
-  }
+  // load reference
+  vtkImageData *reference = loadImage(reference_filename);
+  assert2D(reference);
 
   // 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;
-
+    dose_margin = getMaximum(reference)*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;
-  }
+  // build surface
+  vtkPolyData *data = buildPlane(reference,space_margin,dose_margin);
+  reference->Delete();
+
+  vtkAbstractCellLocator *locator = vtkCellLocator::New();
+  locator->SetDataSet(data);
+  data->Delete();
+  locator->CacheCellBoundsOn();
+  locator->AutomaticOn();
+  locator->BuildLocator();
+
+  // load target
+  vtkImageData *target = loadImage(target_filename);
+  assert2D(target);
+
+  // allocate output
+  vtkImageData *output = vtkImageData::New();
+  output->SetExtent(target->GetWholeExtent());
+  output->SetOrigin(target->GetOrigin());
+  output->SetSpacing(target->GetSpacing());
+  output->SetScalarTypeToFloat();
+  output->AllocateScalars();
+
+  // fill output
+  unsigned long total = 0;
+  unsigned long over_one = 0;
+  for (int kk=0; kk<target->GetNumberOfPoints(); kk++) {
+    double *point = target->GetPoint(kk);
+    double value = target->GetPointData()->GetScalars()->GetTuple1(kk);
+    assert(value>=0);
+    assert(point[2]==0);
+    point[2] = value;
+
+    point[0] /= space_margin;
+    point[1] /= space_margin;
+    point[2] /= dose_margin;
 
-  // 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();
+    double closest_point[3] = {0,0,0};
+    vtkIdType cell_id = 0;
+    int foo = 0;
+    double squared_distance = 0;
 
-    target_scaler->SetInput(target_normalizer->GetOutput());
-    TuneScaler(target_scaler,space_margin);
-    target_scaler->Update();
+    locator->FindClosestPoint(point,closest_point,cell_id,foo,squared_distance);
+
+    double distance = sqrt(squared_distance);
+    output->GetPointData()->GetScalars()->SetTuple1(kk,distance);
+
+    if (value>1) over_one++;
+    total++;
 
-    //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->InsideIsPositiveOn();
-  mapper->SquaredDistanceOff();
-  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);
+  if (verbose) {
+    cout << "total=" << total << endl;
+    cout << "over_one=" << over_one << endl;
+    cout << "ratio=" << static_cast<float>(over_one)/total << endl;
+  }
+
+  locator->Delete();
+  target->Delete();
+
+  saveImage(output,gamma_filename);
+  output->Delete();
 
   return 0;
 }
diff --git a/tools/clitkGammaIndex.ggo b/tools/clitkGammaIndex.ggo
index 1e6870d..20c87f9 100644
--- a/tools/clitkGammaIndex.ggo
+++ b/tools/clitkGammaIndex.ggo
@@ -1,7 +1,7 @@
 #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."
+version "1.1"
+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.\n\nclitkGammaIndex -v -i ref.mhd -j disp.mhd -o gamma.mhd -s 3 -r 0.04\nthis will compute the 3mm 4% gamma index between ref and disp."
 
 option "verbose"     v 	"Verbose"		          flag	   off
 
@@ -12,5 +12,4 @@ 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