[SciPy-User] Dijkstra's algorithm on a lattice
Zachary Pincus
A bit off-topic, but before I write some C or cython to do this, I
thought I'd ask to see if anyone knows of existing code for the task
of finding the shortest (weighted) path between two points on a lattice.
Specifically, I have images with "start" and "end" pixels marked and I
want to find the path through the image with the lowest integrated
intensity. Trivial but tedious to implement, so if anyone has some
good tips I'd be happy to know. (There's already a left-to-right-
shortest-path-finder in the image scikit repository, but that's not
quite what I need.)
Thanks,
Zach
_______________________________________________
SciPy-User mailing list
SciPy...@scipy.org
http://mail.scipy.org/mailman/listinfo/scipy-user
C. Campbell
http://networkx.lanl.gov/
I hope this helps!
Colin
Windows 7: It works the way you want. Learn more.
Stéfan van der Walt
2009/11/19 Zachary Pincus <zachary...@yale.edu>:
> A bit off-topic, but before I write some C or cython to do this, I
> thought I'd ask to see if anyone knows of existing code for the task
> of finding the shortest (weighted) path between two points on a lattice.
This is what the shortest path routine in scikits.image is meant to
do. How can we modify it to make it more useful to you?
Cheers
Stéfan
Zachary Pincus
>> A bit off-topic, but before I write some C or cython to do this, I
>> thought I'd ask to see if anyone knows of existing code for the task
>> of finding the shortest (weighted) path between two points on a
>> lattice.
>
> This is what the shortest path routine in scikits.image is meant to
> do. How can we modify it to make it more useful to you?
Hi Stéfan,
Based on just a rudimentary perusal of that code, I thought it only
found the lowest-cost path from the left to the right of an array...
is this still the case? I'd been needing something to go from an
arbitrary point to any other arbitrary point.
I just started working on some cython code that will compute the
shortest path (under a given length) from any point to any other. It's
not quite Dijkstra (for which one needs to keep a sorted list of the
next pixels to visit), but more like breadth-first search to a given
depth. I'll be happy to send it over if that sort of thing sounds
useful.
Zach
Zachary Pincus
not bad for many purposes.
You pass in a 2D costs array, start- and end-points, and a maximum
number of iterations; the code then keeps track of the minimum
cumulative cost to each pixel from the starting point, as well as the
path thereto. It does this by keeping track of "active" pixels -- any
time a lower cumulative cost to a given pixel is found, that pixel is
made active. Each iteration, all the neighbors of the "active" pixels
are examined to see if their costs can be lowered too. Basically
breadth-first search.
Limitations and oddities:
- Currently, diagonal and vertical/horizontal steps are both allowed.
Easy enough to make this a parameter.
- Paths along the boundary aren't traced out because I didn't want to
deal with an if-check in the inner loop to make sure that the x,y
position plus the current offset wasn't out of bounds. This could be
addressed by (a) padding the input array by one pixel on each side,
(b) putting the if in the inner loop, or (c) having a second pass
through the edge pixels.
- In theory, the code could find the cheapest path from top-left to
bottom-right in a single pass because "active" pixels are marked
immediately as the code iterates through the array. So the max_iters
parameter doesn't guarantee that paths longer than that will not be
found. But it does guarantee that any path found less than that length
is optimal...
Let's say it's BSD licensed, in case anyone finds it of use.
Zach
Gary Ruben
I haven't looked at your code, but your description sounds like you've
got a very nice solution. When you originally asked this, I immediately
thought of Lee's algorithm, or Jarvis's distance-transform based path
planning, which uses a modified distance transform that fixes the start
and goal point costs. I didn't mention them because they don't cover to
your case, but I think your solution is a more general case or theirs -
i.e. you can use yours for navigation/maze solving by setting the
obstacle/wall values to something greater than the maximum distance to
the goal and the floor values to 0.
I think would be a very nice, general routine for scikits.image,
Gary R.
> ------------------------------------------------------------------------
Zachary Pincus
Attached is a simplified version that addresses all of the caveats I
had earlier, plus I added documentation and input-checking.
Boundaries are now handled properly (the bounds-checking is not too
slow, even on huge arrays), and the code now iterates until all paths
have been fully traced. Still BSD licensed; if it might be useful to
scikits.image, please feel free to include it.
Here's a simple "maze" solving example / test.
>>> import numpy
>>> import trace_path
>>> a = numpy.ones((8,8), dtype=numpy.float32)
>>> a[1:-1,1] = 0
>>> a[1,1:-1] = 0
>>> a
array([[ 1., 1., 1., 1., 1., 1., 1., 1.],
[ 1., 0., 0., 0., 0., 0., 0., 1.],
[ 1., 0., 1., 1., 1., 1., 1., 1.],
[ 1., 0., 1., 1., 1., 1., 1., 1.],
[ 1., 0., 1., 1., 1., 1., 1., 1.],
[ 1., 0., 1., 1., 1., 1., 1., 1.],
[ 1., 0., 1., 1., 1., 1., 1., 1.],
[ 1., 1., 1., 1., 1., 1., 1., 1.]], dtype=float32)
>>> trace_path.trace_path(a, (1, 6), [(7, 2)])
(array([[ 1., 1., 1., 1., 1., 1., 1., 1.],
[ 1., 0., 0., 0., 0., 0., 0., 1.],
[ 1., 0., 1., 1., 1., 1., 1., 1.],
[ 1., 0., 1., 2., 2., 2., 2., 2.],
[ 1., 0., 1., 2., 3., 3., 3., 3.],
[ 1., 0., 1., 2., 3., 4., 4., 4.],
[ 1., 0., 1., 2., 3., 4., 5., 5.],
[ 1., 1., 1., 2., 3., 4., 5., 6.]], dtype=float32),
[[(1, 6),
(1, 5),
(1, 4),
(1, 3),
(1, 2),
(2, 1),
(3, 1),
(4, 1),
(5, 1),
(6, 1),
(7, 2)]])
>>> trace_path.trace_path(a, (1, 6), [(7, 2)], diagonal_steps=False)
(array([[ 2., 1., 1., 1., 1., 1., 1., 2.],
[ 1., 0., 0., 0., 0., 0., 0., 1.],
[ 1., 0., 1., 1., 1., 1., 1., 2.],
[ 1., 0., 1., 2., 2., 2., 2., 3.],
[ 1., 0., 1., 2., 3., 3., 3., 4.],
[ 1., 0., 1., 2., 3., 4., 4., 5.],
[ 1., 0., 1., 2., 3., 4., 5., 6.],
[ 2., 1., 2., 3., 4., 5., 6., 7.]], dtype=float32),
[[(1, 6),
(1, 5),
(1, 4),
(1, 3),
(1, 2),
(1, 1),
(2, 1),
(3, 1),
(4, 1),
(5, 1),
(6, 1),
(6, 2),
(7, 2)]])
Zach
Stéfan van der Walt
2009/11/22 Zachary Pincus <zachary...@yale.edu>:
> Boundaries are now handled properly (the bounds-checking is not too
> slow, even on huge arrays), and the code now iterates until all paths
> have been fully traced. Still BSD licensed; if it might be useful to
> scikits.image, please feel free to include it.
This code looks really handy, and I'd love to add it. Would you
consider putting your code in a branch on github?
Simply go to the following URL and click "fork":
http://github.com/stefanv/scikits.image
Add your changes, push back to github and click the button "merge
request", then I'll make sure it gets merged to the main branch.
Thanks!
Stéfan
Stéfan van der Walt
> This code looks really handy, and I'd love to add it. Would you
> consider putting your code in a branch on github?
Actually, don't worry -- I'll add it quickly. Thanks for the contribution!
Cheers
Almar Klein
I have an implementation of the Minimum Cost Path method too. It uses
a binary heap to store the
active pixels (I call it the front). Therefore it does not need to
iterate over the the whole image at
each iteration, but simply pops the pixel with the minimum cumulative
cost from the heap. This
significantly increase speed.
Because I am still changing stuff to my MCP implementation, and I use
it for my research, I am a bit
reluctant to make it publicly available now. I'd be happy to share the
binary heap implementation though.
Looking at the posted code, I think it is incorrect. Each iteration,
you should only check the neighbours
of the pixel that has the minimum cumulative costs. That's why the
binary heap is so important to get
it fast.
I short while ago I made a flash app with some examples. It also
contains the pseudo code, although
I use a slighly different terminology (cost=speed, cumulative cost=time):
http://dl.dropbox.com/u/1463853/mcpExamples.swf
Here's a snipet of how I use the binary heap:
=====
from heap cimport BinaryHeapWithCrossRef
front = BinaryHeapWithCrossRef(frozen_flat)
value, ii= front.Pop() # to get the pixel of the front with the
minimum cumulative cost
front.Push(cumcost_flat[ii], ii) # to insert or update a value
=====
I use flat arrays and scalar indices, so I need to store only one
reference per pixel. This also
makes the implementation work for 3D data (or even higher dimensions
if you wish).
frozen_flat is a flat array, the same size as the imput (cost or
speed) image, that keeps track
whether pixels are frozen (indicating they wont change). A pixel is
frozen right after it is popped
from the heap. I use the same array in the heap to be able to update
values if a pixel is already
in the heap.
I hope this helps a bit. If not, feel free to ask.
Almar
Zachary Pincus
The binary heap code looks extremely useful in general -- thanks for
making it available! Do you have any license you want it under? (BSD
seems preferable if this is to be incorporated into a scikit, e.g.)
It would be great if you would be interested in making your MCP code
available too, even just as a base for others to hack on a bit (rather
than as a finished contribution), but this is of course up to you.
Otherwise I'll probably try to throw together something similar using
the heap code.
> Looking at the posted code, I think it is incorrect. Each iteration,
> you should only check the neighbours of the pixel that has the
> minimum cumulative costs. That's why the binary heap is so important
> to get it fast.
Incorrect means that the code might give a wrong result: is this the
case? I *think* I had satisfied myself that the implementation (while
suboptimal because it does extra work -- a lot in some cases!) would
yield the correct path. (Note that the code doesn't terminate when the
"end" pixel is first assigned a cost, but when no costs are changing
anywhere. Basically, brute-force search instead of Dijkstra's
algorithm. Again, while a lot more than necessary to just find the
minimum cost to a single point, this condition should be sufficient to
ensure that the minimum cost to *every* point in the array has been
found, right? If my analysis is wrong, though, it wouldn't be the
first time!)
> I use flat arrays and scalar indices, so I need to store only one
> reference per pixel. This also makes the implementation work for 3D
> data (or even higher dimensions if you wish).
Do you have code to take the flat index and the shape of the original
array and return the indices to the neighboring pixels, or is there
some other trick with that too?
Anyhow, thanks for your suggestions and contribution! I look forward
to making use of the heap.
Best,
Zach
> <heap.pxd><heap.pyx>_______________________________________________
Almar Klein
> The binary heap code looks extremely useful in general -- thanks for
> making it available! Do you have any license you want it under? (BSD
> seems preferable if this is to be incorporated into a scikit, e.g.)
BSD's fine :)
> It would be great if you would be interested in making your MCP code
> available too, even just as a base for others to hack on a bit (rather
> than as a finished contribution), but this is of course up to you.
> Otherwise I'll probably try to throw together something similar using
> the heap code.
I'll send you my current implementation, you should be able to distill something
usefull from that. The problem is that 1) I make use of another module of
mine to deal with anisotropic data (because my data is anisotropic) 2) I use
the MCP method in a specific way, therefore I needed to make the
implementation more flexible. This makes it less easy to use for other people,
and thus less handy to include in any toolkit as-is.
>> Looking at the posted code, I think it is incorrect. Each iteration,
>> you should only check the neighbours of the pixel that has the
>> minimum cumulative costs. That's why the binary heap is so important
>> to get it fast.
>
> Incorrect means that the code might give a wrong result: is this the
> case? I *think* I had satisfied myself that the implementation (while
> suboptimal because it does extra work -- a lot in some cases!) would
> yield the correct path. (Note that the code doesn't terminate when the
> "end" pixel is first assigned a cost, but when no costs are changing
> anywhere. Basically, brute-force search instead of Dijkstra's
> algorithm. Again, while a lot more than necessary to just find the
> minimum cost to a single point, this condition should be sufficient to
> ensure that the minimum cost to *every* point in the array has been
> found, right? If my analysis is wrong, though, it wouldn't be the
> first time!)
I really mean wrong, sorry. You now select any pixel that is active (meaning
an arbitrary pixel in the front), and from it calculate the cumulative cost
for its neighbours. However, it might be that the cumulative cost of this pixel
is changed later. Therefore you must take the active pixel with the lowest
cumulative cost; so you know it won't be changed.
>> I use flat arrays and scalar indices, so I need to store only one
>> reference per pixel. This also makes the implementation work for 3D
>> data (or even higher dimensions if you wish).
>
>
> Do you have code to take the flat index and the shape of the original
> array and return the indices to the neighboring pixels, or is there
> some other trick with that too?
Yes, it's in the code I'll send you.
Cheers,
Almar
Zachary Pincus
> distill something
> usefull from that. The problem is that 1) I make use of another
> module of
> mine to deal with anisotropic data (because my data is anisotropic)
> 2) I use
> the MCP method in a specific way, therefore I needed to make the
> implementation more flexible. This makes it less easy to use for
> other people,
> and thus less handy to include in any toolkit as-is.
Thanks -- I'll see what I can distill!
Each time a pixel's cumulative cost decreases, I put it back into the
"active" set (i.e. the front), so then the neighbors get re-examined
the next iteration, etc. This should suffice, right? Or am I *still*
missing something? Not that it really matters, because the approach is
rather inefficient for anything except finding the minimum cost to
every single array element (and even then I'm not certain this is
better). But I am curious if I just conceived of the whole problem
wrong. In which case perhaps I'm not the guy you want implementing
this for the scikit!
Zach
Almar Klein
Ah, now I see. I'm sorry. Yes, your code should produce the correct result,
although it will probably evaluate a lot of pixels more than once :)
Almar
Stéfan van der Walt
2009/11/25 Almar Klein <almar...@gmail.com>:
> Ah, now I see. I'm sorry. Yes, your code should produce the correct result,
> although it will probably evaluate a lot of pixels more than once :)
Zach's code was integrated into scikits.image here:
http://github.com/stefanv/scikits.image/blob/master/scikits/image/graph/trace_path.pyx
Any patches or suggestions for improvements are appreciated.
Thanks
Stéfan