1 /*=========================================================================
2 Program: vv http://www.creatis.insa-lyon.fr/rio/vv
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
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.
13 It is distributed under dual licence
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>
33 #include "clitkGammaIndex_ggo.h"
37 #include <vvImageReader.h>
39 vtkImageData* loadImage(const std::string& filename)
41 vvImageReader::Pointer reader = vvImageReader::New();
42 reader->SetInputFilename(filename);
44 vvImage::Pointer vvimage = reader->GetOutput();
45 if (!vvimage) { cerr << "can't load " << filename << endl; return NULL; }
47 vtkImageData *image = vtkImageData::New();
48 image->DeepCopy(vvimage->GetFirstVTKImageData());
52 void saveImage(vtkImageData *image,const std::string &filename) {
53 cout << "saving " << filename << endl;
54 vtkImageWriter *writer = vtkMetaImageWriter::New();
55 writer->SetFileName(filename.c_str());
56 writer->SetInput(image);
62 void insertTriangles(vtkCellArray *cells, vtkPoints *points, const vtkIdType ids[4]) {
63 double p0[3]; points->GetPoint(ids[0],p0);
64 double p1[3]; points->GetPoint(ids[1],p1);
65 double p2[3]; points->GetPoint(ids[2],p2);
66 double p3[3]; points->GetPoint(ids[3],p3);
67 //cout << "----------------------------------" << endl;
68 //cout << "p0=[" << p0[0] << "," << p0[1] << "," << p0[2] << "]" << endl;
69 //cout << "p1=[" << p1[0] << "," << p1[1] << "," << p1[2] << "]" << endl;
70 //cout << "p2=[" << p2[0] << "," << p2[1] << "," << p2[2] << "]" << endl;
71 //cout << "p3=[" << p3[0] << "," << p3[1] << "," << p3[2] << "]" << endl;
73 double center[] = {0,0,0};
74 for (int kk=0; kk<3; kk++) {
82 vtkIdType center_id = points->InsertNextPoint(center);
83 //cout << "center=[" << center[0] << "," << center[1] << "," << center[2] << "]" << endl;
85 cells->InsertNextCell(3);
86 cells->InsertCellPoint(ids[0]);
87 cells->InsertCellPoint(ids[1]);
88 cells->InsertCellPoint(center_id);
90 cells->InsertNextCell(3);
91 cells->InsertCellPoint(ids[1]);
92 cells->InsertCellPoint(ids[3]);
93 cells->InsertCellPoint(center_id);
95 cells->InsertNextCell(3);
96 cells->InsertCellPoint(ids[3]);
97 cells->InsertCellPoint(ids[2]);
98 cells->InsertCellPoint(center_id);
100 cells->InsertNextCell(3);
101 cells->InsertCellPoint(ids[2]);
102 cells->InsertCellPoint(ids[0]);
103 cells->InsertCellPoint(center_id);
106 void insertLine(vtkCellArray *cells, vtkPoints *points, const vtkIdType ids[2]) {
107 cells->InsertNextCell(2);
108 cells->InsertCellPoint(ids[0]);
109 cells->InsertCellPoint(ids[1]);
112 double getMaximum(vtkImageData *image) {
116 vtkPointData* point_data = image->GetPointData();
118 vtkDataArray* scalars = point_data->GetScalars();
121 for (int kk=0; kk<image->GetNumberOfPoints(); kk++) {
122 double value = scalars->GetTuple1(kk);
130 if (maximum<value) maximum = value;
137 vtkPolyData *buildPlane(vtkImageData *image,double spatial_margin,double dose_margin) {
138 vtkPoints *points = vtkPoints::New();
139 for (int kk=0; kk<image->GetNumberOfPoints(); kk++) {
140 double *point = image->GetPoint(kk);
141 double value = image->GetPointData()->GetScalars()->GetTuple1(kk);
146 point[0] /= spatial_margin;
147 point[1] /= spatial_margin;
148 point[2] /= dose_margin;
151 vtkIdType point_id = points->InsertNextPoint(point);
152 assert(kk==point_id);
154 points->InsertNextPoint(point);
158 vtkCellArray *cells = vtkCellArray::New();
159 for (int kk=0; kk<image->GetNumberOfCells(); kk++) {
160 vtkCell *cell = image->GetCell(kk);
162 if (cell->GetNumberOfPoints()==4) {
163 insertTriangles(cells,points,cell->GetPointIds()->GetPointer(0));
167 if (cell->GetNumberOfPoints()==2) {
168 insertLine(cells,points,cell->GetPointIds()->GetPointer(0));
175 vtkPolyData *data = vtkPolyData::New();
176 data->SetPoints(points);
177 data->SetPolys(cells);
184 int main(int argc,char * argv[])
186 clitk::RegisterClitkFactories();
188 args_info_clitkGammaIndex args_info;
190 if (cmdline_parser_clitkGammaIndex(argc, argv, &args_info) != 0)
193 if (!args_info.absolute_dose_margin_given && !args_info.relative_dose_margin_given) {
194 std::cerr << "Specify either relative or absolute dose margin" << endl;
198 bool verbose = args_info.verbose_flag;
200 std::string reference_filename(args_info.reference_arg);
201 std::string target_filename(args_info.target_arg);
202 std::string gamma_filename(args_info.output_arg);
203 double space_margin = args_info.spatial_margin_arg;
204 double dose_rel_margin = args_info.relative_dose_margin_arg;
205 double dose_margin = args_info.absolute_dose_margin_arg;
206 bool use_dose_margin = args_info.absolute_dose_margin_given;
209 cout << "reference_filename=" << reference_filename << endl;
210 cout << "target_filename=" << target_filename << endl;
211 cout << "gamma_filename=" << gamma_filename << endl;
212 cout << "space_margin=" << space_margin << endl;
213 if (use_dose_margin) cout << "dose_margin=" << dose_margin << endl;
214 else cout << "dose_rel_margin=" << dose_rel_margin << endl;
218 vtkImageData* reference = loadImage(reference_filename);
221 // intensity normalisation
222 if (!use_dose_margin) {
223 dose_margin = getMaximum(reference)*dose_rel_margin;
224 if (verbose) cout << "dose_margin=" << dose_margin << endl;
228 vtkPolyData *data = buildPlane(reference,space_margin,dose_margin);
231 vtkAbstractCellLocator *locator = vtkCellLocator::New();
232 locator->SetDataSet(data);
234 locator->CacheCellBoundsOn();
235 locator->AutomaticOn();
236 locator->BuildLocator();
239 vtkImageData *target = loadImage(target_filename);
243 vtkImageData *output = vtkImageData::New();
244 output->SetExtent(target->GetWholeExtent());
245 output->SetOrigin(target->GetOrigin());
246 output->SetSpacing(target->GetSpacing());
247 output->SetScalarTypeToFloat();
248 output->AllocateScalars();
251 unsigned long total = 0;
252 unsigned long over_one = 0;
253 for (int kk=0; kk<target->GetNumberOfPoints(); kk++) {
254 double *point = target->GetPoint(kk);
255 double value = target->GetPointData()->GetScalars()->GetTuple1(kk);
260 point[0] /= space_margin;
261 point[1] /= space_margin;
262 point[2] /= dose_margin;
264 double closest_point[3] = {0,0,0};
265 vtkIdType cell_id = 0;
267 double squared_distance = 0;
269 locator->FindClosestPoint(point,closest_point,cell_id,foo,squared_distance);
271 double distance = sqrt(squared_distance);
272 output->GetPointData()->GetScalars()->SetTuple1(kk,distance);
274 if (value>1) over_one++;
280 cout << "total=" << total << endl;
281 cout << "over_one=" << over_one << endl;
282 cout << "ratio=" << static_cast<float>(over_one)/total << endl;
288 saveImage(output,gamma_filename);