From 585e05997c07c4965a6ad91b7f8d6a7bf17addc6 Mon Sep 17 00:00:00 2001 From: tsteven4 <13596209+tsteven4@users.noreply.github.com> Date: Fri, 26 Nov 2021 07:36:57 -0800 Subject: [PATCH 1/9] Improve kml abstract view. Use welzl's algorithm and nvectors to find the latitude, longitude, and range for the camera. This works everywhere on earth, including the antimeridian and poles. --- CMakeLists.txt | 2 + GPSBabel.pro | 2 + kml.cc | 23 +++--- kml.h | 3 +- src/core/welzl.cc | 179 ++++++++++++++++++++++++++++++++++++++++++++++ src/core/welzl.h | 63 ++++++++++++++++ 6 files changed, 258 insertions(+), 14 deletions(-) create mode 100644 src/core/welzl.cc create mode 100644 src/core/welzl.h diff --git a/CMakeLists.txt b/CMakeLists.txt index 4e12d0013..80652376c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -100,6 +100,7 @@ set(SUPPORT src/core/textstream.cc src/core/usasciicodec.cc src/core/vector3d.cc + src/core/welzl.cc src/core/xmlstreamwriter.cc ) if(${QT_VERSION_MAJOR} EQUAL "6") @@ -193,6 +194,7 @@ set(HEADERS src/core/textstream.h src/core/usasciicodec.h src/core/vector3d.h + src/core/welzl.h src/core/xmlstreamwriter.h src/core/xmltag.h ) diff --git a/GPSBabel.pro b/GPSBabel.pro index 00bc8d03c..924bd3d2c 100644 --- a/GPSBabel.pro +++ b/GPSBabel.pro @@ -111,6 +111,7 @@ SUPPORT = route.cc waypt.cc filter_vecs.cc util.cc vecs.cc mkshort.cc \ src/core/textstream.cc \ src/core/usasciicodec.cc \ src/core/vector3d.cc \ + src/core/welzl.cc \ src/core/xmlstreamwriter.cc versionAtLeast(QT_VERSION, 6.0): SUPPORT += src/core/codecdevice.cc @@ -190,6 +191,7 @@ HEADERS = \ src/core/textstream.h \ src/core/usasciicodec.h \ src/core/vector3d.h \ + src/core/welzl.h \ src/core/xmlstreamwriter.h \ src/core/xmltag.h diff --git a/kml.cc b/kml.cc index 695b30d65..6880f2ee2 100644 --- a/kml.cc +++ b/kml.cc @@ -49,6 +49,7 @@ #include "src/core/datetime.h" // for DateTime #include "src/core/file.h" // for File #include "src/core/logging.h" // for Warning, Fatal +#include "src/core/welzl.h" #include "src/core/xmlstreamwriter.h" // for XmlStreamWriter #include "src/core/xmltag.h" // for xml_findfirst, xml_tag, fs_xml, xml_attribute, xml_findnext #include "units.h" // for fmt_setunits, fmt_speed, fmt_altitude, fmt_distance, units_aviation, units_metric, units_nautical, units_statute @@ -364,7 +365,7 @@ void KmlFormat::rd_deinit() void KmlFormat::wr_init(const QString& fname) { char u = 's'; - waypt_init_bounds(&kml_bounds); + kml_points.clear(); kml_time_min = QDateTime(); kml_time_max = QDateTime(); @@ -814,7 +815,7 @@ void KmlFormat::kml_recompute_time_bounds(const Waypoint* waypointp) void KmlFormat::kml_add_to_bounds(const Waypoint* waypointp) { - waypt_add_to_bounds(&kml_bounds, waypointp); + kml_points.append(gpsbabel::NVector(waypointp->latitude, waypointp->longitude)); kml_recompute_time_bounds(waypointp); } @@ -1682,6 +1683,8 @@ void KmlFormat::kml_write_AbstractView() route_disp_all(nullptr, nullptr, kml_add_to_bounds_lambda); } + gpsbabel::Circle bounds = gpsbabel::Welzl::welzl(kml_points); + writer->writeStartElement(QStringLiteral("LookAt")); if (kml_time_min.isValid() || kml_time_max.isValid()) { @@ -1704,20 +1707,14 @@ void KmlFormat::kml_write_AbstractView() writer->writeEndElement(); // Close gx:TimeSpan tag } -// If our BB spans the antemeridian, flip sign on one. -// This doesn't make our BB optimal, but it at least prevents us from -// zooming to the wrong hemisphere. - if (kml_bounds.min_lon * kml_bounds.max_lon < 0) { - kml_bounds.min_lon = -kml_bounds.max_lon; - } - - writer->writeTextElement(QStringLiteral("longitude"), QString::number((kml_bounds.min_lon + kml_bounds.max_lon) / 2, 'f', precision)); - writer->writeTextElement(QStringLiteral("latitude"), QString::number((kml_bounds.min_lat + kml_bounds.max_lat) / 2, 'f', precision)); + writer->writeTextElement(QStringLiteral("longitude"), QString::number(bounds.center().longitude(), 'f', precision)); + writer->writeTextElement(QStringLiteral("latitude"), QString::number(bounds.center().latitude(), 'f', precision)); // It turns out the length of the diagonal of the bounding box gives us a // reasonable guess for setting the camera altitude. - double bb_size = gcgeodist(kml_bounds.min_lat, kml_bounds.min_lon, - kml_bounds.max_lat, kml_bounds.max_lon); + // TODO: we can easily calculate a range such that the camera sees a bit beyond the points + // for the cases that the range is reasonable. + double bb_size = 2.0 * bounds.radius(); // Clamp bottom zoom level. Otherwise, a single point zooms to grass. if (bb_size < 1000) { bb_size = 1000; diff --git a/kml.h b/kml.h index d324b167f..855d82c98 100644 --- a/kml.h +++ b/kml.h @@ -33,6 +33,7 @@ #include "format.h" #include "src/core/datetime.h" // for DateTime #include "src/core/file.h" // for File +#include "src/core/nvector.h" #include "src/core/xmlstreamwriter.h" // for XmlStreamWriter #include "xmlgeneric.h" // for cb_cdata, cb_end, cb_start, xg_callback, xg_string, xg_cb_type, xml_deinit, xml_ignore_tags, xml_init, xml_read, xg_tag_mapping @@ -216,7 +217,7 @@ class KmlFormat : public Format gpsbabel::XmlStreamWriter* writer{nullptr}; int realtime_positioning{}; - bounds kml_bounds{}; + QVector kml_points; gpsbabel::DateTime kml_time_max; gpsbabel::DateTime kml_time_min; diff --git a/src/core/welzl.cc b/src/core/welzl.cc new file mode 100644 index 000000000..59c1b9c10 --- /dev/null +++ b/src/core/welzl.cc @@ -0,0 +1,179 @@ +/* + Copyright (C) 2021 Robert Lipe, robertlipe+source@gpsbabel.org + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program; if not, write to the Free Software + Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + + This solves the smallest circle problem using Welzl's algorithm. + + Welzl, Emo (1991), "Smallest enclosing disks (balls and ellipsoids)", + in Maurer, H. (ed.), New Results and New Trends in Computer Science, + Lecture Notes in Computer Science, 555, Springer-Verlag, pp. 359–370 + + */ + +#include "src/core/welzl.h" + +#include +#include +#include + +#include + +#include "src/core/nvector.h" +#include "src/core/vector3d.h" + + +namespace gpsbabel +{ + +//using gpsbabel::NVector; +//using gpsbabel::Vector3D; + + +Circle Welzl::b_md(QVector R) +{ + NVector center; + double radius; + + switch (R.size()) { + + case 0: + center = NVector(0.0, 0.0); + radius = 0.0; + break; + + case 1: + center = R.at(0); + radius = 0.0; + break; + + case 2: +#if 0 + double dp_a_b = NVector::dotProduct(R.at(0), R.at(1)); + if (dp_a_b <= (-1.0+8.0*DBL_EPSILON)) { + // The points are so close to be exactly opposite each other that + // there isn't a unique great circle between them. + Vector3D result(nan(""), nan(""), nan("")); + return result; + } +#endif + center = (R.at(0) + R.at(1)).normalize(); + radius = NVector::distance(center, R.at(0)); + //qDebug() << "center2 check" << NVector::distance(center, R.at(0)) << NVector::distance(center, R.at(1)); + break; + + case 3: + auto n_EA_E = R.at(0); + auto n_EB_E = R.at(1); + auto n_EC_E = R.at(2); + + double dp_a_b = NVector::dotProduct(n_EA_E, n_EB_E); + if (dp_a_b >= (1.0-8.0*DBL_EPSILON)) { + // The points are so close we will have trouble constructing a basis. + QVector Rprime{(n_EA_E + n_EB_E).normalize(), n_EC_E}; + return (b_md(Rprime)); + } + double dp_a_c = NVector::dotProduct(n_EA_E, n_EC_E); + if (dp_a_c >= (1.0-8.0*DBL_EPSILON)) { + // The points are so close we will have trouble constructing a basis. + QVector Rprime{(n_EA_E + n_EC_E).normalize(), n_EB_E}; + return (b_md(Rprime)); + } + double dp_b_c = NVector::dotProduct(n_EB_E, n_EC_E); + if (dp_b_c >= (1.0-8.0*DBL_EPSILON)) { + // The points are so close we will have trouble constructing a basis. + QVector Rprime{(n_EB_E + n_EC_E).normalize(), n_EA_E}; + return (b_md(Rprime)); + } + // Form an orthonormal basis with two components being in the plane + // of n_EA_E, n_EB_E and n_EC_E, and the third being perpendicular + // to this plane. + // Call this the W frame. + // The columns of the rotation matrix from E to W are w1, w2 and w3. + Vector3D w1 = (n_EA_E - n_EB_E).normalize(); + Vector3D w3 = NVector::crossProduct(w1, (n_EA_E - n_EC_E).normalize()).normalize(); + Vector3D w2 = NVector::crossProduct(w1, w3); + // Rotate A to the W frame. + Vector3D n_EA_W = Vector3D(w1*n_EA_E, w2*n_EA_E, w3*n_EA_E); + // The z component of n_EA_W, n_EB_W, n_EC_W are all identical. + // The point n_EX_W that is equidistant from n_EA_W, n_EB_W, n_EC_W is + // either (0.0, 0.0, 1.0) or (0.0, 0.0, -1.0). + NVector n_EX_W; + if (n_EA_W.getz() >= 0) { + n_EX_W = Vector3D(0.0, 0.0, 1.0); + } else { + n_EX_W = Vector3D(0.0, 0.0, -1.0); + } + // Translate n_EX_W back to the E frame. + // We need to invert the matrix composed of the column vectors w1, w2, w3. + // Since this matrix is orthogonal it's inverse equals it's transpose. + Vector3D wt1 = Vector3D(w1.getx(), w2.getx(), w3.getx()); + Vector3D wt2 = Vector3D(w1.gety(), w2.gety(), w3.gety()); + Vector3D wt3 = Vector3D(w1.getz(), w2.getz(), w3.getz()); + Vector3D n_EX_E = Vector3D(wt1*n_EX_W, wt2*n_EX_W, wt3*n_EX_W); + + center = n_EX_E; + radius = NVector::distance(center, R.at(0)); + + //qDebug() << R.at(0) << R.at(1) << R.at(2); + //qDebug() << center; + //qDebug() << "center3 check" << NVector::distance(center, R.at(0)) << NVector::distance(center, R.at(1)) << NVector::distance(center, R.at(2)); + break; + } + return Circle(center, radius); +}; + +bool Welzl::outside(Circle D, NVector n) +{ + return NVector::distance(D.center(), n) > D.radius(); +} + +/* This is Welzl's algorithm */ +Circle Welzl::b_mindisk(QVector P, QVector R) +{ + Circle D; + + if (P.isEmpty() || (R.size() == 3)) { + D = b_md(R); + } else { + auto p = P.last(); + auto Pprime(P); + Pprime.removeLast(); + D = b_mindisk(Pprime, R); + if (outside(D, p)) { + auto Pprime(P); + Pprime.removeLast(); + auto Rprime(R); + Rprime.append(p); + D = b_mindisk(Pprime, Rprime); + } + } + return D; +} + +// Return the center and radius of the smallest circle containg the points. +// This works across the antimeridian and at the poles. +Circle Welzl::welzl(QVector Points) +{ + // randomize order of Points + std::random_device rd; + std::mt19937 generator(rd()); + std::shuffle(Points.begin(), Points.end(), generator); + Circle bound = b_mindisk(Points, QVector()); + //qDebug() << bound.center_.latitude() << bound.center_.longitude() << bound.radius_; + return bound; +} + +} // namespace gpsbabel diff --git a/src/core/welzl.h b/src/core/welzl.h new file mode 100644 index 000000000..8a84b23bd --- /dev/null +++ b/src/core/welzl.h @@ -0,0 +1,63 @@ +/* + Copyright (C) 2021 Robert Lipe, robertlipe+source@gpsbabel.org + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program; if not, write to the Free Software + Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + + This solves the smallest circle problem using Welzl's algorithm. + + Welzl, Emo (1991), "Smallest enclosing disks (balls and ellipsoids)", + in Maurer, H. (ed.), New Results and New Trends in Computer Science, + Lecture Notes in Computer Science, 555, Springer-Verlag, pp. 359–370 + + */ + +#ifndef SRC_CORE_WELZL_H +#define SRC_CORE_WELZL_H + +#include +#include "src/core/nvector.h" + + +namespace gpsbabel +{ + +class Circle +{ +public: + Circle() = default; + Circle(NVector center, double radius) : center_(center), radius_(radius) {} + + NVector center() const {return center_;} + double radius() const {return radius_;} + +private: + NVector center_; + double radius_; +}; + +class Welzl +{ +public: + static Circle welzl(QVector); + +private: + static Circle b_md(QVector R); + static bool outside(Circle D, NVector n); + static Circle b_mindisk(QVector P, QVector R); +}; + +} // namespace gpsbabel + +#endif // SRC_CORE_WELZL_H From a727793119e9d6f823b37497dca417a0d7e2b388 Mon Sep 17 00:00:00 2001 From: tsteven4 <13596209+tsteven4@users.noreply.github.com> Date: Fri, 26 Nov 2021 08:58:06 -0700 Subject: [PATCH 2/9] update kml references for abstract view. There is a reproducibility issue with very small changes depeding on the randomization. --- reference/LineStyles.kml | 8 ++++---- reference/bounds-test.kml | 6 +++--- reference/earth-expertgps-track.kml | 6 +++--- reference/earth-expertgps.kml | 6 +++--- reference/earth-gc.kml | 16 +++++++-------- reference/realtime.kml | 2 +- reference/track/bounds-test-track.kml | 6 +++--- reference/track/gtrnctr_power-kml.kml | 24 +++++++++++----------- reference/track/segmented_tracks-track.kml | 4 ++-- reference/track/segmented_tracks.kml | 4 ++-- reference/track/tracks~gpx.kml | 6 +++--- 11 files changed, 44 insertions(+), 44 deletions(-) diff --git a/reference/LineStyles.kml b/reference/LineStyles.kml index 77c81ca37..2c46f6168 100644 --- a/reference/LineStyles.kml +++ b/reference/LineStyles.kml @@ -3,9 +3,9 @@ GPS device - -122.275670 - 37.523278 - 5908.813619 + -122.276195 + 37.523328 + 5678.471387 - + - + - + - Heart Rate + Heart Rate - Cadence + Cadence - Power + Power diff --git a/reference/track/segmented_tracks-track.kml b/reference/track/segmented_tracks-track.kml index 9e3e70ee2..89a7fe5b9 100644 --- a/reference/track/segmented_tracks-track.kml +++ b/reference/track/segmented_tracks-track.kml @@ -7,9 +7,9 @@ 2007-07-27T05:24:05Z 2007-07-27T05:35:00Z - -86.841461 + -86.841519 35.831217 - 1944.793167 + 1901.181635