1
0
mirror of https://github.com/tumic0/GPXSee.git synced 2024-10-06 06:43:22 +02:00

Use the map-provided DEM data for hillshading on IMG maps

This commit is contained in:
Martin Tůma 2024-05-19 16:14:23 +02:00
parent d46ac8435e
commit ff4f3eea60
23 changed files with 1145 additions and 133 deletions

View File

@ -128,6 +128,9 @@ HEADERS += src/common/config.h \
src/map/ENC/objects.h \
src/map/ENC/rastertile.h \
src/map/ENC/style.h \
src/map/IMG/demfile.h \
src/map/IMG/demtile.h \
src/map/IMG/jls.h \
src/map/IMG/section.h \
src/map/IMG/zoom.h \
src/map/conversion.h \
@ -345,6 +348,8 @@ SOURCES += src/main.cpp \
src/map/ENC/mapdata.cpp \
src/map/ENC/rastertile.cpp \
src/map/ENC/style.cpp \
src/map/IMG/demfile.cpp \
src/map/IMG/jls.cpp \
src/map/conversion.cpp \
src/map/encatlas.cpp \
src/map/encmap.cpp \

View File

@ -76,7 +76,6 @@ private:
bool readBytesAligned(quint32 bytes, quint32 &val);
};
template<typename T>
bool BitStream1::read(quint32 bits, T &val)
{

192
src/map/IMG/demfile.cpp Normal file
View File

@ -0,0 +1,192 @@
#include "common/garmin.h"
#include "jls.h"
#include "demfile.h"
using namespace IMG;
using namespace Garmin;
static qint16 limit(const DEMTile *tile, quint16 factor)
{
quint8 f1 = (tile->flags() & 1) != 0;
qint16 l = f1 ? tile->diff() - factor : tile->diff() + 1;
if (tile->flags() > 1) {
for (int i = 1; i < 8; i++) {
if (((tile->flags() >> i) & 1) != 0)
l = (l - 1) - (factor << f1);
}
}
return l;
}
void DEMFile::clear()
{
_levels.clear();
}
bool DEMFile::load(Handle &hdl)
{
quint32 u32, zoomData;
quint16 zooms, zoomDataSize;
if (!(seek(hdl, _gmpOffset + 0x15) && readUInt32(hdl, _flags)
&& readUInt16(hdl, zooms) && readUInt32(hdl, u32)
&& readUInt16(hdl, zoomDataSize) && readUInt32(hdl, zoomData)))
return false;
_levels.reserve(zooms);
for (quint16 i = 0; i < zooms; i++) {
quint32 pixelWidth, pixelHeight, pixelWidth2, pixelHeight2, table, cols,
rows, xr, yr, data;
qint32 lon, lat;
quint16 encoding, size, factor;
qint16 minHeight, maxHeight;
quint8 layer, level;
QList<DEMTile> tiles;
if (!(seek(hdl, zoomData + i * zoomDataSize) && readUInt8(hdl, layer)
&& readUInt8(hdl, level) && readUInt32(hdl, pixelHeight)
&& readUInt32(hdl, pixelWidth) && readUInt32(hdl, pixelHeight2)
&& readUInt32(hdl, pixelWidth2) && readUInt16(hdl, factor)
&& readUInt32(hdl, cols) && readUInt32(hdl, rows)
&& readUInt16(hdl, encoding) && readUInt16(hdl, size)
&& readUInt32(hdl, table) && readUInt32(hdl, data)
&& readInt32(hdl, lon) && readInt32(hdl, lat)
&& readUInt32(hdl, yr) && readUInt32(hdl, xr)
&& readInt16(hdl, minHeight) && readInt16(hdl, maxHeight)))
return false;
if (layer)
continue;
if (!seek(hdl, table))
return false;
tiles.reserve((rows + 1) * (cols + 1));
for (quint32 i = 0; i < rows + 1; i++) {
for (quint32 j = 0; j < cols + 1; j++) {
qint32 x = lon + j * pixelWidth * xr;
qint32 y = lat - i * pixelHeight * yr;
quint32 w = (j == cols) ? (pixelWidth2 + 1) : pixelWidth;
quint32 h = (i == rows) ? (pixelHeight2 + 1) : pixelHeight;
RectC r(Coordinates(toWGS32(x), toWGS32(y)),
Coordinates(toWGS32(x + w * xr), toWGS32(y - h * yr)));
quint32 offset;
qint16 base;
quint16 diff;
quint8 flags = 0;
if (!readVUInt32(hdl, (encoding & 0x3) + 1, offset))
return false;
if (encoding & 0x4) {
if (!readInt16(hdl, base))
return false;
} else {
if (!readInt8(hdl, base))
return false;
}
if (encoding & 0x8) {
if (!readUInt16(hdl, diff))
return false;
} else {
if (!readUInt8(hdl, diff))
return false;
}
if ((encoding & 0x10) && !readUInt8(hdl, flags))
return false;
tiles.append(DEMTile(r, w, h, offset, base, diff, flags));
}
}
_levels.append(Level(RectC(tiles.first().rect().topLeft(),
tiles.last().rect().bottomRight()), toWGS32(xr), toWGS32(yr),
toWGS32(pixelWidth * xr), toWGS32(pixelHeight* yr),
data, rows + 1, cols + 1, factor, level, minHeight, maxHeight, tiles));
}
return !_levels.isEmpty();
}
QList<const DEMTile*> DEMFile::tiles(const RectC &rect, int level) const
{
const Level &lvl = _levels.at(level);
QList<const DEMTile*> ret;
RectC ir(lvl.rect & rect);
double left = (ir.left() - lvl.rect.left()) / lvl.txr;
double top = (lvl.rect.top() - ir.top()) / lvl.tyr;
double right = (ir.right() - lvl.rect.left()) / lvl.txr;
double bottom = (lvl.rect.top() - ir.bottom()) / lvl.tyr;
quint32 t = qMin((quint32)top, lvl.rows - 1);
quint32 l = qMin((quint32)left, lvl.cols - 1);
quint32 b = qMin((quint32)bottom, lvl.rows - 1);
quint32 r = qMin((quint32)right, lvl.cols - 1);
ret.reserve((b - t + 1) * (r - l + 1));
for (quint32 i = t; i <= b; i++)
for (quint32 j = l; j <= r; j++)
ret.append(&lvl.tiles.at(lvl.cols * i + j));
return ret;
}
int DEMFile::level(const Zoom &zoom) const
{
/*
for (int i = 0; i < _levels.size(); i++)
if (_levels.at(i).level >= zoom.level())
return i;
return _levels.size() - 1;
*/
Q_UNUSED(zoom);
return 0;
}
MapData::Elevation *DEMFile::elevations(Handle &hdl, int level,
const DEMTile *tile)
{
const Level &l = _levels.at(level);
MapData::Elevation *ele = new MapData::Elevation();
ele->rect = tile->rect();
ele->xr = l.xr;
ele->yr = l.yr;
if (!tile->diff()) {
ele->m = Matrix<qint16>(tile->h(), tile->w(),
tile->flags() ? -32768 : meters(tile->base()));
return ele;
}
if (!seek(hdl, tile->offset() + l.data))
return ele;
quint16 lim = limit(tile, l.factor);
Matrix<qint16> m(tile->h(), tile->w());
JLS jls(tile->diff(), l.factor, tile->w());
if (jls.decode(this, hdl, m)) {
for (int i = 0; i < m.size(); i++) {
if (m.at(i) >= lim)
m.at(i) = -32768;
else {
m.at(i) += tile->base();
if (m.at(i) < l.minHeight)
m.at(i) = l.minHeight;
if (m.at(i) > l.maxHeight)
m.at(i) = l.maxHeight;
m.at(i) = meters(m.at(i));
}
}
ele->m = m;
}
return ele;
}

57
src/map/IMG/demfile.h Normal file
View File

@ -0,0 +1,57 @@
#ifndef DEMFILE_H
#define DEMFILE_H
#include "common/rtree.h"
#include "subfile.h"
#include "demtile.h"
namespace IMG {
class DEMFile : public SubFile
{
public:
DEMFile(const IMGData *img) : SubFile(img) {}
DEMFile(const QString *path) : SubFile(path) {}
DEMFile(const SubFile *gmp, quint32 offset) : SubFile(gmp, offset) {}
bool load(Handle &hdl);
void clear();
MapData::Elevation *elevations(Handle &hdl, int level, const DEMTile *tile);
int level(const Zoom &zoom) const;
QList<const DEMTile *> tiles(const RectC &rect, int level) const;
private:
struct Level {
Level(const RectC &rect, double xr, double yr, double txr, double tyr,
quint32 data, quint32 rows, quint32 cols, quint16 factor,
quint8 level, qint16 minHeight, qint16 maxHeight,
const QList<DEMTile> &tiles)
: rect(rect), xr(xr), yr(yr), txr(txr), tyr(tyr), data(data),
rows(rows), cols(cols), factor(factor), level(level),
minHeight(minHeight), maxHeight(maxHeight), tiles(tiles) {}
RectC rect;
double xr, yr;
double txr, tyr;
quint32 data;
quint32 rows, cols;
quint16 factor;
quint8 level;
qint16 minHeight;
qint16 maxHeight;
QList<DEMTile> tiles;
};
qint16 meters(qint16 val)
{
return (_flags & 1) ? (qint16)qRound(val * 0.3048) : val;
}
quint32 _flags;
QVector<Level> _levels;
};
}
#endif // DEMFILE_H

36
src/map/IMG/demtile.h Normal file
View File

@ -0,0 +1,36 @@
#ifndef DEMTILE_H
#define DEMTILE_H
#include "common/rectc.h"
namespace IMG {
class DEMTile {
public:
DEMTile(const RectC &rect, quint32 w, quint32 h, quint32 offset,
qint16 base, quint16 diff, quint8 flags)
: _rect(rect), _w(w), _h(h), _offset(offset), _base(base),
_diff(diff), _flags(flags) {}
const RectC &rect() const {return _rect;}
quint32 w() const {return _w;}
quint32 h() const {return _h;}
quint32 offset() const {return _offset;}
qint16 base() const {return _base;}
quint16 diff() const {return _diff;}
quint8 flags() const {return _flags;}
private:
RectC _rect;
quint32 _w, _h;
quint32 _offset;
qint16 _base;
quint16 _diff;
quint8 _flags;
};
}
#endif // DEMTILE_H

View File

@ -7,6 +7,9 @@ using namespace IMG;
static SubFile::Type tileType(const QString &suffix)
{
/* Note: we do not load NOD files from non-NT maps as we have no usage
for them */
if (!suffix.compare("TRE"))
return SubFile::TRE;
else if (!suffix.compare("RGN"))
@ -15,10 +18,12 @@ static SubFile::Type tileType(const QString &suffix)
return SubFile::LBL;
else if (!suffix.compare("TYP"))
return SubFile::TYP;
else if (!suffix.compare("GMP"))
return SubFile::GMP;
else if (!suffix.compare("NET"))
return SubFile::NET;
else if (!suffix.compare("DEM"))
return SubFile::DEM;
else if (!suffix.compare("GMP"))
return SubFile::GMP;
else
return SubFile::Unknown;
}
@ -101,6 +106,7 @@ bool GMAPData::loadTile(const QDir &dir)
_tileTree.Insert(min, max, tile);
_bounds |= tile->bounds();
_hasDEM |= tile->hasDem();
return true;
}

