src/2geom/geom.cpp

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