+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 */
+ /* sw */
+ cols[0] = (gint) floor((east - dem->min_east) / dem->east_scale);
+ rows[0] = (gint) floor((north - dem->min_north) / dem->north_scale);
+ ll[0].lon = (dem->min_east + dem->east_scale*cols[0])/3600;
+ ll[0].lat = (dem->min_north + dem->north_scale*rows[0])/3600;
+ /* nw */
+ cols[1] = cols[0];
+ rows[1] = rows[0] + 1;
+ ll[1].lon = ll[0].lon;
+ ll[1].lat = ll[0].lat + (gdouble)dem->north_scale/3600;
+ /* ne */
+ cols[2] = cols[0] + 1;
+ 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;
+ /* se */
+ cols[3] = cols[0] + 1;
+ 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]);
+ }
+
+#if 0 /* debug */
+ for (i = 0; i < 4; i++)
+ fprintf(stderr, "%f:%f:%d:%d ", ll[i].lat, ll[i].lon, dists[i], elevs[i]);
+ fprintf(stderr, " north_scale=%f\n", dem->north_scale);
+#endif
+
+ 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 )
+{
+ int i;
+ gint16 elevs[4], dists[4];
+ gint16 max_dist;
+ gdouble t = 0.0;
+ gdouble b = 0.0;
+
+ if (!dem_get_ref_points_elev_dist(dem, east, north, elevs, dists))
+ return VIK_DEM_INVALID_ELEVATION;
+
+ max_dist = 0;
+ for (i = 0; i < 4; i++) {
+ if (dists[i] < 1) {
+ return(elevs[i]);
+ }
+ if (dists[i] > max_dist)
+ max_dist = dists[i];
+ }
+
+ gdouble tmp;
+#if 0 /* derived method by Franke & Nielson. Does not seem to work too well here */
+ for (i = 0; i < 4; i++) {
+ tmp = pow((1.0*(max_dist - dists[i])/max_dist*dists[i]), 2);
+ t += tmp*elevs[i];
+ b += tmp;
+ }
+#endif
+
+ for (i = 0; i < 4; i++) {
+ tmp = pow((1.0/dists[i]), 2);
+ t += tmp*elevs[i];
+ b += tmp;
+ }
+
+ // fprintf(stderr, "DEBUG: tmp=%f t=%f b=%f %f\n", tmp, t, b, t/b);
+
+ return(t/b);
+
+}
+