+;--dem-interpolation=value +: Use this option to speciy the method that is used to interpolate data from +hgt raster to the DEM raster. The value bicubic gives the highest precision +but is slower, bilinear is faster but less precise, it tends to smooth the +profile and thus also reduces DEM size compared to bicubic. The value auto +means that bicubic is used where is seems appropriate according to hgt +resolution and dem-dist value, else bilinear is used. The default is auto. +
;--dem-poly=filename : If given, the filename should point to a *.poly file in osmosis polygon file format. The polygon described in the file is used to determine the area Index: resources/help/en/options =================================================================== --- resources/help/en/options (revision 4086) +++ resources/help/en/options (working copy) @@ -463,7 +463,16 @@ Example which should work with levels="0:24, 1:22, 2:20, 3:18": --dem-dists=3312,13248,26512,53024 This was found in a Garmin Demo map for transalpin data created 2009. - + +--dem-interpolation=value + Use this option to speciy the method that is used to interpolate + data from hgt raster to the DEM raster. The value bicubic gives the + highest precision but is slower, bilinear is faster but less precise, + it tends to smooth the profile and thus also reduces DEM size compared to + bicubic. The value auto means that bicubic is used where is seems + appropriate according to hgt resolution and dem-dist value, else bilinear + is used. If not given, the default is auto. + --dem-poly=filename If given, the filename should point to a *.poly file in osmosis polygon file format. The polygon described in the file is used to determine the area Index: src/uk/me/parabola/imgfmt/app/dem/DEMSection.java =================================================================== --- src/uk/me/parabola/imgfmt/app/dem/DEMSection.java (revision 4086) +++ src/uk/me/parabola/imgfmt/app/dem/DEMSection.java (working copy) @@ -59,8 +59,8 @@ pointsDistanceLat = pointDist; pointsDistanceLon = pointDist; - // select bicubic or bilinear interpolation - hgtConverter.setBicubic(zoomLevel, pointDist); + // allow automatic selection of interpolation method for each zoom level + hgtConverter.startNewLevel(pointDist); int []latInfo = getTileInfo(areaHeight, pointsDistanceLat); int []lonInfo = getTileInfo(areaWidth, pointsDistanceLon); @@ -109,9 +109,6 @@ int minBaseHeight = Integer.MAX_VALUE; int maxBaseHeight = Integer.MIN_VALUE; int maxDeltaHeight = Integer.MIN_VALUE; - hgtConverter.setLatDist(pointsDistanceLat); - hgtConverter.setLonDist(pointsDistanceLon); - hgtConverter.clearStat(); for (int m = 0; m < tilesLat; m++) { latOff = top - m * resLat; Index: src/uk/me/parabola/mkgmap/build/MapBuilder.java =================================================================== --- src/uk/me/parabola/mkgmap/build/MapBuilder.java (revision 4086) +++ src/uk/me/parabola/mkgmap/build/MapBuilder.java (working copy) @@ -219,19 +219,20 @@ if (demPolygonFile != null) { demPolygon = Java2DConverter.readPolyFile(demPolygonFile); } - String ipm = props.getProperty("dem-interpolation", "bicubic"); + String ipm = props.getProperty("dem-interpolation", "auto"); switch (ipm) { + case "auto": + demInterpolationMethod = InterpolationMethod.Automatic; + break; case "bicubic": demInterpolationMethod = InterpolationMethod.Bicubic; break; - case "bicubic-spline": - demInterpolationMethod = InterpolationMethod.BicubicSpline; - break; case "bilinear": demInterpolationMethod = InterpolationMethod.Bilinear; break; default: - throw new IllegalArgumentException("invalid argument for option dem-interpolation: '" + ipm + "' supported are 'bilinear', 'bicubic', or 'bicubic-spline'"); + throw new IllegalArgumentException("invalid argument for option dem-interpolation: '" + ipm + + "' supported are 'bilinear', 'bicubic', or 'auto'"); } } Index: src/uk/me/parabola/mkgmap/reader/hgt/Cubic.java =================================================================== --- src/uk/me/parabola/mkgmap/reader/hgt/Cubic.java (revision 4086) +++ src/uk/me/parabola/mkgmap/reader/hgt/Cubic.java (nonexistent) @@ -1,116 +0,0 @@ -package uk.me.parabola.mkgmap.reader.hgt; -/* -Java class to implement cubic and bicubic splines - by Ken Perlin @ NYU, 1998, 2004. - -You have my permission to use freely, as long as you keep the attribution. - Ken Perlin - -What does the class do? - - Cubic spline: If you provide various geometric values for t, then this class creates an object that will interpolate a Cubic spline to give you the value within any value of t between 0 and 1. - - If you want to create a spline path, you can make a one dimensional array of such objects. - - Bicubic spline: If you provide a 4x4 grid of values for geometric quantities in u and v, this class creates an object that will interpolate a Bicubic spline to give you the value within any point of a unit tile in (u,v) space. - - If you want to create a spline surface, you can make a two dimensional array of such objects. - -For a cubic spline the class provides a constructor and a method: - - Cubic(double[] G) - - Given four geometric values over t, calculate cubic coefficients. - - double eval(double t) - - Given a point in the interval t = [0 ... 1], return a value. - - Algorithm: - - f(t) = T M GT, where: - - T = (t3 t2 t 1) , - M is the basis matrix. - - The constructor Cubic(G) calculates the matrix C = M GT - - The method eval(t) calculates the value T C - -For a bicubic spline the class provides a constructor and a method: - - Cubic(double[][] G) - - Given 4×4 geometric values over u×v, calculate bicubic coefficients. - - double eval(double u, double v) - - Given a point in the square [0 ... 1] × [0 ... 1], return a value. - - Algorithm: - - f(u,v) = U M G MT VT , where: - - U = (u3 u2 u 1) , - V = (v3 v2 v 1) , - M is the basis matrix. - - The constructor Cubic(G) calculates the matrix C = M G MT - - The method eval(u,v) calculates the value U C VT - -*/ - -public class Cubic -{ - public static final double[][] BEZIER = { // Bezier basis matrix - {-1 , 3 , -3 , 1 }, - { 3 , -6 , 3 , 0 }, - {-3 , 3 , 0 , 0 }, - { 1 , 0 , 0 , 0 } - }; - public static final double[][] BSPLINE = { // BSpline basis matrix - {-1./6 , 3./6 , -3./6 , 1./6 }, - { 3./6 , -6./6 , 3./6 , 0. }, - {-3./6 , 0. , 3./6 , 0. }, - { 1./6 , 4./6 , 1./6 , 0. } - }; - public static final double[][] CATMULL_ROM = { // 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 } - }; - public static final double[][] HERMITE = { // Hermite basis matrix - { 2 , -2 , 1 , 1 }, - {-3 , 3 , -2 , -1 }, - { 0 , 0 , 1 , 0 }, - { 1 , 0 , 0 , 0 } - }; - - double a, b, c, d; // cubic coefficients vector - - double[][] C = new double[4][4]; // bicubic coefficients matrix - double[][] T = new double[4][4]; // scratch matrix - - Cubic(double[][] M, double[][] G) { - 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[] C3 = C[0], C2 = C[1], C1 = C[2], C0 = C[3]; - - public double eval(double u, double v) { - return u * (u * (u * (v * (v * (v * C3[0] + C3[1]) + C3[2]) + C3[3]) - + (v * (v * (v * C2[0] + C2[1]) + C2[2]) + C2[3])) - + (v * (v * (v * C1[0] + C1[1]) + C1[2]) + C1[3])) - + (v * (v * (v * C0[0] + C0[1]) + C0[2]) + C0[3]); - } -} - Index: src/uk/me/parabola/mkgmap/reader/hgt/HGTConverter.java =================================================================== --- src/uk/me/parabola/mkgmap/reader/hgt/HGTConverter.java (revision 4086) +++ src/uk/me/parabola/mkgmap/reader/hgt/HGTConverter.java (working copy) @@ -48,9 +48,13 @@ private int statRdrRes; private InterpolationMethod interpolationMethod = InterpolationMethod.Bicubic; - public enum InterpolationMethod { - Bilinear, Bicubic, BicubicSpline + /** faster, smoothing, less precise */ + Bilinear , + /** slower, higher precision */ + Bicubic, + /** bicubic for high resolution, else bilinear */ + Automatic }; /** @@ -97,8 +101,17 @@ res = maxRes; // we use the highest available res return; } - + /** + * Allows to change the interpolation method for complex interpolations. + * @param interpolationMethod + */ + public void setInterpolationMethod(InterpolationMethod interpolationMethod) { + this.interpolationMethod = interpolationMethod; + useComplexInterpolation = (interpolationMethod != InterpolationMethod.Bilinear); + } + + /** * Return elevation in meter for a point given in DEM units (32 bit res). * @param lat32 * @param lon32 @@ -105,7 +118,6 @@ * @return height in m or Short.MIN_VALUE if value is invalid */ protected short getElevation(int lat32, int lon32) { - // TODO: maybe calculate the borders in 32 bit res ? int row = (int) ((lat32 - minLat32) * FACTOR); int col = (int) ((lon32 - minLon32) * FACTOR); @@ -132,16 +144,11 @@ short h = HGTReader.UNDEF; statPoints++; - if (useComplexInterpolation && !InterpolationMethod.Bilinear.equals(interpolationMethod)) { - // bicubic or bicubic spline interpolation with 16 points + if (useComplexInterpolation) { + // bicubic (Catmull-Rom) interpolation with 16 points eleArray = fillArray(rdr, row, col, xLeft, yBottom); if (eleArray != null) { - if (InterpolationMethod.Bicubic.equals(interpolationMethod)) { - h = (short) Math.round(bicubicInterpolation(eleArray, qx, qy)); - } else { - Cubic interpolator = new Cubic(Cubic.BSPLINE, eleArray); - h = (short) Math.round(interpolator.eval(qx, qy)); - } + h = (short) Math.round(bicubicInterpolation(eleArray, qx, qy)); statBicubic++; } } @@ -174,8 +181,7 @@ * 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) - { + private double[][] fillArray(HGTReader rdr, int row, int col, int xLeft, int yBottom) { int res = rdr.getRes(); int minX = 0; int minY = 0; @@ -189,8 +195,7 @@ return null; minX = 1; inside = false; - } - else if (xLeft == res - 1) { + } else if (xLeft == res - 1) { if (col >= readers[0].length) return null; maxX = 2; @@ -201,8 +206,7 @@ return null; minY = 1; inside = false; - } - else if (yBottom == res - 1) { + } else if (yBottom == res - 1) { if (row >= readers.length) return null; maxY = 2; @@ -226,22 +230,21 @@ // fill data from adjacent readers, down and up if (xLeft > 0 && xLeft < res - 1) { - if (yBottom == 0) { // bottom edge + if (yBottom == 0) { // bottom edge HGTReader rdrBB = prepReader(res, row - 1, col); if (rdrBB == null) return null; - for (int x = 0; x <= 3; x++ ) { + 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 + } 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++ ) { + for (int x = 0; x <= 3; x++) { h = rdrTT.ele(xLeft + x - 1, 1); if (h == HGTReader.UNDEF) return null; @@ -252,22 +255,21 @@ // fill data from adjacent readers, left and right if (yBottom > 0 && yBottom < res - 1) { - if (xLeft == 0) { // left edgge + if (xLeft == 0) { // left edgge HGTReader rdrLL = prepReader(res, row, col - 1); if (rdrLL == null) return null; - for (int y = 0; y <= 3; y++ ) { + 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 + } 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++ ) { + for (int y = 0; y <= 3; y++) { h = rdrRR.ele(1, yBottom + y - 1); if (h == HGTReader.UNDEF) return null; @@ -278,11 +280,11 @@ // fill data from adjacent readers, corners if (xLeft == 0) { - if (yBottom == 0) { // left bottom corner + 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++ ) { + for (int y = 1; y <= 3; y++) { h = rdrLL.ele(res - 1, yBottom + y - 1); if (h == HGTReader.UNDEF) return null; @@ -292,7 +294,7 @@ HGTReader rdrBB = prepReader(res, row - 1, col); if (rdrBB == null) return null; - for (int x = 1; x <= 3; x++ ) { + for (int x = 1; x <= 3; x++) { h = rdrBB.ele(xLeft + x - 1, res - 1); if (h == HGTReader.UNDEF) return null; @@ -299,7 +301,7 @@ eleArray[x][0] = h; } - HGTReader rdrLB = prepReader(res, row - 1, col -1); + HGTReader rdrLB = prepReader(res, row - 1, col - 1); if (rdrLB == null) return null; h = rdrLB.ele(res - 1, res - 1); @@ -306,12 +308,11 @@ if (h == HGTReader.UNDEF) return null; eleArray[0][0] = h; - } - else if (yBottom == res - 1) { // left top corner + } 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++ ) { + for (int y = 0; y <= 2; y++) { h = rdrLL.ele(res - 1, yBottom + y - 1); if (h == HGTReader.UNDEF) return null; @@ -321,7 +322,7 @@ HGTReader rdrTT = prepReader(res, row + 1, col); if (rdrTT == null) return null; - for (int x = 1; x <= 3; x++ ) { + for (int x = 1; x <= 3; x++) { h = rdrTT.ele(xLeft + x - 1, 1); if (h == HGTReader.UNDEF) return null; @@ -336,13 +337,12 @@ return null; eleArray[0][3] = h; } - } - else if (xLeft == res - 1) { - if (yBottom == 0) { // right bottom corner + } 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++ ) { + for (int y = 1; y <= 3; y++) { h = rdrRR.ele(1, yBottom + y - 1); if (h == HGTReader.UNDEF) return null; @@ -352,7 +352,7 @@ HGTReader rdrBB = prepReader(res, row - 1, col); if (rdrBB == null) return null; - for (int x = 0; x <= 2; x++ ) { + for (int x = 0; x <= 2; x++) { h = rdrBB.ele(xLeft + x - 1, res - 1); if (h == HGTReader.UNDEF) return null; @@ -366,12 +366,11 @@ if (h == HGTReader.UNDEF) return null; eleArray[3][0] = h; - } - else if (yBottom == res - 1) { // right top corner + } 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++ ) { + for (int y = 0; y <= 2; y++) { h = rdrRR.ele(1, yBottom + y - 1); if (h == HGTReader.UNDEF) return null; @@ -381,7 +380,7 @@ HGTReader rdrTT = prepReader(res, row + 1, col); if (rdrTT == null) return null; - for (int x = 0; x <= 2; x++ ) { + for (int x = 0; x <= 2; x++) { h = rdrTT.ele(xLeft + x - 1, 1); if (h == HGTReader.UNDEF) return null; @@ -515,14 +514,14 @@ return (short) Math.round((1.0D - qy) * hxb + qy * hxt); } - public void stats() { - for (int i = 0; i < readers.length; i++) { - for (int j = 0; j < readers[i].length; j++) { - System.out.println(readers[i][j]); - } - - } - } +// public void stats() { +// for (int i = 0; i < readers.length; i++) { +// for (int j = 0; j < readers[i].length; j++) { +// System.out.println(readers[i][j]); +// } +// +// } +// } public int getHighestRes() { return res; @@ -553,34 +552,27 @@ noHeights[0]= outsidePolygonHeight; } - public void setLatDist(int pointsDistance) { - this.pointsDistanceLat = pointsDistance; - } - public void setLonDist(int pointsDistance) { - this.pointsDistanceLon = pointsDistance; - } - /** - * Estimate if bicubic interpolation is useful - * @param zoomLevel number of current DEM layer + * If interpolation method is automatic, decide which to use for a new zoom level. * @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.useComplexInterpolation = true; // DEM resolution is greater than 1/3 HGT resolution - return; + public void startNewLevel(int pointDist) { + clearStat(); + pointsDistanceLat = pointDist; + pointsDistanceLon = pointDist; + if (InterpolationMethod.Automatic.equals(interpolationMethod)) { + if (res > 0) { + 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.useComplexInterpolation = true; // DEM resolution is greater than 1/3 HGT resolution + return; + } } + this.useComplexInterpolation = false; } - this.useComplexInterpolation = false; } - public void clearStat() { + private void clearStat() { statPoints = 0; statBicubic = 0; statBilinear = 0; @@ -589,10 +581,8 @@ statRdrRes = 0; } public void printStat() { - if (useComplexInterpolation) { - log.info("DEM points: " + statPoints + "; bicubic " + statBicubic + ", no HGT " + (statRdrNull + statRdrRes) + + log.info("DEM points: " + statPoints + "; bicubic " + statBicubic + ", no HGT " + (statRdrNull + statRdrRes) + "; bilinear " + statBilinear + ", voids " + statVoid + "; distance " + pointsDistanceLat); - } } /** @@ -605,7 +595,6 @@ * a an array with one value for each point, in the order top -> down and left -> right. */ public short[] getHeights(int lat32, int lon32, int height, int width) { - // TODO Auto-generated method stub short[] realHeights = noHeights; java.awt.geom.Area testArea = null; @@ -650,6 +639,7 @@ /** * Cubic interpolation for 4 points, taken from http://www.paulinternet.nl/?page=bicubic + * Uses Catmull–Rom spline. * @author Paul Breeuwsma */ private static double cubicInterpolation(double[] p, double qx) { @@ -673,7 +663,4 @@ return cubicInterpolation(arr, qx); } - public void setInterpolationMethod(InterpolationMethod interpolationMethod) { - this.interpolationMethod = interpolationMethod; - } }