diff --git a/Modules/ImageStatistics/mitkExtendedLabelStatisticsImageFilter.h b/Modules/ImageStatistics/mitkExtendedLabelStatisticsImageFilter.h index 785b63047d..54d3d46f47 100644 --- a/Modules/ImageStatistics/mitkExtendedLabelStatisticsImageFilter.h +++ b/Modules/ImageStatistics/mitkExtendedLabelStatisticsImageFilter.h @@ -1,353 +1,348 @@ /*=================================================================== 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 #define __mitkExtendedLabelStatisticsImageFilter #include "itkLabelStatisticsImageFilter.h" namespace itk { /** * \class ExtendedLabelStatisticsImageFilter * \brief Extension of the itkLabelStatisticsImageFilter that also calculates the Skewness,Kurtosis,Entropy,Uniformity. * * This class inherits from the itkLabelStatisticsImageFilter and - * uses its results for the calculation of six additional coefficients: - * the Skewness,Kurtosis,Uniformity,UPP,MPP,Entropy + * uses its results for the calculation of seven additional coefficients: + * the Skewness, Kurtosis, Uniformity, UPP, MPP, Entropy and Median * - * 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); /** \class LabelStatistics * \brief Statistics stored per label * \ingroup ITKImageStatistics */ class LabelStatistics { public: // default constructor LabelStatistics() { // initialized to the default values m_Count = NumericTraits< IdentifierType >::ZeroValue(); m_PositivePixelCount = NumericTraits< IdentifierType >::ZeroValue(); m_Sum = NumericTraits< RealType >::ZeroValue(); m_SumOfPositivePixels = NumericTraits< RealType >::ZeroValue(); m_SumOfSquares = NumericTraits< RealType >::ZeroValue(); m_SumOfCubes = NumericTraits< RealType >::ZeroValue(); m_SumOfQuadruples = NumericTraits< RealType >::ZeroValue(); // Set such that the first pixel encountered can be compared m_Minimum = NumericTraits< RealType >::max(); m_Maximum = NumericTraits< RealType >::NonpositiveMin(); // Default these to zero m_Mean = NumericTraits< RealType >::ZeroValue(); m_Sigma = NumericTraits< RealType >::ZeroValue(); m_Variance = NumericTraits< RealType >::ZeroValue(); m_MPP = NumericTraits< RealType >::ZeroValue(); m_Median = NumericTraits< RealType >::ZeroValue(); m_Uniformity = NumericTraits< RealType >::ZeroValue(); m_UPP = NumericTraits< RealType >::ZeroValue(); m_Entropy = NumericTraits< RealType >::ZeroValue(); m_Skewness = NumericTraits< RealType >::ZeroValue(); m_Kurtosis = NumericTraits< RealType >::ZeroValue(); unsigned int imageDimension = itkGetStaticConstMacro(ImageDimension); m_BoundingBox.resize(imageDimension * 2); for ( unsigned int i = 0; i < imageDimension * 2; i += 2 ) { m_BoundingBox[i] = NumericTraits< IndexValueType >::max(); m_BoundingBox[i + 1] = NumericTraits< IndexValueType >::NonpositiveMin(); } m_Histogram = ITK_NULLPTR; } // constructor with histogram enabled LabelStatistics(int size, RealType lowerBound, RealType upperBound) { // initialized to the default values m_Count = NumericTraits< IdentifierType >::ZeroValue(); m_PositivePixelCount = NumericTraits< IdentifierType >::ZeroValue(); m_Sum = NumericTraits< RealType >::ZeroValue(); m_SumOfPositivePixels = NumericTraits< RealType >::ZeroValue(); m_SumOfSquares = NumericTraits< RealType >::ZeroValue(); m_SumOfCubes = NumericTraits< RealType >::ZeroValue(); m_SumOfQuadruples = NumericTraits< RealType >::ZeroValue(); // Set such that the first pixel encountered can be compared m_Minimum = NumericTraits< RealType >::max(); m_Maximum = NumericTraits< RealType >::NonpositiveMin(); // Default these to zero m_Mean = NumericTraits< RealType >::ZeroValue(); m_Sigma = NumericTraits< RealType >::ZeroValue(); m_Variance = NumericTraits< RealType >::ZeroValue(); m_MPP = NumericTraits< RealType >::ZeroValue(); m_Median = NumericTraits< RealType >::ZeroValue(); m_Uniformity = NumericTraits< RealType >::ZeroValue(); m_UPP = NumericTraits< RealType >::ZeroValue(); m_Entropy = NumericTraits< RealType >::ZeroValue(); m_Skewness = NumericTraits< RealType >::ZeroValue(); m_Kurtosis = NumericTraits< RealType >::ZeroValue(); unsigned int imageDimension = itkGetStaticConstMacro(ImageDimension); m_BoundingBox.resize(imageDimension * 2); for ( unsigned int i = 0; i < imageDimension * 2; i += 2 ) { m_BoundingBox[i] = NumericTraits< IndexValueType >::max(); m_BoundingBox[i + 1] = NumericTraits< IndexValueType >::NonpositiveMin(); } // Histogram m_Histogram = HistogramType::New(); typename HistogramType::SizeType hsize; typename HistogramType::MeasurementVectorType lb; typename HistogramType::MeasurementVectorType ub; hsize.SetSize(1); lb.SetSize(1); ub.SetSize(1); m_Histogram->SetMeasurementVectorSize(1); hsize[0] = size; lb[0] = lowerBound; ub[0] = upperBound; m_Histogram->Initialize(hsize, lb, ub); } // need copy constructor because of smart pointer to histogram LabelStatistics(const LabelStatistics & l) { m_Count = l.m_Count; m_Minimum = l.m_Minimum; m_Maximum = l.m_Maximum; m_Mean = l.m_Mean; m_Sum = l.m_Sum; m_SumOfSquares = l.m_SumOfSquares; m_Sigma = l.m_Sigma; m_Variance = l.m_Variance; m_MPP = l.m_MPP; m_Median = l.m_Median; m_Uniformity = l.m_Uniformity; m_UPP = l.m_UPP; m_Entropy = l.m_Entropy; m_Skewness = l.m_Skewness; m_Kurtosis = l.m_Kurtosis; m_BoundingBox = l.m_BoundingBox; m_Histogram = l.m_Histogram; m_SumOfPositivePixels = l.m_SumOfPositivePixels; m_PositivePixelCount = l.m_PositivePixelCount; m_SumOfCubes = l.m_SumOfCubes; m_SumOfQuadruples = l.m_SumOfQuadruples; } // added for completeness LabelStatistics &operator= (const LabelStatistics& l) { if(this != &l) { m_Count = l.m_Count; m_Minimum = l.m_Minimum; m_Maximum = l.m_Maximum; m_Mean = l.m_Mean; m_Sum = l.m_Sum; m_SumOfSquares = l.m_SumOfSquares; m_Sigma = l.m_Sigma; m_Variance = l.m_Variance; m_MPP = l.m_MPP; m_Median = l.m_Median; m_Uniformity = l.m_Uniformity; m_UPP = l.m_UPP; m_Entropy = l.m_Entropy; m_Skewness = l.m_Skewness; m_Kurtosis = l.m_Kurtosis; m_BoundingBox = l.m_BoundingBox; m_Histogram = l.m_Histogram; m_SumOfPositivePixels = l.m_SumOfPositivePixels; m_PositivePixelCount = l.m_PositivePixelCount; m_SumOfCubes = l.m_SumOfCubes; m_SumOfQuadruples = l.m_SumOfQuadruples; } return *this; } IdentifierType m_Count; RealType m_Minimum; RealType m_Maximum; RealType m_Mean; RealType m_Sum; RealType m_SumOfSquares; RealType m_Sigma; RealType m_Variance; RealType m_MPP; RealType m_Median; RealType m_Uniformity; RealType m_UPP; RealType m_Entropy; RealType m_Skewness; RealType m_Kurtosis; IdentifierType m_PositivePixelCount; RealType m_SumOfPositivePixels; RealType m_SumOfCubes; RealType m_SumOfQuadruples; typename Superclass::BoundingBoxType m_BoundingBox; typename HistogramType::Pointer m_Histogram; }; /** Type of the map used to store data per label */ typedef itksys::hash_map< LabelPixelType, LabelStatistics > MapType; typedef typename itksys::hash_map< LabelPixelType, LabelStatistics >::const_iterator StatisticsMapConstIterator; typedef typename itksys::hash_map< LabelPixelType, LabelStatistics >::iterator StatisticsMapIterator; typedef IdentifierType MapSizeType; /** Type of the container used to store valid label values */ typedef std::vector ValidLabelValuesContainerType; /** Return the computed Minimum for a label. */ RealType GetMinimum(LabelPixelType label) const; /** Return the computed Maximum for a label. */ RealType GetMaximum(LabelPixelType label) const; /** Return the computed Mean for a label. */ RealType GetMean(LabelPixelType label) const; /** Return the computed Standard Deviation for a label. */ RealType GetSigma(LabelPixelType label) const; /** Return the computed Variance for a label. */ RealType GetVariance(LabelPixelType label) const; /** Return the computed bounding box for a label. */ typename Superclass::BoundingBoxType GetBoundingBox(LabelPixelType label) const; /** Return the computed region. */ typename Superclass::RegionType GetRegion(LabelPixelType label) const; /** Return the compute Sum for a label. */ RealType GetSum(LabelPixelType label) const; /** Return the number of pixels for a label. */ MapSizeType GetCount(LabelPixelType label) const; /** Return the histogram for a label */ HistogramType::Pointer GetHistogram(LabelPixelType label) const; /*getter method for the new statistics*/ RealType GetSkewness(LabelPixelType label) const; RealType GetKurtosis(LabelPixelType label) const; RealType GetUniformity( LabelPixelType label) const; RealType GetMedian( LabelPixelType label) const; RealType GetEntropy( LabelPixelType label) const; RealType GetMPP( LabelPixelType label) const; RealType GetUPP( LabelPixelType label) const; bool GetMaskingNonEmpty() const; std::list GetRelevantLabels() const; /** specify global Histogram parameters. If the histogram parameters are set with this function, the same min and max value are used for all histograms. */ void SetHistogramParameters(const int numBins, RealType lowerBound, RealType upperBound); /** specify Histogram parameters for each label individually. Labels in the label image that are not represented in the std::maps here will receive global parameters (if available) */ void SetHistogramParametersForLabels(std::map numBins, std::map lowerBound, std::map upperBound); protected: ExtendedLabelStatisticsImageFilter(): m_GlobalHistogramParametersSet(false), m_LabelHistogramParametersSet(false), m_PreferGlobalHistogramParameters(false) { m_NumBins.set_size(1); } virtual ~ExtendedLabelStatisticsImageFilter(){} void AfterThreadedGenerateData(); /** Initialize some accumulators before the threads run. */ void BeforeThreadedGenerateData(); /** Multi-thread version GenerateData. */ void ThreadedGenerateData(const typename TInputImage::RegionType & outputRegionForThread, ThreadIdType threadId); /** Does the specified label exist? Can only be called after a call * a call to Update(). */ bool HasLabel(LabelPixelType label) const { return m_LabelStatistics.find(label) != m_LabelStatistics.end(); } private: std::vector< MapType > m_LabelStatisticsPerThread; MapType m_LabelStatistics; ValidLabelValuesContainerType m_ValidLabelValues; bool m_GlobalHistogramParametersSet; typename HistogramType::SizeType m_NumBins; RealType m_LowerBound; RealType m_UpperBound; bool m_MaskNonEmpty; bool m_LabelHistogramParametersSet; std::map m_LabelMin, m_LabelMax; std::map m_LabelNBins; bool m_PreferGlobalHistogramParameters; }; // end of class } // end namespace itk #ifndef ITK_MANUAL_INSTANTIATION #include "mitkExtendedLabelStatisticsImageFilter.hxx" #endif #endif diff --git a/Modules/ImageStatistics/mitkExtendedStatisticsImageFilter.h b/Modules/ImageStatistics/mitkExtendedStatisticsImageFilter.h index f7518ba18f..912bec98dd 100644 --- a/Modules/ImageStatistics/mitkExtendedStatisticsImageFilter.h +++ b/Modules/ImageStatistics/mitkExtendedStatisticsImageFilter.h @@ -1,216 +1,210 @@ /*=================================================================== 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 #define __mitkExtendedStatisticsImageFilter #include "itkStatisticsImageFilter.h" #include #include namespace itk { /** * \class ExtendedStatisticsImageFilter * \brief Extension of the itkStatisticsImageFilter that also calculates the Skewness and Kurtosis. * * This class inherits from the itkStatisticsImageFilter and * uses its results for the calculation of other additional coefficients: * the Skewness and Kurtosis. * - * As 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. + * StatisticsImageFilter, the skewness, kurtosis, MPP, UPP, Uniformity, Entropy and Median are stored in the outputs + * 7 to 14 by this filter. */ template< class TInputImage > class ExtendedStatisticsImageFilter : public StatisticsImageFilter< TInputImage > { public: /** Standard Self typedef */ typedef ExtendedStatisticsImageFilter Self; typedef StatisticsImageFilter< TInputImage > Superclass; typedef SmartPointer< Self > Pointer; typedef SmartPointer< const Self > ConstPointer; typedef typename Superclass::RealType RealType; typedef typename Superclass::RealObjectType RealObjectType; typedef typename Superclass::PixelType PixelType; /** Histogram-related typedefs */ typedef itk::Statistics::Histogram< RealType > HistogramType; typedef typename HistogramType::Pointer HistogramPointer; itkFactorylessNewMacro( Self ); itkCloneMacro( Self ); itkTypeMacro( ExtendedStatisticsImageFilter, StatisticsImageFilter ); /** * \brief Return the computed Skewness. */ double GetSkewness() const { return this->GetSkewnessOutput()->Get(); } /** * \brief Return the computed Median */ double GetMedian() const { return this->GetMedianOutput()->Get(); } /** * \brief Return the computed Kurtosis. */ double GetKurtosis() const { return this->GetKurtosisOutput()->Get(); } /* \brief Return the computed MPP. */ double GetMPP() const { return this->GetMPPOutput()->Get(); } /** * \brief Return the computed Uniformity. */ double GetUniformity() const { return this->GetUniformityOutput()->Get(); } /** *\brief Return the computed Entropy. */ double GetEntropy() const { return this->GetEntropyOutput()->Get(); } /** * \brief Return the computed UPP. */ double GetUPP() const { return this->GetUPPOutput()->Get(); } /** * \brief Return the computed Histogram. */ const typename HistogramType::Pointer GetHistogram() { if (m_HistogramCalculated) { return m_Histogram; } else { return ITK_NULLPTR; } } /** specify Histogram parameters */ void SetHistogramParameters(const int numBins, RealType lowerBound, RealType upperBound); protected: ExtendedStatisticsImageFilter(); virtual ~ExtendedStatisticsImageFilter(){}; void BeforeThreadedGenerateData(); /** Multi-thread version GenerateData. */ void ThreadedGenerateData(const typename StatisticsImageFilter::RegionType & outputRegionForThread, ThreadIdType threadId); /** * brief Calls AfterThreadedGenerateData() of the superclass and the main methods */ void AfterThreadedGenerateData(); RealObjectType* GetSkewnessOutput(); const RealObjectType* GetSkewnessOutput() const; RealObjectType* GetKurtosisOutput(); const RealObjectType* GetKurtosisOutput() const; RealObjectType* GetMPPOutput(); const RealObjectType* GetMPPOutput() const; RealObjectType* GetEntropyOutput(); const RealObjectType* GetEntropyOutput() const; RealObjectType* GetUniformityOutput(); const RealObjectType* GetUniformityOutput() const; RealObjectType* GetUPPOutput(); const RealObjectType* GetUPPOutput() const; RealObjectType* GetMedianOutput(); const RealObjectType* GetMedianOutput() const; virtual DataObject::Pointer MakeOutput( ProcessObject::DataObjectPointerArraySizeType idx ); private: Array< RealType > m_ThreadSum; Array< RealType > m_SumOfSquares; Array< RealType > m_SumOfCubes; Array< RealType > m_SumOfQuadruples; Array< SizeValueType > m_Count; Array< SizeValueType > m_PositivePixelCount; Array< RealType > m_ThreadSumOfPositivePixels; Array< PixelType > m_ThreadMin; Array< PixelType > m_ThreadMax; std::vector< HistogramPointer > m_HistogramPerThread; HistogramPointer m_Histogram; bool m_UseHistogram; bool m_HistogramCalculated; RealType m_LowerBound, m_UpperBound; int m_NumBins; }; // end of class } // end namespace itk #ifndef ITK_MANUAL_INSTANTIATION #include "mitkExtendedStatisticsImageFilter.hxx" #endif #endif diff --git a/Modules/ImageStatistics/mitkHistogramStatisticsCalculator.h b/Modules/ImageStatistics/mitkHistogramStatisticsCalculator.h index 84eaadec55..a7e3e3d450 100644 --- a/Modules/ImageStatistics/mitkHistogramStatisticsCalculator.h +++ b/Modules/ImageStatistics/mitkHistogramStatisticsCalculator.h @@ -1,40 +1,49 @@ #ifndef MITKHISTOGRAMSTATISTICSCALCULATOR #define MITKHISTOGRAMSTATISTICSCALCULATOR #include #include #include namespace mitk { +/** + * @brief Computes basic histogram statistics such as Uniformity, UPP (Uniformity of positive entries), Entropy and Median (approximation) + */ class MITKIMAGESTATISTICS_EXPORT HistogramStatisticsCalculator { public: typedef double MeasurementType; typedef itk::Statistics::Histogram HistogramType; HistogramStatisticsCalculator(); + /** + * @brief SetHistogram requires a itk::Statistics::Histogram + */ void SetHistogram(HistogramType::Pointer histogram); MeasurementType GetUPP(); MeasurementType GetUniformity(); MeasurementType GetEntropy(); MeasurementType GetMedian(); + /** + * @brief calculate statistics + */ void CalculateStatistics(); protected: private: HistogramType::Pointer m_Histogram; MeasurementType m_Uniformity, m_UPP, m_Entropy, m_Median; bool m_StatisticsCalculated; }; } #endif diff --git a/Modules/ImageStatistics/mitkHotspotMaskGenerator.cpp b/Modules/ImageStatistics/mitkHotspotMaskGenerator.cpp index 3357d72bc1..a641fe15f1 100644 --- a/Modules/ImageStatistics/mitkHotspotMaskGenerator.cpp +++ b/Modules/ImageStatistics/mitkHotspotMaskGenerator.cpp @@ -1,595 +1,595 @@ #include #include #include #include #include #include "mitkImageAccessByItk.h" #include #include #include namespace mitk { HotspotMaskGenerator::HotspotMaskGenerator(): m_HotspotRadiusinMM(6.2035049089940), // radius of a 1cm3 sphere in mm m_HotspotMustBeCompletelyInsideImage(true), m_Label(1) { m_TimeStep = 0; m_InternalMask = mitk::Image::New(); m_InternalMaskUpdateTime = 0; } void HotspotMaskGenerator::SetImage(mitk::Image::Pointer inputImage) { if (inputImage != m_Image) { m_Image = inputImage; m_ConvolutionImageMaxIndex.set_size(inputImage->GetDimension()); m_ConvolutionImageMinIndex.set_size(inputImage->GetDimension()); this->Modified(); } } void HotspotMaskGenerator::SetMask(MaskGenerator::Pointer mask) { if (mask != m_Mask) { m_Mask = mask; this->Modified(); } } HotspotMaskGenerator::~HotspotMaskGenerator() { } void HotspotMaskGenerator::SetHotspotRadiusInMM(double radiusInMillimeter) { if(radiusInMillimeter != m_HotspotRadiusinMM) { m_HotspotRadiusinMM = radiusInMillimeter; this->Modified(); } } const double& HotspotMaskGenerator::GetHotspotRadiusinMM() const { return m_HotspotRadiusinMM; } bool HotspotMaskGenerator::GetHotspotMustBeCompletelyInsideImage() const { return m_HotspotMustBeCompletelyInsideImage; } void HotspotMaskGenerator::SetHotspotMustBeCompletelyInsideImage(bool mustBeCompletelyInImage) { if (m_HotspotMustBeCompletelyInsideImage != mustBeCompletelyInImage) { m_HotspotMustBeCompletelyInsideImage = mustBeCompletelyInImage; this->Modified(); } } mitk::Image::Pointer HotspotMaskGenerator::GetMask() { if (IsUpdateRequired()) { if ( m_Image.IsNull() ) { throw std::runtime_error( "Error: image empty!" ); } if ( m_TimeStep >= m_Image->GetTimeSteps() ) { throw std::runtime_error( "Error: invalid time step!" ); } mitk::ImageTimeSelector::Pointer imageTimeSelector = mitk::ImageTimeSelector::New(); imageTimeSelector->SetInput( m_Image ); imageTimeSelector->SetTimeNr( m_TimeStep ); imageTimeSelector->UpdateLargestPossibleRegion(); mitk::Image::Pointer timeSliceImage = imageTimeSelector->GetOutput(); m_internalImage = timeSliceImage; m_internalMask2D = nullptr; // is this correct when this variable holds a smart pointer? m_internalMask3D = nullptr; if ( m_Mask != nullptr ) { m_Mask->SetTimeStep(m_TimeStep); mitk::Image::Pointer timeSliceMask = m_Mask->GetMask(); if ( m_internalImage->GetDimension() == 3 ) { CastToItkImage(timeSliceMask, m_internalMask3D); AccessFixedDimensionByItk_2(m_internalImage, CalculateHotspotMask, 3, m_internalMask3D, m_Label); } else if ( m_internalImage->GetDimension() == 2 ) { CastToItkImage(timeSliceMask, m_internalMask2D); AccessFixedDimensionByItk_2(m_internalImage, CalculateHotspotMask, 2, m_internalMask2D, m_Label); } else { throw std::runtime_error( "Error: invalid image dimension" ); } } else { if ( m_internalImage->GetDimension() == 3 ) { AccessFixedDimensionByItk_2(m_internalImage, CalculateHotspotMask, 3, m_internalMask3D, m_Label); } else if ( m_internalImage->GetDimension() == 2 ) { AccessFixedDimensionByItk_2(m_internalImage, CalculateHotspotMask, 2, m_internalMask2D, m_Label); } else { throw std::runtime_error( "Error: invalid image dimension" ); } } this->Modified(); } m_InternalMaskUpdateTime = m_InternalMask->GetMTime(); return m_InternalMask; } - void HotspotMaskGenerator::SetLabel(unsigned int label) + void HotspotMaskGenerator::SetLabel(unsigned short label) { if (label != m_Label) { m_Label = label; this->Modified(); } } vnl_vector HotspotMaskGenerator::GetConvolutionImageMinIndex() { this->GetMask(); // make sure we are up to date return m_ConvolutionImageMinIndex; } vnl_vector HotspotMaskGenerator::GetHotspotIndex() { this->GetMask(); // make sure we are up to date return m_ConvolutionImageMaxIndex; } template HotspotMaskGenerator::ImageExtrema HotspotMaskGenerator::CalculateExtremaWorld( const itk::Image* inputImage, typename itk::Image::Pointer 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 HotspotMaskGenerator::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 > HotspotMaskGenerator::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); mitk::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; mitk::Point3D subVoxelIndexPosition; double distanceSquared = 0.0; typedef itk::ContinuousIndex ContinuousIndexType; for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt) { ContinuousIndexType indexPoint(maskIt.GetIndex()); mitk::Point3D voxelPosition; for (unsigned int dimension = 0; dimension < VImageDimension; ++dimension) { voxelPosition[dimension] = indexPoint[dimension]; } maskValue = 0.0; mitk::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 > HotspotMaskGenerator::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 (m_HotspotMustBeCompletelyInsideImage) { // 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 return convolutionImage; } template < typename TPixel, unsigned int VImageDimension> void HotspotMaskGenerator::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; // this is not very smart. I would rather use a 0 initialized mask (not the case here -> blame CalculateHotspotMask) and find the region where I need to iterate over, then iterate only over the small region for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt) { maskIndex = maskIt.GetIndex(); maskImage->TransformIndexToPhysicalPoint(maskIndex, worldPosition); maskIt.Set( worldPosition.EuclideanDistanceTo(sphereCenter) <= sphereRadiusInMM ? 1 : 0 ); } } template void HotspotMaskGenerator::CalculateHotspotMask(itk::Image* inputImage, typename itk::Image::Pointer maskImage, unsigned int label) { 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 = 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()"); } // if mask image is not defined, create an image of the same size as inputImage and fill it with 1's // there is maybe a better way to do this!? if (maskImage == nullptr) { maskImage = MaskImageType::New(); typename MaskImageType::RegionType maskRegion = inputImage->GetLargestPossibleRegion(); typename MaskImageType::SpacingType maskSpacing = inputImage->GetSpacing(); typename MaskImageType::PointType maskOrigin = inputImage->GetOrigin(); typename MaskImageType::DirectionType maskDirection = inputImage->GetDirection(); maskImage->SetRegions(maskRegion); maskImage->Allocate(); maskImage->SetOrigin(maskOrigin); maskImage->SetSpacing(maskSpacing); maskImage->SetDirection(maskDirection); maskImage->FillBuffer(1); label = 1; } // 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); bool isHotspotDefined = convolutionImageInformation.Defined; if (!isHotspotDefined) { MITK_ERROR << "No origin of hotspot-sphere was calculated!"; m_InternalMask = nullptr; } else { // create a binary mask around the "hotspot" region, fill the shape of a sphere around our hotspot center // typename DuplicatorType::Pointer copyMachine = DuplicatorType::New(); // copyMachine->SetInputImage(inputImage); // copyMachine->Update(); // typename CastFilterType::Pointer caster = CastFilterType::New(); // caster->SetInput( copyMachine->GetOutput() ); // caster->Update(); typename MaskImageType::Pointer hotspotMaskITK = MaskImageType::New(); hotspotMaskITK->SetOrigin(inputImage->GetOrigin()); hotspotMaskITK->SetSpacing(inputImage->GetSpacing()); hotspotMaskITK->SetLargestPossibleRegion(inputImage->GetLargestPossibleRegion()); hotspotMaskITK->SetBufferedRegion(inputImage->GetBufferedRegion()); hotspotMaskITK->SetDirection(inputImage->GetDirection()); hotspotMaskITK->SetNumberOfComponentsPerPixel(inputImage->GetNumberOfComponentsPerPixel()); hotspotMaskITK->Allocate(); hotspotMaskITK->FillBuffer(1); 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); FillHotspotMaskPixels(hotspotMaskITK.GetPointer(), maskCenter, m_HotspotRadiusinMM); //obtain mitk::Image::Pointer from itk::Image mitk::Image::Pointer hotspotMaskAsMITKImage = mitk::GrabItkImageMemory(hotspotMaskITK); m_InternalMask = hotspotMaskAsMITKImage; m_ConvolutionImageMaxIndex = convolutionImageInformation.MaxIndex; m_ConvolutionImageMinIndex = convolutionImageInformation.MinIndex; } } bool HotspotMaskGenerator::IsUpdateRequired() const { unsigned long thisClassTimeStamp = this->GetMTime(); unsigned long internalMaskTimeStamp = m_InternalMask->GetMTime(); unsigned long maskGeneratorTimeStamp = m_Mask->GetMTime(); unsigned long inputImageTimeStamp = m_Image->GetMTime(); if (thisClassTimeStamp > m_InternalMaskUpdateTime) // inputs have changed { return true; } if (m_InternalMaskUpdateTime < maskGeneratorTimeStamp || m_InternalMaskUpdateTime < inputImageTimeStamp) // mask image has changed outside of this class { return true; } if (internalMaskTimeStamp > m_InternalMaskUpdateTime) // internal mask has been changed outside of this class { return true; } return false; } } diff --git a/Modules/ImageStatistics/mitkHotspotMaskGenerator.h b/Modules/ImageStatistics/mitkHotspotMaskGenerator.h index ede8abb84b..5e7886c72f 100644 --- a/Modules/ImageStatistics/mitkHotspotMaskGenerator.h +++ b/Modules/ImageStatistics/mitkHotspotMaskGenerator.h @@ -1,132 +1,165 @@ #ifndef MITKHOTSPOTCALCULATOR_H #define MITKHOTSPOTCALCULATOR_H #include #include #include #include #include #include #include #include namespace mitk { +/** + * @brief The HotspotMaskGenerator class is used when a hotspot has to be found in an image. A hotspot is + * the region of the image where the mean intensity is maximal (=brightest spot). It is usually used in PET scans. + * The identification of the hotspot is done as follows: First a cubic (or circular, if image is 2d) + * mask of predefined size is generated. This mask is then convolved with the input image (in fourier domain). + * The maximum value of the convolved image then corresponds to the hotspot. + * If a maskGenerator is set, only the pixels of the convolved image where the corresponding mask is == @a label + * are searched for the maximum value. + */ class MITKIMAGESTATISTICS_EXPORT HotspotMaskGenerator: public MaskGenerator { public: /** Standard Self typedef */ typedef HotspotMaskGenerator Self; typedef MaskGenerator Superclass; typedef itk::SmartPointer< Self > Pointer; typedef itk::SmartPointer< const Self > ConstPointer; /** Method for creation through the object factory. */ itkNewMacro(Self) /** Runtime information support. */ itkTypeMacro(HotspotMaskGenerator, MaskGenerator) + /** + @brief Set the input image + */ void SetImage(mitk::Image::Pointer inputImage); + /** + @brief Set a mask (can be nullptr if no mask is desired) + */ void SetMask(MaskGenerator::Pointer mask); + /** + @brief Set the radius of the hotspot (in MM) + */ void SetHotspotRadiusInMM(double radiusInMillimeter); const double& GetHotspotRadiusinMM() const; + /** + @brief Define whether the hotspot must be completely inside the image. Default is true + */ void SetHotspotMustBeCompletelyInsideImage(bool hotspotCompletelyInsideImage); bool GetHotspotMustBeCompletelyInsideImage() const; - void SetLabel(unsigned int label); + /** + @brief If a maskGenerator is set, this detemines which mask value is used + */ + void SetLabel(unsigned short label); + /** + @brief Computes and returns the hotspot mask. The hotspot mask has the same size as the input image. The hopspot has value 1, the remaining pixels are set to 0 + */ mitk::Image::Pointer GetMask(); + /** + @brief Returns the image index where the hotspot is located + */ vnl_vector GetHotspotIndex(); + /** + @brief Returns the index where the convolution image is minimal (darkest spot in image) + */ vnl_vector GetConvolutionImageMinIndex(); protected: HotspotMaskGenerator(); ~HotspotMaskGenerator(); 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()) { } }; private: /** \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); /** \brief */ template void CalculateHotspotMask(itk::Image* inputImage, typename itk::Image::Pointer maskImage, unsigned int label); template ImageExtrema CalculateExtremaWorld( const itk::Image* inputImage, typename itk::Image::Pointer maskImage, double neccessaryDistanceToImageBorderInMM, unsigned int label); bool IsUpdateRequired() const; HotspotMaskGenerator(const HotspotMaskGenerator &); HotspotMaskGenerator & operator=(const HotspotMaskGenerator &); mitk::Image::Pointer m_Image; MaskGenerator::Pointer m_Mask; mitk::Image::Pointer m_internalImage; itk::Image::Pointer m_internalMask2D; itk::Image::Pointer m_internalMask3D; double m_HotspotRadiusinMM; bool m_HotspotMustBeCompletelyInsideImage; bool m_HotspotParamsChanged; - unsigned int m_Label; + unsigned short m_Label; vnl_vector m_ConvolutionImageMinIndex, m_ConvolutionImageMaxIndex; unsigned long m_InternalMaskUpdateTime; }; } #endif // MITKHOTSPOTCALCULATOR diff --git a/Modules/ImageStatistics/mitkIgnorePixelMaskGenerator.h b/Modules/ImageStatistics/mitkIgnorePixelMaskGenerator.h index 003d65c66d..0f6a5ffb5b 100644 --- a/Modules/ImageStatistics/mitkIgnorePixelMaskGenerator.h +++ b/Modules/ImageStatistics/mitkIgnorePixelMaskGenerator.h @@ -1,61 +1,73 @@ #ifndef MITKIGNOREPIXELMASKGEN_ #define MITKIGNOREPIXELMASKGEN_ #include #include #include #include #include namespace mitk { +/** + * @brief The IgnorePixelMaskGenerator class is used to generate a mask that is zero for specific pixel values in the input image. + */ class MITKIMAGESTATISTICS_EXPORT IgnorePixelMaskGenerator: public MaskGenerator { public: /** Standard Self typedef */ typedef IgnorePixelMaskGenerator Self; typedef MaskGenerator Superclass; typedef itk::SmartPointer< Self > Pointer; typedef itk::SmartPointer< const Self > ConstPointer; typedef double RealType; /** Method for creation through the object factory. */ itkNewMacro(Self) /** Runtime information support. */ itkTypeMacro(IgnorePixelMaskGenerator, MaskGenerator) + /** + * @brief Set input image + */ void SetImage(mitk::Image::Pointer image); + /** + * @brief The mask will be 0 there inputImage==pixelValue and 1 otherwise + */ void SetIgnoredPixelValue(RealType pixelValue); + /** + * @brief Computes and returns the mask + */ mitk::Image::Pointer GetMask(); protected: IgnorePixelMaskGenerator(): m_IgnoredPixelValue(std::numeric_limits::min()) { m_TimeStep = 0; m_InternalMaskUpdateTime = 0; m_InternalMask = mitk::Image::New(); } ~IgnorePixelMaskGenerator(){} template void InternalCalculateMask(typename itk::Image* image); private: bool IsUpdateRequired() const; mitk::Image::Pointer m_Image; RealType m_IgnoredPixelValue; unsigned long m_InternalMaskUpdateTime; }; } #endif diff --git a/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp b/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp index 8228dad6ce..ab20c3e48f 100644 --- a/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp +++ b/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp @@ -1,598 +1,603 @@ #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 "itkImageFileWriter.h" namespace mitk { void ImageStatisticsCalculator::SetInputImage(mitk::Image::Pointer image) { if (image != m_Image) { m_Image = image; m_StatisticsByTimeStep.resize(m_Image->GetTimeSteps()); m_StatisticsUpdateTimePerTimeStep.resize(m_Image->GetTimeSteps()); std::fill(m_StatisticsUpdateTimePerTimeStep.begin(), m_StatisticsUpdateTimePerTimeStep.end(), 0); this->Modified(); } } void ImageStatisticsCalculator::SetMask(mitk::MaskGenerator::Pointer mask) { if (mask != m_MaskGenerator) { m_MaskGenerator = mask; this->Modified(); } } void ImageStatisticsCalculator::SetSecondaryMask(mitk::MaskGenerator::Pointer mask) { if (mask != m_SecondaryMaskGenerator) { m_SecondaryMaskGenerator = mask; this->Modified(); } } void ImageStatisticsCalculator::SetNBinsForHistogramStatistics(unsigned int nBins) { if (nBins != m_nBinsForHistogramStatistics) { m_nBinsForHistogramStatistics = nBins; this->Modified(); this->m_UseBinSizeOverNBins = false; } if (m_UseBinSizeOverNBins) { this->Modified(); this->m_UseBinSizeOverNBins = false; } } unsigned int ImageStatisticsCalculator::GetNBinsForHistogramStatistics() const { return m_nBinsForHistogramStatistics; } void ImageStatisticsCalculator::SetBinSizeForHistogramStatistics(double binSize) { if (binSize != m_binSizeForHistogramStatistics) { m_binSizeForHistogramStatistics = binSize; this->Modified(); this->m_UseBinSizeOverNBins = true; } if (!m_UseBinSizeOverNBins) { this->Modified(); this->m_UseBinSizeOverNBins = true; } } double ImageStatisticsCalculator::GetBinSizeForHistogramStatistics() const { return m_binSizeForHistogramStatistics; } ImageStatisticsCalculator::StatisticsContainer::Pointer ImageStatisticsCalculator::GetStatistics(unsigned int timeStep, unsigned int label) { if (timeStep >= m_StatisticsByTimeStep.size()) { mitkThrow() << "invalid timeStep in ImageStatisticsCalculator_v2::GetStatistics"; } if (m_Image.IsNull()) { mitkThrow() << "no image"; } if (!m_Image->IsInitialized()) { mitkThrow() << "Image not initialized!"; } if (IsUpdateRequired(timeStep)) { if (m_MaskGenerator.IsNotNull()) { m_MaskGenerator->SetTimeStep(timeStep); m_InternalMask = m_MaskGenerator->GetMask(); } if (m_SecondaryMaskGenerator.IsNotNull()) { m_SecondaryMaskGenerator->SetTimeStep(timeStep); m_SecondaryMask = m_SecondaryMaskGenerator->GetMask(); } ImageTimeSelector::Pointer imgTimeSel = ImageTimeSelector::New(); imgTimeSel->SetInput(m_Image); imgTimeSel->SetTimeNr(timeStep); imgTimeSel->UpdateLargestPossibleRegion(); m_ImageTimeSlice = imgTimeSel->GetOutput(); // Calculate statistics with/without mask if (m_MaskGenerator.IsNull() && m_SecondaryMaskGenerator.IsNull()) { // 1) calculate statistics unmasked: AccessByItk_1(m_ImageTimeSlice, InternalCalculateStatisticsUnmasked, timeStep) } else { // 2) calculate statistics masked AccessByItk_1(m_ImageTimeSlice, InternalCalculateStatisticsMasked, timeStep) } this->Modified(); } m_StatisticsUpdateTimePerTimeStep[timeStep] = m_StatisticsByTimeStep[timeStep][m_StatisticsByTimeStep[timeStep].size()-1]->GetMTime(); for (std::vector::iterator it = m_StatisticsByTimeStep[timeStep].begin(); it != m_StatisticsByTimeStep[timeStep].end(); ++it) { StatisticsContainer::Pointer statCont = *it; if (statCont->GetLabel() == label) { return statCont->Clone(); } } // these lines will ony be executed if the requested label could not be found! MITK_WARN << "Invalid label: " << label << " in time step: " << timeStep; return StatisticsContainer::New(); } template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator::InternalCalculateStatisticsUnmasked( typename itk::Image< TPixel, VImageDimension >* image, unsigned int timeStep) { typedef typename itk::Image< TPixel, VImageDimension > ImageType; typedef typename itk::ExtendedStatisticsImageFilter ImageStatisticsFilterType; typedef typename itk::MinMaxImageFilterWithIndex MinMaxFilterType; StatisticsContainer::Pointer statisticsResult = StatisticsContainer::New(); typename ImageStatisticsFilterType::Pointer statisticsFilter = ImageStatisticsFilterType::New(); statisticsFilter->SetInput(image); statisticsFilter->SetCoordinateTolerance(0.001); statisticsFilter->SetDirectionTolerance(0.001); // TODO: this is single threaded. Implement our own image filter that does this multi threaded // typename itk::MinimumMaximumImageCalculator::Pointer imgMinMaxFilter = itk::MinimumMaximumImageCalculator::New(); // imgMinMaxFilter->SetImage(image); // imgMinMaxFilter->Compute(); vnl_vector minIndex, maxIndex; typename MinMaxFilterType::Pointer minMaxFilter = MinMaxFilterType::New(); minMaxFilter->SetInput(image); minMaxFilter->UpdateLargestPossibleRegion(); typename ImageType::PixelType minval = minMaxFilter->GetMin(); typename ImageType::PixelType maxval = minMaxFilter->GetMax(); typename ImageType::IndexType tmpMinIndex = minMaxFilter->GetMinIndex(); typename ImageType::IndexType tmpMaxIndex = minMaxFilter->GetMaxIndex(); // typename ImageType::IndexType tmpMinIndex = imgMinMaxFilter->GetIndexOfMinimum(); // typename ImageType::IndexType tmpMaxIndex = imgMinMaxFilter->GetIndexOfMaximum(); minIndex.set_size(tmpMaxIndex.GetIndexDimension()); maxIndex.set_size(tmpMaxIndex.GetIndexDimension()); for (unsigned int i=0; i < tmpMaxIndex.GetIndexDimension(); i++) { minIndex[i] = tmpMinIndex[i]; maxIndex[i] = tmpMaxIndex[i]; } statisticsResult->SetMinIndex(minIndex); statisticsResult->SetMaxIndex(maxIndex); //convert m_binSize in m_nBins if necessary unsigned int nBinsForHistogram; if (m_UseBinSizeOverNBins) { nBinsForHistogram = std::max(static_cast(std::ceil(maxval - minval)) / m_binSizeForHistogramStatistics, 10.); // do not allow less than 10 bins } else { nBinsForHistogram = m_nBinsForHistogramStatistics; } statisticsFilter->SetHistogramParameters(nBinsForHistogram, minval, maxval); try { statisticsFilter->Update(); } catch (const itk::ExceptionObject& e) { mitkThrow() << "Image statistics calculation failed due to following ITK Exception: \n " << e.what(); } // no mask, therefore just one label = the whole image m_StatisticsByTimeStep[timeStep].resize(1); statisticsResult->SetLabel(1); statisticsResult->SetN(image->GetLargestPossibleRegion().GetNumberOfPixels()); statisticsResult->SetMean(statisticsFilter->GetMean()); statisticsResult->SetMin(statisticsFilter->GetMinimum()); statisticsResult->SetMax(statisticsFilter->GetMaximum()); statisticsResult->SetVariance(statisticsFilter->GetVariance()); statisticsResult->SetStd(statisticsFilter->GetSigma()); statisticsResult->SetSkewness(statisticsFilter->GetSkewness()); statisticsResult->SetKurtosis(statisticsFilter->GetKurtosis()); statisticsResult->SetRMS(std::sqrt(std::pow(statisticsFilter->GetMean(), 2.) + statisticsFilter->GetVariance())); // variance = sigma^2 statisticsResult->SetMPP(statisticsFilter->GetMPP()); statisticsResult->SetEntropy(statisticsFilter->GetEntropy()); statisticsResult->SetMedian(statisticsFilter->GetMedian()); statisticsResult->SetUniformity(statisticsFilter->GetUniformity()); statisticsResult->SetUPP(statisticsFilter->GetUPP()); statisticsResult->SetHistogram(statisticsFilter->GetHistogram()); m_StatisticsByTimeStep[timeStep][0] = statisticsResult; } template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator::InternalCalculateStatisticsMasked( typename itk::Image< TPixel, VImageDimension >* image, unsigned int timeStep) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< MaskPixelType, VImageDimension > MaskType; typedef typename MaskType::PixelType LabelPixelType; typedef itk::ExtendedLabelStatisticsImageFilter< ImageType, MaskType > ImageStatisticsFilterType; typedef MaskUtilities< TPixel, VImageDimension > MaskUtilType; typedef typename itk::MinMaxLabelImageFilterWithIndex MinMaxLabelFilterType; typedef typename ImageType::PixelType InputImgPixelType; // workaround: if m_SecondaryMaskGenerator ist not null but m_MaskGenerator is! (this is the case if we request a 'ignore zuero valued pixels' // mask in the gui but do not define a primary mask) bool swapMasks = false; if (m_SecondaryMask.IsNotNull() && m_InternalMask.IsNull()) { m_InternalMask = m_SecondaryMask; m_SecondaryMask = nullptr; swapMasks = true; } // maskImage has to have the same dimension as image typename MaskType::Pointer maskImage = MaskType::New(); try { // try to access the pixel values directly (no copying or casting). Only works if mask pixels are of pixelType unsigned short maskImage = ImageToItkImage< MaskPixelType, VImageDimension >(m_InternalMask); } catch (itk::ExceptionObject & e) { // if the pixel type of the mask is not short, then we have to make a copy of m_InternalMask (and cast the values) CastToItkImage(m_InternalMask, maskImage); } // if we have a secondary mask (say a ignoreZeroPixelMask) we need to combine the masks (corresponds to AND) if (m_SecondaryMask.IsNotNull()) { typename MaskType::Pointer secondaryMaskImage = MaskType::New(); secondaryMaskImage = ImageToItkImage< MaskPixelType, VImageDimension >(m_SecondaryMask); // secondary mask should be a ignore zero value pixel mask derived from image. it has to be cropped to the mask region (which may be planar or simply smaller) typename MaskUtilities::Pointer secondaryMaskMaskUtil = MaskUtilities::New(); secondaryMaskMaskUtil->SetImage(secondaryMaskImage.GetPointer()); secondaryMaskMaskUtil->SetMask(maskImage.GetPointer()); typename MaskType::Pointer adaptedSecondaryMaskImage = secondaryMaskMaskUtil->ExtractMaskImageRegion(); typename itk::MaskImageFilter2::Pointer maskFilter = itk::MaskImageFilter2::New(); maskFilter->SetInput1(maskImage); maskFilter->SetInput2(adaptedSecondaryMaskImage); maskFilter->SetMaskingValue(1); // all pixels of maskImage where secondaryMaskImage==1 will be kept, all the others are set to 0 maskFilter->UpdateLargestPossibleRegion(); maskImage = maskFilter->GetOutput(); } typename MaskUtilType::Pointer maskUtil = MaskUtilType::New(); maskUtil->SetImage(image); maskUtil->SetMask(maskImage.GetPointer()); // if mask is smaller than image, extract the image region where the mask is typename ImageType::Pointer adaptedImage = ImageType::New(); adaptedImage = maskUtil->ExtractMaskImageRegion(); // this also checks mask sanity // find min, max, minindex and maxindex typename MinMaxLabelFilterType::Pointer minMaxFilter = MinMaxLabelFilterType::New(); minMaxFilter->SetInput(adaptedImage); minMaxFilter->SetLabelInput(maskImage); minMaxFilter->UpdateLargestPossibleRegion(); // set histogram parameters for each label individually (min/max may be different for each label) typedef typename std::map MapType; typedef typename std::pair PairType; std::vector relevantLabels = minMaxFilter->GetRelevantLabels(); MapType minVals; MapType maxVals; std::map nBins; for (LabelPixelType label:relevantLabels) { minVals.insert(PairType(label, minMaxFilter->GetMin(label))); maxVals.insert(PairType(label, minMaxFilter->GetMax(label))); unsigned int nBinsForHistogram; if (m_UseBinSizeOverNBins) { nBinsForHistogram = std::max(static_cast(std::ceil(minMaxFilter->GetMax(label) - minMaxFilter->GetMin(label))) / m_binSizeForHistogramStatistics, 10.); // do not allow less than 10 bins } else { nBinsForHistogram = m_nBinsForHistogramStatistics; } nBins.insert(typename std::pair(label, nBinsForHistogram)); } typename ImageStatisticsFilterType::Pointer imageStatisticsFilter = ImageStatisticsFilterType::New(); imageStatisticsFilter->SetDirectionTolerance(0.001); imageStatisticsFilter->SetCoordinateTolerance(0.001); imageStatisticsFilter->SetInput(adaptedImage); imageStatisticsFilter->SetLabelInput(maskImage); imageStatisticsFilter->SetHistogramParametersForLabels(nBins, minVals, maxVals); imageStatisticsFilter->Update(); std::list labels = imageStatisticsFilter->GetRelevantLabels(); std::list::iterator it = labels.begin(); m_StatisticsByTimeStep[timeStep].resize(0); while(it != labels.end()) { StatisticsContainer::Pointer statisticsResult = StatisticsContainer::New(); // find min, max, minindex and maxindex // make sure to only look in the masked region, use a masker for this vnl_vector minIndex, maxIndex; typename ImageType::IndexType tmpMinIndex = minMaxFilter->GetMinIndex(*it); typename ImageType::IndexType tmpMaxIndex = minMaxFilter->GetMaxIndex(*it); minIndex.set_size(tmpMaxIndex.GetIndexDimension()); maxIndex.set_size(tmpMaxIndex.GetIndexDimension()); for (unsigned int i=0; i < tmpMaxIndex.GetIndexDimension(); i++) { minIndex[i] = tmpMinIndex[i] + (maskImage->GetOrigin()[i] - image->GetOrigin()[i]) / (double) maskImage->GetSpacing()[i]; maxIndex[i] = tmpMaxIndex[i] + (maskImage->GetOrigin()[i] - image->GetOrigin()[i]) / (double) maskImage->GetSpacing()[i]; } statisticsResult->SetMinIndex(minIndex); statisticsResult->SetMaxIndex(maxIndex); // just debug assert(abs(minMaxFilter->GetMax(*it) - imageStatisticsFilter->GetMaximum(*it)) < mitk::eps); assert(abs(minMaxFilter->GetMin(*it) - imageStatisticsFilter->GetMinimum(*it)) < mitk::eps); statisticsResult->SetN(imageStatisticsFilter->GetSum(*it) / (double) imageStatisticsFilter->GetMean(*it)); statisticsResult->SetMean(imageStatisticsFilter->GetMean(*it)); statisticsResult->SetMin(imageStatisticsFilter->GetMinimum(*it)); statisticsResult->SetMax(imageStatisticsFilter->GetMaximum(*it)); statisticsResult->SetVariance(imageStatisticsFilter->GetVariance(*it)); statisticsResult->SetStd(imageStatisticsFilter->GetSigma(*it)); statisticsResult->SetSkewness(imageStatisticsFilter->GetSkewness(*it)); statisticsResult->SetKurtosis(imageStatisticsFilter->GetKurtosis(*it)); statisticsResult->SetRMS(std::sqrt(std::pow(imageStatisticsFilter->GetMean(*it), 2.) + imageStatisticsFilter->GetVariance(*it))); // variance = sigma^2 statisticsResult->SetMPP(imageStatisticsFilter->GetMPP(*it)); statisticsResult->SetLabel(*it); statisticsResult->SetEntropy(imageStatisticsFilter->GetEntropy(*it)); statisticsResult->SetMedian(imageStatisticsFilter->GetMedian(*it)); statisticsResult->SetUniformity(imageStatisticsFilter->GetUniformity(*it)); statisticsResult->SetUPP(imageStatisticsFilter->GetUPP(*it)); statisticsResult->SetHistogram(imageStatisticsFilter->GetHistogram(*it)); m_StatisticsByTimeStep[timeStep].push_back(statisticsResult); ++it; } // swap maskGenerators back if (swapMasks) { m_SecondaryMask = m_InternalMask; m_InternalMask = nullptr; } } bool ImageStatisticsCalculator::IsUpdateRequired(unsigned int timeStep) const { unsigned long thisClassTimeStamp = this->GetMTime(); unsigned long inputImageTimeStamp = m_Image->GetMTime(); unsigned long statisticsTimeStamp = m_StatisticsUpdateTimePerTimeStep[timeStep]; if (thisClassTimeStamp > statisticsTimeStamp) // inputs have changed { return true; } if (inputImageTimeStamp > statisticsTimeStamp) // image has changed { return true; } if (m_MaskGenerator.IsNotNull()) { unsigned long maskGeneratorTimeStamp = m_MaskGenerator->GetMTime(); if (maskGeneratorTimeStamp > statisticsTimeStamp) // there is a mask generator and it has changed { return true; } } if (m_SecondaryMaskGenerator.IsNotNull()) { unsigned long maskGeneratorTimeStamp = m_SecondaryMaskGenerator->GetMTime(); if (maskGeneratorTimeStamp > statisticsTimeStamp) // there is a secondary mask generator and it has changed { return true; } } return false; } ImageStatisticsCalculator::StatisticsContainer::StatisticsContainer(): m_N(std::numeric_limits::max()), m_Mean(std::numeric_limits::max()), m_Min(std::numeric_limits::min()), m_Max(std::numeric_limits::max()), m_Std(std::numeric_limits::max()), m_Variance(std::numeric_limits::max()), m_Skewness(std::numeric_limits::max()), m_Kurtosis(std::numeric_limits::max()), m_RMS(std::numeric_limits::max()), m_MPP(std::numeric_limits::max()), m_Median(std::numeric_limits::max()), m_Uniformity(std::numeric_limits::max()), m_UPP(std::numeric_limits::max()), m_Entropy(std::numeric_limits::max()) { m_minIndex.set_size(0); m_maxIndex.set_size(0); } ImageStatisticsCalculator::statisticsMapType ImageStatisticsCalculator::StatisticsContainer::GetStatisticsAsMap() { ImageStatisticsCalculator::statisticsMapType statisticsAsMap; statisticsAsMap["N"] = m_N; statisticsAsMap["Mean"] = m_Mean; statisticsAsMap["Min"] = m_Min; statisticsAsMap["Max"] = m_Max; statisticsAsMap["StandardDeviation"] = m_Std; statisticsAsMap["Variance"] = m_Variance; statisticsAsMap["Skewness"] = m_Skewness; statisticsAsMap["Kurtosis"] = m_Kurtosis; statisticsAsMap["RMS"] = m_RMS; statisticsAsMap["MPP"] = m_MPP; statisticsAsMap["Median"] = m_Median; statisticsAsMap["Uniformity"] = m_Uniformity; statisticsAsMap["UPP"] = m_UPP; statisticsAsMap["Entropy"] = m_Entropy; + statisticsAsMap["Label"] = m_Label; return statisticsAsMap; } void ImageStatisticsCalculator::StatisticsContainer::Reset() { m_N = std::numeric_limits::max(); m_Mean = std::numeric_limits::max(); m_Min = std::numeric_limits::max(); m_Max = std::numeric_limits::max(); m_Std = std::numeric_limits::max(); m_Variance = std::numeric_limits::max(); m_Skewness = std::numeric_limits::max(); m_Kurtosis = std::numeric_limits::max(); m_RMS = std::numeric_limits::max(); m_MPP = std::numeric_limits::max(); m_Median = std::numeric_limits::max(); m_Uniformity = std::numeric_limits::max(); m_UPP = std::numeric_limits::max(); m_Entropy = std::numeric_limits::max(); + m_Histogram = HistogramType::New(); + m_minIndex.set_size(0); + m_maxIndex.set_size(0); + m_Label = 0; } void ImageStatisticsCalculator::StatisticsContainer::Print() { ImageStatisticsCalculator::statisticsMapType statMap = this->GetStatisticsAsMap(); // print all map key value pairs // const auto& val:statMap for (auto it = statMap.begin(); it != statMap.end(); ++it) { std::cout << it->first << ": " << it->second << std::endl; } // print the min and max index std::cout << "Min Index:" << std::endl; for (auto it = this->GetMinIndex().begin(); it != this->GetMinIndex().end(); it++) { std::cout << *it << " "; } std::cout << std::endl; // print the min and max index std::cout << "Max Index:" << std::endl; for (auto it = this->GetMaxIndex().begin(); it != this->GetMaxIndex().end(); it++) { std::cout << *it << " "; } std::cout << std::endl; } std::string ImageStatisticsCalculator::StatisticsContainer::GetAsString() { std::string res = ""; ImageStatisticsCalculator::statisticsMapType statMap = this->GetStatisticsAsMap(); // print all map key value pairs // const auto& val:statMap for (auto it = statMap.begin(); it != statMap.end(); ++it) { res += std::string(it->first) + ": " + std::to_string(it->second) + "\n"; } // print the min and max index res += "Min Index:" + std::string("\n"); for (auto it = this->GetMinIndex().begin(); it != this->GetMinIndex().end(); it++) { res += std::to_string(*it) + std::string(" "); } res += "\n"; // print the min and max index res += "Max Index:" + std::string("\n"); for (auto it = this->GetMaxIndex().begin(); it != this->GetMaxIndex().end(); it++) { res += std::to_string(*it) + " "; } res += "\n"; return res; } } diff --git a/Modules/ImageStatistics/mitkImageStatisticsCalculator.h b/Modules/ImageStatistics/mitkImageStatisticsCalculator.h index 26d7c1b652..3a6ebda76d 100644 --- a/Modules/ImageStatistics/mitkImageStatisticsCalculator.h +++ b/Modules/ImageStatistics/mitkImageStatisticsCalculator.h @@ -1,397 +1,413 @@ #ifndef MITKIMAGESTATISTICSCALCULATOR #define MITKIMAGESTATISTICSCALCULATOR #include #include #include #include #include #include #include -/*#define mitkGetRefConstMacro(name, type) \ - const type& Get##name() const \ - { \ - return &m_##name; \ - } \ - \ - type& Get##name() \ - { \ - return &m_##name; \ - } \ -*/ - namespace mitk { class MITKIMAGESTATISTICS_EXPORT ImageStatisticsCalculator: public itk::Object { public: /** Standard Self typedef */ typedef ImageStatisticsCalculator Self; typedef itk::Object Superclass; typedef itk::SmartPointer< Self > Pointer; typedef itk::SmartPointer< const Self > ConstPointer; /** Method for creation through the object factory. */ itkNewMacro(Self) /** Runtime information support. */ itkTypeMacro(ImageStatisticsCalculator_v2, itk::Object) typedef double statisticsValueType; typedef std::map statisticsMapType; typedef itk::Statistics::Histogram HistogramType; typedef unsigned short MaskPixelType; + /**Documentation + @brief Container class for storing the computed image statistics. + + Container class for storing the computed image statistics. Stored statistics are: + - N: number of voxels + - Mean + - MPP (Mean of positive pixels) + - Median + - Skewness + - Kurtosis + - Uniformity + - UPP (Uniformity of positive pixels) + - Variance + - Std (Standard Deviation) + - Min + - Max + - RMS (Root Mean Square) + - Label (if applicable, the label (unsigned short) of the mask the statistics belong to) + - Entropy + + It furthermore stores the following: + - MinIndex (Index of Image where the Minimum is located) + - MaxIndex (Index of Image where the Maximum is located) + - Histogram of Pixel Values*/ class StatisticsContainer: public itk::Object { public: /** Standard Self typedef */ typedef StatisticsContainer Self; typedef itk::Object Superclass; typedef itk::SmartPointer< Self > Pointer; typedef itk::SmartPointer< const Self > ConstPointer; /** Method for creation through the object factory. */ itkNewMacro(Self) - //itkCloneMacro(Self) /** Runtime information support. */ itkTypeMacro(StatisticsContainer, itk::Object) typedef double RealType; + /**Documentation + @brief Returns a std::map containing all real valued statistics stored in this class (= all statistics except minIndex, maxIndex and the histogram)*/ statisticsMapType GetStatisticsAsMap(); + /**Documentation + @brief Deletes all stored values*/ void Reset(); void SetN(long n) { m_N = n; } const long& GetN() const { return m_N; } void SetMean(RealType mean) { m_Mean = mean; } const RealType& GetMean() const { return m_Mean; } void SetVariance(RealType variance) { m_Variance = variance; } const RealType& GetVariance() const { return m_Variance; } void SetStd(RealType std) { m_Std = std; } const RealType& GetStd() const { return m_Std; } void SetMin(RealType minVal) { m_Min = minVal; } const RealType& GetMin() const { return m_Min; } void SetMax(RealType maxVal) { m_Max = maxVal; } const RealType& GetMax() const { return m_Max; } void SetRMS(RealType rms) { m_RMS = rms; } const RealType& GetRMS() const { return m_RMS; } void SetSkewness(RealType skewness) { m_Skewness = skewness; } const RealType& GetSkewness() const { return m_Skewness; } void SetKurtosis(RealType kurtosis) { m_Kurtosis = kurtosis; } const RealType& GetKurtosis() const { return m_Kurtosis; } void SetMPP(RealType mpp) { m_MPP = mpp; } const RealType& GetMPP() const { return m_MPP; } void SetLabel(unsigned int label) { m_Label = label; } const unsigned int& GetLabel() const { return m_Label; } void SetMinIndex(vnl_vector minIndex) { m_minIndex = minIndex; } vnl_vector GetMinIndex() const { return m_minIndex; } void SetMaxIndex(vnl_vector maxIndex) { m_maxIndex = maxIndex; } vnl_vector GetMaxIndex() const { return m_maxIndex; } void SetHistogram(HistogramType::Pointer hist) { if (m_Histogram != hist) { m_Histogram = hist; } } const HistogramType::Pointer GetHistogram() const { return m_Histogram; } void SetEntropy(RealType entropy) { m_Entropy = entropy; } const RealType & GetEntropy() const { return m_Entropy; } void SetMedian(RealType median) { m_Median = median; } const RealType & GetMedian() const { return m_Median; } void SetUniformity(RealType uniformity) { m_Uniformity = uniformity; } const RealType & GetUniformity() const { return m_Uniformity; } void SetUPP(RealType upp) { m_UPP = upp; } const RealType & GetUPP() const { return m_UPP; } + /**Documentation + @brief Creates a StatisticsMapType containing all real valued statistics stored in this class (= all statistics except minIndex, maxIndex and the histogram) and prints its contents to std::cout*/ void Print(); + /**Documentation + @brief Generates a string that contains all real valued statistics stored in this class (= all statistics except minIndex, maxIndex and the histogram)*/ std::string GetAsString(); protected: StatisticsContainer(); -// StatisticsContainer(const StatisticsContainer& other) -// { - -// this->SetEntropy(other.GetEntropy()); -// this->SetKurtosis(other.GetKurtosis()); -// this->SetLabel(other.GetLabel()); -// this->SetMax(other.GetMax()); -// this->SetMin(other.GetMin()); -// this->SetMean(other.GetMean()); -// this->SetMedian(other.GetMedian()); -// this->SetMPP(other.GetMPP()); -// this->SetN(other.GetN()); -// this->SetRMS(other.GetRMS()); -// this->SetSkewness(other.GetSkewness()); -// this->SetStd(other.GetStd()); -// this->SetUniformity(other.GetUniformity()); -// this->SetUPP(other.GetUPP()); -// this->SetVariance(other.GetVariance()); -// this->SetHistogram(other.GetHistogram()->Clone()); -// this->SetMinIndex(other.GetMinIndex()); -// this->SetMaxIndex(other.GetMaxIndex()); -// } - - private: itk::LightObject::Pointer InternalClone() const { itk::LightObject::Pointer ioPtr = Superclass::InternalClone(); Self::Pointer rval = dynamic_cast(ioPtr.GetPointer()); if (rval.IsNull()) { itkExceptionMacro(<< "downcast to type " << "StatisticsContainer" << " failed."); } rval->SetEntropy(this->GetEntropy()); rval->SetKurtosis(this->GetKurtosis()); rval->SetLabel(this->GetLabel()); rval->SetMax(this->GetMax()); rval->SetMin(this->GetMin()); rval->SetMean(this->GetMean()); rval->SetMedian(this->GetMedian()); rval->SetMPP(this->GetMPP()); rval->SetN(this->GetN()); rval->SetRMS(this->GetRMS()); rval->SetSkewness(this->GetSkewness()); rval->SetStd(this->GetStd()); rval->SetUniformity(this->GetUniformity()); rval->SetUPP(this->GetUPP()); rval->SetVariance(this->GetVariance()); rval->SetHistogram(this->GetHistogram()); rval->SetMinIndex(this->GetMinIndex()); rval->SetMaxIndex(this->GetMaxIndex()); return ioPtr; } // not pretty, is temporary long m_N; RealType m_Mean, m_Min, m_Max, m_Std, m_Variance; RealType m_Skewness; RealType m_Kurtosis; RealType m_RMS; RealType m_MPP; vnl_vector m_minIndex, m_maxIndex; RealType m_Median; RealType m_Uniformity; RealType m_UPP; RealType m_Entropy; unsigned int m_Label; HistogramType::Pointer m_Histogram; }; + /**Documentation + @brief Set the image for which the statistics are to be computed.*/ void SetInputImage(mitk::Image::Pointer image); + /**Documentation + @brief Set the mask generator that creates the mask which is to be used to calculate statistics. If no more mask is desired simply set @param mask to nullptr*/ void SetMask(mitk::MaskGenerator::Pointer mask); + /**Documentation + @brief Set this if more than one mask should be applied (for instance if a IgnorePixelValueMask were to be used alongside with a segmentation). + Both masks are combined using pixel wise AND operation. The secondary mask does not have to be the same size than the primary but they need to have some overlap*/ void SetSecondaryMask(mitk::MaskGenerator::Pointer mask); + /**Documentation + @brief Set number of bins to be used for histogram statistics. If Bin size is set after number of bins, bin size will be used instead!*/ void SetNBinsForHistogramStatistics(unsigned int nBins); + /**Documentation + @brief Retrieve the number of bins used for histogram statistics. Careful: The return value does not indicate whether NBins or BinSize is used. + That solely depends on which parameter has been set last.*/ unsigned int GetNBinsForHistogramStatistics() const; + /**Documentation + @brief Set bin size to be used for histogram statistics. If nbins is set after bin size, nbins will be used instead!*/ void SetBinSizeForHistogramStatistics(double binSize); + /**Documentation + @brief Retrieve the bin size for histogram statistics. Careful: The return value does not indicate whether NBins or BinSize is used. + That solely depends on which parameter has been set last.*/ double GetBinSizeForHistogramStatistics() const; + /**Documentation + @brief Returns the statistics for label @a label and timeStep @a timeStep. If these requested statistics are not computed yet the computation is done as well. + For performance reasons, statistics for all labels in the image are computed at once. + */ StatisticsContainer::Pointer GetStatistics(unsigned int timeStep=0, unsigned int label=1); protected: ImageStatisticsCalculator(){ m_nBinsForHistogramStatistics = 100; m_binSizeForHistogramStatistics = 10; m_UseBinSizeOverNBins = false; }; private: template < typename TPixel, unsigned int VImageDimension > void InternalCalculateStatisticsUnmasked( typename itk::Image< TPixel, VImageDimension >* image, unsigned int timeStep); template < typename TPixel, unsigned int VImageDimension > typename HistogramType::Pointer InternalCalculateHistogramUnmasked( typename itk::Image< TPixel, VImageDimension >* image, double minVal, double maxVal); template < typename TPixel, unsigned int VImageDimension > void InternalCalculateStatisticsMasked( typename itk::Image< TPixel, VImageDimension >* image, unsigned int timeStep); bool IsUpdateRequired(unsigned int timeStep) const; std::string GetNameOfClass() { return std::string("ImageStatisticsCalculator_v2"); } mitk::Image::Pointer m_Image; mitk::Image::Pointer m_ImageTimeSlice; mitk::MaskGenerator::Pointer m_MaskGenerator; mitk::Image::Pointer m_InternalMask; mitk::MaskGenerator::Pointer m_SecondaryMaskGenerator; mitk::Image::Pointer m_SecondaryMask; unsigned int m_nBinsForHistogramStatistics; double m_binSizeForHistogramStatistics; bool m_UseBinSizeOverNBins; std::vector> m_StatisticsByTimeStep; std::vector m_StatisticsUpdateTimePerTimeStep; }; } #endif // MITKIMAGESTATISTICSCALCULATOR diff --git a/Modules/ImageStatistics/mitkMaskUtilities.h b/Modules/ImageStatistics/mitkMaskUtilities.h index 7a2c4823ed..23d2b6e9fc 100644 --- a/Modules/ImageStatistics/mitkMaskUtilities.h +++ b/Modules/ImageStatistics/mitkMaskUtilities.h @@ -1,51 +1,69 @@ #ifndef MITKMASKUTIL #define MITKMASKUTIL #include #include #include namespace mitk { +/** + * @brief Utility class for mask operations. It checks whether an image and a mask are compatible (spacing, orientation, etc...) + * and it can also crop an image to the LargestPossibleRegion of the Mask + */ template class MITKIMAGESTATISTICS_EXPORT MaskUtilities: public itk::Object { public: /** Standard Self typedef */ typedef MaskUtilities Self; typedef itk::Object Superclass; typedef itk::SmartPointer< Self > Pointer; typedef itk::SmartPointer< const Self > ConstPointer; /** Method for creation through the object factory. */ itkNewMacro(Self) /** Runtime information support. */ itkTypeMacro(MaskUtilities, itk::Object) typedef itk::Image ImageType; typedef itk::Image MaskType; + /** + * @brief Set image + */ void SetImage(ImageType* image); + + /** + * @brief Set mask + */ void SetMask(MaskType* mask); + /** + * @brief Checks whether mask and image are compatible for joint access (as via iterators). + * Spacing and direction must be the same between the two and they must be aligned. Also, the mask must be completely inside the image + */ bool CheckMaskSanity(); + /** + * @brief Crops the image to the LargestPossibleRegion of the mask + */ typename itk::Image::Pointer ExtractMaskImageRegion(); protected: MaskUtilities(){} ~MaskUtilities(){} private: itk::Image* m_Image; itk::Image* m_Mask; }; } #ifndef ITK_MANUAL_INSTANTIATION #include #endif #endif diff --git a/Modules/ImageStatistics/mitkMinMaxLabelmageFilterWithIndex.h b/Modules/ImageStatistics/mitkMinMaxLabelmageFilterWithIndex.h index 5ae6208965..e0c4158777 100644 --- a/Modules/ImageStatistics/mitkMinMaxLabelmageFilterWithIndex.h +++ b/Modules/ImageStatistics/mitkMinMaxLabelmageFilterWithIndex.h @@ -1,171 +1,177 @@ #ifndef MITK_MINMAXLABELIMAGEFILTERWITHINDEX_H #define MITK_MINMAXLABELIMAGEFILTERWITHINDEX_H #include #include #include #include #include "itksys/hash_map.hxx" namespace itk { template class MinMaxLabelImageFilterWithIndex: public itk::ImageToImageFilter { public: /** Standard Self typedef */ typedef MinMaxLabelImageFilterWithIndex Self; typedef ImageToImageFilter< TInputImage, TInputImage > Superclass; typedef SmartPointer< Self > Pointer; typedef SmartPointer< const Self > ConstPointer; /** Method for creation through the object factory. */ itkNewMacro(Self); /** Runtime information support. */ itkTypeMacro(MinMaxLabelImageFilterWithIndex, ImageToImageFilter); typedef typename TInputImage::RegionType RegionType; typedef typename TInputImage::SizeType SizeType; typedef typename TInputImage::IndexType IndexType; typedef typename TInputImage::PixelType PixelType; typedef typename NumericTraits< PixelType >::RealType RealType; typedef typename TLabelImage::RegionType LabelRegionType; typedef typename TLabelImage::SizeType LabelSizeType; typedef typename TLabelImage::IndexType LabelIndexType; typedef typename TLabelImage::PixelType LabelPixelType; + /** + * @brief The LabelExtrema class is just a container for global min/max values and their indices as well as all min and max values (+indices) of the mask labels + */ class LabelExtrema { public: PixelType m_Min, m_Max; IndexType m_MinIndex, m_MaxIndex; LabelExtrema(): m_Min(std::numeric_limits::max()), m_Max(std::numeric_limits::min()) {} }; typedef typename itksys::hash_map ExtremaMapType; typedef typename ExtremaMapType::iterator ExtremaMapTypeIterator; typedef typename ExtremaMapType::const_iterator ExtremaMapTypeConstIterator; typedef typename ExtremaMapType::value_type MapValueType; PixelType GetMin(LabelPixelType label) const { ExtremaMapTypeConstIterator it = m_LabelExtrema.find(label); if (it == m_LabelExtrema.end()) { MITK_ERROR << "invalid label"; } return (*it).second.m_Min; } PixelType GetMax(LabelPixelType label) const { ExtremaMapTypeConstIterator it = m_LabelExtrema.find(label); if (it == m_LabelExtrema.end()) { MITK_ERROR << "invalid label"; } return (*it).second.m_Max; } + /** + * @brief Returns a std::vector containing all labels for which min and max values (and indices) have been computed + */ std::vector GetRelevantLabels() const { std::vector labels; for (auto&& it:m_LabelExtrema) { labels.push_back(it.first); } return labels; } IndexType GetMinIndex(LabelPixelType label) const { ExtremaMapTypeConstIterator it = m_LabelExtrema.find(label); if (it == m_LabelExtrema.end()) { MITK_ERROR << "invalid label"; } return (*it).second.m_MinIndex; } IndexType GetMaxIndex(LabelPixelType label) const { ExtremaMapTypeConstIterator it = m_LabelExtrema.find(label); if (it == m_LabelExtrema.end()) { MITK_ERROR << "invalid label"; } return (*it).second.m_MaxIndex; } PixelType GetGlobalMin() const { return m_GlobalMin; } PixelType GetGlobalMax() const { return m_GlobalMax; } IndexType GetGlobalMinIndex() const { return m_GlobalMinIndex; } IndexType GetGlobalMaxIndex() const { return m_GlobalMaxIndex; } /** Set the label image */ void SetLabelInput(const TLabelImage *input) { // Process object is not const-correct so the const casting is required. this->SetNthInput( 1, const_cast< TLabelImage * >( input ) ); } /** Get the label image */ const TLabelImage * GetLabelInput() const { return itkDynamicCastInDebugMode< TLabelImage * >( const_cast< DataObject * >( this->ProcessObject::GetInput(1) ) ); } protected: void AllocateOutputs(); void ThreadedGenerateData(const RegionType & outputRegionForThread, ThreadIdType threadId); void BeforeThreadedGenerateData(); void AfterThreadedGenerateData(); private: std::vector m_ThreadExtrema; ExtremaMapType m_LabelExtrema; PixelType m_GlobalMin; PixelType m_GlobalMax; IndexType m_GlobalMinIndex, m_GlobalMaxIndex; }; } #include "mitkMinMaxLabelmageFilterWithIndex.hxx" #endif diff --git a/Modules/ImageStatistics/mitkMultiLabelMaskGenerator.h b/Modules/ImageStatistics/mitkMultiLabelMaskGenerator.h index 59a9ccbb0d..0adef4b488 100644 --- a/Modules/ImageStatistics/mitkMultiLabelMaskGenerator.h +++ b/Modules/ImageStatistics/mitkMultiLabelMaskGenerator.h @@ -1,41 +1,43 @@ #ifndef MITKMULTILABELMASKGENERATOR #define MITKMULTILABELMASKGENERATOR #include #include #include #include namespace mitk { - +/** + * @brief The MultiLabelMaskGenerator class NOT IMPLEMENTED YET! + */ class MITKIMAGESTATISTICS_EXPORT MultiLabelMaskGenerator: public MaskGenerator { public: void setLabelSetImage(mitk::LabelSetImage::Pointer labelSetImage); void addLabel(LabelSetImage::PixelType, std::vector::size_type layer=0); void removeLabel(LabelSetImage::PixelType, std::vector::size_type layer=0); void addLabels(std::pair::size_type, std::vector> labelsToAdd); void removeLabels(std::pair::size_type, std::vector> labelsToAdd); void addLabels(std::vector labels, std::vector::size_type layer=0); void removeLabels(std::vector labels, std::vector::size_type layer=0); void removeLayer(std::vector::size_type layer); mitk::Image::Pointer GetMask(); protected: private: mitk::LabelSetImage::Pointer m_LabelSetImage; std::vector> m_selectedLabels; }; } #endif diff --git a/Modules/ImageStatistics/mitkitkMaskImageFilter.h b/Modules/ImageStatistics/mitkitkMaskImageFilter.h index 13a064818a..24fcfd8ecf 100644 --- a/Modules/ImageStatistics/mitkitkMaskImageFilter.h +++ b/Modules/ImageStatistics/mitkitkMaskImageFilter.h @@ -1,290 +1,289 @@ /*========================================================================= * * Copyright Insight Software Consortium * * Licensed under the Apache License, Version 2.0 (the "License"); * you may not use this file except in compliance with the License. * You may obtain a copy of the License at * * http://www.apache.org/licenses/LICENSE-2.0.txt * * Unless required by applicable law or agreed to in writing, software * distributed under the License is distributed on an "AS IS" BASIS, * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. * See the License for the specific language governing permissions and * limitations under the License. * *=========================================================================*/ #ifndef __itkMaskImageFilter2_h #define __itkMaskImageFilter2_h #include "itkBinaryFunctorImageFilter.h" #include "itkNumericTraits.h" #include "itkVariableLengthVector.h" #include namespace itk { namespace Functor { /** - * \class MaskInput + * \class MaskInput2 * \brief * \ingroup ITKImageIntensity */ template< typename TInput, typename TMask, typename TOutput = TInput > class MaskInput2 { public: typedef typename NumericTraits< TInput >::AccumulateType AccumulatorType; MaskInput2() { m_MaskingValue = NumericTraits< TMask >::ZeroValue(); InitializeOutsideValue( static_cast( ITK_NULLPTR ) ); } ~MaskInput2() {} bool operator!=(const MaskInput2 &) const { return false; } bool operator==(const MaskInput2 & other) const { return !( *this != other ); } inline TOutput operator()(const TInput & A, const TMask & B) const { if ( B == m_MaskingValue ) { return static_cast< TOutput >( A ); } else { return m_OutsideValue; } } /** Method to explicitly set the outside value of the mask */ void SetOutsideValue(const TOutput & outsideValue) { m_OutsideValue = outsideValue; } /** Method to get the outside value of the mask */ const TOutput & GetOutsideValue() const { return m_OutsideValue; } /** Method to explicitly set the masking value */ void SetMaskingValue(const TMask & maskingValue) { m_MaskingValue = maskingValue; } /** Method to get the masking value */ const TMask & GetMaskingValue() const { return m_MaskingValue; } private: template < typename TPixelType > void InitializeOutsideValue( TPixelType * ) { this->m_OutsideValue = NumericTraits< TPixelType >::ZeroValue(); } template < typename TValue > void InitializeOutsideValue( VariableLengthVector * ) { // set the outside value to be of zero length this->m_OutsideValue = VariableLengthVector< TValue >(0); } TOutput m_OutsideValue; TMask m_MaskingValue; }; } /** \class MaskImageFilter * \brief Mask an image with a mask. * * This class is templated over the types of the * input image type, the mask image type and the type of the output image. * Numeric conversions (castings) are done by the C++ defaults. * * The pixel type of the input 2 image must have a valid definition of the * operator != with zero. This condition is required because internally this * filter will perform the operation * * \code - * if pixel_from_mask_image != masking_value + * if pixel_from_mask_image == masking_value * pixel_output_image = pixel_input_image * else * pixel_output_image = outside_value * \endcode * * The pixel from the input 1 is cast to the pixel type of the output image. * * Note that the input and the mask images must be of the same size. * - * \warning Any pixel value other than masking value (0 by default) will not be masked out. * * \sa MaskNegatedImageFilter * \ingroup IntensityImageFilters * \ingroup MultiThreaded * \ingroup ITKImageIntensity * * \wiki * \wikiexample{ImageProcessing/MaskImageFilter,Apply a mask to an image} * \endwiki */ template< typename TInputImage, typename TMaskImage, typename TOutputImage = TInputImage > class MITKIMAGESTATISTICS_EXPORT MaskImageFilter2: public BinaryFunctorImageFilter< TInputImage, TMaskImage, TOutputImage, Functor::MaskInput2< typename TInputImage::PixelType, typename TMaskImage::PixelType, typename TOutputImage::PixelType > > { public: /** Standard class typedefs. */ typedef MaskImageFilter2 Self; typedef BinaryFunctorImageFilter< TInputImage, TMaskImage, TOutputImage, Functor::MaskInput2< typename TInputImage::PixelType, typename TMaskImage::PixelType, typename TOutputImage::PixelType > > Superclass; typedef SmartPointer< Self > Pointer; typedef SmartPointer< const Self > ConstPointer; /** Method for creation through the object factory. */ itkNewMacro(Self); /** Runtime information support. */ itkTypeMacro(MaskImageFilter2, BinaryFunctorImageFilter); /** Typedefs **/ typedef TMaskImage MaskImageType; /** Set/Get the mask image. Pixels set in the mask image will retain * the original value of the input image while pixels not set in * the mask will be set to the "OutsideValue". */ void SetMaskImage(const MaskImageType *maskImage) { // Process object is not const-correct so the const casting is required. this->SetNthInput( 1, const_cast< MaskImageType * >( maskImage ) ); } const MaskImageType * GetMaskImage() { return static_cast(this->ProcessObject::GetInput(1)); } /** Method to explicitly set the outside value of the mask. Defaults to 0 */ void SetOutsideValue(const typename TOutputImage::PixelType & outsideValue) { if ( this->GetOutsideValue() != outsideValue ) { this->Modified(); this->GetFunctor().SetOutsideValue(outsideValue); } } const typename TOutputImage::PixelType & GetOutsideValue() const { return this->GetFunctor().GetOutsideValue(); } /** Method to explicitly set the masking value of the mask. Defaults to 0 */ void SetMaskingValue(const typename TMaskImage::PixelType & maskingValue) { if ( this->GetMaskingValue() != maskingValue ) { this->Modified(); this->GetFunctor().SetMaskingValue(maskingValue); } } /** Method to get the masking value of the mask. */ const typename TMaskImage::PixelType & GetMaskingValue() const { return this->GetFunctor().GetMaskingValue(); } void BeforeThreadedGenerateData() { typedef typename TOutputImage::PixelType PixelType; this->CheckOutsideValue( static_cast(ITK_NULLPTR) ); } #ifdef ITK_USE_CONCEPT_CHECKING // Begin concept checking itkConceptMacro( MaskEqualityComparableCheck, ( Concept::EqualityComparable< typename TMaskImage::PixelType > ) ); itkConceptMacro( InputConvertibleToOutputCheck, ( Concept::Convertible< typename TInputImage::PixelType, typename TOutputImage::PixelType > ) ); // End concept checking #endif protected: MaskImageFilter2() {} virtual ~MaskImageFilter2() {} void PrintSelf(std::ostream & os, Indent indent) const { Superclass::PrintSelf(os, indent); os << indent << "OutsideValue: " << this->GetOutsideValue() << std::endl; } private: MaskImageFilter2(const Self &); //purposely not implemented void operator=(const Self &); //purposely not implemented template < typename TPixelType > void CheckOutsideValue( const TPixelType * ) {} template < typename TValue > void CheckOutsideValue( const VariableLengthVector< TValue > * ) { // Check to see if the outside value contains only zeros. If so, // resize it to have the same number of zeros as the output // image. Otherwise, check that the number of components in the // outside value is the same as the number of components in the // output image. If not, throw an exception. VariableLengthVector< TValue > currentValue = this->GetFunctor().GetOutsideValue(); VariableLengthVector< TValue > zeroVector( currentValue.GetSize() ); zeroVector.Fill( NumericTraits< TValue >::ZeroValue() ); if ( currentValue == zeroVector ) { zeroVector.SetSize( this->GetOutput()->GetVectorLength() ); zeroVector.Fill( NumericTraits< TValue >::ZeroValue() ); this->GetFunctor().SetOutsideValue( zeroVector ); } else if ( this->GetFunctor().GetOutsideValue().GetSize() != this->GetOutput()->GetVectorLength() ) { itkExceptionMacro( << "Number of components in OutsideValue: " << this->GetFunctor().GetOutsideValue().GetSize() << " is not the same as the " << "number of components in the image: " << this->GetOutput()->GetVectorLength()); } } }; } // end namespace itk #endif