From 81c40198875923b4c30a2b457d7038c293272e34 Mon Sep 17 00:00:00 2001 From: Stefan van der Walt Date: Sun, 22 Nov 2009 13:09:20 +0200 Subject: [PATCH] Add Zachy Pincus's trace_path routine for finding a low-cos path between a given start and end point. --- scikits/image/graph/__init__.py | 1 + scikits/image/graph/setup.py | 3 + scikits/image/graph/trace_path.pyx | 114 +++++++++++++++++++++++++++++ 3 files changed, 118 insertions(+) create mode 100644 scikits/image/graph/trace_path.pyx diff --git a/scikits/image/graph/__init__.py b/scikits/image/graph/__init__.py index ab753358..da648dff 100644 --- a/scikits/image/graph/__init__.py +++ b/scikits/image/graph/__init__.py @@ -1,5 +1,6 @@ try: from spath import shortest_path + from trace_path import trace_path except ImportError: print """*** The shortest path extension has not been compiled. Run diff --git a/scikits/image/graph/setup.py b/scikits/image/graph/setup.py index 7ff77e6f..6b9950d6 100644 --- a/scikits/image/graph/setup.py +++ b/scikits/image/graph/setup.py @@ -15,9 +15,12 @@ def configuration(parent_package='', top_path=None): # This function tries to create C files from the given .pyx files. If # it fails, we build the checked-in .c files. cython(['spath.pyx'], working_path=base_path) + cython(['trace_path.pyx'], working_path=base_path) config.add_extension('spath', sources=['spath.c'], include_dirs=[get_numpy_include_dirs()]) + config.add_extension('trace_path', sources=['trace_path.c'], + include_dirs=[get_numpy_include_dirs()]) return config diff --git a/scikits/image/graph/trace_path.pyx b/scikits/image/graph/trace_path.pyx new file mode 100644 index 00000000..dc02bee2 --- /dev/null +++ b/scikits/image/graph/trace_path.pyx @@ -0,0 +1,114 @@ +import numpy as numpy +cimport numpy as numpy +cimport cython + +@cython.boundscheck(False) +def trace_path(numpy.ndarray[numpy.float32_t, ndim=2] costs not None, start, ends, diagonal_steps=True): + """Find the lowest-cost path from the start point to each given end point. + + Inputs: 'costs' array; 'start' (x, y) pair; list of 'ends' (x, y) pairs, + and optional 'diagonal_steps' boolean flag. + + Costs are given by the input array: a move onto any given position in the + costs array adds that cost to the path. Paths may be constrained to + vertical and horizontal moves only by passing False for the diagonal_steps + parameter. The costs must be non-negative! + + The array of cumulative costs from the starting point, and a list of paths + from the start to each end point are returned. + + Paths are found by (more or less) breadth-first search outward from the + starting point: each time a lower-cost route to a given pixel is found, that + pixel is marked "active"; the neighbors of all active pixels are then + examined to see if their costs can be lowered as well. This continues until + no pixels are marked active. + """ + if costs.min() < 0: + raise ValueError("All costs must be non-negative.") + try: + a, b = start + except: + raise ValueError("The start point must be an (x, y) pair") + if not (0 <= a < costs.shape[0] and 0 <= b < costs.shape[1]): + raise ValueError("The start point must fall within the array") + for end in ends: + try: + a, b = end + except: + raise ValueError("All end points must be (x, y) pairs") + if not (0 <= a < costs.shape[0] and 0 <= b < costs.shape[1]): + raise ValueError("The end points must fall within the array") + + cdef numpy.ndarray[numpy.float32_t, ndim=2] cumulative_costs = numpy.empty_like(costs) + cumulative_costs.fill(numpy.inf) + cumulative_costs[start] = 0 + costs_shape = (costs.shape[0], costs.shape[1]) + cdef numpy.ndarray[numpy.uint8_t, ndim=2] active_nodes = numpy.zeros(costs_shape, dtype=numpy.uint8) + active_nodes[start] = 1 + cdef numpy.ndarray[numpy.uint8_t, ndim=2] parent_nodes = numpy.empty(costs_shape, dtype=numpy.uint8) + parent_nodes.fill(255) + cdef numpy.ndarray[numpy.int8_t, ndim=2] offsets + if diagonal_steps: + offsets = numpy.array([[-1, -1], [-1, 0], [-1, 1], [0, -1], [0, 1], [1, -1], [1, 0], [1, 1]], dtype=numpy.int8) + else: + offsets = numpy.array([[-1, 0], [0, -1], [0, 1], [1, 0]], dtype=numpy.int8) + + cdef Py_ssize_t x, y, ox, oy, xo, yo, i + cdef Py_ssize_t a_xmax, a_xmin, a_ymax, a_ymin, tmp_xmax, tmp_xmin, tmp_ymax, tmp_ymin + cdef unsigned int xmax, ymax, active, num_steps + xmax = costs.shape[0] - 1 + ymax = costs.shape[1] - 1 + num_steps = 0 + tmp_xmax = tmp_xmin = start[0] + tmp_ymax = tmp_ymin = start[1] + cdef float current_cost, current_cumulative_cost, cumulative_cost, new_cost + + while True: + active = 0 + # iterate over array + for x in range(0, xmax+1): + for y in range(0, ymax+1): + if active_nodes[x, y]: + active_nodes[x, y] = 0 + active = 1 + current_cumulative_cost = cumulative_costs[x, y] + # iterate over offsets + for i in range(8): + ox = offsets[i, 0] + oy = offsets[i, 1] + xo = x + ox + yo = y + oy + if xo < 0 or xo > xmax or yo < 0 or yo > ymax: + continue + current_cost = costs[xo, yo] + new_cost = current_cost + current_cumulative_cost + # if a cheaper path to a given point is found, activate that point + if cumulative_costs[xo, yo] > new_cost: + cumulative_costs[xo, yo] = new_cost + parent_nodes[xo, yo] = i + active_nodes[xo, yo] = 1 + + if not active: + break + + cdef unsigned int startx, starty + startx = start[0] + starty = start[1] + return_paths = [] + # Trace the paths from the endpoints to the start + for end in ends: + path = None + x = end[0] + y = end[1] + if cumulative_costs[x, y] != numpy.inf: + path = [(x, y)] + while not (x == startx and y == starty): + i = parent_nodes[x, y] + ox = offsets[i, 0] + oy = offsets[i, 1] + x -= ox + y -= oy + path.append((x, y)) + path.reverse() + return_paths.append(path) + return cumulative_costs, return_paths \ No newline at end of file