Tuesday, 16 October 2012

c# BiCubic Interpolation

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