- BSD See included LICENSE.txt file
- CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
===========================================================================*/
+
+#include "clitkGammaIndex_ggo.h"
+#include "clitkIO.h"
+#include "clitkDD.h"
+
+#include <iostream>
+#include <cmath>
+#include <cassert>
+
#include <vtkPoints.h>
#include <vtkCellArray.h>
#include <vtkPointData.h>
#include <vtkPNGReader.h>
#include <vtkPolyData.h>
#include <vtkCellLocator.h>
-#include <iostream>
-#include <cmath>
-#include <cassert>
+
+#include <vvImage.h>
+#include <vvImageReader.h>
+#include <vvImageWriter.h>
+
+#include <itkImage.h>
+#include <itkImageRegionIterator.h>
+typedef itk::Image<double> OutputImageType;
+typedef itk::ImageRegionIterator<OutputImageType> OutputImageIterator;
+
using std::endl;
using std::cout;
-#include "clitkGammaIndex_ggo.h"
+vtkImageData* loadImage(const std::string& filename)
+{
+ vvImageReader::Pointer reader = vvImageReader::New();
+ reader->SetInputFilename(filename);
+ reader->Update();
+ vvImage::Pointer vvimage = reader->GetOutput();
+ if (!vvimage) { cerr << "can't load " << filename << endl; return NULL; }
+
+ vtkImageData *image = vtkImageData::New();
+ image->DeepCopy(vvimage->GetFirstVTKImageData());
+ return image;
+}
-#include <vvImage.h>
-#include <vvImageReader.h>
+void saveImage(OutputImageType* image,const std::string &filename) {
+ vvImage::Pointer vvimage = vvImage::New();
+ vvimage->AddItkImage(image);
-void saveImage(vtkImageData *image,const std::string &filename) {
- cout << "saving " << filename << endl;
- vtkImageWriter *writer = vtkMetaImageWriter::New();
- writer->SetFileName(filename.c_str());
- writer->SetInput(image);
- writer->Write();
- writer->Delete();
+ vvImageWriter::Pointer writer = vvImageWriter::New();
+ writer->SetOutputFileName(filename.c_str());
+ writer->SetInput(vvimage);
+ writer->Update();
}
cells->InsertCellPoint(center_id);
}
+void insertLine(vtkCellArray *cells, vtkPoints *points, const vtkIdType ids[2]) {
+ cells->InsertNextCell(2);
+ cells->InsertCellPoint(ids[0]);
+ cells->InsertCellPoint(ids[1]);
+}
+
double getMaximum(vtkImageData *image) {
bool first = true;
double maximum = 0;
+ vtkPointData* point_data = image->GetPointData();
+ assert(point_data);
+ vtkDataArray* scalars = point_data->GetScalars();
+ assert(scalars);
+
for (int kk=0; kk<image->GetNumberOfPoints(); kk++) {
- double value = image->GetPointData()->GetScalars()->GetTuple1(kk);
+ double value = scalars->GetTuple1(kk);
if (first) {
maximum = value;
vtkCellArray *cells = vtkCellArray::New();
for (int kk=0; kk<image->GetNumberOfCells(); kk++) {
vtkCell *cell = image->GetCell(kk);
- assert(cell->GetNumberOfPoints()==4);
- insertTriangles(cells,points,cell->GetPointIds()->GetPointer(0));
+
+ if (cell->GetNumberOfPoints()==4) {
+ insertTriangles(cells,points,cell->GetPointIds()->GetPointer(0));
+ continue;
+ }
+
+ if (cell->GetNumberOfPoints()==2) {
+ insertLine(cells,points,cell->GetPointIds()->GetPointer(0));
+ continue;
+ }
+
+ assert(false);
}
vtkPolyData *data = vtkPolyData::New();
int main(int argc,char * argv[])
{
+ clitk::RegisterClitkFactories();
+
args_info_clitkGammaIndex args_info;
if (cmdline_parser_clitkGammaIndex(argc, argv, &args_info) != 0)
}
// load reference
- vtkImageData* reference = NULL;
+ vtkImageData* reference = loadImage(reference_filename);
+ assert(reference);
+
+ // translate target with arguments values
+ // reference is translated instead of target so that the output space stay the same as target
{
- vvImageReader::Pointer reader = vvImageReader::New();
- reader->SetInputFilename(reference_filename);
- reader->Update();
- vvImage::Pointer vvimage = reader->GetOutput();
- reference = vvimage->GetFirstVTKImageData();
+ double reference_origin[3];
+ reference->GetOrigin(reference_origin);
+ reference_origin[0] -= args_info.translation_x_arg;
+ reference_origin[1] -= args_info.translation_y_arg;
+ reference_origin[2] -= args_info.translation_z_arg;
+ reference->SetOrigin(reference_origin);
}
// intensity normalisation
vtkAbstractCellLocator *locator = vtkCellLocator::New();
locator->SetDataSet(data);
- data->Delete();
+ DD("here");
+ // data->Delete();
locator->CacheCellBoundsOn();
locator->AutomaticOn();
+ DD("BuildLocator");
locator->BuildLocator();
+ DD("end BuildLocator");
// load target
- vtkImageData *target = NULL;
- {
- vvImageReader::Pointer reader = vvImageReader::New();
- reader->SetInputFilename(target_filename);
- reader->Update();
- vvImage::Pointer vvimage = reader->GetOutput();
- target = vvimage->GetFirstVTKImageData();
- }
+ vtkImageData* target = loadImage(target_filename);
+ assert(target);
+
// allocate output
- vtkImageData *output = vtkImageData::New();
- output->SetExtent(target->GetWholeExtent());
- output->SetOrigin(target->GetOrigin());
- output->SetSpacing(target->GetSpacing());
- output->SetScalarTypeToFloat();
- output->AllocateScalars();
+ OutputImageType::Pointer output = OutputImageType::New();
+ {
+ OutputImageType::SizeType::SizeValueType output_array_size[2];
+ output_array_size[0] = target->GetDimensions()[0];
+ output_array_size[1] = target->GetDimensions()[1];
+ OutputImageType::SizeType output_size;
+ output_size.SetSize(output_array_size);
+ output->SetRegions(OutputImageType::RegionType(output_size));
+ output->SetOrigin(target->GetOrigin());
+ output->SetSpacing(target->GetSpacing());
+ output->Allocate();
+ }
// fill output
- unsigned long total = 0;
+ unsigned long kk = 0;
unsigned long over_one = 0;
- for (int kk=0; kk<target->GetNumberOfPoints(); kk++) {
+ OutputImageIterator iter(output,output->GetLargestPossibleRegion());
+ iter.GoToBegin();
+ DD("loop");
+ while (!iter.IsAtEnd()) {
double *point = target->GetPoint(kk);
double value = target->GetPointData()->GetScalars()->GetTuple1(kk);
assert(value>=0);
double squared_distance = 0;
locator->FindClosestPoint(point,closest_point,cell_id,foo,squared_distance);
-
double distance = sqrt(squared_distance);
- output->GetPointData()->GetScalars()->SetTuple1(kk,distance);
-
- if (value>1) over_one++;
- total++;
+ iter.Set(distance);
+ if (distance>1) over_one++;
+ kk++;
+ ++iter;
}
if (verbose) {
- cout << "total=" << total << endl;
+ cout << "total=" << kk << endl;
cout << "over_one=" << over_one << endl;
- cout << "ratio=" << static_cast<float>(over_one)/total << endl;
+ cout << "ratio=" << static_cast<float>(over_one)/kk << endl;
}
locator->Delete();
target->Delete();
saveImage(output,gamma_filename);
- output->Delete();
return 0;
}
-
git://git.creatis.insa-lyon.fr / clitk.git / blobdiff