Fréchet View  1.6.0
A Tool for Exploring Fréchet Distance Algorithms
numeric.h
Go to the documentation of this file.
1 
2 #ifndef FRECHET_VIEW_NUM_ERROR_H
3 #define FRECHET_VIEW_NUM_ERROR_H
4 
5 #include <limits>
6 #include <math.h>
7 #include <qdebug.h>
8 #include <types.h>
9 
10 #include <boost/config.hpp>
11 #include <boost/math/special_functions/next.hpp>
12 
13 #ifdef QUAD_ARITHMETIC
14 /* Quadruple precision floats, supported by the quadmath library.
15  Depends on platform/compiler:
16  Windows+MSC = NO
17  Windows+INTEL = YES
18  Linux+gcc = probably?
19  Linux+INTEL = probably?
20  MacOS+clang = NO
21  MacOS+gcc = NO (but should?)
22  MacOS+INTEL = probably?
23 */
24 # if defined(BOOST_INTEL) || defined(__GNUC__)
25 # include <boost/multiprecision/float128.hpp>
26  typedef boost::multiprecision::float128 float128;
27 #else
28  // fall-back to long double, which in turn may fall back to double
29 # pragma message ("no float128 support. falling back to double precision.")
30  typedef long double float128;
31 # endif
32 #endif
33 
42 typedef long double accurate_float;
43 
44 
45 namespace frechet {
46 
47 using namespace data;
48 
49 namespace numeric {
50 
52 static const constexpr double double_eps = std::numeric_limits<double>::epsilon();
53 
60 template<typename Number>
61 inline Number sq(const Number& x) { return x*x; }
62 
69 template<typename Float>
70 inline Float sqrt(const Float& x);
71 
75 template<>
76 inline double sqrt(const double& x) { return ::sqrt(x); }
77 
81 template<>
82 inline long double sqrt(const long double& x) { return ::sqrtl(x); }
83 
90 template<typename Number>
91 inline Number abs(const Number& x) { return (x<0) ? -x : x; }
92 
93 #if defined(BOOST_MP_USE_FLOAT128) || defined(BOOST_MP_USE_QUAD)
94 
97 template<>
98 inline float128 sqrt(const float128& x) { return sqrt(x); }
99 
100 inline const float128& min(const float128& a, const float128& b)
101 {
102  return (a<b) ? a:b;
103 }
104 #endif
105 
112 inline double round_up(double x, int ulps=1) {
113  return boost::math::float_advance(x,ulps);
114 }
115 
125 template<typename Float>
127  const Float& px, const Float& py,
128  const Float& qx, const Float& qy)
129 {
130  return sq<Float>(px-qx) + sq<Float>(py-qy);
131 }
132 
142 template<typename Float>
143 inline Float euclidean_distance(
144  const Float& px, const Float& py,
145  const Float& qx, const Float& qy)
146 {
147  return sqrt<Float>(euclidean_distance_sq<Float>(px,py,qx,qy));
148 }
149 
157 template<typename Float>
158 inline Float euclidean_distance(const Point& p, const Point& q)
159 {
160  return euclidean_distance<Float>(p.x(),p.y(),q.x(),q.y());
161 }
162 
170 template<typename Float>
171 inline Float euclidean_distance_sq(const Point& p, const Point& q)
172 {
173  return euclidean_distance_sq<Float>(p.x(),p.y(),q.x(),q.y());
174 }
175 
176 
188 template<typename Float>
190  const Float& ax, const Float& ay,
191  const Float& bx, const Float& by,
192  const Float& px, const Float& py)
193 {
194  Float n = (by-ay)*px - by*ax - (bx-ax)*py + bx*ay;
195  Float d = euclidean_distance(ax,ay, bx,by);
196  Q_ASSERT(d!=0.0);
197  return abs(n) / d;
198 }
199 
207 template<typename Float>
208 inline Float euclidean_segment_distance(const QLineF& s, const Point& p)
209 {
210  return euclidean_segment_distance<Float>(
211  s.p1().x(), s.p1().y(),
212  s.p2().x(), s.p2().y(),
213  p.x(),p.y());
214 }
215 
222 inline double min(double a, double b)
223 {
224  return (std::isnan(a) || (b < a)) ? b:a;
225 }
226 
233 inline double max(double a, double b)
234 {
235  return (std::isnan(a) || (b > a)) ? b:a;
236 }
243 inline int lower_floor(double x, int min=0) {
244  int i = (int)x;
245  if (i==x && i > min) --i;
246  return i;
247 }
254 inline int compare(double a, double b)
255 {
256  if (a > b) return +1;
257  if (a < b) return -1;
258  return 0;
259 }
266 {
267  double len = sqrt( sq(p.x()) + sq(p.y()) );
268  return Point(p.x()/len, p.y()/len);
269 }
278 inline Point intersection(Point p1, Point p2, Point p3, Point p4)
279 {
280  double det = (p1.x()-p2.x())*(p3.y()-p4.y()) - (p1.y()-p2.y())*(p3.x()-p4.x());
281  double Px = (p1.x()*p2.y()-p1.y()*p2.x())*(p3.x()-p4.x()) - (p1.x()-p2.x())*(p3.x()*p4.y()-p3.y()*p4.x());
282  double Py = (p1.x()*p2.y()-p1.y()*p2.x())*(p3.y()-p4.y()) - (p1.y()-p2.y())*(p3.x()*p4.y()-p3.y()*p4.x());
283 
284  return Point(Px/det,Py/det);
285 }
295 inline Point between(Point p, Point q, double w)
296 {
297  Q_ASSERT(w>=0 && w<=1.0);
298  return p*(1-w) + q*w;
299 }
306 inline double segment_length(const Curve& c, int i)
307 {
308  return euclidean_distance<double> (c[i],c[i+1]);
309 }
318 inline double max_euclidean_distance(const Curve& P, const Curve& Q)
319 {
320  double max=0.0;
321  for(int i=0; i < P.size(); ++i)
322  {
323  const Point& p = P[i];
324  for(int j=0; j < Q.size(); ++j)
325  {
326  double dsq = euclidean_distance_sq<double> (p,Q[j]);
327  if (dsq > max) max = dsq;
328  }
329  }
330  return sqrt(max);
331 }
332 
333 
334 
335 }} // namespace frechet::numeric
336 
337 #endif //FRECHET_VIEW_NUM_ERROR_H
double round_up(double x, int ulps=1)
rounding function with 'units in the last place' Roughly speaking, a ulp is the smallest representabl...
Definition: numeric.h:112
Point normalized(Point p)
normalize a vector to unit length
Definition: numeric.h:265
global definitions for all algorithms.
long double sqrt(const long double &x)
sqrt() specialisation for "long double"
Definition: numeric.h:82
int compare(double a, double b)
comparator function
Definition: numeric.h:254
long double accurate_float
Definition: numeric.h:42
int lower_floor(double x, int min=0)
floor function with lower limit
Definition: numeric.h:243
QPointF Point
a point in the plane; with double floating point precision. This type is heavily used throughout all ...
Definition: types.h:14
Point intersection(Point p1, Point p2, Point p3, Point p4)
compute the intersection of two lines
Definition: numeric.h:278
Point between(Point p, Point q, double w)
compute a point on a line segment
Definition: numeric.h:295
QPolygonF Curve
a polygonal curve in the plane; with double floating point precision. This type is heavily used throu...
Definition: types.h:20
double min(double a, double b)
minimum function with checks for NAN
Definition: numeric.h:222
double segment_length(const Curve &c, int i)
compute the length of a polygon segment
Definition: numeric.h:306
Float euclidean_distance(const Point &p, const Point &q)
euclidean distance function template for arbitrary floating point types
Definition: numeric.h:158
Float euclidean_distance_sq(const Point &p, const Point &q)
euclidean distance function template for arbitrary floating point types
Definition: numeric.h:171
Number sq(const Number &x)
square function template for arbitrary number types
Definition: numeric.h:61
double max_euclidean_distance(const Curve &P, const Curve &Q)
compute the maximum distance between two polygons This value serves as an upper bound for the Frechet...
Definition: numeric.h:318
Number abs(const Number &x)
abs() function template for arbitrary numerical types
Definition: numeric.h:91
static const constexpr double double_eps
smallest double value
Definition: numeric.h:52
Float euclidean_segment_distance(const QLineF &s, const Point &p)
Euclidean Distance between a point and a segment; function template for arbitrary floating point type...
Definition: numeric.h:208
double max(double a, double b)
maximum function with checks for NAN
Definition: numeric.h:233