diff --git a/Modules/ImageStatistics/Testing/mitkImageStatisticsCalculatorTest.cpp b/Modules/ImageStatistics/Testing/mitkImageStatisticsCalculatorTest.cpp index ca883e8895..7a3b8642b4 100644 --- a/Modules/ImageStatistics/Testing/mitkImageStatisticsCalculatorTest.cpp +++ b/Modules/ImageStatistics/Testing/mitkImageStatisticsCalculatorTest.cpp @@ -1,1133 +1,1133 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #include "mitkImageStatisticsCalculator.h" #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(TestPic3DCroppedNoMask); MITK_TEST(TestPic3DCroppedBinMask); MITK_TEST(TestPic3DCroppedMultilabelMask); MITK_TEST(TestPic3DCroppedPlanarFigure); MITK_TEST(TestUS4DCroppedNoMaskTimeStep1); MITK_TEST(TestUS4DCroppedBinMaskTimeStep1); MITK_TEST(TestUS4DCroppedMultilabelMaskTimeStep1); MITK_TEST(TestUS4DCroppedPlanarFigureTimeStep1); MITK_TEST(TestUS4DCroppedAllTimesteps); MITK_TEST(TestUS4DCropped3DMask); CPPUNIT_TEST_SUITE_END(); public: 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 TestPic3DCroppedNoMask(); void TestPic3DCroppedBinMask(); void TestPic3DCroppedMultilabelMask(); void TestPic3DCroppedPlanarFigure(); void TestUS4DCroppedNoMaskTimeStep1(); void TestUS4DCroppedBinMaskTimeStep1(); void TestUS4DCroppedMultilabelMaskTimeStep1(); void TestUS4DCroppedPlanarFigureTimeStep1(); void TestUS4DCroppedAllTimesteps(); void TestUS4DCropped3DMask(); private: mitk::Image::ConstPointer m_TestImage; mitk::Image::ConstPointer m_Pic3DCroppedImage; mitk::Image::Pointer m_Pic3DCroppedBinMask; mitk::Image::Pointer m_Pic3DCroppedMultilabelMask; mitk::PlanarFigure::Pointer m_Pic3DCroppedPlanarFigure; mitk::Image::ConstPointer m_US4DCroppedImage; mitk::Image::Pointer m_US4DCroppedBinMask; mitk::Image::Pointer m_US4DCroppedMultilabelMask; mitk::Image::Pointer m_US4DCropped3DBinMask; mitk::PlanarFigure::Pointer m_US4DCroppedPlanarFigure; mitk::PlaneGeometry::Pointer m_Geometry; // creates a polygon given a geometry and a vector of 2d points mitk::PlanarPolygon::Pointer GeneratePlanarPolygon(mitk::PlaneGeometry::Pointer geometry, std::vector points); // universal function to calculate statistics const mitk::ImageStatisticsContainer::Pointer mitkImageStatisticsCalculatorTestSuite::ComputeStatistics(mitk::Image::ConstPointer image, mitk::MaskGenerator::Pointer maskGen = nullptr, mitk::MaskGenerator::Pointer secondardMaskGen = nullptr) { 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(); } void VerifyStatistics(mitk::ImageStatisticsContainer::ImageStatisticsObject stats, mitk::ImageStatisticsContainer::RealType testMean, mitk::ImageStatisticsContainer::RealType testSD, mitk::ImageStatisticsContainer::RealType testMedian = 0); // T26098 histogram statistics need to be tested (median, uniformity, UPP, entropy) void VerifyStatistics(mitk::ImageStatisticsContainer::ImageStatisticsObject stats, mitk::ImageStatisticsContainer::VoxelCountType N, mitk::ImageStatisticsContainer::RealType mean, mitk::ImageStatisticsContainer::RealType MPP, mitk::ImageStatisticsContainer::RealType skewness, mitk::ImageStatisticsContainer::RealType kurtosis, mitk::ImageStatisticsContainer::RealType variance, mitk::ImageStatisticsContainer::RealType stdev, mitk::ImageStatisticsContainer::RealType min, mitk::ImageStatisticsContainer::RealType max, mitk::ImageStatisticsContainer::RealType RMS, mitk::ImageStatisticsContainer::IndexType minIndex, mitk::ImageStatisticsContainer::IndexType maxIndex); }; void mitkImageStatisticsCalculatorTestSuite::TestUninitializedImage() { /***************************** * loading uninitialized image to datastorage ******************************/ MITK_INFO << std::endl << "Test uninitialized image: -----------------------------------------------------------------------------------"; mitk::Image::Pointer image = mitk::Image::New(); mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData(image); mitk::ImageStatisticsCalculator::Pointer is = mitk::ImageStatisticsCalculator::New(); CPPUNIT_ASSERT_THROW(is->GetStatistics(), mitk::Exception); } void mitkImageStatisticsCalculatorTestSuite::TestCase1() { /***************************** * one whole white pixel * -> mean of 255 expected ******************************/ MITK_INFO << std::endl << "Test case 1:-----------------------------------------------------------------------------------"; std::string filename = this->GetTestDataFilePath("ImageStatisticsTestData/testimage.nrrd"); m_TestImage = mitk::IOUtil::Load(filename); CPPUNIT_ASSERT_MESSAGE("Failed loading an mitk::Image", m_TestImage.IsNotNull()); m_Geometry = m_TestImage->GetSlicedGeometry()->GetPlaneGeometry(0); CPPUNIT_ASSERT_MESSAGE("Failed getting image geometry", m_Geometry.IsNotNull()); mitk::Point2D pnt1; pnt1[0] = 10.5; pnt1[1] = 3.5; mitk::Point2D pnt2; pnt2[0] = 9.5; pnt2[1] = 3.5; mitk::Point2D pnt3; pnt3[0] = 9.5; pnt3[1] = 4.5; mitk::Point2D pnt4; pnt4[0] = 10.5; pnt4[1] = 4.5; std::vector points{ pnt1,pnt2,pnt3,pnt4 }; auto figure = GeneratePlanarPolygon(m_Geometry, points); mitk::ImageStatisticsContainer::Pointer statisticsContainer; mitk::PlanarFigureMaskGenerator::Pointer planFigMaskGen = mitk::PlanarFigureMaskGenerator::New(); planFigMaskGen->SetInputImage(m_TestImage); planFigMaskGen->SetPlanarFigure(figure.GetPointer()); CPPUNIT_ASSERT_NO_THROW(statisticsContainer = ComputeStatistics(m_TestImage, planFigMaskGen.GetPointer())); auto statisticsObjectTimestep0 = statisticsContainer->GetStatistics(1, 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:-----------------------------------------------------------------------------------"; std::string filename = this->GetTestDataFilePath("ImageStatisticsTestData/testimage.nrrd"); m_TestImage = mitk::IOUtil::Load(filename); m_Geometry = m_TestImage->GetSlicedGeometry()->GetPlaneGeometry(0); mitk::Point2D pnt1; pnt1[0] = 10.0; pnt1[1] = 3.5; mitk::Point2D pnt2; pnt2[0] = 9.5; pnt2[1] = 3.5; mitk::Point2D pnt3; pnt3[0] = 9.5; pnt3[1] = 4.5; mitk::Point2D pnt4; pnt4[0] = 10.0; pnt4[1] = 4.5; std::vector points{ pnt1,pnt2,pnt3,pnt4 }; auto figure = GeneratePlanarPolygon(m_Geometry, points); mitk::PlanarFigureMaskGenerator::Pointer planFigMaskGen = mitk::PlanarFigureMaskGenerator::New(); planFigMaskGen->SetInputImage(m_TestImage); planFigMaskGen->SetPlanarFigure(figure.GetPointer()); mitk::ImageStatisticsContainer::Pointer statisticsContainer; CPPUNIT_ASSERT_NO_THROW(statisticsContainer = ComputeStatistics(m_TestImage, planFigMaskGen.GetPointer())); auto statisticsObjectTimestep0 = statisticsContainer->GetStatistics(1,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:-----------------------------------------------------------------------------------"; std::string filename = this->GetTestDataFilePath("ImageStatisticsTestData/testimage.nrrd"); m_TestImage = mitk::IOUtil::Load(filename); m_Geometry = m_TestImage->GetSlicedGeometry()->GetPlaneGeometry(0); mitk::Point2D pnt1; pnt1[0] = 10.5; pnt1[1] = 3.5; mitk::Point2D pnt2; pnt2[0] = 9.5; pnt2[1] = 3.5; mitk::Point2D pnt3; pnt3[0] = 9.5; pnt3[1] = 4.5; std::vector points{ pnt1,pnt2,pnt3 }; auto figure = GeneratePlanarPolygon(m_Geometry, points); mitk::PlanarFigureMaskGenerator::Pointer planFigMaskGen = mitk::PlanarFigureMaskGenerator::New(); planFigMaskGen->SetInputImage(m_TestImage); planFigMaskGen->SetPlanarFigure(figure.GetPointer()); mitk::ImageStatisticsContainer::Pointer statisticsContainer; CPPUNIT_ASSERT_NO_THROW(statisticsContainer = ComputeStatistics(m_TestImage, planFigMaskGen.GetPointer())); auto statisticsObjectTimestep0 = statisticsContainer->GetStatistics(1, 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:-----------------------------------------------------------------------------------"; std::string filename = this->GetTestDataFilePath("ImageStatisticsTestData/testimage.nrrd"); m_TestImage = mitk::IOUtil::Load(filename); m_Geometry = m_TestImage->GetSlicedGeometry()->GetPlaneGeometry(0); mitk::Point2D pnt1; pnt1[0] = 1.1; pnt1[1] = 1.1; mitk::Point2D pnt2; pnt2[0] = 2.0; pnt2[1] = 2.0; mitk::Point2D pnt3; pnt3[0] = 3.0; pnt3[1] = 1.0; mitk::Point2D pnt4; pnt4[0] = 2.0; pnt4[1] = 0.0; std::vector points{ pnt1,pnt2,pnt3,pnt4 }; auto figure = GeneratePlanarPolygon(m_Geometry, points); mitk::PlanarFigureMaskGenerator::Pointer planFigMaskGen = mitk::PlanarFigureMaskGenerator::New(); planFigMaskGen->SetInputImage(m_TestImage); planFigMaskGen->SetPlanarFigure(figure.GetPointer()); mitk::ImageStatisticsContainer::Pointer statisticsContainer; CPPUNIT_ASSERT_NO_THROW(statisticsContainer = ComputeStatistics(m_TestImage, planFigMaskGen.GetPointer())); auto statisticsObjectTimestep0 = statisticsContainer->GetStatistics(1, 0); this->VerifyStatistics(statisticsObjectTimestep0, 191.25, 127.5, 253.72499847412109); } void mitkImageStatisticsCalculatorTestSuite::TestCase5() { /***************************** * whole pixel (white) + half pixel (gray) in x-direction * -> mean of 191.5 expected ******************************/ MITK_INFO << std::endl << "Test case 5:-----------------------------------------------------------------------------------"; std::string filename = this->GetTestDataFilePath("ImageStatisticsTestData/testimage.nrrd"); m_TestImage = mitk::IOUtil::Load(filename); m_Geometry = m_TestImage->GetSlicedGeometry()->GetPlaneGeometry(0); mitk::Point2D pnt1; pnt1[0] = 11.0; pnt1[1] = 3.5; mitk::Point2D pnt2; pnt2[0] = 9.5; pnt2[1] = 3.5; mitk::Point2D pnt3; pnt3[0] = 9.5; pnt3[1] = 4.5; mitk::Point2D pnt4; pnt4[0] = 11.0; pnt4[1] = 4.5; std::vector points{ pnt1,pnt2,pnt3,pnt4 }; auto figure = GeneratePlanarPolygon(m_Geometry, points); mitk::PlanarFigureMaskGenerator::Pointer planFigMaskGen = mitk::PlanarFigureMaskGenerator::New(); planFigMaskGen->SetInputImage(m_TestImage); planFigMaskGen->SetPlanarFigure(figure.GetPointer()); mitk::ImageStatisticsContainer::Pointer statisticsContainer; CPPUNIT_ASSERT_NO_THROW(statisticsContainer = ComputeStatistics(m_TestImage, planFigMaskGen.GetPointer())); auto statisticsObjectTimestep0 = statisticsContainer->GetStatistics(1, 0); this->VerifyStatistics(statisticsObjectTimestep0, 191.50, 89.802561210691536, 128.63499999046327); } 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:-----------------------------------------------------------------------------------"; std::string filename = this->GetTestDataFilePath("ImageStatisticsTestData/testimage.nrrd"); m_TestImage = mitk::IOUtil::Load(filename); m_Geometry = m_TestImage->GetSlicedGeometry()->GetPlaneGeometry(0); mitk::Point2D pnt1; pnt1[0] = 11.0; pnt1[1] = 3.5; mitk::Point2D pnt2; pnt2[0] = 9.25; pnt2[1] = 3.5; mitk::Point2D pnt3; pnt3[0] = 9.25; pnt3[1] = 4.5; mitk::Point2D pnt4; pnt4[0] = 11.0; pnt4[1] = 4.5; std::vector points{ pnt1,pnt2,pnt3,pnt4 }; auto figure = GeneratePlanarPolygon(m_Geometry, points); mitk::PlanarFigureMaskGenerator::Pointer planFigMaskGen = mitk::PlanarFigureMaskGenerator::New(); planFigMaskGen->SetInputImage(m_TestImage); planFigMaskGen->SetPlanarFigure(figure.GetPointer()); mitk::ImageStatisticsContainer::Pointer statisticsContainer; CPPUNIT_ASSERT_NO_THROW(statisticsContainer = ComputeStatistics(m_TestImage, planFigMaskGen.GetPointer())); auto statisticsObjectTimestep0 = statisticsContainer->GetStatistics(1, 0); this->VerifyStatistics(statisticsObjectTimestep0, 191.5, 89.802561210691536, 128.63499999046327); } 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:-----------------------------------------------------------------------------------"; std::string filename = this->GetTestDataFilePath("ImageStatisticsTestData/testimage.nrrd"); m_TestImage = mitk::IOUtil::Load(filename); m_Geometry = m_TestImage->GetSlicedGeometry()->GetPlaneGeometry(0); mitk::Point2D pnt1; pnt1[0] = 11.0; pnt1[1] = 3.5; mitk::Point2D pnt2; pnt2[0] = 9.0; pnt2[1] = 3.5; mitk::Point2D pnt3; pnt3[0] = 9.0; pnt3[1] = 4.0; mitk::Point2D pnt4; pnt4[0] = 11.0; pnt4[1] = 4.0; std::vector points{ pnt1,pnt2,pnt3,pnt4 }; auto figure = GeneratePlanarPolygon(m_Geometry, points); mitk::PlanarFigureMaskGenerator::Pointer planFigMaskGen = mitk::PlanarFigureMaskGenerator::New(); planFigMaskGen->SetInputImage(m_TestImage); planFigMaskGen->SetPlanarFigure(figure.GetPointer()); mitk::ImageStatisticsContainer::Pointer statisticsContainer; CPPUNIT_ASSERT_NO_THROW(statisticsContainer = ComputeStatistics(m_TestImage, planFigMaskGen.GetPointer())); auto statisticsObjectTimestep0 = statisticsContainer->GetStatistics(1, 0); this->VerifyStatistics(statisticsObjectTimestep0, 127.666666666666667, 127.50032679696680, 128.7750015258789); } void mitkImageStatisticsCalculatorTestSuite::TestCase8() { /***************************** * whole pixel (gray) * -> mean of 128 expected ******************************/ MITK_INFO << std::endl << "Test case 8:-----------------------------------------------------------------------------------"; std::string filename = this->GetTestDataFilePath("ImageStatisticsTestData/testimage.nrrd"); m_TestImage = mitk::IOUtil::Load(filename); m_Geometry = m_TestImage->GetSlicedGeometry()->GetPlaneGeometry(0); mitk::Point2D pnt1; pnt1[0] = 11.5; pnt1[1] = 10.5; mitk::Point2D pnt2; pnt2[0] = 11.5; pnt2[1] = 11.5; mitk::Point2D pnt3; pnt3[0] = 12.5; pnt3[1] = 11.5; mitk::Point2D pnt4; pnt4[0] = 12.5; pnt4[1] = 10.5; std::vector points{ pnt1,pnt2,pnt3,pnt4 }; auto figure = GeneratePlanarPolygon(m_Geometry, points); mitk::PlanarFigureMaskGenerator::Pointer planFigMaskGen = mitk::PlanarFigureMaskGenerator::New(); planFigMaskGen->SetInputImage(m_TestImage); planFigMaskGen->SetPlanarFigure(figure.GetPointer()); mitk::ImageStatisticsContainer::Pointer statisticsContainer; CPPUNIT_ASSERT_NO_THROW(statisticsContainer = ComputeStatistics(m_TestImage, planFigMaskGen.GetPointer())); auto statisticsObjectTimestep0 = statisticsContainer->GetStatistics(1, 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:-----------------------------------------------------------------------------------"; std::string filename = this->GetTestDataFilePath("ImageStatisticsTestData/testimage.nrrd"); m_TestImage = mitk::IOUtil::Load(filename); m_Geometry = m_TestImage->GetSlicedGeometry()->GetPlaneGeometry(0); mitk::Point2D pnt1; pnt1[0] = 11.5; pnt1[1] = 10.5; mitk::Point2D pnt2; pnt2[0] = 11.5; pnt2[1] = 12.0; mitk::Point2D pnt3; pnt3[0] = 12.5; pnt3[1] = 12.0; mitk::Point2D pnt4; pnt4[0] = 12.5; pnt4[1] = 10.5; std::vector points{ pnt1,pnt2,pnt3,pnt4 }; auto figure = GeneratePlanarPolygon(m_Geometry, points); mitk::PlanarFigureMaskGenerator::Pointer planFigMaskGen = mitk::PlanarFigureMaskGenerator::New(); planFigMaskGen->SetInputImage(m_TestImage); planFigMaskGen->SetPlanarFigure(figure.GetPointer()); mitk::ImageStatisticsContainer::Pointer statisticsContainer; CPPUNIT_ASSERT_NO_THROW(statisticsContainer = ComputeStatistics(m_TestImage, planFigMaskGen.GetPointer())); auto statisticsObjectTimestep0 = statisticsContainer->GetStatistics(1, 0); this->VerifyStatistics(statisticsObjectTimestep0, 191.5, 89.802561210691536, 128.63499999046327); } 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:-----------------------------------------------------------------------------------"; std::string filename = this->GetTestDataFilePath("ImageStatisticsTestData/testimage.nrrd"); m_TestImage = mitk::IOUtil::Load(filename); m_Geometry = m_TestImage->GetSlicedGeometry()->GetPlaneGeometry(0); mitk::Point2D pnt1; pnt1[0] = 11.5; pnt1[1] = 10.5; mitk::Point2D pnt2; pnt2[0] = 11.5; pnt2[1] = 13.5; mitk::Point2D pnt3; pnt3[0] = 12.5; pnt3[1] = 13.5; mitk::Point2D pnt4; pnt4[0] = 12.5; pnt4[1] = 10.5; std::vector points{ pnt1,pnt2,pnt3,pnt4 }; auto figure = GeneratePlanarPolygon(m_Geometry, points); mitk::PlanarFigureMaskGenerator::Pointer planFigMaskGen = mitk::PlanarFigureMaskGenerator::New(); planFigMaskGen->SetInputImage(m_TestImage); planFigMaskGen->SetPlanarFigure(figure.GetPointer()); mitk::ImageStatisticsContainer::Pointer statisticsContainer; CPPUNIT_ASSERT_NO_THROW(statisticsContainer = ComputeStatistics(m_TestImage, planFigMaskGen.GetPointer())); auto statisticsObjectTimestep0 = statisticsContainer->GetStatistics(1, 0); this->VerifyStatistics(statisticsObjectTimestep0, 127.666666666666667, 127.50032679696680, 128.7750015258789); } 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:-----------------------------------------------------------------------------------"; std::string filename = this->GetTestDataFilePath("ImageStatisticsTestData/testimage.nrrd"); m_TestImage = mitk::IOUtil::Load(filename); m_Geometry = m_TestImage->GetSlicedGeometry()->GetPlaneGeometry(0); mitk::Point2D pnt1; pnt1[0] = 0.5; pnt1[1] = 0.5; mitk::Point2D pnt2; pnt2[0] = 3.5; pnt2[1] = 3.5; mitk::Point2D pnt3; pnt3[0] = 8.4999; pnt3[1] = 3.5; mitk::Point2D pnt4; pnt4[0] = 5.4999; pnt4[1] = 0.5; std::vector points{ pnt1,pnt2,pnt3,pnt4 }; auto figure = GeneratePlanarPolygon(m_Geometry, points); mitk::PlanarFigureMaskGenerator::Pointer planFigMaskGen = mitk::PlanarFigureMaskGenerator::New(); planFigMaskGen->SetInputImage(m_TestImage); planFigMaskGen->SetPlanarFigure(figure.GetPointer()); mitk::ImageStatisticsContainer::Pointer statisticsContainer; CPPUNIT_ASSERT_NO_THROW(statisticsContainer = ComputeStatistics(m_TestImage, planFigMaskGen.GetPointer())); auto statisticsObjectTimestep0 = statisticsContainer->GetStatistics(1, 0); this->VerifyStatistics(statisticsObjectTimestep0, 204.0, 105.58003057938019, 253.724998474121083); } void mitkImageStatisticsCalculatorTestSuite::TestCase12() { /***************************** * half pixel (white) + whole pixel (white) + half pixel (black) * -> mean of 212.66 expected ******************************/ MITK_INFO << std::endl << "Test case 12:-----------------------------------------------------------------------------------"; std::string filename = this->GetTestDataFilePath("ImageStatisticsTestData/testimage.nrrd"); m_TestImage = mitk::IOUtil::Load(filename); m_Geometry = m_TestImage->GetSlicedGeometry()->GetPlaneGeometry(0); mitk::Point2D pnt1; pnt1[0] = 9.5; pnt1[1] = 0.5; mitk::Point2D pnt2; pnt2[0] = 9.5; pnt2[1] = 2.5; mitk::Point2D pnt3; pnt3[0] = 11.5; pnt3[1] = 2.5; std::vector points{ pnt1,pnt2,pnt3 }; auto figure = GeneratePlanarPolygon(m_Geometry, points); mitk::PlanarFigureMaskGenerator::Pointer planFigMaskGen = mitk::PlanarFigureMaskGenerator::New(); planFigMaskGen->SetInputImage(m_TestImage); planFigMaskGen->SetPlanarFigure(figure.GetPointer()); mitk::ImageStatisticsContainer::Pointer statisticsContainer; CPPUNIT_ASSERT_NO_THROW(statisticsContainer = ComputeStatistics(m_TestImage, planFigMaskGen.GetPointer())); auto statisticsObjectTimestep0 = statisticsContainer->GetStatistics(1, 0); this->VerifyStatistics(statisticsObjectTimestep0, 212.666666666666667, 73.323484187082443, 254.36499786376954); } // T26098 histogram statistics need to be tested (median, uniformity, UPP, entropy) void mitkImageStatisticsCalculatorTestSuite::TestPic3DCroppedNoMask() { MITK_INFO << std::endl << "Test Pic3D cropped without mask:-----------------------------------------------------------------------------------"; std::string Pic3DCroppedFile = this->GetTestDataFilePath("ImageStatisticsTestData/Pic3D_cropped.nrrd"); m_Pic3DCroppedImage = mitk::IOUtil::Load(Pic3DCroppedFile); CPPUNIT_ASSERT_MESSAGE("Failed loading Pic3D_cropped", m_Pic3DCroppedImage.IsNotNull()); //calculated ground truth via script mitk::ImageStatisticsContainer::VoxelCountType expected_N = 27; mitk::ImageStatisticsContainer::RealType expected_mean = -564.1481481481481481; mitk::ImageStatisticsContainer::RealType expected_MPP = 113.66666666666667; //mitk::ImageStatisticsContainer::RealType expected_median = -825; mitk::ImageStatisticsContainer::RealType expected_skewness = 0.7120461106763573; mitk::ImageStatisticsContainer::RealType expected_kurtosis = 1.8794464383714844; mitk::ImageStatisticsContainer::RealType expected_variance = 145946.82336182334; mitk::ImageStatisticsContainer::RealType expected_standarddev = 382.02987234223366; mitk::ImageStatisticsContainer::RealType expected_min = -927; mitk::ImageStatisticsContainer::RealType expected_max = 147; mitk::ImageStatisticsContainer::RealType expected_RMS = 681.32955052662169; mitk::ImageStatisticsContainer::IndexType expected_minIndex; expected_minIndex.set_size(3); expected_minIndex[0] = 2; expected_minIndex[1] = 1; expected_minIndex[2] = 1; mitk::ImageStatisticsContainer::IndexType expected_maxIndex; expected_maxIndex.set_size(3); expected_maxIndex[0] = 0; expected_maxIndex[1] = 1; expected_maxIndex[2] = 2; mitk::ImageStatisticsContainer::Pointer statisticsContainer; CPPUNIT_ASSERT_NO_THROW(statisticsContainer = ComputeStatistics(m_Pic3DCroppedImage)); auto statisticsObject = statisticsContainer->GetStatistics(mitk::ImageStatisticsContainer::NO_MASK_LABEL_VALUE, 0); VerifyStatistics(statisticsObject, expected_N, expected_mean, expected_MPP, expected_skewness, expected_kurtosis, expected_variance, expected_standarddev, expected_min, expected_max, expected_RMS, expected_minIndex, expected_maxIndex); } // T26098 histogram statistics need to be tested (median, uniformity, UPP, entropy) void mitkImageStatisticsCalculatorTestSuite::TestPic3DCroppedBinMask() { MITK_INFO << std::endl << "Test Pic3D cropped binary mask:-----------------------------------------------------------------------------------"; std::string Pic3DCroppedFile = this->GetTestDataFilePath("ImageStatisticsTestData/Pic3D_cropped.nrrd"); m_Pic3DCroppedImage = mitk::IOUtil::Load(Pic3DCroppedFile); CPPUNIT_ASSERT_MESSAGE("Failed loading Pic3D_cropped", m_Pic3DCroppedImage.IsNotNull()); std::string Pic3DCroppedBinMaskFile = this->GetTestDataFilePath("ImageStatisticsTestData/Pic3D_croppedBinMask.nrrd"); m_Pic3DCroppedBinMask = mitk::IOUtil::Load(Pic3DCroppedBinMaskFile); CPPUNIT_ASSERT_MESSAGE("Failed loading Pic3D binary mask", m_Pic3DCroppedBinMask.IsNotNull()); //calculated ground truth via script mitk::ImageStatisticsContainer::RealType expected_kurtosis = 1.0765697398089618; mitk::ImageStatisticsContainer::RealType expected_MPP = -nan(""); mitk::ImageStatisticsContainer::RealType expected_max = -22; mitk::ImageStatisticsContainer::RealType expected_mean = -464; mitk::ImageStatisticsContainer::RealType expected_min = -846; mitk::ImageStatisticsContainer::VoxelCountType expected_N = 4; mitk::ImageStatisticsContainer::RealType expected_RMS = 633.20191618998331; mitk::ImageStatisticsContainer::RealType expected_skewness = 0.0544059290851858; mitk::ImageStatisticsContainer::RealType expected_standarddev = 430.86966320067910; mitk::ImageStatisticsContainer::RealType expected_variance = 185648.66666666663; mitk::ImageStatisticsContainer::IndexType expected_minIndex; expected_minIndex.set_size(3); expected_minIndex[0] = 1; expected_minIndex[1] = 0; expected_minIndex[2] = 0; mitk::ImageStatisticsContainer::IndexType expected_maxIndex; expected_maxIndex.set_size(3); expected_maxIndex[0] = 0; expected_maxIndex[1] = 0; expected_maxIndex[2] = 1; mitk::ImageMaskGenerator::Pointer imgMaskGen = mitk::ImageMaskGenerator::New(); imgMaskGen->SetImageMask(m_Pic3DCroppedBinMask); imgMaskGen->SetInputImage(m_Pic3DCroppedImage); - imgMaskGen->SetTimeStep(0); + imgMaskGen->SetTimePoint(m_Pic3DCroppedImage->GetTimeGeometry()->TimeStepToTimePoint(0)); mitk::ImageStatisticsContainer::Pointer statisticsContainer; CPPUNIT_ASSERT_NO_THROW(statisticsContainer = ComputeStatistics(m_Pic3DCroppedImage, imgMaskGen.GetPointer(), nullptr)); auto statisticsObjectTimestep0 = statisticsContainer->GetStatistics(1, 0); VerifyStatistics(statisticsObjectTimestep0, expected_N, expected_mean, expected_MPP, expected_skewness, expected_kurtosis, expected_variance, expected_standarddev, expected_min, expected_max, expected_RMS, expected_minIndex, expected_maxIndex); } // T26098 histogram statistics need to be tested (median, uniformity, UPP, entropy) void mitkImageStatisticsCalculatorTestSuite::TestPic3DCroppedMultilabelMask() { MITK_INFO << std::endl << "Test Pic3D cropped multi-label mask:-----------------------------------------------------------------------------------"; std::string Pic3DCroppedFile = this->GetTestDataFilePath("ImageStatisticsTestData/Pic3D_cropped.nrrd"); m_Pic3DCroppedImage = mitk::IOUtil::Load(Pic3DCroppedFile); CPPUNIT_ASSERT_MESSAGE("Failed loading Pic3D_cropped", m_Pic3DCroppedImage.IsNotNull()); std::string Pic3DCroppedMultilabelMaskFile = this->GetTestDataFilePath("ImageStatisticsTestData/Pic3D_croppedMultilabelMask.nrrd"); m_Pic3DCroppedMultilabelMask = mitk::IOUtil::Load(Pic3DCroppedMultilabelMaskFile); CPPUNIT_ASSERT_MESSAGE("Failed loading Pic3D multi-label mask", m_Pic3DCroppedMultilabelMask.IsNotNull()); //calculated ground truth via script mitk::ImageStatisticsContainer::RealType expected_kurtosis = 1.5; mitk::ImageStatisticsContainer::RealType expected_MPP = -nan(""); mitk::ImageStatisticsContainer::RealType expected_max = -22; mitk::ImageStatisticsContainer::RealType expected_mean = -586.33333333333333; mitk::ImageStatisticsContainer::RealType expected_min = -916; mitk::ImageStatisticsContainer::VoxelCountType expected_N = 3; mitk::ImageStatisticsContainer::RealType expected_RMS = 764.78566351044469; mitk::ImageStatisticsContainer::RealType expected_skewness = 0.6774469597523700; mitk::ImageStatisticsContainer::RealType expected_standarddev = 491.02987010296363; mitk::ImageStatisticsContainer::RealType expected_variance = 241110.33333333334; mitk::ImageStatisticsContainer::IndexType expected_minIndex; expected_minIndex.set_size(3); expected_minIndex[0] = 2; expected_minIndex[1] = 0; expected_minIndex[2] = 1; mitk::ImageStatisticsContainer::IndexType expected_maxIndex; expected_maxIndex.set_size(3); expected_maxIndex[0] = 0; expected_maxIndex[1] = 0; expected_maxIndex[2] = 1; mitk::ImageMaskGenerator::Pointer imgMaskGen = mitk::ImageMaskGenerator::New(); imgMaskGen->SetImageMask(m_Pic3DCroppedMultilabelMask); imgMaskGen->SetInputImage(m_Pic3DCroppedImage); - imgMaskGen->SetTimeStep(0); + imgMaskGen->SetTimePoint(m_Pic3DCroppedImage->GetTimeGeometry()->TimeStepToTimePoint(0)); mitk::ImageStatisticsContainer::Pointer statisticsContainer; CPPUNIT_ASSERT_NO_THROW(statisticsContainer = ComputeStatistics(m_Pic3DCroppedImage, imgMaskGen.GetPointer(), nullptr)); auto statisticsObjectTimestep0 = statisticsContainer->GetStatistics(2, 0); VerifyStatistics(statisticsObjectTimestep0, expected_N, expected_mean, expected_MPP, expected_skewness, expected_kurtosis, expected_variance, expected_standarddev, expected_min, expected_max, expected_RMS, expected_minIndex, expected_maxIndex); } // T26098 histogram statistics need to be tested (median, uniformity, UPP, entropy) void mitkImageStatisticsCalculatorTestSuite::TestPic3DCroppedPlanarFigure() { MITK_INFO << std::endl << "Test Pic3D cropped planar figure:-----------------------------------------------------------------------------------"; std::string Pic3DCroppedFile = this->GetTestDataFilePath("ImageStatisticsTestData/Pic3D_cropped.nrrd"); m_Pic3DCroppedImage = mitk::IOUtil::Load(Pic3DCroppedFile); CPPUNIT_ASSERT_MESSAGE("Failed loading Pic3D_cropped", m_Pic3DCroppedImage.IsNotNull()); std::string Pic3DCroppedPlanarFigureFile = this->GetTestDataFilePath("ImageStatisticsTestData/Pic3D_croppedPF.pf"); m_Pic3DCroppedPlanarFigure = mitk::IOUtil::Load(Pic3DCroppedPlanarFigureFile); CPPUNIT_ASSERT_MESSAGE("Failed loading Pic3D planar figure", m_Pic3DCroppedPlanarFigure.IsNotNull()); //calculated ground truth via script mitk::ImageStatisticsContainer::RealType expected_kurtosis = 1; mitk::ImageStatisticsContainer::RealType expected_MPP = -nan(""); mitk::ImageStatisticsContainer::RealType expected_max = -67; mitk::ImageStatisticsContainer::RealType expected_mean = -446; mitk::ImageStatisticsContainer::RealType expected_min = -825; mitk::ImageStatisticsContainer::VoxelCountType expected_N = 2; mitk::ImageStatisticsContainer::RealType expected_RMS = 697.27899724572228; mitk::ImageStatisticsContainer::RealType expected_skewness = 0; mitk::ImageStatisticsContainer::RealType expected_standarddev = 535.98694013940303; mitk::ImageStatisticsContainer::RealType expected_variance = 287282.0; mitk::ImageStatisticsContainer::IndexType expected_minIndex; expected_minIndex.set_size(3); expected_minIndex[0] = 1; expected_minIndex[1] = 1; expected_minIndex[2] = 1; mitk::ImageStatisticsContainer::IndexType expected_maxIndex; expected_maxIndex.set_size(3); expected_maxIndex[0] = 0; expected_maxIndex[1] = 1; expected_maxIndex[2] = 1; mitk::PlanarFigureMaskGenerator::Pointer pfMaskGen = mitk::PlanarFigureMaskGenerator::New(); pfMaskGen->SetInputImage(m_Pic3DCroppedImage); pfMaskGen->SetPlanarFigure(m_Pic3DCroppedPlanarFigure); mitk::ImageStatisticsContainer::Pointer statisticsContainer; CPPUNIT_ASSERT_NO_THROW(statisticsContainer = ComputeStatistics(m_Pic3DCroppedImage, pfMaskGen.GetPointer())); auto statisticsObjectTimestep0 = statisticsContainer->GetStatistics(1, 0); VerifyStatistics(statisticsObjectTimestep0, expected_N, expected_mean, expected_MPP, expected_skewness, expected_kurtosis, expected_variance, expected_standarddev, expected_min, expected_max, expected_RMS, expected_minIndex, expected_maxIndex); } // T26098 histogram statistics need to be tested (median, uniformity, UPP, entropy) void mitkImageStatisticsCalculatorTestSuite::TestUS4DCroppedNoMaskTimeStep1() { MITK_INFO << std::endl << "Test US4D cropped without mask time step 1:-----------------------------------------------------------------------------------"; std::string US4DCroppedFile = this->GetTestDataFilePath("ImageStatisticsTestData/US4D_cropped.nrrd"); m_US4DCroppedImage = mitk::IOUtil::Load(US4DCroppedFile); CPPUNIT_ASSERT_MESSAGE("Failed loading US4D_cropped", m_US4DCroppedImage.IsNotNull()); //calculated ground truth via script mitk::ImageStatisticsContainer::RealType expected_kurtosis = 1.5398359155908228; mitk::ImageStatisticsContainer::RealType expected_MPP = 157.74074074074073; mitk::ImageStatisticsContainer::RealType expected_max = 199; mitk::ImageStatisticsContainer::RealType expected_mean = 157.74074074074073; mitk::ImageStatisticsContainer::RealType expected_min = 101; mitk::ImageStatisticsContainer::VoxelCountType expected_N = 27; mitk::ImageStatisticsContainer::RealType expected_RMS = 161.11544579426010; mitk::ImageStatisticsContainer::RealType expected_skewness = 0.0347280313508018; mitk::ImageStatisticsContainer::RealType expected_standarddev = 32.803133753432512; mitk::ImageStatisticsContainer::RealType expected_variance = 1076.0455840455834; mitk::ImageStatisticsContainer::IndexType expected_minIndex; expected_minIndex.set_size(3); expected_minIndex[0] = 0; expected_minIndex[1] = 2; expected_minIndex[2] = 0; mitk::ImageStatisticsContainer::IndexType expected_maxIndex; expected_maxIndex.set_size(3); expected_maxIndex[0] = 0; expected_maxIndex[1] = 0; expected_maxIndex[2] = 1; mitk::ImageStatisticsContainer::Pointer statisticsContainer=mitk::ImageStatisticsContainer::New(); CPPUNIT_ASSERT_NO_THROW(statisticsContainer = ComputeStatistics(m_US4DCroppedImage)); auto statisticsObjectTimestep1 = statisticsContainer->GetStatistics(mitk::ImageStatisticsContainer::NO_MASK_LABEL_VALUE, 1); VerifyStatistics(statisticsObjectTimestep1, expected_N, expected_mean, expected_MPP, expected_skewness, expected_kurtosis, expected_variance, expected_standarddev, expected_min, expected_max, expected_RMS, expected_minIndex, expected_maxIndex); } // T26098 histogram statistics need to be tested (median, uniformity, UPP, entropy) void mitkImageStatisticsCalculatorTestSuite::TestUS4DCroppedBinMaskTimeStep1() { MITK_INFO << std::endl << "Test US4D cropped with binary mask time step 1:-----------------------------------------------------------------------------------"; std::string US4DCroppedFile = this->GetTestDataFilePath("ImageStatisticsTestData/US4D_cropped.nrrd"); m_US4DCroppedImage = mitk::IOUtil::Load(US4DCroppedFile); CPPUNIT_ASSERT_MESSAGE("Failed loading US4D_cropped", m_US4DCroppedImage.IsNotNull()); std::string US4DCroppedBinMaskFile = this->GetTestDataFilePath("ImageStatisticsTestData/US4D_croppedBinMask.nrrd"); m_US4DCroppedBinMask = mitk::IOUtil::Load(US4DCroppedBinMaskFile); CPPUNIT_ASSERT_MESSAGE("Failed loading US4D binary mask", m_US4DCroppedBinMask.IsNotNull()); //calculated ground truth via script mitk::ImageStatisticsContainer::RealType expected_kurtosis = 1.5863739712889191; mitk::ImageStatisticsContainer::RealType expected_MPP = 166.75; mitk::ImageStatisticsContainer::RealType expected_max = 199; mitk::ImageStatisticsContainer::RealType expected_mean = 166.75; mitk::ImageStatisticsContainer::RealType expected_min = 120; mitk::ImageStatisticsContainer::VoxelCountType expected_N = 4; mitk::ImageStatisticsContainer::RealType expected_RMS = 170.68043971117487; mitk::ImageStatisticsContainer::RealType expected_skewness = -0.4285540263894276; mitk::ImageStatisticsContainer::RealType expected_standarddev = 36.417715469260287; mitk::ImageStatisticsContainer::RealType expected_variance = 1326.25; mitk::ImageStatisticsContainer::IndexType expected_minIndex; expected_minIndex.set_size(3); expected_minIndex[0] = 0; expected_minIndex[1] = 0; expected_minIndex[2] = 2; mitk::ImageStatisticsContainer::IndexType expected_maxIndex; expected_maxIndex.set_size(3); expected_maxIndex[0] = 1; expected_maxIndex[1] = 1; expected_maxIndex[2] = 1; mitk::ImageMaskGenerator::Pointer imgMask1 = mitk::ImageMaskGenerator::New(); imgMask1->SetInputImage(m_US4DCroppedImage); imgMask1->SetImageMask(m_US4DCroppedBinMask); mitk::ImageStatisticsContainer::Pointer statisticsContainer=mitk::ImageStatisticsContainer::New(); CPPUNIT_ASSERT_NO_THROW(statisticsContainer = ComputeStatistics(m_US4DCroppedImage, imgMask1.GetPointer(), nullptr)); auto statisticsObjectTimestep1 = statisticsContainer->GetStatistics(1, 1); VerifyStatistics(statisticsObjectTimestep1, expected_N, expected_mean, expected_MPP, expected_skewness, expected_kurtosis, expected_variance, expected_standarddev, expected_min, expected_max, expected_RMS, expected_minIndex, expected_maxIndex); } // T26098 histogram statistics need to be tested (median, uniformity, UPP, entropy) void mitkImageStatisticsCalculatorTestSuite::TestUS4DCroppedMultilabelMaskTimeStep1() { MITK_INFO << std::endl << "Test US4D cropped with multi-label mask time step 1:-----------------------------------------------------------------------------------"; std::string US4DCroppedFile = this->GetTestDataFilePath("ImageStatisticsTestData/US4D_cropped.nrrd"); m_US4DCroppedImage = mitk::IOUtil::Load(US4DCroppedFile); CPPUNIT_ASSERT_MESSAGE("Failed loading US4D_cropped", m_US4DCroppedImage.IsNotNull()); std::string US4DCroppedMultilabelMaskFile = this->GetTestDataFilePath("ImageStatisticsTestData/US4D_croppedMultilabelMask.nrrd"); m_US4DCroppedMultilabelMask = mitk::IOUtil::Load(US4DCroppedMultilabelMaskFile); CPPUNIT_ASSERT_MESSAGE("Failed loading US4D multi-label mask", m_US4DCroppedMultilabelMask.IsNotNull()); //calculated ground truth via script mitk::ImageStatisticsContainer::RealType expected_kurtosis = 1.0432484564918287; mitk::ImageStatisticsContainer::RealType expected_MPP = 159.75; mitk::ImageStatisticsContainer::RealType expected_max = 199; mitk::ImageStatisticsContainer::RealType expected_mean = 159.75; mitk::ImageStatisticsContainer::RealType expected_min = 120; mitk::ImageStatisticsContainer::VoxelCountType expected_N = 4; mitk::ImageStatisticsContainer::RealType expected_RMS = 165.05447333128134; mitk::ImageStatisticsContainer::RealType expected_skewness = -0.004329226115093; mitk::ImageStatisticsContainer::RealType expected_standarddev = 41.508031351374242; mitk::ImageStatisticsContainer::RealType expected_variance = 1722.9166666666670; mitk::ImageStatisticsContainer::IndexType expected_minIndex; expected_minIndex.set_size(3); expected_minIndex[0] = 0; expected_minIndex[1] = 0; expected_minIndex[2] = 2; mitk::ImageStatisticsContainer::IndexType expected_maxIndex; expected_maxIndex.set_size(3); expected_maxIndex[0] = 0; expected_maxIndex[1] = 0; expected_maxIndex[2] = 1; mitk::ImageMaskGenerator::Pointer imgMask1 = mitk::ImageMaskGenerator::New(); imgMask1->SetInputImage(m_US4DCroppedImage); imgMask1->SetImageMask(m_US4DCroppedMultilabelMask); mitk::ImageStatisticsContainer::Pointer statisticsContainer; CPPUNIT_ASSERT_NO_THROW(statisticsContainer = ComputeStatistics(m_US4DCroppedImage, imgMask1.GetPointer(), nullptr)); auto statisticsObjectTimestep1 = statisticsContainer->GetStatistics(1, 1); VerifyStatistics(statisticsObjectTimestep1, expected_N, expected_mean, expected_MPP, expected_skewness, expected_kurtosis, expected_variance, expected_standarddev, expected_min, expected_max, expected_RMS, expected_minIndex, expected_maxIndex); } // T26098 histogram statistics need to be tested (median, uniformity, UPP, entropy) void mitkImageStatisticsCalculatorTestSuite::TestUS4DCroppedPlanarFigureTimeStep1() { MITK_INFO << std::endl << "Test US4D cropped planar figure time step 1:-----------------------------------------------------------------------------------"; std::string US4DCroppedFile = this->GetTestDataFilePath("ImageStatisticsTestData/US4D_cropped.nrrd"); m_US4DCroppedImage = mitk::IOUtil::Load(US4DCroppedFile); CPPUNIT_ASSERT_MESSAGE("Failed loading US4D_cropped", m_US4DCroppedImage.IsNotNull()); std::string US4DCroppedPlanarFigureFile = this->GetTestDataFilePath("ImageStatisticsTestData/US4D_croppedPF.pf"); m_US4DCroppedPlanarFigure = mitk::IOUtil::Load(US4DCroppedPlanarFigureFile); CPPUNIT_ASSERT_MESSAGE("Failed loading US4D planar figure", m_US4DCroppedPlanarFigure.IsNotNull()); //calculated ground truth via script mitk::ImageStatisticsContainer::RealType expected_kurtosis = 1; mitk::ImageStatisticsContainer::RealType expected_MPP = 172.5; mitk::ImageStatisticsContainer::RealType expected_max = 197; mitk::ImageStatisticsContainer::RealType expected_mean = 172.5; mitk::ImageStatisticsContainer::RealType expected_min = 148; mitk::ImageStatisticsContainer::VoxelCountType expected_N = 2; mitk::ImageStatisticsContainer::RealType expected_RMS = 175.94530400098776; mitk::ImageStatisticsContainer::RealType expected_skewness = 0; mitk::ImageStatisticsContainer::RealType expected_standarddev = 34.648232278140831; mitk::ImageStatisticsContainer::RealType expected_variance = 1200.5000000000002; mitk::ImageStatisticsContainer::IndexType expected_minIndex; expected_minIndex.set_size(3); expected_minIndex[0] = 2; expected_minIndex[1] = 2; expected_minIndex[2] = 2; mitk::ImageStatisticsContainer::IndexType expected_maxIndex; expected_maxIndex.set_size(3); expected_maxIndex[0] = 2; expected_maxIndex[1] = 2; expected_maxIndex[2] = 1; mitk::PlanarFigureMaskGenerator::Pointer pfMaskGen = mitk::PlanarFigureMaskGenerator::New(); pfMaskGen->SetInputImage(m_US4DCroppedImage); pfMaskGen->SetPlanarFigure(m_US4DCroppedPlanarFigure); mitk::ImageStatisticsContainer::Pointer statisticsContainer; CPPUNIT_ASSERT_NO_THROW(statisticsContainer = ComputeStatistics(m_US4DCroppedImage, pfMaskGen.GetPointer())); auto statisticsObjectTimestep1 = statisticsContainer->GetStatistics(1, 1); VerifyStatistics(statisticsObjectTimestep1, expected_N, expected_mean, expected_MPP, expected_skewness, expected_kurtosis, expected_variance, expected_standarddev, expected_min, expected_max, expected_RMS, expected_minIndex, expected_maxIndex); } void mitkImageStatisticsCalculatorTestSuite::TestUS4DCroppedAllTimesteps() { MITK_INFO << std::endl << "Test US4D cropped all time steps:-----------------------------------------------------------------------------------"; std::string US4DCroppedFile = this->GetTestDataFilePath("ImageStatisticsTestData/US4D_cropped.nrrd"); m_US4DCroppedImage = mitk::IOUtil::Load(US4DCroppedFile); CPPUNIT_ASSERT_MESSAGE("Failed loading US4D_cropped", m_US4DCroppedImage.IsNotNull()); mitk::ImageStatisticsContainer::Pointer statisticsContainer=mitk::ImageStatisticsContainer::New(); CPPUNIT_ASSERT_NO_THROW(statisticsContainer = ComputeStatistics(m_US4DCroppedImage)); for (int i = 0; i < 4; i++) { CPPUNIT_ASSERT_MESSAGE("Error computing statistics for multiple time steps", statisticsContainer->StatisticsExist(mitk::ImageStatisticsContainer::NO_MASK_LABEL_VALUE, i)); } } void mitkImageStatisticsCalculatorTestSuite::TestUS4DCropped3DMask() { MITK_INFO << std::endl << "Test US4D cropped with 3D binary Mask:-----------------------------------------------------------------------------------"; std::string US4DCroppedFile = this->GetTestDataFilePath("ImageStatisticsTestData/US4D_cropped.nrrd"); m_US4DCroppedImage = mitk::IOUtil::Load(US4DCroppedFile); CPPUNIT_ASSERT_MESSAGE("Failed loading US4D_cropped", m_US4DCroppedImage.IsNotNull()); std::string US4DCropped3DBinMaskFile = this->GetTestDataFilePath("ImageStatisticsTestData/US4D_cropped3DBinMask.nrrd"); m_US4DCropped3DBinMask = mitk::IOUtil::Load(US4DCropped3DBinMaskFile); CPPUNIT_ASSERT_MESSAGE("Failed loading Pic3D binary mask", m_US4DCropped3DBinMask.IsNotNull()); //calculated ground truth via script mitk::ImageStatisticsContainer::RealType expected_kurtosis = 1; mitk::ImageStatisticsContainer::RealType expected_MPP = 198; mitk::ImageStatisticsContainer::RealType expected_max = 199; mitk::ImageStatisticsContainer::RealType expected_mean = 198; mitk::ImageStatisticsContainer::RealType expected_min = 197; mitk::ImageStatisticsContainer::VoxelCountType expected_N = 2; mitk::ImageStatisticsContainer::RealType expected_RMS = 198.00505044063902; mitk::ImageStatisticsContainer::RealType expected_skewness = 0; mitk::ImageStatisticsContainer::RealType expected_standarddev = 1.4142135623730951; mitk::ImageStatisticsContainer::RealType expected_variance = 2; mitk::ImageStatisticsContainer::IndexType expected_minIndex; expected_minIndex.set_size(3); expected_minIndex[0] = 1; expected_minIndex[1] = 2; expected_minIndex[2] = 1; mitk::ImageStatisticsContainer::IndexType expected_maxIndex; expected_maxIndex.set_size(3); expected_maxIndex[0] = 1; expected_maxIndex[1] = 1; expected_maxIndex[2] = 1; mitk::ImageMaskGenerator::Pointer imgMask1 = mitk::ImageMaskGenerator::New(); imgMask1->SetInputImage(m_US4DCroppedImage); imgMask1->SetImageMask(m_US4DCropped3DBinMask); mitk::ImageStatisticsContainer::Pointer statisticsContainer = mitk::ImageStatisticsContainer::New(); CPPUNIT_ASSERT_NO_THROW(statisticsContainer = ComputeStatistics(m_US4DCroppedImage, imgMask1.GetPointer(), nullptr)); auto statisticsObjectTimestep1 = statisticsContainer->GetStatistics(1, 1); VerifyStatistics(statisticsObjectTimestep1, expected_N, expected_mean, expected_MPP, expected_skewness, expected_kurtosis, expected_variance, expected_standarddev, expected_min, expected_max, expected_RMS, expected_minIndex, expected_maxIndex); } mitk::PlanarPolygon::Pointer mitkImageStatisticsCalculatorTestSuite::GeneratePlanarPolygon(mitk::PlaneGeometry::Pointer geometry, std::vector points) { mitk::PlanarPolygon::Pointer figure = mitk::PlanarPolygon::New(); figure->SetPlaneGeometry(geometry); figure->PlaceFigure(points[0]); for (unsigned int i = 1; i < points.size(); i++) { figure->SetControlPoint(i, points[i], true); } return figure; } void mitkImageStatisticsCalculatorTestSuite::VerifyStatistics(mitk::ImageStatisticsContainer::ImageStatisticsObject stats, mitk::ImageStatisticsContainer::RealType testMean, mitk::ImageStatisticsContainer::RealType testSD, mitk::ImageStatisticsContainer::RealType testMedian) { mitk::ImageStatisticsContainer::RealType meanObject = 0; mitk::ImageStatisticsContainer::RealType standardDeviationObject = 0; mitk::ImageStatisticsContainer::RealType medianObject = 0; CPPUNIT_ASSERT_NO_THROW(meanObject = stats.GetValueConverted(mitk::ImageStatisticsConstants::MEAN())); CPPUNIT_ASSERT_NO_THROW(standardDeviationObject = stats.GetValueConverted(mitk::ImageStatisticsConstants::STANDARDDEVIATION())); CPPUNIT_ASSERT_NO_THROW(medianObject = stats.GetValueConverted(mitk::ImageStatisticsConstants::MEDIAN())); CPPUNIT_ASSERT_MESSAGE("Calculated mean gray value is not equal to the desired value.", std::abs(meanObject - testMean) < mitk::eps); CPPUNIT_ASSERT_MESSAGE("Calculated gray value sd is not equal to the desired value.", std::abs(standardDeviationObject - testSD) < mitk::eps); CPPUNIT_ASSERT_MESSAGE("Calculated median gray value is not equal to the desired value.", std::abs(medianObject - testMedian) < mitk::eps); } // T26098 histogram statistics need to be tested (median, uniformity, UPP, entropy) void mitkImageStatisticsCalculatorTestSuite::VerifyStatistics(mitk::ImageStatisticsContainer::ImageStatisticsObject stats, mitk::ImageStatisticsContainer::VoxelCountType N, mitk::ImageStatisticsContainer::RealType mean, mitk::ImageStatisticsContainer::RealType MPP, mitk::ImageStatisticsContainer::RealType skewness, mitk::ImageStatisticsContainer::RealType kurtosis, mitk::ImageStatisticsContainer::RealType variance, mitk::ImageStatisticsContainer::RealType stdev, mitk::ImageStatisticsContainer::RealType min, mitk::ImageStatisticsContainer::RealType max, mitk::ImageStatisticsContainer::RealType RMS, mitk::ImageStatisticsContainer::IndexType minIndex, mitk::ImageStatisticsContainer::IndexType maxIndex) { mitk::ImageStatisticsContainer::VoxelCountType numberOfVoxelsObject; mitk::ImageStatisticsContainer::RealType meanObject = 0; mitk::ImageStatisticsContainer::RealType mppObject = 0; mitk::ImageStatisticsContainer::RealType skewnessObject = 0; mitk::ImageStatisticsContainer::RealType kurtosisObject = 0; mitk::ImageStatisticsContainer::RealType varianceObject = 0; mitk::ImageStatisticsContainer::RealType standardDeviationObject = 0; mitk::ImageStatisticsContainer::RealType minObject = 0; mitk::ImageStatisticsContainer::RealType maxObject = 0; mitk::ImageStatisticsContainer::RealType rmsObject = 0; mitk::ImageStatisticsContainer::IndexType minIndexObject(3,0); mitk::ImageStatisticsContainer::IndexType maxIndexObject(3,0); CPPUNIT_ASSERT_NO_THROW(numberOfVoxelsObject = stats.GetValueConverted(mitk::ImageStatisticsConstants::NUMBEROFVOXELS())); CPPUNIT_ASSERT_NO_THROW(meanObject = stats.GetValueConverted(mitk::ImageStatisticsConstants::MEAN())); CPPUNIT_ASSERT_NO_THROW(mppObject = stats.GetValueConverted(mitk::ImageStatisticsConstants::MPP())); CPPUNIT_ASSERT_NO_THROW(skewnessObject = stats.GetValueConverted(mitk::ImageStatisticsConstants::SKEWNESS())); CPPUNIT_ASSERT_NO_THROW(kurtosisObject = stats.GetValueConverted(mitk::ImageStatisticsConstants::KURTOSIS())); CPPUNIT_ASSERT_NO_THROW(varianceObject = stats.GetValueConverted(mitk::ImageStatisticsConstants::VARIANCE())); CPPUNIT_ASSERT_NO_THROW(standardDeviationObject = stats.GetValueConverted(mitk::ImageStatisticsConstants::STANDARDDEVIATION())); CPPUNIT_ASSERT_NO_THROW(minObject = stats.GetValueConverted(mitk::ImageStatisticsConstants::MINIMUM())); CPPUNIT_ASSERT_NO_THROW(maxObject = stats.GetValueConverted(mitk::ImageStatisticsConstants::MAXIMUM())); CPPUNIT_ASSERT_NO_THROW(rmsObject = stats.GetValueConverted(mitk::ImageStatisticsConstants::RMS())); CPPUNIT_ASSERT_NO_THROW(minIndexObject = stats.GetValueConverted(mitk::ImageStatisticsConstants::MINIMUMPOSITION())); CPPUNIT_ASSERT_NO_THROW(maxIndexObject = stats.GetValueConverted(mitk::ImageStatisticsConstants::MAXIMUMPOSITION())); CPPUNIT_ASSERT_MESSAGE("Calculated value does not fit expected value", numberOfVoxelsObject - N == 0); CPPUNIT_ASSERT_MESSAGE("Calculated value does not fit expected value", std::abs(meanObject - mean) < mitk::eps); // in three test cases MPP is None because the ROI has no positive pixels if (!std::isnan(mppObject)) { CPPUNIT_ASSERT_MESSAGE("Calculated value does not fit expected value", std::abs(mppObject - MPP) < mitk::eps); } CPPUNIT_ASSERT_MESSAGE("Calculated value does not fit expected value", std::abs(skewnessObject - skewness) < mitk::eps); CPPUNIT_ASSERT_MESSAGE("Calculated value does not fit expected value", std::abs(kurtosisObject - kurtosis) < mitk::eps); CPPUNIT_ASSERT_MESSAGE("Calculated value does not fit expected value", std::abs(varianceObject - variance) < mitk::eps); CPPUNIT_ASSERT_MESSAGE("Calculated value does not fit expected value", std::abs(standardDeviationObject - stdev) < mitk::eps); CPPUNIT_ASSERT_MESSAGE("Calculated value does not fit expected value", std::abs(minObject - min) < mitk::eps); CPPUNIT_ASSERT_MESSAGE("Calculated value does not fit expected value", std::abs(maxObject - max) < mitk::eps); CPPUNIT_ASSERT_MESSAGE("Calculated value does not fit expected value", std::abs(rmsObject - RMS) < mitk::eps); for (unsigned int i = 0; i < minIndex.size(); ++i) { CPPUNIT_ASSERT_MESSAGE("Calculated value does not fit expected value", std::abs(minIndexObject[i] - minIndex[i]) < mitk::eps); } for (unsigned int i = 0; i < maxIndex.size(); ++i) { CPPUNIT_ASSERT_MESSAGE("Calculated value does not fit expected value", std::abs(maxIndexObject[i] - maxIndex[i]) < mitk::eps); } } MITK_TEST_SUITE_REGISTRATION(mitkImageStatisticsCalculator) diff --git a/Modules/ImageStatistics/Testing/mitkImageStatisticsContainerTest.cpp b/Modules/ImageStatistics/Testing/mitkImageStatisticsContainerTest.cpp index 886f09458e..8ee6b3fba8 100644 --- a/Modules/ImageStatistics/Testing/mitkImageStatisticsContainerTest.cpp +++ b/Modules/ImageStatistics/Testing/mitkImageStatisticsContainerTest.cpp @@ -1,365 +1,365 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ // Testing #include "mitkTestingMacros.h" #include "mitkTestFixture.h" //MITK includes #include "mitkArbitraryTimeGeometry.h" #include "mitkImageStatisticsContainerNodeHelper.h" #include "mitkImageStatisticsContainerManager.h" #include "mitkImageStatisticsContainer.h" class mitkImageStatisticsContainerTestSuite : public mitk::TestFixture { CPPUNIT_TEST_SUITE(mitkImageStatisticsContainerTestSuite); MITK_TEST(StatisticNamesIO); MITK_TEST(InvalidKey); MITK_TEST(TimeSteps); MITK_TEST(PrintSelf); MITK_TEST(InternalClone); MITK_TEST(StatisticNames); MITK_TEST(OverwriteStatistic); MITK_TEST(AllStatisticNamesForContainer); MITK_TEST(Reset); CPPUNIT_TEST_SUITE_END(); private: mitk::ImageStatisticsContainer::Pointer m_StatisticsContainer; mitk::ImageStatisticsContainer::Pointer m_StatisticsContainer2; mitk::ImageStatisticsContainer::Pointer m_StatisticsContainer3; mitk::ImageStatisticsContainer::ImageStatisticsObject m_StatisticsObject; mitk::ImageStatisticsContainer::ImageStatisticsObject m_StatisticsObject2; mitk::ImageStatisticsContainer::ImageStatisticsObject m_StatisticsObject3; mitk::ArbitraryTimeGeometry::Pointer m_TimeGeometry; mitk::Geometry3D::Pointer m_Geometry1; mitk::Geometry3D::Pointer m_Geometry2; mitk::Geometry3D::Pointer m_Geometry3; mitk::Geometry3D::Pointer m_Geometry3_5; mitk::Geometry3D::Pointer m_Geometry4; mitk::Geometry3D::Pointer m_Geometry5; mitk::TimePointType m_Geometry1MinTP; mitk::TimePointType m_Geometry2MinTP; mitk::TimePointType m_Geometry3MinTP; mitk::TimePointType m_Geometry3_5MinTP; mitk::TimePointType m_Geometry4MinTP; mitk::TimePointType m_Geometry5MinTP; mitk::TimePointType m_Geometry1MaxTP; mitk::TimePointType m_Geometry2MaxTP; mitk::TimePointType m_Geometry3MaxTP; mitk::TimePointType m_Geometry3_5MaxTP; mitk::TimePointType m_Geometry4MaxTP; mitk::TimePointType m_Geometry5MaxTP; std::vector estimatedDefaultStatisticNames { "Mean", "Median", "StandardDeviation", "RMS", "Max", "MaxPosition", "Min", "MinPosition", "#Voxel", "Volume [mm^3]", "Skewness", "Kurtosis", "Uniformity", "Entropy", "MPP", "UPP" }; public: void setUp() override { m_StatisticsContainer = mitk::ImageStatisticsContainer::New(); m_StatisticsContainer2 = mitk::ImageStatisticsContainer::New(); m_StatisticsContainer3 = mitk::ImageStatisticsContainer::New(); m_Geometry1 = mitk::Geometry3D::New(); m_Geometry2 = mitk::Geometry3D::New(); m_Geometry3 = mitk::Geometry3D::New(); m_Geometry3_5 = mitk::Geometry3D::New(); m_Geometry4 = mitk::Geometry3D::New(); m_Geometry5 = mitk::Geometry3D::New(); m_Geometry1MinTP = 1; m_Geometry2MinTP = 2; m_Geometry3MinTP = 3; m_Geometry3_5MinTP = 3.5; m_Geometry4MinTP = 4; m_Geometry5MinTP = 5; m_Geometry1MaxTP = 1.9; m_Geometry2MaxTP = 2.9; m_Geometry3MaxTP = 3.9; m_Geometry3_5MaxTP = 3.9; m_Geometry4MaxTP = 4.9; m_Geometry5MaxTP = 5.9; m_TimeGeometry = mitk::ArbitraryTimeGeometry::New(); m_TimeGeometry->ClearAllGeometries(); m_TimeGeometry->AppendNewTimeStep(m_Geometry1, m_Geometry1MinTP, m_Geometry1MaxTP); m_TimeGeometry->AppendNewTimeStep(m_Geometry2, m_Geometry2MinTP, m_Geometry2MaxTP); m_TimeGeometry->AppendNewTimeStep(m_Geometry3, m_Geometry3MinTP, m_Geometry3MaxTP); m_TimeGeometry->AppendNewTimeStep(m_Geometry4, m_Geometry4MinTP, m_Geometry4MaxTP); m_TimeGeometry->AppendNewTimeStep(m_Geometry5, m_Geometry5MinTP, m_Geometry5MaxTP); } void tearDown() override { m_StatisticsContainer->Reset(); m_StatisticsContainer2->Reset(); m_StatisticsContainer3->Reset(); m_StatisticsObject.Reset(); m_StatisticsObject2.Reset(); m_StatisticsObject3.Reset(); m_TimeGeometry->ClearAllGeometries(); } void StatisticNamesIO() { auto defaultStatisticNames = mitk::ImageStatisticsContainer::ImageStatisticsObject::GetDefaultStatisticNames(); for (size_t i = 0; i < estimatedDefaultStatisticNames.size(); i++) { CPPUNIT_ASSERT_EQUAL_MESSAGE("Estimated default statistics names are not correct.", defaultStatisticNames.at(i), estimatedDefaultStatisticNames.at(i)); } std::string testName = "testName"; double realTypeTest = 3.14; m_StatisticsObject.AddStatistic(testName, realTypeTest); auto customStatisticNames = m_StatisticsObject.GetCustomStatisticNames(); size_t one = 1; CPPUNIT_ASSERT_EQUAL_MESSAGE("Amount of estimated custom statistics names is not correct.", customStatisticNames.size(), one); CPPUNIT_ASSERT_EQUAL_MESSAGE("Estimated custom statistics names are not correct.", customStatisticNames.front(), testName); auto allStatisticNames = m_StatisticsObject.GetAllStatisticNames(); estimatedDefaultStatisticNames.push_back(testName); for (size_t i = 0; i < estimatedDefaultStatisticNames.size(); i++) { CPPUNIT_ASSERT_EQUAL_MESSAGE("Estimated statistics names are not correct.", allStatisticNames.at(i), estimatedDefaultStatisticNames.at(i)); } auto existingStatisticNames = m_StatisticsObject.GetExistingStatisticNames(); CPPUNIT_ASSERT_EQUAL_MESSAGE("Amount of estimated existing statistics names is not correct.", existingStatisticNames.size(), one); CPPUNIT_ASSERT_EQUAL_MESSAGE("Estimated existing statistics names are not correct.", existingStatisticNames.front(), testName); m_StatisticsObject.Reset(); CPPUNIT_ASSERT_MESSAGE("Custom statistics names were not removed correctly.", m_StatisticsObject.GetCustomStatisticNames().empty()); } void InvalidKey() { CPPUNIT_ASSERT_THROW_MESSAGE("Exception should have been thrown because key does not exists.", m_StatisticsObject.GetValueNonConverted("Test"), mitk::Exception); } void TimeSteps() { m_StatisticsContainer->SetTimeGeometry(m_TimeGeometry); m_StatisticsContainer->SetStatistics(1,0, m_StatisticsObject); CPPUNIT_ASSERT_EQUAL_MESSAGE("Previous added time step was not saved correctly.", m_StatisticsContainer->StatisticsExist(1, 0), true); m_StatisticsContainer->SetStatistics(1,1, m_StatisticsObject); CPPUNIT_ASSERT_EQUAL_MESSAGE("Previous added time step was not saved correctly.", m_StatisticsContainer->StatisticsExist(1, 1), true); m_StatisticsContainer->SetStatistics(1,2, m_StatisticsObject); CPPUNIT_ASSERT_EQUAL_MESSAGE("Previous added time step was not saved correctly.", m_StatisticsContainer->StatisticsExist(1, 2), true); m_StatisticsContainer->SetStatistics(1,3, m_StatisticsObject); CPPUNIT_ASSERT_EQUAL_MESSAGE("Previous added time step was not saved correctly.", m_StatisticsContainer->StatisticsExist(1, 3), true); m_StatisticsContainer->SetStatistics(1,4, m_StatisticsObject); CPPUNIT_ASSERT_EQUAL_MESSAGE("Previous added time step was not saved correctly.", m_StatisticsContainer->StatisticsExist(1, 4), true); CPPUNIT_ASSERT_THROW_MESSAGE("Out of timeStep geometry bounds timeStep was added but no exception was thrown.", m_StatisticsContainer->SetStatistics(1,42, m_StatisticsObject), mitk::Exception); CPPUNIT_ASSERT_EQUAL_MESSAGE("Statistics container does not contain the right amount of timeSteps.", static_cast (m_StatisticsContainer->GetExistingTimeSteps(1).size()), 5); CPPUNIT_ASSERT_THROW_MESSAGE("A statistic for a non existing time steps was found.", m_StatisticsContainer->GetStatistics(1, 42), mitk::Exception); } void PrintSelf() { std::stringstream print; // Check print function for empty container CPPUNIT_ASSERT_NO_THROW_MESSAGE("Print function throws an exception.", m_StatisticsContainer->Print(print)); m_StatisticsContainer->SetTimeGeometry(m_TimeGeometry); m_StatisticsObject.AddStatistic("Test", 4.2); m_StatisticsContainer->SetStatistics(1,0, m_StatisticsObject); m_StatisticsContainer->SetStatistics(1,1, m_StatisticsObject); m_StatisticsContainer->SetStatistics(1,2, m_StatisticsObject); m_StatisticsContainer->SetStatistics(1,3, m_StatisticsObject); m_StatisticsContainer->SetStatistics(1,4, m_StatisticsObject); // Check print function for filled container CPPUNIT_ASSERT_NO_THROW_MESSAGE("Print function throws an exception.", m_StatisticsContainer->Print(print)); } void InternalClone() { m_StatisticsContainer->SetTimeGeometry(m_TimeGeometry); m_StatisticsObject.AddStatistic("Test", 4.2); m_StatisticsContainer->SetStatistics(1,0, m_StatisticsObject); m_StatisticsContainer->SetStatistics(1,1, m_StatisticsObject); m_StatisticsContainer->SetStatistics(1,2, m_StatisticsObject); m_StatisticsContainer->SetStatistics(1,3, m_StatisticsObject); m_StatisticsContainer->SetStatistics(1,4, m_StatisticsObject); auto clone = m_StatisticsContainer->Clone(); CPPUNIT_ASSERT_EQUAL_MESSAGE("Internal clone was not cloned correctly.", m_StatisticsContainer->GetExistingTimeSteps(1).size(), clone->GetExistingTimeSteps(1).size()); CPPUNIT_ASSERT_EQUAL_MESSAGE("Internal clone was not cloned correctly.", clone->StatisticsExist(1, 0), true); CPPUNIT_ASSERT_EQUAL_MESSAGE("Internal clone was not cloned correctly.", clone->StatisticsExist(1, 1), true); CPPUNIT_ASSERT_EQUAL_MESSAGE("Internal clone was not cloned correctly.", clone->StatisticsExist(1, 2), true); CPPUNIT_ASSERT_EQUAL_MESSAGE("Internal clone was not cloned correctly.", clone->StatisticsExist(1, 3), true); CPPUNIT_ASSERT_EQUAL_MESSAGE("Internal clone was not cloned correctly.", clone->StatisticsExist(1, 4), true); CPPUNIT_ASSERT_EQUAL_MESSAGE("Internal clone was not cloned correctly.", clone->GetStatistics(1,0).HasStatistic("Test"), true); CPPUNIT_ASSERT_EQUAL_MESSAGE("Internal clone was not cloned correctly.", clone->GetStatistics(1,1).HasStatistic("Test"), true); CPPUNIT_ASSERT_EQUAL_MESSAGE("Internal clone was not cloned correctly.", clone->GetStatistics(1,2).HasStatistic("Test"), true); CPPUNIT_ASSERT_EQUAL_MESSAGE("Internal clone was not cloned correctly.", clone->GetStatistics(1,3).HasStatistic("Test"), true); CPPUNIT_ASSERT_EQUAL_MESSAGE("Internal clone was not cloned correctly.", clone->GetStatistics(1,4).HasStatistic("Test"), true); } void StatisticNames() { m_StatisticsContainer->SetTimeGeometry(m_TimeGeometry); m_StatisticsContainer2->SetTimeGeometry(m_TimeGeometry); m_StatisticsContainer3->SetTimeGeometry(m_TimeGeometry); m_StatisticsObject.AddStatistic("Test", 4.2); m_StatisticsObject2.AddStatistic("Test2", 4.2); m_StatisticsObject3.AddStatistic("Test3", 4.2); // Setting the same statistics for different time steps will only add the statistics name once. m_StatisticsContainer->SetStatistics(1,0, m_StatisticsObject); m_StatisticsContainer->SetStatistics(1,1, m_StatisticsObject); m_StatisticsContainer->SetStatistics(1,2, m_StatisticsObject); m_StatisticsContainer->SetStatistics(1,3, m_StatisticsObject); m_StatisticsContainer->SetStatistics(1,4, m_StatisticsObject); // Setting the same statistics for different time steps will only add the statistics name once. m_StatisticsContainer2->SetStatistics(1,0, m_StatisticsObject2); m_StatisticsContainer2->SetStatistics(1,1, m_StatisticsObject2); m_StatisticsContainer2->SetStatistics(1,2, m_StatisticsObject2); m_StatisticsContainer2->SetStatistics(1,3, m_StatisticsObject2); m_StatisticsContainer2->SetStatistics(1,4, m_StatisticsObject2); // Setting the same statistics for different time steps will only add the statistics name once. m_StatisticsContainer3->SetStatistics(1,0, m_StatisticsObject3); m_StatisticsContainer3->SetStatistics(1,1, m_StatisticsObject3); m_StatisticsContainer3->SetStatistics(1,2, m_StatisticsObject3); m_StatisticsContainer3->SetStatistics(1,3, m_StatisticsObject3); m_StatisticsContainer3->SetStatistics(1,4, m_StatisticsObject3); CPPUNIT_ASSERT_EQUAL_MESSAGE("Estimated statistics names are not correct.", mitk::GetAllStatisticNames(m_StatisticsContainer).size(), estimatedDefaultStatisticNames.size() + 1); CPPUNIT_ASSERT_EQUAL_MESSAGE("Custom statistic name was not saved correctly.", m_StatisticsContainer->GetStatistics(1,0).HasStatistic("Test"), true); std::vector containers; containers.push_back(m_StatisticsContainer2.GetPointer()); containers.push_back(m_StatisticsContainer3.GetPointer()); CPPUNIT_ASSERT_EQUAL_MESSAGE("Estimated statistics names are not correct.", mitk::GetAllStatisticNames(containers).size(), estimatedDefaultStatisticNames.size() + 2); CPPUNIT_ASSERT_EQUAL_MESSAGE("Custom statistic name was not saved correctly.", m_StatisticsContainer2->GetStatistics(1,0).HasStatistic("Test2"), true); CPPUNIT_ASSERT_EQUAL_MESSAGE("Custom statistic name was not saved correctly.", m_StatisticsContainer3->GetStatistics(1,0).HasStatistic("Test3"), true); } void OverwriteStatistic() { m_StatisticsContainer->SetTimeGeometry(m_TimeGeometry); m_StatisticsObject.AddStatistic("Test", 4.2); m_StatisticsContainer->SetStatistics(1,0, m_StatisticsObject); m_StatisticsContainer->SetStatistics(1,1, m_StatisticsObject); m_StatisticsContainer->SetStatistics(1,2, m_StatisticsObject); m_StatisticsContainer->SetStatistics(1,3, m_StatisticsObject); m_StatisticsContainer->SetStatistics(1,4, m_StatisticsObject); CPPUNIT_ASSERT_EQUAL_MESSAGE("Statistic was not correctly added to StatisticsObject.", boost::get (m_StatisticsContainer->GetStatistics(1,0).GetValueNonConverted("Test")), 4.2); // An existing statistic won't be updated by adding another statistic with same name to that object. m_StatisticsObject.AddStatistic("Test", 42.0); CPPUNIT_ASSERT_EQUAL_MESSAGE("Statistic was overwritten.", boost::get(m_StatisticsContainer->GetStatistics(1,0).GetValueNonConverted("Test")), 4.2); } void AllStatisticNamesForContainer() { m_StatisticsContainer->SetTimeGeometry(m_TimeGeometry); m_StatisticsObject.AddStatistic("Test", 4.2); m_StatisticsObject2.AddStatistic("Test2", 4.2); // Setting the same statistics for different time steps will only add the statistics name once. m_StatisticsContainer->SetStatistics(1,0, m_StatisticsObject); m_StatisticsContainer->SetStatistics(1,1, m_StatisticsObject2); m_StatisticsContainer->SetStatistics(1,2, m_StatisticsObject); m_StatisticsContainer->SetStatistics(1,3, m_StatisticsObject2); m_StatisticsContainer->SetStatistics(1,4, m_StatisticsObject); std::vector containers; // Check if empty vector triggers exception. CPPUNIT_ASSERT_NO_THROW_MESSAGE("Empty vector triggered an exception.", mitk::GetAllStatisticNames(containers)); containers.push_back(m_StatisticsContainer2.GetPointer()); containers.push_back(m_StatisticsContainer3.GetPointer()); // Check if vector with empty container triggers exception. CPPUNIT_ASSERT_NO_THROW_MESSAGE("Empty container triggered an exception.", mitk::GetAllStatisticNames(containers)); // Adding filled container to check if data was saved correctly. containers.push_back(m_StatisticsContainer.GetPointer()); CPPUNIT_ASSERT_EQUAL_MESSAGE("Estimated statistics names are not correct.", mitk::GetAllStatisticNames(containers).size(), estimatedDefaultStatisticNames.size() + 2); CPPUNIT_ASSERT_EQUAL_MESSAGE("Custom statistic name was not saved correctly.", m_StatisticsContainer->GetStatistics(1,0).HasStatistic("Test"), true); CPPUNIT_ASSERT_EQUAL_MESSAGE("Custom statistic name was not saved correctly.", m_StatisticsContainer->GetStatistics(1,1).HasStatistic("Test2"), true); } void Reset() { m_StatisticsContainer->SetTimeGeometry(m_TimeGeometry); m_StatisticsObject.AddStatistic("Test", 4.2); m_StatisticsContainer->SetStatistics(1,0, m_StatisticsObject); m_StatisticsContainer->SetStatistics(1,1, m_StatisticsObject); m_StatisticsContainer->SetStatistics(1,2, m_StatisticsObject); m_StatisticsContainer->SetStatistics(1,3, m_StatisticsObject); m_StatisticsContainer->SetStatistics(1,4, m_StatisticsObject); m_StatisticsContainer->Reset(); m_StatisticsObject.Reset(); - CPPUNIT_ASSERT_MESSAGE("Statistics container was reseted correctly.", m_StatisticsContainer->GetExistingLabelValues(false).empty()); + CPPUNIT_ASSERT_MESSAGE("Statistics container was reseted correctly.", m_StatisticsContainer->GetExistingLabelValues().empty()); CPPUNIT_ASSERT_EQUAL_MESSAGE("Statistic was not deleted correctly.", m_StatisticsObject.HasStatistic("Test"), false); } }; MITK_TEST_SUITE_REGISTRATION(mitkImageStatisticsContainer) diff --git a/Modules/ImageStatistics/Testing/mitkImageStatisticsHotspotTest.cpp b/Modules/ImageStatistics/Testing/mitkImageStatisticsHotspotTest.cpp index 0e9198f5ce..3114c2b52c 100644 --- a/Modules/ImageStatistics/Testing/mitkImageStatisticsHotspotTest.cpp +++ b/Modules/ImageStatistics/Testing/mitkImageStatisticsHotspotTest.cpp @@ -1,643 +1,643 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #include "mitkImageStatisticsCalculator.h" #include "itkMultiGaussianImageSource.h" #include "mitkTestingMacros.h" #include #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) { typedef double PixelType; const int Dimension = 3; typedef itk::Image ImageType; 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(); return mitk::GrabItkImageMemory(gaussianGenerator->GetOutput(), nullptr, nullptr, false); } /** \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::ImageStatisticsObject 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); 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]); } } } auto mitkMaskImage = mitk::GrabItkImageMemory(mask, nullptr, nullptr, false); mitk::ImageMaskGenerator::Pointer imgMaskGen = mitk::ImageMaskGenerator::New(); imgMaskGen->SetInputImage(image); 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()->GetStatisticsForTimeStep(0); + return statisticsCalculator->GetStatistics()->GetStatistics(label,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::ImageStatisticsObject 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::ImageStatisticsObject 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/mitkHotspotMaskGenerator.cpp b/Modules/ImageStatistics/mitkHotspotMaskGenerator.cpp index e36b93dce5..76637bfe5d 100644 --- a/Modules/ImageStatistics/mitkHotspotMaskGenerator.cpp +++ b/Modules/ImageStatistics/mitkHotspotMaskGenerator.cpp @@ -1,621 +1,542 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #include #include #include #include #include #include "mitkImageAccessByItk.h" #include #include #include #include namespace mitk { HotspotMaskGenerator::HotspotMaskGenerator(): - m_HotspotRadiusinMM(6.2035049089940), // radius of a 1cm3 sphere in mm + m_HotspotRadiusInMM(6.2035049089940), // radius of a 1cm3 sphere in mm m_HotspotMustBeCompletelyInsideImage(true), m_Label(1) { - m_TimeStep = 0; m_InternalMask = mitk::Image::New(); m_InternalMaskUpdateTime = 0; } - void HotspotMaskGenerator::SetInputImage(mitk::Image::Pointer inputImage) - { - if (inputImage != m_inputImage) - { - m_inputImage = inputImage; - m_ConvolutionImageMaxIndex.set_size(inputImage->GetDimension()); - m_ConvolutionImageMinIndex.set_size(inputImage->GetDimension()); - this->Modified(); - } - } - - void HotspotMaskGenerator::SetMask(MaskGenerator::Pointer mask) - { - if (mask != m_Mask) - { - m_Mask = mask; - this->Modified(); - } - } - HotspotMaskGenerator::~HotspotMaskGenerator() { } - void HotspotMaskGenerator::SetHotspotRadiusInMM(double radiusInMillimeter) + unsigned int HotspotMaskGenerator::GetNumberOfMasks() const { - if(radiusInMillimeter != m_HotspotRadiusinMM) - { - m_HotspotRadiusinMM = radiusInMillimeter; - this->Modified(); - } - } - - const double& HotspotMaskGenerator::GetHotspotRadiusinMM() const - { - return m_HotspotRadiusinMM; + return 1; } - bool HotspotMaskGenerator::GetHotspotMustBeCompletelyInsideImage() const - { - return m_HotspotMustBeCompletelyInsideImage; - } - - void HotspotMaskGenerator::SetHotspotMustBeCompletelyInsideImage(bool mustBeCompletelyInImage) - { - if (m_HotspotMustBeCompletelyInsideImage != mustBeCompletelyInImage) - { - m_HotspotMustBeCompletelyInsideImage = mustBeCompletelyInImage; - this->Modified(); - } - } - - - mitk::Image::ConstPointer HotspotMaskGenerator::GetMask() + mitk::Image::ConstPointer HotspotMaskGenerator::DoGetMask(unsigned int) { if (IsUpdateRequired()) { - if ( m_inputImage.IsNull() ) + if ( m_InputImage.IsNull() ) { throw std::runtime_error( "Error: image empty!" ); } - if ( m_TimeStep >= m_inputImage->GetTimeSteps() ) + if ( m_InputImage->GetTimeGeometry()->IsValidTimePoint(m_TimePoint) ) { throw std::runtime_error( "Error: invalid time step!" ); } - mitk::ImageTimeSelector::Pointer imageTimeSelector = mitk::ImageTimeSelector::New(); - imageTimeSelector->SetInput( m_inputImage ); - imageTimeSelector->SetTimeNr( m_TimeStep ); - imageTimeSelector->UpdateLargestPossibleRegion(); - mitk::Image::Pointer timeSliceImage = imageTimeSelector->GetOutput(); + auto timeSliceImage = SelectImageByTimePoint(m_InputImage, m_TimePoint); - m_internalImage = timeSliceImage; m_internalMask2D = nullptr; // is this correct when this variable holds a smart pointer? m_internalMask3D = nullptr; if ( m_Mask != nullptr ) { - m_Mask->SetTimeStep(m_TimeStep); - mitk::Image::ConstPointer timeSliceMask = m_Mask->GetMask(); + m_Mask->SetTimePoint(m_TimePoint); + mitk::Image::ConstPointer timeSliceMask = m_Mask->GetMask(0); - if ( m_internalImage->GetDimension() == 3 ) + if ( timeSliceImage->GetDimension() == 3 ) { itk::Image::Pointer noneConstMaskImage; //needed to work arround the fact that CastToItkImage currently does not support const itk images. CastToItkImage(timeSliceMask, noneConstMaskImage); m_internalMask3D = noneConstMaskImage; - AccessFixedDimensionByItk_2(m_internalImage, CalculateHotspotMask, 3, m_internalMask3D.GetPointer(), m_Label); + AccessFixedDimensionByItk_2(timeSliceImage, CalculateHotspotMask, 3, m_internalMask3D.GetPointer(), m_Label); } - else if ( m_internalImage->GetDimension() == 2 ) + else if ( timeSliceImage->GetDimension() == 2 ) { itk::Image::Pointer noneConstMaskImage; //needed to work arround the fact that CastToItkImage currently does not support const itk images. CastToItkImage(timeSliceMask, noneConstMaskImage); m_internalMask2D = noneConstMaskImage; - AccessFixedDimensionByItk_2(m_internalImage, CalculateHotspotMask, 2, m_internalMask2D.GetPointer(), m_Label); + AccessFixedDimensionByItk_2(timeSliceImage, CalculateHotspotMask, 2, m_internalMask2D.GetPointer(), m_Label); } else { throw std::runtime_error( "Error: invalid image dimension" ); } } else { - if ( m_internalImage->GetDimension() == 3 ) + if ( timeSliceImage->GetDimension() == 3 ) { - AccessFixedDimensionByItk_2(m_internalImage, CalculateHotspotMask, 3, m_internalMask3D.GetPointer(), m_Label); + AccessFixedDimensionByItk_2(timeSliceImage, CalculateHotspotMask, 3, m_internalMask3D.GetPointer(), m_Label); } - else if ( m_internalImage->GetDimension() == 2 ) + else if ( timeSliceImage->GetDimension() == 2 ) { - AccessFixedDimensionByItk_2(m_internalImage, CalculateHotspotMask, 2, m_internalMask2D.GetPointer(), m_Label); + AccessFixedDimensionByItk_2(timeSliceImage, CalculateHotspotMask, 2, m_internalMask2D.GetPointer(), m_Label); } else { throw std::runtime_error( "Error: invalid image dimension" ); } } this->Modified(); } m_InternalMaskUpdateTime = m_InternalMask->GetMTime(); return m_InternalMask; } - void HotspotMaskGenerator::SetTimeStep(unsigned int timeStep) - { - if (m_TimeStep != timeStep) - { - m_TimeStep = timeStep; - } - } - - void HotspotMaskGenerator::SetLabel(unsigned short label) - { - if (label != m_Label) - { - m_Label = label; - this->Modified(); - } - } - - vnl_vector HotspotMaskGenerator::GetConvolutionImageMinIndex() - { - this->GetMask(); // make sure we are up to date - return m_ConvolutionImageMinIndex; - } - - vnl_vector HotspotMaskGenerator::GetHotspotIndex() - { - this->GetMask(); // make sure we are up to date - return m_ConvolutionImageMaxIndex; - } - template HotspotMaskGenerator::ImageExtrema HotspotMaskGenerator::CalculateExtremaWorld( const itk::Image* inputImage, const itk::Image* maskImage, double neccessaryDistanceToImageBorderInMM, unsigned int label ) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< unsigned short, VImageDimension > MaskImageType; typedef itk::ImageRegionConstIteratorWithIndex MaskImageIteratorType; typedef itk::ImageRegionConstIteratorWithIndex InputImageIndexIteratorType; typename ImageType::SpacingType spacing = inputImage->GetSpacing(); ImageExtrema minMax; minMax.Defined = false; minMax.MaxIndex.set_size(VImageDimension); minMax.MaxIndex.set_size(VImageDimension); typename ImageType::RegionType allowedExtremaRegion = inputImage->GetLargestPossibleRegion(); bool keepDistanceToImageBorders( neccessaryDistanceToImageBorderInMM > 0 ); if (keepDistanceToImageBorders) { itk::IndexValueType distanceInPixels[VImageDimension]; for(unsigned short dimension = 0; dimension < VImageDimension; ++dimension) { // To confirm that the whole hotspot is inside the image we have to keep a specific distance to the image-borders, which is as long as // the radius. To get the amount of indices we divide the radius by spacing and add 0.5 because voxels are center based: // For example with a radius of 2.2 and a spacing of 1 two indices are enough because 2.2 / 1 + 0.5 = 2.7 => 2. // But with a radius of 2.7 we need 3 indices because 2.7 / 1 + 0.5 = 3.2 => 3 distanceInPixels[dimension] = int( neccessaryDistanceToImageBorderInMM / spacing[dimension] + 0.5); } allowedExtremaRegion.ShrinkByRadius(distanceInPixels); } InputImageIndexIteratorType imageIndexIt(inputImage, allowedExtremaRegion); float maxValue = itk::NumericTraits::min(); float minValue = itk::NumericTraits::max(); typename ImageType::IndexType maxIndex; typename ImageType::IndexType minIndex; for(unsigned short i = 0; i < VImageDimension; ++i) { maxIndex[i] = 0; minIndex[i] = 0; } if (maskImage != nullptr) { MaskImageIteratorType maskIt(maskImage, maskImage->GetLargestPossibleRegion()); typename ImageType::IndexType imageIndex; typename ImageType::IndexType maskIndex; for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt) { imageIndex = maskIndex = maskIt.GetIndex(); if(maskIt.Get() == label) { if( allowedExtremaRegion.IsInside(imageIndex) ) { imageIndexIt.SetIndex( imageIndex ); double value = imageIndexIt.Get(); minMax.Defined = true; //Calculate minimum, maximum and corresponding index-values if( value > maxValue ) { maxIndex = imageIndexIt.GetIndex(); maxValue = value; } if(value < minValue ) { minIndex = imageIndexIt.GetIndex(); minValue = value; } } } } } else { for(imageIndexIt.GoToBegin(); !imageIndexIt.IsAtEnd(); ++imageIndexIt) { double value = imageIndexIt.Get(); minMax.Defined = true; //Calculate minimum, maximum and corresponding index-values if( value > maxValue ) { maxIndex = imageIndexIt.GetIndex(); maxValue = value; } if(value < minValue ) { minIndex = imageIndexIt.GetIndex(); minValue = value; } } } minMax.MaxIndex.set_size(VImageDimension); minMax.MinIndex.set_size(VImageDimension); for(unsigned int i = 0; i < minMax.MaxIndex.size(); ++i) { minMax.MaxIndex[i] = maxIndex[i]; } for(unsigned int i = 0; i < minMax.MinIndex.size(); ++i) { minMax.MinIndex[i] = minIndex[i]; } minMax.Max = maxValue; minMax.Min = minValue; return minMax; } template itk::Size HotspotMaskGenerator::CalculateConvolutionKernelSize( double spacing[VImageDimension], double radiusInMM ) { typedef itk::Image< float, VImageDimension > KernelImageType; typedef typename KernelImageType::SizeType SizeType; SizeType maskSize; for(unsigned int i = 0; i < VImageDimension; ++i) { maskSize[i] = static_cast( 2 * radiusInMM / spacing[i]); // We always want an uneven size to have a clear center point in the convolution mask if(maskSize[i] % 2 == 0 ) { ++maskSize[i]; } } return maskSize; } template itk::SmartPointer< itk::Image > HotspotMaskGenerator::GenerateHotspotSearchConvolutionKernel(double mmPerPixel[VImageDimension], double radiusInMM ) { std::stringstream ss; for (unsigned int i = 0; i < VImageDimension; ++i) { ss << mmPerPixel[i]; if (i < VImageDimension -1) ss << ","; } MITK_DEBUG << "Update convolution kernel for spacing (" << ss.str() << ") and radius " << radiusInMM << "mm"; double radiusInMMSquared = radiusInMM * radiusInMM; typedef itk::Image< float, VImageDimension > KernelImageType; typename KernelImageType::Pointer convolutionKernel = KernelImageType::New(); // Calculate size and allocate mask image typedef typename KernelImageType::SizeType SizeType; SizeType maskSize = this->CalculateConvolutionKernelSize(mmPerPixel, radiusInMM); mitk::Point3D convolutionMaskCenterIndex; convolutionMaskCenterIndex.Fill(0.0); for(unsigned int i = 0; i < VImageDimension; ++i) { convolutionMaskCenterIndex[i] = 0.5 * (double)(maskSize[i]-1); } typedef typename KernelImageType::IndexType IndexType; IndexType maskIndex; maskIndex.Fill(0); typedef typename KernelImageType::RegionType RegionType; RegionType maskRegion; maskRegion.SetSize(maskSize); maskRegion.SetIndex(maskIndex); convolutionKernel->SetRegions(maskRegion); convolutionKernel->SetSpacing(mmPerPixel); convolutionKernel->Allocate(); // Fill mask image values by subsampling the image grid typedef itk::ImageRegionIteratorWithIndex MaskIteratorType; MaskIteratorType maskIt(convolutionKernel,maskRegion); int numberOfSubVoxelsPerDimension = 2; // per dimension! int numberOfSubVoxels = ::pow( static_cast(numberOfSubVoxelsPerDimension), static_cast(VImageDimension) ); double subVoxelSizeInPixels = 1.0 / (double)numberOfSubVoxelsPerDimension; double valueOfOneSubVoxel = 1.0 / (double)numberOfSubVoxels; mitk::Point3D subVoxelIndexPosition; double distanceSquared = 0.0; typedef itk::ContinuousIndex ContinuousIndexType; for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt) { ContinuousIndexType indexPoint(maskIt.GetIndex()); mitk::Point3D voxelPosition; for (unsigned int dimension = 0; dimension < VImageDimension; ++dimension) { voxelPosition[dimension] = indexPoint[dimension]; } double maskValue = 0.0; mitk::Vector3D subVoxelOffset; subVoxelOffset.Fill(0.0); // iterate sub-voxels by iterating all possible offsets for (subVoxelOffset[0] = -0.5 + subVoxelSizeInPixels / 2.0; subVoxelOffset[0] < +0.5; subVoxelOffset[0] += subVoxelSizeInPixels) { for (subVoxelOffset[1] = -0.5 + subVoxelSizeInPixels / 2.0; subVoxelOffset[1] < +0.5; subVoxelOffset[1] += subVoxelSizeInPixels) { for (subVoxelOffset[2] = -0.5 + subVoxelSizeInPixels / 2.0; subVoxelOffset[2] < +0.5; subVoxelOffset[2] += subVoxelSizeInPixels) { subVoxelIndexPosition = voxelPosition + subVoxelOffset; // this COULD be integrated into the for-loops if neccessary (add voxelPosition to initializer and end condition) distanceSquared = (subVoxelIndexPosition[0]-convolutionMaskCenterIndex[0]) * mmPerPixel[0] * (subVoxelIndexPosition[0]-convolutionMaskCenterIndex[0]) * mmPerPixel[0] + (subVoxelIndexPosition[1]-convolutionMaskCenterIndex[1]) * mmPerPixel[1] * (subVoxelIndexPosition[1]-convolutionMaskCenterIndex[1]) * mmPerPixel[1] + (subVoxelIndexPosition[2]-convolutionMaskCenterIndex[2]) * mmPerPixel[2] * (subVoxelIndexPosition[2]-convolutionMaskCenterIndex[2]) * mmPerPixel[2]; if (distanceSquared <= radiusInMMSquared) { maskValue += valueOfOneSubVoxel; } } } } maskIt.Set( maskValue ); } return convolutionKernel; } template itk::SmartPointer > HotspotMaskGenerator::GenerateConvolutionImage( const itk::Image* inputImage ) { double mmPerPixel[VImageDimension]; for (unsigned int dimension = 0; dimension < VImageDimension; ++dimension) { mmPerPixel[dimension] = inputImage->GetSpacing()[dimension]; } // update convolution kernel typedef itk::Image< float, VImageDimension > KernelImageType; - typename KernelImageType::Pointer convolutionKernel = this->GenerateHotspotSearchConvolutionKernel(mmPerPixel, m_HotspotRadiusinMM); + typename KernelImageType::Pointer convolutionKernel = this->GenerateHotspotSearchConvolutionKernel(mmPerPixel, m_HotspotRadiusInMM); // update convolution image typedef itk::Image< TPixel, VImageDimension > InputImageType; typedef itk::Image< TPixel, VImageDimension > ConvolutionImageType; typedef itk::FFTConvolutionImageFilter ConvolutionFilterType; itk::VnlFFTImageFilterInitFactory::RegisterFactories(); typename ConvolutionFilterType::Pointer convolutionFilter = ConvolutionFilterType::New(); typedef itk::ConstantBoundaryCondition BoundaryConditionType; BoundaryConditionType boundaryCondition; boundaryCondition.SetConstant(0.0); if (m_HotspotMustBeCompletelyInsideImage) { // overwrite default boundary condition convolutionFilter->SetBoundaryCondition(&boundaryCondition); } convolutionFilter->SetInput(inputImage); convolutionFilter->SetKernelImage(convolutionKernel); convolutionFilter->SetNormalize(true); MITK_DEBUG << "Update Convolution image for hotspot search"; convolutionFilter->UpdateLargestPossibleRegion(); typename ConvolutionImageType::Pointer convolutionImage = convolutionFilter->GetOutput(); convolutionImage->SetSpacing( inputImage->GetSpacing() ); // only workaround because convolution filter seems to ignore spacing of input image return convolutionImage; } template < typename TPixel, unsigned int VImageDimension> void HotspotMaskGenerator::FillHotspotMaskPixels( itk::Image* maskImage, itk::Point sphereCenter, double sphereRadiusInMM ) { typedef itk::Image< TPixel, VImageDimension > MaskImageType; typedef itk::ImageRegionIteratorWithIndex MaskImageIteratorType; MaskImageIteratorType maskIt(maskImage, maskImage->GetLargestPossibleRegion()); typename MaskImageType::IndexType maskIndex; typename MaskImageType::PointType worldPosition; // this is not very smart. I would rather use a 0 initialized mask (not the case here -> blame CalculateHotspotMask) and find the region where I need to iterate over, then iterate only over the small region for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt) { maskIndex = maskIt.GetIndex(); maskImage->TransformIndexToPhysicalPoint(maskIndex, worldPosition); maskIt.Set( worldPosition.EuclideanDistanceTo(sphereCenter) <= sphereRadiusInMM ? 1 : 0 ); } } template void - HotspotMaskGenerator::CalculateHotspotMask(itk::Image* inputImage, + HotspotMaskGenerator::CalculateHotspotMask(const itk::Image* inputImage, const itk::Image* maskImage, unsigned int label) { typedef itk::Image< TPixel, VImageDimension > InputImageType; typedef itk::Image< TPixel, VImageDimension > ConvolutionImageType; typedef itk::Image< unsigned short, VImageDimension > MaskImageType; typename ConvolutionImageType::Pointer convolutionImage = this->GenerateConvolutionImage(inputImage); if (convolutionImage.IsNull()) { MITK_ERROR << "Empty convolution image in CalculateHotspotStatistics(). We should never reach this state (logic error)."; throw std::logic_error("Empty convolution image in CalculateHotspotStatistics()"); } typename MaskImageType::ConstPointer usedMask = maskImage; // if mask image is not defined, create an image of the same size as inputImage and fill it with 1's // there is maybe a better way to do this!? if (maskImage == nullptr) { auto defaultMask = MaskImageType::New(); typename MaskImageType::RegionType maskRegion = inputImage->GetLargestPossibleRegion(); typename MaskImageType::SpacingType maskSpacing = inputImage->GetSpacing(); typename MaskImageType::PointType maskOrigin = inputImage->GetOrigin(); typename MaskImageType::DirectionType maskDirection = inputImage->GetDirection(); defaultMask->SetRegions(maskRegion); defaultMask->Allocate(); defaultMask->SetOrigin(maskOrigin); defaultMask->SetSpacing(maskSpacing); defaultMask->SetDirection(maskDirection); defaultMask->FillBuffer(1); usedMask = defaultMask; label = 1; } // find maximum in convolution image, given the current mask - double requiredDistanceToBorder = m_HotspotMustBeCompletelyInsideImage ? m_HotspotRadiusinMM : -1.0; + double requiredDistanceToBorder = m_HotspotMustBeCompletelyInsideImage ? m_HotspotRadiusInMM : -1.0; ImageExtrema convolutionImageInformation = CalculateExtremaWorld(convolutionImage.GetPointer(), usedMask.GetPointer(), requiredDistanceToBorder, label); bool isHotspotDefined = convolutionImageInformation.Defined; if (!isHotspotDefined) { MITK_ERROR << "No origin of hotspot-sphere was calculated!"; m_InternalMask = nullptr; } else { // create a binary mask around the "hotspot" region, fill the shape of a sphere around our hotspot center // typename DuplicatorType::Pointer copyMachine = DuplicatorType::New(); // copyMachine->SetInputImage(inputImage); // copyMachine->Update(); // typename CastFilterType::Pointer caster = CastFilterType::New(); // caster->SetInput( copyMachine->GetOutput() ); // caster->Update(); typename MaskImageType::Pointer hotspotMaskITK = MaskImageType::New(); hotspotMaskITK->SetOrigin(inputImage->GetOrigin()); hotspotMaskITK->SetSpacing(inputImage->GetSpacing()); hotspotMaskITK->SetLargestPossibleRegion(inputImage->GetLargestPossibleRegion()); hotspotMaskITK->SetBufferedRegion(inputImage->GetBufferedRegion()); hotspotMaskITK->SetDirection(inputImage->GetDirection()); hotspotMaskITK->SetNumberOfComponentsPerPixel(inputImage->GetNumberOfComponentsPerPixel()); hotspotMaskITK->Allocate(); hotspotMaskITK->FillBuffer(1); typedef typename InputImageType::IndexType IndexType; IndexType maskCenterIndex; for (unsigned int d =0; d< VImageDimension;++d) { maskCenterIndex[d]=convolutionImageInformation.MaxIndex[d]; } typename ConvolutionImageType::PointType maskCenter; inputImage->TransformIndexToPhysicalPoint(maskCenterIndex,maskCenter); - FillHotspotMaskPixels(hotspotMaskITK.GetPointer(), maskCenter, m_HotspotRadiusinMM); + FillHotspotMaskPixels(hotspotMaskITK.GetPointer(), maskCenter, m_HotspotRadiusInMM); //obtain mitk::Image::Pointer from itk::Image mitk::Image::Pointer hotspotMaskAsMITKImage = mitk::GrabItkImageMemory(hotspotMaskITK); m_InternalMask = hotspotMaskAsMITKImage; m_ConvolutionImageMaxIndex = convolutionImageInformation.MaxIndex; m_ConvolutionImageMinIndex = convolutionImageInformation.MinIndex; } } bool HotspotMaskGenerator::IsUpdateRequired() const { unsigned long thisClassTimeStamp = this->GetMTime(); unsigned long internalMaskTimeStamp = m_InternalMask->GetMTime(); unsigned long maskGeneratorTimeStamp = m_Mask->GetMTime(); - unsigned long inputImageTimeStamp = m_inputImage->GetMTime(); + unsigned long inputImageTimeStamp = m_InputImage->GetMTime(); if (thisClassTimeStamp > m_InternalMaskUpdateTime) // inputs have changed { return true; } if (m_InternalMaskUpdateTime < maskGeneratorTimeStamp || m_InternalMaskUpdateTime < inputImageTimeStamp) // mask image has changed outside of this class { return true; } if (internalMaskTimeStamp > m_InternalMaskUpdateTime) // internal mask has been changed outside of this class { return true; } return false; } } diff --git a/Modules/ImageStatistics/mitkHotspotMaskGenerator.h b/Modules/ImageStatistics/mitkHotspotMaskGenerator.h index 8a9f05c538..ca388e9bde 100644 --- a/Modules/ImageStatistics/mitkHotspotMaskGenerator.h +++ b/Modules/ImageStatistics/mitkHotspotMaskGenerator.h @@ -1,178 +1,153 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #ifndef mitkHotspotMaskGenerator_h #define mitkHotspotMaskGenerator_h #include #include #include #include #include #include #include #include namespace mitk { /** * @brief The HotspotMaskGenerator class is used when a hotspot has to be found in an image. A hotspot is * the region of the image where the mean intensity is maximal (=brightest spot). It is usually used in PET scans. * The identification of the hotspot is done as follows: First a cubic (or circular, if image is 2d) * mask of predefined size is generated. This mask is then convolved with the input image (in fourier domain). * The maximum value of the convolved image then corresponds to the hotspot. * If a maskGenerator is set, only the pixels of the convolved image where the corresponding mask is == @a label * are searched for the maximum value. */ class MITKIMAGESTATISTICS_EXPORT HotspotMaskGenerator: public MaskGenerator { public: /** Standard Self typedef */ typedef HotspotMaskGenerator Self; typedef MaskGenerator 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(HotspotMaskGenerator, MaskGenerator); - /** - @brief Set the input image. Required for this class - */ - void SetInputImage(mitk::Image::Pointer inputImage); + unsigned int GetNumberOfMasks() const override; /** @brief Set a mask (can be nullptr if no mask is desired) */ - void SetMask(MaskGenerator::Pointer mask); + itkSetObjectMacro(Mask, MaskGenerator); /** @brief Set the radius of the hotspot (in MM) */ - void SetHotspotRadiusInMM(double radiusInMillimeter); - - const double& GetHotspotRadiusinMM() const; + itkGetConstMacro(HotspotRadiusInMM, double); + itkSetMacro(HotspotRadiusInMM, double); /** @brief Define whether the hotspot must be completely inside the image. Default is true */ - void SetHotspotMustBeCompletelyInsideImage(bool hotspotCompletelyInsideImage); - - bool GetHotspotMustBeCompletelyInsideImage() const; + itkGetConstMacro(HotspotMustBeCompletelyInsideImage, bool); + itkSetMacro(HotspotMustBeCompletelyInsideImage, bool); /** @brief If a maskGenerator is set, this detemines which mask value is used */ - void SetLabel(unsigned short label); - - /** - @brief Computes and returns the hotspot mask. The hotspot mask has the same size as the input image. The hopspot has value 1, the remaining pixels are set to 0 - */ - mitk::Image::ConstPointer GetMask() override; - - /** - @brief Returns the image index where the hotspot is located - */ - vnl_vector GetHotspotIndex(); - - /** - @brief Returns the index where the convolution image is minimal (darkest spot in image) - */ - vnl_vector GetConvolutionImageMinIndex(); - - /** - * @brief SetTimeStep is used to set the time step for which the mask is to be generated - * @param timeStep - */ - void SetTimeStep(unsigned int timeStep) override; + itkSetMacro(Label, unsigned short); protected: HotspotMaskGenerator(); ~HotspotMaskGenerator() override; + Image::ConstPointer DoGetMask(unsigned int) override; + class ImageExtrema { public: bool Defined; double Max; double Min; vnl_vector MaxIndex; vnl_vector MinIndex; ImageExtrema() :Defined(false) ,Max(itk::NumericTraits::min()) ,Min(itk::NumericTraits::max()) { } }; private: /** \brief Returns size of convolution kernel depending on spacing and radius. */ template itk::Size CalculateConvolutionKernelSize(double spacing[VImageDimension], double radiusInMM); /** \brief Generates image of kernel which is needed for convolution. */ template itk::SmartPointer< itk::Image > GenerateHotspotSearchConvolutionKernel(double spacing[VImageDimension], double radiusInMM); /** \brief Convolves image with spherical kernel image. Used for hotspot calculation. */ template itk::SmartPointer< itk::Image > GenerateConvolutionImage( const itk::Image* inputImage ); /** \brief Fills pixels of the spherical hotspot mask. */ template < typename TPixel, unsigned int VImageDimension> void FillHotspotMaskPixels( itk::Image* maskImage, itk::Point sphereCenter, double sphereRadiusInMM); /** \brief */ template void - CalculateHotspotMask(itk::Image* inputImage, + CalculateHotspotMask(const itk::Image* inputImage, const itk::Image* maskImage, unsigned int label); template ImageExtrema CalculateExtremaWorld( const itk::Image* inputImage, const itk::Image* maskImage, double neccessaryDistanceToImageBorderInMM, unsigned int label); bool IsUpdateRequired() const; HotspotMaskGenerator(const HotspotMaskGenerator &); HotspotMaskGenerator & operator=(const HotspotMaskGenerator &); MaskGenerator::Pointer m_Mask; mitk::Image::Pointer m_InternalMask; - mitk::Image::Pointer m_internalImage; itk::Image::ConstPointer m_internalMask2D; itk::Image::ConstPointer m_internalMask3D; - double m_HotspotRadiusinMM; + double m_HotspotRadiusInMM; bool m_HotspotMustBeCompletelyInsideImage; unsigned short m_Label; vnl_vector m_ConvolutionImageMinIndex, m_ConvolutionImageMaxIndex; unsigned long m_InternalMaskUpdateTime; }; } #endif diff --git a/Modules/ImageStatistics/mitkIgnorePixelMaskGenerator.cpp b/Modules/ImageStatistics/mitkIgnorePixelMaskGenerator.cpp index ad1787bc84..22f94650d0 100644 --- a/Modules/ImageStatistics/mitkIgnorePixelMaskGenerator.cpp +++ b/Modules/ImageStatistics/mitkIgnorePixelMaskGenerator.cpp @@ -1,132 +1,120 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #include #include #include #include #include #include namespace mitk { void IgnorePixelMaskGenerator::SetIgnoredPixelValue(RealType pixelValue) { if (pixelValue != m_IgnoredPixelValue) { m_IgnoredPixelValue = pixelValue; this->Modified(); } } -void IgnorePixelMaskGenerator::SetTimeStep(unsigned int timeStep) +unsigned int IgnorePixelMaskGenerator::GetNumberOfMasks() const { - if (m_TimeStep != timeStep) - { - m_TimeStep = timeStep; - } + return 1; } -mitk::Image::ConstPointer IgnorePixelMaskGenerator::GetMask() +mitk::Image::ConstPointer IgnorePixelMaskGenerator::DoGetMask(unsigned int maskID) { if (IsUpdateRequired()) { - if (m_inputImage.IsNull()) + if (m_InputImage.IsNull()) { - MITK_ERROR << "Image not set!"; + mitkThrow() << "Image not set!"; } if (m_IgnoredPixelValue == std::numeric_limits::min()) { - MITK_ERROR << "IgnotePixelValue not set!"; - } - - if (m_TimeStep > (m_inputImage->GetTimeSteps() - 1)) - { - MITK_ERROR << "Invalid time step: " << m_TimeStep << ". The image has " << m_inputImage->GetTimeSteps() << " timeSteps!"; + mitkThrow() << "IgnotePixelValue not set!"; } - // extractimage time slice - ImageTimeSelector::Pointer imgTimeSel = ImageTimeSelector::New(); - imgTimeSel->SetInput(m_inputImage); - imgTimeSel->SetTimeNr(m_TimeStep); - imgTimeSel->UpdateLargestPossibleRegion(); + auto timeSliceImage = mitk::SelectImageByTimePoint(m_InputImage, m_TimePoint); - mitk::Image::Pointer timeSliceImage = imgTimeSel->GetOutput(); + if (timeSliceImage.IsNull()) mitkThrow() << "Cannot generate mask. Passed time point is not supported by input image. Invalid time point: "<< m_TimePoint; // update m_InternalMask AccessByItk(timeSliceImage, InternalCalculateMask); m_InternalMask->SetGeometry(timeSliceImage->GetGeometry()); this->Modified(); } m_InternalMaskUpdateTime = m_InternalMask->GetMTime(); return m_InternalMask; } template -void IgnorePixelMaskGenerator::InternalCalculateMask(typename itk::Image* image) +void IgnorePixelMaskGenerator::InternalCalculateMask(typename const itk::Image* image) { typedef itk::Image ImageType; typedef itk::Image MaskType; typename MaskType::Pointer mask = MaskType::New(); mask->SetOrigin(image->GetOrigin()); mask->SetSpacing(image->GetSpacing()); mask->SetLargestPossibleRegion(image->GetLargestPossibleRegion()); mask->SetBufferedRegion(image->GetBufferedRegion()); mask->SetDirection(image->GetDirection()); mask->SetNumberOfComponentsPerPixel(image->GetNumberOfComponentsPerPixel()); mask->Allocate(); mask->FillBuffer(1); // iterate over image and mask and set mask=1 if image=m_IgnorePixelValue itk::ImageRegionConstIterator imageIterator(image, image->GetLargestPossibleRegion()); itk::ImageRegionIterator maskIterator(mask, mask->GetLargestPossibleRegion()); for (imageIterator.GoToBegin(); !imageIterator.IsAtEnd(); ++imageIterator, ++maskIterator) { if (imageIterator.Value() == static_cast(m_IgnoredPixelValue)) { maskIterator.Set(0); } } m_InternalMask = GrabItkImageMemory(mask); } bool IgnorePixelMaskGenerator::IsUpdateRequired() const { unsigned long thisClassTimeStamp = this->GetMTime(); - unsigned long internalMaskTimeStamp = m_InternalMask->GetMTime(); - unsigned long inputImageTimeStamp = m_inputImage->GetMTime(); + unsigned long internalMaskTimeStamp = m_InternalMask.IsNull() ? 0 : m_InternalMask->GetMTime(); + unsigned long inputImageTimeStamp = m_InputImage->GetMTime(); if (thisClassTimeStamp > m_InternalMaskUpdateTime) // inputs have changed { return true; } if (m_InternalMaskUpdateTime < inputImageTimeStamp) // mask image has changed outside of this class { return true; } if (internalMaskTimeStamp > m_InternalMaskUpdateTime) // internal mask has been changed outside of this class { return true; } return false; } } // end namespace diff --git a/Modules/ImageStatistics/mitkIgnorePixelMaskGenerator.h b/Modules/ImageStatistics/mitkIgnorePixelMaskGenerator.h index 4903d30da8..478c443d40 100644 --- a/Modules/ImageStatistics/mitkIgnorePixelMaskGenerator.h +++ b/Modules/ImageStatistics/mitkIgnorePixelMaskGenerator.h @@ -1,84 +1,74 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #ifndef mitkIgnorePixelMaskGenerator_h #define mitkIgnorePixelMaskGenerator_h #include #include #include #include #include namespace mitk { /** * @brief The IgnorePixelMaskGenerator class is used to generate a mask that is zero for specific pixel values in the input image. This class requires an input image. */ class MITKIMAGESTATISTICS_EXPORT IgnorePixelMaskGenerator: public MaskGenerator { public: /** Standard Self typedef */ typedef IgnorePixelMaskGenerator Self; typedef MaskGenerator Superclass; typedef itk::SmartPointer< Self > Pointer; typedef itk::SmartPointer< const Self > ConstPointer; typedef double RealType; /** Method for creation through the object factory. */ itkNewMacro(Self); /** Runtime information support. */ itkTypeMacro(IgnorePixelMaskGenerator, MaskGenerator); /** * @brief The mask will be 0 there inputImage==pixelValue and 1 otherwise */ void SetIgnoredPixelValue(RealType pixelValue); - /** - * @brief Computes and returns the mask - */ - mitk::Image::ConstPointer GetMask() override; - - /** - * @brief SetTimeStep is used to set the time step for which the mask is to be generated - * @param timeStep - */ - void SetTimeStep(unsigned int timeStep) override; + unsigned int GetNumberOfMasks() const override; protected: IgnorePixelMaskGenerator(): - m_IgnoredPixelValue(std::numeric_limits::min()) + m_IgnoredPixelValue(std::numeric_limits::min()), m_InternalMaskUpdateTime(0) { - m_TimeStep = 0; - m_InternalMaskUpdateTime = 0; - m_InternalMask = mitk::Image::New(); } - ~IgnorePixelMaskGenerator() override{} + ~IgnorePixelMaskGenerator() = default; + + mitk::Image::ConstPointer DoGetMask(unsigned int maskID) override; template - void InternalCalculateMask(typename itk::Image* image); + void InternalCalculateMask(typename const itk::Image* image); private: bool IsUpdateRequired() const; mitk::Image::Pointer m_InternalMask; RealType m_IgnoredPixelValue; unsigned long m_InternalMaskUpdateTime; }; } #endif diff --git a/Modules/ImageStatistics/mitkImageMaskGenerator.cpp b/Modules/ImageStatistics/mitkImageMaskGenerator.cpp index abf894aca5..eb0cdcc572 100644 --- a/Modules/ImageStatistics/mitkImageMaskGenerator.cpp +++ b/Modules/ImageStatistics/mitkImageMaskGenerator.cpp @@ -1,99 +1,84 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #include #include #include namespace mitk { -void ImageMaskGenerator::SetImageMask(const Image* maskImage) -{ - if (m_InternalMaskImage != maskImage) - { - m_InternalMaskImage = maskImage; - this->Modified(); - } -} - -void ImageMaskGenerator::SetTimeStep(unsigned int timeStep) -{ - if (timeStep != m_TimeStep) - { - m_TimeStep = timeStep; - UpdateInternalMask(); - } -} - void ImageMaskGenerator::UpdateInternalMask() { - if (this->m_inputImage.IsNull()) + if (this->m_ImageMask.IsNull()) { mitkThrow() << "Cannot update internal mask. Input image is not set."; } - const auto timeGeo = this->m_inputImage->GetTimeGeometry(); - if (!timeGeo->IsValidTimeStep(this->m_TimeStep)) + const auto timeGeo = this->m_InputImage->GetTimeGeometry(); + if (!timeGeo->IsValidTimePoint(this->m_TimePoint)) { mitkThrow() << "Cannot update internal mask. Time step selected that is not supported by input image."; } - auto timePoint = this->m_inputImage->GetTimeGeometry()->TimeStepToTimePoint(this->m_TimeStep); - m_InternalMask = SelectImageByTimePoint(m_InternalMaskImage, timePoint); + m_InternalMask = SelectImageByTimePoint(m_ImageMask, m_TimePoint); if (m_InternalMask.IsNull()) { MITK_WARN << "Warning: time step > number of time steps in mask image, using last time step"; - m_InternalMask = SelectImageByTimeStep(m_InternalMaskImage, m_InternalMaskImage->GetTimeSteps()-1); + m_InternalMask = SelectImageByTimeStep(m_ImageMask, m_ImageMask->GetTimeSteps()-1); } } +unsigned int ImageMaskGenerator::GetNumberOfMasks() const +{ + return 1; +} -mitk::Image::ConstPointer ImageMaskGenerator::GetMask() +mitk::Image::ConstPointer ImageMaskGenerator::DoGetMask(unsigned int) { - if (m_InternalMaskImage.IsNull()) + if (m_ImageMask.IsNull()) { - MITK_ERROR << "Mask Image is nullptr"; + mitkThrow() << "Mask Image is nullptr"; } if (IsUpdateRequired()) { UpdateInternalMask(); } return m_InternalMask; } bool ImageMaskGenerator::IsUpdateRequired() const { - unsigned long internalMaskTimeStamp = m_InternalMask->GetMTime(); - unsigned long maskImageTimeStamp = m_InternalMaskImage->GetMTime(); + unsigned long internalMaskTimeStamp = m_InternalMask.IsNull() ? 0 : m_InternalMask->GetMTime(); + unsigned long maskImageTimeStamp = m_ImageMask->GetMTime(); if (maskImageTimeStamp > internalMaskTimeStamp) // inputs have changed { return true; } if (this->GetMTime() > maskImageTimeStamp) // input has changed { return true; } return false; } } diff --git a/Modules/ImageStatistics/mitkImageMaskGenerator.h b/Modules/ImageStatistics/mitkImageMaskGenerator.h index e475009fd0..e97bc76855 100644 --- a/Modules/ImageStatistics/mitkImageMaskGenerator.h +++ b/Modules/ImageStatistics/mitkImageMaskGenerator.h @@ -1,62 +1,61 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #ifndef mitkImageMaskGenerator_h #define mitkImageMaskGenerator_h #include #include #include #include #include namespace mitk { class MITKIMAGESTATISTICS_EXPORT ImageMaskGenerator: public MaskGenerator { public: /** Standard Self typedef */ typedef ImageMaskGenerator Self; typedef MaskGenerator 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(BinaryImageMaskGenerator, MaskGenerator); - mitk::Image::ConstPointer GetMask() override; + unsigned int GetNumberOfMasks() const override; - void SetTimeStep(unsigned int timeStep) override; - - void SetImageMask(const mitk::Image* maskImage); + itkSetConstObjectMacro(ImageMask, Image) protected: ImageMaskGenerator():Superclass(){ m_InternalMaskUpdateTime = 0; - m_InternalMask = mitk::Image::New(); } + Image::ConstPointer DoGetMask(unsigned int) override; + private: bool IsUpdateRequired() const; void UpdateInternalMask(); - mitk::Image::ConstPointer m_InternalMaskImage; + mitk::Image::ConstPointer m_ImageMask; mitk::Image::ConstPointer m_InternalMask; unsigned long m_InternalMaskUpdateTime; }; } #endif diff --git a/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp b/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp index 25031da9b4..7b61e7ca7e 100644 --- a/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp +++ b/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp @@ -1,510 +1,524 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #include "mitkImageStatisticsCalculator.h" #include #include #include #include #include #include #include #include #include #include #include #include namespace mitk { void ImageStatisticsCalculator::SetInputImage(const mitk::Image *image) { if (image != m_Image) { m_Image = image; this->Modified(); } } void ImageStatisticsCalculator::SetMask(mitk::MaskGenerator *mask) { if (mask != m_MaskGenerator) { m_MaskGenerator = mask; this->Modified(); } } void ImageStatisticsCalculator::SetSecondaryMask(mitk::MaskGenerator *mask) { 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::SetBinSizeForHistogramStatistics(double binSize) { if (binSize != m_binSizeForHistogramStatistics) { m_binSizeForHistogramStatistics = binSize; this->Modified(); this->m_UseBinSizeOverNBins = true; } if (!m_UseBinSizeOverNBins) { this->Modified(); this->m_UseBinSizeOverNBins = true; } } double ImageStatisticsCalculator::GetBinSizeForHistogramStatistics() const { return m_binSizeForHistogramStatistics; } mitk::ImageStatisticsContainer* ImageStatisticsCalculator::GetStatistics() { if (m_Image.IsNull()) { mitkThrow() << "no image"; } if (!m_Image->IsInitialized()) { mitkThrow() << "Image not initialized!"; } if (IsUpdateRequired()) { auto timeGeometry = m_Image->GetTimeGeometry(); m_StatisticContainer = ImageStatisticsContainer::New(); m_StatisticContainer->SetTimeGeometry(timeGeometry->Clone()); // always compute statistics on all timesteps - for (unsigned int timeStep = 0; timeStep < m_Image->GetTimeSteps(); timeStep++) + for (TimeStepType timeStep = 0; timeStep < m_Image->GetTimeSteps(); timeStep++) { - if (m_MaskGenerator.IsNotNull()) + const unsigned int numbersOfMasks = m_MaskGenerator.IsNotNull() ? m_MaskGenerator->GetNumberOfMasks() : 1; + auto timePoint = m_Image->GetTimeGeometry()->TimeStepToTimePoint(timeStep); + //hard coded 0 is + //a workaround , open up a task refer to it that it should be removed + // as soon as the secondary mask is removed anyways and solved by a generator chain. + if (m_SecondaryMaskGenerator.IsNotNull()) + { + if (m_SecondaryMaskGenerator->GetNumberOfMasks() != 1) + mitkThrow() << "Cannot generate secondary mask. ImageStatisticsCalculator does only support secondary mask generators with on mask. Number of masks provided: " + << m_SecondaryMaskGenerator->GetNumberOfMasks(); + m_SecondaryMaskGenerator->SetTimePoint(timePoint); + m_SecondaryMask = m_SecondaryMaskGenerator->GetMask(0); + } + + for (unsigned int maskID = 0; maskID < numbersOfMasks; ++maskID) { - m_MaskGenerator->SetTimeStep(timeStep); - m_InternalMask = m_MaskGenerator->GetMask(); - if (m_MaskGenerator->GetReferenceImage().IsNotNull()) + if (m_MaskGenerator.IsNotNull()) { - m_InternalImageForStatistics = m_MaskGenerator->GetReferenceImage(); + m_MaskGenerator->SetTimePoint(timePoint); + m_InternalMask = m_MaskGenerator->GetMask(maskID); + if (m_MaskGenerator->GetReferenceImage().IsNotNull()) + { + m_InternalImageForStatistics = m_MaskGenerator->GetReferenceImage(); + } + else + { + m_InternalImageForStatistics = m_Image; + } } 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(); - imgTimeSel->Update(); - m_ImageTimeSlice = imgTimeSel->GetOutput(); + m_ImageTimeSlice = SelectImageByTimeStep(m_InternalImageForStatistics, timeStep); - // 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) + // 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) + } } } } return m_StatisticContainer; } template void ImageStatisticsCalculator::InternalCalculateStatisticsUnmasked( - typename itk::Image *image, const TimeGeometry *timeGeometry, TimeStepType timeStep) + typename const itk::Image *image, const TimeGeometry *timeGeometry, TimeStepType timeStep) { typedef typename itk::Image ImageType; typedef typename mitk::StatisticsImageFilter ImageStatisticsFilterType; typedef typename itk::MinMaxImageFilterWithIndex MinMaxFilterType; auto statObj = ImageStatisticsContainer::ImageStatisticsObject(); 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; 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 = imgMinMaxFilter->GetIndexOfMinimum(); // typename ImageType::IndexType tmpMaxIndex = imgMinMaxFilter->GetIndexOfMaximum(); 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]; } 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; } 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(); } 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(ImageStatisticsConstants::NUMBEROFVOXELS(), static_cast(numberOfPixels)); statObj.AddStatistic(ImageStatisticsConstants::VOLUME(), volume); statObj.AddStatistic(ImageStatisticsConstants::MEAN(), statisticsFilter->GetMean()); statObj.AddStatistic(ImageStatisticsConstants::MINIMUM(), static_cast(statisticsFilter->GetMinimum())); statObj.AddStatistic(ImageStatisticsConstants::MAXIMUM(), static_cast(statisticsFilter->GetMaximum())); statObj.AddStatistic(ImageStatisticsConstants::STANDARDDEVIATION(), statisticsFilter->GetSigma()); statObj.AddStatistic(ImageStatisticsConstants::VARIANCE(), variance); statObj.AddStatistic(ImageStatisticsConstants::SKEWNESS(), statisticsFilter->GetSkewness()); statObj.AddStatistic(ImageStatisticsConstants::KURTOSIS(), statisticsFilter->GetKurtosis()); statObj.AddStatistic(ImageStatisticsConstants::RMS(), rms); statObj.AddStatistic(ImageStatisticsConstants::MPP(), statisticsFilter->GetMPP()); statObj.AddStatistic(ImageStatisticsConstants::ENTROPY(), statisticsFilter->GetEntropy()); statObj.AddStatistic(ImageStatisticsConstants::MEDIAN(), statisticsFilter->GetMedian()); statObj.AddStatistic(ImageStatisticsConstants::UNIFORMITY(), statisticsFilter->GetUniformity()); statObj.AddStatistic(ImageStatisticsConstants::UPP(), statisticsFilter->GetUPP()); statObj.m_Histogram = statisticsFilter->GetHistogram(); m_StatisticContainer->SetStatistics(ImageStatisticsContainer::NO_MASK_LABEL_VALUE, timeStep, statObj); } template - double ImageStatisticsCalculator::GetVoxelVolume(typename itk::Image *image) const + double ImageStatisticsCalculator::GetVoxelVolume(typename const 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, + void ImageStatisticsCalculator::InternalCalculateStatisticsMasked(typename const itk::Image *image, const TimeGeometry *timeGeometry, unsigned int timeStep) { typedef itk::Image ImageType; typedef itk::Image MaskType; typedef typename MaskType::PixelType LabelPixelType; typedef LabelStatisticsImageFilter ImageStatisticsFilterType; typedef MaskUtilities MaskUtilType; typedef typename itk::MinMaxLabelImageFilterWithIndex MinMaxLabelFilterType; - // 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) + // workaround: if m_SecondaryMaskGenerator is not null but m_MaskGenerator is! (this is the case if we request a + // 'ignore zero 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; } // maskImage has to have the same dimension as image typename MaskType::ConstPointer 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 &) - { - typename MaskType::Pointer noneConstMaskImage; //needed to work arround the fact that CastToItkImage currently does not support const itk images. + typename MaskType::Pointer noneConstMaskImage; //needed to work around the fact that CastToItkImage currently does not support const itk images. // 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, noneConstMaskImage); maskImage = noneConstMaskImage; } // 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)) { Image::ConstPointer old_img = m_SecondaryMaskGenerator->GetReferenceImage(); m_SecondaryMaskGenerator->SetInputImage(m_MaskGenerator->GetReferenceImage()); - m_SecondaryMask = m_SecondaryMaskGenerator->GetMask(); + m_SecondaryMask = m_SecondaryMaskGenerator->GetMask(0); m_SecondaryMaskGenerator->SetInputImage(old_img); } typename MaskType::ConstPointer 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::ConstPointer 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(); } typename MaskUtilType::Pointer maskUtil = MaskUtilType::New(); maskUtil->SetImage(image); maskUtil->SetMask(maskImage.GetPointer()); // if mask is smaller than image, extract the image region where the mask is typename ImageType::ConstPointer adaptedImage = ImageType::New(); adaptedImage = maskUtil->ExtractMaskImageRegion(); // this also checks mask sanity // find min, max, minindex and maxindex typename MinMaxLabelFilterType::Pointer minMaxFilter = MinMaxLabelFilterType::New(); minMaxFilter->SetInput(adaptedImage); minMaxFilter->SetLabelInput(maskImage); minMaxFilter->UpdateLargestPossibleRegion(); // set histogram parameters for each label individually (min/max may be different for each label) typedef typename std::unordered_map MapType; std::vector relevantLabels = minMaxFilter->GetRelevantLabels(); MapType minVals; MapType maxVals; std::unordered_map nBins; for (LabelPixelType label : relevantLabels) { minVals[label] = static_cast(minMaxFilter->GetMin(label)); maxVals[label] = static_cast(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[label] = nBinsForHistogram; } typename ImageStatisticsFilterType::Pointer imageStatisticsFilter = ImageStatisticsFilterType::New(); imageStatisticsFilter->SetDirectionTolerance(0.001); imageStatisticsFilter->SetCoordinateTolerance(0.001); imageStatisticsFilter->SetInput(adaptedImage); imageStatisticsFilter->SetLabelInput(maskImage); imageStatisticsFilter->SetHistogramParameters(nBins, minVals, maxVals); imageStatisticsFilter->Update(); const auto labels = imageStatisticsFilter->GetValidLabelValues(); for (auto labelValue : labels) { + if (labelValue == ImageStatisticsContainer::NO_MASK_LABEL_VALUE) + { + //we ignore the background of a mask if we compute mask statistics. + continue; + } + ImageStatisticsContainer::ImageStatisticsObject 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; Point3D worldCoordinateMin; Point3D worldCoordinateMax; Point3D indexCoordinateMin; Point3D indexCoordinateMax; m_InternalImageForStatistics->GetGeometry()->IndexToWorld(minMaxFilter->GetMinIndex(labelValue), worldCoordinateMin); m_InternalImageForStatistics->GetGeometry()->IndexToWorld(minMaxFilter->GetMaxIndex(labelValue), 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(ImageStatisticsConstants::MINIMUMPOSITION(), minIndex); statObj.AddStatistic(ImageStatisticsConstants::MAXIMUMPOSITION(), maxIndex); auto voxelVolume = GetVoxelVolume(image); auto numberOfVoxels = static_cast(imageStatisticsFilter->GetCount(labelValue)); auto volume = static_cast(numberOfVoxels) * voxelVolume; auto rms = std::sqrt(std::pow(imageStatisticsFilter->GetMean(labelValue), 2.) + imageStatisticsFilter->GetVariance(labelValue)); // variance = sigma^2 auto variance = imageStatisticsFilter->GetSigma(labelValue) * imageStatisticsFilter->GetSigma(labelValue); statObj.AddStatistic(ImageStatisticsConstants::NUMBEROFVOXELS(), numberOfVoxels); statObj.AddStatistic(ImageStatisticsConstants::VOLUME(), volume); statObj.AddStatistic(ImageStatisticsConstants::MEAN(), imageStatisticsFilter->GetMean(labelValue)); statObj.AddStatistic(ImageStatisticsConstants::MINIMUM(), static_cast(imageStatisticsFilter->GetMinimum(labelValue))); statObj.AddStatistic(ImageStatisticsConstants::MAXIMUM(), static_cast(imageStatisticsFilter->GetMaximum(labelValue))); statObj.AddStatistic(ImageStatisticsConstants::STANDARDDEVIATION(), imageStatisticsFilter->GetSigma(labelValue)); statObj.AddStatistic(ImageStatisticsConstants::VARIANCE(), variance); statObj.AddStatistic(ImageStatisticsConstants::SKEWNESS(), imageStatisticsFilter->GetSkewness(labelValue)); statObj.AddStatistic(ImageStatisticsConstants::KURTOSIS(), imageStatisticsFilter->GetKurtosis(labelValue)); statObj.AddStatistic(ImageStatisticsConstants::RMS(), rms); statObj.AddStatistic(ImageStatisticsConstants::MPP(), imageStatisticsFilter->GetMPP(labelValue)); statObj.AddStatistic(ImageStatisticsConstants::ENTROPY(), imageStatisticsFilter->GetEntropy(labelValue)); statObj.AddStatistic(ImageStatisticsConstants::MEDIAN(), imageStatisticsFilter->GetMedian(labelValue)); statObj.AddStatistic(ImageStatisticsConstants::UNIFORMITY(), imageStatisticsFilter->GetUniformity(labelValue)); statObj.AddStatistic(ImageStatisticsConstants::UPP(), imageStatisticsFilter->GetUPP(labelValue)); statObj.m_Histogram = imageStatisticsFilter->GetHistogram(labelValue); + if (m_StatisticContainer->StatisticsExist(labelValue, timeStep)) + mitkThrow() << "Invalid state/input data. Statistic for a specific label/time step pair was computed more then once. Conflicting label ID: " + << labelValue << " ; conflicting time step: " << timeStep; m_StatisticContainer->SetStatistics(labelValue, timeStep, statObj); } // swap maskGenerators back if (swapMasks) { m_SecondaryMask = m_InternalMask; m_InternalMask = nullptr; } } bool ImageStatisticsCalculator::IsUpdateRequired() const { const auto thisClassTimeStamp = this->GetMTime(); const auto inputImageTimeStamp = m_Image->GetMTime(); if (m_StatisticContainer.IsNull()) { return true; } const auto statisticsTimeStamp = m_StatisticContainer->GetMTime(); if (thisClassTimeStamp > statisticsTimeStamp) // inputs have changed { return true; } if (inputImageTimeStamp > statisticsTimeStamp) // image has changed { return true; } if (m_MaskGenerator.IsNotNull()) { const auto maskGeneratorTimeStamp = m_MaskGenerator->GetMTime(); if (maskGeneratorTimeStamp > statisticsTimeStamp) // there is a mask generator and it has changed { return true; } } if (m_SecondaryMaskGenerator.IsNotNull()) { const auto 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 d8f9c5c7d4..89c55746df 100644 --- a/Modules/ImageStatistics/mitkImageStatisticsCalculator.h +++ b/Modules/ImageStatistics/mitkImageStatisticsCalculator.h @@ -1,118 +1,118 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #ifndef mitkImageStatisticsCalculator_h #define mitkImageStatisticsCalculator_h #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::LabelValueType; /**Documentation @brief Set the image for which the statistics are to be computed.*/ void SetInputImage(const mitk::Image* 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* 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* 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. If these requested statistics are not computed yet the computation is done as well. */ ImageStatisticsContainer* GetStatistics(); 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, const TimeGeometry* timeGeometry, TimeStepType timeStep); + typename const itk::Image< TPixel, VImageDimension >* image, const TimeGeometry* timeGeometry, TimeStepType timeStep); template < typename TPixel, unsigned int VImageDimension > void InternalCalculateStatisticsMasked( - typename itk::Image< TPixel, VImageDimension >* image, const TimeGeometry* timeGeometry, + typename const itk::Image< TPixel, VImageDimension >* image, const TimeGeometry* timeGeometry, unsigned int timeStep); template < typename TPixel, unsigned int VImageDimension > - double GetVoxelVolume(typename itk::Image* image) const; + double GetVoxelVolume(typename const itk::Image* image) const; bool IsUpdateRequired() const; mitk::Image::ConstPointer m_Image; - mitk::Image::Pointer m_ImageTimeSlice; + mitk::Image::ConstPointer m_ImageTimeSlice; mitk::Image::ConstPointer m_InternalImageForStatistics; mitk::MaskGenerator::Pointer m_MaskGenerator; mitk::Image::ConstPointer m_InternalMask; mitk::MaskGenerator::Pointer m_SecondaryMaskGenerator; mitk::Image::ConstPointer m_SecondaryMask; unsigned int m_nBinsForHistogramStatistics; double m_binSizeForHistogramStatistics; bool m_UseBinSizeOverNBins; ImageStatisticsContainer::Pointer m_StatisticContainer; }; } #endif diff --git a/Modules/ImageStatistics/mitkImageStatisticsContainer.cpp b/Modules/ImageStatistics/mitkImageStatisticsContainer.cpp index 90213a565c..17d91eadf7 100644 --- a/Modules/ImageStatistics/mitkImageStatisticsContainer.cpp +++ b/Modules/ImageStatistics/mitkImageStatisticsContainer.cpp @@ -1,306 +1,303 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #include #include #include #include namespace mitk { ImageStatisticsContainer::ImageStatisticsContainer() { this->Reset(); } // The order is derived from the old (<2018) image statistics plugin. const ImageStatisticsContainer::ImageStatisticsObject::StatisticNameVector ImageStatisticsContainer::ImageStatisticsObject::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::ImageStatisticsObject::ImageStatisticsObject() { Reset(); } void ImageStatisticsContainer::ImageStatisticsObject::AddStatistic(const std::string_view key, StatisticsVariantType value) { m_Statistics.emplace(key, value); if (std::find(m_DefaultNames.cbegin(), m_DefaultNames.cend(), key) == m_DefaultNames.cend()) { if (std::find(m_CustomNames.cbegin(), m_CustomNames.cend(), key) == m_CustomNames.cend()) { m_CustomNames.emplace_back(key); } } } const ImageStatisticsContainer::ImageStatisticsObject::StatisticNameVector & ImageStatisticsContainer::ImageStatisticsObject::GetDefaultStatisticNames() { return m_DefaultNames; } const ImageStatisticsContainer::ImageStatisticsObject::StatisticNameVector & ImageStatisticsContainer::ImageStatisticsObject::GetCustomStatisticNames() const { return m_CustomNames; } ImageStatisticsContainer::ImageStatisticsObject::StatisticNameVector ImageStatisticsContainer::ImageStatisticsObject::GetAllStatisticNames() const { StatisticNameVector names = GetDefaultStatisticNames(); names.insert(names.cend(), m_CustomNames.cbegin(), m_CustomNames.cend()); return names; } ImageStatisticsContainer::ImageStatisticsObject::StatisticNameVector ImageStatisticsContainer::ImageStatisticsObject::GetExistingStatisticNames() const { StatisticNameVector names; std::transform(m_Statistics.begin(), m_Statistics.end(), std::back_inserter(names), [](const auto &pair) { return pair.first; }); return names; } bool ImageStatisticsContainer::ImageStatisticsObject::HasStatistic(const std::string_view name) const { return m_Statistics.find(name) != m_Statistics.cend(); } ImageStatisticsContainer::StatisticsVariantType ImageStatisticsContainer::ImageStatisticsObject::GetValueNonConverted( const std::string_view name) const { if (HasStatistic(name)) { return m_Statistics.find(name)->second; } else { mitkThrow() << "invalid statistic key, could not find"; } } void ImageStatisticsContainer::ImageStatisticsObject::Reset() { m_Statistics.clear(); m_CustomNames.clear(); } const static ImageStatisticsContainer::LabelValueType NO_MASK_LABEL_VALUE = 1; bool ImageStatisticsContainer::StatisticsExist(LabelValueType labelValue, TimeStepType timeStep) const { auto labelFinding = m_LabelTimeStep2StatisticsMap.find(labelValue); if (labelFinding == m_LabelTimeStep2StatisticsMap.end()) return false; auto timeFinding = labelFinding->second.find(timeStep); return timeFinding != labelFinding->second.end(); } const ImageStatisticsContainer::HistogramType* ImageStatisticsContainer::GetHistogram(LabelValueType labelValue, TimeStepType timeStep) const { return this->GetStatistics(labelValue, timeStep).m_Histogram; } bool ImageStatisticsContainer::IgnoresZeroVoxel() const { const auto prop = this->GetProperty(mitk::STATS_IGNORE_ZERO_VOXEL_PROPERTY_NAME.c_str()); const auto boolProp = dynamic_cast(prop.GetPointer()); if (nullptr != boolProp) { return boolProp->GetValue(); } return false; } bool ImageStatisticsContainer::IsWIP() const { return this->GetProperty(mitk::STATS_GENERATION_STATUS_PROPERTY_NAME.c_str()).IsNotNull(); } const ImageStatisticsContainer::ImageStatisticsObject &ImageStatisticsContainer::GetStatistics(LabelValueType labelValue, TimeStepType timeStep) const { auto labelFinding = m_LabelTimeStep2StatisticsMap.find(labelValue); if (labelFinding == m_LabelTimeStep2StatisticsMap.end()) mitkThrow() << "Cannot get statistics. Requested label value does not exist. Invalid label:" <second.find(timeStep); if (timeFinding == labelFinding->second.end()) mitkThrow() << "Cannot get statistics. Requested time step does not exist. Invalid time step:" << timeStep; return timeFinding->second; } void ImageStatisticsContainer::SetStatistics(LabelValueType labelValue, TimeStepType timeStep, const ImageStatisticsObject& statistics) { if (!this->GetTimeGeometry()->IsValidTimeStep(timeStep)) mitkThrow() << "Given timeStep " << timeStep << " out of TimeGeometry bounds of the object. TimeSteps in geometry: " << this->GetTimeSteps(); m_LabelTimeStep2StatisticsMap[labelValue][timeStep] = statistics; this->Modified(); } void ImageStatisticsContainer::PrintSelf(std::ostream &os, itk::Indent indent) const { Superclass::PrintSelf(os, indent); for (const auto& [labelValue, timeMap] : m_LabelTimeStep2StatisticsMap) { for (const auto& [timeStep, container] : timeMap) { os << std::endl << indent << "Statistics instance (Label "<< labelValue << ", TimeStep "<< timeStep << "):"; auto statisticsValues = GetStatistics(labelValue,timeStep); auto statisticKeys = statisticsValues.GetExistingStatisticNames(); os << std::endl << indent << "Number of entries: " << statisticKeys.size(); for (const auto& aKey : statisticKeys) { os << std::endl << indent.GetNextIndent() << aKey << ": " << statisticsValues.GetValueNonConverted(aKey); } } } } - ImageStatisticsContainer::LabelValueVectorType ImageStatisticsContainer::GetExistingLabelValues(bool ignoreUnlabeled) const + ImageStatisticsContainer::LabelValueVectorType ImageStatisticsContainer::GetExistingLabelValues() const { LabelValueVectorType result; for (const auto& [labelValue, timeMap] : m_LabelTimeStep2StatisticsMap) { - if (!ignoreUnlabeled || labelValue != NO_MASK_LABEL_VALUE) - { - result.push_back(labelValue); - } + result.push_back(labelValue); } return result; } ImageStatisticsContainer::TimeStepVectorType ImageStatisticsContainer::GetExistingTimeSteps(LabelValueType labelValue) const { TimeStepVectorType result; auto labelFinding = m_LabelTimeStep2StatisticsMap.find(labelValue); if (labelFinding == m_LabelTimeStep2StatisticsMap.end()) mitkThrow() << "Cannot get existing time steps. Requested label value does not exist. Invalid label:" << labelValue; for (const auto& [timestep, stats] : labelFinding->second) { result.push_back(timestep); } return result; } void ImageStatisticsContainer::Reset() { m_LabelTimeStep2StatisticsMap.clear(); this->Modified(); } 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->m_LabelTimeStep2StatisticsMap = m_LabelTimeStep2StatisticsMap; rval->SetTimeGeometry(this->GetTimeGeometry()->Clone()); return ioPtr; } ImageStatisticsContainer::ImageStatisticsObject::StatisticNameVector GetAllStatisticNames( const ImageStatisticsContainer *container) { ImageStatisticsContainer::ImageStatisticsObject::StatisticNameVector names = ImageStatisticsContainer::ImageStatisticsObject::GetDefaultStatisticNames(); if (container) { std::set customKeys; - auto labelValues = container->GetExistingLabelValues(false); + auto labelValues = container->GetExistingLabelValues(); for (const auto labelValue : labelValues) { auto timeSteps = container->GetExistingTimeSteps(labelValue); for (const auto timeStep : timeSteps) { auto statisticKeys = container->GetStatistics(labelValue, timeStep).GetCustomStatisticNames(); customKeys.insert(statisticKeys.cbegin(), statisticKeys.cend()); } } names.insert(names.cend(), customKeys.cbegin(), customKeys.cend()); } return names; } ImageStatisticsContainer::ImageStatisticsObject::StatisticNameVector GetAllStatisticNames( std::vector containers) { ImageStatisticsContainer::ImageStatisticsObject::StatisticNameVector names = ImageStatisticsContainer::ImageStatisticsObject::GetDefaultStatisticNames(); std::set customKeys; for (const auto &container : containers) { std::set customKeys; - auto labelValues = container->GetExistingLabelValues(false); + auto labelValues = container->GetExistingLabelValues(); for (const auto labelValue : labelValues) { auto timeSteps = container->GetExistingTimeSteps(labelValue); for (const auto timeStep : timeSteps) { auto statisticKeys = container->GetStatistics(labelValue, timeStep).GetCustomStatisticNames(); customKeys.insert(statisticKeys.cbegin(), statisticKeys.cend()); } } names.insert(names.cend(), customKeys.cbegin(), customKeys.cend()); } names.insert(names.end(), customKeys.begin(), customKeys.end()); return names; }; } // namespace mitk diff --git a/Modules/ImageStatistics/mitkImageStatisticsContainer.h b/Modules/ImageStatistics/mitkImageStatisticsContainer.h index db315f8460..5d6bc73292 100644 --- a/Modules/ImageStatistics/mitkImageStatisticsContainer.h +++ b/Modules/ImageStatistics/mitkImageStatisticsContainer.h @@ -1,182 +1,182 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #ifndef mitkImageStatisticsContainer_h #define mitkImageStatisticsContainer_h #include #include #include #include #include #include namespace mitk { /** @brief Container class for storing a StatisticsObject for each time step. Stored statistics are: - for the defined statistics, see GetAllStatisticNames - 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 LabelValueType = LabelSetImage::LabelValueType; void SetRequestedRegionToLargestPossibleRegion() override {} bool RequestedRegionIsOutsideOfTheBufferedRegion() override { return false; } bool VerifyRequestedRegion() override { return true; } 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 later on. */ class MITKIMAGESTATISTICS_EXPORT ImageStatisticsObject { public: ImageStatisticsObject(); using RealType = double; using IndexType = vnl_vector; using VoxelCountType = unsigned long; using StatisticsVariantType = boost::variant; /** @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_view 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; StatisticNameVector GetExistingStatisticNames() const; bool HasStatistic(const std::string_view 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_view 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_view name) const; void Reset(); HistogramType::ConstPointer m_Histogram=nullptr; private: using StatisticsMapType = std::map < std::string, StatisticsVariantType, std::less<>>; StatisticsMapType m_Statistics; StatisticNameVector m_CustomNames; static const StatisticNameVector m_DefaultNames; }; using StatisticsVariantType = ImageStatisticsObject::StatisticsVariantType; using RealType = ImageStatisticsObject::RealType; using IndexType = ImageStatisticsObject::IndexType; using VoxelCountType = ImageStatisticsObject::VoxelCountType; using TimeStepVectorType = std::vector; TimeStepVectorType GetExistingTimeSteps(LabelValueType labelValue) const; /** Value that can be used to query for the statistic if no mask was provided.*/ const static LabelValueType NO_MASK_LABEL_VALUE = Label::UNLABELED_VALUE; using LabelValueVectorType = LabelSetImage::LabelValueVectorType; - LabelValueVectorType GetExistingLabelValues(bool ignoreUnlabeled) const; + LabelValueVectorType GetExistingLabelValues() const; /** @brief Deletes all stored values*/ void Reset(); const ImageStatisticsObject& GetStatistics(LabelValueType labelValue, TimeStepType timeStep) const; /** @brief Sets the statisticObject for the given Timestep @pre timeStep must be valid */ void SetStatistics(LabelValueType labelValue, TimeStepType timeStep, const ImageStatisticsObject& statistics); /** @brief Checks if the Time step exists @pre timeStep must be valid */ bool StatisticsExist(LabelValueType labelValue, TimeStepType timeStep) const; /** /brief Returns the histogram of the passed time step. @pre timeStep must be valid*/ const HistogramType* GetHistogram(LabelValueType labelValue, TimeStepType timeStep) const; bool IgnoresZeroVoxel() const; bool IsWIP() const; protected: ImageStatisticsContainer(); void PrintSelf(std::ostream &os, itk::Indent indent) const override; private: itk::LightObject::Pointer InternalClone() const override; using TimeStepMapType = std::map; using LabelMapType = std::map; LabelMapType m_LabelTimeStep2StatisticsMap; }; MITKIMAGESTATISTICS_EXPORT ImageStatisticsContainer::ImageStatisticsObject::StatisticNameVector GetAllStatisticNames(const ImageStatisticsContainer* container); MITKIMAGESTATISTICS_EXPORT ImageStatisticsContainer::ImageStatisticsObject::StatisticNameVector GetAllStatisticNames(std::vector containers); } #endif diff --git a/Modules/ImageStatistics/mitkMaskGenerator.cpp b/Modules/ImageStatistics/mitkMaskGenerator.cpp index 566045014d..121c6469ec 100644 --- a/Modules/ImageStatistics/mitkMaskGenerator.cpp +++ b/Modules/ImageStatistics/mitkMaskGenerator.cpp @@ -1,46 +1,34 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #include namespace mitk { MaskGenerator::MaskGenerator(): - m_TimeStep(0) + m_TimePoint(0) { - m_inputImage = nullptr; } -void MaskGenerator::SetTimeStep(unsigned int timeStep) +mitk::Image::ConstPointer MaskGenerator::GetMask(unsigned int maskID) { - if (timeStep != m_TimeStep) - { - m_TimeStep = timeStep; - this->Modified(); - } -} + if (maskID >= this->GetNumberOfMasks()) mitkThrow() << "Cannot generate and return mask. Passed mask ID is invalid. Invalid ID: " << maskID; -void MaskGenerator::SetInputImage(mitk::Image::ConstPointer inputImg) -{ - if (inputImg != m_inputImage) - { - m_inputImage = inputImg; - this->Modified(); - } + return this->DoGetMask(maskID); } mitk::Image::ConstPointer MaskGenerator::GetReferenceImage() { - return m_inputImage; + return m_InputImage; } } diff --git a/Modules/ImageStatistics/mitkMaskGenerator.h b/Modules/ImageStatistics/mitkMaskGenerator.h index bcca054885..2e5cf11c1c 100644 --- a/Modules/ImageStatistics/mitkMaskGenerator.h +++ b/Modules/ImageStatistics/mitkMaskGenerator.h @@ -1,73 +1,73 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #ifndef mitkMaskGenerator_h #define mitkMaskGenerator_h #include #include #include #include #include namespace mitk { /** * \class MaskGenerator * \brief Base Class for all Mask Generators. Mask generators are classes that provide functionality for the * creation of binary (or unsigned short) masks that can be applied to an image. See dervied classes for more * information. */ class MITKIMAGESTATISTICS_EXPORT MaskGenerator: public itk::Object { public: - /** Standard Self typedef */ - typedef MaskGenerator Self; - typedef itk::Object Superclass; - typedef itk::SmartPointer< Self > Pointer; - typedef itk::SmartPointer< const Self > ConstPointer; - - /** Method for creation through the object factory. */ - itkTypeMacro(MaskGenerator, itk::Object); - - //~MaskGenerator(); + mitkClassMacroItkParent(MaskGenerator, itk::Object); + virtual unsigned int GetNumberOfMasks() const = 0; /** - * @brief GetMask must be overridden by derived classes. - * @return mitk::Image::Pointer of generated mask + * @brief GetMask returns the requested (multi) label mask. + * @param maskID Parameter indicating which mask should be returned. + * @return mitk::Image::Pointer of generated mask requested */ - virtual mitk::Image::ConstPointer GetMask() = 0; + mitk::Image::ConstPointer GetMask(unsigned int maskID); /** * @brief GetReferenceImage per default returns the inputImage (as set by SetInputImage). If no input image is set it will return a nullptr. */ virtual mitk::Image::ConstPointer GetReferenceImage(); /** * @brief SetInputImage is used to set the input image to the mask generator. Some subclasses require an input image, others don't. See the documentation of the specific Mask Generator for more information. */ - void SetInputImage(mitk::Image::ConstPointer inputImg); - - virtual void SetTimeStep(unsigned int timeStep); + itkSetConstObjectMacro(InputImage, mitk::Image); + itkSetMacro(TimePoint, TimePointType); protected: MaskGenerator(); - unsigned int m_TimeStep; - mitk::Image::ConstPointer m_inputImage; + /** + * @brief DoGetMask must be overridden by derived classes. + * @param maskID Parameter indicating which mask should be returned. + * @return mitk::Image::Pointer of generated mask requested + */ + virtual mitk::Image::ConstPointer DoGetMask(unsigned int maskID) = 0; + + + TimePointType m_TimePoint; + mitk::Image::ConstPointer m_InputImage; private: }; } #endif diff --git a/Modules/ImageStatistics/mitkMultiLabelMaskGenerator.cpp b/Modules/ImageStatistics/mitkMultiLabelMaskGenerator.cpp index b615da27de..4a70fef445 100644 --- a/Modules/ImageStatistics/mitkMultiLabelMaskGenerator.cpp +++ b/Modules/ImageStatistics/mitkMultiLabelMaskGenerator.cpp @@ -1,13 +1,28 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #include +#include + +unsigned int mitk::MultiLabelMaskGenerator::GetNumberOfMasks() const +{ + if (m_MultiLabelSegmentation.IsNull()) mitkThrow() << "Invalid state. Cannot get number of masks. MultiLabelSegmentation is not set."; + return m_MultiLabelSegmentation->GetNumberOfLayers(); +} + +mitk::Image::ConstPointer mitk::MultiLabelMaskGenerator::DoGetMask(unsigned int maskID) +{ + if (m_MultiLabelSegmentation.IsNull()) mitkThrow() << "Invalid state. Cannot get number of masks. MultiLabelSegmentation is not set."; + + auto groupImage = m_MultiLabelSegmentation->GetGroupImage(maskID); + return SelectImageByTimePoint(groupImage, m_TimePoint); +} diff --git a/Modules/ImageStatistics/mitkMultiLabelMaskGenerator.h b/Modules/ImageStatistics/mitkMultiLabelMaskGenerator.h index f404255c5c..ed85537c6e 100644 --- a/Modules/ImageStatistics/mitkMultiLabelMaskGenerator.h +++ b/Modules/ImageStatistics/mitkMultiLabelMaskGenerator.h @@ -1,55 +1,49 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #ifndef mitkMultiLabelMaskGenerator_h #define mitkMultiLabelMaskGenerator_h #include #include -#include #include namespace mitk { /** - * @brief The MultiLabelMaskGenerator class NOT IMPLEMENTED YET! + * @brief Class that allows to generate masks (for statistic computation) out of multi label segmentations */ class MITKIMAGESTATISTICS_EXPORT MultiLabelMaskGenerator: public MaskGenerator { public: - /* void setLabelSetImage(mitk::LabelSetImage::Pointer labelSetImage); + /** Standard Self typedef */ + mitkClassMacro(MultiLabelMaskGenerator, MaskGenerator); + itkNewMacro(Self); - void addLabel(LabelSetImage::PixelType, std::vector::size_type layer=0); - void removeLabel(LabelSetImage::PixelType, std::vector::size_type layer=0); + unsigned int GetNumberOfMasks() const override; - void addLabels(std::pair::size_type, std::vector> labelsToAdd); - void removeLabels(std::pair::size_type, std::vector> labelsToAdd); - - void addLabels(std::vector labels, std::vector::size_type layer=0); - void removeLabels(std::vector labels, std::vector::size_type layer=0); - - void removeLayer(std::vector::size_type layer); - - mitk::Image::Pointer GetMask(); + itkSetConstObjectMacro(MultiLabelSegmentation, LabelSetImage); protected: + MultiLabelMaskGenerator() = default; + ~MultiLabelMaskGenerator() = default; -private: - mitk::LabelSetImage::Pointer m_LabelSetImage; - std::vector> m_selectedLabels;*/ + Image::ConstPointer DoGetMask(unsigned int) override; +private: + mitk::LabelSetImage::ConstPointer m_MultiLabelSegmentation; }; } #endif diff --git a/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.cpp b/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.cpp index 6b28bcaab7..6091e37d06 100644 --- a/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.cpp +++ b/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.cpp @@ -1,513 +1,495 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #include #include #include #include "mitkImageAccessByItk.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include namespace mitk { void PlanarFigureMaskGenerator::SetPlanarFigure(mitk::PlanarFigure* planarFigure) { if (nullptr == planarFigure ) { throw std::runtime_error( "Error: planar figure empty!" ); } const PlaneGeometry *planarFigurePlaneGeometry = planarFigure->GetPlaneGeometry(); if ( planarFigurePlaneGeometry == nullptr ) { throw std::runtime_error( "Planar-Figure not yet initialized!" ); } const auto *planarFigureGeometry = dynamic_cast< const PlaneGeometry * >( planarFigurePlaneGeometry ); if ( planarFigureGeometry == nullptr ) { throw std::runtime_error( "Non-planar planar figures not supported!" ); } if (planarFigure != m_PlanarFigure) { this->Modified(); m_PlanarFigure = planarFigure; } } mitk::Image::ConstPointer PlanarFigureMaskGenerator::GetReferenceImage() { if (IsUpdateRequired()) { this->CalculateMask(); } return m_ReferenceImage; } template < typename TPixel, unsigned int VImageDimension > void PlanarFigureMaskGenerator::InternalCalculateMaskFromPlanarFigure( const itk::Image< TPixel, VImageDimension > *image, unsigned int axis ) { typedef itk::Image< unsigned short, 2 > MaskImage2DType; typename MaskImage2DType::Pointer maskImage = MaskImage2DType::New(); maskImage->SetOrigin(image->GetOrigin()); maskImage->SetSpacing(image->GetSpacing()); maskImage->SetLargestPossibleRegion(image->GetLargestPossibleRegion()); maskImage->SetBufferedRegion(image->GetBufferedRegion()); maskImage->SetDirection(image->GetDirection()); maskImage->SetNumberOfComponentsPerPixel(image->GetNumberOfComponentsPerPixel()); maskImage->Allocate(); maskImage->FillBuffer(1); // all PolylinePoints of the PlanarFigure are stored in a vtkPoints object. // These points are used by the vtkLassoStencilSource to create // a vtkImageStencil. const mitk::PlaneGeometry *planarFigurePlaneGeometry = m_PlanarFigure->GetPlaneGeometry(); const typename PlanarFigure::PolyLineType planarFigurePolyline = m_PlanarFigure->GetPolyLine( 0 ); - const mitk::BaseGeometry *imageGeometry3D = m_inputImage->GetGeometry( 0 ); + const mitk::BaseGeometry *imageGeometry3D = m_InputImage->GetGeometry( 0 ); // If there is a second poly line in a closed planar figure, treat it as a hole. PlanarFigure::PolyLineType planarFigureHolePolyline; if (m_PlanarFigure->GetPolyLinesSize() == 2) planarFigureHolePolyline = m_PlanarFigure->GetPolyLine(1); // Determine x- and y-dimensions depending on principal axis // TODO use plane geometry normal to determine that automatically, then check whether the PF is aligned with one of the three principal axis int i0, i1; switch ( axis ) { case 0: i0 = 1; i1 = 2; break; case 1: i0 = 0; i1 = 2; break; case 2: default: i0 = 0; i1 = 1; break; } // store the polyline contour as vtkPoints object vtkSmartPointer points = vtkSmartPointer::New(); for (const auto& point : planarFigurePolyline) { Point3D point3D; // Convert 2D point back to the local index coordinates of the selected image planarFigurePlaneGeometry->Map(point, point3D); imageGeometry3D->WorldToIndex(point3D, point3D); points->InsertNextPoint(point3D[i0], point3D[i1], 0); } vtkSmartPointer holePoints; if (!planarFigureHolePolyline.empty()) { holePoints = vtkSmartPointer::New(); Point3D point3D; for (const auto& point : planarFigureHolePolyline) { planarFigurePlaneGeometry->Map(point, point3D); imageGeometry3D->WorldToIndex(point3D, point3D); holePoints->InsertNextPoint(point3D[i0], point3D[i1], 0); } } // mark a malformed 2D planar figure ( i.e. area = 0 ) as out of bounds // this can happen when all control points of a rectangle lie on the same line = two of the three extents are zero double bounds[6] = {0}; points->GetBounds(bounds); bool extent_x = (fabs(bounds[0] - bounds[1])) < mitk::eps; bool extent_y = (fabs(bounds[2] - bounds[3])) < mitk::eps; bool extent_z = (fabs(bounds[4] - bounds[5])) < mitk::eps; // throw an exception if a closed planar figure is deformed, i.e. has only one non-zero extent if (m_PlanarFigure->IsClosed() && ((extent_x && extent_y) || (extent_x && extent_z) || (extent_y && extent_z))) { mitkThrow() << "Figure has a zero area and cannot be used for masking."; } // create a vtkLassoStencilSource and set the points of the Polygon vtkSmartPointer lassoStencil = vtkSmartPointer::New(); lassoStencil->SetShapeToPolygon(); lassoStencil->SetPoints(points); vtkSmartPointer holeLassoStencil = nullptr; if (holePoints.GetPointer() != nullptr) { holeLassoStencil = vtkSmartPointer::New(); holeLassoStencil->SetShapeToPolygon(); holeLassoStencil->SetPoints(holePoints); } // Export from ITK to VTK (to use a VTK filter) typedef itk::VTKImageImport< MaskImage2DType > ImageImportType; typedef itk::VTKImageExport< MaskImage2DType > ImageExportType; typename ImageExportType::Pointer itkExporter = ImageExportType::New(); itkExporter->SetInput( maskImage ); // itkExporter->SetInput( castFilter->GetOutput() ); vtkSmartPointer vtkImporter = vtkSmartPointer::New(); this->ConnectPipelines( itkExporter, vtkImporter ); // Apply the generated image stencil to the input image vtkSmartPointer imageStencilFilter = vtkSmartPointer::New(); imageStencilFilter->SetInputConnection( vtkImporter->GetOutputPort() ); imageStencilFilter->SetStencilConnection(lassoStencil->GetOutputPort()); imageStencilFilter->ReverseStencilOff(); imageStencilFilter->SetBackgroundValue( 0 ); imageStencilFilter->Update(); vtkSmartPointer holeStencilFilter = nullptr; if (holeLassoStencil.GetPointer() != nullptr) { holeStencilFilter = vtkSmartPointer::New(); holeStencilFilter->SetInputConnection(imageStencilFilter->GetOutputPort()); holeStencilFilter->SetStencilConnection(holeLassoStencil->GetOutputPort()); holeStencilFilter->ReverseStencilOn(); holeStencilFilter->SetBackgroundValue(0); holeStencilFilter->Update(); } // Export from VTK back to ITK vtkSmartPointer vtkExporter = vtkSmartPointer::New(); vtkExporter->SetInputConnection( holeStencilFilter.GetPointer() == nullptr ? imageStencilFilter->GetOutputPort() : holeStencilFilter->GetOutputPort()); vtkExporter->Update(); typename ImageImportType::Pointer itkImporter = ImageImportType::New(); this->ConnectPipelines( vtkExporter, itkImporter ); itkImporter->Update(); typedef itk::ImageDuplicator< ImageImportType::OutputImageType > DuplicatorType; DuplicatorType::Pointer duplicator = DuplicatorType::New(); duplicator->SetInputImage( itkImporter->GetOutput() ); duplicator->Update(); // Store mask m_InternalITKImageMask2D = duplicator->GetOutput(); } template < typename TPixel, unsigned int VImageDimension > void PlanarFigureMaskGenerator::InternalCalculateMaskFromOpenPlanarFigure( const itk::Image< TPixel, VImageDimension > *image, unsigned int axis ) { typedef itk::Image< unsigned short, 2 > MaskImage2DType; typedef itk::LineIterator< MaskImage2DType > LineIteratorType; typedef MaskImage2DType::IndexType IndexType2D; typedef std::vector< IndexType2D > IndexVecType; typename MaskImage2DType::Pointer maskImage = MaskImage2DType::New(); maskImage->SetOrigin(image->GetOrigin()); maskImage->SetSpacing(image->GetSpacing()); maskImage->SetLargestPossibleRegion(image->GetLargestPossibleRegion()); maskImage->SetBufferedRegion(image->GetBufferedRegion()); maskImage->SetDirection(image->GetDirection()); maskImage->SetNumberOfComponentsPerPixel(image->GetNumberOfComponentsPerPixel()); maskImage->Allocate(); maskImage->FillBuffer(0); // all PolylinePoints of the PlanarFigure are stored in a vtkPoints object. const mitk::PlaneGeometry *planarFigurePlaneGeometry = m_PlanarFigure->GetPlaneGeometry(); const typename PlanarFigure::PolyLineType planarFigurePolyline = m_PlanarFigure->GetPolyLine( 0 ); - const mitk::BaseGeometry *imageGeometry3D = m_inputImage->GetGeometry( 0 ); + const mitk::BaseGeometry *imageGeometry3D = m_InputImage->GetGeometry( 0 ); // Determine x- and y-dimensions depending on principal axis // TODO use plane geometry normal to determine that automatically, then check whether the PF is aligned with one of the three principal axis int i0, i1; switch ( axis ) { case 0: i0 = 1; i1 = 2; break; case 1: i0 = 0; i1 = 2; break; case 2: default: i0 = 0; i1 = 1; break; } int numPolyLines = m_PlanarFigure->GetPolyLinesSize(); for ( int lineId = 0; lineId < numPolyLines; ++lineId ) { // store the polyline contour as vtkPoints object IndexVecType pointIndices; for(const auto& point : planarFigurePolyline) { Point3D point3D; planarFigurePlaneGeometry->Map(point, point3D); imageGeometry3D->WorldToIndex(point3D, point3D); IndexType2D index2D; index2D[0] = point3D[i0]; index2D[1] = point3D[i1]; pointIndices.push_back( index2D ); } size_t numLineSegments = pointIndices.size() - 1; for (size_t i = 0; i < numLineSegments; ++i) { LineIteratorType lineIt(maskImage, pointIndices[i], pointIndices[i+1]); while (!lineIt.IsAtEnd()) { lineIt.Set(1); ++lineIt; } } } // Store mask m_InternalITKImageMask2D = maskImage; } bool PlanarFigureMaskGenerator::CheckPlanarFigureIsNotTilted(const PlaneGeometry* planarGeometry, const BaseGeometry *geometry) { if (!planarGeometry) return false; if (!geometry) return false; unsigned int axis; return GetPrincipalAxis(geometry,planarGeometry->GetNormal(), axis); } bool PlanarFigureMaskGenerator::GetPrincipalAxis( const BaseGeometry *geometry, Vector3D vector, unsigned int &axis ) { vector.Normalize(); for ( unsigned int i = 0; i < 3; ++i ) { Vector3D axisVector = geometry->GetAxisVector( i ); axisVector.Normalize(); //normal mitk::eps is to pedantic for this check. See e.g. T27122 //therefore choose a larger epsilon. The value was set a) as small as //possible but b) still allowing to datasets like in (T27122) to pass //when floating rounding errors sum up. const double epsilon = 5e-5; if ( fabs( fabs( axisVector * vector ) - 1.0) < epsilon) { axis = i; return true; } } return false; } void PlanarFigureMaskGenerator::CalculateMask() { - if (m_inputImage.IsNull()) + if (m_InputImage.IsNull()) { - MITK_ERROR << "Image is not set."; + mitkThrow() << "Image is not set."; } if (m_PlanarFigure.IsNull()) { - MITK_ERROR << "PlanarFigure is not set."; + mitkThrow() << "PlanarFigure is not set."; } - if (m_TimeStep != 0) - { - MITK_WARN << "Multiple TimeSteps are not supported in PlanarFigureMaskGenerator (yet)."; - } - - const BaseGeometry *imageGeometry = m_inputImage->GetGeometry(); + const BaseGeometry *imageGeometry = m_InputImage->GetGeometry(); if ( imageGeometry == nullptr ) { - throw std::runtime_error( "Image geometry invalid!" ); + mitkThrow() << "Image geometry invalid!"; } - if (m_inputImage->GetTimeSteps() > 0) - { - mitk::ImageTimeSelector::Pointer imgTimeSel = mitk::ImageTimeSelector::New(); - imgTimeSel->SetInput(m_inputImage); - imgTimeSel->SetTimeNr(m_TimeStep); - imgTimeSel->UpdateLargestPossibleRegion(); - m_InternalTimeSliceImage = imgTimeSel->GetOutput(); - } - else - { - m_InternalTimeSliceImage = m_inputImage; - } + auto timePointImage = SelectImageByTimePoint(m_InputImage, m_TimePoint); + + if (timePointImage.IsNull()) mitkThrow() << "Cannot generate mask. Passed time point is not supported by input image."; m_InternalITKImageMask2D = nullptr; const PlaneGeometry *planarFigurePlaneGeometry = m_PlanarFigure->GetPlaneGeometry(); const auto *planarFigureGeometry = dynamic_cast< const PlaneGeometry * >( planarFigurePlaneGeometry ); //const BaseGeometry *imageGeometry = m_inputImage->GetGeometry(); // Find principal direction of PlanarFigure in input image unsigned int axis; if ( !this->GetPrincipalAxis( imageGeometry, planarFigureGeometry->GetNormal(), axis ) ) { throw std::runtime_error( "Non-aligned planar figures not supported!" ); } m_PlanarFigureAxis = axis; // Find slice number corresponding to PlanarFigure in input image itk::Image< unsigned short, 3 >::IndexType index; imageGeometry->WorldToIndex( planarFigureGeometry->GetOrigin(), index ); unsigned int slice = index[axis]; m_PlanarFigureSlice = slice; // extract image slice which corresponds to the planarFigure and store it in m_InternalImageSlice - mitk::Image::ConstPointer inputImageSlice = extract2DImageSlice(axis, slice); + mitk::Image::ConstPointer inputImageSlice = Extract2DImageSlice(timePointImage, axis, slice); //mitk::IOUtil::Save(inputImageSlice, "/home/fabian/inputSliceImage.nrrd"); // Compute mask from PlanarFigure // rastering for open planar figure: if ( !m_PlanarFigure->IsClosed() ) { AccessFixedDimensionByItk_1(inputImageSlice, InternalCalculateMaskFromOpenPlanarFigure, 2, axis) } else//for closed planar figure { AccessFixedDimensionByItk_1(inputImageSlice, InternalCalculateMaskFromPlanarFigure, 2, axis) } //convert itk mask to mitk::Image::Pointer and return it mitk::Image::Pointer planarFigureMaskImage; planarFigureMaskImage = mitk::GrabItkImageMemory(m_InternalITKImageMask2D); //mitk::IOUtil::Save(planarFigureMaskImage, "/home/fabian/planarFigureMaskImage.nrrd"); //Convert2Dto3DImageFilter::Pointer sliceTo3DImageConverter = Convert2Dto3DImageFilter::New(); //sliceTo3DImageConverter->SetInput(planarFigureMaskImage); //sliceTo3DImageConverter->Update(); //mitk::IOUtil::Save(sliceTo3DImageConverter->GetOutput(), "/home/fabian/3DsliceImage.nrrd"); m_ReferenceImage = inputImageSlice; //mitk::IOUtil::Save(m_ReferenceImage, "/home/fabian/referenceImage.nrrd"); m_InternalMask = planarFigureMaskImage; } -void PlanarFigureMaskGenerator::SetTimeStep(unsigned int timeStep) +unsigned int PlanarFigureMaskGenerator::GetNumberOfMasks() const { - if (timeStep != m_TimeStep) - { - m_TimeStep = timeStep; - this->Modified(); - } + return 1; } -mitk::Image::ConstPointer PlanarFigureMaskGenerator::GetMask() +mitk::Image::ConstPointer PlanarFigureMaskGenerator::DoGetMask(unsigned int) { if (IsUpdateRequired()) { this->CalculateMask(); this->Modified(); } m_InternalMaskUpdateTime = this->GetMTime(); return m_InternalMask; } -mitk::Image::ConstPointer PlanarFigureMaskGenerator::extract2DImageSlice(unsigned int axis, unsigned int slice) +mitk::Image::ConstPointer PlanarFigureMaskGenerator::Extract2DImageSlice(const Image* input, unsigned int axis, unsigned int slice) const { // Extract slice with given position and direction from image - unsigned int dimension = m_InternalTimeSliceImage->GetDimension(); + unsigned int dimension = input->GetDimension(); if (dimension == 3) { ExtractImageFilter::Pointer imageExtractor = ExtractImageFilter::New(); - imageExtractor->SetInput( m_InternalTimeSliceImage ); + imageExtractor->SetInput( input ); imageExtractor->SetSliceDimension( axis ); imageExtractor->SetSliceIndex( slice ); imageExtractor->Update(); return imageExtractor->GetOutput(); } else if(dimension == 2) { - return m_InternalTimeSliceImage; + return input; } else { MITK_ERROR << "Unsupported image dimension. Dimension is: " << dimension << ". Only 2D and 3D images are supported."; return nullptr; } } bool PlanarFigureMaskGenerator::IsUpdateRequired() const { unsigned long thisClassTimeStamp = this->GetMTime(); unsigned long internalMaskTimeStamp = m_InternalMask->GetMTime(); unsigned long planarFigureTimeStamp = m_PlanarFigure->GetMTime(); - unsigned long inputImageTimeStamp = m_inputImage->GetMTime(); + unsigned long inputImageTimeStamp = m_InputImage->GetMTime(); if (thisClassTimeStamp > m_InternalMaskUpdateTime) // inputs have changed { return true; } if (m_InternalMaskUpdateTime < planarFigureTimeStamp || m_InternalMaskUpdateTime < inputImageTimeStamp) // mask image has changed outside of this class { return true; } if (internalMaskTimeStamp > m_InternalMaskUpdateTime) // internal mask has been changed outside of this class { return true; } return false; } } diff --git a/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.h b/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.h index 67a4f6bbf4..1471d0cb0f 100644 --- a/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.h +++ b/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.h @@ -1,146 +1,137 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #ifndef mitkPlanarFigureMaskGenerator_h #define mitkPlanarFigureMaskGenerator_h #include #include #include #include #include #include namespace mitk { /** * \class PlanarFigureMaskGenerator * \brief Derived from MaskGenerator. This class is used to convert a mitk::PlanarFigure into a binary image mask */ class MITKIMAGESTATISTICS_EXPORT PlanarFigureMaskGenerator : public MaskGenerator { public: /** Standard Self typedef */ typedef PlanarFigureMaskGenerator Self; typedef MaskGenerator Superclass; typedef itk::SmartPointer Pointer; typedef itk::SmartPointer ConstPointer; /** Method for creation through the object factory. */ itkNewMacro(Self); /** Runtime information support. */ itkTypeMacro(PlanarFigureMaskGenerator, MaskGenerator); - /** - * @brief GetMask Computes and returns the mask - * @return mitk::Image::Pointer of the generated mask - */ - mitk::Image::ConstPointer GetMask() override; + unsigned int GetNumberOfMasks() const override; void SetPlanarFigure(mitk::PlanarFigure* planarFigure); mitk::Image::ConstPointer GetReferenceImage() override; - /** - * @brief SetTimeStep is used to set the time step for which the mask is to be generated - * @param timeStep - */ - void SetTimeStep(unsigned int timeStep) override; - itkGetConstMacro(PlanarFigureAxis, unsigned int); itkGetConstMacro(PlanarFigureSlice, unsigned int); /** Helper function that indicates if a passed planar geometry is tilted regarding a given geometry and its main axis. *@pre If either planarGeometry or geometry is nullptr it will return false.*/ static bool CheckPlanarFigureIsNotTilted(const PlaneGeometry* planarGeometry, const BaseGeometry *geometry); protected: PlanarFigureMaskGenerator() : Superclass(), m_ReferenceImage(nullptr), m_PlanarFigureAxis(0), m_InternalMaskUpdateTime(0), m_PlanarFigureSlice(0) { m_InternalMask = mitk::Image::New(); } + Image::ConstPointer DoGetMask(unsigned int) override; + private: void CalculateMask(); template void InternalCalculateMaskFromPlanarFigure(const itk::Image *image, unsigned int axis); template void InternalCalculateMaskFromOpenPlanarFigure(const itk::Image *image, unsigned int axis); - mitk::Image::ConstPointer extract2DImageSlice(unsigned int axis, unsigned int slice); + mitk::Image::ConstPointer Extract2DImageSlice(const Image* input, unsigned int axis, unsigned int slice) const; /** Helper function that deduces if the passed vector is equal to one of the primary axis of the geometry.*/ static bool GetPrincipalAxis(const BaseGeometry *geometry, Vector3D vector, unsigned int &axis); /** Connection from ITK to VTK */ template void ConnectPipelines(ITK_Exporter exporter, vtkSmartPointer importer) { importer->SetUpdateInformationCallback(exporter->GetUpdateInformationCallback()); importer->SetPipelineModifiedCallback(exporter->GetPipelineModifiedCallback()); importer->SetWholeExtentCallback(exporter->GetWholeExtentCallback()); importer->SetSpacingCallback(exporter->GetSpacingCallback()); importer->SetOriginCallback(exporter->GetOriginCallback()); importer->SetScalarTypeCallback(exporter->GetScalarTypeCallback()); importer->SetNumberOfComponentsCallback(exporter->GetNumberOfComponentsCallback()); importer->SetPropagateUpdateExtentCallback(exporter->GetPropagateUpdateExtentCallback()); importer->SetUpdateDataCallback(exporter->GetUpdateDataCallback()); importer->SetDataExtentCallback(exporter->GetDataExtentCallback()); importer->SetBufferPointerCallback(exporter->GetBufferPointerCallback()); importer->SetCallbackUserData(exporter->GetCallbackUserData()); } /** Connection from VTK to ITK */ template void ConnectPipelines(vtkSmartPointer exporter, ITK_Importer importer) { importer->SetUpdateInformationCallback(exporter->GetUpdateInformationCallback()); importer->SetPipelineModifiedCallback(exporter->GetPipelineModifiedCallback()); importer->SetWholeExtentCallback(exporter->GetWholeExtentCallback()); importer->SetSpacingCallback(exporter->GetSpacingCallback()); importer->SetOriginCallback(exporter->GetOriginCallback()); importer->SetScalarTypeCallback(exporter->GetScalarTypeCallback()); importer->SetNumberOfComponentsCallback(exporter->GetNumberOfComponentsCallback()); importer->SetPropagateUpdateExtentCallback(exporter->GetPropagateUpdateExtentCallback()); importer->SetUpdateDataCallback(exporter->GetUpdateDataCallback()); importer->SetDataExtentCallback(exporter->GetDataExtentCallback()); importer->SetBufferPointerCallback(exporter->GetBufferPointerCallback()); importer->SetCallbackUserData(exporter->GetCallbackUserData()); } bool IsUpdateRequired() const; mitk::PlanarFigure::Pointer m_PlanarFigure; itk::Image::Pointer m_InternalITKImageMask2D; - mitk::Image::ConstPointer m_InternalTimeSliceImage; mitk::Image::ConstPointer m_ReferenceImage; unsigned int m_PlanarFigureAxis; unsigned long m_InternalMaskUpdateTime; unsigned int m_PlanarFigureSlice; mitk::Image::Pointer m_InternalMask; }; } // namespace mitk #endif diff --git a/Modules/ImageStatisticsUI/Qmitk/QmitkImageStatisticsCalculationRunnable.cpp b/Modules/ImageStatisticsUI/Qmitk/QmitkImageStatisticsCalculationRunnable.cpp index 7d8ad4cbc6..0223ce5987 100644 --- a/Modules/ImageStatisticsUI/Qmitk/QmitkImageStatisticsCalculationRunnable.cpp +++ b/Modules/ImageStatisticsUI/Qmitk/QmitkImageStatisticsCalculationRunnable.cpp @@ -1,183 +1,187 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #include "QmitkImageStatisticsCalculationRunnable.h" #include "mitkImageStatisticsCalculator.h" -#include +#include +#include +#include #include +#include +#include #include #include "mitkStatisticsToImageRelationRule.h" #include "mitkStatisticsToMaskRelationRule.h" #include "mitkImageStatisticsContainerManager.h" #include "mitkProperties.h" QmitkImageStatisticsCalculationRunnable::QmitkImageStatisticsCalculationRunnable() : QmitkDataGenerationJobBase() , m_StatisticsImage(nullptr) - , m_BinaryMask(nullptr) - , m_PlanarFigureMask(nullptr) + , m_MaskData(nullptr) , m_IgnoreZeros(false) , m_HistogramNBins(100) { } QmitkImageStatisticsCalculationRunnable::~QmitkImageStatisticsCalculationRunnable() { } -void QmitkImageStatisticsCalculationRunnable::Initialize(const mitk::Image *image, - const mitk::Image *binaryImage, - const mitk::PlanarFigure *planarFig) +void QmitkImageStatisticsCalculationRunnable::Initialize(const mitk::Image* image, const mitk::BaseData* mask) { + if (nullptr!= mask && + nullptr == dynamic_cast(mask) && + nullptr == dynamic_cast(mask) && + nullptr == dynamic_cast(mask)) + mitkThrow() << "Cannot initialize QmitkImageStatisticsCalculationRunnable. Mask data is not a supported type."; this->m_StatisticsImage = image; - this->m_BinaryMask = binaryImage; - this->m_PlanarFigureMask = planarFig; + this->m_MaskData = mask; } mitk::ImageStatisticsContainer* QmitkImageStatisticsCalculationRunnable::GetStatisticsData() const { return this->m_StatisticsContainer.GetPointer(); } const mitk::Image* QmitkImageStatisticsCalculationRunnable::GetStatisticsImage() const { return this->m_StatisticsImage.GetPointer(); } -const mitk::Image* QmitkImageStatisticsCalculationRunnable::GetMaskImage() const +const mitk::BaseData* QmitkImageStatisticsCalculationRunnable::GetMaskData() const { - return this->m_BinaryMask.GetPointer(); -} - -const mitk::PlanarFigure* QmitkImageStatisticsCalculationRunnable::GetPlanarFigure() const -{ - return this->m_PlanarFigureMask.GetPointer(); + return this->m_MaskData.GetPointer(); } void QmitkImageStatisticsCalculationRunnable::SetIgnoreZeroValueVoxel(bool _arg) { this->m_IgnoreZeros = _arg; } bool QmitkImageStatisticsCalculationRunnable::GetIgnoreZeroValueVoxel() const { return this->m_IgnoreZeros; } void QmitkImageStatisticsCalculationRunnable::SetHistogramNBins(unsigned int nbins) { this->m_HistogramNBins = nbins; } unsigned int QmitkImageStatisticsCalculationRunnable::GetHistogramNBins() const { return this->m_HistogramNBins; } QmitkDataGenerationJobBase::ResultMapType QmitkImageStatisticsCalculationRunnable::GetResults() const { ResultMapType result; result.emplace("statistics", this->GetStatisticsData()); return result; } bool QmitkImageStatisticsCalculationRunnable::RunComputation() { bool statisticCalculationSuccessful = true; mitk::ImageStatisticsCalculator::Pointer calculator = mitk::ImageStatisticsCalculator::New(); if (this->m_StatisticsImage.IsNotNull()) { calculator->SetInputImage(m_StatisticsImage); } 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()) + auto multiLabelMask = dynamic_cast(m_MaskData.GetPointer()); + auto binLabelMask = dynamic_cast(m_MaskData.GetPointer()); + auto pfMask = dynamic_cast(m_MaskData.GetPointer()); + + if (nullptr != multiLabelMask) + { + mitk::MultiLabelMaskGenerator::Pointer imgMask = mitk::MultiLabelMaskGenerator::New(); + imgMask->SetMultiLabelSegmentation(multiLabelMask); + calculator->SetMask(imgMask.GetPointer()); + } + else if (nullptr != binLabelMask) { mitk::ImageMaskGenerator::Pointer imgMask = mitk::ImageMaskGenerator::New(); imgMask->SetInputImage(m_StatisticsImage); - imgMask->SetImageMask(m_BinaryMask); + imgMask->SetImageMask(binLabelMask); calculator->SetMask(imgMask.GetPointer()); } - if (this->m_PlanarFigureMask.IsNotNull()) + else if (nullptr != pfMask) { mitk::PlanarFigureMaskGenerator::Pointer pfMaskGen = mitk::PlanarFigureMaskGenerator::New(); pfMaskGen->SetInputImage(m_StatisticsImage); - pfMaskGen->SetPlanarFigure(m_PlanarFigureMask->Clone()); + pfMaskGen->SetPlanarFigure(pfMask->Clone()); calculator->SetMask(pfMaskGen.GetPointer()); } } catch (const std::exception &e) { MITK_ERROR << "Error while configuring the statistics calculator: " << e.what(); m_LastErrorMessage = e.what(); statisticCalculationSuccessful = false; } if (this->m_IgnoreZeros) { mitk::IgnorePixelMaskGenerator::Pointer ignorePixelValueMaskGen = mitk::IgnorePixelMaskGenerator::New(); ignorePixelValueMaskGen->SetIgnoredPixelValue(0); ignorePixelValueMaskGen->SetInputImage(m_StatisticsImage); calculator->SetSecondaryMask(ignorePixelValueMaskGen.GetPointer()); } else { calculator->SetSecondaryMask(nullptr); } calculator->SetNBinsForHistogramStatistics(m_HistogramNBins); try { calculator->GetStatistics(); } catch (const std::exception &e) { m_LastErrorMessage = "Failure while calculating the statistics: " + std::string(e.what()); MITK_ERROR << m_LastErrorMessage; statisticCalculationSuccessful = false; } if (statisticCalculationSuccessful) { m_StatisticsContainer = calculator->GetStatistics(); auto imageRule = mitk::StatisticsToImageRelationRule::New(); imageRule->Connect(m_StatisticsContainer, m_StatisticsImage); - if (nullptr != m_PlanarFigureMask) - { - auto maskRule = mitk::StatisticsToMaskRelationRule::New(); - maskRule->Connect(m_StatisticsContainer, m_PlanarFigureMask); - } - - if (nullptr != m_BinaryMask) + if (nullptr != m_MaskData) { auto maskRule = mitk::StatisticsToMaskRelationRule::New(); - maskRule->Connect(m_StatisticsContainer, m_BinaryMask); + maskRule->Connect(m_StatisticsContainer, m_MaskData); } m_StatisticsContainer->SetProperty(mitk::STATS_HISTOGRAM_BIN_PROPERTY_NAME.c_str(), mitk::UIntProperty::New(m_HistogramNBins)); m_StatisticsContainer->SetProperty(mitk::STATS_IGNORE_ZERO_VOXEL_PROPERTY_NAME.c_str(), mitk::BoolProperty::New(m_IgnoreZeros)); } return statisticCalculationSuccessful; } diff --git a/Modules/ImageStatisticsUI/Qmitk/QmitkImageStatisticsCalculationRunnable.h b/Modules/ImageStatisticsUI/Qmitk/QmitkImageStatisticsCalculationRunnable.h index 9054dd0676..db41032275 100644 --- a/Modules/ImageStatisticsUI/Qmitk/QmitkImageStatisticsCalculationRunnable.h +++ b/Modules/ImageStatisticsUI/Qmitk/QmitkImageStatisticsCalculationRunnable.h @@ -1,86 +1,82 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #ifndef QmitkImageStatisticsCalculationRunnable_h #define QmitkImageStatisticsCalculationRunnable_h //mitk headers -#include "mitkImage.h" -#include "mitkPlanarFigure.h" -#include "mitkImageStatisticsContainer.h" +#include #include "QmitkDataGenerationJobBase.h" // itk headers #ifndef __itkHistogram_h #include #endif #include /** * /brief This class is executed as background thread for image statistics calculation. * * This class is derived from QRunnable and is intended to be used by QmitkImageStatisticsView * to run the image statistics calculation in a background thread keeping the GUI usable. */ class MITKIMAGESTATISTICSUI_EXPORT QmitkImageStatisticsCalculationRunnable : public QmitkDataGenerationJobBase { Q_OBJECT public: typedef itk::Statistics::Histogram HistogramType; /*! /brief standard constructor. */ QmitkImageStatisticsCalculationRunnable(); /*! /brief standard destructor. */ ~QmitkImageStatisticsCalculationRunnable(); /*! /brief Initializes the object with necessary data. */ - void Initialize(const mitk::Image* image, const mitk::Image* binaryImage, const mitk::PlanarFigure* planarFig); + void Initialize(const mitk::Image* image, const mitk::BaseData* mask); /*! /brief returns the calculated image statistics. */ mitk::ImageStatisticsContainer* GetStatisticsData() const; const mitk::Image* GetStatisticsImage() const; - const mitk::Image* GetMaskImage() const; - const mitk::PlanarFigure* GetPlanarFigure() const; + const mitk::BaseData* GetMaskData() const; /*! /brief Set flag to ignore zero valued voxels */ void SetIgnoreZeroValueVoxel(bool _arg); /*! /brief Get status of zero value voxel ignoring. */ bool GetIgnoreZeroValueVoxel() const; /*! /brief Set bin size for histogram resolution.*/ void SetHistogramNBins(unsigned int nbins); /*! /brief Get bin size for histogram resolution.*/ unsigned int GetHistogramNBins() const; ResultMapType GetResults() const override; protected: bool RunComputation() override; private: mitk::Image::ConstPointer m_StatisticsImage; ///< member variable holds the input image for which the statistics need to be calculated. - mitk::Image::ConstPointer m_BinaryMask; ///< member variable holds the binary mask image for segmentation image statistics calculation. - mitk::PlanarFigure::ConstPointer m_PlanarFigureMask; ///< member variable holds the planar figure for segmentation image statistics calculation. + mitk::BaseData::ConstPointer m_MaskData; ///< member variable holds the data that should be used as mask statistics calculation. mitk::ImageStatisticsContainer::Pointer m_StatisticsContainer; bool m_IgnoreZeros; ///< member variable holds flag to indicate if zero valued voxel should be suppressed unsigned int m_HistogramNBins; ///< member variable holds the bin size for histogram resolution. }; #endif diff --git a/Modules/ImageStatisticsUI/Qmitk/QmitkImageStatisticsDataGenerator.cpp b/Modules/ImageStatisticsUI/Qmitk/QmitkImageStatisticsDataGenerator.cpp index ed2002dc05..d062860977 100644 --- a/Modules/ImageStatisticsUI/Qmitk/QmitkImageStatisticsDataGenerator.cpp +++ b/Modules/ImageStatisticsUI/Qmitk/QmitkImageStatisticsDataGenerator.cpp @@ -1,274 +1,267 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #include "QmitkImageStatisticsDataGenerator.h" #include "mitkImageStatisticsContainer.h" #include "mitkStatisticsToImageRelationRule.h" #include "mitkStatisticsToMaskRelationRule.h" #include "mitkNodePredicateFunction.h" #include "mitkNodePredicateAnd.h" #include "mitkNodePredicateNot.h" #include "mitkNodePredicateDataProperty.h" #include "mitkProperties.h" #include "mitkImageStatisticsContainerManager.h" #include "QmitkImageStatisticsCalculationRunnable.h" void QmitkImageStatisticsDataGenerator::SetIgnoreZeroValueVoxel(bool _arg) { if (m_IgnoreZeroValueVoxel != _arg) { m_IgnoreZeroValueVoxel = _arg; if (m_AutoUpdate) { this->EnsureRecheckingAndGeneration(); } } } bool QmitkImageStatisticsDataGenerator::GetIgnoreZeroValueVoxel() const { return this->m_IgnoreZeroValueVoxel; } void QmitkImageStatisticsDataGenerator::SetHistogramNBins(unsigned int nbins) { if (m_HistogramNBins != nbins) { m_HistogramNBins = nbins; if (m_AutoUpdate) { this->EnsureRecheckingAndGeneration(); } } } unsigned int QmitkImageStatisticsDataGenerator::GetHistogramNBins() const { return this->m_HistogramNBins; } bool QmitkImageStatisticsDataGenerator::ChangedNodeIsRelevant(const mitk::DataNode* changedNode) const { auto result = QmitkImageAndRoiDataGeneratorBase::ChangedNodeIsRelevant(changedNode); if (!result) { if (changedNode->GetProperty(mitk::STATS_GENERATION_STATUS_PROPERTY_NAME.c_str()) == nullptr) { auto stats = dynamic_cast(changedNode->GetData()); result = stats != nullptr && stats->GetMTime()> this->m_GenerationTime.GetMTime(); } } return result; } bool QmitkImageStatisticsDataGenerator::IsValidResultAvailable(const mitk::DataNode* imageNode, const mitk::DataNode* roiNode) const { auto resultNode = this->GetLatestResult(imageNode, roiNode, true, true); return resultNode.IsNotNull(); } mitk::DataNode::Pointer QmitkImageStatisticsDataGenerator::GetLatestResult(const mitk::DataNode* imageNode, const mitk::DataNode* roiNode, bool onlyIfUpToDate, bool noWIP) const { auto storage = m_Storage.Lock(); if (imageNode == nullptr || !imageNode->GetData()) { mitkThrow() << "Image is nullptr"; } const auto image = imageNode->GetData(); const mitk::BaseData* mask = nullptr; if (roiNode) { mask = roiNode->GetData(); } std::lock_guard mutexguard(m_DataMutex); return mitk::ImageStatisticsContainerManager::GetImageStatisticsNode(storage, image, mask, m_IgnoreZeroValueVoxel, m_HistogramNBins, onlyIfUpToDate, noWIP); } void QmitkImageStatisticsDataGenerator::IndicateFutureResults(const mitk::DataNode* imageNode, const mitk::DataNode* roiNode) const { if (imageNode == nullptr || !imageNode->GetData()) { mitkThrow() << "Image node is nullptr"; } auto image = dynamic_cast(imageNode->GetData()); if (!image) { mitkThrow() << "Image node date is nullptr or no image."; } const mitk::BaseData* mask = nullptr; if (roiNode != nullptr) { mask = roiNode->GetData(); } auto resultDataNode = this->GetLatestResult(imageNode, roiNode, true, false); if (resultDataNode.IsNull()) { auto dummyStats = mitk::ImageStatisticsContainer::New(); dummyStats->SetTimeGeometry(image->GetTimeGeometry()->Clone()); auto imageRule = mitk::StatisticsToImageRelationRule::New(); imageRule->Connect(dummyStats, image); if (nullptr != mask) { auto maskRule = mitk::StatisticsToMaskRelationRule::New(); maskRule->Connect(dummyStats, mask); } dummyStats->SetProperty(mitk::STATS_HISTOGRAM_BIN_PROPERTY_NAME.c_str(), mitk::UIntProperty::New(m_HistogramNBins)); dummyStats->SetProperty(mitk::STATS_IGNORE_ZERO_VOXEL_PROPERTY_NAME.c_str(), mitk::BoolProperty::New(m_IgnoreZeroValueVoxel)); auto dummyNode = CreateWIPDataNode(dummyStats, "WIP_"+GenerateStatisticsNodeName(image, mask)); auto storage = m_Storage.Lock(); if (storage != nullptr) { std::lock_guard mutexguard(m_DataMutex); storage->Add(dummyNode); } } } std::pair QmitkImageStatisticsDataGenerator::GetNextMissingGenerationJob(const mitk::DataNode* imageNode, const mitk::DataNode* roiNode) const { auto resultDataNode = this->GetLatestResult(imageNode, roiNode, true, false); std::string status; if (resultDataNode.IsNull() || (resultDataNode->GetStringProperty(mitk::STATS_GENERATION_STATUS_PROPERTY_NAME.c_str(), status) && status == mitk::STATS_GENERATION_STATUS_VALUE_PENDING)) { if (imageNode == nullptr || !imageNode->GetData()) { mitkThrow() << "Image node is nullptr"; } auto image = dynamic_cast(imageNode->GetData()); if (image == nullptr) { mitkThrow() << "Image node date is nullptr or no image."; } - const mitk::Image* mask = nullptr; - const mitk::PlanarFigure* planar = nullptr; + const mitk::BaseData* mask = nullptr; if (roiNode != nullptr) { - mask = dynamic_cast(roiNode->GetData()); - planar = dynamic_cast(roiNode->GetData()); + mask = roiNode->GetData(); } auto newJob = new QmitkImageStatisticsCalculationRunnable; - newJob->Initialize(image, mask, planar); + newJob->Initialize(image, mask); newJob->SetIgnoreZeroValueVoxel(m_IgnoreZeroValueVoxel); newJob->SetHistogramNBins(m_HistogramNBins); return std::pair(newJob, resultDataNode.GetPointer()); } else if (resultDataNode->GetStringProperty(mitk::STATS_GENERATION_STATUS_PROPERTY_NAME.c_str(), status) && status == mitk::STATS_GENERATION_STATUS_VALUE_WORK_IN_PROGRESS) { return std::pair(nullptr, resultDataNode.GetPointer()); } return std::pair(nullptr, nullptr); } void QmitkImageStatisticsDataGenerator::RemoveObsoleteDataNodes(const mitk::DataNode* imageNode, const mitk::DataNode* roiNode) const { if (imageNode == nullptr || !imageNode->GetData()) { mitkThrow() << "Image is nullptr"; } const auto image = imageNode->GetData(); const mitk::BaseData* mask = nullptr; if (roiNode != nullptr) { mask = roiNode->GetData(); } auto lastResult = this->GetLatestResult(imageNode, roiNode, false, false); auto rulePredicate = mitk::ImageStatisticsContainerManager::GetStatisticsPredicateForSources(image, mask); auto notLatestPredicate = mitk::NodePredicateFunction::New([lastResult](const mitk::DataNode* node) { return node != lastResult; }); auto binPredicate = mitk::NodePredicateDataProperty::New(mitk::STATS_HISTOGRAM_BIN_PROPERTY_NAME.c_str(), mitk::UIntProperty::New(m_HistogramNBins)); auto zeroPredicate = mitk::NodePredicateDataProperty::New(mitk::STATS_IGNORE_ZERO_VOXEL_PROPERTY_NAME.c_str(), mitk::BoolProperty::New(m_IgnoreZeroValueVoxel)); mitk::NodePredicateBase::ConstPointer predicate = mitk::NodePredicateAnd::New(rulePredicate, notLatestPredicate).GetPointer(); predicate = mitk::NodePredicateAnd::New(predicate, binPredicate, zeroPredicate).GetPointer(); auto storage = m_Storage.Lock(); if (storage != nullptr) { std::lock_guard mutexguard(m_DataMutex); auto oldStatisticContainerNodes = storage->GetSubset(predicate); storage->Remove(oldStatisticContainerNodes); } } mitk::DataNode::Pointer QmitkImageStatisticsDataGenerator::PrepareResultForStorage(const std::string& /*label*/, mitk::BaseData* result, const QmitkDataGenerationJobBase* job) const { auto statsJob = dynamic_cast(job); if (statsJob != nullptr) { auto resultNode = mitk::DataNode::New(); resultNode->SetProperty("helper object", mitk::BoolProperty::New(true)); resultNode->SetVisibility(false); resultNode->SetData(result); - - const mitk::BaseData* roi = statsJob->GetMaskImage(); - if (roi == nullptr) - { - roi = statsJob->GetPlanarFigure(); - } - resultNode->SetName(this->GenerateStatisticsNodeName(statsJob->GetStatisticsImage(), roi)); + + resultNode->SetName(this->GenerateStatisticsNodeName(statsJob->GetStatisticsImage(), statsJob->GetMaskData())); return resultNode; } return nullptr; } std::string QmitkImageStatisticsDataGenerator::GenerateStatisticsNodeName(const mitk::Image* image, const mitk::BaseData* roi) const { std::stringstream statisticsNodeName; statisticsNodeName << "statistics_bins-" << m_HistogramNBins <<"_"; if (m_IgnoreZeroValueVoxel) { statisticsNodeName << "noZeros_"; } if (image == nullptr) { mitkThrow() << "Image is nullptr"; } statisticsNodeName << image->GetUID(); if (roi != nullptr) { statisticsNodeName << "_" + roi->GetUID(); } return statisticsNodeName.str(); } diff --git a/Modules/ImageStatisticsUI/Qmitk/QmitkImageStatisticsTreeModel.cpp b/Modules/ImageStatisticsUI/Qmitk/QmitkImageStatisticsTreeModel.cpp index 30e45b52c5..25a2d4b05b 100644 --- a/Modules/ImageStatisticsUI/Qmitk/QmitkImageStatisticsTreeModel.cpp +++ b/Modules/ImageStatisticsUI/Qmitk/QmitkImageStatisticsTreeModel.cpp @@ -1,521 +1,516 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #include "QmitkImageStatisticsTreeModel.h" #include "QmitkImageStatisticsTreeItem.h" #include "mitkImageStatisticsContainerManager.h" #include "mitkProportionalTimeGeometry.h" #include "mitkStatisticsToImageRelationRule.h" #include "mitkStatisticsToMaskRelationRule.h" #include "QmitkStyleManager.h" QmitkImageStatisticsTreeModel::QmitkImageStatisticsTreeModel(QObject *parent) : QmitkAbstractDataStorageModel(parent) { m_RootItem = std::make_unique(); } QmitkImageStatisticsTreeModel ::~QmitkImageStatisticsTreeModel() { // set data storage to nullptr so that the event listener gets removed this->SetDataStorage(nullptr); }; void QmitkImageStatisticsTreeModel::DataStorageChanged() { emit beginResetModel(); UpdateByDataStorage(); emit endResetModel(); emit modelChanged(); } void QmitkImageStatisticsTreeModel::NodePredicateChanged() { emit beginResetModel(); UpdateByDataStorage(); emit endResetModel(); emit modelChanged(); } int QmitkImageStatisticsTreeModel::columnCount(const QModelIndex& /*parent*/) const { int columns = m_StatisticNames.size() + 1; return columns; } int QmitkImageStatisticsTreeModel::rowCount(const QModelIndex &parent) const { QmitkImageStatisticsTreeItem *parentItem; if (parent.column() > 0) return 0; if (!parent.isValid()) parentItem = m_RootItem.get(); else parentItem = static_cast(parent.internalPointer()); return parentItem->childCount(); } QVariant QmitkImageStatisticsTreeModel::data(const QModelIndex &index, int role) const { if (!index.isValid()) return QVariant(); QmitkImageStatisticsTreeItem* item = static_cast(index.internalPointer()); if (role == Qt::DisplayRole) { return item->data(index.column()); } else if (role == Qt::DecorationRole && index.column() == 0) { if (item->isWIP() && item->childCount() == 0) return QVariant(QmitkStyleManager::ThemeIcon(QStringLiteral(":/Qmitk/hourglass-half-solid.svg"))); else if (!item->isWIP()) { auto label = item->GetLabelInstance(); if (label.IsNotNull()) { QPixmap pixmap(QSize(20,20)); QColor color(label->GetColor().GetRed() * 255, label->GetColor().GetGreen() * 255, label->GetColor().GetBlue() * 255); pixmap.fill(color); return QVariant(QIcon(pixmap)); } } } return QVariant(); } QModelIndex QmitkImageStatisticsTreeModel::index(int row, int column, const QModelIndex &parent) const { if (!hasIndex(row, column, parent)) return QModelIndex(); QmitkImageStatisticsTreeItem *parentItem; if (!parent.isValid()) parentItem = m_RootItem.get(); else parentItem = static_cast(parent.internalPointer()); QmitkImageStatisticsTreeItem *childItem = parentItem->child(row); if (childItem) return createIndex(row, column, childItem); else return QModelIndex(); } QModelIndex QmitkImageStatisticsTreeModel::parent(const QModelIndex &child) const { if (!child.isValid()) return QModelIndex(); QmitkImageStatisticsTreeItem *childItem = static_cast(child.internalPointer()); QmitkImageStatisticsTreeItem *parentItem = childItem->parentItem(); if (parentItem == m_RootItem.get()) return QModelIndex(); return createIndex(parentItem->row(), 0, parentItem); } Qt::ItemFlags QmitkImageStatisticsTreeModel::flags(const QModelIndex &index) const { if (!index.isValid()) return {}; return QAbstractItemModel::flags(index); } QVariant QmitkImageStatisticsTreeModel::headerData(int section, Qt::Orientation orientation, int role) const { if ((Qt::DisplayRole == role) && (Qt::Horizontal == orientation)) { if (section == 0) { return m_HeaderFirstColumn; } else { return QVariant(m_StatisticNames.at(section - 1).c_str()); } } return QVariant(); } void QmitkImageStatisticsTreeModel::SetImageNodes(const std::vector &nodes) { std::vector> tempNodes; for (const auto &node : nodes) { auto data = node->GetData(); if (data) { auto timeSteps = data->GetTimeSteps(); for (unsigned int i = 0; i < timeSteps; i++) { tempNodes.push_back(std::make_pair(node, i)); } } } emit beginResetModel(); m_TimeStepResolvedImageNodes = std::move(tempNodes); m_ImageNodes = nodes; UpdateByDataStorage(); emit endResetModel(); emit modelChanged(); } void QmitkImageStatisticsTreeModel::SetMaskNodes(const std::vector &nodes) { std::vector> tempNodes; for (const auto &node : nodes) { auto data = node->GetData(); if (data) { auto timeSteps = data->GetTimeSteps(); // special case: apply one mask to each time step of an 4D image if (timeSteps == 1 && m_TimeStepResolvedImageNodes.size() > 1) { timeSteps = m_TimeStepResolvedImageNodes.size(); } for (unsigned int i = 0; i < timeSteps; i++) { tempNodes.push_back(std::make_pair(node, i)); } } } emit beginResetModel(); m_TimeStepResolvedMaskNodes = std::move(tempNodes); m_MaskNodes = nodes; UpdateByDataStorage(); emit endResetModel(); emit modelChanged(); } void QmitkImageStatisticsTreeModel::Clear() { emit beginResetModel(); m_Statistics.clear(); m_ImageNodes.clear(); m_TimeStepResolvedImageNodes.clear(); m_MaskNodes.clear(); m_StatisticNames.clear(); emit endResetModel(); emit modelChanged(); } void QmitkImageStatisticsTreeModel::SetIgnoreZeroValueVoxel(bool _arg) { if (m_IgnoreZeroValueVoxel != _arg) { emit beginResetModel(); m_IgnoreZeroValueVoxel = _arg; UpdateByDataStorage(); emit endResetModel(); emit modelChanged(); } } bool QmitkImageStatisticsTreeModel::GetIgnoreZeroValueVoxel() const { return this->m_IgnoreZeroValueVoxel; } void QmitkImageStatisticsTreeModel::SetHistogramNBins(unsigned int nbins) { if (m_HistogramNBins != nbins) { emit beginResetModel(); m_HistogramNBins = nbins; UpdateByDataStorage(); emit endResetModel(); emit modelChanged(); } } unsigned int QmitkImageStatisticsTreeModel::GetHistogramNBins() const { return this->m_HistogramNBins; } void QmitkImageStatisticsTreeModel::UpdateByDataStorage() { StatisticsContainerVector newStatistics; auto datamanager = m_DataStorage.Lock(); if (datamanager.IsNotNull()) { for (const auto &image : m_ImageNodes) { if (m_MaskNodes.empty()) { auto stats = mitk::ImageStatisticsContainerManager::GetImageStatistics(datamanager, image->GetData(), nullptr, m_IgnoreZeroValueVoxel, m_HistogramNBins, true, false); if (stats.IsNotNull()) { newStatistics.emplace_back(stats); } } else { for (const auto &mask : m_MaskNodes) { auto stats = mitk::ImageStatisticsContainerManager::GetImageStatistics(datamanager, image->GetData(), mask->GetData(), m_IgnoreZeroValueVoxel, m_HistogramNBins, true, false); if (stats.IsNotNull()) { newStatistics.emplace_back(stats); } } } } if (!newStatistics.empty()) { emit dataAvailable(); } } { std::lock_guard locked(m_Mutex); m_Statistics = newStatistics; m_StatisticNames = mitk::GetAllStatisticNames(m_Statistics); BuildHierarchicalModel(); m_BuildTime.Modified(); } } void AddTimeStepTreeItems(const mitk::ImageStatisticsContainer* statistic, const mitk::DataNode* imageNode, const mitk::DataNode* maskNode, mitk::ImageStatisticsContainer::LabelValueType labelValue, const std::vector& statisticNames, bool isWIP, QmitkImageStatisticsTreeItem* parentItem, bool& hasMultipleTimesteps) { // 4. hierarchy level: time steps (optional, only if >1 time step) if (statistic->GetTimeSteps() > 1) { for (unsigned int i = 0; i < statistic->GetTimeSteps(); i++) { QString timeStepLabel = "[" + QString::number(i) + "] " + QString::number(statistic->GetTimeGeometry()->TimeStepToTimePoint(i)) + " ms"; if (statistic->StatisticsExist(labelValue, i)) { auto statisticsItem = new QmitkImageStatisticsTreeItem( statistic->GetStatistics(labelValue,i), statisticNames, timeStepLabel, isWIP, parentItem, imageNode, maskNode); parentItem->appendChild(statisticsItem); } else { auto statisticsItem = new QmitkImageStatisticsTreeItem(statisticNames, QStringLiteral("N/A"), isWIP, parentItem, imageNode, maskNode); } } } hasMultipleTimesteps = hasMultipleTimesteps || (statistic->GetTimeSteps() > 1); } void AddLabelTreeItems(const mitk::ImageStatisticsContainer* statistic, const mitk::DataNode* imageNode, const mitk::DataNode* maskNode, mitk::ImageStatisticsContainer::LabelValueVectorType labelValues, const std::vector& statisticNames, bool isWIP, QmitkImageStatisticsTreeItem* parentItem, bool& hasMultipleTimesteps) { // 3. hierarchy level: labels (optional, only if labels >1) for (const auto labelValue : labelValues) { if (labelValue != mitk::ImageStatisticsContainer::NO_MASK_LABEL_VALUE) { //currently we only show statistics of the labeled pixel if a mask is provided QString labelLabel = QStringLiteral("unnamed label"); const auto multiLabelSeg = dynamic_cast(maskNode->GetData()); const mitk::Label* labelInstance = nullptr; if (nullptr != multiLabelSeg) { labelInstance = multiLabelSeg->GetLabel(labelValue); labelLabel = QString::fromStdString(labelInstance->GetName() + " [" + labelInstance->GetTrackingID() + "]"); } QmitkImageStatisticsTreeItem* labelItem = nullptr; if (statistic->GetTimeSteps() == 1) { // add statistical values directly in this hierarchy level auto statisticsObject = statistic->GetStatistics(labelValue, 0); labelItem = new QmitkImageStatisticsTreeItem(statisticsObject, statisticNames, labelLabel, isWIP, parentItem, imageNode, maskNode, labelInstance); } else { labelItem = new QmitkImageStatisticsTreeItem(statisticNames, labelLabel, isWIP, parentItem, imageNode, maskNode, labelInstance); AddTimeStepTreeItems(statistic, imageNode, maskNode, labelValue, statisticNames, isWIP, labelItem, hasMultipleTimesteps); } parentItem->appendChild(labelItem); } } } void QmitkImageStatisticsTreeModel::BuildHierarchicalModel() { // reset old model m_RootItem.reset(new QmitkImageStatisticsTreeItem()); bool hasMask = false; bool hasMultipleTimesteps = false; std::map dataNodeToTreeItem; for (const auto &statistic : m_Statistics) { bool isWIP = statistic->IsWIP(); // get the connected image data node/mask data node auto imageRule = mitk::StatisticsToImageRelationRule::New(); auto imageOfStatisticsPredicate = imageRule->GetDestinationsDetector(statistic); auto imageFinding = std::find_if(m_ImageNodes.begin(), m_ImageNodes.end(), [&imageOfStatisticsPredicate](const mitk::DataNode::ConstPointer& testNode) { return imageOfStatisticsPredicate->CheckNode(testNode); }); auto maskRule = mitk::StatisticsToMaskRelationRule::New(); auto maskOfStatisticsPredicate = maskRule->GetDestinationsDetector(statistic); auto maskFinding = std::find_if(m_MaskNodes.begin(), m_MaskNodes.end(), [&maskOfStatisticsPredicate](const mitk::DataNode::ConstPointer& testNode) { return maskOfStatisticsPredicate->CheckNode(testNode); }); if (imageFinding == m_ImageNodes.end()) { mitkThrow() << "no image found connected to statistic" << statistic << " Aborting."; } auto& image = *imageFinding; // image: 1. hierarchy level QmitkImageStatisticsTreeItem *imageItem = nullptr; auto search = dataNodeToTreeItem.find(image); if (search != dataNodeToTreeItem.end()) { // the tree item was created previously imageItem = search->second; } else { QString imageLabel = QString::fromStdString(image->GetName()); if (statistic->GetTimeSteps() == 1 && maskFinding == m_MaskNodes.end()) { - bool noZero = statistic->IgnoresZeroVoxel(); - //we have to check if statistics are calculated with no zero voxels, because then - //we have a label/mask (no zero mask), even if no mask is officially defined. - auto labelValue = (!isWIP && noZero) ? statistic->GetExistingLabelValues(true).front() : mitk::ImageStatisticsContainer::NO_MASK_LABEL_VALUE; + auto labelValue = isWIP ? mitk::ImageStatisticsContainer::NO_MASK_LABEL_VALUE : statistic->GetExistingLabelValues().front(); auto statisticsObject = isWIP ? mitk::ImageStatisticsContainer::ImageStatisticsObject() : statistic->GetStatistics(labelValue, 0); // create the final statistics tree item imageItem = new QmitkImageStatisticsTreeItem(statisticsObject, m_StatisticNames, imageLabel, isWIP, m_RootItem.get(), image); } else { imageItem = new QmitkImageStatisticsTreeItem(m_StatisticNames, imageLabel, isWIP, m_RootItem.get(), image); } m_RootItem->appendChild(imageItem); dataNodeToTreeItem.emplace(image, imageItem); } - const auto labelValues = statistic->GetExistingLabelValues(true); //currently we not support showing the statistics for unlabeled pixels if a mask exist if (maskFinding != m_MaskNodes.end()) { + const auto labelValues = statistic->GetExistingLabelValues(); //currently we not support showing the statistics for unlabeled pixels if a mask exist + // mask: 2. hierarchy level exists auto& mask = *maskFinding; QString maskLabel = QString::fromStdString(mask->GetName()); QmitkImageStatisticsTreeItem* maskItem; if (statistic->GetTimeSteps() == 1 && labelValues.size() == 1) { // add statistical values directly in this hierarchy level auto statisticsObject = isWIP ? mitk::ImageStatisticsContainer::ImageStatisticsObject() : statistic->GetStatistics(labelValues.front(), 0); maskItem = new QmitkImageStatisticsTreeItem(statisticsObject, m_StatisticNames, maskLabel, isWIP, imageItem, image, mask); } else { maskItem = new QmitkImageStatisticsTreeItem(m_StatisticNames, maskLabel, isWIP, imageItem, image, mask); } imageItem->appendChild(maskItem); hasMask = true; // 3. hierarchy level: labels (optional, only if more then one label in statistic) if (labelValues.size() > 1) { AddLabelTreeItems(statistic, image, mask, labelValues, m_StatisticNames, isWIP, maskItem, hasMultipleTimesteps); } else { mitk::Label::PixelType labelValue = isWIP ? 0 : labelValues.front(); AddTimeStepTreeItems(statistic, image, mask, labelValue, m_StatisticNames, isWIP, maskItem, hasMultipleTimesteps); } } else { //no mask -> but multi time step - bool noZero = statistic->IgnoresZeroVoxel(); - //we have to check if statistics are calculated with no zero voxels, because then - //we have a label/mask (no zero mask), even if no mask is officially defined. - auto labelValue = (!isWIP && noZero) ? statistic->GetExistingLabelValues(true).front() : mitk::ImageStatisticsContainer::NO_MASK_LABEL_VALUE; + auto labelValue = isWIP ? mitk::ImageStatisticsContainer::NO_MASK_LABEL_VALUE : statistic->GetExistingLabelValues().front(); AddTimeStepTreeItems(statistic, image, nullptr, labelValue, m_StatisticNames, isWIP, imageItem, hasMultipleTimesteps); } } QString headerString = "Images"; if (hasMask) { headerString += "/Masks"; } if (hasMultipleTimesteps) { headerString += "/Timesteps"; } m_HeaderFirstColumn = headerString; } void QmitkImageStatisticsTreeModel::NodeRemoved(const mitk::DataNode* changedNode) { bool isRelevantNode = (nullptr != dynamic_cast(changedNode->GetData())); if (isRelevantNode) { emit beginResetModel(); UpdateByDataStorage(); emit endResetModel(); emit modelChanged(); } } void QmitkImageStatisticsTreeModel::NodeAdded(const mitk::DataNode * changedNode) { bool isRelevantNode = (nullptr != dynamic_cast(changedNode->GetData())); if (isRelevantNode) { emit beginResetModel(); UpdateByDataStorage(); emit endResetModel(); emit modelChanged(); } } void QmitkImageStatisticsTreeModel::NodeChanged(const mitk::DataNode * changedNode) { bool isRelevantNode = m_ImageNodes.end() != std::find(m_ImageNodes.begin(), m_ImageNodes.end(), changedNode); isRelevantNode = isRelevantNode || (m_MaskNodes.end() != std::find(m_MaskNodes.begin(), m_MaskNodes.end(), changedNode)); isRelevantNode = isRelevantNode || (nullptr != dynamic_cast(changedNode->GetData())); if (isRelevantNode) { if (m_BuildTime.GetMTime() < changedNode->GetData()->GetMTime()) { emit beginResetModel(); UpdateByDataStorage(); emit endResetModel(); emit modelChanged(); } } } diff --git a/Plugins/org.mitk.gui.qt.cest/src/internal/QmitkCESTStatisticsView.cpp b/Plugins/org.mitk.gui.qt.cest/src/internal/QmitkCESTStatisticsView.cpp index d41c562a0b..a157d14c16 100644 --- a/Plugins/org.mitk.gui.qt.cest/src/internal/QmitkCESTStatisticsView.cpp +++ b/Plugins/org.mitk.gui.qt.cest/src/internal/QmitkCESTStatisticsView.cpp @@ -1,817 +1,814 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ // itk #include "itksys/SystemTools.hxx" #include #include // Blueberry #include #include // Qmitk #include "QmitkCESTStatisticsView.h" // Qt #include #include #include // qwt #include // mitk #include #include #include #include #include #include #include #include #include #include #include #include #include // boost #include #include // stl #include #include #include #include #include #include namespace { template void GetSortPermutation(std::vector &out, const std::vector &determiningVector, Compare compare = std::less()) { out.resize(determiningVector.size()); std::iota(out.begin(), out.end(), 0); std::sort(out.begin(), out.end(), [&](unsigned i, unsigned j) { return compare(determiningVector[i], determiningVector[j]); }); } template void ApplyPermutation(const std::vector &order, std::vector &vectorToSort) { assert(order.size() == vectorToSort.size()); std::vector tempVector(vectorToSort.size()); for (unsigned i = 0; i < vectorToSort.size(); i++) { tempVector[i] = vectorToSort[order[i]]; } vectorToSort = tempVector; } template void ApplyPermutation(const std::vector &order, std::vector ¤tVector, std::vector &... remainingVectors) { ApplyPermutation(order, currentVector); ApplyPermutation(order, remainingVectors...); } template void SortVectors(const std::vector &orderDeterminingVector, Compare comparison, std::vector &... vectorsToBeSorted) { std::vector order; GetSortPermutation(order, orderDeterminingVector, comparison); ApplyPermutation(order, vectorsToBeSorted...); } } // namespace const std::string QmitkCESTStatisticsView::VIEW_ID = "org.mitk.views.cest.statistics"; QmitkCESTStatisticsView::QmitkCESTStatisticsView(QObject * /*parent*/, const char * /*name*/) { this->m_CalculatorJob = new QmitkImageStatisticsCalculationRunnable(); m_currentSelectedPosition.Fill(0.0); m_currentSelectedTimePoint = 0.; m_CrosshairPointSet = mitk::PointSet::New(); } QmitkCESTStatisticsView::~QmitkCESTStatisticsView() { while (this->m_CalculatorJob->IsRunning()) // wait until thread has finished { itksys::SystemTools::Delay(100); } delete this->m_CalculatorJob; } void QmitkCESTStatisticsView::SetFocus() { m_Controls.threeDimToFourDimPushButton->setFocus(); } void QmitkCESTStatisticsView::CreateQtPartControl(QWidget *parent) { // create GUI widgets from the Qt Designer's .ui file m_Controls.setupUi(parent); connect( m_Controls.threeDimToFourDimPushButton, SIGNAL(clicked()), this, SLOT(OnThreeDimToFourDimPushButtonClicked())); connect((QObject *)this->m_CalculatorJob, SIGNAL(ResultsAvailable()), this, SLOT(OnThreadedStatisticsCalculationEnds()), Qt::QueuedConnection); connect((QObject *)(this->m_Controls.fixedRangeCheckBox), SIGNAL(toggled(bool)), (QObject *)this, SLOT(OnFixedRangeCheckBoxToggled(bool))); connect((QObject *)(this->m_Controls.fixedRangeLowerDoubleSpinBox), SIGNAL(editingFinished()), (QObject *)this, SLOT(OnFixedRangeDoubleSpinBoxChanged())); connect((QObject *)(this->m_Controls.fixedRangeUpperDoubleSpinBox), SIGNAL(editingFinished()), (QObject *)this, SLOT(OnFixedRangeDoubleSpinBoxChanged())); m_Controls.threeDimToFourDimPushButton->setEnabled(false); m_Controls.widget_statistics->SetDataStorage(this->GetDataStorage()); this->m_SliceChangeListener.RenderWindowPartActivated(this->GetRenderWindowPart()); connect(&m_SliceChangeListener, SIGNAL(SliceChanged()), this, SLOT(OnSliceChanged())); } void QmitkCESTStatisticsView::RenderWindowPartActivated(mitk::IRenderWindowPart *renderWindowPart) { this->m_SliceChangeListener.RenderWindowPartActivated(renderWindowPart); } void QmitkCESTStatisticsView::RenderWindowPartDeactivated(mitk::IRenderWindowPart *renderWindowPart) { this->m_SliceChangeListener.RenderWindowPartDeactivated(renderWindowPart); } void QmitkCESTStatisticsView::OnSelectionChanged(berry::IWorkbenchPart::Pointer /*source*/, const QList &nodes) { if (nodes.empty()) { std::stringstream message; message << "Please select an image."; m_Controls.labelWarning->setText(message.str().c_str()); m_Controls.labelWarning->show(); this->Clear(); return; } // iterate all selected objects bool atLeastOneWasCESTImage = false; foreach (mitk::DataNode::Pointer node, nodes) { if (node.IsNull()) { continue; } if (dynamic_cast(node->GetData()) != nullptr) { m_Controls.labelWarning->setVisible(false); bool zSpectrumSet = SetZSpectrum(dynamic_cast( node->GetData()->GetProperty(mitk::CEST_PROPERTY_NAME_OFFSETS().c_str()).GetPointer())); atLeastOneWasCESTImage = atLeastOneWasCESTImage || zSpectrumSet; if (zSpectrumSet) { m_ZImage = dynamic_cast(node->GetData()); m_Controls.widget_statistics->SetImageNodes({node.GetPointer()}); } else { m_MaskImage = dynamic_cast(node->GetData()); m_Controls.widget_statistics->SetMaskNodes({node.GetPointer()}); } } if (dynamic_cast(node->GetData()) != nullptr) { m_MaskPlanarFigure = dynamic_cast(node->GetData()); m_Controls.widget_statistics->SetMaskNodes({node.GetPointer()}); } if (dynamic_cast(node->GetData()) != nullptr) { m_PointSet = dynamic_cast(node->GetData()); } } // We only want to offer normalization or timestep copying if one object is selected if (nodes.size() == 1) { if (dynamic_cast(nodes.front()->GetData())) { m_Controls.threeDimToFourDimPushButton->setDisabled(atLeastOneWasCESTImage); } else { m_Controls.threeDimToFourDimPushButton->setEnabled(false); std::stringstream message; message << "The selected node is not an image."; m_Controls.labelWarning->setText(message.str().c_str()); m_Controls.labelWarning->show(); } this->Clear(); return; } // we always need a mask, either image or planar figure as well as an image for further processing if (nodes.size() != 2) { this->Clear(); return; } m_Controls.threeDimToFourDimPushButton->setEnabled(false); if (!atLeastOneWasCESTImage) { std::stringstream message; message << "None of the selected data nodes contains required CEST meta information"; m_Controls.labelWarning->setText(message.str().c_str()); m_Controls.labelWarning->show(); this->Clear(); return; } bool bothAreImages = (m_ZImage.GetPointer() != nullptr) && (m_MaskImage.GetPointer() != nullptr); if (bothAreImages) { bool geometriesMatch = mitk::Equal(*(m_ZImage->GetTimeGeometry()), *(m_MaskImage->GetTimeGeometry()), mitk::eps, false); if (!geometriesMatch) { std::stringstream message; message << "The selected images have different geometries."; m_Controls.labelWarning->setText(message.str().c_str()); m_Controls.labelWarning->show(); this->Clear(); return; } } if (!this->DataSanityCheck()) { this->Clear(); return; } if (m_PointSet.IsNull()) { // initialize thread and trigger it this->m_CalculatorJob->SetIgnoreZeroValueVoxel(false); - this->m_CalculatorJob->Initialize(m_ZImage.GetPointer(), m_MaskImage.GetPointer(), m_MaskPlanarFigure.GetPointer()); + const mitk::BaseData* maskData = m_MaskImage.GetPointer(); + if (m_MaskImage.IsNull()) maskData = m_MaskPlanarFigure.GetPointer(); + this->m_CalculatorJob->Initialize(m_ZImage.GetPointer(), maskData); std::stringstream message; message << "Calculating statistics..."; m_Controls.labelWarning->setText(message.str().c_str()); m_Controls.labelWarning->show(); try { // Compute statistics QThreadPool::globalInstance()->start(m_CalculatorJob); } catch (const mitk::Exception &e) { std::stringstream message; message << "" << e.GetDescription() << ""; m_Controls.labelWarning->setText(message.str().c_str()); m_Controls.labelWarning->show(); } catch (const std::runtime_error &e) { // In case of exception, print error message on GUI std::stringstream message; message << "" << e.what() << ""; m_Controls.labelWarning->setText(message.str().c_str()); m_Controls.labelWarning->show(); } catch (const std::exception &e) { MITK_ERROR << "Caught exception: " << e.what(); // In case of exception, print error message on GUI std::stringstream message; message << "Error! Unequal Dimensions of Image and Segmentation. No recompute possible "; m_Controls.labelWarning->setText(message.str().c_str()); m_Controls.labelWarning->show(); } while (this->m_CalculatorJob->IsRunning()) // wait until thread has finished { itksys::SystemTools::Delay(100); } } if (m_PointSet.IsNotNull()) { if (m_ZImage->GetDimension() == 4) { AccessFixedDimensionByItk(m_ZImage, PlotPointSet, 4); } else { MITK_WARN << "Expecting a 4D image."; } } } void QmitkCESTStatisticsView::OnThreadedStatisticsCalculationEnds() { this->m_Controls.m_DataViewWidget->SetAxisTitle(QwtPlot::Axis::xBottom, "delta w"); this->m_Controls.m_DataViewWidget->SetAxisTitle(QwtPlot::Axis::yLeft, "z"); if (this->m_CalculatorJob->GetComputationSuccessFlag()) { auto statistics = this->m_CalculatorJob->GetStatisticsData(); std::string statisticsNodeName = "CEST_statistics"; auto statisticsNode = mitk::CreateImageStatisticsNode(statistics, statisticsNodeName); auto imageRule = mitk::StatisticsToImageRelationRule::New(); imageRule->Connect(statistics, m_CalculatorJob->GetStatisticsImage()); - if (m_CalculatorJob->GetMaskImage()) + if (m_CalculatorJob->GetMaskData()) { auto maskRule = mitk::StatisticsToMaskRelationRule::New(); - maskRule->Connect(statistics, m_CalculatorJob->GetMaskImage()); - } - else if (m_CalculatorJob->GetPlanarFigure()) - { - auto planarFigureRule = mitk::StatisticsToMaskRelationRule::New(); - planarFigureRule->Connect(statistics, m_CalculatorJob->GetPlanarFigure()); + maskRule->Connect(statistics, m_CalculatorJob->GetMaskData()); } this->GetDataStorage()->Add(statisticsNode); QmitkPlotWidget::DataVector::size_type numberOfSpectra = this->m_zSpectrum.size(); QmitkPlotWidget::DataVector means(numberOfSpectra); QmitkPlotWidget::DataVector stdevs(numberOfSpectra); for (unsigned int index = 0; index < numberOfSpectra; ++index) { means[index] = statistics->GetStatistics(1,index).GetValueConverted( mitk::ImageStatisticsConstants::MEAN()); stdevs[index] = statistics->GetStatistics(1,index).GetValueConverted( mitk::ImageStatisticsConstants::STANDARDDEVIATION()); } QmitkPlotWidget::DataVector xValues = this->m_zSpectrum; RemoveMZeros(xValues, means, stdevs); ::SortVectors(xValues, std::less(), xValues, means, stdevs); unsigned int curveId = this->m_Controls.m_DataViewWidget->InsertCurve("Spectrum"); this->m_Controls.m_DataViewWidget->SetCurveData(curveId, xValues, means, stdevs, stdevs); this->m_Controls.m_DataViewWidget->SetErrorPen(curveId, QPen(Qt::blue)); QwtSymbol *blueSymbol = new QwtSymbol(QwtSymbol::Rect, QColor(Qt::blue), QColor(Qt::blue), QSize(8, 8)); this->m_Controls.m_DataViewWidget->SetCurveSymbol(curveId, blueSymbol); this->m_Controls.m_DataViewWidget->SetLegendAttribute(curveId, QwtPlotCurve::LegendShowSymbol); QwtLegend *legend = new QwtLegend(); legend->setFrameShape(QFrame::Box); legend->setFrameShadow(QFrame::Sunken); legend->setLineWidth(1); this->m_Controls.m_DataViewWidget->SetLegend(legend, QwtPlot::BottomLegend); m_Controls.m_DataViewWidget->GetPlot() ->axisScaleEngine(QwtPlot::Axis::xBottom) ->setAttributes(QwtScaleEngine::Inverted); this->m_Controls.m_DataViewWidget->Replot(); m_Controls.labelWarning->setVisible(false); m_Controls.m_StatisticsGroupBox->setEnabled(true); m_Controls.m_StatisticsGroupBox->setEnabled(true); if (this->m_Controls.fixedRangeCheckBox->isChecked()) { this->m_Controls.m_DataViewWidget->GetPlot()->setAxisAutoScale(2, false); this->m_Controls.m_DataViewWidget->GetPlot()->setAxisScale( 2, this->m_Controls.fixedRangeLowerDoubleSpinBox->value(), this->m_Controls.fixedRangeUpperDoubleSpinBox->value()); } else { this->m_Controls.m_DataViewWidget->GetPlot()->setAxisAutoScale(2, true); } } else { m_Controls.labelWarning->setText(m_CalculatorJob->GetLastErrorMessage().c_str()); m_Controls.labelWarning->setVisible(true); this->Clear(); } } void QmitkCESTStatisticsView::OnFixedRangeDoubleSpinBoxChanged() { if (this->m_Controls.fixedRangeCheckBox->isChecked()) { this->m_Controls.m_DataViewWidget->GetPlot()->setAxisAutoScale(2, false); this->m_Controls.m_DataViewWidget->GetPlot()->setAxisScale(2, this->m_Controls.fixedRangeLowerDoubleSpinBox->value(), this->m_Controls.fixedRangeUpperDoubleSpinBox->value()); } this->m_Controls.m_DataViewWidget->Replot(); } template void QmitkCESTStatisticsView::PlotPointSet(itk::Image *image) { this->m_Controls.m_DataViewWidget->SetAxisTitle(QwtPlot::Axis::xBottom, "delta w"); this->m_Controls.m_DataViewWidget->SetAxisTitle(QwtPlot::Axis::yLeft, "z"); QmitkPlotWidget::DataVector::size_type numberOfSpectra = this->m_zSpectrum.size(); mitk::PointSet::Pointer internalPointset; if (m_PointSet.IsNotNull()) { internalPointset = m_PointSet; } else { internalPointset = m_CrosshairPointSet; } if (internalPointset.IsNull()) { return; } if (!this->DataSanityCheck()) { m_Controls.labelWarning->setText("Data can not be plotted, internally inconsistent."); m_Controls.labelWarning->show(); return; } auto maxIndex = internalPointset->GetMaxId().Index(); for (std::size_t number = 0; number < maxIndex + 1; ++number) { mitk::PointSet::PointType point; if (!internalPointset->GetPointIfExists(number, &point)) { continue; } if (!this->m_ZImage->GetGeometry()->IsInside(point)) { continue; } itk::Index<3> itkIndex; this->m_ZImage->GetGeometry()->WorldToIndex(point, itkIndex); itk::Index itkIndexTime; itkIndexTime[0] = itkIndex[0]; itkIndexTime[1] = itkIndex[1]; itkIndexTime[2] = itkIndex[2]; QmitkPlotWidget::DataVector values(numberOfSpectra); for (std::size_t step = 0; step < numberOfSpectra; ++step) { if (VImageDimension == 4) { itkIndexTime[3] = step; } values[step] = image->GetPixel(itkIndexTime); } std::stringstream name; name << "Point " << number; // Qcolor enums go from 0 to 19, but 19 is transparent and 0,1 are for bitmaps // 3 is white and thus not visible QColor color(static_cast(number % 17 + 4)); QmitkPlotWidget::DataVector xValues = this->m_zSpectrum; RemoveMZeros(xValues, values); ::SortVectors(xValues, std::less(), xValues, values); unsigned int curveId = this->m_Controls.m_DataViewWidget->InsertCurve(name.str().c_str()); this->m_Controls.m_DataViewWidget->SetCurveData(curveId, xValues, values); this->m_Controls.m_DataViewWidget->SetCurvePen(curveId, QPen(color)); QwtSymbol *symbol = new QwtSymbol(QwtSymbol::Rect, color, color, QSize(8, 8)); this->m_Controls.m_DataViewWidget->SetCurveSymbol(curveId, symbol); this->m_Controls.m_DataViewWidget->SetLegendAttribute(curveId, QwtPlotCurve::LegendShowSymbol); } if (this->m_Controls.fixedRangeCheckBox->isChecked()) { this->m_Controls.m_DataViewWidget->GetPlot()->setAxisAutoScale(2, false); this->m_Controls.m_DataViewWidget->GetPlot()->setAxisScale(2, this->m_Controls.fixedRangeLowerDoubleSpinBox->value(), this->m_Controls.fixedRangeUpperDoubleSpinBox->value()); } else { this->m_Controls.m_DataViewWidget->GetPlot()->setAxisAutoScale(2, true); } QwtLegend *legend = new QwtLegend(); legend->setFrameShape(QFrame::Box); legend->setFrameShadow(QFrame::Sunken); legend->setLineWidth(1); this->m_Controls.m_DataViewWidget->SetLegend(legend, QwtPlot::BottomLegend); m_Controls.m_DataViewWidget->GetPlot() ->axisScaleEngine(QwtPlot::Axis::xBottom) ->setAttributes(QwtScaleEngine::Inverted); this->m_Controls.m_DataViewWidget->Replot(); m_Controls.labelWarning->setVisible(false); } void QmitkCESTStatisticsView::OnFixedRangeCheckBoxToggled(bool state) { this->m_Controls.fixedRangeLowerDoubleSpinBox->setEnabled(state); this->m_Controls.fixedRangeUpperDoubleSpinBox->setEnabled(state); } void QmitkCESTStatisticsView::RemoveMZeros(QmitkPlotWidget::DataVector &xValues, QmitkPlotWidget::DataVector &yValues) { QmitkPlotWidget::DataVector tempX; QmitkPlotWidget::DataVector tempY; for (std::size_t index = 0; index < xValues.size(); ++index) { if ((xValues.at(index) < -299) || (xValues.at(index)) > 299) { // do not include } else { tempX.push_back(xValues.at(index)); tempY.push_back(yValues.at(index)); } } xValues = tempX; yValues = tempY; } void QmitkCESTStatisticsView::RemoveMZeros(QmitkPlotWidget::DataVector &xValues, QmitkPlotWidget::DataVector &yValues, QmitkPlotWidget::DataVector &stdDevs) { QmitkPlotWidget::DataVector tempX; QmitkPlotWidget::DataVector tempY; QmitkPlotWidget::DataVector tempDevs; for (std::size_t index = 0; index < xValues.size(); ++index) { if ((xValues.at(index) < -299) || (xValues.at(index)) > 299) { // do not include } else { tempX.push_back(xValues.at(index)); tempY.push_back(yValues.at(index)); tempDevs.push_back(stdDevs.at(index)); } } xValues = tempX; yValues = tempY; stdDevs = tempDevs; } void QmitkCESTStatisticsView::OnThreeDimToFourDimPushButtonClicked() { QList nodes = this->GetDataManagerSelection(); if (nodes.empty()) return; mitk::DataNode *node = nodes.front(); if (!node) { // Nothing selected. Inform the user and return QMessageBox::information(nullptr, "CEST View", "Please load and select an image before starting image processing."); return; } // here we have a valid mitk::DataNode // a node itself is not very useful, we need its data item (the image) mitk::BaseData *data = node->GetData(); if (data) { // test if this data item is an image or not (could also be a surface or something totally different) mitk::Image *image = dynamic_cast(data); if (image) { if (image->GetDimension() == 4) { AccessFixedDimensionByItk(image, CopyTimesteps, 4); } this->Clear(); } } } template void QmitkCESTStatisticsView::CopyTimesteps(itk::Image *image) { typedef itk::Image ImageType; // typedef itk::PasteImageFilter PasteImageFilterType; unsigned int numberOfTimesteps = image->GetLargestPossibleRegion().GetSize(3); typename ImageType::RegionType sourceRegion = image->GetLargestPossibleRegion(); sourceRegion.SetSize(3, 1); typename ImageType::RegionType targetRegion = image->GetLargestPossibleRegion(); targetRegion.SetSize(3, 1); for (unsigned int timestep = 1; timestep < numberOfTimesteps; ++timestep) { targetRegion.SetIndex(3, timestep); itk::ImageRegionConstIterator sourceIterator(image, sourceRegion); itk::ImageRegionIterator targetIterator(image, targetRegion); while (!sourceIterator.IsAtEnd()) { targetIterator.Set(sourceIterator.Get()); ++sourceIterator; ++targetIterator; } } } bool QmitkCESTStatisticsView::SetZSpectrum(mitk::StringProperty *zSpectrumProperty) { if (zSpectrumProperty == nullptr) { return false; } mitk::LocaleSwitch localeSwitch("C"); std::string zSpectrumString = zSpectrumProperty->GetValueAsString(); std::istringstream iss(zSpectrumString); std::vector zSpectra; std::copy( std::istream_iterator(iss), std::istream_iterator(), std::back_inserter(zSpectra)); m_zSpectrum.clear(); m_zSpectrum.resize(0); for (auto const &spectrumString : zSpectra) { m_zSpectrum.push_back(std::stod(spectrumString)); } return (m_zSpectrum.size() > 0); } bool QmitkCESTStatisticsView::DataSanityCheck() { QmitkPlotWidget::DataVector::size_type numberOfSpectra = m_zSpectrum.size(); // if we do not have a spectrum, the data can not be processed if (numberOfSpectra == 0) { return false; } // if we do not have CEST image data, the data can not be processed if (m_ZImage.IsNull()) { return false; } // if the CEST image data and the meta information do not match, the data can not be processed if (numberOfSpectra != m_ZImage->GetTimeSteps()) { MITK_INFO << "CEST meta information and number of volumes does not match."; return false; } // if we have neither a mask image, a point set nor a mask planar figure, we can not do statistics // statistics on the whole image would not make sense if (m_MaskImage.IsNull() && m_MaskPlanarFigure.IsNull() && m_PointSet.IsNull() && m_CrosshairPointSet->IsEmpty()) { return false; } // if we have a mask image and a mask planar figure, we can not do statistics // we do not know which one to use if (m_MaskImage.IsNotNull() && m_MaskPlanarFigure.IsNotNull()) { return false; } return true; } void QmitkCESTStatisticsView::Clear() { this->m_zSpectrum.clear(); this->m_zSpectrum.resize(0); this->m_ZImage = nullptr; this->m_MaskImage = nullptr; this->m_MaskPlanarFigure = nullptr; this->m_PointSet = nullptr; this->m_Controls.m_DataViewWidget->Clear(); this->m_Controls.m_StatisticsGroupBox->setEnabled(false); this->m_Controls.widget_statistics->SetImageNodes({}); this->m_Controls.widget_statistics->SetMaskNodes({}); } void QmitkCESTStatisticsView::OnSliceChanged() { auto* renderWindowPart = this->GetRenderWindowPart(mitk::WorkbenchUtil::OPEN); auto currentSelectedPosition = renderWindowPart->GetSelectedPosition(nullptr); auto currentSelectedTimePoint = renderWindowPart->GetSelectedTimePoint(); if (m_currentSelectedPosition != currentSelectedPosition || m_currentSelectedTimePoint != currentSelectedTimePoint) //|| m_selectedNodeTime > m_currentPositionTime) { // the current position has been changed or the selected node has been changed since the last position validation -> // check position m_currentSelectedPosition = currentSelectedPosition; m_currentSelectedTimePoint = currentSelectedTimePoint; m_currentPositionTime.Modified(); m_CrosshairPointSet->Clear(); m_CrosshairPointSet->SetPoint(0, m_currentSelectedPosition); QList nodes = this->GetDataManagerSelection(); if (nodes.empty() || nodes.size() > 1) return; mitk::DataNode *node = nodes.front(); if (!node) { return; } if (dynamic_cast(node->GetData()) != nullptr) { m_Controls.labelWarning->setVisible(false); bool zSpectrumSet = SetZSpectrum(dynamic_cast( node->GetData()->GetProperty(mitk::CEST_PROPERTY_NAME_OFFSETS().c_str()).GetPointer())); if (zSpectrumSet) { m_ZImage = dynamic_cast(node->GetData()); } else { return; } } else { return; } this->m_Controls.m_DataViewWidget->Clear(); AccessFixedDimensionByItk(m_ZImage, PlotPointSet, 4); } } diff --git a/Plugins/org.mitk.gui.qt.cest/src/internal/QmitkCESTStatisticsView.h b/Plugins/org.mitk.gui.qt.cest/src/internal/QmitkCESTStatisticsView.h index 13a2bf797c..32ea468492 100644 --- a/Plugins/org.mitk.gui.qt.cest/src/internal/QmitkCESTStatisticsView.h +++ b/Plugins/org.mitk.gui.qt.cest/src/internal/QmitkCESTStatisticsView.h @@ -1,130 +1,131 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #ifndef QmitkCESTStatisticsView_h #define QmitkCESTStatisticsView_h #include #include #include #include "ui_QmitkCESTStatisticsViewControls.h" #include +#include #include #include /** \brief QmitkCESTStatisticsView \warning Basic statistics view for CEST data. \sa QmitkAbstractView \ingroup ${plugin_target}_internal */ class QmitkCESTStatisticsView : public QmitkAbstractView, public mitk::IRenderWindowPartListener { Q_OBJECT public: static const std::string VIEW_ID; /*! \brief default constructor */ QmitkCESTStatisticsView(QObject *parent = nullptr, const char *name = nullptr); /*! \brief default destructor */ ~QmitkCESTStatisticsView() override; protected slots: /// \brief Called when the user clicks the GUI button void OnThreeDimToFourDimPushButtonClicked(); /// \brief takes care of processing the computed data void OnThreadedStatisticsCalculationEnds(); /// \brief Toggle whether or not the plot uses a fixed x range void OnFixedRangeCheckBoxToggled(bool state); /// \brief Adapt axis scale when manual ranges are set void OnFixedRangeDoubleSpinBoxChanged(); /// \brief What to do if the crosshair moves void OnSliceChanged(); protected: void CreateQtPartControl(QWidget *parent) override; void SetFocus() override; void RenderWindowPartActivated(mitk::IRenderWindowPart* renderWindowPart) override; void RenderWindowPartDeactivated(mitk::IRenderWindowPart* renderWindowPart) override; void OnSelectionChanged( berry::IWorkbenchPart::Pointer source, const QList& nodes ) override; /// parse string and set data vector returns true if successful bool SetZSpectrum(mitk::StringProperty* zSpectrumProperty); /** Checks whether the currently set data appears reasonable */ bool DataSanityCheck(); /** Fills the plot based on a point set * * This will only use the first timestep */ template void PlotPointSet(itk::Image* image); /** Deletes all data */ void Clear(); /** Remove MZeros * * Will remove the data for the M0 images from the given input */ void RemoveMZeros(QmitkPlotWidget::DataVector& xValues, QmitkPlotWidget::DataVector& yValues); void RemoveMZeros(QmitkPlotWidget::DataVector& xValues, QmitkPlotWidget::DataVector& yValues, QmitkPlotWidget::DataVector& stdDevs); /** Copies the first timestep of a segmentation to all others */ template void CopyTimesteps(itk::Image* image); Ui::QmitkCESTStatisticsViewControls m_Controls; QmitkImageStatisticsCalculationRunnable* m_CalculatorJob; QmitkPlotWidget::DataVector m_zSpectrum; mitk::Image::Pointer m_ZImage; mitk::Image::Pointer m_MaskImage; mitk::PlanarFigure::Pointer m_MaskPlanarFigure; mitk::PointSet::Pointer m_PointSet; mitk::PointSet::Pointer m_CrosshairPointSet; QmitkSliceNavigationListener m_SliceChangeListener; itk::TimeStamp m_selectedNodeTime; itk::TimeStamp m_currentPositionTime; /** @brief currently valid selected position in the inspector*/ mitk::Point3D m_currentSelectedPosition; mitk::TimePointType m_currentSelectedTimePoint; }; #endif diff --git a/Plugins/org.mitk.gui.qt.measurementtoolbox/src/internal/QmitkImageStatisticsView.cpp b/Plugins/org.mitk.gui.qt.measurementtoolbox/src/internal/QmitkImageStatisticsView.cpp index a59779c8c5..ec0aea40e9 100644 --- a/Plugins/org.mitk.gui.qt.measurementtoolbox/src/internal/QmitkImageStatisticsView.cpp +++ b/Plugins/org.mitk.gui.qt.measurementtoolbox/src/internal/QmitkImageStatisticsView.cpp @@ -1,494 +1,494 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #include "QmitkImageStatisticsView.h" #include // berry includes #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include "mitkPlanarFigureMaskGenerator.h" #include "QmitkImageStatisticsDataGenerator.h" #include "mitkImageStatisticsContainerManager.h" #include const std::string QmitkImageStatisticsView::VIEW_ID = "org.mitk.views.imagestatistics"; namespace { bool CheckPlanarFigureMatchesGeometry(const mitk::PlanarFigure* planarFigure, const mitk::BaseGeometry* imageGeometry) { if (!planarFigure || !imageGeometry) return false; if (!mitk::PlanarFigureMaskGenerator::CheckPlanarFigureIsNotTilted(planarFigure->GetPlaneGeometry(), imageGeometry)) return false; // The rest from here on is only needed until T30279 has been solved. if (planarFigure->IsClosed()) return true; const auto numControlPoints = planarFigure->GetNumberOfControlPoints(); for (unsigned int i = 0; i < numControlPoints; ++i) { if (!imageGeometry->IsInside(planarFigure->GetWorldControlPoint(i))) return false; } return true; } } // unnamed namespace QmitkImageStatisticsView::~QmitkImageStatisticsView() { } void QmitkImageStatisticsView::CreateQtPartControl(QWidget *parent) { m_Controls.setupUi(parent); m_Controls.widget_intensityProfile->SetTheme(GetColorTheme()); m_Controls.groupBox_histogram->setVisible(false); m_Controls.groupBox_intensityProfile->setVisible(false); m_Controls.label_currentlyComputingStatistics->setVisible(false); m_Controls.sliderWidget_histogram->setPrefix("Time: "); m_Controls.sliderWidget_histogram->setDecimals(0); m_Controls.sliderWidget_histogram->setVisible(false); m_Controls.sliderWidget_intensityProfile->setPrefix("Time: "); m_Controls.sliderWidget_intensityProfile->setDecimals(0); m_Controls.sliderWidget_intensityProfile->setVisible(false); ResetGUI(); m_DataGenerator = new QmitkImageStatisticsDataGenerator(parent); m_DataGenerator->SetDataStorage(this->GetDataStorage()); m_DataGenerator->SetAutoUpdate(true); m_Controls.widget_statistics->SetDataStorage(this->GetDataStorage()); m_Controls.imageNodesSelector->SetDataStorage(this->GetDataStorage()); m_Controls.imageNodesSelector->SetNodePredicate(mitk::GetImageStatisticsImagePredicate()); m_Controls.imageNodesSelector->SetSelectionCheckFunction(this->CheckForSameGeometry()); m_Controls.imageNodesSelector->SetSelectionIsOptional(false); m_Controls.imageNodesSelector->SetInvalidInfo(QStringLiteral("Please select images for statistics")); m_Controls.imageNodesSelector->SetPopUpTitel(QStringLiteral("Select input images")); m_Controls.roiNodesSelector->SetPopUpHint(QStringLiteral("You may select multiple images for the statistics computation. But all selected images must have the same geometry.")); m_Controls.roiNodesSelector->SetDataStorage(this->GetDataStorage()); m_Controls.roiNodesSelector->SetNodePredicate(this->GenerateROIPredicate()); m_Controls.roiNodesSelector->SetSelectionIsOptional(true); m_Controls.roiNodesSelector->SetEmptyInfo(QStringLiteral("Please select ROIs")); m_Controls.roiNodesSelector->SetPopUpTitel(QStringLiteral("Select ROIs for statistics computation")); m_Controls.roiNodesSelector->SetPopUpHint(QStringLiteral("You may select ROIs (e.g. planar figures, segmentations) that should be used for the statistics computation. The statistics will only computed for the image parts defined by the ROIs.")); CreateConnections(); this->m_TimePointChangeListener.RenderWindowPartActivated(this->GetRenderWindowPart()); connect(&m_TimePointChangeListener, &QmitkSliceNavigationListener::SelectedTimePointChanged, this, & QmitkImageStatisticsView::OnSelectedTimePointChanged); } void QmitkImageStatisticsView::RenderWindowPartActivated(mitk::IRenderWindowPart* renderWindowPart) { this->m_TimePointChangeListener.RenderWindowPartActivated(renderWindowPart); } void QmitkImageStatisticsView::RenderWindowPartDeactivated(mitk::IRenderWindowPart* renderWindowPart) { this->m_TimePointChangeListener.RenderWindowPartDeactivated(renderWindowPart); } void QmitkImageStatisticsView::CreateConnections() { connect(m_Controls.widget_statistics, &QmitkImageStatisticsWidget::IgnoreZeroValuedVoxelStateChanged, this, &QmitkImageStatisticsView::OnIgnoreZeroValuedVoxelStateChanged); connect(m_Controls.buttonSelection, &QAbstractButton::clicked, this, &QmitkImageStatisticsView::OnButtonSelectionPressed); connect(m_Controls.widget_histogram, &QmitkHistogramVisualizationWidget::RequestHistogramUpdate, this, &QmitkImageStatisticsView::OnRequestHistogramUpdate); connect(m_DataGenerator, &QmitkImageStatisticsDataGenerator::DataGenerationStarted, this, &QmitkImageStatisticsView::OnGenerationStarted); connect(m_DataGenerator, &QmitkImageStatisticsDataGenerator::GenerationFinished, this, &QmitkImageStatisticsView::OnGenerationFinished); connect(m_DataGenerator, &QmitkImageStatisticsDataGenerator::JobError, this, &QmitkImageStatisticsView::OnJobError); connect(m_Controls.imageNodesSelector, &QmitkAbstractNodeSelectionWidget::CurrentSelectionChanged, this, &QmitkImageStatisticsView::OnImageSelectionChanged); connect(m_Controls.roiNodesSelector, &QmitkAbstractNodeSelectionWidget::CurrentSelectionChanged, this, &QmitkImageStatisticsView::OnROISelectionChanged); connect(m_Controls.sliderWidget_intensityProfile, &ctkSliderWidget::valueChanged, this, &QmitkImageStatisticsView::UpdateIntensityProfile); } void QmitkImageStatisticsView::UpdateIntensityProfile() { m_Controls.groupBox_intensityProfile->setVisible(false); const auto selectedImageNodes = m_Controls.imageNodesSelector->GetSelectedNodes(); const auto selectedROINodes = m_Controls.roiNodesSelector->GetSelectedNodes(); if (selectedImageNodes.size()==1 && selectedROINodes.size()==1) { //only supported for one image and roi currently auto image = dynamic_cast(selectedImageNodes.front()->GetData()); auto maskPlanarFigure = dynamic_cast(selectedROINodes.front()->GetData()); if (maskPlanarFigure != nullptr) { if (!maskPlanarFigure->IsClosed()) { mitk::Image::Pointer inputImage; if (image->GetDimension() == 4) { m_Controls.sliderWidget_intensityProfile->setVisible(true); unsigned int maxTimestep = image->GetTimeSteps(); m_Controls.sliderWidget_intensityProfile->setMaximum(maxTimestep - 1); // Intensity profile can only be calculated on 3D, so extract if 4D mitk::ImageTimeSelector::Pointer timeSelector = mitk::ImageTimeSelector::New(); int currentTimestep = static_cast(m_Controls.sliderWidget_intensityProfile->value()); timeSelector->SetInput(image); timeSelector->SetTimeNr(currentTimestep); timeSelector->Update(); inputImage = timeSelector->GetOutput(); } else { m_Controls.sliderWidget_intensityProfile->setVisible(false); inputImage = image; } auto intensityProfile = mitk::ComputeIntensityProfile(inputImage, maskPlanarFigure); m_Controls.groupBox_intensityProfile->setVisible(true); m_Controls.widget_intensityProfile->Reset(); m_Controls.widget_intensityProfile->SetIntensityProfile(intensityProfile.GetPointer(), "Intensity Profile of " + selectedImageNodes.front()->GetName()); } } } } void QmitkImageStatisticsView::UpdateHistogramWidget() { m_Controls.groupBox_histogram->setVisible(false); const auto selectedImageNodes = m_Controls.imageNodesSelector->GetSelectedNodes(); const auto selectedMaskNodes = m_Controls.roiNodesSelector->GetSelectedNodes(); if (selectedImageNodes.size() == 1 && selectedMaskNodes.size()<=1) { //currently only supported for one image and roi due to histogram widget limitations. auto imageNode = selectedImageNodes.front(); const mitk::DataNode* roiNode = nullptr; const mitk::PlanarFigure* planarFigure = nullptr; if (!selectedMaskNodes.empty()) { roiNode = selectedMaskNodes.front(); planarFigure = dynamic_cast(roiNode->GetData()); } if ((planarFigure == nullptr || planarFigure->IsClosed()) && imageNode->GetData()->GetTimeGeometry()->IsValidTimePoint(m_TimePointChangeListener.GetCurrentSelectedTimePoint())) { //if a planar figure is not closed, we show the intensity profile instead of the histogram. auto statisticsNode = m_DataGenerator->GetLatestResult(imageNode, roiNode, true); if (statisticsNode.IsNotNull()) { auto statistics = dynamic_cast(statisticsNode->GetData()); if (statistics && !statistics->IsWIP()) { //currently only supports rois with one label due to histogram widget limitations. - auto labelValues = statistics->GetExistingLabelValues(true); - if (labelValues.size() == 1 || (!(statistics->IgnoresZeroVoxel()) && labelValues.empty())) + auto labelValues = statistics->GetExistingLabelValues(); + if (labelValues.size() == 1) { auto labelValue = labelValues.empty() ? mitk::ImageStatisticsContainer::NO_MASK_LABEL_VALUE : labelValues.front(); const auto timeStep = imageNode->GetData()->GetTimeGeometry()->TimePointToTimeStep(m_TimePointChangeListener.GetCurrentSelectedTimePoint()); if (statistics->StatisticsExist(labelValue, timeStep)) { std::stringstream label; label << imageNode->GetName(); if (imageNode->GetData()->GetTimeSteps() > 1) { label << "[" << timeStep << "]"; } if (roiNode) { label << " with " << roiNode->GetName(); } //Hardcoded labels are currently needed because the current histogram widget (and ChartWidget) //do not allow correct removal or sound update/insertion of several charts. //only thing that works for now is always to update/overwrite the same data label //This is a quick fix for T28223 and T28221 m_Controls.widget_histogram->SetHistogram(statistics->GetHistogram(labelValue, timeStep), "histogram"); m_Controls.groupBox_histogram->setVisible(true); } } } } } } } QmitkChartWidget::ColorTheme QmitkImageStatisticsView::GetColorTheme() const { ctkPluginContext *context = berry::WorkbenchPlugin::GetDefault()->GetPluginContext(); ctkServiceReference styleManagerRef = context->getServiceReference(); if (styleManagerRef) { auto styleManager = context->getService(styleManagerRef); if (styleManager->GetStyle().name == "Dark") { return QmitkChartWidget::ColorTheme::darkstyle; } else { return QmitkChartWidget::ColorTheme::lightstyle; } } return QmitkChartWidget::ColorTheme::darkstyle; } void QmitkImageStatisticsView::ResetGUI() { m_Controls.widget_statistics->Reset(); m_Controls.widget_statistics->setEnabled(false); m_Controls.widget_histogram->Reset(); m_Controls.widget_histogram->setEnabled(false); m_Controls.widget_histogram->SetTheme(GetColorTheme()); } void QmitkImageStatisticsView::OnGenerationStarted(const mitk::DataNode* /*imageNode*/, const mitk::DataNode* /*roiNode*/, const QmitkDataGenerationJobBase* /*job*/) { m_Controls.label_currentlyComputingStatistics->setVisible(true); } void QmitkImageStatisticsView::OnGenerationFinished() { m_Controls.label_currentlyComputingStatistics->setVisible(false); mitk::StatusBar::GetInstance()->Clear(); this->UpdateIntensityProfile(); this->UpdateHistogramWidget(); } void QmitkImageStatisticsView::OnSelectedTimePointChanged(const mitk::TimePointType& /*newTimePoint*/) { this->UpdateHistogramWidget(); } void QmitkImageStatisticsView::OnJobError(QString error, const QmitkDataGenerationJobBase* /*failedJob*/) { mitk::StatusBar::GetInstance()->DisplayErrorText(error.toStdString().c_str()); MITK_WARN << "Error when calculating statistics: " << error; } void QmitkImageStatisticsView::OnRequestHistogramUpdate(unsigned int nbins) { m_Controls.widget_statistics->SetHistogramNBins(nbins); m_DataGenerator->SetHistogramNBins(nbins); this->UpdateIntensityProfile(); this->UpdateHistogramWidget(); } void QmitkImageStatisticsView::OnIgnoreZeroValuedVoxelStateChanged(int state) { auto ignoreZeroValueVoxel = (state == Qt::Unchecked) ? false : true; m_Controls.widget_statistics->SetIgnoreZeroValueVoxel(ignoreZeroValueVoxel); m_DataGenerator->SetIgnoreZeroValueVoxel(ignoreZeroValueVoxel); this->UpdateIntensityProfile(); this->UpdateHistogramWidget(); } void QmitkImageStatisticsView::OnImageSelectionChanged(QmitkAbstractNodeSelectionWidget::NodeList /*nodes*/) { auto images = m_Controls.imageNodesSelector->GetSelectedNodesStdVector(); m_Controls.widget_statistics->SetImageNodes(images); m_Controls.widget_statistics->setEnabled(!images.empty()); m_Controls.roiNodesSelector->SetNodePredicate(this->GenerateROIPredicate()); m_DataGenerator->SetAutoUpdate(false); m_DataGenerator->SetImageNodes(images); m_DataGenerator->Generate(); m_DataGenerator->SetAutoUpdate(true); this->UpdateHistogramWidget(); this->UpdateIntensityProfile(); } void QmitkImageStatisticsView::OnROISelectionChanged(QmitkAbstractNodeSelectionWidget::NodeList /*nodes*/) { auto rois = m_Controls.roiNodesSelector->GetSelectedNodesStdVector(); m_Controls.widget_statistics->SetMaskNodes(rois); m_DataGenerator->SetAutoUpdate(false); m_DataGenerator->SetROINodes(rois); m_DataGenerator->Generate(); m_DataGenerator->SetAutoUpdate(true); this->UpdateHistogramWidget(); this->UpdateIntensityProfile(); } void QmitkImageStatisticsView::OnButtonSelectionPressed() { QmitkNodeSelectionDialog* dialog = new QmitkNodeSelectionDialog(nullptr, "Select input for the statistic","You may select images and ROIs to compute their statistic. ROIs may be segmentations or planar figures."); dialog->SetDataStorage(GetDataStorage()); dialog->SetSelectionCheckFunction(CheckForSameGeometry()); // set predicates auto isPlanarFigurePredicate = mitk::GetImageStatisticsPlanarFigurePredicate(); auto isMaskPredicate = mitk::GetImageStatisticsMaskPredicate(); auto isImagePredicate = mitk::GetImageStatisticsImagePredicate(); auto isMaskOrPlanarFigurePredicate = mitk::NodePredicateOr::New(isPlanarFigurePredicate, isMaskPredicate); auto isImageOrMaskOrPlanarFigurePredicate = mitk::NodePredicateOr::New(isMaskOrPlanarFigurePredicate, isImagePredicate); dialog->SetNodePredicate(isImageOrMaskOrPlanarFigurePredicate); dialog->SetSelectionMode(QAbstractItemView::MultiSelection); dialog->SetCurrentSelection(m_Controls.imageNodesSelector->GetSelectedNodes()+m_Controls.roiNodesSelector->GetSelectedNodes()); if (dialog->exec()) { auto selectedNodeList = dialog->GetSelectedNodes(); m_Controls.imageNodesSelector->SetCurrentSelection(selectedNodeList); m_Controls.roiNodesSelector->SetCurrentSelection(selectedNodeList); } delete dialog; } QmitkNodeSelectionDialog::SelectionCheckFunctionType QmitkImageStatisticsView::CheckForSameGeometry() const { auto isMaskPredicate = mitk::GetImageStatisticsMaskPredicate(); auto lambda = [isMaskPredicate](const QmitkNodeSelectionDialog::NodeList& nodes) { if (nodes.empty()) { return std::string(); } const mitk::Image* imageNodeData = nullptr; for (auto& node : nodes) { auto castedData = dynamic_cast(node->GetData()); if (castedData != nullptr && !isMaskPredicate->CheckNode(node)) { imageNodeData = castedData; break; } } if (imageNodeData == nullptr) { std::stringstream ss; ss << "

