From ff4f3eea60038583f63813aff578cb1918240e40 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Martin=20T=C5=AFma?= <tumic@cbox.cz>
Date: Sun, 19 May 2024 16:14:23 +0200
Subject: [PATCH] Use the map-provided DEM data for hillshading on IMG maps

---
 gpxsee.pro                       |   5 +
 src/map/IMG/bitstream.h          |   1 -
 src/map/IMG/demfile.cpp          | 192 ++++++++++++++++++
 src/map/IMG/demfile.h            |  57 ++++++
 src/map/IMG/demtile.h            |  36 ++++
 src/map/IMG/gmapdata.cpp         |  10 +-
 src/map/IMG/imgdata.cpp          |  10 +-
 src/map/IMG/jls.cpp              | 323 +++++++++++++++++++++++++++++++
 src/map/IMG/jls.h                |  86 ++++++++
 src/map/IMG/mapdata.cpp          |  29 ++-
 src/map/IMG/mapdata.h            |  45 ++++-
 src/map/IMG/rastertile.cpp       | 170 ++++++++++++++--
 src/map/IMG/rastertile.h         |  41 +++-
 src/map/IMG/subfile.h            |  19 +-
 src/map/IMG/vectortile.cpp       |  85 +++++++-
 src/map/IMG/vectortile.h         |  25 ++-
 src/map/hillshading.cpp          |  20 +-
 src/map/hillshading.h            |   2 +-
 src/map/mapsforge/rastertile.cpp |   6 +-
 src/map/mapsforge/rastertile.h   |   2 +-
 src/map/matrix.cpp               |  57 ++----
 src/map/matrix.h                 |  43 ++--
 src/map/transform.cpp            |  14 +-
 23 files changed, 1145 insertions(+), 133 deletions(-)
 create mode 100644 src/map/IMG/demfile.cpp
 create mode 100644 src/map/IMG/demfile.h
 create mode 100644 src/map/IMG/demtile.h
 create mode 100644 src/map/IMG/jls.cpp
 create mode 100644 src/map/IMG/jls.h

diff --git a/gpxsee.pro b/gpxsee.pro
index 15a1bec0..ea0ff7b6 100644
--- a/gpxsee.pro
+++ b/gpxsee.pro
@@ -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 \
diff --git a/src/map/IMG/bitstream.h b/src/map/IMG/bitstream.h
index 93c74e52..43529e96 100644
--- a/src/map/IMG/bitstream.h
+++ b/src/map/IMG/bitstream.h
@@ -76,7 +76,6 @@ private:
 	bool readBytesAligned(quint32 bytes, quint32 &val);
 };
 
-
 template<typename T>
 bool BitStream1::read(quint32 bits, T &val)
 {
diff --git a/src/map/IMG/demfile.cpp b/src/map/IMG/demfile.cpp
new file mode 100644
index 00000000..d172cb6c
--- /dev/null
+++ b/src/map/IMG/demfile.cpp
@@ -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;
+}
diff --git a/src/map/IMG/demfile.h b/src/map/IMG/demfile.h
new file mode 100644
index 00000000..aae83d71
--- /dev/null
+++ b/src/map/IMG/demfile.h
@@ -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
diff --git a/src/map/IMG/demtile.h b/src/map/IMG/demtile.h
new file mode 100644
index 00000000..7568f17d
--- /dev/null
+++ b/src/map/IMG/demtile.h
@@ -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
diff --git a/src/map/IMG/gmapdata.cpp b/src/map/IMG/gmapdata.cpp
index 712eb65c..58b141f0 100644
--- a/src/map/IMG/gmapdata.cpp
+++ b/src/map/IMG/gmapdata.cpp
@@ -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;
 }
