diff --git a/Modules/ImageStatistics/Testing/mitkImageStatisticsCalculatorTest.cpp b/Modules/ImageStatistics/Testing/mitkImageStatisticsCalculatorTest.cpp index 10bd59e1c5..55bbdc70bb 100644 --- a/Modules/ImageStatistics/Testing/mitkImageStatisticsCalculatorTest.cpp +++ b/Modules/ImageStatistics/Testing/mitkImageStatisticsCalculatorTest.cpp @@ -1,1789 +1,1836 @@ /*=================================================================== 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 #include #include #include #include #include #include #include #include #include #include #include /** * \brief Test class for mitkImageStatisticsCalculator * * This test covers: * - instantiation of an ImageStatisticsCalculator class * - correctness of statistics when using PlanarFigures for masking */ class mitkImageStatisticsCalculatorTestSuite : public mitk::TestFixture { CPPUNIT_TEST_SUITE(mitkImageStatisticsCalculatorTestSuite); MITK_TEST(TestUninitializedImage); MITK_TEST(TestCase1); MITK_TEST(TestCase2); MITK_TEST(TestCase3); MITK_TEST(TestCase4); MITK_TEST(TestCase5); MITK_TEST(TestCase6); MITK_TEST(TestCase7); MITK_TEST(TestCase8); MITK_TEST(TestCase9); MITK_TEST(TestCase10); MITK_TEST(TestCase11); MITK_TEST(TestCase12); MITK_TEST(TestImageMaskingEmpty); MITK_TEST(TestImageMaskingNonEmpty); MITK_TEST(TestRecomputeOnModifiedMask); MITK_TEST(TestPic3DStatistics); MITK_TEST(TestPic3DAxialPlanarFigureMaskStatistics); MITK_TEST(TestPic3DSagittalPlanarFigureMaskStatistics); MITK_TEST(TestPic3DCoronalPlanarFigureMaskStatistics); MITK_TEST(TestPic3DImageMaskStatistics_label1); MITK_TEST(TestPic3DImageMaskStatistics_label2); MITK_TEST(TestPic3DIgnorePixelValueMaskStatistics); MITK_TEST(TestPic3DSecondaryMaskStatistics); MITK_TEST(TestUS4DCylStatistics_time1); MITK_TEST(TestUS4DCylAxialPlanarFigureMaskStatistics_time1); MITK_TEST(TestUS4DCylSagittalPlanarFigureMaskStatistics_time1); MITK_TEST(TestUS4DCylCoronalPlanarFigureMaskStatistics_time1); MITK_TEST(TestUS4DCylImageMaskStatistics_time1_label_1); MITK_TEST(TestUS4DCylImageMaskStatistics_time2_label_1); MITK_TEST(TestUS4DCylImageMaskStatistics_time1_label_2); MITK_TEST(TestUS4DCylIgnorePixelValueMaskStatistics_time1); MITK_TEST(TestUS4DCylSecondaryMaskStatistics_time1); CPPUNIT_TEST_SUITE_END(); public: void setUp() override; void tearDown() override; void TestUninitializedImage(); void TestCase1(); void TestCase2(); void TestCase3(); void TestCase4(); void TestCase5(); void TestCase6(); void TestCase7(); void TestCase8(); void TestCase9(); void TestCase10(); void TestCase11(); void TestCase12(); void TestImageMaskingEmpty(); void TestImageMaskingNonEmpty(); void TestRecomputeOnModifiedMask(); void TestPic3DStatistics(); void TestPic3DAxialPlanarFigureMaskStatistics(); void TestPic3DSagittalPlanarFigureMaskStatistics(); void TestPic3DCoronalPlanarFigureMaskStatistics(); void TestPic3DImageMaskStatistics_label1(); void TestPic3DImageMaskStatistics_label2(); void TestPic3DIgnorePixelValueMaskStatistics(); void TestPic3DSecondaryMaskStatistics(); void TestUS4DCylStatistics_time1(); void TestUS4DCylAxialPlanarFigureMaskStatistics_time1(); void TestUS4DCylSagittalPlanarFigureMaskStatistics_time1(); void TestUS4DCylCoronalPlanarFigureMaskStatistics_time1(); void TestUS4DCylImageMaskStatistics_time1_label_1(); void TestUS4DCylImageMaskStatistics_time2_label_1(); void TestUS4DCylImageMaskStatistics_time1_label_2(); void TestUS4DCylIgnorePixelValueMaskStatistics_time1(); void TestUS4DCylSecondaryMaskStatistics_time1(); void TestDifferentNBinsForHistogramStatistics(); void TestDifferentBinSizeForHistogramStatistic(); void TestSwitchFromBinSizeToNBins(); void TestSwitchFromNBinsToBinSize(); private: mitk::Image::Pointer m_TestImage; mitk::Image::Pointer m_Pic3DImage; mitk::Image::Pointer m_Pic3DImageMask; mitk::Image::Pointer m_Pic3DImageMask2; mitk::PlanarFigure::Pointer m_Pic3DPlanarFigureAxial; mitk::PlanarFigure::Pointer m_Pic3DPlanarFigureSagittal; mitk::PlanarFigure::Pointer m_Pic3DPlanarFigureCoronal; mitk::Image::Pointer m_US4DImage; mitk::Image::Pointer m_US4DImageMask; mitk::Image::Pointer m_US4DImageMask2; mitk::PlanarFigure::Pointer m_US4DPlanarFigureAxial; mitk::PlanarFigure::Pointer m_US4DPlanarFigureSagittal; mitk::PlanarFigure::Pointer m_US4DPlanarFigureCoronal; mitk::PlaneGeometry::Pointer m_Geometry; // calculate statistics for the given image and planarpolygon - const mitk::ImageStatisticsContainer::StatisticsObject ComputeStatistics( mitk::Image::Pointer image, + const mitk::ImageStatisticsContainer::Pointer ComputeStatistics( mitk::Image::Pointer image, mitk::PlanarFigure::Pointer polygon ); // calculate statistics for the given image and mask - const mitk::ImageStatisticsContainer::StatisticsObject ComputeStatistics( mitk::Image::Pointer image, + const mitk::ImageStatisticsContainer::Pointer ComputeStatistics(mitk::Image::Pointer image, mitk::Image::Pointer image_mask ); // universal function to calculate statistics - const mitk::ImageStatisticsContainer::StatisticsObject ComputeStatisticsNew(mitk::Image::Pointer image, - int timeStep=0, + const mitk::ImageStatisticsContainer::Pointer ComputeStatisticsNew(mitk::Image::Pointer image, mitk::MaskGenerator::Pointer maskGen=nullptr, mitk::MaskGenerator::Pointer secondardMaskGen=nullptr, unsigned short label=1); void VerifyStatistics(mitk::ImageStatisticsContainer::StatisticsObject stats, double testMean, double testSD, double testMedian=0); void VerifyStatistics(mitk::ImageStatisticsContainer::StatisticsObject stats, unsigned long N, double mean, double MPP, double median, double skewness, double kurtosis, double uniformity, double UPP, double variance, double stdev, double min, double max, double RMS, double entropy, vnl_vector minIndex, vnl_vector maxIndex); }; void mitkImageStatisticsCalculatorTestSuite::tearDown() { m_TestImage = nullptr; m_Pic3DImage = nullptr; m_Pic3DImageMask = nullptr; m_Pic3DImageMask2 = nullptr; m_Pic3DPlanarFigureAxial = nullptr; m_Pic3DPlanarFigureSagittal = nullptr; m_Pic3DPlanarFigureCoronal = nullptr; m_US4DImage = nullptr; m_US4DImageMask = nullptr; m_US4DImageMask2 = nullptr; m_US4DPlanarFigureAxial = nullptr; m_US4DPlanarFigureSagittal = nullptr; m_US4DPlanarFigureCoronal = nullptr; m_Geometry = nullptr; } void mitkImageStatisticsCalculatorTestSuite::setUp() { std::string filename = this->GetTestDataFilePath("ImageStatisticsTestData/testimage.dcm"); std::string Pic3DFile = this->GetTestDataFilePath("Pic3D.nrrd"); std::string Pic3DImageMaskFile = this->GetTestDataFilePath("ImageStatisticsTestData/Pic3D-labels.nrrd"); std::string Pic3DImageMaskFile2 = this->GetTestDataFilePath("ImageStatisticsTestData/Pic3D-labels2.nrrd"); std::string Pic3DAxialPlanarFigureFile = this->GetTestDataFilePath("ImageStatisticsTestData/Pic3DAxialPlanarFigure.pf"); std::string Pic3DSagittalPlanarFigureFile = this->GetTestDataFilePath("ImageStatisticsTestData/Pic3DSagittalPlanarFigure.pf"); std::string Pic3DCoronalPlanarFigureFile = this->GetTestDataFilePath("ImageStatisticsTestData/Pic3DCoronalPlanarFigure.pf"); std::string US4DFile = this->GetTestDataFilePath("US4DCyl.nrrd"); std::string US4DImageMaskFile = this->GetTestDataFilePath("ImageStatisticsTestData/US4D-labels.nrrd"); std::string US4DImageMaskFile2 = this->GetTestDataFilePath("ImageStatisticsTestData/US4D-labels2.nrrd"); std::string US4DAxialPlanarFigureFile = this->GetTestDataFilePath("ImageStatisticsTestData/US4DAxialPlanarFigure.pf"); std::string US4DSagittalPlanarFigureFile = this->GetTestDataFilePath("ImageStatisticsTestData/US4DSagittalPlanarFigure.pf"); std::string US4DCoronalPlanarFigureFile = this->GetTestDataFilePath("ImageStatisticsTestData/US4DCoronalPlanarFigure.pf"); if (filename.empty() || Pic3DFile.empty() || Pic3DImageMaskFile.empty() || Pic3DAxialPlanarFigureFile.empty() || Pic3DSagittalPlanarFigureFile.empty() || Pic3DCoronalPlanarFigureFile.empty() || US4DFile.empty() || US4DImageMaskFile.empty() || US4DAxialPlanarFigureFile.empty() || US4DSagittalPlanarFigureFile.empty() || US4DCoronalPlanarFigureFile.empty()) { MITK_TEST_FAILED_MSG( << "Could not find test file" ) } MITK_TEST_OUTPUT(<< "Loading test image '" << filename << "'") m_TestImage = mitk::IOUtil::Load(filename); MITK_TEST_CONDITION_REQUIRED( m_TestImage.IsNotNull(), "Loaded an mitk::Image" ); m_Geometry = m_TestImage->GetSlicedGeometry()->GetPlaneGeometry(0); MITK_TEST_CONDITION_REQUIRED( m_Geometry.IsNotNull(), "Getting image geometry" ); m_Pic3DImage = mitk::IOUtil::Load(Pic3DFile); MITK_TEST_CONDITION_REQUIRED( m_Pic3DImage.IsNotNull(), "Loaded Pic3D" ); m_Pic3DImageMask = mitk::IOUtil::Load(Pic3DImageMaskFile); MITK_TEST_CONDITION_REQUIRED( m_Pic3DImageMask.IsNotNull(), "Loaded Pic3D image mask" ); m_Pic3DImageMask2 = mitk::IOUtil::Load(Pic3DImageMaskFile2); MITK_TEST_CONDITION_REQUIRED( m_Pic3DImageMask2.IsNotNull(), "Loaded Pic3D image secondary mask" ); m_Pic3DPlanarFigureAxial = mitk::IOUtil::Load(Pic3DAxialPlanarFigureFile); MITK_TEST_CONDITION_REQUIRED( m_Pic3DPlanarFigureAxial.IsNotNull(), "Loaded Pic3D axial planarFigure" ); m_Pic3DPlanarFigureSagittal = mitk::IOUtil::Load(Pic3DSagittalPlanarFigureFile); MITK_TEST_CONDITION_REQUIRED( m_Pic3DPlanarFigureSagittal.IsNotNull(), "Loaded Pic3D sagittal planarFigure" ); m_Pic3DPlanarFigureCoronal = mitk::IOUtil::Load(Pic3DCoronalPlanarFigureFile); MITK_TEST_CONDITION_REQUIRED( m_Pic3DPlanarFigureCoronal.IsNotNull(), "Loaded Pic3D coronal planarFigure" ); m_US4DImage = mitk::IOUtil::Load(US4DFile); MITK_TEST_CONDITION_REQUIRED( m_US4DImage.IsNotNull(), "Loaded US4D" ); m_US4DImageMask = mitk::IOUtil::Load(US4DImageMaskFile); MITK_TEST_CONDITION_REQUIRED( m_US4DImageMask.IsNotNull(), "Loaded US4D image mask" ); m_US4DImageMask2 = mitk::IOUtil::Load(US4DImageMaskFile2); MITK_TEST_CONDITION_REQUIRED( m_US4DImageMask2.IsNotNull(), "Loaded US4D image mask2" ); m_US4DPlanarFigureAxial = mitk::IOUtil::Load(US4DAxialPlanarFigureFile); MITK_TEST_CONDITION_REQUIRED( m_US4DPlanarFigureAxial.IsNotNull(), "Loaded US4D axial planarFigure" ); m_US4DPlanarFigureSagittal = mitk::IOUtil::Load(US4DSagittalPlanarFigureFile); MITK_TEST_CONDITION_REQUIRED( m_US4DPlanarFigureSagittal.IsNotNull(), "Loaded US4D sagittal planarFigure" ); m_US4DPlanarFigureCoronal = mitk::IOUtil::Load(US4DCoronalPlanarFigureFile); MITK_TEST_CONDITION_REQUIRED( m_US4DPlanarFigureCoronal.IsNotNull(), "Loaded US4D coronal planarFigure" ); } void mitkImageStatisticsCalculatorTestSuite::TestCase1() { /***************************** * one whole white pixel * -> mean of 255 expected ******************************/ MITK_INFO << std::endl << "Test case 1:-----------------------------------------------------------------------------------"; mitk::PlanarPolygon::Pointer figure1 = mitk::PlanarPolygon::New(); figure1->SetPlaneGeometry( m_Geometry ); 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); - this->VerifyStatistics(ComputeStatistics(m_TestImage, figure1.GetPointer()), 255.0, 0.0, 255.0); + auto statisticsContainer = ComputeStatistics(m_TestImage, figure1.GetPointer()); + auto statisticsObjectTimestep0 = statisticsContainer->GetStatisticsForTimeStep(0); + + this->VerifyStatistics(statisticsObjectTimestep0, 255.0, 0.0, 255.0); } void mitkImageStatisticsCalculatorTestSuite::TestCase2() { /***************************** * half pixel in x-direction (white) * -> mean of 255 expected ******************************/ MITK_INFO << std::endl << "Test case 2:-----------------------------------------------------------------------------------"; mitk::PlanarPolygon::Pointer figure1 = mitk::PlanarPolygon::New(); figure1->SetPlaneGeometry( m_Geometry ); 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); - this->VerifyStatistics(ComputeStatistics(m_TestImage, figure1.GetPointer()), 255.0, 0.0, 255.0); + auto statisticsContainer = ComputeStatistics(m_TestImage, figure1.GetPointer()); + auto statisticsObjectTimestep0 = statisticsContainer->GetStatisticsForTimeStep(0); + + this->VerifyStatistics(statisticsObjectTimestep0, 255.0, 0.0, 255.0); } void mitkImageStatisticsCalculatorTestSuite::TestCase3() { /***************************** * half pixel in diagonal-direction (white) * -> mean of 255 expected ******************************/ MITK_INFO << std::endl << "Test case 3:-----------------------------------------------------------------------------------"; mitk::PlanarPolygon::Pointer figure1 = mitk::PlanarPolygon::New(); figure1->SetPlaneGeometry( m_Geometry ); 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); - this->VerifyStatistics(ComputeStatistics(m_TestImage, figure1.GetPointer()), 255.0, 0.0, 255.0); + auto statisticsContainer = ComputeStatistics(m_TestImage, figure1.GetPointer()); + auto statisticsObjectTimestep0 = statisticsContainer->GetStatisticsForTimeStep(0); + + this->VerifyStatistics(statisticsObjectTimestep0, 255.0, 0.0, 255.0); } void mitkImageStatisticsCalculatorTestSuite::TestCase4() { /***************************** * one pixel (white) + 2 half pixels (white) + 1 half pixel (black) * -> mean of 191.25 expected ******************************/ MITK_INFO << std::endl << "Test case 4:-----------------------------------------------------------------------------------"; mitk::PlanarPolygon::Pointer figure1 = mitk::PlanarPolygon::New(); figure1->SetPlaneGeometry( m_Geometry ); 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); - this->VerifyStatistics(ComputeStatistics(m_TestImage, figure1.GetPointer()), 191.25, 110.41, 242.250); + auto statisticsContainer = ComputeStatistics(m_TestImage, figure1.GetPointer()); + auto statisticsObjectTimestep0 = statisticsContainer->GetStatisticsForTimeStep(0); + + this->VerifyStatistics(statisticsObjectTimestep0, 191.25, 110.41, 242.250); } void mitkImageStatisticsCalculatorTestSuite::TestCase5() { /***************************** * whole pixel (white) + half pixel (gray) in x-direction * -> mean of 191.5 expected ******************************/ MITK_INFO << std::endl << "Test case 5:-----------------------------------------------------------------------------------"; mitk::PlanarPolygon::Pointer figure1 = mitk::PlanarPolygon::New(); figure1->SetPlaneGeometry( m_Geometry ); 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); - this->VerifyStatistics(ComputeStatistics(m_TestImage, figure1.GetPointer()), 191.50, 63.50, 134.340); + auto statisticsContainer = ComputeStatistics(m_TestImage, figure1.GetPointer()); + auto statisticsObjectTimestep0 = statisticsContainer->GetStatisticsForTimeStep(0); + + this->VerifyStatistics(statisticsObjectTimestep0, 191.50, 63.50, 134.340); } void mitkImageStatisticsCalculatorTestSuite::TestCase6() { /***************************** * quarter pixel (black) + whole pixel (white) + half pixel (gray) in x-direction * -> mean of 191.5 expected ******************************/ MITK_INFO << std::endl << "Test case 6:-----------------------------------------------------------------------------------"; mitk::PlanarPolygon::Pointer figure1 = mitk::PlanarPolygon::New(); figure1->SetPlaneGeometry( m_Geometry ); 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); - this->VerifyStatistics(ComputeStatistics(m_TestImage, figure1.GetPointer()), 191.5, 63.50, 134.340); + auto statisticsContainer = ComputeStatistics(m_TestImage, figure1.GetPointer()); + auto statisticsObjectTimestep0 = statisticsContainer->GetStatisticsForTimeStep(0); + + this->VerifyStatistics(statisticsObjectTimestep0, 191.5, 63.50, 134.340); } void mitkImageStatisticsCalculatorTestSuite::TestCase7() { /***************************** * half pixel (black) + whole pixel (white) + half pixel (gray) in x-direction * -> mean of 127.66 expected ******************************/ MITK_INFO << std::endl << "Test case 7:-----------------------------------------------------------------------------------"; mitk::PlanarPolygon::Pointer figure1 = mitk::PlanarPolygon::New(); figure1->SetPlaneGeometry( m_Geometry ); 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); - this->VerifyStatistics(ComputeStatistics(m_TestImage, figure1.GetPointer()), 127.66, 104.1, 140.250); + auto statisticsContainer = ComputeStatistics(m_TestImage, figure1.GetPointer()); + auto statisticsObjectTimestep0 = statisticsContainer->GetStatisticsForTimeStep(0); + + this->VerifyStatistics(statisticsObjectTimestep0, 127.66, 104.1, 140.250); } void mitkImageStatisticsCalculatorTestSuite::TestCase8() { /***************************** * whole pixel (gray) * -> mean of 128 expected ******************************/ MITK_INFO << std::endl << "Test case 8:-----------------------------------------------------------------------------------"; mitk::PlanarPolygon::Pointer figure2 = mitk::PlanarPolygon::New(); figure2->SetPlaneGeometry( m_Geometry ); 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); - this->VerifyStatistics(ComputeStatistics(m_TestImage, figure2.GetPointer()), 128.0, 0.0, 128.0); + auto statisticsContainer = ComputeStatistics(m_TestImage, figure2.GetPointer()); + auto statisticsObjectTimestep0 = statisticsContainer->GetStatisticsForTimeStep(0); + + this->VerifyStatistics(statisticsObjectTimestep0, 128.0, 0.0, 128.0); } void mitkImageStatisticsCalculatorTestSuite::TestCase9() { /***************************** * whole pixel (gray) + half pixel (white) in y-direction * -> mean of 191.5 expected ******************************/ MITK_INFO << std::endl << "Test case 9:-----------------------------------------------------------------------------------"; mitk::PlanarPolygon::Pointer figure2 = mitk::PlanarPolygon::New(); figure2->SetPlaneGeometry( m_Geometry ); 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); - this->VerifyStatistics(ComputeStatistics(m_TestImage, figure2.GetPointer()), 191.5, 63.50, 134.340); + auto statisticsContainer = ComputeStatistics(m_TestImage, figure2.GetPointer()); + auto statisticsObjectTimestep0 = statisticsContainer->GetStatisticsForTimeStep(0); + + this->VerifyStatistics(statisticsObjectTimestep0, 191.5, 63.50, 134.340); } void mitkImageStatisticsCalculatorTestSuite::TestCase10() { /***************************** * 2 whole pixel (white) + 2 whole pixel (black) in y-direction * -> mean of 127.66 expected ******************************/ MITK_INFO << std::endl << "Test case 10:-----------------------------------------------------------------------------------"; mitk::PlanarPolygon::Pointer figure2 = mitk::PlanarPolygon::New(); figure2->SetPlaneGeometry( m_Geometry ); 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); - this->VerifyStatistics(ComputeStatistics(m_TestImage, figure2.GetPointer()), 127.66, 104.1, 140.250); + auto statisticsContainer = ComputeStatistics(m_TestImage, figure2.GetPointer()); + auto statisticsObjectTimestep0 = statisticsContainer->GetStatisticsForTimeStep(0); + + this->VerifyStatistics(statisticsObjectTimestep0, 127.66, 104.1, 140.250); } void mitkImageStatisticsCalculatorTestSuite::TestCase11() { /***************************** * 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_INFO << std::endl << "Test case 11:-----------------------------------------------------------------------------------"; mitk::PlanarPolygon::Pointer figure2 = mitk::PlanarPolygon::New(); figure2->SetPlaneGeometry( m_Geometry ); 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); - this->VerifyStatistics(ComputeStatistics(m_TestImage, figure2.GetPointer()), 204.0, 102.00, 242.250); + auto statisticsContainer = ComputeStatistics(m_TestImage, figure2.GetPointer()); + auto statisticsObjectTimestep0 = statisticsContainer->GetStatisticsForTimeStep(0); + + this->VerifyStatistics(statisticsObjectTimestep0, 204.0, 102.00, 242.250); } void mitkImageStatisticsCalculatorTestSuite::TestCase12() { /***************************** * half pixel (white) + whole pixel (white) + half pixel (black) * -> mean of 212.66 expected ******************************/ MITK_INFO << std::endl << "Test case 12:-----------------------------------------------------------------------------------"; mitk::PlanarPolygon::Pointer figure2 = mitk::PlanarPolygon::New(); figure2->SetPlaneGeometry( m_Geometry ); 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); - this->VerifyStatistics(ComputeStatistics(m_TestImage, figure2.GetPointer()), 212.66, 59.860, 248.640); + auto statisticsContainer = ComputeStatistics(m_TestImage, figure2.GetPointer()); + auto statisticsObjectTimestep0 = statisticsContainer->GetStatisticsForTimeStep(0); + + this->VerifyStatistics(statisticsObjectTimestep0, 212.66, 59.860, 248.640); } void mitkImageStatisticsCalculatorTestSuite::TestImageMaskingEmpty() { MITK_INFO << std::endl << "TestImageMaskingEmpty:-----------------------------------------------------------------------------------"; mitk::Image::Pointer mask_image = mitk::ImageGenerator::GenerateImageFromReference( m_TestImage, 0 ); - auto statisticsObject = ComputeStatistics( m_TestImage, mask_image ); - // test if empty statisticsContainer - MITK_TEST_CONDITION(statisticsObject.m_Histogram.IsNull(), "histogram is null " << statisticsObject.m_Histogram.IsNull() << " histogram is null: " << true); - MITK_TEST_CONDITION(!statisticsObject.HasStatistic(mitk::ImageStatisticsConstants::MEAN()), "!has mean statistic: " << !statisticsObject.HasStatistic(mitk::ImageStatisticsConstants::MEAN()) << " expected: " << true); - MITK_TEST_CONDITION(!statisticsObject.HasStatistic(mitk::ImageStatisticsConstants::MEDIAN()), "!has median statistic: " << !statisticsObject.HasStatistic(mitk::ImageStatisticsConstants::MEDIAN()) << " expected: " << true); - MITK_TEST_CONDITION(!statisticsObject.HasStatistic(mitk::ImageStatisticsConstants::MINIMUM()), "!has min statistic: " << !statisticsObject.HasStatistic(mitk::ImageStatisticsConstants::MINIMUM()) << " expected: " << true); - MITK_TEST_CONDITION(!statisticsObject.HasStatistic(mitk::ImageStatisticsConstants::MAXIMUM()), "!has max statistic: " << !statisticsObject.HasStatistic(mitk::ImageStatisticsConstants::MAXIMUM()) << " expected: " << true); - MITK_TEST_CONDITION(!statisticsObject.HasStatistic(mitk::ImageStatisticsConstants::STANDARDDEVIATION()), "!has stddev statistic: " << !statisticsObject.HasStatistic(mitk::ImageStatisticsConstants::STANDARDDEVIATION()) << " expected: " << true); - MITK_TEST_CONDITION(!statisticsObject.HasStatistic(mitk::ImageStatisticsConstants::KURTOSIS()), "!has kurtosis statistic: " << !statisticsObject.HasStatistic(mitk::ImageStatisticsConstants::KURTOSIS()) << " expected: " << true); - MITK_TEST_CONDITION(!statisticsObject.HasStatistic(mitk::ImageStatisticsConstants::NUMBEROFVOXELS()), "!has numberOfVoxels statistic: " << !statisticsObject.HasStatistic(mitk::ImageStatisticsConstants::NUMBEROFVOXELS()) << " expected: " << true); + auto statisticsContainer = ComputeStatistics( m_TestImage, mask_image ); + // test if no statisticsContainer for timestep 0 exists + MITK_TEST_CONDITION(!statisticsContainer->TimeStepExists(0), "No statistics for TimeStep 0 does exist."); + MITK_TEST_FOR_EXCEPTION(mitk::Exception, statisticsContainer->GetStatisticsForTimeStep(0)); } void mitkImageStatisticsCalculatorTestSuite::TestImageMaskingNonEmpty() { MITK_INFO << std::endl << "TestImageMaskingNonEmpty:-----------------------------------------------------------------------------------"; mitk::Image::Pointer mask_image = mitk::ImageGenerator::GenerateImageFromReference( m_TestImage, 0 ); // activate voxel in the mask image if (mask_image->GetDimension() == 3) { std::vector< itk::Index<3U> > activated_indices; itk::Index<3U> index = { { 10, 8, 0 } }; activated_indices.push_back(index); index[0] = 9; index[1] = 8; index[2] = 0; activated_indices.push_back(index); index[0] = 9; index[1] = 7; index[2] = 0; activated_indices.push_back(index); index[0] = 10; index[1] = 7; index[2] = 0; activated_indices.push_back(index); std::vector< itk::Index<3U> >::const_iterator indexIter = activated_indices.begin(); mitk::ImagePixelWriteAccessor< unsigned char, 3> writeAccess(mask_image); while (indexIter != activated_indices.end()) { writeAccess.SetPixelByIndex((*indexIter++), 1); } } if (mask_image->GetDimension() == 4) { std::vector< itk::Index<4U> > activated_indices; itk::Index<4U> index = { { 10, 8, 0, 0 } }; activated_indices.push_back(index); index[0] = 9; index[1] = 8; index[2] = 0; index[3] = 0; activated_indices.push_back(index); index[0] = 9; index[1] = 7; index[2] = 0; index[3] = 0; activated_indices.push_back(index); index[0] = 10; index[1] = 7; index[2] = 0; index[3] = 0; activated_indices.push_back(index); std::vector< itk::Index<4U> >::const_iterator indexIter = activated_indices.begin(); mitk::ImagePixelWriteAccessor< unsigned char, 4> writeAccess(mask_image); while (indexIter != activated_indices.end()) { writeAccess.SetPixelByIndex((*indexIter++), 1); } } - this->VerifyStatistics( ComputeStatistics( m_TestImage, mask_image ), 127.5, 127.5, 12.750); + auto statisticsContainer = ComputeStatistics(m_TestImage, mask_image); + auto statisticsObjectTimestep0 = statisticsContainer->GetStatisticsForTimeStep(0); + + this->VerifyStatistics(statisticsObjectTimestep0, 127.5, 127.5, 12.750); } void mitkImageStatisticsCalculatorTestSuite::TestRecomputeOnModifiedMask() { MITK_INFO << std::endl << "TestRecomputeOnModifiedMask:-----------------------------------------------------------------------------------"; mitk::Image::Pointer mask_image = mitk::ImageGenerator::GenerateImageFromReference( m_TestImage, 0 ); mitk::ImageStatisticsCalculator::Pointer statisticsCalculator = mitk::ImageStatisticsCalculator::New(); statisticsCalculator->SetInputImage( m_TestImage ); mitk::ImageMaskGenerator::Pointer imgMaskGen = mitk::ImageMaskGenerator::New(); imgMaskGen->SetImageMask(mask_image); statisticsCalculator->SetMask(imgMaskGen.GetPointer()); - auto statisticsObject = statisticsCalculator->GetStatistics(); + auto statisticsContainer = statisticsCalculator->GetStatistics(); - // test if empty statisticsContainer (mask is empty) - MITK_TEST_CONDITION(statisticsObject.m_Histogram.IsNull(), "histogram is null " << statisticsObject.m_Histogram.IsNull() << " histogram is null: " << true); - MITK_TEST_CONDITION(!statisticsObject.HasStatistic(mitk::ImageStatisticsConstants::MEAN()), "!has mean statistic: " << !statisticsObject.HasStatistic(mitk::ImageStatisticsConstants::MEAN()) << " expected: " << true); - MITK_TEST_CONDITION(!statisticsObject.HasStatistic(mitk::ImageStatisticsConstants::MEDIAN()), "!has median statistic: " << !statisticsObject.HasStatistic(mitk::ImageStatisticsConstants::MEDIAN()) << " expected: " << true); - MITK_TEST_CONDITION(!statisticsObject.HasStatistic(mitk::ImageStatisticsConstants::MINIMUM()), "!has min statistic: " << !statisticsObject.HasStatistic(mitk::ImageStatisticsConstants::MINIMUM()) << " expected: " << true); - MITK_TEST_CONDITION(!statisticsObject.HasStatistic(mitk::ImageStatisticsConstants::MAXIMUM()), "!has max statistic: " << !statisticsObject.HasStatistic(mitk::ImageStatisticsConstants::MAXIMUM()) << " expected: " << true); - MITK_TEST_CONDITION(!statisticsObject.HasStatistic(mitk::ImageStatisticsConstants::STANDARDDEVIATION()), "!has stddev statistic: " << !statisticsObject.HasStatistic(mitk::ImageStatisticsConstants::STANDARDDEVIATION()) << " expected: " << true); - MITK_TEST_CONDITION(!statisticsObject.HasStatistic(mitk::ImageStatisticsConstants::KURTOSIS()), "!has kurtosis statistic: " << !statisticsObject.HasStatistic(mitk::ImageStatisticsConstants::KURTOSIS()) << " expected: " << true); - MITK_TEST_CONDITION(!statisticsObject.HasStatistic(mitk::ImageStatisticsConstants::NUMBEROFVOXELS()), "!has numberOfVoxels statistic: " << !statisticsObject.HasStatistic(mitk::ImageStatisticsConstants::NUMBEROFVOXELS()) << " expected: " << true); + // test if no statisticsContainer for timestep 0 exists + MITK_TEST_CONDITION(!statisticsContainer->TimeStepExists(0), "No statistics for TimeStep 0 does exist."); + MITK_TEST_FOR_EXCEPTION(mitk::Exception, statisticsContainer->GetStatisticsForTimeStep(0)); // activate voxel in the mask image if (mask_image->GetDimension() == 3) { itk::Index<3U> test_index = { { 11, 8, 0 } }; mitk::ImagePixelWriteAccessor< unsigned char, 3> writeAccess(mask_image); writeAccess.SetPixelByIndex(test_index, 1); } if (mask_image->GetDimension() == 4) { itk::Index<4U> test_index = { { 11, 8, 0, 0 } }; mitk::ImagePixelWriteAccessor< unsigned char, 4> writeAccess(mask_image); writeAccess.SetPixelByIndex(test_index, 1); } //Delete if T25625 has been resolved imgMaskGen->Modified(); - mitk::ImageStatisticsContainer::StatisticsObject stat = statisticsCalculator->GetStatistics(); - this->VerifyStatistics( stat, 128.0, 0.0, 128.0); - auto numberOfVoxels = stat.GetValueConverted(mitk::ImageStatisticsConstants::NUMBEROFVOXELS()); + + statisticsContainer = statisticsCalculator->GetStatistics(); + + MITK_TEST_CONDITION(statisticsContainer->TimeStepExists(0), "Statistics for TimeStep 0 does exist."); + auto statisticsObjectTimestep0 = statisticsContainer->GetStatisticsForTimeStep(0); + + this->VerifyStatistics(statisticsObjectTimestep0, 128.0, 0.0, 128.0); + auto numberOfVoxels = statisticsObjectTimestep0.GetValueConverted( + mitk::ImageStatisticsConstants::NUMBEROFVOXELS()); MITK_TEST_CONDITION(numberOfVoxels == 1, "Calculated mask voxel count '" << numberOfVoxels << "' is equal to the desired value '" << 1 << "'" ); } void mitkImageStatisticsCalculatorTestSuite::TestPic3DStatistics() { MITK_INFO << std::endl << "Test plain Pic3D:-----------------------------------------------------------------------------------"; unsigned long expected_N = 3211264; double expected_mean = -365.80015345982144; double expected_MPP = 111.80226129535752; double expected_median = -105.16000366210938; double expected_skewness = -0.26976612134147004; double expected_kurtosis = 1.4655017209571437; double expected_uniformity = 0.06087994379480554; double expected_UPP = 0.011227934437026977; double expected_variance = 224036.80150510342; double expected_standarddev = 473.32525973700518; double expected_min = -1023; double expected_max = 1361; double expected_RMS = 598.20276978323352; double expected_entropy = 4.6727423654570357; vnl_vector expected_minIndex; expected_minIndex.set_size(3); expected_minIndex[0] = 0; expected_minIndex[1] = 0; expected_minIndex[2] = 0; vnl_vector expected_maxIndex; expected_maxIndex.set_size(3); expected_maxIndex[0] = 139; expected_maxIndex[1] = 182; expected_maxIndex[2] = 43; - const mitk::ImageStatisticsContainer::StatisticsObject result = ComputeStatisticsNew(m_Pic3DImage, 0); + auto statisticsContainer = ComputeStatisticsNew(m_Pic3DImage); + auto statisticsObject = statisticsContainer->GetStatisticsForTimeStep(0); - VerifyStatistics(result, + VerifyStatistics(statisticsObject, expected_N, expected_mean, expected_MPP, expected_median, expected_skewness, expected_kurtosis, expected_uniformity, expected_UPP, expected_variance, expected_standarddev, expected_min, expected_max, expected_RMS, expected_entropy, expected_minIndex, expected_maxIndex); } void mitkImageStatisticsCalculatorTestSuite::TestPic3DAxialPlanarFigureMaskStatistics() { MITK_INFO << std::endl << "Test Pic3D axial pf:-----------------------------------------------------------------------------------"; double expected_entropy = 5.6719817476387417; double expected_kurtosis = 5.8846935191205221; double expected_MPP = 230.43933685003768; double expected_max = 1206; double expected_mean = 182.30282131661443; double expected_median = 95.970001220703125; double expected_min = -156; unsigned long expected_N = 3190; double expected_RMS = 301.93844376702253; double expected_skewness = 1.6400489794326298; double expected_standarddev = 240.69172225993557; double expected_UPP = 0.024889790784288681; double expected_uniformity = 0.027579917650180332; double expected_variance = 57932.505164453964; vnl_vector expected_minIndex; expected_minIndex.set_size(3); expected_minIndex[0] = 156; expected_minIndex[1] = 133; expected_minIndex[2] = 24; vnl_vector expected_maxIndex; expected_maxIndex.set_size(3); expected_maxIndex[0] = 125; expected_maxIndex[1] = 167; expected_maxIndex[2] = 24; mitk::PlanarFigureMaskGenerator::Pointer pfMaskGen = mitk::PlanarFigureMaskGenerator::New(); pfMaskGen->SetInputImage(m_Pic3DImage); pfMaskGen->SetPlanarFigure(m_Pic3DPlanarFigureAxial); - const mitk::ImageStatisticsContainer::StatisticsObject result = ComputeStatisticsNew(m_Pic3DImage, 0, pfMaskGen.GetPointer()); + auto statisticsContainer = ComputeStatisticsNew(m_Pic3DImage, pfMaskGen.GetPointer()); + auto statisticsObjectTimestep0 = statisticsContainer->GetStatisticsForTimeStep(0); - VerifyStatistics(result, + VerifyStatistics(statisticsObjectTimestep0, expected_N, expected_mean, expected_MPP, expected_median, expected_skewness, expected_kurtosis, expected_uniformity, expected_UPP, expected_variance, expected_standarddev, expected_min, expected_max, expected_RMS, expected_entropy, expected_minIndex, expected_maxIndex); } void mitkImageStatisticsCalculatorTestSuite::TestPic3DSagittalPlanarFigureMaskStatistics() { MITK_INFO << std::endl << "Test Pic3D sagittal pf:-----------------------------------------------------------------------------------"; double expected_entropy = 5.6051911962074286; double expected_kurtosis = 6.5814062739142338; double expected_MPP = 249.03202846975088; double expected_max = 1240; double expected_mean = 233.93602693602693; double expected_median = 174.9849853515625; double expected_min = -83; unsigned long expected_N = 1188; double expected_RMS = 332.03230188484594; double expected_skewness = 1.7489809015501814; double expected_standarddev = 235.62551813489128; double expected_UPP = 0.026837539253364174; double expected_uniformity = 0.027346982734188126; double expected_variance = 55519.384796335973; vnl_vector expected_minIndex; expected_minIndex.set_size(3); expected_minIndex[0] = 128; expected_minIndex[1] = 119; expected_minIndex[2] = 22; vnl_vector expected_maxIndex; expected_maxIndex.set_size(3); expected_maxIndex[0] = 128; expected_maxIndex[1] = 167; expected_maxIndex[2] = 22; mitk::PlanarFigureMaskGenerator::Pointer pfMaskGen = mitk::PlanarFigureMaskGenerator::New(); pfMaskGen->SetInputImage(m_Pic3DImage); pfMaskGen->SetPlanarFigure(m_Pic3DPlanarFigureSagittal); - const mitk::ImageStatisticsContainer::StatisticsObject result = ComputeStatisticsNew(m_Pic3DImage, 0, pfMaskGen.GetPointer()); + auto statisticsContainer = ComputeStatisticsNew(m_Pic3DImage, pfMaskGen.GetPointer()); + auto statisticsObjectTimestep0 = statisticsContainer->GetStatisticsForTimeStep(0); - VerifyStatistics(result, + VerifyStatistics(statisticsObjectTimestep0, expected_N, expected_mean, expected_MPP, expected_median, expected_skewness, expected_kurtosis, expected_uniformity, expected_UPP, expected_variance, expected_standarddev, expected_min, expected_max, expected_RMS, expected_entropy, expected_minIndex, expected_maxIndex); } void mitkImageStatisticsCalculatorTestSuite::TestPic3DCoronalPlanarFigureMaskStatistics() { MITK_INFO << std::endl << "Test Pic3D coronal pf:-----------------------------------------------------------------------------------"; double expected_entropy = 6.0677398647867449; double expected_kurtosis = 1.6242929941303372; double expected_MPP = 76.649350649350652; double expected_max = 156; double expected_mean = -482.14807692307693; double expected_median = -660.07501220703125; double expected_min = -897; unsigned long expected_N = 520; double expected_RMS = 595.09446729069839; double expected_skewness = 0.51691492278851858; double expected_standarddev = 348.81321207686312; double expected_UPP = 0.0021560650887573964; double expected_uniformity = 0.020295857988165685; double expected_variance = 121670.6569193787; vnl_vector expected_minIndex; expected_minIndex.set_size(3); expected_minIndex[0] = 217; expected_minIndex[1] = 127; expected_minIndex[2] = 43; vnl_vector expected_maxIndex; expected_maxIndex.set_size(3); expected_maxIndex[0] = 209; expected_maxIndex[1] = 127; expected_maxIndex[2] = 39; mitk::PlanarFigureMaskGenerator::Pointer pfMaskGen = mitk::PlanarFigureMaskGenerator::New(); pfMaskGen->SetInputImage(m_Pic3DImage); pfMaskGen->SetPlanarFigure(m_Pic3DPlanarFigureCoronal); - const mitk::ImageStatisticsContainer::StatisticsObject result = ComputeStatisticsNew(m_Pic3DImage, 0, pfMaskGen.GetPointer()); + auto statisticsContainer = ComputeStatisticsNew(m_Pic3DImage, pfMaskGen.GetPointer()); + auto statisticsObjectTimestep0 = statisticsContainer->GetStatisticsForTimeStep(0); - VerifyStatistics(result, + VerifyStatistics(statisticsObjectTimestep0, expected_N, expected_mean, expected_MPP, expected_median, expected_skewness, expected_kurtosis, expected_uniformity, expected_UPP, expected_variance, expected_standarddev, expected_min, expected_max, expected_RMS, expected_entropy, expected_minIndex, expected_maxIndex); } void mitkImageStatisticsCalculatorTestSuite::TestPic3DImageMaskStatistics_label1() { MITK_INFO << std::endl << "Test Pic3D image mask label 1 pf:-----------------------------------------------------------------------------------"; double expected_entropy = 5.695858251095868; double expected_kurtosis = 4.2728827997815717; double expected_MPP = 413.52408256880733; double expected_max = 1206; double expected_mean = 413.52408256880733; double expected_median = 324; double expected_min = 6; unsigned long expected_N = 872; double expected_RMS = 472.02024695145235; double expected_skewness = 1.3396074364415382; double expected_standarddev = 227.59821323493802; double expected_UPP = 0.029758648261930806; double expected_uniformity = 0.029758648261930806; double expected_variance = 51800.946667736309; vnl_vector expected_minIndex; expected_minIndex.set_size(3); expected_minIndex[0] = 135; expected_minIndex[1] = 158; expected_minIndex[2] = 24; vnl_vector expected_maxIndex; expected_maxIndex.set_size(3); expected_maxIndex[0] = 125; expected_maxIndex[1] = 167; expected_maxIndex[2] = 24; mitk::ImageMaskGenerator::Pointer imgMaskGen = mitk::ImageMaskGenerator::New(); imgMaskGen->SetImageMask(m_Pic3DImageMask); imgMaskGen->SetInputImage(m_Pic3DImage); imgMaskGen->SetTimeStep(0); - const mitk::ImageStatisticsContainer::StatisticsObject result = ComputeStatisticsNew(m_Pic3DImage, 0, imgMaskGen.GetPointer(), nullptr, 1); + auto statisticsContainer = ComputeStatisticsNew(m_Pic3DImage, imgMaskGen.GetPointer(), nullptr, 1); + auto statisticsObjectTimestep0 = statisticsContainer->GetStatisticsForTimeStep(0); - VerifyStatistics(result, + VerifyStatistics(statisticsObjectTimestep0, expected_N, expected_mean, expected_MPP, expected_median, expected_skewness, expected_kurtosis, expected_uniformity, expected_UPP, expected_variance, expected_standarddev, expected_min, expected_max, expected_RMS, expected_entropy, expected_minIndex, expected_maxIndex); } void mitkImageStatisticsCalculatorTestSuite::TestPic3DImageMaskStatistics_label2() { MITK_INFO << std::endl << "Test Pic3D image mask label 2 pf:-----------------------------------------------------------------------------------"; double expected_entropy = 4.3685781901212764; double expected_kurtosis = 9.7999112757587934; double expected_MPP = -nan(""); double expected_max = -145; double expected_mean = -897.92833876221493; double expected_median = -969.16499900817871; double expected_min = -1008; unsigned long expected_N = 307; double expected_RMS = 913.01496468179471; double expected_skewness = 2.6658524648889736; double expected_standarddev = 165.29072623903585; double expected_UPP = 0; double expected_uniformity = 0.087544695434434425; double expected_variance = 27321.024180627897; vnl_vector expected_minIndex; expected_minIndex.set_size(3); expected_minIndex[0] = 170; expected_minIndex[1] = 60; expected_minIndex[2] = 24; vnl_vector expected_maxIndex; expected_maxIndex.set_size(3); expected_maxIndex[0] = 173; expected_maxIndex[1] = 57; expected_maxIndex[2] = 24; mitk::ImageMaskGenerator::Pointer imgMaskGen = mitk::ImageMaskGenerator::New(); imgMaskGen->SetImageMask(m_Pic3DImageMask); imgMaskGen->SetInputImage(m_Pic3DImage); imgMaskGen->SetTimeStep(0); - const mitk::ImageStatisticsContainer::StatisticsObject result = ComputeStatisticsNew(m_Pic3DImage, 0, imgMaskGen.GetPointer(), nullptr, 2); + auto statisticsContainer = ComputeStatisticsNew(m_Pic3DImage, imgMaskGen.GetPointer(), nullptr, 2); + auto statisticsObjectTimestep0 = statisticsContainer->GetStatisticsForTimeStep(0); - VerifyStatistics(result, + VerifyStatistics(statisticsObjectTimestep0, expected_N, expected_mean, expected_MPP, expected_median, expected_skewness, expected_kurtosis, expected_uniformity, expected_UPP, expected_variance, expected_standarddev, expected_min, expected_max, expected_RMS, expected_entropy, expected_minIndex, expected_maxIndex); } void mitkImageStatisticsCalculatorTestSuite::TestPic3DIgnorePixelValueMaskStatistics() { MITK_INFO << std::endl << "Test Pic3D ignore zero pixels:-----------------------------------------------------------------------------------"; double expected_entropy = 4.671045011438645; double expected_kurtosis = 1.4638176488404484; double expected_MPP = 111.80226129535752; double expected_max = 1361; double expected_mean = -366.48547402877585; double expected_median = -105.16000366210938; double expected_min = -1023; unsigned long expected_N = 3205259; double expected_RMS = 598.76286909522139; double expected_skewness = -0.26648854845130782; double expected_standarddev = 473.50329537717545; double expected_UPP = 0.011270044547276429; double expected_uniformity = 0.061029773286547614; double expected_variance = 224205.37073304466; vnl_vector expected_minIndex; expected_minIndex.set_size(3); expected_minIndex[0] = 0; expected_minIndex[1] = 0; expected_minIndex[2] = 0; vnl_vector expected_maxIndex; expected_maxIndex.set_size(3); expected_maxIndex[0] = 139; expected_maxIndex[1] = 182; expected_maxIndex[2] = 43; mitk::IgnorePixelMaskGenerator::Pointer ignPixelValMask = mitk::IgnorePixelMaskGenerator::New(); ignPixelValMask->SetInputImage(m_Pic3DImage); ignPixelValMask->SetIgnoredPixelValue(0); ignPixelValMask->SetTimeStep(0); - const mitk::ImageStatisticsContainer::StatisticsObject result = ComputeStatisticsNew(m_Pic3DImage, 0, ignPixelValMask.GetPointer()); + auto statisticsContainer = ComputeStatisticsNew(m_Pic3DImage, ignPixelValMask.GetPointer()); + auto statisticsObjectTimestep0 = statisticsContainer->GetStatisticsForTimeStep(0); - VerifyStatistics(result, + VerifyStatistics(statisticsObjectTimestep0, expected_N, expected_mean, expected_MPP, expected_median, expected_skewness, expected_kurtosis, expected_uniformity, expected_UPP, expected_variance, expected_standarddev, expected_min, expected_max, expected_RMS, expected_entropy, expected_minIndex, expected_maxIndex); } void mitkImageStatisticsCalculatorTestSuite::TestPic3DSecondaryMaskStatistics() { MITK_INFO << std::endl << "Test Pic3D ignore zero pixels AND Image mask 2:-----------------------------------------------------------------------------------"; double expected_entropy = 5.9741637167320176; double expected_kurtosis = 3.490663358061596; double expected_MPP = 332.43534482758622; double expected_max = 1206; double expected_mean = 320.63333333333333; double expected_median = 265.06500244140625; double expected_min = -57; unsigned long expected_N = 720; double expected_RMS = 433.57749531594055; double expected_skewness = 1.1047775627624981; double expected_standarddev = 291.86248474238687; double expected_UPP = 0.020628858024691339; double expected_uniformity = 0.021377314814814797; double expected_variance = 85183.710000000006; vnl_vector expected_minIndex; expected_minIndex.set_size(3); expected_minIndex[0] = 116; expected_minIndex[1] = 170; expected_minIndex[2] = 24; vnl_vector expected_maxIndex; expected_maxIndex.set_size(3); expected_maxIndex[0] = 125; expected_maxIndex[1] = 167; expected_maxIndex[2] = 24; mitk::IgnorePixelMaskGenerator::Pointer ignPixelValMask = mitk::IgnorePixelMaskGenerator::New(); ignPixelValMask->SetInputImage(m_Pic3DImage); ignPixelValMask->SetIgnoredPixelValue(0); ignPixelValMask->SetTimeStep(0); mitk::ImageMaskGenerator::Pointer imgMaskGen2 = mitk::ImageMaskGenerator::New(); imgMaskGen2->SetImageMask(m_Pic3DImageMask2); imgMaskGen2->SetInputImage(m_Pic3DImage); imgMaskGen2->SetTimeStep(0); - const mitk::ImageStatisticsContainer::StatisticsObject result = ComputeStatisticsNew(m_Pic3DImage, 0, imgMaskGen2.GetPointer(), ignPixelValMask.GetPointer()); + auto statisticsContainer = + ComputeStatisticsNew(m_Pic3DImage, imgMaskGen2.GetPointer(), ignPixelValMask.GetPointer()); + auto statisticsObjectTimestep0 = statisticsContainer->GetStatisticsForTimeStep(0); - VerifyStatistics(result, + VerifyStatistics(statisticsObjectTimestep0, expected_N, expected_mean, expected_MPP, expected_median, expected_skewness, expected_kurtosis, expected_uniformity, expected_UPP, expected_variance, expected_standarddev, expected_min, expected_max, expected_RMS, expected_entropy, expected_minIndex, expected_maxIndex); } void mitkImageStatisticsCalculatorTestSuite::TestUS4DCylStatistics_time1() { MITK_INFO << std::endl << "Test plain US4D timeStep1:-----------------------------------------------------------------------------------"; double expected_entropy = 4.8272774900452502; double expected_kurtosis = 6.1336513352934432; double expected_MPP = 53.395358640738536; double expected_max = 199; double expected_mean = 35.771298153622375; double expected_median = 20.894999504089355; double expected_min = 0; unsigned long expected_N = 3409920; double expected_RMS = 59.244523377028408; double expected_skewness = 1.8734292240015058; double expected_standarddev = 47.226346233600559; double expected_UPP = 0.12098731125004937; double expected_uniformity = 0.12098731125004937; double expected_variance = 2230.3277785759178; vnl_vector expected_minIndex; expected_minIndex.set_size(3); expected_minIndex[0] = 0; expected_minIndex[1] = 0; expected_minIndex[2] = 0; vnl_vector expected_maxIndex; expected_maxIndex.set_size(3); expected_maxIndex[0] = 268; expected_maxIndex[1] = 101; expected_maxIndex[2] = 0; - const mitk::ImageStatisticsContainer::StatisticsObject result = ComputeStatisticsNew(m_US4DImage, 1); + auto statisticsContainer = ComputeStatisticsNew(m_US4DImage); + auto statisticsObjectTimestep1 = statisticsContainer->GetStatisticsForTimeStep(1); - VerifyStatistics(result, + VerifyStatistics(statisticsObjectTimestep1, expected_N, expected_mean, expected_MPP, expected_median, expected_skewness, expected_kurtosis, expected_uniformity, expected_UPP, expected_variance, expected_standarddev, expected_min, expected_max, expected_RMS, expected_entropy, expected_minIndex, expected_maxIndex); } void mitkImageStatisticsCalculatorTestSuite::TestUS4DCylAxialPlanarFigureMaskStatistics_time1() { MITK_INFO << std::endl << "Test US4D axial pf timeStep1:-----------------------------------------------------------------------------------"; double expected_entropy = 6.218151288002292; double expected_kurtosis = 1.7322676370242023; double expected_MPP = 121.11663807890223; double expected_max = 199; double expected_mean = 121.11663807890223; double expected_median = 120.14999771118164; double expected_min = 9; unsigned long expected_N = 2332; double expected_RMS = 134.41895158590751; double expected_skewness = -0.1454808104597369; double expected_standarddev = 58.30278317472294; double expected_UPP = 0.021354765820606133; double expected_uniformity = 0.021354765820606133; double expected_variance = 3399.214525918756; vnl_vector expected_minIndex; expected_minIndex.set_size(3); expected_minIndex[0] = 129; expected_minIndex[1] = 131; expected_minIndex[2] = 19; vnl_vector expected_maxIndex; expected_maxIndex.set_size(3); expected_maxIndex[0] = 126; expected_maxIndex[1] = 137; expected_maxIndex[2] = 19; mitk::PlanarFigureMaskGenerator::Pointer pfMaskGen = mitk::PlanarFigureMaskGenerator::New(); pfMaskGen->SetInputImage(m_US4DImage); pfMaskGen->SetPlanarFigure(m_US4DPlanarFigureAxial); - const mitk::ImageStatisticsContainer::StatisticsObject result = ComputeStatisticsNew(m_US4DImage, 1, pfMaskGen.GetPointer()); + auto statisticsContainer = ComputeStatisticsNew(m_US4DImage, pfMaskGen.GetPointer()); + auto statisticsObjectTimestep1 = statisticsContainer->GetStatisticsForTimeStep(1); - VerifyStatistics(result, + VerifyStatistics(statisticsObjectTimestep1, expected_N, expected_mean, expected_MPP, expected_median, expected_skewness, expected_kurtosis, expected_uniformity, expected_UPP, expected_variance, expected_standarddev, expected_min, expected_max, expected_RMS, expected_entropy, expected_minIndex, expected_maxIndex); } void mitkImageStatisticsCalculatorTestSuite::TestUS4DCylSagittalPlanarFigureMaskStatistics_time1() { MITK_INFO << std::endl << "Test US4D sagittal pf timeStep1:-----------------------------------------------------------------------------------"; double expected_entropy = 5.2003987046387508; double expected_kurtosis = 2.7574491062430142; double expected_MPP = 26.212534059945504; double expected_max = 59; double expected_mean = 26.176870748299319; double expected_median = 26.254999160766602; double expected_min = 0; unsigned long expected_N = 735; double expected_RMS = 28.084905283121476; double expected_skewness = 0.18245181360752327; double expected_standarddev = 10.175133541567705; double expected_UPP = 0.032921467906890628; double expected_uniformity = 0.032921467906890628; double expected_variance = 103.53334258873615; vnl_vector expected_minIndex; expected_minIndex.set_size(3); expected_minIndex[0] = 147; expected_minIndex[1] = 94; expected_minIndex[2] = 21; vnl_vector expected_maxIndex; expected_maxIndex.set_size(3); expected_maxIndex[0] = 147; expected_maxIndex[1] = 77; expected_maxIndex[2] = 24; mitk::PlanarFigureMaskGenerator::Pointer pfMaskGen = mitk::PlanarFigureMaskGenerator::New(); pfMaskGen->SetInputImage(m_US4DImage); pfMaskGen->SetPlanarFigure(m_US4DPlanarFigureSagittal); - const mitk::ImageStatisticsContainer::StatisticsObject result = ComputeStatisticsNew(m_US4DImage, 1, pfMaskGen.GetPointer()); + auto statisticsContainer = ComputeStatisticsNew(m_US4DImage, pfMaskGen.GetPointer()); + auto statisticsObjectTimestep1 = statisticsContainer->GetStatisticsForTimeStep(1); - VerifyStatistics(result, + VerifyStatistics(statisticsObjectTimestep1, expected_N, expected_mean, expected_MPP, expected_median, expected_skewness, expected_kurtosis, expected_uniformity, expected_UPP, expected_variance, expected_standarddev, expected_min, expected_max, expected_RMS, expected_entropy, expected_minIndex, expected_maxIndex); } void mitkImageStatisticsCalculatorTestSuite::TestUS4DCylCoronalPlanarFigureMaskStatistics_time1() { MITK_INFO << std::endl << "Test US4D coronal pf timeStep1:-----------------------------------------------------------------------------------"; double expected_entropy = 5.8892941136639161; double expected_kurtosis = 4.6434920707409564; double expected_MPP = 55.486426346239433; double expected_max = 199; double expected_mean = 55.118479221927501; double expected_median = 36.815000534057617; double expected_min = 0; unsigned long expected_N = 2262; double expected_RMS = 71.98149752438627; double expected_skewness = 1.4988288344523237; double expected_standarddev = 46.29567187238105; double expected_UPP = 0.023286748110675673; double expected_uniformity = 0.023286748110675673; double expected_variance = 2143.2892341151742; vnl_vector expected_minIndex; expected_minIndex.set_size(3); expected_minIndex[0] = 214; expected_minIndex[1] = 169; expected_minIndex[2] = 10; vnl_vector expected_maxIndex; expected_maxIndex.set_size(3); expected_maxIndex[0] = 99; expected_maxIndex[1] = 169; expected_maxIndex[2] = 17; mitk::PlanarFigureMaskGenerator::Pointer pfMaskGen = mitk::PlanarFigureMaskGenerator::New(); pfMaskGen->SetInputImage(m_US4DImage); pfMaskGen->SetPlanarFigure(m_US4DPlanarFigureCoronal); - const mitk::ImageStatisticsContainer::StatisticsObject result = ComputeStatisticsNew(m_US4DImage, 1, pfMaskGen.GetPointer()); + auto statisticsContainer = ComputeStatisticsNew(m_US4DImage, pfMaskGen.GetPointer()); + auto statisticsObjectTimestep1 = statisticsContainer->GetStatisticsForTimeStep(1); - VerifyStatistics(result, + VerifyStatistics(statisticsObjectTimestep1, expected_N, expected_mean, expected_MPP, expected_median, expected_skewness, expected_kurtosis, expected_uniformity, expected_UPP, expected_variance, expected_standarddev, expected_min, expected_max, expected_RMS, expected_entropy, expected_minIndex, expected_maxIndex); } void mitkImageStatisticsCalculatorTestSuite::TestUS4DCylImageMaskStatistics_time1_label_1() { MITK_INFO << std::endl << "Test US4D image mask time 1 label 1:-----------------------------------------------------------------------------------"; double expected_entropy = 5.0082903903398677; double expected_kurtosis = 3.6266994778237809; double expected_MPP = 169.58938547486034; double expected_max = 199; double expected_mean = 169.58938547486034; double expected_median = 187.44000244140625; double expected_min = 63; unsigned long expected_N = 716; double expected_RMS = 173.09843164831432; double expected_skewness = -1.2248969838579555; double expected_standarddev = 34.677188083311712; double expected_UPP = 0.076601073624418703; double expected_uniformity = 0.076601073624418703; double expected_variance = 1202.5073733653758; vnl_vector expected_minIndex; expected_minIndex.set_size(3); expected_minIndex[0] = 82; expected_minIndex[1] = 158; expected_minIndex[2] = 19; vnl_vector expected_maxIndex; expected_maxIndex.set_size(3); expected_maxIndex[0] = 126; expected_maxIndex[1] = 140; expected_maxIndex[2] = 19; mitk::ImageMaskGenerator::Pointer imgMask1 = mitk::ImageMaskGenerator::New(); imgMask1->SetInputImage(m_US4DImage); imgMask1->SetImageMask(m_US4DImageMask); - const mitk::ImageStatisticsContainer::StatisticsObject result = ComputeStatisticsNew(m_US4DImage, 1, imgMask1.GetPointer(), nullptr, 1); + auto statisticsContainer = ComputeStatisticsNew(m_US4DImage, imgMask1.GetPointer(), nullptr, 1); + auto statisticsObjectTimestep1 = statisticsContainer->GetStatisticsForTimeStep(1); - VerifyStatistics(result, + VerifyStatistics(statisticsObjectTimestep1, expected_N, expected_mean, expected_MPP, expected_median, expected_skewness, expected_kurtosis, expected_uniformity, expected_UPP, expected_variance, expected_standarddev, expected_min, expected_max, expected_RMS, expected_entropy, expected_minIndex, expected_maxIndex); } void mitkImageStatisticsCalculatorTestSuite::TestUS4DCylImageMaskStatistics_time2_label_1() { MITK_INFO << std::endl << "Test US4D image mask time 2 label 1:-----------------------------------------------------------------------------------"; double expected_entropy = 5.1857604214916506; double expected_kurtosis = 3.0692303858330683; double expected_MPP = 167.97194163860831; double expected_max = 199; double expected_mean = 167.97194163860831; double expected_median = 184.39499664306641; double expected_min = 72; unsigned long expected_N = 891; double expected_RMS = 171.67986611998634; double expected_skewness = -1.1221651136259736; double expected_standarddev = 35.488071983870803; double expected_UPP = 0.063124070232188439; double expected_uniformity = 0.063124070232188439; double expected_variance = 1259.4032531323958; vnl_vector expected_minIndex; expected_minIndex.set_size(3); expected_minIndex[0] = 103; expected_minIndex[1] = 212; expected_minIndex[2] = 19; vnl_vector expected_maxIndex; expected_maxIndex.set_size(3); expected_maxIndex[0] = 102; expected_maxIndex[1] = 168; expected_maxIndex[2] = 19; mitk::ImageMaskGenerator::Pointer imgMask1 = mitk::ImageMaskGenerator::New(); imgMask1->SetInputImage(m_US4DImage); imgMask1->SetImageMask(m_US4DImageMask); - const mitk::ImageStatisticsContainer::StatisticsObject result = ComputeStatisticsNew(m_US4DImage, 2, imgMask1.GetPointer(), nullptr, 1); + auto statisticsContainer = ComputeStatisticsNew(m_US4DImage, imgMask1.GetPointer(), nullptr, 1); + auto statisticsObjectTimestep2 = statisticsContainer->GetStatisticsForTimeStep(2); - VerifyStatistics(result, + VerifyStatistics(statisticsObjectTimestep2, expected_N, expected_mean, expected_MPP, expected_median, expected_skewness, expected_kurtosis, expected_uniformity, expected_UPP, expected_variance, expected_standarddev, expected_min, expected_max, expected_RMS, expected_entropy, expected_minIndex, expected_maxIndex); } void mitkImageStatisticsCalculatorTestSuite::TestUS4DCylImageMaskStatistics_time1_label_2() { MITK_INFO << std::endl << "Test US4D image mask time 1 label 2:-----------------------------------------------------------------------------------"; double expected_entropy = 5.0822234230119001; double expected_kurtosis = 2.4346603343623747; double expected_MPP = 20.733626373626375; double expected_max = 46; double expected_mean = 20.624836029733274; double expected_median = 20.010000228881836; double expected_min = 0; unsigned long expected_N = 2287; double expected_RMS = 22.508347574573804; double expected_skewness = 0.13837218490626488; double expected_standarddev = 9.0134260569684965; double expected_UPP = 0.034783970308787; double expected_uniformity = 0.034783970308787; double expected_variance = 81.241849284438644; vnl_vector expected_minIndex; expected_minIndex.set_size(3); expected_minIndex[0] = 178; expected_minIndex[1] = 76; expected_minIndex[2] = 19; vnl_vector expected_maxIndex; expected_maxIndex.set_size(3); expected_maxIndex[0] = 198; expected_maxIndex[1] = 90; expected_maxIndex[2] = 19; mitk::ImageMaskGenerator::Pointer imgMask1 = mitk::ImageMaskGenerator::New(); imgMask1->SetInputImage(m_US4DImage); imgMask1->SetImageMask(m_US4DImageMask); - const mitk::ImageStatisticsContainer::StatisticsObject result = ComputeStatisticsNew(m_US4DImage, 1, imgMask1.GetPointer(), nullptr, 2); + auto statisticsContainer = ComputeStatisticsNew(m_US4DImage, imgMask1.GetPointer(), nullptr, 2); + auto statisticsObjectTimestep1 = statisticsContainer->GetStatisticsForTimeStep(1); - VerifyStatistics(result, + VerifyStatistics(statisticsObjectTimestep1, expected_N, expected_mean, expected_MPP, expected_median, expected_skewness, expected_kurtosis, expected_uniformity, expected_UPP, expected_variance, expected_standarddev, expected_min, expected_max, expected_RMS, expected_entropy, expected_minIndex, expected_maxIndex); } void mitkImageStatisticsCalculatorTestSuite::TestUS4DCylIgnorePixelValueMaskStatistics_time1() { MITK_INFO << std::endl << "Test US4D ignore zero pixels:-----------------------------------------------------------------------------------"; double expected_entropy = 5.8609813848087962; double expected_kurtosis = 4.7556214582883651; double expected_MPP = 53.395358640738536; double expected_max = 199; double expected_mean = 53.395358640738536; double expected_median = 35.649999618530273; double expected_min = 1; unsigned long expected_N = 2284417; double expected_RMS = 72.382339046507084; double expected_skewness = 1.588289859859108; double expected_standarddev = 48.868585834566694; double expected_UPP = 0.023927063695115193; double expected_uniformity = 0.023927063695115193; double expected_variance = 2388.1386814704128; vnl_vector expected_minIndex; expected_minIndex.set_size(3); expected_minIndex[0] = 187; expected_minIndex[1] = 19; expected_minIndex[2] = 0; vnl_vector expected_maxIndex; expected_maxIndex.set_size(3); expected_maxIndex[0] = 268; expected_maxIndex[1] = 101; expected_maxIndex[2] = 0; mitk::IgnorePixelMaskGenerator::Pointer ignPixelValMask = mitk::IgnorePixelMaskGenerator::New(); ignPixelValMask->SetInputImage(m_US4DImage); ignPixelValMask->SetIgnoredPixelValue(0); ignPixelValMask->SetTimeStep(1); - const mitk::ImageStatisticsContainer::StatisticsObject result = ComputeStatisticsNew(m_US4DImage, 1, ignPixelValMask.GetPointer()); + auto statisticsContainer = ComputeStatisticsNew(m_US4DImage, ignPixelValMask.GetPointer()); + auto statisticsObjectTimestep1 = statisticsContainer->GetStatisticsForTimeStep(1); - VerifyStatistics(result, + VerifyStatistics(statisticsObjectTimestep1, expected_N, expected_mean, expected_MPP, expected_median, expected_skewness, expected_kurtosis, expected_uniformity, expected_UPP, expected_variance, expected_standarddev, expected_min, expected_max, expected_RMS, expected_entropy, expected_minIndex, expected_maxIndex); } void mitkImageStatisticsCalculatorTestSuite::TestUS4DCylSecondaryMaskStatistics_time1() { MITK_INFO << std::endl << "Test US4d ignore zero pixels AND Image mask 2:-----------------------------------------------------------------------------------"; double expected_entropy = 4.9955858614274558; double expected_kurtosis = 17.471042803365179; double expected_MPP = 32.791403286978507; double expected_max = 199; double expected_mean = 32.791403286978507; double expected_median = 25.75; double expected_min = 1; unsigned long expected_N = 17402; double expected_RMS = 42.776697859745241; double expected_skewness = 3.3991813038552596; double expected_standarddev = 27.469433016621732; double expected_UPP = 0.043040554251756687; double expected_uniformity = 0.043040554251756687; double expected_variance = 754.56975025466807; vnl_vector expected_minIndex; expected_minIndex.set_size(3); expected_minIndex[0] = 177; expected_minIndex[1] = 27; expected_minIndex[2] = 36; vnl_vector expected_maxIndex; expected_maxIndex.set_size(3); expected_maxIndex[0] = 109; expected_maxIndex[1] = 116; expected_maxIndex[2] = 36; mitk::IgnorePixelMaskGenerator::Pointer ignPixelValMask = mitk::IgnorePixelMaskGenerator::New(); ignPixelValMask->SetInputImage(m_US4DImage); ignPixelValMask->SetIgnoredPixelValue(0); mitk::ImageMaskGenerator::Pointer imgMaskGen2 = mitk::ImageMaskGenerator::New(); imgMaskGen2->SetImageMask(m_US4DImageMask2); imgMaskGen2->SetInputImage(m_US4DImage); - const mitk::ImageStatisticsContainer::StatisticsObject result = ComputeStatisticsNew(m_US4DImage, 1, imgMaskGen2.GetPointer(), ignPixelValMask.GetPointer()); - //std::cout << result->GetAsString(); + auto statisticsContainer = + ComputeStatisticsNew(m_US4DImage, imgMaskGen2.GetPointer(), ignPixelValMask.GetPointer()); + auto statisticsObjectTimestep1 = statisticsContainer->GetStatisticsForTimeStep(1); - VerifyStatistics(result, + VerifyStatistics(statisticsObjectTimestep1, expected_N, expected_mean, expected_MPP, expected_median, expected_skewness, expected_kurtosis, expected_uniformity, expected_UPP, expected_variance, expected_standarddev, expected_min, expected_max, expected_RMS, expected_entropy, expected_minIndex, expected_maxIndex); } -const mitk::ImageStatisticsContainer::StatisticsObject -mitkImageStatisticsCalculatorTestSuite::ComputeStatistics( mitk::Image::Pointer image, mitk::PlanarFigure::Pointer polygon ) +const mitk::ImageStatisticsContainer::Pointer mitkImageStatisticsCalculatorTestSuite::ComputeStatistics( mitk::Image::Pointer image, mitk::PlanarFigure::Pointer polygon ) { mitk::ImageStatisticsCalculator::Pointer statisticsCalculator = mitk::ImageStatisticsCalculator::New(); statisticsCalculator->SetInputImage( image ); statisticsCalculator->SetNBinsForHistogramStatistics(10); mitk::PlanarFigureMaskGenerator::Pointer planFigMaskGen = mitk::PlanarFigureMaskGenerator::New(); planFigMaskGen->SetInputImage(image); planFigMaskGen->SetPlanarFigure(polygon); statisticsCalculator->SetMask(planFigMaskGen.GetPointer()); try { return statisticsCalculator->GetStatistics(); } catch( ... ) { + return nullptr; } - - return mitk::ImageStatisticsContainer::StatisticsObject{}; } -const mitk::ImageStatisticsContainer::StatisticsObject +const mitk::ImageStatisticsContainer::Pointer mitkImageStatisticsCalculatorTestSuite::ComputeStatistics(mitk::Image::Pointer image, mitk::Image::Pointer image_mask ) { mitk::ImageStatisticsCalculator::Pointer statisticsCalculator = mitk::ImageStatisticsCalculator::New(); statisticsCalculator->SetInputImage(image); statisticsCalculator->SetNBinsForHistogramStatistics(10); mitk::ImageMaskGenerator::Pointer imgMaskGen = mitk::ImageMaskGenerator::New(); imgMaskGen->SetImageMask(image_mask); statisticsCalculator->SetMask(imgMaskGen.GetPointer()); return statisticsCalculator->GetStatistics(); } -const mitk::ImageStatisticsContainer::StatisticsObject +const mitk::ImageStatisticsContainer::Pointer mitkImageStatisticsCalculatorTestSuite::ComputeStatisticsNew(mitk::Image::Pointer image, - int timeStep, mitk::MaskGenerator::Pointer maskGen, mitk::MaskGenerator::Pointer secondardMaskGen, unsigned short label) { mitk::ImageStatisticsCalculator::Pointer imgStatCalc = mitk::ImageStatisticsCalculator::New(); imgStatCalc->SetInputImage(image); if (maskGen.IsNotNull()) { imgStatCalc->SetMask(maskGen.GetPointer()); if (secondardMaskGen.IsNotNull()) { imgStatCalc->SetSecondaryMask(secondardMaskGen.GetPointer()); } } - return imgStatCalc->GetStatistics(timeStep, label); + return imgStatCalc->GetStatistics(label); } void mitkImageStatisticsCalculatorTestSuite::VerifyStatistics(mitk::ImageStatisticsContainer::StatisticsObject stats, double testMean, double testSD, double testMedian) { auto mean = stats.GetValueConverted(mitk::ImageStatisticsConstants::MEAN()); int tmpMean = mean * 100; double calculatedMean = tmpMean / 100.0; MITK_TEST_CONDITION( calculatedMean == testMean, "Calculated mean grayvalue '" << calculatedMean << "' is equal to the desired value '" << testMean << "'" ); auto standardDeviation = stats.GetValueConverted(mitk::ImageStatisticsConstants::STANDARDDEVIATION()); int tmpSD = standardDeviation * 100; double calculatedSD = tmpSD / 100.0; MITK_TEST_CONDITION( calculatedSD == testSD, "Calculated grayvalue sd '" << calculatedSD << "' is equal to the desired value '" << testSD <<"'" ); auto median = stats.GetValueConverted(mitk::ImageStatisticsConstants::MEDIAN()); int tmpMedian = median * 100; double calculatedMedian = tmpMedian / 100.0; MITK_TEST_CONDITION( testMedian == calculatedMedian, "Calculated median grayvalue '" << calculatedMedian << "' is equal to the desired value '" << testMedian << "'"); } void mitkImageStatisticsCalculatorTestSuite::VerifyStatistics(mitk::ImageStatisticsContainer::StatisticsObject stats, unsigned long N, double mean, double MPP, double median, double skewness, double kurtosis, double uniformity, double UPP, double variance, double stdev, double min, double max, double RMS, double entropy, vnl_vector minIndex, vnl_vector maxIndex) { auto numberOfVoxelsObject = stats.GetValueConverted(mitk::ImageStatisticsConstants::NUMBEROFVOXELS()); auto meanObject = stats.GetValueConverted(mitk::ImageStatisticsConstants::MEAN()); auto mppObject = stats.GetValueConverted(mitk::ImageStatisticsConstants::MPP()); auto medianObject = stats.GetValueConverted(mitk::ImageStatisticsConstants::MEDIAN()); auto skewnessObject = stats.GetValueConverted(mitk::ImageStatisticsConstants::SKEWNESS()); auto kurtosisObject = stats.GetValueConverted(mitk::ImageStatisticsConstants::KURTOSIS()); auto uniformityObject = stats.GetValueConverted(mitk::ImageStatisticsConstants::UNIFORMITY()); auto uppObject = stats.GetValueConverted(mitk::ImageStatisticsConstants::UPP()); auto varianceObject = stats.GetValueConverted(mitk::ImageStatisticsConstants::VARIANCE()); auto standardDeviationObject = stats.GetValueConverted(mitk::ImageStatisticsConstants::STANDARDDEVIATION()); auto minObject = stats.GetValueConverted(mitk::ImageStatisticsConstants::MINIMUM()); auto maxObject = stats.GetValueConverted(mitk::ImageStatisticsConstants::MAXIMUM()); auto rmsObject = stats.GetValueConverted(mitk::ImageStatisticsConstants::RMS()); auto entropyObject = stats.GetValueConverted(mitk::ImageStatisticsConstants::ENTROPY()); auto minIndexObject = stats.GetValueConverted(mitk::ImageStatisticsConstants::MINIMUMPOSITION()); auto maxIndexObject = stats.GetValueConverted(mitk::ImageStatisticsConstants::MAXIMUMPOSITION()); MITK_TEST_CONDITION(numberOfVoxelsObject - N == 0, "calculated N: " << numberOfVoxelsObject << " expected N: " << N); MITK_TEST_CONDITION(std::abs(meanObject - mean) < mitk::eps, "calculated mean: " << meanObject << " expected mean: " << mean); // in one test case MPP is None because the roi has no positive pixels if (!std::isnan(mppObject)) { MITK_TEST_CONDITION(std::abs(mppObject - MPP) < mitk::eps, "calculated MPP: " << mppObject << " expected MPP: " << MPP); } MITK_TEST_CONDITION(std::abs(medianObject - median) < mitk::eps, "calculated median: " << medianObject << " expected median: " << median); MITK_TEST_CONDITION(std::abs(skewnessObject - skewness) < mitk::eps, "calculated skewness: " << skewnessObject << " expected skewness: " << skewness); MITK_TEST_CONDITION(std::abs(kurtosisObject - kurtosis) < mitk::eps, "calculated kurtosis: " << kurtosisObject << " expected kurtosis: " << kurtosis); MITK_TEST_CONDITION(std::abs(uniformityObject - uniformity) < mitk::eps, "calculated uniformity: " << uniformityObject << " expected uniformity: " << uniformity); MITK_TEST_CONDITION(std::abs(uppObject - UPP) < mitk::eps, "calculated UPP: " << uppObject << " expected UPP: " << UPP); MITK_TEST_CONDITION(std::abs(varianceObject - variance) < mitk::eps, "calculated variance: " << varianceObject << " expected variance: " << variance); MITK_TEST_CONDITION(std::abs(standardDeviationObject - stdev) < mitk::eps, "calculated stdev: " << standardDeviationObject << " expected stdev: " << stdev); MITK_TEST_CONDITION(std::abs(minObject - min) < mitk::eps, "calculated min: " << minObject << " expected min: " << min); MITK_TEST_CONDITION(std::abs(maxObject - max) < mitk::eps, "calculated max: " << maxObject << " expected max: " << max); MITK_TEST_CONDITION(std::abs(rmsObject - RMS) < mitk::eps, "calculated RMS: " << rmsObject << " expected RMS: " << RMS); MITK_TEST_CONDITION(std::abs(entropyObject - entropy) < mitk::eps, "calculated entropy: " << entropyObject << " expected entropy: " << entropy); for (unsigned int i = 0; i < minIndex.size(); ++i) { MITK_TEST_CONDITION(std::abs(minIndexObject[i] - minIndex[i]) < mitk::eps, "minIndex [" << i << "] = " << minIndexObject[i] << " expected: " << minIndex[i]); } for (unsigned int i = 0; i < maxIndex.size(); ++i) { MITK_TEST_CONDITION(std::abs(maxIndexObject[i] - maxIndex[i]) < mitk::eps, "maxIndex [" << i << "] = " << maxIndexObject[i] << " expected: " << maxIndex[i]); } } void mitkImageStatisticsCalculatorTestSuite::TestUninitializedImage() { /***************************** * loading uninitialized image to datastorage ******************************/ MITK_INFO << std::endl << "Test uninitialized image: -----------------------------------------------------------------------------------"; 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->GetStatistics(); MITK_TEST_FOR_EXCEPTION_END(mitk::Exception) } MITK_TEST_SUITE_REGISTRATION(mitkImageStatisticsCalculator) diff --git a/Modules/ImageStatistics/Testing/mitkImageStatisticsHotspotTest.cpp b/Modules/ImageStatistics/Testing/mitkImageStatisticsHotspotTest.cpp index 82856d604f..925f2f6210 100644 --- a/Modules/ImageStatistics/Testing/mitkImageStatisticsHotspotTest.cpp +++ b/Modules/ImageStatistics/Testing/mitkImageStatisticsHotspotTest.cpp @@ -1,653 +1,653 @@ /*=================================================================== 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 "itkMultiGaussianImageSource.h" #include "mitkTestingMacros.h" #include "mitkImageCast.h" #include #include #include #include #include #include #include /** \section hotspotCalculationTestCases Testcases To see the different Hotspot-Testcases have a look at the \ref hotspottestdoc. Note from an intensive session of checking the test results: - itk::MultiGaussianImageSource needs a review - the test idea is ok, but the combination of XML files for parameters and MultiGaussianImageSource has serious flaws - the XML file should contain exactly the parameters that MultiGaussianImageSource requires - in contrast, now the XML file mentions index coordinates for gaussian centers while the MultiGaussianImageSource expects world coordinates - this requires a transformation (index * spacing assuming no rotation) that was actually broken until recently */ 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 /** \brief XML-Tag "image-rows": size of x-dimension */ int m_ImageRows; /** \brief XML-Tag "image-columns": size of y-dimension */ int m_ImageColumns; /** \brief XML-Tag "image-slices": size of z-dimension */ int m_ImageSlices; /** \brief XML-Tag "numberOfGaussians": number of used gauss-functions */ int m_NumberOfGaussian; /** \brief XML-Tags "spacingX", "spacingY", "spacingZ": spacing of image in every direction */ double m_Spacing[3]; /** \brief XML-Tag "entireHotSpotInImage" */ unsigned int m_EntireHotspotInImage; // XML-Tag /** \brief XML-Tag "centerIndexX: gaussian parameter \warning This parameter READS the centerIndexX parameter from file and is THEN MISUSED to calculate some position in world coordinates, so we require double. */ std::vector m_CenterX; /** \brief XML-Tag "centerIndexY: gaussian parameter \warning This parameter READS the centerIndexX parameter from file and is THEN MISUSED to calculate some position in world coordinates, so we require double. */ std::vector m_CenterY; /** \brief XML-Tag "centerIndexZ: gaussian parameter \warning This parameter READS the centerIndexX parameter from file and is THEN MISUSED to calculate some position in world coordinates, so we require double. */ std::vector m_CenterZ; /** \brief XML-Tag "deviationX: gaussian parameter */ std::vector m_SigmaX; /** \brief XML-Tag "deviationY: gaussian parameter */ std::vector m_SigmaY; /** \brief XML-Tag "deviationZ: gaussian parameter */ std::vector m_SigmaZ; /** \brief XML-Tag "altitude: gaussian parameter */ std::vector m_Altitude; // XML-Tag /** \brief XML-Tag "numberOfLabels": number of different labels which appear in the mask */ unsigned int m_NumberOfLabels; /** \brief XML-Tag "hotspotRadiusInMM": radius of hotspot */ double m_HotspotRadiusInMM; // XML-Tag /** \brief XML-Tag "maximumSizeX": maximum position of ROI in x-dimension */ vnl_vector m_MaxIndexX; /** \brief XML-Tag "minimumSizeX": minimum position of ROI in x-dimension */ vnl_vector m_MinIndexX; /** \brief XML-Tag "maximumSizeX": maximum position of ROI in y-dimension */ vnl_vector m_MaxIndexY; /** \brief XML-Tag "minimumSizeX": minimum position of ROI in y-dimension */ vnl_vector m_MinIndexY; /** \brief XML-Tag "maximumSizeX": maximum position of ROI in z-dimension */ vnl_vector m_MaxIndexZ; /** \brief XML-Tag "minimumSizeX": minimum position of ROI in z-dimension */ vnl_vector m_MinIndexZ; /** \brief XML-Tag "label": value of label */ vnl_vector m_Label; //XML-Tag /** \brief XML-Tag "minimum": minimum inside hotspot */ vnl_vector m_HotspotMin; /** \brief XML-Tag "maximum": maximum inside hotspot */ vnl_vector m_HotspotMax; /** \brief XML-Tag "mean": mean value of hotspot */ vnl_vector m_HotspotMean; /** \brief XML-Tag "maximumIndexX": x-coordinate of maximum-location inside hotspot */ vnl_vector m_HotspotMaxIndexX; /** \brief XML-Tag "maximumIndexX": y-coordinate of maximum-location inside hotspot */ vnl_vector m_HotspotMaxIndexY; /** \brief XML-Tag "maximumIndexX": z-coordinate of maximum-location inside hotspot */ vnl_vector m_HotspotMaxIndexZ; /** \brief XML-Tag "maximumIndexX": x-coordinate of maximum-location inside hotspot */ vnl_vector m_HotspotMinIndexX; /** \brief XML-Tag "maximumIndexX": y-coordinate of maximum-location inside hotspot */ vnl_vector m_HotspotMinIndexY; /** \brief XML-Tag "maximumIndexX": z-coordinate of maximum-location inside hotspot */ vnl_vector m_HotspotMinIndexZ; /** \brief XML-Tag "maximumIndexX": x-coordinate of hotspot-location */ vnl_vector m_HotspotIndexX; /** \brief XML-Tag "maximumIndexX": y-coordinate of hotspot-location */ vnl_vector m_HotspotIndexY; /** \brief XML-Tag "maximumIndexX": z-coordinate of hotspot-location */ vnl_vector m_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 describing 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 four different tags, with all the details described as tag attributes. */ /** \verbatim \endverbatim */ static Parameters ParseParameters(int argc, char* argv[]) { 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; 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"); result.m_EntireHotspotInImage = GetIntegerAttribute( testimage, "entireHotSpotInImage" ); 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] = GetDoubleAttribute(gaussian, "deviationX"); result.m_SigmaY[i] = GetDoubleAttribute(gaussian, "deviationY"); result.m_SigmaZ[i] = GetDoubleAttribute(gaussian, "deviationZ"); result.m_Altitude[i] = GetDoubleAttribute(gaussian, "altitude"); result.m_CenterX[i] = result.m_CenterX[i] * result.m_Spacing[0]; result.m_CenterY[i] = result.m_CenterY[i] * result.m_Spacing[1]; result.m_CenterZ[i] = result.m_CenterZ[i] * result.m_Spacing[2]; result.m_SigmaX[i] = result.m_SigmaX[i] * result.m_Spacing[0]; result.m_SigmaY[i] = result.m_SigmaY[i] * result.m_Spacing[1]; result.m_SigmaZ[i] = result.m_SigmaZ[i] * result.m_Spacing[2]; } NodeList segmentations; domRoot->GetChildren("segmentation", segmentations); MITK_TEST_CONDITION_REQUIRED( segmentations.size() == 1, "One segmentation defined"); itk::DOMNode* segmentation = segmentations[0]; result.m_NumberOfLabels = GetIntegerAttribute(segmentation, "numberOfLabels"); result.m_HotspotRadiusInMM = GetDoubleAttribute(segmentation, "hotspotRadiusInMM"); // read ROI parameters, fill result structure NodeList rois; segmentation->GetChildren("roi", rois); MITK_TEST_CONDITION_REQUIRED( rois.size() >= 1, "At least one ROI defined" ) result.m_MaxIndexX.set_size(result.m_NumberOfLabels); result.m_MinIndexX.set_size(result.m_NumberOfLabels); result.m_MaxIndexY.set_size(result.m_NumberOfLabels); result.m_MinIndexY.set_size(result.m_NumberOfLabels); result.m_MaxIndexZ.set_size(result.m_NumberOfLabels); result.m_MinIndexZ.set_size(result.m_NumberOfLabels); result.m_Label.set_size(result.m_NumberOfLabels); for(unsigned int i = 0; i < rois.size(); ++i) { result.m_MaxIndexX[i] = GetIntegerAttribute(rois[i], "maximumIndexX"); result.m_MinIndexX[i] = GetIntegerAttribute(rois[i], "minimumIndexX"); result.m_MaxIndexY[i] = GetIntegerAttribute(rois[i], "maximumIndexY"); result.m_MinIndexY[i] = GetIntegerAttribute(rois[i], "minimumIndexY"); result.m_MaxIndexZ[i] = GetIntegerAttribute(rois[i], "maximumIndexZ"); result.m_MinIndexZ[i] = GetIntegerAttribute(rois[i], "minimumIndexZ"); result.m_Label[i] = GetIntegerAttribute(rois[i], "label"); } // read statistic parameters, fill result structure NodeList statistics; domRoot->GetChildren("statistic", statistics); MITK_TEST_CONDITION_REQUIRED( statistics.size() >= 1 , "At least one statistic defined" ) MITK_TEST_CONDITION_REQUIRED( statistics.size() == rois.size(), "Same number of rois and corresponding statistics defined"); result.m_HotspotMin.set_size(statistics.size()); result.m_HotspotMax.set_size(statistics.size()); result.m_HotspotMean.set_size(statistics.size()); result.m_HotspotMinIndexX.set_size(statistics.size()); result.m_HotspotMinIndexY.set_size(statistics.size()); result.m_HotspotMinIndexZ.set_size(statistics.size()); result.m_HotspotMaxIndexX.set_size(statistics.size()); result.m_HotspotMaxIndexY.set_size(statistics.size()); result.m_HotspotMaxIndexZ.set_size(statistics.size()); result.m_HotspotIndexX.set_size(statistics.size()); result.m_HotspotIndexY.set_size(statistics.size()); result.m_HotspotIndexZ.set_size(statistics.size()); for(unsigned int i = 0; i < statistics.size(); ++i) { result.m_HotspotMin[i] = GetDoubleAttribute(statistics[i], "minimum"); result.m_HotspotMax[i] = GetDoubleAttribute(statistics[i], "maximum"); result.m_HotspotMean[i] = GetDoubleAttribute(statistics[i], "mean"); result.m_HotspotMinIndexX[i] = GetIntegerAttribute(statistics[i], "minimumIndexX"); result.m_HotspotMinIndexY[i] = GetIntegerAttribute(statistics[i], "minimumIndexY"); result.m_HotspotMinIndexZ[i] = GetIntegerAttribute(statistics[i], "minimumIndexZ"); result.m_HotspotMaxIndexX[i] = GetIntegerAttribute(statistics[i], "maximumIndexX"); result.m_HotspotMaxIndexY[i] = GetIntegerAttribute(statistics[i], "maximumIndexY"); result.m_HotspotMaxIndexZ[i] = GetIntegerAttribute(statistics[i], "maximumIndexZ"); result.m_HotspotIndexX[i] = GetIntegerAttribute(statistics[i], "hotspotIndexX"); result.m_HotspotIndexY[i] = GetIntegerAttribute(statistics[i], "hotspotIndexY"); result.m_HotspotIndexZ[i] = GetIntegerAttribute(statistics[i], "hotspotIndexZ"); } } catch (std::exception& e) { MITK_TEST_CONDITION_REQUIRED(false, "Reading test parameters from XML file. Error message: " << e.what()); } return result; } /** \brief Generate an image that contains a couple of 3D gaussian distributions. Uses the given parameters to produce a test image using class MultiGaussianImageSource. */ static mitk::Image::Pointer BuildTestImage(const Parameters& testParameters) { mitk::Image::Pointer result; typedef double PixelType; const 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->SetRadius(testParameters.m_HotspotRadiusInMM); 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::ImageStatisticsContainer::StatisticsObject CalculateStatistics(mitk::Image* image, const Parameters& testParameters, unsigned int label) { 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(); statisticsCalculator->SetInputImage(image); mitk::Image::Pointer mitkMaskImage; if((testParameters.m_MaxIndexX[label] > testParameters.m_MinIndexX[label] && testParameters.m_MinIndexX[label] >= 0) && (testParameters.m_MaxIndexY[label] > testParameters.m_MinIndexY[label] && testParameters.m_MinIndexY[label] >= 0) && (testParameters.m_MaxIndexZ[label] > testParameters.m_MinIndexZ[label] && testParameters.m_MinIndexZ[label] >= 0)) { for(unsigned 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) { maskIt.Set(0); } for(unsigned int i = 0; i < testParameters.m_NumberOfLabels; ++i) { for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt) { MaskImageType::IndexType index = maskIt.GetIndex(); if((index[0] >= testParameters.m_MinIndexX[i] && index[0] <= testParameters.m_MaxIndexX[i] ) && (index[1] >= testParameters.m_MinIndexY[i] && index[1] <= testParameters.m_MaxIndexY[i] ) && (index[2] >= testParameters.m_MinIndexZ[i] && index[2] <= testParameters.m_MaxIndexZ[i] )) { maskIt.Set(testParameters.m_Label[i]); } } } mitk::CastToMitkImage(mask, mitkMaskImage); mitk::ImageMaskGenerator::Pointer imgMaskGen = mitk::ImageMaskGenerator::New(); imgMaskGen->SetImageMask(mitkMaskImage); mitk::HotspotMaskGenerator::Pointer hotspotMaskGen = mitk::HotspotMaskGenerator::New(); hotspotMaskGen->SetInputImage(image); hotspotMaskGen->SetLabel(testParameters.m_Label[label]); hotspotMaskGen->SetMask(imgMaskGen.GetPointer()); hotspotMaskGen->SetHotspotRadiusInMM(testParameters.m_HotspotRadiusInMM); if(testParameters.m_EntireHotspotInImage == 1) { MITK_INFO << "Hotspot must be completly inside image"; hotspotMaskGen->SetHotspotMustBeCompletelyInsideImage(true); } else { MITK_INFO << "Hotspot must not be completly inside image"; hotspotMaskGen->SetHotspotMustBeCompletelyInsideImage(false); } statisticsCalculator->SetMask(hotspotMaskGen.GetPointer()); MITK_DEBUG << "Masking is set to hotspot+image mask"; } else { mitk::HotspotMaskGenerator::Pointer hotspotMaskGen = mitk::HotspotMaskGenerator::New(); hotspotMaskGen->SetInputImage(image); hotspotMaskGen->SetHotspotRadiusInMM(testParameters.m_HotspotRadiusInMM); if(testParameters.m_EntireHotspotInImage == 1) { MITK_INFO << "Hotspot must be completly inside image"; hotspotMaskGen->SetHotspotMustBeCompletelyInsideImage(true); } else { MITK_INFO << "Hotspot must not be completly inside image"; hotspotMaskGen->SetHotspotMustBeCompletelyInsideImage(false); } MITK_DEBUG << "Masking is set to hotspot only"; } - return statisticsCalculator->GetStatistics(0); + return statisticsCalculator->GetStatistics()->GetStatisticsForTimeStep(0); } static void ValidateStatisticsItem(const std::string& label, double testvalue, double reference, double tolerance) { double diff = ::fabs(reference - testvalue); MITK_TEST_CONDITION( diff < tolerance, "'" << label << "' value close enough to reference value " "(value=" << testvalue << ", reference=" << reference << ", diff=" << diff << ")" ); } static void ValidateStatisticsItem(const std::string& label, const vnl_vector& testvalue, const vnl_vector& reference) { double diffX = ::fabs(double(testvalue[0] - reference[0])); double diffY = ::fabs(double(testvalue[1] - reference[1])); double diffZ = ::fabs(double(testvalue[2] - reference[2])); std::stringstream testPosition; testPosition << testvalue[0] << "," << testvalue[1] << "," << testvalue[2]; std::stringstream referencePosition; referencePosition << reference[0] << "," << reference[1] << "," << reference[2]; MITK_TEST_CONDITION( diffX < mitk::eps && diffY < mitk::eps && diffZ < mitk::eps, "'" << label << "' close enough to reference value " << "(value=[" << testPosition.str() << "]," << " reference=[" << referencePosition.str() << "]"); } /** \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::ImageStatisticsContainer::StatisticsObject hotspotStatistics, const Parameters& testParameters, unsigned int label) { // check all expected test result against actual results double eps = 0.25; // value above the largest tested difference auto mean = hotspotStatistics.GetValueConverted(mitk::ImageStatisticsConstants::MEAN()); auto max = hotspotStatistics.GetValueConverted(mitk::ImageStatisticsConstants::MAXIMUM()); auto min = hotspotStatistics.GetValueConverted(mitk::ImageStatisticsConstants::MINIMUM()); ValidateStatisticsItem("Hotspot mean", mean, testParameters.m_HotspotMean[label], eps); ValidateStatisticsItem("Hotspot maximum", max, testParameters.m_HotspotMax[label], eps); ValidateStatisticsItem("Hotspot minimum", min, testParameters.m_HotspotMin[label], eps); vnl_vector referenceHotspotCenterIndex; referenceHotspotCenterIndex.set_size(3); referenceHotspotCenterIndex[0] = testParameters.m_HotspotIndexX[label]; referenceHotspotCenterIndex[1] = testParameters.m_HotspotIndexY[label]; referenceHotspotCenterIndex[2] = testParameters.m_HotspotIndexZ[label]; // ValidateStatisticsItem("Hotspot center position", statistics.GetHotspotStatistics().GetHotspotIndex(), referenceHotspotCenterIndex); TODO: new image statistics calculator does not give hotspot position // TODO we do not test minimum/maximum positions within the peak/hotspot region, because // these positions are not unique, i.e. there are multiple valid minima/maxima positions. // One solution would be to modify the test cases in order to achive clear positions. // The BETTER/CORRECT solution would be to change the singular position into a set of positions / a region } }; /** \brief Verifies that hotspot statistics part of ImageStatisticsCalculator. The test reads parameters from an XML-file to generate a test-image, calculates the hotspot statistics of the image and checks if the calculated statistics are the same as the specified values of the XML-file. */ int mitkImageStatisticsHotspotTest(int argc, char* argv[]) { MITK_TEST_BEGIN("mitkImageStatisticsHotspotTest") try { mitkImageStatisticsHotspotTestClass::Parameters parameters = mitkImageStatisticsHotspotTestClass::ParseParameters(argc,argv); mitk::Image::Pointer image = mitkImageStatisticsHotspotTestClass::BuildTestImage(parameters); MITK_TEST_CONDITION_REQUIRED( image.IsNotNull(), "Generate test image" ); for(unsigned int label = 0; label < parameters.m_NumberOfLabels; ++label) { mitk::ImageStatisticsContainer::StatisticsObject statistics = mitkImageStatisticsHotspotTestClass::CalculateStatistics(image, parameters, label); mitkImageStatisticsHotspotTestClass::ValidateStatistics(statistics, parameters, label); std::cout << std::endl; } } catch (std::exception& e) { std::cout << "Error: " << e.what() << std::endl; MITK_TEST_CONDITION_REQUIRED( false, "Exception occurred during test execution: " << e.what() ); } catch(...) { MITK_TEST_CONDITION_REQUIRED( false, "Exception occurred during test execution." ); } MITK_TEST_END() } diff --git a/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp b/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp index 5c910668e6..481d709d8c 100644 --- a/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp +++ b/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp @@ -1,542 +1,565 @@ /*=================================================================== 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 +#include #include #include -#include -#include -#include +#include +#include #include +#include +#include #include #include #include -#include -#include -#include namespace mitk { - - void ImageStatisticsCalculator::SetInputImage(mitk::Image::Pointer image) + void ImageStatisticsCalculator::SetInputImage(mitk::Image::Pointer image) + { + if (image != m_Image) { - if (image != m_Image) - { - m_Image = image; - m_StatisticsUpdateTimePerTimeStep.resize(m_Image->GetTimeSteps()); - std::fill(m_StatisticsUpdateTimePerTimeStep.begin(), m_StatisticsUpdateTimePerTimeStep.end(), 0); - this->Modified(); - } + m_Image = image; + this->Modified(); } + } - void ImageStatisticsCalculator::SetMask(mitk::MaskGenerator::Pointer mask) + void ImageStatisticsCalculator::SetMask(mitk::MaskGenerator::Pointer mask) + { + if (mask != m_MaskGenerator) { - - if (mask != m_MaskGenerator) - { - m_MaskGenerator = mask; - this->Modified(); - } - + m_MaskGenerator = mask; + this->Modified(); } + } - void ImageStatisticsCalculator::SetSecondaryMask(mitk::MaskGenerator::Pointer mask) + void ImageStatisticsCalculator::SetSecondaryMask(mitk::MaskGenerator::Pointer mask) + { + if (mask != m_SecondaryMaskGenerator) { + m_SecondaryMaskGenerator = mask; + this->Modified(); + } + } - if (mask != m_SecondaryMaskGenerator) - { - m_SecondaryMaskGenerator = mask; - this->Modified(); - } - + void ImageStatisticsCalculator::SetNBinsForHistogramStatistics(unsigned int nBins) + { + if (nBins != m_nBinsForHistogramStatistics) + { + m_nBinsForHistogramStatistics = nBins; + this->Modified(); + this->m_UseBinSizeOverNBins = false; + } + if (m_UseBinSizeOverNBins) + { + this->Modified(); + this->m_UseBinSizeOverNBins = false; } + } + unsigned int ImageStatisticsCalculator::GetNBinsForHistogramStatistics() const + { + return m_nBinsForHistogramStatistics; + } - void ImageStatisticsCalculator::SetNBinsForHistogramStatistics(unsigned int nBins) + void ImageStatisticsCalculator::SetBinSizeForHistogramStatistics(double binSize) + { + if (binSize != m_binSizeForHistogramStatistics) { - if (nBins != m_nBinsForHistogramStatistics) - { - m_nBinsForHistogramStatistics = nBins; - this->Modified(); - this->m_UseBinSizeOverNBins = false; - } - if (m_UseBinSizeOverNBins) - { - this->Modified(); - this->m_UseBinSizeOverNBins = false; - } + m_binSizeForHistogramStatistics = binSize; + this->Modified(); + this->m_UseBinSizeOverNBins = true; } - - unsigned int ImageStatisticsCalculator::GetNBinsForHistogramStatistics() const + if (!m_UseBinSizeOverNBins) { - return m_nBinsForHistogramStatistics; + this->Modified(); + this->m_UseBinSizeOverNBins = true; } + } + + double ImageStatisticsCalculator::GetBinSizeForHistogramStatistics() const { return m_binSizeForHistogramStatistics; } - void ImageStatisticsCalculator::SetBinSizeForHistogramStatistics(double binSize) + mitk::ImageStatisticsContainer::Pointer ImageStatisticsCalculator::GetStatistics(LabelIndex label) + { + if (m_Image.IsNull()) { - if (binSize != m_binSizeForHistogramStatistics) - { - m_binSizeForHistogramStatistics = binSize; - this->Modified(); - this->m_UseBinSizeOverNBins = true; - } - if (!m_UseBinSizeOverNBins) - { - this->Modified(); - this->m_UseBinSizeOverNBins = true; - } + mitkThrow() << "no image"; } - double ImageStatisticsCalculator::GetBinSizeForHistogramStatistics() const + if (!m_Image->IsInitialized()) { - return m_binSizeForHistogramStatistics; + mitkThrow() << "Image not initialized!"; } - mitk::ImageStatisticsContainer::StatisticsObject ImageStatisticsCalculator::GetStatistics(unsigned int timeStep, unsigned int label) + if (IsUpdateRequired(label)) { - if (!m_StatisticContainers.empty() && timeStep >= m_StatisticContainers.back()->GetTimeSteps()) + // auto aStatisticContainer = ImageStatisticsContainer::New(); + auto timeGeometry = m_Image->GetTimeGeometry(); + // aStatisticContainer->SetTimeGeometry(timeGeometry); + // always compute statistics on all timesteps + for (unsigned int timeStep = 0; timeStep < m_Image->GetTimeSteps(); timeStep++) + { + if (m_MaskGenerator.IsNotNull()) { - mitkThrow() << "invalid timeStep in ImageStatisticsCalculator_v2::GetStatistics"; + m_MaskGenerator->SetTimeStep(timeStep); + //See T25625: otherwise, the mask is not computed again after setting a different time step + m_MaskGenerator->Modified(); + m_InternalMask = m_MaskGenerator->GetMask(); + if (m_MaskGenerator->GetReferenceImage().IsNotNull()) + { + m_InternalImageForStatistics = m_MaskGenerator->GetReferenceImage(); + } + else + { + m_InternalImageForStatistics = m_Image; + } } - - if (m_Image.IsNull()) + else { - mitkThrow() << "no image"; + m_InternalImageForStatistics = m_Image; } - if (!m_Image->IsInitialized()) + if (m_SecondaryMaskGenerator.IsNotNull()) { - mitkThrow() << "Image not initialized!"; + m_SecondaryMaskGenerator->SetTimeStep(timeStep); + m_SecondaryMask = m_SecondaryMaskGenerator->GetMask(); } - auto timeGeometry = m_Image->GetTimeGeometry(); - if (IsUpdateRequired(timeStep)) + ImageTimeSelector::Pointer imgTimeSel = ImageTimeSelector::New(); + imgTimeSel->SetInput(m_InternalImageForStatistics); + imgTimeSel->SetTimeNr(timeStep); + imgTimeSel->UpdateLargestPossibleRegion(); + imgTimeSel->Update(); + m_ImageTimeSlice = imgTimeSel->GetOutput(); + + // Calculate statistics with/without mask + if (m_MaskGenerator.IsNull() && m_SecondaryMaskGenerator.IsNull()) { - if (m_MaskGenerator.IsNotNull()) - { - m_MaskGenerator->SetTimeStep(timeStep); - m_InternalMask = m_MaskGenerator->GetMask(); - if (m_MaskGenerator->GetReferenceImage().IsNotNull()) - { - m_InternalImageForStatistics = m_MaskGenerator->GetReferenceImage(); - } - else - { - m_InternalImageForStatistics = m_Image; - } - } - else - { - m_InternalImageForStatistics = m_Image; - } - - if (m_SecondaryMaskGenerator.IsNotNull()) - { - m_SecondaryMaskGenerator->SetTimeStep(timeStep); - m_SecondaryMask = m_SecondaryMaskGenerator->GetMask(); - } - - ImageTimeSelector::Pointer imgTimeSel = ImageTimeSelector::New(); - imgTimeSel->SetInput(m_InternalImageForStatistics); - imgTimeSel->SetTimeNr(timeStep); - imgTimeSel->UpdateLargestPossibleRegion(); - m_ImageTimeSlice = imgTimeSel->GetOutput(); - - - // Calculate statistics with/without mask - if (m_MaskGenerator.IsNull() && m_SecondaryMaskGenerator.IsNull()) - { - // 1) calculate statistics unmasked: - AccessByItk_2(m_ImageTimeSlice, InternalCalculateStatisticsUnmasked, timeGeometry, timeStep) - - } - else - { - // 2) calculate statistics masked - AccessByItk_2(m_ImageTimeSlice, InternalCalculateStatisticsMasked, timeGeometry, timeStep) - } - - //this->Modified(); + // 1) calculate statistics unmasked: + AccessByItk_2(m_ImageTimeSlice, InternalCalculateStatisticsUnmasked, timeGeometry, timeStep) } - - m_StatisticsUpdateTimePerTimeStep[timeStep] = m_StatisticContainers[m_StatisticContainers.size() - 1]->GetMTime(); - - for (auto it = m_StatisticContainers.begin(); it != m_StatisticContainers.end(); ++it) + else { - ImageStatisticsContainer::Pointer statCont = *it; - ImageStatisticsContainer::StatisticsObject statObject = statCont->GetStatisticsForTimeStep(timeStep); - if (statObject.m_Label == label) - { - return statObject; - } + // 2) calculate statistics masked + AccessByItk_2(m_ImageTimeSlice, InternalCalculateStatisticsMasked, timeGeometry, timeStep) } - // these lines will ony be executed if the requested label could not be found! - MITK_WARN << "Invalid label: " << label << " in time step: " << timeStep; - return ImageStatisticsContainer::StatisticsObject(); + // this->Modified(); + } } - template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator::InternalCalculateStatisticsUnmasked( - typename itk::Image< TPixel, VImageDimension >* image, TimeGeometry* timeGeometry, TimeStepType timeStep) + auto it = m_StatisticContainers.find(label); + if (it != m_StatisticContainers.end()) + { + return it->second; + } + else + { + mitkThrow() << "unknown label"; + return nullptr; + } + } + + template + void ImageStatisticsCalculator::InternalCalculateStatisticsUnmasked( + typename itk::Image *image, TimeGeometry *timeGeometry, TimeStepType timeStep) + { + typedef typename itk::Image ImageType; + typedef typename itk::ExtendedStatisticsImageFilter ImageStatisticsFilterType; + typedef typename itk::MinMaxImageFilterWithIndex MinMaxFilterType; + + // reset statistics container if exists + ImageStatisticsContainer::Pointer statisticContainerForImage; + LabelIndex labelNoMask = 1; + auto it = m_StatisticContainers.find(labelNoMask); + if (it != m_StatisticContainers.end()) { - typedef typename itk::Image< TPixel, VImageDimension > ImageType; - typedef typename itk::ExtendedStatisticsImageFilter ImageStatisticsFilterType; - typedef typename itk::MinMaxImageFilterWithIndex MinMaxFilterType; - + statisticContainerForImage = it->second; + // statisticContainerForImage->Reset(); + // statisticContainerForImage->SetTimeGeometry(timeGeometry); + } + else + { + statisticContainerForImage = ImageStatisticsContainer::New(); + statisticContainerForImage->SetTimeGeometry(timeGeometry); + m_StatisticContainers.emplace(labelNoMask, statisticContainerForImage); + } - // reset statistics container if queried time step exists - if(m_StatisticContainers.empty() || m_StatisticContainers.back()->TimeStepExists(timeStep)) { - m_StatisticContainers.clear(); - ImageStatisticsContainer::Pointer statisticsResult = ImageStatisticsContainer::New(); - statisticsResult->SetTimeGeometry(timeGeometry); - m_StatisticContainers.push_back(statisticsResult); - } + auto statObj = ImageStatisticsContainer::StatisticsObject(); - //TODO for(timeStep) - auto statObj = ImageStatisticsContainer::StatisticsObject(); + typename ImageStatisticsFilterType::Pointer statisticsFilter = ImageStatisticsFilterType::New(); + statisticsFilter->SetInput(image); + statisticsFilter->SetCoordinateTolerance(0.001); + statisticsFilter->SetDirectionTolerance(0.001); - typename ImageStatisticsFilterType::Pointer statisticsFilter = ImageStatisticsFilterType::New(); - statisticsFilter->SetInput(image); - statisticsFilter->SetCoordinateTolerance(0.001); - statisticsFilter->SetDirectionTolerance(0.001); + // TODO: this is single threaded. Implement our own image filter that does this multi threaded + // typename itk::MinimumMaximumImageCalculator::Pointer imgMinMaxFilter = + // itk::MinimumMaximumImageCalculator::New(); imgMinMaxFilter->SetImage(image); + // imgMinMaxFilter->Compute(); + vnl_vector minIndex, maxIndex; - // TODO: this is single threaded. Implement our own image filter that does this multi threaded -// typename itk::MinimumMaximumImageCalculator::Pointer imgMinMaxFilter = itk::MinimumMaximumImageCalculator::New(); -// imgMinMaxFilter->SetImage(image); -// imgMinMaxFilter->Compute(); - vnl_vector minIndex, maxIndex; + typename MinMaxFilterType::Pointer minMaxFilter = MinMaxFilterType::New(); + minMaxFilter->SetInput(image); + minMaxFilter->UpdateLargestPossibleRegion(); + typename ImageType::PixelType minval = minMaxFilter->GetMin(); + typename ImageType::PixelType maxval = minMaxFilter->GetMax(); - typename MinMaxFilterType::Pointer minMaxFilter = MinMaxFilterType::New(); - minMaxFilter->SetInput(image); - minMaxFilter->UpdateLargestPossibleRegion(); - typename ImageType::PixelType minval = minMaxFilter->GetMin(); - typename ImageType::PixelType maxval = minMaxFilter->GetMax(); + typename ImageType::IndexType tmpMinIndex = minMaxFilter->GetMinIndex(); + typename ImageType::IndexType tmpMaxIndex = minMaxFilter->GetMaxIndex(); - typename ImageType::IndexType tmpMinIndex = minMaxFilter->GetMinIndex(); - typename ImageType::IndexType tmpMaxIndex = minMaxFilter->GetMaxIndex(); + // typename ImageType::IndexType tmpMinIndex = imgMinMaxFilter->GetIndexOfMinimum(); + // typename ImageType::IndexType tmpMaxIndex = imgMinMaxFilter->GetIndexOfMaximum(); -// typename ImageType::IndexType tmpMinIndex = imgMinMaxFilter->GetIndexOfMinimum(); -// typename ImageType::IndexType tmpMaxIndex = imgMinMaxFilter->GetIndexOfMaximum(); + minIndex.set_size(tmpMaxIndex.GetIndexDimension()); + maxIndex.set_size(tmpMaxIndex.GetIndexDimension()); - minIndex.set_size(tmpMaxIndex.GetIndexDimension()); - maxIndex.set_size(tmpMaxIndex.GetIndexDimension()); + for (unsigned int i = 0; i < tmpMaxIndex.GetIndexDimension(); i++) + { + minIndex[i] = tmpMinIndex[i]; + maxIndex[i] = tmpMaxIndex[i]; + } - for (unsigned int i=0; i < tmpMaxIndex.GetIndexDimension(); i++) - { - minIndex[i] = tmpMinIndex[i]; - maxIndex[i] = tmpMaxIndex[i]; - } + statObj.AddStatistic(mitk::ImageStatisticsConstants::MINIMUMPOSITION(), minIndex); + statObj.AddStatistic(mitk::ImageStatisticsConstants::MAXIMUMPOSITION(), maxIndex); - statObj.AddStatistic(mitk::ImageStatisticsConstants::MINIMUMPOSITION(), minIndex); - statObj.AddStatistic(mitk::ImageStatisticsConstants::MAXIMUMPOSITION(), maxIndex); + // convert m_binSize in m_nBins if necessary + unsigned int nBinsForHistogram; + if (m_UseBinSizeOverNBins) + { + nBinsForHistogram = std::max(static_cast(std::ceil(maxval - minval)) / m_binSizeForHistogramStatistics, + 10.); // do not allow less than 10 bins + } + else + { + nBinsForHistogram = m_nBinsForHistogramStatistics; + } - //convert m_binSize in m_nBins if necessary - unsigned int nBinsForHistogram; - if (m_UseBinSizeOverNBins) - { - nBinsForHistogram = std::max( - static_cast(std::ceil(maxval - minval)) / m_binSizeForHistogramStatistics, 10.); // do not allow less than 10 bins - } - else - { - nBinsForHistogram = m_nBinsForHistogramStatistics; - } + statisticsFilter->SetHistogramParameters(nBinsForHistogram, minval, maxval); - statisticsFilter->SetHistogramParameters(nBinsForHistogram, minval, maxval); + try + { + statisticsFilter->Update(); + } + catch (const itk::ExceptionObject &e) + { + mitkThrow() << "Image statistics calculation failed due to following ITK Exception: \n " << e.what(); + } - try - { - statisticsFilter->Update(); - } - catch (const itk::ExceptionObject& e) - { - mitkThrow() << "Image statistics calculation failed due to following ITK Exception: \n " << e.what(); - } + auto voxelVolume = GetVoxelVolume(image); + + auto numberOfPixels = image->GetLargestPossibleRegion().GetNumberOfPixels(); + auto volume = static_cast(numberOfPixels) * voxelVolume; + auto variance = statisticsFilter->GetSigma() * statisticsFilter->GetSigma(); + auto rms = + std::sqrt(std::pow(statisticsFilter->GetMean(), 2.) + statisticsFilter->GetVariance()); // variance = sigma^2 + + statObj.AddStatistic(mitk::ImageStatisticsConstants::NUMBEROFVOXELS(), + static_cast(numberOfPixels)); + statObj.AddStatistic(mitk::ImageStatisticsConstants::VOLUME(), volume); + statObj.AddStatistic(mitk::ImageStatisticsConstants::MEAN(), statisticsFilter->GetMean()); + statObj.AddStatistic(mitk::ImageStatisticsConstants::MINIMUM(), + static_cast(statisticsFilter->GetMinimum())); + statObj.AddStatistic(mitk::ImageStatisticsConstants::MAXIMUM(), + static_cast(statisticsFilter->GetMaximum())); + statObj.AddStatistic(mitk::ImageStatisticsConstants::STANDARDDEVIATION(), statisticsFilter->GetSigma()); + statObj.AddStatistic(mitk::ImageStatisticsConstants::VARIANCE(), variance); + statObj.AddStatistic(mitk::ImageStatisticsConstants::SKEWNESS(), statisticsFilter->GetSkewness()); + statObj.AddStatistic(mitk::ImageStatisticsConstants::KURTOSIS(), statisticsFilter->GetKurtosis()); + statObj.AddStatistic(mitk::ImageStatisticsConstants::RMS(), rms); + statObj.AddStatistic(mitk::ImageStatisticsConstants::MPP(), statisticsFilter->GetMPP()); + statObj.AddStatistic(mitk::ImageStatisticsConstants::ENTROPY(), statisticsFilter->GetEntropy()); + statObj.AddStatistic(mitk::ImageStatisticsConstants::MEDIAN(), statisticsFilter->GetMedian()); + statObj.AddStatistic(mitk::ImageStatisticsConstants::UNIFORMITY(), statisticsFilter->GetUniformity()); + statObj.AddStatistic(mitk::ImageStatisticsConstants::UPP(), statisticsFilter->GetUPP()); + statObj.m_Histogram = statisticsFilter->GetHistogram().GetPointer(); + statisticContainerForImage->SetStatisticsForTimeStep(timeStep, statObj); + } + + template + double ImageStatisticsCalculator::GetVoxelVolume(typename itk::Image *image) const + { + auto spacing = image->GetSpacing(); + double voxelVolume = 1.; + for (unsigned int i = 0; i < image->GetImageDimension(); i++) + { + voxelVolume *= spacing[i]; + } + return voxelVolume; + } + + template + void ImageStatisticsCalculator::InternalCalculateStatisticsMasked(typename itk::Image *image, + TimeGeometry *timeGeometry, + unsigned int timeStep) + { + typedef itk::Image ImageType; + typedef itk::Image MaskType; + typedef typename MaskType::PixelType LabelPixelType; + typedef itk::ExtendedLabelStatisticsImageFilter ImageStatisticsFilterType; + typedef MaskUtilities MaskUtilType; + typedef typename itk::MinMaxLabelImageFilterWithIndex MinMaxLabelFilterType; + typedef typename ImageType::PixelType InputImgPixelType; + + // workaround: if m_SecondaryMaskGenerator ist not null but m_MaskGenerator is! (this is the case if we request a + // 'ignore zuero valued pixels' mask in the gui but do not define a primary mask) + bool swapMasks = false; + if (m_SecondaryMask.IsNotNull() && m_InternalMask.IsNull()) + { + m_InternalMask = m_SecondaryMask; + m_SecondaryMask = nullptr; + swapMasks = true; + } - auto voxelVolume = GetVoxelVolume(image); - - // no mask, therefore just one label = the whole image - statObj.m_Label = 1; - auto numberOfPixels = image->GetLargestPossibleRegion().GetNumberOfPixels(); - auto volume = static_cast(numberOfPixels) * voxelVolume; - auto variance = statisticsFilter->GetSigma() * statisticsFilter->GetSigma(); - auto rms = std::sqrt(std::pow(statisticsFilter->GetMean(), 2.) + statisticsFilter->GetVariance()); // variance = sigma^2 - - statObj.AddStatistic(mitk::ImageStatisticsConstants::NUMBEROFVOXELS(), static_cast(numberOfPixels)); - statObj.AddStatistic(mitk::ImageStatisticsConstants::VOLUME(), volume); - statObj.AddStatistic(mitk::ImageStatisticsConstants::MEAN(), statisticsFilter->GetMean()); - statObj.AddStatistic(mitk::ImageStatisticsConstants::MINIMUM(), static_cast(statisticsFilter->GetMinimum())); - statObj.AddStatistic(mitk::ImageStatisticsConstants::MAXIMUM(), static_cast(statisticsFilter->GetMaximum())); - statObj.AddStatistic(mitk::ImageStatisticsConstants::STANDARDDEVIATION(), statisticsFilter->GetSigma()); - statObj.AddStatistic(mitk::ImageStatisticsConstants::VARIANCE(), variance); - statObj.AddStatistic(mitk::ImageStatisticsConstants::SKEWNESS(), statisticsFilter->GetSkewness()); - statObj.AddStatistic(mitk::ImageStatisticsConstants::KURTOSIS(), statisticsFilter->GetKurtosis()); - statObj.AddStatistic(mitk::ImageStatisticsConstants::RMS(), rms); - statObj.AddStatistic(mitk::ImageStatisticsConstants::MPP(), statisticsFilter->GetMPP()); - statObj.AddStatistic(mitk::ImageStatisticsConstants::ENTROPY(), statisticsFilter->GetEntropy()); - statObj.AddStatistic(mitk::ImageStatisticsConstants::MEDIAN(), statisticsFilter->GetMedian()); - statObj.AddStatistic(mitk::ImageStatisticsConstants::UNIFORMITY(), statisticsFilter->GetUniformity()); - statObj.AddStatistic(mitk::ImageStatisticsConstants::UPP(), statisticsFilter->GetUPP()); - statObj.m_Histogram = statisticsFilter->GetHistogram().GetPointer(); - - m_StatisticContainers.back()->SetStatisticsForTimeStep(timeStep, statObj); - } - - template < typename TPixel, unsigned int VImageDimension > - double ImageStatisticsCalculator::GetVoxelVolume(typename itk::Image* image) const - { - auto spacing = image->GetSpacing(); - double voxelVolume = 1.; - for (unsigned int i = 0; i < image->GetImageDimension(); i++) { - voxelVolume *= spacing[i]; - } - return voxelVolume; + // maskImage has to have the same dimension as image + typename MaskType::Pointer maskImage = MaskType::New(); + try + { + // try to access the pixel values directly (no copying or casting). Only works if mask pixels are of pixelType + // unsigned short + maskImage = ImageToItkImage(m_InternalMask); } + catch (const itk::ExceptionObject &) - template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator::InternalCalculateStatisticsMasked( - typename itk::Image< TPixel, VImageDimension >* image, TimeGeometry* timeGeometry, - unsigned int timeStep) { - typedef itk::Image< TPixel, VImageDimension > ImageType; - typedef itk::Image< MaskPixelType, VImageDimension > MaskType; - typedef typename MaskType::PixelType LabelPixelType; - typedef itk::ExtendedLabelStatisticsImageFilter< ImageType, MaskType > ImageStatisticsFilterType; - typedef MaskUtilities< TPixel, VImageDimension > MaskUtilType; - typedef typename itk::MinMaxLabelImageFilterWithIndex MinMaxLabelFilterType; - typedef typename ImageType::PixelType InputImgPixelType; + // if the pixel type of the mask is not short, then we have to make a copy of m_InternalMask (and cast the values) + CastToItkImage(m_InternalMask, maskImage); + } - // workaround: if m_SecondaryMaskGenerator ist not null but m_MaskGenerator is! (this is the case if we request a 'ignore zuero valued pixels' - // mask in the gui but do not define a primary mask) - bool swapMasks = false; - if (m_SecondaryMask.IsNotNull() && m_InternalMask.IsNull()) - { - m_InternalMask = m_SecondaryMask; - m_SecondaryMask = nullptr; - swapMasks = true; - } + // if we have a secondary mask (say a ignoreZeroPixelMask) we need to combine the masks (corresponds to AND) + if (m_SecondaryMask.IsNotNull()) + { + // dirty workaround for a bug when pf mask + any other mask is used in conjunction. We need a proper fix for this + // (Fabian Isensee is responsible and probably working on it!) + if (m_InternalMask->GetDimension() == 2 && + (m_SecondaryMask->GetDimension() == 3 || m_SecondaryMask->GetDimension() == 4)) + { + mitk::Image::Pointer old_img = m_SecondaryMaskGenerator->GetReferenceImage(); + m_SecondaryMaskGenerator->SetInputImage(m_MaskGenerator->GetReferenceImage()); + m_SecondaryMask = m_SecondaryMaskGenerator->GetMask(); + m_SecondaryMaskGenerator->SetInputImage(old_img); + } + typename MaskType::Pointer secondaryMaskImage = MaskType::New(); + secondaryMaskImage = ImageToItkImage(m_SecondaryMask); + + // secondary mask should be a ignore zero value pixel mask derived from image. it has to be cropped to the mask + // region (which may be planar or simply smaller) + typename MaskUtilities::Pointer secondaryMaskMaskUtil = + MaskUtilities::New(); + secondaryMaskMaskUtil->SetImage(secondaryMaskImage.GetPointer()); + secondaryMaskMaskUtil->SetMask(maskImage.GetPointer()); + typename MaskType::Pointer adaptedSecondaryMaskImage = secondaryMaskMaskUtil->ExtractMaskImageRegion(); + + typename itk::MaskImageFilter2::Pointer maskFilter = + itk::MaskImageFilter2::New(); + maskFilter->SetInput1(maskImage); + maskFilter->SetInput2(adaptedSecondaryMaskImage); + maskFilter->SetMaskingValue( + 1); // all pixels of maskImage where secondaryMaskImage==1 will be kept, all the others are set to 0 + maskFilter->UpdateLargestPossibleRegion(); + maskImage = maskFilter->GetOutput(); + } - // maskImage has to have the same dimension as image - typename MaskType::Pointer maskImage = MaskType::New(); - try { - // try to access the pixel values directly (no copying or casting). Only works if mask pixels are of pixelType unsigned short - maskImage = ImageToItkImage< MaskPixelType, VImageDimension >(m_InternalMask); - } - catch (const itk::ExceptionObject &) + typename MaskUtilType::Pointer maskUtil = MaskUtilType::New(); + maskUtil->SetImage(image); + maskUtil->SetMask(maskImage.GetPointer()); - { - // if the pixel type of the mask is not short, then we have to make a copy of m_InternalMask (and cast the values) - CastToItkImage(m_InternalMask, maskImage); - } + // if mask is smaller than image, extract the image region where the mask is + typename ImageType::Pointer adaptedImage = ImageType::New(); - // if we have a secondary mask (say a ignoreZeroPixelMask) we need to combine the masks (corresponds to AND) - if (m_SecondaryMask.IsNotNull()) - { - // dirty workaround for a bug when pf mask + any other mask is used in conjunction. We need a proper fix for this (Fabian Isensee is responsible and probably working on it!) - if (m_InternalMask->GetDimension() == 2 && (m_SecondaryMask->GetDimension() == 3 || m_SecondaryMask->GetDimension() == 4)) - { - mitk::Image::Pointer old_img = m_SecondaryMaskGenerator->GetReferenceImage(); - m_SecondaryMaskGenerator->SetInputImage(m_MaskGenerator->GetReferenceImage()); - m_SecondaryMask = m_SecondaryMaskGenerator->GetMask(); - m_SecondaryMaskGenerator->SetInputImage(old_img); - } - typename MaskType::Pointer secondaryMaskImage = MaskType::New(); - secondaryMaskImage = ImageToItkImage< MaskPixelType, VImageDimension >(m_SecondaryMask); - - // secondary mask should be a ignore zero value pixel mask derived from image. it has to be cropped to the mask region (which may be planar or simply smaller) - typename MaskUtilities::Pointer secondaryMaskMaskUtil = MaskUtilities::New(); - secondaryMaskMaskUtil->SetImage(secondaryMaskImage.GetPointer()); - secondaryMaskMaskUtil->SetMask(maskImage.GetPointer()); - typename MaskType::Pointer adaptedSecondaryMaskImage = secondaryMaskMaskUtil->ExtractMaskImageRegion(); - - typename itk::MaskImageFilter2::Pointer maskFilter = itk::MaskImageFilter2::New(); - maskFilter->SetInput1(maskImage); - maskFilter->SetInput2(adaptedSecondaryMaskImage); - maskFilter->SetMaskingValue(1); // all pixels of maskImage where secondaryMaskImage==1 will be kept, all the others are set to 0 - maskFilter->UpdateLargestPossibleRegion(); - maskImage = maskFilter->GetOutput(); - } + adaptedImage = maskUtil->ExtractMaskImageRegion(); // this also checks mask sanity - typename MaskUtilType::Pointer maskUtil = MaskUtilType::New(); - maskUtil->SetImage(image); - maskUtil->SetMask(maskImage.GetPointer()); + // find min, max, minindex and maxindex + typename MinMaxLabelFilterType::Pointer minMaxFilter = MinMaxLabelFilterType::New(); + minMaxFilter->SetInput(adaptedImage); + minMaxFilter->SetLabelInput(maskImage); + minMaxFilter->UpdateLargestPossibleRegion(); - // if mask is smaller than image, extract the image region where the mask is - typename ImageType::Pointer adaptedImage = ImageType::New(); + // set histogram parameters for each label individually (min/max may be different for each label) + typedef typename std::map MapType; + typedef typename std::pair PairType; - adaptedImage = maskUtil->ExtractMaskImageRegion(); // this also checks mask sanity + std::vector relevantLabels = minMaxFilter->GetRelevantLabels(); + MapType minVals; + MapType maxVals; + std::map nBins; - // find min, max, minindex and maxindex - typename MinMaxLabelFilterType::Pointer minMaxFilter = MinMaxLabelFilterType::New(); - minMaxFilter->SetInput(adaptedImage); - minMaxFilter->SetLabelInput(maskImage); - minMaxFilter->UpdateLargestPossibleRegion(); + for (LabelPixelType label : relevantLabels) + { + minVals.insert(PairType(label, minMaxFilter->GetMin(label))); + maxVals.insert(PairType(label, minMaxFilter->GetMax(label))); + + unsigned int nBinsForHistogram; + if (m_UseBinSizeOverNBins) + { + nBinsForHistogram = + std::max(static_cast(std::ceil(minMaxFilter->GetMax(label) - minMaxFilter->GetMin(label))) / + m_binSizeForHistogramStatistics, + 10.); // do not allow less than 10 bins + } + else + { + nBinsForHistogram = m_nBinsForHistogramStatistics; + } - // set histogram parameters for each label individually (min/max may be different for each label) - typedef typename std::map MapType; - typedef typename std::pair PairType; + nBins.insert(typename std::pair(label, nBinsForHistogram)); + } - std::vector relevantLabels = minMaxFilter->GetRelevantLabels(); - MapType minVals; - MapType maxVals; - std::map nBins; + typename ImageStatisticsFilterType::Pointer imageStatisticsFilter = ImageStatisticsFilterType::New(); + imageStatisticsFilter->SetDirectionTolerance(0.001); + imageStatisticsFilter->SetCoordinateTolerance(0.001); + imageStatisticsFilter->SetInput(adaptedImage); + imageStatisticsFilter->SetLabelInput(maskImage); + imageStatisticsFilter->SetHistogramParametersForLabels(nBins, minVals, maxVals); + imageStatisticsFilter->Update(); - for (LabelPixelType label:relevantLabels) - { - minVals.insert(PairType(label, minMaxFilter->GetMin(label))); - maxVals.insert(PairType(label, minMaxFilter->GetMax(label))); - - unsigned int nBinsForHistogram; - if (m_UseBinSizeOverNBins) - { - nBinsForHistogram = std::max(static_cast( - std::ceil(minMaxFilter->GetMax(label) - minMaxFilter->GetMin(label))) / m_binSizeForHistogramStatistics, 10.); // do not allow less than 10 bins - } - else - { - nBinsForHistogram = m_nBinsForHistogramStatistics; - } - - nBins.insert(typename std::pair(label, nBinsForHistogram)); - } + std::list labels = imageStatisticsFilter->GetRelevantLabels(); + auto it = labels.begin(); - typename ImageStatisticsFilterType::Pointer imageStatisticsFilter = ImageStatisticsFilterType::New(); - imageStatisticsFilter->SetDirectionTolerance(0.001); - imageStatisticsFilter->SetCoordinateTolerance(0.001); - imageStatisticsFilter->SetInput(adaptedImage); - imageStatisticsFilter->SetLabelInput(maskImage); - imageStatisticsFilter->SetHistogramParametersForLabels(nBins, minVals, maxVals); - imageStatisticsFilter->Update(); - - std::list labels = imageStatisticsFilter->GetRelevantLabels(); - auto it = labels.begin(); - - if (m_StatisticContainers.empty() || m_StatisticContainers.back()->TimeStepExists(timeStep)) { - m_StatisticContainers.clear(); - for (size_t i = 0; i < labels.size(); i++) { - ImageStatisticsContainer::Pointer statisticsResult = ImageStatisticsContainer::New(); - statisticsResult->SetTimeGeometry(timeGeometry); - m_StatisticContainers.push_back(statisticsResult); - } - } + while (it != labels.end()) + { + std::cout << "label: " << *it << std::endl; + ImageStatisticsContainer::Pointer statisticContainerForLabelImage; + auto labelIt = m_StatisticContainers.find(*it); + // reset if statisticContainer already exist + if (labelIt != m_StatisticContainers.end()) + { + statisticContainerForLabelImage = labelIt->second; + // statisticContainerForLabelImage->Reset(); + // statisticContainerForLabelImage->SetTimeGeometry(timeGeometry); + } + // create new statisticContainer + else + { + statisticContainerForLabelImage = ImageStatisticsContainer::New(); + statisticContainerForLabelImage->SetTimeGeometry(timeGeometry); + // link label (*it) to statisticContainer + m_StatisticContainers.emplace(*it, statisticContainerForLabelImage); + } - //m_StatisticsByTimeStep[timeStep].resize(0); + ImageStatisticsContainer::StatisticsObject statObj; + + // find min, max, minindex and maxindex + // make sure to only look in the masked region, use a masker for this + + vnl_vector minIndex, maxIndex; + mitk::Point3D worldCoordinateMin; + mitk::Point3D worldCoordinateMax; + mitk::Point3D indexCoordinateMin; + mitk::Point3D indexCoordinateMax; + m_InternalImageForStatistics->GetGeometry()->IndexToWorld(minMaxFilter->GetMinIndex(*it), worldCoordinateMin); + m_InternalImageForStatistics->GetGeometry()->IndexToWorld(minMaxFilter->GetMaxIndex(*it), worldCoordinateMax); + m_Image->GetGeometry()->WorldToIndex(worldCoordinateMin, indexCoordinateMin); + m_Image->GetGeometry()->WorldToIndex(worldCoordinateMax, indexCoordinateMax); + + minIndex.set_size(3); + maxIndex.set_size(3); + + // for (unsigned int i=0; i < tmpMaxIndex.GetIndexDimension(); i++) + for (unsigned int i = 0; i < 3; i++) + { + minIndex[i] = indexCoordinateMin[i]; + maxIndex[i] = indexCoordinateMax[i]; + } - while(it != labels.end()) - { - ImageStatisticsContainer::StatisticsObject statObj = ImageStatisticsContainer::StatisticsObject(); - - // find min, max, minindex and maxindex - // make sure to only look in the masked region, use a masker for this - - vnl_vector minIndex, maxIndex; - mitk::Point3D worldCoordinateMin; - mitk::Point3D worldCoordinateMax; - mitk::Point3D indexCoordinateMin; - mitk::Point3D indexCoordinateMax; - m_InternalImageForStatistics->GetGeometry()->IndexToWorld(minMaxFilter->GetMinIndex(*it), worldCoordinateMin); - m_InternalImageForStatistics->GetGeometry()->IndexToWorld(minMaxFilter->GetMaxIndex(*it), worldCoordinateMax); - m_Image->GetGeometry()->WorldToIndex(worldCoordinateMin, indexCoordinateMin); - m_Image->GetGeometry()->WorldToIndex(worldCoordinateMax, indexCoordinateMax); - - minIndex.set_size(3); - maxIndex.set_size(3); - - //for (unsigned int i=0; i < tmpMaxIndex.GetIndexDimension(); i++) - for (unsigned int i=0; i < 3; i++) - { - minIndex[i] = indexCoordinateMin[i]; - maxIndex[i] = indexCoordinateMax[i]; - } - - statObj.AddStatistic(mitk::ImageStatisticsConstants::MINIMUMPOSITION(), minIndex); - statObj.AddStatistic(mitk::ImageStatisticsConstants::MAXIMUMPOSITION(), maxIndex); - - assert(std::abs(minMaxFilter->GetMax(*it) - imageStatisticsFilter->GetMaximum(*it)) < mitk::eps); - assert(std::abs(minMaxFilter->GetMin(*it) - imageStatisticsFilter->GetMinimum(*it)) < mitk::eps); - - auto voxelVolume = GetVoxelVolume(image); - auto numberOfVoxels = static_cast(imageStatisticsFilter->GetSum(*it) / (double)imageStatisticsFilter->GetMean(*it)); - auto volume = static_cast(numberOfVoxels) * voxelVolume; - auto rms = std::sqrt(std::pow(imageStatisticsFilter->GetMean(*it), 2.) + imageStatisticsFilter->GetVariance(*it)); // variance = sigma^2 - auto variance = imageStatisticsFilter->GetSigma(*it)*imageStatisticsFilter->GetSigma(*it); - - statObj.m_Label = *it; - statObj.AddStatistic(mitk::ImageStatisticsConstants::NUMBEROFVOXELS(), numberOfVoxels); - statObj.AddStatistic(mitk::ImageStatisticsConstants::VOLUME(), volume); - statObj.AddStatistic(mitk::ImageStatisticsConstants::MEAN(), imageStatisticsFilter->GetMean(*it)); - statObj.AddStatistic(mitk::ImageStatisticsConstants::MINIMUM(), imageStatisticsFilter->GetMinimum(*it)); - statObj.AddStatistic(mitk::ImageStatisticsConstants::MAXIMUM(), imageStatisticsFilter->GetMaximum(*it)); - statObj.AddStatistic(mitk::ImageStatisticsConstants::STANDARDDEVIATION(), imageStatisticsFilter->GetSigma(*it)); - statObj.AddStatistic(mitk::ImageStatisticsConstants::VARIANCE(), variance); - statObj.AddStatistic(mitk::ImageStatisticsConstants::SKEWNESS(), imageStatisticsFilter->GetSkewness(*it)); - statObj.AddStatistic(mitk::ImageStatisticsConstants::KURTOSIS(), imageStatisticsFilter->GetKurtosis(*it)); - statObj.AddStatistic(mitk::ImageStatisticsConstants::RMS(), rms); - statObj.AddStatistic(mitk::ImageStatisticsConstants::MPP(), imageStatisticsFilter->GetMPP(*it)); - statObj.AddStatistic(mitk::ImageStatisticsConstants::ENTROPY(), imageStatisticsFilter->GetEntropy(*it)); - statObj.AddStatistic(mitk::ImageStatisticsConstants::MEDIAN(), imageStatisticsFilter->GetMedian(*it)); - statObj.AddStatistic(mitk::ImageStatisticsConstants::UNIFORMITY(), imageStatisticsFilter->GetUniformity(*it)); - statObj.AddStatistic(mitk::ImageStatisticsConstants::UPP(), imageStatisticsFilter->GetUPP(*it)); - statObj.m_Histogram = imageStatisticsFilter->GetHistogram(*it).GetPointer(); - - m_StatisticContainers.at(*it)->SetStatisticsForTimeStep(timeStep, statObj); - ++it; - } + statObj.AddStatistic(mitk::ImageStatisticsConstants::MINIMUMPOSITION(), minIndex); + statObj.AddStatistic(mitk::ImageStatisticsConstants::MAXIMUMPOSITION(), maxIndex); + + assert(std::abs(minMaxFilter->GetMax(*it) - imageStatisticsFilter->GetMaximum(*it)) < mitk::eps); + assert(std::abs(minMaxFilter->GetMin(*it) - imageStatisticsFilter->GetMinimum(*it)) < mitk::eps); + + auto voxelVolume = GetVoxelVolume(image); + auto numberOfVoxels = + static_cast(imageStatisticsFilter->GetSum(*it) / (double)imageStatisticsFilter->GetMean(*it)); + auto volume = static_cast(numberOfVoxels) * voxelVolume; + auto rms = std::sqrt(std::pow(imageStatisticsFilter->GetMean(*it), 2.) + + imageStatisticsFilter->GetVariance(*it)); // variance = sigma^2 + auto variance = imageStatisticsFilter->GetSigma(*it) * imageStatisticsFilter->GetSigma(*it); + + statObj.AddStatistic(mitk::ImageStatisticsConstants::NUMBEROFVOXELS(), numberOfVoxels); + statObj.AddStatistic(mitk::ImageStatisticsConstants::VOLUME(), volume); + statObj.AddStatistic(mitk::ImageStatisticsConstants::MEAN(), imageStatisticsFilter->GetMean(*it)); + statObj.AddStatistic(mitk::ImageStatisticsConstants::MINIMUM(), imageStatisticsFilter->GetMinimum(*it)); + statObj.AddStatistic(mitk::ImageStatisticsConstants::MAXIMUM(), imageStatisticsFilter->GetMaximum(*it)); + statObj.AddStatistic(mitk::ImageStatisticsConstants::STANDARDDEVIATION(), imageStatisticsFilter->GetSigma(*it)); + statObj.AddStatistic(mitk::ImageStatisticsConstants::VARIANCE(), variance); + statObj.AddStatistic(mitk::ImageStatisticsConstants::SKEWNESS(), imageStatisticsFilter->GetSkewness(*it)); + statObj.AddStatistic(mitk::ImageStatisticsConstants::KURTOSIS(), imageStatisticsFilter->GetKurtosis(*it)); + statObj.AddStatistic(mitk::ImageStatisticsConstants::RMS(), rms); + statObj.AddStatistic(mitk::ImageStatisticsConstants::MPP(), imageStatisticsFilter->GetMPP(*it)); + statObj.AddStatistic(mitk::ImageStatisticsConstants::ENTROPY(), imageStatisticsFilter->GetEntropy(*it)); + statObj.AddStatistic(mitk::ImageStatisticsConstants::MEDIAN(), imageStatisticsFilter->GetMedian(*it)); + statObj.AddStatistic(mitk::ImageStatisticsConstants::UNIFORMITY(), imageStatisticsFilter->GetUniformity(*it)); + statObj.AddStatistic(mitk::ImageStatisticsConstants::UPP(), imageStatisticsFilter->GetUPP(*it)); + statObj.m_Histogram = imageStatisticsFilter->GetHistogram(*it).GetPointer(); + + statisticContainerForLabelImage->SetStatisticsForTimeStep(timeStep, statObj); + ++it; + } - // swap maskGenerators back - if (swapMasks) - { - m_SecondaryMask = m_InternalMask; - m_InternalMask = nullptr; - } + // swap maskGenerators back + if (swapMasks) + { + m_SecondaryMask = m_InternalMask; + m_InternalMask = nullptr; } + } - bool ImageStatisticsCalculator::IsUpdateRequired(unsigned int timeStep) const + bool ImageStatisticsCalculator::IsUpdateRequired(LabelIndex label) const + { + unsigned long thisClassTimeStamp = this->GetMTime(); + unsigned long inputImageTimeStamp = m_Image->GetMTime(); + + auto it = m_StatisticContainers.find(label); + if (it == m_StatisticContainers.end()) { - unsigned long thisClassTimeStamp = this->GetMTime(); - unsigned long inputImageTimeStamp = m_Image->GetMTime(); - unsigned long statisticsTimeStamp = m_StatisticsUpdateTimePerTimeStep[timeStep]; + return true; + } - if (thisClassTimeStamp > statisticsTimeStamp) // inputs have changed - { - return true; - } + unsigned long statisticsTimeStamp = it->second->GetMTime(); - if (inputImageTimeStamp > statisticsTimeStamp) // image has changed - { - return true; - } + if (thisClassTimeStamp > statisticsTimeStamp) // inputs have changed + { + return true; + } - if (m_MaskGenerator.IsNotNull()) - { - unsigned long maskGeneratorTimeStamp = m_MaskGenerator->GetMTime(); - if (maskGeneratorTimeStamp > statisticsTimeStamp) // there is a mask generator and it has changed - { - return true; - } - } + if (inputImageTimeStamp > statisticsTimeStamp) // image has changed + { + return true; + } - if (m_SecondaryMaskGenerator.IsNotNull()) - { - unsigned long maskGeneratorTimeStamp = m_SecondaryMaskGenerator->GetMTime(); - if (maskGeneratorTimeStamp > statisticsTimeStamp) // there is a secondary mask generator and it has changed - { - return true; - } - } + if (m_MaskGenerator.IsNotNull()) + { + unsigned long maskGeneratorTimeStamp = m_MaskGenerator->GetMTime(); + if (maskGeneratorTimeStamp > statisticsTimeStamp) // there is a mask generator and it has changed + { + return true; + } + } - return false; + if (m_SecondaryMaskGenerator.IsNotNull()) + { + unsigned long maskGeneratorTimeStamp = m_SecondaryMaskGenerator->GetMTime(); + if (maskGeneratorTimeStamp > statisticsTimeStamp) // there is a secondary mask generator and it has changed + { + return true; + } } -} + + return false; + } +} // namespace mitk diff --git a/Modules/ImageStatistics/mitkImageStatisticsCalculator.h b/Modules/ImageStatistics/mitkImageStatisticsCalculator.h index 1db88d9592..3ac08ea8ac 100644 --- a/Modules/ImageStatistics/mitkImageStatisticsCalculator.h +++ b/Modules/ImageStatistics/mitkImageStatisticsCalculator.h @@ -1,128 +1,128 @@ /*=================================================================== 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 #define MITKIMAGESTATISTICSCALCULATOR #include #include #include #include #include #include namespace mitk { class MITKIMAGESTATISTICS_EXPORT ImageStatisticsCalculator: public itk::Object { public: /** Standard Self typedef */ typedef ImageStatisticsCalculator Self; typedef itk::Object Superclass; typedef itk::SmartPointer< Self > Pointer; typedef itk::SmartPointer< const Self > ConstPointer; /** Method for creation through the object factory. */ itkNewMacro(Self) /** Runtime information support. */ itkTypeMacro(ImageStatisticsCalculator_v2, itk::Object) typedef double statisticsValueType; typedef std::map statisticsMapType; typedef itk::Statistics::Histogram HistogramType; typedef unsigned short MaskPixelType; + using LabelIndex = ImageStatisticsContainer::LabelIndex; /**Documentation @brief Set the image for which the statistics are to be computed.*/ void SetInputImage(mitk::Image::Pointer image); /**Documentation @brief Set the mask generator that creates the mask which is to be used to calculate statistics. If no more mask is desired simply set @param mask to nullptr*/ void SetMask(mitk::MaskGenerator::Pointer mask); /**Documentation @brief Set this if more than one mask should be applied (for instance if a IgnorePixelValueMask were to be used alongside with a segmentation). Both masks are combined using pixel wise AND operation. The secondary mask does not have to be the same size than the primary but they need to have some overlap*/ void SetSecondaryMask(mitk::MaskGenerator::Pointer mask); /**Documentation @brief Set number of bins to be used for histogram statistics. If Bin size is set after number of bins, bin size will be used instead!*/ void SetNBinsForHistogramStatistics(unsigned int nBins); /**Documentation @brief Retrieve the number of bins used for histogram statistics. Careful: The return value does not indicate whether NBins or BinSize is used. That solely depends on which parameter has been set last.*/ unsigned int GetNBinsForHistogramStatistics() const; /**Documentation @brief Set bin size to be used for histogram statistics. If nbins is set after bin size, nbins will be used instead!*/ void SetBinSizeForHistogramStatistics(double binSize); /**Documentation @brief Retrieve the bin size for histogram statistics. Careful: The return value does not indicate whether NBins or BinSize is used. That solely depends on which parameter has been set last.*/ double GetBinSizeForHistogramStatistics() const; /**Documentation - @brief Returns the statistics for label @a label and timeStep @a timeStep. If these requested statistics are not computed yet the computation is done as well. + @brief Returns the statistics for label @a label. If these requested statistics are not computed yet the computation is done as well. For performance reasons, statistics for all labels in the image are computed at once. */ - ImageStatisticsContainer::StatisticsObject GetStatistics(unsigned int timeStep=0, unsigned int label=1); + ImageStatisticsContainer::Pointer GetStatistics(LabelIndex label=1); protected: ImageStatisticsCalculator(){ m_nBinsForHistogramStatistics = 100; m_binSizeForHistogramStatistics = 10; m_UseBinSizeOverNBins = false; }; private: //Calculates statistics for each timestep for image template < typename TPixel, unsigned int VImageDimension > void InternalCalculateStatisticsUnmasked( typename itk::Image< TPixel, VImageDimension >* image, TimeGeometry* timeGeometry, TimeStepType timeStep); template < typename TPixel, unsigned int VImageDimension > void InternalCalculateStatisticsMasked( typename itk::Image< TPixel, VImageDimension >* image, TimeGeometry* timeGeometry, unsigned int timeStep); template < typename TPixel, unsigned int VImageDimension > double GetVoxelVolume(typename itk::Image* image) const; - bool IsUpdateRequired(unsigned int timeStep) const; + bool IsUpdateRequired(LabelIndex label) const; mitk::Image::Pointer m_Image; mitk::Image::Pointer m_ImageTimeSlice; mitk::Image::Pointer m_InternalImageForStatistics; mitk::MaskGenerator::Pointer m_MaskGenerator; mitk::Image::Pointer m_InternalMask; mitk::MaskGenerator::Pointer m_SecondaryMaskGenerator; mitk::Image::Pointer m_SecondaryMask; unsigned int m_nBinsForHistogramStatistics; double m_binSizeForHistogramStatistics; bool m_UseBinSizeOverNBins; - std::vector m_StatisticContainers; - std::vector m_StatisticsUpdateTimePerTimeStep; + std::map m_StatisticContainers; }; } #endif // MITKIMAGESTATISTICSCALCULATOR diff --git a/Modules/ImageStatistics/mitkImageStatisticsContainer.cpp b/Modules/ImageStatistics/mitkImageStatisticsContainer.cpp index d49475aa02..f636fc1f64 100644 --- a/Modules/ImageStatistics/mitkImageStatisticsContainer.cpp +++ b/Modules/ImageStatistics/mitkImageStatisticsContainer.cpp @@ -1,208 +1,210 @@ /*=================================================================== 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 #include namespace mitk { ImageStatisticsContainer::ImageStatisticsContainer() { this->Reset(); } //The order is derived from the old (<2018) image statistics plugin. const ImageStatisticsContainer::StatisticsObject::StatisticNameVector ImageStatisticsContainer::StatisticsObject::m_DefaultNames = { ImageStatisticsConstants::MEAN(), ImageStatisticsConstants::MEDIAN(), ImageStatisticsConstants::STANDARDDEVIATION(), ImageStatisticsConstants::RMS(), ImageStatisticsConstants::MAXIMUM(), ImageStatisticsConstants::MAXIMUMPOSITION(), ImageStatisticsConstants::MINIMUM(), ImageStatisticsConstants::MINIMUMPOSITION(), ImageStatisticsConstants::NUMBEROFVOXELS(), ImageStatisticsConstants::VOLUME(), ImageStatisticsConstants::SKEWNESS(), ImageStatisticsConstants::KURTOSIS(), ImageStatisticsConstants::UNIFORMITY(), ImageStatisticsConstants::ENTROPY(), ImageStatisticsConstants::MPP(), ImageStatisticsConstants::UPP() }; ImageStatisticsContainer::StatisticsObject::StatisticsObject() { Reset(); } void ImageStatisticsContainer::StatisticsObject::AddStatistic(const std::string& key, StatisticsVariantType value) { m_Statistics.emplace(key, value); if (std::find(std::cbegin(m_DefaultNames), std::cend(m_DefaultNames), key)==std::cend(m_DefaultNames)) { if (std::find(std::cbegin(m_CustomNames), std::cend(m_CustomNames), key) == std::cend(m_CustomNames)) { m_CustomNames.emplace_back(key); } } } const ImageStatisticsContainer::StatisticsObject::StatisticNameVector& ImageStatisticsContainer::StatisticsObject::GetDefaultStatisticNames() { return m_DefaultNames; }; const ImageStatisticsContainer::StatisticsObject::StatisticNameVector& ImageStatisticsContainer::StatisticsObject::GetCustomStatisticNames() const { return m_CustomNames; } ImageStatisticsContainer::StatisticsObject::StatisticNameVector ImageStatisticsContainer::StatisticsObject::GetAllStatisticNames() const { StatisticNameVector names = GetDefaultStatisticNames(); names.insert(std::end(names),std::cbegin(m_CustomNames), std::cend(m_CustomNames)); return std::move(names); } bool ImageStatisticsContainer::StatisticsObject::HasStatistic(const std::string& name) const { return m_Statistics.find(name)!=std::cend(m_Statistics); } ImageStatisticsContainer::StatisticsVariantType ImageStatisticsContainer::StatisticsObject::GetValueNonConverted(const std::string& name) const { if (HasStatistic(name)) { return m_Statistics.find(name)->second; } else { mitkThrow() << "invalid statistic key, could not find"; } } void ImageStatisticsContainer::StatisticsObject::Reset() { m_Statistics.clear(); m_CustomNames.clear(); } bool ImageStatisticsContainer::TimeStepExists(TimeStepType timeStep) const { return m_TimeStepMap.find(timeStep) != m_TimeStepMap.end(); } const ImageStatisticsContainer::StatisticsObject& ImageStatisticsContainer::GetStatisticsForTimeStep(TimeStepType timeStep) const { auto it = m_TimeStepMap.find(timeStep); if (it != m_TimeStepMap.end()) { return it->second; } mitkThrow() << "StatisticsObject for timeStep " << timeStep << " not found!"; } void ImageStatisticsContainer::SetStatisticsForTimeStep(TimeStepType timeStep, StatisticsObject statistics) { if (timeStep < this->GetTimeSteps()) { m_TimeStepMap.emplace(timeStep, statistics); + this->Modified(); } else { mitkThrow() << "Given timeStep " << timeStep << " out of timeStep geometry bounds. TimeSteps in geometry: " << this->GetTimeSteps(); } } void ImageStatisticsContainer::PrintSelf(std::ostream &os, itk::Indent indent) const { Superclass::PrintSelf(os, indent); for (unsigned int i = 0; i < this->GetTimeSteps(); i++) { auto statisticsValues = GetStatisticsForTimeStep(i); os << std::endl << indent << "Statistics instance for timeStep " << i << ":"; auto statisticKeys = statisticsValues.GetAllStatisticNames(); + os << std::endl << indent << "Number of entries: " << statisticKeys.size(); for (const auto& aKey : statisticKeys) { os << std::endl << indent.GetNextIndent() << aKey << ": " << statisticsValues.GetValueNonConverted(aKey); } } } unsigned int ImageStatisticsContainer::GetNumberOfTimeSteps()const { return this->GetTimeSteps(); } void ImageStatisticsContainer::Reset() { for (auto iter = m_TimeStepMap.begin(); iter != m_TimeStepMap.end(); iter++) { iter->second.Reset(); } } itk::LightObject::Pointer ImageStatisticsContainer::InternalClone() const { itk::LightObject::Pointer ioPtr = Superclass::InternalClone(); Self::Pointer rval = dynamic_cast(ioPtr.GetPointer()); if (rval.IsNull()) { itkExceptionMacro(<< "downcast to type " << "StatisticsContainer" << " failed."); } rval->SetTimeStepMap(m_TimeStepMap); rval->SetTimeGeometry(this->GetTimeGeometry()->Clone()); return ioPtr; } void ImageStatisticsContainer::SetTimeStepMap(TimeStepMapType map) { m_TimeStepMap = map; } ImageStatisticsContainer::StatisticsObject::StatisticNameVector GetAllStatisticNames(const ImageStatisticsContainer* container) { ImageStatisticsContainer::StatisticsObject::StatisticNameVector names = ImageStatisticsContainer::StatisticsObject::GetDefaultStatisticNames(); if (container) { std::set customKeys; for (unsigned int i = 0; i < container->GetTimeSteps(); i++) { auto statisticKeys = container->GetStatisticsForTimeStep(i).GetCustomStatisticNames(); customKeys.insert(std::cbegin(statisticKeys), std::cend(statisticKeys)); } names.insert(std::end(names), std::begin(customKeys), std::end(customKeys)); } return std::move(names); }; ImageStatisticsContainer::StatisticsObject::StatisticNameVector GetAllStatisticNames(std::vector containers) { ImageStatisticsContainer::StatisticsObject::StatisticNameVector names = ImageStatisticsContainer::StatisticsObject::GetDefaultStatisticNames(); std::set customKeys; for (auto container : containers) { for (unsigned int i = 0; i < container->GetTimeSteps(); i++) { auto statisticKeys = container->GetStatisticsForTimeStep(i).GetCustomStatisticNames(); customKeys.insert(std::cbegin(statisticKeys), std::cend(statisticKeys)); } } names.insert(std::end(names), std::begin(customKeys), std::end(customKeys)); return std::move(names); }; } diff --git a/Modules/ImageStatistics/mitkImageStatisticsContainer.h b/Modules/ImageStatistics/mitkImageStatisticsContainer.h index 9a8542cd86..591928490b 100644 --- a/Modules/ImageStatistics/mitkImageStatisticsContainer.h +++ b/Modules/ImageStatistics/mitkImageStatisticsContainer.h @@ -1,165 +1,163 @@ /*=================================================================== 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 MITKIMAGESTATISTICSCONTAINER #define MITKIMAGESTATISTICSCONTAINER #include #include #include #include #include namespace mitk { /** @brief Container class for storing a StatisticsObject for each timestep. Stored statistics are: - for the defined statistics, see GetAllStatisticNames - - Label (if applicable, the label (unsigned short) of the mask the statistics belong to) - Histogram of Pixel Values*/ class MITKIMAGESTATISTICS_EXPORT ImageStatisticsContainer : public mitk::BaseData { public: mitkClassMacro(ImageStatisticsContainer, mitk::BaseData) itkFactorylessNewMacro(Self) itkCloneMacro(Self) using HistogramType = itk::Statistics::Histogram; using RealType = double; using LabelIndex = unsigned int; using IndexType = vnl_vector; using VoxelCountType = unsigned long; using StatisticsVariantType = boost::variant; using StatisticsMapType = std::map < std::string, StatisticsVariantType>; using StatisticsKeyType = std::string; virtual void SetRequestedRegionToLargestPossibleRegion() override {}; virtual bool RequestedRegionIsOutsideOfTheBufferedRegion() override { return false; }; virtual bool VerifyRequestedRegion() override { return true; }; virtual void SetRequestedRegion(const itk::DataObject*) override {}; /** @brief Container class for storing the computed image statistics. @details The statistics are stored in a map with value as boost::variant. The type used to create the boost::variant is important as only this type can be recovered lateron. */ class MITKIMAGESTATISTICS_EXPORT StatisticsObject { public: StatisticsObject(); /** @brief Adds a statistic to the statistics object @details if already a statistic with that name is included, it is overwritten */ void AddStatistic(const std::string& key, StatisticsVariantType value); using StatisticNameVector = std::vector; /** @brief Returns the names of the default statistics @details The order is derived from the image statistics plugin. */ static const StatisticNameVector& GetDefaultStatisticNames(); /** @brief Returns the names of all custom statistics (defined at runtime and no default names). */ const StatisticNameVector& GetCustomStatisticNames() const; /** @brief Returns the names of all statistics (default and custom defined) Additional custom keys are added at the end in a sorted order. */ StatisticNameVector GetAllStatisticNames() const; bool HasStatistic(const std::string& name) const; /** @brief Converts the requested value to the defined type @param name defined string on creation (AddStatistic) @exception if no statistics with key name was found. */ template TType GetValueConverted(const std::string& name) const { auto value = GetValueNonConverted(name); return boost::get(value); }; /** @brief Returns the requested value @exception if no statistics with key name was found. */ StatisticsVariantType GetValueNonConverted(const std::string& name) const; void Reset(); - LabelIndex m_Label; HistogramType::ConstPointer m_Histogram=nullptr; private: StatisticsMapType m_Statistics; StatisticNameVector m_CustomNames; static const StatisticNameVector m_DefaultNames; }; - typedef std::map TimeStepMapType; + using TimeStepMapType = std::map; unsigned int GetNumberOfTimeSteps() const; /** @brief Deletes all stored values*/ void Reset(); /** @brief Returns the statisticObject for the given Timestep @pre timeStep must be valid */ const StatisticsObject& GetStatisticsForTimeStep(TimeStepType timeStep) const; /** @brief Sets the statisticObject for the given Timestep @pre timeStep must be valid */ void SetStatisticsForTimeStep(TimeStepType timeStep, StatisticsObject statistics); /** @brief Checks if the Time step exists @pre timeStep must be valid */ bool TimeStepExists(TimeStepType timeStep) const; protected: ImageStatisticsContainer(); virtual void PrintSelf(std::ostream &os, itk::Indent indent) const override; private: itk::LightObject::Pointer InternalClone() const override; void SetTimeStepMap(TimeStepMapType map); TimeStepMapType m_TimeStepMap; }; MITKIMAGESTATISTICS_EXPORT ImageStatisticsContainer::StatisticsObject::StatisticNameVector GetAllStatisticNames(const ImageStatisticsContainer* container); MITKIMAGESTATISTICS_EXPORT ImageStatisticsContainer::StatisticsObject::StatisticNameVector GetAllStatisticNames(std::vector containers); } #endif // MITKIMAGESTATISTICSCONTAINER diff --git a/Modules/ImageStatisticsUI/Qmitk/QmitkImageStatisticsCalculationJob.cpp b/Modules/ImageStatisticsUI/Qmitk/QmitkImageStatisticsCalculationJob.cpp index 52d62643a9..ec242732b9 100644 --- a/Modules/ImageStatisticsUI/Qmitk/QmitkImageStatisticsCalculationJob.cpp +++ b/Modules/ImageStatisticsUI/Qmitk/QmitkImageStatisticsCalculationJob.cpp @@ -1,253 +1,246 @@ /*=================================================================== 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 "QmitkImageStatisticsCalculationJob.h" #include "mitkImageStatisticsCalculator.h" #include #include #include QmitkImageStatisticsCalculationJob::QmitkImageStatisticsCalculationJob() : QThread() , m_StatisticsImage(nullptr) , m_BinaryMask(nullptr) , m_PlanarFigureMask(nullptr) , m_TimeStep(0) , m_IgnoreZeros(false) , m_HistogramNBins(100) , m_StatisticChanged(false) , m_CalculationSuccessful(false) { } QmitkImageStatisticsCalculationJob::~QmitkImageStatisticsCalculationJob() { } void QmitkImageStatisticsCalculationJob::Initialize( mitk::Image::ConstPointer image, mitk::Image::ConstPointer binaryImage, mitk::PlanarFigure::ConstPointer planarFig ) { // reset old values if( this->m_StatisticsImage.IsNotNull() ) this->m_StatisticsImage = nullptr; if( this->m_BinaryMask.IsNotNull() ) this->m_BinaryMask = nullptr; if( this->m_PlanarFigureMask.IsNotNull()) this->m_PlanarFigureMask = nullptr; // set new values if passed in if(image.IsNotNull()) this->m_StatisticsImage = image; if(binaryImage.IsNotNull()) this->m_BinaryMask = binaryImage; if(planarFig.IsNotNull()) this->m_PlanarFigureMask = planarFig; } void QmitkImageStatisticsCalculationJob::SetTimeStep( int times ) { this->m_TimeStep = times; } int QmitkImageStatisticsCalculationJob::GetTimeStep() const { return this->m_TimeStep; } mitk::ImageStatisticsContainer::ConstPointer QmitkImageStatisticsCalculationJob::GetStatisticsData() const { mitk::ImageStatisticsContainer::ConstPointer constContainer = this->m_StatisticsContainer.GetPointer(); return constContainer; } mitk::Image::ConstPointer QmitkImageStatisticsCalculationJob::GetStatisticsImage() const { return this->m_StatisticsImage; } mitk::Image::ConstPointer QmitkImageStatisticsCalculationJob::GetMaskImage() const { return this->m_BinaryMask; } mitk::PlanarFigure::ConstPointer QmitkImageStatisticsCalculationJob::GetPlanarFigure() const { return this->m_PlanarFigureMask; } void QmitkImageStatisticsCalculationJob::SetIgnoreZeroValueVoxel(bool _arg) { this->m_IgnoreZeros = _arg; } bool QmitkImageStatisticsCalculationJob::GetIgnoreZeroValueVoxel() const { return this->m_IgnoreZeros; } void QmitkImageStatisticsCalculationJob::SetHistogramNBins(unsigned int nbins) { this->m_HistogramNBins = nbins; } unsigned int QmitkImageStatisticsCalculationJob::GetHistogramNBins() const { return this->m_HistogramNBins; } std::string QmitkImageStatisticsCalculationJob::GetLastErrorMessage() const { return m_message; } QmitkImageStatisticsCalculationJob::HistogramType::ConstPointer QmitkImageStatisticsCalculationJob::GetTimeStepHistogram(unsigned int t) const { if (t >= this->m_HistogramVector.size()) return nullptr; return this->m_HistogramVector[t]; } bool QmitkImageStatisticsCalculationJob::GetStatisticsChangedFlag() const { return m_StatisticChanged; } bool QmitkImageStatisticsCalculationJob::GetStatisticsUpdateSuccessFlag() const { return m_CalculationSuccessful; } void QmitkImageStatisticsCalculationJob::run() { bool statisticCalculationSuccessful = true; mitk::ImageStatisticsCalculator::Pointer calculator = mitk::ImageStatisticsCalculator::New(); if(this->m_StatisticsImage.IsNotNull()) { calculator->SetInputImage(m_StatisticsImage->Clone()); } else { statisticCalculationSuccessful = false; } // Bug 13416 : The ImageStatistics::SetImageMask() method can throw exceptions, i.e. when the dimensionality // of the masked and input image differ, we need to catch them and mark the calculation as failed // the same holds for the ::SetPlanarFigure() try { if(this->m_BinaryMask.IsNotNull()) { mitk::ImageMaskGenerator::Pointer imgMask = mitk::ImageMaskGenerator::New(); imgMask->SetImageMask(m_BinaryMask->Clone()); calculator->SetMask(imgMask.GetPointer()); } if(this->m_PlanarFigureMask.IsNotNull()) { mitk::PlanarFigureMaskGenerator::Pointer pfMaskGen = mitk::PlanarFigureMaskGenerator::New(); pfMaskGen->SetInputImage(m_StatisticsImage->Clone()); pfMaskGen->SetPlanarFigure(m_PlanarFigureMask->Clone()); calculator->SetMask(pfMaskGen.GetPointer()); } } catch (const mitk::Exception& e) { MITK_ERROR << "MITK Exception: " << e.what(); m_message = e.what(); statisticCalculationSuccessful = false; } catch( const itk::ExceptionObject& e) { MITK_ERROR << "ITK Exception:" << e.what(); m_message = e.what(); statisticCalculationSuccessful = false; } catch ( const std::runtime_error &e ) { MITK_ERROR<< "Runtime Exception: " << e.what(); m_message = e.what(); statisticCalculationSuccessful = false; } catch ( const std::exception &e ) { MITK_ERROR<< "Standard Exception: " << e.what(); m_message = e.what(); statisticCalculationSuccessful = false; } bool statisticChanged = false; if (this->m_IgnoreZeros) { mitk::IgnorePixelMaskGenerator::Pointer ignorePixelValueMaskGen = mitk::IgnorePixelMaskGenerator::New(); ignorePixelValueMaskGen->SetIgnoredPixelValue(0); ignorePixelValueMaskGen->SetInputImage(m_StatisticsImage->Clone()); calculator->SetSecondaryMask(ignorePixelValueMaskGen.GetPointer()); } else { calculator->SetSecondaryMask(nullptr); } calculator->SetNBinsForHistogramStatistics(m_HistogramNBins); - for (unsigned int i = 0; i < m_StatisticsImage->GetTimeSteps(); i++) - { try { - calculator->GetStatistics(i); + calculator->GetStatistics(); } catch ( mitk::Exception& e) { m_message = e.GetDescription(); MITK_ERROR<< "MITK Exception: " << e.what(); statisticCalculationSuccessful = false; } catch ( const std::runtime_error &e ) { m_message = "Failure: " + std::string(e.what()); MITK_ERROR<< "Runtime Exception: " << e.what(); statisticCalculationSuccessful = false; } catch ( const std::exception &e ) { m_message = "Failure: " + std::string(e.what()); MITK_ERROR<< "Standard Exception: " << e.what(); statisticCalculationSuccessful = false; } - } this->m_StatisticChanged = statisticChanged; this->m_CalculationSuccessful = statisticCalculationSuccessful; if(statisticCalculationSuccessful) { - // maybe use reset if this line does not work - m_StatisticsContainer = mitk::ImageStatisticsContainer::New(); + m_StatisticsContainer = calculator->GetStatistics(); this->m_HistogramVector.clear(); - mitk::TimeGeometry::Pointer tg = m_StatisticsImage->GetTimeGeometry()->Clone(); - m_StatisticsContainer->SetTimeGeometry(tg); for (unsigned int i = 0; i < m_StatisticsImage->GetTimeSteps(); i++) { - this->m_StatisticsContainer->SetStatisticsForTimeStep(i, calculator->GetStatistics(i)); - this->m_HistogramVector.push_back(calculator->GetStatistics(i).m_Histogram); + this->m_HistogramVector.push_back(calculator->GetStatistics()->GetStatisticsForTimeStep(i).m_Histogram); } } } diff --git a/Modules/Segmentation/Algorithms/mitkFeatureBasedEdgeDetectionFilter.cpp b/Modules/Segmentation/Algorithms/mitkFeatureBasedEdgeDetectionFilter.cpp index bef734bbdb..34dfe97c86 100644 --- a/Modules/Segmentation/Algorithms/mitkFeatureBasedEdgeDetectionFilter.cpp +++ b/Modules/Segmentation/Algorithms/mitkFeatureBasedEdgeDetectionFilter.cpp @@ -1,199 +1,199 @@ /*=================================================================== 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 "mitkFeatureBasedEdgeDetectionFilter.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include mitk::FeatureBasedEdgeDetectionFilter::FeatureBasedEdgeDetectionFilter() { this->SetNumberOfRequiredInputs(1); this->SetNumberOfIndexedOutputs(1); } mitk::FeatureBasedEdgeDetectionFilter::~FeatureBasedEdgeDetectionFilter() { } void mitk::FeatureBasedEdgeDetectionFilter::GenerateData() { mitk::Image::Pointer image = ImageToUnstructuredGridFilter::GetInput(); if (m_SegmentationMask.IsNull()) { MITK_WARN << "Please set a segmentation mask first" << std::endl; return; } // First create a threshold segmentation of the image. The threshold is determined // by the mean +/- stddev of the pixel values that are covered by the segmentation mask // Compute mean and stdDev based on the current segmentation mitk::ImageStatisticsCalculator::Pointer statCalc = mitk::ImageStatisticsCalculator::New(); statCalc->SetInputImage(image); mitk::ImageMaskGenerator::Pointer imgMask = mitk::ImageMaskGenerator::New(); imgMask->SetImageMask(m_SegmentationMask); - auto stats = statCalc->GetStatistics(); + auto stats = statCalc->GetStatistics()->GetStatisticsForTimeStep(0); double mean = stats.GetValueConverted(mitk::ImageStatisticsConstants::MEAN()); double stdDev = stats.GetValueConverted(mitk::ImageStatisticsConstants::STANDARDDEVIATION()); double upperThreshold = mean + stdDev; double lowerThreshold = mean - stdDev; // Perform thresholding mitk::Image::Pointer thresholdImage = mitk::Image::New(); AccessByItk_3(image.GetPointer(), ITKThresholding, lowerThreshold, upperThreshold, thresholdImage) mitk::ProgressBar::GetInstance() ->Progress(2); // Postprocess threshold segmentation // First a closing will be executed mitk::Image::Pointer closedImage = mitk::Image::New(); AccessByItk_1(thresholdImage, ThreadedClosing, closedImage); // Then we will holes that might exist mitk::MorphologicalOperations::FillHoles(closedImage); mitk::ProgressBar::GetInstance()->Progress(); // Extract the binary edges of the resulting segmentation mitk::Image::Pointer edgeImage = mitk::Image::New(); AccessByItk_1(closedImage, ContourSearch, edgeImage); // Convert the edge image into an unstructured grid mitk::ImageToUnstructuredGridFilter::Pointer i2UFilter = mitk::ImageToUnstructuredGridFilter::New(); i2UFilter->SetInput(edgeImage); i2UFilter->SetThreshold(1.0); i2UFilter->Update(); m_PointGrid = this->GetOutput(); if (m_PointGrid.IsNull()) m_PointGrid = mitk::UnstructuredGrid::New(); m_PointGrid->SetVtkUnstructuredGrid(i2UFilter->GetOutput()->GetVtkUnstructuredGrid()); mitk::ProgressBar::GetInstance()->Progress(); } template void mitk::FeatureBasedEdgeDetectionFilter::ThreadedClosing(itk::Image *originalImage, mitk::Image::Pointer &result) { typedef itk::BinaryBallStructuringElement myKernelType; myKernelType ball; ball.SetRadius(1); ball.CreateStructuringElement(); typedef typename itk::Image ImageType; typename itk::DilateObjectMorphologyImageFilter::Pointer dilationFilter = itk::DilateObjectMorphologyImageFilter::New(); dilationFilter->SetInput(originalImage); dilationFilter->SetKernel(ball); dilationFilter->Update(); typename itk::Image::Pointer dilatedImage = dilationFilter->GetOutput(); typename itk::ErodeObjectMorphologyImageFilter::Pointer erodeFilter = itk::ErodeObjectMorphologyImageFilter::New(); erodeFilter->SetInput(dilatedImage); erodeFilter->SetKernel(ball); erodeFilter->Update(); mitk::GrabItkImageMemory(erodeFilter->GetOutput(), result); } template void mitk::FeatureBasedEdgeDetectionFilter::ContourSearch(itk::Image *originalImage, mitk::Image::Pointer &result) { typedef itk::Image ImageType; typedef itk::BinaryContourImageFilter binaryContourImageFilterType; typename binaryContourImageFilterType::Pointer binaryContourFilter = binaryContourImageFilterType::New(); binaryContourFilter->SetInput(originalImage); binaryContourFilter->SetForegroundValue(1); binaryContourFilter->SetBackgroundValue(0); binaryContourFilter->Update(); typename itk::Image::Pointer itkImage = itk::Image::New(); itkImage->Graft(binaryContourFilter->GetOutput()); mitk::GrabItkImageMemory(itkImage, result); } template void mitk::FeatureBasedEdgeDetectionFilter::ITKThresholding(itk::Image *originalImage, double lower, double upper, mitk::Image::Pointer &result) { typedef itk::Image ImageType; typedef itk::Image SegmentationType; typedef itk::BinaryThresholdImageFilter ThresholdFilterType; if (typeid(TPixel) != typeid(float) && typeid(TPixel) != typeid(double)) { // round the thresholds if we have nor a float or double image lower = std::floor(lower + 0.5); upper = std::floor(upper - 0.5); } if (lower >= upper) { upper = lower; } typename ThresholdFilterType::Pointer filter = ThresholdFilterType::New(); filter->SetInput(originalImage); filter->SetLowerThreshold(lower); filter->SetUpperThreshold(upper); filter->SetInsideValue(1); filter->SetOutsideValue(0); filter->Update(); mitk::GrabItkImageMemory(filter->GetOutput(), result); } void mitk::FeatureBasedEdgeDetectionFilter::SetSegmentationMask(mitk::Image::Pointer segmentation) { this->m_SegmentationMask = segmentation; } void mitk::FeatureBasedEdgeDetectionFilter::GenerateOutputInformation() { Superclass::GenerateOutputInformation(); } diff --git a/Modules/SurfaceInterpolation/mitkImageToPointCloudFilter.cpp b/Modules/SurfaceInterpolation/mitkImageToPointCloudFilter.cpp index 9878911009..6b001b3410 100644 --- a/Modules/SurfaceInterpolation/mitkImageToPointCloudFilter.cpp +++ b/Modules/SurfaceInterpolation/mitkImageToPointCloudFilter.cpp @@ -1,167 +1,167 @@ /*=================================================================== 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 "mitkImageToPointCloudFilter.h" #include #include #include #include #include #include #include #include #include #include mitk::ImageToPointCloudFilter::ImageToPointCloudFilter() { m_Method = DetectionMethod(0); this->SetNumberOfRequiredInputs(1); this->SetNumberOfIndexedOutputs(1); } mitk::ImageToPointCloudFilter::~ImageToPointCloudFilter() { } void mitk::ImageToPointCloudFilter::GenerateData() { mitk::Image::ConstPointer image = ImageToUnstructuredGridFilter::GetInput(); m_Geometry = image->GetGeometry(); if (image.IsNull()) { MITK_ERROR << "mitk::ImageToContourFilter: No input available. " "Please set the input!" << std::endl; return; } mitk::Image::Pointer notConstImage = const_cast(image.GetPointer()); switch (m_Method) { case 0: AccessByItk_1(notConstImage.GetPointer(), StdDeviations, 2) break; case 1: AccessByItk_1(notConstImage.GetPointer(), StdDeviations, 3) break; case 2: AccessByItk_1(notConstImage.GetPointer(), StdDeviations, 4) break; default: AccessByItk_1(notConstImage.GetPointer(), StdDeviations, 2) break; } } template void mitk::ImageToPointCloudFilter::StdDeviations(itk::Image *image, int amount) { typedef itk::Image InputImageType; typedef itk::CastImageFilter ImagePTypeToFloatPTypeCasterType; typedef itk::LaplacianImageFilter LaplacianFilterType; typename LaplacianFilterType::Pointer lapFilter = LaplacianFilterType::New(); typename ImagePTypeToFloatPTypeCasterType::Pointer caster = ImagePTypeToFloatPTypeCasterType::New(); caster->SetInput(image); caster->Update(); FloatImageType::Pointer fImage = caster->GetOutput(); lapFilter->SetInput(fImage); lapFilter->UpdateLargestPossibleRegion(); mitk::Image::Pointer edgeImage = mitk::ImportItkImage(lapFilter->GetOutput()); mitk::ImageStatisticsCalculator::Pointer statCalc = mitk::ImageStatisticsCalculator::New(); statCalc->SetInputImage(edgeImage); - auto stats = statCalc->GetStatistics(); + auto stats = statCalc->GetStatistics()->GetStatisticsForTimeStep(0); auto mean = stats.GetValueConverted(mitk::ImageStatisticsConstants::MEAN()); auto stdDev = stats.GetValueConverted(mitk::ImageStatisticsConstants::STANDARDDEVIATION()); double upperThreshold = mean + stdDev * amount; double lowerThreshold = mean - stdDev * amount; typename itk::ImageRegionIterator it(lapFilter->GetOutput(), lapFilter->GetOutput()->GetRequestedRegion()); vtkSmartPointer points = vtkSmartPointer::New(); double greatX = 0, greatY = 0, greatZ = 0; it.GoToBegin(); while (!it.IsAtEnd()) { if (it.Get() > lowerThreshold && it.Get() < upperThreshold) { it.Set(0); } else { it.Set(1); mitk::Point3D imagePoint; mitk::Point3D worldPoint; imagePoint[0] = it.GetIndex()[0]; imagePoint[1] = it.GetIndex()[1]; imagePoint[2] = it.GetIndex()[2]; m_Geometry->IndexToWorld(imagePoint, worldPoint); if (worldPoint[0] > greatX) greatX = worldPoint[0]; if (worldPoint[1] > greatY) greatY = worldPoint[1]; if (worldPoint[2] > greatZ) greatZ = worldPoint[2]; points->InsertNextPoint(worldPoint[0], worldPoint[1], worldPoint[2]); m_NumberOfExtractedPoints++; } ++it; } /*need to build the UnstructuredGrid with at least one vertex otherwise its not visible*/ vtkSmartPointer verts = vtkSmartPointer::New(); verts->GetPointIds()->SetNumberOfIds(m_NumberOfExtractedPoints); for (int i = 0; i < m_NumberOfExtractedPoints; i++) { verts->GetPointIds()->SetId(i, i); } vtkSmartPointer uGrid = vtkSmartPointer::New(); uGrid->Allocate(1); uGrid->InsertNextCell(verts->GetCellType(), verts->GetPointIds()); uGrid->SetPoints(points); mitk::UnstructuredGrid::Pointer outputGrid = mitk::UnstructuredGrid::New(); outputGrid->SetVtkUnstructuredGrid(uGrid); this->SetNthOutput(0, outputGrid); } void mitk::ImageToPointCloudFilter::GenerateOutputInformation() { Superclass::GenerateOutputInformation(); }