diff --git a/Modules/DiffusionImaging/MiniApps/CMakeLists.txt b/Modules/DiffusionImaging/MiniApps/CMakeLists.txt index 21e436f79b..52c0895673 100755 --- a/Modules/DiffusionImaging/MiniApps/CMakeLists.txt +++ b/Modules/DiffusionImaging/MiniApps/CMakeLists.txt @@ -1,121 +1,122 @@ 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 + ImageStatisticsMiniapp_3^^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_3.cpp b/Modules/DiffusionImaging/MiniApps/ImageStatisticsMiniapp_3.cpp new file mode 100644 index 0000000000..9ce6f8a173 --- /dev/null +++ b/Modules/DiffusionImaging/MiniApps/ImageStatisticsMiniapp_3.cpp @@ -0,0 +1,82 @@ +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center, +Division of Medical and Biological Informatics. +All rights reserved. + +This software is distributed WITHOUT ANY WARRANTY; without +even the implied warranty of MERCHANTABILITY or FITNESS FOR +A PARTICULAR PURPOSE. + +See LICENSE.txt or http://www.mitk.org for details. + +===================================================================*/ + +#include "mitkCommandLineParser.h" +#include "mitkImage.h" +#include +#include +#include +#include +#include +#include "mitkIOUtil.h" +#include +#include +#include +#include +#include "mitkImageAccessByItk.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +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; +} + + + +int main( int argc, char* argv[] ) +{ + unsigned int timeStep = 0; + std::string inputImageFile; + inputImageFile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic2DplusT.nrrd"; + // Load image + 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::StatisticsContainer::Pointer result; + + calculator->SetInputImage(inputImage); + calculator->SetNBinsForHistogramStatistics(100); + + for (unsigned int i=0; i < inputImage->GetTimeSteps(); i++) + { + std::cout << "Results for time step " << i << ":" << std::endl; + result = calculator->GetStatistics(i, 1); + result->PrintSelf(); + std::cout << std::endl; + } + + + return EXIT_SUCCESS; +} diff --git a/Modules/DiffusionImaging/MiniApps/ImageStatisticsMiniapp_v2.cpp b/Modules/DiffusionImaging/MiniApps/ImageStatisticsMiniapp_v2.cpp index 57b65c87b4..63395ef046 100644 --- a/Modules/DiffusionImaging/MiniApps/ImageStatisticsMiniapp_v2.cpp +++ b/Modules/DiffusionImaging/MiniApps/ImageStatisticsMiniapp_v2.cpp @@ -1,205 +1,193 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkCommandLineParser.h" #include "mitkImage.h" #include #include #include #include #include #include "mitkIOUtil.h" #include #include #include #include #include "mitkImageAccessByItk.h" #include #include #include #include #include #include #include #include #include #include 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; } std::cout << std::endl; } template < typename TPixel, unsigned int VImageDimension > void get_statistics_boost(itk::Image* itkImage, statistics_res& res){ typedef itk::Image ImageType; itk::ImageRegionConstIterator it(itkImage, itkImage->GetLargestPossibleRegion()); int ctr=0; double sum=0; boost::accumulators::accumulator_set> > acc; for (it.GoToBegin(); !it.IsAtEnd(); ++it) { acc(it.Get()); // currentPixel = it.Get(); // sum+=currentPixel; // ctr+=1; } // res.mean=(double)sum/(double)ctr; res.mean = boost::accumulators::mean(acc); res.variance = boost::accumulators::variance(acc); res.min = boost::accumulators::min(acc); res.max = boost::accumulators::max(acc); res.count = boost::accumulators::count(acc); res.moment = boost::accumulators::moment<2>(acc); std::cout << "sum: " << sum << " N: " << ctr << " mean: " << res.mean << std::endl; } int main( int argc, char* argv[] ) { unsigned int timeStep = 0; std::string inputImageFile; - inputImageFile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic3D.nrrd"; + inputImageFile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic2DplusT.nrrd"; // Load image mitk::Image::Pointer inputImage = mitk::IOUtil::LoadImage(inputImageFile); mitk::Image::Pointer maskImage; std::string maskImageFile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic3D_someSegmentation.nrrd"; maskImage = mitk::IOUtil::LoadImage(maskImageFile); mitk::ImageMaskGenerator::Pointer binaryImageMaskGen = mitk::ImageMaskGenerator::New(); binaryImageMaskGen->SetImageMask(maskImage); mitk::HotspotMaskGenerator::Pointer hotSpotMaskGen = mitk::HotspotMaskGenerator::New(); hotSpotMaskGen->SetImage(inputImage); hotSpotMaskGen->SetHotspotRadiusInMM(10); hotSpotMaskGen->SetTimeStep(0); - std::string pFfile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic3D_rectangle.pf"; + std::string pFfile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic2DplusT_ellipse.pf"; std::vector loadedObjects; loadedObjects = mitk::IOUtil::Load(pFfile); mitk::BaseData::Pointer pf = loadedObjects[0]; mitk::PlanarFigure::Pointer planarFigure = dynamic_cast(pf.GetPointer()); mitk::PlanarFigureMaskGenerator::Pointer planarFigMaskExtr = mitk::PlanarFigureMaskGenerator::New(); planarFigMaskExtr->SetPlanarFigure(planarFigure); planarFigMaskExtr->SetImage(inputImage); // Calculate statistics mitk::ImageStatisticsCalculator_v2::Pointer calculator = mitk::ImageStatisticsCalculator_v2::New(); std::cout << "calculating statistics (unmasked) itk: " << std::endl; - mitk::ImageStatisticsCalculator_v2::statisticsMapType result; + mitk::ImageStatisticsCalculator_v2::StatisticsContainer::Pointer result; calculator->SetInputImage(inputImage); - calculator->SetMask(planarFigMaskExtr.GetPointer()); - calculator->SetDoHistogramStatistics(true); calculator->SetNBinsForHistogramStatistics(5); + std::cout << std::endl; + + std::cout << "unmasked: " << std::endl; + calculator->SetMask(nullptr); + result = calculator->GetStatistics(timeStep); + result->PrintSelf(); + std::cout << std::endl; std::cout << "planarfigure: " << std::endl; - try - { -// result = calculator->GetStatistics(timeStep, 1); - } - catch( const itk::ExceptionObject& e) - { - MITK_ERROR << "Statistic Calculation Failed - ITK Exception:" << e.what(); - return -1; - } + calculator->SetMask(planarFigMaskExtr.GetPointer()); + result = calculator->GetStatistics(timeStep, 1); + result->PrintSelf(); + std::cout << std::endl; -// printMap(result); - std::cout << "image mask: " << std::endl; - calculator->SetMask(binaryImageMaskGen.GetPointer()); - result = calculator->GetStatistics(timeStep, 1); - printMap(result); + + +// std::cout << "image mask: " << std::endl; +// calculator->SetMask(binaryImageMaskGen.GetPointer()); +// result = calculator->GetStatistics(timeStep, 1); +// result->PrintSelf(); // std::cout << "no mask: " << std::endl; // calculator->SetMask(nullptr); // result = calculator->GetStatistics(timeStep); -// printMap(result); +// result->PrintSelf(); // calculator->SetMask(hotSpotMaskGen.GetPointer()); // result = calculator->GetStatistics(timeStep, 1); -// printMap(result); +// result->PrintSelf(); // std::cout << std::endl << std::endl << "calculating statistics old imgstatcalc: " << std::endl; // mitk::ImageStatisticsCalculator::Statistics statisticsStruct; // mitk::ImageStatisticsCalculator::Pointer calculator_old = mitk::ImageStatisticsCalculator::New(); // calculator_old->SetImage(inputImage); // calculator_old->SetImageMask(maskImage); // //calculator_old->SetMaskingModeToImage(); // calculator_old->SetHistogramBinSize(100); // calculator_old->SetPlanarFigure(planarFigure); // calculator_old->SetMaskingModeToPlanarFigure(); // calculator_old->SetMaskingModeToNone(); // calculator_old->SetDoIgnorePixelValue(false); // calculator_old->ComputeStatistics(timeStep); // statisticsStruct = calculator_old->GetStatistics(timeStep); -// std::cout << "Entropy: " << statisticsStruct.GetEntropy() << std::endl; -// std::cout << "Kurtosis: " << statisticsStruct.GetKurtosis() << std::endl; -// std::cout << "MPP: " << statisticsStruct.GetMPP() << std::endl; -// std::cout << "Max: " << statisticsStruct.GetMax() << std::endl; -// std::cout << "Mean: " << statisticsStruct.GetMean() << std::endl; -// std::cout << "Median: " << statisticsStruct.GetMedian() << std::endl; -// std::cout << "Min: " << statisticsStruct.GetMin() << std::endl; -// std::cout << "N: " << statisticsStruct.GetN() << std::endl; -// std::cout << "RMS: " << statisticsStruct.GetRMS() << std::endl; -// std::cout << "Skewness: " << statisticsStruct.GetSkewness() << std::endl; -// std::cout << "Std: " << statisticsStruct.GetSigma() << std::endl; -// std::cout << "UPP: " << statisticsStruct.GetUPP() << std::endl; -// std::cout << "Uniformity: " << statisticsStruct.GetUniformity() << std::endl; -// std::cout << "Variance: " << statisticsStruct.GetVariance() << std::endl << std::endl; + // 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/files.cmake b/Modules/ImageStatistics/files.cmake index d76825a6f1..cefbe9ea3c 100644 --- a/Modules/ImageStatistics/files.cmake +++ b/Modules/ImageStatistics/files.cmake @@ -1,30 +1,31 @@ set(CPP_FILES mitkImageStatisticsCalculator.cpp mitkImageStatisticsCalculator2.cpp mitkPointSetStatisticsCalculator.cpp mitkPointSetDifferenceStatisticsCalculator.cpp mitkIntensityProfile.cpp mitkHotspotMaskGenerator.cpp mitkMaskGenerator.cpp mitkPlanarFigureMaskGenerator.cpp mitkMultiLabelMaskGenerator.cpp mitkImageMaskGenerator.cpp mitkHistogramStatisticsCalculator.cpp mitkMaskUtilities.cpp ) set(H_FILES mitkImageStatisticsCalculator.h mitkImageStatisticsCalculator2.h mitkPointSetDifferenceStatisticsCalculator.h mitkPointSetStatisticsCalculator.h mitkExtendedStatisticsImageFilter.h mitkExtendedLabelStatisticsImageFilter.h mitkHotspotMaskGenerator.h mitkMaskGenerator.h mitkPlanarFigureMaskGenerator.h mitkMultiLabelMaskGenerator.h mitkImageMaskGenerator.h mitkHistogramStatisticsCalculator.h mitkMaskUtilities.h + mitkitkMaskImageFilter.h ) diff --git a/Modules/ImageStatistics/mitkImageStatisticsCalculator2.cpp b/Modules/ImageStatistics/mitkImageStatisticsCalculator2.cpp index 3734371e12..99c69df1ab 100644 --- a/Modules/ImageStatistics/mitkImageStatisticsCalculator2.cpp +++ b/Modules/ImageStatistics/mitkImageStatisticsCalculator2.cpp @@ -1,531 +1,525 @@ #include #include #include #include #include #include #include #include +#include +#include +#include +#include #include #include #include #include #include #include #include +#include #include #include "itkImageFileWriter.h" namespace mitk { void ImageStatisticsCalculator_v2::SetInputImage(mitk::Image::Pointer image) { if (image != m_Image) { m_Image = image; - m_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(); + this->SetAllStatisticsToUpdateRequired(); } } unsigned int ImageStatisticsCalculator_v2::GetNBinsForHistogramStatistics() const { return m_nBinsForHistogramStatistics; } - ImageStatisticsCalculator_v2::statisticsMapType ImageStatisticsCalculator_v2::GetStatistics(unsigned int timeStep, unsigned int label) + ImageStatisticsCalculator_v2::StatisticsContainer::Pointer ImageStatisticsCalculator_v2::GetStatistics(unsigned int timeStep, unsigned int label) { std::cout << "getting statistics..." << std::endl; if (timeStep >= m_StatisticsByTimeStep.size()) { throw std::runtime_error("invalid timeStep in ImageStatisticsCalculator_v2::GetStatistics"); } if (m_Image.IsNull()) { throw std::runtime_error("no image"); } if (m_MaskGenerator.IsNotNull()) { m_MaskGenerator->SetTimeStep(timeStep); m_InternalMask = m_MaskGenerator->GetMask(); } + ImageTimeSelector::Pointer imgTimeSel = ImageTimeSelector::New(); + imgTimeSel->SetInput(m_Image); + imgTimeSel->SetTimeNr(timeStep); + imgTimeSel->UpdateLargestPossibleRegion(); + m_ImageTimeSlice = imgTimeSel->GetOutput(); + if (m_StatisticsUpdateRequiredByTimeStep[timeStep]) { std::cout << "calculating statistics..." << std::endl; // Calculate statistics with/without mask if (m_MaskGenerator.IsNull()) { // 1) calculate statistics unmasked: // plug image into itkstatisticsimagefilter (will be replaced by my awesome filter later on) // retrieve statistics and save them std::cout << "unmasked" << std::endl; - AccessByItk_1(m_Image, InternalCalculateStatisticsUnmasked, timeStep) + AccessByItk_1(m_ImageTimeSlice, InternalCalculateStatisticsUnmasked, timeStep) } else { std::cout << "masked" << std::endl; // 2) calculate statistics masked // extract mask image region // plug mask and image into itklabelstatisticsimagefilter - AccessByItk_1(m_Image, InternalCalculateStatisticsMasked, timeStep) - - } - m_StatisticsUpdateRequiredByTimeStep[timeStep]=false; - } - - - if(m_DoHistogramStatistics && m_HistogramStatisticsUpdateRequiredByTimeStep[timeStep]) - { - std::cout << "histogram statistics..." << std::endl; - - if (m_MaskGenerator.IsNull()) - { - std::cout << "unmasked" << std::endl; - - // 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; + AccessByItk_1(m_ImageTimeSlice, InternalCalculateStatisticsMasked, timeStep) } - else - { - std::cout << "masked" << std::endl; - m_HistogramStatisticsByTimeStep[timeStep].resize(0); - // 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_HistogramStatisticsByTimeStep[timeStep].push_back(histogramStatistics); - } - } - - - m_HistogramStatisticsUpdateRequiredByTimeStep[timeStep] = false; + m_StatisticsUpdateRequiredByTimeStep[timeStep]=false; } - // merge histogram and regular statistics, return result - statisticsMapType statistics; - if (m_MaskGenerator.IsNull()) + for (std::vector::iterator it = m_StatisticsByTimeStep[timeStep].begin(); it != m_StatisticsByTimeStep[timeStep].end(); ++it) { - statistics = m_StatisticsByTimeStep[timeStep][0]->GetStatisticsAsMap(); - if (m_DoHistogramStatistics) + StatisticsContainer::Pointer statCont = *it; + if (statCont->GetLabel() == label) { - statisticsMapType histStatistics = m_HistogramStatisticsByTimeStep[timeStep][0]->GetStatisticsAsMap(); - statistics.insert(histStatistics.begin(), histStatistics.end()); + return statCont; } } - 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; + MITK_ERROR << "Invalid Label: " << label; } 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::InternalCalculateStatisticsUnmasked( typename itk::Image< TPixel, VImageDimension >* image, unsigned int timeStep) { typedef typename itk::Image< TPixel, VImageDimension > ImageType; + StatisticsContainer::Pointer statisticsResult = StatisticsContainer::New(); + typename itk::ExtendedStatisticsImageFilter::Pointer statisticsFilter = itk::ExtendedStatisticsImageFilter::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 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); + // 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() = std::sqrt(std::pow(statisticsFilter->GetMean(), 2.) + statisticsFilter->GetVariance()); // variance = sigma^2 - statisticsResult->GetMPP() = statisticsFilter->GetMPP(); + 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()); + + + // histogram statistics + HistogramType::Pointer histogram = InternalCalculateHistogramUnmasked(image, statisticsResult->GetMin(), statisticsResult->GetMax()); + statisticsResult->SetHistogram(histogram); + + HistogramStatisticsCalculator histStatCalc; + histStatCalc.SetHistogram(histogram); + histStatCalc.CalculateStatistics(); + + statisticsResult->SetEntropy(histStatCalc.GetEntropy()); + statisticsResult->SetMedian(histStatCalc.GetMedian()); + statisticsResult->SetUniformity(histStatCalc.GetUniformity()); + statisticsResult->SetUPP(histStatCalc.GetUPP()); + m_StatisticsByTimeStep[timeStep][0] = statisticsResult; + } - template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator_v2::InternalCalculateHistogramUnmasked( + template < typename TPixel, unsigned int VImageDimension > ImageStatisticsCalculator_v2::HistogramType::Pointer ImageStatisticsCalculator_v2::InternalCalculateHistogramUnmasked( typename itk::Image< TPixel, VImageDimension >* image, double minVal, - double maxVal, - unsigned int timeStep) + double maxVal) { 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(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; + return histogram; } /* 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) { 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; typedef MaskUtilities< TPixel, VImageDimension > MaskUtilType; + typedef itk::MaskImageFilter2 MaskImageFilterType2; + typedef itk::MaskImageFilter MaskImageFilterType; + typedef itk::MinimumMaximumImageCalculator MinMaxCalculatorType; // maskImage has to have the same dimension as image typename MaskType::Pointer maskImage = MaskType::New(); maskImage = ImageToItkImage< unsigned short, VImageDimension >(m_InternalMask); typename MaskUtilType::Pointer maskUtil = MaskUtilType::New(); maskUtil->SetImage(image); maskUtil->SetMask(maskImage.GetPointer()); // if mask is smaller than image, extract the image region where the mask is typename ImageType::Pointer adaptedImage = ImageType::New(); adaptedImage = maskUtil->ExtractMaskImageRegion(); // this also checks mask sanity +// double minVal, maxVal; +// GetMinAndMaxValue(adaptedImage.GetPointer(), minVal, maxVal); + + /** Searching the min and max value is only done over the pixels that are nonzero in the mask. This bevahior differs from the 'old' + * imagestatisticscalculator where min and max were looked for over the entire adaptedImage*/ + // find min, max, minindex and maxindex + // make sure to only look in the masked region, use a masker for this + typename MaskImageFilterType::Pointer maskImageFilter = MaskImageFilterType::New(); + maskImageFilter->SetInput1(adaptedImage); + maskImageFilter->SetInput2(maskImage); + maskImageFilter->SetCoordinateTolerance(0.001); + maskImageFilter->SetDirectionTolerance(0.001); + maskImageFilter->Update(); + + typename MinMaxCalculatorType::Pointer imgMinMaxFilter = MinMaxCalculatorType::New(); + imgMinMaxFilter->SetImage(maskImageFilter->GetOutput()); + imgMinMaxFilter->Compute(); double minVal, maxVal; - GetMinAndMaxValue(adaptedImage.GetPointer(), minVal, maxVal); - - std::string outfile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/adapted_image.nrrd"; - typedef itk::ImageFileWriter< ImageType > WriterType; - typename WriterType::Pointer writer = WriterType::New(); - writer->SetFileName(outfile); - writer->SetInput(adaptedImage); - writer->Update(); + minVal = imgMinMaxFilter->GetMinimum(); + maxVal = imgMinMaxFilter->GetMaximum(); typename ImageStatisticsFilterType::Pointer imageStatisticsFilter = ImageStatisticsFilterType::New(); imageStatisticsFilter->SetDirectionTolerance(0.001); imageStatisticsFilter->SetCoordinateTolerance(0.001); imageStatisticsFilter->SetInput(adaptedImage); imageStatisticsFilter->SetLabelInput(maskImage); imageStatisticsFilter->SetHistogramParameters(m_nBinsForHistogramStatistics, floor(minVal), ceil(maxVal)); imageStatisticsFilter->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() = std::sqrt(std::pow(imageStatisticsFilter->GetMean(*it), 2.) + imageStatisticsFilter->GetVariance(*it)); // variance = sigma^2 - statisticsResult->GetMPP() = imageStatisticsFilter->GetMPP(*it); - statisticsResult->GetLabel() = *it; - m_StatisticsByTimeStep[timeStep].push_back(statisticsResult); - if (m_DoHistogramStatistics) + + // find min, max, minindex and maxindex + // make sure to only look in the masked region, use a masker for this + typename MaskImageFilterType2::Pointer maskImageFilter = MaskImageFilterType2::New(); + maskImageFilter->SetInput1(adaptedImage); + maskImageFilter->SetInput2(maskImage); + maskImageFilter->SetOutsideValue(0); + maskImageFilter->SetMaskingValue(*it); + double tmp_label = *it; + maskImageFilter->SetCoordinateTolerance(0.001); + maskImageFilter->SetDirectionTolerance(0.001); + maskImageFilter->UpdateLargestPossibleRegion(); + +// std::string outfile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic3D_masked_image.nrrd"; +// typedef itk::ImageFileWriter< ImageType > WriterType; +// typename WriterType::Pointer writer = WriterType::New(); +// writer->SetFileName(outfile); +// writer->SetInput(maskImageFilter->GetOutput()); +// writer->Update(); + + // TODO adapt minindex and maxindex if mask was 2D (careful: a mask can be 2d even if the itk::Image is 3D (3d image with a size of 1 in one dimension)) + typename MinMaxCalculatorType::Pointer imgMinMaxFilter = MinMaxCalculatorType::New(); + imgMinMaxFilter->SetImage(maskImageFilter->GetOutput()); + imgMinMaxFilter->Compute(); + vnl_vector minIndex, maxIndex; + 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++) { - HistogramContainer histCont; - histCont.label = *it; - histCont.histogram = imageStatisticsFilter->GetHistogram(*it); - m_HistogramsByTimeStep[timeStep].push_back(histCont); + minIndex[i] = tmpMinIndex[i] + (maskImage->GetOrigin()[i] - image->GetOrigin()[i]) / (double) maskImage->GetSpacing()[i]; + maxIndex[i] = tmpMaxIndex[i] + (maskImage->GetOrigin()[i] - image->GetOrigin()[i]) / (double) maskImage->GetSpacing()[i]; } + + statisticsResult->SetMinIndex(minIndex); + statisticsResult->SetMaxIndex(maxIndex); + + 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); + + HistogramType::Pointer histogram = imageStatisticsFilter->GetHistogram(*it); + statisticsResult->SetHistogram(histogram); + + HistogramStatisticsCalculator histStatCalc; + histStatCalc.SetHistogram(histogram); + histStatCalc.CalculateStatistics(); + + statisticsResult->SetEntropy(histStatCalc.GetEntropy()); + statisticsResult->SetMedian(histStatCalc.GetMedian()); + statisticsResult->SetUniformity(histStatCalc.GetUniformity()); + statisticsResult->SetUPP(histStatCalc.GetUPP()); + + m_StatisticsByTimeStep[timeStep].push_back(statisticsResult); ++it; } } - ImageStatisticsCalculator_v2::HistogramStatisticsContainer::Pointer ImageStatisticsCalculator_v2::InternalCalculateHistogramStatistics(itk::Statistics::Histogram::Pointer histogram) - { - HistogramStatisticsCalculator histStatCalc; - histStatCalc.SetHistogram(histogram); - histStatCalc.CalculateStatistics(); - HistogramStatisticsContainer::Pointer histStatResults = HistogramStatisticsContainer::New(); - histStatResults->GetEntropy() = histStatCalc.GetEntropy(); - histStatResults->GetMedian() = histStatCalc.GetMedian(); - histStatResults->GetUniformity() = histStatCalc.GetUniformity(); - histStatResults->GetUPP() = histStatCalc.GetUPP(); - return histStatResults; - } ImageStatisticsCalculator_v2::StatisticsContainer::StatisticsContainer(): m_N(std::numeric_limits::max()), m_Mean(std::numeric_limits::max()), - m_Min(std::numeric_limits::max()), + m_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_MPP(std::numeric_limits::max()), + m_Median(std::numeric_limits::max()), + m_Uniformity(std::numeric_limits::max()), + m_UPP(std::numeric_limits::max()), + m_Entropy(std::numeric_limits::max()) { - + m_minIndex.set_size(0); + m_maxIndex.set_size(0); } ImageStatisticsCalculator_v2::statisticsMapType ImageStatisticsCalculator_v2::StatisticsContainer::GetStatisticsAsMap() { ImageStatisticsCalculator_v2::statisticsMapType statisticsAsMap; statisticsAsMap["N"] = m_N; statisticsAsMap["Mean"] = m_Mean; statisticsAsMap["Min"] = m_Min; statisticsAsMap["Max"] = m_Max; statisticsAsMap["StandardDeviation"] = m_Std; statisticsAsMap["Variance"] = m_Variance; statisticsAsMap["Skewness"] = m_Skewness; statisticsAsMap["Kurtosis"] = m_Kurtosis; statisticsAsMap["RMS"] = m_RMS; statisticsAsMap["MPP"] = m_MPP; + statisticsAsMap["Median"] = m_Median; + statisticsAsMap["Uniformity"] = m_Uniformity; + statisticsAsMap["UPP"] = m_UPP; + statisticsAsMap["Entropy"] = m_Entropy; return statisticsAsMap; } void ImageStatisticsCalculator_v2::StatisticsContainer::Reset() { m_N = std::numeric_limits::max(); m_Mean = std::numeric_limits::max(); m_Min = std::numeric_limits::max(); m_Max = std::numeric_limits::max(); m_Std = std::numeric_limits::max(); m_Variance = std::numeric_limits::max(); m_Skewness = std::numeric_limits::max(); m_Kurtosis = std::numeric_limits::max(); m_RMS = std::numeric_limits::max(); m_MPP = std::numeric_limits::max(); + m_Median = std::numeric_limits::max(); + m_Uniformity = std::numeric_limits::max(); + m_UPP = std::numeric_limits::max(); + m_Entropy = std::numeric_limits::max(); } - - 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()) + void ImageStatisticsCalculator_v2::StatisticsContainer::PrintSelf() { + ImageStatisticsCalculator_v2::statisticsMapType statMap = this->GetStatisticsAsMap(); + // print all map key value pairs + for (auto it = statMap.begin(); it != statMap.end(); ++it) + { + std::cout << it->first << ": " << it->second << std::endl; + } - } - - 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; + // 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; - return statisticsAsMap; + // 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; } - - 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 8e1c5ffd37..9eecd3d6aa 100644 --- a/Modules/ImageStatistics/mitkImageStatisticsCalculator2.h +++ b/Modules/ImageStatistics/mitkImageStatisticsCalculator2.h @@ -1,357 +1,333 @@ #ifndef MITKIMAGESTATISTICSCALCULATOR2 #define MITKIMAGESTATISTICSCALCULATOR2 #include #include #include #include #include #include #include /*#define mitkGetRefConstMacro(name, type) \ const type& Get##name() const \ { \ return &m_##name; \ } \ \ type& Get##name() \ { \ return &m_##name; \ } \ */ namespace mitk { class MITKIMAGESTATISTICS_EXPORT ImageStatisticsCalculator_v2: public itk::Object { public: /** Standard Self typedef */ typedef ImageStatisticsCalculator_v2 Self; typedef itk::Object Superclass; typedef itk::SmartPointer< Self > Pointer; typedef itk::SmartPointer< const Self > ConstPointer; /** Method for creation through the object factory. */ itkNewMacro(Self) /** Runtime information support. */ itkTypeMacro(ImageStatisticsCalculator_v2, itk::Object) typedef double statisticsValueType; typedef std::map statisticsMapType; typedef itk::Statistics::Histogram HistogramType; - 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=1); - 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() + void SetN(long n) { - return m_N; + m_N = n; } const long& GetN() const { return m_N; } - RealType & GetMean() + void SetMean(RealType mean) { - return m_Mean; + m_Mean = mean; } const RealType& GetMean() const { return m_Mean; } - RealType & GetVariance() + void SetVariance(RealType variance) { - return m_Variance; + m_Variance = variance; } const RealType& GetVariance() const { return m_Variance; } - RealType & GetStd() + void SetStd(RealType std) { - return m_Std; + m_Std = std; } const RealType& GetStd() const { return m_Std; } - RealType & GetMin() + void SetMin(RealType minVal) { - return m_Min; + m_Min = minVal; } const RealType& GetMin() const { return m_Min; } - RealType & GetMax() + void SetMax(RealType maxVal) { - return m_Max; + m_Max = maxVal; } const RealType& GetMax() const { return m_Max; } - RealType & GetRMS() + void SetRMS(RealType rms) { - return m_RMS; + m_RMS = rms; } const RealType& GetRMS() const { return m_RMS; } - RealType & GetSkewness() + void SetSkewness(RealType skewness) { - return m_Skewness; + m_Skewness = skewness; } const RealType& GetSkewness() const { return m_Skewness; } - RealType & GetKurtosis() + void SetKurtosis(RealType kurtosis) { - return m_Kurtosis; + m_Kurtosis = kurtosis; } const RealType& GetKurtosis() const { return m_Kurtosis; } - RealType & GetMPP() + void SetMPP(RealType mpp) { - return m_MPP; + m_MPP = mpp; } const RealType& GetMPP() const { return m_MPP; } - unsigned int & GetLabel() + void SetLabel(unsigned int label) { - return m_Label; + m_Label = label; } const unsigned int& GetLabel() const { return m_Label; } - protected: - StatisticsContainer(); - - - private: - // not pretty, is temporary - long m_N; - RealType m_Mean, m_Min, m_Max, m_Std, m_Variance; - RealType m_Skewness; - RealType m_Kurtosis; - RealType m_RMS; - RealType m_MPP; - 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; - + void SetMinIndex(vnl_vector& minIndex) + { + m_minIndex = minIndex; + } - statisticsMapType GetStatisticsAsMap(); + const vnl_vector& GetMinIndex() const + { + return m_minIndex; + } - void Reset(); + void SetMaxIndex(vnl_vector& maxIndex) + { + m_maxIndex = maxIndex; + } - RealType & GetMedian() + const vnl_vector& GetMaxIndex() const { - return m_Median; + return m_maxIndex; } - const RealType& GetMedian() const + void SetHistogram(HistogramType::Pointer hist) { - return m_Median; + if (m_Histogram != hist) + { + m_Histogram = hist; + } } - RealType & GetUniformity() + const HistogramType::Pointer GetHistogram() const { - return m_Uniformity; + return m_Histogram; } - const RealType& GetUniformity() const + void SetEntropy(RealType entropy) { - return m_Uniformity; + m_Entropy = entropy; } - RealType & GetEntropy() + const RealType & GetEntropy() const { return m_Entropy; } - const RealType& GetEntropy() const + void SetMedian(RealType median) { - return m_Entropy; + m_Median = median; } - RealType & GetUPP() + const RealType & GetMedian() const { - return m_UPP; + return m_Median; } - const RealType& GetUPP() const + void SetUniformity(RealType uniformity) { - return m_UPP; + m_Uniformity = uniformity; } - unsigned int & GetLabel() + const RealType & GetUniformity() const { - return m_Label; + return m_Uniformity; } - const unsigned int& GetLabel() const + void SetUPP(RealType upp) { - return m_Label; + m_UPP = upp; } + const RealType & GetUPP() const + { + return m_UPP; + } + + void PrintSelf(); + protected: - HistogramStatisticsContainer(); + StatisticsContainer(); + private: + // not pretty, is temporary + long m_N; + RealType m_Mean, m_Min, m_Max, m_Std, m_Variance; + RealType m_Skewness; + RealType m_Kurtosis; + RealType m_RMS; + RealType m_MPP; + vnl_vector m_minIndex, m_maxIndex; RealType m_Median; RealType m_Uniformity; RealType m_UPP; RealType m_Entropy; unsigned int m_Label; + HistogramType::Pointer m_Histogram; }; - struct HistogramContainer - { - unsigned int label; - itk::Statistics::Histogram::Pointer histogram; - }; + + + void SetInputImage(mitk::Image::Pointer image); + + void SetMask(mitk::MaskGenerator::Pointer mask); + + void SetNBinsForHistogramStatistics(unsigned int nBins); + + unsigned int GetNBinsForHistogramStatistics() const; + + StatisticsContainer::Pointer GetStatistics(unsigned int timeStep=0, unsigned int label=1); protected: ImageStatisticsCalculator_v2(){ m_nBinsForHistogramStatistics = 100; }; private: void SetAllStatisticsToUpdateRequired(); - void 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 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( + template < typename TPixel, unsigned int VImageDimension > typename HistogramType::Pointer InternalCalculateHistogramUnmasked( typename itk::Image< TPixel, VImageDimension >* image, double minVal, - double maxVal, - unsigned int timeStep); + double maxVal); template < typename TPixel, unsigned int VImageDimension > void InternalCalculateStatisticsMasked( typename itk::Image< TPixel, VImageDimension >* image, unsigned int timeStep); template < typename TPixel, unsigned int VImageDimension > void GetMinAndMaxValue( typename itk::Image< TPixel, VImageDimension >* inputImage, double &minVal, double &maxVal ); std::string GetNameOfClass() { return std::string("ImageStatisticsCalculator_v2"); } mitk::Image::Pointer m_Image; mitk::Image::Pointer m_ImageTimeSlice; mitk::MaskGenerator::Pointer m_MaskGenerator; mitk::Image::Pointer m_InternalMask; - bool m_DoHistogramStatistics; unsigned int m_nBinsForHistogramStatistics; std::vector m_StatisticsUpdateRequiredByTimeStep; // holds which time steps are valid and which ones have to be recalculated std::vector> m_StatisticsByTimeStep; - - std::vector> m_HistogramStatisticsByTimeStep; - std::vector m_HistogramStatisticsUpdateRequiredByTimeStep; - - std::vector> m_HistogramsByTimeStep; }; } #endif // MITKIMAGESTATISTICSCALCULATOR2 diff --git a/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.cpp b/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.cpp index 2df2feb9e8..2c2554d8bd 100644 --- a/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.cpp +++ b/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.cpp @@ -1,388 +1,402 @@ #include #include #include #include "mitkImageAccessByItk.h" #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) { m_Modified = true; m_PlanarFigure = planarFigure; } } void PlanarFigureMaskGenerator::SetImage(mitk::Image::Pointer image) { // 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) { m_Modified = true; m_InternalInputImage = image; } } 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 itk::Image::Pointer maskImage = itk::Image::New(); // maskImage->SetDirection(image->GetDirection()); // maskImage->SetOrigin(image->GetOrigin()); // maskImage->SetSpacing(image->GetSpacing()); // maskImage->SetLargestPossibleRegion(image->GetLargestPossibleRegion()); // 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 ); // 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) + { + mitk::ImageTimeSelector::Pointer imgTimeSel = mitk::ImageTimeSelector::New(); + imgTimeSel->SetInput(m_InternalInputImage); + imgTimeSel->SetTimeNr(m_TimeStep); + imgTimeSel->UpdateLargestPossibleRegion(); + m_InternalTimeSliceImage = imgTimeSel->GetOutput(); + } + else + { + m_InternalTimeSliceImage = m_InternalInputImage; + } + m_InternalITKImageMask2D = nullptr; const PlaneGeometry *planarFigurePlaneGeometry = m_PlanarFigure->GetPlaneGeometry(); const PlaneGeometry *planarFigureGeometry = dynamic_cast< const PlaneGeometry * >( planarFigurePlaneGeometry ); const BaseGeometry *imageGeometry = m_InternalInputImage->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 mitk::Image::Pointer inputImageSlice = extract2DImageSlice(axis, slice); // 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()); Convert2Dto3DImageFilter::Pointer sliceTo3DImageConverter = Convert2Dto3DImageFilter::New(); sliceTo3DImageConverter->SetInput(planarFigureMaskImage); sliceTo3DImageConverter->Update(); m_InternalMask = sliceTo3DImageConverter->GetOutput(); } mitk::Image::Pointer PlanarFigureMaskGenerator::GetMask() { if (m_Modified) { 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(); } m_Modified = false; 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_InternalInputImage->GetDimension(); + unsigned int dimension = m_InternalTimeSliceImage->GetDimension(); mitk::Image::Pointer imageSlice = mitk::Image::New(); if (dimension == 3) { ExtractImageFilter::Pointer imageExtractor = ExtractImageFilter::New(); - imageExtractor->SetInput( m_InternalInputImage ); + imageExtractor->SetInput( m_InternalTimeSliceImage ); imageExtractor->SetSliceDimension( axis ); imageExtractor->SetSliceIndex( slice ); imageExtractor->Update(); imageSlice = imageExtractor->GetOutput(); } else if(dimension == 2) { - imageSlice = m_InternalInputImage; + imageSlice = m_InternalTimeSliceImage; } else { MITK_ERROR << "Unsupported image dimension. Dimension is: " << dimension << ". Only 2D and 3D images are supported."; } return imageSlice; } } diff --git a/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.h b/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.h index d5eb54f483..48e56b7aa1 100644 --- a/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.h +++ b/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.h @@ -1,106 +1,107 @@ #ifndef MITKPLANARFIGUREMASKGENERATOR #define MITKPLANARFIGUREMASKGENERATOR #include #include #include #include #include #include #include #include #include #include namespace mitk { 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) mitk::Image::Pointer GetMask(); void SetPlanarFigure(mitk::PlanarFigure::Pointer planarFigure); /*We need an image because we need its geometry in order to give the points of the planar figure the correct world coordinates*/ void SetImage(const mitk::Image::Pointer inputImage); protected: PlanarFigureMaskGenerator(){} 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()); } mitk::PlanarFigure::Pointer m_PlanarFigure; typename itk::Image::Pointer m_InternalITKImageMask2D; mitk::Image::Pointer m_InternalInputImage; + mitk::Image::Pointer m_InternalTimeSliceImage; unsigned int m_PlanarFigureAxis; }; } #endif // MITKPLANARFIGUREMASKGENERATOR diff --git a/Modules/ImageStatistics/mitkitkMaskImageFilter.h b/Modules/ImageStatistics/mitkitkMaskImageFilter.h new file mode 100644 index 0000000000..13a064818a --- /dev/null +++ b/Modules/ImageStatistics/mitkitkMaskImageFilter.h @@ -0,0 +1,290 @@ +/*========================================================================= + * + * Copyright Insight Software Consortium + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#ifndef __itkMaskImageFilter2_h +#define __itkMaskImageFilter2_h + +#include "itkBinaryFunctorImageFilter.h" +#include "itkNumericTraits.h" +#include "itkVariableLengthVector.h" +#include + +namespace itk +{ +namespace Functor +{ +/** + * \class MaskInput + * \brief + * \ingroup ITKImageIntensity + */ +template< typename TInput, typename TMask, typename TOutput = TInput > +class MaskInput2 +{ +public: + typedef typename NumericTraits< TInput >::AccumulateType AccumulatorType; + + MaskInput2() + { + m_MaskingValue = NumericTraits< TMask >::ZeroValue(); + InitializeOutsideValue( static_cast( ITK_NULLPTR ) ); + } + ~MaskInput2() {} + bool operator!=(const MaskInput2 &) const + { + return false; + } + + bool operator==(const MaskInput2 & other) const + { + return !( *this != other ); + } + + inline TOutput operator()(const TInput & A, const TMask & B) const + { + if ( B == m_MaskingValue ) + { + return static_cast< TOutput >( A ); + } + else + { + return m_OutsideValue; + } + } + + /** Method to explicitly set the outside value of the mask */ + void SetOutsideValue(const TOutput & outsideValue) + { + m_OutsideValue = outsideValue; + } + + /** Method to get the outside value of the mask */ + const TOutput & GetOutsideValue() const + { + return m_OutsideValue; + } + + /** Method to explicitly set the masking value */ + void SetMaskingValue(const TMask & maskingValue) + { + m_MaskingValue = maskingValue; + } + /** Method to get the masking value */ + const TMask & GetMaskingValue() const + { + return m_MaskingValue; + } + +private: + + template < typename TPixelType > + void InitializeOutsideValue( TPixelType * ) + { + this->m_OutsideValue = NumericTraits< TPixelType >::ZeroValue(); + } + + template < typename TValue > + void InitializeOutsideValue( VariableLengthVector * ) + { + // set the outside value to be of zero length + this->m_OutsideValue = VariableLengthVector< TValue >(0); + } + + TOutput m_OutsideValue; + TMask m_MaskingValue; +}; +} +/** \class MaskImageFilter + * \brief Mask an image with a mask. + * + * This class is templated over the types of the + * input image type, the mask image type and the type of the output image. + * Numeric conversions (castings) are done by the C++ defaults. + * + * The pixel type of the input 2 image must have a valid definition of the + * operator != with zero. This condition is required because internally this + * filter will perform the operation + * + * \code + * if pixel_from_mask_image != masking_value + * pixel_output_image = pixel_input_image + * else + * pixel_output_image = outside_value + * \endcode + * + * The pixel from the input 1 is cast to the pixel type of the output image. + * + * Note that the input and the mask images must be of the same size. + * + * \warning Any pixel value other than masking value (0 by default) will not be masked out. + * + * \sa MaskNegatedImageFilter + * \ingroup IntensityImageFilters + * \ingroup MultiThreaded + * \ingroup ITKImageIntensity + * + * \wiki + * \wikiexample{ImageProcessing/MaskImageFilter,Apply a mask to an image} + * \endwiki + */ +template< typename TInputImage, typename TMaskImage, typename TOutputImage = TInputImage > +class MITKIMAGESTATISTICS_EXPORT MaskImageFilter2: + public + BinaryFunctorImageFilter< TInputImage, TMaskImage, TOutputImage, + Functor::MaskInput2< + typename TInputImage::PixelType, + typename TMaskImage::PixelType, + typename TOutputImage::PixelType > > + +{ +public: + /** Standard class typedefs. */ + typedef MaskImageFilter2 Self; + typedef BinaryFunctorImageFilter< TInputImage, TMaskImage, TOutputImage, + Functor::MaskInput2< + typename TInputImage::PixelType, + typename TMaskImage::PixelType, + typename TOutputImage::PixelType > + > Superclass; + + typedef SmartPointer< Self > Pointer; + typedef SmartPointer< const Self > ConstPointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Runtime information support. */ + itkTypeMacro(MaskImageFilter2, + BinaryFunctorImageFilter); + + /** Typedefs **/ + typedef TMaskImage MaskImageType; + + /** Set/Get the mask image. Pixels set in the mask image will retain + * the original value of the input image while pixels not set in + * the mask will be set to the "OutsideValue". + */ + void SetMaskImage(const MaskImageType *maskImage) + { + // Process object is not const-correct so the const casting is required. + this->SetNthInput( 1, const_cast< MaskImageType * >( maskImage ) ); + } + const MaskImageType * GetMaskImage() + { + return static_cast(this->ProcessObject::GetInput(1)); + } + + /** Method to explicitly set the outside value of the mask. Defaults to 0 */ + void SetOutsideValue(const typename TOutputImage::PixelType & outsideValue) + { + if ( this->GetOutsideValue() != outsideValue ) + { + this->Modified(); + this->GetFunctor().SetOutsideValue(outsideValue); + } + } + + const typename TOutputImage::PixelType & GetOutsideValue() const + { + return this->GetFunctor().GetOutsideValue(); + } + + /** Method to explicitly set the masking value of the mask. Defaults to 0 */ + void SetMaskingValue(const typename TMaskImage::PixelType & maskingValue) + { + if ( this->GetMaskingValue() != maskingValue ) + { + this->Modified(); + this->GetFunctor().SetMaskingValue(maskingValue); + } + } + + /** Method to get the masking value of the mask. */ + const typename TMaskImage::PixelType & GetMaskingValue() const + { + return this->GetFunctor().GetMaskingValue(); + } + + void BeforeThreadedGenerateData() + { + typedef typename TOutputImage::PixelType PixelType; + this->CheckOutsideValue( static_cast(ITK_NULLPTR) ); + } + +#ifdef ITK_USE_CONCEPT_CHECKING + // Begin concept checking + itkConceptMacro( MaskEqualityComparableCheck, + ( Concept::EqualityComparable< typename TMaskImage::PixelType > ) ); + itkConceptMacro( InputConvertibleToOutputCheck, + ( Concept::Convertible< typename TInputImage::PixelType, + typename TOutputImage::PixelType > ) ); + // End concept checking +#endif + +protected: + MaskImageFilter2() {} + virtual ~MaskImageFilter2() {} + + void PrintSelf(std::ostream & os, Indent indent) const + { + Superclass::PrintSelf(os, indent); + os << indent << "OutsideValue: " << this->GetOutsideValue() << std::endl; + } + +private: + MaskImageFilter2(const Self &); //purposely not implemented + void operator=(const Self &); //purposely not implemented + + template < typename TPixelType > + void CheckOutsideValue( const TPixelType * ) {} + + template < typename TValue > + void CheckOutsideValue( const VariableLengthVector< TValue > * ) + { + // Check to see if the outside value contains only zeros. If so, + // resize it to have the same number of zeros as the output + // image. Otherwise, check that the number of components in the + // outside value is the same as the number of components in the + // output image. If not, throw an exception. + VariableLengthVector< TValue > currentValue = + this->GetFunctor().GetOutsideValue(); + VariableLengthVector< TValue > zeroVector( currentValue.GetSize() ); + zeroVector.Fill( NumericTraits< TValue >::ZeroValue() ); + + if ( currentValue == zeroVector ) + { + zeroVector.SetSize( this->GetOutput()->GetVectorLength() ); + zeroVector.Fill( NumericTraits< TValue >::ZeroValue() ); + this->GetFunctor().SetOutsideValue( zeroVector ); + } + else if ( this->GetFunctor().GetOutsideValue().GetSize() != + this->GetOutput()->GetVectorLength() ) + { + itkExceptionMacro( + << "Number of components in OutsideValue: " + << this->GetFunctor().GetOutsideValue().GetSize() + << " is not the same as the " + << "number of components in the image: " + << this->GetOutput()->GetVectorLength()); + } + } + +}; +} // end namespace itk + +#endif +