1
0
mirror of https://github.com/tumic0/GPXSee.git synced 2024-11-27 21:24:47 +01:00

Fixed GeoTIFF on platforms where qreal != double

This commit is contained in:
Martin Tůma 2018-04-15 16:27:47 +02:00
parent 4b776c8cc1
commit f2b72ec1b5
27 changed files with 192 additions and 156 deletions

View File

@ -92,6 +92,20 @@ HEADERS += src/config.h \
src/map/transform.h \ src/map/transform.h \
src/map/mapfile.h \ src/map/mapfile.h \
src/map/tifffile.h \ src/map/tifffile.h \
src/map/gcs.h \
src/map/angularunits.h \
src/map/primemeridian.h \
src/map/linearunits.h \
src/map/ct.h \
src/map/mapsource.h \
src/map/tileloader.h \
src/map/wmtsmap.h \
src/map/wmts.h \
src/map/wmsmap.h \
src/map/wms.h \
src/map/crs.h \
src/map/coordinatesystem.h \
src/map/point.h \
src/data/graph.h \ src/data/graph.h \
src/data/poi.h \ src/data/poi.h \
src/data/waypoint.h \ src/data/waypoint.h \
@ -111,19 +125,6 @@ HEADERS += src/config.h \
src/data/fitparser.h \ src/data/fitparser.h \
src/data/igcparser.h \ src/data/igcparser.h \
src/data/nmeaparser.h \ src/data/nmeaparser.h \
src/map/gcs.h \
src/map/angularunits.h \
src/map/primemeridian.h \
src/map/linearunits.h \
src/map/ct.h \
src/map/mapsource.h \
src/map/tileloader.h \
src/map/wmtsmap.h \
src/map/wmts.h \
src/map/wmsmap.h \
src/map/wms.h \
src/map/crs.h \
src/map/coordinatesystem.h \
src/data/oziparsers.h src/data/oziparsers.h
SOURCES += src/main.cpp \ SOURCES += src/main.cpp \
src/common/coordinates.cpp \ src/common/coordinates.cpp \
@ -193,18 +194,6 @@ SOURCES += src/main.cpp \
src/map/transform.cpp \ src/map/transform.cpp \
src/map/mapfile.cpp \ src/map/mapfile.cpp \
src/map/tifffile.cpp \ src/map/tifffile.cpp \
src/data/data.cpp \
src/data/poi.cpp \
src/data/track.cpp \
src/data/route.cpp \
src/data/path.cpp \
src/data/gpxparser.cpp \
src/data/tcxparser.cpp \
src/data/csvparser.cpp \
src/data/kmlparser.cpp \
src/data/fitparser.cpp \
src/data/igcparser.cpp \
src/data/nmeaparser.cpp \
src/map/projection.cpp \ src/map/projection.cpp \
src/map/gcs.cpp \ src/map/gcs.cpp \
src/map/angularunits.cpp \ src/map/angularunits.cpp \
@ -218,6 +207,18 @@ SOURCES += src/main.cpp \
src/map/wms.cpp \ src/map/wms.cpp \
src/map/crs.cpp \ src/map/crs.cpp \
src/map/coordinatesystem.cpp \ src/map/coordinatesystem.cpp \
src/data/data.cpp \
src/data/poi.cpp \
src/data/track.cpp \
src/data/route.cpp \
src/data/path.cpp \
src/data/gpxparser.cpp \
src/data/tcxparser.cpp \
src/data/csvparser.cpp \
src/data/kmlparser.cpp \
src/data/fitparser.cpp \
src/data/igcparser.cpp \
src/data/nmeaparser.cpp \
src/data/oziparsers.cpp src/data/oziparsers.cpp
RESOURCES += gpxsee.qrc RESOURCES += gpxsee.qrc
TRANSLATIONS = lang/gpxsee_cs.ts \ TRANSLATIONS = lang/gpxsee_cs.ts \

View File