View File

@ -8,6 +8,9 @@ using namespace IMG;
static SubFile::Type tileType(const char str[3])
{
/* Note: we do not load NOD files from non-NT maps as we have no usage
for them */
if (!memcmp(str, "TRE", 3))
return SubFile::TRE;
else if (!memcmp(str, "RGN", 3))
@ -16,10 +19,12 @@ static SubFile::Type tileType(const char str[3])
return SubFile::LBL;
else if (!memcmp(str, "TYP", 3))
return SubFile::TYP;
else if (!memcmp(str, "GMP", 3))
return SubFile::GMP;
else if (!memcmp(str, "NET", 3))
return SubFile::NET;
else if (!memcmp(str, "DEM", 3))
return SubFile::DEM;
else if (!memcmp(str, "GMP", 3))
return SubFile::GMP;
else
return SubFile::Unknown;
}
@ -156,6 +161,7 @@ bool IMGData::createTileTree(const TileMap &tileMap)
_tileTree.Insert(min, max, tile);
_bounds |= tile->bounds();
_hasDEM |= tile->hasDem();
}
return (_tileTree.Count() > 0);

323
src/map/IMG/jls.cpp Normal file
View File

@ -0,0 +1,323 @@
#include <cmath>
#include "jls.h"
using namespace IMG;
#define max(a, b) ((a) > (b) ? (a) : (b))
static const quint8 Z[] = {
8, 7, 6, 6, 5, 5, 5, 5, 4, 4, 4, 4, 4, 4, 4, 4,
3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1
};
static const int J[] = {
0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3,
4, 4, 5, 5, 6, 6, 7, 7, 8, 9, 10, 11, 12, 13, 14, 15
};
JLS::JLS(quint16 diff, quint16 factor, quint16 cols)
: _data((cols + 3) * 2)
{
_w = cols;
_maxval = diff;
_near = factor;
_range = ((_maxval + _near * 2) / (_near * 2 + 1)) + 1;
_qbpp = ceil(log2(_range));
quint8 bpp = max(2, ceil(log2(_maxval + 1)));
quint8 LIMIT = 2 * (bpp + max(8, bpp));
_limit = LIMIT - _qbpp - 1;
_runIndex = 0;
_rk = 0;
_rg = 1;
_lrk = 0;
quint16 A = max(2, (_range + 32) / 64);
for (int i = 0; i < 4; i++) {
_a[i] = A;
_b[i] = 0;
_n[i] = 1;
}
_last = _data.data();
_current = _data.data() + (cols + 3);
}
bool JLS::processRunMode(BitStream &bs, quint16 col, quint16 &samples)
{
quint8 z;
quint16 cnt = 0;
while (true) {
if ((qint32)bs.value() < 0) {
z = Z[(bs.value() >> 0x18) ^ 0xff];
for (quint8 i = 0; i < z; i++) {
cnt = cnt + _rg;
if (cnt <= col && _runIndex < 31) {
_runIndex++;
_rk = J[_runIndex];
_rg = 1U << _rk;
}
if (cnt >= col) {
if (!bs.read(i + 1))
return 3;
samples = col;
return true;
}
}
} else
z = 0;
if (z != 8) {
if (!bs.read(z + 1))
return false;
if (_rk) {
samples = (bs.value() >> (32 - _rk)) + cnt;
if (!bs.read(_rk))
return false;
} else
samples = cnt;
_lrk = _rk + 1;
if (_runIndex != 0) {
_runIndex--;
_rk = J[_runIndex];
_rg = 1U << _rk;
}
return true;
}
if (!bs.read(8))
return false;
}
}
bool JLS::decodeError(BitStream &bs, quint8 limit, quint8 k, uint &MErrval)
{
quint8 cnt = 0;
MErrval = 0;
while ((int)bs.value() >= 0) {
cnt = Z[bs.value() >> 0x18];
MErrval += cnt;
if (bs.value() >> 0x18 != 0)
break;
if (!bs.read(8))
return false;
cnt = 0;
}
if (!bs.read(cnt + 1))
return false;
if (MErrval < limit) {
if (k != 0) {
MErrval = (bs.value() >> (0x20 - k)) + (MErrval << k);
if (!bs.read(k))
return false;
}
} else {
MErrval = (bs.value() >> (0x20 - _qbpp)) + 1;
if (!bs.read(_qbpp))
return false;
}
return true;
}
bool JLS::readLine(BitStream &bs)
{
quint8 ictx, rctx;
quint8 k;
uint MErrval;
int Errval;
int Rx = _last[1];
int last0 = _last[0];
int last1 = _last[1];
uint col = 1;
*_current = _last[1];
do {
if (abs(last1 - Rx) > _near) {
int Px = Rx + last1 - last0;
if (Px < 0)
Px = 0;
else if (Px > _maxval)
Px = _maxval;
for (k = 0; _n[1] << k < _a[1]; k++)
;
if (!decodeError(bs, _limit, k, MErrval))
return false;
int mes, meh;
if (MErrval & 1) {
meh = (MErrval + 1) >> 1;
mes = -meh;
} else {
meh = MErrval >> 1;
mes = meh;
}
if ((_near == 0) && (k == 0) && (_b[1] * 2 <= -_n[1])) {
meh = mes + 1;
mes = -mes - 1;
if (MErrval & 1)
meh = mes;
} else
mes = mes * (_near * 2 + 1);
Errval = (Rx < last1) ? mes : -mes;
Rx = Px + Errval;
if (Rx < -_near)
Rx += (_near * 2 + 1) * _range;
else if (Rx > _maxval + _near)
Rx -= (_near * 2 + 1) * _range;
if (Rx < 0)
Rx = 0;
if (Rx > _maxval)
Rx = _maxval;
_a[1] = _a[1] + meh;
_b[1] = _b[1] + mes;
if (_n[1] == 0x40) {
_a[1] = _a[1] >> 1;
_b[1] = _b[1] >> 1;
_n[1] = 0x21;
} else {
_n[1] = _n[1] + 1;
}
if (_b[1] <= -_n[1]) {
_b[1] = _b[1] + _n[1];
if (_b[1] <= -_n[1])
_b[1] = 1 - _n[1];
} else if (_b[1] > 0)
_b[1] = ((_b[1] - _n[1]) >> 0xf) & (_b[1] - _n[1]);
last0 = last1;
last1 = _last[col + 1];
} else {
quint16 samples;
if (!processRunMode(bs, _w - col + 1, samples))
return false;
if (samples != 0) {
for (int i = 0; i < samples; i++) {
if (col > _w)
return false;
_current[col] = Rx;
col++;
}
if (col > _w)
break;
last0 = _last[col];
last1 = _last[col + 1];
} else {
last0 = last1;
last1 = _last[col + 1];
}
rctx = (abs(last0 - Rx) <= _near);
quint16 TEMP = _a[rctx + 2];
if (rctx)
TEMP += _n[rctx + 2] >> 1;
ictx = rctx | 2;
for (k = 0; _n[rctx + 2] << k < TEMP; k++)
;
if (!decodeError(bs, _limit - _lrk, k, MErrval))
return false;
quint16 s = ((k == 0) && (rctx || MErrval)) ?
(_b[ictx] * 2 < _n[ictx]) : 0;
Errval = MErrval + rctx + s;
int evh;
if ((Errval & 1) == 0) {
Errval = Errval / 2;
evh = Errval;
} else {
Errval = s - ((Errval + 1) >> 1);
evh = -Errval;
_b[ictx] = _b[ictx] + 1;
}
Errval *= (_near * 2 + 1);
if (!rctx) {
if (Rx == last0)
return false;
if (Rx < last0)
Rx = last0 + Errval;
else
Rx = last0 - Errval;
} else
Rx = Rx + Errval;
if (Rx < -_near)
Rx += (_near * 2 + 1) * _range;
else if (Rx > _maxval + _near)
Rx -= (_near * 2 + 1) * _range;
if (Rx < 0)
Rx = 0;
if (Rx > _maxval)
Rx = _maxval;
_a[ictx] = _a[ictx] + (evh - rctx);
if (_n[ictx] == 0x40) {
_a[ictx] = _a[ictx] >> 1;
_b[ictx] = _b[ictx] >> 1;
_n[ictx] = 0x21;
} else
_n[ictx] = _n[ictx] + 1;
}
_current[col] = Rx;
col = col + 1;
} while (col <= _w);
return true;
}
bool JLS::decode(const SubFile *file, SubFile::Handle &hdl, Matrix<qint16> &img)
{
BitStream bs(file, hdl);
if (!bs.init())
return false;
for (int i = 0; i < img.h(); i++) {
if (!readLine(bs))
return false;
memcpy(&img.at(i, 0), _current + 1, img.w() * sizeof(quint16));
quint16 *tmp = _last;
_last = _current;
_current = tmp;
}
return true;
}

