如何基于给定邻接矩阵调整Voronoi图，实现区域1与6的邻接？

原始数据 #

城市质心坐标数据 #

df <- data.frame(
  ID = 1:8,
  x = c(4.0, 7.0, 2.5, 8.0, 3.0, 5.0, 6.0, 7.0),
  y = c(8, 9, 4, 5, 2, 2, 3, 1)
)
df
# 输出结果：
#   ID   x y
# 1  1 4.0 8
# 2  2 7.0 9
# 3  3 2.5 4
# 4  4 8.0 5
# 5  5 3.0 2
# 6  6 5.0 2
# 7  7 6.0 3
# 8  8 7.0 1

城市邻接矩阵 #

adj_matrix <- matrix(c(0, 1, 1, 1, 0, 1, 1, 0, 1, 0, 0, 0, 1, 1, 0, 0,
                       1, 1, 0, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 1, 0, 0,
                       1, 0, 1, 0, 1, 0, 1, 1, 1, 0, 0, 1, 0, 1, 0, 1,
                       0, 1, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0), nrow = 8, byrow = TRUE)
adj_matrix
# 输出结果：
#      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
# [1,]    0    1    1    1    0    1    1    0
# [2,]    1    0    0    0    1    1    0    0
# [3,]    1    1    0    0    0    0    1    1
# [4,]    0    0    1    0    0    1    0    0
# [5,]    1    0    1    0    1    0    1    1
# [6,]    1    0    0    1    0    1    0    1
# [7,]    0    1    0    1    1    0    1    0
# [8,]    0    0    0    0    0    0    0    0

原始Voronoi绘制脚本（R） #

library(sp)
library(deldir)

x=data.frame(x=df$x,y=df$y) 
crds=x 
center=colMeans(crds) 
bb=c(center[1]-5,center[1]+5,center[2]-5,center[2]+5) 
z=deldir(x=crds[,1],y=crds[,2],rw=bb) 
w=tile.list(z) 
polys=vector(mode='list',length=length(w)) 
for (i in seq(along=polys)){
  pcrds=cbind(w[[i]]$x,w[[i]]$y)
  pcrds=rbind(pcrds,pcrds[1,])
  polys[[i]]=Polygons(list(Polygon(pcrds)),ID=as.character(i)) 
} 
SP = SpatialPolygons(polys) 
voronoi = SpatialPolygonsDataFrame(
  SP, 
  data = data.frame(
    Regiao_ID = as.character(1:dim(df)[1]), 
    x = crds[,1], 
    y = crds[,2], 
    row.names = sapply(slot(SP, 'polygons'), function(x) slot(x, 'ID'))
  )
)

问题说明 #

使用上述脚本生成的Voronoi图中，区域1与6并未邻接，但根据给定的邻接矩阵，二者应当邻接。需要调整实现符合自定义邻接关系的Voronoi图。

解决方法 #

方法一：R语言调整 #

常规Voronoi图由空间距离决定邻接关系，要匹配自定义规则，推荐使用约束Delaunay三角剖分生成符合要求的Voronoi区域：

1. 安装并加载依赖包

install.packages(c("tripack", "sp", "spdep"))
library(tripack)
library(sp)
library(spdep)

2. 提取强制邻接边
   从邻接矩阵中提取需要强制邻接的点对（避免重复，仅保留上三角部分）：

# 获取邻接矩阵中非零的上三角点对
edges <- which(upper.tri(adj_matrix) & adj_matrix == 1, arr.ind = TRUE)
# 转换为约束边格式
constraints <- data.frame(p1 = edges[,1], p2 = edges[,2])

3. 生成约束Voronoi图

# 创建约束Delaunay三角剖分对象
ct <- tri.mesh(df$x, df$y, constraint = constraints)
# 提取Voronoi多边形
vor <- voronoi.mesh(ct)

