diff --git a/Modules/ImageStatistics/old/mitkExtendedLabelStatisticsImageFilter.h b/Modules/ImageStatistics/old/mitkExtendedLabelStatisticsImageFilter.h deleted file mode 100644 index ae449d14c8..0000000000 --- a/Modules/ImageStatistics/old/mitkExtendedLabelStatisticsImageFilter.h +++ /dev/null @@ -1,130 +0,0 @@ -/*=================================================================== - -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_h -#define __mitkExtendedLabelStatisticsImageFilter_h - -#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 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 itk::Statistics::Histogram HistogramType; - - itkFactorylessNewMacro( Self ); - itkCloneMacro( Self ); - itkTypeMacro(ExtendedLabelStatisticsImageFilter, LabelStatisticsImageFilter); - - /** - * \brief Internal class which stores the calculated coefficients Skewness,Kurtosis,Uniformity,UPP,MPP,Entropy. - */ - class CoefficientsClass - { - public: - - CoefficientsClass() - { - m_Kurtosis = 0.0; - m_Skewness = 0.0; - m_Entropy = -1.0; - m_Uniformity = 0.0; - m_MPP = 0.0; - m_UPP = 0.0; - }; - - ~CoefficientsClass(){}; - - /* the new member coefficients*/ - RealType m_Kurtosis; - RealType m_Skewness; - RealType m_Entropy; - RealType m_Uniformity; - RealType m_MPP; - RealType m_UPP; - }; - - - /*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; - - std::list< int> GetRelevantLabels() const; - bool GetMaskingNonEmpty() const; - - protected: - - typedef std::map< LabelPixelType, CoefficientsClass > CoefficientsMap; - typedef typename CoefficientsMap::const_iterator CoefficientsMapConstIterator; - - ExtendedLabelStatisticsImageFilter(); - - virtual ~ExtendedLabelStatisticsImageFilter(){}; - - /** - * \brief ComputeSkewnessKurtosisAndMPP(),ComputeEntropyUniformityAndUPP() will be called after superclass - * both methods are seprated because one is build up on the pixel values and one is build up a step after on a - * histogram - */ - void ComputeSkewnessKurtosisAndMPP(); - void AfterThreadedGenerateData(); - void ComputeEntropyUniformityAndUPP(); - void CalculateSettingsForLabels(); - - private: - - CoefficientsMap m_LabelStatisticsCoefficients; - std::list< int> m_RelevantLabels; - bool m_MaskNonEmpty; - - }; // end of class - -} // end namespace itk - -#ifndef ITK_MANUAL_INSTANTIATION -#include "mitkExtendedLabelStatisticsImageFilter.hxx" -#endif - -#endif diff --git a/Modules/ImageStatistics/old/mitkExtendedLabelStatisticsImageFilter.hxx b/Modules/ImageStatistics/old/mitkExtendedLabelStatisticsImageFilter.hxx deleted file mode 100644 index b7e26546b7..0000000000 --- a/Modules/ImageStatistics/old/mitkExtendedLabelStatisticsImageFilter.hxx +++ /dev/null @@ -1,323 +0,0 @@ -/*=================================================================== - -The Medical Imaging Interaction Toolkit (MITK) - -Copyright (c) German Cancer Research Center, -Division of Medical and Biological Informatics. -All rights reserved. - -This software is distributed WITHOUT ANY WARRANTY; without -even the implied warranty of MERCHANTABILITY or FITNESS FOR -A PARTICULAR PURPOSE. - -See LICENSE.txt or http://www.mitk.org for details. - -===================================================================*/ -#ifndef _mitkExtendedLabelStatisticsImageFilter_hxx -#define _mitkExtendedLabelStatisticsImageFilter_hxx - -#include "mitkExtendedLabelStatisticsImageFilter.h" - -#include "itkImageRegionConstIteratorWithIndex.h" -#include "itkImageRegionConstIterator.h" -#include -#include -#include "mitkNumericConstants.h" -#include "mitkLogMacros.h" - -namespace itk -{ - template< class TInputImage , class TLabelImage> - ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > - ::ExtendedLabelStatisticsImageFilter() - : LabelStatisticsImageFilter< TInputImage, TLabelImage >() - { - CalculateSettingsForLabels(); - } - - - template< class TInputImage, class TLabelImage > - std::list - ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > - ::GetRelevantLabels() const - { - return m_RelevantLabels; - } - - - template< class TInputImage, class TLabelImage > - bool - ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > - ::GetMaskingNonEmpty() const - { - return m_MaskNonEmpty; - } - - - template< class TInputImage, class TLabelImage > - typename ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >::RealType - ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > - ::GetUniformity(LabelPixelType label) const - { - CoefficientsMapConstIterator mapIt; - - mapIt = m_LabelStatisticsCoefficients.find(label); - if ( mapIt == m_LabelStatisticsCoefficients.end() ) - { - // label does not exist, return a default value - return NumericTraits< PixelType >::Zero; - } - else - { - return ( *mapIt ).second.m_Uniformity; - } - } - - template< class TInputImage, class TLabelImage > - typename ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >::RealType - ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > - ::GetEntropy(LabelPixelType label) const - { - CoefficientsMapConstIterator mapIt; - - mapIt = m_LabelStatisticsCoefficients.find(label); - if ( mapIt == m_LabelStatisticsCoefficients.end() ) - { - // label does not exist, return a default value - return NumericTraits< PixelType >::Zero; - } - else - { - return ( *mapIt ).second.m_Entropy; - } - } - - template< class TInputImage, class TLabelImage > - typename ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >::RealType - ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > - ::GetUPP(LabelPixelType label) const - { - CoefficientsMapConstIterator mapIt; - - mapIt = m_LabelStatisticsCoefficients.find(label); - if ( mapIt == m_LabelStatisticsCoefficients.end() ) - { - // label does not exist, return a default value - return NumericTraits< PixelType >::Zero; - } - else - { - return ( *mapIt ).second.m_UPP; - } - } - - template< class TInputImage, class TLabelImage > - typename ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >::RealType - ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > - ::GetMPP(LabelPixelType label) const - { - CoefficientsMapConstIterator mapIt; - - mapIt = m_LabelStatisticsCoefficients.find(label); - if ( mapIt == m_LabelStatisticsCoefficients.end() ) - { - // label does not exist, return a default value - return NumericTraits< PixelType >::Zero; - } - else - { - return ( *mapIt ).second.m_MPP; - } - } - - template< class TInputImage, class TLabelImage > - typename ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >::RealType - ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > - ::GetKurtosis(LabelPixelType label) const - { - CoefficientsMapConstIterator mapIt; - - mapIt = m_LabelStatisticsCoefficients.find(label); - if ( mapIt == m_LabelStatisticsCoefficients.end() ) - { - // label does not exist, return a default value - return NumericTraits< PixelType >::Zero; - } - else - { - return ( *mapIt ).second.m_Kurtosis; - } - } - - - template< class TInputImage, class TLabelImage > - typename ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >::RealType - ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > - ::GetSkewness(LabelPixelType label) const - { - CoefficientsMapConstIterator mapIt; - - mapIt = m_LabelStatisticsCoefficients.find(label); - if ( mapIt == m_LabelStatisticsCoefficients.end() ) - { - // label does not exist, return a default value - return NumericTraits< PixelType >::Zero; - } - else - { - return ( *mapIt ).second.m_Skewness; - } - } - - - template< class TInputImage, class TLabelImage > - void ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >:: - CalculateSettingsForLabels() - { - LabelPixelType i; - m_MaskNonEmpty = false; - for ( i = 1; i < 4096; ++i ) - { - if ( this->HasLabel( i ) ) - { - m_RelevantLabels.push_back( i ); - m_MaskNonEmpty = true; - m_LabelStatisticsCoefficients.insert( std::make_pair(i, CoefficientsClass()) ); - } - } - } - - - template< class TInputImage, class TLabelImage > - void - ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >:: - ComputeEntropyUniformityAndUPP() - { - RealType partialProbability( 0.0 ); - RealType uniformity( 0.0 ); - RealType entropy( 0.0 ); - RealType upp( 0.0 ); - - LabelPixelType i; - if ( m_MaskNonEmpty ) - { - typename std::list< int >::const_iterator it; - for ( it = m_RelevantLabels.cbegin(), i = 0; - it != m_RelevantLabels.cend(); - ++it, ++i ) - { - HistogramType::Pointer histogramForEntropy = this->GetHistogram(*it); - for (int i = 0; i < histogramForEntropy->Size(); i++) - { - partialProbability = histogramForEntropy->GetFrequency(i,0) / double ( histogramForEntropy->GetTotalFrequency() ) ; - - if( partialProbability != 0) - { - entropy -= partialProbability *( std::log10(partialProbability) / std::log10(2) ) ; - uniformity += std::pow(partialProbability,2); - - - if(histogramForEntropy->GetMeasurement(i,0) > 0) - { - upp += std::pow(partialProbability,2); - } - } - } - m_LabelStatisticsCoefficients[*it].m_Entropy = entropy; - m_LabelStatisticsCoefficients[*it].m_Uniformity = uniformity; - m_LabelStatisticsCoefficients[*it].m_UPP = upp; - } - } - } - - template< class TInputImage, class TLabelImage > - void - ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >:: - ComputeSkewnessKurtosisAndMPP() - { - typename TLabelImage::RegionType Subregion; - - RealType baseOfSkewnessAndCurtosis( 0.0 ); - RealType kurtosis( 0.0 ); - RealType skewness( 0.0 ); - RealType mpp( 0.0 ); - RealType currentPixel( 0.0 ); - - std::list< LabelPixelType> relevantLabels; - LabelPixelType i; - if ( m_MaskNonEmpty ) - { - typename std::list< int >::const_iterator it; - for ( it = m_RelevantLabels.cbegin(), i = 0; - it != m_RelevantLabels.cend(); - ++it ) - { - RealType sigma = this->GetSigma( *it ); - RealType mean = this->GetMean( *it ); - Subregion = Superclass::GetRegion(*it); - - int count( this->GetCount(*it) ); - if ( count == 0 || sigma < mitk::eps) - { - throw std::logic_error( "Empty segmentation" ); - } - - if ( fabs( sigma ) < mitk::sqrteps ) - { - throw std::logic_error( "Sigma == 0" ); - } - - - ImageRegionConstIteratorWithIndex< TInputImage > it1 (this->GetInput(), - Subregion); - ImageRegionConstIterator< TLabelImage > labelIt (this->GetLabelInput(), - Subregion); - - for (it1.GoToBegin(); !it1.IsAtEnd(); ++it1, ++labelIt) - { - if (labelIt.Get() == *it) - { - currentPixel = it1.Get(); - baseOfSkewnessAndCurtosis = (currentPixel -mean) / sigma; - kurtosis += std::pow( baseOfSkewnessAndCurtosis, 4.0 ); - skewness += std::pow( baseOfSkewnessAndCurtosis, 3.0 ); - - if(currentPixel > 0) - { - mpp+= currentPixel; - } - - } - } - m_LabelStatisticsCoefficients[*it].m_Skewness = RealType(skewness/count); - m_LabelStatisticsCoefficients[*it].m_Kurtosis = RealType(kurtosis/count); - m_LabelStatisticsCoefficients[*it].m_MPP = RealType(mpp/count); - } - } - } - - - - template< class TInputImage, class TLabelImage > - void ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >:: - AfterThreadedGenerateData() - { - Superclass::AfterThreadedGenerateData(); - - CalculateSettingsForLabels(); - - ComputeSkewnessKurtosisAndMPP(); - - if(this->GetUseHistograms()) - { - ComputeEntropyUniformityAndUPP(); - } - else - { - MITK_WARN << "Cannot compute coefficients UPP,Entropy,Uniformity because of missing histogram"; - } - } - -} // end namespace itk - -#endif diff --git a/Modules/ImageStatistics/old/mitkExtendedStatisticsImageFilter.h b/Modules/ImageStatistics/old/mitkExtendedStatisticsImageFilter.h deleted file mode 100644 index ea0927d30c..0000000000 --- a/Modules/ImageStatistics/old/mitkExtendedStatisticsImageFilter.h +++ /dev/null @@ -1,205 +0,0 @@ -/*=================================================================== - -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_h -#define __mitkExtendedStatisticsImageFilter_h - -#include "itkStatisticsImageFilter.h" -#include "itkScalarImageToHistogramGenerator.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 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; - - itkFactorylessNewMacro( Self ); - itkCloneMacro( Self ); - itkTypeMacro( ExtendedStatisticsImageFilter, StatisticsImageFilter ); - - - typedef itk::Statistics::ScalarImageToHistogramGenerator< TInputImage > - HistogramGeneratorType; - - /** - * \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 HistogramGeneratorType::HistogramType* - GetHistogram() - { - return m_HistogramGenerator->GetOutput(); - } - - /** - * \brief Set the Binsize for the Histogram. - */ - void SetBinSize(int size); - - protected: - - ExtendedStatisticsImageFilter(); - - virtual ~ExtendedStatisticsImageFilter(){}; - - /** - * brief Calls AfterThreadedGenerateData() of the superclass and the main methods - */ - void AfterThreadedGenerateData(); - void CalculateHistogram(); - - /** - * \brief Compute Entropy,uniformity,MPP,UPP, Median. - * - * The Entropy,uniformity,MPP, Median and UPP will be calculated with the Sigma, Histogram and Mean Value of the - * itkStatisticsImageFilter which comes out of the threadedGenerateData(). - */ - void ComputeSkewnessKurtosisAndMPP(); - void ComputeEntropyUniformityMedianAndUPP(); - - /** - * \brief Histogram. - * - * new members for setting and calculating the hisotgram for those coefficients which depends on this - */ - typename HistogramGeneratorType::Pointer m_HistogramGenerator; - int m_BinSize; - bool m_HistogramCalculated; - - - 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 ); - - - }; // end of class - -} // end namespace itk - -#ifndef ITK_MANUAL_INSTANTIATION -#include "mitkExtendedStatisticsImageFilter.hxx" -#endif - -#endif diff --git a/Modules/ImageStatistics/old/mitkExtendedStatisticsImageFilter.hxx b/Modules/ImageStatistics/old/mitkExtendedStatisticsImageFilter.hxx deleted file mode 100644 index 3f3d810edb..0000000000 --- a/Modules/ImageStatistics/old/mitkExtendedStatisticsImageFilter.hxx +++ /dev/null @@ -1,343 +0,0 @@ -/*=================================================================== - -The Medical Imaging Interaction Toolkit (MITK) - -Copyright (c) German Cancer Research Center, -Division of Medical and Biological Informatics. -All rights reserved. - -This software is distributed WITHOUT ANY WARRANTY; without -even the implied warranty of MERCHANTABILITY or FITNESS FOR -A PARTICULAR PURPOSE. - -See LICENSE.txt or http://www.mitk.org for details. - -===================================================================*/ -#ifndef __mitkExtendedStatisticsImageFilter_hxx -#define __mitkExtendedStatisticsImageFilter_hxx - -#include "mitkExtendedStatisticsImageFilter.h" - - -namespace itk -{ - template< class TInputImage > - ExtendedStatisticsImageFilter< TInputImage >::ExtendedStatisticsImageFilter() - : StatisticsImageFilter< TInputImage >() - { - /* - * add the Skewness,Kurtosis,Entropy,Uniformity,MPP, UPP, Median to the other statistical calculated Values - * of the mitkStatisticsImageFilter as the 7th to the 13th Output - */ - for ( int i = 7; i < 14; ++i ) - { - typename RealObjectType::Pointer output = - static_cast< RealObjectType * >( this->MakeOutput(i).GetPointer() ); - this->ProcessObject::SetNthOutput( i, output.GetPointer() ); - } - - this->GetSkewnessOutput()->Set( 0.0 ); - this->GetKurtosisOutput()->Set( 0.0 ); - this->GetEntropyOutput()->Set( -1.0 ); - this->GetUniformityOutput()->Set( 0.0 ); - this->GetUPPOutput()->Set( 0.0 ); - this->GetMPPOutput()->Set( 0.0 ); - this->GetMedianOutput()->Set( 0.0 ); - - this->m_HistogramCalculated = false; - this->m_HistogramGenerator = HistogramGeneratorType::New(); - } - - template< class TInputImage > - DataObject::Pointer - ExtendedStatisticsImageFilter< TInputImage >::MakeOutput( ProcessObject::DataObjectPointerArraySizeType output) - { - switch ( output ) - { - case 7: - case 8: - case 9: - case 10: - case 11: - case 12: - { - return RealObjectType::New().GetPointer(); - break; - } - case 13: - { - return RealObjectType::New().GetPointer(); - break; - } - default: - { - // might as well make an image - return Superclass::MakeOutput( output ); - break; - } - } - } - - - template< class TInputImage > - void ExtendedStatisticsImageFilter< TInputImage > - ::SetBinSize( int size ) - { - m_BinSize = size; - } - - template< class TInputImage > - typename ExtendedStatisticsImageFilter< TInputImage >::RealObjectType * - ExtendedStatisticsImageFilter< TInputImage > - ::GetSkewnessOutput() - { - return static_cast< RealObjectType * >( this->ProcessObject::GetOutput(7) ); - } - - template< class TInputImage > - const typename ExtendedStatisticsImageFilter< TInputImage >::RealObjectType * - ExtendedStatisticsImageFilter< TInputImage > - ::GetSkewnessOutput() const - { - return static_cast< const RealObjectType * >( this->ProcessObject::GetOutput(7) ); - } - - - template< class TInputImage > - typename ExtendedStatisticsImageFilter< TInputImage >::RealObjectType * - ExtendedStatisticsImageFilter< TInputImage > - ::GetKurtosisOutput() - { - return static_cast< RealObjectType * >( this->ProcessObject::GetOutput(8) ); - } - - template< class TInputImage > - const typename ExtendedStatisticsImageFilter< TInputImage >::RealObjectType * - ExtendedStatisticsImageFilter< TInputImage > - ::GetKurtosisOutput() const - { - return static_cast< const RealObjectType * >( this->ProcessObject::GetOutput(8) ); - } - - template< class TInputImage > - typename ExtendedStatisticsImageFilter< TInputImage >::RealObjectType * - ExtendedStatisticsImageFilter< TInputImage > - ::GetUniformityOutput() - { - return static_cast< RealObjectType * >( this->ProcessObject::GetOutput(9) ); - } - - template< class TInputImage > - const typename ExtendedStatisticsImageFilter< TInputImage >::RealObjectType * - ExtendedStatisticsImageFilter< TInputImage > - ::GetUniformityOutput() const - { - return static_cast< const RealObjectType * >( this->ProcessObject::GetOutput(9) ); - } - - - template< class TInputImage > - typename ExtendedStatisticsImageFilter< TInputImage >::RealObjectType * - ExtendedStatisticsImageFilter< TInputImage > - ::GetEntropyOutput() - { - return static_cast< RealObjectType * >( this->ProcessObject::GetOutput(10) ); - } - - template< class TInputImage > - const typename ExtendedStatisticsImageFilter< TInputImage >::RealObjectType * - ExtendedStatisticsImageFilter< TInputImage > - ::GetEntropyOutput() const - { - return static_cast< const RealObjectType * >( this->ProcessObject::GetOutput(10) ); - } - - template< class TInputImage > - typename ExtendedStatisticsImageFilter< TInputImage >::RealObjectType * - ExtendedStatisticsImageFilter< TInputImage > - ::GetUPPOutput() - { - return static_cast< RealObjectType * >( this->ProcessObject::GetOutput(11) ); - } - - template< class TInputImage > - const typename ExtendedStatisticsImageFilter< TInputImage >::RealObjectType * - ExtendedStatisticsImageFilter< TInputImage > - ::GetUPPOutput() const - { - return static_cast< const RealObjectType * >( this->ProcessObject::GetOutput(11) ); - } - - template< class TInputImage > - typename ExtendedStatisticsImageFilter< TInputImage >::RealObjectType * - ExtendedStatisticsImageFilter< TInputImage > - ::GetMPPOutput() - { - return static_cast< RealObjectType * >( this->ProcessObject::GetOutput(12) ); - } - - template< class TInputImage > - const typename ExtendedStatisticsImageFilter< TInputImage >::RealObjectType * - ExtendedStatisticsImageFilter< TInputImage > - ::GetMPPOutput() const - { - return static_cast< const RealObjectType * >( this->ProcessObject::GetOutput(12) ); - } - - template< class TInputImage > - typename ExtendedStatisticsImageFilter< TInputImage >::RealObjectType * - ExtendedStatisticsImageFilter< TInputImage > - ::GetMedianOutput() - { - return static_cast< RealObjectType * >( this->ProcessObject::GetOutput(13) ); - } - - template< class TInputImage > - const typename ExtendedStatisticsImageFilter< TInputImage >::RealObjectType * - ExtendedStatisticsImageFilter< TInputImage > - ::GetMedianOutput() const - { - return static_cast< const RealObjectType * >( this->ProcessObject::GetOutput(13) ); - } - - - - template< class TInputImage > - void - ExtendedStatisticsImageFilter< TInputImage > - ::AfterThreadedGenerateData() - { - Superclass::AfterThreadedGenerateData(); - - ComputeSkewnessKurtosisAndMPP(); - - CalculateHistogram(); - - if(m_HistogramCalculated == true) - { - ComputeEntropyUniformityMedianAndUPP(); - } - else - { - MITK_WARN << "Cannot compute coefficients UPP,Entropy,Uniformity because of missing histogram"; - } - } - - template< class TInputImage > - void - ExtendedStatisticsImageFilter< TInputImage > - ::ComputeSkewnessKurtosisAndMPP() - { - RealType mean = this->GetMean(); - RealType sigma = this->GetSigma(); - RealType baseOfSkewnessAndKurtosis; - RealType kurtosis(0.0); - RealType skewness(0.0); - RealType mpp(0.0); - RealType currentPixel(0.0); - - if ( sigma < mitk::eps ) - { - throw std::logic_error( "Empty segmentation" ); - } - - - ImageRegionConstIterator< TInputImage > it (this->GetInput(), this->GetInput()->GetLargestPossibleRegion() ); - - int counter = 0; - for (it.GoToBegin(); !it.IsAtEnd(); ++it) - { - currentPixel = it.Get(); - - baseOfSkewnessAndKurtosis = (currentPixel - mean) / sigma ; - kurtosis += std::pow( baseOfSkewnessAndKurtosis, 4.0 ); - skewness += std::pow( baseOfSkewnessAndKurtosis, 3.0 ); - - if(currentPixel > 0) - { - mpp+= currentPixel; - } - counter++; - } - - if ( counter == 0 ) - { - throw std::logic_error( "Empty segmentation" ); - } - - kurtosis = kurtosis / counter; - skewness = skewness / counter; - mpp = mpp/counter; - - this->GetKurtosisOutput()->Set( kurtosis ); - this->GetSkewnessOutput()->Set( skewness ); - this->GetMPPOutput()->Set( mpp ); - } - - - template< class TInputImage > - void - ExtendedStatisticsImageFilter< TInputImage > - ::CalculateHistogram() - { - m_HistogramGenerator->SetInput( this->GetInput() ); - m_HistogramGenerator->SetMarginalScale( 100 ); - m_HistogramGenerator->SetNumberOfBins( m_BinSize ); - m_HistogramGenerator->SetHistogramMin( this->GetMinimum() ); - m_HistogramGenerator->SetHistogramMax( this->GetMaximum() ); - m_HistogramGenerator->Compute(); - - m_HistogramCalculated = true; - } - - - template< class TInputImage > - void - ExtendedStatisticsImageFilter< TInputImage > - ::ComputeEntropyUniformityMedianAndUPP() - { - RealType partialProbability( 0.0 ); - RealType uniformity( 0.0 ); - RealType entropy( 0.0 ); - RealType upp( 0.0 ); - RealType median( 0.0 ); - - const typename HistogramGeneratorType::HistogramType* histogramForEntropy = GetHistogram(); - - double cumulativeProbability = 0.0; - bool medianFound = false; - - for (int i = 0; i < histogramForEntropy->Size(); i++) - { - partialProbability = histogramForEntropy->GetFrequency(i,0) / double ( histogramForEntropy->GetTotalFrequency() ) ; - cumulativeProbability += double ( partialProbability ); - - if( partialProbability != 0) - { - entropy -= partialProbability *( std::log10(partialProbability) / std::log10(2) ) ; - uniformity += std::pow(partialProbability,2); - - if(histogramForEntropy->GetMeasurement(i,0) > 0) - { - upp += std::pow(partialProbability,2); - } - } - - if( cumulativeProbability >= 0.5 && !medianFound ) - { - RealType binMin = histogramForEntropy->GetBinMin( 0, i ); - RealType binMax = histogramForEntropy->GetBinMax( 0, i ); - median = ( binMin + binMax ) / 2.0; - medianFound = true; - } - - } - this->GetEntropyOutput()->Set( entropy ); - this->GetUniformityOutput()->Set( uniformity ); - this->GetUPPOutput()->Set( upp ); - this->GetMedianOutput()->Set( median ); - - } -} - -#endif diff --git a/Modules/ImageStatistics/old/mitkImageStatisticsCalculator.cpp b/Modules/ImageStatistics/old/mitkImageStatisticsCalculator.cpp deleted file mode 100644 index 612aeb71f6..0000000000 --- a/Modules/ImageStatistics/old/mitkImageStatisticsCalculator.cpp +++ /dev/null @@ -1,2246 +0,0 @@ -/*=================================================================== - -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 "mitkImageStatisticsCalculator.h" -#include "mitkImageAccessByItk.h" -#include "mitkImageCast.h" -#include "mitkExtractImageFilter.h" -#include "mitkImageTimeSelector.h" -#include "mitkITKImageImport.h" - -#include -#include - -#include - -#include -#include -#include -#include -#include -#include -#include -#include - -#include - -#include -#include - -#include -#include -#include -#include - -#include -#include -#include -#include -#include - -#include -#include -#include - -#include -#include -#include - -#include - -#include "itkImage.h" - - -//#define DEBUG_HOTSPOTSEARCH - -#define _USE_MATH_DEFINES -#include - -#include "vtkLassoStencilSource.h" - - -namespace mitk -{ - ImageStatisticsCalculator::ImageStatisticsCalculator() - : m_MaskingMode( MASKING_MODE_NONE ), - m_MaskingModeChanged( false ), - m_IgnorePixelValue(0.0), - m_DoIgnorePixelValue(false), - m_IgnorePixelValueChanged(false), - m_PlanarFigureAxis (0), - m_PlanarFigureSlice (0), - m_PlanarFigureCoordinate0 (0), - m_PlanarFigureCoordinate1 (0), - m_HistogramBinSize(1.0), - m_UseDefaultBinSize(true), - m_UseBinSizeBasedOnVOIRegion(false), - m_HotspotRadiusInMM(6.2035049089940), // radius of a 1cm3 sphere in mm - m_CalculateHotspot(false), - m_HotspotRadiusInMMChanged(false), - m_HotspotMustBeCompletelyInsideImage(true) - { - m_EmptyHistogram = HistogramType::New(); - m_EmptyHistogram->SetMeasurementVectorSize(1); - HistogramType::SizeType histogramSize(1); - histogramSize.Fill( 256 ); - m_EmptyHistogram->Initialize( histogramSize ); - - m_EmptyStatistics.Reset(); - } - - ImageStatisticsCalculator::~ImageStatisticsCalculator() - { - } - - - void ImageStatisticsCalculator::SetUseDefaultBinSize(bool useDefault) - { - m_UseDefaultBinSize = useDefault; - } - - - - - ImageStatisticsCalculator::Statistics::Statistics(bool withHotspotStatistics) - :m_HotspotStatistics(withHotspotStatistics ? new Statistics(false) : nullptr) - { - Reset(); - } - - ImageStatisticsCalculator::Statistics::Statistics(const Statistics& other) - :m_HotspotStatistics( nullptr) - { - this->SetLabel( other.GetLabel() ); - this->SetN( other.GetN() ); - this->SetMin( other.GetMin() ); - this->SetMax( other.GetMax() ); - this->SetMedian( other.GetMedian() ); - this->SetMean( other.GetMean() ); - this->SetVariance( other.GetVariance() ); - this->SetKurtosis( other.GetKurtosis() ); - this->SetSkewness( other.GetSkewness() ); - this->SetUniformity( other.GetUniformity() ); - this->SetEntropy( other.GetEntropy() ); - this->SetUPP( other.GetUPP() ); - this->SetMPP( other.GetMPP() ); - this->SetSigma( other.GetSigma() ); - this->SetRMS( other.GetRMS() ); - this->SetMaxIndex( other.GetMaxIndex() ); - this->SetMinIndex( other.GetMinIndex() ); - this->SetHotspotIndex( other.GetHotspotIndex() ); - - if (other.m_HotspotStatistics) - { - this->m_HotspotStatistics = new Statistics(false); - *this->m_HotspotStatistics = *other.m_HotspotStatistics; - } - } - - bool ImageStatisticsCalculator::Statistics::HasHotspotStatistics() const - { - return m_HotspotStatistics != nullptr; - } - - void ImageStatisticsCalculator::Statistics::SetHasHotspotStatistics(bool hasHotspotStatistics) - { - m_HasHotspotStatistics = hasHotspotStatistics; - } - - - ImageStatisticsCalculator::Statistics::~Statistics() - { - delete m_HotspotStatistics; - } - - double ImageStatisticsCalculator::Statistics::GetVariance() const - { - return this->Variance; - } - - void ImageStatisticsCalculator::Statistics::SetVariance( const double value ) - { - if( this->Variance != value ) - { - if( value < 0.0 ) - { - this->Variance = 0.0; // if given value is negative set variance to 0.0 - } - else - { - this->Variance = value; - } - } - } - - double ImageStatisticsCalculator::Statistics::GetSigma() const - { - return this->Sigma; - } - - void ImageStatisticsCalculator::Statistics::SetSigma( const double value ) - { - if( this->Sigma != value ) - { - // for some compiler the value != value works to check for NaN but not for all - // but we can always be sure that the standard deviation is a positive value - if( value != value || value < 0.0 ) - { - // if standard deviation is NaN we just assume 0.0 - this->Sigma = 0.0; - } - else - { - this->Sigma = value; - } - } - } - - void ImageStatisticsCalculator::Statistics::Reset(unsigned int dimension) - { - SetLabel(0); - SetN( 0 ); - SetMin( 0.0 ); - SetMax( 0.0 ); - SetMedian( 0.0 ); - SetVariance( 0.0 ); - SetMean( 0.0 ); - SetSigma( 0.0 ); - SetRMS( 0.0 ); - - SetSkewness( 0.0 ); - SetKurtosis( 0.0 ); - SetUniformity( 0.0 ); - SetEntropy( 0.0 ); - SetMPP( 0.0 ); - SetUPP( 0.0 ); - - vnl_vector zero; - zero.set_size(dimension); - for(unsigned int i = 0; i < dimension; ++i) - { - zero[i] = 0; - } - - SetMaxIndex(zero); - SetMinIndex(zero); - SetHotspotIndex(zero); - - if (m_HotspotStatistics != nullptr) - { - m_HotspotStatistics->Reset(dimension); - } - } - - const ImageStatisticsCalculator::Statistics& - ImageStatisticsCalculator::Statistics::GetHotspotStatistics() const - { - if (m_HotspotStatistics) - { - return *m_HotspotStatistics; - } - else - { - throw std::logic_error("Object has no hostspot statistics, see HasHotspotStatistics()"); - } - } - - ImageStatisticsCalculator::Statistics& - ImageStatisticsCalculator::Statistics::GetHotspotStatistics() - { - if (m_HotspotStatistics) - { - return *m_HotspotStatistics; - } - else - { - throw std::logic_error("Object has no hostspot statistics, see HasHotspotStatistics()"); - } - } - - ImageStatisticsCalculator::Statistics& - ImageStatisticsCalculator::Statistics::operator=(ImageStatisticsCalculator::Statistics const& other) - { - if (this == &other) - return *this; - - this->SetLabel( other.GetLabel() ); - this->SetN( other.GetN() ); - this->SetMin( other.GetMin() ); - this->SetMax( other.GetMax() ); - this->SetMean( other.GetMean() ); - this->SetMedian( other.GetMedian() ); - this->SetVariance( other.GetVariance() ); - this->SetSigma( other.GetSigma() ); - this->SetRMS( other.GetRMS() ); - this->SetMinIndex( other.GetMinIndex() ); - this->SetMaxIndex( other.GetMaxIndex() ); - this->SetHotspotIndex( other.GetHotspotIndex() ); - this->SetSkewness( other.GetSkewness() ); - this->SetKurtosis( other.GetKurtosis() ); - this->SetUniformity( other.GetUniformity() ); - this->SetEntropy( other.GetEntropy() ); - this->SetUPP( other.GetUPP() ); - this->SetMPP( other.GetMPP() ); - - - delete this->m_HotspotStatistics; - this->m_HotspotStatistics = nullptr; - - if (other.m_HotspotStatistics) - { - this->m_HotspotStatistics = new Statistics(false); - *this->m_HotspotStatistics = *other.m_HotspotStatistics; - } - - return *this; - - } - - void ImageStatisticsCalculator::SetImage( const mitk::Image *image ) - { - if ( m_Image != image ) - { - m_Image = image; - this->Modified(); - unsigned int numberOfTimeSteps = image->GetTimeSteps(); - - // Initialize vectors to time-size of this image - m_ImageHistogramVector.resize( numberOfTimeSteps ); - m_MaskedImageHistogramVector.resize( numberOfTimeSteps ); - m_PlanarFigureHistogramVector.resize( numberOfTimeSteps ); - - m_ImageStatisticsVector.resize( numberOfTimeSteps ); - m_MaskedImageStatisticsVector.resize( numberOfTimeSteps ); - m_PlanarFigureStatisticsVector.resize( numberOfTimeSteps ); - - m_ImageStatisticsTimeStampVector.resize( numberOfTimeSteps ); - m_MaskedImageStatisticsTimeStampVector.resize( numberOfTimeSteps ); - m_PlanarFigureStatisticsTimeStampVector.resize( numberOfTimeSteps ); - - m_ImageStatisticsCalculationTriggerVector.resize( numberOfTimeSteps ); - m_MaskedImageStatisticsCalculationTriggerVector.resize( numberOfTimeSteps ); - m_PlanarFigureStatisticsCalculationTriggerVector.resize( numberOfTimeSteps ); - - for ( unsigned int t = 0; t < image->GetTimeSteps(); ++t ) - { - m_ImageStatisticsTimeStampVector[t].Modified(); - m_ImageStatisticsCalculationTriggerVector[t] = true; - } - } - } - - - void ImageStatisticsCalculator::SetImageMask( const mitk::Image *imageMask ) - { - if ( m_Image.IsNull() ) - { - itkExceptionMacro( << "Image needs to be set first!" ); - } - - if ( m_ImageMask != imageMask ) - { - m_ImageMask = imageMask; - this->Modified(); - - for ( unsigned int t = 0; t < m_Image->GetTimeSteps(); ++t ) - { - m_MaskedImageStatisticsTimeStampVector[t].Modified(); - m_MaskedImageStatisticsCalculationTriggerVector[t] = true; - } - } - } - - - void ImageStatisticsCalculator::SetPlanarFigure( mitk::PlanarFigure *planarFigure ) - { - if ( m_Image.IsNull() ) - { - itkExceptionMacro( << "Image needs to be set first!" ); - } - - if ( m_PlanarFigure != planarFigure ) - { - m_PlanarFigure = planarFigure; - this->Modified(); - - for ( unsigned int t = 0; t < m_Image->GetTimeSteps(); ++t ) - { - m_PlanarFigureStatisticsTimeStampVector[t].Modified(); - m_PlanarFigureStatisticsCalculationTriggerVector[t] = true; - } - } - } - - - void ImageStatisticsCalculator::SetMaskingMode( unsigned int mode ) - { - if ( m_MaskingMode != mode ) - { - m_MaskingMode = mode; - m_MaskingModeChanged = true; - this->Modified(); - } - } - - - void ImageStatisticsCalculator::SetMaskingModeToNone() - { - if ( m_MaskingMode != MASKING_MODE_NONE ) - { - m_MaskingMode = MASKING_MODE_NONE; - m_MaskingModeChanged = true; - this->Modified(); - } - } - - - void ImageStatisticsCalculator::SetMaskingModeToImage() - { - if ( m_MaskingMode != MASKING_MODE_IMAGE ) - { - m_MaskingMode = MASKING_MODE_IMAGE; - m_MaskingModeChanged = true; - this->Modified(); - } - } - - - void ImageStatisticsCalculator::SetMaskingModeToPlanarFigure() - { - if ( m_MaskingMode != MASKING_MODE_PLANARFIGURE ) - { - m_MaskingMode = MASKING_MODE_PLANARFIGURE; - m_MaskingModeChanged = true; - this->Modified(); - } - } - - void ImageStatisticsCalculator::SetIgnorePixelValue(double value) - { - if ( m_IgnorePixelValue != value ) - { - m_IgnorePixelValue = value; - if(m_DoIgnorePixelValue) - { - m_IgnorePixelValueChanged = true; - } - this->Modified(); - } - } - - double ImageStatisticsCalculator::GetIgnorePixelValue() - { - return m_IgnorePixelValue; - } - - void ImageStatisticsCalculator::SetDoIgnorePixelValue(bool value) - { - if ( m_DoIgnorePixelValue != value ) - { - m_DoIgnorePixelValue = value; - m_IgnorePixelValueChanged = true; - this->Modified(); - } - } - - bool ImageStatisticsCalculator::GetDoIgnorePixelValue() - { - return m_DoIgnorePixelValue; - } - - void ImageStatisticsCalculator::SetHistogramBinSize(double size) - { - this->m_HistogramBinSize = size; - } - double ImageStatisticsCalculator::GetHistogramBinSize() - { - return this->m_HistogramBinSize; - } - - void ImageStatisticsCalculator::SetHotspotRadiusInMM(double value) - { - if ( m_HotspotRadiusInMM != value ) - { - m_HotspotRadiusInMM = value; - if(m_CalculateHotspot) - { - m_HotspotRadiusInMMChanged = true; - //MITK_INFO <<"Hotspot radius changed, new convolution required"; - } - this->Modified(); - } - } - - double ImageStatisticsCalculator::GetHotspotRadiusInMM() - { - return m_HotspotRadiusInMM; - } - - void ImageStatisticsCalculator::SetCalculateHotspot(bool on) - { - if ( m_CalculateHotspot != on ) - { - m_CalculateHotspot = on; - m_HotspotRadiusInMMChanged = true; - //MITK_INFO <<"Hotspot calculation changed, new convolution required"; - this->Modified(); - } - } - - bool ImageStatisticsCalculator::IsHotspotCalculated() - { - return m_CalculateHotspot; - } - - void ImageStatisticsCalculator::SetHotspotMustBeCompletlyInsideImage(bool hotspotMustBeCompletelyInsideImage, bool warn) - { - m_HotspotMustBeCompletelyInsideImage = hotspotMustBeCompletelyInsideImage; - if (!m_HotspotMustBeCompletelyInsideImage && warn) - { - MITK_WARN << "Hotspot calculation will extrapolate pixels at image borders. Be aware of the consequences for the hotspot location."; - } - } - - bool ImageStatisticsCalculator::GetHotspotMustBeCompletlyInsideImage() const - { - return m_HotspotMustBeCompletelyInsideImage; - } - - - /* Implementation of the min max values for setting the range of the histogram */ - template < typename TPixel, unsigned int VImageDimension > - void ImageStatisticsCalculator::GetMinAndMaxValue( double &min, - double &max, int &counter, double &sigma, - const itk::Image< TPixel, VImageDimension > *InputImage, - itk::Image< unsigned short, VImageDimension > *MaskedImage ) - { - typedef itk::Image< unsigned short, VImageDimension > MaskImageType; - typedef itk::Image< TPixel, VImageDimension > ImageType; - - typedef itk::ImageRegionConstIteratorWithIndex Imageie; - typedef itk::ImageRegionConstIteratorWithIndex Imageie2; - - Imageie2 labelIterator2( MaskedImage, MaskedImage->GetRequestedRegion() ); - Imageie labelIterator3( InputImage, InputImage->GetRequestedRegion() ); - - max = 0; - min = 0; - counter = 0; - sigma = 0; - double SumOfSquares = 0; - double sumSquared = 0; - - double actualPielValue = 0; - int counterOfPixelsInROI = 0; - - - for( labelIterator2.GoToBegin(); !labelIterator2.IsAtEnd(); ++labelIterator2, ++labelIterator3) - { - if( labelIterator2.Value()== 1.0) - { - counter++; - - counterOfPixelsInROI++; - actualPielValue = labelIterator3.Value(); - - sumSquared = sumSquared + actualPielValue; - SumOfSquares = SumOfSquares + std::pow(actualPielValue,2); - - if(counterOfPixelsInROI == 1) - { - max = actualPielValue; - min = actualPielValue; - } - - if(actualPielValue >= max) - { - max = actualPielValue; - } - else if(actualPielValue <= min) - { - min = actualPielValue; - } - - } - } - - if (counter > 1) - { - sigma = ( SumOfSquares - std::pow( sumSquared, 2) / counter ) / ( counter-1 ); - } - else - { - sigma = 0; - } - - } - - - bool ImageStatisticsCalculator::ComputeStatistics( unsigned int timeStep ) - { - - if (m_Image.IsNull() ) - { - mitkThrow() << "Image not set!"; - } - - if (!m_Image->IsInitialized()) - { - mitkThrow() << "Image not initialized!"; - } - - if ( m_Image->GetReferenceCount() == 1 ) - { - // Image no longer valid; we are the only ones to still hold a reference on it - return false; - } - - if ( timeStep >= m_Image->GetTimeSteps() ) - { - throw std::runtime_error( "Error: invalid time step!" ); - } - - // If a mask was set but we are the only ones to still hold a reference on - // it, delete it. - if ( m_ImageMask.IsNotNull() && (m_ImageMask->GetReferenceCount() == 1) ) - { - m_ImageMask = nullptr; - } - - // Check if statistics is already up-to-date - unsigned long imageMTime = m_ImageStatisticsTimeStampVector[timeStep].GetMTime(); - unsigned long maskedImageMTime = m_MaskedImageStatisticsTimeStampVector[timeStep].GetMTime(); - unsigned long planarFigureMTime = m_PlanarFigureStatisticsTimeStampVector[timeStep].GetMTime(); - - bool imageStatisticsCalculationTrigger = m_ImageStatisticsCalculationTriggerVector[timeStep]; - bool maskedImageStatisticsCalculationTrigger = m_MaskedImageStatisticsCalculationTriggerVector[timeStep]; - bool planarFigureStatisticsCalculationTrigger = m_PlanarFigureStatisticsCalculationTriggerVector[timeStep]; - - if ( !m_IgnorePixelValueChanged - && !m_HotspotRadiusInMMChanged - && ((m_MaskingMode != MASKING_MODE_NONE) || (imageMTime > m_Image->GetMTime() && !imageStatisticsCalculationTrigger)) - && ((m_MaskingMode != MASKING_MODE_IMAGE) || (maskedImageMTime > m_ImageMask->GetMTime() && !maskedImageStatisticsCalculationTrigger)) - && ((m_MaskingMode != MASKING_MODE_PLANARFIGURE) || (planarFigureMTime > m_PlanarFigure->GetMTime() && !planarFigureStatisticsCalculationTrigger)) ) - { - // Statistics is up to date! - if ( m_MaskingModeChanged ) - { - m_MaskingModeChanged = false; - } - else - { - return false; - } - } - - // Reset state changed flag - m_MaskingModeChanged = false; - m_IgnorePixelValueChanged = false; - - // Depending on masking mode, extract and/or generate the required image - // and mask data from the user input - this->ExtractImageAndMask( timeStep ); - - - StatisticsContainer *statisticsContainer; - HistogramContainer *histogramContainer; - switch ( m_MaskingMode ) - { - case MASKING_MODE_NONE: - default: - if(!m_DoIgnorePixelValue) - { - statisticsContainer = &m_ImageStatisticsVector[timeStep]; - histogramContainer = &m_ImageHistogramVector[timeStep]; - - m_ImageStatisticsTimeStampVector[timeStep].Modified(); - m_ImageStatisticsCalculationTriggerVector[timeStep] = false; - } - else - { - statisticsContainer = &m_MaskedImageStatisticsVector[timeStep]; - histogramContainer = &m_MaskedImageHistogramVector[timeStep]; - - m_MaskedImageStatisticsTimeStampVector[timeStep].Modified(); - m_MaskedImageStatisticsCalculationTriggerVector[timeStep] = false; - } - break; - - case MASKING_MODE_IMAGE: - statisticsContainer = &m_MaskedImageStatisticsVector[timeStep]; - histogramContainer = &m_MaskedImageHistogramVector[timeStep]; - - m_MaskedImageStatisticsTimeStampVector[timeStep].Modified(); - m_MaskedImageStatisticsCalculationTriggerVector[timeStep] = false; - break; - - case MASKING_MODE_PLANARFIGURE: - statisticsContainer = &m_PlanarFigureStatisticsVector[timeStep]; - histogramContainer = &m_PlanarFigureHistogramVector[timeStep]; - - m_PlanarFigureStatisticsTimeStampVector[timeStep].Modified(); - m_PlanarFigureStatisticsCalculationTriggerVector[timeStep] = false; - break; - } - - // Calculate statistics and histogram(s) - if ( m_InternalImage->GetDimension() == 3 ) - { - if ( m_MaskingMode == MASKING_MODE_NONE && !m_DoIgnorePixelValue ) - { - AccessFixedDimensionByItk_2( - m_InternalImage, - InternalCalculateStatisticsUnmasked, - 3, - statisticsContainer, - histogramContainer ); - } - else - { - AccessFixedDimensionByItk_3( - m_InternalImage, - InternalCalculateStatisticsMasked, - 3, - m_InternalImageMask3D.GetPointer(), - statisticsContainer, - histogramContainer ); - } - } - else if ( m_InternalImage->GetDimension() == 2 ) - { - if ( m_MaskingMode == MASKING_MODE_NONE && !m_DoIgnorePixelValue ) - { - AccessFixedDimensionByItk_2( - m_InternalImage, - InternalCalculateStatisticsUnmasked, - 2, - statisticsContainer, - histogramContainer ); - } - else - { - AccessFixedDimensionByItk_3( - m_InternalImage, - InternalCalculateStatisticsMasked, - 2, - m_InternalImageMask2D.GetPointer(), - statisticsContainer, - histogramContainer ); - } - } - else - { - MITK_ERROR << "ImageStatistics: Image dimension not supported!"; - } - - - // Release unused image smart pointers to free memory - m_InternalImage = mitk::Image::ConstPointer(); - m_InternalImageMask3D = MaskImage3DType::Pointer(); - m_InternalImageMask2D = MaskImage2DType::Pointer(); - - - - return true; - } - - - ImageStatisticsCalculator::BinFrequencyType - ImageStatisticsCalculator::GetBinsAndFreuqencyForHistograms( unsigned int timeStep , unsigned int label ) const - { - const HistogramType *binsAndFrequencyToCalculate = this->GetHistogram(0); - - // ToDo: map should be created on stack not on heap - std::map returnedHistogramMap; - - unsigned int size = binsAndFrequencyToCalculate->Size(); - for( unsigned int bin=0; bin < size; ++bin ) - { - double frequency = binsAndFrequencyToCalculate->GetFrequency( bin, 0 ); - //if( frequency > mitk::eps ) - { - returnedHistogramMap.insert( std::pair(binsAndFrequencyToCalculate->GetMeasurement( bin, 0 ), binsAndFrequencyToCalculate->GetFrequency( bin, 0 ) ) ); - } - } - - return returnedHistogramMap; - } - - const ImageStatisticsCalculator::HistogramType * - ImageStatisticsCalculator::GetHistogram( unsigned int timeStep, unsigned int label ) const - { - if ( m_Image.IsNull() || (timeStep >= m_Image->GetTimeSteps()) ) - { - return nullptr; - } - - switch ( m_MaskingMode ) - { - case MASKING_MODE_NONE: - default: - { - if(m_DoIgnorePixelValue) - return m_MaskedImageHistogramVector[timeStep][label]; - - return m_ImageHistogramVector[timeStep][label]; - } - - case MASKING_MODE_IMAGE: - return m_MaskedImageHistogramVector[timeStep][label]; - - case MASKING_MODE_PLANARFIGURE: - return m_PlanarFigureHistogramVector[timeStep][label]; - } - } - - const ImageStatisticsCalculator::HistogramContainer & - ImageStatisticsCalculator::GetHistogramVector( unsigned int timeStep ) const - { - if ( m_Image.IsNull() || (timeStep >= m_Image->GetTimeSteps()) ) - { - return m_EmptyHistogramContainer; - } - - switch ( m_MaskingMode ) - { - case MASKING_MODE_NONE: - default: - { - if(m_DoIgnorePixelValue) - return m_MaskedImageHistogramVector[timeStep]; - - return m_ImageHistogramVector[timeStep]; - } - - case MASKING_MODE_IMAGE: - return m_MaskedImageHistogramVector[timeStep]; - - case MASKING_MODE_PLANARFIGURE: - return m_PlanarFigureHistogramVector[timeStep]; - } - } - - - const ImageStatisticsCalculator::Statistics & - ImageStatisticsCalculator::GetStatistics( unsigned int timeStep, unsigned int label ) const - { - if ( m_Image.IsNull() || (timeStep >= m_Image->GetTimeSteps()) ) - { - return m_EmptyStatistics; - } - - switch ( m_MaskingMode ) - { - case MASKING_MODE_NONE: - default: - { - if(m_DoIgnorePixelValue) - return m_MaskedImageStatisticsVector[timeStep][label]; - - return m_ImageStatisticsVector[timeStep][label]; - } - case MASKING_MODE_IMAGE: - return m_MaskedImageStatisticsVector[timeStep][label]; - - case MASKING_MODE_PLANARFIGURE: - return m_PlanarFigureStatisticsVector[timeStep][label]; - } - } - - const ImageStatisticsCalculator::StatisticsContainer & - ImageStatisticsCalculator::GetStatisticsVector( unsigned int timeStep ) const - { - if ( m_Image.IsNull() || (timeStep >= m_Image->GetTimeSteps()) ) - { - return m_EmptyStatisticsContainer; - } - - switch ( m_MaskingMode ) - { - case MASKING_MODE_NONE: - default: - { - if(m_DoIgnorePixelValue) - return m_MaskedImageStatisticsVector[timeStep]; - - return m_ImageStatisticsVector[timeStep]; - } - case MASKING_MODE_IMAGE: - return m_MaskedImageStatisticsVector[timeStep]; - - case MASKING_MODE_PLANARFIGURE: - return m_PlanarFigureStatisticsVector[timeStep]; - } - } - - - - void ImageStatisticsCalculator::ExtractImageAndMask( unsigned int timeStep ) - { - if ( m_Image.IsNull() ) - { - throw std::runtime_error( "Error: image empty!" ); - } - - if ( timeStep >= m_Image->GetTimeSteps() ) - { - throw std::runtime_error( "Error: invalid time step!" ); - } - - ImageTimeSelector::Pointer imageTimeSelector = ImageTimeSelector::New(); - imageTimeSelector->SetInput( m_Image ); - imageTimeSelector->SetTimeNr( timeStep ); - imageTimeSelector->UpdateLargestPossibleRegion(); - mitk::Image *timeSliceImage = imageTimeSelector->GetOutput(); - - - switch ( m_MaskingMode ) - { - case MASKING_MODE_NONE: - { - m_InternalImage = timeSliceImage; - m_InternalImageMask2D = nullptr; - m_InternalImageMask3D = nullptr; - if(m_DoIgnorePixelValue) - { - if( m_InternalImage->GetDimension() == 3 ) - { - if(itk::ImageIOBase::USHORT != timeSliceImage->GetPixelType().GetComponentType()) - CastToItkImage( timeSliceImage, m_InternalImageMask3D ); - else - CastToItkImage( timeSliceImage->Clone(), m_InternalImageMask3D ); - m_InternalImageMask3D->FillBuffer(1); - } - if( m_InternalImage->GetDimension() == 2 ) - { - if(itk::ImageIOBase::USHORT != timeSliceImage->GetPixelType().GetComponentType()) - CastToItkImage( timeSliceImage, m_InternalImageMask2D ); - else - CastToItkImage( timeSliceImage->Clone(), m_InternalImageMask2D ); - m_InternalImageMask2D->FillBuffer(1); - } - } - break; - } - - case MASKING_MODE_IMAGE: - { - if ( m_ImageMask.IsNotNull() && (m_ImageMask->GetReferenceCount() > 1) ) - { - if ( timeStep >= m_ImageMask->GetTimeSteps() ) - { - // Use the last mask time step in case the current time step is bigger than the total - // number of mask time steps. - // It makes more sense setting this to the last mask time step than to 0. - // For instance if you have a mask with 2 time steps and an image with 5: - // If time step 0 is selected, the mask will use time step 0. - // If time step 1 is selected, the mask will use time step 1. - // If time step 2+ is selected, the mask will use time step 1. - // If you have a mask with only one time step instead, this will always default to 0. - timeStep = m_ImageMask->GetTimeSteps() - 1; - } - - ImageTimeSelector::Pointer maskedImageTimeSelector = ImageTimeSelector::New(); - maskedImageTimeSelector->SetInput( m_ImageMask ); - maskedImageTimeSelector->SetTimeNr( timeStep ); - maskedImageTimeSelector->UpdateLargestPossibleRegion(); - mitk::Image *timeSliceMaskedImage = maskedImageTimeSelector->GetOutput(); - - m_InternalImage = timeSliceImage; - CastToItkImage( timeSliceMaskedImage, m_InternalImageMask3D ); - } - else - { - throw std::runtime_error( "Error: image mask empty!" ); - } - break; - } - case MASKING_MODE_PLANARFIGURE: - { - m_InternalImageMask2D = nullptr; - - if ( m_PlanarFigure.IsNull() ) - { - throw std::runtime_error( "Error: planar figure empty!" ); - } - if ( !m_PlanarFigure->IsClosed() ) - { - throw std::runtime_error( "Masking not possible for non-closed figures" ); - } - - const BaseGeometry *imageGeometry = timeSliceImage->GetGeometry(); - if ( imageGeometry == nullptr ) - { - throw std::runtime_error( "Image geometry invalid!" ); - } - - const PlaneGeometry *planarFigurePlaneGeometry = m_PlanarFigure->GetPlaneGeometry(); - if ( planarFigurePlaneGeometry == nullptr ) - { - throw std::runtime_error( "Planar-Figure not yet initialized!" ); - } - - const PlaneGeometry *planarFigureGeometry = - dynamic_cast< const PlaneGeometry * >( planarFigurePlaneGeometry ); - if ( planarFigureGeometry == nullptr ) - { - throw std::runtime_error( "Non-planar planar figures not supported!" ); - } - - // Find principal direction of PlanarFigure in input image - unsigned int axis; - if ( !this->GetPrincipalAxis( imageGeometry, - planarFigureGeometry->GetNormal(), axis ) ) - { - throw std::runtime_error( "Non-aligned planar figures not supported!" ); - } - m_PlanarFigureAxis = axis; - - // Find slice number corresponding to PlanarFigure in input image - MaskImage3DType::IndexType index; - imageGeometry->WorldToIndex( planarFigureGeometry->GetOrigin(), index ); - - unsigned int slice = index[axis]; - m_PlanarFigureSlice = slice; - - - // Extract slice with given position and direction from image - unsigned int dimension = timeSliceImage->GetDimension(); - - if (dimension != 2) - { - ExtractImageFilter::Pointer imageExtractor = ExtractImageFilter::New(); - imageExtractor->SetInput( timeSliceImage ); - imageExtractor->SetSliceDimension( axis ); - imageExtractor->SetSliceIndex( slice ); - imageExtractor->Update(); - m_InternalImage = imageExtractor->GetOutput(); - } - else - { - m_InternalImage = timeSliceImage; - } - - // Compute mask from PlanarFigure - AccessFixedDimensionByItk_1( - m_InternalImage, - InternalCalculateMaskFromPlanarFigure, - 2, axis ); - } - } - - if(m_DoIgnorePixelValue) - { - if ( m_InternalImage->GetDimension() == 3 ) - { - AccessFixedDimensionByItk_1( - m_InternalImage, - InternalMaskIgnoredPixels, - 3, - m_InternalImageMask3D.GetPointer() ); - } - else if ( m_InternalImage->GetDimension() == 2 ) - { - AccessFixedDimensionByItk_1( - m_InternalImage, - InternalMaskIgnoredPixels, - 2, - m_InternalImageMask2D.GetPointer() ); - } - } - } - - - bool ImageStatisticsCalculator::GetPrincipalAxis( - const BaseGeometry *geometry, Vector3D vector, - unsigned int &axis ) - { - vector.Normalize(); - for ( unsigned int i = 0; i < 3; ++i ) - { - Vector3D axisVector = geometry->GetAxisVector( i ); - axisVector.Normalize(); - - if ( fabs( fabs( axisVector * vector ) - 1.0) < mitk::eps ) - { - axis = i; - return true; - } - } - - return false; - } - - - unsigned int ImageStatisticsCalculator::calcNumberOfBins(mitk::ScalarType min, mitk::ScalarType max) - { - return std::ceil( ( (max - min ) / m_HistogramBinSize) ); - } - - - template < typename TPixel, unsigned int VImageDimension > - void ImageStatisticsCalculator::InternalCalculateStatisticsUnmasked( - const itk::Image< TPixel, VImageDimension > *image, - StatisticsContainer *statisticsContainer, - HistogramContainer* histogramContainer ) - { - typedef itk::Image< TPixel, VImageDimension > ImageType; - typedef typename ImageType::IndexType IndexType; - - typedef itk::Statistics::ScalarImageToHistogramGenerator< ImageType > - HistogramGeneratorType; - - statisticsContainer->clear(); - histogramContainer->clear(); - - // Progress listening... - typedef itk::SimpleMemberCommand< ImageStatisticsCalculator > ITKCommandType; - ITKCommandType::Pointer progressListener; - progressListener = ITKCommandType::New(); - progressListener->SetCallbackFunction( this, - &ImageStatisticsCalculator::UnmaskedStatisticsProgressUpdate ); - - - // Issue 100 artificial progress events since ScalarIMageToHistogramGenerator - // does not (yet?) support progress reporting - this->InvokeEvent( itk::StartEvent() ); - for ( unsigned int i = 0; i < 100; ++i ) - { - this->UnmaskedStatisticsProgressUpdate(); - } - - // Calculate statistics (separate filter) - typedef itk::ExtendedStatisticsImageFilter< ImageType > StatisticsFilterType; - typename StatisticsFilterType::Pointer statisticsFilter = StatisticsFilterType::New(); - statisticsFilter->SetInput( image ); - statisticsFilter->SetBinSize( 100 ); - statisticsFilter->SetCoordinateTolerance( 0.001 ); - statisticsFilter->SetDirectionTolerance( 0.001 ); - - unsigned long observerTag = statisticsFilter->AddObserver( itk::ProgressEvent(), progressListener ); - try - { - statisticsFilter->Update(); - } - catch (const itk::ExceptionObject& e) - { - mitkThrow() << "Image statistics calculation failed due to following ITK Exception: \n " << e.what(); - } - catch( const std::exception& e ) - { - //mitkThrow() << "Image statistics calculation failed due to following ITK Exception: \n " << e.what(); - } - - statisticsFilter->RemoveObserver( observerTag ); - this->InvokeEvent( itk::EndEvent() ); - - // Calculate minimum and maximum - typedef itk::MinimumMaximumImageCalculator< ImageType > MinMaxFilterType; - typename MinMaxFilterType::Pointer minMaxFilter = MinMaxFilterType::New(); - minMaxFilter->SetImage( image ); - unsigned long observerTag2 = minMaxFilter->AddObserver( itk::ProgressEvent(), progressListener ); - minMaxFilter->Compute(); - minMaxFilter->RemoveObserver( observerTag2 ); - this->InvokeEvent( itk::EndEvent() ); - - Statistics statistics; - statistics.Reset(); - statistics.SetLabel(1); - statistics.SetN(image->GetBufferedRegion().GetNumberOfPixels()); - statistics.SetMin(statisticsFilter->GetMinimum()); - statistics.SetMax(statisticsFilter->GetMaximum()); - statistics.SetMean(statisticsFilter->GetMean()); - statistics.SetMedian(statisticsFilter->GetMedian()); - statistics.SetVariance(statisticsFilter->GetVariance()); - statistics.SetSkewness(statisticsFilter->GetSkewness()); - statistics.SetKurtosis(statisticsFilter->GetKurtosis()); - statistics.SetUniformity( statisticsFilter->GetUniformity()); - statistics.SetEntropy( statisticsFilter->GetEntropy()); - statistics.SetUPP( statisticsFilter->GetUPP()); - statistics.SetMPP( statisticsFilter->GetMPP()); - statistics.SetSigma(statisticsFilter->GetSigma()); - statistics.SetRMS(sqrt( statistics.GetMean() * statistics.GetMean() + statistics.GetSigma() * statistics.GetSigma() )); - - statistics.GetMinIndex().set_size(image->GetImageDimension()); - statistics.GetMaxIndex().set_size(image->GetImageDimension()); - - vnl_vector tmpMaxIndex; - vnl_vector tmpMinIndex; - - tmpMaxIndex.set_size(image->GetImageDimension() ); - tmpMinIndex.set_size(image->GetImageDimension() ); - - for (unsigned int i=0; iGetIndexOfMaximum()[i]; - tmpMinIndex[i] = minMaxFilter->GetIndexOfMinimum()[i]; - } - - statistics.SetMinIndex(tmpMaxIndex); - statistics.SetMinIndex(tmpMinIndex); - - if( IsHotspotCalculated() && VImageDimension == 3 ) - { - typedef itk::Image< unsigned short, VImageDimension > MaskImageType; - typename MaskImageType::Pointer nullMask; - bool isHotspotDefined(false); - Statistics hotspotStatistics = this->CalculateHotspotStatistics(image, nullMask.GetPointer(), m_HotspotRadiusInMM, isHotspotDefined, 0 ); - if (isHotspotDefined) - { - statistics.SetHasHotspotStatistics(true); - statistics.GetHotspotStatistics() = hotspotStatistics; - } - else - { - statistics.SetHasHotspotStatistics(false); - } - - if(statistics.GetHotspotStatistics().HasHotspotStatistics() ) - { - MITK_DEBUG << "Hotspot statistics available"; - statistics.SetHotspotIndex(hotspotStatistics.GetHotspotIndex()); - } - else - { - MITK_ERROR << "No hotspot statistics available!"; - } - } - - statisticsContainer->push_back( statistics ); - - // Calculate histogram - // calculate bin size or number of bins - unsigned int numberOfBins = 200; // default number of bins - if (m_UseDefaultBinSize) - { - m_HistogramBinSize = std::ceil( (statistics.GetMax() - statistics.GetMin() + 1)/numberOfBins ); - } - else - { - numberOfBins = calcNumberOfBins(statistics.GetMin(), statistics.GetMax()); - } - typename HistogramGeneratorType::Pointer histogramGenerator = HistogramGeneratorType::New(); - histogramGenerator->SetInput( image ); - histogramGenerator->SetMarginalScale( 100 ); - histogramGenerator->SetNumberOfBins( numberOfBins ); - histogramGenerator->SetHistogramMin( statistics.GetMin() ); - histogramGenerator->SetHistogramMax( statistics.GetMax() ); - histogramGenerator->Compute(); - histogramContainer->push_back( histogramGenerator->GetOutput() ); - } - - template < typename TPixel, unsigned int VImageDimension > - void ImageStatisticsCalculator::InternalMaskIgnoredPixels( - const itk::Image< TPixel, VImageDimension > *image, - itk::Image< unsigned short, VImageDimension > *maskImage ) - { - typedef itk::Image< TPixel, VImageDimension > ImageType; - typedef itk::Image< unsigned short, VImageDimension > MaskImageType; - - itk::ImageRegionIterator - itmask(maskImage, maskImage->GetLargestPossibleRegion()); - itk::ImageRegionConstIterator - itimage(image, image->GetLargestPossibleRegion()); - - itmask.GoToBegin(); - itimage.GoToBegin(); - - while( !itmask.IsAtEnd() ) - { - if(m_IgnorePixelValue == itimage.Get()) - { - itmask.Set(0); - } - - ++itmask; - ++itimage; - } - } - - template < typename TPixel, unsigned int VImageDimension > - void ImageStatisticsCalculator::InternalCalculateStatisticsMasked( - const itk::Image< TPixel, VImageDimension > *image, - itk::Image< unsigned short, VImageDimension > *maskImage, - StatisticsContainer* statisticsContainer, - HistogramContainer* histogramContainer ) - { - typedef itk::Image< TPixel, VImageDimension > ImageType; - typedef itk::Image< unsigned short, VImageDimension > MaskImageType; - typedef typename ImageType::IndexType IndexType; - typedef typename ImageType::PointType PointType; - typedef typename ImageType::SpacingType SpacingType; - typedef typename ImageType::Pointer ImagePointer; - typedef itk::ExtendedLabelStatisticsImageFilter< ImageType, MaskImageType > LabelStatisticsFilterType; - typedef itk::ChangeInformationImageFilter< MaskImageType > ChangeInformationFilterType; - typedef itk::ExtractImageFilter< ImageType, ImageType > ExtractImageFilterType; - - statisticsContainer->clear(); - histogramContainer->clear(); - - // Make sure that mask is set - if ( maskImage == nullptr ) - { - itkExceptionMacro( << "Mask image needs to be set!" ); - } - - // Make sure that spacing of mask and image are the same - //SpacingType imageSpacing = image->GetSpacing(); - //SpacingType maskSpacing = maskImage->GetSpacing(); - //PointType zeroPoint; zeroPoint.Fill( 0.0 ); - //if ( (zeroPoint + imageSpacing).SquaredEuclideanDistanceTo( (zeroPoint + maskSpacing) ) > mitk::eps ) - //{ - // itkExceptionMacro( << "Mask needs to have same spacing as image! (Image spacing: " << imageSpacing << "; Mask spacing: " << maskSpacing << ")" ); - //} - // Make sure that orientation of mask and image are the same - typedef typename ImageType::DirectionType DirectionType; - DirectionType imageDirection = image->GetDirection(); - DirectionType maskDirection = maskImage->GetDirection(); - for( int i = 0; i < imageDirection.ColumnDimensions; ++i ) - { - for( int j = 0; j < imageDirection.ColumnDimensions; ++j ) - { - double differenceDirection = imageDirection[i][j] - maskDirection[i][j]; - if ( fabs( differenceDirection ) > mitk::eps ) - { - double differenceDirection = imageDirection[i][j] - maskDirection[i][j]; - if ( fabs( differenceDirection ) > 0.001 /*mitk::eps*/ ) // TODO: temp fix (bug 17121) - { - itkExceptionMacro( << "Mask needs to have same direction as image! (Image direction: " << imageDirection << "; Mask direction: " << maskDirection << ")" ); - } - } - } - } - // Make sure that the voxels of mask and image are correctly "aligned", i.e., voxel boundaries are the same in both images - PointType imageOrigin = image->GetOrigin(); - PointType maskOrigin = maskImage->GetOrigin(); - long offset[ImageType::ImageDimension]; - - typedef itk::ContinuousIndex ContinousIndexType; - ContinousIndexType maskOriginContinousIndex, imageOriginContinousIndex; - - image->TransformPhysicalPointToContinuousIndex(maskOrigin, maskOriginContinousIndex); - image->TransformPhysicalPointToContinuousIndex(imageOrigin, imageOriginContinousIndex); - - for ( unsigned int i = 0; i < ImageType::ImageDimension; ++i ) - { - double misalignment = maskOriginContinousIndex[i] - floor( maskOriginContinousIndex[i] + 0.5 ); - if ( fabs( misalignment ) > mitk::eps ) - { - itkWarningMacro( << "Pixels/voxels of mask and image are not sufficiently aligned! (Misalignment: " << misalignment << ")" ); - } - - double indexCoordDistance = maskOriginContinousIndex[i] - imageOriginContinousIndex[i]; - offset[i] = int( indexCoordDistance + image->GetBufferedRegion().GetIndex()[i] + 0.5 ); - } - - // Adapt the origin and region (index/size) of the mask so that the origin of both are the same - typename ChangeInformationFilterType::Pointer adaptMaskFilter; - adaptMaskFilter = ChangeInformationFilterType::New(); - adaptMaskFilter->ChangeOriginOn(); - adaptMaskFilter->ChangeRegionOn(); - adaptMaskFilter->SetInput( maskImage ); - adaptMaskFilter->SetOutputOrigin( image->GetOrigin() ); - adaptMaskFilter->SetOutputOffset( offset ); - adaptMaskFilter->SetCoordinateTolerance( 0.001 ); - adaptMaskFilter->SetDirectionTolerance( 0.001 ); - - - typename MaskImageType::Pointer adaptedMaskImage; - try - { - adaptMaskFilter->Update(); - adaptedMaskImage = adaptMaskFilter->GetOutput(); - } - catch( const itk::ExceptionObject &e) - { - mitkThrow() << "Attempt to adapt shifted origin of the mask image failed due to ITK Exception: \n" << e.what(); - } - catch( const std::exception& e ) - { - //mitkThrow() << "Image statistics calculation failed due to following ITK Exception: \n " << e.what(); - } - - - // Make sure that mask region is contained within image region - if ( adaptedMaskImage.IsNotNull() && - !image->GetLargestPossibleRegion().IsInside( adaptedMaskImage->GetLargestPossibleRegion() ) ) - { - itkWarningMacro( << "Mask region needs to be inside of image region! (Image region: " - << image->GetLargestPossibleRegion() << "; Mask region: " << adaptedMaskImage->GetLargestPossibleRegion() << ")" ); - } - - - // If mask region is smaller than image region, extract the sub-sampled region from the original image - typename ImageType::SizeType imageSize = image->GetBufferedRegion().GetSize(); - typename ImageType::SizeType maskSize = maskImage->GetBufferedRegion().GetSize(); - bool maskSmallerImage = false; - for ( unsigned int i = 0; i < ImageType::ImageDimension; ++i ) - { - if ( maskSize[i] < imageSize[i] ) - { - maskSmallerImage = true; - } - } - - typename ImageType::ConstPointer adaptedImage; - if ( maskSmallerImage ) - { - typename ExtractImageFilterType::Pointer extractImageFilter = ExtractImageFilterType::New(); - extractImageFilter->SetInput( image ); - extractImageFilter->SetExtractionRegion( adaptedMaskImage->GetBufferedRegion() ); - extractImageFilter->SetCoordinateTolerance( 0.001 ); - extractImageFilter->SetDirectionTolerance( 0.001 ); - extractImageFilter->Update(); - adaptedImage = extractImageFilter->GetOutput(); - } - else - { - adaptedImage = image; - } - - // Initialize Filter - typedef itk::StatisticsImageFilter< ImageType > StatisticsFilterType; - typename StatisticsFilterType::Pointer statisticsFilter = StatisticsFilterType::New(); - statisticsFilter->SetInput( adaptedImage ); - - try - { - statisticsFilter->Update(); - } - catch( const itk::ExceptionObject& e) - { - mitkThrow() << "Image statistics initialization computation failed with ITK Exception: \n " << e.what(); - } - catch( const std::exception& e ) - { - //mitkThrow() << "Image statistics calculation failed due to following ITK Exception: \n " << e.what(); - } - - // Calculate bin size or number of bins - unsigned int numberOfBins = m_HistogramBinSize; // default number of bins - double maximum = 0.0; - double minimum = 0.0; - - if (m_UseBinSizeBasedOnVOIRegion) - { - maximum = statisticsFilter->GetMaximum(); - minimum = statisticsFilter->GetMinimum(); - - if (m_UseDefaultBinSize) - { - m_HistogramBinSize = std::ceil( static_cast((statisticsFilter->GetMaximum() - statisticsFilter->GetMinimum() + 1)/numberOfBins) ); - } - else - { - numberOfBins = calcNumberOfBins(statisticsFilter->GetMinimum(), statisticsFilter->GetMaximum()); - } - } - else - { - double sig = 0.0; - int counter = 0; - - //Find the min and max values for the Roi to set the range for the histogram - GetMinAndMaxValue( minimum, maximum, counter, sig, image, maskImage); -// numberOfBins = maximum - minimum; -// if(maximum - minimum <= 10) -// { -// numberOfBins = 100; -// } - } - - - typename LabelStatisticsFilterType::Pointer labelStatisticsFilter = LabelStatisticsFilterType::New(); - labelStatisticsFilter->SetInput( adaptedImage ); - labelStatisticsFilter->SetLabelInput( adaptedMaskImage ); - labelStatisticsFilter->SetCoordinateTolerance( 0.001 ); - labelStatisticsFilter->SetDirectionTolerance( 0.001 ); - labelStatisticsFilter->UseHistogramsOn(); - labelStatisticsFilter->SetHistogramParameters( numberOfBins, floor(minimum), ceil(maximum) ); //statisticsFilter->GetMinimum() statisticsFilter->GetMaximum() - - // Add progress listening - typedef itk::SimpleMemberCommand< ImageStatisticsCalculator > ITKCommandType; - ITKCommandType::Pointer progressListener; - progressListener = ITKCommandType::New(); - progressListener->SetCallbackFunction( this, - &ImageStatisticsCalculator::MaskedStatisticsProgressUpdate ); - unsigned long observerTag = labelStatisticsFilter->AddObserver( - itk::ProgressEvent(), progressListener ); - - // Execute filter - this->InvokeEvent( itk::StartEvent() ); - - // Make sure that only the mask region is considered (otherwise, if the mask region is smaller - // than the image region, the Update() would result in an exception). - labelStatisticsFilter->GetOutput()->SetRequestedRegion( adaptedMaskImage->GetLargestPossibleRegion() ); - - // Execute the filter - try - { - labelStatisticsFilter->Update(); - } - catch( const itk::ExceptionObject& e) - { - mitkThrow() << "Image statistics calculation failed due to following ITK Exception: \n " << e.what(); - } - catch( const std::exception& e ) - { - //mitkThrow() << "Image statistics calculation failed due to following ITK Exception: \n " << e.what(); - } - - this->InvokeEvent( itk::EndEvent() ); - - if( observerTag ) - labelStatisticsFilter->RemoveObserver( observerTag ); - - // Find all relevant labels of mask (other than 0) - std::list< int > relevantLabels = labelStatisticsFilter->GetRelevantLabels(); - unsigned int i; - - if ( labelStatisticsFilter->GetMaskingNonEmpty() ) - { - std::list< int >::iterator it; - for ( it = relevantLabels.begin(), i = 0; - it != relevantLabels.end(); - ++it, ++i ) - { - Statistics statistics; // restore previous code - labelStatisticsFilter->GetHistogram(*it) ; - histogramContainer->push_back( HistogramType::ConstPointer( labelStatisticsFilter->GetHistogram( (*it) ) ) ); - - statistics.SetLabel (*it); - statistics.SetN(labelStatisticsFilter->GetCount( *it )); - statistics.SetMin(labelStatisticsFilter->GetMinimum( *it )); - statistics.SetMax(labelStatisticsFilter->GetMaximum( *it )); - statistics.SetMean(labelStatisticsFilter->GetMean( *it )); - statistics.SetMedian(labelStatisticsFilter->GetMedian( *it)); - statistics.SetMedian(labelStatisticsFilter->GetMedian( *it )); - statistics.SetVariance(labelStatisticsFilter->GetVariance( *it )); - statistics.SetSigma(labelStatisticsFilter->GetSigma( *it )); - statistics.SetSkewness(labelStatisticsFilter->GetSkewness( *it )); - statistics.SetKurtosis(labelStatisticsFilter->GetKurtosis( *it )); - statistics.SetUniformity( labelStatisticsFilter->GetUniformity( *it )); - statistics.SetEntropy( labelStatisticsFilter->GetEntropy( *it )); - statistics.SetUPP( labelStatisticsFilter->GetUPP( *it)); - statistics.SetMPP( labelStatisticsFilter->GetMPP( *it)); - statistics.SetRMS(sqrt( statistics.GetMean() * statistics.GetMean() - + statistics.GetSigma() * statistics.GetSigma() )); - - // restrict image to mask area for min/max index calculation - typedef itk::MaskImageFilter< ImageType, MaskImageType, ImageType > MaskImageFilterType; - typename MaskImageFilterType::Pointer masker = MaskImageFilterType::New(); - bool isMinAndMaxSameValue = (statistics.GetMin() == statistics.GetMax()); - // bug 17962: following is a workaround for the case when min and max are the same, we can probably find a nicer way here - double outsideValue = (isMinAndMaxSameValue ? (statistics.GetMax()/2) : (statistics.GetMin()+statistics.GetMax())/2); - masker->SetOutsideValue( outsideValue ); - masker->SetInput1(adaptedImage); - masker->SetInput2(adaptedMaskImage); - masker->SetCoordinateTolerance( 0.001 ); - masker->SetDirectionTolerance( 0.001 ); - masker->Update(); - // get index of minimum and maximum - typedef itk::MinimumMaximumImageCalculator< ImageType > MinMaxFilterType; - typename MinMaxFilterType::Pointer minMaxFilter = MinMaxFilterType::New(); - minMaxFilter->SetImage( masker->GetOutput() ); - unsigned long observerTag2 = minMaxFilter->AddObserver( itk::ProgressEvent(), progressListener ); - minMaxFilter->Compute(); - minMaxFilter->RemoveObserver( observerTag2 ); - this->InvokeEvent( itk::EndEvent() ); - - typename MinMaxFilterType::IndexType tempMaxIndex = minMaxFilter->GetIndexOfMaximum(); - // bug 17962: following is a workaround for the case when min and max are the same, we can probably find a nicer way here - typename MinMaxFilterType::IndexType tempMinIndex = - (isMinAndMaxSameValue ? minMaxFilter->GetIndexOfMaximum() : minMaxFilter->GetIndexOfMinimum()); - - // FIX BUG 14644 - //If a PlanarFigure is used for segmentation the - //adaptedImage is a single slice (2D). Adding the - // 3. dimension. - - vnl_vector maxIndex; - vnl_vector minIndex; - maxIndex.set_size(m_Image->GetDimension()); - minIndex.set_size(m_Image->GetDimension()); - - if (m_MaskingMode == MASKING_MODE_PLANARFIGURE && m_Image->GetDimension()==3) - { - maxIndex[m_PlanarFigureCoordinate0] = tempMaxIndex[0]; - maxIndex[m_PlanarFigureCoordinate1] = tempMaxIndex[1]; - maxIndex[m_PlanarFigureAxis] = m_PlanarFigureSlice; - - minIndex[m_PlanarFigureCoordinate0] = tempMinIndex[0] ; - minIndex[m_PlanarFigureCoordinate1] = tempMinIndex[1]; - minIndex[m_PlanarFigureAxis] = m_PlanarFigureSlice; - } else - { - for (unsigned int i = 0; ipush_back( statistics ); - } - } - else - { - histogramContainer->push_back( HistogramType::ConstPointer( m_EmptyHistogram ) ); - statisticsContainer->push_back( Statistics() ); - } - } - - - template - ImageStatisticsCalculator::ImageExtrema - ImageStatisticsCalculator::CalculateExtremaWorld( - const itk::Image *inputImage, - itk::Image *maskImage, - double neccessaryDistanceToImageBorderInMM, - unsigned int label) - { - typedef itk::Image< TPixel, VImageDimension > ImageType; - typedef itk::Image< unsigned short, VImageDimension > MaskImageType; - - typedef itk::ImageRegionConstIteratorWithIndex MaskImageIteratorType; - typedef itk::ImageRegionConstIteratorWithIndex InputImageIndexIteratorType; - - typename ImageType::SpacingType spacing = inputImage->GetSpacing(); - - ImageExtrema minMax; - minMax.Defined = false; - minMax.MaxIndex.set_size(VImageDimension); - minMax.MaxIndex.set_size(VImageDimension); - - typename ImageType::RegionType allowedExtremaRegion = inputImage->GetLargestPossibleRegion(); - - bool keepDistanceToImageBorders( neccessaryDistanceToImageBorderInMM > 0 ); - if (keepDistanceToImageBorders) - { - long distanceInPixels[VImageDimension]; - for(unsigned short dimension = 0; dimension < VImageDimension; ++dimension) - { - // To confirm that the whole hotspot is inside the image we have to keep a specific distance to the image-borders, which is as long as - // the radius. To get the amount of indices we divide the radius by spacing and add 0.5 because voxels are center based: - // For example with a radius of 2.2 and a spacing of 1 two indices are enough because 2.2 / 1 + 0.5 = 2.7 => 2. - // But with a radius of 2.7 we need 3 indices because 2.7 / 1 + 0.5 = 3.2 => 3 - distanceInPixels[dimension] = int( neccessaryDistanceToImageBorderInMM / spacing[dimension] + 0.5); - } - - allowedExtremaRegion.ShrinkByRadius(distanceInPixels); - } - - InputImageIndexIteratorType imageIndexIt(inputImage, allowedExtremaRegion); - - float maxValue = itk::NumericTraits::min(); - float minValue = itk::NumericTraits::max(); - - typename ImageType::IndexType maxIndex; - typename ImageType::IndexType minIndex; - - for(unsigned short i = 0; i < VImageDimension; ++i) - { - maxIndex[i] = 0; - minIndex[i] = 0; - } - - if (maskImage != nullptr) - { - MaskImageIteratorType maskIt(maskImage, maskImage->GetLargestPossibleRegion()); - typename ImageType::IndexType imageIndex; - typename ImageType::PointType worldPosition; - typename ImageType::IndexType maskIndex; - - for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt) - { - imageIndex = maskIndex = maskIt.GetIndex(); - - if(maskIt.Get() == label) - { - if( allowedExtremaRegion.IsInside(imageIndex) ) - { - imageIndexIt.SetIndex( imageIndex ); - double value = imageIndexIt.Get(); - minMax.Defined = true; - - //Calculate minimum, maximum and corresponding index-values - if( value > maxValue ) - { - maxIndex = imageIndexIt.GetIndex(); - maxValue = value; - } - - if(value < minValue ) - { - minIndex = imageIndexIt.GetIndex(); - minValue = value; - } - } - } - } - } - else - { - for(imageIndexIt.GoToBegin(); !imageIndexIt.IsAtEnd(); ++imageIndexIt) - { - double value = imageIndexIt.Get(); - minMax.Defined = true; - - //Calculate minimum, maximum and corresponding index-values - if( value > maxValue ) - { - maxIndex = imageIndexIt.GetIndex(); - maxValue = value; - } - - if(value < minValue ) - { - minIndex = imageIndexIt.GetIndex(); - minValue = value; - } - } - } - - minMax.MaxIndex.set_size(VImageDimension); - minMax.MinIndex.set_size(VImageDimension); - - for(unsigned int i = 0; i < minMax.MaxIndex.size(); ++i) - { - minMax.MaxIndex[i] = maxIndex[i]; - } - - for(unsigned int i = 0; i < minMax.MinIndex.size(); ++i) - { - minMax.MinIndex[i] = minIndex[i]; - } - - minMax.Max = maxValue; - minMax.Min = minValue; - - return minMax; - } - - template - itk::Size - ImageStatisticsCalculator - ::CalculateConvolutionKernelSize(double spacing[VImageDimension], double radiusInMM) - { - typedef itk::Image< float, VImageDimension > KernelImageType; - typedef typename KernelImageType::SizeType SizeType; - SizeType maskSize; - - for(unsigned int i = 0; i < VImageDimension; ++i) - { - maskSize[i] = static_cast( 2 * radiusInMM / spacing[i]); - - // We always want an uneven size to have a clear center point in the convolution mask - if(maskSize[i] % 2 == 0 ) - { - ++maskSize[i]; - } - } - return maskSize; - } - - template - itk::SmartPointer< itk::Image > - ImageStatisticsCalculator - ::GenerateHotspotSearchConvolutionKernel(double mmPerPixel[VImageDimension], double radiusInMM) - { - std::stringstream ss; - for (unsigned int i = 0; i < VImageDimension; ++i) - { - ss << mmPerPixel[i]; - if (i < VImageDimension -1) - ss << ","; - } - MITK_DEBUG << "Update convolution kernel for spacing (" << ss.str() << ") and radius " << radiusInMM << "mm"; - - - double radiusInMMSquared = radiusInMM * radiusInMM; - typedef itk::Image< float, VImageDimension > KernelImageType; - typename KernelImageType::Pointer convolutionKernel = KernelImageType::New(); - - // Calculate size and allocate mask image - typedef typename KernelImageType::SizeType SizeType; - SizeType maskSize = this->CalculateConvolutionKernelSize(mmPerPixel, radiusInMM); - - Point3D convolutionMaskCenterIndex; convolutionMaskCenterIndex.Fill(0.0); - for(unsigned int i = 0; i < VImageDimension; ++i) - { - convolutionMaskCenterIndex[i] = 0.5 * (double)(maskSize[i]-1); - } - - typedef typename KernelImageType::IndexType IndexType; - IndexType maskIndex; - maskIndex.Fill(0); - - typedef typename KernelImageType::RegionType RegionType; - RegionType maskRegion; - maskRegion.SetSize(maskSize); - maskRegion.SetIndex(maskIndex); - - convolutionKernel->SetRegions(maskRegion); - convolutionKernel->SetSpacing(mmPerPixel); - convolutionKernel->Allocate(); - - // Fill mask image values by subsampling the image grid - typedef itk::ImageRegionIteratorWithIndex MaskIteratorType; - MaskIteratorType maskIt(convolutionKernel,maskRegion); - - int numberOfSubVoxelsPerDimension = 2; // per dimension! - int numberOfSubVoxels = ::pow( static_cast(numberOfSubVoxelsPerDimension), static_cast(VImageDimension) ); - double subVoxelSizeInPixels = 1.0 / (double)numberOfSubVoxelsPerDimension; - double valueOfOneSubVoxel = 1.0 / (double)numberOfSubVoxels; - double maskValue = 0.0; - Point3D subVoxelIndexPosition; - double distanceSquared = 0.0; - - typedef itk::ContinuousIndex ContinuousIndexType; - for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt) - { - ContinuousIndexType indexPoint(maskIt.GetIndex()); - Point3D voxelPosition; - for (unsigned int dimension = 0; dimension < VImageDimension; ++dimension) - { - voxelPosition[dimension] = indexPoint[dimension]; - } - - maskValue = 0.0; - Vector3D subVoxelOffset; subVoxelOffset.Fill(0.0); - // iterate sub-voxels by iterating all possible offsets - for (subVoxelOffset[0] = -0.5 + subVoxelSizeInPixels / 2.0; - subVoxelOffset[0] < +0.5; - subVoxelOffset[0] += subVoxelSizeInPixels) - { - for (subVoxelOffset[1] = -0.5 + subVoxelSizeInPixels / 2.0; - subVoxelOffset[1] < +0.5; - subVoxelOffset[1] += subVoxelSizeInPixels) - { - for (subVoxelOffset[2] = -0.5 + subVoxelSizeInPixels / 2.0; - subVoxelOffset[2] < +0.5; - subVoxelOffset[2] += subVoxelSizeInPixels) - { - subVoxelIndexPosition = voxelPosition + subVoxelOffset; // this COULD be integrated into the for-loops if neccessary (add voxelPosition to initializer and end condition) - distanceSquared = - (subVoxelIndexPosition[0]-convolutionMaskCenterIndex[0]) * mmPerPixel[0] * (subVoxelIndexPosition[0]-convolutionMaskCenterIndex[0]) * mmPerPixel[0] - + (subVoxelIndexPosition[1]-convolutionMaskCenterIndex[1]) * mmPerPixel[1] * (subVoxelIndexPosition[1]-convolutionMaskCenterIndex[1]) * mmPerPixel[1] - + (subVoxelIndexPosition[2]-convolutionMaskCenterIndex[2]) * mmPerPixel[2] * (subVoxelIndexPosition[2]-convolutionMaskCenterIndex[2]) * mmPerPixel[2]; - - if (distanceSquared <= radiusInMMSquared) - { - maskValue += valueOfOneSubVoxel; - } - } - } - } - maskIt.Set( maskValue ); - } - - return convolutionKernel; - } - - template - itk::SmartPointer > - ImageStatisticsCalculator::GenerateConvolutionImage( const itk::Image* inputImage ) - { - double mmPerPixel[VImageDimension]; - for (unsigned int dimension = 0; dimension < VImageDimension; ++dimension) - { - mmPerPixel[dimension] = inputImage->GetSpacing()[dimension]; - } - - // update convolution kernel - typedef itk::Image< float, VImageDimension > KernelImageType; - typename KernelImageType::Pointer convolutionKernel = this->GenerateHotspotSearchConvolutionKernel(mmPerPixel, m_HotspotRadiusInMM); - - // update convolution image - typedef itk::Image< TPixel, VImageDimension > InputImageType; - typedef itk::Image< TPixel, VImageDimension > ConvolutionImageType; - typedef itk::FFTConvolutionImageFilter ConvolutionFilterType; - - typename ConvolutionFilterType::Pointer convolutionFilter = ConvolutionFilterType::New(); - typedef itk::ConstantBoundaryCondition BoundaryConditionType; - BoundaryConditionType boundaryCondition; - boundaryCondition.SetConstant(0.0); - - if (GetHotspotMustBeCompletlyInsideImage()) - { - // overwrite default boundary condition - convolutionFilter->SetBoundaryCondition(&boundaryCondition); - } - - convolutionFilter->SetInput(inputImage); - convolutionFilter->SetKernelImage(convolutionKernel); - convolutionFilter->SetNormalize(true); - MITK_DEBUG << "Update Convolution image for hotspot search"; - convolutionFilter->UpdateLargestPossibleRegion(); - - typename ConvolutionImageType::Pointer convolutionImage = convolutionFilter->GetOutput(); - convolutionImage->SetSpacing( inputImage->GetSpacing() ); // only workaround because convolution filter seems to ignore spacing of input image - - m_HotspotRadiusInMMChanged = false; - return convolutionImage; - } - - template < typename TPixel, unsigned int VImageDimension> - void - ImageStatisticsCalculator - ::FillHotspotMaskPixels( itk::Image* maskImage, - itk::Point sphereCenter, - double sphereRadiusInMM) - { - typedef itk::Image< TPixel, VImageDimension > MaskImageType; - typedef itk::ImageRegionIteratorWithIndex MaskImageIteratorType; - - MaskImageIteratorType maskIt(maskImage, maskImage->GetLargestPossibleRegion()); - - typename MaskImageType::IndexType maskIndex; - typename MaskImageType::PointType worldPosition; - - for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt) - { - maskIndex = maskIt.GetIndex(); - maskImage->TransformIndexToPhysicalPoint(maskIndex, worldPosition); - maskIt.Set( worldPosition.EuclideanDistanceTo(sphereCenter) <= sphereRadiusInMM ? 1 : 0 ); - } - } - - template < typename TPixel, unsigned int VImageDimension> - ImageStatisticsCalculator::Statistics - ImageStatisticsCalculator::CalculateHotspotStatistics( - const itk::Image* inputImage, - itk::Image* maskImage, - double radiusInMM, - bool& isHotspotDefined, - unsigned int label) - { - // get convolution image (updated in GenerateConvolutionImage()) - typedef itk::Image< TPixel, VImageDimension > InputImageType; - typedef itk::Image< TPixel, VImageDimension > ConvolutionImageType; - typedef itk::Image< float, VImageDimension > KernelImageType; - typedef itk::Image< unsigned short, VImageDimension > MaskImageType; - - //typename ConvolutionImageType::Pointer convolutionImage = dynamic_cast(this->GenerateConvolutionImage(inputImage)); - typename ConvolutionImageType::Pointer convolutionImage = this->GenerateConvolutionImage(inputImage); - - if (convolutionImage.IsNull()) - { - MITK_ERROR << "Empty convolution image in CalculateHotspotStatistics(). We should never reach this state (logic error)."; - throw std::logic_error("Empty convolution image in CalculateHotspotStatistics()"); - } - - // find maximum in convolution image, given the current mask - double requiredDistanceToBorder = m_HotspotMustBeCompletelyInsideImage ? m_HotspotRadiusInMM : -1.0; - ImageExtrema convolutionImageInformation = CalculateExtremaWorld(convolutionImage.GetPointer(), maskImage, requiredDistanceToBorder, label); - - isHotspotDefined = convolutionImageInformation.Defined; - - if (!isHotspotDefined) - { - m_EmptyStatistics.Reset(VImageDimension); - MITK_ERROR << "No origin of hotspot-sphere was calculated! Returning empty statistics"; - return m_EmptyStatistics; - } - else - { - // create a binary mask around the "hotspot" region, fill the shape of a sphere around our hotspot center - typedef itk::ImageDuplicator< InputImageType > DuplicatorType; - typename DuplicatorType::Pointer copyMachine = DuplicatorType::New(); - copyMachine->SetInputImage(inputImage); - copyMachine->Update(); - - typedef itk::CastImageFilter< InputImageType, MaskImageType > CastFilterType; - typename CastFilterType::Pointer caster = CastFilterType::New(); - caster->SetInput( copyMachine->GetOutput() ); - caster->Update(); - typename MaskImageType::Pointer hotspotMaskITK = caster->GetOutput(); - - typedef typename InputImageType::IndexType IndexType; - IndexType maskCenterIndex; - for (unsigned int d =0; d< VImageDimension;++d) maskCenterIndex[d]=convolutionImageInformation.MaxIndex[d]; - typename ConvolutionImageType::PointType maskCenter; - inputImage->TransformIndexToPhysicalPoint(maskCenterIndex,maskCenter); - - this->FillHotspotMaskPixels(hotspotMaskITK.GetPointer(), maskCenter, radiusInMM); - - // calculate statistics within the binary mask - typedef itk::ExtendedLabelStatisticsImageFilter< InputImageType, MaskImageType> LabelStatisticsFilterType; - typename LabelStatisticsFilterType::Pointer labelStatisticsFilter; - labelStatisticsFilter = LabelStatisticsFilterType::New(); - labelStatisticsFilter->SetInput( inputImage ); - labelStatisticsFilter->SetLabelInput( hotspotMaskITK ); - labelStatisticsFilter->SetCoordinateTolerance( 0.001 ); - labelStatisticsFilter->SetDirectionTolerance( 0.001 ); - - labelStatisticsFilter->Update(); - - Statistics hotspotStatistics; - hotspotStatistics.SetHotspotIndex(convolutionImageInformation.MaxIndex); - hotspotStatistics.SetMean(convolutionImageInformation.Max); - - if ( labelStatisticsFilter->HasLabel( 1 ) ) - { - hotspotStatistics.SetLabel (1); - hotspotStatistics.SetN(labelStatisticsFilter->GetCount(1)); - hotspotStatistics.SetMin(labelStatisticsFilter->GetMinimum(1)); - hotspotStatistics.SetMax(labelStatisticsFilter->GetMaximum(1)); - hotspotStatistics.SetMedian(labelStatisticsFilter->GetMedian(1)); - hotspotStatistics.SetVariance(labelStatisticsFilter->GetVariance(1)); - hotspotStatistics.SetSigma(labelStatisticsFilter->GetSigma(1)); - hotspotStatistics.SetRMS(sqrt( hotspotStatistics.GetMean() * hotspotStatistics.GetMean() - + hotspotStatistics.GetSigma() * hotspotStatistics.GetSigma() )); - - MITK_DEBUG << "Statistics for inside hotspot: Mean " << hotspotStatistics.GetMean() - << ", SD " << hotspotStatistics.GetSigma() - << ", Max " << hotspotStatistics.GetMax() - << ", Min " << hotspotStatistics.GetMin(); - } - else - { - MITK_ERROR << "Uh oh! Unable to calculate statistics for hotspot region..."; - return m_EmptyStatistics; - } - - return hotspotStatistics; - } - } - - template < typename TPixel, unsigned int VImageDimension > - void ImageStatisticsCalculator::InternalCalculateMaskFromPlanarFigure( - const itk::Image< TPixel, VImageDimension > *image, unsigned int axis ) - { - typedef itk::Image< TPixel, VImageDimension > ImageType; - - typedef itk::CastImageFilter< ImageType, MaskImage2DType > CastFilterType; - - // Generate mask image as new image with same header as input image and - // initialize with 1. - typename CastFilterType::Pointer castFilter = CastFilterType::New(); - castFilter->SetInput( image ); - castFilter->Update(); - castFilter->GetOutput()->FillBuffer( 1 ); - - // all PolylinePoints of the PlanarFigure are stored in a vtkPoints object. - // These points are used by the vtkLassoStencilSource to create - // a vtkImageStencil. - const mitk::PlaneGeometry *planarFigurePlaneGeometry = m_PlanarFigure->GetPlaneGeometry(); - const typename PlanarFigure::PolyLineType planarFigurePolyline = m_PlanarFigure->GetPolyLine( 0 ); - const mitk::BaseGeometry *imageGeometry3D = m_Image->GetGeometry( 0 ); - // If there is a second poly line in a closed planar figure, treat it as a hole. - PlanarFigure::PolyLineType planarFigureHolePolyline; - - if (m_PlanarFigure->GetPolyLinesSize() == 2) - planarFigureHolePolyline = m_PlanarFigure->GetPolyLine(1); - - - // Determine x- and y-dimensions depending on principal axis - // TODO use plane geometry normal to determine that automatically, then check whether the PF is aligned with one of the three principal axis - int i0, i1; - switch ( axis ) - { - case 0: - i0 = 1; - i1 = 2; - break; - - case 1: - i0 = 0; - i1 = 2; - break; - - case 2: - default: - i0 = 0; - i1 = 1; - break; - } - m_PlanarFigureCoordinate0= i0; - m_PlanarFigureCoordinate1= i1; - - // store the polyline contour as vtkPoints object - bool outOfBounds = false; - vtkSmartPointer points = vtkSmartPointer::New(); - typename PlanarFigure::PolyLineType::const_iterator it; - for ( it = planarFigurePolyline.begin(); - it != planarFigurePolyline.end(); - ++it ) - { - Point3D point3D; - - // Convert 2D point back to the local index coordinates of the selected - // image - // Fabian: From PlaneGeometry documentation: - // Converts a 2D point given in mm (pt2d_mm) relative to the upper-left corner of the geometry into the corresponding world-coordinate (a 3D point in mm, pt3d_mm). - // To convert a 2D point given in units (e.g., pixels in case of an image) into a 2D point given in mm (as required by this method), use IndexToWorld. - planarFigurePlaneGeometry->Map( *it, point3D ); - - // Polygons (partially) outside of the image bounds can not be processed - // further due to a bug in vtkPolyDataToImageStencil - if ( !imageGeometry3D->IsInside( point3D ) ) - { - outOfBounds = true; - } - - // Fabian: Why convert to index coordinates? - imageGeometry3D->WorldToIndex( point3D, point3D ); - - points->InsertNextPoint( point3D[i0], point3D[i1], 0 ); - } - - vtkSmartPointer holePoints = nullptr; - - if (!planarFigureHolePolyline.empty()) - { - holePoints = vtkSmartPointer::New(); - - Point3D point3D; - PlanarFigure::PolyLineType::const_iterator end = planarFigureHolePolyline.end(); - - for (it = planarFigureHolePolyline.begin(); it != end; ++it) - { - // Fabian: same as above - planarFigurePlaneGeometry->Map(*it, point3D); - imageGeometry3D->WorldToIndex(point3D, point3D); - holePoints->InsertNextPoint(point3D[i0], point3D[i1], 0); - } - } - - // mark a malformed 2D planar figure ( i.e. area = 0 ) as out of bounds - // this can happen when all control points of a rectangle lie on the same line = two of the three extents are zero - double bounds[6] = {0, 0, 0, 0, 0, 0}; - points->GetBounds( bounds ); - bool extent_x = (fabs(bounds[0] - bounds[1])) < mitk::eps; - bool extent_y = (fabs(bounds[2] - bounds[3])) < mitk::eps; - bool extent_z = (fabs(bounds[4] - bounds[5])) < mitk::eps; - - // throw an exception if a closed planar figure is deformed, i.e. has only one non-zero extent - if ( m_PlanarFigure->IsClosed() && - ((extent_x && extent_y) || (extent_x && extent_z) || (extent_y && extent_z))) - { - mitkThrow() << "Figure has a zero area and cannot be used for masking."; - } - - if ( outOfBounds ) - { - throw std::runtime_error( "Figure at least partially outside of image bounds!" ); - } - - // create a vtkLassoStencilSource and set the points of the Polygon - vtkSmartPointer lassoStencil = vtkSmartPointer::New(); - lassoStencil->SetShapeToPolygon(); - lassoStencil->SetPoints( points ); - - vtkSmartPointer holeLassoStencil = nullptr; - - if (holePoints.GetPointer() != nullptr) - { - holeLassoStencil = vtkSmartPointer::New(); - holeLassoStencil->SetShapeToPolygon(); - holeLassoStencil->SetPoints(holePoints); - } - - // Export from ITK to VTK (to use a VTK filter) - typedef itk::VTKImageImport< MaskImage2DType > ImageImportType; - typedef itk::VTKImageExport< MaskImage2DType > ImageExportType; - - typename ImageExportType::Pointer itkExporter = ImageExportType::New(); - itkExporter->SetInput( castFilter->GetOutput() ); - - vtkSmartPointer vtkImporter = vtkSmartPointer::New(); - this->ConnectPipelines( itkExporter, vtkImporter ); - - // Apply the generated image stencil to the input image - vtkSmartPointer imageStencilFilter = vtkSmartPointer::New(); - imageStencilFilter->SetInputConnection( vtkImporter->GetOutputPort() ); - imageStencilFilter->SetStencilConnection(lassoStencil->GetOutputPort()); - imageStencilFilter->ReverseStencilOff(); - imageStencilFilter->SetBackgroundValue( 0 ); - imageStencilFilter->Update(); - - vtkSmartPointer holeStencilFilter = nullptr; - - if (holeLassoStencil.GetPointer() != nullptr) - { - holeStencilFilter = vtkSmartPointer::New(); - holeStencilFilter->SetInputConnection(imageStencilFilter->GetOutputPort()); - holeStencilFilter->SetStencilConnection(holeLassoStencil->GetOutputPort()); - holeStencilFilter->ReverseStencilOn(); - holeStencilFilter->SetBackgroundValue(0); - holeStencilFilter->Update(); - } - - // Export from VTK back to ITK - vtkSmartPointer vtkExporter = vtkSmartPointer::New(); - vtkExporter->SetInputConnection( holeStencilFilter.GetPointer() == nullptr - ? imageStencilFilter->GetOutputPort() - : holeStencilFilter->GetOutputPort()); - vtkExporter->Update(); - - typename ImageImportType::Pointer itkImporter = ImageImportType::New(); - this->ConnectPipelines( vtkExporter, itkImporter ); - itkImporter->Update(); - - typedef itk::ImageDuplicator< ImageImportType::OutputImageType > DuplicatorType; - DuplicatorType::Pointer duplicator = DuplicatorType::New(); - duplicator->SetInputImage( itkImporter->GetOutput() ); - duplicator->Update(); - - // Store mask - m_InternalImageMask2D = duplicator->GetOutput(); - } - - - void ImageStatisticsCalculator::UnmaskedStatisticsProgressUpdate() - { - // Need to throw away every second progress event to reach a final count of - // 100 since two consecutive filters are used in this case - static int updateCounter = 0; - if ( updateCounter++ % 2 == 0 ) - { - this->InvokeEvent( itk::ProgressEvent() ); - } - } - - - void ImageStatisticsCalculator::MaskedStatisticsProgressUpdate() - { - this->InvokeEvent( itk::ProgressEvent() ); - } - -} diff --git a/Modules/ImageStatistics/old/mitkImageStatisticsCalculator.h b/Modules/ImageStatistics/old/mitkImageStatisticsCalculator.h deleted file mode 100644 index faaa2cf62c..0000000000 --- a/Modules/ImageStatistics/old/mitkImageStatisticsCalculator.h +++ /dev/null @@ -1,609 +0,0 @@ -/*=================================================================== -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 _MITK_IMAGESTATISTICSCALCULATOR_H -#define _MITK_IMAGESTATISTICSCALCULATOR_H - -#include -#include "MitkImageStatisticsExports.h" -#include -#include -#include "mitkImage.h" -#include "mitkPlanarFigure.h" - -#ifndef __itkHistogram_h -#include -#endif - -#include - - -#include - -// just a helper to unclutter our code -// to be replaced with references to m_Member (when deprecated public members in Statistics are removed) -#define mitkSetGetConstMacro(name, type) \ - virtual type Get##name() const \ - { \ - return this->name; \ - } \ - \ - virtual void Set##name(const type _arg) \ - { \ - if ( this->name != _arg ) \ - { \ - this->name = _arg; \ - } \ - } - -namespace mitk -{ - /** - * \brief Class for calculating statistics and histogram for an (optionally - * masked) image. - * - * Images can be masked by either a label image (of the same dimensions as - * the original image) or by a closed mitk::PlanarFigure, e.g. a circle or - * polygon. When masking with a planar figure, the slice corresponding to the - * plane containing the figure is extracted and then clipped with contour - * defined by the figure. Planar figures need to be aligned along the main axes - * of the image (axial, sagittal, coronal). Planar figures on arbitrary - * rotated planes are not supported. - * - * For each operating mode (no masking, masking by image, masking by planar - * figure), the calculated statistics and histogram are cached so that, when - * switching back and forth between operation modes without modifying mask or - * image, the information doesn't need to be recalculated. - * - * The class also has the possibility to calculate the location and separate - * statistics for a region called "hotspot". The hotspot is a sphere of - * user-defined size and its location is chosen in a way that the average - * pixel value within the sphere is maximized. - * - * \warning Hotspot calculation does not work in case of 2D-images! - * - * Note: currently time-resolved and multi-channel pictures are not properly - * supported. - * - * \section HotspotStatistics_caption Calculation of hotspot statistics - * - * Since calculation of hotspot location and statistics is not - * straight-forward, the following paragraphs will describe it in more detail. - * - * Note: Calculation of hotspot statistics is optional and set to off by default. - * Multilabel-masks are supported. - * - * \subsection HotspotStatistics_description Hotspot Definition - * - * The hotspot of an image is motivated from PET readings. It is defined - * as a spherical region of fixed size which maximizes the average pixel value - * within the region. The following image illustrates the concept: the - * colored areas are different image intensities and the hotspot is located - * in the hottest region of the image. - * - * Note: Only hotspots are calculated for which the whole hotspot-sphere is - * inside the image by default. This behaviour can be changed by - * by calling SetHotspotMustBeCompletlyInsideImage(). - * \warning Note that SetHotspotMustBeCompletlyInsideImage(false) may overrate - * "hot" regions at image borders, because they have a stronger influence on the - * mean value! Think clearly about this fact and make sure this is what you - * want/need in your application, before calling - * SetHotspotMustBeCompletlyInsideImage(false)! - * - * - * \image html hotspotexample.JPG - * - * \subsection HotspotStatistics_calculation Hotspot Calculation - * - * Since only the size of the hotspot is known initially, we need to calculate - * two aspects (both implemented in CalculateHotspotStatistics() ): - * - the hotspot location - * - statistics of the pixels within the hotspot. - * - * Finding the hotspot location requires to calculate the average value at each - * position. This is done by convolution of the image with a sperical kernel - * image which reflects partial volumes (important in the case of low-resolution - * PET images). - * - * Once the hotspot location is known, calculating the actual statistics is a - * simple task which is implemented in CalculateHotspotStatistics() using a second - * instance of the ImageStatisticsCalculator. - * - * Step 1: Finding the hotspot by image convolution - * - * As described above, we use image convolution with a rasterized sphere to - * average the image at each position. To handle coarse resolutions, which would - * normally force us to decide for partially contained voxels whether to count - * them or not, we supersample the kernel image and use non-integer kernel values - * (see GenerateHotspotSearchConvolutionKernel()), which reflect the volume part that is contained in the - * sphere. For example, if three subvoxels are inside the sphere, the corresponding - * kernel voxel gets a value of 0.75 (3 out of 4 subvoxels, see 2D example below). - * - * \image html convolutionkernelsupersampling.jpg - * - * Convolution itself is done by means of the itkFFTConvolutionImageFilter. - * To find the hotspot location, we simply iterate the averaged image and find a - * maximum location (see CalculateExtremaWorld()). In case of images with multiple - * maxima the method returns value and corresponding index of the extrema that is - * found by the iterator first. - * - * Step 2: Computation of hotspot statistics - * - * Once the hotspot location is found, statistics for the region are calculated - * by simply iterating the input image and regarding all pixel centers inside the - * hotspot-sphere for statistics. - * \warning Index positions of maximum/minimum are not provided, because they are not necessarily unique - * \todo If index positions of maximum/minimum are required, output needs to be changed to multiple positions / regions, etc. - * - * \subsection HotspotStatistics_tests Tests - * - * To check the correctness of the hotspot calculation, a special class - * (\ref hotspottestdoc) has been created, which generates images with - * known hotspot location and statistics. A number of unit tests use this class - * to first generate an image of known properites and then verify that - * ImageStatisticsCalculator is able to reproduce the known statistics. - * - */ - class MITKIMAGESTATISTICS_EXPORT ImageStatisticsCalculator : public itk::Object - { - public: - - /** \brief Enum for possible masking modi. */ - enum - { - MASKING_MODE_NONE = 0, - MASKING_MODE_IMAGE = 1, - MASKING_MODE_PLANARFIGURE = 2 - }; - - typedef itk::Statistics::Histogram HistogramType; - typedef HistogramType::ConstIterator HistogramConstIteratorType; - - /** \brief Class for common statistics, includig hotspot properties. */ - class MITKIMAGESTATISTICS_EXPORT Statistics - { - public: - - Statistics(bool withHotspotStatistics = true); - Statistics(const Statistics& other); - - virtual ~Statistics(); - - Statistics& operator=(Statistics const& stats); - - const Statistics& GetHotspotStatistics() const; // real statistics - Statistics& GetHotspotStatistics(); // real statistics - bool HasHotspotStatistics() const; - void SetHasHotspotStatistics(bool hasHotspotStatistics); // set a flag. if set, return empty hotspotstatistics object - - void Reset(unsigned int dimension = 2); - - mitkSetGetConstMacro(Label, unsigned int) - mitkSetGetConstMacro(N, unsigned int) - mitkSetGetConstMacro(Min, double) - mitkSetGetConstMacro(Max, double) - mitkSetGetConstMacro(Mean, double) - mitkSetGetConstMacro(Median, double) - - double GetVariance() const; - /** \brief Set variance - * - * This method checks whether the variance is negative: - * The reason that the variance may be negative is that the underlying itk::LabelStatisticsImageFilter uses a naïve algorithm - * for calculating the variance ( http://en.wikipedia.org/wiki/Algorithms_for_calculating_variance ) which can lead to negative values - * due to rounding errors. - * - * If the variance is negative the value will be set to 0.0, else the given value will be set. - */ - void SetVariance( const double ); - - double GetSigma() const; - /** \brief Set standard deviation (sigma) - * - * This method checks if the given standard deviation is a positive value. This is done because the underlying itk::LabelStatisticsImageFilter uses - * a naïve algorithm to calculate the variance. This may lead to a negative variance and because the square root of the variance is taken it also - * leads to NaN for sigma. - * - * If the given value is not reasonable the value will be set to 0.0, else the given value will be set. - * - * \see SetVariance() - */ - void SetSigma( const double ); - - mitkSetGetConstMacro(RMS, double) - mitkSetGetConstMacro(MinIndex, vnl_vector) - mitkSetGetConstMacro(MaxIndex, vnl_vector) - mitkSetGetConstMacro(HotspotIndex, vnl_vector) - - mitkSetGetConstMacro(Skewness, double) - mitkSetGetConstMacro(Kurtosis, double) - mitkSetGetConstMacro(Uniformity, double) - mitkSetGetConstMacro(Entropy, double) - mitkSetGetConstMacro(UPP, double); - mitkSetGetConstMacro(MPP, double); - - private: - unsigned int Label; - unsigned int N; - double Min; - double Max; - double Mean; - double Median; - double Variance; - double Sigma; - double RMS; - vnl_vector MinIndex; - vnl_vector MaxIndex; - - bool m_SkewnessAndKurtosisSuccesfulCalculated; - bool m_ImageStatisticsSuccesfulCalculated; - - - Statistics* m_HotspotStatistics; - - double Skewness; - double Kurtosis; - double Uniformity; - double Entropy; - double UPP; - double MPP; - - bool m_HasHotspotStatistics; - vnl_vector HotspotIndex; //< index of hotspotsphere origin - }; - - typedef std::vector< HistogramType::ConstPointer > HistogramContainer; - typedef std::vector< Statistics > StatisticsContainer; - - mitkClassMacroItkParent( ImageStatisticsCalculator, itk::Object ); - itkFactorylessNewMacro(Self) - itkCloneMacro(Self) - - /** \brief Automatically calculate bin size to obtain 200 bins. */ - void SetUseDefaultBinSize(bool useDefault); - - /** \brief Set image from which to compute statistics. */ - void SetImage( const mitk::Image *image ); - - /** \brief Set image for masking. */ - void SetImageMask( const mitk::Image *imageMask ); - - /** \brief Set planar figure for masking. */ - void SetPlanarFigure( mitk::PlanarFigure *planarFigure ); - - - /** \brief Set/Get operation mode for masking */ - void SetMaskingMode( unsigned int mode ); - - /** \brief Set/Get operation mode for masking */ - itkGetMacro( MaskingMode, unsigned int ); - - /** \brief Set/Get operation mode for masking */ - void SetMaskingModeToNone(); - - /** \brief Set/Get operation mode for masking */ - void SetMaskingModeToImage(); - - /** \brief Set/Get operation mode for masking */ - void SetMaskingModeToPlanarFigure(); - - /** \brief Set a pixel value for pixels that will be ignored in the statistics */ - void SetIgnorePixelValue(double value); - - /** \brief Get the pixel value for pixels that will be ignored in the statistics */ - double GetIgnorePixelValue(); - - /** \brief Set whether a pixel value should be ignored in the statistics */ - void SetDoIgnorePixelValue(bool doit); - - /** \brief Get whether a pixel value will be ignored in the statistics */ - bool GetDoIgnorePixelValue(); - - /** \brief Set bin size for histogram resolution.*/ - void SetHistogramBinSize( double size); - - /** \brief Get bin size for histogram resolution.*/ - double GetHistogramBinSize(); - - /** \brief Sets the radius for the hotspot */ - void SetHotspotRadiusInMM (double hotspotRadiusInMM); - - /** \brief Returns the radius of the hotspot */ - double GetHotspotRadiusInMM(); - - /** \brief Sets whether the hotspot should be calculated */ - void SetCalculateHotspot(bool calculateHotspot); - - /** \brief Returns true whether the hotspot should be calculated, otherwise false */ - bool IsHotspotCalculated(); - - /** \brief Sets flag whether hotspot is completly inside the image. Please note that if set to false - it can be possible that statistics are calculated for which the whole hotspot is not inside the image! - - \warning regarding positions at the image centers may produce unexpected hotspot locations, please see \ref HotspotStatistics_description - */ - void SetHotspotMustBeCompletlyInsideImage(bool hotspotIsCompletlyInsideImage, bool warn = true); - - /** \brief Returns true if hotspot has to be completly inside the image. */ - bool GetHotspotMustBeCompletlyInsideImage() const; - - /** \brief Compute statistics (together with histogram) for the current - * masking mode. - * - * Computation is not executed if statistics is already up to date. In this - * case, false is returned; otherwise, true.*/ - virtual bool ComputeStatistics( unsigned int timeStep = 0 ); - - - /** \brief Retrieve the histogram depending on the current masking mode. - * - * \param label The label for which to retrieve the histogram in multi-label situations (ascending order). - */ - const HistogramType *GetHistogram( unsigned int timeStep = 0, unsigned int label = 0 ) const; - - /** \brief Retrieve the histogram depending on the current masking mode (for all image labels. */ - const HistogramContainer &GetHistogramVector( unsigned int timeStep = 0 ) const; - - typedef std::map BinFrequencyType; - - /* Returning a map including bin and Frequency*/ - BinFrequencyType GetBinsAndFreuqencyForHistograms( unsigned int timeStep = 0, unsigned int label = 0) const; - - /** \brief Retrieve statistics depending on the current masking mode. - * - - * \param label The label for which to retrieve the statistics in multi-label situations (ascending order). - */ - const Statistics &GetStatistics( unsigned int timeStep = 0, unsigned int label = 0 ) const; - - /** \brief Retrieve statistics depending on the current masking mode (for all image labels). */ - const StatisticsContainer &GetStatisticsVector( unsigned int timeStep = 0 ) const; - - - - - protected: - typedef std::vector< HistogramContainer > HistogramVector; - typedef std::vector< StatisticsContainer > StatisticsVector; - - typedef std::vector< itk::TimeStamp > TimeStampVectorType; - typedef std::vector< bool > BoolVectorType; - - typedef itk::Image< unsigned short, 3 > MaskImage3DType; - typedef itk::Image< unsigned short, 2 > MaskImage2DType; - - ImageStatisticsCalculator(); - - virtual ~ImageStatisticsCalculator(); - - /** \brief Depending on the masking mode, the image and mask from which to - * calculate statistics is extracted from the original input image and mask - * data. - * - * For example, a when using a PlanarFigure as mask, the 2D image slice - * corresponding to the PlanarFigure will be extracted from the original - * image. If masking is disabled, the original image is simply passed - * through. */ - void ExtractImageAndMask( unsigned int timeStep = 0 ); - - /*calculate the min and max value, this is done because we need the min and max value before execution the statistics filter to have the wright range for the histogramm*/ - template < typename TPixel, unsigned int VImageDimension > - void GetMinAndMaxValue(double &minimum, double &maximum, int &counter, double &sigma, const itk::Image< TPixel, VImageDimension > *InputImage, - itk::Image< unsigned short, VImageDimension > *MaskImageType); - - /** \brief If the passed vector matches any of the three principal axes - * of the passed geometry, the ínteger value corresponding to the axis - * is set and true is returned. */ - bool GetPrincipalAxis( const BaseGeometry *geometry, Vector3D vector, - unsigned int &axis ); - - template < typename TPixel, unsigned int VImageDimension > - void InternalCalculateStatisticsUnmasked( - const itk::Image< TPixel, VImageDimension > *image, - StatisticsContainer* statisticsContainer, - HistogramContainer *histogramContainer ); - - template < typename TPixel, unsigned int VImageDimension > - void InternalCalculateStatisticsMasked( - const itk::Image< TPixel, VImageDimension > *image, - itk::Image< unsigned short, VImageDimension > *maskImage, - StatisticsContainer* statisticsContainer, - HistogramContainer* histogramContainer ); - - template < typename TPixel, unsigned int VImageDimension > - void InternalCalculateMaskFromPlanarFigure( - const itk::Image< TPixel, VImageDimension > *image, unsigned int axis ); - - template < typename TPixel, unsigned int VImageDimension > - void InternalMaskIgnoredPixels( - const itk::Image< TPixel, VImageDimension > *image, - itk::Image< unsigned short, VImageDimension > *maskImage ); - - class ImageExtrema - { - public: - bool Defined; - double Max; - double Min; - vnl_vector MaxIndex; - vnl_vector MinIndex; - - ImageExtrema() - :Defined(false) - ,Max(itk::NumericTraits::min()) - ,Min(itk::NumericTraits::max()) - { - } - }; - - - /** \brief Calculates minimum, maximum, mean value and their - * corresponding indices in a given ROI. As input the function - * needs an image and a mask. Returns an ImageExtrema object. */ - template - ImageExtrema CalculateExtremaWorld( - const itk::Image *inputImage, - itk::Image *maskImage, - double neccessaryDistanceToImageBorderInMM, - unsigned int label); - - - /** \brief Calculates the hotspot statistics depending on - * masking mode. Hotspot statistics are calculated for a - * hotspot which is completly located inside the image by default. */ - template < typename TPixel, unsigned int VImageDimension> - Statistics CalculateHotspotStatistics( - const itk::Image *inputImage, - itk::Image *maskImage, - double radiusInMM, - bool& isHotspotDefined, - unsigned int label); - - /** Connection from ITK to VTK */ - template - void ConnectPipelines(ITK_Exporter exporter, vtkSmartPointer importer) - { - importer->SetUpdateInformationCallback(exporter->GetUpdateInformationCallback()); - - importer->SetPipelineModifiedCallback(exporter->GetPipelineModifiedCallback()); - importer->SetWholeExtentCallback(exporter->GetWholeExtentCallback()); - importer->SetSpacingCallback(exporter->GetSpacingCallback()); - importer->SetOriginCallback(exporter->GetOriginCallback()); - importer->SetScalarTypeCallback(exporter->GetScalarTypeCallback()); - - importer->SetNumberOfComponentsCallback(exporter->GetNumberOfComponentsCallback()); - - importer->SetPropagateUpdateExtentCallback(exporter->GetPropagateUpdateExtentCallback()); - importer->SetUpdateDataCallback(exporter->GetUpdateDataCallback()); - importer->SetDataExtentCallback(exporter->GetDataExtentCallback()); - importer->SetBufferPointerCallback(exporter->GetBufferPointerCallback()); - importer->SetCallbackUserData(exporter->GetCallbackUserData()); - } - - /** Connection from VTK to ITK */ - template - void ConnectPipelines(vtkSmartPointer exporter, ITK_Importer importer) - { - importer->SetUpdateInformationCallback(exporter->GetUpdateInformationCallback()); - - importer->SetPipelineModifiedCallback(exporter->GetPipelineModifiedCallback()); - importer->SetWholeExtentCallback(exporter->GetWholeExtentCallback()); - importer->SetSpacingCallback(exporter->GetSpacingCallback()); - importer->SetOriginCallback(exporter->GetOriginCallback()); - importer->SetScalarTypeCallback(exporter->GetScalarTypeCallback()); - - importer->SetNumberOfComponentsCallback(exporter->GetNumberOfComponentsCallback()); - - importer->SetPropagateUpdateExtentCallback(exporter->GetPropagateUpdateExtentCallback()); - importer->SetUpdateDataCallback(exporter->GetUpdateDataCallback()); - importer->SetDataExtentCallback(exporter->GetDataExtentCallback()); - importer->SetBufferPointerCallback(exporter->GetBufferPointerCallback()); - importer->SetCallbackUserData(exporter->GetCallbackUserData()); - } - - - void UnmaskedStatisticsProgressUpdate(); - - void MaskedStatisticsProgressUpdate(); - - /** \brief Returns size of convolution kernel depending on spacing and radius. */ - template - itk::Size - CalculateConvolutionKernelSize(double spacing[VImageDimension], double radiusInMM); - - /** \brief Generates image of kernel which is needed for convolution. */ - template - itk::SmartPointer< itk::Image > - GenerateHotspotSearchConvolutionKernel(double spacing[VImageDimension], double radiusInMM); - - /** \brief Convolves image with spherical kernel image. Used for hotspot calculation. */ - template - itk::SmartPointer< itk::Image > - GenerateConvolutionImage( const itk::Image* inputImage ); - - /** \brief Fills pixels of the spherical hotspot mask. */ - template < typename TPixel, unsigned int VImageDimension> - void - FillHotspotMaskPixels( itk::Image* maskImage, - itk::Point sphereCenter, - double sphereRadiusInMM); - - /** m_Image contains the input image (e.g. 2D, 3D, 3D+t)*/ - mitk::Image::ConstPointer m_Image; - - mitk::Image::ConstPointer m_ImageMask; - - mitk::PlanarFigure::Pointer m_PlanarFigure; - - HistogramVector m_ImageHistogramVector; - HistogramVector m_MaskedImageHistogramVector; - HistogramVector m_PlanarFigureHistogramVector; - - HistogramType::Pointer m_EmptyHistogram; - HistogramContainer m_EmptyHistogramContainer; - - - StatisticsVector m_ImageStatisticsVector; - StatisticsVector m_MaskedImageStatisticsVector; - StatisticsVector m_PlanarFigureStatisticsVector; - StatisticsVector m_MaskedImageHotspotStatisticsVector; - - Statistics m_EmptyStatistics; - StatisticsContainer m_EmptyStatisticsContainer; - - unsigned int m_MaskingMode; - bool m_MaskingModeChanged; - - /** m_InternalImage contains a image volume at one time step (e.g. 2D, 3D)*/ - mitk::Image::ConstPointer m_InternalImage; - MaskImage3DType::Pointer m_InternalImageMask3D; - MaskImage2DType::Pointer m_InternalImageMask2D; - - TimeStampVectorType m_ImageStatisticsTimeStampVector; - TimeStampVectorType m_MaskedImageStatisticsTimeStampVector; - TimeStampVectorType m_PlanarFigureStatisticsTimeStampVector; - - BoolVectorType m_ImageStatisticsCalculationTriggerVector; - BoolVectorType m_MaskedImageStatisticsCalculationTriggerVector; - BoolVectorType m_PlanarFigureStatisticsCalculationTriggerVector; - - double m_IgnorePixelValue; - bool m_DoIgnorePixelValue; - bool m_IgnorePixelValueChanged; - - unsigned int m_PlanarFigureAxis; // Normal axis for PlanarFigure - unsigned int m_PlanarFigureSlice; // Slice which contains PlanarFigure - int m_PlanarFigureCoordinate0; // First plane-axis for PlanarFigure - int m_PlanarFigureCoordinate1; // Second plane-axis for PlanarFigure - - double m_HistogramBinSize; ///Bin size for histogram resoluion. - bool m_UseDefaultBinSize; - bool m_UseBinSizeBasedOnVOIRegion; - double m_HotspotRadiusInMM; - bool m_CalculateHotspot; - bool m_HotspotRadiusInMMChanged; - bool m_HotspotMustBeCompletelyInsideImage; - - - private: - - unsigned int calcNumberOfBins(mitk::ScalarType min, mitk::ScalarType max); - - - }; - -} // namespace - -#endif