[mkgmap-dev] HGT - getElevation()

Fri Jan 19 23:50:19 GMT 2018

Hi Gerd,

I have included bicubic interpolation created by Ken Perlin:
https://github.com/mlevans/pretty-heatmaps/blob/master/DensityArrayCreation/src/DensityArrayCreation/Bicubic.java

Bicubic interpolation is used only for layers with good resolution (3 
arc seconds or better). For other cases bilinear interpolation is used.

Patch is attached.

-- 
Best regards,
Andrzej

-------------- next part --------------
Index: src/uk/me/parabola/imgfmt/app/dem/DEMFile.java
===================================================================

--- src/uk/me/parabola/imgfmt/app/dem/DEMFile.java	(revision 4070)
+++ src/uk/me/parabola/imgfmt/app/dem/DEMFile.java	(working copy)
@@ -52,20 +52,57 @@
 	 * @param outsidePolygonHeight the height value that should be used for points outside of the bounding polygon 
 	 */
 	public void calc(Area area, java.awt.geom.Area demPolygonMapUnits, String pathToHGT, List<Integer> pointDistances, short outsidePolygonHeight) {
-		HGTConverter hgtConverter = new HGTConverter(pathToHGT, area, demPolygonMapUnits);
+		// HGT area is extended by 0.01 degree in each direction
+		HGTConverter hgtConverter = new HGTConverter(pathToHGT, area, demPolygonMapUnits, 0.01D);
 		hgtConverter.setOutsidePolygonHeight(outsidePolygonHeight);
 
+		int top = area.getMaxLat() * 256;
+		int bottom = area.getMinLat() * 256;
+		int left = area.getMinLong() * 256;
+		int right = area.getMaxLong() * 256;
+
 		int zoom = 0;
 		int lastDist = pointDistances.get(pointDistances.size()-1); 
 		for (int pointDist : pointDistances) {
-			DEMSection section = new DEMSection(zoom++, area, hgtConverter, pointDist, pointDist == lastDist);
+			int distance = pointDist;
+			if (distance == -1) {
+				int res = (hgtConverter.getHighestRes() > 0) ? hgtConverter.getHighestRes() : 1200;
+				distance = (int) ((1 << 29) / (res * 45));
+			}
+			// last 4 bits of distance should be 0
+			distance = ((distance + 8)/16)*16;
+
+			int xtop = top;
+			int xleft = left;
+
+			// align DEM to distance raster, if distance not bigger than widening of HGT area
+			if (distance < (int)Math.floor((0.01D/45.0D * (1 << 29)))) {
+				if (xtop >= 0) {
+					xtop -= xtop % distance;
+					if (xtop < 0x3FFFFFFF - distance)
+						xtop += distance;
+				}
+				else {
+					xtop -= xtop % distance;
+				}
+
+				if (xleft >= 0) {
+					xleft -= xleft % distance;
+				}
+				else {
+					xleft -= xleft % distance;
+					if (xleft > Integer.MIN_VALUE + distance)
+						xleft -= distance;
+				}
+			}
+
+			// to improve: DEM secition should get area in 32bit units
+			DEMSection section = new DEMSection(zoom++, xtop, xleft, xtop - bottom, right - xleft, hgtConverter, distance, pointDist == lastDist);
 			demHeader.addSection(section);
 		}
 		return;
 	}
 
