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/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 left = _rect.left() - extend;
-	int right = _rect.right() + 1;
+	int right = _rect.right() + extend;
-	int top = _rect.top() - 1;
+	int top = _rect.top() - extend;
-	int bottom = _rect.bottom() + 1;
+	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)
+	if (_hillShading && _zoom >= 18 && _zoom <= 24) {
-		painter->drawImage(_rect.x(), _rect.y(),
+		MatrixD dem(Filter::blur(elevation(BLUR_RADIUS), BLUR_RADIUS));
-		  HillShading::render(elevation()));
+		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 y = extend; y < m.h() - extend; y++) {
-		for (int x = 1; x < m.w() - 1; x++) {
+		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,
+	static QImage render(const MatrixD &m, int extend, quint8 alpha = 96,
-	  double azimuth = 315, double elevation = 25);
+	  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)
+			if (_hillShading) {
-				painter->drawImage(_rect.x(), _rect.y(),
+				MatrixD dem(Filter::blur(elevation(BLUR_RADIUS), BLUR_RADIUS));
-				  HillShading::render(elevation()));
+				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 left = _rect.left() - extend;
-	int right = _rect.right() + 1;
+	int right = _rect.right() + extend;
-	int top = _rect.top() - 1;
+	int top = _rect.top() - extend;
-	int bottom = _rect.bottom() + 1;
+	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);}