]> Creatis software - clitk.git/blob - tools/clitkGammaIndex.cxx
Update Landmark coordinates with transformation matrix
[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
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 }