diff --git a/Modules/ImageStatistics/Testing/files.cmake b/Modules/ImageStatistics/Testing/files.cmake index 01792788aa..cdfeb842e7 100644 --- a/Modules/ImageStatistics/Testing/files.cmake +++ b/Modules/ImageStatistics/Testing/files.cmake @@ -1,5 +1,6 @@ set(MODULE_TESTS mitkImageStatisticsCalculatorTest.cpp mitkPointSetStatisticsCalculatorTest.cpp mitkPointSetDifferenceStatisticsCalculatorTest.cpp + mitkMultiGaussianTest.cpp ) diff --git a/Modules/ImageStatistics/Testing/itkMultiGaussianImageSource.h b/Modules/ImageStatistics/Testing/itkMultiGaussianImageSource.h new file mode 100644 index 0000000000..64d24106c0 --- /dev/null +++ b/Modules/ImageStatistics/Testing/itkMultiGaussianImageSource.h @@ -0,0 +1,184 @@ +/*========================================================================= + * + * 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 + * + * + * \ingroup DataSources MultiThreaded + * \ingroup ITKTestKernel + * + * \wiki + * \wikiexample{SimpleOperations/RandomImageSource,Produce an image of noise} + * \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; + + /** 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; + + + virtual unsigned int GetNumberOfGaussians() const; + virtual const RadiusType GetRadius() const; + virtual void SetRadius( RadiusType radius ); + virtual const int GetRadiusStepNumber() const; + virtual void SetRadiusStepNumber( unsigned int stepNumber ); + virtual const double GetMaxMeanValue() const; + virtual const IndexType GetSphereMidpoint() const; + virtual const double MultiGaussianFunctionValueAtPoint(double , double, double); + virtual void AddGaussian( VectorType, VectorType, VectorType, VectorType, VectorType, VectorType, VectorType); + virtual void SetNumberOfGausssians( unsigned int ); + 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 + 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 new file mode 100644 index 0000000000..400ca63780 --- /dev/null +++ b/Modules/ImageStatistics/Testing/itkMultiGaussianImageSource.hxx @@ -0,0 +1,438 @@ +/*========================================================================= + * + * 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_NumberOfGaussians = 1; + 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 +MultiGaussianImageSource< TOutputImage > +::GetSphereMidpoint() const +{ + return m_ShpereMidpoint; +} + +template< class TOutputImage > +const double +MultiGaussianImageSource< TOutputImage > +::MultiGaussianFunctionValueAtPoint(double x, double y, double z) +{ + double summand0, summand1, summand2, power, value = 0; + 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; + 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()); + PointType globalCoordinate; + + for(imageIt.GoToBegin(); !imageIt.IsAtEnd(); ++imageIt) + { + valueReal = 0.0; + index = imageIt.GetIndex(); + m_Image->TransformIndexToPhysicalPoint(imageIt.GetIndex(), globalCoordinate); + valueReal = MultiGaussianFunctionValueAtPoint(globalCoordinate[0], globalCoordinate[1] ,globalCoordinate[2]); + imageIt.Set(valueReal); + //progress.CompletedPixel(); + } +} + +//---------------------------------------------------------------------------- +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 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) + { + 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) + { + 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 + 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_MeanValue = meanValueTemp; + m_ShpereMidpoint = index; + } + } +} + +} // end namespace itk +#endif diff --git a/Modules/ImageStatistics/Testing/mitkMultiGaussianTest.cpp b/Modules/ImageStatistics/Testing/mitkMultiGaussianTest.cpp new file mode 100644 index 0000000000..f1b61935d3 --- /dev/null +++ b/Modules/ImageStatistics/Testing/mitkMultiGaussianTest.cpp @@ -0,0 +1,93 @@ +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center, +Division of Medical and Biological Informatics. +All rights reserved. + +This software is distributed WITHOUT ANY WARRANTY; without +even the implied warranty of MERCHANTABILITY or FITNESS FOR +A PARTICULAR PURPOSE. + +See LICENSE.txt or http://www.mitk.org for details. + +===================================================================*/ + +#include "mitkTestingMacros.h" +#include "itkMultiGaussianImageSource.h" + +#include +#include + +int mitkMultiGaussianTest(int, char* []) +{ + // always start with this! + MITK_TEST_BEGIN("mitkMultiGaussianTest") + + typedef double PixelType; + const unsigned int Dimension = 3; + typedef itk::Image ImageType; + typedef itk::MultiGaussianImageSource< ImageType > MultiGaussianImageSource; + itk::MultiGaussianImageSource< ImageType >::VectorType centerXVec, centerYVec, centerZVec, sigmaXVec, sigmaYVec, sigmaZVec, altitudeVec; + MultiGaussianImageSource::Pointer gaussianGenerator = MultiGaussianImageSource::New(); + ImageType::SizeValueType size[3]; + size[0] = 50; + size[1] = 50; + size[2] = 50; + srand (time(NULL)); + unsigned int numberOfGaussians = 7; + unsigned int minWidthOfGaussian = (size[0] + size[1] + size[2]) / 27; // A ninth of the mean image size + unsigned int maxWidthOfGaussian = (size[0] + size[1] + size[2]) / 9; // One-third of the mean image size + unsigned int minAltitudeOfGaussian = 5; + unsigned int maxAltitudeOfGaussian = 200; + double centerX, centerY, centerZ, sigmaX, sigmaY, sigmaZ, altitude; + + gaussianGenerator->SetSize( size ); + gaussianGenerator->SetSpacing( 1 ); + gaussianGenerator->SetRadiusStepNumber(5); + gaussianGenerator->SetRadius(pow(itk::Math::one_over_pi * 0.75 , 1.0 / 3.0) * 10); + gaussianGenerator->SetNumberOfGausssians(numberOfGaussians); + std::ofstream myfile; + myfile.open ("C:/temp/tempParameter3.txt"); + myfile << " CentX \t" << "Y \t" << "Z \t" << "SigX \t" << "Y \t" << "Z \t" << "Altit\n"; + + for( unsigned int i = 0; i < numberOfGaussians; ++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