Chris@16
|
1 // Boost.Geometry (aka GGL, Generic Geometry Library)
|
Chris@16
|
2
|
Chris@16
|
3 // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands.
|
Chris@16
|
4 // Copyright (c) 2008-2012 Bruno Lalande, Paris, France.
|
Chris@16
|
5 // Copyright (c) 2009-2012 Mateusz Loskot, London, UK.
|
Chris@16
|
6
|
Chris@16
|
7 // Parts of Boost.Geometry are redesigned from Geodan's Geographic Library
|
Chris@16
|
8 // (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands.
|
Chris@16
|
9
|
Chris@16
|
10 // Use, modification and distribution is subject to the Boost Software License,
|
Chris@16
|
11 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
|
Chris@16
|
12 // http://www.boost.org/LICENSE_1_0.txt)
|
Chris@16
|
13
|
Chris@16
|
14 #ifndef BOOST_GEOMETRY_UTIL_MATH_HPP
|
Chris@16
|
15 #define BOOST_GEOMETRY_UTIL_MATH_HPP
|
Chris@16
|
16
|
Chris@16
|
17 #include <cmath>
|
Chris@16
|
18 #include <limits>
|
Chris@16
|
19
|
Chris@16
|
20 #include <boost/math/constants/constants.hpp>
|
Chris@16
|
21
|
Chris@16
|
22 #include <boost/geometry/util/select_most_precise.hpp>
|
Chris@16
|
23
|
Chris@16
|
24 namespace boost { namespace geometry
|
Chris@16
|
25 {
|
Chris@16
|
26
|
Chris@16
|
27 namespace math
|
Chris@16
|
28 {
|
Chris@16
|
29
|
Chris@16
|
30 #ifndef DOXYGEN_NO_DETAIL
|
Chris@16
|
31 namespace detail
|
Chris@16
|
32 {
|
Chris@16
|
33
|
Chris@16
|
34
|
Chris@16
|
35 template <typename Type, bool IsFloatingPoint>
|
Chris@16
|
36 struct equals
|
Chris@16
|
37 {
|
Chris@16
|
38 static inline bool apply(Type const& a, Type const& b)
|
Chris@16
|
39 {
|
Chris@16
|
40 return a == b;
|
Chris@16
|
41 }
|
Chris@16
|
42 };
|
Chris@16
|
43
|
Chris@16
|
44 template <typename Type>
|
Chris@16
|
45 struct equals<Type, true>
|
Chris@16
|
46 {
|
Chris@16
|
47 static inline Type get_max(Type const& a, Type const& b, Type const& c)
|
Chris@16
|
48 {
|
Chris@16
|
49 return (std::max)((std::max)(a, b), c);
|
Chris@16
|
50 }
|
Chris@16
|
51
|
Chris@16
|
52 static inline bool apply(Type const& a, Type const& b)
|
Chris@16
|
53 {
|
Chris@16
|
54 if (a == b)
|
Chris@16
|
55 {
|
Chris@16
|
56 return true;
|
Chris@16
|
57 }
|
Chris@16
|
58
|
Chris@16
|
59 // See http://www.parashift.com/c++-faq-lite/newbie.html#faq-29.17,
|
Chris@16
|
60 // FUTURE: replace by some boost tool or boost::test::close_at_tolerance
|
Chris@16
|
61 return std::abs(a - b) <= std::numeric_limits<Type>::epsilon() * get_max(std::abs(a), std::abs(b), 1.0);
|
Chris@16
|
62 }
|
Chris@16
|
63 };
|
Chris@16
|
64
|
Chris@16
|
65 template <typename Type, bool IsFloatingPoint>
|
Chris@16
|
66 struct smaller
|
Chris@16
|
67 {
|
Chris@16
|
68 static inline bool apply(Type const& a, Type const& b)
|
Chris@16
|
69 {
|
Chris@16
|
70 return a < b;
|
Chris@16
|
71 }
|
Chris@16
|
72 };
|
Chris@16
|
73
|
Chris@16
|
74 template <typename Type>
|
Chris@16
|
75 struct smaller<Type, true>
|
Chris@16
|
76 {
|
Chris@16
|
77 static inline bool apply(Type const& a, Type const& b)
|
Chris@16
|
78 {
|
Chris@16
|
79 if (equals<Type, true>::apply(a, b))
|
Chris@16
|
80 {
|
Chris@16
|
81 return false;
|
Chris@16
|
82 }
|
Chris@16
|
83 return a < b;
|
Chris@16
|
84 }
|
Chris@16
|
85 };
|
Chris@16
|
86
|
Chris@16
|
87
|
Chris@16
|
88 template <typename Type, bool IsFloatingPoint>
|
Chris@16
|
89 struct equals_with_epsilon : public equals<Type, IsFloatingPoint> {};
|
Chris@16
|
90
|
Chris@16
|
91
|
Chris@16
|
92 /*!
|
Chris@16
|
93 \brief Short construct to enable partial specialization for PI, currently not possible in Math.
|
Chris@16
|
94 */
|
Chris@16
|
95 template <typename T>
|
Chris@16
|
96 struct define_pi
|
Chris@16
|
97 {
|
Chris@16
|
98 static inline T apply()
|
Chris@16
|
99 {
|
Chris@16
|
100 // Default calls Boost.Math
|
Chris@16
|
101 return boost::math::constants::pi<T>();
|
Chris@16
|
102 }
|
Chris@16
|
103 };
|
Chris@16
|
104
|
Chris@16
|
105 template <typename T>
|
Chris@16
|
106 struct relaxed_epsilon
|
Chris@16
|
107 {
|
Chris@16
|
108 static inline T apply(const T& factor)
|
Chris@16
|
109 {
|
Chris@16
|
110 return factor * std::numeric_limits<T>::epsilon();
|
Chris@16
|
111 }
|
Chris@16
|
112 };
|
Chris@16
|
113
|
Chris@16
|
114
|
Chris@16
|
115 } // namespace detail
|
Chris@16
|
116 #endif
|
Chris@16
|
117
|
Chris@16
|
118
|
Chris@16
|
119 template <typename T>
|
Chris@16
|
120 inline T pi() { return detail::define_pi<T>::apply(); }
|
Chris@16
|
121
|
Chris@16
|
122 template <typename T>
|
Chris@16
|
123 inline T relaxed_epsilon(T const& factor)
|
Chris@16
|
124 {
|
Chris@16
|
125 return detail::relaxed_epsilon<T>::apply(factor);
|
Chris@16
|
126 }
|
Chris@16
|
127
|
Chris@16
|
128
|
Chris@16
|
129 // Maybe replace this by boost equals or boost ublas numeric equals or so
|
Chris@16
|
130
|
Chris@16
|
131 /*!
|
Chris@16
|
132 \brief returns true if both arguments are equal.
|
Chris@16
|
133 \ingroup utility
|
Chris@16
|
134 \param a first argument
|
Chris@16
|
135 \param b second argument
|
Chris@16
|
136 \return true if a == b
|
Chris@16
|
137 \note If both a and b are of an integral type, comparison is done by ==.
|
Chris@16
|
138 If one of the types is floating point, comparison is done by abs and
|
Chris@16
|
139 comparing with epsilon. If one of the types is non-fundamental, it might
|
Chris@16
|
140 be a high-precision number and comparison is done using the == operator
|
Chris@16
|
141 of that class.
|
Chris@16
|
142 */
|
Chris@16
|
143
|
Chris@16
|
144 template <typename T1, typename T2>
|
Chris@16
|
145 inline bool equals(T1 const& a, T2 const& b)
|
Chris@16
|
146 {
|
Chris@16
|
147 typedef typename select_most_precise<T1, T2>::type select_type;
|
Chris@16
|
148 return detail::equals
|
Chris@16
|
149 <
|
Chris@16
|
150 select_type,
|
Chris@16
|
151 boost::is_floating_point<select_type>::type::value
|
Chris@16
|
152 >::apply(a, b);
|
Chris@16
|
153 }
|
Chris@16
|
154
|
Chris@16
|
155 template <typename T1, typename T2>
|
Chris@16
|
156 inline bool equals_with_epsilon(T1 const& a, T2 const& b)
|
Chris@16
|
157 {
|
Chris@16
|
158 typedef typename select_most_precise<T1, T2>::type select_type;
|
Chris@16
|
159 return detail::equals_with_epsilon
|
Chris@16
|
160 <
|
Chris@16
|
161 select_type,
|
Chris@16
|
162 boost::is_floating_point<select_type>::type::value
|
Chris@16
|
163 >::apply(a, b);
|
Chris@16
|
164 }
|
Chris@16
|
165
|
Chris@16
|
166 template <typename T1, typename T2>
|
Chris@16
|
167 inline bool smaller(T1 const& a, T2 const& b)
|
Chris@16
|
168 {
|
Chris@16
|
169 typedef typename select_most_precise<T1, T2>::type select_type;
|
Chris@16
|
170 return detail::smaller
|
Chris@16
|
171 <
|
Chris@16
|
172 select_type,
|
Chris@16
|
173 boost::is_floating_point<select_type>::type::value
|
Chris@16
|
174 >::apply(a, b);
|
Chris@16
|
175 }
|
Chris@16
|
176
|
Chris@16
|
177 template <typename T1, typename T2>
|
Chris@16
|
178 inline bool larger(T1 const& a, T2 const& b)
|
Chris@16
|
179 {
|
Chris@16
|
180 typedef typename select_most_precise<T1, T2>::type select_type;
|
Chris@16
|
181 return detail::smaller
|
Chris@16
|
182 <
|
Chris@16
|
183 select_type,
|
Chris@16
|
184 boost::is_floating_point<select_type>::type::value
|
Chris@16
|
185 >::apply(b, a);
|
Chris@16
|
186 }
|
Chris@16
|
187
|
Chris@16
|
188
|
Chris@16
|
189
|
Chris@16
|
190 double const d2r = geometry::math::pi<double>() / 180.0;
|
Chris@16
|
191 double const r2d = 1.0 / d2r;
|
Chris@16
|
192
|
Chris@16
|
193 /*!
|
Chris@16
|
194 \brief Calculates the haversine of an angle
|
Chris@16
|
195 \ingroup utility
|
Chris@16
|
196 \note See http://en.wikipedia.org/wiki/Haversine_formula
|
Chris@16
|
197 haversin(alpha) = sin2(alpha/2)
|
Chris@16
|
198 */
|
Chris@16
|
199 template <typename T>
|
Chris@16
|
200 inline T hav(T const& theta)
|
Chris@16
|
201 {
|
Chris@16
|
202 T const half = T(0.5);
|
Chris@16
|
203 T const sn = sin(half * theta);
|
Chris@16
|
204 return sn * sn;
|
Chris@16
|
205 }
|
Chris@16
|
206
|
Chris@16
|
207 /*!
|
Chris@16
|
208 \brief Short utility to return the square
|
Chris@16
|
209 \ingroup utility
|
Chris@16
|
210 \param value Value to calculate the square from
|
Chris@16
|
211 \return The squared value
|
Chris@16
|
212 */
|
Chris@16
|
213 template <typename T>
|
Chris@16
|
214 inline T sqr(T const& value)
|
Chris@16
|
215 {
|
Chris@16
|
216 return value * value;
|
Chris@16
|
217 }
|
Chris@16
|
218
|
Chris@16
|
219 /*!
|
Chris@16
|
220 \brief Short utility to workaround gcc/clang problem that abs is converting to integer
|
Chris@16
|
221 and that older versions of MSVC does not support abs of long long...
|
Chris@16
|
222 \ingroup utility
|
Chris@16
|
223 */
|
Chris@16
|
224 template<typename T>
|
Chris@16
|
225 inline T abs(T const& value)
|
Chris@16
|
226 {
|
Chris@16
|
227 T const zero = T();
|
Chris@16
|
228 return value < zero ? -value : value;
|
Chris@16
|
229 }
|
Chris@16
|
230
|
Chris@16
|
231 /*!
|
Chris@16
|
232 \brief Short utility to calculate the sign of a number: -1 (negative), 0 (zero), 1 (positive)
|
Chris@16
|
233 \ingroup utility
|
Chris@16
|
234 */
|
Chris@16
|
235 template <typename T>
|
Chris@16
|
236 static inline int sign(T const& value)
|
Chris@16
|
237 {
|
Chris@16
|
238 T const zero = T();
|
Chris@16
|
239 return value > zero ? 1 : value < zero ? -1 : 0;
|
Chris@16
|
240 }
|
Chris@16
|
241
|
Chris@16
|
242
|
Chris@16
|
243 } // namespace math
|
Chris@16
|
244
|
Chris@16
|
245
|
Chris@16
|
246 }} // namespace boost::geometry
|
Chris@16
|
247
|
Chris@16
|
248 #endif // BOOST_GEOMETRY_UTIL_MATH_HPP
|