Added DEM filtering (bluring)

gpxsee.pro

@@ -139,6 +139,7 @@ HEADERS += src/common/config.h \
     src/map/encjob.h \
     src/map/encmap.h \
     src/map/ENC/iso8211.h \
+    src/map/filter.h \
     src/map/gemfmap.h \
     src/map/gmifile.h \
     src/map/oruxmap.h \
@@ -356,6 +357,7 @@ SOURCES += src/main.cpp \
     src/map/encatlas.cpp \
     src/map/encmap.cpp \
     src/map/ENC/iso8211.cpp \
+    src/map/filter.cpp \
     src/map/gemfmap.cpp \
     src/map/gmifile.cpp \
     src/map/oruxmap.cpp \

src/map/IMG/demfile.cpp

@@ -138,16 +138,11 @@ QList<const DEMTile*> DEMFile::tiles(const RectC &rect, int level) const
 
 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,
@@ -170,7 +165,7 @@ MapData::Elevation *DEMFile::elevations(Handle &hdl, int level,
 
 	quint16 lim = limit(tile, l.factor);
 	Matrix<qint16> m(tile->h(), tile->w());
-	JLS jls(tile->diff(), l.factor, tile->w());
+	JLS jls(tile->diff(), l.factor);
 	if (jls.decode(this, hdl, m)) {
 		for (int i = 0; i < m.size(); i++) {
 			if (m.at(i) >= lim)

src/map/IMG/jls.cpp

@@ -21,10 +21,8 @@ static const int J[] = {
 	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)
+JLS::JLS(quint16 diff, quint16 factor)
 {
-	_w = cols;
 	_maxval = diff;
 	_near = factor;
 
@@ -33,21 +31,6 @@ JLS::JLS(quint16 diff, quint16 factor, quint16 cols)
 	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)
@@ -308,11 +291,28 @@ bool JLS::decode(const SubFile *file, SubFile::Handle &hdl, Matrix<qint16> &img)
 	if (!bs.init())
 		return false;
 
+	_w = img.w();
+	_data = QVector<quint16>((_w + 3) * 2);
+	_last = _data.data();
+	_current = _data.data() + (_w + 3);
+
+	_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;
+	}
+
 	for (int i = 0; i < img.h(); i++) {
 		if (!readLine(bs))
 			return false;
 
-		memcpy(&img.at(i, 0), _current + 1, img.w() * sizeof(quint16));
+		memcpy(&img.at(i, 0), _current + 1, _w * sizeof(quint16));
 
 		quint16 *tmp = _last;
 		_last = _current;

src/map/IMG/jls.h

@@ -10,7 +10,7 @@ namespace IMG {
 class JLS
 {
 public:
-	JLS(quint16 diff, quint16 factor, quint16 cols);
+	JLS(quint16 diff, quint16 factor);
 
 	bool decode(const SubFile *file, SubFile::Handle &hdl, Matrix<qint16> &img);
 
@@ -36,7 +36,7 @@ private:
 			_value <<= bits;
 			_shift += bits;
 
-			while (-1 < (char)_shift) {
+			while (-1 < (qint8)_shift) {
 				if (!_file->readByte(_hdl, &data))
 					return false;
 
@@ -60,10 +60,8 @@ private:
 	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;
@@ -76,6 +74,7 @@ private:
 	qint16 _b[4];
 	quint8 _lrk;
 
+	quint16 _w;
 	QVector<quint16> _data;
 	quint16 *_current;
 	quint16 *_last;

src/map/IMG/rastertile.cpp

@@ -7,6 +7,7 @@
 #include "map/bitmapline.h"
 #include "map/rectd.h"
 #include "map/hillshading.h"
+#include "map/filter.h"
 #include "style.h"
 #include "lblfile.h"
 #include "dem.h"
@@ -28,6 +29,7 @@ using namespace IMG;
 #define ROAD  0
 #define WATER 1
 
+#define BLUR_RADIUS 3
 #define DELTA 0.05 /* DEM3 resolution in degrees */
 
 static const QColor textColor(Qt::black);
@@ -453,15 +455,15 @@ void RasterTile::fetchData(QList<MapData::Poly> &polygons,
 	_data->points(pointRectD.toRectC(_proj, 20), _zoom, &points);
 }
 
-MatrixD RasterTile::elevation() const
+MatrixD RasterTile::elevation(int extend) const
 {
-	MatrixD m(_rect.height() + 2, _rect.width() + 2);
+	MatrixD m(_rect.height() + 2 * extend, _rect.width() + 2 * extend);
 	QVector<Coordinates> ll;
 
-	int left = _rect.left() - 1;
-	int right = _rect.right() + 1;
-	int top = _rect.top() - 1;
-	int bottom = _rect.bottom() + 1;
+	int left = _rect.left() - extend;
+	int right = _rect.right() + extend;
+	int top = _rect.top() - extend;
+	int bottom = _rect.bottom() + extend;
 
 	ll.reserve(m.w() * m.h());
 	for (int y = top; y <= bottom; y++)
@@ -496,9 +498,11 @@ MatrixD RasterTile::elevation() const
 
 void RasterTile::drawHillShading(QPainter *painter) const
 {
-	if (_hillShading && _zoom >= 18 && _zoom <= 24)
-		painter->drawImage(_rect.x(), _rect.y(),
-		  HillShading::render(elevation()));
+	if (_hillShading && _zoom >= 18 && _zoom <= 24) {
+		MatrixD dem(Filter::blur(elevation(BLUR_RADIUS), BLUR_RADIUS));
+		QImage img(HillShading::render(dem, BLUR_RADIUS));
+		painter->drawImage(_rect.x(), _rect.y(), img);
+	}
 }
 
 void RasterTile::render()

src/map/IMG/rastertile.h

@@ -75,7 +75,7 @@ private:
 	const QFont *poiFont(Style::FontSize size = Style::Normal,
 	  int zoom = -1, bool extended = false) const;
 
-	MatrixD elevation() const;
+	MatrixD elevation(int extend) const;
 
 	Projection _proj;
 	Transform _transform;

src/map/filter.cpp

@@ -0,0 +1,114 @@
+#include <cmath>
+#include "filter.h"
+
+static QVector<int> boxesForGauss(double sigma, int n)
+{
+	double wIdeal = sqrt((12 * sigma * sigma / n) + 1);
+	int wl = floor(wIdeal);
+	if (wl % 2 == 0)
+		wl--;
+	int wu = wl + 2;
+
+	double mIdeal = (12 * sigma*sigma - n * wl * wl - 4 * n * wl - 3 * n)
+	  / (-4 * wl - 4);
+	int m = round(mIdeal);
+
+	QVector<int> sizes(n);
+	for (int i = 0; i < n; i++)
+		sizes[i] = i < m ? wl : wu;
+
+	return sizes;
+}
+
+static void boxBlurH4(const MatrixD &src, MatrixD &dst, int r)
+{
+	double iarr = 1.0 / (r + r + 1);
+
+	for (int i = 0; i < src.h(); i++) {
+		int ti = i * src.w(), li = ti, ri = ti + r;
+		double fv = src.at(ti);
+		double lv = src.at(ti + src.w() - 1);
+		double val = (r + 1) * fv;
+
+		for (int j = 0; j < r; j++)
+			val += src.at(ti + j);
+		for (int j = 0; j <= r; j++) {
+			val += src.at(ri++) - fv;
+			dst.at(ti++) = val * iarr;
+		}
+		for (int j = r + 1; j < src.w() - r; j++) {
+			val += src.at(ri++) - src.at(li++);
+			dst.at(ti++) = val * iarr;
+		}
+		for (int j = src.w() - r; j < src.w(); j++) {
+			val += lv - src.at(li++);
+			dst.at(ti++) = val * iarr;
+		}
+	}
+}
+
+static void boxBlurT4(const MatrixD &src, MatrixD &dst, int r)
+{
+	double iarr = 1.0 / (r + r + 1);
+
+	for (int i = 0; i < src.w(); i++) {
+		int ti = i, li = ti, ri = ti + r * src.w();
+		double fv = src.at(ti);
+		double lv = src.at(ti + src.w() * (src.h() - 1));
+		double val = (r + 1) * fv;
+
+		for (int j = 0; j < r; j++)
+			val += src.at(ti + j * src.w());
+		for (int j = 0; j <= r; j++) {
+			val += src.at(ri) - fv;
+			dst.at(ti) = val * iarr;
+			ri += src.w(); ti += src.w();
+		}
+		for (int j = r + 1; j < src.h() - r; j++) {
+			val += src.at(ri) - src.at(li);
+			dst.at(ti) = val * iarr;
+			li += src.w(); ri += src.w(); ti += src.w();
+		}
+		for (int j = src.h() - r; j < src.h(); j++) {
+			val += lv - src.at(li);
+			dst.at(ti) = val * iarr;
+			li += src.w(); ti += src.w();
+		}
+	}
+}
+
+static void boxBlur4(MatrixD &src, MatrixD &dst, int r)
+{
+	for (int i = 0; i < src.size(); i++)
+		dst.at(i) = src.at(i);
+
+	boxBlurH4(dst, src, r);
+	boxBlurT4(src, dst, r);
+}
+
+static void gaussBlur4(MatrixD &src, MatrixD &dst, int r)
+{
+	QVector<int> bxs(boxesForGauss(r, 3));
+
+	boxBlur4(src, dst, (bxs.at(0) - 1) / 2);
+	boxBlur4(dst, src, (bxs.at(1) - 1) / 2);
+	boxBlur4(src, dst, (bxs.at(2) - 1) / 2);
+}
+
+MatrixD Filter::blur(const MatrixD &m, int radius)
+{
+	MatrixD src(m);
+	MatrixD dst(m.h(), m.w());
+
+	for (int i = 0; i < m.size(); i++)
+		if (std::isnan(m.at(i)))
+			src.at(i) = -500;
+
+	gaussBlur4(src, dst, radius);
+
+	for (int i = 0; i < dst.size(); i++)
+		if (std::isnan(m.at(i)))
+			dst.at(i) = NAN;
+
+	return dst;
+}

src/map/filter.h

@@ -0,0 +1,11 @@
+#ifndef FILTER_H
+#define FILTER_H
+
+#include "matrix.h"
+
+namespace Filter
+{
+	MatrixD blur(const MatrixD &m, int radius);
+}
+
+#endif // FILTER_H

src/map/hillshading.cpp

@@ -59,10 +59,11 @@ static void getSubmatrix(int x, int y, const MatrixD &m, SubMatrix &sm)
 	sm.z9 = m.at(bottom, right);
 }
 
-QImage HillShading::render(const MatrixD &m, quint8 alpha, double z,
+QImage HillShading::render(const MatrixD &m, int extend, quint8 alpha, double z,
   double azimuth, double elevation)
 {
-	QImage img(m.w() - 2, m.h() - 2, QImage::Format_ARGB32_Premultiplied);
+	QImage img(m.w() - 2 * extend, m.h() - 2 * extend,
+	  QImage::Format_ARGB32_Premultiplied);
 	uchar *bits = img.bits();
 	int bpl = img.bytesPerLine();
 
@@ -72,8 +73,8 @@ QImage HillShading::render(const MatrixD &m, quint8 alpha, double z,
 
 	getConstants(azimuth, elevation, c);
 
-	for (int y = 1; y < m.h() - 1; y++) {
-		for (int x = 1; x < m.w() - 1; x++) {
+	for (int y = extend; y < m.h() - extend; y++) {
+		for (int x = extend; x < m.w() - extend; x++) {
 			getSubmatrix(x, y, m, sm);
 			getDerivativesHorn(sm, z, d);
 
@@ -89,7 +90,7 @@ QImage HillShading::render(const MatrixD &m, quint8 alpha, double z,
 				pixel = (alpha - val)<<24;
 			}
 
-			*(quint32*)(bits + (y - 1) * bpl + (x - 1) * 4) = pixel;
+			*(quint32*)(bits + (y - extend) * bpl + (x - extend) * 4) = pixel;
 		}
 	}
 

src/map/hillshading.h

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

src/map/mapsforge/rastertile.cpp

@@ -4,6 +4,7 @@
 #include "common/dem.h"
 #include "map/rectd.h"
 #include "map/hillshading.h"
+#include "map/filter.h"
 #include "rastertile.h"
 
 using namespace Mapsforge;
@@ -12,6 +13,8 @@ using namespace Mapsforge;
 #define PATHS_EXTENT 20
 #define SEARCH_EXTENT -0.5
 
+#define BLUR_RADIUS 3
+
 static double LIMIT = cos(deg2rad(170));
 
 static qreal area(const QPainterPath &polygon)
@@ -438,9 +441,11 @@ void RasterTile::drawPaths(QPainter *painter, const QList<MapData::Path> &paths,
 			painter->setBrush(ri->brush());
 			painter->drawEllipse(ll2xy(point->coordinates), radius, radius);
 		} else {
-			if (_hillShading)
-				painter->drawImage(_rect.x(), _rect.y(),
-				  HillShading::render(elevation()));
+			if (_hillShading) {
+				MatrixD dem(Filter::blur(elevation(BLUR_RADIUS), BLUR_RADIUS));
+				QImage img(HillShading::render(dem, BLUR_RADIUS));
+				painter->drawImage(_rect.x(), _rect.y(), img);
+			}
 		}
 	}
 }
@@ -471,14 +476,14 @@ void RasterTile::fetchData(QList<MapData::Path> &paths,
 	_data->points(pointRectD.toRectC(_proj, 20), _zoom, &points);
 }
 
-MatrixD RasterTile::elevation() const
+MatrixD RasterTile::elevation(int extend) const
 {
-	MatrixD m(_rect.height() + 2, _rect.width() + 2);
+	MatrixD m(_rect.height() + 2 * extend, _rect.width() + 2 * extend);
 
-	int left = _rect.left() - 1;
-	int right = _rect.right() + 1;
-	int top = _rect.top() - 1;
-	int bottom = _rect.bottom() + 1;
+	int left = _rect.left() - extend;
+	int right = _rect.right() + extend;
+	int top = _rect.top() - extend;
+	int bottom = _rect.bottom() + extend;
 
 	QVector<Coordinates> ll;
 	ll.reserve(m.w() * m.h());

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);
 
-	MatrixD elevation() const;
+	MatrixD elevation(int extend) const;
 
 	Projection _proj;
 	Transform _transform;

src/map/matrix.h

@@ -15,6 +15,7 @@ public:
 
 	int h() const {return _h;}
 	int w() const {return _w;}
+	const T &at(int n) const {return _m.at(n);}
 	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);}