diff --git a/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp b/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp index bc01620bdb..32eaec9440 100644 --- a/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp +++ b/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp @@ -1,517 +1,536 @@ #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include "itkImageFileWriter.h" namespace mitk { void ImageStatisticsCalculator::SetInputImage(mitk::Image::Pointer image) { if (image != m_Image) { m_Image = image; - m_StatisticsByTimeStep.resize(m_Image->GetTimeSteps()); m_StatisticsUpdateTimePerTimeStep.resize(m_Image->GetTimeSteps()); std::fill(m_StatisticsUpdateTimePerTimeStep.begin(), m_StatisticsUpdateTimePerTimeStep.end(), 0); this->Modified(); } } void ImageStatisticsCalculator::SetMask(mitk::MaskGenerator::Pointer mask) { if (mask != m_MaskGenerator) { m_MaskGenerator = mask; this->Modified(); } } void ImageStatisticsCalculator::SetSecondaryMask(mitk::MaskGenerator::Pointer mask) { if (mask != m_SecondaryMaskGenerator) { m_SecondaryMaskGenerator = mask; this->Modified(); } } void ImageStatisticsCalculator::SetNBinsForHistogramStatistics(unsigned int nBins) { if (nBins != m_nBinsForHistogramStatistics) { m_nBinsForHistogramStatistics = nBins; this->Modified(); this->m_UseBinSizeOverNBins = false; } if (m_UseBinSizeOverNBins) { this->Modified(); this->m_UseBinSizeOverNBins = false; } } unsigned int ImageStatisticsCalculator::GetNBinsForHistogramStatistics() const { return m_nBinsForHistogramStatistics; } void ImageStatisticsCalculator::SetBinSizeForHistogramStatistics(double binSize) { if (binSize != m_binSizeForHistogramStatistics) { m_binSizeForHistogramStatistics = binSize; this->Modified(); this->m_UseBinSizeOverNBins = true; } if (!m_UseBinSizeOverNBins) { this->Modified(); this->m_UseBinSizeOverNBins = true; } } double ImageStatisticsCalculator::GetBinSizeForHistogramStatistics() const { return m_binSizeForHistogramStatistics; } - mitk::StatisticsContainer::Pointer ImageStatisticsCalculator::GetStatistics(unsigned int timeStep, unsigned int label) + mitk::StatisticsContainer::StatisticsObject ImageStatisticsCalculator::GetStatistics(unsigned int timeStep, unsigned int label) { - - if (timeStep >= m_StatisticsByTimeStep.size()) + 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_1(m_ImageTimeSlice, InternalCalculateStatisticsUnmasked, timeStep) + AccessByItk_2(m_ImageTimeSlice, InternalCalculateStatisticsUnmasked, timeGeometry, timeStep) } else { // 2) calculate statistics masked - AccessByItk_1(m_ImageTimeSlice, InternalCalculateStatisticsMasked, timeStep) + AccessByItk_2(m_ImageTimeSlice, InternalCalculateStatisticsMasked, timeGeometry, timeStep) } - //this->Modified(); } - m_StatisticsUpdateTimePerTimeStep[timeStep] = m_StatisticsByTimeStep[timeStep][m_StatisticsByTimeStep[timeStep].size()-1]->GetMTime(); + m_StatisticsUpdateTimePerTimeStep[timeStep] = m_StatisticContainers[m_StatisticContainers.size() - 1]->GetMTime(); - for (auto it = m_StatisticsByTimeStep[timeStep].begin(); it != m_StatisticsByTimeStep[timeStep].end(); ++it) + for (auto it = m_StatisticContainers.begin(); it != m_StatisticContainers.end(); ++it) { StatisticsContainer::Pointer statCont = *it; - if (statCont->GetLabel() == label) + StatisticsContainer::StatisticsObject statObject = statCont->GetStatisticsForTimeStep(timeStep); + if (statObject.m_Label == label) { - return statCont; + 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::New(); + return StatisticsContainer::StatisticsObject(); } template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator::InternalCalculateStatisticsUnmasked( - typename itk::Image< TPixel, VImageDimension >* image, unsigned int timeStep) + 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); + } - StatisticsContainer::Pointer statisticsResult = StatisticsContainer::New(); + //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]; } - statisticsResult->SetMinIndex(minIndex); - statisticsResult->SetMaxIndex(maxIndex); + statObj.m_MinIndex = minIndex; + statObj.m_MaxIndex = 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 - m_StatisticsByTimeStep[timeStep].resize(1); - statisticsResult->SetLabel(1); - statisticsResult->SetN(image->GetLargestPossibleRegion().GetNumberOfPixels()); - statisticsResult->SetVolume(statisticsResult->GetN()*voxelVolume); - statisticsResult->SetMean(statisticsFilter->GetMean()); - statisticsResult->SetMin(statisticsFilter->GetMinimum()); - statisticsResult->SetMax(statisticsFilter->GetMaximum()); - statisticsResult->SetStd(statisticsFilter->GetSigma()); - statisticsResult->SetSkewness(statisticsFilter->GetSkewness()); - statisticsResult->SetKurtosis(statisticsFilter->GetKurtosis()); - statisticsResult->SetRMS(std::sqrt(std::pow(statisticsFilter->GetMean(), 2.) + statisticsFilter->GetVariance())); // variance = sigma^2 - statisticsResult->SetMPP(statisticsFilter->GetMPP()); - - statisticsResult->SetEntropy(statisticsFilter->GetEntropy()); - statisticsResult->SetMedian(statisticsFilter->GetMedian()); - statisticsResult->SetUniformity(statisticsFilter->GetUniformity()); - statisticsResult->SetUPP(statisticsFilter->GetUPP()); - statisticsResult->SetHistogram(statisticsFilter->GetHistogram()); - - m_StatisticsByTimeStep[timeStep][0] = statisticsResult; + statObj.m_Label = 1; + statObj.m_N = image->GetLargestPossibleRegion().GetNumberOfPixels(); + statObj.m_Volume = statObj.m_N * voxelVolume; + statObj.m_Mean = statisticsFilter->GetMean(); + statObj.m_Min = statisticsFilter->GetMinimum(); + statObj.m_Max = statisticsFilter->GetMaximum(); + statObj.m_Std = statisticsFilter->GetSigma(); + statObj.m_Skewness = statisticsFilter->GetSkewness(); + statObj.m_Kurtosis = statisticsFilter->GetKurtosis(); + statObj.m_RMS = std::sqrt(std::pow(statisticsFilter->GetMean(), 2.) + statisticsFilter->GetVariance()); // variance = sigma^2 + statObj.m_MPP = statisticsFilter->GetMPP(); + + statObj.m_Entropy = statisticsFilter->GetEntropy(); + statObj.m_Median = statisticsFilter->GetMedian(); + statObj.m_Uniformity = statisticsFilter->GetUniformity(); + statObj.m_UPP = statisticsFilter->GetUPP(); + statObj.m_Histogram = statisticsFilter->GetHistogram(); + + 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, + 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(); - m_StatisticsByTimeStep[timeStep].resize(0); + + if (m_StatisticContainers.empty() || m_StatisticContainers.back()->TimeStepExists(timeStep)) { + m_StatisticContainers.clear(); + for (int 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::Pointer statisticsResult = StatisticsContainer::New(); + 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]; } - statisticsResult->SetMinIndex(minIndex); - statisticsResult->SetMaxIndex(maxIndex); + statObj.m_MinIndex = minIndex; + statObj.m_MaxIndex = 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); - statisticsResult->SetN(imageStatisticsFilter->GetSum(*it) / (double) imageStatisticsFilter->GetMean(*it)); - statisticsResult->SetVolume(statisticsResult->GetN()*voxelVolume); - statisticsResult->SetMean(imageStatisticsFilter->GetMean(*it)); - statisticsResult->SetMin(imageStatisticsFilter->GetMinimum(*it)); - statisticsResult->SetMax(imageStatisticsFilter->GetMaximum(*it)); - statisticsResult->SetStd(imageStatisticsFilter->GetSigma(*it)); - statisticsResult->SetSkewness(imageStatisticsFilter->GetSkewness(*it)); - statisticsResult->SetKurtosis(imageStatisticsFilter->GetKurtosis(*it)); - statisticsResult->SetRMS(std::sqrt(std::pow(imageStatisticsFilter->GetMean(*it), 2.) + imageStatisticsFilter->GetVariance(*it))); // variance = sigma^2 - statisticsResult->SetMPP(imageStatisticsFilter->GetMPP(*it)); - statisticsResult->SetLabel(*it); - - statisticsResult->SetEntropy(imageStatisticsFilter->GetEntropy(*it)); - statisticsResult->SetMedian(imageStatisticsFilter->GetMedian(*it)); - statisticsResult->SetUniformity(imageStatisticsFilter->GetUniformity(*it)); - statisticsResult->SetUPP(imageStatisticsFilter->GetUPP(*it)); - statisticsResult->SetHistogram(imageStatisticsFilter->GetHistogram(*it)); - - m_StatisticsByTimeStep[timeStep].push_back(statisticsResult); + statObj.m_N = imageStatisticsFilter->GetSum(*it) / (double) imageStatisticsFilter->GetMean(*it); + statObj.m_Volume = statObj.m_N * voxelVolume; + statObj.m_Mean = imageStatisticsFilter->GetMean(*it); + statObj.m_Min = imageStatisticsFilter->GetMinimum(*it); + statObj.m_Max = imageStatisticsFilter->GetMaximum(*it); + statObj.m_Std = imageStatisticsFilter->GetSigma(*it); + statObj.m_Skewness = imageStatisticsFilter->GetSkewness(*it); + statObj.m_Kurtosis = imageStatisticsFilter->GetKurtosis(*it); + statObj.m_RMS = std::sqrt(std::pow(imageStatisticsFilter->GetMean(*it), 2.) + imageStatisticsFilter->GetVariance(*it)); // variance = sigma^2 + statObj.m_MPP = imageStatisticsFilter->GetMPP(*it); + statObj.m_Label = *it; + + statObj.m_Entropy = imageStatisticsFilter->GetEntropy(*it); + statObj.m_Median = imageStatisticsFilter->GetMedian(*it); + statObj.m_Uniformity = imageStatisticsFilter->GetUniformity(*it); + statObj.m_UPP = imageStatisticsFilter->GetUPP(*it); + statObj.m_Histogram = imageStatisticsFilter->GetHistogram(*it); + + m_StatisticContainers.at(*it)->SetStatisticsForTimeStep(timeStep, statObj); + //m_StatisticsByTimeStep[timeStep].push_back(statisticsResult); ++it; } // swap maskGenerators back if (swapMasks) { m_SecondaryMask = m_InternalMask; m_InternalMask = nullptr; } } bool ImageStatisticsCalculator::IsUpdateRequired(unsigned int timeStep) const { unsigned long thisClassTimeStamp = this->GetMTime(); unsigned long inputImageTimeStamp = m_Image->GetMTime(); unsigned long statisticsTimeStamp = m_StatisticsUpdateTimePerTimeStep[timeStep]; if (thisClassTimeStamp > statisticsTimeStamp) // inputs have changed { return true; } if (inputImageTimeStamp > statisticsTimeStamp) // image has changed { return true; } if (m_MaskGenerator.IsNotNull()) { unsigned long maskGeneratorTimeStamp = m_MaskGenerator->GetMTime(); if (maskGeneratorTimeStamp > statisticsTimeStamp) // there is a mask generator and it has changed { return true; } } if (m_SecondaryMaskGenerator.IsNotNull()) { unsigned long maskGeneratorTimeStamp = m_SecondaryMaskGenerator->GetMTime(); if (maskGeneratorTimeStamp > statisticsTimeStamp) // there is a secondary mask generator and it has changed { return true; } } return false; } } diff --git a/Modules/ImageStatistics/mitkImageStatisticsCalculator.h b/Modules/ImageStatistics/mitkImageStatisticsCalculator.h index c1743253c1..469c7202bf 100644 --- a/Modules/ImageStatistics/mitkImageStatisticsCalculator.h +++ b/Modules/ImageStatistics/mitkImageStatisticsCalculator.h @@ -1,114 +1,114 @@ #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::Pointer GetStatistics(unsigned int timeStep=0, unsigned int label=1); + 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, - unsigned int timeStep); + typename itk::Image< TPixel, VImageDimension >* image, TimeGeometry* timeGeometry, TimeStepType timeStep); template < typename TPixel, unsigned int VImageDimension > void InternalCalculateStatisticsMasked( - typename itk::Image< TPixel, VImageDimension >* image, + 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_StatisticsByTimeStep; + std::vector m_StatisticContainers; std::vector m_StatisticsUpdateTimePerTimeStep; }; } #endif // MITKIMAGESTATISTICSCALCULATOR diff --git a/Modules/ImageStatistics/mitkImageStatisticsContainer.cpp b/Modules/ImageStatistics/mitkImageStatisticsContainer.cpp index 3441c6c5cb..04137b1574 100644 --- a/Modules/ImageStatistics/mitkImageStatisticsContainer.cpp +++ b/Modules/ImageStatistics/mitkImageStatisticsContainer.cpp @@ -1,181 +1,126 @@ /*=================================================================== 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 namespace mitk { - StatisticsContainer::StatisticsContainer(): - m_N(0), - m_Volume(nan("")), - m_Mean(nan("")), - m_Min(nan("")), - m_Max(nan("")), - m_Std(nan("")), - m_Skewness(nan("")), - m_Kurtosis(nan("")), - m_RMS(nan("")), - m_MPP(nan("")), - m_Median(nan("")), - m_Uniformity(nan("")), - m_UPP(nan("")), - m_Entropy(nan("")), - m_Label(0) + StatisticsContainer::StatisticsContainer() { - m_MinIndex.set_size(0); - m_MaxIndex.set_size(0); + this->Reset(); } - StatisticsContainer::RealType StatisticsContainer::GetVariance() const - { - return m_Std * m_Std; + bool StatisticsContainer::TimeStepExists(TimeStepType timeStep) const { + return m_TimeStepMap.find(timeStep) != m_TimeStepMap.end(); } - void StatisticsContainer::SetHistogram(HistogramType::Pointer hist) - { - if (m_Histogram != hist) - { - m_Histogram = hist; + StatisticsContainer::StatisticsObject StatisticsContainer::GetStatisticsForTimeStep(TimeStepType timeStep) const { + auto it = m_TimeStepMap.find(timeStep); + if (it != m_TimeStepMap.end()) { + return it->second; + } + mitkThrow() << "StatisticsObject for timeStep " << timeStep << " not found!"; + } + + void StatisticsContainer::SetStatisticsForTimeStep(TimeStepType timeStep, StatisticsObject statistics) { + if (timeStep < this->GetTimeSteps()) { + m_TimeStepMap[timeStep] = statistics; + } + else { + mitkThrow() << "Given timeStep " << timeStep << " out of timeStep geometry bounds. TimeSteps in geometry: " << this->GetTimeSteps(); } } void StatisticsContainer::PrintSelf(std::ostream &os, itk::Indent indent) const { Superclass::PrintSelf(os, indent); - auto statistics = GetStatisticsAsOrderedVector(); - os << std::endl << indent << "Statistics instance:"; - for (const auto& aStatisticValue : statistics) { - os << std::endl << indent.GetNextIndent() << aStatisticValue.first << ": " << aStatisticValue.second; + + for (unsigned int i = 0; i < this->GetTimeSteps(); i++) { + auto statistics = GetStatisticsAsOrderedVector(i); + os << std::endl << indent << "Statistics instance for timeStep " << i << ":"; + for (const auto& aStatisticValue : statistics) { + os << std::endl << indent.GetNextIndent() << aStatisticValue.first << ": " << aStatisticValue.second; + } } } - StatisticsContainer::statisticsOrderedVectorType StatisticsContainer::GetStatisticsAsOrderedVector() const + StatisticsContainer::statisticsOrderedVectorType StatisticsContainer::GetStatisticsAsOrderedVector(TimeStepType timeStep) const { - statisticsOrderedVectorType statisticsAsVector; - statisticsAsVector.emplace_back(std::make_pair("Mean", std::to_string(m_Mean))); - statisticsAsVector.emplace_back(std::make_pair("Median", std::to_string(m_Median))); - statisticsAsVector.emplace_back(std::make_pair("StandardDeviation", std::to_string(m_Std))); - statisticsAsVector.emplace_back(std::make_pair("RMS", std::to_string(m_RMS))); - statisticsAsVector.emplace_back(std::make_pair("Max", std::to_string( m_Max))); - statisticsAsVector.emplace_back(std::make_pair("MaxPosition", convertToString(m_MaxIndex))); - statisticsAsVector.emplace_back(std::make_pair("Min", std::to_string(m_Min))); - statisticsAsVector.emplace_back(std::make_pair("MinPosition", convertToString(m_MinIndex))); - statisticsAsVector.emplace_back(std::make_pair("#Voxel", std::to_string(m_N))); - statisticsAsVector.emplace_back(std::make_pair("Volume [mm^3]", std::to_string(m_Volume))); - statisticsAsVector.emplace_back(std::make_pair("Skewness", std::to_string(m_Skewness))); - statisticsAsVector.emplace_back(std::make_pair("Kurtosis", std::to_string(m_Kurtosis))); - statisticsAsVector.emplace_back(std::make_pair("Uniformity", std::to_string(m_Uniformity))); - statisticsAsVector.emplace_back(std::make_pair("Entropy", std::to_string(m_Entropy))); - statisticsAsVector.emplace_back(std::make_pair("MPP", std::to_string(m_MPP))); - statisticsAsVector.emplace_back(std::make_pair("UPP", std::to_string(m_UPP))); - return statisticsAsVector; + return GetStatisticsForTimeStep(timeStep).GetStatisticsAsOrderedVector(); } + unsigned int StatisticsContainer::GetNumberOfTimeSteps()const { + return this->GetTimeSteps(); + } void StatisticsContainer::Reset() { - m_N = 0; - m_Volume = nan(""); - m_Mean = nan(""); - m_Min = nan(""); - m_Max = nan(""); - m_Std = nan(""); - m_Skewness = nan(""); - m_Kurtosis = nan(""); - m_RMS = nan(""); - m_MPP = nan(""); - m_Median = nan(""); - m_Uniformity = nan(""); - m_UPP = nan(""); - m_Entropy = nan(""); - m_Histogram = HistogramType::New(); - m_MinIndex.set_size(0); - m_MaxIndex.set_size(0); - m_Label = 0; + for (auto iter = m_TimeStepMap.begin(); iter != m_TimeStepMap.end(); iter++) { + iter->second.Reset(); + } } itk::LightObject::Pointer StatisticsContainer::InternalClone() const { itk::LightObject::Pointer ioPtr = Superclass::InternalClone(); Self::Pointer rval = dynamic_cast(ioPtr.GetPointer()); if (rval.IsNull()) { itkExceptionMacro(<< "downcast to type " << "StatisticsContainer" << " failed."); } - rval->SetEntropy(this->GetEntropy()); - rval->SetKurtosis(this->GetKurtosis()); - rval->SetLabel(this->GetLabel()); - rval->SetMax(this->GetMax()); - rval->SetMin(this->GetMin()); - rval->SetMean(this->GetMean()); - rval->SetMedian(this->GetMedian()); - rval->SetMPP(this->GetMPP()); - rval->SetN(this->GetN()); - rval->SetVolume(this->GetVolume()); - rval->SetRMS(this->GetRMS()); - rval->SetSkewness(this->GetSkewness()); - rval->SetStd(this->GetStd()); - rval->SetUniformity(this->GetUniformity()); - rval->SetUPP(this->GetUPP()); - rval->SetHistogram(this->GetHistogram()); - rval->SetMinIndex(this->GetMinIndex()); - rval->SetMaxIndex(this->GetMaxIndex()); + // hopefully this is a deep copy (more sure than before) + rval->SetTimeStepMap(m_TimeStepMap); return ioPtr; } - std::string StatisticsContainer::convertToString(const vnl_vector& index) const - { - std::string result; - for (const auto& aValue : index) { - if (!result.empty()) { - result += "/"; - } - result += std::to_string(aValue); - } - return result; + void StatisticsContainer::SetTimeStepMap(TimeStepMapType map) { + m_TimeStepMap = map; } void StatisticsContainer::Print() { - auto statistics = this->GetStatisticsAsOrderedVector(); + + for (unsigned int i = 0; i < this->GetTimeSteps(); i++) { + auto statistics = this->GetStatisticsAsOrderedVector(i); // print all map key value pairs // const auto& val:statMap for (auto it = statistics.begin(); it != statistics.end(); ++it) { std::cout << it->first << ": " << it->second << std::endl; } // print the min and max index std::cout << "Min Index:" << std::endl; - for (auto it = this->GetMinIndex().begin(); it != this->GetMinIndex().end(); ++it) + for (auto it = m_TimeStepMap[i].m_MinIndex.begin(); it != m_TimeStepMap[i].m_MinIndex.end(); ++it) { std::cout << *it << " "; } std::cout << std::endl; // print the min and max index std::cout << "Max Index:" << std::endl; - for (auto it = this->GetMaxIndex().begin(); it != this->GetMaxIndex().end(); ++it) + for (auto it = m_TimeStepMap[i].m_MaxIndex.begin(); it != m_TimeStepMap[i].m_MaxIndex.end(); ++it) { std::cout << *it << " "; } std::cout << std::endl; + } } } diff --git a/Modules/ImageStatistics/mitkImageStatisticsContainer.h b/Modules/ImageStatistics/mitkImageStatisticsContainer.h index 19fa778f6f..0321a25b3d 100644 --- a/Modules/ImageStatistics/mitkImageStatisticsContainer.h +++ b/Modules/ImageStatistics/mitkImageStatisticsContainer.h @@ -1,161 +1,187 @@ /*=================================================================== 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 MITKIMAGESTATISTICSCONTAINER #define MITKIMAGESTATISTICSCONTAINER #include #include #include namespace mitk { /**Documentation @brief Container class for storing the computed image statistics. Stored statistics are: - N: number of voxels - Mean - MPP (Mean of positive pixels) - Median - Skewness - Kurtosis - Uniformity - UPP (Uniformity of positive pixels) - Std (Standard Deviation) - Min - Max - RMS (Root Mean Square) - Label (if applicable, the label (unsigned short) of the mask the statistics belong to) - Entropy - MinIndex (Index of Image where the Minimum is located) - MaxIndex (Index of Image where the Maximum is located) - Histogram of Pixel Values*/ class MITKIMAGESTATISTICS_EXPORT StatisticsContainer : public mitk::BaseData { public: mitkClassMacro(StatisticsContainer, mitk::BaseData) - /** Method for creation through the object factory. */ - itkNewMacro(Self) + /** Method for creation through the object factory. */ + itkNewMacro(Self) - typedef itk::Statistics::Histogram HistogramType; + + typedef itk::Statistics::Histogram HistogramType; typedef double RealType; typedef std::vector < std::pair > statisticsOrderedVectorType; virtual void SetRequestedRegionToLargestPossibleRegion() override {}; virtual bool RequestedRegionIsOutsideOfTheBufferedRegion() override { return false; }; virtual bool VerifyRequestedRegion() override { return true; }; virtual void SetRequestedRegion(const itk::DataObject*) override {}; + void SetTimeGeometryFromImage() {}; + + struct StatisticsObject { + + StatisticsObject() { + this->Reset(); + } + + long m_N; + RealType m_Volume; + RealType m_Mean, m_Min, m_Max, m_Std; + RealType m_Skewness; + RealType m_Kurtosis; + RealType m_RMS; + RealType m_MPP; + vnl_vector m_MinIndex, m_MaxIndex; + RealType m_Median; + RealType m_Uniformity; + RealType m_UPP; + RealType m_Entropy; + unsigned int m_Label; + + HistogramType::Pointer m_Histogram; + + RealType GetVariance() { + return m_Std * m_Std; + }; + + statisticsOrderedVectorType GetStatisticsAsOrderedVector() { + statisticsOrderedVectorType statisticsAsVector; + statisticsAsVector.emplace_back(std::make_pair("Mean", std::to_string(m_Mean))); + statisticsAsVector.emplace_back(std::make_pair("Median", std::to_string(m_Median))); + statisticsAsVector.emplace_back(std::make_pair("StandardDeviation", std::to_string(m_Std))); + statisticsAsVector.emplace_back(std::make_pair("RMS", std::to_string(m_RMS))); + statisticsAsVector.emplace_back(std::make_pair("Max", std::to_string(m_Max))); + statisticsAsVector.emplace_back(std::make_pair("MaxPosition", convertToString(m_MaxIndex))); + statisticsAsVector.emplace_back(std::make_pair("Min", std::to_string(m_Min))); + statisticsAsVector.emplace_back(std::make_pair("MinPosition", convertToString(m_MinIndex))); + statisticsAsVector.emplace_back(std::make_pair("#Voxel", std::to_string(m_N))); + statisticsAsVector.emplace_back(std::make_pair("Volume [mm^3]", std::to_string(m_Volume))); + statisticsAsVector.emplace_back(std::make_pair("Skewness", std::to_string(m_Skewness))); + statisticsAsVector.emplace_back(std::make_pair("Kurtosis", std::to_string(m_Kurtosis))); + statisticsAsVector.emplace_back(std::make_pair("Uniformity", std::to_string(m_Uniformity))); + statisticsAsVector.emplace_back(std::make_pair("Entropy", std::to_string(m_Entropy))); + statisticsAsVector.emplace_back(std::make_pair("MPP", std::to_string(m_MPP))); + statisticsAsVector.emplace_back(std::make_pair("UPP", std::to_string(m_UPP))); + return statisticsAsVector; + } + + std::string convertToString(const vnl_vector& index) { + std::string result; + for (const auto& aValue : index) { + if (!result.empty()) { + result += "/"; + } + result += std::to_string(aValue); + } + return result; + } + + void Reset() { + m_N = 0; + m_Volume = nan(""); + m_Mean = nan(""); + m_Min = nan(""); + m_Max = nan(""); + m_Std = nan(""); + m_Skewness = nan(""); + m_Kurtosis = nan(""); + m_RMS = nan(""); + m_MPP = nan(""); + m_Median = nan(""); + m_Uniformity = nan(""); + m_UPP = nan(""); + m_Entropy = nan(""); + + m_MinIndex.set_size(0); + m_MaxIndex.set_size(0); + m_Label = 0; + m_Histogram = HistogramType::New(); + }; + + + }; + + typedef std::map TimeStepMapType; + /**Documentation - @brief Returns a std::vector); - itkGetConstMacro(MinIndex, vnl_vector); - - itkSetMacro(MaxIndex, vnl_vector); - itkGetConstMacro(MaxIndex, vnl_vector); - - void SetHistogram(HistogramType::Pointer hist); - - itkGetConstMacro(Histogram, HistogramType::Pointer); - - itkSetMacro(Entropy, RealType); - itkGetConstMacro(Entropy, RealType); - - itkSetMacro(Median, RealType); - itkGetConstMacro(Median, RealType); - - itkSetMacro(Uniformity, RealType); - itkGetConstMacro(Uniformity, RealType); - - itkSetMacro(UPP, RealType); - itkGetConstMacro(UPP, RealType); - /**Documentation @brief Creates a StatisticsMapType containing all real valued statistics stored in this class (= all statistics except minIndex, maxIndex and the histogram) and prints its contents to std::cout*/ void Print(); + StatisticsObject GetStatisticsForTimeStep(TimeStepType timeStep) const; + void SetStatisticsForTimeStep(TimeStepType, StatisticsObject); + bool TimeStepExists(TimeStepType timeStep) const; + protected: StatisticsContainer(); virtual void PrintSelf(std::ostream &os, itk::Indent indent) const override; private: itk::LightObject::Pointer InternalClone() const override; std::string convertToString(const vnl_vector& index) const; - long m_N; - RealType m_Volume; - RealType m_Mean, m_Min, m_Max, m_Std; - RealType m_Skewness; - RealType m_Kurtosis; - RealType m_RMS; - RealType m_MPP; - vnl_vector m_MinIndex, m_MaxIndex; - RealType m_Median; - RealType m_Uniformity; - RealType m_UPP; - RealType m_Entropy; - unsigned int m_Label; - HistogramType::Pointer m_Histogram; + void SetTimeStepMap(TimeStepMapType map); + + TimeStepMapType m_TimeStepMap; }; } #endif // MITKIMAGESTATISTICSCONTAINER diff --git a/Modules/ImageStatistics/mitkIntensityProfile.cpp b/Modules/ImageStatistics/mitkIntensityProfile.cpp index d2fbb9eedf..2d6d89e301 100644 --- a/Modules/ImageStatistics/mitkIntensityProfile.cpp +++ b/Modules/ImageStatistics/mitkIntensityProfile.cpp @@ -1,380 +1,380 @@ /*=================================================================== 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 #include #include #include #include #include #include #include "mitkIntensityProfile.h" using namespace mitk; template static void ReadPixel(const PixelType&, Image::Pointer image, const itk::Index<3>& index, ScalarType* returnValue) { switch (image->GetDimension()) { case 2: { ImagePixelReadAccessor readAccess(image, image->GetSliceData(0)); *returnValue = readAccess.GetPixelByIndex(reinterpret_cast&>(index)); break; } case 3: { ImagePixelReadAccessor readAccess(image, image->GetVolumeData(0)); *returnValue = readAccess.GetPixelByIndex(index); break; } default: *returnValue = 0; break; } } static IntensityProfile::Pointer ComputeIntensityProfile(Image::Pointer image, itk::PolyLineParametricPath<3>::Pointer path) { IntensityProfile::Pointer intensityProfile = IntensityProfile::New(); itk::PolyLineParametricPath<3>::InputType input = path->StartOfInput(); BaseGeometry* imageGeometry = image->GetGeometry(); const PixelType pixelType = image->GetPixelType(); IntensityProfile::MeasurementVectorType measurementVector; itk::PolyLineParametricPath<3>::OffsetType offset; Point3D worldPoint; itk::Index<3> index; do { imageGeometry->IndexToWorld(path->Evaluate(input), worldPoint); imageGeometry->WorldToIndex(worldPoint, index); mitkPixelTypeMultiplex3(ReadPixel, pixelType, image, index, measurementVector.GetDataPointer()); intensityProfile->PushBack(measurementVector); offset = path->IncrementInput(input); } while ((offset[0] | offset[1] | offset[2]) != 0); return intensityProfile; } template static typename itk::InterpolateImageFunction::Pointer CreateInterpolateImageFunction(InterpolateImageFunction::Enum interpolator) { switch (interpolator) { case InterpolateImageFunction::NearestNeighbor: return itk::NearestNeighborInterpolateImageFunction::New().GetPointer(); case InterpolateImageFunction::Linear: return itk::LinearInterpolateImageFunction::New().GetPointer(); case InterpolateImageFunction::WindowedSinc_Blackman_3: return itk::WindowedSincInterpolateImageFunction >::New().GetPointer(); case InterpolateImageFunction::WindowedSinc_Blackman_4: return itk::WindowedSincInterpolateImageFunction >::New().GetPointer(); case InterpolateImageFunction::WindowedSinc_Blackman_5: return itk::WindowedSincInterpolateImageFunction >::New().GetPointer(); case InterpolateImageFunction::WindowedSinc_Cosine_3: return itk::WindowedSincInterpolateImageFunction >::New().GetPointer(); case InterpolateImageFunction::WindowedSinc_Cosine_4: return itk::WindowedSincInterpolateImageFunction >::New().GetPointer(); case InterpolateImageFunction::WindowedSinc_Cosine_5: return itk::WindowedSincInterpolateImageFunction >::New().GetPointer(); case InterpolateImageFunction::WindowedSinc_Hamming_3: return itk::WindowedSincInterpolateImageFunction >::New().GetPointer(); case InterpolateImageFunction::WindowedSinc_Hamming_4: return itk::WindowedSincInterpolateImageFunction >::New().GetPointer(); case InterpolateImageFunction::WindowedSinc_Hamming_5: return itk::WindowedSincInterpolateImageFunction >::New().GetPointer(); case InterpolateImageFunction::WindowedSinc_Lanczos_3: return itk::WindowedSincInterpolateImageFunction >::New().GetPointer(); case InterpolateImageFunction::WindowedSinc_Lanczos_4: return itk::WindowedSincInterpolateImageFunction >::New().GetPointer(); case InterpolateImageFunction::WindowedSinc_Lanczos_5: return itk::WindowedSincInterpolateImageFunction >::New().GetPointer(); case InterpolateImageFunction::WindowedSinc_Welch_3: return itk::WindowedSincInterpolateImageFunction >::New().GetPointer(); case InterpolateImageFunction::WindowedSinc_Welch_4: return itk::WindowedSincInterpolateImageFunction >::New().GetPointer(); case InterpolateImageFunction::WindowedSinc_Welch_5: return itk::WindowedSincInterpolateImageFunction >::New().GetPointer(); default: return itk::NearestNeighborInterpolateImageFunction::New().GetPointer(); } } template static void ComputeIntensityProfile(itk::Image* image, itk::PolyLineParametricPath<3>::Pointer path, unsigned int numSamples, InterpolateImageFunction::Enum interpolator, IntensityProfile::Pointer intensityProfile) { typename itk::InterpolateImageFunction >::Pointer interpolateImageFunction = CreateInterpolateImageFunction >(interpolator); interpolateImageFunction->SetInputImage(image); const itk::PolyLineParametricPath<3>::InputType startOfInput = path->StartOfInput(); const itk::PolyLineParametricPath<3>::InputType delta = 1.0 / (numSamples - 1); IntensityProfile::MeasurementVectorType measurementVector; for (unsigned int i = 0; i < numSamples; ++i) { measurementVector[0] = interpolateImageFunction->EvaluateAtContinuousIndex(path->Evaluate(startOfInput + i * delta)); intensityProfile->PushBack(measurementVector); } } static IntensityProfile::Pointer ComputeIntensityProfile(Image::Pointer image, itk::PolyLineParametricPath<3>::Pointer path, unsigned int numSamples, InterpolateImageFunction::Enum interpolator) { IntensityProfile::Pointer intensityProfile = IntensityProfile::New(); AccessFixedDimensionByItk_n(image, ComputeIntensityProfile, 3, (path, numSamples, interpolator, intensityProfile)); return intensityProfile; } class AddPolyLineElementToPath { public: AddPolyLineElementToPath(const PlaneGeometry* planarFigureGeometry, const BaseGeometry* imageGeometry, itk::PolyLineParametricPath<3>::Pointer path) : m_PlanarFigureGeometry(planarFigureGeometry), m_ImageGeometry(imageGeometry), m_Path(path) { } void operator()(const PlanarFigure::PolyLineElement& polyLineElement) { m_PlanarFigureGeometry->Map(polyLineElement, m_WorldPoint); m_ImageGeometry->WorldToIndex(m_WorldPoint, m_ContinuousIndexPoint); m_Vertex.CastFrom(m_ContinuousIndexPoint); m_Path->AddVertex(m_Vertex); } private: const PlaneGeometry* m_PlanarFigureGeometry; const BaseGeometry* m_ImageGeometry; itk::PolyLineParametricPath<3>::Pointer m_Path; Point3D m_WorldPoint; Point3D m_ContinuousIndexPoint; itk::PolyLineParametricPath<3>::ContinuousIndexType m_Vertex; }; static itk::PolyLineParametricPath<3>::Pointer CreatePathFromPlanarFigure(BaseGeometry* imageGeometry, PlanarFigure* planarFigure) { itk::PolyLineParametricPath<3>::Pointer path = itk::PolyLineParametricPath<3>::New(); const PlanarFigure::PolyLineType polyLine = planarFigure->GetPolyLine(0); std::for_each(polyLine.begin(), polyLine.end(), AddPolyLineElementToPath(planarFigure->GetPlaneGeometry(), imageGeometry, path)); return path; } static void AddPointToPath(const BaseGeometry* imageGeometry, const Point3D& point, itk::PolyLineParametricPath<3>::Pointer path) { Point3D continuousIndexPoint; imageGeometry->WorldToIndex(point, continuousIndexPoint); itk::PolyLineParametricPath<3>::ContinuousIndexType vertex; vertex.CastFrom(continuousIndexPoint); path->AddVertex(vertex); } static itk::PolyLineParametricPath<3>::Pointer CreatePathFromPoints(BaseGeometry* imageGeometry, const Point3D& startPoint, const Point3D& endPoint) { itk::PolyLineParametricPath<3>::Pointer path = itk::PolyLineParametricPath<3>::New(); AddPointToPath(imageGeometry, startPoint, path); AddPointToPath(imageGeometry, endPoint, path); return path; } IntensityProfile::Pointer mitk::ComputeIntensityProfile(Image::Pointer image, PlanarFigure::Pointer planarFigure) { return ::ComputeIntensityProfile(image, CreatePathFromPlanarFigure(image->GetGeometry(), planarFigure)); } IntensityProfile::Pointer mitk::ComputeIntensityProfile(Image::Pointer image, PlanarLine::Pointer planarLine, unsigned int numSamples, InterpolateImageFunction::Enum interpolator) { return ::ComputeIntensityProfile(image, CreatePathFromPlanarFigure(image->GetGeometry(), planarLine.GetPointer()), numSamples, interpolator); } IntensityProfile::Pointer mitk::ComputeIntensityProfile(Image::Pointer image, const Point3D& startPoint, const Point3D& endPoint, unsigned int numSamples, InterpolateImageFunction::Enum interpolator) { return ::ComputeIntensityProfile(image, CreatePathFromPoints(image->GetGeometry(), startPoint, endPoint), numSamples, interpolator); } IntensityProfile::InstanceIdentifier mitk::ComputeGlobalMaximum(IntensityProfile::ConstPointer intensityProfile, IntensityProfile::MeasurementType &max) { max = -vcl_numeric_limits::min(); IntensityProfile::InstanceIdentifier maxIndex = 0; IntensityProfile::ConstIterator end = intensityProfile->End(); IntensityProfile::MeasurementType measurement; for (IntensityProfile::ConstIterator it = intensityProfile->Begin(); it != end; ++it) { measurement = it.GetMeasurementVector()[0]; if (measurement > max) { max = measurement; maxIndex = it.GetInstanceIdentifier(); } } return maxIndex; } IntensityProfile::InstanceIdentifier mitk::ComputeGlobalMinimum(IntensityProfile::ConstPointer intensityProfile, IntensityProfile::MeasurementType &min) { min = vcl_numeric_limits::max(); IntensityProfile::InstanceIdentifier minIndex = 0; IntensityProfile::ConstIterator end = intensityProfile->End(); IntensityProfile::MeasurementType measurement; for (IntensityProfile::ConstIterator it = intensityProfile->Begin(); it != end; ++it) { measurement = it.GetMeasurementVector()[0]; if (measurement < min) { min = measurement; minIndex = it.GetInstanceIdentifier(); } } return minIndex; } IntensityProfile::InstanceIdentifier mitk::ComputeCenterOfMaximumArea(IntensityProfile::ConstPointer intensityProfile, IntensityProfile::InstanceIdentifier radius) { //const IntensityProfile::MeasurementType min = intensityProfile->GetMeasurementVector(ComputeGlobalMinimum(intensityProfile))[0]; IntensityProfile::MeasurementType min; ComputeGlobalMinimum(intensityProfile, min); const IntensityProfile::InstanceIdentifier areaWidth = 1 + 2 * radius; IntensityProfile::MeasurementType maxArea = 0; for (IntensityProfile::InstanceIdentifier i = 0; i < areaWidth; ++i) maxArea += intensityProfile->GetMeasurementVector(i)[0] - min; const IntensityProfile::InstanceIdentifier lastIndex = intensityProfile->Size() - areaWidth; IntensityProfile::InstanceIdentifier centerOfMaxArea = radius; IntensityProfile::MeasurementType area = maxArea; for (IntensityProfile::InstanceIdentifier i = 1; i <= lastIndex; ++i) { area += intensityProfile->GetMeasurementVector(i + areaWidth - 1)[0] - min; area -= intensityProfile->GetMeasurementVector(i - 1)[0] - min; if (area > maxArea) { maxArea = area; centerOfMaxArea = i + radius; // TODO: If multiple areas in the neighborhood have the same intensity chose the middle one instead of the first one. } } return centerOfMaxArea; } std::vector mitk::CreateVectorFromIntensityProfile(IntensityProfile::ConstPointer intensityProfile) { std::vector result; result.reserve(intensityProfile->Size()); IntensityProfile::ConstIterator end = intensityProfile->End(); for (IntensityProfile::ConstIterator it = intensityProfile->Begin(); it != end; ++it) result.push_back(it.GetMeasurementVector()[0]); return result; } IntensityProfile::Pointer mitk::CreateIntensityProfileFromVector(const std::vector& vector) { const IntensityProfile::InstanceIdentifier size = vector.size(); IntensityProfile::Pointer result = IntensityProfile::New(); result->Resize(size); for (IntensityProfile::InstanceIdentifier i = 0; i < size; ++i) result->SetMeasurement(i, 0, vector[i]); return result; } -void mitk::ComputeIntensityProfileStatistics(IntensityProfile::ConstPointer intensityProfile, StatisticsContainer::Pointer stats) +void mitk::ComputeIntensityProfileStatistics(IntensityProfile::ConstPointer intensityProfile, StatisticsContainer::StatisticsObject& stats) { typedef std::vector StatsVecType; StatsVecType statsVec = mitk::CreateVectorFromIntensityProfile( intensityProfile ); IntensityProfile::MeasurementType min; IntensityProfile::MeasurementType max; mitk::ComputeGlobalMinimum( intensityProfile, min ); mitk::ComputeGlobalMaximum( intensityProfile, max ); StatsVecType::size_type numSamples = statsVec.size(); double mean = 0.0; double rms = 0.0; for ( StatsVecType::const_iterator it = statsVec.begin(); it != statsVec.end(); ++it ) { double val = *it; mean += val; rms += val*val; } mean /= numSamples; rms /= numSamples; double var = 0.0; for ( StatsVecType::const_iterator it = statsVec.begin(); it != statsVec.end(); ++it ) { double diff = *it - mean; var += diff*diff; } var /= ( numSamples - 1 ); rms = sqrt( rms ); - stats->SetMin( static_cast( min ) ); - stats->SetMax( static_cast( max ) ); - stats->SetN( numSamples ); - stats->SetMean( mean ); - stats->SetStd( sqrt(var) ); - stats->SetRMS( rms ); + stats.m_Median = static_cast( min ); + stats.m_Max = static_cast( max ); + stats.m_N = numSamples; + stats.m_Mean = mean; + stats.m_Std = sqrt(var); + stats.m_RMS = rms; } diff --git a/Modules/ImageStatistics/mitkIntensityProfile.h b/Modules/ImageStatistics/mitkIntensityProfile.h index cc8e3ce4e3..45003be47d 100644 --- a/Modules/ImageStatistics/mitkIntensityProfile.h +++ b/Modules/ImageStatistics/mitkIntensityProfile.h @@ -1,137 +1,137 @@ /*=================================================================== 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 mitkIntensityProfile_h #define mitkIntensityProfile_h #include #include #include #include #include namespace mitk { typedef itk::Statistics::ListSample::MeasurementVectorType> IntensityProfile; /** \brief Compute intensity profile of an image for each pixel along the first PolyLine of a given planar figure. * * \param[in] image A two or three-dimensional image which consists of single component pixels. * \param[in] planarFigure A planar figure from which the first PolyLine is used to evaluate the intensity profile. * * \return The computed intensity profile. */ MITKIMAGESTATISTICS_EXPORT IntensityProfile::Pointer ComputeIntensityProfile(Image::Pointer image, PlanarFigure::Pointer planarFigure); namespace InterpolateImageFunction { enum Enum { NearestNeighbor, Linear, WindowedSinc_Blackman_3, WindowedSinc_Blackman_4, WindowedSinc_Blackman_5, WindowedSinc_Cosine_3, WindowedSinc_Cosine_4, WindowedSinc_Cosine_5, WindowedSinc_Hamming_3, WindowedSinc_Hamming_4, WindowedSinc_Hamming_5, WindowedSinc_Lanczos_3, WindowedSinc_Lanczos_4, WindowedSinc_Lanczos_5, WindowedSinc_Welch_3, WindowedSinc_Welch_4, WindowedSinc_Welch_5 }; } /** \brief Compute intensity profile of an image for each sample along a planar line. * * \param[in] image A three-dimensional image which consists of single component pixels. * \param[in] planarLine A planar line along which the intensity profile will be evaluated. * \param[in] numSamples Number of samples along the planar line (must be at least 2). * \param[in] interpolator Image interpolation function which is used to read each sample. * * \return The computed intensity profile. */ MITKIMAGESTATISTICS_EXPORT IntensityProfile::Pointer ComputeIntensityProfile(Image::Pointer image, PlanarLine::Pointer planarLine, unsigned int numSamples, InterpolateImageFunction::Enum interpolator = InterpolateImageFunction::NearestNeighbor); /** \brief Compute intensity profile of an image for each sample between two points. * * \param[in] image A three-dimensional image which consists of single component pixels. * \param[in] startPoint A point at which the first sample is to be read. * \param[in] endPoint A point at which the last sample is to be read. * \param[in] numSamples Number of samples between startPoint and endPoint (must be at least 2). * \param[in] interpolator Image interpolation function which is used to read each sample. * * \return The computed intensity profile. */ MITKIMAGESTATISTICS_EXPORT IntensityProfile::Pointer ComputeIntensityProfile(Image::Pointer image, const Point3D& startPoint, const Point3D& endPoint, unsigned int numSamples, InterpolateImageFunction::Enum interpolator = InterpolateImageFunction::NearestNeighbor); /** \brief Compute global maximum of an intensity profile. * * \param[in] intensityProfile An intensity profile. * * \return Index of the global maximum. */ MITKIMAGESTATISTICS_EXPORT IntensityProfile::InstanceIdentifier ComputeGlobalMaximum(IntensityProfile::ConstPointer intensityProfile, IntensityProfile::MeasurementType &max); /** \brief Compute global minimum of an intensity profile. * * \param[in] intensityProfile An intensity profile. * * \return Index of the global minimum. */ MITKIMAGESTATISTICS_EXPORT IntensityProfile::InstanceIdentifier ComputeGlobalMinimum(IntensityProfile::ConstPointer intensityProfile, IntensityProfile::MeasurementType &min); /** \brief Compute statistics of an intensity profile. * * \param[in] intensityProfile An intensity profile. * * \param[in] stats An ImageStatisticsCalculator::Statistics object to hold the calculated statistics. * */ - MITKIMAGESTATISTICS_EXPORT void ComputeIntensityProfileStatistics(IntensityProfile::ConstPointer intensityProfile, StatisticsContainer::Pointer stats); + MITKIMAGESTATISTICS_EXPORT void ComputeIntensityProfileStatistics(IntensityProfile::ConstPointer intensityProfile, StatisticsContainer::StatisticsObject& stats); /** \brief Compute center of maximum area under the curve of an intensity profile. * * \param[in] intensityProfile An intensity profile. * \param[in] radius Radius of the area (width of area equals 1 + 2 * radius). * * \return Index of the maximum area center. */ MITKIMAGESTATISTICS_EXPORT IntensityProfile::InstanceIdentifier ComputeCenterOfMaximumArea(IntensityProfile::ConstPointer intensityProfile, IntensityProfile::InstanceIdentifier radius); /** \brief Convert an intensity profile to a standard library vector. * * \param[in] intensityProfile An intensity profile. * * \return Standard library vector which contains the input intensity profile measurements. */ MITKIMAGESTATISTICS_EXPORT std::vector CreateVectorFromIntensityProfile(IntensityProfile::ConstPointer intensityProfile); /** \brief Convert a standard library vector to an intensity profile. * * \param[in] vector An standard library vector which contains intensity profile measurements. * * \return An intensity profile. */ MITKIMAGESTATISTICS_EXPORT IntensityProfile::Pointer CreateIntensityProfileFromVector(const std::vector& vector); } #endif diff --git a/Modules/ImageStatistics/mitkPointSetDifferenceStatisticsCalculator.cpp b/Modules/ImageStatistics/mitkPointSetDifferenceStatisticsCalculator.cpp index 4b724d0788..4e032d4c45 100644 --- a/Modules/ImageStatistics/mitkPointSetDifferenceStatisticsCalculator.cpp +++ b/Modules/ImageStatistics/mitkPointSetDifferenceStatisticsCalculator.cpp @@ -1,215 +1,215 @@ /*=================================================================== 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 "mitkPointSetDifferenceStatisticsCalculator.h" mitk::PointSetDifferenceStatisticsCalculator::PointSetDifferenceStatisticsCalculator() : m_StatisticsCalculated(false) { - m_Statistics = StatisticsContainer::New(); + m_Statistics = StatisticsContainer::StatisticsObject(); m_PointSet1 = mitk::PointSet::New(); m_PointSet2 = mitk::PointSet::New(); //m_Statistics.Reset(); } mitk::PointSetDifferenceStatisticsCalculator::PointSetDifferenceStatisticsCalculator(mitk::PointSet::Pointer pSet1, mitk::PointSet::Pointer pSet2) { - m_Statistics = StatisticsContainer::New(); + m_Statistics = StatisticsContainer::StatisticsObject(); m_PointSet1 = pSet1; m_PointSet2 = pSet2; m_StatisticsCalculated = false; //m_Statistics.Reset(); } mitk::PointSetDifferenceStatisticsCalculator::~PointSetDifferenceStatisticsCalculator() { } void mitk::PointSetDifferenceStatisticsCalculator::SetPointSets(mitk::PointSet::Pointer pSet1, mitk::PointSet::Pointer pSet2) { if (pSet1.IsNotNull()) { m_PointSet1 = pSet1; } if (pSet2.IsNotNull()) { m_PointSet2 = pSet2; } m_StatisticsCalculated = false; //m_Statistics.Reset(); } std::vector mitk::PointSetDifferenceStatisticsCalculator::GetDifferences() { if (!m_StatisticsCalculated) { this->ComputeStatistics(); } return m_DifferencesVector; } std::vector mitk::PointSetDifferenceStatisticsCalculator::GetSquaredDifferences() { if (!m_StatisticsCalculated) { this->ComputeStatistics(); } return m_SquaredDifferencesVector; } double mitk::PointSetDifferenceStatisticsCalculator::GetMean() { if (!m_StatisticsCalculated) { this->ComputeStatistics(); } - return m_Statistics->GetMean(); + return m_Statistics.m_Mean; } double mitk::PointSetDifferenceStatisticsCalculator::GetSD() { if (!m_StatisticsCalculated) { this->ComputeStatistics(); } - return m_Statistics->GetStd(); + return m_Statistics.m_Std; } double mitk::PointSetDifferenceStatisticsCalculator::GetVariance() { if (!m_StatisticsCalculated) { this->ComputeStatistics(); } - return m_Statistics->GetVariance(); + return m_Statistics.GetVariance(); } double mitk::PointSetDifferenceStatisticsCalculator::GetRMS() { if (!m_StatisticsCalculated) { this->ComputeStatistics(); } - return m_Statistics->GetRMS(); + return m_Statistics.m_RMS; } double mitk::PointSetDifferenceStatisticsCalculator::GetMedian() { if (!m_StatisticsCalculated) { this->ComputeStatistics(); } - return m_Statistics->GetMedian(); + return m_Statistics.m_Median; } double mitk::PointSetDifferenceStatisticsCalculator::GetMax() { if (!m_StatisticsCalculated) { this->ComputeStatistics(); } - return m_Statistics->GetMax(); + return m_Statistics.m_Max; } double mitk::PointSetDifferenceStatisticsCalculator::GetMin() { if (!m_StatisticsCalculated) { this->ComputeStatistics(); } - return m_Statistics->GetMin(); + return m_Statistics.m_Min; } double mitk::PointSetDifferenceStatisticsCalculator::GetNumberOfPoints() { if (!m_StatisticsCalculated) { this->ComputeStatistics(); } - return m_Statistics->GetN(); + return m_Statistics.m_N; } void mitk::PointSetDifferenceStatisticsCalculator::ComputeStatistics() { if ((m_PointSet1.IsNull())||(m_PointSet2.IsNull())) { itkExceptionMacro("Point sets specified are not valid. Please specify correct Point sets"); } else if (m_PointSet1->GetSize()!=m_PointSet2->GetSize()) { itkExceptionMacro("PointSets are not equal. Please make sure that your PointSets have the same size and hold corresponding points."); } else if (m_PointSet1->GetSize()==0) { itkExceptionMacro("There are no points in the PointSets. Please make sure that the PointSets contain points"); } else { double mean = 0.0; double sd = 0.0; double rms= 0.0; std::vector differencesVector; mitk::Point3D point1; mitk::Point3D point2; int numberOfPoints = m_PointSet1->GetSize(); //Iterate over both pointsets in order to compare all points pair-wise mitk::PointSet::PointsIterator end = m_PointSet1->End(); for( mitk::PointSet::PointsIterator pointSetIterator = m_PointSet1->Begin(), pointSetIterator2 = m_PointSet2->Begin(); pointSetIterator != end; ++pointSetIterator, ++pointSetIterator2) //iterate simultaneously over both sets { point1 = pointSetIterator.Value(); point2 = pointSetIterator2.Value(); double squaredDistance = point1.SquaredEuclideanDistanceTo(point2); mean+=sqrt(squaredDistance); rms+=squaredDistance; this->m_SquaredDifferencesVector.push_back(squaredDistance); differencesVector.push_back(sqrt(squaredDistance)); } m_DifferencesVector = differencesVector; mean = mean/numberOfPoints; rms = sqrt(rms/numberOfPoints); for (std::vector::size_type i=0; iSetMean(mean); - m_Statistics->SetStd(sd); - m_Statistics->SetRMS(rms); - m_Statistics->SetMin(differencesVector.at(0)); - m_Statistics->SetMax(differencesVector.at(numberOfPoints-1)); - m_Statistics->SetMedian(median); - m_Statistics->SetN(numberOfPoints); + m_Statistics.m_Mean = mean; + m_Statistics.m_Std = sd; + m_Statistics.m_RMS = rms; + m_Statistics.m_Min = differencesVector.at(0); + m_Statistics.m_Max = differencesVector.at(numberOfPoints-1); + m_Statistics.m_Median = median; + m_Statistics.m_N = numberOfPoints; m_StatisticsCalculated = true; } } diff --git a/Modules/ImageStatistics/mitkPointSetDifferenceStatisticsCalculator.h b/Modules/ImageStatistics/mitkPointSetDifferenceStatisticsCalculator.h index 46e3635937..e37a997e89 100644 --- a/Modules/ImageStatistics/mitkPointSetDifferenceStatisticsCalculator.h +++ b/Modules/ImageStatistics/mitkPointSetDifferenceStatisticsCalculator.h @@ -1,110 +1,110 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef _MITK_PointSetDifferenceStatisticsCalculator_H #define _MITK_PointSetDifferenceStatisticsCalculator_H #include #include #include "mitkImageStatisticsCalculator.h" #include namespace mitk { /** * \brief Class for calculating the difference between two corresponding point sets. * The user can access the single distances between corresponding points as well as a complete statistic (mean, sd, rms, median, max, min) * The point sets must be of equal size! */ class MITKIMAGESTATISTICS_EXPORT PointSetDifferenceStatisticsCalculator : public itk::Object { public: mitkClassMacroItkParent( PointSetDifferenceStatisticsCalculator, itk::Object ); itkFactorylessNewMacro(Self) itkCloneMacro(Self) mitkNewMacro2Param(PointSetDifferenceStatisticsCalculator,mitk::PointSet::Pointer,mitk::PointSet::Pointer); /*! \brief set point sets to be compared */ void SetPointSets(mitk::PointSet::Pointer pSet1, mitk::PointSet::Pointer pSet2); /*! \brief returns a vector holding the differences between the corresponding points of the point sets */ std::vector GetDifferences(); /*! \brief returns a vector holding the squared differences between the corresponding points of the point sets */ std::vector GetSquaredDifferences(); /*! \brief returns the mean distance of all corresponding points of the point sets */ double GetMean(); /*! \brief returns the standard deviation of the distances between all corresponding points of the point sets */ double GetSD(); /*! \brief returns the variance of the distances between all corresponding points of the point sets */ double GetVariance(); /*! \brief returns the root mean squared distance of all corresponding points of the point sets */ double GetRMS(); /*! \brief returns the median distance of all corresponding points of the point sets */ double GetMedian(); /*! \brief returns the maximal distance of all corresponding points of the point sets */ double GetMax(); /*! \brief returns the minimal distance of all corresponding points of the point sets */ double GetMin(); /*! \brief returns the total number of corresponding points of the point sets */ double GetNumberOfPoints(); protected: PointSetDifferenceStatisticsCalculator(); PointSetDifferenceStatisticsCalculator(mitk::PointSet::Pointer,mitk::PointSet::Pointer); ~PointSetDifferenceStatisticsCalculator() override; /*! \brief Method for computing the complete statistics of the differences between the given point sets. */ void ComputeStatistics(); - mitk::StatisticsContainer::Pointer m_Statistics; ///< struct holding the statistics + mitk::StatisticsContainer::StatisticsObject m_Statistics; ///< struct holding the statistics std::vector m_DifferencesVector; ///< vector holding the differences between the corresponding points std::vector m_SquaredDifferencesVector; ///< vector holding the squared differences between the corresponding points mitk::PointSet::Pointer m_PointSet1; ///< first point set used for comparison mitk::PointSet::Pointer m_PointSet2; ///< second point set used for comparison bool m_StatisticsCalculated; ///< flag indicating whether statistics are already calculated or not. }; } #endif // #define _MITK_PointSetDifferenceStatisticsCalculator_H diff --git a/Modules/ImageStatisticsUI/Qmitk/QmitkImageStatisticsCalculationJob.cpp b/Modules/ImageStatisticsUI/Qmitk/QmitkImageStatisticsCalculationJob.cpp index 72692270d1..01283c5569 100644 --- a/Modules/ImageStatisticsUI/Qmitk/QmitkImageStatisticsCalculationJob.cpp +++ b/Modules/ImageStatisticsUI/Qmitk/QmitkImageStatisticsCalculationJob.cpp @@ -1,250 +1,264 @@ /*=================================================================== 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 "QmitkImageStatisticsCalculationJob.h" //QT headers #include #include #include #include #include QmitkImageStatisticsCalculationJob::QmitkImageStatisticsCalculationJob() : QThread() , m_StatisticsImage(nullptr) , m_BinaryMask(nullptr) , m_PlanarFigureMask(nullptr) , m_TimeStep(0) , m_IgnoreZeros(false) , m_HistogramNBins(100) , m_StatisticChanged(false) , m_CalculationSuccessful(false) { } QmitkImageStatisticsCalculationJob::~QmitkImageStatisticsCalculationJob() { } void QmitkImageStatisticsCalculationJob::Initialize( mitk::Image::ConstPointer image, mitk::Image::ConstPointer binaryImage, mitk::PlanarFigure::ConstPointer planarFig ) { // reset old values if( this->m_StatisticsImage.IsNotNull() ) this->m_StatisticsImage = nullptr; if( this->m_BinaryMask.IsNotNull() ) this->m_BinaryMask = nullptr; if( this->m_PlanarFigureMask.IsNotNull()) this->m_PlanarFigureMask = nullptr; // set new values if passed in if(image.IsNotNull()) this->m_StatisticsImage = image; if(binaryImage.IsNotNull()) this->m_BinaryMask = binaryImage; if(planarFig.IsNotNull()) this->m_PlanarFigureMask = planarFig; } void QmitkImageStatisticsCalculationJob::SetTimeStep( int times ) { this->m_TimeStep = times; } int QmitkImageStatisticsCalculationJob::GetTimeStep() const { return this->m_TimeStep; } -std::vector QmitkImageStatisticsCalculationJob::GetStatisticsData() const +mitk::StatisticsContainer::ConstPointer QmitkImageStatisticsCalculationJob::GetStatisticsData() const { - return this->m_StatisticsVector; + mitk::StatisticsContainer::ConstPointer constContainer = this->m_StatisticsContainer; + return constContainer; } mitk::Image::ConstPointer QmitkImageStatisticsCalculationJob::GetStatisticsImage() const { return this->m_StatisticsImage; } mitk::Image::ConstPointer QmitkImageStatisticsCalculationJob::GetMaskImage() const { return this->m_BinaryMask; } mitk::PlanarFigure::ConstPointer QmitkImageStatisticsCalculationJob::GetPlanarFigure() const { return this->m_PlanarFigureMask; } void QmitkImageStatisticsCalculationJob::SetIgnoreZeroValueVoxel(bool _arg) { this->m_IgnoreZeros = _arg; } bool QmitkImageStatisticsCalculationJob::GetIgnoreZeroValueVoxel() const { return this->m_IgnoreZeros; } void QmitkImageStatisticsCalculationJob::SetHistogramNBins(unsigned int nbins) { this->m_HistogramNBins = nbins; } unsigned int QmitkImageStatisticsCalculationJob::GetHistogramNBins() const { return this->m_HistogramNBins; } std::string QmitkImageStatisticsCalculationJob::GetLastErrorMessage() const { return m_message; } QmitkImageStatisticsCalculationJob::HistogramType::ConstPointer QmitkImageStatisticsCalculationJob::GetTimeStepHistogram(unsigned int t) const { if (t >= this->m_HistogramVector.size()) return nullptr; return this->m_HistogramVector[t]; } bool QmitkImageStatisticsCalculationJob::GetStatisticsChangedFlag() const { return m_StatisticChanged; } bool QmitkImageStatisticsCalculationJob::GetStatisticsUpdateSuccessFlag() const { return m_CalculationSuccessful; } void QmitkImageStatisticsCalculationJob::run() { bool statisticCalculationSuccessful = true; mitk::ImageStatisticsCalculator::Pointer calculator = mitk::ImageStatisticsCalculator::New(); if(this->m_StatisticsImage.IsNotNull()) { calculator->SetInputImage(m_StatisticsImage->Clone()); } else { statisticCalculationSuccessful = false; } // Bug 13416 : The ImageStatistics::SetImageMask() method can throw exceptions, i.e. when the dimensionality // of the masked and input image differ, we need to catch them and mark the calculation as failed // the same holds for the ::SetPlanarFigure() try { if(this->m_BinaryMask.IsNotNull()) { mitk::ImageMaskGenerator::Pointer imgMask = mitk::ImageMaskGenerator::New(); imgMask->SetImageMask(m_BinaryMask->Clone()); calculator->SetMask(imgMask.GetPointer()); } if(this->m_PlanarFigureMask.IsNotNull()) { mitk::PlanarFigureMaskGenerator::Pointer pfMaskGen = mitk::PlanarFigureMaskGenerator::New(); pfMaskGen->SetInputImage(m_StatisticsImage->Clone()); pfMaskGen->SetPlanarFigure(m_PlanarFigureMask->Clone()); calculator->SetMask(pfMaskGen.GetPointer()); } } catch (const mitk::Exception& e) { MITK_ERROR << "MITK Exception: " << e.what(); m_message = e.what(); statisticCalculationSuccessful = false; } catch( const itk::ExceptionObject& e) { MITK_ERROR << "ITK Exception:" << e.what(); m_message = e.what(); statisticCalculationSuccessful = false; } catch ( const std::runtime_error &e ) { MITK_ERROR<< "Runtime Exception: " << e.what(); m_message = e.what(); statisticCalculationSuccessful = false; } catch ( const std::exception &e ) { MITK_ERROR<< "Standard Exception: " << e.what(); m_message = e.what(); statisticCalculationSuccessful = false; } bool statisticChanged = false; if (this->m_IgnoreZeros) { mitk::IgnorePixelMaskGenerator::Pointer ignorePixelValueMaskGen = mitk::IgnorePixelMaskGenerator::New(); ignorePixelValueMaskGen->SetIgnoredPixelValue(0); ignorePixelValueMaskGen->SetInputImage(m_StatisticsImage->Clone()); calculator->SetSecondaryMask(ignorePixelValueMaskGen.GetPointer()); } else { calculator->SetSecondaryMask(nullptr); } calculator->SetNBinsForHistogramStatistics(m_HistogramNBins); for (unsigned int i = 0; i < m_StatisticsImage->GetTimeSteps(); i++) { try { calculator->GetStatistics(i); } catch ( mitk::Exception& e) { m_message = e.GetDescription(); MITK_ERROR<< "MITK Exception: " << e.what(); statisticCalculationSuccessful = false; } catch ( const std::runtime_error &e ) { m_message = "Failure: " + std::string(e.what()); MITK_ERROR<< "Runtime Exception: " << e.what(); statisticCalculationSuccessful = false; } catch ( const std::exception &e ) { m_message = "Failure: " + std::string(e.what()); MITK_ERROR<< "Standard Exception: " << e.what(); statisticCalculationSuccessful = false; } } this->m_StatisticChanged = statisticChanged; this->m_CalculationSuccessful = statisticCalculationSuccessful; if(statisticCalculationSuccessful) { - this->m_StatisticsVector.clear(); + //this->m_StatisticsVector.clear(); + //this->m_HistogramVector.clear(); + + //for (unsigned int i = 0; i < m_StatisticsImage->GetTimeSteps(); i++) + //{ + // this->m_StatisticsVector.push_back(calculator->GetStatistics(i).GetPointer()); + // this->m_HistogramVector.push_back((HistogramType*)this->m_StatisticsVector[i]->GetHistogram().GetPointer()); + //} + + // maybe use reset if this line does not work + m_StatisticsContainer = mitk::StatisticsContainer::New(); this->m_HistogramVector.clear(); + mitk::TimeGeometry::Pointer tg = m_StatisticsImage->GetTimeGeometry()->Clone(); + m_StatisticsContainer->SetTimeGeometry(tg); for (unsigned int i = 0; i < m_StatisticsImage->GetTimeSteps(); i++) { - this->m_StatisticsVector.push_back(calculator->GetStatistics(i).GetPointer()); - this->m_HistogramVector.push_back((HistogramType*)this->m_StatisticsVector[i]->GetHistogram().GetPointer()); + this->m_StatisticsContainer->SetStatisticsForTimeStep(i, calculator->GetStatistics(i)); + this->m_HistogramVector.push_back((HistogramType*) calculator->GetStatistics(i).m_Histogram.GetPointer()); } + } } diff --git a/Modules/ImageStatisticsUI/Qmitk/QmitkImageStatisticsCalculationJob.h b/Modules/ImageStatisticsUI/Qmitk/QmitkImageStatisticsCalculationJob.h index 7b87469157..b636d3d80f 100644 --- a/Modules/ImageStatisticsUI/Qmitk/QmitkImageStatisticsCalculationJob.h +++ b/Modules/ImageStatisticsUI/Qmitk/QmitkImageStatisticsCalculationJob.h @@ -1,114 +1,115 @@ /*=================================================================== 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 QMITKIMAGESTATISTICSCALCULATIONTHREAD_H_INCLUDED #define QMITKIMAGESTATISTICSCALCULATIONTHREAD_H_INCLUDED //QT headers #include #include //mitk headers #include "mitkImage.h" #include "mitkPlanarFigure.h" #include "mitkImageStatisticsCalculator.h" #include // itk headers #ifndef __itkHistogram_h #include #endif /** /brief This class is executed as background thread for image statistics calculation. * Documentation: This class is derived from QThread and is intended to be used by QmitkImageStatisticsView to run the image statistics calculation in a background thread keepung the gui usable. */ class MITKIMAGESTATISTICSUI_EXPORT QmitkImageStatisticsCalculationJob : public QThread { Q_OBJECT public: typedef itk::Statistics::Histogram HistogramType; /*! /brief standard constructor. */ QmitkImageStatisticsCalculationJob(); /*! /brief standard destructor. */ ~QmitkImageStatisticsCalculationJob(); /*! /brief Initializes the object with necessary data. */ void Initialize( mitk::Image::ConstPointer image, mitk::Image::ConstPointer binaryImage, mitk::PlanarFigure::ConstPointer planarFig ); /*! /brief returns the calculated image statistics. */ - std::vector GetStatisticsData() const; + mitk::StatisticsContainer::ConstPointer GetStatisticsData() const; mitk::Image::ConstPointer GetStatisticsImage() const; mitk::Image::ConstPointer GetMaskImage() const; mitk::PlanarFigure::ConstPointer GetPlanarFigure() const; /*! /brief Set the time step of the image you want to process. */ void SetTimeStep( int times ); /*! /brief Get the time step of the image you want to process. */ int GetTimeStep() const; /*! /brief Set flag to ignore zero valued voxels */ void SetIgnoreZeroValueVoxel( bool _arg ); /*! /brief Get status of zero value voxel ignoring. */ bool GetIgnoreZeroValueVoxel() const; /*! /brief Set bin size for histogram resolution.*/ void SetHistogramNBins( unsigned int nbins); /*! /brief Get bin size for histogram resolution.*/ unsigned int GetHistogramNBins() const; /*! /brief Returns the histogram of the currently selected time step. */ HistogramType::ConstPointer GetTimeStepHistogram(unsigned int t = 0) const; /*! /brief Returns a flag indicating if the statistics have changed during calculation */ bool GetStatisticsChangedFlag() const; /*! /brief Returns a flag the indicates if the statistics are updated successfully */ bool GetStatisticsUpdateSuccessFlag() const; /*! /brief Method called once the thread is executed. */ void run() override; std::string GetLastErrorMessage() const; private: //member declaration mitk::Image::ConstPointer m_StatisticsImage; ///< member variable holds the input image for which the statistics need to be calculated. mitk::Image::ConstPointer m_BinaryMask; ///< member variable holds the binary mask image for segmentation image statistics calculation. mitk::PlanarFigure::ConstPointer m_PlanarFigureMask; ///< member variable holds the planar figure for segmentation image statistics calculation. - std::vector m_StatisticsVector; ///< member variable holds the result structs. + //std::vector m_StatisticsVector; ///< member variable holds the result structs. + mitk::StatisticsContainer::Pointer m_StatisticsContainer; int m_TimeStep; ///< member variable holds the time step for statistics calculation bool m_IgnoreZeros; ///< member variable holds flag to indicate if zero valued voxel should be suppressed unsigned int m_HistogramNBins; ///< member variable holds the bin size for histogram resolution. bool m_StatisticChanged; ///< flag set if statistics have changed bool m_CalculationSuccessful; ///< flag set if statistics calculation was successful std::vector m_HistogramVector; ///< member holds the histograms of all time steps. std::string m_message; }; #endif // QMITKIMAGESTATISTICSCALCULATIONTHREAD_H_INCLUDED diff --git a/Modules/ImageStatisticsUI/Qmitk/QmitkImageStatisticsTableModel.cpp b/Modules/ImageStatisticsUI/Qmitk/QmitkImageStatisticsTableModel.cpp index 157a19c4fd..36a9d31bd3 100644 --- a/Modules/ImageStatisticsUI/Qmitk/QmitkImageStatisticsTableModel.cpp +++ b/Modules/ImageStatisticsUI/Qmitk/QmitkImageStatisticsTableModel.cpp @@ -1,166 +1,161 @@ /*=================================================================== 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 "QmitkImageStatisticsTableModel.h" +#include QmitkImageStatisticsTableModel::QmitkImageStatisticsTableModel(QObject *parent) : -QAbstractTableModel(parent) + QAbstractTableModel(parent) {} int QmitkImageStatisticsTableModel::rowCount(const QModelIndex &parent) const { - if (parent.isValid()) - { - return 0; - } - - if (!m_statistics.empty()){ - return static_cast(m_statistics.front()->GetStatisticsAsOrderedVector().size()); - } - else - { - return 0; - } + if (parent.isValid()) + { + return 0; + } + mitk::StatisticsContainer::StatisticsObject obj; + return static_cast(obj.GetStatisticsAsOrderedVector().size()); } int QmitkImageStatisticsTableModel::columnCount(const QModelIndex &parent) const { - if (parent.isValid()) - return 0; + if (parent.isValid() || m_statistics.IsNull()) + return 0; - return static_cast(m_statistics.size()); + return static_cast(m_statistics->GetNumberOfTimeSteps()); } QVariant QmitkImageStatisticsTableModel::data(const QModelIndex &index, int role) const { - if (!index.isValid()) - return QVariant(); + if (!index.isValid()) + return QVariant(); - QVariant result; - if (m_viewMode == viewMode::imageXStatistic) { - if (!m_statistics.empty() && index.row() < m_statistics.front()->GetStatisticsAsOrderedVector().size() && index.column() < m_statistics.size()) - { - if (Qt::DisplayRole == role) - { - auto statistics = m_statistics.at(index.column()); - auto statisticsVector = statistics->GetStatisticsAsOrderedVector(); - auto statisticsValueString = statisticsVector.at(index.row()).second; - result = QVariant(QString::fromStdString(statisticsValueString)); - } - else if (Qt::UserRole == role) - { - result = QVariant(index.row()); - } + QVariant result; + if (m_viewMode == viewMode::imageXStatistic) { + if (m_statistics.IsNotNull() && index.row() < rowCount(QModelIndex()) && index.column() < columnCount(QModelIndex())) + { + if (Qt::DisplayRole == role) + { + auto statisticsVector = m_statistics->GetStatisticsAsOrderedVector(index.column()); + auto statisticsValueString = statisticsVector.at(index.row()).second; + result = QVariant(QString::fromStdString(statisticsValueString)); + } + else if (Qt::UserRole == role) + { + result = QVariant(index.row()); } + } + } - return result; + return result; } Qt::ItemFlags QmitkImageStatisticsTableModel::flags(const QModelIndex &index) const { - Qt::ItemFlags flags = QAbstractItemModel::flags(index); + Qt::ItemFlags flags = QAbstractItemModel::flags(index); - return flags; + return flags; } QVariant QmitkImageStatisticsTableModel::headerData(int section, Qt::Orientation orientation, int role) const { - if ((Qt::DisplayRole == role) && (Qt::Horizontal == orientation)) - { - if (!m_imageNodes.empty()) { - if (m_viewMode == viewMode::imageXStatistic) { - std::string maskName; - if (!m_maskNodes.empty() && m_maskNodes.size() == m_imageNodes.size()) { - maskName = " / " + m_maskNodes.at(section)->GetName(); - } - return QVariant((m_imageNodes.at(section)->GetName() + maskName + std::to_string(section)).c_str()); + if ((Qt::DisplayRole == role) && (Qt::Horizontal == orientation)) + { + if (!m_imageNodes.empty()) { + if (m_viewMode == viewMode::imageXStatistic) { + std::string maskName; + if (!m_maskNodes.empty() && m_maskNodes.size() == m_imageNodes.size()) { + maskName = " / " + m_maskNodes.at(section)->GetName(); } + return QVariant((m_imageNodes.at(section)->GetName() + maskName + std::to_string(section)).c_str()); } } - else if ((Qt::DisplayRole == role) && (Qt::Vertical == orientation)){ - if (!m_statistics.empty()) { - if (m_viewMode == viewMode::imageXStatistic) { - return QVariant(m_statisticNames.at(section).c_str()); - } + } + else if ((Qt::DisplayRole == role) && (Qt::Vertical == orientation)) { + if (m_statistics.IsNotNull()) { + if (m_viewMode == viewMode::imageXStatistic) { + return QVariant(m_statisticNames.at(section).c_str()); } } - return QVariant(); + } + return QVariant(); } -void QmitkImageStatisticsTableModel::SetStatistics(const std::vector& statistics) +void QmitkImageStatisticsTableModel::SetStatistics(mitk::StatisticsContainer::ConstPointer statistics) { emit beginResetModel(); m_statistics = statistics; - if (m_statisticNames.empty() && !statistics.empty()) { - auto firstStatisticAsMap = statistics.front()->GetStatisticsAsOrderedVector(); + if (m_statisticNames.empty() && m_statistics.IsNotNull() && m_statistics->TimeStepExists(0)) { + auto firstStatisticAsMap = m_statistics->GetStatisticsAsOrderedVector(0); for (const auto& keyValue : firstStatisticAsMap) { m_statisticNames.push_back(keyValue.first); } } emit endResetModel(); emit dataChanged(QModelIndex(), QModelIndex()); } void QmitkImageStatisticsTableModel::SetImageNodes(const std::vector& nodes) { std::vector tempNodes; for (int i = 0; i < nodes.size(); i++) { int timeSteps = nodes.at(i)->GetData()->GetTimeSteps(); for (int j = 0; j < timeSteps; j++) { tempNodes.push_back(nodes.at(i)); } } m_imageNodes = tempNodes; } void QmitkImageStatisticsTableModel::SetMaskNodes(const std::vector& nodes) { m_maskNodes = nodes; } void QmitkImageStatisticsTableModel::SetViewMode(viewMode m) { m_viewMode = m; } void QmitkImageStatisticsTableModel::SetStatisticsToShow(const std::vector& statisticNames) { m_statisticNamesToShow = statisticNames; } void QmitkImageStatisticsTableModel::SetStatisticsToIgnore(const std::vector& statisticNames) { m_statisticNamesToIgnore = statisticNames; } void QmitkImageStatisticsTableModel::Clear() { emit beginResetModel(); - m_statistics.clear(); + m_statistics = nullptr; m_imageNodes.clear(); m_maskNodes.clear(); m_statisticNamesToIgnore.clear(); m_statisticNamesToShow.clear(); m_statisticNames.clear(); emit endResetModel(); emit dataChanged(QModelIndex(), QModelIndex()); } diff --git a/Modules/ImageStatisticsUI/Qmitk/QmitkImageStatisticsTableModel.h b/Modules/ImageStatisticsUI/Qmitk/QmitkImageStatisticsTableModel.h index bd46a3e727..930f0308e1 100644 --- a/Modules/ImageStatisticsUI/Qmitk/QmitkImageStatisticsTableModel.h +++ b/Modules/ImageStatisticsUI/Qmitk/QmitkImageStatisticsTableModel.h @@ -1,69 +1,69 @@ /*=================================================================== 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 QmitkImageStatisticsTableModel_h #define QmitkImageStatisticsTableModel_h #include //MITK #include #include #include /*! \class QmitkImageStatisticsTableModel Model that takes a mitk::StatisticsContainer and represents it as model in context of the QT view-model-concept. */ class MITKIMAGESTATISTICSUI_EXPORT QmitkImageStatisticsTableModel : public QAbstractTableModel { Q_OBJECT public: enum class viewMode { imageXStatistic, imageXMask }; QmitkImageStatisticsTableModel(QObject *parent = nullptr); virtual ~QmitkImageStatisticsTableModel() {}; - void SetStatistics(const std::vector& statistics); + void SetStatistics(mitk::StatisticsContainer::ConstPointer statistics); void SetImageNodes(const std::vector& nodes); void SetMaskNodes(const std::vector& nodes); void SetViewMode(viewMode m); void SetStatisticsToShow(const std::vector& statisticNames); void SetStatisticsToIgnore(const std::vector& statisticNames); void Clear(); virtual Qt::ItemFlags flags(const QModelIndex &index) const override; virtual QVariant data(const QModelIndex &index, int role) const override; virtual QVariant headerData(int section, Qt::Orientation orientation, int role) const override; virtual int rowCount(const QModelIndex &parent = QModelIndex()) const override; virtual int columnCount(const QModelIndex &parent = QModelIndex()) const override; private: - std::vector m_statistics; + mitk::StatisticsContainer::ConstPointer m_statistics; std::vector m_imageNodes; std::vector m_maskNodes; std::vector m_statisticNamesToShow; std::vector m_statisticNamesToIgnore; std::vector m_statisticNames; viewMode m_viewMode = viewMode::imageXStatistic; }; #endif // mitkQmitkImageStatisticsTableModel_h diff --git a/Modules/ImageStatisticsUI/Qmitk/QmitkImageStatisticsWidget.cpp b/Modules/ImageStatisticsUI/Qmitk/QmitkImageStatisticsWidget.cpp index 7f4c4c42f0..f812e83e03 100644 --- a/Modules/ImageStatisticsUI/Qmitk/QmitkImageStatisticsWidget.cpp +++ b/Modules/ImageStatisticsUI/Qmitk/QmitkImageStatisticsWidget.cpp @@ -1,90 +1,90 @@ /*=================================================================== 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 "QmitkImageStatisticsWidget.h" #include "QmitkTableModelToStringConverter.h" #include #include QmitkImageStatisticsWidget::QmitkImageStatisticsWidget(QWidget* parent) : QWidget(parent) { m_Controls.setupUi(this); m_imageStatisticsModel = new QmitkImageStatisticsTableModel(parent); m_Controls.checkBox4dCompleteTable->setVisible(false); CreateConnections(); + m_ProxyModel = new QSortFilterProxyModel(this); + m_Controls.tableViewStatistics->setModel(m_ProxyModel); + m_ProxyModel->setSourceModel(m_imageStatisticsModel); } -void QmitkImageStatisticsWidget::SetStatistics(const std::vector& sc) +void QmitkImageStatisticsWidget::SetStatistics(mitk::StatisticsContainer::ConstPointer sc) { m_imageStatisticsModel->SetStatistics(sc); - QSortFilterProxyModel *proxyModel = new QSortFilterProxyModel(this); - proxyModel->setSourceModel(m_imageStatisticsModel); - m_Controls.tableViewStatistics->setModel(proxyModel); - m_Controls.tableViewStatistics->resizeColumnsToContents(); + m_Controls.tableViewStatistics->resizeColumnsToContents(); // should call data in table model m_Controls.tableViewStatistics->resizeRowsToContents(); EnableAllGUIElements(); } void QmitkImageStatisticsWidget::SetImageNodes(const std::vector& nodes) { m_imageStatisticsModel->SetImageNodes(nodes); } void QmitkImageStatisticsWidget::SetMaskNodes(const std::vector& nodes) { m_imageStatisticsModel->SetMaskNodes(nodes); } void QmitkImageStatisticsWidget::Reset() { DisableAllGUIElements(); m_imageStatisticsModel->Clear(); } void QmitkImageStatisticsWidget::CreateConnections() { connect(m_Controls.buttonCopyImageStatisticsToClipboard, &QPushButton::clicked, this, &QmitkImageStatisticsWidget::OnClipboardButtonClicked); } void QmitkImageStatisticsWidget::EnableAllGUIElements() { this->setEnabled(true); m_Controls.buttonCopyImageStatisticsToClipboard->setEnabled(true); //temporarily disabled because 4D clipboard functionality is not implemented yet //m_Controls.checkBox4dCompleteTable->setEnabled(true); m_Controls.tableViewStatistics->setEnabled(true); } void QmitkImageStatisticsWidget::DisableAllGUIElements() { this->setEnabled(false); m_Controls.buttonCopyImageStatisticsToClipboard->setEnabled(false); m_Controls.checkBox4dCompleteTable->setEnabled(false); m_Controls.tableViewStatistics->setEnabled(false); } void QmitkImageStatisticsWidget::OnClipboardButtonClicked() { QmitkTableModelToStringConverter converter; converter.SetTableModel(m_imageStatisticsModel); converter.SetIncludeHeaderData(true); converter.SetColumnDelimiter('\t'); QString clipboardAsString = converter.GetString(); QApplication::clipboard()->setText(clipboardAsString, QClipboard::Clipboard); } diff --git a/Modules/ImageStatisticsUI/Qmitk/QmitkImageStatisticsWidget.h b/Modules/ImageStatisticsUI/Qmitk/QmitkImageStatisticsWidget.h index eb5c3bd002..16d5957153 100644 --- a/Modules/ImageStatisticsUI/Qmitk/QmitkImageStatisticsWidget.h +++ b/Modules/ImageStatisticsUI/Qmitk/QmitkImageStatisticsWidget.h @@ -1,53 +1,56 @@ /*=================================================================== 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 QmitkImageStatisticsWidget_H__INCLUDED #define QmitkImageStatisticsWidget_H__INCLUDED #include #include #include #include +class QSortFilterProxyModel; + //Qt #include class MITKIMAGESTATISTICSUI_EXPORT QmitkImageStatisticsWidget : public QWidget { Q_OBJECT public: QmitkImageStatisticsWidget(QWidget* parent = nullptr); - void SetStatistics(const std::vector& sc); + void SetStatistics(mitk::StatisticsContainer::ConstPointer sc); void SetImageNodes(const std::vector& nodes); void SetMaskNodes(const std::vector& nodes); void Reset(); private: void CreateConnections(); void EnableAllGUIElements(); void DisableAllGUIElements(); /** \brief Saves the image statistics to the clipboard */ void OnClipboardButtonClicked(); /** \brief Toogle GUI elements if histogram default bin size checkbox value changed. */ private: Ui::QmitkImageStatisticsControls m_Controls; QmitkImageStatisticsTableModel* m_imageStatisticsModel; + QSortFilterProxyModel* m_ProxyModel; }; #endif // QmitkImageStatisticsWidget_H__INCLUDED diff --git a/Modules/Segmentation/Algorithms/mitkFeatureBasedEdgeDetectionFilter.cpp b/Modules/Segmentation/Algorithms/mitkFeatureBasedEdgeDetectionFilter.cpp index f8875e48d2..408c0288b0 100644 --- a/Modules/Segmentation/Algorithms/mitkFeatureBasedEdgeDetectionFilter.cpp +++ b/Modules/Segmentation/Algorithms/mitkFeatureBasedEdgeDetectionFilter.cpp @@ -1,198 +1,198 @@ /*=================================================================== 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 "mitkFeatureBasedEdgeDetectionFilter.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include mitk::FeatureBasedEdgeDetectionFilter::FeatureBasedEdgeDetectionFilter() { this->SetNumberOfRequiredInputs(1); this->SetNumberOfIndexedOutputs(1); } mitk::FeatureBasedEdgeDetectionFilter::~FeatureBasedEdgeDetectionFilter() { } void mitk::FeatureBasedEdgeDetectionFilter::GenerateData() { mitk::Image::Pointer image = ImageToUnstructuredGridFilter::GetInput(); if (m_SegmentationMask.IsNull()) { MITK_WARN << "Please set a segmentation mask first" << std::endl; return; } // First create a threshold segmentation of the image. The threshold is determined // by the mean +/- stddev of the pixel values that are covered by the segmentation mask // Compute mean and stdDev based on the current segmentation mitk::ImageStatisticsCalculator::Pointer statCalc = mitk::ImageStatisticsCalculator::New(); statCalc->SetInputImage(image); mitk::ImageMaskGenerator::Pointer imgMask = mitk::ImageMaskGenerator::New(); imgMask->SetImageMask(m_SegmentationMask); auto stats = statCalc->GetStatistics(); - double mean = stats->GetMean(); - double stdDev = stats->GetStd(); + double mean = stats.m_Mean; + double stdDev = stats.m_Std; double upperThreshold = mean + stdDev; double lowerThreshold = mean - stdDev; // Perform thresholding mitk::Image::Pointer thresholdImage = mitk::Image::New(); AccessByItk_3(image.GetPointer(), ITKThresholding, lowerThreshold, upperThreshold, thresholdImage) mitk::ProgressBar::GetInstance() ->Progress(2); // Postprocess threshold segmentation // First a closing will be executed mitk::Image::Pointer closedImage = mitk::Image::New(); AccessByItk_1(thresholdImage, ThreadedClosing, closedImage); // Then we will holes that might exist mitk::MorphologicalOperations::FillHoles(closedImage); mitk::ProgressBar::GetInstance()->Progress(); // Extract the binary edges of the resulting segmentation mitk::Image::Pointer edgeImage = mitk::Image::New(); AccessByItk_1(closedImage, ContourSearch, edgeImage); // Convert the edge image into an unstructured grid mitk::ImageToUnstructuredGridFilter::Pointer i2UFilter = mitk::ImageToUnstructuredGridFilter::New(); i2UFilter->SetInput(edgeImage); i2UFilter->SetThreshold(1.0); i2UFilter->Update(); m_PointGrid = this->GetOutput(); if (m_PointGrid.IsNull()) m_PointGrid = mitk::UnstructuredGrid::New(); m_PointGrid->SetVtkUnstructuredGrid(i2UFilter->GetOutput()->GetVtkUnstructuredGrid()); mitk::ProgressBar::GetInstance()->Progress(); } template void mitk::FeatureBasedEdgeDetectionFilter::ThreadedClosing(itk::Image *originalImage, mitk::Image::Pointer &result) { typedef itk::BinaryBallStructuringElement myKernelType; myKernelType ball; ball.SetRadius(1); ball.CreateStructuringElement(); typedef typename itk::Image ImageType; typename itk::DilateObjectMorphologyImageFilter::Pointer dilationFilter = itk::DilateObjectMorphologyImageFilter::New(); dilationFilter->SetInput(originalImage); dilationFilter->SetKernel(ball); dilationFilter->Update(); typename itk::Image::Pointer dilatedImage = dilationFilter->GetOutput(); typename itk::ErodeObjectMorphologyImageFilter::Pointer erodeFilter = itk::ErodeObjectMorphologyImageFilter::New(); erodeFilter->SetInput(dilatedImage); erodeFilter->SetKernel(ball); erodeFilter->Update(); mitk::GrabItkImageMemory(erodeFilter->GetOutput(), result); } template void mitk::FeatureBasedEdgeDetectionFilter::ContourSearch(itk::Image *originalImage, mitk::Image::Pointer &result) { typedef itk::Image ImageType; typedef itk::BinaryContourImageFilter binaryContourImageFilterType; typename binaryContourImageFilterType::Pointer binaryContourFilter = binaryContourImageFilterType::New(); binaryContourFilter->SetInput(originalImage); binaryContourFilter->SetForegroundValue(1); binaryContourFilter->SetBackgroundValue(0); binaryContourFilter->Update(); typename itk::Image::Pointer itkImage = itk::Image::New(); itkImage->Graft(binaryContourFilter->GetOutput()); mitk::GrabItkImageMemory(itkImage, result); } template void mitk::FeatureBasedEdgeDetectionFilter::ITKThresholding(itk::Image *originalImage, double lower, double upper, mitk::Image::Pointer &result) { typedef itk::Image ImageType; typedef itk::Image SegmentationType; typedef itk::BinaryThresholdImageFilter ThresholdFilterType; if (typeid(TPixel) != typeid(float) && typeid(TPixel) != typeid(double)) { // round the thresholds if we have nor a float or double image lower = std::floor(lower + 0.5); upper = std::floor(upper - 0.5); } if (lower >= upper) { upper = lower; } typename ThresholdFilterType::Pointer filter = ThresholdFilterType::New(); filter->SetInput(originalImage); filter->SetLowerThreshold(lower); filter->SetUpperThreshold(upper); filter->SetInsideValue(1); filter->SetOutsideValue(0); filter->Update(); mitk::GrabItkImageMemory(filter->GetOutput(), result); } void mitk::FeatureBasedEdgeDetectionFilter::SetSegmentationMask(mitk::Image::Pointer segmentation) { this->m_SegmentationMask = segmentation; } void mitk::FeatureBasedEdgeDetectionFilter::GenerateOutputInformation() { Superclass::GenerateOutputInformation(); } diff --git a/Modules/SurfaceInterpolation/mitkImageToPointCloudFilter.cpp b/Modules/SurfaceInterpolation/mitkImageToPointCloudFilter.cpp index 5127e87c30..18fad7f8de 100644 --- a/Modules/SurfaceInterpolation/mitkImageToPointCloudFilter.cpp +++ b/Modules/SurfaceInterpolation/mitkImageToPointCloudFilter.cpp @@ -1,166 +1,166 @@ /*=================================================================== 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 "mitkImageToPointCloudFilter.h" #include #include #include #include #include #include #include #include #include mitk::ImageToPointCloudFilter::ImageToPointCloudFilter() { m_Method = DetectionMethod(0); this->SetNumberOfRequiredInputs(1); this->SetNumberOfIndexedOutputs(1); } mitk::ImageToPointCloudFilter::~ImageToPointCloudFilter() { } void mitk::ImageToPointCloudFilter::GenerateData() { mitk::Image::ConstPointer image = ImageToUnstructuredGridFilter::GetInput(); m_Geometry = image->GetGeometry(); if (image.IsNull()) { MITK_ERROR << "mitk::ImageToContourFilter: No input available. " "Please set the input!" << std::endl; return; } mitk::Image::Pointer notConstImage = const_cast(image.GetPointer()); switch (m_Method) { case 0: AccessByItk_1(notConstImage.GetPointer(), StdDeviations, 2) break; case 1: AccessByItk_1(notConstImage.GetPointer(), StdDeviations, 3) break; case 2: AccessByItk_1(notConstImage.GetPointer(), StdDeviations, 4) break; default: AccessByItk_1(notConstImage.GetPointer(), StdDeviations, 2) break; } } template void mitk::ImageToPointCloudFilter::StdDeviations(itk::Image *image, int amount) { typedef itk::Image InputImageType; typedef itk::CastImageFilter ImagePTypeToFloatPTypeCasterType; typedef itk::LaplacianImageFilter LaplacianFilterType; typename LaplacianFilterType::Pointer lapFilter = LaplacianFilterType::New(); typename ImagePTypeToFloatPTypeCasterType::Pointer caster = ImagePTypeToFloatPTypeCasterType::New(); caster->SetInput(image); caster->Update(); FloatImageType::Pointer fImage = caster->GetOutput(); lapFilter->SetInput(fImage); lapFilter->UpdateLargestPossibleRegion(); mitk::Image::Pointer edgeImage = mitk::ImportItkImage(lapFilter->GetOutput()); mitk::ImageStatisticsCalculator::Pointer statCalc = mitk::ImageStatisticsCalculator::New(); statCalc->SetInputImage(edgeImage); auto stats = statCalc->GetStatistics(); - double mean = stats->GetMean(); - double stdDev = stats->GetStd(); + double mean = stats.m_Mean; + double stdDev = stats.m_Std; double upperThreshold = mean + stdDev * amount; double lowerThreshold = mean - stdDev * amount; typename itk::ImageRegionIterator it(lapFilter->GetOutput(), lapFilter->GetOutput()->GetRequestedRegion()); vtkSmartPointer points = vtkSmartPointer::New(); double greatX = 0, greatY = 0, greatZ = 0; it.GoToBegin(); while (!it.IsAtEnd()) { if (it.Get() > lowerThreshold && it.Get() < upperThreshold) { it.Set(0); } else { it.Set(1); mitk::Point3D imagePoint; mitk::Point3D worldPoint; imagePoint[0] = it.GetIndex()[0]; imagePoint[1] = it.GetIndex()[1]; imagePoint[2] = it.GetIndex()[2]; m_Geometry->IndexToWorld(imagePoint, worldPoint); if (worldPoint[0] > greatX) greatX = worldPoint[0]; if (worldPoint[1] > greatY) greatY = worldPoint[1]; if (worldPoint[2] > greatZ) greatZ = worldPoint[2]; points->InsertNextPoint(worldPoint[0], worldPoint[1], worldPoint[2]); m_NumberOfExtractedPoints++; } ++it; } /*need to build the UnstructuredGrid with at least one vertex otherwise its not visible*/ vtkSmartPointer verts = vtkSmartPointer::New(); verts->GetPointIds()->SetNumberOfIds(m_NumberOfExtractedPoints); for (int i = 0; i < m_NumberOfExtractedPoints; i++) { verts->GetPointIds()->SetId(i, i); } vtkSmartPointer uGrid = vtkSmartPointer::New(); uGrid->Allocate(1); uGrid->InsertNextCell(verts->GetCellType(), verts->GetPointIds()); uGrid->SetPoints(points); mitk::UnstructuredGrid::Pointer outputGrid = mitk::UnstructuredGrid::New(); outputGrid->SetVtkUnstructuredGrid(uGrid); this->SetNthOutput(0, outputGrid); } void mitk::ImageToPointCloudFilter::GenerateOutputInformation() { Superclass::GenerateOutputInformation(); } diff --git a/Plugins/org.mitk.gui.qt.measurementtoolbox/src/internal/QmitkImageStatisticsReloadedView.cpp b/Plugins/org.mitk.gui.qt.measurementtoolbox/src/internal/QmitkImageStatisticsReloadedView.cpp index 3de6ca7bca..a533ed4009 100644 --- a/Plugins/org.mitk.gui.qt.measurementtoolbox/src/internal/QmitkImageStatisticsReloadedView.cpp +++ b/Plugins/org.mitk.gui.qt.measurementtoolbox/src/internal/QmitkImageStatisticsReloadedView.cpp @@ -1,398 +1,405 @@ /*=================================================================== 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 "QmitkImageStatisticsReloadedView.h" #include // berry includes #include #include #include #include #include #include #include #include #include #include #include #include #include const std::string QmitkImageStatisticsReloadedView::VIEW_ID = "org.mitk.views.imagestatisticsReloaded"; QmitkImageStatisticsReloadedView::QmitkImageStatisticsReloadedView(QObject* /*parent*/, const char* /*name*/) { this->m_CalculationThread = new QmitkImageStatisticsCalculationJob(); } QmitkImageStatisticsReloadedView::~QmitkImageStatisticsReloadedView() { if (m_selectedPlanarFigure) m_selectedPlanarFigure->RemoveObserver(m_PlanarFigureObserverTag); } void QmitkImageStatisticsReloadedView::CreateQtPartControl(QWidget *parent) { m_Controls.setupUi(parent); m_Controls.widget_histogram->SetTheme(this->GetColorTheme()); m_Controls.widget_intensityProfile->SetTheme(this->GetColorTheme()); m_Controls.groupBox_histogram->setVisible(true); m_Controls.groupBox_intensityProfile->setVisible(false); m_Controls.label_currentlyComputingStatistics->setVisible(false); m_Controls.sliderWidget_histogram->setPrefix("Time: "); m_Controls.sliderWidget_histogram->setDecimals(0); m_Controls.sliderWidget_histogram->setVisible(false); m_statisticsManager = mitk::ImageStatisticsContainerManager::New(); m_statisticsManager->SetDataStorage(this->GetDataStorage().GetPointer()); PrepareDataStorageComboBoxes(); CreateConnections(); } void QmitkImageStatisticsReloadedView::CreateConnections() { connect(this->m_CalculationThread, &QmitkImageStatisticsCalculationJob::finished, this, &QmitkImageStatisticsReloadedView::OnStatisticsCalculationEnds, Qt::QueuedConnection); connect(this->m_Controls.checkBox_ignoreZero, &QCheckBox::stateChanged, this, &QmitkImageStatisticsReloadedView::OnCheckBoxIgnoreZeroStateChanged); connect(this->m_Controls.sliderWidget_histogram, &ctkSliderWidget::valueChanged, this, &QmitkImageStatisticsReloadedView::OnSliderWidgetHistogramChanged); } void QmitkImageStatisticsReloadedView::OnCheckBoxIgnoreZeroStateChanged(int state) { m_ForceRecompute = true; if (state != Qt::Unchecked) { this->m_CalculationThread->SetIgnoreZeroValueVoxel(true); } else { this->m_CalculationThread->SetIgnoreZeroValueVoxel(false); } OnImageOrMaskSelectorChanged(); } void QmitkImageStatisticsReloadedView::OnSliderWidgetHistogramChanged(double value) { unsigned int timeStep = static_cast(value); - HistogramType::ConstPointer histogram = this->m_CalculationThread->GetTimeStepHistogram(timeStep); + auto mask = m_selectedMaskNode ? m_selectedMaskNode->GetData() : nullptr; + auto imageStatistics = m_statisticsManager->GetImageStatistics(m_selectedImageNode->GetData(), mask); + HistogramType::ConstPointer histogram = imageStatistics->GetStatisticsForTimeStep(timeStep).m_Histogram; if (histogram.IsNotNull() && this->m_CalculationThread->GetStatisticsUpdateSuccessFlag()) { this->FillHistogramWidget({ histogram }, { m_selectedImageNode->GetName() }); } } -void QmitkImageStatisticsReloadedView::PartClosed(const berry::IWorkbenchPartReference::Pointer& ) +void QmitkImageStatisticsReloadedView::PartClosed(const berry::IWorkbenchPartReference::Pointer&) { } -void QmitkImageStatisticsReloadedView::FillStatisticsWidget(const std::vector& statistics) +void QmitkImageStatisticsReloadedView::FillStatisticsWidget(mitk::StatisticsContainer::ConstPointer statistics) { m_Controls.widget_statistics->Reset(); m_Controls.widget_statistics->SetStatistics(statistics); m_Controls.widget_statistics->SetImageNodes({ m_selectedImageNode }); if (m_selectedMaskNode) { m_Controls.widget_statistics->SetMaskNodes({ m_selectedMaskNode }); } m_Controls.widget_statistics->setEnabled(true); } void QmitkImageStatisticsReloadedView::FillHistogramWidget(const std::vector& histogram, const std::vector& dataLabels) { m_Controls.groupBox_histogram->setVisible(true); m_Controls.widget_histogram->SetTheme(this->GetColorTheme()); m_Controls.widget_histogram->Reset(); m_Controls.widget_histogram->SetHistogram(histogram.front(), dataLabels.front()); connect(m_Controls.widget_histogram, &QmitkHistogramVisualizationWidget::RequestHistogramUpdate, this, &QmitkImageStatisticsReloadedView::OnRequestHistogramUpdate); } QmitkChartWidget::ChartStyle QmitkImageStatisticsReloadedView::GetColorTheme() const { ctkPluginContext* context = berry::WorkbenchPlugin::GetDefault()->GetPluginContext(); ctkServiceReference styleManagerRef = context->getServiceReference(); if (styleManagerRef) { auto styleManager = context->getService(styleManagerRef); if (styleManager->GetStyle().name == "Dark") { return QmitkChartWidget::ChartStyle::darkstyle; } else { return QmitkChartWidget::ChartStyle::lightstyle; } } return QmitkChartWidget::ChartStyle::darkstyle; } void QmitkImageStatisticsReloadedView::OnImageOrMaskSelectorChanged() { if (this->m_selectedPlanarFigure) { this->m_selectedPlanarFigure->RemoveObserver(this->m_PlanarFigureObserverTag); this->m_selectedPlanarFigure = nullptr; } m_selectedImageNode = m_Controls.imageSelector->GetSelectedNode(); m_selectedMaskNode = m_Controls.maskImageSelector->GetSelectedNode(); m_Controls.groupBox_intensityProfile->setVisible(false); //m_statisticContainerRules.clear(); if (m_selectedImageNode != nullptr) { auto image = dynamic_cast(m_selectedImageNode->GetData()); mitk::Image::Pointer mask = nullptr; mitk::PlanarFigure::Pointer maskPlanarFigure = nullptr; - if (image->GetDimension() == 4){ + if (image->GetDimension() == 4) { m_Controls.sliderWidget_histogram->setVisible(true); unsigned int maxTimestep = image->GetTimeSteps(); m_Controls.sliderWidget_histogram->setMaximum(maxTimestep); } else { m_Controls.sliderWidget_histogram->setVisible(false); } if (m_selectedMaskNode != nullptr) { mask = dynamic_cast(m_selectedMaskNode->GetData()); if (mask == nullptr) { maskPlanarFigure = dynamic_cast(m_selectedMaskNode->GetData()); } } if (mask) { auto imageStatistics = m_statisticsManager->GetImageStatistics(image, mask.GetPointer()); bool imageStatisticsOlderThanInputs = false; if (imageStatistics && (imageStatistics->GetMTime() < image->GetMTime() || imageStatistics->GetMTime() < mask->GetMTime())) { imageStatisticsOlderThanInputs = true; } - //compute statistics with given mask if (!imageStatistics || imageStatisticsOlderThanInputs || m_ForceRecompute) { CalculateStatistics(image, mask.GetPointer()); } //statistics already computed else { this->FillStatisticsWidget({ imageStatistics }); - this->FillHistogramWidget({ imageStatistics->GetHistogram().GetPointer() }, { m_selectedImageNode->GetName() }); + if (imageStatistics->TimeStepExists(0)) { + this->FillHistogramWidget({ imageStatistics->GetStatisticsForTimeStep(0).m_Histogram.GetPointer() }, { m_selectedImageNode->GetName() }); + } } } else if (maskPlanarFigure) { m_selectedPlanarFigure = maskPlanarFigure; ITKCommandType::Pointer changeListener = ITKCommandType::New(); changeListener->SetCallbackFunction(this, &QmitkImageStatisticsReloadedView::OnImageOrMaskSelectorChanged); this->m_PlanarFigureObserverTag = m_selectedPlanarFigure->AddObserver(mitk::EndInteractionPlanarFigureEvent(), changeListener); if (!maskPlanarFigure->IsClosed()) { //compute line profile and display statistics for voxels on line auto intensityProfile = mitk::ComputeIntensityProfile(image, maskPlanarFigure); //Don't show histogram for intensity profiles m_Controls.groupBox_histogram->setVisible(false); m_Controls.groupBox_intensityProfile->setVisible(true); m_Controls.widget_intensityProfile->Reset(); m_Controls.widget_intensityProfile->SetIntensityProfile(intensityProfile.GetPointer(), "Intensity Profile of " + m_selectedImageNode->GetName()); } auto imageStatistics = m_statisticsManager->GetImageStatistics(image, maskPlanarFigure.GetPointer()); bool imageStatisticsOlderThanInputs = false; if (imageStatistics && (imageStatistics->GetMTime() < image->GetMTime() || imageStatistics->GetMTime() < maskPlanarFigure->GetMTime())) { imageStatisticsOlderThanInputs = true; } //for all planar figures: compute statistics with planarFigure as mask if (!imageStatistics || imageStatisticsOlderThanInputs || m_ForceRecompute) { CalculateStatistics(image, nullptr, maskPlanarFigure.GetPointer()); } //statistics already computed else { this->FillStatisticsWidget({ imageStatistics }); if (maskPlanarFigure->IsClosed()) { - this->FillHistogramWidget({ imageStatistics->GetHistogram().GetPointer() }, { m_selectedImageNode->GetName() }); - } + if (imageStatistics->TimeStepExists(0)) { + this->FillHistogramWidget({ imageStatistics->GetStatisticsForTimeStep(0).m_Histogram.GetPointer() }, { m_selectedImageNode->GetName() }); + } } + } } else { auto imageStatistics = m_statisticsManager->GetImageStatistics(image); bool imageStatisticsOlderThanInputs = false; if (imageStatistics && (imageStatistics->GetMTime() < image->GetMTime())) { imageStatisticsOlderThanInputs = true; } //compute statistics with image only if (!imageStatistics || imageStatisticsOlderThanInputs || m_ForceRecompute) { CalculateStatistics(image); } //statistics already computed else { this->FillStatisticsWidget({ imageStatistics }); - this->FillHistogramWidget({ imageStatistics->GetHistogram().GetPointer() }, { m_selectedImageNode->GetName() }); + if (imageStatistics->TimeStepExists(0)) { + this->FillHistogramWidget({ imageStatistics->GetStatisticsForTimeStep(0).m_Histogram.GetPointer() }, { m_selectedImageNode->GetName() }); + } } } } else { ResetGUI(); } m_ForceRecompute = false; } void QmitkImageStatisticsReloadedView::ResetGUI() { m_Controls.widget_statistics->Reset(); m_Controls.widget_statistics->setEnabled(false); m_Controls.widget_histogram->Reset(); m_Controls.widget_histogram->setEnabled(false); } void QmitkImageStatisticsReloadedView::OnStatisticsCalculationEnds() { mitk::StatusBar::GetInstance()->Clear(); if (this->m_CalculationThread->GetStatisticsUpdateSuccessFlag()) { - auto statistics = m_CalculationThread->GetStatisticsData(); - for (auto& statistic : statistics) { + auto statistic = m_CalculationThread->GetStatisticsData(); + //for (auto& statistic : statistics) { mitk::PropertyRelations::RuleResultVectorType rulesForCurrentStatistic; auto statisticNonConst = statistic->Clone(); auto statisticsNodeName = m_selectedImageNode->GetName(); if (m_selectedMaskNode) { statisticsNodeName += "_" + m_selectedMaskNode->GetName(); } statisticsNodeName += "_statistics"; auto statisticsNode = mitk::CreateImageStatisticsNode(statisticNonConst, statisticsNodeName); auto imageRule = mitk::StatisticsToImageRelationRule::New(); imageRule->Connect(statisticNonConst.GetPointer(), m_CalculationThread->GetStatisticsImage().GetPointer()); rulesForCurrentStatistic.push_back(imageRule.GetPointer()); if (m_CalculationThread->GetMaskImage()) { auto maskRule = mitk::StatisticsToMaskRelationRule::New(); maskRule->Connect(statisticNonConst.GetPointer(), m_CalculationThread->GetMaskImage().GetPointer()); rulesForCurrentStatistic.push_back(maskRule.GetPointer()); } else if (m_CalculationThread->GetPlanarFigure()) { auto planarFigureRule = mitk::StatisticsToMaskRelationRule::New(); planarFigureRule->Connect(statisticNonConst.GetPointer(), m_CalculationThread->GetPlanarFigure().GetPointer()); rulesForCurrentStatistic.push_back(planarFigureRule.GetPointer()); } m_statisticContainerRules.push_back(rulesForCurrentStatistic); this->GetDataStorage()->Add(statisticsNode); m_statisticsManager->SetRules(m_statisticContainerRules); - } + //} - this->FillStatisticsWidget(statistics); + this->FillStatisticsWidget(statistic); if (!m_selectedPlanarFigure || m_selectedPlanarFigure->IsClosed()) { this->FillHistogramWidget({ m_CalculationThread->GetTimeStepHistogram() }, { m_selectedImageNode->GetName() }); } } else { mitk::StatusBar::GetInstance()->DisplayErrorText(m_CalculationThread->GetLastErrorMessage().c_str()); m_Controls.widget_histogram->setEnabled(false); m_Controls.widget_statistics->setEnabled(false); } m_Controls.label_currentlyComputingStatistics->setVisible(false); } void QmitkImageStatisticsReloadedView::OnRequestHistogramUpdate(unsigned int nBins) { m_CalculationThread->SetHistogramNBins(nBins); m_CalculationThread->start(); } void QmitkImageStatisticsReloadedView::CalculateStatistics(mitk::Image::ConstPointer image, mitk::Image::ConstPointer mask, mitk::PlanarFigure::ConstPointer maskPlanarFigure) { this->m_StatisticsUpdatePending = true; auto renderPart = this->GetRenderWindowPart(); unsigned int timeStep = renderPart->GetTimeNavigationController()->GetTime()->GetPos(); this->m_CalculationThread->Initialize(image, mask, maskPlanarFigure); this->m_CalculationThread->SetTimeStep(timeStep); try { // Compute statistics this->m_CalculationThread->start(); m_Controls.label_currentlyComputingStatistics->setVisible(true); } catch (const mitk::Exception& e) { mitk::StatusBar::GetInstance()->DisplayErrorText(e.GetDescription()); this->m_StatisticsUpdatePending = false; m_Controls.label_currentlyComputingStatistics->setVisible(false); } catch (const std::runtime_error &e) { mitk::StatusBar::GetInstance()->DisplayErrorText(e.what()); this->m_StatisticsUpdatePending = false; m_Controls.label_currentlyComputingStatistics->setVisible(false); } catch (const std::exception &e) { mitk::StatusBar::GetInstance()->DisplayErrorText(e.what()); this->m_StatisticsUpdatePending = false; m_Controls.label_currentlyComputingStatistics->setVisible(false); } } -void QmitkImageStatisticsReloadedView::OnSelectionChanged( berry::IWorkbenchPart::Pointer part, const QList &nodes ) +void QmitkImageStatisticsReloadedView::OnSelectionChanged(berry::IWorkbenchPart::Pointer part, const QList &nodes) { - Q_UNUSED(part); - Q_UNUSED(nodes); + Q_UNUSED(part); + Q_UNUSED(nodes); } void QmitkImageStatisticsReloadedView::PrepareDataStorageComboBoxes() { auto isPlanarFigurePredicate = mitk::GetImageStatisticsPlanarFigurePredicate(); auto isMaskPredicate = mitk::GetImageStatisticsMaskPredicate(); auto isImagePredicate = mitk::GetImageStatisticsImagePredicate(); auto isMaskOrPlanarFigurePredicate = mitk::NodePredicateOr::New(isPlanarFigurePredicate, isMaskPredicate); m_Controls.imageSelector->SetDataStorage(GetDataStorage()); m_Controls.imageSelector->SetPredicate(isImagePredicate); m_Controls.maskImageSelector->SetDataStorage(GetDataStorage()); m_Controls.maskImageSelector->SetPredicate(isMaskOrPlanarFigurePredicate); m_Controls.maskImageSelector->SetZeroEntryText(""); } void QmitkImageStatisticsReloadedView::Activated() { } void QmitkImageStatisticsReloadedView::Deactivated() { } void QmitkImageStatisticsReloadedView::Visible() { connect(m_Controls.imageSelector, static_cast(&QComboBox::currentIndexChanged), this, &QmitkImageStatisticsReloadedView::OnImageOrMaskSelectorChanged); connect(m_Controls.maskImageSelector, static_cast(&QComboBox::currentIndexChanged), this, &QmitkImageStatisticsReloadedView::OnImageOrMaskSelectorChanged); m_selectedImageNode = m_Controls.imageSelector->GetSelectedNode(); if (m_selectedImageNode) { OnImageOrMaskSelectorChanged(); } else { ResetGUI(); } } void QmitkImageStatisticsReloadedView::Hidden() { m_Controls.imageSelector->disconnect(); m_Controls.maskImageSelector->disconnect(); } void QmitkImageStatisticsReloadedView::SetFocus() { } diff --git a/Plugins/org.mitk.gui.qt.measurementtoolbox/src/internal/QmitkImageStatisticsReloadedView.h b/Plugins/org.mitk.gui.qt.measurementtoolbox/src/internal/QmitkImageStatisticsReloadedView.h index 94d90d1d7c..8c9b7c51ff 100644 --- a/Plugins/org.mitk.gui.qt.measurementtoolbox/src/internal/QmitkImageStatisticsReloadedView.h +++ b/Plugins/org.mitk.gui.qt.measurementtoolbox/src/internal/QmitkImageStatisticsReloadedView.h @@ -1,110 +1,110 @@ /*=================================================================== 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 QmitkImageStatisticsReloadedView_H__INCLUDED #define QmitkImageStatisticsReloadedView_H__INCLUDED #include "ui_QmitkImageStatisticsReloadedViewControls.h" // Qmitk includes #include #include #include #include #include #include #include /*! \brief QmitkImageStatisticsView is a bundle that allows statistics calculation from images. Three modes are supported: 1. Statistics of one image, 2. Statistics of an image and a segmentation, 3. Statistics of an image and a Planar Figure. The statistics calculation is realized in a separate thread to keep the gui accessible during calculation. \ingroup Plugins/org.mitk.gui.qt.measurementtoolbox */ class QmitkImageStatisticsReloadedView : public QmitkAbstractView, public mitk::ILifecycleAwarePart, public berry::IPartListener { Q_OBJECT public: using HistogramType = mitk::StatisticsContainer::HistogramType; /*! \brief default constructor */ QmitkImageStatisticsReloadedView(QObject *parent = nullptr, const char *name = nullptr); /*! \brief default destructor */ virtual ~QmitkImageStatisticsReloadedView(); /*! \brief method for creating the widget containing the application controls, like sliders, buttons etc. */ virtual void CreateQtPartControl(QWidget *parent) override; /*! \brief method for creating the connections of main and control widget */ virtual void CreateConnections(); /*! \brief Is called from the selection mechanism once the data manager selection has changed*/ void OnSelectionChanged(berry::IWorkbenchPart::Pointer part, const QList &selectedNodes) override; void PrepareDataStorageComboBoxes(); static const std::string VIEW_ID; - void FillStatisticsWidget(const std::vector& statistics); + void FillStatisticsWidget(mitk::StatisticsContainer::ConstPointer statistics); void FillHistogramWidget(const std::vector& histogram, const std::vector& dataLabels); QmitkChartWidget::ChartStyle GetColorTheme() const; protected: virtual void Activated() override; virtual void Deactivated() override; virtual void Visible() override; virtual void Hidden() override; virtual void SetFocus() override; /** \brief Is called right before the view closes (before the destructor) */ virtual void PartClosed(const berry::IWorkbenchPartReference::Pointer&) override; /** \brief Required for berry::IPartListener */ virtual Events::Types GetPartEventTypes() const override { return Events::CLOSED; } void OnImageOrMaskSelectorChanged(); void ResetGUI(); void OnStatisticsCalculationEnds(); void OnRequestHistogramUpdate(unsigned int nBins); void OnCheckBoxIgnoreZeroStateChanged(int state); void OnSliderWidgetHistogramChanged(double value); void CalculateStatistics(mitk::Image::ConstPointer image, mitk::Image::ConstPointer mask=nullptr, mitk::PlanarFigure::ConstPointer maskPlanarFigure = nullptr); - + // member variables Ui::QmitkImageStatisticsReloadedViewControls m_Controls; private: typedef itk::SimpleMemberCommand< QmitkImageStatisticsReloadedView > ITKCommandType; QmitkImageStatisticsCalculationJob * m_CalculationThread = nullptr; bool m_StatisticsUpdatePending=false; mitk::DataNode::ConstPointer m_selectedImageNode = nullptr, m_selectedMaskNode = nullptr; std::vector m_statisticContainerRules; mitk::PlanarFigure::Pointer m_selectedPlanarFigure=nullptr; mitk::ImageStatisticsContainerManager::Pointer m_statisticsManager; long m_PlanarFigureObserverTag; bool m_ForceRecompute = false; }; #endif // QmitkImageStatisticsView_H__INCLUDED