86
src/map/IMG/jls.h Normal file
View File

@ -0,0 +1,86 @@
#ifndef IMG_JLS_H
#define IMG_JLS_H
#include <QVector>
#include "bitstream.h"
#include "map/matrix.h"
namespace IMG {
class JLS
{
public:
JLS(quint16 diff, quint16 factor, quint16 cols);
bool decode(const SubFile *file, SubFile::Handle &hdl, Matrix<qint16> &img);
private:
class BitStream
{
public:
BitStream(const SubFile *file, SubFile::Handle &hdl)
: _file(file), _hdl(hdl) {}
bool init()
{
if (!_file->readVUInt32SW(_hdl, 4, _value))
return false;
_shift = (quint8)-8;
return true;
}
bool read(quint8 bits)
{
quint8 data;
_value <<= bits;
_shift += bits;
while (-1 < (char)_shift) {
if (!_file->readByte(_hdl, &data))
return false;
_value |= (quint32)data << _shift;
_shift -= 8;
}
return true;
}
quint32 value() const {return _value;}
private:
const SubFile *_file;
SubFile::Handle &_hdl;
quint32 _value;
quint8 _shift;
};
bool readLine(BitStream &bs);
bool processRunMode(BitStream &bs, quint16 col, quint16 &samples);
bool decodeError(BitStream &bs, quint8 limit, quint8 k, uint &MErrval);
quint16 _w;
quint16 _maxval;
quint16 _near;
quint16 _range;
quint8 _qbpp;
quint8 _limit;
quint8 _runIndex;
quint8 _rk;
quint16 _rg;
quint16 _n[4];
quint16 _a[4];
qint16 _b[4];
quint8 _lrk;
QVector<quint16> _data;
quint16 *_current;
quint16 *_last;
};
}
#endif // IMG_JLS_H

View File

