diff --git a/Modules/ImageStatistics/Testing/mitkMultiGaussianTest.cpp b/Modules/ImageStatistics/Testing/mitkMultiGaussianTest.cpp index 18992600fa..dd0787010e 100644 --- a/Modules/ImageStatistics/Testing/mitkMultiGaussianTest.cpp +++ b/Modules/ImageStatistics/Testing/mitkMultiGaussianTest.cpp @@ -1,308 +1,93 @@ /*=================================================================== 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 "mitkTestingMacros.h" #include "itkMultiGaussianImageSource.h" -#include + #include #include -#include -#include -#include -#include - -//bin\Release\ImageStatisticsTestDriver.exe mitkMultiGaussianTest C:/temp/TestImage/image 1 12 12 10 1 1 5 20 200 "2.5" "2.5" 3 -int mitkMultiGaussianTest(int argc, char* argv[]) +int mitkMultiGaussianTest(int, char* []) { - if ( argc !=14 ) - { - std::cerr << " 14 arguments expected: [outputFilename, numberOfImages, imageSizeX, imageSizeY, imageSizeZ, numberOfGaussians, minWidthOfGaussian, maxWidthOfGaussian, minAltitudeOfGaussian, maxAltitudeOfGaussian, spacingX, spacingY, spacingZ ] " << std::endl; - return EXIT_FAILURE; - } // always start with this! - MITK_TEST_BEGIN("mitkMultiGaussianTest"); - MITK_TEST_CONDITION_REQUIRED(argc == 14, "Test called with 14 parameters"); - - const unsigned int Dimension = 3; - typedef double PixelType; - typedef itk::DOMNode::Pointer DOMNodeType; - typedef itk::Image ImageType; - typedef itk::MultiGaussianImageSource< ImageType > MultiGaussianImageSource; - std::string outputFilename = argv[1], name; - int numberOfImages; - double centerX, centerY, centerZ, sigmaX, sigmaY, sigmaZ, altitude; - unsigned int numberOfGaussians, minWidthOfGaussian, maxWidthOfGaussian, minAltitudeOfGaussian, maxAltitudeOfGaussian; - itk::MultiGaussianImageSource< ImageType >::VectorType centerXVec, centerYVec, centerZVec, sigmaXVec, sigmaYVec, sigmaZVec, altitudeVec; - itk::MultiGaussianImageSource< ImageType >::ItkVectorType regionOfInterestMax, regionOfInterestMin; - itk::MultiGaussianImageSource< ImageType >::IndexType sphereMidpt, maxValueIndexInSphere, minValueIndexInSphere; - itk::MultiGaussianImageSource< ImageType >::ItkVectorType X, Y, Z, Sig, Alt; - MultiGaussianImageSource::Pointer gaussianGenerator; - itk::DOMNodeXMLWriter::Pointer xmlWriter; - itk::MultiGaussianImageSource< ImageType >::SpacingValueArrayType spacing; - DOMNodeType domTestCase, domTestImage, domGaussian, domStatistics, domROI; - ImageType::SizeValueType size[3]; - std::stringstream ss; - double radius = pow(itk::Math::one_over_pi * 0.75 , 1.0 / 3.0) * 10; - char * fileNamePointer; - - if ( ! (std::istringstream(argv[2]) >> numberOfImages) ) numberOfImages = 0; - if ( ! (std::istringstream(argv[3]) >> size[0]) ) size[0] = 0; - if ( ! (std::istringstream(argv[4]) >> size[1]) ) size[1] = 0; - if ( ! (std::istringstream(argv[5]) >> size[2]) ) size[2] = 0; - if ( ! (std::istringstream(argv[6]) >> numberOfGaussians) ) numberOfGaussians = 0; - if ( ! (std::istringstream(argv[7]) >> minWidthOfGaussian) ) minWidthOfGaussian = (size[0] + size[1] + size[2]) / 27; - if ( ! (std::istringstream(argv[8]) >> maxWidthOfGaussian) ) maxWidthOfGaussian = (size[0] + size[1] + size[2]) / 9; - if ( ! (std::istringstream(argv[9]) >> minAltitudeOfGaussian) ) minAltitudeOfGaussian = 5; - if ( ! (std::istringstream(argv[10]) >> maxAltitudeOfGaussian) ) maxAltitudeOfGaussian = 200; - if ( ! (std::istringstream(argv[11]) >> spacing[0]) ) spacing[0] = 1; - if ( ! (std::istringstream(argv[12]) >> spacing[1]) ) spacing[1] = 1; - if ( ! (std::istringstream(argv[13]) >> spacing[2]) ) spacing[2] = 1; - - // Set region of interest in pixels - regionOfInterestMax.SetElement(0, size[0] - static_cast((radius)/spacing[0]) - 1); - regionOfInterestMax.SetElement(1, size[1] - static_cast((radius)/spacing[1]) - 1); - regionOfInterestMax.SetElement(2, size[2] - static_cast((radius)/spacing[2]) - 1); - regionOfInterestMin.SetElement(0, 0 + static_cast((radius+1)/spacing[0])+1); - regionOfInterestMin.SetElement(1, 0 + static_cast((radius+1)/spacing[1])+1); - //regionOfInterestMin.SetElement(2, 0 + static_cast((radius+1)/spacing[2])+1); - regionOfInterestMin.SetElement(0, 0 + static_cast((radius)/spacing[0])+1); - regionOfInterestMin.SetElement(1, 0 + static_cast((radius)/spacing[1])+1); - regionOfInterestMin.SetElement(2, 0 + static_cast((radius)/spacing[2])+1); - + MITK_TEST_BEGIN("mitkMultiGaussianTest") + + typedef double PixelType; + const unsigned int Dimension = 3; + typedef itk::Image ImageType; + typedef itk::MultiGaussianImageSource< ImageType > MultiGaussianImageSource; + itk::MultiGaussianImageSource< ImageType >::VectorType centerXVec, centerYVec, centerZVec, sigmaXVec, sigmaYVec, sigmaZVec, altitudeVec; + MultiGaussianImageSource::Pointer gaussianGenerator = MultiGaussianImageSource::New(); + ImageType::SizeValueType size[3]; + size[0] = 50; + size[1] = 50; + size[2] = 50; srand (time(NULL)); + unsigned int numberOfGaussians = 7; + unsigned int minWidthOfGaussian = (size[0] + size[1] + size[2]) / 27; // A ninth of the mean image size + unsigned int maxWidthOfGaussian = (size[0] + size[1] + size[2]) / 9; // One-third of the mean image size + unsigned int minAltitudeOfGaussian = 5; + unsigned int maxAltitudeOfGaussian = 200; + double centerX, centerY, centerZ, sigmaX, sigmaY, sigmaZ, altitude; + + gaussianGenerator->SetSize( size ); + gaussianGenerator->SetSpacing( 1 ); + gaussianGenerator->SetRadiusStepNumber(5); + gaussianGenerator->SetRadius(pow(itk::Math::one_over_pi * 0.75 , 1.0 / 3.0) * 10); + gaussianGenerator->SetNumberOfGausssians(numberOfGaussians); + // std::ofstream myfile; + // myfile.open ("C:/temp/tempParameter3.txt"); + // myfile << " CentX \t" << "Y \t" << "Z \t" << "SigX \t" << "Y \t" << "Z \t" << "Altit\n"; + int numberAddGaussian = numberOfGaussians; - for(unsigned int k = 1; k <= numberOfImages; ++k) + for( unsigned int i = 0; i < numberAddGaussian; ++i) { - gaussianGenerator = MultiGaussianImageSource::New(); - gaussianGenerator->SetSize( size ); - gaussianGenerator->SetSpacing( spacing ); - gaussianGenerator->SetRadiusStepNumber(8); - gaussianGenerator->SetRadius(radius); - gaussianGenerator->SetNumberOfGausssians(numberOfGaussians); - gaussianGenerator->SetRegionOfInterest(regionOfInterestMin, regionOfInterestMax); - // DOM Node Writer - xmlWriter = itk::DOMNodeXMLWriter::New(); - domTestCase = itk::DOMNode::New(); - domTestCase->SetName("testcase"); - domTestImage = itk::DOMNode::New(); - domTestImage->SetName("testimage"); - ss.str(""); - ss << size[0]; - domTestImage->SetAttribute("image-rows", ss.str()); - ss.str(""); - ss << size[1]; - domTestImage->SetAttribute("image-columns", ss.str()); - ss.str(""); - ss << size[2]; - domTestImage->SetAttribute("image-slices", ss.str()); - ss.str(""); - ss << numberOfGaussians; - domTestImage->SetAttribute("numberOfGaussians", ss.str()); - ss.str(""); - ss << spacing[0]; - domTestImage->SetAttribute("spacingX", ss.str()); - ss.str(""); - ss << spacing[1]; - domTestImage->SetAttribute("spacingY", ss.str()); - ss.str(""); - ss << spacing[2]; - domTestImage->SetAttribute("spacingZ", ss.str()); - domTestCase->AddChildAtBegin(domTestImage); - - for( unsigned int i = 0; i < numberAddGaussian; ++i) - { - - domGaussian = itk::DOMNode::New() ; - domGaussian->SetName("gaussian"); - domTestImage->AddChildAtEnd(domGaussian); - // generate the midpoint and the daviation in mm - centerX = rand() % static_cast(size[0] * spacing[0]); - ss.str(""); - ss << centerX; - domGaussian->SetAttribute("centerIndexX", ss.str()); - - centerY = rand() % static_cast(size[1] * spacing[1]); - ss.str(""); - ss << centerY; - domGaussian->SetAttribute("centerIndexY", ss.str()); - - centerZ = rand() % static_cast(size[2] * spacing[2]); - ss.str(""); - ss << centerZ; - domGaussian->SetAttribute("centerIndexZ", ss.str()); - - sigmaX = minWidthOfGaussian + rand() % (maxWidthOfGaussian - minWidthOfGaussian); - ss.str(""); - ss << sigmaX; - domGaussian->SetAttribute("deviationX", ss.str()); - - sigmaY = minWidthOfGaussian + rand() % (maxWidthOfGaussian - minWidthOfGaussian); - ss.str(""); - ss << sigmaY; - domGaussian->SetAttribute("deviationY", ss.str()); - - sigmaZ = minWidthOfGaussian + rand() % (maxWidthOfGaussian - minWidthOfGaussian); - ss.str(""); - ss << sigmaZ; - domGaussian->SetAttribute("deviationZ", ss.str()); - - altitude = minAltitudeOfGaussian + rand() % (maxAltitudeOfGaussian - minAltitudeOfGaussian); - ss.str(""); - ss << altitude; - domGaussian->SetAttribute("altitude", ss.str()); - - centerXVec.push_back(centerX); - centerYVec.push_back(centerY); - centerZVec.push_back(centerZ); - sigmaXVec.push_back(sigmaX); - sigmaYVec.push_back(sigmaY); - sigmaZVec.push_back(sigmaZ); - altitudeVec.push_back(altitude); - - } - - gaussianGenerator->AddGaussian(centerXVec, centerYVec, centerZVec, sigmaXVec, sigmaYVec, sigmaZVec, altitudeVec); - centerXVec.clear(); - centerYVec.clear(); - centerZVec.clear(); - sigmaXVec.clear(); - sigmaYVec.clear(); - sigmaZVec.clear(); - altitudeVec.clear(); - try { - gaussianGenerator->Update(); - gaussianGenerator->CalculateMidpointAndMeanValue(); - } catch (std::exception& e) - { - std::cout << "Error: " << e.what() << std::endl; - } - - //region of interest - domROI = itk::DOMNode::New(); - domROI->SetName("roi"); - domTestCase->AddChildAtEnd(domROI); - ss.str(""); - ss << regionOfInterestMax[0]; - domROI->SetAttribute("maximumX", ss.str()); - ss.str(""); - ss << regionOfInterestMin[0]; - domROI->SetAttribute("minimumX", ss.str()); - ss.str(""); - ss << regionOfInterestMax[1]; - domROI->SetAttribute("maximumY", ss.str()); - ss.str(""); - ss << regionOfInterestMin[1]; - domROI->SetAttribute("minimumY", ss.str()); - ss.str(""); - ss << regionOfInterestMax[2]; - domROI->SetAttribute("maximumZ", ss.str()); - ss.str(""); - ss << regionOfInterestMin[2]; - domROI->SetAttribute("minimumZ", ss.str()); - - //peak and peak coordinate - domStatistics = itk::DOMNode::New(); - domStatistics->SetName("statistic"); - domTestCase->AddChildAtEnd(domStatistics); - sphereMidpt = gaussianGenerator->GetSphereMidpoint(); - ss.str(""); - ss << sphereMidpt[0]; - domStatistics->SetAttribute("peakIndexX", ss.str()); - ss.str(""); - ss << sphereMidpt[1]; - domStatistics->SetAttribute("peakIndexY", ss.str()); - ss.str(""); - ss << sphereMidpt[2]; - domStatistics->SetAttribute("peakIndexZ", ss.str()); - ss.str(""); - ss << gaussianGenerator->GetMaxMeanValue(); - domStatistics->SetAttribute("peak", ss.str()); - - //optimize the mean value in the sphere - gaussianGenerator->OptimizeMeanValue(); - ss.str(""); - ss << gaussianGenerator->GetMaxMeanValue(); - domStatistics->SetAttribute("peakOptimized", ss.str()); - - - //maximum and maximum coordinate - gaussianGenerator->CalculateMaxAndMinInSphere(); - maxValueIndexInSphere = gaussianGenerator->GetMaxValueIndexInSphere(); - ss.str(""); - ss << maxValueIndexInSphere[0]; - domStatistics->SetAttribute("maximumIndexX", ss.str()); - ss.str(""); - ss << maxValueIndexInSphere[1]; - domStatistics->SetAttribute("maximumIndexY", ss.str()); - ss.str(""); - ss << maxValueIndexInSphere[2]; - domStatistics->SetAttribute("maximumIndexZ", ss.str()); - ss.str(""); - ss << gaussianGenerator->GetMaxValueInSphere(); - domStatistics->SetAttribute("maximum", ss.str()); - - //minimum and minimum coordinate - minValueIndexInSphere = gaussianGenerator->GetMinValueIndexInSphere(); - ss.str(""); - ss << minValueIndexInSphere[0]; - domStatistics->SetAttribute("minimumIndexX", ss.str()); - ss.str(""); - ss << minValueIndexInSphere[1]; - domStatistics->SetAttribute("minimumIndexY", ss.str()); - ss.str(""); - ss << minValueIndexInSphere[2]; - domStatistics->SetAttribute("minimumIndexZ", ss.str()); - ss.str(""); - ss << gaussianGenerator->GetMinValueInSphere(); - domStatistics->SetAttribute("minimum", ss.str()); - - // .xml (Data) - ss.str(""); - if(k < 10){ - ss << outputFilename <<"00"<< k <<".xml"; - }else if(k < 100){ - ss << outputFilename <<"0"<< k <<".xml"; - }else{ ss << outputFilename << k <<".xml";} - name = ss.str(); - fileNamePointer = (char*) name.c_str(); - xmlWriter->SetFileName( fileNamePointer); - xmlWriter->SetInput( domTestCase ); - xmlWriter->Update(); - ImageType::Pointer gaussianImage = gaussianGenerator->GetOutput(); - - //.nrrd (Image) - typedef itk::ImageFileWriter< ImageType > WriterType; - WriterType::Pointer writer = WriterType::New(); - ss.str(""); - if(k < 10) - { - ss << outputFilename <<"00"<< k <<".nrrd"; - }else if(k < 100) - { - ss << outputFilename <<"0"<< k <<".nrrd"; - }else - { - ss << outputFilename << k <<".nrrd"; - } - name = ss.str(); - fileNamePointer = (char*) name.c_str(); - writer->SetFileName( fileNamePointer); - writer->SetInput( gaussianImage ); - writer->Update(); + centerX = rand() % size[0]; + centerY = rand() % size[1]; + centerZ = rand() % size[2]; + sigmaX = minWidthOfGaussian + rand() % (maxWidthOfGaussian - minWidthOfGaussian); + sigmaY = minWidthOfGaussian + rand() % (maxWidthOfGaussian - minWidthOfGaussian); + sigmaZ = minWidthOfGaussian + rand() % (maxWidthOfGaussian - minWidthOfGaussian); + altitude = minAltitudeOfGaussian + rand() % (maxAltitudeOfGaussian - minAltitudeOfGaussian); + //gaussianGenerator->AddGaussian(centerX, centerY, centerZ, sigmaX, sigmaY, sigmaZ, altitude); + centerXVec.push_back(centerX); + centerYVec.push_back(centerY); + centerZVec.push_back(centerZ); + sigmaXVec.push_back(sigmaX); + sigmaYVec.push_back(sigmaY); + sigmaZVec.push_back(sigmaZ); + altitudeVec.push_back(altitude); + // myfile <AddGaussian(centerXVec, centerYVec, centerZVec, sigmaXVec, sigmaYVec, sigmaZVec, altitudeVec); + gaussianGenerator->Update(); + gaussianGenerator->CalculateMidpointAndMeanValue(); + std::cout << "Sphere radius is: " << gaussianGenerator->GetRadius() << std::endl; + std::cout << "Sphere midpoint is: " << gaussianGenerator->GetSphereMidpoint() << std::endl; + std::cout << "Mean value is: " << gaussianGenerator->GetMaxMeanValue() << std::endl; + ImageType::Pointer gaussianImage = gaussianGenerator->GetOutput(); + + //File writer + typedef itk::ImageFileWriter< ImageType > WriterType; + WriterType::Pointer writer = WriterType::New(); + writer->SetFileName( "C:/temp/tempImage33.nrrd" ); + writer->SetInput( gaussianImage ); + writer->Update(); + + MITK_TEST_END() } \ No newline at end of file diff --git a/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp b/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp index a7e206ebf6..20d0462f01 100644 --- a/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp +++ b/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp @@ -1,1672 +1,2003 @@ /*=================================================================== 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 #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 #include #include -#define DEBUG_HOTSPOTSEARCH +//#define DEBUG_HOTSPOTSEARCH #define _USE_MATH_DEFINES #include #if ( ( VTK_MAJOR_VERSION <= 5 ) && ( VTK_MINOR_VERSION<=8) ) -#include "mitkvtkLassoStencilSource.h" + #include "mitkvtkLassoStencilSource.h" #else -#include "vtkLassoStencilSource.h" + #include "vtkLassoStencilSource.h" #endif #include #include 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_EmptyHistogram = HistogramType::New(); - m_EmptyHistogram->SetMeasurementVectorSize(1); - HistogramType::SizeType histogramSize(1); - histogramSize.Fill( 256 ); - m_EmptyHistogram->Initialize( histogramSize ); +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_EmptyHistogram = HistogramType::New(); + m_EmptyHistogram->SetMeasurementVectorSize(1); + HistogramType::SizeType histogramSize(1); + histogramSize.Fill( 256 ); + m_EmptyHistogram->Initialize( histogramSize ); - m_EmptyStatistics.Reset(); - } + m_EmptyStatistics.Reset(); +} - ImageStatisticsCalculator::~ImageStatisticsCalculator() - { - } +ImageStatisticsCalculator::~ImageStatisticsCalculator() +{ +} - void ImageStatisticsCalculator::SetImage( const mitk::Image *image ) +void ImageStatisticsCalculator::SetImage( const mitk::Image *image ) +{ + if ( m_Image != image ) { - if ( m_Image != image ) - { - m_Image = image; - this->Modified(); + m_Image = image; + this->Modified(); - unsigned int numberOfTimeSteps = image->GetTimeSteps(); + 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 ); + // 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_ImageStatisticsVector.resize( numberOfTimeSteps ); + m_MaskedImageStatisticsVector.resize( numberOfTimeSteps ); + m_PlanarFigureStatisticsVector.resize( numberOfTimeSteps ); - m_ImageStatisticsTimeStampVector.resize( numberOfTimeSteps ); - m_MaskedImageStatisticsTimeStampVector.resize( numberOfTimeSteps ); - m_PlanarFigureStatisticsTimeStampVector.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 ); + 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; - } + for ( unsigned int t = 0; t < image->GetTimeSteps(); ++t ) + { + m_ImageStatisticsTimeStampVector[t].Modified(); + m_ImageStatisticsCalculationTriggerVector[t] = true; } } +} - void ImageStatisticsCalculator::SetImageMask( const mitk::Image *imageMask ) +void ImageStatisticsCalculator::SetImageMask( const mitk::Image *imageMask ) +{ + if ( m_Image.IsNull() ) { - if ( m_Image.IsNull() ) - { - itkExceptionMacro( << "Image needs to be set first!" ); - } - - if ( m_Image->GetTimeSteps() != imageMask->GetTimeSteps() ) - { - itkExceptionMacro( << "Image and image mask need to have equal number of time steps!" ); - } - - 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; - } - } + itkExceptionMacro( << "Image needs to be set first!" ); } + if ( m_Image->GetTimeSteps() != imageMask->GetTimeSteps() ) + { + itkExceptionMacro( << "Image and image mask need to have equal number of time steps!" ); + } - void ImageStatisticsCalculator::SetPlanarFigure( mitk::PlanarFigure *planarFigure ) + if ( m_ImageMask != imageMask ) { - if ( m_Image.IsNull() ) - { - itkExceptionMacro( << "Image needs to be set first!" ); - } + m_ImageMask = imageMask; + this->Modified(); - if ( m_PlanarFigure != planarFigure ) + for ( unsigned int t = 0; t < m_Image->GetTimeSteps(); ++t ) { - m_PlanarFigure = planarFigure; - this->Modified(); - - for ( unsigned int t = 0; t < m_Image->GetTimeSteps(); ++t ) - { - m_PlanarFigureStatisticsTimeStampVector[t].Modified(); - m_PlanarFigureStatisticsCalculationTriggerVector[t] = true; - } + m_MaskedImageStatisticsTimeStampVector[t].Modified(); + m_MaskedImageStatisticsCalculationTriggerVector[t] = true; } } +} - void ImageStatisticsCalculator::SetMaskingMode( unsigned int mode ) +void ImageStatisticsCalculator::SetPlanarFigure( mitk::PlanarFigure *planarFigure ) +{ + if ( m_Image.IsNull() ) { - if ( m_MaskingMode != mode ) - { - m_MaskingMode = mode; - m_MaskingModeChanged = true; - this->Modified(); - } + itkExceptionMacro( << "Image needs to be set first!" ); } - - void ImageStatisticsCalculator::SetMaskingModeToNone() + if ( m_PlanarFigure != planarFigure ) { - if ( m_MaskingMode != MASKING_MODE_NONE ) + m_PlanarFigure = planarFigure; + this->Modified(); + + for ( unsigned int t = 0; t < m_Image->GetTimeSteps(); ++t ) { - m_MaskingMode = MASKING_MODE_NONE; - m_MaskingModeChanged = true; - this->Modified(); + m_PlanarFigureStatisticsTimeStampVector[t].Modified(); + m_PlanarFigureStatisticsCalculationTriggerVector[t] = true; } } +} - void ImageStatisticsCalculator::SetMaskingModeToImage() +void ImageStatisticsCalculator::SetMaskingMode( unsigned int mode ) +{ + if ( m_MaskingMode != mode ) { - if ( m_MaskingMode != MASKING_MODE_IMAGE ) - { - m_MaskingMode = MASKING_MODE_IMAGE; - m_MaskingModeChanged = true; - this->Modified(); - } + m_MaskingMode = mode; + m_MaskingModeChanged = true; + this->Modified(); } +} - void ImageStatisticsCalculator::SetMaskingModeToPlanarFigure() +void ImageStatisticsCalculator::SetMaskingModeToNone() +{ + if ( m_MaskingMode != MASKING_MODE_NONE ) { - if ( m_MaskingMode != MASKING_MODE_PLANARFIGURE ) - { - m_MaskingMode = MASKING_MODE_PLANARFIGURE; - m_MaskingModeChanged = true; - this->Modified(); - } + m_MaskingMode = MASKING_MODE_NONE; + m_MaskingModeChanged = true; + this->Modified(); } +} - void ImageStatisticsCalculator::SetIgnorePixelValue(double value) + +void ImageStatisticsCalculator::SetMaskingModeToImage() +{ + if ( m_MaskingMode != MASKING_MODE_IMAGE ) { - if ( m_IgnorePixelValue != value ) - { - m_IgnorePixelValue = value; - if(m_DoIgnorePixelValue) - { - m_IgnorePixelValueChanged = true; - } - this->Modified(); - } + m_MaskingMode = MASKING_MODE_IMAGE; + m_MaskingModeChanged = true; + this->Modified(); } +} + - double ImageStatisticsCalculator::GetIgnorePixelValue() +void ImageStatisticsCalculator::SetMaskingModeToPlanarFigure() +{ + if ( m_MaskingMode != MASKING_MODE_PLANARFIGURE ) { - return m_IgnorePixelValue; + m_MaskingMode = MASKING_MODE_PLANARFIGURE; + m_MaskingModeChanged = true; + this->Modified(); } +} - void ImageStatisticsCalculator::SetDoIgnorePixelValue(bool value) +void ImageStatisticsCalculator::SetIgnorePixelValue(double value) +{ + if ( m_IgnorePixelValue != value ) { - if ( m_DoIgnorePixelValue != value ) + m_IgnorePixelValue = value; + if(m_DoIgnorePixelValue) { - m_DoIgnorePixelValue = value; m_IgnorePixelValueChanged = true; - this->Modified(); } + this->Modified(); } +} - bool ImageStatisticsCalculator::GetDoIgnorePixelValue() +double ImageStatisticsCalculator::GetIgnorePixelValue() +{ + return m_IgnorePixelValue; +} + +void ImageStatisticsCalculator::SetDoIgnorePixelValue(bool value) +{ + if ( m_DoIgnorePixelValue != value ) { - return m_DoIgnorePixelValue; + m_DoIgnorePixelValue = value; + m_IgnorePixelValueChanged = true; + this->Modified(); } +} - void ImageStatisticsCalculator::SetHotspotRadius(double value) +bool ImageStatisticsCalculator::GetDoIgnorePixelValue() +{ + return m_DoIgnorePixelValue; +} + +void ImageStatisticsCalculator::SetHotspotRadius(double value) +{ + m_HotspotRadius = value; +} + +double ImageStatisticsCalculator::GetHotspotRadius() +{ + return m_HotspotRadius; +} + +void ImageStatisticsCalculator::SetCalculateHotspot(bool value) +{ + m_CalculateHotspot = value; +} + +bool ImageStatisticsCalculator::IsHotspotCalculated() +{ + return m_CalculateHotspot; +} + +bool ImageStatisticsCalculator::ComputeStatistics( unsigned int timeStep ) +{ + + if (m_Image.IsNull() ) { - m_HotspotRadius = value; + mitkThrow() << "Image not set!"; } - double ImageStatisticsCalculator::GetHotspotRadius() + if (!m_Image->IsInitialized()) { - return m_HotspotRadius; + mitkThrow() << "Image not initialized!"; } - void ImageStatisticsCalculator::SetCalculateHotspot(bool value) + if ( m_Image->GetReferenceCount() == 1 ) { - m_CalculateHotspot = value; + // Image no longer valid; we are the only ones to still hold a reference on it + return false; } - bool ImageStatisticsCalculator::IsHotspotCalculated() + if ( timeStep >= m_Image->GetTimeSteps() ) { - return m_CalculateHotspot; + throw std::runtime_error( "Error: invalid time step!" ); } - bool ImageStatisticsCalculator::ComputeStatistics( unsigned int timeStep ) + // 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 = NULL; + } - if (m_Image.IsNull() ) - { - mitkThrow() << "Image not set!"; - } + // 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(); - if (!m_Image->IsInitialized()) - { - mitkThrow() << "Image not initialized!"; - } + bool imageStatisticsCalculationTrigger = m_ImageStatisticsCalculationTriggerVector[timeStep]; + bool maskedImageStatisticsCalculationTrigger = m_MaskedImageStatisticsCalculationTriggerVector[timeStep]; + bool planarFigureStatisticsCalculationTrigger = m_PlanarFigureStatisticsCalculationTriggerVector[timeStep]; - if ( m_Image->GetReferenceCount() == 1 ) + if ( !m_IgnorePixelValueChanged + && ((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 ) { - Image no longer valid; we are the only ones to still hold a reference on it - return false; + m_MaskingModeChanged = false; + return true; } - - if ( timeStep >= m_Image->GetTimeSteps() ) + else { - throw std::runtime_error( "Error: invalid time step!" ); + return false; } + } - 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 = NULL; - } + // Reset state changed flag + m_MaskingModeChanged = false; + m_IgnorePixelValueChanged = false; - 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(); + // Depending on masking mode, extract and/or generate the required image + // and mask data from the user input + this->ExtractImageAndMask( timeStep ); - bool imageStatisticsCalculationTrigger = m_ImageStatisticsCalculationTriggerVector[timeStep]; - bool maskedImageStatisticsCalculationTrigger = m_MaskedImageStatisticsCalculationTriggerVector[timeStep]; - bool planarFigureStatisticsCalculationTrigger = m_PlanarFigureStatisticsCalculationTriggerVector[timeStep]; - if ( !m_IgnorePixelValueChanged - && ((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)) ) + StatisticsContainer *statisticsContainer; + HistogramContainer *histogramContainer; + switch ( m_MaskingMode ) + { + case MASKING_MODE_NONE: + default: + if(!m_DoIgnorePixelValue) { - Statistics is up to date! - if ( m_MaskingModeChanged ) - { - m_MaskingModeChanged = false; - return true; - } - 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 = &m_ImageStatisticsVector[timeStep]; + histogramContainer = &m_ImageHistogramVector[timeStep]; - StatisticsContainer *statisticsContainer; - HistogramContainer *histogramContainer; - switch ( m_MaskingMode ) + m_ImageStatisticsTimeStampVector[timeStep].Modified(); + m_ImageStatisticsCalculationTriggerVector[timeStep] = false; + } + else { - 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; + } + break; - case MASKING_MODE_PLANARFIGURE: - statisticsContainer = &m_PlanarFigureStatisticsVector[timeStep]; - histogramContainer = &m_PlanarFigureHistogramVector[timeStep]; + case MASKING_MODE_IMAGE: + statisticsContainer = &m_MaskedImageStatisticsVector[timeStep]; + histogramContainer = &m_MaskedImageHistogramVector[timeStep]; - m_PlanarFigureStatisticsTimeStampVector[timeStep].Modified(); - m_PlanarFigureStatisticsCalculationTriggerVector[timeStep] = false; - break; - } + m_MaskedImageStatisticsTimeStampVector[timeStep].Modified(); + m_MaskedImageStatisticsCalculationTriggerVector[timeStep] = false; + break; - Calculate statistics and histogram(s) - if ( m_InternalImage->GetDimension() == 3 ) + 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 ) { - 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 ); - } + AccessFixedDimensionByItk_2( + m_InternalImage, + InternalCalculateStatisticsUnmasked, + 3, + statisticsContainer, + histogramContainer ); } - else if ( m_InternalImage->GetDimension() == 2 ) + else { - 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 ); - } + 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 { - MITK_ERROR << "ImageStatistics: Image dimension not supported!"; + 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; - } + // Release unused image smart pointers to free memory + m_InternalImage = mitk::Image::ConstPointer(); + m_InternalImageMask3D = MaskImage3DType::Pointer(); + m_InternalImageMask2D = MaskImage2DType::Pointer(); + return true; +} - const ImageStatisticsCalculator::HistogramType * - ImageStatisticsCalculator::GetHistogram( unsigned int timeStep, unsigned int label ) const +const ImageStatisticsCalculator::HistogramType * +ImageStatisticsCalculator::GetHistogram( unsigned int timeStep, unsigned int label ) const +{ + if ( m_Image.IsNull() || (timeStep >= m_Image->GetTimeSteps()) ) { - if ( m_Image.IsNull() || (timeStep >= m_Image->GetTimeSteps()) ) - { - return NULL; - } + return NULL; + } - switch ( m_MaskingMode ) + switch ( m_MaskingMode ) + { + case MASKING_MODE_NONE: + default: { - case MASKING_MODE_NONE: - default: - { - if(m_DoIgnorePixelValue) - return m_MaskedImageHistogramVector[timeStep][label]; + if(m_DoIgnorePixelValue) + return m_MaskedImageHistogramVector[timeStep][label]; - return m_ImageHistogramVector[timeStep][label]; - } + return m_ImageHistogramVector[timeStep][label]; + } - case MASKING_MODE_IMAGE: - return m_MaskedImageHistogramVector[timeStep][label]; + case MASKING_MODE_IMAGE: + return m_MaskedImageHistogramVector[timeStep][label]; - case MASKING_MODE_PLANARFIGURE: - return m_PlanarFigureHistogramVector[timeStep][label]; - } + case MASKING_MODE_PLANARFIGURE: + return m_PlanarFigureHistogramVector[timeStep][label]; } +} - const ImageStatisticsCalculator::HistogramContainer & - ImageStatisticsCalculator::GetHistogramVector( unsigned int timeStep ) const +const ImageStatisticsCalculator::HistogramContainer & +ImageStatisticsCalculator::GetHistogramVector( unsigned int timeStep ) const +{ + if ( m_Image.IsNull() || (timeStep >= m_Image->GetTimeSteps()) ) { - if ( m_Image.IsNull() || (timeStep >= m_Image->GetTimeSteps()) ) - { - return m_EmptyHistogramContainer; - } + return m_EmptyHistogramContainer; + } - switch ( m_MaskingMode ) + switch ( m_MaskingMode ) + { + case MASKING_MODE_NONE: + default: { - case MASKING_MODE_NONE: - default: - { - if(m_DoIgnorePixelValue) - return m_MaskedImageHistogramVector[timeStep]; + if(m_DoIgnorePixelValue) + return m_MaskedImageHistogramVector[timeStep]; - return m_ImageHistogramVector[timeStep]; - } + return m_ImageHistogramVector[timeStep]; + } - case MASKING_MODE_IMAGE: - return m_MaskedImageHistogramVector[timeStep]; + case MASKING_MODE_IMAGE: + return m_MaskedImageHistogramVector[timeStep]; - case MASKING_MODE_PLANARFIGURE: - return m_PlanarFigureHistogramVector[timeStep]; - } + case MASKING_MODE_PLANARFIGURE: + return m_PlanarFigureHistogramVector[timeStep]; } +} - const ImageStatisticsCalculator::Statistics & - ImageStatisticsCalculator::GetStatistics( unsigned int timeStep, unsigned int label ) const +const ImageStatisticsCalculator::Statistics & +ImageStatisticsCalculator::GetStatistics( unsigned int timeStep, unsigned int label ) const +{ + if ( m_Image.IsNull() || (timeStep >= m_Image->GetTimeSteps()) ) { - if ( m_Image.IsNull() || (timeStep >= m_Image->GetTimeSteps()) ) - { - return m_EmptyStatistics; - } + return m_EmptyStatistics; + } - switch ( m_MaskingMode ) + switch ( m_MaskingMode ) + { + case MASKING_MODE_NONE: + default: { - 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]; + if(m_DoIgnorePixelValue) + return m_MaskedImageStatisticsVector[timeStep][label]; - case MASKING_MODE_PLANARFIGURE: - return m_PlanarFigureStatisticsVector[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 +const ImageStatisticsCalculator::StatisticsContainer & +ImageStatisticsCalculator::GetStatisticsVector( unsigned int timeStep ) const +{ + if ( m_Image.IsNull() || (timeStep >= m_Image->GetTimeSteps()) ) { - if ( m_Image.IsNull() || (timeStep >= m_Image->GetTimeSteps()) ) - { - return m_EmptyStatisticsContainer; - } + return m_EmptyStatisticsContainer; + } - switch ( m_MaskingMode ) + switch ( m_MaskingMode ) + { + case MASKING_MODE_NONE: + default: { - case MASKING_MODE_NONE: - default: - { - if(m_DoIgnorePixelValue) - return m_MaskedImageStatisticsVector[timeStep]; - - return m_ImageStatisticsVector[timeStep]; - } - case MASKING_MODE_IMAGE: - return m_MaskedImageStatisticsVector[timeStep]; + if(m_DoIgnorePixelValue) + return m_MaskedImageStatisticsVector[timeStep]; - case MASKING_MODE_PLANARFIGURE: - return m_PlanarFigureStatisticsVector[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 ) +void ImageStatisticsCalculator::ExtractImageAndMask( unsigned int timeStep ) +{ + if ( m_Image.IsNull() ) { - if ( m_Image.IsNull() ) - { - throw std::runtime_error( "Error: image empty!" ); - } + throw std::runtime_error( "Error: image empty!" ); + } - if ( timeStep >= m_Image->GetTimeSteps() ) - { - throw std::runtime_error( "Error: invalid time step!" ); - } + 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(); + ImageTimeSelector::Pointer imageTimeSelector = ImageTimeSelector::New(); + imageTimeSelector->SetInput( m_Image ); + imageTimeSelector->SetTimeNr( timeStep ); + imageTimeSelector->UpdateLargestPossibleRegion(); + mitk::Image *timeSliceImage = imageTimeSelector->GetOutput(); - switch ( m_MaskingMode ) + switch ( m_MaskingMode ) + { + case MASKING_MODE_NONE: { - case MASKING_MODE_NONE: - { - m_InternalImage = timeSliceImage; - m_InternalImageMask2D = NULL; - m_InternalImageMask3D = NULL; + m_InternalImage = timeSliceImage; + m_InternalImageMask2D = NULL; + m_InternalImageMask3D = NULL; - if(m_DoIgnorePixelValue) + if(m_DoIgnorePixelValue) + { + if( m_InternalImage->GetDimension() == 3 ) { - if( m_InternalImage->GetDimension() == 3 ) - { - CastToItkImage( timeSliceImage, m_InternalImageMask3D ); - m_InternalImageMask3D->FillBuffer(1); - } - if( m_InternalImage->GetDimension() == 2 ) - { - CastToItkImage( timeSliceImage, m_InternalImageMask2D ); - m_InternalImageMask2D->FillBuffer(1); - } + CastToItkImage( timeSliceImage, m_InternalImageMask3D ); + m_InternalImageMask3D->FillBuffer(1); + } + if( m_InternalImage->GetDimension() == 2 ) + { + CastToItkImage( timeSliceImage, m_InternalImageMask2D ); + m_InternalImageMask2D->FillBuffer(1); } - break; } + break; + } - case MASKING_MODE_IMAGE: + case MASKING_MODE_IMAGE: + { + if ( m_ImageMask.IsNotNull() && (m_ImageMask->GetReferenceCount() > 1) ) { - if ( m_ImageMask.IsNotNull() && (m_ImageMask->GetReferenceCount() > 1) ) + if ( timeStep < m_ImageMask->GetTimeSteps() ) { - if ( timeStep < m_ImageMask->GetTimeSteps() ) - { - 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 has not enough time steps!" ); - } + 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!" ); + throw std::runtime_error( "Error: image mask has not enough time steps!" ); } - break; } - - case MASKING_MODE_PLANARFIGURE: + else { - m_InternalImageMask2D = NULL; + throw std::runtime_error( "Error: image mask empty!" ); + } + break; + } - 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" ); - } + case MASKING_MODE_PLANARFIGURE: + { + m_InternalImageMask2D = NULL; - const Geometry3D *imageGeometry = timeSliceImage->GetGeometry(); - if ( imageGeometry == NULL ) - { - throw std::runtime_error( "Image geometry invalid!" ); - } + 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 Geometry2D *planarFigureGeometry2D = m_PlanarFigure->GetGeometry2D(); - if ( planarFigureGeometry2D == NULL ) - { - throw std::runtime_error( "Planar-Figure not yet initialized!" ); - } + const Geometry3D *imageGeometry = timeSliceImage->GetGeometry(); + if ( imageGeometry == NULL ) + { + throw std::runtime_error( "Image geometry invalid!" ); + } - const PlaneGeometry *planarFigureGeometry = - dynamic_cast< const PlaneGeometry * >( planarFigureGeometry2D ); - if ( planarFigureGeometry == NULL ) - { - throw std::runtime_error( "Non-planar planar figures not supported!" ); - } + const Geometry2D *planarFigureGeometry2D = m_PlanarFigure->GetGeometry2D(); + if ( planarFigureGeometry2D == NULL ) + { + throw std::runtime_error( "Planar-Figure not yet initialized!" ); + } - 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; + const PlaneGeometry *planarFigureGeometry = + dynamic_cast< const PlaneGeometry * >( planarFigureGeometry2D ); + if ( planarFigureGeometry == NULL ) + { + 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; + // 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; - } + // Extract slice with given position and direction from image + unsigned int dimension = timeSliceImage->GetDimension(); - Compute mask from PlanarFigure - AccessFixedDimensionByItk_1( - m_InternalImage, - InternalCalculateMaskFromPlanarFigure, - 2, axis ); + 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_DoIgnorePixelValue) + { + if ( m_InternalImage->GetDimension() == 3 ) { - if ( m_InternalImage->GetDimension() == 3 ) - { - AccessFixedDimensionByItk_1( + AccessFixedDimensionByItk_1( m_InternalImage, InternalMaskIgnoredPixels, 3, m_InternalImageMask3D.GetPointer() ); - } - else if ( m_InternalImage->GetDimension() == 2 ) - { - AccessFixedDimensionByItk_1( + } + else if ( m_InternalImage->GetDimension() == 2 ) + { + AccessFixedDimensionByItk_1( m_InternalImage, InternalMaskIgnoredPixels, 2, m_InternalImageMask2D.GetPointer() ); - } } } +} - bool ImageStatisticsCalculator::GetPrincipalAxis( - const Geometry3D *geometry, Vector3D vector, - unsigned int &axis ) +bool ImageStatisticsCalculator::GetPrincipalAxis( + const Geometry3D *geometry, Vector3D vector, + unsigned int &axis ) +{ + vector.Normalize(); + for ( unsigned int i = 0; i < 3; ++i ) { - vector.Normalize(); - for ( unsigned int i = 0; i < 3; ++i ) - { - Vector3D axisVector = geometry->GetAxisVector( i ); - axisVector.Normalize(); + Vector3D axisVector = geometry->GetAxisVector( i ); + axisVector.Normalize(); - if ( fabs( fabs( axisVector * vector ) - 1.0) < mitk::eps ) - { - axis = i; - return true; - } + if ( fabs( fabs( axisVector * vector ) - 1.0) < mitk::eps ) + { + axis = i; + return true; } - - return false; } + return false; +} - 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 itk::Image< unsigned short, VImageDimension > MaskImageType; - typedef typename ImageType::IndexType IndexType; - typedef itk::Statistics::ScalarImageToHistogramGenerator< ImageType > - HistogramGeneratorType; +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 itk::Image< unsigned short, VImageDimension > MaskImageType; + typedef typename ImageType::IndexType IndexType; - statisticsContainer->clear(); - histogramContainer->clear(); + typedef itk::Statistics::ScalarImageToHistogramGenerator< ImageType > + HistogramGeneratorType; - Progress listening... - typedef itk::SimpleMemberCommand< ImageStatisticsCalculator > ITKCommandType; - ITKCommandType::Pointer progressListener; - progressListener = ITKCommandType::New(); - progressListener->SetCallbackFunction( this, - &ImageStatisticsCalculator::UnmaskedStatisticsProgressUpdate ); + 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::StatisticsImageFilter< ImageType > StatisticsFilterType; - typename StatisticsFilterType::Pointer statisticsFilter = StatisticsFilterType::New(); - statisticsFilter->SetInput( image ); - unsigned long observerTag = statisticsFilter->AddObserver( itk::ProgressEvent(), progressListener ); - statisticsFilter->Update(); - 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.Label = 1; - statistics.N = image->GetBufferedRegion().GetNumberOfPixels(); - statistics.Min = statisticsFilter->GetMinimum(); - statistics.Max = statisticsFilter->GetMaximum(); - statistics.Mean = statisticsFilter->GetMean(); - statistics.Median = 0.0; - statistics.Sigma = statisticsFilter->GetSigma(); - statistics.RMS = sqrt( statistics.Mean * statistics.Mean + statistics.Sigma * statistics.Sigma ); - - statistics.MinIndex.set_size(image->GetImageDimension()); - statistics.MaxIndex.set_size(image->GetImageDimension()); - for (int i=0; iInvokeEvent( itk::StartEvent() ); + for ( unsigned int i = 0; i < 100; ++i ) + { + this->UnmaskedStatisticsProgressUpdate(); + } + + // Calculate statistics (separate filter) + typedef itk::StatisticsImageFilter< ImageType > StatisticsFilterType; + typename StatisticsFilterType::Pointer statisticsFilter = StatisticsFilterType::New(); + statisticsFilter->SetInput( image ); + unsigned long observerTag = statisticsFilter->AddObserver( itk::ProgressEvent(), progressListener ); + statisticsFilter->Update(); + 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.Label = 1; + statistics.N = image->GetBufferedRegion().GetNumberOfPixels(); + statistics.Min = statisticsFilter->GetMinimum(); + statistics.Max = statisticsFilter->GetMaximum(); + statistics.Mean = statisticsFilter->GetMean(); + statistics.Median = 0.0; + statistics.Sigma = statisticsFilter->GetSigma(); + statistics.RMS = sqrt( statistics.Mean * statistics.Mean + statistics.Sigma * statistics.Sigma ); + + statistics.MinIndex.set_size(image->GetImageDimension()); + statistics.MaxIndex.set_size(image->GetImageDimension()); + for (int i=0; iGetIndexOfMaximum()[i]; statistics.MinIndex[i] = minMaxFilter->GetIndexOfMinimum()[i]; - } + } - statisticsContainer->push_back( statistics ); + statisticsContainer->push_back( statistics ); - Calculate histogram - typename HistogramGeneratorType::Pointer histogramGenerator = HistogramGeneratorType::New(); - histogramGenerator->SetInput( image ); - histogramGenerator->SetMarginalScale( 100 ); - histogramGenerator->SetNumberOfBins( 768 ); - histogramGenerator->SetHistogramMin( statistics.Min ); - histogramGenerator->SetHistogramMax( statistics.Max ); - histogramGenerator->Compute(); + // Calculate histogram + typename HistogramGeneratorType::Pointer histogramGenerator = HistogramGeneratorType::New(); + histogramGenerator->SetInput( image ); + histogramGenerator->SetMarginalScale( 100 ); + histogramGenerator->SetNumberOfBins( 768 ); + histogramGenerator->SetHistogramMin( statistics.Min ); + histogramGenerator->SetHistogramMax( statistics.Max ); + histogramGenerator->Compute(); - histogramContainer->push_back( histogramGenerator->GetOutput() ); - } + 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; +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 + itk::ImageRegionIterator itmask(maskImage, maskImage->GetLargestPossibleRegion()); - itk::ImageRegionConstIterator + itk::ImageRegionConstIterator itimage(image, image->GetLargestPossibleRegion()); - itmask.GoToBegin(); - itimage.GoToBegin(); + itmask.GoToBegin(); + itimage.GoToBegin(); - while( !itmask.IsAtEnd() ) + while( !itmask.IsAtEnd() ) + { + if(m_IgnorePixelValue == itimage.Get()) { - if(m_IgnorePixelValue == itimage.Get()) - { - itmask.Set(0); - } - - ++itmask; - ++itimage; + 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; +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::IndexType IndexType; + typedef typename ImageType::PointType PointType; + typedef typename ImageType::SpacingType SpacingType; - typedef itk::LabelStatisticsImageFilter< ImageType, MaskImageType > LabelStatisticsFilterType; + typedef itk::LabelStatisticsImageFilter< ImageType, MaskImageType > LabelStatisticsFilterType; - typedef itk::ChangeInformationImageFilter< MaskImageType > ChangeInformationFilterType; + typedef itk::ChangeInformationImageFilter< MaskImageType > ChangeInformationFilterType; - typedef itk::ExtractImageFilter< ImageType, ImageType > ExtractImageFilterType; + typedef itk::ExtractImageFilter< ImageType, ImageType > ExtractImageFilterType; - statisticsContainer->clear(); - histogramContainer->clear(); + statisticsContainer->clear(); + histogramContainer->clear(); - Make sure that mask is set - if ( maskImage == NULL ) - { - itkExceptionMacro( << "Mask image needs to be set!" ); - } + // Make sure that mask is set + if ( maskImage == NULL ) + { + 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 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 ) + // 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 ) { - 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 ) > mitk::eps ) - { - itkExceptionMacro( << "Mask needs to have same direction as image! (Image direction: " << imageDirection << "; Mask direction: " << maskDirection << ")" ); - } + 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]; + // 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; + typedef itk::ContinuousIndex ContinousIndexType; + ContinousIndexType maskOriginContinousIndex, imageOriginContinousIndex; - image->TransformPhysicalPointToContinuousIndex(maskOrigin, maskOriginContinousIndex); - image->TransformPhysicalPointToContinuousIndex(imageOrigin, imageOriginContinousIndex); + image->TransformPhysicalPointToContinuousIndex(maskOrigin, maskOriginContinousIndex); + image->TransformPhysicalPointToContinuousIndex(imageOrigin, imageOriginContinousIndex); - for ( unsigned int i = 0; i < ImageType::ImageDimension; ++i ) + for ( unsigned int i = 0; i < ImageType::ImageDimension; ++i ) + { + double misalignment = maskOriginContinousIndex[i] - floor( maskOriginContinousIndex[i] + 0.5 ); + if ( fabs( misalignment ) > mitk::eps ) { - double misalignment = maskOriginContinousIndex[i] - floor( maskOriginContinousIndex[i] + 0.5 ); - if ( fabs( misalignment ) > mitk::eps ) - { - itkExceptionMacro( << "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]; + itkExceptionMacro( << "Pixels/voxels of mask and image are not sufficiently aligned! (Misalignment: " << misalignment << ")" ); } - 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 ); - adaptMaskFilter->Update(); - typename MaskImageType::Pointer adaptedMaskImage = adaptMaskFilter->GetOutput(); + double indexCoordDistance = maskOriginContinousIndex[i] - imageOriginContinousIndex[i]; + offset[i] = (int) indexCoordDistance + image->GetBufferedRegion().GetIndex()[i]; + } + // 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 ); + adaptMaskFilter->Update(); + typename MaskImageType::Pointer adaptedMaskImage = adaptMaskFilter->GetOutput(); - Make sure that mask region is contained within image region - if ( !image->GetLargestPossibleRegion().IsInside( adaptedMaskImage->GetLargestPossibleRegion() ) ) - { - itkExceptionMacro( << "Mask region needs to be inside of image region! (Image region: " - << image->GetLargestPossibleRegion() << "; Mask region: " << adaptedMaskImage->GetLargestPossibleRegion() << ")" ); - } + // Make sure that mask region is contained within image region + if ( !image->GetLargestPossibleRegion().IsInside( adaptedMaskImage->GetLargestPossibleRegion() ) ) + { + itkExceptionMacro( << "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 + // 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] ) { - adaptedImage = image; + 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 ); - - statisticsFilter->Update(); - - int numberOfBins = ( m_DoIgnorePixelValue && (m_MaskingMode == MASKING_MODE_NONE) ) ? 768 : 384; - typename LabelStatisticsFilterType::Pointer labelStatisticsFilter; - labelStatisticsFilter = LabelStatisticsFilterType::New(); - labelStatisticsFilter->SetInput( adaptedImage ); - labelStatisticsFilter->SetLabelInput( adaptedMaskImage ); - labelStatisticsFilter->UseHistogramsOn(); - labelStatisticsFilter->SetHistogramParameters( numberOfBins, 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 - labelStatisticsFilter->Update(); - this->InvokeEvent( itk::EndEvent() ); - labelStatisticsFilter->RemoveObserver( observerTag ); - - Find all relevant labels of mask (other than 0) - std::list< int > relevantLabels; - bool maskNonEmpty = false; - unsigned int i; - for ( i = 1; i < 4096; ++i ) + // Initialize Filter + typedef itk::StatisticsImageFilter< ImageType > StatisticsFilterType; + typename StatisticsFilterType::Pointer statisticsFilter = StatisticsFilterType::New(); + statisticsFilter->SetInput( adaptedImage ); + + statisticsFilter->Update(); + + int numberOfBins = ( m_DoIgnorePixelValue && (m_MaskingMode == MASKING_MODE_NONE) ) ? 768 : 384; + typename LabelStatisticsFilterType::Pointer labelStatisticsFilter; + labelStatisticsFilter = LabelStatisticsFilterType::New(); + labelStatisticsFilter->SetInput( adaptedImage ); + labelStatisticsFilter->SetLabelInput( adaptedMaskImage ); + labelStatisticsFilter->UseHistogramsOn(); + labelStatisticsFilter->SetHistogramParameters( numberOfBins, 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 + labelStatisticsFilter->Update(); + this->InvokeEvent( itk::EndEvent() ); + labelStatisticsFilter->RemoveObserver( observerTag ); + + // Find all relevant labels of mask (other than 0) + std::list< int > relevantLabels; + bool maskNonEmpty = false; + unsigned int i; + for ( i = 1; i < 4096; ++i ) + { + if ( labelStatisticsFilter->HasLabel( i ) ) { - if ( labelStatisticsFilter->HasLabel( i ) ) - { - relevantLabels.push_back( i ); - maskNonEmpty = true; - } + relevantLabels.push_back( i ); + maskNonEmpty = true; } + } - if ( maskNonEmpty ) + if ( maskNonEmpty ) + { + std::list< int >::iterator it; + for ( it = relevantLabels.begin(), i = 0; + it != relevantLabels.end(); + ++it, ++i ) { + histogramContainer->push_back( HistogramType::ConstPointer( labelStatisticsFilter->GetHistogram( (*it) ) ) ); + Statistics statistics; - std::list< int >::iterator it; - for ( it = relevantLabels.begin(), i = 0; - it != relevantLabels.end(); - ++it, ++i ) - { - histogramContainer->push_back( HistogramType::ConstPointer( labelStatisticsFilter->GetHistogram( (*it) ) ) ); - - statistics.Label = (*it); - statistics.N = labelStatisticsFilter->GetCount( *it ); - statistics.Min = labelStatisticsFilter->GetMinimum( *it ); - statistics.Max = labelStatisticsFilter->GetMaximum( *it ); - statistics.Mean = labelStatisticsFilter->GetMean( *it ); - statistics.Median = labelStatisticsFilter->GetMedian( *it ); - statistics.Sigma = labelStatisticsFilter->GetSigma( *it ); - statistics.RMS = sqrt( statistics.Mean * statistics.Mean - + statistics.Sigma * statistics.Sigma ); - - restrict image to mask area for min/max index calculation - typedef itk::MaskImageFilter< ImageType, MaskImageType, ImageType > MaskImageFilterType; - typename MaskImageFilterType::Pointer masker = MaskImageFilterType::New(); - masker->SetOutsideValue( (statistics.Min+statistics.Max)/2 ); - 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() ); + statistics.Label = (*it); + statistics.N = labelStatisticsFilter->GetCount( *it ); + statistics.Min = labelStatisticsFilter->GetMinimum( *it ); + statistics.Max = labelStatisticsFilter->GetMaximum( *it ); + statistics.Mean = labelStatisticsFilter->GetMean( *it ); + statistics.Median = labelStatisticsFilter->GetMedian( *it ); + statistics.Sigma = labelStatisticsFilter->GetSigma( *it ); + statistics.RMS = sqrt( statistics.Mean * statistics.Mean + + statistics.Sigma * statistics.Sigma ); + + // restrict image to mask area for min/max index calculation + typedef itk::MaskImageFilter< ImageType, MaskImageType, ImageType > MaskImageFilterType; + typename MaskImageFilterType::Pointer masker = MaskImageFilterType::New(); + masker->SetOutsideValue( (statistics.Min+statistics.Max)/2 ); + 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() ); statistics.MinIndex.set_size(adaptedImage->GetImageDimension()); statistics.MaxIndex.set_size(adaptedImage->GetImageDimension()); typename MinMaxFilterType::IndexType tempMaxIndex = minMaxFilter->GetIndexOfMaximum(); typename MinMaxFilterType::IndexType tempMinIndex = minMaxFilter->GetIndexOfMinimum(); - FIX BUG 14644 - If a PlanarFigure is used for segmentation the - adaptedImage is a single slice (2D). Adding the - 3. dimension. +// FIX BUG 14644 + //If a PlanarFigure is used for segmentation the + //adaptedImage is a single slice (2D). Adding the + // 3. dimension. if (m_MaskingMode == MASKING_MODE_PLANARFIGURE && m_Image->GetDimension()==3) { - statistics.MaxIndex.set_size(m_Image->GetDimension()); - statistics.MaxIndex[m_PlanarFigureCoordinate0]=tempMaxIndex[0]; - statistics.MaxIndex[m_PlanarFigureCoordinate1]=tempMaxIndex[1]; - statistics.MaxIndex[m_PlanarFigureAxis]=m_PlanarFigureSlice; - - statistics.MinIndex.set_size(m_Image->GetDimension()); - statistics.MinIndex[m_PlanarFigureCoordinate0]=tempMinIndex[0]; - statistics.MinIndex[m_PlanarFigureCoordinate1]=tempMinIndex[1]; - statistics.MinIndex[m_PlanarFigureAxis]=m_PlanarFigureSlice; + statistics.MaxIndex.set_size(m_Image->GetDimension()); + statistics.MaxIndex[m_PlanarFigureCoordinate0]=tempMaxIndex[0]; + statistics.MaxIndex[m_PlanarFigureCoordinate1]=tempMaxIndex[1]; + statistics.MaxIndex[m_PlanarFigureAxis]=m_PlanarFigureSlice; + + statistics.MinIndex.set_size(m_Image->GetDimension()); + statistics.MinIndex[m_PlanarFigureCoordinate0]=tempMinIndex[0]; + statistics.MinIndex[m_PlanarFigureCoordinate1]=tempMinIndex[1]; + statistics.MinIndex[m_PlanarFigureAxis]=m_PlanarFigureSlice; } else { for (int i = 0; i ThreeDimensionalImageType; + //ThreeDimensionalImageType::Pointer image = ThreeDimensionalImageType::New(); + + //ThreeDimensionalImageType::SizeType size; + //ThreeDimensionalImageType::SpacingType spacing; + //ThreeDimensionalImageType::IndexType start; + + //for(int i = 0; i < 3; ++i) + //{ + // size[i] = 50.00; + // spacing[i] = 1.00; + // start[i] = 0.00; + //} + + //ThreeDimensionalImageType::RegionType region; + //region.SetIndex(start); + //region.SetSize(size); + + //image->SetRegions(region); + //image->Allocate(); + + //typedef itk::ImageRegionIteratorWithIndex IteratorType; + //IteratorType imageIt(image, region); + + + //for(imageIt.GoToBegin(); !imageIt.IsAtEnd(); ++imageIt) + //{ + // imageIt.Set(15); + //} + + ///*for(unsigned int x = 70; x <= 80; ++x) + //{ + // for(unsigned int y = 70; y <= 80; ++y) + // { + // for(unsigned int z = 70; z <= 80; ++z) + // { + // ThreeDimensionalImageType::IndexType pixelIndex; + // pixelIndex[0] = x; + // pixelIndex[1] = y; + // pixelIndex[2] = z; + + // image->SetPixel(pixelIndex, 50); + // } + // } + //}*/ + + //ThreeDimensionalImageType::IndexType pixelIndex; + //pixelIndex[0] = 25; + //pixelIndex[1] = 25; + //pixelIndex[2] = 25; + + //image->SetPixel(pixelIndex, 20); + + //pixelIndex[0] = 0; + //pixelIndex[1] = 0; + //pixelIndex[2] = 0; + + //image->SetPixel(pixelIndex, 1000); + + + ///*****************************************************Creating Test Mask**********************************************/ + //typedef itk::Image ThreeDimensionalMaskImageType; + //ThreeDimensionalMaskImageType::Pointer mask = ThreeDimensionalMaskImageType::New(); + + //ThreeDimensionalMaskImageType::SizeType maskSize; + //ThreeDimensionalMaskImageType::SpacingType maskSpacing; + //ThreeDimensionalMaskImageType::IndexType maskStart; + + //for(int i = 0; i < 3; ++i) + //{ + // maskSize[i] = 50; + // maskSpacing[i] = 1.00; + // maskStart[i] = 0.00; + //} + + //ThreeDimensionalMaskImageType::RegionType maskRegion; + //maskRegion.SetIndex(maskStart); + //maskRegion.SetSize(maskSize); + + //mask->SetRegions(maskRegion); + //mask->Allocate(); + + //typedef itk::ImageRegionIteratorWithIndex MaskIteratorType; + //MaskIteratorType maskIt(mask, maskRegion); + + //for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt) + //{ + // maskIt.Set(1); + //} + /*****************************************************Calculate Hotspot Statistics**********************************************/ + const double radiusSUVHotspot = 6.2035049089940; // radius of a 1cm3 sphere in a isotrope image of 1mm spacings + + SetHotspotRadius(radiusSUVHotspot); + SetCalculateHotspot(true); + + if(IsHotspotCalculated()) + { + Statistics hotspotStatistics = CalculateHotspotStatistics (adaptedImage.GetPointer(), adaptedMaskImage.GetPointer(), GetHotspotRadius()); + statistics.HotspotMax = hotspotStatistics.HotspotMax; + statistics.HotspotMin = hotspotStatistics.HotspotMin; + statistics.HotspotPeak = hotspotStatistics.HotspotPeak; + statistics.HotspotMaxIndex = hotspotStatistics.HotspotMaxIndex; + statistics.HotspotMinIndex = hotspotStatistics.HotspotMinIndex; + statistics.HotspotPeakIndex = hotspotStatistics.HotspotPeakIndex; + } + statisticsContainer->push_back( statistics ); + } + } + else + { + histogramContainer->push_back( HistogramType::ConstPointer( m_EmptyHistogram ) ); + statisticsContainer->push_back( Statistics() ); + } +} + +template +ImageStatisticsCalculator::MinMaxIndex ImageStatisticsCalculator::CalculateMinMaxIndex( + const itk::Image *inputImage, + itk::Image *maskImage, + float minBoundary, + float maxBoundary) +{ + typedef itk::Image< TPixel, VImageDimension > ImageType; + typedef itk::Image< unsigned short, VImageDimension > MaskImageType; - /*****************************************************Calculate Hotspot Statistics**********************************************/ + typedef itk::ImageRegionConstIterator MaskImageIteratorType; + typedef itk::ImageRegionConstIteratorWithIndex InputImageIndexIteratorType; - } - const double radiusSUVHotspot = 6.2035049089940; // radius of a 1cm3 sphere in a isotrope image of 1mm spacings + MaskImageIteratorType maskIt(maskImage, maskImage->GetLargestPossibleRegion()); + InputImageIndexIteratorType imageIndexIt(inputImage, inputImage->GetLargestPossibleRegion()); + + float maxValue = -itk::NumericTraits::max(); + float minValue = itk::NumericTraits::max(); + + ImageType::IndexType maxIndex; + ImageType::IndexType minIndex; - SetHotspotRadius(radiusSUVHotspot); - SetCalculateHotspot(true); + double maxBound = itk::NumericTraits::min(); + double minBound = itk::NumericTraits::max(); - if(IsHotspotCalculated()) + bool sameIndex = true; + + for(maskIt.GoToBegin(), imageIndexIt.GoToBegin(); + !maskIt.IsAtEnd() && !imageIndexIt.IsAtEnd(); + ++maskIt, ++imageIndexIt) + { + if(maskIt.Get() > itk::NumericTraits::Zero) + { + double value = imageIndexIt.Get(); + ImageType::IndexType index = imageIndexIt.GetIndex(); + + //Calculate minimum, maximum and corresponding index-values + if( value > maxValue && value < maxBoundary) { - Statistics hotspotStatistics = CalculateHotspotStatistics(adaptedImage.GetPointer(), adaptedMaskImage.GetPointer(), GetHotspotRadius()); - statistics.HotspotMax = hotspotStatistics.HotspotMax; - statistics.HotspotMin = hotspotStatistics.HotspotMin; - statistics.HotspotPeak = hotspotStatistics.HotspotPeak; - statistics.HotspotMaxIndex = hotspotStatistics.HotspotMaxIndex; - statistics.HotspotMinIndex = hotspotStatistics.HotspotMinIndex; - statistics.HotspotPeakIndex = hotspotStatistics.HotspotPeakIndex; + maxIndex = index; + maxValue = value; } - statisticsContainer->push_back( statistics ); + if(value < minValue && value > minBoundary) + { + minIndex = index; + minValue = value; + } + + } + } + + + MinMaxIndex minMax; + + minMax.MinIndex.set_size(inputImage->GetImageDimension()); + minMax.MaxIndex.set_size(inputImage->GetImageDimension()); + + for(int i = 0; i < minMax.MaxIndex.size(); ++i) + minMax.MaxIndex[i] = maxIndex[i]; + + for(int i = 0; i < minMax.MinIndex.size(); ++i) + minMax.MinIndex[i] = minIndex[i]; + + + minMax.Max = maxValue; + minMax.Min = minValue; + + return minMax; +} + +template < typename TPixel, unsigned int VImageDimension> +ImageStatisticsCalculator::Statistics ImageStatisticsCalculator::CalculateHotspotStatistics( + const itk::Image *inputImage, + itk::Image *maskImage, + double radiusInMM) +{ + typedef itk::Image< TPixel, VImageDimension > ImageType; + typedef itk::Image< float, VImageDimension > MaskImageType; + typedef itk::Image SphereMaskImageType; + + typedef itk::ContinuousIndex ContinuousIndexType; + + typedef itk::ImageRegionConstIteratorWithIndex IteratorType; + typedef itk::ImageRegionConstIteratorWithIndex MaskIteratorType; + typedef itk::ImageRegionConstIteratorWithIndex SphereMaskIteratorType; + + typedef typename ImageType::SpacingType SpacingType; + typedef typename ImageType::SizeType SizeType; + typedef typename ImageType::IndexType IndexType; + typedef typename ImageType::PointType PointType; + typedef typename ImageType::RegionType RegionType; + + itk::TimeProbesCollectorBase chronometer; + + MaskImageType::Pointer convolutionMask = MaskImageType::New(); + + SpacingType spacing = inputImage->GetSpacing(); + + IndexType start; + start.Fill(0); + + ContinuousIndexType origin; + SizeType size; + ContinuousIndexType convolutionMaskCenterCoordinate; + /*****************************************************Creating convolution mask**********************************************/ + chronometer.Start("Mask"); + for(unsigned int i = 0; i < VImageDimension; ++i) + { + origin[i] = 0.0; + + double countIndex = 2.0 * radiusInMM / spacing[i]; + + // Rounding up to the next integer by cast + countIndex += 0.9; + int castedIndex = static_cast(countIndex); + + // We always have an uneven number in size to determine a center-point in the convolution mask + if(castedIndex % 2 > 0 ) + { + size[i] = castedIndex; } else { - histogramContainer->push_back( HistogramType::ConstPointer( m_EmptyHistogram ) ); - statisticsContainer->push_back( Statistics() ); + size[i] = castedIndex +1; } + + convolutionMaskCenterCoordinate[i] = (size[i] -1) / 2; } - template - ImageStatisticsCalculator::MinMaxIndex ImageStatisticsCalculator::CalculateMinMaxIndex( - const itk::Image *inputImage, - itk::Image *maskImage) - { - typedef itk::Image< TPixel, VImageDimension > ImageType; - typedef itk::Image< unsigned short, VImageDimension > MaskImageType; - typedef itk::ImageRegionConstIterator MaskImageIteratorType; - typedef itk::ImageRegionConstIteratorWithIndex InputImageIndexIteratorType; + RegionType region; + region.SetSize(size); + region.SetIndex(start); + + convolutionMask->SetRegions(region); + convolutionMask->SetSpacing(spacing); + convolutionMask->Allocate(); - MaskImageIteratorType maskIt(maskImage, maskImage->GetLargestPossibleRegion()); - InputImageIndexIteratorType imageIndexIt(inputImage, inputImage->GetLargestPossibleRegion()); - float maxValue = itk::NumericTraits::min(); - float minValue = itk::NumericTraits::max(); + float subPixelDimensionX = spacing[0]/2; + float subPixelDimensionY = spacing[1]/2; - ImageType::IndexType maxIndex; - ImageType::IndexType minIndex; + double distance = 0.0; + double countSubPixel = 0.0; + double pixelValue = 0.0; - for(maskIt.GoToBegin(), imageIndexIt.GoToBegin(); - !maskIt.IsAtEnd() && !imageIndexIt.IsAtEnd(); - ++maskIt, ++imageIndexIt) + MaskIteratorType maskIt(convolutionMask,region); + + // Generate convolutionMask + for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt) + { + ContinuousIndexType indexPoint(maskIt.GetIndex()); + + for(double x = indexPoint[0] - subPixelDimensionX / 2; + x <= indexPoint[0] + subPixelDimensionX / 2; + x += subPixelDimensionX) { - if(maskIt.Get() > itk::NumericTraits::Zero) + for(double y = indexPoint[1] - subPixelDimensionY /2; + y <= indexPoint[1] + subPixelDimensionY / 2; + y += subPixelDimensionY) { - double value = imageIndexIt.Get(); - - Calculate minimum, maximum and corresponding index-values - if( value > maxValue) + ContinuousIndexType subPixelCenterCoordinate; + subPixelCenterCoordinate[0] = x; + subPixelCenterCoordinate[1] = y; + for( unsigned int i = 2; i < VImageDimension; ++i ) { - maxValue = value; - ImageType::IndexType index = imageIndexIt.GetIndex(); - maxIndex = index; + subPixelCenterCoordinate[i] = indexPoint[i]; } + PointType subPixelCenterCoordinateInPhysicalPoint; + convolutionMask->TransformContinuousIndexToPhysicalPoint(subPixelCenterCoordinate, subPixelCenterCoordinateInPhysicalPoint); - if(value < minValue) - { - minValue = value; - ImageType::IndexType index = imageIndexIt.GetIndex(); - minIndex = index; - } + PointType convolutionMaskCenterCoordinateInPhysicalPoint; + convolutionMask->TransformContinuousIndexToPhysicalPoint(convolutionMaskCenterCoordinate, convolutionMaskCenterCoordinateInPhysicalPoint); + + distance = subPixelCenterCoordinateInPhysicalPoint.EuclideanDistanceTo(convolutionMaskCenterCoordinateInPhysicalPoint); + + if(distance <= radiusInMM) + countSubPixel++; } } - MinMaxIndex minMax; + // pixelValue is the counted SubPixels divided by factor 4 + pixelValue = countSubPixel / 4.00; + convolutionMask->SetPixel(maskIt.GetIndex(),pixelValue); - minMax.MinIndex.set_size(inputImage->GetImageDimension()); - minMax.MaxIndex.set_size(inputImage->GetImageDimension()); + countSubPixel = 0.00; + } - Converting IndexType in vnl_vector - for(int i = 0; i < minMax.MaxIndex.size(); ++i) - minMax.MaxIndex[i] = maxIndex[i]; + chronometer.Stop("Mask"); + /*****************************************************Creating Peak Image**********************************************/ + chronometer.Start("Convolution"); + typedef itk::ConvolutionImageFilter FilterType; + FilterType::Pointer convolutionFilter = FilterType::New(); - for(int i = 0; i < minMax.MinIndex.size(); ++i) - minMax.MinIndex[i] = minIndex[i]; + typedef itk::ConstantBoundaryCondition BoundaryConditionType; + BoundaryConditionType* boundaryCondition = new BoundaryConditionType(); + boundaryCondition->SetConstant(0.0); + convolutionFilter->SetBoundaryCondition(boundaryCondition); + convolutionFilter->SetInput(inputImage); + convolutionFilter->SetKernelImage(convolutionMask); + convolutionFilter->SetNormalize(true); + convolutionFilter->Update(); + chronometer.Stop("Convolution"); + MaskImageType::Pointer peakImage = convolutionFilter->GetOutput(); - minMax.Max = maxValue; - minMax.Min = minValue; + peakImage->SetSpacing( inputImage->GetSpacing() ); // TODO: only workaround because convolution filter seems to ignore spacing of input image - return minMax; - } - template < typename TPixel, unsigned int VImageDimension> - ImageStatisticsCalculator::Statistics ImageStatisticsCalculator::CalculateHotspotStatistics( - const itk::Image *inputImage, - itk::Image *maskImage, - double radiusInMM) - { - typedef itk::Image< TPixel, VImageDimension > ImageType; - typedef itk::Image ConvolutionMaskImageType; - typedef itk::Image SphereMaskImageType; + /*****************************************************Creating Hotspot Sphere**********************************************/ + SphereMaskImageType::Pointer hotspotSphere = SphereMaskImageType::New(); - itk::TimeProbesCollectorBase chronometer; + SphereMaskImageType::SpacingType hotspotSphereSpacing = peakImage->GetSpacing(); + SphereMaskImageType::SizeType hotspotSphereSize; + SphereMaskImageType::RegionType regionSphere; + SphereMaskImageType::IndexType hotspotPeakIndex; + PointType hotspotOrigin; + IndexType offsetInIndex; - typedef typename ImageType::SpacingType SpacingType; - SpacingType spacing = inputImage->GetSpacing(); + typedef itk::EllipseSpatialObject EllipseType; + typedef itk::SpatialObjectToImageFilter SpatialObjectToImageFilter; - typedef typename ConvolutionMaskImageType::IndexType IndexType; - IndexType start; - start.Fill(0); + float maxBoundary = itk::NumericTraits::max(); + float minBoundary = itk::NumericTraits::min(); + double hotspotPeak = itk::NumericTraits::min(); - typedef typename ConvolutionMaskImageType::SizeType SizeType; - SizeType size; + bool isInside = false; - typedef itk::ContinuousIndex ContinuousIndexType; - ContinuousIndexType convolutionMaskCenterCoordinate; + SphereMaskImageType::Pointer croppedRegionMask = SphereMaskImageType::New(); + SphereMaskImageType::SpacingType maskSpacing = peakImage->GetSpacing(); - /*****************************************************Creating convolution mask**********************************************/ - chronometer.Start("Mask"); - for(unsigned int i = 0; i < VImageDimension; ++i) - { - double countIndex = 2.0 * radiusInMM / spacing[i]; + SphereMaskImageType::RegionType maskRegion; + SphereMaskImageType::IndexType maskStart; + maskStart.Fill(0); + SphereMaskImageType::SizeType maskSize = peakImage->GetLargestPossibleRegion().GetSize(); - Rounding up to the next integer by cast - countIndex += 0.99999; - int castedIndex = static_cast(countIndex); + maskRegion.SetIndex(maskStart); + maskRegion.SetSize(peakImage->GetLargestPossibleRegion().GetSize()); + croppedRegionMask->SetRegions(maskRegion); - We always have an uneven number in size to determine a center-point in the convolution mask - if(castedIndex % 2 > 0 ) - { - size[i] = castedIndex; - } - else + int offsetX = static_cast((radiusInMM / maskSpacing[0]) + 0.99999); + int offsetY = static_cast((radiusInMM / maskSpacing[1]) + 0.99999); + int offsetZ = static_cast((radiusInMM / maskSpacing[2]) + 0.99999); + + croppedRegionMask->Allocate(); + + typedef itk::ImageRegionIteratorWithIndex CroppedImageIteratorType; + CroppedImageIteratorType sphereMaskIt(croppedRegionMask, croppedRegionMask->GetLargestPossibleRegion()); + + for(sphereMaskIt.GoToBegin(); !sphereMaskIt.IsAtEnd(); ++sphereMaskIt) + { + sphereMaskIt.Set(0); + } + + for(unsigned int x = offsetX ; x < maskSize[0] - offsetX; ++x) + { + for(unsigned int y = offsetY; y < maskSize[1] - offsetY; ++y) + { + for(unsigned int z = offsetZ; z < maskSize[2] - offsetZ; ++z) { - size[i] = castedIndex +1; + SphereMaskImageType::IndexType pixelIndex; + pixelIndex[0] = x; + pixelIndex[1] = y; + pixelIndex[2] = z; + + croppedRegionMask->SetPixel(pixelIndex, 1); } - convolutionMaskCenterCoordinate[i] = (size[i] -1) / 2; } + } - typedef typename ConvolutionMaskImageType::RegionType RegionType; - RegionType region; - region.SetSize(size); - region.SetIndex(start); - - ConvolutionMaskImageType::Pointer convolutionMask = ConvolutionMaskImageType::New(); + chronometer.Start("Find Peak"); + MinMaxIndex peakInformations = CalculateMinMaxIndex(peakImage.GetPointer(), croppedRegionMask.GetPointer(), minBoundary, maxBoundary); + chronometer.Stop("Find Peak"); - convolutionMask->SetRegions(region); - convolutionMask->SetSpacing(spacing); - convolutionMask->Allocate(); + hotspotPeak = peakInformations.Max; - float subPixelDimensionX = spacing[0]/2; - float subPixelDimensionY = spacing[1]/2; + for(int i = 0; i < VImageDimension; ++i) + hotspotPeakIndex[i] = peakInformations.MaxIndex[i]; - double distance = 0.0; - double countSubPixel = 0.0; - double pixelValue = 0.0; + for(unsigned int i = 0; i < VImageDimension; ++i) + { - typedef typename ConvolutionMaskImageType::PointType PointType; - typedef itk::ImageRegionConstIteratorWithIndex ConvolutionMaskIteratorType; - ConvolutionMaskIteratorType maskIt(convolutionMask,region); + double countIndex = 2.0 * radiusInMM / hotspotSphereSpacing[i]; - Generate convolutionMask + // Rounding up to the next integer by cast + countIndex += 0.9999999; + int castedIndex = static_cast(countIndex); - for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt) + // We always have an uneven number in size to determine a center-point in the convolution mask + if(castedIndex % 2 > 0 ) { - ContinuousIndexType indexPoint(maskIt.GetIndex()); - - for(double x = indexPoint[0] - 0.25; - x <= indexPoint[0] + 0.25; - x += 0.5) - { - for(double y = indexPoint[1] - 0.25; - y <= indexPoint[1] + 0.25; - y += 0.5) - { - ContinuousIndexType subPixelCenterCoordinate; - subPixelCenterCoordinate[0] = x; - subPixelCenterCoordinate[1] = y; - for( unsigned int i = 2; i < VImageDimension; ++i ) - { - subPixelCenterCoordinate[i] = indexPoint[i]; - } + hotspotSphereSize[i] = castedIndex; + } + else + { + hotspotSphereSize[i] = castedIndex +1; + } + } - PointType subPixelCenterCoordinateInPhysicalPoint; - convolutionMask->TransformContinuousIndexToPhysicalPoint(subPixelCenterCoordinate, subPixelCenterCoordinateInPhysicalPoint); + // Initialize SpatialObjectoToImageFilter + itk::SpatialObjectToImageFilter::Pointer spatialObjectToImageFilter + = SpatialObjectToImageFilter::New(); - PointType convolutionMaskCenterCoordinateInPhysicalPoint; - convolutionMask->TransformContinuousIndexToPhysicalPoint(convolutionMaskCenterCoordinate, convolutionMaskCenterCoordinateInPhysicalPoint); + spatialObjectToImageFilter->SetSize(hotspotSphereSize); + spatialObjectToImageFilter->SetSpacing(hotspotSphereSpacing); - distance = subPixelCenterCoordinateInPhysicalPoint.EuclideanDistanceTo(convolutionMaskCenterCoordinateInPhysicalPoint); + // Creating spatial sphere object + EllipseType::Pointer sphere = EllipseType::New(); + sphere->SetRadius(radiusInMM); + typedef EllipseType::TransformType TransformType; + TransformType::Pointer transform = TransformType::New(); - if(distance <= radiusInMM) - countSubPixel++; - } - } - pixelValue is the counted SubPixels divided by factor 4 - pixelValue = countSubPixel / 4.00; - convolutionMask->SetPixel(maskIt.GetIndex(),pixelValue); + transform->SetIdentity(); - countSubPixel = 0.00; - } + TransformType::OutputVectorType translation; - chronometer.Stop("Mask"); - /*****************************************************Creating Peak Image**********************************************/ - chronometer.Start("Convolution"); + // Transform sphere on center-position, set pixelValues inside sphere on 1 and update + for(int i = 0; i < VImageDimension; ++i) + translation[i] = (hotspotSphereSize[i] -1) / 2; - typedef itk::Image< float, VImageDimension > PeakImageType; - typedef itk::FFTConvolutionImageFilter FilterType; - FilterType::Pointer convolutionFilter = FilterType::New(); + transform->Translate(translation, false); - typedef itk::ConstantBoundaryCondition BoundaryConditionType; - BoundaryConditionType* boundaryCondition = new BoundaryConditionType(); - boundaryCondition->SetConstant(0.0); + sphere->SetObjectToParentTransform(transform); - convolutionFilter->SetBoundaryCondition(boundaryCondition); - convolutionFilter->SetInput(inputImage); - convolutionFilter->SetKernelImage(convolutionMask.GetPointer()); - convolutionFilter->SetNormalize(true); - convolutionFilter->Update(); - chronometer.Stop("Convolution"); - PeakImageType::Pointer peakImage = convolutionFilter->GetOutput(); + spatialObjectToImageFilter->SetInput(sphere); - delete boundaryCondition; + sphere->SetDefaultInsideValue(1.00); + sphere->SetDefaultOutsideValue(0.00); - peakImage->SetSpacing( inputImage->GetSpacing() ); // TODO: only workaround because convolution filter seems to ignore spacing of input image + spatialObjectToImageFilter->SetUseObjectValue(true); + spatialObjectToImageFilter->SetOutsideValue(0); + spatialObjectToImageFilter->Update(); + hotspotSphere = spatialObjectToImageFilter->GetOutput(); - /*****************************************************Creating Hotspot Sphere**********************************************/ - typedef itk::Image SphereMaskImageType; - SphereMaskImageType::Pointer hotspotSphere = SphereMaskImageType::New(); + // Calculate new origin for hotspot sphere - SphereMaskImageType::SpacingType hotspotSphereSpacing = peakImage->GetSpacing(); - SphereMaskImageType::SizeType hotspotSphereSize; - SphereMaskImageType::RegionType regionSphere; - SphereMaskImageType::IndexType hotspotPeakIndex; - PointType hotspotOrigin; - IndexType offsetInIndex; + for(int i = 0; i < VImageDimension; ++i) + offsetInIndex[i] = (hotspotSphereSize[i] -1) / 2; - typedef itk::EllipseSpatialObject EllipseType; - typedef itk::SpatialObjectToImageFilter SpatialObjectToImageFilter; + peakImage->TransformIndexToPhysicalPoint(hotspotPeakIndex, hotspotOrigin); - float hotspotPeak = itk::NumericTraits::min(); + PointType offsetInPhysicalPoint; + hotspotSphere->TransformIndexToPhysicalPoint(offsetInIndex, offsetInPhysicalPoint); - SphereMaskImageType::Pointer croppedRegionMask = SphereMaskImageType::New(); - SphereMaskImageType::SpacingType maskSpacing = peakImage->GetSpacing(); + for(int i = 0; i < VImageDimension; ++i) + hotspotOrigin[i] -= offsetInPhysicalPoint[i]; - SphereMaskImageType::IndexType maskStart; - maskStart.Fill(0); - SphereMaskImageType::SizeType maskSize = peakImage->GetLargestPossibleRegion().GetSize(); - SphereMaskImageType::RegionType maskRegion; + hotspotSphere->SetOrigin(hotspotOrigin); + hotspotSphere->Allocate(); - maskRegion.SetIndex(maskStart); - maskRegion.SetSize(peakImage->GetLargestPossibleRegion().GetSize()); +#ifdef DEBUG_HOTSPOTSEARCH - croppedRegionMask->SetRegions(maskRegion); + std::cout << std::endl << std::endl; + std::cout << "hotspotMask: " << std::endl; + unsigned int lastZ = 1000000000; + unsigned int lastY = 1000000000; - int offsetX = static_cast((radiusInMM / maskSpacing[0]) + 0.9999 ); - int offsetY = static_cast((radiusInMM / maskSpacing[1]) + 0.9999 ); - int offsetZ = static_cast((radiusInMM / maskSpacing[2]) + 0.9999 ); + unsigned int hotspotMaskIndexCounter = 0; + SphereMaskIteratorType hotspotMaskIt(hotspotSphere, hotspotSphere->GetLargestPossibleRegion() ); + for(hotspotMaskIt.GoToBegin();!hotspotMaskIt.IsAtEnd();++hotspotMaskIt) + { - croppedRegionMask->Allocate(); + double tmp = hotspotMaskIt.Get(); + if (hotspotMaskIt.GetIndex()[1] != lastY) + { + std::cout << std::endl; + lastY = hotspotMaskIt.GetIndex()[1]; + } + if (hotspotMaskIt.GetIndex()[0] != lastZ) + { + std::cout << tmp << " "; + lastZ = hotspotMaskIt.GetIndex()[0]; + } - typedef itk::ImageRegionIteratorWithIndex CroppedImageIteratorType; - CroppedImageIteratorType sphereMaskIt(croppedRegionMask, croppedRegionMask->GetLargestPossibleRegion()); + hotspotMaskIndexCounter++; - Creating mask which gurantees that sphere is inside the given region in any case - for(sphereMaskIt.GoToBegin(); !sphereMaskIt.IsAtEnd(); ++sphereMaskIt) - { - SphereMaskImageType::IndexType index = sphereMaskIt.GetIndex(); + if(hotspotMaskIndexCounter > 168) { + std::cout << std::endl; + hotspotMaskIndexCounter = 0; + } + } - if((index[0] >= offsetX && index[0] <= maskSize[0] - offsetX -1 ) && - (index[1] >= offsetY && index[1] <= maskSize[1] - offsetY -1) && - (index[2] >= offsetZ && index[2] <= maskSize[2] - offsetZ -1)) - sphereMaskIt.Set(1); - else - sphereMaskIt.Set(0); - } + std::cout << std::endl << std::endl; +#endif - chronometer.Start("Find Peak"); - MinMaxIndex peakInformations = CalculateMinMaxIndex(peakImage.GetPointer(), croppedRegionMask.GetPointer()); - chronometer.Stop("Find Peak"); + /*********************************Creating cropped inputImage for calculation of hotspot statistics****************************/ - hotspotPeak = peakInformations.Max; + IndexType croppedStart; + ImageType::RegionType inputRegion; + peakImage->TransformPhysicalPointToIndex(hotspotOrigin,croppedStart); + ImageType::RegionType::IndexType inputStart = croppedStart; + ImageType::RegionType::SizeType inputSize = hotspotSphere->GetLargestPossibleRegion().GetSize(); + inputRegion.SetIndex(inputStart); + inputRegion.SetSize(inputSize); - for(unsigned int i = 0; i < VImageDimension; ++i) - { - hotspotPeakIndex[i] = peakInformations.MaxIndex[i]; + ImageType::RegionType croppedOutputRegion; + ImageType::IndexType croppedOutputStart; - double countIndex = 2.0 * radiusInMM / hotspotSphereSpacing[i]; + croppedOutputStart.Fill(0); + croppedOutputRegion.SetIndex(croppedOutputStart); + croppedOutputRegion.SetSize(hotspotSphere->GetLargestPossibleRegion().GetSize()); - Rounding up to the next integer by cast - int castedIndex = static_cast(countIndex); + ImageType::Pointer croppedOutputImage = ImageType::New(); + croppedOutputImage->SetRegions(croppedOutputRegion); + croppedOutputImage->Allocate(); - We always have an uneven number in size to determine a center-point in the convolution mask - if(castedIndex % 2 > 0 ) - hotspotSphereSize[i] = castedIndex; - else - hotspotSphereSize[i] = castedIndex +1; - } + typedef itk::ImageRegionConstIterator ImageIteratorType; + ImageIteratorType inputIt(inputImage, inputRegion); - Initialize SpatialObjectoToImageFilter - itk::SpatialObjectToImageFilter::Pointer spatialObjectToImageFilter - = SpatialObjectToImageFilter::New(); + ImageIteratorType croppedOutputImageIt(croppedOutputImage, croppedOutputRegion); + for(inputIt.GoToBegin(), croppedOutputImageIt.GoToBegin(); !inputIt.IsAtEnd(); ++inputIt, ++croppedOutputImageIt) + { + croppedOutputImage->SetPixel(croppedOutputImageIt.GetIndex(), inputIt.Get()); + } - spatialObjectToImageFilter->SetSize(hotspotSphereSize); - spatialObjectToImageFilter->SetSpacing(hotspotSphereSpacing); + // Calculate statistics in Hotspot + maxBoundary = itk::NumericTraits::max(); + minBoundary = itk::NumericTraits::min(); + MinMaxIndex hotspotInformations; + Statistics hotspotStatistics; - Creating spatial sphere object - EllipseType::Pointer sphere = EllipseType::New(); - sphere->SetRadius(radiusInMM); - typedef EllipseType::TransformType TransformType; - TransformType::Pointer transform = TransformType::New(); + chronometer.Start("Hotspot statistics"); + hotspotInformations = CalculateMinMaxIndex(croppedOutputImage.GetPointer(), hotspotSphere.GetPointer(), minBoundary, maxBoundary); + chronometer.Stop("Hotspot statistics"); - transform->SetIdentity(); - TransformType::OutputVectorType translation; + // Add offset to cropped indices + for(int i = 0; i < VImageDimension; ++i) + { + hotspotInformations.MaxIndex[i] += croppedStart[i]; + hotspotInformations.MinIndex[i] += croppedStart[i]; + } - Transform sphere on center-position, set pixelValues inside sphere on 1 and update + hotspotStatistics.HotspotMin = hotspotInformations.Min; + hotspotStatistics.HotspotMinIndex = hotspotInformations.MinIndex; + hotspotStatistics.HotspotMax = hotspotInformations.Max; + hotspotStatistics.HotspotMaxIndex = hotspotInformations.MaxIndex; + hotspotStatistics.HotspotPeak = hotspotPeak; - for(int i = 0; i < VImageDimension; ++i) - translation[i] = ((hotspotSphereSize[i] ) / 2) * hotspotSphereSpacing[i]; + hotspotStatistics.HotspotPeakIndex.set_size(inputImage->GetImageDimension()); + for (int i = 0; i< hotspotStatistics.HotspotPeakIndex.size(); ++i) + { + hotspotStatistics.HotspotPeakIndex[i] = hotspotPeakIndex[i]; + } - transform->Translate(translation, false); + chronometer.Report( std::cout ); + return hotspotStatistics; +} - sphere->SetObjectToParentTransform(transform); +template < typename TPixel, unsigned int VImageDimension> +bool ImageStatisticsCalculator::IsSphereInsideRegion(const itk::Image *inputImage, + itk::Image *sphereImage) +{ + typedef itk::Image< TPixel, VImageDimension > ImageType; + typedef ImageType::RegionType RegionType; + typedef ImageType::SizeType SizeType; + typedef ImageType::PointType PointType; + typedef ImageType::IndexType IndexType; - spatialObjectToImageFilter->SetInput(sphere); + RegionType region(inputImage->GetLargestPossibleRegion().GetIndex(), inputImage->GetLargestPossibleRegion().GetSize() ); - sphere->SetDefaultInsideValue(1.00); - sphere->SetDefaultOutsideValue(0.00); + typedef itk::Image SphereImageType; + typedef SphereImageType::SizeType SphereSizeType; + typedef SphereImageType::IndexType SphereIndexType; - spatialObjectToImageFilter->SetUseObjectValue(true); - spatialObjectToImageFilter->SetOutsideValue(0); + SphereIndexType sphereCenterIndex; + SphereSizeType sphereSize = sphereImage->GetLargestPossibleRegion().GetSize(); - spatialObjectToImageFilter->Update(); - hotspotSphere = spatialObjectToImageFilter->GetOutput(); + double radius = (sphereSize[0] - 1)/2.00; + bool isInside = false; - Calculate new origin for hotspot sphere + for(int i = 0; i < VImageDimension; ++i) + sphereCenterIndex[i] = (sphereSize[i] -1) / 2; - for(int i = 0; i < VImageDimension; ++i) - offsetInIndex[i] = (hotspotSphereSize[i] -1) / 2; - peakImage->TransformIndexToPhysicalPoint(hotspotPeakIndex, hotspotOrigin); + // If the sphere is inside the given region is checked by 14 points which are on the surface of the sphere - PointType offsetInPhysicalPoint; - hotspotSphere->TransformIndexToPhysicalPoint(offsetInIndex, offsetInPhysicalPoint); + /**********************************Calculation of sphere-points which are on axis of coordinates****************************/ - for(int i = 0; i < VImageDimension; ++i) - hotspotOrigin[i] -= offsetInPhysicalPoint[i]; + // First 6 points - hotspotSphere->SetOrigin(hotspotOrigin); - hotspotSphere->Allocate(); + for(int i = 0; i < VImageDimension; ++i) + { + for(double j = sphereCenterIndex[i] - radius; j <= sphereCenterIndex[i]+ radius; j += 2*radius) + { + IndexType sphereIndex; + if(i == 0) + { + sphereIndex[0] = j; + sphereIndex[1] = sphereCenterIndex[1]; + sphereIndex[2] = sphereCenterIndex[2]; + } - /*********************************Creating cropped inputImage for calculation of hotspot statistics****************************/ + if(i == 1) + { + sphereIndex[0] = sphereCenterIndex[0]; + sphereIndex[1] = j; + sphereIndex[2] = sphereCenterIndex[2]; + } - ImageType::IndexType croppedStart; - ImageType::RegionType croppedRegion; + if(i == 2) + { + sphereIndex[0] = sphereCenterIndex[0]; + sphereIndex[1] = sphereCenterIndex[1]; + sphereIndex[2] = j; + } - peakImage->TransformPhysicalPointToIndex(hotspotOrigin,croppedStart); - ImageType::RegionType::SizeType croppedSize = hotspotSphere->GetLargestPossibleRegion().GetSize(); + PointType spherePoint; + sphereImage->TransformIndexToPhysicalPoint(sphereIndex, spherePoint); - croppedRegion.SetIndex(croppedStart); - croppedRegion.SetSize(croppedSize); + IndexType regionIndex; + inputImage->TransformPhysicalPointToIndex(spherePoint, regionIndex); - ImageType::RegionType croppedOutputRegion; - ImageType::IndexType croppedOutputStart; + if(region.IsInside(regionIndex) == true) { + isInside = true; + } + else + { + isInside = false; + return isInside; + } + } + } + /**********************************Calculation of sphere-coordinates which are not on axis of coordinates************************/ + IndexType sphereIndex; + IndexType regionIndex; + PointType spherePoint; - croppedOutputStart.Fill(0); - croppedOutputRegion.SetIndex(croppedOutputStart); - croppedOutputRegion.SetSize(croppedSize); + // First run: summation, second run: subtraction + bool addOrSub = false; - ImageType::Pointer croppedOutputImage = ImageType::New(); - croppedOutputImage->SetRegions(croppedOutputRegion); - croppedOutputImage->Allocate(); + int offsetXDirection = static_cast(radius * sin(45.00 * itk::Math::pi / 180.00) * cos(45.00 * itk::Math::pi / 180.00)); + int offsetYDirection = static_cast(radius * sin(45.00 * itk::Math::pi / 180.00) * sin(45.00 * itk::Math::pi / 180.00)); + int offsetZDirection = static_cast(radius * cos(45.00 * itk::Math::pi / 180.00)); - typedef itk::ImageRegionConstIterator ImageIteratorType; - ImageIteratorType inputIt(inputImage, croppedRegion); + for(int i = 0; i < 2; ++i) + { + if(addOrSub == false) + { + // Point 7 + sphereIndex[0] = sphereCenterIndex[0] + offsetXDirection; + sphereIndex[1] = sphereCenterIndex[1] + offsetYDirection; + sphereIndex[2] = sphereCenterIndex[2] + offsetZDirection; - ImageIteratorType croppedOutputImageIt(croppedOutputImage, croppedOutputRegion); + sphereImage->TransformIndexToPhysicalPoint(sphereIndex, spherePoint); + inputImage->TransformPhysicalPointToIndex(spherePoint, regionIndex); - for(inputIt.GoToBegin(), croppedOutputImageIt.GoToBegin(); !inputIt.IsAtEnd(); ++inputIt, ++croppedOutputImageIt) + if(region.IsInside(regionIndex) == true) { + isInside = true; + } + else + { + isInside = false; + return isInside; + } + } + else { - croppedOutputImage->SetPixel(croppedOutputImageIt.GetIndex(), inputIt.Get()); + // Point 8 + sphereIndex[0] = sphereCenterIndex[0] + offsetXDirection; + sphereIndex[1] = sphereCenterIndex[1] - offsetYDirection; + sphereIndex[2] = sphereCenterIndex[2] + offsetZDirection; + + sphereImage->TransformIndexToPhysicalPoint(sphereIndex, spherePoint); + inputImage->TransformPhysicalPointToIndex(spherePoint, regionIndex); + + if(region.IsInside(regionIndex) == true) + { + isInside = true; + } + else + { + isInside = false; + return isInside; + } } - Calculate statistics in Hotspot - MinMaxIndex hotspotInformations; - Statistics hotspotStatistics; + // Point 9/10 + sphereIndex[2] = sphereCenterIndex[2] - offsetZDirection; - chronometer.Start("Hotspot statistics"); - hotspotInformations = CalculateMinMaxIndex(croppedOutputImage.GetPointer(), hotspotSphere.GetPointer()); - chronometer.Stop("Hotspot statistics"); + sphereImage->TransformIndexToPhysicalPoint(sphereIndex, spherePoint); + inputImage->TransformPhysicalPointToIndex(spherePoint, regionIndex); + if(region.IsInside(regionIndex) == true) + { + isInside = true; + } + else + { + isInside = false; + return isInside; + } - Add offset to cropped indices - for(int i = 0; i < VImageDimension; ++i) - { - hotspotInformations.MaxIndex[i] += croppedStart[i]; - hotspotInformations.MinIndex[i] += croppedStart[i]; - } - hotspotStatistics.HotspotMin = hotspotInformations.Min; - hotspotStatistics.HotspotMinIndex = hotspotInformations.MinIndex; - hotspotStatistics.HotspotMax = hotspotInformations.Max; - hotspotStatistics.HotspotMaxIndex = hotspotInformations.MaxIndex; - hotspotStatistics.HotspotPeak = hotspotPeak; + // Point 11/12 + sphereIndex[0] = sphereCenterIndex[0] - offsetXDirection; - hotspotStatistics.HotspotPeakIndex.set_size(inputImage->GetImageDimension()); - for (int i = 0; i< hotspotStatistics.HotspotPeakIndex.size(); ++i) - { - hotspotStatistics.HotspotPeakIndex[i] = hotspotPeakIndex[i]; - } + sphereImage->TransformIndexToPhysicalPoint(sphereIndex, spherePoint); + inputImage->TransformPhysicalPointToIndex(spherePoint, regionIndex); - chronometer.Report( std::cout ); - return hotspotStatistics; - } + if(region.IsInside(regionIndex) == true) + { + isInside = true; + } + else + { + isInside = false; + return isInside; + } + // Point 13/14 + sphereIndex[2] = sphereCenterIndex[2] + offsetZDirection; + sphereImage->TransformIndexToPhysicalPoint(sphereIndex, spherePoint); + inputImage->TransformPhysicalPointToIndex(spherePoint, regionIndex); - template < typename TPixel, unsigned int VImageDimension > - void ImageStatisticsCalculator::InternalCalculateMaskFromPlanarFigure( - const itk::Image< TPixel, VImageDimension > *image, unsigned int axis ) - { - typedef itk::Image< TPixel, VImageDimension > ImageType; + if(region.IsInside(regionIndex) == true) + { + isInside = true; + } + else + { + isInside = false; + return isInside; + } - typedef itk::CastImageFilter< ImageType, MaskImage2DType > CastFilterType; + addOrSub = true; + } + return isInside; +} - 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::Geometry2D *planarFigureGeometry2D = m_PlanarFigure->GetGeometry2D(); - const typename PlanarFigure::PolyLineType planarFigurePolyline = m_PlanarFigure->GetPolyLine( 0 ); - const mitk::Geometry3D *imageGeometry3D = m_Image->GetGeometry( 0 ); - Determine x- and y-dimensions depending on principal axis - int i0, i1; - switch ( axis ) - { +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::Geometry2D *planarFigureGeometry2D = m_PlanarFigure->GetGeometry2D(); + const typename PlanarFigure::PolyLineType planarFigurePolyline = m_PlanarFigure->GetPolyLine( 0 ); + const mitk::Geometry3D *imageGeometry3D = m_Image->GetGeometry( 0 ); + + // Determine x- and y-dimensions depending on 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 - planarFigureGeometry2D->Map( it->Point, 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 ); + } + 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; - points->InsertNextPoint( point3D[i0], point3D[i1], 0 ); - } + // Convert 2D point back to the local index coordinates of the selected + // image + planarFigureGeometry2D->Map( it->Point, point3D ); - 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))) + // Polygons (partially) outside of the image bounds can not be processed + // further due to a bug in vtkPolyDataToImageStencil + if ( !imageGeometry3D->IsInside( point3D ) ) { - mitkThrow() << "Figure has a zero area and cannot be used for masking."; + outOfBounds = true; } - if ( outOfBounds ) - { - throw std::runtime_error( "Figure at least partially outside of image bounds!" ); - } + imageGeometry3D->WorldToIndex( point3D, point3D ); + + points->InsertNextPoint( point3D[i0], point3D[i1], 0 ); + } - create a vtkLassoStencilSource and set the points of the Polygon - vtkSmartPointer lassoStencil = vtkSmartPointer::New(); - lassoStencil->SetShapeToPolygon(); - lassoStencil->SetPoints( points ); + // 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."; + } - Export from ITK to VTK (to use a VTK filter) - typedef itk::VTKImageImport< MaskImage2DType > ImageImportType; - typedef itk::VTKImageExport< MaskImage2DType > ImageExportType; + if ( outOfBounds ) + { + throw std::runtime_error( "Figure at least partially outside of image bounds!" ); + } - typename ImageExportType::Pointer itkExporter = ImageExportType::New(); - itkExporter->SetInput( castFilter->GetOutput() ); + // create a vtkLassoStencilSource and set the points of the Polygon + vtkSmartPointer lassoStencil = vtkSmartPointer::New(); + lassoStencil->SetShapeToPolygon(); + lassoStencil->SetPoints( points ); - vtkSmartPointer vtkImporter = vtkSmartPointer::New(); - this->ConnectPipelines( itkExporter, vtkImporter ); + // Export from ITK to VTK (to use a VTK filter) + typedef itk::VTKImageImport< MaskImage2DType > ImageImportType; + typedef itk::VTKImageExport< MaskImage2DType > ImageExportType; - Apply the generated image stencil to the input image - vtkSmartPointer imageStencilFilter = vtkSmartPointer::New(); - imageStencilFilter->SetInputConnection( vtkImporter->GetOutputPort() ); - imageStencilFilter->SetStencil( lassoStencil->GetOutput() ); - imageStencilFilter->ReverseStencilOff(); - imageStencilFilter->SetBackgroundValue( 0 ); - imageStencilFilter->Update(); + typename ImageExportType::Pointer itkExporter = ImageExportType::New(); + itkExporter->SetInput( castFilter->GetOutput() ); - Export from VTK back to ITK - vtkSmartPointer vtkExporter = vtkImageExport::New(); // TODO: this is WRONG, should be vtkSmartPointer::New(), but bug # 14455 - vtkExporter->SetInputConnection( imageStencilFilter->GetOutputPort() ); - vtkExporter->Update(); + vtkSmartPointer vtkImporter = vtkSmartPointer::New(); + this->ConnectPipelines( itkExporter, vtkImporter ); - typename ImageImportType::Pointer itkImporter = ImageImportType::New(); - this->ConnectPipelines( vtkExporter, itkImporter ); - itkImporter->Update(); + // Apply the generated image stencil to the input image + vtkSmartPointer imageStencilFilter = vtkSmartPointer::New(); + imageStencilFilter->SetInputConnection( vtkImporter->GetOutputPort() ); + imageStencilFilter->SetStencil( lassoStencil->GetOutput() ); + imageStencilFilter->ReverseStencilOff(); + imageStencilFilter->SetBackgroundValue( 0 ); + imageStencilFilter->Update(); - Store mask - m_InternalImageMask2D = itkImporter->GetOutput(); - } + // Export from VTK back to ITK + vtkSmartPointer vtkExporter = vtkImageExport::New(); // TODO: this is WRONG, should be vtkSmartPointer::New(), but bug # 14455 + vtkExporter->SetInputConnection( imageStencilFilter->GetOutputPort() ); + vtkExporter->Update(); + typename ImageImportType::Pointer itkImporter = ImageImportType::New(); + this->ConnectPipelines( vtkExporter, itkImporter ); + itkImporter->Update(); - 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() ); - } - } + // Store mask + m_InternalImageMask2D = itkImporter->GetOutput(); +} - void ImageStatisticsCalculator::MaskedStatisticsProgressUpdate() +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/mitkImageStatisticsCalculator.h b/Modules/ImageStatistics/mitkImageStatisticsCalculator.h index ccfd992242..976eafa2b1 100644 --- a/Modules/ImageStatistics/mitkImageStatisticsCalculator.h +++ b/Modules/ImageStatistics/mitkImageStatisticsCalculator.h @@ -1,387 +1,403 @@ /*=================================================================== 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. ===================================================================*/ #ifndef _MITK_IMAGESTATISTICSCALCULATOR_H #define _MITK_IMAGESTATISTICSCALCULATOR_H #include #include "ImageStatisticsExports.h" #include #include #ifndef __itkHistogram_h #include #endif #include "mitkImage.h" #include "mitkImageTimeSelector.h" #include "mitkPlanarFigure.h" #include namespace mitk { /** * \brief Class for calculating statistics and histogram for an (optionally * masked) image. * * Images can be masked by either a label image (of the same dimensions as * the original image) or by a closed mitk::PlanarFigure, e.g. a circle or * polygon. When masking with a planar figure, the slice corresponding to the * plane containing the figure is extracted and then clipped with contour * defined by the figure. Planar figures need to be aligned along the main axes * of the image (axial, sagittal, coronal). Planar figures on arbitrary * rotated planes are not supported. * * For each operating mode (no masking, masking by image, masking by planar * figure), the calculated statistics and histogram are cached so that, when * switching back and forth between operation modes without modifying mask or * image, the information doesn't need to be recalculated. * * The class also has the possibility to calculate minimum, maximum, mean * and their corresponding indicies in the hottest spot in a given ROI / VOI. * The size of the hotspot is defined by a sphere with a radius specified by - * the user. + * the user. This procedure is required for the calculation of SUV-statistics + * in PET-images for example. * * Note: currently time-resolved and multi-channel pictures are not properly * supported. */ class ImageStatistics_EXPORT ImageStatisticsCalculator : public itk::Object { public: enum { MASKING_MODE_NONE = 0, MASKING_MODE_IMAGE, MASKING_MODE_PLANARFIGURE }; typedef itk::Statistics::Histogram HistogramType; typedef HistogramType::ConstIterator HistogramConstIteratorType; struct Statistics { int Label; unsigned int N; double Min; double Max; double Mean; double Median; double Variance; double Sigma; double RMS; double HotspotMin; double HotspotMax; double HotspotMean; double HotspotSigma; double HotspotPeak; vnl_vector< int > MinIndex; vnl_vector< int > MaxIndex; vnl_vector HotspotMaxIndex; vnl_vector HotspotMinIndex; vnl_vector HotspotPeakIndex; void Reset() { Label = 0; N = 0; Min = 0.0; Max = 0.0; Mean = 0.0; Median = 0.0; Variance = 0.0; Sigma = 0.0; RMS = 0.0; HotspotMin = 0.0; HotspotMax = 0.0; HotspotMean = 0.0; HotspotPeak = 0.0; HotspotSigma = 0.0; } }; struct MinMaxIndex { double Max; double Min; vnl_vector MaxIndex; vnl_vector MinIndex; }; typedef std::vector< HistogramType::ConstPointer > HistogramContainer; typedef std::vector< Statistics > StatisticsContainer; mitkClassMacro( ImageStatisticsCalculator, itk::Object ); itkNewMacro( ImageStatisticsCalculator ); /** \brief Set image from which to compute statistics. */ void SetImage( const mitk::Image *image ); /** \brief Set image for masking. */ void SetImageMask( const mitk::Image *imageMask ); /** \brief Set planar figure for masking. */ void SetPlanarFigure( mitk::PlanarFigure *planarFigure ); /** \brief Set/Get operation mode for masking */ void SetMaskingMode( unsigned int mode ); /** \brief Set/Get operation mode for masking */ itkGetMacro( MaskingMode, unsigned int ); /** \brief Set/Get operation mode for masking */ void SetMaskingModeToNone(); /** \brief Set/Get operation mode for masking */ void SetMaskingModeToImage(); /** \brief Set/Get operation mode for masking */ void SetMaskingModeToPlanarFigure(); /** \brief Set a pixel value for pixels that will be ignored in the statistics */ void SetIgnorePixelValue(double value); /** \brief Get the pixel value for pixels that will be ignored in the statistics */ double GetIgnorePixelValue(); /** \brief Set whether a pixel value should be ignored in the statistics */ void SetDoIgnorePixelValue(bool doit); /** \brief Get whether a pixel value will be ignored in the statistics */ bool GetDoIgnorePixelValue(); /** \brief Sets the radius for the hotspot */ void SetHotspotRadius (double hotspotRadiusInMM); /** \brief Returns the radius of the hotspot */ double GetHotspotRadius(); /** \brief Sets whether the hotspot should be calculated */ void SetCalculateHotspot(bool calculateHotspot); /** \brief Returns true whether the hotspot should be calculated, otherwise false */ bool IsHotspotCalculated(); /** \brief Compute statistics (together with histogram) for the current * masking mode. * * Computation is not executed if statistics is already up to date. In this * case, false is returned; otherwise, true.*/ virtual bool ComputeStatistics( unsigned int timeStep = 0 ); /** \brief Retrieve the histogram depending on the current masking mode. * * \param label The label for which to retrieve the histogram in multi-label situations (ascending order). */ const HistogramType *GetHistogram( unsigned int timeStep = 0, unsigned int label = 0 ) const; /** \brief Retrieve the histogram depending on the current masking mode (for all image labels. */ const HistogramContainer &GetHistogramVector( unsigned int timeStep = 0 ) const; /** \brief Retrieve statistics depending on the current masking mode. * * \param label The label for which to retrieve the statistics in multi-label situations (ascending order). */ const Statistics &GetStatistics( unsigned int timeStep = 0, unsigned int label = 0 ) const; /** \brief Retrieve statistics depending on the current masking mode (for all image labels). */ const StatisticsContainer &GetStatisticsVector( unsigned int timeStep = 0 ) const; protected: typedef std::vector< HistogramContainer > HistogramVector; typedef std::vector< StatisticsContainer > StatisticsVector; typedef std::vector< itk::TimeStamp > TimeStampVectorType; typedef std::vector< bool > BoolVectorType; typedef itk::Image< unsigned short, 3 > MaskImage3DType; typedef itk::Image< unsigned short, 2 > MaskImage2DType; ImageStatisticsCalculator(); virtual ~ImageStatisticsCalculator(); /** \brief Depending on the masking mode, the image and mask from which to * calculate statistics is extracted from the original input image and mask * data. * * For example, a when using a PlanarFigure as mask, the 2D image slice * corresponding to the PlanarFigure will be extracted from the original * image. If masking is disabled, the original image is simply passed * through. */ void ExtractImageAndMask( unsigned int timeStep = 0 ); /** \brief If the passed vector matches any of the three principal axes * of the passed geometry, the ínteger value corresponding to the axis * is set and true is returned. */ bool GetPrincipalAxis( const Geometry3D *geometry, Vector3D vector, unsigned int &axis ); template < typename TPixel, unsigned int VImageDimension > void InternalCalculateStatisticsUnmasked( const itk::Image< TPixel, VImageDimension > *image, StatisticsContainer* statisticsContainer, HistogramContainer *histogramContainer ); template < typename TPixel, unsigned int VImageDimension > void InternalCalculateStatisticsMasked( const itk::Image< TPixel, VImageDimension > *image, itk::Image< unsigned short, VImageDimension > *maskImage, StatisticsContainer* statisticsContainer, HistogramContainer* histogramContainer ); template < typename TPixel, unsigned int VImageDimension > void InternalCalculateMaskFromPlanarFigure( const itk::Image< TPixel, VImageDimension > *image, unsigned int axis ); template < typename TPixel, unsigned int VImageDimension > void InternalMaskIgnoredPixels( const itk::Image< TPixel, VImageDimension > *image, itk::Image< unsigned short, VImageDimension > *maskImage ); /** \brief Calculates minimum, maximum, mean value and their * corresponding indices in a given ROI. As input the function - * needs an image and a mask It returns a MinMaxIndex object. */ + * needs an image, a mask and for limiting the searched value + * minimum- and maximum boundaries. The boundaries are set + * to the extent of float by default. It returns a + * MinMaxIndex object. */ template MinMaxIndex CalculateMinMaxIndex( const itk::Image *inputImage, - itk::Image *maskImage); + itk::Image *maskImage, + float minBoundary, + float maxBoundary); /** \brief Calculates the hotspot statistics within a given * ROI. As input the function needs an image, a mask which * represents the ROI and a radius which defines the size of * the sphere. The function returns a Statistics object. */ template < typename TPixel, unsigned int VImageDimension> Statistics CalculateHotspotStatistics( const itk::Image *inputImage, itk::Image *maskImage, double radiusInMM); + /** \brief Calculates if a sphere is located within an image. + * As input the function needs an image for whose region should + * be check if the sphere is inside and an image of the sphere. + * To ensure speed only 14 points on the surface of the sphere + * are checked: six on the axis of coordinates and eight more. + * The function returns true if every point is in the region. */ + template < typename TPixel, unsigned int VImageDimension> + bool IsSphereInsideRegion( + const itk::Image *inputImage, + itk::Image *sphereImage); /** 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()); } void UnmaskedStatisticsProgressUpdate(); void MaskedStatisticsProgressUpdate(); /** m_Image contains the input image (e.g. 2D, 3D, 3D+t)*/ mitk::Image::ConstPointer m_Image; mitk::Image::ConstPointer m_ImageMask; mitk::PlanarFigure::Pointer m_PlanarFigure; HistogramVector m_ImageHistogramVector; HistogramVector m_MaskedImageHistogramVector; HistogramVector m_PlanarFigureHistogramVector; HistogramType::Pointer m_EmptyHistogram; HistogramContainer m_EmptyHistogramContainer; StatisticsVector m_ImageStatisticsVector; StatisticsVector m_MaskedImageStatisticsVector; StatisticsVector m_PlanarFigureStatisticsVector; StatisticsVector m_MaskedImageHotspotStatisticsVector; Statistics m_EmptyStatistics; StatisticsContainer m_EmptyStatisticsContainer; unsigned int m_MaskingMode; bool m_MaskingModeChanged; /** m_InternalImage contains a image volume at one time step (e.g. 2D, 3D)*/ mitk::Image::ConstPointer m_InternalImage; MaskImage3DType::Pointer m_InternalImageMask3D; MaskImage2DType::Pointer m_InternalImageMask2D; TimeStampVectorType m_ImageStatisticsTimeStampVector; TimeStampVectorType m_MaskedImageStatisticsTimeStampVector; TimeStampVectorType m_PlanarFigureStatisticsTimeStampVector; BoolVectorType m_ImageStatisticsCalculationTriggerVector; BoolVectorType m_MaskedImageStatisticsCalculationTriggerVector; BoolVectorType m_PlanarFigureStatisticsCalculationTriggerVector; double m_IgnorePixelValue; bool m_DoIgnorePixelValue; bool m_IgnorePixelValueChanged; double m_HotspotRadius; bool m_CalculateHotspot; unsigned int m_PlanarFigureAxis; // Normal axis for PlanarFigure unsigned int m_PlanarFigureSlice; // Slice which contains PlanarFigure int m_PlanarFigureCoordinate0; // First plane-axis for PlanarFigure int m_PlanarFigureCoordinate1; // Second plane-axis for PlanarFigure }; } #endif // #define _MITK_IMAGESTATISTICSCALCULATOR_H