[mkgmap-dev] HGT - getElevation()

Tue Jan 9 18:10:15 GMT 2018

Hi Gerd,

I have prepared a patch with bilinear interpolation, which supports 
missing values. It can restore single missing value, but if more values 
are missing it simply returns the nearest value without interpolation.

-- 
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,42 @@
 	 * @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)
+				return (double)((qx < 0.5D)? ((qy < 0.5D)? hlb: hlt): ((qy < 0.5D)? hrb: hrt));
+			if (qx + qy < 0.4D)
 				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)
+				return (double)((qx < 0.5D)? ((qy < 0.5D)? hlb: hlt): ((qy < 0.5D)? hrb: hrt));
+			if (qx + qy > 1.6D)
 				return HGTReader.UNDEF;
+			hrt = hlt + hrb - hlb;
+		}
+		else if (hrb == HGTReader.UNDEF) {
+			if (hlb == HGTReader.UNDEF || hlt == HGTReader.UNDEF || hrt == HGTReader.UNDEF)
+				return (double)((qx < 0.5D)? ((qy < 0.5D)? hlb: hlt): ((qy < 0.5D)? hrb: hrt));
+			if (qy < qx - 0.4D)
+				return HGTReader.UNDEF;
+			hrb = hlb + hrt - hlt;
+		}
+		else if (hlt == HGTReader.UNDEF) {
+			if (hlb == HGTReader.UNDEF || hrb == HGTReader.UNDEF || hrt == HGTReader.UNDEF)
+				return (double)((qx < 0.5D)? ((qy < 0.5D)? hlb: hlt): ((qy < 0.5D)? hrb: hrt));
+			if (qy > qx + 0.6D)
+				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() {
