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