geom.hh

   1 #ifndef GEOM_HH
   2 #define GEOM_HH
   3
   4 #include <math.h>
   5 #include <vector>
   6 #include <assert.h>
   7 #include "util.hh"
   8
   9 /** Classes for handling the geometry of points and lines. */
  10 namespace geom {
  11
  12   class LineRT; // Forward declaration for Line::Line(const LineRT &line)
  13
  14
  15
  16   //////////////////////////////////////////////////
  17   // Point
  18
  19
  20   /** A class representing a point in a 2-dimensional space. */
  21   struct Point {
  22
  23     /** Create a point. */
  24     Point() : x(0), y(0) { }
  25
  26     /** Create a point (\c x, \c y).
  27      * \param x = the x-coordinate
  28      * \param y = the y-coordinate
  29      */
  30     Point(float x, float y) : x(x), y(y) { }
  31
  32     /** Rotate the point 90 degrees. */
  33     Point& rot90()
  34     {
  35       *this = Point(y, -x);
  36       return *this;
  37     }
  38
  39     /** Return the point rotated 90 degrees.
  40      * \return The original point rotated 90 degrees.
  41      */
  42     Point get_rot90() const
  43     {
  44       return Point(y, -x);
  45     }
  46
  47     /** Compute dot product.
  48      * \param p = the point to compute the dot product with
  49      * \return the dot product
  50     */
  51     float dot(const Point &p) const { return x * p.x + y * p.y; }
  52
  53     /** Add another point to the point.
  54      * \param p = the point to add
  55      * \param scale = the scale for \c p
  56      */
  57     Point& add(const Point &p, float scale = 1)
  58     {
  59       x += scale * p.x;
  60       y += scale * p.y;
  61       return *this;
  62     }
  63
  64     /** Return the result of adding another point from the point.
  65      * \param p = the point to add
  66      * \param scale = the scale for \c p
  67      */
  68     Point get_add(const Point &p, float scale = 1) const
  69     {
  70       return Point(x + scale * p.x, y + scale * p.y);
  71     }
  72
  73     /** Substract another point from the point.
  74      * \param p = the point to substract
  75      * \param scale = the scale for \c p
  76      */
  77     Point& sub(const Point &p, float scale = 1)
  78     {
  79       x -= scale * p.x;
  80       y -= scale * p.y;
  81       return *this;
  82     }
  83
  84     /** Return the result of subtracting another point from the point.
  85      * \param p = the point to substract
  86      * \param scale = the scale for \c p
  87      */
  88     Point get_sub(const Point &p, float scale = 1) const
  89     {
  90       return Point(x - scale * p.x, y - scale * p.y);
  91     }
  92
  93     /** Scale the point.
  94      * \param scale = the scale
  95      */
  96     Point& scale(float scale)
  97     {
  98       x *= scale;
  99       y *= scale;
 100       return *this;
 101     }
 102
 103     /** Return scaled point.
 104      * \param scale = the scale
 105      */
 106     Point get_scale(float scale)
 107     {
 108       return Point(x * scale, y * scale);
 109     }
 110
 111     /** Compute the mean between two points.
 112      * \param point = the point to compute the mean with
 113      */
 114     Point& mean(Point point)
 115     {
 116       add(point).scale(0.5);
 117       return *this;
 118     }
 119
 120     /** Return the mean between two points.
 121      * \param point = the point to compute the mean with
 122      */
 123     Point get_mean(Point point)
 124     {
 125       return get_add(point).get_scale(0.5);
 126     }
 127
 128     /** Compute the angle of the vector in radians.
 129      *
 130      * \note The angle of (0,0) is also defined.  See the man page of
 131      * <EM>atan2</EM> for details.
 132      *
 133      * \return = the angle
 134      */
 135     float angle() const { return atan2(y, x); }
 136
 137     /** Compute the length of the vector.
 138      * \return the length
 139      */
 140     float length() const { return sqrt(x * x + y * y); }
 141
 142     /** Normalize the length of the vector.
 143      * \param length = the new length of the vector
 144      */
 145     Point& normalize(float length = 1)
 146     {
 147       float len = this->length();
 148       x /= len;
 149       y /= len;
 150       return *this;
 151     }
 152
 153     /** Coordinate of the point. */
 154     float x, y;
 155   };
 156
 157
 158
 159   //////////////////////////////////////////////////
 160   // Line
 161
 162
 163   /** A class representing a line as two end points. */
 164   struct Line {
 165
 166     /** Create a line. */
 167     Line() { }
 168
 169     /** Create a line.
 170      * \param a = the first end point
 171      * \param b = the second end point
 172      */
 173     Line(Point a, Point b) : a(a), b(b) { }
 174
 175     /** Create a line.
 176      * \param x1 = the first end point
 177      * \param y1 = the first end point
 178      * \param x2 = the second end point
 179      * \param y2 = the second end point
 180      */
 181     Line(float x1, float y1, float x2, float y2) : a(x1, y1), b(x2, y2) { }
 182
 183     /** Create a line by converting from the polar representation (rho, theta).
 184      *
 185      * \note The end points of the resulting line are guaranteed to
 186      * be on the line, but consider the actual positions to be
 187      * arbitrary.
 188      */
 189     Line(const LineRT &line); // Defined after the definition of LineRT
 190
 191     /** Add a point to the end points of the line.
 192      * \param p = the point to add to the line
 193      * \param scale = the scaling of the point before adding
 194      */
 195     void add(Point p, float scale = 1)
 196     {
 197       a.add(p, scale);
 198       b.add(p, scale);
 199     }
 200
 201     /** Compute the length of the line. */
 202     float length() const
 203     {
 204       return b.get_sub(a).length();
 205     }
 206
 207     /** Compute the normal of the line.
 208      * Note that the length of the normal is not normalized to one.
 209      * \return the normal of the line
 210      */
 211     Point normal() const
 212     {
 213       return b.get_sub(a).rot90();
 214     }
 215
 216     /** Compute the tangent of the line.
 217      * Note that the length of the tangent is not normalized to one.
 218      * \return the tangent of the line
 219      */
 220     Point tangent() const
 221     {
 222       return b.get_sub(a);
 223     }
 224
 225     /** Compute the intersection with \c line.
 226      * \param line = the line to intersect with
 227      * \return the intersection point
 228      */
 229     Point intersection(const Line &line) const
 230     {
 231       Point line_normal = line.normal();
 232       Point tangent = this->tangent();
 233       float k = line.a.get_sub(a).dot(line_normal) / tangent.dot(line_normal);
 234       return a.get_add(tangent, k);
 235     }
 236
 237     /** Cut the line between two lines.
 238      *
 239      * The first point will be the intersection with the first line,
 240      * and the second point will be the intersection with the second
 241      * line.
 242      *
 243      *  \param line1 = the first line
 244      *  \param line2 = the second line
 245      */
 246     Line& cut(const Line &line1, const Line &line2)
 247     {
 248       a = this->intersection(line1);
 249       b = this->intersection(line2);
 250       return *this;
 251     }
 252
 253     /** Compute the cut of the line between two lines.
 254      *
 255      * The first point will be the intersection with the first line,
 256      * and the second point will be the intersection with the second
 257      * line.
 258      *
 259      * \param line1 = the first line
 260      * \param line2 = the second line
 261      */
 262     Line get_cut(const Line &line1, const Line &line2) const
 263     {
 264       return Line(this->intersection(line1), this->intersection(line2));
 265     }
 266
 267     /** Compute the point on the line that is closest to \c point.
 268      * \param point = the point
 269      * \return the point closest to \c point
 270      */
 271     Point closest(const Point &point) const
 272     {
 273       Line line(point, point.get_add(this->normal()));
 274       return this->intersection(line);
 275     }
 276
 277     Point a; //!< The first end point.
 278     Point b; //!< The second end point.
 279
 280     /** Mirror a point with respect to the line.
 281      * \param point = point to mirror
 282      * \param scale = scale for the mirroring
 283      */
 284     Point mirror(const Point &point, float scale = 1) const
 285     {
 286       Point mirror_point = this->closest(point);
 287       return point.get_add(mirror_point.sub(point), 1 + scale);
 288     }
 289
 290     /** Mirror the line with respect to the given line.
 291      * \param mirror = the mirror line
 292      * \param scale = scale for the mirroring
 293      */
 294     Line& mirror(const Line &mirror, float scale = 1)
 295     {
 296       a = mirror.mirror(a, scale);
 297       b = mirror.mirror(b, scale);
 298       return *this;
 299     }
 300
 301     /** Return the line mirrored with respect to the given line.
 302      *  \param mirror = the mirror line
 303      */
 304     Line get_mirror(const Line &mirror) const
 305     {
 306       return Line(mirror.mirror(a), mirror.mirror(b));
 307     }
 308
 309
 310     /** Map the integer points of the line.
 311      *
 312      * The method calls \c obj(x, y) for each integer point of the
 313      * line.
 314      *
 315      * \bug The float coordinates should be handled more elegantly so
 316      * that the line is processed along the longer axis and the other
 317      * coordinate is computed for each middle location.
 318      *
 319      * \param obj = the object to call for each line point
 320      * \return a reference to the object \c obj
 321      */
 322     template <typename T>
 323     T&
 324     map(T &obj)
 325     {
 326       // Check if we have just a point, and process it
 327       if ((int)(a.x + 0.5) == (int)(b.x + 0.5) &&
 328           (int)(a.y + 0.5) == (int)(b.y + 0.5))
 329       {
 330         obj((int)(a.x + 0.5), (int)(b.y + 0.5));
 331         return obj;
 332       }
 333
 334       // Compute the line
 335       float dx = b.x - a.x;
 336       float dy = b.y - a.y;
 337       float div = util::max(util::abs(dx), util::abs(dy));
 338       dx /= div;
 339       dy /= div;
 340
 341       // Process the points of the line
 342       for (int i = 0; i <= div; i++) {
 343         int x = (int)(a.x + i * dx + 0.5);
 344         int y = (int)(a.y + i * dy + 0.5);
 345         obj(x, y);
 346       }
 347
 348       return obj;
 349     }
 350
 351   };
 352
 353
 354
 355   //////////////////////////////////////////////////
 356   // LineRT
 357
 358
 359   /** A class representing a line as the angle and the distance from
 360    * the origin.
 361    *
 362    * In this representation, \c rho is the distance between the line
 363    * and the origin, and \c theta is the angle of the <EM>normal</EM>
 364    * of the line.
 365    *
 366    * \note \c rho can also be negative, and the range of \c theta is
 367    * not restricted.
 368    */
 369   struct LineRT {
 370
 371     /** Create a line. */
 372     LineRT() : rho(0), theta(0) { }
 373
 374     /** Create a line.
 375      * \param rho = the distance to the origin
 376      * \param theta = the angle of the normal
 377      */
 378     LineRT(float rho, float theta) : rho(rho), theta(theta) { }
 379
 380     /** Create a line by converting from an end-point representation.
 381      * \param line = the line to convert from
 382      */
 383     LineRT(const Line &line)
 384     {
 385       Point p = line.closest(Point(0, 0));
 386       rho = p.length();
 387       theta = p.angle();
 388     }
 389
 390     /** Compares lines according to the rho values. */
 391     bool operator<(const LineRT &line) const
 392     {
 393       if (rho < line.rho)
 394         return true;
 395       return false;
 396     }
 397
 398     float rho;   //!< the distance to the origin
 399     float theta; //!< the angle of the normal
 400   };
 401
 402
 403
 404   /** A structure representing a grid as a matrix of points. */
 405   struct Grid {
 406
 407     /** Create a grid of certain size.
 408      * \param width = the width of the grid
 409      * \param height = the height of the grid
 410      */
 411     Grid(int width, int height) : width(width), height(height),
 412                                   points(width * height) { }
 413
 414     /** Create a grid of the intersections between two series of
 415      * lines.
 416      *
 417      * This constructor takes two sets of lines, and for each line in
 418      * the first set, it computes the intersections with each line
 419      * from the second set.  The intersections are stored so that n'th
 420      * row contains the intersections of the n'th line of the first
 421      * set.
 422      *
 423      * \param lines1 = the first set of lines
 424      * \param lines2 = the second set of lines
 425      */
 426     Grid(const std::vector<Line> &lines1, const std::vector<Line> &lines2)
 427       : width(lines2.size()), height(lines1.size()), points(width * height)
 428     {
 429       for (int l1 = 0; l1 < (int)lines1.size(); l1++)
 430         for (int l2 = 0; l2 < (int)lines2.size(); l2++)
 431           (*this)(l2, l1) = lines1[l1].intersection(lines2[l2]);
 432     }
 433
 434     /** Access a point.
 435      *
 436      * \note The function checks the bounds.
 437      *
 438      * \param x = the column of the grid
 439      * \param y = the row of the grid
 440      */
 441     Point &operator()(int x, int y)
 442     {
 443       assert(x >= 0 && x < width);
 444       assert(y >= 0 && y < height);
 445       return points[y * width + x];
 446     }
 447
 448     /** Access the point buffer directly.
 449      *
 450      * \note The function checks the bounds.
 451      *
 452      * \param index = the index of the point in the buffer
 453      */
 454     Point &operator[](int index)
 455     {
 456       assert(index >= 0 && index < (int)points.size());
 457       return points[index];
 458     }
 459
 460
 461     int width; //!< The width of the grid.
 462     int height; //!< The height of the grid.
 463
 464     /** The points in the grid: 1st row, 2nd row, and so on. */
 465     std::vector<Point> points;
 466   };
 467
 468   /** Compute a weighted mean of two points.
 469    * \param p1 = the first points
 470    * \param p2 = the second point
 471    * \param weight = the weight of the first point
 472    */
 473   inline
 474   Point
 475   mean(Point p1, Point p2, float weight = 0.5)
 476   {
 477     p2.sub(p1);
 478     return p1.add(p2, 1.0 - weight);
 479   }
 480
 481   //////////////////////////////////////////////////
 482   // Conversion functions
 483
 484   // Documented above
 485   inline
 486   Line::Line(const LineRT &line)
 487   {
 488     Point normal(cos(line.theta), sin(line.theta));
 489     a = normal.get_scale(line.rho);
 490     b = a.get_add(normal.get_rot90());
 491   }
 492
 493   /** Convert degrees to radians.
 494    * \param angle = the angle in degrees
 495    * \return the angle in radians
 496    */
 497   inline
 498   float rad(float angle) { return angle / 180 * M_PI; }
 499
 500   /** Convert radians to degrees.
 501    * \param angle = the angle in radians
 502    * \return the angle in degrees
 503    */
 504   inline
 505   float deg(float angle) { return angle * 180 / M_PI; }
 506 }
 507
 508 #endif /* GEOM_HH */