|
1
by Jasper van de Gronde
Initial commit |
1 |
#include <2geom/point.h> |
2 |
#include <assert.h> |
|
3 |
#include <2geom/coord.h> |
|
4 |
#include <2geom/isnan.h> //temporary fix for isnan() |
|
5 |
#include <2geom/transforms.h> |
|
6 |
||
7 |
namespace Geom { |
|
8 |
||
9 |
/** \brief Scales this vector to make it a unit vector (within rounding error).
|
|
10 |
*
|
|
11 |
* The current version tries to handle infinite coordinates gracefully,
|
|
12 |
* but it's not clear that any callers need that.
|
|
13 |
*
|
|
14 |
* \pre \f$this \neq (0, 0)\f$
|
|
15 |
* \pre Neither component is NaN.
|
|
16 |
* \post \f$-\epsilon<\left|this\right|-1<\epsilon\f$
|
|
17 |
*/
|
|
18 |
void Point::normalize() { |
|
19 |
double len = hypot(_pt[0], _pt[1]); |
|
20 |
if(len == 0) return; |
|
21 |
if(IS_NAN(len)) return; |
|
22 |
static double const inf = 1e400; |
|
23 |
if(len != inf) { |
|
24 |
*this /= len; |
|
25 |
} else { |
|
26 |
unsigned n_inf_coords = 0; |
|
27 |
/* Delay updating pt in case neither coord is infinite. */
|
|
28 |
Point tmp; |
|
29 |
for ( unsigned i = 0 ; i < 2 ; ++i ) { |
|
30 |
if ( _pt[i] == inf ) { |
|
31 |
++n_inf_coords; |
|
32 |
tmp[i] = 1.0; |
|
33 |
} else if ( _pt[i] == -inf ) { |
|
34 |
++n_inf_coords; |
|
35 |
tmp[i] = -1.0; |
|
36 |
} else { |
|
37 |
tmp[i] = 0.0; |
|
38 |
}
|
|
39 |
}
|
|
40 |
switch (n_inf_coords) { |
|
41 |
case 0: { |
|
42 |
/* Can happen if both coords are near +/-DBL_MAX. */
|
|
43 |
*this /= 4.0; |
|
44 |
len = hypot(_pt[0], _pt[1]); |
|
45 |
assert(len != inf); |
|
46 |
*this /= len; |
|
47 |
break; |
|
48 |
}
|
|
49 |
case 1: { |
|
50 |
*this = tmp; |
|
51 |
break; |
|
52 |
}
|
|
53 |
case 2: { |
|
54 |
*this = tmp * sqrt(0.5); |
|
55 |
break; |
|
56 |
}
|
|
57 |
}
|
|
58 |
}
|
|
59 |
}
|
|
60 |
||
61 |
/** Compute the L1 norm, or manhattan distance, of \a p. */
|
|
62 |
Coord L1(Point const &p) { |
|
63 |
Coord d = 0; |
|
64 |
for ( int i = 0 ; i < 2 ; i++ ) { |
|
65 |
d += fabs(p[i]); |
|
66 |
}
|
|
67 |
return d; |
|
68 |
}
|
|
69 |
||
70 |
/** Compute the L infinity, or maximum, norm of \a p. */
|
|
71 |
Coord LInfty(Point const &p) { |
|
72 |
Coord const a(fabs(p[0])); |
|
73 |
Coord const b(fabs(p[1])); |
|
74 |
return ( a < b || IS_NAN(b) |
|
75 |
? b |
|
76 |
: a ); |
|
77 |
}
|
|
78 |
||
79 |
/** Returns true iff p is a zero vector, i.e.\ Point(0, 0).
|
|
80 |
*
|
|
81 |
* (NaN is considered non-zero.)
|
|
82 |
*/
|
|
83 |
bool
|
|
84 |
is_zero(Point const &p) |
|
85 |
{
|
|
86 |
return ( p[0] == 0 && |
|
87 |
p[1] == 0 ); |
|
88 |
}
|
|
89 |
||
90 |
bool
|
|
91 |
is_unit_vector(Point const &p) |
|
92 |
{
|
|
93 |
return fabs(1.0 - L2(p)) <= 1e-4; |
|
94 |
/* The tolerance of 1e-4 is somewhat arbitrary. Point::normalize is believed to return
|
|
95 |
points well within this tolerance. I'm not aware of any callers that want a small
|
|
96 |
tolerance; most callers would be ok with a tolerance of 0.25. */
|
|
97 |
}
|
|
98 |
||
99 |
Coord atan2(Point const p) { |
|
100 |
return std::atan2(p[Y], p[X]); |
|
101 |
}
|
|
102 |
||
103 |
/** compute the angle turning from a to b. This should give \f$\pi/2\f$ for angle_between(a, rot90(a));
|
|
104 |
* This works by projecting b onto the basis defined by a, rot90(a)
|
|
105 |
*/
|
|
106 |
Coord angle_between(Point const a, Point const b) { |
|
107 |
return std::atan2(cross(b,a), dot(b,a)); |
|
108 |
}
|
|
109 |
||
110 |
||
111 |
||
112 |
/** Returns a version of \a a scaled to be a unit vector (within rounding error).
|
|
113 |
*
|
|
114 |
* The current version tries to handle infinite coordinates gracefully,
|
|
115 |
* but it's not clear that any callers need that.
|
|
116 |
*
|
|
117 |
* \pre a != Point(0, 0).
|
|
118 |
* \pre Neither coordinate is NaN.
|
|
119 |
* \post L2(ret) very near 1.0.
|
|
120 |
*/
|
|
121 |
Point unit_vector(Point const &a) |
|
122 |
{
|
|
123 |
Point ret(a); |
|
124 |
ret.normalize(); |
|
125 |
return ret; |
|
126 |
}
|
|
127 |
||
128 |
Point abs(Point const &b) |
|
129 |
{
|
|
130 |
Point ret; |
|
131 |
for ( int i = 0 ; i < 2 ; i++ ) { |
|
132 |
ret[i] = fabs(b[i]); |
|
133 |
}
|
|
134 |
return ret; |
|
135 |
}
|
|
136 |
||
137 |
Point operator*(Point const &v, Matrix const &m) { |
|
138 |
Point ret; |
|
139 |
for(int i = 0; i < 2; i++) { |
|
140 |
ret[i] = v[X] * m[i] + v[Y] * m[i + 2] + m[i + 4]; |
|
141 |
}
|
|
142 |
return ret; |
|
143 |
}
|
|
144 |
||
145 |
Point operator/(Point const &p, Matrix const &m) { return p * m.inverse(); } |
|
146 |
||
147 |
Point &Point::operator*=(Matrix const &m) |
|
148 |
{
|
|
149 |
*this = *this * m; |
|
150 |
return *this; |
|
151 |
}
|
|
152 |
||
153 |
Point constrain_angle(Point const &A, Point const &B, unsigned int n, Point const &dir) |
|
154 |
{
|
|
155 |
// for special cases we could perhaps use explicit testing (which might be faster)
|
|
156 |
if (n == 0.0) { |
|
157 |
return B; |
|
158 |
}
|
|
159 |
Point diff(B - A); |
|
160 |
double angle = -angle_between(diff, dir); |
|
161 |
double k = round(angle * (double)n / (2.0*M_PI)); |
|
162 |
return A + dir * Rotate(k * 2.0 * M_PI / (double)n) * L2(diff); |
|
163 |
}
|
|
164 |
||
165 |
} //namespace Geom |
|
166 |
||
167 |
/*
|
|
168 |
Local Variables:
|
|
169 |
mode:c++
|
|
170 |
c-file-style:"stroustrup"
|
|
171 |
c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
|
|
172 |
indent-tabs-mode:nil
|
|
173 |
fill-column:99
|
|
174 |
End:
|
|
175 |
*/
|
|
176 |
// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :
|