From 3f71775101ac809a8f3098d302ab447d2db95fb1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Martin=20T=C5=AFma?= Date: Fri, 13 Apr 2018 21:14:12 +0200 Subject: [PATCH] Fixed broken radius-rect computation Always use the double type for coordinate related data --- src/common/coordinates.cpp | 42 ++++---------------------------------- src/common/coordinates.h | 20 ++++++++---------- src/common/rectc.cpp | 39 ++++++++++++++++++++++++++++++++++- src/common/rectc.h | 7 +------ src/data/poi.cpp | 24 ++++++++++------------ src/map/mapfile.cpp | 8 ++++---- src/map/wms.cpp | 2 +- src/map/wmsmap.cpp | 14 ++++++------- src/map/wmsmap.h | 4 ++-- src/map/wmts.h | 8 ++++---- src/map/wmtsmap.cpp | 4 ++-- src/map/wmtsmap.h | 2 +- 12 files changed, 84 insertions(+), 90 deletions(-) diff --git a/src/common/coordinates.cpp b/src/common/coordinates.cpp index 12eeb2b9..b459ec31 100644 --- a/src/common/coordinates.cpp +++ b/src/common/coordinates.cpp @@ -1,50 +1,16 @@ #include "wgs84.h" #include "coordinates.h" -#define MIN_LAT deg2rad(-90.0) -#define MAX_LAT deg2rad(90.0) -#define MIN_LON deg2rad(-180.0) -#define MAX_LON deg2rad(180.0) - -qreal Coordinates::distanceTo(const Coordinates &c) const +double Coordinates::distanceTo(const Coordinates &c) const { - qreal dLat = deg2rad(c.lat() - _lat); - qreal dLon = deg2rad(c.lon() - _lon); - qreal a = pow(sin(dLat / 2.0), 2.0) + double dLat = deg2rad(c.lat() - _lat); + double dLon = deg2rad(c.lon() - _lon); + double a = pow(sin(dLat / 2.0), 2.0) + cos(deg2rad(_lat)) * cos(deg2rad(c.lat())) * pow(sin(dLon / 2.0), 2.0); return (WGS84_RADIUS * (2.0 * atan2(sqrt(a), sqrt(1.0 - a)))); } -QPair Coordinates::boundingRect(qreal distance) const -{ - qreal radDist = distance / WGS84_RADIUS; - - qreal minLat = deg2rad(_lat) - radDist; - qreal maxLat = deg2rad(_lat) + radDist; - - qreal minLon, maxLon; - if (minLat > MIN_LAT && maxLat < MAX_LAT) { - qreal deltaLon = asin(sin(radDist) / cos(_lat)); - minLon = deg2rad(_lon) - deltaLon; - if (minLon < MIN_LON) - minLon += 2.0 * M_PI; - maxLon = deg2rad(_lon) + deltaLon; - if (maxLon > MAX_LON) - maxLon -= 2.0 * M_PI; - } else { - // a pole is within the distance - minLat = qMax(minLat, MIN_LAT); - maxLat = qMin(maxLat, MAX_LAT); - minLon = MIN_LON; - maxLon = MAX_LON; - } - - return QPair(Coordinates(rad2deg(qMin(minLon, - maxLon)), rad2deg(qMin(minLat, maxLat))), Coordinates(rad2deg(qMax(minLon, - maxLon)), rad2deg(qMax(minLat, maxLat)))); -} - #ifndef QT_NO_DEBUG QDebug operator<<(QDebug dbg, const Coordinates &c) { diff --git a/src/common/coordinates.h b/src/common/coordinates.h index fed0390f..e35b314e 100644 --- a/src/common/coordinates.h +++ b/src/common/coordinates.h @@ -2,7 +2,6 @@ #define COORDINATES_H #include -#include #include #ifndef M_PI @@ -15,14 +14,14 @@ class Coordinates { public: Coordinates() {_lon = NAN; _lat = NAN;} - Coordinates(qreal lon, qreal lat) {_lon = lon; _lat = lat;} + Coordinates(double lon, double lat) {_lon = lon; _lat = lat;} - qreal &rlon() {return _lon;} - qreal &rlat() {return _lat;} - void setLon(qreal lon) {_lon = lon;} - void setLat(qreal lat) {_lat = lat;} - qreal lon() const {return _lon;} - qreal lat() const {return _lat;} + double &rlon() {return _lon;} + double &rlat() {return _lat;} + void setLon(double lon) {_lon = lon;} + void setLat(double lat) {_lat = lat;} + double lon() const {return _lon;} + double lat() const {return _lat;} bool isNull() const {return std::isnan(_lon) && std::isnan(_lat);} @@ -30,11 +29,10 @@ public: {return (_lon >= -180.0 && _lon <= 180.0 && _lat >= -90.0 && _lat <= 90.0);} - qreal distanceTo(const Coordinates &c) const; - QPair boundingRect(qreal distance) const; + double distanceTo(const Coordinates &c) const; private: - qreal _lat, _lon; + double _lat, _lon; }; inline bool operator==(const Coordinates &c1, const Coordinates &c2) diff --git a/src/common/rectc.cpp b/src/common/rectc.cpp index 8c0dd430..76a0a4c5 100644 --- a/src/common/rectc.cpp +++ b/src/common/rectc.cpp @@ -1,5 +1,41 @@ +#include "wgs84.h" #include "rectc.h" +#define MIN_LAT deg2rad(-90.0) +#define MAX_LAT deg2rad(90.0) +#define MIN_LON deg2rad(-180.0) +#define MAX_LON deg2rad(180.0) + +RectC::RectC(const Coordinates ¢er, double radius) +{ + double radDist = radius / WGS84_RADIUS; + double radLon = deg2rad(center.lon()); + double radlat = deg2rad(center.lat()); + + double minLat = radlat - radDist; + double maxLat = radlat + radDist; + + double minLon, maxLon; + if (minLat > MIN_LAT && maxLat < MAX_LAT) { + double deltaLon = asin(sin(radDist) / cos(radlat)); + minLon = radLon - deltaLon; + if (minLon < MIN_LON) + minLon += 2.0 * M_PI; + maxLon = radLon + deltaLon; + if (maxLon > MAX_LON) + maxLon -= 2.0 * M_PI; + } else { + // a pole is within the distance + minLat = qMax(minLat, MIN_LAT); + maxLat = qMin(maxLat, MAX_LAT); + minLon = MIN_LON; + maxLon = MAX_LON; + } + + _tl = Coordinates(rad2deg(minLon), rad2deg(maxLat)); + _br = Coordinates(rad2deg(maxLon), rad2deg(minLat)); +} + RectC RectC::operator|(const RectC &r) const { if (isNull()) @@ -134,7 +170,8 @@ void RectC::unite(const Coordinates &c) #ifndef QT_NO_DEBUG QDebug operator<<(QDebug dbg, const RectC &rect) { - dbg.nospace() << "RectC(" << rect.topLeft() << ", " << rect.size() << ")"; + dbg.nospace() << "RectC(" << rect.topLeft() << ", " << rect.bottomRight() + << ")"; return dbg.space(); } #endif // QT_NO_DEBUG diff --git a/src/common/rectc.h b/src/common/rectc.h index 3279aa73..ce61e14b 100644 --- a/src/common/rectc.h +++ b/src/common/rectc.h @@ -2,7 +2,6 @@ #define RECTC_H #include -#include #include "coordinates.h" class RectC @@ -11,6 +10,7 @@ public: RectC() {} RectC(const Coordinates &topLeft, const Coordinates &bottomRight) : _tl(topLeft), _br(bottomRight) {} + RectC(const Coordinates ¢er, double radius); bool isNull() const {return _tl.isNull() && _br.isNull();} @@ -19,14 +19,9 @@ public: Coordinates topLeft() const {return _tl;} Coordinates bottomRight() const {return _br;} - Coordinates center() const {return Coordinates((_tl.lon() + _br.lon()) / 2.0, (_tl.lat() + _br.lat()) / 2.0);} - qreal width() const {return _br.lon() - _tl.lon();} - qreal height() const {return _br.lat() - _tl.lat();} - - QSizeF size() const {return QSizeF(width(), height());} RectC operator|(const RectC &r) const; RectC &operator|=(const RectC &r) {*this = *this | r; return *this;} diff --git a/src/data/poi.cpp b/src/data/poi.cpp index fe2fd726..311916a1 100644 --- a/src/data/poi.cpp +++ b/src/data/poi.cpp @@ -1,5 +1,6 @@ #include #include +#include "common/rectc.h" #include "data.h" #include "poi.h" @@ -101,12 +102,11 @@ QList POI::points(const Path &path) const QSet::const_iterator it; for (int i = 0; i < path.count(); i++) { - const Coordinates &c = path.at(i).coordinates(); - QPair br = c.boundingRect(_radius); - min[0] = br.first.lon(); - min[1] = br.first.lat(); - max[0] = br.second.lon(); - max[1] = br.second.lat(); + RectC br(path.at(i).coordinates(), _radius); + min[0] = br.topLeft().lon(); + min[1] = br.bottomRight().lat(); + max[0] = br.bottomRight().lon(); + max[1] = br.topLeft().lat(); _tree.Search(min, max, cb, &set); } @@ -124,13 +124,11 @@ QList POI::points(const Waypoint &point) const qreal min[2], max[2]; QSet::const_iterator it; - const Coordinates &c = point.coordinates(); - - QPair br = c.boundingRect(_radius); - min[0] = br.first.lon(); - min[1] = br.first.lat(); - max[0] = br.second.lon(); - max[1] = br.second.lat(); + RectC br(point.coordinates(), _radius); + min[0] = br.topLeft().lon(); + min[1] = br.bottomRight().lat(); + max[0] = br.bottomRight().lon(); + max[1] = br.topLeft().lat(); _tree.Search(min, max, cb, &set); diff --git a/src/map/mapfile.cpp b/src/map/mapfile.cpp index abe5eca6..890a24c5 100644 --- a/src/map/mapfile.cpp +++ b/src/map/mapfile.cpp @@ -54,13 +54,13 @@ int MapFile::parse(QIODevice &device, QList &points, int latd = list.at(6).trimmed().toInt(&res); if (!res) ll = false; - qreal latm = list.at(7).trimmed().toFloat(&res); + double latm = list.at(7).trimmed().toDouble(&res); if (!res) ll = false; int lond = list.at(9).trimmed().toInt(&res); if (!res) ll = false; - qreal lonm = list.at(10).trimmed().toFloat(&res); + double lonm = list.at(10).trimmed().toDouble(&res); if (!res) ll = false; if (ll && list.at(8).trimmed() == "S") { @@ -75,10 +75,10 @@ int MapFile::parse(QIODevice &device, QList &points, p.zone = list.at(13).trimmed().toInt(&res); if (!res) p.zone = 0; - qreal ppx = list.at(14).trimmed().toFloat(&res); + double ppx = list.at(14).trimmed().toDouble(&res); if (!res) pp = false; - qreal ppy = list.at(15).trimmed().toFloat(&res); + double ppy = list.at(15).trimmed().toDouble(&res); if (!res) pp = false; if (list.at(16).trimmed() == "S") diff --git a/src/map/wms.cpp b/src/map/wms.cpp index d324cf72..0f94af56 100644 --- a/src/map/wms.cpp +++ b/src/map/wms.cpp @@ -63,7 +63,7 @@ QString WMS::style(QXmlStreamReader &reader) RectC WMS::geographicBoundingBox(QXmlStreamReader &reader) { - qreal left = NAN, top = NAN, right = NAN, bottom = NAN; + double left = NAN, top = NAN, right = NAN, bottom = NAN; while (reader.readNextStartElement()) { if (reader.name() == "westBoundLongitude") diff --git a/src/map/wmsmap.cpp b/src/map/wmsmap.cpp index 62d8a784..71059de4 100644 --- a/src/map/wmsmap.cpp +++ b/src/map/wmsmap.cpp @@ -10,7 +10,7 @@ #define CAPABILITIES_FILE "capabilities.xml" #define TILE_SIZE 256 -qreal WMSMap::sd2res(qreal scaleDenominator) const +double WMSMap::sd2res(double scaleDenominator) const { return scaleDenominator * 0.28e-3 * _projection.units().fromMeters(1.0); } @@ -48,9 +48,9 @@ void WMSMap::computeZooms(const RangeF &scaleDenominator) _zooms.clear(); if (scaleDenominator.size() > 0) { - qreal ld = log2(scaleDenominator.max()) - log2(scaleDenominator.min()); + double ld = log2(scaleDenominator.max()) - log2(scaleDenominator.min()); int cld = ceil(ld); - qreal step = ld / (qreal)cld; + double step = ld / (qreal)cld; qreal lmax = log2(scaleDenominator.max()); for (int i = 0; i <= cld; i++) _zooms.append(pow(2.0, lmax - i * step)); @@ -60,10 +60,10 @@ void WMSMap::computeZooms(const RangeF &scaleDenominator) void WMSMap::updateTransform() { - qreal scaleDenominator = _zooms.at(_zoom); + double scaleDenominator = _zooms.at(_zoom); ReferencePoint tl, br; - qreal pixelSpan = sd2res(scaleDenominator); + double pixelSpan = sd2res(scaleDenominator); if (_projection.isGeographic()) pixelSpan /= deg2rad(WGS84_RADIUS); @@ -150,7 +150,7 @@ void WMSMap::emitLoaded() QRectF WMSMap::bounds() const { - qreal pixelSpan = sd2res(_zooms.at(_zoom)); + double pixelSpan = sd2res(_zooms.at(_zoom)); if (_projection.isGeographic()) pixelSpan /= deg2rad(WGS84_RADIUS); QSizeF size(_boundingBox.width() / pixelSpan, -_boundingBox.height() @@ -176,7 +176,7 @@ int WMSMap::zoomFit(const QSize &size, const RectC &br) QRectF tbr(_projection.ll2xy(br.topLeft()), _projection.ll2xy(br.bottomRight())); QPointF sc(tbr.width() / size.width(), tbr.height() / size.height()); - qreal resolution = qMax(qAbs(sc.x()), qAbs(sc.y())); + double resolution = qMax(qAbs(sc.x()), qAbs(sc.y())); if (_projection.isGeographic()) resolution *= deg2rad(WGS84_RADIUS); diff --git a/src/map/wmsmap.h b/src/map/wmsmap.h index 1e2c2336..cf8a0307 100644 --- a/src/map/wmsmap.h +++ b/src/map/wmsmap.h @@ -46,7 +46,7 @@ private slots: private: QString tileUrl(const QString &version) const; - qreal sd2res(qreal scaleDenominator) const; + double sd2res(double scaleDenominator) const; QString tilesDir() const; void computeZooms(const RangeF &scaleDenominator); void updateTransform(); @@ -62,7 +62,7 @@ private: Projection _projection; Transform _transform; CoordinateSystem _cs; - QVector _zooms; + QVector _zooms; QRectF _boundingBox; int _zoom; bool _block; diff --git a/src/map/wmts.h b/src/map/wmts.h index 728e7cbe..ccdc1661 100644 --- a/src/map/wmts.h +++ b/src/map/wmts.h @@ -54,7 +54,7 @@ public: class Zoom { public: - Zoom(const QString &id, qreal scaleDenominator, const QPointF &topLeft, + Zoom(const QString &id, double scaleDenominator, const QPointF &topLeft, const QSize &tile, const QSize &matrix, const QRect &limits) : _id(id), _scaleDenominator(scaleDenominator), _topLeft(topLeft), _tile(tile), _matrix(matrix), _limits(limits) {} @@ -62,7 +62,7 @@ public: {return _scaleDenominator > other._scaleDenominator;} const QString &id() const {return _id;} - qreal scaleDenominator() const {return _scaleDenominator;} + double scaleDenominator() const {return _scaleDenominator;} const QPointF &topLeft() const {return _topLeft;} const QSize &tile() const {return _tile;} const QSize &matrix() const {return _matrix;} @@ -70,7 +70,7 @@ public: private: QString _id; - qreal _scaleDenominator; + double _scaleDenominator; QPointF _topLeft; QSize _tile; QSize _matrix; @@ -94,7 +94,7 @@ public: private: struct TileMatrix { QString id; - qreal scaleDenominator; + double scaleDenominator; QPointF topLeft; QSize tile; QSize matrix; diff --git a/src/map/wmtsmap.cpp b/src/map/wmtsmap.cpp index 38c06ba9..0b8019ae 100644 --- a/src/map/wmtsmap.cpp +++ b/src/map/wmtsmap.cpp @@ -62,7 +62,7 @@ QString WMTSMap::tilesDir() const return QString(TILES_DIR + "/" + _name); } -qreal WMTSMap::sd2res(qreal scaleDenominator) const +double WMTSMap::sd2res(double scaleDenominator) const { return scaleDenominator * 0.28e-3 * _projection.units().fromMeters(1.0); } @@ -75,7 +75,7 @@ void WMTSMap::updateTransform() QPointF topLeft = (_cs.axisOrder() == CoordinateSystem::YX) ? QPointF(z.topLeft().y(), z.topLeft().x()) : z.topLeft(); - qreal pixelSpan = sd2res(z.scaleDenominator()); + double pixelSpan = sd2res(z.scaleDenominator()); if (_projection.isGeographic()) pixelSpan /= deg2rad(WGS84_RADIUS); QPointF tileSpan(z.tile().width() * pixelSpan, z.tile().height() * pixelSpan); diff --git a/src/map/wmtsmap.h b/src/map/wmtsmap.h index 466bae00..a1c4c969 100644 --- a/src/map/wmtsmap.h +++ b/src/map/wmtsmap.h @@ -46,7 +46,7 @@ private slots: private: bool loadWMTS(); - qreal sd2res(qreal scaleDenominator) const; + double sd2res(double scaleDenominator) const; QString tilesDir() const; void updateTransform();