diff --git a/Modules/DiffusionImaging/MiniApps/HotspotMiniApp.cpp b/Modules/DiffusionImaging/MiniApps/HotspotMiniApp.cpp index 2f76e385e9..c9b3a74f69 100644 --- a/Modules/DiffusionImaging/MiniApps/HotspotMiniApp.cpp +++ b/Modules/DiffusionImaging/MiniApps/HotspotMiniApp.cpp @@ -1,104 +1,103 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkCommandLineParser.h" #include "mitkImage.h" #include "mitkIOUtil.h" #include #include #include #include int main( int argc, char* argv[] ) { mitkCommandLineParser parser; parser.setTitle("Test basic hotspot functionality"); parser.setCategory("Preprocessing Tools"); parser.setDescription(""); parser.setContributor("MBI"); parser.setArgumentPrefix("--", "-"); parser.addArgument("help", "h", mitkCommandLineParser::String, "Help:", "Show this help text"); parser.addArgument("input", "i", mitkCommandLineParser::InputFile, "Input:", "input image", us::Any(),false); parser.addArgument("mask", "m", mitkCommandLineParser::InputFile, "Mask:", "mask image / roi image denotin area on which statistics are calculated", us::Any(),true); parser.addArgument("out", "o", mitkCommandLineParser::OutputFile, "Output", "output file (default: hotspotMiniApp_output.nrrd)", us::Any()); std::cout << "test...." << std::endl; std::map parsedArgs = parser.parseArguments(argc, argv); std::cout << "parsedArgs.size()= " << parsedArgs.size() << std::endl; if (parsedArgs.size()==0 || parsedArgs.count("help") || parsedArgs.count("h")) { // std::cout << "\n\n MiniApp Description: \nCalculates and saves the hotspot of an image." << endl; // std::cout << "Output is written to the designated output file" << endl; // std::cout << parser.helpText(); // return EXIT_SUCCESS; } // Parameters: std::string inputImageFile; if (!parsedArgs.count("i") && !parsedArgs.count("input")) { inputImageFile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic3D.nrrd"; } else { inputImageFile = us::any_cast(parsedArgs["input"]); } mitk::Image::Pointer maskImage; if (parsedArgs.count("mask") || parsedArgs.count("m")) { std::string maskImageFile = us::any_cast(parsedArgs["mask"]); maskImage = mitk::IOUtil::LoadImage(maskImageFile); } std::string outFile; if (parsedArgs.count("out") || parsedArgs.count("o") ) outFile = us::any_cast(parsedArgs["out"]); else outFile = "hotspotMiniApp_output.nrrd"; // Load image and mask mitk::Image::Pointer inputImage = mitk::IOUtil::LoadImage(inputImageFile); // Calculate statistics - mitk::HotspotCalculator hotspotCalc; - hotspotCalc.setInputImage(inputImage); - hotspotCalc.setMaskImage(maskImage); - hotspotCalc.setHotspotRadiusInMM(10); - + mitk::HotspotMaskGenerator::Pointer hotspotCalc = mitk::HotspotMaskGenerator::New(); + hotspotCalc->SetImage(inputImage); + hotspotCalc->SetHotspotRadiusInMM(10); + hotspotCalc->SetTimeStep(0); mitk::Image::Pointer outImage; try { - outImage = hotspotCalc.getHotspotMask(0, 0); + outImage = hotspotCalc->GetMask(); } catch( const itk::ExceptionObject& e) { MITK_ERROR << "Failed - ITK Exception:" << e.what(); return -1; } if (outImage != nullptr) { mitk::IOUtil::SaveImage(outImage, outFile); } return EXIT_SUCCESS; } diff --git a/Modules/DiffusionImaging/MiniApps/ImageStatisticsMiniapp_v2.cpp b/Modules/DiffusionImaging/MiniApps/ImageStatisticsMiniapp_v2.cpp index 4a3c4af7de..5fa25988cb 100644 --- a/Modules/DiffusionImaging/MiniApps/ImageStatisticsMiniapp_v2.cpp +++ b/Modules/DiffusionImaging/MiniApps/ImageStatisticsMiniapp_v2.cpp @@ -1,202 +1,199 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkCommandLineParser.h" #include "mitkImage.h" +#include #include +#include +#include +#include #include "mitkIOUtil.h" #include #include #include #include #include "mitkImageAccessByItk.h" #include #include #include #include #include #include #include #include +#include +#include struct statistics_res{ double mean, variance, min, max, count, moment; }; void printstats(statistics_res s) { std::cout << "mean: " << s.mean << std::endl << "variance: " << s.variance << std::endl << "min: " << s.min << std::endl << "max: " << s.max << std::endl << "count: " << s.count << std::endl << "moment: " << s.moment << std::endl; } template void printMap(std::map input) { for (auto it = input.begin(); it != input.end(); ++it) { std::cout << it->first<< ": " << it->second<< std::endl; } } template < typename TPixel, unsigned int VImageDimension > void get_statistics_boost(itk::Image* itkImage, statistics_res& res){ typedef itk::Image ImageType; itk::ImageRegionConstIterator it(itkImage, itkImage->GetLargestPossibleRegion()); - TPixel currentPixel; int ctr=0; double sum=0; boost::accumulators::accumulator_set> > acc; for (it.GoToBegin(); !it.IsAtEnd(); ++it) { acc(it.Get()); // currentPixel = it.Get(); // sum+=currentPixel; // ctr+=1; } // res.mean=(double)sum/(double)ctr; res.mean = boost::accumulators::mean(acc); res.variance = boost::accumulators::variance(acc); res.min = boost::accumulators::min(acc); res.max = boost::accumulators::max(acc); res.count = boost::accumulators::count(acc); res.moment = boost::accumulators::moment<2>(acc); std::cout << "sum: " << sum << " N: " << ctr << " mean: " << res.mean << std::endl; } int main( int argc, char* argv[] ) { - mitkCommandLineParser parser; - - parser.setTitle("Extract Image Statistics"); - parser.setCategory("Preprocessing Tools"); - parser.setDescription(""); - parser.setContributor("MBI"); - - parser.setArgumentPrefix("--", "-"); - parser.addArgument("help", "h", mitkCommandLineParser::String, "Help:", "Show this help text"); - parser.addArgument("input", "i", mitkCommandLineParser::InputFile, "Input:", "input image", us::Any(),false); - parser.addArgument("mask", "m", mitkCommandLineParser::InputFile, "Mask:", "mask image / roi image denotin area on which statistics are calculated", us::Any(),true); - parser.addArgument("out", "o", mitkCommandLineParser::OutputFile, "Output", "output file (default: filenameOfRoi.nrrd_statistics.txt)", us::Any()); - - std::cout << "test...." << std::endl; - - std::map parsedArgs = parser.parseArguments(argc, argv); - std::cout << "parsedArgs.size()= " << parsedArgs.size() << std::endl; - if (parsedArgs.count("help") || parsedArgs.count("h")) - { - std::cout << "\n\n MiniApp Description: \nCalculates statistics on the supplied image using given mask." << endl; - std::cout << "Output is written to the designated output file in this order:" << endl; - std::cout << "Mean, Standard Deviation, RMS, Max, Min, Number of Voxels, Volume [mm3]" << endl; - std::cout << "\n\n Parameters:"<< endl; - std::cout << parser.helpText(); - return EXIT_SUCCESS; - } - - - // Parameters: unsigned int timeStep = 0; std::string inputImageFile; - if (parsedArgs.count("input") || parsedArgs.count("i")) - { - inputImageFile = us::any_cast(parsedArgs["input"]); - } else - { - inputImageFile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic3D.nrrd"; - } + inputImageFile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic3D.nrrd"; + // Load image + mitk::Image::Pointer inputImage = mitk::IOUtil::LoadImage(inputImageFile); mitk::Image::Pointer maskImage; - if (parsedArgs.count("mask") || parsedArgs.count("m")) - { - std::string maskImageFile = us::any_cast(parsedArgs["mask"]); - maskImage = mitk::IOUtil::LoadImage(maskImageFile); - } - - std::string outFile; - if (parsedArgs.count("out") || parsedArgs.count("o") ) - outFile = us::any_cast(parsedArgs["out"]); - else - outFile = inputImageFile + "_statistics.txt"; - - // Load image and mask - mitk::Image::Pointer inputImage = mitk::IOUtil::LoadImage(inputImageFile); + std::string maskImageFile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic3D_someSegmentation.nrrd"; + maskImage = mitk::IOUtil::LoadImage(maskImageFile); + mitk::BinaryImageMaskGenerator::Pointer binaryImageMaskGen = mitk::BinaryImageMaskGenerator::New(); + binaryImageMaskGen->SetImageMask(maskImage); + + mitk::HotspotMaskGenerator::Pointer hotSpotMaskGen = mitk::HotspotMaskGenerator::New(); + hotSpotMaskGen->SetImage(inputImage); + hotSpotMaskGen->SetHotspotRadiusInMM(10); + hotSpotMaskGen->SetTimeStep(0); + + std::string pFfile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic3D_rectangle.pf"; + std::vector loadedObjects; + loadedObjects = mitk::IOUtil::Load(pFfile); + mitk::BaseData::Pointer pf = loadedObjects[0]; + mitk::PlanarFigure::Pointer planarFigure = dynamic_cast(pf.GetPointer()); + mitk::PlanarFigureMaskGenerator::Pointer planarFigMaskExtr = mitk::PlanarFigureMaskGenerator::New(); + planarFigMaskExtr->SetPlanarFigure(planarFigure); + planarFigMaskExtr->SetImage(inputImage); // Calculate statistics mitk::ImageStatisticsCalculator_v2::Pointer calculator = mitk::ImageStatisticsCalculator_v2::New(); std::cout << "calculating statistics (unmasked) itk: " << std::endl; mitk::ImageStatisticsCalculator_v2::statisticsMapType result; + + calculator->SetInputImage(inputImage); + calculator->SetMask(planarFigMaskExtr.GetPointer()); + calculator->SetDoHistogramStatistics(true); + try { - calculator->SetInputImage(inputImage); - result = calculator->GetStatistics(timeStep); + result = calculator->GetStatistics(timeStep, 1); } catch( const itk::ExceptionObject& e) { MITK_ERROR << "Statistic Calculation Failed - ITK Exception:" << e.what(); return -1; } - // Calculate Volume - double volume = 0; - const mitk::BaseGeometry *geometry = inputImage->GetGeometry(); - if ( geometry != NULL ) - { - const mitk::Vector3D &spacing = inputImage->GetGeometry()->GetSpacing(); - volume = spacing[0] * spacing[1] * spacing[2] * (double) result["N"]; - } - - // Write Results to file - std::ofstream output; - output.open(outFile.c_str()); - output << result["Mean"] << " , "; - output << result["Std"] << " , "; - output << result["RMS"] << " , "; - output << result["Max"] << " , "; - output << result["Min"] << " , "; - output << result["N"] << " , "; - output << volume << "\n"; - - output.flush(); - output.close(); + printMap(result); + calculator->SetMask(binaryImageMaskGen.GetPointer()); + result = calculator->GetStatistics(timeStep, 1); printMap(result); - std::cout << "calculating statistics boost: " << std::endl; + calculator->SetMask(nullptr); + result = calculator->GetStatistics(timeStep); + printMap(result); - statistics_res res; - AccessByItk_n(inputImage, get_statistics_boost, (res)); + calculator->SetMask(hotSpotMaskGen.GetPointer()); + result = calculator->GetStatistics(timeStep, 1); + printMap(result); - printstats(res); +// std::cout << std::endl << std::endl << "calculating statistics old imgstatcalc: " << std::endl; + +// mitk::ImageStatisticsCalculator::Statistics statisticsStruct; +// mitk::ImageStatisticsCalculator::Pointer calculator_old = mitk::ImageStatisticsCalculator::New(); +// calculator_old->SetImage(inputImage); +// calculator_old->SetImageMask(maskImage); +// //calculator_old->SetMaskingModeToImage(); +// calculator_old->SetHistogramBinSize(100); +// calculator_old->SetPlanarFigure(planarFigure); +// calculator_old->SetMaskingModeToPlanarFigure(); +// calculator_old->SetDoIgnorePixelValue(false); +// calculator_old->ComputeStatistics(timeStep); +// statisticsStruct = calculator_old->GetStatistics(timeStep); + +// std::cout << "Entropy: " << statisticsStruct.GetEntropy() << std::endl; +// std::cout << "Kurtosis: " << statisticsStruct.GetKurtosis() << std::endl; +// std::cout << "MPP: " << statisticsStruct.GetMPP() << std::endl; +// std::cout << "Max: " << statisticsStruct.GetMax() << std::endl; +// std::cout << "Mean: " << statisticsStruct.GetMean() << std::endl; +// std::cout << "Median: " << statisticsStruct.GetMedian() << std::endl; +// std::cout << "Min: " << statisticsStruct.GetMin() << std::endl; +// std::cout << "N: " << statisticsStruct.GetN() << std::endl; +// std::cout << "RMS: " << statisticsStruct.GetRMS() << std::endl; +// std::cout << "Skewness: " << statisticsStruct.GetSkewness() << std::endl; +// std::cout << "Std: " << statisticsStruct.GetSigma() << std::endl; +// std::cout << "UPP: " << statisticsStruct.GetUPP() << std::endl; +// std::cout << "Uniformity: " << statisticsStruct.GetUniformity() << std::endl; +// std::cout << "Variance: " << statisticsStruct.GetVariance() << std::endl << std::endl; + + +// std::cout << "calculating statistics boost: " << std::endl; +// statistics_res res; +// AccessByItk_n(inputImage, get_statistics_boost, (res)); +// printstats(res); return EXIT_SUCCESS; } diff --git a/Modules/DiffusionImaging/MiniApps/PlanarFigureMaskCreatorMiniapp.cpp b/Modules/DiffusionImaging/MiniApps/PlanarFigureMaskCreatorMiniapp.cpp index c099bdf1e9..585a20d75a 100644 --- a/Modules/DiffusionImaging/MiniApps/PlanarFigureMaskCreatorMiniapp.cpp +++ b/Modules/DiffusionImaging/MiniApps/PlanarFigureMaskCreatorMiniapp.cpp @@ -1,131 +1,131 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkCommandLineParser.h" #include "mitkImage.h" #include "mitkIOUtil.h" #include #include #include #include #include #include int main( int argc, char* argv[] ) { mitkCommandLineParser parser; parser.setTitle("Test basic hotspot functionality"); parser.setCategory("Preprocessing Tools"); parser.setDescription(""); parser.setContributor("MBI"); parser.setArgumentPrefix("--", "-"); parser.addArgument("help", "h", mitkCommandLineParser::String, "Help:", "Show this help text"); parser.addArgument("input", "i", mitkCommandLineParser::InputFile, "Input:", "input image", us::Any(),false); parser.addArgument("planarfigure", "p", mitkCommandLineParser::InputFile, "PlanarFigure:", "mask image / roi image denotin area on which statistics are calculated", us::Any(),false); parser.addArgument("out", "o", mitkCommandLineParser::OutputFile, "Output", "output file (default: hotspotMiniApp_output.nrrd)", us::Any()); std::cout << "test...." << std::endl; std::map parsedArgs = parser.parseArguments(argc, argv); std::cout << "parsedArgs.size()= " << parsedArgs.size() << std::endl; if (parsedArgs.size()==0 || parsedArgs.count("help") || parsedArgs.count("h")) { // std::cout << "\n\n MiniApp Description: \nCalculates and saves the hotspot of an image." << endl; // std::cout << "Output is written to the designated output file" << endl; // std::cout << parser.helpText(); // return EXIT_SUCCESS; } // Parameters: std::string inputImageFile; if (!parsedArgs.count("i") && !parsedArgs.count("input")) { inputImageFile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic3D.nrrd"; } else { inputImageFile = us::any_cast(parsedArgs["input"]); } std::string inputPFFile; if (parsedArgs.count("planarfigure") || parsedArgs.count("p")) { inputPFFile = us::any_cast(parsedArgs["planarfigure"]); } else { inputPFFile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic3D_rectangle.pf"; } mitk::PlanarFigure::Pointer planarFigure; try { std::vector loadedObjects; // try except loadedObjects = mitk::IOUtil::Load(inputPFFile); mitk::BaseData::Pointer pf = loadedObjects[0]; planarFigure = dynamic_cast(pf.GetPointer()); if (planarFigure.IsNull()) { MITK_ERROR << "something went wrong"; return -1; } } catch (const itk::ExceptionObject& e) { MITK_ERROR << "Failed: " << e.what(); return -1; } std::string outFile; if (parsedArgs.count("out") || parsedArgs.count("o") ) outFile = us::any_cast(parsedArgs["out"]); else outFile = "planarFigureExtraction_output.nrrd"; // Load image and mask mitk::Image::Pointer inputImage = mitk::IOUtil::LoadImage(inputImageFile); // Calculate statistics - mitk::PlanarFigureMaskGenerator planarFigMaskExtr; - planarFigMaskExtr.SetPlanarFigure(planarFigure); - planarFigMaskExtr.SetImage(inputImage); + mitk::PlanarFigureMaskGenerator::Pointer planarFigMaskExtr = mitk::PlanarFigureMaskGenerator::New(); + planarFigMaskExtr->SetPlanarFigure(planarFigure); + planarFigMaskExtr->SetImage(inputImage); mitk::Image::Pointer outImage; try { - outImage = planarFigMaskExtr.GetMask(); + outImage = planarFigMaskExtr->GetMask(); } catch( const itk::ExceptionObject& e) { MITK_ERROR << "Failed - ITK Exception:" << e.what(); return -1; } mitk::BaseGeometry *imageGeo = outImage->GetGeometry(); std::cout << "origin: " << imageGeo->GetOrigin()[0] << " " << imageGeo->GetOrigin()[1] << " " << imageGeo->GetOrigin()[2] << std::endl; if (outImage != nullptr) { mitk::IOUtil::SaveImage(outImage, outFile); } return EXIT_SUCCESS; } diff --git a/Modules/GraphAlgorithms/CMakeLists.txt b/Modules/GraphAlgorithms/CMakeLists.txt index d5f4be6021..64c609e707 100644 --- a/Modules/GraphAlgorithms/CMakeLists.txt +++ b/Modules/GraphAlgorithms/CMakeLists.txt @@ -1,4 +1,4 @@ MITK_CREATE_MODULE( - DEPENDS MitkImageStatistics +# DEPENDS MitkImageStatistics WARNINGS_AS_ERRORS ) diff --git a/Modules/ImageStatistics/mitkBinaryImageMaskGenerator.cpp b/Modules/ImageStatistics/mitkBinaryImageMaskGenerator.cpp index 6713e02d75..84ac7420a9 100644 --- a/Modules/ImageStatistics/mitkBinaryImageMaskGenerator.cpp +++ b/Modules/ImageStatistics/mitkBinaryImageMaskGenerator.cpp @@ -1,43 +1,40 @@ #include #include #include namespace mitk { void BinaryImageMaskGenerator::SetImageMask(Image::Pointer maskImage) { if (m_internalMaskImage != maskImage) { m_internalMaskImage = maskImage; m_Modified = true; } } mitk::Image::Pointer BinaryImageMaskGenerator::GetMask() { if (m_Modified) { // extract timeStep -> is that correct? - if (m_internalMaskImage->GetTimeSteps() > 1) + if (m_TimeStep >= m_internalMaskImage->GetTimeSteps()) { - if (m_TimeStep >= m_internalMaskImage->GetTimeSteps()) - { - throw std::runtime_error("Error: Invalid time step in BinaryImageMaskGenerator!"); - } - ImageTimeSelector::Pointer imageTimeSelector = ImageTimeSelector::New(); - imageTimeSelector->SetInput(m_internalMaskImage); - imageTimeSelector->SetTimeNr(m_TimeStep); - imageTimeSelector->UpdateLargestPossibleRegion(); - - m_InternalMask = mitk::Image::New(); - m_InternalMask = imageTimeSelector->GetOutput(); - m_Modified = false; + throw std::runtime_error("Error: Invalid time step in BinaryImageMaskGenerator!"); } + ImageTimeSelector::Pointer imageTimeSelector = ImageTimeSelector::New(); + imageTimeSelector->SetInput(m_internalMaskImage); + imageTimeSelector->SetTimeNr(m_TimeStep); + imageTimeSelector->UpdateLargestPossibleRegion(); + + m_InternalMask = mitk::Image::New(); + m_InternalMask = imageTimeSelector->GetOutput(); + m_Modified = false; } return m_InternalMask; } } diff --git a/Modules/ImageStatistics/mitkBinaryImageMaskGenerator.h b/Modules/ImageStatistics/mitkBinaryImageMaskGenerator.h index c1fbceaaa7..838377f5e1 100644 --- a/Modules/ImageStatistics/mitkBinaryImageMaskGenerator.h +++ b/Modules/ImageStatistics/mitkBinaryImageMaskGenerator.h @@ -1,27 +1,42 @@ #ifndef mitkBinaryMaskGenerator #define mitkBinaryMaskGenerator #include #include #include +#include +#include namespace mitk { -class MITKIMAGESTATISTICS_EXPORT BinaryImageMaskGenerator: MaskGenerator +class MITKIMAGESTATISTICS_EXPORT BinaryImageMaskGenerator: public MaskGenerator { public: + /** Standard Self typedef */ + typedef BinaryImageMaskGenerator Self; + typedef MaskGenerator Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self) + + /** Runtime information support. */ + itkTypeMacro(BinaryImageMaskGenerator, MaskGenerator) + mitk::Image::Pointer GetMask(); void SetImageMask(mitk::Image::Pointer maskImage); protected: + BinaryImageMaskGenerator():Superclass(){} private: mitk::Image::Pointer m_internalMaskImage; }; } #endif diff --git a/Modules/ImageStatistics/mitkHotspotCalculator.cpp b/Modules/ImageStatistics/mitkHotspotCalculator.cpp index e385b36903..84b8205616 100644 --- a/Modules/ImageStatistics/mitkHotspotCalculator.cpp +++ b/Modules/ImageStatistics/mitkHotspotCalculator.cpp @@ -1,524 +1,535 @@ #include #include #include #include #include #include "mitkImageAccessByItk.h" #include #include #include namespace mitk { - HotspotCalculator::HotspotCalculator(): + HotspotMaskGenerator::HotspotMaskGenerator(): m_HotspotRadiusinMM(6.2035049089940), // radius of a 1cm3 sphere in mm m_HotspotMustBeCompletelyInsideImage(true), - m_HotspotParamsChanged(true) + m_Label(1) { + m_Modified=true; } - void HotspotCalculator::setInputImage(mitk::Image::Pointer inputImage) + void HotspotMaskGenerator::SetImage(mitk::Image::Pointer inputImage) { - m_Image = inputImage; + if (inputImage != m_Image) + { + m_Image = inputImage; + m_Modified = true; + } } - void HotspotCalculator::setMaskImage(mitk::Image::Pointer maskImage) + void HotspotMaskGenerator::SetMask(MaskGenerator::Pointer mask) { - m_Mask = maskImage; + if (mask != m_Mask) + { + m_Mask = mask; + m_Modified = true; + } } - HotspotCalculator::~HotspotCalculator() + HotspotMaskGenerator::~HotspotMaskGenerator() { } - void HotspotCalculator::setHotspotRadiusInMM(double radiusInMillimeter) + void HotspotMaskGenerator::SetHotspotRadiusInMM(double radiusInMillimeter) { if(radiusInMillimeter != m_HotspotRadiusinMM) { m_HotspotRadiusinMM = radiusInMillimeter; - m_HotspotParamsChanged = true; + m_Modified = true; } } - const double& HotspotCalculator::getHotspotRadiusinMM() const + const double& HotspotMaskGenerator::GetHotspotRadiusinMM() const { return m_HotspotRadiusinMM; } - bool HotspotCalculator::getHotspotMustBeCompletelyInsideImage() const + bool HotspotMaskGenerator::GetHotspotMustBeCompletelyInsideImage() const { return m_HotspotMustBeCompletelyInsideImage; } - mitk::Image::Pointer HotspotCalculator::getHotspotMask(const unsigned int timeStep, - unsigned int label) + mitk::Image::Pointer HotspotMaskGenerator::GetMask() { - if ( m_Image.IsNull() ) - { - throw std::runtime_error( "Error: image empty!" ); - } - - if ( timeStep >= m_Image->GetTimeSteps() ) - { - throw std::runtime_error( "Error: invalid time step!" ); - } - - mitk::ImageTimeSelector::Pointer imageTimeSelector = mitk::ImageTimeSelector::New(); - imageTimeSelector->SetInput( m_Image ); - imageTimeSelector->SetTimeNr( timeStep ); - imageTimeSelector->UpdateLargestPossibleRegion(); - mitk::Image::Pointer timeSliceImage = imageTimeSelector->GetOutput(); - - m_internalImage = timeSliceImage; - m_internalMask2D = nullptr; // is this correct when this variable holds a smart pointer? - m_internalMask3D = nullptr; - - if ( m_Mask != nullptr ) + if (m_Modified) { - unsigned int timeStep_mask = timeStep; - if ( timeStep >= m_Mask->GetTimeSteps() ) + if ( m_Image.IsNull() ) { - timeStep_mask = m_Mask->GetTimeSteps(); + throw std::runtime_error( "Error: image empty!" ); } - mitk::ImageTimeSelector::Pointer maskTimeSelector = mitk::ImageTimeSelector::New(); - maskTimeSelector->SetInput( m_Mask ); - maskTimeSelector->SetTimeNr( timeStep_mask ); - maskTimeSelector->UpdateLargestPossibleRegion(); - mitk::Image::Pointer timeSliceMask = maskTimeSelector->GetOutput(); - - if ( m_internalImage->GetDimension() == 3 ) - { - CastToItkImage(timeSliceMask, m_internalMask3D); - AccessFixedDimensionByItk_2(m_internalImage, CalculateHotspotMask, 3, m_internalMask3D, label); - } - else if ( m_internalImage->GetDimension() == 2 ) + if ( m_TimeStep >= m_Image->GetTimeSteps() ) { - CastToItkImage(timeSliceMask, m_internalMask2D); - AccessFixedDimensionByItk_2(m_internalImage, CalculateHotspotMask, 2, m_internalMask2D, label); + throw std::runtime_error( "Error: invalid time step!" ); } - else - { - throw std::runtime_error( "Error: invalid image dimension" ); - } - } - else - { - if ( m_internalImage->GetDimension() == 3 ) - { - AccessFixedDimensionByItk_2(m_internalImage, CalculateHotspotMask, 3, m_internalMask3D, label); - } - else if ( m_internalImage->GetDimension() == 2 ) + mitk::ImageTimeSelector::Pointer imageTimeSelector = mitk::ImageTimeSelector::New(); + imageTimeSelector->SetInput( m_Image ); + imageTimeSelector->SetTimeNr( m_TimeStep ); + imageTimeSelector->UpdateLargestPossibleRegion(); + mitk::Image::Pointer timeSliceImage = imageTimeSelector->GetOutput(); + + m_internalImage = timeSliceImage; + m_internalMask2D = nullptr; // is this correct when this variable holds a smart pointer? + m_internalMask3D = nullptr; + + if ( m_Mask != nullptr ) { - AccessFixedDimensionByItk_2(m_internalImage, CalculateHotspotMask, 2, m_internalMask2D, label); + m_Mask->SetTimeStep(m_TimeStep); + mitk::Image::Pointer timeSliceMask = m_Mask->GetMask(); + + if ( m_internalImage->GetDimension() == 3 ) + { + CastToItkImage(timeSliceMask, m_internalMask3D); + AccessFixedDimensionByItk_2(m_internalImage, CalculateHotspotMask, 3, m_internalMask3D, m_Label); + } + else if ( m_internalImage->GetDimension() == 2 ) + { + CastToItkImage(timeSliceMask, m_internalMask2D); + AccessFixedDimensionByItk_2(m_internalImage, CalculateHotspotMask, 2, m_internalMask2D, m_Label); + } + else + { + throw std::runtime_error( "Error: invalid image dimension" ); + } } else { - throw std::runtime_error( "Error: invalid image dimension" ); + + if ( m_internalImage->GetDimension() == 3 ) + { + AccessFixedDimensionByItk_2(m_internalImage, CalculateHotspotMask, 3, m_internalMask3D, m_Label); + } + else if ( m_internalImage->GetDimension() == 2 ) + { + AccessFixedDimensionByItk_2(m_internalImage, CalculateHotspotMask, 2, m_internalMask2D, m_Label); + } + else + { + throw std::runtime_error( "Error: invalid image dimension" ); + } } + m_Modified = false; } - return m_hotspotMask; + return m_InternalMask; + } + + void HotspotMaskGenerator::SetLabel(unsigned int label) + { + if (label != m_Label) + { + m_Label = label; + m_Modified = true; + } } template - HotspotCalculator::ImageExtrema - HotspotCalculator::CalculateExtremaWorld( const itk::Image* inputImage, + HotspotMaskGenerator::ImageExtrema + HotspotMaskGenerator::CalculateExtremaWorld( const itk::Image* inputImage, typename itk::Image::Pointer maskImage, double neccessaryDistanceToImageBorderInMM, unsigned int label ) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< unsigned short, VImageDimension > MaskImageType; typedef itk::ImageRegionConstIteratorWithIndex MaskImageIteratorType; typedef itk::ImageRegionConstIteratorWithIndex InputImageIndexIteratorType; typename ImageType::SpacingType spacing = inputImage->GetSpacing(); ImageExtrema minMax; minMax.Defined = false; minMax.MaxIndex.set_size(VImageDimension); minMax.MaxIndex.set_size(VImageDimension); typename ImageType::RegionType allowedExtremaRegion = inputImage->GetLargestPossibleRegion(); bool keepDistanceToImageBorders( neccessaryDistanceToImageBorderInMM > 0 ); if (keepDistanceToImageBorders) { long distanceInPixels[VImageDimension]; for(unsigned short dimension = 0; dimension < VImageDimension; ++dimension) { // To confirm that the whole hotspot is inside the image we have to keep a specific distance to the image-borders, which is as long as // the radius. To get the amount of indices we divide the radius by spacing and add 0.5 because voxels are center based: // For example with a radius of 2.2 and a spacing of 1 two indices are enough because 2.2 / 1 + 0.5 = 2.7 => 2. // But with a radius of 2.7 we need 3 indices because 2.7 / 1 + 0.5 = 3.2 => 3 distanceInPixels[dimension] = int( neccessaryDistanceToImageBorderInMM / spacing[dimension] + 0.5); } allowedExtremaRegion.ShrinkByRadius(distanceInPixels); } InputImageIndexIteratorType imageIndexIt(inputImage, allowedExtremaRegion); float maxValue = itk::NumericTraits::min(); float minValue = itk::NumericTraits::max(); typename ImageType::IndexType maxIndex; typename ImageType::IndexType minIndex; for(unsigned short i = 0; i < VImageDimension; ++i) { maxIndex[i] = 0; minIndex[i] = 0; } if (maskImage != nullptr) { MaskImageIteratorType maskIt(maskImage, maskImage->GetLargestPossibleRegion()); typename ImageType::IndexType imageIndex; typename ImageType::PointType worldPosition; typename ImageType::IndexType maskIndex; for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt) { imageIndex = maskIndex = maskIt.GetIndex(); if(maskIt.Get() == label) { if( allowedExtremaRegion.IsInside(imageIndex) ) { imageIndexIt.SetIndex( imageIndex ); double value = imageIndexIt.Get(); minMax.Defined = true; //Calculate minimum, maximum and corresponding index-values if( value > maxValue ) { maxIndex = imageIndexIt.GetIndex(); maxValue = value; } if(value < minValue ) { minIndex = imageIndexIt.GetIndex(); minValue = value; } } } } } else { for(imageIndexIt.GoToBegin(); !imageIndexIt.IsAtEnd(); ++imageIndexIt) { double value = imageIndexIt.Get(); minMax.Defined = true; //Calculate minimum, maximum and corresponding index-values if( value > maxValue ) { maxIndex = imageIndexIt.GetIndex(); maxValue = value; } if(value < minValue ) { minIndex = imageIndexIt.GetIndex(); minValue = value; } } } minMax.MaxIndex.set_size(VImageDimension); minMax.MinIndex.set_size(VImageDimension); for(unsigned int i = 0; i < minMax.MaxIndex.size(); ++i) { minMax.MaxIndex[i] = maxIndex[i]; } for(unsigned int i = 0; i < minMax.MinIndex.size(); ++i) { minMax.MinIndex[i] = minIndex[i]; } minMax.Max = maxValue; minMax.Min = minValue; return minMax; } template itk::Size - HotspotCalculator::CalculateConvolutionKernelSize( double spacing[VImageDimension], + HotspotMaskGenerator::CalculateConvolutionKernelSize( double spacing[VImageDimension], double radiusInMM ) { typedef itk::Image< float, VImageDimension > KernelImageType; typedef typename KernelImageType::SizeType SizeType; SizeType maskSize; for(unsigned int i = 0; i < VImageDimension; ++i) { maskSize[i] = static_cast( 2 * radiusInMM / spacing[i]); // We always want an uneven size to have a clear center point in the convolution mask if(maskSize[i] % 2 == 0 ) { ++maskSize[i]; } } return maskSize; } template itk::SmartPointer< itk::Image > - HotspotCalculator::GenerateHotspotSearchConvolutionKernel(double mmPerPixel[VImageDimension], + HotspotMaskGenerator::GenerateHotspotSearchConvolutionKernel(double mmPerPixel[VImageDimension], double radiusInMM ) { std::stringstream ss; for (unsigned int i = 0; i < VImageDimension; ++i) { ss << mmPerPixel[i]; if (i < VImageDimension -1) ss << ","; } MITK_DEBUG << "Update convolution kernel for spacing (" << ss.str() << ") and radius " << radiusInMM << "mm"; double radiusInMMSquared = radiusInMM * radiusInMM; typedef itk::Image< float, VImageDimension > KernelImageType; typename KernelImageType::Pointer convolutionKernel = KernelImageType::New(); // Calculate size and allocate mask image typedef typename KernelImageType::SizeType SizeType; SizeType maskSize = this->CalculateConvolutionKernelSize(mmPerPixel, radiusInMM); mitk::Point3D convolutionMaskCenterIndex; convolutionMaskCenterIndex.Fill(0.0); for(unsigned int i = 0; i < VImageDimension; ++i) { convolutionMaskCenterIndex[i] = 0.5 * (double)(maskSize[i]-1); } typedef typename KernelImageType::IndexType IndexType; IndexType maskIndex; maskIndex.Fill(0); typedef typename KernelImageType::RegionType RegionType; RegionType maskRegion; maskRegion.SetSize(maskSize); maskRegion.SetIndex(maskIndex); convolutionKernel->SetRegions(maskRegion); convolutionKernel->SetSpacing(mmPerPixel); convolutionKernel->Allocate(); // Fill mask image values by subsampling the image grid typedef itk::ImageRegionIteratorWithIndex MaskIteratorType; MaskIteratorType maskIt(convolutionKernel,maskRegion); int numberOfSubVoxelsPerDimension = 2; // per dimension! int numberOfSubVoxels = ::pow( static_cast(numberOfSubVoxelsPerDimension), static_cast(VImageDimension) ); double subVoxelSizeInPixels = 1.0 / (double)numberOfSubVoxelsPerDimension; double valueOfOneSubVoxel = 1.0 / (double)numberOfSubVoxels; double maskValue = 0.0; mitk::Point3D subVoxelIndexPosition; double distanceSquared = 0.0; typedef itk::ContinuousIndex ContinuousIndexType; for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt) { ContinuousIndexType indexPoint(maskIt.GetIndex()); mitk::Point3D voxelPosition; for (unsigned int dimension = 0; dimension < VImageDimension; ++dimension) { voxelPosition[dimension] = indexPoint[dimension]; } maskValue = 0.0; mitk::Vector3D subVoxelOffset; subVoxelOffset.Fill(0.0); // iterate sub-voxels by iterating all possible offsets for (subVoxelOffset[0] = -0.5 + subVoxelSizeInPixels / 2.0; subVoxelOffset[0] < +0.5; subVoxelOffset[0] += subVoxelSizeInPixels) { for (subVoxelOffset[1] = -0.5 + subVoxelSizeInPixels / 2.0; subVoxelOffset[1] < +0.5; subVoxelOffset[1] += subVoxelSizeInPixels) { for (subVoxelOffset[2] = -0.5 + subVoxelSizeInPixels / 2.0; subVoxelOffset[2] < +0.5; subVoxelOffset[2] += subVoxelSizeInPixels) { subVoxelIndexPosition = voxelPosition + subVoxelOffset; // this COULD be integrated into the for-loops if neccessary (add voxelPosition to initializer and end condition) distanceSquared = (subVoxelIndexPosition[0]-convolutionMaskCenterIndex[0]) * mmPerPixel[0] * (subVoxelIndexPosition[0]-convolutionMaskCenterIndex[0]) * mmPerPixel[0] + (subVoxelIndexPosition[1]-convolutionMaskCenterIndex[1]) * mmPerPixel[1] * (subVoxelIndexPosition[1]-convolutionMaskCenterIndex[1]) * mmPerPixel[1] + (subVoxelIndexPosition[2]-convolutionMaskCenterIndex[2]) * mmPerPixel[2] * (subVoxelIndexPosition[2]-convolutionMaskCenterIndex[2]) * mmPerPixel[2]; if (distanceSquared <= radiusInMMSquared) { maskValue += valueOfOneSubVoxel; } } } } maskIt.Set( maskValue ); } return convolutionKernel; } template itk::SmartPointer > - HotspotCalculator::GenerateConvolutionImage( const itk::Image* inputImage ) + HotspotMaskGenerator::GenerateConvolutionImage( const itk::Image* inputImage ) { double mmPerPixel[VImageDimension]; for (unsigned int dimension = 0; dimension < VImageDimension; ++dimension) { mmPerPixel[dimension] = inputImage->GetSpacing()[dimension]; } // update convolution kernel typedef itk::Image< float, VImageDimension > KernelImageType; typename KernelImageType::Pointer convolutionKernel = this->GenerateHotspotSearchConvolutionKernel(mmPerPixel, m_HotspotRadiusinMM); // update convolution image typedef itk::Image< TPixel, VImageDimension > InputImageType; typedef itk::Image< TPixel, VImageDimension > ConvolutionImageType; typedef itk::FFTConvolutionImageFilter ConvolutionFilterType; typename ConvolutionFilterType::Pointer convolutionFilter = ConvolutionFilterType::New(); typedef itk::ConstantBoundaryCondition BoundaryConditionType; BoundaryConditionType boundaryCondition; boundaryCondition.SetConstant(0.0); if (m_HotspotMustBeCompletelyInsideImage) { // overwrite default boundary condition convolutionFilter->SetBoundaryCondition(&boundaryCondition); } convolutionFilter->SetInput(inputImage); convolutionFilter->SetKernelImage(convolutionKernel); convolutionFilter->SetNormalize(true); MITK_DEBUG << "Update Convolution image for hotspot search"; convolutionFilter->UpdateLargestPossibleRegion(); typename ConvolutionImageType::Pointer convolutionImage = convolutionFilter->GetOutput(); convolutionImage->SetSpacing( inputImage->GetSpacing() ); // only workaround because convolution filter seems to ignore spacing of input image - m_modified = false; return convolutionImage; } template < typename TPixel, unsigned int VImageDimension> void - HotspotCalculator::FillHotspotMaskPixels( itk::Image* maskImage, + HotspotMaskGenerator::FillHotspotMaskPixels( itk::Image* maskImage, itk::Point sphereCenter, double sphereRadiusInMM ) { typedef itk::Image< TPixel, VImageDimension > MaskImageType; typedef itk::ImageRegionIteratorWithIndex MaskImageIteratorType; MaskImageIteratorType maskIt(maskImage, maskImage->GetLargestPossibleRegion()); typename MaskImageType::IndexType maskIndex; typename MaskImageType::PointType worldPosition; // this is not very smart. I would rather use a 0 initialized mask (not the case here -> blame CalculateHotspotMask) and find the region where I need to iterate over, then iterate only over the small region for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt) { maskIndex = maskIt.GetIndex(); maskImage->TransformIndexToPhysicalPoint(maskIndex, worldPosition); maskIt.Set( worldPosition.EuclideanDistanceTo(sphereCenter) <= sphereRadiusInMM ? 1 : 0 ); } } template void - HotspotCalculator::CalculateHotspotMask(itk::Image* inputImage, + HotspotMaskGenerator::CalculateHotspotMask(itk::Image* inputImage, typename itk::Image::Pointer maskImage, unsigned int label) { typedef itk::Image< TPixel, VImageDimension > InputImageType; typedef itk::Image< TPixel, VImageDimension > ConvolutionImageType; typedef itk::Image< float, VImageDimension > KernelImageType; typedef itk::Image< unsigned short, VImageDimension > MaskImageType; typedef itk::ImageDuplicator< InputImageType > DuplicatorType; typedef itk::CastImageFilter< InputImageType, MaskImageType > CastFilterType; typename ConvolutionImageType::Pointer convolutionImage = this->GenerateConvolutionImage(inputImage); if (convolutionImage.IsNull()) { MITK_ERROR << "Empty convolution image in CalculateHotspotStatistics(). We should never reach this state (logic error)."; throw std::logic_error("Empty convolution image in CalculateHotspotStatistics()"); } // if mask image is not defined, create an image of the same size as inputImage and fill it with 1's // there is maybe a better way to do this!? if (maskImage == nullptr) { maskImage = MaskImageType::New(); typename MaskImageType::RegionType maskRegion = inputImage->GetLargestPossibleRegion(); typename MaskImageType::SpacingType maskSpacing = inputImage->GetSpacing(); typename MaskImageType::PointType maskOrigin = inputImage->GetOrigin(); typename MaskImageType::DirectionType maskDirection = inputImage->GetDirection(); maskImage->SetRegions(maskRegion); maskImage->Allocate(); maskImage->SetOrigin(maskOrigin); maskImage->SetSpacing(maskSpacing); maskImage->SetDirection(maskDirection); maskImage->FillBuffer(1); label = 1; } // find maximum in convolution image, given the current mask double requiredDistanceToBorder = m_HotspotMustBeCompletelyInsideImage ? m_HotspotRadiusinMM : -1.0; ImageExtrema convolutionImageInformation = CalculateExtremaWorld(convolutionImage.GetPointer(), maskImage, requiredDistanceToBorder, label); bool isHotspotDefined = convolutionImageInformation.Defined; if (!isHotspotDefined) { MITK_ERROR << "No origin of hotspot-sphere was calculated!"; - m_hotspotMask = nullptr; + m_InternalMask = nullptr; } else { // create a binary mask around the "hotspot" region, fill the shape of a sphere around our hotspot center typename DuplicatorType::Pointer copyMachine = DuplicatorType::New(); copyMachine->SetInputImage(inputImage); copyMachine->Update(); typename CastFilterType::Pointer caster = CastFilterType::New(); caster->SetInput( copyMachine->GetOutput() ); caster->Update(); typename MaskImageType::Pointer hotspotMaskITK = caster->GetOutput(); typedef typename InputImageType::IndexType IndexType; IndexType maskCenterIndex; for (unsigned int d =0; d< VImageDimension;++d) { maskCenterIndex[d]=convolutionImageInformation.MaxIndex[d]; } typename ConvolutionImageType::PointType maskCenter; inputImage->TransformIndexToPhysicalPoint(maskCenterIndex,maskCenter); FillHotspotMaskPixels(hotspotMaskITK.GetPointer(), maskCenter, m_HotspotRadiusinMM); //obtain mitk::Image::Pointer from itk::Image mitk::Image::Pointer hotspotMaskAsMITKImage = mitk::GrabItkImageMemory(hotspotMaskITK); - m_hotspotMask = hotspotMaskAsMITKImage; + m_InternalMask = hotspotMaskAsMITKImage; } } } diff --git a/Modules/ImageStatistics/mitkHotspotCalculator.h b/Modules/ImageStatistics/mitkHotspotCalculator.h index 1dfd7514b6..0d7f6d9827 100644 --- a/Modules/ImageStatistics/mitkHotspotCalculator.h +++ b/Modules/ImageStatistics/mitkHotspotCalculator.h @@ -1,110 +1,124 @@ #ifndef MITKHOTSPOTCALCULATOR_H #define MITKHOTSPOTCALCULATOR_H #include #include #include #include #include #include #include +#include namespace mitk { - class MITKIMAGESTATISTICS_EXPORT HotspotCalculator: itk::Object + class MITKIMAGESTATISTICS_EXPORT HotspotMaskGenerator: public MaskGenerator { public: - HotspotCalculator(); + /** Standard Self typedef */ + typedef HotspotMaskGenerator Self; + typedef MaskGenerator Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; - ~HotspotCalculator(); + /** Method for creation through the object factory. */ + itkNewMacro(Self) - void setInputImage(mitk::Image::Pointer inputImage); + /** Runtime information support. */ + itkTypeMacro(HotspotMaskGenerator, MaskGenerator) - void setMaskImage(mitk::Image::Pointer maskImage); + void SetImage(mitk::Image::Pointer inputImage); - void setHotspotRadiusInMM(double radiusInMillimeter); + void SetMask(MaskGenerator::Pointer mask); - const double& getHotspotRadiusinMM() const; + void SetHotspotRadiusInMM(double radiusInMillimeter); - void setHotspotMustBeCompletelyInsideImage(bool hotspotCompletelyInsideImage); + const double& GetHotspotRadiusinMM() const; - bool getHotspotMustBeCompletelyInsideImage() const; + void SetHotspotMustBeCompletelyInsideImage(bool hotspotCompletelyInsideImage); - mitk::Image::Pointer getHotspotMask(unsigned int timeStep, unsigned int label=0); + bool GetHotspotMustBeCompletelyInsideImage() const; + + void SetLabel(unsigned int label); + + mitk::Image::Pointer GetMask(); protected: + HotspotMaskGenerator(); + + ~HotspotMaskGenerator(); + class ImageExtrema { public: bool Defined; double Max; double Min; vnl_vector MaxIndex; vnl_vector MinIndex; ImageExtrema() :Defined(false) ,Max(itk::NumericTraits::min()) ,Min(itk::NumericTraits::max()) { } }; private: /** \brief Returns size of convolution kernel depending on spacing and radius. */ template itk::Size CalculateConvolutionKernelSize(double spacing[VImageDimension], double radiusInMM); /** \brief Generates image of kernel which is needed for convolution. */ template itk::SmartPointer< itk::Image > GenerateHotspotSearchConvolutionKernel(double spacing[VImageDimension], double radiusInMM); /** \brief Convolves image with spherical kernel image. Used for hotspot calculation. */ template itk::SmartPointer< itk::Image > GenerateConvolutionImage( const itk::Image* inputImage ); /** \brief Fills pixels of the spherical hotspot mask. */ template < typename TPixel, unsigned int VImageDimension> void FillHotspotMaskPixels( itk::Image* maskImage, itk::Point sphereCenter, double sphereRadiusInMM); /** \brief */ template void CalculateHotspotMask(itk::Image* inputImage, typename itk::Image::Pointer maskImage, unsigned int label); template ImageExtrema CalculateExtremaWorld( const itk::Image* inputImage, typename itk::Image::Pointer maskImage, double neccessaryDistanceToImageBorderInMM, unsigned int label); - HotspotCalculator(const HotspotCalculator &); - HotspotCalculator & operator=(const HotspotCalculator &); + HotspotMaskGenerator(const HotspotMaskGenerator &); + HotspotMaskGenerator & operator=(const HotspotMaskGenerator &); mitk::Image::Pointer m_Image; - mitk::Image::Pointer m_Mask = nullptr; + MaskGenerator::Pointer m_Mask; mitk::Image::Pointer m_internalImage; itk::Image::Pointer m_internalMask2D; itk::Image::Pointer m_internalMask3D; - mitk::Image::Pointer m_hotspotMask; double m_HotspotRadiusinMM; bool m_HotspotMustBeCompletelyInsideImage; bool m_HotspotParamsChanged; - bool m_modified = false; + unsigned int m_Label; }; } #endif // MITKHOTSPOTCALCULATOR diff --git a/Modules/ImageStatistics/mitkImageStatisticsCalculator2.cpp b/Modules/ImageStatistics/mitkImageStatisticsCalculator2.cpp index 66bce82291..cc2f1af0c1 100644 --- a/Modules/ImageStatistics/mitkImageStatisticsCalculator2.cpp +++ b/Modules/ImageStatistics/mitkImageStatisticsCalculator2.cpp @@ -1,618 +1,659 @@ #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include +#include "itkImageFileWriter.h" namespace mitk { void ImageStatisticsCalculator_v2::SetInputImage(mitk::Image::Pointer image) { if (image != m_Image) { m_Image = image; m_HistogramStatisticsUpdateRequiredByTimeStep.resize(m_Image->GetTimeSteps()); m_StatisticsUpdateRequiredByTimeStep.resize(m_Image->GetTimeSteps()); m_StatisticsByTimeStep.resize(m_Image->GetTimeSteps()); m_HistogramStatisticsByTimeStep.resize(m_Image->GetTimeSteps()); m_HistogramsByTimeStep.resize(m_Image->GetTimeSteps()); this->SetAllStatisticsToUpdateRequired(); this->SetAllHistogramStatisticsToUpdateRequired(); } } void ImageStatisticsCalculator_v2::SetMask(mitk::MaskGenerator::Pointer mask) { if (mask != m_MaskGenerator) { m_MaskGenerator = mask; this->SetAllStatisticsToUpdateRequired(); this->SetAllHistogramStatisticsToUpdateRequired(); } } void ImageStatisticsCalculator_v2::SetDoHistogramStatistics(bool doHistStatistics) { m_DoHistogramStatistics = doHistStatistics; } bool ImageStatisticsCalculator_v2::GetDoHistogramStatistics() const { return m_DoHistogramStatistics; } void ImageStatisticsCalculator_v2::SetNBinsForHistogramStatistics(unsigned int nBins) { if (nBins != m_nBinsForHistogramStatistics) { m_nBinsForHistogramStatistics = nBins; this->SetAllHistogramStatisticsToUpdateRequired(); } } unsigned int ImageStatisticsCalculator_v2::GetNBinsForHistogramStatistics() const { return m_nBinsForHistogramStatistics; } ImageStatisticsCalculator_v2::statisticsMapType ImageStatisticsCalculator_v2::GetStatistics(unsigned int timeStep, unsigned int label) { + std::cout << "getting statistics..." << std::endl; if (timeStep >= m_StatisticsByTimeStep.size()) { throw std::runtime_error("invalid timeStep in ImageStatisticsCalculator_v2::GetStatistics"); } if (m_Image.IsNull()) { throw std::runtime_error("no image"); } if (m_MaskGenerator.IsNotNull()) { m_MaskGenerator->SetTimeStep(timeStep); m_InternalMask = m_MaskGenerator->GetMask(); } if (m_StatisticsUpdateRequiredByTimeStep[timeStep]) { + std::cout << "calculating statistics..." << std::endl; + // Calculate statistics with/without mask if (m_MaskGenerator.IsNull()) { // 1) calculate statistics unmasked: // plug image into itkstatisticsimagefilter (will be replaced by my awesome filter later on) // retrieve statistics and save them + std::cout << "unmasked" << std::endl; AccessByItk_1(m_Image, InternalCalculateStatisticsUnmasked, timeStep) } else { + std::cout << "masked" << std::endl; + // 2) calculate statistics masked // extract mask image region // plug mask and image into itklabelstatisticsimagefilter AccessByItk_1(m_Image, InternalCalculateStatisticsMasked, timeStep) } m_StatisticsUpdateRequiredByTimeStep[timeStep]=false; } if(m_DoHistogramStatistics && m_HistogramStatisticsUpdateRequiredByTimeStep[timeStep]) { + std::cout << "histogram statistics..." << std::endl; if (m_MaskGenerator.IsNull()) { + std::cout << "unmasked" << std::endl; + // calculate histogram unmasked AccessByItk_3(m_Image, InternalCalculateHistogramUnmasked, m_StatisticsByTimeStep[timeStep][0]->GetMin(), m_StatisticsByTimeStep[timeStep][0]->GetMax(), timeStep); m_HistogramStatisticsByTimeStep[timeStep].resize(1); HistogramStatisticsContainer::Pointer histogramStatistics = this->InternalCalculateHistogramStatistics(m_HistogramsByTimeStep[timeStep][0].histogram); histogramStatistics->GetLabel() = 0; m_HistogramStatisticsByTimeStep[timeStep][0] = histogramStatistics; } else { + std::cout << "masked" << std::endl; + m_HistogramStatisticsByTimeStep[timeStep].resize(0); // the label histograms have already been calculated in InternalCalculateStatisticsMasked for (HistogramContainer histContainer:m_HistogramsByTimeStep[timeStep]) { HistogramStatisticsContainer::Pointer histogramStatistics = this->InternalCalculateHistogramStatistics(histContainer.histogram); histogramStatistics->GetLabel() = histContainer.label; + m_HistogramStatisticsByTimeStep[timeStep].push_back(histogramStatistics); } } m_HistogramStatisticsUpdateRequiredByTimeStep[timeStep] = false; } // merge histogram and regular statistics, return result statisticsMapType statistics; if (m_MaskGenerator.IsNull()) { statistics = m_StatisticsByTimeStep[timeStep][0]->GetStatisticsAsMap(); if (m_DoHistogramStatistics) { statisticsMapType histStatistics = m_HistogramStatisticsByTimeStep[timeStep][0]->GetStatisticsAsMap(); statistics.insert(histStatistics.begin(), histStatistics.end()); } } else { bool statisticsFound(false); for (auto stats: m_StatisticsByTimeStep[timeStep]) { - if (stats->GetLabel() == label) { statistics = stats->GetStatisticsAsMap(); statisticsFound = true; break; } } if (!statisticsFound) { throw std::runtime_error("Invalid label"); } if (m_DoHistogramStatistics) { bool histStatsFound(false); for (auto histStats: m_HistogramStatisticsByTimeStep[timeStep]) { if (histStats->GetLabel() == label) { statisticsMapType histStatsMap = histStats->GetStatisticsAsMap(); statistics.insert(histStatsMap.begin(), histStatsMap.end()); histStatsFound = true; break; } } if (!histStatsFound) { throw std::runtime_error("Invalid label"); } } } return statistics; } void ImageStatisticsCalculator_v2::SetAllStatisticsToUpdateRequired() { for (unsigned int i = 0; i < m_StatisticsUpdateRequiredByTimeStep.size(); i++) { this->SetStatsTimeStepToUpdateRequired(i); } } void ImageStatisticsCalculator_v2::SetStatsTimeStepToUpdateRequired(unsigned int timeStep) { if (timeStep >= m_StatisticsUpdateRequiredByTimeStep.size()) { throw std::runtime_error("invalid timeStep in ImageStatisticsCalculator_v2::SetStatsTimeStepToUpdateRequired"); } m_StatisticsUpdateRequiredByTimeStep[timeStep] = true; m_StatisticsByTimeStep[timeStep].resize(0); } void ImageStatisticsCalculator_v2::SetAllHistogramStatisticsToUpdateRequired() { for (unsigned int i = 0; i < m_HistogramStatisticsUpdateRequiredByTimeStep.size(); i++) { this->SetHistStatsTimeStepToUpdateRequired(i); } } void ImageStatisticsCalculator_v2::SetHistStatsTimeStepToUpdateRequired(unsigned int timeStep) { if (timeStep >= m_HistogramStatisticsUpdateRequiredByTimeStep.size()) { throw std::runtime_error("invalid timeStep in ImageStatisticsCalculator_v2::SetHistStatsTimeStepToUpdateRequired"); } m_HistogramStatisticsUpdateRequiredByTimeStep[timeStep] = true; m_HistogramStatisticsByTimeStep[timeStep].resize(0); } template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator_v2::InternalCheckMaskSanity( typename itk::Image< TPixel, VImageDimension >* image, typename itk::Image< unsigned short, VImageDimension >::Pointer maskImage) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef typename ImageType::PointType PointType; typedef typename ImageType::DirectionType DirectionType; // check direction DirectionType imageDirection = image->GetDirection(); DirectionType maskDirection = maskImage->GetDirection(); for( int i = 0; i < imageDirection.ColumnDimensions; ++i ) { for( int j = 0; j < imageDirection.ColumnDimensions; ++j ) { double differenceDirection = imageDirection[i][j] - maskDirection[i][j]; if ( fabs( differenceDirection ) > mitk::eps ) { double differenceDirection = imageDirection[i][j] - maskDirection[i][j]; if ( fabs( differenceDirection ) > 0.001 /*mitk::eps*/ ) // TODO: temp fix (bug 17121) { itkExceptionMacro( << "Mask needs to have same direction as image! (Image direction: " << imageDirection << "; Mask direction: " << maskDirection << ")" ); } } } } // check spacing PointType imageSpacing = image->GetSpacing(); PointType maskSpacing = maskImage->GetSpacing(); for (unsigned int i = 0; i < VImageDimension; i++) { if ( fabs( maskSpacing[i] - imageSpacing[i] ) > mitk::eps ) { itkExceptionMacro(<< "Spacing of mask and image is not equal. Mask: " << maskSpacing << " image: " << imageSpacing ); } } // check alignment // Make sure that the voxels of mask and image are correctly "aligned", i.e., voxel boundaries are the same in both images PointType imageOrigin = image->GetOrigin(); PointType maskOrigin = maskImage->GetOrigin(); typedef itk::ContinuousIndex ContinousIndexType; ContinousIndexType maskOriginContinousIndex, imageOriginContinousIndex; image->TransformPhysicalPointToContinuousIndex(maskOrigin, maskOriginContinousIndex); image->TransformPhysicalPointToContinuousIndex(imageOrigin, imageOriginContinousIndex); for ( unsigned int i = 0; i < ImageType::ImageDimension; ++i ) { double misalignment = maskOriginContinousIndex[i] - floor( maskOriginContinousIndex[i] + 0.5 ); // misalignment must be a multiple (int) of spacing in that direction if ( fmod(misalignment,imageSpacing[i]) > mitk::eps ) { itkWarningMacro( << "Pixels/voxels of mask and image are not sufficiently aligned! (Misalignment: " << fmod(misalignment,imageSpacing[i]) << ")" ); } } // mask must be completely inside image region // Make sure that mask region is contained within image region if ( maskImage.IsNotNull() && !image->GetLargestPossibleRegion().IsInside( maskImage->GetLargestPossibleRegion() ) ) { itkWarningMacro( << "Mask region needs to be inside of image region! (Image region: " << image->GetLargestPossibleRegion() << "; Mask region: " << maskImage->GetLargestPossibleRegion() << ")" ); } } template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator_v2::InternalCalculateStatisticsUnmasked( typename itk::Image< TPixel, VImageDimension >* image, unsigned int timeStep) { typedef typename itk::Image< TPixel, VImageDimension > ImageType; typename itk::ExtendedStatisticsImageFilter::Pointer statisticsFilter = itk::ExtendedStatisticsImageFilter::New(); statisticsFilter->SetInput(image); statisticsFilter->SetCoordinateTolerance(0.001); statisticsFilter->SetDirectionTolerance(0.001); // TODO: REMOVE THIS, the histogram statistics should be a separate module statisticsFilter->SetBinSize(m_nBinsForHistogramStatistics); try { statisticsFilter->Update(); } catch (const itk::ExceptionObject& e) { mitkThrow() << "Image statistics calculation failed due to following ITK Exception: \n " << e.what(); } // no mask, therefore just one label = the whole image m_StatisticsByTimeStep[timeStep].resize(1); StatisticsContainer::Pointer statisticsResult = StatisticsContainer::New(); statisticsResult->GetLabel() = 0; statisticsResult->GetN() = image->GetLargestPossibleRegion().GetNumberOfPixels(); statisticsResult->GetMean() = statisticsFilter->GetMean(); statisticsResult->GetMin() = statisticsFilter->GetMinimum(); statisticsResult->GetMax() = statisticsFilter->GetMaximum(); statisticsResult->GetVariance() = statisticsFilter->GetVariance(); statisticsResult->GetStd() = statisticsFilter->GetSigma(); statisticsResult->GetSkewness() = statisticsFilter->GetSkewness(); statisticsResult->GetKurtosis() = statisticsFilter->GetKurtosis(); - //statisticsResult->GetRMS() = statisticsFilter->Get(); + statisticsResult->GetRMS() = std::sqrt(std::pow(statisticsFilter->GetMean(), 2.) + statisticsFilter->GetVariance()); // variance = sigma^2 statisticsResult->GetMPP() = statisticsFilter->GetMPP(); m_StatisticsByTimeStep[timeStep][0] = statisticsResult; } template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator_v2::InternalCalculateHistogramUnmasked( typename itk::Image< TPixel, VImageDimension >* image, double minVal, double maxVal, unsigned int timeStep) { typedef typename itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Statistics::ImageToHistogramFilter HistogramFilterType; typedef itk::Statistics::Histogram HistogramType; typedef typename HistogramFilterType::HistogramSizeType SizeType; typename HistogramFilterType::Pointer histogramFilter = HistogramFilterType::New(); histogramFilter->SetInput(image); SizeType histSize(1); histSize[0] = m_nBinsForHistogramStatistics; histogramFilter->SetHistogramSize(histSize); // is that correct?? typename HistogramFilterType::HistogramMeasurementVectorType lowerBound(1); typename HistogramFilterType::HistogramMeasurementVectorType upperBound(1); lowerBound[0] = minVal; upperBound[0] = maxVal; histogramFilter->SetHistogramBinMinimum(lowerBound); histogramFilter->SetHistogramBinMaximum(upperBound); histogramFilter->Update(); typename HistogramType::Pointer histogram = histogramFilter->GetOutput(); HistogramContainer histogramContainer; histogramContainer.label=0; histogramContainer.histogram = histogram; m_HistogramsByTimeStep[timeStep].resize(1); m_HistogramsByTimeStep[timeStep][0] = histogramContainer; } template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator_v2::ExtractMaskImageRegion( typename itk::Image< TPixel, VImageDimension >* image, typename itk::Image< unsigned short, VImageDimension >::Pointer maskImage, - typename itk::Image< TPixel, VImageDimension >::Pointer adaptedImage + typename itk::Image< TPixel, VImageDimension >::Pointer& adaptedImage ) { typedef itk::Image< TPixel, VImageDimension > ImageType; + typedef itk::Image< unsigned short, VImageDimension > MaskImageType; typedef itk::ExtractImageFilter< ImageType, ImageType > ExtractImageFilterType; + typedef itk::ChangeInformationImageFilter< ImageType > ChangeInformationFilterType; typename ImageType::SizeType imageSize = image->GetBufferedRegion().GetSize(); typename ImageType::SizeType maskSize = maskImage->GetBufferedRegion().GetSize(); + + bool maskSmallerImage = false; for ( unsigned int i = 0; i < ImageType::ImageDimension; ++i ) { if ( maskSize[i] < imageSize[i] ) { maskSmallerImage = true; } } if ( maskSmallerImage ) { typename ExtractImageFilterType::Pointer extractImageFilter = ExtractImageFilterType::New(); + typename MaskImageType::PointType maskOrigin = maskImage->GetOrigin(); + typename ImageType::PointType imageOrigin = image->GetOrigin(); + typename MaskImageType::SpacingType maskSpacing = maskImage->GetSpacing(); + typename ImageType::RegionType extractionRegion; + typename ImageType::IndexType extractionRegionIndex; + + + for (unsigned int i=0; i < maskOrigin.GetPointDimension(); i++) + { + extractionRegionIndex[i] = (maskOrigin[i] - imageOrigin[i]) / maskSpacing[i]; + } + + extractionRegion.SetIndex(extractionRegionIndex); + extractionRegion.SetSize(maskImage->GetLargestPossibleRegion().GetSize()); + extractImageFilter->SetInput( image ); - extractImageFilter->SetExtractionRegion( maskImage->GetBufferedRegion() ); + extractImageFilter->SetExtractionRegion( extractionRegion ); extractImageFilter->SetCoordinateTolerance( 0.001 ); extractImageFilter->SetDirectionTolerance( 0.001 ); extractImageFilter->Update(); adaptedImage = extractImageFilter->GetOutput(); + adaptedImage->SetOrigin(maskImage->GetOrigin()); + adaptedImage->SetLargestPossibleRegion(maskImage->GetLargestPossibleRegion()); + adaptedImage->SetBufferedRegion(maskImage->GetBufferedRegion()); + std::cout << "sdgversg" << std::endl; + } else { adaptedImage = image; } } /* Taken from original ImageStatisticsCalculator * This function is not ideal on so many levels 1) it does not support more than one label 2) it assumes that the label id is 1 (which does not have to be the case) 3) single threaded => I suggest writing an itk filter that finds min and max value. This filter could then be called from our new statisticsimagefilter in the beforethreadedgenereatedata (provided that a min/max for the histogram is not yet defined) EDIT: I changed this function so that it simply searches the image min and max and completely disregards the mask. Since we call this function after we cropped the image to the mask region we can run it on the LargestPossibleRegion of the image */ template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator_v2::GetMinAndMaxValue( typename itk::Image< TPixel, VImageDimension >* inputImage, double &minVal, double &maxVal ) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::ImageRegionConstIteratorWithIndex ImageRegionIteratorType; ImageRegionIteratorType imageIterator( inputImage, inputImage->GetLargestPossibleRegion() ); maxVal = std::numeric_limits::min(); minVal = std::numeric_limits::max(); long counter(0); double actualPixelValue(0); for( imageIterator.GoToBegin(); !imageIterator.IsAtEnd(); ++imageIterator) { - if( imageIterator.Value()== 1.0) - { counter++; actualPixelValue = imageIterator.Value(); if(actualPixelValue > maxVal) { maxVal = actualPixelValue; } else if(actualPixelValue < minVal) { minVal = actualPixelValue; } - - } } } template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator_v2::InternalCalculateStatisticsMasked( typename itk::Image< TPixel, VImageDimension >* image, unsigned int timeStep) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< unsigned short, VImageDimension > MaskType; typedef itk::ExtendedLabelStatisticsImageFilter< ImageType, MaskType > ImageStatisticsFilterType; // maskImage has to have the same dimension as image typename MaskType::Pointer maskImage = MaskType::New(); maskImage = ImageToItkImage< unsigned short, VImageDimension >(m_InternalMask); InternalCheckMaskSanity(image, maskImage); // if mask is smaller than image, extract the image region where the mask is typename ImageType::Pointer adaptedImage = ImageType::New(); // uses m_InternalMask this->ExtractMaskImageRegion(image, maskImage, adaptedImage); double minVal, maxVal; GetMinAndMaxValue(adaptedImage.GetPointer(), minVal, maxVal); + std::string outfile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/adapted_image.nrrd"; + typedef itk::ImageFileWriter< ImageType > WriterType; + typename WriterType::Pointer writer = WriterType::New(); + writer->SetFileName(outfile); + writer->SetInput(adaptedImage); + writer->Update(); + typename ImageStatisticsFilterType::Pointer imageStatisticsFilter = ImageStatisticsFilterType::New(); imageStatisticsFilter->SetDirectionTolerance(0.001); imageStatisticsFilter->SetCoordinateTolerance(0.001); - imageStatisticsFilter->SetHistogramParameters(m_nBinsForHistogramStatistics, minVal, maxVal); imageStatisticsFilter->SetInput(adaptedImage); imageStatisticsFilter->SetLabelInput(maskImage); + imageStatisticsFilter->SetHistogramParameters(m_nBinsForHistogramStatistics, floor(minVal), ceil(maxVal)); imageStatisticsFilter->Update(); std::list labels = imageStatisticsFilter->GetRelevantLabels(); std::list::iterator it = labels.begin(); m_StatisticsByTimeStep[timeStep].resize(0); while(it != labels.end()) { StatisticsContainer::Pointer statisticsResult = StatisticsContainer::New(); statisticsResult->GetN() = imageStatisticsFilter->GetSum(*it) / (double) imageStatisticsFilter->GetMean(*it); statisticsResult->GetMean() = imageStatisticsFilter->GetMean(*it); statisticsResult->GetMin() = imageStatisticsFilter->GetMinimum(*it); statisticsResult->GetMax() = imageStatisticsFilter->GetMaximum(*it); statisticsResult->GetVariance() = imageStatisticsFilter->GetVariance(*it); statisticsResult->GetStd() = imageStatisticsFilter->GetSigma(1); statisticsResult->GetSkewness() = imageStatisticsFilter->GetSkewness(*it); statisticsResult->GetKurtosis() = imageStatisticsFilter->GetKurtosis(*it); - //statisticsResult->GetRMS() = imageStatisticsFilter->GetRMS(*it); + statisticsResult->GetRMS() = std::sqrt(std::pow(imageStatisticsFilter->GetMean(*it), 2.) + imageStatisticsFilter->GetVariance(*it)); // variance = sigma^2 statisticsResult->GetMPP() = imageStatisticsFilter->GetMPP(*it); statisticsResult->GetLabel() = *it; + std::cout << imageStatisticsFilter->GetUPP(*it) << " " << imageStatisticsFilter->GetUniformity(*it) << " " << imageStatisticsFilter->GetEntropy(*it) << std::endl; m_StatisticsByTimeStep[timeStep].push_back(statisticsResult); if (m_DoHistogramStatistics) { HistogramContainer histCont; histCont.label = *it; histCont.histogram = imageStatisticsFilter->GetHistogram(*it); m_HistogramsByTimeStep[timeStep].push_back(histCont); } + ++it; } } ImageStatisticsCalculator_v2::HistogramStatisticsContainer::Pointer ImageStatisticsCalculator_v2::InternalCalculateHistogramStatistics(itk::Statistics::Histogram::Pointer histogram) { HistogramStatisticsCalculator histStatCalc; histStatCalc.SetHistogram(histogram); histStatCalc.CalculateStatistics(); HistogramStatisticsContainer::Pointer histStatResults = HistogramStatisticsContainer::New(); histStatResults->GetEntropy() = histStatCalc.GetEntropy(); histStatResults->GetMedian() = histStatCalc.GetMedian(); histStatResults->GetUniformity() = histStatCalc.GetUniformity(); histStatResults->GetUPP() = histStatCalc.GetUPP(); return histStatResults; } ImageStatisticsCalculator_v2::StatisticsContainer::StatisticsContainer(): m_N(std::numeric_limits::max()), m_Mean(std::numeric_limits::max()), m_Min(std::numeric_limits::max()), m_Max(std::numeric_limits::max()), m_Std(std::numeric_limits::max()), m_Variance(std::numeric_limits::max()), m_Skewness(std::numeric_limits::max()), m_Kurtosis(std::numeric_limits::max()), m_RMS(std::numeric_limits::max()), m_MPP(std::numeric_limits::max()) { } ImageStatisticsCalculator_v2::statisticsMapType ImageStatisticsCalculator_v2::StatisticsContainer::GetStatisticsAsMap() { ImageStatisticsCalculator_v2::statisticsMapType statisticsAsMap; statisticsAsMap["N"] = m_N; statisticsAsMap["Mean"] = m_Mean; statisticsAsMap["Min"] = m_Min; statisticsAsMap["Max"] = m_Max; statisticsAsMap["StandardDeviation"] = m_Std; statisticsAsMap["Variance"] = m_Variance; statisticsAsMap["Skewness"] = m_Skewness; statisticsAsMap["Kurtosis"] = m_Kurtosis; statisticsAsMap["RMS"] = m_RMS; statisticsAsMap["MPP"] = m_MPP; return statisticsAsMap; } void ImageStatisticsCalculator_v2::StatisticsContainer::Reset() { m_N = std::numeric_limits::max(); m_Mean = std::numeric_limits::max(); m_Min = std::numeric_limits::max(); m_Max = std::numeric_limits::max(); m_Std = std::numeric_limits::max(); m_Variance = std::numeric_limits::max(); m_Skewness = std::numeric_limits::max(); m_Kurtosis = std::numeric_limits::max(); m_RMS = std::numeric_limits::max(); m_MPP = std::numeric_limits::max(); } ImageStatisticsCalculator_v2::HistogramStatisticsContainer::HistogramStatisticsContainer(): m_Median(std::numeric_limits::max()), m_Uniformity(std::numeric_limits::max()), m_UPP(std::numeric_limits::max()), m_Entropy(std::numeric_limits::max()) { } ImageStatisticsCalculator_v2::statisticsMapType ImageStatisticsCalculator_v2::HistogramStatisticsContainer::GetStatisticsAsMap() { ImageStatisticsCalculator_v2::statisticsMapType statisticsAsMap; statisticsAsMap["Median"] = m_Median; statisticsAsMap["Uniformity"] = m_Uniformity; statisticsAsMap["UPP"] = m_UPP; statisticsAsMap["Entropy"] = m_Entropy; return statisticsAsMap; } void ImageStatisticsCalculator_v2::HistogramStatisticsContainer::Reset() { m_Median = std::numeric_limits::max(); m_Uniformity = std::numeric_limits::max(); m_UPP = std::numeric_limits::max(); m_Entropy = std::numeric_limits::max(); } } diff --git a/Modules/ImageStatistics/mitkImageStatisticsCalculator2.h b/Modules/ImageStatistics/mitkImageStatisticsCalculator2.h index 17522a9d17..64452323bb 100644 --- a/Modules/ImageStatistics/mitkImageStatisticsCalculator2.h +++ b/Modules/ImageStatistics/mitkImageStatisticsCalculator2.h @@ -1,367 +1,367 @@ #ifndef MITKIMAGESTATISTICSCALCULATOR2 #define MITKIMAGESTATISTICSCALCULATOR2 #include #include #include #include #include #include #include /*#define mitkGetRefConstMacro(name, type) \ const type& Get##name() const \ { \ return &m_##name; \ } \ \ type& Get##name() \ { \ return &m_##name; \ } \ */ namespace mitk { class MITKIMAGESTATISTICS_EXPORT ImageStatisticsCalculator_v2: public itk::Object { public: /** Standard Self typedef */ typedef ImageStatisticsCalculator_v2 Self; typedef itk::Object Superclass; typedef itk::SmartPointer< Self > Pointer; typedef itk::SmartPointer< const Self > ConstPointer; /** Method for creation through the object factory. */ itkNewMacro(Self) /** Runtime information support. */ itkTypeMacro(ImageStatisticsCalculator_v2, itk::Object) typedef double statisticsValueType; typedef std::map statisticsMapType; typedef itk::Statistics::Histogram HistogramType; void SetInputImage(mitk::Image::Pointer image); void SetMask(mitk::MaskGenerator::Pointer mask); void SetDoHistogramStatistics(bool doHistStatistics); bool GetDoHistogramStatistics() const; void SetNBinsForHistogramStatistics(unsigned int nBins); unsigned int GetNBinsForHistogramStatistics() const; - statisticsMapType GetStatistics(unsigned int timeStep=0, unsigned int label=0); // not yet implemented + statisticsMapType GetStatistics(unsigned int timeStep=0, unsigned int label=1); class StatisticsContainer: public itk::Object { public: /** Standard Self typedef */ typedef StatisticsContainer Self; typedef itk::Object Superclass; typedef itk::SmartPointer< Self > Pointer; typedef itk::SmartPointer< const Self > ConstPointer; /** Method for creation through the object factory. */ itkNewMacro(Self) /** Runtime information support. */ itkTypeMacro(StatisticsContainer, itk::Object) typedef double RealType; statisticsMapType GetStatisticsAsMap(); void Reset(); long & GetN() { return m_N; } const long& GetN() const { return m_N; } RealType & GetMean() { return m_Mean; } const RealType& GetMean() const { return m_Mean; } RealType & GetVariance() { return m_Variance; } const RealType& GetVariance() const { return m_Variance; } RealType & GetStd() { return m_Std; } const RealType& GetStd() const { return m_Std; } RealType & GetMin() { return m_Min; } const RealType& GetMin() const { return m_Min; } RealType & GetMax() { return m_Max; } const RealType& GetMax() const { return m_Max; } RealType & GetRMS() { return m_RMS; } const RealType& GetRMS() const { return m_RMS; } RealType & GetSkewness() { return m_Skewness; } const RealType& GetSkewness() const { return m_Skewness; } RealType & GetKurtosis() { return m_Kurtosis; } const RealType& GetKurtosis() const { return m_Kurtosis; } RealType & GetMPP() { return m_MPP; } const RealType& GetMPP() const { return m_MPP; } unsigned int & GetLabel() { return m_Label; } const unsigned int& GetLabel() const { return m_Label; } protected: StatisticsContainer(); private: // not pretty, is temporary long m_N; RealType m_Mean, m_Min, m_Max, m_Std, m_Variance; RealType m_Skewness; RealType m_Kurtosis; RealType m_RMS; RealType m_MPP; unsigned int m_Label; }; class HistogramStatisticsContainer: public itk::Object { public: /** Standard Self typedef */ typedef HistogramStatisticsContainer Self; typedef itk::Object Superclass; typedef itk::SmartPointer< Self > Pointer; typedef itk::SmartPointer< const Self > ConstPointer; /** Method for creation through the object factory. */ itkNewMacro(Self) /** Runtime information support. */ itkTypeMacro(HistogramStatisticsContainer, itk::Object) typedef double RealType; statisticsMapType GetStatisticsAsMap(); void Reset(); RealType & GetMedian() { return m_Median; } const RealType& GetMedian() const { return m_Median; } RealType & GetUniformity() { return m_Uniformity; } const RealType& GetUniformity() const { return m_Uniformity; } RealType & GetEntropy() { return m_Entropy; } const RealType& GetEntropy() const { return m_Entropy; } RealType & GetUPP() { return m_UPP; } const RealType& GetUPP() const { return m_UPP; } unsigned int & GetLabel() { return m_Label; } const unsigned int& GetLabel() const { return m_Label; } protected: HistogramStatisticsContainer(); private: RealType m_Median; RealType m_Uniformity; RealType m_UPP; RealType m_Entropy; unsigned int m_Label; }; struct HistogramContainer { unsigned int label; itk::Statistics::Histogram::Pointer histogram; }; protected: ImageStatisticsCalculator_v2(){ m_nBinsForHistogramStatistics = 100; }; private: void SetAllStatisticsToUpdateRequired(); void SetAllHistogramStatisticsToUpdateRequired(); void SetStatsTimeStepToUpdateRequired(unsigned int timeStep); void SetHistStatsTimeStepToUpdateRequired(unsigned int timeStep); void CalculateHistogramStatistics(itk::Statistics::Histogram histogram, HistogramStatisticsContainer* HistStatContainer); template < typename TPixel, unsigned int VImageDimension > void InternalCheckMaskSanity( typename itk::Image< TPixel, VImageDimension >* image, typename itk::Image< unsigned short, VImageDimension >::Pointer maskImage); template < typename TPixel, unsigned int VImageDimension > void InternalCalculateStatisticsUnmasked( typename itk::Image< TPixel, VImageDimension >* image, unsigned int timeStep); HistogramStatisticsContainer::Pointer InternalCalculateHistogramStatistics(itk::Statistics::Histogram::Pointer histogram); template < typename TPixel, unsigned int VImageDimension > void InternalCalculateHistogramUnmasked( typename itk::Image< TPixel, VImageDimension >* image, double minVal, double maxVal, unsigned int timeStep); template < typename TPixel, unsigned int VImageDimension > void InternalCalculateStatisticsMasked( typename itk::Image< TPixel, VImageDimension >* image, unsigned int timeStep); template < typename TPixel, unsigned int VImageDimension > void ExtractMaskImageRegion( typename itk::Image< TPixel, VImageDimension >* image, typename itk::Image< unsigned short, VImageDimension >::Pointer maskImage, - typename itk::Image< TPixel, VImageDimension >::Pointer adaptedImage + typename itk::Image< TPixel, VImageDimension >::Pointer& adaptedImage ); template < typename TPixel, unsigned int VImageDimension > void GetMinAndMaxValue( typename itk::Image< TPixel, VImageDimension >* inputImage, double &minVal, double &maxVal ); std::string GetNameOfClass() { return std::string("ImageStatisticsCalculator_v2"); } mitk::Image::Pointer m_Image; mitk::Image::Pointer m_ImageTimeSlice; mitk::MaskGenerator::Pointer m_MaskGenerator; mitk::Image::Pointer m_InternalMask; bool m_DoHistogramStatistics; unsigned int m_nBinsForHistogramStatistics; std::vector m_StatisticsUpdateRequiredByTimeStep; // holds which time steps are valid and which ones have to be recalculated std::vector> m_StatisticsByTimeStep; std::vector> m_HistogramStatisticsByTimeStep; std::vector m_HistogramStatisticsUpdateRequiredByTimeStep; std::vector> m_HistogramsByTimeStep; }; } #endif // MITKIMAGESTATISTICSCALCULATOR2 diff --git a/Modules/ImageStatistics/mitkMaskGenerator.h b/Modules/ImageStatistics/mitkMaskGenerator.h index 4fbf28eff1..0a863852fd 100644 --- a/Modules/ImageStatistics/mitkMaskGenerator.h +++ b/Modules/ImageStatistics/mitkMaskGenerator.h @@ -1,48 +1,48 @@ #ifndef MITKMASKGENERATOR #define MITKMASKGENERATOR #include #include #include #include #include namespace mitk { class MITKIMAGESTATISTICS_EXPORT MaskGenerator: public itk::Object { public: /** Standard Self typedef */ typedef MaskGenerator Self; typedef itk::Object Superclass; typedef itk::SmartPointer< Self > Pointer; typedef itk::SmartPointer< const Self > ConstPointer; /** Method for creation through the object factory. */ itkNewMacro(Self) /** Runtime information support. */ itkTypeMacro(MaskGenerator, itk::Object) //~MaskGenerator(); - mitk::Image::Pointer GetMask(); + virtual mitk::Image::Pointer GetMask(); void SetTimeStep(unsigned int timeStep); protected: MaskGenerator(); void CheckMaskSanity(mitk::Image::Pointer image); unsigned int m_TimeStep; mitk::Image::Pointer m_InternalMask; bool m_Modified; private: }; } #endif // MITKMASKGENERATOR diff --git a/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.h b/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.h index 38c42f8cac..d5eb54f483 100644 --- a/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.h +++ b/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.h @@ -1,95 +1,106 @@ #ifndef MITKPLANARFIGUREMASKGENERATOR #define MITKPLANARFIGUREMASKGENERATOR #include #include #include #include #include #include #include #include #include #include namespace mitk { class MITKIMAGESTATISTICS_EXPORT PlanarFigureMaskGenerator: public MaskGenerator { public: + /** Standard Self typedef */ + typedef PlanarFigureMaskGenerator Self; + typedef MaskGenerator Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self) + + /** Runtime information support. */ + itkTypeMacro(PlanarFigureMaskGenerator, MaskGenerator) mitk::Image::Pointer GetMask(); void SetPlanarFigure(mitk::PlanarFigure::Pointer planarFigure); /*We need an image because we need its geometry in order to give the points of the planar figure the correct world coordinates*/ void SetImage(const mitk::Image::Pointer inputImage); protected: - + PlanarFigureMaskGenerator(){} private: void CalculateMask(); template < typename TPixel, unsigned int VImageDimension > void InternalCalculateMaskFromPlanarFigure( const itk::Image< TPixel, VImageDimension > *image, unsigned int axis ); mitk::Image::Pointer extract2DImageSlice(unsigned int axis, unsigned int slice); bool GetPrincipalAxis(const BaseGeometry *geometry, Vector3D vector, unsigned int &axis ); /** Connection from ITK to VTK */ template void ConnectPipelines(ITK_Exporter exporter, vtkSmartPointer importer) { importer->SetUpdateInformationCallback(exporter->GetUpdateInformationCallback()); importer->SetPipelineModifiedCallback(exporter->GetPipelineModifiedCallback()); importer->SetWholeExtentCallback(exporter->GetWholeExtentCallback()); importer->SetSpacingCallback(exporter->GetSpacingCallback()); importer->SetOriginCallback(exporter->GetOriginCallback()); importer->SetScalarTypeCallback(exporter->GetScalarTypeCallback()); importer->SetNumberOfComponentsCallback(exporter->GetNumberOfComponentsCallback()); importer->SetPropagateUpdateExtentCallback(exporter->GetPropagateUpdateExtentCallback()); importer->SetUpdateDataCallback(exporter->GetUpdateDataCallback()); importer->SetDataExtentCallback(exporter->GetDataExtentCallback()); importer->SetBufferPointerCallback(exporter->GetBufferPointerCallback()); importer->SetCallbackUserData(exporter->GetCallbackUserData()); } /** Connection from VTK to ITK */ template void ConnectPipelines(vtkSmartPointer exporter, ITK_Importer importer) { importer->SetUpdateInformationCallback(exporter->GetUpdateInformationCallback()); importer->SetPipelineModifiedCallback(exporter->GetPipelineModifiedCallback()); importer->SetWholeExtentCallback(exporter->GetWholeExtentCallback()); importer->SetSpacingCallback(exporter->GetSpacingCallback()); importer->SetOriginCallback(exporter->GetOriginCallback()); importer->SetScalarTypeCallback(exporter->GetScalarTypeCallback()); importer->SetNumberOfComponentsCallback(exporter->GetNumberOfComponentsCallback()); importer->SetPropagateUpdateExtentCallback(exporter->GetPropagateUpdateExtentCallback()); importer->SetUpdateDataCallback(exporter->GetUpdateDataCallback()); importer->SetDataExtentCallback(exporter->GetDataExtentCallback()); importer->SetBufferPointerCallback(exporter->GetBufferPointerCallback()); importer->SetCallbackUserData(exporter->GetCallbackUserData()); } mitk::PlanarFigure::Pointer m_PlanarFigure; typename itk::Image::Pointer m_InternalITKImageMask2D; mitk::Image::Pointer m_InternalInputImage; unsigned int m_PlanarFigureAxis; }; } #endif // MITKPLANARFIGUREMASKGENERATOR