@ -115,7 +115,7 @@ AlbersEqual::AlbersEqual(const Ellipsoid *ellipsoid, double standardParallel1,
_rho0 = (_C < nq0) ? 0 : _a_over_n * sqrt(_C - nq0); _rho0 = (_C < nq0) ? 0 : _a_over_n * sqrt(_C - nq0);
} }
QPointF AlbersEqual::ll2xy(const Coordinates &c) const PointD AlbersEqual::ll2xy(const Coordinates &c) const
{ {
double dlam; double dlam;
double sin_lat; double sin_lat;
@ -139,11 +139,11 @@ QPointF AlbersEqual::ll2xy(const Coordinates &c) const
rho = (_C < nq) ? 0 : _a_over_n * sqrt(_C - nq); rho = (_C < nq) ? 0 : _a_over_n * sqrt(_C - nq);
theta = _n * dlam; theta = _n * dlam;
return QPointF(rho * sin(theta) + _falseEasting, return PointD(rho * sin(theta) + _falseEasting,
_rho0 - rho * cos(theta) + _falseNorthing); _rho0 - rho * cos(theta) + _falseNorthing);
} }
Coordinates AlbersEqual::xy2ll(const QPointF &p) const Coordinates AlbersEqual::xy2ll(const PointD &p) const
{ {
double dy, dx; double dy, dx;
double rho0_minus_dy; double rho0_minus_dy;

View File

@ -14,8 +14,8 @@ public:
virtual CT *clone() const {return new AlbersEqual(*this);} virtual CT *clone() const {return new AlbersEqual(*this);}
virtual QPointF ll2xy(const Coordinates &c) const; virtual PointD ll2xy(const Coordinates &c) const;
virtual Coordinates xy2ll(const QPointF &p) const; virtual Coordinates xy2ll(const PointD &p) const;
private: private:
double _latitudeOrigin; double _latitudeOrigin;

View File

@ -1,8 +1,8 @@
#ifndef CT_H #ifndef CT_H
#define CT_H #define CT_H
#include <QPointF>
#include "common/coordinates.h" #include "common/coordinates.h"
#include "point.h"
class CT { class CT {
public: public:
@ -10,8 +10,8 @@ public:
virtual CT *clone() const = 0; virtual CT *clone() const = 0;
virtual QPointF ll2xy(const Coordinates &c) const = 0; virtual PointD ll2xy(const Coordinates &c) const = 0;
virtual Coordinates xy2ll(const QPointF &p) const = 0; virtual Coordinates xy2ll(const PointD &p) const = 0;
}; };
#endif // CT_H #endif // CT_H

View File

@ -138,7 +138,7 @@ bool GeoTIFF::readIFD(TIFFFile &file, quint32 &offset, Ctx &ctx) const
return true; return true;
} }
bool GeoTIFF::readScale(TIFFFile &file, quint32 offset, QPointF &scale) const bool GeoTIFF::readScale(TIFFFile &file, quint32 offset, PointD &scale) const
{ {
if (!file.seek(offset)) if (!file.seek(offset))
return false; return false;
@ -155,7 +155,7 @@ bool GeoTIFF::readTiepoints(TIFFFile &file, quint32 offset, quint32 count,
QList<ReferencePoint> &points) const QList<ReferencePoint> &points) const
{ {
double z, pz; double z, pz;
QPointF xy, pp; PointD xy, pp;
if (!file.seek(offset)) if (!file.seek(offset))
return false; return false;
@ -175,10 +175,7 @@ bool GeoTIFF::readTiepoints(TIFFFile &file, quint32 offset, quint32 count,
if (!file.readValue(pz)) if (!file.readValue(pz))
return false; return false;
ReferencePoint p; points.append(ReferencePoint(xy, pp));
p.xy = xy.toPoint();
p.pp = pp;
points.append(p);
} }
return true; return true;
@ -450,7 +447,7 @@ GeoTIFF::GeoTIFF(const QString &path)
{ {
quint32 ifd; quint32 ifd;
QList<ReferencePoint> points; QList<ReferencePoint> points;
QPointF scale; PointD scale;
QMap<quint16, Value> kv; QMap<quint16, Value> kv;
Ctx ctx; Ctx ctx;
TIFFFile file; TIFFFile file;
@ -519,12 +516,12 @@ GeoTIFF::GeoTIFF(const QString &path)
else if (ctx.tiepointCount > 1) else if (ctx.tiepointCount > 1)
_transform = Transform(points); _transform = Transform(points);
else if (ctx.matrix) { else if (ctx.matrix) {
double m[16]; double matrix[16];
if (!readMatrix(file, ctx.matrix, m)) { if (!readMatrix(file, ctx.matrix, matrix)) {
_errorString = "Error reading transformation matrix"; _errorString = "Error reading transformation matrix";
return; return;
} }
_transform = Transform(m); _transform = Transform(matrix);
} else { } else {
_errorString = "Incomplete/missing model transformation info"; _errorString = "Incomplete/missing model transformation info";
return; return;

View File

@ -40,7 +40,7 @@ private:
bool readEntry(TIFFFile &file, Ctx &ctx) const; bool readEntry(TIFFFile &file, Ctx &ctx) const;
bool readIFD(TIFFFile &file, quint32 &offset, Ctx &ctx) const; bool readIFD(TIFFFile &file, quint32 &offset, Ctx &ctx) const;
bool readScale(TIFFFile &file, quint32 offset, QPointF &scale) const; bool readScale(TIFFFile &file, quint32 offset, PointD &scale) const;
bool readTiepoints(TIFFFile &file, quint32 offset, quint32 count, bool readTiepoints(TIFFFile &file, quint32 offset, quint32 count,
QList<ReferencePoint> &points) const; QList<ReferencePoint> &points) const;
bool readMatrix(TIFFFile &file, quint32 offset, double matrix[16]) const; bool readMatrix(TIFFFile &file, quint32 offset, double matrix[16]) const;

View File

@ -31,7 +31,7 @@ LambertAzimuthal::LambertAzimuthal(const Ellipsoid *ellipsoid,
_D = _a * (cos(lat0) / sqrt(1.0 - _es2 * sin2(lat0))) / (_Rq * cos(_beta0)); _D = _a * (cos(lat0) / sqrt(1.0 - _es2 * sin2(lat0))) / (_Rq * cos(_beta0));
} }
QPointF LambertAzimuthal::ll2xy(const Coordinates &c) const PointD LambertAzimuthal::ll2xy(const Coordinates &c) const
{ {
double lon = deg2rad(c.lon()); double lon = deg2rad(c.lon());
double lat = deg2rad(c.lat()); double lat = deg2rad(c.lat());
@ -47,10 +47,10 @@ QPointF LambertAzimuthal::ll2xy(const Coordinates &c) const
double y = _falseNorthing + (B / _D) * ((cos(_beta0) * sin(beta)) double y = _falseNorthing + (B / _D) * ((cos(_beta0) * sin(beta))
- (sin(_beta0) * cos(beta) * cos(lon - _lon0))); - (sin(_beta0) * cos(beta) * cos(lon - _lon0)));
return QPointF(x, y); return PointD(x, y);
} }
Coordinates LambertAzimuthal::xy2ll(const QPointF &p) const Coordinates LambertAzimuthal::xy2ll(const PointD &p) const
{ {
double es4 = _es2 * _es2; double es4 = _es2 * _es2;
double es6 = _es2 * es4; double es6 = _es2 * es4;

View File

@ -13,8 +13,8 @@ public:
virtual CT *clone() const {return new LambertAzimuthal(*this);} virtual CT *clone() const {return new LambertAzimuthal(*this);}
virtual QPointF ll2xy(const Coordinates &c) const; virtual PointD ll2xy(const Coordinates &c) const;
virtual Coordinates xy2ll(const QPointF &p) const; virtual Coordinates xy2ll(const PointD &p) const;
private: private:
double _lon0; double _lon0;

View File

@ -93,7 +93,7 @@ LambertConic1::LambertConic1(const Ellipsoid *ellipsoid, double latitudeOrigin,
_rho_olat = _rho0; _rho_olat = _rho0;
} }
QPointF LambertConic1::ll2xy(const Coordinates &c) const PointD LambertConic1::ll2xy(const Coordinates &c) const
{ {
double t; double t;
double rho; double rho;
@ -117,11 +117,11 @@ QPointF LambertConic1::ll2xy(const Coordinates &c) const
theta = _n * dlam; theta = _n * dlam;
return QPointF(rho * sin(theta) + _falseEasting, _rho_olat - rho return PointD(rho * sin(theta) + _falseEasting, _rho_olat - rho
* cos(theta) + _falseNorthing); * cos(theta) + _falseNorthing);
} }
Coordinates LambertConic1::xy2ll(const QPointF &p) const Coordinates LambertConic1::xy2ll(const PointD &p) const
{ {
double dx; double dx;
double dy; double dy;
@ -266,12 +266,12 @@ LambertConic2::LambertConic2(const Ellipsoid *ellipsoid,
falseEasting, falseNorthing); falseEasting, falseNorthing);
} }
QPointF LambertConic2::ll2xy(const Coordinates &c) const PointD LambertConic2::ll2xy(const Coordinates &c) const
{ {
return _lc1.ll2xy(c); return _lc1.ll2xy(c);
} }
Coordinates LambertConic2::xy2ll(const QPointF &p) const Coordinates LambertConic2::xy2ll(const PointD &p) const
{ {
return _lc1.xy2ll(p); return _lc1.xy2ll(p);
} }

View File

@ -15,8 +15,8 @@ public:
virtual CT *clone() const {return new LambertConic1(*this);} virtual CT *clone() const {return new LambertConic1(*this);}
virtual QPointF ll2xy(const Coordinates &c) const; virtual PointD ll2xy(const Coordinates &c) const;
virtual Coordinates xy2ll(const QPointF &p) const; virtual Coordinates xy2ll(const PointD &p) const;
private: private:
double _longitudeOrigin; double _longitudeOrigin;
@ -40,8 +40,8 @@ public:
virtual CT *clone() const {return new LambertConic2(*this);} virtual CT *clone() const {return new LambertConic2(*this);}
virtual QPointF ll2xy(const Coordinates &c) const; virtual PointD ll2xy(const Coordinates &c) const;
virtual Coordinates xy2ll(const QPointF &p) const; virtual Coordinates xy2ll(const PointD &p) const;
private: private:
LambertConic1 _lc1; LambertConic1 _lc1;

View File

@ -11,9 +11,9 @@ public:
virtual CT *clone() const {return new LatLon(*this);} virtual CT *clone() const {return new LatLon(*this);}
virtual QPointF ll2xy(const Coordinates &c) const virtual PointD ll2xy(const Coordinates &c) const
{return QPointF(_au.fromDegrees(c.lon()), _au.fromDegrees(c.lat()));} {return PointD(_au.fromDegrees(c.lon()), _au.fromDegrees(c.lat()));}
virtual Coordinates xy2ll(const QPointF &p) const virtual Coordinates xy2ll(const PointD &p) const
{return Coordinates(_au.toDegrees(p.x()), _au.toDegrees(p.y()));} {return Coordinates(_au.toDegrees(p.x()), _au.toDegrees(p.y()));}
private: private:

View File

@ -3,7 +3,7 @@
#include <cmath> #include <cmath>
#include <QDebug> #include <QDebug>
#include <QPointF> #include "point.h"
class LinearUnits class LinearUnits
{ {
@ -15,11 +15,11 @@ public:
bool isValid() const {return !std::isnan(_f);} bool isValid() const {return !std::isnan(_f);}
double toMeters(double val) const {return val * _f;} double toMeters(double val) const {return val * _f;}
QPointF toMeters(const QPointF &p) const PointD toMeters(const PointD &p) const
{return QPointF(p.x() * _f, p.y() * _f);} {return PointD(p.x() * _f, p.y() * _f);}
double fromMeters(double val) const {return val / _f;} double fromMeters(double val) const {return val / _f;}
QPointF fromMeters(const QPointF &p) const PointD fromMeters(const PointD &p) const
{return QPointF(p.x() / _f, p.y() /_f);} {return PointD(p.x() / _f, p.y() /_f);}
friend bool operator==(const LinearUnits &lu1, const LinearUnits &lu2); friend bool operator==(const LinearUnits &lu1, const LinearUnits &lu2);
#ifndef QT_NO_DEBUG #ifndef QT_NO_DEBUG

View File

@ -84,7 +84,7 @@ int MapFile::parse(QIODevice &device, QList<CalibrationPoint> &points,
if (list.at(16).trimmed() == "S") if (list.at(16).trimmed() == "S")
p.zone = -p.zone; p.zone = -p.zone;
p.rp.xy = QPoint(x, y); p.rp.setXY(PointD(x, y));
if (ll) { if (ll) {
p.ll = Coordinates(lond + lonm/60.0, latd + latm/60.0); p.ll = Coordinates(lond + lonm/60.0, latd + latm/60.0);
if (p.ll.isValid()) if (p.ll.isValid())
@ -92,7 +92,7 @@ int MapFile::parse(QIODevice &device, QList<CalibrationPoint> &points,
else else
return ln; return ln;
} else if (pp) { } else if (pp) {
p.rp.pp = QPointF(ppx, ppy); p.rp.setPP(PointD(ppx, ppy));
points.append(p); points.append(p);
} else } else
return ln; return ln;
@ -238,8 +238,8 @@ bool MapFile::computeTransformation(QList<CalibrationPoint> &points)
QList<ReferencePoint> rp; QList<ReferencePoint> rp;
for (int i = 0; i < points.size(); i++) { for (int i = 0; i < points.size(); i++) {
if (points.at(i).rp.pp.isNull()) if (points.at(i).rp.pp().isNull())
points[i].rp.pp = _projection.ll2xy(points.at(i).ll); points[i].rp.setPP(_projection.ll2xy(points.at(i).ll));
rp.append(points.at(i).rp); rp.append(points.at(i).rp);
} }

View File

@ -3,13 +3,13 @@
#include "common/wgs84.h" #include "common/wgs84.h"
#include "mercator.h" #include "mercator.h"
QPointF Mercator::ll2xy(const Coordinates &c) const PointD Mercator::ll2xy(const Coordinates &c) const
{ {
return QPointF(deg2rad(c.lon()) * WGS84_RADIUS, return PointD(deg2rad(c.lon()) * WGS84_RADIUS,
log(tan(M_PI/4.0 + deg2rad(c.lat())/2.0)) * WGS84_RADIUS); log(tan(M_PI/4.0 + deg2rad(c.lat())/2.0)) * WGS84_RADIUS);
} }
Coordinates Mercator::xy2ll(const QPointF &p) const Coordinates Mercator::xy2ll(const PointD &p) const
{ {
return Coordinates(rad2deg(p.x() / WGS84_RADIUS), return Coordinates(rad2deg(p.x() / WGS84_RADIUS),
rad2deg(2 * atan(exp(p.y() / WGS84_RADIUS)) - M_PI/2)); rad2deg(2 * atan(exp(p.y() / WGS84_RADIUS)) - M_PI/2));

View File

@ -8,8 +8,8 @@ class Mercator : public CT
public: public:
virtual CT *clone() const {return new Mercator(*this);} virtual CT *clone() const {return new Mercator(*this);}
virtual QPointF ll2xy(const Coordinates &c) const; virtual PointD ll2xy(const Coordinates &c) const;
virtual Coordinates xy2ll(const QPointF &p) const; virtual Coordinates xy2ll(const PointD &p) const;
}; };
#endif // MERCATOR_H #endif // MERCATOR_H

View File

@ -325,7 +325,7 @@ void OfflineMap::draw(QPainter *painter, const QRectF &rect)
QPointF OfflineMap::ll2xy(const Coordinates &c) const QPointF OfflineMap::ll2xy(const Coordinates &c) const
{ {
QPointF p(_transform.proj2img(_projection.ll2xy(c))); QPointF p(_transform.proj2img(_projection.ll2xy(c).toPointF()));
return _ozf ? QPointF(p.x() * _scale.x(), p.y() * _scale.y()) : p; return _ozf ? QPointF(p.x() * _scale.x(), p.y() * _scale.y()) : p;
} }
@ -363,8 +363,10 @@ int OfflineMap::zoomFit(const QSize &size, const RectC &rect)
if (!rect.isValid()) if (!rect.isValid())
rescale(0); rescale(0);
else { else {
QPointF tl(_transform.proj2img(_projection.ll2xy(rect.topLeft()))); QPointF tl(_transform.proj2img(_projection.ll2xy(rect.topLeft())
QPointF br(_transform.proj2img(_projection.ll2xy(rect.bottomRight()))); .toPointF()));
QPointF br(_transform.proj2img(_projection.ll2xy(rect.bottomRight())
.toPointF()));
QRect sbr(QRectF(tl, br).toRect().normalized()); QRect sbr(QRectF(tl, br).toRect().normalized());
for (int i = 0; i < _ozf->zooms(); i++) { for (int i = 0; i < _ozf->zooms(); i++) {

View File

@ -42,7 +42,7 @@ public:
QString errorString() const {return _errorString;} QString errorString() const {return _errorString;}
QPointF ll2pp(const Coordinates &c) const QPointF ll2pp(const Coordinates &c) const
{return _projection.ll2xy(c);} {return _projection.ll2xy(c).toPointF();}
QPointF xy2pp(const QPointF &p) const QPointF xy2pp(const QPointF &p) const
{return _transform.img2proj(p);} {return _transform.img2proj(p);}
QPointF pp2xy(const QPointF &p) const QPointF pp2xy(const QPointF &p) const

35
src/map/point.h Normal file
View File

@ -0,0 +1,35 @@
#ifndef POINT_H
#define POINT_H
#include <cmath>
#include <QPointF>
#include <QDebug>
class PointD {
public:
PointD() : _x(NAN), _y(NAN) {}
PointD(double x, double y) : _x(x), _y(y) {}
PointD(const QPointF &p) : _x(p.x()), _y(p.y()) {}
double x() const {return _x;}
double y() const {return _y;}
double &rx() {return _x;}
double &ry() {return _y;}
bool isNull() const {return std::isnan(_x) && std::isnan(_y);}
QPointF toPointF() const {return QPointF(_x, _y);}
private:
double _x, _y;
};
#ifndef QT_NO_DEBUG
inline QDebug operator<<(QDebug dbg, const PointD &p)
{
dbg.nospace() << "PointD(" << p.x() << ", " << p.y() << ")";
return dbg.space();
}
#endif // QT_NO_DEBUG
#endif // POINT_H

View File

@ -103,13 +103,13 @@ Projection &Projection::operator=(const Projection &p)
return *this; return *this;
} }
QPointF Projection::ll2xy(const Coordinates &c) const PointD Projection::ll2xy(const Coordinates &c) const
{ {
return isValid() return isValid()
? _units.fromMeters(_ct->ll2xy(_gcs->fromWGS84(c))) : QPointF(); ? _units.fromMeters(_ct->ll2xy(_gcs->fromWGS84(c))) : PointD();
} }
Coordinates Projection::xy2ll(const QPointF &p) const Coordinates Projection::xy2ll(const PointD &p) const
{ {
return isValid() return isValid()
? _gcs->toWGS84(_ct->xy2ll(_units.toMeters(p))) : Coordinates(); ? _gcs->toWGS84(_ct->xy2ll(_units.toMeters(p))) : Coordinates();

View File

@ -1,9 +1,9 @@
#ifndef PROJECTION_H #ifndef PROJECTION_H
#define PROJECTION_H #define PROJECTION_H
#include <QPointF>
#include <QDebug> #include <QDebug>
#include "common/coordinates.h" #include "common/coordinates.h"
#include "point.h"
#include "linearunits.h" #include "linearunits.h"
#include "coordinatesystem.h" #include "coordinatesystem.h"
@ -83,8 +83,8 @@ public:
bool isValid() const {return !(_gcs == 0 || _ct == 0 || _units.isNull());} bool isValid() const {return !(_gcs == 0 || _ct == 0 || _units.isNull());}
bool isGeographic() const {return _geographic;} bool isGeographic() const {return _geographic;}
QPointF ll2xy(const Coordinates &c) const; PointD ll2xy(const Coordinates &c) const;
Coordinates xy2ll(const QPointF &p) const; Coordinates xy2ll(const PointD &p) const;
const LinearUnits &units() const {return _units;} const LinearUnits &units() const {return _units;}
const CoordinateSystem &coordinateSystem() const {return _cs;} const CoordinateSystem &coordinateSystem() const {return _cs;}

View File

@ -5,21 +5,17 @@
#define NULL_QTRANSFORM 0,0,0,0,0,0,0,0,0 #define NULL_QTRANSFORM 0,0,0,0,0,0,0,0,0
void Transform::simple(const QList<ReferencePoint> &points) void Transform::simple(const ReferencePoint &p1, const ReferencePoint &p2)
{ {
if (points.at(0).xy.x() == points.at(1).xy.x() if (p1.xy().x() == p2.xy().x() || p1.xy().y() == p2.xy().y()) {
|| points.at(0).xy.y() == points.at(1).xy.y()) {
_errorString = "Invalid reference points tuple"; _errorString = "Invalid reference points tuple";
return; return;
} }
qreal sX, sY, dX, dY; double sX = (p1.xy().x() - p2.xy().x()) / (p1.pp().x() - p2.pp().x());
sX = (points.at(0).xy.x() - points.at(1).xy.x()) double sY = (p2.xy().y() - p1.xy().y()) / (p2.pp().y() - p1.pp().y());
/ (points.at(0).pp.x() - points.at(1).pp.x()); double dX = p2.xy().x() - p2.pp().x() * sX;
sY = (points.at(1).xy.y() - points.at(0).xy.y()) double dY = p1.xy().y() - p1.pp().y() * sY;
/ (points.at(1).pp.y() - points.at(0).pp.y());
dX = points.at(1).xy.x() - points.at(1).pp.x() * sX;
dY = points.at(0).xy.y() - points.at(0).pp.y() * sY;
_proj2img = QTransform(sX, 0, 0, sY, dX, dY); _proj2img = QTransform(sX, 0, 0, sY, dX, dY);
_img2proj = _proj2img.inverted(); _img2proj = _proj2img.inverted();
@ -34,11 +30,11 @@ void Transform::affine(const QList<ReferencePoint> &points)
for (int k = 0; k < points.size(); k++) { for (int k = 0; k < points.size(); k++) {
double f[3], t[2]; double f[3], t[2];
f[0] = points.at(k).pp.x(); f[0] = points.at(k).pp().x();
f[1] = points.at(k).pp.y(); f[1] = points.at(k).pp().y();
f[2] = 1.0; f[2] = 1.0;
t[0] = points.at(k).xy.x(); t[0] = points.at(k).xy().x();
t[1] = points.at(k).xy.y(); t[1] = points.at(k).xy().y();
c.m(i,j) += f[i] * t[j]; c.m(i,j) += f[i] * t[j];
} }
} }
@ -49,8 +45,8 @@ void Transform::affine(const QList<ReferencePoint> &points)
for (int qi = 0; qi < points.size(); qi++) { for (int qi = 0; qi < points.size(); qi++) {
double v[3]; double v[3];
v[0] = points.at(qi).pp.x(); v[0] = points.at(qi).pp().x();
v[1] = points.at(qi).pp.y(); v[1] = points.at(qi).pp().y();
v[2] = 1.0; v[2] = 1.0;
for (size_t i = 0; i < Q.h(); i++) for (size_t i = 0; i < Q.h(); i++)
for (size_t j = 0; j < Q.w(); j++) for (size_t j = 0; j < Q.w(); j++)
@ -79,12 +75,17 @@ Transform::Transform(const QList<ReferencePoint> &points)
if (points.count() < 2) if (points.count() < 2)
_errorString = "Insufficient number of reference points"; _errorString = "Insufficient number of reference points";
else if (points.size() == 2) else if (points.size() == 2)
simple(points); simple(points.at(0), points.at(1));
else else
affine(points); affine(points);
} }
Transform::Transform(const ReferencePoint &p, const QPointF &scale) Transform::Transform(const ReferencePoint &p1, const ReferencePoint &p2)
{
simple(p1, p2);
}
Transform::Transform(const ReferencePoint &p, const PointD &scale)
: _proj2img(NULL_QTRANSFORM), _img2proj(NULL_QTRANSFORM) : _proj2img(NULL_QTRANSFORM), _img2proj(NULL_QTRANSFORM)
{ {
if (scale.x() == 0.0 || scale.y() == 0.0) { if (scale.x() == 0.0 || scale.y() == 0.0) {
@ -92,15 +93,16 @@ Transform::Transform(const ReferencePoint &p, const QPointF &scale)
return; return;
} }
_img2proj = QTransform(scale.x(), 0, 0, -scale.y(), p.pp.x() - p.xy.x() _img2proj = QTransform(scale.x(), 0, 0, -scale.y(), p.pp().x() - p.xy().x()
/ scale.x(), p.pp.y() + p.xy.x() / scale.y()); / scale.x(), p.pp().y() + p.xy().x() / scale.y());
_proj2img = _img2proj.inverted(); _proj2img = _img2proj.inverted();
} }
Transform::Transform(double m[16]) Transform::Transform(double matrix[16])
: _proj2img(NULL_QTRANSFORM), _img2proj(NULL_QTRANSFORM) : _proj2img(NULL_QTRANSFORM), _img2proj(NULL_QTRANSFORM)
{ {
_img2proj = QTransform(m[0], m[1], m[4], m[5], m[3], m[7]); _img2proj = QTransform(matrix[0], matrix[1], matrix[4], matrix[5],
matrix[3], matrix[7]);
if (!_img2proj.isInvertible()) if (!_img2proj.isInvertible())
_errorString = "Singular transformation matrix"; _errorString = "Singular transformation matrix";
else else
@ -110,7 +112,7 @@ Transform::Transform(double m[16])
#ifndef QT_NO_DEBUG #ifndef QT_NO_DEBUG
QDebug operator<<(QDebug dbg, const ReferencePoint &p) QDebug operator<<(QDebug dbg, const ReferencePoint &p)
{ {
dbg.nospace() << "ReferencePoint(" << p.xy << ", " << p.pp << ")"; dbg.nospace() << "ReferencePoint(" << p.xy() << ", " << p.pp() << ")";
return dbg.space(); return dbg.space();
} }
#endif // QT_NO_DEBUG #endif // QT_NO_DEBUG

View File

@ -4,10 +4,21 @@
#include <QTransform> #include <QTransform>
#include <QList> #include <QList>
#include <QDebug> #include <QDebug>
#include "point.h"
struct ReferencePoint { class ReferencePoint
QPoint xy; {
QPointF pp; public:
ReferencePoint() {}
ReferencePoint(const PointD &xy, const PointD &pp) : _xy(xy), _pp(pp) {}
const PointD &xy() const {return _xy;}
const PointD &pp() const {return _pp;}
void setXY(const PointD &xy) {_xy = xy;}
void setPP(const PointD &pp) {_pp = pp;}
private:
PointD _xy, _pp;
}; };
class Transform class Transform
@ -15,8 +26,9 @@ class Transform
public: public:
Transform(); Transform();
Transform(const QList<ReferencePoint> &points); Transform(const QList<ReferencePoint> &points);
Transform(const ReferencePoint &p, const QPointF &scale); Transform(const ReferencePoint &p1, const ReferencePoint &p2);
Transform(double m[16]); Transform(const ReferencePoint &p, const PointD &scale);
Transform(double matrix[16]);
QPointF proj2img(const QPointF &p) const {return _proj2img.map(p);} QPointF proj2img(const QPointF &p) const {return _proj2img.map(p);}
QPointF img2proj(const QPointF &p) const {return _img2proj.map(p);} QPointF img2proj(const QPointF &p) const {return _img2proj.map(p);}
@ -26,7 +38,7 @@ public:
const QString &errorString() const {return _errorString;} const QString &errorString() const {return _errorString;}
private: private:
void simple(const QList<ReferencePoint> &points); void simple(const ReferencePoint &p1, const ReferencePoint &p2);
void affine(const QList<ReferencePoint> &points); void affine(const QList<ReferencePoint> &points);
QTransform _proj2img, _img2proj; QTransform _proj2img, _img2proj;

View File

@ -93,7 +93,7 @@ TransverseMercator::TransverseMercator(const Ellipsoid *ellipsoid,
_ep = 315.e0 * _a * (tn4 - tn5) / 512.e0; _ep = 315.e0 * _a * (tn4 - tn5) / 512.e0;
} }
QPointF TransverseMercator::ll2xy(const Coordinates &c) const PointD TransverseMercator::ll2xy(const Coordinates &c) const
{ {
double rl; double rl;
double cl, c2, c3, c5, c7; double cl, c2, c3, c5, c7;
@ -164,10 +164,10 @@ QPointF TransverseMercator::ll2xy(const Coordinates &c) const
x = _falseEasting + dlam * t6 + pow(dlam, 3.e0) * t7 + pow(dlam, 5.e0) x = _falseEasting + dlam * t6 + pow(dlam, 3.e0) * t7 + pow(dlam, 5.e0)
* t8 + pow(dlam, 7.e0) * t9; * t8 + pow(dlam, 7.e0) * t9;
return QPointF(x, y); return PointD(x, y);
} }
Coordinates TransverseMercator::xy2ll(const QPointF &p) const Coordinates TransverseMercator::xy2ll(const PointD &p) const
{ {
double cl; double cl;
double de; double de;

View File

@ -14,8 +14,8 @@ public:
virtual CT *clone() const {return new TransverseMercator(*this);} virtual CT *clone() const {return new TransverseMercator(*this);}
virtual QPointF ll2xy(const Coordinates &c) const; virtual PointD ll2xy(const Coordinates &c) const;
virtual Coordinates xy2ll(const QPointF &p) const; virtual Coordinates xy2ll(const PointD &p) const;
private: private:
double _longitudeOrigin; double _longitudeOrigin;

View File

@ -61,21 +61,14 @@ void WMSMap::computeZooms(const RangeF &scaleDenominator)
void WMSMap::updateTransform() void WMSMap::updateTransform()
{ {
double scaleDenominator = _zooms.at(_zoom); double scaleDenominator = _zooms.at(_zoom);
ReferencePoint tl, br;
double pixelSpan = sd2res(scaleDenominator); double pixelSpan = sd2res(scaleDenominator);
if (_projection.isGeographic()) if (_projection.isGeographic())
pixelSpan /= deg2rad(WGS84_RADIUS); pixelSpan /= deg2rad(WGS84_RADIUS);
tl.xy = QPoint(0, 0); ReferencePoint tl(PointD(0, 0), _boundingBox.topLeft());
tl.pp = _boundingBox.topLeft(); ReferencePoint br(PointD(_boundingBox.width() / pixelSpan,
br.xy = QPoint(_boundingBox.width() / pixelSpan, -_boundingBox.height() -_boundingBox.height() / pixelSpan), _boundingBox.bottomRight());
/ pixelSpan); _transform = Transform(tl, br);
br.pp = _boundingBox.bottomRight();
QList<ReferencePoint> points;
points << tl << br;
_transform = Transform(points);
} }
bool WMSMap::loadWMS() bool WMSMap::loadWMS()
@ -91,8 +84,8 @@ bool WMSMap::loadWMS()
_projection = wms.projection(); _projection = wms.projection();
RectC bb = wms.boundingBox().normalized(); RectC bb = wms.boundingBox().normalized();
_boundingBox = QRectF(_projection.ll2xy(Coordinates(bb.topLeft().lon(), _boundingBox = QRectF(_projection.ll2xy(Coordinates(bb.topLeft().lon(),
bb.bottomRight().lat())), _projection.ll2xy(Coordinates( bb.bottomRight().lat())).toPointF(), _projection.ll2xy(Coordinates(
bb.bottomRight().lon(), bb.topLeft().lat()))); bb.bottomRight().lon(), bb.topLeft().lat())).toPointF());
_tileLoader = TileLoader(tileUrl(wms.version()), tilesDir(), _tileLoader = TileLoader(tileUrl(wms.version()), tilesDir(),
_setup.authorization()); _setup.authorization());
@ -173,8 +166,8 @@ qreal WMSMap::resolution(const QRectF &rect) const
int WMSMap::zoomFit(const QSize &size, const RectC &br) int WMSMap::zoomFit(const QSize &size, const RectC &br)
{ {
if (br.isValid()) { if (br.isValid()) {
QRectF tbr(_projection.ll2xy(br.topLeft()), QRectF tbr(_projection.ll2xy(br.topLeft()).toPointF(),
_projection.ll2xy(br.bottomRight())); _projection.ll2xy(br.bottomRight()).toPointF());
QPointF sc(tbr.width() / size.width(), tbr.height() / size.height()); QPointF sc(tbr.width() / size.width(), tbr.height() / size.height());
double resolution = qMax(qAbs(sc.x()), qAbs(sc.y())); double resolution = qMax(qAbs(sc.x()), qAbs(sc.y()));
if (_projection.isGeographic()) if (_projection.isGeographic())
@ -209,7 +202,7 @@ int WMSMap::zoomOut()
QPointF WMSMap::ll2xy(const Coordinates &c) const QPointF WMSMap::ll2xy(const Coordinates &c) const
{ {
return _transform.proj2img(_projection.ll2xy(c)); return _transform.proj2img(_projection.ll2xy(c).toPointF());
} }
Coordinates WMSMap::xy2ll(const QPointF &p) const Coordinates WMSMap::xy2ll(const QPointF &p) const

View File

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

View File

@ -70,27 +70,21 @@ double WMTSMap::sd2res(double scaleDenominator) const
void WMTSMap::updateTransform() void WMTSMap::updateTransform()
{ {
const WMTS::Zoom &z = _zooms.at(_zoom); const WMTS::Zoom &z = _zooms.at(_zoom);
ReferencePoint tl, br;
QPointF topLeft = (_cs.axisOrder() == CoordinateSystem::YX) PointD topLeft = (_cs.axisOrder() == CoordinateSystem::YX)
? QPointF(z.topLeft().y(), z.topLeft().x()) : z.topLeft(); ? PointD(z.topLeft().y(), z.topLeft().x()) : z.topLeft();
double 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); PointD tileSpan(z.tile().width() * pixelSpan, z.tile().height() * pixelSpan);
QPointF bottomRight(topLeft.x() + tileSpan.x() * z.matrix().width(), PointD bottomRight(topLeft.x() + tileSpan.x() * z.matrix().width(),
topLeft.y() - tileSpan.y() * z.matrix().height()); topLeft.y() - tileSpan.y() * z.matrix().height());
tl.xy = QPoint(0, 0); ReferencePoint tl(PointD(0, 0), topLeft);
tl.pp = topLeft; ReferencePoint br(PointD(z.tile().width() * z.matrix().width(),
br.xy = QPoint(z.tile().width() * z.matrix().width(), z.tile().height() * z.matrix().height()), bottomRight);
z.tile().height() * z.matrix().height()); _transform = Transform(tl, br);
br.pp = bottomRight;
QList<ReferencePoint> points;
points << tl << br;
_transform = Transform(points);
} }
void WMTSMap::load() void WMTSMap::load()
@ -131,8 +125,8 @@ QRectF WMTSMap::bounds() const
int WMTSMap::zoomFit(const QSize &size, const RectC &br) int WMTSMap::zoomFit(const QSize &size, const RectC &br)
{ {
if (br.isValid()) { if (br.isValid()) {
QRectF tbr(_projection.ll2xy(br.topLeft()), QRectF tbr(_projection.ll2xy(br.topLeft()).toPointF(),
_projection.ll2xy(br.bottomRight())); _projection.ll2xy(br.bottomRight()).toPointF());
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())); qreal resolution = qMax(qAbs(sc.x()), qAbs(sc.y()));
if (_projection.isGeographic()) if (_projection.isGeographic())
@ -206,7 +200,7 @@ void WMTSMap::draw(QPainter *painter, const QRectF &rect)
QPointF WMTSMap::ll2xy(const Coordinates &c) const QPointF WMTSMap::ll2xy(const Coordinates &c) const
{ {
return _transform.proj2img(_projection.ll2xy(c)); return _transform.proj2img(_projection.ll2xy(c).toPointF());
} }
Coordinates WMTSMap::xy2ll(const QPointF &p) const Coordinates WMTSMap::xy2ll(const QPointF &p) const