diff --git a/Modules/DiffusionImaging/MiniApps/HotspotMiniApp.cpp b/Modules/DiffusionImaging/MiniApps/HotspotMiniApp.cpp index 903dff8eae..3dae227d76 100644 --- a/Modules/DiffusionImaging/MiniApps/HotspotMiniApp.cpp +++ b/Modules/DiffusionImaging/MiniApps/HotspotMiniApp.cpp @@ -1,103 +1,103 @@ /*=================================================================== 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 "mitkIOUtil.h" #include #include #include #include int main( int argc, char* argv[] ) { mitkCommandLineParser parser; parser.setTitle("Test basic hotspot functionality"); parser.setCategory("Preprocessing Tools"); parser.setDescription(""); parser.setContributor("MBI"); parser.setArgumentPrefix("--", "-"); parser.addArgument("help", "h", mitkCommandLineParser::String, "Help:", "Show this help text"); parser.addArgument("input", "i", mitkCommandLineParser::InputFile, "Input:", "input image", us::Any(),false); parser.addArgument("mask", "m", mitkCommandLineParser::InputFile, "Mask:", "mask image / roi image denotin area on which statistics are calculated", us::Any(),true); parser.addArgument("out", "o", mitkCommandLineParser::OutputFile, "Output", "output file (default: hotspotMiniApp_output.nrrd)", us::Any()); std::cout << "test...." << std::endl; std::map parsedArgs = parser.parseArguments(argc, argv); std::cout << "parsedArgs.size()= " << parsedArgs.size() << std::endl; if (parsedArgs.size()==0 || parsedArgs.count("help") || parsedArgs.count("h")) { // std::cout << "\n\n MiniApp Description: \nCalculates and saves the hotspot of an image." << endl; // std::cout << "Output is written to the designated output file" << endl; // std::cout << parser.helpText(); // return EXIT_SUCCESS; } // Parameters: std::string inputImageFile; if (!parsedArgs.count("i") && !parsedArgs.count("input")) { inputImageFile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic3D.nrrd"; } else { inputImageFile = us::any_cast(parsedArgs["input"]); } mitk::Image::Pointer maskImage; if (parsedArgs.count("mask") || parsedArgs.count("m")) { std::string maskImageFile = us::any_cast(parsedArgs["mask"]); maskImage = mitk::IOUtil::LoadImage(maskImageFile); } std::string outFile; if (parsedArgs.count("out") || parsedArgs.count("o") ) outFile = us::any_cast(parsedArgs["out"]); else outFile = "hotspotMiniApp_output.nrrd"; // Load image and mask mitk::Image::Pointer inputImage = mitk::IOUtil::LoadImage(inputImageFile); // Calculate statistics mitk::HotspotMaskGenerator::Pointer hotspotCalc = mitk::HotspotMaskGenerator::New(); - hotspotCalc->SetImage(inputImage); + hotspotCalc->SetInputImage(inputImage); hotspotCalc->SetHotspotRadiusInMM(10); hotspotCalc->SetTimeStep(0); mitk::Image::Pointer outImage; try { outImage = hotspotCalc->GetMask(); } catch( const itk::ExceptionObject& e) { MITK_ERROR << "Failed - ITK Exception:" << e.what(); return -1; } if (outImage != nullptr) { mitk::IOUtil::SaveImage(outImage, outFile); } return EXIT_SUCCESS; } diff --git a/Modules/DiffusionImaging/MiniApps/ImageStatisticsMiniApp.cpp b/Modules/DiffusionImaging/MiniApps/ImageStatisticsMiniApp.cpp index c8076c8519..f2bef60001 100644 --- a/Modules/DiffusionImaging/MiniApps/ImageStatisticsMiniApp.cpp +++ b/Modules/DiffusionImaging/MiniApps/ImageStatisticsMiniApp.cpp @@ -1,218 +1,218 @@ /*=================================================================== 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 "mitkImageStatisticsCalculator.h" #include "mitkIOUtil.h" #include #include #include #include #include "mitkImageAccessByItk.h" #include #include #include #include #include #include #include #include #include #include struct statistics_res{ double mean, variance, min, max, count, moment; }; void printstats(statistics_res s) { 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; } 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()); TPixel currentPixel; 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[] ) { mitkCommandLineParser parser; parser.setTitle("Extract Image Statistics"); parser.setCategory("Preprocessing Tools"); parser.setDescription(""); parser.setContributor("MBI"); parser.setArgumentPrefix("--", "-"); parser.addArgument("help", "h", mitkCommandLineParser::String, "Help:", "Show this help text"); parser.addArgument("input", "i", mitkCommandLineParser::InputFile, "Input:", "input image", us::Any(),false); parser.addArgument("mask", "m", mitkCommandLineParser::InputFile, "Mask:", "mask image / roi image denotin area on which statistics are calculated", us::Any(),true); parser.addArgument("out", "o", mitkCommandLineParser::OutputFile, "Output", "output file (default: filenameOfRoi.nrrd_statistics.txt)", us::Any()); std::cout << "test...." << std::endl; std::map parsedArgs = parser.parseArguments(argc, argv); std::cout << "parsedArgs.size()= " << parsedArgs.size() << std::endl; if (parsedArgs.size()==0 || parsedArgs.count("help") || parsedArgs.count("h")) { std::cout << "\n\n MiniApp Description: \nCalculates statistics on the supplied image using given mask." << endl; std::cout << "Output is written to the designated output file in this order:" << endl; std::cout << "Mean, Standard Deviation, RMS, Max, Min, Number of Voxels, Volume [mm3]" << endl; std::cout << "\n\n Parameters:"<< endl; std::cout << parser.helpText(); return EXIT_SUCCESS; } // Parameters: bool ignoreZeroValues = false; unsigned int timeStep = 0; std::string inputImageFile = us::any_cast(parsedArgs["input"]); mitk::Image::Pointer maskImage; if (parsedArgs.count("mask") || parsedArgs.count("m")) { std::string maskImageFile = us::any_cast(parsedArgs["mask"]); maskImage = mitk::IOUtil::LoadImage(maskImageFile); } std::string outFile; if (parsedArgs.count("out") || parsedArgs.count("o") ) outFile = us::any_cast(parsedArgs["out"]); else outFile = inputImageFile + "_statistics.txt"; // Load image and mask mitk::Image::Pointer inputImage = mitk::IOUtil::LoadImage(inputImageFile); // Calculate statistics mitk::ImageStatisticsCalculator::StatisticsContainer::Pointer statisticsStruct; mitk::ImageStatisticsCalculator::Pointer calculator = mitk::ImageStatisticsCalculator::New(); try { calculator->SetInputImage(inputImage); if (parsedArgs.count("mask") || parsedArgs.count("m")) { mitk::ImageMaskGenerator::Pointer imgMask = mitk::ImageMaskGenerator::New(); imgMask->SetImageMask(maskImage); imgMask->SetTimeStep(timeStep); calculator->SetMask(imgMask.GetPointer()); } else { calculator->SetMask(nullptr); } } catch( const itk::ExceptionObject& e) { MITK_ERROR << "Statistic Calculation Failed - ITK Exception:" << e.what(); return -1; } if (ignoreZeroValues) { // TODO, cannot have more than one mask, using ignore pixel value will override the image mask :-c mitk::IgnorePixelMaskGenerator::Pointer ignorePixelMask = mitk::IgnorePixelMaskGenerator::New(); - ignorePixelMask->SetImage(inputImage); + ignorePixelMask->SetInputImage(inputImage); ignorePixelMask->SetTimeStep(timeStep); ignorePixelMask->SetIgnoredPixelValue(0); calculator->SetMask(ignorePixelMask.GetPointer()); } std::cout << "calculating statistics itk: " << std::endl; try { statisticsStruct = calculator->GetStatistics(timeStep); } catch ( mitk::Exception& e) { MITK_ERROR<< "MITK Exception: " << e.what(); return -1; } // Calculate Volume double volume = 0; const mitk::BaseGeometry *geometry = inputImage->GetGeometry(); if ( geometry != NULL ) { const mitk::Vector3D &spacing = inputImage->GetGeometry()->GetSpacing(); volume = spacing[0] * spacing[1] * spacing[2] * (double) statisticsStruct->GetN(); } // Write Results to file std::ofstream output; output.open(outFile.c_str()); output << statisticsStruct->GetMean() << " , "; output << statisticsStruct->GetStd() << " , "; output << statisticsStruct->GetRMS() << " , "; output << statisticsStruct->GetMax() << " , "; output << statisticsStruct->GetMin() << " , "; output << statisticsStruct->GetN() << " , "; output << volume << "\n"; output.flush(); output.close(); 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/DiffusionImaging/MiniApps/ImageStatisticsMiniapp_v2.cpp b/Modules/DiffusionImaging/MiniApps/ImageStatisticsMiniapp_v2.cpp index b84b315f73..a81d680524 100644 --- a/Modules/DiffusionImaging/MiniApps/ImageStatisticsMiniapp_v2.cpp +++ b/Modules/DiffusionImaging/MiniApps/ImageStatisticsMiniapp_v2.cpp @@ -1,360 +1,360 @@ /*=================================================================== 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 "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; } //std::string printstats(mitk::ImageStatisticsCalculator::Statistics statisticsStruct) //{ // 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; } 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; mitk::Image::Pointer inputImage = mitk::IOUtil::LoadImage(inputImageFile); mitk::Image::Pointer maskImage; maskImage = mitk::IOUtil::LoadImage(maskImageFile); std::vector loadedObjects; loadedObjects = mitk::IOUtil::Load(pFfile); mitk::BaseData::Pointer pf = loadedObjects[0]; mitk::PlanarFigure::Pointer planarFigure = dynamic_cast(pf.GetPointer()); mitk::ImageStatisticsCalculator::Pointer calculator = mitk::ImageStatisticsCalculator::New(); mitk::ImageStatisticsCalculator::StatisticsContainer::Pointer result; calculator->SetInputImage(inputImage); calculator->SetNBinsForHistogramStatistics(100); long start_time; outFile << "Calculating Statistics on Pic3D" << std::endl << std::endl; outFile << "New Image Statistics Calculator" << 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); + planarFigMaskExtr->SetInputImage(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; start_time = getTimeInMs(); outFile << "3) ignore pixel value mask: " << std::endl; mitk::IgnorePixelMaskGenerator::Pointer ignPixVal = mitk::IgnorePixelMaskGenerator::New(); - ignPixVal->SetImage(inputImage); + ignPixVal->SetInputImage(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); 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->SetInputImage(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(); // ----------------------------------------------------------------------------------------------------------------------------------------------------------- // 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(); // outFile << "1) unmasked: " << std::endl; // calculator_old = mitk::ImageStatisticsCalculator::New(); // calculator_old->SetHistogramBinSize(100); // calculator_old->SetImage(inputImage); // 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->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); // outFile << printstats(statisticsStruct) << std::endl; // outFile << "old image statistics calculator took: " << getTimeInMs()-start_time << " ms." << std::endl << 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); return EXIT_SUCCESS; } diff --git a/Modules/DiffusionImaging/MiniApps/PlanarFigureMaskCreatorMiniapp.cpp b/Modules/DiffusionImaging/MiniApps/PlanarFigureMaskCreatorMiniapp.cpp index 585a20d75a..ec163dc9e1 100644 --- a/Modules/DiffusionImaging/MiniApps/PlanarFigureMaskCreatorMiniapp.cpp +++ b/Modules/DiffusionImaging/MiniApps/PlanarFigureMaskCreatorMiniapp.cpp @@ -1,131 +1,131 @@ /*=================================================================== 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 "mitkIOUtil.h" #include #include #include #include #include #include int main( int argc, char* argv[] ) { mitkCommandLineParser parser; parser.setTitle("Test basic hotspot functionality"); parser.setCategory("Preprocessing Tools"); parser.setDescription(""); parser.setContributor("MBI"); parser.setArgumentPrefix("--", "-"); parser.addArgument("help", "h", mitkCommandLineParser::String, "Help:", "Show this help text"); parser.addArgument("input", "i", mitkCommandLineParser::InputFile, "Input:", "input image", us::Any(),false); parser.addArgument("planarfigure", "p", mitkCommandLineParser::InputFile, "PlanarFigure:", "mask image / roi image denotin area on which statistics are calculated", us::Any(),false); parser.addArgument("out", "o", mitkCommandLineParser::OutputFile, "Output", "output file (default: hotspotMiniApp_output.nrrd)", us::Any()); std::cout << "test...." << std::endl; std::map parsedArgs = parser.parseArguments(argc, argv); std::cout << "parsedArgs.size()= " << parsedArgs.size() << std::endl; if (parsedArgs.size()==0 || parsedArgs.count("help") || parsedArgs.count("h")) { // std::cout << "\n\n MiniApp Description: \nCalculates and saves the hotspot of an image." << endl; // std::cout << "Output is written to the designated output file" << endl; // std::cout << parser.helpText(); // return EXIT_SUCCESS; } // Parameters: std::string inputImageFile; if (!parsedArgs.count("i") && !parsedArgs.count("input")) { inputImageFile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic3D.nrrd"; } else { inputImageFile = us::any_cast(parsedArgs["input"]); } std::string inputPFFile; if (parsedArgs.count("planarfigure") || parsedArgs.count("p")) { inputPFFile = us::any_cast(parsedArgs["planarfigure"]); } else { inputPFFile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic3D_rectangle.pf"; } mitk::PlanarFigure::Pointer planarFigure; try { std::vector loadedObjects; // try except loadedObjects = mitk::IOUtil::Load(inputPFFile); mitk::BaseData::Pointer pf = loadedObjects[0]; planarFigure = dynamic_cast(pf.GetPointer()); if (planarFigure.IsNull()) { MITK_ERROR << "something went wrong"; return -1; } } catch (const itk::ExceptionObject& e) { MITK_ERROR << "Failed: " << e.what(); return -1; } std::string outFile; if (parsedArgs.count("out") || parsedArgs.count("o") ) outFile = us::any_cast(parsedArgs["out"]); else outFile = "planarFigureExtraction_output.nrrd"; // Load image and mask mitk::Image::Pointer inputImage = mitk::IOUtil::LoadImage(inputImageFile); // Calculate statistics mitk::PlanarFigureMaskGenerator::Pointer planarFigMaskExtr = mitk::PlanarFigureMaskGenerator::New(); planarFigMaskExtr->SetPlanarFigure(planarFigure); - planarFigMaskExtr->SetImage(inputImage); + planarFigMaskExtr->SetInputImage(inputImage); mitk::Image::Pointer outImage; try { outImage = planarFigMaskExtr->GetMask(); } catch( const itk::ExceptionObject& e) { MITK_ERROR << "Failed - ITK Exception:" << e.what(); return -1; } mitk::BaseGeometry *imageGeo = outImage->GetGeometry(); std::cout << "origin: " << imageGeo->GetOrigin()[0] << " " << imageGeo->GetOrigin()[1] << " " << imageGeo->GetOrigin()[2] << std::endl; if (outImage != nullptr) { mitk::IOUtil::SaveImage(outImage, outFile); } return EXIT_SUCCESS; } diff --git a/Modules/ImageStatistics/Testing/mitkImageStatisticsCalculatorTest.cpp b/Modules/ImageStatistics/Testing/mitkImageStatisticsCalculatorTest.cpp index 481c456506..b55ebd8dbf 100644 --- a/Modules/ImageStatistics/Testing/mitkImageStatisticsCalculatorTest.cpp +++ b/Modules/ImageStatistics/Testing/mitkImageStatisticsCalculatorTest.cpp @@ -1,542 +1,542 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkImageStatisticsCalculator.h" #include #include #include #include #include #include #include #include #include #include /** * \brief Test class for mitkImageStatisticsCalculator * * This test covers: * - instantiation of an ImageStatisticsCalculator class * - correctness of statistics when using PlanarFigures for masking */ class mitkImageStatisticsCalculatorTestSuite : public mitk::TestFixture { CPPUNIT_TEST_SUITE(mitkImageStatisticsCalculatorTestSuite); MITK_TEST(TestUninitializedImage); MITK_TEST(TestCase1); MITK_TEST(TestCase2); MITK_TEST(TestCase3); MITK_TEST(TestCase4); MITK_TEST(TestCase5); MITK_TEST(TestCase6); MITK_TEST(TestCase7); MITK_TEST(TestCase8); MITK_TEST(TestCase9); MITK_TEST(TestCase10); MITK_TEST(TestCase11); MITK_TEST(TestCase12); MITK_TEST(TestImageMaskingEmpty); MITK_TEST(TestImageMaskingNonEmpty); MITK_TEST(TestRecomputeOnModifiedMask); CPPUNIT_TEST_SUITE_END(); public: void setUp() override; void TestUninitializedImage(); void TestCase1(); void TestCase2(); void TestCase3(); void TestCase4(); void TestCase5(); void TestCase6(); void TestCase7(); void TestCase8(); void TestCase9(); void TestCase10(); void TestCase11(); void TestCase12(); void TestImageMaskingEmpty(); void TestImageMaskingNonEmpty(); void TestRecomputeOnModifiedMask(); private: mitk::Image::Pointer m_Image; mitk::PlaneGeometry::Pointer m_Geometry; // calculate statistics for the given image and planarpolygon const mitk::ImageStatisticsCalculator::StatisticsContainer::Pointer ComputeStatistics( mitk::Image::Pointer image, mitk::PlanarFigure::Pointer polygon ); // calculate statistics for the given image and planarpolygon const mitk::ImageStatisticsCalculator::StatisticsContainer::Pointer ComputeStatistics( mitk::Image::Pointer image, mitk::Image::Pointer image_mask ); void VerifyStatistics(mitk::ImageStatisticsCalculator::StatisticsContainer::Pointer stats, double testMean, double testSD, double testMedian=0); }; void mitkImageStatisticsCalculatorTestSuite::setUp() { std::string filename = this->GetTestDataFilePath("ImageStatistics/testimage.dcm"); if (filename.empty()) { MITK_TEST_FAILED_MSG( << "Could not find test file" ) } MITK_TEST_OUTPUT(<< "Loading test image '" << filename << "'") m_Image = mitk::IOUtil::LoadImage(filename); MITK_TEST_CONDITION_REQUIRED( m_Image.IsNotNull(), "Loaded an mitk::Image" ); m_Geometry = m_Image->GetSlicedGeometry()->GetPlaneGeometry(0); MITK_TEST_CONDITION_REQUIRED( m_Geometry.IsNotNull(), "Getting image geometry" ) } void mitkImageStatisticsCalculatorTestSuite::TestCase1() { /***************************** * one whole white pixel * -> mean of 255 expected ******************************/ MITK_INFO << std::endl << "Test case 1:"; mitk::PlanarPolygon::Pointer figure1 = mitk::PlanarPolygon::New(); figure1->SetPlaneGeometry( m_Geometry ); mitk::Point2D pnt1; pnt1[0] = 10.5 ; pnt1[1] = 3.5; figure1->PlaceFigure( pnt1 ); mitk::Point2D pnt2; pnt2[0] = 9.5; pnt2[1] = 3.5; figure1->SetControlPoint( 1, pnt2, true ); mitk::Point2D pnt3; pnt3[0] = 9.5; pnt3[1] = 4.5; figure1->SetControlPoint( 2, pnt3, true ); mitk::Point2D pnt4; pnt4[0] = 10.5; pnt4[1] = 4.5; figure1->SetControlPoint( 3, pnt4, true ); figure1->GetPolyLine(0); this->VerifyStatistics(ComputeStatistics(m_Image, figure1.GetPointer()), 255.0, 0.0, 255.0); } void mitkImageStatisticsCalculatorTestSuite::TestCase2() { /***************************** * half pixel in x-direction (white) * -> mean of 255 expected ******************************/ MITK_INFO << std::endl << "Test case 2:"; mitk::PlanarPolygon::Pointer figure1 = mitk::PlanarPolygon::New(); figure1->SetPlaneGeometry( m_Geometry ); mitk::Point2D pnt1; pnt1[0] = 10.0 ; pnt1[1] = 3.5; figure1->PlaceFigure( pnt1 ); mitk::Point2D pnt2; pnt2[0] = 9.5; pnt2[1] = 3.5; figure1->SetControlPoint( 1, pnt2, true ); mitk::Point2D pnt3; pnt3[0] = 9.5; pnt3[1] = 4.5; figure1->SetControlPoint( 2, pnt3, true ); mitk::Point2D pnt4; pnt4[0] = 10.0; pnt4[1] = 4.5; figure1->SetControlPoint( 3, pnt4, true ); figure1->GetPolyLine(0); this->VerifyStatistics(ComputeStatistics(m_Image, figure1.GetPointer()), 255.0, 0.0, 255.0); } void mitkImageStatisticsCalculatorTestSuite::TestCase3() { /***************************** * half pixel in diagonal-direction (white) * -> mean of 255 expected ******************************/ MITK_INFO << std::endl << "Test case 3:"; mitk::PlanarPolygon::Pointer figure1 = mitk::PlanarPolygon::New(); figure1->SetPlaneGeometry( m_Geometry ); mitk::Point2D pnt1; pnt1[0] = 10.5 ; pnt1[1] = 3.5; figure1->PlaceFigure( pnt1 ); mitk::Point2D pnt2; pnt2[0] = 9.5; pnt2[1] = 3.5; figure1->SetControlPoint( 1, pnt2, true ); mitk::Point2D pnt3; pnt3[0] = 9.5; pnt3[1] = 4.5; figure1->SetControlPoint( 2, pnt3, true ); figure1->GetPolyLine(0); this->VerifyStatistics(ComputeStatistics(m_Image, figure1.GetPointer()), 255.0, 0.0, 255.0); } void mitkImageStatisticsCalculatorTestSuite::TestCase4() { /***************************** * one pixel (white) + 2 half pixels (white) + 1 half pixel (black) * -> mean of 191.25 expected ******************************/ MITK_INFO << std::endl << "Test case 4:"; mitk::PlanarPolygon::Pointer figure1 = mitk::PlanarPolygon::New(); figure1->SetPlaneGeometry( m_Geometry ); mitk::Point2D pnt1; pnt1[0] = 1.1; pnt1[1] = 1.1; figure1->PlaceFigure( pnt1 ); mitk::Point2D pnt2; pnt2[0] = 2.0; pnt2[1] = 2.0; figure1->SetControlPoint( 1, pnt2, true ); mitk::Point2D pnt3; pnt3[0] = 3.0; pnt3[1] = 1.0; figure1->SetControlPoint( 2, pnt3, true ); mitk::Point2D pnt4; pnt4[0] = 2.0; pnt4[1] = 0.0; figure1->SetControlPoint( 3, pnt4, true ); figure1->GetPolyLine(0); this->VerifyStatistics(ComputeStatistics(m_Image, figure1.GetPointer()), 191.25, 110.41, 242.250); } void mitkImageStatisticsCalculatorTestSuite::TestCase5() { /***************************** * whole pixel (white) + half pixel (gray) in x-direction * -> mean of 191.5 expected ******************************/ MITK_INFO << std::endl << "Test case 5:"; mitk::PlanarPolygon::Pointer figure1 = mitk::PlanarPolygon::New(); figure1->SetPlaneGeometry( m_Geometry ); mitk::Point2D pnt1; pnt1[0] = 11.0; pnt1[1] = 3.5; figure1->PlaceFigure( pnt1 ); mitk::Point2D pnt2; pnt2[0] = 9.5; pnt2[1] = 3.5; figure1->SetControlPoint( 1, pnt2, true ); mitk::Point2D pnt3; pnt3[0] = 9.5; pnt3[1] = 4.5; figure1->SetControlPoint( 2, pnt3, true ); mitk::Point2D pnt4; pnt4[0] = 11.0; pnt4[1] = 4.5; figure1->SetControlPoint( 3, pnt4, true ); figure1->GetPolyLine(0); this->VerifyStatistics(ComputeStatistics(m_Image, figure1.GetPointer()), 191.50, 63.50, 134.340); } void mitkImageStatisticsCalculatorTestSuite::TestCase6() { /***************************** * quarter pixel (black) + whole pixel (white) + half pixel (gray) in x-direction * -> mean of 191.5 expected ******************************/ MITK_INFO << std::endl << "Test case 6:"; mitk::PlanarPolygon::Pointer figure1 = mitk::PlanarPolygon::New(); figure1->SetPlaneGeometry( m_Geometry ); mitk::Point2D pnt1; pnt1[0] = 11.0; pnt1[1] = 3.5; figure1->PlaceFigure( pnt1 ); mitk::Point2D pnt2; pnt2[0] = 9.25; pnt2[1] = 3.5; figure1->SetControlPoint( 1, pnt2, true ); mitk::Point2D pnt3; pnt3[0] = 9.25; pnt3[1] = 4.5; figure1->SetControlPoint( 2, pnt3, true ); mitk::Point2D pnt4; pnt4[0] = 11.0; pnt4[1] = 4.5; figure1->SetControlPoint( 3, pnt4, true ); figure1->GetPolyLine(0); this->VerifyStatistics(ComputeStatistics(m_Image, figure1.GetPointer()), 191.5, 63.50, 134.340); } void mitkImageStatisticsCalculatorTestSuite::TestCase7() { /***************************** * half pixel (black) + whole pixel (white) + half pixel (gray) in x-direction * -> mean of 127.66 expected ******************************/ MITK_INFO << std::endl << "Test case 7:"; mitk::PlanarPolygon::Pointer figure1 = mitk::PlanarPolygon::New(); figure1->SetPlaneGeometry( m_Geometry ); mitk::Point2D pnt1; pnt1[0] = 11.0; pnt1[1] = 3.5; figure1->PlaceFigure( pnt1 ); mitk::Point2D pnt2; pnt2[0] = 9.0; pnt2[1] = 3.5; figure1->SetControlPoint( 1, pnt2, true ); mitk::Point2D pnt3; pnt3[0] = 9.0; pnt3[1] = 4.0; figure1->SetControlPoint( 2, pnt3, true ); mitk::Point2D pnt4; pnt4[0] = 11.0; pnt4[1] = 4.0; figure1->SetControlPoint( 3, pnt4, true ); figure1->GetPolyLine(0); this->VerifyStatistics(ComputeStatistics(m_Image, figure1.GetPointer()), 127.66, 104.1, 140.250); } void mitkImageStatisticsCalculatorTestSuite::TestCase8() { /***************************** * whole pixel (gray) * -> mean of 128 expected ******************************/ MITK_INFO << std::endl << "Test case 8:"; mitk::PlanarPolygon::Pointer figure2 = mitk::PlanarPolygon::New(); figure2->SetPlaneGeometry( m_Geometry ); mitk::Point2D pnt1; pnt1[0] = 11.5; pnt1[1] = 10.5; figure2->PlaceFigure( pnt1 ); mitk::Point2D pnt2; pnt2[0] = 11.5; pnt2[1] = 11.5; figure2->SetControlPoint( 1, pnt2, true ); mitk::Point2D pnt3; pnt3[0] = 12.5; pnt3[1] = 11.5; figure2->SetControlPoint( 2, pnt3, true ); mitk::Point2D pnt4; pnt4[0] = 12.5; pnt4[1] = 10.5; figure2->SetControlPoint( 3, pnt4, true ); figure2->GetPolyLine(0); this->VerifyStatistics(ComputeStatistics(m_Image, figure2.GetPointer()), 128.0, 0.0, 128.0); } void mitkImageStatisticsCalculatorTestSuite::TestCase9() { /***************************** * whole pixel (gray) + half pixel (white) in y-direction * -> mean of 191.5 expected ******************************/ MITK_INFO << std::endl << "Test case 9:"; mitk::PlanarPolygon::Pointer figure2 = mitk::PlanarPolygon::New(); figure2->SetPlaneGeometry( m_Geometry ); mitk::Point2D pnt1; pnt1[0] = 11.5; pnt1[1] = 10.5; figure2->PlaceFigure( pnt1 ); mitk::Point2D pnt2; pnt2[0] = 11.5; pnt2[1] = 12.0; figure2->SetControlPoint( 1, pnt2, true ); mitk::Point2D pnt3; pnt3[0] = 12.5; pnt3[1] = 12.0; figure2->SetControlPoint( 2, pnt3, true ); mitk::Point2D pnt4; pnt4[0] = 12.5; pnt4[1] = 10.5; figure2->SetControlPoint( 3, pnt4, true ); figure2->GetPolyLine(0); this->VerifyStatistics(ComputeStatistics(m_Image, figure2.GetPointer()), 191.5, 63.50, 134.340); } void mitkImageStatisticsCalculatorTestSuite::TestCase10() { /***************************** * 2 whole pixel (white) + 2 whole pixel (black) in y-direction * -> mean of 127.66 expected ******************************/ MITK_INFO << std::endl << "Test case 10:"; mitk::PlanarPolygon::Pointer figure2 = mitk::PlanarPolygon::New(); figure2->SetPlaneGeometry( m_Geometry ); mitk::Point2D pnt1; pnt1[0] = 11.5; pnt1[1] = 10.5; figure2->PlaceFigure( pnt1 ); mitk::Point2D pnt2; pnt2[0] = 11.5; pnt2[1] = 13.5; figure2->SetControlPoint( 1, pnt2, true ); mitk::Point2D pnt3; pnt3[0] = 12.5; pnt3[1] = 13.5; figure2->SetControlPoint( 2, pnt3, true ); mitk::Point2D pnt4; pnt4[0] = 12.5; pnt4[1] = 10.5; figure2->SetControlPoint( 3, pnt4, true ); figure2->GetPolyLine(0); this->VerifyStatistics(ComputeStatistics(m_Image, figure2.GetPointer()), 127.66, 104.1, 140.250); } void mitkImageStatisticsCalculatorTestSuite::TestCase11() { /***************************** * 9 whole pixels (white) + 3 half pixels (white) * + 3 whole pixel (black) [ + 3 slightly less than half pixels (black)] * -> mean of 204.0 expected ******************************/ MITK_INFO << std::endl << "Test case 11:"; mitk::PlanarPolygon::Pointer figure2 = mitk::PlanarPolygon::New(); figure2->SetPlaneGeometry( m_Geometry ); mitk::Point2D pnt1; pnt1[0] = 0.5; pnt1[1] = 0.5; figure2->PlaceFigure( pnt1 ); mitk::Point2D pnt2; pnt2[0] = 3.5; pnt2[1] = 3.5; figure2->SetControlPoint( 1, pnt2, true ); mitk::Point2D pnt3; pnt3[0] = 8.4999; pnt3[1] = 3.5; figure2->SetControlPoint( 2, pnt3, true ); mitk::Point2D pnt4; pnt4[0] = 5.4999; pnt4[1] = 0.5; figure2->SetControlPoint( 3, pnt4, true ); figure2->GetPolyLine(0); this->VerifyStatistics(ComputeStatistics(m_Image, figure2.GetPointer()), 204.0, 102.00, 242.250); } void mitkImageStatisticsCalculatorTestSuite::TestCase12() { /***************************** * half pixel (white) + whole pixel (white) + half pixel (black) * -> mean of 212.66 expected ******************************/ MITK_INFO << std::endl << "Test case 12:"; mitk::PlanarPolygon::Pointer figure2 = mitk::PlanarPolygon::New(); figure2->SetPlaneGeometry( m_Geometry ); mitk::Point2D pnt1; pnt1[0] = 9.5; pnt1[1] = 0.5; figure2->PlaceFigure( pnt1 ); mitk::Point2D pnt2; pnt2[0] = 9.5; pnt2[1] = 2.5; figure2->SetControlPoint( 1, pnt2, true ); mitk::Point2D pnt3; pnt3[0] = 11.5; pnt3[1] = 2.5; figure2->SetControlPoint( 2, pnt3, true ); figure2->GetPolyLine(0); this->VerifyStatistics(ComputeStatistics(m_Image, figure2.GetPointer()), 212.66, 59.860, 248.640); } void mitkImageStatisticsCalculatorTestSuite::TestImageMaskingEmpty() { MITK_INFO << std::endl << "TestImageMaskingEmpty:"; mitk::Image::Pointer mask_image = mitk::ImageGenerator::GenerateImageFromReference( m_Image, 0 ); this->VerifyStatistics( ComputeStatistics( m_Image, mask_image ), -21474836.480, -21474836.480, -21474836.480); // empty statisticsContainer (default values) } void mitkImageStatisticsCalculatorTestSuite::TestImageMaskingNonEmpty() { MITK_INFO << std::endl << "TestImageMaskingNonEmpty:"; mitk::Image::Pointer mask_image = mitk::ImageGenerator::GenerateImageFromReference( m_Image, 0 ); std::vector< itk::Index<3U> > activated_indices; itk::Index<3U> index = {{10, 8, 0}}; activated_indices.push_back( index ); index[0] = 9; index[1] = 8; index[2] = 0; activated_indices.push_back( index ); index[0] = 9; index[1] = 7; index[2] = 0; activated_indices.push_back( index ); index[0] = 10; index[1] = 7; index[2] = 0; activated_indices.push_back( index ); std::vector< itk::Index<3U> >::const_iterator indexIter = activated_indices.begin(); // activate voxel in the mask image mitk::ImagePixelWriteAccessor< unsigned char, 3> writeAccess( mask_image ); while( indexIter != activated_indices.end() ) { writeAccess.SetPixelByIndex( (*indexIter++), 1); } this->VerifyStatistics( ComputeStatistics( m_Image, mask_image ), 127.5, 127.5, 12.750); } void mitkImageStatisticsCalculatorTestSuite::TestRecomputeOnModifiedMask() { MITK_INFO << std::endl << "TestRecomputeOnModifiedMask:"; mitk::Image::Pointer mask_image = mitk::ImageGenerator::GenerateImageFromReference( m_Image, 0 ); mitk::ImageStatisticsCalculator::Pointer statisticsCalculator = mitk::ImageStatisticsCalculator::New(); statisticsCalculator->SetInputImage( m_Image ); mitk::ImageMaskGenerator::Pointer imgMaskGen = mitk::ImageMaskGenerator::New(); imgMaskGen->SetImageMask(mask_image); statisticsCalculator->SetMask(imgMaskGen.GetPointer()); this->VerifyStatistics( statisticsCalculator->GetStatistics(), -21474836.480, -21474836.480, -21474836.480); // activate voxel in the mask image itk::Index<3U> test_index = {11, 8, 0}; mitk::ImagePixelWriteAccessor< unsigned char, 3> writeAccess( mask_image ); writeAccess.SetPixelByIndex( test_index, 1); mask_image->Modified(); mitk::ImageStatisticsCalculator::StatisticsContainer::Pointer stat = statisticsCalculator->GetStatistics(); this->VerifyStatistics( stat, 128.0, 0.0, 128.0); MITK_TEST_CONDITION( stat->GetN() == 1, "Calculated mask voxel count '" << stat->GetN() << "' is equal to the desired value '" << 1 << "'" ); } const mitk::ImageStatisticsCalculator::StatisticsContainer::Pointer mitkImageStatisticsCalculatorTestSuite::ComputeStatistics( mitk::Image::Pointer image, mitk::PlanarFigure::Pointer polygon ) { mitk::ImageStatisticsCalculator::Pointer statisticsCalculator = mitk::ImageStatisticsCalculator::New(); statisticsCalculator->SetInputImage( image ); statisticsCalculator->SetNBinsForHistogramStatistics(10); mitk::PlanarFigureMaskGenerator::Pointer planFigMaskGen = mitk::PlanarFigureMaskGenerator::New(); - planFigMaskGen->SetImage(image); + planFigMaskGen->SetInputImage(image); planFigMaskGen->SetPlanarFigure(polygon); statisticsCalculator->SetMask(planFigMaskGen.GetPointer()); try { return statisticsCalculator->GetStatistics(); } catch( ... ) { } return mitk::ImageStatisticsCalculator::StatisticsContainer::New(); } const mitk::ImageStatisticsCalculator::StatisticsContainer::Pointer mitkImageStatisticsCalculatorTestSuite::ComputeStatistics(mitk::Image::Pointer image, mitk::Image::Pointer image_mask ) { mitk::ImageStatisticsCalculator::Pointer statisticsCalculator = mitk::ImageStatisticsCalculator::New(); statisticsCalculator->SetInputImage(image); statisticsCalculator->SetNBinsForHistogramStatistics(10); mitk::ImageMaskGenerator::Pointer imgMaskGen = mitk::ImageMaskGenerator::New(); imgMaskGen->SetImageMask(image_mask); statisticsCalculator->SetMask(imgMaskGen.GetPointer()); return statisticsCalculator->GetStatistics(); } void mitkImageStatisticsCalculatorTestSuite::VerifyStatistics(mitk::ImageStatisticsCalculator::StatisticsContainer::Pointer stats, double testMean, double testSD, double testMedian) { int tmpMean = stats->GetMean() * 100; double calculatedMean = tmpMean / 100.0; MITK_TEST_CONDITION( calculatedMean == testMean, "Calculated mean grayvalue '" << calculatedMean << "' is equal to the desired value '" << testMean << "'" ); int tmpSD = stats->GetStd() * 100; double calculatedSD = tmpSD / 100.0; MITK_TEST_CONDITION( calculatedSD == testSD, "Calculated grayvalue sd '" << calculatedSD << "' is equal to the desired value '" << testSD <<"'" ); int tmpMedian = stats->GetMedian() * 100; double calculatedMedian = tmpMedian / 100.0; MITK_TEST_CONDITION( testMedian == calculatedMedian, "Calculated median grayvalue '" << calculatedMedian << "' is equal to the desired value '" << testMedian << "'"); } void mitkImageStatisticsCalculatorTestSuite::TestUninitializedImage() { /***************************** * loading uninitialized image to datastorage ******************************/ MITK_TEST_FOR_EXCEPTION_BEGIN(mitk::Exception) mitk::Image::Pointer image = mitk::Image::New(); mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData(image); mitk::ImageStatisticsCalculator::Pointer is = mitk::ImageStatisticsCalculator::New(); is->GetStatistics(); MITK_TEST_FOR_EXCEPTION_END(mitk::Exception) } MITK_TEST_SUITE_REGISTRATION(mitkImageStatisticsCalculator) diff --git a/Modules/ImageStatistics/Testing/mitkImageStatisticsHotspotTest.cpp b/Modules/ImageStatistics/Testing/mitkImageStatisticsHotspotTest.cpp index 0a3d03b00e..5372b3b769 100644 --- a/Modules/ImageStatistics/Testing/mitkImageStatisticsHotspotTest.cpp +++ b/Modules/ImageStatistics/Testing/mitkImageStatisticsHotspotTest.cpp @@ -1,649 +1,649 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkImageStatisticsCalculator.h" #include "itkMultiGaussianImageSource.h" #include "mitkTestingMacros.h" #include "mitkImageCast.h" #include #include #include #include #include #include /** \section hotspotCalculationTestCases Testcases To see the different Hotspot-Testcases have a look at the \ref hotspottestdoc. Note from an intensive session of checking the test results: - itk::MultiGaussianImageSource needs a review - the test idea is ok, but the combination of XML files for parameters and MultiGaussianImageSource has serious flaws - the XML file should contain exactly the parameters that MultiGaussianImageSource requires - in contrast, now the XML file mentions index coordinates for gaussian centers while the MultiGaussianImageSource expects world coordinates - this requires a transformation (index * spacing assuming no rotation) that was actually broken until recently */ struct mitkImageStatisticsHotspotTestClass { /** \brief Test parameters for one test case. Describes all aspects of a single test case: - parameters to generate a test image - parameters of a ROI that describes where to calculate statistics - expected statistics results */ struct Parameters { public: // XML-Tag /** \brief XML-Tag "image-rows": size of x-dimension */ int m_ImageRows; /** \brief XML-Tag "image-columns": size of y-dimension */ int m_ImageColumns; /** \brief XML-Tag "image-slices": size of z-dimension */ int m_ImageSlices; /** \brief XML-Tag "numberOfGaussians": number of used gauss-functions */ int m_NumberOfGaussian; /** \brief XML-Tags "spacingX", "spacingY", "spacingZ": spacing of image in every direction */ double m_Spacing[3]; /** \brief XML-Tag "entireHotSpotInImage" */ unsigned int m_EntireHotspotInImage; // XML-Tag /** \brief XML-Tag "centerIndexX: gaussian parameter \warning This parameter READS the centerIndexX parameter from file and is THEN MISUSED to calculate some position in world coordinates, so we require double. */ std::vector m_CenterX; /** \brief XML-Tag "centerIndexY: gaussian parameter \warning This parameter READS the centerIndexX parameter from file and is THEN MISUSED to calculate some position in world coordinates, so we require double. */ std::vector m_CenterY; /** \brief XML-Tag "centerIndexZ: gaussian parameter \warning This parameter READS the centerIndexX parameter from file and is THEN MISUSED to calculate some position in world coordinates, so we require double. */ std::vector m_CenterZ; /** \brief XML-Tag "deviationX: gaussian parameter */ std::vector m_SigmaX; /** \brief XML-Tag "deviationY: gaussian parameter */ std::vector m_SigmaY; /** \brief XML-Tag "deviationZ: gaussian parameter */ std::vector m_SigmaZ; /** \brief XML-Tag "altitude: gaussian parameter */ std::vector m_Altitude; // XML-Tag /** \brief XML-Tag "numberOfLabels": number of different labels which appear in the mask */ unsigned int m_NumberOfLabels; /** \brief XML-Tag "hotspotRadiusInMM": radius of hotspot */ double m_HotspotRadiusInMM; // XML-Tag /** \brief XML-Tag "maximumSizeX": maximum position of ROI in x-dimension */ vnl_vector m_MaxIndexX; /** \brief XML-Tag "minimumSizeX": minimum position of ROI in x-dimension */ vnl_vector m_MinIndexX; /** \brief XML-Tag "maximumSizeX": maximum position of ROI in y-dimension */ vnl_vector m_MaxIndexY; /** \brief XML-Tag "minimumSizeX": minimum position of ROI in y-dimension */ vnl_vector m_MinIndexY; /** \brief XML-Tag "maximumSizeX": maximum position of ROI in z-dimension */ vnl_vector m_MaxIndexZ; /** \brief XML-Tag "minimumSizeX": minimum position of ROI in z-dimension */ vnl_vector m_MinIndexZ; /** \brief XML-Tag "label": value of label */ vnl_vector m_Label; //XML-Tag /** \brief XML-Tag "minimum": minimum inside hotspot */ vnl_vector m_HotspotMin; /** \brief XML-Tag "maximum": maximum inside hotspot */ vnl_vector m_HotspotMax; /** \brief XML-Tag "mean": mean value of hotspot */ vnl_vector m_HotspotMean; /** \brief XML-Tag "maximumIndexX": x-coordinate of maximum-location inside hotspot */ vnl_vector m_HotspotMaxIndexX; /** \brief XML-Tag "maximumIndexX": y-coordinate of maximum-location inside hotspot */ vnl_vector m_HotspotMaxIndexY; /** \brief XML-Tag "maximumIndexX": z-coordinate of maximum-location inside hotspot */ vnl_vector m_HotspotMaxIndexZ; /** \brief XML-Tag "maximumIndexX": x-coordinate of maximum-location inside hotspot */ vnl_vector m_HotspotMinIndexX; /** \brief XML-Tag "maximumIndexX": y-coordinate of maximum-location inside hotspot */ vnl_vector m_HotspotMinIndexY; /** \brief XML-Tag "maximumIndexX": z-coordinate of maximum-location inside hotspot */ vnl_vector m_HotspotMinIndexZ; /** \brief XML-Tag "maximumIndexX": x-coordinate of hotspot-location */ vnl_vector m_HotspotIndexX; /** \brief XML-Tag "maximumIndexX": y-coordinate of hotspot-location */ vnl_vector m_HotspotIndexY; /** \brief XML-Tag "maximumIndexX": z-coordinate of hotspot-location */ vnl_vector m_HotspotIndexZ; }; /** \brief Find/Convert integer attribute in itk::DOMNode. */ static int GetIntegerAttribute(itk::DOMNode* domNode, const std::string& tag) { assert(domNode); MITK_TEST_CONDITION_REQUIRED( domNode->HasAttribute(tag), "Tag '" << tag << "' is defined in test parameters" ); std::string attributeValue = domNode->GetAttribute(tag); int resultValue; try { //MITK_TEST_OUTPUT( << "Converting tag value '" << attributeValue << "' for tag '" << tag << "' to integer"); std::stringstream(attributeValue) >> resultValue; return resultValue; } catch(std::exception& /*e*/) { MITK_TEST_CONDITION_REQUIRED(false, "Convert tag value '" << attributeValue << "' for tag '" << tag << "' to integer"); return 0; // just to satisfy compiler } } /** \brief Find/Convert double attribute in itk::DOMNode. */ static double GetDoubleAttribute(itk::DOMNode* domNode, const std::string& tag) { assert(domNode); MITK_TEST_CONDITION_REQUIRED( domNode->HasAttribute(tag), "Tag '" << tag << "' is defined in test parameters" ); std::string attributeValue = domNode->GetAttribute(tag); double resultValue; try { //MITK_TEST_OUTPUT( << "Converting tag value '" << attributeValue << "' for tag '" << tag << "' to double"); std::stringstream(attributeValue) >> resultValue; return resultValue; } catch(std::exception& /*e*/) { MITK_TEST_CONDITION_REQUIRED(false, "Convert tag value '" << attributeValue << "' for tag '" << tag << "' to double"); return 0.0; // just to satisfy compiler } } /** \brief Read XML file describing the test parameters. Reads XML file given in first commandline parameter in order to construct a Parameters structure. The XML file should be structurs as the following example, i.e. we describe the three test aspects of Parameters in four different tags, with all the details described as tag attributes. */ /** \verbatim \endverbatim */ static Parameters ParseParameters(int argc, char* argv[]) { MITK_TEST_CONDITION_REQUIRED(argc == 2, "Test is invoked with exactly 1 parameter (XML parameters file)"); MITK_INFO << "Reading parameters from file '" << argv[1] << "'"; std::string filename = argv[1]; Parameters result; itk::DOMNodeXMLReader::Pointer xmlReader = itk::DOMNodeXMLReader::New(); xmlReader->SetFileName( filename ); try { xmlReader->Update(); itk::DOMNode::Pointer domRoot = xmlReader->GetOutput(); typedef std::vector NodeList; NodeList testimages; domRoot->GetChildren("testimage", testimages); MITK_TEST_CONDITION_REQUIRED( testimages.size() == 1, "One test image defined" ) itk::DOMNode* testimage = testimages[0]; result.m_ImageRows = GetIntegerAttribute( testimage, "image-rows" ); result.m_ImageColumns = GetIntegerAttribute( testimage, "image-columns" ); result.m_ImageSlices = GetIntegerAttribute( testimage, "image-slices" ); result.m_NumberOfGaussian = GetIntegerAttribute( testimage, "numberOfGaussians" ); result.m_Spacing[0] = GetDoubleAttribute(testimage, "spacingX"); result.m_Spacing[1] = GetDoubleAttribute(testimage, "spacingY"); result.m_Spacing[2] = GetDoubleAttribute(testimage, "spacingZ"); result.m_EntireHotspotInImage = GetIntegerAttribute( testimage, "entireHotSpotInImage" ); MITK_TEST_OUTPUT( << "Read size parameters (x,y,z): " << result.m_ImageRows << "," << result.m_ImageColumns << "," << result.m_ImageSlices); MITK_TEST_OUTPUT( << "Read spacing parameters (x,y,z): " << result.m_Spacing[0] << "," << result.m_Spacing[1] << "," << result.m_Spacing[2]); NodeList gaussians; testimage->GetChildren("gaussian", gaussians); MITK_TEST_CONDITION_REQUIRED( gaussians.size() >= 1, "At least one gaussian is defined" ) result.m_CenterX.resize(result.m_NumberOfGaussian); result.m_CenterY.resize(result.m_NumberOfGaussian); result.m_CenterZ.resize(result.m_NumberOfGaussian); result.m_SigmaX.resize(result.m_NumberOfGaussian); result.m_SigmaY.resize(result.m_NumberOfGaussian); result.m_SigmaZ.resize(result.m_NumberOfGaussian); result.m_Altitude.resize(result.m_NumberOfGaussian); for(int i = 0; i < result.m_NumberOfGaussian ; ++i) { itk::DOMNode* gaussian = gaussians[i]; result.m_CenterX[i] = GetIntegerAttribute(gaussian, "centerIndexX"); result.m_CenterY[i] = GetIntegerAttribute(gaussian, "centerIndexY"); result.m_CenterZ[i] = GetIntegerAttribute(gaussian, "centerIndexZ"); result.m_SigmaX[i] = GetDoubleAttribute(gaussian, "deviationX"); result.m_SigmaY[i] = GetDoubleAttribute(gaussian, "deviationY"); result.m_SigmaZ[i] = GetDoubleAttribute(gaussian, "deviationZ"); result.m_Altitude[i] = GetDoubleAttribute(gaussian, "altitude"); result.m_CenterX[i] = result.m_CenterX[i] * result.m_Spacing[0]; result.m_CenterY[i] = result.m_CenterY[i] * result.m_Spacing[1]; result.m_CenterZ[i] = result.m_CenterZ[i] * result.m_Spacing[2]; result.m_SigmaX[i] = result.m_SigmaX[i] * result.m_Spacing[0]; result.m_SigmaY[i] = result.m_SigmaY[i] * result.m_Spacing[1]; result.m_SigmaZ[i] = result.m_SigmaZ[i] * result.m_Spacing[2]; } NodeList segmentations; domRoot->GetChildren("segmentation", segmentations); MITK_TEST_CONDITION_REQUIRED( segmentations.size() == 1, "One segmentation defined"); itk::DOMNode* segmentation = segmentations[0]; result.m_NumberOfLabels = GetIntegerAttribute(segmentation, "numberOfLabels"); result.m_HotspotRadiusInMM = GetDoubleAttribute(segmentation, "hotspotRadiusInMM"); // read ROI parameters, fill result structure NodeList rois; segmentation->GetChildren("roi", rois); MITK_TEST_CONDITION_REQUIRED( rois.size() >= 1, "At least one ROI defined" ) result.m_MaxIndexX.set_size(result.m_NumberOfLabels); result.m_MinIndexX.set_size(result.m_NumberOfLabels); result.m_MaxIndexY.set_size(result.m_NumberOfLabels); result.m_MinIndexY.set_size(result.m_NumberOfLabels); result.m_MaxIndexZ.set_size(result.m_NumberOfLabels); result.m_MinIndexZ.set_size(result.m_NumberOfLabels); result.m_Label.set_size(result.m_NumberOfLabels); for(unsigned int i = 0; i < rois.size(); ++i) { result.m_MaxIndexX[i] = GetIntegerAttribute(rois[i], "maximumIndexX"); result.m_MinIndexX[i] = GetIntegerAttribute(rois[i], "minimumIndexX"); result.m_MaxIndexY[i] = GetIntegerAttribute(rois[i], "maximumIndexY"); result.m_MinIndexY[i] = GetIntegerAttribute(rois[i], "minimumIndexY"); result.m_MaxIndexZ[i] = GetIntegerAttribute(rois[i], "maximumIndexZ"); result.m_MinIndexZ[i] = GetIntegerAttribute(rois[i], "minimumIndexZ"); result.m_Label[i] = GetIntegerAttribute(rois[i], "label"); } // read statistic parameters, fill result structure NodeList statistics; domRoot->GetChildren("statistic", statistics); MITK_TEST_CONDITION_REQUIRED( statistics.size() >= 1 , "At least one statistic defined" ) MITK_TEST_CONDITION_REQUIRED( statistics.size() == rois.size(), "Same number of rois and corresponding statistics defined"); result.m_HotspotMin.set_size(statistics.size()); result.m_HotspotMax.set_size(statistics.size()); result.m_HotspotMean.set_size(statistics.size()); result.m_HotspotMinIndexX.set_size(statistics.size()); result.m_HotspotMinIndexY.set_size(statistics.size()); result.m_HotspotMinIndexZ.set_size(statistics.size()); result.m_HotspotMaxIndexX.set_size(statistics.size()); result.m_HotspotMaxIndexY.set_size(statistics.size()); result.m_HotspotMaxIndexZ.set_size(statistics.size()); result.m_HotspotIndexX.set_size(statistics.size()); result.m_HotspotIndexY.set_size(statistics.size()); result.m_HotspotIndexZ.set_size(statistics.size()); for(unsigned int i = 0; i < statistics.size(); ++i) { result.m_HotspotMin[i] = GetDoubleAttribute(statistics[i], "minimum"); result.m_HotspotMax[i] = GetDoubleAttribute(statistics[i], "maximum"); result.m_HotspotMean[i] = GetDoubleAttribute(statistics[i], "mean"); result.m_HotspotMinIndexX[i] = GetIntegerAttribute(statistics[i], "minimumIndexX"); result.m_HotspotMinIndexY[i] = GetIntegerAttribute(statistics[i], "minimumIndexY"); result.m_HotspotMinIndexZ[i] = GetIntegerAttribute(statistics[i], "minimumIndexZ"); result.m_HotspotMaxIndexX[i] = GetIntegerAttribute(statistics[i], "maximumIndexX"); result.m_HotspotMaxIndexY[i] = GetIntegerAttribute(statistics[i], "maximumIndexY"); result.m_HotspotMaxIndexZ[i] = GetIntegerAttribute(statistics[i], "maximumIndexZ"); result.m_HotspotIndexX[i] = GetIntegerAttribute(statistics[i], "hotspotIndexX"); result.m_HotspotIndexY[i] = GetIntegerAttribute(statistics[i], "hotspotIndexY"); result.m_HotspotIndexZ[i] = GetIntegerAttribute(statistics[i], "hotspotIndexZ"); } } catch (std::exception& e) { MITK_TEST_CONDITION_REQUIRED(false, "Reading test parameters from XML file. Error message: " << e.what()); } return result; } /** \brief Generate an image that contains a couple of 3D gaussian distributions. Uses the given parameters to produce a test image using class MultiGaussianImageSource. */ static mitk::Image::Pointer BuildTestImage(const Parameters& testParameters) { mitk::Image::Pointer result; typedef double PixelType; const int Dimension = 3; typedef itk::Image ImageType; ImageType::Pointer image = ImageType::New(); typedef itk::MultiGaussianImageSource< ImageType > MultiGaussianImageSource; MultiGaussianImageSource::Pointer gaussianGenerator = MultiGaussianImageSource::New(); ImageType::SizeValueType size[3]; size[0] = testParameters.m_ImageColumns; size[1] = testParameters.m_ImageRows; size[2] = testParameters.m_ImageSlices; itk::MultiGaussianImageSource::VectorType centerXVec, centerYVec, centerZVec, sigmaXVec, sigmaYVec, sigmaZVec, altitudeVec; for(int i = 0; i < testParameters.m_NumberOfGaussian; ++i) { centerXVec.push_back(testParameters.m_CenterX[i]); centerYVec.push_back(testParameters.m_CenterY[i]); centerZVec.push_back(testParameters.m_CenterZ[i]); sigmaXVec.push_back(testParameters.m_SigmaX[i]); sigmaYVec.push_back(testParameters.m_SigmaY[i]); sigmaZVec.push_back(testParameters.m_SigmaZ[i]); altitudeVec.push_back(testParameters.m_Altitude[i]); } ImageType::SpacingType spacing; for( int i = 0; i < Dimension; ++i ) spacing[i] = testParameters.m_Spacing[i]; gaussianGenerator->SetSize( size ); gaussianGenerator->SetSpacing( spacing ); gaussianGenerator->SetRadius(testParameters.m_HotspotRadiusInMM); gaussianGenerator->SetNumberOfGausssians(testParameters.m_NumberOfGaussian); gaussianGenerator->AddGaussian(centerXVec, centerYVec, centerZVec, sigmaXVec, sigmaYVec, sigmaZVec, altitudeVec); gaussianGenerator->Update(); image = gaussianGenerator->GetOutput(); mitk::CastToMitkImage(image, result); return result; } /** \brief Calculates hotspot statistics for given test image and ROI parameters. Uses ImageStatisticsCalculator to find a hotspot in a defined ROI within the given image. */ static mitk::ImageStatisticsCalculator::StatisticsContainer::Pointer CalculateStatistics(mitk::Image* image, const Parameters& testParameters, unsigned int label) { mitk::ImageStatisticsCalculator::StatisticsContainer::Pointer result; const unsigned int Dimension = 3; typedef itk::Image MaskImageType; MaskImageType::Pointer mask = MaskImageType::New(); MaskImageType::SizeType size; MaskImageType::SpacingType spacing; MaskImageType::IndexType start; mitk::ImageStatisticsCalculator::Pointer statisticsCalculator = mitk::ImageStatisticsCalculator::New(); statisticsCalculator->SetInputImage(image); mitk::Image::Pointer mitkMaskImage; if((testParameters.m_MaxIndexX[label] > testParameters.m_MinIndexX[label] && testParameters.m_MinIndexX[label] >= 0) && (testParameters.m_MaxIndexY[label] > testParameters.m_MinIndexY[label] && testParameters.m_MinIndexY[label] >= 0) && (testParameters.m_MaxIndexZ[label] > testParameters.m_MinIndexZ[label] && testParameters.m_MinIndexZ[label] >= 0)) { for(unsigned int i = 0; i < Dimension; ++i) { start[i] = 0; spacing[i] = testParameters.m_Spacing[i]; } size[0] = testParameters.m_ImageColumns; size[1] = testParameters.m_ImageRows; size[2] = testParameters.m_ImageSlices; MaskImageType::RegionType region; region.SetIndex(start); region.SetSize(size); mask->SetSpacing(spacing); mask->SetRegions(region); mask->Allocate(); typedef itk::ImageRegionIteratorWithIndex MaskImageIteratorType; MaskImageIteratorType maskIt(mask, region); for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt) { maskIt.Set(0); } for(unsigned int i = 0; i < testParameters.m_NumberOfLabels; ++i) { for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt) { MaskImageType::IndexType index = maskIt.GetIndex(); if((index[0] >= testParameters.m_MinIndexX[i] && index[0] <= testParameters.m_MaxIndexX[i] ) && (index[1] >= testParameters.m_MinIndexY[i] && index[1] <= testParameters.m_MaxIndexY[i] ) && (index[2] >= testParameters.m_MinIndexZ[i] && index[2] <= testParameters.m_MaxIndexZ[i] )) { maskIt.Set(testParameters.m_Label[i]); } } } mitk::CastToMitkImage(mask, mitkMaskImage); mitk::ImageMaskGenerator::Pointer imgMaskGen = mitk::ImageMaskGenerator::New(); imgMaskGen->SetImageMask(mitkMaskImage); mitk::HotspotMaskGenerator::Pointer hotspotMaskGen = mitk::HotspotMaskGenerator::New(); - hotspotMaskGen->SetImage(image); + hotspotMaskGen->SetInputImage(image); hotspotMaskGen->SetMask(imgMaskGen.GetPointer()); hotspotMaskGen->SetHotspotRadiusInMM(testParameters.m_HotspotRadiusInMM); if(testParameters.m_EntireHotspotInImage == 1) { MITK_INFO << "Hotspot must be completly inside image"; hotspotMaskGen->SetHotspotMustBeCompletelyInsideImage(true); } else { MITK_INFO << "Hotspot must not be completly inside image"; hotspotMaskGen->SetHotspotMustBeCompletelyInsideImage(false); } statisticsCalculator->SetMask(hotspotMaskGen.GetPointer()); MITK_DEBUG << "Masking is set to hotspot+image mask"; } else { mitk::HotspotMaskGenerator::Pointer hotspotMaskGen = mitk::HotspotMaskGenerator::New(); - hotspotMaskGen->SetImage(image); + hotspotMaskGen->SetInputImage(image); hotspotMaskGen->SetHotspotRadiusInMM(testParameters.m_HotspotRadiusInMM); if(testParameters.m_EntireHotspotInImage == 1) { MITK_INFO << "Hotspot must be completly inside image"; hotspotMaskGen->SetHotspotMustBeCompletelyInsideImage(true); } else { MITK_INFO << "Hotspot must not be completly inside image"; hotspotMaskGen->SetHotspotMustBeCompletelyInsideImage(false); } MITK_DEBUG << "Masking is set to hotspot only"; } result = statisticsCalculator->GetStatistics(0, label); return result; } static void ValidateStatisticsItem(const std::string& label, double testvalue, double reference, double tolerance) { double diff = ::fabs(reference - testvalue); MITK_TEST_CONDITION( diff < tolerance, "'" << label << "' value close enough to reference value " "(value=" << testvalue << ", reference=" << reference << ", diff=" << diff << ")" ); } static void ValidateStatisticsItem(const std::string& label, const vnl_vector& testvalue, const vnl_vector& reference) { double diffX = ::fabs(double(testvalue[0] - reference[0])); double diffY = ::fabs(double(testvalue[1] - reference[1])); double diffZ = ::fabs(double(testvalue[2] - reference[2])); std::stringstream testPosition; testPosition << testvalue[0] << "," << testvalue[1] << "," << testvalue[2]; std::stringstream referencePosition; referencePosition << reference[0] << "," << reference[1] << "," << reference[2]; MITK_TEST_CONDITION( diffX < mitk::eps && diffY < mitk::eps && diffZ < mitk::eps, "'" << label << "' close enough to reference value " << "(value=[" << testPosition.str() << "]," << " reference=[" << referencePosition.str() << "]"); } /** \brief Compares calculated against actual statistics values. Checks validness of all statistics aspects. Lets test fail if any aspect is not sufficiently equal. */ static void ValidateStatistics(const mitk::ImageStatisticsCalculator::StatisticsContainer::Pointer hotspotStatistics, const Parameters& testParameters, unsigned int label) { // check all expected test result against actual results double eps = 0.25; // value above the largest tested difference ValidateStatisticsItem("Hotspot mean", hotspotStatistics->GetMean(), testParameters.m_HotspotMean[label], eps); ValidateStatisticsItem("Hotspot maximum", hotspotStatistics->GetMax(), testParameters.m_HotspotMax[label], eps); ValidateStatisticsItem("Hotspot minimum", hotspotStatistics->GetMin(), testParameters.m_HotspotMin[label], eps); vnl_vector referenceHotspotCenterIndex; referenceHotspotCenterIndex.set_size(3); referenceHotspotCenterIndex[0] = testParameters.m_HotspotIndexX[label]; referenceHotspotCenterIndex[1] = testParameters.m_HotspotIndexY[label]; referenceHotspotCenterIndex[2] = testParameters.m_HotspotIndexZ[label]; // ValidateStatisticsItem("Hotspot center position", statistics.GetHotspotStatistics().GetHotspotIndex(), referenceHotspotCenterIndex); TODO: new image statistics calculator does not give hotspot position // TODO we do not test minimum/maximum positions within the peak/hotspot region, because // these positions are not unique, i.e. there are multiple valid minima/maxima positions. // One solution would be to modify the test cases in order to achive clear positions. // The BETTER/CORRECT solution would be to change the singular position into a set of positions / a region } }; /** \brief Verifies that hotspot statistics part of ImageStatisticsCalculator. The test reads parameters from an XML-file to generate a test-image, calculates the hotspot statistics of the image and checks if the calculated statistics are the same as the specified values of the XML-file. */ int mitkImageStatisticsHotspotTest(int argc, char* argv[]) { MITK_TEST_BEGIN("mitkImageStatisticsHotspotTest") try { mitkImageStatisticsHotspotTestClass::Parameters parameters = mitkImageStatisticsHotspotTestClass::ParseParameters(argc,argv); mitk::Image::Pointer image = mitkImageStatisticsHotspotTestClass::BuildTestImage(parameters); MITK_TEST_CONDITION_REQUIRED( image.IsNotNull(), "Generate test image" ); for(unsigned int label = 0; label < parameters.m_NumberOfLabels; ++label) { mitk::ImageStatisticsCalculator::StatisticsContainer::Pointer statistics = mitkImageStatisticsHotspotTestClass::CalculateStatistics(image, parameters, label); mitkImageStatisticsHotspotTestClass::ValidateStatistics(statistics, parameters, label); std::cout << std::endl; } } catch (std::exception& e) { std::cout << "Error: " << e.what() << std::endl; MITK_TEST_CONDITION_REQUIRED( false, "Exception occurred during test execution: " << e.what() ); } catch(...) { MITK_TEST_CONDITION_REQUIRED( false, "Exception occurred during test execution." ); } MITK_TEST_END() } diff --git a/Modules/ImageStatistics/mitkHotspotMaskGenerator.cpp b/Modules/ImageStatistics/mitkHotspotMaskGenerator.cpp index a641fe15f1..da07560989 100644 --- a/Modules/ImageStatistics/mitkHotspotMaskGenerator.cpp +++ b/Modules/ImageStatistics/mitkHotspotMaskGenerator.cpp @@ -1,595 +1,595 @@ #include #include #include #include #include #include "mitkImageAccessByItk.h" #include #include #include namespace mitk { HotspotMaskGenerator::HotspotMaskGenerator(): m_HotspotRadiusinMM(6.2035049089940), // radius of a 1cm3 sphere in mm m_HotspotMustBeCompletelyInsideImage(true), m_Label(1) { m_TimeStep = 0; m_InternalMask = mitk::Image::New(); m_InternalMaskUpdateTime = 0; } - void HotspotMaskGenerator::SetImage(mitk::Image::Pointer inputImage) + void HotspotMaskGenerator::SetInputImage(mitk::Image::Pointer inputImage) { - if (inputImage != m_Image) + if (inputImage != m_inputImage) { - m_Image = inputImage; + m_inputImage = inputImage; m_ConvolutionImageMaxIndex.set_size(inputImage->GetDimension()); m_ConvolutionImageMinIndex.set_size(inputImage->GetDimension()); this->Modified(); } } void HotspotMaskGenerator::SetMask(MaskGenerator::Pointer mask) { if (mask != m_Mask) { m_Mask = mask; this->Modified(); } } HotspotMaskGenerator::~HotspotMaskGenerator() { } void HotspotMaskGenerator::SetHotspotRadiusInMM(double radiusInMillimeter) { if(radiusInMillimeter != m_HotspotRadiusinMM) { m_HotspotRadiusinMM = radiusInMillimeter; this->Modified(); } } const double& HotspotMaskGenerator::GetHotspotRadiusinMM() const { return m_HotspotRadiusinMM; } bool HotspotMaskGenerator::GetHotspotMustBeCompletelyInsideImage() const { return m_HotspotMustBeCompletelyInsideImage; } void HotspotMaskGenerator::SetHotspotMustBeCompletelyInsideImage(bool mustBeCompletelyInImage) { if (m_HotspotMustBeCompletelyInsideImage != mustBeCompletelyInImage) { m_HotspotMustBeCompletelyInsideImage = mustBeCompletelyInImage; this->Modified(); } } mitk::Image::Pointer HotspotMaskGenerator::GetMask() { if (IsUpdateRequired()) { - if ( m_Image.IsNull() ) + if ( m_inputImage.IsNull() ) { throw std::runtime_error( "Error: image empty!" ); } - if ( m_TimeStep >= m_Image->GetTimeSteps() ) + if ( m_TimeStep >= m_inputImage->GetTimeSteps() ) { throw std::runtime_error( "Error: invalid time step!" ); } mitk::ImageTimeSelector::Pointer imageTimeSelector = mitk::ImageTimeSelector::New(); - imageTimeSelector->SetInput( m_Image ); + imageTimeSelector->SetInput( m_inputImage ); imageTimeSelector->SetTimeNr( m_TimeStep ); imageTimeSelector->UpdateLargestPossibleRegion(); mitk::Image::Pointer timeSliceImage = imageTimeSelector->GetOutput(); m_internalImage = timeSliceImage; m_internalMask2D = nullptr; // is this correct when this variable holds a smart pointer? m_internalMask3D = nullptr; if ( m_Mask != nullptr ) { m_Mask->SetTimeStep(m_TimeStep); mitk::Image::Pointer timeSliceMask = m_Mask->GetMask(); if ( m_internalImage->GetDimension() == 3 ) { CastToItkImage(timeSliceMask, m_internalMask3D); AccessFixedDimensionByItk_2(m_internalImage, CalculateHotspotMask, 3, m_internalMask3D, m_Label); } else if ( m_internalImage->GetDimension() == 2 ) { CastToItkImage(timeSliceMask, m_internalMask2D); AccessFixedDimensionByItk_2(m_internalImage, CalculateHotspotMask, 2, m_internalMask2D, m_Label); } else { throw std::runtime_error( "Error: invalid image dimension" ); } } else { if ( m_internalImage->GetDimension() == 3 ) { AccessFixedDimensionByItk_2(m_internalImage, CalculateHotspotMask, 3, m_internalMask3D, m_Label); } else if ( m_internalImage->GetDimension() == 2 ) { AccessFixedDimensionByItk_2(m_internalImage, CalculateHotspotMask, 2, m_internalMask2D, m_Label); } else { throw std::runtime_error( "Error: invalid image dimension" ); } } this->Modified(); } m_InternalMaskUpdateTime = m_InternalMask->GetMTime(); return m_InternalMask; } void HotspotMaskGenerator::SetLabel(unsigned short label) { if (label != m_Label) { m_Label = label; this->Modified(); } } vnl_vector HotspotMaskGenerator::GetConvolutionImageMinIndex() { this->GetMask(); // make sure we are up to date return m_ConvolutionImageMinIndex; } vnl_vector HotspotMaskGenerator::GetHotspotIndex() { this->GetMask(); // make sure we are up to date return m_ConvolutionImageMaxIndex; } template HotspotMaskGenerator::ImageExtrema HotspotMaskGenerator::CalculateExtremaWorld( const itk::Image* inputImage, typename itk::Image::Pointer maskImage, double neccessaryDistanceToImageBorderInMM, unsigned int label ) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< unsigned short, VImageDimension > MaskImageType; typedef itk::ImageRegionConstIteratorWithIndex MaskImageIteratorType; typedef itk::ImageRegionConstIteratorWithIndex InputImageIndexIteratorType; typename ImageType::SpacingType spacing = inputImage->GetSpacing(); ImageExtrema minMax; minMax.Defined = false; minMax.MaxIndex.set_size(VImageDimension); minMax.MaxIndex.set_size(VImageDimension); typename ImageType::RegionType allowedExtremaRegion = inputImage->GetLargestPossibleRegion(); bool keepDistanceToImageBorders( neccessaryDistanceToImageBorderInMM > 0 ); if (keepDistanceToImageBorders) { long distanceInPixels[VImageDimension]; for(unsigned short dimension = 0; dimension < VImageDimension; ++dimension) { // To confirm that the whole hotspot is inside the image we have to keep a specific distance to the image-borders, which is as long as // the radius. To get the amount of indices we divide the radius by spacing and add 0.5 because voxels are center based: // For example with a radius of 2.2 and a spacing of 1 two indices are enough because 2.2 / 1 + 0.5 = 2.7 => 2. // But with a radius of 2.7 we need 3 indices because 2.7 / 1 + 0.5 = 3.2 => 3 distanceInPixels[dimension] = int( neccessaryDistanceToImageBorderInMM / spacing[dimension] + 0.5); } allowedExtremaRegion.ShrinkByRadius(distanceInPixels); } InputImageIndexIteratorType imageIndexIt(inputImage, allowedExtremaRegion); float maxValue = itk::NumericTraits::min(); float minValue = itk::NumericTraits::max(); typename ImageType::IndexType maxIndex; typename ImageType::IndexType minIndex; for(unsigned short i = 0; i < VImageDimension; ++i) { maxIndex[i] = 0; minIndex[i] = 0; } if (maskImage != nullptr) { MaskImageIteratorType maskIt(maskImage, maskImage->GetLargestPossibleRegion()); typename ImageType::IndexType imageIndex; typename ImageType::PointType worldPosition; typename ImageType::IndexType maskIndex; for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt) { imageIndex = maskIndex = maskIt.GetIndex(); if(maskIt.Get() == label) { if( allowedExtremaRegion.IsInside(imageIndex) ) { imageIndexIt.SetIndex( imageIndex ); double value = imageIndexIt.Get(); minMax.Defined = true; //Calculate minimum, maximum and corresponding index-values if( value > maxValue ) { maxIndex = imageIndexIt.GetIndex(); maxValue = value; } if(value < minValue ) { minIndex = imageIndexIt.GetIndex(); minValue = value; } } } } } else { for(imageIndexIt.GoToBegin(); !imageIndexIt.IsAtEnd(); ++imageIndexIt) { double value = imageIndexIt.Get(); minMax.Defined = true; //Calculate minimum, maximum and corresponding index-values if( value > maxValue ) { maxIndex = imageIndexIt.GetIndex(); maxValue = value; } if(value < minValue ) { minIndex = imageIndexIt.GetIndex(); minValue = value; } } } minMax.MaxIndex.set_size(VImageDimension); minMax.MinIndex.set_size(VImageDimension); for(unsigned int i = 0; i < minMax.MaxIndex.size(); ++i) { minMax.MaxIndex[i] = maxIndex[i]; } for(unsigned int i = 0; i < minMax.MinIndex.size(); ++i) { minMax.MinIndex[i] = minIndex[i]; } minMax.Max = maxValue; minMax.Min = minValue; return minMax; } template itk::Size HotspotMaskGenerator::CalculateConvolutionKernelSize( double spacing[VImageDimension], double radiusInMM ) { typedef itk::Image< float, VImageDimension > KernelImageType; typedef typename KernelImageType::SizeType SizeType; SizeType maskSize; for(unsigned int i = 0; i < VImageDimension; ++i) { maskSize[i] = static_cast( 2 * radiusInMM / spacing[i]); // We always want an uneven size to have a clear center point in the convolution mask if(maskSize[i] % 2 == 0 ) { ++maskSize[i]; } } return maskSize; } template itk::SmartPointer< itk::Image > HotspotMaskGenerator::GenerateHotspotSearchConvolutionKernel(double mmPerPixel[VImageDimension], double radiusInMM ) { std::stringstream ss; for (unsigned int i = 0; i < VImageDimension; ++i) { ss << mmPerPixel[i]; if (i < VImageDimension -1) ss << ","; } MITK_DEBUG << "Update convolution kernel for spacing (" << ss.str() << ") and radius " << radiusInMM << "mm"; double radiusInMMSquared = radiusInMM * radiusInMM; typedef itk::Image< float, VImageDimension > KernelImageType; typename KernelImageType::Pointer convolutionKernel = KernelImageType::New(); // Calculate size and allocate mask image typedef typename KernelImageType::SizeType SizeType; SizeType maskSize = this->CalculateConvolutionKernelSize(mmPerPixel, radiusInMM); mitk::Point3D convolutionMaskCenterIndex; convolutionMaskCenterIndex.Fill(0.0); for(unsigned int i = 0; i < VImageDimension; ++i) { convolutionMaskCenterIndex[i] = 0.5 * (double)(maskSize[i]-1); } typedef typename KernelImageType::IndexType IndexType; IndexType maskIndex; maskIndex.Fill(0); typedef typename KernelImageType::RegionType RegionType; RegionType maskRegion; maskRegion.SetSize(maskSize); maskRegion.SetIndex(maskIndex); convolutionKernel->SetRegions(maskRegion); convolutionKernel->SetSpacing(mmPerPixel); convolutionKernel->Allocate(); // Fill mask image values by subsampling the image grid typedef itk::ImageRegionIteratorWithIndex MaskIteratorType; MaskIteratorType maskIt(convolutionKernel,maskRegion); int numberOfSubVoxelsPerDimension = 2; // per dimension! int numberOfSubVoxels = ::pow( static_cast(numberOfSubVoxelsPerDimension), static_cast(VImageDimension) ); double subVoxelSizeInPixels = 1.0 / (double)numberOfSubVoxelsPerDimension; double valueOfOneSubVoxel = 1.0 / (double)numberOfSubVoxels; double maskValue = 0.0; mitk::Point3D subVoxelIndexPosition; double distanceSquared = 0.0; typedef itk::ContinuousIndex ContinuousIndexType; for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt) { ContinuousIndexType indexPoint(maskIt.GetIndex()); mitk::Point3D voxelPosition; for (unsigned int dimension = 0; dimension < VImageDimension; ++dimension) { voxelPosition[dimension] = indexPoint[dimension]; } maskValue = 0.0; mitk::Vector3D subVoxelOffset; subVoxelOffset.Fill(0.0); // iterate sub-voxels by iterating all possible offsets for (subVoxelOffset[0] = -0.5 + subVoxelSizeInPixels / 2.0; subVoxelOffset[0] < +0.5; subVoxelOffset[0] += subVoxelSizeInPixels) { for (subVoxelOffset[1] = -0.5 + subVoxelSizeInPixels / 2.0; subVoxelOffset[1] < +0.5; subVoxelOffset[1] += subVoxelSizeInPixels) { for (subVoxelOffset[2] = -0.5 + subVoxelSizeInPixels / 2.0; subVoxelOffset[2] < +0.5; subVoxelOffset[2] += subVoxelSizeInPixels) { subVoxelIndexPosition = voxelPosition + subVoxelOffset; // this COULD be integrated into the for-loops if neccessary (add voxelPosition to initializer and end condition) distanceSquared = (subVoxelIndexPosition[0]-convolutionMaskCenterIndex[0]) * mmPerPixel[0] * (subVoxelIndexPosition[0]-convolutionMaskCenterIndex[0]) * mmPerPixel[0] + (subVoxelIndexPosition[1]-convolutionMaskCenterIndex[1]) * mmPerPixel[1] * (subVoxelIndexPosition[1]-convolutionMaskCenterIndex[1]) * mmPerPixel[1] + (subVoxelIndexPosition[2]-convolutionMaskCenterIndex[2]) * mmPerPixel[2] * (subVoxelIndexPosition[2]-convolutionMaskCenterIndex[2]) * mmPerPixel[2]; if (distanceSquared <= radiusInMMSquared) { maskValue += valueOfOneSubVoxel; } } } } maskIt.Set( maskValue ); } return convolutionKernel; } template itk::SmartPointer > HotspotMaskGenerator::GenerateConvolutionImage( const itk::Image* inputImage ) { double mmPerPixel[VImageDimension]; for (unsigned int dimension = 0; dimension < VImageDimension; ++dimension) { mmPerPixel[dimension] = inputImage->GetSpacing()[dimension]; } // update convolution kernel typedef itk::Image< float, VImageDimension > KernelImageType; typename KernelImageType::Pointer convolutionKernel = this->GenerateHotspotSearchConvolutionKernel(mmPerPixel, m_HotspotRadiusinMM); // update convolution image typedef itk::Image< TPixel, VImageDimension > InputImageType; typedef itk::Image< TPixel, VImageDimension > ConvolutionImageType; typedef itk::FFTConvolutionImageFilter ConvolutionFilterType; typename ConvolutionFilterType::Pointer convolutionFilter = ConvolutionFilterType::New(); typedef itk::ConstantBoundaryCondition BoundaryConditionType; BoundaryConditionType boundaryCondition; boundaryCondition.SetConstant(0.0); if (m_HotspotMustBeCompletelyInsideImage) { // overwrite default boundary condition convolutionFilter->SetBoundaryCondition(&boundaryCondition); } convolutionFilter->SetInput(inputImage); convolutionFilter->SetKernelImage(convolutionKernel); convolutionFilter->SetNormalize(true); MITK_DEBUG << "Update Convolution image for hotspot search"; convolutionFilter->UpdateLargestPossibleRegion(); typename ConvolutionImageType::Pointer convolutionImage = convolutionFilter->GetOutput(); convolutionImage->SetSpacing( inputImage->GetSpacing() ); // only workaround because convolution filter seems to ignore spacing of input image return convolutionImage; } template < typename TPixel, unsigned int VImageDimension> void HotspotMaskGenerator::FillHotspotMaskPixels( itk::Image* maskImage, itk::Point sphereCenter, double sphereRadiusInMM ) { typedef itk::Image< TPixel, VImageDimension > MaskImageType; typedef itk::ImageRegionIteratorWithIndex MaskImageIteratorType; MaskImageIteratorType maskIt(maskImage, maskImage->GetLargestPossibleRegion()); typename MaskImageType::IndexType maskIndex; typename MaskImageType::PointType worldPosition; // this is not very smart. I would rather use a 0 initialized mask (not the case here -> blame CalculateHotspotMask) and find the region where I need to iterate over, then iterate only over the small region for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt) { maskIndex = maskIt.GetIndex(); maskImage->TransformIndexToPhysicalPoint(maskIndex, worldPosition); maskIt.Set( worldPosition.EuclideanDistanceTo(sphereCenter) <= sphereRadiusInMM ? 1 : 0 ); } } template void HotspotMaskGenerator::CalculateHotspotMask(itk::Image* inputImage, typename itk::Image::Pointer maskImage, unsigned int label) { typedef itk::Image< TPixel, VImageDimension > InputImageType; typedef itk::Image< TPixel, VImageDimension > ConvolutionImageType; typedef itk::Image< float, VImageDimension > KernelImageType; typedef itk::Image< unsigned short, VImageDimension > MaskImageType; typename ConvolutionImageType::Pointer convolutionImage = this->GenerateConvolutionImage(inputImage); if (convolutionImage.IsNull()) { MITK_ERROR << "Empty convolution image in CalculateHotspotStatistics(). We should never reach this state (logic error)."; throw std::logic_error("Empty convolution image in CalculateHotspotStatistics()"); } // if mask image is not defined, create an image of the same size as inputImage and fill it with 1's // there is maybe a better way to do this!? if (maskImage == nullptr) { maskImage = MaskImageType::New(); typename MaskImageType::RegionType maskRegion = inputImage->GetLargestPossibleRegion(); typename MaskImageType::SpacingType maskSpacing = inputImage->GetSpacing(); typename MaskImageType::PointType maskOrigin = inputImage->GetOrigin(); typename MaskImageType::DirectionType maskDirection = inputImage->GetDirection(); maskImage->SetRegions(maskRegion); maskImage->Allocate(); maskImage->SetOrigin(maskOrigin); maskImage->SetSpacing(maskSpacing); maskImage->SetDirection(maskDirection); maskImage->FillBuffer(1); label = 1; } // find maximum in convolution image, given the current mask double requiredDistanceToBorder = m_HotspotMustBeCompletelyInsideImage ? m_HotspotRadiusinMM : -1.0; ImageExtrema convolutionImageInformation = CalculateExtremaWorld(convolutionImage.GetPointer(), maskImage, requiredDistanceToBorder, label); bool isHotspotDefined = convolutionImageInformation.Defined; if (!isHotspotDefined) { MITK_ERROR << "No origin of hotspot-sphere was calculated!"; m_InternalMask = nullptr; } else { // create a binary mask around the "hotspot" region, fill the shape of a sphere around our hotspot center // typename DuplicatorType::Pointer copyMachine = DuplicatorType::New(); // copyMachine->SetInputImage(inputImage); // copyMachine->Update(); // typename CastFilterType::Pointer caster = CastFilterType::New(); // caster->SetInput( copyMachine->GetOutput() ); // caster->Update(); typename MaskImageType::Pointer hotspotMaskITK = MaskImageType::New(); hotspotMaskITK->SetOrigin(inputImage->GetOrigin()); hotspotMaskITK->SetSpacing(inputImage->GetSpacing()); hotspotMaskITK->SetLargestPossibleRegion(inputImage->GetLargestPossibleRegion()); hotspotMaskITK->SetBufferedRegion(inputImage->GetBufferedRegion()); hotspotMaskITK->SetDirection(inputImage->GetDirection()); hotspotMaskITK->SetNumberOfComponentsPerPixel(inputImage->GetNumberOfComponentsPerPixel()); hotspotMaskITK->Allocate(); hotspotMaskITK->FillBuffer(1); typedef typename InputImageType::IndexType IndexType; IndexType maskCenterIndex; for (unsigned int d =0; d< VImageDimension;++d) { maskCenterIndex[d]=convolutionImageInformation.MaxIndex[d]; } typename ConvolutionImageType::PointType maskCenter; inputImage->TransformIndexToPhysicalPoint(maskCenterIndex,maskCenter); FillHotspotMaskPixels(hotspotMaskITK.GetPointer(), maskCenter, m_HotspotRadiusinMM); //obtain mitk::Image::Pointer from itk::Image mitk::Image::Pointer hotspotMaskAsMITKImage = mitk::GrabItkImageMemory(hotspotMaskITK); m_InternalMask = hotspotMaskAsMITKImage; m_ConvolutionImageMaxIndex = convolutionImageInformation.MaxIndex; m_ConvolutionImageMinIndex = convolutionImageInformation.MinIndex; } } bool HotspotMaskGenerator::IsUpdateRequired() const { unsigned long thisClassTimeStamp = this->GetMTime(); unsigned long internalMaskTimeStamp = m_InternalMask->GetMTime(); unsigned long maskGeneratorTimeStamp = m_Mask->GetMTime(); - unsigned long inputImageTimeStamp = m_Image->GetMTime(); + unsigned long inputImageTimeStamp = m_inputImage->GetMTime(); if (thisClassTimeStamp > m_InternalMaskUpdateTime) // inputs have changed { return true; } if (m_InternalMaskUpdateTime < maskGeneratorTimeStamp || m_InternalMaskUpdateTime < inputImageTimeStamp) // mask image has changed outside of this class { return true; } if (internalMaskTimeStamp > m_InternalMaskUpdateTime) // internal mask has been changed outside of this class { return true; } return false; } } diff --git a/Modules/ImageStatistics/mitkHotspotMaskGenerator.h b/Modules/ImageStatistics/mitkHotspotMaskGenerator.h index 5e7886c72f..909fb6b1a6 100644 --- a/Modules/ImageStatistics/mitkHotspotMaskGenerator.h +++ b/Modules/ImageStatistics/mitkHotspotMaskGenerator.h @@ -1,165 +1,164 @@ #ifndef MITKHOTSPOTCALCULATOR_H #define MITKHOTSPOTCALCULATOR_H #include #include #include #include #include #include #include #include namespace mitk { /** * @brief The HotspotMaskGenerator class is used when a hotspot has to be found in an image. A hotspot is * the region of the image where the mean intensity is maximal (=brightest spot). It is usually used in PET scans. * The identification of the hotspot is done as follows: First a cubic (or circular, if image is 2d) * mask of predefined size is generated. This mask is then convolved with the input image (in fourier domain). * The maximum value of the convolved image then corresponds to the hotspot. * If a maskGenerator is set, only the pixels of the convolved image where the corresponding mask is == @a label * are searched for the maximum value. */ class MITKIMAGESTATISTICS_EXPORT HotspotMaskGenerator: public MaskGenerator { public: /** Standard Self typedef */ typedef HotspotMaskGenerator Self; typedef MaskGenerator 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(HotspotMaskGenerator, MaskGenerator) /** - @brief Set the input image + @brief Set the input image. Required for this class */ - void SetImage(mitk::Image::Pointer inputImage); + void SetInputImage(mitk::Image::Pointer inputImage); /** @brief Set a mask (can be nullptr if no mask is desired) */ void SetMask(MaskGenerator::Pointer mask); /** @brief Set the radius of the hotspot (in MM) */ void SetHotspotRadiusInMM(double radiusInMillimeter); const double& GetHotspotRadiusinMM() const; /** @brief Define whether the hotspot must be completely inside the image. Default is true */ void SetHotspotMustBeCompletelyInsideImage(bool hotspotCompletelyInsideImage); bool GetHotspotMustBeCompletelyInsideImage() const; /** @brief If a maskGenerator is set, this detemines which mask value is used */ void SetLabel(unsigned short label); /** @brief Computes and returns the hotspot mask. The hotspot mask has the same size as the input image. The hopspot has value 1, the remaining pixels are set to 0 */ mitk::Image::Pointer GetMask(); /** @brief Returns the image index where the hotspot is located */ vnl_vector GetHotspotIndex(); /** @brief Returns the index where the convolution image is minimal (darkest spot in image) */ vnl_vector GetConvolutionImageMinIndex(); protected: HotspotMaskGenerator(); ~HotspotMaskGenerator(); class ImageExtrema { public: bool Defined; double Max; double Min; vnl_vector MaxIndex; vnl_vector MinIndex; ImageExtrema() :Defined(false) ,Max(itk::NumericTraits::min()) ,Min(itk::NumericTraits::max()) { } }; private: /** \brief Returns size of convolution kernel depending on spacing and radius. */ template itk::Size CalculateConvolutionKernelSize(double spacing[VImageDimension], double radiusInMM); /** \brief Generates image of kernel which is needed for convolution. */ template itk::SmartPointer< itk::Image > GenerateHotspotSearchConvolutionKernel(double spacing[VImageDimension], double radiusInMM); /** \brief Convolves image with spherical kernel image. Used for hotspot calculation. */ template itk::SmartPointer< itk::Image > GenerateConvolutionImage( const itk::Image* inputImage ); /** \brief Fills pixels of the spherical hotspot mask. */ template < typename TPixel, unsigned int VImageDimension> void FillHotspotMaskPixels( itk::Image* maskImage, itk::Point sphereCenter, double sphereRadiusInMM); /** \brief */ template void CalculateHotspotMask(itk::Image* inputImage, typename itk::Image::Pointer maskImage, unsigned int label); template ImageExtrema CalculateExtremaWorld( const itk::Image* inputImage, typename itk::Image::Pointer maskImage, double neccessaryDistanceToImageBorderInMM, unsigned int label); bool IsUpdateRequired() const; HotspotMaskGenerator(const HotspotMaskGenerator &); HotspotMaskGenerator & operator=(const HotspotMaskGenerator &); - mitk::Image::Pointer m_Image; MaskGenerator::Pointer m_Mask; mitk::Image::Pointer m_internalImage; itk::Image::Pointer m_internalMask2D; itk::Image::Pointer m_internalMask3D; double m_HotspotRadiusinMM; bool m_HotspotMustBeCompletelyInsideImage; bool m_HotspotParamsChanged; unsigned short m_Label; vnl_vector m_ConvolutionImageMinIndex, m_ConvolutionImageMaxIndex; unsigned long m_InternalMaskUpdateTime; }; } #endif // MITKHOTSPOTCALCULATOR diff --git a/Modules/ImageStatistics/mitkIgnorePixelMaskGenerator.cpp b/Modules/ImageStatistics/mitkIgnorePixelMaskGenerator.cpp index 0654055e50..ba068667dd 100644 --- a/Modules/ImageStatistics/mitkIgnorePixelMaskGenerator.cpp +++ b/Modules/ImageStatistics/mitkIgnorePixelMaskGenerator.cpp @@ -1,121 +1,112 @@ #include #include #include #include #include #include namespace mitk { -void IgnorePixelMaskGenerator::SetImage(Image::Pointer image) -{ - if (m_Image != image) - { - m_Image = image; - this->Modified(); - } -} - void IgnorePixelMaskGenerator::SetIgnoredPixelValue(RealType pixelValue) { if (pixelValue != m_IgnoredPixelValue) { m_IgnoredPixelValue = pixelValue; this->Modified(); } } mitk::Image::Pointer IgnorePixelMaskGenerator::GetMask() { if (IsUpdateRequired()) { - if (m_Image.IsNull()) + if (m_inputImage.IsNull()) { MITK_ERROR << "Image not set!"; } if (m_IgnoredPixelValue == std::numeric_limits::min()) { MITK_ERROR << "IgnotePixelValue not set!"; } - if (m_TimeStep > (m_Image->GetTimeSteps() - 1)) + if (m_TimeStep > (m_inputImage->GetTimeSteps() - 1)) { - MITK_ERROR << "Invalid time step: " << m_TimeStep << ". The image has " << m_Image->GetTimeSteps() << " timeSteps!"; + MITK_ERROR << "Invalid time step: " << m_TimeStep << ". The image has " << m_inputImage->GetTimeSteps() << " timeSteps!"; } // extractimage time slice ImageTimeSelector::Pointer imgTimeSel = ImageTimeSelector::New(); - imgTimeSel->SetInput(m_Image); + imgTimeSel->SetInput(m_inputImage); imgTimeSel->SetTimeNr(m_TimeStep); imgTimeSel->UpdateLargestPossibleRegion(); mitk::Image::Pointer timeSliceImage = imgTimeSel->GetOutput(); // update m_InternalMask AccessByItk(timeSliceImage, InternalCalculateMask); m_InternalMask->SetGeometry(timeSliceImage->GetGeometry()); this->Modified(); } m_InternalMaskUpdateTime = m_InternalMask->GetMTime(); return m_InternalMask; } template void IgnorePixelMaskGenerator::InternalCalculateMask(typename itk::Image* image) { typedef itk::Image ImageType; typedef itk::Image MaskType; typename MaskType::Pointer mask = MaskType::New(); mask->SetOrigin(image->GetOrigin()); mask->SetSpacing(image->GetSpacing()); mask->SetLargestPossibleRegion(image->GetLargestPossibleRegion()); mask->SetBufferedRegion(image->GetBufferedRegion()); mask->SetDirection(image->GetDirection()); mask->SetNumberOfComponentsPerPixel(image->GetNumberOfComponentsPerPixel()); mask->Allocate(); mask->FillBuffer(1); // iterate over image and mask and set mask=1 if image=m_IgnorePixelValue itk::ImageRegionConstIterator imageIterator(image, image->GetLargestPossibleRegion()); itk::ImageRegionIterator maskIterator(mask, mask->GetLargestPossibleRegion()); for (imageIterator.GoToBegin(); !imageIterator.IsAtEnd(); ++imageIterator, ++maskIterator) { if (imageIterator.Value() == static_cast(m_IgnoredPixelValue)) { maskIterator.Set(0); } } m_InternalMask = GrabItkImageMemory(mask); } bool IgnorePixelMaskGenerator::IsUpdateRequired() const { unsigned long thisClassTimeStamp = this->GetMTime(); unsigned long internalMaskTimeStamp = m_InternalMask->GetMTime(); - unsigned long inputImageTimeStamp = m_Image->GetMTime(); + unsigned long inputImageTimeStamp = m_inputImage->GetMTime(); if (thisClassTimeStamp > m_InternalMaskUpdateTime) // inputs have changed { return true; } if (m_InternalMaskUpdateTime < inputImageTimeStamp) // mask image has changed outside of this class { return true; } if (internalMaskTimeStamp > m_InternalMaskUpdateTime) // internal mask has been changed outside of this class { return true; } return false; } } // end namespace diff --git a/Modules/ImageStatistics/mitkIgnorePixelMaskGenerator.h b/Modules/ImageStatistics/mitkIgnorePixelMaskGenerator.h index 0f6a5ffb5b..d8eecde982 100644 --- a/Modules/ImageStatistics/mitkIgnorePixelMaskGenerator.h +++ b/Modules/ImageStatistics/mitkIgnorePixelMaskGenerator.h @@ -1,73 +1,67 @@ #ifndef MITKIGNOREPIXELMASKGEN_ #define MITKIGNOREPIXELMASKGEN_ #include #include #include #include #include namespace mitk { /** - * @brief The IgnorePixelMaskGenerator class is used to generate a mask that is zero for specific pixel values in the input image. + * @brief The IgnorePixelMaskGenerator class is used to generate a mask that is zero for specific pixel values in the input image. This class requires an input image. */ class MITKIMAGESTATISTICS_EXPORT IgnorePixelMaskGenerator: public MaskGenerator { public: /** Standard Self typedef */ typedef IgnorePixelMaskGenerator Self; typedef MaskGenerator Superclass; typedef itk::SmartPointer< Self > Pointer; typedef itk::SmartPointer< const Self > ConstPointer; typedef double RealType; /** Method for creation through the object factory. */ itkNewMacro(Self) /** Runtime information support. */ itkTypeMacro(IgnorePixelMaskGenerator, MaskGenerator) - /** - * @brief Set input image - */ - void SetImage(mitk::Image::Pointer image); - /** * @brief The mask will be 0 there inputImage==pixelValue and 1 otherwise */ void SetIgnoredPixelValue(RealType pixelValue); /** * @brief Computes and returns the mask */ mitk::Image::Pointer GetMask(); protected: IgnorePixelMaskGenerator(): m_IgnoredPixelValue(std::numeric_limits::min()) { m_TimeStep = 0; m_InternalMaskUpdateTime = 0; m_InternalMask = mitk::Image::New(); } ~IgnorePixelMaskGenerator(){} template void InternalCalculateMask(typename itk::Image* image); private: bool IsUpdateRequired() const; - mitk::Image::Pointer m_Image; RealType m_IgnoredPixelValue; unsigned long m_InternalMaskUpdateTime; }; } #endif diff --git a/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp b/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp index ab20c3e48f..6c5c204c4f 100644 --- a/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp +++ b/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp @@ -1,603 +1,629 @@ #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; } ImageStatisticsCalculator::StatisticsContainer::Pointer ImageStatisticsCalculator::GetStatistics(unsigned int timeStep, unsigned int label) { if (timeStep >= m_StatisticsByTimeStep.size()) { mitkThrow() << "invalid timeStep in ImageStatisticsCalculator_v2::GetStatistics"; } if (m_Image.IsNull()) { mitkThrow() << "no image"; } if (!m_Image->IsInitialized()) { mitkThrow() << "Image not initialized!"; } 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_Image); + 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) } else { // 2) calculate statistics masked AccessByItk_1(m_ImageTimeSlice, InternalCalculateStatisticsMasked, timeStep) } this->Modified(); } m_StatisticsUpdateTimePerTimeStep[timeStep] = m_StatisticsByTimeStep[timeStep][m_StatisticsByTimeStep[timeStep].size()-1]->GetMTime(); 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->Clone(); } } // 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(); } template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator::InternalCalculateStatisticsUnmasked( typename itk::Image< TPixel, VImageDimension >* image, unsigned int timeStep) { typedef typename itk::Image< TPixel, VImageDimension > ImageType; typedef typename itk::ExtendedStatisticsImageFilter 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); //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(); } // 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::InternalCalculateStatisticsMasked( typename itk::Image< TPixel, VImageDimension >* image, 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 (itk::ExceptionObject & e) { // 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()) { 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(); 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++) + 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); + + //typename ImageType::IndexType tmpMinIndex = minMaxFilter->GetMinIndex(*it); + //typename ImageType::IndexType tmpMaxIndex = minMaxFilter->GetMaxIndex(*it); + + //minIndex.set_size(tmpMaxIndex.GetIndexDimension()); + //maxIndex.set_size(tmpMaxIndex.GetIndexDimension()); + 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] = 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]; + //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]; + minIndex[i] = indexCoordinateMin[i]; + maxIndex[i] = indexCoordinateMax[i]; } statisticsResult->SetMinIndex(minIndex); statisticsResult->SetMaxIndex(maxIndex); // just debug assert(abs(minMaxFilter->GetMax(*it) - imageStatisticsFilter->GetMaximum(*it)) < mitk::eps); assert(abs(minMaxFilter->GetMin(*it) - imageStatisticsFilter->GetMinimum(*it)) < mitk::eps); 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; } // 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; } ImageStatisticsCalculator::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::statisticsMapType ImageStatisticsCalculator::StatisticsContainer::GetStatisticsAsMap() { ImageStatisticsCalculator::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; statisticsAsMap["Label"] = m_Label; return statisticsAsMap; } void ImageStatisticsCalculator::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(); m_Histogram = HistogramType::New(); m_minIndex.set_size(0); m_maxIndex.set_size(0); m_Label = 0; } void ImageStatisticsCalculator::StatisticsContainer::Print() { ImageStatisticsCalculator::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::StatisticsContainer::GetAsString() { std::string res = ""; ImageStatisticsCalculator::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/mitkImageStatisticsCalculator.h b/Modules/ImageStatistics/mitkImageStatisticsCalculator.h index 3a6ebda76d..d21f8a02b1 100644 --- a/Modules/ImageStatistics/mitkImageStatisticsCalculator.h +++ b/Modules/ImageStatistics/mitkImageStatisticsCalculator.h @@ -1,413 +1,414 @@ #ifndef MITKIMAGESTATISTICSCALCULATOR #define MITKIMAGESTATISTICSCALCULATOR #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 Container class for storing the computed image statistics. 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) - Variance - Std (Standard Deviation) - Min - Max - RMS (Root Mean Square) - Label (if applicable, the label (unsigned short) of the mask the statistics belong to) - Entropy It furthermore stores the following: - MinIndex (Index of Image where the Minimum is located) - MaxIndex (Index of Image where the Maximum is located) - Histogram of Pixel Values*/ 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; /**Documentation @brief Returns a std::map containing all real valued statistics stored in this class (= all statistics except minIndex, maxIndex and the histogram)*/ statisticsMapType GetStatisticsAsMap(); /**Documentation @brief Deletes all stored values*/ 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; } vnl_vector GetMinIndex() const { return m_minIndex; } void SetMaxIndex(vnl_vector maxIndex) { m_maxIndex = maxIndex; } 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; } /**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(); /**Documentation @brief Generates a string that contains all real valued statistics stored in this class (= all statistics except minIndex, maxIndex and the histogram)*/ std::string GetAsString(); protected: StatisticsContainer(); private: itk::LightObject::Pointer 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->SetRMS(this->GetRMS()); rval->SetSkewness(this->GetSkewness()); rval->SetStd(this->GetStd()); rval->SetUniformity(this->GetUniformity()); rval->SetUPP(this->GetUPP()); rval->SetVariance(this->GetVariance()); rval->SetHistogram(this->GetHistogram()); rval->SetMinIndex(this->GetMinIndex()); rval->SetMaxIndex(this->GetMaxIndex()); return ioPtr; } // 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; }; /**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); protected: ImageStatisticsCalculator(){ m_nBinsForHistogramStatistics = 100; m_binSizeForHistogramStatistics = 10; m_UseBinSizeOverNBins = false; }; private: 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); bool IsUpdateRequired(unsigned int timeStep) const; std::string GetNameOfClass() { return std::string("ImageStatisticsCalculator_v2"); } 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_StatisticsUpdateTimePerTimeStep; }; } #endif // MITKIMAGESTATISTICSCALCULATOR diff --git a/Modules/ImageStatistics/mitkMaskGenerator.cpp b/Modules/ImageStatistics/mitkMaskGenerator.cpp index ce37102ccd..81c3e3a7aa 100644 --- a/Modules/ImageStatistics/mitkMaskGenerator.cpp +++ b/Modules/ImageStatistics/mitkMaskGenerator.cpp @@ -1,37 +1,51 @@ #include namespace mitk { void MaskGenerator::SetTimeStep(unsigned int timeStep) { if (m_TimeStep != timeStep) { m_TimeStep = timeStep; this->Modified(); } } MaskGenerator::MaskGenerator(): m_TimeStep(0) { - + m_inputImage = nullptr; } mitk::Image::Pointer MaskGenerator::GetMask() { return mitk::Image::New(); } //typename itk::Region<3>::Pointer MaskGenerator::GetImageRegionOfMask(Image::Pointer image) //{ // if (m_InternalMask.IsNull() || m_Modified) // { // MITK_ERROR << "Update MaskGenerator first!"; // } // mitk::BaseGeometry::Pointer imageGeometry = image->GetGeometry(); // mitk::BaseGeometry::Pointer maskGeometry = m_InternalMask->GetGeometry(); //} + +void MaskGenerator::SetInputImage(mitk::Image::Pointer inputImg) +{ + if (inputImg != m_inputImage) + { + m_inputImage = inputImg; + this->Modified(); + } +} + +mitk::Image::Pointer MaskGenerator::GetReferenceImage() +{ + return m_inputImage; +} } diff --git a/Modules/ImageStatistics/mitkMaskGenerator.h b/Modules/ImageStatistics/mitkMaskGenerator.h index cf62263dc5..627c7bc4e1 100644 --- a/Modules/ImageStatistics/mitkMaskGenerator.h +++ b/Modules/ImageStatistics/mitkMaskGenerator.h @@ -1,59 +1,70 @@ #ifndef MITKMASKGENERATOR #define MITKMASKGENERATOR #include #include #include #include #include namespace mitk { /** * \class MaskGenerator * \brief Base Class for all Mask Generators. Mask generators are classes that provide functionality for the * creation of binary (or unsigned short) masks that can be applied to an image. See dervied classes for more * information. */ class MITKIMAGESTATISTICS_EXPORT MaskGenerator: public itk::Object { public: /** Standard Self typedef */ typedef MaskGenerator 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(MaskGenerator, itk::Object) //~MaskGenerator(); /** * @brief GetMask must be overridden by derived classes. * @return mitk::Image::Pointer of generated mask */ virtual mitk::Image::Pointer GetMask(); /** * @brief SetTimeStep is used to set the time step for which the mask is to be generated * @param timeStep */ void SetTimeStep(unsigned int timeStep); + /** + * @brief GetReferenceImage per default returns the inputImage (as set by SetInputImage). If no input image is set it will return a nullptr. + */ + virtual mitk::Image::Pointer GetReferenceImage(); + + /** + * @brief SetInputImage is used to set the input image to the mask generator. Some subclasses require an input image, others don't. See the documentation of the specific Mask Generator for more information. + */ + void SetInputImage(mitk::Image::Pointer inputImg); + protected: MaskGenerator(); unsigned int m_TimeStep; mitk::Image::Pointer m_InternalMask; + mitk::Image::Pointer m_inputImage; private: }; } #endif // MITKMASKGENERATOR diff --git a/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.cpp b/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.cpp index 32815e803e..487e463cf9 100644 --- a/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.cpp +++ b/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.cpp @@ -1,429 +1,418 @@ #include #include #include #include "mitkImageAccessByItk.h" #include #include #include +#include #include #include #include #include #include #include #include #include #include #include #include namespace mitk { void PlanarFigureMaskGenerator::SetPlanarFigure(mitk::PlanarFigure::Pointer planarFigure) { if ( planarFigure.IsNull() ) { throw std::runtime_error( "Error: planar figure empty!" ); } if ( !planarFigure->IsClosed() ) { throw std::runtime_error( "Masking not possible for non-closed figures" ); } const PlaneGeometry *planarFigurePlaneGeometry = planarFigure->GetPlaneGeometry(); if ( planarFigurePlaneGeometry == nullptr ) { throw std::runtime_error( "Planar-Figure not yet initialized!" ); } const PlaneGeometry *planarFigureGeometry = dynamic_cast< const PlaneGeometry * >( planarFigurePlaneGeometry ); if ( planarFigureGeometry == nullptr ) { throw std::runtime_error( "Non-planar planar figures not supported!" ); } if (planarFigure != m_PlanarFigure) { this->Modified(); m_PlanarFigure = planarFigure; } } -void PlanarFigureMaskGenerator::SetImage(mitk::Image::Pointer image) +mitk::Image::Pointer PlanarFigureMaskGenerator::GetReferenceImage() { - // check dimension - unsigned int dimension = image->GetDimension(); - if (dimension > 3) - { - MITK_ERROR << "unsupported image dimension"; - } - - const BaseGeometry *imageGeometry = image->GetGeometry(); - if ( imageGeometry == nullptr ) - { - throw std::runtime_error( "Image geometry invalid!" ); - } - - if (image != m_InternalInputImage) + if (IsUpdateRequired()) { - this->Modified(); - m_InternalInputImage = image; + this->CalculateMask(); } + return m_ReferenceImage; } - template < typename TPixel, unsigned int VImageDimension > void PlanarFigureMaskGenerator::InternalCalculateMaskFromPlanarFigure( const itk::Image< TPixel, VImageDimension > *image, unsigned int axis ) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< unsigned short, 2 > MaskImage2DType; -// typedef itk::CastImageFilter< ImageType, MaskImage2DType > CastFilterType; - - // Generate mask image as new image with same header as input image and - // initialize with 1. -// typename CastFilterType::Pointer castFilter = CastFilterType::New(); -// castFilter->SetInput( image ); -// castFilter->Update(); -// castFilter->GetOutput()->FillBuffer( 1 ); + typename MaskImage2DType::Pointer maskImage = MaskImage2DType::New(); maskImage->SetOrigin(image->GetOrigin()); maskImage->SetSpacing(image->GetSpacing()); maskImage->SetLargestPossibleRegion(image->GetLargestPossibleRegion()); maskImage->SetBufferedRegion(image->GetBufferedRegion()); maskImage->SetDirection(image->GetDirection()); maskImage->SetNumberOfComponentsPerPixel(image->GetNumberOfComponentsPerPixel()); maskImage->Allocate(); maskImage->FillBuffer(1); // all PolylinePoints of the PlanarFigure are stored in a vtkPoints object. // These points are used by the vtkLassoStencilSource to create // a vtkImageStencil. const mitk::PlaneGeometry *planarFigurePlaneGeometry = m_PlanarFigure->GetPlaneGeometry(); const typename PlanarFigure::PolyLineType planarFigurePolyline = m_PlanarFigure->GetPolyLine( 0 ); - const mitk::BaseGeometry *imageGeometry3D = m_InternalInputImage->GetGeometry( 0 ); + const mitk::BaseGeometry *imageGeometry3D = m_inputImage->GetGeometry( 0 ); // If there is a second poly line in a closed planar figure, treat it as a hole. PlanarFigure::PolyLineType planarFigureHolePolyline; if (m_PlanarFigure->GetPolyLinesSize() == 2) planarFigureHolePolyline = m_PlanarFigure->GetPolyLine(1); // Determine x- and y-dimensions depending on principal axis // TODO use plane geometry normal to determine that automatically, then check whether the PF is aligned with one of the three principal axis int i0, i1; switch ( axis ) { case 0: i0 = 1; i1 = 2; break; case 1: i0 = 0; i1 = 2; break; case 2: default: i0 = 0; i1 = 1; break; } // store the polyline contour as vtkPoints object bool outOfBounds = false; vtkSmartPointer points = vtkSmartPointer::New(); typename PlanarFigure::PolyLineType::const_iterator it; for ( it = planarFigurePolyline.begin(); it != planarFigurePolyline.end(); ++it ) { Point3D point3D; // Convert 2D point back to the local index coordinates of the selected // image // Fabian: From PlaneGeometry documentation: // Converts a 2D point given in mm (pt2d_mm) relative to the upper-left corner of the geometry into the corresponding world-coordinate (a 3D point in mm, pt3d_mm). // To convert a 2D point given in units (e.g., pixels in case of an image) into a 2D point given in mm (as required by this method), use IndexToWorld. planarFigurePlaneGeometry->Map( *it, point3D ); // Polygons (partially) outside of the image bounds can not be processed // further due to a bug in vtkPolyDataToImageStencil if ( !imageGeometry3D->IsInside( point3D ) ) { outOfBounds = true; } imageGeometry3D->WorldToIndex( point3D, point3D ); points->InsertNextPoint( point3D[i0], point3D[i1], 0 ); } vtkSmartPointer holePoints = nullptr; if (!planarFigureHolePolyline.empty()) { holePoints = vtkSmartPointer::New(); Point3D point3D; PlanarFigure::PolyLineType::const_iterator end = planarFigureHolePolyline.end(); for (it = planarFigureHolePolyline.begin(); it != end; ++it) { // Fabian: same as above planarFigurePlaneGeometry->Map(*it, point3D); imageGeometry3D->WorldToIndex(point3D, point3D); holePoints->InsertNextPoint(point3D[i0], point3D[i1], 0); } } // mark a malformed 2D planar figure ( i.e. area = 0 ) as out of bounds // this can happen when all control points of a rectangle lie on the same line = two of the three extents are zero double bounds[6] = {0, 0, 0, 0, 0, 0}; points->GetBounds( bounds ); bool extent_x = (fabs(bounds[0] - bounds[1])) < mitk::eps; bool extent_y = (fabs(bounds[2] - bounds[3])) < mitk::eps; bool extent_z = (fabs(bounds[4] - bounds[5])) < mitk::eps; // throw an exception if a closed planar figure is deformed, i.e. has only one non-zero extent if ( m_PlanarFigure->IsClosed() && ((extent_x && extent_y) || (extent_x && extent_z) || (extent_y && extent_z))) { mitkThrow() << "Figure has a zero area and cannot be used for masking."; } if ( outOfBounds ) { throw std::runtime_error( "Figure at least partially outside of image bounds!" ); } // create a vtkLassoStencilSource and set the points of the Polygon vtkSmartPointer lassoStencil = vtkSmartPointer::New(); lassoStencil->SetShapeToPolygon(); lassoStencil->SetPoints( points ); vtkSmartPointer holeLassoStencil = nullptr; if (holePoints.GetPointer() != nullptr) { holeLassoStencil = vtkSmartPointer::New(); holeLassoStencil->SetShapeToPolygon(); holeLassoStencil->SetPoints(holePoints); } // Export from ITK to VTK (to use a VTK filter) typedef itk::VTKImageImport< MaskImage2DType > ImageImportType; typedef itk::VTKImageExport< MaskImage2DType > ImageExportType; typename ImageExportType::Pointer itkExporter = ImageExportType::New(); itkExporter->SetInput( maskImage ); // itkExporter->SetInput( castFilter->GetOutput() ); vtkSmartPointer vtkImporter = vtkSmartPointer::New(); this->ConnectPipelines( itkExporter, vtkImporter ); // Apply the generated image stencil to the input image vtkSmartPointer imageStencilFilter = vtkSmartPointer::New(); imageStencilFilter->SetInputConnection( vtkImporter->GetOutputPort() ); imageStencilFilter->SetStencilConnection(lassoStencil->GetOutputPort()); imageStencilFilter->ReverseStencilOff(); imageStencilFilter->SetBackgroundValue( 0 ); imageStencilFilter->Update(); vtkSmartPointer holeStencilFilter = nullptr; if (holeLassoStencil.GetPointer() != nullptr) { holeStencilFilter = vtkSmartPointer::New(); holeStencilFilter->SetInputConnection(imageStencilFilter->GetOutputPort()); holeStencilFilter->SetStencilConnection(holeLassoStencil->GetOutputPort()); holeStencilFilter->ReverseStencilOn(); holeStencilFilter->SetBackgroundValue(0); holeStencilFilter->Update(); } // Export from VTK back to ITK vtkSmartPointer vtkExporter = vtkSmartPointer::New(); vtkExporter->SetInputConnection( holeStencilFilter.GetPointer() == nullptr ? imageStencilFilter->GetOutputPort() : holeStencilFilter->GetOutputPort()); vtkExporter->Update(); typename ImageImportType::Pointer itkImporter = ImageImportType::New(); this->ConnectPipelines( vtkExporter, itkImporter ); itkImporter->Update(); typedef itk::ImageDuplicator< ImageImportType::OutputImageType > DuplicatorType; DuplicatorType::Pointer duplicator = DuplicatorType::New(); duplicator->SetInputImage( itkImporter->GetOutput() ); duplicator->Update(); // Store mask m_InternalITKImageMask2D = duplicator->GetOutput(); } bool PlanarFigureMaskGenerator::GetPrincipalAxis( const BaseGeometry *geometry, Vector3D vector, unsigned int &axis ) { vector.Normalize(); for ( unsigned int i = 0; i < 3; ++i ) { Vector3D axisVector = geometry->GetAxisVector( i ); axisVector.Normalize(); if ( fabs( fabs( axisVector * vector ) - 1.0) < mitk::eps ) { axis = i; return true; } } return false; } void PlanarFigureMaskGenerator::CalculateMask() { - if (m_InternalInputImage->GetTimeSteps() > 0) + if (m_inputImage.IsNull()) + { + MITK_ERROR << "Image is not set."; + } + + if (m_PlanarFigure.IsNull()) + { + MITK_ERROR << "PlanarFigure is not set."; + } + + if (m_TimeStep != 0) + { + MITK_WARN << "Multiple TimeSteps are not supported in PlanarFigureMaskGenerator (yet)."; + } + + const BaseGeometry *imageGeometry = m_inputImage->GetGeometry(); + if ( imageGeometry == nullptr ) + { + throw std::runtime_error( "Image geometry invalid!" ); + } + + if (m_inputImage->GetTimeSteps() > 0) { mitk::ImageTimeSelector::Pointer imgTimeSel = mitk::ImageTimeSelector::New(); - imgTimeSel->SetInput(m_InternalInputImage); + imgTimeSel->SetInput(m_inputImage); imgTimeSel->SetTimeNr(m_TimeStep); imgTimeSel->UpdateLargestPossibleRegion(); m_InternalTimeSliceImage = imgTimeSel->GetOutput(); } else { - m_InternalTimeSliceImage = m_InternalInputImage; + m_InternalTimeSliceImage = m_inputImage; } m_InternalITKImageMask2D = nullptr; const PlaneGeometry *planarFigurePlaneGeometry = m_PlanarFigure->GetPlaneGeometry(); const PlaneGeometry *planarFigureGeometry = dynamic_cast< const PlaneGeometry * >( planarFigurePlaneGeometry ); - const BaseGeometry *imageGeometry = m_InternalInputImage->GetGeometry(); + //const BaseGeometry *imageGeometry = m_inputImage->GetGeometry(); // Find principal direction of PlanarFigure in input image unsigned int axis; if ( !this->GetPrincipalAxis( imageGeometry, planarFigureGeometry->GetNormal(), axis ) ) { throw std::runtime_error( "Non-aligned planar figures not supported!" ); } m_PlanarFigureAxis = axis; // Find slice number corresponding to PlanarFigure in input image typename itk::Image< unsigned short, 3 >::IndexType index; imageGeometry->WorldToIndex( planarFigureGeometry->GetOrigin(), index ); unsigned int slice = index[axis]; - // extract image slice which corresponds to the planarFigure and sotre it in m_InternalImageSlice + // extract image slice which corresponds to the planarFigure and store it in m_InternalImageSlice mitk::Image::Pointer inputImageSlice = extract2DImageSlice(axis, slice); - + //mitk::IOUtil::SaveImage(inputImageSlice, "/home/fabian/inputSliceImage.nrrd"); // Compute mask from PlanarFigure AccessFixedDimensionByItk_1(inputImageSlice, InternalCalculateMaskFromPlanarFigure, 2, axis) //convert itk mask to mitk::Image::Pointer and return it mitk::Image::Pointer planarFigureMaskImage; planarFigureMaskImage = mitk::GrabItkImageMemory(m_InternalITKImageMask2D); - planarFigureMaskImage->SetGeometry(inputImageSlice->GetGeometry()); + //mitk::IOUtil::SaveImage(planarFigureMaskImage, "/home/fabian/planarFigureMaskImage.nrrd"); - Convert2Dto3DImageFilter::Pointer sliceTo3DImageConverter = Convert2Dto3DImageFilter::New(); - sliceTo3DImageConverter->SetInput(planarFigureMaskImage); - sliceTo3DImageConverter->Update(); + //Convert2Dto3DImageFilter::Pointer sliceTo3DImageConverter = Convert2Dto3DImageFilter::New(); + //sliceTo3DImageConverter->SetInput(planarFigureMaskImage); + //sliceTo3DImageConverter->Update(); + //mitk::IOUtil::SaveImage(sliceTo3DImageConverter->GetOutput(), "/home/fabian/3DsliceImage.nrrd"); - m_InternalMask = sliceTo3DImageConverter->GetOutput(); + m_ReferenceImage = inputImageSlice; + //mitk::IOUtil::SaveImage(m_ReferenceImage, "/home/fabian/referenceImage.nrrd"); + m_InternalMask = planarFigureMaskImage; } mitk::Image::Pointer PlanarFigureMaskGenerator::GetMask() { if (IsUpdateRequired()) { - if (m_InternalInputImage.IsNull()) - { - MITK_ERROR << "Image is not set."; - } - - if (m_PlanarFigure.IsNull()) - { - MITK_ERROR << "PlanarFigure is not set."; - } - - if (m_TimeStep != 0) - { - MITK_WARN << "Multiple TimeSteps are not supported in PlanarFigureMaskGenerator (yet)."; - } - this->CalculateMask(); this->Modified(); } m_InternalMaskUpdateTime = m_InternalMask->GetMTime(); return m_InternalMask; } mitk::Image::Pointer PlanarFigureMaskGenerator::extract2DImageSlice(unsigned int axis, unsigned int slice) { // Extract slice with given position and direction from image unsigned int dimension = m_InternalTimeSliceImage->GetDimension(); mitk::Image::Pointer imageSlice = mitk::Image::New(); if (dimension == 3) { ExtractImageFilter::Pointer imageExtractor = ExtractImageFilter::New(); imageExtractor->SetInput( m_InternalTimeSliceImage ); imageExtractor->SetSliceDimension( axis ); imageExtractor->SetSliceIndex( slice ); imageExtractor->Update(); imageSlice = imageExtractor->GetOutput(); } else if(dimension == 2) { imageSlice = m_InternalTimeSliceImage; } else { MITK_ERROR << "Unsupported image dimension. Dimension is: " << dimension << ". Only 2D and 3D images are supported."; } return imageSlice; } bool PlanarFigureMaskGenerator::IsUpdateRequired() const { unsigned long thisClassTimeStamp = this->GetMTime(); unsigned long internalMaskTimeStamp = m_InternalMask->GetMTime(); unsigned long planarFigureTimeStamp = m_PlanarFigure->GetMTime(); - unsigned long inputImageTimeStamp = m_InternalInputImage->GetMTime(); + unsigned long inputImageTimeStamp = m_inputImage->GetMTime(); if (thisClassTimeStamp > m_InternalMaskUpdateTime) // inputs have changed { return true; } if (m_InternalMaskUpdateTime < planarFigureTimeStamp || m_InternalMaskUpdateTime < inputImageTimeStamp) // mask image has changed outside of this class { return true; } if (internalMaskTimeStamp > m_InternalMaskUpdateTime) // internal mask has been changed outside of this class { return true; } return false; } } diff --git a/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.h b/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.h index 66e713c4e5..03459b476a 100644 --- a/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.h +++ b/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.h @@ -1,123 +1,125 @@ #ifndef MITKPLANARFIGUREMASKGENERATOR #define MITKPLANARFIGUREMASKGENERATOR #include #include #include #include #include #include #include #include #include #include namespace mitk { /** * \class PlanarFigureMaskGenerator * \brief Derived from MaskGenerator. This class is used to convert a mitk::PlanarFigure into a binary image mask */ class MITKIMAGESTATISTICS_EXPORT PlanarFigureMaskGenerator: public MaskGenerator { public: /** Standard Self typedef */ typedef PlanarFigureMaskGenerator Self; typedef MaskGenerator 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(PlanarFigureMaskGenerator, MaskGenerator) /** * @brief GetMask Computes and returns the mask * @return mitk::Image::Pointer of the generated mask */ mitk::Image::Pointer GetMask(); void SetPlanarFigure(mitk::PlanarFigure::Pointer planarFigure); /** * @brief SetImage: We need an image because we need its geometry in order to give the points of the planar figure the correct world coordinates * @param inputImage */ void SetImage(const mitk::Image::Pointer inputImage); + mitk::Image::Pointer GetReferenceImage(); + protected: PlanarFigureMaskGenerator():Superclass(){ m_InternalMaskUpdateTime = 0; m_InternalMask = mitk::Image::New(); + m_ReferenceImage = nullptr; } private: void CalculateMask(); template < typename TPixel, unsigned int VImageDimension > void InternalCalculateMaskFromPlanarFigure( const itk::Image< TPixel, VImageDimension > *image, unsigned int axis ); mitk::Image::Pointer extract2DImageSlice(unsigned int axis, unsigned int slice); bool GetPrincipalAxis(const BaseGeometry *geometry, Vector3D vector, unsigned int &axis ); /** Connection from ITK to VTK */ template void ConnectPipelines(ITK_Exporter exporter, vtkSmartPointer importer) { importer->SetUpdateInformationCallback(exporter->GetUpdateInformationCallback()); importer->SetPipelineModifiedCallback(exporter->GetPipelineModifiedCallback()); importer->SetWholeExtentCallback(exporter->GetWholeExtentCallback()); importer->SetSpacingCallback(exporter->GetSpacingCallback()); importer->SetOriginCallback(exporter->GetOriginCallback()); importer->SetScalarTypeCallback(exporter->GetScalarTypeCallback()); importer->SetNumberOfComponentsCallback(exporter->GetNumberOfComponentsCallback()); importer->SetPropagateUpdateExtentCallback(exporter->GetPropagateUpdateExtentCallback()); importer->SetUpdateDataCallback(exporter->GetUpdateDataCallback()); importer->SetDataExtentCallback(exporter->GetDataExtentCallback()); importer->SetBufferPointerCallback(exporter->GetBufferPointerCallback()); importer->SetCallbackUserData(exporter->GetCallbackUserData()); } /** Connection from VTK to ITK */ template void ConnectPipelines(vtkSmartPointer exporter, ITK_Importer importer) { importer->SetUpdateInformationCallback(exporter->GetUpdateInformationCallback()); importer->SetPipelineModifiedCallback(exporter->GetPipelineModifiedCallback()); importer->SetWholeExtentCallback(exporter->GetWholeExtentCallback()); importer->SetSpacingCallback(exporter->GetSpacingCallback()); importer->SetOriginCallback(exporter->GetOriginCallback()); importer->SetScalarTypeCallback(exporter->GetScalarTypeCallback()); importer->SetNumberOfComponentsCallback(exporter->GetNumberOfComponentsCallback()); importer->SetPropagateUpdateExtentCallback(exporter->GetPropagateUpdateExtentCallback()); importer->SetUpdateDataCallback(exporter->GetUpdateDataCallback()); importer->SetDataExtentCallback(exporter->GetDataExtentCallback()); importer->SetBufferPointerCallback(exporter->GetBufferPointerCallback()); importer->SetCallbackUserData(exporter->GetCallbackUserData()); } bool IsUpdateRequired() const; mitk::PlanarFigure::Pointer m_PlanarFigure; typename itk::Image::Pointer m_InternalITKImageMask2D; - mitk::Image::Pointer m_InternalInputImage; mitk::Image::Pointer m_InternalTimeSliceImage; - + mitk::Image::Pointer m_ReferenceImage; unsigned int m_PlanarFigureAxis; unsigned long m_InternalMaskUpdateTime; }; } #endif // MITKPLANARFIGUREMASKGENERATOR diff --git a/Plugins/org.mitk.gui.qt.measurementtoolbox/src/internal/QmitkImageStatisticsCalculationThread.cpp b/Plugins/org.mitk.gui.qt.measurementtoolbox/src/internal/QmitkImageStatisticsCalculationThread.cpp index 52c5258834..12d321fe1e 100644 --- a/Plugins/org.mitk.gui.qt.measurementtoolbox/src/internal/QmitkImageStatisticsCalculationThread.cpp +++ b/Plugins/org.mitk.gui.qt.measurementtoolbox/src/internal/QmitkImageStatisticsCalculationThread.cpp @@ -1,266 +1,266 @@ /*=================================================================== 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 "QmitkImageStatisticsCalculationThread.h" //QT headers #include #include #include #include #include QmitkImageStatisticsCalculationThread::QmitkImageStatisticsCalculationThread():QThread(), m_StatisticsImage(nullptr), m_BinaryMask(nullptr), m_PlanarFigureMask(nullptr), m_TimeStep(0), m_IgnoreZeros(false), m_CalculationSuccessful(false), m_StatisticChanged(false), m_HistogramBinSize(10.0), m_UseDefaultNBins(true), m_nBinsForHistogramStatistics(100), m_prioritizeNBinsOverBinSize(true) { } QmitkImageStatisticsCalculationThread::~QmitkImageStatisticsCalculationThread() { } void QmitkImageStatisticsCalculationThread::Initialize( mitk::Image::Pointer image, mitk::Image::Pointer binaryImage, mitk::PlanarFigure::Pointer 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->Clone(); if(binaryImage.IsNotNull()) this->m_BinaryMask = binaryImage->Clone(); if(planarFig.IsNotNull()) this->m_PlanarFigureMask = planarFig->Clone(); } void QmitkImageStatisticsCalculationThread::SetUseDefaultNBins(bool useDefault) { m_UseDefaultNBins = useDefault; } void QmitkImageStatisticsCalculationThread::SetTimeStep( int times ) { this->m_TimeStep = times; } int QmitkImageStatisticsCalculationThread::GetTimeStep() { return this->m_TimeStep; } std::vector QmitkImageStatisticsCalculationThread::GetStatisticsData() { return this->m_StatisticsVector; } mitk::Image::Pointer QmitkImageStatisticsCalculationThread::GetStatisticsImage() { return this->m_StatisticsImage; } void QmitkImageStatisticsCalculationThread::SetIgnoreZeroValueVoxel(bool _arg) { this->m_IgnoreZeros = _arg; } bool QmitkImageStatisticsCalculationThread::GetIgnoreZeroValueVoxel() { return this->m_IgnoreZeros; } void QmitkImageStatisticsCalculationThread::SetHistogramBinSize(double size) { this->m_HistogramBinSize = size; this->m_prioritizeNBinsOverBinSize = false; } double QmitkImageStatisticsCalculationThread::GetHistogramBinSize() const { return this->m_HistogramBinSize; } void QmitkImageStatisticsCalculationThread::SetHistogramNBins(double size) { this->m_nBinsForHistogramStatistics = size; this->m_prioritizeNBinsOverBinSize = true; } double QmitkImageStatisticsCalculationThread::GetHistogramNBins() const { return this->m_nBinsForHistogramStatistics; } std::string QmitkImageStatisticsCalculationThread::GetLastErrorMessage() { return m_message; } QmitkImageStatisticsCalculationThread::HistogramType::Pointer QmitkImageStatisticsCalculationThread::GetTimeStepHistogram(unsigned int t) { if (t >= this->m_HistogramVector.size()) return nullptr; return this->m_HistogramVector[t]; } bool QmitkImageStatisticsCalculationThread::GetStatisticsChangedFlag() { return m_StatisticChanged; } bool QmitkImageStatisticsCalculationThread::GetStatisticsUpdateSuccessFlag() { return m_CalculationSuccessful; } void QmitkImageStatisticsCalculationThread::run() { bool statisticCalculationSuccessful = true; mitk::ImageStatisticsCalculator::Pointer calculator = mitk::ImageStatisticsCalculator::New(); if(this->m_StatisticsImage.IsNotNull()) { calculator->SetInputImage(m_StatisticsImage); } 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); calculator->SetMask(imgMask.GetPointer()); } if(this->m_PlanarFigureMask.IsNotNull()) { mitk::PlanarFigureMaskGenerator::Pointer pfMaskGen = mitk::PlanarFigureMaskGenerator::New(); - pfMaskGen->SetImage(m_StatisticsImage); + pfMaskGen->SetInputImage(m_StatisticsImage); pfMaskGen->SetPlanarFigure(m_PlanarFigureMask); calculator->SetMask(pfMaskGen.GetPointer()); } } catch( const itk::ExceptionObject& e) { MITK_ERROR << "ITK Exception:" << e.what(); statisticCalculationSuccessful = false; } catch( const mitk::Exception& e ) { MITK_ERROR<< "MITK Exception: " << e.what(); statisticCalculationSuccessful = false; } catch ( const std::runtime_error &e ) { 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; } bool statisticChanged = false; if (this->m_IgnoreZeros) { mitk::IgnorePixelMaskGenerator::Pointer ignorePixelValueMaskGen = mitk::IgnorePixelMaskGenerator::New(); ignorePixelValueMaskGen->SetIgnoredPixelValue(0); - ignorePixelValueMaskGen->SetImage(m_StatisticsImage); + ignorePixelValueMaskGen->SetInputImage(m_StatisticsImage); calculator->SetSecondaryMask(ignorePixelValueMaskGen.GetPointer()); } else { calculator->SetSecondaryMask(nullptr); } if (m_UseDefaultNBins) { calculator->SetNBinsForHistogramStatistics(100); } else { if (!m_prioritizeNBinsOverBinSize) { calculator->SetBinSizeForHistogramStatistics(m_HistogramBinSize); } else { calculator->SetNBinsForHistogramStatistics(100); } } //calculator->SetHistogramBinSize( m_HistogramBinSize ); //calculator->SetUseDefaultBinSize( m_UseDefaultBinSize ); 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_HistogramVector.clear(); for (unsigned int i = 0; i < m_StatisticsImage->GetTimeSteps(); i++) { this->m_StatisticsVector.push_back(calculator->GetStatistics(i)); this->m_HistogramVector.push_back((HistogramType*)this->m_StatisticsVector[i]->GetHistogram()); } } }