diff --git a/Modules/DiffusionImaging/MiniApps/ImageStatisticsMiniapp_v2.cpp b/Modules/DiffusionImaging/MiniApps/ImageStatisticsMiniapp_v2.cpp index 297b92f63c..fb86a2a91b 100644 --- a/Modules/DiffusionImaging/MiniApps/ImageStatisticsMiniapp_v2.cpp +++ b/Modules/DiffusionImaging/MiniApps/ImageStatisticsMiniapp_v2.cpp @@ -1,229 +1,361 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkCommandLineParser.h" #include "mitkImage.h" #include #include #include #include #include #include "mitkIOUtil.h" #include #include #include #include #include "mitkImageAccessByItk.h" #include #include #include #include #include #include #include #include #include #include #include #include struct statistics_res{ double mean, variance, min, max, count, moment; }; long getTimeInMs() { std::chrono::milliseconds ms = std::chrono::duration_cast< std::chrono::milliseconds > (std::chrono::system_clock::now().time_since_epoch()); long time = ms.count(); return time; } -void printstats(statistics_res s) +std::string printstats(mitk::ImageStatisticsCalculator::Statistics statisticsStruct) { - std::cout << "mean: " << s.mean << std::endl - << "variance: " << s.variance << std::endl - << "min: " << s.min << std::endl - << "max: " << s.max << std::endl - << "count: " << s.count << std::endl - << "moment: " << s.moment << std::endl; + std::string res = ""; + res += std::string("Entropy ") + std::to_string(statisticsStruct.GetEntropy()) + std::string("\n"); + res += std::string("Kurtosis ") + std::to_string(statisticsStruct.GetKurtosis()) + std::string("\n"); + res += std::string("MPP ") + std::to_string(statisticsStruct.GetMPP()) + std::string("\n"); + res += std::string("Max ") + std::to_string(statisticsStruct.GetMax()) + std::string("\n"); + res += std::string("Mean ") + std::to_string(statisticsStruct.GetMean()) + std::string("\n"); + res += std::string("Median ") + std::to_string(statisticsStruct.GetMedian()) + std::string("\n"); + res += std::string("Min ") + std::to_string(statisticsStruct.GetMin()) + std::string("\n"); + res += std::string("N ") + std::to_string(statisticsStruct.GetN()) + std::string("\n"); + res += std::string("RMS ") + std::to_string(statisticsStruct.GetRMS()) + std::string("\n"); + res += std::string("Skewness ") + std::to_string(statisticsStruct.GetSkewness()) + std::string("\n"); + res += std::string("Std ") + std::to_string(statisticsStruct.GetSigma()) + std::string("\n"); + res += std::string("UPP ") + std::to_string(statisticsStruct.GetUPP()) + std::string("\n"); + res += std::string("Uniformity ") + std::to_string(statisticsStruct.GetUniformity()) + std::string("\n"); + res += std::string("Variance ") + std::to_string(statisticsStruct.GetVariance()) + std::string("\n"); + vnl_vector minIndex = statisticsStruct.GetMinIndex(); + vnl_vector maxIndex = statisticsStruct.GetMaxIndex(); + + res += "Min Index: "; + for (auto it = minIndex.begin(); it != minIndex.end(); it++) + { + res += std::to_string(*it) + " "; + } + res += std::string("\n"); + + res += "Max Index: "; + for (auto it = maxIndex.begin(); it != maxIndex.end(); it++) + { + res += std::to_string(*it) + " "; + } + res += std::string("\n"); + return res; } template void printMap(std::map input) { for (auto it = input.begin(); it != input.end(); ++it) { std::cout << it->first<< ": " << it->second<< std::endl; } std::cout << std::endl; } template < typename TPixel, unsigned int VImageDimension > void get_statistics_boost(itk::Image* itkImage, statistics_res& res){ typedef itk::Image ImageType; itk::ImageRegionConstIterator it(itkImage, itkImage->GetLargestPossibleRegion()); int ctr=0; double sum=0; boost::accumulators::accumulator_set> > acc; for (it.GoToBegin(); !it.IsAtEnd(); ++it) { acc(it.Get()); // currentPixel = it.Get(); // sum+=currentPixel; // ctr+=1; } // res.mean=(double)sum/(double)ctr; res.mean = boost::accumulators::mean(acc); res.variance = boost::accumulators::variance(acc); res.min = boost::accumulators::min(acc); res.max = boost::accumulators::max(acc); res.count = boost::accumulators::count(acc); res.moment = boost::accumulators::moment<2>(acc); std::cout << "sum: " << sum << " N: " << ctr << " mean: " << res.mean << std::endl; } -int main( int argc, char* argv[] ) +void compute_statistics(std::string inputImageFile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic3D.nrrd", + std::string outfname = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic3D_statistics_new.txt", + std::string outfname2 = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic3D_statistics_old.txt", + std::string maskImageFile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic3D_someSegmentation.nrrd", + std::string pFfile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic3D_rectangle.pf" ) { + ofstream outFile; + outFile.open(outfname); + unsigned int timeStep = 0; - std::string inputImageFile; - //inputImageFile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic2DplusT.nrrd"; - //inputImageFile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic3D.nrrd"; - inputImageFile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/leber_ct.pic"; - //inputImageFile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/US4DCyl.nrrd"; - // Load image + mitk::Image::Pointer inputImage = mitk::IOUtil::LoadImage(inputImageFile); mitk::Image::Pointer maskImage; - //std::string maskImageFile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic3D_someSegmentation.nrrd"; - std::string maskImageFile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/leber_ct_segmentation.nrrd"; maskImage = mitk::IOUtil::LoadImage(maskImageFile); - mitk::ImageMaskGenerator::Pointer binaryImageMaskGen = mitk::ImageMaskGenerator::New(); - binaryImageMaskGen->SetImageMask(maskImage); - mitk::HotspotMaskGenerator::Pointer hotSpotMaskGen = mitk::HotspotMaskGenerator::New(); - hotSpotMaskGen->SetImage(inputImage); - hotSpotMaskGen->SetHotspotRadiusInMM(3); - hotSpotMaskGen->SetTimeStep(0); - - //std::string pFfile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic2DplusT_ellipse.pf"; - std::string pFfile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic3D_rectangle.pf"; std::vector loadedObjects; loadedObjects = mitk::IOUtil::Load(pFfile); mitk::BaseData::Pointer pf = loadedObjects[0]; mitk::PlanarFigure::Pointer planarFigure = dynamic_cast(pf.GetPointer()); - mitk::PlanarFigureMaskGenerator::Pointer planarFigMaskExtr = mitk::PlanarFigureMaskGenerator::New(); - planarFigMaskExtr->SetPlanarFigure(planarFigure); - planarFigMaskExtr->SetImage(inputImage); mitk::ImageStatisticsCalculator_v2::Pointer calculator = mitk::ImageStatisticsCalculator_v2::New(); - - std::cout << "calculating statistics (unmasked) itk: " << std::endl; mitk::ImageStatisticsCalculator_v2::StatisticsContainer::Pointer result; - calculator->SetInputImage(inputImage); calculator->SetNBinsForHistogramStatistics(100); - std::cout << std::endl; + long start_time; -// std::cout << "unmasked: " << std::endl; -// calculator->SetMask(nullptr); -// result = calculator->GetStatistics(timeStep); -// result->PrintSelf(); -// std::cout << std::endl; + outFile << "Calculating Statistics on Pic3D" << std::endl << std::endl; -// std::cout << "planarfigure: " << std::endl; -// calculator->SetMask(planarFigMaskExtr.GetPointer()); -// result = calculator->GetStatistics(timeStep, 1); -// result->PrintSelf(); -// std::cout << std::endl; + outFile << "New Image Statistics Calculator" << std::endl; -// std::cout << "ignore pixel value mask: " << std::endl; -// mitk::IgnorePixelMaskGenerator::Pointer ignPixVal = mitk::IgnorePixelMaskGenerator::New(); -// ignPixVal->SetImage(inputImage); -// ignPixVal->SetIgnoredPixelValue(0); -// calculator->SetMask(ignPixVal.GetPointer()); -// result = calculator->GetStatistics(timeStep, 1); -// result->PrintSelf(); -// std::cout << std::endl; + start_time = getTimeInMs(); + outFile << "1) unmasked: " << std::endl; + calculator->SetMask(nullptr); + result = calculator->GetStatistics(timeStep); + outFile << result->GetAsString() << std::endl; + outFile << "new image statistics calculator took: " << getTimeInMs()-start_time << " ms." << std::endl << std::endl; + + start_time = getTimeInMs(); + outFile << "2) planarfigure: " << std::endl; + mitk::PlanarFigureMaskGenerator::Pointer planarFigMaskExtr = mitk::PlanarFigureMaskGenerator::New(); + planarFigMaskExtr->SetPlanarFigure(planarFigure); + planarFigMaskExtr->SetImage(inputImage); + calculator->SetMask(planarFigMaskExtr.GetPointer()); + result = calculator->GetStatistics(timeStep, 1); + outFile << result->GetAsString() << std::endl; + outFile << "new image statistics calculator took: " << getTimeInMs()-start_time << " ms." << std::endl << std::endl; - long start_time; start_time = getTimeInMs(); - std::cout << "image mask: " << std::endl; + outFile << "3) ignore pixel value mask: " << std::endl; + mitk::IgnorePixelMaskGenerator::Pointer ignPixVal = mitk::IgnorePixelMaskGenerator::New(); + ignPixVal->SetImage(inputImage); + ignPixVal->SetIgnoredPixelValue(0); + calculator->SetMask(ignPixVal.GetPointer()); + result = calculator->GetStatistics(timeStep, 1); + outFile << result->GetAsString() << std::endl; + outFile << "new image statistics calculator took: " << getTimeInMs()-start_time << " ms." << std::endl << std::endl; + + start_time = getTimeInMs(); + outFile << "4) image mask: " << std::endl; + mitk::ImageMaskGenerator::Pointer binaryImageMaskGen = mitk::ImageMaskGenerator::New(); + binaryImageMaskGen->SetImageMask(maskImage); calculator->SetMask(binaryImageMaskGen.GetPointer()); result = calculator->GetStatistics(timeStep, 1); - result->PrintSelf(); - std::cout << "new image statistics calculator took: " << getTimeInMs()-start_time << " ms." << std::endl; - std::cout << std::endl; + outFile << result->GetAsString() << std::endl; + outFile << "new image statistics calculator took: " << getTimeInMs()-start_time << " ms." << std::endl << std::endl; + + start_time = getTimeInMs(); + outFile << "5) hotspot mask: " << std::endl; + mitk::HotspotMaskGenerator::Pointer hotSpotMaskGen = mitk::HotspotMaskGenerator::New(); + hotSpotMaskGen->SetImage(inputImage); + hotSpotMaskGen->SetHotspotRadiusInMM(10); + hotSpotMaskGen->SetTimeStep(0); + calculator->SetMask(hotSpotMaskGen.GetPointer()); + result = calculator->GetStatistics(timeStep, 1); + outFile << result->GetAsString() << std::endl; + outFile << "new image statistics calculator took: " << getTimeInMs()-start_time << " ms." << std::endl << std::endl; + //mitk::IOUtil::SaveImage(hotSpotMaskGen->GetMask(), "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic2DplusTHotspot.nrrd"); + + outFile << "6) all time steps (image mask): " << std::endl; + if (inputImage->GetTimeSteps() > 1) + { + start_time = getTimeInMs(); + for (unsigned int i=0; i < inputImage->GetTimeSteps(); i++) + { + outFile << "timestep: " << i << std::endl; + calculator->SetMask(binaryImageMaskGen.GetPointer()); + result = calculator->GetStatistics(i, 1); + outFile << result->GetAsString() << std::endl; + outFile << "new image statistics calculator took: " << getTimeInMs()-start_time << " ms." << std::endl << std::endl; + } + } + else + { + outFile << "Input image has only 1 time step!" << std::endl; + } + + outFile.flush(); + outFile.close(); -// for (unsigned int i=0; i < inputImage->GetTimeSteps(); i++) -// { -// std::cout << "no mask: " << std::endl; -// std::cout << "timestep: " << i << std::endl; -// calculator->SetMask(nullptr); -// result = calculator->GetStatistics(i); -// result->PrintSelf(); -// std::cout << std::endl; -// } - -// std::cout << "hotspot mask: " << std::endl; -// calculator->SetMask(hotSpotMaskGen.GetPointer()); -// result = calculator->GetStatistics(timeStep, 1); -// result->PrintSelf(); -// //mitk::IOUtil::SaveImage(hotSpotMaskGen->GetMask(), "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic2DplusTHotspot.nrrd"); -// std::cout << std::endl; + // ----------------------------------------------------------------------------------------------------------------------------------------------------------- + outFile.open(outfname2); std::cout << std::endl << std::endl << "calculating statistics old imgstatcalc: " << std::endl; + mitk::ImageStatisticsCalculator::Statistics statisticsStruct; + mitk::ImageStatisticsCalculator::Pointer calculator_old; + + start_time = getTimeInMs(); - mitk::ImageStatisticsCalculator::Statistics statisticsStruct; - mitk::ImageStatisticsCalculator::Pointer calculator_old = mitk::ImageStatisticsCalculator::New(); + outFile << "1) unmasked: " << std::endl; + calculator_old = mitk::ImageStatisticsCalculator::New(); + calculator_old->SetHistogramBinSize(100); calculator_old->SetImage(inputImage); - calculator_old->SetImageMask(maskImage); - calculator_old->SetMaskingModeToImage(); + calculator_old->SetMaskingModeToNone(); + calculator_old->ComputeStatistics(timeStep); + statisticsStruct = calculator_old->GetStatistics(timeStep); + outFile << printstats(statisticsStruct) << std::endl; + outFile << "old image statistics calculator took: " << getTimeInMs()-start_time << " ms." << std::endl << std::endl; + + start_time = getTimeInMs(); + outFile << "2) planarFigure: " << std::endl; + calculator_old = mitk::ImageStatisticsCalculator::New(); calculator_old->SetHistogramBinSize(100); -// calculator_old->SetPlanarFigure(planarFigure); -// calculator_old->SetMaskingModeToPlanarFigure(); - //calculator_old->SetMaskingModeToNone(); + calculator_old->SetImage(inputImage); + calculator_old->SetPlanarFigure(planarFigure); + calculator_old->SetMaskingModeToPlanarFigure(); + calculator_old->ComputeStatistics(timeStep); + statisticsStruct = calculator_old->GetStatistics(timeStep); + outFile << printstats(statisticsStruct) << std::endl; + outFile << "old image statistics calculator took: " << getTimeInMs()-start_time << " ms." << std::endl << std::endl; + + start_time = getTimeInMs(); + outFile << "3) IgnorePixelValue: " << std::endl; + calculator_old = mitk::ImageStatisticsCalculator::New(); + calculator_old->SetHistogramBinSize(100); + calculator_old->SetImage(inputImage); + calculator_old->SetMaskingModeToNone(); + calculator_old->SetDoIgnorePixelValue(true); + calculator_old->SetIgnorePixelValue(0); + calculator_old->ComputeStatistics(timeStep); + statisticsStruct = calculator_old->GetStatistics(timeStep); + outFile << printstats(statisticsStruct) << std::endl; + outFile << "old image statistics calculator took: " << getTimeInMs()-start_time << " ms." << std::endl << std::endl; calculator_old->SetDoIgnorePixelValue(false); + + start_time = getTimeInMs(); + outFile << "4) Image Mask: " << std::endl; + calculator_old = mitk::ImageStatisticsCalculator::New(); + calculator_old->SetHistogramBinSize(100); + calculator_old->SetImage(inputImage); + calculator_old->SetImageMask(maskImage); + calculator_old->SetMaskingModeToImage(); calculator_old->ComputeStatistics(timeStep); statisticsStruct = calculator_old->GetStatistics(timeStep); - std::cout << "old image statistics calculator took: " << getTimeInMs()-start_time << " ms." << std::endl; + outFile << printstats(statisticsStruct) << std::endl; + outFile << "old image statistics calculator took: " << getTimeInMs()-start_time << " ms." << std::endl << std::endl; - std::cout << "entropy " << statisticsStruct.GetEntropy() << std::endl; - std::cout << "uniformity " << statisticsStruct.GetUniformity() << std::endl; - std::cout << "upp " << statisticsStruct.GetUPP() << std::endl; - std::cout << "median " << statisticsStruct.GetMedian() << std::endl; + start_time = getTimeInMs(); + outFile << "5) Hotspot Mask: " << std::endl; + calculator_old = mitk::ImageStatisticsCalculator::New(); + calculator_old->SetHistogramBinSize(100); + calculator_old->SetImage(inputImage); + calculator_old->SetMaskingModeToNone(); + calculator_old->SetCalculateHotspot(true); + calculator_old->SetHotspotRadiusInMM(10); + calculator_old->ComputeStatistics(timeStep); + statisticsStruct = calculator_old->GetStatistics(timeStep).GetHotspotStatistics(); + outFile << printstats(statisticsStruct) << std::endl; + outFile << "old image statistics calculator took: " << getTimeInMs()-start_time << " ms." << std::endl << std::endl; + calculator_old->SetCalculateHotspot(false); + + outFile << "6) all time steps (image mask): " << std::endl; + if (inputImage->GetTimeSteps() > 1) + { + start_time = getTimeInMs(); + calculator_old = mitk::ImageStatisticsCalculator::New(); + calculator_old->SetHistogramBinSize(100); + calculator_old->SetImage(inputImage); + calculator_old->SetImageMask(maskImage); + calculator_old->SetMaskingModeToImage(); + for (unsigned int i=0; i < inputImage->GetTimeSteps(); i++) + { + calculator_old->ComputeStatistics(i); + statisticsStruct = calculator_old->GetStatistics(i); + outFile << printstats(statisticsStruct) << std::endl; + outFile << "old image statistics calculator took: " << getTimeInMs()-start_time << " ms." << std::endl << std::endl; + } + } + else + { + outFile << "Input image has only 1 time step!" << std::endl; + } + + outFile.flush(); + outFile.close(); +} + +int main( int argc, char* argv[] ) +{ + std::string inputImageFile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic3D.nrrd"; + std::string outfname = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic3D_statistics_new.txt"; + std::string outfname2 = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic3D_statistics_old.txt"; + std::string maskImageFile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic3D_someSegmentation.nrrd"; + std::string pFfile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic3D_rectangle.pf"; + compute_statistics(inputImageFile, outfname, outfname2, maskImageFile, pFfile); + + inputImageFile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic2DplusT.nrrd"; + outfname = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic2DplusT_statistics_new.txt"; + outfname2 = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic2DplusT_statistics_old.txt"; + maskImageFile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic2DplusT_someSegmentation.nrrd"; + pFfile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic2DplusT_ellipse.pf"; + compute_statistics(inputImageFile, outfname, outfname2, maskImageFile, pFfile); + + inputImageFile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/leber_ct.pic"; + outfname = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/leber_ct_statistics_new.txt"; + outfname2 = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/leber_ct_statistics_old.txt"; + maskImageFile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/leber_ct_segmentation.nrrd"; + pFfile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/leber_ct_PF.pf"; + compute_statistics(inputImageFile, outfname, outfname2, maskImageFile, pFfile); -// std::cout << "calculating statistics boost: " << std::endl; -// statistics_res res; -// AccessByItk_n(inputImage, get_statistics_boost, (res)); -// printstats(res); return EXIT_SUCCESS; } diff --git a/Modules/ImageStatistics/mitkImageStatisticsCalculator2.cpp b/Modules/ImageStatistics/mitkImageStatisticsCalculator2.cpp index 1c018701aa..d1a4e1de3f 100644 --- a/Modules/ImageStatistics/mitkImageStatisticsCalculator2.cpp +++ b/Modules/ImageStatistics/mitkImageStatisticsCalculator2.cpp @@ -1,454 +1,477 @@ #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include "itkImageFileWriter.h" namespace mitk { void ImageStatisticsCalculator_v2::SetInputImage(mitk::Image::Pointer image) { if (image != m_Image) { m_Image = image; m_StatisticsUpdateRequiredByTimeStep.resize(m_Image->GetTimeSteps()); m_StatisticsByTimeStep.resize(m_Image->GetTimeSteps()); this->SetAllStatisticsToUpdateRequired(); } } void ImageStatisticsCalculator_v2::SetMask(mitk::MaskGenerator::Pointer mask) { if (mask != m_MaskGenerator) { m_MaskGenerator = mask; this->SetAllStatisticsToUpdateRequired(); } } void ImageStatisticsCalculator_v2::SetNBinsForHistogramStatistics(unsigned int nBins) { if (nBins != m_nBinsForHistogramStatistics) { m_nBinsForHistogramStatistics = nBins; this->SetAllStatisticsToUpdateRequired(); } } unsigned int ImageStatisticsCalculator_v2::GetNBinsForHistogramStatistics() const { return m_nBinsForHistogramStatistics; } ImageStatisticsCalculator_v2::StatisticsContainer::Pointer ImageStatisticsCalculator_v2::GetStatistics(unsigned int timeStep, unsigned int label) { - std::cout << "getting statistics..." << std::endl; if (timeStep >= m_StatisticsByTimeStep.size()) { throw std::runtime_error("invalid timeStep in ImageStatisticsCalculator_v2::GetStatistics"); } if (m_Image.IsNull()) { throw std::runtime_error("no image"); } if (m_MaskGenerator.IsNotNull()) { m_MaskGenerator->SetTimeStep(timeStep); m_InternalMask = m_MaskGenerator->GetMask(); } ImageTimeSelector::Pointer imgTimeSel = ImageTimeSelector::New(); imgTimeSel->SetInput(m_Image); imgTimeSel->SetTimeNr(timeStep); imgTimeSel->UpdateLargestPossibleRegion(); m_ImageTimeSlice = imgTimeSel->GetOutput(); if (m_StatisticsUpdateRequiredByTimeStep[timeStep]) { - std::cout << "calculating statistics..." << std::endl; - // Calculate statistics with/without mask if (m_MaskGenerator.IsNull()) { // 1) calculate statistics unmasked: // plug image into itkstatisticsimagefilter (will be replaced by my awesome filter later on) // retrieve statistics and save them - std::cout << "unmasked" << std::endl; - AccessByItk_1(m_ImageTimeSlice, InternalCalculateStatisticsUnmasked, timeStep) } else { - std::cout << "masked" << std::endl; - // 2) calculate statistics masked // extract mask image region // plug mask and image into itklabelstatisticsimagefilter AccessByItk_1(m_ImageTimeSlice, InternalCalculateStatisticsMasked, timeStep) } m_StatisticsUpdateRequiredByTimeStep[timeStep]=false; } for (std::vector::iterator it = m_StatisticsByTimeStep[timeStep].begin(); it != m_StatisticsByTimeStep[timeStep].end(); ++it) { StatisticsContainer::Pointer statCont = *it; if (statCont->GetLabel() == label) { return statCont; } } - MITK_ERROR << "Invalid Label: " << label; + MITK_WARN << "Invalid label: " << label << " in time step: " << timeStep; + return StatisticsContainer::New(); } void ImageStatisticsCalculator_v2::SetAllStatisticsToUpdateRequired() { for (unsigned int i = 0; i < m_StatisticsUpdateRequiredByTimeStep.size(); i++) { this->SetStatsTimeStepToUpdateRequired(i); } } void ImageStatisticsCalculator_v2::SetStatsTimeStepToUpdateRequired(unsigned int timeStep) { if (timeStep >= m_StatisticsUpdateRequiredByTimeStep.size()) { throw std::runtime_error("invalid timeStep in ImageStatisticsCalculator_v2::SetStatsTimeStepToUpdateRequired"); } m_StatisticsUpdateRequiredByTimeStep[timeStep] = true; m_StatisticsByTimeStep[timeStep].resize(0); } template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator_v2::InternalCalculateStatisticsUnmasked( typename itk::Image< TPixel, VImageDimension >* image, unsigned int timeStep) { typedef typename itk::Image< TPixel, VImageDimension > ImageType; typedef typename itk::ExtendedStatisticsImageFilter2 ImageStatisticsFilterType; typedef typename itk::MinMaxImageFilterWithIndex MinMaxFilterType; StatisticsContainer::Pointer statisticsResult = StatisticsContainer::New(); typename 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); //statisticsFilter->SetHistogramParameters(m_nBinsForHistogramStatistics, imgMinMaxFilter->GetMinimum(), imgMinMaxFilter->GetMaximum()); statisticsFilter->SetHistogramParameters(m_nBinsForHistogramStatistics, minval, maxval); try { statisticsFilter->Update(); } catch (const itk::ExceptionObject& e) { mitkThrow() << "Image statistics calculation failed due to following ITK Exception: \n " << e.what(); } // no mask, therefore just one label = the whole image m_StatisticsByTimeStep[timeStep].resize(1); statisticsResult->SetLabel(1); statisticsResult->SetN(image->GetLargestPossibleRegion().GetNumberOfPixels()); statisticsResult->SetMean(statisticsFilter->GetMean()); statisticsResult->SetMin(statisticsFilter->GetMinimum()); statisticsResult->SetMax(statisticsFilter->GetMaximum()); statisticsResult->SetVariance(statisticsFilter->GetVariance()); statisticsResult->SetStd(statisticsFilter->GetSigma()); statisticsResult->SetSkewness(statisticsFilter->GetSkewness()); statisticsResult->SetKurtosis(statisticsFilter->GetKurtosis()); statisticsResult->SetRMS(std::sqrt(std::pow(statisticsFilter->GetMean(), 2.) + statisticsFilter->GetVariance())); // variance = sigma^2 statisticsResult->SetMPP(statisticsFilter->GetMPP()); statisticsResult->SetEntropy(statisticsFilter->GetEntropy()); statisticsResult->SetMedian(statisticsFilter->GetMedian()); statisticsResult->SetUniformity(statisticsFilter->GetUniformity()); statisticsResult->SetUPP(statisticsFilter->GetUPP()); statisticsResult->SetHistogram(statisticsFilter->GetHistogram()); m_StatisticsByTimeStep[timeStep][0] = statisticsResult; } template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator_v2::InternalCalculateStatisticsMasked( typename itk::Image< TPixel, VImageDimension >* image, unsigned int timeStep) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< unsigned short, VImageDimension > MaskType; typedef typename MaskType::PixelType LabelPixelType; typedef itk::ExtendedLabelStatisticsImageFilter2< ImageType, MaskType > ImageStatisticsFilterType; typedef MaskUtilities< TPixel, VImageDimension > MaskUtilType; typedef typename itk::MinMaxLabelImageFilterWithIndex MinMaxLabelFilterType; typedef typename ImageType::PixelType InputImgPixelType; // maskImage has to have the same dimension as image typename MaskType::Pointer maskImage = MaskType::New(); maskImage = ImageToItkImage< unsigned short, VImageDimension >(m_InternalMask); typename MaskUtilType::Pointer maskUtil = MaskUtilType::New(); maskUtil->SetImage(image); maskUtil->SetMask(maskImage.GetPointer()); // if mask is smaller than image, extract the image region where the mask is typename ImageType::Pointer adaptedImage = ImageType::New(); adaptedImage = maskUtil->ExtractMaskImageRegion(); // this also checks mask sanity // typename itk::ImageFileWriter::Pointer imgFileWriter = itk::ImageFileWriter::New(); // imgFileWriter->SetInput(adaptedImage); // imgFileWriter->SetFileName("/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/ct_leber_adapted_image.nrrd"); // imgFileWriter->Update(); // typename itk::ImageFileWriter::Pointer imgFileWriter2 = itk::ImageFileWriter::New(); // imgFileWriter2->SetInput(maskImage); // imgFileWriter2->SetFileName("/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/ct_leber_adapted_image_mask.nrrd"); // imgFileWriter2->Update(); // 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 typedef typename std::map MapType; typedef typename std::pair PairType; std::vector relevantLabels = minMaxFilter->GetRelevantLabels(); MapType minVals; MapType maxVals; std::map nBins; for (LabelPixelType label:relevantLabels) { minVals.insert(PairType(label, minMaxFilter->GetMin(label))); maxVals.insert(PairType(label, minMaxFilter->GetMax(label))); nBins.insert(typename std::pair(label, m_nBinsForHistogramStatistics)); } // minVal = minMaxFilter->GetGlobalMin(); // maxVal = minMaxFilter->GetGlobalMax(); typename ImageStatisticsFilterType::Pointer imageStatisticsFilter = ImageStatisticsFilterType::New(); imageStatisticsFilter->SetDirectionTolerance(0.001); imageStatisticsFilter->SetCoordinateTolerance(0.001); imageStatisticsFilter->SetInput(adaptedImage); imageStatisticsFilter->SetLabelInput(maskImage); //imageStatisticsFilter->SetHistogramParameters(m_nBinsForHistogramStatistics, floor(minVal), ceil(maxVal)); imageStatisticsFilter->SetHistogramParametersForLabels(nBins, minVals, maxVals); imageStatisticsFilter->Update(); std::list labels = imageStatisticsFilter->GetRelevantLabels(); std::list::iterator it = labels.begin(); m_StatisticsByTimeStep[timeStep].resize(0); while(it != labels.end()) { StatisticsContainer::Pointer statisticsResult = StatisticsContainer::New(); // find min, max, minindex and maxindex // make sure to only look in the masked region, use a masker for this vnl_vector minIndex, maxIndex; typename ImageType::IndexType tmpMinIndex = minMaxFilter->GetMinIndex(*it); typename ImageType::IndexType tmpMaxIndex = minMaxFilter->GetMaxIndex(*it); minIndex.set_size(tmpMaxIndex.GetIndexDimension()); maxIndex.set_size(tmpMaxIndex.GetIndexDimension()); for (unsigned int i=0; i < tmpMaxIndex.GetIndexDimension(); i++) { minIndex[i] = tmpMinIndex[i] + (maskImage->GetOrigin()[i] - image->GetOrigin()[i]) / (double) maskImage->GetSpacing()[i]; maxIndex[i] = tmpMaxIndex[i] + (maskImage->GetOrigin()[i] - image->GetOrigin()[i]) / (double) maskImage->GetSpacing()[i]; } statisticsResult->SetMinIndex(minIndex); statisticsResult->SetMaxIndex(maxIndex); // just debug assert(minMaxFilter->GetMax(*it) == imageStatisticsFilter->GetMaximum(*it)); assert(minMaxFilter->GetMin(*it) == imageStatisticsFilter->GetMinimum(*it)); statisticsResult->SetN(imageStatisticsFilter->GetSum(*it) / (double) imageStatisticsFilter->GetMean(*it)); statisticsResult->SetMean(imageStatisticsFilter->GetMean(*it)); statisticsResult->SetMin(imageStatisticsFilter->GetMinimum(*it)); statisticsResult->SetMax(imageStatisticsFilter->GetMaximum(*it)); statisticsResult->SetVariance(imageStatisticsFilter->GetVariance(*it)); statisticsResult->SetStd(imageStatisticsFilter->GetSigma(*it)); statisticsResult->SetSkewness(imageStatisticsFilter->GetSkewness(*it)); statisticsResult->SetKurtosis(imageStatisticsFilter->GetKurtosis(*it)); statisticsResult->SetRMS(std::sqrt(std::pow(imageStatisticsFilter->GetMean(*it), 2.) + imageStatisticsFilter->GetVariance(*it))); // variance = sigma^2 statisticsResult->SetMPP(imageStatisticsFilter->GetMPP(*it)); statisticsResult->SetLabel(*it); statisticsResult->SetEntropy(imageStatisticsFilter->GetEntropy(*it)); statisticsResult->SetMedian(imageStatisticsFilter->GetMedian(*it)); statisticsResult->SetUniformity(imageStatisticsFilter->GetUniformity(*it)); statisticsResult->SetUPP(imageStatisticsFilter->GetUPP(*it)); statisticsResult->SetHistogram(imageStatisticsFilter->GetHistogram(*it)); m_StatisticsByTimeStep[timeStep].push_back(statisticsResult); ++it; } } ImageStatisticsCalculator_v2::StatisticsContainer::StatisticsContainer(): m_N(std::numeric_limits::max()), m_Mean(std::numeric_limits::max()), m_Min(std::numeric_limits::min()), m_Max(std::numeric_limits::max()), m_Std(std::numeric_limits::max()), m_Variance(std::numeric_limits::max()), m_Skewness(std::numeric_limits::max()), m_Kurtosis(std::numeric_limits::max()), m_RMS(std::numeric_limits::max()), m_MPP(std::numeric_limits::max()), m_Median(std::numeric_limits::max()), m_Uniformity(std::numeric_limits::max()), m_UPP(std::numeric_limits::max()), m_Entropy(std::numeric_limits::max()) { m_minIndex.set_size(0); m_maxIndex.set_size(0); } ImageStatisticsCalculator_v2::statisticsMapType ImageStatisticsCalculator_v2::StatisticsContainer::GetStatisticsAsMap() { ImageStatisticsCalculator_v2::statisticsMapType statisticsAsMap; statisticsAsMap["N"] = m_N; statisticsAsMap["Mean"] = m_Mean; statisticsAsMap["Min"] = m_Min; statisticsAsMap["Max"] = m_Max; statisticsAsMap["StandardDeviation"] = m_Std; statisticsAsMap["Variance"] = m_Variance; statisticsAsMap["Skewness"] = m_Skewness; statisticsAsMap["Kurtosis"] = m_Kurtosis; statisticsAsMap["RMS"] = m_RMS; statisticsAsMap["MPP"] = m_MPP; statisticsAsMap["Median"] = m_Median; statisticsAsMap["Uniformity"] = m_Uniformity; statisticsAsMap["UPP"] = m_UPP; statisticsAsMap["Entropy"] = m_Entropy; return statisticsAsMap; } void ImageStatisticsCalculator_v2::StatisticsContainer::Reset() { m_N = std::numeric_limits::max(); m_Mean = std::numeric_limits::max(); m_Min = std::numeric_limits::max(); m_Max = std::numeric_limits::max(); m_Std = std::numeric_limits::max(); m_Variance = std::numeric_limits::max(); m_Skewness = std::numeric_limits::max(); m_Kurtosis = std::numeric_limits::max(); m_RMS = std::numeric_limits::max(); m_MPP = std::numeric_limits::max(); m_Median = std::numeric_limits::max(); m_Uniformity = std::numeric_limits::max(); m_UPP = std::numeric_limits::max(); m_Entropy = std::numeric_limits::max(); } void ImageStatisticsCalculator_v2::StatisticsContainer::PrintSelf() { ImageStatisticsCalculator_v2::statisticsMapType statMap = this->GetStatisticsAsMap(); // print all map key value pairs // const auto& val:statMap for (auto it = statMap.begin(); it != statMap.end(); ++it) { std::cout << it->first << ": " << it->second << std::endl; } // print the min and max index std::cout << "Min Index:" << std::endl; for (auto it = this->GetMinIndex().begin(); it != this->GetMinIndex().end(); it++) { std::cout << *it << " "; } std::cout << std::endl; // print the min and max index std::cout << "Max Index:" << std::endl; for (auto it = this->GetMaxIndex().begin(); it != this->GetMaxIndex().end(); it++) { std::cout << *it << " "; } std::cout << std::endl; } + std::string ImageStatisticsCalculator_v2::StatisticsContainer::GetAsString() + { + std::string res = ""; + ImageStatisticsCalculator_v2::statisticsMapType statMap = this->GetStatisticsAsMap(); + // print all map key value pairs + // const auto& val:statMap + for (auto it = statMap.begin(); it != statMap.end(); ++it) + { + res += std::string(it->first) + ": " + std::to_string(it->second) + "\n"; + } + + // print the min and max index + res += "Min Index:" + std::string("\n"); + for (auto it = this->GetMinIndex().begin(); it != this->GetMinIndex().end(); it++) + { + res += std::to_string(*it) + std::string(" "); + } + res += "\n"; + + // print the min and max index + res += "Max Index:" + std::string("\n"); + for (auto it = this->GetMaxIndex().begin(); it != this->GetMaxIndex().end(); it++) + { + res += std::to_string(*it) + " "; + } + res += "\n"; + return res; + } + } diff --git a/Modules/ImageStatistics/mitkImageStatisticsCalculator2.h b/Modules/ImageStatistics/mitkImageStatisticsCalculator2.h index fdd8af3207..f06b0109e7 100644 --- a/Modules/ImageStatistics/mitkImageStatisticsCalculator2.h +++ b/Modules/ImageStatistics/mitkImageStatisticsCalculator2.h @@ -1,333 +1,336 @@ #ifndef MITKIMAGESTATISTICSCALCULATOR2 #define MITKIMAGESTATISTICSCALCULATOR2 #include #include #include #include #include #include #include /*#define mitkGetRefConstMacro(name, type) \ const type& Get##name() const \ { \ return &m_##name; \ } \ \ type& Get##name() \ { \ return &m_##name; \ } \ */ namespace mitk { class MITKIMAGESTATISTICS_EXPORT ImageStatisticsCalculator_v2: public itk::Object { public: /** Standard Self typedef */ typedef ImageStatisticsCalculator_v2 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; class StatisticsContainer: public itk::Object { public: /** Standard Self typedef */ typedef StatisticsContainer 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(StatisticsContainer, itk::Object) typedef double RealType; statisticsMapType GetStatisticsAsMap(); void Reset(); void SetN(long n) { m_N = n; } const long& GetN() const { return m_N; } void SetMean(RealType mean) { m_Mean = mean; } const RealType& GetMean() const { return m_Mean; } void SetVariance(RealType variance) { m_Variance = variance; } const RealType& GetVariance() const { return m_Variance; } void SetStd(RealType std) { m_Std = std; } const RealType& GetStd() const { return m_Std; } void SetMin(RealType minVal) { m_Min = minVal; } const RealType& GetMin() const { return m_Min; } void SetMax(RealType maxVal) { m_Max = maxVal; } const RealType& GetMax() const { return m_Max; } void SetRMS(RealType rms) { m_RMS = rms; } const RealType& GetRMS() const { return m_RMS; } void SetSkewness(RealType skewness) { m_Skewness = skewness; } const RealType& GetSkewness() const { return m_Skewness; } void SetKurtosis(RealType kurtosis) { m_Kurtosis = kurtosis; } const RealType& GetKurtosis() const { return m_Kurtosis; } void SetMPP(RealType mpp) { m_MPP = mpp; } const RealType& GetMPP() const { return m_MPP; } void SetLabel(unsigned int label) { m_Label = label; } const unsigned int& GetLabel() const { return m_Label; } void SetMinIndex(vnl_vector& minIndex) { m_minIndex = minIndex; } const vnl_vector& GetMinIndex() const { return m_minIndex; } void SetMaxIndex(vnl_vector& maxIndex) { m_maxIndex = maxIndex; } const vnl_vector& GetMaxIndex() const { return m_maxIndex; } void SetHistogram(HistogramType::Pointer hist) { if (m_Histogram != hist) { m_Histogram = hist; } } const HistogramType::Pointer GetHistogram() const { return m_Histogram; } void SetEntropy(RealType entropy) { m_Entropy = entropy; } const RealType & GetEntropy() const { return m_Entropy; } void SetMedian(RealType median) { m_Median = median; } const RealType & GetMedian() const { return m_Median; } void SetUniformity(RealType uniformity) { m_Uniformity = uniformity; } const RealType & GetUniformity() const { return m_Uniformity; } void SetUPP(RealType upp) { m_UPP = upp; } const RealType & GetUPP() const { return m_UPP; } void PrintSelf(); + std::string GetAsString(); + + protected: StatisticsContainer(); private: // not pretty, is temporary long m_N; RealType m_Mean, m_Min, m_Max, m_Std, m_Variance; 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 SetInputImage(mitk::Image::Pointer image); void SetMask(mitk::MaskGenerator::Pointer mask); void SetNBinsForHistogramStatistics(unsigned int nBins); unsigned int GetNBinsForHistogramStatistics() const; StatisticsContainer::Pointer GetStatistics(unsigned int timeStep=0, unsigned int label=1); protected: ImageStatisticsCalculator_v2(){ m_nBinsForHistogramStatistics = 100; }; private: void SetAllStatisticsToUpdateRequired(); void SetStatsTimeStepToUpdateRequired(unsigned int timeStep); template < typename TPixel, unsigned int VImageDimension > void InternalCalculateStatisticsUnmasked( typename itk::Image< TPixel, VImageDimension >* image, unsigned int timeStep); template < typename TPixel, unsigned int VImageDimension > typename HistogramType::Pointer InternalCalculateHistogramUnmasked( typename itk::Image< TPixel, VImageDimension >* image, double minVal, double maxVal); template < typename TPixel, unsigned int VImageDimension > void InternalCalculateStatisticsMasked( typename itk::Image< TPixel, VImageDimension >* image, unsigned int timeStep); // template < typename TPixel, unsigned int VImageDimension > void GetMinAndMaxValue( // typename itk::Image< TPixel, VImageDimension >* inputImage, // double &minVal, // double &maxVal // ); std::string GetNameOfClass() { return std::string("ImageStatisticsCalculator_v2"); } mitk::Image::Pointer m_Image; mitk::Image::Pointer m_ImageTimeSlice; mitk::MaskGenerator::Pointer m_MaskGenerator; mitk::Image::Pointer m_InternalMask; unsigned int m_nBinsForHistogramStatistics; std::vector m_StatisticsUpdateRequiredByTimeStep; // holds which time steps are valid and which ones have to be recalculated std::vector> m_StatisticsByTimeStep; }; } #endif // MITKIMAGESTATISTICSCALCULATOR2