Feature Request: Considering 'Knight's move' in find_cost/route_through_array

[marhofmann]

I found a paper which shows the impact of considering the "Knight's move" in shortest path analysis. I'll try to summerize the findings regarding this issue as follows:

"...To convert a raster grid into a raster network, cells are most commonly connected to their neighbors according to a specified radius R (...). R = 0 denotes connecting cells to their orthogonal neighbors (rook’s move), R = 1 denotes connecting cells to their orthogonal and diagonal neighbors (queen’s move), and R = 2 denotes additionally connecting cells via knight’s moves."

and further

"... Goodchild (...) calculated the worst-case geometric elongation for a shortest path when traversing a raster network of uniform cost, and found the R = 0 network imparts a 41.4% elongation error, the R = 1 network imparts an 8.2% elongation error, and the R = 2 network imparts a 2.79% elongation error. Higher radius values could be used to further reduce the elongation error, but Huber and Church (...) found that the R = 2 network provides the best trade-off between accuracy and computational burden on real geographic data."

After deeper analysis coming to the conclusion:

"... The R = 2 paths do sometimes align with knight’s move directions, but overall display the least amount of visible geometric distortion due to having the fewest route alignment restrictions. While discretizing continuous space will always result in some level of geometric distortion, it is clear that R = 2 connectivity dramatically reduces geometric distortion with minimal additional complexity as compared to the R = 1 connectivity that most GIS software currently uses."

"... it was found that R = 2 networks provided more alternatives, with less orientation error, and that their solutions consistently had lower objective costs than the R = 1 network paths."

Source

I find this very useful and promising. Would it be possible to implement it to find even "cheaper" paths through a raster dataset with scikit-image?

[jni]

Thanks for the issue @marhofmann!

Typically, we don't implement state of the art, "promising" approaches in scikit-image, favouring rather well-established algorithms that have already had high adoption in the scientific community.

Having said that, route_through_array could use some improvement. In particular, we should make it easier to specify arbitrary, custom neighbourhoods. Actually, MCP, the class that route_through_array is based on, does allow this:

https://scikit-image.org/docs/dev/api/skimage.graph.html#skimage.graph.MCP

See the description of offsets.

You can actually find the route through array source code and see that it's very simple. So you can implement your own using MCP and offsets. Perhaps an example for the gallery for implementing this approach would be a welcome contribution! Would you like to try your hand at that?

In the meantime, I'll definitely keep this open as a reminder to improve our path-finding API.

[marhofmann]

Thank your very much for your answer @jni! So if i understood you correctly, adding the offsets for "Knight's moves" to the MCP constructor should do the trick? But a "Knight's moves" respects movements by 2 in any direction eg: [2,1], [-1,2] etc. (total length of sqrt(5))

From the documentation of the MCP_Geometric class which should be more suitable for this case:

On the other hand, a move from (1, 1) to (2, 2) is along the diagonal and is sqrt(2) in length. Half of this move is within the pixel (1, 1) and the other half in (2, 2), so the cost of this move is calculated as (sqrt(2)/2)*costs[1,1] + (sqrt(2)/2)*costs[2,2]. These calculations don't make a lot of sense with offsets of magnitude greater than 1. Use the sampling argument in order to deal with anisotropic data.

So when i try to pass an offset with a magnitude greater 2 i receive the following (even when passing sampling = (1,1))

ValueError: all offset components must be 0, 1, or -1

So as far as i have seen, in _mcp.pyx in line 748 & 749

if np.absolute(self.offsets).max() > 1: raise ValueError('all offset components must be 0, 1, or -1')

should be adaped to this usecase?!

EDIT: So after looking a bit deeper into the code I noticed, that it might not be so straight forward. The move of length = sqrt(5) should be split up into 4 pixels! So if (0, 0) is the current pixel and a Knight's move to (2,1) should be processed, the cost of the pixels (0,0) (1,0), (1,1), and (2,1) should be considered in a fraction of sqrt(5). I'll think about how the cost of each pixel should be considered in this case.

[marhofmann]

Since i have never done anything with cython before and building the scikit-image extensions didn't look so easy to me, i tried a workaround by subclassing the MCP_Flexible from pure python.

I wanted to overload travel_cost in a way, that "Knight's moves" can be considered. Sadly i couldn't access index and new_index of find_cost from pure python, which is necessary to get the cost of all fields on the way of a Knight's move.

Example:

• Current position= (0, 0)
• New position after Knight's move= (2, 1)
• Fields/pixels which need to be considered by calculating cost: (0,0), (1,0), (1,1), (2,1)

Simply passing index and new_index to _travel_cost and make it accessible for MCP_Flexible here (without any usage in all other classes) would solve my problem and additionally would make the MCP_Flexible class even more flexible. Other kinds of moves could be considered as well!