[mkgmap-dev] HGT - getElevation()
From Andrzej Popowski popej at poczta.onet.pl on Sun Jan 21 17:38:05 GMT 2018
Hi, I'm attaching a next patch for bicubic interpolation. It include following changes: - some errors corrected, - more statistic of interpolation, - changed condition for applying bicubic interpolation. In this version bicubic interpolation is applied, when required DEM resolution is not lower than 1/3 of HGT resolution. I assume, that for low DEM resolution (big dem-dists values), there is no reason to make precise interpolation. I have uploaded compiled mkgmap: http://files.mkgmap.org.uk/download/404/mkgmap-dem-tdb-r4071-bicubic-6.zip -- 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 4071) +++ src/uk/me/parabola/imgfmt/app/dem/DEMFile.java (working copy) @@ -52,20 +52,56 @@ * @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) Math.round((1 << 29) / (res * 45.0D)); + } + // 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; + } + } + + 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 4071) +++ src/uk/me/parabola/imgfmt/app/dem/DEMSection.java (working copy) @@ -48,24 +48,23 @@ 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 pointsDistanceLat = pointDist; pointsDistanceLon = pointDist; - - int []latInfo = getTileInfo(bbox.getHeight() * 256, pointsDistanceLat); - int []lonInfo = getTileInfo(bbox.getWidth() * 256, pointsDistanceLon); + + // select bicubic or bilinear interpolation + hgtConverter.setBicubic(zoomLevel, 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 +112,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 +150,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 4071) +++ src/uk/me/parabola/mkgmap/reader/hgt/HGTConverter.java (working copy) @@ -38,7 +38,16 @@ private int lastRow = -1; private int pointsDistanceLat; private int pointsDistanceLon; - + private boolean useBicubic; + + private int statPoints; + private int statBicubic; + private int statBilinear; + private int statVoid; + private int statRdrNull; + private int statRdrRes; + + /** * Class to extract elevation information from SRTM files in hgt format. * @param path a comma separated list of directories which may contain *.hgt files. @@ -46,14 +55,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 +113,43 @@ 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); + statBilinear++; + if (h == HGTReader.UNDEF) statVoid++; + } + if (h == HGTReader.UNDEF && log.isLoggable(Level.WARNING)) { double lon = lon32 * FACTOR; double lat = lat32 * FACTOR; @@ -129,7 +158,263 @@ } 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; + + // check borders + if (xLeft == 0) { + if (col <= 0) + return null; + minX = 1; + inside = false; + } + else if (xLeft == res - 1) { + if (col >= readers[0].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.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 again + 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) + return 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 if (yBottom == res - 1) { // top edge + HGTReader rdrTT = prepReader(res, row + 1, col); + if (rdrTT == null) + return 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; + } + } + } + + // 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) + return 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 if (xLeft == res - 1) { // right edge + HGTReader rdrRR = prepReader(res, row, col + 1); + if (rdrRR == null) + return 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; + } + } + } + + // 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) + return 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; + } + + HGTReader rdrBB = prepReader(res, row - 1, col); + if (rdrBB == null) + return 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; + } + + HGTReader rdrLB = prepReader(res, row - 1, col -1); + if (rdrLB == null) + return null; + h = rdrLB.ele(res - 1, res - 1); + if (h == HGTReader.UNDEF) + return null; + eleArray[0][0] = h; + } + else if (yBottom == res - 1) { // left top corner + HGTReader rdrLL = prepReader(res, row, col - 1); + if (rdrLL == null) + return 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; + } + + HGTReader rdrTT = prepReader(res, row + 1, col); + if (rdrTT == null) + return 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; + } + + HGTReader rdrLT = prepReader(res, row + 1, col - 1); + if (rdrLT == null) + return null; + h = rdrLT.ele(res - 1, 1); + if (h == HGTReader.UNDEF) + return null; + eleArray[0][3] = h; + } + } + else if (xLeft == res - 1) { + if (yBottom == 0) { // right bottom corner + HGTReader rdrRR = prepReader(res, row, col + 1); + if (rdrRR == null) + return 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; + } + + HGTReader rdrBB = prepReader(res, row - 1, col); + if (rdrBB == null) + return 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; + } + + HGTReader rdrRB = prepReader(res, row - 1, col + 1); + if (rdrRB == null) + return null; + h = rdrRB.ele(1, res - 1); + if (h == HGTReader.UNDEF) + return null; + eleArray[3][0] = h; + } + else if (yBottom == res - 1) { // right top corner + HGTReader rdrRR = prepReader(res, row, col + 1); + if (rdrRR == null) + return 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; + } + + HGTReader rdrTT = prepReader(res, row + 1, col); + if (rdrTT == null) + return 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; + } + + HGTReader rdrRT = prepReader(res, row + 1, col + 1); + if (rdrRT == null) + return null; + h = rdrRT.ele(1, 1); + if (h == HGTReader.UNDEF) + return null; + eleArray[3][3] = h; + } + } + + // all 16 values present + return eleArray; + } + + /** + * + */ + private HGTReader prepReader(int res, int row, int col) { + HGTReader rdr = readers[row][col]; + + if (rdr == null) { + statRdrNull++; + return null; + } + + // do not use if different resolution + if (res != rdr.getRes()) { + statRdrRes++; + return null; + } + + rdr.prepRead(); + if (row > lastRow) + lastRow = row; + return rdr; + } + /** * Try to free java heap memory allocated by the readers. */ @@ -264,7 +549,41 @@ this.pointsDistanceLon = pointsDistance; } + /** + * Estimate if bicubic interpolatin is useful + * @param zoomLevel number of current DEM layer + * @param pointDist distance between DEM points + */ + public void setBicubic(int zoomLevel, int pointDist) { + if (res > 0) { + //if (zoomLevel == 0) { + // this.useBicubic = true; // always use bicubic for most detailed layer, for overview too + // return; + //} + int distHGTx3 = (1 << 29)/(45 / 3 * res); // 3 * HGT points distance in DEM units + if (distHGTx3 + 20 > pointDist) { // account for rounding of pointDist and distHGTx3 + this.useBicubic = true; // DEM resolution is greater than 1/3 HGT resolution + return; + } + } + this.useBicubic = false; + } + public void clearStat() { + statPoints = 0; + statBicubic = 0; + statBilinear = 0; + statVoid = 0; + statRdrNull = 0; + statRdrRes = 0; + } + public void printStat() { + if (useBicubic) { + System.out.println("DEM points: " + statPoints + "; bicubic " + statBicubic + ", no HGT " + (statRdrNull + statRdrRes) + + "; bilinear " + statBilinear + ", voids " + statVoid + "; distance " + pointsDistanceLat); + } + } + /** * Fill array with real height values for a given upper left corner of a rectangle. * @param lat32 latitude of upper left corner
- Previous message: [mkgmap-dev] HGT - getElevation()
- Next message: [mkgmap-dev] HGT - getElevation()
- Messages sorted by: [ date ] [ thread ] [ subject ] [ author ]
More information about the mkgmap-dev mailing list