2 #include <vtkCellArray.h>
3 #include <vtkPointData.h>
4 #include <vtkImageData.h>
5 #include <vtkMetaImageReader.h>
6 #include <vtkMetaImageWriter.h>
7 #include <vtkPNGReader.h>
8 #include <vtkPolyData.h>
9 #include <vtkCellLocator.h>
16 #include "clitkGammaIndex_ggo.h"
18 vtkImageData *loadImage(const std::string &filename) {
19 vtkImageReader2 *reader = vtkMetaImageReader::New();
20 //vtkImageReader2 *reader = vtkPNGReader::New();
21 reader->SetFileName(filename.c_str());
24 vtkImageData *image = reader->GetOutput();
25 image->Register(NULL);
31 void saveImage(vtkImageData *image,const std::string &filename) {
32 cout << "saving " << filename << endl;
33 vtkImageWriter *writer = vtkMetaImageWriter::New();
34 writer->SetFileName(filename.c_str());
35 writer->SetInput(image);
40 void insertTriangles(vtkCellArray *cells, vtkPoints *points, const vtkIdType ids[4]) {
41 double p0[3]; points->GetPoint(ids[0],p0);
42 double p1[3]; points->GetPoint(ids[1],p1);
43 double p2[3]; points->GetPoint(ids[2],p2);
44 double p3[3]; points->GetPoint(ids[3],p3);
45 //cout << "----------------------------------" << endl;
46 //cout << "p0=[" << p0[0] << "," << p0[1] << "," << p0[2] << "]" << endl;
47 //cout << "p1=[" << p1[0] << "," << p1[1] << "," << p1[2] << "]" << endl;
48 //cout << "p2=[" << p2[0] << "," << p2[1] << "," << p2[2] << "]" << endl;
49 //cout << "p3=[" << p3[0] << "," << p3[1] << "," << p3[2] << "]" << endl;
51 double center[] = {0,0,0};
52 for (int kk=0; kk<3; kk++) {
60 vtkIdType center_id = points->InsertNextPoint(center);
61 //cout << "center=[" << center[0] << "," << center[1] << "," << center[2] << "]" << endl;
63 cells->InsertNextCell(3);
64 cells->InsertCellPoint(ids[0]);
65 cells->InsertCellPoint(ids[1]);
66 cells->InsertCellPoint(center_id);
68 cells->InsertNextCell(3);
69 cells->InsertCellPoint(ids[1]);
70 cells->InsertCellPoint(ids[3]);
71 cells->InsertCellPoint(center_id);
73 cells->InsertNextCell(3);
74 cells->InsertCellPoint(ids[3]);
75 cells->InsertCellPoint(ids[2]);
76 cells->InsertCellPoint(center_id);
78 cells->InsertNextCell(3);
79 cells->InsertCellPoint(ids[2]);
80 cells->InsertCellPoint(ids[0]);
81 cells->InsertCellPoint(center_id);
84 double getMaximum(vtkImageData *image) {
88 for (int kk=0; kk<image->GetNumberOfPoints(); kk++) {
89 double value = image->GetPointData()->GetScalars()->GetTuple1(kk);
97 if (maximum<value) maximum = value;
103 vtkPolyData *buildPlane(vtkImageData *image,double spatial_margin,double dose_margin) {
104 vtkPoints *points = vtkPoints::New();
105 for (int kk=0; kk<image->GetNumberOfPoints(); kk++) {
106 double *point = image->GetPoint(kk);
107 double value = image->GetPointData()->GetScalars()->GetTuple1(kk);
112 point[0] /= spatial_margin;
113 point[1] /= spatial_margin;
114 point[2] /= dose_margin;
117 vtkIdType point_id = points->InsertNextPoint(point);
118 assert(kk==point_id);
120 points->InsertNextPoint(point);
124 vtkCellArray *cells = vtkCellArray::New();
125 for (int kk=0; kk<image->GetNumberOfCells(); kk++) {
126 vtkCell *cell = image->GetCell(kk);
127 assert(cell->GetNumberOfPoints()==4);
128 insertTriangles(cells,points,cell->GetPointIds()->GetPointer(0));
131 vtkPolyData *data = vtkPolyData::New();
132 data->SetPoints(points);
133 data->SetPolys(cells);
140 void assert2D(vtkImageData *image) {
142 int *extent = image->GetWholeExtent();
143 assert(extent[4]==0);
144 assert(extent[5]==0);
148 int main(int argc,char * argv[])
150 args_info_clitkGammaIndex args_info;
152 if (cmdline_parser_clitkGammaIndex(argc, argv, &args_info) != 0)
155 if (!args_info.absolute_dose_margin_given && !args_info.relative_dose_margin_given) {
156 std::cerr << "Specify either relative or absolute dose margin" << endl;
160 bool verbose = args_info.verbose_flag;
162 std::string reference_filename(args_info.reference_arg);
163 std::string target_filename(args_info.target_arg);
164 std::string gamma_filename(args_info.output_arg);
165 double space_margin = args_info.spatial_margin_arg;
166 double dose_rel_margin = args_info.relative_dose_margin_arg;
167 double dose_margin = args_info.absolute_dose_margin_arg;
168 bool use_dose_margin = args_info.absolute_dose_margin_given;
171 cout << "reference_filename=" << reference_filename << endl;
172 cout << "target_filename=" << target_filename << endl;
173 cout << "gamma_filename=" << gamma_filename << endl;
174 cout << "space_margin=" << space_margin << endl;
175 if (use_dose_margin) cout << "dose_margin=" << dose_margin << endl;
176 else cout << "dose_rel_margin=" << dose_rel_margin << endl;
180 vtkImageData *reference = loadImage(reference_filename);
183 // intensity normalisation
184 if (!use_dose_margin) {
185 dose_margin = getMaximum(reference)*dose_rel_margin;
186 if (verbose) cout << "dose_margin=" << dose_margin << endl;
190 vtkPolyData *data = buildPlane(reference,space_margin,dose_margin);
193 vtkAbstractCellLocator *locator = vtkCellLocator::New();
194 locator->SetDataSet(data);
196 locator->CacheCellBoundsOn();
197 locator->AutomaticOn();
198 locator->BuildLocator();
201 vtkImageData *target = loadImage(target_filename);
205 vtkImageData *output = vtkImageData::New();
206 output->SetExtent(target->GetWholeExtent());
207 output->SetOrigin(target->GetOrigin());
208 output->SetSpacing(target->GetSpacing());
209 output->SetScalarTypeToFloat();
210 output->AllocateScalars();
213 unsigned long total = 0;
214 unsigned long over_one = 0;
215 for (int kk=0; kk<target->GetNumberOfPoints(); kk++) {
216 double *point = target->GetPoint(kk);
217 double value = target->GetPointData()->GetScalars()->GetTuple1(kk);
222 point[0] /= space_margin;
223 point[1] /= space_margin;
224 point[2] /= dose_margin;
226 double closest_point[3] = {0,0,0};
227 vtkIdType cell_id = 0;
229 double squared_distance = 0;
231 locator->FindClosestPoint(point,closest_point,cell_id,foo,squared_distance);
233 double distance = sqrt(squared_distance);
234 output->GetPointData()->GetScalars()->SetTuple1(kk,distance);
236 if (value>1) over_one++;
242 cout << "total=" << total << endl;
243 cout << "over_one=" << over_one << endl;
244 cout << "ratio=" << static_cast<float>(over_one)/total << endl;
250 saveImage(output,gamma_filename);