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 itkFrustumSpatialFunction_hxx
19 #define itkFrustumSpatialFunction_hxx
20
21 #include "itkFrustumSpatialFunction.h"
22
23 namespace itk
24 {
25 template< unsigned int VDimension, typename TInput >
FrustumSpatialFunction()26 FrustumSpatialFunction< VDimension, TInput >::FrustumSpatialFunction() :
27 m_RotationPlane( RotateInXZPlane )
28 {
29 m_Apex.Fill( 0.0f );
30 }
31
32 template< unsigned int VDimension, typename TInput >
33 typename FrustumSpatialFunction< VDimension, TInput >::OutputType
34 FrustumSpatialFunction< VDimension, TInput >
Evaluate(const InputType & position) const35 ::Evaluate(const InputType & position) const
36 {
37 using PointType = InputType;
38 using VectorType = typename PointType::VectorType;
39
40 VectorType relativePosition = position - m_Apex;
41 const double distanceToApex = relativePosition.GetNorm();
42
43 // Check Top and Bottom planes. If the angle is negative, the top plane
44 // value may be less than the bottom plane, which is still correct.
45 if ( m_TopPlane <= m_BottomPlane )
46 {
47 if ( distanceToApex < m_TopPlane
48 || distanceToApex > m_BottomPlane )
49 {
50 return 0;
51 }
52 }
53 else
54 {
55 if ( distanceToApex > m_TopPlane
56 || distanceToApex < m_BottomPlane )
57 {
58 return 0;
59 }
60 }
61
62 if ( m_RotationPlane == RotateInXZPlane )
63 {
64 const double dx = relativePosition[0];
65 const double dy = relativePosition[1];
66 const double dz = relativePosition[2];
67
68 const double distanceXZ = std::sqrt(dx * dx + dz * dz);
69
70 const double deg2rad = std::atan(1.0) / 45.0;
71
72 // Check planes along Y
73 const double angleY = std::atan2(dy, distanceXZ);
74 if ( std::fabs(angleY) > m_ApertureAngleY * deg2rad )
75 {
76 return 0;
77 }
78
79 // Check planes along X
80 const double angleX = std::atan2(dx, dz);
81
82 if ( std::cos(angleX + ( 180.0 + m_AngleZ ) * deg2rad) <
83 std::cos(deg2rad * m_ApertureAngleX) )
84 {
85 return 0;
86 }
87
88 return 1;
89 }
90 else if ( m_RotationPlane == RotateInYZPlane )
91 {
92 const double dx = relativePosition[0];
93 const double dy = relativePosition[1];
94 const double dz = relativePosition[2];
95
96 const double distanceYZ = std::sqrt(dy * dy + dz * dz);
97
98 const double deg2rad = std::atan(1.0) / 45.0;
99
100 // Check planes along X
101 const double angleX = std::atan2(dx, distanceYZ);
102 if ( std::fabs(angleX) > m_ApertureAngleX * deg2rad )
103 {
104 return 0;
105 }
106
107 // Check planes along Y
108 const double angleY = std::atan2(dy, dz);
109
110 if ( std::cos(angleY + ( 180.0 + m_AngleZ ) * deg2rad) <
111 std::cos(deg2rad * m_ApertureAngleY) )
112 {
113 return 0;
114 }
115
116 return 1;
117 }
118 else
119 {
120 itkExceptionMacro(<< "Rotation plane not set or set to an unsupported value!");
121 }
122 }
123
124 template< unsigned int VDimension, typename TInput >
125 void
PrintSelf(std::ostream & os,Indent indent) const126 FrustumSpatialFunction< VDimension, TInput >::PrintSelf(std::ostream & os, Indent indent) const
127 {
128 Superclass::PrintSelf(os, indent);
129
130 os << indent << "Apex: " << m_Apex << std::endl;
131 os << indent << "AngleZ: " << m_AngleZ << std::endl;
132 os << indent << "ApertureAngleX: " << m_ApertureAngleX << std::endl;
133 os << indent << "ApertureAngleY: " << m_ApertureAngleY << std::endl;
134 os << indent << "TopPlane: " << m_TopPlane << std::endl;
135 os << indent << "BottomPlane: " << m_BottomPlane << std::endl;
136 os << indent << "RotationPlane: " << m_RotationPlane << std::endl;
137 }
138 } // end namespace itk
139
140 #endif
141