Coordinates of the shortest path between two points in a 2D grid

I have a square grid of a specific size, say 50 x 50. There are two points, a source, say (1,1), and a destination, say (26,35). Is there a function I can utilize to get a column vector with the coordinates of the shortest path that an agent might take going from the source to the destination. The agent can travel horizontally, vertically, and diagonally.

One of the ways I can see of doing this would be to create a graph from the grid and use shortestpath function. However, the node number would have to be related to the coordinate with a function, and ideally I would like to avoid going into that complication. Any help is appreciated.

[Matt J]

However, the node number would have to be related to the coordinate with a function.

I don't see why that would be true. Why not a look-up table instead?

[Tejas]

I'm not sure how I could describe the connections between each pair of points in the grid to make a graph of it. As far as I understand it, for each point, I'd have to check whether its neighboring points are within the grid (all permutations of x+1 and y+1), and if they are, assign the corresponding node as being connected. In your opinion, would this be easier to implement than trying to find the shortest path using grid coordinates?

Answer 1

[Matt J]

I'm not sure how I could describe the connections between each pair of points in the grid to make a graph of it.

Creating the graph for a given set of grid dimensions M,N is quite easy (example below).

M=10;N=10; %grid dimensions

[I,J]=ndgrid(1:M,1:N);
IJ=[I(:),J(:)];

D=pdist2(IJ,IJ);

A=(D>=0.001 & D<=sqrt(2)*1.001).*D; %adjacency matrix

G=graph(A);

[Matt J]

Incidentally, my spatialgraph2D class

https://www.mathworks.com/matlabcentral/fileexchange/73630-spatialgraph2d

can be useful for visualizing the graph generated by the code above.

obj=spatialgraph2D(G,IJ(:,1),IJ(:,2));
plot(obj)

[Matt J]

Here is a sample shortest path caluclation between two points P1 and P2:

P1=[3,1]; P2=[4,3];

nodenum=@(P) sub2ind([M,N],P(1),P(2));

idx=shortestpath(G,nodenum(P1),nodenum(P2));

path_points = [I(idx),J(idx)]

resulting in,

path_points =

     3     1
     3     2
     4     3

[Tejas]

This is very helpful. I imagined creating a matrix where each (x,y) location would correspond to the node number of that point, but sub2ind does exactly that. Thank you!

[Tejas]

Is there an efficient way to do this step? I have realized that the variable D is huge if I increase M and N beyond 100. I suppose we don't really need to find the distance of every point with every other point in the grid, but I'm not sure of another way to build the adjacency matrix.

[Image Analyst]

Did you even see my answer below? There is a function bwdistgeodesic(). Did you read Steve's blog I referred you to, or perhaps even try to use/adapt it?

[Matt J]

Is there an efficient way to do this step? I have realized that the variable D is huge if I increase M and N beyond 100.

There are other input syntaxes to pdist2 that let you calculate only specific distances. One way to modify the original answer would be,

M=100;N=100; %grid dimensions

tic;

    [I,J]=ndgrid(1:M,1:N);
    IJ=[I(:),J(:)];
    
    [W,K]=pdist2(IJ,IJ,'euclidean','Smallest',9);
    
    W(1,:)=[];
    
    s=repmat(K(1,:),8,1);
    t=K(2:end,:);
    
    idx = W>=0.001 & W<=sqrt(2)*1.001;
    W=W(idx); s=s(idx);t=t(idx);
    
    st=sort([s(:),t(:)],2);
    A=unique( [st,W] ,'rows');
    
    
    G=graph(A(:,1),A(:,2),A(:,3));

toc; % Elapsed time is 0.304730 seconds.

[Tejas]

@Image Analyst, Steve's blog is far too difficult for me to include in my code. I need to understand what exactly is happening and be able to explain it to someone who will be using the simulation as well. Since Matt's code already works for me and is something I understand, I'd be happier with a more efficient manipulation of it than a completely different route. I tried some things and failed, which is why I asked here.

[Matt J]

Steve's blog is far too difficult for me to include in my code. I need to understand what exactly is happening

Well, I'm not sure what the difference is between what Steve blog's proposes and what Walter proposes, but the idea of using bwdistgeodesic() should be conceptually quite simple. You are framing the path traversal as the movement between two points in a binary image while stepping only on the white pixels. The main challenge is setting up an appropriate image (I showed how to do that here) and also locating a subset of the pixels which correspond to your points.

Answer 2

[Walter Roberson]

For the case where each location has the same "cost" to move over, and not all locations are passible, then:

https://www.mathworks.com/help/images/ref/bwdistgeodesic.html

Create a binary grid with 1 for the grid points that can be travelled through, 0 for impassible.

Call bwdistgeodesic passing in the binary image and the row and the column of the destination.

Now for any given source point, examine the corresponding point in the transform. If it is inf then you cannot reach the destination.

Otherwise repeat: move in the direction that decreases the count the most; if there are multiple places with the same count then there are multiple shortest path and you can choose any of them. When you get to 0 you have arrived.

The same transform matrix can be used for multiple source points to the same destination.

If you have the opposite situation, one source point and multiple destinations to investigate, just reverse source and destination for the above algorithm.

If all locations are passible (within the grid) and all of them have the same cost, then there are other algorithms that are much more simple.

[Tejas]

This certainly seems helpful for the cases where the agent is essentially travelling on a 3D landscape, where the height of each point corresponds to the cost of reaching it, and some locations cannot be travelled through. For now though, my simulation is restricted to a plane with equal cost of each step and no obstacles. Thank you anyway!

[Matt J]

For now though, my simulation is restricted to a plane with equal cost of each step and no obstacles.

That's not entirely true. Because your path planning is constrained to purely horizontal, vertical, and diagonal steps, it is equivalent to maneuvering through an image like that below, where the black regions are obstacles and the start/end points are the red sub-lattice points. Bear in mind also that my solution, as currently code, treats diagonal paths as having sqrt(2) times higher cost than vertical or horizontal steps. You can adjust that, if you wish, of course.

N=11;
BW=eye(N)+fliplr(eye(N))>0; 
BW(1,:)=1; BW(:,1)=1;
BW(end,:)=[];BW(:,end)=[];
BW=repmat(BW,5,5);
BW(end+1,:)=1; BW(:,end+1)=1;

imagesc(BW); colormap(gray(256))

[I,J]=ndgrid(1:N-1:size(BW,1));
hold on; scatter(I(:),J(:),'ro','MarkerFaceColor','r','SizeData',60); hold off

[George Abrahams]

I'll expand upon @Walter Roberson's answer, as I suspect people are arriving at this post with that in mind.

If there are obstacles and distance is all that you care about, bwdistgeodesic is probably the most efficient solution.

Your pixels might be weighted depending on some preference variable, however. For example, you might want to find the shortest walking path between 2 points on a map while minimizing increases in elevation (going up hills). Or, you might want to follow the veins in a retinal scan by following the path of lowest intensity, like in the figure below.

In this case, you'll need an undirected graph with an algorithm like Dijkstra's, implemented by the built-in shortestpath.

As you have figured, you'll need to therefore find someway of converting the image into a graph, with edges between connected pixels being assigned the relevent weights. For this, you could use my bwgraph function on File Exchange, along with the NodeWeights name-value argument. I posted an example of that here.

Answer 3

[Image Analyst]

See Steve's 5-part blog for a thorough discussion.

http://blogs.mathworks.com/steve/2011/11/01/exploring-shortest-paths-part-1/