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