[mkgmap-dev] HGT - getElevation()

Tue Jan 9 20:04:20 GMT 2018

Hi,

Gerd, your additional interpolation would be redundant.
Side note: comparison of double and integer doesn't look safe. Maybe 
would be better if interpolateHeight() returned integer?

Frank, I understand your idea, but please note, that your choice of 
triangle is arbitrary. There are 2 ways of dividing a square and if you 
would choose the other way, results would be different and not only for 
the middle point. This suggest that the solution isn't complete. Classic 
bilinear interpolation seems better, but differences are small and 
actual map looks more or less the same.

I have attached second patch, which apply interpolation when only 2 
values on any edge are present.

Another topic. I observe some DEM artifacts, when displaying area near 
the border of dem-poly. There are small rectangles without shading 
within a bigger tile. They appear at random, when scrolling or zooming 
map in BaseCamp or Mapsource. Ctrl-G doesn't help. See attached img.

I haven't noticed the problem, when compiling a map without dem-poly. 
Maybe it is a result of slanted edge of DEM clipped by poly.

-- 
Best regards,
Andrzej
-------------- next part --------------
Index: src/uk/me/parabola/mkgmap/reader/hgt/HGTConverter.java
===================================================================

--- src/uk/me/parabola/mkgmap/reader/hgt/HGTConverter.java	(revision 4040)
+++ src/uk/me/parabola/mkgmap/reader/hgt/HGTConverter.java	(working copy)
@@ -109,7 +109,7 @@
 		int hLB = rdr.ele(xLeft, yBottom);
 		int hRB = rdr.ele(xRight, yBottom);
 		lastRow = row;
-		double rc = interpolatedHeightInNormatedRectangle(x1-xLeft, y1-yBottom, hLT, hRT, hRB, hLB);
+		double rc = interpolatedHeight(x1-xLeft, y1-yBottom, hLT, hRT, hRB, hLB);
 		if (rc == HGTReader.UNDEF) {
 			int sum = 0;
 			int valid = 0;
@@ -147,10 +147,10 @@
 			}
 		}
 	}
-	
+
 	/**
 	 * Interpolate the height of point p from the 4 closest values in the hgt matrix.
-	 * Code is a copy from Frank Stinners program BuildDEMFile (Hgtreader.cs) 
+	 * Bilinear interpolation with single node restore
 	 * @param qx value from 0 .. 1 gives relative x position in matrix 
 	 * @param qy value from 0 .. 1 gives relative y position in matrix
 	 * @param hlt height left top
@@ -159,42 +159,70 @@
 	 * @param hlb height left bottom
 	 * @return the interpolated height
 	 */