@ -7,29 +7,37 @@
using namespace IMG;
#define CACHED_SUBDIVS_COUNT 2048 // ~32MB for both caches together
#define CACHED_SUBDIVS_COUNT 2048 // ~32MB for both caches together
#define CACHED_DEMTILES_COUNT 1024 // ~32MB
bool MapData::polyCb(VectorTile *tile, void *context)
{
PolyCTX *ctx = (PolyCTX*)context;
tile->polys(ctx->rect, ctx->zoom, ctx->polygons, ctx->lines,
ctx->polyCache, ctx->lock);
ctx->cache, ctx->lock);
return true;
}
bool MapData::pointCb(VectorTile *tile, void *context)
{
PointCTX *ctx = (PointCTX*)context;
tile->points(ctx->rect, ctx->zoom, ctx->points, ctx->pointCache, ctx->lock);
tile->points(ctx->rect, ctx->zoom, ctx->points, ctx->cache, ctx->lock);
return true;
}
bool MapData::elevationCb(VectorTile *tile, void *context)
{
ElevationCTX *ctx = (ElevationCTX*)context;
tile->elevations(ctx->rect, ctx->zoom, ctx->elevations, ctx->cache, ctx->lock);
return true;
}
MapData::MapData(const QString &fileName)
: _fileName(fileName), _typ(0), _style(0), _valid(false)
: _fileName(fileName), _typ(0), _style(0), _hasDEM(false), _valid(false)
{
_polyCache.setMaxCost(CACHED_SUBDIVS_COUNT);
_pointCache.setMaxCost(CACHED_SUBDIVS_COUNT);
_demCache.setMaxCost(CACHED_DEMTILES_COUNT);
}
MapData::~MapData()
@ -69,6 +77,19 @@ void MapData::points(const RectC &rect, int bits, QList<Point> *points)
_tileTree.Search(min, max, pointCb, &ctx);
}
void MapData::elevations(const RectC &rect, int bits, QList<Elevation> *elevations)
{
ElevationCTX ctx(rect, zoom(bits), elevations, &_demCache, &_demLock);
double min[2], max[2];
min[0] = rect.left();
min[1] = rect.bottom();
max[0] = rect.right();
max[1] = rect.top();
_tileTree.Search(min, max, elevationCb, &ctx);
}
void MapData::load(qreal ratio)
{
Q_ASSERT(!_style);

View File

@ -10,6 +10,7 @@
#include "common/rtree.h"
#include "common/range.h"
#include "common/hash.h"
#include "map/matrix.h"
#include "label.h"
#include "raster.h"
#include "zoom.h"
@ -20,6 +21,7 @@ class Style;
class SubDiv;
class SubFile;
class VectorTile;
class DEMTile;
class MapData
{
@ -55,6 +57,13 @@ public:
{return id < other.id;}
};
struct Elevation {
Matrix<qint16> m;
RectC rect;
double xr;
double yr;
};
MapData(const QString &fileName);
virtual ~MapData();
@ -65,10 +74,13 @@ public:
void polys(const RectC &rect, int bits, QList<Poly> *polygons,
QList<Poly> *lines);
void points(const RectC &rect, int bits, QList<Point> *points);
void elevations(const RectC &rect, int bits, QList<Elevation> *elevations);
void load(qreal ratio);
void clear();
bool hasDEM() const {return _hasDEM;}
const QString &fileName() const {return _fileName;}
bool isValid() const {return _valid;}
@ -87,6 +99,7 @@ protected:
TileTree _tileTree;
QList<Zoom> _zooms;
Range _zoomLevels;
bool _hasDEM;
bool _valid;
QString _errorString;
@ -103,34 +116,48 @@ private:
typedef QCache<const SubDiv*, Polys> PolyCache;
typedef QCache<const SubDiv*, QList<Point> > PointCache;
typedef QCache<const DEMTile*, Elevation> ElevationCache;
struct PolyCTX
{
PolyCTX(const RectC &rect, const Zoom &zoom,
QList<MapData::Poly> *polygons, QList<MapData::Poly> *lines,
PolyCache *polyCache, QMutex *lock)
PolyCache *cache, QMutex *lock)
: rect(rect), zoom(zoom), polygons(polygons), lines(lines),
polyCache(polyCache), lock(lock) {}
cache(cache), lock(lock) {}
const RectC &rect;
const Zoom &zoom;
QList<MapData::Poly> *polygons;
QList<MapData::Poly> *lines;
PolyCache *polyCache;
PolyCache *cache;
QMutex *lock;
};
struct PointCTX
{
PointCTX(const RectC &rect, const Zoom &zoom,
QList<MapData::Point> *points, PointCache *pointCache, QMutex *lock)
: rect(rect), zoom(zoom), points(points), pointCache(pointCache),
lock(lock) {}
QList<MapData::Point> *points, PointCache *cache, QMutex *lock)
: rect(rect), zoom(zoom), points(points), cache(cache), lock(lock) {}
const RectC &rect;
const Zoom &zoom;
QList<MapData::Point> *points;
PointCache *pointCache;
PointCache *cache;
QMutex *lock;
};
struct ElevationCTX
{
ElevationCTX(const RectC &rect, const Zoom &zoom,
QList<Elevation> *elevations, ElevationCache *cache, QMutex *lock)
: rect(rect), zoom(zoom), elevations(elevations), cache(cache),
lock(lock) {}
const RectC &rect;
const Zoom &zoom;
QList<Elevation> *elevations;
ElevationCache *cache;
QMutex *lock;
};
@ -138,10 +165,12 @@ private:
static bool polyCb(VectorTile *tile, void *context);
static bool pointCb(VectorTile *tile, void *context);
static bool elevationCb(VectorTile *tile, void *context);
PolyCache _polyCache;
PointCache _pointCache;
QMutex _lock;
ElevationCache _demCache;
QMutex _lock, _demLock;
friend class VectorTile;
};

View File

