1
0
mirror of https://github.com/tumic0/GPXSee.git synced 2024-11-30 22:51:16 +01:00

Fixed broken radius-rect computation

Always use the double type for coordinate related data
This commit is contained in:
Martin Tůma 2018-04-13 21:14:12 +02:00
parent f8009b0151
commit 3f71775101
12 changed files with 84 additions and 90 deletions

View File

@ -1,50 +1,16 @@
#include "wgs84.h" #include "wgs84.h"
#include "coordinates.h" #include "coordinates.h"
#define MIN_LAT deg2rad(-90.0) double Coordinates::distanceTo(const Coordinates &c) const
#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
{ {
qreal dLat = deg2rad(c.lat() - _lat); double dLat = deg2rad(c.lat() - _lat);
qreal dLon = deg2rad(c.lon() - _lon); double dLon = deg2rad(c.lon() - _lon);
qreal a = pow(sin(dLat / 2.0), 2.0) double a = pow(sin(dLat / 2.0), 2.0)
+ cos(deg2rad(_lat)) * cos(deg2rad(c.lat())) * pow(sin(dLon / 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)))); return (WGS84_RADIUS * (2.0 * atan2(sqrt(a), sqrt(1.0 - a))));
} }
QPair<Coordinates, Coordinates> 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, Coordinates>(Coordinates(rad2deg(qMin(minLon,
maxLon)), rad2deg(qMin(minLat, maxLat))), Coordinates(rad2deg(qMax(minLon,
maxLon)), rad2deg(qMax(minLat, maxLat))));
}
#ifndef QT_NO_DEBUG #ifndef QT_NO_DEBUG
QDebug operator<<(QDebug dbg, const Coordinates &c) QDebug operator<<(QDebug dbg, const Coordinates &c)
{ {

View File

@ -2,7 +2,6 @@
#define COORDINATES_H #define COORDINATES_H
#include <cmath> #include <cmath>
#include <QPair>
#include <QDebug> #include <QDebug>
#ifndef M_PI #ifndef M_PI
@ -15,14 +14,14 @@ class Coordinates
{ {
public: public:
Coordinates() {_lon = NAN; _lat = NAN;} 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;} double &rlon() {return _lon;}
qreal &rlat() {return _lat;} double &rlat() {return _lat;}
void setLon(qreal lon) {_lon = lon;} void setLon(double lon) {_lon = lon;}
void setLat(qreal lat) {_lat = lat;} void setLat(double lat) {_lat = lat;}
qreal lon() const {return _lon;} double lon() const {return _lon;}
qreal lat() const {return _lat;} double lat() const {return _lat;}
bool isNull() const bool isNull() const
{return std::isnan(_lon) && std::isnan(_lat);} {return std::isnan(_lon) && std::isnan(_lat);}
@ -30,11 +29,10 @@ public:
{return (_lon >= -180.0 && _lon <= 180.0 {return (_lon >= -180.0 && _lon <= 180.0
&& _lat >= -90.0 && _lat <= 90.0);} && _lat >= -90.0 && _lat <= 90.0);}
qreal distanceTo(const Coordinates &c) const; double distanceTo(const Coordinates &c) const;
QPair<Coordinates, Coordinates> boundingRect(qreal distance) const;
private: private:
qreal _lat, _lon; double _lat, _lon;
}; };
inline bool operator==(const Coordinates &c1, const Coordinates &c2) inline bool operator==(const Coordinates &c1, const Coordinates &c2)

View File

@ -1,5 +1,41 @@
#include "wgs84.h"
#include "rectc.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 &center, 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 RectC RectC::operator|(const RectC &r) const
{ {
if (isNull()) if (isNull())
@ -134,7 +170,8 @@ void RectC::unite(const Coordinates &c)
#ifndef QT_NO_DEBUG #ifndef QT_NO_DEBUG
QDebug operator<<(QDebug dbg, const RectC &rect) QDebug operator<<(QDebug dbg, const RectC &rect)
{ {
dbg.nospace() << "RectC(" << rect.topLeft() << ", " << rect.size() << ")"; dbg.nospace() << "RectC(" << rect.topLeft() << ", " << rect.bottomRight()
<< ")";
return dbg.space(); return dbg.space();
} }
#endif // QT_NO_DEBUG #endif // QT_NO_DEBUG

View File

@ -2,7 +2,6 @@
#define RECTC_H #define RECTC_H
#include <QDebug> #include <QDebug>
#include <QSizeF>
#include "coordinates.h" #include "coordinates.h"
class RectC class RectC
@ -11,6 +10,7 @@ public:
RectC() {} RectC() {}
RectC(const Coordinates &topLeft, const Coordinates &bottomRight) RectC(const Coordinates &topLeft, const Coordinates &bottomRight)
: _tl(topLeft), _br(bottomRight) {} : _tl(topLeft), _br(bottomRight) {}
RectC(const Coordinates &center, double radius);
bool isNull() const bool isNull() const
{return _tl.isNull() && _br.isNull();} {return _tl.isNull() && _br.isNull();}
@ -19,14 +19,9 @@ public:
Coordinates topLeft() const {return _tl;} Coordinates topLeft() const {return _tl;}
Coordinates bottomRight() const {return _br;} Coordinates bottomRight() const {return _br;}
Coordinates center() const Coordinates center() const
{return Coordinates((_tl.lon() + _br.lon()) / 2.0, {return Coordinates((_tl.lon() + _br.lon()) / 2.0,
(_tl.lat() + _br.lat()) / 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) const;
RectC &operator|=(const RectC &r) {*this = *this | r; return *this;} RectC &operator|=(const RectC &r) {*this = *this | r; return *this;}

View File

@ -1,5 +1,6 @@
#include <QFile> #include <QFile>
#include <QDir> #include <QDir>
#include "common/rectc.h"
#include "data.h" #include "data.h"
#include "poi.h" #include "poi.h"
@ -101,12 +102,11 @@ QList<Waypoint> POI::points(const Path &path) const
QSet<int>::const_iterator it; QSet<int>::const_iterator it;
for (int i = 0; i < path.count(); i++) { for (int i = 0; i < path.count(); i++) {
const Coordinates &c = path.at(i).coordinates(); RectC br(path.at(i).coordinates(), _radius);
QPair<Coordinates, Coordinates> br = c.boundingRect(_radius); min[0] = br.topLeft().lon();
min[0] = br.first.lon(); min[1] = br.bottomRight().lat();
min[1] = br.first.lat(); max[0] = br.bottomRight().lon();
max[0] = br.second.lon(); max[1] = br.topLeft().lat();
max[1] = br.second.lat();
_tree.Search(min, max, cb, &set); _tree.Search(min, max, cb, &set);
} }
@ -124,13 +124,11 @@ QList<Waypoint> POI::points(const Waypoint &point) const
qreal min[2], max[2]; qreal min[2], max[2];
QSet<int>::const_iterator it; QSet<int>::const_iterator it;
const Coordinates &c = point.coordinates(); RectC br(point.coordinates(), _radius);
min[0] = br.topLeft().lon();
QPair<Coordinates, Coordinates> br = c.boundingRect(_radius); min[1] = br.bottomRight().lat();
min[0] = br.first.lon(); max[0] = br.bottomRight().lon();
min[1] = br.first.lat(); max[1] = br.topLeft().lat();
max[0] = br.second.lon();
max[1] = br.second.lat();
_tree.Search(min, max, cb, &set); _tree.Search(min, max, cb, &set);

View File

@ -54,13 +54,13 @@ int MapFile::parse(QIODevice &device, QList<CalibrationPoint> &points,
int latd = list.at(6).trimmed().toInt(&res); int latd = list.at(6).trimmed().toInt(&res);
if (!res) if (!res)
ll = false; ll = false;
qreal latm = list.at(7).trimmed().toFloat(&res); double latm = list.at(7).trimmed().toDouble(&res);
if (!res) if (!res)
ll = false; ll = false;
int lond = list.at(9).trimmed().toInt(&res); int lond = list.at(9).trimmed().toInt(&res);
if (!res) if (!res)
ll = false; ll = false;
qreal lonm = list.at(10).trimmed().toFloat(&res); double lonm = list.at(10).trimmed().toDouble(&res);
if (!res) if (!res)
ll = false; ll = false;
if (ll && list.at(8).trimmed() == "S") { if (ll && list.at(8).trimmed() == "S") {
@ -75,10 +75,10 @@ int MapFile::parse(QIODevice &device, QList<CalibrationPoint> &points,
p.zone = list.at(13).trimmed().toInt(&res); p.zone = list.at(13).trimmed().toInt(&res);
if (!res) if (!res)
p.zone = 0; p.zone = 0;
qreal ppx = list.at(14).trimmed().toFloat(&res); double ppx = list.at(14).trimmed().toDouble(&res);
if (!res) if (!res)
pp = false; pp = false;
qreal ppy = list.at(15).trimmed().toFloat(&res); double ppy = list.at(15).trimmed().toDouble(&res);
if (!res) if (!res)
pp = false; pp = false;
if (list.at(16).trimmed() == "S") if (list.at(16).trimmed() == "S")

View File

@ -63,7 +63,7 @@ QString WMS::style(QXmlStreamReader &reader)
RectC WMS::geographicBoundingBox(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()) { while (reader.readNextStartElement()) {
if (reader.name() == "westBoundLongitude") if (reader.name() == "westBoundLongitude")

View File

@ -10,7 +10,7 @@
#define CAPABILITIES_FILE "capabilities.xml" #define CAPABILITIES_FILE "capabilities.xml"
#define TILE_SIZE 256 #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); return scaleDenominator * 0.28e-3 * _projection.units().fromMeters(1.0);
} }
@ -48,9 +48,9 @@ void WMSMap::computeZooms(const RangeF &scaleDenominator)
_zooms.clear(); _zooms.clear();
if (scaleDenominator.size() > 0) { if (scaleDenominator.size() > 0) {
qreal ld = log2(scaleDenominator.max()) - log2(scaleDenominator.min()); double ld = log2(scaleDenominator.max()) - log2(scaleDenominator.min());
int cld = ceil(ld); int cld = ceil(ld);
qreal step = ld / (qreal)cld; double step = ld / (qreal)cld;
qreal lmax = log2(scaleDenominator.max()); qreal lmax = log2(scaleDenominator.max());
for (int i = 0; i <= cld; i++) for (int i = 0; i <= cld; i++)
_zooms.append(pow(2.0, lmax - i * step)); _zooms.append(pow(2.0, lmax - i * step));
@ -60,10 +60,10 @@ void WMSMap::computeZooms(const RangeF &scaleDenominator)
void WMSMap::updateTransform() void WMSMap::updateTransform()
{ {
qreal scaleDenominator = _zooms.at(_zoom); double scaleDenominator = _zooms.at(_zoom);
ReferencePoint tl, br; ReferencePoint tl, br;
qreal pixelSpan = sd2res(scaleDenominator); double pixelSpan = sd2res(scaleDenominator);
if (_projection.isGeographic()) if (_projection.isGeographic())
pixelSpan /= deg2rad(WGS84_RADIUS); pixelSpan /= deg2rad(WGS84_RADIUS);
@ -150,7 +150,7 @@ void WMSMap::emitLoaded()
QRectF WMSMap::bounds() const QRectF WMSMap::bounds() const
{ {
qreal pixelSpan = sd2res(_zooms.at(_zoom)); double pixelSpan = sd2res(_zooms.at(_zoom));
if (_projection.isGeographic()) if (_projection.isGeographic())
pixelSpan /= deg2rad(WGS84_RADIUS); pixelSpan /= deg2rad(WGS84_RADIUS);
QSizeF size(_boundingBox.width() / pixelSpan, -_boundingBox.height() 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()), QRectF tbr(_projection.ll2xy(br.topLeft()),
_projection.ll2xy(br.bottomRight())); _projection.ll2xy(br.bottomRight()));
QPointF sc(tbr.width() / size.width(), tbr.height() / size.height()); 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()) if (_projection.isGeographic())
resolution *= deg2rad(WGS84_RADIUS); resolution *= deg2rad(WGS84_RADIUS);

View File

@ -46,7 +46,7 @@ private slots:
private: private:
QString tileUrl(const QString &version) const; QString tileUrl(const QString &version) const;
qreal sd2res(qreal scaleDenominator) const; double sd2res(double scaleDenominator) const;
QString tilesDir() const; QString tilesDir() const;
void computeZooms(const RangeF &scaleDenominator); void computeZooms(const RangeF &scaleDenominator);
void updateTransform(); void updateTransform();
@ -62,7 +62,7 @@ private:
Projection _projection; Projection _projection;
Transform _transform; Transform _transform;
CoordinateSystem _cs; CoordinateSystem _cs;
QVector<qreal> _zooms; QVector<double> _zooms;
QRectF _boundingBox; QRectF _boundingBox;
int _zoom; int _zoom;
bool _block; bool _block;

View File

@ -54,7 +54,7 @@ public:
class Zoom class Zoom
{ {
public: 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) : const QSize &tile, const QSize &matrix, const QRect &limits) :
_id(id), _scaleDenominator(scaleDenominator), _topLeft(topLeft), _id(id), _scaleDenominator(scaleDenominator), _topLeft(topLeft),
_tile(tile), _matrix(matrix), _limits(limits) {} _tile(tile), _matrix(matrix), _limits(limits) {}
@ -62,7 +62,7 @@ public:
{return _scaleDenominator > other._scaleDenominator;} {return _scaleDenominator > other._scaleDenominator;}
const QString &id() const {return _id;} const QString &id() const {return _id;}
qreal scaleDenominator() const {return _scaleDenominator;} double scaleDenominator() const {return _scaleDenominator;}
const QPointF &topLeft() const {return _topLeft;} const QPointF &topLeft() const {return _topLeft;}
const QSize &tile() const {return _tile;} const QSize &tile() const {return _tile;}
const QSize &matrix() const {return _matrix;} const QSize &matrix() const {return _matrix;}
@ -70,7 +70,7 @@ public:
private: private:
QString _id; QString _id;
qreal _scaleDenominator; double _scaleDenominator;
QPointF _topLeft; QPointF _topLeft;
QSize _tile; QSize _tile;
QSize _matrix; QSize _matrix;
@ -94,7 +94,7 @@ public:
private: private:
struct TileMatrix { struct TileMatrix {
QString id; QString id;
qreal scaleDenominator; double scaleDenominator;
QPointF topLeft; QPointF topLeft;
QSize tile; QSize tile;
QSize matrix; QSize matrix;

View File

@ -62,7 +62,7 @@ QString WMTSMap::tilesDir() const
return QString(TILES_DIR + "/" + _name); 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); return scaleDenominator * 0.28e-3 * _projection.units().fromMeters(1.0);
} }
@ -75,7 +75,7 @@ void WMTSMap::updateTransform()
QPointF topLeft = (_cs.axisOrder() == CoordinateSystem::YX) QPointF topLeft = (_cs.axisOrder() == CoordinateSystem::YX)
? QPointF(z.topLeft().y(), z.topLeft().x()) : z.topLeft(); ? QPointF(z.topLeft().y(), z.topLeft().x()) : z.topLeft();
qreal pixelSpan = sd2res(z.scaleDenominator()); double pixelSpan = sd2res(z.scaleDenominator());
if (_projection.isGeographic()) if (_projection.isGeographic())
pixelSpan /= deg2rad(WGS84_RADIUS); pixelSpan /= deg2rad(WGS84_RADIUS);
QPointF tileSpan(z.tile().width() * pixelSpan, z.tile().height() * pixelSpan); QPointF tileSpan(z.tile().width() * pixelSpan, z.tile().height() * pixelSpan);

View File

@ -46,7 +46,7 @@ private slots:
private: private:
bool loadWMTS(); bool loadWMTS();
qreal sd2res(qreal scaleDenominator) const; double sd2res(double scaleDenominator) const;
QString tilesDir() const; QString tilesDir() const;
void updateTransform(); void updateTransform();