/usr/include/GeographicLib/PolygonArea.hpp is in libgeographiclib-dev 1.21-1ubuntu1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 | /**
* \file PolygonArea.hpp
* \brief Header for GeographicLib::PolygonArea class
*
* Copyright (c) Charles Karney (2010, 2011) <charles@karney.com> and licensed
* under the MIT/X11 License. For more information, see
* http://geographiclib.sourceforge.net/
**********************************************************************/
#if !defined(GEOGRAPHICLIB_POLYGONAREA_HPP)
#define GEOGRAPHICLIB_POLYGONAREA_HPP \
"$Id: 7a339f312a9c977b9fccad3c0c8bfa9009d863e2 $"
#include <GeographicLib/Geodesic.hpp>
#include <GeographicLib/Constants.hpp>
#include <GeographicLib/Accumulator.hpp>
namespace GeographicLib {
/**
* \brief Polygon Areas.
*
* This computes the area of a geodesic polygon using the method given
* Section 15 of
* - C. F. F. Karney,
* <a href="http://arxiv.org/abs/1102.1215v1">Geodesics
* on an ellipsoid of revolution</a>,
* Feb. 2011;
* preprint
* <a href="http://arxiv.org/abs/1102.1215v1">arxiv:1102.1215v1</a>.
* .
* See also Section 6 of
* - C. F. F. Karney,
* <a href="http://arxiv.org/abs/1109.4448">Algorithms for geodesics</a>,
* Sept. 2011;
* preprint
* <a href="http://arxiv.org/abs/1109.4448">arxiv:1109.4448</a>.
*
* This class lets you add vertices one at a time to the polygon. The area
* and perimeter are accumulated in two times the standard floating point
* precision to guard against the loss of accuracy with many-sided polygons.
* At any point you can ask for the perimeter and area so far. There's an
* option to treat the points as defining a polyline instead of a polygon; in
* that case, only the perimeter is computed.
*
* Example of use:
* \include example-PolygonArea.cpp
*
* <a href="Planimeter.1.html">Planimeter</a> is a command-line utility
* providing access to the functionality of PolygonArea.
**********************************************************************/
class GEOGRAPHIC_EXPORT PolygonArea {
private:
typedef Math::real real;
Geodesic _earth;
real _area0; // Full ellipsoid area
bool _polyline; // Assume polyline (don't close and skip area)
unsigned _mask;
unsigned _num;
int _crossings;
Accumulator<real> _areasum, _perimetersum;
real _lat0, _lon0, _lat1, _lon1;
// Copied from Geodesic class
static inline real AngNormalize(real x) throw() {
// Place angle in [-180, 180). Assumes x is in [-540, 540).
//
// g++ 4.4.4 holds a temporary in an extended register causing an error
// with the triangle 89,0.1;89,90.1;89,-179.9. The volatile declaration
// fixes this. (The bug probably triggered because transit and
// AngNormalize are inline functions. So don't port this change over to
// Geodesic.hpp.)
volatile real y = x;
return y >= 180 ? y - 360 : (y < -180 ? y + 360 : y);
}
static inline int transit(real lon1, real lon2) {
// Return 1 or -1 if crossing prime meridian in east or west direction.
// Otherwise return zero.
lon1 = AngNormalize(lon1);
lon2 = AngNormalize(lon2);
// treat lon12 = -180 as an eastward geodesic, so convert to 180.
real lon12 = -AngNormalize(lon1 - lon2); // In (-180, 180]
int cross =
lon1 < 0 && lon2 >= 0 && lon12 > 0 ? 1 :
(lon2 < 0 && lon1 >= 0 && lon12 < 0 ? -1 : 0);
return cross;
}
public:
/**
* Constructor for PolygonArea.
*
* @param[in] earth the Geodesic object to use for geodesic calculations.
* By default this uses the WGS84 ellipsoid.
* @param[in] polyline if true that treat the points as defining a polyline
* instead of a polygon (default = false).
**********************************************************************/
PolygonArea(const Geodesic& earth, bool polyline = false) throw()
: _earth(earth)
, _area0(_earth.EllipsoidArea())
, _polyline(polyline)
, _mask(Geodesic::DISTANCE | (_polyline ? 0 : Geodesic::AREA))
{
Clear();
}
/**
* Clear PolygonArea, allowing a new polygon to be started.
**********************************************************************/
void Clear() throw() {
_num = 0;
_crossings = 0;
_areasum = 0;
_perimetersum = 0;
_lat0 = _lon0 = _lat1 = _lon1 = 0;
}
/**
* Add a point to the polygon or polyline.
*
* @param[in] lat the latitude of the point (degrees).
* @param[in] lon the latitude of the point (degrees).
*
* \e lat should be in the range [-90, 90] and \e lon should be in the
* range [-180, 360].
**********************************************************************/
void AddPoint(real lat, real lon) throw();
/**
* Return the results so far.
*
* @param[in] reverse if true then clockwise (instead of counter-clockwise)
* traversal counts as a positive area.
* @param[in] sign if true 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] perimeter the perimeter of the polygon or length of the
* polyline (meters).
* @param[out] area the area of the polygon (meters^2); only set if
* polyline is false in the constructor.
* @return the number of points.
**********************************************************************/
unsigned Compute(bool reverse, bool sign,
real& perimeter, real& area) const throw();
/**
* 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 AddPoint and Compute are used.
*
* @param[in] lat the latitude of the test point (degrees).
* @param[in] lon the longitude of the test point (degrees).
* @param[in] reverse if true then clockwise (instead of counter-clockwise)
* traversal counts as a positive area.
* @param[in] sign if true 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] perimeter the approximate perimeter of the polygon or length
* of the polyline (meters).
* @param[out] area the approximate area of the polygon (meters^2); only
* set if polyline is false in the constructor.
* @return the number of points.
*
* \e lat should be in the range [-90, 90] and \e lon should be in the
* range [-180, 360].
**********************************************************************/
unsigned TestCompute(real lat, real lon, bool reverse, bool sign,
real& perimeter, real& area) const throw();
/** \name Inspector functions
**********************************************************************/
///@{
/**
* @return \e a the equatorial radius of the ellipsoid (meters). This is
* the value inherited from the Geodesic object used in the constructor.
**********************************************************************/
Math::real MajorRadius() const throw() { return _earth.MajorRadius(); }
/**
* @return \e f the flattening of the ellipsoid. This is the value
* inherited from the Geodesic object used in the constructor.
**********************************************************************/
Math::real Flattening() const throw() { return _earth.Flattening(); }
///@}
};
} // namespace GeographicLib
#endif // GEOGRAPHICLIB_POLYGONAREA_HPP
|