diff --git a/Modules/ImageStatistics/mitkExtendedLabelStatisticsImageFilter.h b/Modules/ImageStatistics/mitkExtendedLabelStatisticsImageFilter.h index b3b99841c6..ec06faaf9b 100644 --- a/Modules/ImageStatistics/mitkExtendedLabelStatisticsImageFilter.h +++ b/Modules/ImageStatistics/mitkExtendedLabelStatisticsImageFilter.h @@ -1,350 +1,351 @@ /*=================================================================== 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 seven additional coefficients: * the Skewness, Kurtosis, Uniformity, UPP, MPP, Entropy and Median * * */ template< class TInputImage, class TLabelImage > class ExtendedLabelStatisticsImageFilter : public LabelStatisticsImageFilter< TInputImage, TLabelImage > { public: typedef ExtendedLabelStatisticsImageFilter Self; typedef LabelStatisticsImageFilter < TInputImage, TLabelImage > Superclass; typedef SmartPointer< Self > Pointer; typedef SmartPointer< const Self > ConstPointer; typedef typename Superclass::LabelPixelType LabelPixelType; typedef typename Superclass::RealType RealType; typedef typename Superclass::PixelType PixelType; typedef typename Superclass::MapIterator MapIterator; typedef typename Superclass::BoundingBoxType BoundingBoxType; typedef typename Superclass::RegionType RegionType; typedef itk::Statistics::Histogram HistogramType; itkFactorylessNewMacro( Self ); itkCloneMacro( Self ); itkTypeMacro(ExtendedLabelStatisticsImageFilter, LabelStatisticsImageFilter); /** \class LabelStatistics * \brief Statistics stored per label * \ingroup ITKImageStatistics */ class LabelStatistics { public: // default constructor LabelStatistics() { // initialized to the default values m_Count = NumericTraits< IdentifierType >::ZeroValue(); m_PositivePixelCount = NumericTraits< IdentifierType >::ZeroValue(); m_Sum = NumericTraits< RealType >::ZeroValue(); m_SumOfPositivePixels = NumericTraits< RealType >::ZeroValue(); m_SumOfSquares = NumericTraits< RealType >::ZeroValue(); m_SumOfCubes = NumericTraits< RealType >::ZeroValue(); m_SumOfQuadruples = NumericTraits< RealType >::ZeroValue(); // Set such that the first pixel encountered can be compared m_Minimum = NumericTraits< RealType >::max(); m_Maximum = NumericTraits< RealType >::NonpositiveMin(); // Default these to zero m_Mean = NumericTraits< RealType >::ZeroValue(); m_Sigma = NumericTraits< RealType >::ZeroValue(); m_Variance = NumericTraits< RealType >::ZeroValue(); m_MPP = NumericTraits< RealType >::ZeroValue(); m_Median = NumericTraits< RealType >::ZeroValue(); m_Uniformity = NumericTraits< RealType >::ZeroValue(); m_UPP = NumericTraits< RealType >::ZeroValue(); m_Entropy = NumericTraits< RealType >::ZeroValue(); m_Skewness = NumericTraits< RealType >::ZeroValue(); m_Kurtosis = NumericTraits< RealType >::ZeroValue(); unsigned int imageDimension = itkGetStaticConstMacro(ImageDimension); m_BoundingBox.resize(imageDimension * 2); for ( unsigned int i = 0; i < imageDimension * 2; i += 2 ) { m_BoundingBox[i] = NumericTraits< IndexValueType >::max(); m_BoundingBox[i + 1] = NumericTraits< IndexValueType >::NonpositiveMin(); } m_Histogram = nullptr; } // constructor with histogram enabled LabelStatistics(int size, RealType lowerBound, RealType upperBound) { // initialized to the default values m_Count = NumericTraits< IdentifierType >::ZeroValue(); m_PositivePixelCount = NumericTraits< IdentifierType >::ZeroValue(); m_Sum = NumericTraits< RealType >::ZeroValue(); m_SumOfPositivePixels = NumericTraits< RealType >::ZeroValue(); m_SumOfSquares = NumericTraits< RealType >::ZeroValue(); m_SumOfCubes = NumericTraits< RealType >::ZeroValue(); m_SumOfQuadruples = NumericTraits< RealType >::ZeroValue(); // Set such that the first pixel encountered can be compared m_Minimum = NumericTraits< RealType >::max(); m_Maximum = NumericTraits< RealType >::NonpositiveMin(); // Default these to zero m_Mean = NumericTraits< RealType >::ZeroValue(); m_Sigma = NumericTraits< RealType >::ZeroValue(); m_Variance = NumericTraits< RealType >::ZeroValue(); m_MPP = NumericTraits< RealType >::ZeroValue(); m_Median = NumericTraits< RealType >::ZeroValue(); m_Uniformity = NumericTraits< RealType >::ZeroValue(); m_UPP = NumericTraits< RealType >::ZeroValue(); m_Entropy = NumericTraits< RealType >::ZeroValue(); m_Skewness = NumericTraits< RealType >::ZeroValue(); m_Kurtosis = NumericTraits< RealType >::ZeroValue(); unsigned int imageDimension = itkGetStaticConstMacro(ImageDimension); m_BoundingBox.resize(imageDimension * 2); for ( unsigned int i = 0; i < imageDimension * 2; i += 2 ) { m_BoundingBox[i] = NumericTraits< IndexValueType >::max(); m_BoundingBox[i + 1] = NumericTraits< IndexValueType >::NonpositiveMin(); } // Histogram m_Histogram = HistogramType::New(); typename HistogramType::SizeType hsize; typename HistogramType::MeasurementVectorType lb; typename HistogramType::MeasurementVectorType ub; hsize.SetSize(1); lb.SetSize(1); ub.SetSize(1); m_Histogram->SetMeasurementVectorSize(1); hsize[0] = size; lb[0] = lowerBound; ub[0] = upperBound; m_Histogram->Initialize(hsize, lb, ub); } // need copy constructor because of smart pointer to histogram LabelStatistics(const LabelStatistics & l) { m_Count = l.m_Count; m_Minimum = l.m_Minimum; m_Maximum = l.m_Maximum; m_Mean = l.m_Mean; m_Sum = l.m_Sum; m_SumOfSquares = l.m_SumOfSquares; m_Sigma = l.m_Sigma; m_Variance = l.m_Variance; m_MPP = l.m_MPP; m_Median = l.m_Median; m_Uniformity = l.m_Uniformity; m_UPP = l.m_UPP; m_Entropy = l.m_Entropy; m_Skewness = l.m_Skewness; m_Kurtosis = l.m_Kurtosis; m_BoundingBox = l.m_BoundingBox; m_Histogram = l.m_Histogram; m_SumOfPositivePixels = l.m_SumOfPositivePixels; m_PositivePixelCount = l.m_PositivePixelCount; m_SumOfCubes = l.m_SumOfCubes; m_SumOfQuadruples = l.m_SumOfQuadruples; } // added for completeness LabelStatistics &operator= (const LabelStatistics& l) { if(this != &l) { m_Count = l.m_Count; m_Minimum = l.m_Minimum; m_Maximum = l.m_Maximum; m_Mean = l.m_Mean; m_Sum = l.m_Sum; m_SumOfSquares = l.m_SumOfSquares; m_Sigma = l.m_Sigma; m_Variance = l.m_Variance; m_MPP = l.m_MPP; m_Median = l.m_Median; m_Uniformity = l.m_Uniformity; m_UPP = l.m_UPP; m_Entropy = l.m_Entropy; m_Skewness = l.m_Skewness; m_Kurtosis = l.m_Kurtosis; m_BoundingBox = l.m_BoundingBox; m_Histogram = l.m_Histogram; m_SumOfPositivePixels = l.m_SumOfPositivePixels; m_PositivePixelCount = l.m_PositivePixelCount; m_SumOfCubes = l.m_SumOfCubes; m_SumOfQuadruples = l.m_SumOfQuadruples; } return *this; } IdentifierType m_Count; RealType m_Minimum; RealType m_Maximum; RealType m_Mean; RealType m_Sum; RealType m_SumOfSquares; RealType m_Sigma; RealType m_Variance; RealType m_MPP; RealType m_Median; RealType m_Uniformity; RealType m_UPP; RealType m_Entropy; RealType m_Skewness; RealType m_Kurtosis; IdentifierType m_PositivePixelCount; RealType m_SumOfPositivePixels; RealType m_SumOfCubes; RealType m_SumOfQuadruples; typename Superclass::BoundingBoxType m_BoundingBox; typename HistogramType::Pointer m_Histogram; }; /** Type of the map used to store data per label */ typedef itksys::hash_map< LabelPixelType, LabelStatistics > MapType; typedef typename itksys::hash_map< LabelPixelType, LabelStatistics >::const_iterator StatisticsMapConstIterator; typedef typename itksys::hash_map< LabelPixelType, LabelStatistics >::iterator StatisticsMapIterator; typedef IdentifierType MapSizeType; /** Type of the container used to store valid label values */ typedef std::vector ValidLabelValuesContainerType; /** Return the computed Minimum for a label. */ RealType GetMinimum(LabelPixelType label) const; /** Return the computed Maximum for a label. */ RealType GetMaximum(LabelPixelType label) const; /** Return the computed Mean for a label. */ RealType GetMean(LabelPixelType label) const; /** Return the computed Standard Deviation for a label. */ RealType GetSigma(LabelPixelType label) const; /** Return the computed Variance for a label. */ RealType GetVariance(LabelPixelType label) const; /** Return the computed bounding box for a label. */ BoundingBoxType GetBoundingBox(LabelPixelType label) const; /** Return the computed region. */ RegionType GetRegion(LabelPixelType label) const; /** Return the compute Sum for a label. */ RealType GetSum(LabelPixelType label) const; /** Return the number of pixels for a label. */ MapSizeType GetCount(LabelPixelType label) const; /** Return the histogram for a label */ HistogramType::Pointer GetHistogram(LabelPixelType label) const; /*getter method for the new statistics*/ RealType GetSkewness(LabelPixelType label) const; RealType GetKurtosis(LabelPixelType label) const; RealType GetUniformity( LabelPixelType label) const; RealType GetMedian( LabelPixelType label) const; RealType GetEntropy( LabelPixelType label) const; RealType GetMPP( LabelPixelType label) const; RealType GetUPP( LabelPixelType label) const; bool GetMaskingNonEmpty() const; std::list GetRelevantLabels() const; /** specify global Histogram parameters. If the histogram parameters are set with this function, the same min and max value are used for all histograms. */ void SetHistogramParameters(const int numBins, RealType lowerBound, RealType upperBound); /** specify Histogram parameters for each label individually. Labels in the label image that are not represented in the std::maps here will receive global parameters (if available) */ void SetHistogramParametersForLabels(std::map numBins, std::map lowerBound, std::map upperBound); protected: ExtendedLabelStatisticsImageFilter(): m_GlobalHistogramParametersSet(false), + m_MaskNonEmpty(false), m_LabelHistogramParametersSet(false), m_PreferGlobalHistogramParameters(false) { m_NumBins.set_size(1); } ~ExtendedLabelStatisticsImageFilter() override{} void AfterThreadedGenerateData() override; /** Initialize some accumulators before the threads run. */ void BeforeThreadedGenerateData() override; /** Multi-thread version GenerateData. */ void ThreadedGenerateData(const typename TInputImage::RegionType & outputRegionForThread, ThreadIdType threadId) override; /** Does the specified label exist? Can only be called after a call * a call to Update(). */ bool HasLabel(LabelPixelType label) const { return m_LabelStatistics.find(label) != m_LabelStatistics.end(); } private: std::vector< MapType > m_LabelStatisticsPerThread; MapType m_LabelStatistics; ValidLabelValuesContainerType m_ValidLabelValues; bool m_GlobalHistogramParametersSet; typename HistogramType::SizeType m_NumBins; RealType m_LowerBound; RealType m_UpperBound; bool m_MaskNonEmpty; bool m_LabelHistogramParametersSet; std::map m_LabelMin, m_LabelMax; std::map m_LabelNBins; bool m_PreferGlobalHistogramParameters; }; // end of class } // end namespace itk #ifndef ITK_MANUAL_INSTANTIATION #include "mitkExtendedLabelStatisticsImageFilter.hxx" #endif #endif diff --git a/Modules/ImageStatistics/mitkExtendedStatisticsImageFilter.hxx b/Modules/ImageStatistics/mitkExtendedStatisticsImageFilter.hxx index 9ea49c069b..bb8e3b7e3f 100644 --- a/Modules/ImageStatistics/mitkExtendedStatisticsImageFilter.hxx +++ b/Modules/ImageStatistics/mitkExtendedStatisticsImageFilter.hxx @@ -1,477 +1,474 @@ /*=================================================================== 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" #include namespace itk { template< class TInputImage > ExtendedStatisticsImageFilter< TInputImage >::ExtendedStatisticsImageFilter() : StatisticsImageFilter< TInputImage >(), m_ThreadSum(1), m_SumOfSquares(1), m_SumOfCubes(1), m_SumOfQuadruples(1), m_Count(1), m_PositivePixelCount(1), m_ThreadSumOfPositivePixels(1), m_ThreadMin(1), m_ThreadMax(1), m_UseHistogram(false), m_HistogramCalculated(false) { /* * add the Skewness,Kurtosis,Entropy,Uniformity,MPP, UPP, Median to the other statistical calculated Values * of the mitkStatisticsImageFilter as the 7th to the 13th Output */ for ( int i = 7; i < 14; ++i ) { typename RealObjectType::Pointer output = static_cast< RealObjectType * >( this->MakeOutput(i).GetPointer() ); this->ProcessObject::SetNthOutput( i, output.GetPointer() ); } this->GetSkewnessOutput()->Set( 0.0 ); this->GetKurtosisOutput()->Set( 0.0 ); this->GetEntropyOutput()->Set( -1.0 ); this->GetUniformityOutput()->Set( 0.0 ); this->GetUPPOutput()->Set( 0.0 ); this->GetMPPOutput()->Set( 0.0 ); this->GetMedianOutput()->Set( 0.0 ); m_Histogram = ITK_NULLPTR; m_HistogramPerThread.resize(0); } template< class TInputImage > DataObject::Pointer ExtendedStatisticsImageFilter< TInputImage >::MakeOutput( ProcessObject::DataObjectPointerArraySizeType output) { switch ( output ) { case 7: case 8: case 9: case 10: case 11: case 12: { return RealObjectType::New().GetPointer(); - break; } case 13: { return RealObjectType::New().GetPointer(); - break; } default: { // might as well make an image return Superclass::MakeOutput( output ); - break; } } } template< class TInputImage > typename ExtendedStatisticsImageFilter< TInputImage >::RealObjectType * ExtendedStatisticsImageFilter< TInputImage > ::GetSkewnessOutput() { return static_cast< RealObjectType * >( this->ProcessObject::GetOutput(7) ); } template< class TInputImage > const typename ExtendedStatisticsImageFilter< TInputImage >::RealObjectType * ExtendedStatisticsImageFilter< TInputImage > ::GetSkewnessOutput() const { return static_cast< const RealObjectType * >( this->ProcessObject::GetOutput(7) ); } template< class TInputImage > typename ExtendedStatisticsImageFilter< TInputImage >::RealObjectType * ExtendedStatisticsImageFilter< TInputImage > ::GetKurtosisOutput() { return static_cast< RealObjectType * >( this->ProcessObject::GetOutput(8) ); } template< class TInputImage > const typename ExtendedStatisticsImageFilter< TInputImage >::RealObjectType * ExtendedStatisticsImageFilter< TInputImage > ::GetKurtosisOutput() const { return static_cast< const RealObjectType * >( this->ProcessObject::GetOutput(8) ); } template< class TInputImage > typename ExtendedStatisticsImageFilter< TInputImage >::RealObjectType * ExtendedStatisticsImageFilter< TInputImage > ::GetUniformityOutput() { return static_cast< RealObjectType * >( this->ProcessObject::GetOutput(9) ); } template< class TInputImage > const typename ExtendedStatisticsImageFilter< TInputImage >::RealObjectType * ExtendedStatisticsImageFilter< TInputImage > ::GetUniformityOutput() const { return static_cast< const RealObjectType * >( this->ProcessObject::GetOutput(9) ); } template< class TInputImage > typename ExtendedStatisticsImageFilter< TInputImage >::RealObjectType * ExtendedStatisticsImageFilter< TInputImage > ::GetEntropyOutput() { return static_cast< RealObjectType * >( this->ProcessObject::GetOutput(10) ); } template< class TInputImage > const typename ExtendedStatisticsImageFilter< TInputImage >::RealObjectType * ExtendedStatisticsImageFilter< TInputImage > ::GetEntropyOutput() const { return static_cast< const RealObjectType * >( this->ProcessObject::GetOutput(10) ); } template< class TInputImage > typename ExtendedStatisticsImageFilter< TInputImage >::RealObjectType * ExtendedStatisticsImageFilter< TInputImage > ::GetUPPOutput() { return static_cast< RealObjectType * >( this->ProcessObject::GetOutput(11) ); } template< class TInputImage > const typename ExtendedStatisticsImageFilter< TInputImage >::RealObjectType * ExtendedStatisticsImageFilter< TInputImage > ::GetUPPOutput() const { return static_cast< const RealObjectType * >( this->ProcessObject::GetOutput(11) ); } template< class TInputImage > typename ExtendedStatisticsImageFilter< TInputImage >::RealObjectType * ExtendedStatisticsImageFilter< TInputImage > ::GetMPPOutput() { return static_cast< RealObjectType * >( this->ProcessObject::GetOutput(12) ); } template< class TInputImage > const typename ExtendedStatisticsImageFilter< TInputImage >::RealObjectType * ExtendedStatisticsImageFilter< TInputImage > ::GetMPPOutput() const { return static_cast< const RealObjectType * >( this->ProcessObject::GetOutput(12) ); } template< class TInputImage > typename ExtendedStatisticsImageFilter< TInputImage >::RealObjectType * ExtendedStatisticsImageFilter< TInputImage > ::GetMedianOutput() { return static_cast< RealObjectType * >( this->ProcessObject::GetOutput(13) ); } template< class TInputImage > const typename ExtendedStatisticsImageFilter< TInputImage >::RealObjectType * ExtendedStatisticsImageFilter< TInputImage > ::GetMedianOutput() const { return static_cast< const RealObjectType * >( this->ProcessObject::GetOutput(13) ); } template< typename TInputImage > void ExtendedStatisticsImageFilter< TInputImage > ::BeforeThreadedGenerateData() { ThreadIdType numberOfThreads = this->GetNumberOfThreads(); if (m_UseHistogram) { // Histogram m_Histogram = HistogramType::New(); typename HistogramType::SizeType hsize; typename HistogramType::MeasurementVectorType lb; typename HistogramType::MeasurementVectorType ub; hsize.SetSize(1); lb.SetSize(1); ub.SetSize(1); m_Histogram->SetMeasurementVectorSize(1); hsize[0] = m_NumBins; lb[0] = m_LowerBound; ub[0] = m_UpperBound; m_Histogram->Initialize(hsize, lb, ub); m_HistogramPerThread.resize(numberOfThreads); for (unsigned int i = 0; i < numberOfThreads; i++) { typename HistogramType::Pointer hist = HistogramType::New(); typename HistogramType::SizeType hsize; typename HistogramType::MeasurementVectorType lb; typename HistogramType::MeasurementVectorType ub; hsize.SetSize(1); lb.SetSize(1); ub.SetSize(1); hist->SetMeasurementVectorSize(1); hsize[0] = m_NumBins; lb[0] = m_LowerBound; ub[0] = m_UpperBound; hist->Initialize(hsize, lb, ub); m_HistogramPerThread[i] = hist; } } // Resize the thread temporaries m_Count.SetSize(numberOfThreads); m_SumOfSquares.SetSize(numberOfThreads); m_SumOfCubes.SetSize(numberOfThreads); m_SumOfQuadruples.SetSize(numberOfThreads); m_ThreadSum.SetSize(numberOfThreads); m_ThreadMin.SetSize(numberOfThreads); m_ThreadMax.SetSize(numberOfThreads); m_PositivePixelCount.SetSize(numberOfThreads); m_ThreadSumOfPositivePixels.SetSize(numberOfThreads); // Initialize the temporaries m_Count.Fill(NumericTraits< SizeValueType >::Zero); m_ThreadSum.Fill(NumericTraits< RealType >::Zero); m_SumOfSquares.Fill(NumericTraits< RealType >::Zero); m_SumOfCubes.Fill(NumericTraits< RealType >::Zero); m_SumOfQuadruples.Fill(NumericTraits< RealType >::Zero); m_ThreadMin.Fill( NumericTraits< PixelType >::max() ); m_ThreadMax.Fill( NumericTraits< PixelType >::NonpositiveMin() ); m_PositivePixelCount.Fill(NumericTraits< SizeValueType >::Zero); m_ThreadSumOfPositivePixels.Fill(NumericTraits< SizeValueType >::Zero); } template< typename TInputImage > void ExtendedStatisticsImageFilter< TInputImage > ::ThreadedGenerateData(const typename StatisticsImageFilter::RegionType & outputRegionForThread, ThreadIdType threadId) { const SizeValueType size0 = outputRegionForThread.GetSize(0); if( size0 == 0) { return; } RealType realValue; PixelType value; typename HistogramType::IndexType histogramIndex(1); typename HistogramType::MeasurementVectorType histogramMeasurement(1); RealType sum = NumericTraits< RealType >::ZeroValue(); RealType sumOfPositivePixels = NumericTraits< RealType >::ZeroValue(); RealType sumOfSquares = NumericTraits< RealType >::ZeroValue(); RealType sumOfCubes = NumericTraits< RealType >::ZeroValue(); RealType sumOfQuadruples = NumericTraits< RealType >::ZeroValue(); SizeValueType count = NumericTraits< SizeValueType >::ZeroValue(); SizeValueType countOfPositivePixels = NumericTraits< SizeValueType >::ZeroValue(); PixelType min = NumericTraits< PixelType >::max(); PixelType max = NumericTraits< PixelType >::NonpositiveMin(); ImageScanlineConstIterator< TInputImage > it (this->GetInput(), outputRegionForThread); // support progress methods/callbacks const size_t numberOfLinesToProcess = outputRegionForThread.GetNumberOfPixels() / size0; ProgressReporter progress( this, threadId, numberOfLinesToProcess ); // do the work while ( !it.IsAtEnd() ) { while ( !it.IsAtEndOfLine() ) { value = it.Get(); realValue = static_cast< RealType >( value ); if (m_UseHistogram) { histogramMeasurement[0] = value; m_HistogramPerThread[threadId]->GetIndex(histogramMeasurement, histogramIndex); m_HistogramPerThread[threadId]->IncreaseFrequencyOfIndex(histogramIndex, 1); } if ( value < min ) { min = value; } if ( value > max ) { max = value; } if (value > 0) { sumOfPositivePixels += realValue; ++countOfPositivePixels; } sum += realValue; sumOfSquares += ( realValue * realValue ); sumOfCubes += std::pow(realValue, 3.); sumOfQuadruples += std::pow(realValue, 4.); ++count; ++it; } it.NextLine(); progress.CompletedPixel(); } m_ThreadSum[threadId] = sum; m_SumOfSquares[threadId] = sumOfSquares; m_Count[threadId] = count; m_ThreadMin[threadId] = min; m_ThreadMax[threadId] = max; m_ThreadSumOfPositivePixels[threadId] = sumOfPositivePixels; m_PositivePixelCount[threadId] = countOfPositivePixels; m_SumOfCubes[threadId] = sumOfCubes; m_SumOfQuadruples[threadId] = sumOfQuadruples; } template< class TInputImage > void ExtendedStatisticsImageFilter< TInputImage > ::AfterThreadedGenerateData() { ThreadIdType i; SizeValueType count = 0; SizeValueType countOfPositivePixels = 0; RealType sumOfSquares = 0; RealType sumOfCubes = 0; RealType sumOfQuadruples = 0; RealType sumOfPositivePixels = 0; ThreadIdType numberOfThreads = this->GetNumberOfThreads(); PixelType minimum; PixelType maximum; RealType mean; RealType meanOfPositivePixels; RealType sigma; RealType variance; RealType skewness; RealType kurtosis; RealType sum; sum = sumOfSquares = sumOfCubes = sumOfQuadruples = NumericTraits< RealType >::ZeroValue(); count = countOfPositivePixels = 0; // Find the min/max over all threads and accumulate count, sum and // sum of squares minimum = NumericTraits< PixelType >::max(); maximum = NumericTraits< PixelType >::NonpositiveMin(); for ( i = 0; i < numberOfThreads; i++ ) { count += m_Count[i]; sum += m_ThreadSum[i]; sumOfSquares += m_SumOfSquares[i]; sumOfCubes += m_SumOfCubes[i]; sumOfQuadruples += m_SumOfQuadruples[i]; sumOfPositivePixels += m_ThreadSumOfPositivePixels[i]; countOfPositivePixels += m_PositivePixelCount[i]; if ( m_ThreadMin[i] < minimum ) { minimum = m_ThreadMin[i]; } if ( m_ThreadMax[i] > maximum ) { maximum = m_ThreadMax[i]; } // if enabled, update the histogram for this label if ( m_UseHistogram ) { typename HistogramType::IndexType index; index.SetSize(1); for ( int bin = 0; bin < m_NumBins; ++bin ) { index[0] = bin; m_Histogram->IncreaseFrequency( bin, m_HistogramPerThread[i]->GetFrequency(bin) ); } } } // compute statistics mean = sum / static_cast< RealType >( count ); meanOfPositivePixels = sumOfPositivePixels / static_cast< RealType >( countOfPositivePixels ); // unbiased estimate variance = ( sumOfSquares - ( sum * sum / static_cast< RealType >( count ) ) ) / ( static_cast< RealType >( count ) ); RealType secondMoment, thirdMoment, fourthMoment; secondMoment = sumOfSquares / static_cast< RealType >( count ); thirdMoment = sumOfCubes / static_cast< RealType >( count ); fourthMoment = sumOfQuadruples / static_cast< RealType >( count ); skewness = (thirdMoment - 3. * secondMoment * mean + 2. * std::pow(mean, 3.)) / std::pow(secondMoment - std::pow(mean, 2.), 1.5); // see http://www.boost.org/doc/libs/1_51_0/doc/html/boost/accumulators/impl/skewness_impl.html kurtosis = (fourthMoment - 4. * thirdMoment * mean + 6. * secondMoment * std::pow(mean, 2.) - 3. * std::pow(mean, 4.)) / std::pow(secondMoment - std::pow(mean, 2.), 2.); // see http://www.boost.org/doc/libs/1_46_1/doc/html/boost/accumulators/impl/kurtosis_impl.html, dropped -3 sigma = std::sqrt(variance); // Set the outputs this->GetMinimumOutput()->Set(minimum); this->GetMaximumOutput()->Set(maximum); this->GetMeanOutput()->Set(mean); this->GetSigmaOutput()->Set(sigma); this->GetVarianceOutput()->Set(variance); this->GetSumOutput()->Set(sum); this->GetKurtosisOutput()->Set(kurtosis); this->GetSkewnessOutput()->Set(skewness); this->GetMPPOutput()->Set(meanOfPositivePixels); if (m_UseHistogram) { mitk::HistogramStatisticsCalculator histStatCalc; histStatCalc.SetHistogram(m_Histogram); histStatCalc.CalculateStatistics(); this->GetUniformityOutput()->Set(histStatCalc.GetUniformity()); this->GetUPPOutput()->Set(histStatCalc.GetUPP()); this->GetEntropyOutput()->Set(histStatCalc.GetEntropy()); this->GetMedianOutput()->Set(histStatCalc.GetMedian()); } m_HistogramCalculated = true; } template< typename TInputImage > void ExtendedStatisticsImageFilter< TInputImage > ::SetHistogramParameters(const int numBins, RealType lowerBound, RealType upperBound) { m_NumBins = numBins; m_LowerBound = lowerBound; m_UpperBound = upperBound; m_UseHistogram = true; } } #endif diff --git a/Modules/ImageStatistics/mitkHistogramStatisticsCalculator.cpp b/Modules/ImageStatistics/mitkHistogramStatisticsCalculator.cpp index b7af6a479d..be86834261 100644 --- a/Modules/ImageStatistics/mitkHistogramStatisticsCalculator.cpp +++ b/Modules/ImageStatistics/mitkHistogramStatisticsCalculator.cpp @@ -1,108 +1,107 @@ #include #include namespace mitk { HistogramStatisticsCalculator::HistogramStatisticsCalculator(): - m_StatisticsCalculated(false) + m_Uniformity(0), m_UPP(0), m_Entropy(0), m_Median(0), m_StatisticsCalculated(false) { } void HistogramStatisticsCalculator::SetHistogram(HistogramType::Pointer histogram) { if (m_Histogram != histogram) { m_Histogram = histogram; m_StatisticsCalculated = false; } } HistogramStatisticsCalculator::MeasurementType HistogramStatisticsCalculator::GetEntropy() { if (!m_StatisticsCalculated) { MITK_WARN("Statistics have not yet been calculated, running calculation now..."); CalculateStatistics(); } return m_Entropy; } HistogramStatisticsCalculator::MeasurementType HistogramStatisticsCalculator::GetMedian() { if (!m_StatisticsCalculated) { MITK_WARN("Statistics have not yet been calculated, running calculation now..."); CalculateStatistics(); } return m_Median; } HistogramStatisticsCalculator::MeasurementType HistogramStatisticsCalculator::GetUniformity() { if (!m_StatisticsCalculated) { MITK_WARN("Statistics have not yet been calculated, running calculation now..."); CalculateStatistics(); } return m_Uniformity; } HistogramStatisticsCalculator::MeasurementType HistogramStatisticsCalculator::GetUPP() { if (!m_StatisticsCalculated) { MITK_WARN("Statistics have not yet been calculated, running calculation now..."); CalculateStatistics(); } return m_UPP; } void HistogramStatisticsCalculator::CalculateStatistics() { if (m_Histogram.IsNull()) { throw std::runtime_error("Histogram not set in HistogramStatisticsCalculator::CalculateStatistics()"); } unsigned int nBins = m_Histogram->GetSize()[0]; m_Uniformity = 0; m_Entropy = 0; m_UPP = 0; m_Median = 0; MeasurementType cumulativeProbability = 0.0; - MeasurementType partialProbability; bool medianFound(false); for (unsigned int i = 0; i < nBins; i++) { - partialProbability = m_Histogram->GetFrequency(i, 0) / double( m_Histogram->GetTotalFrequency() ); + MeasurementType partialProbability = m_Histogram->GetFrequency(i, 0) / double( m_Histogram->GetTotalFrequency() ); cumulativeProbability += partialProbability; if (partialProbability != 0) { m_Entropy -= partialProbability * (std::log2( partialProbability )); m_Uniformity += std::pow(partialProbability, 2); if (m_Histogram->GetMeasurement(i, 0) > 0) { m_UPP += std::pow(partialProbability, 2); } } if (cumulativeProbability >= 0.5 && !medianFound) { MeasurementType binMin = m_Histogram->GetBinMin(0, i); MeasurementType binMax = m_Histogram->GetBinMax(0, i); m_Median = (binMax + binMin) / 2.0; medianFound = true; } } m_StatisticsCalculated = true; } } diff --git a/Modules/ImageStatistics/mitkHotspotMaskGenerator.cpp b/Modules/ImageStatistics/mitkHotspotMaskGenerator.cpp index 5f8a2f6844..9e594ba377 100644 --- a/Modules/ImageStatistics/mitkHotspotMaskGenerator.cpp +++ b/Modules/ImageStatistics/mitkHotspotMaskGenerator.cpp @@ -1,602 +1,601 @@ #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::SetInputImage(mitk::Image::Pointer inputImage) { if (inputImage != m_inputImage) { m_inputImage = 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_inputImage.IsNull() ) { throw std::runtime_error( "Error: image empty!" ); } if ( m_TimeStep >= m_inputImage->GetTimeSteps() ) { throw std::runtime_error( "Error: invalid time step!" ); } mitk::ImageTimeSelector::Pointer imageTimeSelector = mitk::ImageTimeSelector::New(); imageTimeSelector->SetInput( m_inputImage ); 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::SetTimeStep(unsigned int timeStep) { if (m_TimeStep != timeStep) { m_TimeStep = timeStep; } } 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) { itk::IndexValueType 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; + double 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< 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_inputImage->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 d416012f45..88ec126506 100644 --- a/Modules/ImageStatistics/mitkHotspotMaskGenerator.h +++ b/Modules/ImageStatistics/mitkHotspotMaskGenerator.h @@ -1,170 +1,169 @@ #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. Required for this class */ void SetInputImage(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; /** @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() override; /** @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(); /** * @brief SetTimeStep is used to set the time step for which the mask is to be generated * @param timeStep */ void SetTimeStep(unsigned int timeStep) override; protected: HotspotMaskGenerator(); ~HotspotMaskGenerator() override; 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 &); 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 short m_Label; vnl_vector m_ConvolutionImageMinIndex, m_ConvolutionImageMaxIndex; unsigned long m_InternalMaskUpdateTime; }; } #endif // MITKHOTSPOTCALCULATOR diff --git a/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp b/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp index 5ee2941914..3a5dc97fdd 100644 --- a/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp +++ b/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp @@ -1,629 +1,630 @@ #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_MaskGenerator->GetReferenceImage().IsNotNull()) { m_InternalImageForStatistics = m_MaskGenerator->GetReferenceImage(); } else { m_InternalImageForStatistics = m_Image; } } else { m_InternalImageForStatistics = m_Image; } if (m_SecondaryMaskGenerator.IsNotNull()) { m_SecondaryMaskGenerator->SetTimeStep(timeStep); m_SecondaryMask = m_SecondaryMaskGenerator->GetMask(); } ImageTimeSelector::Pointer imgTimeSel = ImageTimeSelector::New(); imgTimeSel->SetInput(m_InternalImageForStatistics); imgTimeSel->SetTimeNr(timeStep); imgTimeSel->UpdateLargestPossibleRegion(); 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 (auto 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 (const itk::ExceptionObject &) { // if the pixel type of the mask is not short, then we have to make a copy of m_InternalMask (and cast the values) CastToItkImage(m_InternalMask, maskImage); } // if we have a secondary mask (say a ignoreZeroPixelMask) we need to combine the masks (corresponds to AND) if (m_SecondaryMask.IsNotNull()) { // dirty workaround for a bug when pf mask + any other mask is used in conjunction. We need a proper fix for this (Fabian Isensee is responsible and probably working on it!) if (m_InternalMask->GetDimension() == 2 && (m_SecondaryMask->GetDimension() == 3 || m_SecondaryMask->GetDimension() == 4)) { mitk::Image::Pointer old_img = m_SecondaryMaskGenerator->GetReferenceImage(); m_SecondaryMaskGenerator->SetInputImage(m_MaskGenerator->GetReferenceImage()); m_SecondaryMask = m_SecondaryMaskGenerator->GetMask(); m_SecondaryMaskGenerator->SetInputImage(old_img); } typename MaskType::Pointer secondaryMaskImage = MaskType::New(); secondaryMaskImage = ImageToItkImage< 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(); auto 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; mitk::Point3D worldCoordinateMin; mitk::Point3D worldCoordinateMax; mitk::Point3D indexCoordinateMin; mitk::Point3D indexCoordinateMax; m_InternalImageForStatistics->GetGeometry()->IndexToWorld(minMaxFilter->GetMinIndex(*it), worldCoordinateMin); m_InternalImageForStatistics->GetGeometry()->IndexToWorld(minMaxFilter->GetMaxIndex(*it), worldCoordinateMax); m_Image->GetGeometry()->WorldToIndex(worldCoordinateMin, indexCoordinateMin); m_Image->GetGeometry()->WorldToIndex(worldCoordinateMax, indexCoordinateMax); minIndex.set_size(3); maxIndex.set_size(3); //for (unsigned int i=0; i < tmpMaxIndex.GetIndexDimension(); i++) for (unsigned int i=0; i < 3; i++) { minIndex[i] = indexCoordinateMin[i]; maxIndex[i] = indexCoordinateMax[i]; } statisticsResult->SetMinIndex(minIndex); statisticsResult->SetMaxIndex(maxIndex); assert(std::abs(minMaxFilter->GetMax(*it) - imageStatisticsFilter->GetMaximum(*it)) < mitk::eps); assert(std::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(0), m_Mean(nan("")), m_Min(nan("")), m_Max(nan("")), m_Std(nan("")), m_Variance(nan("")), m_Skewness(nan("")), m_Kurtosis(nan("")), m_RMS(nan("")), m_MPP(nan("")), m_Median(nan("")), m_Uniformity(nan("")), m_UPP(nan("")), - m_Entropy(nan("")) + m_Entropy(nan("")), + m_Label(0) { 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 = 0; m_Mean = nan(""); m_Min = nan(""); m_Max = nan(""); m_Std = nan(""); m_Variance = nan(""); m_Skewness = nan(""); m_Kurtosis = nan(""); m_RMS = nan(""); m_MPP = nan(""); m_Median = nan(""); m_Uniformity = nan(""); m_UPP = nan(""); m_Entropy = nan(""); 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/mitkMaskUtilities.h b/Modules/ImageStatistics/mitkMaskUtilities.h index 73bd79c187..79d93ab86b 100644 --- a/Modules/ImageStatistics/mitkMaskUtilities.h +++ b/Modules/ImageStatistics/mitkMaskUtilities.h @@ -1,69 +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(): m_Image(nullptr), m_Mask(nullptr){} ~MaskUtilities() override{} private: itk::Image* m_Image; itk::Image* m_Mask; }; } #ifndef ITK_MANUAL_INSTANTIATION #include #endif #endif diff --git a/Modules/ImageStatistics/mitkMinMaxImageFilterWithIndex.hxx b/Modules/ImageStatistics/mitkMinMaxImageFilterWithIndex.hxx index c925cb5aa4..77b4510f4f 100644 --- a/Modules/ImageStatistics/mitkMinMaxImageFilterWithIndex.hxx +++ b/Modules/ImageStatistics/mitkMinMaxImageFilterWithIndex.hxx @@ -1,106 +1,109 @@ #ifndef MITK_MinMaxImageFilterWithIndex_HXX #define MITK_MinMaxImageFilterWithIndex_HXX #include #include namespace itk { template< typename TInputImage > void MinMaxImageFilterWithIndex< TInputImage >::AllocateOutputs() { // Pass the input through as the output typename TInputImage::Pointer image = const_cast< TInputImage * >( this->GetInput() ); this->GraftOutput(image); // Nothing that needs to be allocated for the remaining outputs } template< typename TInputImage > void MinMaxImageFilterWithIndex< TInputImage >::ThreadedGenerateData(const RegionType & outputRegionForThread, ThreadIdType threadId) { const SizeValueType size0 = outputRegionForThread.GetSize(0); if( size0 == 0) { return; } PixelType value; PixelType threadMin, threadMax; - IndexType threadMinIndex, threadMaxIndex; + IndexType threadMinIndex; + IndexType threadMaxIndex; + threadMinIndex.Fill(0); + threadMaxIndex.Fill(0); threadMin = std::numeric_limits::max(); threadMax = std::numeric_limits::min(); ImageRegionConstIteratorWithIndex< TInputImage > it (this->GetInput(), outputRegionForThread); // do the work while ( !it.IsAtEnd() ) { value = it.Get(); if (value < threadMin) { threadMin = value; threadMinIndex = it.GetIndex(); } if (value > threadMax) { threadMax = value; threadMaxIndex = it.GetIndex(); } ++it; } m_ThreadMax[threadId] = threadMax; m_ThreadMin[threadId] = threadMin; m_ThreadMaxIndex[threadId] = threadMaxIndex; m_ThreadMinIndex[threadId] = threadMinIndex; } template< typename TInputImage > void MinMaxImageFilterWithIndex< TInputImage >::BeforeThreadedGenerateData() { ThreadIdType numberOfThreads = this->GetNumberOfThreads(); m_ThreadMin.resize(numberOfThreads); m_ThreadMax.resize(numberOfThreads); m_ThreadMinIndex.resize(numberOfThreads); m_ThreadMaxIndex.resize(numberOfThreads); for (unsigned int i =0; i < numberOfThreads; i++) { m_ThreadMin[i] = std::numeric_limits::max(); m_ThreadMax[i] = std::numeric_limits::min(); } m_Min = std::numeric_limits::max(); m_Max = std::numeric_limits::min(); } template< typename TInputImage > void MinMaxImageFilterWithIndex< TInputImage >::AfterThreadedGenerateData() { ThreadIdType numberOfThreads = this->GetNumberOfThreads(); for (ThreadIdType i = 0; i < numberOfThreads; i++) { if (m_ThreadMin[i] < m_Min) { m_Min = m_ThreadMin[i]; m_MinIndex = m_ThreadMinIndex[i]; } if (m_ThreadMax[i] > m_Max) { m_Max = m_ThreadMax[i]; m_MaxIndex = m_ThreadMaxIndex[i]; } } } } #endif diff --git a/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.h b/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.h index f28688f3e3..dcc30b9e70 100644 --- a/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.h +++ b/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.h @@ -1,131 +1,132 @@ #ifndef MITKPLANARFIGUREMASKGENERATOR #define MITKPLANARFIGUREMASKGENERATOR #include #include #include #include #include #include #include #include #include #include namespace mitk { /** * \class PlanarFigureMaskGenerator * \brief Derived from MaskGenerator. This class is used to convert a mitk::PlanarFigure into a binary image mask */ class MITKIMAGESTATISTICS_EXPORT PlanarFigureMaskGenerator: public MaskGenerator { public: /** Standard Self typedef */ typedef PlanarFigureMaskGenerator 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(PlanarFigureMaskGenerator, MaskGenerator) /** * @brief GetMask Computes and returns the mask * @return mitk::Image::Pointer of the generated mask */ mitk::Image::Pointer GetMask() override; void SetPlanarFigure(mitk::PlanarFigure::Pointer planarFigure); mitk::Image::Pointer GetReferenceImage() override; /** * @brief SetTimeStep is used to set the time step for which the mask is to be generated * @param timeStep */ void SetTimeStep(unsigned int timeStep) override; protected: PlanarFigureMaskGenerator():Superclass(){ m_InternalMaskUpdateTime = 0; m_InternalMask = mitk::Image::New(); m_ReferenceImage = nullptr; + m_PlanarFigureAxis = 0; } private: void CalculateMask(); template < typename TPixel, unsigned int VImageDimension > void InternalCalculateMaskFromPlanarFigure( const itk::Image< TPixel, VImageDimension > *image, unsigned int axis ); template < typename TPixel, unsigned int VImageDimension > void InternalCalculateMaskFromOpenPlanarFigure( const itk::Image< TPixel, VImageDimension > *image, unsigned int axis ); mitk::Image::Pointer extract2DImageSlice(unsigned int axis, unsigned int slice); bool GetPrincipalAxis(const BaseGeometry *geometry, Vector3D vector, unsigned int &axis ); /** 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()); } bool IsUpdateRequired() const; mitk::PlanarFigure::Pointer m_PlanarFigure; itk::Image::Pointer m_InternalITKImageMask2D; mitk::Image::Pointer m_InternalTimeSliceImage; mitk::Image::Pointer m_ReferenceImage; unsigned int m_PlanarFigureAxis; unsigned long m_InternalMaskUpdateTime; }; } #endif // MITKPLANARFIGUREMASKGENERATOR diff --git a/Modules/ImageStatistics/mitkPointSetStatisticsCalculator.h b/Modules/ImageStatistics/mitkPointSetStatisticsCalculator.h index c612d2e152..9495c501b7 100644 --- a/Modules/ImageStatistics/mitkPointSetStatisticsCalculator.h +++ b/Modules/ImageStatistics/mitkPointSetStatisticsCalculator.h @@ -1,122 +1,122 @@ /*=================================================================== 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_PointSetSTATISTICSCALCULATOR_H #define _MITK_PointSetSTATISTICSCALCULATOR_H #include #include #include namespace mitk { /** * \brief Class for calculating statistics (like standard derivation, RMS, mean, etc.) for a PointSet. */ class MITKIMAGESTATISTICS_EXPORT PointSetStatisticsCalculator : public itk::Object { public: mitkClassMacroItkParent( PointSetStatisticsCalculator, itk::Object ); itkFactorylessNewMacro(Self) itkCloneMacro(Self) mitkNewMacro1Param(PointSetStatisticsCalculator,mitk::PointSet::Pointer); /** @brief Sets the point set which will be analysed. */ void SetPointSet(mitk::PointSet::Pointer pSet); /** @return Returns the mean position of the analysed point set (only valid navigation data). Returns [0;0;0] if there is no valid navigation data.*/ mitk::Point3D GetPositionMean(); /** @return Returns the standard derivation of each component (x, y and z) of the analysed point set (only valid navigation data). Returns [0;0;0] if there is no valid navigation data.*/ mitk::Vector3D GetPositionStandardDeviation(); /** @return Returns the sample standard derivation of each component (x, y and z) of the analysed point set (only valid navigation data). Returns [0;0;0] if there is no valid navigation data.*/ mitk::Vector3D GetPositionSampleStandardDeviation(); /** @return Returns the mean distance to the mean postion (=mean error) of the analysed point set (only valid navigation data). Returns 0 if there is no valid navigation data. */ double GetPositionErrorMean(); /** @return Returns the standard derivation of the errors of all positions of the analysed point set (only valid navigation data). Returns 0 if there is no valid navigation data. */ double GetPositionErrorStandardDeviation(); /** @return Returns the sample standard derivation of the errors of all positions of the analysed point set (only valid navigation data). Returns 0 if there is no valid navigation data. */ double GetPositionErrorSampleStandardDeviation(); /** @return Returns the RMS of the errors of all positions of the analysed point set (only valid navigation data). Returns 0 if there is no valid navigation data. */ double GetPositionErrorRMS(); /** @return Returns the median of the errors of all positions of the analysed point set (only valid navigation data). Returns 0 if there is no valid navigation data. */ double GetPositionErrorMedian(); /** @return Returns the maximum of the errors of all positions of the analysed point set (only valid navigation data). Returns 0 if there is no valid navigation data. */ double GetPositionErrorMax(); /** @return Returns the minimum of the errors of all positions of the analysed point set (only valid navigation data). Returns 0 if there is no valid navigation data. */ double GetPositionErrorMin(); //##################################################################################################### //this both methods are used by another class an so they are public... perhaps we want to move them //out of this class because they have nothing to do with point sets. /** @return returns the standard derivation of the given list (NOT of the point set).*/ double GetStabw(std::vector list); /** @return returns the sample standard derivation of the given list (NOT of the point set).*/ double GetSampleStabw(std::vector list); //##################################################################################################### protected: PointSetStatisticsCalculator(); - PointSetStatisticsCalculator(mitk::PointSet::Pointer); + explicit PointSetStatisticsCalculator(mitk::PointSet::Pointer); ~PointSetStatisticsCalculator() override; // TODO: Remove the std::vector data structure and use mitk::PointSet everywhere /** @return Returns a list with the distances to the mean of the list */ std::vector GetErrorList(std::vector list); /** @return Returns the mean of the point list. Returns [0;0;0] if the list is empty. */ mitk::Point3D GetMean(std::vector list); /** @brief Converts a point set to a vector of Point3D. */ std::vector PointSetToVector(mitk::PointSet::Pointer pSet); /** @return Returns true if all positions in the evaluated points set are equal. False if not. */ bool CheckIfAllPositionsAreEqual(); mitk::PointSet::Pointer m_PointSet; double GetMean(std::vector list); double GetMedian(std::vector list); double GetMax(std::vector list); double GetMin(std::vector list); }; } #endif // #define _MITK_PointSetSTATISTICSCALCULATOR_H