From 6f39426485bdcfc768656f042224470de519e401 Mon Sep 17 00:00:00 2001 From: Zachary Pincus Date: Tue, 29 Dec 2009 13:05:53 -0500 Subject: [PATCH] Added classes for minimum-cost-pathfinding through n-d arrays. mcp.py presents the interface and a simple usage example; the classes of interest are in _mcp.pyx. Basic tests are included. --- scikits/image/graph/_mcp.pxd | 29 ++ scikits/image/graph/_mcp.pyx | 557 ++++++++++++++++++++++++++ scikits/image/graph/mcp.py | 41 ++ scikits/image/graph/tests/test_mcp.py | 113 ++++++ 4 files changed, 740 insertions(+) create mode 100644 scikits/image/graph/_mcp.pxd create mode 100644 scikits/image/graph/_mcp.pyx create mode 100644 scikits/image/graph/mcp.py create mode 100644 scikits/image/graph/tests/test_mcp.py diff --git a/scikits/image/graph/_mcp.pxd b/scikits/image/graph/_mcp.pxd new file mode 100644 index 00000000..04b24a1a --- /dev/null +++ b/scikits/image/graph/_mcp.pxd @@ -0,0 +1,29 @@ +""" This is the definition file for mcp.pyx. +It contains the definitions of the mcp class, such that +other cython modules can "cimport mcp" and subclass it. +""" + +cimport numpy as np +cimport heap + +# determine datatypes for mcp +ctypedef np.float64_t FLOAT_T +ctypedef double FLOAT_C + +cdef class MCP: + cdef heap.FastUpdateBinaryHeap costs_heap + cdef object costs_shape + cdef int dim + cdef object flat_costs + cdef object flat_cumulative_costs + cdef object traceback_offsets + cdef object flat_edge_map + cdef readonly object offsets + cdef object flat_offsets + cdef object offset_lengths + cdef int dirty + cdef int use_start_cost + # if use_start_cost is true, the cost of the starting element is added to + # the cost of the path. Set to true by default in the base class... + + cdef FLOAT_C _travel_cost(self, FLOAT_C old_cost, FLOAT_C new_cost, FLOAT_C offset_length) diff --git a/scikits/image/graph/_mcp.pyx b/scikits/image/graph/_mcp.pyx new file mode 100644 index 00000000..83df06a9 --- /dev/null +++ b/scikits/image/graph/_mcp.pyx @@ -0,0 +1,557 @@ +"""Cython implementation of Dijkstra's minimum cost path algorithm, +for use with data on a n-dimensional lattice. + +Original author: Zachary Pincus +Inspired by code from Almar Klein + +License: BSD + +Copyright 2009 Zachary Pincus + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions +are met: + +1. Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. +2. Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the distribution. + +THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR +IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES +OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. +IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, +INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT +NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF +THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +""" + + +import cython +cimport numpy as np +import numpy as np +cimport heap +import heap + +FLOAT = np.float64 + +def _get_edge_map(shape): + """_get_edge_map(shape) + + Return an array with edge points/lines/planes/hyperplanes marked. + + Given a shape (of length n), return an edge_map array with a shape of + original_shape + (n,), where, for each dimension, edge_map[...,dim] will + have zeros at indices not along an edge in that dimension, -1s at indices + along the lower boundary, and +1s on the upper boundary. + + This allows one to, given an nd index, calculate not only if the index is + at the edge of the array, but if so, which edge(s) it lies along. + """ + d = len(shape) + edges = np.zeros(shape+(d,), order='F', dtype=np.int8) + for i in range(d): + slices = [slice(None)] * (d+1) + slices[d] = i + slices[i] = 0 + edges[tuple(slices)] = -1 + slices[i] = -1 + edges[tuple(slices)] = 1 + return edges + +def _offset_edge_map(shape, offsets): + """_offset_edge_map(shape, offsets) + Return an array with positions marked where offsets will step out of bounds. + + Given a shape (of length n) and a list of n-d offsets, return a shape + (n,) + sized edge_map, where, for each dimension edge_map[...,dim] has zeros at + indices at which none of the given offsets (in that dimension) will step + out of bounds. If the value is nonzero, it gives the largest offset (in + terms of absolute value) that will step out of bounds in that direction. + + An example will be explanatory: + >>> offsets = [[-2,0], [1,1], [0,2]] + >>> edge_map = _offset_edge_map((4,4), offsets) + >>> edge_map[...,0] + array([[-2, -2, -2, -2], + [-2, -2, -2, -2], + [ 0, 0, 0, 0], + [ 1, 1, 1, 1]], dtype=int8) + + >>> edge_map[...,1] + array([[0, 0, 2, 1], + [0, 0, 2, 1], + [0, 0, 2, 1], + [0, 0, 2, 1]], dtype=int8) + + """ + d = len(shape) + edges = np.zeros(shape+(d,), order='F', dtype=np.int8) + offsets = np.asarray(offsets) + for i in range(d): + slices = [slice(None)] * (d+1) + slices[d] = i + distinct_offsets = set(offsets[:,i]) + if 0 in distinct_offsets: + distinct_offsets.remove(0) + for offset in sorted(distinct_offsets, key=np.absolute, reverse=True): + # process offsets with larger absolute values first, so that smaller + # offsets will overwrite the correct region of the offsets array. + slice_stop = -offset + if offset > 0: + slice_stop -= 1 + slice_step = -np.sign(offset) + slices[i] = slice(None, slice_stop, slice_step) + edges[tuple(slices)] = offset + return edges + +def make_offsets(d, fully_connected): + """make_offsets(d, fully_connected) + + Make a list of offsets from a center point defining a n-dim neighborhood. + + Parameters + ---------- + d : int + dimension of the offsets to produce + fully_connected : bool + whether the neighborhood should be singly- of fully-connected + + Returns + ------- + offsets : list of tuples of length `d` + + Example + ------- + + The singly-connected 2-d neighborhood is four offsets: + + >>> make_offsets(2, False) + [(-1,0), (1,0), (0,-1), (0,1)] + + While the fully-connected 2-d neighborhood is the full cartesian product + of {-1, 0, 1} (less the origin (0,0)). + """ + if fully_connected: + mask = np.ones([3]*d, dtype=np.uint8) + mask[tuple([1]*d)] = 0 + else: + mask = np.zeros([3]*d, dtype=np.uint8) + for i in range(d): + indices = [1]*d + indices[i] = (0, -1) + mask[tuple(indices)] = 1 + offsets = [] + for indices, value in np.ndenumerate(mask): + if value == 1: + indices = np.array(indices) - 1 + offsets.append(indices) + return offsets + +def _unravel_index_fortran(flat_indices, shape): + """_unravel_index_fortran(flat_indices, shape) + + Given a flat index into an n-d fortran-strided array, return an index tuple. + """ + strides = np.multiply.accumulate([1] + list(shape[:-1])) + indices = [tuple(idx/strides % shape) for idx in flat_indices] + return indices + +def _ravel_index_fortran(indices, shape): + """_ravel_index_fortran(flat_indices, shape) + + Given an index tuple into an n-d fortran-strided array, return a flat index. + """ + strides = np.multiply.accumulate([1] + list(shape[:-1])) + flat_indices = [np.sum(strides * idx) for idx in indices] + return flat_indices + +def _normalize_indices(indices, shape): + """_normalize_indices(indices, shape) + + Make all indices positive. If an index is out-of-bounds, return None. + """ + new_indices = [] + for index in indices: + if len(index) != len(shape): + return None + new_index = [] + for i, s in zip(index, shape): + if i < 0: + i = s+i + if not (0 <= i < s): + return None + new_index.append(i) + new_indices.append(new_index) + return new_indices + +cdef class MCP: + """MCP(costs, offsets=None, fully_connected=True) + + A class for finding the minimum cost path through a given n-d costs array. + + Given an n-d costs array, this class can be used to find the minimum-cost + path through that array from any set of points to any other set of points. + Basic usage is to initialize the class and call find_costs() with a one + or more starting indices (and an optional list of end indices). After + that, call traceback() one or more times to find the path from any given + end-position to the closest starting index. New paths through the same + costs array can be found by calling find_costs() repeatedly. + + The cost of a path is calculated simply as the sum of the values of the + `costs` array at each point on the path. The class MCP_Geometric, on the + other hand, accounts for the fact that diagonal vs. axial moves are of + different lengths, and weights the path cost accordingly. + + Parameters + ---------- + costs : ndarray + offsets : iterable, optional + A list of offset tuples: each offset specifies a valid move from a + given n-d position. + If not provided, offsets corresponding to a singly- or fully-connected + n-d neighborhood will be constructed with make_offsets(), using the + `fully_connected` parameter value. + fully_connected : bool, optional + If no `offsets` are provided, this determines the connectivity of the + generated neighborhood. If true, the path may go along diagonals + between elements of the `costs` array; otherwise only axial moves are + permitted. + + Attributes + ---------- + offsets : ndarray + Equivalent to the `offsets` provided to the constructor, or if none + were so provided, the offsets created for the requested n-d + neighborhood. These are useful for interpreting the `traceback` array + returned by the find_costs() method. + """ + def __init__(self, costs, offsets=None, fully_connected=True): + """__init__(costs, offsets=None, fully_connected=True) + + See class documentation. + """ + costs = np.asarray(costs) + if costs.min() < 0: + raise ValueError('minimum cost must be positive') + if not np.can_cast(costs.dtype, FLOAT): + raise TypeError('cannot cast costs array to ' + str(FLOAT)) + # We use flat, fortran-style indexing here (could use C-style, but this is + # my code and I like fortran-style! Also, it's faster when working with + # image arrays, which are often already fortran-strided.) + self.flat_costs = costs.astype(FLOAT).flatten('F') + size = self.flat_costs.shape[0] + self.flat_cumulative_costs = np.empty(size, dtype=FLOAT) + self.flat_cumulative_costs.fill(np.inf) + self.dim = len(costs.shape) + self.costs_shape = costs.shape + self.costs_heap = heap.FastUpdateBinaryHeap(initial_capacity=size, max_reference=size-1) + + # This array stores, for each point, the index into the offset array (see + # below) that leads to that point from the predecessor point. + self.traceback_offsets = np.empty(size, dtype=np.int16) + self.traceback_offsets.fill(-1) + + # The offsets are a list of relative offsets from a central point to each + # point in the relevant neighborhood. (e.g. (-1, 0) might be a 2d offset). + # These offsets are raveled to provide flat, 1d offsets that can be used + # in the same way for flat indices to move to neighboring points. + if offsets is None: + offsets = make_offsets(self.dim, fully_connected) + self.offsets = np.array(offsets, dtype=np.int8) + self.flat_offsets = np.array(_ravel_index_fortran(self.offsets, self.costs_shape), dtype=np.int32) + + # Instead of unraveling each index during the pathfinding algorithm, we + # will use a pre-computed "edge map" that specifies for each dimension + # whether a given index is on a lower or upper boundary (or none at all) + # Flatten this map to get something that can be indexed as by the same + # flat indices as elsewhere. + # The edge map stores more than a boolean "on some edge" flag so as to + # allow us to examine the non-out-of-bounds neighbors for a given edge + # point while excluding the neighbors which are outside the array. + self.flat_edge_map = _offset_edge_map(costs.shape, self.offsets).reshape((size, self.dim), order='F') + + # The offset lengths are the distances traveled along each offset + self.offset_lengths = np.sqrt(np.sum(self.offsets**2, axis=1)).astype(FLOAT) + self.dirty = 0 + self.use_start_cost = 1 + + def _reset(self): + """_reset() + + Clears paths found by find_costs(). + """ + self.costs_heap.reset() + self.traceback_offsets.fill(-1) + self.flat_cumulative_costs.fill(np.inf) + self.dirty = 0 + + cdef FLOAT_C _travel_cost(self, FLOAT_C old_cost, FLOAT_C new_cost, FLOAT_C offset_length): + return new_cost + + @cython.boundscheck(False) + def find_costs(self, starts, ends=None, find_all_ends=True): + """find_costs(starts, ends=None, find_all_ends=True) + + Find the minimum-cost path from the given starting points. + + This method finds the minimum-cost path to the specified ending + indices from any one of the specified starting indices. If no end + positions are given, then the minimum-cost path to every position in + the costs array will be found. + + Parameters + ---------- + starts : iterable + A list of n-d starting indices (where n is the dimension of the + `costs` array). The minimum cost path to the closest/cheapest + starting point will be found. + ends : iterable, optional + A list of n-d ending indices. + find_all_ends : bool, optional + If 'True' (default), the minimum-cost-path to every specified + end-position will be found; otherwise the algorithm will stop when + a a path is found to any end-position. (If no `ends` were + specified, then this parameter has no effect.) + + Returns + ------- + cumulative_costs : ndarray + Same shape as the `costs` array; this array records the minimum + cost path from the nearest/cheapest starting index to each index + considered. (If `ends` were specified, not all elements in the + array will necessarily be considered: positions not evaluated will + have a cumulative cost of inf. If `find_all_ends` is 'False', only + one of the specified end-positions will have a finite cumulative + cost.) + traceback : ndarray + Same shape as the `costs` array; this array contains the offset to + any given index from its predecessor index. The offset indices + index into the `offsets` attribute, which is a array of n-d + offsets. In the 2-d case, if offsets[traceback[x, y]] is (-1, -1), + that means that the predecessor of [x, y] in the minimum cost path + to some start position is [x+1, y+1]. Note that if the + offset_index is -1, then the given index was not considered. + """ + + # basic variables to use for end-finding; also fix up the start and end + # lists + cdef int use_ends = 0 + cdef int num_ends + cdef int all_ends = find_all_ends + cdef np.ndarray[np.uint32_t, ndim=1] flat_ends + starts = _normalize_indices(starts, self.costs_shape) + if starts is None: + raise ValueError('start points must all be within the costs array') + if ends is not None: + ends = _normalize_indices(ends, self.costs_shape) + if ends is None: + raise ValueError('end points must all be within the costs array') + use_ends = 1 + num_ends = len(ends) + flat_ends = np.array(_ravel_index_fortran(ends, self.costs_shape), dtype=np.uint32) + + if self.dirty: + self._reset() + + # lookup and array-ify object attributes for fast use + cdef heap.FastUpdateBinaryHeap costs_heap = self.costs_heap + cdef np.ndarray[FLOAT_T, ndim=1] flat_costs = self.flat_costs + cdef np.ndarray[FLOAT_T, ndim=1] flat_cumulative_costs = self.flat_cumulative_costs + cdef np.ndarray[np.int16_t, ndim=1] traceback_offsets = self.traceback_offsets + cdef np.ndarray[np.int8_t, ndim=2] flat_edge_map = self.flat_edge_map + cdef np.ndarray[np.int8_t, ndim=2] offsets = self.offsets + cdef np.ndarray[np.int32_t, ndim=1] flat_offsets = self.flat_offsets + cdef np.ndarray[FLOAT_T, ndim=1] offset_lengths = self.offset_lengths + + cdef int dim = self.dim + cdef int num_offsets = len(flat_offsets) + + # push each start point into the heap. Note that we use flat indexing! + for start in _ravel_index_fortran(starts, self.costs_shape): + if self.use_start_cost: + costs_heap.push_fast(flat_costs[start], start) + else: + costs_heap.push_fast(0, start) + + cdef double cost, new_cost + cdef unsigned int index, new_index + cdef int is_at_edge, use_offset + cdef unsigned int d, i + cdef int offset, edge_val + cdef int num_ends_found = 0 + cdef double inf = np.inf + cdef double travel_cost + while 1: + # Find the point with the minimum cost in the heap. Once popped, this + # point's minimum cost path has been found. + if costs_heap.count == 0: + # nothing in the heap: we've found paths to every point in the array + break + + cost = costs_heap.pop_fast() + index = costs_heap._popped_ref + + # Record the cost we found to this point + flat_cumulative_costs[index] = cost + + if use_ends: + # If we're only tracing out a path to one or more endpoints, check to + # see if this is an endpoint, and if so, if we're done pathfinding. + for i in range(num_ends): + if index == flat_ends[i]: + num_ends_found += 1 + break + if (num_ends_found and not all_ends) or num_ends_found == num_ends: + # if we've found one or all of the end points (as requested), stop searching + break + + # Look into the edge map to see if this point is at an edge along any axis + is_at_edge = 0 + for d in range(dim): + if flat_edge_map[index, d] != 0: + is_at_edge = 1 + break + + # Now examine the points neighboring the given point + for i in range(num_offsets): + # First, if we're at some edge, scrutinize the offset to ensure that + # it won't put us out-of-bounds. If, for example, the edge_map at + # (x, y) is (-1, 0) -- though of course we use flat indexing below -- + # that means that (x, y) is along the lower edge of the array; thus + # offsets with -1 or more negative in the x-dimension should not be used! + use_offset = 1 + if is_at_edge: + for d in range(dim): + offset = offsets[i, d] + edge_val = flat_edge_map[index, d] + if ( (offset < 0 and edge_val < 0 and offset <= edge_val) + or (offset > 0 and edge_val > 0 and offset >= edge_val)): + use_offset = 0 + break + # If not at an edge, or the specific offset doesn't push over + # the edge, then we go on. + if not use_offset: + continue + + # using the flat offsets, calculate the new flat index + new_index = index + flat_offsets[i] + + # If we have already found the best path here then ignore this point + if flat_cumulative_costs[new_index] != inf: + continue + + # Now we ask the heap to append or update the cost to this new point, + # but only if that point isn't already in the heap, or it is but the + # new cost is lower. + travel_cost = self._travel_cost(flat_costs[index], flat_costs[new_index], offset_lengths[i]) + costs_heap.push_if_lower_fast(cost + travel_cost, new_index) + # If we did perform an append or update, we should record the offset + # from the predecessor to this new point + if costs_heap._pushed: + traceback_offsets[new_index] = i + + # Un-flatten the costs and traceback arrays for human consumption. + cumulative_costs = flat_cumulative_costs.reshape(self.costs_shape, order='F') + traceback = traceback_offsets.reshape(self.costs_shape, order='F') + self.dirty = 1 + return cumulative_costs, traceback + + @cython.boundscheck(False) + def traceback(self, end): + """traceback(end) + + Trace a minimum cost path through the pre-calculated traceback array. + + This convenience function reconstructs the the minimum cost path to a + given end position from one of the starting indices provided to + find_costs(), which must have been called previously. This function + can be called as many times as desired after find_costs() has been + run. + + Parameters + ---------- + end : iterable + An n-d index into the `costs` array. + + Returns + ------- + traceback : list of n-d tuples + A list of indices into the `costs` array, starting with one of + the start positions passed to find_costs(), and ending with the + given `end` index. These indices specify the minimum-cost path + from any given start index to the `end` index. (The total cost + of that path can be read out from the `cumulative_costs` array + returned by find_costs().) + """ + if not self.dirty: + raise Exception('find_costs() must be run before traceback()') + ends = _normalize_indices([end], self.costs_shape) + if ends is None: + raise ValueError('the specified end point must be within the costs array') + traceback = [tuple(ends[0])] + + cdef unsigned int flat_position = _ravel_index_fortran(ends, self.costs_shape)[0] + if self.flat_cumulative_costs[flat_position] == np.inf: + raise ValueError('no minimum-cost path was found to the specified end point') + + cdef np.ndarray[np.int32_t, ndim=1] position = np.array(ends[0], dtype=np.int32) + cdef np.ndarray[np.int16_t, ndim=1] traceback_offsets = self.traceback_offsets + cdef np.ndarray[np.int8_t, ndim=2] offsets = self.offsets + cdef np.ndarray[np.int32_t, ndim=1] flat_offsets = self.flat_offsets + + cdef unsigned int offset, d + cdef int dim = self.dim + while 1: + offset = traceback_offsets[flat_position] + if offset == -1: + # At a point where we can go no further: probably a start point + break + flat_position -= flat_offsets[offset] + for d in range(dim): + position[d] -= offsets[offset, d] + traceback.append(tuple(position)) + return traceback[::-1] + +cdef class MCP_Geometric(MCP): + """MCP_Geometric(costs, offsets=None, fully_connected=True) + + Find distance-weighted minimum cost 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 costs array contains at each position + the "cost" of a unit distance of travel through that position. For + example, a move (in 2-d) from (1, 1) to (1, 2) is assumed to originate in + the center of the pixel (1, 1) and terminate in the center of (1, 2). The + entire move is of distance 1, half through (1, 1) and half through (1, 2); + thus the cost of that move is `(1/2)*costs[1,1] + (1/2)*costs[1,2]`. + + 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. + """ + + def __init__(self, costs, offsets=None, fully_connected=True): + """__init__(costs, offsets=None, fully_connected=True) + + See class documentation. + """ + MCP.__init__(self, costs, offsets, fully_connected) + if np.absolute(self.offsets).max() > 1: + raise ValueError('all offset components must be 0, 1, or -1') + self.use_start_cost = 0 + + cdef FLOAT_C _travel_cost(self, FLOAT_C old_cost, FLOAT_C new_cost, FLOAT_C offset_length): + return offset_length * 0.5 * (old_cost + new_cost) + \ No newline at end of file diff --git a/scikits/image/graph/mcp.py b/scikits/image/graph/mcp.py new file mode 100644 index 00000000..523fb6c4 --- /dev/null +++ b/scikits/image/graph/mcp.py @@ -0,0 +1,41 @@ +from _mcp import MCP, MCP_Geometric, make_offsets + +def route_through_array(array, start, end, fully_connected=True, geometric=True): + """route_through_array(array, start, end, fully_connected=True, geometric=True) + + Simple example of how to use the MCP and MCP_Geometric classes. + + See the MCP and MCP_Geometric class documentation for explanation of the + path-finding algorithm. + + Parameters + ---------- + array : ndarray + Array of costs. + start : iterable + n-d index into `array` defining the starting point + end : iterable + n-d index into `array` defining the end point + fully_connected : bool (optional) + If True, diagonal moves are permitted, if False, only axial moves. + geometric : bool (optional) + If True, the MCP_Geometric class is used to calculate costs, if False, + the MCP base class is used. See the class documentation for + an explanation of the differences between MCP and MCP_Geometric. + + Returns + ------- + path : list + List of n-d index tuples defining the path from `start` to `end`. + cost : float + Cost of the path. + """ + start, end = tuple(start), tuple(end) + if geometric: + mcp_class = MCP_Geometric + else: + mcp_class = MCP + m = mcp_class(array, fully_connected=fully_connected) + costs, traceback_array = m.find_costs([start], [end]) + return m.traceback(end), costs[end] + diff --git a/scikits/image/graph/tests/test_mcp.py b/scikits/image/graph/tests/test_mcp.py new file mode 100644 index 00000000..81018665 --- /dev/null +++ b/scikits/image/graph/tests/test_mcp.py @@ -0,0 +1,113 @@ +import numpy as np +from numpy.testing import * + +import scikits.image.graph.mcp as mcp + +a = np.ones((8,8), dtype=np.float32) +a[1:-1, 1] = 0 +a[1, 1:-1] = 0 + +## 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) + +def test_basic(): + m = mcp.MCP(a, fully_connected=True) + costs, traceback = m.find_costs([(1,6)]) + return_path = m.traceback((7, 2)) + assert_array_equal(costs, + [[ 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.]]) + + assert_array_equal(return_path, + [(1, 6), + (1, 5), + (1, 4), + (1, 3), + (1, 2), + (2, 1), + (3, 1), + (4, 1), + (5, 1), + (6, 1), + (7, 2)]) + +def test_route(): + return_path, cost = mcp.route_through_array(a, (1,6), (7,2), geometric=True) + assert_almost_equal(cost, np.sqrt(2)/2) + assert_array_equal(return_path, + [(1, 6), + (1, 5), + (1, 4), + (1, 3), + (1, 2), + (2, 1), + (3, 1), + (4, 1), + (5, 1), + (6, 1), + (7, 2)]) + +def test_no_diagonal(): + m = mcp.MCP(a, fully_connected=False) + costs, traceback = m.find_costs([(1,6)]) + return_path = m.traceback((7, 2)) + assert_array_equal(costs, + [[ 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.]]) + assert_array_equal(return_path, + [(1, 6), + (1, 5), + (1, 4), + (1, 3), + (1, 2), + (1, 1), + (2, 1), + (3, 1), + (4, 1), + (5, 1), + (6, 1), + (7, 1), + (7, 2)]) + + +def test_crashing(): + _test_random((1000,1000)) + _test_random((10,20,30,40)) + +def _test_random(shape): + # Just tests for crashing -- not for correctness. + np.random.seed(0) + a = np.random.random(shape).astype(np.float32) + starts = [[0]*len(shape), [-1]*len(shape), (np.random.random(len(shape))*shape).astype(int)] + ends = [(np.random.random(len(shape))*shape).astype(int) for i in range(4)] + m = mcp.MCP(a, fully_connected=True) + costs, offsets = m.find_costs(starts) + for point in [(np.random.random(len(shape))*shape).astype(int) for i in range(4)]: + m.traceback(point) + m.reset() + m.find_costs(starts, ends) + for end in ends: + m.traceback(end) + return a, costs, offsets + + +if __name__ == "__main__": + run_module_suite()