diff --git a/Modules/DiffusionImaging/MiniApps/ImageStatisticsMiniapp_v2.cpp b/Modules/DiffusionImaging/MiniApps/ImageStatisticsMiniapp_v2.cpp index e4a44899cd..57b65c87b4 100644 --- a/Modules/DiffusionImaging/MiniApps/ImageStatisticsMiniapp_v2.cpp +++ b/Modules/DiffusionImaging/MiniApps/ImageStatisticsMiniapp_v2.cpp @@ -1,200 +1,205 @@ /*=================================================================== 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; } + std::cout << 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::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); + calculator->SetNBinsForHistogramStatistics(5); + + std::cout << "planarfigure: " << std::endl; try { // result = calculator->GetStatistics(timeStep, 1); } catch( const itk::ExceptionObject& e) { MITK_ERROR << "Statistic Calculation Failed - ITK Exception:" << e.what(); return -1; } // printMap(result); -// calculator->SetMask(binaryImageMaskGen.GetPointer()); -// result = calculator->GetStatistics(timeStep, 1); -// printMap(result); - - calculator->SetMask(nullptr); - result = calculator->GetStatistics(timeStep); + std::cout << "image mask: " << std::endl; + calculator->SetMask(binaryImageMaskGen.GetPointer()); + result = calculator->GetStatistics(timeStep, 1); printMap(result); +// std::cout << "no mask: " << std::endl; +// calculator->SetMask(nullptr); +// result = calculator->GetStatistics(timeStep); +// 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->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; +// 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 5c078b083a..d76825a6f1 100644 --- a/Modules/ImageStatistics/files.cmake +++ b/Modules/ImageStatistics/files.cmake @@ -1,28 +1,30 @@ set(CPP_FILES mitkImageStatisticsCalculator.cpp mitkImageStatisticsCalculator2.cpp mitkPointSetStatisticsCalculator.cpp mitkPointSetDifferenceStatisticsCalculator.cpp mitkIntensityProfile.cpp mitkHotspotMaskGenerator.cpp mitkMaskGenerator.cpp mitkPlanarFigureMaskGenerator.cpp mitkMultiLabelMaskGenerator.cpp mitkImageMaskGenerator.cpp mitkHistogramStatisticsCalculator.cpp + mitkMaskUtilities.cpp ) set(H_FILES mitkImageStatisticsCalculator.h mitkImageStatisticsCalculator2.h mitkPointSetDifferenceStatisticsCalculator.h mitkPointSetStatisticsCalculator.h mitkExtendedStatisticsImageFilter.h mitkExtendedLabelStatisticsImageFilter.h mitkHotspotMaskGenerator.h mitkMaskGenerator.h mitkPlanarFigureMaskGenerator.h mitkMultiLabelMaskGenerator.h mitkImageMaskGenerator.h mitkHistogramStatisticsCalculator.h + mitkMaskUtilities.h ) diff --git a/Modules/ImageStatistics/mitkExtendedLabelStatisticsImageFilter.hxx b/Modules/ImageStatistics/mitkExtendedLabelStatisticsImageFilter.hxx index 42865e425c..b7e26546b7 100644 --- a/Modules/ImageStatistics/mitkExtendedLabelStatisticsImageFilter.hxx +++ b/Modules/ImageStatistics/mitkExtendedLabelStatisticsImageFilter.hxx @@ -1,324 +1,323 @@ /*=================================================================== 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. ===================================================================*/ #ifndef _mitkExtendedLabelStatisticsImageFilter_hxx #define _mitkExtendedLabelStatisticsImageFilter_hxx #include "mitkExtendedLabelStatisticsImageFilter.h" #include "itkImageRegionConstIteratorWithIndex.h" #include "itkImageRegionConstIterator.h" #include #include #include "mitkNumericConstants.h" #include "mitkLogMacros.h" namespace itk { template< class TInputImage , class TLabelImage> ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > ::ExtendedLabelStatisticsImageFilter() : LabelStatisticsImageFilter< TInputImage, TLabelImage >() { CalculateSettingsForLabels(); } template< class TInputImage, class TLabelImage > std::list ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > ::GetRelevantLabels() const { return m_RelevantLabels; } template< class TInputImage, class TLabelImage > bool ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > ::GetMaskingNonEmpty() const { return m_MaskNonEmpty; } template< class TInputImage, class TLabelImage > typename ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >::RealType ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > ::GetUniformity(LabelPixelType label) const { CoefficientsMapConstIterator mapIt; mapIt = m_LabelStatisticsCoefficients.find(label); if ( mapIt == m_LabelStatisticsCoefficients.end() ) { // label does not exist, return a default value return NumericTraits< PixelType >::Zero; } else { return ( *mapIt ).second.m_Uniformity; } } template< class TInputImage, class TLabelImage > typename ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >::RealType ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > ::GetEntropy(LabelPixelType label) const { CoefficientsMapConstIterator mapIt; mapIt = m_LabelStatisticsCoefficients.find(label); if ( mapIt == m_LabelStatisticsCoefficients.end() ) { // label does not exist, return a default value return NumericTraits< PixelType >::Zero; } else { return ( *mapIt ).second.m_Entropy; } } template< class TInputImage, class TLabelImage > typename ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >::RealType ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > ::GetUPP(LabelPixelType label) const { CoefficientsMapConstIterator mapIt; mapIt = m_LabelStatisticsCoefficients.find(label); if ( mapIt == m_LabelStatisticsCoefficients.end() ) { // label does not exist, return a default value return NumericTraits< PixelType >::Zero; } else { return ( *mapIt ).second.m_UPP; } } template< class TInputImage, class TLabelImage > typename ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >::RealType ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > ::GetMPP(LabelPixelType label) const { CoefficientsMapConstIterator mapIt; mapIt = m_LabelStatisticsCoefficients.find(label); if ( mapIt == m_LabelStatisticsCoefficients.end() ) { // label does not exist, return a default value return NumericTraits< PixelType >::Zero; } else { return ( *mapIt ).second.m_MPP; } } template< class TInputImage, class TLabelImage > typename ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >::RealType ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > ::GetKurtosis(LabelPixelType label) const { CoefficientsMapConstIterator mapIt; mapIt = m_LabelStatisticsCoefficients.find(label); if ( mapIt == m_LabelStatisticsCoefficients.end() ) { // label does not exist, return a default value return NumericTraits< PixelType >::Zero; } else { return ( *mapIt ).second.m_Kurtosis; } } template< class TInputImage, class TLabelImage > typename ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >::RealType ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > ::GetSkewness(LabelPixelType label) const { CoefficientsMapConstIterator mapIt; mapIt = m_LabelStatisticsCoefficients.find(label); if ( mapIt == m_LabelStatisticsCoefficients.end() ) { // label does not exist, return a default value return NumericTraits< PixelType >::Zero; } else { return ( *mapIt ).second.m_Skewness; } } template< class TInputImage, class TLabelImage > void ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >:: CalculateSettingsForLabels() { LabelPixelType i; m_MaskNonEmpty = false; for ( i = 1; i < 4096; ++i ) { if ( this->HasLabel( i ) ) { m_RelevantLabels.push_back( i ); m_MaskNonEmpty = true; m_LabelStatisticsCoefficients.insert( std::make_pair(i, CoefficientsClass()) ); } } } template< class TInputImage, class TLabelImage > void ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >:: ComputeEntropyUniformityAndUPP() { - double baseChange = std::log10(2); RealType partialProbability( 0.0 ); RealType uniformity( 0.0 ); RealType entropy( 0.0 ); RealType upp( 0.0 ); LabelPixelType i; if ( m_MaskNonEmpty ) { typename std::list< int >::const_iterator it; for ( it = m_RelevantLabels.cbegin(), i = 0; it != m_RelevantLabels.cend(); ++it, ++i ) { HistogramType::Pointer histogramForEntropy = this->GetHistogram(*it); for (int i = 0; i < histogramForEntropy->Size(); i++) { partialProbability = histogramForEntropy->GetFrequency(i,0) / double ( histogramForEntropy->GetTotalFrequency() ) ; if( partialProbability != 0) { entropy -= partialProbability *( std::log10(partialProbability) / std::log10(2) ) ; uniformity += std::pow(partialProbability,2); if(histogramForEntropy->GetMeasurement(i,0) > 0) { upp += std::pow(partialProbability,2); } } } m_LabelStatisticsCoefficients[*it].m_Entropy = entropy; m_LabelStatisticsCoefficients[*it].m_Uniformity = uniformity; m_LabelStatisticsCoefficients[*it].m_UPP = upp; } } } template< class TInputImage, class TLabelImage > void ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >:: ComputeSkewnessKurtosisAndMPP() { typename TLabelImage::RegionType Subregion; RealType baseOfSkewnessAndCurtosis( 0.0 ); RealType kurtosis( 0.0 ); RealType skewness( 0.0 ); RealType mpp( 0.0 ); RealType currentPixel( 0.0 ); std::list< LabelPixelType> relevantLabels; LabelPixelType i; if ( m_MaskNonEmpty ) { typename std::list< int >::const_iterator it; for ( it = m_RelevantLabels.cbegin(), i = 0; it != m_RelevantLabels.cend(); ++it ) { RealType sigma = this->GetSigma( *it ); RealType mean = this->GetMean( *it ); Subregion = Superclass::GetRegion(*it); int count( this->GetCount(*it) ); if ( count == 0 || sigma < mitk::eps) { throw std::logic_error( "Empty segmentation" ); } if ( fabs( sigma ) < mitk::sqrteps ) { throw std::logic_error( "Sigma == 0" ); } ImageRegionConstIteratorWithIndex< TInputImage > it1 (this->GetInput(), Subregion); ImageRegionConstIterator< TLabelImage > labelIt (this->GetLabelInput(), Subregion); for (it1.GoToBegin(); !it1.IsAtEnd(); ++it1, ++labelIt) { if (labelIt.Get() == *it) { currentPixel = it1.Get(); baseOfSkewnessAndCurtosis = (currentPixel -mean) / sigma; kurtosis += std::pow( baseOfSkewnessAndCurtosis, 4.0 ); skewness += std::pow( baseOfSkewnessAndCurtosis, 3.0 ); if(currentPixel > 0) { mpp+= currentPixel; } } } m_LabelStatisticsCoefficients[*it].m_Skewness = RealType(skewness/count); m_LabelStatisticsCoefficients[*it].m_Kurtosis = RealType(kurtosis/count); m_LabelStatisticsCoefficients[*it].m_MPP = RealType(mpp/count); } } } template< class TInputImage, class TLabelImage > void ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >:: AfterThreadedGenerateData() { Superclass::AfterThreadedGenerateData(); CalculateSettingsForLabels(); ComputeSkewnessKurtosisAndMPP(); if(this->GetUseHistograms()) { ComputeEntropyUniformityAndUPP(); } else { MITK_WARN << "Cannot compute coefficients UPP,Entropy,Uniformity because of missing histogram"; } } } // end namespace itk #endif diff --git a/Modules/ImageStatistics/mitkExtendedStatisticsImageFilter.hxx b/Modules/ImageStatistics/mitkExtendedStatisticsImageFilter.hxx index 004cf1a90b..3f3d810edb 100644 --- a/Modules/ImageStatistics/mitkExtendedStatisticsImageFilter.hxx +++ b/Modules/ImageStatistics/mitkExtendedStatisticsImageFilter.hxx @@ -1,344 +1,343 @@ /*=================================================================== 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. ===================================================================*/ #ifndef __mitkExtendedStatisticsImageFilter_hxx #define __mitkExtendedStatisticsImageFilter_hxx #include "mitkExtendedStatisticsImageFilter.h" namespace itk { template< class TInputImage > ExtendedStatisticsImageFilter< TInputImage >::ExtendedStatisticsImageFilter() : StatisticsImageFilter< TInputImage >() { /* * add the Skewness,Kurtosis,Entropy,Uniformity,MPP, UPP, Median to the other statistical calculated Values * of the mitkStatisticsImageFilter as the 7th to the 13th Output */ for ( int i = 7; i < 14; ++i ) { typename RealObjectType::Pointer output = static_cast< RealObjectType * >( this->MakeOutput(i).GetPointer() ); this->ProcessObject::SetNthOutput( i, output.GetPointer() ); } this->GetSkewnessOutput()->Set( 0.0 ); this->GetKurtosisOutput()->Set( 0.0 ); this->GetEntropyOutput()->Set( -1.0 ); this->GetUniformityOutput()->Set( 0.0 ); this->GetUPPOutput()->Set( 0.0 ); this->GetMPPOutput()->Set( 0.0 ); this->GetMedianOutput()->Set( 0.0 ); this->m_HistogramCalculated = false; this->m_HistogramGenerator = HistogramGeneratorType::New(); } template< class TInputImage > DataObject::Pointer ExtendedStatisticsImageFilter< TInputImage >::MakeOutput( ProcessObject::DataObjectPointerArraySizeType output) { switch ( output ) { case 7: case 8: case 9: case 10: case 11: case 12: { return RealObjectType::New().GetPointer(); break; } case 13: { return RealObjectType::New().GetPointer(); break; } default: { // might as well make an image return Superclass::MakeOutput( output ); break; } } } template< class TInputImage > void ExtendedStatisticsImageFilter< TInputImage > ::SetBinSize( int size ) { m_BinSize = size; } template< class TInputImage > typename ExtendedStatisticsImageFilter< TInputImage >::RealObjectType * ExtendedStatisticsImageFilter< TInputImage > ::GetSkewnessOutput() { return static_cast< RealObjectType * >( this->ProcessObject::GetOutput(7) ); } template< class TInputImage > const typename ExtendedStatisticsImageFilter< TInputImage >::RealObjectType * ExtendedStatisticsImageFilter< TInputImage > ::GetSkewnessOutput() const { return static_cast< const RealObjectType * >( this->ProcessObject::GetOutput(7) ); } template< class TInputImage > typename ExtendedStatisticsImageFilter< TInputImage >::RealObjectType * ExtendedStatisticsImageFilter< TInputImage > ::GetKurtosisOutput() { return static_cast< RealObjectType * >( this->ProcessObject::GetOutput(8) ); } template< class TInputImage > const typename ExtendedStatisticsImageFilter< TInputImage >::RealObjectType * ExtendedStatisticsImageFilter< TInputImage > ::GetKurtosisOutput() const { return static_cast< const RealObjectType * >( this->ProcessObject::GetOutput(8) ); } template< class TInputImage > typename ExtendedStatisticsImageFilter< TInputImage >::RealObjectType * ExtendedStatisticsImageFilter< TInputImage > ::GetUniformityOutput() { return static_cast< RealObjectType * >( this->ProcessObject::GetOutput(9) ); } template< class TInputImage > const typename ExtendedStatisticsImageFilter< TInputImage >::RealObjectType * ExtendedStatisticsImageFilter< TInputImage > ::GetUniformityOutput() const { return static_cast< const RealObjectType * >( this->ProcessObject::GetOutput(9) ); } template< class TInputImage > typename ExtendedStatisticsImageFilter< TInputImage >::RealObjectType * ExtendedStatisticsImageFilter< TInputImage > ::GetEntropyOutput() { return static_cast< RealObjectType * >( this->ProcessObject::GetOutput(10) ); } template< class TInputImage > const typename ExtendedStatisticsImageFilter< TInputImage >::RealObjectType * ExtendedStatisticsImageFilter< TInputImage > ::GetEntropyOutput() const { return static_cast< const RealObjectType * >( this->ProcessObject::GetOutput(10) ); } template< class TInputImage > typename ExtendedStatisticsImageFilter< TInputImage >::RealObjectType * ExtendedStatisticsImageFilter< TInputImage > ::GetUPPOutput() { return static_cast< RealObjectType * >( this->ProcessObject::GetOutput(11) ); } template< class TInputImage > const typename ExtendedStatisticsImageFilter< TInputImage >::RealObjectType * ExtendedStatisticsImageFilter< TInputImage > ::GetUPPOutput() const { return static_cast< const RealObjectType * >( this->ProcessObject::GetOutput(11) ); } template< class TInputImage > typename ExtendedStatisticsImageFilter< TInputImage >::RealObjectType * ExtendedStatisticsImageFilter< TInputImage > ::GetMPPOutput() { return static_cast< RealObjectType * >( this->ProcessObject::GetOutput(12) ); } template< class TInputImage > const typename ExtendedStatisticsImageFilter< TInputImage >::RealObjectType * ExtendedStatisticsImageFilter< TInputImage > ::GetMPPOutput() const { return static_cast< const RealObjectType * >( this->ProcessObject::GetOutput(12) ); } template< class TInputImage > typename ExtendedStatisticsImageFilter< TInputImage >::RealObjectType * ExtendedStatisticsImageFilter< TInputImage > ::GetMedianOutput() { return static_cast< RealObjectType * >( this->ProcessObject::GetOutput(13) ); } template< class TInputImage > const typename ExtendedStatisticsImageFilter< TInputImage >::RealObjectType * ExtendedStatisticsImageFilter< TInputImage > ::GetMedianOutput() const { return static_cast< const RealObjectType * >( this->ProcessObject::GetOutput(13) ); } template< class TInputImage > void ExtendedStatisticsImageFilter< TInputImage > ::AfterThreadedGenerateData() { Superclass::AfterThreadedGenerateData(); ComputeSkewnessKurtosisAndMPP(); CalculateHistogram(); if(m_HistogramCalculated == true) { ComputeEntropyUniformityMedianAndUPP(); } else { MITK_WARN << "Cannot compute coefficients UPP,Entropy,Uniformity because of missing histogram"; } } template< class TInputImage > void ExtendedStatisticsImageFilter< TInputImage > ::ComputeSkewnessKurtosisAndMPP() { RealType mean = this->GetMean(); RealType sigma = this->GetSigma(); RealType baseOfSkewnessAndKurtosis; RealType kurtosis(0.0); RealType skewness(0.0); RealType mpp(0.0); RealType currentPixel(0.0); if ( sigma < mitk::eps ) { throw std::logic_error( "Empty segmentation" ); } ImageRegionConstIterator< TInputImage > it (this->GetInput(), this->GetInput()->GetLargestPossibleRegion() ); int counter = 0; for (it.GoToBegin(); !it.IsAtEnd(); ++it) { currentPixel = it.Get(); baseOfSkewnessAndKurtosis = (currentPixel - mean) / sigma ; kurtosis += std::pow( baseOfSkewnessAndKurtosis, 4.0 ); skewness += std::pow( baseOfSkewnessAndKurtosis, 3.0 ); if(currentPixel > 0) { mpp+= currentPixel; } counter++; } if ( counter == 0 ) { throw std::logic_error( "Empty segmentation" ); } kurtosis = kurtosis / counter; skewness = skewness / counter; mpp = mpp/counter; this->GetKurtosisOutput()->Set( kurtosis ); this->GetSkewnessOutput()->Set( skewness ); this->GetMPPOutput()->Set( mpp ); } template< class TInputImage > void ExtendedStatisticsImageFilter< TInputImage > ::CalculateHistogram() { m_HistogramGenerator->SetInput( this->GetInput() ); m_HistogramGenerator->SetMarginalScale( 100 ); m_HistogramGenerator->SetNumberOfBins( m_BinSize ); m_HistogramGenerator->SetHistogramMin( this->GetMinimum() ); m_HistogramGenerator->SetHistogramMax( this->GetMaximum() ); m_HistogramGenerator->Compute(); m_HistogramCalculated = true; } template< class TInputImage > void ExtendedStatisticsImageFilter< TInputImage > ::ComputeEntropyUniformityMedianAndUPP() { - double baseChange = std::log10(2); RealType partialProbability( 0.0 ); RealType uniformity( 0.0 ); RealType entropy( 0.0 ); RealType upp( 0.0 ); RealType median( 0.0 ); const typename HistogramGeneratorType::HistogramType* histogramForEntropy = GetHistogram(); double cumulativeProbability = 0.0; bool medianFound = false; for (int i = 0; i < histogramForEntropy->Size(); i++) { partialProbability = histogramForEntropy->GetFrequency(i,0) / double ( histogramForEntropy->GetTotalFrequency() ) ; cumulativeProbability += double ( partialProbability ); if( partialProbability != 0) { entropy -= partialProbability *( std::log10(partialProbability) / std::log10(2) ) ; uniformity += std::pow(partialProbability,2); if(histogramForEntropy->GetMeasurement(i,0) > 0) { upp += std::pow(partialProbability,2); } } if( cumulativeProbability >= 0.5 && !medianFound ) { RealType binMin = histogramForEntropy->GetBinMin( 0, i ); RealType binMax = histogramForEntropy->GetBinMax( 0, i ); median = ( binMin + binMax ) / 2.0; medianFound = true; } } this->GetEntropyOutput()->Set( entropy ); this->GetUniformityOutput()->Set( uniformity ); this->GetUPPOutput()->Set( upp ); this->GetMedianOutput()->Set( median ); } } #endif diff --git a/Modules/ImageStatistics/mitkImageStatisticsCalculator2.cpp b/Modules/ImageStatistics/mitkImageStatisticsCalculator2.cpp index 26e596b992..3734371e12 100644 --- a/Modules/ImageStatistics/mitkImageStatisticsCalculator2.cpp +++ b/Modules/ImageStatistics/mitkImageStatisticsCalculator2.cpp @@ -1,657 +1,531 @@ #include #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()); - - } - 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; + typedef MaskUtilities< TPixel, VImageDimension > MaskUtilType; // 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); + typename MaskUtilType::Pointer maskUtil = MaskUtilType::New(); + maskUtil->SetImage(image); + maskUtil->SetMask(maskImage.GetPointer()); // 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); + adaptedImage = maskUtil->ExtractMaskImageRegion(); // this also checks mask sanity 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; 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 64452323bb..8e1c5ffd37 100644 --- a/Modules/ImageStatistics/mitkImageStatisticsCalculator2.h +++ b/Modules/ImageStatistics/mitkImageStatisticsCalculator2.h @@ -1,367 +1,357 @@ #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=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 - ); - 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/mitkMaskUtilities.cpp b/Modules/ImageStatistics/mitkMaskUtilities.cpp new file mode 100644 index 0000000000..38f82022ac --- /dev/null +++ b/Modules/ImageStatistics/mitkMaskUtilities.cpp @@ -0,0 +1,189 @@ +#ifndef MITKMASKUTIL_CPP +#define MITKMASKUTIL_CPP + +#include +//#include +#include +#include +#include +#include + +namespace mitk +{ + template + void MaskUtilities::SetImage(ImageType* image) + { + if (image != m_Image) + { + m_Image = image; + } + } + + template + void MaskUtilities::SetMask(MaskType* mask) + { + if (mask != m_Mask) + { + m_Mask = mask; + } + } + + template + bool MaskUtilities::CheckMaskSanity() + { + if (m_Mask==nullptr || m_Image==nullptr) + { + MITK_ERROR << "Set an image and a mask first"; + } + + typedef itk::Image< TPixel, VImageDimension > ImageType; + typedef typename ImageType::PointType PointType; + typedef typename ImageType::DirectionType DirectionType; + + bool maskSanity = true; + + + if (m_Mask==nullptr) + { + MITK_ERROR << "Something went wrong when casting the mitk mask image to an itk mask image. Do the mask and the input image have the same dimension?"; + // note to self: We could try to convert say a 2d mask to a 3d mask if the image is 3d. (mask and image dimension have to match.) + } + // check direction + DirectionType imageDirection = m_Image->GetDirection(); + DirectionType maskDirection = m_Mask->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) + { + maskSanity = false; + MITK_INFO << "Mask needs to have same direction as image! (Image direction: " << imageDirection << "; Mask direction: " << maskDirection << ")"; + } + } + } + } + + // check spacing + PointType imageSpacing = m_Image->GetSpacing(); + PointType maskSpacing = m_Mask->GetSpacing(); + for (unsigned int i = 0; i < VImageDimension; i++) + { + if ( fabs( maskSpacing[i] - imageSpacing[i] ) > mitk::eps ) + { + maskSanity = false; + MITK_INFO << "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 = m_Image->GetOrigin(); + PointType maskOrigin = m_Mask->GetOrigin(); + + typedef itk::ContinuousIndex ContinousIndexType; + ContinousIndexType maskOriginContinousIndex, imageOriginContinousIndex; + + m_Image->TransformPhysicalPointToContinuousIndex(maskOrigin, maskOriginContinousIndex); + m_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 ) + { + maskSanity = false; + MITK_INFO << "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 ( m_Mask!=nullptr && + !m_Image->GetLargestPossibleRegion().IsInside( m_Mask->GetLargestPossibleRegion() ) ) + { + maskSanity = false; + MITK_INFO << "Mask region needs to be inside of image region! (Image region: " + << m_Image->GetLargestPossibleRegion() << "; Mask region: " << m_Mask->GetLargestPossibleRegion() << ")"; + } + return maskSanity; + } + + template + typename itk::Image::Pointer MaskUtilities::ExtractMaskImageRegion() + { + if (m_Mask==nullptr || m_Image==nullptr) + { + MITK_ERROR << "Set an image and a mask first"; + } + + bool maskSanity = CheckMaskSanity(); + + if (!maskSanity) + { + MITK_ERROR << "Mask and image are not compatible"; + } + + typedef itk::Image< TPixel, VImageDimension > ImageType; + typedef itk::Image< unsigned short, VImageDimension > MaskType; + typedef itk::ExtractImageFilter< ImageType, ImageType > ExtractImageFilterType; + + typename ImageType::SizeType imageSize = m_Image->GetBufferedRegion().GetSize(); + typename ImageType::SizeType maskSize = m_Mask->GetBufferedRegion().GetSize(); + + typename itk::Image::Pointer extractedImg = itk::Image::New(); + + 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 MaskType::PointType maskOrigin = m_Mask->GetOrigin(); + typename ImageType::PointType imageOrigin = m_Image->GetOrigin(); + typename MaskType::SpacingType maskSpacing = m_Mask->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(m_Mask->GetLargestPossibleRegion().GetSize()); + + extractImageFilter->SetInput( m_Image ); + extractImageFilter->SetExtractionRegion( extractionRegion ); + extractImageFilter->SetCoordinateTolerance( 0.001 ); + extractImageFilter->SetDirectionTolerance( 0.001 ); + extractImageFilter->Update(); + extractedImg = extractImageFilter->GetOutput(); + extractedImg->SetOrigin(m_Mask->GetOrigin()); + extractedImg->SetLargestPossibleRegion(m_Mask->GetLargestPossibleRegion()); + extractedImg->SetBufferedRegion(m_Mask->GetBufferedRegion()); + + } + else + { + extractedImg = m_Image; + } + + return extractedImg; + } + +} + +#endif diff --git a/Modules/ImageStatistics/mitkMaskUtilities.h b/Modules/ImageStatistics/mitkMaskUtilities.h new file mode 100644 index 0000000000..7a2c4823ed --- /dev/null +++ b/Modules/ImageStatistics/mitkMaskUtilities.h @@ -0,0 +1,51 @@ +#ifndef MITKMASKUTIL +#define MITKMASKUTIL + +#include +#include +#include + +namespace mitk +{ +template +class MITKIMAGESTATISTICS_EXPORT MaskUtilities: public itk::Object + { + public: + /** Standard Self typedef */ + typedef MaskUtilities 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(MaskUtilities, itk::Object) + + typedef itk::Image ImageType; + typedef itk::Image MaskType; + + void SetImage(ImageType* image); + void SetMask(MaskType* mask); + + bool CheckMaskSanity(); + + typename itk::Image::Pointer ExtractMaskImageRegion(); + + protected: + MaskUtilities(){} + + ~MaskUtilities(){} + + private: + itk::Image* m_Image; + itk::Image* m_Mask; + }; +} + +#ifndef ITK_MANUAL_INSTANTIATION +#include +#endif + +#endif