From 82d0db04054ac02ddbd921427ddc2edf73d475a0 Mon Sep 17 00:00:00 2001 From: dsarrut Date: Fri, 25 Mar 2011 09:30:51 +0000 Subject: [PATCH] clitkGammaIndex uses vtkCellLocator --- tools/CMakeLists.txt | 2 +- tools/clitkGammaIndex.cxx | 420 ++++++++++++++++---------------------- tools/clitkGammaIndex.ggo | 5 +- 3 files changed, 184 insertions(+), 243 deletions(-) diff --git a/tools/CMakeLists.txt b/tools/CMakeLists.txt index b9d2441..789241b 100644 --- a/tools/CMakeLists.txt +++ b/tools/CMakeLists.txt @@ -98,7 +98,7 @@ IF (CLITK_BUILD_TOOLS) WRAP_GGO(clitkGammaIndex_GGO_C clitkGammaIndex.ggo) ADD_EXECUTABLE(clitkGammaIndex clitkGammaIndex.cxx ${clitkGammaIndex_GGO_C}) - TARGET_LINK_LIBRARIES(clitkGammaIndex ITKBasicFilters ITKCommon ITKIO) + TARGET_LINK_LIBRARIES(clitkGammaIndex vtkCommon vtkFiltering vtkGraphics vtkIO vtkImaging) ADD_EXECUTABLE(clitkImageArithm clitkImageArithm.cxx ${clitkImageArithm_GGO_C}) TARGET_LINK_LIBRARIES(clitkImageArithm clitkImageArithmImageLib clitkCommon ITKIO ) diff --git a/tools/clitkGammaIndex.cxx b/tools/clitkGammaIndex.cxx index 39cbf6f..c784122 100644 --- a/tools/clitkGammaIndex.cxx +++ b/tools/clitkGammaIndex.cxx @@ -1,195 +1,148 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include #include -using std::cout; +#include +#include using std::endl; - -#include -#include -#include -#include -#include -#include -#include +using std::cout; #include "clitkGammaIndex_ggo.h" -const unsigned int image_dim = 2; - -typedef itk::Image Image; -typedef itk::ImageRegionIteratorWithIndex ImageIterator; -typedef itk::ImageRegionConstIteratorWithIndex ImageConstIterator; -typedef itk::ImageFileReader Reader; -typedef itk::MinimumMaximumImageCalculator MinMaxer; -typedef itk::ShiftScaleImageFilter Normalizer; -typedef itk::ChangeInformationImageFilter Scaler; - -typedef itk::Image ImageBin; -typedef itk::Image ImageMap; -typedef itk::SignedMaurerDistanceMapImageFilter Mapper; - -template -void SaveImage(const ImageType *image, const std::string &filename) { - typedef typename itk::ImageFileWriter Writer; - typename Writer::Pointer writer = Writer::New(); - writer->SetFileName(filename); +vtkImageData *loadImage(const std::string &filename) { + vtkImageReader2 *reader = vtkMetaImageReader::New(); + //vtkImageReader2 *reader = vtkPNGReader::New(); + reader->SetFileName(filename.c_str()); + reader->Update(); + + vtkImageData *image = reader->GetOutput(); + image->Register(NULL); + reader->Delete(); + + return 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->Update(); + writer->Write(); + writer->Delete(); } -ImageBin::Pointer AllocateImageBin(const Image *reference, const Image *target, unsigned int dose_size, float dose_max) { - const Image::RegionType reference_region = reference->GetLargestPossibleRegion(); - const Image::RegionType target_region = target->GetLargestPossibleRegion(); - const Image::SpacingType reference_spacing = reference->GetSpacing(); - const Image::SpacingType target_spacing = target->GetSpacing(); - const Image::PointType reference_origin = reference->GetOrigin(); - const Image::PointType target_origin = target->GetOrigin(); - - assert(reference_region == target_region); - assert(reference_spacing == target_spacing); - assert(reference_origin == target_origin); - - ImageBin::RegionType region; - ImageBin::SpacingType spacing; - ImageBin::PointType origin; - for (unsigned long kk=0; kkGetPoint(ids[0],p0); + double p1[3]; points->GetPoint(ids[1],p1); + double p2[3]; points->GetPoint(ids[2],p2); + double p3[3]; points->GetPoint(ids[3],p3); + //cout << "----------------------------------" << endl; + //cout << "p0=[" << p0[0] << "," << p0[1] << "," << p0[2] << "]" << endl; + //cout << "p1=[" << p1[0] << "," << p1[1] << "," << p1[2] << "]" << endl; + //cout << "p2=[" << p2[0] << "," << p2[1] << "," << p2[2] << "]" << endl; + //cout << "p3=[" << p3[0] << "," << p3[1] << "," << p3[2] << "]" << endl; + + double center[] = {0,0,0}; + for (int kk=0; kk<3; kk++) { + center[kk] += p0[kk]; + center[kk] += p1[kk]; + center[kk] += p2[kk]; + center[kk] += p3[kk]; + center[kk] /= 4.; } - const unsigned long dose_dimension = Image::GetImageDimension(); - region.SetSize(dose_dimension,dose_size); - region.SetIndex(dose_dimension,0); - spacing[dose_dimension] = dose_max/(dose_size-1); - origin[dose_dimension] = 0; - - ImageBin::Pointer map = ImageBin::New(); - map->SetRegions(region); - map->SetSpacing(spacing); - map->SetOrigin(origin); - map->Allocate(); - map->FillBuffer(0); - - return map; -} -Image::Pointer AllocateImageGamma(const Image *target) { - const Image::RegionType target_region = target->GetLargestPossibleRegion(); - const Image::SpacingType target_spacing = target->GetSpacing(); - const Image::PointType target_origin = target->GetOrigin(); + vtkIdType center_id = points->InsertNextPoint(center); + //cout << "center=[" << center[0] << "," << center[1] << "," << center[2] << "]" << endl; + + cells->InsertNextCell(3); + cells->InsertCellPoint(ids[0]); + cells->InsertCellPoint(ids[1]); + cells->InsertCellPoint(center_id); - Image::Pointer gamma = Image::New(); - gamma->SetRegions(target_region); - gamma->SetSpacing(target_spacing); - gamma->SetOrigin(target_origin); - gamma->Allocate(); - gamma->FillBuffer(0); + cells->InsertNextCell(3); + cells->InsertCellPoint(ids[1]); + cells->InsertCellPoint(ids[3]); + cells->InsertCellPoint(center_id); - return gamma; + cells->InsertNextCell(3); + cells->InsertCellPoint(ids[3]); + cells->InsertCellPoint(ids[2]); + cells->InsertCellPoint(center_id); + + cells->InsertNextCell(3); + cells->InsertCellPoint(ids[2]); + cells->InsertCellPoint(ids[0]); + cells->InsertCellPoint(center_id); } -void FillImageBin(ImageBin *image_bin, const Image *reference) { - ImageConstIterator iterator(reference,reference->GetLargestPossibleRegion()); - iterator.GoToBegin(); - while (!iterator.IsAtEnd()) { - Image::PixelType value = iterator.Get(); - Image::PointType point; - reference->TransformIndexToPhysicalPoint(iterator.GetIndex(),point); - - ImageBin::PointType point_bin; - for (unsigned long kk=0; kkTransformPhysicalPointToIndex(point_bin,index_bin); -#else - bool found = image_bin->TransformPhysicalPointToIndex(point_bin,index_bin); - assert(found); -#endif + for (int kk=0; kkGetNumberOfPoints(); kk++) { + double value = image->GetPointData()->GetScalars()->GetTuple1(kk); - while (index_bin[Image::GetImageDimension()] >=0 ) { - image_bin->SetPixel(index_bin,255); - index_bin[Image::GetImageDimension()]--; + if (first) { + maximum = value; + first = false; + continue; } - - ++iterator; + + if (maximumGetLargestPossibleRegion()); - ImageConstIterator target_iterator(target,target->GetLargestPossibleRegion()); - gamma_iterator.GoToBegin(); - target_iterator.GoToBegin(); - while (!target_iterator.IsAtEnd()) { - assert(!gamma_iterator.IsAtEnd()); - - Image::PixelType value = target_iterator.Get(); - Image::PointType point; - target->TransformIndexToPhysicalPoint(target_iterator.GetIndex(),point); - - ImageMap::PointType point_map; - point_map.Fill(0); - for (unsigned long kk=0; kkTransformPhysicalPointToIndex(point_map,index_map); +vtkPolyData *buildPlane(vtkImageData *image,double spatial_margin,double dose_margin) { + vtkPoints *points = vtkPoints::New(); + for (int kk=0; kkGetNumberOfPoints(); kk++) { + double *point = image->GetPoint(kk); + double value = image->GetPointData()->GetScalars()->GetTuple1(kk); + assert(value>=0); + assert(point[2]==0); + point[2] = value; + + point[0] /= spatial_margin; + point[1] /= spatial_margin; + point[2] /= dose_margin; + +#ifndef NDEBUG + vtkIdType point_id = points->InsertNextPoint(point); + assert(kk==point_id); #else - bool found = distance->TransformPhysicalPointToIndex(point_map,index_map); - assert(found); + points->InsertNextPoint(point); #endif - - gamma_iterator.Set(fabsf(distance->GetPixel(index_map))); - - ++gamma_iterator; - ++target_iterator; } -} -void TuneScaler(Scaler *scaler,float space_margin) { - Scaler::InputImageType::PointType origin = scaler->GetInput()->GetOrigin(); - Scaler::InputImageType::SpacingType spacing = scaler->GetInput()->GetSpacing(); - for (unsigned int kk=0; kkGetNumberOfCells(); kk++) { + vtkCell *cell = image->GetCell(kk); + assert(cell->GetNumberOfPoints()==4); + insertTriangles(cells,points,cell->GetPointIds()->GetPointer(0)); } - scaler->SetOutputSpacing(spacing); - scaler->SetOutputOrigin(origin); - scaler->ChangeSpacingOn(); - scaler->ChangeOriginOn(); -} + vtkPolyData *data = vtkPolyData::New(); + data->SetPoints(points); + data->SetPolys(cells); + points->Delete(); + cells->Delete(); -Image::PixelType GetImageMaximum(const Image *image) { - MinMaxer::Pointer minmaxer = MinMaxer::New(); - minmaxer->SetImage(image); - minmaxer->ComputeMaximum(); - return minmaxer->GetMaximum(); + return data; } -void ComputeGammaRatio(const Image *image) { - ImageConstIterator iterator(image,image->GetLargestPossibleRegion()); - iterator.GoToBegin(); - unsigned long total = 0; - unsigned long over_one = 0; - while (!iterator.IsAtEnd()) { - Image::PixelType value = iterator.Get(); - - if (value>1) over_one++; - total++; - - ++iterator; - } - - cout << "total=" << total << endl; - cout << "over_one=" << over_one << endl; - cout << "ratio=" << static_cast(over_one)/total << endl; +void assert2D(vtkImageData *image) { +#ifndef NDEBUG + int *extent = image->GetWholeExtent(); + assert(extent[4]==0); + assert(extent[5]==0); +#endif } int main(int argc,char * argv[]) @@ -204,109 +157,98 @@ int main(int argc,char * argv[]) exit(1); } - if (args_info.isodose_number_arg <= 0) { - std::cerr << "Specify a valid isodose number (>0)" << endl; - exit(1); - } - bool verbose = args_info.verbose_flag; std::string reference_filename(args_info.reference_arg); std::string target_filename(args_info.target_arg); std::string gamma_filename(args_info.output_arg); - float space_margin = args_info.spatial_margin_arg; - float dose_rel_margin = args_info.relative_dose_margin_arg; - float dose_margin = args_info.absolute_dose_margin_arg; + double space_margin = args_info.spatial_margin_arg; + double dose_rel_margin = args_info.relative_dose_margin_arg; + double dose_margin = args_info.absolute_dose_margin_arg; bool use_dose_margin = args_info.absolute_dose_margin_given; - unsigned int dose_size = args_info.isodose_number_arg; if (verbose) { cout << "reference_filename=" << reference_filename << endl; cout << "target_filename=" << target_filename << endl; cout << "gamma_filename=" << gamma_filename << endl; - cout << "dose_size=" << dose_size << endl; cout << "space_margin=" << space_margin << endl; if (use_dose_margin) cout << "dose_margin=" << dose_margin << endl; else cout << "dose_rel_margin=" << dose_rel_margin << endl; } - // load images - Reader::Pointer reference_reader = Reader::New(); - Reader::Pointer target_reader = Reader::New(); - { - reference_reader->SetFileName(reference_filename); - target_reader->SetFileName(target_filename); - reference_reader->Update(); - target_reader->Update(); - } + // load reference + vtkImageData *reference = loadImage(reference_filename); + assert2D(reference); // intensity normalisation if (!use_dose_margin) { - MinMaxer::PixelType reference_max = GetImageMaximum(reference_reader->GetOutput()); - //MinMaxer::PixelType target_max = GetImageMaximum(target_reader->GetOutput()); - - dose_margin = reference_max*dose_rel_margin; - + dose_margin = getMaximum(reference)*dose_rel_margin; if (verbose) cout << "dose_margin=" << dose_margin << endl; } - // 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; kkGetNumberOfPoints(); 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; - // normalize space coordinates - Scaler::Pointer reference_scaler = Scaler::New(); - Scaler::Pointer target_scaler = Scaler::New(); - { - reference_scaler->SetInput(reference_normalizer->GetOutput()); - TuneScaler(reference_scaler,space_margin); - reference_scaler->Update(); + double closest_point[3] = {0,0,0}; + vtkIdType cell_id = 0; + int foo = 0; + double squared_distance = 0; - target_scaler->SetInput(target_normalizer->GetOutput()); - TuneScaler(target_scaler,space_margin); - target_scaler->Update(); + 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++; - //SaveImage(reference_scaler->GetOutput(),"norm_reference.mhd"); - //SaveImage(reference_scaler->GetOutput(),"norm_target.mhd"); } - // compute hyper surface plane - float reference_dose_scaled_max = GetImageMaximum(reference_scaler->GetOutput()); - float target_dose_scaled_max = GetImageMaximum(target_scaler->GetOutput()); - float dose_scaled_max = reference_dose_scaled_max > target_dose_scaled_max ? reference_dose_scaled_max : target_dose_scaled_max; - ImageBin::Pointer image_bin = AllocateImageBin(reference_scaler->GetOutput(),target_scaler->GetOutput(),dose_size,dose_scaled_max); - FillImageBin(image_bin,reference_scaler->GetOutput()); - //SaveImage(image_bin.GetPointer(),"surface.mhd"); - - // compute distance map - Mapper::Pointer mapper = Mapper::New(); - mapper->InsideIsPositiveOn(); - mapper->SquaredDistanceOff(); - mapper->UseImageSpacingOn(); - mapper->SetInput(image_bin); - mapper->Update(); - //SaveImage(mapper->GetOutput(),"distance.mhd"); - - // extract gamma index from distance map - Image::Pointer image_gamma = AllocateImageGamma(target_normalizer->GetOutput()); - FillImageGamma(image_gamma,target_scaler->GetOutput(),mapper->GetOutput()); - SaveImage(image_gamma.GetPointer(),gamma_filename); - - if (verbose) ComputeGammaRatio(image_gamma); + if (verbose) { + cout << "total=" << total << endl; + cout << "over_one=" << over_one << endl; + cout << "ratio=" << static_cast(over_one)/total << endl; + } + + locator->Delete(); + target->Delete(); + + saveImage(output,gamma_filename); + output->Delete(); return 0; } diff --git a/tools/clitkGammaIndex.ggo b/tools/clitkGammaIndex.ggo index 1e6870d..20c87f9 100644 --- a/tools/clitkGammaIndex.ggo +++ b/tools/clitkGammaIndex.ggo @@ -1,7 +1,7 @@ #File clitkGammaIndex.ggo package "clitkGammaIndex" -version "1.0" -purpose "Compute the gamma index map for target dose map vs reference dose map. If verbose is given, the ratio of gamma>1 pixels of the total pixel in the image is computed. Absolute dose margin has priority over relative dose margin. The relative dose margin is relative to the maximum dose in the reference image. Reference and target images must share the same origin, the same spacing and the same size." +version "1.1" +purpose "Compute the gamma index map for target dose map vs reference dose map. If verbose is given, the ratio of gamma>1 pixels of the total pixel in the image is computed. Absolute dose margin has priority over relative dose margin. The relative dose margin is relative to the maximum dose in the reference image.\n\nclitkGammaIndex -v -i ref.mhd -j disp.mhd -o gamma.mhd -s 3 -r 0.04\nthis will compute the 3mm 4% gamma index between ref and disp." option "verbose" v "Verbose" flag off @@ -12,5 +12,4 @@ option "output" o "Output image filename" string yes option "spatial-margin" s "Spatial margin [mm]" double yes option "relative-dose-margin" r "Dose margin relative to max dose in reference [%]" double no option "absolute-dose-margin" d "Absolute dose margin [Gray]" double no -option "isodose-number" n "Number of isodose computed in distance map" int default="128" no -- 2.46.0