src/2geom/geom.cpp

   1 /**
   2  *  \file geom.cpp
   3  *  \brief Various geometrical calculations.
   4  */
   5
   6 #ifdef HAVE_CONFIG_H
   7 # include <config.h>
   8 #endif
   9 #include <2geom/geom.h>
  10 #include <2geom/point.h>
  11 #include <algorithm>
  12
  13 namespace Geom {
  14
  15 /**
  16  * Finds the intersection of the two (infinite) lines
  17  * defined by the points p such that dot(n0, p) == d0 and dot(n1, p) == d1.
  18  *
  19  * If the two lines intersect, then \a result becomes their point of
  20  * intersection; otherwise, \a result remains unchanged.
  21  *
  22  * This function finds the intersection of the two lines (infinite)
  23  * defined by n0.X = d0 and x1.X = d1.  The algorithm is as follows:
  24  * To compute the intersection point use kramer's rule:
  25  * \verbatim
  26  * convert lines to form
  27  * ax + by = c
  28  * dx + ey = f
  29  *
  30  * (
  31  *  e.g. a = (x2 - x1), b = (y2 - y1), c = (x2 - x1)*x1 + (y2 - y1)*y1
  32  * )
  33  *
  34  * In our case we use:
  35  *   a = n0.x     d = n1.x
  36  *   b = n0.y     e = n1.y
  37  *   c = d0        f = d1
  38  *
  39  * so:
  40  *
  41  * adx + bdy = cd
  42  * adx + aey = af
  43  *
  44  * bdy - aey = cd - af
  45  * (bd - ae)y = cd - af
  46  *
  47  * y = (cd - af)/(bd - ae)
  48  *
  49  * repeat for x and you get:
  50  *
  51  * x = (fb - ce)/(bd - ae)                \endverbatim
  52  *
  53  * If the denominator (bd-ae) is 0 then the lines are parallel, if the
  54  * numerators are 0 then the lines coincide.
  55  *
  56  * \todo Why not use existing but outcommented code below
  57  * (HAVE_NEW_INTERSECTOR_CODE)?
  58  */
  59 IntersectorKind
  60 line_intersection(Geom::Point const &n0, double const d0,
  61                   Geom::Point const &n1, double const d1,
  62                   Geom::Point &result)
  63 {
  64     double denominator = dot(Geom::rot90(n0), n1);
  65     double X = n1[Geom::Y] * d0 -
  66         n0[Geom::Y] * d1;
  67     /* X = (-d1, d0) dot (n0[Y], n1[Y]) */
  68
  69     if (denominator == 0) {
  70         if ( X == 0 ) {
  71             return coincident;
  72         } else {
  73             return parallel;
  74         }
  75     }
  76
  77     double Y = n0[Geom::X] * d1 -
  78         n1[Geom::X] * d0;
  79
  80     result = Geom::Point(X, Y) / denominator;
  81
  82     return intersects;
  83 }
  84
  85
  86
  87
  88 /* ccw exists as a building block */
  89 int
  90 intersector_ccw(const Geom::Point& p0, const Geom::Point& p1,
  91         const Geom::Point& p2)
  92 /* Determine which way a set of three points winds. */
  93 {
  94     Geom::Point d1 = p1 - p0;
  95     Geom::Point d2 = p2 - p0;
  96     /* compare slopes but avoid division operation */
  97     double c = dot(Geom::rot90(d1), d2);
  98     if(c > 0)
  99         return +1; // ccw - do these match def'n in header?
 100     if(c < 0)
 101         return -1; // cw
 102
 103     /* Colinear [or NaN].  Decide the order. */
 104     if ( ( d1[0] * d2[0] < 0 )  ||
 105          ( d1[1] * d2[1] < 0 ) ) {
 106         return -1; // p2  <  p0 < p1
 107     } else if ( dot(d1,d1) < dot(d2,d2) ) {
 108         return +1; // p0 <= p1  <  p2
 109     } else {
 110         return 0; // p0 <= p2 <= p1
 111     }
 112 }
 113
 114 /** Determine whether the line segment from p00 to p01 intersects the
 115     infinite line passing through p10 and p11. This doesn't find the
 116     point of intersection, use the line_intersect function above,
 117     or the segment_intersection interface below.
 118
 119     \pre neither segment is zero-length; i.e. p00 != p01 and p10 != p11.
 120 */
 121 bool
 122 line_segment_intersectp(Geom::Point const &p00, Geom::Point const &p01,
 123                         Geom::Point const &p10, Geom::Point const &p11)
 124 {
 125     if(p00 == p01) return false;
 126     if(p10 == p11) return false;
 127
 128     return ((intersector_ccw(p00, p01, p10) * intersector_ccw(p00, p01, p11)) <= 0 );
 129 }
 130
 131
 132 /** Determine whether two line segments intersect.  This doesn't find
 133     the point of intersection, use the line_intersect function above,
 134     or the segment_intersection interface below.
 135
 136     \pre neither segment is zero-length; i.e. p00 != p01 and p10 != p11.
 137 */
 138 bool
 139 segment_intersectp(Geom::Point const &p00, Geom::Point const &p01,
 140                    Geom::Point const &p10, Geom::Point const &p11)
 141 {
 142     if(p00 == p01) return false;
 143     if(p10 == p11) return false;
 144
 145     /* true iff (    (the p1 segment straddles the p0 infinite line)
 146      *           and (the p0 segment straddles the p1 infinite line) ). */
 147     return (line_segment_intersectp(p00, p01, p10, p11) &&
 148             line_segment_intersectp(p10, p11, p00, p01));
 149 }
 150
 151 /** Determine whether \& where a line segments intersects an (infinite) line.
 152
 153 If there is no intersection, then \a result remains unchanged.
 154
 155 \pre neither segment is zero-length; i.e. p00 != p01 and p10 != p11.
 156 **/
 157 IntersectorKind
 158 line_segment_intersect(Geom::Point const &p00, Geom::Point const &p01,
 159                        Geom::Point const &p10, Geom::Point const &p11,
 160                        Geom::Point &result)
 161 {
 162     if(line_segment_intersectp(p00, p01, p10, p11)) {
 163         Geom::Point n0 = (p01 - p00).ccw();
 164         double d0 = dot(n0,p00);
 165
 166         Geom::Point n1 = (p11 - p10).ccw();
 167         double d1 = dot(n1,p10);
 168         return line_intersection(n0, d0, n1, d1, result);
 169     } else {
 170         return no_intersection;
 171     }
 172 }
 173
 174
 175 /** Determine whether \& where two line segments intersect.
 176
 177 If the two segments don't intersect, then \a result remains unchanged.
 178
 179 \pre neither segment is zero-length; i.e. p00 != p01 and p10 != p11.
 180 **/
 181 IntersectorKind
 182 segment_intersect(Geom::Point const &p00, Geom::Point const &p01,
 183                   Geom::Point const &p10, Geom::Point const &p11,
 184                   Geom::Point &result)
 185 {
 186     if(segment_intersectp(p00, p01, p10, p11)) {
 187         Geom::Point n0 = (p01 - p00).ccw();
 188         double d0 = dot(n0,p00);
 189
 190         Geom::Point n1 = (p11 - p10).ccw();
 191         double d1 = dot(n1,p10);
 192         return line_intersection(n0, d0, n1, d1, result);
 193     } else {
 194         return no_intersection;
 195     }
 196 }
 197
 198 /** Determine whether \& where two line segments intersect.
 199
 200 If the two segments don't intersect, then \a result remains unchanged.
 201
 202 \pre neither segment is zero-length; i.e. p00 != p01 and p10 != p11.
 203 **/
 204 IntersectorKind
 205 line_twopoint_intersect(Geom::Point const &p00, Geom::Point const &p01,
 206                         Geom::Point const &p10, Geom::Point const &p11,
 207                         Geom::Point &result)
 208 {
 209     Geom::Point n0 = (p01 - p00).ccw();
 210     double d0 = dot(n0,p00);
 211
 212     Geom::Point n1 = (p11 - p10).ccw();
 213     double d1 = dot(n1,p10);
 214     return line_intersection(n0, d0, n1, d1, result);
 215 }
 216
 217 // this is used to compare points for std::sort below
 218 static bool
 219 is_less(Point const &A, Point const &B)
 220 {
 221     if (A[X] < B[X]) {
 222         return true;
 223     } else if (A[X] == B[X] && A[Y] < B[Y]) {
 224         return true;
 225     } else {
 226         return false;
 227     }
 228 }
 229
 230 // TODO: this can doubtlessly be improved
 231 static void
 232 eliminate_duplicates_p(std::vector<Point> &pts)
 233 {
 234     unsigned int size = pts.size();
 235
 236     if (size < 2)
 237         return;
 238
 239     if (size == 2) {
 240         if (pts[0] == pts[1]) {
 241             pts.pop_back();
 242         }
 243     } else {
 244         std::sort(pts.begin(), pts.end(), &is_less);
 245         if (size == 3) {
 246             if (pts[0] == pts[1]) {
 247                 pts.erase(pts.begin());
 248             } else if (pts[1] == pts[2]) {
 249                 pts.pop_back();
 250             }
 251         } else {
 252             // we have size == 4
 253             if (pts[2] == pts[3]) {
 254                 pts.pop_back();
 255             }
 256             if (pts[0] == pts[1]) {
 257                 pts.erase(pts.begin());
 258             }
 259         }
 260     }
 261 }
 262
 263 /** Determine whether \& where an (infinite) line intersects a rectangle.
 264  *
 265  * \a c0, \a c1 are diagonal corners of the rectangle and
 266  * \a p1, \a p1 are distinct points on the line
 267  *
 268  * \return A list (possibly empty) of points of intersection.  If two such points (say \a r0 and \a
 269  * r1) then it is guaranteed that the order of \a r0, \a r1 along the line is the same as the that
 270  * of \a c0, \a c1 (i.e., the vectors \a r1 - \a r0 and \a p1 - \a p0 point into the same
 271  * direction).
 272  */
 273 std::vector<Geom::Point>
 274 rect_line_intersect(Geom::Point const &c0, Geom::Point const &c1,
 275                     Geom::Point const &p0, Geom::Point const &p1)
 276 {
 277     using namespace Geom;
 278
 279     std::vector<Point> results;
 280
 281     Point A(c0);
 282     Point C(c1);
 283
 284     Point B(A[X], C[Y]);
 285     Point D(C[X], A[Y]);
 286
 287     Point res;
 288
 289     if (line_segment_intersect(p0, p1, A, B, res) == intersects) {
 290         results.push_back(res);
 291     }
 292     if (line_segment_intersect(p0, p1, B, C, res) == intersects) {
 293         results.push_back(res);
 294     }
 295     if (line_segment_intersect(p0, p1, C, D, res) == intersects) {
 296         results.push_back(res);
 297     }
 298     if (line_segment_intersect(p0, p1, D, A, res) == intersects) {
 299         results.push_back(res);
 300     }
 301
 302     eliminate_duplicates_p(results);
 303
 304     if (results.size() == 2) {
 305         // sort the results so that the order is the same as that of p0 and p1
 306         Point dir1 (results[1] - results[0]);
 307         Point dir2 (p1 - p0);
 308         if (dot(dir1, dir2) < 0) {
 309             std::swap(results[0], results[1]);
 310         }
 311     }
 312
 313     return results;
 314 }
 315
 316 /**
 317  * polyCentroid: Calculates the centroid (xCentroid, yCentroid) and area of a polygon, given its
 318  * vertices (x[0], y[0]) ... (x[n-1], y[n-1]). It is assumed that the contour is closed, i.e., that
 319  * the vertex following (x[n-1], y[n-1]) is (x[0], y[0]).  The algebraic sign of the area is
 320  * positive for counterclockwise ordering of vertices in x-y plane; otherwise negative.
 321
 322  * Returned values:
 323     0 for normal execution;
 324     1 if the polygon is degenerate (number of vertices < 3);
 325     2 if area = 0 (and the centroid is undefined).
 326
 327     * for now we require the path to be a polyline and assume it is closed.
 328 **/
 329
 330 int centroid(std::vector<Geom::Point> p, Geom::Point& centroid, double &area) {
 331     const unsigned n = p.size();
 332     if (n < 3)
 333         return 1;
 334     Geom::Point centroid_tmp(0,0);
 335     double atmp = 0;
 336     for (unsigned i = n-1, j = 0; j < n; i = j, j++) {
 337         const double ai = -cross(p[j], p[i]);
 338         atmp += ai;
 339         centroid_tmp += (p[j] + p[i])*ai; // first moment.
 340     }
 341     area = atmp / 2;
 342     if (atmp != 0) {
 343         centroid = centroid_tmp / (3 * atmp);
 344         return 0;
 345     }
 346     return 2;
 347 }
 348
 349 }
 350
 351 /*
 352   Local Variables:
 353   mode:c++
 354   c-file-style:"stroustrup"
 355   c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
 356   indent-tabs-mode:nil
 357   fill-column:99
 358   End:
 359 */
 360 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :