X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=tools%2FclitkGammaIndex.cxx;h=be0f8072218330f4b82a99ddc7c1b7f50419698e;hb=f2547abd3ebf99109ffd9ad4a7bc7b2e73ea91b5;hp=d6dfaa78ce0905fc0a5e16ff1f0008b1b739a747;hpb=18fd13ca69bbe5a6d91e070fef22ca9c64c92063;p=clitk.git

diff --git a/tools/clitkGammaIndex.cxx b/tools/clitkGammaIndex.cxx
index d6dfaa7..be0f807 100644
--- a/tools/clitkGammaIndex.cxx
+++ b/tools/clitkGammaIndex.cxx
@@ -1,308 +1,306 @@
+/*=========================================================================
+  Program:   vv                     http://www.creatis.insa-lyon.fr/rio/vv
+
+  Authors belong to:
+  - University of LYON              http://www.universite-lyon.fr/
+  - Léon Bérard cancer center       http://www.centreleonberard.fr
+  - CREATIS CNRS laboratory         http://www.creatis.insa-lyon.fr
+
+  This software is distributed WITHOUT ANY WARRANTY; without even
+  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+  PURPOSE.  See the copyright notices for more information.
+
+  It is distributed under dual licence
+
+  - BSD        See included LICENSE.txt file
+  - CeCILL-B   http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
+===========================================================================*/
+#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 <itkDanielssonDistanceMapImageFilter.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::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();
-}
+#include <clitkIO.h>
+#include <vvImage.h>
+#include <vvImageReader.h>
+#include <vvImageWriter.h>
 
-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;
-}
+#include <itkImage.h>
+#include <itkImageRegionIterator.h>
+typedef itk::Image<double> OutputImageType;
+typedef itk::ImageRegionIterator<OutputImageType> OutputImageIterator;
 
-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();
+vtkImageData* loadImage(const std::string& filename)
+{
+    vvImageReader::Pointer reader = vvImageReader::New();
+    reader->SetInputFilename(filename);
+    reader->Update();
+    vvImage::Pointer vvimage = reader->GetOutput();
+    if (!vvimage) { cerr << "can't load " << filename << endl; return NULL; }
+
+    vtkImageData *image = vtkImageData::New();
+    image->DeepCopy(vvimage->GetFirstVTKImageData());
+    return image;
+}
 
-  Image::Pointer gamma = Image::New();
-  gamma->SetRegions(target_region);
-  gamma->SetSpacing(target_spacing);
-  gamma->SetOrigin(target_origin);
-  gamma->Allocate();
-  gamma->FillBuffer(0);
+void saveImage(OutputImageType* image,const std::string &filename) {
+    vvImage::Pointer vvimage = vvImage::New();
+    vvimage->AddItkImage(image);
 
-  return gamma;
+    vvImageWriter::Pointer writer = vvImageWriter::New();
+    writer->SetOutputFileName(filename.c_str());
+    writer->SetInput(vvimage);
+    writer->Update();
 }
 
-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];
+
+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.;
     }
-    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
+    vtkIdType center_id = points->InsertNextPoint(center);
+    //cout << "center=[" << center[0] << "," << center[1] << "," << center[2] << "]" << endl;
 
-    image_bin->SetPixel(index_bin,255);
-    
-    ++iterator;
-  }
-}
+    cells->InsertNextCell(3);
+    cells->InsertCellPoint(ids[0]);
+    cells->InsertCellPoint(ids[1]);
+    cells->InsertCellPoint(center_id);
 
-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;
+    cells->InsertNextCell(3);
+    cells->InsertCellPoint(ids[1]);
+    cells->InsertCellPoint(ids[3]);
+    cells->InsertCellPoint(center_id);
 
-    ImageMap::IndexType index_map;
-#ifdef NDEBUG
-    distance->TransformPhysicalPointToIndex(point_map,index_map);
-#else
-    bool found = distance->TransformPhysicalPointToIndex(point_map,index_map);
-    assert(found);
-#endif
+    cells->InsertNextCell(3);
+    cells->InsertCellPoint(ids[3]);
+    cells->InsertCellPoint(ids[2]);
+    cells->InsertCellPoint(center_id);
 
-    gamma_iterator.Set(distance->GetPixel(index_map));
-    
-    ++gamma_iterator;
-    ++target_iterator;
-  }
+    cells->InsertNextCell(3);
+    cells->InsertCellPoint(ids[2]);
+    cells->InsertCellPoint(ids[0]);
+    cells->InsertCellPoint(center_id);
 }
 