Select at least one image.

"; return ss.str(); } auto imageGeoPredicate = mitk::NodePredicateGeometry::New(imageNodeData->GetGeometry()); auto maskGeoPredicate = mitk::NodePredicateSubGeometry::New(imageNodeData->GetGeometry()); for (auto& rightNode : nodes) { if (imageNodeData != rightNode->GetData()) { bool validGeometry = true; if (isMaskPredicate->CheckNode(rightNode)) { validGeometry = maskGeoPredicate->CheckNode(rightNode); } else if (dynamic_cast(rightNode->GetData())) { validGeometry = imageGeoPredicate->CheckNode(rightNode); } else { const auto planarFigure = dynamic_cast(rightNode->GetData()); const auto imageGeometry = imageNodeData->GetGeometry(); validGeometry = CheckPlanarFigureMatchesGeometry(planarFigure, imageGeometry); } if (!validGeometry) { std::stringstream ss; ss << "

Invalid selection: All selected nodes must have the same geometry.

Differing node i.a.: \""; ss << rightNode->GetName() <<"\"

"; return ss.str(); } } } return std::string(); }; return lambda; } mitk::NodePredicateBase::Pointer QmitkImageStatisticsView::GenerateROIPredicate() const { auto isPlanarFigurePredicate = mitk::GetImageStatisticsPlanarFigurePredicate(); auto isMaskPredicate = mitk::GetImageStatisticsMaskPredicate(); auto isMaskOrPlanarFigurePredicate = mitk::NodePredicateOr::New(isPlanarFigurePredicate, isMaskPredicate); mitk::NodePredicateBase::Pointer result = isMaskOrPlanarFigurePredicate.GetPointer(); if(!m_Controls.imageNodesSelector->GetSelectedNodes().empty()) { auto image = m_Controls.imageNodesSelector->GetSelectedNodes().front()->GetData(); auto imageGeoPredicate = mitk::NodePredicateSubGeometry::New(image->GetGeometry()); auto lambda = [image, imageGeoPredicate](const mitk::DataNode* node) { bool sameGeometry = true; if (dynamic_cast(node->GetData()) != nullptr) { sameGeometry = imageGeoPredicate->CheckNode(node); } else { const auto planarFigure = dynamic_cast(node->GetData()); sameGeometry = CheckPlanarFigureMatchesGeometry(planarFigure, image->GetGeometry()); } return sameGeometry; }; result = mitk::NodePredicateAnd::New(isMaskOrPlanarFigurePredicate, mitk::NodePredicateFunction::New(lambda)).GetPointer(); } return result; }