diff --git a/Modules/ImageStatistics/Testing/mitkImageStatisticsTextureAnalysisTest.cpp b/Modules/ImageStatistics/Testing/mitkImageStatisticsTextureAnalysisTest.cpp index 11f751fa84..039239e967 100644 --- a/Modules/ImageStatistics/Testing/mitkImageStatisticsTextureAnalysisTest.cpp +++ b/Modules/ImageStatistics/Testing/mitkImageStatisticsTextureAnalysisTest.cpp @@ -1,235 +1,234 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkTestingMacros.h" #include "itkImage.h" #include "mitkExtendedLabelStatisticsImageFilter.h" #include "mitkExtendedStatisticsImageFilter.h" #include "mitkNumericConstants.h" /** \section Testing of Skewness and Kurtosis * This test class is for testing the added coefficients Skewness and Kurtosis * for the mitkExtendedLabelStatisticsImageFilter (Masked Images) and for the * mitkExtendedStatisticsImageFilter (Unmasked Images). Both filter will be tested * against two pictures. */ class mitkImageStatisticsTextureAnalysisTestClass { /** * \brief Explanation of the mitkImageStatisticsTextureAnalysisTestClass test class * * this test class produce test images and masking images with the method CreatingTestImageForDifferentLabelSize. * TestInstanceFortheMaskedStatisticsFilter and TestInstanceFortheUnmaskedStatisticsFilter are the two Instances * for the filters of masking and unmasking images. * TestofSkewnessAndKurtosisForMaskedImagesand TestofSkewnessAndKurtosisForUnmaskedImages are the correlated test * for the checking the values. */ public: typedef itk::Image< int,3 >ImageType; typedef ImageType::Pointer PointerOfImage; typedef itk::ExtendedLabelStatisticsImageFilter< ImageType, ImageType > LabelStatisticsFilterType; typedef LabelStatisticsFilterType::Pointer labelStatisticsFilterPointer; typedef itk::ExtendedStatisticsImageFilter< ImageType > StatisticsFilterType; typedef StatisticsFilterType::Pointer StatisticsFilterPointer; ImageType::Pointer CreatingTestImageForDifferentLabelSize( int factorOfDividingThePicture, int bufferValue, int labelValue) { ImageType::Pointer image = ImageType :: New(); ImageType::IndexType start; ImageType::SizeType size; start[0] = 0; start[1] = 0; start[2] = 0; size[0] = 100; size[1] = 100; size[2] = 100; ImageType:: RegionType region; region.SetSize( size ); region.SetIndex( start ); image->SetRegions(region); image->Allocate(); image->FillBuffer(bufferValue); - int i = 0; for(unsigned int r = 0; r < 50; r++) { for(int c = 0; c < factorOfDividingThePicture; c++) { for(unsigned int l = 0; l < 100; l++) { ImageType::IndexType pixelIndex; pixelIndex[0] = r; pixelIndex[1] = c; pixelIndex[2] = l; image->SetPixel(pixelIndex, labelValue); } } } return image; } LabelStatisticsFilterType::Pointer TestInstanceFortheMaskedStatisticsFilter(ImageType::Pointer image, ImageType::Pointer maskImage) { LabelStatisticsFilterType::Pointer labelStatisticsFilter; labelStatisticsFilter = LabelStatisticsFilterType::New(); labelStatisticsFilter->SetInput( image ); labelStatisticsFilter->UseHistogramsOn(); labelStatisticsFilter->SetHistogramParameters( 20, -10, 10); labelStatisticsFilter->SetLabelInput( maskImage ); labelStatisticsFilter->Update(); return labelStatisticsFilter; } StatisticsFilterType::Pointer TestInstanceFortheUnmaskedStatisticsFilter(ImageType::Pointer image ) { StatisticsFilterType::Pointer StatisticsFilter; StatisticsFilter = StatisticsFilterType::New(); StatisticsFilter->SetInput( image ); StatisticsFilter->SetHistogramParameters( 20, -10, 10 ); StatisticsFilter->Update(); return StatisticsFilter; } //test for Skewness,Kurtosis and MPP for masked Images void TestofSkewnessKurtosisAndMPPForMaskedImages(LabelStatisticsFilterType::Pointer labelStatisticsFilter, double expectedSkewness, double expectedKurtosis, double expectedMPP) { // let's create an object of our class bool isSkewsnessLowerlimitCorrect = labelStatisticsFilter->GetSkewness( 1 )- expectedKurtosis+ std::pow(10,-3) <= expectedSkewness; bool isSkewsnessUpperlimitCorrect = labelStatisticsFilter->GetSkewness( 1 )+ expectedKurtosis+ std::pow(10,-3) >= expectedSkewness; MITK_TEST_CONDITION( isSkewsnessLowerlimitCorrect && isSkewsnessUpperlimitCorrect,"expectedSkewness: " << expectedSkewness << " actual Value: " << labelStatisticsFilter->GetSkewness( 1 ) ); bool isKurtosisUpperlimitCorrect = labelStatisticsFilter->GetKurtosis( 1 ) <= expectedKurtosis+ std::pow(10,-3); bool isKurtosisLowerlimitCorrect = expectedKurtosis- std::pow(10,-3) <= labelStatisticsFilter->GetKurtosis( 1 ); MITK_TEST_CONDITION( isKurtosisUpperlimitCorrect && isKurtosisLowerlimitCorrect,"expectedKurtosis: " << expectedKurtosis << " actual Value: " << labelStatisticsFilter->GetKurtosis( 1 ) ); MITK_TEST_CONDITION( ( expectedMPP - labelStatisticsFilter->GetMPP( 1 ) ) < 1, "expected MPP: " << expectedMPP << " actual Value: " << labelStatisticsFilter->GetMPP( 1 ) ); } //test for Entropy,Uniformity and UPP for masked Images void TestofEntropyUniformityAndUppForMaskedImages(LabelStatisticsFilterType::Pointer labelStatisticsFilter, double expectedEntropy, double expectedUniformity, double expectedUPP) { bool calculatedEntropyLowerLimit = labelStatisticsFilter->GetEntropy( 1 ) >= expectedEntropy - std::pow(10,-3); bool calculatedUniformityLowerLimit = labelStatisticsFilter->GetUniformity( 1 ) >= expectedUniformity - std::pow(10,-3); bool calculatedUppLowerLimit = labelStatisticsFilter->GetUPP( 1 ) >= expectedUPP - std::pow(10,-3); bool calculatedEntropyUpperLimit = labelStatisticsFilter->GetEntropy( 1 ) <= expectedEntropy + std::pow(10,-3); bool calculatedUniformityUpperLimit = labelStatisticsFilter->GetUniformity( 1 ) <= expectedUniformity + std::pow(10,-3); bool calculatedUppUpperLimit = labelStatisticsFilter->GetUPP( 1 ) <= expectedUPP + std::pow(10,-3); MITK_TEST_CONDITION( calculatedEntropyLowerLimit && calculatedEntropyUpperLimit, "expected Entropy: " << expectedEntropy << " actual Value: " << labelStatisticsFilter->GetEntropy( 1 ) ); MITK_TEST_CONDITION( calculatedUniformityLowerLimit && calculatedUniformityUpperLimit, "expected Uniformity: " << expectedUniformity << " actual Value: " << labelStatisticsFilter->GetUniformity( 1 ) ); MITK_TEST_CONDITION( calculatedUppLowerLimit && calculatedUppUpperLimit, "expected UPP: " << expectedUPP << " actual Value: " << labelStatisticsFilter->GetUPP( 1 ) ); } //test for Skewness,Kurtosis and MPP for unmasked Images void TestofSkewnessKurtosisAndMPPForUnmaskedImages(StatisticsFilterType::Pointer StatisticsFilter, double expectedSkewness, double expectedKurtosis, double expectedMPP) { // let's create an object of our class bool isSkewsnessLowerlimitCorrect = StatisticsFilter->GetSkewness()- expectedKurtosis+ std::pow(10,-3) <= expectedSkewness; bool isSkewsnessUpperlimitCorrect = StatisticsFilter->GetSkewness()+ expectedKurtosis+ std::pow(10,-3) >= expectedSkewness; MITK_TEST_CONDITION( isSkewsnessLowerlimitCorrect && isSkewsnessUpperlimitCorrect,"expectedSkewness: " << expectedSkewness << " actual Value: " << StatisticsFilter->GetSkewness() ); bool isKurtosisUpperlimitCorrect = StatisticsFilter->GetKurtosis() <= expectedKurtosis+ std::pow(10,-3); bool isKurtosisLowerlimitCorrect = expectedKurtosis- std::pow(10,-3) <= StatisticsFilter->GetKurtosis(); MITK_TEST_CONDITION( isKurtosisUpperlimitCorrect && isKurtosisLowerlimitCorrect,"expectedKurtosis: " << expectedKurtosis << " actual Value: " << StatisticsFilter->GetKurtosis() ); MITK_TEST_CONDITION( ( expectedMPP - StatisticsFilter->GetMPP() ) < mitk::eps, "expected MPP: " << expectedMPP << " actual Value: " << StatisticsFilter->GetMPP() ); } //test for Entropy,Uniformity and UPP for unmasked Images void TestofEntropyUniformityAndUppForUnmaskedImages(StatisticsFilterType::Pointer StatisticsFilter, double expectedEntropy, double expectedUniformity, double expectedUPP) { bool calculatedEntropyLowerLimit = StatisticsFilter->GetEntropy() >= expectedEntropy - std::pow(10,-3); bool calculatedUniformityLowerLimit = StatisticsFilter->GetUniformity() >= expectedUniformity - std::pow(10,-3); bool calculatedUppLowerLimit = StatisticsFilter->GetUPP() >= expectedUPP - std::pow(10,-3); bool calculatedEntropyUpperLimit = StatisticsFilter->GetEntropy() <= expectedEntropy + std::pow(10,-3); bool calculatedUniformityUpperLimit = StatisticsFilter->GetUniformity() <= expectedUniformity + std::pow(10,-3); bool calculatedUppUpperLimit = StatisticsFilter->GetUPP() <= expectedUPP + std::pow(10,-3); MITK_TEST_CONDITION( calculatedEntropyLowerLimit && calculatedEntropyUpperLimit, "expected Entropy: " << expectedEntropy << " actual Value: " << StatisticsFilter->GetEntropy() ); MITK_TEST_CONDITION( calculatedUniformityLowerLimit && calculatedUniformityUpperLimit, "expected Uniformity: " << expectedUniformity << " actual Value: " << StatisticsFilter->GetUniformity() ); MITK_TEST_CONDITION( calculatedUppLowerLimit && calculatedUppUpperLimit, "expected UPP: " << expectedUPP << " actual Value: " << StatisticsFilter->GetUPP() ); } }; int mitkImageStatisticsTextureAnalysisTest(int, char* []) { // always start with this! MITK_TEST_BEGIN("mitkImageStatisticsTextureAnalysisTest") mitkImageStatisticsTextureAnalysisTestClass testclassInstance; mitkImageStatisticsTextureAnalysisTestClass::PointerOfImage labelImage = testclassInstance.CreatingTestImageForDifferentLabelSize(100, 1, 1); mitkImageStatisticsTextureAnalysisTestClass::PointerOfImage image = testclassInstance.CreatingTestImageForDifferentLabelSize(100, 3, 2); mitkImageStatisticsTextureAnalysisTestClass::PointerOfImage image2 = testclassInstance.CreatingTestImageForDifferentLabelSize(50, 3, 2); //test for masked images mitkImageStatisticsTextureAnalysisTestClass::labelStatisticsFilterPointer mitkLabelFilter= testclassInstance.TestInstanceFortheMaskedStatisticsFilter( image,labelImage); testclassInstance.TestofSkewnessKurtosisAndMPPForMaskedImages(mitkLabelFilter, 0, 0.999998, 2.5); testclassInstance.TestofEntropyUniformityAndUppForMaskedImages(mitkLabelFilter, 1, 0.5, 0.5); mitkImageStatisticsTextureAnalysisTestClass::labelStatisticsFilterPointer mitkLabelFilter2= testclassInstance.TestInstanceFortheMaskedStatisticsFilter( image2,labelImage); testclassInstance.TestofSkewnessKurtosisAndMPPForMaskedImages(mitkLabelFilter2, -1.1547, 2.33333, 2.75); testclassInstance.TestofEntropyUniformityAndUppForMaskedImages(mitkLabelFilter2, 0.811278, 0.625, 0.625); //test for unmasked images mitkImageStatisticsTextureAnalysisTestClass::StatisticsFilterPointer mitkFilter= testclassInstance.TestInstanceFortheUnmaskedStatisticsFilter( image); testclassInstance.TestofSkewnessKurtosisAndMPPForUnmaskedImages(mitkFilter, 0, 0.999998, 2.5); testclassInstance.TestofEntropyUniformityAndUppForUnmaskedImages(mitkFilter, 1, 0.5, 0.5); mitkImageStatisticsTextureAnalysisTestClass::StatisticsFilterPointer mitkFilter2= testclassInstance.TestInstanceFortheUnmaskedStatisticsFilter( image2); testclassInstance.TestofSkewnessKurtosisAndMPPForUnmaskedImages(mitkFilter2, -1.1547, 2.33333, 2.75); testclassInstance.TestofEntropyUniformityAndUppForUnmaskedImages( mitkFilter2, 0.811278, 0.625, 0.625); MITK_TEST_END() } diff --git a/Modules/ImageStatistics/itkMultiGaussianImageSource.hxx b/Modules/ImageStatistics/itkMultiGaussianImageSource.hxx index 91d77316c1..71a4da4114 100644 --- a/Modules/ImageStatistics/itkMultiGaussianImageSource.hxx +++ b/Modules/ImageStatistics/itkMultiGaussianImageSource.hxx @@ -1,1004 +1,1003 @@ /*========================================================================= * * Copyright Insight Software Consortium * * Licensed under the Apache License, Version 2.0 (the "License"); * you may not use this file except in compliance with the License. * You may obtain a copy of the License at * * http://www.apache.org/licenses/LICENSE-2.0.txt * * Unless required by applicable law or agreed to in writing, software * distributed under the License is distributed on an "AS IS" BASIS, * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. * See the License for the specific language governing permissions and * limitations under the License. * *=========================================================================*/ /*========================================================================= * * Portions of this file are subject to the VTK Toolkit Version 3 copyright. * * Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen * * For complete copyright, license and disclaimer of warranty information * please refer to the NOTICE file at the top of the ITK source tree. * *=========================================================================*/ #ifndef __itkMultiGaussianImageSource_hxx #define __itkMultiGaussianImageSource_hxx #include #include #include #include "itkMultiGaussianImageSource.h" #include "itkImageRegionIterator.h" #include "itkObjectFactory.h" #include "itkProgressReporter.h" #include "itkDOMNodeXMLWriter.h" #include "stdlib.h" namespace itk { /** * */ template< class TOutputImage > MultiGaussianImageSource< TOutputImage > ::MultiGaussianImageSource() { //Initial image is 100 wide in each direction. for ( unsigned int i = 0; i < TOutputImage::GetImageDimension(); i++ ) { m_Size[i] = 100; m_Spacing[i] = 1.0; m_Origin[i] = 0.0; m_SphereMidpoint[i] = 0; } m_NumberOfGaussians = 0; m_Radius = 1; m_RadiusStepNumber = 5; m_MeanValue = 0; m_Min = NumericTraits< OutputImagePixelType >::NonpositiveMin(); m_Max = NumericTraits< OutputImagePixelType >::max(); } //----------------------------------------------------------------------------------------------------------------------- template< class TOutputImage > MultiGaussianImageSource< TOutputImage > ::~MultiGaussianImageSource() {} //----------------------------------------------------------------------------------------------------------------------- template< class TOutputImage > void MultiGaussianImageSource< TOutputImage > ::SetSize(SizeValueArrayType sizeArray) { const unsigned int count = TOutputImage::ImageDimension; unsigned int i; for ( i = 0; i < count; i++ ) { if ( sizeArray[i] != this->m_Size[i] ) { break; } } if ( i < count ) { this->Modified(); for ( i = 0; i < count; i++ ) { this->m_Size[i] = sizeArray[i]; } } } template< class TOutputImage > const typename MultiGaussianImageSource< TOutputImage >::SizeValueType * MultiGaussianImageSource< TOutputImage > ::GetSize() const { return this->m_Size.GetSize(); } //----------------------------------------------------------------------------------------------------------------------- template< class TOutputImage > void MultiGaussianImageSource< TOutputImage > ::SetSpacing(SpacingValueArrayType spacingArray) { const unsigned int count = TOutputImage::ImageDimension; unsigned int i; for ( i = 0; i < count; i++ ) { if ( spacingArray[i] != this->m_Spacing[i] ) { break; } } if ( i < count ) { this->Modified(); for ( i = 0; i < count; i++ ) { this->m_Spacing[i] = spacingArray[i]; } } } //----------------------------------------------------------------------------------------------------------------------- template< class TOutputImage > void MultiGaussianImageSource< TOutputImage > ::SetOrigin(PointValueArrayType originArray) { const unsigned int count = TOutputImage::ImageDimension; unsigned int i; for ( i = 0; i < count; i++ ) { if ( originArray[i] != this->m_Origin[i] ) { break; } } if ( i < count ) { this->Modified(); for ( i = 0; i < count; i++ ) { this->m_Origin[i] = originArray[i]; } } } //----------------------------------------------------------------------------------------------------------------------- template< class TOutputImage > const typename MultiGaussianImageSource< TOutputImage >::PointValueType * MultiGaussianImageSource< TOutputImage > ::GetOrigin() const { for ( unsigned int i = 0; i < TOutputImage::ImageDimension; i++ ) { this->m_OriginArray[i] = this->m_Origin[i]; } return this->m_OriginArray; } template< class TOutputImage > const typename MultiGaussianImageSource< TOutputImage >::SpacingValueType * MultiGaussianImageSource< TOutputImage > ::GetSpacing() const { for ( unsigned int i = 0; i < TOutputImage::ImageDimension; i++ ) { this->m_SpacingArray[i] = this->m_Spacing[i]; } return this->m_SpacingArray; } //----------------------------------------------------------------------------------------------------------------------- /** * */ template< class TOutputImage > void MultiGaussianImageSource< TOutputImage > ::PrintSelf(std::ostream & os, Indent indent) const { Superclass::PrintSelf(os, indent); os << indent << "Max: " << static_cast< typename NumericTraits< OutputImagePixelType >::PrintType >( m_Max ) << std::endl; os << indent << "Min: " << static_cast< typename NumericTraits< OutputImagePixelType >::PrintType >( m_Min ) << std::endl; os << indent << "Origin: ["; unsigned int ii = 0; while( ii < TOutputImage::ImageDimension - 1 ) { os << m_Origin[ii] << ", "; ++ii; } os << m_Origin[ii] << "]" << std::endl; os << indent << "Spacing: ["; ii = 0; while( ii < TOutputImage::ImageDimension - 1 ) { os << m_Spacing[ii] << ", "; ++ii; } os << m_Spacing[ii] << "]" << std::endl; os << indent << "Size: ["; ii = 0; while( ii < TOutputImage::ImageDimension - 1 ) { os << m_Size[ii] << ", "; ++ii; } os << m_Size[ii] << "]" << std::endl; } //----------------------------------------------------------------------------------------------------------------------- template< class TOutputImage > unsigned int MultiGaussianImageSource< TOutputImage > ::GetNumberOfGaussians() const { return this->m_NumberOfGaussians; } //----------------------------------------------------------------------------------------------------------------------- template< class TOutputImage > typename MultiGaussianImageSource< TOutputImage >::RadiusType MultiGaussianImageSource< TOutputImage > ::GetRadius() const { return this->m_Radius; } //----------------------------------------------------------------------------------------------------------------------- template< class TOutputImage > void MultiGaussianImageSource< TOutputImage > ::SetRadius( RadiusType radius ) { this->m_Radius = radius; } //----------------------------------------------------------------------------------------------------------------------- template< class TOutputImage > void MultiGaussianImageSource< TOutputImage > ::SetNumberOfGausssians( unsigned int n ) { this->m_NumberOfGaussians = n; } //----------------------------------------------------------------------------------------------------------------------- template< class TOutputImage > void MultiGaussianImageSource< TOutputImage > ::SetRegionOfInterest( ItkVectorType roiMin, ItkVectorType roiMax ) { m_RegionOfInterestMax = roiMax; m_RegionOfInterestMin = roiMin; } //----------------------------------------------------------------------------------------------------------------------- template< class TOutputImage > const typename MultiGaussianImageSource< TOutputImage >::OutputImagePixelType MultiGaussianImageSource< TOutputImage > ::GetMaxMeanValue() const { return m_MeanValue; } //----------------------------------------------------------------------------------------------------------------------- template< class TOutputImage > const typename MultiGaussianImageSource< TOutputImage >::OutputImagePixelType MultiGaussianImageSource< TOutputImage > ::GetMaxValueInSphere() const { return m_MaxValueInSphere; } //----------------------------------------------------------------------------------------------------------------------- template< class TOutputImage > const typename MultiGaussianImageSource< TOutputImage >::IndexType MultiGaussianImageSource< TOutputImage > ::GetMaxValueIndexInSphere() const { return m_MaxValueIndexInSphere; } //----------------------------------------------------------------------------------------------------------------------- template< class TOutputImage > const typename MultiGaussianImageSource< TOutputImage >::OutputImagePixelType MultiGaussianImageSource< TOutputImage > ::GetMinValueInSphere() const { return m_MinValueInSphere; } //----------------------------------------------------------------------------------------------------------------------- template< class TOutputImage > const typename MultiGaussianImageSource< TOutputImage >::IndexType MultiGaussianImageSource< TOutputImage > ::GetMinValueIndexInSphere() const { return m_MinValueIndexInSphere; } //----------------------------------------------------------------------------------------------------------------------- template< class TOutputImage > const typename MultiGaussianImageSource< TOutputImage >::IndexType MultiGaussianImageSource< TOutputImage > ::GetSphereMidpoint() const { return m_SphereMidpoint; } //----------------------------------------------------------------------------------------------------------------------- /* Calculate and return value of the integral of the gaussian in a cuboid region with the dimension 3: in the x-axis between xMin and xMax and in the y-axis between yMin and yMax and in the z-axis also between zMin and zMax. */ template< class TOutputImage > double MultiGaussianImageSource< TOutputImage > ::MultiGaussianFunctionValueAtCuboid(double xMin, double xMax, double yMin, double yMax, double zMin, double zMax) { double mean = 0; double summand0, summand1, summand2, value, factor; for(unsigned int n = 0; n < m_NumberOfGaussians; ++n) { summand0 = FunctionPhi((xMax - m_CenterX[n]) / m_SigmaX[n] ) - FunctionPhi((xMin - m_CenterX[n]) / m_SigmaX[n] ); summand1 = FunctionPhi((yMax - m_CenterY[n]) / m_SigmaY[n] ) - FunctionPhi((yMin - m_CenterY[n]) / m_SigmaY[n] ); summand2 = FunctionPhi((zMax - m_CenterZ[n]) / m_SigmaZ[n] ) - FunctionPhi((zMin - m_CenterZ[n]) / m_SigmaZ[n] ); value = summand0 * summand1 * summand2; factor = (m_SigmaX[n] * m_SigmaY[n] * m_SigmaZ[n] ) * pow(2.0 * itk::Math::pi, 1.5 ); mean = mean + factor * value * m_Altitude[n]; } return mean; } //--------------------------------------------------------------------------------------------------------------------- /* Returns the linear interpolation of the values of the standard normal distribution function. This values could be seen in the vector m_NormalDistValues. */ template< class TOutputImage > double MultiGaussianImageSource< TOutputImage > ::FunctionPhi(double value) { double phiAtValue; //linear interpolation between the values int indexValue = static_cast( 100 * value); if( indexValue > 409 ) { return phiAtValue = 1.0; } else if( indexValue < -409 ) { return phiAtValue = 0.0; } else if( indexValue == 409 ) { return phiAtValue = m_NormalDistValues[409]; } else if( indexValue == -409 ) { return phiAtValue = 1.0 - m_NormalDistValues[409]; } else { bool negative = false; if (indexValue < 0.0) { negative = true; value = -value; } int indexUp = static_cast( 100 * value) + 1; int indexDown = static_cast( 100 * value); double alpha = 100.0 * value - static_cast(indexDown); phiAtValue = (1.0 - alpha) * m_NormalDistValues[indexDown] + alpha * m_NormalDistValues[indexUp] ; if(negative) { phiAtValue = 1.0 - phiAtValue; } return phiAtValue; } } //---------------------------------------------------------------------------------------------------------------------- /* Set the midpoint of the cuboid in a vector m_Midpoints. This cuboids discretise the sphere with the octree method. Set the radius of the cuboid ( = 0.5 * length of the side of the cuboid ) in the vector m_RadiusCuboid. */ template< class TOutputImage > void MultiGaussianImageSource< TOutputImage > ::InsertPoints( PointType globalCoordinateMidpointCuboid, double cuboidRadius) { typename MapContainerPoints::ElementIdentifier id = m_Midpoints.Size(); m_Midpoints.InsertElement(id, globalCoordinateMidpointCuboid); m_RadiusCuboid.InsertElement(id, cuboidRadius); } //---------------------------------------------------------------------------------------------------------------------- /* This recursive method realise the octree method. It subdivide a cuboid in eight cuboids, when this cuboid crosses the boundary of sphere. If the cuboid is inside the sphere, it calculates the integral and, if uncommented, insert the midpoint of the cuboid in the m_Midpoints vector. */ template< class TOutputImage > void MultiGaussianImageSource< TOutputImage > ::CalculateEdgesInSphere( PointType globalCoordinateMidpointCuboid, PointType globalCoordinateMidpointSphere, double cuboidRadius, int level ) { double xMin, xMax, yMin, yMax, zMin, zMax; double cuboidRadiusRecursion = cuboidRadius; PointType newMidpoint; int intersect = IntersectTheSphere( globalCoordinateMidpointCuboid, globalCoordinateMidpointSphere, cuboidRadiusRecursion); if( intersect == 1 ) { if (level < 4) { // cuboid intersect the sphere -> call CalculateEdgesInSphere (this method) for the subdivided cuboid cuboidRadiusRecursion = cuboidRadiusRecursion / 2.0; for(int i = -1; i < 2; i+=2) { for(int k = -1; k < 2; k+=2) { for(int j = -1; j < 2; j+=2) { newMidpoint[0] = globalCoordinateMidpointCuboid[0] + static_cast(i) * cuboidRadiusRecursion; newMidpoint[1] = globalCoordinateMidpointCuboid[1] + static_cast(k) * cuboidRadiusRecursion; newMidpoint[2] = globalCoordinateMidpointCuboid[2] + static_cast(j) * cuboidRadiusRecursion; this->CalculateEdgesInSphere( newMidpoint, globalCoordinateMidpointSphere, cuboidRadiusRecursion, level + 1 ); } } } } // last step of recursion -> on the boundary else { // Calculate the integral and take the half of it (because we are on the boundary) xMin = globalCoordinateMidpointCuboid[0] - cuboidRadius; xMax = globalCoordinateMidpointCuboid[0] + cuboidRadius; yMin = globalCoordinateMidpointCuboid[1] - cuboidRadius; yMax = globalCoordinateMidpointCuboid[1] + cuboidRadius; zMin = globalCoordinateMidpointCuboid[2] - cuboidRadius; zMax = globalCoordinateMidpointCuboid[2] + cuboidRadius; // size is in index coordinate -> multiply by the spacing -> global coordinate of the image boundary // yz Plane bool yzPlaneNotCrossYSection = xMin >= m_Origin[0] && xMax <= m_Size[0] * m_Spacing[0]; // xz Plane bool xzPlaneNotCrossYSection = yMin >= m_Origin[1] && yMax <= m_Size[1] * m_Spacing[1]; // xy Plane bool xyPlaneNotCrossZSection = zMin >= m_Origin[2] && zMax <= m_Size[2] * m_Spacing[2]; //check if the boundary of the integral is inside the image if( yzPlaneNotCrossYSection && xzPlaneNotCrossYSection && xyPlaneNotCrossZSection) { //double temp = this->MultiGaussianFunctionValueAtCuboid( xMin, xMax, yMin, yMax, zMin, zMax ) * 0.5; m_meanValueTemp = m_meanValueTemp + this->MultiGaussianFunctionValueAtCuboid( xMin, xMax, yMin, yMax, zMin, zMax ) * 0.5; m_Volume = m_Volume + pow( 2.0 * cuboidRadius, 3.0) / 2.0; } } } else if(intersect == 2) { // cuboid in the sphere // To insert the midpoint and the radius of the cuboid in a vector, so that we can visualise the midpoints, uncomment the next line and the line WriteXMLToTestTheCuboidInsideTheSphere(); // InsertPoints(globalCoordinateMidpointCuboid, cuboidRadius); // Calculate the integral boundary xMin = globalCoordinateMidpointCuboid[0] - cuboidRadius; xMax = globalCoordinateMidpointCuboid[0] + cuboidRadius; yMin = globalCoordinateMidpointCuboid[1] - cuboidRadius; yMax = globalCoordinateMidpointCuboid[1] + cuboidRadius; zMin = globalCoordinateMidpointCuboid[2] - cuboidRadius; zMax = globalCoordinateMidpointCuboid[2] + cuboidRadius; // size is in index coordinate -> multiply by the spacing -> global coordinate of the image boundary // yz Plane // bool yzPlaneAtOriginCrossXSection = xMin <= m_Origin[0] && xMax >= m_Origin[0]; // bool yzPlaneAtImageBorderCrossXSection = xMin <= m_Size[0] * m_Spacing[0] && xMax >= m_Size[0] * m_Spacing[0]; bool yzPlaneNotCrossYSection = xMin >= m_Origin[0] && xMax <= m_Size[0] * m_Spacing[0]; // xz Plane // bool xzPlaneAtOriginCrossYSection = yMin <= m_Origin[1] && yMax >= m_Origin[1]; // bool xzPlaneAtImageBorderCrossYSection = yMin <= m_Size[1] * m_Spacing[1] && yMax >= m_Size[1] * m_Spacing[1]; bool xzPlaneNotCrossYSection = yMin >= m_Origin[1] && yMax <= m_Size[1] * m_Spacing[1]; // xy Plane // bool xyPlaneAtOriginCrossZSection = zMin <= m_Origin[2] && zMax >= m_Origin[2]; // bool xyPlaneAtImageBorderCrossZSection = zMin <= m_Size[2] * m_Spacing[2] && zMax >= m_Size[2] * m_Spacing[2]; bool xyPlaneNotCrossZSection = zMin >= m_Origin[2] && zMax <= m_Size[2] * m_Spacing[2]; if( yzPlaneNotCrossYSection && xzPlaneNotCrossYSection && xyPlaneNotCrossZSection) { m_meanValueTemp = m_meanValueTemp + this->MultiGaussianFunctionValueAtCuboid( xMin, xMax, yMin, yMax, zMin, zMax ); m_Volume = m_Volume + pow( 2.0 * cuboidRadius, 3.0); } // check if the boundary of the image intersect the cuboid and if yes, change the limits of the cuboid to be only inside the image; therefor we cut the sphere and neglect the part of it outside the image else /* if( // one plane crosses the cuboid ( (yzPlaneAtOriginCrossXSection && xzPlaneNotCrossYSection && xyPlaneNotCrossZSection) || (yzPlaneAtImageBorderCrossXSection && xzPlaneNotCrossYSection && xyPlaneNotCrossZSection) || (yzPlaneNotCrossYSection && xzPlaneAtOriginCrossYSection && xyPlaneNotCrossZSection) || (yzPlaneNotCrossYSection && xzPlaneAtImageBorderCrossYSection && xyPlaneNotCrossZSection) || (yzPlaneNotCrossYSection && xzPlaneNotCrossYSection && xyPlaneAtOriginCrossZSection) || (yzPlaneNotCrossYSection && xzPlaneNotCrossYSection && xyPlaneAtImageBorderCrossZSection) ) || // two plane cross the cuboid (on the image edges possible) ( (yzPlaneAtOriginCrossXSection && xzPlaneAtOriginCrossYSection && xyPlaneNotCrossZSection) || ( (yzPlaneAtOriginCrossXSection && xzPlaneNotCrossYSection && xyPlaneAtOriginCrossZSection) || (yzPlaneAtImageBorderCrossXSection && xzPlaneAtImageBorderCrossYSection && xyPlaneNotCrossZSection) || (yzPlaneAtImageBorderCrossXSection && xzPlaneNotCrossYSection && xyPlaneAtImageBorderCrossZSection) || (yzPlaneNotCrossYSection && xzPlaneAtOriginCrossYSection && xyPlaneNotCrossZSection) || (yzPlaneNotCrossYSection && xzPlaneAtOriginCrossYSection && xyPlaneAtOriginCrossZSection) || (yzPlaneAtImageBorderCrossXSection && xzPlaneAtImageBorderCrossYSection && xyPlaneNotCrossZSection) || (yzPlaneNotCrossYSection && xzPlaneAtImageBorderCrossYSection && xyPlaneAtImageBorderCrossZSection) || (yzPlaneNotCrossYSection && xzPlaneAtImageBorderCrossYSection && xyPlaneAtOriginCrossZSection) || (yzPlaneAtOriginCrossXSection && xzPlaneNotCrossYSection && xyPlaneAtOriginCrossZSection) || (yzPlaneNotCrossYSection && xzPlaneAtImageBorderCrossYSection && xyPlaneAtImageBorderCrossZSection)) || (yzPlaneAtImageBorderCrossXSection && xzPlaneNotCrossYSection && xyPlaneAtImageBorderCrossZSection) ) ) */ { // x-Axis if(xMin <= m_Origin[0] && xMax >= m_Origin[0]) { xMin = m_Origin[0]; }else if(xMin <= m_Size[0] * m_Spacing[0] && xMax >= m_Size[0] * m_Spacing[0]) { xMax = m_Size[0] * m_Spacing[0]; } // y-Axis if(yMin <= m_Origin[1] && yMax >= m_Origin[1]) { yMin = m_Origin[1]; }else if(yMin <= m_Size[1] * m_Spacing[1] && yMax >= m_Size[1] * m_Spacing[1]) { yMax = m_Size[1] * m_Spacing[1]; } // z-Axis if(zMin <= m_Origin[2] && zMax >= m_Origin[2]) { zMin = m_Origin[2]; }else if(zMin <= m_Size[2] * m_Spacing[2] && zMax >= m_Size[2] * m_Spacing[2]) { zMax = m_Size[2] * m_Spacing[2]; } m_meanValueTemp = m_meanValueTemp + this->MultiGaussianFunctionValueAtCuboid( xMin, xMax, yMin, yMax, zMin, zMax ); m_Volume = m_Volume + (xMax - xMin) * (yMax - yMin) * (zMax - zMin) ; } } } //----------------------------------------------------------------------------------------------------------------------- /* Start the octree recursion in eigth directions for the sphere with midpoint globalCoordinateMidpointSphere and, if uncommented, write this in a file, so that we can visualise it. */ template< class TOutputImage > void MultiGaussianImageSource< TOutputImage > ::GenerateCuboidSegmentationInSphere(PointType globalCoordinateMidpointSphere) { double cuboidRadius = m_Radius * 0.5; PointType newMidpoint; for(int i = -1; i < 2; i+=2) { for(int k = -1; k < 2; k+=2) { for(int j = -1; j < 2; j+=2) { newMidpoint[0] = globalCoordinateMidpointSphere[0] + static_cast(i) * cuboidRadius; newMidpoint[1] = globalCoordinateMidpointSphere[1] + static_cast(k) * cuboidRadius; newMidpoint[2] = globalCoordinateMidpointSphere[2] + static_cast(j) * cuboidRadius; CalculateEdgesInSphere( newMidpoint, globalCoordinateMidpointSphere, cuboidRadius, 0); } } } if(m_WriteMPS) { m_WriteMPS = 0; // uncomment to write an .mps file to visualise the midpoints // std::cout << "Wrote .xml to visualise the midpoints." << std::endl; // WriteXMLToTestTheCuboidInsideTheSphere(); } } //---------------------------------------------------------------------------------------------------------------------- /* This class allows by the method CalculateTheMidpointAndTheMeanValueWithOctree() to find a sphere with a specified radius that has a maximal mean value over all sphere with that radius with midpoint inside or at the boundary of the image. We approximaze the sphere with the octree recursiv method. CalculateTheMidpointAndTheMeanValueWithOctree works as follows: 1. The for-loops traverse the region of interest and assume the current point to be the wanted sphere midpoint. 2. Calculate the mean value for that sphere. 3. Compare with the until-now-found-maximum and take the bigger one. */ template< class TOutputImage > void MultiGaussianImageSource< TOutputImage > ::CalculateTheMidpointAndTheMeanValueWithOctree() { m_MeanValue = 0.0; double meanValueTemp; PointType midpoint; //typename MapContainerPoints::ElementIdentifier cuboidNumber = m_Midpoints.Size(); // SetNormalDistributionValues(); //double radius; //double xMin, xMax, yMin, yMax, zMin, zMax; m_WriteMPS = 1; PointType globalCoordinateMidpointSphere; IndexType index; OutputImageRegionType regionOfInterest; IndexType indexR; indexR.SetElement( 0, m_RegionOfInterestMin[0] ); indexR.SetElement( 1, m_RegionOfInterestMin[1] ); indexR.SetElement( 2, m_RegionOfInterestMin[2] ); regionOfInterest.SetIndex(indexR); SizeType sizeROI; sizeROI.SetElement( 0, m_RegionOfInterestMax[0] - m_RegionOfInterestMin[0] + 1); sizeROI.SetElement( 1, m_RegionOfInterestMax[1] - m_RegionOfInterestMin[1] + 1); sizeROI.SetElement( 2, m_RegionOfInterestMax[2] - m_RegionOfInterestMin[2] + 1); regionOfInterest.SetSize(sizeROI); typename TOutputImage::Pointer image = this->GetOutput(0); IteratorType regionOfInterestIterator(image, regionOfInterest); for(regionOfInterestIterator.GoToBegin(); !regionOfInterestIterator.IsAtEnd(); ++regionOfInterestIterator) { index = regionOfInterestIterator.GetIndex(); image->TransformIndexToPhysicalPoint(index, globalCoordinateMidpointSphere); m_Volume = 0.0; m_meanValueTemp = 0.0; this->GenerateCuboidSegmentationInSphere(globalCoordinateMidpointSphere); meanValueTemp = m_meanValueTemp / m_Volume; // std::cout << "index: " << index <<" meanvalue: " << meanValueTemp << std::endl; if(meanValueTemp > m_MeanValue) { m_GlobalCoordinate = globalCoordinateMidpointSphere; m_MeanValue = meanValueTemp; m_SphereMidpoint = index; } } } //---------------------------------------------------------------------------------------------------------------------- /* Check if a cuboid intersect the sphere boundary. Returns 2, if the cuboid is inside the sphere; returns 1, if the cuboid intersects the sphere boundary and 0, if the cuboid is out of the sphere. */ template< class TOutputImage > unsigned int MultiGaussianImageSource< TOutputImage > ::IntersectTheSphere( PointType globalCoordinateMidpointCuboid, PointType globalCoordinateMidpointSphere, double cuboidRadius) { unsigned int intersect = 1; PointType cuboidEdge; int count = 0; for(int i = -1; i < 2; i+=2) { for(int k = -1; k < 2; k+=2) { for(int j = -1; j < 2; j+=2) { cuboidEdge[0] = globalCoordinateMidpointCuboid[0] + static_cast(i) * cuboidRadius; cuboidEdge[1] = globalCoordinateMidpointCuboid[1] + static_cast(k) * cuboidRadius; cuboidEdge[2] = globalCoordinateMidpointCuboid[2] + static_cast(j) * cuboidRadius; if (globalCoordinateMidpointSphere.SquaredEuclideanDistanceTo(cuboidEdge) <= m_Radius * m_Radius) { ++count; } } } } if ( count == 0 ) { // cuboid not in the sphere intersect = 0; } if (count == 8 ) { // cuboid in the sphere intersect = 2; } return intersect; } //---------------------------------------------------------------------------------------------------------------------- template< class TOutputImage > double MultiGaussianImageSource< TOutputImage > ::MultiGaussianFunctionValueAtPoint(double x, double y, double z) { //this claculate the mean value in the voxel //integrate over the voxel with midpoint [x, y, z] double summand0, summand1, summand2/*, power*/, value = 0.0, factor; double xMin, xMax, yMin, yMax, zMin, zMax, mean; mean = 0.0; // the for-loop represent the sum of the gaussian function xMin = x - m_Spacing[0] / 2.0; xMax = x + m_Spacing[0] / 2.0; yMin = y - m_Spacing[1] / 2.0; yMax = y + m_Spacing[1] / 2.0; zMin = z - m_Spacing[2] / 2.0; zMax = z + m_Spacing[2] / 2.0; for( unsigned int n = 0; n < m_NumberOfGaussians; ++n ) { summand0 = FunctionPhi( (xMax - m_CenterX[n]) / m_SigmaX[n] ) - FunctionPhi( (xMin - m_CenterX[n]) / m_SigmaX[n] ); summand1 = FunctionPhi( (yMax - m_CenterY[n]) / m_SigmaY[n] ) - FunctionPhi( (yMin - m_CenterY[n]) / m_SigmaY[n] ); summand2 = FunctionPhi( (zMax - m_CenterZ[n]) / m_SigmaZ[n] ) - FunctionPhi( (zMin - m_CenterZ[n]) / m_SigmaZ[n] ); value = summand0 * summand1 * summand2; factor = ( m_SigmaX[n] * m_SigmaY[n] * m_SigmaZ[n] ) * pow( 2.0 * itk::Math::pi, 1.5 ); mean = mean + factor * value * m_Altitude[n]; } value = mean / (m_Spacing[0] * m_Spacing[1] * m_Spacing[2] ); /* //this calculate the value of the gaussian at the midpoint of the voxel: double summand0, summand1, summand2, power, value = 0.0; // the for-loop represent the sum of the gaussian function for(unsigned int n =0; n < m_NumberOfGaussians; ++n) { summand0 = ( x - m_CenterX[n] ) / m_SigmaX[n]; summand1 = ( y - m_CenterY[n] ) / m_SigmaY[n]; summand2 = ( z - m_CenterZ[n] ) / m_SigmaZ[n]; power = summand0 * summand0 + summand1 * summand1 + summand2 * summand2; value = value + m_Altitude[n] * pow(itk::Math::e, -0.5 * power); } */ // std::cout << "X: " << xMin << " " << x << " "<< xMax << " Y: "<< yMin << " " << y << " " << yMax << " Z: "<< zMin << " "<< z << " " << zMax << " value: " << value << std::endl; return value; } //---------------------------------------------------------------------------------------------------------------------- template< class TOutputImage > void MultiGaussianImageSource< TOutputImage > ::AddGaussian( VectorType x, VectorType y, VectorType z, VectorType sx, VectorType sy, VectorType sz, VectorType altitude) { for(unsigned int i = 0; i < x.size(); ++i) { m_CenterX.push_back( x[i] ); m_CenterY.push_back( y[i] ); m_CenterZ.push_back( z[i] ); m_SigmaX.push_back( sx[i] ); m_SigmaY.push_back( sy[i] ); m_SigmaZ.push_back( sz[i] ); m_Altitude.push_back(altitude[i]); } } //----------------------------------------------------------------------------------------------------------------------- template< typename TOutputImage > void MultiGaussianImageSource< TOutputImage > ::GenerateOutputInformation() { TOutputImage *output; IndexType index; index.Fill(0); output = this->GetOutput(0); typename TOutputImage::RegionType largestPossibleRegion; largestPossibleRegion.SetSize(this->m_Size); largestPossibleRegion.SetIndex(index); output->SetLargestPossibleRegion(largestPossibleRegion); output->SetSpacing(m_Spacing); output->SetOrigin(m_Origin); } //----------------------------------------------------------------------------------------------------------------------- template< typename TOutputImage > void MultiGaussianImageSource< TOutputImage > ::GenerateData() { itkDebugMacro(<< "Generating a image of scalars "); double valueReal; IndexType index; - typedef typename TOutputImage::PixelType scalarType; typename TOutputImage::Pointer image = this->GetOutput(0); image = this->GetOutput(0); image->SetBufferedRegion( image->GetRequestedRegion() ); image->Allocate(); IteratorType imageIt(image, image->GetLargestPossibleRegion()); PointType globalCoordinate; this->SetNormalDistributionValues(); for(imageIt.GoToBegin(); !imageIt.IsAtEnd(); ++imageIt) { valueReal = 0.0; index = imageIt.GetIndex(); image->TransformIndexToPhysicalPoint(imageIt.GetIndex(), globalCoordinate); valueReal = MultiGaussianFunctionValueAtPoint(globalCoordinate[0], globalCoordinate[1] ,globalCoordinate[2]); imageIt.Set(valueReal); } } //----------------------------------------------------------------------------------------------------------------------- /*This method is used to write a .mps file, so that we can visualize the midpoints of the approximated sphere as a scatterplot (for example with MITK Workbench). */ template< typename TOutputImage > void MultiGaussianImageSource< TOutputImage > ::WriteXMLToTestTheCuboidInsideTheSphere() { std::stringstream ss; int numberSummand = 1.0; //write an .mps test file itk::DOMNodeXMLWriter::Pointer xmlWriter; typedef itk::DOMNode::Pointer DOMNodeType; DOMNodeType domXML, domPointSetFile, domFileVersion, domPointSet, domPoint, domId, domX, domY, domZ; xmlWriter = itk::DOMNodeXMLWriter::New(); domXML = itk::DOMNode::New(); domXML->SetName("?xml"); domPointSetFile = itk::DOMNode::New(); domPointSetFile->SetName("point_set_file"); //domFileVersion = itk::DOMNode::New(); //domFileVersion->SetName("file_version"); domPointSet = itk::DOMNode::New(); domPointSet->SetName("point_set"); ss.str(""); ss << 1.0; domXML->SetAttribute("version", ss.str()); domXML->AddChildAtBegin(domPointSetFile); //domPointSetFile -> AddChildAtBegin(domFileVersion); domPointSetFile -> AddChildAtBegin(domPointSet); unsigned int cap = m_Midpoints.Size(); for(unsigned int iter = 0 ; iter < cap; ++iter) { domPoint = itk::DOMNode::New(); domPoint->SetName("point"); domX = itk::DOMNode::New(); domX->SetName("x"); domY = itk::DOMNode::New(); domY->SetName("y"); domZ = itk::DOMNode::New(); domZ->SetName("z"); domId = itk::DOMNode::New(); domId->SetName("id"); ss.str(""); ss << numberSummand; domId->AddTextChildAtBegin(ss.str()); domPoint -> AddChildAtEnd(domId); double scaleFactor = 10.0; PointType point = m_Midpoints.GetElement( numberSummand - 1 ); ss.str(""); ss << point[0] * scaleFactor; domX->AddTextChildAtBegin(ss.str()); domPoint -> AddChildAtEnd(domX); ss.str(""); ss << point[1] * scaleFactor; domY->AddTextChildAtBegin(ss.str()); domPoint -> AddChildAtEnd(domY); ss.str(""); ss << point[2] * scaleFactor; domZ->AddTextChildAtBegin(ss.str()); domPoint -> AddChildAtEnd(domZ); domPointSet -> AddChildAtEnd(domPoint); numberSummand += 1.0; } // .mps (Data) ss.str(""); ss << "C:/temp/CuboidsInTheSphere.mps"; std::string name = ss.str(); char * fileNamePointer = (char*) name.c_str(); xmlWriter->SetFileName( fileNamePointer); xmlWriter->SetInput( domXML ); xmlWriter->Update(); } //----------------------------------------------------------------------------------------------------------------------- template< typename TOutputImage > void MultiGaussianImageSource< TOutputImage > ::CalculateMaxAndMinInSphere() { IndexType index; typename MultiGaussianImageSource< TOutputImage >::OutputImagePixelType value; m_MaxValueInSphere = std::numeric_limits::min(); m_MinValueInSphere = std::numeric_limits::max(); int radInt, sizeRegion; OutputImageRegionType cuboidRegion; IndexType indexR; SizeType sizeR; int indexRegion, originAsIndex; for( unsigned int i = 0; i < 3; ++i ) { radInt = static_cast(m_Radius/m_Spacing[i]); indexRegion = m_SphereMidpoint[i] - radInt; originAsIndex = static_cast(m_Origin[i]/m_Spacing[i]); if( originAsIndex > indexRegion ) { indexR.SetElement(i, originAsIndex ); } else { indexR.SetElement(i, indexRegion ); } sizeRegion = 2 *radInt + 1; int sizeOutputImage = m_Size[i]; if( (indexR[i] + sizeRegion) > (originAsIndex + sizeOutputImage) ) { std::cout << "Not the entire sphere is in the image!" << std::endl; sizeR.SetElement(i, m_Size[i] - indexRegion ); } else { sizeR.SetElement(i, sizeRegion ); } } cuboidRegion.SetIndex(indexR); cuboidRegion.SetSize(sizeR); typename TOutputImage::Pointer image = this->GetOutput(0); IteratorType cuboidRegionOfInterestIterator(image, cuboidRegion); PointType globalCoordinate; for(cuboidRegionOfInterestIterator.GoToBegin(); !cuboidRegionOfInterestIterator.IsAtEnd(); ++cuboidRegionOfInterestIterator) { index = cuboidRegionOfInterestIterator.GetIndex(); if(IsInImage(index)) { image->TransformIndexToPhysicalPoint(cuboidRegionOfInterestIterator.GetIndex(), globalCoordinate); if( m_GlobalCoordinate.EuclideanDistanceTo(globalCoordinate) <= m_Radius ) { value = cuboidRegionOfInterestIterator.Get(); if(m_MaxValueInSphere < value) { m_MaxValueInSphere = value; m_MaxValueIndexInSphere = index; } if(m_MinValueInSphere > value) { m_MinValueInSphere = value; m_MinValueIndexInSphere = index; } } } } } //----------------------------------------------------------------------------------------------------------------------- template< typename TOutputImage > bool MultiGaussianImageSource< TOutputImage > ::IsInImage(IndexType index) { bool isInImage = true; int originAsIndex; for( unsigned int i = 0; i < 3; ++i ) { originAsIndex = static_cast(m_Origin[i]/m_Spacing[i]); int sizeOfOutputImage = m_Size[i]; if( index[i] < originAsIndex || index[i] > (originAsIndex + sizeOfOutputImage) ) return false; } return isInImage; } //----------------------------------------------------------------------------------------------------------------------- template< typename TOutputImage > void MultiGaussianImageSource< TOutputImage > ::SetNormalDistributionValues() { double temp[410] = { 0.5 , 0.50399 , 0.50798, 0.51197, 0.51595, 0.51994, 0.52392, 0.5279, 0.53188, 0.53586, 0.53983, 0.5438, 0.54776, 0.55172, 0.55567, 0.55962, 0.56356, 0.56749, 0.57142, 0.57535, 0.57926, 0.58317, 0.58706, 0.59095, 0.59483 , 0.59871, 0.60257, 0.60642, 0.61026, 0.61409, 0.61791, 0.62172, 0.62552, 0.6293, 0.63307, 0.63683, 0.64058, 0.64431, 0.64803, 0.65173, 0.65542, 0.6591, 0.66276, 0.6664, 0.67003, 0.67364, 0.67724, 0.68082, 0.68439, 0.68793, 0.69146, 0.69497, 0.69847, 0.70194, 0.7054, 0.70884, 0.71226, 0.71566, 0.71904, 0.7224, 0.72575, 0.72907, 0.73237, 0.73565, 0.73891, 0.74215, 0.74537, 0.74857, 0.75175, 0.7549, 0.75804, 0.76115, 0.76424, 0.7673, 0.77035, 0.77337, 0.77637, 0.77935, 0.7823, 0.78524, 0.78814, 0.79103, 0.79389, 0.79673, 0.79955, 0.80234, 0.80511, 0.80785, 0.81057, 0.81327, 0.81594, 0.81859, 0.82121, 0.82381, 0.82639, 0.82894, 0.83147, 0.83398, 0.83646, 0.83891, 0.84134, 0.84375, 0.84614, 0.84849, 0.85083, 0.85314, 0.85543, 0.85769, 0.85993, 0.86214, 0.86433, 0.8665, 0.86864, 0.87076, 0.87286, 0.87493, 0.87698, 0.879, 0.881, 0.88298, 0.88493, 0.88686, 0.88877, 0.89065, 0.89251, 0.89435, 0.89617, 0.89796, 0.89973, 0.90147, 0.9032, 0.9049, 0.90658, 0.90824, 0.90988, 0.91149, 0.91309, 0.91466, 0.91621, 0.91774, 0.91924, 0.92073, 0.9222, 0.92364, 0.92507, 0.92647, 0.92785, 0.92922, 0.93056, 0.93189, 0.93319, 0.93448, 0.93574, 0.93699, 0.93822, 0.93943, 0.94062, 0.94179, 0.94295, 0.94408, 0.9452, 0.9463, 0.94738, 0.94845, 0.9495, 0.95053, 0.95154, 0.95254, 0.95352, 0.95449, 0.95543, 0.95637, 0.95728, 0.95818, 0.95907, 0.95994, 0.9608, 0.96164, 0.96246, 0.96327, 0.96407, 0.96485, 0.96562, 0.96638, 0.96712, 0.96784, 0.96856, 0.96926, 0.96995, 0.97062, 0.97128, 0.97193, 0.97257, 0.9732, 0.97381, 0.97441, 0.975, 0.97558, 0.97615, 0.9767, 0.97725, 0.97778, 0.97831, 0.97882, 0.97932, 0.97982, 0.9803, 0.98077, 0.98124, 0.98169, 0.98214, 0.98257, 0.983, 0.98341, 0.98382, 0.98422, 0.98461, 0.985, 0.98537, 0.98574, 0.9861, 0.98645, 0.98679, 0.98713, 0.98745, 0.98778, 0.98809, 0.9884, 0.9887, 0.98899, 0.98928, 0.98956, 0.98983, 0.9901, 0.99036, 0.99061, 0.99086, 0.99111, 0.99134, 0.99158, 0.9918, 0.99202, 0.99224, 0.99245, 0.99266, 0.99286, 0.99305, 0.99324, 0.99343, 0.99361, 0.99379, 0.99396, 0.99413, 0.9943, 0.99446, 0.99461, 0.99477, 0.99492, 0.99506, 0.9952, 0.99534, 0.99547, 0.9956, 0.99573, 0.99585, 0.99598, 0.99609, 0.99621, 0.99632, 0.99643, 0.99653, 0.99664, 0.99674, 0.99683, 0.99693, 0.99702, 0.99711, 0.9972, 0.99728, 0.99736, 0.99744, 0.99752, 0.9976, 0.99767, 0.99774, 0.99781, 0.99788, 0.99795, 0.99801, 0.99807, 0.99813, 0.99819, 0.99825, 0.99831, 0.99836, 0.99841, 0.99846, 0.99851, 0.99856, 0.99861, 0.99865, 0.99869, 0.99874, 0.99878, 0.99882, 0.99886, 0.99889, 0.99893, 0.99896, 0.999, 0.99903, 0.99906, 0.9991, 0.99913, 0.99916, 0.99918, 0.99921, 0.99924, 0.99926, 0.99929, 0.99931, 0.99934, 0.99936, 0.99938, 0.9994, 0.99942, 0.99944, 0.99946, 0.99948, 0.9995, 0.99952, 0.99953, 0.99955, 0.99957, 0.99958, 0.9996, 0.99961, 0.99962, 0.99964, 0.99965, 0.99966, 0.99968, 0.99969, 0.9997, 0.99971, 0.99972, 0.99973, 0.99974, 0.99975, 0.99976, 0.99977, 0.99978, 0.99978, 0.99979, 0.9998, 0.99981, 0.99981, 0.99982, 0.99983, 0.99983, 0.99984, 0.99985, 0.99985, 0.99986, 0.99986, 0.99987, 0.99987, 0.99988, 0.99988, 0.99989, 0.99989, 0.9999, 0.9999, 0.9999, 0.99991, 0.99991, 0.99992, 0.99992, 0.99992, 0.99992, 0.99993, 0.99993, 0.99993, 0.99994, 0.99994, 0.99994, 0.99994, 0.99995, 0.99995, 0.99995, 0.99995, 0.99995, 0.99996, 0.99996, 0.99996, 0.99996, 0.99996, 0.99996, 0.99997, 0.99997, 0.99997, 0.99997, 0.99997, 0.99997, 0.99997, 0.99997, 0.99998, 0.99998, 0.99998, 0.99998 }; for(int i=0; i < 410; i++) { m_NormalDistValues[i] = temp[i]; } } } // end namespace itk #endif diff --git a/Modules/ImageStatistics/mitkExtendedLabelStatisticsImageFilter.hxx b/Modules/ImageStatistics/mitkExtendedLabelStatisticsImageFilter.hxx index 57e4dda221..403574bbda 100644 --- a/Modules/ImageStatistics/mitkExtendedLabelStatisticsImageFilter.hxx +++ b/Modules/ImageStatistics/mitkExtendedLabelStatisticsImageFilter.hxx @@ -1,765 +1,760 @@ /*=================================================================== 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 _mitkExtendedLabelStatisticsImageFilter_hxx #define _mitkExtendedLabelStatisticsImageFilter_hxx #include "mitkExtendedLabelStatisticsImageFilter.h" #include "itkImageRegionConstIteratorWithIndex.h" #include "itkImageRegionConstIterator.h" #include #include #include "mitkNumericConstants.h" #include "mitkLogMacros.h" #include namespace itk { template< class TInputImage, class TLabelImage > bool ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > ::GetMaskingNonEmpty() const { return m_MaskNonEmpty; } template< typename TInputImage, typename TLabelImage > void ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > ::SetHistogramParameters(const int numBins, RealType lowerBound, RealType upperBound) { m_NumBins[0] = numBins; m_LowerBound = lowerBound; m_UpperBound = upperBound; m_GlobalHistogramParametersSet = true; m_PreferGlobalHistogramParameters = true; this->Modified(); } template< typename TInputImage, typename TLabelImage > void ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > ::SetHistogramParametersForLabels(std::map numBins, std::map lowerBound, std::map upperBound) { m_LabelMin = lowerBound; m_LabelMax = upperBound; m_LabelNBins = numBins; m_LabelHistogramParametersSet = true; m_PreferGlobalHistogramParameters = false; this->Modified(); } template< class TInputImage, class TLabelImage > typename ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >::RealType ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > ::GetUniformity(LabelPixelType label) const { StatisticsMapConstIterator mapIt; mapIt = m_LabelStatistics.find(label); if ( mapIt == m_LabelStatistics.end() ) { // label does not exist, return a default value return NumericTraits< PixelType >::Zero; } else { return ( *mapIt ).second.m_Uniformity; } } template< class TInputImage, class TLabelImage > typename ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >::RealType ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > ::GetMedian(LabelPixelType label) const { StatisticsMapConstIterator mapIt; mapIt = m_LabelStatistics.find(label); if ( mapIt == m_LabelStatistics.end() ) { // label does not exist, return a default value return NumericTraits< PixelType >::Zero; } else { return ( *mapIt ).second.m_Median; } } template< class TInputImage, class TLabelImage > typename ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >::RealType ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > ::GetEntropy(LabelPixelType label) const { StatisticsMapConstIterator mapIt; mapIt = m_LabelStatistics.find(label); if ( mapIt == m_LabelStatistics.end() ) { // label does not exist, return a default value return NumericTraits< PixelType >::Zero; } else { return ( *mapIt ).second.m_Entropy; } } template< class TInputImage, class TLabelImage > typename ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >::RealType ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > ::GetUPP(LabelPixelType label) const { StatisticsMapConstIterator mapIt; mapIt = m_LabelStatistics.find(label); if ( mapIt == m_LabelStatistics.end() ) { // label does not exist, return a default value return NumericTraits< PixelType >::Zero; } else { return ( *mapIt ).second.m_UPP; } } template< class TInputImage, class TLabelImage > typename ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >::RealType ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > ::GetMPP(LabelPixelType label) const { StatisticsMapConstIterator mapIt; mapIt = m_LabelStatistics.find(label); if ( mapIt == m_LabelStatistics.end() ) { // label does not exist, return a default value return NumericTraits< PixelType >::Zero; } else { return ( *mapIt ).second.m_MPP; } } template< class TInputImage, class TLabelImage > typename ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >::RealType ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > ::GetKurtosis(LabelPixelType label) const { StatisticsMapConstIterator mapIt; mapIt = m_LabelStatistics.find(label); if ( mapIt == m_LabelStatistics.end() ) { // label does not exist, return a default value return NumericTraits< PixelType >::Zero; } else { return ( *mapIt ).second.m_Kurtosis; } } template< class TInputImage, class TLabelImage > typename ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >::RealType ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > ::GetSkewness(LabelPixelType label) const { StatisticsMapConstIterator mapIt; mapIt = m_LabelStatistics.find(label); if ( mapIt == m_LabelStatistics.end() ) { // label does not exist, return a default value return NumericTraits< PixelType >::Zero; } else { return ( *mapIt ).second.m_Skewness; } } template< typename TInputImage, typename TLabelImage > typename ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >::RealType ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > ::GetMinimum(LabelPixelType label) const { StatisticsMapConstIterator mapIt; mapIt = m_LabelStatistics.find(label); if ( mapIt == m_LabelStatistics.end() ) { // label does not exist, return a default value return NumericTraits< PixelType >::max(); } else { return ( *mapIt ).second.m_Minimum; } } template< typename TInputImage, typename TLabelImage > typename ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >::RealType ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > ::GetMaximum(LabelPixelType label) const { StatisticsMapConstIterator mapIt; mapIt = m_LabelStatistics.find(label); if ( mapIt == m_LabelStatistics.end() ) { // label does not exist, return a default value return NumericTraits< PixelType >::NonpositiveMin(); } else { return ( *mapIt ).second.m_Maximum; } } template< typename TInputImage, typename TLabelImage > typename ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >::RealType ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > ::GetMean(LabelPixelType label) const { StatisticsMapConstIterator mapIt; mapIt = m_LabelStatistics.find(label); if ( mapIt == m_LabelStatistics.end() ) { // label does not exist, return a default value return NumericTraits< PixelType >::ZeroValue(); } else { return ( *mapIt ).second.m_Mean; } } template< typename TInputImage, typename TLabelImage > typename ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >::RealType ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > ::GetSum(LabelPixelType label) const { StatisticsMapConstIterator mapIt; mapIt = m_LabelStatistics.find(label); if ( mapIt == m_LabelStatistics.end() ) { // label does not exist, return a default value return NumericTraits< PixelType >::ZeroValue(); } else { return ( *mapIt ).second.m_Sum; } } template< typename TInputImage, typename TLabelImage > typename ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >::RealType ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > ::GetSigma(LabelPixelType label) const { StatisticsMapConstIterator mapIt; mapIt = m_LabelStatistics.find(label); if ( mapIt == m_LabelStatistics.end() ) { // label does not exist, return a default value return NumericTraits< PixelType >::ZeroValue(); } else { return ( *mapIt ).second.m_Sigma; } } template< typename TInputImage, typename TLabelImage > typename ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >::RealType ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > ::GetVariance(LabelPixelType label) const { StatisticsMapConstIterator mapIt; mapIt = m_LabelStatistics.find(label); if ( mapIt == m_LabelStatistics.end() ) { // label does not exist, return a default value return NumericTraits< PixelType >::ZeroValue(); } else { return ( *mapIt ).second.m_Variance; } } template< typename TInputImage, typename TLabelImage > typename ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >::BoundingBoxType ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > ::GetBoundingBox(LabelPixelType label) const { StatisticsMapConstIterator mapIt; mapIt = m_LabelStatistics.find(label); if ( mapIt == m_LabelStatistics.end() ) { typename Superclass::BoundingBoxType emptyBox; // label does not exist, return a default value return emptyBox; } else { return ( *mapIt ).second.m_BoundingBox; } } template< typename TInputImage, typename TLabelImage > typename ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >::RegionType ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > ::GetRegion(LabelPixelType label) const { StatisticsMapConstIterator mapIt; mapIt = m_LabelStatistics.find(label); if ( mapIt == m_LabelStatistics.end() ) { typename Superclass::RegionType emptyRegion; // label does not exist, return a default value return emptyRegion; } else { typename Superclass::BoundingBoxType bbox = this->GetBoundingBox(label); typename Superclass::IndexType index; typename Superclass::SizeType size; unsigned int dimension = bbox.size() / 2; for ( unsigned int i = 0; i < dimension; i++ ) { index[i] = bbox[2 * i]; size[i] = bbox[2 * i + 1] - bbox[2 * i] + 1; } typename Superclass::RegionType region; region.SetSize(size); region.SetIndex(index); return region; } } template< typename TInputImage, typename TLabelImage > typename ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >::MapSizeType ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > ::GetCount(LabelPixelType label) const { StatisticsMapConstIterator mapIt; mapIt = m_LabelStatistics.find(label); if ( mapIt == m_LabelStatistics.end() ) { // label does not exist, return a default value return 0; } else { return ( *mapIt ).second.m_Count; } } template< typename TInputImage, typename TLabelImage > typename ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >::HistogramType::Pointer ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > ::GetHistogram(LabelPixelType label) const { StatisticsMapConstIterator mapIt; mapIt = m_LabelStatistics.find(label); if ( mapIt == m_LabelStatistics.end() ) { // label does not exist, return a default value return ITK_NULLPTR; } else { // this will be zero if histograms have not been enabled return ( *mapIt ).second.m_Histogram; } } template< typename TInputImage, typename TLabelImage > void ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > ::BeforeThreadedGenerateData() { ThreadIdType numberOfThreads = this->GetNumberOfThreads(); // Resize the thread temporaries m_LabelStatisticsPerThread.resize(numberOfThreads); // Initialize the temporaries for ( ThreadIdType i = 0; i < numberOfThreads; ++i ) { m_LabelStatisticsPerThread[i].clear(); } // Initialize the final map m_LabelStatistics.clear(); } template< typename TInputImage, typename TLabelImage > std::list ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > ::GetRelevantLabels() const { std::list< int> relevantLabels; for (int i = 0; i < 4096; ++i ) { if ( this->HasLabel( i ) ) { relevantLabels.push_back( i ); } } return relevantLabels; } template< typename TInputImage, typename TLabelImage > void ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage > ::ThreadedGenerateData(const typename TInputImage::RegionType & outputRegionForThread, ThreadIdType threadId) { typename HistogramType::IndexType histogramIndex(1); typename HistogramType::MeasurementVectorType histogramMeasurement(1); const SizeValueType size0 = outputRegionForThread.GetSize(0); if( size0 == 0) { return; } ImageLinearConstIteratorWithIndex< TInputImage > it (this->GetInput(), outputRegionForThread); ImageScanlineConstIterator< TLabelImage > labelIt (this->GetLabelInput(), outputRegionForThread); StatisticsMapIterator mapIt; // support progress methods/callbacks const size_t numberOfLinesToProcess = outputRegionForThread.GetNumberOfPixels() / size0; ProgressReporter progress( this, threadId, numberOfLinesToProcess ); typedef typename MapType::value_type MapValueType; // do the work while ( !it.IsAtEnd() ) { while ( !it.IsAtEndOfLine() ) { const RealType & value = static_cast< RealType >( it.Get() ); const LabelPixelType & label = labelIt.Get(); // is the label already in this thread? mapIt = m_LabelStatisticsPerThread[threadId].find(label); if ( mapIt == m_LabelStatisticsPerThread[threadId].end() ) { // if global histogram parameters are set and preferred then use them if ( m_PreferGlobalHistogramParameters && m_GlobalHistogramParametersSet ) { mapIt = m_LabelStatisticsPerThread[threadId].insert( MapValueType( label, LabelStatistics(m_NumBins[0], m_LowerBound, m_UpperBound) ) ).first; } // if we have label histogram parameters then use them. If we encounter a label that has no parameters then use global settings if available else if(!m_PreferGlobalHistogramParameters && m_LabelHistogramParametersSet) { typename std::map::iterator lbIt, ubIt; typename std::map::iterator nbIt; lbIt = m_LabelMin.find(label); ubIt = m_LabelMax.find(label); nbIt = m_LabelNBins.find(label); // if any of the parameters is lacking for the current label but global histogram params are available, use the global parameters if ((lbIt == m_LabelMin.end() || ubIt == m_LabelMax.end() || nbIt == m_LabelNBins.end()) && m_GlobalHistogramParametersSet) { mapIt = m_LabelStatisticsPerThread[threadId].insert( MapValueType( label, LabelStatistics(m_NumBins[0], m_LowerBound, m_UpperBound) ) ).first; } // if any of the parameters is lacking for the current label and global histogram params are not available, dont use histograms for this label else if ((lbIt == m_LabelMin.end() || ubIt == m_LabelMax.end() || nbIt == m_LabelNBins.end()) && !m_GlobalHistogramParametersSet) { mapIt = m_LabelStatisticsPerThread[threadId].insert( MapValueType( label, LabelStatistics() ) ).first; } // label histogram parameters are available, use them! else { PixelType lowerBound, upperBound; unsigned int nBins; lowerBound = (*lbIt).second; upperBound = (*ubIt).second; nBins = (*nbIt).second; mapIt = m_LabelStatisticsPerThread[threadId].insert( MapValueType( label, LabelStatistics(nBins, lowerBound, upperBound) ) ).first; } } // neither global nor label specific histogram parameters are set -> don't use histograms else { mapIt = m_LabelStatisticsPerThread[threadId].insert( MapValueType( label, LabelStatistics() ) ).first; } } typename MapType::mapped_type &labelStats = ( *mapIt ).second; // update the values for this label and this thread if ( value < labelStats.m_Minimum ) { labelStats.m_Minimum = value; } if ( value > labelStats.m_Maximum ) { labelStats.m_Maximum = value; } // bounding box is min,max pairs for ( unsigned int i = 0; i < ( 2 * TInputImage::ImageDimension ); i += 2 ) { const typename TInputImage::IndexType & index = it.GetIndex(); if ( labelStats.m_BoundingBox[i] > index[i / 2] ) { labelStats.m_BoundingBox[i] = index[i / 2]; } if ( labelStats.m_BoundingBox[i + 1] < index[i / 2] ) { labelStats.m_BoundingBox[i + 1] = index[i / 2]; } } labelStats.m_Sum += value; labelStats.m_SumOfSquares += ( value * value ); labelStats.m_Count++; labelStats.m_SumOfCubes += std::pow(value, 3.); labelStats.m_SumOfQuadruples += std::pow(value, 4.); if (value > 0) { labelStats.m_PositivePixelCount++; labelStats.m_SumOfPositivePixels += value; } // if enabled, update the histogram for this label if ( labelStats.m_Histogram.IsNotNull() ) - { + { histogramMeasurement[0] = value; labelStats.m_Histogram->GetIndex(histogramMeasurement, histogramIndex); labelStats.m_Histogram->IncreaseFrequencyOfIndex(histogramIndex, 1); - } - else - { - int x = 0; } - ++labelIt; ++it; } labelIt.NextLine(); it.NextLine(); progress.CompletedPixel(); } } template< class TInputImage, class TLabelImage > void ExtendedLabelStatisticsImageFilter< TInputImage, TLabelImage >:: AfterThreadedGenerateData() { StatisticsMapIterator mapIt; StatisticsMapConstIterator threadIt; ThreadIdType i; ThreadIdType numberOfThreads = this->GetNumberOfThreads(); // Run through the map for each thread and accumulate the count, // sum, and sumofsquares for ( i = 0; i < numberOfThreads; i++ ) { // iterate over the map for this thread for ( threadIt = m_LabelStatisticsPerThread[i].begin(); threadIt != m_LabelStatisticsPerThread[i].end(); ++threadIt ) { // does this label exist in the cumulative structure yet? mapIt = m_LabelStatistics.find( ( *threadIt ).first ); if ( mapIt == m_LabelStatistics.end() ) { // create a new entry typedef typename MapType::value_type MapValueType; if ( m_GlobalHistogramParametersSet || m_LabelHistogramParametersSet ) { // mapIt = m_LabelStatistics.insert( MapValueType( ( *threadIt ).first, // LabelStatistics(m_NumBins[0], m_LowerBound, // m_UpperBound) ) ).first; mapIt = m_LabelStatistics.insert( MapValueType( *threadIt ) ).first; continue; } else { mapIt = m_LabelStatistics.insert( MapValueType( ( *threadIt ).first, LabelStatistics() ) ).first; } } typename MapType::mapped_type &labelStats = ( *mapIt ).second; // accumulate the information from this thread labelStats.m_Count += ( *threadIt ).second.m_Count; labelStats.m_Sum += ( *threadIt ).second.m_Sum; labelStats.m_SumOfSquares += ( *threadIt ).second.m_SumOfSquares; labelStats.m_SumOfPositivePixels += ( *threadIt ).second.m_SumOfPositivePixels; labelStats.m_PositivePixelCount += ( *threadIt ).second.m_PositivePixelCount; labelStats.m_SumOfCubes += ( *threadIt ).second.m_SumOfCubes; labelStats.m_SumOfQuadruples += ( *threadIt ).second.m_SumOfQuadruples; if ( labelStats.m_Minimum > ( *threadIt ).second.m_Minimum ) { labelStats.m_Minimum = ( *threadIt ).second.m_Minimum; } if ( labelStats.m_Maximum < ( *threadIt ).second.m_Maximum ) { labelStats.m_Maximum = ( *threadIt ).second.m_Maximum; } //bounding box is min,max pairs int dimension = labelStats.m_BoundingBox.size() / 2; for ( int ii = 0; ii < ( dimension * 2 ); ii += 2 ) { if ( labelStats.m_BoundingBox[ii] > ( *threadIt ).second.m_BoundingBox[ii] ) { labelStats.m_BoundingBox[ii] = ( *threadIt ).second.m_BoundingBox[ii]; } if ( labelStats.m_BoundingBox[ii + 1] < ( *threadIt ).second.m_BoundingBox[ii + 1] ) { labelStats.m_BoundingBox[ii + 1] = ( *threadIt ).second.m_BoundingBox[ii + 1]; } } // if enabled, update the histogram for this label if ( m_GlobalHistogramParametersSet || m_LabelHistogramParametersSet ) { typename HistogramType::IndexType index; index.SetSize(1); for ( unsigned int bin = 0; bin < labelStats.m_Histogram->Size(); bin++ ) { index[0] = bin; labelStats.m_Histogram->IncreaseFrequency( bin, ( *threadIt ).second.m_Histogram->GetFrequency(bin) ); } } } // end of thread map iterator loop } // end of thread loop // compute the remainder of the statistics for ( mapIt = m_LabelStatistics.begin(); mapIt != m_LabelStatistics.end(); ++mapIt ) { typename MapType::mapped_type &labelStats = ( *mapIt ).second; // mean labelStats.m_Mean = labelStats.m_Sum / static_cast< RealType >( labelStats.m_Count ); // MPP labelStats.m_MPP = labelStats.m_SumOfPositivePixels / static_cast< RealType >( labelStats.m_PositivePixelCount ); // variance if ( labelStats.m_Count > 0 ) { // unbiased estimate of variance LabelStatistics & ls = mapIt->second; const RealType sumSquared = ls.m_Sum * ls.m_Sum; const RealType count = static_cast< RealType >( ls.m_Count ); ls.m_Variance = ( ls.m_SumOfSquares - sumSquared / count ) / ( count ); RealType secondMoment = ls.m_SumOfSquares / count; RealType thirdMoment = ls.m_SumOfCubes / count; RealType fourthMoment = ls.m_SumOfQuadruples / count; ls.m_Skewness = (thirdMoment - 3. * secondMoment * ls.m_Mean + 2. * std::pow(ls.m_Mean, 3.)) / std::pow(secondMoment - std::pow(ls.m_Mean, 2.), 1.5); // see http://www.boost.org/doc/libs/1_51_0/doc/html/boost/accumulators/impl/skewness_impl.html ls.m_Kurtosis = (fourthMoment - 4. * thirdMoment * ls.m_Mean + 6. * secondMoment * std::pow(ls.m_Mean, 2.) - 3. * std::pow(ls.m_Mean, 4.)) / std::pow(secondMoment - std::pow(ls.m_Mean, 2.), 2.); // see http://www.boost.org/doc/libs/1_51_0/doc/html/boost/accumulators/impl/kurtosis_impl.html, dropped -3 } else { labelStats.m_Variance = NumericTraits< RealType >::ZeroValue(); labelStats.m_Skewness = NumericTraits< RealType >::ZeroValue(); labelStats.m_Kurtosis = NumericTraits< RealType >::ZeroValue(); } // sigma labelStats.m_Sigma = std::sqrt( labelStats.m_Variance ); // histogram statistics if (labelStats.m_Histogram.IsNotNull()) { mitk::HistogramStatisticsCalculator histStatCalc; histStatCalc.SetHistogram(labelStats.m_Histogram); histStatCalc.CalculateStatistics(); labelStats.m_Median = histStatCalc.GetMedian(); labelStats.m_Entropy = histStatCalc.GetEntropy(); labelStats.m_Uniformity = histStatCalc.GetUniformity(); labelStats.m_UPP = histStatCalc.GetUPP(); } } { //Now update the cached vector of valid labels. m_ValidLabelValues.resize(0); m_ValidLabelValues.reserve(m_LabelStatistics.size()); for ( mapIt = m_LabelStatistics.begin(); mapIt != m_LabelStatistics.end(); ++mapIt ) { m_ValidLabelValues.push_back(mapIt->first); } } } } // end namespace itk #endif diff --git a/Modules/ImageStatistics/mitkExtendedStatisticsImageFilter.h b/Modules/ImageStatistics/mitkExtendedStatisticsImageFilter.h index fc4ea75878..75977cf508 100644 --- a/Modules/ImageStatistics/mitkExtendedStatisticsImageFilter.h +++ b/Modules/ImageStatistics/mitkExtendedStatisticsImageFilter.h @@ -1,210 +1,211 @@ /*=================================================================== 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 __mitkExtendedStatisticsImageFilter #define __mitkExtendedStatisticsImageFilter #include "itkStatisticsImageFilter.h" #include #include namespace itk { /** * \class ExtendedStatisticsImageFilter * \brief Extension of the itkStatisticsImageFilter that also calculates the Skewness and Kurtosis. * * This class inherits from the itkStatisticsImageFilter and * uses its results for the calculation of other additional coefficients: * the Skewness and Kurtosis. * * As the StatisticsImageFilter stores the statistics in the outputs 1 to 6 by the * StatisticsImageFilter, the skewness, kurtosis, MPP, UPP, Uniformity, Entropy and Median are stored in the outputs * 7 to 14 by this filter. */ template< class TInputImage > class ExtendedStatisticsImageFilter : public StatisticsImageFilter< TInputImage > { public: /** Standard Self typedef */ typedef ExtendedStatisticsImageFilter Self; typedef StatisticsImageFilter< TInputImage > Superclass; typedef SmartPointer< Self > Pointer; typedef SmartPointer< const Self > ConstPointer; typedef typename Superclass::RealType RealType; typedef typename Superclass::RealObjectType RealObjectType; typedef typename Superclass::PixelType PixelType; /** Histogram-related typedefs */ typedef itk::Statistics::Histogram< RealType > HistogramType; typedef typename HistogramType::Pointer HistogramPointer; itkFactorylessNewMacro( Self ); itkCloneMacro( Self ); itkTypeMacro( ExtendedStatisticsImageFilter, StatisticsImageFilter ); /** * \brief Return the computed Skewness. */ double GetSkewness() const { return this->GetSkewnessOutput()->Get(); } /** * \brief Return the computed Median */ double GetMedian() const { return this->GetMedianOutput()->Get(); } /** * \brief Return the computed Kurtosis. */ double GetKurtosis() const { return this->GetKurtosisOutput()->Get(); } /* \brief Return the computed MPP. */ double GetMPP() const { return this->GetMPPOutput()->Get(); } /** * \brief Return the computed Uniformity. */ double GetUniformity() const { return this->GetUniformityOutput()->Get(); } /** *\brief Return the computed Entropy. */ double GetEntropy() const { return this->GetEntropyOutput()->Get(); } /** * \brief Return the computed UPP. */ double GetUPP() const { return this->GetUPPOutput()->Get(); } /** * \brief Return the computed Histogram. */ const typename HistogramType::Pointer GetHistogram() { if (m_HistogramCalculated) { return m_Histogram; } else { return nullptr; } } /** specify Histogram parameters */ void SetHistogramParameters(const int numBins, RealType lowerBound, RealType upperBound); protected: ExtendedStatisticsImageFilter(); virtual ~ExtendedStatisticsImageFilter(){}; void BeforeThreadedGenerateData(); /** Multi-thread version GenerateData. */ void ThreadedGenerateData(const typename StatisticsImageFilter::RegionType & outputRegionForThread, ThreadIdType threadId); /** * brief Calls AfterThreadedGenerateData() of the superclass and the main methods */ void AfterThreadedGenerateData(); RealObjectType* GetSkewnessOutput(); const RealObjectType* GetSkewnessOutput() const; RealObjectType* GetKurtosisOutput(); const RealObjectType* GetKurtosisOutput() const; RealObjectType* GetMPPOutput(); const RealObjectType* GetMPPOutput() const; RealObjectType* GetEntropyOutput(); const RealObjectType* GetEntropyOutput() const; RealObjectType* GetUniformityOutput(); const RealObjectType* GetUniformityOutput() const; RealObjectType* GetUPPOutput(); const RealObjectType* GetUPPOutput() const; RealObjectType* GetMedianOutput(); const RealObjectType* GetMedianOutput() const; - virtual DataObject::Pointer MakeOutput( ProcessObject::DataObjectPointerArraySizeType idx ); + using Superclass::MakeOutput; + virtual DataObject::Pointer MakeOutput( ProcessObject::DataObjectPointerArraySizeType idx ) override; private: Array< RealType > m_ThreadSum; Array< RealType > m_SumOfSquares; Array< RealType > m_SumOfCubes; Array< RealType > m_SumOfQuadruples; Array< SizeValueType > m_Count; Array< SizeValueType > m_PositivePixelCount; Array< RealType > m_ThreadSumOfPositivePixels; Array< PixelType > m_ThreadMin; Array< PixelType > m_ThreadMax; std::vector< HistogramPointer > m_HistogramPerThread; HistogramPointer m_Histogram; bool m_UseHistogram; bool m_HistogramCalculated; RealType m_LowerBound, m_UpperBound; int m_NumBins; }; // end of class } // end namespace itk #ifndef ITK_MANUAL_INSTANTIATION #include "mitkExtendedStatisticsImageFilter.hxx" #endif #endif diff --git a/Modules/ImageStatistics/mitkExtendedStatisticsImageFilter.hxx b/Modules/ImageStatistics/mitkExtendedStatisticsImageFilter.hxx index 9413e33c0c..9ea49c069b 100644 --- a/Modules/ImageStatistics/mitkExtendedStatisticsImageFilter.hxx +++ b/Modules/ImageStatistics/mitkExtendedStatisticsImageFilter.hxx @@ -1,477 +1,477 @@ /*=================================================================== 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 __mitkExtendedStatisticsImageFilter_hxx #define __mitkExtendedStatisticsImageFilter_hxx #include "mitkExtendedStatisticsImageFilter.h" #include namespace itk { template< class TInputImage > ExtendedStatisticsImageFilter< TInputImage >::ExtendedStatisticsImageFilter() : StatisticsImageFilter< TInputImage >(), m_ThreadSum(1), m_SumOfSquares(1), m_SumOfCubes(1), m_SumOfQuadruples(1), m_Count(1), - m_ThreadMin(1), - m_ThreadMax(1), m_PositivePixelCount(1), m_ThreadSumOfPositivePixels(1), + m_ThreadMin(1), + m_ThreadMax(1), m_UseHistogram(false), m_HistogramCalculated(false) { /* * add the Skewness,Kurtosis,Entropy,Uniformity,MPP, UPP, Median to the other statistical calculated Values * of the mitkStatisticsImageFilter as the 7th to the 13th Output */ for ( int i = 7; i < 14; ++i ) { typename RealObjectType::Pointer output = static_cast< RealObjectType * >( this->MakeOutput(i).GetPointer() ); this->ProcessObject::SetNthOutput( i, output.GetPointer() ); } this->GetSkewnessOutput()->Set( 0.0 ); this->GetKurtosisOutput()->Set( 0.0 ); this->GetEntropyOutput()->Set( -1.0 ); this->GetUniformityOutput()->Set( 0.0 ); this->GetUPPOutput()->Set( 0.0 ); this->GetMPPOutput()->Set( 0.0 ); this->GetMedianOutput()->Set( 0.0 ); m_Histogram = ITK_NULLPTR; m_HistogramPerThread.resize(0); } template< class TInputImage > DataObject::Pointer ExtendedStatisticsImageFilter< TInputImage >::MakeOutput( ProcessObject::DataObjectPointerArraySizeType output) { switch ( output ) { case 7: case 8: case 9: case 10: case 11: case 12: { return RealObjectType::New().GetPointer(); break; } case 13: { return RealObjectType::New().GetPointer(); break; } default: { // might as well make an image return Superclass::MakeOutput( output ); break; } } } template< class TInputImage > typename ExtendedStatisticsImageFilter< TInputImage >::RealObjectType * ExtendedStatisticsImageFilter< TInputImage > ::GetSkewnessOutput() { return static_cast< RealObjectType * >( this->ProcessObject::GetOutput(7) ); } template< class TInputImage > const typename ExtendedStatisticsImageFilter< TInputImage >::RealObjectType * ExtendedStatisticsImageFilter< TInputImage > ::GetSkewnessOutput() const { return static_cast< const RealObjectType * >( this->ProcessObject::GetOutput(7) ); } template< class TInputImage > typename ExtendedStatisticsImageFilter< TInputImage >::RealObjectType * ExtendedStatisticsImageFilter< TInputImage > ::GetKurtosisOutput() { return static_cast< RealObjectType * >( this->ProcessObject::GetOutput(8) ); } template< class TInputImage > const typename ExtendedStatisticsImageFilter< TInputImage >::RealObjectType * ExtendedStatisticsImageFilter< TInputImage > ::GetKurtosisOutput() const { return static_cast< const RealObjectType * >( this->ProcessObject::GetOutput(8) ); } template< class TInputImage > typename ExtendedStatisticsImageFilter< TInputImage >::RealObjectType * ExtendedStatisticsImageFilter< TInputImage > ::GetUniformityOutput() { return static_cast< RealObjectType * >( this->ProcessObject::GetOutput(9) ); } template< class TInputImage > const typename ExtendedStatisticsImageFilter< TInputImage >::RealObjectType * ExtendedStatisticsImageFilter< TInputImage > ::GetUniformityOutput() const { return static_cast< const RealObjectType * >( this->ProcessObject::GetOutput(9) ); } template< class TInputImage > typename ExtendedStatisticsImageFilter< TInputImage >::RealObjectType * ExtendedStatisticsImageFilter< TInputImage > ::GetEntropyOutput() { return static_cast< RealObjectType * >( this->ProcessObject::GetOutput(10) ); } template< class TInputImage > const typename ExtendedStatisticsImageFilter< TInputImage >::RealObjectType * ExtendedStatisticsImageFilter< TInputImage > ::GetEntropyOutput() const { return static_cast< const RealObjectType * >( this->ProcessObject::GetOutput(10) ); } template< class TInputImage > typename ExtendedStatisticsImageFilter< TInputImage >::RealObjectType * ExtendedStatisticsImageFilter< TInputImage > ::GetUPPOutput() { return static_cast< RealObjectType * >( this->ProcessObject::GetOutput(11) ); } template< class TInputImage > const typename ExtendedStatisticsImageFilter< TInputImage >::RealObjectType * ExtendedStatisticsImageFilter< TInputImage > ::GetUPPOutput() const { return static_cast< const RealObjectType * >( this->ProcessObject::GetOutput(11) ); } template< class TInputImage > typename ExtendedStatisticsImageFilter< TInputImage >::RealObjectType * ExtendedStatisticsImageFilter< TInputImage > ::GetMPPOutput() { return static_cast< RealObjectType * >( this->ProcessObject::GetOutput(12) ); } template< class TInputImage > const typename ExtendedStatisticsImageFilter< TInputImage >::RealObjectType * ExtendedStatisticsImageFilter< TInputImage > ::GetMPPOutput() const { return static_cast< const RealObjectType * >( this->ProcessObject::GetOutput(12) ); } template< class TInputImage > typename ExtendedStatisticsImageFilter< TInputImage >::RealObjectType * ExtendedStatisticsImageFilter< TInputImage > ::GetMedianOutput() { return static_cast< RealObjectType * >( this->ProcessObject::GetOutput(13) ); } template< class TInputImage > const typename ExtendedStatisticsImageFilter< TInputImage >::RealObjectType * ExtendedStatisticsImageFilter< TInputImage > ::GetMedianOutput() const { return static_cast< const RealObjectType * >( this->ProcessObject::GetOutput(13) ); } template< typename TInputImage > void ExtendedStatisticsImageFilter< TInputImage > ::BeforeThreadedGenerateData() { ThreadIdType numberOfThreads = this->GetNumberOfThreads(); if (m_UseHistogram) { // Histogram m_Histogram = HistogramType::New(); typename HistogramType::SizeType hsize; typename HistogramType::MeasurementVectorType lb; typename HistogramType::MeasurementVectorType ub; hsize.SetSize(1); lb.SetSize(1); ub.SetSize(1); m_Histogram->SetMeasurementVectorSize(1); hsize[0] = m_NumBins; lb[0] = m_LowerBound; ub[0] = m_UpperBound; m_Histogram->Initialize(hsize, lb, ub); m_HistogramPerThread.resize(numberOfThreads); for (unsigned int i = 0; i < numberOfThreads; i++) { typename HistogramType::Pointer hist = HistogramType::New(); typename HistogramType::SizeType hsize; typename HistogramType::MeasurementVectorType lb; typename HistogramType::MeasurementVectorType ub; hsize.SetSize(1); lb.SetSize(1); ub.SetSize(1); hist->SetMeasurementVectorSize(1); hsize[0] = m_NumBins; lb[0] = m_LowerBound; ub[0] = m_UpperBound; hist->Initialize(hsize, lb, ub); m_HistogramPerThread[i] = hist; } } // Resize the thread temporaries m_Count.SetSize(numberOfThreads); m_SumOfSquares.SetSize(numberOfThreads); m_SumOfCubes.SetSize(numberOfThreads); m_SumOfQuadruples.SetSize(numberOfThreads); m_ThreadSum.SetSize(numberOfThreads); m_ThreadMin.SetSize(numberOfThreads); m_ThreadMax.SetSize(numberOfThreads); m_PositivePixelCount.SetSize(numberOfThreads); m_ThreadSumOfPositivePixels.SetSize(numberOfThreads); // Initialize the temporaries m_Count.Fill(NumericTraits< SizeValueType >::Zero); m_ThreadSum.Fill(NumericTraits< RealType >::Zero); m_SumOfSquares.Fill(NumericTraits< RealType >::Zero); m_SumOfCubes.Fill(NumericTraits< RealType >::Zero); m_SumOfQuadruples.Fill(NumericTraits< RealType >::Zero); m_ThreadMin.Fill( NumericTraits< PixelType >::max() ); m_ThreadMax.Fill( NumericTraits< PixelType >::NonpositiveMin() ); m_PositivePixelCount.Fill(NumericTraits< SizeValueType >::Zero); m_ThreadSumOfPositivePixels.Fill(NumericTraits< SizeValueType >::Zero); } template< typename TInputImage > void ExtendedStatisticsImageFilter< TInputImage > ::ThreadedGenerateData(const typename StatisticsImageFilter::RegionType & outputRegionForThread, ThreadIdType threadId) { const SizeValueType size0 = outputRegionForThread.GetSize(0); if( size0 == 0) { return; } RealType realValue; PixelType value; typename HistogramType::IndexType histogramIndex(1); typename HistogramType::MeasurementVectorType histogramMeasurement(1); RealType sum = NumericTraits< RealType >::ZeroValue(); RealType sumOfPositivePixels = NumericTraits< RealType >::ZeroValue(); RealType sumOfSquares = NumericTraits< RealType >::ZeroValue(); RealType sumOfCubes = NumericTraits< RealType >::ZeroValue(); RealType sumOfQuadruples = NumericTraits< RealType >::ZeroValue(); SizeValueType count = NumericTraits< SizeValueType >::ZeroValue(); SizeValueType countOfPositivePixels = NumericTraits< SizeValueType >::ZeroValue(); PixelType min = NumericTraits< PixelType >::max(); PixelType max = NumericTraits< PixelType >::NonpositiveMin(); ImageScanlineConstIterator< TInputImage > it (this->GetInput(), outputRegionForThread); // support progress methods/callbacks const size_t numberOfLinesToProcess = outputRegionForThread.GetNumberOfPixels() / size0; ProgressReporter progress( this, threadId, numberOfLinesToProcess ); // do the work while ( !it.IsAtEnd() ) { while ( !it.IsAtEndOfLine() ) { value = it.Get(); realValue = static_cast< RealType >( value ); if (m_UseHistogram) { histogramMeasurement[0] = value; m_HistogramPerThread[threadId]->GetIndex(histogramMeasurement, histogramIndex); m_HistogramPerThread[threadId]->IncreaseFrequencyOfIndex(histogramIndex, 1); } if ( value < min ) { min = value; } if ( value > max ) { max = value; } if (value > 0) { sumOfPositivePixels += realValue; ++countOfPositivePixels; } sum += realValue; sumOfSquares += ( realValue * realValue ); sumOfCubes += std::pow(realValue, 3.); sumOfQuadruples += std::pow(realValue, 4.); ++count; ++it; } it.NextLine(); progress.CompletedPixel(); } m_ThreadSum[threadId] = sum; m_SumOfSquares[threadId] = sumOfSquares; m_Count[threadId] = count; m_ThreadMin[threadId] = min; m_ThreadMax[threadId] = max; m_ThreadSumOfPositivePixels[threadId] = sumOfPositivePixels; m_PositivePixelCount[threadId] = countOfPositivePixels; m_SumOfCubes[threadId] = sumOfCubes; m_SumOfQuadruples[threadId] = sumOfQuadruples; } template< class TInputImage > void ExtendedStatisticsImageFilter< TInputImage > ::AfterThreadedGenerateData() { ThreadIdType i; SizeValueType count = 0; SizeValueType countOfPositivePixels = 0; RealType sumOfSquares = 0; RealType sumOfCubes = 0; RealType sumOfQuadruples = 0; RealType sumOfPositivePixels = 0; ThreadIdType numberOfThreads = this->GetNumberOfThreads(); PixelType minimum; PixelType maximum; RealType mean; RealType meanOfPositivePixels; RealType sigma; RealType variance; RealType skewness; RealType kurtosis; RealType sum; sum = sumOfSquares = sumOfCubes = sumOfQuadruples = NumericTraits< RealType >::ZeroValue(); count = countOfPositivePixels = 0; // Find the min/max over all threads and accumulate count, sum and // sum of squares minimum = NumericTraits< PixelType >::max(); maximum = NumericTraits< PixelType >::NonpositiveMin(); for ( i = 0; i < numberOfThreads; i++ ) { count += m_Count[i]; sum += m_ThreadSum[i]; sumOfSquares += m_SumOfSquares[i]; sumOfCubes += m_SumOfCubes[i]; sumOfQuadruples += m_SumOfQuadruples[i]; sumOfPositivePixels += m_ThreadSumOfPositivePixels[i]; countOfPositivePixels += m_PositivePixelCount[i]; if ( m_ThreadMin[i] < minimum ) { minimum = m_ThreadMin[i]; } if ( m_ThreadMax[i] > maximum ) { maximum = m_ThreadMax[i]; } // if enabled, update the histogram for this label if ( m_UseHistogram ) { typename HistogramType::IndexType index; index.SetSize(1); for ( int bin = 0; bin < m_NumBins; ++bin ) { index[0] = bin; m_Histogram->IncreaseFrequency( bin, m_HistogramPerThread[i]->GetFrequency(bin) ); } } } // compute statistics mean = sum / static_cast< RealType >( count ); meanOfPositivePixels = sumOfPositivePixels / static_cast< RealType >( countOfPositivePixels ); // unbiased estimate variance = ( sumOfSquares - ( sum * sum / static_cast< RealType >( count ) ) ) / ( static_cast< RealType >( count ) ); RealType secondMoment, thirdMoment, fourthMoment; secondMoment = sumOfSquares / static_cast< RealType >( count ); thirdMoment = sumOfCubes / static_cast< RealType >( count ); fourthMoment = sumOfQuadruples / static_cast< RealType >( count ); skewness = (thirdMoment - 3. * secondMoment * mean + 2. * std::pow(mean, 3.)) / std::pow(secondMoment - std::pow(mean, 2.), 1.5); // see http://www.boost.org/doc/libs/1_51_0/doc/html/boost/accumulators/impl/skewness_impl.html kurtosis = (fourthMoment - 4. * thirdMoment * mean + 6. * secondMoment * std::pow(mean, 2.) - 3. * std::pow(mean, 4.)) / std::pow(secondMoment - std::pow(mean, 2.), 2.); // see http://www.boost.org/doc/libs/1_46_1/doc/html/boost/accumulators/impl/kurtosis_impl.html, dropped -3 sigma = std::sqrt(variance); // Set the outputs this->GetMinimumOutput()->Set(minimum); this->GetMaximumOutput()->Set(maximum); this->GetMeanOutput()->Set(mean); this->GetSigmaOutput()->Set(sigma); this->GetVarianceOutput()->Set(variance); this->GetSumOutput()->Set(sum); this->GetKurtosisOutput()->Set(kurtosis); this->GetSkewnessOutput()->Set(skewness); this->GetMPPOutput()->Set(meanOfPositivePixels); if (m_UseHistogram) { mitk::HistogramStatisticsCalculator histStatCalc; histStatCalc.SetHistogram(m_Histogram); histStatCalc.CalculateStatistics(); this->GetUniformityOutput()->Set(histStatCalc.GetUniformity()); this->GetUPPOutput()->Set(histStatCalc.GetUPP()); this->GetEntropyOutput()->Set(histStatCalc.GetEntropy()); this->GetMedianOutput()->Set(histStatCalc.GetMedian()); } m_HistogramCalculated = true; } template< typename TInputImage > void ExtendedStatisticsImageFilter< TInputImage > ::SetHistogramParameters(const int numBins, RealType lowerBound, RealType upperBound) { m_NumBins = numBins; m_LowerBound = lowerBound; m_UpperBound = upperBound; m_UseHistogram = true; } } #endif diff --git a/Modules/ImageStatistics/mitkHotspotMaskGenerator.cpp b/Modules/ImageStatistics/mitkHotspotMaskGenerator.cpp index c1b74388eb..5f8a2f6844 100644 --- a/Modules/ImageStatistics/mitkHotspotMaskGenerator.cpp +++ b/Modules/ImageStatistics/mitkHotspotMaskGenerator.cpp @@ -1,603 +1,602 @@ #include #include #include #include #include #include "mitkImageAccessByItk.h" #include #include #include namespace mitk { HotspotMaskGenerator::HotspotMaskGenerator(): m_HotspotRadiusinMM(6.2035049089940), // radius of a 1cm3 sphere in mm m_HotspotMustBeCompletelyInsideImage(true), m_Label(1) { m_TimeStep = 0; m_InternalMask = mitk::Image::New(); m_InternalMaskUpdateTime = 0; } void HotspotMaskGenerator::SetInputImage(mitk::Image::Pointer inputImage) { if (inputImage != m_inputImage) { m_inputImage = inputImage; m_ConvolutionImageMaxIndex.set_size(inputImage->GetDimension()); m_ConvolutionImageMinIndex.set_size(inputImage->GetDimension()); this->Modified(); } } void HotspotMaskGenerator::SetMask(MaskGenerator::Pointer mask) { if (mask != m_Mask) { m_Mask = mask; this->Modified(); } } HotspotMaskGenerator::~HotspotMaskGenerator() { } void HotspotMaskGenerator::SetHotspotRadiusInMM(double radiusInMillimeter) { if(radiusInMillimeter != m_HotspotRadiusinMM) { m_HotspotRadiusinMM = radiusInMillimeter; this->Modified(); } } const double& HotspotMaskGenerator::GetHotspotRadiusinMM() const { return m_HotspotRadiusinMM; } bool HotspotMaskGenerator::GetHotspotMustBeCompletelyInsideImage() const { return m_HotspotMustBeCompletelyInsideImage; } void HotspotMaskGenerator::SetHotspotMustBeCompletelyInsideImage(bool mustBeCompletelyInImage) { if (m_HotspotMustBeCompletelyInsideImage != mustBeCompletelyInImage) { m_HotspotMustBeCompletelyInsideImage = mustBeCompletelyInImage; this->Modified(); } } mitk::Image::Pointer HotspotMaskGenerator::GetMask() { if (IsUpdateRequired()) { if ( m_inputImage.IsNull() ) { throw std::runtime_error( "Error: image empty!" ); } if ( m_TimeStep >= m_inputImage->GetTimeSteps() ) { throw std::runtime_error( "Error: invalid time step!" ); } mitk::ImageTimeSelector::Pointer imageTimeSelector = mitk::ImageTimeSelector::New(); imageTimeSelector->SetInput( m_inputImage ); imageTimeSelector->SetTimeNr( m_TimeStep ); imageTimeSelector->UpdateLargestPossibleRegion(); mitk::Image::Pointer timeSliceImage = imageTimeSelector->GetOutput(); m_internalImage = timeSliceImage; m_internalMask2D = nullptr; // is this correct when this variable holds a smart pointer? m_internalMask3D = nullptr; if ( m_Mask != nullptr ) { m_Mask->SetTimeStep(m_TimeStep); mitk::Image::Pointer timeSliceMask = m_Mask->GetMask(); if ( m_internalImage->GetDimension() == 3 ) { CastToItkImage(timeSliceMask, m_internalMask3D); AccessFixedDimensionByItk_2(m_internalImage, CalculateHotspotMask, 3, m_internalMask3D, m_Label); } else if ( m_internalImage->GetDimension() == 2 ) { CastToItkImage(timeSliceMask, m_internalMask2D); AccessFixedDimensionByItk_2(m_internalImage, CalculateHotspotMask, 2, m_internalMask2D, m_Label); } else { throw std::runtime_error( "Error: invalid image dimension" ); } } else { if ( m_internalImage->GetDimension() == 3 ) { AccessFixedDimensionByItk_2(m_internalImage, CalculateHotspotMask, 3, m_internalMask3D, m_Label); } else if ( m_internalImage->GetDimension() == 2 ) { AccessFixedDimensionByItk_2(m_internalImage, CalculateHotspotMask, 2, m_internalMask2D, m_Label); } else { throw std::runtime_error( "Error: invalid image dimension" ); } } this->Modified(); } m_InternalMaskUpdateTime = m_InternalMask->GetMTime(); return m_InternalMask; } void HotspotMaskGenerator::SetTimeStep(unsigned int timeStep) { if (m_TimeStep != timeStep) { m_TimeStep = timeStep; } } void HotspotMaskGenerator::SetLabel(unsigned short label) { if (label != m_Label) { m_Label = label; this->Modified(); } } vnl_vector HotspotMaskGenerator::GetConvolutionImageMinIndex() { this->GetMask(); // make sure we are up to date return m_ConvolutionImageMinIndex; } vnl_vector HotspotMaskGenerator::GetHotspotIndex() { this->GetMask(); // make sure we are up to date return m_ConvolutionImageMaxIndex; } template HotspotMaskGenerator::ImageExtrema HotspotMaskGenerator::CalculateExtremaWorld( const itk::Image* inputImage, typename itk::Image::Pointer maskImage, double neccessaryDistanceToImageBorderInMM, unsigned int label ) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< unsigned short, VImageDimension > MaskImageType; typedef itk::ImageRegionConstIteratorWithIndex MaskImageIteratorType; typedef itk::ImageRegionConstIteratorWithIndex InputImageIndexIteratorType; typename ImageType::SpacingType spacing = inputImage->GetSpacing(); ImageExtrema minMax; minMax.Defined = false; minMax.MaxIndex.set_size(VImageDimension); minMax.MaxIndex.set_size(VImageDimension); typename ImageType::RegionType allowedExtremaRegion = inputImage->GetLargestPossibleRegion(); bool keepDistanceToImageBorders( neccessaryDistanceToImageBorderInMM > 0 ); if (keepDistanceToImageBorders) { itk::IndexValueType distanceInPixels[VImageDimension]; for(unsigned short dimension = 0; dimension < VImageDimension; ++dimension) { // To confirm that the whole hotspot is inside the image we have to keep a specific distance to the image-borders, which is as long as // the radius. To get the amount of indices we divide the radius by spacing and add 0.5 because voxels are center based: // For example with a radius of 2.2 and a spacing of 1 two indices are enough because 2.2 / 1 + 0.5 = 2.7 => 2. // But with a radius of 2.7 we need 3 indices because 2.7 / 1 + 0.5 = 3.2 => 3 distanceInPixels[dimension] = int( neccessaryDistanceToImageBorderInMM / spacing[dimension] + 0.5); } allowedExtremaRegion.ShrinkByRadius(distanceInPixels); } InputImageIndexIteratorType imageIndexIt(inputImage, allowedExtremaRegion); float maxValue = itk::NumericTraits::min(); float minValue = itk::NumericTraits::max(); typename ImageType::IndexType maxIndex; typename ImageType::IndexType minIndex; for(unsigned short i = 0; i < VImageDimension; ++i) { maxIndex[i] = 0; minIndex[i] = 0; } if (maskImage != nullptr) { MaskImageIteratorType maskIt(maskImage, maskImage->GetLargestPossibleRegion()); typename ImageType::IndexType imageIndex; typename ImageType::PointType worldPosition; typename ImageType::IndexType maskIndex; for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt) { imageIndex = maskIndex = maskIt.GetIndex(); if(maskIt.Get() == label) { if( allowedExtremaRegion.IsInside(imageIndex) ) { imageIndexIt.SetIndex( imageIndex ); double value = imageIndexIt.Get(); minMax.Defined = true; //Calculate minimum, maximum and corresponding index-values if( value > maxValue ) { maxIndex = imageIndexIt.GetIndex(); maxValue = value; } if(value < minValue ) { minIndex = imageIndexIt.GetIndex(); minValue = value; } } } } } else { for(imageIndexIt.GoToBegin(); !imageIndexIt.IsAtEnd(); ++imageIndexIt) { double value = imageIndexIt.Get(); minMax.Defined = true; //Calculate minimum, maximum and corresponding index-values if( value > maxValue ) { maxIndex = imageIndexIt.GetIndex(); maxValue = value; } if(value < minValue ) { minIndex = imageIndexIt.GetIndex(); minValue = value; } } } minMax.MaxIndex.set_size(VImageDimension); minMax.MinIndex.set_size(VImageDimension); for(unsigned int i = 0; i < minMax.MaxIndex.size(); ++i) { minMax.MaxIndex[i] = maxIndex[i]; } for(unsigned int i = 0; i < minMax.MinIndex.size(); ++i) { minMax.MinIndex[i] = minIndex[i]; } minMax.Max = maxValue; minMax.Min = minValue; return minMax; } template itk::Size HotspotMaskGenerator::CalculateConvolutionKernelSize( double spacing[VImageDimension], double radiusInMM ) { typedef itk::Image< float, VImageDimension > KernelImageType; typedef typename KernelImageType::SizeType SizeType; SizeType maskSize; for(unsigned int i = 0; i < VImageDimension; ++i) { maskSize[i] = static_cast( 2 * radiusInMM / spacing[i]); // We always want an uneven size to have a clear center point in the convolution mask if(maskSize[i] % 2 == 0 ) { ++maskSize[i]; } } return maskSize; } template itk::SmartPointer< itk::Image > HotspotMaskGenerator::GenerateHotspotSearchConvolutionKernel(double mmPerPixel[VImageDimension], double radiusInMM ) { std::stringstream ss; for (unsigned int i = 0; i < VImageDimension; ++i) { ss << mmPerPixel[i]; if (i < VImageDimension -1) ss << ","; } MITK_DEBUG << "Update convolution kernel for spacing (" << ss.str() << ") and radius " << radiusInMM << "mm"; double radiusInMMSquared = radiusInMM * radiusInMM; typedef itk::Image< float, VImageDimension > KernelImageType; typename KernelImageType::Pointer convolutionKernel = KernelImageType::New(); // Calculate size and allocate mask image typedef typename KernelImageType::SizeType SizeType; SizeType maskSize = this->CalculateConvolutionKernelSize(mmPerPixel, radiusInMM); mitk::Point3D convolutionMaskCenterIndex; convolutionMaskCenterIndex.Fill(0.0); for(unsigned int i = 0; i < VImageDimension; ++i) { convolutionMaskCenterIndex[i] = 0.5 * (double)(maskSize[i]-1); } typedef typename KernelImageType::IndexType IndexType; IndexType maskIndex; maskIndex.Fill(0); typedef typename KernelImageType::RegionType RegionType; RegionType maskRegion; maskRegion.SetSize(maskSize); maskRegion.SetIndex(maskIndex); convolutionKernel->SetRegions(maskRegion); convolutionKernel->SetSpacing(mmPerPixel); convolutionKernel->Allocate(); // Fill mask image values by subsampling the image grid typedef itk::ImageRegionIteratorWithIndex MaskIteratorType; MaskIteratorType maskIt(convolutionKernel,maskRegion); int numberOfSubVoxelsPerDimension = 2; // per dimension! int numberOfSubVoxels = ::pow( static_cast(numberOfSubVoxelsPerDimension), static_cast(VImageDimension) ); double subVoxelSizeInPixels = 1.0 / (double)numberOfSubVoxelsPerDimension; double valueOfOneSubVoxel = 1.0 / (double)numberOfSubVoxels; double maskValue = 0.0; mitk::Point3D subVoxelIndexPosition; double distanceSquared = 0.0; typedef itk::ContinuousIndex ContinuousIndexType; for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt) { ContinuousIndexType indexPoint(maskIt.GetIndex()); mitk::Point3D voxelPosition; for (unsigned int dimension = 0; dimension < VImageDimension; ++dimension) { voxelPosition[dimension] = indexPoint[dimension]; } maskValue = 0.0; mitk::Vector3D subVoxelOffset; subVoxelOffset.Fill(0.0); // iterate sub-voxels by iterating all possible offsets for (subVoxelOffset[0] = -0.5 + subVoxelSizeInPixels / 2.0; subVoxelOffset[0] < +0.5; subVoxelOffset[0] += subVoxelSizeInPixels) { for (subVoxelOffset[1] = -0.5 + subVoxelSizeInPixels / 2.0; subVoxelOffset[1] < +0.5; subVoxelOffset[1] += subVoxelSizeInPixels) { for (subVoxelOffset[2] = -0.5 + subVoxelSizeInPixels / 2.0; subVoxelOffset[2] < +0.5; subVoxelOffset[2] += subVoxelSizeInPixels) { subVoxelIndexPosition = voxelPosition + subVoxelOffset; // this COULD be integrated into the for-loops if neccessary (add voxelPosition to initializer and end condition) distanceSquared = (subVoxelIndexPosition[0]-convolutionMaskCenterIndex[0]) * mmPerPixel[0] * (subVoxelIndexPosition[0]-convolutionMaskCenterIndex[0]) * mmPerPixel[0] + (subVoxelIndexPosition[1]-convolutionMaskCenterIndex[1]) * mmPerPixel[1] * (subVoxelIndexPosition[1]-convolutionMaskCenterIndex[1]) * mmPerPixel[1] + (subVoxelIndexPosition[2]-convolutionMaskCenterIndex[2]) * mmPerPixel[2] * (subVoxelIndexPosition[2]-convolutionMaskCenterIndex[2]) * mmPerPixel[2]; if (distanceSquared <= radiusInMMSquared) { maskValue += valueOfOneSubVoxel; } } } } maskIt.Set( maskValue ); } return convolutionKernel; } template itk::SmartPointer > HotspotMaskGenerator::GenerateConvolutionImage( const itk::Image* inputImage ) { double mmPerPixel[VImageDimension]; for (unsigned int dimension = 0; dimension < VImageDimension; ++dimension) { mmPerPixel[dimension] = inputImage->GetSpacing()[dimension]; } // update convolution kernel typedef itk::Image< float, VImageDimension > KernelImageType; typename KernelImageType::Pointer convolutionKernel = this->GenerateHotspotSearchConvolutionKernel(mmPerPixel, m_HotspotRadiusinMM); // update convolution image typedef itk::Image< TPixel, VImageDimension > InputImageType; typedef itk::Image< TPixel, VImageDimension > ConvolutionImageType; typedef itk::FFTConvolutionImageFilter ConvolutionFilterType; typename ConvolutionFilterType::Pointer convolutionFilter = ConvolutionFilterType::New(); typedef itk::ConstantBoundaryCondition BoundaryConditionType; BoundaryConditionType boundaryCondition; boundaryCondition.SetConstant(0.0); if (m_HotspotMustBeCompletelyInsideImage) { // overwrite default boundary condition convolutionFilter->SetBoundaryCondition(&boundaryCondition); } convolutionFilter->SetInput(inputImage); convolutionFilter->SetKernelImage(convolutionKernel); convolutionFilter->SetNormalize(true); MITK_DEBUG << "Update Convolution image for hotspot search"; convolutionFilter->UpdateLargestPossibleRegion(); typename ConvolutionImageType::Pointer convolutionImage = convolutionFilter->GetOutput(); convolutionImage->SetSpacing( inputImage->GetSpacing() ); // only workaround because convolution filter seems to ignore spacing of input image return convolutionImage; } template < typename TPixel, unsigned int VImageDimension> void HotspotMaskGenerator::FillHotspotMaskPixels( itk::Image* maskImage, itk::Point sphereCenter, double sphereRadiusInMM ) { typedef itk::Image< TPixel, VImageDimension > MaskImageType; typedef itk::ImageRegionIteratorWithIndex MaskImageIteratorType; MaskImageIteratorType maskIt(maskImage, maskImage->GetLargestPossibleRegion()); typename MaskImageType::IndexType maskIndex; typename MaskImageType::PointType worldPosition; // this is not very smart. I would rather use a 0 initialized mask (not the case here -> blame CalculateHotspotMask) and find the region where I need to iterate over, then iterate only over the small region for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt) { maskIndex = maskIt.GetIndex(); maskImage->TransformIndexToPhysicalPoint(maskIndex, worldPosition); maskIt.Set( worldPosition.EuclideanDistanceTo(sphereCenter) <= sphereRadiusInMM ? 1 : 0 ); } } template void HotspotMaskGenerator::CalculateHotspotMask(itk::Image* inputImage, typename itk::Image::Pointer maskImage, unsigned int label) { typedef itk::Image< TPixel, VImageDimension > InputImageType; typedef itk::Image< TPixel, VImageDimension > ConvolutionImageType; - typedef itk::Image< float, VImageDimension > KernelImageType; typedef itk::Image< unsigned short, VImageDimension > MaskImageType; typename ConvolutionImageType::Pointer convolutionImage = this->GenerateConvolutionImage(inputImage); if (convolutionImage.IsNull()) { MITK_ERROR << "Empty convolution image in CalculateHotspotStatistics(). We should never reach this state (logic error)."; throw std::logic_error("Empty convolution image in CalculateHotspotStatistics()"); } // if mask image is not defined, create an image of the same size as inputImage and fill it with 1's // there is maybe a better way to do this!? if (maskImage == nullptr) { maskImage = MaskImageType::New(); typename MaskImageType::RegionType maskRegion = inputImage->GetLargestPossibleRegion(); typename MaskImageType::SpacingType maskSpacing = inputImage->GetSpacing(); typename MaskImageType::PointType maskOrigin = inputImage->GetOrigin(); typename MaskImageType::DirectionType maskDirection = inputImage->GetDirection(); maskImage->SetRegions(maskRegion); maskImage->Allocate(); maskImage->SetOrigin(maskOrigin); maskImage->SetSpacing(maskSpacing); maskImage->SetDirection(maskDirection); maskImage->FillBuffer(1); label = 1; } // find maximum in convolution image, given the current mask double requiredDistanceToBorder = m_HotspotMustBeCompletelyInsideImage ? m_HotspotRadiusinMM : -1.0; ImageExtrema convolutionImageInformation = CalculateExtremaWorld(convolutionImage.GetPointer(), maskImage, requiredDistanceToBorder, label); bool isHotspotDefined = convolutionImageInformation.Defined; if (!isHotspotDefined) { MITK_ERROR << "No origin of hotspot-sphere was calculated!"; m_InternalMask = nullptr; } else { // create a binary mask around the "hotspot" region, fill the shape of a sphere around our hotspot center // typename DuplicatorType::Pointer copyMachine = DuplicatorType::New(); // copyMachine->SetInputImage(inputImage); // copyMachine->Update(); // typename CastFilterType::Pointer caster = CastFilterType::New(); // caster->SetInput( copyMachine->GetOutput() ); // caster->Update(); typename MaskImageType::Pointer hotspotMaskITK = MaskImageType::New(); hotspotMaskITK->SetOrigin(inputImage->GetOrigin()); hotspotMaskITK->SetSpacing(inputImage->GetSpacing()); hotspotMaskITK->SetLargestPossibleRegion(inputImage->GetLargestPossibleRegion()); hotspotMaskITK->SetBufferedRegion(inputImage->GetBufferedRegion()); hotspotMaskITK->SetDirection(inputImage->GetDirection()); hotspotMaskITK->SetNumberOfComponentsPerPixel(inputImage->GetNumberOfComponentsPerPixel()); hotspotMaskITK->Allocate(); hotspotMaskITK->FillBuffer(1); typedef typename InputImageType::IndexType IndexType; IndexType maskCenterIndex; for (unsigned int d =0; d< VImageDimension;++d) { maskCenterIndex[d]=convolutionImageInformation.MaxIndex[d]; } typename ConvolutionImageType::PointType maskCenter; inputImage->TransformIndexToPhysicalPoint(maskCenterIndex,maskCenter); FillHotspotMaskPixels(hotspotMaskITK.GetPointer(), maskCenter, m_HotspotRadiusinMM); //obtain mitk::Image::Pointer from itk::Image mitk::Image::Pointer hotspotMaskAsMITKImage = mitk::GrabItkImageMemory(hotspotMaskITK); m_InternalMask = hotspotMaskAsMITKImage; m_ConvolutionImageMaxIndex = convolutionImageInformation.MaxIndex; m_ConvolutionImageMinIndex = convolutionImageInformation.MinIndex; } } bool HotspotMaskGenerator::IsUpdateRequired() const { unsigned long thisClassTimeStamp = this->GetMTime(); unsigned long internalMaskTimeStamp = m_InternalMask->GetMTime(); unsigned long maskGeneratorTimeStamp = m_Mask->GetMTime(); unsigned long inputImageTimeStamp = m_inputImage->GetMTime(); if (thisClassTimeStamp > m_InternalMaskUpdateTime) // inputs have changed { return true; } if (m_InternalMaskUpdateTime < maskGeneratorTimeStamp || m_InternalMaskUpdateTime < inputImageTimeStamp) // mask image has changed outside of this class { return true; } if (internalMaskTimeStamp > m_InternalMaskUpdateTime) // internal mask has been changed outside of this class { return true; } return false; } } diff --git a/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp b/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp index c659d8f075..71034570bd 100644 --- a/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp +++ b/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp @@ -1,642 +1,629 @@ #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include "itkImageFileWriter.h" namespace mitk { void ImageStatisticsCalculator::SetInputImage(mitk::Image::Pointer image) { if (image != m_Image) { m_Image = image; m_StatisticsByTimeStep.resize(m_Image->GetTimeSteps()); m_StatisticsUpdateTimePerTimeStep.resize(m_Image->GetTimeSteps()); std::fill(m_StatisticsUpdateTimePerTimeStep.begin(), m_StatisticsUpdateTimePerTimeStep.end(), 0); this->Modified(); } } void ImageStatisticsCalculator::SetMask(mitk::MaskGenerator::Pointer mask) { if (mask != m_MaskGenerator) { m_MaskGenerator = mask; this->Modified(); } } void ImageStatisticsCalculator::SetSecondaryMask(mitk::MaskGenerator::Pointer mask) { if (mask != m_SecondaryMaskGenerator) { m_SecondaryMaskGenerator = mask; this->Modified(); } } void ImageStatisticsCalculator::SetNBinsForHistogramStatistics(unsigned int nBins) { if (nBins != m_nBinsForHistogramStatistics) { m_nBinsForHistogramStatistics = nBins; this->Modified(); this->m_UseBinSizeOverNBins = false; } if (m_UseBinSizeOverNBins) { this->Modified(); this->m_UseBinSizeOverNBins = false; } } unsigned int ImageStatisticsCalculator::GetNBinsForHistogramStatistics() const { return m_nBinsForHistogramStatistics; } void ImageStatisticsCalculator::SetBinSizeForHistogramStatistics(double binSize) { if (binSize != m_binSizeForHistogramStatistics) { m_binSizeForHistogramStatistics = binSize; this->Modified(); this->m_UseBinSizeOverNBins = true; } if (!m_UseBinSizeOverNBins) { this->Modified(); this->m_UseBinSizeOverNBins = true; } } double ImageStatisticsCalculator::GetBinSizeForHistogramStatistics() const { return m_binSizeForHistogramStatistics; } ImageStatisticsCalculator::StatisticsContainer::Pointer ImageStatisticsCalculator::GetStatistics(unsigned int timeStep, unsigned int label) { if (timeStep >= m_StatisticsByTimeStep.size()) { mitkThrow() << "invalid timeStep in ImageStatisticsCalculator_v2::GetStatistics"; } if (m_Image.IsNull()) { mitkThrow() << "no image"; } if (!m_Image->IsInitialized()) { mitkThrow() << "Image not initialized!"; } if (IsUpdateRequired(timeStep)) { if (m_MaskGenerator.IsNotNull()) { m_MaskGenerator->SetTimeStep(timeStep); m_InternalMask = m_MaskGenerator->GetMask(); if (m_MaskGenerator->GetReferenceImage().IsNotNull()) { m_InternalImageForStatistics = m_MaskGenerator->GetReferenceImage(); } else { m_InternalImageForStatistics = m_Image; } } else { m_InternalImageForStatistics = m_Image; } if (m_SecondaryMaskGenerator.IsNotNull()) { m_SecondaryMaskGenerator->SetTimeStep(timeStep); m_SecondaryMask = m_SecondaryMaskGenerator->GetMask(); } ImageTimeSelector::Pointer imgTimeSel = ImageTimeSelector::New(); imgTimeSel->SetInput(m_InternalImageForStatistics); imgTimeSel->SetTimeNr(timeStep); imgTimeSel->UpdateLargestPossibleRegion(); m_ImageTimeSlice = imgTimeSel->GetOutput(); // Calculate statistics with/without mask if (m_MaskGenerator.IsNull() && m_SecondaryMaskGenerator.IsNull()) { // 1) calculate statistics unmasked: AccessByItk_1(m_ImageTimeSlice, InternalCalculateStatisticsUnmasked, timeStep) } else { // 2) calculate statistics masked AccessByItk_1(m_ImageTimeSlice, InternalCalculateStatisticsMasked, timeStep) } //this->Modified(); } m_StatisticsUpdateTimePerTimeStep[timeStep] = m_StatisticsByTimeStep[timeStep][m_StatisticsByTimeStep[timeStep].size()-1]->GetMTime(); for (std::vector::iterator it = m_StatisticsByTimeStep[timeStep].begin(); it != m_StatisticsByTimeStep[timeStep].end(); ++it) { StatisticsContainer::Pointer statCont = *it; if (statCont->GetLabel() == label) { return statCont->Clone(); } } // these lines will ony be executed if the requested label could not be found! MITK_WARN << "Invalid label: " << label << " in time step: " << timeStep; return StatisticsContainer::New(); } template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator::InternalCalculateStatisticsUnmasked( typename itk::Image< TPixel, VImageDimension >* image, unsigned int timeStep) { typedef typename itk::Image< TPixel, VImageDimension > ImageType; typedef typename itk::ExtendedStatisticsImageFilter ImageStatisticsFilterType; typedef typename itk::MinMaxImageFilterWithIndex MinMaxFilterType; StatisticsContainer::Pointer statisticsResult = StatisticsContainer::New(); typename ImageStatisticsFilterType::Pointer statisticsFilter = ImageStatisticsFilterType::New(); statisticsFilter->SetInput(image); statisticsFilter->SetCoordinateTolerance(0.001); statisticsFilter->SetDirectionTolerance(0.001); // TODO: this is single threaded. Implement our own image filter that does this multi threaded // typename itk::MinimumMaximumImageCalculator::Pointer imgMinMaxFilter = itk::MinimumMaximumImageCalculator::New(); // imgMinMaxFilter->SetImage(image); // imgMinMaxFilter->Compute(); vnl_vector minIndex, maxIndex; typename MinMaxFilterType::Pointer minMaxFilter = MinMaxFilterType::New(); minMaxFilter->SetInput(image); minMaxFilter->UpdateLargestPossibleRegion(); typename ImageType::PixelType minval = minMaxFilter->GetMin(); typename ImageType::PixelType maxval = minMaxFilter->GetMax(); typename ImageType::IndexType tmpMinIndex = minMaxFilter->GetMinIndex(); typename ImageType::IndexType tmpMaxIndex = minMaxFilter->GetMaxIndex(); // typename ImageType::IndexType tmpMinIndex = imgMinMaxFilter->GetIndexOfMinimum(); // typename ImageType::IndexType tmpMaxIndex = imgMinMaxFilter->GetIndexOfMaximum(); minIndex.set_size(tmpMaxIndex.GetIndexDimension()); maxIndex.set_size(tmpMaxIndex.GetIndexDimension()); for (unsigned int i=0; i < tmpMaxIndex.GetIndexDimension(); i++) { minIndex[i] = tmpMinIndex[i]; maxIndex[i] = tmpMaxIndex[i]; } statisticsResult->SetMinIndex(minIndex); statisticsResult->SetMaxIndex(maxIndex); //convert m_binSize in m_nBins if necessary unsigned int nBinsForHistogram; if (m_UseBinSizeOverNBins) { nBinsForHistogram = std::max(static_cast(std::ceil(maxval - minval)) / m_binSizeForHistogramStatistics, 10.); // do not allow less than 10 bins } else { nBinsForHistogram = m_nBinsForHistogramStatistics; } statisticsFilter->SetHistogramParameters(nBinsForHistogram, minval, maxval); try { statisticsFilter->Update(); } catch (const itk::ExceptionObject& e) { mitkThrow() << "Image statistics calculation failed due to following ITK Exception: \n " << e.what(); } // no mask, therefore just one label = the whole image m_StatisticsByTimeStep[timeStep].resize(1); statisticsResult->SetLabel(1); statisticsResult->SetN(image->GetLargestPossibleRegion().GetNumberOfPixels()); statisticsResult->SetMean(statisticsFilter->GetMean()); statisticsResult->SetMin(statisticsFilter->GetMinimum()); statisticsResult->SetMax(statisticsFilter->GetMaximum()); statisticsResult->SetVariance(statisticsFilter->GetVariance()); statisticsResult->SetStd(statisticsFilter->GetSigma()); statisticsResult->SetSkewness(statisticsFilter->GetSkewness()); statisticsResult->SetKurtosis(statisticsFilter->GetKurtosis()); statisticsResult->SetRMS(std::sqrt(std::pow(statisticsFilter->GetMean(), 2.) + statisticsFilter->GetVariance())); // variance = sigma^2 statisticsResult->SetMPP(statisticsFilter->GetMPP()); statisticsResult->SetEntropy(statisticsFilter->GetEntropy()); statisticsResult->SetMedian(statisticsFilter->GetMedian()); statisticsResult->SetUniformity(statisticsFilter->GetUniformity()); statisticsResult->SetUPP(statisticsFilter->GetUPP()); statisticsResult->SetHistogram(statisticsFilter->GetHistogram()); m_StatisticsByTimeStep[timeStep][0] = statisticsResult; } template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator::InternalCalculateStatisticsMasked( typename itk::Image< TPixel, VImageDimension >* image, unsigned int timeStep) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< MaskPixelType, VImageDimension > MaskType; typedef typename MaskType::PixelType LabelPixelType; typedef itk::ExtendedLabelStatisticsImageFilter< ImageType, MaskType > ImageStatisticsFilterType; typedef MaskUtilities< TPixel, VImageDimension > MaskUtilType; typedef typename itk::MinMaxLabelImageFilterWithIndex MinMaxLabelFilterType; typedef typename ImageType::PixelType InputImgPixelType; // workaround: if m_SecondaryMaskGenerator ist not null but m_MaskGenerator is! (this is the case if we request a 'ignore zuero valued pixels' // mask in the gui but do not define a primary mask) bool swapMasks = false; if (m_SecondaryMask.IsNotNull() && m_InternalMask.IsNull()) { m_InternalMask = m_SecondaryMask; m_SecondaryMask = nullptr; swapMasks = true; } // maskImage has to have the same dimension as image typename MaskType::Pointer maskImage = MaskType::New(); try { // try to access the pixel values directly (no copying or casting). Only works if mask pixels are of pixelType unsigned short maskImage = ImageToItkImage< MaskPixelType, VImageDimension >(m_InternalMask); } catch (const itk::ExceptionObject &) { // if the pixel type of the mask is not short, then we have to make a copy of m_InternalMask (and cast the values) CastToItkImage(m_InternalMask, maskImage); } // if we have a secondary mask (say a ignoreZeroPixelMask) we need to combine the masks (corresponds to AND) if (m_SecondaryMask.IsNotNull()) { // dirty workaround for a bug when pf mask + any other mask is used in conjunction. We need a proper fix for this (Fabian Isensee is responsible and probably working on it!) if (m_InternalMask->GetDimension() == 2 && (m_SecondaryMask->GetDimension() == 3 || m_SecondaryMask->GetDimension() == 4)) { mitk::Image::Pointer old_img = m_SecondaryMaskGenerator->GetReferenceImage(); m_SecondaryMaskGenerator->SetInputImage(m_MaskGenerator->GetReferenceImage()); m_SecondaryMask = m_SecondaryMaskGenerator->GetMask(); m_SecondaryMaskGenerator->SetInputImage(old_img); } typename MaskType::Pointer secondaryMaskImage = MaskType::New(); secondaryMaskImage = ImageToItkImage< MaskPixelType, VImageDimension >(m_SecondaryMask); // secondary mask should be a ignore zero value pixel mask derived from image. it has to be cropped to the mask region (which may be planar or simply smaller) typename MaskUtilities::Pointer secondaryMaskMaskUtil = MaskUtilities::New(); secondaryMaskMaskUtil->SetImage(secondaryMaskImage.GetPointer()); secondaryMaskMaskUtil->SetMask(maskImage.GetPointer()); typename MaskType::Pointer adaptedSecondaryMaskImage = secondaryMaskMaskUtil->ExtractMaskImageRegion(); typename itk::MaskImageFilter2::Pointer maskFilter = itk::MaskImageFilter2::New(); maskFilter->SetInput1(maskImage); maskFilter->SetInput2(adaptedSecondaryMaskImage); maskFilter->SetMaskingValue(1); // all pixels of maskImage where secondaryMaskImage==1 will be kept, all the others are set to 0 maskFilter->UpdateLargestPossibleRegion(); maskImage = maskFilter->GetOutput(); } typename MaskUtilType::Pointer maskUtil = MaskUtilType::New(); maskUtil->SetImage(image); maskUtil->SetMask(maskImage.GetPointer()); // if mask is smaller than image, extract the image region where the mask is typename ImageType::Pointer adaptedImage = ImageType::New(); adaptedImage = maskUtil->ExtractMaskImageRegion(); // this also checks mask sanity // find min, max, minindex and maxindex typename MinMaxLabelFilterType::Pointer minMaxFilter = MinMaxLabelFilterType::New(); minMaxFilter->SetInput(adaptedImage); minMaxFilter->SetLabelInput(maskImage); minMaxFilter->UpdateLargestPossibleRegion(); // set histogram parameters for each label individually (min/max may be different for each label) typedef typename std::map MapType; typedef typename std::pair PairType; std::vector relevantLabels = minMaxFilter->GetRelevantLabels(); MapType minVals; MapType maxVals; std::map nBins; for (LabelPixelType label:relevantLabels) { minVals.insert(PairType(label, minMaxFilter->GetMin(label))); maxVals.insert(PairType(label, minMaxFilter->GetMax(label))); unsigned int nBinsForHistogram; if (m_UseBinSizeOverNBins) { nBinsForHistogram = std::max(static_cast(std::ceil(minMaxFilter->GetMax(label) - minMaxFilter->GetMin(label))) / m_binSizeForHistogramStatistics, 10.); // do not allow less than 10 bins } else { nBinsForHistogram = m_nBinsForHistogramStatistics; } nBins.insert(typename std::pair(label, nBinsForHistogram)); } typename ImageStatisticsFilterType::Pointer imageStatisticsFilter = ImageStatisticsFilterType::New(); imageStatisticsFilter->SetDirectionTolerance(0.001); imageStatisticsFilter->SetCoordinateTolerance(0.001); imageStatisticsFilter->SetInput(adaptedImage); imageStatisticsFilter->SetLabelInput(maskImage); imageStatisticsFilter->SetHistogramParametersForLabels(nBins, minVals, maxVals); imageStatisticsFilter->Update(); std::list labels = imageStatisticsFilter->GetRelevantLabels(); std::list::iterator it = labels.begin(); m_StatisticsByTimeStep[timeStep].resize(0); while(it != labels.end()) { StatisticsContainer::Pointer statisticsResult = StatisticsContainer::New(); // find min, max, minindex and maxindex // make sure to only look in the masked region, use a masker for this vnl_vector minIndex, maxIndex; mitk::Point3D worldCoordinateMin; mitk::Point3D worldCoordinateMax; mitk::Point3D indexCoordinateMin; mitk::Point3D indexCoordinateMax; m_InternalImageForStatistics->GetGeometry()->IndexToWorld(minMaxFilter->GetMinIndex(*it), worldCoordinateMin); m_InternalImageForStatistics->GetGeometry()->IndexToWorld(minMaxFilter->GetMaxIndex(*it), worldCoordinateMax); m_Image->GetGeometry()->WorldToIndex(worldCoordinateMin, indexCoordinateMin); m_Image->GetGeometry()->WorldToIndex(worldCoordinateMax, indexCoordinateMax); - //typename ImageType::IndexType tmpMinIndex = minMaxFilter->GetMinIndex(*it); - //typename ImageType::IndexType tmpMaxIndex = minMaxFilter->GetMaxIndex(*it); - - //minIndex.set_size(tmpMaxIndex.GetIndexDimension()); - //maxIndex.set_size(tmpMaxIndex.GetIndexDimension()); minIndex.set_size(3); maxIndex.set_size(3); //for (unsigned int i=0; i < tmpMaxIndex.GetIndexDimension(); i++) for (unsigned int i=0; i < 3; i++) { - //minIndex[i] = tmpMinIndex[i] + (maskImage->GetOrigin()[i] - image->GetOrigin()[i]) / (double) maskImage->GetSpacing()[i]; - //maxIndex[i] = tmpMaxIndex[i] + (maskImage->GetOrigin()[i] - image->GetOrigin()[i]) / (double) maskImage->GetSpacing()[i]; minIndex[i] = indexCoordinateMin[i]; maxIndex[i] = indexCoordinateMax[i]; } statisticsResult->SetMinIndex(minIndex); statisticsResult->SetMaxIndex(maxIndex); - // just debug - TPixel min_Filter = minMaxFilter->GetMin(*it); - TPixel max_Filter = minMaxFilter->GetMax(*it); - TPixel min_Itk = imageStatisticsFilter->GetMinimum(*it); - TPixel max_Itk = imageStatisticsFilter->GetMaximum(*it); - assert(abs(minMaxFilter->GetMax(*it) - imageStatisticsFilter->GetMaximum(*it)) < mitk::eps); assert(abs(minMaxFilter->GetMin(*it) - imageStatisticsFilter->GetMinimum(*it)) < mitk::eps); statisticsResult->SetN(imageStatisticsFilter->GetSum(*it) / (double) imageStatisticsFilter->GetMean(*it)); statisticsResult->SetMean(imageStatisticsFilter->GetMean(*it)); statisticsResult->SetMin(imageStatisticsFilter->GetMinimum(*it)); statisticsResult->SetMax(imageStatisticsFilter->GetMaximum(*it)); statisticsResult->SetVariance(imageStatisticsFilter->GetVariance(*it)); statisticsResult->SetStd(imageStatisticsFilter->GetSigma(*it)); statisticsResult->SetSkewness(imageStatisticsFilter->GetSkewness(*it)); statisticsResult->SetKurtosis(imageStatisticsFilter->GetKurtosis(*it)); statisticsResult->SetRMS(std::sqrt(std::pow(imageStatisticsFilter->GetMean(*it), 2.) + imageStatisticsFilter->GetVariance(*it))); // variance = sigma^2 statisticsResult->SetMPP(imageStatisticsFilter->GetMPP(*it)); statisticsResult->SetLabel(*it); statisticsResult->SetEntropy(imageStatisticsFilter->GetEntropy(*it)); statisticsResult->SetMedian(imageStatisticsFilter->GetMedian(*it)); statisticsResult->SetUniformity(imageStatisticsFilter->GetUniformity(*it)); statisticsResult->SetUPP(imageStatisticsFilter->GetUPP(*it)); statisticsResult->SetHistogram(imageStatisticsFilter->GetHistogram(*it)); m_StatisticsByTimeStep[timeStep].push_back(statisticsResult); ++it; } // swap maskGenerators back if (swapMasks) { m_SecondaryMask = m_InternalMask; m_InternalMask = nullptr; } } bool ImageStatisticsCalculator::IsUpdateRequired(unsigned int timeStep) const { unsigned long thisClassTimeStamp = this->GetMTime(); unsigned long inputImageTimeStamp = m_Image->GetMTime(); unsigned long statisticsTimeStamp = m_StatisticsUpdateTimePerTimeStep[timeStep]; if (thisClassTimeStamp > statisticsTimeStamp) // inputs have changed { return true; } if (inputImageTimeStamp > statisticsTimeStamp) // image has changed { return true; } if (m_MaskGenerator.IsNotNull()) { unsigned long maskGeneratorTimeStamp = m_MaskGenerator->GetMTime(); if (maskGeneratorTimeStamp > statisticsTimeStamp) // there is a mask generator and it has changed { return true; } } if (m_SecondaryMaskGenerator.IsNotNull()) { unsigned long maskGeneratorTimeStamp = m_SecondaryMaskGenerator->GetMTime(); if (maskGeneratorTimeStamp > statisticsTimeStamp) // there is a secondary mask generator and it has changed { return true; } } return false; } ImageStatisticsCalculator::StatisticsContainer::StatisticsContainer(): - m_N(nan("")), + m_N(0), m_Mean(nan("")), m_Min(nan("")), m_Max(nan("")), m_Std(nan("")), m_Variance(nan("")), m_Skewness(nan("")), m_Kurtosis(nan("")), m_RMS(nan("")), m_MPP(nan("")), m_Median(nan("")), m_Uniformity(nan("")), m_UPP(nan("")), m_Entropy(nan("")) { m_minIndex.set_size(0); m_maxIndex.set_size(0); } ImageStatisticsCalculator::statisticsMapType ImageStatisticsCalculator::StatisticsContainer::GetStatisticsAsMap() { ImageStatisticsCalculator::statisticsMapType statisticsAsMap; statisticsAsMap["N"] = m_N; statisticsAsMap["Mean"] = m_Mean; statisticsAsMap["Min"] = m_Min; statisticsAsMap["Max"] = m_Max; statisticsAsMap["StandardDeviation"] = m_Std; statisticsAsMap["Variance"] = m_Variance; statisticsAsMap["Skewness"] = m_Skewness; statisticsAsMap["Kurtosis"] = m_Kurtosis; statisticsAsMap["RMS"] = m_RMS; statisticsAsMap["MPP"] = m_MPP; statisticsAsMap["Median"] = m_Median; statisticsAsMap["Uniformity"] = m_Uniformity; statisticsAsMap["UPP"] = m_UPP; statisticsAsMap["Entropy"] = m_Entropy; statisticsAsMap["Label"] = m_Label; return statisticsAsMap; } void ImageStatisticsCalculator::StatisticsContainer::Reset() { - m_N = nan(""); + m_N = 0; m_Mean = nan(""); m_Min = nan(""); m_Max = nan(""); m_Std = nan(""); m_Variance = nan(""); m_Skewness = nan(""); m_Kurtosis = nan(""); m_RMS = nan(""); m_MPP = nan(""); m_Median = nan(""); m_Uniformity = nan(""); m_UPP = nan(""); m_Entropy = nan(""); m_Histogram = HistogramType::New(); m_minIndex.set_size(0); m_maxIndex.set_size(0); m_Label = 0; } void ImageStatisticsCalculator::StatisticsContainer::Print() { ImageStatisticsCalculator::statisticsMapType statMap = this->GetStatisticsAsMap(); // print all map key value pairs // const auto& val:statMap for (auto it = statMap.begin(); it != statMap.end(); ++it) { std::cout << it->first << ": " << it->second << std::endl; } // print the min and max index std::cout << "Min Index:" << std::endl; for (auto it = this->GetMinIndex().begin(); it != this->GetMinIndex().end(); ++it) { std::cout << *it << " "; } std::cout << std::endl; // print the min and max index std::cout << "Max Index:" << std::endl; for (auto it = this->GetMaxIndex().begin(); it != this->GetMaxIndex().end(); ++it) { std::cout << *it << " "; } std::cout << std::endl; } std::string ImageStatisticsCalculator::StatisticsContainer::GetAsString() { std::string res = ""; ImageStatisticsCalculator::statisticsMapType statMap = this->GetStatisticsAsMap(); // print all map key value pairs // const auto& val:statMap for (auto it = statMap.begin(); it != statMap.end(); ++it) { res += std::string(it->first) + ": " + std::to_string(it->second) + "\n"; } // print the min and max index res += "Min Index:" + std::string("\n"); for (auto it = this->GetMinIndex().begin(); it != this->GetMinIndex().end(); it++) { res += std::to_string(*it) + std::string(" "); } res += "\n"; // print the min and max index res += "Max Index:" + std::string("\n"); for (auto it = this->GetMaxIndex().begin(); it != this->GetMaxIndex().end(); it++) { res += std::to_string(*it) + " "; } res += "\n"; return res; } } diff --git a/Modules/ImageStatistics/mitkMaskUtilities.cpp b/Modules/ImageStatistics/mitkMaskUtilities.cpp index 38f82022ac..9da8e3b2aa 100644 --- a/Modules/ImageStatistics/mitkMaskUtilities.cpp +++ b/Modules/ImageStatistics/mitkMaskUtilities.cpp @@ -1,189 +1,189 @@ #ifndef MITKMASKUTIL_CPP #define MITKMASKUTIL_CPP #include //#include #include #include #include #include namespace mitk { template void MaskUtilities::SetImage(ImageType* image) { if (image != m_Image) { m_Image = image; } } template void MaskUtilities::SetMask(MaskType* mask) { if (mask != m_Mask) { m_Mask = mask; } } template bool MaskUtilities::CheckMaskSanity() { if (m_Mask==nullptr || m_Image==nullptr) { MITK_ERROR << "Set an image and a mask first"; } typedef itk::Image< TPixel, VImageDimension > ImageType; typedef typename ImageType::PointType PointType; typedef typename ImageType::DirectionType DirectionType; bool maskSanity = true; if (m_Mask==nullptr) { MITK_ERROR << "Something went wrong when casting the mitk mask image to an itk mask image. Do the mask and the input image have the same dimension?"; // note to self: We could try to convert say a 2d mask to a 3d mask if the image is 3d. (mask and image dimension have to match.) } // check direction DirectionType imageDirection = m_Image->GetDirection(); DirectionType maskDirection = m_Mask->GetDirection(); - for( int i = 0; i < imageDirection.ColumnDimensions; ++i ) + for(unsigned int i = 0; i < imageDirection.ColumnDimensions; ++i ) { - for( int j = 0; j < imageDirection.ColumnDimensions; ++j ) + for(unsigned int j = 0; j < imageDirection.ColumnDimensions; ++j ) { double differenceDirection = imageDirection[i][j] - maskDirection[i][j]; if ( fabs( differenceDirection ) > mitk::eps ) { double differenceDirection = imageDirection[i][j] - maskDirection[i][j]; if ( fabs( differenceDirection ) > 0.001 /*mitk::eps*/ ) // TODO: temp fix (bug 17121) { maskSanity = false; MITK_INFO << "Mask needs to have same direction as image! (Image direction: " << imageDirection << "; Mask direction: " << maskDirection << ")"; } } } } // check spacing PointType imageSpacing = m_Image->GetSpacing(); PointType maskSpacing = m_Mask->GetSpacing(); for (unsigned int i = 0; i < VImageDimension; i++) { if ( fabs( maskSpacing[i] - imageSpacing[i] ) > mitk::eps ) { maskSanity = false; MITK_INFO << "Spacing of mask and image is not equal. Mask: " << maskSpacing << " image: " << imageSpacing; } } // check alignment // Make sure that the voxels of mask and image are correctly "aligned", i.e., voxel boundaries are the same in both images PointType imageOrigin = m_Image->GetOrigin(); PointType maskOrigin = m_Mask->GetOrigin(); typedef itk::ContinuousIndex ContinousIndexType; ContinousIndexType maskOriginContinousIndex, imageOriginContinousIndex; m_Image->TransformPhysicalPointToContinuousIndex(maskOrigin, maskOriginContinousIndex); m_Image->TransformPhysicalPointToContinuousIndex(imageOrigin, imageOriginContinousIndex); for ( unsigned int i = 0; i < ImageType::ImageDimension; ++i ) { double misalignment = maskOriginContinousIndex[i] - floor( maskOriginContinousIndex[i] + 0.5 ); // misalignment must be a multiple (int) of spacing in that direction if ( fmod(misalignment,imageSpacing[i]) > mitk::eps ) { maskSanity = false; MITK_INFO << "Pixels/voxels of mask and image are not sufficiently aligned! (Misalignment: " << fmod(misalignment,imageSpacing[i]) << ")"; } } // mask must be completely inside image region // Make sure that mask region is contained within image region if ( m_Mask!=nullptr && !m_Image->GetLargestPossibleRegion().IsInside( m_Mask->GetLargestPossibleRegion() ) ) { maskSanity = false; MITK_INFO << "Mask region needs to be inside of image region! (Image region: " << m_Image->GetLargestPossibleRegion() << "; Mask region: " << m_Mask->GetLargestPossibleRegion() << ")"; } return maskSanity; } template typename itk::Image::Pointer MaskUtilities::ExtractMaskImageRegion() { if (m_Mask==nullptr || m_Image==nullptr) { MITK_ERROR << "Set an image and a mask first"; } bool maskSanity = CheckMaskSanity(); if (!maskSanity) { MITK_ERROR << "Mask and image are not compatible"; } typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< unsigned short, VImageDimension > MaskType; typedef itk::ExtractImageFilter< ImageType, ImageType > ExtractImageFilterType; typename ImageType::SizeType imageSize = m_Image->GetBufferedRegion().GetSize(); typename ImageType::SizeType maskSize = m_Mask->GetBufferedRegion().GetSize(); typename itk::Image::Pointer extractedImg = itk::Image::New(); bool maskSmallerImage = false; for ( unsigned int i = 0; i < ImageType::ImageDimension; ++i ) { if ( maskSize[i] < imageSize[i] ) { maskSmallerImage = true; } } if ( maskSmallerImage ) { typename ExtractImageFilterType::Pointer extractImageFilter = ExtractImageFilterType::New(); typename MaskType::PointType maskOrigin = m_Mask->GetOrigin(); typename ImageType::PointType imageOrigin = m_Image->GetOrigin(); typename MaskType::SpacingType maskSpacing = m_Mask->GetSpacing(); typename ImageType::RegionType extractionRegion; typename ImageType::IndexType extractionRegionIndex; for (unsigned int i=0; i < maskOrigin.GetPointDimension(); i++) { extractionRegionIndex[i] = (maskOrigin[i] - imageOrigin[i]) / maskSpacing[i]; } extractionRegion.SetIndex(extractionRegionIndex); extractionRegion.SetSize(m_Mask->GetLargestPossibleRegion().GetSize()); extractImageFilter->SetInput( m_Image ); extractImageFilter->SetExtractionRegion( extractionRegion ); extractImageFilter->SetCoordinateTolerance( 0.001 ); extractImageFilter->SetDirectionTolerance( 0.001 ); extractImageFilter->Update(); extractedImg = extractImageFilter->GetOutput(); extractedImg->SetOrigin(m_Mask->GetOrigin()); extractedImg->SetLargestPossibleRegion(m_Mask->GetLargestPossibleRegion()); extractedImg->SetBufferedRegion(m_Mask->GetBufferedRegion()); } else { extractedImg = m_Image; } return extractedImg; } } #endif diff --git a/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.cpp b/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.cpp index 4f7b19990d..d6d4cddcd6 100644 --- a/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.cpp +++ b/Modules/ImageStatistics/mitkPlanarFigureMaskGenerator.cpp @@ -1,426 +1,425 @@ #include #include #include #include "mitkImageAccessByItk.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include namespace mitk { void PlanarFigureMaskGenerator::SetPlanarFigure(mitk::PlanarFigure::Pointer planarFigure) { if ( planarFigure.IsNull() ) { throw std::runtime_error( "Error: planar figure empty!" ); } if ( !planarFigure->IsClosed() ) { throw std::runtime_error( "Masking not possible for non-closed figures" ); } const PlaneGeometry *planarFigurePlaneGeometry = planarFigure->GetPlaneGeometry(); if ( planarFigurePlaneGeometry == nullptr ) { throw std::runtime_error( "Planar-Figure not yet initialized!" ); } const PlaneGeometry *planarFigureGeometry = dynamic_cast< const PlaneGeometry * >( planarFigurePlaneGeometry ); if ( planarFigureGeometry == nullptr ) { throw std::runtime_error( "Non-planar planar figures not supported!" ); } if (planarFigure != m_PlanarFigure) { this->Modified(); m_PlanarFigure = planarFigure; } } mitk::Image::Pointer PlanarFigureMaskGenerator::GetReferenceImage() { if (IsUpdateRequired()) { this->CalculateMask(); } return m_ReferenceImage; } template < typename TPixel, unsigned int VImageDimension > void PlanarFigureMaskGenerator::InternalCalculateMaskFromPlanarFigure( const itk::Image< TPixel, VImageDimension > *image, unsigned int axis ) { - typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< unsigned short, 2 > MaskImage2DType; typename MaskImage2DType::Pointer maskImage = MaskImage2DType::New(); maskImage->SetOrigin(image->GetOrigin()); maskImage->SetSpacing(image->GetSpacing()); maskImage->SetLargestPossibleRegion(image->GetLargestPossibleRegion()); maskImage->SetBufferedRegion(image->GetBufferedRegion()); maskImage->SetDirection(image->GetDirection()); maskImage->SetNumberOfComponentsPerPixel(image->GetNumberOfComponentsPerPixel()); maskImage->Allocate(); maskImage->FillBuffer(1); // all PolylinePoints of the PlanarFigure are stored in a vtkPoints object. // These points are used by the vtkLassoStencilSource to create // a vtkImageStencil. const mitk::PlaneGeometry *planarFigurePlaneGeometry = m_PlanarFigure->GetPlaneGeometry(); const typename PlanarFigure::PolyLineType planarFigurePolyline = m_PlanarFigure->GetPolyLine( 0 ); const mitk::BaseGeometry *imageGeometry3D = m_inputImage->GetGeometry( 0 ); // If there is a second poly line in a closed planar figure, treat it as a hole. PlanarFigure::PolyLineType planarFigureHolePolyline; if (m_PlanarFigure->GetPolyLinesSize() == 2) planarFigureHolePolyline = m_PlanarFigure->GetPolyLine(1); // Determine x- and y-dimensions depending on principal axis // TODO use plane geometry normal to determine that automatically, then check whether the PF is aligned with one of the three principal axis int i0, i1; switch ( axis ) { case 0: i0 = 1; i1 = 2; break; case 1: i0 = 0; i1 = 2; break; case 2: default: i0 = 0; i1 = 1; break; } // store the polyline contour as vtkPoints object 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 // Fabian: From PlaneGeometry documentation: // Converts a 2D point given in mm (pt2d_mm) relative to the upper-left corner of the geometry into the corresponding world-coordinate (a 3D point in mm, pt3d_mm). // To convert a 2D point given in units (e.g., pixels in case of an image) into a 2D point given in mm (as required by this method), use IndexToWorld. planarFigurePlaneGeometry->Map( *it, 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 ); } vtkSmartPointer holePoints = nullptr; if (!planarFigureHolePolyline.empty()) { holePoints = vtkSmartPointer::New(); Point3D point3D; PlanarFigure::PolyLineType::const_iterator end = planarFigureHolePolyline.end(); for (it = planarFigureHolePolyline.begin(); it != end; ++it) { // Fabian: same as above planarFigurePlaneGeometry->Map(*it, point3D); imageGeometry3D->WorldToIndex(point3D, point3D); holePoints->InsertNextPoint(point3D[i0], point3D[i1], 0); } } // mark a malformed 2D planar figure ( i.e. area = 0 ) as out of bounds // this can happen when all control points of a rectangle lie on the same line = two of the three extents are zero double bounds[6] = {0, 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 ); vtkSmartPointer holeLassoStencil = nullptr; if (holePoints.GetPointer() != nullptr) { holeLassoStencil = vtkSmartPointer::New(); holeLassoStencil->SetShapeToPolygon(); holeLassoStencil->SetPoints(holePoints); } // Export from ITK to VTK (to use a VTK filter) typedef itk::VTKImageImport< MaskImage2DType > ImageImportType; typedef itk::VTKImageExport< MaskImage2DType > ImageExportType; typename ImageExportType::Pointer itkExporter = ImageExportType::New(); itkExporter->SetInput( maskImage ); // itkExporter->SetInput( castFilter->GetOutput() ); vtkSmartPointer vtkImporter = vtkSmartPointer::New(); this->ConnectPipelines( itkExporter, vtkImporter ); // Apply the generated image stencil to the input image vtkSmartPointer imageStencilFilter = vtkSmartPointer::New(); imageStencilFilter->SetInputConnection( vtkImporter->GetOutputPort() ); imageStencilFilter->SetStencilConnection(lassoStencil->GetOutputPort()); imageStencilFilter->ReverseStencilOff(); imageStencilFilter->SetBackgroundValue( 0 ); imageStencilFilter->Update(); vtkSmartPointer holeStencilFilter = nullptr; if (holeLassoStencil.GetPointer() != nullptr) { holeStencilFilter = vtkSmartPointer::New(); holeStencilFilter->SetInputConnection(imageStencilFilter->GetOutputPort()); holeStencilFilter->SetStencilConnection(holeLassoStencil->GetOutputPort()); holeStencilFilter->ReverseStencilOn(); holeStencilFilter->SetBackgroundValue(0); holeStencilFilter->Update(); } // Export from VTK back to ITK vtkSmartPointer vtkExporter = vtkSmartPointer::New(); vtkExporter->SetInputConnection( holeStencilFilter.GetPointer() == nullptr ? imageStencilFilter->GetOutputPort() : holeStencilFilter->GetOutputPort()); vtkExporter->Update(); typename ImageImportType::Pointer itkImporter = ImageImportType::New(); this->ConnectPipelines( vtkExporter, itkImporter ); itkImporter->Update(); typedef itk::ImageDuplicator< ImageImportType::OutputImageType > DuplicatorType; DuplicatorType::Pointer duplicator = DuplicatorType::New(); duplicator->SetInputImage( itkImporter->GetOutput() ); duplicator->Update(); // Store mask m_InternalITKImageMask2D = duplicator->GetOutput(); } bool PlanarFigureMaskGenerator::GetPrincipalAxis( const BaseGeometry *geometry, Vector3D vector, unsigned int &axis ) { vector.Normalize(); for ( unsigned int i = 0; i < 3; ++i ) { Vector3D axisVector = geometry->GetAxisVector( i ); axisVector.Normalize(); if ( fabs( fabs( axisVector * vector ) - 1.0) < mitk::eps ) { axis = i; return true; } } return false; } void PlanarFigureMaskGenerator::CalculateMask() { if (m_inputImage.IsNull()) { MITK_ERROR << "Image is not set."; } if (m_PlanarFigure.IsNull()) { MITK_ERROR << "PlanarFigure is not set."; } if (m_TimeStep != 0) { MITK_WARN << "Multiple TimeSteps are not supported in PlanarFigureMaskGenerator (yet)."; } const BaseGeometry *imageGeometry = m_inputImage->GetGeometry(); if ( imageGeometry == nullptr ) { throw std::runtime_error( "Image geometry invalid!" ); } if (m_inputImage->GetTimeSteps() > 0) { mitk::ImageTimeSelector::Pointer imgTimeSel = mitk::ImageTimeSelector::New(); imgTimeSel->SetInput(m_inputImage); imgTimeSel->SetTimeNr(m_TimeStep); imgTimeSel->UpdateLargestPossibleRegion(); m_InternalTimeSliceImage = imgTimeSel->GetOutput(); } else { m_InternalTimeSliceImage = m_inputImage; } m_InternalITKImageMask2D = nullptr; const PlaneGeometry *planarFigurePlaneGeometry = m_PlanarFigure->GetPlaneGeometry(); const PlaneGeometry *planarFigureGeometry = dynamic_cast< const PlaneGeometry * >( planarFigurePlaneGeometry ); //const BaseGeometry *imageGeometry = m_inputImage->GetGeometry(); // Find principal direction of PlanarFigure in input image unsigned int axis; if ( !this->GetPrincipalAxis( imageGeometry, planarFigureGeometry->GetNormal(), axis ) ) { throw std::runtime_error( "Non-aligned planar figures not supported!" ); } m_PlanarFigureAxis = axis; // Find slice number corresponding to PlanarFigure in input image itk::Image< unsigned short, 3 >::IndexType index; imageGeometry->WorldToIndex( planarFigureGeometry->GetOrigin(), index ); unsigned int slice = index[axis]; // extract image slice which corresponds to the planarFigure and store it in m_InternalImageSlice mitk::Image::Pointer inputImageSlice = extract2DImageSlice(axis, slice); //mitk::IOUtil::SaveImage(inputImageSlice, "/home/fabian/inputSliceImage.nrrd"); // Compute mask from PlanarFigure AccessFixedDimensionByItk_1(inputImageSlice, InternalCalculateMaskFromPlanarFigure, 2, axis) //convert itk mask to mitk::Image::Pointer and return it mitk::Image::Pointer planarFigureMaskImage; planarFigureMaskImage = mitk::GrabItkImageMemory(m_InternalITKImageMask2D); //mitk::IOUtil::SaveImage(planarFigureMaskImage, "/home/fabian/planarFigureMaskImage.nrrd"); //Convert2Dto3DImageFilter::Pointer sliceTo3DImageConverter = Convert2Dto3DImageFilter::New(); //sliceTo3DImageConverter->SetInput(planarFigureMaskImage); //sliceTo3DImageConverter->Update(); //mitk::IOUtil::SaveImage(sliceTo3DImageConverter->GetOutput(), "/home/fabian/3DsliceImage.nrrd"); m_ReferenceImage = inputImageSlice; //mitk::IOUtil::SaveImage(m_ReferenceImage, "/home/fabian/referenceImage.nrrd"); m_InternalMask = planarFigureMaskImage; } void PlanarFigureMaskGenerator::SetTimeStep(unsigned int timeStep) { if (timeStep != m_TimeStep) { m_TimeStep = timeStep; } } mitk::Image::Pointer PlanarFigureMaskGenerator::GetMask() { if (IsUpdateRequired()) { this->CalculateMask(); this->Modified(); } m_InternalMaskUpdateTime = this->GetMTime(); return m_InternalMask; } mitk::Image::Pointer PlanarFigureMaskGenerator::extract2DImageSlice(unsigned int axis, unsigned int slice) { // Extract slice with given position and direction from image unsigned int dimension = m_InternalTimeSliceImage->GetDimension(); mitk::Image::Pointer imageSlice = mitk::Image::New(); if (dimension == 3) { ExtractImageFilter::Pointer imageExtractor = ExtractImageFilter::New(); imageExtractor->SetInput( m_InternalTimeSliceImage ); imageExtractor->SetSliceDimension( axis ); imageExtractor->SetSliceIndex( slice ); imageExtractor->Update(); imageSlice = imageExtractor->GetOutput(); } else if(dimension == 2) { imageSlice = m_InternalTimeSliceImage; } else { MITK_ERROR << "Unsupported image dimension. Dimension is: " << dimension << ". Only 2D and 3D images are supported."; } return imageSlice; } bool PlanarFigureMaskGenerator::IsUpdateRequired() const { unsigned long thisClassTimeStamp = this->GetMTime(); unsigned long internalMaskTimeStamp = m_InternalMask->GetMTime(); unsigned long planarFigureTimeStamp = m_PlanarFigure->GetMTime(); unsigned long inputImageTimeStamp = m_inputImage->GetMTime(); if (thisClassTimeStamp > m_InternalMaskUpdateTime) // inputs have changed { return true; } if (m_InternalMaskUpdateTime < planarFigureTimeStamp || m_InternalMaskUpdateTime < inputImageTimeStamp) // mask image has changed outside of this class { return true; } if (internalMaskTimeStamp > m_InternalMaskUpdateTime) // internal mask has been changed outside of this class { return true; } return false; } }