]> Creatis software - clitk.git/blob - tools/clitkGammaIndex.cxx
Rolled back to middle slice on first load
[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 #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