diff --git a/Modules/ImageStatistics/Testing/mitkImageStatisticsCalculatorTest.cpp b/Modules/ImageStatistics/Testing/mitkImageStatisticsCalculatorTest.cpp index db48266dc9..0e3c62345e 100644 --- a/Modules/ImageStatistics/Testing/mitkImageStatisticsCalculatorTest.cpp +++ b/Modules/ImageStatistics/Testing/mitkImageStatisticsCalculatorTest.cpp @@ -1,623 +1,496 @@ /*=================================================================== 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 "mitkStandardFileLocations.h" #include "mitkDicomSeriesReader.h" #include "mitkTestingMacros.h" #include "mitkImageStatisticsCalculator.h" #include "mitkPlanarPolygon.h" #include "mitkDicomSeriesReader.h" #include #include "vtkStreamingDemandDrivenPipeline.h" #include #include -#include "itkImage.h" -#include -#include -#include - //#include /** * \brief Test class for mitkImageStatisticsCalculator * * This test covers: * - instantiation of an ImageStatisticsCalculator class * - correctness of statistics when using PlanarFigures for masking */ class mitkImageStatisticsCalculatorTestClass { public: struct testCase { int id; mitk::PlanarFigure::Pointer figure; double mean; double sd; }; // calculate statistics for the given image and planarpolygon static const mitk::ImageStatisticsCalculator::Statistics TestStatistics( mitk::Image::Pointer image, mitk::PlanarFigure::Pointer polygon ) { mitk::ImageStatisticsCalculator::Pointer statisticsCalculator = mitk::ImageStatisticsCalculator::New(); statisticsCalculator->SetImage( image ); statisticsCalculator->SetMaskingModeToPlanarFigure(); statisticsCalculator->SetPlanarFigure( polygon ); statisticsCalculator->ComputeStatistics(); return statisticsCalculator->GetStatistics(); } - // calculate statistics for the given image and mask - static const mitk::ImageStatisticsCalculator::Statistics TestStatisticsWithMask(mitk::Image::Pointer image, mitk::Image::Pointer mask) - { - mitk::ImageStatisticsCalculator::Pointer statisticsCalculator = mitk::ImageStatisticsCalculator::New(); - statisticsCalculator->SetImage(image); - statisticsCalculator->SetMaskingModeToImage(); - statisticsCalculator->SetImageMask(mask); - statisticsCalculator->ComputeStatistics(); - return statisticsCalculator->GetStatistics(); - } // returns a vector of defined test-cases static std::vector InitializeTestCases( mitk::Geometry2D::Pointer geom ) { std::vector testCases; { /***************************** * one whole white pixel * -> mean of 255 expected ******************************/ mitk::PlanarPolygon::Pointer figure1 = mitk::PlanarPolygon::New(); figure1->SetGeometry2D( geom ); mitk::Point2D pnt1; pnt1[0] = 10.5 ; pnt1[1] = 3.5; figure1->PlaceFigure( pnt1 ); mitk::Point2D pnt2; pnt2[0] = 9.5; pnt2[1] = 3.5; figure1->SetControlPoint( 1, pnt2, true ); mitk::Point2D pnt3; pnt3[0] = 9.5; pnt3[1] = 4.5; figure1->SetControlPoint( 2, pnt3, true ); mitk::Point2D pnt4; pnt4[0] = 10.5; pnt4[1] = 4.5; figure1->SetControlPoint( 3, pnt4, true ); figure1->GetPolyLine(0); testCase test; test.id = testCases.size(); test.figure = figure1; test.mean = 255.0; test.sd = 0.0; testCases.push_back( test ); } { /***************************** * half pixel in x-direction (white) * -> mean of 255 expected ******************************/ mitk::PlanarPolygon::Pointer figure1 = mitk::PlanarPolygon::New(); figure1->SetGeometry2D( geom ); mitk::Point2D pnt1; pnt1[0] = 10.0 ; pnt1[1] = 3.5; figure1->PlaceFigure( pnt1 ); mitk::Point2D pnt2; pnt2[0] = 9.5; pnt2[1] = 3.5; figure1->SetControlPoint( 1, pnt2, true ); mitk::Point2D pnt3; pnt3[0] = 9.5; pnt3[1] = 4.5; figure1->SetControlPoint( 2, pnt3, true ); mitk::Point2D pnt4; pnt4[0] = 10.0; pnt4[1] = 4.5; figure1->SetControlPoint( 3, pnt4, true ); figure1->GetPolyLine(0); testCase test; test.id = testCases.size(); test.figure = figure1; test.mean = 255.0; test.sd = 0.0; testCases.push_back( test ); } { /***************************** * half pixel in diagonal-direction (white) * -> mean of 255 expected ******************************/ mitk::PlanarPolygon::Pointer figure1 = mitk::PlanarPolygon::New(); figure1->SetGeometry2D( geom ); mitk::Point2D pnt1; pnt1[0] = 10.5 ; pnt1[1] = 3.5; figure1->PlaceFigure( pnt1 ); mitk::Point2D pnt2; pnt2[0] = 9.5; pnt2[1] = 3.5; figure1->SetControlPoint( 1, pnt2, true ); mitk::Point2D pnt3; pnt3[0] = 9.5; pnt3[1] = 4.5; figure1->SetControlPoint( 2, pnt3, true ); figure1->GetPolyLine(0); testCase test; test.id = testCases.size(); test.figure = figure1; test.mean = 255.0; test.sd = 0.0; testCases.push_back( test ); } { /***************************** * one pixel (white) + 2 half pixels (white) + 1 half pixel (black) * -> mean of 191.25 expected ******************************/ mitk::PlanarPolygon::Pointer figure1 = mitk::PlanarPolygon::New(); figure1->SetGeometry2D( geom ); mitk::Point2D pnt1; pnt1[0] = 1.1; pnt1[1] = 1.1; figure1->PlaceFigure( pnt1 ); mitk::Point2D pnt2; pnt2[0] = 2.0; pnt2[1] = 2.0; figure1->SetControlPoint( 1, pnt2, true ); mitk::Point2D pnt3; pnt3[0] = 3.0; pnt3[1] = 1.0; figure1->SetControlPoint( 2, pnt3, true ); mitk::Point2D pnt4; pnt4[0] = 2.0; pnt4[1] = 0.0; figure1->SetControlPoint( 3, pnt4, true ); figure1->GetPolyLine(0); testCase test; test.id = testCases.size(); test.figure = figure1; test.mean = 191.25; test.sd = 127.5; testCases.push_back( test ); } { /***************************** * whole pixel (white) + half pixel (gray) in x-direction * -> mean of 191.5 expected ******************************/ mitk::PlanarPolygon::Pointer figure1 = mitk::PlanarPolygon::New(); figure1->SetGeometry2D( geom ); mitk::Point2D pnt1; pnt1[0] = 11.0; pnt1[1] = 3.5; figure1->PlaceFigure( pnt1 ); mitk::Point2D pnt2; pnt2[0] = 9.5; pnt2[1] = 3.5; figure1->SetControlPoint( 1, pnt2, true ); mitk::Point2D pnt3; pnt3[0] = 9.5; pnt3[1] = 4.5; figure1->SetControlPoint( 2, pnt3, true ); mitk::Point2D pnt4; pnt4[0] = 11.0; pnt4[1] = 4.5; figure1->SetControlPoint( 3, pnt4, true ); figure1->GetPolyLine(0); testCase test; test.id = testCases.size(); test.figure = figure1; test.mean = 191.50; test.sd = 89.80; testCases.push_back( test ); } { /***************************** * quarter pixel (black) + whole pixel (white) + half pixel (gray) in x-direction * -> mean of 191.5 expected ******************************/ mitk::PlanarPolygon::Pointer figure1 = mitk::PlanarPolygon::New(); figure1->SetGeometry2D( geom ); mitk::Point2D pnt1; pnt1[0] = 11.0; pnt1[1] = 3.5; figure1->PlaceFigure( pnt1 ); mitk::Point2D pnt2; pnt2[0] = 9.25; pnt2[1] = 3.5; figure1->SetControlPoint( 1, pnt2, true ); mitk::Point2D pnt3; pnt3[0] = 9.25; pnt3[1] = 4.5; figure1->SetControlPoint( 2, pnt3, true ); mitk::Point2D pnt4; pnt4[0] = 11.0; pnt4[1] = 4.5; figure1->SetControlPoint( 3, pnt4, true ); figure1->GetPolyLine(0); testCase test; test.id = testCases.size(); test.figure = figure1; test.mean = 191.5; test.sd = 89.80; testCases.push_back( test ); } { /***************************** * half pixel (black) + whole pixel (white) + half pixel (gray) in x-direction * -> mean of 127.66 expected ******************************/ mitk::PlanarPolygon::Pointer figure1 = mitk::PlanarPolygon::New(); figure1->SetGeometry2D( geom ); mitk::Point2D pnt1; pnt1[0] = 11.0; pnt1[1] = 3.5; figure1->PlaceFigure( pnt1 ); mitk::Point2D pnt2; pnt2[0] = 9.0; pnt2[1] = 3.5; figure1->SetControlPoint( 1, pnt2, true ); mitk::Point2D pnt3; pnt3[0] = 9.0; pnt3[1] = 4.0; figure1->SetControlPoint( 2, pnt3, true ); mitk::Point2D pnt4; pnt4[0] = 11.0; pnt4[1] = 4.0; figure1->SetControlPoint( 3, pnt4, true ); figure1->GetPolyLine(0); testCase test; test.id = testCases.size(); test.figure = figure1; test.mean = 127.66; test.sd = 127.5; testCases.push_back( test ); } { /***************************** * whole pixel (gray) * -> mean of 128 expected ******************************/ mitk::PlanarPolygon::Pointer figure2 = mitk::PlanarPolygon::New(); figure2->SetGeometry2D( geom ); mitk::Point2D pnt1; pnt1[0] = 11.5; pnt1[1] = 10.5; figure2->PlaceFigure( pnt1 ); mitk::Point2D pnt2; pnt2[0] = 11.5; pnt2[1] = 11.5; figure2->SetControlPoint( 1, pnt2, true ); mitk::Point2D pnt3; pnt3[0] = 12.5; pnt3[1] = 11.5; figure2->SetControlPoint( 2, pnt3, true ); mitk::Point2D pnt4; pnt4[0] = 12.5; pnt4[1] = 10.5; figure2->SetControlPoint( 3, pnt4, true ); figure2->GetPolyLine(0); testCase test; test.id = testCases.size(); test.figure = figure2; test.mean = 128.0; test.sd = 0.0; testCases.push_back( test ); } { /***************************** * whole pixel (gray) + half pixel (white) in y-direction * -> mean of 191.5 expected ******************************/ mitk::PlanarPolygon::Pointer figure2 = mitk::PlanarPolygon::New(); figure2->SetGeometry2D( geom ); mitk::Point2D pnt1; pnt1[0] = 11.5; pnt1[1] = 10.5; figure2->PlaceFigure( pnt1 ); mitk::Point2D pnt2; pnt2[0] = 11.5; pnt2[1] = 12.0; figure2->SetControlPoint( 1, pnt2, true ); mitk::Point2D pnt3; pnt3[0] = 12.5; pnt3[1] = 12.0; figure2->SetControlPoint( 2, pnt3, true ); mitk::Point2D pnt4; pnt4[0] = 12.5; pnt4[1] = 10.5; figure2->SetControlPoint( 3, pnt4, true ); figure2->GetPolyLine(0); testCase test; test.id = testCases.size(); test.figure = figure2; test.mean = 191.5; test.sd = 89.80; testCases.push_back( test ); } { /***************************** * 2 whole pixel (white) + 2 whole pixel (black) in y-direction * -> mean of 127.66 expected ******************************/ mitk::PlanarPolygon::Pointer figure2 = mitk::PlanarPolygon::New(); figure2->SetGeometry2D( geom ); mitk::Point2D pnt1; pnt1[0] = 11.5; pnt1[1] = 10.5; figure2->PlaceFigure( pnt1 ); mitk::Point2D pnt2; pnt2[0] = 11.5; pnt2[1] = 13.5; figure2->SetControlPoint( 1, pnt2, true ); mitk::Point2D pnt3; pnt3[0] = 12.5; pnt3[1] = 13.5; figure2->SetControlPoint( 2, pnt3, true ); mitk::Point2D pnt4; pnt4[0] = 12.5; pnt4[1] = 10.5; figure2->SetControlPoint( 3, pnt4, true ); figure2->GetPolyLine(0); testCase test; test.id = testCases.size(); test.figure = figure2; test.mean = 127.66; test.sd = 127.5; testCases.push_back( test ); } { /***************************** * 9 whole pixels (white) + 3 half pixels (white) * + 3 whole pixel (black) [ + 3 slightly less than half pixels (black)] * -> mean of 204.0 expected ******************************/ mitk::PlanarPolygon::Pointer figure2 = mitk::PlanarPolygon::New(); figure2->SetGeometry2D( geom ); mitk::Point2D pnt1; pnt1[0] = 0.5; pnt1[1] = 0.5; figure2->PlaceFigure( pnt1 ); mitk::Point2D pnt2; pnt2[0] = 3.5; pnt2[1] = 3.5; figure2->SetControlPoint( 1, pnt2, true ); mitk::Point2D pnt3; pnt3[0] = 8.4999; pnt3[1] = 3.5; figure2->SetControlPoint( 2, pnt3, true ); mitk::Point2D pnt4; pnt4[0] = 5.4999; pnt4[1] = 0.5; figure2->SetControlPoint( 3, pnt4, true ); figure2->GetPolyLine(0); testCase test; test.id = testCases.size(); test.figure = figure2; test.mean = 204.0; test.sd = 105.58; testCases.push_back( test ); } { /***************************** * half pixel (white) + whole pixel (white) + half pixel (black) * -> mean of 212.66 expected ******************************/ mitk::PlanarPolygon::Pointer figure2 = mitk::PlanarPolygon::New(); figure2->SetGeometry2D( geom ); mitk::Point2D pnt1; pnt1[0] = 9.5; pnt1[1] = 0.5; figure2->PlaceFigure( pnt1 ); mitk::Point2D pnt2; pnt2[0] = 9.5; pnt2[1] = 2.5; figure2->SetControlPoint( 1, pnt2, true ); mitk::Point2D pnt3; pnt3[0] = 11.5; pnt3[1] = 2.5; figure2->SetControlPoint( 2, pnt3, true ); figure2->GetPolyLine(0); testCase test; test.id = testCases.size(); test.figure = figure2; test.mean = 212.66; test.sd = 73.32; testCases.push_back( test ); } return testCases; } - - // loads the test image - static mitk::Image::Pointer GetTestImage(const char *nameOfFile, const char *path) + static mitk::Image::Pointer GetTestImage() { mitk::StandardFileLocations::Pointer locator = mitk::StandardFileLocations::GetInstance(); - const std::string filename = locator->FindFile(nameOfFile, path); + const std::string filename = locator->FindFile("testimage.dcm", "Modules/MitkExt/Testing/Data"); if (filename.empty()) { MITK_ERROR << "Could not find test file"; return NULL; } else { - MITK_INFO << "Found " << nameOfFile; + MITK_INFO << "Found testimage.dcm"; } itk::FilenamesContainer file; file.push_back( filename ); mitk::DicomSeriesReader* reader = new mitk::DicomSeriesReader; mitk::DataNode::Pointer node = reader->LoadDicomSeries( file, false, false ); mitk::Image::Pointer image = dynamic_cast( node->GetData() ); return image; } }; int TestUnitilizedImage() { /***************************** * loading uninitialized image to datastorage ******************************/ std::cout << "Testing loading uninitialized image to datastorage:"; try{ MITK_TEST_FOR_EXCEPTION_BEGIN(mitk::Exception) mitk::Image::Pointer image = mitk::Image::New(); mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData(image); mitk::ImageStatisticsCalculator::Pointer is = mitk::ImageStatisticsCalculator::New(); is->ComputeStatistics(); MITK_TEST_FOR_EXCEPTION_END(mitk::Exception) } catch (const mitk::Exception&) { std::cout << "Success: Loading uninitialized image has thrown an exception." << std::endl; } return 0; } -// static method TestHotspotSearch(image): - // -> formuliert verschiedene TestFälle: Polygon und Masken (Segmentierung) - static void TestHotspotSearch(mitk::Image::Pointer image, unsigned int testCase) - { - // testCase for PlanarFigure-1 - if(testCase == 0) - { - - } - - // testCase for Mask-1 - if(testCase == 1) - { - - } - } - - static void CreateTestImage() - { - - typedef itk::Image ThreeDimensionalImageType; - ThreeDimensionalImageType::Pointer image = ThreeDimensionalImageType::New(); - - ThreeDimensionalImageType::SizeType size; - ThreeDimensionalImageType::SpacingType spacing; - ThreeDimensionalImageType::IndexType start; - - for(int i = 0; i < 3; ++i) - { - size[i] = 10; - 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(10); - } - - for(unsigned int x = 2; x <= 8; ++x) - { - for(unsigned int y = 2; y <= 8; ++y) - { - for(unsigned int z = 2; z <= 8; ++z) - { - ThreeDimensionalImageType::IndexType pixelIndex; - pixelIndex[0] = x; - pixelIndex[1] = y; - pixelIndex[2] = z; - - image->SetPixel(pixelIndex, 200); - } - } - } - - ThreeDimensionalImageType::IndexType pixelIndex; - pixelIndex[0] = 5; - pixelIndex[1] = 5; - pixelIndex[2] = 5; - - image->SetPixel(pixelIndex, 255); - - 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] = 10; - 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); - } - - /*typedef itk::ImageFileWriter WriterType; - WriterType::Pointer writer = WriterType::New(); - writer->SetInput(ThreeDimage); - writer->SetFileName("D:\\dreiD.dcm"); - writer->Update();*/ - } - int mitkImageStatisticsCalculatorTest(int, char* []) { // always start with this! MITK_TEST_BEGIN("mitkImageStatisticsCalculatorTest") //QCoreApplication app(argc, argv); - mitk::Image::Pointer image = mitkImageStatisticsCalculatorTestClass::GetTestImage(' ', ' '); + mitk::Image::Pointer image = mitkImageStatisticsCalculatorTestClass::GetTestImage(); MITK_TEST_CONDITION_REQUIRED( image.IsNotNull(), "Loading test image" ); mitk::Geometry2D::Pointer geom = image->GetSlicedGeometry()->GetGeometry2D(0); std::vector allTestCases = mitkImageStatisticsCalculatorTestClass::InitializeTestCases( geom ); for ( std::vector::size_type i=0; i + #include + #include #include struct mitkImageStatisticsHotspotTestClass { /** \brief Test parameters for one test case. Describes all aspects of a single test case: - parameters to generate a test image - parameters of a ROI that describes where to calculate statistics - expected statistics results */ struct Parameters { - int m_ExpectedValueU; - int m_ExpectedValueV; - - int m_ParameterX; - int m_ParameterY; + // XML-Tag + int m_ImageRows; // XML-Tag "image-rows" + int m_ImageColumns; // XML-Tag "image-columns" + int m_ImageSlices; // XML-Tag "image-slices" + int m_NumberOfGaussian; // XML-Tag "numberOfGaussians" + + // XML-Tag + std::vector m_CenterX; // XML-Tag "centerIndexX" + std::vector m_CenterY; // XML-Tag "centerIndexY" + std::vector m_CenterZ; // XML-Tag "centerIndexZ" + + std::vector m_SigmaX; // XML-Tag "deviationX" + std::vector m_SigmaY; // XML-Tag "deviationY" + std::vector m_SigmaZ; // XML-Tag "deviationZ" + + std::vector m_Altitude; // XML-Tag "altitude" + + // XML-Tag + float m_HotspotMinimum; // XML-Tag "minimum" + float m_HotspotMaximum; // XML-Tag "maximum" + float m_HotspotPeak; // XML-Tag "peak" + + std::vector m_HotspotMaximumIndex; // XML-Tag "maximumIndexX", XML-Tag "maximumIndexY", XML-Tag "maximumIndexZ" + std::vector m_HotspotMinimumIndex; // XML-Tag "minimumIndexX", XML-Tag "minimumIndexY", XML-Tag "minimumIndexZ" + std::vector m_HotspotPeakIndex; // XML-Tag "peakIndexX", XML-Tag "peakIndexY", XML-Tag "peakIndexZ" }; /** \brief Find/Convert integer attribute in itk::DOMNode. */ static int GetIntegerAttribute(itk::DOMNode* domNode, const std::string& tag) { assert(domNode); MITK_TEST_CONDITION_REQUIRED( domNode->HasAttribute(tag), "Tag '" << tag << "' is defined in test parameters" ); std::string attributeValue = domNode->GetAttribute(tag); int resultValue; try { //MITK_TEST_OUTPUT( << "Converting tag value '" << attributeValue << "' for tag '" << tag << "' to integer"); std::stringstream(attributeValue) >> resultValue; return resultValue; } catch(std::exception& e) { MITK_TEST_CONDITION_REQUIRED(false, "Convert tag value '" << attributeValue << "' for tag '" << tag << "' to integer"); return 0; // just to satisfy compiler } } - /** \brief Find/Convert double attribute in itk::DOMNode. */ static double GetDoubleAttribute(itk::DOMNode* domNode, const std::string& tag) { assert(domNode); MITK_TEST_CONDITION_REQUIRED( domNode->HasAttribute(tag), "Tag '" << tag << "' is defined in test parameters" ); std::string attributeValue = domNode->GetAttribute(tag); double resultValue; try { //MITK_TEST_OUTPUT( << "Converting tag value '" << attributeValue << "' for tag '" << tag << "' to double"); std::stringstream(attributeValue) >> resultValue; return resultValue; } catch(std::exception& e) { MITK_TEST_CONDITION_REQUIRED(false, "Convert tag value '" << attributeValue << "' for tag '" << tag << "' to double"); return 0.0; // just to satisfy compiler } } /** \brief Read XML file describging the test parameters. Reads XML file given in first commandline parameter in order to construct a Parameters structure. The XML file should be structurs as the following example, i.e. we describe the three test aspects of Parameters in three different tags, with all the details described as tag attributes. \verbatim \verbatim The different parameters are interpreted as follows: ... TODO TODO TODO ... */ static Parameters ParseParameters(int argc, char* argv[]) { // - parse parameters // - fill ALL values of result structure // - if necessary, provide c'tor and default parameters to Parameters MITK_TEST_CONDITION_REQUIRED(argc == 2, "Test is invoked with exactly 1 parameter (XML parameters file)"); MITK_INFO << "Reading parameters from file '" << argv[1] << "'"; std::string filename = argv[1]; Parameters result; itk::DOMNodeXMLReader::Pointer xmlReader = itk::DOMNodeXMLReader::New(); xmlReader->SetFileName( filename ); try { xmlReader->Update(); itk::DOMNode::Pointer domRoot = xmlReader->GetOutput(); - + typedef std::vector NodeList; // read test image parameters, fill result structure - std::vector testimages; + NodeList testimages; domRoot->GetChildren("testimage", testimages); MITK_TEST_CONDITION_REQUIRED( testimages.size() == 1, "One test image defined" ) itk::DOMNode* testimage = testimages[0]; - int x = GetIntegerAttribute( testimage, "x" ); - int y = GetIntegerAttribute( testimage, "y" ); + result.m_ImageRows = GetIntegerAttribute( testimage, "image-rows" ); + result.m_ImageColumns = GetIntegerAttribute( testimage, "image-columns" ); + result.m_ImageSlices = GetIntegerAttribute( testimage, "image-slices" ); + + result.m_NumberOfGaussian = GetIntegerAttribute( testimage, "numberOfGaussians" ); + + MITK_TEST_OUTPUT( << "Read size parameters (x,y,z): " << result.m_ImageRows << "," << result.m_ImageColumns << "," << result.m_ImageSlices); + + NodeList gaussians; + testimage->GetChildren("gaussian", gaussians); + MITK_TEST_CONDITION_REQUIRED( gaussians.size() >= 1, "At least one gaussian is defined" ) - MITK_TEST_OUTPUT( << "Read integer parameters (x,y)=" << x << "," << y); + std::vector tmpCenterX(result.m_NumberOfGaussian,0); + std::vector tmpCenterY(result.m_NumberOfGaussian,0); + std::vector tmpCenterZ(result.m_NumberOfGaussian,0); - double xd = GetDoubleAttribute( testimage, "x" ); - double yd = GetDoubleAttribute( testimage, "y" ); + std::vector tmpSigmaX(result.m_NumberOfGaussian,0); + std::vector tmpSigmaY(result.m_NumberOfGaussian,0); + std::vector tmpSigmaZ(result.m_NumberOfGaussian,0); - MITK_TEST_OUTPUT( << "Read double parameters (x,y)=" << xd << "," << yd); + std::vector tmpAltitude(result.m_NumberOfGaussian,0); + for(int i = 0; i < result.m_NumberOfGaussian ; ++i) + { + itk::DOMNode* gaussian = gaussians[i]; + + tmpCenterX[i] = GetIntegerAttribute(gaussian, "centerIndexX"); + tmpCenterY[i] = GetIntegerAttribute(gaussian, "centerIndexY"); + tmpCenterZ[i] = GetIntegerAttribute(gaussian, "centerIndexZ"); + + tmpSigmaX[i] = GetIntegerAttribute(gaussian, "deviationX"); + tmpSigmaY[i] = GetIntegerAttribute(gaussian, "deviationY"); + tmpSigmaZ[i] = GetIntegerAttribute(gaussian, "deviationZ"); + + tmpAltitude[i] = GetIntegerAttribute(gaussian, "altitude"); + } + + result.m_CenterX = tmpCenterX; + result.m_CenterY = tmpCenterY; + result.m_CenterZ = tmpCenterZ; + + result.m_SigmaX = tmpSigmaX; + result.m_SigmaY = tmpSigmaY; + result.m_SigmaZ = tmpSigmaZ; + + result.m_Altitude = tmpAltitude; // read ROI parameters, fill result structure - std::vector rois; + NodeList rois; domRoot->GetChildren("roi", rois); MITK_TEST_CONDITION_REQUIRED( rois.size() == 1, "One ROI defined" ) itk::DOMNode* roi = rois[0]; - double roiU = GetDoubleAttribute( roi, "u" ); - double roiV = GetDoubleAttribute( roi, "v" ); + result.m_RoiMaximumX = GetIntegerAttribute(roi, "maximumX"); + result.m_RoiMinimumX = GetIntegerAttribute(roi, "minimumX"); + result.m_RoiMaximumY = GetIntegerAttribute(roi, "maximumY"); + result.m_RoiMinimumY = GetIntegerAttribute(roi, "minimumY"); + result.m_RoiMaximumZ = GetIntegerAttribute(roi, "maximumZ"); + result.m_RoiMinimumZ = GetIntegerAttribute(roi, "minimumZ"); + + // read statistic parameters, fill result structure + NodeList statistics; + domRoot->GetChildren("statistic", statistics); + MITK_TEST_CONDITION_REQUIRED( statistics.size() == 1, "One statistic defined" ) + itk::DOMNode* statistic = statistics[0]; + + result.m_HotspotMinimum = GetDoubleAttribute(statistic, "minimum"); + result.m_HotspotMaximum = GetDoubleAttribute(statistic, "maximum"); + result.m_HotspotPeak = GetDoubleAttribute(statistic, "peak"); + + std::vector tmpMinimumIndex(3,0); - double roiUndefined = GetDoubleAttribute( roi, "imWagenVorMir" ); + tmpMinimumIndex[0] = GetIntegerAttribute(statistic, "minimumIndexX"); + tmpMinimumIndex[1] = GetIntegerAttribute(statistic, "minimumIndexY"); + tmpMinimumIndex[2] = GetIntegerAttribute(statistic, "minimumIndexZ"); + + result.m_HotspotMinimumIndex = tmpMinimumIndex; + + + std::vector tmpMaximumIndex(3,0); + + tmpMaximumIndex[0] = GetIntegerAttribute(statistic, "maximumIndexX"); + tmpMaximumIndex[1] = GetIntegerAttribute(statistic, "maximumIndexY"); + tmpMaximumIndex[2] = GetIntegerAttribute(statistic, "maximumIndexZ"); + + result.m_HotspotMaximumIndex = tmpMaximumIndex; + + std::vector tmpPeakIndex(3,0); + + tmpPeakIndex[0] = GetIntegerAttribute(statistic, "peakIndexX"); + tmpPeakIndex[1] = GetIntegerAttribute(statistic, "peakIndexY"); + tmpPeakIndex[2] = GetIntegerAttribute(statistic, "peakIndexZ"); + + result.m_HotspotPeakIndex = tmpPeakIndex; + + return result; } catch (std::exception& e) { MITK_TEST_CONDITION_REQUIRED(false, "Reading test parameters from XML file. Error message: " << e.what()); } if (false /* and all parameters nicely found */) { return result; } else { throw std::invalid_argument("Test called with invalid parameters.."); // TODO provide details if possible } } /** \brief Generate an image that contains a couple of 3D gaussian distributions. Uses the given parameters to produce a test image using class TODO... bla */ + static mitk::Image::Pointer BuildTestImage(const Parameters& testParameters) { + // evaluate parameters, create corresponding image mitk::Image::Pointer result; - // evaluate parameters, create corresponding image + typedef double PixelType; + const unsigned int Dimension = 3; + typedef itk::Image ImageType; + ImageType::Pointer image = ImageType::New(); + typedef itk::MultiGaussianImageSource< ImageType > MultiGaussianImageSource; + MultiGaussianImageSource::Pointer gaussianGenerator = MultiGaussianImageSource::New(); + ImageType::SizeValueType size[3]; + size[0] = testParameters.m_ImageColumns; + size[1] = testParameters.m_ImageRows; + size[2] = testParameters.m_ImageSlices; + + itk::MultiGaussianImageSource::VectorType centerXVec, centerYVec, centerZVec, sigmaXVec, sigmaYVec, sigmaZVec, altitudeVec; + + for(int i = 0; i < testParameters.m_NumberOfGaussian; ++i) + { + centerXVec.push_back(testParameters.m_CenterX[i]); + centerYVec.push_back(testParameters.m_CenterY[i]); + centerZVec.push_back(testParameters.m_CenterZ[i]); + + sigmaXVec.push_back(testParameters.m_SigmaX[i]); + sigmaYVec.push_back(testParameters.m_SigmaY[i]); + sigmaZVec.push_back(testParameters.m_SigmaZ[i]); + + altitudeVec.push_back(testParameters.m_Altitude[i]); + } + + 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(testParameters.m_NumberOfGaussian); + + gaussianGenerator->AddGaussian(centerXVec, centerYVec, centerZVec, + sigmaXVec, sigmaYVec, sigmaZVec, altitudeVec); + + gaussianGenerator->Update(); + + image = gaussianGenerator->GetOutput(); + + mitk::CastToMitkImage(image, result); return result; } /** \brief Calculates hotspot statistics for given test image and ROI parameters. Uses ImageStatisticsCalculator to find a hotspot in a defined ROI within the given image. */ static mitk::ImageStatisticsCalculator::Statistics CalculateStatistics(mitk::Image* image, const Parameters& testParameters) { mitk::ImageStatisticsCalculator::Statistics result; + const unsigned int Dimension = 3; + typedef itk::Image MaskImageType; + MaskImageType::Pointer mask = MaskImageType::New(); + + MaskImageType::SizeType size; + MaskImageType::SpacingType spacing; + MaskImageType::IndexType start; + + mitk::ImageStatisticsCalculator::Pointer statisticsCalculator = mitk::ImageStatisticsCalculator::New(); + + for(int i = 0; i < 3; ++i) + { + spacing[i] = 1.00; + start[i] = 0.00; + } + + size[0] = testParameters.m_ImageColumns; + size[1] = testParameters.m_ImageRows; + size[2] = testParameters.m_ImageSlices; + + MaskImageType::RegionType region; + region.SetIndex(start); + region.SetSize(size); + + mask->SetRegions(region); + mask->Allocate(); + + for(int x = testParameters.m_RoiMinimumX; x < testParameters.m_RoiMaximumX; ++x) + { + for(int y = testParameters.m_RoiMinimumY; y < testParameters.m_RoiMaximumY; ++y) + { + for(int z = testParameters.m_RoiMinimumZ; z < testParameters.m_RoiMaximumZ; ++z) + { + MaskImageType::IndexType pixelIndex; + pixelIndex[0] = x; + pixelIndex[1] = y; + pixelIndex[2] = z; + + mask->SetPixel(pixelIndex, 1.00); + } + } + } + + mitk::Image::Pointer mitkMaskImage; + mitk::CastToMitkImage(mask, mitkMaskImage); + + statisticsCalculator->SetImage(image); + statisticsCalculator->SetImageMask(mitkMaskImage); + statisticsCalculator->SetMaskingModeToImage(); + statisticsCalculator->ComputeStatistics(); + result = statisticsCalculator->GetStatistics(); // create calculator object // fill parameters (mask, planar figure, etc.) // execute calculation // retrieve result and return from function // handle errors w/o crash! return result; } /** \brief Compares calculated against actual statistics values. Checks validness of all statistics aspects. Lets test fail if any aspect is not sufficiently equal. */ static void ValidateStatistics(const mitk::ImageStatisticsCalculator::Statistics& statistics, const Parameters& testParameters) { // check all expected test result against actual results - double actualValue = 0.0; - double expectedValue = testParameters.m_ExpectedValueU; + double actualPeakValue = testParameters.m_HotspotPeak; + double expectedPeakValue = statistics.HotspotPeak; + + double actualMaxValue = testParameters.m_HotspotMaximum; + double expectedMaxValue = statistics.HotspotMax; + + double actualMinValue = testParameters.m_HotspotMinimum; + double expectedMinValue = statistics.HotspotMin; + + + //Peak Index + std::vector actualPeakIndex = testParameters.m_HotspotPeakIndex; + vnl_vector expectedVnlPeakIndex; + expectedVnlPeakIndex = statistics.HotspotPeakIndex; + + vnl_vector actualVnlPeakIndex; + actualVnlPeakIndex.set_size(3); + + for(int i = 0; i < actualVnlPeakIndex.size(); ++i) + actualVnlPeakIndex[i] = actualPeakIndex[i]; + + // MaxIndex + std::vector actualMaxIndex = testParameters.m_HotspotMaximumIndex; + vnl_vector expectedVnlMaxIndex; + expectedVnlMaxIndex = statistics.HotspotMaxIndex; + + vnl_vector actualVnlMaxIndex; + actualVnlMaxIndex.set_size(3); + + for(int i = 0; i < actualVnlMaxIndex.size(); ++i) + actualVnlMaxIndex[i] = actualMaxIndex[i]; + + //MinIndex + std::vector actualMinIndex = testParameters.m_HotspotMinimumIndex; + vnl_vector expectedVnlMinIndex; + expectedVnlMinIndex = statistics.HotspotMinIndex; + + vnl_vector actualVnlMinIndex; + actualVnlMinIndex.set_size(3); + + for(int i = 0; i < actualVnlMinIndex.size(); ++i) + actualVnlMinIndex[i] = actualMinIndex[i]; + + double eps = 0.001; // float comparisons, allow tiny differences - MITK_TEST_CONDITION( ::fabs(actualValue - expectedValue) < mitk::eps, "Calculated statistics value " << actualValue << " (expected " << expectedValue << ")" ); + MITK_TEST_CONDITION( ::fabs(actualPeakValue - expectedPeakValue) < eps, "Actual hotspotPeak value " << actualPeakValue << " (expected " << expectedPeakValue << ")" ); + MITK_TEST_CONDITION( ::fabs(actualMaxValue - expectedMaxValue) < eps, "Actual hotspotMax value " << actualMaxValue << " (expected " << expectedMaxValue << ")" ); + MITK_TEST_CONDITION( ::fabs(actualMinValue - expectedMinValue) < eps, "Actual hotspotMin value " << actualMinValue << " (expected " << expectedMinValue << ")" ); + + MITK_TEST_CONDITION( expectedVnlPeakIndex == actualVnlPeakIndex, "Actual hotspotPeakIndex " << actualVnlPeakIndex << " (expected " << expectedVnlPeakIndex << ")" ); + MITK_TEST_CONDITION( expectedVnlMaxIndex == actualVnlMaxIndex, "Actual hotspotMaxIndex " << actualVnlMaxIndex << " (expected " << expectedVnlMaxIndex << ")" ); + MITK_TEST_CONDITION( expectedVnlMinIndex == actualVnlMinIndex, "Actual hotspotMinIndex " << actualVnlMinIndex << " (expected " << expectedVnlMinIndex << ")" ); } }; /** \brief Verifies that TODO hotspot statistics part of ImageStatisticsCalculator. bla... */ int mitkImageStatisticsHotspotTest(int argc, char* argv[]) { MITK_TEST_BEGIN("mitkImageStatisticsHotspotTest") // parse commandline parameters (see CMakeLists.txt) mitkImageStatisticsHotspotTestClass::Parameters parameters = mitkImageStatisticsHotspotTestClass::ParseParameters(argc,argv); // build a test image as described in parameters mitk::Image::Pointer image = mitkImageStatisticsHotspotTestClass::BuildTestImage(parameters); MITK_TEST_CONDITION_REQUIRED( image.IsNotNull(), "Generate test image" ); // calculate statistics for this image (potentially use parameters for statistics ROI) mitk::ImageStatisticsCalculator::Statistics statistics = mitkImageStatisticsHotspotTestClass::CalculateStatistics(image, parameters); // compare statistics against stored expected values mitkImageStatisticsHotspotTestClass::ValidateStatistics(statistics, parameters); MITK_TEST_END() } diff --git a/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp b/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp index a50222ad61..63386c8de7 100644 --- a/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp +++ b/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp @@ -1,2064 +1,2055 @@ /*=================================================================== 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 //#define DEBUG_HOTSPOTSEARCH #define _USE_MATH_DEFINES #include #if ( ( VTK_MAJOR_VERSION <= 5 ) && ( VTK_MINOR_VERSION<=8) ) #include "mitkvtkLassoStencilSource.h" #else #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 ); m_EmptyStatistics.Reset(); } ImageStatisticsCalculator::~ImageStatisticsCalculator() { } void ImageStatisticsCalculator::SetImage( const mitk::Image *image ) { if ( m_Image != image ) { m_Image = image; this->Modified(); unsigned int numberOfTimeSteps = image->GetTimeSteps(); // Initialize vectors to time-size of this image m_ImageHistogramVector.resize( numberOfTimeSteps ); m_MaskedImageHistogramVector.resize( numberOfTimeSteps ); m_PlanarFigureHistogramVector.resize( numberOfTimeSteps ); m_ImageStatisticsVector.resize( numberOfTimeSteps ); m_MaskedImageStatisticsVector.resize( numberOfTimeSteps ); m_PlanarFigureStatisticsVector.resize( numberOfTimeSteps ); m_ImageStatisticsTimeStampVector.resize( numberOfTimeSteps ); m_MaskedImageStatisticsTimeStampVector.resize( numberOfTimeSteps ); m_PlanarFigureStatisticsTimeStampVector.resize( numberOfTimeSteps ); m_ImageStatisticsCalculationTriggerVector.resize( numberOfTimeSteps ); m_MaskedImageStatisticsCalculationTriggerVector.resize( numberOfTimeSteps ); m_PlanarFigureStatisticsCalculationTriggerVector.resize( numberOfTimeSteps ); for ( unsigned int t = 0; t < image->GetTimeSteps(); ++t ) { m_ImageStatisticsTimeStampVector[t].Modified(); m_ImageStatisticsCalculationTriggerVector[t] = true; } } } void ImageStatisticsCalculator::SetImageMask( const mitk::Image *imageMask ) { if ( m_Image.IsNull() ) { itkExceptionMacro( << "Image needs to be set first!" ); } if ( m_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; } } } void ImageStatisticsCalculator::SetPlanarFigure( mitk::PlanarFigure *planarFigure ) { if ( m_Image.IsNull() ) { itkExceptionMacro( << "Image needs to be set first!" ); } if ( m_PlanarFigure != planarFigure ) { m_PlanarFigure = planarFigure; this->Modified(); for ( unsigned int t = 0; t < m_Image->GetTimeSteps(); ++t ) { m_PlanarFigureStatisticsTimeStampVector[t].Modified(); m_PlanarFigureStatisticsCalculationTriggerVector[t] = true; } } } void ImageStatisticsCalculator::SetMaskingMode( unsigned int mode ) { if ( m_MaskingMode != mode ) { m_MaskingMode = mode; m_MaskingModeChanged = true; this->Modified(); } } void ImageStatisticsCalculator::SetMaskingModeToNone() { if ( m_MaskingMode != MASKING_MODE_NONE ) { m_MaskingMode = MASKING_MODE_NONE; m_MaskingModeChanged = true; this->Modified(); } } void ImageStatisticsCalculator::SetMaskingModeToImage() { if ( m_MaskingMode != MASKING_MODE_IMAGE ) { m_MaskingMode = MASKING_MODE_IMAGE; m_MaskingModeChanged = true; this->Modified(); } } void ImageStatisticsCalculator::SetMaskingModeToPlanarFigure() { if ( m_MaskingMode != MASKING_MODE_PLANARFIGURE ) { m_MaskingMode = MASKING_MODE_PLANARFIGURE; m_MaskingModeChanged = true; this->Modified(); } } void ImageStatisticsCalculator::SetIgnorePixelValue(double value) { if ( m_IgnorePixelValue != value ) { m_IgnorePixelValue = value; if(m_DoIgnorePixelValue) { m_IgnorePixelValueChanged = true; } this->Modified(); } } double ImageStatisticsCalculator::GetIgnorePixelValue() { return m_IgnorePixelValue; } void ImageStatisticsCalculator::SetDoIgnorePixelValue(bool value) { if ( m_DoIgnorePixelValue != value ) { m_DoIgnorePixelValue = value; m_IgnorePixelValueChanged = true; this->Modified(); } } bool ImageStatisticsCalculator::GetDoIgnorePixelValue() { return m_DoIgnorePixelValue; } void ImageStatisticsCalculator::SetHotspotSize(double value) { m_HotspotSize = value; } double ImageStatisticsCalculator::GetHotspotSize() { return m_HotspotSize; } void ImageStatisticsCalculator::SetCalculateHotspot(bool value) { m_CalculateHotspot = value; } bool ImageStatisticsCalculator::IsHotspotCalculated() { return m_CalculateHotspot; } bool ImageStatisticsCalculator::ComputeStatistics( unsigned int timeStep ) { if (m_Image.IsNull() ) { mitkThrow() << "Image not set!"; } if (!m_Image->IsInitialized()) { mitkThrow() << "Image not initialized!"; } if ( m_Image->GetReferenceCount() == 1 ) { // Image no longer valid; we are the only ones to still hold a reference on it return false; } if ( timeStep >= m_Image->GetTimeSteps() ) { throw std::runtime_error( "Error: invalid time step!" ); } // If a mask was set but we are the only ones to still hold a reference on // it, delete it. if ( m_ImageMask.IsNotNull() && (m_ImageMask->GetReferenceCount() == 1) ) { m_ImageMask = NULL; } // Check if statistics is already up-to-date unsigned long imageMTime = m_ImageStatisticsTimeStampVector[timeStep].GetMTime(); unsigned long maskedImageMTime = m_MaskedImageStatisticsTimeStampVector[timeStep].GetMTime(); unsigned long planarFigureMTime = m_PlanarFigureStatisticsTimeStampVector[timeStep].GetMTime(); bool imageStatisticsCalculationTrigger = m_ImageStatisticsCalculationTriggerVector[timeStep]; bool maskedImageStatisticsCalculationTrigger = m_MaskedImageStatisticsCalculationTriggerVector[timeStep]; bool planarFigureStatisticsCalculationTrigger = m_PlanarFigureStatisticsCalculationTriggerVector[timeStep]; if ( !m_IgnorePixelValueChanged && ((m_MaskingMode != MASKING_MODE_NONE) || (imageMTime > m_Image->GetMTime() && !imageStatisticsCalculationTrigger)) && ((m_MaskingMode != MASKING_MODE_IMAGE) || (maskedImageMTime > m_ImageMask->GetMTime() && !maskedImageStatisticsCalculationTrigger)) && ((m_MaskingMode != MASKING_MODE_PLANARFIGURE) || (planarFigureMTime > m_PlanarFigure->GetMTime() && !planarFigureStatisticsCalculationTrigger)) ) { // Statistics is up to date! if ( m_MaskingModeChanged ) { m_MaskingModeChanged = false; 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 *statisticsContainer; HistogramContainer *histogramContainer; switch ( m_MaskingMode ) { case MASKING_MODE_NONE: default: if(!m_DoIgnorePixelValue) { statisticsContainer = &m_ImageStatisticsVector[timeStep]; histogramContainer = &m_ImageHistogramVector[timeStep]; m_ImageStatisticsTimeStampVector[timeStep].Modified(); m_ImageStatisticsCalculationTriggerVector[timeStep] = false; } else { statisticsContainer = &m_MaskedImageStatisticsVector[timeStep]; histogramContainer = &m_MaskedImageHistogramVector[timeStep]; m_MaskedImageStatisticsTimeStampVector[timeStep].Modified(); m_MaskedImageStatisticsCalculationTriggerVector[timeStep] = false; } break; case MASKING_MODE_IMAGE: statisticsContainer = &m_MaskedImageStatisticsVector[timeStep]; histogramContainer = &m_MaskedImageHistogramVector[timeStep]; m_MaskedImageStatisticsTimeStampVector[timeStep].Modified(); m_MaskedImageStatisticsCalculationTriggerVector[timeStep] = false; break; case MASKING_MODE_PLANARFIGURE: statisticsContainer = &m_PlanarFigureStatisticsVector[timeStep]; histogramContainer = &m_PlanarFigureHistogramVector[timeStep]; m_PlanarFigureStatisticsTimeStampVector[timeStep].Modified(); m_PlanarFigureStatisticsCalculationTriggerVector[timeStep] = false; break; } // Calculate statistics and histogram(s) if ( m_InternalImage->GetDimension() == 3 ) { if ( m_MaskingMode == MASKING_MODE_NONE && !m_DoIgnorePixelValue ) { AccessFixedDimensionByItk_2( m_InternalImage, InternalCalculateStatisticsUnmasked, 3, statisticsContainer, histogramContainer ); } else { AccessFixedDimensionByItk_3( m_InternalImage, InternalCalculateStatisticsMasked, 3, m_InternalImageMask3D.GetPointer(), statisticsContainer, histogramContainer ); } } else if ( m_InternalImage->GetDimension() == 2 ) { if ( m_MaskingMode == MASKING_MODE_NONE && !m_DoIgnorePixelValue ) { AccessFixedDimensionByItk_2( m_InternalImage, InternalCalculateStatisticsUnmasked, 2, statisticsContainer, histogramContainer ); } else { AccessFixedDimensionByItk_3( m_InternalImage, InternalCalculateStatisticsMasked, 2, m_InternalImageMask2D.GetPointer(), statisticsContainer, histogramContainer ); } } else { MITK_ERROR << "ImageStatistics: Image dimension not supported!"; } // Release unused image smart pointers to free memory m_InternalImage = mitk::Image::ConstPointer(); m_InternalImageMask3D = MaskImage3DType::Pointer(); m_InternalImageMask2D = MaskImage2DType::Pointer(); return true; } const ImageStatisticsCalculator::HistogramType * ImageStatisticsCalculator::GetHistogram( unsigned int timeStep, unsigned int label ) const { if ( m_Image.IsNull() || (timeStep >= m_Image->GetTimeSteps()) ) { return NULL; } switch ( m_MaskingMode ) { case MASKING_MODE_NONE: default: { if(m_DoIgnorePixelValue) return m_MaskedImageHistogramVector[timeStep][label]; return m_ImageHistogramVector[timeStep][label]; } case MASKING_MODE_IMAGE: return m_MaskedImageHistogramVector[timeStep][label]; case MASKING_MODE_PLANARFIGURE: return m_PlanarFigureHistogramVector[timeStep][label]; } } const ImageStatisticsCalculator::HistogramContainer & ImageStatisticsCalculator::GetHistogramVector( unsigned int timeStep ) const { if ( m_Image.IsNull() || (timeStep >= m_Image->GetTimeSteps()) ) { return m_EmptyHistogramContainer; } switch ( m_MaskingMode ) { case MASKING_MODE_NONE: default: { if(m_DoIgnorePixelValue) return m_MaskedImageHistogramVector[timeStep]; return m_ImageHistogramVector[timeStep]; } case MASKING_MODE_IMAGE: return m_MaskedImageHistogramVector[timeStep]; case MASKING_MODE_PLANARFIGURE: return m_PlanarFigureHistogramVector[timeStep]; } } const ImageStatisticsCalculator::Statistics & ImageStatisticsCalculator::GetStatistics( unsigned int timeStep, unsigned int label ) const { if ( m_Image.IsNull() || (timeStep >= m_Image->GetTimeSteps()) ) { return m_EmptyStatistics; } switch ( m_MaskingMode ) { case MASKING_MODE_NONE: default: { if(m_DoIgnorePixelValue) return m_MaskedImageStatisticsVector[timeStep][label]; return m_ImageStatisticsVector[timeStep][label]; } case MASKING_MODE_IMAGE: return m_MaskedImageStatisticsVector[timeStep][label]; case MASKING_MODE_PLANARFIGURE: return m_PlanarFigureStatisticsVector[timeStep][label]; } } const ImageStatisticsCalculator::Statistics & ImageStatisticsCalculator::GetHotspotStatistics( unsigned int timeStep, unsigned int label) const { if(m_Image.IsNull() || (timeStep >= m_Image->GetTimeSteps()) ) return m_EmptyStatistics; return m_MaskedImageHotspotStatisticsVector[timeStep][label]; } const ImageStatisticsCalculator::StatisticsContainer & ImageStatisticsCalculator::GetStatisticsVector( unsigned int timeStep ) const { if ( m_Image.IsNull() || (timeStep >= m_Image->GetTimeSteps()) ) { return m_EmptyStatisticsContainer; } switch ( m_MaskingMode ) { case MASKING_MODE_NONE: default: { if(m_DoIgnorePixelValue) return m_MaskedImageStatisticsVector[timeStep]; return m_ImageStatisticsVector[timeStep]; } case MASKING_MODE_IMAGE: return m_MaskedImageStatisticsVector[timeStep]; case MASKING_MODE_PLANARFIGURE: return m_PlanarFigureStatisticsVector[timeStep]; } } void ImageStatisticsCalculator::ExtractImageAndMask( unsigned int timeStep ) { if ( m_Image.IsNull() ) { throw std::runtime_error( "Error: image empty!" ); } if ( timeStep >= m_Image->GetTimeSteps() ) { throw std::runtime_error( "Error: invalid time step!" ); } ImageTimeSelector::Pointer imageTimeSelector = ImageTimeSelector::New(); imageTimeSelector->SetInput( m_Image ); imageTimeSelector->SetTimeNr( timeStep ); imageTimeSelector->UpdateLargestPossibleRegion(); mitk::Image *timeSliceImage = imageTimeSelector->GetOutput(); switch ( m_MaskingMode ) { case MASKING_MODE_NONE: { m_InternalImage = timeSliceImage; m_InternalImageMask2D = NULL; m_InternalImageMask3D = NULL; if(m_DoIgnorePixelValue) { 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); } } break; } case MASKING_MODE_IMAGE: { if ( m_ImageMask.IsNotNull() && (m_ImageMask->GetReferenceCount() > 1) ) { 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!" ); } } else { throw std::runtime_error( "Error: image mask empty!" ); } break; } case MASKING_MODE_PLANARFIGURE: { m_InternalImageMask2D = NULL; 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 Geometry3D *imageGeometry = timeSliceImage->GetGeometry(); if ( imageGeometry == NULL ) { throw std::runtime_error( "Image geometry invalid!" ); } const Geometry2D *planarFigureGeometry2D = m_PlanarFigure->GetGeometry2D(); if ( planarFigureGeometry2D == NULL ) { throw std::runtime_error( "Planar-Figure not yet initialized!" ); } 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; // Extract slice with given position and direction from image unsigned int dimension = timeSliceImage->GetDimension(); if (dimension != 2) { ExtractImageFilter::Pointer imageExtractor = ExtractImageFilter::New(); imageExtractor->SetInput( timeSliceImage ); imageExtractor->SetSliceDimension( axis ); imageExtractor->SetSliceIndex( slice ); imageExtractor->Update(); m_InternalImage = imageExtractor->GetOutput(); } else { m_InternalImage = timeSliceImage; } // Compute mask from PlanarFigure AccessFixedDimensionByItk_1( m_InternalImage, InternalCalculateMaskFromPlanarFigure, 2, axis ); } } if(m_DoIgnorePixelValue) { if ( m_InternalImage->GetDimension() == 3 ) { AccessFixedDimensionByItk_1( m_InternalImage, InternalMaskIgnoredPixels, 3, m_InternalImageMask3D.GetPointer() ); } else if ( m_InternalImage->GetDimension() == 2 ) { AccessFixedDimensionByItk_1( m_InternalImage, InternalMaskIgnoredPixels, 2, m_InternalImageMask2D.GetPointer() ); } } } bool ImageStatisticsCalculator::GetPrincipalAxis( const Geometry3D *geometry, Vector3D vector, unsigned int &axis ) { vector.Normalize(); for ( unsigned int i = 0; i < 3; ++i ) { Vector3D axisVector = geometry->GetAxisVector( i ); axisVector.Normalize(); if ( fabs( fabs( axisVector * vector ) - 1.0) < mitk::eps ) { axis = i; return true; } } return false; } 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; 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; iGetIndexOfMaximum()[i]; statistics.MinIndex[i] = minMaxFilter->GetIndexOfMinimum()[i]; } 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(); histogramContainer->push_back( histogramGenerator->GetOutput() ); } template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator::InternalMaskIgnoredPixels( const itk::Image< TPixel, VImageDimension > *image, itk::Image< unsigned short, VImageDimension > *maskImage ) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< unsigned short, VImageDimension > MaskImageType; itk::ImageRegionIterator itmask(maskImage, maskImage->GetLargestPossibleRegion()); itk::ImageRegionConstIterator itimage(image, image->GetLargestPossibleRegion()); itmask.GoToBegin(); itimage.GoToBegin(); while( !itmask.IsAtEnd() ) { if(m_IgnorePixelValue == itimage.Get()) { itmask.Set(0); } ++itmask; ++itimage; } } template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator::InternalCalculateStatisticsMasked( const itk::Image< TPixel, VImageDimension > *image, itk::Image< unsigned short, VImageDimension > *maskImage, StatisticsContainer* statisticsContainer, HistogramContainer* histogramContainer ) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< unsigned short, VImageDimension > MaskImageType; typedef typename ImageType::IndexType IndexType; typedef typename ImageType::PointType PointType; typedef typename ImageType::SpacingType SpacingType; typedef itk::LabelStatisticsImageFilter< ImageType, MaskImageType > LabelStatisticsFilterType; typedef itk::ChangeInformationImageFilter< MaskImageType > ChangeInformationFilterType; typedef itk::ExtractImageFilter< ImageType, ImageType > ExtractImageFilterType; statisticsContainer->clear(); histogramContainer->clear(); // Make sure that mask is set if ( maskImage == 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 orientation of mask and image are the same typedef typename ImageType::DirectionType DirectionType; DirectionType imageDirection = image->GetDirection(); DirectionType maskDirection = maskImage->GetDirection(); for( int i = 0; i < imageDirection.ColumnDimensions; ++i ) { for( int j = 0; j < imageDirection.ColumnDimensions; ++j ) { double differenceDirection = imageDirection[i][j] - maskDirection[i][j]; if ( fabs( differenceDirection ) > mitk::eps ) { itkExceptionMacro( << "Mask needs to have same direction as image! (Image direction: " << imageDirection << "; Mask direction: " << maskDirection << ")" ); } } } // Make sure that the voxels of mask and image are correctly "aligned", i.e., voxel boundaries are the same in both images PointType imageOrigin = image->GetOrigin(); PointType maskOrigin = maskImage->GetOrigin(); long offset[ImageType::ImageDimension]; typedef itk::ContinuousIndex ContinousIndexType; ContinousIndexType maskOriginContinousIndex, imageOriginContinousIndex; image->TransformPhysicalPointToContinuousIndex(maskOrigin, maskOriginContinousIndex); image->TransformPhysicalPointToContinuousIndex(imageOrigin, imageOriginContinousIndex); for ( unsigned int i = 0; i < ImageType::ImageDimension; ++i ) { double misalignment = maskOriginContinousIndex[i] - floor( maskOriginContinousIndex[i] + 0.5 ); if ( fabs( misalignment ) > mitk::eps ) { 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]; } // 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() << ")" ); } // If mask region is smaller than image region, extract the sub-sampled region from the original image typename ImageType::SizeType imageSize = image->GetBufferedRegion().GetSize(); typename ImageType::SizeType maskSize = maskImage->GetBufferedRegion().GetSize(); bool maskSmallerImage = false; for ( unsigned int i = 0; i < ImageType::ImageDimension; ++i ) { if ( maskSize[i] < imageSize[i] ) { maskSmallerImage = true; } } typename ImageType::ConstPointer adaptedImage; if ( maskSmallerImage ) { typename ExtractImageFilterType::Pointer extractImageFilter = ExtractImageFilterType::New(); extractImageFilter->SetInput( image ); extractImageFilter->SetExtractionRegion( adaptedMaskImage->GetBufferedRegion() ); extractImageFilter->Update(); adaptedImage = extractImageFilter->GetOutput(); } else { adaptedImage = image; } // Initialize Filter typedef itk::StatisticsImageFilter< ImageType > StatisticsFilterType; typename StatisticsFilterType::Pointer statisticsFilter = StatisticsFilterType::New(); statisticsFilter->SetInput( adaptedImage ); 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 ) ) { relevantLabels.push_back( i ); maskNonEmpty = true; } } 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; 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. 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; } 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**********************************************/ - Statistics hotspotStatistics = CalculateHotspotStatistics (adaptedImage.GetPointer(), adaptedMaskImage.GetPointer(), 3.00); + Statistics hotspotStatistics = CalculateHotspotStatistics (adaptedImage.GetPointer(), adaptedMaskImage.GetPointer(), 6.2035049089940); 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; typedef itk::ImageRegionConstIterator MaskImageIteratorType; typedef itk::ImageRegionConstIteratorWithIndex InputImageIndexIteratorType; MaskImageIteratorType maskIt(maskImage, maskImage->GetLargestPossibleRegion()); InputImageIndexIteratorType imageIndexIt(inputImage, inputImage->GetLargestPossibleRegion()); float maxValue = itk::NumericTraits::min(); float minValue = itk::NumericTraits::max(); ImageType::IndexType maxIndex; ImageType::IndexType minIndex; double maxBound = itk::NumericTraits::min(); double minBound = itk::NumericTraits::max(); std::list< ImageType::IndexType > maxIndexList; std::list< ImageType::IndexType > minIndexList; 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) { maxIndexList.clear(); - maxIndexList.push_back( index ); + maxIndexList.push_back(index); maxValue = value; } else if ( value == maxValue ) { maxIndexList.push_back( index ); } if(value < minValue && value > minBoundary) { minIndexList.clear(); - minIndexList.push_back( index ); + minIndexList.push_back(index); minValue = value; } else if ( value == minValue ) { minIndexList.push_back( index ); } } } std::list> tmpList; + MinMaxIndex minMax; + + // Convert Max-IndexType-list into vnl_vector-list for(std::list< ImageType::IndexType >::iterator it = maxIndexList.begin(); it != maxIndexList.end(); ++it) { vnl_vector indexVector; indexVector.set_size(inputImage->GetImageDimension()); ImageType::IndexType index = *it; for(int i = 0; i < VImageDimension ; ++i) indexVector[i] = index[i]; tmpList.push_back(indexVector); } - MinMaxIndex minMax; minMax.MaxIndexList = tmpList; - tmpList.clear(); + // Convert Min-IndexType-list into vnl_vector-list for(std::list< ImageType::IndexType >::iterator it = minIndexList.begin(); it != minIndexList.end(); ++it) { vnl_vector indexVector; indexVector.set_size(inputImage->GetImageDimension()); ImageType::IndexType index = *it; for(int i = 0; i < indexVector.size(); ++i) indexVector[i] = index[i]; tmpList.push_back(indexVector); } minMax.MinIndexList = tmpList; minMax.MinIndex.set_size(inputImage->GetImageDimension()); minMax.MaxIndex.set_size(inputImage->GetImageDimension()); 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; MaskImageType::Pointer convolutionMask = MaskImageType::New(); SpacingType spacing = inputImage->GetSpacing(); IndexType start; start.Fill(0); ContinuousIndexType origin; SizeType size; ContinuousIndexType convolutionMaskCenterCoordinate; /*****************************************************Creating convolution 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 { size[i] = castedIndex +1; } convolutionMaskCenterCoordinate[i] = (size[i] -1) / 2; } RegionType region; region.SetSize(size); region.SetIndex(start); convolutionMask->SetRegions(region); convolutionMask->SetSpacing(spacing); convolutionMask->Allocate(); float subPixelDimensionX = spacing[0]/2; float subPixelDimensionY = spacing[1]/2; double distance = 0.0; double countSubPixel = 0.0; double pixelValue = 0.0; 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) { for(double y = indexPoint[1] - subPixelDimensionY /2; y <= indexPoint[1] + subPixelDimensionY / 2; y += subPixelDimensionY) { ContinuousIndexType subPixelCenterCoordinate; subPixelCenterCoordinate[0] = x; subPixelCenterCoordinate[1] = y; for( unsigned int i = 2; i < VImageDimension; ++i ) { subPixelCenterCoordinate[i] = indexPoint[i]; } distance = subPixelCenterCoordinate.EuclideanDistanceTo(convolutionMaskCenterCoordinate); if(distance <= RadiusInMM) countSubPixel++; } } // pixelValue is the counted SubPixels divided by factor 4 pixelValue = countSubPixel / 4.00; convolutionMask->SetPixel(maskIt.GetIndex(),pixelValue); countSubPixel = 0.00; } /*****************************************************Creating Peak Image**********************************************/ typedef itk::ConvolutionImageFilter FilterType; FilterType::Pointer convolutionFilter = FilterType::New(); convolutionFilter->SetInput(inputImage); convolutionFilter->SetKernelImage(convolutionMask); convolutionFilter->SetNormalize(true); convolutionFilter->Update(); MaskImageType::Pointer peakImage = convolutionFilter->GetOutput(); /*std::cout << std::endl << std::endl; std::cout << "PeakImage: " << std::endl; unsigned int lastZ = 1000000000; unsigned int lastY = 1000000000; unsigned int hotspotMaskIndexCounter = 0; MaskIteratorType hotspotMaskIt(peakImage, peakImage->GetLargestPossibleRegion() ); for(hotspotMaskIt.GoToBegin();!hotspotMaskIt.IsAtEnd();++hotspotMaskIt) { 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]; } hotspotMaskIndexCounter++; if(hotspotMaskIndexCounter > 899) { std::cout << std::endl; hotspotMaskIndexCounter = 0; } } std::cout << std::endl << std::endl;*/ /*****************************************************Creating Hotspot Sphere**********************************************/ SphereMaskImageType::Pointer hotspotSphere = SphereMaskImageType::New(); SphereMaskImageType::SpacingType hotspotSphereSpacing = peakImage->GetSpacing(); SphereMaskImageType::SizeType hotspotSphereSize; SphereMaskImageType::RegionType regionSphere; SphereMaskImageType::IndexType hotspotPeakIndex; PointType hotspotOrigin; IndexType offsetInIndex; typedef itk::EllipseSpatialObject EllipseType; typedef itk::SpatialObjectToImageFilter SpatialObjectToImageFilter; const double radiusSUVHotspot = 6.2035049089940; float maxBoundary = itk::NumericTraits::max(); float minBoundary = itk::NumericTraits::min(); double hotspotPeak = itk::NumericTraits::min(); bool isInside = false; std::list> peakIndexList; - /* vnl_vector headIndex; - headIndex.set_size(VImageDimension); - - peakIndexList.push_back(head);*/ - int i = 0; while (isInside == false) { MinMaxIndex peakInformations = CalculateMinMaxIndex(peakImage.GetPointer(), maskImage, minBoundary, maxBoundary); peakIndexList = peakInformations.MaxIndexList; - for(std::list>::iterator it = peakIndexList.begin(); it != peakIndexList.end(); ++it) - { - vnl_vector tmp; - tmp = *it; - std::cout << "List[" << i << "]: " << tmp << std::endl; - } - - if(maxBoundary == peakInformations.Max) - { - std::cout << "Abortion" << std::endl; - break; - } - - if(!peakIndexList.empty()) + for(std::list>::iterator iterator = peakIndexList.begin(); iterator != peakIndexList.end(); ++iterator) { - for(std::list>::iterator iterator = peakIndexList.begin(); iterator != peakIndexList.end(); ++iterator) + /*if(maxBoundary == peakInformations.Max) { - vnl_vector tmpIndex = *iterator; - SphereMaskImageType::IndexType peakIndex; + std::cout << "Abortion" << std::endl; + break; + }*/ - for(unsigned int i = 0; i < VImageDimension; ++i) - { - // Converting vnl_vector in IndexType for origin of sphere - peakIndex[i] = tmpIndex[i]; + vnl_vector tmpIndex = *iterator; + SphereMaskImageType::IndexType peakIndex; - double countIndex = 2.0 * radiusSUVHotspot / hotspotSphereSpacing[i]; + for(unsigned int i = 0; i < VImageDimension; ++i) + { + // Converting vnl_vector in IndexType for origin of sphere + peakIndex[i] = tmpIndex[i]; - // Rounding up to the next integer by cast - countIndex += 0.9; - int castedIndex = static_cast(countIndex); + double countIndex = 2.0 * radiusSUVHotspot / hotspotSphereSpacing[i]; - // 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; - } + // 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 ) + { + hotspotSphereSize[i] = castedIndex; + } + else + { + hotspotSphereSize[i] = castedIndex +1; } + } - // Initialize SpatialObjectoToImageFilter - itk::SpatialObjectToImageFilter::Pointer spatialObjectToImageFilter - = SpatialObjectToImageFilter::New(); - spatialObjectToImageFilter->SetSize(hotspotSphereSize); - spatialObjectToImageFilter->SetSpacing(hotspotSphereSpacing); - // Creating spatial sphere object - EllipseType::Pointer sphere = EllipseType::New(); - sphere->SetRadius(radiusSUVHotspot); + // Initialize SpatialObjectoToImageFilter + itk::SpatialObjectToImageFilter::Pointer spatialObjectToImageFilter + = SpatialObjectToImageFilter::New(); - typedef EllipseType::TransformType TransformType; - TransformType::Pointer transform = TransformType::New(); + spatialObjectToImageFilter->SetSize(hotspotSphereSize); + spatialObjectToImageFilter->SetSpacing(hotspotSphereSpacing); - transform->SetIdentity(); + // Creating spatial sphere object + EllipseType::Pointer sphere = EllipseType::New(); + sphere->SetRadius(radiusSUVHotspot); + typedef EllipseType::TransformType TransformType; + TransformType::Pointer transform = TransformType::New(); - TransformType::OutputVectorType translation; + transform->SetIdentity(); - // 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; + TransformType::OutputVectorType translation; + + // 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; - transform->Translate(translation, false); + transform->Translate(translation, false); - sphere->SetObjectToParentTransform(transform); + sphere->SetObjectToParentTransform(transform); - spatialObjectToImageFilter->SetInput(sphere); + spatialObjectToImageFilter->SetInput(sphere); - sphere->SetDefaultInsideValue(1.00); - sphere->SetDefaultOutsideValue(0.00); + sphere->SetDefaultInsideValue(1.00); + sphere->SetDefaultOutsideValue(0.00); - spatialObjectToImageFilter->SetUseObjectValue(true); - spatialObjectToImageFilter->SetOutsideValue(0); + spatialObjectToImageFilter->SetUseObjectValue(true); + spatialObjectToImageFilter->SetOutsideValue(0); - spatialObjectToImageFilter->Update(); - hotspotSphere = spatialObjectToImageFilter->GetOutput(); + spatialObjectToImageFilter->Update(); + hotspotSphere = spatialObjectToImageFilter->GetOutput(); #ifdef DEBUG_HOTSPOTSEARCH - std::cout << std::endl << std::endl; - std::cout << "hotspotMask: " << std::endl; - unsigned int lastZ = 1000000000; - unsigned int lastY = 1000000000; + std::cout << std::endl << std::endl; + std::cout << "hotspotMask: " << std::endl; + unsigned int lastZ = 1000000000; + unsigned int lastY = 1000000000; - unsigned int hotspotMaskIndexCounter = 0; - SphereMaskIteratorType hotspotMaskIt(hotspotSphere, hotspotSphere->GetLargestPossibleRegion() ); - for(hotspotMaskIt.GoToBegin();!hotspotMaskIt.IsAtEnd();++hotspotMaskIt) - { + unsigned int hotspotMaskIndexCounter = 0; + SphereMaskIteratorType hotspotMaskIt(hotspotSphere, hotspotSphere->GetLargestPossibleRegion() ); + for(hotspotMaskIt.GoToBegin();!hotspotMaskIt.IsAtEnd();++hotspotMaskIt) + { - 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]; - } + 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]; + } - hotspotMaskIndexCounter++; + hotspotMaskIndexCounter++; - if(hotspotMaskIndexCounter > 168) { - std::cout << std::endl; - hotspotMaskIndexCounter = 0; - } + if(hotspotMaskIndexCounter > 168) { + std::cout << std::endl; + hotspotMaskIndexCounter = 0; } + } - std::cout << std::endl << std::endl; + std::cout << std::endl << std::endl; #endif - // Calculate new origin for hotspot sphere + // Calculate new origin for hotspot sphere - for(int i = 0; i < VImageDimension; ++i) - offsetInIndex[i] = (hotspotSphereSize[i] -1) / 2; + for(int i = 0; i < VImageDimension; ++i) + offsetInIndex[i] = (hotspotSphereSize[i] -1) / 2; - peakImage->TransformIndexToPhysicalPoint(peakIndex, hotspotOrigin); + peakImage->TransformIndexToPhysicalPoint(peakIndex, hotspotOrigin); - PointType offsetInPhysicalPoint; - hotspotSphere->TransformIndexToPhysicalPoint(offsetInIndex, offsetInPhysicalPoint); + PointType offsetInPhysicalPoint; + hotspotSphere->TransformIndexToPhysicalPoint(offsetInIndex, offsetInPhysicalPoint); - for(int i = 0; i < VImageDimension; ++i) - hotspotOrigin[i] -= offsetInPhysicalPoint[i]; + for(int i = 0; i < VImageDimension; ++i) + hotspotOrigin[i] -= offsetInPhysicalPoint[i]; - hotspotSphere->SetOrigin(hotspotOrigin); - hotspotSphere->Allocate(); + hotspotSphere->SetOrigin(hotspotOrigin); + hotspotSphere->Allocate(); - // If sphere is not in peakImage we need a new peakValue - bool sphereInRegion = IsSphereInsideRegion(peakImage.GetPointer(), hotspotSphere.GetPointer() ); + /*vnl_vector tmp; + tmp = *iterator; + std::cout << "PeakValue: " << peakInformations.Max << " List[" << i << "]: " << tmp << std::endl;*/ - if(sphereInRegion == false) - { - maxBoundary = peakInformations.Max; - minBoundary = peakInformations.Min; - } - else - { - hotspotPeak = peakInformations.Max; - - for(int i = 0; i < VImageDimension; ++i) - hotspotPeakIndex[i] = peakIndex[i]; + // If sphere is not in peakImage we need a new peakValue + bool sphereInRegion = IsSphereInsideRegion(peakImage.GetPointer(), hotspotSphere.GetPointer() ); - isInside = true; - break; - } + if(sphereInRegion == false) + { + maxBoundary = peakInformations.Max; + minBoundary = peakInformations.Min; } + else + { + hotspotPeak = peakInformations.Max; + + for(int i = 0; i < VImageDimension; ++i) + hotspotPeakIndex[i] = peakIndex[i]; - i++; + isInside = true; + break; + } } } /*********************************Creating cropped inputImage for calculation of hotspot statistics****************************/ IndexType croppedStart; for(int i = 0; i < 3; ++i) croppedStart[i] = offsetInIndex[i] - hotspotSphereSize[i]; 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); ImageType::RegionType outputRegion; ImageType::IndexType outputStart; outputStart.Fill(0); outputRegion.SetIndex(outputStart); outputRegion.SetSize(hotspotSphere->GetLargestPossibleRegion().GetSize()); ImageType::Pointer croppedInputImage = ImageType::New(); croppedInputImage->SetRegions(outputRegion); croppedInputImage->Allocate(); typedef itk::ImageRegionConstIterator InputImageIteratorType; InputImageIteratorType inputIt(inputImage, inputRegion); InputImageIteratorType croppedInputImageIt(croppedInputImage, outputRegion); for(inputIt.GoToBegin(), croppedInputImageIt.GoToBegin(); !inputIt.IsAtEnd(); ++inputIt, ++croppedInputImageIt) { croppedInputImage->SetPixel(croppedInputImageIt.GetIndex(), inputIt.Get()); - //croppedInputImageIt.Set(inputIt.Get() ); } // Calculate statistics in Hotspot maxBoundary = itk::NumericTraits::max(); minBoundary = itk::NumericTraits::min(); MinMaxIndex hotspotInformations; Statistics hotspotStatistics; std::list> maxIndexList; hotspotInformations = CalculateMinMaxIndex(croppedInputImage.GetPointer(), hotspotSphere.GetPointer(), minBoundary, maxBoundary); // Min and max index is the first item in the corresponding list hotspotInformations.MaxIndex = hotspotInformations.MaxIndexList.front(); hotspotInformations.MinIndex = hotspotInformations.MinIndexList.front(); // 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; hotspotStatistics.HotspotPeakIndex.set_size(inputImage->GetImageDimension()); for (int i = 0; i< hotspotStatistics.HotspotPeakIndex.size(); ++i) { hotspotStatistics.HotspotPeakIndex[i] = hotspotPeakIndex[i]; } - return hotspotStatistics; } 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; RegionType region(inputImage->GetLargestPossibleRegion().GetIndex(), inputImage->GetLargestPossibleRegion().GetSize() ); typedef itk::Image SphereImageType; typedef SphereImageType::SizeType SphereSizeType; typedef SphereImageType::IndexType SphereIndexType; SphereIndexType sphereCenterIndex; SphereSizeType sphereSize = sphereImage->GetLargestPossibleRegion().GetSize(); double radius = (sphereSize[0] - 1)/2.00; bool isInside = false; for(int i = 0; i < VImageDimension; ++i) sphereCenterIndex[i] = (sphereSize[i] -1) / 2; /**********************************Calculation of sphere-coordinates which are on axis of coordinates****************************/ 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]; } if(i == 1) { sphereIndex[0] = sphereCenterIndex[0]; sphereIndex[1] = j; sphereIndex[2] = sphereCenterIndex[2]; } if(i == 2) { sphereIndex[0] = sphereCenterIndex[0]; sphereIndex[1] = sphereCenterIndex[1]; sphereIndex[2] = j; } PointType spherePoint; sphereImage->TransformIndexToPhysicalPoint(sphereIndex, spherePoint); IndexType regionIndex; inputImage->TransformPhysicalPointToIndex(spherePoint, regionIndex); 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; // First run: summation, second run: subtraction bool addOrSub = false; + 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)); + for(int i = 0; i < 2; ++i) { if(addOrSub == false) { // Point 1 - sphereIndex[0] = sphereCenterIndex[0] + radius * sin(45.00 * itk::Math::pi / 180.00) * cos(45.00 * itk::Math::pi / 180.00); - sphereIndex[1] = sphereCenterIndex[1] + radius * sin(45.00 * itk::Math::pi / 180.00) * sin(45.00 * itk::Math::pi / 180.00); - sphereIndex[2] = sphereCenterIndex[2] + radius * cos(45.00 * itk::Math::pi / 180.00); + 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; } } else { // Point 2 - sphereIndex[0] = sphereCenterIndex[0] + radius * sin(45.00 * itk::Math::pi / 180.00) * cos(45.00 * itk::Math::pi / 180.00); - sphereIndex[1] = -(sphereCenterIndex[1] + radius * sin(45.00 * itk::Math::pi / 180.00) * sin(45.00 * itk::Math::pi / 180.00)); - sphereIndex[2] = sphereCenterIndex[2] + radius * cos(45.00 * itk::Math::pi / 180.00); + 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; } } // Point 3/4 - sphereIndex[2] = -(sphereCenterIndex[2] + radius * cos(45.00 * itk::Math::pi / 180.00)); + 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; } // Point 5/6 - sphereIndex[0] = -(sphereCenterIndex[0] + radius * sin(45.00 * itk::Math::pi / 180.00) * cos(45.00 * itk::Math::pi / 180.00)); + sphereIndex[0] = sphereCenterIndex[0] - offsetXDirection; sphereImage->TransformIndexToPhysicalPoint(sphereIndex, spherePoint); inputImage->TransformPhysicalPointToIndex(spherePoint, regionIndex); if(region.IsInside(regionIndex) == true) { isInside = true; } else { isInside = false; return isInside; } // Point 7/8 - sphereIndex[2] = sphereCenterIndex[2] + radius * cos(45.00 * itk::Math::pi / 180.00); + 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; } addOrSub = true; } return isInside; } 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 ); points->InsertNextPoint( point3D[i0], point3D[i1], 0 ); } // mark a malformed 2D planar figure ( i.e. area = 0 ) as out of bounds // this can happen when all control points of a rectangle lie on the same line = two of the three extents are zero double bounds[6] = {0, 0, 0, 0, 0, 0}; points->GetBounds( bounds ); bool extent_x = (fabs(bounds[0] - bounds[1])) < mitk::eps; bool extent_y = (fabs(bounds[2] - bounds[3])) < mitk::eps; bool extent_z = (fabs(bounds[4] - bounds[5])) < mitk::eps; // throw an exception if a closed planar figure is deformed, i.e. has only one non-zero extent if ( m_PlanarFigure->IsClosed() && ((extent_x && extent_y) || (extent_x && extent_z) || (extent_y && extent_z))) { mitkThrow() << "Figure has a zero area and cannot be used for masking."; } if ( outOfBounds ) { throw std::runtime_error( "Figure at least partially outside of image bounds!" ); } // create a vtkLassoStencilSource and set the points of the Polygon vtkSmartPointer lassoStencil = vtkSmartPointer::New(); lassoStencil->SetShapeToPolygon(); lassoStencil->SetPoints( points ); // Export from ITK to VTK (to use a VTK filter) typedef itk::VTKImageImport< MaskImage2DType > ImageImportType; typedef itk::VTKImageExport< MaskImage2DType > ImageExportType; typename ImageExportType::Pointer itkExporter = ImageExportType::New(); itkExporter->SetInput( castFilter->GetOutput() ); vtkSmartPointer vtkImporter = vtkSmartPointer::New(); this->ConnectPipelines( itkExporter, vtkImporter ); // Apply the generated image stencil to the input image vtkSmartPointer imageStencilFilter = vtkSmartPointer::New(); imageStencilFilter->SetInputConnection( vtkImporter->GetOutputPort() ); imageStencilFilter->SetStencil( lassoStencil->GetOutput() ); imageStencilFilter->ReverseStencilOff(); imageStencilFilter->SetBackgroundValue( 0 ); imageStencilFilter->Update(); // 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(); // Store mask m_InternalImageMask2D = itkImporter->GetOutput(); } void ImageStatisticsCalculator::UnmaskedStatisticsProgressUpdate() { // Need to throw away every second progress event to reach a final count of // 100 since two consecutive filters are used in this case static int updateCounter = 0; if ( updateCounter++ % 2 == 0 ) { this->InvokeEvent( itk::ProgressEvent() ); } } void ImageStatisticsCalculator::MaskedStatisticsProgressUpdate() { this->InvokeEvent( itk::ProgressEvent() ); } }