Interpolating DEM data.

- DEM data can now be interpolated using ione of 3 methods:
  o no_interpolation: Take the nearest south-west data.
  o simple_interpolation: Interpolated using a simple method.
  o best_interpolation: Interpolated using Shepard method.

- "Apply DEM data" for a track now uses "best_interpolation".

- Elevation-distance graph now uses "best_interpolation".

- Elevation information displayed on status bar now uses most suitable
  interpolation methods according to the current zoom level.

Not much testing has been made to justify performance/accuracy trade-off
of the interpolating methods.

Signed-off-by: Quy Tonthat <qtonthat@gmail.com>

diff --git a/ChangeLog b/ChangeLog
index e3c1a4f918846a76090663d125a3ba50b6233b01..7d8263d6acd3e3757306752ca3842356612af324 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,9 @@
+2007-10-16
+Quy Tonthat <qtonthat@gmail.com>:
+       * Interpolating DEM data in 3 different methods. Track data,
+       elevation-distance graph and elevation info on status bar now make use
+       of interpolation.
+
 2007-10-15
 Quy Tonthat <qtonthat@gmail.com>:
        * Add more room to the top of elevation-distance graph.

diff --git a/src/dem.c b/src/dem.c
index 6fa862fc82452fc1b6ac573115a26f7fd43fa35d..a3f13ee3e4bd487093bf70968aa1e1f6c02a1f64 100644 (file)
--- a/src/dem.c
+++ b/src/dem.c
@@ -485,6 +485,77 @@ gint16 vik_dem_get_east_north ( VikDEM *dem, gdouble east, gdouble north )
   return vik_dem_get_xy ( dem, col, row );
 }
 
