Hi,
Gerd, your additional interpolation would be redundant.
Side note: comparison of double and integer doesn't look safe. Maybe
would be better if interpolateHeight() returned integer?
Frank, I understand your idea, but please note, that your choice of
triangle is arbitrary. There are 2 ways of dividing a square and if you
would choose the other way, results would be different and not only for
the middle point. This suggest that the solution isn't complete. Classic
bilinear interpolation seems better, but differences are small and
actual map looks more or less the same.
I have attached second patch, which apply interpolation when only 2
values on any edge are present.
Another topic. I observe some DEM artifacts, when displaying area near
the border of dem-poly. There are small rectangles without shading
within a bigger tile. They appear at random, when scrolling or zooming
map in BaseCamp or Mapsource. Ctrl-G doesn't help. See attached img.
I haven't noticed the problem, when compiling a map without dem-poly.
Maybe it is a result of slanted edge of DEM clipped by poly.
--
Best regards,
Andrzej
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,70 @@
* @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) {
+ if (hrt != HGTReader.UNDEF && hlt !=
HGTReader.UNDEF && qy > 0.5D) //top edge
+ return (1.0D - qx)*hlt + qx*hrt;
+ if (hrt != HGTReader.UNDEF && hrb !=
HGTReader.UNDEF && qx > 0.5D) //right edge
+ return (1.0D - qy)*hrb + qy*hrt;
+ //if (hlt != HGTReader.UNDEF && hrb !=
HGTReader.UNDEF && qx + qy > 0.5D && gx + qy < 1.5D) //diagonal
+ // nearest value
+ return (double)((qx < 0.5D)? ((qy < 0.5D)? hlb:
hlt): ((qy < 0.5D)? hrb: hrt));
+ }
+ if (qx + qy < 0.4D) // point is near missing value
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) {
+ if (hlb != HGTReader.UNDEF && hrb !=
HGTReader.UNDEF && qy < 0.5D) //lower edge
+ return (1.0D - qx)*hlb + qx*hrb;
+ if (hlb != HGTReader.UNDEF && hlt !=
HGTReader.UNDEF && qx < 0.5D) //left edge
+ return (1.0D - qy)*hlb + qy*hlt;
+ //if (hlt != HGTReader.UNDEF && hrb !=
HGTReader.UNDEF && qx + qy > 0.5D && gx + qy < 1.5D) //diagonal
+ // nearest value
+ return (double)((qx < 0.5D)? ((qy < 0.5D)? hlb:
hlt): ((qy < 0.5D)? hrb: hrt));
+ }
+ if (qx + qy > 1.6D) // point is near missing value
return HGTReader.UNDEF;
+ hrt = hlt + hrb - hlb;
+ }
+ else if (hrb == HGTReader.UNDEF) {
+ if (hlb == HGTReader.UNDEF || hlt == HGTReader.UNDEF ||
hrt == HGTReader.UNDEF) {
+ if (hlt != HGTReader.UNDEF && hrt !=
HGTReader.UNDEF && qy > 0.5D) //top edge
+ return (1.0D - qx)*hlt + qx*hrt;
+ if (hlt != HGTReader.UNDEF && hlb !=
HGTReader.UNDEF && qx < 0.5D) //left edge
+ return (1.0D - qy)*hlb + qy*hlt;
+ //if (hlb != HGTReader.UNDEF && hrt !=
HGTReader.UNDEF && qy > qx - 0.5D && qy < qx + 0.5D) //diagonal
+ // nearest value
+ return (double)((qx < 0.5D)? ((qy < 0.5D)? hlb:
hlt): ((qy < 0.5D)? hrb: hrt));
+ }
+ if (qy < qx - 0.4D) // point is near missing value
+ return HGTReader.UNDEF;
+ hrb = hlb + hrt - hlt;
+ }
+ else if (hlt == HGTReader.UNDEF) {
+ if (hlb == HGTReader.UNDEF || hrb == HGTReader.UNDEF ||
hrt == HGTReader.UNDEF) {
+ if (hrb != HGTReader.UNDEF && hlb !=
HGTReader.UNDEF && qy < 0.5D) //lower edge
+ return (1.0D - qx)*hlb + qx*hrb;
+ if (hrb != HGTReader.UNDEF && hrt !=
HGTReader.UNDEF && qx > 0.5D) //right edge
+ return (1.0D - qy)*hrb + qy*hrt;
+ //if (hlb != HGTReader.UNDEF && hrt !=
HGTReader.UNDEF && qy > qx - 0.5D && qy < qx + 0.5D) //diagonal
+ // nearest value
+ return (double)((qx < 0.5D)? ((qy < 0.5D)? hlb:
hlt): ((qy < 0.5D)? hrb: hrt));
+ }
+ if (qy > qx + 0.6D) // point is near missing value
+ 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() {
_______________________________________________
mkgmap-dev mailing list
mkgmap-dev@lists.mkgmap.org.uk
http://www.mkgmap.org.uk/mailman/listinfo/mkgmap-dev
