From 18fd13ca69bbe5a6d91e070fef22ca9c64c92063 Mon Sep 17 00:00:00 2001 From: delmon Date: Wed, 23 Mar 2011 13:38:26 +0000 Subject: [PATCH] pgueth: + Add clitkGammaIndex computation. --- tools/CMakeLists.txt | 4 + tools/clitkGammaIndex.cxx | 308 ++++++++++++++++++++++++++++++++++++++ tools/clitkGammaIndex.ggo | 16 ++ 3 files changed, 328 insertions(+) create mode 100644 tools/clitkGammaIndex.cxx create mode 100644 tools/clitkGammaIndex.ggo diff --git a/tools/CMakeLists.txt b/tools/CMakeLists.txt index ac74780..b9d2441 100644 --- a/tools/CMakeLists.txt +++ b/tools/CMakeLists.txt @@ -96,6 +96,10 @@ IF (CLITK_BUILD_TOOLS) ADD_EXECUTABLE(clitkGuerreroVentilation clitkGuerreroVentilation.cxx clitkGuerreroVentilationGenericFilter.cxx ${clitkGuerreroVentilation_GGO_C}) TARGET_LINK_LIBRARIES(clitkGuerreroVentilation clitkCommon ITKIO) + WRAP_GGO(clitkGammaIndex_GGO_C clitkGammaIndex.ggo) + ADD_EXECUTABLE(clitkGammaIndex clitkGammaIndex.cxx ${clitkGammaIndex_GGO_C}) + TARGET_LINK_LIBRARIES(clitkGammaIndex ITKBasicFilters ITKCommon ITKIO) + 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 new file mode 100644 index 0000000..d6dfaa7 --- /dev/null +++ b/tools/clitkGammaIndex.cxx @@ -0,0 +1,308 @@ +#include +using std::cout; +using std::endl; + +#include +#include +#include +#include +#include +#include +#include + +#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::DanielssonDistanceMapImageFilter 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); + writer->SetInput(image); + writer->Update(); +} + +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; kkSetRegions(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(); + + Image::Pointer gamma = Image::New(); + gamma->SetRegions(target_region); + gamma->SetSpacing(target_spacing); + gamma->SetOrigin(target_origin); + gamma->Allocate(); + gamma->FillBuffer(0); + + return gamma; +} + +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 + + image_bin->SetPixel(index_bin,255); + + ++iterator; + } +} + +void FillImageGamma(Image *gamma, const Image *target, const ImageMap *distance) { + ImageIterator gamma_iterator(gamma,gamma->GetLargestPossibleRegion()); + 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); +#else + bool found = distance->TransformPhysicalPointToIndex(point_map,index_map); + assert(found); +#endif + + gamma_iterator.Set(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; kkSetOutputSpacing(spacing); + scaler->SetOutputOrigin(origin); + scaler->ChangeSpacingOn(); + scaler->ChangeOriginOn(); +} + +Image::PixelType GetImageMaximum(const Image *image) { + MinMaxer::Pointer minmaxer = MinMaxer::New(); + minmaxer->SetImage(image); + minmaxer->ComputeMaximum(); + return minmaxer->GetMaximum(); +} + +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; +} + +int main(int argc,char * argv[]) +{ + args_info_clitkGammaIndex args_info; + + if (cmdline_parser_clitkGammaIndex(argc, argv, &args_info) != 0) + exit(1); + + if (!args_info.absolute_dose_margin_given && !args_info.relative_dose_margin_given) { + std::cerr << "Specify either relative or absolute dose margin" << endl; + 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; + 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(); + } + + // 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; + + 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; + } + + // 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(); + + target_scaler->SetInput(target_normalizer->GetOutput()); + TuneScaler(target_scaler,space_margin); + target_scaler->Update(); + + //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->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); + + return 0; +} + diff --git a/tools/clitkGammaIndex.ggo b/tools/clitkGammaIndex.ggo new file mode 100644 index 0000000..1e6870d --- /dev/null +++ b/tools/clitkGammaIndex.ggo @@ -0,0 +1,16 @@ +#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." + +option "verbose" v "Verbose" flag off + +option "reference" i "Reference dose map" string yes +option "target" j "Target dose map" string yes +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.47.1