diff --git a/Modules/ImageStatistics/Testing/CMakeLists.txt b/Modules/ImageStatistics/Testing/CMakeLists.txt index 153cd81e2e..012a6a7a63 100644 --- a/Modules/ImageStatistics/Testing/CMakeLists.txt +++ b/Modules/ImageStatistics/Testing/CMakeLists.txt @@ -1 +1,3 @@ MITK_CREATE_MODULE_TESTS() + +mitkAddCustomModuleTest(mitkImageStatisticsHotspotTest_Case1 mitkImageStatisticsHotspotTest ${CMAKE_CURRENT_SOURCE_DIR}/Data/Hotspot_Case1.xml) diff --git a/Modules/ImageStatistics/Testing/Data/Hotspot_Case1.xml b/Modules/ImageStatistics/Testing/Data/Hotspot_Case1.xml new file mode 100644 index 0000000000..eebe8211ed --- /dev/null +++ b/Modules/ImageStatistics/Testing/Data/Hotspot_Case1.xml @@ -0,0 +1,19 @@ + + + + + + + + + + + + + + + + + + + diff --git a/Modules/ImageStatistics/Testing/Data/Hotspot_Case10.xml b/Modules/ImageStatistics/Testing/Data/Hotspot_Case10.xml new file mode 100644 index 0000000000..cad93554d2 --- /dev/null +++ b/Modules/ImageStatistics/Testing/Data/Hotspot_Case10.xml @@ -0,0 +1,19 @@ + + + + + + + + + + + + + + + + + + + diff --git a/Modules/ImageStatistics/Testing/Data/Hotspot_Case11.xml b/Modules/ImageStatistics/Testing/Data/Hotspot_Case11.xml new file mode 100644 index 0000000000..36c0fdc14a --- /dev/null +++ b/Modules/ImageStatistics/Testing/Data/Hotspot_Case11.xml @@ -0,0 +1,19 @@ + + + + + + + + + + + + + + + + + + + diff --git a/Modules/ImageStatistics/Testing/Data/Hotspot_Case12.xml b/Modules/ImageStatistics/Testing/Data/Hotspot_Case12.xml new file mode 100644 index 0000000000..96ee21f963 --- /dev/null +++ b/Modules/ImageStatistics/Testing/Data/Hotspot_Case12.xml @@ -0,0 +1,19 @@ + + + + + + + + + + + + + + + + + + + diff --git a/Modules/ImageStatistics/Testing/Data/Hotspot_Case13.xml b/Modules/ImageStatistics/Testing/Data/Hotspot_Case13.xml new file mode 100644 index 0000000000..c170b9e7a2 --- /dev/null +++ b/Modules/ImageStatistics/Testing/Data/Hotspot_Case13.xml @@ -0,0 +1,19 @@ + + + + + + + + + + + + + + + + + + + diff --git a/Modules/ImageStatistics/Testing/Data/Hotspot_Case14.xml b/Modules/ImageStatistics/Testing/Data/Hotspot_Case14.xml new file mode 100644 index 0000000000..f1323b4a24 --- /dev/null +++ b/Modules/ImageStatistics/Testing/Data/Hotspot_Case14.xml @@ -0,0 +1,19 @@ + + + + + + + + + + + + + + + + + + + diff --git a/Modules/ImageStatistics/Testing/Data/Hotspot_Case15.xml b/Modules/ImageStatistics/Testing/Data/Hotspot_Case15.xml new file mode 100644 index 0000000000..03375e359a --- /dev/null +++ b/Modules/ImageStatistics/Testing/Data/Hotspot_Case15.xml @@ -0,0 +1,19 @@ + + + + + + + + + + + + + + + + + + + diff --git a/Modules/ImageStatistics/Testing/Data/Hotspot_Case16.xml b/Modules/ImageStatistics/Testing/Data/Hotspot_Case16.xml new file mode 100644 index 0000000000..86d47ce303 --- /dev/null +++ b/Modules/ImageStatistics/Testing/Data/Hotspot_Case16.xml @@ -0,0 +1,19 @@ + + + + + + + + + + + + + + + + + + + diff --git a/Modules/ImageStatistics/Testing/Data/Hotspot_Case17.xml b/Modules/ImageStatistics/Testing/Data/Hotspot_Case17.xml new file mode 100644 index 0000000000..c4c82b8a5a --- /dev/null +++ b/Modules/ImageStatistics/Testing/Data/Hotspot_Case17.xml @@ -0,0 +1,19 @@ + + + + + + + + + + + + + + + + + + + diff --git a/Modules/ImageStatistics/Testing/Data/Hotspot_Case18.xml b/Modules/ImageStatistics/Testing/Data/Hotspot_Case18.xml new file mode 100644 index 0000000000..5beb23c2a0 --- /dev/null +++ b/Modules/ImageStatistics/Testing/Data/Hotspot_Case18.xml @@ -0,0 +1,19 @@ + + + + + + + + + + + + + + + + + + + diff --git a/Modules/ImageStatistics/Testing/Data/Hotspot_Case19.xml b/Modules/ImageStatistics/Testing/Data/Hotspot_Case19.xml new file mode 100644 index 0000000000..058cfd8731 --- /dev/null +++ b/Modules/ImageStatistics/Testing/Data/Hotspot_Case19.xml @@ -0,0 +1,19 @@ + + + + + + + + + + + + + + + + + + + diff --git a/Modules/ImageStatistics/Testing/Data/Hotspot_Case2.xml b/Modules/ImageStatistics/Testing/Data/Hotspot_Case2.xml new file mode 100644 index 0000000000..823baa7bb8 --- /dev/null +++ b/Modules/ImageStatistics/Testing/Data/Hotspot_Case2.xml @@ -0,0 +1,19 @@ + + + + + + + + + + + + + + + + + + + diff --git a/Modules/ImageStatistics/Testing/Data/Hotspot_Case20.xml b/Modules/ImageStatistics/Testing/Data/Hotspot_Case20.xml new file mode 100644 index 0000000000..81308b6585 --- /dev/null +++ b/Modules/ImageStatistics/Testing/Data/Hotspot_Case20.xml @@ -0,0 +1,19 @@ + + + + + + + + + + + + + + + + + + + diff --git a/Modules/ImageStatistics/Testing/Data/Hotspot_Case21.xml b/Modules/ImageStatistics/Testing/Data/Hotspot_Case21.xml new file mode 100644 index 0000000000..ef484fb0dd --- /dev/null +++ b/Modules/ImageStatistics/Testing/Data/Hotspot_Case21.xml @@ -0,0 +1,19 @@ + + + + + + + + + + + + + + + + + + + diff --git a/Modules/ImageStatistics/Testing/Data/Hotspot_Case3.xml b/Modules/ImageStatistics/Testing/Data/Hotspot_Case3.xml new file mode 100644 index 0000000000..033da8710a --- /dev/null +++ b/Modules/ImageStatistics/Testing/Data/Hotspot_Case3.xml @@ -0,0 +1,19 @@ + + + + + + + + + + + + + + + + + + + diff --git a/Modules/ImageStatistics/Testing/Data/Hotspot_Case4.xml b/Modules/ImageStatistics/Testing/Data/Hotspot_Case4.xml new file mode 100644 index 0000000000..96180b2922 --- /dev/null +++ b/Modules/ImageStatistics/Testing/Data/Hotspot_Case4.xml @@ -0,0 +1,19 @@ + + + + + + + + + + + + + + + + + + + diff --git a/Modules/ImageStatistics/Testing/Data/Hotspot_Case5.xml b/Modules/ImageStatistics/Testing/Data/Hotspot_Case5.xml new file mode 100644 index 0000000000..b88b176040 --- /dev/null +++ b/Modules/ImageStatistics/Testing/Data/Hotspot_Case5.xml @@ -0,0 +1,19 @@ + + + + + + + + + + + + + + + + + + + diff --git a/Modules/ImageStatistics/Testing/Data/Hotspot_Case6.xml b/Modules/ImageStatistics/Testing/Data/Hotspot_Case6.xml new file mode 100644 index 0000000000..1948c59ba4 --- /dev/null +++ b/Modules/ImageStatistics/Testing/Data/Hotspot_Case6.xml @@ -0,0 +1,19 @@ + + + + + + + + + + + + + + + + + + + diff --git a/Modules/ImageStatistics/Testing/Data/Hotspot_Case7.xml b/Modules/ImageStatistics/Testing/Data/Hotspot_Case7.xml new file mode 100644 index 0000000000..1addb3b5dc --- /dev/null +++ b/Modules/ImageStatistics/Testing/Data/Hotspot_Case7.xml @@ -0,0 +1,19 @@ + + + + + + + + + + + + + + + + + + + diff --git a/Modules/ImageStatistics/Testing/Data/Hotspot_Case8.xml b/Modules/ImageStatistics/Testing/Data/Hotspot_Case8.xml new file mode 100644 index 0000000000..5a8f17e349 --- /dev/null +++ b/Modules/ImageStatistics/Testing/Data/Hotspot_Case8.xml @@ -0,0 +1,19 @@ + + + + + + + + + + + + + + + + + + + diff --git a/Modules/ImageStatistics/Testing/Data/Hotspot_Case9.xml b/Modules/ImageStatistics/Testing/Data/Hotspot_Case9.xml new file mode 100644 index 0000000000..f786bffda3 --- /dev/null +++ b/Modules/ImageStatistics/Testing/Data/Hotspot_Case9.xml @@ -0,0 +1,19 @@ + + + + + + + + + + + + + + + + + + + diff --git a/Modules/ImageStatistics/Testing/Data/neue Bilder/Hotspot_Case1.xml b/Modules/ImageStatistics/Testing/Data/neue Bilder/Hotspot_Case1.xml new file mode 100644 index 0000000000..e57655f2d1 --- /dev/null +++ b/Modules/ImageStatistics/Testing/Data/neue Bilder/Hotspot_Case1.xml @@ -0,0 +1,7 @@ + + + + + + + diff --git a/Modules/ImageStatistics/Testing/Data/neue Bilder/Hotspot_Case2.xml b/Modules/ImageStatistics/Testing/Data/neue Bilder/Hotspot_Case2.xml new file mode 100644 index 0000000000..088d464625 --- /dev/null +++ b/Modules/ImageStatistics/Testing/Data/neue Bilder/Hotspot_Case2.xml @@ -0,0 +1,7 @@ + + + + + + + diff --git a/Modules/ImageStatistics/Testing/Data/neue Bilder/Hotspot_Case3.xml b/Modules/ImageStatistics/Testing/Data/neue Bilder/Hotspot_Case3.xml new file mode 100644 index 0000000000..8bd11ad0c4 --- /dev/null +++ b/Modules/ImageStatistics/Testing/Data/neue Bilder/Hotspot_Case3.xml @@ -0,0 +1,7 @@ + + + + + + + diff --git a/Modules/ImageStatistics/Testing/Data/neue Bilder/Hotspot_Case4.xml b/Modules/ImageStatistics/Testing/Data/neue Bilder/Hotspot_Case4.xml new file mode 100644 index 0000000000..438ae152c1 --- /dev/null +++ b/Modules/ImageStatistics/Testing/Data/neue Bilder/Hotspot_Case4.xml @@ -0,0 +1,7 @@ + + + + + + + diff --git a/Modules/ImageStatistics/Testing/Data/neue Bilder/Hotspot_Case5.xml b/Modules/ImageStatistics/Testing/Data/neue Bilder/Hotspot_Case5.xml new file mode 100644 index 0000000000..3657c60064 --- /dev/null +++ b/Modules/ImageStatistics/Testing/Data/neue Bilder/Hotspot_Case5.xml @@ -0,0 +1,7 @@ + + + + + + + diff --git a/Modules/ImageStatistics/Testing/Data/neue Bilder/Hotspot_Case6.xml b/Modules/ImageStatistics/Testing/Data/neue Bilder/Hotspot_Case6.xml new file mode 100644 index 0000000000..ecc58a523b --- /dev/null +++ b/Modules/ImageStatistics/Testing/Data/neue Bilder/Hotspot_Case6.xml @@ -0,0 +1,7 @@ + + + + + + + diff --git a/Modules/ImageStatistics/Testing/Data/neue Bilder/Hotspot_Case7.xml b/Modules/ImageStatistics/Testing/Data/neue Bilder/Hotspot_Case7.xml new file mode 100644 index 0000000000..0ac9a6c591 --- /dev/null +++ b/Modules/ImageStatistics/Testing/Data/neue Bilder/Hotspot_Case7.xml @@ -0,0 +1,7 @@ + + + + + + + diff --git a/Modules/ImageStatistics/Testing/Data/neue Bilder/Hotspot_Case8.xml b/Modules/ImageStatistics/Testing/Data/neue Bilder/Hotspot_Case8.xml new file mode 100644 index 0000000000..37d5cf1cec --- /dev/null +++ b/Modules/ImageStatistics/Testing/Data/neue Bilder/Hotspot_Case8.xml @@ -0,0 +1,7 @@ + + + + + + + diff --git a/Modules/ImageStatistics/Testing/Data/neue Bilder/Hotspot_Case9.xml b/Modules/ImageStatistics/Testing/Data/neue Bilder/Hotspot_Case9.xml new file mode 100644 index 0000000000..f234557dc1 --- /dev/null +++ b/Modules/ImageStatistics/Testing/Data/neue Bilder/Hotspot_Case9.xml @@ -0,0 +1,7 @@ + + + + + + + diff --git a/Modules/ImageStatistics/Testing/files.cmake b/Modules/ImageStatistics/Testing/files.cmake index 01792788aa..5fb127cbdf 100644 --- a/Modules/ImageStatistics/Testing/files.cmake +++ b/Modules/ImageStatistics/Testing/files.cmake @@ -1,5 +1,10 @@ set(MODULE_TESTS mitkImageStatisticsCalculatorTest.cpp mitkPointSetStatisticsCalculatorTest.cpp mitkPointSetDifferenceStatisticsCalculatorTest.cpp + mitkMultiGaussianTest.cpp +) + +set(MODULE_CUSTOM_TESTS + mitkImageStatisticsHotspotTest.cpp ) diff --git a/Modules/ImageStatistics/Testing/itkMultiGaussianImageSource.h b/Modules/ImageStatistics/Testing/itkMultiGaussianImageSource.h new file mode 100644 index 0000000000..98a236789c --- /dev/null +++ b/Modules/ImageStatistics/Testing/itkMultiGaussianImageSource.h @@ -0,0 +1,200 @@ +/*========================================================================= + * + * 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_h +#define __itkMultiGaussianImageSource_h + +#include "itkImageSource.h" +#include "itkNumericTraits.h" +#include "itkImageRegionIteratorWithIndex.h" +#include "itkImageFileWriter.h" + +namespace itk +{ +/** \class MultiGaussianImageSource + * \brief Generate an 3-dimensional multigaussian image. + * This class defines an 3-dimensional Image, in which the value at one voxel equals the value of a multigaussian function evaluated at the voxel's * coordinates. The multigaussian function is build as a sum of N gaussian function. This is defined by the following parameters: + * 1. CenterX, CenterY, CenterZ - vectors of the size of N determining the expectancy value at the x-, y- and the z-axis. That means: The i-th + * gaussian bell curve takes its maximal value at the voxel with index [CenterX(i); CenterY(i); Centerz(i)]. + * 2. SigmaX, SigmaY, SigmaZ - vectors of the size of N determining the deviation at the x-, y- and the z-axis. That means: The width of the i-th + * gaussian bell curve in the x-axis is SigmaX(i), in the y-axis is SigmaY(i) and in the z-axis is SigmaZ(i). + * 3. Altitude - vector of the size of N determining the altitude: the i-th gaussian bell curve has a height of Altitude(i). + * This class allows by the method CalculateMidpointAndMeanValue() 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. + * + * \ingroup DataSources + * \ingroup ITKTestKernel + * + * \wiki + * \wikiexample{} + * \endwiki + */ +template< typename TOutputImage > +class ITK_EXPORT MultiGaussianImageSource:public ImageSource< TOutputImage > +{ +public: + /** Standard class typedefs. */ + typedef MultiGaussianImageSource Self; + typedef ImageSource< TOutputImage > Superclass; + typedef SmartPointer< Self > Pointer; + typedef SmartPointer< const Self > ConstPointer; + + /** Typedef for the output image PixelType. */ + typedef typename TOutputImage::PixelType OutputImagePixelType; + + /** Typedef to describe the output image region type. */ + typedef typename TOutputImage::RegionType OutputImageRegionType; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Basic types from the OutputImageType */ + typedef typename TOutputImage::SizeType SizeType; + typedef typename TOutputImage::IndexType IndexType; + typedef typename TOutputImage::SpacingType SpacingType; + typedef typename TOutputImage::PointType PointType; + typedef typename SizeType::SizeValueType SizeValueType; + typedef SizeValueType SizeValueArrayType[TOutputImage::ImageDimension]; + typedef typename TOutputImage::SpacingValueType SpacingValueType; + typedef SpacingValueType SpacingValueArrayType[TOutputImage::ImageDimension]; + typedef typename TOutputImage::PointValueType PointValueType; + typedef PointValueType PointValueArrayType[TOutputImage::ImageDimension]; + /** Typedef to describe the sphere radius type. */ + typedef double RadiusType; + /** Typedef to describe the standard vector type. */ + typedef std::vector VectorType; + /** Typedef to describe the itk vector type. */ + typedef itk::Vector ItkVectorType; + /** Typedef to describe the ImageRegionIteratorWithIndex type. */ + typedef ImageRegionIteratorWithIndex IteratorType; + /** Typedef to describe the Poiner type at the output image. */ + typedef typename TOutputImage::Pointer ImageType; + + /** Set/Get size of the output image */ + itkSetMacro(Size, SizeType); + virtual void SetSize(SizeValueArrayType sizeArray); + + virtual const SizeValueType * GetSize() const; + + /** Set/Get spacing of the output image */ + itkSetMacro(Spacing, SpacingType); + virtual void SetSpacing(SpacingValueArrayType spacingArray); + virtual const SpacingValueType * GetSpacing() const; + + /** Set/Get origin of the output image */ + itkSetMacro(Origin, PointType); + virtual void SetOrigin(PointValueArrayType originArray); + virtual const PointValueType * GetOrigin() const; + + + /** Get the number of gaussian functions in the output image */ + virtual unsigned int GetNumberOfGaussians() const; + + virtual void SetNumberOfGausssians( unsigned int ); + /** Set/Get the radius of the sphere */ + virtual const RadiusType GetRadius() const; + virtual void SetRadius( RadiusType radius ); + /** Set/Get the number of steps to traverse the radius of the sphere */ + virtual const int GetRadiusStepNumber() const; + virtual void SetRadiusStepNumber( unsigned int stepNumber ); + /** Get the maximal mean value in a sphere over all possible spheres with midpoint in the image */ + virtual const double GetMaxMeanValue() const; + /** Get the index of the midpoint of a sphere with the maximal mean value */ + virtual const ItkVectorType GetSphereMidpoint() const; + /** Calculates the value of the multigaussian function at a Point given by its coordinates [x;y;z] */ + virtual const double MultiGaussianFunctionValueAtPoint(double , double, double); + /** Adds a multigaussian defined by the parameter: CenterX, CenterY, CenterZ, SigmaX, SigmaY, SigmaZ, Altitude. + All parameters should have the same size, which determinates the number of the gaussian added. */ + virtual void AddGaussian( VectorType, VectorType, VectorType, VectorType, VectorType, VectorType, VectorType); + /** Calculates and set the index of the midpoint of the sphere with the maximal mean value as well as the mean value*/ + virtual void CalculateMidpointAndMeanValue(); + + /** Set the minimum possible pixel value. By default, it is + * NumericTraits::min(). */ + itkSetClampMacro( Min, OutputImagePixelType, + NumericTraits< OutputImagePixelType >::NonpositiveMin(), + NumericTraits< OutputImagePixelType >::max() ); + + /** Get the minimum possible pixel value. */ + itkGetConstMacro(Min, OutputImagePixelType); + + /** Set the maximum possible pixel value. By default, it is + * NumericTraits::max(). */ + itkSetClampMacro( Max, OutputImagePixelType, + NumericTraits< OutputImagePixelType >::NonpositiveMin(), + NumericTraits< OutputImagePixelType >::max() ); + + /** Get the maximum possible pixel value. */ + itkGetConstMacro(Max, OutputImagePixelType); + + + + +protected: + MultiGaussianImageSource(); + ~MultiGaussianImageSource(); + void PrintSelf(std::ostream & os, Indent indent) const; + + virtual void GenerateData(); + virtual void GenerateOutputInformation(); + +private: + MultiGaussianImageSource(const MultiGaussianImageSource &); //purposely not implemented + void operator=(const MultiGaussianImageSource &); //purposely not implemented + + SizeType m_Size; //size of the output image + SpacingType m_Spacing; //spacing + PointType m_Origin; //origin + + + unsigned int m_NumberOfGaussians; //number of Gaussians + RadiusType m_Radius; //radius of the sphere + unsigned int m_RadiusStepNumber; //number of steps to traverse the sphere radius + OutputImagePixelType m_MeanValue; //mean value in the wanted sphere + ItkVectorType m_SphereMidpoint; //midpoint of the wanted sphere + VectorType m_SigmaX; //deviation in the x-axis + VectorType m_SigmaY; //deviation in the y-axis + VectorType m_SigmaZ; //deviation in the z-axis + VectorType m_CenterX; //x-coordinate of the mean value of Gaussians + VectorType m_CenterY; //y-coordinate of the mean value of Gaussians + VectorType m_CenterZ; //z-coordinate of the mean value of Gaussians + VectorType m_Altitude; //amplitude + + typename TOutputImage::PixelType m_Min; //minimum possible value + typename TOutputImage::PixelType m_Max; //maximum possible value + + // The following variables are deprecated, and provided here just for + // backward compatibility. It use is discouraged. + mutable PointValueArrayType m_OriginArray; + mutable SpacingValueArrayType m_SpacingArray; +}; +} // end namespace itk + +#ifndef ITK_MANUAL_INSTANTIATION +#include "itkMultiGaussianImageSource.hxx" +#endif + +#endif diff --git a/Modules/ImageStatistics/Testing/itkMultiGaussianImageSource.hxx b/Modules/ImageStatistics/Testing/itkMultiGaussianImageSource.hxx new file mode 100644 index 0000000000..f39a2b7122 --- /dev/null +++ b/Modules/ImageStatistics/Testing/itkMultiGaussianImageSource.hxx @@ -0,0 +1,446 @@ +/*========================================================================= + * + * 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 "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 > +const 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 > +const int +MultiGaussianImageSource< TOutputImage > +::GetRadiusStepNumber() const +{ + return this->m_RadiusStepNumber; +} + +template< class TOutputImage > +void +MultiGaussianImageSource< TOutputImage > +::SetRadiusStepNumber( unsigned int stepNumber ) +{ + this->m_RadiusStepNumber = stepNumber; +} + +template< class TOutputImage > +void +MultiGaussianImageSource< TOutputImage > +::SetNumberOfGausssians( unsigned int n ) +{ + this->m_NumberOfGaussians = n; +} + +template< class TOutputImage > +const double +MultiGaussianImageSource< TOutputImage > +::GetMaxMeanValue() const +{ + return m_MeanValue; +} + + +//---------------------------------------------------------------------------- +template< class TOutputImage > +const typename MultiGaussianImageSource< TOutputImage >::ItkVectorType +MultiGaussianImageSource< TOutputImage > +::GetSphereMidpoint() const +{ + return m_SphereMidpoint; +} + +template< class TOutputImage > +const double +MultiGaussianImageSource< TOutputImage > +::MultiGaussianFunctionValueAtPoint(double x, double y, double z) +{ + double summand0, summand1, summand2, power, value = 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); + } + 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; + + 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 class allows by the method CalculateMidpointAndMeanValue() 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. The parameter RadiusStepNumber controls the accuracy of that calculation (the higher the value the higher the exactness). + The algorithm works as follows: + 1. the first three for-loops traverse the image and assume the current point to be the wanted sphere midpoint + 2. calculate the mean value for that sphere (use sphere coordinates): + 2.1. traverse the radius of the sphere with step size Radius divided by RadiusStepNumber (the for-loop with index i) + 2.2. define a variable dist, which gives a approximately distance between the points at the sphere surface + (here we take such a distance, that on the smallest sphere are located 8 points) + 2.3. calculate the angles so that the points are equally spaced on the surface (for-loops with indexes j and k) + 2.3. for all radius length add the values at the points on the sphere and divide by the number of points added + (the values at each point is calculate with the method MultiGaussianFunctionValueAtPoint()) + 3. Compare with the until-now-found-maximum and take the bigger one +*/ +template< typename TOutputImage > +void +MultiGaussianImageSource< TOutputImage > +::CalculateMidpointAndMeanValue() +{ + itkDebugMacro(<< "Generating a image of scalars "); + int numberSummand = 0, angleStepNumberOverTwo; + double meanValueTemp, valueAdd, value, x, y, z, temp; + double riStep, fijStep, psikStep, ri, fij, psik; + double dist = itk::Math::pi * m_Radius / (2 * m_RadiusStepNumber); + m_MeanValue = 0; + riStep = m_Radius / m_RadiusStepNumber; + for(unsigned int index0 = 0; index0 < m_Size[0]; ++index0) + { + for(unsigned int index1 = 0; index1 < m_Size[1]; ++index1) + { + for(unsigned int index2 = 0; index2 < m_Size[2]; ++index2) + { + numberSummand = 0; + value = 0.0; + ri = riStep; + for(unsigned int i = 0; i < m_RadiusStepNumber; ++i) + { + angleStepNumberOverTwo = static_cast( itk::Math::pi * ri / dist); + fij = 0; + fijStep = itk::Math::pi / angleStepNumberOverTwo; + for(unsigned int j = 0; j <= angleStepNumberOverTwo; ++j) // from 0 to pi + { + z = ri * cos(fij); + psikStep = 2.0 * itk::Math::pi / (2.0 * angleStepNumberOverTwo); + psik = -itk::Math::pi + psikStep; + + temp = ri * sin(fij); + for(unsigned int k = 0; k < 2 * angleStepNumberOverTwo; ++k) // from -pi to pi + { + x = temp * cos(psik); + y = temp * sin(psik); + numberSummand++; + valueAdd = MultiGaussianFunctionValueAtPoint(x + index0, y + index1, z + index2); + value = value + valueAdd; + psik = psik + psikStep; + } + fij = fij + fijStep; + } + ri = ri + riStep; + } + + meanValueTemp = value / numberSummand; + if(meanValueTemp > m_MeanValue) + { + m_MeanValue = meanValueTemp; + m_SphereMidpoint.SetElement( 0, index0 ); + m_SphereMidpoint.SetElement( 1, index1 ); + m_SphereMidpoint.SetElement( 2, index2 ); + } + } + } + } +} +} // end namespace itk +#endif diff --git a/Modules/ImageStatistics/Testing/mitkImageStatisticsHotspotTest.cpp b/Modules/ImageStatistics/Testing/mitkImageStatisticsHotspotTest.cpp new file mode 100644 index 0000000000..5d6b277e8a --- /dev/null +++ b/Modules/ImageStatistics/Testing/mitkImageStatisticsHotspotTest.cpp @@ -0,0 +1,544 @@ +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center, +Division of Medical and Biological Informatics. +All rights reserved. + +This software is distributed WITHOUT ANY WARRANTY; without +even the implied warranty of MERCHANTABILITY or FITNESS FOR +A PARTICULAR PURPOSE. + +See LICENSE.txt or http://www.mitk.org for details. + +===================================================================*/ + +#include "mitkImageStatisticsCalculator.h" +#include "itkMultiGaussianImageSource.h" +#include "mitkTestingMacros.h" + +#include + +#include + + +#include +#include + +struct mitkImageStatisticsHotspotTestClass +{ + + /** + \brief Test parameters for one test case. + + Describes all aspects of a single test case: + - parameters to generate a test image + - parameters of a ROI that describes where to calculate statistics + - expected statistics results + */ + struct Parameters + { + // XML-Tag + int m_ImageRows; // XML-Tag "image-rows" + int m_ImageColumns; // XML-Tag "image-columns" + int m_ImageSlices; // XML-Tag "image-slices" + + int m_NumberOfGaussian; // XML-Tag "numberOfGaussians" + + std::vector m_Spacing; // XML-Tag "spacingX", "spacingY", "spacingZ" + + // XML-Tag + std::vector m_CenterX; // XML-Tag "centerIndexX" + std::vector m_CenterY; // XML-Tag "centerIndexY" + std::vector m_CenterZ; // XML-Tag "centerIndexZ" + + std::vector m_SigmaX; // XML-Tag "deviationX" + std::vector m_SigmaY; // XML-Tag "deviationY" + std::vector m_SigmaZ; // XML-Tag "deviationZ" + + std::vector m_Altitude; // XML-Tag "altitude" + + // XML-Tag + float m_HotspotMinimum; // XML-Tag "minimum" + float m_HotspotMaximum; // XML-Tag "maximum" + float m_HotspotPeak; // XML-Tag "peak" + + std::vector m_HotspotMaximumIndex; // XML-Tag "maximumIndexX", XML-Tag "maximumIndexY", XML-Tag "maximumIndexZ" + std::vector m_HotspotMinimumIndex; // XML-Tag "minimumIndexX", XML-Tag "minimumIndexY", XML-Tag "minimumIndexZ" + std::vector m_HotspotPeakIndex; // XML-Tag "peakIndexX", XML-Tag "peakIndexY", XML-Tag "peakIndexZ" + }; + + /** + \brief Find/Convert integer attribute in itk::DOMNode. + */ + static int GetIntegerAttribute(itk::DOMNode* domNode, const std::string& tag) + { + assert(domNode); + MITK_TEST_CONDITION_REQUIRED( domNode->HasAttribute(tag), "Tag '" << tag << "' is defined in test parameters" ); + std::string attributeValue = domNode->GetAttribute(tag); + + int resultValue; + try + { + //MITK_TEST_OUTPUT( << "Converting tag value '" << attributeValue << "' for tag '" << tag << "' to integer"); + std::stringstream(attributeValue) >> resultValue; + return resultValue; + } + catch(std::exception& e) + { + MITK_TEST_CONDITION_REQUIRED(false, "Convert tag value '" << attributeValue << "' for tag '" << tag << "' to integer"); + return 0; // just to satisfy compiler + } + } + /** + \brief Find/Convert double attribute in itk::DOMNode. + */ + static double GetDoubleAttribute(itk::DOMNode* domNode, const std::string& tag) + { + assert(domNode); + MITK_TEST_CONDITION_REQUIRED( domNode->HasAttribute(tag), "Tag '" << tag << "' is defined in test parameters" ); + std::string attributeValue = domNode->GetAttribute(tag); + + double resultValue; + try + { + //MITK_TEST_OUTPUT( << "Converting tag value '" << attributeValue << "' for tag '" << tag << "' to double"); + std::stringstream(attributeValue) >> resultValue; + return resultValue; + } + catch(std::exception& e) + { + MITK_TEST_CONDITION_REQUIRED(false, "Convert tag value '" << attributeValue << "' for tag '" << tag << "' to double"); + return 0.0; // just to satisfy compiler + } + } + + + /** + \brief Read XML file describging the test parameters. + + Reads XML file given in first commandline parameter in order + to construct a Parameters structure. The XML file should be + structurs as the following example, i.e. we describe the + three test aspects of Parameters in three different tags, + with all the details described as tag attributes. + +\verbatim + + + + + + + + + +\verbatim + + The different parameters are interpreted as follows: + ... TODO TODO TODO ... + + */ + static Parameters ParseParameters(int argc, char* argv[]) + { + // - parse parameters + // - fill ALL values of result structure + // - if necessary, provide c'tor and default parameters to Parameters + + MITK_TEST_CONDITION_REQUIRED(argc == 2, "Test is invoked with exactly 1 parameter (XML parameters file)"); + MITK_INFO << "Reading parameters from file '" << argv[1] << "'"; + std::string filename = argv[1]; + + Parameters result; + + itk::DOMNodeXMLReader::Pointer xmlReader = itk::DOMNodeXMLReader::New(); + xmlReader->SetFileName( filename ); + try + { + xmlReader->Update(); + itk::DOMNode::Pointer domRoot = xmlReader->GetOutput(); + typedef std::vector NodeList; + // read test image parameters, fill result structure + NodeList testimages; + domRoot->GetChildren("testimage", testimages); + MITK_TEST_CONDITION_REQUIRED( testimages.size() == 1, "One test image defined" ) + itk::DOMNode* testimage = testimages[0]; + + result.m_ImageRows = GetIntegerAttribute( testimage, "image-rows" ); + result.m_ImageColumns = GetIntegerAttribute( testimage, "image-columns" ); + result.m_ImageSlices = GetIntegerAttribute( testimage, "image-slices" ); + + result.m_NumberOfGaussian = GetIntegerAttribute( testimage, "numberOfGaussians" ); + + std::vector tmpSpacing(3,0); + + tmpSpacing[0] = GetDoubleAttribute(testimage, "spacingX"); + tmpSpacing[1] = GetDoubleAttribute(testimage, "spacingY"); + tmpSpacing[2] = GetDoubleAttribute(testimage, "spacingZ"); + + result.m_Spacing = tmpSpacing; + + MITK_TEST_OUTPUT( << "Read size parameters (x,y,z): " << result.m_ImageRows << "," << result.m_ImageColumns << "," << result.m_ImageSlices); + MITK_TEST_OUTPUT( << "Read spacing parameters (x,y,z): " << result.m_Spacing[0] << "," << result.m_Spacing[1] << "," << result.m_Spacing[2]); + + NodeList gaussians; + testimage->GetChildren("gaussian", gaussians); + MITK_TEST_CONDITION_REQUIRED( gaussians.size() >= 1, "At least one gaussian is defined" ) + + std::vector tmpCenterX(result.m_NumberOfGaussian,0); + std::vector tmpCenterY(result.m_NumberOfGaussian,0); + std::vector tmpCenterZ(result.m_NumberOfGaussian,0); + + std::vector tmpSigmaX(result.m_NumberOfGaussian,0); + std::vector tmpSigmaY(result.m_NumberOfGaussian,0); + std::vector tmpSigmaZ(result.m_NumberOfGaussian,0); + + std::vector tmpAltitude(result.m_NumberOfGaussian,0); + + for(int i = 0; i < result.m_NumberOfGaussian ; ++i) + { + itk::DOMNode* gaussian = gaussians[i]; + + tmpCenterX[i] = GetIntegerAttribute(gaussian, "centerIndexX"); + tmpCenterY[i] = GetIntegerAttribute(gaussian, "centerIndexY"); + tmpCenterZ[i] = GetIntegerAttribute(gaussian, "centerIndexZ"); + + tmpSigmaX[i] = GetIntegerAttribute(gaussian, "deviationX"); + tmpSigmaY[i] = GetIntegerAttribute(gaussian, "deviationY"); + tmpSigmaZ[i] = GetIntegerAttribute(gaussian, "deviationZ"); + + tmpAltitude[i] = GetIntegerAttribute(gaussian, "altitude"); + } + + result.m_CenterX = tmpCenterX; + result.m_CenterY = tmpCenterY; + result.m_CenterZ = tmpCenterZ; + + result.m_SigmaX = tmpSigmaX; + result.m_SigmaY = tmpSigmaY; + result.m_SigmaZ = tmpSigmaZ; + + result.m_Altitude = tmpAltitude; + + // read ROI parameters, fill result structure + NodeList rois; + domRoot->GetChildren("roi", rois); + MITK_TEST_CONDITION_REQUIRED( rois.size() == 1, "One ROI defined" ) + itk::DOMNode* roi = rois[0]; + + result.m_RoiMaximumX = GetIntegerAttribute(roi, "maximumX"); + result.m_RoiMinimumX = GetIntegerAttribute(roi, "minimumX"); + result.m_RoiMaximumY = GetIntegerAttribute(roi, "maximumY"); + result.m_RoiMinimumY = GetIntegerAttribute(roi, "minimumY"); + result.m_RoiMaximumZ = GetIntegerAttribute(roi, "maximumZ"); + result.m_RoiMinimumZ = GetIntegerAttribute(roi, "minimumZ"); + + // read statistic parameters, fill result structure + NodeList statistics; + domRoot->GetChildren("statistic", statistics); + MITK_TEST_CONDITION_REQUIRED( statistics.size() == 1, "One statistic defined" ) + itk::DOMNode* statistic = statistics[0]; + + result.m_HotspotMinimum = GetDoubleAttribute(statistic, "minimum"); + result.m_HotspotMaximum = GetDoubleAttribute(statistic, "maximum"); + result.m_HotspotPeak = GetDoubleAttribute(statistic, "peakOptimized"); + + std::vector tmpMinimumIndex(3,0); + + tmpMinimumIndex[0] = GetIntegerAttribute(statistic, "minimumIndexX"); + tmpMinimumIndex[1] = GetIntegerAttribute(statistic, "minimumIndexY"); + tmpMinimumIndex[2] = GetIntegerAttribute(statistic, "minimumIndexZ"); + + result.m_HotspotMinimumIndex = tmpMinimumIndex; + + + std::vector tmpMaximumIndex(3,0); + + tmpMaximumIndex[0] = GetIntegerAttribute(statistic, "maximumIndexX"); + tmpMaximumIndex[1] = GetIntegerAttribute(statistic, "maximumIndexY"); + tmpMaximumIndex[2] = GetIntegerAttribute(statistic, "maximumIndexZ"); + + result.m_HotspotMaximumIndex = tmpMaximumIndex; + + std::vector tmpPeakIndex(3,0); + + tmpPeakIndex[0] = GetIntegerAttribute(statistic, "peakIndexX"); + tmpPeakIndex[1] = GetIntegerAttribute(statistic, "peakIndexY"); + tmpPeakIndex[2] = GetIntegerAttribute(statistic, "peakIndexZ"); + + result.m_HotspotPeakIndex = tmpPeakIndex; + + return result; + } + catch (std::exception& e) + { + MITK_TEST_CONDITION_REQUIRED(false, "Reading test parameters from XML file. Error message: " << e.what()); + } + + if (false /* and all parameters nicely found */) + { + return result; + } + else + { + throw std::invalid_argument("Test called with invalid parameters.."); // TODO provide details if possible + } + } + + /** + \brief Generate an image that contains a couple of 3D gaussian distributions. + + Uses the given parameters to produce a test image using class TODO... bla + */ + + static mitk::Image::Pointer BuildTestImage(const Parameters& testParameters) + { + // evaluate parameters, create corresponding image + mitk::Image::Pointer result; + + typedef double PixelType; + const unsigned int Dimension = 3; + typedef itk::Image ImageType; + ImageType::Pointer image = ImageType::New(); + typedef itk::MultiGaussianImageSource< ImageType > MultiGaussianImageSource; + MultiGaussianImageSource::Pointer gaussianGenerator = MultiGaussianImageSource::New(); + ImageType::SizeValueType size[3]; + size[0] = testParameters.m_ImageColumns; + size[1] = testParameters.m_ImageRows; + size[2] = testParameters.m_ImageSlices; + + itk::MultiGaussianImageSource::VectorType centerXVec, centerYVec, centerZVec, sigmaXVec, sigmaYVec, sigmaZVec, altitudeVec; + + for(int i = 0; i < testParameters.m_NumberOfGaussian; ++i) + { + centerXVec.push_back(testParameters.m_CenterX[i]); + centerYVec.push_back(testParameters.m_CenterY[i]); + centerZVec.push_back(testParameters.m_CenterZ[i]); + + sigmaXVec.push_back(testParameters.m_SigmaX[i]); + sigmaYVec.push_back(testParameters.m_SigmaY[i]); + sigmaZVec.push_back(testParameters.m_SigmaZ[i]); + + altitudeVec.push_back(testParameters.m_Altitude[i]); + } + + ImageType::SpacingType spacing; + + for(int i = 0; i < Dimension; ++i) + spacing[i] = testParameters.m_Spacing[i]; + + gaussianGenerator->SetSize( size ); + gaussianGenerator->SetSpacing( spacing ); + gaussianGenerator->SetRadiusStepNumber(5); + gaussianGenerator->SetRadius(pow(itk::Math::one_over_pi * 0.75 , 1.0 / 3.0) * 10); + gaussianGenerator->SetNumberOfGausssians(testParameters.m_NumberOfGaussian); + + gaussianGenerator->AddGaussian(centerXVec, centerYVec, centerZVec, + sigmaXVec, sigmaYVec, sigmaZVec, altitudeVec); + + gaussianGenerator->Update(); + + image = gaussianGenerator->GetOutput(); + + mitk::CastToMitkImage(image, result); + + return result; + } + + /** + \brief Calculates hotspot statistics for given test image and ROI parameters. + + Uses ImageStatisticsCalculator to find a hotspot in a defined ROI within the given image. + */ + static mitk::ImageStatisticsCalculator::Statistics CalculateStatistics(mitk::Image* image, const Parameters& testParameters) + { + mitk::ImageStatisticsCalculator::Statistics result; + const unsigned int Dimension = 3; + typedef itk::Image MaskImageType; + MaskImageType::Pointer mask = MaskImageType::New(); + + MaskImageType::SizeType size; + MaskImageType::SpacingType spacing; + MaskImageType::IndexType start; + + mitk::ImageStatisticsCalculator::Pointer statisticsCalculator = mitk::ImageStatisticsCalculator::New(); + + for(int i = 0; i < Dimension; ++i) + { + start[i] = 0.00; + spacing[i] = testParameters.m_Spacing[i]; + } + + + + size[0] = testParameters.m_ImageColumns; + size[1] = testParameters.m_ImageRows; + size[2] = testParameters.m_ImageSlices; + + MaskImageType::RegionType region; + region.SetIndex(start); + region.SetSize(size); + + mask->SetSpacing(spacing); + mask->SetRegions(region); + mask->Allocate(); + + for(int x = testParameters.m_RoiMinimumX; x < testParameters.m_RoiMaximumX; ++x) + { + for(int y = testParameters.m_RoiMinimumY; y < testParameters.m_RoiMaximumY; ++y) + { + for(int z = testParameters.m_RoiMinimumZ; z < testParameters.m_RoiMaximumZ; ++z) + { + MaskImageType::IndexType pixelIndex; + pixelIndex[0] = x; + pixelIndex[1] = y; + pixelIndex[2] = z; + + mask->SetPixel(pixelIndex, 1.00); + } + } + } + + mitk::Image::Pointer mitkMaskImage; + mitk::CastToMitkImage(mask, mitkMaskImage); + + statisticsCalculator->SetImage(image); + statisticsCalculator->SetImageMask(mitkMaskImage); + statisticsCalculator->SetMaskingModeToImage(); + statisticsCalculator->ComputeStatistics(); + result = statisticsCalculator->GetStatistics(); + + // create calculator object + // fill parameters (mask, planar figure, etc.) + // execute calculation + // retrieve result and return from function + // handle errors w/o crash! + + return result; + } + + /** + \brief Compares calculated against actual statistics values. + + Checks validness of all statistics aspects. Lets test fail if any aspect is not sufficiently equal. + */ + static void ValidateStatistics(const mitk::ImageStatisticsCalculator::Statistics& statistics, const Parameters& testParameters) + { + // check all expected test result against actual results + + double actualPeakValue = testParameters.m_HotspotPeak; + double expectedPeakValue = statistics.HotspotPeak; + + double actualMaxValue = testParameters.m_HotspotMaximum; + double expectedMaxValue = statistics.HotspotMax; + + double actualMinValue = testParameters.m_HotspotMinimum; + double expectedMinValue = statistics.HotspotMin; + + + //Peak Index + std::vector actualPeakIndex = testParameters.m_HotspotPeakIndex; + vnl_vector expectedVnlPeakIndex; + expectedVnlPeakIndex = statistics.HotspotPeakIndex; + + vnl_vector actualVnlPeakIndex; + actualVnlPeakIndex.set_size(3); + + for(int i = 0; i < actualVnlPeakIndex.size(); ++i) + actualVnlPeakIndex[i] = actualPeakIndex[i]; + + // MaxIndex + std::vector actualMaxIndex = testParameters.m_HotspotMaximumIndex; + vnl_vector expectedVnlMaxIndex; + expectedVnlMaxIndex = statistics.HotspotMaxIndex; + + vnl_vector actualVnlMaxIndex; + actualVnlMaxIndex.set_size(3); + + for(int i = 0; i < actualVnlMaxIndex.size(); ++i) + actualVnlMaxIndex[i] = actualMaxIndex[i]; + + //MinIndex + std::vector actualMinIndex = testParameters.m_HotspotMinimumIndex; + vnl_vector expectedVnlMinIndex; + expectedVnlMinIndex = statistics.HotspotMinIndex; + + vnl_vector actualVnlMinIndex; + actualVnlMinIndex.set_size(3); + + for(int i = 0; i < actualVnlMinIndex.size(); ++i) + actualVnlMinIndex[i] = actualMinIndex[i]; + + double eps = 0.001; + + // float comparisons, allow tiny differences + MITK_TEST_CONDITION( ::fabs(actualPeakValue - expectedPeakValue) < eps, "Actual hotspotPeak value " << actualPeakValue << " (expected " << expectedPeakValue << ")" ); + MITK_TEST_CONDITION( ::fabs(actualMaxValue - expectedMaxValue) < eps, "Actual hotspotMax value " << actualMaxValue << " (expected " << expectedMaxValue << ")" ); + MITK_TEST_CONDITION( ::fabs(actualMinValue - expectedMinValue) < eps, "Actual hotspotMin value " << actualMinValue << " (expected " << expectedMinValue << ")" ); + + MITK_TEST_CONDITION( expectedVnlPeakIndex == actualVnlPeakIndex, "Actual hotspotPeakIndex " << actualVnlPeakIndex << " (expected " << expectedVnlPeakIndex << ")" ); + MITK_TEST_CONDITION( expectedVnlMaxIndex == actualVnlMaxIndex, "Actual hotspotMaxIndex " << actualVnlMaxIndex << " (expected " << expectedVnlMaxIndex << ")" ); + MITK_TEST_CONDITION( expectedVnlMinIndex == actualVnlMinIndex, "Actual hotspotMinIndex " << actualVnlMinIndex << " (expected " << expectedVnlMinIndex << ")" ); + } +}; + + +#include + +/** + \brief Verifies that TODO hotspot statistics part of ImageStatisticsCalculator. + + bla... +*/ +int mitkImageStatisticsHotspotTest(int argc, char* argv[]) +{ + MITK_TEST_BEGIN("mitkImageStatisticsHotspotTest") + try { + // parse commandline parameters (see CMakeLists.txt) + mitkImageStatisticsHotspotTestClass::Parameters parameters = mitkImageStatisticsHotspotTestClass::ParseParameters(argc,argv); + + // build a test image as described in parameters + mitk::Image::Pointer image = mitkImageStatisticsHotspotTestClass::BuildTestImage(parameters); + MITK_TEST_CONDITION_REQUIRED( image.IsNotNull(), "Generate test image" ); + + itk::TimeProbe clock; + clock.Start(); + + // calculate statistics for this image (potentially use parameters for statistics ROI) + mitk::ImageStatisticsCalculator::Statistics statistics = mitkImageStatisticsHotspotTestClass::CalculateStatistics(image, parameters); + + clock.Stop(); + std::cout << "Statistics time consumed: " << clock.GetTotal() << std::endl; + // compare statistics against stored expected values + mitkImageStatisticsHotspotTestClass::ValidateStatistics(statistics, parameters); + + } + catch (std::exception& e) + { + std::cout << "Error: " << e.what() << std::endl; + } + + + MITK_TEST_END() +} diff --git a/Modules/ImageStatistics/Testing/mitkMultiGaussianTest.cpp b/Modules/ImageStatistics/Testing/mitkMultiGaussianTest.cpp new file mode 100644 index 0000000000..dd0787010e --- /dev/null +++ b/Modules/ImageStatistics/Testing/mitkMultiGaussianTest.cpp @@ -0,0 +1,93 @@ +/*=================================================================== + +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 "itkMultiGaussianImageSource.h" + +#include +#include + +int mitkMultiGaussianTest(int, char* []) +{ + // always start with this! + MITK_TEST_BEGIN("mitkMultiGaussianTest") + + typedef double PixelType; + const unsigned int Dimension = 3; + typedef itk::Image ImageType; + typedef itk::MultiGaussianImageSource< ImageType > MultiGaussianImageSource; + itk::MultiGaussianImageSource< ImageType >::VectorType centerXVec, centerYVec, centerZVec, sigmaXVec, sigmaYVec, sigmaZVec, altitudeVec; + MultiGaussianImageSource::Pointer gaussianGenerator = MultiGaussianImageSource::New(); + ImageType::SizeValueType size[3]; + size[0] = 50; + size[1] = 50; + size[2] = 50; + srand (time(NULL)); + unsigned int numberOfGaussians = 7; + unsigned int minWidthOfGaussian = (size[0] + size[1] + size[2]) / 27; // A ninth of the mean image size + unsigned int maxWidthOfGaussian = (size[0] + size[1] + size[2]) / 9; // One-third of the mean image size + unsigned int minAltitudeOfGaussian = 5; + unsigned int maxAltitudeOfGaussian = 200; + double centerX, centerY, centerZ, sigmaX, sigmaY, sigmaZ, altitude; + + gaussianGenerator->SetSize( size ); + gaussianGenerator->SetSpacing( 1 ); + gaussianGenerator->SetRadiusStepNumber(5); + gaussianGenerator->SetRadius(pow(itk::Math::one_over_pi * 0.75 , 1.0 / 3.0) * 10); + gaussianGenerator->SetNumberOfGausssians(numberOfGaussians); + // std::ofstream myfile; + // myfile.open ("C:/temp/tempParameter3.txt"); + // myfile << " CentX \t" << "Y \t" << "Z \t" << "SigX \t" << "Y \t" << "Z \t" << "Altit\n"; + + int numberAddGaussian = numberOfGaussians; + for( unsigned int i = 0; i < numberAddGaussian; ++i) + { + centerX = rand() % size[0]; + centerY = rand() % size[1]; + centerZ = rand() % size[2]; + sigmaX = minWidthOfGaussian + rand() % (maxWidthOfGaussian - minWidthOfGaussian); + sigmaY = minWidthOfGaussian + rand() % (maxWidthOfGaussian - minWidthOfGaussian); + sigmaZ = minWidthOfGaussian + rand() % (maxWidthOfGaussian - minWidthOfGaussian); + altitude = minAltitudeOfGaussian + rand() % (maxAltitudeOfGaussian - minAltitudeOfGaussian); + //gaussianGenerator->AddGaussian(centerX, centerY, centerZ, sigmaX, sigmaY, sigmaZ, altitude); + centerXVec.push_back(centerX); + centerYVec.push_back(centerY); + centerZVec.push_back(centerZ); + sigmaXVec.push_back(sigmaX); + sigmaYVec.push_back(sigmaY); + sigmaZVec.push_back(sigmaZ); + altitudeVec.push_back(altitude); + // myfile <AddGaussian(centerXVec, centerYVec, centerZVec, sigmaXVec, sigmaYVec, sigmaZVec, altitudeVec); + gaussianGenerator->Update(); + gaussianGenerator->CalculateMidpointAndMeanValue(); + std::cout << "Sphere radius is: " << gaussianGenerator->GetRadius() << std::endl; + std::cout << "Sphere midpoint is: " << gaussianGenerator->GetSphereMidpoint() << std::endl; + std::cout << "Mean value is: " << gaussianGenerator->GetMaxMeanValue() << std::endl; + ImageType::Pointer gaussianImage = gaussianGenerator->GetOutput(); + + //File writer + typedef itk::ImageFileWriter< ImageType > WriterType; + WriterType::Pointer writer = WriterType::New(); + writer->SetFileName( "C:/temp/tempImage33.nrrd" ); + writer->SetInput( gaussianImage ); + writer->Update(); + + MITK_TEST_END() +} \ No newline at end of file diff --git a/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp b/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp index 31dabfb6ef..a1526b6c26 100644 --- a/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp +++ b/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp @@ -1,1217 +1,1755 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkImageStatisticsCalculator.h" #include "mitkImageAccessByItk.h" #include "mitkImageCast.h" #include "mitkExtractImageFilter.h" +#include "mitkImageTimeSelector.h" #include #include #include #include #include #include #include +#include +#include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include - +#include +#include +#include +#include +#include +#include + +#include +#include #include +//#define DEBUG_HOTSPOTSEARCH + +#define _USE_MATH_DEFINES +#include + #if ( ( VTK_MAJOR_VERSION <= 5 ) && ( VTK_MINOR_VERSION<=8) ) #include "mitkvtkLassoStencilSource.h" #else #include "vtkLassoStencilSource.h" #endif #include #include +// TODO DM: sort includes, check if they are really needed + namespace mitk { ImageStatisticsCalculator::ImageStatisticsCalculator() : m_MaskingMode( MASKING_MODE_NONE ), m_MaskingModeChanged( false ), m_IgnorePixelValue(0.0), m_DoIgnorePixelValue(false), m_IgnorePixelValueChanged(false), m_PlanarFigureAxis (0), m_PlanarFigureSlice (0), m_PlanarFigureCoordinate0 (0), - m_PlanarFigureCoordinate1 (0) + m_PlanarFigureCoordinate1 (0), // TODO DM: check order of variable initialization + m_HotspotRadiusInMM(6.2035049089940), // radius of a 1cm3 sphere in mm + m_CalculateHotspot(false) { m_EmptyHistogram = HistogramType::New(); m_EmptyHistogram->SetMeasurementVectorSize(1); HistogramType::SizeType histogramSize(1); histogramSize.Fill( 256 ); m_EmptyHistogram->Initialize( histogramSize ); m_EmptyStatistics.Reset(); } ImageStatisticsCalculator::~ImageStatisticsCalculator() { } void ImageStatisticsCalculator::SetImage( const mitk::Image *image ) { if ( m_Image != image ) { m_Image = image; this->Modified(); unsigned int numberOfTimeSteps = image->GetTimeSteps(); // Initialize vectors to time-size of this image m_ImageHistogramVector.resize( numberOfTimeSteps ); m_MaskedImageHistogramVector.resize( numberOfTimeSteps ); m_PlanarFigureHistogramVector.resize( numberOfTimeSteps ); m_ImageStatisticsVector.resize( numberOfTimeSteps ); m_MaskedImageStatisticsVector.resize( numberOfTimeSteps ); m_PlanarFigureStatisticsVector.resize( numberOfTimeSteps ); m_ImageStatisticsTimeStampVector.resize( numberOfTimeSteps ); m_MaskedImageStatisticsTimeStampVector.resize( numberOfTimeSteps ); m_PlanarFigureStatisticsTimeStampVector.resize( numberOfTimeSteps ); m_ImageStatisticsCalculationTriggerVector.resize( numberOfTimeSteps ); m_MaskedImageStatisticsCalculationTriggerVector.resize( numberOfTimeSteps ); m_PlanarFigureStatisticsCalculationTriggerVector.resize( numberOfTimeSteps ); for ( unsigned int t = 0; t < image->GetTimeSteps(); ++t ) { m_ImageStatisticsTimeStampVector[t].Modified(); m_ImageStatisticsCalculationTriggerVector[t] = true; } } } void ImageStatisticsCalculator::SetImageMask( const mitk::Image *imageMask ) { if ( m_Image.IsNull() ) { itkExceptionMacro( << "Image needs to be set first!" ); } if ( m_Image->GetTimeSteps() != imageMask->GetTimeSteps() ) { itkExceptionMacro( << "Image and image mask need to have equal number of time steps!" ); } if ( m_ImageMask != imageMask ) { m_ImageMask = imageMask; this->Modified(); for ( unsigned int t = 0; t < m_Image->GetTimeSteps(); ++t ) { m_MaskedImageStatisticsTimeStampVector[t].Modified(); m_MaskedImageStatisticsCalculationTriggerVector[t] = true; } } } void ImageStatisticsCalculator::SetPlanarFigure( mitk::PlanarFigure *planarFigure ) { if ( m_Image.IsNull() ) { itkExceptionMacro( << "Image needs to be set first!" ); } if ( m_PlanarFigure != planarFigure ) { m_PlanarFigure = planarFigure; this->Modified(); for ( unsigned int t = 0; t < m_Image->GetTimeSteps(); ++t ) { m_PlanarFigureStatisticsTimeStampVector[t].Modified(); m_PlanarFigureStatisticsCalculationTriggerVector[t] = true; } } } void ImageStatisticsCalculator::SetMaskingMode( unsigned int mode ) { if ( m_MaskingMode != mode ) { m_MaskingMode = mode; m_MaskingModeChanged = true; this->Modified(); } } void ImageStatisticsCalculator::SetMaskingModeToNone() { if ( m_MaskingMode != MASKING_MODE_NONE ) { m_MaskingMode = MASKING_MODE_NONE; m_MaskingModeChanged = true; this->Modified(); } } void ImageStatisticsCalculator::SetMaskingModeToImage() { if ( m_MaskingMode != MASKING_MODE_IMAGE ) { m_MaskingMode = MASKING_MODE_IMAGE; m_MaskingModeChanged = true; this->Modified(); } } void ImageStatisticsCalculator::SetMaskingModeToPlanarFigure() { if ( m_MaskingMode != MASKING_MODE_PLANARFIGURE ) { m_MaskingMode = MASKING_MODE_PLANARFIGURE; m_MaskingModeChanged = true; this->Modified(); } } void ImageStatisticsCalculator::SetIgnorePixelValue(double value) { if ( m_IgnorePixelValue != value ) { m_IgnorePixelValue = value; if(m_DoIgnorePixelValue) { m_IgnorePixelValueChanged = true; } this->Modified(); } } double ImageStatisticsCalculator::GetIgnorePixelValue() { return m_IgnorePixelValue; } void ImageStatisticsCalculator::SetDoIgnorePixelValue(bool value) { if ( m_DoIgnorePixelValue != value ) { m_DoIgnorePixelValue = value; m_IgnorePixelValueChanged = true; this->Modified(); } } bool ImageStatisticsCalculator::GetDoIgnorePixelValue() { return m_DoIgnorePixelValue; } +void ImageStatisticsCalculator::SetHotspotRadius(double value) +{ + m_HotspotRadiusInMM = value; +} + +double ImageStatisticsCalculator::GetHotspotRadius() +{ + return m_HotspotRadiusInMM; +} + +void ImageStatisticsCalculator::SetCalculateHotspot(bool value) +{ + m_CalculateHotspot = value; +} + +bool ImageStatisticsCalculator::IsHotspotCalculated() +{ + return m_CalculateHotspot; +} + bool ImageStatisticsCalculator::ComputeStatistics( unsigned int timeStep ) { if (m_Image.IsNull() ) { mitkThrow() << "Image not set!"; } if (!m_Image->IsInitialized()) { mitkThrow() << "Image not initialized!"; } if ( m_Image->GetReferenceCount() == 1 ) { // Image no longer valid; we are the only ones to still hold a reference on it return false; } if ( timeStep >= m_Image->GetTimeSteps() ) { throw std::runtime_error( "Error: invalid time step!" ); } // If a mask was set but we are the only ones to still hold a reference on // it, delete it. if ( m_ImageMask.IsNotNull() && (m_ImageMask->GetReferenceCount() == 1) ) { m_ImageMask = NULL; } - // Check if statistics is already up-to-date unsigned long imageMTime = m_ImageStatisticsTimeStampVector[timeStep].GetMTime(); unsigned long maskedImageMTime = m_MaskedImageStatisticsTimeStampVector[timeStep].GetMTime(); unsigned long planarFigureMTime = m_PlanarFigureStatisticsTimeStampVector[timeStep].GetMTime(); bool imageStatisticsCalculationTrigger = m_ImageStatisticsCalculationTriggerVector[timeStep]; bool maskedImageStatisticsCalculationTrigger = m_MaskedImageStatisticsCalculationTriggerVector[timeStep]; bool planarFigureStatisticsCalculationTrigger = m_PlanarFigureStatisticsCalculationTriggerVector[timeStep]; if ( !m_IgnorePixelValueChanged && ((m_MaskingMode != MASKING_MODE_NONE) || (imageMTime > m_Image->GetMTime() && !imageStatisticsCalculationTrigger)) && ((m_MaskingMode != MASKING_MODE_IMAGE) || (maskedImageMTime > m_ImageMask->GetMTime() && !maskedImageStatisticsCalculationTrigger)) && ((m_MaskingMode != MASKING_MODE_PLANARFIGURE) || (planarFigureMTime > m_PlanarFigure->GetMTime() && !planarFigureStatisticsCalculationTrigger)) ) { // Statistics is up to date! if ( m_MaskingModeChanged ) { m_MaskingModeChanged = false; return true; } else { return false; } } // Reset state changed flag m_MaskingModeChanged = false; m_IgnorePixelValueChanged = false; // Depending on masking mode, extract and/or generate the required image // and mask data from the user input this->ExtractImageAndMask( timeStep ); StatisticsContainer *statisticsContainer; HistogramContainer *histogramContainer; switch ( m_MaskingMode ) { case MASKING_MODE_NONE: default: if(!m_DoIgnorePixelValue) { statisticsContainer = &m_ImageStatisticsVector[timeStep]; histogramContainer = &m_ImageHistogramVector[timeStep]; m_ImageStatisticsTimeStampVector[timeStep].Modified(); m_ImageStatisticsCalculationTriggerVector[timeStep] = false; } else { statisticsContainer = &m_MaskedImageStatisticsVector[timeStep]; histogramContainer = &m_MaskedImageHistogramVector[timeStep]; m_MaskedImageStatisticsTimeStampVector[timeStep].Modified(); m_MaskedImageStatisticsCalculationTriggerVector[timeStep] = false; } break; case MASKING_MODE_IMAGE: statisticsContainer = &m_MaskedImageStatisticsVector[timeStep]; histogramContainer = &m_MaskedImageHistogramVector[timeStep]; m_MaskedImageStatisticsTimeStampVector[timeStep].Modified(); m_MaskedImageStatisticsCalculationTriggerVector[timeStep] = false; break; case MASKING_MODE_PLANARFIGURE: statisticsContainer = &m_PlanarFigureStatisticsVector[timeStep]; histogramContainer = &m_PlanarFigureHistogramVector[timeStep]; m_PlanarFigureStatisticsTimeStampVector[timeStep].Modified(); m_PlanarFigureStatisticsCalculationTriggerVector[timeStep] = false; break; } // Calculate statistics and histogram(s) if ( m_InternalImage->GetDimension() == 3 ) { if ( m_MaskingMode == MASKING_MODE_NONE && !m_DoIgnorePixelValue ) { AccessFixedDimensionByItk_2( m_InternalImage, InternalCalculateStatisticsUnmasked, 3, statisticsContainer, histogramContainer ); } else { AccessFixedDimensionByItk_3( m_InternalImage, InternalCalculateStatisticsMasked, 3, m_InternalImageMask3D.GetPointer(), statisticsContainer, histogramContainer ); } } else if ( m_InternalImage->GetDimension() == 2 ) { if ( m_MaskingMode == MASKING_MODE_NONE && !m_DoIgnorePixelValue ) { AccessFixedDimensionByItk_2( m_InternalImage, InternalCalculateStatisticsUnmasked, 2, statisticsContainer, histogramContainer ); } else { AccessFixedDimensionByItk_3( m_InternalImage, InternalCalculateStatisticsMasked, 2, m_InternalImageMask2D.GetPointer(), statisticsContainer, histogramContainer ); } } else { MITK_ERROR << "ImageStatistics: Image dimension not supported!"; } // Release unused image smart pointers to free memory m_InternalImage = mitk::Image::ConstPointer(); m_InternalImageMask3D = MaskImage3DType::Pointer(); m_InternalImageMask2D = MaskImage2DType::Pointer(); return true; } const ImageStatisticsCalculator::HistogramType * ImageStatisticsCalculator::GetHistogram( unsigned int timeStep, unsigned int label ) const { if ( m_Image.IsNull() || (timeStep >= m_Image->GetTimeSteps()) ) { return NULL; } switch ( m_MaskingMode ) { case MASKING_MODE_NONE: default: { if(m_DoIgnorePixelValue) return m_MaskedImageHistogramVector[timeStep][label]; return m_ImageHistogramVector[timeStep][label]; } case MASKING_MODE_IMAGE: return m_MaskedImageHistogramVector[timeStep][label]; case MASKING_MODE_PLANARFIGURE: return m_PlanarFigureHistogramVector[timeStep][label]; } } const ImageStatisticsCalculator::HistogramContainer & ImageStatisticsCalculator::GetHistogramVector( unsigned int timeStep ) const { if ( m_Image.IsNull() || (timeStep >= m_Image->GetTimeSteps()) ) { return m_EmptyHistogramContainer; } switch ( m_MaskingMode ) { case MASKING_MODE_NONE: default: { if(m_DoIgnorePixelValue) return m_MaskedImageHistogramVector[timeStep]; return m_ImageHistogramVector[timeStep]; } case MASKING_MODE_IMAGE: return m_MaskedImageHistogramVector[timeStep]; case MASKING_MODE_PLANARFIGURE: return m_PlanarFigureHistogramVector[timeStep]; } } const ImageStatisticsCalculator::Statistics & ImageStatisticsCalculator::GetStatistics( unsigned int timeStep, unsigned int label ) const { if ( m_Image.IsNull() || (timeStep >= m_Image->GetTimeSteps()) ) { return m_EmptyStatistics; } switch ( m_MaskingMode ) { case MASKING_MODE_NONE: default: { if(m_DoIgnorePixelValue) return m_MaskedImageStatisticsVector[timeStep][label]; return m_ImageStatisticsVector[timeStep][label]; } case MASKING_MODE_IMAGE: return m_MaskedImageStatisticsVector[timeStep][label]; case MASKING_MODE_PLANARFIGURE: return m_PlanarFigureStatisticsVector[timeStep][label]; } } const ImageStatisticsCalculator::StatisticsContainer & ImageStatisticsCalculator::GetStatisticsVector( unsigned int timeStep ) const { if ( m_Image.IsNull() || (timeStep >= m_Image->GetTimeSteps()) ) { return m_EmptyStatisticsContainer; } switch ( m_MaskingMode ) { case MASKING_MODE_NONE: default: { if(m_DoIgnorePixelValue) return m_MaskedImageStatisticsVector[timeStep]; return m_ImageStatisticsVector[timeStep]; } case MASKING_MODE_IMAGE: return m_MaskedImageStatisticsVector[timeStep]; case MASKING_MODE_PLANARFIGURE: return m_PlanarFigureStatisticsVector[timeStep]; } } void ImageStatisticsCalculator::ExtractImageAndMask( unsigned int timeStep ) { if ( m_Image.IsNull() ) { throw std::runtime_error( "Error: image empty!" ); } if ( timeStep >= m_Image->GetTimeSteps() ) { throw std::runtime_error( "Error: invalid time step!" ); } ImageTimeSelector::Pointer imageTimeSelector = ImageTimeSelector::New(); imageTimeSelector->SetInput( m_Image ); imageTimeSelector->SetTimeNr( timeStep ); imageTimeSelector->UpdateLargestPossibleRegion(); mitk::Image *timeSliceImage = imageTimeSelector->GetOutput(); switch ( m_MaskingMode ) { case MASKING_MODE_NONE: { m_InternalImage = timeSliceImage; m_InternalImageMask2D = NULL; m_InternalImageMask3D = NULL; if(m_DoIgnorePixelValue) { if( m_InternalImage->GetDimension() == 3 ) { CastToItkImage( timeSliceImage, m_InternalImageMask3D ); m_InternalImageMask3D->FillBuffer(1); } if( m_InternalImage->GetDimension() == 2 ) { CastToItkImage( timeSliceImage, m_InternalImageMask2D ); m_InternalImageMask2D->FillBuffer(1); } } break; } case MASKING_MODE_IMAGE: { if ( m_ImageMask.IsNotNull() && (m_ImageMask->GetReferenceCount() > 1) ) { if ( timeStep < m_ImageMask->GetTimeSteps() ) { ImageTimeSelector::Pointer maskedImageTimeSelector = ImageTimeSelector::New(); maskedImageTimeSelector->SetInput( m_ImageMask ); maskedImageTimeSelector->SetTimeNr( timeStep ); maskedImageTimeSelector->UpdateLargestPossibleRegion(); mitk::Image *timeSliceMaskedImage = maskedImageTimeSelector->GetOutput(); m_InternalImage = timeSliceImage; CastToItkImage( timeSliceMaskedImage, m_InternalImageMask3D ); } else { throw std::runtime_error( "Error: image mask has not enough time steps!" ); } } else { throw std::runtime_error( "Error: image mask empty!" ); } break; } case MASKING_MODE_PLANARFIGURE: { m_InternalImageMask2D = NULL; if ( m_PlanarFigure.IsNull() ) { throw std::runtime_error( "Error: planar figure empty!" ); } if ( !m_PlanarFigure->IsClosed() ) { throw std::runtime_error( "Masking not possible for non-closed figures" ); } const Geometry3D *imageGeometry = timeSliceImage->GetGeometry(); if ( imageGeometry == NULL ) { throw std::runtime_error( "Image geometry invalid!" ); } const Geometry2D *planarFigureGeometry2D = m_PlanarFigure->GetGeometry2D(); if ( planarFigureGeometry2D == NULL ) { throw std::runtime_error( "Planar-Figure not yet initialized!" ); } const PlaneGeometry *planarFigureGeometry = dynamic_cast< const PlaneGeometry * >( planarFigureGeometry2D ); if ( planarFigureGeometry == NULL ) { throw std::runtime_error( "Non-planar planar figures not supported!" ); } // Find principal direction of PlanarFigure in input image unsigned int axis; if ( !this->GetPrincipalAxis( imageGeometry, planarFigureGeometry->GetNormal(), axis ) ) { throw std::runtime_error( "Non-aligned planar figures not supported!" ); } m_PlanarFigureAxis = axis; // Find slice number corresponding to PlanarFigure in input image MaskImage3DType::IndexType index; imageGeometry->WorldToIndex( planarFigureGeometry->GetOrigin(), index ); unsigned int slice = index[axis]; m_PlanarFigureSlice = slice; // Extract slice with given position and direction from image - ExtractImageFilter::Pointer imageExtractor = ExtractImageFilter::New(); - imageExtractor->SetInput( timeSliceImage ); - imageExtractor->SetSliceDimension( axis ); - imageExtractor->SetSliceIndex( slice ); - imageExtractor->Update(); - m_InternalImage = imageExtractor->GetOutput(); + unsigned int dimension = timeSliceImage->GetDimension(); + if (dimension != 2) + { + ExtractImageFilter::Pointer imageExtractor = ExtractImageFilter::New(); + imageExtractor->SetInput( timeSliceImage ); + imageExtractor->SetSliceDimension( axis ); + imageExtractor->SetSliceIndex( slice ); + imageExtractor->Update(); + m_InternalImage = imageExtractor->GetOutput(); + } + else + { + m_InternalImage = timeSliceImage; + } // Compute mask from PlanarFigure AccessFixedDimensionByItk_1( m_InternalImage, InternalCalculateMaskFromPlanarFigure, 2, axis ); } } if(m_DoIgnorePixelValue) { if ( m_InternalImage->GetDimension() == 3 ) { AccessFixedDimensionByItk_1( m_InternalImage, InternalMaskIgnoredPixels, 3, m_InternalImageMask3D.GetPointer() ); } else if ( m_InternalImage->GetDimension() == 2 ) { AccessFixedDimensionByItk_1( m_InternalImage, InternalMaskIgnoredPixels, 2, m_InternalImageMask2D.GetPointer() ); } } } bool ImageStatisticsCalculator::GetPrincipalAxis( const Geometry3D *geometry, Vector3D vector, unsigned int &axis ) { vector.Normalize(); for ( unsigned int i = 0; i < 3; ++i ) { Vector3D axisVector = geometry->GetAxisVector( i ); axisVector.Normalize(); if ( fabs( fabs( axisVector * vector ) - 1.0) < mitk::eps ) { axis = i; return true; } } return false; } template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator::InternalCalculateStatisticsUnmasked( const itk::Image< TPixel, VImageDimension > *image, StatisticsContainer *statisticsContainer, HistogramContainer* histogramContainer ) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< unsigned short, VImageDimension > MaskImageType; typedef typename ImageType::IndexType IndexType; typedef itk::Statistics::ScalarImageToHistogramGenerator< ImageType > HistogramGeneratorType; statisticsContainer->clear(); histogramContainer->clear(); // Progress listening... typedef itk::SimpleMemberCommand< ImageStatisticsCalculator > ITKCommandType; ITKCommandType::Pointer progressListener; progressListener = ITKCommandType::New(); progressListener->SetCallbackFunction( this, &ImageStatisticsCalculator::UnmaskedStatisticsProgressUpdate ); // Issue 100 artificial progress events since ScalarIMageToHistogramGenerator // does not (yet?) support progress reporting this->InvokeEvent( itk::StartEvent() ); for ( unsigned int i = 0; i < 100; ++i ) { this->UnmaskedStatisticsProgressUpdate(); } // Calculate statistics (separate filter) typedef itk::StatisticsImageFilter< ImageType > StatisticsFilterType; typename StatisticsFilterType::Pointer statisticsFilter = StatisticsFilterType::New(); statisticsFilter->SetInput( image ); unsigned long observerTag = statisticsFilter->AddObserver( itk::ProgressEvent(), progressListener ); statisticsFilter->Update(); statisticsFilter->RemoveObserver( observerTag ); this->InvokeEvent( itk::EndEvent() ); // Calculate minimum and maximum typedef itk::MinimumMaximumImageCalculator< ImageType > MinMaxFilterType; typename MinMaxFilterType::Pointer minMaxFilter = MinMaxFilterType::New(); minMaxFilter->SetImage( image ); unsigned long observerTag2 = minMaxFilter->AddObserver( itk::ProgressEvent(), progressListener ); minMaxFilter->Compute(); minMaxFilter->RemoveObserver( observerTag2 ); this->InvokeEvent( itk::EndEvent() ); Statistics statistics; statistics.Reset(); statistics.Label = 1; statistics.N = image->GetBufferedRegion().GetNumberOfPixels(); statistics.Min = statisticsFilter->GetMinimum(); statistics.Max = statisticsFilter->GetMaximum(); statistics.Mean = statisticsFilter->GetMean(); statistics.Median = 0.0; statistics.Sigma = statisticsFilter->GetSigma(); statistics.RMS = sqrt( statistics.Mean * statistics.Mean + statistics.Sigma * statistics.Sigma ); statistics.MinIndex.set_size(image->GetImageDimension()); statistics.MaxIndex.set_size(image->GetImageDimension()); - for (int i=0; iGetIndexOfMaximum()[i]; statistics.MinIndex[i] = minMaxFilter->GetIndexOfMinimum()[i]; } statisticsContainer->push_back( statistics ); // Calculate histogram typename HistogramGeneratorType::Pointer histogramGenerator = HistogramGeneratorType::New(); histogramGenerator->SetInput( image ); histogramGenerator->SetMarginalScale( 100 ); histogramGenerator->SetNumberOfBins( 768 ); histogramGenerator->SetHistogramMin( statistics.Min ); histogramGenerator->SetHistogramMax( statistics.Max ); histogramGenerator->Compute(); + // TODO DM: add hotspot search here! + histogramContainer->push_back( histogramGenerator->GetOutput() ); } template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator::InternalMaskIgnoredPixels( const itk::Image< TPixel, VImageDimension > *image, itk::Image< unsigned short, VImageDimension > *maskImage ) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< unsigned short, VImageDimension > MaskImageType; itk::ImageRegionIterator itmask(maskImage, maskImage->GetLargestPossibleRegion()); itk::ImageRegionConstIterator itimage(image, image->GetLargestPossibleRegion()); itmask.GoToBegin(); itimage.GoToBegin(); while( !itmask.IsAtEnd() ) { if(m_IgnorePixelValue == itimage.Get()) { itmask.Set(0); } ++itmask; ++itimage; } } template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator::InternalCalculateStatisticsMasked( const itk::Image< TPixel, VImageDimension > *image, itk::Image< unsigned short, VImageDimension > *maskImage, StatisticsContainer* statisticsContainer, HistogramContainer* histogramContainer ) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< unsigned short, VImageDimension > MaskImageType; typedef typename ImageType::IndexType IndexType; typedef typename ImageType::PointType PointType; typedef typename ImageType::SpacingType SpacingType; typedef itk::LabelStatisticsImageFilter< ImageType, MaskImageType > LabelStatisticsFilterType; typedef itk::ChangeInformationImageFilter< MaskImageType > ChangeInformationFilterType; typedef itk::ExtractImageFilter< ImageType, ImageType > ExtractImageFilterType; statisticsContainer->clear(); histogramContainer->clear(); - // Make sure that mask is set if ( maskImage == NULL ) { itkExceptionMacro( << "Mask image needs to be set!" ); } // Make sure that spacing of mask and image are the same SpacingType imageSpacing = image->GetSpacing(); SpacingType maskSpacing = maskImage->GetSpacing(); PointType zeroPoint; zeroPoint.Fill( 0.0 ); if ( (zeroPoint + imageSpacing).SquaredEuclideanDistanceTo( (zeroPoint + maskSpacing) ) > mitk::eps ) { itkExceptionMacro( << "Mask needs to have same spacing as image! (Image spacing: " << imageSpacing << "; Mask spacing: " << maskSpacing << ")" ); } // Make sure that orientation of mask and image are the same typedef typename ImageType::DirectionType DirectionType; DirectionType imageDirection = image->GetDirection(); DirectionType maskDirection = maskImage->GetDirection(); for( int i = 0; i < imageDirection.ColumnDimensions; ++i ) { for( int j = 0; j < imageDirection.ColumnDimensions; ++j ) { double differenceDirection = imageDirection[i][j] - maskDirection[i][j]; if ( fabs( differenceDirection ) > mitk::eps ) { itkExceptionMacro( << "Mask needs to have same direction as image! (Image direction: " << imageDirection << "; Mask direction: " << maskDirection << ")" ); } } } // Make sure that the voxels of mask and image are correctly "aligned", i.e., voxel boundaries are the same in both images PointType imageOrigin = image->GetOrigin(); PointType maskOrigin = maskImage->GetOrigin(); long offset[ImageType::ImageDimension]; typedef itk::ContinuousIndex ContinousIndexType; ContinousIndexType maskOriginContinousIndex, imageOriginContinousIndex; image->TransformPhysicalPointToContinuousIndex(maskOrigin, maskOriginContinousIndex); image->TransformPhysicalPointToContinuousIndex(imageOrigin, imageOriginContinousIndex); for ( unsigned int i = 0; i < ImageType::ImageDimension; ++i ) { double misalignment = maskOriginContinousIndex[i] - floor( maskOriginContinousIndex[i] + 0.5 ); if ( fabs( misalignment ) > mitk::eps ) { itkExceptionMacro( << "Pixels/voxels of mask and image are not sufficiently aligned! (Misalignment: " << misalignment << ")" ); } double indexCoordDistance = maskOriginContinousIndex[i] - imageOriginContinousIndex[i]; offset[i] = (int) indexCoordDistance + image->GetBufferedRegion().GetIndex()[i]; } // Adapt the origin and region (index/size) of the mask so that the origin of both are the same typename ChangeInformationFilterType::Pointer adaptMaskFilter; adaptMaskFilter = ChangeInformationFilterType::New(); adaptMaskFilter->ChangeOriginOn(); adaptMaskFilter->ChangeRegionOn(); adaptMaskFilter->SetInput( maskImage ); adaptMaskFilter->SetOutputOrigin( image->GetOrigin() ); adaptMaskFilter->SetOutputOffset( offset ); adaptMaskFilter->Update(); typename MaskImageType::Pointer adaptedMaskImage = adaptMaskFilter->GetOutput(); // Make sure that mask region is contained within image region if ( !image->GetLargestPossibleRegion().IsInside( adaptedMaskImage->GetLargestPossibleRegion() ) ) { itkExceptionMacro( << "Mask region needs to be inside of image region! (Image region: " << image->GetLargestPossibleRegion() << "; Mask region: " << adaptedMaskImage->GetLargestPossibleRegion() << ")" ); } // If mask region is smaller than image region, extract the sub-sampled region from the original image typename ImageType::SizeType imageSize = image->GetBufferedRegion().GetSize(); typename ImageType::SizeType maskSize = maskImage->GetBufferedRegion().GetSize(); bool maskSmallerImage = false; for ( unsigned int i = 0; i < ImageType::ImageDimension; ++i ) { if ( maskSize[i] < imageSize[i] ) { maskSmallerImage = true; } } typename ImageType::ConstPointer adaptedImage; if ( maskSmallerImage ) { typename ExtractImageFilterType::Pointer extractImageFilter = ExtractImageFilterType::New(); extractImageFilter->SetInput( image ); extractImageFilter->SetExtractionRegion( adaptedMaskImage->GetBufferedRegion() ); extractImageFilter->Update(); adaptedImage = extractImageFilter->GetOutput(); } else { adaptedImage = image; } // Initialize Filter typedef itk::StatisticsImageFilter< ImageType > StatisticsFilterType; typename StatisticsFilterType::Pointer statisticsFilter = StatisticsFilterType::New(); statisticsFilter->SetInput( adaptedImage ); statisticsFilter->Update(); int numberOfBins = ( m_DoIgnorePixelValue && (m_MaskingMode == MASKING_MODE_NONE) ) ? 768 : 384; typename LabelStatisticsFilterType::Pointer labelStatisticsFilter; labelStatisticsFilter = LabelStatisticsFilterType::New(); labelStatisticsFilter->SetInput( adaptedImage ); labelStatisticsFilter->SetLabelInput( adaptedMaskImage ); labelStatisticsFilter->UseHistogramsOn(); labelStatisticsFilter->SetHistogramParameters( numberOfBins, statisticsFilter->GetMinimum(), statisticsFilter->GetMaximum() ); // Add progress listening typedef itk::SimpleMemberCommand< ImageStatisticsCalculator > ITKCommandType; ITKCommandType::Pointer progressListener; progressListener = ITKCommandType::New(); progressListener->SetCallbackFunction( this, &ImageStatisticsCalculator::MaskedStatisticsProgressUpdate ); unsigned long observerTag = labelStatisticsFilter->AddObserver( itk::ProgressEvent(), progressListener ); // Execute filter this->InvokeEvent( itk::StartEvent() ); // Make sure that only the mask region is considered (otherwise, if the mask region is smaller // than the image region, the Update() would result in an exception). labelStatisticsFilter->GetOutput()->SetRequestedRegion( adaptedMaskImage->GetLargestPossibleRegion() ); // Execute the filter labelStatisticsFilter->Update(); this->InvokeEvent( itk::EndEvent() ); labelStatisticsFilter->RemoveObserver( observerTag ); // Find all relevant labels of mask (other than 0) std::list< int > relevantLabels; bool maskNonEmpty = false; unsigned int i; for ( i = 1; i < 4096; ++i ) { if ( labelStatisticsFilter->HasLabel( i ) ) { relevantLabels.push_back( i ); maskNonEmpty = true; } } + if ( maskNonEmpty ) { + Statistics statistics; std::list< int >::iterator it; for ( it = relevantLabels.begin(), i = 0; it != relevantLabels.end(); ++it, ++i ) { histogramContainer->push_back( HistogramType::ConstPointer( labelStatisticsFilter->GetHistogram( (*it) ) ) ); - Statistics statistics; statistics.Label = (*it); statistics.N = labelStatisticsFilter->GetCount( *it ); statistics.Min = labelStatisticsFilter->GetMinimum( *it ); statistics.Max = labelStatisticsFilter->GetMaximum( *it ); statistics.Mean = labelStatisticsFilter->GetMean( *it ); statistics.Median = labelStatisticsFilter->GetMedian( *it ); statistics.Sigma = labelStatisticsFilter->GetSigma( *it ); statistics.RMS = sqrt( statistics.Mean * statistics.Mean + statistics.Sigma * statistics.Sigma ); // restrict image to mask area for min/max index calculation typedef itk::MaskImageFilter< ImageType, MaskImageType, ImageType > MaskImageFilterType; typename MaskImageFilterType::Pointer masker = MaskImageFilterType::New(); masker->SetOutsideValue( (statistics.Min+statistics.Max)/2 ); masker->SetInput1(adaptedImage); masker->SetInput2(adaptedMaskImage); masker->Update(); // get index of minimum and maximum typedef itk::MinimumMaximumImageCalculator< ImageType > MinMaxFilterType; typename MinMaxFilterType::Pointer minMaxFilter = MinMaxFilterType::New(); minMaxFilter->SetImage( masker->GetOutput() ); unsigned long observerTag2 = minMaxFilter->AddObserver( itk::ProgressEvent(), progressListener ); minMaxFilter->Compute(); minMaxFilter->RemoveObserver( observerTag2 ); this->InvokeEvent( itk::EndEvent() ); statistics.MinIndex.set_size(adaptedImage->GetImageDimension()); statistics.MaxIndex.set_size(adaptedImage->GetImageDimension()); typename MinMaxFilterType::IndexType tempMaxIndex = minMaxFilter->GetIndexOfMaximum(); typename MinMaxFilterType::IndexType tempMinIndex = minMaxFilter->GetIndexOfMinimum(); // FIX BUG 14644 //If a PlanarFigure is used for segmentation the //adaptedImage is a single slice (2D). Adding the // 3. dimension. if (m_MaskingMode == MASKING_MODE_PLANARFIGURE && m_Image->GetDimension()==3) { statistics.MaxIndex.set_size(m_Image->GetDimension()); statistics.MaxIndex[m_PlanarFigureCoordinate0]=tempMaxIndex[0]; statistics.MaxIndex[m_PlanarFigureCoordinate1]=tempMaxIndex[1]; statistics.MaxIndex[m_PlanarFigureAxis]=m_PlanarFigureSlice; statistics.MinIndex.set_size(m_Image->GetDimension()); statistics.MinIndex[m_PlanarFigureCoordinate0]=tempMinIndex[0]; statistics.MinIndex[m_PlanarFigureCoordinate1]=tempMinIndex[1]; statistics.MinIndex[m_PlanarFigureAxis]=m_PlanarFigureSlice; } else { - for (int i = 0; ipush_back( statistics ); + // TODO DM: what about different label values? ImageStatisticsCalculator usually calculates statistics sets for EACH label in the given mask + // TODO DM: it would be more consistent if we calculate hotspot statistics for EACH label, not only for the "unequal 0" label (after all other TODOs) + /*****************************************************Calculate Hotspot Statistics**********************************************/ + + if(IsHotspotCalculated()) + { + // TODO DM: CalculateHotspotStatistics should + // 1. regard mask + // 2. calculate a hotspot (and its statistics) per mask label/value + // 3. use LabelStatisticsImageFilter where possible + Statistics hotspotStatistics = CalculateHotspotStatistics (adaptedImage.GetPointer(), adaptedMaskImage.GetPointer(), GetHotspotRadius()); + statistics.HotspotMax = hotspotStatistics.HotspotMax; + statistics.HotspotMin = hotspotStatistics.HotspotMin; + statistics.HotspotMean = hotspotStatistics.HotspotMean; + statistics.HotspotMaxIndex = hotspotStatistics.HotspotMaxIndex; + statistics.HotspotMinIndex = hotspotStatistics.HotspotMinIndex; + statistics.HotspotIndex = hotspotStatistics.HotspotIndex; + // TODO DM: add other statistics: N, RMS, ... ; clear role of peak/mean } + statisticsContainer->push_back( statistics ); } else { histogramContainer->push_back( HistogramType::ConstPointer( m_EmptyHistogram ) ); - statisticsContainer->push_back( Statistics() );; + statisticsContainer->push_back( Statistics() ); // TODO DM: this is uninitialized! (refactor into real class!) } } +// TODO DM: needs to be modified to calculate a specific or multiple(!) labels +template +ImageStatisticsCalculator::MinMaxIndex ImageStatisticsCalculator::CalculateMinMaxIndex( + const itk::Image *inputImage, + itk::Image *maskImage) +{ + typedef itk::Image< TPixel, VImageDimension > ImageType; + typedef itk::Image< unsigned short, VImageDimension > MaskImageType; + + typedef itk::ImageRegionConstIterator MaskImageIteratorType; + typedef itk::ImageRegionConstIteratorWithIndex InputImageIndexIteratorType; + + MaskImageIteratorType maskIt(maskImage, maskImage->GetLargestPossibleRegion()); // TODO DM: we should use the same regions here + InputImageIndexIteratorType imageIndexIt(inputImage, inputImage->GetLargestPossibleRegion()); + + float maxValue = itk::NumericTraits::min(); // TODO DM: I DID correct this before: use named functions instead of using - + float minValue = itk::NumericTraits::max(); + + typename ImageType::IndexType maxIndex; + typename ImageType::IndexType minIndex; + + for(maskIt.GoToBegin(), imageIndexIt.GoToBegin(); + !maskIt.IsAtEnd() && !imageIndexIt.IsAtEnd(); + ++maskIt, ++imageIndexIt) + { + if(maskIt.Get() > itk::NumericTraits::Zero) // TODO DM: this is where multiple mask values could be used + { + double value = imageIndexIt.Get(); + + //Calculate minimum, maximum and corresponding index-values + if( value > maxValue ) + { + maxIndex = imageIndexIt.GetIndex(); + maxValue = value; + } + + if(value < minValue ) + { + minIndex = imageIndexIt.GetIndex(); + minValue = value; + } + } + } + + MinMaxIndex minMax; + + minMax.MinIndex.set_size(inputImage->GetImageDimension()); + minMax.MaxIndex.set_size(inputImage->GetImageDimension()); + + 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::SmartPointer< itk::Image > +ImageStatisticsCalculator +::GenerateHotspotSearchConvolutionMask(double spacing[VImageDimension], double radiusInMM) +{ + double radiusInMMSquared = radiusInMM * radiusInMM; + typedef itk::Image< float, VImageDimension > MaskImageType; + typename MaskImageType::Pointer convolutionMask = MaskImageType::New(); + + // Calculate size and allocate mask image + typedef typename MaskImageType::IndexType IndexType; + IndexType maskIndex; + maskIndex.Fill(0); + + typedef typename MaskImageType::SizeType SizeType; + SizeType maskSize; + + Point3D convolutionMaskCenter; convolutionMaskCenter.Fill(0.0); + for(unsigned int i = 0; i < VImageDimension; ++i) + { + maskSize[i] = ::ceil( 2.0 * radiusInMM / spacing[i] ); + + // We always need an uneven size to determine a clear center point in the convolution mask // TODO DM: actually this is not true, is it? I don't see a reason + if(maskSize[i] % 2 == 0 ) + { + ++maskSize[i]; + } + + // TODO DM: center is wrong (below is the corrected calculation) + // convolutionMaskCenterCoordinate[i] = (maskSize[i] -1) / 2; + // coordinates are center based: with 1 pixel the center is at 0.0 + // coordinates are center based: with 2 pixels the center is at 0.5 + // coordinates are center based: with 3 pixels the center is at 1.0 + // coordinates are center based: with 4 pixels the center is at 1.5 + // etc. + convolutionMaskCenter[i] = 0.5 * (double)(maskSize[i]-1); + } + + typedef typename MaskImageType::RegionType RegionType; + RegionType maskRegion; + maskRegion.SetSize(maskSize); + maskRegion.SetIndex(maskIndex); + + convolutionMask->SetRegions(maskRegion); + convolutionMask->SetSpacing(spacing); + convolutionMask->Allocate(); + + // Fill mask image values by subsampling the image grid + typedef itk::ImageRegionIteratorWithIndex MaskIteratorType; + MaskIteratorType maskIt(convolutionMask,maskRegion); + + int numberOfSubVoxelsPerDimension = 2; // per dimension! + int numberOfSubVoxels = ::pow( static_cast(numberOfSubVoxelsPerDimension), static_cast(VImageDimension) ); + double subVoxelSize = 1.0 / (double)numberOfSubVoxelsPerDimension; //(double)numberOfSubVoxels; + double valueOfOneSubVoxel = 1.0 / (double)numberOfSubVoxels; + double maskValue = 0.0; + Point3D subVoxelPosition; + double distanceSquared = 0.0; + + typedef itk::ContinuousIndex ContinuousIndexType; + for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt) + { + ContinuousIndexType indexPoint(maskIt.GetIndex()); + Point3D voxelPosition; + for (unsigned int dimension = 0; dimension < VImageDimension; ++dimension) + { + voxelPosition[dimension] = indexPoint[dimension]; + } + + // TODO DM: regard all dimensions, including z! (former code used only x/y) + // TODO DM: generalize: not x, y, z but a for loop over dimension + // TODO DM: this could be done by calling a recursive method, handing over the "remaining number of dimensions to iterate" + + maskValue = 0.0; + Vector3D subVoxelOffset; subVoxelOffset.Fill(0.0); + // iterate sub-voxels by iterating all possible offsets + for (subVoxelOffset[0] = -0.5 + subVoxelSize / 2.0; + subVoxelOffset[0] < +0.5; + subVoxelOffset[0] += subVoxelSize) + { + for (subVoxelOffset[1] = -0.5 + subVoxelSize / 2.0; + subVoxelOffset[1] < +0.5; + subVoxelOffset[1] += subVoxelSize) + { + for (subVoxelOffset[2] = -0.5 + subVoxelSize / 2.0; + subVoxelOffset[2] < +0.5; + subVoxelOffset[2] += subVoxelSize) + { + subVoxelPosition = voxelPosition + subVoxelOffset; // TODO DM: this COULD be integrated into the for-loops if neccessary (add voxelPosition to initializer and end condition) + //if ( subVoxelPosition.EuclideanDistanceTo( convolutionMaskCenter ) < radiusInMM ) // TODO DM: this is too much matrix operations, we calculate ourselves, check if this time is relevant + distanceSquared = (subVoxelPosition[0]-convolutionMaskCenter[0]) / spacing[0] * (subVoxelPosition[0]-convolutionMaskCenter[0]) / spacing[0] + + (subVoxelPosition[1]-convolutionMaskCenter[1]) / spacing[1] * (subVoxelPosition[1]-convolutionMaskCenter[1]) / spacing[1] + + (subVoxelPosition[2]-convolutionMaskCenter[2]) / spacing[2] * (subVoxelPosition[2]-convolutionMaskCenter[2]) / spacing[2]; + + if (distanceSquared <= radiusInMMSquared) + { + maskValue += valueOfOneSubVoxel; + } + } + } + } + maskIt.Set( maskValue ); + } + return convolutionMask; +} + + +// TODO DM: should be refactored into multiple smaller methosd. This one is too large +template < typename TPixel, unsigned int VImageDimension> +ImageStatisticsCalculator::Statistics ImageStatisticsCalculator::CalculateHotspotStatistics( + const itk::Image* inputImage, + itk::Image* maskImage, // TODO DM: this parameter is completely ignored, although the method is currently ONLY called in the masked input case + double radiusInMM) +{ + typedef itk::Image< TPixel, VImageDimension > InputImageType; + typedef itk::Image< float, VImageDimension > MaskImageType; + + double spacing[VImageDimension]; + for (unsigned int dimension = 0; dimension < VImageDimension; ++dimension) + { + spacing[dimension] = inputImage->GetSpacing()[dimension]; + } + + typename MaskImageType::Pointer convolutionMask = this->GenerateHotspotSearchConvolutionMask(spacing, radiusInMM); + + typedef typename InputImageType::IndexType IndexType; + typedef typename InputImageType::SizeType SizeType; + typedef typename MaskImageType::PointType PointType; + + // Convolution of spherical mask and input image + + typedef itk::Image< float, VImageDimension > ConvolutionImageType; + typedef itk::FFTConvolutionImageFilter ConvolutionFilterType; // TODO DM: this line said ConvolutionImageFilter before: why?? + typedef itk::ConstantBoundaryCondition BoundaryConditionType; + BoundaryConditionType boundaryCondition; + boundaryCondition.SetConstant(0.0); + + typename ConvolutionFilterType::Pointer convolutionFilter = ConvolutionFilterType::New(); + convolutionFilter->SetBoundaryCondition(&boundaryCondition); + convolutionFilter->SetInput(inputImage); + convolutionFilter->SetKernelImage(convolutionMask); + convolutionFilter->SetNormalize(true); + convolutionFilter->Update(); + // TODO DM: above Update will calculate the convolution image for ALL of the input image + // in cases where we have a masked image (always in the first application use case) + // this is too much! it would be enough to calculate the minimum and maximum index of the mask (in each dimension), + // in order to define a region for convolutionFilter. This cold save significant time (perhaps enough to fall back to ConvolutionFilterType instead of FFTConvolutionImageFilter) + // TODO: performance analysis after these changes! + + typename ConvolutionImageType::Pointer hotspotImage = convolutionFilter->GetOutput(); + hotspotImage->SetSpacing( inputImage->GetSpacing() ); // TODO: only workaround because convolution filter seems to ignore spacing of input image + + + // TODO DM: why a spatial object? Objective here should be to 1. find position and value of maximum value in convolution image + /*****************************************************Creating Hotspot Sphere**********************************************/ + typedef itk::Image SphereMaskImageType; + typename SphereMaskImageType::Pointer hotspotSphere = SphereMaskImageType::New(); + + typedef itk::EllipseSpatialObject EllipseType; + typedef itk::SpatialObjectToImageFilter SpatialObjectToImageFilter; + + double hotspotMean = itk::NumericTraits::min(); + + typename SphereMaskImageType::Pointer croppedRegionMask = SphereMaskImageType::New(); + + typename SphereMaskImageType::IndexType peakStart; + peakStart.Fill(0); + typename SphereMaskImageType::SizeType sphereMaskSize = hotspotImage->GetLargestPossibleRegion().GetSize(); + + // TODO DM: this creates an image of the input image size! + typename SphereMaskImageType::RegionType peakRegion; + peakRegion.SetIndex(peakStart); + peakRegion.SetSize(hotspotImage->GetLargestPossibleRegion().GetSize()); + + croppedRegionMask->SetRegions(peakRegion); + croppedRegionMask->Allocate(); + + int offsetX = static_cast((radiusInMM / spacing[0]) + 0.99999); + int offsetY = static_cast((radiusInMM / spacing[1]) + 0.99999); + int offsetZ = static_cast((radiusInMM / spacing[2]) + 0.99999); + + typedef itk::ImageRegionIteratorWithIndex CroppedImageIteratorType; + CroppedImageIteratorType sphereMaskIt(croppedRegionMask, peakRegion); + + for(sphereMaskIt.GoToBegin(); !sphereMaskIt.IsAtEnd(); ++sphereMaskIt) + { + IndexType index = sphereMaskIt.GetIndex(); + + if((index[0] >= offsetX && index[0] <= sphereMaskSize[0] - offsetX -1) && + (index[1] >= offsetY && index[1] <= sphereMaskSize[1] - offsetY -1) && + (index[2] >= offsetZ && index[2] <= sphereMaskSize[2] - offsetZ -1)) + sphereMaskIt.Set(1); + else + sphereMaskIt.Set(0); + } + + typedef typename itk::Image InputMaskImageType; + typedef itk::ImageRegionIteratorWithIndex MaskImageIteratorType; + MaskImageIteratorType inputMaskIt(maskImage, maskImage->GetLargestPossibleRegion()); + CroppedImageIteratorType sphereMaskIterator(croppedRegionMask, croppedRegionMask->GetLargestPossibleRegion()); + + for(inputMaskIt.GoToBegin(), sphereMaskIterator.GoToBegin(); + !inputMaskIt.IsAtEnd() &&!sphereMaskIterator.IsAtEnd(); + ++inputMaskIt, ++sphereMaskIterator) + { + unsigned int maskValue = inputMaskIt.Get(); + unsigned int sphereMaskValue = sphereMaskIterator.Get(); + + if(maskValue > 0 && sphereMaskValue > 0) + sphereMaskIterator.Set(1); + else + sphereMaskIterator.Set(0); + } + + // TODO DM: sphereMaskIt seems to define a box region where a sphere could fit inside the input image + // this seems to come from an idea that Hannes mentioned and what I commented on in line 1244 + // CONVOLUTION should be restricted to an area where we can possibly find result values (i.e. regions inside the mask) + // in addition, if we require the sphere to be completely contained inside the input image (talk to Mathias/Danial for definition) + // THEN we should reduce the mask image before working with it (and prior to using it as a bounding region for convolution) + // + // Besides the comment above, a spatial object is not useful here. A simple itk::ImageRegion would be enough! (and it would fit into the iterator initialization) + MinMaxIndex peakInformations = CalculateMinMaxIndex(hotspotImage.GetPointer(), croppedRegionMask.GetPointer()); + + hotspotMean = peakInformations.Max; + typename SphereMaskImageType::IndexType hotspotIndex; + for(int i = 0; i < VImageDimension; ++i) + hotspotIndex[i] = peakInformations.MaxIndex[i]; + + typename SphereMaskImageType::SizeType hotspotSphereSize; + typename SphereMaskImageType::SpacingType hotspotSphereSpacing = inputImage->GetSpacing(); // TODO DM: we don't need a third spacing definition; all our calculations are for one and the same image with just one spacing in variable "spacing" + + // TODO DM: remove this and use previously calculated mask size! This is redundant + for(unsigned int i = 0; i < VImageDimension; ++i) + { + + double countIndex = 2.0 * radiusInMM / hotspotSphereSpacing[i]; + + // Rounding up to the next integer by cast + countIndex += 0.9999999; + int castedIndex = static_cast(countIndex); + + // We always have an uneven number in size to determine a center-point in the convolution mask + if(castedIndex % 2 > 0 ) + { + hotspotSphereSize[i] = castedIndex; + } + else + { + hotspotSphereSize[i] = castedIndex +1; + } + } + + // Initialize SpatialObjectoToImageFilter + typename itk::SpatialObjectToImageFilter::Pointer spatialObjectToImageFilter + = SpatialObjectToImageFilter::New(); + + spatialObjectToImageFilter->SetSize(hotspotSphereSize); + spatialObjectToImageFilter->SetSpacing(hotspotSphereSpacing); + + // Creating spatial sphere object + typename EllipseType::Pointer sphere = EllipseType::New(); + sphere->SetRadius(radiusInMM); + typedef typename EllipseType::TransformType TransformType; + typename TransformType::Pointer transform = TransformType::New(); + + transform->SetIdentity(); + + typename TransformType::OutputVectorType translation; + + // Transform sphere on center-position, set pixelValues inside sphere on 1 and update + for(int i = 0; i < VImageDimension; ++i) + translation[i] = static_cast((hotspotSphereSize[i] -1) * hotspotSphereSpacing[i] / 2); + + transform->Translate(translation, false); + + sphere->SetObjectToParentTransform(transform); + + spatialObjectToImageFilter->SetInput(sphere); + + sphere->SetDefaultInsideValue(1.00); + sphere->SetDefaultOutsideValue(0.00); + + spatialObjectToImageFilter->SetUseObjectValue(true); + spatialObjectToImageFilter->SetOutsideValue(0); + + spatialObjectToImageFilter->Update(); + hotspotSphere = spatialObjectToImageFilter->GetOutput(); + + // Calculate new origin for hotspot sphere + + IndexType offsetInIndex; + + for(int i = 0; i < VImageDimension; ++i) + offsetInIndex[i] = hotspotSphereSize[i] / 2; + + typename ConvolutionImageType::PointType hotspotOrigin; + hotspotImage->TransformIndexToPhysicalPoint(hotspotIndex, hotspotOrigin); + + PointType offsetInPhysicalPoint; + hotspotSphere->TransformIndexToPhysicalPoint(offsetInIndex, offsetInPhysicalPoint); + + for(int i = 0; i < VImageDimension; ++i) + hotspotOrigin[i] -= offsetInPhysicalPoint[i]; + + hotspotSphere->SetOrigin(hotspotOrigin); + hotspotSphere->Allocate(); + + /* TODO DM: you don't need all of the above "spatial object sphere" code. + It should be possible to replace all of the below code with a single call + to your CalculateMinMaxIndex method. + */ + +#ifdef DEBUG_HOTSPOTSEARCH + + std::cout << std::endl << std::endl; + std::cout << "hotspotMask: " << std::endl; + unsigned int lastZ = 1000000000; + unsigned int lastY = 1000000000; + + unsigned int hotspotMaskIndexCounter = 0; + + typedef itk::ImageRegionConstIteratorWithIndex SphereMaskIteratorType; + SphereMaskIteratorType hotspotMaskIt(hotspotSphere, hotspotSphere->GetLargestPossibleRegion() ); + + for(hotspotMaskIt.GoToBegin();!hotspotMaskIt.IsAtEnd();++hotspotMaskIt) + { + + double tmp = hotspotMaskIt.Get(); + if (hotspotMaskIt.GetIndex()[1] != lastY) + { + std::cout << std::endl; + lastY = hotspotMaskIt.GetIndex()[1]; + } + if (hotspotMaskIt.GetIndex()[0] != lastZ) + { + std::cout << tmp << " "; + lastZ = hotspotMaskIt.GetIndex()[0]; + } + + hotspotMaskIndexCounter++; + + if(hotspotMaskIndexCounter > hotspotSphereSize[0] * hotspotSphereSize[1] -1) { + std::cout << std::endl; + hotspotMaskIndexCounter = 0; + } + } + + std::cout << std::endl << std::endl; +#endif + + /*********************************Creating cropped inputImage for calculation of hotspot statistics****************************/ + + typename InputImageType::IndexType croppedStart; + hotspotImage->TransformPhysicalPointToIndex(hotspotOrigin,croppedStart); + + typename InputImageType::RegionType::SizeType croppedSize = hotspotSphere->GetLargestPossibleRegion().GetSize(); + typename InputImageType::RegionType inputRegion; + inputRegion.SetIndex(croppedStart); + inputRegion.SetSize(croppedSize); + + typename InputImageType::IndexType croppedOutputStart; + croppedOutputStart.Fill(0); + + typename InputImageType::RegionType croppedOutputRegion; + croppedOutputRegion.SetIndex(croppedOutputStart); + croppedOutputRegion.SetSize(hotspotSphere->GetLargestPossibleRegion().GetSize()); + + typename InputImageType::Pointer croppedOutputImage = InputImageType::New(); + croppedOutputImage->SetRegions(croppedOutputRegion); + croppedOutputImage->Allocate(); + + typedef itk::ImageRegionConstIterator ImageIteratorType; + ImageIteratorType inputIt(inputImage, inputRegion); + + ImageIteratorType croppedOutputImageIt(croppedOutputImage, croppedOutputRegion); + + for(inputIt.GoToBegin(), croppedOutputImageIt.GoToBegin(); !inputIt.IsAtEnd(); ++inputIt, ++croppedOutputImageIt) + { + croppedOutputImage->SetPixel(croppedOutputImageIt.GetIndex(), inputIt.Get()); + } + + // Calculate statistics in Hotspot + MinMaxIndex hotspotInformations; + Statistics hotspotStatistics; + + hotspotInformations = CalculateMinMaxIndex(croppedOutputImage.GetPointer(), hotspotSphere.GetPointer()); + + // Add offset to cropped indices + for(int i = 0; i < VImageDimension; ++i) + { + hotspotInformations.MaxIndex[i] += croppedStart[i]; + hotspotInformations.MinIndex[i] += croppedStart[i]; + } + + hotspotStatistics.HotspotMin = hotspotInformations.Min; + hotspotStatistics.HotspotMinIndex = hotspotInformations.MinIndex; + hotspotStatistics.HotspotMax = hotspotInformations.Max; + hotspotStatistics.HotspotMaxIndex = hotspotInformations.MaxIndex; + hotspotStatistics.HotspotMean = hotspotMean; + + hotspotStatistics.HotspotIndex.set_size(inputImage->GetImageDimension()); + for (int i = 0; i< hotspotStatistics.HotspotIndex.size(); ++i) + { + hotspotStatistics.HotspotIndex[i] = hotspotIndex[i]; + } + + return hotspotStatistics; +} template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator::InternalCalculateMaskFromPlanarFigure( const itk::Image< TPixel, VImageDimension > *image, unsigned int axis ) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::CastImageFilter< ImageType, MaskImage2DType > CastFilterType; // Generate mask image as new image with same header as input image and // initialize with "1". typename CastFilterType::Pointer castFilter = CastFilterType::New(); castFilter->SetInput( image ); castFilter->Update(); castFilter->GetOutput()->FillBuffer( 1 ); // all PolylinePoints of the PlanarFigure are stored in a vtkPoints object. // These points are used by the vtkLassoStencilSource to create // a vtkImageStencil. const mitk::Geometry2D *planarFigureGeometry2D = m_PlanarFigure->GetGeometry2D(); const typename PlanarFigure::PolyLineType planarFigurePolyline = m_PlanarFigure->GetPolyLine( 0 ); const mitk::Geometry3D *imageGeometry3D = m_Image->GetGeometry( 0 ); // Determine x- and y-dimensions depending on principal axis int i0, i1; switch ( axis ) { case 0: i0 = 1; i1 = 2; break; case 1: i0 = 0; i1 = 2; break; case 2: default: i0 = 0; i1 = 1; break; } m_PlanarFigureCoordinate0= i0; m_PlanarFigureCoordinate1= i1; // store the polyline contour as vtkPoints object bool outOfBounds = false; vtkSmartPointer points = vtkSmartPointer::New(); typename PlanarFigure::PolyLineType::const_iterator it; for ( it = planarFigurePolyline.begin(); it != planarFigurePolyline.end(); ++it ) { Point3D point3D; // Convert 2D point back to the local index coordinates of the selected // image planarFigureGeometry2D->Map( it->Point, point3D ); // Polygons (partially) outside of the image bounds can not be processed // further due to a bug in vtkPolyDataToImageStencil if ( !imageGeometry3D->IsInside( point3D ) ) { outOfBounds = true; } imageGeometry3D->WorldToIndex( point3D, point3D ); points->InsertNextPoint( point3D[i0], point3D[i1], 0 ); } // mark a malformed 2D planar figure ( i.e. area = 0 ) as out of bounds // this can happen when all control points of a rectangle lie on the same line = two of the three extents are zero double bounds[6] = {0, 0, 0, 0, 0, 0}; points->GetBounds( bounds ); bool extent_x = (fabs(bounds[0] - bounds[1])) < mitk::eps; bool extent_y = (fabs(bounds[2] - bounds[3])) < mitk::eps; bool extent_z = (fabs(bounds[4] - bounds[5])) < mitk::eps; // throw an exception if a closed planar figure is deformed, i.e. has only one non-zero extent if ( m_PlanarFigure->IsClosed() && ((extent_x && extent_y) || (extent_x && extent_z) || (extent_y && extent_z))) { mitkThrow() << "Figure has a zero area and cannot be used for masking."; } if ( outOfBounds ) { throw std::runtime_error( "Figure at least partially outside of image bounds!" ); } // create a vtkLassoStencilSource and set the points of the Polygon vtkSmartPointer lassoStencil = vtkSmartPointer::New(); lassoStencil->SetShapeToPolygon(); lassoStencil->SetPoints( points ); // Export from ITK to VTK (to use a VTK filter) typedef itk::VTKImageImport< MaskImage2DType > ImageImportType; typedef itk::VTKImageExport< MaskImage2DType > ImageExportType; typename ImageExportType::Pointer itkExporter = ImageExportType::New(); itkExporter->SetInput( castFilter->GetOutput() ); vtkSmartPointer vtkImporter = vtkSmartPointer::New(); this->ConnectPipelines( itkExporter, vtkImporter ); // Apply the generated image stencil to the input image vtkSmartPointer imageStencilFilter = vtkSmartPointer::New(); imageStencilFilter->SetInputConnection( vtkImporter->GetOutputPort() ); imageStencilFilter->SetStencil( lassoStencil->GetOutput() ); imageStencilFilter->ReverseStencilOff(); imageStencilFilter->SetBackgroundValue( 0 ); imageStencilFilter->Update(); // Export from VTK back to ITK vtkSmartPointer vtkExporter = vtkImageExport::New(); // TODO: this is WRONG, should be vtkSmartPointer::New(), but bug # 14455 vtkExporter->SetInputConnection( imageStencilFilter->GetOutputPort() ); vtkExporter->Update(); typename ImageImportType::Pointer itkImporter = ImageImportType::New(); this->ConnectPipelines( vtkExporter, itkImporter ); itkImporter->Update(); // Store mask m_InternalImageMask2D = itkImporter->GetOutput(); } void ImageStatisticsCalculator::UnmaskedStatisticsProgressUpdate() { // Need to throw away every second progress event to reach a final count of // 100 since two consecutive filters are used in this case static int updateCounter = 0; if ( updateCounter++ % 2 == 0 ) { this->InvokeEvent( itk::ProgressEvent() ); } } void ImageStatisticsCalculator::MaskedStatisticsProgressUpdate() { this->InvokeEvent( itk::ProgressEvent() ); } } diff --git a/Modules/ImageStatistics/mitkImageStatisticsCalculator.h b/Modules/ImageStatistics/mitkImageStatisticsCalculator.h index 34288e8ca7..6b48fd18a7 100644 --- a/Modules/ImageStatistics/mitkImageStatisticsCalculator.h +++ b/Modules/ImageStatistics/mitkImageStatisticsCalculator.h @@ -1,328 +1,397 @@ /*=================================================================== - The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ +#ifndef mitkImageStatisticsCalculator_h +#define mitkImageStatisticsCalculator_h -#ifndef _MITK_IMAGESTATISTICSCALCULATOR_H -#define _MITK_IMAGESTATISTICSCALCULATOR_H - -#include -#include "ImageStatisticsExports.h" -#include -#include +#include "mitkImage.h" +#include "mitkPlanarFigure.h" +// TODO DM: why the ifndef? #ifndef __itkHistogram_h #include #endif +#include -#include "mitkImage.h" -#include "mitkImageTimeSelector.h" -#include "mitkPlanarFigure.h" #include +#include "ImageStatisticsExports.h" + namespace mitk { /** * \brief Class for calculating statistics and histogram for an (optionally * masked) image. * * Images can be masked by either a label image (of the same dimensions as * the original image) or by a closed mitk::PlanarFigure, e.g. a circle or * polygon. When masking with a planar figure, the slice corresponding to the * plane containing the figure is extracted and then clipped with contour * defined by the figure. Planar figures need to be aligned along the main axes * of the image (axial, sagittal, coronal). Planar figures on arbitrary * rotated planes are not supported. * * For each operating mode (no masking, masking by image, masking by planar * figure), the calculated statistics and histogram are cached so that, when * switching back and forth between operation modes without modifying mask or * image, the information doesn't need to be recalculated. * + * The class also has the possibility to calculate minimum, maximum, mean + * and their corresponding indicies in the hottest spot in a given ROI / VOI. + * The size of the hotspot is defined by a sphere with a radius specified by + * the user. This procedure is required for the calculation of SUV-statistics + * in PET-images for example. + * * Note: currently time-resolved and multi-channel pictures are not properly * supported. */ class ImageStatistics_EXPORT ImageStatisticsCalculator : public itk::Object { public: + /** + TODO DM: document + */ enum { MASKING_MODE_NONE = 0, - MASKING_MODE_IMAGE, - MASKING_MODE_PLANARFIGURE + MASKING_MODE_IMAGE = 1, + MASKING_MODE_PLANARFIGURE = 2 }; typedef itk::Statistics::Histogram HistogramType; typedef HistogramType::ConstIterator HistogramConstIteratorType; + /** + TODO DM: document + */ struct Statistics { int Label; - unsigned int N; - double Min; - double Max; - double Mean; - double Median; - double Variance; - double Sigma; - double RMS; + unsigned int N; //< number of voxels + double Min; //< mimimum value + double Max; //< maximum value + double Mean; //< mean value + double Median; //< median value + double Variance; //< variance of values // TODO DM: remove, was never filled with values ; check if any calling code within MITK used this member! + double Sigma; //< standard deviation of values (== square root of variance) + double RMS; //< root means square (TODO DM: check mesning) + double HotspotMin; //< mimimum value inside hotspot + double HotspotMax; //< maximum value inside hotspot + double HotspotMean; //< mean value inside hotspot + double HotspotSigma; //< standard deviation of values inside hotspot + //TODO DM: where is variance? does not make much sense, but should be consistent with usual statistics + //TODO DM: same goes for N + //TODO DM: same goes for RMS + double HotspotPeak; //< TODO DM: should this not replace "mean" the two values could be irritating vnl_vector< int > MinIndex; vnl_vector< int > MaxIndex; + vnl_vector HotspotMaxIndex; + vnl_vector HotspotMinIndex; + vnl_vector HotspotIndex; //< TODO DM: couldn't this be named "hotspot index"? We need to clear naming of hotspotmean, hotspotpeakindex, and hotspotpeak - void Reset() + // TODO DM: make this struct a real class and put this into a constructor + void Reset() // TODO DM: move to .cpp file (mitk::ImageStatisticsCalculator::Statistics::Reset() {...}) { Label = 0; N = 0; Min = 0.0; Max = 0.0; Mean = 0.0; Median = 0.0; Variance = 0.0; Sigma = 0.0; RMS = 0.0; + HotspotMin = 0.0; + HotspotMax = 0.0; + HotspotMean = 0.0; + HotspotPeak = 0.0; + HotspotSigma = 0.0; // TODO DM: also reset index values! Check that everything is initialized } }; + struct MinMaxIndex // TODO DM: why this structure? could at least be private + { + double Max; + double Min; + vnl_vector MaxIndex; + vnl_vector MinIndex; + }; + typedef std::vector< HistogramType::ConstPointer > HistogramContainer; typedef std::vector< Statistics > StatisticsContainer; - mitkClassMacro( ImageStatisticsCalculator, itk::Object ); itkNewMacro( ImageStatisticsCalculator ); /** \brief Set image from which to compute statistics. */ void SetImage( const mitk::Image *image ); /** \brief Set image for masking. */ void SetImageMask( const mitk::Image *imageMask ); /** \brief Set planar figure for masking. */ void SetPlanarFigure( mitk::PlanarFigure *planarFigure ); /** \brief Set/Get operation mode for masking */ void SetMaskingMode( unsigned int mode ); /** \brief Set/Get operation mode for masking */ itkGetMacro( MaskingMode, unsigned int ); /** \brief Set/Get operation mode for masking */ void SetMaskingModeToNone(); /** \brief Set/Get operation mode for masking */ void SetMaskingModeToImage(); /** \brief Set/Get operation mode for masking */ void SetMaskingModeToPlanarFigure(); /** \brief Set a pixel value for pixels that will be ignored in the statistics */ void SetIgnorePixelValue(double value); /** \brief Get the pixel value for pixels that will be ignored in the statistics */ double GetIgnorePixelValue(); - /** \brief Set wether a pixel value should be ignored in the statistics */ + /** \brief Set whether a pixel value should be ignored in the statistics */ void SetDoIgnorePixelValue(bool doit); - /** \brief Get wether a pixel value will be ignored in the statistics */ + /** \brief Get whether a pixel value will be ignored in the statistics */ bool GetDoIgnorePixelValue(); + /** \brief Sets the radius for the hotspot */ + void SetHotspotRadius (double hotspotRadiusInMM); // TODO in mm + + /** \brief Returns the radius of the hotspot */ + double GetHotspotRadius(); // TODO in mm + + /** \brief Sets whether the hotspot should be calculated */ + void SetCalculateHotspot(bool calculateHotspot); + + /** \brief Returns true whether the hotspot should be calculated, otherwise false */ + bool IsHotspotCalculated(); + /** \brief Compute statistics (together with histogram) for the current * masking mode. * * Computation is not executed if statistics is already up to date. In this * case, false is returned; otherwise, true.*/ virtual bool ComputeStatistics( unsigned int timeStep = 0 ); /** \brief Retrieve the histogram depending on the current masking mode. * * \param label The label for which to retrieve the histogram in multi-label situations (ascending order). */ const HistogramType *GetHistogram( unsigned int timeStep = 0, unsigned int label = 0 ) const; /** \brief Retrieve the histogram depending on the current masking mode (for all image labels. */ const HistogramContainer &GetHistogramVector( unsigned int timeStep = 0 ) const; /** \brief Retrieve statistics depending on the current masking mode. * * \param label The label for which to retrieve the statistics in multi-label situations (ascending order). */ const Statistics &GetStatistics( unsigned int timeStep = 0, unsigned int label = 0 ) const; + /** \brief Retrieve statistics depending on the current masking mode (for all image labels). */ const StatisticsContainer &GetStatisticsVector( unsigned int timeStep = 0 ) const; protected: typedef std::vector< HistogramContainer > HistogramVector; typedef std::vector< StatisticsContainer > StatisticsVector; typedef std::vector< itk::TimeStamp > TimeStampVectorType; typedef std::vector< bool > BoolVectorType; - - typedef itk::Image< unsigned short, 3 > MaskImage3DType; typedef itk::Image< unsigned short, 2 > MaskImage2DType; ImageStatisticsCalculator(); virtual ~ImageStatisticsCalculator(); /** \brief Depending on the masking mode, the image and mask from which to * calculate statistics is extracted from the original input image and mask * data. * * For example, a when using a PlanarFigure as mask, the 2D image slice * corresponding to the PlanarFigure will be extracted from the original * image. If masking is disabled, the original image is simply passed * through. */ void ExtractImageAndMask( unsigned int timeStep = 0 ); /** \brief If the passed vector matches any of the three principal axes * of the passed geometry, the ínteger value corresponding to the axis * is set and true is returned. */ bool GetPrincipalAxis( const Geometry3D *geometry, Vector3D vector, unsigned int &axis ); template < typename TPixel, unsigned int VImageDimension > void InternalCalculateStatisticsUnmasked( const itk::Image< TPixel, VImageDimension > *image, StatisticsContainer* statisticsContainer, HistogramContainer *histogramContainer ); template < typename TPixel, unsigned int VImageDimension > void InternalCalculateStatisticsMasked( const itk::Image< TPixel, VImageDimension > *image, itk::Image< unsigned short, VImageDimension > *maskImage, StatisticsContainer* statisticsContainer, HistogramContainer* histogramContainer ); template < typename TPixel, unsigned int VImageDimension > void InternalCalculateMaskFromPlanarFigure( const itk::Image< TPixel, VImageDimension > *image, unsigned int axis ); template < typename TPixel, unsigned int VImageDimension > void InternalMaskIgnoredPixels( const itk::Image< TPixel, VImageDimension > *image, itk::Image< unsigned short, VImageDimension > *maskImage ); + /** \brief Calculates minimum, maximum, mean value and their + * corresponding indices in a given ROI. As input the function + * needs an image and a mask. It returns a MinMaxIndex object. */ + template + MinMaxIndex CalculateMinMaxIndex( + const itk::Image *inputImage, + itk::Image *maskImage); + + /** \brief Calculates the hotspot statistics within a given + * ROI. As input the function needs an image, a mask which + * represents the ROI and a radius which defines the size of + * the sphere. The function returns a Statistics object. */ + template < typename TPixel, unsigned int VImageDimension> + Statistics CalculateHotspotStatistics( + const itk::Image *inputImage, + itk::Image *maskImage, + double radiusInMM); + /** Connection from ITK to VTK */ template void ConnectPipelines(ITK_Exporter exporter, vtkSmartPointer importer) { importer->SetUpdateInformationCallback(exporter->GetUpdateInformationCallback()); importer->SetPipelineModifiedCallback(exporter->GetPipelineModifiedCallback()); importer->SetWholeExtentCallback(exporter->GetWholeExtentCallback()); importer->SetSpacingCallback(exporter->GetSpacingCallback()); importer->SetOriginCallback(exporter->GetOriginCallback()); importer->SetScalarTypeCallback(exporter->GetScalarTypeCallback()); importer->SetNumberOfComponentsCallback(exporter->GetNumberOfComponentsCallback()); importer->SetPropagateUpdateExtentCallback(exporter->GetPropagateUpdateExtentCallback()); importer->SetUpdateDataCallback(exporter->GetUpdateDataCallback()); importer->SetDataExtentCallback(exporter->GetDataExtentCallback()); importer->SetBufferPointerCallback(exporter->GetBufferPointerCallback()); importer->SetCallbackUserData(exporter->GetCallbackUserData()); } /** Connection from VTK to ITK */ template void ConnectPipelines(vtkSmartPointer exporter, ITK_Importer importer) { importer->SetUpdateInformationCallback(exporter->GetUpdateInformationCallback()); importer->SetPipelineModifiedCallback(exporter->GetPipelineModifiedCallback()); importer->SetWholeExtentCallback(exporter->GetWholeExtentCallback()); importer->SetSpacingCallback(exporter->GetSpacingCallback()); importer->SetOriginCallback(exporter->GetOriginCallback()); importer->SetScalarTypeCallback(exporter->GetScalarTypeCallback()); importer->SetNumberOfComponentsCallback(exporter->GetNumberOfComponentsCallback()); importer->SetPropagateUpdateExtentCallback(exporter->GetPropagateUpdateExtentCallback()); importer->SetUpdateDataCallback(exporter->GetUpdateDataCallback()); importer->SetDataExtentCallback(exporter->GetDataExtentCallback()); importer->SetBufferPointerCallback(exporter->GetBufferPointerCallback()); importer->SetCallbackUserData(exporter->GetCallbackUserData()); } void UnmaskedStatisticsProgressUpdate(); void MaskedStatisticsProgressUpdate(); + template + itk::SmartPointer< itk::Image > + GenerateHotspotSearchConvolutionMask(double spacing[VImageDimension], double radiusInMM); /** m_Image contains the input image (e.g. 2D, 3D, 3D+t)*/ mitk::Image::ConstPointer m_Image; mitk::Image::ConstPointer m_ImageMask; mitk::PlanarFigure::Pointer m_PlanarFigure; HistogramVector m_ImageHistogramVector; HistogramVector m_MaskedImageHistogramVector; HistogramVector m_PlanarFigureHistogramVector; HistogramType::Pointer m_EmptyHistogram; HistogramContainer m_EmptyHistogramContainer; StatisticsVector m_ImageStatisticsVector; StatisticsVector m_MaskedImageStatisticsVector; StatisticsVector m_PlanarFigureStatisticsVector; + StatisticsVector m_MaskedImageHotspotStatisticsVector; Statistics m_EmptyStatistics; StatisticsContainer m_EmptyStatisticsContainer; unsigned int m_MaskingMode; bool m_MaskingModeChanged; /** m_InternalImage contains a image volume at one time step (e.g. 2D, 3D)*/ mitk::Image::ConstPointer m_InternalImage; MaskImage3DType::Pointer m_InternalImageMask3D; MaskImage2DType::Pointer m_InternalImageMask2D; TimeStampVectorType m_ImageStatisticsTimeStampVector; TimeStampVectorType m_MaskedImageStatisticsTimeStampVector; TimeStampVectorType m_PlanarFigureStatisticsTimeStampVector; BoolVectorType m_ImageStatisticsCalculationTriggerVector; BoolVectorType m_MaskedImageStatisticsCalculationTriggerVector; BoolVectorType m_PlanarFigureStatisticsCalculationTriggerVector; double m_IgnorePixelValue; bool m_DoIgnorePixelValue; bool m_IgnorePixelValueChanged; + double m_HotspotRadiusInMM; + bool m_CalculateHotspot; + unsigned int m_PlanarFigureAxis; // Normal axis for PlanarFigure unsigned int m_PlanarFigureSlice; // Slice which contains PlanarFigure int m_PlanarFigureCoordinate0; // First plane-axis for PlanarFigure int m_PlanarFigureCoordinate1; // Second plane-axis for PlanarFigure }; -} +} // namespace -#endif // #define _MITK_IMAGESTATISTICSCALCULATOR_H +#endif diff --git a/Plugins/org.mitk.gui.qt.measurementtoolbox/src/internal/QmitkImageStatisticsCalculationThread.cpp b/Plugins/org.mitk.gui.qt.measurementtoolbox/src/internal/QmitkImageStatisticsCalculationThread.cpp index 3014fee6ee..da6ffdab75 100644 --- a/Plugins/org.mitk.gui.qt.measurementtoolbox/src/internal/QmitkImageStatisticsCalculationThread.cpp +++ b/Plugins/org.mitk.gui.qt.measurementtoolbox/src/internal/QmitkImageStatisticsCalculationThread.cpp @@ -1,181 +1,182 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "QmitkImageStatisticsCalculationThread.h" //QT headers #include #include QmitkImageStatisticsCalculationThread::QmitkImageStatisticsCalculationThread():QThread(), m_StatisticsImage(NULL), m_BinaryMask(NULL), m_PlanarFigureMask(NULL), m_TimeStep(0), m_IgnoreZeros(false), m_CalculationSuccessful(false), m_StatisticChanged(false) { } QmitkImageStatisticsCalculationThread::~QmitkImageStatisticsCalculationThread() { } void QmitkImageStatisticsCalculationThread::Initialize( mitk::Image::Pointer image, mitk::Image::Pointer binaryImage, mitk::PlanarFigure::Pointer planarFig ) { // reset old values if( this->m_StatisticsImage.IsNotNull() ) this->m_StatisticsImage = 0; if( this->m_BinaryMask.IsNotNull() ) this->m_BinaryMask = 0; if( this->m_PlanarFigureMask.IsNotNull()) this->m_PlanarFigureMask = 0; // set new values if passed in if(image.IsNotNull()) this->m_StatisticsImage = image->Clone(); if(binaryImage.IsNotNull()) this->m_BinaryMask = binaryImage->Clone(); if(planarFig.IsNotNull()) this->m_PlanarFigureMask = dynamic_cast(planarFig.GetPointer()); // once clone methods for planar figures are implemented, copy the data here! } void QmitkImageStatisticsCalculationThread::SetTimeStep( int times ) { this->m_TimeStep = times; } int QmitkImageStatisticsCalculationThread::GetTimeStep() { return this->m_TimeStep; } mitk::ImageStatisticsCalculator::Statistics QmitkImageStatisticsCalculationThread::GetStatisticsData() { return this->m_StatisticsStruct; } mitk::Image::Pointer QmitkImageStatisticsCalculationThread::GetStatisticsImage() { return this->m_StatisticsImage; } void QmitkImageStatisticsCalculationThread::SetIgnoreZeroValueVoxel(bool _arg) { this->m_IgnoreZeros = _arg; } bool QmitkImageStatisticsCalculationThread::GetIgnoreZeroValueVoxel() { return this->m_IgnoreZeros; } std::string QmitkImageStatisticsCalculationThread::GetLastErrorMessage() { return m_message; } QmitkImageStatisticsCalculationThread::HistogramType::Pointer QmitkImageStatisticsCalculationThread::GetTimeStepHistogram() { return this->m_TimeStepHistogram; } bool QmitkImageStatisticsCalculationThread::GetStatisticsChangedFlag() { return m_StatisticChanged; } bool QmitkImageStatisticsCalculationThread::GetStatisticsUpdateSuccessFlag() { return m_CalculationSuccessful; } void QmitkImageStatisticsCalculationThread::run() { bool statisticCalculationSuccessful = true; mitk::ImageStatisticsCalculator::Pointer calculator = mitk::ImageStatisticsCalculator::New(); if(this->m_StatisticsImage.IsNotNull()) { calculator->SetImage(m_StatisticsImage); calculator->SetMaskingModeToNone(); } else { statisticCalculationSuccessful = false; } // Bug 13416 : The ImageStatistics::SetImageMask() method can throw exceptions, i.e. when the dimensionality // of the masked and input image differ, we need to catch them and mark the calculation as failed // the same holds for the ::SetPlanarFigure() try { if(this->m_BinaryMask.IsNotNull()) { calculator->SetImageMask(m_BinaryMask); calculator->SetMaskingModeToImage(); } if(this->m_PlanarFigureMask.IsNotNull()) { calculator->SetPlanarFigure(m_PlanarFigureMask); calculator->SetMaskingModeToPlanarFigure(); } } catch( const itk::ExceptionObject& e) { MITK_ERROR << "ITK Exception:" << e.what(); statisticCalculationSuccessful = false; } bool statisticChanged = false; calculator->SetDoIgnorePixelValue(this->m_IgnoreZeros); calculator->SetIgnorePixelValue(0); try { + //calculator->SetCalculateHotspot(true); statisticChanged = calculator->ComputeStatistics(m_TimeStep); } catch ( mitk::Exception& e) { //m_message = e.GetDescription(); MITK_ERROR<< "MITK Exception: " << e.what(); statisticCalculationSuccessful = false; } catch ( const std::runtime_error &e ) { //m_message = "Failure: " + std::string(e.what()); MITK_ERROR<< "Runtime Exception: " << e.what(); statisticCalculationSuccessful = false; } catch ( const std::exception &e ) { //m_message = "Failure: " + std::string(e.what()); MITK_ERROR<< "Standard Exception: " << e.what(); statisticCalculationSuccessful = false; } this->m_StatisticChanged = statisticChanged; this->m_CalculationSuccessful = statisticCalculationSuccessful; if(statisticCalculationSuccessful) { this->m_StatisticsStruct = calculator->GetStatistics(m_TimeStep); if(this->m_TimeStepHistogram.IsNotNull()) { this->m_TimeStepHistogram = NULL; } this->m_TimeStepHistogram = (HistogramType*) calculator->GetHistogram(m_TimeStep); } } diff --git a/Plugins/org.mitk.gui.qt.measurementtoolbox/src/internal/QmitkImageStatisticsView.cpp b/Plugins/org.mitk.gui.qt.measurementtoolbox/src/internal/QmitkImageStatisticsView.cpp index 263d732b64..a00c945b3e 100644 --- a/Plugins/org.mitk.gui.qt.measurementtoolbox/src/internal/QmitkImageStatisticsView.cpp +++ b/Plugins/org.mitk.gui.qt.measurementtoolbox/src/internal/QmitkImageStatisticsView.cpp @@ -1,702 +1,743 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "QmitkImageStatisticsView.h" // Qt includes #include // berry includes #include // mitk includes #include "mitkNodePredicateDataType.h" #include "mitkPlanarFigureInteractor.h" // itk includes #include "itksys/SystemTools.hxx" #include #include const std::string QmitkImageStatisticsView::VIEW_ID = "org.mitk.views.imagestatistics"; QmitkImageStatisticsView::QmitkImageStatisticsView(QObject* /*parent*/, const char* /*name*/) : m_Controls( NULL ), m_TimeStepperAdapter( NULL ), m_SelectedImage( NULL ), m_SelectedImageMask( NULL ), m_SelectedPlanarFigure( NULL ), m_ImageObserverTag( -1 ), m_ImageMaskObserverTag( -1 ), m_PlanarFigureObserverTag( -1 ), m_CurrentStatisticsValid( false ), m_StatisticsUpdatePending( false ), m_DataNodeSelectionChanged ( false ), m_Visible(false) { this->m_CalculationThread = new QmitkImageStatisticsCalculationThread; } QmitkImageStatisticsView::~QmitkImageStatisticsView() { if ( m_SelectedImage != NULL ) m_SelectedImage->RemoveObserver( m_ImageObserverTag ); if ( m_SelectedImageMask != NULL ) m_SelectedImageMask->RemoveObserver( m_ImageMaskObserverTag ); if ( m_SelectedPlanarFigure != NULL ) m_SelectedPlanarFigure->RemoveObserver( m_PlanarFigureObserverTag ); while(this->m_CalculationThread->isRunning()) // wait until thread has finished { itksys::SystemTools::Delay(100); } delete this->m_CalculationThread; } void QmitkImageStatisticsView::CreateQtPartControl(QWidget *parent) { if (m_Controls == NULL) { m_Controls = new Ui::QmitkImageStatisticsViewControls; m_Controls->setupUi(parent); this->CreateConnections(); m_Controls->m_ErrorMessageLabel->hide(); m_Controls->m_StatisticsWidgetStack->setCurrentIndex( 0 ); m_Controls->m_LineProfileWidget->SetPathModeToPlanarFigure(); } } void QmitkImageStatisticsView::CreateConnections() { if ( m_Controls ) { connect( (QObject*)(this->m_Controls->m_ButtonCopyHistogramToClipboard), SIGNAL(clicked()),(QObject*) this, SLOT(OnClipboardHistogramButtonClicked()) ); connect( (QObject*)(this->m_Controls->m_ButtonCopyStatisticsToClipboard), SIGNAL(clicked()),(QObject*) this, SLOT(OnClipboardStatisticsButtonClicked()) ); connect( (QObject*)(this->m_Controls->m_IgnoreZerosCheckbox), SIGNAL(clicked()),(QObject*) this, SLOT(OnIgnoreZerosCheckboxClicked()) ); connect( (QObject*) this->m_CalculationThread, SIGNAL(finished()),this, SLOT( OnThreadedStatisticsCalculationEnds()),Qt::QueuedConnection); connect( (QObject*) this, SIGNAL(StatisticsUpdate()),this, SLOT( RequestStatisticsUpdate()), Qt::QueuedConnection); connect( (QObject*) this->m_Controls->m_StatisticsTable, SIGNAL(cellDoubleClicked(int,int)),this, SLOT( JumpToCoordinates(int,int)) ); connect( (QObject*) (this->m_Controls->m_barRadioButton), SIGNAL(clicked()), (QObject*) (this->m_Controls->m_JSHistogram), SLOT(OnBarRadioButtonSelected())); connect( (QObject*) (this->m_Controls->m_lineRadioButton), SIGNAL(clicked()), (QObject*) (this->m_Controls->m_JSHistogram), SLOT(OnLineRadioButtonSelected())); } } void QmitkImageStatisticsView::JumpToCoordinates(int row ,int col) { mitk::Point3D world; if (row==4) world = m_WorldMin; else if (row==3) world = m_WorldMax; else return; mitk::IRenderWindowPart* part = this->GetRenderWindowPart(); if (part) { part->GetQmitkRenderWindow("axial")->GetSliceNavigationController()->SelectSliceByPoint(world); part->GetQmitkRenderWindow("sagittal")->GetSliceNavigationController()->SelectSliceByPoint(world); part->GetQmitkRenderWindow("coronal")->GetSliceNavigationController()->SelectSliceByPoint(world); } } void QmitkImageStatisticsView::OnIgnoreZerosCheckboxClicked() { emit StatisticsUpdate(); } void QmitkImageStatisticsView::OnClipboardHistogramButtonClicked() { if ( m_CurrentStatisticsValid ) { typedef mitk::ImageStatisticsCalculator::HistogramType HistogramType; const HistogramType *histogram = this->m_CalculationThread->GetTimeStepHistogram().GetPointer(); QString clipboard( "Measurement \t Frequency\n" ); for ( HistogramType::ConstIterator it = histogram->Begin(); it != histogram->End(); ++it ) { clipboard = clipboard.append( "%L1 \t %L2\n" ) .arg( it.GetMeasurementVector()[0], 0, 'f', 2 ) .arg( it.GetFrequency() ); } QApplication::clipboard()->setText( clipboard, QClipboard::Clipboard ); } else { QApplication::clipboard()->clear(); } } void QmitkImageStatisticsView::OnClipboardStatisticsButtonClicked() { if ( this->m_CurrentStatisticsValid ) { const mitk::ImageStatisticsCalculator::Statistics &statistics = this->m_CalculationThread->GetStatisticsData(); // Copy statistics to clipboard ("%Ln" will use the default locale for // number formatting) QString clipboard( "Mean \t StdDev \t RMS \t Max \t Min \t N \t V (mm³)\n" ); clipboard = clipboard.append( "%L1 \t %L2 \t %L3 \t %L4 \t %L5 \t %L6 \t %L7" ) .arg( statistics.Mean, 0, 'f', 10 ) .arg( statistics.Sigma, 0, 'f', 10 ) .arg( statistics.RMS, 0, 'f', 10 ) .arg( statistics.Max, 0, 'f', 10 ) .arg( statistics.Min, 0, 'f', 10 ) .arg( statistics.N ) .arg( m_Controls->m_StatisticsTable->item( 0, 6 )->text() ); QApplication::clipboard()->setText( clipboard, QClipboard::Clipboard ); } else { QApplication::clipboard()->clear(); } } void QmitkImageStatisticsView::OnSelectionChanged( berry::IWorkbenchPart::Pointer /*part*/, const QList &selectedNodes ) { if (this->m_Visible) { this->SelectionChanged( selectedNodes ); } else { this->m_DataNodeSelectionChanged = true; } } void QmitkImageStatisticsView::SelectionChanged(const QList &selectedNodes) { if( this->m_StatisticsUpdatePending ) { this->m_DataNodeSelectionChanged = true; return; // not ready for new data now! } if (selectedNodes.size() == this->m_SelectedDataNodes.size()) { int i = 0; for (; i < selectedNodes.size(); ++i) { if (selectedNodes.at(i) != this->m_SelectedDataNodes.at(i)) { break; } } // node selection did not change if (i == selectedNodes.size()) return; } this->ReinitData(); if (!selectedNodes.size()) { m_Controls->m_JSHistogram->ClearHistogram(); m_Controls->m_lineRadioButton->setEnabled(true); m_Controls->m_barRadioButton->setEnabled(true); m_Controls->m_InfoLabel->setText(QString("")); } if(selectedNodes.size() == 1 || selectedNodes.size() == 2) { bool isBinary = false; selectedNodes.value(0)->GetBoolProperty("binary",isBinary); if(isBinary) { m_Controls->m_JSHistogram->ClearHistogram(); m_Controls->m_lineRadioButton->setEnabled(true); m_Controls->m_barRadioButton->setEnabled(true); m_Controls->m_InfoLabel->setText(QString("")); } for (int i= 0; i< selectedNodes.size(); ++i) { this->m_SelectedDataNodes.push_back(selectedNodes.at(i)); } this->m_DataNodeSelectionChanged = false; this->m_Controls->m_ErrorMessageLabel->setText( "" ); this->m_Controls->m_ErrorMessageLabel->hide(); emit StatisticsUpdate(); } else { this->m_DataNodeSelectionChanged = false; } } void QmitkImageStatisticsView::ReinitData() { while( this->m_CalculationThread->isRunning()) // wait until thread has finished { itksys::SystemTools::Delay(100); } if(this->m_SelectedImage != NULL) { this->m_SelectedImage->RemoveObserver( this->m_ImageObserverTag); this->m_SelectedImage = NULL; } if(this->m_SelectedImageMask != NULL) { this->m_SelectedImageMask->RemoveObserver( this->m_ImageMaskObserverTag); this->m_SelectedImageMask = NULL; } if(this->m_SelectedPlanarFigure != NULL) { this->m_SelectedPlanarFigure->RemoveObserver( this->m_PlanarFigureObserverTag); this->m_SelectedPlanarFigure = NULL; } this->m_SelectedDataNodes.clear(); this->m_StatisticsUpdatePending = false; m_Controls->m_ErrorMessageLabel->setText( "" ); m_Controls->m_ErrorMessageLabel->hide(); this->InvalidateStatisticsTableView(); m_Controls->m_StatisticsWidgetStack->setCurrentIndex( 0 ); } void QmitkImageStatisticsView::OnThreadedStatisticsCalculationEnds() { std::stringstream message; message << ""; m_Controls->m_ErrorMessageLabel->setText( message.str().c_str() ); m_Controls->m_ErrorMessageLabel->hide(); this->WriteStatisticsToGUI(); } void QmitkImageStatisticsView::UpdateStatistics() { mitk::IRenderWindowPart* renderPart = this->GetRenderWindowPart(); if ( renderPart == NULL ) { this->m_StatisticsUpdatePending = false; return; } m_WorldMin.Fill(-1); m_WorldMax.Fill(-1); // classify selected nodes mitk::NodePredicateDataType::Pointer imagePredicate = mitk::NodePredicateDataType::New("Image"); std::string maskName = std::string(); std::string maskType = std::string(); unsigned int maskDimension = 0; // reset data from last run ITKCommandType::Pointer changeListener = ITKCommandType::New(); changeListener->SetCallbackFunction( this, &QmitkImageStatisticsView::SelectedDataModified ); mitk::DataNode::Pointer planarFigureNode; for( int i= 0 ; i < this->m_SelectedDataNodes.size(); ++i) { mitk::PlanarFigure::Pointer planarFig = dynamic_cast(this->m_SelectedDataNodes.at(i)->GetData()); if( imagePredicate->CheckNode(this->m_SelectedDataNodes.at(i)) ) { bool isMask = false; this->m_SelectedDataNodes.at(i)->GetPropertyValue("binary", isMask); if( this->m_SelectedImageMask == NULL && isMask) { this->m_SelectedImageMask = dynamic_cast(this->m_SelectedDataNodes.at(i)->GetData()); this->m_ImageMaskObserverTag = this->m_SelectedImageMask->AddObserver(itk::ModifiedEvent(), changeListener); maskName = this->m_SelectedDataNodes.at(i)->GetName(); maskType = m_SelectedImageMask->GetNameOfClass(); maskDimension = 3; } else if( !isMask ) { if(this->m_SelectedImage == NULL) { this->m_SelectedImage = static_cast(this->m_SelectedDataNodes.at(i)->GetData()); this->m_ImageObserverTag = this->m_SelectedImage->AddObserver(itk::ModifiedEvent(), changeListener); } } } else if (planarFig.IsNotNull()) { if(this->m_SelectedPlanarFigure == NULL) { this->m_SelectedPlanarFigure = planarFig; this->m_PlanarFigureObserverTag = this->m_SelectedPlanarFigure->AddObserver(mitk::EndInteractionPlanarFigureEvent(), changeListener); maskName = this->m_SelectedDataNodes.at(i)->GetName(); maskType = this->m_SelectedPlanarFigure->GetNameOfClass(); maskDimension = 2; planarFigureNode = m_SelectedDataNodes.at(i); } } else { std::stringstream message; message << "" << "Invalid data node type!" << ""; m_Controls->m_ErrorMessageLabel->setText( message.str().c_str() ); m_Controls->m_ErrorMessageLabel->show(); } } if(maskName == "") { maskName = "None"; maskType = ""; maskDimension = 0; } if (m_SelectedPlanarFigure != NULL && m_SelectedImage == NULL) { mitk::DataStorage::SetOfObjects::ConstPointer parentSet = this->GetDataStorage()->GetSources(planarFigureNode); for (int i=0; iSize(); i++) { mitk::DataNode::Pointer node = parentSet->ElementAt(i); if( imagePredicate->CheckNode(node) ) { bool isMask = false; node->GetPropertyValue("binary", isMask); if( !isMask ) { if(this->m_SelectedImage == NULL) { this->m_SelectedImage = static_cast(node->GetData()); this->m_ImageObserverTag = this->m_SelectedImage->AddObserver(itk::ModifiedEvent(), changeListener); } } } } } unsigned int timeStep = renderPart->GetTimeNavigationController()->GetTime()->GetPos(); if ( m_SelectedImage != NULL && m_SelectedImage->IsInitialized()) { // Check if a the selected image is a multi-channel image. If yes, statistics // cannot be calculated currently. if ( m_SelectedImage->GetPixelType().GetNumberOfComponents() > 1 ) { std::stringstream message; message << "Multi-component images not supported."; m_Controls->m_ErrorMessageLabel->setText( message.str().c_str() ); m_Controls->m_ErrorMessageLabel->show(); this->InvalidateStatisticsTableView(); m_Controls->m_StatisticsWidgetStack->setCurrentIndex( 0 ); m_Controls->m_JSHistogram->ClearHistogram(); m_CurrentStatisticsValid = false; this->m_StatisticsUpdatePending = false; m_Controls->m_lineRadioButton->setEnabled(true); m_Controls->m_barRadioButton->setEnabled(true); m_Controls->m_InfoLabel->setText(QString("")); return; } std::stringstream maskLabel; maskLabel << maskName; if ( maskDimension > 0 ) { maskLabel << " [" << maskDimension << "D " << maskType << "]"; } m_Controls->m_SelectedMaskLabel->setText( maskLabel.str().c_str() ); // check time step validity if(m_SelectedImage->GetDimension() <= 3 && timeStep > m_SelectedImage->GetDimension(3)-1) { timeStep = m_SelectedImage->GetDimension(3)-1; } //// initialize thread and trigger it this->m_CalculationThread->SetIgnoreZeroValueVoxel( m_Controls->m_IgnoreZerosCheckbox->isChecked() ); this->m_CalculationThread->Initialize( m_SelectedImage, m_SelectedImageMask, m_SelectedPlanarFigure ); this->m_CalculationThread->SetTimeStep( timeStep ); std::stringstream message; message << "Calculating statistics..."; m_Controls->m_ErrorMessageLabel->setText( message.str().c_str() ); m_Controls->m_ErrorMessageLabel->show(); try { // Compute statistics this->m_CalculationThread->start(); } catch ( const mitk::Exception& e) { std::stringstream message; message << "" << e.GetDescription() << ""; m_Controls->m_ErrorMessageLabel->setText( message.str().c_str() ); m_Controls->m_ErrorMessageLabel->show(); this->m_StatisticsUpdatePending = false; } catch ( const std::runtime_error &e ) { // In case of exception, print error message on GUI std::stringstream message; message << "" << e.what() << ""; m_Controls->m_ErrorMessageLabel->setText( message.str().c_str() ); m_Controls->m_ErrorMessageLabel->show(); this->m_StatisticsUpdatePending = false; } catch ( const std::exception &e ) { MITK_ERROR << "Caught exception: " << e.what(); // In case of exception, print error message on GUI std::stringstream message; message << "Error! Unequal Dimensions of Image and Segmentation. No recompute possible "; m_Controls->m_ErrorMessageLabel->setText( message.str().c_str() ); m_Controls->m_ErrorMessageLabel->show(); this->m_StatisticsUpdatePending = false; } } else { this->m_StatisticsUpdatePending = false; } } void QmitkImageStatisticsView::SelectedDataModified() { if( !m_StatisticsUpdatePending ) { emit StatisticsUpdate(); } } void QmitkImageStatisticsView::NodeRemoved(const mitk::DataNode *node) { while(this->m_CalculationThread->isRunning()) // wait until thread has finished { itksys::SystemTools::Delay(100); } if (node->GetData() == m_SelectedImage) { m_SelectedImage = NULL; } } void QmitkImageStatisticsView::RequestStatisticsUpdate() { if ( !m_StatisticsUpdatePending ) { if(this->m_DataNodeSelectionChanged) { this->SelectionChanged(this->GetCurrentSelection()); } else { this->m_StatisticsUpdatePending = true; this->UpdateStatistics(); } } if (this->GetRenderWindowPart()) this->GetRenderWindowPart()->RequestUpdate(); } void QmitkImageStatisticsView::WriteStatisticsToGUI() { m_Controls->m_lineRadioButton->setEnabled(true); m_Controls->m_barRadioButton->setEnabled(true); m_Controls->m_InfoLabel->setText(QString("")); if(m_DataNodeSelectionChanged) { this->m_StatisticsUpdatePending = false; this->RequestStatisticsUpdate(); return; // stop visualization of results and calculate statistics of new selection } if ( this->m_CalculationThread->GetStatisticsUpdateSuccessFlag()) { if ( this->m_CalculationThread->GetStatisticsChangedFlag() ) { // Do not show any error messages m_Controls->m_ErrorMessageLabel->hide(); m_CurrentStatisticsValid = true; } if (m_Controls->m_barRadioButton->isChecked()) { m_Controls->m_JSHistogram->OnBarRadioButtonSelected(); } m_Controls->m_StatisticsWidgetStack->setCurrentIndex( 0 ); m_Controls->m_JSHistogram->ComputeHistogram( this->m_CalculationThread->GetTimeStepHistogram().GetPointer() ); this->FillStatisticsTableView( this->m_CalculationThread->GetStatisticsData(), this->m_CalculationThread->GetStatisticsImage()); } else { m_Controls->m_SelectedMaskLabel->setText( "None" ); m_Controls->m_ErrorMessageLabel->setText( m_CalculationThread->GetLastErrorMessage().c_str() ); m_Controls->m_ErrorMessageLabel->show(); // Clear statistics and histogram this->InvalidateStatisticsTableView(); m_Controls->m_StatisticsWidgetStack->setCurrentIndex( 0 ); //m_Controls->m_JSHistogram->clearHistogram(); m_CurrentStatisticsValid = false; // If a (non-closed) PlanarFigure is selected, display a line profile widget if ( m_SelectedPlanarFigure != NULL ) { // check whether PlanarFigure is initialized const mitk::Geometry2D *planarFigureGeometry2D = m_SelectedPlanarFigure->GetGeometry2D(); if ( planarFigureGeometry2D == NULL ) { // Clear statistics, histogram, and GUI this->InvalidateStatisticsTableView(); m_Controls->m_StatisticsWidgetStack->setCurrentIndex( 0 ); m_Controls->m_JSHistogram->ClearHistogram(); m_CurrentStatisticsValid = false; m_Controls->m_ErrorMessageLabel->hide(); m_Controls->m_SelectedMaskLabel->setText( "None" ); this->m_StatisticsUpdatePending = false; m_Controls->m_lineRadioButton->setEnabled(true); m_Controls->m_barRadioButton->setEnabled(true); m_Controls->m_InfoLabel->setText(QString("")); return; } unsigned int timeStep = this->GetRenderWindowPart()->GetTimeNavigationController()->GetTime()->GetPos(); m_Controls->m_JSHistogram->SetImage(this->m_CalculationThread->GetStatisticsImage()); m_Controls->m_JSHistogram->SetPlanarFigure(m_SelectedPlanarFigure); m_Controls->m_JSHistogram->ComputeIntensityProfile(timeStep); m_Controls->m_lineRadioButton->setEnabled(false); m_Controls->m_barRadioButton->setEnabled(false); std::stringstream message; message << "Only linegraph available for an intesityprofile!"; m_Controls->m_InfoLabel->setText(message.str().c_str()); } } this->m_StatisticsUpdatePending = false; } void QmitkImageStatisticsView::FillStatisticsTableView( const mitk::ImageStatisticsCalculator::Statistics &s, const mitk::Image *image ) { if (s.MaxIndex.size()==3) { mitk::Point3D index; index[0] = s.MaxIndex[0]; index[1] = s.MaxIndex[1]; index[2] = s.MaxIndex[2]; m_SelectedImage->GetGeometry()->IndexToWorld(index, m_WorldMax); index[0] = s.MinIndex[0]; index[1] = s.MinIndex[1]; index[2] = s.MinIndex[2]; m_SelectedImage->GetGeometry()->IndexToWorld(index, m_WorldMin); } int decimals = 2; mitk::PixelType doublePix = mitk::MakeScalarPixelType< double >(); mitk::PixelType floatPix = mitk::MakeScalarPixelType< float >(); if (image->GetPixelType()==doublePix || image->GetPixelType()==floatPix) decimals = 5; this->m_Controls->m_StatisticsTable->setItem( 0, 0, new QTableWidgetItem( QString("%1").arg(s.Mean, 0, 'f', decimals) ) ); this->m_Controls->m_StatisticsTable->setItem( 0, 1, new QTableWidgetItem( QString("%1").arg(s.Sigma, 0, 'f', decimals) ) ); this->m_Controls->m_StatisticsTable->setItem( 0, 2, new QTableWidgetItem( QString("%1").arg(s.RMS, 0, 'f', decimals) ) ); QString max; max.append(QString("%1").arg(s.Max, 0, 'f', decimals)); max += " ("; for (int i=0; im_Controls->m_StatisticsTable->setItem( 0, 3, new QTableWidgetItem( max ) ); QString min; min.append(QString("%1").arg(s.Min, 0, 'f', decimals)); min += " ("; for (int i=0; im_Controls->m_StatisticsTable->setItem( 0, 4, new QTableWidgetItem( min ) ); this->m_Controls->m_StatisticsTable->setItem( 0, 5, new QTableWidgetItem( QString("%1").arg(s.N) ) ); const mitk::Geometry3D *geometry = image->GetGeometry(); if ( geometry != NULL ) { const mitk::Vector3D &spacing = image->GetGeometry()->GetSpacing(); double volume = spacing[0] * spacing[1] * spacing[2] * (double) s.N; this->m_Controls->m_StatisticsTable->setItem( 0, 6, new QTableWidgetItem( QString("%1").arg(volume, 0, 'f', decimals) ) ); } else { this->m_Controls->m_StatisticsTable->setItem( 0, 6, new QTableWidgetItem( "NA" ) ); } + + + QString hotspotPeak; hotspotPeak.append(QString("%1").arg(s.HotspotPeak, 0, 'f', decimals)); + hotspotPeak += " ("; + for (int i=0; im_Controls->m_StatisticsTable->setItem( 0, 7, new QTableWidgetItem( hotspotPeak ) ); + + + QString hotspotMax; hotspotMax.append(QString("%1").arg(s.HotspotMax, 0, 'f', decimals)); + hotspotMax += " ("; + for (int i=0; im_Controls->m_StatisticsTable->setItem( 0, 8, new QTableWidgetItem( hotspotMax ) ); + + + QString hotspotMin; hotspotMin.append(QString("%1").arg(s.HotspotMin, 0, 'f', decimals)); + hotspotMin += " ("; + for (int i=0; im_Controls->m_StatisticsTable->setItem( 0, 9, new QTableWidgetItem( hotspotMin ) ); + + } void QmitkImageStatisticsView::InvalidateStatisticsTableView() { for ( unsigned int i = 0; i < 7; ++i ) { this->m_Controls->m_StatisticsTable->setItem( 0, i, new QTableWidgetItem( "NA" ) ); } } void QmitkImageStatisticsView::Activated() { } void QmitkImageStatisticsView::Deactivated() { } void QmitkImageStatisticsView::Visible() { m_Visible = true; if (m_DataNodeSelectionChanged) { if (this->IsCurrentSelectionValid()) { this->SelectionChanged(this->GetCurrentSelection()); } else { this->SelectionChanged(this->GetDataManagerSelection()); } m_DataNodeSelectionChanged = false; } } void QmitkImageStatisticsView::Hidden() { m_Visible = false; } void QmitkImageStatisticsView::SetFocus() { } diff --git a/Plugins/org.mitk.gui.qt.measurementtoolbox/src/internal/QmitkImageStatisticsViewControls.ui b/Plugins/org.mitk.gui.qt.measurementtoolbox/src/internal/QmitkImageStatisticsViewControls.ui index 711c668e07..6a16df6cd2 100644 --- a/Plugins/org.mitk.gui.qt.measurementtoolbox/src/internal/QmitkImageStatisticsViewControls.ui +++ b/Plugins/org.mitk.gui.qt.measurementtoolbox/src/internal/QmitkImageStatisticsViewControls.ui @@ -1,477 +1,480 @@ QmitkImageStatisticsViewControls true 0 0 465 800 Form 0 0 Mask: None 2 Qt::Horizontal 40 20 color: rgb(255, 0, 0); Error Message Qt::AutoText Ignore zero-valued voxels Statistics 9 9 9 - - - 0 - 0 - - 100 175 - - - 16777215 - 170 - - Qt::ScrollBarAlwaysOff Qt::ScrollBarAsNeeded true true true Qt::DotLine true false false 80 true 80 false true true false 25 25 false false Mean StdDev RMS Max Min N V (mm³) + + + Hotspot Peak + + + + + Hotspot Max + + + + + Hotspot Min + + Component 1 0 0 Copy to Clipboard Qt::Horizontal 40 20 150 160 Histogram false 0 0 0 0 0 0 0 Copy to Clipboard Qt::Horizontal 40 20 0 0 0 50 16777215 16777215 Plot 0 20 411 31 QLayout::SetDefaultConstraint Qt::Horizontal QSizePolicy::Preferred 10 0 0 0 Use barchart true 0 0 0 0 Use linegraph Qt::Horizontal 40 20 Qt::Vertical 20 40 QmitkVtkHistogramWidget QWidget
QmitkVtkHistogramWidget.h
1 QmitkVtkLineProfileWidget QWidget
QmitkVtkLineProfileWidget.h
1 QmitkHistogramJSWidget QWidget
QmitkHistogramJSWidget.h
1