- // scale intensity
- Normalizer::Pointer reference_normalizer = Normalizer::New();
- Normalizer::Pointer target_normalizer = Normalizer::New();
- {
- reference_normalizer->SetShift(0);
- reference_normalizer->SetScale(1/dose_margin);
- reference_normalizer->SetInput(reference_reader->GetOutput());
- reference_normalizer->Update();
-
- target_normalizer->SetShift(0);
- target_normalizer->SetScale(1/dose_margin);
- target_normalizer->SetInput(target_reader->GetOutput());
- target_normalizer->Update();
-
- //cout << "scale=" << reference_normalizer->GetScale() << "/" << target_normalizer->GetScale() << endl;
- //cout << "shift=" << reference_normalizer->GetShift() << "/" << target_normalizer->GetShift() << endl;
- }
+ // build surface
+ vtkPolyData *data = buildPlane(reference,space_margin,dose_margin);
+ reference->Delete();
+
+ vtkAbstractCellLocator *locator = vtkCellLocator::New();
+ locator->SetDataSet(data);
+ data->Delete();
+ locator->CacheCellBoundsOn();
+ locator->AutomaticOn();
+ locator->BuildLocator();
+
+ // load target
+ vtkImageData *target = loadImage(target_filename);
+ assert2D(target);
+
+ // allocate output
+ vtkImageData *output = vtkImageData::New();
+ output->SetExtent(target->GetWholeExtent());
+ output->SetOrigin(target->GetOrigin());
+ output->SetSpacing(target->GetSpacing());
+ output->SetScalarTypeToFloat();
+ output->AllocateScalars();
+
+ // fill output
+ unsigned long total = 0;
+ unsigned long over_one = 0;
+ for (int kk=0; kk<target->GetNumberOfPoints(); kk++) {
+ double *point = target->GetPoint(kk);
+ double value = target->GetPointData()->GetScalars()->GetTuple1(kk);
+ assert(value>=0);
+ assert(point[2]==0);
+ point[2] = value;
+
+ point[0] /= space_margin;
+ point[1] /= space_margin;
+ point[2] /= dose_margin;
+
+ double closest_point[3] = {0,0,0};
+ vtkIdType cell_id = 0;
+ int foo = 0;
+ double squared_distance = 0;
+
+ locator->FindClosestPoint(point,closest_point,cell_id,foo,squared_distance);