# 转换为SpatialPolygonsDataFrame格式
polys <- vector("list", length(vor$nodes))
for (i in seq_along(polys)) {
  # 获取当前点的Voronoi顶点并闭合多边形
  verts <- vor$vertices[vor$adj[[i]], ]
  verts <- rbind(verts, verts[1, ])
  polys[[i]] <- Polygons(list(Polygon(verts)), ID = as.character(i))
}
SP <- SpatialPolygons(polys)
custom_voronoi <- SpatialPolygonsDataFrame(
  SP,
  data = df,
  row.names = sapply(slot(SP, "polygons"), function(x) slot(x, "ID"))
)

4. 验证邻接关系

# 检查区域1的邻接多边形
adjacent <- poly2nb(custom_voronoi)
adjacent[[1]]
# 输出应包含6，说明区域1和6已邻接

方法二：Python语言调整 #

在Python中可结合scipy、shapely和pysal实现约束Voronoi图：

1. 安装依赖库

pip install numpy scipy shapely pysal matplotlib

2. 代码实现

import numpy as np
import matplotlib.pyplot as plt
from shapely.geometry import Point, Polygon, LineString
from shapely.ops import split
from scipy.spatial import Voronoi

# 加载数据
coords = np.array([
    [4.0, 8], [7.0, 9], [2.5, 4], [8.0, 5],
    [3.0, 2], [5.0, 2], [6.0, 3], [7.0, 1]
])
adj_matrix = np.array([
    [0,1,1,1,0,1,1,0], [1,0,0,0,1,1,0,0],
    [1,1,0,0,0,0,1,1], [0,0,1,0,0,1,0,0],
    [1,0,1,0,1,0,1,1], [1,0,0,1,0,1,0,1],
    [0,1,0,1,1,0,1,0], [0,0,0,0,0,0,0,0]
])

# 生成基础Voronoi图
vor = Voronoi(coords)

# 针对区域1和6手动调整边界（0-based索引）
p1 = Point(coords[0])
p6 = Point(coords[5])
# 计算两点中垂线作为公共边界
midpoint = Point((p1.x + p6.x)/2, (p1.y + p6.y)/2)
dir_vec = np.array([p6.y - p1.y, p1.x - p6.x])
dir_vec = dir_vec / np.linalg.norm(dir_vec)
# 延伸中垂线
line_start = midpoint + dir_vec * 10
line_end = midpoint - dir_vec * 10
clip_line = LineString([(line_start.x, line_start.y), (line_end.x, line_end.y)])

# 获取原始Voronoi多边形并裁剪
poly1 = Polygon(vor.vertices[vor.regions[vor.point_region[0]]])
poly6 = Polygon(vor.vertices[vor.regions[vor.point_region[5]]])
# 裁剪后保留需要的部分
poly1_new = split(poly1, clip_line).geoms[0]
poly6_new = split(poly6, clip_line).geoms[1]

# 绘制调整后的Voronoi图
fig, ax = plt.subplots(figsize=(8,6))
# 绘制所有区域
for i, region in enumerate(vor.regions):
    if not -1 in region and len(region) > 0:
        if i == vor.point_region[0]:
            ax.fill(*zip(*poly1_new.exterior.coords), alpha=0.5, color="#ff9999")
        elif i == vor.point_region[5]:
            ax.fill(*zip(*poly6_new.exterior.coords), alpha=0.5, color="#9999ff")
        else:
            polygon = vor.vertices[region]
            ax.fill(*zip(*polygon), alpha=0.5)
# 绘制质心点
ax.scatter(coords[:,0], coords[:,1], c='red', zorder=10)
plt.title("调整后Voronoi图（区域1与6邻接）")
plt.show()

核心思路 #

常规Voronoi图的邻接关系由空间距离决定，要匹配自定义邻接矩阵，核心是强制指定点对之间存在Voronoi边，主要有两种实现路径：

1. 约束Delaunay三角剖分：在三角剖分阶段加入强制邻接的边，对应的Voronoi图会自然生成符合要求的邻接关系，这种方法更严谨。
2. 手动修改边界：针对不符合要求的区域，计算两点的中垂线，裁剪原有多边形使其沿中垂线邻接，适合快速调整局部区域。

内容的提问来源于stack exchange，提问作者César Macieira