Rewrite spath.shortest_path to use new mcp classes

spath.shortest_path presents the same interface, but now uses the new mcp
classes internally, which means that it is now fully n-d, and can also find
the shortest path along any given axis. The old output format has been
retained, but is a bit weird for non-2d arrays. A more sensible n-d output
format is available, but not enabled by default.
This commit is contained in:
Zachary Pincus
2009-12-29 13:17:18 -05:00
parent 6f39426485
commit 8b2a8b5031
3 changed files with 130 additions and 29 deletions
+29
View File
@@ -0,0 +1,29 @@
import _mcp
cimport _mcp
cdef extern from "math.h":
double fabs(double f)
cdef class MCP_Diff(_mcp.MCP):
"""MCP_Diff(costs, offsets=None, fully_connected=True)
Find minimum-difference paths through an n-d costs array.
See the documentation for MCP for full details. This class differs from
MCP in that the cost of a path is not simply the sum of the costs along
that path.
This class instead assumes that the cost of moving from one point to
another is the absolute value of the difference in the costs between the
two points
"""
def __init__(self, costs, offsets=None, fully_connected=True):
"""__init__(costs, offsets=None, fully_connected=True)
See class documentation.
"""
_mcp.MCP.__init__(self, costs, offsets, fully_connected)
self.use_start_cost = 0
cdef _mcp.FLOAT_C _travel_cost(self, _mcp.FLOAT_C old_cost, _mcp.FLOAT_C new_cost, _mcp.FLOAT_C offset_length):
return fabs(old_cost - new_cost)
+73
View File
@@ -0,0 +1,73 @@
import numpy as np
import _spath
def shortest_path(arr, reach=1, axis=-1, output_indexlist=False):
"""Find the shortest path through an n-d array from one side to another.
Parameters
----------
arr : ndarray of float64
reach : int, optional
By default (``reach = 1``), the shortest path can only move
one row up or down for every step it moves forward (i.e.,
the path gradient is limited to 1). `reach` defines the
number of elements that can be skipped along each non-axis
dimension at each step.
axis : int, optional
The axis along which the path must always move forward (default -1)
output_indexlist: bool, optional
See return value `p` for explanation.
Returns
-------
p : iterable of int
For each step along `axis`, the coordinate of the shortest path.
If `output_indexlist` is True, then the path is returned as a list of
n-d tuples that index into `arr`. If False, then the path is returned
as an array listing the coordinates of the path along the non-axis
dimensions for each step along the axis dimension. That is,
`p.shape == (arr.shape[axis], arr.ndim-1)` except that p is squeezed
before returning so if `arr.ndim == 2`, then
`p.shape == (arr.shape[axis],)`
cost : float
Cost of path. This is the absolute sum of all the
differences along the path.
"""
# First: calculate the valid moves from any given position. Basically,
# always move +1 along the given axis, and then can move anywhere within
# a grid defined by the reach.
if axis < 0:
axis += arr.ndim
offset_ind_shape = (2*reach + 1,) * (arr.ndim - 1)
offset_indices = np.indices(offset_ind_shape) - reach
offset_indices = np.insert(offset_indices, axis, np.ones(offset_ind_shape), axis=0)
offset_size = np.multiply.reduce(offset_ind_shape)
offsets = np.reshape(offset_indices, (arr.ndim, offset_size), order='F').T
# Valid starting positions are anywhere on the hyperplane defined by
# position 0 on the given axis. Ending positions are anywhere on the
# hyperplane at position -1 along the same.
non_axis_shape = arr.shape[:axis] + arr.shape[axis+1:]
non_axis_indices = np.indices(non_axis_shape)
non_axis_size = np.multiply.reduce(non_axis_shape)
start_indices = np.insert(non_axis_indices, axis, np.zeros(non_axis_shape), axis=0)
starts = np.reshape(start_indices, (arr.ndim, non_axis_size), order='F').T
end_indices = np.insert(non_axis_indices, axis, -np.ones(non_axis_shape), axis=0)
ends = np.reshape(end_indices, (arr.ndim, non_axis_size), order='F').T
# Find the minimum-cost path to one of the end-points
m = _spath.MCP_Diff(arr, offsets=offsets)
costs, traceback = m.find_costs(starts, ends, find_all_ends=False)
# Figure out which end-point was found
for end in ends:
cost = costs[tuple(end)]
if cost != np.inf:
break
traceback = m.traceback(end)
if not output_indexlist:
traceback = np.array(traceback)
traceback = np.concatenate([traceback[:,:axis], traceback[:,axis+1:]], axis=1)
traceback = np.squeeze(traceback)
return traceback, cost
+28 -29
View File
@@ -1,34 +1,33 @@
import numpy as np
from numpy.testing import *
import scikits.image.graph.spath as spath
def test_basic():
x = np.array([[1, 1, 3],
[0, 2, 0],
[4, 3, 1]])
path, cost = spath.shortest_path(x)
assert_array_equal(path, [0, 0, 1])
assert_equal(cost, 1)
from scikits.image.graph import shortest_path
class TestShortestPath:
def test_basic(self):
x = np.array([[1, 1, 3],
[0, 2, 0],
[4, 3, 1]])
path, cost = shortest_path(x)
assert_array_equal(path, [0, 0, 1])
assert_equal(cost, 1)
def test_reach(self):
x = np.array([[1, 1, 3],
[0, 2, 0],
[4, 3, 1]])
path, cost = shortest_path(x, reach=2)
assert_array_equal(path, [0, 0, 2])
assert_equal(cost, 0)
def test_non_square(self):
x = np.array([[1, 1, 1, 1, 5, 5, 5],
[5, 0, 0, 5, 9, 1, 1],
[0, 5, 1, 0, 5, 5, 0],
[6, 1, 1, 5, 0, 0, 1]])
path, cost = shortest_path(x, reach=2)
assert_array_equal(path, [2, 1, 1, 2, 3, 3, 2])
assert_equal(cost, 0)
def test_reach():
x = np.array([[1, 1, 3],
[0, 2, 0],
[4, 3, 1]])
path, cost = spath.shortest_path(x, reach=2)
assert_array_equal(path, [0, 0, 2])
assert_equal(cost, 0)
def test_non_square():
x = np.array([[1, 1, 1, 1, 5, 5, 5],
[5, 0, 0, 5, 9, 1, 1],
[0, 5, 1, 0, 5, 5, 0],
[6, 1, 1, 5, 0, 0, 1]])
path, cost = spath.shortest_path(x, reach=2)
assert_array_equal(path, [2, 1, 1, 2, 3, 3, 2])
assert_equal(cost, 0)
if __name__ == "__main__":
run_module_suite()
run_module_suite()