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