00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015 #include "GeographicLib/Geodesic.hpp"
00016 #include "GeographicLib/DMS.hpp"
00017 #include "GeographicLib/GeoCoords.hpp"
00018 #include <iostream>
00019
00020 #include "Planimeter.usage"
00021
00022 int main(int argc, char* argv[]) {
00023 using namespace GeographicLib;
00024 typedef Math::real real;
00025
00026 class Accumulator {
00027
00028
00029
00030 private:
00031
00032
00033 real _s, _t;
00034
00035
00036 static inline real sum(real u, real v, real& t) {
00037 volatile real s = u + v;
00038 volatile real up = s - v;
00039 volatile real vpp = s - up;
00040 up -= u;
00041 vpp -= v;
00042 t = -(up + vpp);
00043
00044
00045 return s;
00046 }
00047 public:
00048 Accumulator() throw() : _s(0), _t(0) {};
00049 void Clear() throw() { _s = 0; _t = 0; }
00050
00051 void Add(real y) throw() {
00052 _s = sum(_s, y, y);
00053 _t += y;
00054 }
00055 void Negate() throw() { _s *= -1; _t *= -1; }
00056
00057 real Sum(real y) const throw() {
00058 real s = sum(_s, y, y);
00059 return s + (_t + y);
00060 }
00061
00062 real Sum() const throw() { return _s + _t; }
00063 };
00064
00065 class GeodesicPolygon {
00066 private:
00067 const Geodesic& _g;
00068 const real _area0;
00069 unsigned _num;
00070 int _crossings;
00071 Accumulator _area, _perimeter;
00072 real _lat0, _lon0, _lat1, _lon1;
00073
00074 static inline real AngNormalize(real x) throw() {
00075
00076 return x >= 180 ? x - 360 : x < -180 ? x + 360 : x;
00077 }
00078 static inline int transit(real lon1, real lon2) {
00079
00080
00081 lon1 = AngNormalize(lon1);
00082 lon2 = AngNormalize(lon2);
00083
00084 real lon12 = -AngNormalize(lon1 - lon2);
00085 int cross =
00086 lon1 < 0 && lon2 >= 0 && lon12 > 0 ? 1 :
00087 lon2 < 0 && lon1 >= 0 && lon12 < 0 ? -1 : 0;
00088 return cross;
00089 }
00090 public:
00091 GeodesicPolygon(const Geodesic& g) throw()
00092 : _g(g)
00093 , _area0(_g.EllipsoidArea())
00094 {
00095 Clear();
00096 }
00097 void Clear() throw() {
00098 _num = 0;
00099 _crossings = 0;
00100 _area.Clear();
00101 _perimeter.Clear();
00102 _lat0 = _lon0 = _lat1 = _lon1 = 0;
00103 }
00104 void AddPoint(real lat, real lon) throw() {
00105 if (_num == 0) {
00106 _lat0 = _lat1 = lat;
00107 _lon0 = _lon1 = lon;
00108 } else {
00109 real s12, S12, t;
00110 _g.GenInverse(_lat1, _lon1, lat, lon,
00111 Geodesic::DISTANCE | Geodesic::AREA,
00112 s12, t, t, t, t, t, S12);
00113 _perimeter.Add(s12);
00114 _area.Add(S12);
00115 if (_area.Sum() > _area0/2)
00116 _area.Add(-_area0);
00117 else if (_area.Sum() <= -_area0/2)
00118 _area.Add(_area0);
00119 _crossings += transit(_lon1, lon);
00120 _lat1 = lat;
00121 _lon1 = lon;
00122 }
00123 ++_num;
00124 }
00125 unsigned Compute(bool reverse, bool sign,
00126 real& perimeter, real& area) const throw() {
00127 real s12, S12, t;
00128 if (_num < 2) {
00129 perimeter = area = 0;
00130 return _num;
00131 }
00132 _g.GenInverse(_lat1, _lon1, _lat0, _lon0,
00133 Geodesic::DISTANCE | Geodesic::AREA,
00134 s12, t, t, t, t, t, S12);
00135 perimeter = _perimeter.Sum(s12);
00136 Accumulator area1(_area);
00137 area1.Add(S12);
00138 if (area1.Sum() > _area0/2)
00139 area1.Add(-_area0);
00140 else if (area1.Sum() <= -_area0/2)
00141 area1.Add(_area0);
00142 int crossings = _crossings + transit(_lon1, _lon0);
00143 if (crossings & 1) {
00144 if (area1.Sum() < 0)
00145 area1.Sum(_area0/2);
00146 else
00147 area1.Sum(-_area0/2);
00148 }
00149
00150
00151 if (!reverse)
00152 area1.Negate();
00153
00154 if (sign) {
00155 if (area1.Sum() > _area0/2)
00156 area1.Add(-_area0);
00157 else if (area1.Sum() <= -_area0/2)
00158 area1.Add(_area0);
00159 } else {
00160 if (area1.Sum() >= _area0)
00161 area1.Add(-_area0);
00162 else if (area1.Sum() < 0)
00163 area1.Add(_area0);
00164 }
00165 area = area1.Sum();
00166 return _num;
00167 }
00168 };
00169
00170 real
00171 a = Constants::WGS84_a<real>(),
00172 r = Constants::WGS84_r<real>();
00173 bool reverse = false, sign = false;
00174 for (int m = 1; m < argc; ++m) {
00175 std::string arg(argv[m]);
00176 if (arg == "-r")
00177 reverse = !reverse;
00178 else if (arg == "-s")
00179 sign = !sign;
00180 else if (arg == "-e") {
00181 if (m + 2 >= argc) return usage(1, true);
00182 try {
00183 a = DMS::Decode(std::string(argv[m + 1]));
00184 r = DMS::Decode(std::string(argv[m + 2]));
00185 }
00186 catch (const std::exception& e) {
00187 std::cerr << "Error decoding arguments of -e: " << e.what() << "\n";
00188 return 1;
00189 }
00190 m += 2;
00191 } else if (arg == "--version") {
00192 std::cout
00193 << PROGRAM_NAME
00194 << ": $Id: Planimeter.cpp 6978 2011-02-21 22:42:11Z karney $\n"
00195 << "GeographicLib version " << GEOGRAPHICLIB_VERSION << "\n";
00196 return 0;
00197 } else
00198 return usage(!(arg == "-h" || arg == "--help"), arg != "--help");
00199 }
00200
00201 const Geodesic geod(a, r);
00202 GeodesicPolygon poly(geod);
00203 GeoCoords p;
00204
00205 std::string s;
00206 real perimeter, area;
00207 unsigned num;
00208 while (std::getline(std::cin, s)) {
00209 try {
00210 p.Reset(s);
00211 if (p.Latitude() != p.Latitude() || p.Longitude() != p.Longitude())
00212 throw GeographicErr("NAN");
00213 }
00214 catch (const GeographicErr&) {
00215 num = poly.Compute(reverse, sign, perimeter, area);
00216 if (num > 0)
00217 std::cout << num << " "
00218 << DMS::Encode(perimeter, 8, DMS::NUMBER) << " "
00219 << DMS::Encode(area, 4, DMS::NUMBER) << "\n";
00220 poly.Clear();
00221 continue;
00222 }
00223 poly.AddPoint(p.Latitude(), p.Longitude());
00224 }
00225 num = poly.Compute(reverse, sign, perimeter, area);
00226 if (num > 0)
00227 std::cout << num << " "
00228 << DMS::Encode(perimeter, 8, DMS::NUMBER) << " "
00229 << DMS::Encode(area, 4, DMS::NUMBER) << "\n";
00230 poly.Clear();
00231 return 0;
00232 }