diff --git a/Modules/DiffusionImaging/MiniApps/ImageStatisticsMiniapp_v2.cpp b/Modules/DiffusionImaging/MiniApps/ImageStatisticsMiniapp_v2.cpp index 9735413843..297b92f63c 100644 --- a/Modules/DiffusionImaging/MiniApps/ImageStatisticsMiniapp_v2.cpp +++ b/Modules/DiffusionImaging/MiniApps/ImageStatisticsMiniapp_v2.cpp @@ -1,210 +1,229 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkCommandLineParser.h" #include "mitkImage.h" #include #include #include #include #include #include "mitkIOUtil.h" #include #include #include #include #include "mitkImageAccessByItk.h" #include #include #include #include #include #include #include #include #include #include #include +#include struct statistics_res{ double mean, variance, min, max, count, moment; }; +long getTimeInMs() +{ + std::chrono::milliseconds ms = std::chrono::duration_cast< std::chrono::milliseconds > (std::chrono::system_clock::now().time_since_epoch()); + long time = ms.count(); + return time; +} + void printstats(statistics_res s) { std::cout << "mean: " << s.mean << std::endl << "variance: " << s.variance << std::endl << "min: " << s.min << std::endl << "max: " << s.max << std::endl << "count: " << s.count << std::endl << "moment: " << s.moment << std::endl; } template void printMap(std::map input) { for (auto it = input.begin(); it != input.end(); ++it) { std::cout << it->first<< ": " << it->second<< std::endl; } std::cout << std::endl; } template < typename TPixel, unsigned int VImageDimension > void get_statistics_boost(itk::Image* itkImage, statistics_res& res){ typedef itk::Image ImageType; itk::ImageRegionConstIterator it(itkImage, itkImage->GetLargestPossibleRegion()); int ctr=0; double sum=0; boost::accumulators::accumulator_set> > acc; for (it.GoToBegin(); !it.IsAtEnd(); ++it) { acc(it.Get()); // currentPixel = it.Get(); // sum+=currentPixel; // ctr+=1; } // res.mean=(double)sum/(double)ctr; res.mean = boost::accumulators::mean(acc); res.variance = boost::accumulators::variance(acc); res.min = boost::accumulators::min(acc); res.max = boost::accumulators::max(acc); res.count = boost::accumulators::count(acc); res.moment = boost::accumulators::moment<2>(acc); std::cout << "sum: " << sum << " N: " << ctr << " mean: " << res.mean << std::endl; } int main( int argc, char* argv[] ) { unsigned int timeStep = 0; std::string inputImageFile; //inputImageFile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic2DplusT.nrrd"; //inputImageFile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic3D.nrrd"; - inputImageFile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/US4DCyl.nrrd"; + inputImageFile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/leber_ct.pic"; + //inputImageFile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/US4DCyl.nrrd"; // Load image mitk::Image::Pointer inputImage = mitk::IOUtil::LoadImage(inputImageFile); mitk::Image::Pointer maskImage; - std::string maskImageFile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic3D_someSegmentation.nrrd"; + //std::string maskImageFile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic3D_someSegmentation.nrrd"; + std::string maskImageFile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/leber_ct_segmentation.nrrd"; maskImage = mitk::IOUtil::LoadImage(maskImageFile); mitk::ImageMaskGenerator::Pointer binaryImageMaskGen = mitk::ImageMaskGenerator::New(); binaryImageMaskGen->SetImageMask(maskImage); mitk::HotspotMaskGenerator::Pointer hotSpotMaskGen = mitk::HotspotMaskGenerator::New(); hotSpotMaskGen->SetImage(inputImage); hotSpotMaskGen->SetHotspotRadiusInMM(3); hotSpotMaskGen->SetTimeStep(0); - std::string pFfile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic2DplusT_ellipse.pf"; + //std::string pFfile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic2DplusT_ellipse.pf"; + std::string pFfile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic3D_rectangle.pf"; std::vector loadedObjects; loadedObjects = mitk::IOUtil::Load(pFfile); mitk::BaseData::Pointer pf = loadedObjects[0]; mitk::PlanarFigure::Pointer planarFigure = dynamic_cast(pf.GetPointer()); mitk::PlanarFigureMaskGenerator::Pointer planarFigMaskExtr = mitk::PlanarFigureMaskGenerator::New(); planarFigMaskExtr->SetPlanarFigure(planarFigure); planarFigMaskExtr->SetImage(inputImage); - // Calculate statistics mitk::ImageStatisticsCalculator_v2::Pointer calculator = mitk::ImageStatisticsCalculator_v2::New(); std::cout << "calculating statistics (unmasked) itk: " << std::endl; mitk::ImageStatisticsCalculator_v2::StatisticsContainer::Pointer result; calculator->SetInputImage(inputImage); calculator->SetNBinsForHistogramStatistics(100); std::cout << std::endl; + // std::cout << "unmasked: " << std::endl; // calculator->SetMask(nullptr); // result = calculator->GetStatistics(timeStep); // result->PrintSelf(); // std::cout << std::endl; // std::cout << "planarfigure: " << std::endl; // calculator->SetMask(planarFigMaskExtr.GetPointer()); // result = calculator->GetStatistics(timeStep, 1); // result->PrintSelf(); // std::cout << std::endl; // std::cout << "ignore pixel value mask: " << std::endl; // mitk::IgnorePixelMaskGenerator::Pointer ignPixVal = mitk::IgnorePixelMaskGenerator::New(); // ignPixVal->SetImage(inputImage); // ignPixVal->SetIgnoredPixelValue(0); // calculator->SetMask(ignPixVal.GetPointer()); // result = calculator->GetStatistics(timeStep, 1); // result->PrintSelf(); // std::cout << std::endl; + long start_time; + start_time = getTimeInMs(); + std::cout << "image mask: " << std::endl; + calculator->SetMask(binaryImageMaskGen.GetPointer()); + result = calculator->GetStatistics(timeStep, 1); + result->PrintSelf(); + std::cout << "new image statistics calculator took: " << getTimeInMs()-start_time << " ms." << std::endl; + std::cout << std::endl; -// std::cout << "image mask: " << std::endl; -// calculator->SetMask(binaryImageMaskGen.GetPointer()); -// result = calculator->GetStatistics(timeStep, 1); -// result->PrintSelf(); - - for (unsigned int i=0; i < inputImage->GetTimeSteps(); i++) - { - std::cout << "no mask: " << std::endl; - std::cout << "timestep: " << i << std::endl; - calculator->SetMask(nullptr); - result = calculator->GetStatistics(i); - result->PrintSelf(); - std::cout << std::endl; - } +// for (unsigned int i=0; i < inputImage->GetTimeSteps(); i++) +// { +// std::cout << "no mask: " << std::endl; +// std::cout << "timestep: " << i << std::endl; +// calculator->SetMask(nullptr); +// result = calculator->GetStatistics(i); +// result->PrintSelf(); +// std::cout << std::endl; +// } // std::cout << "hotspot mask: " << std::endl; // calculator->SetMask(hotSpotMaskGen.GetPointer()); // result = calculator->GetStatistics(timeStep, 1); // result->PrintSelf(); -// mitk::IOUtil::SaveImage(hotSpotMaskGen->GetMask(), "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic2DplusTHotspot.nrrd"); +// //mitk::IOUtil::SaveImage(hotSpotMaskGen->GetMask(), "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic2DplusTHotspot.nrrd"); +// std::cout << std::endl; -// std::cout << std::endl << std::endl << "calculating statistics old imgstatcalc: " << std::endl; + std::cout << std::endl << std::endl << "calculating statistics old imgstatcalc: " << std::endl; -// mitk::ImageStatisticsCalculator::Statistics statisticsStruct; -// mitk::ImageStatisticsCalculator::Pointer calculator_old = mitk::ImageStatisticsCalculator::New(); -// calculator_old->SetImage(inputImage); -// calculator_old->SetImageMask(maskImage); -// //calculator_old->SetMaskingModeToImage(); -// calculator_old->SetHistogramBinSize(100); + start_time = getTimeInMs(); + mitk::ImageStatisticsCalculator::Statistics statisticsStruct; + mitk::ImageStatisticsCalculator::Pointer calculator_old = mitk::ImageStatisticsCalculator::New(); + calculator_old->SetImage(inputImage); + calculator_old->SetImageMask(maskImage); + calculator_old->SetMaskingModeToImage(); + calculator_old->SetHistogramBinSize(100); // calculator_old->SetPlanarFigure(planarFigure); // calculator_old->SetMaskingModeToPlanarFigure(); -// calculator_old->SetMaskingModeToNone(); -// calculator_old->SetDoIgnorePixelValue(false); -// calculator_old->ComputeStatistics(timeStep); -// statisticsStruct = calculator_old->GetStatistics(timeStep); - - - + //calculator_old->SetMaskingModeToNone(); + calculator_old->SetDoIgnorePixelValue(false); + calculator_old->ComputeStatistics(timeStep); + statisticsStruct = calculator_old->GetStatistics(timeStep); + std::cout << "old image statistics calculator took: " << getTimeInMs()-start_time << " ms." << std::endl; + + std::cout << "entropy " << statisticsStruct.GetEntropy() << std::endl; + std::cout << "uniformity " << statisticsStruct.GetUniformity() << std::endl; + std::cout << "upp " << statisticsStruct.GetUPP() << std::endl; + std::cout << "median " << statisticsStruct.GetMedian() << std::endl; // std::cout << "calculating statistics boost: " << std::endl; // statistics_res res; // AccessByItk_n(inputImage, get_statistics_boost, (res)); // printstats(res); return EXIT_SUCCESS; } diff --git a/Modules/ImageStatistics/files.cmake b/Modules/ImageStatistics/files.cmake index c77d7d9ee7..b25b226b74 100644 --- a/Modules/ImageStatistics/files.cmake +++ b/Modules/ImageStatistics/files.cmake @@ -1,35 +1,37 @@ set(CPP_FILES mitkImageStatisticsCalculator.cpp mitkImageStatisticsCalculator2.cpp mitkPointSetStatisticsCalculator.cpp mitkPointSetDifferenceStatisticsCalculator.cpp mitkIntensityProfile.cpp mitkHotspotMaskGenerator.cpp mitkMaskGenerator.cpp mitkPlanarFigureMaskGenerator.cpp mitkMultiLabelMaskGenerator.cpp mitkImageMaskGenerator.cpp mitkHistogramStatisticsCalculator.cpp mitkMaskUtilities.cpp mitkIgnorePixelMaskGenerator.cpp ) set(H_FILES mitkImageStatisticsCalculator.h mitkImageStatisticsCalculator2.h mitkPointSetDifferenceStatisticsCalculator.h mitkPointSetStatisticsCalculator.h mitkExtendedStatisticsImageFilter.h mitkExtendedStatisticsImageFilter_2.h mitkExtendedLabelStatisticsImageFilter.h mitkExtendedLabelStatisticsImageFilter_2.h mitkHotspotMaskGenerator.h mitkMaskGenerator.h mitkPlanarFigureMaskGenerator.h mitkMultiLabelMaskGenerator.h mitkImageMaskGenerator.h mitkHistogramStatisticsCalculator.h mitkMaskUtilities.h mitkitkMaskImageFilter.h mitkIgnorePixelMaskGenerator.h + mitkMinMaxImageFilterWithIndex.h + mitkMinMaxLabelmageFilterWithIndex.h ) diff --git a/Modules/ImageStatistics/mitkExtendedLabelStatisticsImageFilter_2.h b/Modules/ImageStatistics/mitkExtendedLabelStatisticsImageFilter_2.h index ee7b7cd92a..bc2a6432cd 100644 --- a/Modules/ImageStatistics/mitkExtendedLabelStatisticsImageFilter_2.h +++ b/Modules/ImageStatistics/mitkExtendedLabelStatisticsImageFilter_2.h @@ -1,293 +1,353 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef __mitkExtendedLabelStatisticsImageFilter_h2 #define __mitkExtendedLabelStatisticsImageFilter_h2 #include "itkLabelStatisticsImageFilter.h" namespace itk { /** * \class ExtendedLabelStatisticsImageFilter * \brief Extension of the itkLabelStatisticsImageFilter that also calculates the Skewness,Kurtosis,Entropy,Uniformity. * * This class inherits from the itkLabelStatisticsImageFilter and * uses its results for the calculation of six additional coefficients: * the Skewness,Kurtosis,Uniformity,UPP,MPP,Entropy * * As these coefficient are based on the mean and the sigma which are both calculated * by the LabelStatisticsImageFilter, the method AfterThreadedGenerateData() is overwritten * and calls ComputeSkewnessKurtosisAndMPP() and ComputeEntropyUniformityAndUPP after the AfterThreadedGenerateData() * while the second coefficient Method is only called when the the method useHistogram is on!!! * implementation of the superclass is called. * */ template< class TInputImage, class TLabelImage > class ExtendedLabelStatisticsImageFilter2 : public LabelStatisticsImageFilter< TInputImage, TLabelImage > { public: typedef ExtendedLabelStatisticsImageFilter2 Self; typedef LabelStatisticsImageFilter < TInputImage, TLabelImage > Superclass; typedef SmartPointer< Self > Pointer; typedef SmartPointer< const Self > ConstPointer; typedef typename Superclass::LabelPixelType LabelPixelType; typedef typename Superclass::RealType RealType; typedef typename Superclass::PixelType PixelType; typedef typename Superclass::MapIterator MapIterator; typedef itk::Statistics::Histogram HistogramType; itkFactorylessNewMacro( Self ); itkCloneMacro( Self ); itkTypeMacro(ExtendedLabelStatisticsImageFilter2, LabelStatisticsImageFilter); /** \class LabelStatistics * \brief Statistics stored per label * \ingroup ITKImageStatistics */ class LabelStatistics { public: // default constructor LabelStatistics() { // initialized to the default values m_Count = NumericTraits< IdentifierType >::ZeroValue(); m_PositivePixelCount = NumericTraits< IdentifierType >::ZeroValue(); m_Sum = NumericTraits< RealType >::ZeroValue(); m_SumOfPositivePixels = NumericTraits< RealType >::ZeroValue(); m_SumOfSquares = NumericTraits< RealType >::ZeroValue(); m_SumOfCubes = NumericTraits< RealType >::ZeroValue(); m_SumOfQuadruples = NumericTraits< RealType >::ZeroValue(); // Set such that the first pixel encountered can be compared m_Minimum = NumericTraits< RealType >::max(); m_Maximum = NumericTraits< RealType >::NonpositiveMin(); // Default these to zero m_Mean = NumericTraits< RealType >::ZeroValue(); m_Sigma = NumericTraits< RealType >::ZeroValue(); m_Variance = NumericTraits< RealType >::ZeroValue(); m_MPP = NumericTraits< RealType >::ZeroValue(); m_Median = NumericTraits< RealType >::ZeroValue(); m_Uniformity = NumericTraits< RealType >::ZeroValue(); m_UPP = NumericTraits< RealType >::ZeroValue(); m_Entropy = NumericTraits< RealType >::ZeroValue(); m_Skewness = NumericTraits< RealType >::ZeroValue(); m_Kurtosis = NumericTraits< RealType >::ZeroValue(); unsigned int imageDimension = itkGetStaticConstMacro(ImageDimension); m_BoundingBox.resize(imageDimension * 2); for ( unsigned int i = 0; i < imageDimension * 2; i += 2 ) { m_BoundingBox[i] = NumericTraits< IndexValueType >::max(); m_BoundingBox[i + 1] = NumericTraits< IndexValueType >::NonpositiveMin(); } m_Histogram = ITK_NULLPTR; } // constructor with histogram enabled LabelStatistics(int size, RealType lowerBound, RealType upperBound) { // initialized to the default values m_Count = NumericTraits< IdentifierType >::ZeroValue(); m_PositivePixelCount = NumericTraits< IdentifierType >::ZeroValue(); m_Sum = NumericTraits< RealType >::ZeroValue(); m_SumOfPositivePixels = NumericTraits< RealType >::ZeroValue(); m_SumOfSquares = NumericTraits< RealType >::ZeroValue(); m_SumOfCubes = NumericTraits< RealType >::ZeroValue(); m_SumOfQuadruples = NumericTraits< RealType >::ZeroValue(); // Set such that the first pixel encountered can be compared m_Minimum = NumericTraits< RealType >::max(); m_Maximum = NumericTraits< RealType >::NonpositiveMin(); // Default these to zero m_Mean = NumericTraits< RealType >::ZeroValue(); m_Sigma = NumericTraits< RealType >::ZeroValue(); m_Variance = NumericTraits< RealType >::ZeroValue(); m_MPP = NumericTraits< RealType >::ZeroValue(); m_Median = NumericTraits< RealType >::ZeroValue(); m_Uniformity = NumericTraits< RealType >::ZeroValue(); m_UPP = NumericTraits< RealType >::ZeroValue(); m_Entropy = NumericTraits< RealType >::ZeroValue(); m_Skewness = NumericTraits< RealType >::ZeroValue(); m_Kurtosis = NumericTraits< RealType >::ZeroValue(); unsigned int imageDimension = itkGetStaticConstMacro(ImageDimension); m_BoundingBox.resize(imageDimension * 2); for ( unsigned int i = 0; i < imageDimension * 2; i += 2 ) { m_BoundingBox[i] = NumericTraits< IndexValueType >::max(); m_BoundingBox[i + 1] = NumericTraits< IndexValueType >::NonpositiveMin(); } // Histogram m_Histogram = HistogramType::New(); typename HistogramType::SizeType hsize; typename HistogramType::MeasurementVectorType lb; typename HistogramType::MeasurementVectorType ub; hsize.SetSize(1); lb.SetSize(1); ub.SetSize(1); m_Histogram->SetMeasurementVectorSize(1); hsize[0] = size; lb[0] = lowerBound; ub[0] = upperBound; m_Histogram->Initialize(hsize, lb, ub); } // need copy constructor because of smart pointer to histogram LabelStatistics(const LabelStatistics & l) { m_Count = l.m_Count; m_Minimum = l.m_Minimum; m_Maximum = l.m_Maximum; m_Mean = l.m_Mean; m_Sum = l.m_Sum; m_SumOfSquares = l.m_SumOfSquares; m_Sigma = l.m_Sigma; m_Variance = l.m_Variance; m_MPP = l.m_MPP; m_Median = l.m_Median; m_Uniformity = l.m_Uniformity; m_UPP = l.m_UPP; m_Entropy = l.m_Entropy; m_Skewness = l.m_Skewness; m_Kurtosis = l.m_Kurtosis; m_BoundingBox = l.m_BoundingBox; m_Histogram = l.m_Histogram; m_SumOfPositivePixels = l.m_SumOfPositivePixels; m_PositivePixelCount = l.m_PositivePixelCount; m_SumOfCubes = l.m_SumOfCubes; m_SumOfQuadruples = l.m_SumOfQuadruples; } // added for completeness LabelStatistics &operator= (const LabelStatistics& l) { if(this != &l) { m_Count = l.m_Count; m_Minimum = l.m_Minimum; m_Maximum = l.m_Maximum; m_Mean = l.m_Mean; m_Sum = l.m_Sum; m_SumOfSquares = l.m_SumOfSquares; m_Sigma = l.m_Sigma; m_Variance = l.m_Variance; m_MPP = l.m_MPP; m_Median = l.m_Median; m_Uniformity = l.m_Uniformity; m_UPP = l.m_UPP; m_Entropy = l.m_Entropy; m_Skewness = l.m_Skewness; m_Kurtosis = l.m_Kurtosis; m_BoundingBox = l.m_BoundingBox; m_Histogram = l.m_Histogram; m_SumOfPositivePixels = l.m_SumOfPositivePixels; m_PositivePixelCount = l.m_PositivePixelCount; m_SumOfCubes = l.m_SumOfCubes; m_SumOfQuadruples = l.m_SumOfQuadruples; } return *this; } IdentifierType m_Count; RealType m_Minimum; RealType m_Maximum; RealType m_Mean; RealType m_Sum; RealType m_SumOfSquares; RealType m_Sigma; RealType m_Variance; RealType m_MPP; RealType m_Median; RealType m_Uniformity; RealType m_UPP; RealType m_Entropy; RealType m_Skewness; RealType m_Kurtosis; IdentifierType m_PositivePixelCount; RealType m_SumOfPositivePixels; RealType m_SumOfCubes; RealType m_SumOfQuadruples; typename Superclass::BoundingBoxType m_BoundingBox; typename HistogramType::Pointer m_Histogram; }; /** Type of the map used to store data per label */ typedef itksys::hash_map< LabelPixelType, LabelStatistics > MapType; - typedef typename itksys::hash_map< LabelPixelType, LabelStatistics >::const_iterator MapConstIterator; + typedef 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: - ExtendedLabelStatisticsImageFilter2(); + ExtendedLabelStatisticsImageFilter2(): + m_GlobalHistogramParametersSet(false), + m_LabelHistogramParametersSet(false), + m_PreferGlobalHistogramParameters(false) + { + m_NumBins.set_size(1); + } virtual ~ExtendedLabelStatisticsImageFilter2(){} void AfterThreadedGenerateData(); /** Initialize some accumulators before the threads run. */ void BeforeThreadedGenerateData(); /** Multi-thread version GenerateData. */ void ThreadedGenerateData(const typename TInputImage::RegionType & outputRegionForThread, - typename TInputImage::ThreadIdType threadId); + 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_UseHistograms; + 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_2.hxx" #endif #endif diff --git a/Modules/ImageStatistics/mitkExtendedLabelStatisticsImageFilter_2.hxx b/Modules/ImageStatistics/mitkExtendedLabelStatisticsImageFilter_2.hxx index a38ce4d79d..39dbe5285c 100644 --- a/Modules/ImageStatistics/mitkExtendedLabelStatisticsImageFilter_2.hxx +++ b/Modules/ImageStatistics/mitkExtendedLabelStatisticsImageFilter_2.hxx @@ -1,464 +1,765 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef _mitkExtendedLabelStatisticsImageFilter2_hxx #define _mitkExtendedLabelStatisticsImageFilter2_hxx #include "mitkExtendedLabelStatisticsImageFilter_2.h" #include "itkImageRegionConstIteratorWithIndex.h" #include "itkImageRegionConstIterator.h" #include #include #include "mitkNumericConstants.h" #include "mitkLogMacros.h" #include namespace itk { template< class TInputImage, class TLabelImage > bool ExtendedLabelStatisticsImageFilter2< TInputImage, TLabelImage > ::GetMaskingNonEmpty() const { return m_MaskNonEmpty; } + template< typename TInputImage, typename TLabelImage > + void + ExtendedLabelStatisticsImageFilter2< TInputImage, TLabelImage > + ::SetHistogramParameters(const int numBins, RealType lowerBound, RealType upperBound) + { + m_NumBins[0] = numBins; + m_LowerBound = lowerBound; + m_UpperBound = upperBound; + m_GlobalHistogramParametersSet = true; + m_PreferGlobalHistogramParameters = true; + this->Modified(); + } + + template< typename TInputImage, typename TLabelImage > + void + ExtendedLabelStatisticsImageFilter2< TInputImage, TLabelImage > + ::SetHistogramParametersForLabels(std::map numBins, std::map lowerBound, + std::map upperBound) + { + m_LabelMin = lowerBound; + m_LabelMax = upperBound; + m_LabelNBins = numBins; + m_LabelHistogramParametersSet = true; + m_PreferGlobalHistogramParameters = false; + this->Modified(); + } template< class TInputImage, class TLabelImage > typename ExtendedLabelStatisticsImageFilter2< TInputImage, TLabelImage >::RealType ExtendedLabelStatisticsImageFilter2< TInputImage, TLabelImage > ::GetUniformity(LabelPixelType label) const { - MapIterator mapIt; + StatisticsMapConstIterator mapIt; mapIt = m_LabelStatistics.find(label); if ( mapIt == m_LabelStatistics.end() ) { // label does not exist, return a default value return NumericTraits< PixelType >::Zero; } else { return ( *mapIt ).second.m_Uniformity; } } template< class TInputImage, class TLabelImage > typename ExtendedLabelStatisticsImageFilter2< TInputImage, TLabelImage >::RealType ExtendedLabelStatisticsImageFilter2< TInputImage, TLabelImage > ::GetMedian(LabelPixelType label) const { - MapIterator mapIt; + StatisticsMapConstIterator mapIt; mapIt = m_LabelStatistics.find(label); if ( mapIt == m_LabelStatistics.end() ) { // label does not exist, return a default value return NumericTraits< PixelType >::Zero; } else { return ( *mapIt ).second.m_Median; } } template< class TInputImage, class TLabelImage > typename ExtendedLabelStatisticsImageFilter2< TInputImage, TLabelImage >::RealType ExtendedLabelStatisticsImageFilter2< TInputImage, TLabelImage > ::GetEntropy(LabelPixelType label) const { - MapIterator mapIt; + StatisticsMapConstIterator mapIt; mapIt = m_LabelStatistics.find(label); if ( mapIt == m_LabelStatistics.end() ) { // label does not exist, return a default value return NumericTraits< PixelType >::Zero; } else { return ( *mapIt ).second.m_Entropy; } } template< class TInputImage, class TLabelImage > typename ExtendedLabelStatisticsImageFilter2< TInputImage, TLabelImage >::RealType ExtendedLabelStatisticsImageFilter2< TInputImage, TLabelImage > ::GetUPP(LabelPixelType label) const { - MapIterator mapIt; + StatisticsMapConstIterator mapIt; mapIt = m_LabelStatistics.find(label); if ( mapIt == m_LabelStatistics.end() ) { // label does not exist, return a default value return NumericTraits< PixelType >::Zero; } else { return ( *mapIt ).second.m_UPP; } } template< class TInputImage, class TLabelImage > typename ExtendedLabelStatisticsImageFilter2< TInputImage, TLabelImage >::RealType ExtendedLabelStatisticsImageFilter2< TInputImage, TLabelImage > ::GetMPP(LabelPixelType label) const { - MapIterator mapIt; + StatisticsMapConstIterator mapIt; mapIt = m_LabelStatistics.find(label); if ( mapIt == m_LabelStatistics.end() ) { // label does not exist, return a default value return NumericTraits< PixelType >::Zero; } else { return ( *mapIt ).second.m_MPP; } } template< class TInputImage, class TLabelImage > typename ExtendedLabelStatisticsImageFilter2< TInputImage, TLabelImage >::RealType ExtendedLabelStatisticsImageFilter2< TInputImage, TLabelImage > ::GetKurtosis(LabelPixelType label) const { - MapIterator mapIt; + StatisticsMapConstIterator mapIt; mapIt = m_LabelStatistics.find(label); if ( mapIt == m_LabelStatistics.end() ) { // label does not exist, return a default value return NumericTraits< PixelType >::Zero; } else { return ( *mapIt ).second.m_Kurtosis; } } template< class TInputImage, class TLabelImage > typename ExtendedLabelStatisticsImageFilter2< TInputImage, TLabelImage >::RealType ExtendedLabelStatisticsImageFilter2< TInputImage, TLabelImage > ::GetSkewness(LabelPixelType label) const { - MapIterator mapIt; + StatisticsMapConstIterator mapIt; mapIt = m_LabelStatistics.find(label); if ( mapIt == m_LabelStatistics.end() ) { // label does not exist, return a default value return NumericTraits< PixelType >::Zero; } else { return ( *mapIt ).second.m_Skewness; } } + + template< typename TInputImage, typename TLabelImage > + typename ExtendedLabelStatisticsImageFilter2< TInputImage, TLabelImage >::RealType + ExtendedLabelStatisticsImageFilter2< TInputImage, TLabelImage > + ::GetMinimum(LabelPixelType label) const + { + StatisticsMapConstIterator mapIt; + + mapIt = m_LabelStatistics.find(label); + if ( mapIt == m_LabelStatistics.end() ) + { + // label does not exist, return a default value + return NumericTraits< PixelType >::max(); + } + else + { + return ( *mapIt ).second.m_Minimum; + } + } + + template< typename TInputImage, typename TLabelImage > + typename ExtendedLabelStatisticsImageFilter2< TInputImage, TLabelImage >::RealType + ExtendedLabelStatisticsImageFilter2< TInputImage, TLabelImage > + ::GetMaximum(LabelPixelType label) const + { + StatisticsMapConstIterator mapIt; + + mapIt = m_LabelStatistics.find(label); + if ( mapIt == m_LabelStatistics.end() ) + { + // label does not exist, return a default value + return NumericTraits< PixelType >::NonpositiveMin(); + } + else + { + return ( *mapIt ).second.m_Maximum; + } + } + + template< typename TInputImage, typename TLabelImage > + typename ExtendedLabelStatisticsImageFilter2< TInputImage, TLabelImage >::RealType + ExtendedLabelStatisticsImageFilter2< TInputImage, TLabelImage > + ::GetMean(LabelPixelType label) const + { + StatisticsMapConstIterator mapIt; + + mapIt = m_LabelStatistics.find(label); + if ( mapIt == m_LabelStatistics.end() ) + { + // label does not exist, return a default value + return NumericTraits< PixelType >::ZeroValue(); + } + else + { + return ( *mapIt ).second.m_Mean; + } + } + + template< typename TInputImage, typename TLabelImage > + typename ExtendedLabelStatisticsImageFilter2< TInputImage, TLabelImage >::RealType + ExtendedLabelStatisticsImageFilter2< TInputImage, TLabelImage > + ::GetSum(LabelPixelType label) const + { + StatisticsMapConstIterator mapIt; + + mapIt = m_LabelStatistics.find(label); + if ( mapIt == m_LabelStatistics.end() ) + { + // label does not exist, return a default value + return NumericTraits< PixelType >::ZeroValue(); + } + else + { + return ( *mapIt ).second.m_Sum; + } + } + + template< typename TInputImage, typename TLabelImage > + typename ExtendedLabelStatisticsImageFilter2< TInputImage, TLabelImage >::RealType + ExtendedLabelStatisticsImageFilter2< TInputImage, TLabelImage > + ::GetSigma(LabelPixelType label) const + { + StatisticsMapConstIterator mapIt; + + mapIt = m_LabelStatistics.find(label); + if ( mapIt == m_LabelStatistics.end() ) + { + // label does not exist, return a default value + return NumericTraits< PixelType >::ZeroValue(); + } + else + { + return ( *mapIt ).second.m_Sigma; + } + } + + template< typename TInputImage, typename TLabelImage > + typename ExtendedLabelStatisticsImageFilter2< TInputImage, TLabelImage >::RealType + ExtendedLabelStatisticsImageFilter2< TInputImage, TLabelImage > + ::GetVariance(LabelPixelType label) const + { + StatisticsMapConstIterator mapIt; + + mapIt = m_LabelStatistics.find(label); + if ( mapIt == m_LabelStatistics.end() ) + { + // label does not exist, return a default value + return NumericTraits< PixelType >::ZeroValue(); + } + else + { + return ( *mapIt ).second.m_Variance; + } + } + + template< typename TInputImage, typename TLabelImage > + typename ExtendedLabelStatisticsImageFilter2< TInputImage, TLabelImage >::Superclass::BoundingBoxType + ExtendedLabelStatisticsImageFilter2< TInputImage, TLabelImage > + ::GetBoundingBox(LabelPixelType label) const + { + StatisticsMapConstIterator mapIt; + + mapIt = m_LabelStatistics.find(label); + if ( mapIt == m_LabelStatistics.end() ) + { + typename Superclass::BoundingBoxType emptyBox; + // label does not exist, return a default value + return emptyBox; + } + else + { + return ( *mapIt ).second.m_BoundingBox; + } + } + + template< typename TInputImage, typename TLabelImage > + typename ExtendedLabelStatisticsImageFilter2< TInputImage, TLabelImage >::Superclass::RegionType + ExtendedLabelStatisticsImageFilter2< TInputImage, TLabelImage > + ::GetRegion(LabelPixelType label) const + { + StatisticsMapConstIterator mapIt; + + mapIt = m_LabelStatistics.find(label); + + if ( mapIt == m_LabelStatistics.end() ) + { + typename Superclass::RegionType emptyRegion; + // label does not exist, return a default value + return emptyRegion; + } + else + { + typename Superclass::BoundingBoxType bbox = this->GetBoundingBox(label); + typename Superclass::IndexType index; + typename Superclass::SizeType size; + + unsigned int dimension = bbox.size() / 2; + + for ( unsigned int i = 0; i < dimension; i++ ) + { + index[i] = bbox[2 * i]; + size[i] = bbox[2 * i + 1] - bbox[2 * i] + 1; + } + typename Superclass::RegionType region; + region.SetSize(size); + region.SetIndex(index); + + return region; + } + } + + template< typename TInputImage, typename TLabelImage > + typename ExtendedLabelStatisticsImageFilter2< TInputImage, TLabelImage >::MapSizeType + ExtendedLabelStatisticsImageFilter2< TInputImage, TLabelImage > + ::GetCount(LabelPixelType label) const + { + StatisticsMapConstIterator mapIt; + + mapIt = m_LabelStatistics.find(label); + if ( mapIt == m_LabelStatistics.end() ) + { + // label does not exist, return a default value + return 0; + } + else + { + return ( *mapIt ).second.m_Count; + } + } + + + template< typename TInputImage, typename TLabelImage > + typename ExtendedLabelStatisticsImageFilter2< TInputImage, TLabelImage >::HistogramType::Pointer + ExtendedLabelStatisticsImageFilter2< TInputImage, TLabelImage > + ::GetHistogram(LabelPixelType label) const + { + StatisticsMapConstIterator mapIt; + + mapIt = m_LabelStatistics.find(label); + if ( mapIt == m_LabelStatistics.end() ) + { + // label does not exist, return a default value + return ITK_NULLPTR; + } + else + { + // this will be zero if histograms have not been enabled + return ( *mapIt ).second.m_Histogram; + } + } + + + template< typename TInputImage, typename TLabelImage > void ExtendedLabelStatisticsImageFilter2< TInputImage, TLabelImage > ::BeforeThreadedGenerateData() { ThreadIdType numberOfThreads = this->GetNumberOfThreads(); // Resize the thread temporaries m_LabelStatisticsPerThread.resize(numberOfThreads); // Initialize the temporaries for ( ThreadIdType i = 0; i < numberOfThreads; ++i ) { m_LabelStatisticsPerThread[i].clear(); } // Initialize the final map m_LabelStatistics.clear(); } + template< typename TInputImage, typename TLabelImage > + std::list + ExtendedLabelStatisticsImageFilter2< TInputImage, TLabelImage > + ::GetRelevantLabels() const + { + std::list< int> relevantLabels; + for (int i = 0; i < 4096; ++i ) + { + if ( this->HasLabel( i ) ) + { + relevantLabels.push_back( i ); + } + } + return relevantLabels; + } + template< typename TInputImage, typename TLabelImage > void ExtendedLabelStatisticsImageFilter2< TInputImage, TLabelImage > ::ThreadedGenerateData(const typename TInputImage::RegionType & outputRegionForThread, - typename TInputImage::ThreadIdType threadId) + ThreadIdType threadId) { typename HistogramType::IndexType histogramIndex(1); typename HistogramType::MeasurementVectorType histogramMeasurement(1); const SizeValueType size0 = outputRegionForThread.GetSize(0); if( size0 == 0) { return; } ImageLinearConstIteratorWithIndex< TInputImage > it (this->GetInput(), outputRegionForThread); ImageScanlineConstIterator< TLabelImage > labelIt (this->GetLabelInput(), outputRegionForThread); - MapIterator mapIt; + StatisticsMapIterator mapIt; // support progress methods/callbacks const size_t numberOfLinesToProcess = outputRegionForThread.GetNumberOfPixels() / size0; ProgressReporter progress( this, threadId, numberOfLinesToProcess ); + typedef typename MapType::value_type MapValueType; + // do the work while ( !it.IsAtEnd() ) { while ( !it.IsAtEndOfLine() ) { const RealType & value = static_cast< RealType >( it.Get() ); const LabelPixelType & label = labelIt.Get(); // is the label already in this thread? mapIt = m_LabelStatisticsPerThread[threadId].find(label); if ( mapIt == m_LabelStatisticsPerThread[threadId].end() ) { - // create a new statistics object - typedef typename MapType::value_type MapValueType; - if ( m_UseHistograms ) + // if global histogram parameters are set and preferred then use them + if ( m_PreferGlobalHistogramParameters && m_GlobalHistogramParametersSet ) { mapIt = m_LabelStatisticsPerThread[threadId].insert( MapValueType( label, LabelStatistics(m_NumBins[0], m_LowerBound, m_UpperBound) ) ).first; } + // if we have label histogram parameters then use them. If we encounter a label that has no parameters then use global settings if available + else if(!m_PreferGlobalHistogramParameters && m_LabelHistogramParametersSet) + { + typename std::map::iterator lbIt, ubIt; + typename std::map::iterator nbIt; + + lbIt = m_LabelMin.find(label); + ubIt = m_LabelMax.find(label); + nbIt = m_LabelNBins.find(label); + + // if any of the parameters is lacking for the current label but global histogram params are available, use the global parameters + if ((lbIt == m_LabelMin.end() || ubIt == m_LabelMax.end() || nbIt == m_LabelNBins.end()) && m_GlobalHistogramParametersSet) + { + mapIt = m_LabelStatisticsPerThread[threadId].insert( MapValueType( label, + LabelStatistics(m_NumBins[0], m_LowerBound, + m_UpperBound) ) ).first; + } + // if any of the parameters is lacking for the current label and global histogram params are not available, dont use histograms for this label + else if ((lbIt == m_LabelMin.end() || ubIt == m_LabelMax.end() || nbIt == m_LabelNBins.end()) && !m_GlobalHistogramParametersSet) + { + mapIt = m_LabelStatisticsPerThread[threadId].insert( MapValueType( label, + LabelStatistics() ) ).first; + } + // label histogram parameters are available, use them! + else + { + PixelType lowerBound, upperBound; + unsigned int nBins; + lowerBound = (*lbIt).second; + upperBound = (*ubIt).second; + nBins = (*nbIt).second; + mapIt = m_LabelStatisticsPerThread[threadId].insert( MapValueType( label, + LabelStatistics(nBins, lowerBound, upperBound) ) ).first; + } + } + // neither global nor label specific histogram parameters are set -> don't use histograms else { mapIt = m_LabelStatisticsPerThread[threadId].insert( MapValueType( label, LabelStatistics() ) ).first; } } typename MapType::mapped_type &labelStats = ( *mapIt ).second; // update the values for this label and this thread if ( value < labelStats.m_Minimum ) { labelStats.m_Minimum = value; } if ( value > labelStats.m_Maximum ) { labelStats.m_Maximum = value; } // bounding box is min,max pairs for ( unsigned int i = 0; i < ( 2 * TInputImage::ImageDimension ); i += 2 ) { const typename TInputImage::IndexType & index = it.GetIndex(); if ( labelStats.m_BoundingBox[i] > index[i / 2] ) { labelStats.m_BoundingBox[i] = index[i / 2]; } if ( labelStats.m_BoundingBox[i + 1] < index[i / 2] ) { labelStats.m_BoundingBox[i + 1] = index[i / 2]; } } labelStats.m_Sum += value; labelStats.m_SumOfSquares += ( value * value ); labelStats.m_Count++; labelStats.m_SumOfCubes += std::pow(value, 3.); labelStats.m_SumOfQuadruples += std::pow(value, 4.); if (value > 0) { labelStats.m_PositivePixelCount++; - labelStats.m_SumOfPisitivePixels += value; + labelStats.m_SumOfPositivePixels += value; } // if enabled, update the histogram for this label - if ( m_UseHistograms ) + if ( labelStats.m_Histogram.IsNotNull() ) { histogramMeasurement[0] = value; labelStats.m_Histogram->GetIndex(histogramMeasurement, histogramIndex); labelStats.m_Histogram->IncreaseFrequencyOfIndex(histogramIndex, 1); } + else + { + int x = 0; + } ++labelIt; ++it; } labelIt.NextLine(); it.NextLine(); progress.CompletedPixel(); } } template< class TInputImage, class TLabelImage > void ExtendedLabelStatisticsImageFilter2< TInputImage, TLabelImage >:: AfterThreadedGenerateData() { - MapIterator mapIt; - MapConstIterator threadIt; - ThreadIdType i; + StatisticsMapIterator mapIt; + StatisticsMapConstIterator threadIt; + ThreadIdType i; ThreadIdType numberOfThreads = this->GetNumberOfThreads(); // Run through the map for each thread and accumulate the count, // sum, and sumofsquares for ( i = 0; i < numberOfThreads; i++ ) { // iterate over the map for this thread for ( threadIt = m_LabelStatisticsPerThread[i].begin(); threadIt != m_LabelStatisticsPerThread[i].end(); ++threadIt ) { // does this label exist in the cumulative structure yet? mapIt = m_LabelStatistics.find( ( *threadIt ).first ); if ( mapIt == m_LabelStatistics.end() ) { // create a new entry typedef typename MapType::value_type MapValueType; - if ( m_UseHistograms ) + if ( m_GlobalHistogramParametersSet || m_LabelHistogramParametersSet ) { - mapIt = m_LabelStatistics.insert( MapValueType( ( *threadIt ).first, - LabelStatistics(m_NumBins[0], m_LowerBound, - m_UpperBound) ) ).first; +// mapIt = m_LabelStatistics.insert( MapValueType( ( *threadIt ).first, +// LabelStatistics(m_NumBins[0], m_LowerBound, +// m_UpperBound) ) ).first; + mapIt = m_LabelStatistics.insert( MapValueType( *threadIt ) ).first; + continue; } else { mapIt = m_LabelStatistics.insert( MapValueType( ( *threadIt ).first, LabelStatistics() ) ).first; } } typename MapType::mapped_type &labelStats = ( *mapIt ).second; // accumulate the information from this thread labelStats.m_Count += ( *threadIt ).second.m_Count; labelStats.m_Sum += ( *threadIt ).second.m_Sum; labelStats.m_SumOfSquares += ( *threadIt ).second.m_SumOfSquares; labelStats.m_SumOfPositivePixels += ( *threadIt ).second.m_SumOfPositivePixels; labelStats.m_PositivePixelCount += ( *threadIt ).second.m_PositivePixelCount; labelStats.m_SumOfCubes += ( *threadIt ).second.m_SumOfCubes; labelStats.m_SumOfQuadruples += ( *threadIt ).second.m_SumOfQuadruples; if ( labelStats.m_Minimum > ( *threadIt ).second.m_Minimum ) { labelStats.m_Minimum = ( *threadIt ).second.m_Minimum; } if ( labelStats.m_Maximum < ( *threadIt ).second.m_Maximum ) { labelStats.m_Maximum = ( *threadIt ).second.m_Maximum; } //bounding box is min,max pairs int dimension = labelStats.m_BoundingBox.size() / 2; for ( int ii = 0; ii < ( dimension * 2 ); ii += 2 ) { if ( labelStats.m_BoundingBox[ii] > ( *threadIt ).second.m_BoundingBox[ii] ) { labelStats.m_BoundingBox[ii] = ( *threadIt ).second.m_BoundingBox[ii]; } if ( labelStats.m_BoundingBox[ii + 1] < ( *threadIt ).second.m_BoundingBox[ii + 1] ) { labelStats.m_BoundingBox[ii + 1] = ( *threadIt ).second.m_BoundingBox[ii + 1]; } } // if enabled, update the histogram for this label - if ( m_UseHistograms ) + if ( m_GlobalHistogramParametersSet || m_LabelHistogramParametersSet ) { typename HistogramType::IndexType index; index.SetSize(1); - for ( unsigned int bin = 0; bin < m_NumBins[0]; bin++ ) + for ( unsigned int bin = 0; bin < labelStats.m_Histogram->Size(); bin++ ) { index[0] = bin; labelStats.m_Histogram->IncreaseFrequency( bin, ( *threadIt ).second.m_Histogram->GetFrequency(bin) ); } } } // end of thread map iterator loop } // end of thread loop // compute the remainder of the statistics for ( mapIt = m_LabelStatistics.begin(); mapIt != m_LabelStatistics.end(); ++mapIt ) { typename MapType::mapped_type &labelStats = ( *mapIt ).second; // mean labelStats.m_Mean = labelStats.m_Sum / static_cast< RealType >( labelStats.m_Count ); // MPP labelStats.m_MPP = labelStats.m_SumOfPositivePixels / static_cast< RealType >( labelStats.m_PositivePixelCount ); // variance if ( labelStats.m_Count > 0 ) { // unbiased estimate of variance LabelStatistics & ls = mapIt->second; const RealType sumSquared = ls.m_Sum * ls.m_Sum; const RealType count = static_cast< RealType >( ls.m_Count ); ls.m_Variance = ( ls.m_SumOfSquares - sumSquared / count ) / ( count ); RealType secondMoment = ls.m_SumOfSquares / count; RealType thirdMoment = ls.m_SumOfCubes / count; RealType fourthMoment = ls.m_SumOfQuadruples / count; ls.m_Skewness = (thirdMoment - 3. * secondMoment * ls.m_Mean + 2. * std::pow(ls.m_Mean, 3.)) / std::pow(secondMoment - std::pow(ls.m_Mean, 2.), 1.5); // see http://www.boost.org/doc/libs/1_51_0/doc/html/boost/accumulators/impl/skewness_impl.html - ls.m_Kurtosis = (fourthMoment - 4. * thirdMoment * ls.m_Mean + 6. * secondMoment * std::pow(ls.m_Mean, 2.) - 3. * std::pow(ls.m_Mean, 4.)) / std::pow(secondMoment - std::pow(ls.m_Mean, 2.), 1.5) - 3.; // see http://www.boost.org/doc/libs/1_51_0/doc/html/boost/accumulators/impl/kurtosis_impl.html + ls.m_Kurtosis = (fourthMoment - 4. * thirdMoment * ls.m_Mean + 6. * secondMoment * std::pow(ls.m_Mean, 2.) - 3. * std::pow(ls.m_Mean, 4.)) / std::pow(secondMoment - std::pow(ls.m_Mean, 2.), 2.); // see http://www.boost.org/doc/libs/1_51_0/doc/html/boost/accumulators/impl/kurtosis_impl.html, dropped -3 } else { labelStats.m_Variance = NumericTraits< RealType >::ZeroValue(); labelStats.m_Skewness = NumericTraits< RealType >::ZeroValue(); labelStats.m_Kurtosis = NumericTraits< RealType >::ZeroValue(); } // sigma labelStats.m_Sigma = std::sqrt( labelStats.m_Variance ); // histogram statistics - mitk::HistogramStatisticsCalculator histStatCalc; - histStatCalc.SetHistogram(labelStats.m_Histogram); - histStatCalc.CalculateStatistics(); - labelStats.m_Median = histStatCalc.GetMedian(); - labelStats.m_Entropy = histStatCalc.GetEntropy(); - labelStats.m_Uniformity = histStatCalc.GetUniformity(); - labelStats.m_UPP = histStatCalc.GetUPP(); + if (labelStats.m_Histogram.IsNotNull()) + { + mitk::HistogramStatisticsCalculator histStatCalc; + histStatCalc.SetHistogram(labelStats.m_Histogram); + histStatCalc.CalculateStatistics(); + labelStats.m_Median = histStatCalc.GetMedian(); + labelStats.m_Entropy = histStatCalc.GetEntropy(); + labelStats.m_Uniformity = histStatCalc.GetUniformity(); + labelStats.m_UPP = histStatCalc.GetUPP(); + } } { //Now update the cached vector of valid labels. m_ValidLabelValues.resize(0); m_ValidLabelValues.reserve(m_LabelStatistics.size()); for ( mapIt = m_LabelStatistics.begin(); mapIt != m_LabelStatistics.end(); ++mapIt ) { m_ValidLabelValues.push_back(mapIt->first); } } } } // end namespace itk #endif diff --git a/Modules/ImageStatistics/mitkExtendedStatisticsImageFilter_2.h b/Modules/ImageStatistics/mitkExtendedStatisticsImageFilter_2.h index 50382ab31d..ab829ba0d4 100644 --- a/Modules/ImageStatistics/mitkExtendedStatisticsImageFilter_2.h +++ b/Modules/ImageStatistics/mitkExtendedStatisticsImageFilter_2.h @@ -1,216 +1,216 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef __mitkExtendedStatisticsImageFilter_h2 #define __mitkExtendedStatisticsImageFilter_h2 #include "itkStatisticsImageFilter.h" #include #include namespace itk { /** * \class ExtendedStatisticsImageFilter * \brief Extension of the itkStatisticsImageFilter that also calculates the Skewness and Kurtosis. * * This class inherits from the itkStatisticsImageFilter and * uses its results for the calculation of other additional coefficients: * the Skewness and Kurtosis. * * As these coefficient are based on the mean and the sigma which are both calculated * by the StatisticsImageFilter, the method AfterThreadedGenerateData() is overwritten * and calls ComputeSkewnessKurtosisAndMPP() and ComputeEntropyUniformityAndUPP * after the AfterThreadedGenerateData() * implementation of the superclass is called. * * As the StatisticsImageFilter stores the statistics in the outputs 1 to 6 by the * StatisticsImageFilter, the skewness, kurtosis,MPP,UPP,Uniformity and Entropy are stored in the outputs * 7 to 13 by this filter. */ template< class TInputImage > class ExtendedStatisticsImageFilter2 : public StatisticsImageFilter< TInputImage > { public: /** Standard Self typedef */ typedef ExtendedStatisticsImageFilter2 Self; typedef StatisticsImageFilter< TInputImage > Superclass; typedef SmartPointer< Self > Pointer; typedef SmartPointer< const Self > ConstPointer; typedef typename Superclass::RealType RealType; typedef typename Superclass::RealObjectType RealObjectType; typedef typename Superclass::PixelType PixelType; /** Histogram-related typedefs */ typedef itk::Statistics::Histogram< RealType > HistogramType; typedef typename HistogramType::Pointer HistogramPointer; itkFactorylessNewMacro( Self ); itkCloneMacro( Self ); itkTypeMacro( ExtendedStatisticsImageFilter2, StatisticsImageFilter ); /** * \brief Return the computed Skewness. */ double GetSkewness() const { return this->GetSkewnessOutput()->Get(); } /** * \brief Return the computed Median */ double GetMedian() const { return this->GetMedianOutput()->Get(); } /** * \brief Return the computed Kurtosis. */ double GetKurtosis() const { return this->GetKurtosisOutput()->Get(); } /* \brief Return the computed MPP. */ double GetMPP() const { return this->GetMPPOutput()->Get(); } /** * \brief Return the computed Uniformity. */ double GetUniformity() const { return this->GetUniformityOutput()->Get(); } /** *\brief Return the computed Entropy. */ double GetEntropy() const { return this->GetEntropyOutput()->Get(); } /** * \brief Return the computed UPP. */ double GetUPP() const { return this->GetUPPOutput()->Get(); } /** * \brief Return the computed Histogram. */ const typename HistogramType::Pointer GetHistogram() { if (m_HistogramCalculated) { return m_Histogram; } else { return ITK_NULLPTR; } } /** specify Histogram parameters */ void SetHistogramParameters(const int numBins, RealType lowerBound, RealType upperBound); protected: ExtendedStatisticsImageFilter2(); virtual ~ExtendedStatisticsImageFilter2(){}; void BeforeThreadedGenerateData(); /** Multi-thread version GenerateData. */ void ThreadedGenerateData(const typename StatisticsImageFilter::RegionType & outputRegionForThread, - typename StatisticsImageFilter::ThreadIdType threadId); + ThreadIdType threadId); /** * brief Calls AfterThreadedGenerateData() of the superclass and the main methods */ void AfterThreadedGenerateData(); RealObjectType* GetSkewnessOutput(); const RealObjectType* GetSkewnessOutput() const; RealObjectType* GetKurtosisOutput(); const RealObjectType* GetKurtosisOutput() const; RealObjectType* GetMPPOutput(); const RealObjectType* GetMPPOutput() const; RealObjectType* GetEntropyOutput(); const RealObjectType* GetEntropyOutput() const; RealObjectType* GetUniformityOutput(); const RealObjectType* GetUniformityOutput() const; RealObjectType* GetUPPOutput(); const RealObjectType* GetUPPOutput() const; RealObjectType* GetMedianOutput(); const RealObjectType* GetMedianOutput() const; virtual DataObject::Pointer MakeOutput( ProcessObject::DataObjectPointerArraySizeType idx ); private: Array< RealType > m_ThreadSum; Array< RealType > m_SumOfSquares; Array< RealType > m_SumOfCubes; Array< RealType > m_SumOfQuadruples; Array< SizeValueType > m_Count; Array< SizeValueType > m_PositivePixelCount; Array< RealType > m_ThreadSumOfPositivePixels; Array< PixelType > m_ThreadMin; Array< PixelType > m_ThreadMax; - Array< HistogramPointer > m_HistogramPerThread; + std::vector< HistogramPointer > m_HistogramPerThread; HistogramPointer m_Histogram; bool m_UseHistogram; bool m_HistogramCalculated; RealType m_LowerBound, m_UpperBound; - typename HistogramType::SizeType m_NumBins; + int m_NumBins; }; // end of class } // end namespace itk #ifndef ITK_MANUAL_INSTANTIATION #include "mitkExtendedStatisticsImageFilter_2.hxx" #endif #endif diff --git a/Modules/ImageStatistics/mitkExtendedStatisticsImageFilter_2.hxx b/Modules/ImageStatistics/mitkExtendedStatisticsImageFilter_2.hxx index b669e93de5..3b9db1fc08 100644 --- a/Modules/ImageStatistics/mitkExtendedStatisticsImageFilter_2.hxx +++ b/Modules/ImageStatistics/mitkExtendedStatisticsImageFilter_2.hxx @@ -1,475 +1,477 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef __mitkExtendedStatisticsImageFilter2_hxx #define __mitkExtendedStatisticsImageFilter2_hxx #include "mitkExtendedStatisticsImageFilter_2.h" #include namespace itk { template< class TInputImage > ExtendedStatisticsImageFilter2< TInputImage >::ExtendedStatisticsImageFilter2() : StatisticsImageFilter< TInputImage >(), m_ThreadSum(1), m_SumOfSquares(1), m_SumOfCubes(1), m_SumOfQuadruples(1), m_Count(1), m_ThreadMin(1), m_ThreadMax(1), m_PositivePixelCount(1), m_ThreadSumOfPositivePixels(1), - m_HistogramPerThread(1), m_UseHistogram(false), m_HistogramCalculated(false) { /* * add the Skewness,Kurtosis,Entropy,Uniformity,MPP, UPP, Median to the other statistical calculated Values * of the mitkStatisticsImageFilter as the 7th to the 13th Output */ for ( int i = 7; i < 14; ++i ) { typename RealObjectType::Pointer output = static_cast< RealObjectType * >( this->MakeOutput(i).GetPointer() ); this->ProcessObject::SetNthOutput( i, output.GetPointer() ); } this->GetSkewnessOutput()->Set( 0.0 ); this->GetKurtosisOutput()->Set( 0.0 ); this->GetEntropyOutput()->Set( -1.0 ); this->GetUniformityOutput()->Set( 0.0 ); this->GetUPPOutput()->Set( 0.0 ); this->GetMPPOutput()->Set( 0.0 ); this->GetMedianOutput()->Set( 0.0 ); + + m_Histogram = ITK_NULLPTR; + m_HistogramPerThread.resize(0); } template< class TInputImage > DataObject::Pointer ExtendedStatisticsImageFilter2< TInputImage >::MakeOutput( ProcessObject::DataObjectPointerArraySizeType output) { switch ( output ) { case 7: case 8: case 9: case 10: case 11: case 12: { return RealObjectType::New().GetPointer(); break; } case 13: { return RealObjectType::New().GetPointer(); break; } default: { // might as well make an image return Superclass::MakeOutput( output ); break; } } } template< class TInputImage > typename ExtendedStatisticsImageFilter2< TInputImage >::RealObjectType * ExtendedStatisticsImageFilter2< TInputImage > ::GetSkewnessOutput() { return static_cast< RealObjectType * >( this->ProcessObject::GetOutput(7) ); } template< class TInputImage > const typename ExtendedStatisticsImageFilter2< TInputImage >::RealObjectType * ExtendedStatisticsImageFilter2< TInputImage > ::GetSkewnessOutput() const { return static_cast< const RealObjectType * >( this->ProcessObject::GetOutput(7) ); } template< class TInputImage > typename ExtendedStatisticsImageFilter2< TInputImage >::RealObjectType * ExtendedStatisticsImageFilter2< TInputImage > ::GetKurtosisOutput() { return static_cast< RealObjectType * >( this->ProcessObject::GetOutput(8) ); } template< class TInputImage > const typename ExtendedStatisticsImageFilter2< TInputImage >::RealObjectType * ExtendedStatisticsImageFilter2< TInputImage > ::GetKurtosisOutput() const { return static_cast< const RealObjectType * >( this->ProcessObject::GetOutput(8) ); } template< class TInputImage > typename ExtendedStatisticsImageFilter2< TInputImage >::RealObjectType * ExtendedStatisticsImageFilter2< TInputImage > ::GetUniformityOutput() { return static_cast< RealObjectType * >( this->ProcessObject::GetOutput(9) ); } template< class TInputImage > const typename ExtendedStatisticsImageFilter2< TInputImage >::RealObjectType * ExtendedStatisticsImageFilter2< TInputImage > ::GetUniformityOutput() const { return static_cast< const RealObjectType * >( this->ProcessObject::GetOutput(9) ); } template< class TInputImage > typename ExtendedStatisticsImageFilter2< TInputImage >::RealObjectType * ExtendedStatisticsImageFilter2< TInputImage > ::GetEntropyOutput() { return static_cast< RealObjectType * >( this->ProcessObject::GetOutput(10) ); } template< class TInputImage > const typename ExtendedStatisticsImageFilter2< TInputImage >::RealObjectType * ExtendedStatisticsImageFilter2< TInputImage > ::GetEntropyOutput() const { return static_cast< const RealObjectType * >( this->ProcessObject::GetOutput(10) ); } template< class TInputImage > typename ExtendedStatisticsImageFilter2< TInputImage >::RealObjectType * ExtendedStatisticsImageFilter2< TInputImage > ::GetUPPOutput() { return static_cast< RealObjectType * >( this->ProcessObject::GetOutput(11) ); } template< class TInputImage > const typename ExtendedStatisticsImageFilter2< TInputImage >::RealObjectType * ExtendedStatisticsImageFilter2< TInputImage > ::GetUPPOutput() const { return static_cast< const RealObjectType * >( this->ProcessObject::GetOutput(11) ); } template< class TInputImage > typename ExtendedStatisticsImageFilter2< TInputImage >::RealObjectType * ExtendedStatisticsImageFilter2< TInputImage > ::GetMPPOutput() { return static_cast< RealObjectType * >( this->ProcessObject::GetOutput(12) ); } template< class TInputImage > const typename ExtendedStatisticsImageFilter2< TInputImage >::RealObjectType * ExtendedStatisticsImageFilter2< TInputImage > ::GetMPPOutput() const { return static_cast< const RealObjectType * >( this->ProcessObject::GetOutput(12) ); } template< class TInputImage > typename ExtendedStatisticsImageFilter2< TInputImage >::RealObjectType * ExtendedStatisticsImageFilter2< TInputImage > ::GetMedianOutput() { return static_cast< RealObjectType * >( this->ProcessObject::GetOutput(13) ); } template< class TInputImage > const typename ExtendedStatisticsImageFilter2< TInputImage >::RealObjectType * ExtendedStatisticsImageFilter2< TInputImage > ::GetMedianOutput() const { return static_cast< const RealObjectType * >( this->ProcessObject::GetOutput(13) ); } template< typename TInputImage > void ExtendedStatisticsImageFilter2< TInputImage > ::BeforeThreadedGenerateData() { ThreadIdType numberOfThreads = this->GetNumberOfThreads(); if (m_UseHistogram) { // Histogram m_Histogram = HistogramType::New(); typename HistogramType::SizeType hsize; typename HistogramType::MeasurementVectorType lb; typename HistogramType::MeasurementVectorType ub; hsize.SetSize(1); lb.SetSize(1); ub.SetSize(1); m_Histogram->SetMeasurementVectorSize(1); hsize[0] = m_NumBins; lb[0] = m_LowerBound; ub[0] = m_UpperBound; m_Histogram->Initialize(hsize, lb, ub); - m_HistogramPerThread.set_size(numberOfThreads); + 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); - m_Histogram->SetMeasurementVectorSize(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 ExtendedStatisticsImageFilter2< TInputImage > ::ThreadedGenerateData(const typename StatisticsImageFilter::RegionType & outputRegionForThread, - typename StatisticsImageFilter::ThreadIdType threadId) + ThreadIdType threadId) { const SizeValueType size0 = outputRegionForThread.GetSize(0); if( size0 == 0) { return; } RealType realValue; PixelType value; typename HistogramType::IndexType histogramIndex(1); typename HistogramType::MeasurementVectorType histogramMeasurement(1); RealType sum = NumericTraits< RealType >::ZeroValue(); RealType sumOfPositivePixels = NumericTraits< RealType >::ZeroValue(); RealType sumOfSquares = NumericTraits< RealType >::ZeroValue(); RealType sumOfCubes = NumericTraits< RealType >::ZeroValue(); RealType sumOfQuadruples = NumericTraits< RealType >::ZeroValue(); SizeValueType count = NumericTraits< SizeValueType >::ZeroValue(); SizeValueType countOfPositivePixels = NumericTraits< SizeValueType >::ZeroValue(); PixelType min = NumericTraits< PixelType >::max(); PixelType max = NumericTraits< PixelType >::NonpositiveMin(); ImageScanlineConstIterator< TInputImage > it (this->GetInput(), outputRegionForThread); // support progress methods/callbacks const size_t numberOfLinesToProcess = outputRegionForThread.GetNumberOfPixels() / size0; ProgressReporter progress( this, threadId, numberOfLinesToProcess ); // do the work while ( !it.IsAtEnd() ) { while ( !it.IsAtEndOfLine() ) { value = it.Get(); realValue = static_cast< RealType >( value ); if (m_UseHistogram) { histogramMeasurement[0] = value; m_HistogramPerThread[threadId]->GetIndex(histogramMeasurement, histogramIndex); m_HistogramPerThread[threadId]->IncreaseFrequencyOfIndex(histogramIndex, 1); } if ( value < min ) { min = value; } if ( value > max ) { max = value; } if (value > 0) { sumOfPositivePixels += realValue; ++countOfPositivePixels; } sum += realValue; sumOfSquares += ( realValue * realValue ); sumOfCubes += std::pow(realValue, 3.); sumOfQuadruples += std::pow(realValue, 4.); ++count; ++it; } it.NextLine(); progress.CompletedPixel(); } m_ThreadSum[threadId] = sum; m_SumOfSquares[threadId] = sumOfSquares; m_Count[threadId] = count; m_ThreadMin[threadId] = min; m_ThreadMax[threadId] = max; m_ThreadSumOfPositivePixels[threadId] = sumOfPositivePixels; m_PositivePixelCount[threadId] = countOfPositivePixels; m_SumOfCubes[threadId] = sumOfCubes; m_SumOfQuadruples[threadId] = sumOfQuadruples; } template< class TInputImage > void ExtendedStatisticsImageFilter2< TInputImage > ::AfterThreadedGenerateData() { ThreadIdType i; SizeValueType count; SizeValueType countOfPositivePixels; RealType sumOfSquares; RealType sumOfCubes; RealType sumOfQuadruples; RealType sumOfPositivePixels; ThreadIdType numberOfThreads = this->GetNumberOfThreads(); PixelType minimum; PixelType maximum; RealType mean; RealType meanOfPositivePixels; RealType sigma; RealType variance; RealType skewness; RealType kurtosis; RealType sum; sum = sumOfSquares = sumOfCubes = sumOfQuadruples = NumericTraits< RealType >::ZeroValue(); count = countOfPositivePixels = 0; // Find the min/max over all threads and accumulate count, sum and // sum of squares minimum = NumericTraits< PixelType >::max(); maximum = NumericTraits< PixelType >::NonpositiveMin(); for ( i = 0; i < numberOfThreads; i++ ) { count += m_Count[i]; sum += m_ThreadSum[i]; sumOfSquares += m_SumOfSquares[i]; sumOfCubes += m_SumOfCubes[i]; sumOfQuadruples += m_SumOfQuadruples[i]; sumOfPositivePixels += m_ThreadSumOfPositivePixels[i]; countOfPositivePixels += m_PositivePixelCount[i]; if ( m_ThreadMin[i] < minimum ) { minimum = m_ThreadMin[i]; } if ( m_ThreadMax[i] > maximum ) { maximum = m_ThreadMax[i]; } // if enabled, update the histogram for this label if ( m_UseHistogram ) { typename HistogramType::IndexType index; index.SetSize(1); - for ( unsigned int bin = 0; bin < m_NumBins[0]; bin++ ) + for ( unsigned 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.) - 3; // see http://www.boost.org/doc/libs/1_46_1/doc/html/boost/accumulators/impl/kurtosis_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 ExtendedStatisticsImageFilter2< TInputImage > ::SetHistogramParameters(const int numBins, RealType lowerBound, RealType upperBound) { - m_NumBins[0] = numBins; + m_NumBins = numBins; m_LowerBound = lowerBound; m_UpperBound = upperBound; m_UseHistogram = true; } } #endif diff --git a/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp b/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp index 03128441fc..612aeb71f6 100644 --- a/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp +++ b/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp @@ -1,2246 +1,2246 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkImageStatisticsCalculator.h" #include "mitkImageAccessByItk.h" #include "mitkImageCast.h" #include "mitkExtractImageFilter.h" #include "mitkImageTimeSelector.h" #include "mitkITKImageImport.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include "itkImage.h" //#define DEBUG_HOTSPOTSEARCH #define _USE_MATH_DEFINES #include #include "vtkLassoStencilSource.h" namespace mitk { ImageStatisticsCalculator::ImageStatisticsCalculator() : m_MaskingMode( MASKING_MODE_NONE ), m_MaskingModeChanged( false ), m_IgnorePixelValue(0.0), m_DoIgnorePixelValue(false), m_IgnorePixelValueChanged(false), m_PlanarFigureAxis (0), m_PlanarFigureSlice (0), m_PlanarFigureCoordinate0 (0), m_PlanarFigureCoordinate1 (0), m_HistogramBinSize(1.0), m_UseDefaultBinSize(true), m_UseBinSizeBasedOnVOIRegion(false), m_HotspotRadiusInMM(6.2035049089940), // radius of a 1cm3 sphere in mm m_CalculateHotspot(false), m_HotspotRadiusInMMChanged(false), m_HotspotMustBeCompletelyInsideImage(true) { m_EmptyHistogram = HistogramType::New(); m_EmptyHistogram->SetMeasurementVectorSize(1); HistogramType::SizeType histogramSize(1); histogramSize.Fill( 256 ); m_EmptyHistogram->Initialize( histogramSize ); m_EmptyStatistics.Reset(); } ImageStatisticsCalculator::~ImageStatisticsCalculator() { } void ImageStatisticsCalculator::SetUseDefaultBinSize(bool useDefault) { m_UseDefaultBinSize = useDefault; } ImageStatisticsCalculator::Statistics::Statistics(bool withHotspotStatistics) :m_HotspotStatistics(withHotspotStatistics ? new Statistics(false) : nullptr) { Reset(); } ImageStatisticsCalculator::Statistics::Statistics(const Statistics& other) :m_HotspotStatistics( nullptr) { this->SetLabel( other.GetLabel() ); this->SetN( other.GetN() ); this->SetMin( other.GetMin() ); this->SetMax( other.GetMax() ); this->SetMedian( other.GetMedian() ); this->SetMean( other.GetMean() ); this->SetVariance( other.GetVariance() ); this->SetKurtosis( other.GetKurtosis() ); this->SetSkewness( other.GetSkewness() ); this->SetUniformity( other.GetUniformity() ); this->SetEntropy( other.GetEntropy() ); this->SetUPP( other.GetUPP() ); this->SetMPP( other.GetMPP() ); this->SetSigma( other.GetSigma() ); this->SetRMS( other.GetRMS() ); this->SetMaxIndex( other.GetMaxIndex() ); this->SetMinIndex( other.GetMinIndex() ); this->SetHotspotIndex( other.GetHotspotIndex() ); if (other.m_HotspotStatistics) { this->m_HotspotStatistics = new Statistics(false); *this->m_HotspotStatistics = *other.m_HotspotStatistics; } } bool ImageStatisticsCalculator::Statistics::HasHotspotStatistics() const { return m_HotspotStatistics != nullptr; } void ImageStatisticsCalculator::Statistics::SetHasHotspotStatistics(bool hasHotspotStatistics) { m_HasHotspotStatistics = hasHotspotStatistics; } ImageStatisticsCalculator::Statistics::~Statistics() { delete m_HotspotStatistics; } double ImageStatisticsCalculator::Statistics::GetVariance() const { return this->Variance; } void ImageStatisticsCalculator::Statistics::SetVariance( const double value ) { if( this->Variance != value ) { if( value < 0.0 ) { this->Variance = 0.0; // if given value is negative set variance to 0.0 } else { this->Variance = value; } } } double ImageStatisticsCalculator::Statistics::GetSigma() const { return this->Sigma; } void ImageStatisticsCalculator::Statistics::SetSigma( const double value ) { if( this->Sigma != value ) { // for some compiler the value != value works to check for NaN but not for all // but we can always be sure that the standard deviation is a positive value if( value != value || value < 0.0 ) { // if standard deviation is NaN we just assume 0.0 this->Sigma = 0.0; } else { this->Sigma = value; } } } void ImageStatisticsCalculator::Statistics::Reset(unsigned int dimension) { SetLabel(0); SetN( 0 ); SetMin( 0.0 ); SetMax( 0.0 ); SetMedian( 0.0 ); SetVariance( 0.0 ); SetMean( 0.0 ); SetSigma( 0.0 ); SetRMS( 0.0 ); SetSkewness( 0.0 ); SetKurtosis( 0.0 ); SetUniformity( 0.0 ); SetEntropy( 0.0 ); SetMPP( 0.0 ); SetUPP( 0.0 ); vnl_vector zero; zero.set_size(dimension); for(unsigned int i = 0; i < dimension; ++i) { zero[i] = 0; } SetMaxIndex(zero); SetMinIndex(zero); SetHotspotIndex(zero); if (m_HotspotStatistics != nullptr) { m_HotspotStatistics->Reset(dimension); } } const ImageStatisticsCalculator::Statistics& ImageStatisticsCalculator::Statistics::GetHotspotStatistics() const { if (m_HotspotStatistics) { return *m_HotspotStatistics; } else { throw std::logic_error("Object has no hostspot statistics, see HasHotspotStatistics()"); } } ImageStatisticsCalculator::Statistics& ImageStatisticsCalculator::Statistics::GetHotspotStatistics() { if (m_HotspotStatistics) { return *m_HotspotStatistics; } else { throw std::logic_error("Object has no hostspot statistics, see HasHotspotStatistics()"); } } ImageStatisticsCalculator::Statistics& ImageStatisticsCalculator::Statistics::operator=(ImageStatisticsCalculator::Statistics const& other) { if (this == &other) return *this; this->SetLabel( other.GetLabel() ); this->SetN( other.GetN() ); this->SetMin( other.GetMin() ); this->SetMax( other.GetMax() ); this->SetMean( other.GetMean() ); this->SetMedian( other.GetMedian() ); this->SetVariance( other.GetVariance() ); this->SetSigma( other.GetSigma() ); this->SetRMS( other.GetRMS() ); this->SetMinIndex( other.GetMinIndex() ); this->SetMaxIndex( other.GetMaxIndex() ); this->SetHotspotIndex( other.GetHotspotIndex() ); this->SetSkewness( other.GetSkewness() ); this->SetKurtosis( other.GetKurtosis() ); this->SetUniformity( other.GetUniformity() ); this->SetEntropy( other.GetEntropy() ); this->SetUPP( other.GetUPP() ); this->SetMPP( other.GetMPP() ); delete this->m_HotspotStatistics; this->m_HotspotStatistics = nullptr; if (other.m_HotspotStatistics) { this->m_HotspotStatistics = new Statistics(false); *this->m_HotspotStatistics = *other.m_HotspotStatistics; } return *this; } void ImageStatisticsCalculator::SetImage( const mitk::Image *image ) { if ( m_Image != image ) { m_Image = image; this->Modified(); unsigned int numberOfTimeSteps = image->GetTimeSteps(); // Initialize vectors to time-size of this image m_ImageHistogramVector.resize( numberOfTimeSteps ); m_MaskedImageHistogramVector.resize( numberOfTimeSteps ); m_PlanarFigureHistogramVector.resize( numberOfTimeSteps ); m_ImageStatisticsVector.resize( numberOfTimeSteps ); m_MaskedImageStatisticsVector.resize( numberOfTimeSteps ); m_PlanarFigureStatisticsVector.resize( numberOfTimeSteps ); m_ImageStatisticsTimeStampVector.resize( numberOfTimeSteps ); m_MaskedImageStatisticsTimeStampVector.resize( numberOfTimeSteps ); m_PlanarFigureStatisticsTimeStampVector.resize( numberOfTimeSteps ); m_ImageStatisticsCalculationTriggerVector.resize( numberOfTimeSteps ); m_MaskedImageStatisticsCalculationTriggerVector.resize( numberOfTimeSteps ); m_PlanarFigureStatisticsCalculationTriggerVector.resize( numberOfTimeSteps ); for ( unsigned int t = 0; t < image->GetTimeSteps(); ++t ) { m_ImageStatisticsTimeStampVector[t].Modified(); m_ImageStatisticsCalculationTriggerVector[t] = true; } } } void ImageStatisticsCalculator::SetImageMask( const mitk::Image *imageMask ) { if ( m_Image.IsNull() ) { itkExceptionMacro( << "Image needs to be set first!" ); } if ( m_ImageMask != imageMask ) { m_ImageMask = imageMask; this->Modified(); for ( unsigned int t = 0; t < m_Image->GetTimeSteps(); ++t ) { m_MaskedImageStatisticsTimeStampVector[t].Modified(); m_MaskedImageStatisticsCalculationTriggerVector[t] = true; } } } void ImageStatisticsCalculator::SetPlanarFigure( mitk::PlanarFigure *planarFigure ) { if ( m_Image.IsNull() ) { itkExceptionMacro( << "Image needs to be set first!" ); } if ( m_PlanarFigure != planarFigure ) { m_PlanarFigure = planarFigure; this->Modified(); for ( unsigned int t = 0; t < m_Image->GetTimeSteps(); ++t ) { m_PlanarFigureStatisticsTimeStampVector[t].Modified(); m_PlanarFigureStatisticsCalculationTriggerVector[t] = true; } } } void ImageStatisticsCalculator::SetMaskingMode( unsigned int mode ) { if ( m_MaskingMode != mode ) { m_MaskingMode = mode; m_MaskingModeChanged = true; this->Modified(); } } void ImageStatisticsCalculator::SetMaskingModeToNone() { if ( m_MaskingMode != MASKING_MODE_NONE ) { m_MaskingMode = MASKING_MODE_NONE; m_MaskingModeChanged = true; this->Modified(); } } void ImageStatisticsCalculator::SetMaskingModeToImage() { if ( m_MaskingMode != MASKING_MODE_IMAGE ) { m_MaskingMode = MASKING_MODE_IMAGE; m_MaskingModeChanged = true; this->Modified(); } } void ImageStatisticsCalculator::SetMaskingModeToPlanarFigure() { if ( m_MaskingMode != MASKING_MODE_PLANARFIGURE ) { m_MaskingMode = MASKING_MODE_PLANARFIGURE; m_MaskingModeChanged = true; this->Modified(); } } void ImageStatisticsCalculator::SetIgnorePixelValue(double value) { if ( m_IgnorePixelValue != value ) { m_IgnorePixelValue = value; if(m_DoIgnorePixelValue) { m_IgnorePixelValueChanged = true; } this->Modified(); } } double ImageStatisticsCalculator::GetIgnorePixelValue() { return m_IgnorePixelValue; } void ImageStatisticsCalculator::SetDoIgnorePixelValue(bool value) { if ( m_DoIgnorePixelValue != value ) { m_DoIgnorePixelValue = value; m_IgnorePixelValueChanged = true; this->Modified(); } } bool ImageStatisticsCalculator::GetDoIgnorePixelValue() { return m_DoIgnorePixelValue; } void ImageStatisticsCalculator::SetHistogramBinSize(double size) { this->m_HistogramBinSize = size; } double ImageStatisticsCalculator::GetHistogramBinSize() { return this->m_HistogramBinSize; } void ImageStatisticsCalculator::SetHotspotRadiusInMM(double value) { if ( m_HotspotRadiusInMM != value ) { m_HotspotRadiusInMM = value; if(m_CalculateHotspot) { m_HotspotRadiusInMMChanged = true; //MITK_INFO <<"Hotspot radius changed, new convolution required"; } this->Modified(); } } double ImageStatisticsCalculator::GetHotspotRadiusInMM() { return m_HotspotRadiusInMM; } void ImageStatisticsCalculator::SetCalculateHotspot(bool on) { if ( m_CalculateHotspot != on ) { m_CalculateHotspot = on; m_HotspotRadiusInMMChanged = true; //MITK_INFO <<"Hotspot calculation changed, new convolution required"; this->Modified(); } } bool ImageStatisticsCalculator::IsHotspotCalculated() { return m_CalculateHotspot; } void ImageStatisticsCalculator::SetHotspotMustBeCompletlyInsideImage(bool hotspotMustBeCompletelyInsideImage, bool warn) { m_HotspotMustBeCompletelyInsideImage = hotspotMustBeCompletelyInsideImage; if (!m_HotspotMustBeCompletelyInsideImage && warn) { MITK_WARN << "Hotspot calculation will extrapolate pixels at image borders. Be aware of the consequences for the hotspot location."; } } bool ImageStatisticsCalculator::GetHotspotMustBeCompletlyInsideImage() const { return m_HotspotMustBeCompletelyInsideImage; } /* Implementation of the min max values for setting the range of the histogram */ template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator::GetMinAndMaxValue( double &min, double &max, int &counter, double &sigma, const itk::Image< TPixel, VImageDimension > *InputImage, itk::Image< unsigned short, VImageDimension > *MaskedImage ) { typedef itk::Image< unsigned short, VImageDimension > MaskImageType; typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::ImageRegionConstIteratorWithIndex Imageie; typedef itk::ImageRegionConstIteratorWithIndex Imageie2; Imageie2 labelIterator2( MaskedImage, MaskedImage->GetRequestedRegion() ); Imageie labelIterator3( InputImage, InputImage->GetRequestedRegion() ); max = 0; min = 0; counter = 0; sigma = 0; double SumOfSquares = 0; double sumSquared = 0; double actualPielValue = 0; int counterOfPixelsInROI = 0; for( labelIterator2.GoToBegin(); !labelIterator2.IsAtEnd(); ++labelIterator2, ++labelIterator3) { if( labelIterator2.Value()== 1.0) { counter++; counterOfPixelsInROI++; actualPielValue = labelIterator3.Value(); sumSquared = sumSquared + actualPielValue; SumOfSquares = SumOfSquares + std::pow(actualPielValue,2); if(counterOfPixelsInROI == 1) { max = actualPielValue; min = actualPielValue; } if(actualPielValue >= max) { max = actualPielValue; } else if(actualPielValue <= min) { min = actualPielValue; } } } if (counter > 1) { sigma = ( SumOfSquares - std::pow( sumSquared, 2) / counter ) / ( counter-1 ); } else { sigma = 0; } } bool ImageStatisticsCalculator::ComputeStatistics( unsigned int timeStep ) { if (m_Image.IsNull() ) { mitkThrow() << "Image not set!"; } if (!m_Image->IsInitialized()) { mitkThrow() << "Image not initialized!"; } if ( m_Image->GetReferenceCount() == 1 ) { // Image no longer valid; we are the only ones to still hold a reference on it return false; } if ( timeStep >= m_Image->GetTimeSteps() ) { throw std::runtime_error( "Error: invalid time step!" ); } // If a mask was set but we are the only ones to still hold a reference on // it, delete it. if ( m_ImageMask.IsNotNull() && (m_ImageMask->GetReferenceCount() == 1) ) { m_ImageMask = nullptr; } // Check if statistics is already up-to-date unsigned long imageMTime = m_ImageStatisticsTimeStampVector[timeStep].GetMTime(); unsigned long maskedImageMTime = m_MaskedImageStatisticsTimeStampVector[timeStep].GetMTime(); unsigned long planarFigureMTime = m_PlanarFigureStatisticsTimeStampVector[timeStep].GetMTime(); bool imageStatisticsCalculationTrigger = m_ImageStatisticsCalculationTriggerVector[timeStep]; bool maskedImageStatisticsCalculationTrigger = m_MaskedImageStatisticsCalculationTriggerVector[timeStep]; bool planarFigureStatisticsCalculationTrigger = m_PlanarFigureStatisticsCalculationTriggerVector[timeStep]; if ( !m_IgnorePixelValueChanged && !m_HotspotRadiusInMMChanged && ((m_MaskingMode != MASKING_MODE_NONE) || (imageMTime > m_Image->GetMTime() && !imageStatisticsCalculationTrigger)) && ((m_MaskingMode != MASKING_MODE_IMAGE) || (maskedImageMTime > m_ImageMask->GetMTime() && !maskedImageStatisticsCalculationTrigger)) && ((m_MaskingMode != MASKING_MODE_PLANARFIGURE) || (planarFigureMTime > m_PlanarFigure->GetMTime() && !planarFigureStatisticsCalculationTrigger)) ) { // Statistics is up to date! if ( m_MaskingModeChanged ) { m_MaskingModeChanged = false; } else { return false; } } // Reset state changed flag m_MaskingModeChanged = false; m_IgnorePixelValueChanged = false; // Depending on masking mode, extract and/or generate the required image // and mask data from the user input this->ExtractImageAndMask( timeStep ); StatisticsContainer *statisticsContainer; HistogramContainer *histogramContainer; switch ( m_MaskingMode ) { case MASKING_MODE_NONE: default: if(!m_DoIgnorePixelValue) { statisticsContainer = &m_ImageStatisticsVector[timeStep]; histogramContainer = &m_ImageHistogramVector[timeStep]; m_ImageStatisticsTimeStampVector[timeStep].Modified(); m_ImageStatisticsCalculationTriggerVector[timeStep] = false; } else { statisticsContainer = &m_MaskedImageStatisticsVector[timeStep]; histogramContainer = &m_MaskedImageHistogramVector[timeStep]; m_MaskedImageStatisticsTimeStampVector[timeStep].Modified(); m_MaskedImageStatisticsCalculationTriggerVector[timeStep] = false; } break; case MASKING_MODE_IMAGE: statisticsContainer = &m_MaskedImageStatisticsVector[timeStep]; histogramContainer = &m_MaskedImageHistogramVector[timeStep]; m_MaskedImageStatisticsTimeStampVector[timeStep].Modified(); m_MaskedImageStatisticsCalculationTriggerVector[timeStep] = false; break; case MASKING_MODE_PLANARFIGURE: statisticsContainer = &m_PlanarFigureStatisticsVector[timeStep]; histogramContainer = &m_PlanarFigureHistogramVector[timeStep]; m_PlanarFigureStatisticsTimeStampVector[timeStep].Modified(); m_PlanarFigureStatisticsCalculationTriggerVector[timeStep] = false; break; } // Calculate statistics and histogram(s) if ( m_InternalImage->GetDimension() == 3 ) { if ( m_MaskingMode == MASKING_MODE_NONE && !m_DoIgnorePixelValue ) { AccessFixedDimensionByItk_2( m_InternalImage, InternalCalculateStatisticsUnmasked, 3, statisticsContainer, histogramContainer ); } else { AccessFixedDimensionByItk_3( m_InternalImage, InternalCalculateStatisticsMasked, 3, m_InternalImageMask3D.GetPointer(), statisticsContainer, histogramContainer ); } } else if ( m_InternalImage->GetDimension() == 2 ) { if ( m_MaskingMode == MASKING_MODE_NONE && !m_DoIgnorePixelValue ) { AccessFixedDimensionByItk_2( m_InternalImage, InternalCalculateStatisticsUnmasked, 2, statisticsContainer, histogramContainer ); } else { AccessFixedDimensionByItk_3( m_InternalImage, InternalCalculateStatisticsMasked, 2, m_InternalImageMask2D.GetPointer(), statisticsContainer, histogramContainer ); } } else { MITK_ERROR << "ImageStatistics: Image dimension not supported!"; } // Release unused image smart pointers to free memory m_InternalImage = mitk::Image::ConstPointer(); m_InternalImageMask3D = MaskImage3DType::Pointer(); m_InternalImageMask2D = MaskImage2DType::Pointer(); return true; } ImageStatisticsCalculator::BinFrequencyType ImageStatisticsCalculator::GetBinsAndFreuqencyForHistograms( unsigned int timeStep , unsigned int label ) const { const HistogramType *binsAndFrequencyToCalculate = this->GetHistogram(0); // ToDo: map should be created on stack not on heap std::map returnedHistogramMap; unsigned int size = binsAndFrequencyToCalculate->Size(); for( unsigned int bin=0; bin < size; ++bin ) { double frequency = binsAndFrequencyToCalculate->GetFrequency( bin, 0 ); //if( frequency > mitk::eps ) { returnedHistogramMap.insert( std::pair(binsAndFrequencyToCalculate->GetMeasurement( bin, 0 ), binsAndFrequencyToCalculate->GetFrequency( bin, 0 ) ) ); } } return returnedHistogramMap; } const ImageStatisticsCalculator::HistogramType * ImageStatisticsCalculator::GetHistogram( unsigned int timeStep, unsigned int label ) const { if ( m_Image.IsNull() || (timeStep >= m_Image->GetTimeSteps()) ) { return nullptr; } switch ( m_MaskingMode ) { case MASKING_MODE_NONE: default: { if(m_DoIgnorePixelValue) return m_MaskedImageHistogramVector[timeStep][label]; return m_ImageHistogramVector[timeStep][label]; } case MASKING_MODE_IMAGE: return m_MaskedImageHistogramVector[timeStep][label]; case MASKING_MODE_PLANARFIGURE: return m_PlanarFigureHistogramVector[timeStep][label]; } } const ImageStatisticsCalculator::HistogramContainer & ImageStatisticsCalculator::GetHistogramVector( unsigned int timeStep ) const { if ( m_Image.IsNull() || (timeStep >= m_Image->GetTimeSteps()) ) { return m_EmptyHistogramContainer; } switch ( m_MaskingMode ) { case MASKING_MODE_NONE: default: { if(m_DoIgnorePixelValue) return m_MaskedImageHistogramVector[timeStep]; return m_ImageHistogramVector[timeStep]; } case MASKING_MODE_IMAGE: return m_MaskedImageHistogramVector[timeStep]; case MASKING_MODE_PLANARFIGURE: return m_PlanarFigureHistogramVector[timeStep]; } } const ImageStatisticsCalculator::Statistics & ImageStatisticsCalculator::GetStatistics( unsigned int timeStep, unsigned int label ) const { if ( m_Image.IsNull() || (timeStep >= m_Image->GetTimeSteps()) ) { return m_EmptyStatistics; } switch ( m_MaskingMode ) { case MASKING_MODE_NONE: default: { if(m_DoIgnorePixelValue) return m_MaskedImageStatisticsVector[timeStep][label]; return m_ImageStatisticsVector[timeStep][label]; } case MASKING_MODE_IMAGE: return m_MaskedImageStatisticsVector[timeStep][label]; case MASKING_MODE_PLANARFIGURE: return m_PlanarFigureStatisticsVector[timeStep][label]; } } const ImageStatisticsCalculator::StatisticsContainer & ImageStatisticsCalculator::GetStatisticsVector( unsigned int timeStep ) const { if ( m_Image.IsNull() || (timeStep >= m_Image->GetTimeSteps()) ) { return m_EmptyStatisticsContainer; } switch ( m_MaskingMode ) { case MASKING_MODE_NONE: default: { if(m_DoIgnorePixelValue) return m_MaskedImageStatisticsVector[timeStep]; return m_ImageStatisticsVector[timeStep]; } case MASKING_MODE_IMAGE: return m_MaskedImageStatisticsVector[timeStep]; case MASKING_MODE_PLANARFIGURE: return m_PlanarFigureStatisticsVector[timeStep]; } } void ImageStatisticsCalculator::ExtractImageAndMask( unsigned int timeStep ) { if ( m_Image.IsNull() ) { throw std::runtime_error( "Error: image empty!" ); } if ( timeStep >= m_Image->GetTimeSteps() ) { throw std::runtime_error( "Error: invalid time step!" ); } ImageTimeSelector::Pointer imageTimeSelector = ImageTimeSelector::New(); imageTimeSelector->SetInput( m_Image ); imageTimeSelector->SetTimeNr( timeStep ); imageTimeSelector->UpdateLargestPossibleRegion(); mitk::Image *timeSliceImage = imageTimeSelector->GetOutput(); switch ( m_MaskingMode ) { case MASKING_MODE_NONE: { m_InternalImage = timeSliceImage; m_InternalImageMask2D = nullptr; m_InternalImageMask3D = nullptr; if(m_DoIgnorePixelValue) { if( m_InternalImage->GetDimension() == 3 ) { if(itk::ImageIOBase::USHORT != timeSliceImage->GetPixelType().GetComponentType()) CastToItkImage( timeSliceImage, m_InternalImageMask3D ); else CastToItkImage( timeSliceImage->Clone(), m_InternalImageMask3D ); m_InternalImageMask3D->FillBuffer(1); } if( m_InternalImage->GetDimension() == 2 ) { if(itk::ImageIOBase::USHORT != timeSliceImage->GetPixelType().GetComponentType()) CastToItkImage( timeSliceImage, m_InternalImageMask2D ); else CastToItkImage( timeSliceImage->Clone(), m_InternalImageMask2D ); m_InternalImageMask2D->FillBuffer(1); } } break; } case MASKING_MODE_IMAGE: { if ( m_ImageMask.IsNotNull() && (m_ImageMask->GetReferenceCount() > 1) ) { if ( timeStep >= m_ImageMask->GetTimeSteps() ) { // Use the last mask time step in case the current time step is bigger than the total // number of mask time steps. // It makes more sense setting this to the last mask time step than to 0. // For instance if you have a mask with 2 time steps and an image with 5: // If time step 0 is selected, the mask will use time step 0. // If time step 1 is selected, the mask will use time step 1. // If time step 2+ is selected, the mask will use time step 1. // If you have a mask with only one time step instead, this will always default to 0. timeStep = m_ImageMask->GetTimeSteps() - 1; } ImageTimeSelector::Pointer maskedImageTimeSelector = ImageTimeSelector::New(); maskedImageTimeSelector->SetInput( m_ImageMask ); maskedImageTimeSelector->SetTimeNr( timeStep ); maskedImageTimeSelector->UpdateLargestPossibleRegion(); mitk::Image *timeSliceMaskedImage = maskedImageTimeSelector->GetOutput(); m_InternalImage = timeSliceImage; CastToItkImage( timeSliceMaskedImage, m_InternalImageMask3D ); } else { throw std::runtime_error( "Error: image mask empty!" ); } break; } case MASKING_MODE_PLANARFIGURE: { m_InternalImageMask2D = nullptr; if ( m_PlanarFigure.IsNull() ) { throw std::runtime_error( "Error: planar figure empty!" ); } if ( !m_PlanarFigure->IsClosed() ) { throw std::runtime_error( "Masking not possible for non-closed figures" ); } const BaseGeometry *imageGeometry = timeSliceImage->GetGeometry(); if ( imageGeometry == nullptr ) { throw std::runtime_error( "Image geometry invalid!" ); } const PlaneGeometry *planarFigurePlaneGeometry = m_PlanarFigure->GetPlaneGeometry(); if ( planarFigurePlaneGeometry == nullptr ) { throw std::runtime_error( "Planar-Figure not yet initialized!" ); } const PlaneGeometry *planarFigureGeometry = dynamic_cast< const PlaneGeometry * >( planarFigurePlaneGeometry ); if ( planarFigureGeometry == nullptr ) { throw std::runtime_error( "Non-planar planar figures not supported!" ); } // Find principal direction of PlanarFigure in input image unsigned int axis; if ( !this->GetPrincipalAxis( imageGeometry, planarFigureGeometry->GetNormal(), axis ) ) { throw std::runtime_error( "Non-aligned planar figures not supported!" ); } m_PlanarFigureAxis = axis; // Find slice number corresponding to PlanarFigure in input image MaskImage3DType::IndexType index; imageGeometry->WorldToIndex( planarFigureGeometry->GetOrigin(), index ); unsigned int slice = index[axis]; m_PlanarFigureSlice = slice; // Extract slice with given position and direction from image unsigned int dimension = timeSliceImage->GetDimension(); if (dimension != 2) { ExtractImageFilter::Pointer imageExtractor = ExtractImageFilter::New(); imageExtractor->SetInput( timeSliceImage ); imageExtractor->SetSliceDimension( axis ); imageExtractor->SetSliceIndex( slice ); imageExtractor->Update(); m_InternalImage = imageExtractor->GetOutput(); } else { m_InternalImage = timeSliceImage; } // Compute mask from PlanarFigure AccessFixedDimensionByItk_1( m_InternalImage, InternalCalculateMaskFromPlanarFigure, 2, axis ); } } if(m_DoIgnorePixelValue) { if ( m_InternalImage->GetDimension() == 3 ) { AccessFixedDimensionByItk_1( m_InternalImage, InternalMaskIgnoredPixels, 3, m_InternalImageMask3D.GetPointer() ); } else if ( m_InternalImage->GetDimension() == 2 ) { AccessFixedDimensionByItk_1( m_InternalImage, InternalMaskIgnoredPixels, 2, m_InternalImageMask2D.GetPointer() ); } } } bool ImageStatisticsCalculator::GetPrincipalAxis( const BaseGeometry *geometry, Vector3D vector, unsigned int &axis ) { vector.Normalize(); for ( unsigned int i = 0; i < 3; ++i ) { Vector3D axisVector = geometry->GetAxisVector( i ); axisVector.Normalize(); if ( fabs( fabs( axisVector * vector ) - 1.0) < mitk::eps ) { axis = i; return true; } } return false; } unsigned int ImageStatisticsCalculator::calcNumberOfBins(mitk::ScalarType min, mitk::ScalarType max) { return std::ceil( ( (max - min ) / m_HistogramBinSize) ); } template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator::InternalCalculateStatisticsUnmasked( const itk::Image< TPixel, VImageDimension > *image, StatisticsContainer *statisticsContainer, HistogramContainer* histogramContainer ) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef typename ImageType::IndexType IndexType; typedef itk::Statistics::ScalarImageToHistogramGenerator< ImageType > HistogramGeneratorType; statisticsContainer->clear(); histogramContainer->clear(); // Progress listening... typedef itk::SimpleMemberCommand< ImageStatisticsCalculator > ITKCommandType; ITKCommandType::Pointer progressListener; progressListener = ITKCommandType::New(); progressListener->SetCallbackFunction( this, &ImageStatisticsCalculator::UnmaskedStatisticsProgressUpdate ); // Issue 100 artificial progress events since ScalarIMageToHistogramGenerator // does not (yet?) support progress reporting this->InvokeEvent( itk::StartEvent() ); for ( unsigned int i = 0; i < 100; ++i ) { this->UnmaskedStatisticsProgressUpdate(); } // Calculate statistics (separate filter) typedef itk::ExtendedStatisticsImageFilter< ImageType > StatisticsFilterType; typename StatisticsFilterType::Pointer statisticsFilter = StatisticsFilterType::New(); statisticsFilter->SetInput( image ); statisticsFilter->SetBinSize( 100 ); statisticsFilter->SetCoordinateTolerance( 0.001 ); statisticsFilter->SetDirectionTolerance( 0.001 ); unsigned long observerTag = statisticsFilter->AddObserver( itk::ProgressEvent(), progressListener ); try { statisticsFilter->Update(); } catch (const itk::ExceptionObject& e) { mitkThrow() << "Image statistics calculation failed due to following ITK Exception: \n " << e.what(); } catch( const std::exception& e ) { //mitkThrow() << "Image statistics calculation failed due to following ITK Exception: \n " << e.what(); } statisticsFilter->RemoveObserver( observerTag ); this->InvokeEvent( itk::EndEvent() ); // Calculate minimum and maximum typedef itk::MinimumMaximumImageCalculator< ImageType > MinMaxFilterType; typename MinMaxFilterType::Pointer minMaxFilter = MinMaxFilterType::New(); minMaxFilter->SetImage( image ); unsigned long observerTag2 = minMaxFilter->AddObserver( itk::ProgressEvent(), progressListener ); minMaxFilter->Compute(); minMaxFilter->RemoveObserver( observerTag2 ); this->InvokeEvent( itk::EndEvent() ); Statistics statistics; statistics.Reset(); statistics.SetLabel(1); statistics.SetN(image->GetBufferedRegion().GetNumberOfPixels()); statistics.SetMin(statisticsFilter->GetMinimum()); statistics.SetMax(statisticsFilter->GetMaximum()); statistics.SetMean(statisticsFilter->GetMean()); statistics.SetMedian(statisticsFilter->GetMedian()); statistics.SetVariance(statisticsFilter->GetVariance()); statistics.SetSkewness(statisticsFilter->GetSkewness()); statistics.SetKurtosis(statisticsFilter->GetKurtosis()); statistics.SetUniformity( statisticsFilter->GetUniformity()); statistics.SetEntropy( statisticsFilter->GetEntropy()); statistics.SetUPP( statisticsFilter->GetUPP()); statistics.SetMPP( statisticsFilter->GetMPP()); statistics.SetSigma(statisticsFilter->GetSigma()); statistics.SetRMS(sqrt( statistics.GetMean() * statistics.GetMean() + statistics.GetSigma() * statistics.GetSigma() )); statistics.GetMinIndex().set_size(image->GetImageDimension()); statistics.GetMaxIndex().set_size(image->GetImageDimension()); vnl_vector tmpMaxIndex; vnl_vector tmpMinIndex; tmpMaxIndex.set_size(image->GetImageDimension() ); tmpMinIndex.set_size(image->GetImageDimension() ); for (unsigned int i=0; iGetIndexOfMaximum()[i]; tmpMinIndex[i] = minMaxFilter->GetIndexOfMinimum()[i]; } statistics.SetMinIndex(tmpMaxIndex); statistics.SetMinIndex(tmpMinIndex); if( IsHotspotCalculated() && VImageDimension == 3 ) { typedef itk::Image< unsigned short, VImageDimension > MaskImageType; typename MaskImageType::Pointer nullMask; bool isHotspotDefined(false); Statistics hotspotStatistics = this->CalculateHotspotStatistics(image, nullMask.GetPointer(), m_HotspotRadiusInMM, isHotspotDefined, 0 ); if (isHotspotDefined) { statistics.SetHasHotspotStatistics(true); statistics.GetHotspotStatistics() = hotspotStatistics; } else { statistics.SetHasHotspotStatistics(false); } if(statistics.GetHotspotStatistics().HasHotspotStatistics() ) { MITK_DEBUG << "Hotspot statistics available"; statistics.SetHotspotIndex(hotspotStatistics.GetHotspotIndex()); } else { MITK_ERROR << "No hotspot statistics available!"; } } statisticsContainer->push_back( statistics ); // Calculate histogram // calculate bin size or number of bins unsigned int numberOfBins = 200; // default number of bins if (m_UseDefaultBinSize) { m_HistogramBinSize = std::ceil( (statistics.GetMax() - statistics.GetMin() + 1)/numberOfBins ); } else { numberOfBins = calcNumberOfBins(statistics.GetMin(), statistics.GetMax()); } typename HistogramGeneratorType::Pointer histogramGenerator = HistogramGeneratorType::New(); histogramGenerator->SetInput( image ); histogramGenerator->SetMarginalScale( 100 ); histogramGenerator->SetNumberOfBins( numberOfBins ); histogramGenerator->SetHistogramMin( statistics.GetMin() ); histogramGenerator->SetHistogramMax( statistics.GetMax() ); histogramGenerator->Compute(); histogramContainer->push_back( histogramGenerator->GetOutput() ); } template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator::InternalMaskIgnoredPixels( const itk::Image< TPixel, VImageDimension > *image, itk::Image< unsigned short, VImageDimension > *maskImage ) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< unsigned short, VImageDimension > MaskImageType; itk::ImageRegionIterator itmask(maskImage, maskImage->GetLargestPossibleRegion()); itk::ImageRegionConstIterator itimage(image, image->GetLargestPossibleRegion()); itmask.GoToBegin(); itimage.GoToBegin(); while( !itmask.IsAtEnd() ) { if(m_IgnorePixelValue == itimage.Get()) { itmask.Set(0); } ++itmask; ++itimage; } } template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator::InternalCalculateStatisticsMasked( const itk::Image< TPixel, VImageDimension > *image, itk::Image< unsigned short, VImageDimension > *maskImage, StatisticsContainer* statisticsContainer, HistogramContainer* histogramContainer ) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< unsigned short, VImageDimension > MaskImageType; typedef typename ImageType::IndexType IndexType; typedef typename ImageType::PointType PointType; typedef typename ImageType::SpacingType SpacingType; typedef typename ImageType::Pointer ImagePointer; typedef itk::ExtendedLabelStatisticsImageFilter< ImageType, MaskImageType > LabelStatisticsFilterType; typedef itk::ChangeInformationImageFilter< MaskImageType > ChangeInformationFilterType; typedef itk::ExtractImageFilter< ImageType, ImageType > ExtractImageFilterType; statisticsContainer->clear(); histogramContainer->clear(); // Make sure that mask is set if ( maskImage == nullptr ) { itkExceptionMacro( << "Mask image needs to be set!" ); } // Make sure that spacing of mask and image are the same //SpacingType imageSpacing = image->GetSpacing(); //SpacingType maskSpacing = maskImage->GetSpacing(); //PointType zeroPoint; zeroPoint.Fill( 0.0 ); //if ( (zeroPoint + imageSpacing).SquaredEuclideanDistanceTo( (zeroPoint + maskSpacing) ) > mitk::eps ) //{ // itkExceptionMacro( << "Mask needs to have same spacing as image! (Image spacing: " << imageSpacing << "; Mask spacing: " << maskSpacing << ")" ); //} // Make sure that orientation of mask and image are the same typedef typename ImageType::DirectionType DirectionType; DirectionType imageDirection = image->GetDirection(); DirectionType maskDirection = maskImage->GetDirection(); for( int i = 0; i < imageDirection.ColumnDimensions; ++i ) { for( int j = 0; j < imageDirection.ColumnDimensions; ++j ) { double differenceDirection = imageDirection[i][j] - maskDirection[i][j]; if ( fabs( differenceDirection ) > mitk::eps ) { double differenceDirection = imageDirection[i][j] - maskDirection[i][j]; if ( fabs( differenceDirection ) > 0.001 /*mitk::eps*/ ) // TODO: temp fix (bug 17121) { itkExceptionMacro( << "Mask needs to have same direction as image! (Image direction: " << imageDirection << "; Mask direction: " << maskDirection << ")" ); } } } } // Make sure that the voxels of mask and image are correctly "aligned", i.e., voxel boundaries are the same in both images PointType imageOrigin = image->GetOrigin(); PointType maskOrigin = maskImage->GetOrigin(); long offset[ImageType::ImageDimension]; typedef itk::ContinuousIndex ContinousIndexType; ContinousIndexType maskOriginContinousIndex, imageOriginContinousIndex; image->TransformPhysicalPointToContinuousIndex(maskOrigin, maskOriginContinousIndex); image->TransformPhysicalPointToContinuousIndex(imageOrigin, imageOriginContinousIndex); for ( unsigned int i = 0; i < ImageType::ImageDimension; ++i ) { double misalignment = maskOriginContinousIndex[i] - floor( maskOriginContinousIndex[i] + 0.5 ); if ( fabs( misalignment ) > mitk::eps ) { itkWarningMacro( << "Pixels/voxels of mask and image are not sufficiently aligned! (Misalignment: " << misalignment << ")" ); } double indexCoordDistance = maskOriginContinousIndex[i] - imageOriginContinousIndex[i]; offset[i] = int( indexCoordDistance + image->GetBufferedRegion().GetIndex()[i] + 0.5 ); } // Adapt the origin and region (index/size) of the mask so that the origin of both are the same typename ChangeInformationFilterType::Pointer adaptMaskFilter; adaptMaskFilter = ChangeInformationFilterType::New(); adaptMaskFilter->ChangeOriginOn(); adaptMaskFilter->ChangeRegionOn(); adaptMaskFilter->SetInput( maskImage ); adaptMaskFilter->SetOutputOrigin( image->GetOrigin() ); adaptMaskFilter->SetOutputOffset( offset ); adaptMaskFilter->SetCoordinateTolerance( 0.001 ); adaptMaskFilter->SetDirectionTolerance( 0.001 ); typename MaskImageType::Pointer adaptedMaskImage; try { adaptMaskFilter->Update(); adaptedMaskImage = adaptMaskFilter->GetOutput(); } catch( const itk::ExceptionObject &e) { mitkThrow() << "Attempt to adapt shifted origin of the mask image failed due to ITK Exception: \n" << e.what(); } catch( const std::exception& e ) { //mitkThrow() << "Image statistics calculation failed due to following ITK Exception: \n " << e.what(); } // Make sure that mask region is contained within image region if ( adaptedMaskImage.IsNotNull() && !image->GetLargestPossibleRegion().IsInside( adaptedMaskImage->GetLargestPossibleRegion() ) ) { itkWarningMacro( << "Mask region needs to be inside of image region! (Image region: " << image->GetLargestPossibleRegion() << "; Mask region: " << adaptedMaskImage->GetLargestPossibleRegion() << ")" ); } // If mask region is smaller than image region, extract the sub-sampled region from the original image typename ImageType::SizeType imageSize = image->GetBufferedRegion().GetSize(); typename ImageType::SizeType maskSize = maskImage->GetBufferedRegion().GetSize(); bool maskSmallerImage = false; for ( unsigned int i = 0; i < ImageType::ImageDimension; ++i ) { if ( maskSize[i] < imageSize[i] ) { maskSmallerImage = true; } } typename ImageType::ConstPointer adaptedImage; if ( maskSmallerImage ) { typename ExtractImageFilterType::Pointer extractImageFilter = ExtractImageFilterType::New(); extractImageFilter->SetInput( image ); extractImageFilter->SetExtractionRegion( adaptedMaskImage->GetBufferedRegion() ); extractImageFilter->SetCoordinateTolerance( 0.001 ); extractImageFilter->SetDirectionTolerance( 0.001 ); extractImageFilter->Update(); adaptedImage = extractImageFilter->GetOutput(); } else { adaptedImage = image; } // Initialize Filter typedef itk::StatisticsImageFilter< ImageType > StatisticsFilterType; typename StatisticsFilterType::Pointer statisticsFilter = StatisticsFilterType::New(); statisticsFilter->SetInput( adaptedImage ); try { statisticsFilter->Update(); } catch( const itk::ExceptionObject& e) { mitkThrow() << "Image statistics initialization computation failed with ITK Exception: \n " << e.what(); } catch( const std::exception& e ) { //mitkThrow() << "Image statistics calculation failed due to following ITK Exception: \n " << e.what(); } // Calculate bin size or number of bins - unsigned int numberOfBins = 200; // default number of bins + unsigned int numberOfBins = m_HistogramBinSize; // default number of bins double maximum = 0.0; double minimum = 0.0; if (m_UseBinSizeBasedOnVOIRegion) { maximum = statisticsFilter->GetMaximum(); minimum = statisticsFilter->GetMinimum(); if (m_UseDefaultBinSize) { m_HistogramBinSize = std::ceil( static_cast((statisticsFilter->GetMaximum() - statisticsFilter->GetMinimum() + 1)/numberOfBins) ); } else { numberOfBins = calcNumberOfBins(statisticsFilter->GetMinimum(), statisticsFilter->GetMaximum()); } } else { double sig = 0.0; int counter = 0; //Find the min and max values for the Roi to set the range for the histogram GetMinAndMaxValue( minimum, maximum, counter, sig, image, maskImage); - numberOfBins = maximum - minimum; - if(maximum - minimum <= 10) - { - numberOfBins = 100; - } +// numberOfBins = maximum - minimum; +// if(maximum - minimum <= 10) +// { +// numberOfBins = 100; +// } } typename LabelStatisticsFilterType::Pointer labelStatisticsFilter = LabelStatisticsFilterType::New(); labelStatisticsFilter->SetInput( adaptedImage ); labelStatisticsFilter->SetLabelInput( adaptedMaskImage ); labelStatisticsFilter->SetCoordinateTolerance( 0.001 ); labelStatisticsFilter->SetDirectionTolerance( 0.001 ); labelStatisticsFilter->UseHistogramsOn(); labelStatisticsFilter->SetHistogramParameters( numberOfBins, floor(minimum), ceil(maximum) ); //statisticsFilter->GetMinimum() statisticsFilter->GetMaximum() // Add progress listening typedef itk::SimpleMemberCommand< ImageStatisticsCalculator > ITKCommandType; ITKCommandType::Pointer progressListener; progressListener = ITKCommandType::New(); progressListener->SetCallbackFunction( this, &ImageStatisticsCalculator::MaskedStatisticsProgressUpdate ); unsigned long observerTag = labelStatisticsFilter->AddObserver( itk::ProgressEvent(), progressListener ); // Execute filter this->InvokeEvent( itk::StartEvent() ); // Make sure that only the mask region is considered (otherwise, if the mask region is smaller // than the image region, the Update() would result in an exception). labelStatisticsFilter->GetOutput()->SetRequestedRegion( adaptedMaskImage->GetLargestPossibleRegion() ); // Execute the filter try { labelStatisticsFilter->Update(); } catch( const itk::ExceptionObject& e) { mitkThrow() << "Image statistics calculation failed due to following ITK Exception: \n " << e.what(); } catch( const std::exception& e ) { //mitkThrow() << "Image statistics calculation failed due to following ITK Exception: \n " << e.what(); } this->InvokeEvent( itk::EndEvent() ); if( observerTag ) labelStatisticsFilter->RemoveObserver( observerTag ); // Find all relevant labels of mask (other than 0) std::list< int > relevantLabels = labelStatisticsFilter->GetRelevantLabels(); unsigned int i; if ( labelStatisticsFilter->GetMaskingNonEmpty() ) { std::list< int >::iterator it; for ( it = relevantLabels.begin(), i = 0; it != relevantLabels.end(); ++it, ++i ) { Statistics statistics; // restore previous code labelStatisticsFilter->GetHistogram(*it) ; histogramContainer->push_back( HistogramType::ConstPointer( labelStatisticsFilter->GetHistogram( (*it) ) ) ); statistics.SetLabel (*it); statistics.SetN(labelStatisticsFilter->GetCount( *it )); statistics.SetMin(labelStatisticsFilter->GetMinimum( *it )); statistics.SetMax(labelStatisticsFilter->GetMaximum( *it )); statistics.SetMean(labelStatisticsFilter->GetMean( *it )); statistics.SetMedian(labelStatisticsFilter->GetMedian( *it)); statistics.SetMedian(labelStatisticsFilter->GetMedian( *it )); statistics.SetVariance(labelStatisticsFilter->GetVariance( *it )); statistics.SetSigma(labelStatisticsFilter->GetSigma( *it )); statistics.SetSkewness(labelStatisticsFilter->GetSkewness( *it )); statistics.SetKurtosis(labelStatisticsFilter->GetKurtosis( *it )); statistics.SetUniformity( labelStatisticsFilter->GetUniformity( *it )); statistics.SetEntropy( labelStatisticsFilter->GetEntropy( *it )); statistics.SetUPP( labelStatisticsFilter->GetUPP( *it)); statistics.SetMPP( labelStatisticsFilter->GetMPP( *it)); statistics.SetRMS(sqrt( statistics.GetMean() * statistics.GetMean() + statistics.GetSigma() * statistics.GetSigma() )); // restrict image to mask area for min/max index calculation typedef itk::MaskImageFilter< ImageType, MaskImageType, ImageType > MaskImageFilterType; typename MaskImageFilterType::Pointer masker = MaskImageFilterType::New(); bool isMinAndMaxSameValue = (statistics.GetMin() == statistics.GetMax()); // bug 17962: following is a workaround for the case when min and max are the same, we can probably find a nicer way here double outsideValue = (isMinAndMaxSameValue ? (statistics.GetMax()/2) : (statistics.GetMin()+statistics.GetMax())/2); masker->SetOutsideValue( outsideValue ); masker->SetInput1(adaptedImage); masker->SetInput2(adaptedMaskImage); masker->SetCoordinateTolerance( 0.001 ); masker->SetDirectionTolerance( 0.001 ); masker->Update(); // get index of minimum and maximum typedef itk::MinimumMaximumImageCalculator< ImageType > MinMaxFilterType; typename MinMaxFilterType::Pointer minMaxFilter = MinMaxFilterType::New(); minMaxFilter->SetImage( masker->GetOutput() ); unsigned long observerTag2 = minMaxFilter->AddObserver( itk::ProgressEvent(), progressListener ); minMaxFilter->Compute(); minMaxFilter->RemoveObserver( observerTag2 ); this->InvokeEvent( itk::EndEvent() ); typename MinMaxFilterType::IndexType tempMaxIndex = minMaxFilter->GetIndexOfMaximum(); // bug 17962: following is a workaround for the case when min and max are the same, we can probably find a nicer way here typename MinMaxFilterType::IndexType tempMinIndex = (isMinAndMaxSameValue ? minMaxFilter->GetIndexOfMaximum() : minMaxFilter->GetIndexOfMinimum()); // FIX BUG 14644 //If a PlanarFigure is used for segmentation the //adaptedImage is a single slice (2D). Adding the // 3. dimension. vnl_vector maxIndex; vnl_vector minIndex; maxIndex.set_size(m_Image->GetDimension()); minIndex.set_size(m_Image->GetDimension()); if (m_MaskingMode == MASKING_MODE_PLANARFIGURE && m_Image->GetDimension()==3) { maxIndex[m_PlanarFigureCoordinate0] = tempMaxIndex[0]; maxIndex[m_PlanarFigureCoordinate1] = tempMaxIndex[1]; maxIndex[m_PlanarFigureAxis] = m_PlanarFigureSlice; minIndex[m_PlanarFigureCoordinate0] = tempMinIndex[0] ; minIndex[m_PlanarFigureCoordinate1] = tempMinIndex[1]; minIndex[m_PlanarFigureAxis] = m_PlanarFigureSlice; } else { for (unsigned int i = 0; ipush_back( statistics ); } } else { histogramContainer->push_back( HistogramType::ConstPointer( m_EmptyHistogram ) ); statisticsContainer->push_back( Statistics() ); } } template ImageStatisticsCalculator::ImageExtrema ImageStatisticsCalculator::CalculateExtremaWorld( const itk::Image *inputImage, itk::Image *maskImage, double neccessaryDistanceToImageBorderInMM, unsigned int label) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< unsigned short, VImageDimension > MaskImageType; typedef itk::ImageRegionConstIteratorWithIndex MaskImageIteratorType; typedef itk::ImageRegionConstIteratorWithIndex InputImageIndexIteratorType; typename ImageType::SpacingType spacing = inputImage->GetSpacing(); ImageExtrema minMax; minMax.Defined = false; minMax.MaxIndex.set_size(VImageDimension); minMax.MaxIndex.set_size(VImageDimension); typename ImageType::RegionType allowedExtremaRegion = inputImage->GetLargestPossibleRegion(); bool keepDistanceToImageBorders( neccessaryDistanceToImageBorderInMM > 0 ); if (keepDistanceToImageBorders) { long distanceInPixels[VImageDimension]; for(unsigned short dimension = 0; dimension < VImageDimension; ++dimension) { // To confirm that the whole hotspot is inside the image we have to keep a specific distance to the image-borders, which is as long as // the radius. To get the amount of indices we divide the radius by spacing and add 0.5 because voxels are center based: // For example with a radius of 2.2 and a spacing of 1 two indices are enough because 2.2 / 1 + 0.5 = 2.7 => 2. // But with a radius of 2.7 we need 3 indices because 2.7 / 1 + 0.5 = 3.2 => 3 distanceInPixels[dimension] = int( neccessaryDistanceToImageBorderInMM / spacing[dimension] + 0.5); } allowedExtremaRegion.ShrinkByRadius(distanceInPixels); } InputImageIndexIteratorType imageIndexIt(inputImage, allowedExtremaRegion); float maxValue = itk::NumericTraits::min(); float minValue = itk::NumericTraits::max(); typename ImageType::IndexType maxIndex; typename ImageType::IndexType minIndex; for(unsigned short i = 0; i < VImageDimension; ++i) { maxIndex[i] = 0; minIndex[i] = 0; } if (maskImage != nullptr) { MaskImageIteratorType maskIt(maskImage, maskImage->GetLargestPossibleRegion()); typename ImageType::IndexType imageIndex; typename ImageType::PointType worldPosition; typename ImageType::IndexType maskIndex; for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt) { imageIndex = maskIndex = maskIt.GetIndex(); if(maskIt.Get() == label) { if( allowedExtremaRegion.IsInside(imageIndex) ) { imageIndexIt.SetIndex( imageIndex ); double value = imageIndexIt.Get(); minMax.Defined = true; //Calculate minimum, maximum and corresponding index-values if( value > maxValue ) { maxIndex = imageIndexIt.GetIndex(); maxValue = value; } if(value < minValue ) { minIndex = imageIndexIt.GetIndex(); minValue = value; } } } } } else { for(imageIndexIt.GoToBegin(); !imageIndexIt.IsAtEnd(); ++imageIndexIt) { double value = imageIndexIt.Get(); minMax.Defined = true; //Calculate minimum, maximum and corresponding index-values if( value > maxValue ) { maxIndex = imageIndexIt.GetIndex(); maxValue = value; } if(value < minValue ) { minIndex = imageIndexIt.GetIndex(); minValue = value; } } } minMax.MaxIndex.set_size(VImageDimension); minMax.MinIndex.set_size(VImageDimension); for(unsigned int i = 0; i < minMax.MaxIndex.size(); ++i) { minMax.MaxIndex[i] = maxIndex[i]; } for(unsigned int i = 0; i < minMax.MinIndex.size(); ++i) { minMax.MinIndex[i] = minIndex[i]; } minMax.Max = maxValue; minMax.Min = minValue; return minMax; } template itk::Size ImageStatisticsCalculator ::CalculateConvolutionKernelSize(double spacing[VImageDimension], double radiusInMM) { typedef itk::Image< float, VImageDimension > KernelImageType; typedef typename KernelImageType::SizeType SizeType; SizeType maskSize; for(unsigned int i = 0; i < VImageDimension; ++i) { maskSize[i] = static_cast( 2 * radiusInMM / spacing[i]); // We always want an uneven size to have a clear center point in the convolution mask if(maskSize[i] % 2 == 0 ) { ++maskSize[i]; } } return maskSize; } template itk::SmartPointer< itk::Image > ImageStatisticsCalculator ::GenerateHotspotSearchConvolutionKernel(double mmPerPixel[VImageDimension], double radiusInMM) { std::stringstream ss; for (unsigned int i = 0; i < VImageDimension; ++i) { ss << mmPerPixel[i]; if (i < VImageDimension -1) ss << ","; } MITK_DEBUG << "Update convolution kernel for spacing (" << ss.str() << ") and radius " << radiusInMM << "mm"; double radiusInMMSquared = radiusInMM * radiusInMM; typedef itk::Image< float, VImageDimension > KernelImageType; typename KernelImageType::Pointer convolutionKernel = KernelImageType::New(); // Calculate size and allocate mask image typedef typename KernelImageType::SizeType SizeType; SizeType maskSize = this->CalculateConvolutionKernelSize(mmPerPixel, radiusInMM); Point3D convolutionMaskCenterIndex; convolutionMaskCenterIndex.Fill(0.0); for(unsigned int i = 0; i < VImageDimension; ++i) { convolutionMaskCenterIndex[i] = 0.5 * (double)(maskSize[i]-1); } typedef typename KernelImageType::IndexType IndexType; IndexType maskIndex; maskIndex.Fill(0); typedef typename KernelImageType::RegionType RegionType; RegionType maskRegion; maskRegion.SetSize(maskSize); maskRegion.SetIndex(maskIndex); convolutionKernel->SetRegions(maskRegion); convolutionKernel->SetSpacing(mmPerPixel); convolutionKernel->Allocate(); // Fill mask image values by subsampling the image grid typedef itk::ImageRegionIteratorWithIndex MaskIteratorType; MaskIteratorType maskIt(convolutionKernel,maskRegion); int numberOfSubVoxelsPerDimension = 2; // per dimension! int numberOfSubVoxels = ::pow( static_cast(numberOfSubVoxelsPerDimension), static_cast(VImageDimension) ); double subVoxelSizeInPixels = 1.0 / (double)numberOfSubVoxelsPerDimension; double valueOfOneSubVoxel = 1.0 / (double)numberOfSubVoxels; double maskValue = 0.0; Point3D subVoxelIndexPosition; double distanceSquared = 0.0; typedef itk::ContinuousIndex ContinuousIndexType; for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt) { ContinuousIndexType indexPoint(maskIt.GetIndex()); Point3D voxelPosition; for (unsigned int dimension = 0; dimension < VImageDimension; ++dimension) { voxelPosition[dimension] = indexPoint[dimension]; } maskValue = 0.0; Vector3D subVoxelOffset; subVoxelOffset.Fill(0.0); // iterate sub-voxels by iterating all possible offsets for (subVoxelOffset[0] = -0.5 + subVoxelSizeInPixels / 2.0; subVoxelOffset[0] < +0.5; subVoxelOffset[0] += subVoxelSizeInPixels) { for (subVoxelOffset[1] = -0.5 + subVoxelSizeInPixels / 2.0; subVoxelOffset[1] < +0.5; subVoxelOffset[1] += subVoxelSizeInPixels) { for (subVoxelOffset[2] = -0.5 + subVoxelSizeInPixels / 2.0; subVoxelOffset[2] < +0.5; subVoxelOffset[2] += subVoxelSizeInPixels) { subVoxelIndexPosition = voxelPosition + subVoxelOffset; // this COULD be integrated into the for-loops if neccessary (add voxelPosition to initializer and end condition) distanceSquared = (subVoxelIndexPosition[0]-convolutionMaskCenterIndex[0]) * mmPerPixel[0] * (subVoxelIndexPosition[0]-convolutionMaskCenterIndex[0]) * mmPerPixel[0] + (subVoxelIndexPosition[1]-convolutionMaskCenterIndex[1]) * mmPerPixel[1] * (subVoxelIndexPosition[1]-convolutionMaskCenterIndex[1]) * mmPerPixel[1] + (subVoxelIndexPosition[2]-convolutionMaskCenterIndex[2]) * mmPerPixel[2] * (subVoxelIndexPosition[2]-convolutionMaskCenterIndex[2]) * mmPerPixel[2]; if (distanceSquared <= radiusInMMSquared) { maskValue += valueOfOneSubVoxel; } } } } maskIt.Set( maskValue ); } return convolutionKernel; } template itk::SmartPointer > ImageStatisticsCalculator::GenerateConvolutionImage( const itk::Image* inputImage ) { double mmPerPixel[VImageDimension]; for (unsigned int dimension = 0; dimension < VImageDimension; ++dimension) { mmPerPixel[dimension] = inputImage->GetSpacing()[dimension]; } // update convolution kernel typedef itk::Image< float, VImageDimension > KernelImageType; typename KernelImageType::Pointer convolutionKernel = this->GenerateHotspotSearchConvolutionKernel(mmPerPixel, m_HotspotRadiusInMM); // update convolution image typedef itk::Image< TPixel, VImageDimension > InputImageType; typedef itk::Image< TPixel, VImageDimension > ConvolutionImageType; typedef itk::FFTConvolutionImageFilter ConvolutionFilterType; typename ConvolutionFilterType::Pointer convolutionFilter = ConvolutionFilterType::New(); typedef itk::ConstantBoundaryCondition BoundaryConditionType; BoundaryConditionType boundaryCondition; boundaryCondition.SetConstant(0.0); if (GetHotspotMustBeCompletlyInsideImage()) { // overwrite default boundary condition convolutionFilter->SetBoundaryCondition(&boundaryCondition); } convolutionFilter->SetInput(inputImage); convolutionFilter->SetKernelImage(convolutionKernel); convolutionFilter->SetNormalize(true); MITK_DEBUG << "Update Convolution image for hotspot search"; convolutionFilter->UpdateLargestPossibleRegion(); typename ConvolutionImageType::Pointer convolutionImage = convolutionFilter->GetOutput(); convolutionImage->SetSpacing( inputImage->GetSpacing() ); // only workaround because convolution filter seems to ignore spacing of input image m_HotspotRadiusInMMChanged = false; return convolutionImage; } template < typename TPixel, unsigned int VImageDimension> void ImageStatisticsCalculator ::FillHotspotMaskPixels( itk::Image* maskImage, itk::Point sphereCenter, double sphereRadiusInMM) { typedef itk::Image< TPixel, VImageDimension > MaskImageType; typedef itk::ImageRegionIteratorWithIndex MaskImageIteratorType; MaskImageIteratorType maskIt(maskImage, maskImage->GetLargestPossibleRegion()); typename MaskImageType::IndexType maskIndex; typename MaskImageType::PointType worldPosition; for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt) { maskIndex = maskIt.GetIndex(); maskImage->TransformIndexToPhysicalPoint(maskIndex, worldPosition); maskIt.Set( worldPosition.EuclideanDistanceTo(sphereCenter) <= sphereRadiusInMM ? 1 : 0 ); } } template < typename TPixel, unsigned int VImageDimension> ImageStatisticsCalculator::Statistics ImageStatisticsCalculator::CalculateHotspotStatistics( const itk::Image* inputImage, itk::Image* maskImage, double radiusInMM, bool& isHotspotDefined, unsigned int label) { // get convolution image (updated in GenerateConvolutionImage()) typedef itk::Image< TPixel, VImageDimension > InputImageType; typedef itk::Image< TPixel, VImageDimension > ConvolutionImageType; typedef itk::Image< float, VImageDimension > KernelImageType; typedef itk::Image< unsigned short, VImageDimension > MaskImageType; //typename ConvolutionImageType::Pointer convolutionImage = dynamic_cast(this->GenerateConvolutionImage(inputImage)); typename ConvolutionImageType::Pointer convolutionImage = this->GenerateConvolutionImage(inputImage); if (convolutionImage.IsNull()) { MITK_ERROR << "Empty convolution image in CalculateHotspotStatistics(). We should never reach this state (logic error)."; throw std::logic_error("Empty convolution image in CalculateHotspotStatistics()"); } // find maximum in convolution image, given the current mask double requiredDistanceToBorder = m_HotspotMustBeCompletelyInsideImage ? m_HotspotRadiusInMM : -1.0; ImageExtrema convolutionImageInformation = CalculateExtremaWorld(convolutionImage.GetPointer(), maskImage, requiredDistanceToBorder, label); isHotspotDefined = convolutionImageInformation.Defined; if (!isHotspotDefined) { m_EmptyStatistics.Reset(VImageDimension); MITK_ERROR << "No origin of hotspot-sphere was calculated! Returning empty statistics"; return m_EmptyStatistics; } else { // create a binary mask around the "hotspot" region, fill the shape of a sphere around our hotspot center typedef itk::ImageDuplicator< InputImageType > DuplicatorType; typename DuplicatorType::Pointer copyMachine = DuplicatorType::New(); copyMachine->SetInputImage(inputImage); copyMachine->Update(); typedef itk::CastImageFilter< InputImageType, MaskImageType > CastFilterType; typename CastFilterType::Pointer caster = CastFilterType::New(); caster->SetInput( copyMachine->GetOutput() ); caster->Update(); typename MaskImageType::Pointer hotspotMaskITK = caster->GetOutput(); typedef typename InputImageType::IndexType IndexType; IndexType maskCenterIndex; for (unsigned int d =0; d< VImageDimension;++d) maskCenterIndex[d]=convolutionImageInformation.MaxIndex[d]; typename ConvolutionImageType::PointType maskCenter; inputImage->TransformIndexToPhysicalPoint(maskCenterIndex,maskCenter); this->FillHotspotMaskPixels(hotspotMaskITK.GetPointer(), maskCenter, radiusInMM); // calculate statistics within the binary mask typedef itk::ExtendedLabelStatisticsImageFilter< InputImageType, MaskImageType> LabelStatisticsFilterType; typename LabelStatisticsFilterType::Pointer labelStatisticsFilter; labelStatisticsFilter = LabelStatisticsFilterType::New(); labelStatisticsFilter->SetInput( inputImage ); labelStatisticsFilter->SetLabelInput( hotspotMaskITK ); labelStatisticsFilter->SetCoordinateTolerance( 0.001 ); labelStatisticsFilter->SetDirectionTolerance( 0.001 ); labelStatisticsFilter->Update(); Statistics hotspotStatistics; hotspotStatistics.SetHotspotIndex(convolutionImageInformation.MaxIndex); hotspotStatistics.SetMean(convolutionImageInformation.Max); if ( labelStatisticsFilter->HasLabel( 1 ) ) { hotspotStatistics.SetLabel (1); hotspotStatistics.SetN(labelStatisticsFilter->GetCount(1)); hotspotStatistics.SetMin(labelStatisticsFilter->GetMinimum(1)); hotspotStatistics.SetMax(labelStatisticsFilter->GetMaximum(1)); hotspotStatistics.SetMedian(labelStatisticsFilter->GetMedian(1)); hotspotStatistics.SetVariance(labelStatisticsFilter->GetVariance(1)); hotspotStatistics.SetSigma(labelStatisticsFilter->GetSigma(1)); hotspotStatistics.SetRMS(sqrt( hotspotStatistics.GetMean() * hotspotStatistics.GetMean() + hotspotStatistics.GetSigma() * hotspotStatistics.GetSigma() )); MITK_DEBUG << "Statistics for inside hotspot: Mean " << hotspotStatistics.GetMean() << ", SD " << hotspotStatistics.GetSigma() << ", Max " << hotspotStatistics.GetMax() << ", Min " << hotspotStatistics.GetMin(); } else { MITK_ERROR << "Uh oh! Unable to calculate statistics for hotspot region..."; return m_EmptyStatistics; } return hotspotStatistics; } } template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator::InternalCalculateMaskFromPlanarFigure( const itk::Image< TPixel, VImageDimension > *image, unsigned int axis ) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::CastImageFilter< ImageType, MaskImage2DType > CastFilterType; // Generate mask image as new image with same header as input image and // initialize with 1. typename CastFilterType::Pointer castFilter = CastFilterType::New(); castFilter->SetInput( image ); castFilter->Update(); castFilter->GetOutput()->FillBuffer( 1 ); // all PolylinePoints of the PlanarFigure are stored in a vtkPoints object. // These points are used by the vtkLassoStencilSource to create // a vtkImageStencil. const mitk::PlaneGeometry *planarFigurePlaneGeometry = m_PlanarFigure->GetPlaneGeometry(); const typename PlanarFigure::PolyLineType planarFigurePolyline = m_PlanarFigure->GetPolyLine( 0 ); const mitk::BaseGeometry *imageGeometry3D = m_Image->GetGeometry( 0 ); // If there is a second poly line in a closed planar figure, treat it as a hole. PlanarFigure::PolyLineType planarFigureHolePolyline; if (m_PlanarFigure->GetPolyLinesSize() == 2) planarFigureHolePolyline = m_PlanarFigure->GetPolyLine(1); // Determine x- and y-dimensions depending on principal axis // TODO use plane geometry normal to determine that automatically, then check whether the PF is aligned with one of the three principal axis int i0, i1; switch ( axis ) { case 0: i0 = 1; i1 = 2; break; case 1: i0 = 0; i1 = 2; break; case 2: default: i0 = 0; i1 = 1; break; } m_PlanarFigureCoordinate0= i0; m_PlanarFigureCoordinate1= i1; // store the polyline contour as vtkPoints object bool outOfBounds = false; vtkSmartPointer points = vtkSmartPointer::New(); typename PlanarFigure::PolyLineType::const_iterator it; for ( it = planarFigurePolyline.begin(); it != planarFigurePolyline.end(); ++it ) { Point3D point3D; // Convert 2D point back to the local index coordinates of the selected // image // Fabian: From PlaneGeometry documentation: // Converts a 2D point given in mm (pt2d_mm) relative to the upper-left corner of the geometry into the corresponding world-coordinate (a 3D point in mm, pt3d_mm). // To convert a 2D point given in units (e.g., pixels in case of an image) into a 2D point given in mm (as required by this method), use IndexToWorld. planarFigurePlaneGeometry->Map( *it, point3D ); // Polygons (partially) outside of the image bounds can not be processed // further due to a bug in vtkPolyDataToImageStencil if ( !imageGeometry3D->IsInside( point3D ) ) { outOfBounds = true; } // Fabian: Why convert to index coordinates? imageGeometry3D->WorldToIndex( point3D, point3D ); points->InsertNextPoint( point3D[i0], point3D[i1], 0 ); } vtkSmartPointer holePoints = nullptr; if (!planarFigureHolePolyline.empty()) { holePoints = vtkSmartPointer::New(); Point3D point3D; PlanarFigure::PolyLineType::const_iterator end = planarFigureHolePolyline.end(); for (it = planarFigureHolePolyline.begin(); it != end; ++it) { // Fabian: same as above planarFigurePlaneGeometry->Map(*it, point3D); imageGeometry3D->WorldToIndex(point3D, point3D); holePoints->InsertNextPoint(point3D[i0], point3D[i1], 0); } } // mark a malformed 2D planar figure ( i.e. area = 0 ) as out of bounds // this can happen when all control points of a rectangle lie on the same line = two of the three extents are zero double bounds[6] = {0, 0, 0, 0, 0, 0}; points->GetBounds( bounds ); bool extent_x = (fabs(bounds[0] - bounds[1])) < mitk::eps; bool extent_y = (fabs(bounds[2] - bounds[3])) < mitk::eps; bool extent_z = (fabs(bounds[4] - bounds[5])) < mitk::eps; // throw an exception if a closed planar figure is deformed, i.e. has only one non-zero extent if ( m_PlanarFigure->IsClosed() && ((extent_x && extent_y) || (extent_x && extent_z) || (extent_y && extent_z))) { mitkThrow() << "Figure has a zero area and cannot be used for masking."; } if ( outOfBounds ) { throw std::runtime_error( "Figure at least partially outside of image bounds!" ); } // create a vtkLassoStencilSource and set the points of the Polygon vtkSmartPointer lassoStencil = vtkSmartPointer::New(); lassoStencil->SetShapeToPolygon(); lassoStencil->SetPoints( points ); vtkSmartPointer holeLassoStencil = nullptr; if (holePoints.GetPointer() != nullptr) { holeLassoStencil = vtkSmartPointer::New(); holeLassoStencil->SetShapeToPolygon(); holeLassoStencil->SetPoints(holePoints); } // Export from ITK to VTK (to use a VTK filter) typedef itk::VTKImageImport< MaskImage2DType > ImageImportType; typedef itk::VTKImageExport< MaskImage2DType > ImageExportType; typename ImageExportType::Pointer itkExporter = ImageExportType::New(); itkExporter->SetInput( castFilter->GetOutput() ); vtkSmartPointer vtkImporter = vtkSmartPointer::New(); this->ConnectPipelines( itkExporter, vtkImporter ); // Apply the generated image stencil to the input image vtkSmartPointer imageStencilFilter = vtkSmartPointer::New(); imageStencilFilter->SetInputConnection( vtkImporter->GetOutputPort() ); imageStencilFilter->SetStencilConnection(lassoStencil->GetOutputPort()); imageStencilFilter->ReverseStencilOff(); imageStencilFilter->SetBackgroundValue( 0 ); imageStencilFilter->Update(); vtkSmartPointer holeStencilFilter = nullptr; if (holeLassoStencil.GetPointer() != nullptr) { holeStencilFilter = vtkSmartPointer::New(); holeStencilFilter->SetInputConnection(imageStencilFilter->GetOutputPort()); holeStencilFilter->SetStencilConnection(holeLassoStencil->GetOutputPort()); holeStencilFilter->ReverseStencilOn(); holeStencilFilter->SetBackgroundValue(0); holeStencilFilter->Update(); } // Export from VTK back to ITK vtkSmartPointer vtkExporter = vtkSmartPointer::New(); vtkExporter->SetInputConnection( holeStencilFilter.GetPointer() == nullptr ? imageStencilFilter->GetOutputPort() : holeStencilFilter->GetOutputPort()); vtkExporter->Update(); typename ImageImportType::Pointer itkImporter = ImageImportType::New(); this->ConnectPipelines( vtkExporter, itkImporter ); itkImporter->Update(); typedef itk::ImageDuplicator< ImageImportType::OutputImageType > DuplicatorType; DuplicatorType::Pointer duplicator = DuplicatorType::New(); duplicator->SetInputImage( itkImporter->GetOutput() ); duplicator->Update(); // Store mask m_InternalImageMask2D = duplicator->GetOutput(); } void ImageStatisticsCalculator::UnmaskedStatisticsProgressUpdate() { // Need to throw away every second progress event to reach a final count of // 100 since two consecutive filters are used in this case static int updateCounter = 0; if ( updateCounter++ % 2 == 0 ) { this->InvokeEvent( itk::ProgressEvent() ); } } void ImageStatisticsCalculator::MaskedStatisticsProgressUpdate() { this->InvokeEvent( itk::ProgressEvent() ); } } diff --git a/Modules/ImageStatistics/mitkImageStatisticsCalculator2.cpp b/Modules/ImageStatistics/mitkImageStatisticsCalculator2.cpp index 24dcf08723..1c018701aa 100644 --- a/Modules/ImageStatistics/mitkImageStatisticsCalculator2.cpp +++ b/Modules/ImageStatistics/mitkImageStatisticsCalculator2.cpp @@ -1,519 +1,454 @@ #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_v2::SetInputImage(mitk::Image::Pointer image) { if (image != m_Image) { m_Image = image; m_StatisticsUpdateRequiredByTimeStep.resize(m_Image->GetTimeSteps()); m_StatisticsByTimeStep.resize(m_Image->GetTimeSteps()); this->SetAllStatisticsToUpdateRequired(); } } void ImageStatisticsCalculator_v2::SetMask(mitk::MaskGenerator::Pointer mask) { if (mask != m_MaskGenerator) { m_MaskGenerator = mask; this->SetAllStatisticsToUpdateRequired(); } } void ImageStatisticsCalculator_v2::SetNBinsForHistogramStatistics(unsigned int nBins) { if (nBins != m_nBinsForHistogramStatistics) { m_nBinsForHistogramStatistics = nBins; this->SetAllStatisticsToUpdateRequired(); } } unsigned int ImageStatisticsCalculator_v2::GetNBinsForHistogramStatistics() const { return m_nBinsForHistogramStatistics; } ImageStatisticsCalculator_v2::StatisticsContainer::Pointer ImageStatisticsCalculator_v2::GetStatistics(unsigned int timeStep, unsigned int label) { std::cout << "getting statistics..." << std::endl; if (timeStep >= m_StatisticsByTimeStep.size()) { throw std::runtime_error("invalid timeStep in ImageStatisticsCalculator_v2::GetStatistics"); } if (m_Image.IsNull()) { throw std::runtime_error("no image"); } if (m_MaskGenerator.IsNotNull()) { m_MaskGenerator->SetTimeStep(timeStep); m_InternalMask = m_MaskGenerator->GetMask(); } ImageTimeSelector::Pointer imgTimeSel = ImageTimeSelector::New(); imgTimeSel->SetInput(m_Image); imgTimeSel->SetTimeNr(timeStep); imgTimeSel->UpdateLargestPossibleRegion(); m_ImageTimeSlice = imgTimeSel->GetOutput(); if (m_StatisticsUpdateRequiredByTimeStep[timeStep]) { std::cout << "calculating statistics..." << std::endl; // Calculate statistics with/without mask if (m_MaskGenerator.IsNull()) { // 1) calculate statistics unmasked: // plug image into itkstatisticsimagefilter (will be replaced by my awesome filter later on) // retrieve statistics and save them std::cout << "unmasked" << std::endl; AccessByItk_1(m_ImageTimeSlice, InternalCalculateStatisticsUnmasked, timeStep) } else { std::cout << "masked" << std::endl; // 2) calculate statistics masked // extract mask image region // plug mask and image into itklabelstatisticsimagefilter AccessByItk_1(m_ImageTimeSlice, InternalCalculateStatisticsMasked, timeStep) } m_StatisticsUpdateRequiredByTimeStep[timeStep]=false; } for (std::vector::iterator it = m_StatisticsByTimeStep[timeStep].begin(); it != m_StatisticsByTimeStep[timeStep].end(); ++it) { StatisticsContainer::Pointer statCont = *it; if (statCont->GetLabel() == label) { return statCont; } } MITK_ERROR << "Invalid Label: " << label; } void ImageStatisticsCalculator_v2::SetAllStatisticsToUpdateRequired() { for (unsigned int i = 0; i < m_StatisticsUpdateRequiredByTimeStep.size(); i++) { this->SetStatsTimeStepToUpdateRequired(i); } } void ImageStatisticsCalculator_v2::SetStatsTimeStepToUpdateRequired(unsigned int timeStep) { if (timeStep >= m_StatisticsUpdateRequiredByTimeStep.size()) { throw std::runtime_error("invalid timeStep in ImageStatisticsCalculator_v2::SetStatsTimeStepToUpdateRequired"); } m_StatisticsUpdateRequiredByTimeStep[timeStep] = true; m_StatisticsByTimeStep[timeStep].resize(0); } template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator_v2::InternalCalculateStatisticsUnmasked( typename itk::Image< TPixel, VImageDimension >* image, unsigned int timeStep) { typedef typename itk::Image< TPixel, VImageDimension > ImageType; + typedef typename itk::ExtendedStatisticsImageFilter2 ImageStatisticsFilterType; + typedef typename itk::MinMaxImageFilterWithIndex MinMaxFilterType; StatisticsContainer::Pointer statisticsResult = StatisticsContainer::New(); - typename itk::ExtendedStatisticsImageFilter::Pointer statisticsFilter = itk::ExtendedStatisticsImageFilter::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(); +// typename itk::MinimumMaximumImageCalculator::Pointer imgMinMaxFilter = itk::MinimumMaximumImageCalculator::New(); +// imgMinMaxFilter->SetImage(image); +// imgMinMaxFilter->Compute(); vnl_vector minIndex, maxIndex; - typename ImageType::IndexType tmpMinIndex = imgMinMaxFilter->GetIndexOfMinimum(); - typename ImageType::IndexType tmpMaxIndex = imgMinMaxFilter->GetIndexOfMaximum(); + 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); - // TODO: REMOVE THIS, the histogram statistics should be a separate module - statisticsFilter->SetBinSize(m_nBinsForHistogramStatistics); + //statisticsFilter->SetHistogramParameters(m_nBinsForHistogramStatistics, imgMinMaxFilter->GetMinimum(), imgMinMaxFilter->GetMaximum()); + statisticsFilter->SetHistogramParameters(m_nBinsForHistogramStatistics, 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()); - - // histogram statistics - HistogramType::Pointer histogram = InternalCalculateHistogramUnmasked(image, statisticsResult->GetMin(), statisticsResult->GetMax()); - statisticsResult->SetHistogram(histogram); - - HistogramStatisticsCalculator histStatCalc; - histStatCalc.SetHistogram(histogram); - histStatCalc.CalculateStatistics(); - - statisticsResult->SetEntropy(histStatCalc.GetEntropy()); - statisticsResult->SetMedian(histStatCalc.GetMedian()); - statisticsResult->SetUniformity(histStatCalc.GetUniformity()); - statisticsResult->SetUPP(histStatCalc.GetUPP()); + 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 > ImageStatisticsCalculator_v2::HistogramType::Pointer ImageStatisticsCalculator_v2::InternalCalculateHistogramUnmasked( - typename itk::Image< TPixel, VImageDimension >* image, - double minVal, - double maxVal) - { - typedef typename itk::Image< TPixel, VImageDimension > ImageType; - typedef itk::Statistics::ImageToHistogramFilter HistogramFilterType; - typedef itk::Statistics::Histogram HistogramType; - typedef typename HistogramFilterType::HistogramSizeType SizeType; - - typename HistogramFilterType::Pointer histogramFilter = HistogramFilterType::New(); - histogramFilter->SetInput(image); - SizeType histSize(1); - histSize[0] = m_nBinsForHistogramStatistics; - histogramFilter->SetHistogramSize(histSize); // is that correct?? - typename HistogramFilterType::HistogramMeasurementVectorType lowerBound(1); - typename HistogramFilterType::HistogramMeasurementVectorType upperBound(1); - lowerBound[0] = minVal; - upperBound[0] = maxVal; - histogramFilter->SetHistogramBinMinimum(lowerBound); - histogramFilter->SetHistogramBinMaximum(upperBound); - - histogramFilter->Update(); - typename HistogramType::Pointer histogram = histogramFilter->GetOutput(); - - return histogram; - } - - - - /* the whole min and max calculation thing is completely bad... need a new itk filter for that */ -// template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator_v2::GetMinAndMaxValue( -// typename itk::Image< TPixel, VImageDimension >* inputImage, -// double &minVal, -// double &maxVal -// ) -// { -// typedef itk::Image< TPixel, VImageDimension > ImageType; - -// typedef itk::ImageRegionConstIteratorWithIndex ImageRegionIteratorType; - -// ImageRegionIteratorType imageIterator( inputImage, inputImage->GetLargestPossibleRegion() ); - -// maxVal = std::numeric_limits::min(); -// minVal = std::numeric_limits::max(); -// long counter(0); -// double actualPixelValue(0); - -// for( imageIterator.GoToBegin(); !imageIterator.IsAtEnd(); ++imageIterator) -// { -// counter++; - -// actualPixelValue = imageIterator.Value(); - -// if(actualPixelValue > maxVal) -// { -// maxVal = actualPixelValue; -// } -// else if(actualPixelValue < minVal) -// { -// minVal = actualPixelValue; -// } -// } -// } - template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator_v2::InternalCalculateStatisticsMasked( typename itk::Image< TPixel, VImageDimension >* image, unsigned int timeStep) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< unsigned short, VImageDimension > MaskType; - typedef itk::ExtendedLabelStatisticsImageFilter< ImageType, MaskType > ImageStatisticsFilterType; + typedef typename MaskType::PixelType LabelPixelType; + typedef itk::ExtendedLabelStatisticsImageFilter2< ImageType, MaskType > ImageStatisticsFilterType; typedef MaskUtilities< TPixel, VImageDimension > MaskUtilType; - typedef itk::MaskImageFilter2 MaskImageFilterType2; - typedef itk::MaskImageFilter MaskImageFilterType; - typedef itk::MinimumMaximumImageCalculator MinMaxCalculatorType; + typedef typename itk::MinMaxLabelImageFilterWithIndex MinMaxLabelFilterType; + typedef typename ImageType::PixelType InputImgPixelType; // maskImage has to have the same dimension as image typename MaskType::Pointer maskImage = MaskType::New(); maskImage = ImageToItkImage< unsigned short, VImageDimension >(m_InternalMask); typename MaskUtilType::Pointer maskUtil = MaskUtilType::New(); maskUtil->SetImage(image); maskUtil->SetMask(maskImage.GetPointer()); // if mask is smaller than image, extract the image region where the mask is typename ImageType::Pointer adaptedImage = ImageType::New(); adaptedImage = maskUtil->ExtractMaskImageRegion(); // this also checks mask sanity -// double minVal, maxVal; -// GetMinAndMaxValue(adaptedImage.GetPointer(), minVal, maxVal); +// typename itk::ImageFileWriter::Pointer imgFileWriter = itk::ImageFileWriter::New(); +// imgFileWriter->SetInput(adaptedImage); +// imgFileWriter->SetFileName("/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/ct_leber_adapted_image.nrrd"); +// imgFileWriter->Update(); + +// typename itk::ImageFileWriter::Pointer imgFileWriter2 = itk::ImageFileWriter::New(); +// imgFileWriter2->SetInput(maskImage); +// imgFileWriter2->SetFileName("/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/ct_leber_adapted_image_mask.nrrd"); +// imgFileWriter2->Update(); - /** Searching the min and max value is only done over the pixels that are nonzero in the mask. This bevahior differs from the 'old' - * imagestatisticscalculator where min and max were looked for over the entire adaptedImage*/ // find min, max, minindex and maxindex - // make sure to only look in the masked region, use a masker for this - typename MaskImageFilterType::Pointer maskImageFilter = MaskImageFilterType::New(); - maskImageFilter->SetInput1(adaptedImage); - maskImageFilter->SetInput2(maskImage); - maskImageFilter->SetCoordinateTolerance(0.001); - maskImageFilter->SetDirectionTolerance(0.001); - maskImageFilter->Update(); - - typename MinMaxCalculatorType::Pointer imgMinMaxFilter = MinMaxCalculatorType::New(); - imgMinMaxFilter->SetImage(maskImageFilter->GetOutput()); - imgMinMaxFilter->Compute(); - double minVal, maxVal; - minVal = imgMinMaxFilter->GetMinimum(); - maxVal = imgMinMaxFilter->GetMaximum(); + typename MinMaxLabelFilterType::Pointer minMaxFilter = MinMaxLabelFilterType::New(); + minMaxFilter->SetInput(adaptedImage); + minMaxFilter->SetLabelInput(maskImage); + minMaxFilter->UpdateLargestPossibleRegion(); + + // set histogram parameters for each label individually + 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))); + nBins.insert(typename std::pair(label, m_nBinsForHistogramStatistics)); + } + +// minVal = minMaxFilter->GetGlobalMin(); +// maxVal = minMaxFilter->GetGlobalMax(); typename ImageStatisticsFilterType::Pointer imageStatisticsFilter = ImageStatisticsFilterType::New(); imageStatisticsFilter->SetDirectionTolerance(0.001); imageStatisticsFilter->SetCoordinateTolerance(0.001); imageStatisticsFilter->SetInput(adaptedImage); imageStatisticsFilter->SetLabelInput(maskImage); - imageStatisticsFilter->SetHistogramParameters(m_nBinsForHistogramStatistics, floor(minVal), ceil(maxVal)); + //imageStatisticsFilter->SetHistogramParameters(m_nBinsForHistogramStatistics, floor(minVal), ceil(maxVal)); + 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 - typename MaskImageFilterType2::Pointer maskImageFilter = MaskImageFilterType2::New(); - maskImageFilter->SetInput1(adaptedImage); - maskImageFilter->SetInput2(maskImage); - maskImageFilter->SetOutsideValue(0); - maskImageFilter->SetMaskingValue(*it); - maskImageFilter->SetCoordinateTolerance(0.001); - maskImageFilter->SetDirectionTolerance(0.001); - maskImageFilter->UpdateLargestPossibleRegion(); - -// std::string outfile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic3D_masked_image.nrrd"; -// typedef itk::ImageFileWriter< ImageType > WriterType; -// typename WriterType::Pointer writer = WriterType::New(); -// writer->SetFileName(outfile); -// writer->SetInput(maskImageFilter->GetOutput()); -// writer->Update(); - - // TODO adapt minindex and maxindex if mask was 2D (careful: a mask can be 2d even if the itk::Image is 3D (3d image with a size of 1 in one dimension)) - typename MinMaxCalculatorType::Pointer imgMinMaxFilter = MinMaxCalculatorType::New(); - imgMinMaxFilter->SetImage(maskImageFilter->GetOutput()); - imgMinMaxFilter->Compute(); + vnl_vector minIndex, maxIndex; - typename ImageType::IndexType tmpMinIndex = imgMinMaxFilter->GetIndexOfMinimum(); - typename ImageType::IndexType tmpMaxIndex = imgMinMaxFilter->GetIndexOfMaximum(); + 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(minMaxFilter->GetMax(*it) == imageStatisticsFilter->GetMaximum(*it)); + assert(minMaxFilter->GetMin(*it) == imageStatisticsFilter->GetMinimum(*it)); + + statisticsResult->SetN(imageStatisticsFilter->GetSum(*it) / (double) imageStatisticsFilter->GetMean(*it)); statisticsResult->SetMean(imageStatisticsFilter->GetMean(*it)); statisticsResult->SetMin(imageStatisticsFilter->GetMinimum(*it)); statisticsResult->SetMax(imageStatisticsFilter->GetMaximum(*it)); statisticsResult->SetVariance(imageStatisticsFilter->GetVariance(*it)); statisticsResult->SetStd(imageStatisticsFilter->GetSigma(*it)); statisticsResult->SetSkewness(imageStatisticsFilter->GetSkewness(*it)); statisticsResult->SetKurtosis(imageStatisticsFilter->GetKurtosis(*it)); statisticsResult->SetRMS(std::sqrt(std::pow(imageStatisticsFilter->GetMean(*it), 2.) + imageStatisticsFilter->GetVariance(*it))); // variance = sigma^2 statisticsResult->SetMPP(imageStatisticsFilter->GetMPP(*it)); statisticsResult->SetLabel(*it); - HistogramType::Pointer histogram = imageStatisticsFilter->GetHistogram(*it); - statisticsResult->SetHistogram(histogram); - - HistogramStatisticsCalculator histStatCalc; - histStatCalc.SetHistogram(histogram); - histStatCalc.CalculateStatistics(); - - statisticsResult->SetEntropy(histStatCalc.GetEntropy()); - statisticsResult->SetMedian(histStatCalc.GetMedian()); - statisticsResult->SetUniformity(histStatCalc.GetUniformity()); - statisticsResult->SetUPP(histStatCalc.GetUPP()); + 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; } } ImageStatisticsCalculator_v2::StatisticsContainer::StatisticsContainer(): m_N(std::numeric_limits::max()), m_Mean(std::numeric_limits::max()), m_Min(std::numeric_limits::min()), m_Max(std::numeric_limits::max()), m_Std(std::numeric_limits::max()), m_Variance(std::numeric_limits::max()), m_Skewness(std::numeric_limits::max()), m_Kurtosis(std::numeric_limits::max()), m_RMS(std::numeric_limits::max()), m_MPP(std::numeric_limits::max()), m_Median(std::numeric_limits::max()), m_Uniformity(std::numeric_limits::max()), m_UPP(std::numeric_limits::max()), m_Entropy(std::numeric_limits::max()) { m_minIndex.set_size(0); m_maxIndex.set_size(0); } ImageStatisticsCalculator_v2::statisticsMapType ImageStatisticsCalculator_v2::StatisticsContainer::GetStatisticsAsMap() { ImageStatisticsCalculator_v2::statisticsMapType statisticsAsMap; statisticsAsMap["N"] = m_N; statisticsAsMap["Mean"] = m_Mean; statisticsAsMap["Min"] = m_Min; statisticsAsMap["Max"] = m_Max; statisticsAsMap["StandardDeviation"] = m_Std; statisticsAsMap["Variance"] = m_Variance; statisticsAsMap["Skewness"] = m_Skewness; statisticsAsMap["Kurtosis"] = m_Kurtosis; statisticsAsMap["RMS"] = m_RMS; statisticsAsMap["MPP"] = m_MPP; statisticsAsMap["Median"] = m_Median; statisticsAsMap["Uniformity"] = m_Uniformity; statisticsAsMap["UPP"] = m_UPP; statisticsAsMap["Entropy"] = m_Entropy; return statisticsAsMap; } void ImageStatisticsCalculator_v2::StatisticsContainer::Reset() { m_N = std::numeric_limits::max(); m_Mean = std::numeric_limits::max(); m_Min = std::numeric_limits::max(); m_Max = std::numeric_limits::max(); m_Std = std::numeric_limits::max(); m_Variance = std::numeric_limits::max(); m_Skewness = std::numeric_limits::max(); m_Kurtosis = std::numeric_limits::max(); m_RMS = std::numeric_limits::max(); m_MPP = std::numeric_limits::max(); m_Median = std::numeric_limits::max(); m_Uniformity = std::numeric_limits::max(); m_UPP = std::numeric_limits::max(); m_Entropy = std::numeric_limits::max(); } void ImageStatisticsCalculator_v2::StatisticsContainer::PrintSelf() { ImageStatisticsCalculator_v2::statisticsMapType statMap = this->GetStatisticsAsMap(); // print all map key value pairs // const auto& val:statMap for (auto it = statMap.begin(); it != statMap.end(); ++it) { std::cout << it->first << ": " << it->second << std::endl; } // print the min and max index std::cout << "Min Index:" << std::endl; for (auto it = this->GetMinIndex().begin(); it != this->GetMinIndex().end(); it++) { std::cout << *it << " "; } std::cout << std::endl; // print the min and max index std::cout << "Max Index:" << std::endl; for (auto it = this->GetMaxIndex().begin(); it != this->GetMaxIndex().end(); it++) { std::cout << *it << " "; } std::cout << std::endl; } } diff --git a/Modules/ImageStatistics/mitkMinMaxImageFilterWithIndex.h b/Modules/ImageStatistics/mitkMinMaxImageFilterWithIndex.h new file mode 100644 index 0000000000..b5768f4cd4 --- /dev/null +++ b/Modules/ImageStatistics/mitkMinMaxImageFilterWithIndex.h @@ -0,0 +1,84 @@ +#ifndef MITK_MINMAXIMAGEFILTERWITHINDEX_H +#define MITK_MINMAXIMAGEFILTERWITHINDEX_H + +#include + +#include +#include +#include + + + +namespace itk +{ +template +class MinMaxImageFilterWithIndex: public itk::ImageToImageFilter +{ +public: + /** Standard Self typedef */ + typedef MinMaxImageFilterWithIndex 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(MinMaxImageFilterWithIndex, 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; + + + RealType GetMin() const + { + return m_Min; + } + + RealType GetMax() const + { + return m_Max; + } + + IndexType GetMinIndex() const + { + return m_MinIndex; + } + + IndexType GetMaxIndex() const + { + return m_MaxIndex; + } + +protected: + void AllocateOutputs(); + + void ThreadedGenerateData(const RegionType & + outputRegionForThread, + ThreadIdType threadId); + + void BeforeThreadedGenerateData(); + + void AfterThreadedGenerateData(); + +private: + std::vector m_ThreadMin; + std::vector m_ThreadMax; + std::vector m_ThreadMinIndex; + std::vector m_ThreadMaxIndex; + + PixelType m_Min; + PixelType m_Max; + IndexType m_MinIndex; + IndexType m_MaxIndex; +}; +} + +#include "mitkMinMaxImageFilterWithIndex.hxx" + + +#endif diff --git a/Modules/ImageStatistics/mitkMinMaxImageFilterWithIndex.hxx b/Modules/ImageStatistics/mitkMinMaxImageFilterWithIndex.hxx new file mode 100644 index 0000000000..c925cb5aa4 --- /dev/null +++ b/Modules/ImageStatistics/mitkMinMaxImageFilterWithIndex.hxx @@ -0,0 +1,106 @@ +#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; + + 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/mitkMinMaxLabelmageFilterWithIndex.h b/Modules/ImageStatistics/mitkMinMaxLabelmageFilterWithIndex.h new file mode 100644 index 0000000000..5ae6208965 --- /dev/null +++ b/Modules/ImageStatistics/mitkMinMaxLabelmageFilterWithIndex.h @@ -0,0 +1,171 @@ +#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; + + 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; + } + + 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/mitkMinMaxLabelmageFilterWithIndex.hxx b/Modules/ImageStatistics/mitkMinMaxLabelmageFilterWithIndex.hxx new file mode 100644 index 0000000000..c1a0af7ae2 --- /dev/null +++ b/Modules/ImageStatistics/mitkMinMaxLabelmageFilterWithIndex.hxx @@ -0,0 +1,130 @@ +#ifndef MITK_MinMaxLabelImageFilterWithIndex_HXX +#define MITK_MinMaxLabelImageFilterWithIndex_HXX + +#include +#include + +namespace itk +{ + + +template< typename TInputImage, typename TLabelImage > +void MinMaxLabelImageFilterWithIndex< TInputImage, TLabelImage >::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, typename TLabelImage > +void MinMaxLabelImageFilterWithIndex< TInputImage, TLabelImage >::ThreadedGenerateData(const RegionType & + outputRegionForThread, + ThreadIdType threadId) +{ + const SizeValueType size0 = outputRegionForThread.GetSize(0); + if( size0 == 0) + { + return; + } + PixelType value; + LabelPixelType label; + + ExtremaMapType threadExtrema; + ExtremaMapTypeIterator threadExtremaIt; + + ImageRegionConstIteratorWithIndex< TInputImage > it (this->GetInput(), outputRegionForThread); + ImageRegionConstIteratorWithIndex< TLabelImage > labelit (this->GetLabelInput(), outputRegionForThread); + + // do the work + while ( !it.IsAtEnd() ) + { + value = it.Get(); + label = labelit.Get(); + + threadExtremaIt = threadExtrema.find(label); + + // if label does not exist yet, create a new entry in the map. + if (threadExtremaIt == threadExtrema.end()) + { + threadExtremaIt = threadExtrema.insert( MapValueType(label, LabelExtrema()) ).first; + } + + if (value < (*threadExtremaIt).second.m_Min) + { + (*threadExtremaIt).second.m_Min = value; + (*threadExtremaIt).second.m_MinIndex = it.GetIndex(); + } + if (value > (*threadExtremaIt).second.m_Max) + { + (*threadExtremaIt).second.m_Max = value; + (*threadExtremaIt).second.m_MaxIndex = it.GetIndex(); + } + ++it; + ++labelit; + } + + m_ThreadExtrema[threadId] = threadExtrema; +} + +template< typename TInputImage, typename TLabelImage > +void MinMaxLabelImageFilterWithIndex< TInputImage, TLabelImage >::BeforeThreadedGenerateData() +{ + ThreadIdType numberOfThreads = this->GetNumberOfThreads(); + m_ThreadExtrema.resize(numberOfThreads); + + for (unsigned int i =0; i < numberOfThreads; i++) + { + m_ThreadExtrema[i] = ExtremaMapType(); + } +} + +template< typename TInputImage, typename TLabelImage > +void MinMaxLabelImageFilterWithIndex< TInputImage, TLabelImage >::AfterThreadedGenerateData() +{ + ThreadIdType numberOfThreads = this->GetNumberOfThreads(); + + m_GlobalMin = std::numeric_limits::max(); + m_GlobalMax = std::numeric_limits::min(); + + ExtremaMapTypeIterator it; + + for (ThreadIdType i = 0; i < numberOfThreads; i++) + { + for (auto&& it2 : m_ThreadExtrema[i]) + { + it = m_LabelExtrema.find(it2.first); + if (it == m_LabelExtrema.end()) + { + it = m_LabelExtrema.insert( MapValueType(it2.first, LabelExtrema()) ).first; + } + + if (it2.second.m_Min < (*it).second.m_Min) + { + (*it).second.m_Min = it2.second.m_Min; + (*it).second.m_MinIndex = it2.second.m_MinIndex; + if (it2.second.m_Min < m_GlobalMin) + { + m_GlobalMin = it2.second.m_Min; + m_GlobalMinIndex = it2.second.m_MinIndex; + } + } + + if (it2.second.m_Max > (*it).second.m_Max) + { + (*it).second.m_Max = it2.second.m_Max; + (*it).second.m_MaxIndex = it2.second.m_MaxIndex; + if (it2.second.m_Max > m_GlobalMax) + { + m_GlobalMax = it2.second.m_Max; + m_GlobalMaxIndex = it2.second.m_MaxIndex; + } + } + } + } +} +} +#endif