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 itkSymmetricEllipsoidInteriorExteriorSpatialFunction_hxx
19 #define itkSymmetricEllipsoidInteriorExteriorSpatialFunction_hxx
20
21 #include "itkSymmetricEllipsoidInteriorExteriorSpatialFunction.h"
22 #include <cmath>
23
24 namespace itk
25 {
26 template< unsigned int VDimension, typename TInput >
27 SymmetricEllipsoidInteriorExteriorSpatialFunction< VDimension, TInput >
SymmetricEllipsoidInteriorExteriorSpatialFunction()28 ::SymmetricEllipsoidInteriorExteriorSpatialFunction()
29 {
30 m_Center.Fill(0.0); // Origin of ellipsoid
31 m_Orientation.Fill(1.0); // Orientation of unique axis
32 }
33
34 template< unsigned int VDimension, typename TInput >
35 typename SymmetricEllipsoidInteriorExteriorSpatialFunction< VDimension, TInput >::OutputType
36 SymmetricEllipsoidInteriorExteriorSpatialFunction< VDimension, TInput >
Evaluate(const InputType & position) const37 ::Evaluate(const InputType & position) const
38 {
39 double uniqueTerm; // Term in ellipsoid equation for unique axis
40 double symmetricTerm; // Term in ellipsoid equation for symmetric axes
41
42 Vector< double, VDimension > pointVector;
43 Vector< double, VDimension > symmetricVector;
44
45 // Project the position onto the major axis, normalize by axis length,
46 // and determine whether position is inside ellipsoid.
47 for ( unsigned int i = 0; i < VDimension; i++ )
48 {
49 pointVector[i] = position[i] - m_Center[i];
50 }
51
52 uniqueTerm = std::pow( static_cast< double >( ( ( pointVector * m_Orientation ) / ( .5 * m_UniqueAxis ) ) ),
53 static_cast< double >( 2 ) );
54 symmetricVector = pointVector - ( m_Orientation * ( pointVector * m_Orientation ) );
55 symmetricTerm = std::pow(
56 static_cast< double >( ( ( symmetricVector.GetNorm() ) / ( .5 * m_SymmetricAxes ) ) ), static_cast< double >( 2 ) );
57
58 if ( ( uniqueTerm + symmetricTerm ) >= 0 && ( uniqueTerm + symmetricTerm ) <= 1 )
59 {
60 return 1; // Inside the ellipsoid.
61 }
62 //Default return value assumes outside the ellipsoid
63 return 0; // Outside the ellipsoid.
64 }
65
66 template< unsigned int VDimension, typename TInput >
67 void SymmetricEllipsoidInteriorExteriorSpatialFunction< VDimension, TInput >
PrintSelf(std::ostream & os,Indent indent) const68 ::PrintSelf(std::ostream & os, Indent indent) const
69 {
70 Superclass::PrintSelf(os, indent);
71
72 os << indent << "Origin of Ellipsoid: ";
73 os << m_Center << std::endl;
74 os << indent << "Unique Axis Orientation: ";
75 os << m_Orientation << std::endl;
76 os << indent << "Unique Axis Length: ";
77 os << m_UniqueAxis << std::endl;
78 os << indent << "Symmetric Axis Length: ";
79 os << m_SymmetricAxes << std::endl;
80 }
81
82 template< unsigned int VDimension, typename TInput >
83 void SymmetricEllipsoidInteriorExteriorSpatialFunction< VDimension, TInput >
SetOrientation(VectorType orientation,double uniqueAxis,double symmetricAxes)84 ::SetOrientation(VectorType orientation, double uniqueAxis, double symmetricAxes)
85 {
86 m_Orientation = orientation; // Orientation of unique axis of ellipsoid
87 m_SymmetricAxes = symmetricAxes; // Length of symmetric axes
88 m_UniqueAxis = uniqueAxis; // Length of unique axis
89 }
90 } // end namespace itk
91
92 #endif
93