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 "geom.h"
  10 #include "point.h"
  11
  12 namespace Geom {
  13 /**
  14  * Finds the intersection of the two (infinite) lines
  15  * defined by the points p such that dot(n0, p) == d0 and dot(n1, p) == d1.
  16  *
  17  * If the two lines intersect, then \a result becomes their point of
  18  * intersection; otherwise, \a result remains unchanged.
  19  *
  20  * This function finds the intersection of the two lines (infinite)
  21  * defined by n0.X = d0 and x1.X = d1.  The algorithm is as follows:
  22  * To compute the intersection point use kramer's rule:
  23  * \verbatim
  24  * convert lines to form
  25  * ax + by = c
  26  * dx + ey = f
  27  *
  28  * (
  29  *  e.g. a = (x2 - x1), b = (y2 - y1), c = (x2 - x1)*x1 + (y2 - y1)*y1
  30  * )
  31  *
  32  * In our case we use:
  33  *   a = n0.x     d = n1.x
  34  *   b = n0.y     e = n1.y
  35  *   c = d0        f = d1
  36  *
  37  * so:
  38  *
  39  * adx + bdy = cd
  40  * adx + aey = af
  41  *
  42  * bdy - aey = cd - af
  43  * (bd - ae)y = cd - af
  44  *
  45  * y = (cd - af)/(bd - ae)
  46  *
  47  * repeat for x and you get:
  48  *
  49  * x = (fb - ce)/(bd - ae)                \endverbatim
  50  *
  51  * If the denominator (bd-ae) is 0 then the lines are parallel, if the
  52  * numerators are then 0 then the lines coincide.
  53  *
  54  * \todo Why not use existing but outcommented code below
  55  * (HAVE_NEW_INTERSECTOR_CODE)?
  56  */
  57 IntersectorKind
  58 line_intersection(Geom::Point const &n0, double const d0,
  59                   Geom::Point const &n1, double const d1,
  60                   Geom::Point &result)
  61 {
  62     double denominator = dot(Geom::rot90(n0), n1);
  63     double X = n1[Geom::Y] * d0 -
  64         n0[Geom::Y] * d1;
  65     /* X = (-d1, d0) dot (n0[Y], n1[Y]) */
  66
  67     if (denominator == 0) {
  68         if ( X == 0 ) {
  69             return coincident;
  70         } else {
  71             return parallel;
  72         }
  73     }
  74
  75     double Y = n0[Geom::X] * d1 -
  76         n1[Geom::X] * d0;
  77
  78     result = Geom::Point(X, Y) / denominator;
  79
  80     return intersects;
  81 }
  82
  83
  84
  85
  86 /* ccw exists as a building block */
  87 int
  88 intersector_ccw(const Geom::Point& p0, const Geom::Point& p1,
  89         const Geom::Point& p2)
  90 /* Determine which way a set of three points winds. */
  91 {
  92     Geom::Point d1 = p1 - p0;
  93     Geom::Point d2 = p2 - p0;
  94     /* compare slopes but avoid division operation */
  95     double c = dot(Geom::rot90(d1), d2);
  96     if(c > 0)
  97         return +1; // ccw - do these match def'n in header?
  98     if(c < 0)
  99         return -1; // cw
 100
 101     /* Colinear [or NaN].  Decide the order. */
 102     if ( ( d1[0] * d2[0] < 0 )  ||
 103          ( d1[1] * d2[1] < 0 ) ) {
 104         return -1; // p2  <  p0 < p1
 105     } else if ( dot(d1,d1) < dot(d2,d2) ) {
 106         return +1; // p0 <= p1  <  p2
 107     } else {
 108         return 0; // p0 <= p2 <= p1
 109     }
 110 }
 111
 112 /** Determine whether two line segments intersect.  This doesn't find
 113     the point of intersection, use the line_intersect function above,
 114     or the segment_intersection interface below.
 115
 116     \pre neither segment is zero-length; i.e. p00 != p01 and p10 != p11.
 117 */
 118 static bool
 119 segment_intersectp(Geom::Point const &p00, Geom::Point const &p01,
 120                                Geom::Point const &p10, Geom::Point const &p11)
 121 {
 122     if(p00 == p01) return false;
 123     if(p10 == p11) return false;
 124
 125     /* true iff (    (the p1 segment straddles the p0 infinite line)
 126      *           and (the p0 segment straddles the p1 infinite line) ). */
 127     return ((intersector_ccw(p00,p01, p10)
 128              *intersector_ccw(p00, p01, p11)) <=0 )
 129         &&
 130         ((intersector_ccw(p10,p11, p00)
 131           *intersector_ccw(p10, p11, p01)) <=0 );
 132 }
 133
 134
 135 /** Determine whether \& where two line segments intersect.
 136
 137 If the two segments don't intersect, then \a result remains unchanged.
 138
 139 \pre neither segment is zero-length; i.e. p00 != p01 and p10 != p11.
 140 **/
 141 IntersectorKind
 142 segment_intersect(Geom::Point const &p00, Geom::Point const &p01,
 143                               Geom::Point const &p10, Geom::Point const &p11,
 144                               Geom::Point &result)
 145 {
 146     if(segment_intersectp(p00, p01, p10, p11)) {
 147         Geom::Point n0 = (p01 - p00).ccw();
 148         double d0 = dot(n0,p00);
 149
 150         Geom::Point n1 = (p11 - p10).ccw();
 151         double d1 = dot(n1,p10);
 152         return line_intersection(n0, d0, n1, d1, result);
 153     } else {
 154         return no_intersection;
 155     }
 156 }
 157
 158 /** Determine whether \& where two line segments intersect.
 159
 160 If the two segments don't intersect, then \a result remains unchanged.
 161
 162 \pre neither segment is zero-length; i.e. p00 != p01 and p10 != p11.
 163 **/
 164 IntersectorKind
 165 line_twopoint_intersect(Geom::Point const &p00, Geom::Point const &p01,
 166                         Geom::Point const &p10, Geom::Point const &p11,
 167                         Geom::Point &result)
 168 {
 169     Geom::Point n0 = (p01 - p00).ccw();
 170     double d0 = dot(n0,p00);
 171
 172     Geom::Point n1 = (p11 - p10).ccw();
 173     double d1 = dot(n1,p10);
 174     return line_intersection(n0, d0, n1, d1, result);
 175 }
 176
 177 /**
 178  * polyCentroid: Calculates the centroid (xCentroid, yCentroid) and area of a polygon, given its
 179  * vertices (x[0], y[0]) ... (x[n-1], y[n-1]). It is assumed that the contour is closed, i.e., that
 180  * the vertex following (x[n-1], y[n-1]) is (x[0], y[0]).  The algebraic sign of the area is
 181  * positive for counterclockwise ordering of vertices in x-y plane; otherwise negative.
 182
 183  * Returned values:
 184     0 for normal execution;
 185     1 if the polygon is degenerate (number of vertices < 3);
 186     2 if area = 0 (and the centroid is undefined).
 187
 188     * for now we require the path to be a polyline and assume it is closed.
 189 **/
 190
 191 int centroid(std::vector<Geom::Point> p, Geom::Point& centroid, double &area) {
 192     const unsigned n = p.size();
 193     if (n < 3)
 194         return 1;
 195     Geom::Point centroid_tmp(0,0);
 196     double atmp = 0;
 197     for (unsigned i = n-1, j = 0; j < n; i = j, j++) {
 198         const double ai = -cross(p[j], p[i]);
 199         atmp += ai;
 200         centroid_tmp += (p[j] + p[i])*ai; // first moment.
 201     }
 202     area = atmp / 2;
 203     if (atmp != 0) {
 204         centroid = centroid_tmp / (3 * atmp);
 205         return 0;
 206     }
 207     return 2;
 208 }
 209
 210 }
 211
 212 /*
 213   Local Variables:
 214   mode:c++
 215   c-file-style:"stroustrup"
 216   c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
 217   indent-tabs-mode:nil
 218   fill-column:99
 219   End:
 220 */
 221 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :