1 /** @file 2 * @brief Cartesian point / 2D vector and related operations 3 *//* 4 * Authors: 5 * Michael G. Sloan <mgsloan@gmail.com> 6 * Nathan Hurst <njh@njhurst.com> 7 * Krzysztof Kosiński <tweenk.pl@gmail.com> 8 * 9 * Copyright (C) 2006-2009 Authors 10 * 11 * This library is free software; you can redistribute it and/or 12 * modify it either under the terms of the GNU Lesser General Public 13 * License version 2.1 as published by the Free Software Foundation 14 * (the "LGPL") or, at your option, under the terms of the Mozilla 15 * Public License Version 1.1 (the "MPL"). If you do not alter this 16 * notice, a recipient may use your version of this file under either 17 * the MPL or the LGPL. 18 * 19 * You should have received a copy of the LGPL along with this library 20 * in the file COPYING-LGPL-2.1; if not, write to the Free Software 21 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 22 * You should have received a copy of the MPL along with this library 23 * in the file COPYING-MPL-1.1 24 * 25 * The contents of this file are subject to the Mozilla Public License 26 * Version 1.1 (the "License"); you may not use this file except in 27 * compliance with the License. You may obtain a copy of the License at 28 * http://www.mozilla.org/MPL/ 29 * 30 * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY 31 * OF ANY KIND, either express or implied. See the LGPL or the MPL for 32 * the specific language governing rights and limitations. 33 */ 34 35 #ifndef LIB2GEOM_SEEN_POINT_H 36 #define LIB2GEOM_SEEN_POINT_H 37 38 #include <iostream> 39 #include <iterator> 40 #include <boost/operators.hpp> 41 #include <2geom/forward.h> 42 #include <2geom/coord.h> 43 #include <2geom/int-point.h> 44 #include <2geom/math-utils.h> 45 #include <2geom/utils.h> 46 47 namespace Geom { 48 49 class Point 50 : boost::additive< Point 51 , boost::totally_ordered< Point 52 , boost::multiplicative< Point, Coord 53 , MultipliableNoncommutative< Point, Affine 54 , MultipliableNoncommutative< Point, Translate 55 , MultipliableNoncommutative< Point, Rotate 56 , MultipliableNoncommutative< Point, Scale 57 , MultipliableNoncommutative< Point, HShear 58 , MultipliableNoncommutative< Point, VShear 59 , MultipliableNoncommutative< Point, Zoom 60 > > > > > > > > > > // base class chaining, see documentation for Boost.Operator 61 { 62 Coord _pt[2]; 63 public: 64 typedef Coord D1Value; 65 typedef Coord &D1Reference; 66 typedef Coord const &D1ConstReference; 67 68 /// @name Create points 69 /// @{ 70 /** Construct a point on the origin. */ Point()71 Point() 72 { _pt[X] = _pt[Y] = 0; } 73 74 /** Construct a point from its coordinates. */ Point(Coord x,Coord y)75 Point(Coord x, Coord y) { 76 _pt[X] = x; _pt[Y] = y; 77 } 78 /** Construct from integer point. */ Point(IntPoint const & p)79 Point(IntPoint const &p) { 80 _pt[X] = p[X]; 81 _pt[Y] = p[Y]; 82 } 83 /** @brief Construct a point from its polar coordinates. 84 * The angle is specified in radians, in the mathematical convention (increasing 85 * counter-clockwise from +X). */ polar(Coord angle,Coord radius)86 static Point polar(Coord angle, Coord radius) { 87 Point ret(polar(angle)); 88 ret *= radius; 89 return ret; 90 } 91 /** @brief Construct an unit vector from its angle. 92 * The angle is specified in radians, in the mathematical convention (increasing 93 * counter-clockwise from +X). */ 94 static Point polar(Coord angle); 95 /// @} 96 97 /// @name Access the coordinates of a point 98 /// @{ 99 Coord operator[](unsigned i) const { return _pt[i]; } 100 Coord &operator[](unsigned i) { return _pt[i]; } 101 102 Coord operator[](Dim2 d) const noexcept { return _pt[d]; } 103 Coord &operator[](Dim2 d) noexcept { return _pt[d]; } 104 x()105 Coord x() const noexcept { return _pt[X]; } x()106 Coord &x() noexcept { return _pt[X]; } y()107 Coord y() const noexcept { return _pt[Y]; } y()108 Coord &y() noexcept { return _pt[Y]; } 109 /// @} 110 111 /// @name Vector operations 112 /// @{ 113 /** @brief Compute the distance from origin. 114 * @return Length of the vector from origin to this point */ length()115 Coord length() const { return hypot(_pt[0], _pt[1]); } 116 void normalize(); normalized()117 Point normalized() const { 118 Point ret(*this); 119 ret.normalize(); 120 return ret; 121 } 122 123 /** @brief Return a point like this point but rotated -90 degrees. 124 * If the y axis grows downwards and the x axis grows to the 125 * right, then this is 90 degrees counter-clockwise. */ ccw()126 Point ccw() const { 127 return Point(_pt[Y], -_pt[X]); 128 } 129 130 /** @brief Return a point like this point but rotated +90 degrees. 131 * If the y axis grows downwards and the x axis grows to the 132 * right, then this is 90 degrees clockwise. */ cw()133 Point cw() const { 134 return Point(-_pt[Y], _pt[X]); 135 } 136 /// @} 137 138 /// @name Vector-like arithmetic operations 139 /// @{ 140 Point operator-() const { 141 return Point(-_pt[X], -_pt[Y]); 142 } 143 Point &operator+=(Point const &o) { 144 for ( unsigned i = 0 ; i < 2 ; ++i ) { 145 _pt[i] += o._pt[i]; 146 } 147 return *this; 148 } 149 Point &operator-=(Point const &o) { 150 for ( unsigned i = 0 ; i < 2 ; ++i ) { 151 _pt[i] -= o._pt[i]; 152 } 153 return *this; 154 } 155 Point &operator*=(Coord s) { 156 for (double & i : _pt) i *= s; 157 return *this; 158 } 159 Point &operator/=(Coord s) { 160 //TODO: s == 0? 161 for (double & i : _pt) i /= s; 162 return *this; 163 } 164 /// @} 165 166 /// @name Affine transformations 167 /// @{ 168 Point &operator*=(Affine const &m); 169 // implemented in transforms.cpp 170 Point &operator*=(Translate const &t); 171 Point &operator*=(Scale const &s); 172 Point &operator*=(Rotate const &r); 173 Point &operator*=(HShear const &s); 174 Point &operator*=(VShear const &s); 175 Point &operator*=(Zoom const &z); 176 /// @} 177 178 /// @name Conversion to integer points 179 /// @{ 180 /** @brief Round to nearest integer coordinates. */ round()181 IntPoint round() const { 182 IntPoint ret(::round(_pt[X]), ::round(_pt[Y])); 183 return ret; 184 } 185 /** @brief Round coordinates downwards. */ floor()186 IntPoint floor() const { 187 IntPoint ret(::floor(_pt[X]), ::floor(_pt[Y])); 188 return ret; 189 } 190 /** @brief Round coordinates upwards. */ ceil()191 IntPoint ceil() const { 192 IntPoint ret(::ceil(_pt[X]), ::ceil(_pt[Y])); 193 return ret; 194 } 195 /// @} 196 197 /// @name Various utilities 198 /// @{ 199 /** @brief Check whether both coordinates are finite. */ isFinite()200 bool isFinite() const { 201 for (double i : _pt) { 202 if(!std::isfinite(i)) return false; 203 } 204 return true; 205 } 206 /** @brief Check whether both coordinates are zero. */ isZero()207 bool isZero() const { 208 return _pt[X] == 0 && _pt[Y] == 0; 209 } 210 /** @brief Check whether the length of the vector is close to 1. */ 211 bool isNormalized(Coord eps=EPSILON) const { 212 return are_near(length(), 1.0, eps); 213 } 214 /** @brief Equality operator. 215 * This tests for exact identity (as opposed to are_near()). Note that due to numerical 216 * errors, this test might return false even if the points should be identical. */ 217 bool operator==(const Point &in_pnt) const { 218 return (_pt[X] == in_pnt[X]) && (_pt[Y] == in_pnt[Y]); 219 } 220 /** @brief Lexicographical ordering for points. 221 * Y coordinate is regarded as more significant. When sorting according to this 222 * ordering, the points will be sorted according to the Y coordinate, and within 223 * points with the same Y coordinate according to the X coordinate. */ 224 bool operator<(const Point &p) const { 225 return _pt[Y] < p[Y] || (_pt[Y] == p[Y] && _pt[X] < p[X]); 226 } 227 /// @} 228 229 /** @brief Lexicographical ordering functor. 230 * @param d The dimension with higher significance */ 231 template <Dim2 DIM> struct LexLess; 232 template <Dim2 DIM> struct LexGreater; 233 //template <Dim2 DIM, typename First = std::less<Coord>, typename Second = std::less<Coord> > LexOrder; 234 /** @brief Lexicographical ordering functor with runtime dimension. */ 235 struct LexLessRt { LexLessRtLexLessRt236 LexLessRt(Dim2 d) : dim(d) {} 237 inline bool operator()(Point const &a, Point const &b) const; 238 private: 239 Dim2 dim; 240 }; 241 struct LexGreaterRt { LexGreaterRtLexGreaterRt242 LexGreaterRt(Dim2 d) : dim(d) {} 243 inline bool operator()(Point const &a, Point const &b) const; 244 private: 245 Dim2 dim; 246 }; 247 //template <typename First = std::less<Coord>, typename Second = std::less<Coord> > LexOrder 248 }; 249 250 /** @brief Output operator for points. 251 * Prints out the coordinates. 252 * @relates Point */ 253 std::ostream &operator<<(std::ostream &out, const Geom::Point &p); 254 255 template<> struct Point::LexLess<X> { 256 typedef std::less<Coord> Primary; 257 typedef std::less<Coord> Secondary; 258 typedef std::less<Coord> XOrder; 259 typedef std::less<Coord> YOrder; 260 bool operator()(Point const &a, Point const &b) const { 261 return a[X] < b[X] || (a[X] == b[X] && a[Y] < b[Y]); 262 } 263 }; 264 template<> struct Point::LexLess<Y> { 265 typedef std::less<Coord> Primary; 266 typedef std::less<Coord> Secondary; 267 typedef std::less<Coord> XOrder; 268 typedef std::less<Coord> YOrder; 269 bool operator()(Point const &a, Point const &b) const { 270 return a[Y] < b[Y] || (a[Y] == b[Y] && a[X] < b[X]); 271 } 272 }; 273 template<> struct Point::LexGreater<X> { 274 typedef std::greater<Coord> Primary; 275 typedef std::greater<Coord> Secondary; 276 typedef std::greater<Coord> XOrder; 277 typedef std::greater<Coord> YOrder; 278 bool operator()(Point const &a, Point const &b) const { 279 return a[X] > b[X] || (a[X] == b[X] && a[Y] > b[Y]); 280 } 281 }; 282 template<> struct Point::LexGreater<Y> { 283 typedef std::greater<Coord> Primary; 284 typedef std::greater<Coord> Secondary; 285 typedef std::greater<Coord> XOrder; 286 typedef std::greater<Coord> YOrder; 287 bool operator()(Point const &a, Point const &b) const { 288 return a[Y] > b[Y] || (a[Y] == b[Y] && a[X] > b[X]); 289 } 290 }; 291 inline bool Point::LexLessRt::operator()(Point const &a, Point const &b) const { 292 return dim ? Point::LexLess<Y>()(a, b) : Point::LexLess<X>()(a, b); 293 } 294 inline bool Point::LexGreaterRt::operator()(Point const &a, Point const &b) const { 295 return dim ? Point::LexGreater<Y>()(a, b) : Point::LexGreater<X>()(a, b); 296 } 297 298 /** @brief Compute the second (Euclidean) norm of @a p. 299 * This corresponds to the length of @a p. The result will not overflow even if 300 * \f$p_X^2 + p_Y^2\f$ is larger that the maximum value that can be stored 301 * in a <code>double</code>. 302 * @return \f$\sqrt{p_X^2 + p_Y^2}\f$ 303 * @relates Point */ 304 inline Coord L2(Point const &p) { 305 return p.length(); 306 } 307 308 /** @brief Compute the square of the Euclidean norm of @a p. 309 * Warning: this can overflow where L2 won't. 310 * @return \f$p_X^2 + p_Y^2\f$ 311 * @relates Point */ 312 inline Coord L2sq(Point const &p) { 313 return p[0]*p[0] + p[1]*p[1]; 314 } 315 316 /** @brief Returns p * Geom::rotate_degrees(90), but more efficient. 317 * 318 * Angle direction in 2Geom: If you use the traditional mathematics convention that y 319 * increases upwards, then positive angles are anticlockwise as per the mathematics convention. If 320 * you take the common non-mathematical convention that y increases downwards, then positive angles 321 * are clockwise, as is common outside of mathematics. 322 * 323 * There is no function to rotate by -90 degrees: use -rot90(p) instead. 324 * @relates Point */ 325 inline Point rot90(Point const &p) { 326 return Point(-p[Y], p[X]); 327 } 328 329 /** @brief Linear interpolation between two points. 330 * @param t Time value 331 * @param a First point 332 * @param b Second point 333 * @return Point on a line between a and b. The ratio of its distance from a 334 * and the distance between a and b will be equal to t. 335 * @relates Point */ 336 inline Point lerp(Coord t, Point const &a, Point const &b) { 337 return (1 - t) * a + t * b; 338 } 339 340 /** @brief Return a point halfway between the specified ones. 341 * @relates Point */ 342 inline Point middle_point(Point const &p1, Point const &p2) { 343 return lerp(0.5, p1, p2); 344 } 345 346 /** @brief Compute the dot product of a and b. 347 * Dot product can be interpreted as a measure of how parallel the vectors are. 348 * For perpendicular vectors, it is zero. For parallel ones, its absolute value is highest, 349 * and the sign depends on whether they point in the same direction (+) or opposite ones (-). 350 * @return \f$a \cdot b = a_X b_X + a_Y b_Y\f$. 351 * @relates Point */ 352 inline Coord dot(Point const &a, Point const &b) { 353 return a[X] * b[X] + a[Y] * b[Y]; 354 } 355 356 /** @brief Compute the 2D cross product. 357 * This is also known as "perp dot product". It will be zero for parallel vectors, 358 * and the absolute value will be highest for perpendicular vectors. 359 * @return \f$a \times b = a_X b_Y - a_Y b_X\f$. 360 * @relates Point*/ 361 inline Coord cross(Point const &a, Point const &b) 362 { 363 // equivalent implementation: 364 // return dot(a, b.ccw()); 365 return a[X] * b[Y] - a[Y] * b[X]; 366 } 367 368 /// Compute the (Euclidean) distance between points. 369 /// @relates Point 370 inline Coord distance (Point const &a, Point const &b) { 371 return (a - b).length(); 372 } 373 374 /// Compute the square of the distance between points. 375 /// @relates Point 376 inline Coord distanceSq (Point const &a, Point const &b) { 377 return L2sq(a - b); 378 } 379 380 //IMPL: NearConcept 381 /// Test whether two points are no further apart than some threshold. 382 /// @relates Point 383 inline bool are_near(Point const &a, Point const &b, double eps = EPSILON) { 384 // do not use an unqualified calls to distance before the empty 385 // specialization of iterator_traits is defined - see end of file 386 return are_near((a - b).length(), 0, eps); 387 } 388 389 /// Test whether the relative distance between two points is less than some threshold. 390 inline bool are_near_rel(Point const &a, Point const &b, double eps = EPSILON) { 391 return (a - b).length() <= eps * (a.length() + b.length()) / 2; 392 } 393 394 /// Test whether three points lie approximately on the same line. 395 /// @relates Point 396 inline bool are_collinear(Point const& p1, Point const& p2, Point const& p3, 397 double eps = EPSILON) 398 { 399 return are_near( cross(p3, p2) - cross(p3, p1) + cross(p2, p1), 0, eps); 400 } 401 402 Point unit_vector(Point const &a); 403 Coord L1(Point const &p); 404 Coord LInfty(Point const &p); 405 bool is_zero(Point const &p); 406 bool is_unit_vector(Point const &p, Coord eps = EPSILON); 407 double atan2(Point const &p); 408 double angle_between(Point const &a, Point const &b); 409 Point abs(Point const &b); 410 Point constrain_angle(Point const &A, Point const &B, unsigned int n = 4, Geom::Point const &dir = Geom::Point(1,0)); 411 412 } // end namespace Geom 413 414 // This is required to fix a bug in GCC 4.3.3 (and probably others) that causes the compiler 415 // to try to instantiate the iterator_traits template and fail. Probably it thinks that Point 416 // is an iterator and tries to use std::distance instead of Geom::distance. 417 namespace std { 418 template <> class iterator_traits<Geom::Point> {}; 419 } 420 421 #endif // LIB2GEOM_SEEN_POINT_H 422 423 /* 424 Local Variables: 425 mode:c++ 426 c-file-style:"stroustrup" 427 c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +)) 428 indent-tabs-mode:nil 429 fill-column:99 430 End: 431 */ 432 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 : 433