diff --git a/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp b/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp index f126ae1dc5..6a3db779e7 100644 --- a/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp +++ b/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp @@ -1,546 +1,542 @@ +/*=================================================================== -#include -#include - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include +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 -#include #include #include #include #include #include #include #include #include #include #include - - #include -#include "itkImageFileWriter.h" - namespace mitk { void ImageStatisticsCalculator::SetInputImage(mitk::Image::Pointer image) { if (image != m_Image) { m_Image = image; m_StatisticsUpdateTimePerTimeStep.resize(m_Image->GetTimeSteps()); std::fill(m_StatisticsUpdateTimePerTimeStep.begin(), m_StatisticsUpdateTimePerTimeStep.end(), 0); this->Modified(); } } void ImageStatisticsCalculator::SetMask(mitk::MaskGenerator::Pointer mask) { if (mask != m_MaskGenerator) { m_MaskGenerator = mask; this->Modified(); } } void ImageStatisticsCalculator::SetSecondaryMask(mitk::MaskGenerator::Pointer mask) { if (mask != m_SecondaryMaskGenerator) { m_SecondaryMaskGenerator = mask; this->Modified(); } } void ImageStatisticsCalculator::SetNBinsForHistogramStatistics(unsigned int nBins) { if (nBins != m_nBinsForHistogramStatistics) { m_nBinsForHistogramStatistics = nBins; this->Modified(); this->m_UseBinSizeOverNBins = false; } if (m_UseBinSizeOverNBins) { this->Modified(); this->m_UseBinSizeOverNBins = false; } } unsigned int ImageStatisticsCalculator::GetNBinsForHistogramStatistics() const { return m_nBinsForHistogramStatistics; } void ImageStatisticsCalculator::SetBinSizeForHistogramStatistics(double binSize) { if (binSize != m_binSizeForHistogramStatistics) { m_binSizeForHistogramStatistics = binSize; this->Modified(); this->m_UseBinSizeOverNBins = true; } if (!m_UseBinSizeOverNBins) { this->Modified(); this->m_UseBinSizeOverNBins = true; } } double ImageStatisticsCalculator::GetBinSizeForHistogramStatistics() const { return m_binSizeForHistogramStatistics; } mitk::StatisticsContainer::StatisticsObject ImageStatisticsCalculator::GetStatistics(unsigned int timeStep, unsigned int label) { if (!m_StatisticContainers.empty() && timeStep >= m_StatisticContainers.back()->GetTimeSteps()) { mitkThrow() << "invalid timeStep in ImageStatisticsCalculator_v2::GetStatistics"; } if (m_Image.IsNull()) { mitkThrow() << "no image"; } if (!m_Image->IsInitialized()) { mitkThrow() << "Image not initialized!"; } auto timeGeometry = m_Image->GetTimeGeometry(); if (IsUpdateRequired(timeStep)) { if (m_MaskGenerator.IsNotNull()) { m_MaskGenerator->SetTimeStep(timeStep); m_InternalMask = m_MaskGenerator->GetMask(); if (m_MaskGenerator->GetReferenceImage().IsNotNull()) { m_InternalImageForStatistics = m_MaskGenerator->GetReferenceImage(); } else { m_InternalImageForStatistics = m_Image; } } else { m_InternalImageForStatistics = m_Image; } if (m_SecondaryMaskGenerator.IsNotNull()) { m_SecondaryMaskGenerator->SetTimeStep(timeStep); m_SecondaryMask = m_SecondaryMaskGenerator->GetMask(); } ImageTimeSelector::Pointer imgTimeSel = ImageTimeSelector::New(); imgTimeSel->SetInput(m_InternalImageForStatistics); imgTimeSel->SetTimeNr(timeStep); imgTimeSel->UpdateLargestPossibleRegion(); m_ImageTimeSlice = imgTimeSel->GetOutput(); // Calculate statistics with/without mask if (m_MaskGenerator.IsNull() && m_SecondaryMaskGenerator.IsNull()) { // 1) calculate statistics unmasked: AccessByItk_2(m_ImageTimeSlice, InternalCalculateStatisticsUnmasked, timeGeometry, timeStep) } else { // 2) calculate statistics masked AccessByItk_2(m_ImageTimeSlice, InternalCalculateStatisticsMasked, timeGeometry, timeStep) } //this->Modified(); } m_StatisticsUpdateTimePerTimeStep[timeStep] = m_StatisticContainers[m_StatisticContainers.size() - 1]->GetMTime(); for (auto it = m_StatisticContainers.begin(); it != m_StatisticContainers.end(); ++it) { StatisticsContainer::Pointer statCont = *it; StatisticsContainer::StatisticsObject statObject = statCont->GetStatisticsForTimeStep(timeStep); if (statObject.m_Label == label) { return statObject; } } // these lines will ony be executed if the requested label could not be found! MITK_WARN << "Invalid label: " << label << " in time step: " << timeStep; return StatisticsContainer::StatisticsObject(); } template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator::InternalCalculateStatisticsUnmasked( typename itk::Image< TPixel, VImageDimension >* image, TimeGeometry* timeGeometry, TimeStepType timeStep) { typedef typename itk::Image< TPixel, VImageDimension > ImageType; typedef typename itk::ExtendedStatisticsImageFilter ImageStatisticsFilterType; typedef typename itk::MinMaxImageFilterWithIndex MinMaxFilterType; // reset statistics container if queried time step exists if(m_StatisticContainers.empty() || m_StatisticContainers.back()->TimeStepExists(timeStep)) { m_StatisticContainers.clear(); StatisticsContainer::Pointer statisticsResult = StatisticsContainer::New(); statisticsResult->SetTimeGeometry(timeGeometry); m_StatisticContainers.push_back(statisticsResult); } //TODO for(timeStep) auto statObj = StatisticsContainer::StatisticsObject(); 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); // no mask, therefore just one label = the whole image statObj.m_Label = 1; 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(); m_StatisticContainers.back()->SetStatisticsForTimeStep(timeStep, statObj); } template < typename TPixel, unsigned int VImageDimension > 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 < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator::InternalCalculateStatisticsMasked( typename itk::Image< TPixel, VImageDimension >* image, TimeGeometry* timeGeometry, unsigned int timeStep) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< MaskPixelType, VImageDimension > MaskType; typedef typename MaskType::PixelType LabelPixelType; typedef itk::ExtendedLabelStatisticsImageFilter< ImageType, MaskType > ImageStatisticsFilterType; typedef MaskUtilities< TPixel, VImageDimension > MaskUtilType; typedef typename itk::MinMaxLabelImageFilterWithIndex MinMaxLabelFilterType; typedef typename ImageType::PixelType InputImgPixelType; // workaround: if m_SecondaryMaskGenerator ist not null but m_MaskGenerator is! (this is the case if we request a 'ignore zuero valued pixels' // mask in the gui but do not define a primary mask) bool swapMasks = false; if (m_SecondaryMask.IsNotNull() && m_InternalMask.IsNull()) { m_InternalMask = m_SecondaryMask; m_SecondaryMask = nullptr; swapMasks = true; } // maskImage has to have the same dimension as image typename MaskType::Pointer maskImage = MaskType::New(); try { // try to access the pixel values directly (no copying or casting). Only works if mask pixels are of pixelType unsigned short maskImage = ImageToItkImage< MaskPixelType, VImageDimension >(m_InternalMask); } catch (const itk::ExceptionObject &) { // if the pixel type of the mask is not short, then we have to make a copy of m_InternalMask (and cast the values) CastToItkImage(m_InternalMask, maskImage); } // if we have a secondary mask (say a ignoreZeroPixelMask) we need to combine the masks (corresponds to AND) if (m_SecondaryMask.IsNotNull()) { // dirty workaround for a bug when pf mask + any other mask is used in conjunction. We need a proper fix for this (Fabian Isensee is responsible and probably working on it!) if (m_InternalMask->GetDimension() == 2 && (m_SecondaryMask->GetDimension() == 3 || m_SecondaryMask->GetDimension() == 4)) { mitk::Image::Pointer old_img = m_SecondaryMaskGenerator->GetReferenceImage(); m_SecondaryMaskGenerator->SetInputImage(m_MaskGenerator->GetReferenceImage()); m_SecondaryMask = m_SecondaryMaskGenerator->GetMask(); m_SecondaryMaskGenerator->SetInputImage(old_img); } typename MaskType::Pointer secondaryMaskImage = MaskType::New(); secondaryMaskImage = ImageToItkImage< MaskPixelType, VImageDimension >(m_SecondaryMask); // secondary mask should be a ignore zero value pixel mask derived from image. it has to be cropped to the mask region (which may be planar or simply smaller) typename MaskUtilities::Pointer secondaryMaskMaskUtil = MaskUtilities::New(); secondaryMaskMaskUtil->SetImage(secondaryMaskImage.GetPointer()); secondaryMaskMaskUtil->SetMask(maskImage.GetPointer()); typename MaskType::Pointer adaptedSecondaryMaskImage = secondaryMaskMaskUtil->ExtractMaskImageRegion(); typename itk::MaskImageFilter2::Pointer maskFilter = itk::MaskImageFilter2::New(); maskFilter->SetInput1(maskImage); maskFilter->SetInput2(adaptedSecondaryMaskImage); maskFilter->SetMaskingValue(1); // all pixels of maskImage where secondaryMaskImage==1 will be kept, all the others are set to 0 maskFilter->UpdateLargestPossibleRegion(); maskImage = maskFilter->GetOutput(); } typename MaskUtilType::Pointer maskUtil = MaskUtilType::New(); maskUtil->SetImage(image); maskUtil->SetMask(maskImage.GetPointer()); // if mask is smaller than image, extract the image region where the mask is typename ImageType::Pointer adaptedImage = ImageType::New(); adaptedImage = maskUtil->ExtractMaskImageRegion(); // this also checks mask sanity // find min, max, minindex and maxindex typename MinMaxLabelFilterType::Pointer minMaxFilter = MinMaxLabelFilterType::New(); minMaxFilter->SetInput(adaptedImage); minMaxFilter->SetLabelInput(maskImage); minMaxFilter->UpdateLargestPossibleRegion(); // set histogram parameters for each label individually (min/max may be different for each label) typedef typename std::map MapType; typedef typename std::pair PairType; std::vector relevantLabels = minMaxFilter->GetRelevantLabels(); MapType minVals; MapType maxVals; std::map nBins; for (LabelPixelType label:relevantLabels) { minVals.insert(PairType(label, minMaxFilter->GetMin(label))); maxVals.insert(PairType(label, minMaxFilter->GetMax(label))); unsigned int nBinsForHistogram; if (m_UseBinSizeOverNBins) { nBinsForHistogram = std::max(static_cast( std::ceil(minMaxFilter->GetMax(label) - minMaxFilter->GetMin(label))) / m_binSizeForHistogramStatistics, 10.); // do not allow less than 10 bins } else { nBinsForHistogram = m_nBinsForHistogramStatistics; } nBins.insert(typename std::pair(label, nBinsForHistogram)); } typename ImageStatisticsFilterType::Pointer imageStatisticsFilter = ImageStatisticsFilterType::New(); imageStatisticsFilter->SetDirectionTolerance(0.001); imageStatisticsFilter->SetCoordinateTolerance(0.001); imageStatisticsFilter->SetInput(adaptedImage); imageStatisticsFilter->SetLabelInput(maskImage); imageStatisticsFilter->SetHistogramParametersForLabels(nBins, minVals, maxVals); imageStatisticsFilter->Update(); std::list labels = imageStatisticsFilter->GetRelevantLabels(); auto it = labels.begin(); if (m_StatisticContainers.empty() || m_StatisticContainers.back()->TimeStepExists(timeStep)) { m_StatisticContainers.clear(); for (size_t i = 0; i < labels.size(); i++) { StatisticsContainer::Pointer statisticsResult = StatisticsContainer::New(); statisticsResult->SetTimeGeometry(timeGeometry); m_StatisticContainers.push_back(statisticsResult); } } //m_StatisticsByTimeStep[timeStep].resize(0); while(it != labels.end()) { StatisticsContainer::StatisticsObject statObj = StatisticsContainer::StatisticsObject(); // 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->GetSum(*it) / (double)imageStatisticsFilter->GetMean(*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.m_Label = *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(); m_StatisticContainers.at(*it)->SetStatisticsForTimeStep(timeStep, statObj); ++it; } // swap maskGenerators back if (swapMasks) { m_SecondaryMask = m_InternalMask; m_InternalMask = nullptr; } } bool ImageStatisticsCalculator::IsUpdateRequired(unsigned int timeStep) const { unsigned long thisClassTimeStamp = this->GetMTime(); unsigned long inputImageTimeStamp = m_Image->GetMTime(); unsigned long statisticsTimeStamp = m_StatisticsUpdateTimePerTimeStep[timeStep]; if (thisClassTimeStamp > statisticsTimeStamp) // inputs have changed { return true; } if (inputImageTimeStamp > statisticsTimeStamp) // image has changed { return true; } if (m_MaskGenerator.IsNotNull()) { unsigned long maskGeneratorTimeStamp = m_MaskGenerator->GetMTime(); if (maskGeneratorTimeStamp > statisticsTimeStamp) // there is a mask generator and it has changed { return true; } } if (m_SecondaryMaskGenerator.IsNotNull()) { unsigned long maskGeneratorTimeStamp = m_SecondaryMaskGenerator->GetMTime(); if (maskGeneratorTimeStamp > statisticsTimeStamp) // there is a secondary mask generator and it has changed { return true; } } return false; } } diff --git a/Modules/ImageStatistics/mitkImageStatisticsCalculator.h b/Modules/ImageStatistics/mitkImageStatisticsCalculator.h index 469c7202bf..fc12245b0e 100644 --- a/Modules/ImageStatistics/mitkImageStatisticsCalculator.h +++ b/Modules/ImageStatistics/mitkImageStatisticsCalculator.h @@ -1,114 +1,128 @@ +/*=================================================================== + +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 MITKIMAGESTATISTICSCALCULATOR #define MITKIMAGESTATISTICSCALCULATOR #include #include #include #include #include -#include #include -#include namespace mitk { class MITKIMAGESTATISTICS_EXPORT ImageStatisticsCalculator: public itk::Object { public: /** Standard Self typedef */ typedef ImageStatisticsCalculator Self; typedef itk::Object Superclass; typedef itk::SmartPointer< Self > Pointer; typedef itk::SmartPointer< const Self > ConstPointer; /** Method for creation through the object factory. */ itkNewMacro(Self) /** Runtime information support. */ itkTypeMacro(ImageStatisticsCalculator_v2, itk::Object) typedef double statisticsValueType; typedef std::map statisticsMapType; typedef itk::Statistics::Histogram HistogramType; typedef unsigned short MaskPixelType; /**Documentation @brief Set the image for which the statistics are to be computed.*/ void SetInputImage(mitk::Image::Pointer image); /**Documentation @brief Set the mask generator that creates the mask which is to be used to calculate statistics. If no more mask is desired simply set @param mask to nullptr*/ void SetMask(mitk::MaskGenerator::Pointer mask); /**Documentation @brief Set this if more than one mask should be applied (for instance if a IgnorePixelValueMask were to be used alongside with a segmentation). Both masks are combined using pixel wise AND operation. The secondary mask does not have to be the same size than the primary but they need to have some overlap*/ void SetSecondaryMask(mitk::MaskGenerator::Pointer mask); /**Documentation @brief Set number of bins to be used for histogram statistics. If Bin size is set after number of bins, bin size will be used instead!*/ void SetNBinsForHistogramStatistics(unsigned int nBins); /**Documentation @brief Retrieve the number of bins used for histogram statistics. Careful: The return value does not indicate whether NBins or BinSize is used. That solely depends on which parameter has been set last.*/ unsigned int GetNBinsForHistogramStatistics() const; /**Documentation @brief Set bin size to be used for histogram statistics. If nbins is set after bin size, nbins will be used instead!*/ void SetBinSizeForHistogramStatistics(double binSize); /**Documentation @brief Retrieve the bin size for histogram statistics. Careful: The return value does not indicate whether NBins or BinSize is used. That solely depends on which parameter has been set last.*/ double GetBinSizeForHistogramStatistics() const; /**Documentation @brief Returns the statistics for label @a label and timeStep @a timeStep. If these requested statistics are not computed yet the computation is done as well. For performance reasons, statistics for all labels in the image are computed at once. */ StatisticsContainer::StatisticsObject GetStatistics(unsigned int timeStep=0, unsigned int label=1); protected: ImageStatisticsCalculator(){ m_nBinsForHistogramStatistics = 100; m_binSizeForHistogramStatistics = 10; m_UseBinSizeOverNBins = false; }; private: //Calculates statistics for each timestep for image template < typename TPixel, unsigned int VImageDimension > void InternalCalculateStatisticsUnmasked( typename itk::Image< TPixel, VImageDimension >* image, TimeGeometry* timeGeometry, TimeStepType timeStep); template < typename TPixel, unsigned int VImageDimension > void InternalCalculateStatisticsMasked( typename itk::Image< TPixel, VImageDimension >* image, TimeGeometry* timeGeometry, unsigned int timeStep); template < typename TPixel, unsigned int VImageDimension > double GetVoxelVolume(typename itk::Image* image) const; bool IsUpdateRequired(unsigned int timeStep) const; mitk::Image::Pointer m_Image; mitk::Image::Pointer m_ImageTimeSlice; mitk::Image::Pointer m_InternalImageForStatistics; mitk::MaskGenerator::Pointer m_MaskGenerator; mitk::Image::Pointer m_InternalMask; mitk::MaskGenerator::Pointer m_SecondaryMaskGenerator; mitk::Image::Pointer m_SecondaryMask; unsigned int m_nBinsForHistogramStatistics; double m_binSizeForHistogramStatistics; bool m_UseBinSizeOverNBins; std::vector m_StatisticContainers; std::vector m_StatisticsUpdateTimePerTimeStep; }; } #endif // MITKIMAGESTATISTICSCALCULATOR