diff --git a/Modules/ImageStatistics/Testing/itkMultiGaussianImageSource.h b/Modules/ImageStatistics/Testing/itkMultiGaussianImageSource.h index 64d24106c0..98a236789c 100644 --- a/Modules/ImageStatistics/Testing/itkMultiGaussianImageSource.h +++ b/Modules/ImageStatistics/Testing/itkMultiGaussianImageSource.h @@ -1,184 +1,200 @@ /*========================================================================= * * Copyright Insight Software Consortium * * Licensed under the Apache License, Version 2.0 (the "License"); * you may not use this file except in compliance with the License. * You may obtain a copy of the License at * * http://www.apache.org/licenses/LICENSE-2.0.txt * * Unless required by applicable law or agreed to in writing, software * distributed under the License is distributed on an "AS IS" BASIS, * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. * See the License for the specific language governing permissions and * limitations under the License. * *=========================================================================*/ /*========================================================================= * * Portions of this file are subject to the VTK Toolkit Version 3 copyright. * * Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen * * For complete copyright, license and disclaimer of warranty information * please refer to the NOTICE file at the top of the ITK source tree. * *=========================================================================*/ #ifndef __itkMultiGaussianImageSource_h #define __itkMultiGaussianImageSource_h #include "itkImageSource.h" #include "itkNumericTraits.h" #include "itkImageRegionIteratorWithIndex.h" #include "itkImageFileWriter.h" namespace itk { /** \class MultiGaussianImageSource - * \brief Generate an n-dimensional image TODO + * \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 MultiThreaded + * \ingroup DataSources * \ingroup ITKTestKernel * * \wiki - * \wikiexample{SimpleOperations/RandomImageSource,Produce an image of noise} + * \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; - /** Run-time type information (and related methods). */ - //itkTypeMacro(RandomImageSource, ImageSource); - /** 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 double RadiusType; - typedef std::vector VectorType; - typedef ImageRegionIteratorWithIndex IteratorType; - typedef typename TOutputImage::Pointer ImageType; + typedef typename TOutputImage::SizeType SizeType; + typedef typename TOutputImage::IndexType IndexType; + typedef typename TOutputImage::SpacingType SpacingType; + typedef typename TOutputImage::PointType PointType; + typedef typename SizeType::SizeValueType SizeValueType; + typedef SizeValueType SizeValueArrayType[TOutputImage::ImageDimension]; + typedef typename TOutputImage::SpacingValueType SpacingValueType; + typedef SpacingValueType SpacingValueArrayType[TOutputImage::ImageDimension]; + typedef typename TOutputImage::PointValueType PointValueType; + typedef PointValueType PointValueArrayType[TOutputImage::ImageDimension]; + /** Typedef to describe the sphere radius type. */ + typedef double RadiusType; + /** Typedef to describe the standard vector type. */ + typedef std::vector VectorType; + /** Typedef to describe the itk vector type. */ + typedef itk::Vector ItkVectorType; + /** Typedef to describe the ImageRegionIteratorWithIndex type. */ + typedef ImageRegionIteratorWithIndex IteratorType; + /** Typedef to describe the Poiner type at the output image. */ + typedef typename TOutputImage::Pointer ImageType; /** Set/Get size of the output image */ itkSetMacro(Size, SizeType); virtual void SetSize(SizeValueArrayType sizeArray); virtual const SizeValueType * GetSize() const; /** Set/Get spacing of the output image */ itkSetMacro(Spacing, SpacingType); virtual void SetSpacing(SpacingValueArrayType spacingArray); - virtual const SpacingValueType * GetSpacing() const; /** Set/Get origin of the output image */ itkSetMacro(Origin, PointType); virtual void SetOrigin(PointValueArrayType originArray); - virtual const PointValueType * GetOrigin() const; + /** Get the number of gaussian functions in the output image */ virtual unsigned int GetNumberOfGaussians() const; + + virtual void SetNumberOfGausssians( unsigned int ); + /** Set/Get the radius of the sphere */ virtual const RadiusType GetRadius() const; virtual void SetRadius( RadiusType radius ); + /** Set/Get the number of steps to traverse the radius of the sphere */ virtual const int GetRadiusStepNumber() const; virtual void SetRadiusStepNumber( unsigned int stepNumber ); + /** Get the maximal mean value in a sphere over all possible spheres with midpoint in the image */ virtual const double GetMaxMeanValue() const; - virtual const IndexType GetSphereMidpoint() 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); - virtual void SetNumberOfGausssians( unsigned int ); + /** Calculates and set the index of the midpoint of the sphere with the maximal mean value as well as the mean value*/ virtual void CalculateMidpointAndMeanValue(); /** Set the minimum possible pixel value. By default, it is * NumericTraits::min(). */ itkSetClampMacro( Min, OutputImagePixelType, NumericTraits< OutputImagePixelType >::NonpositiveMin(), NumericTraits< OutputImagePixelType >::max() ); /** Get the minimum possible pixel value. */ itkGetConstMacro(Min, OutputImagePixelType); /** Set the maximum possible pixel value. By default, it is * NumericTraits::max(). */ itkSetClampMacro( Max, OutputImagePixelType, NumericTraits< OutputImagePixelType >::NonpositiveMin(), NumericTraits< OutputImagePixelType >::max() ); /** Get the maximum possible pixel value. */ itkGetConstMacro(Max, OutputImagePixelType); protected: MultiGaussianImageSource(); ~MultiGaussianImageSource(); void PrintSelf(std::ostream & os, Indent indent) const; virtual void GenerateData(); virtual void GenerateOutputInformation(); private: MultiGaussianImageSource(const MultiGaussianImageSource &); //purposely not implemented void operator=(const MultiGaussianImageSource &); //purposely not implemented SizeType m_Size; //size of the output image SpacingType m_Spacing; //spacing PointType m_Origin; //origin unsigned int m_NumberOfGaussians; //number of Gaussians - RadiusType m_Radius; //radius of the shpere - unsigned int m_RadiusStepNumber; //number of steps to treverse the sphere radius - OutputImagePixelType m_MeanValue; //mean value in the wanted shpere - IndexType m_ShpereMidpoint; //midpoint of the wanted shpere + 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 - ImageType m_Image; 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 index 400ca63780..f39a2b7122 100644 --- a/Modules/ImageStatistics/Testing/itkMultiGaussianImageSource.hxx +++ b/Modules/ImageStatistics/Testing/itkMultiGaussianImageSource.hxx @@ -1,438 +1,446 @@ /*========================================================================= * * Copyright Insight Software Consortium * * Licensed under the Apache License, Version 2.0 (the "License"); * you may not use this file except in compliance with the License. * You may obtain a copy of the License at * * http://www.apache.org/licenses/LICENSE-2.0.txt * * Unless required by applicable law or agreed to in writing, software * distributed under the License is distributed on an "AS IS" BASIS, * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. * See the License for the specific language governing permissions and * limitations under the License. * *=========================================================================*/ /*========================================================================= * * Portions of this file are subject to the VTK Toolkit Version 3 copyright. * * Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen * * For complete copyright, license and disclaimer of warranty information * please refer to the NOTICE file at the top of the ITK source tree. * *=========================================================================*/ #ifndef __itkMultiGaussianImageSource_hxx #define __itkMultiGaussianImageSource_hxx - - #include #include #include #include "itkMultiGaussianImageSource.h" #include "itkImageRegionIterator.h" #include "itkObjectFactory.h" #include "itkProgressReporter.h" #include "stdlib.h" namespace itk { /** * */ template< class TOutputImage > MultiGaussianImageSource< TOutputImage > ::MultiGaussianImageSource() { //Initial image is 100 wide in each direction. for ( unsigned int i = 0; i < TOutputImage::GetImageDimension(); i++ ) { m_Size[i] = 100; m_Spacing[i] = 1.0; m_Origin[i] = 0.0; - m_ShpereMidpoint[i] = 0; + m_SphereMidpoint[i] = 0; } - m_NumberOfGaussians = 1; + m_NumberOfGaussians = 0; m_Radius = 1; m_RadiusStepNumber = 5; m_MeanValue = 0; m_Min = NumericTraits< OutputImagePixelType >::NonpositiveMin(); m_Max = NumericTraits< OutputImagePixelType >::max(); } template< class TOutputImage > MultiGaussianImageSource< TOutputImage > ::~MultiGaussianImageSource() {} template< class TOutputImage > void MultiGaussianImageSource< TOutputImage > ::SetSize(SizeValueArrayType sizeArray) { const unsigned int count = TOutputImage::ImageDimension; unsigned int i; for ( i = 0; i < count; i++ ) { if ( sizeArray[i] != this->m_Size[i] ) { break; } } if ( i < count ) { this->Modified(); for ( i = 0; i < count; i++ ) { this->m_Size[i] = sizeArray[i]; } } } template< class TOutputImage > const typename MultiGaussianImageSource< TOutputImage >::SizeValueType * MultiGaussianImageSource< TOutputImage > ::GetSize() const { return this->m_Size.GetSize(); } template< class TOutputImage > void MultiGaussianImageSource< TOutputImage > ::SetSpacing(SpacingValueArrayType spacingArray) { const unsigned int count = TOutputImage::ImageDimension; unsigned int i; for ( i = 0; i < count; i++ ) { if ( spacingArray[i] != this->m_Spacing[i] ) { break; } } if ( i < count ) { this->Modified(); for ( i = 0; i < count; i++ ) { this->m_Spacing[i] = spacingArray[i]; } } } template< class TOutputImage > void MultiGaussianImageSource< TOutputImage > ::SetOrigin(PointValueArrayType originArray) { const unsigned int count = TOutputImage::ImageDimension; unsigned int i; for ( i = 0; i < count; i++ ) { if ( originArray[i] != this->m_Origin[i] ) { break; } } if ( i < count ) { this->Modified(); for ( i = 0; i < count; i++ ) { this->m_Origin[i] = originArray[i]; } } } template< class TOutputImage > const typename MultiGaussianImageSource< TOutputImage >::PointValueType * MultiGaussianImageSource< TOutputImage > ::GetOrigin() const { for ( unsigned int i = 0; i < TOutputImage::ImageDimension; i++ ) { this->m_OriginArray[i] = this->m_Origin[i]; } return this->m_OriginArray; } template< class TOutputImage > const typename MultiGaussianImageSource< TOutputImage >::SpacingValueType * MultiGaussianImageSource< TOutputImage > ::GetSpacing() const { for ( unsigned int i = 0; i < TOutputImage::ImageDimension; i++ ) { this->m_SpacingArray[i] = this->m_Spacing[i]; } return this->m_SpacingArray; } /** * */ template< class TOutputImage > void MultiGaussianImageSource< TOutputImage > ::PrintSelf(std::ostream & os, Indent indent) const { Superclass::PrintSelf(os, indent); os << indent << "Max: " << static_cast< typename NumericTraits< OutputImagePixelType >::PrintType >( m_Max ) << std::endl; os << indent << "Min: " << static_cast< typename NumericTraits< OutputImagePixelType >::PrintType >( m_Min ) << std::endl; os << indent << "Origin: ["; unsigned int ii = 0; while( ii < TOutputImage::ImageDimension - 1 ) { os << m_Origin[ii] << ", "; ++ii; } os << m_Origin[ii] << "]" << std::endl; os << indent << "Spacing: ["; ii = 0; while( ii < TOutputImage::ImageDimension - 1 ) { os << m_Spacing[ii] << ", "; ++ii; } os << m_Spacing[ii] << "]" << std::endl; os << indent << "Size: ["; ii = 0; while( ii < TOutputImage::ImageDimension - 1 ) { os << m_Size[ii] << ", "; ++ii; } os << m_Size[ii] << "]" << std::endl; } template< class TOutputImage > unsigned int MultiGaussianImageSource< TOutputImage > ::GetNumberOfGaussians() const { return this->m_NumberOfGaussians; } template< class TOutputImage > const typename MultiGaussianImageSource< TOutputImage >::RadiusType MultiGaussianImageSource< TOutputImage > ::GetRadius() const { return this->m_Radius; } template< class TOutputImage > void MultiGaussianImageSource< TOutputImage > ::SetRadius( RadiusType radius ) { this->m_Radius = radius; } template< class TOutputImage > const int MultiGaussianImageSource< TOutputImage > ::GetRadiusStepNumber() const { return this->m_RadiusStepNumber; } template< class TOutputImage > void MultiGaussianImageSource< TOutputImage > ::SetRadiusStepNumber( unsigned int stepNumber ) { this->m_RadiusStepNumber = stepNumber; } template< class TOutputImage > void MultiGaussianImageSource< TOutputImage > ::SetNumberOfGausssians( unsigned int n ) { this->m_NumberOfGaussians = n; } template< class TOutputImage > const double MultiGaussianImageSource< TOutputImage > ::GetMaxMeanValue() const { return m_MeanValue; } //---------------------------------------------------------------------------- template< class TOutputImage > -const typename MultiGaussianImageSource< TOutputImage >::IndexType +const typename MultiGaussianImageSource< TOutputImage >::ItkVectorType MultiGaussianImageSource< TOutputImage > ::GetSphereMidpoint() const { - return m_ShpereMidpoint; + return m_SphereMidpoint; } template< class TOutputImage > const double MultiGaussianImageSource< TOutputImage > ::MultiGaussianFunctionValueAtPoint(double x, double y, double z) { double summand0, summand1, summand2, power, value = 0; + // the for-loop represent the sum of the gaussian function for(unsigned int n =0; n < m_NumberOfGaussians; ++n) { summand0 = (x - m_CenterX[n]) / m_SigmaX[n]; summand1 = (y - m_CenterY[n]) / m_SigmaY[n]; summand2 = (z - m_CenterZ[n]) / m_SigmaZ[n]; power = summand0 * summand0 + summand1 * summand1 + summand2 * summand2; value = value + m_Altitude[n] * pow(itk::Math::e, -0.5 * power); } return value; } template< class TOutputImage > void MultiGaussianImageSource< TOutputImage > ::AddGaussian( VectorType x, VectorType y, VectorType z, VectorType sx, VectorType sy, VectorType sz, VectorType altitude) { for(unsigned int i = 0; i < x.size(); ++i) { m_CenterX.push_back(x[i]); m_CenterY.push_back(y[i]); m_CenterZ.push_back(z[i]); m_SigmaX.push_back(sx[i]); m_SigmaY.push_back(sy[i]); m_SigmaZ.push_back(sz[i]); m_Altitude.push_back(altitude[i]); } } //---------------------------------------------------------------------------- template< typename TOutputImage > void MultiGaussianImageSource< TOutputImage > ::GenerateOutputInformation() { TOutputImage *output; IndexType index; index.Fill(0); output = this->GetOutput(0); typename TOutputImage::RegionType largestPossibleRegion; largestPossibleRegion.SetSize(this->m_Size); largestPossibleRegion.SetIndex(index); output->SetLargestPossibleRegion(largestPossibleRegion); output->SetSpacing(m_Spacing); output->SetOrigin(m_Origin); } //---------------------------------------------------------------------------- template< typename TOutputImage > void MultiGaussianImageSource< TOutputImage > ::GenerateData() { itkDebugMacro(<< "Generating a image of scalars "); - double valueReal, value; + double valueReal; IndexType index; - // Support progress methods/callbacks - //ProgressReporter progress( this, threadId, outputRegionForThread.GetNumberOfPixels() ); - - // typedef typename TOutputImage::PixelType scalarType; - //typename TOutputImage::Pointer image = this->GetOutput(0); - m_Image = this->GetOutput(0); - m_Image->SetBufferedRegion( m_Image->GetRequestedRegion() ); - m_Image->Allocate(); - IteratorType imageIt(m_Image, m_Image->GetLargestPossibleRegion()); + 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(); - m_Image->TransformIndexToPhysicalPoint(imageIt.GetIndex(), globalCoordinate); + image->TransformIndexToPhysicalPoint(imageIt.GetIndex(), globalCoordinate); valueReal = MultiGaussianFunctionValueAtPoint(globalCoordinate[0], globalCoordinate[1] ,globalCoordinate[2]); imageIt.Set(valueReal); - //progress.CompletedPixel(); } } //---------------------------------------------------------------------------- +/* + This class allows by the method CalculateMidpointAndMeanValue() to find a sphere with a specified radius that has a maximal mean value over all sphere with that radius with midpoint inside or at the boundary of the image. The parameter RadiusStepNumber controls the accuracy of that calculation (the higher the value the higher the exactness). + The algorithm works as follows: + 1. the first three for-loops traverse the image and assume the current point to be the wanted sphere midpoint + 2. calculate the mean value for that sphere (use sphere coordinates): + 2.1. traverse the radius of the sphere with step size Radius divided by RadiusStepNumber (the for-loop with index i) + 2.2. define a variable dist, which gives a approximately distance between the points at the sphere surface + (here we take such a distance, that on the smallest sphere are located 8 points) + 2.3. calculate the angles so that the points are equally spaced on the surface (for-loops with indexes j and k) + 2.3. for all radius length add the values at the points on the sphere and divide by the number of points added + (the values at each point is calculate with the method MultiGaussianFunctionValueAtPoint()) + 3. Compare with the until-now-found-maximum and take the bigger one +*/ template< typename TOutputImage > void MultiGaussianImageSource< TOutputImage > ::CalculateMidpointAndMeanValue() { itkDebugMacro(<< "Generating a image of scalars "); int numberSummand = 0, angleStepNumberOverTwo; - double valueReal, meanValueTemp, valueAdd, value, x, y, z, temp; + 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; - IteratorType imageIt(m_Image, m_Image->GetLargestPossibleRegion()); - PointType globalCoordinate; - IndexType index; - - for(imageIt.GoToBegin(); !imageIt.IsAtEnd(); ++imageIt) + for(unsigned int index0 = 0; index0 < m_Size[0]; ++index0) { - numberSummand = 1; - index = imageIt.GetIndex(); - m_Image->TransformIndexToPhysicalPoint(imageIt.GetIndex(), globalCoordinate); - value = imageIt.Value(); - ri = riStep; - for(unsigned int i = 0; i < m_RadiusStepNumber; ++i) + for(unsigned int index1 = 0; index1 < m_Size[1]; ++index1) { - angleStepNumberOverTwo = static_cast( itk::Math::pi * ri / dist); - fij = 0; - fijStep = itk::Math::pi / angleStepNumberOverTwo; - for(unsigned int j = 0; j <= angleStepNumberOverTwo; ++j) // from 0 to pi + for(unsigned int index2 = 0; index2 < m_Size[2]; ++index2) { - 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 + numberSummand = 0; + value = 0.0; + ri = riStep; + for(unsigned int i = 0; i < m_RadiusStepNumber; ++i) { - x = temp * cos(psik); - y = temp * sin(psik); - numberSummand++; - valueAdd = MultiGaussianFunctionValueAtPoint(x + globalCoordinate[0], y + globalCoordinate[1], z + globalCoordinate[2]); - value = value + valueAdd; + angleStepNumberOverTwo = static_cast( itk::Math::pi * ri / dist); + fij = 0; + fijStep = itk::Math::pi / angleStepNumberOverTwo; + for(unsigned int j = 0; j <= angleStepNumberOverTwo; ++j) // from 0 to pi + { + z = ri * cos(fij); + psikStep = 2.0 * itk::Math::pi / (2.0 * angleStepNumberOverTwo); + psik = -itk::Math::pi + psikStep; + + temp = ri * sin(fij); + for(unsigned int k = 0; k < 2 * angleStepNumberOverTwo; ++k) // from -pi to pi + { + x = temp * cos(psik); + y = temp * sin(psik); + numberSummand++; + valueAdd = MultiGaussianFunctionValueAtPoint(x + index0, y + index1, z + index2); + value = value + valueAdd; + psik = psik + psikStep; + } + fij = fij + fijStep; + } + ri = ri + riStep; + } - psik = psik + psikStep; + 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 ); } - fij = fij + fijStep ; } - ri = ri + riStep; - } - - meanValueTemp = value / numberSummand; - if(meanValueTemp > m_MeanValue) - { - m_MeanValue = meanValueTemp; - m_ShpereMidpoint = index; } } } - } // end namespace itk #endif diff --git a/Modules/ImageStatistics/Testing/mitkMultiGaussianTest.cpp b/Modules/ImageStatistics/Testing/mitkMultiGaussianTest.cpp index f1b61935d3..dd0787010e 100644 --- a/Modules/ImageStatistics/Testing/mitkMultiGaussianTest.cpp +++ b/Modules/ImageStatistics/Testing/mitkMultiGaussianTest.cpp @@ -1,93 +1,93 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkTestingMacros.h" #include "itkMultiGaussianImageSource.h" #include #include int mitkMultiGaussianTest(int, char* []) { // always start with this! MITK_TEST_BEGIN("mitkMultiGaussianTest") typedef double PixelType; const unsigned int Dimension = 3; typedef itk::Image ImageType; typedef itk::MultiGaussianImageSource< ImageType > MultiGaussianImageSource; itk::MultiGaussianImageSource< ImageType >::VectorType centerXVec, centerYVec, centerZVec, sigmaXVec, sigmaYVec, sigmaZVec, altitudeVec; MultiGaussianImageSource::Pointer gaussianGenerator = MultiGaussianImageSource::New(); ImageType::SizeValueType size[3]; size[0] = 50; size[1] = 50; size[2] = 50; srand (time(NULL)); unsigned int numberOfGaussians = 7; unsigned int minWidthOfGaussian = (size[0] + size[1] + size[2]) / 27; // A ninth of the mean image size unsigned int maxWidthOfGaussian = (size[0] + size[1] + size[2]) / 9; // One-third of the mean image size unsigned int minAltitudeOfGaussian = 5; unsigned int maxAltitudeOfGaussian = 200; double centerX, centerY, centerZ, sigmaX, sigmaY, sigmaZ, altitude; gaussianGenerator->SetSize( size ); gaussianGenerator->SetSpacing( 1 ); gaussianGenerator->SetRadiusStepNumber(5); gaussianGenerator->SetRadius(pow(itk::Math::one_over_pi * 0.75 , 1.0 / 3.0) * 10); gaussianGenerator->SetNumberOfGausssians(numberOfGaussians); - std::ofstream myfile; - myfile.open ("C:/temp/tempParameter3.txt"); - myfile << " CentX \t" << "Y \t" << "Z \t" << "SigX \t" << "Y \t" << "Z \t" << "Altit\n"; + // std::ofstream myfile; + // myfile.open ("C:/temp/tempParameter3.txt"); + // myfile << " CentX \t" << "Y \t" << "Z \t" << "SigX \t" << "Y \t" << "Z \t" << "Altit\n"; - for( unsigned int i = 0; i < numberOfGaussians; ++i) + int numberAddGaussian = numberOfGaussians; + for( unsigned int i = 0; i < numberAddGaussian; ++i) { centerX = rand() % size[0]; centerY = rand() % size[1]; centerZ = rand() % size[2]; sigmaX = minWidthOfGaussian + rand() % (maxWidthOfGaussian - minWidthOfGaussian); sigmaY = minWidthOfGaussian + rand() % (maxWidthOfGaussian - minWidthOfGaussian); sigmaZ = minWidthOfGaussian + rand() % (maxWidthOfGaussian - minWidthOfGaussian); altitude = minAltitudeOfGaussian + rand() % (maxAltitudeOfGaussian - minAltitudeOfGaussian); //gaussianGenerator->AddGaussian(centerX, centerY, centerZ, sigmaX, sigmaY, sigmaZ, altitude); centerXVec.push_back(centerX); centerYVec.push_back(centerY); centerZVec.push_back(centerZ); sigmaXVec.push_back(sigmaX); sigmaYVec.push_back(sigmaY); sigmaZVec.push_back(sigmaZ); altitudeVec.push_back(altitude); - - myfile <AddGaussian(centerXVec, centerYVec, centerZVec, sigmaXVec, sigmaYVec, sigmaZVec, altitudeVec); gaussianGenerator->Update(); gaussianGenerator->CalculateMidpointAndMeanValue(); std::cout << "Sphere radius is: " << gaussianGenerator->GetRadius() << std::endl; std::cout << "Sphere midpoint is: " << gaussianGenerator->GetSphereMidpoint() << std::endl; std::cout << "Mean value is: " << gaussianGenerator->GetMaxMeanValue() << std::endl; ImageType::Pointer gaussianImage = gaussianGenerator->GetOutput(); //File writer typedef itk::ImageFileWriter< ImageType > WriterType; WriterType::Pointer writer = WriterType::New(); writer->SetFileName( "C:/temp/tempImage33.nrrd" ); writer->SetInput( gaussianImage ); writer->Update(); MITK_TEST_END() } \ No newline at end of file