GRSHORTESTPATH Find a shortest path between two nodes in the graph.
PATH = grShortestPath(NODES, EDGES, NODE1, NODE2, WEIGHTS)
NODES and EDGES defines the graph, NODE1 and NODE2 are indices of the
node extremities, and WEIGHTS is the set of weights associated to each
edge.
The function returns a set of node indices.
PATH = grShortestPath(NODES, EDGES, NODEINDS, WEIGHTS)
Specifies two or more points that must be traversed by the path, in the
specified order.
% Create a simple tree graph, and compute shortest path
[x y] = meshgrid([10 20 30], [10 20 30]);
nodes = [x(:) y(:)];
edges = [1 2;2 3;2 5;3 6; 4 5;4 7;5 8; 8 9];
drawGraph(nodes, edges)
axis equal; axis([0 40 0 40]);
drawNodeLabels(nodes, 1:9)
path = grShortestPath(nodes, edges, 1, 9);
% same as:
path = grShortestPath(nodes, edges, [1 9]);
See also
graphs, grFindMaximalLengthPath0001 function [nodePath, edgePath] = grShortestPath(nodes, edges, ind0, ind1, edgeWeights) 0002 %GRSHORTESTPATH Find a shortest path between two nodes in the graph. 0003 % 0004 % PATH = grShortestPath(NODES, EDGES, NODE1, NODE2, WEIGHTS) 0005 % NODES and EDGES defines the graph, NODE1 and NODE2 are indices of the 0006 % node extremities, and WEIGHTS is the set of weights associated to each 0007 % edge. 0008 % The function returns a set of node indices. 0009 % 0010 % PATH = grShortestPath(NODES, EDGES, NODEINDS, WEIGHTS) 0011 % Specifies two or more points that must be traversed by the path, in the 0012 % specified order. 0013 % 0014 % % Create a simple tree graph, and compute shortest path 0015 % [x y] = meshgrid([10 20 30], [10 20 30]); 0016 % nodes = [x(:) y(:)]; 0017 % edges = [1 2;2 3;2 5;3 6; 4 5;4 7;5 8; 8 9]; 0018 % drawGraph(nodes, edges) 0019 % axis equal; axis([0 40 0 40]); 0020 % drawNodeLabels(nodes, 1:9) 0021 % path = grShortestPath(nodes, edges, 1, 9); 0022 % % same as: 0023 % path = grShortestPath(nodes, edges, [1 9]); 0024 % 0025 % See also 0026 % graphs, grFindMaximalLengthPath 0027 % 0028 0029 % ------ 0030 % Author: David Legland 0031 % E-mail: david.legland@inrae.fr 0032 % Created: 2011-05-22, using Matlab 7.9.0.529 (R2009b) 0033 % Copyright 2011-2024 INRA - Cepia Software Platform 0034 0035 % process the case where several node indices are specified in first arg. 0036 if length(ind0) > 1 0037 % following optional argument is edge values 0038 if exist('ind1', 'var') 0039 edgeWeights = ind1; 0040 else 0041 edgeWeights = ones(size(edges, 1), 1); 0042 end 0043 0044 % concatenate path pieces 0045 nodePath = ind0(1); 0046 edgePath = size(0, 2); 0047 for i = 2:length(ind0) 0048 [node0, edge0] = grShortestPath(nodes, edges, ind0(i-1), ind0(i), edgeWeights); 0049 nodePath = [nodePath ; node0(2:end)]; %#ok<AGROW> 0050 edgePath = [edgePath ; edge0]; %#ok<AGROW> 0051 end 0052 0053 return; 0054 end 0055 0056 % ensure weights are defined 0057 if ~exist('edgeWeights', 'var') 0058 edgeWeights = ones(size(edges, 1), 1); 0059 end 0060 0061 0062 % check indices limits 0063 nNodes = size(nodes, 1); 0064 if max(ind0) > nNodes 0065 error('Start index exceed number of nodes in the graph'); 0066 end 0067 if max(ind1) > nNodes 0068 error('End index exceed number of nodes in the graph'); 0069 end 0070 0071 0072 %% Initialisations 0073 0074 % consider infinite distance for all nodes 0075 dists = inf * ones(nNodes, 1); 0076 0077 % initial node is at distance 0 by definition 0078 dists(ind0) = 0; 0079 0080 % create an array of indices for the predessor of each node 0081 preds = zeros(nNodes, 1); 0082 preds(ind0) = 0; 0083 0084 % allocate memory to store the subgraph, which is a tree 0085 edges2 = zeros(nNodes-1, 2); 0086 0087 % create an array of node flags 0088 unprocessedNodeInds = 1:nNodes; 0089 0090 edgeCount = 0; 0091 0092 0093 %% Iterate until all nodes are flagged to 1 0094 0095 while ~isempty(unprocessedNodeInds) 0096 % choose unprocessed node with lowest distance 0097 [tmp, ind] = min(dists(unprocessedNodeInds)); %#ok<ASGLU> 0098 ind = unprocessedNodeInds(ind); 0099 ind = ind(1); 0100 0101 % mark current node as processed 0102 unprocessedNodeInds(unprocessedNodeInds == ind) = []; 0103 0104 % if the current node is the target, it is not necessary to continue... 0105 if ind == ind1 0106 break; 0107 end 0108 0109 % find the index list of edges incident to current node 0110 adjEdgeInds = grAdjacentEdges(edges, ind); 0111 0112 % select edges whose opposite node has not yet been processed 0113 inds2 = sum(ismember(edges(adjEdgeInds, :), unprocessedNodeInds), 2) > 0; 0114 adjEdgeInds = adjEdgeInds(inds2); 0115 0116 % iterate over incident edges to update adjacent nodes 0117 for iNeigh = 1:length(adjEdgeInds) 0118 % extract index of current adjacent node 0119 edge = edges(adjEdgeInds(iNeigh), :); 0120 adjNodeInd = edge(~ismember(edge, ind)); 0121 0122 newDist = dists(ind) + edgeWeights(adjEdgeInds(iNeigh)); 0123 % dists(adjNodeInd) = min(dists(adjNodeInd), newDist); 0124 if newDist < dists(adjNodeInd) 0125 dists(adjNodeInd) = newDist; 0126 preds(adjNodeInd) = ind; 0127 end 0128 end 0129 0130 % find the list of adjacent nodes 0131 adjNodeInds = unique(edges(adjEdgeInds, :)); 0132 adjNodeInds(adjNodeInds == ind) = []; 0133 0134 % choose the adjacent nodes with lowest distance, and add the 0135 % corresponding edges to the subgraph 0136 if ~isempty(adjNodeInds) 0137 minDist = min(dists(adjNodeInds)); 0138 closestNodeInds = adjNodeInds(dists(adjNodeInds) <= minDist); 0139 for iCloseNode = 1:length(closestNodeInds) 0140 edgeCount = edgeCount + 1; 0141 edges2(edgeCount, :) = [ind closestNodeInds(iCloseNode)]; 0142 end 0143 end 0144 end 0145 0146 0147 %% Path creation 0148 0149 % create the path: start from end index, and identify successive set of 0150 % neighbor edges and nodes 0151 0152 nodeInd = ind1; 0153 nodePath = nodeInd; 0154 edgePath = []; 0155 0156 % find predecessors until we reach ind0 node 0157 while nodeInd ~= ind0 0158 newNodeInd = preds(nodeInd); 0159 nodePath = [nodePath ; newNodeInd]; %#ok<AGROW> 0160 0161 % search the edge (both directions) in the list of edges 0162 e_tmp = [nodeInd newNodeInd]; 0163 [~,edgeInd] = ismember ([e_tmp; e_tmp(end:-1:1)], edges, 'rows'); 0164 edgeInd(edgeInd == 0) = []; % erase the one that isn't there 0165 edgePath = [edgePath ; edgeInd]; %#ok<AGROW> 0166 0167 nodeInd = newNodeInd; 0168 end 0169 0170 % reverse the path 0171 nodePath = nodePath(end:-1:1); 0172 edgePath = edgePath(end:-1:1); 0173