tools/clitkGammaIndex.cxx

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