Google Groups no longer supports new Usenet posts or subscriptions. Historical content remains viewable.
Dismiss

3D surface curve fit and interpolation?

2,222 views
Skip to first unread message

Leland C. Scott

unread,
Oct 25, 2008, 12:33:29 AM10/25/08
to
I have a collection of data points (x,y,z) for a surface, not on a regular
grid, that I want to generate a surface approximation and then use it to
interpolate "z" values based on the "x" and "y" values supplied. Ideally the
data values would be in a 3 column wide matrix read in from disk at run
time. Each column vector would correspond to the "x", "y" and "z" data
values. The interpolation I want to do is going to be a function that
returns a value for "z" based on the supplied values for "x" and "y". I
would prefer cubic spline interpolation between the data points with linear
interpolation used at the edges of the surface. I've looked through the
scilab help file and several functions look like they may do what I want,
but I'm not sure which one, and the documentation isn't very extensive. I
can do this in MathCAD but I want to dump that package for something that is
open source like scilab.

Any advise? Also most useful would be some real code examples to look at,
scilab files, PDF files, on the web etc.

--
Regards,
Leland C. Scott
KC8LDO


Richard Owlett

unread,
Oct 25, 2008, 10:05:58 AM10/25/08
to

You might want to ask over on comp.dsp
This sounds like a signal processing type problem, though I'm no expert.
The term "non uniform sampling" comes to mind.
HTH


bruno

unread,
Oct 26, 2008, 4:46:51 AM10/26/08
to

You can try the cshep2d (build the interpolant)
and eval_cshep2d (evaluation of the interpolant
on (x,y) points) functions. The cshep2d help
page provides a short example.

hth
Bruno

Richard Owlett

unread,
Oct 26, 2008, 7:28:30 AM10/26/08
to

Could you point me to a web reference of the underlying math good for
someone who had 2 years of college math >40 years ago. My Google search
wasn't very productive - poor choice of search terms probably.
Thank you.

Enrico Segre

unread,
Oct 26, 2008, 10:12:33 AM10/26/08
to
I think I have already answered this question a few times on the ng in
the past, if you scan backward. Sparse point interpolation in N-d is
not a trival task, regardless of the form you store on file the input
data. Sophisticate routines usually first create a meshing of the
support points (for instance a Delaynay triangulation, but there are
other discussible choices), then represent the surface on each cell
according to some rules - for instance n-continuity - and weightings
of the neighbouring data, and define a prolongation method for points
out of the mesh. Scilab by itself just doesn't do any of this, besides
a primitive mesh2d triangulator, which is unusable beyond some tens of
points. AFAIK, what comes closest to it is either the function
grid2Ddata of my toolbox (somewhat primitive, inverse distance based).
CGAL (the full one) most likely also includes all the ingredients for
building something more sophisticate (if not a sparse interpolator
itself), but AFAIR the port CGLAB just doesn't address the problem,
and development there is apparently dead.
hth, Enrico

Enrico Segre

unread,
Oct 28, 2008, 3:25:27 AM10/28/08
to

right, I forgot - and that we were in touch at that time!
Enrico

pat....@laposte.net

unread,
Oct 28, 2008, 5:50:13 AM10/28/08
to

Possible solution I have used to deal with results from DOE. I need to
find a function Y=f(X1,X2,X3).

Hereafter is the code :

// Results to fit
X1=[2000,2000,4000,4000,4000,2000,2000,4000,3000,2000,4000,3000];
X2=[1,10,10,10,1,10,5.5,5.5,10,1,1,5.5];
X3=[1.5,1,1,1.5,1.5,1.5,1.5,1,1.5,1,1,1];
X=[X1;X2;X3];
Y=[199,191,174,171,176,189,190,174,177,204,181,181];
Z=[Y;X];

// Fitting function
function y=func_fit(x,p)
x1=x(1,:)
x2=x(2,:)
x3=x(3,:)

y=p(1)+p(2)*x1+p(3)*x2+p(4)*x3+p(5)*x1.*x2+p(6)*x2.*x3+p(7)*x3.*x1+p(8)*x1.*x1+p(9)*x2.*x2
endfunction

//The Error function
function e=Err(p,m),
e=(Y-func_fit(X,p)),
endfunction

//Solve the problem
p0=1:9;
[pfit,v]=lsqrsolve(p0',Err,size(X,2));

bruno

unread,
Oct 28, 2008, 7:29:11 AM10/28/08
to

> Could you point me to a web reference of the underlying math good for
> someone who had 2 years of college math >40 years ago. My Google search
> wasn't very productive - poor choice of search terms probably.

I don't know if there is any tuturial on the web (cubic shepard
interpolation in google seems to give only papers from math
journals). The underlying code is from R. J. Renka and described
in the paper:

