representing terrain with bezier surfaces

[Derek Hecht]

Hello all,

At my job one of the things we do is determine elevation values for
map positions by examining sample elevation grids. I was wondering if
there is any source code/papers/discussions on taking a sample grid and
representing it with a bezier surface or some sort of quadratic surface
equation to speed up elevation determination. Instead of interpolating
elevation based on the grid, just plug in the x,y to the correct
equation and get the result. Has this been done at all, and is it fast
enough to be worth doing.

Any help/directions/pointers would be appreciated,

Thanks,

Derek

--
( BURP! )
@..@ / Derek Hecht, Software Engineer
(----) Hughes Training, Link Division
( >__< ) E-MAIL: dhe...@link.com
^^ ~~ ^^ PHONE : (817) 695-8722

[Daniel Phillips]

Derek Hecht says:
>
> At my job one of the things we do is determine elevation values for
>map positions by examining sample elevation grids. I was wondering if
>there is any source code/papers/discussions on taking a sample grid and
>representing it with a bezier surface or some sort of quadratic surface
>equation to speed up elevation determination. Instead of interpolating
>elevation based on the grid, just plug in the x,y to the correct
>equation and get the result. Has this been done at all, and is it fast
>enough to be worth doing.

It works fine. You use a bezier function, preferably cubic, with two
inputs (X, Z) and one output (Y). (I'm sure some will take issue
with my choice of coordinate names, but to each his own.) It's not
easy to get a high speed VR type rendering happening, though.
Good luck.

---------------
Daniel Phillips
phil...@dowcow.com

[Kekoa Proudfoot]

In article <322DC4...@link.com>, Derek Hecht <dhe...@link.com> wrote:
> At my job one of the things we do is determine elevation values for
>map positions by examining sample elevation grids. I was wondering if
>there is any source code/papers/discussions on taking a sample grid and
>representing it with a bezier surface or some sort of quadratic surface
>equation to speed up elevation determination. Instead of interpolating
>elevation based on the grid, just plug in the x,y to the correct
>equation and get the result. Has this been done at all, and is it fast
>enough to be worth doing.

It's not clear to me exactly what you want to do here (do you want to
replace quadrilateral regions of your map with Bezier patches?) ... in any
event, one of the people in the lab I'm working in wrote a Siggraph paper
which covers fitting B-splines to dense polygonal meshes in 3D. The
paper's web page is at:

http://www-graphics.stanford.edu/papers/surfacefitting/

Although the information in this paper is a bit more general than that
which I'm guessing you might need, the paper (and it's list of references)
might be a good jumping-off point for what you're trying to do.

As for how fast it is to get a z value from an x,y pair ... I can't say.

Hope this helps!

-Kekoa
ke...@graphics.stanford.edu

[Daniel Phillips]

Kekoa Proudfoot says:
>Derek Hecht wrote:
>>...I was wondering if there is any source code/papers/discussions

>>on taking a sample grid and representing it with a bezier surface
>>or some sort of quadratic surface equation to speed up elevation
>>determination. Instead of interpolating elevation based on the
>>grid, just plug in the x,y to the correct equation and get the
>>result. Has this been done at all, and is it fast enough to be
>>worth doing.
>

>...one of the people in the lab I'm working in wrote a Siggraph paper

>which covers fitting B-splines to dense polygonal meshes in 3D. The
>paper's web page is at:
>
>http://www-graphics.stanford.edu/papers/surfacefitting/

Right, I glossed over the detail of fitting the surface last time,
didn't I? Here is a good starting point:

Graphics Gems II: Least-Squares Approximation to Bezier Curves
and Surfaces
Doug Moore and Joe Warren

---------------
Daniel Phillips
phil...@dowcow.com

[Thatcher Ulrich]

Daniel Phillips wrote:

>
> Derek Hecht says:
> >
> > At my job one of the things we do is determine elevation values for

> >map positions by examining sample elevation grids. I was wondering if

> >there is any source code/papers/discussions on taking a sample grid and
> >representing it with a bezier surface or some sort of quadratic surface
> >equation to speed up elevation determination. Instead of interpolating
> >elevation based on the grid, just plug in the x,y to the correct
> >equation and get the result. Has this been done at all, and is it fast
> >enough to be worth doing.
>

> It works fine. You use a bezier function, preferably cubic, with two
> inputs (X, Z) and one output (Y). (I'm sure some will take issue
> with my choice of coordinate names, but to each his own.) It's not
> easy to get a high speed VR type rendering happening, though.
> Good luck.
>
> ---------------
> Daniel Phillips
> phil...@dowcow.com

Here's some working code. It refers to a couple of functions not listed here,
but you should be able to fill in the blanks easily. My dim memory is that
the code uses cubic patches, and ensures that the second derivative at grid
edges is continuous, and that the surface passes through the grid points.

// The following functions are for computing points on a splined surface
// between grid points. Our method uses bicubic interpolation splines to
// get a smooth surface which follows our grid of altitude points.

void mult_mm(frac result[4][4], frac a[4][4], frac b[4][4])
// 4x4 Matrix multiply.
{
int i, j, k;

for (i=0; i<4; i++) {
for (j=0; j<4; j++) {
result[i][j] = 0;
for (k=0; k<4; k++) {
result[i][j] += a[i][k]*b[k][j];
}
}
}
}

void mult_mv(frac result[4], frac a[4][4], frac b[4])
// 4x4 times 1x4...
{
int i, j;

for (i=0; i<4; i++) {
result[i] = 0;
for (j=0; j<4; j++) {
result[i] += a[i][j]*b[j];
}
}
}

frac mult_vv(frac a[4], frac b[4])
// vector multiply.
{
int i;
frac result;

result = 0;
for (i=0; i<4; i++) {
result += a[i]*b[i];
}
return(result);
}

frac mb[4][4]={ {-1, 3, -3, 1},
{3, -6, 3, 0},
{-3, 3, 0, 0},
{1, 0, 0, 0} };

frac mi[4][4]={ {-4.5, 13.5, -13.5, 4.5},
{9, -22.5, 18, -4.5},
{-5.5, 9, -4.5, 1},
{1, 0, 0, 0} };
frac mit[4][4]={{-4.5, 9, -5.5, 1},
{13.5, -22.5, 9, 0},
{-13.5, 18, -4.5, 0},
{4.5, -4.5, 1, 0}};

frac model_GetAltitude(const vector& ref)
// Returns the altitude at a given spot in the universe. Ignores the y element
// of the ref vector; just uses the x & z.
{
frac temp[4][4], ci[4][4], p[4][4];
frac x0, x1, y0, y1, z0, z1;
int iX, iZ;

x0 = ref.gelem(0)*fOneOverGridRes;
z0 = ref.gelem(2)*fOneOverGridRes;
iX = (int) floor(x0);
iZ = (int) floor(z0);

// Fill p matrix with grid point altitudes...
for (int i=0; i<4; i++) {
for (int j=0; j<4; j++) {
// Queries elevation value from grid.
p[i][j] = model_GetAltitude(iX+i-1, iZ+j-1);
}
}

// Make ci matrix for finding points on interpolated bicubic surface...
mult_mm(temp, p,mit);
mult_mm(ci, mi, temp);

frac vs[4] = {1,1,1,1}, vt[4]={1,1,1,1}, tempv[4];
frac s, t;

s = (x0-iX+1)/3; t = (z0-iZ+1)/3;
vs[2] = s; vs[1] = s*s; vs[0] = s*s*s;
vt[2] = t; vt[1] = t*t; vt[0] = t*t*t;

mult_mv(tempv, ci, vt);
y0 = mult_vv(vs, tempv);

return(y0);
}