到目前为止我发现的：一种常用的方法，例如在有限元模拟中，是网格划分或三角剖分。
以下是一组点的“Delaunay三角剖分”的gnuplot实现的尝试。 https://en.wikipedia.org/wiki/Delaunay_triangulation
但是，我知道 gnuplot 并不是真正完成此类任务的工具。
所以，我可能不知道有更好的解决方案。我很想了解它们。
德劳内三角剖分：
下面的代码肯定不是最有效的三角测量方法，欢迎改进
程序（短版）：
- 通过增加 x 对 N 个数据点进行排序，如果 x 相同，则增加 y
- 找到前 m>=3 个不共线的点
- 从 m 到 N 的循环点
3.1) 查找与点 m+1 的连接不与任何当前船体段相交的所有船体点
3.2) 将这些船体点连接到点 m+1 并相应地修改船体
- 循环所有内边缘
4.1) 找到包含当前边的 2 个三角形。这些形成一个四边形
4.2）如果四边形是凸的，检查对角线是否需要翻转（“Lawson-flip”）
4.3) 从 4) 重新开始，直到不再需要翻转为止
为了给三角形着色
- 以质心为第四个点将每个三角形分成 3 个四边形
- 根据各个数据点的 z 值为 3 个子四边形着色
评论：
- gnuplot 没有原生排序功能（尤其是按 >=2 列排序），因此您必须使用
sort（Linux 上已经包含，在 Windows 上您必须安装，例如来自 GnuWin 的 CoreUtils。
- 翻转边缘需要一些时间。我猜，它与
O(n^2) 缩放。因此，超过 100 个数据点就变得不切实际，因为它只会花费太长时间。但似乎有些算法应该在 O(n log n) 中运行。
- 欢迎改进，或者甚至在 gnuplot 中的实现会很棒;-)
代码：
### Delaunay triangulation (gnuplot implementation attempt, requires gnuplot 5.4)
reset session
# get some test data
Random=0 # set to 0 to read data from existing FILE
if (Random) {
FILE = "tbTriangulationRandom.dat"
set print FILE
do for [i=0:99] {
print sprintf("%g %g %g",x=invnorm(rand(0))*10,y=invnorm(rand(0))*10,x*y)
}
set print
}
else {
FILE = "tbTriangulationTest.dat"
}
# sort data by x increasing values and if x is identical by increasing y values
set table $Data
plot '<sort -n -k 1 -k 2 '.FILE u 1:2:3 w table
unset table
# definition of quite a few variables and functions
colX = 1 # x column
colY = 2 # y column
colZ = 3 # z column
N = |$Data| # number of points
EDGES = '' # list of (inner) edges
HULL = '' # list of hull segments
TRIANGLES = '' # list of triangles
HULLPOINTS = '' # list of hullpoints
array Px[N] # set point array size
array Py[N] # set point array size
array Pz[N] # set point array size
newEdge(p1,p2) = sprintf(" %d %d ",p1,p2)
Edge(n) = sprintf(" %s %s ",word(EDGES,2*n-1),word(EDGES,2*n))
EdgeP(n,p) = int(word(EDGES,2*n-2+p))
changeEdge(n,p1,p2) = (_edge=Edge(n), _pos = strstrt(EDGES,_edge), _pos ? \
EDGES[1:_pos-1].newEdge(p1,p2). \
EDGES[_pos+strlen(_edge):strlen(EDGES)] : EDGES)
TriangleCount(n) = words(TRIANGLES)/3
TriangleN(n) = sprintf(" %s %s %s ", \
word(TRIANGLES,3*n-2),word(TRIANGLES,3*n-1),word(TRIANGLES,3*n))
newTriangle(p1,p2,p3) = p1<p2 && p2<p3 ? sprintf(" %d %d %d ",p1,p2,p3) : \
p1<p3 && p3<p2 ? sprintf(" %d %d %d ",p1,p3,p2) : \
p2<p1 && p1<p3 ? sprintf(" %d %d %d ",p2,p1,p3) : \
p2<p3 && p3<p1 ? sprintf(" %d %d %d ",p2,p3,p1) : \
p3<p1 && p1<p2 ? sprintf(" %d %d %d ",p3,p1,p2) : \
sprintf(" %d %d %d ",p3,p2,p1)
changeTA(n,p1,p2,p3) = (TA=TriangleN(n), _pos = strstrt(TRIANGLES,TA), _pos ? \
TRIANGLES[1:_pos-1].newTriangle(p1,p2,p3). \
TRIANGLES[_pos+strlen(TA):strlen(TRIANGLES)] : TRIANGLES)
TAp(n,p) = int(word(TRIANGLES,3*n-3+p))
TAx(n,p) = Px[TAp(n,p)] # x-coordinate of point p of triangle n
TAy(n,p) = Py[TAp(n,p)] # y-coordinate of point p of triangle n
HullP(n,p) = int(word(HULL,2*n-2+p)) # hull segment point number
HScount(n) = int(words(HULL))/2 # number of hull segments
getHullPoints(n) = (_tmp = '', sum [_i=1:words(HULL)] ((_s=' '.word(HULL,_i).' ', _tmp = \
strstrt(_tmp,_s) ? _tmp : _tmp._s ), 0), _tmp)
removeFromHull(seg) = (seg, _pos = strstrt(HULL,seg), _pos ? \
HULL[1:_pos-1].HULL[_pos+strlen(seg):strlen(HULL)] : HULL)
# orientation of 3 points, either -1=clockwise, 0=collinear, 1=counterclockwise
Orientation(p1,p2,p3) = sgn((Px[p2]-Px[p1])*(Py[p3]-Py[p1]) - (Px[p3]-Px[p1])*(Py[p2]-Py[p1]))
# check for intersection of segment p1-p2 with segment p3-p4, 0=no intersection, 1=intersection
IntersectCheck(p1,p2,p3,p4) = (Orientation(p1,p3,p2)==Orientation(p1,p4,p2)) || \
(Orientation(p3,p1,p4)==Orientation(p3,p2,p4)) ? 0 : 1
Sinus(p1,p2) = (_dx=Px[p2]-Px[p1], _dy=Py[p2]-Py[p1], _dy/sqrt(_dx**2 + _dy**2))
### Macros for later use
# Creating inner edges datablock macro
CreateEdgeDatablock = '\
set print $EDGES; \
do for [i=1:words(EDGES)/2] { \
p1 = int(word(EDGES,2*i-1)); \
p2 = int(word(EDGES,2*i)); \
print sprintf("%g %g %g %g %d %d",Px[p1],Py[p1],Px[p2]-Px[p1],Py[p2]-Py[p1],p1,p2) \
}; \
set print '
# Creating hull datablock macro
CreateHullDatablock = '\
set print $HULL; \
do for [i=1:words(HULL)/2] { \
p1 = int(word(HULL,2*i-1)); \
p2 = int(word(HULL,2*i)); \
print sprintf("%g %g %g %g %d %d",Px[p1],Py[p1],Px[p2]-Px[p1],Py[p2]-Py[p1],p1,p2) \
}; \
set print '
# plotting everything
PlotEverything = '\
plot $EDGES u 1:2:3:4 w vec lw 1.0 lc "red" nohead, \
$HULL u 1:2:3:4 w vec lw 1.5 lc "blue" nohead, \
$Data u 1:2 w p pt 7 ps 0.5 lc "black", \
$Data u 1:2:($0+1) w labels offset 0.5,0.5 '
# put datpoints into arrays
set table $Dummy
plot $Data u (Px[int($0)+1]=column(colX),Py[int($0)+1]=column(colY),Pz[int($0)+1]=column(colZ),'') w table
unset table
# get first m>=3 points which are not all collinear
HULL = HULL.newEdge(1,2) # add first 2 points to hull in any case
do for [p=3:N] {
if (Orientation(p-2,p-1,p)==0) { # orientation==0 if collinear
HULL = HULL.newEdge(p-1,p)
}
else { break } # stop if first >=3 non-collinear points found
}
HPcountInit = words(getHullPoints(0)) # get initial number of hull points
# actual plotting starts from here
set offset 1,1,1,1
set key noautotitle
set palette rgb 33,13,10
set rmargin screen 0.8
plot $Data u 1:2 w p pt 7 ps 0.5 lc "black", \
'' u 1:2:($0+1) w labels offset 0.5,0.5
set label 1 at graph 0.02,0.97 "Adding points... "
# loop all data points
do for [p=HPcountInit+1:N] {
print sprintf("### Adding P%d",p)
HPlist = getHullPoints(0)
HPcount = words(HPlist)
set print $NewConnections # initalize/empty datablock for new connections
print ""
set print
do for [hpt in HPlist] { # loop and check all hull points
hp = int(hpt)
# print sprintf("Check hull point P%d", hp)
c = 0
do for [hs=1:HScount(0)] { # loop all hull segments
hp1 = HullP(hs,1)
hp2 = HullP(hs,2)
# print sprintf("Check %d-%d with %d-%d", hp1, hp2, hp, p)
if (p!=hp1 && p!=hp2 && hp!=hp1 && hp!=hp2) {
c = c || IntersectCheck(hp1,hp2,hp,p)
if (c) { break }
}
}
if (c==0) { # if no intersections with hull
set print $NewConnections append # add new connections to datablock
print sprintf("%g %g", hp, Sinus(p,hp))
set print
}
}
# sort datablock clockwise (a bit cumbersome in gnuplot)
set table $ConnectSorted
plot $NewConnections u 1:2:2 smooth zsort # requires gnuplot 5.4.0
set table $Dummy
plot Connect='' $ConnectSorted u (Connect=Connect.sprintf(" %d",$1),'') w table
unset table
# add new edges
Ccount = int(words(Connect))
do for [i=1:Ccount-1] {
p1 = int(word(Connect,i))
p2 = int(word(Connect,i+1))
TRIANGLES = TRIANGLES.sprintf(" %d %d %d ", p1<p2?p1:p2, p2<p1?p1:p2, p) # numbers in ascending order
if (i==1) { HULL=HULL.newEdge(p1,p) }
if (i>1 && i<Ccount) { EDGES = EDGES.newEdge(p1,p) }
if (i==Ccount-1) {
HULL = HULL.newEdge(p2,p)
}
if (p!=HPcountInit+1) { # remove hull segments, except initial ones
NewEdge = p1<p2 ? sprintf(" %d %d ",p1,p2) : sprintf(" %d %d ",p2,p1)
# print sprintf("Remove %s",NewEdge)
HULL = removeFromHull(NewEdge)
EDGES = EDGES.NewEdge
@CreateEdgeDatablock
@CreateHullDatablock
@PlotEverything
}
}
}
# flip diagonal of a quadrangle if Det(p1,p2,p3,p4) and Orientation(p1,p2,p3) have the same sign
#
m11(p1,p4) = Px[p1]-Px[p4]
m21(p2,p4) = Px[p2]-Px[p4]
m31(p3,p4) = Px[p3]-Px[p4]
m12(p1,p4) = Py[p1]-Py[p4]
m22(p2,p4) = Py[p2]-Py[p4]
m32(p3,p4) = Py[p3]-Py[p4]
m13(p1,p4) = (Px[p1]-Px[p4])**2 + (Py[p1]-Py[p4])**2
m23(p2,p4) = (Px[p2]-Px[p4])**2 + (Py[p2]-Py[p4])**2
m33(p3,p4) = (Px[p3]-Px[p4])**2 + (Py[p3]-Py[p4])**2
Det(p1,p2,p3,p4) = m11(p1,p4)*(m22(p2,p4)*m33(p3,p4) - m32(p3,p4)*m23(p2,p4)) + \
m12(p1,p4)*(m23(p2,p4)*m31(p3,p4) - m33(p3,p4)*m21(p2,p4)) + \
m13(p1,p4)*(m21(p2,p4)*m32(p3,p4) - m31(p3,p4)*m22(p2,p4))
# create triangle data
set print $Triangles
do for [i=1:TriangleCount(0)] {
print sprintf("%g %g",TAx(i,1),TAy(i,1))
print sprintf("%g %g",TAx(i,2),TAy(i,2))
print sprintf("%g %g",TAx(i,3),TAy(i,3))
print sprintf("%g %g",TAx(i,1),TAy(i,1))
print ""
}
unset print
@CreateEdgeDatablock
@CreateHullDatablock
@PlotEverything
set label 1 "Flipping diagonals... "
###
# loop edges and check if need to flip. If on edge was flipped, start over again.
flip = 0
flipCount = 0
flippedAtLeastOnce = 1
while (flippedAtLeastOnce) {
flippedAtLeastOnce=0
do for [e=1:words(EDGES)/2] { # loop all inner edges
# find the 2 triangles with this edge
p1 = EdgeP(e,1)
p2 = EdgeP(e,2)
found = 0
do for [t=1:TriangleCount(0)] { # loop all triangles
tap1 = TAp(t,1)
tap2 = TAp(t,2)
tap3 = TAp(t,3)
p = p1==tap1 ? p2==tap2 ? tap3 : p2==tap3 ? tap2 : 0 : p1==tap2 ? p2==tap3 ? tap1 : 0 : 0
# print sprintf("%d %d %d: %d",tap1,tap2,tap3,p)
if (p!=0) {
if (found==1) {
ta2=t; p4=p;
flip = sgn(Det(p1,p2,p3,p4))==Orientation(p1,p2,p3)
flippedAtLeastOnce = flippedAtLeastOnce || flip
if (flip) {
flipCount = flipCount+1
print sprintf("Flip % 3d: %d-%d with %d-%d",flipCount,p1,p2,p3,p4)
EDGES = changeEdge(e,p3,p4)
TRIANGLES = changeTA(ta1,p1,p3,p4)
TRIANGLES = changeTA(ta2,p2,p3,p4)
@CreateEdgeDatablock
@CreateHullDatablock
@PlotEverything
break # start over again
}
}
if (found==0) { ta1=t; p3=p; found=1}
}
}
}
}
###
# create quadrangles datablock
Center2x(p1,p2) = (Px[p1]+Px[p2])/2. # x-center of 2 points
Center2y(p1,p2) = (Py[p1]+Py[p2])/2. # y-center of 2 points
Center3x(p1,p2,p3) = (Px[p1]+Px[p2]+Px[p3])/3. # x-center between 3 points
Center3y(p1,p2,p3) = (Py[p1]+Py[p2]+Py[p3])/3. # x-center between 3 points
set print $QUADRANGLES
do for [i=1:TriangleCount(0)] {
do for [p=0:2] {
z = Pz[TAp(i,p+1)]
tap1 = TAp(i,p+1)
tap2 = TAp(i,(p+1)%3+1)
tap3 = TAp(i,(p+2)%3+1)
print sprintf("%g %g %g", Px[tap1], Py[tap1], z)
print sprintf("%g %g %g", Center2x(tap1,tap2), Center2y(tap1,tap2), z)
print sprintf("%g %g %g", Center3x(tap1,tap2,tap3), Center3y(tap1,tap2,tap3), z)
print sprintf("%g %g %g", Center2x(tap1,tap3), Center2y(tap1,tap3), z)
print sprintf("%g %g %g", Px[tap1], Py[tap1], z)
print ''
}
}
set print
set label 1 "Coloring areas."
plot $QUADRANGLES u 1:2:3 w filledcurves closed lc palette, \
$EDGES u 1:2:3:4 w vec lw 1.0 lc "grey" nohead, \
$HULL u 1:2:3:4 w vec lw 1.5 lc "blue" nohead, \
$Data u 1:2:3 w p pt 7 ps 0.5 lc palette
### end of code
结果：
动画：（实际上，在我的旧笔记本电脑上大约需要 3 分钟）