diff --git a/Modules/ImageStatistics/files.cmake b/Modules/ImageStatistics/files.cmake index 23d6c94783..3d89557ffe 100644 --- a/Modules/ImageStatistics/files.cmake +++ b/Modules/ImageStatistics/files.cmake @@ -1,50 +1,50 @@ set(CPP_FILES mitkImageStatisticsCalculator.cpp mitkImageStatisticsContainer.cpp mitkPointSetStatisticsCalculator.cpp mitkPointSetDifferenceStatisticsCalculator.cpp mitkIntensityProfile.cpp mitkHotspotMaskGenerator.cpp mitkMaskGenerator.cpp mitkPlanarFigureMaskGenerator.cpp mitkMultiLabelMaskGenerator.cpp mitkImageMaskGenerator.cpp mitkHistogramStatisticsCalculator.cpp mitkIgnorePixelMaskGenerator.cpp mitkImageStatisticsPredicateHelper.cpp mitkImageStatisticsContainerNodeHelper.cpp mitkImageStatisticsContainerManager.cpp mitkStatisticsToImageRelationRule.cpp mitkStatisticsToMaskRelationRule.cpp mitkImageStatisticsConstants.cpp ) set(H_FILES mitkImageStatisticsCalculator.h mitkImageStatisticsContainer.h mitkPointSetDifferenceStatisticsCalculator.h mitkPointSetStatisticsCalculator.h mitkStatisticsImageFilter.h - mitkExtendedLabelStatisticsImageFilter.h + mitkLabelStatisticsImageFilter.h mitkHotspotMaskGenerator.h mitkMaskGenerator.h mitkPlanarFigureMaskGenerator.h mitkMultiLabelMaskGenerator.h mitkImageMaskGenerator.h mitkHistogramStatisticsCalculator.h mitkMaskUtilities.h mitkitkMaskImageFilter.h mitkIgnorePixelMaskGenerator.h mitkMinMaxImageFilterWithIndex.h mitkMinMaxLabelmageFilterWithIndex.h mitkImageStatisticsPredicateHelper.h mitkImageStatisticsContainerNodeHelper.h mitkImageStatisticsContainerManager.h mitkStatisticsToImageRelationRule.h mitkStatisticsToMaskRelationRule.h mitkImageStatisticsConstants.h ) set(TPP_FILES mitkMaskUtilities.tpp ) diff --git a/Modules/ImageStatistics/mitkExtendedLabelStatisticsImageFilter.h b/Modules/ImageStatistics/mitkExtendedLabelStatisticsImageFilter.h deleted file mode 100644 index 23081fa423..0000000000 --- a/Modules/ImageStatistics/mitkExtendedLabelStatisticsImageFilter.h +++ /dev/null @@ -1,345 +0,0 @@ -/*============================================================================ - -The Medical Imaging Interaction Toolkit (MITK) - -Copyright (c) German Cancer Research Center (DKFZ) -All rights reserved. - -Use of this source code is governed by a 3-clause BSD license that can be -found in the LICENSE file. - -============================================================================*/ -#ifndef __mitkExtendedLabelStatisticsImageFilter -#define __mitkExtendedLabelStatisticsImageFilter - -#include "itkLabelStatisticsImageFilter.h" - -namespace itk -{ - /** - * \class ExtendedLabelStatisticsImageFilter - * \brief Extension of the itkLabelStatisticsImageFilter that also calculates the Skewness,Kurtosis,Entropy,Uniformity. - * - * This class inherits from the itkLabelStatisticsImageFilter and - * uses its results for the calculation of seven additional coefficients: - * the Skewness, Kurtosis, Uniformity, UPP, MPP, Entropy and Median - * - * - */ - template< class TInputImage, class TLabelImage > - class ExtendedLabelStatisticsImageFilter : public LabelStatisticsImageFilter< TInputImage, TLabelImage > - { - public: - - typedef ExtendedLabelStatisticsImageFilter Self; - typedef LabelStatisticsImageFilter < TInputImage, TLabelImage > Superclass; - typedef SmartPointer< Self > Pointer; - typedef SmartPointer< const Self > ConstPointer; - typedef typename Superclass::LabelPixelType LabelPixelType; - typedef typename Superclass::RealType RealType; - typedef typename Superclass::PixelType PixelType; - typedef typename Superclass::MapIterator MapIterator; - typedef typename Superclass::BoundingBoxType BoundingBoxType; - typedef typename Superclass::RegionType RegionType; - typedef itk::Statistics::Histogram HistogramType; - - itkFactorylessNewMacro( Self ); - itkCloneMacro( Self ); - itkTypeMacro(ExtendedLabelStatisticsImageFilter, LabelStatisticsImageFilter); - - - /** \class LabelStatistics - * \brief Statistics stored per label - * \ingroup ITKImageStatistics - */ - class LabelStatistics - { - public: - - // default constructor - LabelStatistics() - { - // initialized to the default values - m_Count = NumericTraits< IdentifierType >::ZeroValue(); - m_PositivePixelCount = NumericTraits< IdentifierType >::ZeroValue(); - m_Sum = NumericTraits< RealType >::ZeroValue(); - m_SumOfPositivePixels = NumericTraits< RealType >::ZeroValue(); - - m_SumOfSquares = NumericTraits< RealType >::ZeroValue(); - m_SumOfCubes = NumericTraits< RealType >::ZeroValue(); - m_SumOfQuadruples = NumericTraits< RealType >::ZeroValue(); - - // Set such that the first pixel encountered can be compared - m_Minimum = NumericTraits< RealType >::max(); - m_Maximum = NumericTraits< RealType >::NonpositiveMin(); - - // Default these to zero - m_Mean = NumericTraits< RealType >::ZeroValue(); - m_Sigma = NumericTraits< RealType >::ZeroValue(); - m_Variance = NumericTraits< RealType >::ZeroValue(); - m_MPP = NumericTraits< RealType >::ZeroValue(); - m_Median = NumericTraits< RealType >::ZeroValue(); - m_Uniformity = NumericTraits< RealType >::ZeroValue(); - m_UPP = NumericTraits< RealType >::ZeroValue(); - m_Entropy = NumericTraits< RealType >::ZeroValue(); - m_Skewness = NumericTraits< RealType >::ZeroValue(); - m_Kurtosis = NumericTraits< RealType >::ZeroValue(); - - unsigned int imageDimension = itkGetStaticConstMacro(ImageDimension); - m_BoundingBox.resize(imageDimension * 2); - for ( unsigned int i = 0; i < imageDimension * 2; i += 2 ) - { - m_BoundingBox[i] = NumericTraits< IndexValueType >::max(); - m_BoundingBox[i + 1] = NumericTraits< IndexValueType >::NonpositiveMin(); - } - m_Histogram = nullptr; - } - - // constructor with histogram enabled - LabelStatistics(int size, RealType lowerBound, RealType upperBound) - { - // initialized to the default values - m_Count = NumericTraits< IdentifierType >::ZeroValue(); - m_PositivePixelCount = NumericTraits< IdentifierType >::ZeroValue(); - m_Sum = NumericTraits< RealType >::ZeroValue(); - m_SumOfPositivePixels = NumericTraits< RealType >::ZeroValue(); - - m_SumOfSquares = NumericTraits< RealType >::ZeroValue(); - m_SumOfCubes = NumericTraits< RealType >::ZeroValue(); - m_SumOfQuadruples = NumericTraits< RealType >::ZeroValue(); - - // Set such that the first pixel encountered can be compared - m_Minimum = NumericTraits< RealType >::max(); - m_Maximum = NumericTraits< RealType >::NonpositiveMin(); - - // Default these to zero - m_Mean = NumericTraits< RealType >::ZeroValue(); - m_Sigma = NumericTraits< RealType >::ZeroValue(); - m_Variance = NumericTraits< RealType >::ZeroValue(); - m_MPP = NumericTraits< RealType >::ZeroValue(); - m_Median = NumericTraits< RealType >::ZeroValue(); - m_Uniformity = NumericTraits< RealType >::ZeroValue(); - m_UPP = NumericTraits< RealType >::ZeroValue(); - m_Entropy = NumericTraits< RealType >::ZeroValue(); - m_Skewness = NumericTraits< RealType >::ZeroValue(); - m_Kurtosis = NumericTraits< RealType >::ZeroValue(); - - - unsigned int imageDimension = itkGetStaticConstMacro(ImageDimension); - m_BoundingBox.resize(imageDimension * 2); - for ( unsigned int i = 0; i < imageDimension * 2; i += 2 ) - { - m_BoundingBox[i] = NumericTraits< IndexValueType >::max(); - m_BoundingBox[i + 1] = NumericTraits< IndexValueType >::NonpositiveMin(); - } - - // Histogram - m_Histogram = HistogramType::New(); - typename HistogramType::SizeType hsize; - typename HistogramType::MeasurementVectorType lb; - typename HistogramType::MeasurementVectorType ub; - hsize.SetSize(1); - lb.SetSize(1); - ub.SetSize(1); - m_Histogram->SetMeasurementVectorSize(1); - hsize[0] = size; - lb[0] = lowerBound; - ub[0] = upperBound; - m_Histogram->Initialize(hsize, lb, ub); - } - - // need copy constructor because of smart pointer to histogram - LabelStatistics(const LabelStatistics & l) - { - m_Count = l.m_Count; - m_Minimum = l.m_Minimum; - m_Maximum = l.m_Maximum; - m_Mean = l.m_Mean; - m_Sum = l.m_Sum; - m_SumOfSquares = l.m_SumOfSquares; - m_Sigma = l.m_Sigma; - m_Variance = l.m_Variance; - m_MPP = l.m_MPP; - m_Median = l.m_Median; - m_Uniformity = l.m_Uniformity; - m_UPP = l.m_UPP; - m_Entropy = l.m_Entropy; - m_Skewness = l.m_Skewness; - m_Kurtosis = l.m_Kurtosis; - m_BoundingBox = l.m_BoundingBox; - m_Histogram = l.m_Histogram; - m_SumOfPositivePixels = l.m_SumOfPositivePixels; - m_PositivePixelCount = l.m_PositivePixelCount; - m_SumOfCubes = l.m_SumOfCubes; - m_SumOfQuadruples = l.m_SumOfQuadruples; - } - - // added for completeness - LabelStatistics &operator= (const LabelStatistics& l) - { - if(this != &l) - { - m_Count = l.m_Count; - m_Minimum = l.m_Minimum; - m_Maximum = l.m_Maximum; - m_Mean = l.m_Mean; - m_Sum = l.m_Sum; - m_SumOfSquares = l.m_SumOfSquares; - m_Sigma = l.m_Sigma; - m_Variance = l.m_Variance; - m_MPP = l.m_MPP; - m_Median = l.m_Median; - m_Uniformity = l.m_Uniformity; - m_UPP = l.m_UPP; - m_Entropy = l.m_Entropy; - m_Skewness = l.m_Skewness; - m_Kurtosis = l.m_Kurtosis; - m_BoundingBox = l.m_BoundingBox; - m_Histogram = l.m_Histogram; - m_SumOfPositivePixels = l.m_SumOfPositivePixels; - m_PositivePixelCount = l.m_PositivePixelCount; - m_SumOfCubes = l.m_SumOfCubes; - m_SumOfQuadruples = l.m_SumOfQuadruples; - } - return *this; - } - - IdentifierType m_Count; - RealType m_Minimum; - RealType m_Maximum; - RealType m_Mean; - RealType m_Sum; - RealType m_SumOfSquares; - RealType m_Sigma; - RealType m_Variance; - RealType m_MPP; - RealType m_Median; - RealType m_Uniformity; - RealType m_UPP; - RealType m_Entropy; - RealType m_Skewness; - RealType m_Kurtosis; - IdentifierType m_PositivePixelCount; - RealType m_SumOfPositivePixels; - RealType m_SumOfCubes; - RealType m_SumOfQuadruples; - typename Superclass::BoundingBoxType m_BoundingBox; - typename HistogramType::Pointer m_Histogram; - }; - - /** Type of the map used to store data per label */ - typedef std::unordered_map< LabelPixelType, LabelStatistics > MapType; - typedef typename std::unordered_map< LabelPixelType, LabelStatistics >::const_iterator StatisticsMapConstIterator; - typedef typename std::unordered_map< LabelPixelType, LabelStatistics >::iterator StatisticsMapIterator; - typedef IdentifierType MapSizeType; - - /** Type of the container used to store valid label values */ - typedef std::vector ValidLabelValuesContainerType; - - /** Return the computed Minimum for a label. */ - RealType GetMinimum(LabelPixelType label) const; - - /** Return the computed Maximum for a label. */ - RealType GetMaximum(LabelPixelType label) const; - - /** Return the computed Mean for a label. */ - RealType GetMean(LabelPixelType label) const; - - /** Return the computed Standard Deviation for a label. */ - RealType GetSigma(LabelPixelType label) const; - - /** Return the computed Variance for a label. */ - RealType GetVariance(LabelPixelType label) const; - - /** Return the computed bounding box for a label. */ - BoundingBoxType GetBoundingBox(LabelPixelType label) const; - - /** Return the computed region. */ - RegionType GetRegion(LabelPixelType label) const; - - /** Return the compute Sum for a label. */ - RealType GetSum(LabelPixelType label) const; - - /** Return the number of pixels for a label. */ - MapSizeType GetCount(LabelPixelType label) const; - - /** Return the histogram for a label */ - HistogramType::Pointer GetHistogram(LabelPixelType label) const; - - /*getter method for the new statistics*/ - RealType GetSkewness(LabelPixelType label) const; - RealType GetKurtosis(LabelPixelType label) const; - RealType GetUniformity( LabelPixelType label) const; - RealType GetMedian( LabelPixelType label) const; - RealType GetEntropy( LabelPixelType label) const; - RealType GetMPP( LabelPixelType label) const; - RealType GetUPP( LabelPixelType label) const; - - bool GetMaskingNonEmpty() const; - - std::list GetRelevantLabels() const; - - - /** specify global Histogram parameters. If the histogram parameters are set with this function, the same min and max value are used for all histograms. */ - void SetHistogramParameters(const int numBins, RealType lowerBound, - RealType upperBound); - - /** specify Histogram parameters for each label individually. Labels in the label image that are not represented in the std::maps here will receive global parameters (if available) */ - void SetHistogramParametersForLabels(std::map numBins, std::map lowerBound, - std::map upperBound); - - protected: - ExtendedLabelStatisticsImageFilter(): - m_GlobalHistogramParametersSet(false), - m_MaskNonEmpty(false), - m_LabelHistogramParametersSet(false), - m_PreferGlobalHistogramParameters(false) - { - m_NumBins.set_size(1); - } - - ~ExtendedLabelStatisticsImageFilter() override{} - - void AfterStreamedGenerateData() override; - - /** Initialize some accumulators before the threads run. */ - void BeforeStreamedGenerateData() override; - - /** Multi-thread version GenerateData. */ - void ThreadedStreamedGenerateData(const RegionType& region) override; - - /** Does the specified label exist? Can only be called after a call - * a call to Update(). */ - bool HasLabel(LabelPixelType label) const - { - return m_LabelStatistics.find(label) != m_LabelStatistics.end(); - } - - private: - std::vector< MapType > m_LabelStatisticsPerThread; - MapType m_LabelStatistics; - ValidLabelValuesContainerType m_ValidLabelValues; - - bool m_GlobalHistogramParametersSet; - - typename HistogramType::SizeType m_NumBins; - - RealType m_LowerBound; - RealType m_UpperBound; - - bool m_MaskNonEmpty; - - bool m_LabelHistogramParametersSet; - std::map m_LabelMin, m_LabelMax; - std::map m_LabelNBins; - bool m_PreferGlobalHistogramParameters; - - }; // end of class - -} // end namespace itk - -#ifndef ITK_MANUAL_INSTANTIATION -#include "mitkExtendedLabelStatisticsImageFilter.hxx" -#endif - -#endif diff --git a/Modules/ImageStatistics/mitkExtendedLabelStatisticsImageFilter.hxx b/Modules/ImageStatistics/mitkExtendedLabelStatisticsImageFilter.hxx deleted file mode 100644 index f2ecd243d6..0000000000 --- a/Modules/ImageStatistics/mitkExtendedLabelStatisticsImageFilter.hxx +++ /dev/null @@ -1,734 +0,0 @@ -/*============================================================================ - -The Medical Imaging Interaction Toolkit (MITK) - -Copyright (c) German Cancer Research Center (DKFZ) -All rights reserved. - -Use of this source code is governed by a 3-clause BSD license that can be -found in the LICENSE file. - -============================================================================*/ -#ifndef _mitkExtendedLabelStatisticsImageFilter_hxx -#define _mitkExtendedLabelStatisticsImageFilter_hxx - -#include "mitkExtendedLabelStatisticsImageFilter.h" - -#include "itkImageRegionConstIteratorWithIndex.h" -#include "itkImageRegionConstIterator.h" -#include -#include -#include "mitkNumericConstants.h" -#include "mitkLogMacros.h" -#include - -namespace itk -{ - template< class TInputImage, class TLabelImage > - bool - ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > - ::GetMaskingNonEmpty() const - { - return m_MaskNonEmpty; - } - - template< typename TInputImage, typename TLabelImage > - void - ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > - ::SetHistogramParameters(const int numBins, RealType lowerBound, RealType upperBound) - { - m_NumBins[0] = numBins; - m_LowerBound = lowerBound; - m_UpperBound = upperBound; - m_GlobalHistogramParametersSet = true; - m_PreferGlobalHistogramParameters = true; - this->Modified(); - } - - template< typename TInputImage, typename TLabelImage > - void - ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > - ::SetHistogramParametersForLabels(std::map numBins, std::map lowerBound, - std::map upperBound) - { - m_LabelMin = lowerBound; - m_LabelMax = upperBound; - m_LabelNBins = numBins; - m_LabelHistogramParametersSet = true; - m_PreferGlobalHistogramParameters = false; - this->Modified(); - } - - template< class TInputImage, class TLabelImage > - typename ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >::RealType - ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > - ::GetUniformity(LabelPixelType label) const - { - StatisticsMapConstIterator mapIt; - - mapIt = m_LabelStatistics.find(label); - if ( mapIt == m_LabelStatistics.end() ) - { - mitkThrow() << "Label does not exist"; - } - else - { - return ( *mapIt ).second.m_Uniformity; - } - } - - - template< class TInputImage, class TLabelImage > - typename ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >::RealType - ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > - ::GetMedian(LabelPixelType label) const - { - StatisticsMapConstIterator mapIt; - - mapIt = m_LabelStatistics.find(label); - if ( mapIt == m_LabelStatistics.end() ) - { - mitkThrow() << "Label does not exist"; - } - else - { - return ( *mapIt ).second.m_Median; - } - } - - template< class TInputImage, class TLabelImage > - typename ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >::RealType - ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > - ::GetEntropy(LabelPixelType label) const - { - StatisticsMapConstIterator mapIt; - - mapIt = m_LabelStatistics.find(label); - if ( mapIt == m_LabelStatistics.end() ) - { - mitkThrow() << "Label does not exist"; - } - else - { - return ( *mapIt ).second.m_Entropy; - } - } - - template< class TInputImage, class TLabelImage > - typename ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >::RealType - ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > - ::GetUPP(LabelPixelType label) const - { - StatisticsMapConstIterator mapIt; - - mapIt = m_LabelStatistics.find(label); - if ( mapIt == m_LabelStatistics.end() ) - { - mitkThrow() << "Label does not exist"; - } - else - { - return ( *mapIt ).second.m_UPP; - } - } - - template< class TInputImage, class TLabelImage > - typename ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >::RealType - ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > - ::GetMPP(LabelPixelType label) const - { - StatisticsMapConstIterator mapIt; - - mapIt = m_LabelStatistics.find(label); - if ( mapIt == m_LabelStatistics.end() ) - { - mitkThrow() << "Label does not exist"; - } - else - { - return ( *mapIt ).second.m_MPP; - } - } - - template< class TInputImage, class TLabelImage > - typename ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >::RealType - ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > - ::GetKurtosis(LabelPixelType label) const - { - StatisticsMapConstIterator mapIt; - - mapIt = m_LabelStatistics.find(label); - if ( mapIt == m_LabelStatistics.end() ) - { - mitkThrow() << "Label does not exist"; - } - else - { - return ( *mapIt ).second.m_Kurtosis; - } - } - - - template< class TInputImage, class TLabelImage > - typename ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >::RealType - ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > - ::GetSkewness(LabelPixelType label) const - { - StatisticsMapConstIterator mapIt; - - mapIt = m_LabelStatistics.find(label); - if ( mapIt == m_LabelStatistics.end() ) - { - mitkThrow() << "Label does not exist"; - } - else - { - return ( *mapIt ).second.m_Skewness; - } - } - - - - template< typename TInputImage, typename TLabelImage > - typename ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >::RealType - ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > - ::GetMinimum(LabelPixelType label) const - { - StatisticsMapConstIterator mapIt; - - mapIt = m_LabelStatistics.find(label); - if ( mapIt == m_LabelStatistics.end() ) - { - mitkThrow() << "Label does not exist"; - } - else - { - return ( *mapIt ).second.m_Minimum; - } - } - - template< typename TInputImage, typename TLabelImage > - typename ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >::RealType - ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > - ::GetMaximum(LabelPixelType label) const - { - StatisticsMapConstIterator mapIt; - - mapIt = m_LabelStatistics.find(label); - if ( mapIt == m_LabelStatistics.end() ) - { - mitkThrow() << "Label does not exist"; - } - else - { - return ( *mapIt ).second.m_Maximum; - } - } - - template< typename TInputImage, typename TLabelImage > - typename ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >::RealType - ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > - ::GetMean(LabelPixelType label) const - { - StatisticsMapConstIterator mapIt; - - mapIt = m_LabelStatistics.find(label); - if ( mapIt == m_LabelStatistics.end() ) - { - mitkThrow() << "Label does not exist"; - } - else - { - return ( *mapIt ).second.m_Mean; - } - } - - template< typename TInputImage, typename TLabelImage > - typename ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >::RealType - ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > - ::GetSum(LabelPixelType label) const - { - StatisticsMapConstIterator mapIt; - - mapIt = m_LabelStatistics.find(label); - if ( mapIt == m_LabelStatistics.end() ) - { - mitkThrow() << "Label does not exist"; - } - else - { - return ( *mapIt ).second.m_Sum; - } - } - - template< typename TInputImage, typename TLabelImage > - typename ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >::RealType - ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > - ::GetSigma(LabelPixelType label) const - { - StatisticsMapConstIterator mapIt; - - mapIt = m_LabelStatistics.find(label); - if ( mapIt == m_LabelStatistics.end() ) - { - mitkThrow() << "Label does not exist"; - } - else - { - return ( *mapIt ).second.m_Sigma; - } - } - - template< typename TInputImage, typename TLabelImage > - typename ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >::RealType - ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > - ::GetVariance(LabelPixelType label) const - { - StatisticsMapConstIterator mapIt; - - mapIt = m_LabelStatistics.find(label); - if ( mapIt == m_LabelStatistics.end() ) - { - mitkThrow() << "Label does not exist"; - } - else - { - return ( *mapIt ).second.m_Variance; - } - } - - template< typename TInputImage, typename TLabelImage > - typename ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >::BoundingBoxType - ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > - ::GetBoundingBox(LabelPixelType label) const - { - StatisticsMapConstIterator mapIt; - - mapIt = m_LabelStatistics.find(label); - if ( mapIt == m_LabelStatistics.end() ) - { - mitkThrow() << "Label does not exist"; - } - else - { - return ( *mapIt ).second.m_BoundingBox; - } - } - - template< typename TInputImage, typename TLabelImage > - typename ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >::RegionType - ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > - ::GetRegion(LabelPixelType label) const - { - StatisticsMapConstIterator mapIt; - - mapIt = m_LabelStatistics.find(label); - - if ( mapIt == m_LabelStatistics.end() ) - { - mitkThrow() << "Label does not exist"; - } - else - { - typename Superclass::BoundingBoxType bbox = this->GetBoundingBox(label); - typename Superclass::IndexType index; - typename Superclass::SizeType size; - - unsigned int dimension = bbox.size() / 2; - - for ( unsigned int i = 0; i < dimension; i++ ) - { - index[i] = bbox[2 * i]; - size[i] = bbox[2 * i + 1] - bbox[2 * i] + 1; - } - typename Superclass::RegionType region; - region.SetSize(size); - region.SetIndex(index); - - return region; - } - } - - template< typename TInputImage, typename TLabelImage > - typename ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >::MapSizeType - ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > - ::GetCount(LabelPixelType label) const - { - StatisticsMapConstIterator mapIt; - - mapIt = m_LabelStatistics.find(label); - if ( mapIt == m_LabelStatistics.end() ) - { - mitkThrow() << "Label does not exist"; - } - else - { - return ( *mapIt ).second.m_Count; - } - } - - - template< typename TInputImage, typename TLabelImage > - typename ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >::HistogramType::Pointer - ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > - ::GetHistogram(LabelPixelType label) const - { - StatisticsMapConstIterator mapIt; - - mapIt = m_LabelStatistics.find(label); - if ( mapIt == m_LabelStatistics.end() ) - { - mitkThrow() << "Label does not exist"; - } - else - { - // this will be zero if histograms have not been enabled - return ( *mapIt ).second.m_Histogram; - } - } - - - - template< typename TInputImage, typename TLabelImage > - void - ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > - ::BeforeStreamedGenerateData() - { - ThreadIdType numberOfThreads = this->GetNumberOfThreads(); - - // Resize the thread temporaries - m_LabelStatisticsPerThread.resize(numberOfThreads); - - // Initialize the temporaries - for ( ThreadIdType i = 0; i < numberOfThreads; ++i ) - { - m_LabelStatisticsPerThread[i].clear(); - } - - // Initialize the final map - m_LabelStatistics.clear(); - } - - template< typename TInputImage, typename TLabelImage > - std::list - ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > - ::GetRelevantLabels() const - { - std::list< int> relevantLabels; - for (int i = 0; i < 4096; ++i ) - { - if ( this->HasLabel( i ) ) - { - relevantLabels.push_back( i ); - } - } - return relevantLabels; - } - - template - void ExtendedLabelStatisticsImageFilter::ThreadedStreamedGenerateData(const RegionType& region) - { - - typename HistogramType::IndexType histogramIndex(1); - typename HistogramType::MeasurementVectorType histogramMeasurement(1); - - const SizeValueType size0 = outputRegionForThread.GetSize(0); - if( size0 == 0) - { - return; - } - - ImageLinearConstIteratorWithIndex< TInputImage > it (this->GetInput(), - outputRegionForThread); - - ImageScanlineConstIterator< TLabelImage > labelIt (this->GetLabelInput(), - outputRegionForThread); - - StatisticsMapIterator mapIt; - - // support progress methods/callbacks - const size_t numberOfLinesToProcess = outputRegionForThread.GetNumberOfPixels() / size0; - ProgressReporter progress( this, threadId, numberOfLinesToProcess ); - - typedef typename MapType::value_type MapValueType; - - // do the work - while ( !it.IsAtEnd() ) - { - while ( !it.IsAtEndOfLine() ) - { - const RealType & value = static_cast< RealType >( it.Get() ); - - const LabelPixelType & label = labelIt.Get(); - - // is the label already in this thread? - mapIt = m_LabelStatisticsPerThread[threadId].find(label); - if ( mapIt == m_LabelStatisticsPerThread[threadId].end() ) - { - // if global histogram parameters are set and preferred then use them - if ( m_PreferGlobalHistogramParameters && m_GlobalHistogramParametersSet ) - { - mapIt = m_LabelStatisticsPerThread[threadId].insert( MapValueType( label, - LabelStatistics(m_NumBins[0], m_LowerBound, - m_UpperBound) ) ).first; - } - // if we have label histogram parameters then use them. If we encounter a label that has no parameters then use global settings if available - else if(!m_PreferGlobalHistogramParameters && m_LabelHistogramParametersSet) - { - typename std::map::iterator lbIt, ubIt; - typename std::map::iterator nbIt; - - lbIt = m_LabelMin.find(label); - ubIt = m_LabelMax.find(label); - nbIt = m_LabelNBins.find(label); - - // if any of the parameters is lacking for the current label but global histogram params are available, use the global parameters - if ((lbIt == m_LabelMin.end() || ubIt == m_LabelMax.end() || nbIt == m_LabelNBins.end()) && m_GlobalHistogramParametersSet) - { - mapIt = m_LabelStatisticsPerThread[threadId].insert( MapValueType( label, - LabelStatistics(m_NumBins[0], m_LowerBound, - m_UpperBound) ) ).first; - } - // if any of the parameters is lacking for the current label and global histogram params are not available, dont use histograms for this label - else if ((lbIt == m_LabelMin.end() || ubIt == m_LabelMax.end() || nbIt == m_LabelNBins.end()) && !m_GlobalHistogramParametersSet) - { - mapIt = m_LabelStatisticsPerThread[threadId].insert( MapValueType( label, - LabelStatistics() ) ).first; - } - // label histogram parameters are available, use them! - else - { - PixelType lowerBound, upperBound; - unsigned int nBins; - lowerBound = (*lbIt).second; - upperBound = (*ubIt).second; - nBins = (*nbIt).second; - mapIt = m_LabelStatisticsPerThread[threadId].insert( MapValueType( label, - LabelStatistics(nBins, lowerBound, upperBound) ) ).first; - } - } - // neither global nor label specific histogram parameters are set -> don't use histograms - else - { - mapIt = m_LabelStatisticsPerThread[threadId].insert( MapValueType( label, - LabelStatistics() ) ).first; - } - } - - typename MapType::mapped_type &labelStats = ( *mapIt ).second; - - // update the values for this label and this thread - if ( value < labelStats.m_Minimum ) - { - labelStats.m_Minimum = value; - } - if ( value > labelStats.m_Maximum ) - { - labelStats.m_Maximum = value; - } - - // bounding box is min,max pairs - for ( unsigned int i = 0; i < ( 2 * TInputImage::ImageDimension ); i += 2 ) - { - const typename TInputImage::IndexType & index = it.GetIndex(); - if ( labelStats.m_BoundingBox[i] > index[i / 2] ) - { - labelStats.m_BoundingBox[i] = index[i / 2]; - } - if ( labelStats.m_BoundingBox[i + 1] < index[i / 2] ) - { - labelStats.m_BoundingBox[i + 1] = index[i / 2]; - } - } - - labelStats.m_Sum += value; - labelStats.m_SumOfSquares += ( value * value ); - labelStats.m_Count++; - labelStats.m_SumOfCubes += std::pow(value, 3.); - labelStats.m_SumOfQuadruples += std::pow(value, 4.); - - if (value > 0) - { - labelStats.m_PositivePixelCount++; - labelStats.m_SumOfPositivePixels += value; - } - - // if enabled, update the histogram for this label - if ( labelStats.m_Histogram.IsNotNull() ) - { - histogramMeasurement[0] = value; - labelStats.m_Histogram->GetIndex(histogramMeasurement, histogramIndex); - labelStats.m_Histogram->IncreaseFrequencyOfIndex(histogramIndex, 1); - } - - ++labelIt; - ++it; - } - labelIt.NextLine(); - it.NextLine(); - progress.CompletedPixel(); - } - - } - - - template< class TInputImage, class TLabelImage > - void ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >:: - AfterStreamedGenerateData() - { - StatisticsMapIterator mapIt; - StatisticsMapConstIterator threadIt; - ThreadIdType i; - ThreadIdType numberOfThreads = this->GetNumberOfThreads(); - - // Run through the map for each thread and accumulate the count, - // sum, and sumofsquares - for ( i = 0; i < numberOfThreads; i++ ) - { - // iterate over the map for this thread - for ( threadIt = m_LabelStatisticsPerThread[i].begin(); - threadIt != m_LabelStatisticsPerThread[i].end(); - ++threadIt ) - { - // does this label exist in the cumulative structure yet? - mapIt = m_LabelStatistics.find( ( *threadIt ).first ); - if ( mapIt == m_LabelStatistics.end() ) - { - // create a new entry - typedef typename MapType::value_type MapValueType; - if ( m_GlobalHistogramParametersSet || m_LabelHistogramParametersSet ) - { -// mapIt = m_LabelStatistics.insert( MapValueType( ( *threadIt ).first, -// LabelStatistics(m_NumBins[0], m_LowerBound, -// m_UpperBound) ) ).first; - mapIt = m_LabelStatistics.insert( MapValueType( *threadIt ) ).first; - continue; - } - else - { - mapIt = m_LabelStatistics.insert( MapValueType( ( *threadIt ).first, - LabelStatistics() ) ).first; - } - } - - - typename MapType::mapped_type &labelStats = ( *mapIt ).second; - - // accumulate the information from this thread - labelStats.m_Count += ( *threadIt ).second.m_Count; - labelStats.m_Sum += ( *threadIt ).second.m_Sum; - labelStats.m_SumOfSquares += ( *threadIt ).second.m_SumOfSquares; - labelStats.m_SumOfPositivePixels += ( *threadIt ).second.m_SumOfPositivePixels; - labelStats.m_PositivePixelCount += ( *threadIt ).second.m_PositivePixelCount; - labelStats.m_SumOfCubes += ( *threadIt ).second.m_SumOfCubes; - labelStats.m_SumOfQuadruples += ( *threadIt ).second.m_SumOfQuadruples; - - if ( labelStats.m_Minimum > ( *threadIt ).second.m_Minimum ) - { - labelStats.m_Minimum = ( *threadIt ).second.m_Minimum; - } - if ( labelStats.m_Maximum < ( *threadIt ).second.m_Maximum ) - { - labelStats.m_Maximum = ( *threadIt ).second.m_Maximum; - } - - //bounding box is min,max pairs - int dimension = labelStats.m_BoundingBox.size() / 2; - for ( int ii = 0; ii < ( dimension * 2 ); ii += 2 ) - { - if ( labelStats.m_BoundingBox[ii] > ( *threadIt ).second.m_BoundingBox[ii] ) - { - labelStats.m_BoundingBox[ii] = ( *threadIt ).second.m_BoundingBox[ii]; - } - if ( labelStats.m_BoundingBox[ii + 1] < ( *threadIt ).second.m_BoundingBox[ii + 1] ) - { - labelStats.m_BoundingBox[ii + 1] = ( *threadIt ).second.m_BoundingBox[ii + 1]; - } - } - - // if enabled, update the histogram for this label - if ( m_GlobalHistogramParametersSet || m_LabelHistogramParametersSet ) - { - typename HistogramType::IndexType index; - index.SetSize(1); - for ( unsigned int bin = 0; bin < labelStats.m_Histogram->Size(); bin++ ) - { - index[0] = bin; - labelStats.m_Histogram->IncreaseFrequency( bin, ( *threadIt ).second.m_Histogram->GetFrequency(bin) ); - } - } - } // end of thread map iterator loop - } // end of thread loop - - // compute the remainder of the statistics - for ( mapIt = m_LabelStatistics.begin(); - mapIt != m_LabelStatistics.end(); - ++mapIt ) - { - typename MapType::mapped_type &labelStats = ( *mapIt ).second; - - // mean - labelStats.m_Mean = labelStats.m_Sum - / static_cast< RealType >( labelStats.m_Count ); - - // MPP - labelStats.m_MPP = labelStats.m_SumOfPositivePixels - / static_cast< RealType >( labelStats.m_PositivePixelCount ); - - // variance - if ( labelStats.m_Count > 0 ) - { - // unbiased estimate of variance - LabelStatistics & ls = mapIt->second; - const RealType sumSquared = ls.m_Sum * ls.m_Sum; - const RealType count = static_cast< RealType >( ls.m_Count ); - - ls.m_Variance = ( ls.m_SumOfSquares - sumSquared / count ) / ( count ); - - RealType secondMoment = ls.m_SumOfSquares / count; - RealType thirdMoment = ls.m_SumOfCubes / count; - RealType fourthMoment = ls.m_SumOfQuadruples / count; - - ls.m_Skewness = (thirdMoment - 3. * secondMoment * ls.m_Mean + 2. * std::pow(ls.m_Mean, 3.)) / std::pow(secondMoment - std::pow(ls.m_Mean, 2.), 1.5); // see http://www.boost.org/doc/libs/1_51_0/doc/html/boost/accumulators/impl/skewness_impl.html - ls.m_Kurtosis = (fourthMoment - 4. * thirdMoment * ls.m_Mean + 6. * secondMoment * std::pow(ls.m_Mean, 2.) - 3. * std::pow(ls.m_Mean, 4.)) / std::pow(secondMoment - std::pow(ls.m_Mean, 2.), 2.); // see http://www.boost.org/doc/libs/1_51_0/doc/html/boost/accumulators/impl/kurtosis_impl.html, dropped -3 - } - else - { - labelStats.m_Variance = NumericTraits< RealType >::ZeroValue(); - labelStats.m_Skewness = NumericTraits< RealType >::ZeroValue(); - labelStats.m_Kurtosis = NumericTraits< RealType >::ZeroValue(); - } - - // sigma - labelStats.m_Sigma = std::sqrt( labelStats.m_Variance ); - - // histogram statistics - if (labelStats.m_Histogram.IsNotNull()) - { - mitk::HistogramStatisticsCalculator histStatCalc; - histStatCalc.SetHistogram(labelStats.m_Histogram); - histStatCalc.CalculateStatistics(); - labelStats.m_Median = histStatCalc.GetMedian(); - labelStats.m_Entropy = histStatCalc.GetEntropy(); - labelStats.m_Uniformity = histStatCalc.GetUniformity(); - labelStats.m_UPP = histStatCalc.GetUPP(); - } - - } - - { - //Now update the cached vector of valid labels. - m_ValidLabelValues.resize(0); - m_ValidLabelValues.reserve(m_LabelStatistics.size()); - for ( mapIt = m_LabelStatistics.begin(); - mapIt != m_LabelStatistics.end(); - ++mapIt ) - { - m_ValidLabelValues.push_back(mapIt->first); - } - } - } - -} // end namespace itk - -#endif diff --git a/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp b/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp index 65c9962af9..c35e9d86db 100644 --- a/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp +++ b/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp @@ -1,552 +1,549 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #include "mitkImageStatisticsCalculator.h" -#include +#include #include #include #include #include #include #include #include #include #include #include #include namespace mitk { void ImageStatisticsCalculator::SetInputImage(const mitk::Image *image) { if (image != m_Image) { m_Image = image; this->Modified(); } } void ImageStatisticsCalculator::SetMask(mitk::MaskGenerator *mask) { if (mask != m_MaskGenerator) { m_MaskGenerator = mask; this->Modified(); } } void ImageStatisticsCalculator::SetSecondaryMask(mitk::MaskGenerator *mask) { if (mask != m_SecondaryMaskGenerator) { m_SecondaryMaskGenerator = mask; this->Modified(); } } void ImageStatisticsCalculator::SetNBinsForHistogramStatistics(unsigned int nBins) { if (nBins != m_nBinsForHistogramStatistics) { m_nBinsForHistogramStatistics = nBins; this->Modified(); this->m_UseBinSizeOverNBins = false; } if (m_UseBinSizeOverNBins) { this->Modified(); this->m_UseBinSizeOverNBins = false; } } unsigned int ImageStatisticsCalculator::GetNBinsForHistogramStatistics() const { return m_nBinsForHistogramStatistics; } void ImageStatisticsCalculator::SetBinSizeForHistogramStatistics(double binSize) { if (binSize != m_binSizeForHistogramStatistics) { m_binSizeForHistogramStatistics = binSize; this->Modified(); this->m_UseBinSizeOverNBins = true; } if (!m_UseBinSizeOverNBins) { this->Modified(); this->m_UseBinSizeOverNBins = true; } } double ImageStatisticsCalculator::GetBinSizeForHistogramStatistics() const { return m_binSizeForHistogramStatistics; } mitk::ImageStatisticsContainer* ImageStatisticsCalculator::GetStatistics(LabelIndex label) { if (m_Image.IsNull()) { mitkThrow() << "no image"; } if (!m_Image->IsInitialized()) { mitkThrow() << "Image not initialized!"; } if (IsUpdateRequired(label)) { auto timeGeometry = m_Image->GetTimeGeometry(); // always compute statistics on all timesteps for (unsigned int timeStep = 0; timeStep < m_Image->GetTimeSteps(); timeStep++) { if (m_MaskGenerator.IsNotNull()) { m_MaskGenerator->SetTimeStep(timeStep); //See T25625: otherwise, the mask is not computed again after setting a different time step m_MaskGenerator->Modified(); m_InternalMask = m_MaskGenerator->GetMask(); if (m_MaskGenerator->GetReferenceImage().IsNotNull()) { m_InternalImageForStatistics = m_MaskGenerator->GetReferenceImage(); } else { m_InternalImageForStatistics = m_Image; } } else { m_InternalImageForStatistics = m_Image; } if (m_SecondaryMaskGenerator.IsNotNull()) { m_SecondaryMaskGenerator->SetTimeStep(timeStep); m_SecondaryMask = m_SecondaryMaskGenerator->GetMask(); } ImageTimeSelector::Pointer imgTimeSel = ImageTimeSelector::New(); imgTimeSel->SetInput(m_InternalImageForStatistics); imgTimeSel->SetTimeNr(timeStep); imgTimeSel->UpdateLargestPossibleRegion(); imgTimeSel->Update(); m_ImageTimeSlice = imgTimeSel->GetOutput(); // Calculate statistics with/without mask if (m_MaskGenerator.IsNull() && m_SecondaryMaskGenerator.IsNull()) { // 1) calculate statistics unmasked: AccessByItk_2(m_ImageTimeSlice, InternalCalculateStatisticsUnmasked, timeGeometry, timeStep) } else { // 2) calculate statistics masked AccessByItk_2(m_ImageTimeSlice, InternalCalculateStatisticsMasked, timeGeometry, timeStep) } } } auto it = m_StatisticContainers.find(label); if (it != m_StatisticContainers.end()) { return (it->second).GetPointer(); } else { mitkThrow() << "unknown label"; return nullptr; } } template void ImageStatisticsCalculator::InternalCalculateStatisticsUnmasked( typename itk::Image *image, const TimeGeometry *timeGeometry, TimeStepType timeStep) { typedef typename itk::Image ImageType; typedef typename StatisticsImageFilter ImageStatisticsFilterType; typedef typename itk::MinMaxImageFilterWithIndex MinMaxFilterType; // reset statistics container if exists ImageStatisticsContainer::Pointer statisticContainerForImage; LabelIndex labelNoMask = 1; auto it = m_StatisticContainers.find(labelNoMask); if (it != m_StatisticContainers.end()) { statisticContainerForImage = it->second; } else { statisticContainerForImage = ImageStatisticsContainer::New(); statisticContainerForImage->SetTimeGeometry(const_cast(timeGeometry)); m_StatisticContainers.emplace(labelNoMask, statisticContainerForImage); } auto statObj = ImageStatisticsContainer::ImageStatisticsObject(); typename ImageStatisticsFilterType::Pointer statisticsFilter = ImageStatisticsFilterType::New(); statisticsFilter->SetInput(image); statisticsFilter->SetCoordinateTolerance(0.001); statisticsFilter->SetDirectionTolerance(0.001); // TODO: this is single threaded. Implement our own image filter that does this multi threaded // typename itk::MinimumMaximumImageCalculator::Pointer imgMinMaxFilter = // itk::MinimumMaximumImageCalculator::New(); imgMinMaxFilter->SetImage(image); // imgMinMaxFilter->Compute(); vnl_vector minIndex, maxIndex; typename MinMaxFilterType::Pointer minMaxFilter = MinMaxFilterType::New(); minMaxFilter->SetInput(image); minMaxFilter->UpdateLargestPossibleRegion(); typename ImageType::PixelType minval = minMaxFilter->GetMin(); typename ImageType::PixelType maxval = minMaxFilter->GetMax(); typename ImageType::IndexType tmpMinIndex = minMaxFilter->GetMinIndex(); typename ImageType::IndexType tmpMaxIndex = minMaxFilter->GetMaxIndex(); // typename ImageType::IndexType tmpMinIndex = imgMinMaxFilter->GetIndexOfMinimum(); // typename ImageType::IndexType tmpMaxIndex = imgMinMaxFilter->GetIndexOfMaximum(); minIndex.set_size(tmpMaxIndex.GetIndexDimension()); maxIndex.set_size(tmpMaxIndex.GetIndexDimension()); for (unsigned int i = 0; i < tmpMaxIndex.GetIndexDimension(); i++) { minIndex[i] = tmpMinIndex[i]; maxIndex[i] = tmpMaxIndex[i]; } statObj.AddStatistic(mitk::ImageStatisticsConstants::MINIMUMPOSITION(), minIndex); statObj.AddStatistic(mitk::ImageStatisticsConstants::MAXIMUMPOSITION(), maxIndex); // convert m_binSize in m_nBins if necessary unsigned int nBinsForHistogram; if (m_UseBinSizeOverNBins) { nBinsForHistogram = std::max(static_cast(std::ceil(maxval - minval)) / m_binSizeForHistogramStatistics, 10.); // do not allow less than 10 bins } else { nBinsForHistogram = m_nBinsForHistogramStatistics; } statisticsFilter->SetHistogramParameters(nBinsForHistogram, minval, maxval); try { statisticsFilter->Update(); } catch (const itk::ExceptionObject &e) { mitkThrow() << "Image statistics calculation failed due to following ITK Exception: \n " << e.what(); } auto voxelVolume = GetVoxelVolume(image); auto numberOfPixels = image->GetLargestPossibleRegion().GetNumberOfPixels(); auto volume = static_cast(numberOfPixels) * voxelVolume; auto variance = statisticsFilter->GetSigma() * statisticsFilter->GetSigma(); auto rms = std::sqrt(std::pow(statisticsFilter->GetMean(), 2.) + statisticsFilter->GetVariance()); // variance = sigma^2 statObj.AddStatistic(mitk::ImageStatisticsConstants::NUMBEROFVOXELS(), static_cast(numberOfPixels)); statObj.AddStatistic(mitk::ImageStatisticsConstants::VOLUME(), volume); statObj.AddStatistic(mitk::ImageStatisticsConstants::MEAN(), statisticsFilter->GetMean()); statObj.AddStatistic(mitk::ImageStatisticsConstants::MINIMUM(), static_cast(statisticsFilter->GetMinimum())); statObj.AddStatistic(mitk::ImageStatisticsConstants::MAXIMUM(), static_cast(statisticsFilter->GetMaximum())); statObj.AddStatistic(mitk::ImageStatisticsConstants::STANDARDDEVIATION(), statisticsFilter->GetSigma()); statObj.AddStatistic(mitk::ImageStatisticsConstants::VARIANCE(), variance); statObj.AddStatistic(mitk::ImageStatisticsConstants::SKEWNESS(), statisticsFilter->GetSkewness()); statObj.AddStatistic(mitk::ImageStatisticsConstants::KURTOSIS(), statisticsFilter->GetKurtosis()); statObj.AddStatistic(mitk::ImageStatisticsConstants::RMS(), rms); statObj.AddStatistic(mitk::ImageStatisticsConstants::MPP(), statisticsFilter->GetMPP()); statObj.AddStatistic(mitk::ImageStatisticsConstants::ENTROPY(), statisticsFilter->GetEntropy()); statObj.AddStatistic(mitk::ImageStatisticsConstants::MEDIAN(), statisticsFilter->GetMedian()); statObj.AddStatistic(mitk::ImageStatisticsConstants::UNIFORMITY(), statisticsFilter->GetUniformity()); statObj.AddStatistic(mitk::ImageStatisticsConstants::UPP(), statisticsFilter->GetUPP()); statObj.m_Histogram = statisticsFilter->GetHistogram(); statisticContainerForImage->SetStatisticsForTimeStep(timeStep, statObj); } template double ImageStatisticsCalculator::GetVoxelVolume(typename itk::Image *image) const { auto spacing = image->GetSpacing(); double voxelVolume = 1.; for (unsigned int i = 0; i < image->GetImageDimension(); i++) { voxelVolume *= spacing[i]; } return voxelVolume; } template void ImageStatisticsCalculator::InternalCalculateStatisticsMasked(typename itk::Image *image, const TimeGeometry *timeGeometry, unsigned int timeStep) { typedef itk::Image ImageType; typedef itk::Image MaskType; typedef typename MaskType::PixelType LabelPixelType; - typedef itk::ExtendedLabelStatisticsImageFilter ImageStatisticsFilterType; + typedef LabelStatisticsImageFilter ImageStatisticsFilterType; typedef MaskUtilities MaskUtilType; typedef typename itk::MinMaxLabelImageFilterWithIndex MinMaxLabelFilterType; typedef typename ImageType::PixelType InputImgPixelType; // workaround: if m_SecondaryMaskGenerator ist not null but m_MaskGenerator is! (this is the case if we request a // 'ignore zuero valued pixels' mask in the gui but do not define a primary mask) bool swapMasks = false; if (m_SecondaryMask.IsNotNull() && m_InternalMask.IsNull()) { m_InternalMask = m_SecondaryMask; m_SecondaryMask = nullptr; swapMasks = true; } // maskImage has to have the same dimension as image typename MaskType::Pointer maskImage = MaskType::New(); try { // try to access the pixel values directly (no copying or casting). Only works if mask pixels are of pixelType // unsigned short maskImage = ImageToItkImage(m_InternalMask); } catch (const itk::ExceptionObject &) { // if the pixel type of the mask is not short, then we have to make a copy of m_InternalMask (and cast the values) CastToItkImage(m_InternalMask, maskImage); } // if we have a secondary mask (say a ignoreZeroPixelMask) we need to combine the masks (corresponds to AND) if (m_SecondaryMask.IsNotNull()) { // dirty workaround for a bug when pf mask + any other mask is used in conjunction. We need a proper fix for this // (Fabian Isensee is responsible and probably working on it!) if (m_InternalMask->GetDimension() == 2 && (m_SecondaryMask->GetDimension() == 3 || m_SecondaryMask->GetDimension() == 4)) { mitk::Image::ConstPointer old_img = m_SecondaryMaskGenerator->GetReferenceImage(); m_SecondaryMaskGenerator->SetInputImage(m_MaskGenerator->GetReferenceImage()); m_SecondaryMask = m_SecondaryMaskGenerator->GetMask(); m_SecondaryMaskGenerator->SetInputImage(old_img); } typename MaskType::Pointer secondaryMaskImage = MaskType::New(); secondaryMaskImage = ImageToItkImage(m_SecondaryMask); // secondary mask should be a ignore zero value pixel mask derived from image. it has to be cropped to the mask // region (which may be planar or simply smaller) typename MaskUtilities::Pointer secondaryMaskMaskUtil = MaskUtilities::New(); secondaryMaskMaskUtil->SetImage(secondaryMaskImage.GetPointer()); secondaryMaskMaskUtil->SetMask(maskImage.GetPointer()); typename MaskType::Pointer adaptedSecondaryMaskImage = secondaryMaskMaskUtil->ExtractMaskImageRegion(); typename itk::MaskImageFilter2::Pointer maskFilter = itk::MaskImageFilter2::New(); maskFilter->SetInput1(maskImage); maskFilter->SetInput2(adaptedSecondaryMaskImage); maskFilter->SetMaskingValue( 1); // all pixels of maskImage where secondaryMaskImage==1 will be kept, all the others are set to 0 maskFilter->UpdateLargestPossibleRegion(); maskImage = maskFilter->GetOutput(); } 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(); adaptedImage = maskUtil->ExtractMaskImageRegion(); // this also checks mask sanity // find min, max, minindex and maxindex typename MinMaxLabelFilterType::Pointer minMaxFilter = MinMaxLabelFilterType::New(); minMaxFilter->SetInput(adaptedImage); minMaxFilter->SetLabelInput(maskImage); minMaxFilter->UpdateLargestPossibleRegion(); // set histogram parameters for each label individually (min/max may be different for each label) - typedef typename std::map MapType; - typedef typename std::pair PairType; + typedef typename std::unordered_map MapType; std::vector relevantLabels = minMaxFilter->GetRelevantLabels(); MapType minVals; MapType maxVals; - std::map nBins; + std::unordered_map nBins; for (LabelPixelType label : relevantLabels) { - minVals.insert(PairType(label, minMaxFilter->GetMin(label))); - maxVals.insert(PairType(label, minMaxFilter->GetMax(label))); + minVals[label] = static_cast(minMaxFilter->GetMin(label)); + maxVals[label] = static_cast(minMaxFilter->GetMax(label)); unsigned int nBinsForHistogram; if (m_UseBinSizeOverNBins) { nBinsForHistogram = std::max(static_cast(std::ceil(minMaxFilter->GetMax(label) - minMaxFilter->GetMin(label))) / m_binSizeForHistogramStatistics, 10.); // do not allow less than 10 bins } else { nBinsForHistogram = m_nBinsForHistogramStatistics; } - nBins.insert(typename std::pair(label, nBinsForHistogram)); + nBins[label] = nBinsForHistogram; } typename ImageStatisticsFilterType::Pointer imageStatisticsFilter = ImageStatisticsFilterType::New(); imageStatisticsFilter->SetDirectionTolerance(0.001); imageStatisticsFilter->SetCoordinateTolerance(0.001); imageStatisticsFilter->SetInput(adaptedImage); imageStatisticsFilter->SetLabelInput(maskImage); - imageStatisticsFilter->SetHistogramParametersForLabels(nBins, minVals, maxVals); + imageStatisticsFilter->SetHistogramParameters(nBins, minVals, maxVals); imageStatisticsFilter->Update(); - std::list labels = imageStatisticsFilter->GetRelevantLabels(); + auto labels = imageStatisticsFilter->GetValidLabelValues(); auto it = labels.begin(); while (it != labels.end()) { ImageStatisticsContainer::Pointer statisticContainerForLabelImage; auto labelIt = m_StatisticContainers.find(*it); // reset if statisticContainer already exist if (labelIt != m_StatisticContainers.end()) { statisticContainerForLabelImage = labelIt->second; } // create new statisticContainer else { statisticContainerForLabelImage = ImageStatisticsContainer::New(); statisticContainerForLabelImage->SetTimeGeometry(const_cast(timeGeometry)); // link label (*it) to statisticContainer m_StatisticContainers.emplace(*it, statisticContainerForLabelImage); } ImageStatisticsContainer::ImageStatisticsObject statObj; // find min, max, minindex and maxindex // make sure to only look in the masked region, use a masker for this vnl_vector minIndex, maxIndex; mitk::Point3D worldCoordinateMin; mitk::Point3D worldCoordinateMax; mitk::Point3D indexCoordinateMin; mitk::Point3D indexCoordinateMax; m_InternalImageForStatistics->GetGeometry()->IndexToWorld(minMaxFilter->GetMinIndex(*it), worldCoordinateMin); m_InternalImageForStatistics->GetGeometry()->IndexToWorld(minMaxFilter->GetMaxIndex(*it), worldCoordinateMax); m_Image->GetGeometry()->WorldToIndex(worldCoordinateMin, indexCoordinateMin); m_Image->GetGeometry()->WorldToIndex(worldCoordinateMax, indexCoordinateMax); minIndex.set_size(3); maxIndex.set_size(3); // for (unsigned int i=0; i < tmpMaxIndex.GetIndexDimension(); i++) for (unsigned int i = 0; i < 3; i++) { minIndex[i] = indexCoordinateMin[i]; maxIndex[i] = indexCoordinateMax[i]; } statObj.AddStatistic(mitk::ImageStatisticsConstants::MINIMUMPOSITION(), minIndex); statObj.AddStatistic(mitk::ImageStatisticsConstants::MAXIMUMPOSITION(), maxIndex); - assert(std::abs(minMaxFilter->GetMax(*it) - imageStatisticsFilter->GetMaximum(*it)) < mitk::eps); - assert(std::abs(minMaxFilter->GetMin(*it) - imageStatisticsFilter->GetMinimum(*it)) < mitk::eps); - auto voxelVolume = GetVoxelVolume(image); auto numberOfVoxels = static_cast(imageStatisticsFilter->GetCount(*it)); auto volume = static_cast(numberOfVoxels) * voxelVolume; auto rms = std::sqrt(std::pow(imageStatisticsFilter->GetMean(*it), 2.) + imageStatisticsFilter->GetVariance(*it)); // variance = sigma^2 auto variance = imageStatisticsFilter->GetSigma(*it) * imageStatisticsFilter->GetSigma(*it); statObj.AddStatistic(mitk::ImageStatisticsConstants::NUMBEROFVOXELS(), numberOfVoxels); statObj.AddStatistic(mitk::ImageStatisticsConstants::VOLUME(), volume); statObj.AddStatistic(mitk::ImageStatisticsConstants::MEAN(), imageStatisticsFilter->GetMean(*it)); - statObj.AddStatistic(mitk::ImageStatisticsConstants::MINIMUM(), imageStatisticsFilter->GetMinimum(*it)); - statObj.AddStatistic(mitk::ImageStatisticsConstants::MAXIMUM(), imageStatisticsFilter->GetMaximum(*it)); + statObj.AddStatistic(mitk::ImageStatisticsConstants::MINIMUM(), + static_cast(imageStatisticsFilter->GetMinimum(*it))); + statObj.AddStatistic(mitk::ImageStatisticsConstants::MAXIMUM(), + static_cast(imageStatisticsFilter->GetMaximum(*it))); statObj.AddStatistic(mitk::ImageStatisticsConstants::STANDARDDEVIATION(), imageStatisticsFilter->GetSigma(*it)); statObj.AddStatistic(mitk::ImageStatisticsConstants::VARIANCE(), variance); statObj.AddStatistic(mitk::ImageStatisticsConstants::SKEWNESS(), imageStatisticsFilter->GetSkewness(*it)); statObj.AddStatistic(mitk::ImageStatisticsConstants::KURTOSIS(), imageStatisticsFilter->GetKurtosis(*it)); statObj.AddStatistic(mitk::ImageStatisticsConstants::RMS(), rms); statObj.AddStatistic(mitk::ImageStatisticsConstants::MPP(), imageStatisticsFilter->GetMPP(*it)); statObj.AddStatistic(mitk::ImageStatisticsConstants::ENTROPY(), imageStatisticsFilter->GetEntropy(*it)); statObj.AddStatistic(mitk::ImageStatisticsConstants::MEDIAN(), imageStatisticsFilter->GetMedian(*it)); statObj.AddStatistic(mitk::ImageStatisticsConstants::UNIFORMITY(), imageStatisticsFilter->GetUniformity(*it)); statObj.AddStatistic(mitk::ImageStatisticsConstants::UPP(), imageStatisticsFilter->GetUPP(*it)); - statObj.m_Histogram = imageStatisticsFilter->GetHistogram(*it).GetPointer(); - + statObj.m_Histogram = imageStatisticsFilter->GetHistogram(*it); statisticContainerForLabelImage->SetStatisticsForTimeStep(timeStep, statObj); ++it; } // swap maskGenerators back if (swapMasks) { m_SecondaryMask = m_InternalMask; m_InternalMask = nullptr; } } bool ImageStatisticsCalculator::IsUpdateRequired(LabelIndex label) const { unsigned long thisClassTimeStamp = this->GetMTime(); unsigned long inputImageTimeStamp = m_Image->GetMTime(); auto it = m_StatisticContainers.find(label); if (it == m_StatisticContainers.end()) { return true; } unsigned long statisticsTimeStamp = it->second->GetMTime(); if (thisClassTimeStamp > statisticsTimeStamp) // inputs have changed { return true; } if (inputImageTimeStamp > statisticsTimeStamp) // image has changed { return true; } if (m_MaskGenerator.IsNotNull()) { unsigned long maskGeneratorTimeStamp = m_MaskGenerator->GetMTime(); if (maskGeneratorTimeStamp > statisticsTimeStamp) // there is a mask generator and it has changed { return true; } } if (m_SecondaryMaskGenerator.IsNotNull()) { unsigned long maskGeneratorTimeStamp = m_SecondaryMaskGenerator->GetMTime(); if (maskGeneratorTimeStamp > statisticsTimeStamp) // there is a secondary mask generator and it has changed { return true; } } return false; } } // namespace mitk diff --git a/Modules/ImageStatistics/mitkLabelStatisticsImageFilter.h b/Modules/ImageStatistics/mitkLabelStatisticsImageFilter.h new file mode 100644 index 0000000000..fe9beb7000 --- /dev/null +++ b/Modules/ImageStatistics/mitkLabelStatisticsImageFilter.h @@ -0,0 +1,168 @@ +/*============================================================================ + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center (DKFZ) +All rights reserved. + +Use of this source code is governed by a 3-clause BSD license that can be +found in the LICENSE file. + +============================================================================*/ + +// This file is based on ITK's itkLabelStatisticsImageFilter.h + +#ifndef mitkLabelStatisticsImageFilter_h +#define mitkLabelStatisticsImageFilter_h + +#include +#include +#include +#include +#include + +#include +#include +#include + +#include + +namespace mitk +{ + template + class LabelStatisticsImageFilter : public itk::ImageSink + { + public: + using Self = LabelStatisticsImageFilter; + using Superclass = itk::ImageSink; + using Pointer = itk::SmartPointer; + using ConstPointer = itk::SmartPointer; + + itkFactorylessNewMacro(Self); + + itkTypeMacro(LabelStatisticsImageFilter, itk::ImageSink); + + using IndexType = typename TInputImage::IndexType; + using SizeType = typename TInputImage::SizeType; + using RegionType = typename TInputImage::RegionType; + using PixelType = typename TInputImage::PixelType; + using LabelPixelType = typename mitk::Label::PixelType; + + static constexpr unsigned int ImageDimension = TInputImage::ImageDimension; + + using RealType = typename itk::NumericTraits::RealType; + + using DataObjectPointer = typename itk::DataObject::Pointer; + + using RealObjectType = itk::SimpleDataObjectDecorator; + + using BoundingBoxType = std::vector; + + using HistogramType = itk::Statistics::Histogram; + using HistogramPointer = typename HistogramType::Pointer; + + class LabelStatistics + { + public: + LabelStatistics(); + LabelStatistics(unsigned int size, RealType lowerBound, RealType upperBound); + ~LabelStatistics(); + + itk::SizeValueType m_Count; + itk::SizeValueType m_CountOfPositivePixels; + RealType m_Min; + RealType m_Max; + RealType m_Mean; + itk::CompensatedSummation m_Sum; + itk::CompensatedSummation m_SumOfPositivePixels; + itk::CompensatedSummation m_SumOfSquares; + itk::CompensatedSummation m_SumOfCubes; + itk::CompensatedSummation m_SumOfQuadruples; + RealType m_Sigma; + RealType m_Variance; + RealType m_MPP; + RealType m_Median; + RealType m_Uniformity; + RealType m_UPP; + RealType m_Entropy; + RealType m_Skewness; + RealType m_Kurtosis; + BoundingBoxType m_BoundingBox; + HistogramPointer m_Histogram; + }; + + using MapType = std::unordered_map; + using MapIterator = typename MapType::iterator; + using MapConstIterator = typename MapType::const_iterator; + + using ValidLabelValuesContainerType = std::vector; + const ValidLabelValuesContainerType& GetValidLabelValues() const; + + void SetHistogramParameters( + const std::unordered_map& sizes, + const std::unordered_map& lowerBounds, + const std::unordered_map& upperBounds); + + using LabelImageType = itk::Image; + + itkSetInputMacro(LabelInput, LabelImageType); + itkGetInputMacro(LabelInput, LabelImageType); + + bool HasLabel(LabelPixelType label) const; + unsigned int GetNumberOfObjects() const; + unsigned int GetNumberOfLabels() const; + + PixelType GetMinimum(LabelPixelType label) const; + PixelType GetMaximum(LabelPixelType label) const; + RealType GetMean(LabelPixelType label) const; + RealType GetSigma(LabelPixelType label) const; + RealType GetVariance(LabelPixelType label) const; + BoundingBoxType GetBoundingBox(LabelPixelType label) const; + RegionType GetRegion(LabelPixelType label) const; + RealType GetSum(LabelPixelType label) const; + RealType GetSumOfSquares(LabelPixelType label) const; + RealType GetSumOfCubes(LabelPixelType label) const; + RealType GetSumOfQuadruples(LabelPixelType label) const; + RealType GetSkewness(LabelPixelType label) const; + RealType GetKurtosis(LabelPixelType label) const; + RealType GetMPP(LabelPixelType label) const; + itk::SizeValueType GetCount(LabelPixelType label) const; + HistogramPointer GetHistogram(LabelPixelType label) const; + RealType GetEntropy(LabelPixelType label) const; + RealType GetUniformity(LabelPixelType label) const; + RealType GetUPP(LabelPixelType label) const; + RealType GetMedian(LabelPixelType label) const; + + protected: + LabelStatisticsImageFilter(); + ~LabelStatisticsImageFilter(); + + void BeforeStreamedGenerateData() override; + void ThreadedStreamedGenerateData(const RegionType&) override; + void AfterStreamedGenerateData() override; + + void PrintSelf(std::ostream& os, itk::Indent indent) const override; + + private: + const LabelStatistics& GetLabelStatistics(LabelPixelType label) const; + const LabelStatistics& GetLabelHistogramStatistics(LabelPixelType label) const; + + void MergeMap(MapType& map1, MapType& map2) const; + + MapType m_LabelStatistics; + ValidLabelValuesContainerType m_ValidLabelValues; + + bool m_ComputeHistograms; + std::unordered_map m_HistogramSizes; + std::unordered_map m_HistogramLowerBounds; + std::unordered_map m_HistogramUpperBounds; + + std::mutex m_Mutex; + }; +} + +#ifndef ITK_MANUAL_INSTANTIATION +#include +#endif + +#endif diff --git a/Modules/ImageStatistics/mitkLabelStatisticsImageFilter.hxx b/Modules/ImageStatistics/mitkLabelStatisticsImageFilter.hxx new file mode 100644 index 0000000000..367fd3f5c2 --- /dev/null +++ b/Modules/ImageStatistics/mitkLabelStatisticsImageFilter.hxx @@ -0,0 +1,546 @@ +/*============================================================================ + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center (DKFZ) +All rights reserved. + +Use of this source code is governed by a 3-clause BSD license that can be +found in the LICENSE file. + +============================================================================*/ + +#ifndef mitkLabelStatisticsImageFilter_hxx +#define mitkLabelStatisticsImageFilter_hxx + +#include "mitkLabelStatisticsImageFilter.h" + +#include +#include + +template +mitk::LabelStatisticsImageFilter::LabelStatistics::LabelStatistics() + : m_Count(0), + m_CountOfPositivePixels(0), + m_Min(itk::NumericTraits::max()), + m_Max(itk::NumericTraits::NonpositiveMin()), + m_Mean(0), + m_Sum(0), + m_SumOfPositivePixels(0), + m_SumOfSquares(0), + m_SumOfCubes(0), + m_SumOfQuadruples(0), + m_Sigma(0), + m_Variance(0), + m_MPP(0), + m_Median(0), + m_Uniformity(0), + m_UPP(0), + m_Entropy(0), + m_Skewness(0), + m_Kurtosis(0) +{ + m_BoundingBox.resize(ImageDimension * 2); + + for (std::remove_const_t i = 0; i < ImageDimension * 2; i += 2) + { + m_BoundingBox[i] = itk::NumericTraits::max(); + m_BoundingBox[i + 1] = itk::NumericTraits::NonpositiveMin(); + } +} + +template +mitk::LabelStatisticsImageFilter::LabelStatistics::LabelStatistics(unsigned int size, RealType lowerBound, RealType upperBound) + : LabelStatistics() +{ + typename HistogramType::SizeType histogramSize; + histogramSize.SetSize(1); + histogramSize[0] = size; + + typename HistogramType::MeasurementVectorType histogramLowerBound; + histogramLowerBound.SetSize(1); + histogramLowerBound[0] = lowerBound; + + typename HistogramType::MeasurementVectorType histogramUpperBound; + histogramUpperBound.SetSize(1); + histogramUpperBound[0] = upperBound; + + m_Histogram = HistogramType::New(); + m_Histogram->SetMeasurementVectorSize(1); + m_Histogram->Initialize(histogramSize, histogramLowerBound, histogramUpperBound); +} + +template +mitk::LabelStatisticsImageFilter::LabelStatistics::~LabelStatistics() +{ +} + +template +mitk::LabelStatisticsImageFilter::LabelStatisticsImageFilter() + : m_ComputeHistograms(false) +{ + this->AddRequiredInputName("LabelInput"); +} + +template +mitk::LabelStatisticsImageFilter::~LabelStatisticsImageFilter() +{ +} + +template +auto mitk::LabelStatisticsImageFilter::BeforeStreamedGenerateData() -> void +{ + this->AllocateOutputs(); + m_LabelStatistics.clear(); +} + +template +auto mitk::LabelStatisticsImageFilter::ThreadedStreamedGenerateData(const RegionType& region) -> void +{ + if (0 == region.GetSize(0)) + return; + + MapType localStats; + + typename HistogramType::MeasurementVectorType histogramMeasurement(1); + typename HistogramType::IndexType histogramIndex(1); + + using TLabelImage = itk::Image; + + itk::ImageLinearConstIteratorWithIndex it(this->GetInput(), region); + itk::ImageScanlineConstIterator labelIt(this->GetLabelInput(), region); + + auto mapIt = localStats.end(); + + while (!it.IsAtEnd()) + { + while (!it.IsAtEndOfLine()) + { + const auto& value = static_cast(it.Get()); + const auto& index = it.GetIndex(); + const auto& label = labelIt.Get(); + + mapIt = localStats.find(label); + + if (mapIt == localStats.end()) + { + mapIt = m_ComputeHistograms + ? localStats.emplace(label, LabelStatistics(m_HistogramSizes[label], m_HistogramLowerBounds[label], m_HistogramUpperBounds[label])).first + : localStats.emplace(label, LabelStatistics()).first; + } + + auto& labelStats = mapIt->second; + + labelStats.m_Min = std::min(labelStats.m_Min, value); + labelStats.m_Max = std::max(labelStats.m_Max, value); + labelStats.m_Sum += value; + auto squareValue = value * value; + labelStats.m_SumOfSquares += squareValue; + labelStats.m_SumOfCubes += squareValue * value; + labelStats.m_SumOfQuadruples += squareValue * squareValue; + ++labelStats.m_Count; + + if (0 < value) + { + labelStats.m_SumOfPositivePixels += value; + ++labelStats.m_CountOfPositivePixels; + } + + for (unsigned int i = 0; i < (ImageDimension * 2); i += 2) + { + labelStats.m_BoundingBox[i] = std::min(labelStats.m_BoundingBox[i], index[i / 2]); + labelStats.m_BoundingBox[i + 1] = std::max(labelStats.m_BoundingBox[i + 1], index[i / 2]); + } + + if (m_ComputeHistograms) + { + histogramMeasurement[0] = value; + labelStats.m_Histogram->GetIndex(histogramMeasurement, histogramIndex); + labelStats.m_Histogram->IncreaseFrequencyOfIndex(histogramIndex, 1); + } + + ++labelIt; + ++it; + } + + labelIt.NextLine(); + it.NextLine(); + } + + // Merge localStats and m_LabelStatistics concurrently safe in a local copy + + while (true) + { + std::unique_lock lock(m_Mutex); + + if (m_LabelStatistics.empty()) + { + std::swap(m_LabelStatistics, localStats); + break; + } + else + { + MapType toMerge; + std::swap(m_LabelStatistics, toMerge); + lock.unlock(); + this->MergeMap(localStats, toMerge); + } + } +} + +template +auto mitk::LabelStatisticsImageFilter::AfterStreamedGenerateData() -> void +{ + Superclass::AfterStreamedGenerateData(); + + m_ValidLabelValues.clear(); + m_ValidLabelValues.reserve(m_LabelStatistics.size()); + + for (auto& val : m_LabelStatistics) + { + m_ValidLabelValues.push_back(val.first); + auto& stats = val.second; + + const auto& sum = stats.m_Sum.GetSum(); + const auto& sumOfSquares = stats.m_SumOfSquares.GetSum(); + const auto& sumOfCubes = stats.m_SumOfCubes.GetSum(); + const auto& sumOfQuadruples = stats.m_SumOfQuadruples.GetSum(); + const auto& sumOfPositivePixels = stats.m_SumOfPositivePixels.GetSum(); + + const RealType count = stats.m_Count; + const RealType countOfPositivePixels = stats.m_CountOfPositivePixels; + + stats.m_Mean = sum / count; + const auto& mean = stats.m_Mean; + + if (count > 1) + { + auto sumSquared = sum * sum; + stats.m_Variance = (sumOfSquares - sumSquared / count) / (count - 1.0); + } + else + { + stats.m_Variance = 0.0; + } + + stats.m_Sigma = std::sqrt(stats.m_Variance); + + const auto secondMoment = sumOfSquares / count; + const auto thirdMoment = sumOfCubes / count; + const auto fourthMoment = sumOfQuadruples / count; + + stats.m_Skewness = (thirdMoment - 3 * secondMoment * mean + 2 * std::pow(mean, 3)) / std::pow(secondMoment - std::pow(mean, 2), 1.5); + stats.m_Kurtosis = (fourthMoment - 4 * thirdMoment * mean + 6 * secondMoment * std::pow(mean, 2) - 3 * std::pow(mean, 4)) / std::pow(secondMoment - std::pow(mean, 2), 2); + stats.m_MPP = sumOfPositivePixels / countOfPositivePixels; + + if (m_ComputeHistograms) + { + mitk::HistogramStatisticsCalculator histogramStatisticsCalculator; + histogramStatisticsCalculator.SetHistogram(stats.m_Histogram); + histogramStatisticsCalculator.CalculateStatistics(); + + stats.m_Entropy = histogramStatisticsCalculator.GetEntropy(); + stats.m_Uniformity = histogramStatisticsCalculator.GetUniformity(); + stats.m_UPP = histogramStatisticsCalculator.GetUPP(); + stats.m_Median = histogramStatisticsCalculator.GetMedian(); + } + } +} + +template +auto mitk::LabelStatisticsImageFilter::SetHistogramParameters( + const std::unordered_map& sizes, + const std::unordered_map& lowerBounds, + const std::unordered_map& upperBounds) -> void +{ + bool modified = false; + + if (m_HistogramSizes != sizes) + { + m_HistogramSizes = sizes; + modified = true; + } + + if (m_HistogramLowerBounds != lowerBounds) + { + m_HistogramLowerBounds = lowerBounds; + modified = true; + } + + if (m_HistogramUpperBounds != upperBounds) + { + m_HistogramUpperBounds = upperBounds; + modified = true; + } + + m_ComputeHistograms = true; + + if (modified) + this->Modified(); +} + +template +auto mitk::LabelStatisticsImageFilter::MergeMap(MapType& map1, MapType& map2) const -> void +{ + for (auto& elem2 : map2) + { + auto iter1 = map1.find(elem2.first); + + if (map1.end() == iter1) + { + map1.emplace(elem2.first, std::move(elem2.second)); + } + else + { + const auto label = iter1->first; + auto& stats1 = iter1->second; + auto& stats2 = elem2.second; + + stats1.m_Count += stats2.m_Count; + stats1.m_Sum += stats2.m_Sum; + stats1.m_SumOfSquares += stats2.m_SumOfSquares; + + stats1.m_Min = std::min(stats1.m_Min, stats2.m_Min); + stats1.m_Max = std::max(stats1.m_Max, stats2.m_Max); + + stats1.m_Sum += stats2.m_Sum; + stats1.m_SumOfSquares += stats2.m_SumOfSquares; + stats1.m_SumOfCubes += stats2.m_SumOfCubes; + stats1.m_SumOfQuadruples += stats2.m_SumOfQuadruples; + stats1.m_Count += stats2.m_Count; + stats1.m_SumOfPositivePixels += stats2.m_SumOfPositivePixels; + stats1.m_CountOfPositivePixels += stats2.m_CountOfPositivePixels; + + for (unsigned int i = 0; i < (ImageDimension * 2); i += 2) + { + stats1.m_BoundingBox[i] = std::min(stats1.m_BoundingBox[i], stats2.m_BoundingBox[i]); + stats1.m_BoundingBox[i + 1] = std::max(stats1.m_BoundingBox[i + 1], stats2.m_BoundingBox[i + 1]); + } + + if (m_ComputeHistograms) + { + typename HistogramType::IndexType index; + index.SetSize(1); + + const auto histogramSize = m_HistogramSizes.at(label); + + for (unsigned int bin = 0; bin < histogramSize; ++bin) + { + index[0] = bin; + stats1.m_Histogram->IncreaseFrequency(bin, stats2.m_Histogram->GetFrequency(bin)); + } + } + } + } +} + +template +auto mitk::LabelStatisticsImageFilter::HasLabel(LabelPixelType label) const -> bool +{ + return m_LabelStatistics.find(label) != m_LabelStatistics.end(); +} + +template +auto mitk::LabelStatisticsImageFilter::GetNumberOfObjects() const -> unsigned int +{ + return static_cast(m_LabelStatistics.size()); +} + +template +auto mitk::LabelStatisticsImageFilter::GetNumberOfLabels() const -> unsigned int +{ + return this->GetNumberOfObjects(); +} + +template +auto mitk::LabelStatisticsImageFilter::GetValidLabelValues() const -> const ValidLabelValuesContainerType& +{ + return m_ValidLabelValues; +} + +template +auto mitk::LabelStatisticsImageFilter::GetLabelStatistics(LabelPixelType label) const -> const LabelStatistics& +{ + MapConstIterator it = m_LabelStatistics.find(label); + + if (m_LabelStatistics.end() != it) + return it->second; + + mitkThrow() << "Label " << label << " does not exist"; +} + +template +auto mitk::LabelStatisticsImageFilter::GetLabelHistogramStatistics(LabelPixelType label) const -> const LabelStatistics& +{ + const auto& labelStatistics = this->GetLabelStatistics(label); + + if (m_ComputeHistograms && labelStatistics.m_Histogram.IsNotNull()) + return labelStatistics; + + mitkThrow() << "Histogram was not computed for label " << label; +} + +template +auto mitk::LabelStatisticsImageFilter::GetMinimum(LabelPixelType label) const -> PixelType +{ + const auto& labelStatistics = this->GetLabelStatistics(label); + return labelStatistics.m_Min; +} + +template +auto mitk::LabelStatisticsImageFilter::GetMaximum(LabelPixelType label) const -> PixelType +{ + const auto& labelStatistics = this->GetLabelStatistics(label); + return labelStatistics.m_Max; +} + +template +auto mitk::LabelStatisticsImageFilter::GetMean(LabelPixelType label) const -> RealType +{ + const auto& labelStatistics = this->GetLabelStatistics(label); + return labelStatistics.m_Mean; +} + +template +auto mitk::LabelStatisticsImageFilter::GetSigma(LabelPixelType label) const -> RealType +{ + const auto& labelStatistics = this->GetLabelStatistics(label); + return labelStatistics.m_Sigma; +} + +template +auto mitk::LabelStatisticsImageFilter::GetVariance(LabelPixelType label) const -> RealType +{ + const auto& labelStatistics = this->GetLabelStatistics(label); + return labelStatistics.m_Variance; +} + +template +auto mitk::LabelStatisticsImageFilter::GetSkewness(LabelPixelType label) const ->RealType +{ + const auto& labelStatistics = this->GetLabelStatistics(label); + return labelStatistics.m_Skewness; +} + +template +auto mitk::LabelStatisticsImageFilter::GetKurtosis(LabelPixelType label) const -> RealType +{ + const auto& labelStatistics = this->GetLabelStatistics(label); + return labelStatistics.m_Kurtosis; +} + +template +auto mitk::LabelStatisticsImageFilter::GetMPP(LabelPixelType label) const -> RealType +{ + const auto& labelStatistics = this->GetLabelStatistics(label); + return labelStatistics.m_MPP; +} + +template +auto mitk::LabelStatisticsImageFilter::GetCount(LabelPixelType label) const-> itk::SizeValueType +{ + const auto& labelStatistics = this->GetLabelStatistics(label); + return labelStatistics.m_Count; +} + +template +auto mitk::LabelStatisticsImageFilter::GetHistogram(LabelPixelType label) const -> HistogramPointer +{ + const auto& labelStatistics = this->GetLabelHistogramStatistics(label); + return labelStatistics.m_Histogram; +} + +template +auto mitk::LabelStatisticsImageFilter::GetEntropy(LabelPixelType label) const -> RealType +{ + const auto& labelStatistics = this->GetLabelHistogramStatistics(label); + return labelStatistics.m_Entropy; +} + +template +auto mitk::LabelStatisticsImageFilter::GetUniformity(LabelPixelType label) const -> RealType +{ + const auto& labelStatistics = this->GetLabelHistogramStatistics(label); + return labelStatistics.m_Uniformity; +} + +template +auto mitk::LabelStatisticsImageFilter::GetUPP(LabelPixelType label) const -> RealType +{ + const auto& labelStatistics = this->GetLabelHistogramStatistics(label); + return labelStatistics.m_UPP; +} + +template +auto mitk::LabelStatisticsImageFilter::GetMedian(LabelPixelType label) const -> RealType +{ + const auto& labelStatistics = this->GetLabelHistogramStatistics(label); + return labelStatistics.m_Median; +} + +template +auto mitk::LabelStatisticsImageFilter::GetSum(LabelPixelType label) const -> RealType +{ + const auto& labelStatistics = this->GetLabelStatistics(label); + return labelStatistics.m_Sum.GetSum(); +} + +template +auto mitk::LabelStatisticsImageFilter::GetSumOfSquares(LabelPixelType label) const -> RealType +{ + const auto& labelStatistics = this->GetLabelStatistics(label); + return labelStatistics.m_SumOfSquares.GetSum(); +} + +template +auto mitk::LabelStatisticsImageFilter::GetSumOfCubes(LabelPixelType label) const -> RealType +{ + const auto& labelStatistics = this->GetLabelStatistics(label); + return labelStatistics.m_SumOfCubes.GetSum(); +} + +template +auto mitk::LabelStatisticsImageFilter::GetSumOfQuadruples(LabelPixelType label) const -> RealType +{ + const auto& labelStatistics = this->GetLabelStatistics(label); + return labelStatistics.m_SumOfQuadruples.GetSum(); +} + +template +auto mitk::LabelStatisticsImageFilter::GetBoundingBox(LabelPixelType label) const -> BoundingBoxType +{ + const auto& labelStatistics = this->GetLabelStatistics(label); + return labelStatistics.m_BoundingBox; +} + +template +auto mitk::LabelStatisticsImageFilter::GetRegion(LabelPixelType label) const -> RegionType +{ + const auto& labelStatistics = this->GetLabelStatistics(label); + + IndexType index; + SizeType size; + + for (unsigned int i = 0; i < ImageDimension; ++i) + { + index[i] = labelStatistics.m_BoundingBox[2 * i]; + size[i] = labelStatistics.m_BoundingBox[2 * i + 1] - labelStatistics.m_BoundingBox[2 * i] + 1; + } + + RegionType region; + region.SetIndex(index); + region.SetSize(size); + + return region; +} + +template +void mitk::LabelStatisticsImageFilter::PrintSelf(std::ostream& os, itk::Indent indent) const +{ + Superclass::PrintSelf(os, indent); + + os << indent << "Number of labels: " << m_LabelStatistics.size() << std::endl; + os << indent << "Compute histograms: " << m_ComputeHistograms << std::endl; +} + +#endif diff --git a/Modules/ImageStatistics/mitkStatisticsImageFilter.h b/Modules/ImageStatistics/mitkStatisticsImageFilter.h index 18302d41d5..d971f0ea50 100644 --- a/Modules/ImageStatistics/mitkStatisticsImageFilter.h +++ b/Modules/ImageStatistics/mitkStatisticsImageFilter.h @@ -1,191 +1,142 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ // This file is based on ITK's itkStatisticsImageFilter.h #ifndef mitkStatisticsImageFilter #define mitkStatisticsImageFilter #include #include #include #include #include #include #include #include namespace mitk { template class StatisticsImageFilter : public itk::ImageSink { public: using Self = StatisticsImageFilter; using Superclass = itk::ImageSink; using Pointer = itk::SmartPointer; using ConstPointer = itk::SmartPointer; itkFactorylessNewMacro(Self); itkTypeMacro(StatisticsImageFilter, itk::ImageSink); - using InputImagePointer = typename TInputImage::Pointer; - using RegionType = typename TInputImage::RegionType; - using SizeType = typename TInputImage::SizeType; - using IndexType = typename TInputImage::IndexType; using PixelType = typename TInputImage::PixelType; - static constexpr unsigned int ImageDimension = TInputImage::ImageDimension; - using RealType = typename itk::NumericTraits::RealType; using HistogramType = typename itk::Statistics::Histogram; using HistogramPointer = itk::SmartPointer; using DataObjectPointer = typename itk::DataObject::Pointer; template using SimpleDataObjectDecorator = itk::SimpleDataObjectDecorator; using RealObjectType = SimpleDataObjectDecorator; using PixelObjectType = SimpleDataObjectDecorator; - /** Return the computed Minimum. */ itkGetDecoratedOutputMacro(Minimum, PixelType); - - /** Return the computed Maximum. */ itkGetDecoratedOutputMacro(Maximum, PixelType); - - /** Return the computed Mean. */ itkGetDecoratedOutputMacro(Mean, RealType); - - /** Return the computed Standard Deviation. */ itkGetDecoratedOutputMacro(Sigma, RealType); - - /** Return the computed Variance. */ itkGetDecoratedOutputMacro(Variance, RealType); - - /** Return the computed Sum. */ itkGetDecoratedOutputMacro(Sum, RealType); - - /** Return the computed Sum of Squares. */ itkGetDecoratedOutputMacro(SumOfSquares, RealType); - - /** Return the computed Sum of Cubes. */ itkGetDecoratedOutputMacro(SumOfCubes, RealType); - - /** Return the computed Sum of Quadruples. */ itkGetDecoratedOutputMacro(SumOfQuadruples, RealType); - - /** Return the computed Skewness. */ itkGetDecoratedOutputMacro(Skewness, RealType); - - /** Return the computed Kurtosis. */ itkGetDecoratedOutputMacro(Kurtosis, RealType); - - /** Return the computed Mean of Positive Pixels. */ itkGetDecoratedOutputMacro(MPP, RealType); - - /** Return the computed Histogram. */ itkGetDecoratedOutputMacro(Histogram, HistogramPointer); - - /** Return the computed Entropy. */ itkGetDecoratedOutputMacro(Entropy, RealType); - - /** Return the computed Uniformity. */ itkGetDecoratedOutputMacro(Uniformity, RealType); - - /** Return the computed UPP. */ itkGetDecoratedOutputMacro(UPP, RealType); - - /** Return the computed Median. */ itkGetDecoratedOutputMacro(Median, RealType); void SetHistogramParameters(unsigned int size, RealType lowerBound, RealType upperBound); using DataObjectIdentifierType = itk::ProcessObject::DataObjectIdentifierType; using Superclass::MakeOutput; /** Make a DataObject of the correct type to be used as the specified output. */ DataObjectPointer MakeOutput(const DataObjectIdentifierType& name) override; -#ifdef ITK_USE_CONCEPT_CHECKING - // Begin concept checking - itkConceptMacro(InputHasNumericTraitsCheck, (itk::Concept::HasNumericTraits)); - // End concept checking -#endif - protected: StatisticsImageFilter(); ~StatisticsImageFilter(); itkSetDecoratedOutputMacro(Minimum, PixelType); itkSetDecoratedOutputMacro(Maximum, PixelType); itkSetDecoratedOutputMacro(Mean, RealType); itkSetDecoratedOutputMacro(Sigma, RealType); itkSetDecoratedOutputMacro(Variance, RealType); itkSetDecoratedOutputMacro(Sum, RealType); itkSetDecoratedOutputMacro(SumOfSquares, RealType); itkSetDecoratedOutputMacro(SumOfCubes, RealType); itkSetDecoratedOutputMacro(SumOfQuadruples, RealType); itkSetDecoratedOutputMacro(Skewness, RealType); itkSetDecoratedOutputMacro(Kurtosis, RealType); itkSetDecoratedOutputMacro(MPP, RealType); itkSetDecoratedOutputMacro(Histogram, HistogramPointer); itkSetDecoratedOutputMacro(Entropy, RealType); itkSetDecoratedOutputMacro(Uniformity, RealType); itkSetDecoratedOutputMacro(UPP, RealType); itkSetDecoratedOutputMacro(Median, RealType); - /** Initialize some accumulators before the threads run. */ void BeforeStreamedGenerateData() override; - void ThreadedStreamedGenerateData(const RegionType&) override; - - /** Set outputs to computed values from all regions. */ void AfterStreamedGenerateData() override; void PrintSelf(std::ostream& os, itk::Indent indent) const override; private: HistogramPointer CreateInitializedHistogram() const; - bool m_CalculateHistogram; + bool m_ComputeHistogram; unsigned int m_HistogramSize; - PixelType m_HistogramLowerBound; - PixelType m_HistogramUpperBound; + RealType m_HistogramLowerBound; + RealType m_HistogramUpperBound; HistogramPointer m_Histogram; itk::CompensatedSummation m_Sum; itk::CompensatedSummation m_SumOfPositivePixels; itk::CompensatedSummation m_SumOfSquares; itk::CompensatedSummation m_SumOfCubes; itk::CompensatedSummation m_SumOfQuadruples; itk::SizeValueType m_Count; itk::SizeValueType m_CountOfPositivePixels; PixelType m_Min; PixelType m_Max; std::mutex m_Mutex; }; } #ifndef ITK_MANUAL_INSTANTIATION #include #endif #endif diff --git a/Modules/ImageStatistics/mitkStatisticsImageFilter.hxx b/Modules/ImageStatistics/mitkStatisticsImageFilter.hxx index b0ea7f5486..8c44f722dc 100644 --- a/Modules/ImageStatistics/mitkStatisticsImageFilter.hxx +++ b/Modules/ImageStatistics/mitkStatisticsImageFilter.hxx @@ -1,325 +1,325 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #ifndef mitkStatisticsImageFilter_hxx #define mitkStatisticsImageFilter_hxx #include #include #include template mitk::StatisticsImageFilter::StatisticsImageFilter() - : m_CalculateHistogram(false), + : m_ComputeHistogram(false), m_HistogramSize(0), m_HistogramLowerBound(itk::NumericTraits::NonpositiveMin()), m_HistogramUpperBound(itk::NumericTraits::max()), m_Sum(1), m_SumOfPositivePixels(1), m_SumOfSquares(1), m_SumOfCubes(1), m_SumOfQuadruples(1), m_Count(1), m_CountOfPositivePixels(1), m_Min(1), m_Max(1) { this->SetNumberOfRequiredInputs(1); this->SetMinimum(itk::NumericTraits::max()); this->SetMaximum(itk::NumericTraits::NonpositiveMin()); this->SetMean(itk::NumericTraits::max()); this->SetSigma(itk::NumericTraits::max()); this->SetVariance(itk::NumericTraits::max()); this->SetSum(0); this->SetSumOfSquares(0); this->SetSumOfCubes(0); this->SetSumOfQuadruples(0); this->SetSkewness(0); this->SetKurtosis(0); this->SetMPP(0); this->SetEntropy(-1.0); this->SetUniformity(0); this->SetUPP(0); this->SetMedian(0); } template mitk::StatisticsImageFilter::~StatisticsImageFilter() { } template typename mitk::StatisticsImageFilter::DataObjectPointer mitk::StatisticsImageFilter::MakeOutput(const DataObjectIdentifierType& name) { if (name == "Minimum" || name == "Maximum") { return PixelObjectType::New(); } if (name == "Mean" || name == "Sigma" || name == "Variance" || name == "Sum" || name == "SumOfSquares" || name == "SumOfCubes" || name == "SumOfQuadruples" || name == "Skewness" || name == "Kurtosis" || name == "MPP" || name == "Entropy" || name == "Uniformity" || name == "UPP" || name == "Median") { return RealObjectType::New(); } if (name == "Histogram") { return HistogramType::New(); } return Superclass::MakeOutput(name); } template void mitk::StatisticsImageFilter::SetHistogramParameters(unsigned int size, RealType lowerBound, RealType upperBound) { bool modified = false; if (m_HistogramSize != size) { m_HistogramSize = size; modified = true; } if (m_HistogramLowerBound != lowerBound) { m_HistogramLowerBound = lowerBound; modified = true; } if (m_HistogramUpperBound != upperBound) { m_HistogramUpperBound = upperBound; modified = true; } - m_CalculateHistogram = true; + m_ComputeHistogram = true; if (modified) this->Modified(); } template auto mitk::StatisticsImageFilter::CreateInitializedHistogram() const -> HistogramPointer { typename HistogramType::SizeType size; size.SetSize(1); size.Fill(m_HistogramSize); typename HistogramType::MeasurementVectorType lowerBound; lowerBound.SetSize(1); lowerBound.Fill(m_HistogramLowerBound); typename HistogramType::MeasurementVectorType upperBound; upperBound.SetSize(1); upperBound.Fill(m_HistogramUpperBound); auto histogram = HistogramType::New(); histogram->SetMeasurementVectorSize(1); histogram->Initialize(size, lowerBound, upperBound); return histogram; } template void mitk::StatisticsImageFilter::BeforeStreamedGenerateData() { Superclass::BeforeStreamedGenerateData(); m_Sum = 0; m_SumOfPositivePixels = 0; m_SumOfSquares = 0; m_SumOfCubes = 0; m_SumOfQuadruples = 0; m_Count = 0; m_CountOfPositivePixels = 0; m_Min = itk::NumericTraits::max(); m_Max = itk::NumericTraits::NonpositiveMin(); - if (m_CalculateHistogram) + if (m_ComputeHistogram) m_Histogram = this->CreateInitializedHistogram(); } template void mitk::StatisticsImageFilter::ThreadedStreamedGenerateData(const RegionType& regionForThread) { itk::CompensatedSummation sum = 0; itk::CompensatedSummation sumOfPositivePixels = 0; itk::CompensatedSummation sumOfSquares = 0; itk::CompensatedSummation sumOfCubes = 0; itk::CompensatedSummation sumOfQuadruples = 0; itk::SizeValueType count = 0; itk::SizeValueType countOfPositivePixels = 0; auto min = itk::NumericTraits::max(); auto max = itk::NumericTraits::NonpositiveMin(); RealType realValue = 0; RealType squareValue = 0; HistogramPointer histogram; typename HistogramType::MeasurementVectorType histogramMeasurement; typename HistogramType::IndexType histogramIndex; - if (m_CalculateHistogram) // Initialize histogram + if (m_ComputeHistogram) // Initialize histogram { histogram = this->CreateInitializedHistogram(); histogramMeasurement.SetSize(1); } itk::ImageScanlineConstIterator it(this->GetInput(), regionForThread); while (!it.IsAtEnd()) { while (!it.IsAtEndOfLine()) { const auto& value = it.Get(); realValue = static_cast(value); - if (m_CalculateHistogram) // Compute histogram + if (m_ComputeHistogram) // Compute histogram { histogramMeasurement[0] = realValue; histogram->GetIndex(histogramMeasurement, histogramIndex); histogram->IncreaseFrequencyOfIndex(histogramIndex, 1); } min = std::min(min, value); max = std::max(max, value); squareValue = realValue * realValue; sum += realValue; sumOfSquares += squareValue; sumOfCubes += squareValue * realValue; sumOfQuadruples += squareValue * squareValue; ++count; if (0 < realValue) { sumOfPositivePixels += realValue; ++countOfPositivePixels; } ++it; } it.NextLine(); } std::lock_guard mutexHolder(m_Mutex); - if (m_CalculateHistogram) // Merge histograms + if (m_ComputeHistogram) // Merge histograms { typename HistogramType::ConstIterator histogramIt = histogram->Begin(); typename HistogramType::ConstIterator histogramEnd = histogram->End(); while (histogramIt != histogramEnd) { m_Histogram->GetIndex(histogramIt.GetMeasurementVector(), histogramIndex); m_Histogram->IncreaseFrequencyOfIndex(histogramIndex, histogramIt.GetFrequency()); ++histogramIt; } } m_Sum += sum; m_SumOfPositivePixels += sumOfPositivePixels; m_SumOfSquares += sumOfSquares; m_SumOfCubes += sumOfCubes; m_SumOfQuadruples += sumOfQuadruples; m_Count += count; m_CountOfPositivePixels += countOfPositivePixels; m_Min = std::min(min, m_Min); m_Max = std::max(max, m_Max); } template void mitk::StatisticsImageFilter::AfterStreamedGenerateData() { Superclass::AfterStreamedGenerateData(); const RealType sum = m_Sum.GetSum(); const RealType sumOfPositivePixels = m_SumOfPositivePixels.GetSum(); const RealType sumOfSquares = m_SumOfSquares.GetSum(); const RealType sumOfCubes = m_SumOfCubes.GetSum(); const RealType sumOfQuadruples = m_SumOfQuadruples.GetSum(); const itk::SizeValueType count = m_Count; const itk::SizeValueType countOfPositivePixels = m_CountOfPositivePixels; const PixelType minimum = m_Min; const PixelType maximum = m_Max; const RealType mean = sum / static_cast(count); const RealType variance = (sumOfSquares - (sum * sum / static_cast(count))) / (static_cast(count) - 1); const RealType sigma = std::sqrt(variance); const RealType secondMoment = sumOfSquares / static_cast(count); const RealType thirdMoment = sumOfCubes / static_cast(count); const RealType fourthMoment = sumOfQuadruples / static_cast(count); const RealType skewness = (thirdMoment - 3 * secondMoment * mean + 2 * std::pow(mean, 3)) / std::pow(secondMoment - std::pow(mean, 2), 1.5); const RealType kurtosis = (fourthMoment - 4 * thirdMoment * mean + 6 * secondMoment * std::pow(mean, 2) - 3 * std::pow(mean, 4)) / std::pow(secondMoment - std::pow(mean, 2), 2); const RealType meanOfPositivePixels = sumOfPositivePixels / static_cast(countOfPositivePixels); this->SetMinimum(minimum); this->SetMaximum(maximum); this->SetMean(mean); this->SetSigma(sigma); this->SetVariance(variance); this->SetSum(sum); this->SetSumOfSquares(sumOfSquares); this->SetSumOfCubes(sumOfCubes); this->SetSumOfQuadruples(sumOfQuadruples); this->SetSkewness(skewness); this->SetKurtosis(kurtosis); this->SetMPP(meanOfPositivePixels); - if (m_CalculateHistogram) + if (m_ComputeHistogram) { this->SetHistogram(m_Histogram); mitk::HistogramStatisticsCalculator histogramStatisticsCalculator; histogramStatisticsCalculator.SetHistogram(m_Histogram); histogramStatisticsCalculator.CalculateStatistics(); this->SetEntropy(histogramStatisticsCalculator.GetEntropy()); this->SetUniformity(histogramStatisticsCalculator.GetUniformity()); this->SetUPP(histogramStatisticsCalculator.GetUPP()); this->SetMedian(histogramStatisticsCalculator.GetMedian()); } } template void mitk::StatisticsImageFilter::PrintSelf(std::ostream& os, itk::Indent indent) const { Superclass::PrintSelf(os, indent); os << indent << "SumOfCubes: " << this->GetSumOfCubes() << std::endl; os << indent << "SumOfQuadruples: " << this->GetSumOfQuadruples() << std::endl; os << indent << "Skewness: " << this->GetSkewness() << std::endl; os << indent << "Kurtosis: " << this->GetKurtosis() << std::endl; os << indent << "MPP: " << this->GetMPP() << std::endl; os << indent << "Entropy: " << this->GetEntropy() << std::endl; os << indent << "Uniformity: " << this->GetUniformity() << std::endl; os << indent << "UPP: " << this->GetUPP() << std::endl; os << indent << "Median: " << this->GetMedian() << std::endl; } #endif