So today at work i had to implement a bicubic interpolation implementation to interpolate a geoid for height calculation, here is what i came up with, very noobish, For those interested you can download the example project https://www.dropbox.com/s/bk0cf1a32lms0ai/Interp.zip
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
namespace Interp
{
class Program
{
static Importer importer = new Importer();
static void Main(string[] args)
{
importer.load_tol("ausgeoid09.tol");
importer.load_tor("ausgeoid09.tor");
// E = X = Long = Col
// N = Y = Lat = Lines
double rett = get_height(132.32, -29.92); // 4.2
//rett = get_height(-8, 108);
//rett = get_height(-8.000001, 130.000001);
double hight = get_geoid(132.32,-29.92); // 77.935
}
static double get_height(double X, double Y)
{
int xx = 0;
int yy = 0;
xx = (int)((X - importer.UP_LEFT_E) / 0.01666667);
yy = (int)(( importer.UP_LEFT_N - Y) / 0.01666667);
return importer.tor_data[ yy * importer.COLUMNS + xx];
}
static double Interpolate(double z, double[] height)
{
return height[1] + 0.5 * z
* (height[2] - height[0] + z
* (2.0 * height[0] - 5.0 * height[1] + 4.0 * height[2] - height[3] + z
* (3.0 * (height[1] - height[2]) + height[3] - height[0])));
}
static double get_geoid(double X, double Y)
{
int xx = 0;
int yy = 0;
xx = (int)Math.Floor((X - importer.UP_LEFT_E) / 0.01666667) -1 ;
yy = (int)Math.Floor((importer.UP_LEFT_N - Y) / 0.01666667) -1 ;
double tX = X % 0.01666667;
double tY = Y % 0.01666667;
tX /= 0.01666667;
tY /= -0.01666667;
double[] four_line = new double[4];
double[] four_col = new double[4];
for (int m = 0; m < 4; m++)
{
for (int k = 0; k < 4; k++)
{
four_line[k] = importer.tor_data[(yy + m) * importer.COLUMNS + xx + k];
}
// calclulate the cubic interpolation
four_col[m] = Interpolate(tY, four_line);
}
double result1 = Interpolate(tX, four_col);
return result1;
}
}
}
No comments:
Post a Comment