Interpolation with Polynomials and Splines

• Download source files - 37.03 KB
• Download demo project - 7.52 KB

Introduction

First of all, take a look here, where Tobias von Petersdorff has published his interpolation form written in Java.

Have you tried to find some sample code written in C# that helps you to draw a Spline? I've tried with no success.

This is a port of his sample project.

switch (mode) { 
case(InterpolationMode.POLY):
// use divided difference algorithm to compute Newton form coefficients
for (int k=1; k<=np-1; k++) 
{
for (int i=0; i<=np-1-k; i++) 
{
yCoords[i]=(yCoords[i+1]-yCoords[i])/(xCoords[i+k]-xCoords[i]);
}
}

// for equidistant points along x-axis evaluate polynomial and draw line
float dt = ((float) this.graphPanel.Width-1)/npp;
for (int k=0; k<=npp; k++) 
{
x=k*dt;
// evaluate polynomial at t
y = yCoords[0];
for (int i=1; i<=np-1; i++) 
{ 
y=y*(x-xCoords[i])+yCoords[i];
}
// draw line

try
{
g.DrawLine(linePen, (int)oldx,(int)oldy,(int)x,(int)y);
}
catch(Exception ex)
{
Console.WriteLine(ex.Message + " x1={0} y1={1} x2={2} y2={3}", oldx, oldy, x, y);
}
oldx=x;
oldy=y;
}
break;
case(InterpolationMode.NAT_SPL): 
if (np>1)
{ 
float[] a = new float[np];
float x1;
float x2;
float[] h = new float[np];
for (int i=1; i<=np-1; i++)
{
h[i] = xCoords[i] - xCoords[i-1];
}
if (np>2)
{
float[] sub = new float[np-1];
float[] diag = new float[np-1];
float[] sup = new float[np-1];
for (int i=1; i<=np-2; i++)
{
diag[i] = (h[i] + h[i+1])/3;
sup[i] = h[i+1]/6;
sub[i] = h[i]/6;
a[i] = (yCoords[i+1]-yCoords[i])/h[i+1]-(yCoords[i]-yCoords[i-1])/h[i];
}
solveTridiag(sub,diag,sup,ref a,np-2);
}
// note that a[0]=a[np-1]=0
// draw
oldx=xCoords[0];
oldy=yCoords[0];
try
{
g.DrawLine(linePen, (int)oldx,(int)oldy,(int)oldx,(int)oldy);
}
catch(Exception ex)
{
Console.WriteLine(ex.Message + " x1={0} y1={1} x2={0} y2={1}", oldx, oldy);
}

for (int i=1; i<=np-1; i++) 
{ // loop over intervals between nodes
for (int j=1; j<=precision; j++)
{
x1 = (h[i]*j)/precision;
x2 = h[i] - x1;
y = ((-a[i-1]/6*(x2+h[i])*x1+yCoords[i-1])*x2 +
(-a[i]/6*(x1+h[i])*x2+yCoords[i])*x1)/h[i];
x=xCoords[i-1]+x1;
if ((lastPen == null) || (lastPen == linePen2))
linePen = linePen1;
else
linePen = linePen2;
lastPen = linePen;

try
{
g.DrawLine(linePen, (int)oldx,(int)oldy,(int)x,(int)y);
}
catch(Exception ex)
{
Console.WriteLine(ex.Message + " x1={0} y1={1} x2={2} y2={3}", oldx, oldy, x, y);
}

oldx=x;
oldy=y;
}
}
}
break;
}

History

• 20^th October, 2005: Initial post