+static gboolean dem_get_ref_points_elev_dist(VikDEM *dem,
+    gdouble east, gdouble north, /* in seconds */
+    gint16 *elevs, gint16 *dists)
+{
+  int i;
+  int cols[4], rows[4];
+  struct LatLon ll[4];
+  struct LatLon pos;
+
+  if ( east > dem->max_east || east < dem->min_east ||
+      north > dem->max_north || north < dem->min_north )
+    return FALSE;  /* got nothing */
+
+  pos.lon = east/3600;
+  pos.lat = north/3600;
+
+  /* order of the data: sw, nw, ne, se */
+  cols[0] = (gint) floor((east - dem->min_east) / dem->east_scale);  /* sw */
+  rows[0] = (gint) floor((north - dem->min_north) / dem->north_scale);
+  ll[0].lon = (east - fabs(fmod(east, dem->east_scale)))/3600;
+  ll[0].lat = (north - fabs(fmod(north, dem->north_scale)))/3600;
+
+  cols[1] = cols[0];         /*nw*/
+  rows[1] = rows[0] + 1;
+  ll[1].lon = ll[0].lon;
+  ll[1].lat = ll[0].lat + (gdouble)dem->north_scale/3600;
+
+  cols[2] = cols[0] + 1;     /*ne*/
+  rows[2] = rows[0] + 1;
+  ll[2].lon = ll[0].lon + (gdouble)dem->east_scale/3600;
+  ll[2].lat = ll[0].lat + (gdouble)dem->north_scale/3600;
+
+  cols[3] = cols[0] + 1;     /*se*/
+  rows[3] = rows[0];
+  ll[3].lon = ll[0].lon + (gdouble)dem->east_scale/3600;
+  ll[3].lat = ll[0].lat;
+
+  for (i = 0; i < 4; i++) {
+    if ((elevs[i] = vik_dem_get_xy(dem, cols[i], rows[i])) == VIK_DEM_INVALID_ELEVATION)
+      return FALSE;
+    dists[i] = a_coords_latlon_diff(&pos, &ll[i]);
+  }
+
+  return TRUE;  /* all OK */
+}
+
+gint16 vik_dem_get_simple_interpol ( VikDEM *dem, gdouble east, gdouble north )
+{
+  int i;
+  gint16 elevs[4], dists[4];
+
+  if (!dem_get_ref_points_elev_dist(dem, east, north, elevs, dists))
+    return VIK_DEM_INVALID_ELEVATION;
+
+  for (i = 0; i < 4; i++) {
+    if (dists[i] < 1) {
+      return(elevs[i]);
+    }
+  }
+
+  gdouble t = (gdouble)elevs[0]/dists[0] + (gdouble)elevs[1]/dists[1] + (gdouble)elevs[2]/dists[2] + (gdouble)elevs[3]/dists[3];
+  gdouble b = 1.0/dists[0] + 1.0/dists[1] + 1.0/dists[2] + 1.0/dists[3];
+
+  return(t/b);
+}
+
+gint16 vik_dem_get_shepard_interpol ( VikDEM *dem, gdouble east, gdouble north )
+{
+  return vik_dem_get_simple_interpol(dem, east, north);
+}
+
 void vik_dem_east_north_to_xy ( VikDEM *dem, gdouble east, gdouble north, guint *col, guint *row )
 {
   *col = (gint) floor((east - dem->min_east) / dem->east_scale);

diff --git a/src/dem.h b/src/dem.h
index 998ab088f1cdfa1a341314fa7783d9519efb4b3c..85964d1ce9d3a3fa1b285f94e0d1a92a5be5619c 100644 (file)
--- a/src/dem.h
+++ b/src/dem.h
@@ -45,6 +45,8 @@ void vik_dem_free ( VikDEM *dem );
 gint16 vik_dem_get_xy ( VikDEM *dem, guint x, guint y );
 
 gint16 vik_dem_get_east_north ( VikDEM *dem, gdouble east, gdouble north );
+gint16 vik_dem_get_simple_interpol ( VikDEM *dem, gdouble east, gdouble north );
+gint16 vik_dem_get_best_interpol ( VikDEM *dem, gdouble east, gdouble north );
 
 void vik_dem_east_north_to_xy ( VikDEM *dem, gdouble east, gdouble north, guint *col, guint *row );
 

diff --git a/src/dems.c b/src/dems.c
index b278e5183b2a09ab72489829d0dbeddd5a8275c6..da268311318f7b227117649d3ee93a78090c3e11 100644 (file)
--- a/src/dems.c
+++ b/src/dems.c
@@ -161,32 +161,46 @@ gint16 a_dems_list_get_elev_by_coord ( GList *dems, const VikCoord *coord )
 
 typedef struct {
   const VikCoord *coord;
+  VikDemInterpol method;
   gint elev;
 } CoordElev;
 
 static gboolean get_elev_by_coord(gpointer key, LoadedDEM *ldem, CoordElev *ce)
 {
   VikDEM *dem = ldem->dem;
+  gdouble lat, lon;
 
   if ( dem->horiz_units == VIK_DEM_HORIZ_LL_ARCSECONDS ) {
     struct LatLon ll_tmp;
     vik_coord_to_latlon (ce->coord, &ll_tmp );
-    ll_tmp.lat *= 3600;
-    ll_tmp.lon *= 3600;
-    ce->elev = vik_dem_get_east_north(dem, ll_tmp.lon, ll_tmp.lat);
-    return (ce->elev != VIK_DEM_INVALID_ELEVATION);
+    lat = ll_tmp.lat * 3600;
+    lon = ll_tmp.lon * 3600;
   } else if (dem->horiz_units == VIK_DEM_HORIZ_UTM_METERS) {
     static struct UTM utm_tmp;
+    if (utm_tmp.zone != dem->utm_zone)
+      return FALSE;
     vik_coord_to_utm (ce->coord, &utm_tmp);
-    if ( utm_tmp.zone == dem->utm_zone &&
-             (ce->elev = vik_dem_get_east_north(dem, utm_tmp.easting, utm_tmp.northing)) != VIK_DEM_INVALID_ELEVATION )
-      return TRUE;
+    lat = utm_tmp.northing;
+    lon = utm_tmp.easting;
+  } else
+    return FALSE;
+
+  switch (ce->method) {
+    case VIK_DEM_INTERPOL_NONE:
+      ce->elev = vik_dem_get_east_north(dem, lon, lat);
+      break;
+    case VIK_DEM_INTERPOL_SIMPLE:
+      ce->elev = vik_dem_get_simple_interpol(dem, lon, lat);
+      break;
+    case VIK_DEM_INTERPOL_BEST:
+      ce->elev = vik_dem_get_shepard_interpol(dem, lon, lat);
+      break;
   }
-  return FALSE;
+  return (ce->elev != VIK_DEM_INVALID_ELEVATION);
 }
 
 /* TODO: keep a (sorted) linked list of DEMs and select the best resolution one */
-gint16 a_dems_get_elev_by_coord ( const VikCoord *coord )
+gint16 a_dems_get_elev_by_coord ( const VikCoord *coord, VikDemInterpol method )
 {
   CoordElev ce;
 
@@ -194,6 +208,7 @@ gint16 a_dems_get_elev_by_coord ( const VikCoord *coord )
     return VIK_DEM_INVALID_ELEVATION;
 
   ce.coord = coord;
+  ce.method = method;
   ce.elev = VIK_DEM_INVALID_ELEVATION;
 
   if(!g_hash_table_find(loaded_dems, get_elev_by_coord, &ce))

diff --git a/src/dems.h b/src/dems.h
index cf8b1a9f436274049148db9f2d261641742f4e09..0d93244d9bf15a71fd6ce9402e7a6f3c43945c71 100644 (file)
--- a/src/dems.h
+++ b/src/dems.h
@@ -4,6 +4,12 @@
 #include "dem.h"
 #include "vikcoord.h"
 
+typedef enum {
+  VIK_DEM_INTERPOL_NONE = 0,
+  VIK_DEM_INTERPOL_SIMPLE,
+  VIK_DEM_INTERPOL_BEST,
+} VikDemInterpol;
+
 void a_dems_uninit ();
 VikDEM *a_dems_load(const gchar *filename);
 void a_dems_unref(const gchar *filename);
@@ -12,7 +18,7 @@ void a_dems_load_list ( GList **dems );
 void a_dems_list_free ( GList *dems );
 GList *a_dems_list_copy ( GList *dems );
 gint16 a_dems_list_get_elev_by_coord ( GList *dems, const VikCoord *coord );
-gint16 a_dems_get_elev_by_coord ( const VikCoord *coord );
+gint16 a_dems_get_elev_by_coord ( const VikCoord *coord, VikDemInterpol method);
 
 #endif
 #include <glib.h>

diff --git a/src/viktrack.c b/src/viktrack.c
index a5aa269e7238a79e097c52213c76fe290c0278ce..e2f455e4fe3401415bf427fd7f0e515ae30a03f2 100644 (file)
--- a/src/viktrack.c
+++ b/src/viktrack.c
@@ -809,7 +809,7 @@ void vik_track_apply_dem_data ( VikTrack *tr )
     /* TODO: of the 4 possible choices we have for choosing an elevation
      * (trackpoint in between samples), choose the one with the least elevation change
      * as the last */
-    elev = a_dems_get_elev_by_coord ( &(VIK_TRACKPOINT(tp_iter->data)->coord) );
+    elev = a_dems_get_elev_by_coord ( &(VIK_TRACKPOINT(tp_iter->data)->coord), VIK_DEM_INTERPOL_BEST );
     if ( elev != VIK_DEM_INVALID_ELEVATION )
       VIK_TRACKPOINT(tp_iter->data)->altitude = elev;
     tp_iter = tp_iter->next;

diff --git a/src/viktrwlayer_propwin.c b/src/viktrwlayer_propwin.c
index 769a1c60ae7a8f17fd7c6782f7d1cbc1a1da24d2..8e0225c0d31c84f3a5167f6ade02a09ac2e7bf2d 100644 (file)
--- a/src/viktrwlayer_propwin.c
+++ b/src/viktrwlayer_propwin.c
@@ -169,7 +169,7 @@ static void draw_dem_alt_speed_dist(VikTrack *tr, GdkDrawable *pix, GdkGC *alt_g
 
   for (iter = tr->trackpoints->next; iter; iter = iter->next) {
     int x, y_alt, y_speed;
-    gint16 elev = a_dems_get_elev_by_coord(&(VIK_TRACKPOINT(iter->data)->coord));
+    gint16 elev = a_dems_get_elev_by_coord(&(VIK_TRACKPOINT(iter->data)->coord), VIK_DEM_INTERPOL_SIMPLE);
     elev -= alt_offset;
     dist += vik_coord_diff ( &(VIK_TRACKPOINT(iter->data)->coord),
       &(VIK_TRACKPOINT(iter->prev->data)->coord) );

diff --git a/src/vikwindow.c b/src/vikwindow.c
index 28129741671f4f5a0e371c9be514ffb8fd8d5f37..ca1db1db94bcaf7e566fb5333174d4fbb40f6ed4 100644 (file)
--- a/src/vikwindow.c
+++ b/src/vikwindow.c
@@ -474,7 +474,8 @@ static void draw_mouse_motion (VikWindow *vw, GdkEventMotion *event)
   vik_viewport_screen_to_coord ( vw->viking_vvp, event->x, event->y, &coord );
   vik_coord_to_utm ( &coord, &utm );
   a_coords_utm_to_latlon ( &utm, &ll );
-  if ((alt = a_dems_get_elev_by_coord(&coord)) != VIK_DEM_INVALID_ELEVATION)
+  /* TODO: Change interpolate method according to scale */
+  if ((alt = a_dems_get_elev_by_coord(&coord, VIK_DEM_INTERPOL_SIMPLE)) != VIK_DEM_INVALID_ELEVATION)
     g_snprintf ( pointer_buf, 36, "Cursor: %f %f %dm", ll.lat, ll.lon, alt );
   else
     g_snprintf ( pointer_buf, 36, "Cursor: %f %f", ll.lat, ll.lon );