diff --git a/Modules/ImageStatistics/Testing/itkMultiGaussianImageSource.h b/Modules/ImageStatistics/Testing/itkMultiGaussianImageSource.h deleted file mode 100644 index b869d50bc8..0000000000 --- a/Modules/ImageStatistics/Testing/itkMultiGaussianImageSource.h +++ /dev/null @@ -1,219 +0,0 @@ -/*========================================================================= - * - * 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; - /** Set the number of gaussian function*/ - 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; - /** Set the number of steps to traverse the radius */ - 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(); - /** Calculates and set the index an the value of maximulm and minimum in the wanted sphere*/ - virtual void CalculateMaxAndMinInSphere(); - /** Get the index in the sphere with maximal value*/ - virtual const ItkVectorType GetMaxValueIndexInSphere() const; - /** Get the maximal value in the sphere*/ - virtual const double GetMaxValueInSphere() const; - /** Get the index in the sphere with minimal value*/ - virtual const ItkVectorType GetMinValueIndexInSphere() const; - /** Get the minimal value in the sphere*/ - virtual const double GetMinValueInSphere() const; - /** Set the region of interest */ - virtual void SetRegionOfInterest(ItkVectorType, ItkVectorType); - /** Optimize the mean value in the wanted sphere*/ - virtual void OptimizeMeanValue(); - /** Check whether the point is in the region of interest */ - virtual const bool IsInRegionOfInterest(unsigned int, unsigned int, unsigned int); - /** 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 - OutputImagePixelType m_MaxValueInSphere; //maximal value in the wanted sphere - ItkVectorType m_MaxValueIndexInSphere; //index of the maximal value in the wanted sphere - OutputImagePixelType m_MinValueInSphere; //minimal value in the wanted sphere - ItkVectorType m_MinValueIndexInSphere; //index of the minimal value in the wanted sphere - 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 - ItkVectorType m_RegionOfInterestMax; // maximal values for the coordinates in the region of interest - ItkVectorType m_RegionOfInterestMin; // minimal values for the coordinates in the region of interest - 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 deleted file mode 100644 index 9caaafe472..0000000000 --- a/Modules/ImageStatistics/Testing/itkMultiGaussianImageSource.hxx +++ /dev/null @@ -1,585 +0,0 @@ -/*========================================================================= - * - * 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 > -void -MultiGaussianImageSource< TOutputImage > -::SetRegionOfInterest( ItkVectorType roiMin, ItkVectorType roiMax ) -{ - m_RegionOfInterestMax.operator=(roiMax); - m_RegionOfInterestMin.operator=(roiMin); -} - - -template< class TOutputImage > -const double -MultiGaussianImageSource< TOutputImage > -::GetMaxMeanValue() const -{ - return m_MeanValue; -} - -template< class TOutputImage > -const double -MultiGaussianImageSource< TOutputImage > -::GetMaxValueInSphere() const -{ - return m_MaxValueInSphere; -} - -template< class TOutputImage > -const typename MultiGaussianImageSource< TOutputImage >::ItkVectorType -MultiGaussianImageSource< TOutputImage > -::GetMaxValueIndexInSphere() const -{ - return m_MaxValueIndexInSphere; -} - -template< class TOutputImage > -const double -MultiGaussianImageSource< TOutputImage > -::GetMinValueInSphere() const -{ - return m_MinValueInSphere; -} - -template< class TOutputImage > -const typename MultiGaussianImageSource< TOutputImage >::ItkVectorType -MultiGaussianImageSource< TOutputImage > -::GetMinValueIndexInSphere() const -{ - return m_MinValueIndexInSphere; -} -//---------------------------------------------------------------------------- -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< class TOutputImage > -const bool -MultiGaussianImageSource< TOutputImage > -::IsInRegionOfInterest(unsigned int ind0, unsigned int ind1, unsigned int ind2 ) -{ - bool isInROI = 0; - if(ind0 >= m_RegionOfInterestMin[0] && ind0 <=m_RegionOfInterestMax[0] && - ind1 >= m_RegionOfInterestMin[1] && ind1 <=m_RegionOfInterestMax[1] && - ind2 >= m_RegionOfInterestMin[2] && ind2 <=m_RegionOfInterestMax[2]) - { - isInROI = 1; - } - return isInROI; -} -//---------------------------------------------------------------------------- -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 region of interest 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) - { - if(IsInRegionOfInterest(index0, index1, index2)) - { - m_MaxValueInSphere = 0.0; - m_MinValueInSphere = 1000.0; - 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.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 ); - } - } - } - } - } -} - -template< typename TOutputImage > -void -MultiGaussianImageSource< TOutputImage > -::OptimizeMeanValue() -{ - int radiusStepNumberOptimized = m_RadiusStepNumber * 4; - int numberSummand = 0, angleStepNumberOverTwo; - double valueAdd, value, x, y, z, temp; - double riStep, fijStep, psikStep, ri, fij, psik; - double dist = itk::Math::pi * m_Radius / (2 * radiusStepNumberOptimized); - m_MeanValue = 0; - riStep = m_Radius / radiusStepNumberOptimized; - - int index0 = m_SphereMidpoint[0]; - int index1 = m_SphereMidpoint[1]; - int index2 = m_SphereMidpoint[2]; - ri = riStep; - value = MultiGaussianFunctionValueAtPoint(index0, index1, index2); - for(unsigned int i = 0; i < radiusStepNumberOptimized; ++i) - { - angleStepNumberOverTwo = static_cast( itk::Math::pi * ri / dist); - fij = 0.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; - } - m_MeanValue = value / numberSummand; - -} - -template< typename TOutputImage > -void -MultiGaussianImageSource< TOutputImage > -::CalculateMaxAndMinInSphere() -{ - double value; - m_MaxValueInSphere = 0.0; - m_MinValueInSphere = 100 * m_MeanValue; - int radInt = static_cast(m_Radius); - for(int index0 = -radInt; index0 <= radInt; ++index0 ){ - for(int index1 = -radInt; index1 <= radInt; ++index1 ){ - for(int index2 = -radInt; index2 <= radInt; ++index2 ){ - if( index0 * index0 + index1 * index1 + index2 * index2 < m_Radius * m_Radius ) - { - value = MultiGaussianFunctionValueAtPoint( m_SphereMidpoint[0] + index0, m_SphereMidpoint[1] + index1, m_SphereMidpoint[2] + index2); - if(m_MaxValueInSphere < value){ - m_MaxValueInSphere = value; - m_MaxValueIndexInSphere.SetNthComponent( 0, m_SphereMidpoint[0] + index0 ); - m_MaxValueIndexInSphere.SetNthComponent( 1, m_SphereMidpoint[1] + index1 ); - m_MaxValueIndexInSphere.SetNthComponent( 2, m_SphereMidpoint[2] + index2 ); - } - if(m_MinValueInSphere > value){ - m_MinValueInSphere = value; - m_MinValueIndexInSphere.SetNthComponent( 0, m_SphereMidpoint[0] + index0 ); - m_MinValueIndexInSphere.SetNthComponent( 1, m_SphereMidpoint[1] + index1 ); - m_MinValueIndexInSphere.SetNthComponent( 2, m_SphereMidpoint[2] + index2 ); - } - } - } - } - } -} -} // end namespace itk -#endif diff --git a/Modules/ImageStatistics/Testing/mitkMultiGaussianTest.cpp b/Modules/ImageStatistics/Testing/mitkMultiGaussianTest.cpp index 96440f5769..00f21a4086 100644 --- a/Modules/ImageStatistics/Testing/mitkMultiGaussianTest.cpp +++ b/Modules/ImageStatistics/Testing/mitkMultiGaussianTest.cpp @@ -1,305 +1,900 @@ /*=================================================================== 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 #include #include #include #include +#include #include +// Commandline:(for exmaple) mitkMultiGaussianTest C:/temp/output C:/temp/inputFile.xml +// +//For Example: input.xml +// +// > +// +// +// +// +// +//For Example: output.xml +// +// +// +// +// +// +// +// +// //bin\Release\ImageStatisticsTestDriver.exe mitkMultiGaussianTest C:/temp/TestImage/image 1 12 12 10 1 1 5 20 200 "2.5" "2.5" 3 +// +// bin\Release\ImageStatisticsTestDriver.exe mitkMultiGaussianTest C:/temp/HotSpotTestImage/ImageNeu30 5 30 30 15 1 5 15 20 200 1 1 2 int mitkMultiGaussianTest(int argc, char* argv[]) { - if ( argc !=14 ) + //if ( argc != 14 || argc != 3 ) + //{ + // std::cerr << " 14 arguments expected: [outputFilename, numberOfImages, imageSizeX, imageSizeY, imageSizeZ, numberOfGaussians, minWidthOfGaussian, maxWidthOfGaussian, minAltitudeOfGaussian, maxAltitudeOfGaussian, spacingX, spacingY, spacingZ ] \n OR \n 3 argument expected : [outputFilename, input.xml] \n " << std::endl; + // return EXIT_FAILURE; + //} + //// always start with this! + //MITK_TEST_BEGIN("mitkMultiGaussianTest"); + //MITK_TEST_CONDITION_REQUIRED(argc == 14, "Test called with 14 parameters"); + + if( argc == 14) + { + + + const unsigned int Dimension = 3; + typedef double PixelType; + typedef itk::DOMNode::Pointer DOMNodeType; + typedef itk::Image ImageType; + typedef itk::MultiGaussianImageSource< ImageType > MultiGaussianImageSource; + std::string outputFilename = argv[1], name; + int numberOfImages; + double centerX, centerY, centerZ, sigmaX, sigmaY, sigmaZ, altitude; + unsigned int numberOfGaussians, minWidthOfGaussian, maxWidthOfGaussian, minAltitudeOfGaussian, maxAltitudeOfGaussian; + itk::MultiGaussianImageSource< ImageType >::VectorType centerXVec, centerYVec, centerZVec, sigmaXVec, sigmaYVec, sigmaZVec, altitudeVec; + itk::MultiGaussianImageSource< ImageType >::ItkVectorType regionOfInterestMax, regionOfInterestMin; + itk::MultiGaussianImageSource< ImageType >::IndexType sphereMidpt, maxValueIndexInSphere, minValueIndexInSphere; + itk::MultiGaussianImageSource< ImageType >::ItkVectorType X, Y, Z, Sig, Alt; + MultiGaussianImageSource::Pointer gaussianGenerator; + itk::DOMNodeXMLWriter::Pointer xmlWriter; + itk::MultiGaussianImageSource< ImageType >::SpacingValueArrayType spacing; + DOMNodeType domTestCase, domTestImage, domGaussian, domStatistics, domROI; + ImageType::SizeValueType size[3]; + std::stringstream ss; + double radius = pow(itk::Math::one_over_pi * 0.75 , 1.0 / 3.0) * 10; + char * fileNamePointer; + + if ( ! (std::istringstream(argv[2]) >> numberOfImages) ) numberOfImages = 0; + if ( ! (std::istringstream(argv[3]) >> size[0]) ) size[0] = 0; + if ( ! (std::istringstream(argv[4]) >> size[1]) ) size[1] = 0; + if ( ! (std::istringstream(argv[5]) >> size[2]) ) size[2] = 0; + if ( ! (std::istringstream(argv[6]) >> numberOfGaussians) ) numberOfGaussians = 0; + if ( ! (std::istringstream(argv[7]) >> minWidthOfGaussian) ) minWidthOfGaussian = (size[0] + size[1] + size[2]) / 27; + if ( ! (std::istringstream(argv[8]) >> maxWidthOfGaussian) ) maxWidthOfGaussian = (size[0] + size[1] + size[2]) / 9; + if ( ! (std::istringstream(argv[9]) >> minAltitudeOfGaussian) ) minAltitudeOfGaussian = 5; + if ( ! (std::istringstream(argv[10]) >> maxAltitudeOfGaussian) ) maxAltitudeOfGaussian = 200; + if ( ! (std::istringstream(argv[11]) >> spacing[0]) ) spacing[0] = 1; + if ( ! (std::istringstream(argv[12]) >> spacing[1]) ) spacing[1] = 1; + if ( ! (std::istringstream(argv[13]) >> spacing[2]) ) spacing[2] = 1; + + // Set region of interest in pixels + regionOfInterestMax.SetElement(0, size[0] - static_cast((radius)/spacing[0] + 0.9999) - 1); + regionOfInterestMax.SetElement(1, size[1] - static_cast((radius)/spacing[1] + 0.9999) - 1); + regionOfInterestMax.SetElement(2, size[2] - static_cast((radius)/spacing[2] + 0.9999) - 1); + regionOfInterestMin.SetElement(0, 0 + static_cast((radius)/spacing[0]+ 0.9999)); + regionOfInterestMin.SetElement(1, 0 + static_cast((radius)/spacing[1]+ 0.9999)); + regionOfInterestMin.SetElement(2, 0 + static_cast((radius)/spacing[2]+ 0.9999)); + + srand (time(NULL)); + int numberAddGaussian = numberOfGaussians; + unsigned int k = 0; + unsigned int count = 0; + while(k != numberOfImages && count < 500) + // for(unsigned int k = 1; k <= numberOfImages; ++k) { - std::cerr << " 14 arguments expected: [outputFilename, numberOfImages, imageSizeX, imageSizeY, imageSizeZ, numberOfGaussians, minWidthOfGaussian, maxWidthOfGaussian, minAltitudeOfGaussian, maxAltitudeOfGaussian, spacingX, spacingY, spacingZ ] " << std::endl; - return EXIT_FAILURE; + ++count; + gaussianGenerator = MultiGaussianImageSource::New(); + gaussianGenerator->SetSize( size ); + gaussianGenerator->SetSpacing( spacing ); + gaussianGenerator->SetRadiusStepNumber(8); + gaussianGenerator->SetRadius(radius); + gaussianGenerator->SetNumberOfGausssians(numberOfGaussians); + gaussianGenerator->SetRegionOfInterest(regionOfInterestMin, regionOfInterestMax); + // DOM Node Writer + xmlWriter = itk::DOMNodeXMLWriter::New(); + domTestCase = itk::DOMNode::New(); + domTestCase->SetName("testcase"); + domTestImage = itk::DOMNode::New(); + domTestImage->SetName("testimage"); + ss.str(""); + ss << size[0]; + domTestImage->SetAttribute("image-rows", ss.str()); + ss.str(""); + ss << size[1]; + domTestImage->SetAttribute("image-columns", ss.str()); + ss.str(""); + ss << size[2]; + domTestImage->SetAttribute("image-slices", ss.str()); + ss.str(""); + ss << numberOfGaussians; + domTestImage->SetAttribute("numberOfGaussians", ss.str()); + ss.str(""); + ss << spacing[0]; + domTestImage->SetAttribute("spacingX", ss.str()); + ss.str(""); + ss << spacing[1]; + domTestImage->SetAttribute("spacingY", ss.str()); + ss.str(""); + ss << spacing[2]; + domTestImage->SetAttribute("spacingZ", ss.str()); + domTestCase->AddChildAtBegin(domTestImage); + + for( unsigned int i = 0; i < numberAddGaussian; ++i) + { + + domGaussian = itk::DOMNode::New() ; + domGaussian->SetName("gaussian"); + domTestImage->AddChildAtEnd(domGaussian); + // generate the midpoint and the daviation in mm + centerX = rand() % static_cast(size[0] * spacing[0]); + ss.str(""); + ss << centerX; + domGaussian->SetAttribute("centerIndexX", ss.str()); + + centerY = rand() % static_cast(size[1] * spacing[1]); + ss.str(""); + ss << centerY; + domGaussian->SetAttribute("centerIndexY", ss.str()); + + centerZ = rand() % static_cast(size[2] * spacing[2]); + ss.str(""); + ss << centerZ; + domGaussian->SetAttribute("centerIndexZ", ss.str()); + + sigmaX = minWidthOfGaussian + rand() % (maxWidthOfGaussian - minWidthOfGaussian); + ss.str(""); + ss << sigmaX; + domGaussian->SetAttribute("deviationX", ss.str()); + + sigmaY = minWidthOfGaussian + rand() % (maxWidthOfGaussian - minWidthOfGaussian); + ss.str(""); + ss << sigmaY; + domGaussian->SetAttribute("deviationY", ss.str()); + + sigmaZ = minWidthOfGaussian + rand() % (maxWidthOfGaussian - minWidthOfGaussian); + ss.str(""); + ss << sigmaZ; + domGaussian->SetAttribute("deviationZ", ss.str()); + + altitude = minAltitudeOfGaussian + rand() % (maxAltitudeOfGaussian - minAltitudeOfGaussian); + ss.str(""); + ss << altitude; + domGaussian->SetAttribute("altitude", ss.str()); + + 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); + + } + + gaussianGenerator->AddGaussian(centerXVec, centerYVec, centerZVec, sigmaXVec, sigmaYVec, sigmaZVec, altitudeVec); + centerXVec.clear(); + centerYVec.clear(); + centerZVec.clear(); + sigmaXVec.clear(); + sigmaYVec.clear(); + sigmaZVec.clear(); + altitudeVec.clear(); + try { + gaussianGenerator->Update(); + gaussianGenerator->CalculateMidpointAndMeanValue(); + } catch (std::exception& e) + { + std::cout << "Error: " << e.what() << std::endl; + } + + //region of interest + domROI = itk::DOMNode::New(); + domROI->SetName("roi"); + domTestCase->AddChildAtEnd(domROI); + ss.str(""); + ss << radius; + domROI->SetAttribute("hotspotRadiusInMM", ss.str()); + ss.str(""); + ss << regionOfInterestMax[0]; + domROI->SetAttribute("maximumSizeX", ss.str()); + ss.str(""); + ss << regionOfInterestMin[0]; + domROI->SetAttribute("minimumSizeX", ss.str()); + ss.str(""); + ss << regionOfInterestMax[1]; + domROI->SetAttribute("maximumSizeY", ss.str()); + ss.str(""); + ss << regionOfInterestMin[1]; + domROI->SetAttribute("minimumSizeY", ss.str()); + ss.str(""); + ss << regionOfInterestMax[2]; + domROI->SetAttribute("maximumSizeZ", ss.str()); + ss.str(""); + ss << regionOfInterestMin[2]; + domROI->SetAttribute("minimumSizeZ", ss.str()); + + //peak and peak coordinate + domStatistics = itk::DOMNode::New(); + domStatistics->SetName("statistic"); + domTestCase->AddChildAtEnd(domStatistics); + sphereMidpt = gaussianGenerator->GetSphereMidpoint(); + ss.str(""); + ss << sphereMidpt[0]; + domStatistics->SetAttribute("hotspotIndexX", ss.str()); + ss.str(""); + ss << sphereMidpt[1]; + domStatistics->SetAttribute("hotspotIndexY", ss.str()); + ss.str(""); + ss << sphereMidpt[2]; + domStatistics->SetAttribute("hotspotIndexZ", ss.str()); + ss.str(""); + ss << gaussianGenerator->GetMaxMeanValue(); + domStatistics->SetAttribute("peak", ss.str()); + + //optimize the mean value in the sphere + // gaussianGenerator->OptimizeMeanValue(); + ss.str(""); + ss << gaussianGenerator->GetMaxMeanValue(); + domStatistics->SetAttribute("mean", ss.str()); + + + //maximum and maximum coordinate + gaussianGenerator->CalculateMaxAndMinInSphere(); + maxValueIndexInSphere = gaussianGenerator->GetMaxValueIndexInSphere(); + ss.str(""); + ss << maxValueIndexInSphere[0]; + domStatistics->SetAttribute("maximumIndexX", ss.str()); + ss.str(""); + ss << maxValueIndexInSphere[1]; + domStatistics->SetAttribute("maximumIndexY", ss.str()); + ss.str(""); + ss << maxValueIndexInSphere[2]; + domStatistics->SetAttribute("maximumIndexZ", ss.str()); + ss.str(""); + ss << gaussianGenerator->GetMaxValueInSphere(); + domStatistics->SetAttribute("maximum", ss.str()); + + //minimum and minimum coordinate + minValueIndexInSphere = gaussianGenerator->GetMinValueIndexInSphere(); + ss.str(""); + ss << minValueIndexInSphere[0]; + domStatistics->SetAttribute("minimumIndexX", ss.str()); + ss.str(""); + ss << minValueIndexInSphere[1]; + domStatistics->SetAttribute("minimumIndexY", ss.str()); + ss.str(""); + ss << minValueIndexInSphere[2]; + domStatistics->SetAttribute("minimumIndexZ", ss.str()); + ss.str(""); + ss << gaussianGenerator->GetMinValueInSphere(); + domStatistics->SetAttribute("minimum", ss.str()); + + + // test if the midpoint of the sphere is not the same as the maximum's index and saves only such a case + if( sphereMidpt[0]!= maxValueIndexInSphere[0] && sphereMidpt[1]!= maxValueIndexInSphere[1] && sphereMidpt[2]!= maxValueIndexInSphere[2]) + { + k++; + // .xml (Data) + ss.str(""); + if(k < 10){ + ss << outputFilename <<"00"<< k <<".xml"; + }else if(k < 100){ + ss << outputFilename <<"0"<< k <<".xml"; + }else{ ss << outputFilename << k <<".xml";} + name = ss.str(); + fileNamePointer = (char*) name.c_str(); + xmlWriter->SetFileName( fileNamePointer); + xmlWriter->SetInput( domTestCase ); + xmlWriter->Update(); + ImageType::Pointer gaussianImage = gaussianGenerator->GetOutput(); + + //.nrrd (Image) + typedef itk::ImageFileWriter< ImageType > WriterType; + WriterType::Pointer writer = WriterType::New(); + ss.str(""); + if(k < 10) + { + ss << outputFilename <<"00"<< k <<".nrrd"; + }else if(k < 100) + { + ss << outputFilename <<"0"<< k <<".nrrd"; + }else + { + ss << outputFilename << k <<".nrrd"; + } + name = ss.str(); + fileNamePointer = (char*) name.c_str(); + writer->SetFileName( fileNamePointer); + writer->SetInput( gaussianImage ); + writer->Update(); + } } - // always start with this! - MITK_TEST_BEGIN("mitkMultiGaussianTest"); - MITK_TEST_CONDITION_REQUIRED(argc == 14, "Test called with 14 parameters"); - - const unsigned int Dimension = 3; - typedef double PixelType; - typedef itk::DOMNode::Pointer DOMNodeType; - typedef itk::Image ImageType; - typedef itk::MultiGaussianImageSource< ImageType > MultiGaussianImageSource; - std::string outputFilename = argv[1], name; - int numberOfImages; - double centerX, centerY, centerZ, sigmaX, sigmaY, sigmaZ, altitude; - unsigned int numberOfGaussians, minWidthOfGaussian, maxWidthOfGaussian, minAltitudeOfGaussian, maxAltitudeOfGaussian; - itk::MultiGaussianImageSource< ImageType >::VectorType centerXVec, centerYVec, centerZVec, sigmaXVec, sigmaYVec, sigmaZVec, altitudeVec; - itk::MultiGaussianImageSource< ImageType >::ItkVectorType regionOfInterestMax, regionOfInterestMin; - itk::MultiGaussianImageSource< ImageType >::IndexType sphereMidpt, maxValueIndexInSphere, minValueIndexInSphere; - itk::MultiGaussianImageSource< ImageType >::ItkVectorType X, Y, Z, Sig, Alt; - MultiGaussianImageSource::Pointer gaussianGenerator; - itk::DOMNodeXMLWriter::Pointer xmlWriter; - itk::MultiGaussianImageSource< ImageType >::SpacingValueArrayType spacing; - DOMNodeType domTestCase, domTestImage, domGaussian, domStatistics, domROI; - ImageType::SizeValueType size[3]; - std::stringstream ss; - double radius = pow(itk::Math::one_over_pi * 0.75 , 1.0 / 3.0) * 10; - char * fileNamePointer; - - if ( ! (std::istringstream(argv[2]) >> numberOfImages) ) numberOfImages = 0; - if ( ! (std::istringstream(argv[3]) >> size[0]) ) size[0] = 0; - if ( ! (std::istringstream(argv[4]) >> size[1]) ) size[1] = 0; - if ( ! (std::istringstream(argv[5]) >> size[2]) ) size[2] = 0; - if ( ! (std::istringstream(argv[6]) >> numberOfGaussians) ) numberOfGaussians = 0; - if ( ! (std::istringstream(argv[7]) >> minWidthOfGaussian) ) minWidthOfGaussian = (size[0] + size[1] + size[2]) / 27; - if ( ! (std::istringstream(argv[8]) >> maxWidthOfGaussian) ) maxWidthOfGaussian = (size[0] + size[1] + size[2]) / 9; - if ( ! (std::istringstream(argv[9]) >> minAltitudeOfGaussian) ) minAltitudeOfGaussian = 5; - if ( ! (std::istringstream(argv[10]) >> maxAltitudeOfGaussian) ) maxAltitudeOfGaussian = 200; - if ( ! (std::istringstream(argv[11]) >> spacing[0]) ) spacing[0] = 1; - if ( ! (std::istringstream(argv[12]) >> spacing[1]) ) spacing[1] = 1; - if ( ! (std::istringstream(argv[13]) >> spacing[2]) ) spacing[2] = 1; - - // Set region of interest in pixels - regionOfInterestMax.SetElement(0, size[0] - static_cast((radius)/spacing[0] + 0.9999) - 1); - regionOfInterestMax.SetElement(1, size[1] - static_cast((radius)/spacing[1] + 0.9999) - 1); - regionOfInterestMax.SetElement(2, size[2] - static_cast((radius)/spacing[2] + 0.9999) - 1); - regionOfInterestMin.SetElement(0, 0 + static_cast((radius)/spacing[0]+ 0.9999)); - regionOfInterestMin.SetElement(1, 0 + static_cast((radius)/spacing[1]+ 0.9999)); - regionOfInterestMin.SetElement(2, 0 + static_cast((radius)/spacing[2]+ 0.9999)); - - srand (time(NULL)); - int numberAddGaussian = numberOfGaussians; - for(unsigned int k = 1; k <= numberOfImages; ++k) + } + + + //-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- + // In the inputFile.xml we find the characteristics of the Gaussian and the ROI's. Hier we can have more then one ROI -> we find the hot spot for each of the ROI's. + else { - gaussianGenerator = MultiGaussianImageSource::New(); - gaussianGenerator->SetSize( size ); - gaussianGenerator->SetSpacing( spacing ); - gaussianGenerator->SetRadiusStepNumber(8); - gaussianGenerator->SetRadius(radius); - gaussianGenerator->SetNumberOfGausssians(numberOfGaussians); - gaussianGenerator->SetRegionOfInterest(regionOfInterestMin, regionOfInterestMax); - // DOM Node Writer + + const unsigned int Dimension = 3; + typedef double PixelType; + typedef itk::DOMNode::Pointer DOMNodeType; + typedef itk::Image ImageType; + typedef itk::MultiGaussianImageSource< ImageType > MultiGaussianImageSource; + std::string outputFilename = argv[1], name; + int numberOfImages; + double centerX, centerY, centerZ, sigmaX, sigmaY, sigmaZ, altitude, hotSpotRadiusInMM; + unsigned int numberOfGaussians, minWidthOfGaussian, maxWidthOfGaussian, minAltitudeOfGaussian, maxAltitudeOfGaussian, numberOfLabels; + itk::MultiGaussianImageSource< ImageType >::VectorType centerXVec, centerYVec, centerZVec, sigmaXVec, sigmaYVec, sigmaZVec, altitudeVec, ROImaxSizeX, ROIminSizeX, ROImaxSizeY, ROIminSizeY, ROImaxSizeZ, ROIminSizeZ, label; + itk::MultiGaussianImageSource< ImageType >::ItkVectorType regionOfInterestMax, regionOfInterestMin; + itk::MultiGaussianImageSource< ImageType >::IndexType sphereMidpt, maxValueIndexInSphere, minValueIndexInSphere; + itk::MultiGaussianImageSource< ImageType >::ItkVectorType X, Y, Z, Sig, Alt; + MultiGaussianImageSource::Pointer gaussianGenerator; + itk::DOMNodeXMLWriter::Pointer xmlWriter; + itk::MultiGaussianImageSource< ImageType >::SpacingValueArrayType spacing; + DOMNodeType domTestCase, domTestImage, domGaussian, domSegmentation, domStatistics, domROI; + ImageType::SizeValueType size[3]; + std::stringstream ss; + double radius = pow(itk::Math::one_over_pi * 0.75 , 1.0 / 3.0) * 10; + char * fileNamePointer; + std::string attributeValue; + double value; + bool entireHotSpotInROI; + int maxSize, minSize; + + std::string filename = argv[2]; + itk::DOMNodeXMLReader::Pointer xmlReader = itk::DOMNodeXMLReader::New(); + xmlReader->SetFileName( filename ); + 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]; + + attributeValue = testimage->GetAttribute("image-rows"); + std::stringstream(attributeValue) >> size[0]; + attributeValue = testimage->GetAttribute("image-columns"); + std::stringstream(attributeValue) >> size[1]; + attributeValue = testimage->GetAttribute("image-slices"); + std::stringstream(attributeValue) >> size[2]; + + attributeValue = testimage->GetAttribute( "numberOfGaussians" ); + std::stringstream(attributeValue) >> numberOfGaussians; + + attributeValue = testimage->GetAttribute( "spacingX" ); + std::stringstream(attributeValue) >> spacing[0]; + attributeValue = testimage->GetAttribute( "spacingY" ); + std::stringstream(attributeValue) >> spacing[1]; + attributeValue = testimage->GetAttribute( "spacingZ" ); + std::stringstream(attributeValue) >> spacing[2]; + attributeValue = testimage->GetAttribute( "entireHotSpotInROI" ); + std::stringstream(attributeValue) >> entireHotSpotInROI; + + std::cout << "Read size parameters (x,y,z): " << size[0] << "," << size[1] << "," << size[2] << "\n" << std::endl; + std::cout << "Read spacing parameters (x,y,z): " << spacing[0] << "," << spacing[1] << "," << spacing[2]<< "\n" << std::endl; + + NodeList gaussians; + testimage->GetChildren("gaussian", gaussians); + itk::DOMNode* gaussian; + for(int i = 0; i < numberOfGaussians ; ++i) + { + gaussian = gaussians[i]; + + attributeValue = gaussian->GetAttribute( "centerIndexX" ); + std::stringstream(attributeValue) >> value; + centerXVec.push_back(value); + + attributeValue = gaussian->GetAttribute( "centerIndexY" ); + std::stringstream(attributeValue) >> value; + centerYVec.push_back(value); + + attributeValue = gaussian->GetAttribute( "centerIndexZ" ); + std::stringstream(attributeValue) >> value; + centerZVec.push_back(value); + + std::cout << "Read center of Gaussian (x,y,z): " << centerXVec[i] << "," << centerYVec[i] << "," << centerZVec[i] << "\n" << std::endl; + + attributeValue = gaussian->GetAttribute( "deviationX" ); + std::stringstream(attributeValue) >> value; + sigmaXVec.push_back(value); + + attributeValue = gaussian->GetAttribute( "deviationY" ); + std::stringstream(attributeValue) >> value; + sigmaYVec.push_back(value); + + attributeValue = gaussian->GetAttribute( "deviationZ" ); + std::stringstream(attributeValue) >> value; + sigmaZVec.push_back(value); + + std::cout << "Read deviation of Gaussian (x,y,z): " << sigmaXVec[i] << "," << sigmaYVec[i] << "," << sigmaZVec[i] << "\n" << std::endl; + + attributeValue = gaussian->GetAttribute( "altitude" ); + std::stringstream(attributeValue) >> value; + altitudeVec.push_back(value); + std::cout << "Read altitude: " << altitudeVec[i] << "\n" << std::endl; + } + + + + // read ROI's parameter + NodeList segmentations; + domRoot->GetChildren("segmentation", segmentations); + MITK_TEST_CONDITION_REQUIRED( segmentations.size() == 1, "One ROI image defined" ) + itk::DOMNode* segmentation = segmentations[0]; + + attributeValue = segmentation->GetAttribute("numberOfLabels"); + std::stringstream(attributeValue) >> numberOfLabels; + attributeValue = segmentation->GetAttribute("hotspotRadiusInMM"); + std::stringstream(attributeValue) >> hotSpotRadiusInMM; + + + std::cout << "Read number of labels: " << numberOfLabels << "\n" << std::endl; + std::cout << "Read radius in mm : " << hotSpotRadiusInMM << "\n" << std::endl; + + NodeList rois; + segmentation->GetChildren("roi", rois); + itk::DOMNode* roi; + for(int i = 0; i < numberOfLabels ; ++i) + { + roi = rois[i]; + + attributeValue = roi->GetAttribute( "label" ); + std::stringstream(attributeValue) >> value; + label.push_back(value); + + attributeValue = roi->GetAttribute( "maximumSizeX" ); + std::stringstream(attributeValue) >> value; + ROImaxSizeX.push_back(value); + + attributeValue = roi->GetAttribute( "minimumSizeX" ); + std::stringstream(attributeValue) >> value; + ROIminSizeX.push_back(value); + + attributeValue = roi->GetAttribute( "maximumSizeY" ); + std::stringstream(attributeValue) >> value; + ROImaxSizeY.push_back(value); + + attributeValue = roi->GetAttribute( "minimumSizeY" ); + std::stringstream(attributeValue) >> value; + ROIminSizeY.push_back(value); + + attributeValue = roi->GetAttribute( "maximumSizeZ" ); + std::stringstream(attributeValue) >> value; + ROImaxSizeZ.push_back(value); + + attributeValue = roi->GetAttribute( "minimumSizeZ" ); + std::stringstream(attributeValue) >> value; + ROIminSizeZ.push_back(value); + + std::cout << "Read ROI with label nummber: " << label[i] << "with min and max values in the x-, y-, z-Achse: [" << ROIminSizeX[i] << ROImaxSizeX[i] <<"], [" << ROIminSizeY[i] << ROImaxSizeY[i] <<"], [" << ROIminSizeZ[i] << ROImaxSizeZ[i] <<"]\n" << std::endl; + } + + //write test image parameter + xmlWriter = itk::DOMNodeXMLWriter::New(); domTestCase = itk::DOMNode::New(); domTestCase->SetName("testcase"); domTestImage = itk::DOMNode::New(); domTestImage->SetName("testimage"); ss.str(""); ss << size[0]; domTestImage->SetAttribute("image-rows", ss.str()); ss.str(""); ss << size[1]; domTestImage->SetAttribute("image-columns", ss.str()); ss.str(""); ss << size[2]; domTestImage->SetAttribute("image-slices", ss.str()); ss.str(""); ss << numberOfGaussians; domTestImage->SetAttribute("numberOfGaussians", ss.str()); ss.str(""); ss << spacing[0]; domTestImage->SetAttribute("spacingX", ss.str()); ss.str(""); ss << spacing[1]; domTestImage->SetAttribute("spacingY", ss.str()); ss.str(""); ss << spacing[2]; domTestImage->SetAttribute("spacingZ", ss.str()); + ss.str(""); + ss << entireHotSpotInROI; + domTestImage->SetAttribute("entireHotSpotInROI", ss.str()); + domTestCase->AddChildAtBegin(domTestImage); + int numberAddGaussian = numberOfGaussians; for( unsigned int i = 0; i < numberAddGaussian; ++i) { domGaussian = itk::DOMNode::New() ; domGaussian->SetName("gaussian"); domTestImage->AddChildAtEnd(domGaussian); // generate the midpoint and the daviation in mm - centerX = rand() % static_cast(size[0] * spacing[0]); + centerX = centerXVec[i]; ss.str(""); ss << centerX; domGaussian->SetAttribute("centerIndexX", ss.str()); - centerY = rand() % static_cast(size[1] * spacing[1]); + centerY = centerYVec[i]; ss.str(""); ss << centerY; domGaussian->SetAttribute("centerIndexY", ss.str()); - centerZ = rand() % static_cast(size[2] * spacing[2]); + centerZ = centerZVec[i]; ss.str(""); ss << centerZ; domGaussian->SetAttribute("centerIndexZ", ss.str()); - sigmaX = minWidthOfGaussian + rand() % (maxWidthOfGaussian - minWidthOfGaussian); + sigmaX = sigmaXVec[i]; ss.str(""); ss << sigmaX; domGaussian->SetAttribute("deviationX", ss.str()); - sigmaY = minWidthOfGaussian + rand() % (maxWidthOfGaussian - minWidthOfGaussian); + sigmaY = sigmaYVec[i]; ss.str(""); ss << sigmaY; domGaussian->SetAttribute("deviationY", ss.str()); - sigmaZ = minWidthOfGaussian + rand() % (maxWidthOfGaussian - minWidthOfGaussian); + sigmaZ = sigmaZVec[i]; ss.str(""); ss << sigmaZ; domGaussian->SetAttribute("deviationZ", ss.str()); - altitude = minAltitudeOfGaussian + rand() % (maxAltitudeOfGaussian - minAltitudeOfGaussian); + altitude = altitudeVec[i]; ss.str(""); ss << altitude; domGaussian->SetAttribute("altitude", ss.str()); - 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); - } + radius = hotSpotRadiusInMM; + + gaussianGenerator = MultiGaussianImageSource::New(); + gaussianGenerator->SetSize( size ); + gaussianGenerator->SetSpacing( spacing ); + gaussianGenerator->SetRadiusStepNumber(8); + gaussianGenerator->SetRadius(radius); + gaussianGenerator->SetNumberOfGausssians(numberOfGaussians); gaussianGenerator->AddGaussian(centerXVec, centerYVec, centerZVec, sigmaXVec, sigmaYVec, sigmaZVec, altitudeVec); - centerXVec.clear(); - centerYVec.clear(); - centerZVec.clear(); - sigmaXVec.clear(); - sigmaYVec.clear(); - sigmaZVec.clear(); - altitudeVec.clear(); - try { - gaussianGenerator->Update(); - gaussianGenerator->CalculateMidpointAndMeanValue(); - } catch (std::exception& e) - { - std::cout << "Error: " << e.what() << std::endl; - } - //region of interest - domROI = itk::DOMNode::New(); - domROI->SetName("roi"); - domTestCase->AddChildAtEnd(domROI); - ss.str(""); - ss << regionOfInterestMax[0]; - domROI->SetAttribute("maximumX", ss.str()); - ss.str(""); - ss << regionOfInterestMin[0]; - domROI->SetAttribute("minimumX", ss.str()); - ss.str(""); - ss << regionOfInterestMax[1]; - domROI->SetAttribute("maximumY", ss.str()); - ss.str(""); - ss << regionOfInterestMin[1]; - domROI->SetAttribute("minimumY", ss.str()); - ss.str(""); - ss << regionOfInterestMax[2]; - domROI->SetAttribute("maximumZ", ss.str()); - ss.str(""); - ss << regionOfInterestMin[2]; - domROI->SetAttribute("minimumZ", ss.str()); - - //peak and peak coordinate - domStatistics = itk::DOMNode::New(); - domStatistics->SetName("statistic"); - domTestCase->AddChildAtEnd(domStatistics); - sphereMidpt = gaussianGenerator->GetSphereMidpoint(); - ss.str(""); - ss << sphereMidpt[0]; - domStatistics->SetAttribute("peakIndexX", ss.str()); - ss.str(""); - ss << sphereMidpt[1]; - domStatistics->SetAttribute("peakIndexY", ss.str()); + domSegmentation = itk::DOMNode::New(); + domSegmentation->SetName("segmentation"); + domTestCase->AddChildAtEnd(domSegmentation); + + ss.str(""); - ss << sphereMidpt[2]; - domStatistics->SetAttribute("peakIndexZ", ss.str()); + ss << numberOfLabels; + domSegmentation->SetAttribute("numberOfLabels", ss.str()); ss.str(""); - ss << gaussianGenerator->GetMaxMeanValue(); - domStatistics->SetAttribute("peak", ss.str()); + ss << hotSpotRadiusInMM; + domSegmentation->SetAttribute("hotspotRadiusInMM", ss.str()); - //optimize the mean value in the sphere - gaussianGenerator->OptimizeMeanValue(); - ss.str(""); - ss << gaussianGenerator->GetMaxMeanValue(); - domStatistics->SetAttribute("peakOptimized", ss.str()); + // Set region of interest in index values + if(entireHotSpotInROI) + { + // set the region of interest for each label i + for (unsigned int i = 0; i < numberOfLabels; ++i) + { + maxSize = ROImaxSizeX[i] - static_cast((radius)/spacing[0] + 0.9999) - 1; + minSize = ROIminSizeX[i] + static_cast((radius)/spacing[0]+ 0.9999); + if(minSize > maxSize) + { + std::cout << "The sphere is larger then the region of interest! Set the roi to be the whole image " << std::endl; + maxSize = size[0]; + minSize = 0.0; + } + // the maximum in the x-Axis + regionOfInterestMax.SetElement( 0, (maxSize < size[0]) ? maxSize : size[0] ); + // the minimum in the x-Axis + regionOfInterestMin.SetElement(0, (minSize > 0.0) ? minSize : 0.0 ); + + + maxSize = ROImaxSizeY[i] - static_cast((radius)/spacing[1] + 0.9999) - 1; + minSize = ROIminSizeY[i] + static_cast((radius)/spacing[1]+ 0.9999); + if(minSize > maxSize) + { + std::cout << "The sphere is larger then the region of interest! Set the roi to be the whole image " << std::endl; + maxSize = size[1]; + minSize = 0.0; + } + // the maximum in the y-Axis + regionOfInterestMax.SetElement(1, (maxSize < size[1]) ? maxSize : size[1] ); + // the minimum in the y-Axis + regionOfInterestMin.SetElement(1, (minSize > 0.0) ? minSize : 0.0 ); + + maxSize = ROImaxSizeZ[i] - static_cast((radius)/spacing[2] + 0.9999) - 1; + minSize = ROIminSizeZ[i] + static_cast((radius)/spacing[2]+ 0.9999); + if(minSize > maxSize) + { + std::cout << "The sphere is larger then the region of interest! Set the roi to be the whole image " << std::endl; + maxSize = size[2]; + minSize = 0.0; + } + // the maximum in the z-Axis + regionOfInterestMax.SetElement(2, (maxSize < size[2]) ? maxSize : size[2] ); + // the minimum in the z-Axis + regionOfInterestMin.SetElement(2, (minSize > 0.0) ? minSize : 0.0 ); + + + + + gaussianGenerator->SetRegionOfInterest(regionOfInterestMin, regionOfInterestMax); + // TODO !! Kann ich die Update Methode ausserhalb der Schleife machen, nachdem ich Calculate Midpoint and Mean Value aufgerufen habe? Wird es korrekt berechnt? Letztendlich sind diese beide von einander unabhaengig. Unten auch! + gaussianGenerator->Update(); + // gaussianGenerator->CalculateMidpointAndMeanValue(); + // gaussianGenerator->CalculateMidpoint(); + + //TODO + gaussianGenerator->GenerateCuboidSegmentationInSphere(); + + + + //region of interest + domROI = itk::DOMNode::New(); + domROI->SetName("roi"); + domSegmentation->AddChildAtEnd(domROI); + + ss.str(""); + ss << label[i]; + domROI->SetAttribute("label", ss.str()); + ss.str(""); + ss << ROImaxSizeX[i]; //regionOfInterestMax[0]; + domROI->SetAttribute("maximumSizeX", ss.str()); + ss.str(""); + ss << ROIminSizeX[i]; //regionOfInterestMin[0]; + domROI->SetAttribute("minimumSizeX", ss.str()); + ss.str(""); + ss << ROImaxSizeY[i]; //regionOfInterestMax[1]; + domROI->SetAttribute("maximumSizeY", ss.str()); + ss.str(""); + ss << ROIminSizeY[i]; //regionOfInterestMin[1]; + domROI->SetAttribute("minimumSizeY", ss.str()); + ss.str(""); + ss << ROImaxSizeZ[i]; //regionOfInterestMax[2]; + domROI->SetAttribute("maximumSizeZ", ss.str()); + ss.str(""); + ss << ROIminSizeZ[i]; //regionOfInterestMin[2]; + domROI->SetAttribute("minimumSizeZ", ss.str()); + + //peak and peak coordinate + domStatistics = itk::DOMNode::New(); + domStatistics->SetName("statistic"); + domTestCase->AddChildAtEnd(domStatistics); + sphereMidpt = gaussianGenerator->GetSphereMidpoint(); + ss.str(""); + ss << sphereMidpt[0]; + domStatistics->SetAttribute("hotspotIndexX", ss.str()); + ss.str(""); + ss << sphereMidpt[1]; + domStatistics->SetAttribute("hotspotIndexY", ss.str()); + ss.str(""); + ss << sphereMidpt[2]; + domStatistics->SetAttribute("hotspotIndexZ", ss.str()); + ss.str(""); + ss << gaussianGenerator->GetMaxMeanValue(); + domStatistics->SetAttribute("peak", ss.str()); + + //optimize the mean value in the sphere + // gaussianGenerator->OptimizeMeanValue(); + ss.str(""); + ss << gaussianGenerator->GetMaxMeanValue(); + domStatistics->SetAttribute("mean", ss.str()); + + + //maximum and maximum coordinate + gaussianGenerator->CalculateMaxAndMinInSphere(); + maxValueIndexInSphere = gaussianGenerator->GetMaxValueIndexInSphere(); + ss.str(""); + ss << maxValueIndexInSphere[0]; + domStatistics->SetAttribute("maximumIndexX", ss.str()); + ss.str(""); + ss << maxValueIndexInSphere[1]; + domStatistics->SetAttribute("maximumIndexY", ss.str()); + ss.str(""); + ss << maxValueIndexInSphere[2]; + domStatistics->SetAttribute("maximumIndexZ", ss.str()); + ss.str(""); + ss << gaussianGenerator->GetMaxValueInSphere(); + domStatistics->SetAttribute("maximum", ss.str()); + + //minimum and minimum coordinate + minValueIndexInSphere = gaussianGenerator->GetMinValueIndexInSphere(); + ss.str(""); + ss << minValueIndexInSphere[0]; + domStatistics->SetAttribute("minimumIndexX", ss.str()); + ss.str(""); + ss << minValueIndexInSphere[1]; + domStatistics->SetAttribute("minimumIndexY", ss.str()); + ss.str(""); + ss << minValueIndexInSphere[2]; + domStatistics->SetAttribute("minimumIndexZ", ss.str()); + ss.str(""); + ss << gaussianGenerator->GetMinValueInSphere(); + domStatistics->SetAttribute("minimum", ss.str()); + + } + } + else + { + // set the region of interest for each label i + for (unsigned int i = 0; i < numberOfLabels; ++i) + { + for (unsigned int k = 0 ; k < 3; ++k) + { + + maxSize = ROImaxSizeX[i]; + minSize = ROIminSizeX[i]; + if(minSize > maxSize) + { + std::cout << "The sphere is larger then the region of interest! Set the roi to be the whole image " << std::endl; + maxSize = size[k]; + minSize = 0.0; + } + // the maximum in the k-Axis + regionOfInterestMax.SetElement(k, (maxSize < size[k]) ? maxSize : size[k] ); + // the minimum in the k-Axis + regionOfInterestMin.SetElement(k, (minSize > 0.0) ? minSize : 0.0 ); + } + + + gaussianGenerator->SetRegionOfInterest(regionOfInterestMin, regionOfInterestMax); + gaussianGenerator->Update(); + gaussianGenerator->CalculateMidpointAndMeanValue(); + + + //region of interest + domROI = itk::DOMNode::New(); + domROI->SetName("roi"); + domSegmentation->AddChildAtEnd(domROI); + + ss.str(""); + ss << label[i]; + domROI->SetAttribute("label", ss.str()); + ss.str(""); + ss << regionOfInterestMax[0]; + domROI->SetAttribute("maximumSizeX", ss.str()); + ss.str(""); + ss << regionOfInterestMin[0]; + domROI->SetAttribute("minimumSizeX", ss.str()); + ss.str(""); + ss << regionOfInterestMax[1]; + domROI->SetAttribute("maximumSizeY", ss.str()); + ss.str(""); + ss << regionOfInterestMin[1]; + domROI->SetAttribute("minimumSizeY", ss.str()); + ss.str(""); + ss << regionOfInterestMax[2]; + domROI->SetAttribute("maximumSizeZ", ss.str()); + ss.str(""); + ss << regionOfInterestMin[2]; + domROI->SetAttribute("minimumSizeZ", ss.str()); + + //peak and peak coordinate + domStatistics = itk::DOMNode::New(); + domStatistics->SetName("statistic"); + domTestCase->AddChildAtEnd(domStatistics); + sphereMidpt = gaussianGenerator->GetSphereMidpoint(); + ss.str(""); + ss << sphereMidpt[0]; + domStatistics->SetAttribute("hotspotIndexX", ss.str()); + ss.str(""); + ss << sphereMidpt[1]; + domStatistics->SetAttribute("hotspotIndexY", ss.str()); + ss.str(""); + ss << sphereMidpt[2]; + domStatistics->SetAttribute("hotspotIndexZ", ss.str()); + ss.str(""); + ss << gaussianGenerator->GetMaxMeanValue(); + domStatistics->SetAttribute("peak", ss.str()); + + //optimize the mean value in the sphere + // gaussianGenerator->OptimizeMeanValue(); + ss.str(""); + ss << gaussianGenerator->GetMaxMeanValue(); + domStatistics->SetAttribute("mean", ss.str()); + + + //maximum and maximum coordinate + gaussianGenerator->CalculateMaxAndMinInSphere(); + maxValueIndexInSphere = gaussianGenerator->GetMaxValueIndexInSphere(); + ss.str(""); + ss << maxValueIndexInSphere[0]; + domStatistics->SetAttribute("maximumIndexX", ss.str()); + ss.str(""); + ss << maxValueIndexInSphere[1]; + domStatistics->SetAttribute("maximumIndexY", ss.str()); + ss.str(""); + ss << maxValueIndexInSphere[2]; + domStatistics->SetAttribute("maximumIndexZ", ss.str()); + ss.str(""); + ss << gaussianGenerator->GetMaxValueInSphere(); + domStatistics->SetAttribute("maximum", ss.str()); + + //minimum and minimum coordinate + minValueIndexInSphere = gaussianGenerator->GetMinValueIndexInSphere(); + ss.str(""); + ss << minValueIndexInSphere[0]; + domStatistics->SetAttribute("minimumIndexX", ss.str()); + ss.str(""); + ss << minValueIndexInSphere[1]; + domStatistics->SetAttribute("minimumIndexY", ss.str()); + ss.str(""); + ss << minValueIndexInSphere[2]; + domStatistics->SetAttribute("minimumIndexZ", ss.str()); + ss.str(""); + ss << gaussianGenerator->GetMinValueInSphere(); + domStatistics->SetAttribute("minimum", ss.str()); + + } + } - //maximum and maximum coordinate - gaussianGenerator->CalculateMaxAndMinInSphere(); - maxValueIndexInSphere = gaussianGenerator->GetMaxValueIndexInSphere(); - ss.str(""); - ss << maxValueIndexInSphere[0]; - domStatistics->SetAttribute("maximumIndexX", ss.str()); - ss.str(""); - ss << maxValueIndexInSphere[1]; - domStatistics->SetAttribute("maximumIndexY", ss.str()); - ss.str(""); - ss << maxValueIndexInSphere[2]; - domStatistics->SetAttribute("maximumIndexZ", ss.str()); - ss.str(""); - ss << gaussianGenerator->GetMaxValueInSphere(); - domStatistics->SetAttribute("maximum", ss.str()); - //minimum and minimum coordinate - minValueIndexInSphere = gaussianGenerator->GetMinValueIndexInSphere(); - ss.str(""); - ss << minValueIndexInSphere[0]; - domStatistics->SetAttribute("minimumIndexX", ss.str()); - ss.str(""); - ss << minValueIndexInSphere[1]; - domStatistics->SetAttribute("minimumIndexY", ss.str()); - ss.str(""); - ss << minValueIndexInSphere[2]; - domStatistics->SetAttribute("minimumIndexZ", ss.str()); - ss.str(""); - ss << gaussianGenerator->GetMinValueInSphere(); - domStatistics->SetAttribute("minimum", ss.str()); // .xml (Data) ss.str(""); - if(k < 10){ - ss << outputFilename <<"00"<< k <<".xml"; - }else if(k < 100){ - ss << outputFilename <<"0"<< k <<".xml"; - }else{ ss << outputFilename << k <<".xml";} + ss << outputFilename << ".xml"; name = ss.str(); fileNamePointer = (char*) name.c_str(); xmlWriter->SetFileName( fileNamePointer); xmlWriter->SetInput( domTestCase ); xmlWriter->Update(); ImageType::Pointer gaussianImage = gaussianGenerator->GetOutput(); //.nrrd (Image) typedef itk::ImageFileWriter< ImageType > WriterType; WriterType::Pointer writer = WriterType::New(); ss.str(""); - if(k < 10) - { - ss << outputFilename <<"00"<< k <<".nrrd"; - }else if(k < 100) - { - ss << outputFilename <<"0"<< k <<".nrrd"; - }else - { - ss << outputFilename << k <<".nrrd"; - } + ss << outputFilename << ".nrrd"; name = ss.str(); fileNamePointer = (char*) name.c_str(); writer->SetFileName( fileNamePointer); writer->SetInput( gaussianImage ); writer->Update(); + + // gaussianGenerator -> WriteXMLToTestTheCuboid(); + + + } - MITK_TEST_END() -} \ No newline at end of file + // MITK_TEST_END() + +} diff --git a/Modules/ImageStatistics/itkMultiGaussianImageSource.h b/Modules/ImageStatistics/itkMultiGaussianImageSource.h index 71fd32c17a..41970adb52 100644 --- a/Modules/ImageStatistics/itkMultiGaussianImageSource.h +++ b/Modules/ImageStatistics/itkMultiGaussianImageSource.h @@ -1,291 +1,323 @@ /*========================================================================= * * 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" +#include + 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-the gaussian bell curve deviation 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. Furthermore it can calculate the maximal und minimal pixel intensities and whose indices in the founded sphere. \section algorithm_description Algorithm description \subsection gaus_generation Generation of a multigauss image A multigauss function consists of the sum of \f$N\f$ gauss function. The \f$i\f$-th \linebreak (\f$1 \leq i \leq N\f$) gaussian is described with the following seven parameters: 1. \f$x_0^{(i)}\f$ is the expectancy value in the \f$x\f$-Axis 2. \f$y_0^{(i)}\f$ is the expectancy value in the \f$y\f$-Axis 3. \f$z_0^{(i)}\f$ is the expectancy value in the \f$z\f$-Axis 4. \f$\sigma_x^{(i)}\f$ is the deviation in the \f$x\f$-Axis 5. \f$\sigma_y^{(i)}\f$ is the deviation in the \f$y\f$-Axis 6. \f$\sigma_z^{(i)}\f$ is the deviation in the \f$z\f$-Axis 7. \f$a^{(i)}\f$ is the altitude of the gaussian. A gauss function has the following form: \f{eqnarray*}{ f^{(i)}(x,y,z) = a^{(i)} exp \left[ - \left( \frac{(x - x_0^{(i)})^2}{2 (\sigma_x^{(i)})^2} + \frac{(y - y_0^{(i)})^2}{2 (\sigma_y^{(i)})^2} + \frac{(z - z_0^{(i)})^2}{2 (\sigma_z^{(i)})^2} \right) \right]. \f} A multigauss function has then the form: \f{align*}{ f_N(x,y,z) =& \sum_{i = 1}^{N}f^{(i)}(x,y,z)\\ =&\sum_{i=1}^{N} a^{(i)} exp \left[ - \left( \frac{(x - x_0^{(i)})^2}{2 (\sigma_x^{(i)})^2} + \frac{(y - y_0^{(i)})^2}{2 (\sigma_y^{(i)})^2} + \frac{(z - z_0^{(i)})^2}{2 (\sigma_z^{(i)})^2} \right) \right]. \f} The multigauss function \f$f_N\f$ will be evaluated at each voxel coordinate to become the voxel intensity. \subsection calculating_statistics Algorithm for calculating statiscic in a sphere This section explains how we can find a sphere region which has a maximal mean value over all sphere regions with a fixed radius. Furthermore we want to calculate the maximal and minimal value in the wanted sphere. To calculate the mean value in a sphere we use the simple definition of a mean value of a function in a bounded discrete space: the sum of the values of the function divided by the number of the summand. We discretize the spherical region of interest and use the spherical coordinate system to describe each point \linebreak [\f$x,y,z\f$] in this region: \f{eqnarray*}{ x = r \sin(\phi) \cos(\psi)\hspace{20pt} y = r \sin(\phi) \sin(\psi)\hspace{20pt} z = r \cos(\phi), \f} where \f$r \in [0;R]\f$, \f$\phi \in [0; \pi]\f$, \f$ \psi \in [-\pi; \pi)\f$. We decided to have only one discretizing parameter T, which is equal to the number of steps we make to traverse the radius \f$R\f$ of the sphere. So the higher the value of T is, the more accurate the mean value is. The two angle \f$\psi\f$ and \f$\phi\f$ have an equal step size, which is calculate as follows: Define a constant distance (\f$dist\f$) between the points on the sphere with the smallest radius. In our class we set \f$dist\f$ such a value, that on the equator of the smallest sphere only four points lie. Then we calculate for each radius the angle step size in such a way, that the points on the equator of the particular sphere have always the same distance and this should be approximately the same as \f$dist\f$ ( exactly the same distance as \f$dist\f$ is not always possible, because the length of the equator is generally not a multiple of a natural number). With such a discretization we almost achieve a uniformly distribution of the points in the "hot spot" region. So the following term equals the mean value in the sphere with midpoint [\f$x_l, y_l, z_l\f$] for a multigauss function: \f{align*}{ m_l = & \frac{1}{\kappa} \sum_{i = 1}^{T} \sum_{j = 0}^{\eta(i)} \sum_{k = 0}^{2\eta(i)} \sum_{i = 1}^{N} a^{(i)} exp \Big[ - \Big( \frac{(r_i \sin(\phi_j) \cos(\psi_k) - x_0^{(i)} + x_l)^2}{2 (\sigma_x^{(i)})^2} \\ &+\frac{(r_i \sin(\phi_j) \sin(\psi_k) - y_0^{(i)} + y_l)^2}{2 (\sigma_y^{(i)})^2} +\frac{(r_i \cos(\phi_j) - z_0^{(i)} + z_l)^2}{2 (\sigma_z^{(i)})^2}\Big)\Big] \f} Here denotes \f$\kappa\f$ the number of summand; the radius of the \f$i\f$-th sphere equals its index multiplied by the radius step size and the angle \f$\phi_j\f$ equals its index multiplied by the angle step size (analogical \f$\psi_j\f$): \f{align*}{ r_i = i \frac{R}{T}, \hspace{20pt} \phi_j = j \frac{\pi}{\eta(i)}, \hspace{20pt} \psi_k = \pi + k \frac{\pi}{\eta(k)}, \f} where \f$\eta(i)\f$ is the number of angle steps: \f$\eta(i)= \left \lfloor \frac{\pi r_i}{dist} \right \rfloor\f$ and \f$dist\f$ the fixed distance between two points: \f$dist = \frac{\pi R}{2T}\f$. We specify the volume of the "hot spot" to be \f$1\f$ cm\f$^3\f$. The radius of the sphere is there for: \f{align*}{ R = \left(\frac{3}{4\pi}\right)^{1/3}. \f} Now we know how to calculate the mean value in a sphere for given midpoint and radius of the wanted region. So we assume each voxel in the given image to be the sphere midpoint and we calculate the mean value as described above. If the new mean value is greater then the "maximum-until-now", we take the new value to be the "maximum-until-now". Then we go to the next voxel and make the same calculation and so on. At the same time we save the coordinates of the midpoint voxel. After we found the midpoint and the maximal mean value, we can calculate the maximum and the minimum in the sphere: we just traverse all the voxels in the region and take the maximum and minimum value and the respective coordinates. We can also optimize the calculation of the mean value in the founded sphere just by increasing the number of radius steps and calculate this value one more time. We set the new step number to be four time greater then T. */ 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 typename itk::ImageRegion<3> ::SizeValueType SizeRegionType; /** 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 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; + typedef MapContainer MapContainerPoints; + typedef MapContainer MapContainerEdges; + typedef MapContainer MapContainerRadius; + /** 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; /** Set the number of gaussian function*/ 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; /** Set the number of steps to traverse the radius */ 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 OutputImagePixelType GetMaxMeanValue() const; /** Get the index of the midpoint of a sphere with the maximal mean value */ virtual const IndexType 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(); /** Calculates and set the index an the value of maximulm and minimum in the wanted sphere*/ virtual void CalculateMaxAndMinInSphere(); /** Get the index in the sphere with maximal value*/ virtual const IndexType GetMaxValueIndexInSphere() const; /** Get the maximal value in the sphere*/ virtual const OutputImagePixelType GetMaxValueInSphere() const; /** Get the index in the sphere with minimal value*/ virtual const IndexType GetMinValueIndexInSphere() const; /** Get the minimal value in the sphere*/ virtual const OutputImagePixelType GetMinValueInSphere() const; /** Set the region of interest */ virtual void SetRegionOfInterest(ItkVectorType, ItkVectorType); /** Optimize the mean value in the wanted sphere*/ virtual void OptimizeMeanValue(); + /** Write a .mps file to visualise the point in the sphere*/ + virtual void WriteXMLToTest(); + virtual void WriteXMLToTestTheCuboid(); + virtual void WriteXMLToTestTheCuboidInsideTheSphere(); + virtual void CalculateTheMidPointMeanValueInCuboid(); + virtual void CalculateEdgesInSphere( PointType globalCoordinateMidpointCuboid, PointType globalCoordinateMidpointSphere, double cuboidRadius); + virtual double MultiGaussianFunctionValueAtCuboid(double xMin, double xMax, double yMin, double yMax, double zMin, double zMax); + virtual void InsertPoints( PointType globalCoordinateMidpointCuboid, double cuboidRadius); + virtual void GenerateCuboidSegmentationInSphere(); + virtual double FunctionPhi(double value); + + virtual void CalculateMidpoint(); + + virtual unsigned int IntesectTheSphere( PointType globalCoordinateMidpointCuboid, PointType globalCoordinateMidpointSphere, double sideLength); + void SetNormalDistributionValues(); + + double Quadtrees( PointType globalCoordinateMidpointCuboid, PointType globalCoordinateMidpointSphere, double sideLength, double meanValueTemp); + + /** Set the minimum possible pixel value. By default, it is * NumericTraits::min(). */ itkSetClampMacro( Min, OutputImagePixelType, NumericTraits< OutputImagePixelType >::NonpositiveMin(), NumericTraits< OutputImagePixelType >::max() ); + /** Check if a index is inside the image*/ + bool IsInImage(IndexType index); /** 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 OutputImagePixelType m_MaxValueInSphere; //maximal value in the wanted sphere IndexType m_MaxValueIndexInSphere; //index of the maximal value in the wanted sphere OutputImagePixelType m_MinValueInSphere; //minimal value in the wanted sphere IndexType m_MinValueIndexInSphere; //index of the minimal value in the wanted sphere 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 OutputImagePixelType m_ValueAtMidpoint; //value at the midpoint of the wanted sphere IndexType 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 ItkVectorType m_RegionOfInterestMax; // maximal values for the coordinates in the region of interest ItkVectorType m_RegionOfInterestMin; // minimal values for the coordinates in the region of interest typename TOutputImage::PixelType m_Min; //minimum possible value typename TOutputImage::PixelType m_Max; //maximum possible value PointType m_GlobalCoordinate; // physical coordiante of the sphere midpoint + + + MapContainerEdges m_Edges; + MapContainerPoints m_Midpoints, m_EdgePoints; + MapContainerRadius m_RadiusCuboid; + + VectorType m_XCoordToTest, m_YCoordToTest, m_ZCoordToTest; + double m_NormalDistValues [410]; // 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/itkMultiGaussianImageSource.hxx b/Modules/ImageStatistics/itkMultiGaussianImageSource.hxx index a0deb01ebd..2bee070eeb 100644 --- a/Modules/ImageStatistics/itkMultiGaussianImageSource.hxx +++ b/Modules/ImageStatistics/itkMultiGaussianImageSource.hxx @@ -1,589 +1,1640 @@ /*========================================================================= - * - * 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. - * - *=========================================================================*/ +* +* 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. - * - *=========================================================================*/ +* +* 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++ ) + /** + * + */ + 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_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_NumberOfGaussians = 0; + m_Radius = 1; + m_RadiusStepNumber = 5; + m_MeanValue = 0; - m_Min = NumericTraits< OutputImagePixelType >::NonpositiveMin(); - m_Max = NumericTraits< OutputImagePixelType >::max(); -} + m_Min = NumericTraits< OutputImagePixelType >::NonpositiveMin(); + m_Max = NumericTraits< OutputImagePixelType >::max(); + } -template< class TOutputImage > -MultiGaussianImageSource< TOutputImage > -::~MultiGaussianImageSource() -{} + template< class TOutputImage > + MultiGaussianImageSource< TOutputImage > + ::~MultiGaussianImageSource() + {} -template< class TOutputImage > -void -MultiGaussianImageSource< TOutputImage > -::SetSize(SizeValueArrayType sizeArray) -{ - const unsigned int count = TOutputImage::ImageDimension; - unsigned int i; + template< class TOutputImage > + void + MultiGaussianImageSource< TOutputImage > + ::SetSize(SizeValueArrayType sizeArray) + { + const unsigned int count = TOutputImage::ImageDimension; + unsigned int i; - for ( i = 0; i < count; i++ ) + for ( i = 0; i < count; i++ ) { - if ( sizeArray[i] != this->m_Size[i] ) + if ( sizeArray[i] != this->m_Size[i] ) { - break; + break; } } - if ( i < count ) + if ( i < count ) { - this->Modified(); - for ( i = 0; i < count; i++ ) + this->Modified(); + for ( i = 0; i < count; i++ ) { - this->m_Size[i] = sizeArray[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 > + 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; + template< class TOutputImage > + void + MultiGaussianImageSource< TOutputImage > + ::SetSpacing(SpacingValueArrayType spacingArray) + { + const unsigned int count = TOutputImage::ImageDimension; + unsigned int i; - for ( i = 0; i < count; i++ ) + for ( i = 0; i < count; i++ ) { - if ( spacingArray[i] != this->m_Spacing[i] ) + if ( spacingArray[i] != this->m_Spacing[i] ) { - break; + break; } } - if ( i < count ) + if ( i < count ) { - this->Modified(); - for ( i = 0; i < count; i++ ) + this->Modified(); + for ( i = 0; i < count; i++ ) { - this->m_Spacing[i] = spacingArray[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; + template< class TOutputImage > + void + MultiGaussianImageSource< TOutputImage > + ::SetOrigin(PointValueArrayType originArray) + { + const unsigned int count = TOutputImage::ImageDimension; + unsigned int i; - for ( i = 0; i < count; i++ ) + for ( i = 0; i < count; i++ ) { - if ( originArray[i] != this->m_Origin[i] ) + if ( originArray[i] != this->m_Origin[i] ) { - break; + break; } } - if ( i < count ) + if ( i < count ) { - this->Modified(); - for ( i = 0; i < count; i++ ) + this->Modified(); + for ( i = 0; i < count; i++ ) { - this->m_Origin[i] = originArray[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++ ) + 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]; + this->m_OriginArray[i] = this->m_Origin[i]; } - return this->m_OriginArray; -} + 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++ ) + 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]; + this->m_SpacingArray[i] = this->m_Spacing[i]; } - return this->m_SpacingArray; -} + 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; + /** + * + */ + 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 << indent << "Origin: ["; + unsigned int ii = 0; + while( ii < TOutputImage::ImageDimension - 1 ) { - os << m_Origin[ii] << ", "; - ++ii; + os << m_Origin[ii] << ", "; + ++ii; } - os << m_Origin[ii] << "]" << std::endl; + os << m_Origin[ii] << "]" << std::endl; - os << indent << "Spacing: ["; - ii = 0; - while( ii < TOutputImage::ImageDimension - 1 ) + os << indent << "Spacing: ["; + ii = 0; + while( ii < TOutputImage::ImageDimension - 1 ) { - os << m_Spacing[ii] << ", "; - ++ii; + os << m_Spacing[ii] << ", "; + ++ii; } - os << m_Spacing[ii] << "]" << std::endl; + os << m_Spacing[ii] << "]" << std::endl; - os << indent << "Size: ["; - ii = 0; - while( ii < TOutputImage::ImageDimension - 1 ) + os << indent << "Size: ["; + ii = 0; + while( ii < TOutputImage::ImageDimension - 1 ) { - os << m_Size[ii] << ", "; - ++ii; + os << m_Size[ii] << ", "; + ++ii; } - os << m_Size[ii] << "]" << std::endl; -} + os << m_Size[ii] << "]" << std::endl; + } -template< class TOutputImage > -unsigned int -MultiGaussianImageSource< TOutputImage > -::GetNumberOfGaussians() const -{ - return this->m_NumberOfGaussians; -} + 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 > + 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 > + 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 > + 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 > + ::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 > + void + MultiGaussianImageSource< TOutputImage > + ::SetNumberOfGausssians( unsigned int n ) + { + this->m_NumberOfGaussians = n; + } -template< class TOutputImage > -void -MultiGaussianImageSource< TOutputImage > -::SetRegionOfInterest( ItkVectorType roiMin, ItkVectorType roiMax ) -{ - m_RegionOfInterestMax.operator=(roiMax); - m_RegionOfInterestMin.operator=(roiMin); -} + template< class TOutputImage > + void + MultiGaussianImageSource< TOutputImage > + ::SetRegionOfInterest( ItkVectorType roiMin, ItkVectorType roiMax ) + { + m_RegionOfInterestMax.operator=(roiMax); + m_RegionOfInterestMin.operator=(roiMin); + } -template< class TOutputImage > -const typename MultiGaussianImageSource< TOutputImage >::OutputImagePixelType -MultiGaussianImageSource< TOutputImage > -::GetMaxMeanValue() const -{ - return m_MeanValue; -} + template< class TOutputImage > + const typename MultiGaussianImageSource< TOutputImage >::OutputImagePixelType + MultiGaussianImageSource< TOutputImage > + ::GetMaxMeanValue() const + { + return m_MeanValue; + } -template< class TOutputImage > -const typename MultiGaussianImageSource< TOutputImage >::OutputImagePixelType -MultiGaussianImageSource< TOutputImage > -::GetMaxValueInSphere() const -{ - return m_MaxValueInSphere; -} + template< class TOutputImage > + const typename MultiGaussianImageSource< TOutputImage >::OutputImagePixelType + MultiGaussianImageSource< TOutputImage > + ::GetMaxValueInSphere() const + { + return m_MaxValueInSphere; + } -template< class TOutputImage > -const typename MultiGaussianImageSource< TOutputImage >::IndexType -MultiGaussianImageSource< TOutputImage > -::GetMaxValueIndexInSphere() const -{ - return m_MaxValueIndexInSphere; -} + template< class TOutputImage > + const typename MultiGaussianImageSource< TOutputImage >::IndexType + MultiGaussianImageSource< TOutputImage > + ::GetMaxValueIndexInSphere() const + { + return m_MaxValueIndexInSphere; + } -template< class TOutputImage > -const typename MultiGaussianImageSource< TOutputImage >::OutputImagePixelType -MultiGaussianImageSource< TOutputImage > -::GetMinValueInSphere() const -{ - return m_MinValueInSphere; -} + template< class TOutputImage > + const typename MultiGaussianImageSource< TOutputImage >::OutputImagePixelType + MultiGaussianImageSource< TOutputImage > + ::GetMinValueInSphere() const + { + return m_MinValueInSphere; + } -template< class TOutputImage > -const typename MultiGaussianImageSource< TOutputImage >::IndexType -MultiGaussianImageSource< TOutputImage > -::GetMinValueIndexInSphere() const -{ - return m_MinValueIndexInSphere; -} -//---------------------------------------------------------------------------- -template< class TOutputImage > -const typename MultiGaussianImageSource< TOutputImage >::IndexType -MultiGaussianImageSource< TOutputImage > -::GetSphereMidpoint() const -{ - return m_SphereMidpoint; -} + template< class TOutputImage > + const typename MultiGaussianImageSource< TOutputImage >::IndexType + MultiGaussianImageSource< TOutputImage > + ::GetMinValueIndexInSphere() const + { + return m_MinValueIndexInSphere; + } + //---------------------------------------------------------------------------- + template< class TOutputImage > + const typename MultiGaussianImageSource< TOutputImage >::IndexType + MultiGaussianImageSource< TOutputImage > + ::GetSphereMidpoint() const + { + return m_SphereMidpoint; + } -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) + //---------------------------------------------------------------------------------------------------------------------- + + template< class TOutputImage > + double + MultiGaussianImageSource< TOutputImage > + ::MultiGaussianFunctionValueAtCuboid(double xMin, double xMax, double yMin, double yMax, double zMin, double zMax) { - summand0 = (x - m_CenterX[n]) / m_SigmaX[n]; - summand1 = (y - m_CenterY[n]) / m_SigmaY[n]; - summand2 = (z - m_CenterZ[n]) / m_SigmaZ[n]; + double mean = 0; + double summand0, summand1, summand2, value, factor; - power = summand0 * summand0 + summand1 * summand1 + summand2 * summand2; - value = value + m_Altitude[n] * pow(itk::Math::e, -0.5 * power); + for(unsigned int n =0; n < m_NumberOfGaussians; ++n) + { + summand0 = FunctionPhi( (xMax - m_CenterX[n]) / m_SigmaX[n] ) + FunctionPhi( (xMin - m_CenterX[n]) / m_SigmaX[n] ); + summand1 = FunctionPhi( (yMax - m_CenterY[n]) / m_SigmaY[n] ) + FunctionPhi( (yMin - m_CenterY[n]) / m_SigmaY[n] ); + summand2 = FunctionPhi( (zMax - m_CenterZ[n]) / m_SigmaZ[n] ) + FunctionPhi( (zMin - m_CenterZ[n]) / m_SigmaZ[n] ); + + value = summand0 * summand1 * summand2; + factor = (m_SigmaX[n] * m_SigmaY[n] * m_SigmaZ[n] * pow(2.0 * itk::Math::pi, 1.5 )); // ((xMax - xMin) * (yMax - yMin) * (zMax - zMin)); + mean = mean + m_Altitude[n] * factor * value; // * m_Altitude[n] ???? + + } + return mean; } - 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: + + template< class TOutputImage > + double + MultiGaussianImageSource< TOutputImage > + ::FunctionPhi(double value) + { + double phiAtValue; + + // TODO qubische interpolation + int indexValue = static_cast( 100 * value); + if (value < 0.0) + { + phiAtValue = 1.0 - m_NormalDistValues[- indexValue + 1]; + } + if(indexValue > 410) + { + phiAtValue = 1; + } + phiAtValue = m_NormalDistValues[indexValue]; + return phiAtValue; + } + //---------------------------------------------------------------------------------------------------------------------- + template< class TOutputImage > + void + MultiGaussianImageSource< TOutputImage > + ::InsertPoints( PointType globalCoordinateMidpointCuboid, double cuboidRadius) + { + //TODO letzter Punkt entfernen.....! + + PointType edgePoint; + MapContainerPoints::ElementIdentifier identifier; + VectorType edgesIds; + MapContainerPoints::ElementIdentifier id = m_Edges.Size() + 1; + MapContainerPoints::ElementIdentifier id1 = m_Midpoints.Size() + 1; + for(int i = -1; i < 2; i+=2) + { + for(int k = -1; k < 2; k+=2) + { + for(int j = -1; j < 2; j+=2) + { + edgePoint[0] = globalCoordinateMidpointCuboid[0] + i * cuboidRadius; + edgePoint[1] = globalCoordinateMidpointCuboid[1] + k * cuboidRadius; + edgePoint[2] = globalCoordinateMidpointCuboid[2] + j * cuboidRadius; + + identifier = m_EdgePoints.Size() + 1; + m_EdgePoints.InsertElement(identifier, globalCoordinateMidpointCuboid); + edgesIds.push_back(identifier); + } + } + } + + m_Midpoints.InsertElement(id1, globalCoordinateMidpointCuboid); + m_RadiusCuboid.InsertElement(id1, cuboidRadius); + m_Edges.InsertElement(id, edgesIds); + + } + + //---------------------------------------------------------------------------------------------------------------------- + template< class TOutputImage > + void + MultiGaussianImageSource< TOutputImage > + ::CalculateEdgesInSphere( PointType globalCoordinateMidpointCuboid, PointType globalCoordinateMidpointSphere, double cuboidRadius) + { + double cuboidRadiusRecursion = cuboidRadius; + double minCuboidRadius = m_Spacing[0] / 5.0; + bool exit = 0; + while(cuboidRadiusRecursion > minCuboidRadius && !exit) + { + int intersect = IntesectTheSphere( globalCoordinateMidpointCuboid, globalCoordinateMidpointSphere, cuboidRadiusRecursion); + + if( intersect == 1) + { + cuboidRadiusRecursion = cuboidRadiusRecursion / 2.0; + for(int i = -1; i < 2; i+=2) + { + for(int k = -1; k < 2; k+=2) + { + for(int j = -1; j < 2; j+=2) + { + // cuboid intersect the sphere + PointType newMidpoint; + + newMidpoint[0] = globalCoordinateMidpointCuboid[0] + static_cast(i) * cuboidRadiusRecursion; + newMidpoint[1] = globalCoordinateMidpointCuboid[1] + static_cast(k) * cuboidRadiusRecursion; + newMidpoint[2] = globalCoordinateMidpointCuboid[2] + static_cast(j) * cuboidRadiusRecursion; + + CalculateEdgesInSphere( newMidpoint, globalCoordinateMidpointSphere, cuboidRadiusRecursion); + } + } + } + } + else if(intersect == 2) + { + // cuboid in the sphere + // insert the edge points of the cuboid + InsertPoints(globalCoordinateMidpointCuboid, cuboidRadiusRecursion); + exit = 1; + + // WriteXMLToTestTheCuboidInsideTheSphere(); + + } + else if (intersect == 0) + { + // cuboid not in the sphere + exit = 1; + } + + } + + } + + + //----------------------------------------------------------------------------------------------------------------------- + + template< class TOutputImage > + void + MultiGaussianImageSource< TOutputImage > + ::GenerateCuboidSegmentationInSphere() + { + + double cuboidRadius; + + + PointType globalCoordinateMidpointCuboid; + PointType globalCoordinateMidpointSphere; + IndexType index; + OutputImageRegionType regionOfInterest; + IndexType indexR; + indexR.SetElement(0, m_RegionOfInterestMin[0]); + indexR.SetElement(1, m_RegionOfInterestMin[1]); + indexR.SetElement(2, m_RegionOfInterestMin[2]); + regionOfInterest.SetIndex(indexR); + SizeType sizeROI; + sizeROI.SetElement(0, m_RegionOfInterestMax[0] - m_RegionOfInterestMin[0] +1); + sizeROI.SetElement(1, m_RegionOfInterestMax[1] - m_RegionOfInterestMin[1] +1); + sizeROI.SetElement(2, m_RegionOfInterestMax[2] - m_RegionOfInterestMin[2] +1); + regionOfInterest.SetSize(sizeROI); + typename TOutputImage::Pointer image = this->GetOutput(0); + IteratorType regionOfInterestIterator(image, regionOfInterest); + + for(regionOfInterestIterator.GoToBegin(); !regionOfInterestIterator.IsAtEnd(); ++regionOfInterestIterator) + { + cuboidRadius = m_Radius / 2.0; + index = regionOfInterestIterator.GetIndex(); + image->TransformIndexToPhysicalPoint(index, globalCoordinateMidpointSphere); + for(int i = -1; i < 2; i++) + { + for(int k = -1; k < 2; k++) + { + for(int j = -1; j < 2; j++) + { + PointType newMidpoint; + newMidpoint[0] = globalCoordinateMidpointSphere[0] + static_cast(i) * cuboidRadius; + newMidpoint[1] = globalCoordinateMidpointSphere[1] + static_cast(k) * cuboidRadius; + newMidpoint[2] = globalCoordinateMidpointSphere[2] + static_cast(j) * cuboidRadius; + CalculateEdgesInSphere( newMidpoint, globalCoordinateMidpointSphere, cuboidRadius); + } + } + } + + + + + WriteXMLToTestTheCuboidInsideTheSphere(); + //to test only one step!!! TODO + regionOfInterestIterator.GoToReverseBegin(); + } + + } + + + //---------------------------------------------------------------------------------------------------------------------- + + template< class TOutputImage > + void + MultiGaussianImageSource< TOutputImage > + ::CalculateTheMidPointMeanValueInCuboid() + { + + m_MeanValue = 0.0; + double meanValueTemp; + PointType midpoint; + MapContainerPoints::ElementIdentifier cuboidNumber = m_Midpoints.Size(); + SetNormalDistributionValues(); + double radius; + double xMin, xMax, yMin, yMax, zMin, zMax; + + + PointType globalCoordinateMidpointSphere; + IndexType index; + OutputImageRegionType regionOfInterest; + IndexType indexR; + indexR.SetElement(0, m_RegionOfInterestMin[0]); + indexR.SetElement(1, m_RegionOfInterestMin[1]); + indexR.SetElement(2, m_RegionOfInterestMin[2]); + regionOfInterest.SetIndex(indexR); + SizeType sizeROI; + sizeROI.SetElement(0, m_RegionOfInterestMax[0] - m_RegionOfInterestMin[0] +1); + sizeROI.SetElement(1, m_RegionOfInterestMax[1] - m_RegionOfInterestMin[1] +1); + sizeROI.SetElement(2, m_RegionOfInterestMax[2] - m_RegionOfInterestMin[2] +1); + regionOfInterest.SetSize(sizeROI); + typename TOutputImage::Pointer image = this->GetOutput(0); + IteratorType regionOfInterestIterator(image, regionOfInterest); + + for(regionOfInterestIterator.GoToBegin(); !regionOfInterestIterator.IsAtEnd(); ++regionOfInterestIterator) + { + + index = regionOfInterestIterator.GetIndex(); + image->TransformIndexToPhysicalPoint(index, globalCoordinateMidpointSphere); + + meanValueTemp = 0.0; + for(unsigned int i = 0; i < cuboidNumber; i++ ) + { + + midpoint = m_Midpoints.ElementAt(i); + midpoint[0] = midpoint[0] + globalCoordinateMidpointSphere[0]; + midpoint[2] = midpoint[2] + globalCoordinateMidpointSphere[2]; + midpoint[1] = midpoint[1] + globalCoordinateMidpointSphere[1]; + radius= m_RadiusCuboid.ElementAt(i); + xMin = midpoint[0] - radius; + xMax = midpoint[0] + radius; + yMin = midpoint[1] - radius; + yMax = midpoint[1] + radius; + zMin = midpoint[2] - radius; + zMax = midpoint[2] + radius; + meanValueTemp = meanValueTemp + MultiGaussianFunctionValueAtCuboid( xMin, xMax, yMin, yMax, zMin, zMax ); + } + + meanValueTemp = meanValueTemp / (m_Radius * itk::Math::pi * itk::Math::pi); + if(meanValueTemp > m_MeanValue) + { + m_GlobalCoordinate = globalCoordinateMidpointSphere; + m_MeanValue = meanValueTemp; + m_SphereMidpoint = index; + } + + } + + } + + //----------------------------------------------------------------------------------------------------------------------- + + template< class TOutputImage > + void + MultiGaussianImageSource< TOutputImage > + ::CalculateMidpoint() + { + double valueMean = 0; + m_MeanValue = 0.0; + double sideLength = m_Radius; + double cuboidRadius; + double sideLengthNew = m_Radius; + double meanValueTemp; + bool intersect; + PointType globalCoordinateMidpointCuboid; + PointType globalCoordinateMidpointSphere; + IndexType index; + OutputImageRegionType regionOfInterest; + IndexType indexR; + indexR.SetElement(0, m_RegionOfInterestMin[0]); + indexR.SetElement(1, m_RegionOfInterestMin[1]); + indexR.SetElement(2, m_RegionOfInterestMin[2]); + regionOfInterest.SetIndex(indexR); + SizeType sizeROI; + sizeROI.SetElement(0, m_RegionOfInterestMax[0] - m_RegionOfInterestMin[0] +1); + sizeROI.SetElement(1, m_RegionOfInterestMax[1] - m_RegionOfInterestMin[1] +1); + sizeROI.SetElement(2, m_RegionOfInterestMax[2] - m_RegionOfInterestMin[2] +1); + regionOfInterest.SetSize(sizeROI); + typename TOutputImage::Pointer image = this->GetOutput(0); + IteratorType regionOfInterestIterator(image, regionOfInterest); + SetNormalDistributionValues(); + + for(regionOfInterestIterator.GoToBegin(); !regionOfInterestIterator.IsAtEnd(); ++regionOfInterestIterator) + { + cuboidRadius = sideLength / 2.0; + meanValueTemp = 0.0; + index = regionOfInterestIterator.GetIndex(); + image->TransformIndexToPhysicalPoint(index, globalCoordinateMidpointSphere); + + globalCoordinateMidpointCuboid[0] = globalCoordinateMidpointSphere[0] - cuboidRadius; + globalCoordinateMidpointCuboid[1] = globalCoordinateMidpointSphere[1] - cuboidRadius; + globalCoordinateMidpointCuboid[2] = globalCoordinateMidpointSphere[2] - cuboidRadius; + meanValueTemp = meanValueTemp + Quadtrees(globalCoordinateMidpointCuboid, globalCoordinateMidpointSphere, cuboidRadius, meanValueTemp); + + globalCoordinateMidpointCuboid[0] = globalCoordinateMidpointSphere[0] + cuboidRadius; + globalCoordinateMidpointCuboid[1] = globalCoordinateMidpointSphere[1] - cuboidRadius; + globalCoordinateMidpointCuboid[2] = globalCoordinateMidpointSphere[2] - cuboidRadius; + meanValueTemp = meanValueTemp + Quadtrees(globalCoordinateMidpointCuboid, globalCoordinateMidpointSphere, cuboidRadius, meanValueTemp); + + globalCoordinateMidpointCuboid[0] = globalCoordinateMidpointSphere[0] + cuboidRadius; + globalCoordinateMidpointCuboid[1] = globalCoordinateMidpointSphere[1] + cuboidRadius; + globalCoordinateMidpointCuboid[2] = globalCoordinateMidpointSphere[2] - cuboidRadius; + meanValueTemp = meanValueTemp + Quadtrees(globalCoordinateMidpointCuboid, globalCoordinateMidpointSphere, cuboidRadius, meanValueTemp); + + globalCoordinateMidpointCuboid[0] = globalCoordinateMidpointSphere[0] + cuboidRadius; + globalCoordinateMidpointCuboid[1] = globalCoordinateMidpointSphere[1] + cuboidRadius; + globalCoordinateMidpointCuboid[2] = globalCoordinateMidpointSphere[2] + cuboidRadius; + meanValueTemp = meanValueTemp + Quadtrees(globalCoordinateMidpointCuboid, globalCoordinateMidpointSphere, cuboidRadius, meanValueTemp); + + globalCoordinateMidpointCuboid[0] = globalCoordinateMidpointSphere[0] + cuboidRadius; + globalCoordinateMidpointCuboid[1] = globalCoordinateMidpointSphere[1] - cuboidRadius; + globalCoordinateMidpointCuboid[2] = globalCoordinateMidpointSphere[2] + cuboidRadius; + meanValueTemp = meanValueTemp + Quadtrees(globalCoordinateMidpointCuboid, globalCoordinateMidpointSphere, + cuboidRadius, meanValueTemp); + + globalCoordinateMidpointCuboid[0] = globalCoordinateMidpointSphere[0] - cuboidRadius; + globalCoordinateMidpointCuboid[1] = globalCoordinateMidpointSphere[1] + cuboidRadius; + globalCoordinateMidpointCuboid[2] = globalCoordinateMidpointSphere[2] + cuboidRadius; + meanValueTemp = meanValueTemp + Quadtrees(globalCoordinateMidpointCuboid, globalCoordinateMidpointSphere, cuboidRadius, meanValueTemp); + + globalCoordinateMidpointCuboid[0] = globalCoordinateMidpointSphere[0] - cuboidRadius; + globalCoordinateMidpointCuboid[1] = globalCoordinateMidpointSphere[1] + cuboidRadius; + globalCoordinateMidpointCuboid[2] = globalCoordinateMidpointSphere[2] - cuboidRadius; + meanValueTemp = meanValueTemp + Quadtrees(globalCoordinateMidpointCuboid, globalCoordinateMidpointSphere, cuboidRadius, meanValueTemp); + + globalCoordinateMidpointCuboid[0] = globalCoordinateMidpointSphere[0] - cuboidRadius; + globalCoordinateMidpointCuboid[1] = globalCoordinateMidpointSphere[1] - cuboidRadius; + globalCoordinateMidpointCuboid[2] = globalCoordinateMidpointSphere[2] + cuboidRadius; + meanValueTemp = meanValueTemp + Quadtrees(globalCoordinateMidpointCuboid, globalCoordinateMidpointSphere, + cuboidRadius, meanValueTemp); + + + + meanValueTemp = meanValueTemp / (m_Radius * itk::Math::pi * itk::Math::pi); + if(meanValueTemp > m_MeanValue) + { + m_GlobalCoordinate = globalCoordinateMidpointSphere; + m_MeanValue = meanValueTemp; + m_SphereMidpoint = index; + } + + + //to test only one step!!! TODO + regionOfInterestIterator.GoToReverseBegin(); + + } + } + //---------------------------------------------------------------------------------------------------------------------- + + template< class TOutputImage > + double + MultiGaussianImageSource< TOutputImage > + ::Quadtrees( PointType globalCoordinateMidpointCuboid, PointType globalCoordinateMidpointSphere, double cuboidRadius, double meanValueTemp) + { + double xMin, xMax, yMin, yMax, zMin, zMax; + double minDistance = m_Spacing[0] / 2.0; + unsigned int intersect; + + while(cuboidRadius > minDistance) + { + intersect = IntesectTheSphere( globalCoordinateMidpointCuboid, globalCoordinateMidpointSphere, cuboidRadius); + if( intersect == 1 ) + { + cuboidRadius = cuboidRadius / 2.0; + globalCoordinateMidpointCuboid[0] = globalCoordinateMidpointCuboid[0] - cuboidRadius; + globalCoordinateMidpointCuboid[1] = globalCoordinateMidpointCuboid[1] - cuboidRadius; + globalCoordinateMidpointCuboid[2] = globalCoordinateMidpointCuboid[2] - cuboidRadius; + meanValueTemp = meanValueTemp + Quadtrees( globalCoordinateMidpointCuboid, globalCoordinateMidpointSphere, cuboidRadius, meanValueTemp); + + globalCoordinateMidpointCuboid[0] = globalCoordinateMidpointCuboid[0] - cuboidRadius; + globalCoordinateMidpointCuboid[1] = globalCoordinateMidpointCuboid[1] + cuboidRadius; + globalCoordinateMidpointCuboid[2] = globalCoordinateMidpointCuboid[2] + cuboidRadius; + meanValueTemp = meanValueTemp + Quadtrees( globalCoordinateMidpointCuboid, globalCoordinateMidpointSphere, cuboidRadius, meanValueTemp); + + globalCoordinateMidpointCuboid[0] = globalCoordinateMidpointCuboid[0] + cuboidRadius; + globalCoordinateMidpointCuboid[1] = globalCoordinateMidpointCuboid[1] - cuboidRadius; + globalCoordinateMidpointCuboid[2] = globalCoordinateMidpointCuboid[2] - cuboidRadius; + meanValueTemp = meanValueTemp + Quadtrees( globalCoordinateMidpointCuboid, globalCoordinateMidpointSphere, cuboidRadius, meanValueTemp); + + globalCoordinateMidpointCuboid[0] = globalCoordinateMidpointCuboid[0] + cuboidRadius; + globalCoordinateMidpointCuboid[1] = globalCoordinateMidpointCuboid[1] + cuboidRadius; + globalCoordinateMidpointCuboid[2] = globalCoordinateMidpointCuboid[2] - cuboidRadius; + meanValueTemp = meanValueTemp + Quadtrees( globalCoordinateMidpointCuboid, globalCoordinateMidpointSphere, cuboidRadius, meanValueTemp); + + globalCoordinateMidpointCuboid[0] = globalCoordinateMidpointCuboid[0] + cuboidRadius; + globalCoordinateMidpointCuboid[1] = globalCoordinateMidpointCuboid[1] - cuboidRadius; + globalCoordinateMidpointCuboid[2] = globalCoordinateMidpointCuboid[2] + cuboidRadius; + meanValueTemp = meanValueTemp + Quadtrees( globalCoordinateMidpointCuboid, globalCoordinateMidpointSphere, cuboidRadius, meanValueTemp); + + globalCoordinateMidpointCuboid[0] = globalCoordinateMidpointCuboid[0] - cuboidRadius; + globalCoordinateMidpointCuboid[1] = globalCoordinateMidpointCuboid[1] + cuboidRadius; + globalCoordinateMidpointCuboid[2] = globalCoordinateMidpointCuboid[2] - cuboidRadius; + meanValueTemp = meanValueTemp + Quadtrees( globalCoordinateMidpointCuboid, globalCoordinateMidpointSphere, cuboidRadius, meanValueTemp); + + globalCoordinateMidpointCuboid[0] = globalCoordinateMidpointCuboid[0] - cuboidRadius; + globalCoordinateMidpointCuboid[1] = globalCoordinateMidpointCuboid[1] - cuboidRadius; + globalCoordinateMidpointCuboid[2] = globalCoordinateMidpointCuboid[2] + cuboidRadius; + meanValueTemp = meanValueTemp + Quadtrees( globalCoordinateMidpointCuboid, globalCoordinateMidpointSphere, cuboidRadius, meanValueTemp); + + globalCoordinateMidpointCuboid[0] = globalCoordinateMidpointCuboid[0] + cuboidRadius; + globalCoordinateMidpointCuboid[1] = globalCoordinateMidpointCuboid[1] + cuboidRadius; + globalCoordinateMidpointCuboid[2] = globalCoordinateMidpointCuboid[2] + cuboidRadius; + meanValueTemp = meanValueTemp + Quadtrees( globalCoordinateMidpointCuboid, globalCoordinateMidpointSphere, cuboidRadius, meanValueTemp); + + + } + else if( intersect == 2 ) + { + + xMin = globalCoordinateMidpointCuboid[0] - cuboidRadius; + xMax = globalCoordinateMidpointCuboid[0] + cuboidRadius; + yMin = globalCoordinateMidpointCuboid[1] - cuboidRadius; + yMax = globalCoordinateMidpointCuboid[1] + cuboidRadius; + zMin = globalCoordinateMidpointCuboid[2] - cuboidRadius; + zMax = globalCoordinateMidpointCuboid[2] + cuboidRadius; + + meanValueTemp = meanValueTemp + MultiGaussianFunctionValueAtCuboid( xMin, xMax, yMin, yMax, zMin, zMax); + cuboidRadius = minDistance; + } + else if( intersect == 0 ) + { + cuboidRadius = minDistance; + } + } + return meanValueTemp; + } + //---------------------------------------------------------------------------------------------------------------------- + + template< class TOutputImage > + unsigned int + MultiGaussianImageSource< TOutputImage > + ::IntesectTheSphere( PointType globalCoordinateMidpointCuboid, PointType globalCoordinateMidpointSphere, double cuboidRadius) + { + + unsigned int intersect = 1; + PointType cuboidEdge1, cuboidEdge2, cuboidEdge3, cuboidEdge4, cuboidEdge5, cuboidEdge6, cuboidEdge7, cuboidEdge8; + int count = 0, index = 0; + + cuboidEdge1[0] = globalCoordinateMidpointCuboid[0] - cuboidRadius; + cuboidEdge1[1] = globalCoordinateMidpointCuboid[1] - cuboidRadius; + cuboidEdge1[2] = globalCoordinateMidpointCuboid[2] - cuboidRadius; + + cuboidEdge2[0] = globalCoordinateMidpointCuboid[0] + cuboidRadius; + cuboidEdge2[1] = globalCoordinateMidpointCuboid[1] - cuboidRadius; + cuboidEdge2[2] = globalCoordinateMidpointCuboid[2] - cuboidRadius; + + cuboidEdge3[0] = globalCoordinateMidpointCuboid[0] + cuboidRadius; + cuboidEdge3[1] = globalCoordinateMidpointCuboid[1] + cuboidRadius; + cuboidEdge3[2] = globalCoordinateMidpointCuboid[2] - cuboidRadius; + + cuboidEdge4[0] = globalCoordinateMidpointCuboid[0] + cuboidRadius; + cuboidEdge4[1] = globalCoordinateMidpointCuboid[1] - cuboidRadius; + cuboidEdge4[2] = globalCoordinateMidpointCuboid[2] + cuboidRadius; + + cuboidEdge5[0] = globalCoordinateMidpointCuboid[0] - cuboidRadius; + cuboidEdge5[1] = globalCoordinateMidpointCuboid[1] + cuboidRadius; + cuboidEdge5[2] = globalCoordinateMidpointCuboid[2] - cuboidRadius; + + cuboidEdge6[0] = globalCoordinateMidpointCuboid[0] - cuboidRadius; + cuboidEdge6[1] = globalCoordinateMidpointCuboid[1] + cuboidRadius; + cuboidEdge6[2] = globalCoordinateMidpointCuboid[2] + cuboidRadius; + + cuboidEdge7[0] = globalCoordinateMidpointCuboid[0] - cuboidRadius; + cuboidEdge7[1] = globalCoordinateMidpointCuboid[1] - cuboidRadius; + cuboidEdge7[2] = globalCoordinateMidpointCuboid[2] + cuboidRadius; + + cuboidEdge8[0] = globalCoordinateMidpointCuboid[0] + cuboidRadius; + cuboidEdge8[1] = globalCoordinateMidpointCuboid[1] + cuboidRadius; + cuboidEdge8[2] = globalCoordinateMidpointCuboid[2] + cuboidRadius; + + + + count = (globalCoordinateMidpointSphere.EuclideanDistanceTo(cuboidEdge1) < m_Radius) + ( globalCoordinateMidpointSphere.EuclideanDistanceTo(cuboidEdge2) < m_Radius) + ( globalCoordinateMidpointSphere.EuclideanDistanceTo(cuboidEdge3) < m_Radius) + ( globalCoordinateMidpointSphere.EuclideanDistanceTo(cuboidEdge4) < m_Radius) + ( globalCoordinateMidpointSphere.EuclideanDistanceTo(cuboidEdge5) < m_Radius) + ( globalCoordinateMidpointSphere.EuclideanDistanceTo(cuboidEdge6) < m_Radius) + + ( globalCoordinateMidpointSphere.EuclideanDistanceTo(cuboidEdge7) < m_Radius) + ( globalCoordinateMidpointSphere.EuclideanDistanceTo(cuboidEdge8) < m_Radius); + if ( count == 0 ) + { + intersect = 0; + } + if (count == 8 ) + { + intersect = 2; + } + + //if(count != 0 && count != 8) + //{ + // //to write a point set + // m_XCoordToTest.push_back(cuboidEdge1[0]); + // m_YCoordToTest.push_back(cuboidEdge1[1]); + // m_ZCoordToTest.push_back(cuboidEdge1[2]); + // //to write a point set + // m_XCoordToTest.push_back(cuboidEdge2[0]); + // m_YCoordToTest.push_back(cuboidEdge2[1]); + // m_ZCoordToTest.push_back(cuboidEdge2[2]); + // //to write a point set + // m_XCoordToTest.push_back(cuboidEdge3[0]); + // m_YCoordToTest.push_back(cuboidEdge3[1]); + // m_ZCoordToTest.push_back(cuboidEdge3[2]); + // //to write a point set + // m_XCoordToTest.push_back(cuboidEdge4[0]); + // m_YCoordToTest.push_back(cuboidEdge4[1]); + // m_ZCoordToTest.push_back(cuboidEdge4[2]); + // //to write a point set + // m_XCoordToTest.push_back(cuboidEdge5[0]); + // m_YCoordToTest.push_back(cuboidEdge5[1]); + // m_ZCoordToTest.push_back(cuboidEdge5[2]); + // //to write a point set + // m_XCoordToTest.push_back(cuboidEdge6[0]); + // m_YCoordToTest.push_back(cuboidEdge6[1]); + // m_ZCoordToTest.push_back(cuboidEdge6[2]); + // //to write a point set + // m_XCoordToTest.push_back(cuboidEdge7[0]); + // m_YCoordToTest.push_back(cuboidEdge7[1]); + // m_ZCoordToTest.push_back(cuboidEdge7[2]); + // //to write a point set + // m_XCoordToTest.push_back(cuboidEdge8[0]); + // m_YCoordToTest.push_back(cuboidEdge8[1]); + // m_ZCoordToTest.push_back(cuboidEdge8[2]); + // WriteXMLToTestTheCuboid(); + + //} + + return intersect; + } + + //---------------------------------------------------------------------------------------------------------------------- + + 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 region of interest 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()) + 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 x, y, z, temp; - MultiGaussianImageSource< TOutputImage >::OutputImagePixelType value, meanValueTemp, valueAdd; - PointType globalCoordinate; - IndexType index; - 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; - OutputImageRegionType regionOfInterest; - IndexType indexR; - indexR.SetElement(0, m_RegionOfInterestMin[0]); - indexR.SetElement(1, m_RegionOfInterestMin[1]); - indexR.SetElement(2, m_RegionOfInterestMin[2]); - regionOfInterest.SetIndex(indexR); - SizeType sizeROI; - sizeROI.SetElement(0, m_RegionOfInterestMax[0] - m_RegionOfInterestMin[0] +1); - sizeROI.SetElement(1, m_RegionOfInterestMax[1] - m_RegionOfInterestMin[1] +1); - sizeROI.SetElement(2, m_RegionOfInterestMax[2] - m_RegionOfInterestMin[2] +1); - regionOfInterest.SetSize(sizeROI); - typename TOutputImage::Pointer image = this->GetOutput(0); - IteratorType regionOfInterestIterator(image, regionOfInterest); - - for(regionOfInterestIterator.GoToBegin(); !regionOfInterestIterator.IsAtEnd(); ++regionOfInterestIterator) - { - numberSummand = 1; - value = regionOfInterestIterator.Get(); - m_ValueAtMidpoint = value; - index = regionOfInterestIterator.GetIndex(); - image->TransformIndexToPhysicalPoint(regionOfInterestIterator.GetIndex(), globalCoordinate); + */ + template< typename TOutputImage > + void + MultiGaussianImageSource< TOutputImage > + ::CalculateMidpointAndMeanValue() + { + itkDebugMacro(<< "Generating a image of scalars "); + + double numberSummand = 0.0; + int fiStepNumber, psiStepNumber; + double x, y, z, temp, abst; + MultiGaussianImageSource< TOutputImage >::OutputImagePixelType value, meanValueTemp, valueAdd; + PointType globalCoordinate; + IndexType index; + double riStep, fijStep, psikStep, ri, fij, psik, rNew; + double dist = itk::Math::pi * m_Radius / static_cast< double >(2 * m_RadiusStepNumber); + m_MeanValue = 0.0; + riStep = m_Radius / static_cast< double >( m_RadiusStepNumber ); + OutputImageRegionType regionOfInterest; + IndexType indexR; + indexR.SetElement(0, m_RegionOfInterestMin[0]); + indexR.SetElement(1, m_RegionOfInterestMin[1]); + indexR.SetElement(2, m_RegionOfInterestMin[2]); + regionOfInterest.SetIndex(indexR); + SizeType sizeROI; + sizeROI.SetElement(0, m_RegionOfInterestMax[0] - m_RegionOfInterestMin[0] +1); + sizeROI.SetElement(1, m_RegionOfInterestMax[1] - m_RegionOfInterestMin[1] +1); + sizeROI.SetElement(2, m_RegionOfInterestMax[2] - m_RegionOfInterestMin[2] +1); + regionOfInterest.SetSize(sizeROI); + typename TOutputImage::Pointer image = this->GetOutput(0); + IteratorType regionOfInterestIterator(image, regionOfInterest); + + for(regionOfInterestIterator.GoToBegin(); !regionOfInterestIterator.IsAtEnd(); ++regionOfInterestIterator) + { + numberSummand = 1.0; + value = regionOfInterestIterator.Get(); + m_ValueAtMidpoint = value; + index = regionOfInterestIterator.GetIndex(); + + + image->TransformIndexToPhysicalPoint(index, globalCoordinate); + ri = riStep; + + //for(unsigned int i = 0; i < m_RadiusStepNumber; ++i) + //{ + // angleStepNumberOverTwo = static_cast( itk::Math::pi * ri / dist); + // //fij = 0.0; + // fijStep = itk::Math::pi / static_cast< double >( angleStepNumberOverTwo ); + // fij = fijStep; + // std::cout << "i = " << i << " r = " << ri << "\n" << std::endl; + // for(unsigned int j = 0; j < angleStepNumberOverTwo -1 ; ++j) // from 0 to pi j<= angleStep... + // { + // z = ri * cos(fij); + // psikStep = 2.0 * itk::Math::pi / (2.0 * static_cast< double >( angleStepNumberOverTwo ) ); + // psik = -itk::Math::pi + psikStep; + // temp = ri * sin(fij); + // std::cout << " j = " << j << " fi = " << fij << "\n" << std::endl; + // for(unsigned int k = 0; k < 2 * angleStepNumberOverTwo; ++k) // from -pi to pi + // { + // std::cout << " k = " << k << " psi = " << psik << "\n" << std::endl; + // x = temp * cos(psik); + // y = temp * sin(psik); + // numberSummand += 1.0; + // // std::cout << "x, y, z: " << x <<" " << y << " " << z << "\n" < xMax ) xMax = tempX; + // //if ( tempY > yMax ) yMax = tempY; + // //if ( tempZ > zMax ) zMax = tempZ; + // + // } + // fij = fij + fijStep; + // } + // ri = ri + riStep; + //} + // second implementation ... corrected! + // for(unsigned int i = 0; i < m_RadiusStepNumber; ++i) + //{ + // angleStepNumberOverTwo = static_cast( itk::Math::pi * ri / dist); + // //fij = 0.0; + // fijStep = itk::Math::pi / static_cast< double >( angleStepNumberOverTwo ); + // fij = fijStep; + // // std::cout << "i = " << i << " r = " << ri << "\n" << std::endl; + // for(unsigned int j = 0; j < angleStepNumberOverTwo -1 ; ++j) // from 0 to pi j<= angleStep... + // { + // z = ri * cos(fij); + // /* double l = 2.0 * ri / static_cast< double >( angleStepNumberOverTwo - 1); + // double tempVar = (j + 1) * l; + // abst = (tempVar < ri ) ? tempVar : (2.0 * ri - tempVar); + // double testVar = abst * (ri - abst); + // rNew = std::sqrt(abst * (ri - abst)); */ + // if(fij < itk::Math::pi / 2.0) + // { + // rNew = sin(fij) * ri; + // } + // else + // { + // rNew = sin(itk::Math::pi - fij) * ri; + // } + // psiStepNumber = static_cast( 2.0 * itk::Math::pi * rNew / dist); + // psikStep = 2.0 * itk::Math::pi / psiStepNumber ; + // psik = -itk::Math::pi + psikStep; + // temp = ri * sin(fij); + // // std::cout << " j = " << j << " fi = " << fij << "\n" << std::endl; + // for(unsigned int k = 0; k < psiStepNumber; ++k) // from -pi to pi + // { + // // std::cout << " k = " << k << " psi = " << psik << "\n" << std::endl; + // x = temp * cos(psik); + // y = temp * sin(psik); + // numberSummand += 1.0; + // valueAdd = MultiGaussianFunctionValueAtPoint(x + globalCoordinate[0], y + globalCoordinate[1], z + globalCoordinate[2]); + // value = value + valueAdd; + // psik = psik + psikStep; + // } + // fij = fij + fijStep; + // } + // ri = ri + riStep; + //} + + + for(unsigned int i = 0; i < m_RadiusStepNumber; ++i) + { + fiStepNumber = static_cast( (itk::Math::pi * ri / dist)); + //fij = 0.0; + fijStep = itk::Math::pi / (static_cast< double >( fiStepNumber ) * 2.0); + fij = fijStep; + // std::cout << "i = " << i << " r = " << ri << "\n" << std::endl; + for(unsigned int j = 0; j < fiStepNumber -2 ; ++j) // from 0 to pi/2 j<= angleStep... + { + z = ri * cos(fij); + if(fij < itk::Math::pi / 2.0) + { + rNew = sin(fij) * ri; + } + else + { + rNew = sin(itk::Math::pi - fij) * ri; + } + + psiStepNumber = static_cast(( 2.0 * itk::Math::pi * rNew / dist) / 4.0); + psikStep = itk::Math::pi / (static_cast< double >( psiStepNumber ) * 2.0); + psik = psikStep; + temp = ri * sin(fij); + // std::cout << " j = " << j << " fi = " << fij << "\n" << std::endl; + if(psiStepNumber > 2){ + for(unsigned int k = 0; k < psiStepNumber -2; ++k) // from 0 to pi/2 + { + // std::cout << " k = " << k << " psi = " << psik << "\n" << std::endl; + x = temp * cos(psik); + y = temp * sin(psik); + numberSummand += 8.0; + + + valueAdd = MultiGaussianFunctionValueAtPoint(x + globalCoordinate[0], y + globalCoordinate[1], z + globalCoordinate[2]); + value = value + valueAdd; + valueAdd = MultiGaussianFunctionValueAtPoint( - x + globalCoordinate[0], y + globalCoordinate[1], z + globalCoordinate[2]); + value = value + valueAdd; + valueAdd = MultiGaussianFunctionValueAtPoint( - x + globalCoordinate[0], - y + globalCoordinate[1], - z + globalCoordinate[2]); + value = value + valueAdd; + valueAdd = MultiGaussianFunctionValueAtPoint(x + globalCoordinate[0], - y + globalCoordinate[1], z + globalCoordinate[2]); + value = value + valueAdd; + valueAdd = MultiGaussianFunctionValueAtPoint(x + globalCoordinate[0], - y + globalCoordinate[1], - z + globalCoordinate[2]); + value = value + valueAdd; + valueAdd = MultiGaussianFunctionValueAtPoint(x + globalCoordinate[0], y + globalCoordinate[1], - z + globalCoordinate[2]); + value = value + valueAdd; + valueAdd = MultiGaussianFunctionValueAtPoint(- x + globalCoordinate[0], y + globalCoordinate[1], - z + globalCoordinate[2]); + value = value + valueAdd; + valueAdd = MultiGaussianFunctionValueAtPoint(- x + globalCoordinate[0], - y + globalCoordinate[1], z + globalCoordinate[2]); + value = value + valueAdd; + + + + + psik = psik + psikStep; + + } + } + fij = fij + fijStep; + } + ri = ri + riStep; + } + + + meanValueTemp = value / numberSummand; + if(meanValueTemp > m_MeanValue) + { + m_GlobalCoordinate = globalCoordinate; + m_MeanValue = meanValueTemp; + m_SphereMidpoint = index; + } + } + } + + //------------------------------------------------------------------------------------------------------------------------- + template< typename TOutputImage > + void + MultiGaussianImageSource< TOutputImage > + ::WriteXMLToTestTheCuboidInsideTheSphere() + { + std::stringstream ss; + + + int numberSummand = 1.0; + + //write an .mps test file + itk::DOMNodeXMLWriter::Pointer xmlWriter; + typedef itk::DOMNode::Pointer DOMNodeType; + DOMNodeType domXML, domPointSetFile, domFileVersion, domPointSet, domPoint, domId, domX, domY, domZ; + + xmlWriter = itk::DOMNodeXMLWriter::New(); + domXML = itk::DOMNode::New(); + domXML->SetName("?xml"); + domPointSetFile = itk::DOMNode::New(); + domPointSetFile->SetName("point_set_file"); + //domFileVersion = itk::DOMNode::New(); + //domFileVersion->SetName("file_version"); + domPointSet = itk::DOMNode::New(); + domPointSet->SetName("point_set"); + + ss.str(""); + ss << 1.0; + domXML->SetAttribute("version", ss.str()); + domXML->AddChildAtBegin(domPointSetFile); + //domPointSetFile -> AddChildAtBegin(domFileVersion); + domPointSetFile -> AddChildAtBegin(domPointSet); + + + int cap = m_Midpoints.Size(); + for(unsigned int iter = 0 ; iter < cap; ++iter) + { + numberSummand += 1.0; + + domPoint = itk::DOMNode::New(); + domPoint->SetName("point"); + domX = itk::DOMNode::New(); + domX->SetName("x"); + domY = itk::DOMNode::New(); + domY->SetName("y"); + domZ = itk::DOMNode::New(); + domZ->SetName("z"); + domId = itk::DOMNode::New(); + domId->SetName("id"); + + + ss.str(""); + ss << numberSummand; + domId->AddTextChildAtBegin(ss.str()); + domPoint -> AddChildAtEnd(domId); + + double scaleFactor = 10.0; + PointType point = m_Midpoints.GetElement(numberSummand); + ss.str(""); + ss << point[0] * scaleFactor; + domX->AddTextChildAtBegin(ss.str()); + domPoint -> AddChildAtEnd(domX); + ss.str(""); + ss << point[1] * scaleFactor; + domY->AddTextChildAtBegin(ss.str()); + domPoint -> AddChildAtEnd(domY); + ss.str(""); + ss << point[2] * scaleFactor; + domZ->AddTextChildAtBegin(ss.str()); + domPoint -> AddChildAtEnd(domZ); + domPointSet -> AddChildAtEnd(domPoint); + + } + + + // .mps (Data) + ss.str(""); + ss << "C:/temp/CuboidsInTheSphere.mps"; + std::string name = ss.str(); + char * fileNamePointer = (char*) name.c_str(); + xmlWriter->SetFileName( fileNamePointer); + xmlWriter->SetInput( domXML ); + xmlWriter->Update(); + + + + + } + + //--------------------------------------------------------------------------------------------------------------------- + + template< typename TOutputImage > + void + MultiGaussianImageSource< TOutputImage > + ::WriteXMLToTestTheCuboid() + { + // to create a point set + + + std::stringstream ss; + int fiStepNumber, psiStepNumber; + double x, y, z, temp; + MultiGaussianImageSource< TOutputImage >::OutputImagePixelType meanValueTemp; + PointType globalCoordinate; + double riStep, fijStep, psikStep, ri, fij, psik, rNew; + + + int angleStepNumber = 4; + double dist = 2.0 * itk::Math::pi * m_Radius / static_cast< double >(angleStepNumber * m_RadiusStepNumber); + m_MeanValue = 0.0; + riStep = m_Radius / static_cast< double >( m_RadiusStepNumber); + int numberSummand = 1.0; + + //write an .mps test file + itk::DOMNodeXMLWriter::Pointer xmlWriter; + typedef itk::DOMNode::Pointer DOMNodeType; + DOMNodeType domXML, domPointSetFile, domFileVersion, domPointSet, domPoint, domId, domX, domY, domZ; + + xmlWriter = itk::DOMNodeXMLWriter::New(); + domXML = itk::DOMNode::New(); + domXML->SetName("?xml"); + domPointSetFile = itk::DOMNode::New(); + domPointSetFile->SetName("point_set_file"); + //domFileVersion = itk::DOMNode::New(); + //domFileVersion->SetName("file_version"); + domPointSet = itk::DOMNode::New(); + domPointSet->SetName("point_set"); + + ss.str(""); + ss << 1.0; + domXML->SetAttribute("version", ss.str()); + domXML->AddChildAtBegin(domPointSetFile); + //domPointSetFile -> AddChildAtBegin(domFileVersion); + domPointSetFile -> AddChildAtBegin(domPointSet); + + + int cap = m_XCoordToTest.size(); + for(unsigned int iter = 0 ; iter < cap; ++iter) + { + numberSummand += 1.0; + + domPoint = itk::DOMNode::New(); + domPoint->SetName("point"); + domX = itk::DOMNode::New(); + domX->SetName("x"); + domY = itk::DOMNode::New(); + domY->SetName("y"); + domZ = itk::DOMNode::New(); + domZ->SetName("z"); + domId = itk::DOMNode::New(); + domId->SetName("id"); + + + ss.str(""); + ss << numberSummand; + domId->AddTextChildAtBegin(ss.str()); + domPoint -> AddChildAtEnd(domId); + + double scaleFactor = 10.0; + + ss.str(""); + ss << m_XCoordToTest[iter] * scaleFactor; + domX->AddTextChildAtBegin(ss.str()); + domPoint -> AddChildAtEnd(domX); + ss.str(""); + ss << m_YCoordToTest[iter] * scaleFactor; + domY->AddTextChildAtBegin(ss.str()); + domPoint -> AddChildAtEnd(domY); + ss.str(""); + ss << m_ZCoordToTest[iter] * scaleFactor; + domZ->AddTextChildAtBegin(ss.str()); + domPoint -> AddChildAtEnd(domZ); + domPointSet -> AddChildAtEnd(domPoint); + + + + } + + + // .mps (Data) + ss.str(""); + ss << "C:/temp/CuboidMidpointsSet.mps"; + std::string name = ss.str(); + char * fileNamePointer = (char*) name.c_str(); + xmlWriter->SetFileName( fileNamePointer); + xmlWriter->SetInput( domXML ); + xmlWriter->Update(); + + + } + //------------------------------------------------------------------------------------------------------------------ + template< typename TOutputImage > + void + MultiGaussianImageSource< TOutputImage > + ::WriteXMLToTest() + { + // to create a point set + + double numberSummand = 0.0; + std::stringstream ss; + int fiStepNumber, psiStepNumber; + double x, y, z, temp; + MultiGaussianImageSource< TOutputImage >::OutputImagePixelType meanValueTemp; + PointType globalCoordinate; + double riStep, fijStep, psikStep, ri, fij, psik, rNew; + + + int angleStepNumber = 4; + double dist = 2.0 * itk::Math::pi * m_Radius / static_cast< double >(angleStepNumber * m_RadiusStepNumber); + m_MeanValue = 0.0; + riStep = m_Radius / static_cast< double >( m_RadiusStepNumber); + numberSummand = 1.0; + + //write an .mps test file + itk::DOMNodeXMLWriter::Pointer xmlWriter; + typedef itk::DOMNode::Pointer DOMNodeType; + DOMNodeType domXML, domPointSetFile, domFileVersion, domPointSet, domPoint, domId, domX, domY, domZ; + + xmlWriter = itk::DOMNodeXMLWriter::New(); + domXML = itk::DOMNode::New(); + domXML->SetName("?xml"); + domPointSetFile = itk::DOMNode::New(); + domPointSetFile->SetName("point_set_file"); + //domFileVersion = itk::DOMNode::New(); + //domFileVersion->SetName("file_version"); + domPointSet = itk::DOMNode::New(); + domPointSet->SetName("point_set"); + + ss.str(""); + ss << 1.0; + domXML->SetAttribute("version", ss.str()); + domXML->AddChildAtBegin(domPointSetFile); + //domPointSetFile -> AddChildAtBegin(domFileVersion); + domPointSetFile -> AddChildAtBegin(domPointSet); + + + typename TOutputImage::Pointer image = this->GetOutput(0); + image->TransformIndexToPhysicalPoint(m_SphereMidpoint, globalCoordinate); ri = riStep; for(unsigned int i = 0; i < m_RadiusStepNumber; ++i) + { + fiStepNumber = static_cast( (itk::Math::pi * ri / dist) / 2.0); + //fij = 0.0; + fijStep = itk::Math::pi / (static_cast< double >( fiStepNumber ) * 2.0); + fij = fijStep; + // std::cout << "i = " << i << " r = " << ri << "\n" << std::endl; + for(unsigned int j = 0; j < fiStepNumber - 2 ; ++j) // from 0 to pi/2 j<= angleStep... + { + z = ri * cos(fij); + if(fij < itk::Math::pi / 2.0) + { + rNew = sin(fij) * ri; + } + else + { + rNew = sin(itk::Math::pi - fij) * ri; + } + + psiStepNumber = static_cast(( 2.0 * itk::Math::pi * rNew / dist) / 4.0); + psikStep = itk::Math::pi / (static_cast< double >( psiStepNumber ) * 2.0); + psik = psikStep; + temp = ri * sin(fij); + + for(unsigned int k = 0; k < psiStepNumber - 2; ++k) // from 0 to pi/2 + { + x = temp * cos(psik); + y = temp * sin(psik); + numberSummand += 1.0; + + domPoint = itk::DOMNode::New(); + domPoint->SetName("point"); + domX = itk::DOMNode::New(); + domX->SetName("x"); + domY = itk::DOMNode::New(); + domY->SetName("y"); + domZ = itk::DOMNode::New(); + domZ->SetName("z"); + domId = itk::DOMNode::New(); + domId->SetName("id"); + + + ss.str(""); + ss << numberSummand; + domId->AddTextChildAtBegin(ss.str()); + domPoint -> AddChildAtEnd(domId); + + double scaleFactor = 40.0; + + ss.str(""); + ss << (x + globalCoordinate[0]) * scaleFactor; + domX->AddTextChildAtBegin(ss.str()); + domPoint -> AddChildAtEnd(domX); + ss.str(""); + ss << (y + globalCoordinate[1]) * scaleFactor; + domY->AddTextChildAtBegin(ss.str()); + domPoint -> AddChildAtEnd(domY); + ss.str(""); + ss << (z + globalCoordinate[2]) * scaleFactor; + domZ->AddTextChildAtBegin(ss.str()); + domPoint -> AddChildAtEnd(domZ); + domPointSet -> AddChildAtEnd(domPoint); + + + psik = psik + psikStep; + + } + fij = fij + fijStep; + } + ri = ri + riStep; + } + + + // .mps (Data) + ss.str(""); + ss << "C:/temp/SpherePointSet.mps"; + std::string name = ss.str(); + char * fileNamePointer = (char*) name.c_str(); + xmlWriter->SetFileName( fileNamePointer); + xmlWriter->SetInput( domXML ); + xmlWriter->Update(); + + + } + + //------------------------------------------------------------------------------------------------------------------------------------------- + template< typename TOutputImage > + void + MultiGaussianImageSource< TOutputImage > + ::OptimizeMeanValue() + { + int radiusStepNumberOptimized = m_RadiusStepNumber * 4; + int numberSummand = 1, angleStepNumberOverTwo; + double x, y, z, temp; + double riStep, fijStep, psikStep, ri, fij, psik; + double dist = itk::Math::pi * m_Radius / (2 * radiusStepNumberOptimized); + MultiGaussianImageSource< TOutputImage >::OutputImagePixelType valueAdd, value, + m_MeanValue = 0; + riStep = m_Radius / radiusStepNumberOptimized; + ri = riStep; + value = m_ValueAtMidpoint; + for(unsigned int i = 0; i < radiusStepNumberOptimized; ++i) { angleStepNumberOverTwo = static_cast( itk::Math::pi * ri / dist); fij = 0.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 + globalCoordinate[0], y + globalCoordinate[1], z + globalCoordinate[2]); + valueAdd = MultiGaussianFunctionValueAtPoint(x + m_GlobalCoordinate[0] , y + m_GlobalCoordinate[1], z + m_GlobalCoordinate[2]); value = value + valueAdd; psik = psik + psikStep; } fij = fij + fijStep; } ri = ri + riStep; } - meanValueTemp = value / numberSummand; - if(meanValueTemp > m_MeanValue) - { - m_GlobalCoordinate = globalCoordinate; - m_MeanValue = meanValueTemp; - m_SphereMidpoint = index; - } + m_MeanValue = value / numberSummand; } -} -//------------------------------------------------------------------------------------------------------------------------------------------- -template< typename TOutputImage > -void -MultiGaussianImageSource< TOutputImage > -::OptimizeMeanValue() -{ - int radiusStepNumberOptimized = m_RadiusStepNumber * 4; - int numberSummand = 1, angleStepNumberOverTwo; - double x, y, z, temp; - double riStep, fijStep, psikStep, ri, fij, psik; - double dist = itk::Math::pi * m_Radius / (2 * radiusStepNumberOptimized); - MultiGaussianImageSource< TOutputImage >::OutputImagePixelType valueAdd, value, - m_MeanValue = 0; - riStep = m_Radius / radiusStepNumberOptimized; - ri = riStep; - value = m_ValueAtMidpoint; - for(unsigned int i = 0; i < radiusStepNumberOptimized; ++i) - { - angleStepNumberOverTwo = static_cast( itk::Math::pi * ri / dist); - fij = 0.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 + //------------------------------------------------------------------------------------------------------------------------------------------- + template< typename TOutputImage > + void + MultiGaussianImageSource< TOutputImage > + ::CalculateMaxAndMinInSphere() + { + IndexType index; + MultiGaussianImageSource< TOutputImage >::OutputImagePixelType value; + m_MaxValueInSphere = std::numeric_limits::min(); + m_MinValueInSphere = std::numeric_limits::max(); + int radInt; + OutputImageRegionType cuboidRegion; + IndexType indexR; + SizeType sizeR; + int indexRegion, sizeRegion; + for( unsigned int i = 0; i < 3; ++i ) + { + radInt = static_cast(m_Radius/m_Spacing[i]); + indexRegion = m_SphereMidpoint[i] - radInt; + if( m_Origin[i] > indexRegion ) { - x = temp * cos(psik); - y = temp * sin(psik); - numberSummand++; - valueAdd = MultiGaussianFunctionValueAtPoint(x + m_GlobalCoordinate[0] , y + m_GlobalCoordinate[1], z + m_GlobalCoordinate[2]); - value = value + valueAdd; - psik = psik + psikStep; + indexR.SetElement(i, m_Origin[i] ); } - fij = fij + fijStep; - } - ri = ri + riStep; - } - m_MeanValue = value / numberSummand; -} - + else + { + indexR.SetElement(i, indexRegion ); + } + sizeRegion = 2 *radInt + 1; -//------------------------------------------------------------------------------------------------------------------------------------------- -template< typename TOutputImage > -void -MultiGaussianImageSource< TOutputImage > -::CalculateMaxAndMinInSphere() -{ - IndexType index; - MultiGaussianImageSource< TOutputImage >::OutputImagePixelType value; - m_MaxValueInSphere = std::numeric_limits::min(); - m_MinValueInSphere = std::numeric_limits::max(); - int radInt0 = static_cast(m_Radius/m_Spacing[0]); - int radInt1 = static_cast(m_Radius/m_Spacing[1]); - int radInt2 = static_cast(m_Radius/m_Spacing[2]); - OutputImageRegionType cuboidRegion; - IndexType indexR; - indexR.SetElement(0, m_SphereMidpoint[0] - radInt0); - indexR.SetElement(1, m_SphereMidpoint[1] - radInt1); - indexR.SetElement(2, m_SphereMidpoint[2] - radInt2); - cuboidRegion.SetIndex(indexR); - SizeType size; - size.SetElement(0, 2 *radInt0 + 1); - size.SetElement(1, 2 *radInt1 + 1); - size.SetElement(2, 2* radInt2 + 1); - cuboidRegion.SetSize(size); - typename TOutputImage::Pointer image = this->GetOutput(0); - IteratorType cuboidRegionOfInterestIterator(image, cuboidRegion); - PointType globalCoordinate; - for(cuboidRegionOfInterestIterator.GoToBegin(); !cuboidRegionOfInterestIterator.IsAtEnd(); ++cuboidRegionOfInterestIterator) - { - index = cuboidRegionOfInterestIterator.GetIndex(); - image->TransformIndexToPhysicalPoint(cuboidRegionOfInterestIterator.GetIndex(), globalCoordinate); - if( m_GlobalCoordinate.EuclideanDistanceTo(globalCoordinate) <= m_Radius ) - { - value = cuboidRegionOfInterestIterator.Get(); - if(m_MaxValueInSphere < value) + if(indexR[i] + sizeRegion > m_Size[i]) { - m_MaxValueInSphere = value; - m_MaxValueIndexInSphere = index; + sizeR.SetElement(i, m_Size[i] ); } - if(m_MinValueInSphere > value) + else { - m_MinValueInSphere = value; - m_MinValueIndexInSphere = index; + sizeR.SetElement(i, sizeRegion ); + } + + } + cuboidRegion.SetIndex(indexR); + cuboidRegion.SetSize(sizeR); + typename TOutputImage::Pointer image = this->GetOutput(0); + IteratorType cuboidRegionOfInterestIterator(image, cuboidRegion); + PointType globalCoordinate; + for(cuboidRegionOfInterestIterator.GoToBegin(); !cuboidRegionOfInterestIterator.IsAtEnd(); ++cuboidRegionOfInterestIterator) + { + index = cuboidRegionOfInterestIterator.GetIndex(); + if(IsInImage(index)) + { + image->TransformIndexToPhysicalPoint(cuboidRegionOfInterestIterator.GetIndex(), globalCoordinate); + if( m_GlobalCoordinate.EuclideanDistanceTo(globalCoordinate) <= m_Radius ) + { + value = cuboidRegionOfInterestIterator.Get(); + if(m_MaxValueInSphere < value) + { + m_MaxValueInSphere = value; + m_MaxValueIndexInSphere = index; + } + if(m_MinValueInSphere > value) + { + m_MinValueInSphere = value; + m_MinValueIndexInSphere = index; + } + } } } } -} + //-------------------------------------------------------------------------------------------------------------------------- + + template< typename TOutputImage > + bool + MultiGaussianImageSource< TOutputImage > + ::IsInImage(IndexType index) + { + bool isInImage = 0; + if(index[0] >= m_Origin[0] && index[0] <= m_Size[0] && index[1] >= m_Origin[1] && index[1] <= m_Size[1] && index[2] >= m_Origin[2] && index[2] <= m_Size[2]) + { + isInImage = 1; + } + return isInImage; + } + + + //----------------------------------------------------------------------------------------------------------------------- + + + template< typename TOutputImage > + void + MultiGaussianImageSource< TOutputImage > + ::SetNormalDistributionValues() + { + + double temp[410] = { 0.5 , 0.50399 , 0.50798, 0.51197, 0.51595, 0.51994, 0.52392, 0.5279, 0.53188, 0.53586, 0.53983, 0.5438, 0.54776, 0.55172, 0.55567, 0.55962, 0.56356, 0.56749, 0.57142, 0.57535, 0.57926, 0.58317, 0.58706, 0.59095, 0.59483 , 0.59871, 0.60257, 0.60642, 0.61026, 0.61409, 0.61791, 0.62172, 0.62552, 0.6293, 0.63307, 0.63683, 0.64058, 0.64431, 0.64803, 0.65173, 0.65542, 0.6591, 0.66276, 0.6664, 0.67003, 0.67364, 0.67724, 0.68082, 0.68439, 0.68793, 0.69146, 0.69497, 0.69847, 0.70194, 0.7054, 0.70884, 0.71226, 0.71566, 0.71904, 0.7224, 0.72575, 0.72907, 0.73237, 0.73565, 0.73891, 0.74215, 0.74537, 0.74857, 0.75175, 0.7549, 0.75804, 0.76115, 0.76424, 0.7673, 0.77035, 0.77337, 0.77637, 0.77935, 0.7823, 0.78524, 0.78814, 0.79103, 0.79389, 0.79673, 0.79955, 0.80234, 0.80511, 0.80785, 0.81057, 0.81327, 0.81594, 0.81859, 0.82121, 0.82381, 0.82639, 0.82894, 0.83147, 0.83398, 0.83646, 0.83891, 0.84134, 0.84375, 0.84614, 0.84849, 0.85083, 0.85314, 0.85543, 0.85769, 0.85993, 0.86214, 0.86433, 0.8665, 0.86864, 0.87076, 0.87286, 0.87493, 0.87698, 0.879, 0.881, 0.88298, 0.88493, 0.88686, 0.88877, 0.89065, 0.89251, 0.89435, 0.89617, 0.89796, 0.89973, 0.90147, 0.9032, 0.9049, 0.90658, 0.90824, 0.90988, 0.91149, 0.91309, 0.91466, 0.91621, 0.91774, 0.91924, 0.92073, 0.9222, 0.92364, 0.92507, 0.92647, 0.92785, 0.92922, 0.93056, 0.93189, 0.93319, 0.93448, 0.93574, 0.93699, 0.93822, 0.93943, 0.94062, 0.94179, 0.94295, 0.94408, 0.9452, 0.9463, 0.94738, 0.94845, 0.9495, 0.95053, 0.95154, 0.95254, 0.95352, 0.95449, 0.95543, 0.95637, 0.95728, 0.95818, 0.95907, 0.95994, 0.9608, 0.96164, 0.96246, 0.96327, 0.96407, 0.96485, 0.96562, 0.96638, 0.96712, 0.96784, 0.96856, 0.96926, 0.96995, 0.97062, 0.97128, 0.97193, 0.97257, 0.9732, 0.97381, 0.97441, 0.975, 0.97558, 0.97615, 0.9767, 0.97725, 0.97778, 0.97831, 0.97882, 0.97932, 0.97982, 0.9803, 0.98077, 0.98124, 0.98169, 0.98214, 0.98257, 0.983, 0.98341, 0.98382, 0.98422, 0.98461, 0.985, 0.98537, 0.98574, 0.9861, 0.98645, 0.98679, 0.98713, 0.98745, 0.98778, 0.98809, 0.9884, 0.9887, 0.98899, 0.98928, 0.98956, 0.98983, 0.9901, 0.99036, 0.99061, 0.99086, 0.99111, 0.99134, 0.99158, 0.9918, 0.99202, 0.99224, 0.99245, 0.99266, 0.99286, 0.99305, 0.99324, 0.99343, 0.99361, 0.99379, 0.99396, 0.99413, 0.9943, 0.99446, 0.99461, 0.99477, 0.99492, 0.99506, 0.9952, 0.99534, 0.99547, 0.9956, 0.99573, 0.99585, 0.99598, 0.99609, 0.99621, 0.99632, 0.99643, 0.99653, 0.99664, 0.99674, 0.99683, 0.99693, 0.99702, 0.99711, 0.9972, 0.99728, 0.99736, 0.99744, 0.99752, 0.9976, 0.99767, 0.99774, 0.99781, 0.99788, 0.99795, 0.99801, 0.99807, 0.99813, 0.99819, 0.99825, 0.99831, 0.99836, 0.99841, 0.99846, 0.99851, 0.99856, 0.99861, 0.99865, 0.99869, 0.99874, 0.99878, 0.99882, 0.99886, 0.99889, 0.99893, 0.99896, 0.999, 0.99903, 0.99906, 0.9991, 0.99913, 0.99916, 0.99918, 0.99921, 0.99924, 0.99926, 0.99929, 0.99931, 0.99934, 0.99936, 0.99938, 0.9994, 0.99942, 0.99944, 0.99946, 0.99948, 0.9995, 0.99952, 0.99953, 0.99955, 0.99957, 0.99958, 0.9996, 0.99961, 0.99962, 0.99964, 0.99965, 0.99966, 0.99968, 0.99969, 0.9997, 0.99971, 0.99972, 0.99973, 0.99974, 0.99975, 0.99976, 0.99977, 0.99978, 0.99978, 0.99979, 0.9998, 0.99981, 0.99981, 0.99982, 0.99983, 0.99983, 0.99984, 0.99985, 0.99985, 0.99986, 0.99986, 0.99987, 0.99987, 0.99988, 0.99988, 0.99989, 0.99989, 0.9999, 0.9999, 0.9999, 0.99991, 0.99991, 0.99992, 0.99992, 0.99992, 0.99992, 0.99993, 0.99993, 0.99993, 0.99994, 0.99994, 0.99994, 0.99994, 0.99995, 0.99995, 0.99995, 0.99995, 0.99995, 0.99996, 0.99996, 0.99996, 0.99996, 0.99996, 0.99996, 0.99997, 0.99997, 0.99997, 0.99997, 0.99997, 0.99997, 0.99997, 0.99997, 0.99998, 0.99998, 0.99998, 0.99998 }; + + for(int i=0; i < 410; i++) + { + m_NormalDistValues[i] = temp[i]; + } + + } + + } // end namespace itk #endif + +