mirror of
https://github.com/wassname/scikit-image.git
synced 2026-07-20 12:40:31 +08:00
Merge pull request #854 from almarklein/update-mcp
Updates to MCP algorithm. Improves speed, adds anisotropy, more flexibility in general, and a new MCP class to find connections between seed points.
This commit is contained in:
Stefan van der Walt
@@ -1,7 +1,9 @@
|
||||
from .spath import shortest_path
|
||||
from .mcp import MCP, MCP_Geometric, route_through_array
|
||||
from .mcp import MCP, MCP_Geometric, MCP_Connect, MCP_Flexible, route_through_array
|
||||
|
||||
__all__ = ['shortest_path',
|
||||
'MCP',
|
||||
'MCP_Geometric',
|
||||
'MCP_Connect',
|
||||
'MCP_Flexible',
|
||||
'route_through_array']
|
||||
+23
-9
@@ -9,22 +9,36 @@ cimport numpy as cnp
|
||||
ctypedef heap.BOOL_T BOOL_T
|
||||
ctypedef unsigned char DIM_T
|
||||
ctypedef cnp.float64_t FLOAT_T
|
||||
ctypedef cnp.intp_t INDEX_T
|
||||
ctypedef cnp.int8_t EDGE_T
|
||||
ctypedef cnp.int8_t OFFSET_T
|
||||
ctypedef cnp.int16_t OFFSETS_INDEX_T
|
||||
|
||||
|
||||
cdef class MCP:
|
||||
cdef heap.FastUpdateBinaryHeap costs_heap
|
||||
cdef object costs_shape
|
||||
cdef object _starts
|
||||
cdef object _ends
|
||||
cdef DIM_T dim
|
||||
cdef object flat_costs
|
||||
cdef object flat_cumulative_costs
|
||||
cdef object traceback_offsets
|
||||
cdef object flat_pos_edge_map
|
||||
cdef object flat_neg_edge_map
|
||||
cdef readonly object offsets
|
||||
cdef object flat_offsets
|
||||
cdef object offset_lengths
|
||||
cdef BOOL_T dirty
|
||||
cdef BOOL_T 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...
|
||||
|
||||
|
||||
# Arrays used during front propagation
|
||||
cdef FLOAT_T [:] flat_costs
|
||||
cdef FLOAT_T [:] flat_cumulative_costs
|
||||
cdef OFFSETS_INDEX_T [:] traceback_offsets
|
||||
cdef EDGE_T [:,:] flat_pos_edge_map
|
||||
cdef EDGE_T [:,:] flat_neg_edge_map
|
||||
cdef OFFSET_T [:,:] offsets
|
||||
cdef INDEX_T [:] flat_offsets
|
||||
cdef FLOAT_T [:] offset_lengths
|
||||
|
||||
# Methods
|
||||
cpdef int goal_reached(self, INDEX_T index, FLOAT_T cumcost)
|
||||
cdef FLOAT_T _travel_cost(self, FLOAT_T old_cost, FLOAT_T new_cost, FLOAT_T offset_length)
|
||||
cdef void _examine_neighbor(self, INDEX_T index, INDEX_T new_index, FLOAT_T offset_length)
|
||||
cdef void _update_node(self, INDEX_T index, INDEX_T new_index, FLOAT_T offset_length)
|
||||
|
||||
+380
-99
@@ -5,6 +5,7 @@ for use with data on a n-dimensional lattice.
|
||||
|
||||
Original author: Zachary Pincus
|
||||
Inspired by code from Almar Klein
|
||||
Later modifications by Almar Klein (Dec 2013)
|
||||
|
||||
License: BSD
|
||||
|
||||
@@ -39,13 +40,9 @@ import heap
|
||||
cimport numpy as cnp
|
||||
cimport heap
|
||||
|
||||
ctypedef cnp.int8_t OFFSET_T
|
||||
OFFSET_D = np.int8
|
||||
ctypedef cnp.int16_t OFFSETS_INDEX_T
|
||||
OFFSETS_INDEX_D = np.int16
|
||||
ctypedef cnp.int8_t EDGE_T
|
||||
EDGE_D = np.int8
|
||||
ctypedef cnp.intp_t INDEX_T
|
||||
INDEX_D = np.intp
|
||||
FLOAT_D = np.float64
|
||||
|
||||
@@ -106,7 +103,7 @@ def _offset_edge_map(shape, offsets):
|
||||
[0, 0, 2, 1]], dtype=int8)
|
||||
|
||||
"""
|
||||
indices = np.indices(shape) # indices.shape = (n,)+shape
|
||||
indices = np.indices(shape) # indices.shape = (n,)+shape
|
||||
|
||||
#get the distance from each index to the upper or lower edge in each dim
|
||||
pos_edges = (shape - indices.T).T
|
||||
@@ -172,11 +169,12 @@ 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.
|
||||
|
||||
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]
|
||||
indices = [tuple((idx // strides) % shape) for idx in flat_indices]
|
||||
return indices
|
||||
|
||||
|
||||
@@ -185,7 +183,8 @@ 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.
|
||||
Given an index tuple into an n-d fortran-strided array, return a
|
||||
flat index.
|
||||
|
||||
"""
|
||||
strides = np.multiply.accumulate([1] + list(shape[:-1]))
|
||||
@@ -209,7 +208,7 @@ def _normalize_indices(indices, shape):
|
||||
for i, s in zip(index, shape):
|
||||
i = int(i)
|
||||
if i < 0:
|
||||
i = s+i
|
||||
i = s + i
|
||||
if not (0 <= i < s):
|
||||
return None
|
||||
new_index.append(i)
|
||||
@@ -220,14 +219,15 @@ def _normalize_indices(indices, shape):
|
||||
@cython.boundscheck(True)
|
||||
@cython.wraparound(True)
|
||||
def _reverse(arr):
|
||||
"""Reverse index an array safely, with bounds/wraparound checks on."""
|
||||
"""Reverse index an array safely, with bounds/wraparound checks on.
|
||||
"""
|
||||
return arr[::-1]
|
||||
|
||||
|
||||
@cython.boundscheck(False)
|
||||
@cython.wraparound(False)
|
||||
cdef class MCP:
|
||||
"""MCP(costs, offsets=None, fully_connected=True)
|
||||
"""MCP(costs, offsets=None, fully_connected=True, sampling=None)
|
||||
|
||||
A class for finding the minimum cost path through a given n-d costs array.
|
||||
|
||||
@@ -261,7 +261,10 @@ 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.
|
||||
|
||||
sampling : tuple, optional
|
||||
For each dimension, specifies the distance between two cells/voxels.
|
||||
If not given or None, the distance is assumed unit.
|
||||
|
||||
Attributes
|
||||
----------
|
||||
offsets : ndarray
|
||||
@@ -271,15 +274,27 @@ cdef class MCP:
|
||||
returned by the find_costs() method.
|
||||
|
||||
"""
|
||||
def __init__(self, costs, offsets=None, fully_connected=True):
|
||||
"""__init__(costs, offsets=None, fully_connected=True)
|
||||
|
||||
def __init__(self, costs, offsets=None, fully_connected=True,
|
||||
sampling=None):
|
||||
"""__init__(costs, offsets=None, fully_connected=True, sampling=None)
|
||||
|
||||
See class documentation.
|
||||
"""
|
||||
costs = np.asarray(costs)
|
||||
if not np.can_cast(costs.dtype, FLOAT_D):
|
||||
raise TypeError('cannot cast costs array to ' + str(FLOAT_D))
|
||||
|
||||
|
||||
# Check sampling
|
||||
if sampling is None:
|
||||
sampling = np.array([1.0 for s in costs.shape], FLOAT_D)
|
||||
elif isinstance(sampling, (list, tuple)):
|
||||
sampling = np.array(sampling, FLOAT_D)
|
||||
if sampling.ndim != 1 or len(sampling) != costs.ndim:
|
||||
raise ValueError('Need one sampling element per dimension.')
|
||||
else:
|
||||
raise ValueError('Invalid type for sampling: %r.' % type(sampling))
|
||||
|
||||
# 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
|
||||
@@ -287,23 +302,21 @@ cdef class MCP:
|
||||
self.flat_costs = costs.astype(FLOAT_D).flatten('F')
|
||||
size = self.flat_costs.shape[0]
|
||||
self.flat_cumulative_costs = np.empty(size, dtype=FLOAT_D)
|
||||
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,
|
||||
self.costs_heap = heap.FastUpdateBinaryHeap(initial_capacity=128,
|
||||
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=OFFSETS_INDEX_D)
|
||||
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.
|
||||
# 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=OFFSET_D)
|
||||
@@ -313,9 +326,9 @@ cdef class MCP:
|
||||
|
||||
# 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.
|
||||
# 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.
|
||||
@@ -325,26 +338,80 @@ cdef class MCP:
|
||||
|
||||
|
||||
# The offset lengths are the distances traveled along each offset
|
||||
self.offset_lengths = np.sqrt(
|
||||
np.sum(self.offsets**2, axis=1)).astype(FLOAT_D)
|
||||
self.offset_lengths = np.sqrt(np.sum((sampling * self.offsets)**2,
|
||||
axis=1)).astype(FLOAT_D)
|
||||
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.traceback_offsets[...] = -2 # -2 is not reached, -1 is start
|
||||
self.flat_cumulative_costs[...] = np.inf
|
||||
self.dirty = 0
|
||||
|
||||
|
||||
# Get starts and ends
|
||||
# We do not pass them in as arguments for backwards compat
|
||||
starts, ends = self._starts, self._ends
|
||||
|
||||
# push each start point into the heap. Note that we use flat indexing!
|
||||
for start in _ravel_index_fortran(starts, self.costs_shape):
|
||||
self.traceback_offsets[start] = -1
|
||||
if self.use_start_cost:
|
||||
self.costs_heap.push_fast(self.flat_costs[start], start)
|
||||
else:
|
||||
self.costs_heap.push_fast(0, start)
|
||||
|
||||
|
||||
cdef FLOAT_T _travel_cost(self, FLOAT_T old_cost,
|
||||
FLOAT_T new_cost, FLOAT_T offset_length):
|
||||
""" float _travel_cost(float old_cost, float new_cost,
|
||||
float offset_length)
|
||||
The travel cost for going from the current node to the next.
|
||||
Default is simply the cost of the next node.
|
||||
"""
|
||||
return new_cost
|
||||
|
||||
def find_costs(self, starts, ends=None, find_all_ends=True):
|
||||
|
||||
|
||||
cpdef int goal_reached(self, INDEX_T index, FLOAT_T cumcost):
|
||||
""" int goal_reached(int index, float cumcost)
|
||||
This method is called each iteration after popping an index
|
||||
from the heap, before examining the neighbours.
|
||||
|
||||
This method can be overloaded to modify the behavior of the MCP
|
||||
algorithm. An example might be to stop the algorithm when a
|
||||
certain cumulative cost is reached, or when the front is a
|
||||
certain distance away from the seed point.
|
||||
|
||||
This method should return 1 if the algorithm should not check
|
||||
the current point's neighbours and 2 if the algorithm is now
|
||||
done.
|
||||
"""
|
||||
return 0
|
||||
|
||||
|
||||
cdef void _examine_neighbor(self, INDEX_T index, INDEX_T new_index,
|
||||
FLOAT_T offset_length):
|
||||
""" _examine_neighbor(int index, int new_index, float offset_length)
|
||||
This method is called once for every pair of neighboring nodes,
|
||||
as soon as both nodes become frozen.
|
||||
"""
|
||||
pass
|
||||
|
||||
|
||||
cdef void _update_node(self, INDEX_T index, INDEX_T new_index,
|
||||
FLOAT_T offset_length):
|
||||
""" _update_node(int index, int new_index, float offset_length)
|
||||
This method is called when a node is updated.
|
||||
"""
|
||||
pass
|
||||
|
||||
|
||||
def find_costs(self, starts, ends=None, find_all_ends=True,
|
||||
max_coverage=1.0, max_cumulative_cost=None, max_cost=None):
|
||||
"""
|
||||
Find the minimum-cost path from the given starting points.
|
||||
|
||||
@@ -392,10 +459,13 @@ cdef class MCP:
|
||||
cdef BOOL_T use_ends = 0
|
||||
cdef INDEX_T num_ends
|
||||
cdef BOOL_T all_ends = find_all_ends
|
||||
cdef cnp.ndarray[INDEX_T, ndim=1] flat_ends
|
||||
cdef INDEX_T [:] flat_ends
|
||||
starts = _normalize_indices(starts, self.costs_shape)
|
||||
if starts is None:
|
||||
raise ValueError('start points must all be within the costs array')
|
||||
elif not starts:
|
||||
raise ValueError('no valid start points to start front' +
|
||||
'propagation')
|
||||
if ends is not None:
|
||||
ends = _normalize_indices(ends, self.costs_shape)
|
||||
if ends is None:
|
||||
@@ -406,57 +476,66 @@ cdef class MCP:
|
||||
flat_ends = np.array(_ravel_index_fortran(
|
||||
ends, self.costs_shape), dtype=INDEX_D)
|
||||
|
||||
if self.dirty:
|
||||
self._reset()
|
||||
|
||||
# lookup and array-ify object attributes for fast use
|
||||
# Always perform a reset to (re)initialize our arrays and start
|
||||
# positions
|
||||
self._starts, self._ends = starts, ends
|
||||
self._reset()
|
||||
|
||||
# Get shorter names for arrays
|
||||
cdef FLOAT_T [:] flat_costs = self.flat_costs
|
||||
cdef FLOAT_T [:] flat_cumulative_costs = self.flat_cumulative_costs
|
||||
cdef OFFSETS_INDEX_T [:] traceback_offsets = self.traceback_offsets
|
||||
cdef EDGE_T [:,:] flat_pos_edge_map = self.flat_pos_edge_map
|
||||
cdef EDGE_T [:,:] flat_neg_edge_map = self.flat_neg_edge_map
|
||||
cdef OFFSET_T [:,:] offsets = self.offsets
|
||||
cdef INDEX_T [:] flat_offsets = self.flat_offsets
|
||||
cdef FLOAT_T [:] offset_lengths = self.offset_lengths
|
||||
|
||||
# Short names for other attributes
|
||||
cdef heap.FastUpdateBinaryHeap costs_heap = self.costs_heap
|
||||
cdef cnp.ndarray[FLOAT_T, ndim=1] flat_costs = self.flat_costs
|
||||
cdef cnp.ndarray[FLOAT_T, ndim=1] flat_cumulative_costs = \
|
||||
self.flat_cumulative_costs
|
||||
cdef cnp.ndarray[OFFSETS_INDEX_T, ndim=1] traceback_offsets = \
|
||||
self.traceback_offsets
|
||||
cdef cnp.ndarray[EDGE_T, ndim=2] flat_pos_edge_map = \
|
||||
self.flat_pos_edge_map
|
||||
cdef cnp.ndarray[EDGE_T, ndim=2] flat_neg_edge_map = \
|
||||
self.flat_neg_edge_map
|
||||
cdef cnp.ndarray[OFFSET_T, ndim=2] offsets = self.offsets
|
||||
cdef cnp.ndarray[INDEX_T, ndim=1] flat_offsets = self.flat_offsets
|
||||
cdef cnp.ndarray[FLOAT_T, ndim=1] offset_lengths = self.offset_lengths
|
||||
|
||||
cdef DIM_T 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 FLOAT_T cost, new_cost
|
||||
|
||||
# Variables used during front propagation
|
||||
cdef FLOAT_T cost, new_cost, cumcost, new_cumcost, offset_length
|
||||
cdef INDEX_T index, new_index
|
||||
cdef BOOL_T is_at_edge, use_offset
|
||||
cdef INDEX_T d, i
|
||||
cdef INDEX_T d, i, iter
|
||||
cdef OFFSET_T offset
|
||||
cdef EDGE_T pos_edge_val, neg_edge_val
|
||||
cdef int num_ends_found = 0
|
||||
cdef FLOAT_T inf = np.inf
|
||||
cdef FLOAT_T travel_cost
|
||||
while 1:
|
||||
cdef int goal_reached
|
||||
|
||||
cdef INDEX_T maxiter = int(max_coverage * flat_costs.size)
|
||||
|
||||
for iter in range(maxiter):
|
||||
|
||||
# This is rather like a while loop, except we are guaranteed to
|
||||
# exit, which is nice during developing to prevent eternal loops.
|
||||
|
||||
# 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()
|
||||
|
||||
# Get current cumulative cost and index from the heap
|
||||
cumcost = costs_heap.pop_fast()
|
||||
index = costs_heap._popped_ref
|
||||
|
||||
# Record the cost we found to this point
|
||||
flat_cumulative_costs[index] = cost
|
||||
|
||||
flat_cumulative_costs[index] = cumcost
|
||||
|
||||
# Check if goal is reached
|
||||
goal_reached = self.goal_reached(index, cumcost)
|
||||
if goal_reached > 0:
|
||||
if goal_reached == 1:
|
||||
continue # Skip neighbours
|
||||
else:
|
||||
break # Done completely
|
||||
|
||||
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
|
||||
@@ -470,7 +549,7 @@ cdef class MCP:
|
||||
# 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
|
||||
@@ -504,58 +583,75 @@ cdef class MCP:
|
||||
# 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]
|
||||
|
||||
|
||||
# Get offset length
|
||||
offset_length = offset_lengths[i]
|
||||
|
||||
# If we have already found the best path here then
|
||||
# ignore this point
|
||||
if flat_cumulative_costs[new_index] != inf:
|
||||
# Give subclass the oportunity to examine these two nodes
|
||||
# Note that only when both nodes are "frozen" their
|
||||
# cumulative cost is set. By doing the check here, each
|
||||
# pair of nodes is checked exactly once.
|
||||
self._examine_neighbor(index, new_index, offset_length)
|
||||
continue
|
||||
# If the cost at this point is negative or infinite, ignore it
|
||||
|
||||
# Get cost and new cost
|
||||
cost = flat_costs[index]
|
||||
new_cost = flat_costs[new_index]
|
||||
|
||||
# If the cost at this point is negative or infinite, ignore it
|
||||
if new_cost < 0 or new_cost == inf:
|
||||
continue
|
||||
|
||||
|
||||
# Calculate new cumulative cost
|
||||
new_cumcost = cumcost + self._travel_cost(cost, new_cost,
|
||||
offset_length)
|
||||
|
||||
# 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],
|
||||
new_cost,
|
||||
offset_lengths[i])
|
||||
# don't push infs into the heap though!
|
||||
new_cost = cost + travel_cost
|
||||
if new_cost != inf:
|
||||
costs_heap.push_if_lower_fast(new_cost, new_index)
|
||||
if new_cumcost != inf:
|
||||
costs_heap.push_if_lower_fast(new_cumcost, 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
|
||||
|
||||
self._update_node(index, new_index, offset_length)
|
||||
|
||||
|
||||
# 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')
|
||||
cumulative_costs = np.asarray(flat_cumulative_costs)
|
||||
cumulative_costs = cumulative_costs.reshape(self.costs_shape,
|
||||
order='F')
|
||||
traceback = np.asarray(traceback_offsets)
|
||||
traceback = traceback.reshape(self.costs_shape, order='F')
|
||||
self.dirty = 1
|
||||
return cumulative_costs, traceback
|
||||
|
||||
|
||||
|
||||
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
|
||||
@@ -579,14 +675,14 @@ cdef class MCP:
|
||||
if self.flat_cumulative_costs[flat_position] == np.inf:
|
||||
raise ValueError('no minimum-cost path was found '
|
||||
'to the specified end point')
|
||||
|
||||
cdef cnp.ndarray[INDEX_T, ndim=1] position = \
|
||||
np.array(ends[0], dtype=INDEX_D)
|
||||
cdef cnp.ndarray[OFFSETS_INDEX_T, ndim=1] traceback_offsets = \
|
||||
self.traceback_offsets
|
||||
cdef cnp.ndarray[OFFSET_T, ndim=2] offsets = self.offsets
|
||||
cdef cnp.ndarray[INDEX_T, ndim=1] flat_offsets = self.flat_offsets
|
||||
|
||||
|
||||
# Short names for arrays
|
||||
cdef OFFSETS_INDEX_T [:] traceback_offsets = self.traceback_offsets
|
||||
cdef OFFSET_T [:,:] offsets = self.offsets
|
||||
cdef INDEX_T [:] flat_offsets = self.flat_offsets
|
||||
# New array
|
||||
cdef INDEX_T [:] position = np.array(ends[0], dtype=INDEX_D)
|
||||
|
||||
cdef OFFSETS_INDEX_T offset
|
||||
cdef DIM_T d
|
||||
cdef DIM_T dim = self.dim
|
||||
@@ -602,6 +698,7 @@ cdef class MCP:
|
||||
return _reverse(traceback)
|
||||
|
||||
|
||||
|
||||
@cython.boundscheck(False)
|
||||
@cython.wraparound(False)
|
||||
cdef class MCP_Geometric(MCP):
|
||||
@@ -626,15 +723,17 @@ cdef class MCP_Geometric(MCP):
|
||||
`(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.
|
||||
greater than 1. Use the `sampling` argument in order to deal with
|
||||
anisotropic data.
|
||||
"""
|
||||
|
||||
def __init__(self, costs, offsets=None, fully_connected=True):
|
||||
"""__init__(costs, offsets=None, fully_connected=True)
|
||||
def __init__(self, costs, offsets=None, fully_connected=True,
|
||||
sampling=None):
|
||||
"""__init__(costs, offsets=None, fully_connected=True, sampling=None)
|
||||
|
||||
See class documentation.
|
||||
"""
|
||||
MCP.__init__(self, costs, offsets, fully_connected)
|
||||
MCP.__init__(self, costs, offsets, fully_connected, sampling)
|
||||
if np.absolute(self.offsets).max() > 1:
|
||||
raise ValueError('all offset components must be 0, 1, or -1')
|
||||
self.use_start_cost = 0
|
||||
@@ -642,3 +741,185 @@ cdef class MCP_Geometric(MCP):
|
||||
cdef FLOAT_T _travel_cost(self, FLOAT_T old_cost, FLOAT_T new_cost,
|
||||
FLOAT_T offset_length):
|
||||
return offset_length * 0.5 * (old_cost + new_cost)
|
||||
|
||||
|
||||
|
||||
@cython.boundscheck(True)
|
||||
@cython.wraparound(True)
|
||||
cdef class MCP_Connect(MCP):
|
||||
"""MCP_Connect(costs, offsets=None, fully_connected=True)
|
||||
|
||||
Connect source points using the distance-weighted minimum cost function.
|
||||
|
||||
A front is grown from each seed point simultaneously, while the
|
||||
origin of the front is tracked as well. When two fronts meet,
|
||||
create_connection() is called. This method must be overloaded to
|
||||
deal with the found edges in a way that is appropriate for the
|
||||
application.
|
||||
"""
|
||||
|
||||
cdef INDEX_T [:] flat_idmap
|
||||
|
||||
|
||||
def __init__(self, costs, offsets=None, fully_connected=True,
|
||||
sampling=None):
|
||||
MCP.__init__(self, costs, offsets, fully_connected, sampling)
|
||||
|
||||
# Create id map to keep track of origin of nodes
|
||||
self.flat_idmap = np.zeros(self.costs_shape, INDEX_D).ravel('F')
|
||||
|
||||
|
||||
def _reset(self):
|
||||
""" Reset the id map.
|
||||
"""
|
||||
MCP._reset(self)
|
||||
starts, ends = self._starts, self._ends
|
||||
|
||||
# Reset idmap
|
||||
self.flat_idmap[...] = -1
|
||||
id = 0
|
||||
for start in _ravel_index_fortran(starts, self.costs_shape):
|
||||
self.flat_idmap[start] = id
|
||||
id += 1
|
||||
|
||||
|
||||
cdef FLOAT_T _travel_cost(self, FLOAT_T old_cost, FLOAT_T new_cost,
|
||||
FLOAT_T offset_length):
|
||||
""" Equivalent to MCP_Geometric.
|
||||
"""
|
||||
return offset_length * 0.5 * (old_cost + new_cost)
|
||||
|
||||
|
||||
cdef void _examine_neighbor(self, INDEX_T index, INDEX_T new_index,
|
||||
FLOAT_T offset_length):
|
||||
""" Check whether two fronts are meeting. If so, the flat_traceback
|
||||
is obtained and a connection is created.
|
||||
"""
|
||||
|
||||
# Short names
|
||||
cdef INDEX_T [:] flat_idmap = self.flat_idmap
|
||||
cdef FLOAT_T [:] flat_cumulative_costs = self.flat_cumulative_costs
|
||||
|
||||
# Get ids
|
||||
cdef INDEX_T id1 = flat_idmap[index]
|
||||
cdef INDEX_T id2 = flat_idmap[new_index]
|
||||
|
||||
if id2 < 0 or id1 < 0:
|
||||
pass
|
||||
elif id2 != id1:
|
||||
# We reached the 'front' of another seed point!
|
||||
# Get position/coordinates
|
||||
pos1, pos2 = _unravel_index_fortran([index, new_index],
|
||||
self.costs_shape)
|
||||
# Also get the costs, so we can keep the path with the least cost
|
||||
cost1 = flat_cumulative_costs[index]
|
||||
cost2 = flat_cumulative_costs[new_index]
|
||||
# Create connection
|
||||
self.create_connection(id1, id2, pos1, pos2, cost1, cost2)
|
||||
|
||||
|
||||
def create_connection(self, id1, id2, tb1, tb2, cost1, cost2):
|
||||
""" create_connection id1, id2, pos1, pos2, cost1, cost2)
|
||||
|
||||
Overload this method to keep track of the connections that are
|
||||
found during MCP processing. Note that a connection with the
|
||||
same ids can be found multiple times (but with different
|
||||
positions and costs).
|
||||
|
||||
At the time that this method is called, both points are "frozen"
|
||||
and will not be visited again by the MCP algorithm.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
id1 : int
|
||||
The seed point id where the first neighbor originated from.
|
||||
id2 : int
|
||||
The seed point id where the second neighbor originated from.
|
||||
pos1 : tuple
|
||||
The index of of the first neighbour in the connection.
|
||||
pos2 : tuple
|
||||
The index of of the second neighbour in the connection.
|
||||
cost1 : float
|
||||
The cumulative cost at `pos1`.
|
||||
cost2 : float
|
||||
The cumulative costs at `pos2`.
|
||||
"""
|
||||
pass
|
||||
|
||||
|
||||
cdef void _update_node(self, INDEX_T index, INDEX_T new_index,
|
||||
FLOAT_T offset_length):
|
||||
""" Keep track of the id map so that we know which seed point
|
||||
a certain front originates from.
|
||||
"""
|
||||
self.flat_idmap[new_index] = self.flat_idmap[index]
|
||||
|
||||
|
||||
|
||||
@cython.boundscheck(False)
|
||||
@cython.wraparound(False)
|
||||
cdef class MCP_Flexible(MCP):
|
||||
"""MCP_Flexible(costs, offsets=None, fully_connected=True)
|
||||
|
||||
Find minimum cost paths through an N-d costs array.
|
||||
|
||||
See the documentation for MCP for full details. This class differs from
|
||||
MCP in that several methods can be overloaded (from pure Python) to
|
||||
modify the behavior of the algorithm and/or create custom algorithms
|
||||
based on MCP. Note that goal_reached can also be overloaded in the
|
||||
MCP class.
|
||||
|
||||
"""
|
||||
|
||||
def travel_cost(self, FLOAT_T old_cost, FLOAT_T new_cost,
|
||||
FLOAT_T offset_length):
|
||||
""" travel_cost(old_cost, new_cost, offset_length)
|
||||
This method calculates the travel cost for going from the
|
||||
current node to the next. The default implementation returns
|
||||
new_cost. Overload this method to adapt the behaviour of the
|
||||
algorithm.
|
||||
"""
|
||||
return new_cost
|
||||
|
||||
|
||||
def examine_neighbor(self, INDEX_T index, INDEX_T new_index,
|
||||
FLOAT_T offset_length):
|
||||
""" examine_neighbor(index, new_index, offset_length)
|
||||
This method is called once for every pair of neighboring nodes,
|
||||
as soon as both nodes are frozen.
|
||||
|
||||
This method can be overloaded to obtain information about
|
||||
neightboring nodes, and/or to modify the behavior of the MCP
|
||||
algorithm. One example is the MCP_Connect class, which checks
|
||||
for meeting fronts using this hook.
|
||||
"""
|
||||
pass
|
||||
|
||||
|
||||
def update_node(self, INDEX_T index, INDEX_T new_index,
|
||||
FLOAT_T offset_length):
|
||||
""" update_node(index, new_index, offset_length)
|
||||
This method is called when a node is updated, right after
|
||||
new_index is pushed onto the heap and the traceback map is
|
||||
updated.
|
||||
|
||||
This method can be overloaded to keep track of other arrays
|
||||
that are used by a specific implementation of the algorithm.
|
||||
For instance the MCP_Connect class uses it to update an id map.
|
||||
"""
|
||||
pass
|
||||
|
||||
|
||||
cdef FLOAT_T _travel_cost(self, FLOAT_T old_cost, FLOAT_T new_cost,
|
||||
FLOAT_T offset_length):
|
||||
return self.travel_cost(old_cost, new_cost, offset_length)
|
||||
|
||||
|
||||
cdef void _examine_neighbor(self, INDEX_T index, INDEX_T new_index,
|
||||
FLOAT_T offset_length):
|
||||
self.examine_neighbor(index, new_index, offset_length)
|
||||
|
||||
|
||||
cdef void _update_node(self, INDEX_T index, INDEX_T new_index,
|
||||
FLOAT_T offset_length):
|
||||
self.update_node(index, new_index, offset_length)
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
from ._mcp import MCP, MCP_Geometric
|
||||
from ._mcp import MCP, MCP_Geometric, MCP_Connect, MCP_Flexible
|
||||
|
||||
|
||||
def route_through_array(array, start, end, fully_connected=True,
|
||||
|
||||
@@ -0,0 +1,53 @@
|
||||
import skimage.graph.mcp as mcp
|
||||
from numpy.testing import (assert_array_equal,
|
||||
assert_almost_equal,
|
||||
)
|
||||
|
||||
import numpy as np
|
||||
|
||||
a = np.ones((8, 8), dtype=np.float32)
|
||||
|
||||
|
||||
horizontal_ramp = np.array([[ 0., 1., 2., 3., 4., 5., 6., 7.,],
|
||||
[ 0., 1., 2., 3., 4., 5., 6., 7.,],
|
||||
[ 0., 1., 2., 3., 4., 5., 6., 7.,],
|
||||
[ 0., 1., 2., 3., 4., 5., 6., 7.,],
|
||||
[ 0., 1., 2., 3., 4., 5., 6., 7.,],
|
||||
[ 0., 1., 2., 3., 4., 5., 6., 7.,],
|
||||
[ 0., 1., 2., 3., 4., 5., 6., 7.,],
|
||||
[ 0., 1., 2., 3., 4., 5., 6., 7.,]])
|
||||
|
||||
vertical_ramp = np.array( [[ 0., 0., 0., 0., 0., 0., 0., 0.,],
|
||||
[ 1., 1., 1., 1., 1., 1., 1., 1.,],
|
||||
[ 2., 2., 2., 2., 2., 2., 2., 2.,],
|
||||
[ 3., 3., 3., 3., 3., 3., 3., 3.,],
|
||||
[ 4., 4., 4., 4., 4., 4., 4., 4.,],
|
||||
[ 5., 5., 5., 5., 5., 5., 5., 5.,],
|
||||
[ 6., 6., 6., 6., 6., 6., 6., 6.,],
|
||||
[ 7., 7., 7., 7., 7., 7., 7., 7.,]])
|
||||
|
||||
|
||||
def test_anisotropy():
|
||||
|
||||
# Create seeds; vertical seeds create a horizonral ramp
|
||||
seeds_for_horizontal = [(i, 0) for i in range(8) ]
|
||||
seeds_for_vertcal = [(0, i) for i in range(8) ]
|
||||
|
||||
|
||||
for sy in range(1, 5):
|
||||
for sx in range(1,5):
|
||||
sampling = sy, sx
|
||||
# Trace horizontally
|
||||
m1 = mcp.MCP_Geometric(a, sampling=sampling, fully_connected=True)
|
||||
costs1, traceback = m1.find_costs(seeds_for_horizontal)
|
||||
# Trace vertically
|
||||
m2 = mcp.MCP_Geometric(a, sampling=sampling, fully_connected=True)
|
||||
costs2, traceback = m2.find_costs(seeds_for_vertcal)
|
||||
|
||||
# Check
|
||||
assert_array_equal(costs1, horizontal_ramp * sx)
|
||||
assert_array_equal(costs2, vertical_ramp * sy)
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
np.testing.run_module_suite()
|
||||
@@ -0,0 +1,83 @@
|
||||
import skimage.graph.mcp as mcp
|
||||
# import stentseg.graph._mcp as mcp
|
||||
from numpy.testing import (assert_array_equal,
|
||||
assert_almost_equal,
|
||||
)
|
||||
|
||||
import numpy as np
|
||||
|
||||
|
||||
a = np.ones((8, 8), dtype=np.float32)
|
||||
|
||||
count = 0
|
||||
class MCP(mcp.MCP_Connect):
|
||||
|
||||
def _reset(self):
|
||||
""" Reset the id map.
|
||||
"""
|
||||
mcp.MCP_Connect._reset(self)
|
||||
self._conn = {}
|
||||
self._bestconn = {}
|
||||
|
||||
|
||||
def create_connection(self, id1, id2, pos1, pos2, cost1, cost2):
|
||||
# Process data
|
||||
hash = min(id1, id2), max(id1, id2)
|
||||
val = min(pos1, pos2), max(pos1, pos2)
|
||||
cost = min(cost1, cost2)
|
||||
# Add to total list
|
||||
self._conn.setdefault(hash, []).append(val)
|
||||
# Keep track of connection with lowest cost
|
||||
curcost = self._bestconn.get(hash, (np.inf,))[0]
|
||||
if cost < curcost:
|
||||
self._bestconn[hash] = (cost,) + val
|
||||
|
||||
|
||||
def test_connections():
|
||||
|
||||
# Create MCP object with three seed points
|
||||
mcp = MCP(a)
|
||||
costs, traceback = mcp.find_costs([ (1,1), (7,7), (1,7) ])
|
||||
|
||||
# Test that all three seed points are connected
|
||||
connections = set(mcp._conn.keys())
|
||||
assert (0, 1) in connections
|
||||
assert (1, 2) in connections
|
||||
assert (0, 2) in connections
|
||||
|
||||
# Test that any two neighbors have only been connected once
|
||||
for position_tuples in mcp._conn.values():
|
||||
n1 = len(position_tuples)
|
||||
n2 = len(set(position_tuples))
|
||||
assert n1 == n2
|
||||
|
||||
# For seed 0 and 1
|
||||
cost, pos1, pos2 = mcp._bestconn[(0,1)]
|
||||
# Test meeting points
|
||||
assert (pos1, pos2) == ( (3,3), (4,4) )
|
||||
# Test the whole path
|
||||
path = mcp.traceback(pos1) + list(reversed(mcp.traceback(pos2)))
|
||||
assert_array_equal(path,
|
||||
[(1, 1), (2, 2), (3, 3), (4, 4), (5, 5), (6, 6), (7, 7)])
|
||||
|
||||
# For seed 1 and 2
|
||||
cost, pos1, pos2 = mcp._bestconn[(1,2)]
|
||||
# Test meeting points
|
||||
assert (pos1, pos2) == ( (3,7), (4,7) )
|
||||
# Test the whole path
|
||||
path = mcp.traceback(pos1) + list(reversed(mcp.traceback(pos2)))
|
||||
assert_array_equal(path,
|
||||
[(1, 7), (2, 7), (3, 7), (4, 7), (5, 7), (6, 7), (7, 7)])
|
||||
|
||||
# For seed 0 and 2
|
||||
cost, pos1, pos2 = mcp._bestconn[(0,2)]
|
||||
# Test meeting points
|
||||
assert (pos1, pos2) == ( (1,3), (1,4) )
|
||||
# Test the whole path
|
||||
path = mcp.traceback(pos1) + list(reversed(mcp.traceback(pos2)))
|
||||
assert_array_equal(path,
|
||||
[(1, 1), (1, 2), (1, 3), (1, 4), (1, 5), (1, 6), (1, 7)])
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
np.testing.run_module_suite()
|
||||
@@ -0,0 +1,61 @@
|
||||
import skimage.graph.mcp as mcp
|
||||
from numpy.testing import (assert_array_equal,
|
||||
assert_almost_equal,
|
||||
)
|
||||
|
||||
import numpy as np
|
||||
|
||||
a = np.ones((8, 8), dtype=np.float32)
|
||||
a[1::2] *= 2.0
|
||||
|
||||
|
||||
class TestFlexibleMCP(mcp.MCP_Flexible):
|
||||
""" Simple MCP subclass that allows the front to travel
|
||||
a certain distance from the seed point, and uses a constant
|
||||
cost factor that is independant of the cost array.
|
||||
"""
|
||||
|
||||
def _reset(self):
|
||||
mcp.MCP_Flexible._reset(self)
|
||||
self._distance = np.zeros((8, 8), dtype=np.float32).ravel()
|
||||
|
||||
def goal_reached(self, index, cumcost):
|
||||
if self._distance[index] > 4:
|
||||
return 2
|
||||
else:
|
||||
return 0
|
||||
|
||||
def travel_cost(self, index, new_index, offset_length):
|
||||
return 1.0 # fixed cost
|
||||
|
||||
def examine_neighbor(self, index, new_index, offset_length):
|
||||
pass # We do not test this
|
||||
|
||||
def update_node(self, index, new_index, offset_length):
|
||||
self._distance[new_index] = self._distance[index] + 1
|
||||
|
||||
|
||||
def test_flexible():
|
||||
|
||||
# Create MCP and do a traceback
|
||||
mcp = TestFlexibleMCP(a)
|
||||
costs, traceback = mcp.find_costs([(0, 0)])
|
||||
|
||||
# Check that inner part is correct. This basically
|
||||
# tests whether travel_cost works.
|
||||
assert_array_equal(costs[:4,:4], [[1, 2, 3, 4],
|
||||
[2, 2, 3, 4],
|
||||
[3, 3, 3, 4],
|
||||
[4, 4, 4, 4]])
|
||||
|
||||
# Test that the algorithm stopped at the right distance.
|
||||
# Note that some of the costs are filled in but not yet frozen,
|
||||
# so we take a bit of margin
|
||||
assert np.all( costs[-2:,:] == np.inf )
|
||||
assert np.all( costs[:,-2:] == np.inf )
|
||||
|
||||
#print(costs)
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
np.testing.run_module_suite()
|
||||
@@ -116,8 +116,8 @@ def test_offsets():
|
||||
m = mcp.MCP(a, offsets=offsets)
|
||||
costs, traceback = m.find_costs([(1, 6)])
|
||||
assert_array_equal(traceback,
|
||||
[[-1, -1, -1, -1, -1, -1, -1, -1],
|
||||
[-1, -1, -1, -1, -1, -1, -1, -1],
|
||||
[[-2, -2, -2, -2, -2, -2, -2, -2],
|
||||
[-2, -2, -2, -2, -2, -2, -1, -2],
|
||||
[15, 14, 13, 12, 11, 10, 0, 1],
|
||||
[10, 0, 1, 2, 3, 4, 5, 6],
|
||||
[10, 0, 1, 2, 3, 4, 5, 6],
|
||||
@@ -151,5 +151,6 @@ def _test_random(shape):
|
||||
return a, costs, offsets
|
||||
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
np.testing.run_module_suite()
|
||||
|
||||
Reference in New Issue
Block a user