diff --git a/scikits/image/graph/_spath.pyx b/scikits/image/graph/_spath.pyx new file mode 100644 index 00000000..d4624ef2 --- /dev/null +++ b/scikits/image/graph/_spath.pyx @@ -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) diff --git a/scikits/image/graph/spath.py b/scikits/image/graph/spath.py new file mode 100644 index 00000000..92bd012e --- /dev/null +++ b/scikits/image/graph/spath.py @@ -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 + diff --git a/scikits/image/graph/tests/test_spath.py b/scikits/image/graph/tests/test_spath.py index e126a424..b28d1f26 100644 --- a/scikits/image/graph/tests/test_spath.py +++ b/scikits/image/graph/tests/test_spath.py @@ -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() \ No newline at end of file