diff --git a/Modules/ImageStatistics/Testing/mitkImageStatisticsCalculatorTest.cpp b/Modules/ImageStatistics/Testing/mitkImageStatisticsCalculatorTest.cpp index 0e3c62345e..780d8f130c 100644 --- a/Modules/ImageStatistics/Testing/mitkImageStatisticsCalculatorTest.cpp +++ b/Modules/ImageStatistics/Testing/mitkImageStatisticsCalculatorTest.cpp @@ -1,496 +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 /** * \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(); } // 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() { mitk::StandardFileLocations::Pointer locator = mitk::StandardFileLocations::GetInstance(); 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 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; } int mitkImageStatisticsCalculatorTest(int, char* []) { // always start with this! MITK_TEST_BEGIN("mitkImageStatisticsCalculatorTest") //QCoreApplication app(argc, argv); 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 { + public: // 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" float m_Spacing[3]; // XML-Tag "spacingX", "spacingY", "spacingZ" // 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 int m_MaxSizeX; // XML-Tag "maximumSizeX" int m_MinSizeX; // XML-Tag "minimumSizeX" int m_MaxSizeY; // XML-Tag "maximumSizeY" int m_MinSizeY; // XML-Tag "minimumSizeY" int m_MaxSizeZ; // XML-Tag "maximumSizeZ" int m_MinSizeZ; // XML-Tag "minimumSizeZ" float m_hotspotRadiusInMM; // XML-Tag "hotspotRadiusInMM" //XML-Tag float m_HotspotMin; // XML-Tag "minimum" float m_HotspotMax; // XML-Tag "maximum" float m_HotspotMean; // XML-Tag "mean" vnl_vector m_HotspotMaxIndex; // XML-Tag "maximumIndexX", XML-Tag "maximumIndexY", XML-Tag "maximumIndexZ" vnl_vector m_HotspotMinIndex; // XML-Tag "minimumIndexX", XML-Tag "minimumIndexY", XML-Tag "minimumIndexZ" vnl_vector m_HotspotIndex; // XML-Tag "hotspotIndexX", XML-Tag "hotspotIndexY", XML-Tag "hotspotIndexZ" }; /** \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 NodeList testimages; domRoot->GetChildren("testimage", testimages); MITK_TEST_CONDITION_REQUIRED( testimages.size() == 1, "One test image defined" ) itk::DOMNode* testimage = testimages[0]; 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" ); result.m_Spacing[0] = GetDoubleAttribute(testimage, "spacingX"); result.m_Spacing[1] = GetDoubleAttribute(testimage, "spacingY"); result.m_Spacing[2] = GetDoubleAttribute(testimage, "spacingZ"); MITK_TEST_OUTPUT( << "Read size parameters (x,y,z): " << result.m_ImageRows << "," << result.m_ImageColumns << "," << result.m_ImageSlices); MITK_TEST_OUTPUT( << "Read spacing parameters (x,y,z): " << result.m_Spacing[0] << "," << result.m_Spacing[1] << "," << result.m_Spacing[2]); NodeList gaussians; testimage->GetChildren("gaussian", gaussians); MITK_TEST_CONDITION_REQUIRED( gaussians.size() >= 1, "At least one gaussian is defined" ) result.m_CenterX.resize(result.m_NumberOfGaussian); result.m_CenterY.resize(result.m_NumberOfGaussian); result.m_CenterZ.resize(result.m_NumberOfGaussian); result.m_SigmaX.resize(result.m_NumberOfGaussian); result.m_SigmaY.resize(result.m_NumberOfGaussian); result.m_SigmaZ.resize(result.m_NumberOfGaussian); result.m_Altitude.resize(result.m_NumberOfGaussian); for(int i = 0; i < result.m_NumberOfGaussian ; ++i) { itk::DOMNode* gaussian = gaussians[i]; result.m_CenterX[i] = GetIntegerAttribute(gaussian, "centerIndexX"); result.m_CenterY[i] = GetIntegerAttribute(gaussian, "centerIndexY"); result.m_CenterZ[i] = GetIntegerAttribute(gaussian, "centerIndexZ"); result.m_SigmaX[i] = GetIntegerAttribute(gaussian, "deviationX"); result.m_SigmaY[i] = GetIntegerAttribute(gaussian, "deviationY"); result.m_SigmaZ[i] = GetIntegerAttribute(gaussian, "deviationZ"); result.m_Altitude[i] = GetIntegerAttribute(gaussian, "altitude"); } // read ROI parameters, fill result structure NodeList rois; domRoot->GetChildren("roi", rois); MITK_TEST_CONDITION_REQUIRED( rois.size() == 1, "One ROI defined" ) itk::DOMNode* roi = rois[0]; result.m_MaxSizeX = GetIntegerAttribute(roi, "maximumSizeX"); result.m_MinSizeX = GetIntegerAttribute(roi, "minimumSizeX"); result.m_MaxSizeY = GetIntegerAttribute(roi, "maximumSizeY"); result.m_MinSizeY = GetIntegerAttribute(roi, "minimumSizeY"); result.m_MaxSizeZ = GetIntegerAttribute(roi, "maximumSizeZ"); result.m_MinSizeZ = GetIntegerAttribute(roi, "minimumSizeZ"); result.m_hotspotRadiusInMM = GetDoubleAttribute(roi, "hotspotRadiusInMM"); // 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_HotspotMin = GetDoubleAttribute(statistic, "minimum"); result.m_HotspotMax = GetDoubleAttribute(statistic, "maximum"); result.m_HotspotMean = GetDoubleAttribute(statistic, "mean"); result.m_HotspotMinIndex.set_size(3); result.m_HotspotMinIndex[0] = GetIntegerAttribute(statistic, "minimumIndexX"); result.m_HotspotMinIndex[1] = GetIntegerAttribute(statistic, "minimumIndexY"); result.m_HotspotMinIndex[2] = GetIntegerAttribute(statistic, "minimumIndexZ"); result.m_HotspotMaxIndex.set_size(3); result.m_HotspotMaxIndex[0] = GetIntegerAttribute(statistic, "maximumIndexX"); result.m_HotspotMaxIndex[1] = GetIntegerAttribute(statistic, "maximumIndexY"); result.m_HotspotMaxIndex[2] = GetIntegerAttribute(statistic, "maximumIndexZ"); result.m_HotspotIndex.set_size(3); result.m_HotspotIndex[0] = GetIntegerAttribute(statistic, "hotspotIndexX"); result.m_HotspotIndex[1] = GetIntegerAttribute(statistic, "hotspotIndexY"); result.m_HotspotIndex[2] = GetIntegerAttribute(statistic, "hotspotIndexZ"); 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; 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]); } ImageType::SpacingType spacing; for(int i = 0; i < Dimension; ++i) spacing[i] = testParameters.m_Spacing[i]; gaussianGenerator->SetSize( size ); gaussianGenerator->SetSpacing( spacing ); 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 < Dimension; ++i) { start[i] = 0; spacing[i] = testParameters.m_Spacing[i]; } 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->SetSpacing(spacing); mask->SetRegions(region); mask->Allocate(); typedef itk::ImageRegionIteratorWithIndex MaskImageIteratorType; MaskImageIteratorType maskIt(mask, region); for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt) { MaskImageType::IndexType index = maskIt.GetIndex(); if((index[0] >= testParameters.m_MinSizeX && index[0] <= testParameters.m_MaxSizeX) && (index[1] >= testParameters.m_MinSizeY && index[1] <= testParameters.m_MaxSizeY) && (index[2] >= testParameters.m_MinSizeZ && index[2] <= testParameters.m_MaxSizeZ)) maskIt.Set(1); else maskIt.Set(0); } mitk::Image::Pointer mitkMaskImage; mitk::CastToMitkImage(mask, mitkMaskImage); statisticsCalculator->SetImage(image); statisticsCalculator->SetImageMask(mitkMaskImage); statisticsCalculator->SetMaskingModeToImage(); - statisticsCalculator->SetHotspotRadius(testParameters.m_hotspotRadiusInMM); + statisticsCalculator->SetHotspotRadiusInMM(testParameters.m_hotspotRadiusInMM); statisticsCalculator->SetCalculateHotspot(true); 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 eps = 0.001; // float comparisons, allow tiny differences - MITK_TEST_CONDITION( ::fabs(testParameters.m_HotspotMean - statistics.HotspotMean) < eps, "Mean value of hotspot in XML-File: " << testParameters.m_HotspotMean << " (Mean value of hotspot calculated in mitkImageStatisticsCalculator: " << statistics.HotspotMean << ")" ); - MITK_TEST_CONDITION( ::fabs(testParameters.m_HotspotMax- statistics.HotspotMax) < eps, "Maximum of hotspot in XML-File: " << testParameters.m_HotspotMax << " (Maximum of hotspot calculated in mitkImageStatisticsCalculator: " << statistics.HotspotMax << ")" ); - MITK_TEST_CONDITION( ::fabs(testParameters.m_HotspotMin - statistics.HotspotMin) < eps, "Minimum of hotspot in XML-File: " << testParameters.m_HotspotMin << " (Minimum of hotspot calculated in mitkImageStatisticsCalculator: " << statistics.HotspotMin << ")" ); + MITK_TEST_CONDITION( ::fabs(testParameters.m_HotspotMean - statistics.GetHotspotMean() ) < eps, "Mean value of hotspot in XML-File: " << testParameters.m_HotspotMean << " (Mean value of hotspot calculated in mitkImageStatisticsCalculator: " << statistics.GetHotspotMean() << ")" ); + MITK_TEST_CONDITION( ::fabs(testParameters.m_HotspotMax- statistics.GetHotspotMax() ) < eps, "Maximum of hotspot in XML-File: " << testParameters.m_HotspotMax << " (Maximum of hotspot calculated in mitkImageStatisticsCalculator: " << statistics.GetHotspotMax() << ")" ); + MITK_TEST_CONDITION( ::fabs(testParameters.m_HotspotMin - statistics.GetHotspotMin() ) < eps, "Minimum of hotspot in XML-File: " << testParameters.m_HotspotMin << " (Minimum of hotspot calculated in mitkImageStatisticsCalculator: " << statistics.GetHotspotMin() << ")" ); - MITK_TEST_CONDITION( statistics.HotspotIndex == testParameters.m_HotspotIndex, "Index of hotspot in XML-File: " << testParameters.m_HotspotIndex << " (Index of hotspot calculated in mitkImageStatisticsCalculator: " << statistics.HotspotIndex << ")" ); - MITK_TEST_CONDITION( statistics.HotspotMaxIndex == testParameters.m_HotspotMaxIndex, "Index of HotspotMaximum in XML-File: " << testParameters.m_HotspotMaxIndex << " (Index of HotspotMaximum calculated in mitkImageStatisticsCalculator: " << statistics.HotspotMaxIndex << ")" ); - MITK_TEST_CONDITION( statistics.HotspotMinIndex == testParameters.m_HotspotMinIndex, "Index of HotspotMinimum in XML-File " << testParameters.m_HotspotMinIndex << " (Index of HotspotMinimum calculated in mitkImageStatisticsCalculator: " << statistics.HotspotMinIndex << ")" ); + MITK_TEST_CONDITION( statistics.GetHotspotIndex() == testParameters.m_HotspotIndex, "Index of hotspot in XML-File: " << testParameters.m_HotspotIndex << " (Index of hotspot calculated in mitkImageStatisticsCalculator: " << statistics.GetHotspotIndex() << ")" ); + MITK_TEST_CONDITION( statistics.GetHotspotMaxIndex() == testParameters.m_HotspotMaxIndex, "Index of HotspotMaximum in XML-File: " << testParameters.m_HotspotMaxIndex << " (Index of HotspotMaximum calculated in mitkImageStatisticsCalculator: " << statistics.GetHotspotMaxIndex() << ")" ); + MITK_TEST_CONDITION( statistics.GetHotspotMinIndex() == testParameters.m_HotspotMinIndex, "Index of HotspotMinimum in XML-File " << testParameters.m_HotspotMinIndex << " (Index of HotspotMinimum calculated in mitkImageStatisticsCalculator: " << statistics.GetHotspotMinIndex() << ")" ); } }; #include /** \brief Verifies that TODO hotspot statistics part of ImageStatisticsCalculator. bla... */ int mitkImageStatisticsHotspotTest(int argc, char* argv[]) { MITK_TEST_BEGIN("mitkImageStatisticsHotspotTest") try { // 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" ); itk::TimeProbe clock; clock.Start(); // calculate statistics for this image (potentially use parameters for statistics ROI) mitk::ImageStatisticsCalculator::Statistics statistics = mitkImageStatisticsHotspotTestClass::CalculateStatistics(image, parameters); clock.Stop(); std::cout << "Statistics time consumed: " << clock.GetTotal() << std::endl; // compare statistics against stored expected values mitkImageStatisticsHotspotTestClass::ValidateStatistics(statistics, parameters); } catch (std::exception& e) { std::cout << "Error: " << e.what() << std::endl; } MITK_TEST_END() } diff --git a/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp b/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp index eac745b4b3..fa218be244 100644 --- a/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp +++ b/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp @@ -1,1792 +1,1848 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkImageStatisticsCalculator.h" #include "mitkImageAccessByItk.h" #include "mitkImageCast.h" #include "mitkExtractImageFilter.h" #include "mitkImageTimeSelector.h" #include "mitkITKImageImport.h" #include #include #include #include #include #include #include #include #include #include -#include #include #include - #include #include #include #include #include #include #include #include #include -#include -#include -#include #include -#include -#include #include -#include #include #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 - // TODO DM: sort includes, check if they are really needed 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), // TODO DM: check order of variable initialization + m_PlanarFigureCoordinate1 (0), // TODO DM: check order of variable initialization m_HotspotRadiusInMM(6.2035049089940), // radius of a 1cm3 sphere in mm m_CalculateHotspot(false), m_HotspotRadiusInMMChanged(false) { m_EmptyHistogram = HistogramType::New(); m_EmptyHistogram->SetMeasurementVectorSize(1); HistogramType::SizeType histogramSize(1); histogramSize.Fill( 256 ); m_EmptyHistogram->Initialize( histogramSize ); - m_EmptyStatistics.Reset(); -} + // m_EmptyStatistics.Reset(); +} ImageStatisticsCalculator::~ImageStatisticsCalculator() { } +ImageStatisticsCalculator::Statistics::Statistics() +: m_Label(0), + m_N(0), + m_Min(0.0), + m_Max(0.0), + m_Median(0.0), + m_Mean(0.0), + m_Sigma(0.0), + m_RMS(0.0), + m_MaxIndex(1), //TODO: Dimension anpassen + m_MinIndex(1), //TODO: Dimension anpassen + m_HotspotN(0), + m_HotspotMin(0.0), + m_HotspotMax(0.0), + m_HotspotMedian(0.0), + m_HotspotMean(0.0), + m_HotspotSigma(0.0), + m_HotspotRMS(0.0), + m_HotspotIndex(1), //TODO: Dimension anpassen + m_HotspotMaxIndex(1), //TODO: Dimension anpassen + m_HotspotMinIndex(1) //TODO: Dimension anpassen +{ +} + + 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::SetHotspotRadius(double value) +void ImageStatisticsCalculator::SetHotspotRadiusInMM(double value) { if ( m_HotspotRadiusInMM != value ) { m_HotspotRadiusInMM = value; if(m_CalculateHotspot) { m_HotspotRadiusInMMChanged = true; MITK_INFO <<"Hotspot radius changed, new convolution required"; } this->Modified(); } } -double ImageStatisticsCalculator::GetHotspotRadius() +double ImageStatisticsCalculator::GetHotspotRadiusInMM() { return m_HotspotRadiusInMM; } void ImageStatisticsCalculator::SetCalculateHotspot(bool on) { if ( m_CalculateHotspot != on ) { m_CalculateHotspot = on; m_HotspotRadiusInMMChanged = true; MITK_INFO <<"Hotspot calculation changed, new convolution required"; this->Modified(); } } bool ImageStatisticsCalculator::IsHotspotCalculated() { return m_CalculateHotspot; } 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_HotspotRadiusInMMChanged && ((m_MaskingMode != MASKING_MODE_NONE) || (imageMTime > m_Image->GetMTime() && !imageStatisticsCalculationTrigger)) && ((m_MaskingMode != MASKING_MODE_IMAGE) || (maskedImageMTime > m_ImageMask->GetMTime() && !maskedImageStatisticsCalculationTrigger)) && ((m_MaskingMode != MASKING_MODE_PLANARFIGURE) || (planarFigureMTime > m_PlanarFigure->GetMTime() && !planarFigureStatisticsCalculationTrigger)) ) { // Statistics is up to date! if ( m_MaskingModeChanged ) { m_MaskingModeChanged = false; } else { return false; } } // Reset state changed flag m_MaskingModeChanged = false; m_IgnorePixelValueChanged = false; // Depending on masking mode, extract and/or generate the required image // and mask data from the user input this->ExtractImageAndMask( timeStep ); StatisticsContainer *statisticsContainer; HistogramContainer *histogramContainer; switch ( m_MaskingMode ) { case MASKING_MODE_NONE: default: if(!m_DoIgnorePixelValue) { statisticsContainer = &m_ImageStatisticsVector[timeStep]; histogramContainer = &m_ImageHistogramVector[timeStep]; m_ImageStatisticsTimeStampVector[timeStep].Modified(); m_ImageStatisticsCalculationTriggerVector[timeStep] = false; } else { statisticsContainer = &m_MaskedImageStatisticsVector[timeStep]; histogramContainer = &m_MaskedImageHistogramVector[timeStep]; m_MaskedImageStatisticsTimeStampVector[timeStep].Modified(); m_MaskedImageStatisticsCalculationTriggerVector[timeStep] = false; } break; case MASKING_MODE_IMAGE: statisticsContainer = &m_MaskedImageStatisticsVector[timeStep]; histogramContainer = &m_MaskedImageHistogramVector[timeStep]; m_MaskedImageStatisticsTimeStampVector[timeStep].Modified(); m_MaskedImageStatisticsCalculationTriggerVector[timeStep] = false; break; case MASKING_MODE_PLANARFIGURE: statisticsContainer = &m_PlanarFigureStatisticsVector[timeStep]; histogramContainer = &m_PlanarFigureHistogramVector[timeStep]; m_PlanarFigureStatisticsTimeStampVector[timeStep].Modified(); m_PlanarFigureStatisticsCalculationTriggerVector[timeStep] = false; break; } // Calculate statistics and histogram(s) if ( m_InternalImage->GetDimension() == 3 ) { if ( m_MaskingMode == MASKING_MODE_NONE && !m_DoIgnorePixelValue ) { AccessFixedDimensionByItk_2( m_InternalImage, InternalCalculateStatisticsUnmasked, 3, statisticsContainer, histogramContainer ); } else { AccessFixedDimensionByItk_3( m_InternalImage, InternalCalculateStatisticsMasked, 3, m_InternalImageMask3D.GetPointer(), statisticsContainer, histogramContainer ); } } else if ( m_InternalImage->GetDimension() == 2 ) { if ( m_MaskingMode == MASKING_MODE_NONE && !m_DoIgnorePixelValue ) { AccessFixedDimensionByItk_2( m_InternalImage, InternalCalculateStatisticsUnmasked, 2, statisticsContainer, histogramContainer ); } else { AccessFixedDimensionByItk_3( m_InternalImage, InternalCalculateStatisticsMasked, 2, m_InternalImageMask2D.GetPointer(), statisticsContainer, histogramContainer ); } } else { MITK_ERROR << "ImageStatistics: Image dimension not supported!"; } // Release unused image smart pointers to free memory m_InternalImage = mitk::Image::ConstPointer(); m_InternalImageMask3D = MaskImage3DType::Pointer(); m_InternalImageMask2D = MaskImage2DType::Pointer(); return true; } 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::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() ); } } MITK_DEBUG << "Update of convolution image required?\n m_CalculateHotspot: " << m_CalculateHotspot << "\n m_HotspotSearchConvolutionImage: " << (void*) m_HotspotSearchConvolutionImage.GetPointer() << "\n m_ImageStatisticsCalculationTriggerVector["<GetMTime() << "\n ImageStatistics::MTime: " << this->GetMTime() << "\n m_Image->GetMTime(): " << m_Image->GetMTime(); if( m_CalculateHotspot && ( m_HotspotSearchConvolutionImage.IsNull() || m_Image->GetMTime() > this->GetMTime() // TODO check when m_InternalImage 'really' changes; depends on timeStep || m_HotspotRadiusInMMChanged == true ) ) { MITK_DEBUG <<" --> Update required."; if ( m_InternalImage->GetDimension() == 3 ) { AccessFixedDimensionByItk( m_InternalImage, InternalUpdateConvolutionImage, 3 ); } else if ( m_InternalImage->GetDimension() == 2 ) { AccessFixedDimensionByItk( m_InternalImage, InternalUpdateConvolutionImage, 2 ); } } else { MITK_DEBUG <<" --> Update required."; } } 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 statistics; //statistics.Reset(); + statistics.SetLabel(1); + statistics.SetN(image->GetBufferedRegion().GetNumberOfPixels()); + statistics.SetMin(statisticsFilter->GetMinimum()); + statistics.SetMax(statisticsFilter->GetMaximum()); + statistics.SetMean(statisticsFilter->GetMean()); + statistics.SetMedian(0.0); + statistics.SetSigma(statisticsFilter->GetSigma()); + statistics.SetRMS(sqrt( statistics.GetMean() * statistics.GetMean() + statistics.GetSigma() * statistics.GetSigma() )); + + statistics.GetMinIndex().set_size(image->GetImageDimension()); + statistics.GetMaxIndex().set_size(image->GetImageDimension()); - statistics.MinIndex.set_size(image->GetImageDimension()); - statistics.MaxIndex.set_size(image->GetImageDimension()); - for (unsigned int i=0; i tmpMaxIndex; + vnl_vector tmpMinIndex; + + for (unsigned int i=0; iGetIndexOfMaximum()[i]; - statistics.MinIndex[i] = minMaxFilter->GetIndexOfMinimum()[i]; + tmpMaxIndex[i] = minMaxFilter->GetIndexOfMaximum()[i]; + tmpMinIndex[i] = minMaxFilter->GetIndexOfMinimum()[i]; + } + + statistics.SetMinIndex(tmpMaxIndex); + statistics.SetMinIndex(tmpMinIndex); + + if( IsHotspotCalculated() ) + { + typedef itk::Image< unsigned short, VImageDimension > MaskImageType; + typename MaskImageType::Pointer maskImageITK = MaskImageType::New(); + + maskImageITK->SetRegions ( image->GetLargestPossibleRegion() ); + maskImageITK->Allocate(); + typedef itk::ImageRegionIteratorWithIndex< MaskImageType > MaskImageIteratorType; + + MaskImageIteratorType maskIt(maskImageITK, image->GetLargestPossibleRegion()); + + for(maskIt.GoToBegin();!maskIt.IsAtEnd(); ++maskIt) + maskIt.Set(1); + + Image::Pointer maskImageMITK = ImportItkImage( maskImageITK ); + + ImageStatisticsCalculator::Pointer hotspotCalculator = ImageStatisticsCalculator::New(); + hotspotCalculator->SetImage( this->m_Image ); + hotspotCalculator->SetMaskingModeToImage(); + hotspotCalculator->SetImageMask( maskImageMITK ); + hotspotCalculator->SetCalculateHotspot( true ); + hotspotCalculator->SetHotspotRadiusInMM( GetHotspotRadiusInMM() ); + hotspotCalculator->ComputeStatistics(0); // TODO: timestep + + Statistics hotspotStatistics = hotspotCalculator->GetStatistics(); + ImageExtrema hotspotExtrema = CalculateExtremaWorld(image, maskImageITK.GetPointer()); + + statistics.SetHotspotN(hotspotStatistics.GetHotspotN()); + statistics.SetHotspotMax (hotspotStatistics.GetHotspotMax()); + statistics.SetHotspotMin(hotspotStatistics.GetHotspotMin()); + statistics.SetHotspotMean(hotspotStatistics.GetHotspotMean()); + statistics.SetHotspotMedian(0.0); + statistics.SetHotspotSigma (hotspotStatistics.GetHotspotSigma()); + statistics.SetHotspotRMS (sqrt (statistics.GetHotspotMean() * statistics.GetHotspotMean() + + statistics.GetHotspotSigma() * statistics.GetHotspotSigma() )); + statistics.SetHotspotMaxIndex(hotspotStatistics.GetHotspotMaxIndex()); + statistics.SetHotspotMinIndex(hotspotStatistics.GetHotspotMinIndex()); + statistics.SetHotspotIndex(hotspotExtrema.MaxIndex); } 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->SetHistogramMin( statistics.GetMin() ); + histogramGenerator->SetHistogramMax( statistics.GetMax() ); histogramGenerator->Compute(); // TODO DM: add hotspot search here! 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; } } MITK_INFO <<"Mask has pixels? " << maskNonEmpty; if ( maskNonEmpty ) { Statistics statistics; std::list< int >::iterator it; for ( it = relevantLabels.begin(), i = 0; it != relevantLabels.end(); ++it, ++i ) { histogramContainer->push_back( HistogramType::ConstPointer( labelStatisticsFilter->GetHistogram( (*it) ) ) ); - statistics.Label = (*it); - statistics.N = labelStatisticsFilter->GetCount( *it ); - statistics.Min = labelStatisticsFilter->GetMinimum( *it ); - statistics.Max = labelStatisticsFilter->GetMaximum( *it ); - statistics.Mean = labelStatisticsFilter->GetMean( *it ); - statistics.Median = labelStatisticsFilter->GetMedian( *it ); - statistics.Sigma = labelStatisticsFilter->GetSigma( *it ); - statistics.RMS = sqrt( statistics.Mean * statistics.Mean - + statistics.Sigma * statistics.Sigma ); + statistics.SetLabel (*it); + statistics.SetN(labelStatisticsFilter->GetCount( *it )); + statistics.SetMin(labelStatisticsFilter->GetMinimum( *it )); + statistics.SetMax(labelStatisticsFilter->GetMaximum( *it )); + statistics.SetMean(labelStatisticsFilter->GetMean( *it )); + statistics.SetMedian(labelStatisticsFilter->GetMedian( *it )); + statistics.SetSigma(labelStatisticsFilter->GetSigma( *it )); + statistics.SetRMS(sqrt( statistics.GetMean() * statistics.GetMean() + + statistics.GetSigma() * statistics.GetSigma() )); // restrict image to mask area for min/max index calculation typedef itk::MaskImageFilter< ImageType, MaskImageType, ImageType > MaskImageFilterType; typename MaskImageFilterType::Pointer masker = MaskImageFilterType::New(); - masker->SetOutsideValue( (statistics.Min+statistics.Max)/2 ); + masker->SetOutsideValue( (statistics.GetMin()+statistics.GetMax())/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. + + vnl_vector tmpMaxIndex; + vnl_vector tmpMinIndex; + tmpMaxIndex.set_size(m_Image->GetDimension()); + tmpMinIndex.set_size(m_Image->GetDimension()); + 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; + tmpMaxIndex[m_PlanarFigureCoordinate0] = tempMaxIndex[0]; + tmpMaxIndex[m_PlanarFigureCoordinate1] = tempMaxIndex[1]; + tmpMaxIndex[m_PlanarFigureAxis] = m_PlanarFigureSlice; + + tmpMinIndex[m_PlanarFigureCoordinate0] = tempMinIndex[0] ; + tmpMinIndex[m_PlanarFigureCoordinate1] = tempMinIndex[1]; + tmpMinIndex[m_PlanarFigureAxis] = m_PlanarFigureSlice; } else { - for (unsigned int i = 0; ipush_back( statistics ); } else { histogramContainer->push_back( HistogramType::ConstPointer( m_EmptyHistogram ) ); statisticsContainer->push_back( Statistics() ); // TODO DM: this is uninitialized! (refactor into real class!) } } // TODO DM: needs to be modified to calculate a specific or multiple(!) labels template ImageStatisticsCalculator::ImageExtrema ImageStatisticsCalculator::CalculateExtremaWorld( const itk::Image *inputImage, itk::Image *maskImage) { // TODO iterate onlye the region of the mask (in the inputImage), which is usually smaller (not sure if this is possible) typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< unsigned short, VImageDimension > MaskImageType; typedef itk::ImageRegionConstIteratorWithIndex MaskImageIteratorType; typedef itk::ImageRegionConstIteratorWithIndex InputImageIndexIteratorType; - MaskImageIteratorType maskIt(maskImage, maskImage->GetLargestPossibleRegion()); // TODO DM: we should use the same regions here + MaskImageIteratorType maskIt(maskImage, inputImage->GetLargestPossibleRegion()); InputImageIndexIteratorType imageIndexIt(inputImage, inputImage->GetLargestPossibleRegion()); - float maxValue = itk::NumericTraits::min(); // TODO DM: I DID correct this before: use named functions instead of using - + float maxValue = itk::NumericTraits::min(); float minValue = itk::NumericTraits::max(); typename ImageType::IndexType maxIndex; typename ImageType::IndexType minIndex; typename ImageType::IndexType imageIndex; typename ImageType::PointType worldPosition; typename ImageType::IndexType maskIndex; for(imageIndexIt.GoToBegin(); !imageIndexIt.IsAtEnd(); ++imageIndexIt) { imageIndex = imageIndexIt.GetIndex(); inputImage->TransformIndexToPhysicalPoint(imageIndex, worldPosition); maskImage->TransformPhysicalPointToIndex(worldPosition, maskIndex); maskIt.SetIndex( maskIndex ); if(maskIt.Get() > 0) { double value = imageIndexIt.Get(); //Calculate minimum, maximum and corresponding index-values if( value > maxValue ) { maxIndex = imageIndexIt.GetIndex(); maxValue = value; } if(value < minValue ) { minIndex = imageIndexIt.GetIndex(); minValue = value; } } } ImageExtrema minMax; minMax.MinIndex.set_size(inputImage->GetImageDimension()); minMax.MaxIndex.set_size(inputImage->GetImageDimension()); for(unsigned int i = 0; i < minMax.MaxIndex.size(); ++i) minMax.MaxIndex[i] = maxIndex[i]; for(unsigned int i = 0; i < minMax.MinIndex.size(); ++i) minMax.MinIndex[i] = minIndex[i]; minMax.Max = maxValue; minMax.Min = minValue; return minMax; } - -// TODO DM: needs to be modified to calculate a specific or multiple(!) labels -template -ImageStatisticsCalculator::ImageExtrema ImageStatisticsCalculator::CalculateExtrema( - const itk::Image *inputImage, - itk::Image *maskImage) -{ - typedef itk::Image< TPixel, VImageDimension > ImageType; - typedef itk::Image< unsigned short, VImageDimension > MaskImageType; - - typedef itk::ImageRegionConstIterator MaskImageIteratorType; - typedef itk::ImageRegionConstIteratorWithIndex InputImageIndexIteratorType; - - MaskImageIteratorType maskIt(maskImage, maskImage->GetLargestPossibleRegion()); // TODO DM: we should use the same regions here - InputImageIndexIteratorType imageIndexIt(inputImage, inputImage->GetLargestPossibleRegion()); - - float maxValue = itk::NumericTraits::min(); // TODO DM: I DID correct this before: use named functions instead of using - - float minValue = itk::NumericTraits::max(); - - typename ImageType::IndexType maxIndex; - typename ImageType::IndexType minIndex; - - for(maskIt.GoToBegin(), imageIndexIt.GoToBegin(); - !maskIt.IsAtEnd() && !imageIndexIt.IsAtEnd(); - ++maskIt, ++imageIndexIt) - { - if(maskIt.Get() > itk::NumericTraits::Zero) // TODO DM: this is where multiple mask values could be used - { - double value = imageIndexIt.Get(); - - //Calculate minimum, maximum and corresponding index-values - if( value > maxValue ) - { - maxIndex = imageIndexIt.GetIndex(); - maxValue = value; - } - - if(value < minValue ) - { - minIndex = imageIndexIt.GetIndex(); - minValue = value; - } - } - } - - ImageExtrema minMax; - - minMax.MinIndex.set_size(inputImage->GetImageDimension()); - minMax.MaxIndex.set_size(inputImage->GetImageDimension()); - - for(unsigned int i = 0; i < minMax.MaxIndex.size(); ++i) - minMax.MaxIndex[i] = maxIndex[i]; - - for(unsigned int i = 0; i < minMax.MinIndex.size(); ++i) - minMax.MinIndex[i] = minIndex[i]; - - - minMax.Max = maxValue; - minMax.Min = minValue; - - return minMax; -} - template itk::Size ImageStatisticsCalculator ::CalculateConvolutionKernelSize(double spacing[VImageDimension], double radiusInMM) { typedef itk::Image< float, VImageDimension > KernelImageType; typedef typename KernelImageType::SizeType SizeType; SizeType maskSize; for(unsigned int i = 0; i < VImageDimension; ++i) { maskSize[i] = ::ceil( 2.0 * radiusInMM / spacing[i] ); // We always want an uneven size to have a clear center point in the convolution mask if(maskSize[i] % 2 == 0 ) { ++maskSize[i]; } } return maskSize; } template itk::SmartPointer< itk::Image > ImageStatisticsCalculator ::GenerateHotspotSearchConvolutionKernel(double spacing[VImageDimension], double radiusInMM) { std::stringstream ss; for (unsigned int i = 0; i < VImageDimension; ++i) { ss << spacing[i]; if (i < VImageDimension -1) ss << ","; } MITK_DEBUG << "Update convolution kernel for spacing (" << ss.str() << ") and radius " << radiusInMM << "mm"; double radiusInMMSquared = radiusInMM * radiusInMM; typedef itk::Image< float, VImageDimension > KernelImageType; typename KernelImageType::Pointer convolutionKernel = KernelImageType::New(); // Calculate size and allocate mask image typedef typename KernelImageType::IndexType IndexType; IndexType maskIndex; maskIndex.Fill(0); typedef typename KernelImageType::SizeType SizeType; SizeType maskSize = this->CalculateConvolutionKernelSize(spacing, radiusInMM); Point3D convolutionMaskCenter; convolutionMaskCenter.Fill(0.0); for(unsigned int i = 0; i < VImageDimension; ++i) { convolutionMaskCenter[i] = 0.5 * (double)(maskSize[i]-1); } typedef typename KernelImageType::RegionType RegionType; RegionType maskRegion; maskRegion.SetSize(maskSize); maskRegion.SetIndex(maskIndex); convolutionKernel->SetRegions(maskRegion); convolutionKernel->SetSpacing(spacing); convolutionKernel->Allocate(); // Fill mask image values by subsampling the image grid typedef itk::ImageRegionIteratorWithIndex MaskIteratorType; MaskIteratorType maskIt(convolutionKernel,maskRegion); int numberOfSubVoxelsPerDimension = 2; // per dimension! int numberOfSubVoxels = ::pow( static_cast(numberOfSubVoxelsPerDimension), static_cast(VImageDimension) ); - double subVoxelSize = 1.0 / (double)numberOfSubVoxelsPerDimension; //(double)numberOfSubVoxels; + double subVoxelSize = 1.0 / (double)numberOfSubVoxelsPerDimension; double valueOfOneSubVoxel = 1.0 / (double)numberOfSubVoxels; double maskValue = 0.0; Point3D subVoxelPosition; double distanceSquared = 0.0; typedef itk::ContinuousIndex ContinuousIndexType; for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt) { ContinuousIndexType indexPoint(maskIt.GetIndex()); Point3D voxelPosition; for (unsigned int dimension = 0; dimension < VImageDimension; ++dimension) { voxelPosition[dimension] = indexPoint[dimension]; } // TODO DM: regard all dimensions, including z! (former code used only x/y) // TODO DM: generalize: not x, y, z but a for loop over dimension // TODO DM: this could be done by calling a recursive method, handing over the "remaining number of dimensions to iterate" maskValue = 0.0; Vector3D subVoxelOffset; subVoxelOffset.Fill(0.0); // iterate sub-voxels by iterating all possible offsets for (subVoxelOffset[0] = -0.5 + subVoxelSize / 2.0; subVoxelOffset[0] < +0.5; subVoxelOffset[0] += subVoxelSize) { for (subVoxelOffset[1] = -0.5 + subVoxelSize / 2.0; subVoxelOffset[1] < +0.5; subVoxelOffset[1] += subVoxelSize) { for (subVoxelOffset[2] = -0.5 + subVoxelSize / 2.0; subVoxelOffset[2] < +0.5; subVoxelOffset[2] += subVoxelSize) { subVoxelPosition = voxelPosition + subVoxelOffset; // TODO DM: this COULD be integrated into the for-loops if neccessary (add voxelPosition to initializer and end condition) //if ( subVoxelPosition.EuclideanDistanceTo( convolutionMaskCenter ) < radiusInMM ) // TODO DM: this is too much matrix operations, we calculate ourselves, check if this time is relevant - distanceSquared = (subVoxelPosition[0]-convolutionMaskCenter[0]) / spacing[0] * (subVoxelPosition[0]-convolutionMaskCenter[0]) / spacing[0] + // maskValue += valueOfOneSubVoxel; + distanceSquared = (subVoxelPosition[0]-convolutionMaskCenter[0]) / spacing[0] * (subVoxelPosition[0]-convolutionMaskCenter[0]) / spacing[0] + (subVoxelPosition[1]-convolutionMaskCenter[1]) / spacing[1] * (subVoxelPosition[1]-convolutionMaskCenter[1]) / spacing[1] + (subVoxelPosition[2]-convolutionMaskCenter[2]) / spacing[2] * (subVoxelPosition[2]-convolutionMaskCenter[2]) / spacing[2]; if (distanceSquared <= radiusInMMSquared) { maskValue += valueOfOneSubVoxel; } } } } maskIt.Set( maskValue ); } return convolutionKernel; } template void ImageStatisticsCalculator::InternalUpdateConvolutionImage( itk::Image* inputImage ) { double spacing[VImageDimension]; for (unsigned int dimension = 0; dimension < VImageDimension; ++dimension) { spacing[dimension] = inputImage->GetSpacing()[dimension]; } // update convolution kernel typedef itk::Image< float, VImageDimension > KernelImageType; typename KernelImageType::Pointer convolutionKernel = this->GenerateHotspotSearchConvolutionKernel(spacing, m_HotspotRadiusInMM); // TODO: if GenerateHotspotSearchConvolutionKernel() consumes relevant time, we need an additional condition // update convolution image typedef itk::Image< TPixel, VImageDimension > InputImageType; typedef itk::Image< TPixel, VImageDimension > ConvolutionImageType; typedef itk::FFTConvolutionImageFilter ConvolutionFilterType; typedef itk::ConstantBoundaryCondition BoundaryConditionType; BoundaryConditionType boundaryCondition; boundaryCondition.SetConstant(0.0); typename ConvolutionFilterType::Pointer convolutionFilter = ConvolutionFilterType::New(); convolutionFilter->SetBoundaryCondition(&boundaryCondition); convolutionFilter->SetInput(inputImage); convolutionFilter->SetKernelImage(convolutionKernel); convolutionFilter->SetNormalize(true); //MITK_INFO << "Update Convolution image for hotspot search"; convolutionFilter->UpdateLargestPossibleRegion(); // TODO check if we could benefit from restricting this for a region typename ConvolutionImageType::Pointer convolutionImage = convolutionFilter->GetOutput(); convolutionImage->SetSpacing( inputImage->GetSpacing() ); // TODO: only workaround because convolution filter seems to ignore spacing of input image m_HotspotSearchConvolutionImage = convolutionImage.GetPointer(); m_HotspotRadiusInMMChanged = false; } template < typename TPixel, unsigned int VImageDimension> void ImageStatisticsCalculator ::FillHotspotMaskPixels( itk::Image* maskImage, itk::Point sphereCenter, double sphereRadiusInMM) { typedef itk::Image< TPixel, VImageDimension > MaskImageType; typedef itk::ImageRegionIteratorWithIndex MaskImageIteratorType; MaskImageIteratorType maskIt(maskImage, maskImage->GetLargestPossibleRegion()); // TODO DM: we should use the same regions here typename MaskImageType::IndexType maskIndex; typename MaskImageType::PointType worldPosition; for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt) { maskIndex = maskIt.GetIndex(); maskImage->TransformIndexToPhysicalPoint(maskIndex, worldPosition); maskIt.Set( worldPosition.EuclideanDistanceTo(sphereCenter) <= sphereRadiusInMM ? 1 : 0 ); } } template < typename TPixel, unsigned int VImageDimension> ImageStatisticsCalculator::Statistics ImageStatisticsCalculator::CalculateHotspotStatistics( const itk::Image* inputImage, itk::Image* maskImage, double radiusInMM) { //MITK_INFO << "CalculateHotspotStatistics()"; // get convolution image (updated in InternalUpdateConvolutionImage()) typedef itk::Image< TPixel, VImageDimension > ConvolutionImageType; typedef itk::Image< float, VImageDimension > KernelImageType; typedef itk::Image< unsigned short, VImageDimension > MaskImageType; typename ConvolutionImageType::Pointer convolutionImage = dynamic_cast(m_HotspotSearchConvolutionImage.GetPointer()); if (convolutionImage.IsNull()) { MITK_ERROR << "Empty convolution image in CalculateHotspotStatistics(). We should never reach this state (logic error)."; throw std::logic_error("Empty convolution image in CalculateHotspotStatistics()"); } // find maximum in convolution image, given the current mask (this might change over time, while we assume the input fixed (TODO wrong assumption)) ImageExtrema pi = CalculateExtremaWorld(convolutionImage.GetPointer(), maskImage); /* MITK_INFO <<"Initial search for hotspot: " "\n Index: " << pi.MaxIndex[0] << "," << pi.MaxIndex[1] << "," << pi.MaxIndex[2] << "\n Value(hotspot): " << pi.Max<< "\n Index(min): " << pi.MinIndex[0] << "," << pi.MinIndex[1] << "," << pi.MinIndex[2] << "\n Value(min): " << pi.Min;*/ // create mask corresponding to hotspot region // mask is defined on the inputImage grid and is // dimensioned as the hotspot convolution kernel (the sphere) double spacing[VImageDimension]; for (unsigned int dimension = 0; dimension < VImageDimension; ++dimension) { spacing[dimension] = inputImage->GetSpacing()[dimension]; } typedef typename ConvolutionImageType::SizeType SizeType; SizeType maskSize = this->CalculateConvolutionKernelSize(spacing, m_HotspotRadiusInMM); typedef typename ConvolutionImageType::IndexType IndexType; IndexType maskIndex; maskIndex.Fill(0); //MITK_INFO << "Hotspot statistics mask started with size ["< shift -2 required if (maskIndex[dimension] < 0) { maskIndex[dimension] = 0; } /* TODO if (maskIndex[dimension] + maskSize[dimension] > inputImage->GetBufferedRegion().GetSize()[dimension] ) { maskSize[dimension] = inputImage->GetBufferedRegionLargestPossibleRegion().GetSize()[dimension] - maskIndex[dimension]; } */ } //MITK_INFO << "Hotspot statistics mask corrected as region of size ["<CopyInformation( inputImage ); // type not optimal, but image grid is good typedef typename ConvolutionImageType::RegionType RegionType; RegionType hotspotMaskRegion; IndexType mi; mi.Fill(0); hotspotMaskRegion.SetIndex( mi ); hotspotMaskRegion.SetSize( maskSize ); hotspotMaskITK->SetRegions( hotspotMaskRegion ); hotspotMaskITK->Allocate(); typename ConvolutionImageType::PointType maskOrigin; inputImage->TransformIndexToPhysicalPoint(maskIndex,maskOrigin); //MITK_INFO << "Mask origin at: " << maskOrigin; hotspotMaskITK->SetOrigin(maskOrigin); IndexType maskCenterIndex; for (unsigned int d =0; d< VImageDimension;++d) maskCenterIndex[d]=pi.MaxIndex[d]; typename ConvolutionImageType::PointType maskCenter; inputImage->TransformIndexToPhysicalPoint(maskCenterIndex,maskCenter); //MITK_INFO << "Mask center in input image: " << maskCenter; FillHotspotMaskPixels(hotspotMaskITK.GetPointer(), maskCenter, m_HotspotRadiusInMM); Image::Pointer hotspotMaskMITK = ImportItkImage( hotspotMaskITK ); Point3D maskCenterPosition = hotspotMaskMITK->GetGeometry()->GetCenter(); //MITK_INFO << "Mask center: " << maskCenterPosition; // use second instance of ImageStatisticsCalculator to calculate hotspot statistics ImageStatisticsCalculator::Pointer calculator = ImageStatisticsCalculator::New(); calculator->SetImage( this->m_Image ); calculator->SetMaskingModeToImage(); calculator->SetImageMask( hotspotMaskMITK ); calculator->SetCalculateHotspot( false ); calculator->ComputeStatistics(0); // TODO: timestep Statistics hotspotMaskStatistics = calculator->GetStatistics(0); Statistics hotspotStatistics; - hotspotStatistics.HotspotMin = hotspotMaskStatistics.Min; - hotspotStatistics.HotspotMinIndex = hotspotMaskStatistics.MinIndex; - hotspotStatistics.HotspotMax = hotspotMaskStatistics.Max; - hotspotStatistics.HotspotMaxIndex = hotspotMaskStatistics.MaxIndex; - hotspotStatistics.HotspotMean = hotspotMaskStatistics.Mean; - hotspotStatistics.HotspotMedian = hotspotMaskStatistics.Median; - hotspotStatistics.HotspotIndex = pi.MaxIndex; - hotspotStatistics.HotspotSigma = hotspotMaskStatistics.Sigma; - hotspotStatistics.HotspotRMS = sqrt( hotspotMaskStatistics.Mean * hotspotMaskStatistics.Mean - + hotspotMaskStatistics.Sigma * hotspotMaskStatistics.Sigma ); - hotspotStatistics.HotspotN = hotspotMaskStatistics.N; + hotspotStatistics.SetHotspotMin(hotspotMaskStatistics.GetMin()); + hotspotStatistics.SetHotspotMinIndex(hotspotMaskStatistics.GetMinIndex()); + hotspotStatistics.SetHotspotMax(hotspotMaskStatistics.GetMax()); + hotspotStatistics.SetHotspotMaxIndex(hotspotMaskStatistics.GetMaxIndex()); + hotspotStatistics.SetHotspotMean(hotspotMaskStatistics.GetMean()); + hotspotStatistics.SetHotspotMedian(hotspotMaskStatistics.GetMedian()); + hotspotStatistics.SetHotspotIndex(pi.MaxIndex); + hotspotStatistics.SetHotspotSigma(hotspotMaskStatistics.GetSigma()); + hotspotStatistics.SetHotspotRMS(sqrt( hotspotMaskStatistics.GetMean() * hotspotMaskStatistics.GetMean() + + hotspotMaskStatistics.GetSigma() * hotspotMaskStatistics.GetSigma() )); + hotspotStatistics.SetHotspotN(hotspotMaskStatistics.GetN()); /*MITK_INFO << "----- Hotspot search results:" "\n Index: " << hotspotStatistics.HotspotIndex[0] << "," << hotspotStatistics.HotspotIndex[1] << "," << hotspotStatistics.HotspotIndex[2] << "\n Value: " << hotspotStatistics.HotspotMean << "\n Max Index: " << hotspotStatistics.HotspotMaxIndex[0] << "," << hotspotStatistics.HotspotMaxIndex[1] << "," << hotspotStatistics.HotspotMaxIndex[2] << "\n Max Value: " << hotspotStatistics.HotspotMax << "\n Min Index: " << hotspotStatistics.HotspotMinIndex[0] << "," << hotspotStatistics.HotspotMinIndex[1] << "," << hotspotStatistics.HotspotMinIndex[2] << "\n Min Value: " << hotspotStatistics.HotspotMin;*/ return hotspotStatistics; } template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator::InternalCalculateMaskFromPlanarFigure( const itk::Image< TPixel, VImageDimension > *image, unsigned int axis ) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::CastImageFilter< ImageType, MaskImage2DType > CastFilterType; // Generate mask image as new image with same header as input image and // initialize with 1. typename CastFilterType::Pointer castFilter = CastFilterType::New(); castFilter->SetInput( image ); castFilter->Update(); castFilter->GetOutput()->FillBuffer( 1 ); // all PolylinePoints of the PlanarFigure are stored in a vtkPoints object. // These points are used by the vtkLassoStencilSource to create // a vtkImageStencil. const mitk::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() ); } +// Get-functions statistics +unsigned int ImageStatisticsCalculator::Statistics::GetN() const { return m_N; } +double ImageStatisticsCalculator::Statistics::GetMean() const { return m_Mean; } +double ImageStatisticsCalculator::Statistics::GetMin() const { return m_Min; } +double ImageStatisticsCalculator::Statistics::GetMax() const { return m_Max; } +double ImageStatisticsCalculator::Statistics::GetMedian() const { return m_Median; } +double ImageStatisticsCalculator::Statistics::GetVariance() const { return m_Variance; } +double ImageStatisticsCalculator::Statistics::GetSigma() const { return m_Sigma; } +double ImageStatisticsCalculator::Statistics::GetRMS() const { return m_RMS; } +vnl_vector ImageStatisticsCalculator::Statistics::GetMaxIndex() const { return m_MaxIndex; } +vnl_vector ImageStatisticsCalculator::Statistics::GetMinIndex() const { return m_MinIndex; } + +// Set-fucntions statistics +void ImageStatisticsCalculator::Statistics::SetLabel(unsigned int label) { m_Label = label; } +void ImageStatisticsCalculator::Statistics::SetN(unsigned int n) { m_N = n; } +void ImageStatisticsCalculator::Statistics::SetMean(double mean) { m_Mean = mean; } +void ImageStatisticsCalculator::Statistics::SetMin(double min) { m_Min = min; } +void ImageStatisticsCalculator::Statistics::SetMax(double max) { m_Max = max; } +void ImageStatisticsCalculator::Statistics::SetMedian(double median) { m_Median = median; } +void ImageStatisticsCalculator::Statistics::SetVariance(double variance) { m_Variance = variance; } +void ImageStatisticsCalculator::Statistics::SetSigma(double sigma) { m_Sigma = sigma; } +void ImageStatisticsCalculator::Statistics::SetRMS(double rms) { m_RMS = rms; } +void ImageStatisticsCalculator::Statistics::SetMaxIndex(vnl_vector maxIndex) { m_MaxIndex = maxIndex; } +void ImageStatisticsCalculator::Statistics::SetMinIndex(vnl_vector minIndex) { m_MinIndex = minIndex; } + +// Get-fucntions hotspot-statistics +unsigned int ImageStatisticsCalculator::Statistics::GetHotspotN() const { return m_HotspotN; } +double ImageStatisticsCalculator::Statistics::GetHotspotMean() const { return m_HotspotMean; } +double ImageStatisticsCalculator::Statistics::GetHotspotMin() const { return m_HotspotMin; } +double ImageStatisticsCalculator::Statistics::GetHotspotMax() const { return m_HotspotMax; } +double ImageStatisticsCalculator::Statistics::GetHotspotMedian() const { return m_HotspotMedian; } +double ImageStatisticsCalculator::Statistics::GetHotspotVariance() const { return m_HotspotVariance; } +double ImageStatisticsCalculator::Statistics::GetHotspotSigma() const { return m_HotspotSigma; } +double ImageStatisticsCalculator::Statistics::GetHotspotRMS() const { return m_HotspotRMS; } +vnl_vector ImageStatisticsCalculator::Statistics::GetHotspotMaxIndex() const { return m_HotspotMaxIndex; } +vnl_vector ImageStatisticsCalculator::Statistics::GetHotspotMinIndex() const { return m_HotspotMinIndex; } +vnl_vector ImageStatisticsCalculator::Statistics::GetHotspotIndex() const {return m_HotspotIndex;} + +// Set-functions hotspot-statistics +void ImageStatisticsCalculator::Statistics::SetHotspotN(unsigned int n) { m_HotspotN = n; } +void ImageStatisticsCalculator::Statistics::SetHotspotMean(double mean) { m_HotspotMean = mean; } +void ImageStatisticsCalculator::Statistics::SetHotspotMin(double min) { m_HotspotMin = min; } +void ImageStatisticsCalculator::Statistics::SetHotspotMax(double max) { m_HotspotMax = max; } +void ImageStatisticsCalculator::Statistics::SetHotspotMedian(double median) { m_HotspotMedian = median; } +void ImageStatisticsCalculator::Statistics::SetHotspotVariance(double variance) { m_HotspotVariance = variance; } +void ImageStatisticsCalculator::Statistics::SetHotspotSigma(double sigma) { m_HotspotSigma = sigma; } +void ImageStatisticsCalculator::Statistics::SetHotspotRMS(double rms) { m_HotspotRMS = rms; } +void ImageStatisticsCalculator::Statistics::SetHotspotMaxIndex(vnl_vector maxIndex) { m_HotspotMaxIndex = maxIndex; } +void ImageStatisticsCalculator::Statistics::SetHotspotMinIndex(vnl_vector minIndex) { m_HotspotMinIndex = minIndex; } +void ImageStatisticsCalculator::Statistics::SetHotspotIndex(vnl_vector index) { m_HotspotIndex = index; } } diff --git a/Modules/ImageStatistics/mitkImageStatisticsCalculator.h b/Modules/ImageStatistics/mitkImageStatisticsCalculator.h index deef51ec3c..73e7ca90ab 100644 --- a/Modules/ImageStatistics/mitkImageStatisticsCalculator.h +++ b/Modules/ImageStatistics/mitkImageStatisticsCalculator.h @@ -1,426 +1,524 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef mitkImageStatisticsCalculator_h #define mitkImageStatisticsCalculator_h #include "mitkImage.h" #include "mitkPlanarFigure.h" // TODO DM: why the ifndef? #ifndef __itkHistogram_h #include #endif #include #include #include "ImageStatisticsExports.h" namespace mitk { /** * \brief Class for calculating statistics and histogram for an (optionally * masked) image. * * Images can be masked by either a label image (of the same dimensions as * the original image) or by a closed mitk::PlanarFigure, e.g. a circle or * polygon. When masking with a planar figure, the slice corresponding to the * plane containing the figure is extracted and then clipped with contour * defined by the figure. Planar figures need to be aligned along the main axes * of the image (axial, sagittal, coronal). Planar figures on arbitrary * rotated planes are not supported. * * For each operating mode (no masking, masking by image, masking by planar * figure), the calculated statistics and histogram are cached so that, when * switching back and forth between operation modes without modifying mask or * image, the information doesn't need to be recalculated. * - * The class also has the possibility to calculate minimum, maximum, mean - * and their corresponding indicies in the hottest spot in a given ROI / VOI. - * The size of the hotspot is defined by a sphere with a radius specified by - * the user. This procedure is required for the calculation of SUV-statistics - * in PET-images for example. + * The class also has the possibility to calculate the location and separate + * statistics for a region called “hotspot”. The hotspot is a sphere of + * user-defined size and its location is chosen in a way that the average + * pixel value within the sphere is maximized. * * Note: currently time-resolved and multi-channel pictures are not properly * supported. */ + + /** + * \section HotspotStatistics_caption Calculation of hotspot statistics + * + * Since calculation of hotspot location and statistics is not + * straight-forward, the following paragraphs will describe it in more detail. + * + * Note: Calculation of hotspot statistics is optional. + * + * \subsection HotspotStatistics_description Hotspot Definition + * + * The “hotspot” of an image is motivated from PET readings. It is defined + * as a spherical region of fixed size which maximizes the average pixel value + * within the region. The following image illustrates the concept: the + * colored areas are different image intensities and the hotspot is located + * in the “hottest” region of the image. + * + * [image with pixelvalue scale] + * + * \subsection HotspotStatistics_calculation Hotspot Calculation + * + * Since only the size of the hotspot is known initially, we need to calculate + * two aspects (both implemented in ComputeHotspotStatistics): + * - the hotspot location + * - statistics of the pixels within the hotspot + * Finding the hotspot location requires to calculate the average value at each + * position. This is done by convolution of the image with a sperical kernel + * image which reflects partial volumes (important in the case of low-resolution + * PET images). + * + * Once the hotspot location is known, calculating the actual statistics is a + * simple task which is implemented in ...TODO... + * + * \b Step 1: Finding the hotspot by image convolution + * + * As described above, we use image convolution with a rasterized sphere to + * average the image at each position. To handle coarse resolutions, which would + * normally force us to decide for partially contained voxels whether to count + * them or not, we supersample the kernel image and use non-integer kernel values + * (see TODO-MethodName), which reflect the volume part that is contained in the + * sphere. For example, if three subvoxels are inside the sphere, the corresponding + * kernel voxel gets a value of 0.75 (3 out of 4 subvoxels, see 2D example below). + * + * [image with subsampled pixel] + * + * Convolution itself is done by means of the itk::FFTConvolutionImageFilter. + * To find the hotspot location, we simply iterate the averaged image and find a + * maximum location (see CalculateExtrema()). + * TODO? Discuss cases with multiple maxima!? + * -> "In case of images with multiple maxima the method returns value and corresponding + * index of the extrema that is found by the iterator first" + * + * \b Step 2: Computation of hotspot statistics + * + * Once the hotspot location is found, statistics for the region are calculated + * by simply iterating the input image and regarding all pixel centers inside the + * hotspot-sphere for statistics. + * + * \subsection HotspotStatistics_tests Tests + * + * To check the correctness of the hotspot calculation, a special class + * (mitkImageStatisticsHotspotTest) has been created, which generates images with + * known hotspot location and statistics. A number of unit tests use this class + * to first generate an image of known properites and then verify that + * ImageStatisticsCalculator is able to reproduce the known statistics. + * +*/ class ImageStatistics_EXPORT ImageStatisticsCalculator : public itk::Object { public: /** TODO DM: document */ enum { MASKING_MODE_NONE = 0, MASKING_MODE_IMAGE = 1, MASKING_MODE_PLANARFIGURE = 2 }; typedef itk::Statistics::Histogram HistogramType; typedef HistogramType::ConstIterator HistogramConstIteratorType; /** TODO DM: document */ - struct Statistics + class ImageStatistics_EXPORT Statistics { - int Label; - unsigned int N; //< number of voxels - double Min; //< mimimum value - double Max; //< maximum value - double Mean; //< mean value - double Median; //< median value - double Variance; //< variance of values // TODO DM: remove, was never filled with values ; check if any calling code within MITK used this member! - double Sigma; //< standard deviation of values (== square root of variance) - double RMS; //< root means square (TODO DM: check mesning) - - unsigned int HotspotN; - double HotspotMin; //< mimimum value inside hotspot - double HotspotMax; //< maximum value inside hotspot - double HotspotMean; //< mean value inside hotspot - double HotspotMedian; - double HotspotSigma; //< standard deviation of values inside hotspot - double HotspotRMS; - - vnl_vector< int > MinIndex; - vnl_vector< int > MaxIndex; - vnl_vector HotspotMaxIndex; - vnl_vector HotspotMinIndex; - vnl_vector HotspotIndex; //< TODO DM: couldn't this be named "hotspot index"? We need to clear naming of hotspotmean, hotspotpeakindex, and hotspotpeak - - Statistics() - { - this->Reset(); - } - - // TODO DM: make this struct a real class and put this into a constructor - void Reset() // TODO DM: move to .cpp file (mitk::ImageStatisticsCalculator::Statistics::Reset() {...}) - { - Label = 0; - N = 0; - Min = 0.0; - Max = 0.0; - Mean = 0.0; - Median = 0.0; - Variance = 0.0; - Sigma = 0.0; - RMS = 0.0; - HotspotMin = 0.0; - HotspotMax = 0.0; - HotspotMean = 0.0; - HotspotSigma = 0.0; // TODO DM: also reset index values! Check that everything is initialized - } + public: + + Statistics(); + + /** \brief Initialize every member of Statistics to zero. */ + //void Reset() + + unsigned int GetN() const; + double GetMin() const; + double GetMax() const; + double GetMean() const; + double GetMedian() const; + double GetVariance() const; + double GetSigma() const; + double GetRMS() const; + vnl_vector GetMaxIndex() const; + vnl_vector GetMinIndex() const; + + void SetLabel(unsigned int label); + void SetN(unsigned int n); + void SetMin(double min); + void SetMax(double max); + void SetMean(double mean); + void SetMedian(double median); + void SetVariance(double variance); + void SetSigma(double sigma); + void SetRMS(double rms); + void SetMaxIndex(vnl_vector maxIndex); + void SetMinIndex(vnl_vector minIndex); + + unsigned int GetHotspotN() const; + double GetHotspotMin() const; + double GetHotspotMax() const; + double GetHotspotMean() const; + double GetHotspotMedian() const; + double GetHotspotVariance() const; + double GetHotspotSigma() const; + double GetHotspotRMS() const; + vnl_vector GetHotspotMaxIndex() const; + vnl_vector GetHotspotMinIndex() const; + vnl_vector GetHotspotIndex() const; + + void SetHotspotLabel(unsigned int label); + void SetHotspotN(unsigned int n); + void SetHotspotMin(double min); + void SetHotspotMax(double max); + void SetHotspotMean(double mean); + void SetHotspotMedian(double median); + void SetHotspotVariance(double variance); + void SetHotspotSigma(double sigma); + void SetHotspotRMS(double rms); + void SetHotspotMaxIndex(vnl_vector hotspotMaxIndex); + void SetHotspotMinIndex(vnl_vector hotspotMinIndex); + void SetHotspotIndex(vnl_vector hotspotIndex); + + private: + + int m_Label; + unsigned int m_N; //< number of voxels + double m_Min; //< mimimum value + double m_Max; //< maximum value + double m_Mean; //< mean value + double m_Median; //< median value + double m_Variance; + double m_Sigma; //< standard deviation of values (== square root of variance) + double m_RMS; //< root mean square + + vnl_vector< int > m_MinIndex; //< index of minimum value + vnl_vector< int > m_MaxIndex; //< index of maximum value + + unsigned int m_HotspotN; //< number of voxels inside hotspot + double m_HotspotMin; //< mimimum value inside hotspot + double m_HotspotMax; //< maximum value inside hotspot + double m_HotspotMean; //< mean value of hotspot + double m_HotspotMedian; //< median value of hotspot + double m_HotspotVariance; + double m_HotspotSigma; //< standard deviation of values inside hotspot + double m_HotspotRMS; //< root mean square of hotspot + + vnl_vector m_HotspotMinIndex; //< index of minimum value inside hotspot + vnl_vector m_HotspotMaxIndex; //< index of maximum value inside hotspot + vnl_vector m_HotspotIndex; //< index of hotspot origin }; typedef std::vector< HistogramType::ConstPointer > HistogramContainer; typedef std::vector< Statistics > StatisticsContainer; mitkClassMacro( ImageStatisticsCalculator, itk::Object ); itkNewMacro( ImageStatisticsCalculator ); /** \brief Set image from which to compute statistics. */ void SetImage( const mitk::Image *image ); /** \brief Set image for masking. */ void SetImageMask( const mitk::Image *imageMask ); /** \brief Set planar figure for masking. */ void SetPlanarFigure( mitk::PlanarFigure *planarFigure ); /** \brief Set/Get operation mode for masking */ void SetMaskingMode( unsigned int mode ); /** \brief Set/Get operation mode for masking */ itkGetMacro( MaskingMode, unsigned int ); /** \brief Set/Get operation mode for masking */ void SetMaskingModeToNone(); /** \brief Set/Get operation mode for masking */ void SetMaskingModeToImage(); /** \brief Set/Get operation mode for masking */ void SetMaskingModeToPlanarFigure(); /** \brief Set a pixel value for pixels that will be ignored in the statistics */ void SetIgnorePixelValue(double value); /** \brief Get the pixel value for pixels that will be ignored in the statistics */ double GetIgnorePixelValue(); /** \brief Set whether a pixel value should be ignored in the statistics */ void SetDoIgnorePixelValue(bool doit); /** \brief Get whether a pixel value will be ignored in the statistics */ bool GetDoIgnorePixelValue(); /** \brief Sets the radius for the hotspot */ - void SetHotspotRadius (double hotspotRadiusInMM); // TODO in mm + void SetHotspotRadiusInMM (double hotspotRadiusInMM); /** \brief Returns the radius of the hotspot */ - double GetHotspotRadius(); // TODO in mm + double GetHotspotRadiusInMM(); /** \brief Sets whether the hotspot should be calculated */ void SetCalculateHotspot(bool calculateHotspot); /** \brief Returns true whether the hotspot should be calculated, otherwise false */ bool IsHotspotCalculated(); /** \brief Compute statistics (together with histogram) for the current * masking mode. * * Computation is not executed if statistics is already up to date. In this * case, false is returned; otherwise, true.*/ virtual bool ComputeStatistics( unsigned int timeStep = 0 ); /** \brief Retrieve the histogram depending on the current masking mode. * * \param label The label for which to retrieve the histogram in multi-label situations (ascending order). */ const HistogramType *GetHistogram( unsigned int timeStep = 0, unsigned int label = 0 ) const; /** \brief Retrieve the histogram depending on the current masking mode (for all image labels. */ const HistogramContainer &GetHistogramVector( unsigned int timeStep = 0 ) const; /** \brief Retrieve statistics depending on the current masking mode. * * \param label The label for which to retrieve the statistics in multi-label situations (ascending order). */ const Statistics &GetStatistics( unsigned int timeStep = 0, unsigned int label = 0 ) const; /** \brief Retrieve statistics depending on the current masking mode (for all image labels). */ const StatisticsContainer &GetStatisticsVector( unsigned int timeStep = 0 ) const; protected: typedef std::vector< HistogramContainer > HistogramVector; typedef std::vector< StatisticsContainer > StatisticsVector; typedef std::vector< itk::TimeStamp > TimeStampVectorType; typedef std::vector< bool > BoolVectorType; typedef itk::Image< unsigned short, 3 > MaskImage3DType; typedef itk::Image< unsigned short, 2 > MaskImage2DType; ImageStatisticsCalculator(); virtual ~ImageStatisticsCalculator(); /** \brief Depending on the masking mode, the image and mask from which to * calculate statistics is extracted from the original input image and mask * data. * * For example, a when using a PlanarFigure as mask, the 2D image slice * corresponding to the PlanarFigure will be extracted from the original * image. If masking is disabled, the original image is simply passed * through. */ void ExtractImageAndMask( unsigned int timeStep = 0 ); /** \brief If the passed vector matches any of the three principal axes * of the passed geometry, the ínteger value corresponding to the axis * is set and true is returned. */ bool GetPrincipalAxis( const Geometry3D *geometry, Vector3D vector, unsigned int &axis ); template < typename TPixel, unsigned int VImageDimension > void InternalCalculateStatisticsUnmasked( const itk::Image< TPixel, VImageDimension > *image, StatisticsContainer* statisticsContainer, HistogramContainer *histogramContainer ); template < typename TPixel, unsigned int VImageDimension > void InternalCalculateStatisticsMasked( const itk::Image< TPixel, VImageDimension > *image, itk::Image< unsigned short, VImageDimension > *maskImage, StatisticsContainer* statisticsContainer, HistogramContainer* histogramContainer ); template < typename TPixel, unsigned int VImageDimension > void InternalCalculateMaskFromPlanarFigure( const itk::Image< TPixel, VImageDimension > *image, unsigned int axis ); template < typename TPixel, unsigned int VImageDimension > void InternalMaskIgnoredPixels( const itk::Image< TPixel, VImageDimension > *image, itk::Image< unsigned short, VImageDimension > *maskImage ); struct ImageExtrema { double Max; double Min; vnl_vector MaxIndex; vnl_vector MinIndex; }; - /** \brief Calculates minimum, maximum, mean value and their + + /** \brief Calculates minimum, maximum, mean value and their * corresponding indices in a given ROI. As input the function * needs an image and a mask. It returns a ImageExtrema object. */ - template - ImageExtrema CalculateExtrema( - const itk::Image *inputImage, - itk::Image *maskImage); - template ImageExtrema CalculateExtremaWorld( const itk::Image *inputImage, itk::Image *maskImage); /** \brief Calculates the hotspot statistics within a given * ROI. As input the function needs an image, a mask which * represents the ROI and a radius which defines the size of * the sphere. The function returns a Statistics object. */ template < typename TPixel, unsigned int VImageDimension> Statistics CalculateHotspotStatistics( const itk::Image *inputImage, itk::Image *maskImage, double radiusInMM); /** Connection from ITK to VTK */ template void ConnectPipelines(ITK_Exporter exporter, vtkSmartPointer importer) { importer->SetUpdateInformationCallback(exporter->GetUpdateInformationCallback()); importer->SetPipelineModifiedCallback(exporter->GetPipelineModifiedCallback()); importer->SetWholeExtentCallback(exporter->GetWholeExtentCallback()); importer->SetSpacingCallback(exporter->GetSpacingCallback()); importer->SetOriginCallback(exporter->GetOriginCallback()); importer->SetScalarTypeCallback(exporter->GetScalarTypeCallback()); importer->SetNumberOfComponentsCallback(exporter->GetNumberOfComponentsCallback()); importer->SetPropagateUpdateExtentCallback(exporter->GetPropagateUpdateExtentCallback()); importer->SetUpdateDataCallback(exporter->GetUpdateDataCallback()); importer->SetDataExtentCallback(exporter->GetDataExtentCallback()); importer->SetBufferPointerCallback(exporter->GetBufferPointerCallback()); importer->SetCallbackUserData(exporter->GetCallbackUserData()); } /** Connection from VTK to ITK */ template void ConnectPipelines(vtkSmartPointer exporter, ITK_Importer importer) { importer->SetUpdateInformationCallback(exporter->GetUpdateInformationCallback()); importer->SetPipelineModifiedCallback(exporter->GetPipelineModifiedCallback()); importer->SetWholeExtentCallback(exporter->GetWholeExtentCallback()); importer->SetSpacingCallback(exporter->GetSpacingCallback()); importer->SetOriginCallback(exporter->GetOriginCallback()); importer->SetScalarTypeCallback(exporter->GetScalarTypeCallback()); importer->SetNumberOfComponentsCallback(exporter->GetNumberOfComponentsCallback()); importer->SetPropagateUpdateExtentCallback(exporter->GetPropagateUpdateExtentCallback()); importer->SetUpdateDataCallback(exporter->GetUpdateDataCallback()); importer->SetDataExtentCallback(exporter->GetDataExtentCallback()); importer->SetBufferPointerCallback(exporter->GetBufferPointerCallback()); importer->SetCallbackUserData(exporter->GetCallbackUserData()); } void UnmaskedStatisticsProgressUpdate(); void MaskedStatisticsProgressUpdate(); template itk::Size CalculateConvolutionKernelSize(double spacing[VImageDimension], double radiusInMM); template itk::SmartPointer< itk::Image > GenerateHotspotSearchConvolutionKernel(double spacing[VImageDimension], double radiusInMM); /** Uses members m_HotspotRadiusInMM */ template void InternalUpdateConvolutionImage( itk::Image* inputImage ); template < typename TPixel, unsigned int VImageDimension> void FillHotspotMaskPixels( itk::Image* maskImage, itk::Point sphereCenter, double sphereRadiusInMM); /** m_Image contains the input image (e.g. 2D, 3D, 3D+t)*/ mitk::Image::ConstPointer m_Image; mitk::Image::ConstPointer m_ImageMask; mitk::PlanarFigure::Pointer m_PlanarFigure; HistogramVector m_ImageHistogramVector; HistogramVector m_MaskedImageHistogramVector; HistogramVector m_PlanarFigureHistogramVector; HistogramType::Pointer m_EmptyHistogram; HistogramContainer m_EmptyHistogramContainer; StatisticsVector m_ImageStatisticsVector; StatisticsVector m_MaskedImageStatisticsVector; StatisticsVector m_PlanarFigureStatisticsVector; StatisticsVector m_MaskedImageHotspotStatisticsVector; Statistics m_EmptyStatistics; StatisticsContainer m_EmptyStatisticsContainer; unsigned int m_MaskingMode; bool m_MaskingModeChanged; /** m_InternalImage contains a image volume at one time step (e.g. 2D, 3D)*/ mitk::Image::ConstPointer m_InternalImage; MaskImage3DType::Pointer m_InternalImageMask3D; MaskImage2DType::Pointer m_InternalImageMask2D; TimeStampVectorType m_ImageStatisticsTimeStampVector; TimeStampVectorType m_MaskedImageStatisticsTimeStampVector; TimeStampVectorType m_PlanarFigureStatisticsTimeStampVector; BoolVectorType m_ImageStatisticsCalculationTriggerVector; BoolVectorType m_MaskedImageStatisticsCalculationTriggerVector; BoolVectorType m_PlanarFigureStatisticsCalculationTriggerVector; double m_IgnorePixelValue; bool m_DoIgnorePixelValue; bool m_IgnorePixelValueChanged; itk::Object::Pointer m_HotspotSearchConvolutionImage; // itk::Image - double m_HotspotRadiusInMM; - bool m_CalculateHotspot; - bool m_HotspotRadiusInMMChanged; - unsigned int m_PlanarFigureAxis; // Normal axis for PlanarFigure unsigned int m_PlanarFigureSlice; // Slice which contains PlanarFigure int m_PlanarFigureCoordinate0; // First plane-axis for PlanarFigure int m_PlanarFigureCoordinate1; // Second plane-axis for PlanarFigure + double m_HotspotRadiusInMM; + bool m_CalculateHotspot; + bool m_HotspotRadiusInMMChanged; + }; } // namespace #endif diff --git a/Modules/ImageStatistics/mitkPointSetDifferenceStatisticsCalculator.cpp b/Modules/ImageStatistics/mitkPointSetDifferenceStatisticsCalculator.cpp index 8d2173d8df..8e5d173598 100644 --- a/Modules/ImageStatistics/mitkPointSetDifferenceStatisticsCalculator.cpp +++ b/Modules/ImageStatistics/mitkPointSetDifferenceStatisticsCalculator.cpp @@ -1,214 +1,214 @@ /*=================================================================== 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 "mitkPointSetDifferenceStatisticsCalculator.h" mitk::PointSetDifferenceStatisticsCalculator::PointSetDifferenceStatisticsCalculator() : m_StatisticsCalculated(false) { m_PointSet1 = mitk::PointSet::New(); m_PointSet2 = mitk::PointSet::New(); - m_Statistics.Reset(); + //m_Statistics.Reset(); } mitk::PointSetDifferenceStatisticsCalculator::PointSetDifferenceStatisticsCalculator(mitk::PointSet::Pointer pSet1, mitk::PointSet::Pointer pSet2) { m_PointSet1 = pSet1; m_PointSet2 = pSet2; m_StatisticsCalculated = false; - m_Statistics.Reset(); + //m_Statistics.Reset(); } mitk::PointSetDifferenceStatisticsCalculator::~PointSetDifferenceStatisticsCalculator() { } void mitk::PointSetDifferenceStatisticsCalculator::SetPointSets(mitk::PointSet::Pointer pSet1, mitk::PointSet::Pointer pSet2) { if (pSet1.IsNotNull()) { m_PointSet1 = pSet1; } if (pSet2.IsNotNull()) { m_PointSet2 = pSet2; } m_StatisticsCalculated = false; - m_Statistics.Reset(); + //m_Statistics.Reset(); } std::vector mitk::PointSetDifferenceStatisticsCalculator::GetDifferences() { if (!m_StatisticsCalculated) { this->ComputeStatistics(); } return m_DifferencesVector; } std::vector mitk::PointSetDifferenceStatisticsCalculator::GetSquaredDifferences() { if (!m_StatisticsCalculated) { this->ComputeStatistics(); } return m_SquaredDifferencesVector; } double mitk::PointSetDifferenceStatisticsCalculator::GetMean() { if (!m_StatisticsCalculated) { this->ComputeStatistics(); } - return m_Statistics.Mean; + return m_Statistics.GetMean(); } double mitk::PointSetDifferenceStatisticsCalculator::GetSD() { if (!m_StatisticsCalculated) { this->ComputeStatistics(); } - return m_Statistics.Sigma; + return m_Statistics.GetSigma(); } double mitk::PointSetDifferenceStatisticsCalculator::GetVariance() { if (!m_StatisticsCalculated) { this->ComputeStatistics(); } - return m_Statistics.Variance; + return m_Statistics.GetVariance(); } double mitk::PointSetDifferenceStatisticsCalculator::GetRMS() { if (!m_StatisticsCalculated) { this->ComputeStatistics(); } - return m_Statistics.RMS; + return m_Statistics.GetRMS(); } double mitk::PointSetDifferenceStatisticsCalculator::GetMedian() { if (!m_StatisticsCalculated) { this->ComputeStatistics(); } - return m_Statistics.Median; + return m_Statistics.GetMedian(); } double mitk::PointSetDifferenceStatisticsCalculator::GetMax() { if (!m_StatisticsCalculated) { this->ComputeStatistics(); } - return m_Statistics.Max; + return m_Statistics.GetMax(); } double mitk::PointSetDifferenceStatisticsCalculator::GetMin() { if (!m_StatisticsCalculated) { this->ComputeStatistics(); } - return m_Statistics.Min; + return m_Statistics.GetMin(); } double mitk::PointSetDifferenceStatisticsCalculator::GetNumberOfPoints() { if (!m_StatisticsCalculated) { this->ComputeStatistics(); } - return m_Statistics.N; + return m_Statistics.GetN(); } void mitk::PointSetDifferenceStatisticsCalculator::ComputeStatistics() { if ((m_PointSet1.IsNull())||(m_PointSet2.IsNull())) { itkExceptionMacro("Point sets specified are not valid. Please specify correct Point sets"); } else if (m_PointSet1->GetSize()!=m_PointSet2->GetSize()) { itkExceptionMacro("PointSets are not equal. Please make sure that your PointSets have the same size and hold corresponding points."); } else if (m_PointSet1->GetSize()==0) { itkExceptionMacro("There are no points in the PointSets. Please make sure that the PointSets contain points"); } else { double mean = 0.0; double sd = 0.0; double rms= 0.0; std::vector differencesVector; mitk::Point3D point1; mitk::Point3D point2; int numberOfPoints = m_PointSet1->GetSize(); //Iterate over both pointsets in order to compare all points pair-wise mitk::PointSet::PointsIterator end = m_PointSet1->End(); for( mitk::PointSet::PointsIterator pointSetIterator = m_PointSet1->Begin(), pointSetIterator2 = m_PointSet2->Begin(); pointSetIterator != end; ++pointSetIterator, ++pointSetIterator2) //iterate simultaneously over both sets { point1 = pointSetIterator.Value(); point2 = pointSetIterator2.Value(); double squaredDistance = point1.SquaredEuclideanDistanceTo(point2); mean+=sqrt(squaredDistance); rms+=squaredDistance; this->m_SquaredDifferencesVector.push_back(squaredDistance); differencesVector.push_back(sqrt(squaredDistance)); } m_DifferencesVector = differencesVector; mean = mean/numberOfPoints; rms = sqrt(rms/numberOfPoints); for (std::vector::size_type i=0; i