在 Matlab 中优化 Voronoi 图的唯一边的计算

Question

我正在执行一些计算，我需要评估以某些节点为中心的 Voronoi 多边形之间的通量。为此，我需要找到多对多边形之间的共同边，例如。 V1 & V2 如下图所示。每条边只应评估一次。

为此，我获取节点的 x 和 y 坐标并执行 Delaunay triangulation 以查找相邻节点。然后我运行一个循环来找出哪些节点有共同的顶点。然后我计算 Voronoi 多边形并创建一个数组 (istr)，其索引为“边缘”编号，值为中心 voronoi 多边形。然后，“neigh”数组将索引作为“edge”编号和所有相邻多边形值。在检查以确保我不会对每条边重复此评估后，我然后计算边（即每个多边形之间共享的顶点）。

我可以使用下面的代码来计算边，但是在 for...循环中计算 nodneigh 存在很大的瓶颈，因为需要迭代地访问单元数组的组件。需要更多时间的是使用单元函数计算 edge 以访问 Delaunay 三角剖分/Voronoi 多边形的输出。

我的问题是如何加快这两个瓶颈？虽然我很欣赏 Matlab 中单元阵列的灵活性，但我觉得当我不需要它时，它确实会减慢一切。我尝试用 NaN 填充单元格数组，将其转换为矩阵并执行逐行相交，但这并没有那么成功：arrayfun 需要更长的时间，而且我似乎无法使用 intersect with GPU 计算。

% Create dummy data
nstr    = 1000;         % number of particles
x       = rand(nstr,1);   % particle x coordinates
y       = rand(nstr,1);   % particle y coordinates

% Delaunay triangulation
DT      = delaunayTriangulation(x,y);

% Determine node neighbors of the original nodes
nodneigh    = cell(nstr,1);
numtotneigh = 0;                        % initialise total # of neighbors
bla         = DT.vertexAttachments;     % Get the particle/triangle IDs

% BOTTLENECK 1: Find out which particles/triangles are neighbours
for istr = 1:nstr
    nodneigh{istr} = setdiff(unique(DT.ConnectivityList(bla{istr},:)),istr);
    numtotneigh = numtotneigh+length(nodneigh{istr});
end

% Construct Thiessen polygons by Voronoi tessalation
[voro_V,voro_R] = DT.voronoiDiagram;

% Bookkeeping - create an index of edges with associated voronoi regions
cellsz = cellfun(@size,nodneigh,'uni',false);
cellsz = cell2mat(cellsz);
cellsz = cellsz(:,1);
temp = [1:nstr];
idx([cumsum([1 cellsz(cellsz>0)'])]) = 1;

istr    = temp(cumsum(idx(1:find(idx,1,'last')-1)))'; % Region number
neigh   = vertcat(nodneigh{:});                       % Region neighbours 
neigh_m = mod(neigh,nstr);   

% Make sure neighbourship has not already been evaluated
idx             = neigh_m == 0;
neigh_m(idx,:)  = nstr;

neigh    = vertcat(nodneigh{:});  

% BOTTLENECK 2:
% Determine which edges are common to both central and neighbour regions
edge     = cellfun(@intersect,voro_R(istr),voro_R(neigh),...
                'UniformOutput',false);
edge     = cell2mat(edge);

Answer 1

您可以遍历边缘并计算从边缘中点到所有站点的距离。然后按升序对距离进行排序，对于内部 voronoi 多边形，选择第一个和第二个。对于外部多边形选择第一个。基本上是一条边分开/划分 2 个多边形。

它应该更快（瓶颈 #2）。您不需要为 voronoi 图计算（所有）三角剖分的邻居。当您遍历所有三角形边缘并检查重影（双边缘）（瓶颈＃1）时，它会起作用。

Answer 2

好的，终于得到了我满意的东西。通过利用DelaunayTriangulation 类的edges 方法，我可以解决第一个瓶颈。我认为现在的计算略有不同，但我猜它们现在更严格了。我已经设法通过将单元阵列移动到矩阵然后循环来加速第二个瓶颈。我使用ismember 而不是相交，部分原因是我发现this helpful post 关于如何在matlab 中使用intersect 和setdiff 提高性能。

 % Create dummy data
nstr    = 1000;         % number of particles
x       = rand(nstr,1);   % particle x coordinates
y       = rand(nstr,1);   % particle y coordinates

% Delaunay triangulation
DT      = delaunayTriangulation(x,y);

 % construct Thiessen polygons by Voronoi tessalation
[voro_V,voro_R] = DT.voronoiDiagram;
dt_ed = DT.edges;
istr = dt_ed(dt_ed(:,1)<=nstr,1);
neigh = dt_ed(dt_ed(:,1)<=nstr,2);

% Determine cross-sectional area and Dtrans of all nodes                  
neigh_m = mod(neigh,nstr);                    

% if neigh_m == 0, neigh_m = nstr; end
% if the index of the neighboring streamline is smaller than
% that of the current streamline, the relationship has already
% been evaluated
idx             = neigh_m == 0;
neigh_m(idx,:)  = nstr;

temp_is = nan(numel(istr),40);
temp_ne = nan(numel(istr),40);
edge = nan(numel(istr),2);
for index = 1:numel(istr)
   temp_len1     = length(voro_R{istr(index)});
   temp_len2     = length(voro_R{neigh(index)});
   temp_is(index,1:temp_len1) = voro_R{istr(index)};
   temp_ne(index,1:temp_len2) = voro_R{neigh(index)};
   edge(index,:) = temp_is(index,ismember(temp_is(index,:),temp_ne(index,:)));
end

这些更改使我的代码运行速度比原始代码快 3-4 倍。如果有人有更好的想法（也许基于 GPU 的东西会有所帮助？）我愿意接受更好的答案。