diff --git a/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp b/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp index 723e7702db..65c9962af9 100644 --- a/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp +++ b/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp @@ -1,552 +1,552 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #include "mitkImageStatisticsCalculator.h" #include #include #include #include #include #include #include #include #include #include #include #include namespace mitk { void ImageStatisticsCalculator::SetInputImage(const mitk::Image *image) { if (image != m_Image) { m_Image = image; this->Modified(); } } void ImageStatisticsCalculator::SetMask(mitk::MaskGenerator *mask) { if (mask != m_MaskGenerator) { m_MaskGenerator = mask; this->Modified(); } } void ImageStatisticsCalculator::SetSecondaryMask(mitk::MaskGenerator *mask) { if (mask != m_SecondaryMaskGenerator) { m_SecondaryMaskGenerator = mask; this->Modified(); } } void ImageStatisticsCalculator::SetNBinsForHistogramStatistics(unsigned int nBins) { if (nBins != m_nBinsForHistogramStatistics) { m_nBinsForHistogramStatistics = nBins; this->Modified(); this->m_UseBinSizeOverNBins = false; } if (m_UseBinSizeOverNBins) { this->Modified(); this->m_UseBinSizeOverNBins = false; } } unsigned int ImageStatisticsCalculator::GetNBinsForHistogramStatistics() const { return m_nBinsForHistogramStatistics; } void ImageStatisticsCalculator::SetBinSizeForHistogramStatistics(double binSize) { if (binSize != m_binSizeForHistogramStatistics) { m_binSizeForHistogramStatistics = binSize; this->Modified(); this->m_UseBinSizeOverNBins = true; } if (!m_UseBinSizeOverNBins) { this->Modified(); this->m_UseBinSizeOverNBins = true; } } double ImageStatisticsCalculator::GetBinSizeForHistogramStatistics() const { return m_binSizeForHistogramStatistics; } mitk::ImageStatisticsContainer* ImageStatisticsCalculator::GetStatistics(LabelIndex label) { if (m_Image.IsNull()) { mitkThrow() << "no image"; } if (!m_Image->IsInitialized()) { mitkThrow() << "Image not initialized!"; } if (IsUpdateRequired(label)) { auto timeGeometry = m_Image->GetTimeGeometry(); // always compute statistics on all timesteps for (unsigned int timeStep = 0; timeStep < m_Image->GetTimeSteps(); timeStep++) { if (m_MaskGenerator.IsNotNull()) { m_MaskGenerator->SetTimeStep(timeStep); //See T25625: otherwise, the mask is not computed again after setting a different time step m_MaskGenerator->Modified(); m_InternalMask = m_MaskGenerator->GetMask(); if (m_MaskGenerator->GetReferenceImage().IsNotNull()) { m_InternalImageForStatistics = m_MaskGenerator->GetReferenceImage(); } else { m_InternalImageForStatistics = m_Image; } } else { m_InternalImageForStatistics = m_Image; } if (m_SecondaryMaskGenerator.IsNotNull()) { m_SecondaryMaskGenerator->SetTimeStep(timeStep); m_SecondaryMask = m_SecondaryMaskGenerator->GetMask(); } ImageTimeSelector::Pointer imgTimeSel = ImageTimeSelector::New(); imgTimeSel->SetInput(m_InternalImageForStatistics); imgTimeSel->SetTimeNr(timeStep); imgTimeSel->UpdateLargestPossibleRegion(); imgTimeSel->Update(); m_ImageTimeSlice = imgTimeSel->GetOutput(); // Calculate statistics with/without mask if (m_MaskGenerator.IsNull() && m_SecondaryMaskGenerator.IsNull()) { // 1) calculate statistics unmasked: AccessByItk_2(m_ImageTimeSlice, InternalCalculateStatisticsUnmasked, timeGeometry, timeStep) } else { // 2) calculate statistics masked AccessByItk_2(m_ImageTimeSlice, InternalCalculateStatisticsMasked, timeGeometry, timeStep) } } } auto it = m_StatisticContainers.find(label); if (it != m_StatisticContainers.end()) { return (it->second).GetPointer(); } else { mitkThrow() << "unknown label"; return nullptr; } } template void ImageStatisticsCalculator::InternalCalculateStatisticsUnmasked( typename itk::Image *image, const TimeGeometry *timeGeometry, TimeStepType timeStep) { typedef typename itk::Image ImageType; typedef typename StatisticsImageFilter ImageStatisticsFilterType; typedef typename itk::MinMaxImageFilterWithIndex MinMaxFilterType; // reset statistics container if exists ImageStatisticsContainer::Pointer statisticContainerForImage; LabelIndex labelNoMask = 1; auto it = m_StatisticContainers.find(labelNoMask); if (it != m_StatisticContainers.end()) { statisticContainerForImage = it->second; } else { statisticContainerForImage = ImageStatisticsContainer::New(); statisticContainerForImage->SetTimeGeometry(const_cast(timeGeometry)); m_StatisticContainers.emplace(labelNoMask, statisticContainerForImage); } auto statObj = ImageStatisticsContainer::ImageStatisticsObject(); typename ImageStatisticsFilterType::Pointer statisticsFilter = ImageStatisticsFilterType::New(); statisticsFilter->SetInput(image); statisticsFilter->SetCoordinateTolerance(0.001); statisticsFilter->SetDirectionTolerance(0.001); // TODO: this is single threaded. Implement our own image filter that does this multi threaded // typename itk::MinimumMaximumImageCalculator::Pointer imgMinMaxFilter = // itk::MinimumMaximumImageCalculator::New(); imgMinMaxFilter->SetImage(image); // imgMinMaxFilter->Compute(); vnl_vector minIndex, maxIndex; typename MinMaxFilterType::Pointer minMaxFilter = MinMaxFilterType::New(); minMaxFilter->SetInput(image); minMaxFilter->UpdateLargestPossibleRegion(); typename ImageType::PixelType minval = minMaxFilter->GetMin(); typename ImageType::PixelType maxval = minMaxFilter->GetMax(); typename ImageType::IndexType tmpMinIndex = minMaxFilter->GetMinIndex(); typename ImageType::IndexType tmpMaxIndex = minMaxFilter->GetMaxIndex(); // typename ImageType::IndexType tmpMinIndex = imgMinMaxFilter->GetIndexOfMinimum(); // typename ImageType::IndexType tmpMaxIndex = imgMinMaxFilter->GetIndexOfMaximum(); minIndex.set_size(tmpMaxIndex.GetIndexDimension()); maxIndex.set_size(tmpMaxIndex.GetIndexDimension()); for (unsigned int i = 0; i < tmpMaxIndex.GetIndexDimension(); i++) { minIndex[i] = tmpMinIndex[i]; maxIndex[i] = tmpMaxIndex[i]; } statObj.AddStatistic(mitk::ImageStatisticsConstants::MINIMUMPOSITION(), minIndex); statObj.AddStatistic(mitk::ImageStatisticsConstants::MAXIMUMPOSITION(), maxIndex); // convert m_binSize in m_nBins if necessary unsigned int nBinsForHistogram; if (m_UseBinSizeOverNBins) { nBinsForHistogram = std::max(static_cast(std::ceil(maxval - minval)) / m_binSizeForHistogramStatistics, 10.); // do not allow less than 10 bins } else { nBinsForHistogram = m_nBinsForHistogramStatistics; } statisticsFilter->SetHistogramParameters(nBinsForHistogram, minval, maxval); try { statisticsFilter->Update(); } catch (const itk::ExceptionObject &e) { mitkThrow() << "Image statistics calculation failed due to following ITK Exception: \n " << e.what(); } auto voxelVolume = GetVoxelVolume(image); auto numberOfPixels = image->GetLargestPossibleRegion().GetNumberOfPixels(); auto volume = static_cast(numberOfPixels) * voxelVolume; auto variance = statisticsFilter->GetSigma() * statisticsFilter->GetSigma(); auto rms = std::sqrt(std::pow(statisticsFilter->GetMean(), 2.) + statisticsFilter->GetVariance()); // variance = sigma^2 statObj.AddStatistic(mitk::ImageStatisticsConstants::NUMBEROFVOXELS(), static_cast(numberOfPixels)); statObj.AddStatistic(mitk::ImageStatisticsConstants::VOLUME(), volume); statObj.AddStatistic(mitk::ImageStatisticsConstants::MEAN(), statisticsFilter->GetMean()); statObj.AddStatistic(mitk::ImageStatisticsConstants::MINIMUM(), static_cast(statisticsFilter->GetMinimum())); statObj.AddStatistic(mitk::ImageStatisticsConstants::MAXIMUM(), static_cast(statisticsFilter->GetMaximum())); statObj.AddStatistic(mitk::ImageStatisticsConstants::STANDARDDEVIATION(), statisticsFilter->GetSigma()); statObj.AddStatistic(mitk::ImageStatisticsConstants::VARIANCE(), variance); statObj.AddStatistic(mitk::ImageStatisticsConstants::SKEWNESS(), statisticsFilter->GetSkewness()); statObj.AddStatistic(mitk::ImageStatisticsConstants::KURTOSIS(), statisticsFilter->GetKurtosis()); statObj.AddStatistic(mitk::ImageStatisticsConstants::RMS(), rms); statObj.AddStatistic(mitk::ImageStatisticsConstants::MPP(), statisticsFilter->GetMPP()); statObj.AddStatistic(mitk::ImageStatisticsConstants::ENTROPY(), statisticsFilter->GetEntropy()); statObj.AddStatistic(mitk::ImageStatisticsConstants::MEDIAN(), statisticsFilter->GetMedian()); statObj.AddStatistic(mitk::ImageStatisticsConstants::UNIFORMITY(), statisticsFilter->GetUniformity()); statObj.AddStatistic(mitk::ImageStatisticsConstants::UPP(), statisticsFilter->GetUPP()); - statObj.m_Histogram = statisticsFilter->GetHistogram().GetPointer(); + statObj.m_Histogram = statisticsFilter->GetHistogram(); statisticContainerForImage->SetStatisticsForTimeStep(timeStep, statObj); } template double ImageStatisticsCalculator::GetVoxelVolume(typename itk::Image *image) const { auto spacing = image->GetSpacing(); double voxelVolume = 1.; for (unsigned int i = 0; i < image->GetImageDimension(); i++) { voxelVolume *= spacing[i]; } return voxelVolume; } template void ImageStatisticsCalculator::InternalCalculateStatisticsMasked(typename itk::Image *image, const TimeGeometry *timeGeometry, unsigned int timeStep) { typedef itk::Image ImageType; typedef itk::Image MaskType; typedef typename MaskType::PixelType LabelPixelType; typedef itk::ExtendedLabelStatisticsImageFilter ImageStatisticsFilterType; typedef MaskUtilities MaskUtilType; typedef typename itk::MinMaxLabelImageFilterWithIndex MinMaxLabelFilterType; typedef typename ImageType::PixelType InputImgPixelType; // workaround: if m_SecondaryMaskGenerator ist not null but m_MaskGenerator is! (this is the case if we request a // 'ignore zuero valued pixels' mask in the gui but do not define a primary mask) bool swapMasks = false; if (m_SecondaryMask.IsNotNull() && m_InternalMask.IsNull()) { m_InternalMask = m_SecondaryMask; m_SecondaryMask = nullptr; swapMasks = true; } // maskImage has to have the same dimension as image typename MaskType::Pointer maskImage = MaskType::New(); try { // try to access the pixel values directly (no copying or casting). Only works if mask pixels are of pixelType // unsigned short maskImage = ImageToItkImage(m_InternalMask); } catch (const itk::ExceptionObject &) { // if the pixel type of the mask is not short, then we have to make a copy of m_InternalMask (and cast the values) CastToItkImage(m_InternalMask, maskImage); } // if we have a secondary mask (say a ignoreZeroPixelMask) we need to combine the masks (corresponds to AND) if (m_SecondaryMask.IsNotNull()) { // dirty workaround for a bug when pf mask + any other mask is used in conjunction. We need a proper fix for this // (Fabian Isensee is responsible and probably working on it!) if (m_InternalMask->GetDimension() == 2 && (m_SecondaryMask->GetDimension() == 3 || m_SecondaryMask->GetDimension() == 4)) { mitk::Image::ConstPointer old_img = m_SecondaryMaskGenerator->GetReferenceImage(); m_SecondaryMaskGenerator->SetInputImage(m_MaskGenerator->GetReferenceImage()); m_SecondaryMask = m_SecondaryMaskGenerator->GetMask(); m_SecondaryMaskGenerator->SetInputImage(old_img); } typename MaskType::Pointer secondaryMaskImage = MaskType::New(); secondaryMaskImage = ImageToItkImage(m_SecondaryMask); // secondary mask should be a ignore zero value pixel mask derived from image. it has to be cropped to the mask // region (which may be planar or simply smaller) typename MaskUtilities::Pointer secondaryMaskMaskUtil = MaskUtilities::New(); secondaryMaskMaskUtil->SetImage(secondaryMaskImage.GetPointer()); secondaryMaskMaskUtil->SetMask(maskImage.GetPointer()); typename MaskType::Pointer adaptedSecondaryMaskImage = secondaryMaskMaskUtil->ExtractMaskImageRegion(); typename itk::MaskImageFilter2::Pointer maskFilter = itk::MaskImageFilter2::New(); maskFilter->SetInput1(maskImage); maskFilter->SetInput2(adaptedSecondaryMaskImage); maskFilter->SetMaskingValue( 1); // all pixels of maskImage where secondaryMaskImage==1 will be kept, all the others are set to 0 maskFilter->UpdateLargestPossibleRegion(); maskImage = maskFilter->GetOutput(); } typename MaskUtilType::Pointer maskUtil = MaskUtilType::New(); maskUtil->SetImage(image); maskUtil->SetMask(maskImage.GetPointer()); // if mask is smaller than image, extract the image region where the mask is typename ImageType::Pointer adaptedImage = ImageType::New(); adaptedImage = maskUtil->ExtractMaskImageRegion(); // this also checks mask sanity // find min, max, minindex and maxindex typename MinMaxLabelFilterType::Pointer minMaxFilter = MinMaxLabelFilterType::New(); minMaxFilter->SetInput(adaptedImage); minMaxFilter->SetLabelInput(maskImage); minMaxFilter->UpdateLargestPossibleRegion(); // set histogram parameters for each label individually (min/max may be different for each label) typedef typename std::map MapType; typedef typename std::pair PairType; std::vector relevantLabels = minMaxFilter->GetRelevantLabels(); MapType minVals; MapType maxVals; std::map nBins; for (LabelPixelType label : relevantLabels) { minVals.insert(PairType(label, minMaxFilter->GetMin(label))); maxVals.insert(PairType(label, minMaxFilter->GetMax(label))); unsigned int nBinsForHistogram; if (m_UseBinSizeOverNBins) { nBinsForHistogram = std::max(static_cast(std::ceil(minMaxFilter->GetMax(label) - minMaxFilter->GetMin(label))) / m_binSizeForHistogramStatistics, 10.); // do not allow less than 10 bins } else { nBinsForHistogram = m_nBinsForHistogramStatistics; } nBins.insert(typename std::pair(label, nBinsForHistogram)); } typename ImageStatisticsFilterType::Pointer imageStatisticsFilter = ImageStatisticsFilterType::New(); imageStatisticsFilter->SetDirectionTolerance(0.001); imageStatisticsFilter->SetCoordinateTolerance(0.001); imageStatisticsFilter->SetInput(adaptedImage); imageStatisticsFilter->SetLabelInput(maskImage); imageStatisticsFilter->SetHistogramParametersForLabels(nBins, minVals, maxVals); imageStatisticsFilter->Update(); std::list labels = imageStatisticsFilter->GetRelevantLabels(); auto it = labels.begin(); while (it != labels.end()) { ImageStatisticsContainer::Pointer statisticContainerForLabelImage; auto labelIt = m_StatisticContainers.find(*it); // reset if statisticContainer already exist if (labelIt != m_StatisticContainers.end()) { statisticContainerForLabelImage = labelIt->second; } // create new statisticContainer else { statisticContainerForLabelImage = ImageStatisticsContainer::New(); statisticContainerForLabelImage->SetTimeGeometry(const_cast(timeGeometry)); // link label (*it) to statisticContainer m_StatisticContainers.emplace(*it, statisticContainerForLabelImage); } ImageStatisticsContainer::ImageStatisticsObject statObj; // find min, max, minindex and maxindex // make sure to only look in the masked region, use a masker for this vnl_vector minIndex, maxIndex; mitk::Point3D worldCoordinateMin; mitk::Point3D worldCoordinateMax; mitk::Point3D indexCoordinateMin; mitk::Point3D indexCoordinateMax; m_InternalImageForStatistics->GetGeometry()->IndexToWorld(minMaxFilter->GetMinIndex(*it), worldCoordinateMin); m_InternalImageForStatistics->GetGeometry()->IndexToWorld(minMaxFilter->GetMaxIndex(*it), worldCoordinateMax); m_Image->GetGeometry()->WorldToIndex(worldCoordinateMin, indexCoordinateMin); m_Image->GetGeometry()->WorldToIndex(worldCoordinateMax, indexCoordinateMax); minIndex.set_size(3); maxIndex.set_size(3); // for (unsigned int i=0; i < tmpMaxIndex.GetIndexDimension(); i++) for (unsigned int i = 0; i < 3; i++) { minIndex[i] = indexCoordinateMin[i]; maxIndex[i] = indexCoordinateMax[i]; } statObj.AddStatistic(mitk::ImageStatisticsConstants::MINIMUMPOSITION(), minIndex); statObj.AddStatistic(mitk::ImageStatisticsConstants::MAXIMUMPOSITION(), maxIndex); assert(std::abs(minMaxFilter->GetMax(*it) - imageStatisticsFilter->GetMaximum(*it)) < mitk::eps); assert(std::abs(minMaxFilter->GetMin(*it) - imageStatisticsFilter->GetMinimum(*it)) < mitk::eps); auto voxelVolume = GetVoxelVolume(image); auto numberOfVoxels = static_cast(imageStatisticsFilter->GetCount(*it)); auto volume = static_cast(numberOfVoxels) * voxelVolume; auto rms = std::sqrt(std::pow(imageStatisticsFilter->GetMean(*it), 2.) + imageStatisticsFilter->GetVariance(*it)); // variance = sigma^2 auto variance = imageStatisticsFilter->GetSigma(*it) * imageStatisticsFilter->GetSigma(*it); statObj.AddStatistic(mitk::ImageStatisticsConstants::NUMBEROFVOXELS(), numberOfVoxels); statObj.AddStatistic(mitk::ImageStatisticsConstants::VOLUME(), volume); statObj.AddStatistic(mitk::ImageStatisticsConstants::MEAN(), imageStatisticsFilter->GetMean(*it)); statObj.AddStatistic(mitk::ImageStatisticsConstants::MINIMUM(), imageStatisticsFilter->GetMinimum(*it)); statObj.AddStatistic(mitk::ImageStatisticsConstants::MAXIMUM(), imageStatisticsFilter->GetMaximum(*it)); statObj.AddStatistic(mitk::ImageStatisticsConstants::STANDARDDEVIATION(), imageStatisticsFilter->GetSigma(*it)); statObj.AddStatistic(mitk::ImageStatisticsConstants::VARIANCE(), variance); statObj.AddStatistic(mitk::ImageStatisticsConstants::SKEWNESS(), imageStatisticsFilter->GetSkewness(*it)); statObj.AddStatistic(mitk::ImageStatisticsConstants::KURTOSIS(), imageStatisticsFilter->GetKurtosis(*it)); statObj.AddStatistic(mitk::ImageStatisticsConstants::RMS(), rms); statObj.AddStatistic(mitk::ImageStatisticsConstants::MPP(), imageStatisticsFilter->GetMPP(*it)); statObj.AddStatistic(mitk::ImageStatisticsConstants::ENTROPY(), imageStatisticsFilter->GetEntropy(*it)); statObj.AddStatistic(mitk::ImageStatisticsConstants::MEDIAN(), imageStatisticsFilter->GetMedian(*it)); statObj.AddStatistic(mitk::ImageStatisticsConstants::UNIFORMITY(), imageStatisticsFilter->GetUniformity(*it)); statObj.AddStatistic(mitk::ImageStatisticsConstants::UPP(), imageStatisticsFilter->GetUPP(*it)); statObj.m_Histogram = imageStatisticsFilter->GetHistogram(*it).GetPointer(); statisticContainerForLabelImage->SetStatisticsForTimeStep(timeStep, statObj); ++it; } // swap maskGenerators back if (swapMasks) { m_SecondaryMask = m_InternalMask; m_InternalMask = nullptr; } } bool ImageStatisticsCalculator::IsUpdateRequired(LabelIndex label) const { unsigned long thisClassTimeStamp = this->GetMTime(); unsigned long inputImageTimeStamp = m_Image->GetMTime(); auto it = m_StatisticContainers.find(label); if (it == m_StatisticContainers.end()) { return true; } unsigned long statisticsTimeStamp = it->second->GetMTime(); if (thisClassTimeStamp > statisticsTimeStamp) // inputs have changed { return true; } if (inputImageTimeStamp > statisticsTimeStamp) // image has changed { return true; } if (m_MaskGenerator.IsNotNull()) { unsigned long maskGeneratorTimeStamp = m_MaskGenerator->GetMTime(); if (maskGeneratorTimeStamp > statisticsTimeStamp) // there is a mask generator and it has changed { return true; } } if (m_SecondaryMaskGenerator.IsNotNull()) { unsigned long maskGeneratorTimeStamp = m_SecondaryMaskGenerator->GetMTime(); if (maskGeneratorTimeStamp > statisticsTimeStamp) // there is a secondary mask generator and it has changed { return true; } } return false; } } // namespace mitk diff --git a/Modules/ImageStatistics/mitkStatisticsImageFilter.h b/Modules/ImageStatistics/mitkStatisticsImageFilter.h index 22c3db0225..18302d41d5 100644 --- a/Modules/ImageStatistics/mitkStatisticsImageFilter.h +++ b/Modules/ImageStatistics/mitkStatisticsImageFilter.h @@ -1,181 +1,191 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ // This file is based on ITK's itkStatisticsImageFilter.h #ifndef mitkStatisticsImageFilter #define mitkStatisticsImageFilter #include #include #include +#include #include #include #include #include namespace mitk { template class StatisticsImageFilter : public itk::ImageSink { public: - /** Standard Self type alias */ using Self = StatisticsImageFilter; using Superclass = itk::ImageSink; using Pointer = itk::SmartPointer; using ConstPointer = itk::SmartPointer; - /** Method for creation bypassing the object factory. */ itkFactorylessNewMacro(Self); - /** Runtime information support. */ itkTypeMacro(StatisticsImageFilter, itk::ImageSink); - /** Image related type alias. */ using InputImagePointer = typename TInputImage::Pointer; using RegionType = typename TInputImage::RegionType; using SizeType = typename TInputImage::SizeType; using IndexType = typename TInputImage::IndexType; using PixelType = typename TInputImage::PixelType; - /** Image related type alias. */ static constexpr unsigned int ImageDimension = TInputImage::ImageDimension; - /** Type to use for computations. */ using RealType = typename itk::NumericTraits::RealType; + + using HistogramType = typename itk::Statistics::Histogram; + using HistogramPointer = itk::SmartPointer; - /** Smart Pointer type to a DataObject. */ using DataObjectPointer = typename itk::DataObject::Pointer; template using SimpleDataObjectDecorator = itk::SimpleDataObjectDecorator; - /** Type of DataObjects used for scalar outputs */ using RealObjectType = SimpleDataObjectDecorator; using PixelObjectType = SimpleDataObjectDecorator; /** Return the computed Minimum. */ itkGetDecoratedOutputMacro(Minimum, PixelType); /** Return the computed Maximum. */ itkGetDecoratedOutputMacro(Maximum, PixelType); /** Return the computed Mean. */ itkGetDecoratedOutputMacro(Mean, RealType); /** Return the computed Standard Deviation. */ itkGetDecoratedOutputMacro(Sigma, RealType); /** Return the computed Variance. */ itkGetDecoratedOutputMacro(Variance, RealType); /** Return the computed Sum. */ itkGetDecoratedOutputMacro(Sum, RealType); /** Return the computed Sum of Squares. */ itkGetDecoratedOutputMacro(SumOfSquares, RealType); /** Return the computed Sum of Cubes. */ itkGetDecoratedOutputMacro(SumOfCubes, RealType); /** Return the computed Sum of Quadruples. */ itkGetDecoratedOutputMacro(SumOfQuadruples, RealType); /** Return the computed Skewness. */ itkGetDecoratedOutputMacro(Skewness, RealType); /** Return the computed Kurtosis. */ itkGetDecoratedOutputMacro(Kurtosis, RealType); /** Return the computed Mean of Positive Pixels. */ itkGetDecoratedOutputMacro(MPP, RealType); + /** Return the computed Histogram. */ + itkGetDecoratedOutputMacro(Histogram, HistogramPointer); + /** Return the computed Entropy. */ itkGetDecoratedOutputMacro(Entropy, RealType); /** Return the computed Uniformity. */ itkGetDecoratedOutputMacro(Uniformity, RealType); /** Return the computed UPP. */ itkGetDecoratedOutputMacro(UPP, RealType); /** Return the computed Median. */ itkGetDecoratedOutputMacro(Median, RealType); + void SetHistogramParameters(unsigned int size, RealType lowerBound, RealType upperBound); + using DataObjectIdentifierType = itk::ProcessObject::DataObjectIdentifierType; using Superclass::MakeOutput; /** Make a DataObject of the correct type to be used as the specified output. */ DataObjectPointer MakeOutput(const DataObjectIdentifierType& name) override; #ifdef ITK_USE_CONCEPT_CHECKING // Begin concept checking itkConceptMacro(InputHasNumericTraitsCheck, (itk::Concept::HasNumericTraits)); // End concept checking #endif protected: StatisticsImageFilter(); - ~StatisticsImageFilter() override = default; + ~StatisticsImageFilter(); itkSetDecoratedOutputMacro(Minimum, PixelType); itkSetDecoratedOutputMacro(Maximum, PixelType); itkSetDecoratedOutputMacro(Mean, RealType); itkSetDecoratedOutputMacro(Sigma, RealType); itkSetDecoratedOutputMacro(Variance, RealType); itkSetDecoratedOutputMacro(Sum, RealType); itkSetDecoratedOutputMacro(SumOfSquares, RealType); itkSetDecoratedOutputMacro(SumOfCubes, RealType); itkSetDecoratedOutputMacro(SumOfQuadruples, RealType); itkSetDecoratedOutputMacro(Skewness, RealType); itkSetDecoratedOutputMacro(Kurtosis, RealType); itkSetDecoratedOutputMacro(MPP, RealType); + itkSetDecoratedOutputMacro(Histogram, HistogramPointer); itkSetDecoratedOutputMacro(Entropy, RealType); itkSetDecoratedOutputMacro(Uniformity, RealType); itkSetDecoratedOutputMacro(UPP, RealType); itkSetDecoratedOutputMacro(Median, RealType); /** Initialize some accumulators before the threads run. */ void BeforeStreamedGenerateData() override; void ThreadedStreamedGenerateData(const RegionType&) override; /** Set outputs to computed values from all regions. */ void AfterStreamedGenerateData() override; void PrintSelf(std::ostream& os, itk::Indent indent) const override; private: - itk::CompensatedSummation m_ThreadSum; + HistogramPointer CreateInitializedHistogram() const; + + bool m_CalculateHistogram; + unsigned int m_HistogramSize; + PixelType m_HistogramLowerBound; + PixelType m_HistogramUpperBound; + HistogramPointer m_Histogram; + + itk::CompensatedSummation m_Sum; itk::CompensatedSummation m_SumOfPositivePixels; itk::CompensatedSummation m_SumOfSquares; itk::CompensatedSummation m_SumOfCubes; itk::CompensatedSummation m_SumOfQuadruples; itk::SizeValueType m_Count; itk::SizeValueType m_CountOfPositivePixels; - PixelType m_ThreadMin; - PixelType m_ThreadMax; + PixelType m_Min; + PixelType m_Max; std::mutex m_Mutex; }; } #ifndef ITK_MANUAL_INSTANTIATION #include #endif #endif diff --git a/Modules/ImageStatistics/mitkStatisticsImageFilter.hxx b/Modules/ImageStatistics/mitkStatisticsImageFilter.hxx index 0b97f4a8c7..215f8983d7 100644 --- a/Modules/ImageStatistics/mitkStatisticsImageFilter.hxx +++ b/Modules/ImageStatistics/mitkStatisticsImageFilter.hxx @@ -1,214 +1,323 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #ifndef mitkStatisticsImageFilter_hxx #define mitkStatisticsImageFilter_hxx #include #include template mitk::StatisticsImageFilter::StatisticsImageFilter() - : m_ThreadSum(1), + : m_CalculateHistogram(false), + m_HistogramSize(0), + m_HistogramLowerBound(itk::NumericTraits::NonpositiveMin()), + m_HistogramUpperBound(itk::NumericTraits::max()), + m_Sum(1), m_SumOfPositivePixels(1), m_SumOfSquares(1), m_SumOfCubes(1), m_SumOfQuadruples(1), m_Count(1), m_CountOfPositivePixels(1), - m_ThreadMin(1), - m_ThreadMax(1) + m_Min(1), + m_Max(1) { this->SetNumberOfRequiredInputs(1); this->SetMinimum(itk::NumericTraits::max()); this->SetMaximum(itk::NumericTraits::NonpositiveMin()); this->SetMean(itk::NumericTraits::max()); this->SetSigma(itk::NumericTraits::max()); this->SetVariance(itk::NumericTraits::max()); this->SetSum(0); this->SetSumOfSquares(0); this->SetSumOfCubes(0); this->SetSumOfQuadruples(0); this->SetSkewness(0); this->SetKurtosis(0); this->SetMPP(0); this->SetEntropy(-1.0); this->SetUniformity(0); this->SetUPP(0); this->SetMedian(0); } +template +mitk::StatisticsImageFilter::~StatisticsImageFilter() +{ +} + template typename mitk::StatisticsImageFilter::DataObjectPointer mitk::StatisticsImageFilter::MakeOutput(const DataObjectIdentifierType& name) { if (name == "Minimum" || name == "Maximum") { return PixelObjectType::New(); } if (name == "Mean" || name == "Sigma" || name == "Variance" || name == "Sum" || name == "SumOfSquares" || name == "SumOfCubes" || name == "SumOfQuadruples" || name == "Skewness" || name == "Kurtosis" || name == "MPP" || name == "Entropy" || name == "Uniformity" || name == "UPP" || name == "Median") { return RealObjectType::New(); } + if (name == "Histogram") + { + return HistogramType::New(); + } + return Superclass::MakeOutput(name); } +template +void mitk::StatisticsImageFilter::SetHistogramParameters(unsigned int size, RealType lowerBound, RealType upperBound) +{ + bool modified = false; + + if (m_HistogramSize != size) + { + m_HistogramSize = size; + modified = true; + } + + if (m_HistogramLowerBound != lowerBound) + { + m_HistogramLowerBound = lowerBound; + modified = true; + } + + if (m_HistogramUpperBound != upperBound) + { + m_HistogramUpperBound = upperBound; + modified = true; + } + + m_CalculateHistogram = true; + + if (modified) + this->Modified(); +} + +template +auto mitk::StatisticsImageFilter::CreateInitializedHistogram() const -> HistogramPointer +{ + typename HistogramType::SizeType size; + size.SetSize(1); + size[0] = m_HistogramSize; + + typename HistogramType::MeasurementVectorType lowerBound; + lowerBound.SetSize(1); + lowerBound[0] = m_HistogramLowerBound; + + typename HistogramType::MeasurementVectorType upperBound; + upperBound.SetSize(1); + upperBound.[0] = m_HistogramUpperBound; + + auto histogram = HistogramType::New(); + histogram->Initialize(size, lowerBound, upperBound); + + return histogram; +} + template void mitk::StatisticsImageFilter::BeforeStreamedGenerateData() { Superclass::BeforeStreamedGenerateData(); - m_ThreadSum = 0; + m_Sum = 0; m_SumOfPositivePixels = 0; m_SumOfSquares = 0; m_SumOfCubes = 0; m_SumOfQuadruples = 0; m_Count = 0; m_CountOfPositivePixels = 0; - m_ThreadMin = itk::NumericTraits::max(); - m_ThreadMax = itk::NumericTraits::NonpositiveMin(); + m_Min = itk::NumericTraits::max(); + m_Max = itk::NumericTraits::NonpositiveMin(); + + if (m_CalculateHistogram) + m_Histogram = this->CreateInitializedHistogram(); } template void mitk::StatisticsImageFilter::ThreadedStreamedGenerateData(const RegionType& regionForThread) { itk::CompensatedSummation sum = 0; itk::CompensatedSummation sumOfPositivePixels = 0; itk::CompensatedSummation sumOfSquares = 0; itk::CompensatedSummation sumOfCubes = 0; itk::CompensatedSummation sumOfQuadruples = 0; itk::SizeValueType count = 0; itk::SizeValueType countOfPositivePixels = 0; auto min = itk::NumericTraits::max(); auto max = itk::NumericTraits::NonpositiveMin(); + RealType realValue = 0; RealType squareValue = 0; + HistogramPointer histogram; + typename HistogramType::MeasurementVectorType histogramMeasurement; + typename HistogramType::IndexType histogramIndex; + + if (m_CalculateHistogram) // Initialize histogram + { + histogram = this->CreateInitializedHistogram(); + histogramMeasurement.SetSize(1); + } + itk::ImageScanlineConstIterator it(this->GetInput(), regionForThread); while (!it.IsAtEnd()) { while (!it.IsAtEndOfLine()) { const auto& value = it.Get(); - const auto realValue = static_cast(value); + realValue = static_cast(value); + + if (m_CalculateHistogram) // Compute histogram + { + histogramMeasurement[0] = realValue; + histogram->GetIndex(histogramMeasurement, histogramIndex); + histogram->IncreaseFrequencyOfIndex(histogramIndex, 1); + } + min = std::min(min, value); max = std::max(max, value); squareValue = realValue * realValue; sum += realValue; sumOfSquares += squareValue; sumOfCubes += squareValue * realValue; sumOfQuadruples += squareValue * squareValue; ++count; if (0 < realValue) { sumOfPositivePixels += realValue; ++countOfPositivePixels; } ++it; } it.NextLine(); } std::lock_guard mutexHolder(m_Mutex); - m_ThreadSum += sum; + if (m_CalculateHistogram) // Merge histograms + { + typename HistogramType::ConstIterator histogramIt = histogram->Begin(); + typename HistogramType::ConstIterator histogramEnd = histogram->End(); + + while (histogramIt != histogramEnd) + { + m_Histogram->GetIndex(histogramIt.GetMeasurementVector(), histogramIndex); + m_Histogram->IncreaseFrequencyOfIndex(histogramIndex, histogramIt.GetFrequency()); + ++histogramIt; + } + } + + m_Sum += sum; m_SumOfPositivePixels += sumOfPositivePixels; m_SumOfSquares += sumOfSquares; m_SumOfCubes += sumOfCubes; m_SumOfQuadruples += sumOfQuadruples; m_Count += count; m_CountOfPositivePixels += countOfPositivePixels; - m_ThreadMin = std::min(min, m_ThreadMin); - m_ThreadMax = std::max(max, m_ThreadMax); + m_Min = std::min(min, m_Min); + m_Max = std::max(max, m_Max); } template void mitk::StatisticsImageFilter::AfterStreamedGenerateData() { Superclass::AfterStreamedGenerateData(); - const RealType sum = m_ThreadSum; - const RealType sumOfPositivePixels = m_SumOfPositivePixels; - const RealType sumOfSquares = m_SumOfSquares; - const RealType sumOfCubes = m_SumOfCubes; - const RealType sumOfQuadruples = m_SumOfQuadruples; + const RealType sum = m_Sum.GetSum(); + const RealType sumOfPositivePixels = m_SumOfPositivePixels.GetSum(); + const RealType sumOfSquares = m_SumOfSquares.GetSum(); + const RealType sumOfCubes = m_SumOfCubes.GetSum(); + const RealType sumOfQuadruples = m_SumOfQuadruples.GetSum(); const itk::SizeValueType count = m_Count; const itk::SizeValueType countOfPositivePixels = m_CountOfPositivePixels; - const PixelType minimum = m_ThreadMin; - const PixelType maximum = m_ThreadMax; + const PixelType minimum = m_Min; + const PixelType maximum = m_Max; const RealType mean = sum / static_cast(count); const RealType variance = (sumOfSquares - (sum * sum / static_cast(count))) / (static_cast(count) - 1); const RealType sigma = std::sqrt(variance); const RealType secondMoment = sumOfSquares / static_cast(count); const RealType thirdMoment = sumOfCubes / static_cast(count); const RealType fourthMoment = sumOfQuadruples / static_cast(count); const RealType skewness = (thirdMoment - 3 * secondMoment * mean + 2 * std::pow(mean, 3)) / std::pow(secondMoment - std::pow(mean, 2), 1.5); const RealType kurtosis = (fourthMoment - 4 * thirdMoment * mean + 6 * secondMoment * std::pow(mean, 2) - 3 * std::pow(mean, 4)) / std::pow(secondMoment - std::pow(mean, 2), 2); const RealType meanOfPositivePixels = sumOfPositivePixels / static_cast(countOfPositivePixels); this->SetMinimum(minimum); this->SetMaximum(maximum); this->SetMean(mean); this->SetSigma(sigma); this->SetVariance(variance); this->SetSum(sum); this->SetSumOfSquares(sumOfSquares); this->SetSumOfCubes(sumOfCubes); this->SetSumOfQuadruples(sumOfQuadruples); this->SetSkewness(skewness); this->SetKurtosis(kurtosis); this->SetMPP(meanOfPositivePixels); - this->SetEntropy(entropy); - this->SetUniformity(uniformity); - this->SetUPP(upp); - this->SetMedian(median); + + if (m_CalculateHistogram) + { + this->SetHistogram(m_Histogram); + + mitk::HistogramStatisticsCalculator histogramStatisticsCalculator; + histogramStatisticsCalculator.SetHistogram(m_Histogram); + histogramStatisticsCalculator.CalculateStatistics(); + + this->SetEntropy(histogramStatisticsCalculator.GetEntropy()); + this->SetUniformity(histogramStatisticsCalculator.GetUniformity()); + this->SetUPP(histogramStatisticsCalculator.GetUPP()); + this->SetMedian(histogramStatisticsCalculator.GetMedian()); + } } template void mitk::StatisticsImageFilter::PrintSelf(std::ostream& os, itk::Indent indent) const { Superclass::PrintSelf(os, indent); os << indent << "SumOfCubes: " << this->GetSumOfCubes() << std::endl; os << indent << "SumOfQuadruples: " << this->GetSumOfQuadruples() << std::endl; os << indent << "Skewness: " << this->GetSkewness() << std::endl; os << indent << "Kurtosis: " << this->GetKurtosis() << std::endl; os << indent << "MPP: " << this->GetMPP() << std::endl; os << indent << "Entropy: " << this->GetEntropy() << std::endl; os << indent << "Uniformity: " << this->GetUniformity() << std::endl; os << indent << "UPP: " << this->GetUPP() << std::endl; os << indent << "Median: " << this->GetMedian() << std::endl; } #endif