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
|
// Example of using the GeographicLib::NearestNeighbor class. WARNING: this
// creates a file, pointset.xml or pointset.txt, in the current directory.
// Read lat/lon locations from locations.txt and lat/lon queries from
// queries.txt. For each query print to standard output: the index for the
// closest location and the distance to it. Print statistics to standard error
// at the end.
#include <iostream>
#include <exception>
#include <vector>
#include <fstream>
#include <string>
#if !defined(GEOGRAPHICLIB_HAVE_BOOST_SERIALIZATION)
#define GEOGRAPHICLIB_HAVE_BOOST_SERIALIZATION 0
#endif
#if GEOGRAPHICLIB_HAVE_BOOST_SERIALIZATION
// If Boost serialization is available, use it.
#include <boost/archive/xml_iarchive.hpp>
#include <boost/archive/xml_oarchive.hpp>
#endif
#include <GeographicLib/NearestNeighbor.hpp>
#include <GeographicLib/Geodesic.hpp>
#include <GeographicLib/DMS.hpp>
using namespace std;
using namespace GeographicLib;
// A structure to hold a geographic coordinate.
struct pos {
double _lat, _lon;
pos(double lat = 0, double lon = 0) : _lat(lat), _lon(lon) {}
};
// A class to compute the distance between 2 positions.
class DistanceCalculator {
private:
Geodesic _geod;
public:
explicit DistanceCalculator(const Geodesic& geod) : _geod(geod) {}
double operator() (const pos& a, const pos& b) const {
double d;
_geod.Inverse(a._lat, a._lon, b._lat, b._lon, d);
if ( !(d >= 0) )
// Catch illegal positions which result in d = NaN
throw GeographicErr("distance doesn't satisfy d >= 0");
return d;
}
};
int main() {
try {
// Read in locations
vector<pos> locs;
double lat, lon;
string sa, sb;
{
ifstream is("locations.txt");
if (!is.good())
throw GeographicErr("locations.txt not readable");
while (is >> sa >> sb) {
DMS::DecodeLatLon(sa, sb, lat, lon);
locs.push_back(pos(lat, lon));
}
if (locs.size() == 0)
throw GeographicErr("need at least one location");
}
// Define a distance function object
DistanceCalculator distance(Geodesic::WGS84());
// Create NearestNeighbor object
NearestNeighbor<double, pos, DistanceCalculator> pointset;
{
// Used saved object if it is available
#if GEOGRAPHICLIB_HAVE_BOOST_SERIALIZATION
ifstream is("pointset.xml");
if (is.good()) {
boost::archive::xml_iarchive ia(is);
ia >> BOOST_SERIALIZATION_NVP(pointset);
}
#else
ifstream is("pointset.txt");
if (is.good())
is >> pointset;
#endif
}
// Is the saved pointset up-to-date?
if (pointset.NumPoints() != int(locs.size())) {
// else initialize it
pointset.Initialize(locs, distance);
// and save it
#if GEOGRAPHICLIB_HAVE_BOOST_SERIALIZATION
ofstream os("pointset.xml");
if (!os.good())
throw GeographicErr("cannot write to pointset.xml");
boost::archive::xml_oarchive oa(os);
oa << BOOST_SERIALIZATION_NVP(pointset);
#else
ofstream os("pointset.txt");
if (!os.good())
throw GeographicErr("cannot write to pointset.txt");
os << pointset << "\n";
#endif
}
ifstream is("queries.txt");
double d;
int count = 0;
vector<int> k;
while (is >> sa >> sb) {
++count;
DMS::DecodeLatLon(sa, sb, lat, lon);
d = pointset.Search(locs, distance, pos(lat, lon), k);
if (k.size() != 1)
throw GeographicErr("unexpected number of results");
cout << k[0] << " " << d << "\n";
}
int setupcost, numsearches, searchcost, mincost, maxcost;
double mean, sd;
pointset.Statistics(setupcost, numsearches, searchcost,
mincost, maxcost, mean, sd);
int
totcost = setupcost + searchcost,
exhaustivecost = count * pointset.NumPoints();
cerr
<< "Number of distance calculations = " << totcost << "\n"
<< "With an exhaustive search = " << exhaustivecost << "\n"
<< "Ratio = " << double(totcost) / exhaustivecost << "\n"
<< "Efficiency improvement = "
<< 100 * (1 - double(totcost) / exhaustivecost) << "%\n";
}
catch (const exception& e) {
cerr << "Caught exception: " << e.what() << "\n";
return 1;
}
}
|
