（算法高手请进）如何将已知的一些无规则三维数据(x,y,z)进行网格化? (300分)

石油大学有本《计算机地质制图》上面有网格化的源代码。是FORTRAN的。
有本《计算机图形学》上面有几种插值方法和曲面拟合的源代码，是PASCAL的。

to:wk_knife
大佬呀, 你说的这几本书我一本都没找到！麻烦你好人做到底就把那本书的地址贴出来给我啦
免得我象无头苍蝇般乱撞....

TO:zwx,
你可不可以把你的源码发一份给我呀！ [] 到时等卖出去赚到钱了咱哥俩就一齐去玩如何呀

{没有梦的人生是无味的...}

估计是没有电子书，SORRY
我是建议你去找找书。

to chji:
留下Email,我发一个网格化的源代码给你，是Turbo C 2.0的，改到Delphi不会太难的。
说明一下，我测试过，控制点较少的情况下，应用效果不好，有时甚至与实际情况相差
十万八千里！

TO:zwx,
我的邮箱: need@21cn.com 或 index@163.com

其实我已找到几个内插值的算法，但使用起来与你说的情况一样:
"控制点较少的情况下，应用效果不好，" 真不知WinSurf的网格化功能是如何做到的!

继续关注

呵，不错啊，居然有人响应
请问你只是需要作图还是需要计算？
我以前是用opengl做的
iie@163.net

哈哈，很高兴还有这么多同志，好像有许多是石油行业的吗？
to chji:
我觉得可以使用surfer在后台网格化，然后，用Noswing的程序绘制等值线
to NoSwing ：
你的程序可不可以发给我一份
To iie:
如果我两方面都要呢？
to all:如果友好结果，我也重酬！

我以前做过，不过两方面都很差：）
图形显示是用OpenGL，计算是用矩形积分。

准确地说，你要做的是TIN(不规则三角网)，黄杏元有一本书有多种TIN生成的算法，另外，
有好多的OpenGL和DX程序都有此功能。

to:NoSwing,
请给我一份等值线源代码好么？
谢谢。
zhangsn98@sohu.com

up

我怕我会喜欢你！（你可是有......）

高手，可以教我学Delphi吗？我知道你很厉害的

前几天翻出大学时图书馆漏还的书，叫《逼近论》（我居然还有这么一本书，我记忆中一直以为是一本讲时间序列的书）
，讲的就是如何科学的插值——科学一词有点多余，不过事实如此。

用超星找找吧，很可能有。

我现在对那Delphi还是半懂不懂，看书也是老样子，没有进步一样，怎么办？？？？？？？？

这个问题我曾经作过。当然最好的网格化方法的是克立金方法，
这是数学地质中一种著名的插值方法，其算法本身而言，Surf的帮助中有较详细的说明。
但是因为涉及到一些统计量的计算，因此这种方法比较难于实现。还有一种简单的方法叫
反比半径法，这种方法的原理是已知点对估计点的影响随其距离的增加而单调递减。这种方法
非常简单，我就是用这种方法做的，但是其效果不是很好，在数据比较少的情况下。Surfur
中也有这种算法计算网格数据。当然还有许多别的算法，但效果都不如克立金。
以前我曾经下载过一个Fortran语言编制的网格数据生成的动态链接库，
可惜他用文件进行数据交换，我现在一时找不到了。过一段时间我会重新编制这方面的程序
，我将使用FORTRAN商业算法集IMSL中现有的函数SURF/DSURF 函数，其说明如下：
SURF/DSURF (Single/Double precision)
Compute a smooth bivariate interpolant to scattered data that is locally a quintic
polynomial in two variables.
Usage
CALL SURF (NDATA, XYDATA, FDATA, NXOUT, NYOUT, XOUT, YOUT,
SUR, LDSUR)
Arguments
NDATA — Number of data points. (Input)
NDATA must be at least four.
XYDATA — A 2 by NDATA array containing the coordinates of the interpolation
points. (Input)
These points must be distinct. The x-coordinate of the I-th data point is stored in
XYDATA(1, I) and the y-coordinate of the I-th data point is stored in XYDATA(2,
I).
FDATA — Array of length NDATA containing the interpolation values. (Input)
FDATA(I) contains the value at (XYDATA(1, I), XYDATA(2, I)).
NXOUT — The number of elements in XOUT. (Input)
NYOUT — The number of elements in YOUT. (Input)
XOUT — Array of length NXOUT containing an increasing sequence of points.
(Input)
These points are the x-coordinates of a grid on which the interpolated surface is to
be evaluated.
YOUT — Array of length NYOUT containing an increasing sequence of points.
(Input)
These points are the y-coordinates of a grid on which the interpolated surface is to
be evaluated.
SUR — Matrix of size NXOUT by NYOUT. (Output)
This matrix contains the values of the surface on the XOUT by YOUT grid, i.e.
SUR(I, J) contains the interpolated value at (XOUT(I), YOUT(J)).
LDSUR — Leading dimension of SUR exactly as specified in the dimension
statement of the calling program. (Input)
LDSUR must be at least as large as NXOUT.
Comments
1. Automatic workspace usage is
SURF 31 * NDATA + NXOUT * NYOUT + 6 * NDATA units, or
DSURF 31 * NDATA+ NXOUT * NYOUT + 12 * NDATA units.
Workspace may be explicitly provided, if desired, by use of
S2RF/DS2RF. The reference is
CALL S2RF (NDATA, XYDATA, FDATA, NXOUT, NYOUT, XOUT,
YOUT, SUR, LDSUR, IWK, WK)
The additional arguments are as follows:
IWK — Work array of length 31 * NDATA + NXOUT * NYOUT.
WK — Work array of length 6 * NDATA.
2. Informational errors
Type Code
4 5 The data point values must be distinct.
4 6 The XOUT values must be strictly increasing.
4 7 The YOUT values must be strictly increasing.
3. This method of interpolation reproduces linear functions.
Algorithm
This routine is designed to compute a C  interpolant to scattered data in the
plane. Given the data points
x y f in i i i i
N , , 1 6 < A=1
SURF returns (in SUR, the user-specified grid) the values of the interpolant s. The
computation of s is as follows: First the Delaunay triangulation of the points
x y i i i
N , 1 6 < A=1
is computed. On each triangle T in this triangulation, s has the form
s x y c x y x y T
m n
mn
T m n , , 05= &amp;aring; " &amp;Icirc;
+ &amp;pound;5
Thus, s is a bivariate quintic polynomial on each triangle of the triangulation. In
addition, we have
s(xL, yL) = fL for i = 1, &amp;frac14;, N
and s is continuously differentiable across the boundaries of neighboring
triangles. These conditions do not exhaust the freedom implied by the above
representation. This additional freedom is exploited in an attempt to produce an
interpolant that is faithful to the global shape properties implied by the data. For
more information on this routine, we refer the reader to the article by Akima
(1978). The grid is specified by the two integer variables NXOUT, NYOUT that
represent, respectively, the number of grid points in the first (second) variable
and by two real vectors that represent, respectively, the first (second) coordinates
of the grid.
从说明来看这个算法好性不是克立金法。也不知道效果到底怎么样