diff --git a/Modules/DiffusionImaging/MiniApps/HotspotMiniApp.cpp b/Modules/DiffusionImaging/MiniApps/HotspotMiniApp.cpp index c9b3a74f69..903dff8eae 100644 --- a/Modules/DiffusionImaging/MiniApps/HotspotMiniApp.cpp +++ b/Modules/DiffusionImaging/MiniApps/HotspotMiniApp.cpp @@ -1,103 +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 +#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::HotspotMaskGenerator::Pointer hotspotCalc = mitk::HotspotMaskGenerator::New(); hotspotCalc->SetImage(inputImage); hotspotCalc->SetHotspotRadiusInMM(10); hotspotCalc->SetTimeStep(0); mitk::Image::Pointer outImage; try { 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 5fa25988cb..e4a44899cd 100644 --- a/Modules/DiffusionImaging/MiniApps/ImageStatisticsMiniapp_v2.cpp +++ b/Modules/DiffusionImaging/MiniApps/ImageStatisticsMiniapp_v2.cpp @@ -1,199 +1,200 @@ /*=================================================================== 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 +#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()); 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[] ) { unsigned int timeStep = 0; std::string inputImageFile; 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; 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(); + mitk::ImageMaskGenerator::Pointer binaryImageMaskGen = mitk::ImageMaskGenerator::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); - + calculator->SetNBinsForHistogramStatistics(100); try { - result = calculator->GetStatistics(timeStep, 1); +// result = calculator->GetStatistics(timeStep, 1); } catch( const itk::ExceptionObject& e) { MITK_ERROR << "Statistic Calculation Failed - ITK Exception:" << e.what(); return -1; } - printMap(result); +// printMap(result); - calculator->SetMask(binaryImageMaskGen.GetPointer()); - result = calculator->GetStatistics(timeStep, 1); - printMap(result); +// calculator->SetMask(binaryImageMaskGen.GetPointer()); +// result = calculator->GetStatistics(timeStep, 1); +// printMap(result); calculator->SetMask(nullptr); result = calculator->GetStatistics(timeStep); printMap(result); - calculator->SetMask(hotSpotMaskGen.GetPointer()); - result = calculator->GetStatistics(timeStep, 1); - printMap(result); +// calculator->SetMask(hotSpotMaskGen.GetPointer()); +// result = calculator->GetStatistics(timeStep, 1); +// printMap(result); // 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; + 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->SetMaskingModeToNone(); + 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/ImageStatistics/files.cmake b/Modules/ImageStatistics/files.cmake index 80a9d5d335..5c078b083a 100644 --- a/Modules/ImageStatistics/files.cmake +++ b/Modules/ImageStatistics/files.cmake @@ -1,28 +1,28 @@ set(CPP_FILES mitkImageStatisticsCalculator.cpp mitkImageStatisticsCalculator2.cpp mitkPointSetStatisticsCalculator.cpp mitkPointSetDifferenceStatisticsCalculator.cpp mitkIntensityProfile.cpp - mitkHotspotCalculator.cpp + mitkHotspotMaskGenerator.cpp mitkMaskGenerator.cpp mitkPlanarFigureMaskGenerator.cpp mitkMultiLabelMaskGenerator.cpp - mitkBinaryImageMaskGenerator.cpp + mitkImageMaskGenerator.cpp mitkHistogramStatisticsCalculator.cpp ) set(H_FILES mitkImageStatisticsCalculator.h mitkImageStatisticsCalculator2.h mitkPointSetDifferenceStatisticsCalculator.h mitkPointSetStatisticsCalculator.h mitkExtendedStatisticsImageFilter.h mitkExtendedLabelStatisticsImageFilter.h - mitkHotspotCalculator.h + mitkHotspotMaskGenerator.h mitkMaskGenerator.h mitkPlanarFigureMaskGenerator.h mitkMultiLabelMaskGenerator.h - mitkBinaryImageMaskGenerator.h + mitkImageMaskGenerator.h mitkHistogramStatisticsCalculator.h ) diff --git a/Modules/ImageStatistics/mitkHistogramStatisticsCalculator.cpp b/Modules/ImageStatistics/mitkHistogramStatisticsCalculator.cpp index 59e4f9febd..b7af6a479d 100644 --- a/Modules/ImageStatistics/mitkHistogramStatisticsCalculator.cpp +++ b/Modules/ImageStatistics/mitkHistogramStatisticsCalculator.cpp @@ -1,108 +1,108 @@ #include #include namespace mitk { HistogramStatisticsCalculator::HistogramStatisticsCalculator(): m_StatisticsCalculated(false) { } void HistogramStatisticsCalculator::SetHistogram(HistogramType::Pointer histogram) { if (m_Histogram != histogram) { m_Histogram = histogram; m_StatisticsCalculated = false; } } HistogramStatisticsCalculator::MeasurementType HistogramStatisticsCalculator::GetEntropy() { if (!m_StatisticsCalculated) { MITK_WARN("Statistics have not yet been calculated, running calculation now..."); CalculateStatistics(); } return m_Entropy; } HistogramStatisticsCalculator::MeasurementType HistogramStatisticsCalculator::GetMedian() { if (!m_StatisticsCalculated) { MITK_WARN("Statistics have not yet been calculated, running calculation now..."); CalculateStatistics(); } return m_Median; } HistogramStatisticsCalculator::MeasurementType HistogramStatisticsCalculator::GetUniformity() { if (!m_StatisticsCalculated) { MITK_WARN("Statistics have not yet been calculated, running calculation now..."); CalculateStatistics(); } return m_Uniformity; } HistogramStatisticsCalculator::MeasurementType HistogramStatisticsCalculator::GetUPP() { if (!m_StatisticsCalculated) { MITK_WARN("Statistics have not yet been calculated, running calculation now..."); CalculateStatistics(); } return m_UPP; } void HistogramStatisticsCalculator::CalculateStatistics() { if (m_Histogram.IsNull()) { throw std::runtime_error("Histogram not set in HistogramStatisticsCalculator::CalculateStatistics()"); } unsigned int nBins = m_Histogram->GetSize()[0]; m_Uniformity = 0; m_Entropy = 0; m_UPP = 0; m_Median = 0; MeasurementType cumulativeProbability = 0.0; MeasurementType partialProbability; bool medianFound(false); for (unsigned int i = 0; i < nBins; i++) { - partialProbability = m_Histogram->GetFrequency(i) / double( m_Histogram->GetTotalFrequency() ); + partialProbability = m_Histogram->GetFrequency(i, 0) / double( m_Histogram->GetTotalFrequency() ); cumulativeProbability += partialProbability; if (partialProbability != 0) { m_Entropy -= partialProbability * (std::log2( partialProbability )); m_Uniformity += std::pow(partialProbability, 2); if (m_Histogram->GetMeasurement(i, 0) > 0) { m_UPP += std::pow(partialProbability, 2); } } if (cumulativeProbability >= 0.5 && !medianFound) { MeasurementType binMin = m_Histogram->GetBinMin(0, i); MeasurementType binMax = m_Histogram->GetBinMax(0, i); m_Median = (binMax + binMin) / 2.0; medianFound = true; } } m_StatisticsCalculated = true; } } diff --git a/Modules/ImageStatistics/mitkHotspotCalculator.cpp b/Modules/ImageStatistics/mitkHotspotMaskGenerator.cpp similarity index 99% rename from Modules/ImageStatistics/mitkHotspotCalculator.cpp rename to Modules/ImageStatistics/mitkHotspotMaskGenerator.cpp index 84b8205616..bc14376d56 100644 --- a/Modules/ImageStatistics/mitkHotspotCalculator.cpp +++ b/Modules/ImageStatistics/mitkHotspotMaskGenerator.cpp @@ -1,535 +1,535 @@ -#include +#include #include #include #include #include #include "mitkImageAccessByItk.h" #include #include #include namespace mitk { HotspotMaskGenerator::HotspotMaskGenerator(): m_HotspotRadiusinMM(6.2035049089940), // radius of a 1cm3 sphere in mm m_HotspotMustBeCompletelyInsideImage(true), m_Label(1) { m_Modified=true; } void HotspotMaskGenerator::SetImage(mitk::Image::Pointer inputImage) { if (inputImage != m_Image) { m_Image = inputImage; m_Modified = true; } } void HotspotMaskGenerator::SetMask(MaskGenerator::Pointer mask) { if (mask != m_Mask) { m_Mask = mask; m_Modified = true; } } HotspotMaskGenerator::~HotspotMaskGenerator() { } void HotspotMaskGenerator::SetHotspotRadiusInMM(double radiusInMillimeter) { if(radiusInMillimeter != m_HotspotRadiusinMM) { m_HotspotRadiusinMM = radiusInMillimeter; m_Modified = true; } } const double& HotspotMaskGenerator::GetHotspotRadiusinMM() const { return m_HotspotRadiusinMM; } bool HotspotMaskGenerator::GetHotspotMustBeCompletelyInsideImage() const { return m_HotspotMustBeCompletelyInsideImage; } mitk::Image::Pointer HotspotMaskGenerator::GetMask() { if (m_Modified) { if ( m_Image.IsNull() ) { throw std::runtime_error( "Error: image empty!" ); } if ( m_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( 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 ) { 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 { 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_InternalMask; } void HotspotMaskGenerator::SetLabel(unsigned int label) { if (label != m_Label) { m_Label = label; m_Modified = true; } } template 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 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 > 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 > 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 return convolutionImage; } template < typename TPixel, unsigned int VImageDimension> void 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 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_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_InternalMask = hotspotMaskAsMITKImage; } } } diff --git a/Modules/ImageStatistics/mitkHotspotCalculator.h b/Modules/ImageStatistics/mitkHotspotMaskGenerator.h similarity index 100% rename from Modules/ImageStatistics/mitkHotspotCalculator.h rename to Modules/ImageStatistics/mitkHotspotMaskGenerator.h diff --git a/Modules/ImageStatistics/mitkBinaryImageMaskGenerator.cpp b/Modules/ImageStatistics/mitkImageMaskGenerator.cpp similarity index 57% rename from Modules/ImageStatistics/mitkBinaryImageMaskGenerator.cpp rename to Modules/ImageStatistics/mitkImageMaskGenerator.cpp index 84ac7420a9..57b4735f79 100644 --- a/Modules/ImageStatistics/mitkBinaryImageMaskGenerator.cpp +++ b/Modules/ImageStatistics/mitkImageMaskGenerator.cpp @@ -1,40 +1,46 @@ -#include +#include #include #include namespace mitk { -void BinaryImageMaskGenerator::SetImageMask(Image::Pointer maskImage) +void ImageMaskGenerator::SetImageMask(Image::Pointer maskImage) { if (m_internalMaskImage != maskImage) { m_internalMaskImage = maskImage; m_Modified = true; } } -mitk::Image::Pointer BinaryImageMaskGenerator::GetMask() +mitk::Image::Pointer ImageMaskGenerator::GetMask() { if (m_Modified) { - // extract timeStep -> is that correct? + unsigned int timeStepForExtraction; + if (m_TimeStep >= m_internalMaskImage->GetTimeSteps()) { - throw std::runtime_error("Error: Invalid time step in BinaryImageMaskGenerator!"); + MITK_WARN << "Warning: time step > number of time steps in mask image"; + timeStepForExtraction = m_internalMaskImage->GetTimeSteps() - 1; + } + else + { + timeStepForExtraction = m_TimeStep; } ImageTimeSelector::Pointer imageTimeSelector = ImageTimeSelector::New(); imageTimeSelector->SetInput(m_internalMaskImage); - imageTimeSelector->SetTimeNr(m_TimeStep); + imageTimeSelector->SetTimeNr(timeStepForExtraction); 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/mitkImageMaskGenerator.h similarity index 81% rename from Modules/ImageStatistics/mitkBinaryImageMaskGenerator.h rename to Modules/ImageStatistics/mitkImageMaskGenerator.h index 838377f5e1..8fd381d636 100644 --- a/Modules/ImageStatistics/mitkBinaryImageMaskGenerator.h +++ b/Modules/ImageStatistics/mitkImageMaskGenerator.h @@ -1,42 +1,42 @@ #ifndef mitkBinaryMaskGenerator #define mitkBinaryMaskGenerator #include #include #include #include #include namespace mitk { -class MITKIMAGESTATISTICS_EXPORT BinaryImageMaskGenerator: public MaskGenerator +class MITKIMAGESTATISTICS_EXPORT ImageMaskGenerator: public MaskGenerator { public: /** Standard Self typedef */ - typedef BinaryImageMaskGenerator Self; + typedef ImageMaskGenerator 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(){} + ImageMaskGenerator():Superclass(){} private: mitk::Image::Pointer m_internalMaskImage; }; } #endif diff --git a/Modules/ImageStatistics/mitkImageStatisticsCalculator2.cpp b/Modules/ImageStatistics/mitkImageStatisticsCalculator2.cpp index cc2f1af0c1..26e596b992 100644 --- a/Modules/ImageStatistics/mitkImageStatisticsCalculator2.cpp +++ b/Modules/ImageStatistics/mitkImageStatisticsCalculator2.cpp @@ -1,659 +1,657 @@ #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() = 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 ) { 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( 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) { 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->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() = 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(); } }