R. J. Renka, "Algorithm 790. CSHEP2D: Cubic Shepard Method for
Bivariate Interpolation of Scattered Data", ACM Trans. Math.
Software, Vol. 25, No. 1, Mar. 1999, pp. 70-73.

In the source code there are a few explanations, see:

http://viewvc.scilab.org/bin/cgi/viewvc.cgi/trunk/scilab/modules/interpolation/src/fortran/cshep2d.f?revision=13046&view=markup

hth
Bruno

Richard Owlett

unread,
Oct 28, 2008, 8:26:12 AM10/28/08
to


Thank you. If I read carefully I'll get what I need.

Leland C. Scott

unread,
Nov 6, 2008, 12:23:49 AM11/6/08
to

"bruno" <bruno....@iecn.u-nancy.fr> wrote in message
news:d93e1032-270e-4277...@k30g2000hse.googlegroups.com...

I tried it and it does work, sort of anyway. I tried using values between
the contour curves on the graph and over some ranges the eval function
return value is way off. The MathCAD equivalent function has a way to limit
the area around the point being interpolated to improve accuracy. I don't
see anything like that for the cshep2d function in scilab.

bruno

unread,
Nov 6, 2008, 2:58:41 PM11/6/08
to
On 6 nov, 06:23, "Leland C. Scott" <kc8...@arrl.net> wrote:
> "bruno" <bruno.pin...@iecn.u-nancy.fr> wrote in message

I don't understand exactly your problem with cshep2d. I know that if
you evaluate the interpolant at some point (x,y) enough far away the
interpolation
points region, the interpolant is set to 0 (see the example at the
eval_cshep2d
help page). Did you mean you met this problem ?

hth
Bruno

Leland C. Scott

unread,
Nov 6, 2008, 7:48:29 PM11/6/08
to
Bruno;

No, the problem is the interpolation is just plain wrong. For example I have
two contour lines, one at a z = -1 and the other at z = -0.5. Picking a
point between the two should result in an interpolated value of
around -0.75, but instead I get a -0.1 something, clearly out of line with
what I expected. Other regions seem to work OK but over a region on the
graph I'm looking at the interpolated results just don't look right.

I've used MathCAD's multivariate polynomial function "loess()" with a small
"span" value works like I expected. See their help file blurb below.

____________________________________________________________________________
____________________________________________________________________________
Multivariate polynomial regression, localized (MathCAD Professional only)
loess(Mxy,vz,span) Returns a vector which interp uses to find a set of
second order polynomial surfaces that best fit the x, y and z data values
specified in Mxy and vz within a neighborhood of the x and y data values.
span controls the size of this neighborhood.

interp(vs,Mxy,vz,v) Returns the z value corresponding to the x and y values
you supply in v.

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


Arguments:
Mxy is an m x 2 matrix containing x and y coordinates.

vz is an m element vector containing the z coordinates corresponding to the
m points specified in Mxy.

vs is the vector generated by either loess or regress.

span is a positive real number for specifying how big a neighborhood you
want to use. Use larger values of span when the data behaves very
differently over different ranges of x. A good default value is span=0.75.

v is a two element vector whose elements are the x and y values at which
you want to evaluate the corresponding z coordinate on the regression
surface.

Notes:
a.. You can use loess when you have a set of measured z values
corresponding to x and y values and you want to fit a polynomial surface
through those z values.

b.. The loess function is restricted to at most four independent
variables.

c.. You should always use the interp function after using the loess
function.

d.. loess tries to fit different second order polynomials depending on
where you are on the curve. It does this by examining the data in a small
neighborhood of the point you are interested in. The span argument controls
the size of this neighborhood. As span gets larger, loess becomes equivalent
to regress with n=2

_____________________________________________________________________________

_____________________________________________________________________________

I was hoping that the scilab "cshep2d" function would work like the above,
but it doesn't seem to do very well.

The MathCAD reference manual cites the following as the basis for the
function:

"local polynomial estimation (Cleveland and Devlin, 1988)"

And I found the following on the web with more information:

http://en.wikipedia.org/wiki/Local_regression

--
Regards,
Leland C. Scott
KC8LDO

"bruno" <bruno....@iecn.u-nancy.fr> wrote in message
news:4eaec8ed-3d67-4d32...@w39g2000prb.googlegroups.com...

b.gif

poorn...@gmail.com

unread,
Jul 3, 2014, 11:33:36 PM7/3/14
to

Hi Leland,

I'm wondering if you could help me out?

I trying to create a 3D scatter plot of x,y,z data that has been sampled at non-uniform intervals.

Is there a way to do this Mathcad?

Thanks & regards,
Poornaka
0 new messages