一个古老的未解答的问题,我觉得很有趣。当然,您始终可以手动绘制箭头并填充没有意义的区域,尤其是对于大量点而言。我不知道 gnuplot 中有这样的 Voronoi 功能。当使用不规则网格中的数据点对地图进行着色时,这可能特别有用。我知道 gnuplot 有一些插值功能(检查help dgrid3d
)。
过去,我尝试过Delaunay 三角测量的 gnuplot 实现 https://stackoverflow.com/q/68507660/7295599,但是,速度相当慢。显然,下面的 Voronoi 方法似乎要快得多。
这个程序似乎与运行O(n^2)
,即对于 100、200、400 个数据点,我的旧笔记本电脑分别需要大约 3、11 和 40 秒。
它能做什么:
- 确定给定点
p0
与所有其他点的垂直平分线并将它们存储为向量 (x,y,a=角度)$Vectors
.
- 此外,添加边界(
xmin,ymax,ymin,ymax
) 也作为向量。
- 确定到点距离最小的向量(边界除外)
p0
(这将在p1
).
- 计算该向量与所有其他向量(除了它本身和前一个向量)的所有交集。逆时针方向找到最近的交叉路口(将在
p2
).
- 继续直到你得到第一个
p1
again.
我的解释和脚本可能不太容易理解。
两者当然都可以改进。欢迎提出建议!
例如,也许使用数组可能比写入数据块更快。
评论:脚本中似乎还有一个小错误。有时,可能会发生某些区域(尤其是边界处)未正确填充的情况。
Script:(需要 gnuplot 5.2.0,2017 年 9 月,因为数据块的索引)
### plot Voronoi-diagram from given set of points
reset session
time0 = time(0.0)
FILE = "SO40883823.dat"
# create some random test data
set table $Data
set samples 100
plot '+' u (sprintf("%g %g %g",x0=rand(0), y0=rand(0), z0=(int(rand(0)*0xffffff)))) w table
# plot FILE u 1:2:3 w table
unset table
set size square
set key noautotitle
xmin = 0 # adjust the borders
xmax = 1
ymin = 0
ymax = 1
set xrange[xmin:xmax]
set yrange[ymin:ymax]
set tics out
set offset 0.01,0.01,0.01,0.01
set angle degrees
j = {0,1}
x(i) = real(word($Data[i],1))
y(i) = real(word($Data[i],2))
z(i) = real(word($Data[i],3))
getDmin(colX,colY,colA,colL) = (L=column(colL), \
column(colA)==column(colA) || column(0)<3 ? dmin!=dmin ? \
(idx1=column(0),x1=column(colX),y1=column(colY),a1=column(colA),dmin=L) : L<dmin ? \
(idx1=column(0),x1=column(colX),y1=column(colY),a1=column(colA),dmin=L) : 0 : NaN)
M(x,y,a) = ((x-x0)*sin(a0) - (y-y0)*cos(a0))/(sin(a)*cos(a0)-sin(a0)*cos(a))
xs(x,y,a) = M(x,y,a)==0 ? NaN : x + cos(a)*M(x,y,a)
ys(x,y,a) = M(x,y,a)==0 ? NaN : y + sin(a)*M(x,y,a)
# orientation of 3 2D-points p0,p1,p2: -1=clockwise, 0=linear, +1=counterclockwise
Orientation(p0,p1,p2) = sgn((real(p1)-real(p0))*(imag(p2)-imag(p0)) - \
(real(p2)-real(p0))*(imag(p1)-imag(p0)))
colX = 1
colY = 2
colA = 3
colL = 4
set print $Outlines
print "" # initialize datablock
set print
do for [i=1:|$Data|] {
set print $Vectors
print sprintf("%g %g %g %g %g", x(i), ymin, 0, y(i)-ymin, 0) # bottom border
print sprintf("%g %g %g %g %g", xmax, y(i), 90, xmax-x(i), 1) # right
print sprintf("%g %g %g %g %g", x(i), ymax, 180, ymax-y(i), 2) # top
print sprintf("%g %g %g %g %g", xmin, y(i), 270, x(i)-xmin, 3) # left
x00 = x(i)
y00 = y(i)
z00 = z(i)
p0 = x00 + j*y00
do for [k=1:|$Data|] {
x0 = (x00 + x(k))/2.
y0 = (y00 + y(k))/2.
dx = x(k)-x00
dy = y(k)-y00
a0 = dx==0 && dy==0 ? NaN : atan2(dy,dx)+90
dL = sqrt(dx**2 + dy**2)
print sprintf("%g %g %g %g %g", x0, y0, a0, dL, k+3)
}
set print
# get vector with minimal distance from p0
dmin = NaN
stats [*:*][*:*] $Vectors u (getDmin(colX,colY,colA,colL)) nooutput # get x1,y1,a1,idx1
p1 = x1 + j*y1
set print $Outlines append
print sprintf("# Outline for point %g: %g,%g", i, x00,y00)
set print
idxC = idx00 = idx1
k = 0
while (idxC!=idx00 || k<2) {
k=k+1
x0 = x1; y0=y1; a0=a1; idx0 = idx1; idx2=idx1;
dmin = p2m = NaN
stats [*:*][*:*] $Vectors u ($0==idx0 || $0==idxC ? NaN : \
sprintf("%g %g %g %g %g", xc=xs($1,$2,$3), \
yc=ys($1,$2,$3), (p2=xc+j*yc, ox=Orientation(p0,p1,p2)), \
(dL = (x0-xc)**2+(y0-yc)**2), \
(ox==1 ? dmin!=dmin ? \
(x1=xc, y1=yc, a1=$3, idx1=$0, p2m=p2, dmin=dL) : dL<dmin ? \
(x1=xc, y1=yc, a1=$3, idx1=$0, p2m=p2, dmin=dL) : dmin : dmin,$0))) nooutput
set print $Outlines append
print sprintf("%g %g %g % 2d",x1, y1, z00, idx0)
set print
p1 = p2m
idx2 = idx1
idxC = idx0
}
set print $Outlines append
print ""; print ""
set print
print sprintf("Point: % 3d",i)
}
plot for [i=0:*] $Outlines u 1:2:3 index i w filledcurves lc rgb var fs solid 0.5, \
$Data u 1:2 w p pt 7 ps 0.5 lc "black"
print sprintf("Elapsed time: %.3f sec", time(0.0)-time0)
### end of script
Result:
OP文件中的数据为SO40883823.dat
第三列有一些颜色。为了创建这个,请评论line 10
并取消注释line 11
在上面的脚本中。下图有一些额外的点编号(不在脚本中)。
0 0 0xff0000
100 0 0x00ff00
0 100 0x0000ff
100 100 0xff00ff
10 20 0xee82ee
25 60 0xffff00
35 80 0xffa544
40 30 0x00ffff
生成一些随机测试数据:100 分
...以及 400 分: