diff --git a/Modules/ImageStatistics/files.cmake b/Modules/ImageStatistics/files.cmake index ffca98b49e..c77d7d9ee7 100644 --- a/Modules/ImageStatistics/files.cmake +++ b/Modules/ImageStatistics/files.cmake @@ -1,33 +1,35 @@ set(CPP_FILES mitkImageStatisticsCalculator.cpp mitkImageStatisticsCalculator2.cpp mitkPointSetStatisticsCalculator.cpp mitkPointSetDifferenceStatisticsCalculator.cpp mitkIntensityProfile.cpp mitkHotspotMaskGenerator.cpp mitkMaskGenerator.cpp mitkPlanarFigureMaskGenerator.cpp mitkMultiLabelMaskGenerator.cpp mitkImageMaskGenerator.cpp mitkHistogramStatisticsCalculator.cpp mitkMaskUtilities.cpp mitkIgnorePixelMaskGenerator.cpp ) set(H_FILES mitkImageStatisticsCalculator.h mitkImageStatisticsCalculator2.h mitkPointSetDifferenceStatisticsCalculator.h mitkPointSetStatisticsCalculator.h mitkExtendedStatisticsImageFilter.h + mitkExtendedStatisticsImageFilter_2.h mitkExtendedLabelStatisticsImageFilter.h + mitkExtendedLabelStatisticsImageFilter_2.h mitkHotspotMaskGenerator.h mitkMaskGenerator.h mitkPlanarFigureMaskGenerator.h mitkMultiLabelMaskGenerator.h mitkImageMaskGenerator.h mitkHistogramStatisticsCalculator.h mitkMaskUtilities.h mitkitkMaskImageFilter.h mitkIgnorePixelMaskGenerator.h ) diff --git a/Modules/ImageStatistics/mitkExtendedLabelStatisticsImageFilter_2.h b/Modules/ImageStatistics/mitkExtendedLabelStatisticsImageFilter_2.h new file mode 100644 index 0000000000..5931cbc24f --- /dev/null +++ b/Modules/ImageStatistics/mitkExtendedLabelStatisticsImageFilter_2.h @@ -0,0 +1,292 @@ +/*=================================================================== + +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*/ + RealType GetSkewness(LabelPixelType label) const; + RealType GetKurtosis(LabelPixelType label) const; + RealType GetUniformity( 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 new file mode 100644 index 0000000000..b05e02f46a --- /dev/null +++ b/Modules/ImageStatistics/mitkExtendedLabelStatisticsImageFilter_2.hxx @@ -0,0 +1,444 @@ +/*=================================================================== + +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 > + ::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 diff --git a/Modules/ImageStatistics/mitkExtendedStatisticsImageFilter_2.h b/Modules/ImageStatistics/mitkExtendedStatisticsImageFilter_2.h new file mode 100644 index 0000000000..50382ab31d --- /dev/null +++ b/Modules/ImageStatistics/mitkExtendedStatisticsImageFilter_2.h @@ -0,0 +1,216 @@ +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center, +Division of Medical and Biological Informatics. +All rights reserved. + +This software is distributed WITHOUT ANY WARRANTY; without +even the implied warranty of MERCHANTABILITY or FITNESS FOR +A PARTICULAR PURPOSE. + +See LICENSE.txt or http://www.mitk.org for details. + +===================================================================*/ +#ifndef __mitkExtendedStatisticsImageFilter_h2 +#define __mitkExtendedStatisticsImageFilter_h2 + +#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 these coefficient are based on the mean and the sigma which are both calculated + * by the StatisticsImageFilter, the method AfterThreadedGenerateData() is overwritten + * and calls ComputeSkewnessKurtosisAndMPP() and ComputeEntropyUniformityAndUPP + * after the AfterThreadedGenerateData() + * implementation of the superclass is called. + * + * As the StatisticsImageFilter stores the statistics in the outputs 1 to 6 by the + * StatisticsImageFilter, the skewness, kurtosis,MPP,UPP,Uniformity and Entropy are stored in the outputs + * 7 to 13 by this filter. + */ + template< class TInputImage > + class ExtendedStatisticsImageFilter2 : public StatisticsImageFilter< TInputImage > + { + public: + /** Standard Self typedef */ + typedef ExtendedStatisticsImageFilter2 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( ExtendedStatisticsImageFilter2, 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 ITK_NULLPTR; + } + } + + + /** specify Histogram parameters */ + void SetHistogramParameters(const int numBins, RealType lowerBound, + RealType upperBound); + + protected: + + ExtendedStatisticsImageFilter2(); + + virtual ~ExtendedStatisticsImageFilter2(){}; + + void BeforeThreadedGenerateData(); + + /** Multi-thread version GenerateData. */ + void ThreadedGenerateData(const typename StatisticsImageFilter::RegionType & + outputRegionForThread, + typename StatisticsImageFilter::ThreadIdType threadId); + + /** + * brief Calls AfterThreadedGenerateData() of the superclass and the main methods + */ + void AfterThreadedGenerateData(); + + + 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; + + virtual DataObject::Pointer MakeOutput( ProcessObject::DataObjectPointerArraySizeType idx ); + +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; + Array< HistogramPointer > m_HistogramPerThread; + HistogramPointer m_Histogram; + bool m_UseHistogram; + bool m_HistogramCalculated; + RealType m_LowerBound, m_UpperBound; + typename HistogramType::SizeType m_NumBins; + + + }; // end of class + +} // end namespace itk + +#ifndef ITK_MANUAL_INSTANTIATION +#include "mitkExtendedStatisticsImageFilter_2.hxx" +#endif + +#endif diff --git a/Modules/ImageStatistics/mitkExtendedStatisticsImageFilter_2.hxx b/Modules/ImageStatistics/mitkExtendedStatisticsImageFilter_2.hxx new file mode 100644 index 0000000000..b669e93de5 --- /dev/null +++ b/Modules/ImageStatistics/mitkExtendedStatisticsImageFilter_2.hxx @@ -0,0 +1,475 @@ +/*=================================================================== + +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 __mitkExtendedStatisticsImageFilter2_hxx +#define __mitkExtendedStatisticsImageFilter2_hxx + +#include "mitkExtendedStatisticsImageFilter_2.h" +#include + + +namespace itk +{ + template< class TInputImage > + ExtendedStatisticsImageFilter2< TInputImage >::ExtendedStatisticsImageFilter2() + : StatisticsImageFilter< TInputImage >(), + m_ThreadSum(1), + m_SumOfSquares(1), + m_SumOfCubes(1), + m_SumOfQuadruples(1), + m_Count(1), + m_ThreadMin(1), + m_ThreadMax(1), + m_PositivePixelCount(1), + m_ThreadSumOfPositivePixels(1), + m_HistogramPerThread(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 ); + } + + template< class TInputImage > + DataObject::Pointer + ExtendedStatisticsImageFilter2< TInputImage >::MakeOutput( ProcessObject::DataObjectPointerArraySizeType output) + { + switch ( output ) + { + case 7: + case 8: + case 9: + case 10: + case 11: + case 12: + { + return RealObjectType::New().GetPointer(); + break; + } + case 13: + { + return RealObjectType::New().GetPointer(); + break; + } + default: + { + // might as well make an image + return Superclass::MakeOutput( output ); + break; + } + } + } + + template< class TInputImage > + typename ExtendedStatisticsImageFilter2< TInputImage >::RealObjectType * + ExtendedStatisticsImageFilter2< TInputImage > + ::GetSkewnessOutput() + { + return static_cast< RealObjectType * >( this->ProcessObject::GetOutput(7) ); + } + + template< class TInputImage > + const typename ExtendedStatisticsImageFilter2< TInputImage >::RealObjectType * + ExtendedStatisticsImageFilter2< TInputImage > + ::GetSkewnessOutput() const + { + return static_cast< const RealObjectType * >( this->ProcessObject::GetOutput(7) ); + } + + + template< class TInputImage > + typename ExtendedStatisticsImageFilter2< TInputImage >::RealObjectType * + ExtendedStatisticsImageFilter2< TInputImage > + ::GetKurtosisOutput() + { + return static_cast< RealObjectType * >( this->ProcessObject::GetOutput(8) ); + } + + template< class TInputImage > + const typename ExtendedStatisticsImageFilter2< TInputImage >::RealObjectType * + ExtendedStatisticsImageFilter2< TInputImage > + ::GetKurtosisOutput() const + { + return static_cast< const RealObjectType * >( this->ProcessObject::GetOutput(8) ); + } + + template< class TInputImage > + typename ExtendedStatisticsImageFilter2< TInputImage >::RealObjectType * + ExtendedStatisticsImageFilter2< TInputImage > + ::GetUniformityOutput() + { + return static_cast< RealObjectType * >( this->ProcessObject::GetOutput(9) ); + } + + template< class TInputImage > + const typename ExtendedStatisticsImageFilter2< TInputImage >::RealObjectType * + ExtendedStatisticsImageFilter2< TInputImage > + ::GetUniformityOutput() const + { + return static_cast< const RealObjectType * >( this->ProcessObject::GetOutput(9) ); + } + + + template< class TInputImage > + typename ExtendedStatisticsImageFilter2< TInputImage >::RealObjectType * + ExtendedStatisticsImageFilter2< TInputImage > + ::GetEntropyOutput() + { + return static_cast< RealObjectType * >( this->ProcessObject::GetOutput(10) ); + } + + template< class TInputImage > + const typename ExtendedStatisticsImageFilter2< TInputImage >::RealObjectType * + ExtendedStatisticsImageFilter2< TInputImage > + ::GetEntropyOutput() const + { + return static_cast< const RealObjectType * >( this->ProcessObject::GetOutput(10) ); + } + + template< class TInputImage > + typename ExtendedStatisticsImageFilter2< TInputImage >::RealObjectType * + ExtendedStatisticsImageFilter2< TInputImage > + ::GetUPPOutput() + { + return static_cast< RealObjectType * >( this->ProcessObject::GetOutput(11) ); + } + + template< class TInputImage > + const typename ExtendedStatisticsImageFilter2< TInputImage >::RealObjectType * + ExtendedStatisticsImageFilter2< TInputImage > + ::GetUPPOutput() const + { + return static_cast< const RealObjectType * >( this->ProcessObject::GetOutput(11) ); + } + + template< class TInputImage > + typename ExtendedStatisticsImageFilter2< TInputImage >::RealObjectType * + ExtendedStatisticsImageFilter2< TInputImage > + ::GetMPPOutput() + { + return static_cast< RealObjectType * >( this->ProcessObject::GetOutput(12) ); + } + + template< class TInputImage > + const typename ExtendedStatisticsImageFilter2< TInputImage >::RealObjectType * + ExtendedStatisticsImageFilter2< TInputImage > + ::GetMPPOutput() const + { + return static_cast< const RealObjectType * >( this->ProcessObject::GetOutput(12) ); + } + + template< class TInputImage > + typename ExtendedStatisticsImageFilter2< TInputImage >::RealObjectType * + ExtendedStatisticsImageFilter2< TInputImage > + ::GetMedianOutput() + { + return static_cast< RealObjectType * >( this->ProcessObject::GetOutput(13) ); + } + + template< class TInputImage > + const typename ExtendedStatisticsImageFilter2< TInputImage >::RealObjectType * + ExtendedStatisticsImageFilter2< TInputImage > + ::GetMedianOutput() const + { + return static_cast< const RealObjectType * >( this->ProcessObject::GetOutput(13) ); + } + + + template< typename TInputImage > + void + ExtendedStatisticsImageFilter2< 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.set_size(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); + m_Histogram->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 + ExtendedStatisticsImageFilter2< TInputImage > + ::ThreadedGenerateData(const typename StatisticsImageFilter::RegionType & + outputRegionForThread, + typename StatisticsImageFilter::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 + ExtendedStatisticsImageFilter2< TInputImage > + ::AfterThreadedGenerateData() + { + ThreadIdType i; + SizeValueType count; + SizeValueType countOfPositivePixels; + RealType sumOfSquares; + RealType sumOfCubes; + RealType sumOfQuadruples; + RealType sumOfPositivePixels; + + 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 ( unsigned int bin = 0; bin < m_NumBins[0]; 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.) - 3; // see http://www.boost.org/doc/libs/1_46_1/doc/html/boost/accumulators/impl/kurtosis_impl.html + 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 + ExtendedStatisticsImageFilter2< TInputImage > + ::SetHistogramParameters(const int numBins, RealType lowerBound, RealType upperBound) + { + m_NumBins[0] = numBins; + m_LowerBound = lowerBound; + m_UpperBound = upperBound; + m_UseHistogram = true; + } + +} + +#endif diff --git a/Modules/ImageStatistics/mitkImageStatisticsCalculator2.cpp b/Modules/ImageStatistics/mitkImageStatisticsCalculator2.cpp index 17fdd8bf26..24dcf08723 100644 --- a/Modules/ImageStatistics/mitkImageStatisticsCalculator2.cpp +++ b/Modules/ImageStatistics/mitkImageStatisticsCalculator2.cpp @@ -1,516 +1,519 @@ #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include +#include #include +#include #include #include #include "itkImageFileWriter.h" namespace mitk { void ImageStatisticsCalculator_v2::SetInputImage(mitk::Image::Pointer image) { if (image != m_Image) { m_Image = image; m_StatisticsUpdateRequiredByTimeStep.resize(m_Image->GetTimeSteps()); m_StatisticsByTimeStep.resize(m_Image->GetTimeSteps()); this->SetAllStatisticsToUpdateRequired(); } } void ImageStatisticsCalculator_v2::SetMask(mitk::MaskGenerator::Pointer mask) { if (mask != m_MaskGenerator) { m_MaskGenerator = mask; this->SetAllStatisticsToUpdateRequired(); } } void ImageStatisticsCalculator_v2::SetNBinsForHistogramStatistics(unsigned int nBins) { if (nBins != m_nBinsForHistogramStatistics) { m_nBinsForHistogramStatistics = nBins; this->SetAllStatisticsToUpdateRequired(); } } unsigned int ImageStatisticsCalculator_v2::GetNBinsForHistogramStatistics() const { return m_nBinsForHistogramStatistics; } ImageStatisticsCalculator_v2::StatisticsContainer::Pointer ImageStatisticsCalculator_v2::GetStatistics(unsigned int timeStep, unsigned int label) { std::cout << "getting statistics..." << std::endl; if (timeStep >= m_StatisticsByTimeStep.size()) { throw std::runtime_error("invalid timeStep in ImageStatisticsCalculator_v2::GetStatistics"); } if (m_Image.IsNull()) { throw std::runtime_error("no image"); } if (m_MaskGenerator.IsNotNull()) { m_MaskGenerator->SetTimeStep(timeStep); m_InternalMask = m_MaskGenerator->GetMask(); } ImageTimeSelector::Pointer imgTimeSel = ImageTimeSelector::New(); imgTimeSel->SetInput(m_Image); imgTimeSel->SetTimeNr(timeStep); imgTimeSel->UpdateLargestPossibleRegion(); m_ImageTimeSlice = imgTimeSel->GetOutput(); if (m_StatisticsUpdateRequiredByTimeStep[timeStep]) { std::cout << "calculating statistics..." << std::endl; // Calculate statistics with/without mask if (m_MaskGenerator.IsNull()) { // 1) calculate statistics unmasked: // plug image into itkstatisticsimagefilter (will be replaced by my awesome filter later on) // retrieve statistics and save them std::cout << "unmasked" << std::endl; AccessByItk_1(m_ImageTimeSlice, InternalCalculateStatisticsUnmasked, timeStep) } else { std::cout << "masked" << std::endl; // 2) calculate statistics masked // extract mask image region // plug mask and image into itklabelstatisticsimagefilter AccessByItk_1(m_ImageTimeSlice, InternalCalculateStatisticsMasked, timeStep) } m_StatisticsUpdateRequiredByTimeStep[timeStep]=false; } for (std::vector::iterator it = m_StatisticsByTimeStep[timeStep].begin(); it != m_StatisticsByTimeStep[timeStep].end(); ++it) { StatisticsContainer::Pointer statCont = *it; if (statCont->GetLabel() == label) { return statCont; } } MITK_ERROR << "Invalid Label: " << label; } void ImageStatisticsCalculator_v2::SetAllStatisticsToUpdateRequired() { for (unsigned int i = 0; i < m_StatisticsUpdateRequiredByTimeStep.size(); i++) { this->SetStatsTimeStepToUpdateRequired(i); } } void ImageStatisticsCalculator_v2::SetStatsTimeStepToUpdateRequired(unsigned int timeStep) { if (timeStep >= m_StatisticsUpdateRequiredByTimeStep.size()) { throw std::runtime_error("invalid timeStep in ImageStatisticsCalculator_v2::SetStatsTimeStepToUpdateRequired"); } m_StatisticsUpdateRequiredByTimeStep[timeStep] = true; m_StatisticsByTimeStep[timeStep].resize(0); } template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator_v2::InternalCalculateStatisticsUnmasked( typename itk::Image< TPixel, VImageDimension >* image, unsigned int timeStep) { typedef typename itk::Image< TPixel, VImageDimension > ImageType; StatisticsContainer::Pointer statisticsResult = StatisticsContainer::New(); typename itk::ExtendedStatisticsImageFilter::Pointer statisticsFilter = itk::ExtendedStatisticsImageFilter::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 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]; } statisticsResult->SetMinIndex(minIndex); statisticsResult->SetMaxIndex(maxIndex); // TODO: REMOVE THIS, the histogram statistics should be a separate module statisticsFilter->SetBinSize(m_nBinsForHistogramStatistics); try { statisticsFilter->Update(); } catch (const itk::ExceptionObject& e) { mitkThrow() << "Image statistics calculation failed due to following ITK Exception: \n " << e.what(); } // no mask, therefore just one label = the whole image m_StatisticsByTimeStep[timeStep].resize(1); statisticsResult->SetLabel(1); statisticsResult->SetN(image->GetLargestPossibleRegion().GetNumberOfPixels()); statisticsResult->SetMean(statisticsFilter->GetMean()); statisticsResult->SetMin(statisticsFilter->GetMinimum()); statisticsResult->SetMax(statisticsFilter->GetMaximum()); statisticsResult->SetVariance(statisticsFilter->GetVariance()); statisticsResult->SetStd(statisticsFilter->GetSigma()); statisticsResult->SetSkewness(statisticsFilter->GetSkewness()); statisticsResult->SetKurtosis(statisticsFilter->GetKurtosis()); statisticsResult->SetRMS(std::sqrt(std::pow(statisticsFilter->GetMean(), 2.) + statisticsFilter->GetVariance())); // variance = sigma^2 statisticsResult->SetMPP(statisticsFilter->GetMPP()); // histogram statistics HistogramType::Pointer histogram = InternalCalculateHistogramUnmasked(image, statisticsResult->GetMin(), statisticsResult->GetMax()); statisticsResult->SetHistogram(histogram); HistogramStatisticsCalculator histStatCalc; histStatCalc.SetHistogram(histogram); histStatCalc.CalculateStatistics(); statisticsResult->SetEntropy(histStatCalc.GetEntropy()); statisticsResult->SetMedian(histStatCalc.GetMedian()); statisticsResult->SetUniformity(histStatCalc.GetUniformity()); statisticsResult->SetUPP(histStatCalc.GetUPP()); m_StatisticsByTimeStep[timeStep][0] = statisticsResult; } template < typename TPixel, unsigned int VImageDimension > ImageStatisticsCalculator_v2::HistogramType::Pointer ImageStatisticsCalculator_v2::InternalCalculateHistogramUnmasked( typename itk::Image< TPixel, VImageDimension >* image, double minVal, double maxVal) { typedef typename itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Statistics::ImageToHistogramFilter HistogramFilterType; typedef itk::Statistics::Histogram HistogramType; typedef typename HistogramFilterType::HistogramSizeType SizeType; typename HistogramFilterType::Pointer histogramFilter = HistogramFilterType::New(); histogramFilter->SetInput(image); SizeType histSize(1); histSize[0] = m_nBinsForHistogramStatistics; histogramFilter->SetHistogramSize(histSize); // is that correct?? typename HistogramFilterType::HistogramMeasurementVectorType lowerBound(1); typename HistogramFilterType::HistogramMeasurementVectorType upperBound(1); lowerBound[0] = minVal; upperBound[0] = maxVal; histogramFilter->SetHistogramBinMinimum(lowerBound); histogramFilter->SetHistogramBinMaximum(upperBound); histogramFilter->Update(); typename HistogramType::Pointer histogram = histogramFilter->GetOutput(); return histogram; } /* the whole min and max calculation thing is completely bad... need a new itk filter for that */ // template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator_v2::GetMinAndMaxValue( // typename itk::Image< TPixel, VImageDimension >* inputImage, // double &minVal, // double &maxVal // ) // { // typedef itk::Image< TPixel, VImageDimension > ImageType; // typedef itk::ImageRegionConstIteratorWithIndex ImageRegionIteratorType; // ImageRegionIteratorType imageIterator( inputImage, inputImage->GetLargestPossibleRegion() ); // maxVal = std::numeric_limits::min(); // minVal = std::numeric_limits::max(); // long counter(0); // double actualPixelValue(0); // for( imageIterator.GoToBegin(); !imageIterator.IsAtEnd(); ++imageIterator) // { // counter++; // actualPixelValue = imageIterator.Value(); // if(actualPixelValue > maxVal) // { // maxVal = actualPixelValue; // } // else if(actualPixelValue < minVal) // { // minVal = actualPixelValue; // } // } // } template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator_v2::InternalCalculateStatisticsMasked( typename itk::Image< TPixel, VImageDimension >* image, unsigned int timeStep) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< unsigned short, VImageDimension > MaskType; typedef itk::ExtendedLabelStatisticsImageFilter< ImageType, MaskType > ImageStatisticsFilterType; typedef MaskUtilities< TPixel, VImageDimension > MaskUtilType; typedef itk::MaskImageFilter2 MaskImageFilterType2; typedef itk::MaskImageFilter MaskImageFilterType; typedef itk::MinimumMaximumImageCalculator MinMaxCalculatorType; // maskImage has to have the same dimension as image typename MaskType::Pointer maskImage = MaskType::New(); maskImage = ImageToItkImage< unsigned short, VImageDimension >(m_InternalMask); 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 // double minVal, maxVal; // GetMinAndMaxValue(adaptedImage.GetPointer(), minVal, maxVal); /** Searching the min and max value is only done over the pixels that are nonzero in the mask. This bevahior differs from the 'old' * imagestatisticscalculator where min and max were looked for over the entire adaptedImage*/ // find min, max, minindex and maxindex // make sure to only look in the masked region, use a masker for this typename MaskImageFilterType::Pointer maskImageFilter = MaskImageFilterType::New(); maskImageFilter->SetInput1(adaptedImage); maskImageFilter->SetInput2(maskImage); maskImageFilter->SetCoordinateTolerance(0.001); maskImageFilter->SetDirectionTolerance(0.001); maskImageFilter->Update(); typename MinMaxCalculatorType::Pointer imgMinMaxFilter = MinMaxCalculatorType::New(); imgMinMaxFilter->SetImage(maskImageFilter->GetOutput()); imgMinMaxFilter->Compute(); double minVal, maxVal; minVal = imgMinMaxFilter->GetMinimum(); maxVal = imgMinMaxFilter->GetMaximum(); typename ImageStatisticsFilterType::Pointer imageStatisticsFilter = ImageStatisticsFilterType::New(); imageStatisticsFilter->SetDirectionTolerance(0.001); imageStatisticsFilter->SetCoordinateTolerance(0.001); imageStatisticsFilter->SetInput(adaptedImage); imageStatisticsFilter->SetLabelInput(maskImage); imageStatisticsFilter->SetHistogramParameters(m_nBinsForHistogramStatistics, floor(minVal), ceil(maxVal)); imageStatisticsFilter->Update(); std::list labels = imageStatisticsFilter->GetRelevantLabels(); std::list::iterator it = labels.begin(); m_StatisticsByTimeStep[timeStep].resize(0); while(it != labels.end()) { StatisticsContainer::Pointer statisticsResult = StatisticsContainer::New(); // find min, max, minindex and maxindex // make sure to only look in the masked region, use a masker for this typename MaskImageFilterType2::Pointer maskImageFilter = MaskImageFilterType2::New(); maskImageFilter->SetInput1(adaptedImage); maskImageFilter->SetInput2(maskImage); maskImageFilter->SetOutsideValue(0); maskImageFilter->SetMaskingValue(*it); maskImageFilter->SetCoordinateTolerance(0.001); maskImageFilter->SetDirectionTolerance(0.001); maskImageFilter->UpdateLargestPossibleRegion(); // std::string outfile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic3D_masked_image.nrrd"; // typedef itk::ImageFileWriter< ImageType > WriterType; // typename WriterType::Pointer writer = WriterType::New(); // writer->SetFileName(outfile); // writer->SetInput(maskImageFilter->GetOutput()); // writer->Update(); // TODO adapt minindex and maxindex if mask was 2D (careful: a mask can be 2d even if the itk::Image is 3D (3d image with a size of 1 in one dimension)) typename MinMaxCalculatorType::Pointer imgMinMaxFilter = MinMaxCalculatorType::New(); imgMinMaxFilter->SetImage(maskImageFilter->GetOutput()); imgMinMaxFilter->Compute(); vnl_vector minIndex, maxIndex; 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] + (maskImage->GetOrigin()[i] - image->GetOrigin()[i]) / (double) maskImage->GetSpacing()[i]; maxIndex[i] = tmpMaxIndex[i] + (maskImage->GetOrigin()[i] - image->GetOrigin()[i]) / (double) maskImage->GetSpacing()[i]; } statisticsResult->SetMinIndex(minIndex); statisticsResult->SetMaxIndex(maxIndex); statisticsResult->SetN(imageStatisticsFilter->GetSum(*it) / (double) imageStatisticsFilter->GetMean(*it)); statisticsResult->SetMean(imageStatisticsFilter->GetMean(*it)); statisticsResult->SetMin(imageStatisticsFilter->GetMinimum(*it)); statisticsResult->SetMax(imageStatisticsFilter->GetMaximum(*it)); statisticsResult->SetVariance(imageStatisticsFilter->GetVariance(*it)); statisticsResult->SetStd(imageStatisticsFilter->GetSigma(*it)); statisticsResult->SetSkewness(imageStatisticsFilter->GetSkewness(*it)); statisticsResult->SetKurtosis(imageStatisticsFilter->GetKurtosis(*it)); statisticsResult->SetRMS(std::sqrt(std::pow(imageStatisticsFilter->GetMean(*it), 2.) + imageStatisticsFilter->GetVariance(*it))); // variance = sigma^2 statisticsResult->SetMPP(imageStatisticsFilter->GetMPP(*it)); statisticsResult->SetLabel(*it); HistogramType::Pointer histogram = imageStatisticsFilter->GetHistogram(*it); statisticsResult->SetHistogram(histogram); HistogramStatisticsCalculator histStatCalc; histStatCalc.SetHistogram(histogram); histStatCalc.CalculateStatistics(); statisticsResult->SetEntropy(histStatCalc.GetEntropy()); statisticsResult->SetMedian(histStatCalc.GetMedian()); statisticsResult->SetUniformity(histStatCalc.GetUniformity()); statisticsResult->SetUPP(histStatCalc.GetUPP()); m_StatisticsByTimeStep[timeStep].push_back(statisticsResult); ++it; } } ImageStatisticsCalculator_v2::StatisticsContainer::StatisticsContainer(): m_N(std::numeric_limits::max()), m_Mean(std::numeric_limits::max()), m_Min(std::numeric_limits::min()), m_Max(std::numeric_limits::max()), m_Std(std::numeric_limits::max()), m_Variance(std::numeric_limits::max()), m_Skewness(std::numeric_limits::max()), m_Kurtosis(std::numeric_limits::max()), m_RMS(std::numeric_limits::max()), m_MPP(std::numeric_limits::max()), m_Median(std::numeric_limits::max()), m_Uniformity(std::numeric_limits::max()), m_UPP(std::numeric_limits::max()), m_Entropy(std::numeric_limits::max()) { m_minIndex.set_size(0); m_maxIndex.set_size(0); } ImageStatisticsCalculator_v2::statisticsMapType ImageStatisticsCalculator_v2::StatisticsContainer::GetStatisticsAsMap() { ImageStatisticsCalculator_v2::statisticsMapType statisticsAsMap; statisticsAsMap["N"] = m_N; statisticsAsMap["Mean"] = m_Mean; statisticsAsMap["Min"] = m_Min; statisticsAsMap["Max"] = m_Max; statisticsAsMap["StandardDeviation"] = m_Std; statisticsAsMap["Variance"] = m_Variance; statisticsAsMap["Skewness"] = m_Skewness; statisticsAsMap["Kurtosis"] = m_Kurtosis; statisticsAsMap["RMS"] = m_RMS; statisticsAsMap["MPP"] = m_MPP; statisticsAsMap["Median"] = m_Median; statisticsAsMap["Uniformity"] = m_Uniformity; statisticsAsMap["UPP"] = m_UPP; statisticsAsMap["Entropy"] = m_Entropy; return statisticsAsMap; } void ImageStatisticsCalculator_v2::StatisticsContainer::Reset() { m_N = std::numeric_limits::max(); m_Mean = std::numeric_limits::max(); m_Min = std::numeric_limits::max(); m_Max = std::numeric_limits::max(); m_Std = std::numeric_limits::max(); m_Variance = std::numeric_limits::max(); m_Skewness = std::numeric_limits::max(); m_Kurtosis = std::numeric_limits::max(); m_RMS = std::numeric_limits::max(); m_MPP = std::numeric_limits::max(); m_Median = std::numeric_limits::max(); m_Uniformity = std::numeric_limits::max(); m_UPP = std::numeric_limits::max(); m_Entropy = std::numeric_limits::max(); } void ImageStatisticsCalculator_v2::StatisticsContainer::PrintSelf() { ImageStatisticsCalculator_v2::statisticsMapType statMap = this->GetStatisticsAsMap(); // print all map key value pairs + // const auto& val:statMap for (auto it = statMap.begin(); it != statMap.end(); ++it) { std::cout << it->first << ": " << it->second << std::endl; } // print the min and max index std::cout << "Min Index:" << std::endl; for (auto it = this->GetMinIndex().begin(); it != this->GetMinIndex().end(); it++) { std::cout << *it << " "; } std::cout << std::endl; // print the min and max index std::cout << "Max Index:" << std::endl; for (auto it = this->GetMaxIndex().begin(); it != this->GetMaxIndex().end(); it++) { std::cout << *it << " "; } std::cout << std::endl; } }