I have a 2D Delaunay triangulation which contains a

scalar value at each vertex. I've created a terrain

map that uses bilinear interpolation.

/******************************************************

STRUCTURES : Related to my terrain map.

******************************************************/

typedef struct {

double x;

double y;

} Pnt2dd;

typedef Pnt2dd *Pnt2ddptr;

typedef struct {

double A;

double B;

double C;

} EqPlane;

typedef struct vrtx2 {

Pnt2dd point;

long z;

} Vrtx2dd;

typedef Vrtx2dd *Vrtx2ddptr;

typedef struct edge2 {

Tria2ddptr adjface[2];

Vrtx2ddptr endpts[2];

} Edge2dd;

typedef Edge2dd *Edge2ddptr;

typedef struct tria2 {

Edge2ddptr edge[3];

Vrtx2ddptr vertex[3];

EqPlane eq;

} Tria2dd;

typedef Tria2dd *Tria2ddptr;

/******************************************************

FUNCTION : void CalculateABC(Tria2ddptr face)

INPUTS : Pointer to a face of the "triangulation".

OUTPUT : Populated facet constants (A, B & C) for

member "eq" of the triangle structure.

PURPOSE : Populate the bilinear "eq" constants.

******************************************************/

void CalculateABC(Tria2ddptr face)

{

double x0 = face->vertex[0]->point.x;

double x1 = face->vertex[1]->point.x;

double x2 = face->vertex[2]->point.x;

double y0 = face->vertex[0]->point.y;

double y1 = face->vertex[1]->point.y;

double y2 = face->vertex[2]->point.y;

double z0 = face->vertex[0]->z;

double z1 = face->vertex[1]->z;

double z2 = face->vertex[2]->z;

double D, Da, Db, Dc;

D = x0*(y1-y2) - x1*(y0-y2) + x2*(y0-y1);

if (fabs(D) > 0.001)

{ /* Cramer's Rule (Aggies drool) */

Da = z0*(y1-y2) - z1*(y0-y2) + z2*(y0-y1);

Db = x0*(z1-z2) - x1*(z0-z2) + x2*(z0-z1);

Dc = x0*(y1*z2 - y2*z1) -

x1*(y0*z2 - y2*z0) +

x2*(y0*z1 - y1*z0);

face->eq.A = Da/D;

face->eq.B = Db/D;

face->eq.C = Dc/D;

}

else

/* Handle Erroneous Data */

}

For a point "p" inside a triangle I then determine a "z"

by substituting that triangle's stored constants within

this equation:

zp = face->eq.A * p.x +

face->eq.B * p.y +

face->eq.C;

---------------------------------------------------------

I'd like to enhance this by having it do some sort of

spline interpolation. I've looked in Numerical Recipe's

offering for bicubic splines, but it assumes a Cartesian

grid of data elements. My triangles are not only

irregularly shaped, but they vary radically in size.

Any suggestions about the following will be welcome:

1) How should I turn my "EqPlane" into an "EqSpline"?

2) Algorithm for a "CalculateSplineConsts" function?

3) Equation for actually perfroming the interpolation?

Let's not get *. (I'm a simple person with a simple

mind.)

Thanks,

Bob

Sent via Deja.com http://www.veryComputer.com/

Share what you know. Learn what you don't.