[mkgmap-dev] HGT - getElevation()

From Andrzej Popowski popej at poczta.onet.pl on Tue Jan 9 20:04:20 GMT 2018

```Hi,

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

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 --------------
===================================================================
@@ -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);
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) {
-			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)
-
+	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 (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
-
-			// 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
-
+			hlb = hlt + hrb - hrt;
+		}
+		else if (hrt == 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
+			hrt = hlt + hrb - hlb;
+		}
+		else if (hrb == 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
+			hrb = hlb + hrt - hlt;
+		}
+		else if (hlt == 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
+			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>
```