Optimising the calculation of unique edges of a Voronoi diagram in Matlab

Question

I am performing some calculations where I need to evaluate fluxes between Voronoi Polygons that are centred around some nodes. To do this I need to find the common edges between pairs of polygons eg. V1 & V2 as shown in the figure below. Each edge should be evaluated once only.

To do this I take my x- & y-coordinates of the nodes and do a Delaunay triangulation to find the neighbour nodes. I then run a loop to figure out which nodes have common vertices. I then calculate Voronoi polygons and create an array (istr) with the index of the 'edge' number and the value is the central voronoi polygon. The 'neigh' array then has the index as the 'edge' number and all the neighbouring polygon values. After a check to make sure I don't repeat this evaluation for each edge I then calculate the edges (i.e. the vertices shared between each polygon).

I can calculate the edges using the code below, however there is a massive bottleneck in the for...loop for calculating nodneigh as components of a cell array need to be accessed iteratively. What takes up even more time is the calculation of edge using a cell function to access the outputs of the Delaunay triangulation/Voronoi polygons.

My question is how can I speed up these two bottlenecks? Whilst I appreciate the flexibility of the cell array in Matlab, I feel like it is really slowing everything down when I don't need it to. I have tried padding out the cell array with NaNs, converting it to a matrix and performing a row-wise intersect but this hasn't been so successful: arrayfun takes much, much longer and I can't seem to use intersect with GPU computing.

% 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

You can loop through the edges and compute the distance from midpoint of the edge to all sites. Then sort the distances in ascending order and for inner voronoi polygons pick the first and the second. For outer polygons pick the first. Basically an edge separate/divide 2 polygons.

It should be faster (bottleneck #2). You don't need to compute (all) the neighbors of the triangulation for the voronoi diagram. It works when you walk all triangles edges and check for ghosts (double edges) (bottleneck #1).

Answer 2

Ok, finally got something I am happy with. By utilising the edges method of the 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 in matlab.

 % 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

These changes have made my code run 3-4 times faster than my original code. If anyone has any better ideas (perhaps something GPU-based would help?) I am open to nicer answers.