在 Matlab 中优化 Voronoi 图的唯一边的计算
Optimising the calculation of unique edges of a Voronoi diagram in Matlab
我正在执行一些计算，我需要评估以某些节点为中心的 Voronoi 多边形之间的通量。为此，我需要找到多边形对之间的公共边，例如。 V1&V2如下图所示。每条边应仅评估一次。
为此，我获取节点的 x 和 y 坐标并执行 Delaunay triangulation 以查找相邻节点。然后我 运行 一个循环来找出哪些节点具有公共顶点。然后我计算 Voronoi 多边形并创建一个数组 (istr)，索引为 'edge' 数字，值是中央 voronoi 多边形。 'neigh' 数组然后将索引作为 'edge' 数字和所有相邻的多边形值。在检查以确保我不会对每条边重复此评估之后，我然后计算边（即每个多边形之间共享的顶点）。
我可以使用下面的代码计算边缘，但是在计算 nodneigh 的 for... 循环中存在巨大瓶颈，因为元胞数组的组件需要迭代访问。占用更多时间的是 edge 的计算，使用单元函数访问 Delaunay triangulation/Voronoi 多边形的输出。
我的问题是如何加速这两个瓶颈？虽然我很欣赏 Matlab 中元胞数组的灵活性，但我觉得当我不需要它时它真的会减慢一切。我尝试用 NaN 填充元胞数组，将其转换为矩阵并执行逐行相交，但这并不是那么成功：arrayfun 需要很多很多时间，而且我似乎无法使用与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);
您可以遍历边缘并计算从边缘中点到所有站点的距离。然后按升序对距离进行排序，对于内部 voronoi 多边形，选择第一个和第二个。对于外部多边形，选择第一个。基本上是一条边 separate/divide 2 个多边形。
它应该更快（瓶颈 #2）。您不需要为 voronoi 图计算（所有）三角剖分的邻居。当您遍历所有三角形边并检查重影（双边）（瓶颈 #1）时，它会起作用。
好的，终于得到了令我满意的东西。利用matlab中DelaunayTriangulation class I could resolve the first bottleneck. I think there is a slight difference in the calculations now, but I am guessing they are now a bit more rigorous. I have managed to speed up the second bottleneck by moving the cell array to a matrix and then looping through this. Instead of intersect I am using ismember, partly because I found this helpful post about how to improve performance with intersect and setdiff的edges方法。
% 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 的东西会有所帮助？）我愿意接受更好的答案。
我正在执行一些计算，我需要评估以某些节点为中心的 Voronoi 多边形之间的通量。为此，我需要找到多边形对之间的公共边，例如。 V1&V2如下图所示。每条边应仅评估一次。
为此，我获取节点的 x 和 y 坐标并执行 Delaunay triangulation 以查找相邻节点。然后我 运行 一个循环来找出哪些节点具有公共顶点。然后我计算 Voronoi 多边形并创建一个数组 (istr)，索引为 'edge' 数字，值是中央 voronoi 多边形。 'neigh' 数组然后将索引作为 'edge' 数字和所有相邻的多边形值。在检查以确保我不会对每条边重复此评估之后，我然后计算边（即每个多边形之间共享的顶点）。
我可以使用下面的代码计算边缘，但是在计算 nodneigh 的 for... 循环中存在巨大瓶颈，因为元胞数组的组件需要迭代访问。占用更多时间的是 edge 的计算，使用单元函数访问 Delaunay triangulation/Voronoi 多边形的输出。
我的问题是如何加速这两个瓶颈？虽然我很欣赏 Matlab 中元胞数组的灵活性，但我觉得当我不需要它时它真的会减慢一切。我尝试用 NaN 填充元胞数组，将其转换为矩阵并执行逐行相交，但这并不是那么成功：arrayfun 需要很多很多时间，而且我似乎无法使用与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);
您可以遍历边缘并计算从边缘中点到所有站点的距离。然后按升序对距离进行排序，对于内部 voronoi 多边形，选择第一个和第二个。对于外部多边形，选择第一个。基本上是一条边 separate/divide 2 个多边形。
它应该更快（瓶颈 #2）。您不需要为 voronoi 图计算（所有）三角剖分的邻居。当您遍历所有三角形边并检查重影（双边）（瓶颈 #1）时，它会起作用。
好的，终于得到了令我满意的东西。利用matlab中DelaunayTriangulation class I could resolve the first bottleneck. I think there is a slight difference in the calculations now, but I am guessing they are now a bit more rigorous. I have managed to speed up the second bottleneck by moving the cell array to a matrix and then looping through this. Instead of intersect I am using ismember, partly because I found this helpful post about how to improve performance with intersect and setdiff的edges方法。
% 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 的东西会有所帮助？）我愿意接受更好的答案。