example-NearestNeighbor.cpp
Go to the documentation of this file.
1 // Example of using the GeographicLib::NearestNeighbor class. WARNING: this
2 // creates a file, pointset.xml or pointset.txt, in the current directory.
3 // Read lat/lon locations from locations.txt and lat/lon queries from
4 // queries.txt. For each query print to standard output: the index for the
5 // closest location and the distance to it. Print statistics to standard error
6 // at the end.
7 
8 #include <iostream>
9 #include <exception>
10 #include <vector>
11 #include <fstream>
12 #include <string>
13 
14 #if !defined(GEOGRAPHICLIB_HAVE_BOOST_SERIALIZATION)
15 #define GEOGRAPHICLIB_HAVE_BOOST_SERIALIZATION 0
16 #endif
17 
18 #if GEOGRAPHICLIB_HAVE_BOOST_SERIALIZATION
19 // If Boost serialization is available, use it.
20 #include <boost/archive/xml_iarchive.hpp>
21 #include <boost/archive/xml_oarchive.hpp>
22 #endif
23 
26 #include <GeographicLib/DMS.hpp>
27 
28 using namespace std;
29 using namespace GeographicLib;
30 
31 // A structure to hold a geographic coordinate.
32 struct pos {
33  double _lat, _lon;
34  pos(double lat = 0, double lon = 0) : _lat(lat), _lon(lon) {}
35 };
36 
37 // A class to compute the distance between 2 positions.
39 private:
41 public:
42  explicit DistanceCalculator(const Geodesic& geod) : _geod(geod) {}
43  double operator() (const pos& a, const pos& b) const {
44  double d;
45  _geod.Inverse(a._lat, a._lon, b._lat, b._lon, d);
46  if ( !(d >= 0) )
47  // Catch illegal positions which result in d = NaN
48  throw GeographicErr("distance doesn't satisfy d >= 0");
49  return d;
50  }
51 };
52 
53 int main() {
54  try {
55  // Read in locations
56  vector<pos> locs;
57  double lat, lon;
58  string sa, sb;
59  {
60  ifstream is("locations.txt");
61  if (!is.good())
62  throw GeographicErr("locations.txt not readable");
63  while (is >> sa >> sb) {
64  DMS::DecodeLatLon(sa, sb, lat, lon);
65  locs.push_back(pos(lat, lon));
66  }
67  if (locs.size() == 0)
68  throw GeographicErr("need at least one location");
69  }
70 
71  // Define a distance function object
72  DistanceCalculator distance(Geodesic::WGS84());
73 
74  // Create NearestNeighbor object
76 
77  {
78  // Used saved object if it is available
79 #if GEOGRAPHICLIB_HAVE_BOOST_SERIALIZATION
80  ifstream is("pointset.xml");
81  if (is.good()) {
82  boost::archive::xml_iarchive ia(is);
83  ia >> BOOST_SERIALIZATION_NVP(pointset);
84  }
85 #else
86  ifstream is("pointset.txt");
87  if (is.good())
88  is >> pointset;
89 #endif
90  }
91  // Is the saved pointset up-to-date?
92  if (pointset.NumPoints() != int(locs.size())) {
93  // else initialize it
94  pointset.Initialize(locs, distance);
95  // and save it
96 #if GEOGRAPHICLIB_HAVE_BOOST_SERIALIZATION
97  ofstream os("pointset.xml");
98  if (!os.good())
99  throw GeographicErr("cannot write to pointset.xml");
100  boost::archive::xml_oarchive oa(os);
101  oa << BOOST_SERIALIZATION_NVP(pointset);
102 #else
103  ofstream os("pointset.txt");
104  if (!os.good())
105  throw GeographicErr("cannot write to pointset.txt");
106  os << pointset << "\n";
107 #endif
108  }
109 
110  ifstream is("queries.txt");
111  double d;
112  int count = 0;
113  vector<int> k;
114  while (is >> sa >> sb) {
115  ++count;
116  DMS::DecodeLatLon(sa, sb, lat, lon);
117  d = pointset.Search(locs, distance, pos(lat, lon), k);
118  if (k.size() != 1)
119  throw GeographicErr("unexpected number of results");
120  cout << k[0] << " " << d << "\n";
121  }
122  int setupcost, numsearches, searchcost, mincost, maxcost;
123  double mean, sd;
124  pointset.Statistics(setupcost, numsearches, searchcost,
125  mincost, maxcost, mean, sd);
126  int
127  totcost = setupcost + searchcost,
128  exhaustivecost = count * pointset.NumPoints();
129  cerr
130  << "Number of distance calculations = " << totcost << "\n"
131  << "With an exhaustive search = " << exhaustivecost << "\n"
132  << "Ratio = " << double(totcost) / exhaustivecost << "\n"
133  << "Efficiency improvement = "
134  << 100 * (1 - double(totcost) / exhaustivecost) << "%\n";
135  }
136  catch (const exception& e) {
137  cerr << "Caught exception: " << e.what() << "\n";
138  return 1;
139  }
140 }
DistanceCalculator(const Geodesic &geod)
Header for GeographicLib::NearestNeighbor class.
Scalar * b
Definition: benchVecAdd.cpp:17
static const double lat
pos(double lat=0, double lon=0)
Double_ distance(const OrientedPlane3_ &p)
Definition: BFloat16.h:88
static const Vector2 mean(20, 40)
Header for GeographicLib::Geodesic class.
void Statistics(int &setupcost, int &numsearches, int &searchcost, int &mincost, int &maxcost, double &mean, double &sd) const
dist_t Search(const std::vector< pos_t > &pts, const distfun_t &dist, const pos_t &query, std::vector< int > &ind, int k=1, dist_t maxdist=std::numeric_limits< dist_t >::max(), dist_t mindist=-1, bool exhaustive=true, dist_t tol=0) const
Namespace for GeographicLib.
Array< double, 1, 3 > e(1./3., 0.5, 2.)
Exception handling for GeographicLib.
Definition: Constants.hpp:389
ofstream os("timeSchurFactors.csv")
void Initialize(const std::vector< pos_t > &pts, const distfun_t &dist, int bucket=4)
static const double lon
Nearest-neighbor calculations.
internal::enable_if< internal::valid_indexed_view_overload< RowIndices, ColIndices >::value &&internal::traits< typename EIGEN_INDEXED_VIEW_METHOD_TYPE< RowIndices, ColIndices >::type >::ReturnAsIndexedView, typename EIGEN_INDEXED_VIEW_METHOD_TYPE< RowIndices, ColIndices >::type >::type operator()(const RowIndices &rowIndices, const ColIndices &colIndices) EIGEN_INDEXED_VIEW_METHOD_CONST
Geodesic calculations
Definition: Geodesic.hpp:172
Math::real Inverse(real lat1, real lon1, real lat2, real lon2, real &s12, real &azi1, real &azi2, real &m12, real &M12, real &M21, real &S12) const
Definition: Geodesic.hpp:674
Header for GeographicLib::DMS class.


gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:34:13