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