diff --git a/Modules/ImageStatistics/files.cmake b/Modules/ImageStatistics/files.cmake index be094a5f7a..80a9d5d335 100644 --- a/Modules/ImageStatistics/files.cmake +++ b/Modules/ImageStatistics/files.cmake @@ -1,24 +1,28 @@ set(CPP_FILES mitkImageStatisticsCalculator.cpp mitkImageStatisticsCalculator2.cpp mitkPointSetStatisticsCalculator.cpp mitkPointSetDifferenceStatisticsCalculator.cpp mitkIntensityProfile.cpp mitkHotspotCalculator.cpp mitkMaskGenerator.cpp mitkPlanarFigureMaskGenerator.cpp mitkMultiLabelMaskGenerator.cpp + mitkBinaryImageMaskGenerator.cpp + mitkHistogramStatisticsCalculator.cpp ) set(H_FILES mitkImageStatisticsCalculator.h mitkImageStatisticsCalculator2.h mitkPointSetDifferenceStatisticsCalculator.h mitkPointSetStatisticsCalculator.h mitkExtendedStatisticsImageFilter.h mitkExtendedLabelStatisticsImageFilter.h mitkHotspotCalculator.h mitkMaskGenerator.h mitkPlanarFigureMaskGenerator.h mitkMultiLabelMaskGenerator.h + mitkBinaryImageMaskGenerator.h + mitkHistogramStatisticsCalculator.h ) diff --git a/Modules/ImageStatistics/mitkBinaryImageMaskGenerator.cpp b/Modules/ImageStatistics/mitkBinaryImageMaskGenerator.cpp new file mode 100644 index 0000000000..6713e02d75 --- /dev/null +++ b/Modules/ImageStatistics/mitkBinaryImageMaskGenerator.cpp @@ -0,0 +1,43 @@ + + +#include +#include +#include + +namespace mitk { + +void BinaryImageMaskGenerator::SetImageMask(Image::Pointer maskImage) +{ + if (m_internalMaskImage != maskImage) + { + m_internalMaskImage = maskImage; + m_Modified = true; + } +} + +mitk::Image::Pointer BinaryImageMaskGenerator::GetMask() +{ + if (m_Modified) + { + // extract timeStep -> is that correct? + if (m_internalMaskImage->GetTimeSteps() > 1) + { + if (m_TimeStep >= m_internalMaskImage->GetTimeSteps()) + { + throw std::runtime_error("Error: Invalid time step in BinaryImageMaskGenerator!"); + } + ImageTimeSelector::Pointer imageTimeSelector = ImageTimeSelector::New(); + imageTimeSelector->SetInput(m_internalMaskImage); + imageTimeSelector->SetTimeNr(m_TimeStep); + imageTimeSelector->UpdateLargestPossibleRegion(); + + m_InternalMask = mitk::Image::New(); + m_InternalMask = imageTimeSelector->GetOutput(); + m_Modified = false; + } + } + + return m_InternalMask; +} + +} diff --git a/Modules/ImageStatistics/mitkBinaryImageMaskGenerator.h b/Modules/ImageStatistics/mitkBinaryImageMaskGenerator.h new file mode 100644 index 0000000000..c1fbceaaa7 --- /dev/null +++ b/Modules/ImageStatistics/mitkBinaryImageMaskGenerator.h @@ -0,0 +1,27 @@ +#ifndef mitkBinaryMaskGenerator +#define mitkBinaryMaskGenerator + +#include +#include +#include + +namespace mitk +{ +class MITKIMAGESTATISTICS_EXPORT BinaryImageMaskGenerator: MaskGenerator +{ +public: + mitk::Image::Pointer GetMask(); + + void SetImageMask(mitk::Image::Pointer maskImage); + +protected: + +private: + mitk::Image::Pointer m_internalMaskImage; + +}; + + +} + +#endif diff --git a/Modules/ImageStatistics/mitkHistogramStatisticsCalculator.cpp b/Modules/ImageStatistics/mitkHistogramStatisticsCalculator.cpp new file mode 100644 index 0000000000..f2d2744608 --- /dev/null +++ b/Modules/ImageStatistics/mitkHistogramStatisticsCalculator.cpp @@ -0,0 +1,108 @@ +#include +#include + +namespace mitk { + +HistogramStatisticsCalculator::HistogramStatisticsCalculator(): + m_StatisticsCalculated(false) +{ + +} + +void HistogramStatisticsCalculator::SetHistogram(HistogramType::Pointer histogram) +{ + if (m_Histogram != histogram) + { + m_Histogram = histogram; + m_StatisticsCalculated = false; + } +} + +HistogramStatisticsCalculator::MeasurementType HistogramStatisticsCalculator::GetEntropy() const +{ + if (!m_StatisticsCalculated) + { + MITK_WARN("Statistics have not yet been calculated, running calculation now..."); + this->CalculateStatistics(); + } + return m_Entropy; +} + +HistogramStatisticsCalculator::MeasurementType HistogramStatisticsCalculator::GetMedian() const +{ + if (!m_StatisticsCalculated) + { + MITK_WARN("Statistics have not yet been calculated, running calculation now..."); + this->CalculateStatistics(); + } + return m_Median; +} + +HistogramStatisticsCalculator::MeasurementType HistogramStatisticsCalculator::GetUniformity() const +{ + if (!m_StatisticsCalculated) + { + MITK_WARN("Statistics have not yet been calculated, running calculation now..."); + this->CalculateStatistics(); + } + return m_Uniformity; +} + +HistogramStatisticsCalculator::MeasurementType HistogramStatisticsCalculator::GetUPP() const +{ + if (!m_StatisticsCalculated) + { + MITK_WARN("Statistics have not yet been calculated, running calculation now..."); + this->CalculateStatistics(); + } + return m_UPP; +} + +void HistogramStatisticsCalculator::CalculateStatistics() +{ + if (m_Histogram.IsNull()) + { + throw std::runtime_error("Histogram not set in HistogramStatisticsCalculator::CalculateStatistics()"); + } + + unsigned int nBins = m_Histogram->GetSize(); + m_Uniformity = 0; + m_Entropy = 0; + m_UPP = 0; + m_Median = 0; + + MeasurementType cumulativeProbability = 0.0; + MeasurementType partialProbability; + bool medianFound(false); + + for (unsigned int i = 0; i < nBins; i++) + { + partialProbability = m_Histogram->GetFrequency(i) / double( m_Histogram->GetTotalFrequency() ); + cumulativeProbability += partialProbability; + + if (partialProbability != 0) + { + m_Entropy -= partialProbability * (std::log2( partialProbability )); + m_Uniformity += std::pow(partialProbability, 2); + + if (m_Histogram->GetMeasurement(i, 0) > 0) + { + m_UPP += std::pow(partialProbability, 2); + } + + } + + if (cumulativeProbability >= 0.5 && !medianFound) + { + MeasurementType binMin = m_Histogram->GetBinMin(0, i); + MeasurementType binMax = m_Histogram->GetBinMax(0, i); + m_Median = (binMax + binMin) / 2.0; + medianFound = true; + } + + } + m_StatisticsCalculated = true; +} + +} + diff --git a/Modules/ImageStatistics/mitkHistogramStatisticsCalculator.h b/Modules/ImageStatistics/mitkHistogramStatisticsCalculator.h new file mode 100644 index 0000000000..a9c6cfc76c --- /dev/null +++ b/Modules/ImageStatistics/mitkHistogramStatisticsCalculator.h @@ -0,0 +1,40 @@ +#ifndef MITKHISTOGRAMSTATISTICSCALCULATOR +#define MITKHISTOGRAMSTATISTICSCALCULATOR + +#include +#include +#include + + +namespace mitk +{ + class MITKIMAGESTATISTICS_EXPORT HistogramStatisticsCalculator + { + public: + typedef double MeasurementType; + typedef itk::Statistics::Histogram HistogramType; + + HistogramStatisticsCalculator(); + + void SetHistogram(HistogramType::Pointer histogram); + + MeasurementType GetUPP() const; + + MeasurementType GetUniformity() const; + + MeasurementType GetEntropy() const; + + MeasurementType GetMedian() const; + + void CalculateStatistics(); + + protected: + + private: + HistogramType::Pointer m_Histogram; + MeasurementType m_Uniformity, m_UPP, m_Entropy, m_Median; + bool m_StatisticsCalculated; + }; +} + +#endif diff --git a/Modules/ImageStatistics/mitkImageStatisticsCalculator2.cpp b/Modules/ImageStatistics/mitkImageStatisticsCalculator2.cpp index 10e5f1e4cb..c46d6f8476 100644 --- a/Modules/ImageStatistics/mitkImageStatisticsCalculator2.cpp +++ b/Modules/ImageStatistics/mitkImageStatisticsCalculator2.cpp @@ -1,5 +1,564 @@ -#ifndef MITKIMAGESTATISTICSCALCULATOR2 -#define MITKIMAGESTATISTICSCALCULATOR2 -#endif // MITKIMAGESTATISTICSCALCULATOR2 +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +namespace mitk +{ + + void ImageStatisticsCalculator_v2::SetInputImage(mitk::Image::Pointer image) + { + if (image != m_Image) + { + m_Image = image; + m_HistogramStatisticsUpdateRequiredByTimeStep.resize(m_Image->GetTimeSteps()); + m_StatisticsUpdateRequiredByTimeStep.resize(m_Image->GetTimeSteps()); + m_HistogramsByTimeStep.resize(m_Image->GetTimeSteps()); + StatisticsContainer* stats = new StatisticsContainer(); + this->SetAllStatisticsToUpdateRequired(); + this->SetAllHistogramStatisticsToUpdateRequired(); + } + } + + void ImageStatisticsCalculator_v2::SetMask(mitk::MaskGenerator mask) + { + if (mask != m_MaskGenerator) + { + m_MaskGenerator = mask; + this->SetAllStatisticsToUpdateRequired(); + this->SetAllHistogramStatisticsToUpdateRequired(); + } + } + + + void ImageStatisticsCalculator_v2::SetDoHistogramStatistics(bool doHistStatistics) + { + m_DoHistogramStatistics = doHistStatistics; + } + + bool ImageStatisticsCalculator_v2::GetDoHistogramStatistics() const + { + return m_DoHistogramStatistics; + } + + void ImageStatisticsCalculator_v2::SetNBinsForHistogramStatistics(unsigned int nBins) + { + if (nBins != m_nBinsForHistogramStatistics) + { + m_nBinsForHistogramStatistics = nBins; + this->SetAllHistogramStatisticsToUpdateRequired(); + } + } + + unsigned int ImageStatisticsCalculator_v2::GetNBinsForHistogramStatistics() const + { + return m_nBinsForHistogramStatistics; + } + + ImageStatisticsCalculator_v2::statisticsMapType ImageStatisticsCalculator_v2::GetStatistics(unsigned int timeStep, unsigned int label) + { + 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 != nullptr) + { + m_MaskGenerator.SetTimeStep(timeStep); + m_InternalMask = m_MaskGenerator.GetMask(); + } + + if (m_StatisticsUpdateRequiredByTimeStep[timeStep]) + { + // Calculate statistics with/without mask + + if (m_MaskGenerator == nullptr) + { + // 1) calculate statistics unmasked: + // plug image into itkstatisticsimagefilter (will be replaced by my awesome filter later on) + // retrieve statistics and save them + + AccessByItk_1(m_Image, InternalCalculateStatisticsUnmasked, timeStep) + + } + else + { + // 2) calculate statistics masked + // extract mask image region + // plug mask and image into itklabelstatisticsimagefilter + AccessByItk_1(m_Image, InternalCalculateStatisticsMasked, timeStep) + + } + m_StatisticsUpdateRequiredByTimeStep[timeStep]=false; + } + + + if(m_DoHistogramStatistics && m_HistogramStatisticsUpdateRequiredByTimeStep[timeStep]) + { + + if (m_MaskGenerator == nullptr) + { + // calculate histogram unmasked + AccessByItk_1(m_Image, InternalCalculateHistogramUnmasked, timeStep); + m_HistogramStatisticsByTimeStep[timeStep].resize(1); + HistogramStatisticsContainer::Pointer histogramStatistics = this->InternalCalculateHistogramStatistics(m_HistogramsByTimeStep[timeStep][0].histogram); + histogramStatistics->SetLabel(0); + m_HistogramStatisticsByTimeStep[timeStep][0] = histogramStatistics; + } + else + { + // the label histograms have already been calculated in InternalCalculateStatisticsMasked + for (HistogramContainer histContainer:m_HistogramsByTimeStep[timeStep]) + { + HistogramStatisticsContainer::Pointer histogramStatistics = this->InternalCalculateHistogramStatistics(histContainer.histogram); + histogramStatistics->SetLabel(histContainer.label); + } + } + + + m_HistogramStatisticsUpdateRequiredByTimeStep[timeStep] = false; + + } + + + // merge histogram and regular statistics, return result + + } + + + 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]->Reset(); + } + + void ImageStatisticsCalculator_v2::SetAllHistogramStatisticsToUpdateRequired() + { + for (unsigned int i = 0; i < m_HistogramStatisticsUpdateRequiredByTimeStep.size(); i++) + { + this->SetHistStatsTimeStepToUpdateRequired(i); + } + } + + void ImageStatisticsCalculator_v2::SetHistStatsTimeStepToUpdateRequired(unsigned int timeStep) + { + if (timeStep >= m_HistogramStatisticsUpdateRequiredByTimeStep.size()) + { + throw std::runtime_error("invalid timeStep in ImageStatisticsCalculator_v2::SetHistStatsTimeStepToUpdateRequired"); + } + m_HistogramStatisticsUpdateRequiredByTimeStep[timeStep] = true; + m_HistogramStatisticsByTimeStep[timeStep]->Reset(); + } + + void ImageStatisticsCalculator_v2::CalculateHistogramStatistics(itk::Statistics::Histogram histogram, HistogramStatisticsContainer *HistStatContainer) + { + typedef itk::Statistics::Histogram HistogramType; + + } + + template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator_v2::InternalCheckMaskSanity( + typename itk::Image< TPixel, VImageDimension >::Pointer image, + typename itk::Image< unsigned int, VImageDimension >::Pointer maskImage) + { + 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; + typedef typename ImageType::DirectionType DirectionType; + + + // check direction + 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 << ")" ); + } + } + } + } + + // check spacing + PointType imageSpacing = image->GetSpacing(); + PointType maskSpacing = mask->GetSpacing(); + for (unsigned int i = 0; i < VImageDimension; i++) + { + if ( fabs( maskSpacing[i] - imageSpacing[i] ) > mitk::eps ) + { + itkExceptionMacro(<< "Spacing of mask and image is not equal. Mask: " << maskSpacing << " image: " << imageSpacing ); + } + } + + // check alignment + // 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(); + + 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 ); + // misalignment must be a multiple (int) of spacing in that direction + if ( misalignment%imageSpacing[i] > mitk::eps ) + { + itkWarningMacro( << "Pixels/voxels of mask and image are not sufficiently aligned! (Misalignment: " << misalignment%imageSpacing[i] << ")" ); + } + } + + // mask must be completely inside image region + // Make sure that mask region is contained within image region + if ( maskImage.IsNotNull() && + !image->GetLargestPossibleRegion().IsInside( maskImage->GetLargestPossibleRegion() ) ) + { + itkWarningMacro( << "Mask region needs to be inside of image region! (Image region: " + << image->GetLargestPossibleRegion() << "; Mask region: " << maskImage->GetLargestPossibleRegion() << ")" ); + } + + } + + template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator_v2::InternalCalculateStatisticsUnmasked( + typename itk::Image< TPixel, VImageDimension >::Pointer image, unsigned int timeStep) + { + typedef typename itk::Image< TPixel, VImageDimension > ImageType; + + typename itk::ExtendedStatisticsImageFilter::Pointer statisticsFilter = itk::ExtendedStatisticsImageFilter::New(); + statisticsFilter->SetInput(image); + statisticsFilter->SetCoordinateTolerance(0.001); + statisticsFilter->SetDirectionTolerance(0.001); + + 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_StatisticsUpdateRequiredByTimeStep[timeStep].resize(1); + StatisticsContainer::Pointer statisticsResult = StatisticsContainer::New(); + statisticsResult->GetLabel() = 0; + statisticsResult->GetN() = image->GetLargestPossibleRegion().GetNumberOfPixels(); + statisticsResult->GetMean() = statisticsFilter->GetMean(); + statisticsResult->GetMin() = statisticsFilter->GetMinimum(); + statisticsResult->GetMax() = statisticsFilter->GetMaximum(); + statisticsResult->GetVariance() = statisticsFilter->GetVariance(); + statisticsResult->GetStd() = statisticsFilter->GetSigma(); + statisticsResult->GetSkewness() = statisticsFilter->GetSkewness(); + statisticsResult->GetKurtosis() = statisticsFilter->GetKurtosis(); + statisticsResult->GetRMS() = statisticsFilter->GetRMS(); + statisticsResult->GetMPP() = statisticsFilter->GetMPP(); + m_StatisticsUpdateRequiredByTimeStep[timeStep][0] = statisticsResult; + + } + + template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator_v2::InternalCalculateHistogramUnmasked( + typename itk::Image< TPixel, VImageDimension >::Pointer image, + unsigned int timeStep) + { + typedef typename itk::Image< TPixel, VImageDimension > ImageType; + typedef itk::Statistics::ImageToHistogramFilter HistogramFilterType; + typedef itk::Statistics::Histogram HistogramType; + + typename HistogramFilterType::Pointer histogramFilter = HistogramFilterType::New(); + histogramFilter->SetInput(image); + histogramFilter->SetHistogramSize(m_nBinsForHistogramStatistics); // is that correct?? + histogramFilter->SetHistogramBinMaximum(statisticsResult->GetMax()); + histogramFilter->SetHistogramBinMinimum(statisticsResult->GetMin()); + + histogramFilter->Update(); + typename HistogramType::Pointer histogram = histogramFilter->GetOutput(); + + HistogramContainer histogramContainer; + histogramContainer.label=0; + histogramContainer.histogram = histogram; + m_HistogramsByTimeStep[timeStep].resize(1); + m_HistogramsByTimeStep[timeStep][0] = histogramContainer; + } + + + template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator_v2::ExtractMaskImageRegion( + typename itk::Image< TPixel, VImageDimension >::Pointer image, + typename itk::Image< unsigned int, VImageDimension >::Pointer maskImage, + typename itk::Image< TPixel, VImageDimension >::Pointer adaptedImage + ) + { + typedef itk::Image< TPixel, VImageDimension > ImageType; + typedef itk::Image< unsigned short, VImageDimension > MaskType; + typedef itk::ExtractImageFilter< ImageType, ImageType > ExtractImageFilterType; + + 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; + } + } + + if ( maskSmallerImage ) + { + typename ExtractImageFilterType::Pointer extractImageFilter = ExtractImageFilterType::New(); + extractImageFilter->SetInput( image ); + extractImageFilter->SetExtractionRegion( maskImage->GetBufferedRegion() ); + extractImageFilter->SetCoordinateTolerance( 0.001 ); + extractImageFilter->SetDirectionTolerance( 0.001 ); + extractImageFilter->Update(); + adaptedImage = extractImageFilter->GetOutput(); + } + else + { + adaptedImage = image; + } + } + + + /* Taken from original ImageStatisticsCalculator + * This function is not ideal on so many levels + 1) it does not support more than one label + 2) it assumes that the label id is 1 (which does not have to be the case) + 3) single threaded + => I suggest writing an itk filter that finds min and max value. This filter could then be called from our new + statisticsimagefilter in the beforethreadedgenereatedata (provided that a min/max for the histogram is not yet defined) + EDIT: I changed this function so that it simply searches the image min and max and completely disregards the mask. + Since we call this function after we cropped the image to the mask region we can run it on the LargestPossibleRegion of the image */ + template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator_v2::GetMinAndMaxValue( + typename itk::Image< TPixel, VImageDimension >::ConstPointer 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) + { + if( maskIterator.Value()== 1.0) + { + counter++; + + actualPixelValue = imageIterator.Value(); + + if(actualPixelValue >= maxVal) + { + maxVal = actualPixelValue; + } + else if(actualPielValue <= minVal) + { + minVal = actualPixelValue; + } + + } + } + } + + + template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator_v2::InternalCalculateStatisticsMasked( + typename itk::Image< TPixel, VImageDimension >::Pointer image, + unsigned int timeStep) + { + typedef itk::Image< TPixel, VImageDimension > ImageType; + typedef itk::Image< unsigned short, VImageDimension > MaskType; + typedef itk::ExtractImageFilter< ImageType, ImageType > ExtractImageFilterType; + typedef itk::ExtendedLabelStatisticsImageFilter< ImageType, MaskType > ImageStatisticsFilterType; + + // maskImage has to have the same dimension as image + MaskType::Pointer maskImage = MaskType::New(); + maskImage = ImageToItkImage(m_InternalMask); + + InternalCheckMaskSanity(image, maskImage); + + // if mask is smaller than image, extract the image region where the mask is + typename ImageType::Pointer adaptedImage = ImageType::New(); + + // uses m_InternalMask + this->ExtractMaskImageRegion(image, maskImage, adaptedImage); + + double minVal, maxVal; + GetMinAndMaxValue(adaptedImage, minVal, maxVal); + + ImageStatisticsFilterType::Pointer imageStatisticsFilter = ImageStatisticsFilterType::New(); + imageStatisticsFilter->SetDirectionTolerance(0.001); + imageStatisticsFilter->SetCoordinateTolerance(0.001); + imageStatisticsFilter->SetHistogramParameters(m_nBinsForHistogramStatistics, minVal, maxVal); + imageStatisticsFilter->SetInput(adaptedImage); + imageStatisticsFilter->SetLabelInput(maskImage); + 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(); + statisticsResult->GetN() = imageStatisticsFilter->GetSum(*it) / (double) imageStatisticsFilter->GetMean(*it); + statisticsResult->GetMean() = imageStatisticsFilter->GetMean(*it); + statisticsResult->GetMin() = imageStatisticsFilter->GetMinimum(*it); + statisticsResult->GetMax() = imageStatisticsFilter->GetMaximum(*it); + statisticsResult->GetVariance() = imageStatisticsFilter->GetVariance(*it); + statisticsResult->GetStd() = imageStatisticsFilter->GetSigma(1); + statisticsResult->GetSkewness() = imageStatisticsFilter->GetSkewness(*it); + statisticsResult->GetKurtosis() = imageStatisticsFilter->GetKurtosis(*it); + statisticsResult->GetRMS() = imageStatisticsFilter->GetRMS(*it); + statisticsResult->GetMPP() = imageStatisticsFilter->GetMPP(*it); + statisticsResult->GetLabel() = *it; + m_StatisticsByTimeStep[timeStep].push_back(statisticsResult); + if (m_DoHistogramStatistics) + { + m_HistogramsByTimeStep[timeStep].push_back(imageStatisticsFilter->GetHistogram(*it)); + } + } + } + + ImageStatisticsCalculator_v2::HistogramStatisticsContainer::Pointer ImageStatisticsCalculator_v2::InternalCalculateHistogramStatistics(itk::Statistics::Histogram::Pointer histogram) + { + HistogramStatisticsCalculator histStatCalc; + histStatCalc.SetHistogram(histogram); + histStatCalc.CalculateStatistics(); + HistogramStatisticsContainer::Pointer histStatResults = HistogramStatisticsContainer::New(); + histStatResults->GetEntropy() = histStatCalc.GetEntropy(); + histStatResults->GetMedian() = histStatCalc.GetMedian(); + histStatResults->GetUniformity() = histStatCalc.GetUniformity(); + histStatResults->GetUPP() = histStatCalc.GetUPP(); + return histStatResults; + } + + ImageStatisticsCalculator_v2::StatisticsContainer::StatisticsContainer(): + 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()), + { + + } + + 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; + + 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(); + } + + + ImageStatisticsCalculator_v2::HistogramStatisticsContainer::HistogramStatisticsContainer(): + m_Median(std::numeric_limits::max()), + m_Uniformity(std::numeric_limits::max()), + m_UPP(std::numeric_limits::max()), + m_Entropy(std::numeric_limits::max()), + { + + } + + ImageStatisticsCalculator_v2::statisticsMapType ImageStatisticsCalculator_v2::HistogramStatisticsContainer::GetStatisticsAsMap() + { + ImageStatisticsCalculator_v2::statisticsMapType statisticsAsMap; + + statisticsAsMap['Median'] = m_Median; + statisticsAsMap['Uniformity'] = m_Uniformity; + statisticsAsMap['UPP'] = m_UPP; + statisticsAsMap['Entropy'] = m_Entropy; + + return statisticsAsMap; + } + + + void ImageStatisticsCalculator_v2::HistogramStatisticsContainer::Reset() + { + m_Median = std::numeric_limits::max(); + m_Uniformity = std::numeric_limits::max(); + m_UPP = std::numeric_limits::max(); + m_Entropy = std::numeric_limits::max(); + + } +} diff --git a/Modules/ImageStatistics/mitkImageStatisticsCalculator2.h b/Modules/ImageStatistics/mitkImageStatisticsCalculator2.h index 10e5f1e4cb..bf81958697 100644 --- a/Modules/ImageStatistics/mitkImageStatisticsCalculator2.h +++ b/Modules/ImageStatistics/mitkImageStatisticsCalculator2.h @@ -1,5 +1,178 @@ #ifndef MITKIMAGESTATISTICSCALCULATOR2 #define MITKIMAGESTATISTICSCALCULATOR2 +#include +#include +#include +#include +#include +#include + + +#define mitkGetRefConstMacro(name, type) \ + virtual const type& Get##name() const \ + { \ + return &m_##name; \ + } \ + \ + virtual type& Get##name() \ + { \ + return &m_##name; \ + } \ + + +namespace mitk +{ + class MITKIMAGESTATISTICS_EXPORT ImageStatisticsCalculator_v2 + { + public: + typedef double statisticsValueType; + typedef std::map statisticsMapType; + typedef itk::Statistics::Histogram HistogramType; + + void SetInputImage(mitk::Image::Pointer image); + + void SetMask(mitk::MaskGenerator mask); + + void SetDoHistogramStatistics(bool doHistStatistics); + + bool GetDoHistogramStatistics() const; + + void SetNBinsForHistogramStatistics(unsigned int nBins); + + unsigned int GetNBinsForHistogramStatistics() const; + + statisticsMapType GetStatistics(unsigned int timeStep=0, unsigned int label=0); // not yet implemented + + class StatisticsContainer:itk::Object + { + public: + typedef double RealType; + + statisticsMapType GetStatisticsAsMap() const; + + void Reset(); + + mitkGetRefConstMacro(N, long) + mitkGetRefConstMacro(Mean, RealType) + mitkGetRefConstMacro(Variance, RealType) + mitkGetRefConstMacro(Std, RealType) + mitkGetRefConstMacro(Min, RealType) + mitkGetRefConstMacro(Max, RealType) + mitkGetRefConstMacro(RMS, RealType) + mitkGetRefConstMacro(Skewness, RealType) + mitkGetRefConstMacro(Kurtosis, RealType) + mitkGetRefConstMacro(MPP, RealType) + mitkGetRefConstMacro(Label, unsigned int) + + protected: + StatisticsContainer(); + + + private: + // not pretty, is temporary + RealType m_Mean, m_Variance, m_Std, m_Min, m_Max, m_RMS; + RealType m_Skewness; + RealType m_Kurtosis; + RealType m_MPP; + long m_N; + unsigned int m_Label; + + }; + + class HistogramStatisticsContainer:itk::Object + { + public: + typedef double RealType; + + + statisticsMapType GetStatisticsAsMap() const; + + void Reset(); + + mitkGetRefConstMacro(Median, RealType) + mitkGetRefConstMacro(Uniformity, RealType) + mitkGetRefConstMacro(Entropy, RealType) + mitkGetRefConstMacro(UPP, RealType) + mitkGetRefConstMacro(Label, unsigned int) + + protected: + HistogramStatisticsContainer(); + + private: + RealType m_Median; + RealType m_Uniformity; + RealType m_Entropy; + RealType m_UPP; + unsigned int m_Label; + + }; + + struct HistogramContainer + { + unsigned int label; + itk::Statistics::Histogram::Pointer histogram; + }; + + protected: + + private: + void SetAllStatisticsToUpdateRequired(); + void SetAllHistogramStatisticsToUpdateRequired(); + + void SetStatsTimeStepToUpdateRequired(unsigned int timeStep); + void SetHistStatsTimeStepToUpdateRequired(unsigned int timeStep); + + void CalculateHistogramStatistics(itk::Statistics::Histogram histogram, HistogramStatisticsContainer* HistStatContainer); + + template < typename TPixel, unsigned int VImageDimension > void InternalCheckMaskSanity( + typename itk::Image< TPixel, VImageDimension >::Pointer image, + typename itk::Image< unsigned int, VImageDimension >::Pointer maskImage); + + template < typename TPixel, unsigned int VImageDimension > void InternalCalculateStatisticsUnmasked( + typename itk::Image< TPixel, VImageDimension >::Pointer image, + unsigned int timeStep); + + HistogramStatisticsContainer::Pointer InternalCalculateHistogramStatistics(itk::Statistics::Histogram::Pointer histogram); + + template < typename TPixel, unsigned int VImageDimension > void InternalCalculateHistogramUnmasked( + typename itk::Image< TPixel, VImageDimension >::Pointer image, + unsigned int timeStep); + + template < typename TPixel, unsigned int VImageDimension > void InternalCalculateStatisticsMasked( + typename itk::Image< TPixel, VImageDimension >::Pointer image, + unsigned int timeStep); + + template < typename TPixel, unsigned int VImageDimension > void ExtractMaskImageRegion( + typename itk::Image< TPixel, VImageDimension >::Pointer image, + typename itk::Image< unsigned int, VImageDimension >::Pointer maskImage, + typename itk::Image< TPixel, VImageDimension >::Pointer adaptedImage + ); + + template < typename TPixel, unsigned int VImageDimension > void GetMinAndMaxValue( + typename itk::Image< TPixel, VImageDimension >::ConstPointer inputImage, + double &minVal, + double &maxVal + ); + + mitk::Image::Pointer m_Image; + mitk::Image::Pointer m_ImageTimeSlice; + + mitk::MaskGenerator m_MaskGenerator; + mitk::Image::Pointer m_InternalMask; + + bool m_DoHistogramStatistics; + unsigned int m_nBinsForHistogramStatistics; + + std::vector m_StatisticsUpdateRequiredByTimeStep; // holds which time steps are valid and which ones have to be recalculated + std::vector> m_StatisticsByTimeStep; + + std::vector> m_HistogramStatisticsByTimeStep; + std::vector m_HistogramStatisticsUpdateRequiredByTimeStep; + + std::vector> m_HistogramsByTimeStep; + }; + +} #endif // MITKIMAGESTATISTICSCALCULATOR2 diff --git a/Modules/ImageStatistics/mitkMaskGenerator.h b/Modules/ImageStatistics/mitkMaskGenerator.h index 7c19bf6b6b..c747417633 100644 --- a/Modules/ImageStatistics/mitkMaskGenerator.h +++ b/Modules/ImageStatistics/mitkMaskGenerator.h @@ -1,33 +1,34 @@ #ifndef MITKMASKGENERATOR #define MITKMASKGENERATOR #include #include #include namespace mitk { class MITKIMAGESTATISTICS_EXPORT MaskGenerator { public: MaskGenerator(); //~MaskGenerator(); virtual mitk::Image::Pointer GetMask() = 0; void SetTimeStep(unsigned int timeStep); protected: + void CheckMaskSanity(mitk::Image::Pointer image); + unsigned int m_TimeStep; mitk::Image::Pointer m_InternalMask; bool m_Modified; private: - void CheckMaskSanity(mitk::Image::Pointer); }; } #endif // MITKMASKGENERATOR diff --git a/Modules/ImageStatistics/mitkMultiLabelMaskGenerator.h b/Modules/ImageStatistics/mitkMultiLabelMaskGenerator.h index 4a0a0be26f..cc84356104 100644 --- a/Modules/ImageStatistics/mitkMultiLabelMaskGenerator.h +++ b/Modules/ImageStatistics/mitkMultiLabelMaskGenerator.h @@ -1,25 +1,41 @@ #ifndef MITKMULTILABELMASKGENERATOR #define MITKMULTILABELMASKGENERATOR #include #include #include #include + namespace mitk { class MITKIMAGESTATISTICS_EXPORT MultiLabelMaskGenerator: public MaskGenerator { public: - void SetLabelSetImage(mitk::LabelSetImage::Pointer labelSetImage); + void setLabelSetImage(mitk::LabelSetImage::Pointer labelSetImage); + + void addLabel(LabelSetImage::PixelType, std::vector::size_type layer=0); + void removeLabel(LabelSetImage::PixelType, std::vector::size_type layer=0); + + void addLabels(std::pair> labelsToAdd); + void removeLabels(std::pair> labelsToAdd); + + void addLabels(std::vector labels, std::vector::size_type layer=0); + void removeLabels(std::vector labels, std::vector::size_type layer=0); + + void removeLayer(std::vector::size_type layer); + + mitk::Image::Pointer GetMask(); protected: private: + mitk::LabelSetImage::Pointer m_LabelSetImage; + std::vector> m_selectedLabels; }; } #endif diff --git a/Modules/MaskGenerator/CMakeLists.txt b/Modules/MaskGenerator/CMakeLists.txt new file mode 100644 index 0000000000..d073a4481d --- /dev/null +++ b/Modules/MaskGenerator/CMakeLists.txt @@ -0,0 +1,13 @@ +MITK_CREATE_MODULE( + DEPENDS MitkImageStatistics + #PACKAGE_DEPENDS + # PUBLIC ITK|ITKIOXML + # PRIVATE ITK|ITKVTK+ITKConvolution + # WARNINGS_AS_ERRORS +) + +if(BUILD_TESTING) + + # add_subdirectory(Testing) + +endif(BUILD_TESTING) diff --git a/Modules/MaskGenerator/files.cmake b/Modules/MaskGenerator/files.cmake new file mode 100644 index 0000000000..36e7a3e616 --- /dev/null +++ b/Modules/MaskGenerator/files.cmake @@ -0,0 +1,7 @@ +set(CPP_FILES + +) + +set(H_FILES + +) diff --git a/Modules/MaskGenerator/mitkDeleteMe.cpp b/Modules/MaskGenerator/mitkDeleteMe.cpp new file mode 100644 index 0000000000..8f5253af72 --- /dev/null +++ b/Modules/MaskGenerator/mitkDeleteMe.cpp @@ -0,0 +1,6 @@ +#include + +void deleteMe::poebelMalHerum() +{ + std::cout << "was fuer ein dummes Modul" << std::endl; +} \ No newline at end of file diff --git a/Modules/MaskGenerator/mitkDeleteMe.h b/Modules/MaskGenerator/mitkDeleteMe.h new file mode 100644 index 0000000000..e2b9473478 --- /dev/null +++ b/Modules/MaskGenerator/mitkDeleteMe.h @@ -0,0 +1,7 @@ +#include + +class deleteMe +{ +public: +void poebelMalHerum(); +}; diff --git a/Modules/ModuleList.cmake b/Modules/ModuleList.cmake index 63e0a822c9..8b28890951 100644 --- a/Modules/ModuleList.cmake +++ b/Modules/ModuleList.cmake @@ -1,78 +1,78 @@ # The entries in the mitk_modules list must be # ordered according to their dependencies. set(mitk_modules Core CommandLine AppUtil DCMTesting RDF LegacyIO DataTypesExt Overlays LegacyGL AlgorithmsExt MapperExt DICOMReader DICOMReaderServices DICOMTesting Qt4Qt5TestModule SceneSerializationBase PlanarFigure ImageDenoising ImageExtraction - ImageStatistics LegacyAdaptors SceneSerialization Gizmo GraphAlgorithms Multilabel + ImageStatistics ContourModel SurfaceInterpolation Segmentation PlanarFigureSegmentation OpenViewCore QmlItems QtWidgets QtWidgetsExt SegmentationUI DiffusionImaging GPGPU OpenIGTLink IGTBase IGT CameraCalibration RigidRegistration RigidRegistrationUI DeformableRegistration DeformableRegistrationUI OpenCL OpenCVVideoSupport QtOverlays ToFHardware ToFProcessing ToFUI US USUI DicomUI Simulation Remeshing Python QtPython Persistence OpenIGTLinkUI IGTUI VtkShaders DicomRT RTUI IOExt XNAT TubeGraph BiophotonicsHardware Classification TumorInvasionAnalysis ) if(MITK_ENABLE_PIC_READER) list(APPEND mitk_modules IpPicSupportIO) endif()