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
  19 #include "clitkGammaIndex_ggo.h"
  20 #include "clitkIO.h"
  21 #include "clitkDD.h"
  22
  23 #include <iostream>
  24 #include <cmath>
  25 #include <cassert>
  26
  27 #include <vtkPoints.h>
  28 #include <vtkCellArray.h>
  29 #include <vtkPointData.h>
  30 #include <vtkImageData.h>
  31 #include <vtkMetaImageReader.h>
  32 #include <vtkMetaImageWriter.h>
  33 #include <vtkPNGReader.h>
  34 #include <vtkPolyData.h>
  35 #include <vtkCellLocator.h>
  36
  37 #include <vvImage.h>
  38 #include <vvImageReader.h>
  39 #include <vvImageWriter.h>
  40
  41 #include <itkImage.h>
  42 #include <itkImageRegionIterator.h>
  43 typedef itk::Image<double> OutputImageType;
  44 typedef itk::ImageRegionIterator<OutputImageType> OutputImageIterator;
  45
  46 using std::endl;
  47 using std::cout;
  48
  49 vtkImageData* loadImage(const std::string& filename)
  50 {
  51     vvImageReader::Pointer reader = vvImageReader::New();
  52     reader->SetInputFilename(filename);
  53     reader->Update();
  54     vvImage::Pointer vvimage = reader->GetOutput();
  55     if (!vvimage) { cerr << "can't load " << filename << endl; return NULL; }
  56
  57     vtkImageData *image = vtkImageData::New();
  58     image->DeepCopy(vvimage->GetFirstVTKImageData());
  59     return image;
  60 }
  61
  62 void saveImage(OutputImageType* image,const std::string &filename) {
  63     vvImage::Pointer vvimage = vvImage::New();
  64     vvimage->AddItkImage(image);
  65
  66     vvImageWriter::Pointer writer = vvImageWriter::New();
  67     writer->SetOutputFileName(filename.c_str());
  68     writer->SetInput(vvimage);
  69     writer->Update();
  70 }
  71
  72
  73 void insertTriangles(vtkCellArray *cells, vtkPoints *points, const vtkIdType ids[4]) {
  74     double p0[3]; points->GetPoint(ids[0],p0);
  75     double p1[3]; points->GetPoint(ids[1],p1);
  76     double p2[3]; points->GetPoint(ids[2],p2);
  77     double p3[3]; points->GetPoint(ids[3],p3);
  78     //cout << "----------------------------------" << endl;
  79     //cout << "p0=[" << p0[0] << "," << p0[1] << "," << p0[2] << "]" << endl;
  80     //cout << "p1=[" << p1[0] << "," << p1[1] << "," << p1[2] << "]" << endl;
  81     //cout << "p2=[" << p2[0] << "," << p2[1] << "," << p2[2] << "]" << endl;
  82     //cout << "p3=[" << p3[0] << "," << p3[1] << "," << p3[2] << "]" << endl;
  83
  84     double center[] = {0,0,0};
  85     for (int kk=0; kk<3; kk++) {
  86         center[kk] += p0[kk];
  87         center[kk] += p1[kk];
  88         center[kk] += p2[kk];
  89         center[kk] += p3[kk];
  90         center[kk] /= 4.;
  91     }
  92
  93     vtkIdType center_id = points->InsertNextPoint(center);
  94     //cout << "center=[" << center[0] << "," << center[1] << "," << center[2] << "]" << endl;
  95
  96     cells->InsertNextCell(3);
  97     cells->InsertCellPoint(ids[0]);
  98     cells->InsertCellPoint(ids[1]);
  99     cells->InsertCellPoint(center_id);
 100
 101     cells->InsertNextCell(3);
 102     cells->InsertCellPoint(ids[1]);
 103     cells->InsertCellPoint(ids[3]);
 104     cells->InsertCellPoint(center_id);
 105
 106     cells->InsertNextCell(3);
 107     cells->InsertCellPoint(ids[3]);
 108     cells->InsertCellPoint(ids[2]);
 109     cells->InsertCellPoint(center_id);
 110
 111     cells->InsertNextCell(3);
 112     cells->InsertCellPoint(ids[2]);
 113     cells->InsertCellPoint(ids[0]);
 114     cells->InsertCellPoint(center_id);
 115 }
 116
 117 void insertLine(vtkCellArray *cells, vtkPoints *points, const vtkIdType ids[2]) {
 118     cells->InsertNextCell(2);
 119     cells->InsertCellPoint(ids[0]);
 120     cells->InsertCellPoint(ids[1]);
 121 }
 122
 123 double getMaximum(vtkImageData *image) {
 124     bool first = true;
 125     double maximum = 0;
 126
 127     vtkPointData* point_data = image->GetPointData();
 128     assert(point_data);
 129     vtkDataArray* scalars = point_data->GetScalars();
 130     assert(scalars);
 131
 132     for (int kk=0; kk<image->GetNumberOfPoints(); kk++) {
 133         double value = scalars->GetTuple1(kk);
 134
 135         if (first) {
 136             maximum = value;
 137             first = false;
 138             continue;
 139         }
 140
 141         if (maximum<value) maximum = value;
 142     }
 143
 144     return maximum;
 145 }
 146
 147
 148 vtkPolyData *buildPlane(vtkImageData *image,double spatial_margin,double dose_margin) {
 149     vtkPoints *points = vtkPoints::New();
 150     for (int kk=0; kk<image->GetNumberOfPoints(); kk++) {
 151         double *point = image->GetPoint(kk);
 152         double value = image->GetPointData()->GetScalars()->GetTuple1(kk);
 153         assert(value>=0);
 154         assert(point[2]==0);
 155         point[2] = value;
 156
 157         point[0] /= spatial_margin;
 158         point[1] /= spatial_margin;
 159         point[2] /= dose_margin;
 160
 161 #ifndef NDEBUG
 162         vtkIdType point_id = points->InsertNextPoint(point);
 163         assert(kk==point_id);
 164 #else
 165         points->InsertNextPoint(point);
 166 #endif
 167     }
 168
 169     vtkCellArray *cells = vtkCellArray::New();
 170     for (int kk=0; kk<image->GetNumberOfCells(); kk++) {
 171         vtkCell *cell = image->GetCell(kk);
 172
 173         if (cell->GetNumberOfPoints()==4) {
 174             insertTriangles(cells,points,cell->GetPointIds()->GetPointer(0));
 175             continue;
 176         }
 177
 178         if (cell->GetNumberOfPoints()==2) {
 179             insertLine(cells,points,cell->GetPointIds()->GetPointer(0));
 180             continue;
 181         }
 182
 183         assert(false);
 184     }
 185
 186     vtkPolyData *data = vtkPolyData::New();
 187     data->SetPoints(points);
 188     data->SetPolys(cells);
 189     points->Delete();
 190     cells->Delete();
 191
 192     return data;
 193 }
 194
 195 int main(int argc,char * argv[])
 196 {
 197     clitk::RegisterClitkFactories();
 198
 199     args_info_clitkGammaIndex args_info;
 200
 201     if (cmdline_parser_clitkGammaIndex(argc, argv, &args_info) != 0)
 202         exit(1);
 203
 204     if (!args_info.absolute_dose_margin_given && !args_info.relative_dose_margin_given) {
 205         std::cerr << "Specify either relative or absolute dose margin" << endl;
 206         exit(1);
 207     }
 208
 209     bool verbose = args_info.verbose_flag;
 210
 211     std::string reference_filename(args_info.reference_arg);
 212     std::string target_filename(args_info.target_arg);
 213     std::string gamma_filename(args_info.output_arg);
 214     double space_margin = args_info.spatial_margin_arg;
 215     double dose_rel_margin = args_info.relative_dose_margin_arg;
 216     double dose_margin = args_info.absolute_dose_margin_arg;
 217     bool use_dose_margin = args_info.absolute_dose_margin_given;
 218
 219     if (verbose) {
 220         cout << "reference_filename=" << reference_filename << endl;
 221         cout << "target_filename=" << target_filename << endl;
 222         cout << "gamma_filename=" << gamma_filename << endl;
 223         cout << "space_margin=" << space_margin << endl;
 224         if (use_dose_margin) cout << "dose_margin=" << dose_margin << endl;
 225         else cout << "dose_rel_margin=" << dose_rel_margin << endl;
 226     }
 227
 228     // load reference
 229     vtkImageData* reference = loadImage(reference_filename);
 230     assert(reference);
 231
 232     // translate target with arguments values
 233     // reference is translated instead of target so that the output space stay the same as target
 234     {
 235         double reference_origin[3];
 236         reference->GetOrigin(reference_origin);
 237         reference_origin[0] -= args_info.translation_x_arg;
 238         reference_origin[1] -= args_info.translation_y_arg;
 239         reference_origin[2] -= args_info.translation_z_arg;
 240         reference->SetOrigin(reference_origin);
 241     }
 242
 243     // intensity normalisation
 244     if (!use_dose_margin) {
 245         dose_margin = getMaximum(reference)*dose_rel_margin;
 246         if (verbose) cout << "dose_margin=" << dose_margin << endl;
 247     }
 248
 249     // build surface
 250     vtkPolyData *data = buildPlane(reference,space_margin,dose_margin);
 251     reference->Delete();
 252
 253     vtkAbstractCellLocator *locator = vtkCellLocator::New();
 254     locator->SetDataSet(data);
 255     DD("here");
 256     //    data->Delete();
 257     locator->CacheCellBoundsOn();
 258     locator->AutomaticOn();
 259     DD("BuildLocator");
 260     locator->BuildLocator();
 261     DD("end BuildLocator");
 262
 263     // load target
 264     vtkImageData* target = loadImage(target_filename);
 265     assert(target);
 266
 267
 268     // allocate output
 269     OutputImageType::Pointer output = OutputImageType::New();
 270     {
 271         OutputImageType::SizeType::SizeValueType output_array_size[2];
 272         output_array_size[0] = target->GetDimensions()[0];
 273         output_array_size[1] = target->GetDimensions()[1];
 274         OutputImageType::SizeType output_size;
 275         output_size.SetSize(output_array_size);
 276         output->SetRegions(OutputImageType::RegionType(output_size));
 277         output->SetOrigin(target->GetOrigin());
 278         output->SetSpacing(target->GetSpacing());
 279         output->Allocate();
 280     }
 281
 282     // fill output
 283     unsigned long  kk = 0;
 284     unsigned long over_one = 0;
 285     OutputImageIterator iter(output,output->GetLargestPossibleRegion());
 286     iter.GoToBegin();
 287     DD("loop");
 288     while (!iter.IsAtEnd()) {
 289         double *point = target->GetPoint(kk);
 290         double value = target->GetPointData()->GetScalars()->GetTuple1(kk);
 291         assert(value>=0);
 292         assert(point[2]==0);
 293         point[2] = value;
 294
 295         point[0] /= space_margin;
 296         point[1] /= space_margin;
 297         point[2] /= dose_margin;
 298
 299         double closest_point[3] = {0,0,0};
 300         vtkIdType cell_id = 0;
 301         int foo = 0;
 302         double squared_distance = 0;
 303
 304         locator->FindClosestPoint(point,closest_point,cell_id,foo,squared_distance);
 305         double distance = sqrt(squared_distance);
 306         iter.Set(distance);
 307
 308         if (distance>1) over_one++;
 309         kk++;
 310         ++iter;
 311     }
 312
 313     if (verbose) {
 314         cout << "total=" << kk << endl;
 315         cout << "over_one=" << over_one << endl;
 316         cout << "ratio=" << static_cast<float>(over_one)/kk << endl;
 317     }
 318
 319     locator->Delete();
 320     target->Delete();
 321
 322     saveImage(output,gamma_filename);
 323
 324     return 0;
 325 }