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