diff --git a/Modules/ImageStatistics/Testing/mitkImageStatisticsCalculatorTest.cpp b/Modules/ImageStatistics/Testing/mitkImageStatisticsCalculatorTest.cpp index f16f0683b3..db48266dc9 100644 --- a/Modules/ImageStatistics/Testing/mitkImageStatisticsCalculatorTest.cpp +++ b/Modules/ImageStatistics/Testing/mitkImageStatisticsCalculatorTest.cpp @@ -1,722 +1,623 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkStandardFileLocations.h" #include "mitkDicomSeriesReader.h" #include "mitkTestingMacros.h" #include "mitkImageStatisticsCalculator.h" #include "mitkPlanarPolygon.h" #include "mitkDicomSeriesReader.h" #include #include "vtkStreamingDemandDrivenPipeline.h" #include #include #include "itkImage.h" -#include +#include #include #include //#include /** * \brief Test class for mitkImageStatisticsCalculator * * This test covers: * - instantiation of an ImageStatisticsCalculator class * - correctness of statistics when using PlanarFigures for masking */ class mitkImageStatisticsCalculatorTestClass { public: struct testCase { int id; mitk::PlanarFigure::Pointer figure; double mean; double sd; }; // calculate statistics for the given image and planarpolygon static const mitk::ImageStatisticsCalculator::Statistics TestStatistics( mitk::Image::Pointer image, mitk::PlanarFigure::Pointer polygon ) { mitk::ImageStatisticsCalculator::Pointer statisticsCalculator = mitk::ImageStatisticsCalculator::New(); statisticsCalculator->SetImage( image ); statisticsCalculator->SetMaskingModeToPlanarFigure(); statisticsCalculator->SetPlanarFigure( polygon ); statisticsCalculator->ComputeStatistics(); return statisticsCalculator->GetStatistics(); } + // calculate statistics for the given image and mask + static const mitk::ImageStatisticsCalculator::Statistics TestStatisticsWithMask(mitk::Image::Pointer image, mitk::Image::Pointer mask) + { + mitk::ImageStatisticsCalculator::Pointer statisticsCalculator = mitk::ImageStatisticsCalculator::New(); + statisticsCalculator->SetImage(image); + statisticsCalculator->SetMaskingModeToImage(); + statisticsCalculator->SetImageMask(mask); + statisticsCalculator->ComputeStatistics(); + return statisticsCalculator->GetStatistics(); + } // returns a vector of defined test-cases static std::vector InitializeTestCases( mitk::Geometry2D::Pointer geom ) { std::vector testCases; { /***************************** * one whole white pixel * -> mean of 255 expected ******************************/ mitk::PlanarPolygon::Pointer figure1 = mitk::PlanarPolygon::New(); figure1->SetGeometry2D( geom ); mitk::Point2D pnt1; pnt1[0] = 10.5 ; pnt1[1] = 3.5; figure1->PlaceFigure( pnt1 ); mitk::Point2D pnt2; pnt2[0] = 9.5; pnt2[1] = 3.5; figure1->SetControlPoint( 1, pnt2, true ); mitk::Point2D pnt3; pnt3[0] = 9.5; pnt3[1] = 4.5; figure1->SetControlPoint( 2, pnt3, true ); mitk::Point2D pnt4; pnt4[0] = 10.5; pnt4[1] = 4.5; figure1->SetControlPoint( 3, pnt4, true ); figure1->GetPolyLine(0); testCase test; test.id = testCases.size(); test.figure = figure1; test.mean = 255.0; test.sd = 0.0; testCases.push_back( test ); } { /***************************** * half pixel in x-direction (white) * -> mean of 255 expected ******************************/ mitk::PlanarPolygon::Pointer figure1 = mitk::PlanarPolygon::New(); figure1->SetGeometry2D( geom ); mitk::Point2D pnt1; pnt1[0] = 10.0 ; pnt1[1] = 3.5; figure1->PlaceFigure( pnt1 ); mitk::Point2D pnt2; pnt2[0] = 9.5; pnt2[1] = 3.5; figure1->SetControlPoint( 1, pnt2, true ); mitk::Point2D pnt3; pnt3[0] = 9.5; pnt3[1] = 4.5; figure1->SetControlPoint( 2, pnt3, true ); mitk::Point2D pnt4; pnt4[0] = 10.0; pnt4[1] = 4.5; figure1->SetControlPoint( 3, pnt4, true ); figure1->GetPolyLine(0); testCase test; test.id = testCases.size(); test.figure = figure1; test.mean = 255.0; test.sd = 0.0; testCases.push_back( test ); } { /***************************** * half pixel in diagonal-direction (white) * -> mean of 255 expected ******************************/ mitk::PlanarPolygon::Pointer figure1 = mitk::PlanarPolygon::New(); figure1->SetGeometry2D( geom ); mitk::Point2D pnt1; pnt1[0] = 10.5 ; pnt1[1] = 3.5; figure1->PlaceFigure( pnt1 ); mitk::Point2D pnt2; pnt2[0] = 9.5; pnt2[1] = 3.5; figure1->SetControlPoint( 1, pnt2, true ); mitk::Point2D pnt3; pnt3[0] = 9.5; pnt3[1] = 4.5; figure1->SetControlPoint( 2, pnt3, true ); figure1->GetPolyLine(0); testCase test; test.id = testCases.size(); test.figure = figure1; test.mean = 255.0; test.sd = 0.0; testCases.push_back( test ); } { /***************************** * one pixel (white) + 2 half pixels (white) + 1 half pixel (black) * -> mean of 191.25 expected ******************************/ mitk::PlanarPolygon::Pointer figure1 = mitk::PlanarPolygon::New(); figure1->SetGeometry2D( geom ); mitk::Point2D pnt1; pnt1[0] = 1.1; pnt1[1] = 1.1; figure1->PlaceFigure( pnt1 ); mitk::Point2D pnt2; pnt2[0] = 2.0; pnt2[1] = 2.0; figure1->SetControlPoint( 1, pnt2, true ); mitk::Point2D pnt3; pnt3[0] = 3.0; pnt3[1] = 1.0; figure1->SetControlPoint( 2, pnt3, true ); mitk::Point2D pnt4; pnt4[0] = 2.0; pnt4[1] = 0.0; figure1->SetControlPoint( 3, pnt4, true ); figure1->GetPolyLine(0); testCase test; test.id = testCases.size(); test.figure = figure1; test.mean = 191.25; test.sd = 127.5; testCases.push_back( test ); } { /***************************** * whole pixel (white) + half pixel (gray) in x-direction * -> mean of 191.5 expected ******************************/ mitk::PlanarPolygon::Pointer figure1 = mitk::PlanarPolygon::New(); figure1->SetGeometry2D( geom ); mitk::Point2D pnt1; pnt1[0] = 11.0; pnt1[1] = 3.5; figure1->PlaceFigure( pnt1 ); mitk::Point2D pnt2; pnt2[0] = 9.5; pnt2[1] = 3.5; figure1->SetControlPoint( 1, pnt2, true ); mitk::Point2D pnt3; pnt3[0] = 9.5; pnt3[1] = 4.5; figure1->SetControlPoint( 2, pnt3, true ); mitk::Point2D pnt4; pnt4[0] = 11.0; pnt4[1] = 4.5; figure1->SetControlPoint( 3, pnt4, true ); figure1->GetPolyLine(0); testCase test; test.id = testCases.size(); test.figure = figure1; test.mean = 191.50; test.sd = 89.80; testCases.push_back( test ); } { /***************************** * quarter pixel (black) + whole pixel (white) + half pixel (gray) in x-direction * -> mean of 191.5 expected ******************************/ mitk::PlanarPolygon::Pointer figure1 = mitk::PlanarPolygon::New(); figure1->SetGeometry2D( geom ); mitk::Point2D pnt1; pnt1[0] = 11.0; pnt1[1] = 3.5; figure1->PlaceFigure( pnt1 ); mitk::Point2D pnt2; pnt2[0] = 9.25; pnt2[1] = 3.5; figure1->SetControlPoint( 1, pnt2, true ); mitk::Point2D pnt3; pnt3[0] = 9.25; pnt3[1] = 4.5; figure1->SetControlPoint( 2, pnt3, true ); mitk::Point2D pnt4; pnt4[0] = 11.0; pnt4[1] = 4.5; figure1->SetControlPoint( 3, pnt4, true ); figure1->GetPolyLine(0); testCase test; test.id = testCases.size(); test.figure = figure1; test.mean = 191.5; test.sd = 89.80; testCases.push_back( test ); } { /***************************** * half pixel (black) + whole pixel (white) + half pixel (gray) in x-direction * -> mean of 127.66 expected ******************************/ mitk::PlanarPolygon::Pointer figure1 = mitk::PlanarPolygon::New(); figure1->SetGeometry2D( geom ); mitk::Point2D pnt1; pnt1[0] = 11.0; pnt1[1] = 3.5; figure1->PlaceFigure( pnt1 ); mitk::Point2D pnt2; pnt2[0] = 9.0; pnt2[1] = 3.5; figure1->SetControlPoint( 1, pnt2, true ); mitk::Point2D pnt3; pnt3[0] = 9.0; pnt3[1] = 4.0; figure1->SetControlPoint( 2, pnt3, true ); mitk::Point2D pnt4; pnt4[0] = 11.0; pnt4[1] = 4.0; figure1->SetControlPoint( 3, pnt4, true ); figure1->GetPolyLine(0); testCase test; test.id = testCases.size(); test.figure = figure1; test.mean = 127.66; test.sd = 127.5; testCases.push_back( test ); } { /***************************** * whole pixel (gray) * -> mean of 128 expected ******************************/ mitk::PlanarPolygon::Pointer figure2 = mitk::PlanarPolygon::New(); figure2->SetGeometry2D( geom ); mitk::Point2D pnt1; pnt1[0] = 11.5; pnt1[1] = 10.5; figure2->PlaceFigure( pnt1 ); mitk::Point2D pnt2; pnt2[0] = 11.5; pnt2[1] = 11.5; figure2->SetControlPoint( 1, pnt2, true ); mitk::Point2D pnt3; pnt3[0] = 12.5; pnt3[1] = 11.5; figure2->SetControlPoint( 2, pnt3, true ); mitk::Point2D pnt4; pnt4[0] = 12.5; pnt4[1] = 10.5; figure2->SetControlPoint( 3, pnt4, true ); figure2->GetPolyLine(0); testCase test; test.id = testCases.size(); test.figure = figure2; test.mean = 128.0; test.sd = 0.0; testCases.push_back( test ); } { /***************************** * whole pixel (gray) + half pixel (white) in y-direction * -> mean of 191.5 expected ******************************/ mitk::PlanarPolygon::Pointer figure2 = mitk::PlanarPolygon::New(); figure2->SetGeometry2D( geom ); mitk::Point2D pnt1; pnt1[0] = 11.5; pnt1[1] = 10.5; figure2->PlaceFigure( pnt1 ); mitk::Point2D pnt2; pnt2[0] = 11.5; pnt2[1] = 12.0; figure2->SetControlPoint( 1, pnt2, true ); mitk::Point2D pnt3; pnt3[0] = 12.5; pnt3[1] = 12.0; figure2->SetControlPoint( 2, pnt3, true ); mitk::Point2D pnt4; pnt4[0] = 12.5; pnt4[1] = 10.5; figure2->SetControlPoint( 3, pnt4, true ); figure2->GetPolyLine(0); testCase test; test.id = testCases.size(); test.figure = figure2; test.mean = 191.5; test.sd = 89.80; testCases.push_back( test ); } { /***************************** * 2 whole pixel (white) + 2 whole pixel (black) in y-direction * -> mean of 127.66 expected ******************************/ mitk::PlanarPolygon::Pointer figure2 = mitk::PlanarPolygon::New(); figure2->SetGeometry2D( geom ); mitk::Point2D pnt1; pnt1[0] = 11.5; pnt1[1] = 10.5; figure2->PlaceFigure( pnt1 ); mitk::Point2D pnt2; pnt2[0] = 11.5; pnt2[1] = 13.5; figure2->SetControlPoint( 1, pnt2, true ); mitk::Point2D pnt3; pnt3[0] = 12.5; pnt3[1] = 13.5; figure2->SetControlPoint( 2, pnt3, true ); mitk::Point2D pnt4; pnt4[0] = 12.5; pnt4[1] = 10.5; figure2->SetControlPoint( 3, pnt4, true ); figure2->GetPolyLine(0); testCase test; test.id = testCases.size(); test.figure = figure2; test.mean = 127.66; test.sd = 127.5; testCases.push_back( test ); } { /***************************** * 9 whole pixels (white) + 3 half pixels (white) * + 3 whole pixel (black) [ + 3 slightly less than half pixels (black)] * -> mean of 204.0 expected ******************************/ mitk::PlanarPolygon::Pointer figure2 = mitk::PlanarPolygon::New(); figure2->SetGeometry2D( geom ); mitk::Point2D pnt1; pnt1[0] = 0.5; pnt1[1] = 0.5; figure2->PlaceFigure( pnt1 ); mitk::Point2D pnt2; pnt2[0] = 3.5; pnt2[1] = 3.5; figure2->SetControlPoint( 1, pnt2, true ); mitk::Point2D pnt3; pnt3[0] = 8.4999; pnt3[1] = 3.5; figure2->SetControlPoint( 2, pnt3, true ); mitk::Point2D pnt4; pnt4[0] = 5.4999; pnt4[1] = 0.5; figure2->SetControlPoint( 3, pnt4, true ); figure2->GetPolyLine(0); testCase test; test.id = testCases.size(); test.figure = figure2; test.mean = 204.0; test.sd = 105.58; testCases.push_back( test ); } { /***************************** * half pixel (white) + whole pixel (white) + half pixel (black) * -> mean of 212.66 expected ******************************/ mitk::PlanarPolygon::Pointer figure2 = mitk::PlanarPolygon::New(); figure2->SetGeometry2D( geom ); mitk::Point2D pnt1; pnt1[0] = 9.5; pnt1[1] = 0.5; figure2->PlaceFigure( pnt1 ); mitk::Point2D pnt2; pnt2[0] = 9.5; pnt2[1] = 2.5; figure2->SetControlPoint( 1, pnt2, true ); mitk::Point2D pnt3; pnt3[0] = 11.5; pnt3[1] = 2.5; figure2->SetControlPoint( 2, pnt3, true ); figure2->GetPolyLine(0); testCase test; test.id = testCases.size(); test.figure = figure2; test.mean = 212.66; test.sd = 73.32; testCases.push_back( test ); } return testCases; } - static void GetTestConvolutionMask(itk::Image< float, 2 > *mask, - double RadiusInMM) - { - /***************************** - * radius of 2.5 mm - * -> expected ConvolutionMask: - * 0,25 0,75 1 0,75 0,25 - * 0,75 1 1 1 0,75 - * 1 1 1 1 1 - * 0,75 1 1 1 0,75 - * 0,25 0,75 1 0,75 0,25 - ******************************/ - - typedef itk::Image ImageType; - - ImageType::Pointer convolutionMask = ImageType::New(); - - ImageType::RegionType region; - ImageType::IndexType start; - ImageType::SizeType size; - ImageType::IndexType pixelIndex; - ImageType::SpacingType spacing; - ImageType::IndexType origin; - - - for(unsigned int i = 0; i < 2; ++i) - { - spacing[i] = 1.0; - start[i] = 0.0; - - double countIndex = 2.0 * RadiusInMM / spacing[i]; - - //Rounding up to the next integer by cast - countIndex += 0.9; - int castedIndex = static_cast(countIndex); - - //We always have an uneven number in size to determine a center-point in the convolution mask - if(castedIndex % 2 > 0 ) - { - size[i] = castedIndex; - } - else - { - size[i] = castedIndex +1; - } - } - - region.SetSize(size); - region.SetIndex(start); - - convolutionMask->SetRegions(region); - convolutionMask->SetSpacing(spacing); - convolutionMask->Allocate(); - - //y = 0 - pixelIndex[0] = 0; - pixelIndex[1] = 0; - convolutionMask->SetPixel(pixelIndex, 0.25); - pixelIndex[0] = 1; - pixelIndex[1] = 0; - convolutionMask->SetPixel(pixelIndex, 0.75); - - pixelIndex[0] = 2; - pixelIndex[1] = 0; - convolutionMask->SetPixel(pixelIndex, 1.0); - - pixelIndex[0] = 3; - pixelIndex[1] = 0; - convolutionMask->SetPixel(pixelIndex, 0.75); - - pixelIndex[0] = 4; - pixelIndex[1] = 0; - convolutionMask->SetPixel(pixelIndex, 0.25); - - //y = 1 - pixelIndex[1] = 0; - pixelIndex[1] = 1; - convolutionMask->SetPixel(pixelIndex, 0.75); - - pixelIndex[0] = 1; - pixelIndex[1] = 1; - convolutionMask->SetPixel(pixelIndex, 1.0); - - pixelIndex[0] = 2; - pixelIndex[1] = 1; - convolutionMask->SetPixel(pixelIndex, 1.0); - - pixelIndex[0] = 3; - pixelIndex[1] = 1; - convolutionMask->SetPixel(pixelIndex, 1.0); - - pixelIndex[0] = 4; - pixelIndex[1] = 1; - convolutionMask->SetPixel(pixelIndex, 0.75); - - //y = 2 - pixelIndex[1] = 0; - pixelIndex[1] = 2; - convolutionMask->SetPixel(pixelIndex, 1.0); - - pixelIndex[0] = 1; - pixelIndex[1] = 2; - convolutionMask->SetPixel(pixelIndex, 1.0); - - pixelIndex[0] = 2; - pixelIndex[1] = 2; - convolutionMask->SetPixel(pixelIndex, 1.0); - - pixelIndex[0] = 3; - pixelIndex[1] = 2; - convolutionMask->SetPixel(pixelIndex, 1.0); - - pixelIndex[0] = 4; - pixelIndex[1] = 2; - convolutionMask->SetPixel(pixelIndex, 1.0); - - //y = 3 - pixelIndex[1] = 0; - pixelIndex[1] = 3; - convolutionMask->SetPixel(pixelIndex, 0.75); - - pixelIndex[0] = 1; - pixelIndex[1] = 3; - convolutionMask->SetPixel(pixelIndex, 1.0); - - pixelIndex[0] = 2; - pixelIndex[1] = 3; - convolutionMask->SetPixel(pixelIndex, 1.0); - - pixelIndex[0] = 3; - pixelIndex[1] = 3; - convolutionMask->SetPixel(pixelIndex, 1.0); - - pixelIndex[0] = 4; - pixelIndex[1] = 3; - convolutionMask->SetPixel(pixelIndex, 0.75); - - //y = 4 - pixelIndex[1] = 0; - pixelIndex[1] = 4; - convolutionMask->SetPixel(pixelIndex, 0.25); - - pixelIndex[0] = 1; - pixelIndex[1] = 4; - convolutionMask->SetPixel(pixelIndex, 0.75); - - pixelIndex[0] = 2; - pixelIndex[1] = 4; - convolutionMask->SetPixel(pixelIndex, 1.0); - - pixelIndex[0] = 3; - pixelIndex[1] = 4; - convolutionMask->SetPixel(pixelIndex, 0.75); - - pixelIndex[0] = 4; - pixelIndex[1] = 4; - convolutionMask->SetPixel(pixelIndex, 0.25); - - typedef itk::ImageRegionConstIteratorWithIndex IteratorType; - IteratorType it(convolutionMask, convolutionMask->GetLargestPossibleRegion()); - - std::cout << std::endl << std::endl; - std::cout << "Selbsterzeugtes Bild: "<< std::endl; - - //#ifdef DEBUG_SEGMENTATION_CORRECTION - // Debug Ausgabe: - unsigned int lastZ = 1000000000; - unsigned int lastY = 1000000000; - for(it.GoToBegin(); !it.IsAtEnd(); ++it) - { - double x = it.Get(); - if (it.GetIndex()[1] != lastY) - { - std::cout << std::endl; - lastY = it.GetIndex()[1]; - } - if (it.GetIndex()[0] != lastZ) - { - std::cout << x << " "; - lastZ = it.GetIndex()[0]; - } - } - //#endif - - typedef itk::SimilarityIndexImageFilter FilterType; - FilterType::Pointer similiarityFilter = FilterType::New(); - - similiarityFilter->SetInput1(mask); - similiarityFilter->SetInput2(convolutionMask); - similiarityFilter->Update(); - - ImageType::Pointer similarityImage = similiarityFilter->GetOutput(); - IteratorType iterator(similarityImage, similarityImage->GetLargestPossibleRegion()); - - std::cout << "Bild nach Filter: "<< std::endl; - //#ifdef DEBUG_SEGMENTATION_CORRECTION - // Debug Ausgabe: - lastZ = 1000000000; - lastY = 1000000000; - for(iterator.GoToBegin(); !iterator.IsAtEnd(); ++iterator) - { - double x = iterator.Get(); - if (iterator.GetIndex()[1] != lastY) - { - std::cout << std::endl; - lastY = iterator.GetIndex()[1]; - } - if (iterator.GetIndex()[0] != lastZ) - { - std::cout << x << " "; - lastZ = iterator.GetIndex()[0]; - } - } - //#endif - - } // loads the test image - static mitk::Image::Pointer GetTestImage() + static mitk::Image::Pointer GetTestImage(const char *nameOfFile, const char *path) { mitk::StandardFileLocations::Pointer locator = mitk::StandardFileLocations::GetInstance(); - const std::string filename = locator->FindFile("testimage.dcm", "Modules/MitkExt/Testing/Data"); + const std::string filename = locator->FindFile(nameOfFile, path); if (filename.empty()) { MITK_ERROR << "Could not find test file"; return NULL; } else { - MITK_INFO << "Found testimage.dcm"; + MITK_INFO << "Found " << nameOfFile; } itk::FilenamesContainer file; file.push_back( filename ); mitk::DicomSeriesReader* reader = new mitk::DicomSeriesReader; mitk::DataNode::Pointer node = reader->LoadDicomSeries( file, false, false ); mitk::Image::Pointer image = dynamic_cast( node->GetData() ); return image; } }; int TestUnitilizedImage() { /***************************** * loading uninitialized image to datastorage ******************************/ std::cout << "Testing loading uninitialized image to datastorage:"; try{ MITK_TEST_FOR_EXCEPTION_BEGIN(mitk::Exception) mitk::Image::Pointer image = mitk::Image::New(); mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData(image); mitk::ImageStatisticsCalculator::Pointer is = mitk::ImageStatisticsCalculator::New(); is->ComputeStatistics(); MITK_TEST_FOR_EXCEPTION_END(mitk::Exception) } catch (const mitk::Exception&) { std::cout << "Success: Loading uninitialized image has thrown an exception." << std::endl; } return 0; } +// static method TestHotspotSearch(image): + // -> formuliert verschiedene TestFälle: Polygon und Masken (Segmentierung) + static void TestHotspotSearch(mitk::Image::Pointer image, unsigned int testCase) + { + // testCase for PlanarFigure-1 + if(testCase == 0) + { + + } + + // testCase for Mask-1 + if(testCase == 1) + { + + } + } + + static void CreateTestImage() + { + + typedef itk::Image ThreeDimensionalImageType; + ThreeDimensionalImageType::Pointer image = ThreeDimensionalImageType::New(); + + ThreeDimensionalImageType::SizeType size; + ThreeDimensionalImageType::SpacingType spacing; + ThreeDimensionalImageType::IndexType start; + + for(int i = 0; i < 3; ++i) + { + size[i] = 10; + spacing[i] = 1.00; + start[i] = 0.00; + } + + ThreeDimensionalImageType::RegionType region; + region.SetIndex(start); + region.SetSize(size); + + image->SetRegions(region); + image->Allocate(); + + typedef itk::ImageRegionIteratorWithIndex IteratorType; + IteratorType imageIt(image, region); + + + for(imageIt.GoToBegin(); !imageIt.IsAtEnd(); ++imageIt) + { + imageIt.Set(10); + } + + for(unsigned int x = 2; x <= 8; ++x) + { + for(unsigned int y = 2; y <= 8; ++y) + { + for(unsigned int z = 2; z <= 8; ++z) + { + ThreeDimensionalImageType::IndexType pixelIndex; + pixelIndex[0] = x; + pixelIndex[1] = y; + pixelIndex[2] = z; + + image->SetPixel(pixelIndex, 200); + } + } + } + + ThreeDimensionalImageType::IndexType pixelIndex; + pixelIndex[0] = 5; + pixelIndex[1] = 5; + pixelIndex[2] = 5; + + image->SetPixel(pixelIndex, 255); + + typedef itk::Image ThreeDimensionalMaskImageType; + ThreeDimensionalMaskImageType::Pointer mask = ThreeDimensionalMaskImageType::New(); + + ThreeDimensionalMaskImageType::SizeType maskSize; + ThreeDimensionalMaskImageType::SpacingType maskSpacing; + ThreeDimensionalMaskImageType::IndexType maskStart; + + for(int i = 0; i < 3; ++i) + { + maskSize[i] = 10; + maskSpacing[i] = 1.00; + maskStart[i] = 0.00; + } + + ThreeDimensionalMaskImageType::RegionType maskRegion; + maskRegion.SetIndex(maskStart); + maskRegion.SetSize(maskSize); + + mask->SetRegions(maskRegion); + mask->Allocate(); + + typedef itk::ImageRegionIteratorWithIndex MaskIteratorType; + MaskIteratorType maskIt(mask, maskRegion); + + + for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt) + { + maskIt.Set(1); + } + + /*typedef itk::ImageFileWriter WriterType; + WriterType::Pointer writer = WriterType::New(); + writer->SetInput(ThreeDimage); + writer->SetFileName("D:\\dreiD.dcm"); + writer->Update();*/ + } + int mitkImageStatisticsCalculatorTest(int, char* []) { // always start with this! MITK_TEST_BEGIN("mitkImageStatisticsCalculatorTest") //QCoreApplication app(argc, argv); - mitk::Image::Pointer image = mitkImageStatisticsCalculatorTestClass::GetTestImage(); + mitk::Image::Pointer image = mitkImageStatisticsCalculatorTestClass::GetTestImage(' ', ' '); MITK_TEST_CONDITION_REQUIRED( image.IsNotNull(), "Loading test image" ); mitk::Geometry2D::Pointer geom = image->GetSlicedGeometry()->GetGeometry2D(0); std::vector allTestCases = mitkImageStatisticsCalculatorTestClass::InitializeTestCases( geom ); for ( std::vector::size_type i=0; i formuliert verschiedene TestFälle: Polygon und Masken (Segmentierung) - TestUnitilizedImage(); MITK_TEST_END() } diff --git a/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp b/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp index bc1bf783c8..7a058d23b9 100644 --- a/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp +++ b/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp @@ -1,1942 +1,2056 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkImageStatisticsCalculator.h" #include "mitkImageAccessByItk.h" #include "mitkImageCast.h" #include "mitkExtractImageFilter.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include + //#define DEBUG_HOTSPOTSEARCH #define _USE_MATH_DEFINES #include #if ( ( VTK_MAJOR_VERSION <= 5 ) && ( VTK_MINOR_VERSION<=8) ) #include "mitkvtkLassoStencilSource.h" #else #include "vtkLassoStencilSource.h" #endif #include #include namespace mitk { ImageStatisticsCalculator::ImageStatisticsCalculator() : m_MaskingMode( MASKING_MODE_NONE ), m_MaskingModeChanged( false ), m_IgnorePixelValue(0.0), m_DoIgnorePixelValue(false), m_IgnorePixelValueChanged(false), m_PlanarFigureAxis (0), m_PlanarFigureSlice (0), m_PlanarFigureCoordinate0 (0), m_PlanarFigureCoordinate1 (0) { m_EmptyHistogram = HistogramType::New(); m_EmptyHistogram->SetMeasurementVectorSize(1); HistogramType::SizeType histogramSize(1); histogramSize.Fill( 256 ); m_EmptyHistogram->Initialize( histogramSize ); m_EmptyStatistics.Reset(); } ImageStatisticsCalculator::~ImageStatisticsCalculator() { } void ImageStatisticsCalculator::SetImage( const mitk::Image *image ) { if ( m_Image != image ) { m_Image = image; this->Modified(); unsigned int numberOfTimeSteps = image->GetTimeSteps(); // Initialize vectors to time-size of this image m_ImageHistogramVector.resize( numberOfTimeSteps ); m_MaskedImageHistogramVector.resize( numberOfTimeSteps ); m_PlanarFigureHistogramVector.resize( numberOfTimeSteps ); m_ImageStatisticsVector.resize( numberOfTimeSteps ); m_MaskedImageStatisticsVector.resize( numberOfTimeSteps ); m_PlanarFigureStatisticsVector.resize( numberOfTimeSteps ); m_ImageStatisticsTimeStampVector.resize( numberOfTimeSteps ); m_MaskedImageStatisticsTimeStampVector.resize( numberOfTimeSteps ); m_PlanarFigureStatisticsTimeStampVector.resize( numberOfTimeSteps ); m_ImageStatisticsCalculationTriggerVector.resize( numberOfTimeSteps ); m_MaskedImageStatisticsCalculationTriggerVector.resize( numberOfTimeSteps ); m_PlanarFigureStatisticsCalculationTriggerVector.resize( numberOfTimeSteps ); for ( unsigned int t = 0; t < image->GetTimeSteps(); ++t ) { m_ImageStatisticsTimeStampVector[t].Modified(); m_ImageStatisticsCalculationTriggerVector[t] = true; } } } void ImageStatisticsCalculator::SetImageMask( const mitk::Image *imageMask ) { if ( m_Image.IsNull() ) { itkExceptionMacro( << "Image needs to be set first!" ); } if ( m_Image->GetTimeSteps() != imageMask->GetTimeSteps() ) { itkExceptionMacro( << "Image and image mask need to have equal number of time steps!" ); } if ( m_ImageMask != imageMask ) { m_ImageMask = imageMask; this->Modified(); for ( unsigned int t = 0; t < m_Image->GetTimeSteps(); ++t ) { m_MaskedImageStatisticsTimeStampVector[t].Modified(); m_MaskedImageStatisticsCalculationTriggerVector[t] = true; } } } void ImageStatisticsCalculator::SetPlanarFigure( mitk::PlanarFigure *planarFigure ) { if ( m_Image.IsNull() ) { itkExceptionMacro( << "Image needs to be set first!" ); } if ( m_PlanarFigure != planarFigure ) { m_PlanarFigure = planarFigure; this->Modified(); for ( unsigned int t = 0; t < m_Image->GetTimeSteps(); ++t ) { m_PlanarFigureStatisticsTimeStampVector[t].Modified(); m_PlanarFigureStatisticsCalculationTriggerVector[t] = true; } } } void ImageStatisticsCalculator::SetMaskingMode( unsigned int mode ) { if ( m_MaskingMode != mode ) { m_MaskingMode = mode; m_MaskingModeChanged = true; this->Modified(); } } void ImageStatisticsCalculator::SetMaskingModeToNone() { if ( m_MaskingMode != MASKING_MODE_NONE ) { m_MaskingMode = MASKING_MODE_NONE; m_MaskingModeChanged = true; this->Modified(); } } void ImageStatisticsCalculator::SetMaskingModeToImage() { if ( m_MaskingMode != MASKING_MODE_IMAGE ) { m_MaskingMode = MASKING_MODE_IMAGE; m_MaskingModeChanged = true; this->Modified(); } } void ImageStatisticsCalculator::SetMaskingModeToPlanarFigure() { if ( m_MaskingMode != MASKING_MODE_PLANARFIGURE ) { m_MaskingMode = MASKING_MODE_PLANARFIGURE; m_MaskingModeChanged = true; this->Modified(); } } void ImageStatisticsCalculator::SetIgnorePixelValue(double value) { if ( m_IgnorePixelValue != value ) { m_IgnorePixelValue = value; if(m_DoIgnorePixelValue) { m_IgnorePixelValueChanged = true; } this->Modified(); } } double ImageStatisticsCalculator::GetIgnorePixelValue() { return m_IgnorePixelValue; } void ImageStatisticsCalculator::SetDoIgnorePixelValue(bool value) { if ( m_DoIgnorePixelValue != value ) { m_DoIgnorePixelValue = value; m_IgnorePixelValueChanged = true; this->Modified(); } } bool ImageStatisticsCalculator::GetDoIgnorePixelValue() { return m_DoIgnorePixelValue; } void ImageStatisticsCalculator::SetHotspotSize(double value) { m_HotspotSize = value; } double ImageStatisticsCalculator::GetHotspotSize() { return m_HotspotSize; } void ImageStatisticsCalculator::SetCalculateHotspot(bool value) { m_CalculateHotspot = value; } bool ImageStatisticsCalculator::IsHotspotCalculated() { return m_CalculateHotspot; } bool ImageStatisticsCalculator::ComputeStatistics( unsigned int timeStep ) { if (m_Image.IsNull() ) { mitkThrow() << "Image not set!"; } if (!m_Image->IsInitialized()) { mitkThrow() << "Image not initialized!"; } if ( m_Image->GetReferenceCount() == 1 ) { // Image no longer valid; we are the only ones to still hold a reference on it return false; } if ( timeStep >= m_Image->GetTimeSteps() ) { throw std::runtime_error( "Error: invalid time step!" ); } // If a mask was set but we are the only ones to still hold a reference on // it, delete it. if ( m_ImageMask.IsNotNull() && (m_ImageMask->GetReferenceCount() == 1) ) { m_ImageMask = NULL; } // Check if statistics is already up-to-date unsigned long imageMTime = m_ImageStatisticsTimeStampVector[timeStep].GetMTime(); unsigned long maskedImageMTime = m_MaskedImageStatisticsTimeStampVector[timeStep].GetMTime(); unsigned long planarFigureMTime = m_PlanarFigureStatisticsTimeStampVector[timeStep].GetMTime(); bool imageStatisticsCalculationTrigger = m_ImageStatisticsCalculationTriggerVector[timeStep]; bool maskedImageStatisticsCalculationTrigger = m_MaskedImageStatisticsCalculationTriggerVector[timeStep]; bool planarFigureStatisticsCalculationTrigger = m_PlanarFigureStatisticsCalculationTriggerVector[timeStep]; if ( !m_IgnorePixelValueChanged && ((m_MaskingMode != MASKING_MODE_NONE) || (imageMTime > m_Image->GetMTime() && !imageStatisticsCalculationTrigger)) && ((m_MaskingMode != MASKING_MODE_IMAGE) || (maskedImageMTime > m_ImageMask->GetMTime() && !maskedImageStatisticsCalculationTrigger)) && ((m_MaskingMode != MASKING_MODE_PLANARFIGURE) || (planarFigureMTime > m_PlanarFigure->GetMTime() && !planarFigureStatisticsCalculationTrigger)) ) { // Statistics is up to date! if ( m_MaskingModeChanged ) { m_MaskingModeChanged = false; return true; } else { return false; } } // Reset state changed flag m_MaskingModeChanged = false; m_IgnorePixelValueChanged = false; // Depending on masking mode, extract and/or generate the required image // and mask data from the user input this->ExtractImageAndMask( timeStep ); StatisticsContainer *statisticsContainer; HistogramContainer *histogramContainer; switch ( m_MaskingMode ) { case MASKING_MODE_NONE: default: if(!m_DoIgnorePixelValue) { statisticsContainer = &m_ImageStatisticsVector[timeStep]; histogramContainer = &m_ImageHistogramVector[timeStep]; m_ImageStatisticsTimeStampVector[timeStep].Modified(); m_ImageStatisticsCalculationTriggerVector[timeStep] = false; } else { statisticsContainer = &m_MaskedImageStatisticsVector[timeStep]; histogramContainer = &m_MaskedImageHistogramVector[timeStep]; m_MaskedImageStatisticsTimeStampVector[timeStep].Modified(); m_MaskedImageStatisticsCalculationTriggerVector[timeStep] = false; } break; case MASKING_MODE_IMAGE: statisticsContainer = &m_MaskedImageStatisticsVector[timeStep]; histogramContainer = &m_MaskedImageHistogramVector[timeStep]; m_MaskedImageStatisticsTimeStampVector[timeStep].Modified(); m_MaskedImageStatisticsCalculationTriggerVector[timeStep] = false; break; case MASKING_MODE_PLANARFIGURE: statisticsContainer = &m_PlanarFigureStatisticsVector[timeStep]; histogramContainer = &m_PlanarFigureHistogramVector[timeStep]; m_PlanarFigureStatisticsTimeStampVector[timeStep].Modified(); m_PlanarFigureStatisticsCalculationTriggerVector[timeStep] = false; break; } // Calculate statistics and histogram(s) if ( m_InternalImage->GetDimension() == 3 ) { if ( m_MaskingMode == MASKING_MODE_NONE && !m_DoIgnorePixelValue ) { AccessFixedDimensionByItk_2( m_InternalImage, InternalCalculateStatisticsUnmasked, 3, statisticsContainer, histogramContainer ); } else { AccessFixedDimensionByItk_3( m_InternalImage, InternalCalculateStatisticsMasked, 3, m_InternalImageMask3D.GetPointer(), statisticsContainer, histogramContainer ); } } else if ( m_InternalImage->GetDimension() == 2 ) { if ( m_MaskingMode == MASKING_MODE_NONE && !m_DoIgnorePixelValue ) { AccessFixedDimensionByItk_2( m_InternalImage, InternalCalculateStatisticsUnmasked, 2, statisticsContainer, histogramContainer ); } else { AccessFixedDimensionByItk_3( m_InternalImage, InternalCalculateStatisticsMasked, 2, m_InternalImageMask2D.GetPointer(), statisticsContainer, histogramContainer ); } } else { MITK_ERROR << "ImageStatistics: Image dimension not supported!"; } // Release unused image smart pointers to free memory m_InternalImage = mitk::Image::ConstPointer(); m_InternalImageMask3D = MaskImage3DType::Pointer(); m_InternalImageMask2D = MaskImage2DType::Pointer(); return true; } const ImageStatisticsCalculator::HistogramType * ImageStatisticsCalculator::GetHistogram( unsigned int timeStep, unsigned int label ) const { if ( m_Image.IsNull() || (timeStep >= m_Image->GetTimeSteps()) ) { return NULL; } switch ( m_MaskingMode ) { case MASKING_MODE_NONE: default: { if(m_DoIgnorePixelValue) return m_MaskedImageHistogramVector[timeStep][label]; return m_ImageHistogramVector[timeStep][label]; } case MASKING_MODE_IMAGE: return m_MaskedImageHistogramVector[timeStep][label]; case MASKING_MODE_PLANARFIGURE: return m_PlanarFigureHistogramVector[timeStep][label]; } } const ImageStatisticsCalculator::HistogramContainer & ImageStatisticsCalculator::GetHistogramVector( unsigned int timeStep ) const { if ( m_Image.IsNull() || (timeStep >= m_Image->GetTimeSteps()) ) { return m_EmptyHistogramContainer; } switch ( m_MaskingMode ) { case MASKING_MODE_NONE: default: { if(m_DoIgnorePixelValue) return m_MaskedImageHistogramVector[timeStep]; return m_ImageHistogramVector[timeStep]; } case MASKING_MODE_IMAGE: return m_MaskedImageHistogramVector[timeStep]; case MASKING_MODE_PLANARFIGURE: return m_PlanarFigureHistogramVector[timeStep]; } } const ImageStatisticsCalculator::Statistics & ImageStatisticsCalculator::GetStatistics( unsigned int timeStep, unsigned int label ) const { if ( m_Image.IsNull() || (timeStep >= m_Image->GetTimeSteps()) ) { return m_EmptyStatistics; } switch ( m_MaskingMode ) { case MASKING_MODE_NONE: default: { if(m_DoIgnorePixelValue) return m_MaskedImageStatisticsVector[timeStep][label]; return m_ImageStatisticsVector[timeStep][label]; } case MASKING_MODE_IMAGE: return m_MaskedImageStatisticsVector[timeStep][label]; case MASKING_MODE_PLANARFIGURE: return m_PlanarFigureStatisticsVector[timeStep][label]; } } const ImageStatisticsCalculator::Statistics & ImageStatisticsCalculator::GetHotspotStatistics( unsigned int timeStep, unsigned int label) const { if(m_Image.IsNull() || (timeStep >= m_Image->GetTimeSteps()) ) return m_EmptyStatistics; return m_MaskedImageHotspotStatisticsVector[timeStep][label]; } const ImageStatisticsCalculator::StatisticsContainer & ImageStatisticsCalculator::GetStatisticsVector( unsigned int timeStep ) const { if ( m_Image.IsNull() || (timeStep >= m_Image->GetTimeSteps()) ) { return m_EmptyStatisticsContainer; } switch ( m_MaskingMode ) { case MASKING_MODE_NONE: default: { if(m_DoIgnorePixelValue) return m_MaskedImageStatisticsVector[timeStep]; return m_ImageStatisticsVector[timeStep]; } case MASKING_MODE_IMAGE: return m_MaskedImageStatisticsVector[timeStep]; case MASKING_MODE_PLANARFIGURE: return m_PlanarFigureStatisticsVector[timeStep]; } } void ImageStatisticsCalculator::ExtractImageAndMask( unsigned int timeStep ) { if ( m_Image.IsNull() ) { throw std::runtime_error( "Error: image empty!" ); } if ( timeStep >= m_Image->GetTimeSteps() ) { throw std::runtime_error( "Error: invalid time step!" ); } ImageTimeSelector::Pointer imageTimeSelector = ImageTimeSelector::New(); imageTimeSelector->SetInput( m_Image ); imageTimeSelector->SetTimeNr( timeStep ); imageTimeSelector->UpdateLargestPossibleRegion(); mitk::Image *timeSliceImage = imageTimeSelector->GetOutput(); switch ( m_MaskingMode ) { case MASKING_MODE_NONE: { m_InternalImage = timeSliceImage; m_InternalImageMask2D = NULL; m_InternalImageMask3D = NULL; if(m_DoIgnorePixelValue) { if( m_InternalImage->GetDimension() == 3 ) { CastToItkImage( timeSliceImage, m_InternalImageMask3D ); m_InternalImageMask3D->FillBuffer(1); } if( m_InternalImage->GetDimension() == 2 ) { CastToItkImage( timeSliceImage, m_InternalImageMask2D ); m_InternalImageMask2D->FillBuffer(1); } } break; } case MASKING_MODE_IMAGE: { if ( m_ImageMask.IsNotNull() && (m_ImageMask->GetReferenceCount() > 1) ) { if ( timeStep < m_ImageMask->GetTimeSteps() ) { ImageTimeSelector::Pointer maskedImageTimeSelector = ImageTimeSelector::New(); maskedImageTimeSelector->SetInput( m_ImageMask ); maskedImageTimeSelector->SetTimeNr( timeStep ); maskedImageTimeSelector->UpdateLargestPossibleRegion(); mitk::Image *timeSliceMaskedImage = maskedImageTimeSelector->GetOutput(); m_InternalImage = timeSliceImage; CastToItkImage( timeSliceMaskedImage, m_InternalImageMask3D ); } else { throw std::runtime_error( "Error: image mask has not enough time steps!" ); } } else { throw std::runtime_error( "Error: image mask empty!" ); } break; } case MASKING_MODE_PLANARFIGURE: { m_InternalImageMask2D = NULL; if ( m_PlanarFigure.IsNull() ) { throw std::runtime_error( "Error: planar figure empty!" ); } if ( !m_PlanarFigure->IsClosed() ) { throw std::runtime_error( "Masking not possible for non-closed figures" ); } const Geometry3D *imageGeometry = timeSliceImage->GetGeometry(); if ( imageGeometry == NULL ) { throw std::runtime_error( "Image geometry invalid!" ); } const Geometry2D *planarFigureGeometry2D = m_PlanarFigure->GetGeometry2D(); if ( planarFigureGeometry2D == NULL ) { throw std::runtime_error( "Planar-Figure not yet initialized!" ); } const PlaneGeometry *planarFigureGeometry = dynamic_cast< const PlaneGeometry * >( planarFigureGeometry2D ); if ( planarFigureGeometry == NULL ) { throw std::runtime_error( "Non-planar planar figures not supported!" ); } // Find principal direction of PlanarFigure in input image unsigned int axis; if ( !this->GetPrincipalAxis( imageGeometry, planarFigureGeometry->GetNormal(), axis ) ) { throw std::runtime_error( "Non-aligned planar figures not supported!" ); } m_PlanarFigureAxis = axis; // Find slice number corresponding to PlanarFigure in input image MaskImage3DType::IndexType index; imageGeometry->WorldToIndex( planarFigureGeometry->GetOrigin(), index ); unsigned int slice = index[axis]; m_PlanarFigureSlice = slice; // Extract slice with given position and direction from image ExtractImageFilter::Pointer imageExtractor = ExtractImageFilter::New(); imageExtractor->SetInput( timeSliceImage ); imageExtractor->SetSliceDimension( axis ); imageExtractor->SetSliceIndex( slice ); imageExtractor->Update(); m_InternalImage = imageExtractor->GetOutput(); // Compute mask from PlanarFigure AccessFixedDimensionByItk_1( m_InternalImage, InternalCalculateMaskFromPlanarFigure, 2, axis ); } } if(m_DoIgnorePixelValue) { if ( m_InternalImage->GetDimension() == 3 ) { AccessFixedDimensionByItk_1( m_InternalImage, InternalMaskIgnoredPixels, 3, m_InternalImageMask3D.GetPointer() ); } else if ( m_InternalImage->GetDimension() == 2 ) { AccessFixedDimensionByItk_1( m_InternalImage, InternalMaskIgnoredPixels, 2, m_InternalImageMask2D.GetPointer() ); } } } bool ImageStatisticsCalculator::GetPrincipalAxis( const Geometry3D *geometry, Vector3D vector, unsigned int &axis ) { vector.Normalize(); for ( unsigned int i = 0; i < 3; ++i ) { Vector3D axisVector = geometry->GetAxisVector( i ); axisVector.Normalize(); if ( fabs( fabs( axisVector * vector ) - 1.0) < mitk::eps ) { axis = i; return true; } } return false; } template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator::InternalCalculateStatisticsUnmasked( const itk::Image< TPixel, VImageDimension > *image, StatisticsContainer *statisticsContainer, HistogramContainer* histogramContainer ) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< unsigned short, VImageDimension > MaskImageType; typedef typename ImageType::IndexType IndexType; typedef itk::Statistics::ScalarImageToHistogramGenerator< ImageType > HistogramGeneratorType; statisticsContainer->clear(); histogramContainer->clear(); // Progress listening... typedef itk::SimpleMemberCommand< ImageStatisticsCalculator > ITKCommandType; ITKCommandType::Pointer progressListener; progressListener = ITKCommandType::New(); progressListener->SetCallbackFunction( this, &ImageStatisticsCalculator::UnmaskedStatisticsProgressUpdate ); // Issue 100 artificial progress events since ScalarIMageToHistogramGenerator // does not (yet?) support progress reporting this->InvokeEvent( itk::StartEvent() ); for ( unsigned int i = 0; i < 100; ++i ) { this->UnmaskedStatisticsProgressUpdate(); } // Calculate statistics (separate filter) typedef itk::StatisticsImageFilter< ImageType > StatisticsFilterType; typename StatisticsFilterType::Pointer statisticsFilter = StatisticsFilterType::New(); statisticsFilter->SetInput( image ); unsigned long observerTag = statisticsFilter->AddObserver( itk::ProgressEvent(), progressListener ); statisticsFilter->Update(); statisticsFilter->RemoveObserver( observerTag ); this->InvokeEvent( itk::EndEvent() ); // Calculate minimum and maximum typedef itk::MinimumMaximumImageCalculator< ImageType > MinMaxFilterType; typename MinMaxFilterType::Pointer minMaxFilter = MinMaxFilterType::New(); minMaxFilter->SetImage( image ); unsigned long observerTag2 = minMaxFilter->AddObserver( itk::ProgressEvent(), progressListener ); minMaxFilter->Compute(); minMaxFilter->RemoveObserver( observerTag2 ); this->InvokeEvent( itk::EndEvent() ); Statistics statistics; statistics.Reset(); statistics.Label = 1; statistics.N = image->GetBufferedRegion().GetNumberOfPixels(); statistics.Min = statisticsFilter->GetMinimum(); statistics.Max = statisticsFilter->GetMaximum(); statistics.Mean = statisticsFilter->GetMean(); statistics.Median = 0.0; statistics.Sigma = statisticsFilter->GetSigma(); statistics.RMS = sqrt( statistics.Mean * statistics.Mean + statistics.Sigma * statistics.Sigma ); statistics.MinIndex.set_size(image->GetImageDimension()); statistics.MaxIndex.set_size(image->GetImageDimension()); for (int i=0; iGetIndexOfMaximum()[i]; statistics.MinIndex[i] = minMaxFilter->GetIndexOfMinimum()[i]; } statisticsContainer->push_back( statistics ); // Calculate histogram typename HistogramGeneratorType::Pointer histogramGenerator = HistogramGeneratorType::New(); histogramGenerator->SetInput( image ); histogramGenerator->SetMarginalScale( 100 ); histogramGenerator->SetNumberOfBins( 768 ); histogramGenerator->SetHistogramMin( statistics.Min ); histogramGenerator->SetHistogramMax( statistics.Max ); histogramGenerator->Compute(); histogramContainer->push_back( histogramGenerator->GetOutput() ); } template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator::InternalMaskIgnoredPixels( const itk::Image< TPixel, VImageDimension > *image, itk::Image< unsigned short, VImageDimension > *maskImage ) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< unsigned short, VImageDimension > MaskImageType; itk::ImageRegionIterator itmask(maskImage, maskImage->GetLargestPossibleRegion()); itk::ImageRegionConstIterator itimage(image, image->GetLargestPossibleRegion()); itmask.GoToBegin(); itimage.GoToBegin(); while( !itmask.IsAtEnd() ) { if(m_IgnorePixelValue == itimage.Get()) { itmask.Set(0); } ++itmask; ++itimage; } } template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator::InternalCalculateStatisticsMasked( const itk::Image< TPixel, VImageDimension > *image, itk::Image< unsigned short, VImageDimension > *maskImage, StatisticsContainer* statisticsContainer, HistogramContainer* histogramContainer ) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< unsigned short, VImageDimension > MaskImageType; typedef typename ImageType::IndexType IndexType; typedef typename ImageType::PointType PointType; typedef typename ImageType::SpacingType SpacingType; typedef itk::LabelStatisticsImageFilter< ImageType, MaskImageType > LabelStatisticsFilterType; typedef itk::ChangeInformationImageFilter< MaskImageType > ChangeInformationFilterType; typedef itk::ExtractImageFilter< ImageType, ImageType > ExtractImageFilterType; statisticsContainer->clear(); histogramContainer->clear(); // Make sure that mask is set if ( maskImage == NULL ) { itkExceptionMacro( << "Mask image needs to be set!" ); } // Make sure that spacing of mask and image are the same SpacingType imageSpacing = image->GetSpacing(); SpacingType maskSpacing = maskImage->GetSpacing(); PointType zeroPoint; zeroPoint.Fill( 0.0 ); if ( (zeroPoint + imageSpacing).SquaredEuclideanDistanceTo( (zeroPoint + maskSpacing) ) > mitk::eps ) { itkExceptionMacro( << "Mask needs to have same spacing as image! (Image spacing: " << imageSpacing << "; Mask spacing: " << maskSpacing << ")" ); } // Make sure that orientation of mask and image are the same typedef typename ImageType::DirectionType DirectionType; DirectionType imageDirection = image->GetDirection(); DirectionType maskDirection = maskImage->GetDirection(); for( int i = 0; i < imageDirection.ColumnDimensions; ++i ) { for( int j = 0; j < imageDirection.ColumnDimensions; ++j ) { double differenceDirection = imageDirection[i][j] - maskDirection[i][j]; if ( fabs( differenceDirection ) > mitk::eps ) { itkExceptionMacro( << "Mask needs to have same direction as image! (Image direction: " << imageDirection << "; Mask direction: " << maskDirection << ")" ); } } } // Make sure that the voxels of mask and image are correctly "aligned", i.e., voxel boundaries are the same in both images PointType imageOrigin = image->GetOrigin(); PointType maskOrigin = maskImage->GetOrigin(); long offset[ImageType::ImageDimension]; typedef itk::ContinuousIndex ContinousIndexType; ContinousIndexType maskOriginContinousIndex, imageOriginContinousIndex; image->TransformPhysicalPointToContinuousIndex(maskOrigin, maskOriginContinousIndex); image->TransformPhysicalPointToContinuousIndex(imageOrigin, imageOriginContinousIndex); for ( unsigned int i = 0; i < ImageType::ImageDimension; ++i ) { double misalignment = maskOriginContinousIndex[i] - floor( maskOriginContinousIndex[i] + 0.5 ); if ( fabs( misalignment ) > mitk::eps ) { itkExceptionMacro( << "Pixels/voxels of mask and image are not sufficiently aligned! (Misalignment: " << misalignment << ")" ); } double indexCoordDistance = maskOriginContinousIndex[i] - imageOriginContinousIndex[i]; offset[i] = (int) indexCoordDistance + image->GetBufferedRegion().GetIndex()[i]; } // Adapt the origin and region (index/size) of the mask so that the origin of both are the same typename ChangeInformationFilterType::Pointer adaptMaskFilter; adaptMaskFilter = ChangeInformationFilterType::New(); adaptMaskFilter->ChangeOriginOn(); adaptMaskFilter->ChangeRegionOn(); adaptMaskFilter->SetInput( maskImage ); adaptMaskFilter->SetOutputOrigin( image->GetOrigin() ); adaptMaskFilter->SetOutputOffset( offset ); adaptMaskFilter->Update(); typename MaskImageType::Pointer adaptedMaskImage = adaptMaskFilter->GetOutput(); // Make sure that mask region is contained within image region if ( !image->GetLargestPossibleRegion().IsInside( adaptedMaskImage->GetLargestPossibleRegion() ) ) { itkExceptionMacro( << "Mask region needs to be inside of image region! (Image region: " << image->GetLargestPossibleRegion() << "; Mask region: " << adaptedMaskImage->GetLargestPossibleRegion() << ")" ); } // If mask region is smaller than image region, extract the sub-sampled region from the original image typename ImageType::SizeType imageSize = image->GetBufferedRegion().GetSize(); typename ImageType::SizeType maskSize = maskImage->GetBufferedRegion().GetSize(); bool maskSmallerImage = false; for ( unsigned int i = 0; i < ImageType::ImageDimension; ++i ) { if ( maskSize[i] < imageSize[i] ) { maskSmallerImage = true; } } typename ImageType::ConstPointer adaptedImage; if ( maskSmallerImage ) { typename ExtractImageFilterType::Pointer extractImageFilter = ExtractImageFilterType::New(); extractImageFilter->SetInput( image ); extractImageFilter->SetExtractionRegion( adaptedMaskImage->GetBufferedRegion() ); extractImageFilter->Update(); adaptedImage = extractImageFilter->GetOutput(); } else { adaptedImage = image; } // Initialize Filter typedef itk::StatisticsImageFilter< ImageType > StatisticsFilterType; typename StatisticsFilterType::Pointer statisticsFilter = StatisticsFilterType::New(); statisticsFilter->SetInput( adaptedImage ); statisticsFilter->Update(); int numberOfBins = ( m_DoIgnorePixelValue && (m_MaskingMode == MASKING_MODE_NONE) ) ? 768 : 384; typename LabelStatisticsFilterType::Pointer labelStatisticsFilter; labelStatisticsFilter = LabelStatisticsFilterType::New(); labelStatisticsFilter->SetInput( adaptedImage ); labelStatisticsFilter->SetLabelInput( adaptedMaskImage ); labelStatisticsFilter->UseHistogramsOn(); labelStatisticsFilter->SetHistogramParameters( numberOfBins, statisticsFilter->GetMinimum(), statisticsFilter->GetMaximum() ); // Add progress listening typedef itk::SimpleMemberCommand< ImageStatisticsCalculator > ITKCommandType; ITKCommandType::Pointer progressListener; progressListener = ITKCommandType::New(); progressListener->SetCallbackFunction( this, &ImageStatisticsCalculator::MaskedStatisticsProgressUpdate ); unsigned long observerTag = labelStatisticsFilter->AddObserver( itk::ProgressEvent(), progressListener ); // Execute filter this->InvokeEvent( itk::StartEvent() ); // Make sure that only the mask region is considered (otherwise, if the mask region is smaller // than the image region, the Update() would result in an exception). labelStatisticsFilter->GetOutput()->SetRequestedRegion( adaptedMaskImage->GetLargestPossibleRegion() ); // Execute the filter labelStatisticsFilter->Update(); this->InvokeEvent( itk::EndEvent() ); labelStatisticsFilter->RemoveObserver( observerTag ); // Find all relevant labels of mask (other than 0) std::list< int > relevantLabels; bool maskNonEmpty = false; unsigned int i; for ( i = 1; i < 4096; ++i ) { if ( labelStatisticsFilter->HasLabel( i ) ) { relevantLabels.push_back( i ); maskNonEmpty = true; } } if ( maskNonEmpty ) { std::list< int >::iterator it; for ( it = relevantLabels.begin(), i = 0; it != relevantLabels.end(); ++it, ++i ) { histogramContainer->push_back( HistogramType::ConstPointer( labelStatisticsFilter->GetHistogram( (*it) ) ) ); Statistics statistics; statistics.Label = (*it); statistics.N = labelStatisticsFilter->GetCount( *it ); statistics.Min = labelStatisticsFilter->GetMinimum( *it ); statistics.Max = labelStatisticsFilter->GetMaximum( *it ); statistics.Mean = labelStatisticsFilter->GetMean( *it ); statistics.Median = labelStatisticsFilter->GetMedian( *it ); statistics.Sigma = labelStatisticsFilter->GetSigma( *it ); statistics.RMS = sqrt( statistics.Mean * statistics.Mean + statistics.Sigma * statistics.Sigma ); // restrict image to mask area for min/max index calculation typedef itk::MaskImageFilter< ImageType, MaskImageType, ImageType > MaskImageFilterType; typename MaskImageFilterType::Pointer masker = MaskImageFilterType::New(); masker->SetOutsideValue( (statistics.Min+statistics.Max)/2 ); masker->SetInput1(adaptedImage); masker->SetInput2(adaptedMaskImage); masker->Update(); // get index of minimum and maximum typedef itk::MinimumMaximumImageCalculator< ImageType > MinMaxFilterType; typename MinMaxFilterType::Pointer minMaxFilter = MinMaxFilterType::New(); minMaxFilter->SetImage( masker->GetOutput() ); unsigned long observerTag2 = minMaxFilter->AddObserver( itk::ProgressEvent(), progressListener ); minMaxFilter->Compute(); minMaxFilter->RemoveObserver( observerTag2 ); this->InvokeEvent( itk::EndEvent() ); statistics.MinIndex.set_size(adaptedImage->GetImageDimension()); statistics.MaxIndex.set_size(adaptedImage->GetImageDimension()); typename MinMaxFilterType::IndexType tempMaxIndex = minMaxFilter->GetIndexOfMaximum(); typename MinMaxFilterType::IndexType tempMinIndex = minMaxFilter->GetIndexOfMinimum(); // FIX BUG 14644 //If a PlanarFigure is used for segmentation the //adaptedImage is a single slice (2D). Adding the // 3. dimension. if (m_MaskingMode == MASKING_MODE_PLANARFIGURE && m_Image->GetDimension()==3) { statistics.MaxIndex.set_size(m_Image->GetDimension()); statistics.MaxIndex[m_PlanarFigureCoordinate0]=tempMaxIndex[0]; statistics.MaxIndex[m_PlanarFigureCoordinate1]=tempMaxIndex[1]; statistics.MaxIndex[m_PlanarFigureAxis]=m_PlanarFigureSlice; statistics.MinIndex.set_size(m_Image->GetDimension()); statistics.MinIndex[m_PlanarFigureCoordinate0]=tempMinIndex[0]; statistics.MinIndex[m_PlanarFigureCoordinate1]=tempMinIndex[1]; statistics.MinIndex[m_PlanarFigureAxis]=m_PlanarFigureSlice; } else { for (int i = 0; i ThreeDimensionalImageType; - ThreeDimensionalImageType::Pointer image = ThreeDimensionalImageType::New(); + typedef itk::Image ThreeDimensionalImageType; + ThreeDimensionalImageType::Pointer image = ThreeDimensionalImageType::New(); ThreeDimensionalImageType::SizeType size; ThreeDimensionalImageType::SpacingType spacing; ThreeDimensionalImageType::IndexType start; for(int i = 0; i < 3; ++i) { - size[i] = 100.00; + size[i] = 50.00; spacing[i] = 1.00; start[i] = 0.00; } ThreeDimensionalImageType::RegionType region; region.SetIndex(start); region.SetSize(size); image->SetRegions(region); image->Allocate(); typedef itk::ImageRegionIteratorWithIndex IteratorType; IteratorType imageIt(image, region); for(imageIt.GoToBegin(); !imageIt.IsAtEnd(); ++imageIt) { imageIt.Set(15); } - for(unsigned int x = 70; x <= 80; ++x) + /*for(unsigned int x = 70; x <= 80; ++x) { for(unsigned int y = 70; y <= 80; ++y) { for(unsigned int z = 70; z <= 80; ++z) { ThreeDimensionalImageType::IndexType pixelIndex; pixelIndex[0] = x; pixelIndex[1] = y; pixelIndex[2] = z; image->SetPixel(pixelIndex, 50); } } - } + }*/ ThreeDimensionalImageType::IndexType pixelIndex; - pixelIndex[0] = 75; - pixelIndex[1] = 75; - pixelIndex[2] = 75; + pixelIndex[0] = 25; + pixelIndex[1] = 25; + pixelIndex[2] = 25; - image->SetPixel(pixelIndex, 1500); + image->SetPixel(pixelIndex, 20); pixelIndex[0] = 0; pixelIndex[1] = 0; pixelIndex[2] = 0; - image->SetPixel(pixelIndex, 10000000); + image->SetPixel(pixelIndex, 1000); /*****************************************************Creating Test Mask**********************************************/ - typedef itk::Image ThreeDimensionalMaskImageType; + typedef itk::Image ThreeDimensionalMaskImageType; ThreeDimensionalMaskImageType::Pointer mask = ThreeDimensionalMaskImageType::New(); ThreeDimensionalMaskImageType::SizeType maskSize; ThreeDimensionalMaskImageType::SpacingType maskSpacing; ThreeDimensionalMaskImageType::IndexType maskStart; for(int i = 0; i < 3; ++i) { - maskSize[i] = 100; + maskSize[i] = 50; maskSpacing[i] = 1.00; maskStart[i] = 0.00; } ThreeDimensionalMaskImageType::RegionType maskRegion; maskRegion.SetIndex(maskStart); maskRegion.SetSize(maskSize); mask->SetRegions(maskRegion); mask->Allocate(); typedef itk::ImageRegionIteratorWithIndex MaskIteratorType; MaskIteratorType maskIt(mask, maskRegion); for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt) { maskIt.Set(1); } /*****************************************************Calculate Hotspot Statistics**********************************************/ - Statistics hotspotStatistics = CalculateHotspotStatistics(image.GetPointer(), mask.GetPointer(), 3.00);; + Statistics hotspotStatistics = CalculateHotspotStatistics (adaptedImage.GetPointer(), adaptedMaskImage.GetPointer(), 3.00); statistics.HotspotMax = hotspotStatistics.HotspotMax; statistics.HotspotMin = hotspotStatistics.HotspotMin; statistics.HotspotPeak = hotspotStatistics.HotspotPeak; statistics.HotspotMaxIndex = hotspotStatistics.HotspotMaxIndex; statistics.HotspotMinIndex = hotspotStatistics.HotspotMinIndex; statistics.HotspotPeakIndex = hotspotStatistics.HotspotPeakIndex; statisticsContainer->push_back( statistics ); } } else { histogramContainer->push_back( HistogramType::ConstPointer( m_EmptyHistogram ) ); statisticsContainer->push_back( Statistics() ); } } template ImageStatisticsCalculator::MinMaxIndex ImageStatisticsCalculator::CalculateMinMaxIndex( const itk::Image *inputImage, itk::Image *maskImage, float minBoundary, float maxBoundary) { - // TODO: check if inputImage and maskImage have same dimensions typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< unsigned short, VImageDimension > MaskImageType; typedef itk::ImageRegionConstIterator MaskImageIteratorType; typedef itk::ImageRegionConstIteratorWithIndex InputImageIndexIteratorType; MaskImageIteratorType maskIt(maskImage, maskImage->GetLargestPossibleRegion()); InputImageIndexIteratorType imageIndexIt(inputImage, inputImage->GetLargestPossibleRegion()); float maxValue = itk::NumericTraits::min(); float minValue = itk::NumericTraits::max(); ImageType::IndexType maxIndex; ImageType::IndexType minIndex; + double maxBound = itk::NumericTraits::min(); + double minBound = itk::NumericTraits::max(); + + std::list< ImageType::IndexType > maxIndexList; + std::list< ImageType::IndexType > minIndexList; + + bool sameIndex = true; + for(maskIt.GoToBegin(), imageIndexIt.GoToBegin(); - !maskIt.IsAtEnd() && !imageIndexIt.IsAtEnd(); - ++maskIt, ++imageIndexIt) + !maskIt.IsAtEnd() && !imageIndexIt.IsAtEnd(); + ++maskIt, ++imageIndexIt) { if(maskIt.Get() > itk::NumericTraits::Zero) { - double indexValue = imageIndexIt.Get(); - // Calculate minimum, maximum and corresponding index-values - if(indexValue > maxValue && indexValue < maxBoundary) + double value = imageIndexIt.Get(); + ImageType::IndexType index = imageIndexIt.GetIndex(); + //Calculate minimum, maximum and corresponding index-values + if( value > maxValue && value < maxBoundary) + { + maxIndexList.clear(); + maxIndexList.push_back( index ); + maxValue = value; + } + else if ( value == maxValue ) { - maxIndex = imageIndexIt.GetIndex(); - maxValue = indexValue; + maxIndexList.push_back( index ); } - if(indexValue < minValue && indexValue > minBoundary) + if(value < minValue && value > minBoundary) + { + minIndexList.clear(); + minIndexList.push_back( index ); + minValue = value; + } + else if ( value == minValue ) { - minIndex = imageIndexIt.GetIndex(); - minValue = indexValue; + minIndexList.push_back( index ); } } } - MinMaxIndex minMax; - minMax.MinIndex.set_size(inputImage->GetImageDimension()); - minMax.MaxIndex.set_size(inputImage->GetImageDimension()); + std::list> tmpList; - for (int i = 0; i::iterator it = maxIndexList.begin(); it != maxIndexList.end(); ++it) { - minMax.MinIndex[i] = minIndex[i]; + vnl_vector indexVector; + indexVector.set_size(inputImage->GetImageDimension()); + + ImageType::IndexType index = *it; + + for(int i = 0; i < VImageDimension ; ++i) + indexVector[i] = index[i]; + + tmpList.push_back(indexVector); } - for (int i = 0; i::iterator it = minIndexList.begin(); it != minIndexList.end(); ++it) { - minMax.MaxIndex[i] = maxIndex[i]; + vnl_vector indexVector; + indexVector.set_size(inputImage->GetImageDimension()); + + ImageType::IndexType index = *it; + + for(int i = 0; i < indexVector.size(); ++i) + indexVector[i] = index[i]; + + tmpList.push_back(indexVector); } + minMax.MinIndexList = tmpList; + + minMax.MinIndex.set_size(inputImage->GetImageDimension()); + minMax.MaxIndex.set_size(inputImage->GetImageDimension()); + minMax.Max = maxValue; minMax.Min = minValue; return minMax; } template < typename TPixel, unsigned int VImageDimension> ImageStatisticsCalculator::Statistics ImageStatisticsCalculator::CalculateHotspotStatistics( - itk::Image *inputImage, + const itk::Image *inputImage, itk::Image *maskImage, double RadiusInMM) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< float, VImageDimension > MaskImageType; typedef itk::Image SphereMaskImageType; typedef itk::ContinuousIndex ContinuousIndexType; typedef itk::ImageRegionConstIteratorWithIndex IteratorType; typedef itk::ImageRegionConstIteratorWithIndex MaskIteratorType; - typedef itk::ImageRegionIteratorWithIndex SphereMaskIteratorType; + typedef itk::ImageRegionConstIteratorWithIndex SphereMaskIteratorType; typedef typename ImageType::SpacingType SpacingType; typedef typename ImageType::SizeType SizeType; typedef typename ImageType::IndexType IndexType; typedef typename ImageType::PointType PointType; typedef typename ImageType::RegionType RegionType; MaskImageType::Pointer convolutionMask = MaskImageType::New(); SpacingType spacing = inputImage->GetSpacing(); IndexType start; start.Fill(0); ContinuousIndexType origin; SizeType size; ContinuousIndexType convolutionMaskCenterCoordinate; /*****************************************************Creating convolution mask**********************************************/ for(unsigned int i = 0; i < VImageDimension; ++i) { origin[i] = 0.0; double countIndex = 2.0 * RadiusInMM / spacing[i]; // Rounding up to the next integer by cast countIndex += 0.9; int castedIndex = static_cast(countIndex); // We always have an uneven number in size to determine a center-point in the convolution mask if(castedIndex % 2 > 0 ) { size[i] = castedIndex; } else { size[i] = castedIndex +1; } convolutionMaskCenterCoordinate[i] = (size[i] -1) / 2; } RegionType region; region.SetSize(size); region.SetIndex(start); convolutionMask->SetRegions(region); convolutionMask->SetSpacing(spacing); convolutionMask->Allocate(); float subPixelDimensionX = spacing[0]/2; float subPixelDimensionY = spacing[1]/2; double distance = 0.0; double countSubPixel = 0.0; double pixelValue = 0.0; MaskIteratorType maskIt(convolutionMask,region); // Generate convolutionMask for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt) { ContinuousIndexType indexPoint(maskIt.GetIndex()); for(double x = indexPoint[0] - subPixelDimensionX / 2; x <= indexPoint[0] + subPixelDimensionX / 2; x += subPixelDimensionX) { for(double y = indexPoint[1] - subPixelDimensionY /2; y <= indexPoint[1] + subPixelDimensionY / 2; y += subPixelDimensionY) { ContinuousIndexType subPixelCenterCoordinate; subPixelCenterCoordinate[0] = x; subPixelCenterCoordinate[1] = y; for( unsigned int i = 2; i < VImageDimension; ++i ) { subPixelCenterCoordinate[i] = indexPoint[i]; } distance = subPixelCenterCoordinate.EuclideanDistanceTo(convolutionMaskCenterCoordinate); if(distance <= RadiusInMM) countSubPixel++; } } // pixelValue is the counted SubPixels divided by factor 4 pixelValue = countSubPixel / 4.00; convolutionMask->SetPixel(maskIt.GetIndex(),pixelValue); countSubPixel = 0.00; } /*****************************************************Creating Peak Image**********************************************/ typedef itk::ConvolutionImageFilter FilterType; FilterType::Pointer convolutionFilter = FilterType::New(); convolutionFilter->SetInput(inputImage); convolutionFilter->SetKernelImage(convolutionMask); convolutionFilter->SetNormalize(true); convolutionFilter->Update(); MaskImageType::Pointer peakImage = convolutionFilter->GetOutput(); + /*std::cout << std::endl << std::endl; + std::cout << "PeakImage: " << std::endl; + unsigned int lastZ = 1000000000; + unsigned int lastY = 1000000000; + + unsigned int hotspotMaskIndexCounter = 0; + MaskIteratorType hotspotMaskIt(peakImage, peakImage->GetLargestPossibleRegion() ); + for(hotspotMaskIt.GoToBegin();!hotspotMaskIt.IsAtEnd();++hotspotMaskIt) + { + + double tmp = hotspotMaskIt.Get(); + if (hotspotMaskIt.GetIndex()[1] != lastY) + { + std::cout << std::endl; + lastY = hotspotMaskIt.GetIndex()[1]; + } + if (hotspotMaskIt.GetIndex()[0] != lastZ) + { + std::cout << tmp << " "; + lastZ = hotspotMaskIt.GetIndex()[0]; + } + + hotspotMaskIndexCounter++; + + if(hotspotMaskIndexCounter > 899) { + std::cout << std::endl; + hotspotMaskIndexCounter = 0; + } + } + + std::cout << std::endl << std::endl;*/ + /*****************************************************Creating Hotspot Sphere**********************************************/ SphereMaskImageType::Pointer hotspotSphere = SphereMaskImageType::New(); SphereMaskImageType::SpacingType hotspotSphereSpacing = peakImage->GetSpacing(); SphereMaskImageType::SizeType hotspotSphereSize; SphereMaskImageType::RegionType regionSphere; SphereMaskImageType::IndexType hotspotPeakIndex; PointType hotspotOrigin; IndexType offsetInIndex; typedef itk::EllipseSpatialObject EllipseType; typedef itk::SpatialObjectToImageFilter SpatialObjectToImageFilter; const double radiusSUVHotspot = 6.2035049089940; float maxBoundary = itk::NumericTraits::max(); float minBoundary = itk::NumericTraits::min(); double hotspotPeak = itk::NumericTraits::min(); bool isInside = false; + std::list> peakIndexList; + + /* vnl_vector headIndex; + headIndex.set_size(VImageDimension); + + peakIndexList.push_back(head);*/ + + int i = 0; + while (isInside == false) { + MinMaxIndex peakInformations = CalculateMinMaxIndex(peakImage.GetPointer(), maskImage, minBoundary, maxBoundary); + peakIndexList = peakInformations.MaxIndexList; + + for(std::list>::iterator it = peakIndexList.begin(); it != peakIndexList.end(); ++it) + { + vnl_vector tmp; + tmp = *it; + std::cout << "List[" << i << "]: " << tmp << std::endl; + } if(maxBoundary == peakInformations.Max) { std::cout << "Abortion" << std::endl; break; } - SphereMaskImageType::IndexType peakIndex; - - for(unsigned int i = 0; i < VImageDimension; ++i) + if(!peakIndexList.empty()) { - double countIndex = 2.0 * radiusSUVHotspot / hotspotSphereSpacing[i]; + for(std::list>::iterator iterator = peakIndexList.begin(); iterator != peakIndexList.end(); ++iterator) + { + vnl_vector tmpIndex = *iterator; + SphereMaskImageType::IndexType peakIndex; - // Rounding up to the next integer by cast - countIndex += 0.9; - int castedIndex = static_cast(countIndex); + for(unsigned int i = 0; i < VImageDimension; ++i) + { + // Converting vnl_vector in IndexType for origin of sphere + peakIndex[i] = tmpIndex[i]; - // We always have an uneven number in size to determine a center-point in the convolution mask - if(castedIndex % 2 > 0 ) - { - hotspotSphereSize[i] = castedIndex; - } - else - { - hotspotSphereSize[i] = castedIndex +1; - } - // Calculating peakIndex for Origin - peakIndex[i] = peakInformations.MaxIndex[i]; - } + double countIndex = 2.0 * radiusSUVHotspot / hotspotSphereSpacing[i]; + + // Rounding up to the next integer by cast + countIndex += 0.9; + int castedIndex = static_cast(countIndex); + + // We always have an uneven number in size to determine a center-point in the convolution mask + if(castedIndex % 2 > 0 ) + { + hotspotSphereSize[i] = castedIndex; + } + else + { + hotspotSphereSize[i] = castedIndex +1; + } + + } - // Initialize SpatialObjectoToImageFilter - itk::SpatialObjectToImageFilter::Pointer spatialObjectToImageFilter - = SpatialObjectToImageFilter::New(); + // Initialize SpatialObjectoToImageFilter + itk::SpatialObjectToImageFilter::Pointer spatialObjectToImageFilter + = SpatialObjectToImageFilter::New(); - spatialObjectToImageFilter->SetSize(hotspotSphereSize); - spatialObjectToImageFilter->SetSpacing(hotspotSphereSpacing); + spatialObjectToImageFilter->SetSize(hotspotSphereSize); + spatialObjectToImageFilter->SetSpacing(hotspotSphereSpacing); - // Creating spatial sphere object - EllipseType::Pointer sphere = EllipseType::New(); - sphere->SetRadius(radiusSUVHotspot); + // Creating spatial sphere object + EllipseType::Pointer sphere = EllipseType::New(); + sphere->SetRadius(radiusSUVHotspot); - typedef EllipseType::TransformType TransformType; - TransformType::Pointer transform = TransformType::New(); + typedef EllipseType::TransformType TransformType; + TransformType::Pointer transform = TransformType::New(); - transform->SetIdentity(); + transform->SetIdentity(); - TransformType::OutputVectorType translation; + TransformType::OutputVectorType translation; - // Transform sphere on center-position, set pixelValues inside sphere on 1 and update - for(int i = 0; i < VImageDimension; ++i) - translation[i] = (hotspotSphereSize[i] -1) / 2; + // Transform sphere on center-position, set pixelValues inside sphere on 1 and update + for(int i = 0; i < VImageDimension; ++i) + translation[i] = (hotspotSphereSize[i] -1) / 2; - transform->Translate(translation, false); + transform->Translate(translation, false); - sphere->SetObjectToParentTransform(transform); + sphere->SetObjectToParentTransform(transform); - spatialObjectToImageFilter->SetInput(sphere); + spatialObjectToImageFilter->SetInput(sphere); - sphere->SetDefaultInsideValue(1.00); - sphere->SetDefaultOutsideValue(0.00); + sphere->SetDefaultInsideValue(1.00); + sphere->SetDefaultOutsideValue(0.00); - spatialObjectToImageFilter->SetUseObjectValue(true); - spatialObjectToImageFilter->SetOutsideValue(0); + spatialObjectToImageFilter->SetUseObjectValue(true); + spatialObjectToImageFilter->SetOutsideValue(0); - spatialObjectToImageFilter->Update(); - hotspotSphere = spatialObjectToImageFilter->GetOutput(); + spatialObjectToImageFilter->Update(); + hotspotSphere = spatialObjectToImageFilter->GetOutput(); #ifdef DEBUG_HOTSPOTSEARCH - std::cout << std::endl << std::endl; - std::cout << "hotspotMask: " << std::endl; - unsigned int lastZ = 1000000000; - unsigned int lastY = 1000000000; + std::cout << std::endl << std::endl; + std::cout << "hotspotMask: " << std::endl; + unsigned int lastZ = 1000000000; + unsigned int lastY = 1000000000; - unsigned int hotspotMaskIndexCounter = 0; - SphereMaskIteratorType hotspotMaskIt(hotspotSphere, hotspotSphere->GetLargestPossibleRegion() ); - for(hotspotMaskIt.GoToBegin();!hotspotMaskIt.IsAtEnd();++hotspotMaskIt) - { + unsigned int hotspotMaskIndexCounter = 0; + SphereMaskIteratorType hotspotMaskIt(hotspotSphere, hotspotSphere->GetLargestPossibleRegion() ); + for(hotspotMaskIt.GoToBegin();!hotspotMaskIt.IsAtEnd();++hotspotMaskIt) + { - double tmp = hotspotMaskIt.Get(); - if (hotspotMaskIt.GetIndex()[1] != lastY) - { - std::cout << std::endl; - lastY = hotspotMaskIt.GetIndex()[1]; - } - if (hotspotMaskIt.GetIndex()[0] != lastZ) - { - std::cout << tmp << " "; - lastZ = hotspotMaskIt.GetIndex()[0]; - } + double tmp = hotspotMaskIt.Get(); + if (hotspotMaskIt.GetIndex()[1] != lastY) + { + std::cout << std::endl; + lastY = hotspotMaskIt.GetIndex()[1]; + } + if (hotspotMaskIt.GetIndex()[0] != lastZ) + { + std::cout << tmp << " "; + lastZ = hotspotMaskIt.GetIndex()[0]; + } - hotspotMaskIndexCounter++; + hotspotMaskIndexCounter++; - if(hotspotMaskIndexCounter > 168) { - std::cout << std::endl; - hotspotMaskIndexCounter = 0; - } - } + if(hotspotMaskIndexCounter > 168) { + std::cout << std::endl; + hotspotMaskIndexCounter = 0; + } + } - std::cout << std::endl << std::endl; + std::cout << std::endl << std::endl; #endif - // Calculate new origin for hotspot sphere + // Calculate new origin for hotspot sphere - for(int i = 0; i < VImageDimension; ++i) - offsetInIndex[i] = (hotspotSphereSize[i] -1) / 2; + for(int i = 0; i < VImageDimension; ++i) + offsetInIndex[i] = (hotspotSphereSize[i] -1) / 2; - peakImage->TransformIndexToPhysicalPoint(peakIndex, hotspotOrigin); + peakImage->TransformIndexToPhysicalPoint(peakIndex, hotspotOrigin); - PointType offsetInPhysicalPoint; - hotspotSphere->TransformIndexToPhysicalPoint(offsetInIndex, offsetInPhysicalPoint); + PointType offsetInPhysicalPoint; + hotspotSphere->TransformIndexToPhysicalPoint(offsetInIndex, offsetInPhysicalPoint); - for(int i = 0; i < VImageDimension; ++i) - hotspotOrigin[i] -= offsetInPhysicalPoint[i]; + for(int i = 0; i < VImageDimension; ++i) + hotspotOrigin[i] -= offsetInPhysicalPoint[i]; - hotspotSphere->SetOrigin(hotspotOrigin); - hotspotSphere->Allocate(); + hotspotSphere->SetOrigin(hotspotOrigin); + hotspotSphere->Allocate(); - // If sphere is not in peakImage we need a new peakValue - bool sphereInRegion = IsSphereInsideRegion(peakImage.GetPointer(), hotspotSphere.GetPointer() ); + // If sphere is not in peakImage we need a new peakValue + bool sphereInRegion = IsSphereInsideRegion(peakImage.GetPointer(), hotspotSphere.GetPointer() ); - if(sphereInRegion == false) - { - maxBoundary = peakInformations.Max; - minBoundary = peakInformations.Min; + if(sphereInRegion == false) + { + maxBoundary = peakInformations.Max; + minBoundary = peakInformations.Min; + } + else + { + hotspotPeak = peakInformations.Max; - } - else - { - hotspotPeak = peakInformations.Max; + for(int i = 0; i < VImageDimension; ++i) + hotspotPeakIndex[i] = peakIndex[i]; - for(int i = 0; i < VImageDimension; ++i) - hotspotPeakIndex[i] = peakInformations.MaxIndex[i]; + isInside = true; + break; + } + } - isInside = true; + i++; } } + /*********************************Creating cropped inputImage for calculation of hotspot statistics****************************/ IndexType croppedStart; for(int i = 0; i < 3; ++i) croppedStart[i] = offsetInIndex[i] - hotspotSphereSize[i]; ImageType::RegionType inputRegion; peakImage->TransformPhysicalPointToIndex(hotspotOrigin,croppedStart); ImageType::RegionType::IndexType inputStart = croppedStart; ImageType::RegionType::SizeType inputSize = hotspotSphere->GetLargestPossibleRegion().GetSize(); inputRegion.SetIndex(inputStart); inputRegion.SetSize(inputSize); ImageType::RegionType outputRegion; ImageType::IndexType outputStart; outputStart.Fill(0); outputRegion.SetIndex(outputStart); outputRegion.SetSize(hotspotSphere->GetLargestPossibleRegion().GetSize()); ImageType::Pointer croppedInputImage = ImageType::New(); croppedInputImage->SetRegions(outputRegion); croppedInputImage->Allocate(); - typedef itk::ImageRegionIterator InputImageIteratorType; + typedef itk::ImageRegionConstIterator InputImageIteratorType; InputImageIteratorType inputIt(inputImage, inputRegion); - SphereMaskIteratorType croppedInputImageIt(croppedInputImage, outputRegion); - for(inputIt.GoToBegin(), croppedInputImageIt.GoToBegin(); !inputIt.IsAtEnd(); ++inputIt, ++croppedInputImageIt) - { - croppedInputImageIt.Set(inputIt.Get() ); - } + InputImageIteratorType croppedInputImageIt(croppedInputImage, outputRegion); + + for(inputIt.GoToBegin(), croppedInputImageIt.GoToBegin(); !inputIt.IsAtEnd(); ++inputIt, ++croppedInputImageIt) + { + croppedInputImage->SetPixel(croppedInputImageIt.GetIndex(), inputIt.Get()); + //croppedInputImageIt.Set(inputIt.Get() ); + } // Calculate statistics in Hotspot maxBoundary = itk::NumericTraits::max(); minBoundary = itk::NumericTraits::min(); MinMaxIndex hotspotInformations; Statistics hotspotStatistics; + std::list> maxIndexList; + hotspotInformations = CalculateMinMaxIndex(croppedInputImage.GetPointer(), hotspotSphere.GetPointer(), minBoundary, maxBoundary); + // Min and max index is the first item in the corresponding list + hotspotInformations.MaxIndex = hotspotInformations.MaxIndexList.front(); + hotspotInformations.MinIndex = hotspotInformations.MinIndexList.front(); + // Add offset to cropped indices for(int i = 0; i < VImageDimension; ++i) { hotspotInformations.MaxIndex[i] += croppedStart[i]; hotspotInformations.MinIndex[i] += croppedStart[i]; } hotspotStatistics.HotspotMin = hotspotInformations.Min; hotspotStatistics.HotspotMinIndex = hotspotInformations.MinIndex; hotspotStatistics.HotspotMax = hotspotInformations.Max; hotspotStatistics.HotspotMaxIndex = hotspotInformations.MaxIndex; hotspotStatistics.HotspotPeak = hotspotPeak; hotspotStatistics.HotspotPeakIndex.set_size(inputImage->GetImageDimension()); for (int i = 0; i< hotspotStatistics.HotspotPeakIndex.size(); ++i) { hotspotStatistics.HotspotPeakIndex[i] = hotspotPeakIndex[i]; } + return hotspotStatistics; } template < typename TPixel, unsigned int VImageDimension> bool ImageStatisticsCalculator::IsSphereInsideRegion(const itk::Image *inputImage, itk::Image *sphereImage) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef ImageType::RegionType RegionType; typedef ImageType::SizeType SizeType; typedef ImageType::PointType PointType; typedef ImageType::IndexType IndexType; RegionType region(inputImage->GetLargestPossibleRegion().GetIndex(), inputImage->GetLargestPossibleRegion().GetSize() ); typedef itk::Image SphereImageType; typedef SphereImageType::SizeType SphereSizeType; typedef SphereImageType::IndexType SphereIndexType; SphereIndexType sphereCenterIndex; SphereSizeType sphereSize = sphereImage->GetLargestPossibleRegion().GetSize(); double radius = (sphereSize[0] - 1)/2.00; bool isInside = false; for(int i = 0; i < VImageDimension; ++i) sphereCenterIndex[i] = (sphereSize[i] -1) / 2; /**********************************Calculation of sphere-coordinates which are on axis of coordinates****************************/ for(int i = 0; i < VImageDimension; ++i) { for(double j = sphereCenterIndex[i] - radius; j <= sphereCenterIndex[i]+ radius; j += 2*radius) { IndexType sphereIndex; if(i == 0) { sphereIndex[0] = j; sphereIndex[1] = sphereCenterIndex[1]; sphereIndex[2] = sphereCenterIndex[2]; } if(i == 1) { sphereIndex[0] = sphereCenterIndex[0]; sphereIndex[1] = j; sphereIndex[2] = sphereCenterIndex[2]; } if(i == 2) { sphereIndex[0] = sphereCenterIndex[0]; sphereIndex[1] = sphereCenterIndex[1]; sphereIndex[2] = j; } PointType spherePoint; sphereImage->TransformIndexToPhysicalPoint(sphereIndex, spherePoint); IndexType regionIndex; inputImage->TransformPhysicalPointToIndex(spherePoint, regionIndex); if(region.IsInside(regionIndex) == true) { isInside = true; } else { isInside = false; return isInside; } } } /**********************************Calculation of sphere-coordinates which are not on axis of coordinates************************/ IndexType sphereIndex; IndexType regionIndex; PointType spherePoint; // First run: summation, second run: subtraction bool addOrSub = false; for(int i = 0; i < 2; ++i) { if(addOrSub == false) { // Point 1 sphereIndex[0] = sphereCenterIndex[0] + radius * sin(45.00 * itk::Math::pi / 180.00) * cos(45.00 * itk::Math::pi / 180.00); sphereIndex[1] = sphereCenterIndex[1] + radius * sin(45.00 * itk::Math::pi / 180.00) * sin(45.00 * itk::Math::pi / 180.00); sphereIndex[2] = sphereCenterIndex[2] + radius * cos(45.00 * itk::Math::pi / 180.00); sphereImage->TransformIndexToPhysicalPoint(sphereIndex, spherePoint); inputImage->TransformPhysicalPointToIndex(spherePoint, regionIndex); if(region.IsInside(regionIndex) == true) { isInside = true; } else { isInside = false; return isInside; } } else { // Point 2 sphereIndex[0] = sphereCenterIndex[0] + radius * sin(45.00 * itk::Math::pi / 180.00) * cos(45.00 * itk::Math::pi / 180.00); sphereIndex[1] = -(sphereCenterIndex[1] + radius * sin(45.00 * itk::Math::pi / 180.00) * sin(45.00 * itk::Math::pi / 180.00)); sphereIndex[2] = sphereCenterIndex[2] + radius * cos(45.00 * itk::Math::pi / 180.00); sphereImage->TransformIndexToPhysicalPoint(sphereIndex, spherePoint); inputImage->TransformPhysicalPointToIndex(spherePoint, regionIndex); if(region.IsInside(regionIndex) == true) { isInside = true; } else { isInside = false; return isInside; } } // Point 3/4 sphereIndex[2] = -(sphereCenterIndex[2] + radius * cos(45.00 * itk::Math::pi / 180.00)); sphereImage->TransformIndexToPhysicalPoint(sphereIndex, spherePoint); inputImage->TransformPhysicalPointToIndex(spherePoint, regionIndex); if(region.IsInside(regionIndex) == true) { isInside = true; } else { isInside = false; return isInside; } // Point 5/6 sphereIndex[0] = -(sphereCenterIndex[0] + radius * sin(45.00 * itk::Math::pi / 180.00) * cos(45.00 * itk::Math::pi / 180.00)); sphereImage->TransformIndexToPhysicalPoint(sphereIndex, spherePoint); inputImage->TransformPhysicalPointToIndex(spherePoint, regionIndex); if(region.IsInside(regionIndex) == true) { isInside = true; } else { isInside = false; return isInside; } // Point 7/8 sphereIndex[2] = sphereCenterIndex[2] + radius * cos(45.00 * itk::Math::pi / 180.00); sphereImage->TransformIndexToPhysicalPoint(sphereIndex, spherePoint); inputImage->TransformPhysicalPointToIndex(spherePoint, regionIndex); if(region.IsInside(regionIndex) == true) { isInside = true; } else { isInside = false; return isInside; } addOrSub = true; } return isInside; } template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator::InternalCalculateMaskFromPlanarFigure( const itk::Image< TPixel, VImageDimension > *image, unsigned int axis ) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::CastImageFilter< ImageType, MaskImage2DType > CastFilterType; // Generate mask image as new image with same header as input image and // initialize with "1". typename CastFilterType::Pointer castFilter = CastFilterType::New(); castFilter->SetInput( image ); castFilter->Update(); castFilter->GetOutput()->FillBuffer( 1 ); // all PolylinePoints of the PlanarFigure are stored in a vtkPoints object. // These points are used by the vtkLassoStencilSource to create // a vtkImageStencil. const mitk::Geometry2D *planarFigureGeometry2D = m_PlanarFigure->GetGeometry2D(); const typename PlanarFigure::PolyLineType planarFigurePolyline = m_PlanarFigure->GetPolyLine( 0 ); const mitk::Geometry3D *imageGeometry3D = m_Image->GetGeometry( 0 ); // Determine x- and y-dimensions depending on principal axis int i0, i1; switch ( axis ) { case 0: i0 = 1; i1 = 2; break; case 1: i0 = 0; i1 = 2; break; case 2: default: i0 = 0; i1 = 1; break; } m_PlanarFigureCoordinate0= i0; m_PlanarFigureCoordinate1= i1; // store the polyline contour as vtkPoints object bool outOfBounds = false; vtkSmartPointer points = vtkSmartPointer::New(); typename PlanarFigure::PolyLineType::const_iterator it; for ( it = planarFigurePolyline.begin(); it != planarFigurePolyline.end(); ++it ) { Point3D point3D; // Convert 2D point back to the local index coordinates of the selected // image planarFigureGeometry2D->Map( it->Point, point3D ); // Polygons (partially) outside of the image bounds can not be processed // further due to a bug in vtkPolyDataToImageStencil if ( !imageGeometry3D->IsInside( point3D ) ) { outOfBounds = true; } imageGeometry3D->WorldToIndex( point3D, point3D ); points->InsertNextPoint( point3D[i0], point3D[i1], 0 ); } // mark a malformed 2D planar figure ( i.e. area = 0 ) as out of bounds // this can happen when all control points of a rectangle lie on the same line = two of the three extents are zero double bounds[6] = {0, 0, 0, 0, 0, 0}; points->GetBounds( bounds ); bool extent_x = (fabs(bounds[0] - bounds[1])) < mitk::eps; bool extent_y = (fabs(bounds[2] - bounds[3])) < mitk::eps; bool extent_z = (fabs(bounds[4] - bounds[5])) < mitk::eps; // throw an exception if a closed planar figure is deformed, i.e. has only one non-zero extent if ( m_PlanarFigure->IsClosed() && ((extent_x && extent_y) || (extent_x && extent_z) || (extent_y && extent_z))) { mitkThrow() << "Figure has a zero area and cannot be used for masking."; } if ( outOfBounds ) { throw std::runtime_error( "Figure at least partially outside of image bounds!" ); } // create a vtkLassoStencilSource and set the points of the Polygon vtkSmartPointer lassoStencil = vtkSmartPointer::New(); lassoStencil->SetShapeToPolygon(); lassoStencil->SetPoints( points ); // Export from ITK to VTK (to use a VTK filter) typedef itk::VTKImageImport< MaskImage2DType > ImageImportType; typedef itk::VTKImageExport< MaskImage2DType > ImageExportType; typename ImageExportType::Pointer itkExporter = ImageExportType::New(); itkExporter->SetInput( castFilter->GetOutput() ); vtkSmartPointer vtkImporter = vtkSmartPointer::New(); this->ConnectPipelines( itkExporter, vtkImporter ); // Apply the generated image stencil to the input image vtkSmartPointer imageStencilFilter = vtkSmartPointer::New(); imageStencilFilter->SetInputConnection( vtkImporter->GetOutputPort() ); imageStencilFilter->SetStencil( lassoStencil->GetOutput() ); imageStencilFilter->ReverseStencilOff(); imageStencilFilter->SetBackgroundValue( 0 ); imageStencilFilter->Update(); // Export from VTK back to ITK vtkSmartPointer vtkExporter = vtkImageExport::New(); // TODO: this is WRONG, should be vtkSmartPointer::New(), but bug # 14455 vtkExporter->SetInputConnection( imageStencilFilter->GetOutputPort() ); vtkExporter->Update(); typename ImageImportType::Pointer itkImporter = ImageImportType::New(); this->ConnectPipelines( vtkExporter, itkImporter ); itkImporter->Update(); // Store mask m_InternalImageMask2D = itkImporter->GetOutput(); } void ImageStatisticsCalculator::UnmaskedStatisticsProgressUpdate() { // Need to throw away every second progress event to reach a final count of // 100 since two consecutive filters are used in this case static int updateCounter = 0; if ( updateCounter++ % 2 == 0 ) { this->InvokeEvent( itk::ProgressEvent() ); } } void ImageStatisticsCalculator::MaskedStatisticsProgressUpdate() { this->InvokeEvent( itk::ProgressEvent() ); } } diff --git a/Modules/ImageStatistics/mitkImageStatisticsCalculator.h b/Modules/ImageStatistics/mitkImageStatisticsCalculator.h index 0c0b8adc70..1c26326416 100644 --- a/Modules/ImageStatistics/mitkImageStatisticsCalculator.h +++ b/Modules/ImageStatistics/mitkImageStatisticsCalculator.h @@ -1,384 +1,386 @@ /*=================================================================== 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 _MITK_IMAGESTATISTICSCALCULATOR_H #define _MITK_IMAGESTATISTICSCALCULATOR_H #include #include "ImageStatisticsExports.h" #include #include #ifndef __itkHistogram_h #include #endif #include "mitkImage.h" #include "mitkImageTimeSelector.h" #include "mitkPlanarFigure.h" #include namespace mitk { /** * \brief Class for calculating statistics and histogram for an (optionally * masked) image. * * Images can be masked by either a label image (of the same dimensions as * the original image) or by a closed mitk::PlanarFigure, e.g. a circle or * polygon. When masking with a planar figure, the slice corresponding to the * plane containing the figure is extracted and then clipped with contour * defined by the figure. Planar figures need to be aligned along the main axes * of the image (axial, sagittal, coronal). Planar figures on arbitrary * rotated planes are not supported. * * For each operating mode (no masking, masking by image, masking by planar * figure), the calculated statistics and histogram are cached so that, when * switching back and forth between operation modes without modifying mask or * image, the information doesn't need to be recalculated. * * Note: currently time-resolved and multi-channel pictures are not properly * supported. */ class ImageStatistics_EXPORT ImageStatisticsCalculator : public itk::Object { public: enum { MASKING_MODE_NONE = 0, MASKING_MODE_IMAGE, MASKING_MODE_PLANARFIGURE }; typedef itk::Statistics::Histogram HistogramType; typedef HistogramType::ConstIterator HistogramConstIteratorType; struct Statistics { int Label; unsigned int N; double Min; double Max; double Mean; double Median; double Variance; double Sigma; double RMS; double HotspotMin; double HotspotMax; double HotspotMean; double HotspotSigma; double HotspotPeak; vnl_vector< int > MinIndex; vnl_vector< int > MaxIndex; vnl_vector HotspotMaxIndex; vnl_vector HotspotMinIndex; vnl_vector HotspotPeakIndex; void Reset() { Label = 0; N = 0; Min = 0.0; Max = 0.0; Mean = 0.0; Median = 0.0; Variance = 0.0; Sigma = 0.0; RMS = 0.0; HotspotMin = 0.0; HotspotMax = 0.0; HotspotMean = 0.0; HotspotPeak = 0.0; HotspotSigma = 0.0; } }; struct MinMaxIndex { double Max; double Min; vnl_vector MaxIndex; vnl_vector MinIndex; + std::list< vnl_vector > MaxIndexList; + std::list< vnl_vector > MinIndexList; }; typedef std::vector< HistogramType::ConstPointer > HistogramContainer; typedef std::vector< Statistics > StatisticsContainer; mitkClassMacro( ImageStatisticsCalculator, itk::Object ); itkNewMacro( ImageStatisticsCalculator ); /** \brief Set image from which to compute statistics. */ void SetImage( const mitk::Image *image ); /** \brief Set image for masking. */ void SetImageMask( const mitk::Image *imageMask ); /** \brief Set planar figure for masking. */ void SetPlanarFigure( mitk::PlanarFigure *planarFigure ); /** \brief Set/Get operation mode for masking */ void SetMaskingMode( unsigned int mode ); /** \brief Set/Get operation mode for masking */ itkGetMacro( MaskingMode, unsigned int ); /** \brief Set/Get operation mode for masking */ void SetMaskingModeToNone(); /** \brief Set/Get operation mode for masking */ void SetMaskingModeToImage(); /** \brief Set/Get operation mode for masking */ void SetMaskingModeToPlanarFigure(); /** \brief Set a pixel value for pixels that will be ignored in the statistics */ void SetIgnorePixelValue(double value); /** \brief Get the pixel value for pixels that will be ignored in the statistics */ double GetIgnorePixelValue(); /** \brief Set wether a pixel value should be ignored in the statistics */ void SetDoIgnorePixelValue(bool doit); /** \brief Get wether a pixel value will be ignored in the statistics */ bool GetDoIgnorePixelValue(); void SetHotspotSize (double hotspotRadiusInMM); double GetHotspotSize(); void SetCalculateHotspot(bool calculateHotspot); bool IsHotspotCalculated(); /** \brief Compute statistics (together with histogram) for the current * masking mode. * * Computation is not executed if statistics is already up to date. In this * case, false is returned; otherwise, true.*/ virtual bool ComputeStatistics( unsigned int timeStep = 0 ); /** \brief Retrieve the histogram depending on the current masking mode. * * \param label The label for which to retrieve the histogram in multi-label situations (ascending order). */ const HistogramType *GetHistogram( unsigned int timeStep = 0, unsigned int label = 0 ) const; /** \brief Retrieve the histogram depending on the current masking mode (for all image labels. */ const HistogramContainer &GetHistogramVector( unsigned int timeStep = 0 ) const; /** \brief Retrieve statistics depending on the current masking mode. * * \param label The label for which to retrieve the statistics in multi-label situations (ascending order). */ const Statistics &GetStatistics( unsigned int timeStep = 0, unsigned int label = 0 ) const; /** \brief Retrieve statistics depending on the current masking mode. * TODO: Kommentare anpassen! * \param label The label for which to retrieve the statistics in multi-label situations (ascending order). */ const Statistics &GetHotspotStatistics( unsigned int timeStep = 0, unsigned int label = 0 ) const; /** \brief Retrieve statistics depending on the current masking mode (for all image labels). */ const StatisticsContainer &GetStatisticsVector( unsigned int timeStep = 0 ) const; protected: typedef std::vector< HistogramContainer > HistogramVector; typedef std::vector< StatisticsContainer > StatisticsVector; typedef std::vector< itk::TimeStamp > TimeStampVectorType; typedef std::vector< bool > BoolVectorType; typedef itk::Image< unsigned short, 3 > MaskImage3DType; typedef itk::Image< unsigned short, 2 > MaskImage2DType; ImageStatisticsCalculator(); virtual ~ImageStatisticsCalculator(); /** \brief Depending on the masking mode, the image and mask from which to * calculate statistics is extracted from the original input image and mask * data. * * For example, a when using a PlanarFigure as mask, the 2D image slice * corresponding to the PlanarFigure will be extracted from the original * image. If masking is disabled, the original image is simply passed * through. */ void ExtractImageAndMask( unsigned int timeStep = 0 ); /** \brief If the passed vector matches any of the three principal axes * of the passed geometry, the ínteger value corresponding to the axis * is set and true is returned. */ bool GetPrincipalAxis( const Geometry3D *geometry, Vector3D vector, unsigned int &axis ); template < typename TPixel, unsigned int VImageDimension > void InternalCalculateStatisticsUnmasked( const itk::Image< TPixel, VImageDimension > *image, StatisticsContainer* statisticsContainer, HistogramContainer *histogramContainer ); template < typename TPixel, unsigned int VImageDimension > void InternalCalculateStatisticsMasked( const itk::Image< TPixel, VImageDimension > *image, itk::Image< unsigned short, VImageDimension > *maskImage, StatisticsContainer* statisticsContainer, HistogramContainer* histogramContainer ); template < typename TPixel, unsigned int VImageDimension > void InternalCalculateMaskFromPlanarFigure( const itk::Image< TPixel, VImageDimension > *image, unsigned int axis ); template < typename TPixel, unsigned int VImageDimension > void InternalMaskIgnoredPixels( const itk::Image< TPixel, VImageDimension > *image, itk::Image< unsigned short, VImageDimension > *maskImage ); template MinMaxIndex CalculateMinMaxIndex( const itk::Image *inputImage, itk::Image *maskImage, float minBoundary, float maxBoundary); template < typename TPixel, unsigned int VImageDimension> Statistics CalculateHotspotStatistics( - itk::Image *inputImage, + const itk::Image *inputImage, itk::Image *maskImage, double RadiusInMM); template < typename TPixel, unsigned int VImageDimension> bool IsSphereInsideRegion( const itk::Image *inputImage, itk::Image *sphereImage); /** Connection from ITK to VTK */ template void ConnectPipelines(ITK_Exporter exporter, vtkSmartPointer importer) { importer->SetUpdateInformationCallback(exporter->GetUpdateInformationCallback()); importer->SetPipelineModifiedCallback(exporter->GetPipelineModifiedCallback()); importer->SetWholeExtentCallback(exporter->GetWholeExtentCallback()); importer->SetSpacingCallback(exporter->GetSpacingCallback()); importer->SetOriginCallback(exporter->GetOriginCallback()); importer->SetScalarTypeCallback(exporter->GetScalarTypeCallback()); importer->SetNumberOfComponentsCallback(exporter->GetNumberOfComponentsCallback()); importer->SetPropagateUpdateExtentCallback(exporter->GetPropagateUpdateExtentCallback()); importer->SetUpdateDataCallback(exporter->GetUpdateDataCallback()); importer->SetDataExtentCallback(exporter->GetDataExtentCallback()); importer->SetBufferPointerCallback(exporter->GetBufferPointerCallback()); importer->SetCallbackUserData(exporter->GetCallbackUserData()); } /** Connection from VTK to ITK */ template void ConnectPipelines(vtkSmartPointer exporter, ITK_Importer importer) { importer->SetUpdateInformationCallback(exporter->GetUpdateInformationCallback()); importer->SetPipelineModifiedCallback(exporter->GetPipelineModifiedCallback()); importer->SetWholeExtentCallback(exporter->GetWholeExtentCallback()); importer->SetSpacingCallback(exporter->GetSpacingCallback()); importer->SetOriginCallback(exporter->GetOriginCallback()); importer->SetScalarTypeCallback(exporter->GetScalarTypeCallback()); importer->SetNumberOfComponentsCallback(exporter->GetNumberOfComponentsCallback()); importer->SetPropagateUpdateExtentCallback(exporter->GetPropagateUpdateExtentCallback()); importer->SetUpdateDataCallback(exporter->GetUpdateDataCallback()); importer->SetDataExtentCallback(exporter->GetDataExtentCallback()); importer->SetBufferPointerCallback(exporter->GetBufferPointerCallback()); importer->SetCallbackUserData(exporter->GetCallbackUserData()); } void UnmaskedStatisticsProgressUpdate(); void MaskedStatisticsProgressUpdate(); /** m_Image contains the input image (e.g. 2D, 3D, 3D+t)*/ mitk::Image::ConstPointer m_Image; mitk::Image::ConstPointer m_ImageMask; mitk::PlanarFigure::Pointer m_PlanarFigure; HistogramVector m_ImageHistogramVector; HistogramVector m_MaskedImageHistogramVector; HistogramVector m_PlanarFigureHistogramVector; HistogramType::Pointer m_EmptyHistogram; HistogramContainer m_EmptyHistogramContainer; StatisticsVector m_ImageStatisticsVector; StatisticsVector m_MaskedImageStatisticsVector; StatisticsVector m_PlanarFigureStatisticsVector; StatisticsVector m_MaskedImageHotspotStatisticsVector; Statistics m_EmptyStatistics; StatisticsContainer m_EmptyStatisticsContainer; unsigned int m_MaskingMode; bool m_MaskingModeChanged; /** m_InternalImage contains a image volume at one time step (e.g. 2D, 3D)*/ mitk::Image::ConstPointer m_InternalImage; MaskImage3DType::Pointer m_InternalImageMask3D; MaskImage2DType::Pointer m_InternalImageMask2D; TimeStampVectorType m_ImageStatisticsTimeStampVector; TimeStampVectorType m_MaskedImageStatisticsTimeStampVector; TimeStampVectorType m_PlanarFigureStatisticsTimeStampVector; BoolVectorType m_ImageStatisticsCalculationTriggerVector; BoolVectorType m_MaskedImageStatisticsCalculationTriggerVector; BoolVectorType m_PlanarFigureStatisticsCalculationTriggerVector; double m_IgnorePixelValue; bool m_DoIgnorePixelValue; bool m_IgnorePixelValueChanged; double m_HotspotSize; bool m_CalculateHotspot; unsigned int m_PlanarFigureAxis; // Normal axis for PlanarFigure unsigned int m_PlanarFigureSlice; // Slice which contains PlanarFigure int m_PlanarFigureCoordinate0; // First plane-axis for PlanarFigure int m_PlanarFigureCoordinate1; // Second plane-axis for PlanarFigure }; } #endif // #define _MITK_IMAGESTATISTICSCALCULATOR_H