diff --git a/Modules/DiffusionImaging/MiniApps/CMakeLists.txt b/Modules/DiffusionImaging/MiniApps/CMakeLists.txt index 91ee4ac0a0..e4564aff4f 100755 --- a/Modules/DiffusionImaging/MiniApps/CMakeLists.txt +++ b/Modules/DiffusionImaging/MiniApps/CMakeLists.txt @@ -1,118 +1,120 @@ 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 + 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/HotspotMiniApp.cpp b/Modules/DiffusionImaging/MiniApps/HotspotMiniApp.cpp new file mode 100644 index 0000000000..2f76e385e9 --- /dev/null +++ b/Modules/DiffusionImaging/MiniApps/HotspotMiniApp.cpp @@ -0,0 +1,104 @@ +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center, +Division of Medical and Biological Informatics. +All rights reserved. + +This software is distributed WITHOUT ANY WARRANTY; without +even the implied warranty of MERCHANTABILITY or FITNESS FOR +A PARTICULAR PURPOSE. + +See LICENSE.txt or http://www.mitk.org for details. + +===================================================================*/ + +#include "mitkCommandLineParser.h" +#include "mitkImage.h" +#include "mitkIOUtil.h" +#include +#include +#include +#include + + +int main( int argc, char* argv[] ) +{ + mitkCommandLineParser parser; + + parser.setTitle("Test basic hotspot functionality"); + parser.setCategory("Preprocessing Tools"); + parser.setDescription(""); + parser.setContributor("MBI"); + + parser.setArgumentPrefix("--", "-"); + parser.addArgument("help", "h", mitkCommandLineParser::String, "Help:", "Show this help text"); + parser.addArgument("input", "i", mitkCommandLineParser::InputFile, "Input:", "input image", us::Any(),false); + parser.addArgument("mask", "m", mitkCommandLineParser::InputFile, "Mask:", "mask image / roi image denotin area on which statistics are calculated", us::Any(),true); + parser.addArgument("out", "o", mitkCommandLineParser::OutputFile, "Output", "output file (default: hotspotMiniApp_output.nrrd)", us::Any()); + + std::cout << "test...." << std::endl; + + std::map parsedArgs = parser.parseArguments(argc, argv); + std::cout << "parsedArgs.size()= " << parsedArgs.size() << std::endl; + if (parsedArgs.size()==0 || parsedArgs.count("help") || parsedArgs.count("h")) + { +// std::cout << "\n\n MiniApp Description: \nCalculates and saves the hotspot of an image." << endl; +// std::cout << "Output is written to the designated output file" << endl; +// std::cout << parser.helpText(); +// return EXIT_SUCCESS; + } + + // Parameters: + std::string inputImageFile; + if (!parsedArgs.count("i") && !parsedArgs.count("input")) + { + inputImageFile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic3D.nrrd"; + } + else + { + inputImageFile = us::any_cast(parsedArgs["input"]); + } + + mitk::Image::Pointer maskImage; + if (parsedArgs.count("mask") || parsedArgs.count("m")) + { + std::string maskImageFile = us::any_cast(parsedArgs["mask"]); + maskImage = mitk::IOUtil::LoadImage(maskImageFile); + } + + std::string outFile; + if (parsedArgs.count("out") || parsedArgs.count("o") ) + outFile = us::any_cast(parsedArgs["out"]); + else + outFile = "hotspotMiniApp_output.nrrd"; + + // Load image and mask + mitk::Image::Pointer inputImage = mitk::IOUtil::LoadImage(inputImageFile); + + // Calculate statistics + mitk::HotspotCalculator hotspotCalc; + hotspotCalc.setInputImage(inputImage); + hotspotCalc.setMaskImage(maskImage); + hotspotCalc.setHotspotRadiusInMM(10); + + mitk::Image::Pointer outImage; + try + { + outImage = hotspotCalc.getHotspotMask(0, 0); + } + catch( const itk::ExceptionObject& e) + { + MITK_ERROR << "Failed - ITK Exception:" << e.what(); + return -1; + } + + if (outImage != nullptr) + { + mitk::IOUtil::SaveImage(outImage, outFile); + } + + + return EXIT_SUCCESS; +} diff --git a/Modules/DiffusionImaging/MiniApps/PlanarFigureMaskCreatorMiniapp.cpp b/Modules/DiffusionImaging/MiniApps/PlanarFigureMaskCreatorMiniapp.cpp new file mode 100644 index 0000000000..c099bdf1e9 --- /dev/null +++ b/Modules/DiffusionImaging/MiniApps/PlanarFigureMaskCreatorMiniapp.cpp @@ -0,0 +1,131 @@ +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center, +Division of Medical and Biological Informatics. +All rights reserved. + +This software is distributed WITHOUT ANY WARRANTY; without +even the implied warranty of MERCHANTABILITY or FITNESS FOR +A PARTICULAR PURPOSE. + +See LICENSE.txt or http://www.mitk.org for details. + +===================================================================*/ + +#include "mitkCommandLineParser.h" +#include "mitkImage.h" +#include "mitkIOUtil.h" +#include +#include +#include +#include +#include +#include + + +int main( int argc, char* argv[] ) +{ + mitkCommandLineParser parser; + + parser.setTitle("Test basic hotspot functionality"); + parser.setCategory("Preprocessing Tools"); + parser.setDescription(""); + parser.setContributor("MBI"); + + parser.setArgumentPrefix("--", "-"); + parser.addArgument("help", "h", mitkCommandLineParser::String, "Help:", "Show this help text"); + parser.addArgument("input", "i", mitkCommandLineParser::InputFile, "Input:", "input image", us::Any(),false); + parser.addArgument("planarfigure", "p", mitkCommandLineParser::InputFile, "PlanarFigure:", "mask image / roi image denotin area on which statistics are calculated", us::Any(),false); + parser.addArgument("out", "o", mitkCommandLineParser::OutputFile, "Output", "output file (default: hotspotMiniApp_output.nrrd)", us::Any()); + + std::cout << "test...." << std::endl; + + std::map parsedArgs = parser.parseArguments(argc, argv); + std::cout << "parsedArgs.size()= " << parsedArgs.size() << std::endl; + if (parsedArgs.size()==0 || parsedArgs.count("help") || parsedArgs.count("h")) + { +// std::cout << "\n\n MiniApp Description: \nCalculates and saves the hotspot of an image." << endl; +// std::cout << "Output is written to the designated output file" << endl; +// std::cout << parser.helpText(); +// return EXIT_SUCCESS; + } + + // Parameters: + std::string inputImageFile; + if (!parsedArgs.count("i") && !parsedArgs.count("input")) + { + inputImageFile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic3D.nrrd"; + } + else + { + inputImageFile = us::any_cast(parsedArgs["input"]); + } + + std::string inputPFFile; + if (parsedArgs.count("planarfigure") || parsedArgs.count("p")) + { + inputPFFile = us::any_cast(parsedArgs["planarfigure"]); + } + else + { + inputPFFile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic3D_rectangle.pf"; + } + + + mitk::PlanarFigure::Pointer planarFigure; + try + { + std::vector loadedObjects; + // try except + loadedObjects = mitk::IOUtil::Load(inputPFFile); + mitk::BaseData::Pointer pf = loadedObjects[0]; + planarFigure = dynamic_cast(pf.GetPointer()); + if (planarFigure.IsNull()) + { + MITK_ERROR << "something went wrong"; + return -1; + } + } + catch (const itk::ExceptionObject& e) + { + MITK_ERROR << "Failed: " << e.what(); + return -1; + } + + std::string outFile; + if (parsedArgs.count("out") || parsedArgs.count("o") ) + outFile = us::any_cast(parsedArgs["out"]); + else + outFile = "planarFigureExtraction_output.nrrd"; + + // Load image and mask + mitk::Image::Pointer inputImage = mitk::IOUtil::LoadImage(inputImageFile); + + // Calculate statistics + mitk::PlanarFigureMaskGenerator planarFigMaskExtr; + planarFigMaskExtr.SetPlanarFigure(planarFigure); + planarFigMaskExtr.SetImage(inputImage); + + mitk::Image::Pointer outImage; + try + { + outImage = planarFigMaskExtr.GetMask(); + } + catch( const itk::ExceptionObject& e) + { + MITK_ERROR << "Failed - ITK Exception:" << e.what(); + return -1; + } + + mitk::BaseGeometry *imageGeo = outImage->GetGeometry(); + std::cout << "origin: " << imageGeo->GetOrigin()[0] << " " << imageGeo->GetOrigin()[1] << " " << imageGeo->GetOrigin()[2] << std::endl; + if (outImage != nullptr) + { + mitk::IOUtil::SaveImage(outImage, outFile); + } + + + return EXIT_SUCCESS; +} diff --git a/Modules/ImageStatistics/CMakeLists.txt b/Modules/ImageStatistics/CMakeLists.txt index 5264c1c320..fa0a7cf973 100644 --- a/Modules/ImageStatistics/CMakeLists.txt +++ b/Modules/ImageStatistics/CMakeLists.txt @@ -1,13 +1,13 @@ MITK_CREATE_MODULE( - DEPENDS MitkImageExtraction MitkPlanarFigure + DEPENDS MitkImageExtraction MitkPlanarFigure MitkMultilabel PACKAGE_DEPENDS PUBLIC ITK|ITKIOXML PRIVATE ITK|ITKVTK+ITKConvolution # WARNINGS_AS_ERRORS ) if(BUILD_TESTING) add_subdirectory(Testing) endif(BUILD_TESTING) diff --git a/Modules/ImageStatistics/files.cmake b/Modules/ImageStatistics/files.cmake index 29d3c3547a..be094a5f7a 100644 --- a/Modules/ImageStatistics/files.cmake +++ b/Modules/ImageStatistics/files.cmake @@ -1,14 +1,24 @@ set(CPP_FILES mitkImageStatisticsCalculator.cpp + mitkImageStatisticsCalculator2.cpp mitkPointSetStatisticsCalculator.cpp mitkPointSetDifferenceStatisticsCalculator.cpp mitkIntensityProfile.cpp + mitkHotspotCalculator.cpp + mitkMaskGenerator.cpp + mitkPlanarFigureMaskGenerator.cpp + mitkMultiLabelMaskGenerator.cpp ) set(H_FILES mitkImageStatisticsCalculator.h + mitkImageStatisticsCalculator2.h mitkPointSetDifferenceStatisticsCalculator.h mitkPointSetStatisticsCalculator.h mitkExtendedStatisticsImageFilter.h mitkExtendedLabelStatisticsImageFilter.h + mitkHotspotCalculator.h + mitkMaskGenerator.h + mitkPlanarFigureMaskGenerator.h + mitkMultiLabelMaskGenerator.h ) diff --git a/Modules/ImageStatistics/mitkHotspotCalculator.cpp b/Modules/ImageStatistics/mitkHotspotCalculator.cpp new file mode 100644 index 0000000000..e385b36903 --- /dev/null +++ b/Modules/ImageStatistics/mitkHotspotCalculator.cpp @@ -0,0 +1,524 @@ +#include +#include +#include +#include +#include +#include "mitkImageAccessByItk.h" +#include +#include +#include + +namespace mitk +{ + HotspotCalculator::HotspotCalculator(): + m_HotspotRadiusinMM(6.2035049089940), // radius of a 1cm3 sphere in mm + m_HotspotMustBeCompletelyInsideImage(true), + m_HotspotParamsChanged(true) + { + } + + void HotspotCalculator::setInputImage(mitk::Image::Pointer inputImage) + { + m_Image = inputImage; + } + + void HotspotCalculator::setMaskImage(mitk::Image::Pointer maskImage) + { + m_Mask = maskImage; + } + + HotspotCalculator::~HotspotCalculator() + { + } + + void HotspotCalculator::setHotspotRadiusInMM(double radiusInMillimeter) + { + if(radiusInMillimeter != m_HotspotRadiusinMM) + { + m_HotspotRadiusinMM = radiusInMillimeter; + m_HotspotParamsChanged = true; + } + } + + const double& HotspotCalculator::getHotspotRadiusinMM() const + { + return m_HotspotRadiusinMM; + } + + bool HotspotCalculator::getHotspotMustBeCompletelyInsideImage() const + { + return m_HotspotMustBeCompletelyInsideImage; + } + + mitk::Image::Pointer HotspotCalculator::getHotspotMask(const unsigned int timeStep, + unsigned int label) + { + if ( m_Image.IsNull() ) + { + throw std::runtime_error( "Error: image empty!" ); + } + + if ( timeStep >= m_Image->GetTimeSteps() ) + { + throw std::runtime_error( "Error: invalid time step!" ); + } + + mitk::ImageTimeSelector::Pointer imageTimeSelector = mitk::ImageTimeSelector::New(); + imageTimeSelector->SetInput( m_Image ); + imageTimeSelector->SetTimeNr( timeStep ); + imageTimeSelector->UpdateLargestPossibleRegion(); + mitk::Image::Pointer timeSliceImage = imageTimeSelector->GetOutput(); + + m_internalImage = timeSliceImage; + m_internalMask2D = nullptr; // is this correct when this variable holds a smart pointer? + m_internalMask3D = nullptr; + + if ( m_Mask != nullptr ) + { + unsigned int timeStep_mask = timeStep; + if ( timeStep >= m_Mask->GetTimeSteps() ) + { + timeStep_mask = m_Mask->GetTimeSteps(); + } + + mitk::ImageTimeSelector::Pointer maskTimeSelector = mitk::ImageTimeSelector::New(); + maskTimeSelector->SetInput( m_Mask ); + maskTimeSelector->SetTimeNr( timeStep_mask ); + maskTimeSelector->UpdateLargestPossibleRegion(); + mitk::Image::Pointer timeSliceMask = maskTimeSelector->GetOutput(); + + if ( m_internalImage->GetDimension() == 3 ) + { + CastToItkImage(timeSliceMask, m_internalMask3D); + AccessFixedDimensionByItk_2(m_internalImage, CalculateHotspotMask, 3, m_internalMask3D, label); + } + else if ( m_internalImage->GetDimension() == 2 ) + { + CastToItkImage(timeSliceMask, m_internalMask2D); + AccessFixedDimensionByItk_2(m_internalImage, CalculateHotspotMask, 2, m_internalMask2D, label); + } + else + { + throw std::runtime_error( "Error: invalid image dimension" ); + } + } + else + { + + if ( m_internalImage->GetDimension() == 3 ) + { + AccessFixedDimensionByItk_2(m_internalImage, CalculateHotspotMask, 3, m_internalMask3D, label); + } + else if ( m_internalImage->GetDimension() == 2 ) + { + AccessFixedDimensionByItk_2(m_internalImage, CalculateHotspotMask, 2, m_internalMask2D, label); + } + else + { + throw std::runtime_error( "Error: invalid image dimension" ); + } + } + + return m_hotspotMask; + } + + template + HotspotCalculator::ImageExtrema + HotspotCalculator::CalculateExtremaWorld( const itk::Image* inputImage, + typename itk::Image::Pointer maskImage, + double neccessaryDistanceToImageBorderInMM, + unsigned int label ) + { + typedef itk::Image< TPixel, VImageDimension > ImageType; + typedef itk::Image< unsigned short, VImageDimension > MaskImageType; + + typedef itk::ImageRegionConstIteratorWithIndex MaskImageIteratorType; + typedef itk::ImageRegionConstIteratorWithIndex InputImageIndexIteratorType; + + typename ImageType::SpacingType spacing = inputImage->GetSpacing(); + + ImageExtrema minMax; + minMax.Defined = false; + minMax.MaxIndex.set_size(VImageDimension); + minMax.MaxIndex.set_size(VImageDimension); + + typename ImageType::RegionType allowedExtremaRegion = inputImage->GetLargestPossibleRegion(); + + bool keepDistanceToImageBorders( neccessaryDistanceToImageBorderInMM > 0 ); + if (keepDistanceToImageBorders) + { + long distanceInPixels[VImageDimension]; + for(unsigned short dimension = 0; dimension < VImageDimension; ++dimension) + { + // To confirm that the whole hotspot is inside the image we have to keep a specific distance to the image-borders, which is as long as + // the radius. To get the amount of indices we divide the radius by spacing and add 0.5 because voxels are center based: + // For example with a radius of 2.2 and a spacing of 1 two indices are enough because 2.2 / 1 + 0.5 = 2.7 => 2. + // But with a radius of 2.7 we need 3 indices because 2.7 / 1 + 0.5 = 3.2 => 3 + distanceInPixels[dimension] = int( neccessaryDistanceToImageBorderInMM / spacing[dimension] + 0.5); + } + + allowedExtremaRegion.ShrinkByRadius(distanceInPixels); + } + + InputImageIndexIteratorType imageIndexIt(inputImage, allowedExtremaRegion); + + float maxValue = itk::NumericTraits::min(); + float minValue = itk::NumericTraits::max(); + + typename ImageType::IndexType maxIndex; + typename ImageType::IndexType minIndex; + + for(unsigned short i = 0; i < VImageDimension; ++i) + { + maxIndex[i] = 0; + minIndex[i] = 0; + } + + if (maskImage != nullptr) + { + MaskImageIteratorType maskIt(maskImage, maskImage->GetLargestPossibleRegion()); + typename ImageType::IndexType imageIndex; + typename ImageType::PointType worldPosition; + typename ImageType::IndexType maskIndex; + + for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt) + { + imageIndex = maskIndex = maskIt.GetIndex(); + + if(maskIt.Get() == label) + { + if( allowedExtremaRegion.IsInside(imageIndex) ) + { + imageIndexIt.SetIndex( imageIndex ); + double value = imageIndexIt.Get(); + minMax.Defined = true; + + //Calculate minimum, maximum and corresponding index-values + if( value > maxValue ) + { + maxIndex = imageIndexIt.GetIndex(); + maxValue = value; + } + + if(value < minValue ) + { + minIndex = imageIndexIt.GetIndex(); + minValue = value; + } + } + } + } + } + else + { + for(imageIndexIt.GoToBegin(); !imageIndexIt.IsAtEnd(); ++imageIndexIt) + { + double value = imageIndexIt.Get(); + minMax.Defined = true; + + //Calculate minimum, maximum and corresponding index-values + if( value > maxValue ) + { + maxIndex = imageIndexIt.GetIndex(); + maxValue = value; + } + + if(value < minValue ) + { + minIndex = imageIndexIt.GetIndex(); + minValue = value; + } + } + } + + minMax.MaxIndex.set_size(VImageDimension); + minMax.MinIndex.set_size(VImageDimension); + + for(unsigned int i = 0; i < minMax.MaxIndex.size(); ++i) + { + minMax.MaxIndex[i] = maxIndex[i]; + } + + for(unsigned int i = 0; i < minMax.MinIndex.size(); ++i) + { + minMax.MinIndex[i] = minIndex[i]; + } + + minMax.Max = maxValue; + minMax.Min = minValue; + + return minMax; + } + + template + itk::Size + HotspotCalculator::CalculateConvolutionKernelSize( double spacing[VImageDimension], + double radiusInMM ) + { + typedef itk::Image< float, VImageDimension > KernelImageType; + typedef typename KernelImageType::SizeType SizeType; + SizeType maskSize; + + for(unsigned int i = 0; i < VImageDimension; ++i) + { + maskSize[i] = static_cast( 2 * radiusInMM / spacing[i]); + + // We always want an uneven size to have a clear center point in the convolution mask + if(maskSize[i] % 2 == 0 ) + { + ++maskSize[i]; + } + } + return maskSize; + } + + template + itk::SmartPointer< itk::Image > + HotspotCalculator::GenerateHotspotSearchConvolutionKernel(double mmPerPixel[VImageDimension], + double radiusInMM ) + { + std::stringstream ss; + for (unsigned int i = 0; i < VImageDimension; ++i) + { + ss << mmPerPixel[i]; + if (i < VImageDimension -1) + ss << ","; + } + MITK_DEBUG << "Update convolution kernel for spacing (" << ss.str() << ") and radius " << radiusInMM << "mm"; + + + double radiusInMMSquared = radiusInMM * radiusInMM; + typedef itk::Image< float, VImageDimension > KernelImageType; + typename KernelImageType::Pointer convolutionKernel = KernelImageType::New(); + + // Calculate size and allocate mask image + typedef typename KernelImageType::SizeType SizeType; + SizeType maskSize = this->CalculateConvolutionKernelSize(mmPerPixel, radiusInMM); + + mitk::Point3D convolutionMaskCenterIndex; + convolutionMaskCenterIndex.Fill(0.0); + for(unsigned int i = 0; i < VImageDimension; ++i) + { + convolutionMaskCenterIndex[i] = 0.5 * (double)(maskSize[i]-1); + } + + typedef typename KernelImageType::IndexType IndexType; + IndexType maskIndex; + maskIndex.Fill(0); + + typedef typename KernelImageType::RegionType RegionType; + RegionType maskRegion; + maskRegion.SetSize(maskSize); + maskRegion.SetIndex(maskIndex); + + convolutionKernel->SetRegions(maskRegion); + convolutionKernel->SetSpacing(mmPerPixel); + convolutionKernel->Allocate(); + + // Fill mask image values by subsampling the image grid + typedef itk::ImageRegionIteratorWithIndex MaskIteratorType; + MaskIteratorType maskIt(convolutionKernel,maskRegion); + + int numberOfSubVoxelsPerDimension = 2; // per dimension! + int numberOfSubVoxels = ::pow( static_cast(numberOfSubVoxelsPerDimension), static_cast(VImageDimension) ); + double subVoxelSizeInPixels = 1.0 / (double)numberOfSubVoxelsPerDimension; + double valueOfOneSubVoxel = 1.0 / (double)numberOfSubVoxels; + double maskValue = 0.0; + mitk::Point3D subVoxelIndexPosition; + double distanceSquared = 0.0; + + typedef itk::ContinuousIndex ContinuousIndexType; + for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt) + { + ContinuousIndexType indexPoint(maskIt.GetIndex()); + mitk::Point3D voxelPosition; + for (unsigned int dimension = 0; dimension < VImageDimension; ++dimension) + { + voxelPosition[dimension] = indexPoint[dimension]; + } + + maskValue = 0.0; + mitk::Vector3D subVoxelOffset; subVoxelOffset.Fill(0.0); + // iterate sub-voxels by iterating all possible offsets + for (subVoxelOffset[0] = -0.5 + subVoxelSizeInPixels / 2.0; + subVoxelOffset[0] < +0.5; + subVoxelOffset[0] += subVoxelSizeInPixels) + { + for (subVoxelOffset[1] = -0.5 + subVoxelSizeInPixels / 2.0; + subVoxelOffset[1] < +0.5; + subVoxelOffset[1] += subVoxelSizeInPixels) + { + for (subVoxelOffset[2] = -0.5 + subVoxelSizeInPixels / 2.0; + subVoxelOffset[2] < +0.5; + subVoxelOffset[2] += subVoxelSizeInPixels) + { + subVoxelIndexPosition = voxelPosition + subVoxelOffset; // this COULD be integrated into the for-loops if neccessary (add voxelPosition to initializer and end condition) + distanceSquared = + (subVoxelIndexPosition[0]-convolutionMaskCenterIndex[0]) * mmPerPixel[0] * (subVoxelIndexPosition[0]-convolutionMaskCenterIndex[0]) * mmPerPixel[0] + + (subVoxelIndexPosition[1]-convolutionMaskCenterIndex[1]) * mmPerPixel[1] * (subVoxelIndexPosition[1]-convolutionMaskCenterIndex[1]) * mmPerPixel[1] + + (subVoxelIndexPosition[2]-convolutionMaskCenterIndex[2]) * mmPerPixel[2] * (subVoxelIndexPosition[2]-convolutionMaskCenterIndex[2]) * mmPerPixel[2]; + + if (distanceSquared <= radiusInMMSquared) + { + maskValue += valueOfOneSubVoxel; + } + } + } + } + maskIt.Set( maskValue ); + } + + return convolutionKernel; + } + + template + itk::SmartPointer > + HotspotCalculator::GenerateConvolutionImage( const itk::Image* inputImage ) + { + double mmPerPixel[VImageDimension]; + for (unsigned int dimension = 0; dimension < VImageDimension; ++dimension) + { + mmPerPixel[dimension] = inputImage->GetSpacing()[dimension]; + } + + // update convolution kernel + typedef itk::Image< float, VImageDimension > KernelImageType; + typename KernelImageType::Pointer convolutionKernel = this->GenerateHotspotSearchConvolutionKernel(mmPerPixel, m_HotspotRadiusinMM); + + // update convolution image + typedef itk::Image< TPixel, VImageDimension > InputImageType; + typedef itk::Image< TPixel, VImageDimension > ConvolutionImageType; + typedef itk::FFTConvolutionImageFilter ConvolutionFilterType; + + typename ConvolutionFilterType::Pointer convolutionFilter = ConvolutionFilterType::New(); + typedef itk::ConstantBoundaryCondition BoundaryConditionType; + BoundaryConditionType boundaryCondition; + boundaryCondition.SetConstant(0.0); + + if (m_HotspotMustBeCompletelyInsideImage) + { + // overwrite default boundary condition + convolutionFilter->SetBoundaryCondition(&boundaryCondition); + } + + convolutionFilter->SetInput(inputImage); + convolutionFilter->SetKernelImage(convolutionKernel); + convolutionFilter->SetNormalize(true); + MITK_DEBUG << "Update Convolution image for hotspot search"; + convolutionFilter->UpdateLargestPossibleRegion(); + + typename ConvolutionImageType::Pointer convolutionImage = convolutionFilter->GetOutput(); + convolutionImage->SetSpacing( inputImage->GetSpacing() ); // only workaround because convolution filter seems to ignore spacing of input image + + m_modified = false; + return convolutionImage; + } + + template < typename TPixel, unsigned int VImageDimension> + void + HotspotCalculator::FillHotspotMaskPixels( itk::Image* maskImage, + itk::Point sphereCenter, + double sphereRadiusInMM ) + { + typedef itk::Image< TPixel, VImageDimension > MaskImageType; + typedef itk::ImageRegionIteratorWithIndex MaskImageIteratorType; + + MaskImageIteratorType maskIt(maskImage, maskImage->GetLargestPossibleRegion()); + + typename MaskImageType::IndexType maskIndex; + typename MaskImageType::PointType worldPosition; + + // this is not very smart. I would rather use a 0 initialized mask (not the case here -> blame CalculateHotspotMask) and find the region where I need to iterate over, then iterate only over the small region + for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt) + { + maskIndex = maskIt.GetIndex(); + maskImage->TransformIndexToPhysicalPoint(maskIndex, worldPosition); + maskIt.Set( worldPosition.EuclideanDistanceTo(sphereCenter) <= sphereRadiusInMM ? 1 : 0 ); + } + } + + template + void + HotspotCalculator::CalculateHotspotMask(itk::Image* inputImage, + typename itk::Image::Pointer maskImage, + unsigned int label) + { + typedef itk::Image< TPixel, VImageDimension > InputImageType; + typedef itk::Image< TPixel, VImageDimension > ConvolutionImageType; + typedef itk::Image< float, VImageDimension > KernelImageType; + typedef itk::Image< unsigned short, VImageDimension > MaskImageType; + typedef itk::ImageDuplicator< InputImageType > DuplicatorType; + typedef itk::CastImageFilter< InputImageType, MaskImageType > CastFilterType; + + typename ConvolutionImageType::Pointer convolutionImage = this->GenerateConvolutionImage(inputImage); + + if (convolutionImage.IsNull()) + { + MITK_ERROR << "Empty convolution image in CalculateHotspotStatistics(). We should never reach this state (logic error)."; + throw std::logic_error("Empty convolution image in CalculateHotspotStatistics()"); + } + + // if mask image is not defined, create an image of the same size as inputImage and fill it with 1's + // there is maybe a better way to do this!? + if (maskImage == nullptr) + { + maskImage = MaskImageType::New(); + typename MaskImageType::RegionType maskRegion = inputImage->GetLargestPossibleRegion(); + typename MaskImageType::SpacingType maskSpacing = inputImage->GetSpacing(); + typename MaskImageType::PointType maskOrigin = inputImage->GetOrigin(); + typename MaskImageType::DirectionType maskDirection = inputImage->GetDirection(); + maskImage->SetRegions(maskRegion); + maskImage->Allocate(); + maskImage->SetOrigin(maskOrigin); + maskImage->SetSpacing(maskSpacing); + maskImage->SetDirection(maskDirection); + + maskImage->FillBuffer(1); + + label = 1; + } + + // find maximum in convolution image, given the current mask + double requiredDistanceToBorder = m_HotspotMustBeCompletelyInsideImage ? m_HotspotRadiusinMM : -1.0; + ImageExtrema convolutionImageInformation = CalculateExtremaWorld(convolutionImage.GetPointer(), maskImage, requiredDistanceToBorder, label); + + bool isHotspotDefined = convolutionImageInformation.Defined; + + if (!isHotspotDefined) + { + MITK_ERROR << "No origin of hotspot-sphere was calculated!"; + m_hotspotMask = nullptr; + } + else + { + // create a binary mask around the "hotspot" region, fill the shape of a sphere around our hotspot center + typename DuplicatorType::Pointer copyMachine = DuplicatorType::New(); + copyMachine->SetInputImage(inputImage); + copyMachine->Update(); + + typename CastFilterType::Pointer caster = CastFilterType::New(); + caster->SetInput( copyMachine->GetOutput() ); + caster->Update(); + typename MaskImageType::Pointer hotspotMaskITK = caster->GetOutput(); + + typedef typename InputImageType::IndexType IndexType; + IndexType maskCenterIndex; + for (unsigned int d =0; d< VImageDimension;++d) + { + maskCenterIndex[d]=convolutionImageInformation.MaxIndex[d]; + } + + typename ConvolutionImageType::PointType maskCenter; + inputImage->TransformIndexToPhysicalPoint(maskCenterIndex,maskCenter); + + FillHotspotMaskPixels(hotspotMaskITK.GetPointer(), maskCenter, m_HotspotRadiusinMM); + + //obtain mitk::Image::Pointer from itk::Image + mitk::Image::Pointer hotspotMaskAsMITKImage = mitk::GrabItkImageMemory(hotspotMaskITK); + + m_hotspotMask = hotspotMaskAsMITKImage; + } + } +} diff --git a/Modules/ImageStatistics/mitkHotspotCalculator.h b/Modules/ImageStatistics/mitkHotspotCalculator.h new file mode 100644 index 0000000000..1dfd7514b6 --- /dev/null +++ b/Modules/ImageStatistics/mitkHotspotCalculator.h @@ -0,0 +1,110 @@ +#ifndef MITKHOTSPOTCALCULATOR_H +#define MITKHOTSPOTCALCULATOR_H + +#include +#include +#include +#include +#include +#include +#include + + +namespace mitk +{ + class MITKIMAGESTATISTICS_EXPORT HotspotCalculator: itk::Object + { + public: + HotspotCalculator(); + + ~HotspotCalculator(); + + void setInputImage(mitk::Image::Pointer inputImage); + + void setMaskImage(mitk::Image::Pointer maskImage); + + void setHotspotRadiusInMM(double radiusInMillimeter); + + const double& getHotspotRadiusinMM() const; + + void setHotspotMustBeCompletelyInsideImage(bool hotspotCompletelyInsideImage); + + bool getHotspotMustBeCompletelyInsideImage() const; + + mitk::Image::Pointer getHotspotMask(unsigned int timeStep, unsigned int label=0); + + protected: + class ImageExtrema + { + public: + bool Defined; + double Max; + double Min; + vnl_vector MaxIndex; + vnl_vector MinIndex; + + ImageExtrema() + :Defined(false) + ,Max(itk::NumericTraits::min()) + ,Min(itk::NumericTraits::max()) + { + } + }; + + private: + /** \brief Returns size of convolution kernel depending on spacing and radius. */ + template + itk::Size + CalculateConvolutionKernelSize(double spacing[VImageDimension], double radiusInMM); + + /** \brief Generates image of kernel which is needed for convolution. */ + template + itk::SmartPointer< itk::Image > + GenerateHotspotSearchConvolutionKernel(double spacing[VImageDimension], double radiusInMM); + + /** \brief Convolves image with spherical kernel image. Used for hotspot calculation. */ + template + itk::SmartPointer< itk::Image > + GenerateConvolutionImage( const itk::Image* inputImage ); + + + /** \brief Fills pixels of the spherical hotspot mask. */ + template < typename TPixel, unsigned int VImageDimension> + void + FillHotspotMaskPixels( itk::Image* maskImage, + itk::Point sphereCenter, + double sphereRadiusInMM); + + + /** \brief */ + template + void + CalculateHotspotMask(itk::Image* inputImage, + typename itk::Image::Pointer maskImage, + unsigned int label); + + + template + ImageExtrema CalculateExtremaWorld( const itk::Image* inputImage, + typename itk::Image::Pointer maskImage, + double neccessaryDistanceToImageBorderInMM, + unsigned int label); + + HotspotCalculator(const HotspotCalculator &); + HotspotCalculator & operator=(const HotspotCalculator &); + + mitk::Image::Pointer m_Image; + mitk::Image::Pointer m_Mask = nullptr; + mitk::Image::Pointer m_internalImage; + itk::Image::Pointer m_internalMask2D; + itk::Image::Pointer m_internalMask3D; + mitk::Image::Pointer m_hotspotMask; + double m_HotspotRadiusinMM; + bool m_HotspotMustBeCompletelyInsideImage; + bool m_HotspotParamsChanged; + bool m_modified = false; + }; +} +#endif // MITKHOTSPOTCALCULATOR + + diff --git a/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp b/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp index 03bd6e601e..de7f4e94b1 100644 --- a/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp +++ b/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp @@ -1,2226 +1,2232 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkImageStatisticsCalculator.h" #include "mitkImageAccessByItk.h" #include "mitkImageCast.h" #include "mitkExtractImageFilter.h" #include "mitkImageTimeSelector.h" #include "mitkITKImageImport.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include "itkImage.h" //#define DEBUG_HOTSPOTSEARCH #define _USE_MATH_DEFINES #include #include "vtkLassoStencilSource.h" namespace mitk { ImageStatisticsCalculator::ImageStatisticsCalculator() : m_MaskingMode( MASKING_MODE_NONE ), m_MaskingModeChanged( false ), m_IgnorePixelValue(0.0), m_DoIgnorePixelValue(false), m_IgnorePixelValueChanged(false), m_PlanarFigureAxis (0), m_PlanarFigureSlice (0), m_PlanarFigureCoordinate0 (0), m_PlanarFigureCoordinate1 (0), m_HistogramBinSize(1.0), m_UseDefaultBinSize(true), m_UseBinSizeBasedOnVOIRegion(false), m_HotspotRadiusInMM(6.2035049089940), // radius of a 1cm3 sphere in mm m_CalculateHotspot(false), m_HotspotRadiusInMMChanged(false), m_HotspotMustBeCompletelyInsideImage(true) { m_EmptyHistogram = HistogramType::New(); m_EmptyHistogram->SetMeasurementVectorSize(1); HistogramType::SizeType histogramSize(1); histogramSize.Fill( 256 ); m_EmptyHistogram->Initialize( histogramSize ); m_EmptyStatistics.Reset(); } ImageStatisticsCalculator::~ImageStatisticsCalculator() { } void ImageStatisticsCalculator::SetUseDefaultBinSize(bool useDefault) { m_UseDefaultBinSize = useDefault; } ImageStatisticsCalculator::Statistics::Statistics(bool withHotspotStatistics) :m_HotspotStatistics(withHotspotStatistics ? new Statistics(false) : nullptr) { Reset(); } ImageStatisticsCalculator::Statistics::Statistics(const Statistics& other) :m_HotspotStatistics( nullptr) { this->SetLabel( other.GetLabel() ); this->SetN( other.GetN() ); this->SetMin( other.GetMin() ); this->SetMax( other.GetMax() ); this->SetMedian( other.GetMedian() ); this->SetMean( other.GetMean() ); this->SetVariance( other.GetVariance() ); this->SetKurtosis( other.GetKurtosis() ); this->SetSkewness( other.GetSkewness() ); this->SetUniformity( other.GetUniformity() ); this->SetEntropy( other.GetEntropy() ); this->SetUPP( other.GetUPP() ); this->SetMPP( other.GetMPP() ); this->SetSigma( other.GetSigma() ); this->SetRMS( other.GetRMS() ); this->SetMaxIndex( other.GetMaxIndex() ); this->SetMinIndex( other.GetMinIndex() ); this->SetHotspotIndex( other.GetHotspotIndex() ); if (other.m_HotspotStatistics) { this->m_HotspotStatistics = new Statistics(false); *this->m_HotspotStatistics = *other.m_HotspotStatistics; } } bool ImageStatisticsCalculator::Statistics::HasHotspotStatistics() const { return m_HotspotStatistics != nullptr; } void ImageStatisticsCalculator::Statistics::SetHasHotspotStatistics(bool hasHotspotStatistics) { m_HasHotspotStatistics = hasHotspotStatistics; } ImageStatisticsCalculator::Statistics::~Statistics() { delete m_HotspotStatistics; } double ImageStatisticsCalculator::Statistics::GetVariance() const { return this->Variance; } void ImageStatisticsCalculator::Statistics::SetVariance( const double value ) { if( this->Variance != value ) { if( value < 0.0 ) { this->Variance = 0.0; // if given value is negative set variance to 0.0 } else { this->Variance = value; } } } double ImageStatisticsCalculator::Statistics::GetSigma() const { return this->Sigma; } void ImageStatisticsCalculator::Statistics::SetSigma( const double value ) { if( this->Sigma != value ) { // for some compiler the value != value works to check for NaN but not for all // but we can always be sure that the standard deviation is a positive value if( value != value || value < 0.0 ) { // if standard deviation is NaN we just assume 0.0 this->Sigma = 0.0; } else { this->Sigma = value; } } } void ImageStatisticsCalculator::Statistics::Reset(unsigned int dimension) { SetLabel(0); SetN( 0 ); SetMin( 0.0 ); SetMax( 0.0 ); SetMedian( 0.0 ); SetVariance( 0.0 ); SetMean( 0.0 ); SetSigma( 0.0 ); SetRMS( 0.0 ); SetSkewness( 0.0 ); SetKurtosis( 0.0 ); SetUniformity( 0.0 ); SetEntropy( 0.0 ); SetMPP( 0.0 ); SetUPP( 0.0 ); vnl_vector zero; zero.set_size(dimension); for(unsigned int i = 0; i < dimension; ++i) { zero[i] = 0; } SetMaxIndex(zero); SetMinIndex(zero); SetHotspotIndex(zero); if (m_HotspotStatistics != nullptr) { m_HotspotStatistics->Reset(dimension); } } const ImageStatisticsCalculator::Statistics& ImageStatisticsCalculator::Statistics::GetHotspotStatistics() const { if (m_HotspotStatistics) { return *m_HotspotStatistics; } else { throw std::logic_error("Object has no hostspot statistics, see HasHotspotStatistics()"); } } ImageStatisticsCalculator::Statistics& ImageStatisticsCalculator::Statistics::GetHotspotStatistics() { if (m_HotspotStatistics) { return *m_HotspotStatistics; } else { throw std::logic_error("Object has no hostspot statistics, see HasHotspotStatistics()"); } } ImageStatisticsCalculator::Statistics& ImageStatisticsCalculator::Statistics::operator=(ImageStatisticsCalculator::Statistics const& other) { if (this == &other) return *this; this->SetLabel( other.GetLabel() ); this->SetN( other.GetN() ); this->SetMin( other.GetMin() ); this->SetMax( other.GetMax() ); this->SetMean( other.GetMean() ); this->SetMedian( other.GetMedian() ); this->SetVariance( other.GetVariance() ); this->SetSigma( other.GetSigma() ); this->SetRMS( other.GetRMS() ); this->SetMinIndex( other.GetMinIndex() ); this->SetMaxIndex( other.GetMaxIndex() ); this->SetHotspotIndex( other.GetHotspotIndex() ); this->SetSkewness( other.GetSkewness() ); this->SetKurtosis( other.GetKurtosis() ); this->SetUniformity( other.GetUniformity() ); this->SetEntropy( other.GetEntropy() ); this->SetUPP( other.GetUPP() ); this->SetMPP( other.GetMPP() ); delete this->m_HotspotStatistics; this->m_HotspotStatistics = nullptr; if (other.m_HotspotStatistics) { this->m_HotspotStatistics = new Statistics(false); *this->m_HotspotStatistics = *other.m_HotspotStatistics; } return *this; } void ImageStatisticsCalculator::SetImage( const mitk::Image *image ) { if ( m_Image != image ) { m_Image = image; this->Modified(); unsigned int numberOfTimeSteps = image->GetTimeSteps(); // Initialize vectors to time-size of this image m_ImageHistogramVector.resize( numberOfTimeSteps ); m_MaskedImageHistogramVector.resize( numberOfTimeSteps ); m_PlanarFigureHistogramVector.resize( numberOfTimeSteps ); m_ImageStatisticsVector.resize( numberOfTimeSteps ); m_MaskedImageStatisticsVector.resize( numberOfTimeSteps ); m_PlanarFigureStatisticsVector.resize( numberOfTimeSteps ); m_ImageStatisticsTimeStampVector.resize( numberOfTimeSteps ); m_MaskedImageStatisticsTimeStampVector.resize( numberOfTimeSteps ); m_PlanarFigureStatisticsTimeStampVector.resize( numberOfTimeSteps ); m_ImageStatisticsCalculationTriggerVector.resize( numberOfTimeSteps ); m_MaskedImageStatisticsCalculationTriggerVector.resize( numberOfTimeSteps ); m_PlanarFigureStatisticsCalculationTriggerVector.resize( numberOfTimeSteps ); for ( unsigned int t = 0; t < image->GetTimeSteps(); ++t ) { m_ImageStatisticsTimeStampVector[t].Modified(); m_ImageStatisticsCalculationTriggerVector[t] = true; } } } void ImageStatisticsCalculator::SetImageMask( const mitk::Image *imageMask ) { if ( m_Image.IsNull() ) { itkExceptionMacro( << "Image needs to be set first!" ); } if ( m_ImageMask != imageMask ) { m_ImageMask = imageMask; this->Modified(); for ( unsigned int t = 0; t < m_Image->GetTimeSteps(); ++t ) { m_MaskedImageStatisticsTimeStampVector[t].Modified(); m_MaskedImageStatisticsCalculationTriggerVector[t] = true; } } } void ImageStatisticsCalculator::SetPlanarFigure( mitk::PlanarFigure *planarFigure ) { if ( m_Image.IsNull() ) { itkExceptionMacro( << "Image needs to be set first!" ); } if ( m_PlanarFigure != planarFigure ) { m_PlanarFigure = planarFigure; this->Modified(); for ( unsigned int t = 0; t < m_Image->GetTimeSteps(); ++t ) { m_PlanarFigureStatisticsTimeStampVector[t].Modified(); m_PlanarFigureStatisticsCalculationTriggerVector[t] = true; } } } void ImageStatisticsCalculator::SetMaskingMode( unsigned int mode ) { if ( m_MaskingMode != mode ) { m_MaskingMode = mode; m_MaskingModeChanged = true; this->Modified(); } } void ImageStatisticsCalculator::SetMaskingModeToNone() { if ( m_MaskingMode != MASKING_MODE_NONE ) { m_MaskingMode = MASKING_MODE_NONE; m_MaskingModeChanged = true; this->Modified(); } } void ImageStatisticsCalculator::SetMaskingModeToImage() { if ( m_MaskingMode != MASKING_MODE_IMAGE ) { m_MaskingMode = MASKING_MODE_IMAGE; m_MaskingModeChanged = true; this->Modified(); } } void ImageStatisticsCalculator::SetMaskingModeToPlanarFigure() { if ( m_MaskingMode != MASKING_MODE_PLANARFIGURE ) { m_MaskingMode = MASKING_MODE_PLANARFIGURE; m_MaskingModeChanged = true; this->Modified(); } } void ImageStatisticsCalculator::SetIgnorePixelValue(double value) { if ( m_IgnorePixelValue != value ) { m_IgnorePixelValue = value; if(m_DoIgnorePixelValue) { m_IgnorePixelValueChanged = true; } this->Modified(); } } double ImageStatisticsCalculator::GetIgnorePixelValue() { return m_IgnorePixelValue; } void ImageStatisticsCalculator::SetDoIgnorePixelValue(bool value) { if ( m_DoIgnorePixelValue != value ) { m_DoIgnorePixelValue = value; m_IgnorePixelValueChanged = true; this->Modified(); } } bool ImageStatisticsCalculator::GetDoIgnorePixelValue() { return m_DoIgnorePixelValue; } void ImageStatisticsCalculator::SetHistogramBinSize(double size) { this->m_HistogramBinSize = size; } double ImageStatisticsCalculator::GetHistogramBinSize() { return this->m_HistogramBinSize; } void ImageStatisticsCalculator::SetHotspotRadiusInMM(double value) { if ( m_HotspotRadiusInMM != value ) { m_HotspotRadiusInMM = value; if(m_CalculateHotspot) { m_HotspotRadiusInMMChanged = true; //MITK_INFO <<"Hotspot radius changed, new convolution required"; } this->Modified(); } } double ImageStatisticsCalculator::GetHotspotRadiusInMM() { return m_HotspotRadiusInMM; } void ImageStatisticsCalculator::SetCalculateHotspot(bool on) { if ( m_CalculateHotspot != on ) { m_CalculateHotspot = on; m_HotspotRadiusInMMChanged = true; //MITK_INFO <<"Hotspot calculation changed, new convolution required"; this->Modified(); } } bool ImageStatisticsCalculator::IsHotspotCalculated() { return m_CalculateHotspot; } void ImageStatisticsCalculator::SetHotspotMustBeCompletlyInsideImage(bool hotspotMustBeCompletelyInsideImage, bool warn) { m_HotspotMustBeCompletelyInsideImage = hotspotMustBeCompletelyInsideImage; if (!m_HotspotMustBeCompletelyInsideImage && warn) { MITK_WARN << "Hotspot calculation will extrapolate pixels at image borders. Be aware of the consequences for the hotspot location."; } } bool ImageStatisticsCalculator::GetHotspotMustBeCompletlyInsideImage() const { return m_HotspotMustBeCompletelyInsideImage; } /* Implementation of the min max values for setting the range of the histogram */ template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator::GetMinAndMaxValue( double &min, double &max, int &counter, double &sigma, const itk::Image< TPixel, VImageDimension > *InputImage, itk::Image< unsigned short, VImageDimension > *MaskedImage ) { typedef itk::Image< unsigned short, VImageDimension > MaskImageType; typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::ImageRegionConstIteratorWithIndex Imageie; typedef itk::ImageRegionConstIteratorWithIndex Imageie2; Imageie2 labelIterator2( MaskedImage, MaskedImage->GetRequestedRegion() ); Imageie labelIterator3( InputImage, InputImage->GetRequestedRegion() ); max = 0; min = 0; counter = 0; sigma = 0; double SumOfSquares = 0; double sumSquared = 0; double actualPielValue = 0; int counterOfPixelsInROI = 0; for( labelIterator2.GoToBegin(); !labelIterator2.IsAtEnd(); ++labelIterator2, ++labelIterator3) { if( labelIterator2.Value()== 1.0) { counter++; counterOfPixelsInROI++; actualPielValue = labelIterator3.Value(); sumSquared = sumSquared + actualPielValue; SumOfSquares = SumOfSquares + std::pow(actualPielValue,2); if(counterOfPixelsInROI == 1) { max = actualPielValue; min = actualPielValue; } if(actualPielValue >= max) { max = actualPielValue; } else if(actualPielValue <= min) { min = actualPielValue; } } } if (counter > 1) { sigma = ( SumOfSquares - std::pow( sumSquared, 2) / counter ) / ( counter-1 ); } else { sigma = 0; } } bool ImageStatisticsCalculator::ComputeStatistics( unsigned int timeStep ) { if (m_Image.IsNull() ) { mitkThrow() << "Image not set!"; } if (!m_Image->IsInitialized()) { mitkThrow() << "Image not initialized!"; } if ( m_Image->GetReferenceCount() == 1 ) { // Image no longer valid; we are the only ones to still hold a reference on it return false; } if ( timeStep >= m_Image->GetTimeSteps() ) { throw std::runtime_error( "Error: invalid time step!" ); } // If a mask was set but we are the only ones to still hold a reference on // it, delete it. if ( m_ImageMask.IsNotNull() && (m_ImageMask->GetReferenceCount() == 1) ) { m_ImageMask = nullptr; } // Check if statistics is already up-to-date unsigned long imageMTime = m_ImageStatisticsTimeStampVector[timeStep].GetMTime(); unsigned long maskedImageMTime = m_MaskedImageStatisticsTimeStampVector[timeStep].GetMTime(); unsigned long planarFigureMTime = m_PlanarFigureStatisticsTimeStampVector[timeStep].GetMTime(); bool imageStatisticsCalculationTrigger = m_ImageStatisticsCalculationTriggerVector[timeStep]; bool maskedImageStatisticsCalculationTrigger = m_MaskedImageStatisticsCalculationTriggerVector[timeStep]; bool planarFigureStatisticsCalculationTrigger = m_PlanarFigureStatisticsCalculationTriggerVector[timeStep]; if ( !m_IgnorePixelValueChanged && !m_HotspotRadiusInMMChanged && ((m_MaskingMode != MASKING_MODE_NONE) || (imageMTime > m_Image->GetMTime() && !imageStatisticsCalculationTrigger)) && ((m_MaskingMode != MASKING_MODE_IMAGE) || (maskedImageMTime > m_ImageMask->GetMTime() && !maskedImageStatisticsCalculationTrigger)) && ((m_MaskingMode != MASKING_MODE_PLANARFIGURE) || (planarFigureMTime > m_PlanarFigure->GetMTime() && !planarFigureStatisticsCalculationTrigger)) ) { // Statistics is up to date! if ( m_MaskingModeChanged ) { m_MaskingModeChanged = false; } else { return false; } } // Reset state changed flag m_MaskingModeChanged = false; m_IgnorePixelValueChanged = false; // Depending on masking mode, extract and/or generate the required image // and mask data from the user input this->ExtractImageAndMask( timeStep ); StatisticsContainer *statisticsContainer; HistogramContainer *histogramContainer; switch ( m_MaskingMode ) { case MASKING_MODE_NONE: default: if(!m_DoIgnorePixelValue) { statisticsContainer = &m_ImageStatisticsVector[timeStep]; histogramContainer = &m_ImageHistogramVector[timeStep]; m_ImageStatisticsTimeStampVector[timeStep].Modified(); m_ImageStatisticsCalculationTriggerVector[timeStep] = false; } else { statisticsContainer = &m_MaskedImageStatisticsVector[timeStep]; histogramContainer = &m_MaskedImageHistogramVector[timeStep]; m_MaskedImageStatisticsTimeStampVector[timeStep].Modified(); m_MaskedImageStatisticsCalculationTriggerVector[timeStep] = false; } break; case MASKING_MODE_IMAGE: statisticsContainer = &m_MaskedImageStatisticsVector[timeStep]; histogramContainer = &m_MaskedImageHistogramVector[timeStep]; m_MaskedImageStatisticsTimeStampVector[timeStep].Modified(); m_MaskedImageStatisticsCalculationTriggerVector[timeStep] = false; break; case MASKING_MODE_PLANARFIGURE: statisticsContainer = &m_PlanarFigureStatisticsVector[timeStep]; histogramContainer = &m_PlanarFigureHistogramVector[timeStep]; m_PlanarFigureStatisticsTimeStampVector[timeStep].Modified(); m_PlanarFigureStatisticsCalculationTriggerVector[timeStep] = false; break; } // Calculate statistics and histogram(s) if ( m_InternalImage->GetDimension() == 3 ) { if ( m_MaskingMode == MASKING_MODE_NONE && !m_DoIgnorePixelValue ) { AccessFixedDimensionByItk_2( m_InternalImage, InternalCalculateStatisticsUnmasked, 3, statisticsContainer, histogramContainer ); } else { AccessFixedDimensionByItk_3( m_InternalImage, InternalCalculateStatisticsMasked, 3, m_InternalImageMask3D.GetPointer(), statisticsContainer, histogramContainer ); } } else if ( m_InternalImage->GetDimension() == 2 ) { if ( m_MaskingMode == MASKING_MODE_NONE && !m_DoIgnorePixelValue ) { AccessFixedDimensionByItk_2( m_InternalImage, InternalCalculateStatisticsUnmasked, 2, statisticsContainer, histogramContainer ); } else { AccessFixedDimensionByItk_3( m_InternalImage, InternalCalculateStatisticsMasked, 2, m_InternalImageMask2D.GetPointer(), statisticsContainer, histogramContainer ); } } else { MITK_ERROR << "ImageStatistics: Image dimension not supported!"; } // Release unused image smart pointers to free memory m_InternalImage = mitk::Image::ConstPointer(); m_InternalImageMask3D = MaskImage3DType::Pointer(); m_InternalImageMask2D = MaskImage2DType::Pointer(); return true; } ImageStatisticsCalculator::BinFrequencyType ImageStatisticsCalculator::GetBinsAndFreuqencyForHistograms( unsigned int timeStep , unsigned int label ) const { const HistogramType *binsAndFrequencyToCalculate = this->GetHistogram(0); // ToDo: map should be created on stack not on heap std::map returnedHistogramMap; unsigned int size = binsAndFrequencyToCalculate->Size(); for( unsigned int bin=0; bin < size; ++bin ) { double frequency = binsAndFrequencyToCalculate->GetFrequency( bin, 0 ); //if( frequency > mitk::eps ) { returnedHistogramMap.insert( std::pair(binsAndFrequencyToCalculate->GetMeasurement( bin, 0 ), binsAndFrequencyToCalculate->GetFrequency( bin, 0 ) ) ); } } return returnedHistogramMap; } const ImageStatisticsCalculator::HistogramType * ImageStatisticsCalculator::GetHistogram( unsigned int timeStep, unsigned int label ) const { if ( m_Image.IsNull() || (timeStep >= m_Image->GetTimeSteps()) ) { return nullptr; } switch ( m_MaskingMode ) { case MASKING_MODE_NONE: default: { if(m_DoIgnorePixelValue) return m_MaskedImageHistogramVector[timeStep][label]; return m_ImageHistogramVector[timeStep][label]; } case MASKING_MODE_IMAGE: return m_MaskedImageHistogramVector[timeStep][label]; case MASKING_MODE_PLANARFIGURE: return m_PlanarFigureHistogramVector[timeStep][label]; } } const ImageStatisticsCalculator::HistogramContainer & ImageStatisticsCalculator::GetHistogramVector( unsigned int timeStep ) const { if ( m_Image.IsNull() || (timeStep >= m_Image->GetTimeSteps()) ) { return m_EmptyHistogramContainer; } switch ( m_MaskingMode ) { case MASKING_MODE_NONE: default: { if(m_DoIgnorePixelValue) return m_MaskedImageHistogramVector[timeStep]; return m_ImageHistogramVector[timeStep]; } case MASKING_MODE_IMAGE: return m_MaskedImageHistogramVector[timeStep]; case MASKING_MODE_PLANARFIGURE: return m_PlanarFigureHistogramVector[timeStep]; } } const ImageStatisticsCalculator::Statistics & ImageStatisticsCalculator::GetStatistics( unsigned int timeStep, unsigned int label ) const { if ( m_Image.IsNull() || (timeStep >= m_Image->GetTimeSteps()) ) { return m_EmptyStatistics; } switch ( m_MaskingMode ) { case MASKING_MODE_NONE: default: { if(m_DoIgnorePixelValue) return m_MaskedImageStatisticsVector[timeStep][label]; return m_ImageStatisticsVector[timeStep][label]; } case MASKING_MODE_IMAGE: return m_MaskedImageStatisticsVector[timeStep][label]; case MASKING_MODE_PLANARFIGURE: return m_PlanarFigureStatisticsVector[timeStep][label]; } } const ImageStatisticsCalculator::StatisticsContainer & ImageStatisticsCalculator::GetStatisticsVector( unsigned int timeStep ) const { if ( m_Image.IsNull() || (timeStep >= m_Image->GetTimeSteps()) ) { return m_EmptyStatisticsContainer; } switch ( m_MaskingMode ) { case MASKING_MODE_NONE: default: { if(m_DoIgnorePixelValue) return m_MaskedImageStatisticsVector[timeStep]; return m_ImageStatisticsVector[timeStep]; } case MASKING_MODE_IMAGE: return m_MaskedImageStatisticsVector[timeStep]; case MASKING_MODE_PLANARFIGURE: return m_PlanarFigureStatisticsVector[timeStep]; } } void ImageStatisticsCalculator::ExtractImageAndMask( unsigned int timeStep ) { if ( m_Image.IsNull() ) { throw std::runtime_error( "Error: image empty!" ); } if ( timeStep >= m_Image->GetTimeSteps() ) { throw std::runtime_error( "Error: invalid time step!" ); } ImageTimeSelector::Pointer imageTimeSelector = ImageTimeSelector::New(); imageTimeSelector->SetInput( m_Image ); imageTimeSelector->SetTimeNr( timeStep ); imageTimeSelector->UpdateLargestPossibleRegion(); mitk::Image *timeSliceImage = imageTimeSelector->GetOutput(); switch ( m_MaskingMode ) { case MASKING_MODE_NONE: { m_InternalImage = timeSliceImage; m_InternalImageMask2D = nullptr; m_InternalImageMask3D = nullptr; if(m_DoIgnorePixelValue) { if( m_InternalImage->GetDimension() == 3 ) { if(itk::ImageIOBase::USHORT != timeSliceImage->GetPixelType().GetComponentType()) CastToItkImage( timeSliceImage, m_InternalImageMask3D ); else CastToItkImage( timeSliceImage->Clone(), m_InternalImageMask3D ); m_InternalImageMask3D->FillBuffer(1); } if( m_InternalImage->GetDimension() == 2 ) { if(itk::ImageIOBase::USHORT != timeSliceImage->GetPixelType().GetComponentType()) CastToItkImage( timeSliceImage, m_InternalImageMask2D ); else CastToItkImage( timeSliceImage->Clone(), m_InternalImageMask2D ); m_InternalImageMask2D->FillBuffer(1); } } break; } case MASKING_MODE_IMAGE: { if ( m_ImageMask.IsNotNull() && (m_ImageMask->GetReferenceCount() > 1) ) { if ( timeStep >= m_ImageMask->GetTimeSteps() ) { // Use the last mask time step in case the current time step is bigger than the total // number of mask time steps. // It makes more sense setting this to the last mask time step than to 0. // For instance if you have a mask with 2 time steps and an image with 5: // If time step 0 is selected, the mask will use time step 0. // If time step 1 is selected, the mask will use time step 1. // If time step 2+ is selected, the mask will use time step 1. // If you have a mask with only one time step instead, this will always default to 0. timeStep = m_ImageMask->GetTimeSteps() - 1; } ImageTimeSelector::Pointer maskedImageTimeSelector = ImageTimeSelector::New(); maskedImageTimeSelector->SetInput( m_ImageMask ); maskedImageTimeSelector->SetTimeNr( timeStep ); maskedImageTimeSelector->UpdateLargestPossibleRegion(); mitk::Image *timeSliceMaskedImage = maskedImageTimeSelector->GetOutput(); m_InternalImage = timeSliceImage; CastToItkImage( timeSliceMaskedImage, m_InternalImageMask3D ); } else { throw std::runtime_error( "Error: image mask empty!" ); } break; } case MASKING_MODE_PLANARFIGURE: { m_InternalImageMask2D = nullptr; if ( m_PlanarFigure.IsNull() ) { throw std::runtime_error( "Error: planar figure empty!" ); } if ( !m_PlanarFigure->IsClosed() ) { throw std::runtime_error( "Masking not possible for non-closed figures" ); } const BaseGeometry *imageGeometry = timeSliceImage->GetGeometry(); if ( imageGeometry == nullptr ) { throw std::runtime_error( "Image geometry invalid!" ); } const PlaneGeometry *planarFigurePlaneGeometry = m_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!" ); } // 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 MaskImage3DType::IndexType index; imageGeometry->WorldToIndex( planarFigureGeometry->GetOrigin(), index ); unsigned int slice = index[axis]; m_PlanarFigureSlice = slice; // Extract slice with given position and direction from image unsigned int dimension = timeSliceImage->GetDimension(); if (dimension != 2) { ExtractImageFilter::Pointer imageExtractor = ExtractImageFilter::New(); imageExtractor->SetInput( timeSliceImage ); imageExtractor->SetSliceDimension( axis ); imageExtractor->SetSliceIndex( slice ); imageExtractor->Update(); m_InternalImage = imageExtractor->GetOutput(); } else { m_InternalImage = timeSliceImage; } // Compute mask from PlanarFigure AccessFixedDimensionByItk_1( m_InternalImage, InternalCalculateMaskFromPlanarFigure, 2, axis ); } } if(m_DoIgnorePixelValue) { if ( m_InternalImage->GetDimension() == 3 ) { AccessFixedDimensionByItk_1( m_InternalImage, InternalMaskIgnoredPixels, 3, m_InternalImageMask3D.GetPointer() ); } else if ( m_InternalImage->GetDimension() == 2 ) { AccessFixedDimensionByItk_1( m_InternalImage, InternalMaskIgnoredPixels, 2, m_InternalImageMask2D.GetPointer() ); } } } bool ImageStatisticsCalculator::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; } unsigned int ImageStatisticsCalculator::calcNumberOfBins(mitk::ScalarType min, mitk::ScalarType max) { return std::ceil( ( (max - min ) / m_HistogramBinSize) ); } template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator::InternalCalculateStatisticsUnmasked( const itk::Image< TPixel, VImageDimension > *image, StatisticsContainer *statisticsContainer, HistogramContainer* histogramContainer ) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef typename ImageType::IndexType IndexType; typedef itk::Statistics::ScalarImageToHistogramGenerator< ImageType > HistogramGeneratorType; statisticsContainer->clear(); histogramContainer->clear(); // Progress listening... typedef itk::SimpleMemberCommand< ImageStatisticsCalculator > ITKCommandType; ITKCommandType::Pointer progressListener; progressListener = ITKCommandType::New(); progressListener->SetCallbackFunction( this, &ImageStatisticsCalculator::UnmaskedStatisticsProgressUpdate ); // Issue 100 artificial progress events since ScalarIMageToHistogramGenerator // does not (yet?) support progress reporting this->InvokeEvent( itk::StartEvent() ); for ( unsigned int i = 0; i < 100; ++i ) { this->UnmaskedStatisticsProgressUpdate(); } // Calculate statistics (separate filter) typedef itk::ExtendedStatisticsImageFilter< ImageType > StatisticsFilterType; typename StatisticsFilterType::Pointer statisticsFilter = StatisticsFilterType::New(); statisticsFilter->SetInput( image ); statisticsFilter->SetBinSize( 100 ); unsigned long observerTag = statisticsFilter->AddObserver( itk::ProgressEvent(), progressListener ); try { statisticsFilter->Update(); } catch (const itk::ExceptionObject& e) { mitkThrow() << "Image statistics calculation failed due to following ITK Exception: \n " << e.what(); } catch( const std::exception& e ) { //mitkThrow() << "Image statistics calculation failed due to following ITK Exception: \n " << e.what(); } statisticsFilter->RemoveObserver( observerTag ); this->InvokeEvent( itk::EndEvent() ); // Calculate minimum and maximum typedef itk::MinimumMaximumImageCalculator< ImageType > MinMaxFilterType; typename MinMaxFilterType::Pointer minMaxFilter = MinMaxFilterType::New(); minMaxFilter->SetImage( image ); unsigned long observerTag2 = minMaxFilter->AddObserver( itk::ProgressEvent(), progressListener ); minMaxFilter->Compute(); minMaxFilter->RemoveObserver( observerTag2 ); this->InvokeEvent( itk::EndEvent() ); Statistics statistics; statistics.Reset(); statistics.SetLabel(1); statistics.SetN(image->GetBufferedRegion().GetNumberOfPixels()); statistics.SetMin(statisticsFilter->GetMinimum()); statistics.SetMax(statisticsFilter->GetMaximum()); statistics.SetMean(statisticsFilter->GetMean()); statistics.SetMedian(statisticsFilter->GetMedian()); statistics.SetVariance(statisticsFilter->GetVariance()); statistics.SetSkewness(statisticsFilter->GetSkewness()); statistics.SetKurtosis(statisticsFilter->GetKurtosis()); statistics.SetUniformity( statisticsFilter->GetUniformity()); statistics.SetEntropy( statisticsFilter->GetEntropy()); statistics.SetUPP( statisticsFilter->GetUPP()); statistics.SetMPP( statisticsFilter->GetMPP()); statistics.SetSigma(statisticsFilter->GetSigma()); statistics.SetRMS(sqrt( statistics.GetMean() * statistics.GetMean() + statistics.GetSigma() * statistics.GetSigma() )); statistics.GetMinIndex().set_size(image->GetImageDimension()); statistics.GetMaxIndex().set_size(image->GetImageDimension()); vnl_vector tmpMaxIndex; vnl_vector tmpMinIndex; tmpMaxIndex.set_size(image->GetImageDimension() ); tmpMinIndex.set_size(image->GetImageDimension() ); for (unsigned int i=0; iGetIndexOfMaximum()[i]; tmpMinIndex[i] = minMaxFilter->GetIndexOfMinimum()[i]; } statistics.SetMinIndex(tmpMaxIndex); statistics.SetMinIndex(tmpMinIndex); if( IsHotspotCalculated() && VImageDimension == 3 ) { typedef itk::Image< unsigned short, VImageDimension > MaskImageType; typename MaskImageType::Pointer nullMask; bool isHotspotDefined(false); Statistics hotspotStatistics = this->CalculateHotspotStatistics(image, nullMask.GetPointer(), m_HotspotRadiusInMM, isHotspotDefined, 0 ); if (isHotspotDefined) { statistics.SetHasHotspotStatistics(true); statistics.GetHotspotStatistics() = hotspotStatistics; } else { statistics.SetHasHotspotStatistics(false); } if(statistics.GetHotspotStatistics().HasHotspotStatistics() ) { MITK_DEBUG << "Hotspot statistics available"; statistics.SetHotspotIndex(hotspotStatistics.GetHotspotIndex()); } else { MITK_ERROR << "No hotspot statistics available!"; } } statisticsContainer->push_back( statistics ); // Calculate histogram // calculate bin size or number of bins unsigned int numberOfBins = 200; // default number of bins if (m_UseDefaultBinSize) { m_HistogramBinSize = std::ceil( (statistics.GetMax() - statistics.GetMin() + 1)/numberOfBins ); } else { numberOfBins = calcNumberOfBins(statistics.GetMin(), statistics.GetMax()); } typename HistogramGeneratorType::Pointer histogramGenerator = HistogramGeneratorType::New(); histogramGenerator->SetInput( image ); histogramGenerator->SetMarginalScale( 100 ); histogramGenerator->SetNumberOfBins( numberOfBins ); histogramGenerator->SetHistogramMin( statistics.GetMin() ); histogramGenerator->SetHistogramMax( statistics.GetMax() ); histogramGenerator->Compute(); histogramContainer->push_back( histogramGenerator->GetOutput() ); } template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator::InternalMaskIgnoredPixels( const itk::Image< TPixel, VImageDimension > *image, itk::Image< unsigned short, VImageDimension > *maskImage ) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< unsigned short, VImageDimension > MaskImageType; itk::ImageRegionIterator itmask(maskImage, maskImage->GetLargestPossibleRegion()); itk::ImageRegionConstIterator itimage(image, image->GetLargestPossibleRegion()); itmask.GoToBegin(); itimage.GoToBegin(); while( !itmask.IsAtEnd() ) { if(m_IgnorePixelValue == itimage.Get()) { itmask.Set(0); } ++itmask; ++itimage; } } template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator::InternalCalculateStatisticsMasked( const itk::Image< TPixel, VImageDimension > *image, itk::Image< unsigned short, VImageDimension > *maskImage, StatisticsContainer* statisticsContainer, HistogramContainer* histogramContainer ) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< unsigned short, VImageDimension > MaskImageType; typedef typename ImageType::IndexType IndexType; typedef typename ImageType::PointType PointType; typedef typename ImageType::SpacingType SpacingType; typedef typename ImageType::Pointer ImagePointer; typedef itk::ExtendedLabelStatisticsImageFilter< ImageType, MaskImageType > LabelStatisticsFilterType; typedef itk::ChangeInformationImageFilter< MaskImageType > ChangeInformationFilterType; typedef itk::ExtractImageFilter< ImageType, ImageType > ExtractImageFilterType; statisticsContainer->clear(); histogramContainer->clear(); // Make sure that mask is set if ( maskImage == nullptr ) { itkExceptionMacro( << "Mask image needs to be set!" ); } // Make sure that spacing of mask and image are the same //SpacingType imageSpacing = image->GetSpacing(); //SpacingType maskSpacing = maskImage->GetSpacing(); //PointType zeroPoint; zeroPoint.Fill( 0.0 ); //if ( (zeroPoint + imageSpacing).SquaredEuclideanDistanceTo( (zeroPoint + maskSpacing) ) > mitk::eps ) //{ // itkExceptionMacro( << "Mask needs to have same spacing as image! (Image spacing: " << imageSpacing << "; Mask spacing: " << maskSpacing << ")" ); //} // Make sure that orientation of mask and image are the same typedef typename ImageType::DirectionType DirectionType; 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 << ")" ); } } } } // 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(); long offset[ImageType::ImageDimension]; 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 ); if ( fabs( misalignment ) > mitk::eps ) { itkWarningMacro( << "Pixels/voxels of mask and image are not sufficiently aligned! (Misalignment: " << misalignment << ")" ); } double indexCoordDistance = maskOriginContinousIndex[i] - imageOriginContinousIndex[i]; offset[i] = int( indexCoordDistance + image->GetBufferedRegion().GetIndex()[i] + 0.5 ); } // Adapt the origin and region (index/size) of the mask so that the origin of both are the same typename ChangeInformationFilterType::Pointer adaptMaskFilter; adaptMaskFilter = ChangeInformationFilterType::New(); adaptMaskFilter->ChangeOriginOn(); adaptMaskFilter->ChangeRegionOn(); adaptMaskFilter->SetInput( maskImage ); adaptMaskFilter->SetOutputOrigin( image->GetOrigin() ); adaptMaskFilter->SetOutputOffset( offset ); typename MaskImageType::Pointer adaptedMaskImage; try { adaptMaskFilter->Update(); adaptedMaskImage = adaptMaskFilter->GetOutput(); } catch( const itk::ExceptionObject &e) { mitkThrow() << "Attempt to adapt shifted origin of the mask image failed due to ITK Exception: \n" << e.what(); } catch( const std::exception& e ) { //mitkThrow() << "Image statistics calculation failed due to following ITK Exception: \n " << e.what(); } // Make sure that mask region is contained within image region if ( adaptedMaskImage.IsNotNull() && !image->GetLargestPossibleRegion().IsInside( adaptedMaskImage->GetLargestPossibleRegion() ) ) { itkWarningMacro( << "Mask region needs to be inside of image region! (Image region: " << image->GetLargestPossibleRegion() << "; Mask region: " << adaptedMaskImage->GetLargestPossibleRegion() << ")" ); } // If mask region is smaller than image region, extract the sub-sampled region from the original image 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; } } typename ImageType::ConstPointer adaptedImage; if ( maskSmallerImage ) { typename ExtractImageFilterType::Pointer extractImageFilter = ExtractImageFilterType::New(); extractImageFilter->SetInput( image ); extractImageFilter->SetExtractionRegion( adaptedMaskImage->GetBufferedRegion() ); extractImageFilter->Update(); adaptedImage = extractImageFilter->GetOutput(); } else { adaptedImage = image; } // Initialize Filter typedef itk::StatisticsImageFilter< ImageType > StatisticsFilterType; typename StatisticsFilterType::Pointer statisticsFilter = StatisticsFilterType::New(); statisticsFilter->SetInput( adaptedImage ); try { statisticsFilter->Update(); } catch( const itk::ExceptionObject& e) { mitkThrow() << "Image statistics initialization computation failed with ITK Exception: \n " << e.what(); } catch( const std::exception& e ) { //mitkThrow() << "Image statistics calculation failed due to following ITK Exception: \n " << e.what(); } // Calculate bin size or number of bins unsigned int numberOfBins = 200; // default number of bins double maximum = 0.0; double minimum = 0.0; if (m_UseBinSizeBasedOnVOIRegion) { maximum = statisticsFilter->GetMaximum(); minimum = statisticsFilter->GetMinimum(); if (m_UseDefaultBinSize) { m_HistogramBinSize = std::ceil( static_cast((statisticsFilter->GetMaximum() - statisticsFilter->GetMinimum() + 1)/numberOfBins) ); } else { numberOfBins = calcNumberOfBins(statisticsFilter->GetMinimum(), statisticsFilter->GetMaximum()); } } else { double sig = 0.0; int counter = 0; //Find the min and max values for the Roi to set the range for the histogram GetMinAndMaxValue( minimum, maximum, counter, sig, image, maskImage); numberOfBins = maximum - minimum; if(maximum - minimum <= 10) { numberOfBins = 100; } } typename LabelStatisticsFilterType::Pointer labelStatisticsFilter = LabelStatisticsFilterType::New(); labelStatisticsFilter->SetInput( adaptedImage ); labelStatisticsFilter->SetLabelInput( adaptedMaskImage ); labelStatisticsFilter->UseHistogramsOn(); labelStatisticsFilter->SetHistogramParameters( numberOfBins, floor(minimum), ceil(maximum) ); //statisticsFilter->GetMinimum() statisticsFilter->GetMaximum() // Add progress listening typedef itk::SimpleMemberCommand< ImageStatisticsCalculator > ITKCommandType; ITKCommandType::Pointer progressListener; progressListener = ITKCommandType::New(); progressListener->SetCallbackFunction( this, &ImageStatisticsCalculator::MaskedStatisticsProgressUpdate ); unsigned long observerTag = labelStatisticsFilter->AddObserver( itk::ProgressEvent(), progressListener ); // Execute filter this->InvokeEvent( itk::StartEvent() ); // Make sure that only the mask region is considered (otherwise, if the mask region is smaller // than the image region, the Update() would result in an exception). labelStatisticsFilter->GetOutput()->SetRequestedRegion( adaptedMaskImage->GetLargestPossibleRegion() ); // Execute the filter try { labelStatisticsFilter->Update(); } catch( const itk::ExceptionObject& e) { mitkThrow() << "Image statistics calculation failed due to following ITK Exception: \n " << e.what(); } catch( const std::exception& e ) { //mitkThrow() << "Image statistics calculation failed due to following ITK Exception: \n " << e.what(); } this->InvokeEvent( itk::EndEvent() ); if( observerTag ) labelStatisticsFilter->RemoveObserver( observerTag ); // Find all relevant labels of mask (other than 0) std::list< int > relevantLabels = labelStatisticsFilter->GetRelevantLabels(); unsigned int i; if ( labelStatisticsFilter->GetMaskingNonEmpty() ) { std::list< int >::iterator it; for ( it = relevantLabels.begin(), i = 0; it != relevantLabels.end(); ++it, ++i ) { Statistics statistics; // restore previous code labelStatisticsFilter->GetHistogram(*it) ; histogramContainer->push_back( HistogramType::ConstPointer( labelStatisticsFilter->GetHistogram( (*it) ) ) ); statistics.SetLabel (*it); statistics.SetN(labelStatisticsFilter->GetCount( *it )); statistics.SetMin(labelStatisticsFilter->GetMinimum( *it )); statistics.SetMax(labelStatisticsFilter->GetMaximum( *it )); statistics.SetMean(labelStatisticsFilter->GetMean( *it )); statistics.SetMedian(labelStatisticsFilter->GetMedian( *it)); statistics.SetMedian(labelStatisticsFilter->GetMedian( *it )); statistics.SetVariance(labelStatisticsFilter->GetVariance( *it )); statistics.SetSigma(labelStatisticsFilter->GetSigma( *it )); statistics.SetSkewness(labelStatisticsFilter->GetSkewness( *it )); statistics.SetKurtosis(labelStatisticsFilter->GetKurtosis( *it )); statistics.SetUniformity( labelStatisticsFilter->GetUniformity( *it )); statistics.SetEntropy( labelStatisticsFilter->GetEntropy( *it )); statistics.SetUPP( labelStatisticsFilter->GetUPP( *it)); statistics.SetMPP( labelStatisticsFilter->GetMPP( *it)); statistics.SetRMS(sqrt( statistics.GetMean() * statistics.GetMean() + statistics.GetSigma() * statistics.GetSigma() )); // restrict image to mask area for min/max index calculation typedef itk::MaskImageFilter< ImageType, MaskImageType, ImageType > MaskImageFilterType; typename MaskImageFilterType::Pointer masker = MaskImageFilterType::New(); bool isMinAndMaxSameValue = (statistics.GetMin() == statistics.GetMax()); // bug 17962: following is a workaround for the case when min and max are the same, we can probably find a nicer way here double outsideValue = (isMinAndMaxSameValue ? (statistics.GetMax()/2) : (statistics.GetMin()+statistics.GetMax())/2); masker->SetOutsideValue( outsideValue ); masker->SetInput1(adaptedImage); masker->SetInput2(adaptedMaskImage); masker->Update(); // get index of minimum and maximum typedef itk::MinimumMaximumImageCalculator< ImageType > MinMaxFilterType; typename MinMaxFilterType::Pointer minMaxFilter = MinMaxFilterType::New(); minMaxFilter->SetImage( masker->GetOutput() ); unsigned long observerTag2 = minMaxFilter->AddObserver( itk::ProgressEvent(), progressListener ); minMaxFilter->Compute(); minMaxFilter->RemoveObserver( observerTag2 ); this->InvokeEvent( itk::EndEvent() ); typename MinMaxFilterType::IndexType tempMaxIndex = minMaxFilter->GetIndexOfMaximum(); // bug 17962: following is a workaround for the case when min and max are the same, we can probably find a nicer way here typename MinMaxFilterType::IndexType tempMinIndex = (isMinAndMaxSameValue ? minMaxFilter->GetIndexOfMaximum() : minMaxFilter->GetIndexOfMinimum()); // FIX BUG 14644 //If a PlanarFigure is used for segmentation the //adaptedImage is a single slice (2D). Adding the // 3. dimension. vnl_vector maxIndex; vnl_vector minIndex; maxIndex.set_size(m_Image->GetDimension()); minIndex.set_size(m_Image->GetDimension()); if (m_MaskingMode == MASKING_MODE_PLANARFIGURE && m_Image->GetDimension()==3) { maxIndex[m_PlanarFigureCoordinate0] = tempMaxIndex[0]; maxIndex[m_PlanarFigureCoordinate1] = tempMaxIndex[1]; maxIndex[m_PlanarFigureAxis] = m_PlanarFigureSlice; minIndex[m_PlanarFigureCoordinate0] = tempMinIndex[0] ; minIndex[m_PlanarFigureCoordinate1] = tempMinIndex[1]; minIndex[m_PlanarFigureAxis] = m_PlanarFigureSlice; } else { for (unsigned int i = 0; ipush_back( statistics ); } } else { histogramContainer->push_back( HistogramType::ConstPointer( m_EmptyHistogram ) ); statisticsContainer->push_back( Statistics() ); } } template ImageStatisticsCalculator::ImageExtrema ImageStatisticsCalculator::CalculateExtremaWorld( const itk::Image *inputImage, itk::Image *maskImage, double neccessaryDistanceToImageBorderInMM, unsigned int label) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< unsigned short, VImageDimension > MaskImageType; typedef itk::ImageRegionConstIteratorWithIndex MaskImageIteratorType; typedef itk::ImageRegionConstIteratorWithIndex InputImageIndexIteratorType; typename ImageType::SpacingType spacing = inputImage->GetSpacing(); ImageExtrema minMax; minMax.Defined = false; minMax.MaxIndex.set_size(VImageDimension); minMax.MaxIndex.set_size(VImageDimension); typename ImageType::RegionType allowedExtremaRegion = inputImage->GetLargestPossibleRegion(); bool keepDistanceToImageBorders( neccessaryDistanceToImageBorderInMM > 0 ); if (keepDistanceToImageBorders) { long distanceInPixels[VImageDimension]; for(unsigned short dimension = 0; dimension < VImageDimension; ++dimension) { // To confirm that the whole hotspot is inside the image we have to keep a specific distance to the image-borders, which is as long as // the radius. To get the amount of indices we divide the radius by spacing and add 0.5 because voxels are center based: // For example with a radius of 2.2 and a spacing of 1 two indices are enough because 2.2 / 1 + 0.5 = 2.7 => 2. // But with a radius of 2.7 we need 3 indices because 2.7 / 1 + 0.5 = 3.2 => 3 distanceInPixels[dimension] = int( neccessaryDistanceToImageBorderInMM / spacing[dimension] + 0.5); } allowedExtremaRegion.ShrinkByRadius(distanceInPixels); } InputImageIndexIteratorType imageIndexIt(inputImage, allowedExtremaRegion); float maxValue = itk::NumericTraits::min(); float minValue = itk::NumericTraits::max(); typename ImageType::IndexType maxIndex; typename ImageType::IndexType minIndex; for(unsigned short i = 0; i < VImageDimension; ++i) { maxIndex[i] = 0; minIndex[i] = 0; } if (maskImage != nullptr) { MaskImageIteratorType maskIt(maskImage, maskImage->GetLargestPossibleRegion()); typename ImageType::IndexType imageIndex; typename ImageType::PointType worldPosition; typename ImageType::IndexType maskIndex; for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt) { imageIndex = maskIndex = maskIt.GetIndex(); if(maskIt.Get() == label) { if( allowedExtremaRegion.IsInside(imageIndex) ) { imageIndexIt.SetIndex( imageIndex ); double value = imageIndexIt.Get(); minMax.Defined = true; //Calculate minimum, maximum and corresponding index-values if( value > maxValue ) { maxIndex = imageIndexIt.GetIndex(); maxValue = value; } if(value < minValue ) { minIndex = imageIndexIt.GetIndex(); minValue = value; } } } } } else { for(imageIndexIt.GoToBegin(); !imageIndexIt.IsAtEnd(); ++imageIndexIt) { double value = imageIndexIt.Get(); minMax.Defined = true; //Calculate minimum, maximum and corresponding index-values if( value > maxValue ) { maxIndex = imageIndexIt.GetIndex(); maxValue = value; } if(value < minValue ) { minIndex = imageIndexIt.GetIndex(); minValue = value; } } } minMax.MaxIndex.set_size(VImageDimension); minMax.MinIndex.set_size(VImageDimension); for(unsigned int i = 0; i < minMax.MaxIndex.size(); ++i) { minMax.MaxIndex[i] = maxIndex[i]; } for(unsigned int i = 0; i < minMax.MinIndex.size(); ++i) { minMax.MinIndex[i] = minIndex[i]; } minMax.Max = maxValue; minMax.Min = minValue; return minMax; } template itk::Size ImageStatisticsCalculator ::CalculateConvolutionKernelSize(double spacing[VImageDimension], double radiusInMM) { typedef itk::Image< float, VImageDimension > KernelImageType; typedef typename KernelImageType::SizeType SizeType; SizeType maskSize; for(unsigned int i = 0; i < VImageDimension; ++i) { maskSize[i] = static_cast( 2 * radiusInMM / spacing[i]); // We always want an uneven size to have a clear center point in the convolution mask if(maskSize[i] % 2 == 0 ) { ++maskSize[i]; } } return maskSize; } template itk::SmartPointer< itk::Image > ImageStatisticsCalculator ::GenerateHotspotSearchConvolutionKernel(double mmPerPixel[VImageDimension], double radiusInMM) { std::stringstream ss; for (unsigned int i = 0; i < VImageDimension; ++i) { ss << mmPerPixel[i]; if (i < VImageDimension -1) ss << ","; } MITK_DEBUG << "Update convolution kernel for spacing (" << ss.str() << ") and radius " << radiusInMM << "mm"; double radiusInMMSquared = radiusInMM * radiusInMM; typedef itk::Image< float, VImageDimension > KernelImageType; typename KernelImageType::Pointer convolutionKernel = KernelImageType::New(); // Calculate size and allocate mask image typedef typename KernelImageType::SizeType SizeType; SizeType maskSize = this->CalculateConvolutionKernelSize(mmPerPixel, radiusInMM); Point3D convolutionMaskCenterIndex; convolutionMaskCenterIndex.Fill(0.0); for(unsigned int i = 0; i < VImageDimension; ++i) { convolutionMaskCenterIndex[i] = 0.5 * (double)(maskSize[i]-1); } typedef typename KernelImageType::IndexType IndexType; IndexType maskIndex; maskIndex.Fill(0); typedef typename KernelImageType::RegionType RegionType; RegionType maskRegion; maskRegion.SetSize(maskSize); maskRegion.SetIndex(maskIndex); convolutionKernel->SetRegions(maskRegion); convolutionKernel->SetSpacing(mmPerPixel); convolutionKernel->Allocate(); // Fill mask image values by subsampling the image grid typedef itk::ImageRegionIteratorWithIndex MaskIteratorType; MaskIteratorType maskIt(convolutionKernel,maskRegion); int numberOfSubVoxelsPerDimension = 2; // per dimension! int numberOfSubVoxels = ::pow( static_cast(numberOfSubVoxelsPerDimension), static_cast(VImageDimension) ); double subVoxelSizeInPixels = 1.0 / (double)numberOfSubVoxelsPerDimension; double valueOfOneSubVoxel = 1.0 / (double)numberOfSubVoxels; double maskValue = 0.0; Point3D subVoxelIndexPosition; double distanceSquared = 0.0; typedef itk::ContinuousIndex ContinuousIndexType; for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt) { ContinuousIndexType indexPoint(maskIt.GetIndex()); Point3D voxelPosition; for (unsigned int dimension = 0; dimension < VImageDimension; ++dimension) { voxelPosition[dimension] = indexPoint[dimension]; } maskValue = 0.0; Vector3D subVoxelOffset; subVoxelOffset.Fill(0.0); // iterate sub-voxels by iterating all possible offsets for (subVoxelOffset[0] = -0.5 + subVoxelSizeInPixels / 2.0; subVoxelOffset[0] < +0.5; subVoxelOffset[0] += subVoxelSizeInPixels) { for (subVoxelOffset[1] = -0.5 + subVoxelSizeInPixels / 2.0; subVoxelOffset[1] < +0.5; subVoxelOffset[1] += subVoxelSizeInPixels) { for (subVoxelOffset[2] = -0.5 + subVoxelSizeInPixels / 2.0; subVoxelOffset[2] < +0.5; subVoxelOffset[2] += subVoxelSizeInPixels) { subVoxelIndexPosition = voxelPosition + subVoxelOffset; // this COULD be integrated into the for-loops if neccessary (add voxelPosition to initializer and end condition) distanceSquared = (subVoxelIndexPosition[0]-convolutionMaskCenterIndex[0]) * mmPerPixel[0] * (subVoxelIndexPosition[0]-convolutionMaskCenterIndex[0]) * mmPerPixel[0] + (subVoxelIndexPosition[1]-convolutionMaskCenterIndex[1]) * mmPerPixel[1] * (subVoxelIndexPosition[1]-convolutionMaskCenterIndex[1]) * mmPerPixel[1] + (subVoxelIndexPosition[2]-convolutionMaskCenterIndex[2]) * mmPerPixel[2] * (subVoxelIndexPosition[2]-convolutionMaskCenterIndex[2]) * mmPerPixel[2]; if (distanceSquared <= radiusInMMSquared) { maskValue += valueOfOneSubVoxel; } } } } maskIt.Set( maskValue ); } return convolutionKernel; } template itk::SmartPointer > ImageStatisticsCalculator::GenerateConvolutionImage( const itk::Image* inputImage ) { double mmPerPixel[VImageDimension]; for (unsigned int dimension = 0; dimension < VImageDimension; ++dimension) { mmPerPixel[dimension] = inputImage->GetSpacing()[dimension]; } // update convolution kernel typedef itk::Image< float, VImageDimension > KernelImageType; typename KernelImageType::Pointer convolutionKernel = this->GenerateHotspotSearchConvolutionKernel(mmPerPixel, m_HotspotRadiusInMM); // update convolution image typedef itk::Image< TPixel, VImageDimension > InputImageType; typedef itk::Image< TPixel, VImageDimension > ConvolutionImageType; typedef itk::FFTConvolutionImageFilter ConvolutionFilterType; typename ConvolutionFilterType::Pointer convolutionFilter = ConvolutionFilterType::New(); typedef itk::ConstantBoundaryCondition BoundaryConditionType; BoundaryConditionType boundaryCondition; boundaryCondition.SetConstant(0.0); if (GetHotspotMustBeCompletlyInsideImage()) { // overwrite default boundary condition convolutionFilter->SetBoundaryCondition(&boundaryCondition); } convolutionFilter->SetInput(inputImage); convolutionFilter->SetKernelImage(convolutionKernel); convolutionFilter->SetNormalize(true); MITK_DEBUG << "Update Convolution image for hotspot search"; convolutionFilter->UpdateLargestPossibleRegion(); typename ConvolutionImageType::Pointer convolutionImage = convolutionFilter->GetOutput(); convolutionImage->SetSpacing( inputImage->GetSpacing() ); // only workaround because convolution filter seems to ignore spacing of input image m_HotspotRadiusInMMChanged = false; return convolutionImage; } template < typename TPixel, unsigned int VImageDimension> void ImageStatisticsCalculator ::FillHotspotMaskPixels( itk::Image* maskImage, itk::Point sphereCenter, double sphereRadiusInMM) { typedef itk::Image< TPixel, VImageDimension > MaskImageType; typedef itk::ImageRegionIteratorWithIndex MaskImageIteratorType; MaskImageIteratorType maskIt(maskImage, maskImage->GetLargestPossibleRegion()); typename MaskImageType::IndexType maskIndex; typename MaskImageType::PointType worldPosition; for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt) { maskIndex = maskIt.GetIndex(); maskImage->TransformIndexToPhysicalPoint(maskIndex, worldPosition); maskIt.Set( worldPosition.EuclideanDistanceTo(sphereCenter) <= sphereRadiusInMM ? 1 : 0 ); } } template < typename TPixel, unsigned int VImageDimension> ImageStatisticsCalculator::Statistics ImageStatisticsCalculator::CalculateHotspotStatistics( const itk::Image* inputImage, itk::Image* maskImage, double radiusInMM, bool& isHotspotDefined, unsigned int label) { // get convolution image (updated in GenerateConvolutionImage()) typedef itk::Image< TPixel, VImageDimension > InputImageType; typedef itk::Image< TPixel, VImageDimension > ConvolutionImageType; typedef itk::Image< float, VImageDimension > KernelImageType; typedef itk::Image< unsigned short, VImageDimension > MaskImageType; //typename ConvolutionImageType::Pointer convolutionImage = dynamic_cast(this->GenerateConvolutionImage(inputImage)); typename ConvolutionImageType::Pointer convolutionImage = this->GenerateConvolutionImage(inputImage); if (convolutionImage.IsNull()) { MITK_ERROR << "Empty convolution image in CalculateHotspotStatistics(). We should never reach this state (logic error)."; throw std::logic_error("Empty convolution image in CalculateHotspotStatistics()"); } // find maximum in convolution image, given the current mask double requiredDistanceToBorder = m_HotspotMustBeCompletelyInsideImage ? m_HotspotRadiusInMM : -1.0; ImageExtrema convolutionImageInformation = CalculateExtremaWorld(convolutionImage.GetPointer(), maskImage, requiredDistanceToBorder, label); isHotspotDefined = convolutionImageInformation.Defined; if (!isHotspotDefined) { m_EmptyStatistics.Reset(VImageDimension); MITK_ERROR << "No origin of hotspot-sphere was calculated! Returning empty statistics"; return m_EmptyStatistics; } else { // create a binary mask around the "hotspot" region, fill the shape of a sphere around our hotspot center typedef itk::ImageDuplicator< InputImageType > DuplicatorType; typename DuplicatorType::Pointer copyMachine = DuplicatorType::New(); copyMachine->SetInputImage(inputImage); copyMachine->Update(); typedef itk::CastImageFilter< InputImageType, MaskImageType > CastFilterType; typename CastFilterType::Pointer caster = CastFilterType::New(); caster->SetInput( copyMachine->GetOutput() ); caster->Update(); typename MaskImageType::Pointer hotspotMaskITK = caster->GetOutput(); typedef typename InputImageType::IndexType IndexType; IndexType maskCenterIndex; for (unsigned int d =0; d< VImageDimension;++d) maskCenterIndex[d]=convolutionImageInformation.MaxIndex[d]; typename ConvolutionImageType::PointType maskCenter; inputImage->TransformIndexToPhysicalPoint(maskCenterIndex,maskCenter); this->FillHotspotMaskPixels(hotspotMaskITK.GetPointer(), maskCenter, radiusInMM); // calculate statistics within the binary mask typedef itk::ExtendedLabelStatisticsImageFilter< InputImageType, MaskImageType> LabelStatisticsFilterType; typename LabelStatisticsFilterType::Pointer labelStatisticsFilter; labelStatisticsFilter = LabelStatisticsFilterType::New(); labelStatisticsFilter->SetInput( inputImage ); labelStatisticsFilter->SetLabelInput( hotspotMaskITK ); labelStatisticsFilter->Update(); Statistics hotspotStatistics; hotspotStatistics.SetHotspotIndex(convolutionImageInformation.MaxIndex); hotspotStatistics.SetMean(convolutionImageInformation.Max); if ( labelStatisticsFilter->HasLabel( 1 ) ) { hotspotStatistics.SetLabel (1); hotspotStatistics.SetN(labelStatisticsFilter->GetCount(1)); hotspotStatistics.SetMin(labelStatisticsFilter->GetMinimum(1)); hotspotStatistics.SetMax(labelStatisticsFilter->GetMaximum(1)); hotspotStatistics.SetMedian(labelStatisticsFilter->GetMedian(1)); hotspotStatistics.SetVariance(labelStatisticsFilter->GetVariance(1)); hotspotStatistics.SetSigma(labelStatisticsFilter->GetSigma(1)); hotspotStatistics.SetRMS(sqrt( hotspotStatistics.GetMean() * hotspotStatistics.GetMean() + hotspotStatistics.GetSigma() * hotspotStatistics.GetSigma() )); MITK_DEBUG << "Statistics for inside hotspot: Mean " << hotspotStatistics.GetMean() << ", SD " << hotspotStatistics.GetSigma() << ", Max " << hotspotStatistics.GetMax() << ", Min " << hotspotStatistics.GetMin(); } else { MITK_ERROR << "Uh oh! Unable to calculate statistics for hotspot region..."; return m_EmptyStatistics; } return hotspotStatistics; } } template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator::InternalCalculateMaskFromPlanarFigure( const itk::Image< TPixel, VImageDimension > *image, unsigned int axis ) { typedef itk::Image< TPixel, VImageDimension > ImageType; 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 ); // 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_Image->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; } m_PlanarFigureCoordinate0= i0; m_PlanarFigureCoordinate1= i1; // 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; } + // Fabian: Why convert to index coordinates? 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( 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_InternalImageMask2D = duplicator->GetOutput(); } void ImageStatisticsCalculator::UnmaskedStatisticsProgressUpdate() { // Need to throw away every second progress event to reach a final count of // 100 since two consecutive filters are used in this case static int updateCounter = 0; if ( updateCounter++ % 2 == 0 ) { this->InvokeEvent( itk::ProgressEvent() ); } } void ImageStatisticsCalculator::MaskedStatisticsProgressUpdate() { this->InvokeEvent( itk::ProgressEvent() ); } } diff --git a/Modules/ImageStatistics/mitkImageStatisticsCalculator2.cpp b/Modules/ImageStatistics/mitkImageStatisticsCalculator2.cpp new file mode 100644 index 0000000000..10e5f1e4cb --- /dev/null +++ b/Modules/ImageStatistics/mitkImageStatisticsCalculator2.cpp @@ -0,0 +1,5 @@ +#ifndef MITKIMAGESTATISTICSCALCULATOR2 +#define MITKIMAGESTATISTICSCALCULATOR2 + +#endif // MITKIMAGESTATISTICSCALCULATOR2 + diff --git a/Modules/ImageStatistics/mitkImageStatisticsCalculator2.h b/Modules/ImageStatistics/mitkImageStatisticsCalculator2.h new file mode 100644 index 0000000000..10e5f1e4cb --- /dev/null +++ b/Modules/ImageStatistics/mitkImageStatisticsCalculator2.h @@ -0,0 +1,5 @@ +#ifndef MITKIMAGESTATISTICSCALCULATOR2 +#define MITKIMAGESTATISTICSCALCULATOR2 + +#endif // MITKIMAGESTATISTICSCALCULATOR2 + diff --git a/Modules/ImageStatistics/mitkMaskGenerator.cpp b/Modules/ImageStatistics/mitkMaskGenerator.cpp new file mode 100644 index 0000000000..6f34324163 --- /dev/null +++ b/Modules/ImageStatistics/mitkMaskGenerator.cpp @@ -0,0 +1,33 @@ +#include + +namespace mitk +{ +void MaskGenerator::SetTimeStep(unsigned int timeStep) +{ + if (m_TimeStep != timeStep) + { + m_TimeStep = timeStep; + m_Modified = true; + } +} + +MaskGenerator::MaskGenerator(): + m_TimeStep(0), + m_Modified(false) +{ + +} + +//typename itk::Region<3>::Pointer MaskGenerator::GetImageRegionOfMask(Image::Pointer image) +//{ +// if (m_InternalMask.IsNull() || m_Modified) +// { +// MITK_ERROR << "Update MaskGenerator first!"; +// } + +// mitk::BaseGeometry::Pointer imageGeometry = image->GetGeometry(); +// mitk::BaseGeometry::Pointer maskGeometry = m_InternalMask->GetGeometry(); + + +//} +} diff --git a/Modules/ImageStatistics/mitkMaskGenerator.h b/Modules/ImageStatistics/mitkMaskGenerator.h new file mode 100644 index 0000000000..7c19bf6b6b --- /dev/null +++ b/Modules/ImageStatistics/mitkMaskGenerator.h @@ -0,0 +1,33 @@ +#ifndef MITKMASKGENERATOR +#define MITKMASKGENERATOR + +#include +#include +#include + +namespace mitk +{ +class MITKIMAGESTATISTICS_EXPORT MaskGenerator +{ +public: + MaskGenerator(); + + //~MaskGenerator(); + + virtual mitk::Image::Pointer GetMask() = 0; + + void SetTimeStep(unsigned int timeStep); + +protected: + unsigned int m_TimeStep; + mitk::Image::Pointer m_InternalMask; + bool m_Modified; + +private: + void CheckMaskSanity(mitk::Image::Pointer); + +}; +} + +#endif // MITKMASKGENERATOR + diff --git a/Modules/ImageStatistics/mitkMultiLabelMaskGenerator.cpp b/Modules/ImageStatistics/mitkMultiLabelMaskGenerator.cpp new file mode 100644 index 0000000000..e3d467c857 --- /dev/null +++ b/Modules/ImageStatistics/mitkMultiLabelMaskGenerator.cpp @@ -0,0 +1 @@ +#include diff --git a/Modules/ImageStatistics/mitkMultiLabelMaskGenerator.h b/Modules/ImageStatistics/mitkMultiLabelMaskGenerator.h new file mode 100644 index 0000000000..4a0a0be26f --- /dev/null +++ b/Modules/ImageStatistics/mitkMultiLabelMaskGenerator.h @@ -0,0 +1,25 @@ +#ifndef MITKMULTILABELMASKGENERATOR +#define MITKMULTILABELMASKGENERATOR + +#include +#include +#include +#include + +namespace mitk +{ + +class MITKIMAGESTATISTICS_EXPORT MultiLabelMaskGenerator: public MaskGenerator +{ +public: + void SetLabelSetImage(mitk::LabelSetImage::Pointer labelSetImage); + +protected: + +private: + +}; + +} + +#endif diff --git a/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.cpp b/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.cpp new file mode 100644 index 0000000000..2df2feb9e8 --- /dev/null +++ b/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.cpp @@ -0,0 +1,388 @@ +#include +#include +#include +#include "mitkImageAccessByItk.h" +#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() +{ + 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(); + mitk::Image::Pointer imageSlice = mitk::Image::New(); + + if (dimension == 3) + { + ExtractImageFilter::Pointer imageExtractor = ExtractImageFilter::New(); + imageExtractor->SetInput( m_InternalInputImage ); + imageExtractor->SetSliceDimension( axis ); + imageExtractor->SetSliceIndex( slice ); + imageExtractor->Update(); + imageSlice = imageExtractor->GetOutput(); + } + else if(dimension == 2) + { + imageSlice = m_InternalInputImage; + } + 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 new file mode 100644 index 0000000000..38c42f8cac --- /dev/null +++ b/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.h @@ -0,0 +1,95 @@ +#ifndef MITKPLANARFIGUREMASKGENERATOR +#define MITKPLANARFIGUREMASKGENERATOR + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace mitk +{ +class MITKIMAGESTATISTICS_EXPORT PlanarFigureMaskGenerator: public MaskGenerator + { + public: + + 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: + + + 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; + + unsigned int m_PlanarFigureAxis; + }; +} + +#endif // MITKPLANARFIGUREMASKGENERATOR + diff --git a/Modules/SurfaceInterpolation/mitkReduceContourSetFilter.cpp b/Modules/SurfaceInterpolation/mitkReduceContourSetFilter.cpp index 13119b9af1..cdb59bb323 100644 --- a/Modules/SurfaceInterpolation/mitkReduceContourSetFilter.cpp +++ b/Modules/SurfaceInterpolation/mitkReduceContourSetFilter.cpp @@ -1,502 +1,503 @@ /*=================================================================== 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 "mitkReduceContourSetFilter.h" mitk::ReduceContourSetFilter::ReduceContourSetFilter() { m_MaxSegmentLenght = 0; m_StepSize = 10; m_Tolerance = -1; m_ReductionType = DOUGLAS_PEUCKER; m_MaxSpacing = -1; m_MinSpacing = -1; this->m_UseProgressBar = false; this->m_ProgressStepSize = 1; m_NumberOfPointsAfterReduction = 0; mitk::Surface::Pointer output = mitk::Surface::New(); this->SetNthOutput(0, output.GetPointer()); } mitk::ReduceContourSetFilter::~ReduceContourSetFilter() { } void mitk::ReduceContourSetFilter::GenerateData() { unsigned int numberOfInputs = this->GetNumberOfIndexedInputs(); unsigned int numberOfOutputs (0); vtkSmartPointer newPolyData; vtkSmartPointer newPolygons; vtkSmartPointer newPoints; //For the purpose of evaluation // unsigned int numberOfPointsBefore (0); m_NumberOfPointsAfterReduction=0; for(unsigned int i = 0; i < numberOfInputs; i++) { mitk::Surface* currentSurface = const_cast( this->GetInput(i) ); vtkSmartPointer polyData = currentSurface->GetVtkPolyData(); newPolyData = vtkSmartPointer::New(); newPolygons = vtkSmartPointer::New(); newPoints = vtkSmartPointer::New(); vtkSmartPointer existingPolys = polyData->GetPolys(); vtkSmartPointer existingPoints = polyData->GetPoints(); existingPolys->InitTraversal(); vtkIdType* cell (nullptr); vtkIdType cellSize (0); for( existingPolys->InitTraversal(); existingPolys->GetNextCell(cellSize, cell);) { bool incorporatePolygon = this->CheckForIntersection(cell,cellSize,existingPoints, /*numberOfIntersections, intersectionPoints, */i); if ( !incorporatePolygon ) continue; vtkSmartPointer newPolygon = vtkSmartPointer::New(); if(m_ReductionType == NTH_POINT) { this->ReduceNumberOfPointsByNthPoint(cellSize, cell, existingPoints, newPolygon, newPoints); if (newPolygon->GetPointIds()->GetNumberOfIds() != 0) { newPolygons->InsertNextCell(newPolygon); } } else if (m_ReductionType == DOUGLAS_PEUCKER) { this->ReduceNumberOfPointsByDouglasPeucker(cellSize, cell, existingPoints, newPolygon, newPoints); if (newPolygon->GetPointIds()->GetNumberOfIds() > 3) { newPolygons->InsertNextCell(newPolygon); } } //Again for evaluation // numberOfPointsBefore += cellSize; m_NumberOfPointsAfterReduction += newPolygon->GetPointIds()->GetNumberOfIds(); } if (newPolygons->GetNumberOfCells() != 0) { newPolyData->SetPolys(newPolygons); newPolyData->SetPoints(newPoints); newPolyData->BuildLinks(); this->SetNumberOfIndexedOutputs(numberOfOutputs + 1); mitk::Surface::Pointer surface = mitk::Surface::New(); this->SetNthOutput(numberOfOutputs, surface.GetPointer()); surface->SetVtkPolyData(newPolyData); numberOfOutputs++; } } // MITK_INFO<<"Points before: "<SetNumberOfIndexedOutputs(numberOfOutputs); //Setting progressbar if (this->m_UseProgressBar) mitk::ProgressBar::GetInstance()->Progress(this->m_ProgressStepSize); } void mitk::ReduceContourSetFilter::ReduceNumberOfPointsByNthPoint (vtkIdType cellSize, vtkIdType* cell, vtkPoints* points, vtkPolygon* reducedPolygon, vtkPoints* reducedPoints) { unsigned int newNumberOfPoints (0); unsigned int mod = cellSize%m_StepSize; if(mod == 0) { newNumberOfPoints = cellSize/m_StepSize; } else { newNumberOfPoints = ( (cellSize-mod)/m_StepSize )+1; } if (newNumberOfPoints <= 3) { return; } reducedPolygon->GetPointIds()->SetNumberOfIds(newNumberOfPoints); reducedPolygon->GetPoints()->SetNumberOfPoints(newNumberOfPoints); for (vtkIdType i = 0; i < cellSize; i++) { if (i%m_StepSize == 0) { double point[3]; points->GetPoint(cell[i], point); vtkIdType id = reducedPoints->InsertNextPoint(point); reducedPolygon->GetPointIds()->SetId(i/m_StepSize, id); } } vtkIdType id = cell[0]; double point[3]; points->GetPoint(id, point); id = reducedPoints->InsertNextPoint(point); reducedPolygon->GetPointIds()->SetId(newNumberOfPoints-1, id); } void mitk::ReduceContourSetFilter::ReduceNumberOfPointsByDouglasPeucker(vtkIdType cellSize, vtkIdType* cell, vtkPoints* points, vtkPolygon* reducedPolygon, vtkPoints* reducedPoints) { //If the cell is too small to obtain a reduced polygon with the given stepsize return if (cellSize <= static_cast(m_StepSize*3))return; /* What we do now is (see the Douglas Peucker Algorithm): 1. Divide the current contour in two line segments (start - middle; middle - end), put them into the stack 2. Fetch first line segment and create the following vectors: - v1 = (start;end) - v2 = (start;currentPoint) -> for each point of the current line segment! 3. Calculate the distance from the currentPoint to v1: a. Determine the length of the orthogonal projection of v2 to v1 by: l = v2 * (normalized v1) b. There a three possibilities for the distance then: d = sqrt(lenght(v2)^2 - l^2) if l > 0 and l < length(v1) d = lenght(v2-v1) if l > 0 and l > lenght(v1) d = length(v2) if l < 0 because v2 is then pointing in a different direction than v1 4. Memorize the point with the biggest distance and create two new line segments with it at the end of the iteration and put it into the stack 5. If the distance value D <= m_Tolerance, then add the start and end index and the corresponding points to the reduced ones */ //First of all set tolerance if none is specified if(m_Tolerance < 0) { if(m_MaxSpacing > 0) { m_Tolerance = m_MinSpacing; } else { m_Tolerance = 1.5; } } std::stack lineSegments; //1. Divide in line segments LineSegment ls2; ls2.StartIndex = cell[cellSize/2]; ls2.EndIndex = cell[cellSize-1]; lineSegments.push(ls2); LineSegment ls1; ls1.StartIndex = cell[0]; ls1.EndIndex = cell[cellSize/2]; lineSegments.push(ls1); LineSegment currentSegment; double v1[3]; double v2[3]; double tempV[3]; double lenghtV1; double currentMaxDistance (0); vtkIdType currentMaxDistanceIndex (0); double l; double d; vtkIdType pointId (0); //Add the start index to the reduced points. From now on just the end indices will be added pointId = reducedPoints->InsertNextPoint(points->GetPoint(cell[0])); reducedPolygon->GetPointIds()->InsertNextId(pointId); while (!lineSegments.empty()) { currentSegment = lineSegments.top(); lineSegments.pop(); //2. Create vectors points->GetPoint(currentSegment.EndIndex, tempV); points->GetPoint(currentSegment.StartIndex, v1); v1[0] = tempV[0]-v1[0]; v1[1] = tempV[1]-v1[1]; v1[2] = tempV[2]-v1[2]; lenghtV1 = vtkMath::Norm(v1); vtkMath::Normalize(v1); int range = currentSegment.EndIndex - currentSegment.StartIndex; for (int i = 1; i < abs(range); ++i) { points->GetPoint(currentSegment.StartIndex+i, tempV); points->GetPoint(currentSegment.StartIndex, v2); v2[0] = tempV[0]-v2[0]; v2[1] = tempV[1]-v2[1]; v2[2] = tempV[2]-v2[2]; //3. Calculate the distance l = vtkMath::Dot(v2, v1); d = vtkMath::Norm(v2); if (l > 0 && l < lenghtV1) { d = sqrt((d*d-l*l)); } else if (l > 0 && l > lenghtV1) { tempV[0] = lenghtV1*v1[0] - v2[0]; tempV[1] = lenghtV1*v1[1] - v2[1]; tempV[2] = lenghtV1*v1[2] - v2[2]; d = vtkMath::Norm(tempV); } //4. Memorize maximum distance if (d > currentMaxDistance) { currentMaxDistance = d; currentMaxDistanceIndex = currentSegment.StartIndex+i; } } //4. & 5. if (currentMaxDistance <= m_Tolerance) { //double temp[3]; int segmentLenght = currentSegment.EndIndex - currentSegment.StartIndex; if (segmentLenght > (int)m_MaxSegmentLenght) { m_MaxSegmentLenght = (unsigned int)segmentLenght; } // MITK_INFO<<"Lenght: "< 25) { unsigned int newLenght(segmentLenght); while (newLenght > 25) { newLenght = newLenght*0.5; } unsigned int divisions = abs(segmentLenght)/newLenght; // MITK_INFO<<"Divisions: "<InsertNextPoint(points->GetPoint(currentSegment.StartIndex + newLenght*i)); reducedPolygon->GetPointIds()->InsertNextId(pointId); } } // MITK_INFO<<"Inserting END: "<InsertNextPoint(points->GetPoint(currentSegment.EndIndex)); reducedPolygon->GetPointIds()->InsertNextId(pointId); } else { ls2.StartIndex = currentMaxDistanceIndex; ls2.EndIndex = currentSegment.EndIndex; lineSegments.push(ls2); ls1.StartIndex = currentSegment.StartIndex; ls1.EndIndex = currentMaxDistanceIndex; lineSegments.push(ls1); } currentMaxDistance = 0; } } bool mitk::ReduceContourSetFilter::CheckForIntersection (vtkIdType* currentCell, vtkIdType currentCellSize, vtkPoints* currentPoints,/* vtkIdType numberOfIntersections, vtkIdType* intersectionPoints,*/ unsigned int currentInputIndex) { /* If we check the current cell for intersections then we have to consider three possibilies: 1. There is another cell among all the other input surfaces which intersects the current polygon: - That means we have to save the intersection points because these points should not be eliminated 2. There current polygon exists just because of an intersection of another polygon with the current plane defined by the current polygon - That means the current polygon should not be incorporated and all of its points should be eliminated 3. There is no intersection - That mean we can just reduce the current polygons points without considering any intersections */ for (unsigned int i = 0; i < this->GetNumberOfIndexedInputs(); i++) { //Don't check for intersection with the polygon itself if (i == currentInputIndex) continue; //Get the next polydata to check for intersection vtkSmartPointer poly = const_cast( this->GetInput(i) )->GetVtkPolyData(); vtkSmartPointer polygonArray = poly->GetPolys(); polygonArray->InitTraversal(); vtkIdType anotherInputPolygonSize (0); vtkIdType* anotherInputPolygonIDs(nullptr); /* The procedure is: - Create the equation of the plane, defined by the points of next input - Calculate the distance of each point of the current polygon to the plane - If the maximum distance is not bigger than 1.5 of the maximum spacing AND the minimal distance is not bigger than 0.5 of the minimum spacing then the current contour is an intersection contour */ for( polygonArray->InitTraversal(); polygonArray->GetNextCell(anotherInputPolygonSize, anotherInputPolygonIDs);) { //Choosing three plane points to calculate the plane vectors double p1[3]; double p2[3]; double p3[3]; //The plane vectors double v1[3]; double v2[3] = { 0 }; //The plane normal double normal[3]; //Create first Vector poly->GetPoint(anotherInputPolygonIDs[0], p1); poly->GetPoint(anotherInputPolygonIDs[1], p2); v1[0] = p2[0]-p1[0]; v1[1] = p2[1]-p1[1]; v1[2] = p2[2]-p1[2]; //Find 3rd point for 2nd vector (The angle between the two plane vectors should be bigger than 30 degrees) double maxDistance (0); double minDistance (10000); for (vtkIdType j = 2; j < anotherInputPolygonSize; j++) { poly->GetPoint(anotherInputPolygonIDs[j], p3); v2[0] = p3[0]-p1[0]; v2[1] = p3[1]-p1[1]; v2[2] = p3[2]-p1[2]; //Calculate the angle between the two vector for the current point double dotV1V2 = vtkMath::Dot(v1,v2); double absV1 = sqrt(vtkMath::Dot(v1,v1)); double absV2 = sqrt(vtkMath::Dot(v2,v2)); double cosV1V2 = dotV1V2/(absV1*absV2); double arccos = acos(cosV1V2); double degree = vtkMath::DegreesFromRadians(arccos); //If angle is bigger than 30 degrees break if (degree > 30) break; }//for (to find 3rd point) //Calculate normal of the plane by taking the cross product of the two vectors vtkMath::Cross(v1,v2,normal); vtkMath::Normalize(normal); //Determine position of the plane double lambda = vtkMath::Dot(normal, p1); /* Calculate the distance to the plane for each point of the current polygon If the distance is zero then save the currentPoint as intersection point */ for (vtkIdType k = 0; k < currentCellSize; k++) { double currentPoint[3]; currentPoints->GetPoint(currentCell[k], currentPoint); double tempPoint[3]; tempPoint[0] = normal[0]*currentPoint[0]; tempPoint[1] = normal[1]*currentPoint[1]; tempPoint[2] = normal[2]*currentPoint[2]; double temp = tempPoint[0]+tempPoint[1]+tempPoint[2]-lambda; double distance = fabs(temp); if (distance > maxDistance) { maxDistance = distance; } if (distance < minDistance) { minDistance = distance; } }//for (to calculate distance and intersections with currentPolygon) if (maxDistance < 1.5*m_MaxSpacing && minDistance < 0.5*m_MinSpacing) { return false; } //Because we are considering the plane defined by the acual input polygon only one iteration is sufficient //We do not need to consider each cell of the plane break; }//for (to traverse through all cells of actualInputPolyData) }//for (to iterate through all inputs) return true; } void mitk::ReduceContourSetFilter::GenerateOutputInformation() { Superclass::GenerateOutputInformation(); } void mitk::ReduceContourSetFilter::Reset() { for (unsigned int i = 0; i < this->GetNumberOfIndexedInputs(); i++) { this->PopBackInput(); } this->SetNumberOfIndexedInputs(0); this->SetNumberOfIndexedOutputs(0); + // BUG XXXXX Fix mitk::Surface::Pointer output = mitk::Surface::New(); this->SetNthOutput(0, output.GetPointer()); m_NumberOfPointsAfterReduction = 0; } void mitk::ReduceContourSetFilter::SetUseProgressBar(bool status) { this->m_UseProgressBar = status; } void mitk::ReduceContourSetFilter::SetProgressStepSize(unsigned int stepSize) { this->m_ProgressStepSize = stepSize; }