-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();
+void insertLine(vtkCellArray *cells, vtkPoints *points, const vtkIdType ids[2]) {
+    cells->InsertNextCell(2);
+    cells->InsertCellPoint(ids[0]);
+    cells->InsertCellPoint(ids[1]);
 }
 
-Image::PixelType GetImageMaximum(const Image *image) {
-    MinMaxer::Pointer minmaxer = MinMaxer::New();
-    minmaxer->SetImage(image);
-    minmaxer->ComputeMaximum();
-    return minmaxer->GetMaximum();
+double getMaximum(vtkImageData *image) {
+    bool first = true;
+    double maximum = 0;
+
+    vtkPointData* point_data = image->GetPointData();
+    assert(point_data);
+    vtkDataArray* scalars = point_data->GetScalars();
+    assert(scalars);
+
+    for (int kk=0; kk<image->GetNumberOfPoints(); kk++) {
+	double value = scalars->GetTuple1(kk);
+
+	if (first) {
+	    maximum = value;
+	    first = false;
+	    continue;
+	}
+
+	if (maximum<value) maximum = value;
+    }
+
+    return maximum;
 }
 
-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++;
+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
+	points->InsertNextPoint(point);
+#endif
+    }
+
+    vtkCellArray *cells = vtkCellArray::New();
+    for (int kk=0; kk<image->GetNumberOfCells(); kk++) {
+	vtkCell *cell = image->GetCell(kk);
+
+	if (cell->GetNumberOfPoints()==4) {
+	    insertTriangles(cells,points,cell->GetPointIds()->GetPointer(0));
+	    continue;
+	}
+
+	if (cell->GetNumberOfPoints()==2) {
+	    insertLine(cells,points,cell->GetPointIds()->GetPointer(0));
+	    continue;
+	}
+
+	assert(false);
+    }
 
-    ++iterator;
-  }
+    vtkPolyData *data = vtkPolyData::New();
+    data->SetPoints(points);
+    data->SetPolys(cells);
+    points->Delete();
+    cells->Delete();
 
-  cout << "total=" << total << endl;
-  cout << "over_one=" << over_one << endl;
-  cout << "ratio=" << static_cast<float>(over_one)/total << endl;
+    return data;
 }
 
 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;
+    clitk::RegisterClitkFactories();
+
+    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);
+    }
+
+    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);
+    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;
+
+    if (verbose) {
+	cout << "reference_filename=" << reference_filename << endl;
+	cout << "target_filename=" << target_filename << endl;
+	cout << "gamma_filename=" << gamma_filename << 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 reference
+    vtkImageData* reference = loadImage(reference_filename);
+    assert(reference);
+
+    // intensity normalisation
+    if (!use_dose_margin) {
+	dose_margin = getMaximum(reference)*dose_rel_margin;
+	if (verbose) cout << "dose_margin=" << dose_margin << 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);
+    assert(target);
+
+    // allocate output
+    OutputImageType::Pointer output = OutputImageType::New();
+    {
+	OutputImageType::SizeType::SizeValueType output_array_size[2];
+	output_array_size[0] = target->GetDimensions()[0];
+	output_array_size[1] = target->GetDimensions()[1];
+	OutputImageType::SizeType output_size;
+	output_size.SetSize(output_array_size);
+	output->SetRegions(OutputImageType::RegionType(output_size));
+	output->SetOrigin(target->GetOrigin());
+	output->SetSpacing(target->GetSpacing());
+	output->Allocate();
+    }
+
+    // fill output
+    unsigned long  kk = 0;
+    unsigned long over_one = 0;
+    OutputImageIterator iter(output,output->GetLargestPossibleRegion());
+    iter.GoToBegin();
+    while (!iter.IsAtEnd()) {
+	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;
+
+	double closest_point[3] = {0,0,0};
+	vtkIdType cell_id = 0;
+	int foo = 0;
+	double squared_distance = 0;
+
+	locator->FindClosestPoint(point,closest_point,cell_id,foo,squared_distance);
+	double distance = sqrt(squared_distance);
+	iter.Set(distance);
+
+	if (distance>1) over_one++;
+	kk++;
+	++iter;
+    }
+
+    if (verbose) {
+	cout << "total=" << kk << endl;
+	cout << "over_one=" << over_one << endl;
+	cout << "ratio=" << static_cast<float>(over_one)/kk << endl;
+    }
+
+    locator->Delete();
+    target->Delete();
+
+    saveImage(output,gamma_filename);
+
+    return 0;
 }