@ -27,6 +27,7 @@ using namespace IMG;
#define ROAD 0
#define WATER 1
static const QColor textColor(Qt::black);
static const QColor haloColor(Qt::white);
static const QColor shieldColor(Qt::white);
@ -107,7 +108,22 @@ static bool rectNearPolygon(const QPolygonF &polygon, const QRectF &rect)
|| polygon.containsPoint(rect.bottomRight(), Qt::OddEvenFill)));
}
const QFont *RasterTile::poiFont(Style::FontSize size, int zoom, bool extended)
static double interpolate(double dx, double dy, double p0, double p1, double p2,
double p3)
{
return p0 * (1.0 - dx) * (1.0 - dy) + p1 * dx * (1.0 - dy)
+ p2 * dy * (1.0 - dx) + p3 * dx * dy;
}
static double val(const Matrix<qint16> &m, int row, int col)
{
qint16 v = m.at(row, col);
return (v == -32768) ? NAN : (double)v;
}
const QFont *RasterTile::poiFont(Style::FontSize size, int zoom,
bool extended) const
{
if (zoom > 25)
size = Style::Normal;
@ -122,7 +138,7 @@ const QFont *RasterTile::poiFont(Style::FontSize size, int zoom, bool extended)
}
}
void RasterTile::ll2xy(QList<MapData::Poly> &polys)
void RasterTile::ll2xy(QList<MapData::Poly> &polys) const
{
for (int i = 0; i < polys.size(); i++) {
MapData::Poly &poly = polys[i];
@ -133,7 +149,7 @@ void RasterTile::ll2xy(QList<MapData::Poly> &polys)
}
}
void RasterTile::ll2xy(QList<MapData::Point> &points)
void RasterTile::ll2xy(QList<MapData::Point> &points) const
{
for (int i = 0; i < points.size(); i++) {
QPointF p(ll2xy(points.at(i).coordinates));
@ -142,7 +158,7 @@ void RasterTile::ll2xy(QList<MapData::Point> &points)
}
void RasterTile::drawPolygons(QPainter *painter,
const QList<MapData::Poly> &polygons)
const QList<MapData::Poly> &polygons) const
{
QCache<const LBLFile *, SubFile::Handle> hc(16);
@ -177,7 +193,9 @@ void RasterTile::drawPolygons(QPainter *painter,
//painter->setPen(Qt::blue);
//painter->setBrush(Qt::NoBrush);
//painter->setRenderHint(QPainter::Antialiasing, false);
//painter->drawRect(QRectF(tl, br));
//painter->setRenderHint(QPainter::Antialiasing);
} else {
const Style::Polygon &style = _data->style()->polygon(poly.type);
@ -189,7 +207,8 @@ void RasterTile::drawPolygons(QPainter *painter,
}
}
void RasterTile::drawLines(QPainter *painter, const QList<MapData::Poly> &lines)
void RasterTile::drawLines(QPainter *painter,
const QList<MapData::Poly> &lines) const
{
painter->setBrush(Qt::NoBrush);
@ -218,7 +237,7 @@ void RasterTile::drawLines(QPainter *painter, const QList<MapData::Poly> &lines)
}
void RasterTile::drawTextItems(QPainter *painter,
const QList<TextItem*> &textItems)
const QList<TextItem*> &textItems) const
{
for (int i = 0; i < textItems.size(); i++)
textItems.at(i)->paint(painter);
@ -446,30 +465,147 @@ void RasterTile::fetchData(QList<MapData::Poly> &polygons,
_data->points(pointRectD.toRectC(_proj, 20), _zoom, &points);
}
Matrix RasterTile::elevation() const
bool RasterTile::edgeCb(const MapData::Elevation *e, void *context)
{
Matrix m(_rect.height() + 2, _rect.width() + 2);
EdgeCTX *ctx = (EdgeCTX*)context;
double x = (ctx->c.lon() - e->rect.left()) / e->xr;
double y = (e->rect.top() - ctx->c.lat()) / e->yr;
int row = qMin((int)y, e->m.h() - 1);
int col = qMin((int)x, e->m.w() - 1);
ctx->ele = val(e->m, row, col);
return std::isnan(ctx->ele);
}
double RasterTile::edge(const DEMTRee &tree, const Coordinates &c)
{
double min[2], max[2];
double ele = NAN;
EdgeCTX ctx(c, ele);
min[0] = c.lon();
min[1] = c.lat();
max[0] = c.lon();
max[1] = c.lat();
tree.Search(min, max, edgeCb, &ctx);
return ele;
}
double RasterTile::elevation(const DEMTRee &tree, const MapData::Elevation *e,
const Coordinates &c)
{
double x = (c.lon() - e->rect.left()) / e->xr;
double y = (e->rect.top() - c.lat()) / e->yr;
int row = qMin((int)y, e->m.h() - 1);
int col = qMin((int)x, e->m.w() - 1);
double p0 = val(e->m, row, col);
double p1 = (col == e->m.w() - 1)
? edge(tree, Coordinates(e->rect.left() + (col + 1) * e->xr + e->xr/2,
e->rect.top() - row * e->yr))
: val(e->m, row, col + 1);
double p2 = (row == e->m.h() - 1)
? edge(tree, Coordinates(e->rect.left() + col * e->xr,
e->rect.top() - (row + 1) * e->yr - e->yr/2))
: val(e->m, row + 1, col);
double p3 = ((col == e->m.w() - 1) || (row == e->m.h() - 1))
? edge(tree, Coordinates(e->rect.left() + (col + 1) * e->xr + e->xr/2,
e->rect.top() - (row + 1) * e->yr - e->yr/2))
: val(e->m, row + 1, col + 1);
return interpolate(x - col, y - row, p0, p1, p2, p3);
}
void RasterTile::buildTree(const QList<MapData::Elevation> &tiles, DEMTRee &tree)
{
double min[2], max[2];
for (int i = 0; i < tiles.size(); i++) {
const MapData::Elevation &e = tiles.at(i);
min[0] = e.rect.left();
min[1] = e.rect.bottom();
max[0] = e.rect.right();
max[1] = e.rect.top();
tree.Insert(min, max, &e);
}
}
bool RasterTile::elevationCb(const MapData::Elevation *e, void *context)
{
ElevationCTX *ctx = (ElevationCTX*)context;
ctx->ele = elevation(ctx->tree, e, ctx->c);
return std::isnan(ctx->ele);
}
void RasterTile::searchTree(const DEMTRee &tree, const Coordinates &c,
double &ele)
{
double min[2], max[2];
ElevationCTX ctx(tree, c, ele);
min[0] = c.lon();
min[1] = c.lat();
max[0] = c.lon();
max[1] = c.lat();
ele = NAN;
tree.Search(min, max, elevationCb, &ctx);
}
MatrixD RasterTile::elevation() const
{
MatrixD m(_rect.height() + 2, _rect.width() + 2);
QVector<Coordinates> ll;
int left = _rect.left() - 1;
int right = _rect.right() + 1;
int top = _rect.top() - 1;
int bottom = _rect.bottom() + 1;
QVector<Coordinates> ll;
ll.reserve(m.w() * m.h());
for (int y = top; y <= bottom; y++) {
for (int y = top; y <= bottom; y++)
for (int x = left; x <= right; x++)
ll.append(xy2ll(QPointF(x, y)));
if (_data->hasDEM()) {
RectC rect;
QList<MapData::Elevation> tiles;
DEMTRee tree;
for (int i = 0; i < ll.size(); i++)
rect = rect.united(ll.at(i));
// Extra margin for edge()
rect = rect.united(xy2ll(QPointF(right + 1, bottom + 1)));
_data->elevations(rect, _zoom, &tiles);
buildTree(tiles, tree);
for (int i = 0; i < ll.size(); i++)
searchTree(tree, ll.at(i), m.at(i));
} else {
DEM::lock();
for (int i = 0; i < ll.size(); i++)
m.at(i) = DEM::elevation(ll.at(i));
DEM::unlock();
}
DEM::lock();
for (int i = 0; i < ll.size(); i++)
m.m(i) = DEM::elevation(ll.at(i));
DEM::unlock();
return m;
}
void RasterTile::drawHillShading(QPainter *painter) const
{
if (_hillShading && _zoom >= 18 && _zoom <= 24)
painter->drawImage(_rect.x(), _rect.y(),
HillShading::render(elevation()));
}
void RasterTile::render()
{
QImage img(_rect.width() * _ratio, _rect.height() * _ratio,
@ -501,14 +637,14 @@ void RasterTile::render()
painter.translate(-_rect.x(), -_rect.y());
drawPolygons(&painter, polygons);
if (_hillShading && _zoom >= 18 && _zoom <= 24)
painter.drawImage(_rect.x(), _rect.y(), HillShading::render(elevation()));
drawHillShading(&painter);
drawLines(&painter, lines);
drawTextItems(&painter, textItems);
qDeleteAll(textItems);
//painter.setPen(Qt::red);
//painter.setBrush(Qt::NoBrush);
//painter.setRenderHint(QPainter::Antialiasing, false);
//painter.drawRect(_rect);

View File

@ -30,18 +30,36 @@ public:
void render();
private:
typedef RTree<const MapData::Elevation*, double, 2> DEMTRee;
struct ElevationCTX {
ElevationCTX(const DEMTRee &tree, const Coordinates &c, double &ele)
: tree(tree), c(c), ele(ele) {}
const DEMTRee &tree;
const Coordinates &c;
double &ele;
};
struct EdgeCTX {
EdgeCTX(const Coordinates &c, double &ele) : c(c), ele(ele) {}
const Coordinates &c;
double &ele;
};
void fetchData(QList<MapData::Poly> &polygons, QList<MapData::Poly> &lines,
QList<MapData::Point> &points);
QPointF ll2xy(const Coordinates &c) const
{return _transform.proj2img(_proj.ll2xy(c));}
Coordinates xy2ll(const QPointF &p) const
{return _proj.xy2ll(_transform.img2proj(p));}
void ll2xy(QList<MapData::Poly> &polys);
void ll2xy(QList<MapData::Point> &points);
void ll2xy(QList<MapData::Poly> &polys) const;
void ll2xy(QList<MapData::Point> &points) const;
void drawPolygons(QPainter *painter, const QList<MapData::Poly> &polygons);
void drawLines(QPainter *painter, const QList<MapData::Poly> &lines);
void drawTextItems(QPainter *painter, const QList<TextItem*> &textItems);
void drawPolygons(QPainter *painter, const QList<MapData::Poly> &polygons) const;
void drawLines(QPainter *painter, const QList<MapData::Poly> &lines) const;
void drawTextItems(QPainter *painter, const QList<TextItem*> &textItems) const;
void drawHillShading(QPainter *painter) const;
void processPolygons(const QList<MapData::Poly> &polygons,
QList<TextItem *> &textItems);
@ -55,9 +73,18 @@ private:
QList<TextItem*> &textItems, const QImage (&arrows)[2]);
const QFont *poiFont(Style::FontSize size = Style::Normal,
int zoom = -1, bool extended = false);
int zoom = -1, bool extended = false) const;
Matrix elevation() const;
MatrixD elevation() const;
static double edge(const DEMTRee &tree, const Coordinates &c);
static double elevation(const DEMTRee &tree, const MapData::Elevation *e,
const Coordinates &c);
static void buildTree(const QList<MapData::Elevation> &tiles, DEMTRee &tree);
static void searchTree(const DEMTRee &tree, const Coordinates &c,
double &ele);
static bool elevationCb(const MapData::Elevation *e, void *context);
static bool edgeCb(const MapData::Elevation *e, void *context);
Projection _proj;
Transform _transform;

View File

@ -13,7 +13,7 @@ namespace IMG {
class SubFile
{
public:
enum Type {Unknown, TRE, RGN, LBL, NET, NOD, TYP, GMP};
enum Type {Unknown, TRE, RGN, LBL, NET, NOD, DEM, TYP, GMP};
class Handle
{
@ -76,6 +76,16 @@ public:
return true;
}
template<typename T>
bool readInt8(Handle &handle, T &val) const
{
quint8 b;
if (!readByte(handle, &b))
return false;
val = (qint8)b;
return true;
}
template<typename T>
bool readUInt16(Handle &handle, T &val) const
{
@ -90,7 +100,7 @@ public:
{
if (!readUInt16(handle, (quint16&)val))
return false;
if((quint16)val > 0x7FFF)
if ((quint16)val > 0x7FFF)
val = (val & 0x7FFF) - 0x8000;
return true;
}
@ -125,6 +135,11 @@ public:
return true;
}
bool readInt32(Handle &handle, qint32 &val) const
{
return readUInt32(handle, (quint32&)val);
}
bool readVUInt32SW(Handle &hdl, quint32 bytes, quint32 &val) const
{
quint8 b;

View File

@ -18,7 +18,6 @@ static void copyPoints(const RectC &rect, QList<MapData::Point> *src,
dst->append(src->at(j));
}
SubFile *VectorTile::file(SubFile::Type type)
{
switch (type) {
@ -32,6 +31,8 @@ SubFile *VectorTile::file(SubFile::Type type)
return _net;
case SubFile::NOD:
return _nod;
case SubFile::DEM:
return _dem;
case SubFile::GMP:
return _gmp;
default:
@ -53,11 +54,12 @@ bool VectorTile::init()
bool VectorTile::initGMP()
{
SubFile::Handle hdl(_gmp);
quint32 tre, rgn, lbl, net, nod;
quint32 tre, rgn, lbl, net, nod, dem;
if (!(_gmp->seek(hdl, 0x19) && _gmp->readUInt32(hdl, tre)
&& _gmp->readUInt32(hdl, rgn) && _gmp->readUInt32(hdl, lbl)
&& _gmp->readUInt32(hdl, net) && _gmp->readUInt32(hdl, nod)))
&& _gmp->readUInt32(hdl, net) && _gmp->readUInt32(hdl, nod)
&& _gmp->readUInt32(hdl, dem)))
return false;
_tre = tre ? new TREFile(_gmp, tre) : 0;
@ -65,6 +67,7 @@ bool VectorTile::initGMP()
_lbl = lbl ? new LBLFile(_gmp, lbl) : 0;
_net = net ? new NETFile(_gmp, net) : 0;
_nod = nod ? new NODFile(_gmp, nod) : 0;
_dem = dem ? new DEMFile(_gmp, dem) : 0;
return true;
}
@ -88,6 +91,18 @@ bool VectorTile::load(SubFile::Handle &rgnHdl, SubFile::Handle &lblHdl,
return true;
}
bool VectorTile::loadDem(SubFile::Handle &hdl)
{
_demLoaded = -1;
if (!_dem || !_dem->load(hdl))
return false;
_demLoaded = 1;
return true;
}
void VectorTile::clear()
{
_tre->clear();
@ -96,13 +111,16 @@ void VectorTile::clear()
_lbl->clear();
if (_net)
_net->clear();
if (_dem)
_dem->clear();
_loaded = 0;
_demLoaded = 0;
}
void VectorTile::polys(const RectC &rect, const Zoom &zoom,
QList<MapData::Poly> *polygons, QList<MapData::Poly> *lines,
MapData::PolyCache *polyCache, QMutex *lock)
MapData::PolyCache *cache, QMutex *lock)
{
SubFile::Handle *rgnHdl = 0, *lblHdl = 0, *netHdl = 0, *nodHdl = 0,
*nodHdl2 = 0;
@ -131,7 +149,7 @@ void VectorTile::polys(const RectC &rect, const Zoom &zoom,
for (int i = 0; i < subdivs.size(); i++) {
SubDiv *subdiv = subdivs.at(i);
MapData::Polys *polys = polyCache->object(subdiv);
MapData::Polys *polys = cache->object(subdiv);
if (!polys) {
quint32 shift = _tre->shift(subdiv->bits());
QList<MapData::Poly> p, l;
@ -165,7 +183,7 @@ void VectorTile::polys(const RectC &rect, const Zoom &zoom,
copyPolys(rect, &p, polygons);
copyPolys(rect, &l, lines);
polyCache->insert(subdiv, new MapData::Polys(p, l));
cache->insert(subdiv, new MapData::Polys(p, l));
} else {
copyPolys(rect, &(polys->polygons), polygons);
copyPolys(rect, &(polys->lines), lines);
@ -178,8 +196,7 @@ void VectorTile::polys(const RectC &rect, const Zoom &zoom,
}
void VectorTile::points(const RectC &rect, const Zoom &zoom,
QList<MapData::Point> *points, QCache<const SubDiv *,
QList<MapData::Point> > *pointCache, QMutex *lock)
QList<MapData::Point> *points, MapData::PointCache *cache, QMutex *lock)
{
SubFile::Handle *rgnHdl = 0, *lblHdl = 0;
@ -207,7 +224,7 @@ void VectorTile::points(const RectC &rect, const Zoom &zoom,
for (int i = 0; i < subdivs.size(); i++) {
SubDiv *subdiv = subdivs.at(i);
QList<MapData::Point> *pl = pointCache->object(subdiv);
QList<MapData::Point> *pl = cache->object(subdiv);
if (!pl) {
QList<MapData::Point> p;
@ -226,7 +243,7 @@ void VectorTile::points(const RectC &rect, const Zoom &zoom,
_rgn->extPointObjects(*rgnHdl, subdiv, _lbl, *lblHdl, &p);
copyPoints(rect, &p, points);
pointCache->insert(subdiv, new QList<MapData::Point>(p));
cache->insert(subdiv, new QList<MapData::Point>(p));
} else
copyPoints(rect, pl, points);
}
@ -236,6 +253,54 @@ void VectorTile::points(const RectC &rect, const Zoom &zoom,
delete rgnHdl; delete lblHdl;
}
void VectorTile::elevations(const RectC &rect, const Zoom &zoom,
QList<MapData::Elevation> *elevations, MapData::ElevationCache *cache,
QMutex *lock)
{
SubFile::Handle *hdl = 0;
lock->lock();
if (_demLoaded < 0) {
lock->unlock();
return;
}
if (!_demLoaded) {
hdl = new SubFile::Handle(_dem);
if (!loadDem(*hdl)) {
lock->unlock();
delete hdl;
return;
}
}
int level = _dem->level(zoom);
QList<const DEMTile*> tiles(_dem->tiles(rect, level));
for (int i = 0; i < tiles.size(); i++) {
const DEMTile *tile = tiles.at(i);
MapData::Elevation *el = cache->object(tile);
if (!el) {
if (!hdl)
hdl = new SubFile::Handle(_dem);
el = _dem->elevations(*hdl, level, tile);
if (!el->m.isNull())
elevations->append(*el);
cache->insert(tile, el);
} else {
if (!el->m.isNull())
elevations->append(*el);
}
}
lock->unlock();
delete hdl;
}
#ifndef QT_NO_DEBUG
QDebug operator<<(QDebug dbg, const VectorTile &tile)
{

View File

@ -6,17 +6,19 @@
#include "lblfile.h"
#include "netfile.h"
#include "nodfile.h"
#include "demfile.h"
namespace IMG {
class VectorTile {
public:
VectorTile()
: _tre(0), _rgn(0), _lbl(0), _net(0), _nod(0), _gmp(0), _loaded(0) {}
: _tre(0), _rgn(0), _lbl(0), _net(0), _nod(0), _dem(0), _gmp(0),
_loaded(0), _demLoaded(0) {}
~VectorTile()
{
delete _tre; delete _rgn; delete _lbl; delete _net; delete _nod;
delete _gmp;
delete _dem; delete _gmp;
}
bool init();
@ -24,21 +26,25 @@ public:
const RectC &bounds() const {return _tre->bounds();}
QVector<Zoom> zooms() const {return _tre->zooms();}
bool hasDem() const {return _dem != 0;}
SubFile *file(SubFile::Type type);
void polys(const RectC &rect, const Zoom &zoom,
QList<MapData::Poly> *polygons, QList<MapData::Poly> *lines,
MapData::PolyCache *polyCache, QMutex *lock);
MapData::PolyCache *cache, QMutex *lock);
void points(const RectC &rect, const Zoom &zoom,
QList<MapData::Point> *points, QCache<const SubDiv*,
QList<MapData::Point> > *pointCache, QMutex *lock);
QList<MapData::Point> *points, MapData::PointCache *cache, QMutex *lock);
void elevations(const RectC &rect, const Zoom &zoom,
QList<MapData::Elevation> *elevations, MapData::ElevationCache *cache,
QMutex *lock);
static bool isTileFile(SubFile::Type type)
{
return (type == SubFile::TRE || type == SubFile::LBL
|| type == SubFile::RGN || type == SubFile::NET
|| type == SubFile::NOD || type == SubFile::GMP);
|| type == SubFile::NOD || type == SubFile::DEM
|| type == SubFile::GMP);
}
template<typename T>
@ -60,6 +66,9 @@ public:
case SubFile::NOD:
_nod = new NODFile(container);
return _nod;
case SubFile::DEM:
_dem = new DEMFile(container);
return _dem;
case SubFile::GMP:
_gmp = new SubFile(container);
return _gmp;
@ -72,15 +81,17 @@ private:
bool initGMP();
bool load(SubFile::Handle &rgnHdl, SubFile::Handle &lblHdl,
SubFile::Handle &netHdl, SubFile::Handle &nodHdl);
bool loadDem(SubFile::Handle &demHdl);
TREFile *_tre;
RGNFile *_rgn;
LBLFile *_lbl;
NETFile *_net;
NODFile *_nod;
DEMFile *_dem;
SubFile *_gmp;
int _loaded;
int _loaded, _demLoaded;
};
}

View File

@ -42,24 +42,24 @@ static void getDerivativesHorn(const SubMatrix &sm, double z, Derivatives &d)
d.dzdy = (z * (sm.z1 + 2 * sm.z2 + sm.z3 - sm.z7 - 2 * sm.z8 - sm.z9)) / 8;
}
static void getSubmatrix(int x, int y, const Matrix &m, SubMatrix &sm)
static void getSubmatrix(int x, int y, const MatrixD &m, SubMatrix &sm)
{
int left = x - 1;
int right = x + 1;
int top = y - 1;
int bottom = y + 1;
sm.z1 = m.m(top, left);
sm.z2 = m.m(top, x);
sm.z3 = m.m(top, right);
sm.z4 = m.m(y, left);
sm.z6 = m.m(y, right);
sm.z7 = m.m(bottom, left);
sm.z8 = m.m(bottom, x);
sm.z9 = m.m(bottom, right);
sm.z1 = m.at(top, left);
sm.z2 = m.at(top, x);
sm.z3 = m.at(top, right);
sm.z4 = m.at(y, left);
sm.z6 = m.at(y, right);
sm.z7 = m.at(bottom, left);
sm.z8 = m.at(bottom, x);
sm.z9 = m.at(bottom, right);
}
QImage HillShading::render(const Matrix &m, quint8 alpha, double z,
QImage HillShading::render(const MatrixD &m, quint8 alpha, double z,
double azimuth, double elevation)
{
QImage img(m.w() - 2, m.h() - 2, QImage::Format_ARGB32_Premultiplied);

View File

@ -7,7 +7,7 @@
class HillShading
{
public:
static QImage render(const Matrix &m, quint8 alpha = 96, double z = 0.3,
static QImage render(const MatrixD &m, quint8 alpha = 96, double z = 0.3,
double azimuth = 315, double elevation = 25);
};

View File

@ -471,9 +471,9 @@ void RasterTile::fetchData(QList<MapData::Path> &paths,
_data->points(pointRectD.toRectC(_proj, 20), _zoom, &points);
}
Matrix RasterTile::elevation() const
MatrixD RasterTile::elevation() const
{
Matrix m(_rect.height() + 2, _rect.width() + 2);
MatrixD m(_rect.height() + 2, _rect.width() + 2);
int left = _rect.left() - 1;
int right = _rect.right() + 1;
@ -489,7 +489,7 @@ Matrix RasterTile::elevation() const
DEM::lock();
for (int i = 0; i < ll.size(); i++)
m.m(i) = DEM::elevation(ll.at(i));
m.at(i) = DEM::elevation(ll.at(i));
DEM::unlock();
return m;

View File

@ -215,7 +215,7 @@ private:
void drawPaths(QPainter *painter, const QList<MapData::Path> &paths,
const QList<MapData::Point> &points, QVector<PainterPath> &painterPaths);
Matrix elevation() const;
MatrixD elevation() const;
Projection _proj;
Transform _transform;

View File

@ -1,79 +1,56 @@
#include "matrix.h"
#define abs(x) ((x)<0 ? -(x) : (x))
Matrix::Matrix(int h, int w)
{
_h = h;
_w = w;
_m.resize(_h * _w);
}
bool Matrix::eliminate(double epsilon)
bool MatrixD::eliminate(double epsilon)
{
double temp;
for (int i = 0; i < _h; i++) {
int maxrow = i;
for (int j = i+1; j < _h; j++)
if (abs(m(j, i)) > abs(m(maxrow, i)))
if (abs(at(j, i)) > abs(at(maxrow, i)))
maxrow = j;
for (int j = 0; j < _w; j++) {
temp = m(i, j);
m(i, j) = m(maxrow, j);
m(maxrow, j) = temp;
temp = at(i, j);
at(i, j) = at(maxrow, j);
at(maxrow, j) = temp;
}
if (abs(m(i, i)) <= epsilon)
if (abs(at(i, i)) <= epsilon)
return false;
for (int j = i+1; j<_h; j++) {
temp = m(j, i) / m(i, i);
temp = at(j, i) / at(i, i);
for (int k = i; k < _w; k++)
m(j, k) -= m(i, k) * temp;
at(j, k) -= at(i, k) * temp;
}
}
for (int i = _h-1; i >= 0; i--) {
temp = m(i, i);
temp = at(i, i);
for (int j = 0; j < i; j++)
for (int k = _w-1; k >= i; k--)
m(j, k) -= m(i, k) * m(j, i) / temp;
m(i, i) /= temp;
at(j, k) -= at(i, k) * at(j, i) / temp;
at(i, i) /= temp;
for (int j = _h; j < _w; j++)
m(i, j) /= temp;
at(i, j) /= temp;
}
return true;
}
Matrix Matrix::augemented(const Matrix &M) const
MatrixD MatrixD::augemented(const MatrixD &M) const
{
if (_h != M._h)
return Matrix();
return MatrixD();
Matrix A(_h, _w + M._w);
MatrixD A(_h, _w + M._w);
for (int i = 0; i < _h; i++)
for (int j = 0; j < _w; j++)
A.m(i, j) = m(i, j);
A.at(i, j) = at(i, j);
for (int i = 0; i < _h; i++)
for (int j = _w; j < A._w; j++)
A.m(i, j) = M.m(i, j-_w);
A.at(i, j) = M.at(i, j-_w);
return A;
}
#ifndef QT_NO_DEBUG
QDebug operator<<(QDebug dbg, const Matrix &matrix)
{
dbg.nospace() << "Matrix(" << "\n";
for (int i = 0; i < matrix.h(); i++) {
for (int j = 0; j < matrix.w(); j++)
dbg << "\t" << matrix.m(i, j);
dbg << "\n";
}
dbg << ")";
return dbg.space();
}
#endif // QT_NO_DEBUG

View File

@ -5,31 +5,52 @@
#include <QVector>
#include <QDebug>
template <class T>
class Matrix
{
public:
Matrix() {_h = 0; _w = 0;}
Matrix(int h, int w);
Matrix(int h, int w) : _h(h), _w(w) {_m.resize(_h * _w);}
Matrix(int h, int w, const T &val) : _m(h * w, val), _h(h), _w(w) {}
int h() const {return _h;}
int w() const {return _w;}
double &m(int n) {return _m[n];}
double &m(int i, int j) {return _m[_w * i + j];}
double const &m(int i, int j) const {return _m.at(_w * i + j);}
T &at(int n) {return _m[n];}
T &at(int i, int j) {return _m[_w * i + j];}
T const &at(int i, int j) const {return _m.at(_w * i + j);}
bool isNull() const {return (_h == 0 || _w == 0);}
int size() const {return _m.size();}
protected:
QVector<T> _m;
int _h, _w;
};
class MatrixD : public Matrix<double>
{
public:
MatrixD() : Matrix<double>() {}
MatrixD(int h, int w) : Matrix<double>(h, w) {}
bool eliminate(double epsilon = DBL_EPSILON);
Matrix augemented(const Matrix &M) const;
private:
QVector<double> _m;
int _h;
int _w;
MatrixD augemented(const MatrixD &M) const;
};
#ifndef QT_NO_DEBUG
QDebug operator<<(QDebug dbg, const Matrix &matrix);
template <class T>
inline QDebug operator<<(QDebug dbg, const Matrix<T> &matrix)
{
dbg.nospace() << "Matrix(" << "\n";
for (int i = 0; i < matrix.h(); i++) {
for (int j = 0; j < matrix.w(); j++)
dbg << "\t" << matrix.at(i, j);
dbg << "\n";
}
dbg << ")";
return dbg.space();
}
#endif // QT_NO_DEBUG
#endif // MATRIX_H

View File

@ -22,7 +22,7 @@ void Transform::simple(const ReferencePoint &p1, const ReferencePoint &p2)
void Transform::affine(const QList<ReferencePoint> &points)
{
Matrix c(3, 2);
MatrixD c(3, 2);
for (int i = 0; i < c.h(); i++) {
for (int j = 0; j < c.w(); j++) {
for (int k = 0; k < points.size(); k++) {
@ -33,12 +33,12 @@ void Transform::affine(const QList<ReferencePoint> &points)
f[2] = 1.0;
t[0] = points.at(k).xy().x();
t[1] = points.at(k).xy().y();
c.m(i,j) += f[i] * t[j];
c.at(i,j) += f[i] * t[j];
}
}
}
Matrix Q(3, 3);
MatrixD Q(3, 3);
for (int qi = 0; qi < points.size(); qi++) {
double v[3];
@ -47,17 +47,17 @@ void Transform::affine(const QList<ReferencePoint> &points)
v[2] = 1.0;
for (int i = 0; i < Q.h(); i++)
for (int j = 0; j < Q.w(); j++)
Q.m(i,j) += v[i] * v[j];
Q.at(i,j) += v[i] * v[j];
}
Matrix M(Q.augemented(c));
MatrixD M(Q.augemented(c));
if (!M.eliminate()) {
_errorString = "Singular transformation matrix";
return;
}
_proj2img = QTransform(M.m(0,3), M.m(0,4), M.m(1,3), M.m(1,4), M.m(2,3),
M.m(2,4));
_proj2img = QTransform(M.at(0,3), M.at(0,4), M.at(1,3), M.at(1,4), M.at(2,3),
M.at(2,4));
_img2proj = _proj2img.inverted();
}