diff --git a/Modules/DiffusionImaging/MiniApps/ImageStatisticsMiniapp_3.cpp b/Modules/DiffusionImaging/MiniApps/ImageStatisticsMiniapp_3.cpp index 9ce6f8a173..0e42234a9a 100644 --- a/Modules/DiffusionImaging/MiniApps/ImageStatisticsMiniapp_3.cpp +++ b/Modules/DiffusionImaging/MiniApps/ImageStatisticsMiniapp_3.cpp @@ -1,82 +1,71 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkCommandLineParser.h" #include "mitkImage.h" #include #include #include #include #include #include "mitkIOUtil.h" #include #include #include #include #include "mitkImageAccessByItk.h" #include #include #include #include #include #include #include #include #include #include - - -template -void printMap(std::map input) -{ - for (auto it = input.begin(); it != input.end(); ++it) - { - std::cout << it->first<< ": " << it->second<< std::endl; - } - std::cout << std::endl; -} - +#include int main( int argc, char* argv[] ) { unsigned int timeStep = 0; std::string inputImageFile; inputImageFile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic2DplusT.nrrd"; // Load image mitk::Image::Pointer inputImage = mitk::IOUtil::LoadImage(inputImageFile); // Calculate statistics mitk::ImageStatisticsCalculator_v2::Pointer calculator = mitk::ImageStatisticsCalculator_v2::New(); std::cout << "calculating statistics (unmasked) itk: " << std::endl; mitk::ImageStatisticsCalculator_v2::StatisticsContainer::Pointer result; calculator->SetInputImage(inputImage); calculator->SetNBinsForHistogramStatistics(100); for (unsigned int i=0; i < inputImage->GetTimeSteps(); i++) { std::cout << "Results for time step " << i << ":" << std::endl; result = calculator->GetStatistics(i, 1); result->PrintSelf(); std::cout << std::endl; } return EXIT_SUCCESS; } diff --git a/Modules/ImageStatistics/mitkExtendedLabelStatisticsImageFilter_2.h b/Modules/ImageStatistics/mitkExtendedLabelStatisticsImageFilter_2.h index 5931cbc24f..ee7b7cd92a 100644 --- a/Modules/ImageStatistics/mitkExtendedLabelStatisticsImageFilter_2.h +++ b/Modules/ImageStatistics/mitkExtendedLabelStatisticsImageFilter_2.h @@ -1,292 +1,293 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef __mitkExtendedLabelStatisticsImageFilter_h2 #define __mitkExtendedLabelStatisticsImageFilter_h2 #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 six additional coefficients: * the Skewness,Kurtosis,Uniformity,UPP,MPP,Entropy * * As these coefficient are based on the mean and the sigma which are both calculated * by the LabelStatisticsImageFilter, the method AfterThreadedGenerateData() is overwritten * and calls ComputeSkewnessKurtosisAndMPP() and ComputeEntropyUniformityAndUPP after the AfterThreadedGenerateData() * while the second coefficient Method is only called when the the method useHistogram is on!!! * implementation of the superclass is called. * */ template< class TInputImage, class TLabelImage > class ExtendedLabelStatisticsImageFilter2 : public LabelStatisticsImageFilter< TInputImage, TLabelImage > { public: typedef ExtendedLabelStatisticsImageFilter2 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 itk::Statistics::Histogram HistogramType; itkFactorylessNewMacro( Self ); itkCloneMacro( Self ); itkTypeMacro(ExtendedLabelStatisticsImageFilter2, 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 = ITK_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 MapConstIterator; typedef IdentifierType MapSizeType; /** Type of the container used to store valid label values */ typedef std::vector ValidLabelValuesContainerType; - /*getter method for the new coefficients*/ + /*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; protected: ExtendedLabelStatisticsImageFilter2(); virtual ~ExtendedLabelStatisticsImageFilter2(){} void AfterThreadedGenerateData(); /** Initialize some accumulators before the threads run. */ void BeforeThreadedGenerateData(); /** Multi-thread version GenerateData. */ void ThreadedGenerateData(const typename TInputImage::RegionType & outputRegionForThread, typename TInputImage::ThreadIdType threadId); private: std::vector< MapType > m_LabelStatisticsPerThread; MapType m_LabelStatistics; ValidLabelValuesContainerType m_ValidLabelValues; bool m_UseHistograms; typename HistogramType::SizeType m_NumBins; RealType m_LowerBound; RealType m_UpperBound; bool m_MaskNonEmpty; }; // end of class } // end namespace itk #ifndef ITK_MANUAL_INSTANTIATION #include "mitkExtendedLabelStatisticsImageFilter_2.hxx" #endif #endif diff --git a/Modules/ImageStatistics/mitkExtendedLabelStatisticsImageFilter_2.hxx b/Modules/ImageStatistics/mitkExtendedLabelStatisticsImageFilter_2.hxx index b05e02f46a..a38ce4d79d 100644 --- a/Modules/ImageStatistics/mitkExtendedLabelStatisticsImageFilter_2.hxx +++ b/Modules/ImageStatistics/mitkExtendedLabelStatisticsImageFilter_2.hxx @@ -1,444 +1,464 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef _mitkExtendedLabelStatisticsImageFilter2_hxx #define _mitkExtendedLabelStatisticsImageFilter2_hxx #include "mitkExtendedLabelStatisticsImageFilter_2.h" #include "itkImageRegionConstIteratorWithIndex.h" #include "itkImageRegionConstIterator.h" #include #include #include "mitkNumericConstants.h" #include "mitkLogMacros.h" #include namespace itk { template< class TInputImage, class TLabelImage > bool ExtendedLabelStatisticsImageFilter2< TInputImage, TLabelImage > ::GetMaskingNonEmpty() const { return m_MaskNonEmpty; } template< class TInputImage, class TLabelImage > typename ExtendedLabelStatisticsImageFilter2< TInputImage, TLabelImage >::RealType ExtendedLabelStatisticsImageFilter2< TInputImage, TLabelImage > ::GetUniformity(LabelPixelType label) const { MapIterator mapIt; mapIt = m_LabelStatistics.find(label); if ( mapIt == m_LabelStatistics.end() ) { // label does not exist, return a default value return NumericTraits< PixelType >::Zero; } else { return ( *mapIt ).second.m_Uniformity; } } + + template< class TInputImage, class TLabelImage > + typename ExtendedLabelStatisticsImageFilter2< TInputImage, TLabelImage >::RealType + ExtendedLabelStatisticsImageFilter2< TInputImage, TLabelImage > + ::GetMedian(LabelPixelType label) const + { + MapIterator mapIt; + + mapIt = m_LabelStatistics.find(label); + if ( mapIt == m_LabelStatistics.end() ) + { + // label does not exist, return a default value + return NumericTraits< PixelType >::Zero; + } + else + { + return ( *mapIt ).second.m_Median; + } + } + template< class TInputImage, class TLabelImage > typename ExtendedLabelStatisticsImageFilter2< TInputImage, TLabelImage >::RealType ExtendedLabelStatisticsImageFilter2< TInputImage, TLabelImage > ::GetEntropy(LabelPixelType label) const { MapIterator mapIt; mapIt = m_LabelStatistics.find(label); if ( mapIt == m_LabelStatistics.end() ) { // label does not exist, return a default value return NumericTraits< PixelType >::Zero; } else { return ( *mapIt ).second.m_Entropy; } } template< class TInputImage, class TLabelImage > typename ExtendedLabelStatisticsImageFilter2< TInputImage, TLabelImage >::RealType ExtendedLabelStatisticsImageFilter2< TInputImage, TLabelImage > ::GetUPP(LabelPixelType label) const { MapIterator mapIt; mapIt = m_LabelStatistics.find(label); if ( mapIt == m_LabelStatistics.end() ) { // label does not exist, return a default value return NumericTraits< PixelType >::Zero; } else { return ( *mapIt ).second.m_UPP; } } template< class TInputImage, class TLabelImage > typename ExtendedLabelStatisticsImageFilter2< TInputImage, TLabelImage >::RealType ExtendedLabelStatisticsImageFilter2< TInputImage, TLabelImage > ::GetMPP(LabelPixelType label) const { MapIterator mapIt; mapIt = m_LabelStatistics.find(label); if ( mapIt == m_LabelStatistics.end() ) { // label does not exist, return a default value return NumericTraits< PixelType >::Zero; } else { return ( *mapIt ).second.m_MPP; } } template< class TInputImage, class TLabelImage > typename ExtendedLabelStatisticsImageFilter2< TInputImage, TLabelImage >::RealType ExtendedLabelStatisticsImageFilter2< TInputImage, TLabelImage > ::GetKurtosis(LabelPixelType label) const { MapIterator mapIt; mapIt = m_LabelStatistics.find(label); if ( mapIt == m_LabelStatistics.end() ) { // label does not exist, return a default value return NumericTraits< PixelType >::Zero; } else { return ( *mapIt ).second.m_Kurtosis; } } template< class TInputImage, class TLabelImage > typename ExtendedLabelStatisticsImageFilter2< TInputImage, TLabelImage >::RealType ExtendedLabelStatisticsImageFilter2< TInputImage, TLabelImage > ::GetSkewness(LabelPixelType label) const { MapIterator mapIt; mapIt = m_LabelStatistics.find(label); if ( mapIt == m_LabelStatistics.end() ) { // label does not exist, return a default value return NumericTraits< PixelType >::Zero; } else { return ( *mapIt ).second.m_Skewness; } } template< typename TInputImage, typename TLabelImage > void ExtendedLabelStatisticsImageFilter2< TInputImage, TLabelImage > ::BeforeThreadedGenerateData() { 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 > void ExtendedLabelStatisticsImageFilter2< TInputImage, TLabelImage > ::ThreadedGenerateData(const typename TInputImage::RegionType & outputRegionForThread, typename TInputImage::ThreadIdType threadId) { 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); MapIterator mapIt; // 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() ) { 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() ) { // create a new statistics object typedef typename MapType::value_type MapValueType; if ( m_UseHistograms ) { mapIt = m_LabelStatisticsPerThread[threadId].insert( MapValueType( label, LabelStatistics(m_NumBins[0], m_LowerBound, m_UpperBound) ) ).first; } 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_SumOfPisitivePixels += value; } // if enabled, update the histogram for this label if ( m_UseHistograms ) { 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 ExtendedLabelStatisticsImageFilter2< TInputImage, TLabelImage >:: AfterThreadedGenerateData() { MapIterator mapIt; MapConstIterator 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_UseHistograms ) { mapIt = m_LabelStatistics.insert( MapValueType( ( *threadIt ).first, LabelStatistics(m_NumBins[0], m_LowerBound, m_UpperBound) ) ).first; } 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_UseHistograms ) { typename HistogramType::IndexType index; index.SetSize(1); for ( unsigned int bin = 0; bin < m_NumBins[0]; 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.), 1.5) - 3.; // see http://www.boost.org/doc/libs/1_51_0/doc/html/boost/accumulators/impl/kurtosis_impl.html } 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 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