diff --git a/skimage/graph/__init__.py b/skimage/graph/__init__.py index eb817c77..a335971d 100644 --- a/skimage/graph/__init__.py +++ b/skimage/graph/__init__.py @@ -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'] \ No newline at end of file diff --git a/skimage/graph/_mcp.pxd b/skimage/graph/_mcp.pxd index b2e2a548..1cfdbc9a 100644 --- a/skimage/graph/_mcp.pxd +++ b/skimage/graph/_mcp.pxd @@ -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) + \ No newline at end of file diff --git a/skimage/graph/_mcp.pyx b/skimage/graph/_mcp.pyx index 34904dcc..c65922db 100644 --- a/skimage/graph/_mcp.pyx +++ b/skimage/graph/_mcp.pyx @@ -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) diff --git a/skimage/graph/mcp.py b/skimage/graph/mcp.py index dc584226..bf693a45 100644 --- a/skimage/graph/mcp.py +++ b/skimage/graph/mcp.py @@ -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, diff --git a/skimage/graph/tests/test_anisotropy.py b/skimage/graph/tests/test_anisotropy.py new file mode 100644 index 00000000..199d2e73 --- /dev/null +++ b/skimage/graph/tests/test_anisotropy.py @@ -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() diff --git a/skimage/graph/tests/test_connect.py b/skimage/graph/tests/test_connect.py new file mode 100644 index 00000000..619ac3b4 --- /dev/null +++ b/skimage/graph/tests/test_connect.py @@ -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() \ No newline at end of file diff --git a/skimage/graph/tests/test_flexible.py b/skimage/graph/tests/test_flexible.py new file mode 100644 index 00000000..9ae50759 --- /dev/null +++ b/skimage/graph/tests/test_flexible.py @@ -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() diff --git a/skimage/graph/tests/test_mcp.py b/skimage/graph/tests/test_mcp.py index 560f19d0..5c44abad 100644 --- a/skimage/graph/tests/test_mcp.py +++ b/skimage/graph/tests/test_mcp.py @@ -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()