diff --git a/Modules/DiffusionImaging/MiniApps/CMakeLists.txt b/Modules/DiffusionImaging/MiniApps/CMakeLists.txt index e4564aff4f..21e436f79b 100755 --- a/Modules/DiffusionImaging/MiniApps/CMakeLists.txt +++ b/Modules/DiffusionImaging/MiniApps/CMakeLists.txt @@ -1,120 +1,121 @@ option(BUILD_DiffusionMiniApps "Build commandline tools for diffusion" OFF) if(BUILD_DiffusionMiniApps OR MITK_BUILD_ALL_APPS) # needed include directories include_directories( ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_BINARY_DIR} ) # list of diffusion miniapps # if an app requires additional dependencies # they are added after a "^^" and separated by "_" set( diffusionminiapps DwiDenoising^^ ImageResampler^^ NetworkCreation^^MitkFiberTracking_MitkConnectomics NetworkStatistics^^MitkConnectomics ExportShImage^^ ImageStatisticsMiniApp^^MitkImageStatistics + ImageStatisticsMiniapp_v2^^MitkImageStatistics PlanarFigureMaskCreatorMiniapp^^MitkImageStatistics HotspotMiniApp^^MitkImageStatistics Fiberfox^^MitkFiberTracking MultishellMethods^^MitkFiberTracking PeaksAngularError^^MitkFiberTracking PeakExtraction^^MitkFiberTracking FiberExtraction^^MitkFiberTracking FiberProcessing^^MitkFiberTracking FiberDirectionExtraction^^MitkFiberTracking LocalDirectionalFiberPlausibility^^MitkFiberTracking StreamlineTracking^^MitkFiberTracking GibbsTracking^^MitkFiberTracking CopyGeometry^^ DiffusionIndices^^ TractometerMetrics^^MitkFiberTracking QballReconstruction^^ Registration^^ FileFormatConverter^^MitkFiberTracking TensorReconstruction^^ TensorDerivedMapsExtraction^^ DiffusionDICOMLoader^^ DFTraining^^MitkFiberTracking DFTracking^^MitkFiberTracking ) foreach(diffusionminiapp ${diffusionminiapps}) # extract mini app name and dependencies string(REPLACE "^^" "\\;" miniapp_info ${diffusionminiapp}) set(miniapp_info_list ${miniapp_info}) list(GET miniapp_info_list 0 appname) list(GET miniapp_info_list 1 raw_dependencies) string(REPLACE "_" "\\;" dependencies "${raw_dependencies}") set(dependencies_list ${dependencies}) mitk_create_executable(${appname} DEPENDS MitkCore MitkDiffusionCore MitkCommandLine ${dependencies_list} PACKAGE_DEPENDS ITK CPP_FILES ${appname}.cpp ) if(EXECUTABLE_IS_ENABLED) # On Linux, create a shell script to start a relocatable application if(UNIX AND NOT APPLE) install(PROGRAMS "${MITK_SOURCE_DIR}/CMake/RunInstalledApp.sh" DESTINATION "." RENAME ${EXECUTABLE_TARGET}.sh) endif() get_target_property(_is_bundle ${EXECUTABLE_TARGET} MACOSX_BUNDLE) if(APPLE) if(_is_bundle) set(_target_locations ${EXECUTABLE_TARGET}.app) set(${_target_locations}_qt_plugins_install_dir ${EXECUTABLE_TARGET}.app/Contents/MacOS) set(_bundle_dest_dir ${EXECUTABLE_TARGET}.app/Contents/MacOS) set(_qt_plugins_for_current_bundle ${EXECUTABLE_TARGET}.app/Contents/MacOS) set(_qt_conf_install_dirs ${EXECUTABLE_TARGET}.app/Contents/Resources) install(TARGETS ${EXECUTABLE_TARGET} BUNDLE DESTINATION . ) else() if(NOT MACOSX_BUNDLE_NAMES) set(_qt_conf_install_dirs bin) set(_target_locations bin/${EXECUTABLE_TARGET}) set(${_target_locations}_qt_plugins_install_dir bin) install(TARGETS ${EXECUTABLE_TARGET} RUNTIME DESTINATION bin) else() foreach(bundle_name ${MACOSX_BUNDLE_NAMES}) list(APPEND _qt_conf_install_dirs ${bundle_name}.app/Contents/Resources) set(_current_target_location ${bundle_name}.app/Contents/MacOS/${EXECUTABLE_TARGET}) list(APPEND _target_locations ${_current_target_location}) set(${_current_target_location}_qt_plugins_install_dir ${bundle_name}.app/Contents/MacOS) message( " set(${_current_target_location}_qt_plugins_install_dir ${bundle_name}.app/Contents/MacOS) ") install(TARGETS ${EXECUTABLE_TARGET} RUNTIME DESTINATION ${bundle_name}.app/Contents/MacOS/) endforeach() endif() endif() else() set(_target_locations bin/${EXECUTABLE_TARGET}${CMAKE_EXECUTABLE_SUFFIX}) set(${_target_locations}_qt_plugins_install_dir bin) set(_qt_conf_install_dirs bin) install(TARGETS ${EXECUTABLE_TARGET} RUNTIME DESTINATION bin) endif() endif() endforeach() # This mini app does not depend on mitkDiffusionImaging at all mitk_create_executable(Dicom2Nrrd DEPENDS MitkCore MitkCommandLine CPP_FILES Dicom2Nrrd.cpp ) # On Linux, create a shell script to start a relocatable application if(UNIX AND NOT APPLE) install(PROGRAMS "${MITK_SOURCE_DIR}/CMake/RunInstalledApp.sh" DESTINATION "." RENAME ${EXECUTABLE_TARGET}.sh) endif() if(EXECUTABLE_IS_ENABLED) MITK_INSTALL_TARGETS(EXECUTABLES ${EXECUTABLE_TARGET}) endif() endif() diff --git a/Modules/DiffusionImaging/MiniApps/ImageStatisticsMiniapp_v2.cpp b/Modules/DiffusionImaging/MiniApps/ImageStatisticsMiniapp_v2.cpp new file mode 100644 index 0000000000..4a3c4af7de --- /dev/null +++ b/Modules/DiffusionImaging/MiniApps/ImageStatisticsMiniapp_v2.cpp @@ -0,0 +1,202 @@ +/*=================================================================== + +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 "mitkIOUtil.h" +#include +#include +#include +#include +#include "mitkImageAccessByItk.h" +#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 +void printMap(std::map input) +{ + for (auto it = input.begin(); it != input.end(); ++it) + { + std::cout << it->first<< ": " << it->second<< 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.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: + unsigned int timeStep = 0; + std::string inputImageFile; + if (parsedArgs.count("input") || parsedArgs.count("i")) + { + inputImageFile = us::any_cast(parsedArgs["input"]); + } else + { + inputImageFile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic3D.nrrd"; + } + + 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_v2::Pointer calculator = mitk::ImageStatisticsCalculator_v2::New(); + + std::cout << "calculating statistics (unmasked) itk: " << std::endl; + mitk::ImageStatisticsCalculator_v2::statisticsMapType result; + try + { + calculator->SetInputImage(inputImage); + result = calculator->GetStatistics(timeStep); + } + catch( const itk::ExceptionObject& e) + { + MITK_ERROR << "Statistic Calculation Failed - ITK 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) result["N"]; + } + + // Write Results to file + std::ofstream output; + output.open(outFile.c_str()); + output << result["Mean"] << " , "; + output << result["Std"] << " , "; + output << result["RMS"] << " , "; + output << result["Max"] << " , "; + output << result["Min"] << " , "; + output << result["N"] << " , "; + output << volume << "\n"; + + output.flush(); + output.close(); + + printMap(result); + + std::cout << "calculating statistics boost: " << std::endl; + + statistics_res res; + AccessByItk_n(inputImage, get_statistics_boost, (res)); + + printstats(res); + + return EXIT_SUCCESS; +} diff --git a/Modules/ImageStatistics/mitkImageStatisticsCalculator2.cpp b/Modules/ImageStatistics/mitkImageStatisticsCalculator2.cpp index 03fc78c17d..66bce82291 100644 --- a/Modules/ImageStatistics/mitkImageStatisticsCalculator2.cpp +++ b/Modules/ImageStatistics/mitkImageStatisticsCalculator2.cpp @@ -1,565 +1,618 @@ #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include namespace mitk { void ImageStatisticsCalculator_v2::SetInputImage(mitk::Image::Pointer image) { if (image != m_Image) { m_Image = image; m_HistogramStatisticsUpdateRequiredByTimeStep.resize(m_Image->GetTimeSteps()); m_StatisticsUpdateRequiredByTimeStep.resize(m_Image->GetTimeSteps()); + m_StatisticsByTimeStep.resize(m_Image->GetTimeSteps()); + m_HistogramStatisticsByTimeStep.resize(m_Image->GetTimeSteps()); m_HistogramsByTimeStep.resize(m_Image->GetTimeSteps()); this->SetAllStatisticsToUpdateRequired(); this->SetAllHistogramStatisticsToUpdateRequired(); } } void ImageStatisticsCalculator_v2::SetMask(mitk::MaskGenerator::Pointer mask) { if (mask != m_MaskGenerator) { m_MaskGenerator = mask; this->SetAllStatisticsToUpdateRequired(); this->SetAllHistogramStatisticsToUpdateRequired(); } } void ImageStatisticsCalculator_v2::SetDoHistogramStatistics(bool doHistStatistics) { m_DoHistogramStatistics = doHistStatistics; } bool ImageStatisticsCalculator_v2::GetDoHistogramStatistics() const { return m_DoHistogramStatistics; } void ImageStatisticsCalculator_v2::SetNBinsForHistogramStatistics(unsigned int nBins) { if (nBins != m_nBinsForHistogramStatistics) { m_nBinsForHistogramStatistics = nBins; this->SetAllHistogramStatisticsToUpdateRequired(); } } unsigned int ImageStatisticsCalculator_v2::GetNBinsForHistogramStatistics() const { return m_nBinsForHistogramStatistics; } ImageStatisticsCalculator_v2::statisticsMapType ImageStatisticsCalculator_v2::GetStatistics(unsigned int timeStep, unsigned int label) { if (timeStep >= m_StatisticsByTimeStep.size()) { throw std::runtime_error("invalid timeStep in ImageStatisticsCalculator_v2::GetStatistics"); } if (m_Image.IsNull()) { throw std::runtime_error("no image"); } if (m_MaskGenerator.IsNotNull()) { m_MaskGenerator->SetTimeStep(timeStep); m_InternalMask = m_MaskGenerator->GetMask(); } if (m_StatisticsUpdateRequiredByTimeStep[timeStep]) { // Calculate statistics with/without mask if (m_MaskGenerator.IsNull()) { // 1) calculate statistics unmasked: // plug image into itkstatisticsimagefilter (will be replaced by my awesome filter later on) // retrieve statistics and save them AccessByItk_1(m_Image, InternalCalculateStatisticsUnmasked, timeStep) } else { // 2) calculate statistics masked // extract mask image region // plug mask and image into itklabelstatisticsimagefilter AccessByItk_1(m_Image, InternalCalculateStatisticsMasked, timeStep) } m_StatisticsUpdateRequiredByTimeStep[timeStep]=false; } if(m_DoHistogramStatistics && m_HistogramStatisticsUpdateRequiredByTimeStep[timeStep]) { if (m_MaskGenerator.IsNull()) { // calculate histogram unmasked AccessByItk_3(m_Image, InternalCalculateHistogramUnmasked, m_StatisticsByTimeStep[timeStep][0]->GetMin(), m_StatisticsByTimeStep[timeStep][0]->GetMax(), timeStep); m_HistogramStatisticsByTimeStep[timeStep].resize(1); HistogramStatisticsContainer::Pointer histogramStatistics = this->InternalCalculateHistogramStatistics(m_HistogramsByTimeStep[timeStep][0].histogram); histogramStatistics->GetLabel() = 0; m_HistogramStatisticsByTimeStep[timeStep][0] = histogramStatistics; } else { // the label histograms have already been calculated in InternalCalculateStatisticsMasked for (HistogramContainer histContainer:m_HistogramsByTimeStep[timeStep]) { HistogramStatisticsContainer::Pointer histogramStatistics = this->InternalCalculateHistogramStatistics(histContainer.histogram); histogramStatistics->GetLabel() = histContainer.label; } } m_HistogramStatisticsUpdateRequiredByTimeStep[timeStep] = false; } - // merge histogram and regular statistics, return result + statisticsMapType statistics; + if (m_MaskGenerator.IsNull()) + { + statistics = m_StatisticsByTimeStep[timeStep][0]->GetStatisticsAsMap(); + if (m_DoHistogramStatistics) + { + statisticsMapType histStatistics = m_HistogramStatisticsByTimeStep[timeStep][0]->GetStatisticsAsMap(); + statistics.insert(histStatistics.begin(), histStatistics.end()); + } + } + else + { + bool statisticsFound(false); + for (auto stats: m_StatisticsByTimeStep[timeStep]) + { + + if (stats->GetLabel() == label) + { + statistics = stats->GetStatisticsAsMap(); + statisticsFound = true; + break; + } + } + if (!statisticsFound) + { + throw std::runtime_error("Invalid label"); + } + + if (m_DoHistogramStatistics) + { + bool histStatsFound(false); + for (auto histStats: m_HistogramStatisticsByTimeStep[timeStep]) + { + if (histStats->GetLabel() == label) + { + statisticsMapType histStatsMap = histStats->GetStatisticsAsMap(); + statistics.insert(histStatsMap.begin(), histStatsMap.end()); + histStatsFound = true; + break; + } + } + if (!histStatsFound) + { + throw std::runtime_error("Invalid label"); + } + } + } + return statistics; } void ImageStatisticsCalculator_v2::SetAllStatisticsToUpdateRequired() { for (unsigned int i = 0; i < m_StatisticsUpdateRequiredByTimeStep.size(); i++) { this->SetStatsTimeStepToUpdateRequired(i); } } void ImageStatisticsCalculator_v2::SetStatsTimeStepToUpdateRequired(unsigned int timeStep) { if (timeStep >= m_StatisticsUpdateRequiredByTimeStep.size()) { throw std::runtime_error("invalid timeStep in ImageStatisticsCalculator_v2::SetStatsTimeStepToUpdateRequired"); } m_StatisticsUpdateRequiredByTimeStep[timeStep] = true; m_StatisticsByTimeStep[timeStep].resize(0); } void ImageStatisticsCalculator_v2::SetAllHistogramStatisticsToUpdateRequired() { for (unsigned int i = 0; i < m_HistogramStatisticsUpdateRequiredByTimeStep.size(); i++) { this->SetHistStatsTimeStepToUpdateRequired(i); } } void ImageStatisticsCalculator_v2::SetHistStatsTimeStepToUpdateRequired(unsigned int timeStep) { if (timeStep >= m_HistogramStatisticsUpdateRequiredByTimeStep.size()) { throw std::runtime_error("invalid timeStep in ImageStatisticsCalculator_v2::SetHistStatsTimeStepToUpdateRequired"); } m_HistogramStatisticsUpdateRequiredByTimeStep[timeStep] = true; m_HistogramStatisticsByTimeStep[timeStep].resize(0); } template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator_v2::InternalCheckMaskSanity( typename itk::Image< TPixel, VImageDimension >* image, typename itk::Image< unsigned short, VImageDimension >::Pointer maskImage) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef typename ImageType::PointType PointType; typedef typename ImageType::DirectionType DirectionType; // check direction DirectionType imageDirection = image->GetDirection(); DirectionType maskDirection = maskImage->GetDirection(); for( int i = 0; i < imageDirection.ColumnDimensions; ++i ) { for( int j = 0; j < imageDirection.ColumnDimensions; ++j ) { double differenceDirection = imageDirection[i][j] - maskDirection[i][j]; if ( fabs( differenceDirection ) > mitk::eps ) { double differenceDirection = imageDirection[i][j] - maskDirection[i][j]; if ( fabs( differenceDirection ) > 0.001 /*mitk::eps*/ ) // TODO: temp fix (bug 17121) { itkExceptionMacro( << "Mask needs to have same direction as image! (Image direction: " << imageDirection << "; Mask direction: " << maskDirection << ")" ); } } } } // check spacing PointType imageSpacing = image->GetSpacing(); PointType maskSpacing = maskImage->GetSpacing(); for (unsigned int i = 0; i < VImageDimension; i++) { if ( fabs( maskSpacing[i] - imageSpacing[i] ) > mitk::eps ) { itkExceptionMacro(<< "Spacing of mask and image is not equal. Mask: " << maskSpacing << " image: " << imageSpacing ); } } // check alignment // Make sure that the voxels of mask and image are correctly "aligned", i.e., voxel boundaries are the same in both images PointType imageOrigin = image->GetOrigin(); PointType maskOrigin = maskImage->GetOrigin(); typedef itk::ContinuousIndex ContinousIndexType; ContinousIndexType maskOriginContinousIndex, imageOriginContinousIndex; image->TransformPhysicalPointToContinuousIndex(maskOrigin, maskOriginContinousIndex); image->TransformPhysicalPointToContinuousIndex(imageOrigin, imageOriginContinousIndex); for ( unsigned int i = 0; i < ImageType::ImageDimension; ++i ) { double misalignment = maskOriginContinousIndex[i] - floor( maskOriginContinousIndex[i] + 0.5 ); // misalignment must be a multiple (int) of spacing in that direction if ( fmod(misalignment,imageSpacing[i]) > mitk::eps ) { itkWarningMacro( << "Pixels/voxels of mask and image are not sufficiently aligned! (Misalignment: " << fmod(misalignment,imageSpacing[i]) << ")" ); } } // mask must be completely inside image region // Make sure that mask region is contained within image region if ( maskImage.IsNotNull() && !image->GetLargestPossibleRegion().IsInside( maskImage->GetLargestPossibleRegion() ) ) { itkWarningMacro( << "Mask region needs to be inside of image region! (Image region: " << image->GetLargestPossibleRegion() << "; Mask region: " << maskImage->GetLargestPossibleRegion() << ")" ); } } template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator_v2::InternalCalculateStatisticsUnmasked( typename itk::Image< TPixel, VImageDimension >* image, unsigned int timeStep) { typedef typename itk::Image< TPixel, VImageDimension > ImageType; typename itk::ExtendedStatisticsImageFilter::Pointer statisticsFilter = itk::ExtendedStatisticsImageFilter::New(); statisticsFilter->SetInput(image); statisticsFilter->SetCoordinateTolerance(0.001); statisticsFilter->SetDirectionTolerance(0.001); + // TODO: REMOVE THIS, the histogram statistics should be a separate module + statisticsFilter->SetBinSize(m_nBinsForHistogramStatistics); + 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); StatisticsContainer::Pointer statisticsResult = StatisticsContainer::New(); statisticsResult->GetLabel() = 0; statisticsResult->GetN() = image->GetLargestPossibleRegion().GetNumberOfPixels(); statisticsResult->GetMean() = statisticsFilter->GetMean(); statisticsResult->GetMin() = statisticsFilter->GetMinimum(); statisticsResult->GetMax() = statisticsFilter->GetMaximum(); statisticsResult->GetVariance() = statisticsFilter->GetVariance(); statisticsResult->GetStd() = statisticsFilter->GetSigma(); statisticsResult->GetSkewness() = statisticsFilter->GetSkewness(); statisticsResult->GetKurtosis() = statisticsFilter->GetKurtosis(); //statisticsResult->GetRMS() = statisticsFilter->Get(); statisticsResult->GetMPP() = statisticsFilter->GetMPP(); m_StatisticsByTimeStep[timeStep][0] = statisticsResult; } template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator_v2::InternalCalculateHistogramUnmasked( typename itk::Image< TPixel, VImageDimension >* image, double minVal, double maxVal, unsigned int timeStep) { typedef typename itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Statistics::ImageToHistogramFilter HistogramFilterType; typedef itk::Statistics::Histogram HistogramType; typedef typename HistogramFilterType::HistogramSizeType SizeType; typename HistogramFilterType::Pointer histogramFilter = HistogramFilterType::New(); histogramFilter->SetInput(image); - SizeType histSize = static_cast(m_nBinsForHistogramStatistics); + SizeType histSize(1); + histSize[0] = m_nBinsForHistogramStatistics; histogramFilter->SetHistogramSize(histSize); // is that correct?? typename HistogramFilterType::HistogramMeasurementVectorType lowerBound(1); typename HistogramFilterType::HistogramMeasurementVectorType upperBound(1); lowerBound[0] = minVal; upperBound[0] = maxVal; histogramFilter->SetHistogramBinMinimum(lowerBound); histogramFilter->SetHistogramBinMaximum(upperBound); histogramFilter->Update(); typename HistogramType::Pointer histogram = histogramFilter->GetOutput(); HistogramContainer histogramContainer; histogramContainer.label=0; histogramContainer.histogram = histogram; m_HistogramsByTimeStep[timeStep].resize(1); m_HistogramsByTimeStep[timeStep][0] = histogramContainer; } template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator_v2::ExtractMaskImageRegion( typename itk::Image< TPixel, VImageDimension >* image, typename itk::Image< unsigned short, VImageDimension >::Pointer maskImage, typename itk::Image< TPixel, VImageDimension >::Pointer adaptedImage ) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::ExtractImageFilter< ImageType, ImageType > ExtractImageFilterType; typename ImageType::SizeType imageSize = image->GetBufferedRegion().GetSize(); typename ImageType::SizeType maskSize = maskImage->GetBufferedRegion().GetSize(); bool maskSmallerImage = false; for ( unsigned int i = 0; i < ImageType::ImageDimension; ++i ) { if ( maskSize[i] < imageSize[i] ) { maskSmallerImage = true; } } if ( maskSmallerImage ) { typename ExtractImageFilterType::Pointer extractImageFilter = ExtractImageFilterType::New(); extractImageFilter->SetInput( image ); extractImageFilter->SetExtractionRegion( maskImage->GetBufferedRegion() ); extractImageFilter->SetCoordinateTolerance( 0.001 ); extractImageFilter->SetDirectionTolerance( 0.001 ); extractImageFilter->Update(); adaptedImage = extractImageFilter->GetOutput(); } else { adaptedImage = image; } } /* Taken from original ImageStatisticsCalculator * This function is not ideal on so many levels 1) it does not support more than one label 2) it assumes that the label id is 1 (which does not have to be the case) 3) single threaded => I suggest writing an itk filter that finds min and max value. This filter could then be called from our new statisticsimagefilter in the beforethreadedgenereatedata (provided that a min/max for the histogram is not yet defined) EDIT: I changed this function so that it simply searches the image min and max and completely disregards the mask. Since we call this function after we cropped the image to the mask region we can run it on the LargestPossibleRegion of the image */ template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator_v2::GetMinAndMaxValue( typename itk::Image< TPixel, VImageDimension >* inputImage, double &minVal, double &maxVal ) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::ImageRegionConstIteratorWithIndex ImageRegionIteratorType; ImageRegionIteratorType imageIterator( inputImage, inputImage->GetLargestPossibleRegion() ); maxVal = std::numeric_limits::min(); minVal = std::numeric_limits::max(); long counter(0); double actualPixelValue(0); for( imageIterator.GoToBegin(); !imageIterator.IsAtEnd(); ++imageIterator) { if( imageIterator.Value()== 1.0) { counter++; actualPixelValue = imageIterator.Value(); if(actualPixelValue > maxVal) { maxVal = actualPixelValue; } else if(actualPixelValue < minVal) { minVal = actualPixelValue; } } } } template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator_v2::InternalCalculateStatisticsMasked( typename itk::Image< TPixel, VImageDimension >* image, unsigned int timeStep) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< unsigned short, VImageDimension > MaskType; typedef itk::ExtendedLabelStatisticsImageFilter< ImageType, MaskType > ImageStatisticsFilterType; // maskImage has to have the same dimension as image typename MaskType::Pointer maskImage = MaskType::New(); maskImage = ImageToItkImage< unsigned short, VImageDimension >(m_InternalMask); InternalCheckMaskSanity(image, maskImage); // if mask is smaller than image, extract the image region where the mask is typename ImageType::Pointer adaptedImage = ImageType::New(); // uses m_InternalMask this->ExtractMaskImageRegion(image, maskImage, adaptedImage); double minVal, maxVal; GetMinAndMaxValue(adaptedImage.GetPointer(), minVal, maxVal); typename ImageStatisticsFilterType::Pointer imageStatisticsFilter = ImageStatisticsFilterType::New(); imageStatisticsFilter->SetDirectionTolerance(0.001); imageStatisticsFilter->SetCoordinateTolerance(0.001); imageStatisticsFilter->SetHistogramParameters(m_nBinsForHistogramStatistics, minVal, maxVal); imageStatisticsFilter->SetInput(adaptedImage); imageStatisticsFilter->SetLabelInput(maskImage); imageStatisticsFilter->Update(); std::list labels = imageStatisticsFilter->GetRelevantLabels(); std::list::iterator it = labels.begin(); m_StatisticsByTimeStep[timeStep].resize(0); while(it != labels.end()) { StatisticsContainer::Pointer statisticsResult = StatisticsContainer::New(); statisticsResult->GetN() = imageStatisticsFilter->GetSum(*it) / (double) imageStatisticsFilter->GetMean(*it); statisticsResult->GetMean() = imageStatisticsFilter->GetMean(*it); statisticsResult->GetMin() = imageStatisticsFilter->GetMinimum(*it); statisticsResult->GetMax() = imageStatisticsFilter->GetMaximum(*it); statisticsResult->GetVariance() = imageStatisticsFilter->GetVariance(*it); statisticsResult->GetStd() = imageStatisticsFilter->GetSigma(1); statisticsResult->GetSkewness() = imageStatisticsFilter->GetSkewness(*it); statisticsResult->GetKurtosis() = imageStatisticsFilter->GetKurtosis(*it); //statisticsResult->GetRMS() = imageStatisticsFilter->GetRMS(*it); statisticsResult->GetMPP() = imageStatisticsFilter->GetMPP(*it); statisticsResult->GetLabel() = *it; m_StatisticsByTimeStep[timeStep].push_back(statisticsResult); if (m_DoHistogramStatistics) { HistogramContainer histCont; histCont.label = *it; histCont.histogram = imageStatisticsFilter->GetHistogram(*it); m_HistogramsByTimeStep[timeStep].push_back(histCont); } } } ImageStatisticsCalculator_v2::HistogramStatisticsContainer::Pointer ImageStatisticsCalculator_v2::InternalCalculateHistogramStatistics(itk::Statistics::Histogram::Pointer histogram) { HistogramStatisticsCalculator histStatCalc; histStatCalc.SetHistogram(histogram); histStatCalc.CalculateStatistics(); HistogramStatisticsContainer::Pointer histStatResults = HistogramStatisticsContainer::New(); histStatResults->GetEntropy() = histStatCalc.GetEntropy(); histStatResults->GetMedian() = histStatCalc.GetMedian(); histStatResults->GetUniformity() = histStatCalc.GetUniformity(); histStatResults->GetUPP() = histStatCalc.GetUPP(); return histStatResults; } ImageStatisticsCalculator_v2::StatisticsContainer::StatisticsContainer(): m_N(std::numeric_limits::max()), m_Mean(std::numeric_limits::max()), m_Min(std::numeric_limits::max()), m_Max(std::numeric_limits::max()), m_Std(std::numeric_limits::max()), m_Variance(std::numeric_limits::max()), m_Skewness(std::numeric_limits::max()), m_Kurtosis(std::numeric_limits::max()), m_RMS(std::numeric_limits::max()), m_MPP(std::numeric_limits::max()) { } ImageStatisticsCalculator_v2::statisticsMapType ImageStatisticsCalculator_v2::StatisticsContainer::GetStatisticsAsMap() { ImageStatisticsCalculator_v2::statisticsMapType statisticsAsMap; statisticsAsMap["N"] = m_N; statisticsAsMap["Mean"] = m_Mean; statisticsAsMap["Min"] = m_Min; statisticsAsMap["Max"] = m_Max; statisticsAsMap["StandardDeviation"] = m_Std; statisticsAsMap["Variance"] = m_Variance; statisticsAsMap["Skewness"] = m_Skewness; statisticsAsMap["Kurtosis"] = m_Kurtosis; statisticsAsMap["RMS"] = m_RMS; statisticsAsMap["MPP"] = m_MPP; return statisticsAsMap; } void ImageStatisticsCalculator_v2::StatisticsContainer::Reset() { m_N = std::numeric_limits::max(); m_Mean = std::numeric_limits::max(); m_Min = std::numeric_limits::max(); m_Max = std::numeric_limits::max(); m_Std = std::numeric_limits::max(); m_Variance = std::numeric_limits::max(); m_Skewness = std::numeric_limits::max(); m_Kurtosis = std::numeric_limits::max(); m_RMS = std::numeric_limits::max(); m_MPP = std::numeric_limits::max(); } ImageStatisticsCalculator_v2::HistogramStatisticsContainer::HistogramStatisticsContainer(): m_Median(std::numeric_limits::max()), m_Uniformity(std::numeric_limits::max()), m_UPP(std::numeric_limits::max()), m_Entropy(std::numeric_limits::max()) { } ImageStatisticsCalculator_v2::statisticsMapType ImageStatisticsCalculator_v2::HistogramStatisticsContainer::GetStatisticsAsMap() { ImageStatisticsCalculator_v2::statisticsMapType statisticsAsMap; statisticsAsMap["Median"] = m_Median; statisticsAsMap["Uniformity"] = m_Uniformity; statisticsAsMap["UPP"] = m_UPP; statisticsAsMap["Entropy"] = m_Entropy; return statisticsAsMap; } void ImageStatisticsCalculator_v2::HistogramStatisticsContainer::Reset() { m_Median = std::numeric_limits::max(); m_Uniformity = std::numeric_limits::max(); m_UPP = std::numeric_limits::max(); m_Entropy = std::numeric_limits::max(); } } diff --git a/Modules/ImageStatistics/mitkImageStatisticsCalculator2.h b/Modules/ImageStatistics/mitkImageStatisticsCalculator2.h index 80e2607e42..17522a9d17 100644 --- a/Modules/ImageStatistics/mitkImageStatisticsCalculator2.h +++ b/Modules/ImageStatistics/mitkImageStatisticsCalculator2.h @@ -1,352 +1,367 @@ #ifndef MITKIMAGESTATISTICSCALCULATOR2 #define MITKIMAGESTATISTICSCALCULATOR2 #include #include #include #include #include #include #include - /*#define mitkGetRefConstMacro(name, type) \ const type& Get##name() const \ { \ return &m_##name; \ } \ \ type& Get##name() \ { \ return &m_##name; \ } \ */ namespace mitk { - class MITKIMAGESTATISTICS_EXPORT ImageStatisticsCalculator_v2 + class MITKIMAGESTATISTICS_EXPORT ImageStatisticsCalculator_v2: public itk::Object { public: + /** Standard Self typedef */ + typedef ImageStatisticsCalculator_v2 Self; + typedef itk::Object Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self) + + /** Runtime information support. */ + itkTypeMacro(ImageStatisticsCalculator_v2, itk::Object) + typedef double statisticsValueType; typedef std::map statisticsMapType; typedef itk::Statistics::Histogram HistogramType; void SetInputImage(mitk::Image::Pointer image); void SetMask(mitk::MaskGenerator::Pointer mask); void SetDoHistogramStatistics(bool doHistStatistics); bool GetDoHistogramStatistics() const; void SetNBinsForHistogramStatistics(unsigned int nBins); unsigned int GetNBinsForHistogramStatistics() const; statisticsMapType GetStatistics(unsigned int timeStep=0, unsigned int label=0); // not yet implemented class StatisticsContainer: public itk::Object { public: /** Standard Self typedef */ typedef StatisticsContainer Self; typedef itk::Object Superclass; typedef itk::SmartPointer< Self > Pointer; typedef itk::SmartPointer< const Self > ConstPointer; /** Method for creation through the object factory. */ itkNewMacro(Self) /** Runtime information support. */ itkTypeMacro(StatisticsContainer, itk::Object) typedef double RealType; statisticsMapType GetStatisticsAsMap(); void Reset(); long & GetN() { return m_N; } const long& GetN() const { return m_N; } RealType & GetMean() { return m_Mean; } const RealType& GetMean() const { return m_Mean; } RealType & GetVariance() { return m_Variance; } const RealType& GetVariance() const { return m_Variance; } RealType & GetStd() { return m_Std; } const RealType& GetStd() const { return m_Std; } RealType & GetMin() { return m_Min; } const RealType& GetMin() const { return m_Min; } RealType & GetMax() { return m_Max; } const RealType& GetMax() const { return m_Max; } RealType & GetRMS() { return m_RMS; } const RealType& GetRMS() const { return m_RMS; } RealType & GetSkewness() { return m_Skewness; } const RealType& GetSkewness() const { return m_Skewness; } RealType & GetKurtosis() { return m_Kurtosis; } const RealType& GetKurtosis() const { return m_Kurtosis; } RealType & GetMPP() { return m_MPP; } const RealType& GetMPP() const { return m_MPP; } unsigned int & GetLabel() { return m_Label; } const unsigned int& GetLabel() const { return m_Label; } protected: StatisticsContainer(); private: // not pretty, is temporary - RealType m_Mean, m_Variance, m_Std, m_Min, m_Max, m_RMS; + 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; - long m_N; unsigned int m_Label; }; class HistogramStatisticsContainer: public itk::Object { public: /** Standard Self typedef */ typedef HistogramStatisticsContainer 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(HistogramStatisticsContainer, itk::Object) typedef double RealType; statisticsMapType GetStatisticsAsMap(); void Reset(); RealType & GetMedian() { return m_Median; } const RealType& GetMedian() const { return m_Median; } RealType & GetUniformity() { return m_Uniformity; } const RealType& GetUniformity() const { return m_Uniformity; } RealType & GetEntropy() { return m_Entropy; } const RealType& GetEntropy() const { return m_Entropy; } RealType & GetUPP() { return m_UPP; } const RealType& GetUPP() const { return m_UPP; } unsigned int & GetLabel() { return m_Label; } const unsigned int& GetLabel() const { return m_Label; } protected: HistogramStatisticsContainer(); private: RealType m_Median; RealType m_Uniformity; - RealType m_Entropy; RealType m_UPP; + RealType m_Entropy; unsigned int m_Label; }; struct HistogramContainer { unsigned int label; itk::Statistics::Histogram::Pointer histogram; }; protected: + ImageStatisticsCalculator_v2(){ + m_nBinsForHistogramStatistics = 100; + }; private: void SetAllStatisticsToUpdateRequired(); void SetAllHistogramStatisticsToUpdateRequired(); void SetStatsTimeStepToUpdateRequired(unsigned int timeStep); void SetHistStatsTimeStepToUpdateRequired(unsigned int timeStep); void CalculateHistogramStatistics(itk::Statistics::Histogram histogram, HistogramStatisticsContainer* HistStatContainer); template < typename TPixel, unsigned int VImageDimension > void InternalCheckMaskSanity( typename itk::Image< TPixel, VImageDimension >* image, typename itk::Image< unsigned short, VImageDimension >::Pointer maskImage); template < typename TPixel, unsigned int VImageDimension > void InternalCalculateStatisticsUnmasked( typename itk::Image< TPixel, VImageDimension >* image, unsigned int timeStep); HistogramStatisticsContainer::Pointer InternalCalculateHistogramStatistics(itk::Statistics::Histogram::Pointer histogram); template < typename TPixel, unsigned int VImageDimension > void InternalCalculateHistogramUnmasked( typename itk::Image< TPixel, VImageDimension >* image, double minVal, double maxVal, unsigned int timeStep); template < typename TPixel, unsigned int VImageDimension > void InternalCalculateStatisticsMasked( typename itk::Image< TPixel, VImageDimension >* image, unsigned int timeStep); template < typename TPixel, unsigned int VImageDimension > void ExtractMaskImageRegion( typename itk::Image< TPixel, VImageDimension >* image, typename itk::Image< unsigned short, VImageDimension >::Pointer maskImage, typename itk::Image< TPixel, VImageDimension >::Pointer adaptedImage ); template < typename TPixel, unsigned int VImageDimension > void GetMinAndMaxValue( typename itk::Image< TPixel, VImageDimension >* inputImage, double &minVal, double &maxVal ); std::string GetNameOfClass() { return std::string("ImageStatisticsCalculator_v2"); } mitk::Image::Pointer m_Image; mitk::Image::Pointer m_ImageTimeSlice; mitk::MaskGenerator::Pointer m_MaskGenerator; mitk::Image::Pointer m_InternalMask; bool m_DoHistogramStatistics; unsigned int m_nBinsForHistogramStatistics; std::vector m_StatisticsUpdateRequiredByTimeStep; // holds which time steps are valid and which ones have to be recalculated std::vector> m_StatisticsByTimeStep; std::vector> m_HistogramStatisticsByTimeStep; std::vector m_HistogramStatisticsUpdateRequiredByTimeStep; std::vector> m_HistogramsByTimeStep; }; } #endif // MITKIMAGESTATISTICSCALCULATOR2