allow polygon vertices to be specified incrementally for geodesic area (#221).

git-svn-id: http://svn.osgeo.org/metacrs/proj/trunk@2377 4e78687f-474d-0410-85f9-8d5e500ac6b2

2 files changed, 495 insertions, 75 deletions

diff --git a/src/geodesic.c b/src/geodesic.c
index 957f2144..bd9fc960 100644
--- a/src/geodesic.c
+++ b/src/geodesic.c
@@ -1,8 +1,12 @@
 /**
  * \file geodesic.c
  * \brief Implementation of the geodesic routines in C
+ *
+ * For the full documentation see geodesic.h.
  **********************************************************************/
 
+/** @cond SKIP */
+
 /*
  * This is a C implementation of the geodesic algorithms described in
  *
@@ -20,8 +24,6 @@
  */
 
 #include "geodesic.h"
-
-/** @cond SKIP */
 #include <math.h>
 
 #define GEOGRAPHICLIB_GEODESIC_ORDER 6
@@ -221,9 +223,10 @@ static real A2m1f(real eps);
 static void C2f(real eps, real c[]);
 static int transit(real lon1, real lon2);
 static void accini(real s[]);
+static void acccopy(const real s[], real t[]);
 static void accadd(real s[], real y);
-
-/** @endcond */
+static real accsum(const real s[], real y);
+static void accneg(real s[]);
 
 void geod_init(struct geod_geodesic* g, real a, real f) {
   if (!init) Init();
@@ -925,16 +928,12 @@ real geod_geninverse(const struct geod_geodesic* g,
 }
 
 void geod_inverse(const struct geod_geodesic* g,
-                  double lat1, double lon1, double lat2, double lon2,
-                  double* ps12, double* pazi1, double* pazi2) {
+                  real lat1, real lon1, real lat2, real lon2,
+                  real* ps12, real* pazi1, real* pazi2) {
   geod_geninverse(g, lat1, lon1, lat2, lon2, ps12, pazi1, pazi2, 0, 0, 0, 0);
 }
 
-/** @cond SKIP */
-
-real SinCosSeries(boolx sinp,
-                  real sinx, real cosx,
-                  const real c[], int n) {
+real SinCosSeries(boolx sinp, real sinx, real cosx, const real c[], int n) {
   /* Evaluate
    * y = sinp ? sum(c[i] * sin( 2*i    * x), i, 1, n) :
    *            sum(c[i] * cos((2*i+1) * x), i, 0, n-1)
@@ -1522,6 +1521,11 @@ void accini(real s[]) {
   s[0] = s[1] = 0;
 }
 
+void acccopy(const real s[], real t[]) {
+  /* Copy an accumulator; t = s. */
+  t[0] = s[0]; t[1] = s[1];
+}
+
 void accadd(real s[], real y) {
   /* Add y to an accumulator. */
   real u, z = sumx(y, s[1], &u);
@@ -1532,30 +1536,235 @@ void accadd(real s[], real y) {
     s[1] = s[1] + u;
 }
 
-/** @endcond */
+real accsum(const real s[], real y) {
+  /* Return accumulator + y (but don't add to accumulator). */
+  real t[2];
+  acccopy(s, t);
+  accadd(t, y);
+  return t[0];
+}
 
-void geod_polygonarea(const struct geod_geodesic* g,
-                      real lats[], real lons[], int n,
-                      real* pA, real* pP) {
-  int i, crossings = 0;
-  real area0 = 4 * pi * g->c2, A[2], P[2];
-  accini(A); accini(P);
-  for (i = 0; i < n; ++i) {
+void accneg(real s[]) {
+  /* Negate an accumulator. */
+  s[0] = -s[0]; s[1] = -s[1];
+}
+
+void geod_polygon_init(struct geod_polygon* p, boolx polylinep) {
+  p->lat0 = p->lon0 = p->lat = p->lon = NaN;
+  p->polyline = (polylinep != 0);
+  accini(p->P);
+  accini(p->A);
+  p->num = p->crossings = 0;
+}
+
+void geod_polygon_addpoint(const struct geod_geodesic* g,
+                           struct geod_polygon* p,
+                           real lat, real lon) {
+  lon = AngNormalize(lon);
+  if (p->num == 0) {
+    p->lat0 = p->lat = lat;
+    p->lon0 = p->lon = lon;
+  } else {
+    real s12, S12;
+    geod_geninverse(g, p->lat, p->lon, lat, lon,
+                    &s12, 0, 0, 0, 0, 0, p->polyline ? 0 : &S12);
+    accadd(p->P, s12);
+    if (!p->polyline) {
+      accadd(p->A, S12);
+      p->crossings += transit(p->lon, lon);
+    }
+    p->lat = lat; p->lon = lon;
+  }
+  ++p->num;
+}
+
+void geod_polygon_addedge(const struct geod_geodesic* g,
+                          struct geod_polygon* p,
+                          real azi, real s) {
+  if (p->num) {                 /* Do nothing is num is zero */
+    real lat, lon, S12;
+    geod_gendirect(g, p->lat, p->lon, azi, FALSE, s,
+                   &lat, &lon, 0,
+                   0, 0, 0, 0, p->polyline ? 0 : &S12);
+    accadd(p->P, s);
+    if (!p->polyline) {
+      accadd(p->A, S12);
+      p->crossings += transit(p->lon, lon);
+    }
+    p->lat = lat; p->lon = lon;
+    ++p->num;
+  }
+}
+
+unsigned geod_polygon_compute(const struct geod_geodesic* g,
+                              const struct geod_polygon* p,
+                              boolx reverse, boolx sign,
+                              real* pA, real* pP) {
+  real s12, S12, t[2], area0;
+  int crossings;
+  if (p->num < 2) {
+    if (pP) *pP = 0;
+    if (!p->polyline && pA) *pA = 0;
+    return p->num;
+  }
+  if (p->polyline) {
+    if (pP) *pP = p->P[0];
+    return p->num;
+  }
+  geod_geninverse(g, p->lat, p->lon, p->lat0, p->lon0,
+                  &s12, 0, 0, 0, 0, 0, &S12);
+  if (pP) *pP = accsum(p->P, s12);
+  acccopy(p->A, t);
+  accadd(t, S12);
+  crossings = p->crossings + transit(p->lon, p->lon0);
+  area0 = 4 * pi * g->c2;
+  if (crossings & 1)
+    accadd(t, (t[0] < 0 ? 1 : -1) * area0/2);
+  /* area is with the clockwise sense.  If !reverse convert to
+   * counter-clockwise convention. */
+  if (!reverse)
+    accneg(t);
+  /* If sign put area in (-area0/2, area0/2], else put area in [0, area0) */
+  if (sign) {
+    if (t[0] > area0/2)
+      accadd(t, -area0);
+    else if (t[0] <= -area0/2)
+      accadd(t, +area0);
+  } else {
+    if (t[0] >= area0)
+      accadd(t, -area0);
+    else if (t[0] < 0)
+      accadd(t, +area0);
+  }
+  if (pA) *pA = 0 + t[0];
+  return p->num;
+}
+
+unsigned geod_polygon_testpoint(const struct geod_geodesic* g,
+                                const struct geod_polygon* p,
+                                real lat, real lon,
+                                boolx reverse, boolx sign,
+                                real* pA, real* pP) {
+  real perimeter, tempsum, area0;
+  int crossings, i;
+  unsigned num = p->num + 1;
+  if (num == 1) {
+    if (pP) *pP = 0;
+    if (!p->polyline && pA) *pA = 0;
+    return num;
+  }
+  perimeter = p->P[0];
+  tempsum = p->polyline ? 0 : p->A[0];
+  crossings = p->crossings;
+  for (i = 0; i < (p->polyline ? 1 : 2); ++i) {
     real s12, S12;
-    geod_geninverse(g, lats[i], lons[i], lats[(i + 1) % n], lons[(i + 1) % n],
-               &s12, 0, 0, 0, 0, 0, &S12);
-    accadd(P, s12);
-    /* The minus sign is due to the counter-clockwise convention */
-    accadd(A, -S12);
-    crossings += transit(lons[i], lons[(i + 1) % n]);
+    geod_geninverse(g,
+                    i == 0 ? p->lat  : lat, i == 0 ? p->lon  : lon,
+                    i != 0 ? p->lat0 : lat, i != 0 ? p->lon0 : lon,
+                    &s12, 0, 0, 0, 0, 0, p->polyline ? 0 : &S12);
+    perimeter += s12;
+    if (!p->polyline) {
+      tempsum += S12;
+      crossings += transit(i == 0 ? p->lon  : lon,
+                           i != 0 ? p->lon0 : lon);
+    }
   }
+
+  if (pP) *pP = perimeter;
+  if (p->polyline)
+    return num;
+
+  area0 = 4 * pi * g->c2;
   if (crossings & 1)
-    accadd(A, (A[0] < 0 ? 1 : -1) * area0/2);
-  /* Put area in (-area0/2, area0/2] */
-  if (A[0] > area0/2)
-    accadd(A, -area0);
-  else if (A[0] <= -area0/2)
-    accadd(A, +area0);
-  if (pA) *pA = A[0];
-  if (pP) *pP = P[0];
+    tempsum += (tempsum < 0 ? 1 : -1) * area0/2;
+  /* area is with the clockwise sense.  If !reverse convert to
+   * counter-clockwise convention. */
+  if (!reverse)
+    tempsum *= -1;
+  /* If sign put area in (-area0/2, area0/2], else put area in [0, area0) */
+  if (sign) {
+    if (tempsum > area0/2)
+      tempsum -= area0;
+    else if (tempsum <= -area0/2)
+      tempsum += area0;
+  } else {
+    if (tempsum >= area0)
+      tempsum -= area0;
+    else if (tempsum < 0)
+      tempsum += area0;
+  }
+  if (pA) *pA = 0 + tempsum;
+  return num;
 }
+
+unsigned geod_polygon_testedge(const struct geod_geodesic* g,
+                               const struct geod_polygon* p,
+                               real azi, real s,
+                               boolx reverse, boolx sign,
+                               real* pA, real* pP) {
+  real perimeter, tempsum, area0;
+  int crossings;
+  unsigned num = p->num + 1;
+  if (num == 1) {               /* we don't have a starting point! */
+    if (pP) *pP = NaN;
+    if (!p->polyline && pA) *pA = NaN;
+    return 0;
+  }
+  perimeter = p->P[0] + s;
+  if (p->polyline) {
+    if (pP) *pP = perimeter;
+    return num;
+  }
+
+  tempsum = p->A[0];
+  crossings = p->crossings;
+  {
+    real lat, lon, s12, S12;
+    geod_gendirect(g, p->lat, p->lon, azi, FALSE, s,
+                   &lat, &lon, 0,
+                   0, 0, 0, 0, &S12);
+    tempsum += S12;
+    crossings += transit(p->lon, lon);
+    geod_geninverse(g, lat,  lon, p->lat0,  p->lon0,
+                    &s12, 0, 0, 0, 0, 0, &S12);
+    perimeter += s12;
+    tempsum += S12;
+    crossings += transit(lon, p->lon0);
+  }
+
+  area0 = 4 * pi * g->c2;
+  if (crossings & 1)
+    tempsum += (tempsum < 0 ? 1 : -1) * area0/2;
+  /* area is with the clockwise sense.  If !reverse convert to
+   * counter-clockwise convention. */
+  if (!reverse)
+    tempsum *= -1;
+  /* If sign put area in (-area0/2, area0/2], else put area in [0, area0) */
+  if (sign) {
+    if (tempsum > area0/2)
+      tempsum -= area0;
+    else if (tempsum <= -area0/2)
+      tempsum += area0;
+  } else {
+    if (tempsum >= area0)
+      tempsum -= area0;
+    else if (tempsum < 0)
+      tempsum += area0;
+  }
+  if (pP) *pP = perimeter;
+  if (pA) *pA = 0 + tempsum;
+  return num;
+}
+
+void geod_polygonarea(const struct geod_geodesic* g,
+                      real lats[], real lons[], int n,
+                      real* pA, real* pP) {
+  int i;
+  struct geod_polygon p;
+  geod_polygon_init(&p, FALSE);
+  for (i = 0; i < n; ++i)
+    geod_polygon_addpoint(g, &p, lats[i], lons[i]);
+  geod_polygon_compute(g, &p, FALSE, TRUE, pA, pP);
+}
+
+/** @endcond */
diff --git a/src/geodesic.h b/src/geodesic.h
index 87d087a5..69985009 100644
--- a/src/geodesic.h
+++ b/src/geodesic.h
@@ -27,8 +27,8 @@
  * - the direct problem -- given \e lat1, \e lon1, \e s12, and \e azi1,
  *   determine \e lat2, \e lon2, and \e azi2.  This is solved by the function
  *   geod_direct().
- * - the inverse problem -- given \e lat1, \e lon1, \e lat2, \e lon2, determine
- *   \e s12, \e azi1, and \e azi2.  This is solved by the function
+ * - the inverse problem -- given \e lat1, \e lon1, and \e lat2, \e lon2,
+ *   determine \e s12, \e azi1, and \e azi2.  This is solved by the function
  *   geod_inverse().
  *
  * The ellipsoid is specified by its equatorial radius \e a (typically in
@@ -97,31 +97,44 @@
  *   azi2] + [\e d, \e d], for arbitrary \e d.
  *
  * These routines are a simple transcription of the corresponding C++ classes
- * in <a href="http://geographiclib.sf.net"> GeographicLib</a>.  The
- * "class data" is represented by the structs geod_geodesic and
- * geod_geodesicline and these are passed as initial arguments to the member
- * functions.  (But note that the PolygonArea class has been replaced by a
- * simple function interface and the area computation does use the accurate
- * summing providing by the Accumulator class.)  Most of the internal comments
- * have been retained.  However, in the process of transcription some
- * documentation has been lost and the documentation for the C++ classes,
- * GeographicLib::Geodesic, GeographicLib::GeodesicLine, and
- * GeographicLib::PolygonArea, should be consulted.  The C++ code remains the
- * "reference implementation".  Think twice about restructuring the internals
- * of the C code since this may make porting fixes from the C++ code more
- * difficult.
+ * in <a href="http://geographiclib.sf.net"> GeographicLib</a>.  The "class
+ * data" is represented by the structs geod_geodesic, geod_geodesicline,
+ * geod_polygon and pointers to these objects are passed as initial arguments
+ * to the member functions.  Most of the internal comments have been retained.
+ * However, in the process of transcription some documentation has been lost
+ * and the documentation for the C++ classes, GeographicLib::Geodesic,
+ * GeographicLib::GeodesicLine, and GeographicLib::PolygonArea, should be
+ * consulted.  The C++ code remains the "reference implementation".  Think
+ * twice about restructuring the internals of the C code since this may make
+ * porting fixes from the C++ code more difficult.
  *
  * Copyright (c) Charles Karney (2012-2013) <charles@karney.com> and licensed
  * under the MIT/X11 License.  For more information, see
  * http://geographiclib.sourceforge.net/
  *
  * This library was distributed with
- * <a href="../index.html">GeographicLib</a> 1.31.
+ * <a href="../index.html">GeographicLib</a> 1.32.
  **********************************************************************/
 
 #if !defined(GEODESIC_H)
 #define GEODESIC_H 1
 
+/**
+ * The major version of the geodesic library.  (This tracks the version of
+ * GeographicLib.)
+ **********************************************************************/
+#define GEODESIC_VERSION_MAJOR 1
+/**
+ * The minor version of the geodesic library.  (This tracks the version of
+ * GeographicLib.)
+ **********************************************************************/
+#define GEODESIC_VERSION_MINOR 32
+/**
+ * The patch level of the geodesic library.  (This tracks the version of
+ * GeographicLib.)
+ **********************************************************************/
+#define GEODESIC_VERSION_PATCH 0
+
 #if defined(__cplusplus)
 extern "C" {
 #endif
@@ -159,30 +172,50 @@ extern "C" {
   };
 
   /**
-   * Initialize a geod_geodesic object
+   * The struct for accumulating information about a geodesic polygon.  This is
+   * used for computing the perimeter and area of a polygon.  This must be
+   * initialized by geod_polygon_init() before use.
+   **********************************************************************/
+  struct geod_polygon {
+    double lat;                 /**< the current latitude */
+    double lon;                 /**< the current longitude */
+    /**< @cond SKIP */
+    double lat0;
+    double lon0;
+    double A[2];
+    double P[2];
+    int polyline;
+    int crossings;
+    /**< @endcond */
+    unsigned num;               /**< the number of points so far */
+  };
+
+  /**
+   * Initialize a geod_geodesic object.
    *
-   * @param[out] g the object to be initialized.
+   * @param[out] g a pointer to the object to be initialized.
    * @param[in] a the equatorial radius (meters).
    * @param[in] f the flattening.
    **********************************************************************/
   void geod_init(struct geod_geodesic* g, double a, double f);
 
   /**
-   * Initialize a geod_geodesicline object
+   * Initialize a geod_geodesicline object.
    *
-   * @param[out] l the object to be initialized.
-   * @param[in] g the geod_geodesic object specifying the ellipsoid.
+   * @param[out] l a pointer to the object to be initialized.
+   * @param[in] g a pointer to the geod_geodesic object specifying the
+   *   ellipsoid.
    * @param[in] lat1 latitude of point 1 (degrees).
    * @param[in] lon1 longitude of point 1 (degrees).
    * @param[in] azi1 azimuth at point 1 (degrees).
-   * @param[in] caps bitor'ed combination of geod_mask values specifying the
+   * @param[in] caps bitor'ed combination of geod_mask() values specifying the
    *   capabilities the geod_geodesicline object should possess, i.e., which
    *   quantities can be returned in calls to geod_position() and
    *   geod_genposition().
    *
    * \e g must have been initialized with a call to geod_init().  \e lat1
-   * should be in the range [&minus;90&deg;, 90&deg;]; \e azi1 should be in the
-   * range [&minus;540&deg;, 540&deg;).
+   * should be in the range [&minus;90&deg;, 90&deg;]; \e lon1 and \e azi1
+   * should be in the range [&minus;540&deg;, 540&deg;).
    *
    * The geod_mask values are [see geod_mask()]:
    * - \e caps |= GEOD_LATITUDE for the latitude \e lat2; this is
@@ -210,7 +243,8 @@ extern "C" {
   /**
    * Solve the direct geodesic problem.
    *
-   * @param[in] g the geod_geodesic object specifying the ellipsoid.
+   * @param[in] g a pointer to the geod_geodesic object specifying the
+   *   ellipsoid.
    * @param[in] lat1 latitude of point 1 (degrees).
    * @param[in] lon1 longitude of point 1 (degrees).
    * @param[in] azi1 azimuth at point 1 (degrees).
@@ -224,7 +258,7 @@ extern "C" {
    * should be in the range [&minus;90&deg;, 90&deg;]; \e lon1 and \e azi1
    * should be in the range [&minus;540&deg;, 540&deg;).  The values of \e lon2
    * and \e azi2 returned are in the range [&minus;180&deg;, 180&deg;).  Any of
-   * the "return" arguments plat2, etc., may be replaced by 0, if you do not
+   * the "return" arguments \e plat2, etc., may be replaced by 0, if you do not
    * need some quantities computed.
    *
    * If either point is at a pole, the azimuth is defined by keeping the
@@ -250,7 +284,8 @@ extern "C" {
   /**
    * Solve the inverse geodesic problem.
    *
-   * @param[in] g the geod_geodesic object specifying the ellipsoid.
+   * @param[in] g a pointer to the geod_geodesic object specifying the
+   *   ellipsoid.
    * @param[in] lat1 latitude of point 1 (degrees).
    * @param[in] lon1 longitude of point 1 (degrees).
    * @param[in] lat2 latitude of point 2 (degrees).
@@ -264,8 +299,8 @@ extern "C" {
    * and \e lat2 should be in the range [&minus;90&deg;, 90&deg;]; \e lon1 and
    * \e lon2 should be in the range [&minus;540&deg;, 540&deg;).  The values of
    * \e azi1 and \e azi2 returned are in the range [&minus;180&deg;, 180&deg;).
-   * Any of the "return" arguments ps12, etc., may be replaced by 0, if you do
-   * not need some quantities computed.
+   * Any of the "return" arguments \e ps12, etc., may be replaced by 0, if you
+   * do not need some quantities computed.
    *
    * If either point is at a pole, the azimuth is defined by keeping the
    * longitude fixed, writing \e lat = &plusmn;(90&deg; &minus; &epsilon;), and
@@ -292,18 +327,19 @@ extern "C" {
   /**
    * Compute the position along a geod_geodesicline.
    *
-   * @param[in] l the geod_geodesicline object specifying the geodesic line.
+   * @param[in] l a pointer to the geod_geodesicline object specifying the
+   *   geodesic line.
    * @param[in] s12 distance between point 1 and point 2 (meters); it can be
    *   negative.
    * @param[out] plat2 pointer to the latitude of point 2 (degrees).
    * @param[out] plon2 pointer to the longitude of point 2 (degrees); requires
-   *   that the \e was initialized with \e caps |= GEOD_LONGITUDE.
+   *   that \e l was initialized with \e caps |= GEOD_LONGITUDE.
    * @param[out] pazi2 pointer to the (forward) azimuth at point 2 (degrees).
    *
    * \e l must have been initialized with a call to geod_lineinit() with \e
    * caps |= GEOD_DISTANCE_IN.  The values of \e lon2 and \e azi2 returned are
    * in the range [&minus;180&deg;, 180&deg;).  Any of the "return" arguments
-   * plat2, etc., may be replaced by 0, if you do not need some quantities
+   * \e plat2, etc., may be replaced by 0, if you do not need some quantities
    * computed.
    *
    * Example, compute way points between JFK and Singapore Changi Airport
@@ -340,7 +376,8 @@ extern "C" {
   /**
    * The general direct geodesic problem.
    *
-   * @param[in] g the geod_geodesic object specifying the ellipsoid.
+   * @param[in] g a pointer to the geod_geodesic object specifying the
+   *   ellipsoid.
    * @param[in] lat1 latitude of point 1 (degrees).
    * @param[in] lon1 longitude of point 1 (degrees).
    * @param[in] azi1 azimuth at point 1 (degrees).
@@ -367,7 +404,7 @@ extern "C" {
    * should be in the range [&minus;90&deg;, 90&deg;]; \e lon1 and \e azi1
    * should be in the range [&minus;540&deg;, 540&deg;).  The function value \e
    * a12 equals \e s12_a12 is \e arcmode is non-zero.  Any of the "return"
-   * arguments plat2, etc., may be replaced by 0, if you do not need some
+   * arguments \e plat2, etc., may be replaced by 0, if you do not need some
    * quantities computed.
    **********************************************************************/
   double geod_gendirect(const struct geod_geodesic* g,
@@ -380,7 +417,8 @@ extern "C" {
   /**
    * The general inverse geodesic calculation.
    *
-   * @param[in] g the geod_geodesic object specifying the ellipsoid.
+   * @param[in] g a pointer to the geod_geodesic object specifying the
+   *   ellipsoid.
    * @param[in] lat1 latitude of point 1 (degrees).
    * @param[in] lon1 longitude of point 1 (degrees).
    * @param[in] lat2 latitude of point 2 (degrees).
@@ -401,7 +439,7 @@ extern "C" {
    * \e g must have been initialized with a call to geod_init().  \e lat1
    * and \e lat2 should be in the range [&minus;90&deg;, 90&deg;]; \e lon1 and
    * \e lon2 should be in the range [&minus;540&deg;, 540&deg;).  Any of the
-   * "return" arguments ps12, etc., may be replaced by 0, if you do not need
+   * "return" arguments \e ps12, etc., may be replaced by 0, if you do not need
    * some quantities computed.
    **********************************************************************/
   double geod_geninverse(const struct geod_geodesic* g,
@@ -413,7 +451,8 @@ extern "C" {
   /**
    * The general position function.
    *
-   * @param[in] l the geod_geodesicline object specifying the geodesic line.
+   * @param[in] l a pointer to the geod_geodesicline object specifying the
+   *   geodesic line.
    * @param[in] arcmode flag determining the meaning of the second parameter;
    *   if arcmode is 0, then \e l must have been initialized with \e caps |=
    *   GEOD_DISTANCE_IN.
@@ -443,11 +482,11 @@ extern "C" {
    * \e l must have been initialized with a call to geod_lineinit() with \e
    * caps |= GEOD_DISTANCE_IN.  The values of \e lon2 and \e azi2 returned are
    * in the range [&minus;180&deg;, 180&deg;).  Any of the "return" arguments
-   * plat2, etc., may be replaced by 0, if you do not need some quantities
+   * \e plat2, etc., may be replaced by 0, if you do not need some quantities
    * computed.  Requesting a value which \e l is not capable of computing is
    * not an error; the corresponding argument will not be altered.
    *
-   * Example, computer way points between JFK and Singapore Changi Airport
+   * Example, compute way points between JFK and Singapore Changi Airport
    * using geod_genposition().  In this example, the points are evenly space in
    * arc length (and so only approximately equally space in distance).  This is
    * faster than using geod_position() would be appropriate if drawing the path
@@ -476,9 +515,181 @@ extern "C" {
                           double* pS12);
 
   /**
-   * The area of a geodesic polygon.
+   * Initialize a geod_polygon object.
+   *
+   * @param[out] p a pointer to the object to be initialized.
+   * @param[in] polylinep non-zero if a polyline instead of a polygon.
+   *
+   * If \e polylinep is zero, then the sequence of vertices and edges added by
+   * geod_polygon_addpoint() and geod_polygon_addedge() define a polygon and
+   * the perimeter and area are returned by geod_polygon_compute().  If \e
+   * polylinep is non-zero, then the vertices and edges define a polyline and
+   * only the perimeter is returned by geod_polygon_compute().
+   *
+   * An example of the use of this function is given in the documentation for
+   * geod_polygon_compute().
+   **********************************************************************/
+  void geod_polygon_init(struct geod_polygon* p, int polylinep);
+
+  /**
+   * Add a point to the polygon or polyline.
+   *
+   * @param[in] g a pointer to the geod_geodesic object specifying the
+   *   ellipsoid.
+   * @param[in,out] p a pointer to the geod_polygon object specifying the
+   *   polygon.
+   * @param[in] lat the latitude of the point (degrees).
+   * @param[in] lon the longitude of the point (degrees).
+   *
+   * \e g and \e p must have been initialized with calls to geod_init() and
+   * geod_polygon_init(), respectively.  The same \e g must be used for all the
+   * points and edges in a polygon.  \e lat should be in the range
+   * [&minus;90&deg;, 90&deg;] and \e lon should be in the range
+   * [&minus;540&deg;, 540&deg;).
+   *
+   * An example of the use of this function is given in the documentation for
+   * geod_polygon_compute().
+   **********************************************************************/
+  void geod_polygon_addpoint(const struct geod_geodesic* g,
+                             struct geod_polygon* p,
+                             double lat, double lon);
+
+  /**
+   * Add an edge to the polygon or polyline.
+   *
+   * @param[in] g a pointer to the geod_geodesic object specifying the
+   *   ellipsoid.
+   * @param[in,out] p a pointer to the geod_polygon object specifying the
+   *   polygon.
+   * @param[in] azi azimuth at current point (degrees).
+   * @param[in] s distance from current point to next point (meters).
+   *
+   * \e g and \e p must have been initialized with calls to geod_init() and
+   * geod_polygon_init(), respectively.  The same \e g must be used for all the
+   * points and edges in a polygon.  \e azi should be in the range
+   * [&minus;540&deg;, 540&deg;).  This does nothing if no points have been
+   * added yet.  The \e lat and \e lon fields of \e p give the location of
+   * the new vertex.
+   **********************************************************************/
+  void geod_polygon_addedge(const struct geod_geodesic* g,
+                            struct geod_polygon* p,
+                            double azi, double s);
+
+  /**
+   * Return the results for a polygon.
+   *
+   * @param[in] g a pointer to the geod_geodesic object specifying the
+   *   ellipsoid.
+   * @param[in] p a pointer to the geod_polygon object specifying the polygon.
+   * @param[in] reverse if non-zero then clockwise (instead of
+   *   counter-clockwise) traversal counts as a positive area.
+   * @param[in] sign if non-zero then return a signed result for the area if
+   *   the polygon is traversed in the "wrong" direction instead of returning
+   *   the area for the rest of the earth.
+   * @param[out] pA pointer to the area of the polygon (meters<sup>2</sup>);
+   *   only set if \e polyline is non-zero in the call to geod_polygon_init().
+   * @param[out] pP pointer to the perimeter of the polygon or length of the
+   *   polyline (meters).
+   * @return the number of points.
+   *
+   * Only simple polygons (which are not self-intersecting) are allowed.
+   * There's no need to "close" the polygon by repeating the first vertex.  Set
+   * \e pA or \e pP to zero, if you do not want the corresponding quantity
+   * returned.
+   *
+   * Example, compute the perimeter and area of the geodesic triangle with
+   * vertices (0&deg;N,0&deg;E), (0&deg;N,90&deg;E), (90&deg;N,0&deg;E).
+   @code
+   double A, P;
+   int n;
+   struct geod_geodesic g;
+   struct geod_polygon p;
+   geod_init(&g, 6378137, 1/298.257223563);
+   geod_polygon_init(&p, 0);
+
+   geod_polygon_addpoint(&g, &p,  0,  0);
+   geod_polygon_addpoint(&g, &p,  0, 90);
+   geod_polygon_addpoint(&g, &p, 90,  0);
+   n = geod_polygon_compute(&g, &p, 0, 1, &A, &P);
+   printf("%d %.8f %.3f\n", n, P, A);
+   @endcode
+   **********************************************************************/
+  unsigned geod_polygon_compute(const struct geod_geodesic* g,
+                                const struct geod_polygon* p,
+                                int reverse, int sign,
+                                double* pA, double* pP);
+
+  /**
+   * Return the results assuming a tentative final test point is added;
+   * however, the data for the test point is not saved.  This lets you report a
+   * running result for the perimeter and area as the user moves the mouse
+   * cursor.  Ordinary floating point arithmetic is used to accumulate the data
+   * for the test point; thus the area and perimeter returned are less accurate
+   * than if geod_polygon_addpoint() and geod_polygon_compute() are used.
+   *
+   * @param[in] g a pointer to the geod_geodesic object specifying the
+   *   ellipsoid.
+   * @param[in] p a pointer to the geod_polygon object specifying the polygon.
+   * @param[in] lat the latitude of the test point (degrees).
+   * @param[in] lon the longitude of the test point (degrees).
+   * @param[in] reverse if non-zero then clockwise (instead of
+   *   counter-clockwise) traversal counts as a positive area.
+   * @param[in] sign if non-zero then return a signed result for the area if
+   *   the polygon is traversed in the "wrong" direction instead of returning
+   *   the area for the rest of the earth.
+   * @param[out] pA pointer to the area of the polygon (meters<sup>2</sup>);
+   *   only set if \e polyline is non-zero in the call to geod_polygon_init().
+   * @param[out] pP pointer to the perimeter of the polygon or length of the
+   *   polyline (meters).
+   * @return the number of points.
+   *
+   * \e lat should be in the range [&minus;90&deg;, 90&deg;] and \e
+   * lon should be in the range [&minus;540&deg;, 540&deg;).
+   **********************************************************************/
+  unsigned geod_polygon_testpoint(const struct geod_geodesic* g,
+                                  const struct geod_polygon* p,
+                                  double lat, double lon,
+                                  int reverse, int sign,
+                                  double* pA, double* pP);
+
+  /**
+   * Return the results assuming a tentative final test point is added via an
+   * azimuth and distance; however, the data for the test point is not saved.
+   * This lets you report a running result for the perimeter and area as the
+   * user moves the mouse cursor.  Ordinary floating point arithmetic is used
+   * to accumulate the data for the test point; thus the area and perimeter
+   * returned are less accurate than if geod_polygon_addedge() and
+   * geod_polygon_compute() are used.
+   *
+   * @param[in] g a pointer to the geod_geodesic object specifying the
+   *   ellipsoid.
+   * @param[in] p a pointer to the geod_polygon object specifying the polygon.
+   * @param[in] azi azimuth at current point (degrees).
+   * @param[in] s distance from current point to final test point (meters).
+   * @param[in] reverse if non-zero then clockwise (instead of
+   *   counter-clockwise) traversal counts as a positive area.
+   * @param[in] sign if non-zero then return a signed result for the area if
+   *   the polygon is traversed in the "wrong" direction instead of returning
+   *   the area for the rest of the earth.
+   * @param[out] pA pointer to the area of the polygon (meters<sup>2</sup>);
+   *   only set if \e polyline is non-zero in the call to geod_polygon_init().
+   * @param[out] pP pointer to the perimeter of the polygon or length of the
+   *   polyline (meters).
+   * @return the number of points.
+   *
+   * \e azi should be in the range [&minus;540&deg;, 540&deg;).
+   **********************************************************************/
+  unsigned geod_polygon_testedge(const struct geod_geodesic* g,
+                                 const struct geod_polygon* p,
+                                 double azi, double s,
+                                 int reverse, int sign,
+                                 double* pA, double* pP);
+
+  /**
+   * A simple interface for computing the area of a geodesic polygon.
    *
-   * @param[in] g the geod_geodesic object specifying the ellipsoid.
+   * @param[in] g a pointer to the geod_geodesic object specifying the
+   *   ellipsoid.
    * @param[in] lats an array of latitudes of the polygon vertices (degrees).
    * @param[in] lons an array of longitudes of the polygon vertices (degrees).
    * @param[in] n the number of vertices.
@@ -512,7 +723,7 @@ extern "C" {
                         double* pA, double* pP);
 
   /**
-   * mask values for geod_geodesicline capabilities.
+   * mask values for the the \e caps argument to geod_lineinit().
    **********************************************************************/
   enum geod_mask {
     GEOD_NONE         = 0U,                     /**< Calculate nothing */