X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=tools%2FclitkGammaIndex.cxx;h=1f5fd03f3a78c905ecd666dc2376569b674e9fa3;hb=1512c37902327dde2650268e1fbd7fbc085b7a89;hp=e2369d0dbc0c2b9663480eae9714f7b8c1f4d9b8;hpb=c3358d7e514585e32f26e91d1c450d1c129929bc;p=clitk.git diff --git a/tools/clitkGammaIndex.cxx b/tools/clitkGammaIndex.cxx index e2369d0..1f5fd03 100644 --- a/tools/clitkGammaIndex.cxx +++ b/tools/clitkGammaIndex.cxx @@ -15,6 +15,15 @@ - 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 +#include +#include + #include #include #include @@ -24,24 +33,40 @@ #include #include #include -#include -#include -#include + +#include +#include +#include + +#include +#include +typedef itk::Image OutputImageType; +typedef itk::ImageRegionIterator 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 -#include +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(); } @@ -89,12 +114,23 @@ void insertTriangles(vtkCellArray *cells, vtkPoints *points, const vtkIdType ids 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; kkGetNumberOfPoints(); kk++) { - double value = image->GetPointData()->GetScalars()->GetTuple1(kk); + double value = scalars->GetTuple1(kk); if (first) { maximum = value; @@ -133,8 +169,18 @@ vtkPolyData *buildPlane(vtkImageData *image,double spatial_margin,double dose_ma vtkCellArray *cells = vtkCellArray::New(); for (int kk=0; kkGetNumberOfCells(); 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(); @@ -148,6 +194,8 @@ vtkPolyData *buildPlane(vtkImageData *image,double spatial_margin,double dose_ma int main(int argc,char * argv[]) { + clitk::RegisterClitkFactories(); + args_info_clitkGammaIndex args_info; if (cmdline_parser_clitkGammaIndex(argc, argv, &args_info) != 0) @@ -178,13 +226,18 @@ int main(int argc,char * argv[]) } // 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 @@ -199,33 +252,40 @@ int main(int argc,char * argv[]) 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; kkGetNumberOfPoints(); 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); @@ -242,27 +302,24 @@ int main(int argc,char * argv[]) 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(over_one)/total << endl; + cout << "ratio=" << static_cast(over_one)/kk << endl; } locator->Delete(); target->Delete(); saveImage(output,gamma_filename); - output->Delete(); return 0; } -