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
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
_______________________________________________
mkgmap-dev mailing list
mkgmap-dev@lists.mkgmap.org.uk
http://www.mkgmap.org.uk/mailman/listinfo/mkgmap-dev