-	double interpolatedHeightInNormatedRectangle(double qx, double qy, int hlt, int hrt, int hrb, int hlb) {
-		if (hlb == HGTReader.UNDEF || hrt == HGTReader.UNDEF)
-			return HGTReader.UNDEF; // keine Berechnung m??glich
-
-		/*
-		 * In welchem Dreieck liegt der Punkt? oben +-/ |/
-		 * 
-		 * unten /| /-+
-		 */
-		if (qy >= qx) { // oberes Dreieck aus hlb, hrt und hlt (Anstieg py/px
-						// ist gr???Ÿer als height/width)
-
-			if (hlt == HGTReader.UNDEF)
+	double interpolatedHeight(double qx, double qy, int hlt, int hrt, int hrb, int hlb) {
+		// extrapolate single node height if requested point is not near
+		// for multiple missing nodes, return the height of the neares node
+		if (hlb == HGTReader.UNDEF) {
+			if (hrb == HGTReader.UNDEF || hlt == HGTReader.UNDEF || hrt == HGTReader.UNDEF) {
+				if (hrt != HGTReader.UNDEF && hlt != HGTReader.UNDEF && qy > 0.5D)	//top edge
+					return (1.0D - qx)*hlt + qx*hrt;
+				if (hrt != HGTReader.UNDEF && hrb != HGTReader.UNDEF && qx > 0.5D)	//right edge
+					return (1.0D - qy)*hrb + qy*hrt;
+				//if (hlt != HGTReader.UNDEF && hrb != HGTReader.UNDEF && qx + qy > 0.5D && gx + qy < 1.5D)	//diagonal
+				// nearest value
+				return (double)((qx < 0.5D)? ((qy < 0.5D)? hlb: hlt): ((qy < 0.5D)? hrb: hrt));
+			}
+			if (qx + qy < 0.4D)	// point is near missing value
 				return HGTReader.UNDEF;
-
-			// hlt als Koordinatenursprung normieren; mit hrt und hlb 3 Punkte
-			// einer Ebene (3-Punkt-Gleichung)
-			hrt -= hlt;
-			hlb -= hlt;
-			qy -= 1;
-
-			return hlt + qx * hrt - qy * hlb;
-
-		} else { // unteres Dreieck aus hlb, hrb und hrt
-
-			if (hrb == HGTReader.UNDEF)
+			hlb = hlt + hrb - hrt;
+		}
+		else if (hrt == HGTReader.UNDEF) {
+			if (hlb == HGTReader.UNDEF || hrb == HGTReader.UNDEF || hlt == HGTReader.UNDEF) {
+				if (hlb != HGTReader.UNDEF && hrb != HGTReader.UNDEF && qy < 0.5D)	//lower edge
+					return (1.0D - qx)*hlb + qx*hrb;
+				if (hlb != HGTReader.UNDEF && hlt != HGTReader.UNDEF && qx < 0.5D)	//left edge
+					return (1.0D - qy)*hlb + qy*hlt;
+				//if (hlt != HGTReader.UNDEF && hrb != HGTReader.UNDEF && qx + qy > 0.5D && gx + qy < 1.5D)	//diagonal
+				// nearest value
+				return (double)((qx < 0.5D)? ((qy < 0.5D)? hlb: hlt): ((qy < 0.5D)? hrb: hrt));
+			}
+			if (qx + qy > 1.6D)	// point is near missing value
 				return HGTReader.UNDEF;
+			hrt = hlt + hrb - hlb;
+		}
+		else if (hrb == HGTReader.UNDEF) {
+			if (hlb == HGTReader.UNDEF || hlt == HGTReader.UNDEF || hrt == HGTReader.UNDEF) {
+				if (hlt != HGTReader.UNDEF && hrt != HGTReader.UNDEF && qy > 0.5D)	//top edge
+					return (1.0D - qx)*hlt + qx*hrt;
+				if (hlt != HGTReader.UNDEF && hlb != HGTReader.UNDEF && qx < 0.5D)	//left edge
+					return (1.0D - qy)*hlb + qy*hlt;
+				//if (hlb != HGTReader.UNDEF && hrt != HGTReader.UNDEF && qy > qx - 0.5D && qy < qx + 0.5D)	//diagonal
+				// nearest value
+				return (double)((qx < 0.5D)? ((qy < 0.5D)? hlb: hlt): ((qy < 0.5D)? hrb: hrt));
+			}
+			if (qy < qx - 0.4D)	// point is near missing value
+				return HGTReader.UNDEF;
+			hrb = hlb + hrt - hlt;
+		}
+		else if (hlt == HGTReader.UNDEF) {
+			if (hlb == HGTReader.UNDEF || hrb == HGTReader.UNDEF || hrt == HGTReader.UNDEF) {
+				if (hrb != HGTReader.UNDEF && hlb != HGTReader.UNDEF && qy < 0.5D)	//lower edge
+					return (1.0D - qx)*hlb + qx*hrb;
+				if (hrb != HGTReader.UNDEF && hrt != HGTReader.UNDEF && qx > 0.5D)	//right edge
+					return (1.0D - qy)*hrb + qy*hrt;
+				//if (hlb != HGTReader.UNDEF && hrt != HGTReader.UNDEF && qy > qx - 0.5D && qy < qx + 0.5D)	//diagonal
+				// nearest value
+				return (double)((qx < 0.5D)? ((qy < 0.5D)? hlb: hlt): ((qy < 0.5D)? hrb: hrt));
+			}
+			if (qy > qx + 0.6D)	// point is near missing value
+				return HGTReader.UNDEF;
+			hlt = hlb + hrt - hrb;
+		}
 
-			// hrb als Koordinatenursprung normieren; mit hrt und hlb 3 Punkte
-			// einer Ebene (3-Punkt-Gleichung)
-			hrt -= hrb;
-			hlb -= hrb;
-			qx -= 1;
-
-			return hrb - qx * hlb + qy * hrt;
-		}
+		// bilinera interpolation
+		double hxt = (1.0D - qx)*hlt + qx*hrt;
+		double hxb = (1.0D - qx)*hlb + qx*hrb;
+		return (1.0D - qy)*hxb + qy*hxt;
 	}
 
 	public void stats() {
-------------- next part --------------
A non-text attachment was scrubbed...
Name: dem-poly.png
Type: image/png
Size: 43740 bytes
Desc: not available
URL: <http://www.mkgmap.org.uk/pipermail/mkgmap-dev/attachments/20180109/a04eda6a/attachment-0001.png>