diff --git a/src/map/IMG/imgdata.cpp b/src/map/IMG/imgdata.cpp
index 81672b3a..a180b570 100644
--- a/src/map/IMG/imgdata.cpp
+++ b/src/map/IMG/imgdata.cpp
@@ -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);
diff --git a/src/map/IMG/jls.cpp b/src/map/IMG/jls.cpp
new file mode 100644
index 00000000..74fa3e95
--- /dev/null
+++ b/src/map/IMG/jls.cpp
@@ -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;
+}
diff --git a/src/map/IMG/jls.h b/src/map/IMG/jls.h
new file mode 100644
index 00000000..c501df4b
--- /dev/null
+++ b/src/map/IMG/jls.h
@@ -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
diff --git a/src/map/IMG/mapdata.cpp b/src/map/IMG/mapdata.cpp
index e4dfe84c..6fb797d5 100644
--- a/src/map/IMG/mapdata.cpp
+++ b/src/map/IMG/mapdata.cpp
@@ -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);
diff --git a/src/map/IMG/mapdata.h b/src/map/IMG/mapdata.h
index b850fe34..14dbb8db 100644
--- a/src/map/IMG/mapdata.h
+++ b/src/map/IMG/mapdata.h
@@ -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;
 };
diff --git a/src/map/IMG/rastertile.cpp b/src/map/IMG/rastertile.cpp
index fd5c9f03..98925f31 100644
--- a/src/map/IMG/rastertile.cpp
+++ b/src/map/IMG/rastertile.cpp
@@ -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);
 
diff --git a/src/map/IMG/rastertile.h b/src/map/IMG/rastertile.h
index d6afe09d..97687d63 100644
--- a/src/map/IMG/rastertile.h
+++ b/src/map/IMG/rastertile.h
@@ -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;
diff --git a/src/map/IMG/subfile.h b/src/map/IMG/subfile.h
index bdfcdf9d..d0eeaee8 100644
--- a/src/map/IMG/subfile.h
+++ b/src/map/IMG/subfile.h
@@ -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;
diff --git a/src/map/IMG/vectortile.cpp b/src/map/IMG/vectortile.cpp
index 814b6a6a..186fe374 100644
--- a/src/map/IMG/vectortile.cpp
+++ b/src/map/IMG/vectortile.cpp
@@ -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)
 {
diff --git a/src/map/IMG/vectortile.h b/src/map/IMG/vectortile.h
index 9f0a4214..ca663460 100644
--- a/src/map/IMG/vectortile.h
+++ b/src/map/IMG/vectortile.h
@@ -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;
 };
 
 }
diff --git a/src/map/hillshading.cpp b/src/map/hillshading.cpp
index 414fe8a4..f50cd43d 100644
--- a/src/map/hillshading.cpp
+++ b/src/map/hillshading.cpp
@@ -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);
diff --git a/src/map/hillshading.h b/src/map/hillshading.h
index affc4e55..3ab08beb 100644
--- a/src/map/hillshading.h
+++ b/src/map/hillshading.h
@@ -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);
 };
 
diff --git a/src/map/mapsforge/rastertile.cpp b/src/map/mapsforge/rastertile.cpp
index 0ee12397..9ba9cee3 100644
--- a/src/map/mapsforge/rastertile.cpp
+++ b/src/map/mapsforge/rastertile.cpp
@@ -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;
diff --git a/src/map/mapsforge/rastertile.h b/src/map/mapsforge/rastertile.h
index e0a33028..47613931 100644
--- a/src/map/mapsforge/rastertile.h
+++ b/src/map/mapsforge/rastertile.h
@@ -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;
diff --git a/src/map/matrix.cpp b/src/map/matrix.cpp
index fd20c52b..50c94a79 100644
--- a/src/map/matrix.cpp
+++ b/src/map/matrix.cpp
@@ -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
diff --git a/src/map/matrix.h b/src/map/matrix.h
index 8e8c0c3b..942bdb40 100644
--- a/src/map/matrix.h
+++ b/src/map/matrix.h
@@ -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
diff --git a/src/map/transform.cpp b/src/map/transform.cpp
index d91acb34..a57ac974 100644
--- a/src/map/transform.cpp
+++ b/src/map/transform.cpp
@@ -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();
 }