mirror of
https://github.com/wassname/scikit-image.git
synced 2026-07-30 12:31:08 +08:00
graph: mcp: Clean up long lines. Fix tests.
This commit is contained in:
+171
-131
@@ -1,3 +1,5 @@
|
||||
# -*- python -*-
|
||||
|
||||
"""Cython implementation of Dijkstra's minimum cost path algorithm,
|
||||
for use with data on a n-dimensional lattice.
|
||||
|
||||
@@ -30,8 +32,7 @@ THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
|
||||
THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
|
||||
"""
|
||||
|
||||
|
||||
import cython
|
||||
import cython
|
||||
cimport numpy as np
|
||||
import numpy as np
|
||||
cimport heap
|
||||
@@ -40,17 +41,16 @@ import heap
|
||||
FLOAT = np.float64
|
||||
|
||||
def _get_edge_map(shape):
|
||||
"""_get_edge_map(shape)
|
||||
|
||||
Return an array with edge points/lines/planes/hyperplanes marked.
|
||||
|
||||
"""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)
|
||||
@@ -64,8 +64,8 @@ def _get_edge_map(shape):
|
||||
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.
|
||||
"""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
|
||||
@@ -105,36 +105,36 @@ def _offset_edge_map(shape, offsets):
|
||||
if offset > 0:
|
||||
slice_stop -= 1
|
||||
slice_step = -np.sign(offset)
|
||||
slices[i] = slice(None, slice_stop, slice_step)
|
||||
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.
|
||||
|
||||
"""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)
|
||||
@@ -154,8 +154,9 @@ def make_offsets(d, fully_connected):
|
||||
|
||||
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]
|
||||
@@ -163,8 +164,9 @@ def _unravel_index_fortran(flat_indices, shape):
|
||||
|
||||
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]
|
||||
@@ -172,8 +174,9 @@ def _ravel_index_fortran(indices, shape):
|
||||
|
||||
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:
|
||||
@@ -191,9 +194,9 @@ def _normalize_indices(indices, shape):
|
||||
|
||||
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
|
||||
@@ -201,12 +204,12 @@ cdef class MCP:
|
||||
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
|
||||
`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.
|
||||
|
||||
different lengths, and weights the path cost accordingly.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
costs : ndarray
|
||||
@@ -221,7 +224,7 @@ cdef class MCP:
|
||||
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
|
||||
@@ -229,10 +232,11 @@ cdef class MCP:
|
||||
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)
|
||||
@@ -240,75 +244,85 @@ cdef class MCP:
|
||||
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.)
|
||||
|
||||
# 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.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).
|
||||
|
||||
# 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)
|
||||
|
||||
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
|
||||
# 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')
|
||||
|
||||
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.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):
|
||||
|
||||
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
|
||||
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.
|
||||
@@ -317,9 +331,9 @@ cdef class MCP:
|
||||
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
|
||||
@@ -328,7 +342,7 @@ cdef class MCP:
|
||||
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
|
||||
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
|
||||
@@ -336,8 +350,8 @@ cdef class MCP:
|
||||
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
|
||||
@@ -350,34 +364,38 @@ cdef class MCP:
|
||||
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')
|
||||
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)
|
||||
|
||||
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[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
|
||||
@@ -387,76 +405,93 @@ cdef class MCP:
|
||||
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.
|
||||
# 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
|
||||
# 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.
|
||||
# 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
|
||||
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
|
||||
|
||||
# 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
|
||||
|
||||
# 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!
|
||||
# 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)):
|
||||
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 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 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])
|
||||
|
||||
# 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 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')
|
||||
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
|
||||
@@ -464,27 +499,27 @@ cdef class MCP:
|
||||
@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
|
||||
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().)
|
||||
"""
|
||||
@@ -492,18 +527,23 @@ cdef class MCP:
|
||||
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')
|
||||
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]
|
||||
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
|
||||
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:
|
||||
@@ -519,39 +559,39 @@ cdef class MCP:
|
||||
|
||||
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
|
||||
|
||||
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
|
||||
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
|
||||
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
|
||||
|
||||
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):
|
||||
|
||||
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)
|
||||
|
||||
@@ -1,13 +1,11 @@
|
||||
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.
|
||||
|
||||
"""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
|
||||
@@ -20,9 +18,9 @@ def route_through_array(array, start, end, fully_connected=True, geometric=True)
|
||||
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
|
||||
the MCP base class is used. See the class documentation for
|
||||
an explanation of the differences between MCP and MCP_Geometric.
|
||||
|
||||
|
||||
Returns
|
||||
-------
|
||||
path : list
|
||||
@@ -38,4 +36,3 @@ def route_through_array(array, start, end, fully_connected=True, geometric=True)
|
||||
m = mcp_class(array, fully_connected=fully_connected)
|
||||
costs, traceback_array = m.find_costs([start], [end])
|
||||
return m.traceback(end), costs[end]
|
||||
|
||||
|
||||
@@ -12,7 +12,7 @@ def configuration(parent_package='', top_path=None):
|
||||
config.add_data_dir('tests')
|
||||
|
||||
# This function tries to create C files from the given .pyx files. If
|
||||
# it fails, we build the checked-in .c files.
|
||||
# it fails, try to build with pre-generated .c files.
|
||||
cython(['_spath.pyx'], working_path=base_path)
|
||||
cython(['_mcp.pyx'], working_path=base_path)
|
||||
cython(['heap.pyx'], working_path=base_path)
|
||||
@@ -34,4 +34,4 @@ if __name__ == '__main__':
|
||||
url = 'http://stefanv.github.com/scikits.image/',
|
||||
license = 'Modified BSD',
|
||||
**(configuration(top_path='').todict())
|
||||
)
|
||||
)
|
||||
|
||||
@@ -1,12 +1,12 @@
|
||||
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.],
|
||||
@@ -15,7 +15,7 @@ a[1, 1:-1] = 0
|
||||
## [ 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)])
|
||||
@@ -29,7 +29,7 @@ def test_basic():
|
||||
[ 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),
|
||||
@@ -58,7 +58,7 @@ def test_route():
|
||||
(5, 1),
|
||||
(6, 1),
|
||||
(7, 2)])
|
||||
|
||||
|
||||
def test_no_diagonal():
|
||||
m = mcp.MCP(a, fully_connected=False)
|
||||
costs, traceback = m.find_costs([(1,6)])
|
||||
@@ -86,23 +86,25 @@ def test_no_diagonal():
|
||||
(6, 1),
|
||||
(7, 1),
|
||||
(7, 2)])
|
||||
|
||||
|
||||
|
||||
def test_crashing():
|
||||
_test_random((1000,1000))
|
||||
_test_random((10,20,30,40))
|
||||
for shape in [(100, 100), (5, 8, 13, 17)]:
|
||||
yield _test_random, shape
|
||||
|
||||
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)]
|
||||
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)]:
|
||||
for point in [(np.random.random(len(shape))*shape).astype(int)
|
||||
for i in range(4)]:
|
||||
m.traceback(point)
|
||||
m.reset()
|
||||
m._reset()
|
||||
m.find_costs(starts, ends)
|
||||
for end in ends:
|
||||
m.traceback(end)
|
||||
|
||||
@@ -1,8 +1,8 @@
|
||||
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],
|
||||
@@ -27,7 +27,7 @@ def test_non_square():
|
||||
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()
|
||||
|
||||
Reference in New Issue
Block a user