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 itkEllipsoidInteriorExteriorSpatialFunction_hxx
19 #define itkEllipsoidInteriorExteriorSpatialFunction_hxx
20
21 #include "itkEllipsoidInteriorExteriorSpatialFunction.h"
22 #include <cmath>
23
24 namespace itk
25 {
26 template< unsigned int VDimension, typename TInput >
27 EllipsoidInteriorExteriorSpatialFunction< VDimension, TInput >
EllipsoidInteriorExteriorSpatialFunction()28 ::EllipsoidInteriorExteriorSpatialFunction()
29 {
30 m_Orientations = nullptr;
31 m_Axes.Fill(1.0f); // Lengths of ellipsoid axes.
32 m_Center.Fill(0.0f); // Origin of ellipsoid
33 }
34
35 template< unsigned int VDimension, typename TInput >
36 EllipsoidInteriorExteriorSpatialFunction< VDimension, TInput >
~EllipsoidInteriorExteriorSpatialFunction()37 ::~EllipsoidInteriorExteriorSpatialFunction()
38 {
39 unsigned int i;
40
41 if ( m_Orientations )
42 {
43 for ( i = 0; i < VDimension; i++ )
44 {
45 delete[] m_Orientations[i];
46 }
47 delete[] m_Orientations;
48 }
49 }
50
51 template< unsigned int VDimension, typename TInput >
52 typename EllipsoidInteriorExteriorSpatialFunction< VDimension, TInput >::OutputType
53 EllipsoidInteriorExteriorSpatialFunction< VDimension, TInput >
Evaluate(const InputType & position) const54 ::Evaluate(const InputType & position) const
55 {
56 double distanceSquared = 0;
57
58 Vector< double, VDimension > orientationVector;
59 Vector< double, VDimension > pointVector;
60
61 // Project the position onto each of the axes, normalize by axis length,
62 // and determine whether position is inside ellipsoid. The length of axis0,
63 // m_Axis[0] is orientated in the direction of m_Orientations[0].
64 for ( unsigned int i = 0; i < VDimension; i++ )
65 {
66 pointVector[i] = position[i] - m_Center[i];
67 }
68
69 for ( unsigned int i = 0; i < VDimension; i++ )
70 {
71 for ( unsigned int j = 0; j < VDimension; j++ )
72 {
73 orientationVector[j] = m_Orientations[i][j];
74 }
75 distanceSquared += std::pow( static_cast< double >( ( orientationVector * pointVector ) / ( .5 * m_Axes[i] ) ),
76 static_cast< double >( 2 ) );
77 }
78
79 if ( distanceSquared <= 1 )
80 {
81 return 1; // Inside the ellipsoid.
82 }
83 //Default return value assumes outside the ellipsoid
84 return 0; // Outside the ellipsoid.
85 }
86
87 template< unsigned int VDimension, typename TInput >
88 void EllipsoidInteriorExteriorSpatialFunction< VDimension, TInput >
SetOrientations(const OrientationType & orientations)89 ::SetOrientations(const OrientationType & orientations)
90 {
91 unsigned int i, j;
92
93 // Initialize orientation vectors.
94 if ( m_Orientations )
95 {
96 for ( i = 0; i < VDimension; i++ )
97 {
98 delete[] m_Orientations[i];
99 }
100 delete[] m_Orientations;
101 }
102 m_Orientations = new double *[VDimension];
103 for ( i = 0; i < VDimension; i++ )
104 {
105 m_Orientations[i] = new double[VDimension];
106 }
107
108 // Set orientation vectors (must be orthogonal).
109 for ( i = 0; i < VDimension; i++ )
110 {
111 for ( j = 0; j < VDimension; j++ )
112 {
113 m_Orientations[i][j] = orientations[i][j];
114 }
115 }
116 }
117
118 template< unsigned int VDimension, typename TInput >
119 void EllipsoidInteriorExteriorSpatialFunction< VDimension, TInput >
PrintSelf(std::ostream & os,Indent indent) const120 ::PrintSelf(std::ostream & os, Indent indent) const
121 {
122 unsigned int i, j;
123
124 Superclass::PrintSelf(os, indent);
125
126 os << indent << "Lengths of Ellipsoid Axes: " << m_Axes << std::endl;
127 os << indent << "Origin of Ellipsoid: " << m_Center << std::endl;
128 if ( m_Orientations )
129 {
130 os << indent << "Orientations: " << std::endl;
131 for ( i = 0; i < VDimension; i++ )
132 {
133 for ( j = 0; j < VDimension; j++ )
134 {
135 os << indent << indent << m_Orientations[i][j] << " ";
136 }
137 os << std::endl;
138 }
139 }
140 }
141 } // end namespace itk
142
143 #endif
144