tools/clitkGammaIndex.cxx

   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 #include <vtkPoints.h>
  19 #include <vtkCellArray.h>
  20 #include <vtkPointData.h>
  21 #include <vtkImageData.h>
  22 #include <vtkMetaImageReader.h>
  23 #include <vtkMetaImageWriter.h>
  24 #include <vtkPNGReader.h>
  25 #include <vtkPolyData.h>
  26 #include <vtkCellLocator.h>
  27 #include <iostream>
  28 #include <cmath>
  29 #include <cassert>
  30 using std::endl;
  31 using std::cout;
  32
  33 #include "clitkGammaIndex_ggo.h"
  34
  35 vtkImageData *loadImage(const std::string &filename) {
  36   vtkImageReader2 *reader = vtkMetaImageReader::New();
  37   //vtkImageReader2 *reader = vtkPNGReader::New();
  38   reader->SetFileName(filename.c_str());
  39   reader->Update();
  40
  41   vtkImageData *image = reader->GetOutput();
  42   image->Register(NULL);
  43   reader->Delete();
  44
  45   return image;
  46 }
  47
  48 void saveImage(vtkImageData *image,const std::string &filename) {
  49   cout << "saving " << filename << endl;
  50   vtkImageWriter *writer = vtkMetaImageWriter::New();
  51   writer->SetFileName(filename.c_str());
  52   writer->SetInput(image);
  53   writer->Write();
  54   writer->Delete();
  55 }
  56
  57 void insertTriangles(vtkCellArray *cells, vtkPoints *points, const vtkIdType ids[4]) {
  58   double p0[3]; points->GetPoint(ids[0],p0);
  59   double p1[3]; points->GetPoint(ids[1],p1);
  60   double p2[3]; points->GetPoint(ids[2],p2);
  61   double p3[3]; points->GetPoint(ids[3],p3);
  62   //cout << "----------------------------------" << endl;
  63   //cout << "p0=[" << p0[0] << "," << p0[1] << "," << p0[2] << "]" << endl;
  64   //cout << "p1=[" << p1[0] << "," << p1[1] << "," << p1[2] << "]" << endl;
  65   //cout << "p2=[" << p2[0] << "," << p2[1] << "," << p2[2] << "]" << endl;
  66   //cout << "p3=[" << p3[0] << "," << p3[1] << "," << p3[2] << "]" << endl;
  67
  68   double center[] = {0,0,0};
  69   for (int kk=0; kk<3; kk++) {
  70     center[kk] += p0[kk];
  71     center[kk] += p1[kk];
  72     center[kk] += p2[kk];
  73     center[kk] += p3[kk];
  74     center[kk] /= 4.;
  75   }
  76
  77   vtkIdType center_id = points->InsertNextPoint(center);
  78   //cout << "center=[" << center[0] << "," << center[1] << "," << center[2] << "]" << endl;
  79
  80   cells->InsertNextCell(3);
  81   cells->InsertCellPoint(ids[0]);
  82   cells->InsertCellPoint(ids[1]);
  83   cells->InsertCellPoint(center_id);
  84
  85   cells->InsertNextCell(3);
  86   cells->InsertCellPoint(ids[1]);
  87   cells->InsertCellPoint(ids[3]);
  88   cells->InsertCellPoint(center_id);
  89
  90   cells->InsertNextCell(3);
  91   cells->InsertCellPoint(ids[3]);
  92   cells->InsertCellPoint(ids[2]);
  93   cells->InsertCellPoint(center_id);
  94
  95   cells->InsertNextCell(3);
  96   cells->InsertCellPoint(ids[2]);
  97   cells->InsertCellPoint(ids[0]);
  98   cells->InsertCellPoint(center_id);
  99 }
 100
 101 double getMaximum(vtkImageData *image) {
 102   bool first = true;
 103   double maximum = 0;
 104
 105   for (int kk=0; kk<image->GetNumberOfPoints(); kk++) {
 106     double value = image->GetPointData()->GetScalars()->GetTuple1(kk);
 107
 108     if (first) {
 109       maximum = value;
 110       first = false;
 111       continue;
 112     }
 113
 114     if (maximum<value) maximum = value;
 115   }
 116
 117   return maximum;
 118 }
 119
 120 vtkPolyData *buildPlane(vtkImageData *image,double spatial_margin,double dose_margin) {
 121   vtkPoints *points = vtkPoints::New();
 122   for (int kk=0; kk<image->GetNumberOfPoints(); kk++) {
 123     double *point = image->GetPoint(kk);
 124     double value = image->GetPointData()->GetScalars()->GetTuple1(kk);
 125     assert(value>=0);
 126     assert(point[2]==0);
 127     point[2] = value;
 128
 129     point[0] /= spatial_margin;
 130     point[1] /= spatial_margin;
 131     point[2] /= dose_margin;
 132
 133 #ifndef NDEBUG
 134     vtkIdType point_id = points->InsertNextPoint(point);
 135     assert(kk==point_id);
 136 #else
 137     points->InsertNextPoint(point);
 138 #endif
 139   }
 140
 141   vtkCellArray *cells = vtkCellArray::New();
 142   for (int kk=0; kk<image->GetNumberOfCells(); kk++) {
 143     vtkCell *cell = image->GetCell(kk);
 144     assert(cell->GetNumberOfPoints()==4);
 145     insertTriangles(cells,points,cell->GetPointIds()->GetPointer(0));
 146   }
 147
 148   vtkPolyData *data = vtkPolyData::New();
 149   data->SetPoints(points);
 150   data->SetPolys(cells);
 151   points->Delete();
 152   cells->Delete();
 153
 154   return data;
 155 }
 156
 157 void assert2D(vtkImageData *image) {
 158 #ifndef NDEBUG
 159   int *extent = image->GetWholeExtent();
 160   assert(extent[4]==0);
 161   assert(extent[5]==0);
 162 #endif
 163 }
 164
 165 int main(int argc,char * argv[])
 166 {
 167   args_info_clitkGammaIndex args_info;
 168
 169   if (cmdline_parser_clitkGammaIndex(argc, argv, &args_info) != 0)
 170     exit(1);
 171
 172   if (!args_info.absolute_dose_margin_given && !args_info.relative_dose_margin_given) {
 173     std::cerr << "Specify either relative or absolute dose margin" << endl;
 174     exit(1);
 175   }
 176
 177   bool verbose = args_info.verbose_flag;
 178
 179   std::string reference_filename(args_info.reference_arg);
 180   std::string target_filename(args_info.target_arg);
 181   std::string gamma_filename(args_info.output_arg);
 182   double space_margin = args_info.spatial_margin_arg;
 183   double dose_rel_margin = args_info.relative_dose_margin_arg;
 184   double dose_margin = args_info.absolute_dose_margin_arg;
 185   bool use_dose_margin = args_info.absolute_dose_margin_given;
 186
 187   if (verbose) {
 188     cout << "reference_filename=" << reference_filename << endl;
 189     cout << "target_filename=" << target_filename << endl;
 190     cout << "gamma_filename=" << gamma_filename << endl;
 191     cout << "space_margin=" << space_margin << endl;
 192     if (use_dose_margin) cout << "dose_margin=" << dose_margin << endl;
 193     else cout << "dose_rel_margin=" << dose_rel_margin << endl;
 194   }
 195
 196   // load reference
 197   vtkImageData *reference = loadImage(reference_filename);
 198   assert2D(reference);
 199
 200   // intensity normalisation
 201   if (!use_dose_margin) {
 202     dose_margin = getMaximum(reference)*dose_rel_margin;
 203     if (verbose) cout << "dose_margin=" << dose_margin << endl;
 204   }
 205
 206   // build surface
 207   vtkPolyData *data = buildPlane(reference,space_margin,dose_margin);
 208   reference->Delete();
 209
 210   vtkAbstractCellLocator *locator = vtkCellLocator::New();
 211   locator->SetDataSet(data);
 212   data->Delete();
 213   locator->CacheCellBoundsOn();
 214   locator->AutomaticOn();
 215   locator->BuildLocator();
 216
 217   // load target
 218   vtkImageData *target = loadImage(target_filename);
 219   assert2D(target);
 220
 221   // allocate output
 222   vtkImageData *output = vtkImageData::New();
 223   output->SetExtent(target->GetWholeExtent());
 224   output->SetOrigin(target->GetOrigin());
 225   output->SetSpacing(target->GetSpacing());
 226   output->SetScalarTypeToFloat();
 227   output->AllocateScalars();
 228
 229   // fill output
 230   unsigned long total = 0;
 231   unsigned long over_one = 0;
 232   for (int kk=0; kk<target->GetNumberOfPoints(); kk++) {
 233     double *point = target->GetPoint(kk);
 234     double value = target->GetPointData()->GetScalars()->GetTuple1(kk);
 235     assert(value>=0);
 236     assert(point[2]==0);
 237     point[2] = value;
 238
 239     point[0] /= space_margin;
 240     point[1] /= space_margin;
 241     point[2] /= dose_margin;
 242
 243     double closest_point[3] = {0,0,0};
 244     vtkIdType cell_id = 0;
 245     int foo = 0;
 246     double squared_distance = 0;
 247
 248     locator->FindClosestPoint(point,closest_point,cell_id,foo,squared_distance);
 249
 250     double distance = sqrt(squared_distance);
 251     output->GetPointData()->GetScalars()->SetTuple1(kk,distance);
 252
 253     if (value>1) over_one++;
 254     total++;
 255
 256   }
 257
 258   if (verbose) {
 259     cout << "total=" << total << endl;
 260     cout << "over_one=" << over_one << endl;
 261     cout << "ratio=" << static_cast<float>(over_one)/total << endl;
 262   }
 263
 264   locator->Delete();
 265   target->Delete();
 266
 267   saveImage(output,gamma_filename);
 268   output->Delete();
 269
 270   return 0;
 271 }
 272