diff --git a/Modules/DiffusionImaging/MiniApps/ImageStatisticsMiniapp_3.cpp b/Modules/DiffusionImaging/MiniApps/ImageStatisticsMiniapp_3.cpp index a2c397c6e0..a9e73c9710 100644 --- a/Modules/DiffusionImaging/MiniApps/ImageStatisticsMiniapp_3.cpp +++ b/Modules/DiffusionImaging/MiniApps/ImageStatisticsMiniapp_3.cpp @@ -1,69 +1,72 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkCommandLineParser.h" #include "mitkImage.h" #include #include #include #include #include "mitkIOUtil.h" #include #include #include #include #include "mitkImageAccessByItk.h" #include #include #include #include #include #include #include #include #include #include +#include int main( int argc, char* argv[] ) { unsigned int timeStep = 0; std::string inputImageFile; inputImageFile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic2DplusT.nrrd"; // Load image mitk::Image::Pointer inputImage = mitk::IOUtil::LoadImage(inputImageFile); // Calculate statistics mitk::ImageStatisticsCalculator::Pointer calculator = mitk::ImageStatisticsCalculator::New(); std::cout << "calculating statistics (unmasked) itk: " << std::endl; mitk::ImageStatisticsCalculator::StatisticsContainer::Pointer result; calculator->SetInputImage(inputImage); calculator->SetNBinsForHistogramStatistics(100); for (unsigned int i=0; i < inputImage->GetTimeSteps(); i++) { std::cout << "Results for time step " << i << ":" << std::endl; result = calculator->GetStatistics(i, 1); result->Print(); std::cout << std::endl; } + mitk::Image::Pointer tmp = mitk::ImageGenerator::GenerateImageFromReference(inputImage, 0); + return EXIT_SUCCESS; } diff --git a/Modules/ImageStatistics/Testing/files.cmake b/Modules/ImageStatistics/Testing/files.cmake index b01d4dbfb1..84b7c4cf8c 100644 --- a/Modules/ImageStatistics/Testing/files.cmake +++ b/Modules/ImageStatistics/Testing/files.cmake @@ -1,11 +1,11 @@ set(MODULE_TESTS - # mitkImageStatisticsCalculatorTest.cpp - # mitkPointSetStatisticsCalculatorTest.cpp - # mitkPointSetDifferenceStatisticsCalculatorTest.cpp - # mitkImageStatisticsTextureAnalysisTest.cpp + mitkImageStatisticsCalculatorTest.cpp + mitkPointSetStatisticsCalculatorTest.cpp + mitkPointSetDifferenceStatisticsCalculatorTest.cpp + mitkImageStatisticsTextureAnalysisTest.cpp ) set(MODULE_CUSTOM_TESTS # mitkImageStatisticsHotspotTest.cpp # mitkMultiGaussianTest.cpp # TODO: activate test to generate new test cases for mitkImageStatisticsHotspotTest ) diff --git a/Modules/ImageStatistics/Testing/mitkImageStatisticsCalculatorTest.cpp b/Modules/ImageStatistics/Testing/mitkImageStatisticsCalculatorTest.cpp index 1cec357532..481c456506 100644 --- a/Modules/ImageStatistics/Testing/mitkImageStatisticsCalculatorTest.cpp +++ b/Modules/ImageStatistics/Testing/mitkImageStatisticsCalculatorTest.cpp @@ -1,507 +1,542 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkImageStatisticsCalculator.h" #include #include #include #include #include #include #include #include +#include +#include + /** * \brief Test class for mitkImageStatisticsCalculator * * This test covers: * - instantiation of an ImageStatisticsCalculator class * - correctness of statistics when using PlanarFigures for masking */ class mitkImageStatisticsCalculatorTestSuite : public mitk::TestFixture { CPPUNIT_TEST_SUITE(mitkImageStatisticsCalculatorTestSuite); MITK_TEST(TestUninitializedImage); MITK_TEST(TestCase1); MITK_TEST(TestCase2); MITK_TEST(TestCase3); MITK_TEST(TestCase4); MITK_TEST(TestCase5); MITK_TEST(TestCase6); MITK_TEST(TestCase7); MITK_TEST(TestCase8); MITK_TEST(TestCase9); MITK_TEST(TestCase10); MITK_TEST(TestCase11); MITK_TEST(TestCase12); MITK_TEST(TestImageMaskingEmpty); MITK_TEST(TestImageMaskingNonEmpty); MITK_TEST(TestRecomputeOnModifiedMask); CPPUNIT_TEST_SUITE_END(); public: void setUp() override; void TestUninitializedImage(); void TestCase1(); void TestCase2(); void TestCase3(); void TestCase4(); void TestCase5(); void TestCase6(); void TestCase7(); void TestCase8(); void TestCase9(); void TestCase10(); void TestCase11(); void TestCase12(); void TestImageMaskingEmpty(); void TestImageMaskingNonEmpty(); void TestRecomputeOnModifiedMask(); private: mitk::Image::Pointer m_Image; mitk::PlaneGeometry::Pointer m_Geometry; // calculate statistics for the given image and planarpolygon - const mitk::ImageStatisticsCalculator::Statistics ComputeStatistics( mitk::Image::Pointer image, + const mitk::ImageStatisticsCalculator::StatisticsContainer::Pointer ComputeStatistics( mitk::Image::Pointer image, mitk::PlanarFigure::Pointer polygon ); // calculate statistics for the given image and planarpolygon - const mitk::ImageStatisticsCalculator::Statistics ComputeStatistics( mitk::Image::Pointer image, + const mitk::ImageStatisticsCalculator::StatisticsContainer::Pointer ComputeStatistics( mitk::Image::Pointer image, mitk::Image::Pointer image_mask ); - void VerifyStatistics(const mitk::ImageStatisticsCalculator::Statistics& stats, + void VerifyStatistics(mitk::ImageStatisticsCalculator::StatisticsContainer::Pointer stats, double testMean, double testSD, double testMedian=0); }; void mitkImageStatisticsCalculatorTestSuite::setUp() { std::string filename = this->GetTestDataFilePath("ImageStatistics/testimage.dcm"); if (filename.empty()) { MITK_TEST_FAILED_MSG( << "Could not find test file" ) } MITK_TEST_OUTPUT(<< "Loading test image '" << filename << "'") m_Image = mitk::IOUtil::LoadImage(filename); MITK_TEST_CONDITION_REQUIRED( m_Image.IsNotNull(), "Loaded an mitk::Image" ); m_Geometry = m_Image->GetSlicedGeometry()->GetPlaneGeometry(0); MITK_TEST_CONDITION_REQUIRED( m_Geometry.IsNotNull(), "Getting image geometry" ) } void mitkImageStatisticsCalculatorTestSuite::TestCase1() { /***************************** * one whole white pixel * -> mean of 255 expected ******************************/ + MITK_INFO << std::endl << "Test case 1:"; mitk::PlanarPolygon::Pointer figure1 = mitk::PlanarPolygon::New(); figure1->SetPlaneGeometry( m_Geometry ); mitk::Point2D pnt1; pnt1[0] = 10.5 ; pnt1[1] = 3.5; figure1->PlaceFigure( pnt1 ); mitk::Point2D pnt2; pnt2[0] = 9.5; pnt2[1] = 3.5; figure1->SetControlPoint( 1, pnt2, true ); mitk::Point2D pnt3; pnt3[0] = 9.5; pnt3[1] = 4.5; figure1->SetControlPoint( 2, pnt3, true ); mitk::Point2D pnt4; pnt4[0] = 10.5; pnt4[1] = 4.5; figure1->SetControlPoint( 3, pnt4, true ); figure1->GetPolyLine(0); this->VerifyStatistics(ComputeStatistics(m_Image, figure1.GetPointer()), 255.0, 0.0, 255.0); } void mitkImageStatisticsCalculatorTestSuite::TestCase2() { /***************************** * half pixel in x-direction (white) * -> mean of 255 expected ******************************/ + MITK_INFO << std::endl << "Test case 2:"; + mitk::PlanarPolygon::Pointer figure1 = mitk::PlanarPolygon::New(); figure1->SetPlaneGeometry( m_Geometry ); mitk::Point2D pnt1; pnt1[0] = 10.0 ; pnt1[1] = 3.5; figure1->PlaceFigure( pnt1 ); mitk::Point2D pnt2; pnt2[0] = 9.5; pnt2[1] = 3.5; figure1->SetControlPoint( 1, pnt2, true ); mitk::Point2D pnt3; pnt3[0] = 9.5; pnt3[1] = 4.5; figure1->SetControlPoint( 2, pnt3, true ); mitk::Point2D pnt4; pnt4[0] = 10.0; pnt4[1] = 4.5; figure1->SetControlPoint( 3, pnt4, true ); figure1->GetPolyLine(0); this->VerifyStatistics(ComputeStatistics(m_Image, figure1.GetPointer()), 255.0, 0.0, 255.0); } void mitkImageStatisticsCalculatorTestSuite::TestCase3() { /***************************** * half pixel in diagonal-direction (white) * -> mean of 255 expected ******************************/ + MITK_INFO << std::endl << "Test case 3:"; mitk::PlanarPolygon::Pointer figure1 = mitk::PlanarPolygon::New(); figure1->SetPlaneGeometry( m_Geometry ); mitk::Point2D pnt1; pnt1[0] = 10.5 ; pnt1[1] = 3.5; figure1->PlaceFigure( pnt1 ); mitk::Point2D pnt2; pnt2[0] = 9.5; pnt2[1] = 3.5; figure1->SetControlPoint( 1, pnt2, true ); mitk::Point2D pnt3; pnt3[0] = 9.5; pnt3[1] = 4.5; figure1->SetControlPoint( 2, pnt3, true ); figure1->GetPolyLine(0); this->VerifyStatistics(ComputeStatistics(m_Image, figure1.GetPointer()), 255.0, 0.0, 255.0); } void mitkImageStatisticsCalculatorTestSuite::TestCase4() { /***************************** * one pixel (white) + 2 half pixels (white) + 1 half pixel (black) * -> mean of 191.25 expected ******************************/ + MITK_INFO << std::endl << "Test case 4:"; mitk::PlanarPolygon::Pointer figure1 = mitk::PlanarPolygon::New(); figure1->SetPlaneGeometry( m_Geometry ); mitk::Point2D pnt1; pnt1[0] = 1.1; pnt1[1] = 1.1; figure1->PlaceFigure( pnt1 ); mitk::Point2D pnt2; pnt2[0] = 2.0; pnt2[1] = 2.0; figure1->SetControlPoint( 1, pnt2, true ); mitk::Point2D pnt3; pnt3[0] = 3.0; pnt3[1] = 1.0; figure1->SetControlPoint( 2, pnt3, true ); mitk::Point2D pnt4; pnt4[0] = 2.0; pnt4[1] = 0.0; figure1->SetControlPoint( 3, pnt4, true ); figure1->GetPolyLine(0); - this->VerifyStatistics(ComputeStatistics(m_Image, figure1.GetPointer()), 191.25, 127.5, 254.50); + this->VerifyStatistics(ComputeStatistics(m_Image, figure1.GetPointer()), 191.25, 110.41, 242.250); } void mitkImageStatisticsCalculatorTestSuite::TestCase5() { /***************************** * whole pixel (white) + half pixel (gray) in x-direction * -> mean of 191.5 expected ******************************/ + MITK_INFO << std::endl << "Test case 5:"; mitk::PlanarPolygon::Pointer figure1 = mitk::PlanarPolygon::New(); figure1->SetPlaneGeometry( m_Geometry ); mitk::Point2D pnt1; pnt1[0] = 11.0; pnt1[1] = 3.5; figure1->PlaceFigure( pnt1 ); mitk::Point2D pnt2; pnt2[0] = 9.5; pnt2[1] = 3.5; figure1->SetControlPoint( 1, pnt2, true ); mitk::Point2D pnt3; pnt3[0] = 9.5; pnt3[1] = 4.5; figure1->SetControlPoint( 2, pnt3, true ); mitk::Point2D pnt4; pnt4[0] = 11.0; pnt4[1] = 4.5; figure1->SetControlPoint( 3, pnt4, true ); figure1->GetPolyLine(0); - this->VerifyStatistics(ComputeStatistics(m_Image, figure1.GetPointer()), 191.50, 89.80, 254.5); + this->VerifyStatistics(ComputeStatistics(m_Image, figure1.GetPointer()), 191.50, 63.50, 134.340); } void mitkImageStatisticsCalculatorTestSuite::TestCase6() { /***************************** * quarter pixel (black) + whole pixel (white) + half pixel (gray) in x-direction * -> mean of 191.5 expected ******************************/ + MITK_INFO << std::endl << "Test case 6:"; mitk::PlanarPolygon::Pointer figure1 = mitk::PlanarPolygon::New(); figure1->SetPlaneGeometry( m_Geometry ); mitk::Point2D pnt1; pnt1[0] = 11.0; pnt1[1] = 3.5; figure1->PlaceFigure( pnt1 ); mitk::Point2D pnt2; pnt2[0] = 9.25; pnt2[1] = 3.5; figure1->SetControlPoint( 1, pnt2, true ); mitk::Point2D pnt3; pnt3[0] = 9.25; pnt3[1] = 4.5; figure1->SetControlPoint( 2, pnt3, true ); mitk::Point2D pnt4; pnt4[0] = 11.0; pnt4[1] = 4.5; figure1->SetControlPoint( 3, pnt4, true ); figure1->GetPolyLine(0); - this->VerifyStatistics(ComputeStatistics(m_Image, figure1.GetPointer()), 191.5, 89.80, 254.5); + this->VerifyStatistics(ComputeStatistics(m_Image, figure1.GetPointer()), 191.5, 63.50, 134.340); } void mitkImageStatisticsCalculatorTestSuite::TestCase7() { /***************************** * half pixel (black) + whole pixel (white) + half pixel (gray) in x-direction * -> mean of 127.66 expected ******************************/ + MITK_INFO << std::endl << "Test case 7:"; + mitk::PlanarPolygon::Pointer figure1 = mitk::PlanarPolygon::New(); figure1->SetPlaneGeometry( m_Geometry ); mitk::Point2D pnt1; pnt1[0] = 11.0; pnt1[1] = 3.5; figure1->PlaceFigure( pnt1 ); mitk::Point2D pnt2; pnt2[0] = 9.0; pnt2[1] = 3.5; figure1->SetControlPoint( 1, pnt2, true ); mitk::Point2D pnt3; pnt3[0] = 9.0; pnt3[1] = 4.0; figure1->SetControlPoint( 2, pnt3, true ); mitk::Point2D pnt4; pnt4[0] = 11.0; pnt4[1] = 4.0; figure1->SetControlPoint( 3, pnt4, true ); figure1->GetPolyLine(0); - this->VerifyStatistics(ComputeStatistics(m_Image, figure1.GetPointer()), 127.66, 127.5, 128.5); + this->VerifyStatistics(ComputeStatistics(m_Image, figure1.GetPointer()), 127.66, 104.1, 140.250); } void mitkImageStatisticsCalculatorTestSuite::TestCase8() { /***************************** * whole pixel (gray) * -> mean of 128 expected ******************************/ + MITK_INFO << std::endl << "Test case 8:"; + mitk::PlanarPolygon::Pointer figure2 = mitk::PlanarPolygon::New(); figure2->SetPlaneGeometry( m_Geometry ); mitk::Point2D pnt1; pnt1[0] = 11.5; pnt1[1] = 10.5; figure2->PlaceFigure( pnt1 ); mitk::Point2D pnt2; pnt2[0] = 11.5; pnt2[1] = 11.5; figure2->SetControlPoint( 1, pnt2, true ); mitk::Point2D pnt3; pnt3[0] = 12.5; pnt3[1] = 11.5; figure2->SetControlPoint( 2, pnt3, true ); mitk::Point2D pnt4; pnt4[0] = 12.5; pnt4[1] = 10.5; figure2->SetControlPoint( 3, pnt4, true ); figure2->GetPolyLine(0); this->VerifyStatistics(ComputeStatistics(m_Image, figure2.GetPointer()), 128.0, 0.0, 128.0); } void mitkImageStatisticsCalculatorTestSuite::TestCase9() { /***************************** * whole pixel (gray) + half pixel (white) in y-direction * -> mean of 191.5 expected ******************************/ + MITK_INFO << std::endl << "Test case 9:"; + mitk::PlanarPolygon::Pointer figure2 = mitk::PlanarPolygon::New(); figure2->SetPlaneGeometry( m_Geometry ); mitk::Point2D pnt1; pnt1[0] = 11.5; pnt1[1] = 10.5; figure2->PlaceFigure( pnt1 ); mitk::Point2D pnt2; pnt2[0] = 11.5; pnt2[1] = 12.0; figure2->SetControlPoint( 1, pnt2, true ); mitk::Point2D pnt3; pnt3[0] = 12.5; pnt3[1] = 12.0; figure2->SetControlPoint( 2, pnt3, true ); mitk::Point2D pnt4; pnt4[0] = 12.5; pnt4[1] = 10.5; figure2->SetControlPoint( 3, pnt4, true ); figure2->GetPolyLine(0); - this->VerifyStatistics(ComputeStatistics(m_Image, figure2.GetPointer()), 191.5, 89.80, 254.5); + this->VerifyStatistics(ComputeStatistics(m_Image, figure2.GetPointer()), 191.5, 63.50, 134.340); } void mitkImageStatisticsCalculatorTestSuite::TestCase10() { /***************************** * 2 whole pixel (white) + 2 whole pixel (black) in y-direction * -> mean of 127.66 expected ******************************/ + MITK_INFO << std::endl << "Test case 10:"; + mitk::PlanarPolygon::Pointer figure2 = mitk::PlanarPolygon::New(); figure2->SetPlaneGeometry( m_Geometry ); mitk::Point2D pnt1; pnt1[0] = 11.5; pnt1[1] = 10.5; figure2->PlaceFigure( pnt1 ); mitk::Point2D pnt2; pnt2[0] = 11.5; pnt2[1] = 13.5; figure2->SetControlPoint( 1, pnt2, true ); mitk::Point2D pnt3; pnt3[0] = 12.5; pnt3[1] = 13.5; figure2->SetControlPoint( 2, pnt3, true ); mitk::Point2D pnt4; pnt4[0] = 12.5; pnt4[1] = 10.5; figure2->SetControlPoint( 3, pnt4, true ); figure2->GetPolyLine(0); - this->VerifyStatistics(ComputeStatistics(m_Image, figure2.GetPointer()), 127.66, 127.5, 128.5); + this->VerifyStatistics(ComputeStatistics(m_Image, figure2.GetPointer()), 127.66, 104.1, 140.250); } void mitkImageStatisticsCalculatorTestSuite::TestCase11() { /***************************** * 9 whole pixels (white) + 3 half pixels (white) * + 3 whole pixel (black) [ + 3 slightly less than half pixels (black)] * -> mean of 204.0 expected ******************************/ + MITK_INFO << std::endl << "Test case 11:"; + mitk::PlanarPolygon::Pointer figure2 = mitk::PlanarPolygon::New(); figure2->SetPlaneGeometry( m_Geometry ); mitk::Point2D pnt1; pnt1[0] = 0.5; pnt1[1] = 0.5; figure2->PlaceFigure( pnt1 ); mitk::Point2D pnt2; pnt2[0] = 3.5; pnt2[1] = 3.5; figure2->SetControlPoint( 1, pnt2, true ); mitk::Point2D pnt3; pnt3[0] = 8.4999; pnt3[1] = 3.5; figure2->SetControlPoint( 2, pnt3, true ); mitk::Point2D pnt4; pnt4[0] = 5.4999; pnt4[1] = 0.5; figure2->SetControlPoint( 3, pnt4, true ); figure2->GetPolyLine(0); - this->VerifyStatistics(ComputeStatistics(m_Image, figure2.GetPointer()), 204.0, 105.58, 254.5); + this->VerifyStatistics(ComputeStatistics(m_Image, figure2.GetPointer()), 204.0, 102.00, 242.250); } void mitkImageStatisticsCalculatorTestSuite::TestCase12() { /***************************** * half pixel (white) + whole pixel (white) + half pixel (black) * -> mean of 212.66 expected ******************************/ + MITK_INFO << std::endl << "Test case 12:"; + mitk::PlanarPolygon::Pointer figure2 = mitk::PlanarPolygon::New(); figure2->SetPlaneGeometry( m_Geometry ); mitk::Point2D pnt1; pnt1[0] = 9.5; pnt1[1] = 0.5; figure2->PlaceFigure( pnt1 ); mitk::Point2D pnt2; pnt2[0] = 9.5; pnt2[1] = 2.5; figure2->SetControlPoint( 1, pnt2, true ); mitk::Point2D pnt3; pnt3[0] = 11.5; pnt3[1] = 2.5; figure2->SetControlPoint( 2, pnt3, true ); figure2->GetPolyLine(0); - this->VerifyStatistics(ComputeStatistics(m_Image, figure2.GetPointer()), 212.66, 73.32, 254.5); + this->VerifyStatistics(ComputeStatistics(m_Image, figure2.GetPointer()), 212.66, 59.860, 248.640); } void mitkImageStatisticsCalculatorTestSuite::TestImageMaskingEmpty() { + MITK_INFO << std::endl << "TestImageMaskingEmpty:"; + mitk::Image::Pointer mask_image = mitk::ImageGenerator::GenerateImageFromReference( m_Image, 0 ); - this->VerifyStatistics( ComputeStatistics( m_Image, mask_image ), 0.0, 0.0, 0.0); + this->VerifyStatistics( ComputeStatistics( m_Image, mask_image ), -21474836.480, -21474836.480, -21474836.480); // empty statisticsContainer (default values) } void mitkImageStatisticsCalculatorTestSuite::TestImageMaskingNonEmpty() { + MITK_INFO << std::endl << "TestImageMaskingNonEmpty:"; + mitk::Image::Pointer mask_image = mitk::ImageGenerator::GenerateImageFromReference( m_Image, 0 ); std::vector< itk::Index<3U> > activated_indices; itk::Index<3U> index = {{10, 8, 0}}; activated_indices.push_back( index ); index[0] = 9; index[1] = 8; index[2] = 0; activated_indices.push_back( index ); index[0] = 9; index[1] = 7; index[2] = 0; activated_indices.push_back( index ); index[0] = 10; index[1] = 7; index[2] = 0; activated_indices.push_back( index ); std::vector< itk::Index<3U> >::const_iterator indexIter = activated_indices.begin(); // activate voxel in the mask image mitk::ImagePixelWriteAccessor< unsigned char, 3> writeAccess( mask_image ); while( indexIter != activated_indices.end() ) { writeAccess.SetPixelByIndex( (*indexIter++), 1); } - this->VerifyStatistics( ComputeStatistics( m_Image, mask_image ), 127.5, 147.22, 254.5); + this->VerifyStatistics( ComputeStatistics( m_Image, mask_image ), 127.5, 127.5, 12.750); } void mitkImageStatisticsCalculatorTestSuite::TestRecomputeOnModifiedMask() { + MITK_INFO << std::endl << "TestRecomputeOnModifiedMask:"; mitk::Image::Pointer mask_image = mitk::ImageGenerator::GenerateImageFromReference( m_Image, 0 ); mitk::ImageStatisticsCalculator::Pointer statisticsCalculator = mitk::ImageStatisticsCalculator::New(); - statisticsCalculator->SetImage( m_Image ); - statisticsCalculator->SetImageMask( mask_image ); - statisticsCalculator->SetMaskingModeToImage(); + statisticsCalculator->SetInputImage( m_Image ); + + mitk::ImageMaskGenerator::Pointer imgMaskGen = mitk::ImageMaskGenerator::New(); + imgMaskGen->SetImageMask(mask_image); + + statisticsCalculator->SetMask(imgMaskGen.GetPointer()); - statisticsCalculator->ComputeStatistics(); - this->VerifyStatistics( statisticsCalculator->GetStatistics(), 0.0, 0.0, 0.0); + this->VerifyStatistics( statisticsCalculator->GetStatistics(), -21474836.480, -21474836.480, -21474836.480); // activate voxel in the mask image itk::Index<3U> test_index = {11, 8, 0}; mitk::ImagePixelWriteAccessor< unsigned char, 3> writeAccess( mask_image ); writeAccess.SetPixelByIndex( test_index, 1); mask_image->Modified(); - statisticsCalculator->ComputeStatistics(); - const mitk::ImageStatisticsCalculator::Statistics stat = statisticsCalculator->GetStatistics(); + mitk::ImageStatisticsCalculator::StatisticsContainer::Pointer stat = statisticsCalculator->GetStatistics(); this->VerifyStatistics( stat, 128.0, 0.0, 128.0); - MITK_TEST_CONDITION( stat.GetN() == 1, "Calculated mask voxel count '" << stat.GetN() << "' is equal to the desired value '" << 1 << "'" ); + MITK_TEST_CONDITION( stat->GetN() == 1, "Calculated mask voxel count '" << stat->GetN() << "' is equal to the desired value '" << 1 << "'" ); } -const mitk::ImageStatisticsCalculator::Statistics +const mitk::ImageStatisticsCalculator::StatisticsContainer::Pointer mitkImageStatisticsCalculatorTestSuite::ComputeStatistics( mitk::Image::Pointer image, mitk::PlanarFigure::Pointer polygon ) { mitk::ImageStatisticsCalculator::Pointer statisticsCalculator = mitk::ImageStatisticsCalculator::New(); - statisticsCalculator->SetImage( image ); - statisticsCalculator->SetMaskingModeToPlanarFigure(); - statisticsCalculator->SetPlanarFigure( polygon ); + statisticsCalculator->SetInputImage( image ); + + statisticsCalculator->SetNBinsForHistogramStatistics(10); + + mitk::PlanarFigureMaskGenerator::Pointer planFigMaskGen = mitk::PlanarFigureMaskGenerator::New(); + planFigMaskGen->SetImage(image); + planFigMaskGen->SetPlanarFigure(polygon); + + statisticsCalculator->SetMask(planFigMaskGen.GetPointer()); try { - statisticsCalculator->ComputeStatistics(); return statisticsCalculator->GetStatistics(); } catch( ... ) { } - return mitk::ImageStatisticsCalculator::Statistics(); + return mitk::ImageStatisticsCalculator::StatisticsContainer::New(); } -const mitk::ImageStatisticsCalculator::Statistics +const mitk::ImageStatisticsCalculator::StatisticsContainer::Pointer mitkImageStatisticsCalculatorTestSuite::ComputeStatistics(mitk::Image::Pointer image, mitk::Image::Pointer image_mask ) { mitk::ImageStatisticsCalculator::Pointer statisticsCalculator = mitk::ImageStatisticsCalculator::New(); - statisticsCalculator->SetImage( image ); - statisticsCalculator->SetImageMask( image_mask ); - statisticsCalculator->SetMaskingModeToImage(); + statisticsCalculator->SetInputImage(image); + + statisticsCalculator->SetNBinsForHistogramStatistics(10); - statisticsCalculator->ComputeStatistics(); + mitk::ImageMaskGenerator::Pointer imgMaskGen = mitk::ImageMaskGenerator::New(); + imgMaskGen->SetImageMask(image_mask); + statisticsCalculator->SetMask(imgMaskGen.GetPointer()); return statisticsCalculator->GetStatistics(); } -void mitkImageStatisticsCalculatorTestSuite::VerifyStatistics(const mitk::ImageStatisticsCalculator::Statistics& stats, +void mitkImageStatisticsCalculatorTestSuite::VerifyStatistics(mitk::ImageStatisticsCalculator::StatisticsContainer::Pointer stats, double testMean, double testSD, double testMedian) { - int tmpMean = stats.GetMean() * 100; + int tmpMean = stats->GetMean() * 100; double calculatedMean = tmpMean / 100.0; MITK_TEST_CONDITION( calculatedMean == testMean, "Calculated mean grayvalue '" << calculatedMean << "' is equal to the desired value '" << testMean << "'" ); - int tmpSD = stats.GetSigma() * 100; + int tmpSD = stats->GetStd() * 100; double calculatedSD = tmpSD / 100.0; MITK_TEST_CONDITION( calculatedSD == testSD, "Calculated grayvalue sd '" << calculatedSD << "' is equal to the desired value '" << testSD <<"'" ); - int tmpMedian = stats.GetMedian() * 100; + int tmpMedian = stats->GetMedian() * 100; double calculatedMedian = tmpMedian / 100.0; MITK_TEST_CONDITION( testMedian == calculatedMedian, "Calculated median grayvalue '" << calculatedMedian << "' is equal to the desired value '" << testMedian << "'"); } void mitkImageStatisticsCalculatorTestSuite::TestUninitializedImage() { /***************************** * loading uninitialized image to datastorage ******************************/ MITK_TEST_FOR_EXCEPTION_BEGIN(mitk::Exception) mitk::Image::Pointer image = mitk::Image::New(); mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData(image); mitk::ImageStatisticsCalculator::Pointer is = mitk::ImageStatisticsCalculator::New(); - is->ComputeStatistics(); + is->GetStatistics(); MITK_TEST_FOR_EXCEPTION_END(mitk::Exception) } MITK_TEST_SUITE_REGISTRATION(mitkImageStatisticsCalculator) diff --git a/Modules/ImageStatistics/Testing/mitkImageStatisticsHotspotTest.cpp b/Modules/ImageStatistics/Testing/mitkImageStatisticsHotspotTest.cpp index baa2df1120..c1fae82f5c 100644 --- a/Modules/ImageStatistics/Testing/mitkImageStatisticsHotspotTest.cpp +++ b/Modules/ImageStatistics/Testing/mitkImageStatisticsHotspotTest.cpp @@ -1,632 +1,632 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkImageStatisticsCalculator.h" #include "itkMultiGaussianImageSource.h" #include "mitkTestingMacros.h" #include "mitkImageCast.h" #include #include #include #include /** \section hotspotCalculationTestCases Testcases To see the different Hotspot-Testcases have a look at the \ref hotspottestdoc. Note from an intensive session of checking the test results: - itk::MultiGaussianImageSource needs a review - the test idea is ok, but the combination of XML files for parameters and MultiGaussianImageSource has serious flaws - the XML file should contain exactly the parameters that MultiGaussianImageSource requires - in contrast, now the XML file mentions index coordinates for gaussian centers while the MultiGaussianImageSource expects world coordinates - this requires a transformation (index * spacing assuming no rotation) that was actually broken until recently */ struct mitkImageStatisticsHotspotTestClass { /** \brief Test parameters for one test case. Describes all aspects of a single test case: - parameters to generate a test image - parameters of a ROI that describes where to calculate statistics - expected statistics results */ struct Parameters { public: // XML-Tag /** \brief XML-Tag "image-rows": size of x-dimension */ int m_ImageRows; /** \brief XML-Tag "image-columns": size of y-dimension */ int m_ImageColumns; /** \brief XML-Tag "image-slices": size of z-dimension */ int m_ImageSlices; /** \brief XML-Tag "numberOfGaussians": number of used gauss-functions */ int m_NumberOfGaussian; /** \brief XML-Tags "spacingX", "spacingY", "spacingZ": spacing of image in every direction */ double m_Spacing[3]; /** \brief XML-Tag "entireHotSpotInImage" */ unsigned int m_EntireHotspotInImage; // XML-Tag /** \brief XML-Tag "centerIndexX: gaussian parameter \warning This parameter READS the centerIndexX parameter from file and is THEN MISUSED to calculate some position in world coordinates, so we require double. */ std::vector m_CenterX; /** \brief XML-Tag "centerIndexY: gaussian parameter \warning This parameter READS the centerIndexX parameter from file and is THEN MISUSED to calculate some position in world coordinates, so we require double. */ std::vector m_CenterY; /** \brief XML-Tag "centerIndexZ: gaussian parameter \warning This parameter READS the centerIndexX parameter from file and is THEN MISUSED to calculate some position in world coordinates, so we require double. */ std::vector m_CenterZ; /** \brief XML-Tag "deviationX: gaussian parameter */ std::vector m_SigmaX; /** \brief XML-Tag "deviationY: gaussian parameter */ std::vector m_SigmaY; /** \brief XML-Tag "deviationZ: gaussian parameter */ std::vector m_SigmaZ; /** \brief XML-Tag "altitude: gaussian parameter */ std::vector m_Altitude; // XML-Tag /** \brief XML-Tag "numberOfLabels": number of different labels which appear in the mask */ unsigned int m_NumberOfLabels; /** \brief XML-Tag "hotspotRadiusInMM": radius of hotspot */ double m_HotspotRadiusInMM; // XML-Tag /** \brief XML-Tag "maximumSizeX": maximum position of ROI in x-dimension */ vnl_vector m_MaxIndexX; /** \brief XML-Tag "minimumSizeX": minimum position of ROI in x-dimension */ vnl_vector m_MinIndexX; /** \brief XML-Tag "maximumSizeX": maximum position of ROI in y-dimension */ vnl_vector m_MaxIndexY; /** \brief XML-Tag "minimumSizeX": minimum position of ROI in y-dimension */ vnl_vector m_MinIndexY; /** \brief XML-Tag "maximumSizeX": maximum position of ROI in z-dimension */ vnl_vector m_MaxIndexZ; /** \brief XML-Tag "minimumSizeX": minimum position of ROI in z-dimension */ vnl_vector m_MinIndexZ; /** \brief XML-Tag "label": value of label */ vnl_vector m_Label; //XML-Tag /** \brief XML-Tag "minimum": minimum inside hotspot */ vnl_vector m_HotspotMin; /** \brief XML-Tag "maximum": maximum inside hotspot */ vnl_vector m_HotspotMax; /** \brief XML-Tag "mean": mean value of hotspot */ vnl_vector m_HotspotMean; /** \brief XML-Tag "maximumIndexX": x-coordinate of maximum-location inside hotspot */ vnl_vector m_HotspotMaxIndexX; /** \brief XML-Tag "maximumIndexX": y-coordinate of maximum-location inside hotspot */ vnl_vector m_HotspotMaxIndexY; /** \brief XML-Tag "maximumIndexX": z-coordinate of maximum-location inside hotspot */ vnl_vector m_HotspotMaxIndexZ; /** \brief XML-Tag "maximumIndexX": x-coordinate of maximum-location inside hotspot */ vnl_vector m_HotspotMinIndexX; /** \brief XML-Tag "maximumIndexX": y-coordinate of maximum-location inside hotspot */ vnl_vector m_HotspotMinIndexY; /** \brief XML-Tag "maximumIndexX": z-coordinate of maximum-location inside hotspot */ vnl_vector m_HotspotMinIndexZ; /** \brief XML-Tag "maximumIndexX": x-coordinate of hotspot-location */ vnl_vector m_HotspotIndexX; /** \brief XML-Tag "maximumIndexX": y-coordinate of hotspot-location */ vnl_vector m_HotspotIndexY; /** \brief XML-Tag "maximumIndexX": z-coordinate of hotspot-location */ vnl_vector m_HotspotIndexZ; }; /** \brief Find/Convert integer attribute in itk::DOMNode. */ static int GetIntegerAttribute(itk::DOMNode* domNode, const std::string& tag) { assert(domNode); MITK_TEST_CONDITION_REQUIRED( domNode->HasAttribute(tag), "Tag '" << tag << "' is defined in test parameters" ); std::string attributeValue = domNode->GetAttribute(tag); int resultValue; try { //MITK_TEST_OUTPUT( << "Converting tag value '" << attributeValue << "' for tag '" << tag << "' to integer"); std::stringstream(attributeValue) >> resultValue; return resultValue; } catch(std::exception& /*e*/) { MITK_TEST_CONDITION_REQUIRED(false, "Convert tag value '" << attributeValue << "' for tag '" << tag << "' to integer"); return 0; // just to satisfy compiler } } /** \brief Find/Convert double attribute in itk::DOMNode. */ static double GetDoubleAttribute(itk::DOMNode* domNode, const std::string& tag) { assert(domNode); MITK_TEST_CONDITION_REQUIRED( domNode->HasAttribute(tag), "Tag '" << tag << "' is defined in test parameters" ); std::string attributeValue = domNode->GetAttribute(tag); double resultValue; try { //MITK_TEST_OUTPUT( << "Converting tag value '" << attributeValue << "' for tag '" << tag << "' to double"); std::stringstream(attributeValue) >> resultValue; return resultValue; } catch(std::exception& /*e*/) { MITK_TEST_CONDITION_REQUIRED(false, "Convert tag value '" << attributeValue << "' for tag '" << tag << "' to double"); return 0.0; // just to satisfy compiler } } /** \brief Read XML file describing the test parameters. Reads XML file given in first commandline parameter in order to construct a Parameters structure. The XML file should be structurs as the following example, i.e. we describe the three test aspects of Parameters in four different tags, with all the details described as tag attributes. */ /** \verbatim \endverbatim */ static Parameters ParseParameters(int argc, char* argv[]) { MITK_TEST_CONDITION_REQUIRED(argc == 2, "Test is invoked with exactly 1 parameter (XML parameters file)"); MITK_INFO << "Reading parameters from file '" << argv[1] << "'"; std::string filename = argv[1]; Parameters result; itk::DOMNodeXMLReader::Pointer xmlReader = itk::DOMNodeXMLReader::New(); xmlReader->SetFileName( filename ); try { xmlReader->Update(); itk::DOMNode::Pointer domRoot = xmlReader->GetOutput(); typedef std::vector NodeList; NodeList testimages; domRoot->GetChildren("testimage", testimages); MITK_TEST_CONDITION_REQUIRED( testimages.size() == 1, "One test image defined" ) itk::DOMNode* testimage = testimages[0]; result.m_ImageRows = GetIntegerAttribute( testimage, "image-rows" ); result.m_ImageColumns = GetIntegerAttribute( testimage, "image-columns" ); result.m_ImageSlices = GetIntegerAttribute( testimage, "image-slices" ); result.m_NumberOfGaussian = GetIntegerAttribute( testimage, "numberOfGaussians" ); result.m_Spacing[0] = GetDoubleAttribute(testimage, "spacingX"); result.m_Spacing[1] = GetDoubleAttribute(testimage, "spacingY"); result.m_Spacing[2] = GetDoubleAttribute(testimage, "spacingZ"); result.m_EntireHotspotInImage = GetIntegerAttribute( testimage, "entireHotSpotInImage" ); MITK_TEST_OUTPUT( << "Read size parameters (x,y,z): " << result.m_ImageRows << "," << result.m_ImageColumns << "," << result.m_ImageSlices); MITK_TEST_OUTPUT( << "Read spacing parameters (x,y,z): " << result.m_Spacing[0] << "," << result.m_Spacing[1] << "," << result.m_Spacing[2]); NodeList gaussians; testimage->GetChildren("gaussian", gaussians); MITK_TEST_CONDITION_REQUIRED( gaussians.size() >= 1, "At least one gaussian is defined" ) result.m_CenterX.resize(result.m_NumberOfGaussian); result.m_CenterY.resize(result.m_NumberOfGaussian); result.m_CenterZ.resize(result.m_NumberOfGaussian); result.m_SigmaX.resize(result.m_NumberOfGaussian); result.m_SigmaY.resize(result.m_NumberOfGaussian); result.m_SigmaZ.resize(result.m_NumberOfGaussian); result.m_Altitude.resize(result.m_NumberOfGaussian); for(int i = 0; i < result.m_NumberOfGaussian ; ++i) { itk::DOMNode* gaussian = gaussians[i]; result.m_CenterX[i] = GetIntegerAttribute(gaussian, "centerIndexX"); result.m_CenterY[i] = GetIntegerAttribute(gaussian, "centerIndexY"); result.m_CenterZ[i] = GetIntegerAttribute(gaussian, "centerIndexZ"); result.m_SigmaX[i] = GetDoubleAttribute(gaussian, "deviationX"); result.m_SigmaY[i] = GetDoubleAttribute(gaussian, "deviationY"); result.m_SigmaZ[i] = GetDoubleAttribute(gaussian, "deviationZ"); result.m_Altitude[i] = GetDoubleAttribute(gaussian, "altitude"); result.m_CenterX[i] = result.m_CenterX[i] * result.m_Spacing[0]; result.m_CenterY[i] = result.m_CenterY[i] * result.m_Spacing[1]; result.m_CenterZ[i] = result.m_CenterZ[i] * result.m_Spacing[2]; result.m_SigmaX[i] = result.m_SigmaX[i] * result.m_Spacing[0]; result.m_SigmaY[i] = result.m_SigmaY[i] * result.m_Spacing[1]; result.m_SigmaZ[i] = result.m_SigmaZ[i] * result.m_Spacing[2]; } NodeList segmentations; domRoot->GetChildren("segmentation", segmentations); MITK_TEST_CONDITION_REQUIRED( segmentations.size() == 1, "One segmentation defined"); itk::DOMNode* segmentation = segmentations[0]; result.m_NumberOfLabels = GetIntegerAttribute(segmentation, "numberOfLabels"); result.m_HotspotRadiusInMM = GetDoubleAttribute(segmentation, "hotspotRadiusInMM"); // read ROI parameters, fill result structure NodeList rois; segmentation->GetChildren("roi", rois); MITK_TEST_CONDITION_REQUIRED( rois.size() >= 1, "At least one ROI defined" ) result.m_MaxIndexX.set_size(result.m_NumberOfLabels); result.m_MinIndexX.set_size(result.m_NumberOfLabels); result.m_MaxIndexY.set_size(result.m_NumberOfLabels); result.m_MinIndexY.set_size(result.m_NumberOfLabels); result.m_MaxIndexZ.set_size(result.m_NumberOfLabels); result.m_MinIndexZ.set_size(result.m_NumberOfLabels); result.m_Label.set_size(result.m_NumberOfLabels); for(unsigned int i = 0; i < rois.size(); ++i) { result.m_MaxIndexX[i] = GetIntegerAttribute(rois[i], "maximumIndexX"); result.m_MinIndexX[i] = GetIntegerAttribute(rois[i], "minimumIndexX"); result.m_MaxIndexY[i] = GetIntegerAttribute(rois[i], "maximumIndexY"); result.m_MinIndexY[i] = GetIntegerAttribute(rois[i], "minimumIndexY"); result.m_MaxIndexZ[i] = GetIntegerAttribute(rois[i], "maximumIndexZ"); result.m_MinIndexZ[i] = GetIntegerAttribute(rois[i], "minimumIndexZ"); result.m_Label[i] = GetIntegerAttribute(rois[i], "label"); } // read statistic parameters, fill result structure NodeList statistics; domRoot->GetChildren("statistic", statistics); MITK_TEST_CONDITION_REQUIRED( statistics.size() >= 1 , "At least one statistic defined" ) MITK_TEST_CONDITION_REQUIRED( statistics.size() == rois.size(), "Same number of rois and corresponding statistics defined"); result.m_HotspotMin.set_size(statistics.size()); result.m_HotspotMax.set_size(statistics.size()); result.m_HotspotMean.set_size(statistics.size()); result.m_HotspotMinIndexX.set_size(statistics.size()); result.m_HotspotMinIndexY.set_size(statistics.size()); result.m_HotspotMinIndexZ.set_size(statistics.size()); result.m_HotspotMaxIndexX.set_size(statistics.size()); result.m_HotspotMaxIndexY.set_size(statistics.size()); result.m_HotspotMaxIndexZ.set_size(statistics.size()); result.m_HotspotIndexX.set_size(statistics.size()); result.m_HotspotIndexY.set_size(statistics.size()); result.m_HotspotIndexZ.set_size(statistics.size()); for(unsigned int i = 0; i < statistics.size(); ++i) { result.m_HotspotMin[i] = GetDoubleAttribute(statistics[i], "minimum"); result.m_HotspotMax[i] = GetDoubleAttribute(statistics[i], "maximum"); result.m_HotspotMean[i] = GetDoubleAttribute(statistics[i], "mean"); result.m_HotspotMinIndexX[i] = GetIntegerAttribute(statistics[i], "minimumIndexX"); result.m_HotspotMinIndexY[i] = GetIntegerAttribute(statistics[i], "minimumIndexY"); result.m_HotspotMinIndexZ[i] = GetIntegerAttribute(statistics[i], "minimumIndexZ"); result.m_HotspotMaxIndexX[i] = GetIntegerAttribute(statistics[i], "maximumIndexX"); result.m_HotspotMaxIndexY[i] = GetIntegerAttribute(statistics[i], "maximumIndexY"); result.m_HotspotMaxIndexZ[i] = GetIntegerAttribute(statistics[i], "maximumIndexZ"); result.m_HotspotIndexX[i] = GetIntegerAttribute(statistics[i], "hotspotIndexX"); result.m_HotspotIndexY[i] = GetIntegerAttribute(statistics[i], "hotspotIndexY"); result.m_HotspotIndexZ[i] = GetIntegerAttribute(statistics[i], "hotspotIndexZ"); } } catch (std::exception& e) { MITK_TEST_CONDITION_REQUIRED(false, "Reading test parameters from XML file. Error message: " << e.what()); } return result; } /** \brief Generate an image that contains a couple of 3D gaussian distributions. Uses the given parameters to produce a test image using class MultiGaussianImageSource. */ static mitk::Image::Pointer BuildTestImage(const Parameters& testParameters) { mitk::Image::Pointer result; typedef double PixelType; const int Dimension = 3; typedef itk::Image ImageType; ImageType::Pointer image = ImageType::New(); typedef itk::MultiGaussianImageSource< ImageType > MultiGaussianImageSource; MultiGaussianImageSource::Pointer gaussianGenerator = MultiGaussianImageSource::New(); ImageType::SizeValueType size[3]; size[0] = testParameters.m_ImageColumns; size[1] = testParameters.m_ImageRows; size[2] = testParameters.m_ImageSlices; itk::MultiGaussianImageSource::VectorType centerXVec, centerYVec, centerZVec, sigmaXVec, sigmaYVec, sigmaZVec, altitudeVec; for(int i = 0; i < testParameters.m_NumberOfGaussian; ++i) { centerXVec.push_back(testParameters.m_CenterX[i]); centerYVec.push_back(testParameters.m_CenterY[i]); centerZVec.push_back(testParameters.m_CenterZ[i]); sigmaXVec.push_back(testParameters.m_SigmaX[i]); sigmaYVec.push_back(testParameters.m_SigmaY[i]); sigmaZVec.push_back(testParameters.m_SigmaZ[i]); altitudeVec.push_back(testParameters.m_Altitude[i]); } ImageType::SpacingType spacing; for( int i = 0; i < Dimension; ++i ) spacing[i] = testParameters.m_Spacing[i]; gaussianGenerator->SetSize( size ); gaussianGenerator->SetSpacing( spacing ); gaussianGenerator->SetRadius(testParameters.m_HotspotRadiusInMM); gaussianGenerator->SetNumberOfGausssians(testParameters.m_NumberOfGaussian); gaussianGenerator->AddGaussian(centerXVec, centerYVec, centerZVec, sigmaXVec, sigmaYVec, sigmaZVec, altitudeVec); gaussianGenerator->Update(); image = gaussianGenerator->GetOutput(); mitk::CastToMitkImage(image, result); return result; } /** \brief Calculates hotspot statistics for given test image and ROI parameters. Uses ImageStatisticsCalculator to find a hotspot in a defined ROI within the given image. */ - static mitk::ImageStatisticsCalculator::Statistics CalculateStatistics(mitk::Image* image, const Parameters& testParameters, unsigned int label) + static mitk::ImageStatisticsCalculator::StatisticsContainer::Pointer CalculateStatistics(mitk::Image* image, const Parameters& testParameters, unsigned int label) { - mitk::ImageStatisticsCalculator::Statistics result; + mitk::ImageStatisticsCalculator::StatisticsContainer::Pointer result; const unsigned int Dimension = 3; typedef itk::Image MaskImageType; MaskImageType::Pointer mask = MaskImageType::New(); MaskImageType::SizeType size; MaskImageType::SpacingType spacing; MaskImageType::IndexType start; mitk::ImageStatisticsCalculator::Pointer statisticsCalculator = mitk::ImageStatisticsCalculator::New(); - statisticsCalculator->SetImage(image); + statisticsCalculator->SetInputImage(image); mitk::Image::Pointer mitkMaskImage; if((testParameters.m_MaxIndexX[label] > testParameters.m_MinIndexX[label] && testParameters.m_MinIndexX[label] >= 0) && (testParameters.m_MaxIndexY[label] > testParameters.m_MinIndexY[label] && testParameters.m_MinIndexY[label] >= 0) && (testParameters.m_MaxIndexZ[label] > testParameters.m_MinIndexZ[label] && testParameters.m_MinIndexZ[label] >= 0)) { for(unsigned int i = 0; i < Dimension; ++i) { start[i] = 0; spacing[i] = testParameters.m_Spacing[i]; } size[0] = testParameters.m_ImageColumns; size[1] = testParameters.m_ImageRows; size[2] = testParameters.m_ImageSlices; MaskImageType::RegionType region; region.SetIndex(start); region.SetSize(size); mask->SetSpacing(spacing); mask->SetRegions(region); mask->Allocate(); typedef itk::ImageRegionIteratorWithIndex MaskImageIteratorType; MaskImageIteratorType maskIt(mask, region); for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt) { maskIt.Set(0); } for(unsigned int i = 0; i < testParameters.m_NumberOfLabels; ++i) { for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt) { MaskImageType::IndexType index = maskIt.GetIndex(); if((index[0] >= testParameters.m_MinIndexX[i] && index[0] <= testParameters.m_MaxIndexX[i] ) && (index[1] >= testParameters.m_MinIndexY[i] && index[1] <= testParameters.m_MaxIndexY[i] ) && (index[2] >= testParameters.m_MinIndexZ[i] && index[2] <= testParameters.m_MaxIndexZ[i] )) { maskIt.Set(testParameters.m_Label[i]); } } } MITK_DEBUG << "Masking mode has set to image"; mitk::CastToMitkImage(mask, mitkMaskImage); statisticsCalculator->SetImageMask(mitkMaskImage); statisticsCalculator->SetMaskingModeToImage(); } else { MITK_DEBUG << "Masking mode has set to none"; statisticsCalculator->SetMaskingModeToNone(); } statisticsCalculator->SetHotspotRadiusInMM(testParameters.m_HotspotRadiusInMM); statisticsCalculator->SetCalculateHotspot(true); if(testParameters.m_EntireHotspotInImage == 1) { MITK_INFO << "Hotspot must be completly inside image"; statisticsCalculator->SetHotspotMustBeCompletlyInsideImage(true); } else { MITK_INFO << "Hotspot must not be completly inside image"; statisticsCalculator->SetHotspotMustBeCompletlyInsideImage(false); } statisticsCalculator->ComputeStatistics(); result = statisticsCalculator->GetStatistics(0, label); return result; } 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::ImageStatisticsCalculator::Statistics& statistics, const Parameters& testParameters, unsigned int label) + static void ValidateHotspotStatistics(const mitk::ImageStatisticsCalculator::StatisticsContainer::Pointer 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 - ValidateStatisticsItem("Hotspot mean", statistics.GetHotspotStatistics().GetMean(), testParameters.m_HotspotMean[label], eps); - ValidateStatisticsItem("Hotspot maximum", statistics.GetHotspotStatistics().GetMax(), testParameters.m_HotspotMax[label], eps); - ValidateStatisticsItem("Hotspot minimum", statistics.GetHotspotStatistics().GetMin(), testParameters.m_HotspotMin[label], eps); + ValidateStatisticsItem("Hotspot mean", hotspotStatistics->GetMean(), testParameters.m_HotspotMean[label], eps); + ValidateStatisticsItem("Hotspot maximum", hotspotStatistics->GetMax(), testParameters.m_HotspotMax[label], eps); + ValidateStatisticsItem("Hotspot minimum", hotspotStatistics->GetMin(), 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); + // 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::ImageStatisticsCalculator::Statistics statistics = mitkImageStatisticsHotspotTestClass::CalculateStatistics(image, parameters, label); + mitk::ImageStatisticsCalculator::StatisticsContainer::Pointer 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/Testing/mitkImageStatisticsTextureAnalysisTest.cpp b/Modules/ImageStatistics/Testing/mitkImageStatisticsTextureAnalysisTest.cpp index 1907147299..d5d888d586 100644 --- a/Modules/ImageStatistics/Testing/mitkImageStatisticsTextureAnalysisTest.cpp +++ b/Modules/ImageStatistics/Testing/mitkImageStatisticsTextureAnalysisTest.cpp @@ -1,235 +1,235 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkTestingMacros.h" #include "itkImage.h" #include "mitkExtendedLabelStatisticsImageFilter.h" #include "mitkExtendedStatisticsImageFilter.h" #include "mitkNumericConstants.h" /** \section Testing of Skewness and Kurtosis * This test class is for testing the added coefficients Skewness and Kurtosis * for the mitkExtendedLabelStatisticsImageFilter (Masked Images) and for the * mitkExtendedStatisticsImageFilter (Unmasked Images). Both filter will be tested * against two pictures. */ class mitkImageStatisticsTextureAnalysisTestClass { /** * \brief Explanation of the mitkImageStatisticsTextureAnalysisTestClass test class * * this test class produce test images and masking images with the method CreatingTestImageForDifferentLabelSize. * TestInstanceFortheMaskedStatisticsFilter and TestInstanceFortheUnmaskedStatisticsFilter are the two Instances * for the filters of masking and unmasking images. * TestofSkewnessAndKurtosisForMaskedImagesand TestofSkewnessAndKurtosisForUnmaskedImages are the correlated test * for the checking the values. */ public: typedef itk::Image< int,3 >ImageType; typedef ImageType::Pointer PointerOfImage; typedef itk::ExtendedLabelStatisticsImageFilter< ImageType, ImageType > LabelStatisticsFilterType; typedef LabelStatisticsFilterType::Pointer labelStatisticsFilterPointer; typedef itk::ExtendedStatisticsImageFilter< ImageType > StatisticsFilterType; typedef StatisticsFilterType::Pointer StatisticsFilterPointer; ImageType::Pointer CreatingTestImageForDifferentLabelSize( int factorOfDividingThePicture, int bufferValue, int labelValue) { ImageType::Pointer image = ImageType :: New(); ImageType::IndexType start; ImageType::SizeType size; start[0] = 0; start[1] = 0; start[2] = 0; size[0] = 100; size[1] = 100; size[2] = 100; ImageType:: RegionType region; region.SetSize( size ); region.SetIndex( start ); image->SetRegions(region); image->Allocate(); image->FillBuffer(bufferValue); int i = 0; for(unsigned int r = 0; r < 50; r++) { for(unsigned int c = 0; c < factorOfDividingThePicture; c++) { for(unsigned int l = 0; l < 100; l++) { ImageType::IndexType pixelIndex; pixelIndex[0] = r; pixelIndex[1] = c; pixelIndex[2] = l; image->SetPixel(pixelIndex, labelValue); } } } return image; } LabelStatisticsFilterType::Pointer TestInstanceFortheMaskedStatisticsFilter(ImageType::Pointer image, ImageType::Pointer maskImage) { LabelStatisticsFilterType::Pointer labelStatisticsFilter; labelStatisticsFilter = LabelStatisticsFilterType::New(); labelStatisticsFilter->SetInput( image ); labelStatisticsFilter->UseHistogramsOn(); labelStatisticsFilter->SetHistogramParameters( 20, -10, 10); labelStatisticsFilter->SetLabelInput( maskImage ); labelStatisticsFilter->Update(); return labelStatisticsFilter; } StatisticsFilterType::Pointer TestInstanceFortheUnmaskedStatisticsFilter(ImageType::Pointer image ) { StatisticsFilterType::Pointer StatisticsFilter; StatisticsFilter = StatisticsFilterType::New(); StatisticsFilter->SetInput( image ); - StatisticsFilter->SetBinSize( 20 ); + StatisticsFilter->SetHistogramParameters( 20, -10, 10 ); StatisticsFilter->Update(); return StatisticsFilter; } //test for Skewness,Kurtosis and MPP for masked Images void TestofSkewnessKurtosisAndMPPForMaskedImages(LabelStatisticsFilterType::Pointer labelStatisticsFilter, double expectedSkewness, double expectedKurtosis, double expectedMPP) { // let's create an object of our class bool isSkewsnessLowerlimitCorrect = labelStatisticsFilter->GetSkewness( 1 )- expectedKurtosis+ std::pow(10,-3) <= expectedSkewness; bool isSkewsnessUpperlimitCorrect = labelStatisticsFilter->GetSkewness( 1 )+ expectedKurtosis+ std::pow(10,-3) >= expectedSkewness; MITK_TEST_CONDITION( isSkewsnessLowerlimitCorrect && isSkewsnessUpperlimitCorrect,"expectedSkewness: " << expectedSkewness << " actual Value: " << labelStatisticsFilter->GetSkewness( 1 ) ); bool isKurtosisUpperlimitCorrect = labelStatisticsFilter->GetKurtosis( 1 ) <= expectedKurtosis+ std::pow(10,-3); bool isKurtosisLowerlimitCorrect = expectedKurtosis- std::pow(10,-3) <= labelStatisticsFilter->GetKurtosis( 1 ); MITK_TEST_CONDITION( isKurtosisUpperlimitCorrect && isKurtosisLowerlimitCorrect,"expectedKurtosis: " << expectedKurtosis << " actual Value: " << labelStatisticsFilter->GetKurtosis( 1 ) ); MITK_TEST_CONDITION( ( expectedMPP - labelStatisticsFilter->GetMPP( 1 ) ) < 1, "expected MPP: " << expectedMPP << " actual Value: " << labelStatisticsFilter->GetMPP( 1 ) ); } //test for Entropy,Uniformity and UPP for masked Images void TestofEntropyUniformityAndUppForMaskedImages(LabelStatisticsFilterType::Pointer labelStatisticsFilter, double expectedEntropy, double expectedUniformity, double expectedUPP) { bool calculatedEntropyLowerLimit = labelStatisticsFilter->GetEntropy( 1 ) >= expectedEntropy - std::pow(10,-3); bool calculatedUniformityLowerLimit = labelStatisticsFilter->GetUniformity( 1 ) >= expectedUniformity - std::pow(10,-3); bool calculatedUppLowerLimit = labelStatisticsFilter->GetUPP( 1 ) >= expectedUPP - std::pow(10,-3); bool calculatedEntropyUpperLimit = labelStatisticsFilter->GetEntropy( 1 ) <= expectedEntropy + std::pow(10,-3); bool calculatedUniformityUpperLimit = labelStatisticsFilter->GetUniformity( 1 ) <= expectedUniformity + std::pow(10,-3); bool calculatedUppUpperLimit = labelStatisticsFilter->GetUPP( 1 ) <= expectedUPP + std::pow(10,-3); MITK_TEST_CONDITION( calculatedEntropyLowerLimit && calculatedEntropyUpperLimit, "expected Entropy: " << expectedEntropy << " actual Value: " << labelStatisticsFilter->GetEntropy( 1 ) ); MITK_TEST_CONDITION( calculatedUniformityLowerLimit && calculatedUniformityUpperLimit, "expected Uniformity: " << expectedUniformity << " actual Value: " << labelStatisticsFilter->GetUniformity( 1 ) ); MITK_TEST_CONDITION( calculatedUppLowerLimit && calculatedUppUpperLimit, "expected UPP: " << expectedUPP << " actual Value: " << labelStatisticsFilter->GetUPP( 1 ) ); } //test for Skewness,Kurtosis and MPP for unmasked Images void TestofSkewnessKurtosisAndMPPForUnmaskedImages(StatisticsFilterType::Pointer StatisticsFilter, double expectedSkewness, double expectedKurtosis, double expectedMPP) { // let's create an object of our class bool isSkewsnessLowerlimitCorrect = StatisticsFilter->GetSkewness()- expectedKurtosis+ std::pow(10,-3) <= expectedSkewness; bool isSkewsnessUpperlimitCorrect = StatisticsFilter->GetSkewness()+ expectedKurtosis+ std::pow(10,-3) >= expectedSkewness; MITK_TEST_CONDITION( isSkewsnessLowerlimitCorrect && isSkewsnessUpperlimitCorrect,"expectedSkewness: " << expectedSkewness << " actual Value: " << StatisticsFilter->GetSkewness() ); bool isKurtosisUpperlimitCorrect = StatisticsFilter->GetKurtosis() <= expectedKurtosis+ std::pow(10,-3); bool isKurtosisLowerlimitCorrect = expectedKurtosis- std::pow(10,-3) <= StatisticsFilter->GetKurtosis(); MITK_TEST_CONDITION( isKurtosisUpperlimitCorrect && isKurtosisLowerlimitCorrect,"expectedKurtosis: " << expectedKurtosis << " actual Value: " << StatisticsFilter->GetKurtosis() ); MITK_TEST_CONDITION( ( expectedMPP - StatisticsFilter->GetMPP() ) < mitk::eps, "expected MPP: " << expectedMPP << " actual Value: " << StatisticsFilter->GetMPP() ); } //test for Entropy,Uniformity and UPP for unmasked Images void TestofEntropyUniformityAndUppForUnmaskedImages(StatisticsFilterType::Pointer StatisticsFilter, double expectedEntropy, double expectedUniformity, double expectedUPP) { bool calculatedEntropyLowerLimit = StatisticsFilter->GetEntropy() >= expectedEntropy - std::pow(10,-3); bool calculatedUniformityLowerLimit = StatisticsFilter->GetUniformity() >= expectedUniformity - std::pow(10,-3); bool calculatedUppLowerLimit = StatisticsFilter->GetUPP() >= expectedUPP - std::pow(10,-3); bool calculatedEntropyUpperLimit = StatisticsFilter->GetEntropy() <= expectedEntropy + std::pow(10,-3); bool calculatedUniformityUpperLimit = StatisticsFilter->GetUniformity() <= expectedUniformity + std::pow(10,-3); bool calculatedUppUpperLimit = StatisticsFilter->GetUPP() <= expectedUPP + std::pow(10,-3); MITK_TEST_CONDITION( calculatedEntropyLowerLimit && calculatedEntropyUpperLimit, "expected Entropy: " << expectedEntropy << " actual Value: " << StatisticsFilter->GetEntropy() ); MITK_TEST_CONDITION( calculatedUniformityLowerLimit && calculatedUniformityUpperLimit, "expected Uniformity: " << expectedUniformity << " actual Value: " << StatisticsFilter->GetUniformity() ); MITK_TEST_CONDITION( calculatedUppLowerLimit && calculatedUppUpperLimit, "expected UPP: " << expectedUPP << " actual Value: " << StatisticsFilter->GetUPP() ); } }; int mitkImageStatisticsTextureAnalysisTest(int, char* []) { // always start with this! MITK_TEST_BEGIN("mitkImageStatisticsTextureAnalysisTest") mitkImageStatisticsTextureAnalysisTestClass testclassInstance; mitkImageStatisticsTextureAnalysisTestClass::PointerOfImage labelImage = testclassInstance.CreatingTestImageForDifferentLabelSize(100, 1, 1); mitkImageStatisticsTextureAnalysisTestClass::PointerOfImage image = testclassInstance.CreatingTestImageForDifferentLabelSize(100, 3, 2); mitkImageStatisticsTextureAnalysisTestClass::PointerOfImage image2 = testclassInstance.CreatingTestImageForDifferentLabelSize(50, 3, 2); //test for masked images mitkImageStatisticsTextureAnalysisTestClass::labelStatisticsFilterPointer mitkLabelFilter= testclassInstance.TestInstanceFortheMaskedStatisticsFilter( image,labelImage); testclassInstance.TestofSkewnessKurtosisAndMPPForMaskedImages(mitkLabelFilter, 0, 0.999998, 2.5); testclassInstance.TestofEntropyUniformityAndUppForMaskedImages(mitkLabelFilter, 1, 0.5, 0.5); mitkImageStatisticsTextureAnalysisTestClass::labelStatisticsFilterPointer mitkLabelFilter2= testclassInstance.TestInstanceFortheMaskedStatisticsFilter( image2,labelImage); testclassInstance.TestofSkewnessKurtosisAndMPPForMaskedImages(mitkLabelFilter2, -1.1547, 2.33333, 2.75); testclassInstance.TestofEntropyUniformityAndUppForMaskedImages(mitkLabelFilter2, 0.811278, 0.625, 0.625); //test for unmasked images mitkImageStatisticsTextureAnalysisTestClass::StatisticsFilterPointer mitkFilter= testclassInstance.TestInstanceFortheUnmaskedStatisticsFilter( image); testclassInstance.TestofSkewnessKurtosisAndMPPForUnmaskedImages(mitkFilter, 0, 0.999998, 2.5); testclassInstance.TestofEntropyUniformityAndUppForUnmaskedImages(mitkFilter, 1, 0.5, 0.5); mitkImageStatisticsTextureAnalysisTestClass::StatisticsFilterPointer mitkFilter2= testclassInstance.TestInstanceFortheUnmaskedStatisticsFilter( image2); testclassInstance.TestofSkewnessKurtosisAndMPPForUnmaskedImages(mitkFilter2, -1.1547, 2.33333, 2.75); testclassInstance.TestofEntropyUniformityAndUppForUnmaskedImages( mitkFilter2, 0.811278, 0.625, 0.625); MITK_TEST_END() } diff --git a/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp b/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp index 0386beb5e7..beaa1fb7a5 100644 --- a/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp +++ b/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp @@ -1,515 +1,576 @@ #include #include #include #include #include #include #include #include #include #include #include -#include +#include #include #include #include #include #include #include #include #include #include #include +#include +#include + #include #include "itkImageFileWriter.h" namespace mitk { void ImageStatisticsCalculator::SetInputImage(mitk::Image::Pointer image) { if (image != m_Image) { m_Image = image; m_StatisticsUpdateRequiredByTimeStep.resize(m_Image->GetTimeSteps()); m_StatisticsByTimeStep.resize(m_Image->GetTimeSteps()); this->SetAllStatisticsToUpdateRequired(); } } void ImageStatisticsCalculator::SetMask(mitk::MaskGenerator::Pointer mask) { if (mask != m_MaskGenerator) { m_MaskGenerator = mask; this->SetAllStatisticsToUpdateRequired(); } } void ImageStatisticsCalculator::SetSecondaryMask(mitk::MaskGenerator::Pointer mask) { if (mask != m_SecondaryMaskGenerator) { m_SecondaryMaskGenerator = mask; this->SetAllStatisticsToUpdateRequired(); } } void ImageStatisticsCalculator::SetNBinsForHistogramStatistics(unsigned int nBins) { if (nBins != m_nBinsForHistogramStatistics) { m_nBinsForHistogramStatistics = nBins; this->SetAllStatisticsToUpdateRequired(); + this->m_UseBinSizeOverNBins = false; + } + if (m_UseBinSizeOverNBins) + { + this->SetAllStatisticsToUpdateRequired(); + 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->SetAllStatisticsToUpdateRequired(); + this->m_UseBinSizeOverNBins = true; + } + if (!m_UseBinSizeOverNBins) + { + this->SetAllStatisticsToUpdateRequired(); + this->m_UseBinSizeOverNBins = true; + } + } + + double ImageStatisticsCalculator::GetBinSizeForHistogramStatistics() const + { + return m_binSizeForHistogramStatistics; + } + ImageStatisticsCalculator::StatisticsContainer::Pointer ImageStatisticsCalculator::GetStatistics(unsigned int timeStep, unsigned int label) { + if (timeStep >= m_StatisticsByTimeStep.size()) { - throw std::runtime_error("invalid timeStep in ImageStatisticsCalculator_v2::GetStatistics"); + mitkThrow() << "invalid timeStep in ImageStatisticsCalculator_v2::GetStatistics"; } if (m_Image.IsNull()) { - throw std::runtime_error("no image"); + mitkThrow() << "no image"; + } + + if (!m_Image->IsInitialized()) + { + mitkThrow() << "Image not initialized!"; } if (m_MaskGenerator.IsNotNull()) { m_MaskGenerator->SetTimeStep(timeStep); m_InternalMask = m_MaskGenerator->GetMask(); } if (m_SecondaryMaskGenerator.IsNotNull()) { m_SecondaryMaskGenerator->SetTimeStep(timeStep); m_SecondaryMask = m_SecondaryMaskGenerator->GetMask(); } ImageTimeSelector::Pointer imgTimeSel = ImageTimeSelector::New(); imgTimeSel->SetInput(m_Image); imgTimeSel->SetTimeNr(timeStep); imgTimeSel->UpdateLargestPossibleRegion(); m_ImageTimeSlice = imgTimeSel->GetOutput(); if (m_StatisticsUpdateRequiredByTimeStep[timeStep]) { // Calculate statistics with/without mask if (m_MaskGenerator.IsNull() && m_SecondaryMaskGenerator.IsNull()) { // 1) calculate statistics unmasked: AccessByItk_1(m_ImageTimeSlice, InternalCalculateStatisticsUnmasked, timeStep) } else { // 2) calculate statistics masked AccessByItk_1(m_ImageTimeSlice, InternalCalculateStatisticsMasked, timeStep) } m_StatisticsUpdateRequiredByTimeStep[timeStep]=false; } for (std::vector::iterator it = m_StatisticsByTimeStep[timeStep].begin(); it != m_StatisticsByTimeStep[timeStep].end(); ++it) { StatisticsContainer::Pointer statCont = *it; if (statCont->GetLabel() == label) { return statCont; } } // these lines will ony be executed if the requested label could not be found! MITK_WARN << "Invalid label: " << label << " in time step: " << timeStep; return StatisticsContainer::New(); } void ImageStatisticsCalculator::SetAllStatisticsToUpdateRequired() { for (unsigned int i = 0; i < m_StatisticsUpdateRequiredByTimeStep.size(); i++) { this->SetStatsTimeStepToUpdateRequired(i); } } void ImageStatisticsCalculator::SetStatsTimeStepToUpdateRequired(unsigned int timeStep) { if (timeStep >= m_StatisticsUpdateRequiredByTimeStep.size()) { throw std::runtime_error("invalid timeStep in ImageStatisticsCalculator_v2::SetStatsTimeStepToUpdateRequired"); } m_StatisticsUpdateRequiredByTimeStep[timeStep] = true; m_StatisticsByTimeStep[timeStep].resize(0); } template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator::InternalCalculateStatisticsUnmasked( typename itk::Image< TPixel, VImageDimension >* image, unsigned int timeStep) { typedef typename itk::Image< TPixel, VImageDimension > ImageType; typedef typename itk::ExtendedStatisticsImageFilter ImageStatisticsFilterType; typedef typename itk::MinMaxImageFilterWithIndex MinMaxFilterType; StatisticsContainer::Pointer statisticsResult = StatisticsContainer::New(); 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]; } statisticsResult->SetMinIndex(minIndex); statisticsResult->SetMaxIndex(maxIndex); - //statisticsFilter->SetHistogramParameters(m_nBinsForHistogramStatistics, imgMinMaxFilter->GetMinimum(), imgMinMaxFilter->GetMaximum()); - statisticsFilter->SetHistogramParameters(m_nBinsForHistogramStatistics, minval, maxval); + //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(); } // no mask, therefore just one label = the whole image m_StatisticsByTimeStep[timeStep].resize(1); statisticsResult->SetLabel(1); statisticsResult->SetN(image->GetLargestPossibleRegion().GetNumberOfPixels()); statisticsResult->SetMean(statisticsFilter->GetMean()); statisticsResult->SetMin(statisticsFilter->GetMinimum()); statisticsResult->SetMax(statisticsFilter->GetMaximum()); statisticsResult->SetVariance(statisticsFilter->GetVariance()); statisticsResult->SetStd(statisticsFilter->GetSigma()); statisticsResult->SetSkewness(statisticsFilter->GetSkewness()); statisticsResult->SetKurtosis(statisticsFilter->GetKurtosis()); statisticsResult->SetRMS(std::sqrt(std::pow(statisticsFilter->GetMean(), 2.) + statisticsFilter->GetVariance())); // variance = sigma^2 statisticsResult->SetMPP(statisticsFilter->GetMPP()); statisticsResult->SetEntropy(statisticsFilter->GetEntropy()); statisticsResult->SetMedian(statisticsFilter->GetMedian()); statisticsResult->SetUniformity(statisticsFilter->GetUniformity()); statisticsResult->SetUPP(statisticsFilter->GetUPP()); statisticsResult->SetHistogram(statisticsFilter->GetHistogram()); m_StatisticsByTimeStep[timeStep][0] = statisticsResult; } template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator::InternalCalculateStatisticsMasked( typename itk::Image< TPixel, VImageDimension >* image, unsigned int timeStep) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< unsigned short, VImageDimension > MaskType; typedef typename MaskType::PixelType LabelPixelType; typedef itk::ExtendedLabelStatisticsImageFilter< ImageType, MaskType > ImageStatisticsFilterType; typedef MaskUtilities< TPixel, VImageDimension > MaskUtilType; typedef typename itk::MinMaxLabelImageFilterWithIndex MinMaxLabelFilterType; typedef typename ImageType::PixelType InputImgPixelType; // workaround: if m_SecondaryMaskGenerator ist not null but m_MaskGenerator is! (this is the case if we request a 'ignore zuero valued pixels' // mask in the gui but do not define a primary mask) bool swapMasks = false; if (m_SecondaryMask.IsNotNull() && m_InternalMask.IsNull()) { m_InternalMask = m_SecondaryMask; m_SecondaryMask = nullptr; swapMasks = true; } // maskImage has to have the same dimension as image typename MaskType::Pointer maskImage = MaskType::New(); - maskImage = ImageToItkImage< unsigned short, VImageDimension >(m_InternalMask); + try { + // try to access the pixel values directly (no copying or casting). Only works if mask pixels are of pixelType unsigned short + maskImage = ImageToItkImage< unsigned short, VImageDimension >(m_InternalMask); + } + catch (itk::ExceptionObject & e) + { + // if the pixel type of the mask is not short, then we have to make a copy of m_InternalMask (and cast the values) + CastToItkImage(m_InternalMask, maskImage); + } // if we have a secondary mask (say a ignoreZeroPixelMask) we need to combine the masks (corresponds to AND) if (m_SecondaryMask.IsNotNull()) { typename MaskType::Pointer secondaryMaskImage = MaskType::New(); secondaryMaskImage = ImageToItkImage< unsigned short, VImageDimension >(m_SecondaryMask); // secondary mask should be a ignore zero value pixel mask derived from image. it has to be cropped to the mask region (which may be planar or simply smaller) typename MaskUtilities::Pointer secondaryMaskMaskUtil = MaskUtilities::New(); secondaryMaskMaskUtil->SetImage(secondaryMaskImage.GetPointer()); secondaryMaskMaskUtil->SetMask(maskImage.GetPointer()); typename MaskType::Pointer adaptedSecondaryMaskImage = secondaryMaskMaskUtil->ExtractMaskImageRegion(); typename itk::MaskImageFilter2::Pointer maskFilter = itk::MaskImageFilter2::New(); maskFilter->SetInput1(maskImage); maskFilter->SetInput2(adaptedSecondaryMaskImage); maskFilter->SetMaskingValue(1); // all pixels of maskImage where secondaryMaskImage==1 will be kept, all the others are set to 0 maskFilter->UpdateLargestPossibleRegion(); maskImage = maskFilter->GetOutput(); } 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::Pointer 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 + // set histogram parameters for each label individually (min/max may be different for each label) typedef typename std::map MapType; typedef typename std::pair PairType; std::vector relevantLabels = minMaxFilter->GetRelevantLabels(); MapType minVals; MapType maxVals; std::map nBins; for (LabelPixelType label:relevantLabels) { minVals.insert(PairType(label, minMaxFilter->GetMin(label))); maxVals.insert(PairType(label, minMaxFilter->GetMax(label))); - nBins.insert(typename std::pair(label, m_nBinsForHistogramStatistics)); - } -// minVal = minMaxFilter->GetGlobalMin(); -// maxVal = minMaxFilter->GetGlobalMax(); + unsigned int nBinsForHistogram; + if (m_UseBinSizeOverNBins) + { + nBinsForHistogram = std::max(static_cast(std::ceil(minMaxFilter->GetMax(label) - minMaxFilter->GetMin(label))) / m_binSizeForHistogramStatistics, 10.); // do not allow less than 10 bins + } + else + { + nBinsForHistogram = m_nBinsForHistogramStatistics; + } + + nBins.insert(typename std::pair(label, nBinsForHistogram)); + } typename ImageStatisticsFilterType::Pointer imageStatisticsFilter = ImageStatisticsFilterType::New(); imageStatisticsFilter->SetDirectionTolerance(0.001); imageStatisticsFilter->SetCoordinateTolerance(0.001); imageStatisticsFilter->SetInput(adaptedImage); imageStatisticsFilter->SetLabelInput(maskImage); - //imageStatisticsFilter->SetHistogramParameters(m_nBinsForHistogramStatistics, floor(minVal), ceil(maxVal)); imageStatisticsFilter->SetHistogramParametersForLabels(nBins, minVals, maxVals); imageStatisticsFilter->Update(); std::list labels = imageStatisticsFilter->GetRelevantLabels(); std::list::iterator it = labels.begin(); m_StatisticsByTimeStep[timeStep].resize(0); while(it != labels.end()) { StatisticsContainer::Pointer statisticsResult = StatisticsContainer::New(); // find min, max, minindex and maxindex // make sure to only look in the masked region, use a masker for this vnl_vector minIndex, maxIndex; typename ImageType::IndexType tmpMinIndex = minMaxFilter->GetMinIndex(*it); typename ImageType::IndexType tmpMaxIndex = minMaxFilter->GetMaxIndex(*it); minIndex.set_size(tmpMaxIndex.GetIndexDimension()); maxIndex.set_size(tmpMaxIndex.GetIndexDimension()); for (unsigned int i=0; i < tmpMaxIndex.GetIndexDimension(); i++) { minIndex[i] = tmpMinIndex[i] + (maskImage->GetOrigin()[i] - image->GetOrigin()[i]) / (double) maskImage->GetSpacing()[i]; maxIndex[i] = tmpMaxIndex[i] + (maskImage->GetOrigin()[i] - image->GetOrigin()[i]) / (double) maskImage->GetSpacing()[i]; } statisticsResult->SetMinIndex(minIndex); statisticsResult->SetMaxIndex(maxIndex); // just debug assert(minMaxFilter->GetMax(*it) == imageStatisticsFilter->GetMaximum(*it)); assert(minMaxFilter->GetMin(*it) == imageStatisticsFilter->GetMinimum(*it)); statisticsResult->SetN(imageStatisticsFilter->GetSum(*it) / (double) imageStatisticsFilter->GetMean(*it)); statisticsResult->SetMean(imageStatisticsFilter->GetMean(*it)); statisticsResult->SetMin(imageStatisticsFilter->GetMinimum(*it)); statisticsResult->SetMax(imageStatisticsFilter->GetMaximum(*it)); statisticsResult->SetVariance(imageStatisticsFilter->GetVariance(*it)); statisticsResult->SetStd(imageStatisticsFilter->GetSigma(*it)); statisticsResult->SetSkewness(imageStatisticsFilter->GetSkewness(*it)); statisticsResult->SetKurtosis(imageStatisticsFilter->GetKurtosis(*it)); statisticsResult->SetRMS(std::sqrt(std::pow(imageStatisticsFilter->GetMean(*it), 2.) + imageStatisticsFilter->GetVariance(*it))); // variance = sigma^2 statisticsResult->SetMPP(imageStatisticsFilter->GetMPP(*it)); statisticsResult->SetLabel(*it); statisticsResult->SetEntropy(imageStatisticsFilter->GetEntropy(*it)); statisticsResult->SetMedian(imageStatisticsFilter->GetMedian(*it)); statisticsResult->SetUniformity(imageStatisticsFilter->GetUniformity(*it)); statisticsResult->SetUPP(imageStatisticsFilter->GetUPP(*it)); statisticsResult->SetHistogram(imageStatisticsFilter->GetHistogram(*it)); m_StatisticsByTimeStep[timeStep].push_back(statisticsResult); ++it; } // swap maskGenerators back if (swapMasks) { m_SecondaryMask = m_InternalMask; m_InternalMask = nullptr; } } ImageStatisticsCalculator::StatisticsContainer::StatisticsContainer(): m_N(std::numeric_limits::max()), m_Mean(std::numeric_limits::max()), m_Min(std::numeric_limits::min()), m_Max(std::numeric_limits::max()), m_Std(std::numeric_limits::max()), m_Variance(std::numeric_limits::max()), m_Skewness(std::numeric_limits::max()), m_Kurtosis(std::numeric_limits::max()), m_RMS(std::numeric_limits::max()), m_MPP(std::numeric_limits::max()), m_Median(std::numeric_limits::max()), m_Uniformity(std::numeric_limits::max()), m_UPP(std::numeric_limits::max()), m_Entropy(std::numeric_limits::max()) { m_minIndex.set_size(0); m_maxIndex.set_size(0); } ImageStatisticsCalculator::statisticsMapType ImageStatisticsCalculator::StatisticsContainer::GetStatisticsAsMap() { ImageStatisticsCalculator::statisticsMapType statisticsAsMap; statisticsAsMap["N"] = m_N; statisticsAsMap["Mean"] = m_Mean; statisticsAsMap["Min"] = m_Min; statisticsAsMap["Max"] = m_Max; statisticsAsMap["StandardDeviation"] = m_Std; statisticsAsMap["Variance"] = m_Variance; statisticsAsMap["Skewness"] = m_Skewness; statisticsAsMap["Kurtosis"] = m_Kurtosis; statisticsAsMap["RMS"] = m_RMS; statisticsAsMap["MPP"] = m_MPP; statisticsAsMap["Median"] = m_Median; statisticsAsMap["Uniformity"] = m_Uniformity; statisticsAsMap["UPP"] = m_UPP; statisticsAsMap["Entropy"] = m_Entropy; return statisticsAsMap; } void ImageStatisticsCalculator::StatisticsContainer::Reset() { m_N = std::numeric_limits::max(); m_Mean = std::numeric_limits::max(); m_Min = std::numeric_limits::max(); m_Max = std::numeric_limits::max(); m_Std = std::numeric_limits::max(); m_Variance = std::numeric_limits::max(); m_Skewness = std::numeric_limits::max(); m_Kurtosis = std::numeric_limits::max(); m_RMS = std::numeric_limits::max(); m_MPP = std::numeric_limits::max(); m_Median = std::numeric_limits::max(); m_Uniformity = std::numeric_limits::max(); m_UPP = std::numeric_limits::max(); m_Entropy = std::numeric_limits::max(); } void ImageStatisticsCalculator::StatisticsContainer::Print() { ImageStatisticsCalculator::statisticsMapType statMap = this->GetStatisticsAsMap(); // print all map key value pairs // const auto& val:statMap for (auto it = statMap.begin(); it != statMap.end(); ++it) { std::cout << it->first << ": " << it->second << std::endl; } // print the min and max index std::cout << "Min Index:" << std::endl; for (auto it = this->GetMinIndex().begin(); it != this->GetMinIndex().end(); it++) { std::cout << *it << " "; } std::cout << std::endl; // print the min and max index std::cout << "Max Index:" << std::endl; for (auto it = this->GetMaxIndex().begin(); it != this->GetMaxIndex().end(); it++) { std::cout << *it << " "; } std::cout << std::endl; } std::string ImageStatisticsCalculator::StatisticsContainer::GetAsString() { std::string res = ""; ImageStatisticsCalculator::statisticsMapType statMap = this->GetStatisticsAsMap(); // print all map key value pairs // const auto& val:statMap for (auto it = statMap.begin(); it != statMap.end(); ++it) { res += std::string(it->first) + ": " + std::to_string(it->second) + "\n"; } // print the min and max index res += "Min Index:" + std::string("\n"); for (auto it = this->GetMinIndex().begin(); it != this->GetMinIndex().end(); it++) { res += std::to_string(*it) + std::string(" "); } res += "\n"; // print the min and max index res += "Max Index:" + std::string("\n"); for (auto it = this->GetMaxIndex().begin(); it != this->GetMaxIndex().end(); it++) { res += std::to_string(*it) + " "; } res += "\n"; return res; } + } diff --git a/Modules/ImageStatistics/mitkImageStatisticsCalculator.h b/Modules/ImageStatistics/mitkImageStatisticsCalculator.h index a502e815cf..9540464d4e 100644 --- a/Modules/ImageStatistics/mitkImageStatisticsCalculator.h +++ b/Modules/ImageStatistics/mitkImageStatisticsCalculator.h @@ -1,341 +1,345 @@ #ifndef MITKIMAGESTATISTICSCALCULATOR #define MITKIMAGESTATISTICSCALCULATOR #include #include #include #include #include #include #include /*#define mitkGetRefConstMacro(name, type) \ const type& Get##name() const \ { \ return &m_##name; \ } \ \ type& Get##name() \ { \ return &m_##name; \ } \ */ 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; class StatisticsContainer: public itk::Object { public: /** Standard Self typedef */ typedef StatisticsContainer 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(StatisticsContainer, itk::Object) typedef double RealType; statisticsMapType GetStatisticsAsMap(); void Reset(); void SetN(long n) { m_N = n; } const long& GetN() const { return m_N; } void SetMean(RealType mean) { m_Mean = mean; } const RealType& GetMean() const { return m_Mean; } void SetVariance(RealType variance) { m_Variance = variance; } const RealType& GetVariance() const { return m_Variance; } void SetStd(RealType std) { m_Std = std; } const RealType& GetStd() const { return m_Std; } void SetMin(RealType minVal) { m_Min = minVal; } const RealType& GetMin() const { return m_Min; } void SetMax(RealType maxVal) { m_Max = maxVal; } const RealType& GetMax() const { return m_Max; } void SetRMS(RealType rms) { m_RMS = rms; } const RealType& GetRMS() const { return m_RMS; } void SetSkewness(RealType skewness) { m_Skewness = skewness; } const RealType& GetSkewness() const { return m_Skewness; } void SetKurtosis(RealType kurtosis) { m_Kurtosis = kurtosis; } const RealType& GetKurtosis() const { return m_Kurtosis; } void SetMPP(RealType mpp) { m_MPP = mpp; } const RealType& GetMPP() const { return m_MPP; } void SetLabel(unsigned int label) { m_Label = label; } const unsigned int& GetLabel() const { return m_Label; } void SetMinIndex(vnl_vector& minIndex) { m_minIndex = minIndex; } const vnl_vector& GetMinIndex() const { return m_minIndex; } void SetMaxIndex(vnl_vector& maxIndex) { m_maxIndex = maxIndex; } const vnl_vector& GetMaxIndex() const { return m_maxIndex; } void SetHistogram(HistogramType::Pointer hist) { if (m_Histogram != hist) { m_Histogram = hist; } } const HistogramType::Pointer GetHistogram() const { return m_Histogram; } void SetEntropy(RealType entropy) { m_Entropy = entropy; } const RealType & GetEntropy() const { return m_Entropy; } void SetMedian(RealType median) { m_Median = median; } const RealType & GetMedian() const { return m_Median; } void SetUniformity(RealType uniformity) { m_Uniformity = uniformity; } const RealType & GetUniformity() const { return m_Uniformity; } void SetUPP(RealType upp) { m_UPP = upp; } const RealType & GetUPP() const { return m_UPP; } void Print(); std::string GetAsString(); protected: StatisticsContainer(); private: // not pretty, is temporary long m_N; RealType m_Mean, m_Min, m_Max, m_Std, m_Variance; RealType m_Skewness; RealType m_Kurtosis; RealType m_RMS; RealType m_MPP; vnl_vector m_minIndex, m_maxIndex; RealType m_Median; RealType m_Uniformity; RealType m_UPP; RealType m_Entropy; unsigned int m_Label; HistogramType::Pointer m_Histogram; }; void SetInputImage(mitk::Image::Pointer image); void SetMask(mitk::MaskGenerator::Pointer mask); void SetSecondaryMask(mitk::MaskGenerator::Pointer mask); void SetNBinsForHistogramStatistics(unsigned int nBins); unsigned int GetNBinsForHistogramStatistics() const; + void SetBinSizeForHistogramStatistics(double binSize); + + double GetBinSizeForHistogramStatistics() const; + StatisticsContainer::Pointer GetStatistics(unsigned int timeStep=0, unsigned int label=1); protected: ImageStatisticsCalculator(){ m_nBinsForHistogramStatistics = 100; + m_binSizeForHistogramStatistics = 10; + m_UseBinSizeOverNBins = false; }; + private: void SetAllStatisticsToUpdateRequired(); void SetStatsTimeStepToUpdateRequired(unsigned int timeStep); template < typename TPixel, unsigned int VImageDimension > void InternalCalculateStatisticsUnmasked( typename itk::Image< TPixel, VImageDimension >* image, unsigned int timeStep); template < typename TPixel, unsigned int VImageDimension > typename HistogramType::Pointer InternalCalculateHistogramUnmasked( typename itk::Image< TPixel, VImageDimension >* image, double minVal, double maxVal); template < typename TPixel, unsigned int VImageDimension > void InternalCalculateStatisticsMasked( typename itk::Image< TPixel, VImageDimension >* image, unsigned int timeStep); -// template < typename TPixel, unsigned int VImageDimension > void GetMinAndMaxValue( -// typename itk::Image< TPixel, VImageDimension >* inputImage, -// double &minVal, -// double &maxVal -// ); std::string GetNameOfClass() { return std::string("ImageStatisticsCalculator_v2"); } mitk::Image::Pointer m_Image; mitk::Image::Pointer m_ImageTimeSlice; mitk::MaskGenerator::Pointer m_MaskGenerator; mitk::Image::Pointer m_InternalMask; mitk::MaskGenerator::Pointer m_SecondaryMaskGenerator; mitk::Image::Pointer m_SecondaryMask; unsigned int m_nBinsForHistogramStatistics; + double m_binSizeForHistogramStatistics; + bool m_UseBinSizeOverNBins; std::vector m_StatisticsUpdateRequiredByTimeStep; // holds which time steps are valid and which ones have to be recalculated std::vector> m_StatisticsByTimeStep; }; } #endif // MITKIMAGESTATISTICSCALCULATOR diff --git a/Plugins/org.mitk.gui.qt.measurementtoolbox/src/internal/QmitkImageStatisticsCalculationThread.cpp b/Plugins/org.mitk.gui.qt.measurementtoolbox/src/internal/QmitkImageStatisticsCalculationThread.cpp index a163340ca4..52c5258834 100644 --- a/Plugins/org.mitk.gui.qt.measurementtoolbox/src/internal/QmitkImageStatisticsCalculationThread.cpp +++ b/Plugins/org.mitk.gui.qt.measurementtoolbox/src/internal/QmitkImageStatisticsCalculationThread.cpp @@ -1,240 +1,266 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "QmitkImageStatisticsCalculationThread.h" //QT headers #include #include #include #include #include QmitkImageStatisticsCalculationThread::QmitkImageStatisticsCalculationThread():QThread(), m_StatisticsImage(nullptr), m_BinaryMask(nullptr), m_PlanarFigureMask(nullptr), m_TimeStep(0), - m_IgnoreZeros(false), m_CalculationSuccessful(false), m_StatisticChanged(false), m_HistogramBinSize(1.0), m_UseDefaultBinSize(true) + m_IgnoreZeros(false), m_CalculationSuccessful(false), m_StatisticChanged(false), m_HistogramBinSize(10.0), m_UseDefaultNBins(true), m_nBinsForHistogramStatistics(100), m_prioritizeNBinsOverBinSize(true) { } QmitkImageStatisticsCalculationThread::~QmitkImageStatisticsCalculationThread() { } void QmitkImageStatisticsCalculationThread::Initialize( mitk::Image::Pointer image, mitk::Image::Pointer binaryImage, mitk::PlanarFigure::Pointer planarFig ) { // reset old values if( this->m_StatisticsImage.IsNotNull() ) this->m_StatisticsImage = nullptr; if( this->m_BinaryMask.IsNotNull() ) this->m_BinaryMask = nullptr; if( this->m_PlanarFigureMask.IsNotNull()) this->m_PlanarFigureMask = nullptr; // set new values if passed in if(image.IsNotNull()) this->m_StatisticsImage = image->Clone(); if(binaryImage.IsNotNull()) this->m_BinaryMask = binaryImage->Clone(); if(planarFig.IsNotNull()) this->m_PlanarFigureMask = planarFig->Clone(); } -void QmitkImageStatisticsCalculationThread::SetUseDefaultBinSize(bool useDefault) +void QmitkImageStatisticsCalculationThread::SetUseDefaultNBins(bool useDefault) { - m_UseDefaultBinSize = useDefault; + m_UseDefaultNBins = useDefault; } void QmitkImageStatisticsCalculationThread::SetTimeStep( int times ) { this->m_TimeStep = times; } int QmitkImageStatisticsCalculationThread::GetTimeStep() { return this->m_TimeStep; } std::vector QmitkImageStatisticsCalculationThread::GetStatisticsData() { return this->m_StatisticsVector; } mitk::Image::Pointer QmitkImageStatisticsCalculationThread::GetStatisticsImage() { return this->m_StatisticsImage; } void QmitkImageStatisticsCalculationThread::SetIgnoreZeroValueVoxel(bool _arg) { this->m_IgnoreZeros = _arg; } bool QmitkImageStatisticsCalculationThread::GetIgnoreZeroValueVoxel() { return this->m_IgnoreZeros; } void QmitkImageStatisticsCalculationThread::SetHistogramBinSize(double size) { this->m_HistogramBinSize = size; + this->m_prioritizeNBinsOverBinSize = false; } -double QmitkImageStatisticsCalculationThread::GetHistogramBinSize() +double QmitkImageStatisticsCalculationThread::GetHistogramBinSize() const { return this->m_HistogramBinSize; } +void QmitkImageStatisticsCalculationThread::SetHistogramNBins(double size) +{ + this->m_nBinsForHistogramStatistics = size; + this->m_prioritizeNBinsOverBinSize = true; +} + +double QmitkImageStatisticsCalculationThread::GetHistogramNBins() const +{ + return this->m_nBinsForHistogramStatistics; +} + std::string QmitkImageStatisticsCalculationThread::GetLastErrorMessage() { return m_message; } QmitkImageStatisticsCalculationThread::HistogramType::Pointer QmitkImageStatisticsCalculationThread::GetTimeStepHistogram(unsigned int t) { if (t >= this->m_HistogramVector.size()) return nullptr; return this->m_HistogramVector[t]; } bool QmitkImageStatisticsCalculationThread::GetStatisticsChangedFlag() { return m_StatisticChanged; } bool QmitkImageStatisticsCalculationThread::GetStatisticsUpdateSuccessFlag() { return m_CalculationSuccessful; } void QmitkImageStatisticsCalculationThread::run() { 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()) { mitk::ImageMaskGenerator::Pointer imgMask = mitk::ImageMaskGenerator::New(); imgMask->SetImageMask(m_BinaryMask); calculator->SetMask(imgMask.GetPointer()); } if(this->m_PlanarFigureMask.IsNotNull()) { mitk::PlanarFigureMaskGenerator::Pointer pfMaskGen = mitk::PlanarFigureMaskGenerator::New(); pfMaskGen->SetImage(m_StatisticsImage); pfMaskGen->SetPlanarFigure(m_PlanarFigureMask); calculator->SetMask(pfMaskGen.GetPointer()); } } catch( const itk::ExceptionObject& e) { MITK_ERROR << "ITK Exception:" << e.what(); statisticCalculationSuccessful = false; } catch( const mitk::Exception& e ) { MITK_ERROR<< "MITK Exception: " << e.what(); statisticCalculationSuccessful = false; } catch ( const std::runtime_error &e ) { MITK_ERROR<< "Runtime Exception: " << e.what(); statisticCalculationSuccessful = false; } catch ( const std::exception &e ) { //m_message = "Failure: " + std::string(e.what()); MITK_ERROR<< "Standard Exception: " << e.what(); statisticCalculationSuccessful = false; } bool statisticChanged = false; if (this->m_IgnoreZeros) { mitk::IgnorePixelMaskGenerator::Pointer ignorePixelValueMaskGen = mitk::IgnorePixelMaskGenerator::New(); ignorePixelValueMaskGen->SetIgnoredPixelValue(0); ignorePixelValueMaskGen->SetImage(m_StatisticsImage); calculator->SetSecondaryMask(ignorePixelValueMaskGen.GetPointer()); } else { calculator->SetSecondaryMask(nullptr); } - calculator->SetNBinsForHistogramStatistics(m_HistogramBinSize); + + if (m_UseDefaultNBins) + { + calculator->SetNBinsForHistogramStatistics(100); + } + else + { + if (!m_prioritizeNBinsOverBinSize) + { + calculator->SetBinSizeForHistogramStatistics(m_HistogramBinSize); + } + else + { + calculator->SetNBinsForHistogramStatistics(100); + } + } + //calculator->SetHistogramBinSize( m_HistogramBinSize ); //calculator->SetUseDefaultBinSize( m_UseDefaultBinSize ); for (unsigned int i = 0; i < m_StatisticsImage->GetTimeSteps(); i++) { try { calculator->GetStatistics(i); } catch ( mitk::Exception& e) { //m_message = e.GetDescription(); MITK_ERROR<< "MITK Exception: " << e.what(); statisticCalculationSuccessful = false; } catch ( const std::runtime_error &e ) { //m_message = "Failure: " + std::string(e.what()); MITK_ERROR<< "Runtime Exception: " << e.what(); statisticCalculationSuccessful = false; } catch ( const std::exception &e ) { //m_message = "Failure: " + std::string(e.what()); MITK_ERROR<< "Standard Exception: " << e.what(); statisticCalculationSuccessful = false; } } this->m_StatisticChanged = statisticChanged; this->m_CalculationSuccessful = statisticCalculationSuccessful; if(statisticCalculationSuccessful) { this->m_StatisticsVector.clear(); this->m_HistogramVector.clear(); for (unsigned int i = 0; i < m_StatisticsImage->GetTimeSteps(); i++) { this->m_StatisticsVector.push_back(calculator->GetStatistics(i)); this->m_HistogramVector.push_back((HistogramType*)this->m_StatisticsVector[i]->GetHistogram()); } } - - m_HistogramBinSize = calculator->GetNBinsForHistogramStatistics(); } diff --git a/Plugins/org.mitk.gui.qt.measurementtoolbox/src/internal/QmitkImageStatisticsCalculationThread.h b/Plugins/org.mitk.gui.qt.measurementtoolbox/src/internal/QmitkImageStatisticsCalculationThread.h index 55b8ead215..a18d6fd287 100644 --- a/Plugins/org.mitk.gui.qt.measurementtoolbox/src/internal/QmitkImageStatisticsCalculationThread.h +++ b/Plugins/org.mitk.gui.qt.measurementtoolbox/src/internal/QmitkImageStatisticsCalculationThread.h @@ -1,115 +1,123 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef QMITKIMAGESTATISTICSCALCULATIONTHREAD_H_INCLUDED #define QMITKIMAGESTATISTICSCALCULATIONTHREAD_H_INCLUDED //QT headers #include #include //mitk headers #include "mitkImage.h" #include "mitkPlanarFigure.h" #include "mitkImageStatisticsCalculator.h" // itk headers #ifndef __itkHistogram_h #include #endif /** /brief This class is executed as background thread for image statistics calculation. * Documentation: This class is derived from QThread and is intended to be used by QmitkImageStatisticsView to run the image statistics calculation in a background thread keepung the gui usable. * \ingroup Plugins/MeasurementToolbox */ class QmitkImageStatisticsCalculationThread : public QThread { Q_OBJECT public: typedef itk::Statistics::Histogram HistogramType; /*! /brief standard constructor. */ QmitkImageStatisticsCalculationThread(); /*! /brief standard destructor. */ ~QmitkImageStatisticsCalculationThread(); /*! *\brief Automatically calculate bin size to obtain 200 bins. */ - void SetUseDefaultBinSize(bool useDefault); + void SetUseDefaultNBins(bool useDefault); /*! /brief Initializes the object with necessary data. */ void Initialize( mitk::Image::Pointer image, mitk::Image::Pointer binaryImage, mitk::PlanarFigure::Pointer planarFig ); /*! /brief returns the calculated image statistics. */ std::vector GetStatisticsData(); /*! /brief */ mitk::Image::Pointer GetStatisticsImage(); /*! /brief Set the time step of the image you want to process. */ void SetTimeStep( int times ); /*! /brief Get the time step of the image you want to process. */ int GetTimeStep(); /*! /brief Set flag to ignore zero valued voxels */ void SetIgnoreZeroValueVoxel( bool _arg ); /*! /brief Get status of zero value voxel ignoring. */ bool GetIgnoreZeroValueVoxel(); /*! /brief Set bin size for histogram resolution.*/ void SetHistogramBinSize( double size); /*! /brief Get bin size for histogram resolution.*/ - double GetHistogramBinSize(); + double GetHistogramBinSize() const; + /*! + /brief Set bin size for histogram resolution.*/ + void SetHistogramNBins( double size); + /*! + /brief Get bin size for histogram resolution.*/ + double GetHistogramNBins() const; /*! /brief Returns the histogram of the currently selected time step. */ HistogramType::Pointer GetTimeStepHistogram(unsigned int t = 0); /*! /brief Returns a flag indicating if the statistics have changed during calculation */ bool GetStatisticsChangedFlag(); /*! /brief Returns a flag the indicates if the statistics are updated successfully */ bool GetStatisticsUpdateSuccessFlag(); /*! /brief Method called once the thread is executed. */ void run() override; std::string GetLastErrorMessage(); private: //member declaration mitk::Image::Pointer m_StatisticsImage; ///< member variable holds the input image for which the statistics need to be calculated. mitk::Image::Pointer m_BinaryMask; ///< member variable holds the binary mask image for segmentation image statistics calculation. mitk::PlanarFigure::Pointer m_PlanarFigureMask; ///< member variable holds the planar figure for segmentation image statistics calculation. std::vector m_StatisticsVector; ///< member variable holds the result structs. int m_TimeStep; ///< member variable holds the time step for statistics calculation bool m_IgnoreZeros; ///< member variable holds flag to indicate if zero valued voxel should be suppressed double m_HistogramBinSize; ///< member variable holds the bin size for histogram resolution. bool m_StatisticChanged; ///< flag set if statistics have changed bool m_CalculationSuccessful; ///< flag set if statistics calculation was successful std::vector m_HistogramVector; ///< member holds the histograms of all time steps. std::string m_message; - bool m_UseDefaultBinSize; + bool m_UseDefaultNBins; + unsigned int m_nBinsForHistogramStatistics; + bool m_prioritizeNBinsOverBinSize; }; #endif // QMITKIMAGESTATISTICSCALCULATIONTHREAD_H_INCLUDED 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 79b2f9f0a2..276c45d13b 100644 --- a/Plugins/org.mitk.gui.qt.measurementtoolbox/src/internal/QmitkImageStatisticsView.cpp +++ b/Plugins/org.mitk.gui.qt.measurementtoolbox/src/internal/QmitkImageStatisticsView.cpp @@ -1,1317 +1,1325 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "QmitkImageStatisticsView.h" // Qt includes #include #include #include // berry includes #include // mitk includes #include "mitkNodePredicateDataType.h" #include "mitkNodePredicateOr.h" #include "mitkPlanarFigureInteractor.h" // itk includes #include "itksys/SystemTools.hxx" #include #include #include const std::string QmitkImageStatisticsView::VIEW_ID = "org.mitk.views.imagestatistics"; const int QmitkImageStatisticsView::STAT_TABLE_BASE_HEIGHT = 180; QmitkImageStatisticsView::QmitkImageStatisticsView(QObject* /*parent*/, const char* /*name*/) : m_Controls( NULL ), m_TimeStepperAdapter( NULL ), m_SelectedImage( NULL ), m_SelectedImageMask( NULL ), m_SelectedPlanarFigure( NULL ), m_ImageObserverTag( -1 ), m_ImageMaskObserverTag( -1 ), m_PlanarFigureObserverTag( -1 ), m_TimeObserverTag( -1 ), m_CurrentStatisticsValid( false ), m_StatisticsUpdatePending( false ), m_DataNodeSelectionChanged ( false ), m_Visible(false) { this->m_CalculationThread = new QmitkImageStatisticsCalculationThread; } QmitkImageStatisticsView::~QmitkImageStatisticsView() { if ( m_SelectedImage != NULL ) m_SelectedImage->RemoveObserver( m_ImageObserverTag ); if ( m_SelectedImageMask != NULL ) m_SelectedImageMask->RemoveObserver( m_ImageMaskObserverTag ); if ( m_SelectedPlanarFigure != NULL ) m_SelectedPlanarFigure->RemoveObserver( m_PlanarFigureObserverTag ); while(this->m_CalculationThread->isRunning()) // wait until thread has finished { itksys::SystemTools::Delay(100); } delete this->m_CalculationThread; } void QmitkImageStatisticsView::CreateQtPartControl(QWidget *parent) { if (m_Controls == NULL) { m_Controls = new Ui::QmitkImageStatisticsViewControls; m_Controls->setupUi(parent); CreateConnections(); m_Controls->m_ErrorMessageLabel->hide(); m_Controls->m_StatisticsWidgetStack->setCurrentIndex( 0 ); m_Controls->m_BinSizeFrame->setVisible(false); } } void QmitkImageStatisticsView::CreateConnections() { if ( m_Controls ) { connect( (QObject*)(this->m_Controls->m_ButtonCopyHistogramToClipboard), SIGNAL(clicked()),(QObject*) this, SLOT(OnClipboardHistogramButtonClicked()) ); connect( (QObject*)(this->m_Controls->m_ButtonCopyStatisticsToClipboard), SIGNAL(clicked()),(QObject*) this, SLOT(OnClipboardStatisticsButtonClicked()) ); connect( (QObject*)(this->m_Controls->m_IgnoreZerosCheckbox), SIGNAL(clicked()),(QObject*) this, SLOT(OnIgnoreZerosCheckboxClicked()) ); connect( (QObject*) this->m_CalculationThread, SIGNAL(finished()),this, SLOT( OnThreadedStatisticsCalculationEnds()),Qt::QueuedConnection); connect( (QObject*) this, SIGNAL(StatisticsUpdate()),this, SLOT( RequestStatisticsUpdate()), Qt::QueuedConnection); connect( (QObject*) this->m_Controls->m_StatisticsTable, SIGNAL(cellDoubleClicked(int,int)),this, SLOT( JumpToCoordinates(int,int)) ); connect( (QObject*) (this->m_Controls->m_barRadioButton), SIGNAL(clicked()), (QObject*) (this->m_Controls->m_JSHistogram), SLOT(OnBarRadioButtonSelected())); connect( (QObject*) (this->m_Controls->m_lineRadioButton), SIGNAL(clicked()), (QObject*) (this->m_Controls->m_JSHistogram), SLOT(OnLineRadioButtonSelected())); connect( (QObject*) (this->m_Controls->m_HistogramBinSizeSpinbox), SIGNAL(editingFinished()), this, SLOT(OnHistogramBinSizeBoxValueChanged())); connect( (QObject*)(this->m_Controls->m_UseDefaultBinSizeBox), SIGNAL(clicked()),(QObject*) this, SLOT(OnDefaultBinSizeBoxChanged()) ); } } void QmitkImageStatisticsView::OnDefaultBinSizeBoxChanged() { - if (m_CalculationThread!=NULL) - m_Controls->m_HistogramBinSizeSpinbox->setValue(m_CalculationThread->GetHistogramBinSize()); - if (m_Controls->m_UseDefaultBinSizeBox->isChecked()) - m_Controls->m_BinSizeFrame->setVisible(false); - else - m_Controls->m_BinSizeFrame->setVisible(true); +// if (m_CalculationThread!=NULL) +// { +// m_Controls->m_HistogramBinSizeSpinbox->setValue(m_CalculationThread->GetHistogramBinSize()); +// } + + m_Controls->m_BinSizeFrame->setVisible(!m_Controls->m_UseDefaultBinSizeBox->isChecked()); + m_CalculationThread->SetUseDefaultNBins(m_Controls->m_UseDefaultBinSizeBox->isChecked()); + + this->UpdateStatistics(); } void QmitkImageStatisticsView::PartClosed(const berry::IWorkbenchPartReference::Pointer& ) { } void QmitkImageStatisticsView::OnTimeChanged(const itk::EventObject& e) { if (this->m_SelectedDataNodes.isEmpty() || this->m_SelectedImage == NULL) return; const mitk::SliceNavigationController::GeometryTimeEvent* timeEvent = dynamic_cast(&e); assert(timeEvent != NULL); unsigned int timestep = timeEvent->GetPos(); if (this->m_SelectedImage->GetTimeSteps() > 1) { for (int x = 0; x < this->m_Controls->m_StatisticsTable->columnCount(); x++) { for (int y = 0; y < this->m_Controls->m_StatisticsTable->rowCount(); y++) { QTableWidgetItem* item = this->m_Controls->m_StatisticsTable->item(y, x); if (item == NULL) break; if (x == timestep) { item->setBackgroundColor(Qt::yellow); } else { if (y % 2 == 0) item->setBackground(this->m_Controls->m_StatisticsTable->palette().base()); else item->setBackground(this->m_Controls->m_StatisticsTable->palette().alternateBase()); } } } this->m_Controls->m_StatisticsTable->viewport()->update(); } if ((this->m_SelectedImage->GetTimeSteps() == 1 && timestep == 0) || this->m_SelectedImage->GetTimeSteps() > 1) { // display histogram for selected timestep this->m_Controls->m_JSHistogram->ClearHistogram(); QmitkImageStatisticsCalculationThread::HistogramType::Pointer histogram = this->m_CalculationThread->GetTimeStepHistogram(timestep); if (histogram.IsNotNull()) { bool closedFigure = this->m_CalculationThread->GetStatisticsUpdateSuccessFlag(); if ( closedFigure ) { this->m_Controls->m_JSHistogram->ComputeHistogram(histogram.GetPointer()); } //this->m_Controls->m_JSHistogram->ComputeHistogram(histogram.GetPointer()); /*else { m_Controls->m_JSHistogram->ComputeIntensityProfile(timestep, true); }*/ // this->m_Controls->m_JSHistogram->SignalGraphChanged(); // hacky way to make sure the protected SignalGraphChanged() is called if (this->m_Controls->m_JSHistogram->GetUseLineGraph()) { this->m_Controls->m_JSHistogram->OnBarRadioButtonSelected(); this->m_Controls->m_JSHistogram->OnLineRadioButtonSelected(); } else { this->m_Controls->m_JSHistogram->OnLineRadioButtonSelected(); this->m_Controls->m_JSHistogram->OnBarRadioButtonSelected(); } } } } void QmitkImageStatisticsView::JumpToCoordinates(int row ,int col) { if(m_SelectedDataNodes.isEmpty()) { MITK_WARN("QmitkImageStatisticsView") << "No data node selected for statistics calculation." ; return; } mitk::Point3D world; if (row==4 && !m_WorldMinList.empty()) world = m_WorldMinList[col]; else if (row==3 && !m_WorldMaxList.empty()) world = m_WorldMaxList[col]; else return; mitk::IRenderWindowPart* part = this->GetRenderWindowPart(); if (part) { part->GetQmitkRenderWindow("axial")->GetSliceNavigationController()->SelectSliceByPoint(world); part->GetQmitkRenderWindow("sagittal")->GetSliceNavigationController()->SelectSliceByPoint(world); part->GetQmitkRenderWindow("coronal")->GetSliceNavigationController()->SelectSliceByPoint(world); mitk::SliceNavigationController::GeometryTimeEvent timeEvent(this->m_SelectedImage->GetTimeGeometry(), col); part->GetQmitkRenderWindow("axial")->GetSliceNavigationController()->SetGeometryTime(timeEvent); } } void QmitkImageStatisticsView::OnIgnoreZerosCheckboxClicked() { emit StatisticsUpdate(); } void QmitkImageStatisticsView::OnClipboardHistogramButtonClicked() { if ( m_CurrentStatisticsValid && !( m_SelectedPlanarFigure != NULL)) { const unsigned int t = this->GetRenderWindowPart()->GetTimeNavigationController()->GetTime()->GetPos(); typedef mitk::ImageStatisticsCalculator::HistogramType HistogramType; const HistogramType *histogram = this->m_CalculationThread->GetTimeStepHistogram(t).GetPointer(); QString clipboard( "Measurement \t Frequency\n" ); for ( HistogramType::ConstIterator it = histogram->Begin(); it != histogram->End(); ++it ) { if( m_Controls->m_HistogramBinSizeSpinbox->value() == 1.0) { clipboard = clipboard.append( "%L1 \t %L2\n" ) .arg( it.GetMeasurementVector()[0], 0, 'f', 0 ) .arg( it.GetFrequency() ); } else { clipboard = clipboard.append( "%L1 \t %L2\n" ) .arg( it.GetMeasurementVector()[0], 0, 'f', 2 ) .arg( it.GetFrequency() ); } } QApplication::clipboard()->setText( clipboard, QClipboard::Clipboard ); } // If a (non-closed) PlanarFigure is selected, display a line profile widget else if ( m_CurrentStatisticsValid && (m_SelectedPlanarFigure != NULL )) { auto intensity = m_Controls->m_JSHistogram->GetFrequency(); auto pixel = m_Controls->m_JSHistogram->GetMeasurement(); QString clipboard( "Pixel \t Intensity\n" ); auto j = pixel.begin(); for (auto i = intensity.begin(); i < intensity.end(); i++) { assert(j != pixel.end()); clipboard = clipboard.append( "%L1 \t %L2\n" ) .arg( (*j).toString()) .arg( (*i).toString()); j++; } QApplication::clipboard()->setText( clipboard, QClipboard::Clipboard ); } else { QApplication::clipboard()->clear(); } } void QmitkImageStatisticsView::OnClipboardStatisticsButtonClicked() { QLocale tempLocal; QLocale::setDefault(QLocale(QLocale::English, QLocale::UnitedStates)); if ( m_CurrentStatisticsValid && !( m_SelectedPlanarFigure != NULL)) { const std::vector &statistics = this->m_CalculationThread->GetStatisticsData(); // Set time borders for for loop ;) unsigned int startT, endT; if(this->m_Controls->m_CheckBox4dCompleteTable->checkState()==Qt::CheckState::Unchecked) { startT = this->GetRenderWindowPart()->GetTimeNavigationController()->GetTime()-> GetPos(); endT = startT+1; } else { startT = 0; endT = statistics.size(); } QVector< QVector > statisticsTable; QStringList headline; // Create Headline headline << " " << "Mean" << "Median" << "StdDev" << "RMS" << "Max" << "Min" << "NumberOfVoxels" << "Skewness" << "Kurtosis" << "Uniformity" << "Entropy" << "MPP" << "UPP" << "V [mm³]"; for(int i=0;i row; row.append(headline.at(i)); statisticsTable.append(row); } // Fill Table for(unsigned int t=startT;tGetMean()) << QString::number(statistics[t]->GetMedian()) << QString::number(statistics[t]->GetStd()) << QString::number(statistics[t]->GetRMS()) << QString::number(statistics[t]->GetMax()) << QString::number(statistics[t]->GetMin()) << QString::number(statistics[t]->GetN()) << QString::number(statistics[t]->GetSkewness()) << QString::number(statistics[t]->GetKurtosis()) << QString::number(statistics[t]->GetUniformity()) << QString::number(statistics[t]->GetEntropy()) << QString::number(statistics[t]->GetMPP()) << QString::number(statistics[t]->GetUPP()) << QString::number(m_Controls->m_StatisticsTable->item(7, 0)->data(Qt::DisplayRole).toDouble()); for(int z=0;zsetText(clipboard, QClipboard::Clipboard); } else { QApplication::clipboard()->clear(); } QLocale::setDefault(tempLocal); } void QmitkImageStatisticsView::OnSelectionChanged( berry::IWorkbenchPart::Pointer /*part*/, const QList &selectedNodes ) { if (this->m_Visible) { this->SelectionChanged( selectedNodes ); } else { this->m_DataNodeSelectionChanged = true; } } void QmitkImageStatisticsView::SelectionChanged(const QList &selectedNodes) { if( this->m_StatisticsUpdatePending ) { this->m_DataNodeSelectionChanged = true; return; // not ready for new data now! } if (selectedNodes.size() == this->m_SelectedDataNodes.size()) { int i = 0; for (; i < selectedNodes.size(); ++i) { if (selectedNodes.at(i) != this->m_SelectedDataNodes.at(i)) { break; } } // node selection did not change if (i == selectedNodes.size()) return; } //reset the feature image and image mask field m_Controls->m_SelectedFeatureImageLabel->setText("None"); m_Controls->m_SelectedMaskLabel->setText("None"); this->ReinitData(); if (selectedNodes.isEmpty()) { m_Controls->m_JSHistogram->ClearHistogram(); m_Controls->m_lineRadioButton->setEnabled(true); m_Controls->m_barRadioButton->setEnabled(true); m_Controls->m_HistogramBinSizeSpinbox->setEnabled(true); m_Controls->m_HistogramBinSizeCaptionLabel->setEnabled(true); // m_Controls->m_HistogramBinSizeLabel->setEnabled(true); m_Controls->m_InfoLabel->setText(QString("")); // m_Controls->horizontalLayout_3->setEnabled(false); m_Controls->groupBox->setEnabled(false); m_Controls->groupBox_3->setEnabled(false); } else { // m_Controls->horizontalLayout_3->setEnabled(true); m_Controls->groupBox->setEnabled(true); m_Controls->groupBox_3->setEnabled(true); } if(selectedNodes.size() == 1 || selectedNodes.size() == 2) { bool isBinary = false; selectedNodes.value(0)->GetBoolProperty("binary",isBinary); mitk::NodePredicateDataType::Pointer isLabelSet = mitk::NodePredicateDataType::New("LabelSetImage"); isBinary |= isLabelSet->CheckNode(selectedNodes.value(0)); if(isBinary) { m_Controls->m_JSHistogram->ClearHistogram(); m_Controls->m_lineRadioButton->setEnabled(true); m_Controls->m_barRadioButton->setEnabled(true); m_Controls->m_HistogramBinSizeSpinbox->setEnabled(true); m_Controls->m_HistogramBinSizeCaptionLabel->setEnabled(true); // m_Controls->m_HistogramBinSizeLabel->setEnabled(true); m_Controls->m_InfoLabel->setText(QString("")); } for (int i= 0; i< selectedNodes.size(); ++i) { this->m_SelectedDataNodes.push_back(selectedNodes.at(i)); } this->m_DataNodeSelectionChanged = false; this->m_Controls->m_ErrorMessageLabel->setText( "" ); this->m_Controls->m_ErrorMessageLabel->hide(); emit StatisticsUpdate(); } else { this->m_DataNodeSelectionChanged = false; } } void QmitkImageStatisticsView::ReinitData() { while( this->m_CalculationThread->isRunning()) // wait until thread has finished { itksys::SystemTools::Delay(100); } if(this->m_SelectedImage != NULL) { this->m_SelectedImage->RemoveObserver( this->m_ImageObserverTag); this->m_SelectedImage = NULL; } if(this->m_SelectedImageMask != NULL) { this->m_SelectedImageMask->RemoveObserver( this->m_ImageMaskObserverTag); this->m_SelectedImageMask = NULL; } if(this->m_SelectedPlanarFigure != NULL) { this->m_SelectedPlanarFigure->RemoveObserver( this->m_PlanarFigureObserverTag); this->m_SelectedPlanarFigure = NULL; } this->m_SelectedDataNodes.clear(); this->m_StatisticsUpdatePending = false; m_Controls->m_ErrorMessageLabel->setText( "" ); m_Controls->m_ErrorMessageLabel->hide(); this->InvalidateStatisticsTableView(); m_Controls->m_JSHistogram->ClearHistogram(); m_Controls->m_StatisticsWidgetStack->setCurrentIndex( 0 ); } void QmitkImageStatisticsView::OnThreadedStatisticsCalculationEnds() { std::stringstream message; message << ""; m_Controls->m_ErrorMessageLabel->setText( message.str().c_str() ); m_Controls->m_ErrorMessageLabel->hide(); this->WriteStatisticsToGUI(); } void QmitkImageStatisticsView::UpdateStatistics() { mitk::IRenderWindowPart* renderPart = this->GetRenderWindowPart(); if ( renderPart == NULL ) { this->m_StatisticsUpdatePending = false; return; } m_WorldMinList.clear(); m_WorldMaxList.clear(); // classify selected nodes mitk::NodePredicateDataType::Pointer isImage = mitk::NodePredicateDataType::New("Image"); mitk::NodePredicateDataType::Pointer isLabelSet = mitk::NodePredicateDataType::New("LabelSetImage"); mitk::NodePredicateOr::Pointer imagePredicate = mitk::NodePredicateOr::New(isImage, isLabelSet); std::string maskName = std::string(); std::string maskType = std::string(); std::string featureImageName = std::string(); unsigned int maskDimension = 0; // reset data from last run ITKCommandType::Pointer changeListener = ITKCommandType::New(); changeListener->SetCallbackFunction( this, &QmitkImageStatisticsView::SelectedDataModified ); mitk::DataNode::Pointer planarFigureNode; for( int i= 0 ; i < this->m_SelectedDataNodes.size(); ++i) { mitk::PlanarFigure::Pointer planarFig = dynamic_cast(this->m_SelectedDataNodes.at(i)->GetData()); if( imagePredicate->CheckNode(this->m_SelectedDataNodes.at(i)) ) { bool isMask = false; this->m_SelectedDataNodes.at(i)->GetPropertyValue("binary", isMask); isMask |= isLabelSet->CheckNode(this->m_SelectedDataNodes.at(i)); if( this->m_SelectedImageMask == NULL && isMask) { this->m_SelectedImageMask = dynamic_cast(this->m_SelectedDataNodes.at(i)->GetData()); this->m_ImageMaskObserverTag = this->m_SelectedImageMask->AddObserver(itk::ModifiedEvent(), changeListener); maskName = this->m_SelectedDataNodes.at(i)->GetName(); maskType = m_SelectedImageMask->GetNameOfClass(); maskDimension = 3; } else if( !isMask ) { if(this->m_SelectedImage == NULL) { this->m_SelectedImage = static_cast(this->m_SelectedDataNodes.at(i)->GetData()); this->m_ImageObserverTag = this->m_SelectedImage->AddObserver(itk::ModifiedEvent(), changeListener); } featureImageName = this->m_SelectedDataNodes.at(i)->GetName(); } } else if (planarFig.IsNotNull()) { if(this->m_SelectedPlanarFigure == NULL) { this->m_SelectedPlanarFigure = planarFig; this->m_PlanarFigureObserverTag = this->m_SelectedPlanarFigure->AddObserver(mitk::EndInteractionPlanarFigureEvent(), changeListener); maskName = this->m_SelectedDataNodes.at(i)->GetName(); maskType = this->m_SelectedPlanarFigure->GetNameOfClass(); maskDimension = 2; planarFigureNode = m_SelectedDataNodes.at(i); } } else { std::stringstream message; message << "" << "Invalid data node type!" << ""; m_Controls->m_ErrorMessageLabel->setText( message.str().c_str() ); m_Controls->m_ErrorMessageLabel->show(); } } if(maskName == "") { maskName = "None"; maskType = ""; maskDimension = 0; } if(featureImageName == "") { featureImageName = "None"; } if (m_SelectedPlanarFigure != NULL && m_SelectedImage == NULL) { mitk::DataStorage::SetOfObjects::ConstPointer parentSet = this->GetDataStorage()->GetSources(planarFigureNode); for (int i=0; iSize(); i++) { mitk::DataNode::Pointer node = parentSet->ElementAt(i); if( imagePredicate->CheckNode(node) ) { bool isMask = false; node->GetPropertyValue("binary", isMask); isMask |= isLabelSet->CheckNode(node); if( !isMask ) { if(this->m_SelectedImage == NULL) { this->m_SelectedImage = static_cast(node->GetData()); this->m_ImageObserverTag = this->m_SelectedImage->AddObserver(itk::ModifiedEvent(), changeListener); } } } } } unsigned int timeStep = renderPart->GetTimeNavigationController()->GetTime()->GetPos(); if ( m_SelectedImage != NULL && m_SelectedImage->IsInitialized()) { // Check if a the selected image is a multi-channel image. If yes, statistics // cannot be calculated currently. if ( m_SelectedImage->GetPixelType().GetNumberOfComponents() > 1 ) { std::stringstream message; message << "Multi-component images not supported."; m_Controls->m_ErrorMessageLabel->setText( message.str().c_str() ); m_Controls->m_ErrorMessageLabel->show(); this->InvalidateStatisticsTableView(); m_Controls->m_StatisticsWidgetStack->setCurrentIndex( 0 ); m_Controls->m_JSHistogram->ClearHistogram(); m_CurrentStatisticsValid = false; this->m_StatisticsUpdatePending = false; m_Controls->m_lineRadioButton->setEnabled(true); m_Controls->m_barRadioButton->setEnabled(true); m_Controls->m_HistogramBinSizeSpinbox->setEnabled(true); m_Controls->m_HistogramBinSizeCaptionLabel->setEnabled(true); // m_Controls->m_HistogramBinSizeLabel->setEnabled(true); m_Controls->m_InfoLabel->setText(QString("")); return; } std::stringstream maskLabel; maskLabel << maskName; if ( maskDimension > 0 ) { maskLabel << " [" << maskDimension << "D " << maskType << "]"; } m_Controls->m_SelectedMaskLabel->setText( maskLabel.str().c_str() ); m_Controls->m_SelectedFeatureImageLabel->setText(featureImageName.c_str()); // check time step validity if(m_SelectedImage->GetDimension() <= 3 && timeStep > m_SelectedImage->GetDimension(3)-1) { timeStep = m_SelectedImage->GetDimension(3)-1; } // Add the used mask time step to the mask label so the user knows which mask time step was used // if the image time step is bigger than the total number of mask time steps (see // ImageStatisticsCalculator::ExtractImageAndMask) if (m_SelectedImageMask != NULL) { unsigned int maskTimeStep = timeStep; if (maskTimeStep >= m_SelectedImageMask->GetTimeSteps()) { maskTimeStep = m_SelectedImageMask->GetTimeSteps() - 1; } m_Controls->m_SelectedMaskLabel->setText(m_Controls->m_SelectedMaskLabel->text() + QString(" (t=") + QString::number(maskTimeStep) + QString(")")); } //// initialize thread and trigger it this->m_CalculationThread->SetIgnoreZeroValueVoxel( m_Controls->m_IgnoreZerosCheckbox->isChecked() ); this->m_CalculationThread->Initialize( m_SelectedImage, m_SelectedImageMask, m_SelectedPlanarFigure ); this->m_CalculationThread->SetTimeStep( timeStep ); - this->m_CalculationThread->SetHistogramBinSize(m_Controls->m_HistogramBinSizeSpinbox->value()); + std::stringstream message; message << "Calculating statistics..."; m_Controls->m_ErrorMessageLabel->setText( message.str().c_str() ); m_Controls->m_ErrorMessageLabel->show(); try { // Compute statistics - this->m_CalculationThread->SetUseDefaultBinSize(m_Controls->m_UseDefaultBinSizeBox->isChecked()); + // this->m_CalculationThread->SetUseDefaultBinSize(m_Controls->m_UseDefaultBinSizeBox->isChecked()); this->m_CalculationThread->start(); } catch ( const mitk::Exception& e) { std::stringstream message; message << "" << e.GetDescription() << ""; m_Controls->m_ErrorMessageLabel->setText( message.str().c_str() ); m_Controls->m_ErrorMessageLabel->show(); this->m_StatisticsUpdatePending = false; } catch ( const std::runtime_error &e ) { // In case of exception, print error message on GUI std::stringstream message; message << "" << e.what() << ""; m_Controls->m_ErrorMessageLabel->setText( message.str().c_str() ); m_Controls->m_ErrorMessageLabel->show(); this->m_StatisticsUpdatePending = false; } 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->m_ErrorMessageLabel->setText( message.str().c_str() ); m_Controls->m_ErrorMessageLabel->show(); this->m_StatisticsUpdatePending = false; } } else { this->m_StatisticsUpdatePending = false; } } void QmitkImageStatisticsView::SelectedDataModified() { if( !m_StatisticsUpdatePending ) { emit StatisticsUpdate(); } } void QmitkImageStatisticsView::NodeRemoved(const mitk::DataNode *node) { while(this->m_CalculationThread->isRunning()) // wait until thread has finished { itksys::SystemTools::Delay(100); } if (node->GetData() == m_SelectedImage) { m_SelectedImage = NULL; } } void QmitkImageStatisticsView::RequestStatisticsUpdate() { if ( !m_StatisticsUpdatePending ) { if(this->m_DataNodeSelectionChanged) { this->SelectionChanged(this->GetCurrentSelection()); } else { this->m_StatisticsUpdatePending = true; this->UpdateStatistics(); } } if (this->GetRenderWindowPart()) this->GetRenderWindowPart()->RequestUpdate(); } void QmitkImageStatisticsView::OnHistogramBinSizeBoxValueChanged() { - this->UpdateStatistics(); + if (m_Controls->m_HistogramBinSizeSpinbox->value() != m_HistogramBinSize) + { + m_HistogramBinSize = m_Controls->m_HistogramBinSizeSpinbox->value(); + this->m_CalculationThread->SetHistogramBinSize(m_Controls->m_HistogramBinSizeSpinbox->value()); + this->UpdateStatistics(); + } } void QmitkImageStatisticsView::WriteStatisticsToGUI() { m_Controls->m_lineRadioButton->setEnabled(true); m_Controls->m_barRadioButton->setEnabled(true); m_Controls->m_HistogramBinSizeSpinbox->setEnabled(true); m_Controls->m_HistogramBinSizeCaptionLabel->setEnabled(true); // m_Controls->m_HistogramBinSizeLabel->setEnabled(true); m_Controls->m_InfoLabel->setText(QString("")); if(m_DataNodeSelectionChanged) { this->m_StatisticsUpdatePending = false; this->RequestStatisticsUpdate(); return; // stop visualization of results and calculate statistics of new selection } if ( this->m_CalculationThread->GetStatisticsUpdateSuccessFlag()) { if ( this->m_CalculationThread->GetStatisticsChangedFlag() ) { // Do not show any error messages m_Controls->m_ErrorMessageLabel->hide(); m_CurrentStatisticsValid = true; } if (m_Controls->m_barRadioButton->isChecked()) { m_Controls->m_JSHistogram->OnBarRadioButtonSelected(); } m_Controls->m_StatisticsWidgetStack->setCurrentIndex( 0 ); m_Controls->m_HistogramBinSizeSpinbox->setValue( this->m_CalculationThread->GetHistogramBinSize() ); //m_Controls->m_JSHistogram->ComputeHistogram( this->m_CalculationThread->GetTimeStepHistogram(this->m_CalculationThread->GetTimeStep()).GetPointer() ); this->FillStatisticsTableView( this->m_CalculationThread->GetStatisticsData(), this->m_CalculationThread->GetStatisticsImage()); } else { m_Controls->m_SelectedMaskLabel->setText( "None" ); m_Controls->m_ErrorMessageLabel->setText( m_CalculationThread->GetLastErrorMessage().c_str() ); m_Controls->m_ErrorMessageLabel->show(); // Clear statistics and histogram this->InvalidateStatisticsTableView(); m_Controls->m_StatisticsWidgetStack->setCurrentIndex( 0 ); //m_Controls->m_JSHistogram->clearHistogram(); m_CurrentStatisticsValid = false; // If a (non-closed) PlanarFigure is selected, display a line profile widget if ( m_SelectedPlanarFigure != NULL ) { // Check if the (closed) planar figure is out of bounds and so no image mask could be calculated--> Intensity Profile can not be calculated bool outOfBounds = false; if ( m_SelectedPlanarFigure->IsClosed() && m_SelectedImageMask == NULL) { outOfBounds = true; std::stringstream message; message << "Planar figure is on a rotated image plane or outside the image bounds."; m_Controls->m_InfoLabel->setText(message.str().c_str()); } // check whether PlanarFigure is initialized const mitk::PlaneGeometry *planarFigurePlaneGeometry = m_SelectedPlanarFigure->GetPlaneGeometry(); if ( !(planarFigurePlaneGeometry == NULL || outOfBounds)) { unsigned int timeStep = this->GetRenderWindowPart()->GetTimeNavigationController()->GetTime()->GetPos(); m_Controls->m_JSHistogram->SetImage(this->m_CalculationThread->GetStatisticsImage()); m_Controls->m_JSHistogram->SetPlanarFigure(m_SelectedPlanarFigure); m_Controls->m_JSHistogram->ComputeIntensityProfile(timeStep, true); //m_Controls->m_JSHistogram->ComputeIntensityProfile(timeStep); m_Controls->m_lineRadioButton->setEnabled(false); m_Controls->m_barRadioButton->setEnabled(false); m_Controls->m_HistogramBinSizeSpinbox->setEnabled(false); m_Controls->m_HistogramBinSizeCaptionLabel->setEnabled(false); // m_Controls->m_HistogramBinSizeLabel->setEnabled(false); this->FillLinearProfileStatisticsTableView( this->m_CalculationThread->GetStatisticsImage() ); std::stringstream message; message << "Only linegraph available for an intensity profile!"; m_Controls->m_InfoLabel->setText(message.str().c_str()); m_CurrentStatisticsValid = true; } else { // Clear statistics, histogram, and GUI this->InvalidateStatisticsTableView(); m_Controls->m_StatisticsWidgetStack->setCurrentIndex( 0 ); m_Controls->m_JSHistogram->ClearHistogram(); m_CurrentStatisticsValid = false; m_Controls->m_ErrorMessageLabel->hide(); m_Controls->m_SelectedMaskLabel->setText( "None" ); this->m_StatisticsUpdatePending = false; m_Controls->m_lineRadioButton->setEnabled(true); m_Controls->m_barRadioButton->setEnabled(true); m_Controls->m_HistogramBinSizeSpinbox->setEnabled(true); m_Controls->m_HistogramBinSizeCaptionLabel->setEnabled(true); // m_Controls->m_HistogramBinSizeLabel->setEnabled(true); if (!outOfBounds) m_Controls->m_InfoLabel->setText(QString("")); return; // Sebastian Wirkert: would suggest to remove this return, since it is an artifact of previous // code architecture. However, removing it will cause m_StatisticsUpdatePending to be set to false // in case of invalid statistics which it previously was not. } } } this->m_StatisticsUpdatePending = false; } void QmitkImageStatisticsView::FillStatisticsTableView( const std::vector &s, const mitk::Image *image ) { this->m_Controls->m_StatisticsTable->setColumnCount(image->GetTimeSteps()); this->m_Controls->m_StatisticsTable->horizontalHeader()->setVisible(image->GetTimeSteps() > 1); // Set Checkbox for complete copy of statistic table if(image->GetTimeSteps()>1) { this->m_Controls->m_CheckBox4dCompleteTable->setEnabled(true); } else { this->m_Controls->m_CheckBox4dCompleteTable->setEnabled(false); this->m_Controls->m_CheckBox4dCompleteTable->setChecked(false); } int decimals = 2; mitk::PixelType doublePix = mitk::MakeScalarPixelType< double >(); mitk::PixelType floatPix = mitk::MakeScalarPixelType< float >(); if (image->GetPixelType()==doublePix || image->GetPixelType()==floatPix) { decimals = 5; } for (unsigned int t = 0; t < image->GetTimeSteps(); t++) { this->m_Controls->m_StatisticsTable->setHorizontalHeaderItem(t, new QTableWidgetItem(QString::number(t))); if (s[t]->GetMaxIndex().size()==3) { mitk::Point3D index, max, min; index[0] = s[t]->GetMaxIndex()[0]; index[1] = s[t]->GetMaxIndex()[1]; index[2] = s[t]->GetMaxIndex()[2]; m_SelectedImage->GetGeometry()->IndexToWorld(index, max); this->m_WorldMaxList.push_back(max); index[0] = s[t]->GetMinIndex()[0]; index[1] = s[t]->GetMinIndex()[1]; index[2] = s[t]->GetMinIndex()[2]; m_SelectedImage->GetGeometry()->IndexToWorld(index, min); this->m_WorldMinList.push_back(min); } typedef mitk::ImageStatisticsCalculator::StatisticsContainer::RealType RealType; RealType maxVal = std::numeric_limits::max(); this->m_Controls->m_StatisticsTable->setItem( 0, t, new QTableWidgetItem( QString("%1").arg(s[t]->GetMean(), 0, 'f', decimals) ) ); this->m_Controls->m_StatisticsTable->setItem( 1, t, new QTableWidgetItem( QString("%1").arg(s[t]->GetMedian(), 0, 'f', decimals) ) ); this->m_Controls->m_StatisticsTable->setItem( 2, t, new QTableWidgetItem( QString("%1").arg(s[t]->GetStd(), 0, 'f', decimals) ) ); this->m_Controls->m_StatisticsTable->setItem( 3, t, new QTableWidgetItem( QString("%1").arg(s[t]->GetRMS(), 0, 'f', decimals) ) ); QString max; max.append(QString("%1").arg(s[t]->GetMax(), 0, 'f', decimals)); max += " ("; for (int i=0; iGetMaxIndex().size(); i++) { max += QString::number(s[t]->GetMaxIndex()[i]); if (iGetMaxIndex().size()-1) max += ","; } max += ")"; this->m_Controls->m_StatisticsTable->setItem( 4, t, new QTableWidgetItem( max ) ); QString min; min.append(QString("%1").arg(s[t]->GetMin(), 0, 'f', decimals)); min += " ("; for (int i=0; iGetMinIndex().size(); i++) { min += QString::number(s[t]->GetMinIndex()[i]); if (iGetMinIndex().size()-1) min += ","; } min += ")"; this->m_Controls->m_StatisticsTable->setItem( 5, t, new QTableWidgetItem( min ) ); this->m_Controls->m_StatisticsTable->setItem( 6, t, new QTableWidgetItem( QString("%1").arg(s[t]->GetN()) ) ); const mitk::BaseGeometry *geometry = image->GetGeometry(); if ( geometry != NULL ) { const mitk::Vector3D &spacing = image->GetGeometry()->GetSpacing(); double volume = spacing[0] * spacing[1] * spacing[2] * (double) s[t]->GetN(); this->m_Controls->m_StatisticsTable->setItem( 7, t, new QTableWidgetItem( QString("%1").arg(volume, 0, 'f', decimals) ) ); } else { this->m_Controls->m_StatisticsTable->setItem( 7, t, new QTableWidgetItem( "NA" ) ); } //statistics of higher order should have 5 decimal places because they used to be very small this->m_Controls->m_StatisticsTable->setItem( 8, t, new QTableWidgetItem( QString("%1").arg(s[t]->GetSkewness(), 0, 'f', 5) ) ); this->m_Controls->m_StatisticsTable->setItem( 9, t, new QTableWidgetItem( QString("%1").arg(s[t]->GetKurtosis(), 0, 'f', 5) ) ); this->m_Controls->m_StatisticsTable->setItem( 10, t, new QTableWidgetItem( QString("%1").arg(s[t]->GetUniformity(), 0, 'f', 5) ) ); this->m_Controls->m_StatisticsTable->setItem( 11, t, new QTableWidgetItem( QString("%1").arg(s[t]->GetEntropy(), 0, 'f', 5) ) ); this->m_Controls->m_StatisticsTable->setItem( 12, t, new QTableWidgetItem( QString("%1").arg(s[t]->GetMPP(), 0, 'f', decimals) ) ); this->m_Controls->m_StatisticsTable->setItem( 13, t, new QTableWidgetItem( QString("%1").arg(s[t]->GetUPP(), 0, 'f', 5) ) ); } this->m_Controls->m_StatisticsTable->resizeColumnsToContents(); int height = STAT_TABLE_BASE_HEIGHT; if (this->m_Controls->m_StatisticsTable->horizontalHeader()->isVisible()) height += this->m_Controls->m_StatisticsTable->horizontalHeader()->height(); if (this->m_Controls->m_StatisticsTable->horizontalScrollBar()->isVisible()) height += this->m_Controls->m_StatisticsTable->horizontalScrollBar()->height(); this->m_Controls->m_StatisticsTable->setMinimumHeight(height); // make sure the current timestep's column is highlighted (and the correct histogram is displayed) unsigned int t = this->GetRenderWindowPart()->GetTimeNavigationController()->GetTime()-> GetPos(); mitk::SliceNavigationController::GeometryTimeEvent timeEvent(this->m_SelectedImage->GetTimeGeometry(), t); this->OnTimeChanged(timeEvent); t = std::min(image->GetTimeSteps() - 1, t); // See bug 18340 /*QString hotspotMean; hotspotMean.append(QString("%1").arg(s[t].GetHotspotStatistics().GetMean(), 0, 'f', decimals)); hotspotMean += " ("; for (int i=0; im_Controls->m_StatisticsTable->setItem( 7, t, new QTableWidgetItem( hotspotMean ) ); QString hotspotMax; hotspotMax.append(QString("%1").arg(s[t].GetHotspotStatistics().GetMax(), 0, 'f', decimals)); hotspotMax += " ("; for (int i=0; im_Controls->m_StatisticsTable->setItem( 8, t, new QTableWidgetItem( hotspotMax ) ); QString hotspotMin; hotspotMin.append(QString("%1").arg(s[t].GetHotspotStatistics().GetMin(), 0, 'f', decimals)); hotspotMin += " ("; for (int i=0; im_Controls->m_StatisticsTable->setItem( 9, t, new QTableWidgetItem( hotspotMin ) );*/ } std::vector QmitkImageStatisticsView::CalculateStatisticsForPlanarFigure( const mitk::Image *image) { std::vector result; int decimals = 2; mitk::PixelType doublePix = mitk::MakeScalarPixelType< double >(); mitk::PixelType floatPix = mitk::MakeScalarPixelType< float >(); if (image->GetPixelType()==doublePix || image->GetPixelType()==floatPix) { decimals = 5; } mitk::ImageStatisticsCalculator::StatisticsContainer::Pointer stats = m_Controls->m_JSHistogram->GetStatistics(); typedef mitk::ImageStatisticsCalculator::StatisticsContainer::RealType RealType; RealType maxVal = std::numeric_limits::max(); if (stats->GetMean() == maxVal) { result.push_back(QString("NA")); } else { result.push_back(QString("%1").arg(stats->GetMean(), 0, 'f', decimals)); } if (stats->GetMedian() == maxVal) { result.push_back(QString("NA")); } else { result.push_back(QString("%1").arg(stats->GetMedian(), 0, 'f', decimals)); } if (stats->GetStd() == maxVal) { result.push_back(QString("NA")); } else { result.push_back( QString("%1").arg( stats->GetStd(), 0, 'f', decimals)); } if (stats->GetRMS() == maxVal) { result.push_back(QString("NA")); } else { result.push_back(QString("%1").arg( stats->GetRMS(), 0, 'f', decimals)); } if (stats->GetMax() == maxVal) { result.push_back(QString("NA")); } else { QString max; max.append(QString("%1").arg(stats->GetMax(), 0, 'f', decimals)); result.push_back(max); } if (stats->GetMin() == maxVal) { result.push_back(QString("NA")); } else { QString min; min.append(QString("%1").arg(stats->GetMin(), 0, 'f', decimals)); result.push_back(min); } if (stats->GetN() == maxVal) { result.push_back(QString("NA")); } else { result.push_back(QString("%1").arg(stats->GetN())); } result.push_back(QString("NA")); //statistics of higher order should have 5 decimal places because they used to be very small if (stats->GetSkewness() == maxVal) { result.push_back(QString("NA")); } else { result.push_back(QString("%1").arg(stats->GetSkewness(), 0, 'f', 5 )); } if (stats->GetKurtosis() == maxVal) { result.push_back(QString("NA")); } else { result.push_back(QString("%1").arg(stats->GetKurtosis(), 0, 'f', 5) ); } if (stats->GetUniformity() == maxVal) { result.push_back(QString("NA")); } else { result.push_back(QString("%1").arg(stats->GetUniformity(), 0, 'f', 5) ); } if (stats->GetEntropy() == maxVal) { result.push_back(QString("NA")); } else { result.push_back(QString("%1").arg(stats->GetEntropy(), 0, 'f', 5) ); } if (stats->GetMPP() == maxVal) { result.push_back(QString("NA")); } else { result.push_back(QString("%1").arg(stats->GetMPP(), 0, 'f', decimals) ); } if (stats->GetUPP() == maxVal) { result.push_back(QString("NA")); } else { result.push_back(QString("%1").arg(stats->GetUPP(), 0, 'f', 5) ); } return result; } void QmitkImageStatisticsView::FillLinearProfileStatisticsTableView( const mitk::Image *image ) { this->m_Controls->m_StatisticsTable->setColumnCount(1); this->m_Controls->m_StatisticsTable->horizontalHeader()->setVisible(false); m_PlanarFigureStatistics = this->CalculateStatisticsForPlanarFigure(image); for (int i = 0; i< m_PlanarFigureStatistics.size(); i++) { this->m_Controls->m_StatisticsTable->setItem( i, 0, new QTableWidgetItem(m_PlanarFigureStatistics[i] )); } this->m_Controls->m_StatisticsTable->resizeColumnsToContents(); int height = STAT_TABLE_BASE_HEIGHT; if (this->m_Controls->m_StatisticsTable->horizontalHeader()->isVisible()) height += this->m_Controls->m_StatisticsTable->horizontalHeader()->height(); if (this->m_Controls->m_StatisticsTable->horizontalScrollBar()->isVisible()) height += this->m_Controls->m_StatisticsTable->horizontalScrollBar()->height(); this->m_Controls->m_StatisticsTable->setMinimumHeight(height); } void QmitkImageStatisticsView::InvalidateStatisticsTableView() { this->m_Controls->m_StatisticsTable->horizontalHeader()->setVisible(false); this->m_Controls->m_StatisticsTable->setColumnCount(1); for ( unsigned int i = 0; i < this->m_Controls->m_StatisticsTable->rowCount(); ++i ) { { this->m_Controls->m_StatisticsTable->setItem( i, 0, new QTableWidgetItem( "NA" ) ); } } this->m_Controls->m_StatisticsTable->setMinimumHeight(STAT_TABLE_BASE_HEIGHT); } void QmitkImageStatisticsView::Activated() { } void QmitkImageStatisticsView::Deactivated() { } void QmitkImageStatisticsView::Visible() { m_Visible = true; mitk::IRenderWindowPart* renderWindow = GetRenderWindowPart(); if (renderWindow) { itk::ReceptorMemberCommand::Pointer cmdTimeEvent = itk::ReceptorMemberCommand::New(); cmdTimeEvent->SetCallbackFunction(this, &QmitkImageStatisticsView::OnTimeChanged); // It is sufficient to add the observer to the axial render window since the GeometryTimeEvent // is always triggered by all views. m_TimeObserverTag = renderWindow->GetQmitkRenderWindow("axial")-> GetSliceNavigationController()-> AddObserver(mitk::SliceNavigationController::GeometryTimeEvent(NULL, 0), cmdTimeEvent); } if (m_DataNodeSelectionChanged) { if (this->IsCurrentSelectionValid()) { this->SelectionChanged(this->GetCurrentSelection()); } else { this->SelectionChanged(this->GetDataManagerSelection()); } m_DataNodeSelectionChanged = false; } } void QmitkImageStatisticsView::Hidden() { m_Visible = false; // The slice navigation controller observer is removed here instead of in the destructor. // If it was called in the destructor, the application would freeze because the view's // destructor gets called after the render windows have been destructed. if ( m_TimeObserverTag != NULL ) { mitk::IRenderWindowPart* renderWindow = GetRenderWindowPart(); if (renderWindow) { renderWindow->GetQmitkRenderWindow("axial")->GetSliceNavigationController()-> RemoveObserver( m_TimeObserverTag ); } m_TimeObserverTag = NULL; } } void QmitkImageStatisticsView::SetFocus() { } diff --git a/Plugins/org.mitk.gui.qt.measurementtoolbox/src/internal/QmitkImageStatisticsView.h b/Plugins/org.mitk.gui.qt.measurementtoolbox/src/internal/QmitkImageStatisticsView.h index 7da9a2eb3c..0fb88d3013 100644 --- a/Plugins/org.mitk.gui.qt.measurementtoolbox/src/internal/QmitkImageStatisticsView.h +++ b/Plugins/org.mitk.gui.qt.measurementtoolbox/src/internal/QmitkImageStatisticsView.h @@ -1,189 +1,191 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef QmitkImageStatisticsView_H__INCLUDED #define QmitkImageStatisticsView_H__INCLUDED #include "ui_QmitkImageStatisticsViewControls.h" // Qmitk includes #include #include "QmitkStepperAdapter.h" #include "QmitkImageStatisticsCalculationThread.h" #include // mitk includes #include "mitkImageStatisticsCalculator.h" #include "mitkILifecycleAwarePart.h" #include "mitkPlanarLine.h" /*! \brief QmitkImageStatisticsView is a bundle that allows statistics calculation from images. Three modes are supported: 1. Statistics of one image, 2. Statistics of an image and a segmentation, 3. Statistics of an image and a Planar Figure. The statistics calculation is realized in a seperate thread to keep the gui accessable during calculation. \ingroup Plugins/org.mitk.gui.qt.measurementtoolbox */ class QmitkImageStatisticsView : public QmitkAbstractView, public mitk::ILifecycleAwarePart, public berry::IPartListener { Q_OBJECT private: /*! \ Convenient typedefs */ typedef mitk::DataStorage::SetOfObjects ConstVector; typedef ConstVector::ConstPointer ConstVectorPointer; typedef ConstVector::ConstIterator ConstVectorIterator; typedef std::map< mitk::Image *, mitk::ImageStatisticsCalculator::Pointer > ImageStatisticsMapType; typedef QList SelectedDataNodeVectorType; typedef itk::SimpleMemberCommand< QmitkImageStatisticsView > ITKCommandType; public: /*! \brief default constructor */ QmitkImageStatisticsView(QObject *parent=nullptr, const char *name=nullptr); /*! \brief default destructor */ virtual ~QmitkImageStatisticsView(); /*! \brief method for creating the widget containing the application controls, like sliders, buttons etc. */ virtual void CreateQtPartControl(QWidget *parent) override; /*! \brief method for creating the connections of main and control widget */ virtual void CreateConnections(); /*! \brief not implemented*/ //bool IsExclusiveFunctionality() const; /*! \brief Is called from the selection mechanism once the data manager selection has changed*/ void OnSelectionChanged( berry::IWorkbenchPart::Pointer part, const QList &nodes ) override; static const std::string VIEW_ID; static const int STAT_TABLE_BASE_HEIGHT; public slots: /** \brief Called when the statistics update is finished, sets the results to GUI.*/ void OnThreadedStatisticsCalculationEnds(); /** \brief Update bin size for histogram resolution. */ void OnHistogramBinSizeBoxValueChanged(); protected slots: /** \brief Saves the histogram to the clipboard */ void OnClipboardHistogramButtonClicked(); /** \brief Saves the statistics to the clipboard */ void OnClipboardStatisticsButtonClicked(); /** \brief Indicates if zeros should be excluded from statistics calculation */ void OnIgnoreZerosCheckboxClicked( ); /** \brief Checks if update is possible and calls StatisticsUpdate() possible */ void RequestStatisticsUpdate(); /** \brief Jump to coordinates stored in the double clicked cell */ void JumpToCoordinates(int row, int col); /** \brief Toogle GUI elements if histogram default bin size checkbox value changed. */ void OnDefaultBinSizeBoxChanged(); signals: /** \brief Method to set the data to the member and start the threaded statistics update */ void StatisticsUpdate(); protected: /** \brief Writes the calculated statistics to the GUI */ void FillStatisticsTableView(const std::vector &s, const mitk::Image *image ); std::vector CalculateStatisticsForPlanarFigure( const mitk::Image *image); void FillLinearProfileStatisticsTableView( const mitk::Image *image ); /** \brief Removes statistics from the GUI */ void InvalidateStatisticsTableView(); /** \brief Recalculate statistics for currently selected image and mask and * update the GUI. */ void UpdateStatistics(); /** \brief Listener for progress events to update progress bar. */ void UpdateProgressBar(); /** \brief Removes any cached images which are no longer referenced elsewhere. */ void RemoveOrphanImages(); /** \brief Computes an Intensity Profile along line and updates the histogram widget with it. */ void ComputeIntensityProfile( mitk::PlanarLine* line ); /** \brief Removes all Observers to images, masks and planar figures and sets corresponding members to zero */ void ClearObservers(); void Activated() override; void Deactivated() override; void Visible() override; void Hidden() override; void SetFocus() override; /** \brief Method called when itkModifiedEvent is called by selected data. */ void SelectedDataModified(); /** \brief Method called when the data manager selection changes */ void SelectionChanged(const QList &selectedNodes); /** \brief Method called to remove old selection when a new selection is present */ void ReinitData(); /** \brief writes the statistics to the gui*/ void WriteStatisticsToGUI(); void NodeRemoved(const mitk::DataNode *node) override; /** \brief Is called right before the view closes (before the destructor) */ virtual void PartClosed(const berry::IWorkbenchPartReference::Pointer& ) override; /** \brief Is called from the image navigator once the time step has changed */ void OnTimeChanged( const itk::EventObject& ); /** \brief Required for berry::IPartListener */ virtual Events::Types GetPartEventTypes() const override { return Events::CLOSED; } // member variables Ui::QmitkImageStatisticsViewControls *m_Controls; // if you have a planar figure selected, the statistics values will be saved in this one. std::vector m_PlanarFigureStatistics; QmitkImageStatisticsCalculationThread* m_CalculationThread; QmitkStepperAdapter* m_TimeStepperAdapter; unsigned int m_CurrentTime; QString m_Clipboard; // Image and mask data mitk::Image* m_SelectedImage; mitk::Image* m_SelectedImageMask; mitk::PlanarFigure* m_SelectedPlanarFigure; // observer tags long m_ImageObserverTag; long m_ImageMaskObserverTag; long m_PlanarFigureObserverTag; long m_TimeObserverTag; SelectedDataNodeVectorType m_SelectedDataNodes; bool m_CurrentStatisticsValid; bool m_StatisticsUpdatePending; bool m_StatisticsIntegrationPending; bool m_DataNodeSelectionChanged; bool m_Visible; + double m_HistogramBinSize; + std::vector m_WorldMinList; std::vector m_WorldMaxList; }; #endif // QmitkImageStatisticsView_H__INCLUDED