1 /*=========================================================================
2 *
3 * Copyright Insight Software Consortium
4 *
5 * Licensed under the Apache License, Version 2.0 (the "License");
6 * you may not use this file except in compliance with the License.
7 * You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0.txt
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 *
17 *=========================================================================*/
18 #ifndef itkGaussianSpatialFunction_hxx
19 #define itkGaussianSpatialFunction_hxx
20
21 #include <cmath>
22 #include "itkMath.h"
23 #include "itkGaussianSpatialFunction.h"
24
25 namespace itk
26 {
27 template< typename TOutput, unsigned int VImageDimension, typename TInput >
28 GaussianSpatialFunction< TOutput, VImageDimension, TInput >
GaussianSpatialFunction()29 ::GaussianSpatialFunction()
30
31 {
32 m_Mean = ArrayType::Filled(10.0);
33 m_Sigma = ArrayType::Filled(5.0);
34 }
35
36 template< typename TOutput, unsigned int VImageDimension, typename TInput >
37 typename GaussianSpatialFunction< TOutput, VImageDimension, TInput >::OutputType
38 GaussianSpatialFunction< TOutput, VImageDimension, TInput >
Evaluate(const TInput & position) const39 ::Evaluate(const TInput & position) const
40 {
41 // We have to compute the Gaussian in several stages, because of the
42 // n-dimensional generalization
43
44 // Normalizing the Gaussian is important for statistical applications
45 // but is generally not desirable for creating images because of the
46 // very small numbers involved (would need to use doubles)
47 double prefixDenom = 1.0;
48
49 if ( m_Normalized )
50 {
51 const double squareRootOfTwoPi = std::sqrt(2.0 * itk::Math::pi);
52
53 for ( unsigned int i = 0; i < VImageDimension; ++i )
54 {
55 prefixDenom *= m_Sigma[i] * squareRootOfTwoPi;
56 }
57 }
58
59 double suffixExp = 0;
60
61 for ( unsigned int i = 0; i < VImageDimension; ++i )
62 {
63 suffixExp += ( position[i] - m_Mean[i] ) * ( position[i] - m_Mean[i] )
64 / ( 2 * m_Sigma[i] * m_Sigma[i] );
65 }
66
67 const double value = m_Scale * ( 1 / prefixDenom ) * std::exp(-1 * suffixExp);
68
69 return static_cast< TOutput >( value );
70 }
71
72 template< typename TOutput, unsigned int VImageDimension, typename TInput >
73 void
74 GaussianSpatialFunction< TOutput, VImageDimension, TInput >
PrintSelf(std::ostream & os,Indent indent) const75 ::PrintSelf(std::ostream & os, Indent indent) const
76 {
77 Superclass::PrintSelf(os, indent);
78
79 os << indent << "Sigma: " << m_Sigma << std::endl;
80 os << indent << "Mean: " << m_Mean << std::endl;
81 os << indent << "Scale: " << m_Scale << std::endl;
82 os << indent << "Normalized?: " << m_Normalized << std::endl;
83 }
84 } // end namespace itk
85
86 #endif
87