diff --git a/Modules/ImageStatistics/Testing/mitkImageStatisticsTextureAnalysisTest.cpp b/Modules/ImageStatistics/Testing/mitkImageStatisticsTextureAnalysisTest.cpp index f55425c872..4b3ad225f1 100644 --- a/Modules/ImageStatistics/Testing/mitkImageStatisticsTextureAnalysisTest.cpp +++ b/Modules/ImageStatistics/Testing/mitkImageStatisticsTextureAnalysisTest.cpp @@ -1,230 +1,230 @@ /*============================================================================ 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 "mitkTestingMacros.h" #include "itkImage.h" #include "mitkExtendedLabelStatisticsImageFilter.h" -#include "mitkExtendedStatisticsImageFilter.h" +#include "mitkStatisticsImageFilter.h" #include "mitkNumericConstants.h" /** \section Testing of Skewness and Kurtosis * This test class is for testing the added coefficients Skewness and Kurtosis * for the mitkExtendedLabelStatisticsImageFilter (Masked Images) and for the * mitkExtendedStatisticsImageFilter (Unmasked Images). Both filter will be tested * against two pictures. */ class mitkImageStatisticsTextureAnalysisTestClass { /** * \brief Explanation of the mitkImageStatisticsTextureAnalysisTestClass test class * * this test class produce test images and masking images with the method CreatingTestImageForDifferentLabelSize. * TestInstanceFortheMaskedStatisticsFilter and TestInstanceFortheUnmaskedStatisticsFilter are the two Instances * for the filters of masking and unmasking images. * TestofSkewnessAndKurtosisForMaskedImagesand TestofSkewnessAndKurtosisForUnmaskedImages are the correlated test * for the checking the values. */ public: typedef itk::Image< int,3 >ImageType; typedef ImageType::Pointer PointerOfImage; typedef itk::ExtendedLabelStatisticsImageFilter< ImageType, ImageType > LabelStatisticsFilterType; typedef LabelStatisticsFilterType::Pointer labelStatisticsFilterPointer; - typedef itk::ExtendedStatisticsImageFilter< ImageType > StatisticsFilterType; + typedef mitk::StatisticsImageFilter< ImageType > StatisticsFilterType; typedef StatisticsFilterType::Pointer StatisticsFilterPointer; ImageType::Pointer CreatingTestImageForDifferentLabelSize( int factorOfDividingThePicture, int bufferValue, int labelValue) { ImageType::Pointer image = ImageType :: New(); ImageType::IndexType start; ImageType::SizeType size; start[0] = 0; start[1] = 0; start[2] = 0; size[0] = 100; size[1] = 100; size[2] = 100; ImageType:: RegionType region; region.SetSize( size ); region.SetIndex( start ); image->SetRegions(region); image->Allocate(); image->FillBuffer(bufferValue); for(unsigned int r = 0; r < 50; r++) { for(int c = 0; c < factorOfDividingThePicture; c++) { for(unsigned int l = 0; l < 100; l++) { ImageType::IndexType pixelIndex; pixelIndex[0] = r; pixelIndex[1] = c; pixelIndex[2] = l; image->SetPixel(pixelIndex, labelValue); } } } return image; } LabelStatisticsFilterType::Pointer TestInstanceFortheMaskedStatisticsFilter(ImageType::Pointer image, ImageType::Pointer maskImage) { LabelStatisticsFilterType::Pointer labelStatisticsFilter; labelStatisticsFilter = LabelStatisticsFilterType::New(); labelStatisticsFilter->SetInput( image ); labelStatisticsFilter->UseHistogramsOn(); labelStatisticsFilter->SetHistogramParameters( 20, -10, 10); labelStatisticsFilter->SetLabelInput( maskImage ); labelStatisticsFilter->Update(); return labelStatisticsFilter; } StatisticsFilterType::Pointer TestInstanceFortheUnmaskedStatisticsFilter(ImageType::Pointer image ) { StatisticsFilterType::Pointer StatisticsFilter; StatisticsFilter = StatisticsFilterType::New(); StatisticsFilter->SetInput( image ); StatisticsFilter->SetHistogramParameters( 20, -10, 10 ); StatisticsFilter->Update(); return StatisticsFilter; } //test for Skewness,Kurtosis and MPP for masked Images void TestofSkewnessKurtosisAndMPPForMaskedImages(LabelStatisticsFilterType::Pointer labelStatisticsFilter, double expectedSkewness, double expectedKurtosis, double expectedMPP) { // let's create an object of our class bool isSkewsnessLowerlimitCorrect = labelStatisticsFilter->GetSkewness( 1 )- expectedKurtosis+ std::pow(10,-3) <= expectedSkewness; bool isSkewsnessUpperlimitCorrect = labelStatisticsFilter->GetSkewness( 1 )+ expectedKurtosis+ std::pow(10,-3) >= expectedSkewness; MITK_TEST_CONDITION( isSkewsnessLowerlimitCorrect && isSkewsnessUpperlimitCorrect,"expectedSkewness: " << expectedSkewness << " actual Value: " << labelStatisticsFilter->GetSkewness( 1 ) ); bool isKurtosisUpperlimitCorrect = labelStatisticsFilter->GetKurtosis( 1 ) <= expectedKurtosis+ std::pow(10,-3); bool isKurtosisLowerlimitCorrect = expectedKurtosis- std::pow(10,-3) <= labelStatisticsFilter->GetKurtosis( 1 ); MITK_TEST_CONDITION( isKurtosisUpperlimitCorrect && isKurtosisLowerlimitCorrect,"expectedKurtosis: " << expectedKurtosis << " actual Value: " << labelStatisticsFilter->GetKurtosis( 1 ) ); MITK_TEST_CONDITION( ( expectedMPP - labelStatisticsFilter->GetMPP( 1 ) ) < 1, "expected MPP: " << expectedMPP << " actual Value: " << labelStatisticsFilter->GetMPP( 1 ) ); } //test for Entropy,Uniformity and UPP for masked Images void TestofEntropyUniformityAndUppForMaskedImages(LabelStatisticsFilterType::Pointer labelStatisticsFilter, double expectedEntropy, double expectedUniformity, double expectedUPP) { bool calculatedEntropyLowerLimit = labelStatisticsFilter->GetEntropy( 1 ) >= expectedEntropy - std::pow(10,-3); bool calculatedUniformityLowerLimit = labelStatisticsFilter->GetUniformity( 1 ) >= expectedUniformity - std::pow(10,-3); bool calculatedUppLowerLimit = labelStatisticsFilter->GetUPP( 1 ) >= expectedUPP - std::pow(10,-3); bool calculatedEntropyUpperLimit = labelStatisticsFilter->GetEntropy( 1 ) <= expectedEntropy + std::pow(10,-3); bool calculatedUniformityUpperLimit = labelStatisticsFilter->GetUniformity( 1 ) <= expectedUniformity + std::pow(10,-3); bool calculatedUppUpperLimit = labelStatisticsFilter->GetUPP( 1 ) <= expectedUPP + std::pow(10,-3); MITK_TEST_CONDITION( calculatedEntropyLowerLimit && calculatedEntropyUpperLimit, "expected Entropy: " << expectedEntropy << " actual Value: " << labelStatisticsFilter->GetEntropy( 1 ) ); MITK_TEST_CONDITION( calculatedUniformityLowerLimit && calculatedUniformityUpperLimit, "expected Uniformity: " << expectedUniformity << " actual Value: " << labelStatisticsFilter->GetUniformity( 1 ) ); MITK_TEST_CONDITION( calculatedUppLowerLimit && calculatedUppUpperLimit, "expected UPP: " << expectedUPP << " actual Value: " << labelStatisticsFilter->GetUPP( 1 ) ); } //test for Skewness,Kurtosis and MPP for unmasked Images void TestofSkewnessKurtosisAndMPPForUnmaskedImages(StatisticsFilterType::Pointer StatisticsFilter, double expectedSkewness, double expectedKurtosis, double expectedMPP) { // let's create an object of our class bool isSkewsnessLowerlimitCorrect = StatisticsFilter->GetSkewness()- expectedKurtosis+ std::pow(10,-3) <= expectedSkewness; bool isSkewsnessUpperlimitCorrect = StatisticsFilter->GetSkewness()+ expectedKurtosis+ std::pow(10,-3) >= expectedSkewness; MITK_TEST_CONDITION( isSkewsnessLowerlimitCorrect && isSkewsnessUpperlimitCorrect,"expectedSkewness: " << expectedSkewness << " actual Value: " << StatisticsFilter->GetSkewness() ); bool isKurtosisUpperlimitCorrect = StatisticsFilter->GetKurtosis() <= expectedKurtosis+ std::pow(10,-3); bool isKurtosisLowerlimitCorrect = expectedKurtosis- std::pow(10,-3) <= StatisticsFilter->GetKurtosis(); MITK_TEST_CONDITION( isKurtosisUpperlimitCorrect && isKurtosisLowerlimitCorrect,"expectedKurtosis: " << expectedKurtosis << " actual Value: " << StatisticsFilter->GetKurtosis() ); MITK_TEST_CONDITION( ( expectedMPP - StatisticsFilter->GetMPP() ) < mitk::eps, "expected MPP: " << expectedMPP << " actual Value: " << StatisticsFilter->GetMPP() ); } //test for Entropy,Uniformity and UPP for unmasked Images void TestofEntropyUniformityAndUppForUnmaskedImages(StatisticsFilterType::Pointer StatisticsFilter, double expectedEntropy, double expectedUniformity, double expectedUPP) { bool calculatedEntropyLowerLimit = StatisticsFilter->GetEntropy() >= expectedEntropy - std::pow(10,-3); bool calculatedUniformityLowerLimit = StatisticsFilter->GetUniformity() >= expectedUniformity - std::pow(10,-3); bool calculatedUppLowerLimit = StatisticsFilter->GetUPP() >= expectedUPP - std::pow(10,-3); bool calculatedEntropyUpperLimit = StatisticsFilter->GetEntropy() <= expectedEntropy + std::pow(10,-3); bool calculatedUniformityUpperLimit = StatisticsFilter->GetUniformity() <= expectedUniformity + std::pow(10,-3); bool calculatedUppUpperLimit = StatisticsFilter->GetUPP() <= expectedUPP + std::pow(10,-3); MITK_TEST_CONDITION( calculatedEntropyLowerLimit && calculatedEntropyUpperLimit, "expected Entropy: " << expectedEntropy << " actual Value: " << StatisticsFilter->GetEntropy() ); MITK_TEST_CONDITION( calculatedUniformityLowerLimit && calculatedUniformityUpperLimit, "expected Uniformity: " << expectedUniformity << " actual Value: " << StatisticsFilter->GetUniformity() ); MITK_TEST_CONDITION( calculatedUppLowerLimit && calculatedUppUpperLimit, "expected UPP: " << expectedUPP << " actual Value: " << StatisticsFilter->GetUPP() ); } }; int mitkImageStatisticsTextureAnalysisTest(int, char* []) { // always start with this! MITK_TEST_BEGIN("mitkImageStatisticsTextureAnalysisTest") mitkImageStatisticsTextureAnalysisTestClass testclassInstance; mitkImageStatisticsTextureAnalysisTestClass::PointerOfImage labelImage = testclassInstance.CreatingTestImageForDifferentLabelSize(100, 1, 1); mitkImageStatisticsTextureAnalysisTestClass::PointerOfImage image = testclassInstance.CreatingTestImageForDifferentLabelSize(100, 3, 2); mitkImageStatisticsTextureAnalysisTestClass::PointerOfImage image2 = testclassInstance.CreatingTestImageForDifferentLabelSize(50, 3, 2); //test for masked images mitkImageStatisticsTextureAnalysisTestClass::labelStatisticsFilterPointer mitkLabelFilter= testclassInstance.TestInstanceFortheMaskedStatisticsFilter( image,labelImage); testclassInstance.TestofSkewnessKurtosisAndMPPForMaskedImages(mitkLabelFilter, 0, 0.999998, 2.5); testclassInstance.TestofEntropyUniformityAndUppForMaskedImages(mitkLabelFilter, 1, 0.5, 0.5); mitkImageStatisticsTextureAnalysisTestClass::labelStatisticsFilterPointer mitkLabelFilter2= testclassInstance.TestInstanceFortheMaskedStatisticsFilter( image2,labelImage); testclassInstance.TestofSkewnessKurtosisAndMPPForMaskedImages(mitkLabelFilter2, -1.1547, 2.33333, 2.75); testclassInstance.TestofEntropyUniformityAndUppForMaskedImages(mitkLabelFilter2, 0.811278, 0.625, 0.625); //test for unmasked images mitkImageStatisticsTextureAnalysisTestClass::StatisticsFilterPointer mitkFilter= testclassInstance.TestInstanceFortheUnmaskedStatisticsFilter( image); testclassInstance.TestofSkewnessKurtosisAndMPPForUnmaskedImages(mitkFilter, 0, 0.999998, 2.5); testclassInstance.TestofEntropyUniformityAndUppForUnmaskedImages(mitkFilter, 1, 0.5, 0.5); mitkImageStatisticsTextureAnalysisTestClass::StatisticsFilterPointer mitkFilter2= testclassInstance.TestInstanceFortheUnmaskedStatisticsFilter( image2); testclassInstance.TestofSkewnessKurtosisAndMPPForUnmaskedImages(mitkFilter2, -1.1547, 2.33333, 2.75); testclassInstance.TestofEntropyUniformityAndUppForUnmaskedImages( mitkFilter2, 0.811278, 0.625, 0.625); MITK_TEST_END() } diff --git a/Modules/ImageStatistics/files.cmake b/Modules/ImageStatistics/files.cmake index 454e2941ce..23d6c94783 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 - mitkExtendedStatisticsImageFilter.h + mitkStatisticsImageFilter.h mitkExtendedLabelStatisticsImageFilter.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 -) \ No newline at end of file +) diff --git a/Modules/ImageStatistics/mitkExtendedLabelStatisticsImageFilter.h b/Modules/ImageStatistics/mitkExtendedLabelStatisticsImageFilter.h index 37352663b7..23081fa423 100644 --- a/Modules/ImageStatistics/mitkExtendedLabelStatisticsImageFilter.h +++ b/Modules/ImageStatistics/mitkExtendedLabelStatisticsImageFilter.h @@ -1,347 +1,345 @@ /*============================================================================ 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 itksys::hash_map< LabelPixelType, LabelStatistics > MapType; - typedef typename itksys::hash_map< LabelPixelType, LabelStatistics >::const_iterator StatisticsMapConstIterator; - typedef typename itksys::hash_map< LabelPixelType, LabelStatistics >::iterator StatisticsMapIterator; - typedef IdentifierType MapSizeType; + 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 AfterThreadedGenerateData() override; + void AfterStreamedGenerateData() override; /** Initialize some accumulators before the threads run. */ - void BeforeThreadedGenerateData() override; + void BeforeStreamedGenerateData() override; /** Multi-thread version GenerateData. */ - void ThreadedGenerateData(const typename TInputImage::RegionType & - outputRegionForThread, - ThreadIdType threadId) override; + 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 index 05223f2105..f2ecd243d6 100644 --- a/Modules/ImageStatistics/mitkExtendedLabelStatisticsImageFilter.hxx +++ b/Modules/ImageStatistics/mitkExtendedLabelStatisticsImageFilter.hxx @@ -1,737 +1,734 @@ /*============================================================================ 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 > - ::BeforeThreadedGenerateData() + ::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< typename TInputImage, typename TLabelImage > - void - ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > - ::ThreadedGenerateData(const typename TInputImage::RegionType & outputRegionForThread, - ThreadIdType threadId) + 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 >:: - AfterThreadedGenerateData() + 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/mitkExtendedStatisticsImageFilter.h b/Modules/ImageStatistics/mitkExtendedStatisticsImageFilter.h deleted file mode 100644 index 00236bce16..0000000000 --- a/Modules/ImageStatistics/mitkExtendedStatisticsImageFilter.h +++ /dev/null @@ -1,207 +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 __mitkExtendedStatisticsImageFilter -#define __mitkExtendedStatisticsImageFilter - -#include "itkStatisticsImageFilter.h" -#include -#include - -namespace itk -{ - /** - * \class ExtendedStatisticsImageFilter - * \brief Extension of the itkStatisticsImageFilter that also calculates the Skewness and Kurtosis. - * - * This class inherits from the itkStatisticsImageFilter and - * uses its results for the calculation of other additional coefficients: - * the Skewness and Kurtosis. - * - * As the StatisticsImageFilter stores the statistics in the outputs 1 to 6 by the - * StatisticsImageFilter, the skewness, kurtosis, MPP, UPP, Uniformity, Entropy and Median are stored in the outputs - * 7 to 14 by this filter. - */ - template< class TInputImage > - class ExtendedStatisticsImageFilter : public StatisticsImageFilter< TInputImage > - { - public: - /** Standard Self typedef */ - typedef ExtendedStatisticsImageFilter Self; - typedef StatisticsImageFilter< TInputImage > Superclass; - typedef SmartPointer< Self > Pointer; - typedef SmartPointer< const Self > ConstPointer; - typedef typename Superclass::RealType RealType; - typedef typename Superclass::RealObjectType RealObjectType; - typedef typename Superclass::PixelType PixelType; - - /** Histogram-related typedefs */ - typedef itk::Statistics::Histogram< RealType > HistogramType; - typedef typename HistogramType::Pointer HistogramPointer; - - itkFactorylessNewMacro( Self ); - itkCloneMacro( Self ); - itkTypeMacro( ExtendedStatisticsImageFilter, StatisticsImageFilter ); - - /** - * \brief Return the computed Skewness. - */ - double GetSkewness() const - { - return this->GetSkewnessOutput()->Get(); - } - - /** - * \brief Return the computed Median - */ - double GetMedian() const - { - return this->GetMedianOutput()->Get(); - } - - /** - * \brief Return the computed Kurtosis. - */ - double GetKurtosis() const - { - return this->GetKurtosisOutput()->Get(); - } - - /* \brief Return the computed MPP. - */ - double GetMPP() const - { - return this->GetMPPOutput()->Get(); - } - - /** - * \brief Return the computed Uniformity. - */ - double GetUniformity() const - { - return this->GetUniformityOutput()->Get(); - } - - /** - *\brief Return the computed Entropy. - */ - double GetEntropy() const - { - return this->GetEntropyOutput()->Get(); - } - - /** - * \brief Return the computed UPP. - */ - double GetUPP() const - { - return this->GetUPPOutput()->Get(); - } - - - /** - * \brief Return the computed Histogram. - */ - const typename HistogramType::Pointer - GetHistogram() - { - if (m_HistogramCalculated) - { - return m_Histogram; - } - else - { - return nullptr; - } - } - - - /** specify Histogram parameters */ - void SetHistogramParameters(const int numBins, RealType lowerBound, - RealType upperBound); - - protected: - - ExtendedStatisticsImageFilter(); - - ~ExtendedStatisticsImageFilter() override{}; - - void BeforeThreadedGenerateData() override; - - /** Multi-thread version GenerateData. */ - void ThreadedGenerateData(const typename StatisticsImageFilter::RegionType & - outputRegionForThread, - ThreadIdType threadId) override; - - /** - * brief Calls AfterThreadedGenerateData() of the superclass and the main methods - */ - void AfterThreadedGenerateData() override; - - - RealObjectType* GetSkewnessOutput(); - - const RealObjectType* GetSkewnessOutput() const; - - RealObjectType* GetKurtosisOutput(); - - const RealObjectType* GetKurtosisOutput() const; - - RealObjectType* GetMPPOutput(); - - const RealObjectType* GetMPPOutput() const; - - RealObjectType* GetEntropyOutput(); - - const RealObjectType* GetEntropyOutput() const; - - RealObjectType* GetUniformityOutput(); - - const RealObjectType* GetUniformityOutput() const; - - RealObjectType* GetUPPOutput(); - - const RealObjectType* GetUPPOutput() const; - - RealObjectType* GetMedianOutput(); - - const RealObjectType* GetMedianOutput() const; - - using Superclass::MakeOutput; - DataObject::Pointer MakeOutput( ProcessObject::DataObjectPointerArraySizeType idx ) override; - -private: - Array< RealType > m_ThreadSum; - Array< RealType > m_SumOfSquares; - Array< RealType > m_SumOfCubes; - Array< RealType > m_SumOfQuadruples; - Array< SizeValueType > m_Count; - Array< SizeValueType > m_PositivePixelCount; - Array< RealType > m_ThreadSumOfPositivePixels; - Array< PixelType > m_ThreadMin; - Array< PixelType > m_ThreadMax; - std::vector< HistogramPointer > m_HistogramPerThread; - HistogramPointer m_Histogram; - bool m_UseHistogram; - bool m_HistogramCalculated; - RealType m_LowerBound, m_UpperBound; - int m_NumBins; - - - }; // end of class - -} // end namespace itk - -#ifndef ITK_MANUAL_INSTANTIATION -#include "mitkExtendedStatisticsImageFilter.hxx" -#endif - -#endif diff --git a/Modules/ImageStatistics/mitkExtendedStatisticsImageFilter.hxx b/Modules/ImageStatistics/mitkExtendedStatisticsImageFilter.hxx deleted file mode 100644 index 7ade1af4e9..0000000000 --- a/Modules/ImageStatistics/mitkExtendedStatisticsImageFilter.hxx +++ /dev/null @@ -1,470 +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 __mitkExtendedStatisticsImageFilter_hxx -#define __mitkExtendedStatisticsImageFilter_hxx - -#include "mitkExtendedStatisticsImageFilter.h" -#include - - -namespace itk -{ - template< class TInputImage > - ExtendedStatisticsImageFilter< TInputImage >::ExtendedStatisticsImageFilter() - : StatisticsImageFilter< TInputImage >(), - m_ThreadSum(1), - m_SumOfSquares(1), - m_SumOfCubes(1), - m_SumOfQuadruples(1), - m_Count(1), - m_PositivePixelCount(1), - m_ThreadSumOfPositivePixels(1), - m_ThreadMin(1), - m_ThreadMax(1), - m_UseHistogram(false), - m_HistogramCalculated(false) - { - /* - * 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 ); - - m_Histogram = ITK_NULLPTR; - m_HistogramPerThread.resize(0); - } - - 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(); - } - case 13: - { - return RealObjectType::New().GetPointer(); - } - default: - { - // might as well make an image - return Superclass::MakeOutput( output ); - } - } - } - - 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< typename TInputImage > - void - ExtendedStatisticsImageFilter< TInputImage > - ::BeforeThreadedGenerateData() - { - ThreadIdType numberOfThreads = this->GetNumberOfThreads(); - - if (m_UseHistogram) - { - // 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] = m_NumBins; - lb[0] = m_LowerBound; - ub[0] = m_UpperBound; - m_Histogram->Initialize(hsize, lb, ub); - - m_HistogramPerThread.resize(numberOfThreads); - - for (unsigned int i = 0; i < numberOfThreads; i++) - { - typename HistogramType::Pointer hist = HistogramType::New(); - typename HistogramType::SizeType hsize; - typename HistogramType::MeasurementVectorType lb; - typename HistogramType::MeasurementVectorType ub; - hsize.SetSize(1); - lb.SetSize(1); - ub.SetSize(1); - hist->SetMeasurementVectorSize(1); - hsize[0] = m_NumBins; - lb[0] = m_LowerBound; - ub[0] = m_UpperBound; - hist->Initialize(hsize, lb, ub); - m_HistogramPerThread[i] = hist; - } - } - - // Resize the thread temporaries - m_Count.SetSize(numberOfThreads); - m_SumOfSquares.SetSize(numberOfThreads); - m_SumOfCubes.SetSize(numberOfThreads); - m_SumOfQuadruples.SetSize(numberOfThreads); - m_ThreadSum.SetSize(numberOfThreads); - m_ThreadMin.SetSize(numberOfThreads); - m_ThreadMax.SetSize(numberOfThreads); - m_PositivePixelCount.SetSize(numberOfThreads); - m_ThreadSumOfPositivePixels.SetSize(numberOfThreads); - - // Initialize the temporaries - m_Count.Fill(NumericTraits< SizeValueType >::Zero); - m_ThreadSum.Fill(NumericTraits< RealType >::Zero); - m_SumOfSquares.Fill(NumericTraits< RealType >::Zero); - m_SumOfCubes.Fill(NumericTraits< RealType >::Zero); - m_SumOfQuadruples.Fill(NumericTraits< RealType >::Zero); - m_ThreadMin.Fill( NumericTraits< PixelType >::max() ); - m_ThreadMax.Fill( NumericTraits< PixelType >::NonpositiveMin() ); - m_PositivePixelCount.Fill(NumericTraits< SizeValueType >::Zero); - m_ThreadSumOfPositivePixels.Fill(NumericTraits< SizeValueType >::Zero); - } - - - template< typename TInputImage > - void - ExtendedStatisticsImageFilter< TInputImage > - ::ThreadedGenerateData(const typename StatisticsImageFilter::RegionType & - outputRegionForThread, - ThreadIdType threadId) - { - const SizeValueType size0 = outputRegionForThread.GetSize(0); - if( size0 == 0) - { - return; - } - RealType realValue; - PixelType value; - - typename HistogramType::IndexType histogramIndex(1); - typename HistogramType::MeasurementVectorType histogramMeasurement(1); - - RealType sum = NumericTraits< RealType >::ZeroValue(); - RealType sumOfPositivePixels = NumericTraits< RealType >::ZeroValue(); - RealType sumOfSquares = NumericTraits< RealType >::ZeroValue(); - RealType sumOfCubes = NumericTraits< RealType >::ZeroValue(); - RealType sumOfQuadruples = NumericTraits< RealType >::ZeroValue(); - SizeValueType count = NumericTraits< SizeValueType >::ZeroValue(); - SizeValueType countOfPositivePixels = NumericTraits< SizeValueType >::ZeroValue(); - PixelType min = NumericTraits< PixelType >::max(); - PixelType max = NumericTraits< PixelType >::NonpositiveMin(); - - ImageScanlineConstIterator< TInputImage > it (this->GetInput(), outputRegionForThread); - - // support progress methods/callbacks - const size_t numberOfLinesToProcess = outputRegionForThread.GetNumberOfPixels() / size0; - ProgressReporter progress( this, threadId, numberOfLinesToProcess ); - - // do the work - while ( !it.IsAtEnd() ) - { - while ( !it.IsAtEndOfLine() ) - { - value = it.Get(); - realValue = static_cast< RealType >( value ); - - if (m_UseHistogram) - { - histogramMeasurement[0] = value; - m_HistogramPerThread[threadId]->GetIndex(histogramMeasurement, histogramIndex); - m_HistogramPerThread[threadId]->IncreaseFrequencyOfIndex(histogramIndex, 1); - } - - if ( value < min ) - { - min = value; - } - if ( value > max ) - { - max = value; - } - if (value > 0) - { - sumOfPositivePixels += realValue; - ++countOfPositivePixels; - } - - sum += realValue; - sumOfSquares += ( realValue * realValue ); - sumOfCubes += std::pow(realValue, 3.); - sumOfQuadruples += std::pow(realValue, 4.); - ++count; - ++it; - } - it.NextLine(); - progress.CompletedPixel(); - } - - m_ThreadSum[threadId] = sum; - m_SumOfSquares[threadId] = sumOfSquares; - m_Count[threadId] = count; - m_ThreadMin[threadId] = min; - m_ThreadMax[threadId] = max; - m_ThreadSumOfPositivePixels[threadId] = sumOfPositivePixels; - m_PositivePixelCount[threadId] = countOfPositivePixels; - m_SumOfCubes[threadId] = sumOfCubes; - m_SumOfQuadruples[threadId] = sumOfQuadruples; - } - - - template< class TInputImage > - void - ExtendedStatisticsImageFilter< TInputImage > - ::AfterThreadedGenerateData() - { - ThreadIdType i; - SizeValueType count = 0; - SizeValueType countOfPositivePixels = 0; - RealType sumOfSquares = 0; - RealType sumOfCubes = 0; - RealType sumOfQuadruples = 0; - RealType sumOfPositivePixels = 0; - - ThreadIdType numberOfThreads = this->GetNumberOfThreads(); - - PixelType minimum; - PixelType maximum; - RealType mean; - RealType meanOfPositivePixels; - RealType sigma; - RealType variance; - RealType skewness; - RealType kurtosis; - RealType sum; - - sum = sumOfSquares = sumOfCubes = sumOfQuadruples = NumericTraits< RealType >::ZeroValue(); - count = countOfPositivePixels = 0; - - // Find the min/max over all threads and accumulate count, sum and - // sum of squares - minimum = NumericTraits< PixelType >::max(); - maximum = NumericTraits< PixelType >::NonpositiveMin(); - for ( i = 0; i < numberOfThreads; i++ ) - { - count += m_Count[i]; - sum += m_ThreadSum[i]; - sumOfSquares += m_SumOfSquares[i]; - sumOfCubes += m_SumOfCubes[i]; - sumOfQuadruples += m_SumOfQuadruples[i]; - sumOfPositivePixels += m_ThreadSumOfPositivePixels[i]; - countOfPositivePixels += m_PositivePixelCount[i]; - - if ( m_ThreadMin[i] < minimum ) - { - minimum = m_ThreadMin[i]; - } - if ( m_ThreadMax[i] > maximum ) - { - maximum = m_ThreadMax[i]; - } - // if enabled, update the histogram for this label - if ( m_UseHistogram ) - { - typename HistogramType::IndexType index; - index.SetSize(1); - for ( int bin = 0; bin < m_NumBins; ++bin ) - { - index[0] = bin; - m_Histogram->IncreaseFrequency( bin, m_HistogramPerThread[i]->GetFrequency(bin) ); - } - } - } - // compute statistics - mean = sum / static_cast< RealType >( count ); - meanOfPositivePixels = sumOfPositivePixels / static_cast< RealType >( countOfPositivePixels ); - - // unbiased estimate - variance = ( sumOfSquares - ( sum * sum / static_cast< RealType >( count ) ) ) - / ( static_cast< RealType >( count ) ); - RealType secondMoment, thirdMoment, fourthMoment; - secondMoment = sumOfSquares / static_cast< RealType >( count ); - thirdMoment = sumOfCubes / static_cast< RealType >( count ); - fourthMoment = sumOfQuadruples / static_cast< RealType >( count ); - - skewness = (thirdMoment - 3. * secondMoment * mean + 2. * std::pow(mean, 3.)) / std::pow(secondMoment - std::pow(mean, 2.), 1.5); // see http://www.boost.org/doc/libs/1_51_0/doc/html/boost/accumulators/impl/skewness_impl.html - kurtosis = (fourthMoment - 4. * thirdMoment * mean + 6. * secondMoment * std::pow(mean, 2.) - 3. * std::pow(mean, 4.)) / std::pow(secondMoment - std::pow(mean, 2.), 2.); // see http://www.boost.org/doc/libs/1_46_1/doc/html/boost/accumulators/impl/kurtosis_impl.html, dropped -3 - sigma = std::sqrt(variance); - - // Set the outputs - this->GetMinimumOutput()->Set(minimum); - this->GetMaximumOutput()->Set(maximum); - this->GetMeanOutput()->Set(mean); - this->GetSigmaOutput()->Set(sigma); - this->GetVarianceOutput()->Set(variance); - this->GetSumOutput()->Set(sum); - this->GetKurtosisOutput()->Set(kurtosis); - this->GetSkewnessOutput()->Set(skewness); - this->GetMPPOutput()->Set(meanOfPositivePixels); - - if (m_UseHistogram) - { - mitk::HistogramStatisticsCalculator histStatCalc; - histStatCalc.SetHistogram(m_Histogram); - histStatCalc.CalculateStatistics(); - this->GetUniformityOutput()->Set(histStatCalc.GetUniformity()); - this->GetUPPOutput()->Set(histStatCalc.GetUPP()); - this->GetEntropyOutput()->Set(histStatCalc.GetEntropy()); - this->GetMedianOutput()->Set(histStatCalc.GetMedian()); - - } - m_HistogramCalculated = true; - } - - template< typename TInputImage > - void - ExtendedStatisticsImageFilter< TInputImage > - ::SetHistogramParameters(const int numBins, RealType lowerBound, RealType upperBound) - { - m_NumBins = numBins; - m_LowerBound = lowerBound; - m_UpperBound = upperBound; - m_UseHistogram = true; - } - -} - -#endif diff --git a/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp b/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp index bfa5041c2f..723e7702db 100644 --- a/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp +++ b/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp @@ -1,552 +1,552 @@ /*============================================================================ 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 itk::ExtendedStatisticsImageFilter ImageStatisticsFilterType; + 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().GetPointer(); 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 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; std::vector relevantLabels = minMaxFilter->GetRelevantLabels(); MapType minVals; MapType maxVals; std::map nBins; for (LabelPixelType label : relevantLabels) { minVals.insert(PairType(label, minMaxFilter->GetMin(label))); maxVals.insert(PairType(label, 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)); } 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->Update(); std::list labels = imageStatisticsFilter->GetRelevantLabels(); 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::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(); 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/mitkStatisticsImageFilter.h b/Modules/ImageStatistics/mitkStatisticsImageFilter.h new file mode 100644 index 0000000000..22c3db0225 --- /dev/null +++ b/Modules/ImageStatistics/mitkStatisticsImageFilter.h @@ -0,0 +1,181 @@ +/*============================================================================ + +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 + +namespace mitk +{ + template + class StatisticsImageFilter : public itk::ImageSink + { + public: + /** Standard Self type alias */ + using Self = StatisticsImageFilter; + using Superclass = itk::ImageSink; + using Pointer = itk::SmartPointer; + using ConstPointer = itk::SmartPointer; + + /** Method for creation bypassing the object factory. */ + itkFactorylessNewMacro(Self); + + /** Runtime information support. */ + itkTypeMacro(StatisticsImageFilter, itk::ImageSink); + + /** Image related type alias. */ + 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; + + /** Image related type alias. */ + static constexpr unsigned int ImageDimension = TInputImage::ImageDimension; + + /** Type to use for computations. */ + using RealType = typename itk::NumericTraits::RealType; + + /** Smart Pointer type to a DataObject. */ + using DataObjectPointer = typename itk::DataObject::Pointer; + + template + using SimpleDataObjectDecorator = itk::SimpleDataObjectDecorator; + + /** Type of DataObjects used for scalar outputs */ + 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 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); + + 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() override = default; + + 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(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: + itk::CompensatedSummation m_ThreadSum; + 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_ThreadMin; + PixelType m_ThreadMax; + + std::mutex m_Mutex; + }; +} + +#ifndef ITK_MANUAL_INSTANTIATION +#include +#endif + +#endif diff --git a/Modules/ImageStatistics/mitkStatisticsImageFilter.hxx b/Modules/ImageStatistics/mitkStatisticsImageFilter.hxx new file mode 100644 index 0000000000..0b97f4a8c7 --- /dev/null +++ b/Modules/ImageStatistics/mitkStatisticsImageFilter.hxx @@ -0,0 +1,214 @@ +/*============================================================================ + +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 + +template +mitk::StatisticsImageFilter::StatisticsImageFilter() + : m_ThreadSum(1), + m_SumOfPositivePixels(1), + m_SumOfSquares(1), + m_SumOfCubes(1), + m_SumOfQuadruples(1), + m_Count(1), + m_CountOfPositivePixels(1), + m_ThreadMin(1), + m_ThreadMax(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 +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(); + } + + return Superclass::MakeOutput(name); +} + +template +void mitk::StatisticsImageFilter::BeforeStreamedGenerateData() +{ + Superclass::BeforeStreamedGenerateData(); + + m_ThreadSum = 0; + m_SumOfPositivePixels = 0; + m_SumOfSquares = 0; + m_SumOfCubes = 0; + m_SumOfQuadruples = 0; + m_Count = 0; + m_CountOfPositivePixels = 0; + m_ThreadMin = itk::NumericTraits::max(); + m_ThreadMax = itk::NumericTraits::NonpositiveMin(); +} + +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 squareValue = 0; + + itk::ImageScanlineConstIterator it(this->GetInput(), regionForThread); + + while (!it.IsAtEnd()) + { + while (!it.IsAtEndOfLine()) + { + const auto& value = it.Get(); + const auto realValue = static_cast(value); + 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); + + m_ThreadSum += sum; + m_SumOfPositivePixels += sumOfPositivePixels; + m_SumOfSquares += sumOfSquares; + m_SumOfCubes += sumOfCubes; + m_SumOfQuadruples += sumOfQuadruples; + m_Count += count; + m_CountOfPositivePixels += countOfPositivePixels; + m_ThreadMin = std::min(min, m_ThreadMin); + m_ThreadMax = std::max(max, m_ThreadMax); +} + +template +void mitk::StatisticsImageFilter::AfterStreamedGenerateData() +{ + Superclass::AfterStreamedGenerateData(); + + const RealType sum = m_ThreadSum; + const RealType sumOfPositivePixels = m_SumOfPositivePixels; + const RealType sumOfSquares = m_SumOfSquares; + const RealType sumOfCubes = m_SumOfCubes; + const RealType sumOfQuadruples = m_SumOfQuadruples; + const itk::SizeValueType count = m_Count; + const itk::SizeValueType countOfPositivePixels = m_CountOfPositivePixels; + const PixelType minimum = m_ThreadMin; + const PixelType maximum = m_ThreadMax; + + 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); + this->SetEntropy(entropy); + this->SetUniformity(uniformity); + this->SetUPP(upp); + this->SetMedian(median); +} + +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