-	
-	
 	public void write() {
 		ImgFileWriter w = getWriter();
 		if (w instanceof BufferedImgFileWriter) {
Index: src/uk/me/parabola/imgfmt/app/dem/DEMSection.java
===================================================================
--- src/uk/me/parabola/imgfmt/app/dem/DEMSection.java	(revision 4070)
+++ src/uk/me/parabola/imgfmt/app/dem/DEMSection.java	(working copy)
@@ -48,16 +48,12 @@
 	private int maxHeight = Integer.MIN_VALUE;
 	private List<DEMTile> tiles = new ArrayList<>();
 	
-	public DEMSection(int zoomLevel, Area bbox, HGTConverter hgtConverter, int pointDist, boolean lastLevel) {
+	public DEMSection(int zoomLevel, int areaTop, int areaLeft, int areaHeight, int areaWidth, HGTConverter hgtConverter, int pointDist, boolean lastLevel) {
 		this.zoomLevel = zoomLevel;
 		this.lastLevel = lastLevel;
 		
-		if (pointDist == -1) {
-			int res = (hgtConverter.getHighestRes() > 0) ? hgtConverter.getHighestRes() : 1200;
-			pointDist = (int) ((1 << 29) / (res * 45));
-		}
-		this.top = bbox.getMaxLat() * 256;
-		this.left = bbox.getMinLong() * 256;
+		this.top = areaTop;
+		this.left = areaLeft;
 
 		// calculate raster that starts at top left corner
 		// last row and right column have non-standard height / row values 
@@ -64,8 +60,12 @@
 		pointsDistanceLat = pointDist; 
 		pointsDistanceLon = pointDist;
 		
-		int []latInfo = getTileInfo(bbox.getHeight() * 256, pointsDistanceLat);
-		int []lonInfo = getTileInfo(bbox.getWidth() * 256, pointsDistanceLon);
+		//minimal resolution for bicubic interpolation is 3 arc seconds
+		hgtConverter.setBicubic(pointDist <= 9942);
+		//if (pointDist <= 9942) System.out.println("using bicubic interpolatin for point distance " + pointDist);
+
+		int []latInfo = getTileInfo(areaHeight, pointsDistanceLat);
+		int []lonInfo = getTileInfo(areaWidth, pointsDistanceLon);
 		// store the values written to the header
 		tilesLat = latInfo[0];
 		tilesLon = lonInfo[0];
@@ -113,6 +113,7 @@
 		int maxDeltaHeight = Integer.MIN_VALUE;
 		hgtConverter.setLatDist(pointsDistanceLat);
 		hgtConverter.setLonDist(pointsDistanceLon);
+		hgtConverter.clearStat();
 		for (int m = 0; m < tilesLat; m++) {
 			latOff = top - m * resLat;
 			
@@ -150,6 +151,8 @@
 			}
 		}
 		
+		hgtConverter.printStat();
+
 		if (dataLen > 0) {
 			minHeight = minBaseHeight;
 		} else { 
Index: src/uk/me/parabola/mkgmap/reader/hgt/Bicubic.java
===================================================================
--- src/uk/me/parabola/mkgmap/reader/hgt/Bicubic.java	(nonexistent)
+++ src/uk/me/parabola/mkgmap/reader/hgt/Bicubic.java	(working copy)
@@ -0,0 +1,79 @@
+/*
+Java class to implement a bicubic interpolator
+      by Ken Perlin @ NYU, 1998.
+
+You have my permission to use freely, as long as you keep the attribution. - Ken Perlin
+
+Note: this Bicubic.html file also works as a legal Bicubic.java file. If you save the source under that name, you can just run javac on it.
+Why does this class exist?
+
+    I created this class because I needed to evaluate a fairly expensive function over x,y many times. It turned out that I was able to rewrite the function as a sum of a "low frequency" part that only varied slowly, and a much less expensive "high frequency" part that just contained the details.
+
+    I used the bicubic patches as follows: First I presampled the low frequency part on a coarse grid. During the inner loop of the evaluation, I used this bicubic approximation, and then added in the high frequency part. Through this approach I was able to speed up the computation by a large factor.
+
+What does the class do?
+
+    If you provide a 4x4 grid of values (ie: all the corners in a 3x3 arrangement of square tiles), this class creates an object that will interpolate a Catmull-Rom spline to give you the value within any point of the center tile.
+
+    If you want to create a spline surface, you can make a two dimensional array of such objects, arranged so that adjoining objects overlap by one grid point.
+
+The class provides a constructor and a method:
+
+    Bicubic(double[][] G)
+
+        Given 16 values on the grid [-1,0,1,2] X [-1,0,1,2], calculate bicubic coefficients.
+
+    double eval(double x, double y)
+
+        Given a point in the square [0,1] X [0,1], return a value.
+
+Algorithm:
+
+    f(x,y) = X * M * G * MT * YT , where:
+
+        X = (x3 x2 x 1) ,
+        Y = (y3 y2 y 1) ,
+        M is the Catmull-Rom basis matrix.
+
+    The constructor Bicubic() calculates the matrix C = M * G * MT
+
+    The method eval(x,y) calculates the value X * C * YT
+
+*/
+
+package uk.me.parabola.mkgmap.reader.hgt;
+
+public class Bicubic
+{
+   static double[][] M = {             // Catmull-Rom basis matrix
+      {-0.5,  1.5, -1.5, 0.5},
+      { 1  , -2.5,  2  ,-0.5},
+      {-0.5,  0  ,  0.5, 0  },
+      { 0  ,  1  ,  0  , 0  }
+   };
+
+   double[][] C = new double[4][4];    // coefficients matrix
+
+   Bicubic(double[][] G) {
+      double[][] T = new double[4][4];
+
+      for (int i = 0 ; i < 4 ; i++)    // T = G * MT
+      for (int j = 0 ; j < 4 ; j++)
+      for (int k = 0 ; k < 4 ; k++)
+	 T[i][j] += G[i][k] * M[j][k];
+
+      for (int i = 0 ; i < 4 ; i++)    // C = M * T
+      for (int j = 0 ; j < 4 ; j++)
+      for (int k = 0 ; k < 4 ; k++)
+	 C[i][j] += M[i][k] * T[k][j];
+   }
+
+   double[] X3 = C[0], X2 = C[1], X1 = C[2], X0 = C[3];
+
+   public double eval(double x, double y) {
+      return x*(x*(x*(y * (y * (y * X3[0] + X3[1]) + X3[2]) + X3[3])
+                   + (y * (y * (y * X2[0] + X2[1]) + X2[2]) + X2[3]))
+                +    (y * (y * (y * X1[0] + X1[1]) + X1[2]) + X1[3]))
+             +       (y * (y * (y * X0[0] + X0[1]) + X0[2]) + X0[3]);
+   }
+}
Index: src/uk/me/parabola/mkgmap/reader/hgt/HGTConverter.java
===================================================================
--- src/uk/me/parabola/mkgmap/reader/hgt/HGTConverter.java	(revision 4070)
+++ src/uk/me/parabola/mkgmap/reader/hgt/HGTConverter.java	(working copy)
@@ -38,6 +38,10 @@
 	private int lastRow = -1;
 	private int pointsDistanceLat;
 	private int pointsDistanceLon;
+	private boolean useBicubic;
+
+	private int statPoints;
+	private int statBicubic;
 	
 	/**
 	 * Class to extract elevation information from SRTM files in hgt format.
@@ -46,14 +50,18 @@
 	 * @param demPolygonMapUnits optional bounding polygon which describes the area for 
 	 * which elevation should be read from hgt files.
 	 */
-	public HGTConverter(String path, Area bbox, java.awt.geom.Area demPolygonMapUnits) {
-		int minLat = (int) Math.floor(Utils.toDegrees(bbox.getMinLat()));
-		int minLon = (int) Math.floor(Utils.toDegrees(bbox.getMinLong()));
-		// the bbox is slightly enlarged so that we dont fail when rounding has no effect
-		// e.g. if getMaxLat() returns 0
-		int maxLat = (int) Math.ceil(Utils.toDegrees(bbox.getMaxLat() + 1));
-		int maxLon = (int) Math.ceil(Utils.toDegrees(bbox.getMaxLong() + 1));
+	public HGTConverter(String path, Area bbox, java.awt.geom.Area demPolygonMapUnits, double extra) {
+		// make bigger box for interpolation or aligning of areas
+		int minLat = (int) Math.floor(Utils.toDegrees(bbox.getMinLat()) - extra);
+		int minLon = (int) Math.floor(Utils.toDegrees(bbox.getMinLong()) - extra);
+		int maxLat = (int) Math.ceil(Utils.toDegrees(bbox.getMaxLat()) + extra);
+		int maxLon = (int) Math.ceil(Utils.toDegrees(bbox.getMaxLong()) + extra);
 
+		if (minLat < -90) minLat = -90;
+		if (maxLat > 90) maxLat = 90;
+		if (minLon < -180) minLon = -180;
+		if (maxLon > 180) maxLon = 180;
+
 		minLat32 = Utils.toMapUnit(minLat) * 256; 
 		minLon32 = Utils.toMapUnit(minLon) * 256; 
 		// create matrix of readers 
@@ -100,27 +108,41 @@
 		rdr.prepRead();
 		if (res <= 0)
 			return 0; // assumed to be an area in the ocean
+		lastRow = row;
+
 		double scale  = res * FACTOR;
 		
-		
 		double y1 = (lat32 - minLat32) * scale - row * res;
 		double x1 = (lon32 - minLon32) * scale - col * res;
 		int xLeft = (int) x1;
 		int yBottom = (int) y1;
-		int xRight = xLeft + 1;
-		int yTop = yBottom + 1;
-		
-		int hLT = rdr.ele(xLeft, yTop);
-		int hRT = rdr.ele(xRight, yTop);
-		int hLB = rdr.ele(xLeft, yBottom);
-		int hRB = rdr.ele(xRight, yBottom);
-		lastRow = row;
-		
-		short h = interpolatedHeight(x1-xLeft, y1-yBottom, hLT, hRT, hRB, hLB);
-//		System.out.println(String.format("%.7f %.7f: (%.2f) %d", lat,lon,rc,rc2));
-//		if (Math.round(rc2) != (short) rc2) {
-//			long dd = 4;
-//		}
+
+		short h = HGTReader.UNDEF;
+
+		statPoints++;
+		if (useBicubic) {
+			// bicubic interpolation with 16 points
+			double[][] eleArray = fillArray(rdr, row, col, xLeft, yBottom);
+			if (eleArray != null) {
+				Bicubic interpolator = new Bicubic(eleArray);
+				h = (short) Math.round(interpolator.eval(x1-xLeft, y1-yBottom));
+				statBicubic++;
+			}
+		}
+
+		if (h == HGTReader.UNDEF) {
+			// use bilinear interpolation if bicubic not available
+			int xRight = xLeft + 1;
+			int yTop = yBottom + 1;
+
+			int hLT = rdr.ele(xLeft, yTop);
+			int hRT = rdr.ele(xRight, yTop);
+			int hLB = rdr.ele(xLeft, yBottom);
+			int hRB = rdr.ele(xRight, yBottom);
+			
+			h = interpolatedHeight(x1-xLeft, y1-yBottom, hLT, hRT, hRB, hLB);
+		}
+
 		if (h == HGTReader.UNDEF && log.isLoggable(Level.WARNING)) {
 			double lon = lon32 * FACTOR;
 			double lat = lat32 * FACTOR;
@@ -129,7 +151,275 @@
 		}
 		return h;
 	}
+
+	/**
+	 * Fill 16 values of HGT near required coordinates
+	 * can use HGTreaders near the current one
+	 */
+	private double[][] fillArray(HGTReader rdr, int row, int col, int xLeft, int yBottom)
+	{
+		int res = rdr.getRes();
+		int minX = 0;
+		int minY = 0;
+		int maxX = 3;
+		int maxY = 3;
+		boolean inside = true;
+
+		// checking borders
+		if (xLeft == 0) {
+			if (col == 0)
+				return null;
+			minX = 1;
+			inside = false;
+		}
+		else if (xLeft == res - 1) {
+			if (col == readers.length)
+				return null;
+			maxX = 2;
+			inside = false;
+		}
+		if (yBottom == 0) {
+			if (row == 0)
+				return null;
+			minY = 1;
+			inside = false;
+		}
+		else if (yBottom == res - 1) {
+			if (row == readers[0].length)
+				return null;
+			maxY = 2;
+			inside = false;
+		}
+
+		// fill data from current reader
+		double[][] eleArray = new double[4][4];
+		short h;
+		for (int x = minX; x <= maxX; x++) {
+			for (int y = minY; y <= maxY; y++) {
+				h = rdr.ele(xLeft + x - 1, yBottom + y - 1);
+				if (h == HGTReader.UNDEF)
+					return null;
+				eleArray[x][y] = h;
+			}
+		}
+
+		if (inside)		// no need to check borders
+			return eleArray;
+
+		// fill data from adjacent readers, down and up
+		if (xLeft > 0 && xLeft < res - 1) {
+			if (yBottom == 0) {	// bottom edge
+				HGTReader rdrBB = prepReader(res, row - 1, col);
+				if (rdrBB != null) {
+					for (int x = 0; x <= 3; x++ ) {
+						h = rdrBB.ele(xLeft + x - 1, res - 1);
+						if (h == HGTReader.UNDEF)
+							return null;
+						eleArray[x][0] = h;
+					}
+				}
+				else return null;
+			}
+			else if (yBottom == res - 1) {	// top edge
+				HGTReader rdrTT = prepReader(res, row + 1, col);
+				if (rdrTT != null) {
+					for (int x = 0; x <= 3; x++ ) {
+						h = rdrTT.ele(xLeft + x - 1, 1);
+						if (h == HGTReader.UNDEF)
+							return null;
+						eleArray[x][3] = h;
+					}
+				}
+				else return null;
+			}
+		}
+
+		// fill data from adjacent readers, left and right
+		if (yBottom > 0 && yBottom < res - 1) {
+			if (xLeft == 0) {	// left edgge
+				HGTReader rdrLL = prepReader(res, row, col - 1);
+				if (rdrLL != null) {
+					for (int y = 0; y <= 3; y++ ) {
+						h = rdrLL.ele(res - 1, yBottom + y - 1);
+						if (h == HGTReader.UNDEF)
+							return null;
+						eleArray[0][y] = h;
+					}
+				}
+				else return null;
+			}
+			else if (xLeft == res - 1) {	// right edge
+				HGTReader rdrRR = prepReader(res, row, col + 1);
+				if (rdrRR != null) {
+					for (int y = 0; y <= 3; y++ ) {
+						h = rdrRR.ele(1, yBottom + y - 1);
+						if (h == HGTReader.UNDEF)
+							return null;
+						eleArray[3][y] = h;
+					}
+				}
+				else return null;
+			}
+		}
+
+		// fill data from adjacent readers, corners
+		if (xLeft == 0) {
+			if (yBottom == 0) {	// left bottom corner
+				HGTReader rdrLL = prepReader(res, row, col - 1);
+				if (rdrLL != null) {
+					for (int y = 1; y <= 3; y++ ) {
+						h = rdrLL.ele(res - 1, yBottom + y - 1);
+						if (h == HGTReader.UNDEF)
+							return null;
+						eleArray[0][y] = h;
+					}
+				}
+				else return null;
+
+				HGTReader rdrBB = prepReader(res, row - 1, col);
+				if (rdrBB != null) {
+					for (int x = 1; x <= 3; x++ ) {
+						h = rdrBB.ele(xLeft + x - 1, res - 1);
+						if (h == HGTReader.UNDEF)
+							return null;
+						eleArray[x][0] = h;
+					}
+				}
+				else return null;
+
+				HGTReader rdrLB = prepReader(res, row - 1, col -1);
+				if (rdrLB != null) {
+					h = rdrLB.ele(res - 1, res - 1);
+					if (h == HGTReader.UNDEF)
+						return null;
+					eleArray[0][0] = h;
+				}
+				else return null;
+			}
+			else if (yBottom == res - 1) {	// left top corner
+				HGTReader rdrLL = prepReader(res, row, col - 1);
+				if (rdrLL != null) {
+					for (int y = 0; y <= 2; y++ ) {
+						h = rdrLL.ele(res - 1, yBottom + y - 1);
+						if (h == HGTReader.UNDEF)
+							return null;
+						eleArray[0][y] = h;
+					}
+				}
+				else return null;
+
+				HGTReader rdrTT = prepReader(res, row + 1, col);
+				if (rdrTT != null) {
+					for (int x = 1; x <= 3; x++ ) {
+						h = rdrTT.ele(xLeft + x - 1, 1);
+						if (h == HGTReader.UNDEF)
+							return null;
+						eleArray[x][3] = h;
+					}
+				}
+				else return null;
+
+				HGTReader rdrLT = prepReader(res, row + 1, col - 1);
+				if (rdrLT != null) {
+					h = rdrLT.ele(res - 1, 1);
+					if (h == HGTReader.UNDEF)
+						return null;
+					eleArray[0][3] = h;
+				}
+				else return null;
+			}
+		}
+		else if (xLeft == res - 1) {
+			if (yBottom == 0) {	// right bottom corner
+				HGTReader rdrRR = prepReader(res, row, col + 1);
+				if (rdrRR != null) {
+					for (int y = 1; y <= 3; y++ ) {
+						h = rdrRR.ele(1, yBottom + y - 1);
+						if (h == HGTReader.UNDEF)
+							return null;
+						eleArray[3][y] = h;
+					}
+				}
+				else return null;
+
+				HGTReader rdrBB = prepReader(res, row - 1, col);
+				if (rdrBB != null) {
+					for (int x = 0; x <= 2; x++ ) {
+						h = rdrBB.ele(xLeft + x - 1, res - 1);
+						if (h == HGTReader.UNDEF)
+							return null;
+						eleArray[x][0] = h;
+					}
+				}
+				else return null;
+
+				HGTReader rdrRB = prepReader(res, row - 1, col + 1);
+				if (rdrRB != null) {
+					h = rdrRB.ele(1, res - 1);
+					if (h == HGTReader.UNDEF)
+						return null;
+					eleArray[3][0] = h;
+				}
+				else return null;
+			}
+			else if (yBottom == res - 1) {	// right top corner
+				HGTReader rdrRR = prepReader(res, row, col + 1);
+				if (rdrRR != null) {
+					for (int y = 0; y <= 2; y++ ) {
+						h = rdrRR.ele(1, yBottom + y - 1);
+						if (h == HGTReader.UNDEF)
+							return null;
+						eleArray[3][y] = h;
+					}
+				}
+				else return null;
+
+				HGTReader rdrTT = prepReader(res, row + 1, col);
+				if (rdrTT != null) {
+					for (int x = 0; x <= 2; x++ ) {
+						h = rdrTT.ele(xLeft + x - 1, 1);
+						if (h == HGTReader.UNDEF)
+							return null;
+						eleArray[x][3] = h;
+					}
+				}
+				else return null;
+
+				HGTReader rdrRT = prepReader(res, row + 1, col + 1);
+				if (rdrRT != null) {
+					h = rdrRT.ele(1, 1);
+					if (h == HGTReader.UNDEF)
+						return null;
+					eleArray[3][3] = h;
+				}
+				else return null;
+			}
+		}
+
+		// all 16 values present
+		return eleArray;
+	}
+
+	/**
+	 * 
+	 */
+	private HGTReader prepReader(int res, int row, int col) {
+		HGTReader rdr = readers[row][col];
+
+		if (rdr == null)
+			return null;
+
+		// do not use if different resolution
+		if (res != rdr.getRes())
+			return null;
+
+		rdr.prepRead();
+		if (row > lastRow)
+			lastRow = row;
 	
+		return rdr;
+	}
+
 	/**
 	 * Try to free java heap memory allocated by the readers.
 	 */
@@ -263,7 +553,20 @@
 	public void setLonDist(int pointsDistance) {
 		this.pointsDistanceLon = pointsDistance;
 	}
+	public void setBicubic(boolean bicubic) {
+		this.useBicubic = bicubic;
+	}
 
+	public void clearStat() {
+		statPoints = 0;
+		statBicubic = 0;
+	}
+	public void printStat() {
+		if (useBicubic) {
+			System.out.println("DEM points calculated bicubic " + statBicubic + ", bilinear " + (statPoints - statBicubic) +
+				", distances " + pointsDistanceLat + " " + pointsDistanceLon);
+		}
+	}
 
 	/**
 	 * Fill array with real height values for a given upper left corner of a rectangle.
