src/2geom/geom.cpp

   1 /**
   2  *  \file src/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
  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 then 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 two line segments intersect.  This doesn't find
 114     the point of intersection, use the line_intersect function above,
 115     or the segment_intersection interface below.
 116
 117     \pre neither segment is zero-length; i.e. p00 != p01 and p10 != p11.
 118 */
 119 static bool
 120 segment_intersectp(Geom::Point const &p00, Geom::Point const &p01,
 121                                Geom::Point const &p10, Geom::Point const &p11)
 122 {
 123     if(p00 == p01) return false;
 124     if(p10 == p11) return false;
 125
 126     /* true iff (    (the p1 segment straddles the p0 infinite line)
 127      *           and (the p0 segment straddles the p1 infinite line) ). */
 128     return ((intersector_ccw(p00,p01, p10)
 129              *intersector_ccw(p00, p01, p11)) <=0 )
 130         &&
 131         ((intersector_ccw(p10,p11, p00)
 132           *intersector_ccw(p10, p11, p01)) <=0 );
 133 }
 134
 135
 136 /** Determine whether \& where two line segments intersect.
 137
 138 If the two segments don't intersect, then \a result remains unchanged.
 139
 140 \pre neither segment is zero-length; i.e. p00 != p01 and p10 != p11.
 141 **/
 142 IntersectorKind
 143 segment_intersect(Geom::Point const &p00, Geom::Point const &p01,
 144                               Geom::Point const &p10, Geom::Point const &p11,
 145                               Geom::Point &result)
 146 {
 147     if(segment_intersectp(p00, p01, p10, p11)) {
 148         Geom::Point n0 = (p01 - p00).ccw();
 149         double d0 = dot(n0,p00);
 150
 151         Geom::Point n1 = (p11 - p10).ccw();
 152         double d1 = dot(n1,p10);
 153         return line_intersection(n0, d0, n1, d1, result);
 154     } else {
 155         return no_intersection;
 156     }
 157 }
 158
 159 /** Determine whether \& where two line segments intersect.
 160
 161 If the two segments don't intersect, then \a result remains unchanged.
 162
 163 \pre neither segment is zero-length; i.e. p00 != p01 and p10 != p11.
 164 **/
 165 IntersectorKind
 166 line_twopoint_intersect(Geom::Point const &p00, Geom::Point const &p01,
 167                         Geom::Point const &p10, Geom::Point const &p11,
 168                         Geom::Point &result)
 169 {
 170     Geom::Point n0 = (p01 - p00).ccw();
 171     double d0 = dot(n0,p00);
 172
 173     Geom::Point n1 = (p11 - p10).ccw();
 174     double d1 = dot(n1,p10);
 175     return line_intersection(n0, d0, n1, d1, result);
 176 }
 177
 178 /**
 179  * polyCentroid: Calculates the centroid (xCentroid, yCentroid) and area of a polygon, given its
 180  * vertices (x[0], y[0]) ... (x[n-1], y[n-1]). It is assumed that the contour is closed, i.e., that
 181  * the vertex following (x[n-1], y[n-1]) is (x[0], y[0]).  The algebraic sign of the area is
 182  * positive for counterclockwise ordering of vertices in x-y plane; otherwise negative.
 183
 184  * Returned values:
 185     0 for normal execution;
 186     1 if the polygon is degenerate (number of vertices < 3);
 187     2 if area = 0 (and the centroid is undefined).
 188
 189     * for now we require the path to be a polyline and assume it is closed.
 190 **/
 191
 192 int centroid(std::vector<Geom::Point> p, Geom::Point& centroid, double &area) {
 193     const unsigned n = p.size();
 194     if (n < 3)
 195         return 1;
 196     Geom::Point centroid_tmp(0,0);
 197     double atmp = 0;
 198     for (unsigned i = n-1, j = 0; j < n; i = j, j++) {
 199         const double ai = -cross(p[j], p[i]);
 200         atmp += ai;
 201         centroid_tmp += (p[j] + p[i])*ai; // first moment.
 202     }
 203     area = atmp / 2;
 204     if (atmp != 0) {
 205         centroid = centroid_tmp / (3 * atmp);
 206         return 0;
 207     }
 208     return 2;
 209 }
 210
 211 }
 212
 213 /*
 214   Local Variables:
 215   mode:c++
 216   c-file-style:"stroustrup"
 217   c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
 218   indent-tabs-mode:nil
 219   fill-column:99
 220   End:
 221 */
 222 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :