From 4b2136bc7d17426458a518aded70120a408bf476 Mon Sep 17 00:00:00 2001 From: Almar Klein Date: Fri, 20 Sep 2013 11:25:14 +0200 Subject: [PATCH 01/16] Make mcp evolution algorithm a bit more readable. In particular, distinguish between cost and cumcost, the same variable name was used for this. --- skimage/graph/_mcp.pyx | 33 +++++++++++++++++++-------------- 1 file changed, 19 insertions(+), 14 deletions(-) diff --git a/skimage/graph/_mcp.pyx b/skimage/graph/_mcp.pyx index 34904dcc..42985a37 100644 --- a/skimage/graph/_mcp.pyx +++ b/skimage/graph/_mcp.pyx @@ -434,7 +434,7 @@ cdef class MCP: else: costs_heap.push_fast(0, start) - cdef FLOAT_T cost, new_cost + cdef FLOAT_T cost, new_cost, cumcost, new_cumcost cdef INDEX_T index, new_index cdef BOOL_T is_at_edge, use_offset cdef INDEX_T d, i @@ -442,7 +442,7 @@ cdef class MCP: 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: # Find the point with the minimum cost in the heap. Once # popped, this point's minimum cost path has been found. @@ -450,12 +450,13 @@ cdef class MCP: # 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 if use_ends: # If we're only tracing out a path to one or more @@ -512,27 +513,31 @@ cdef class MCP: # ignore this point if flat_cumulative_costs[new_index] != inf: 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_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], - 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 - + # Un-flatten the costs and traceback arrays for human consumption. cumulative_costs = flat_cumulative_costs.reshape(self.costs_shape, order='F') From 8d09300155385e16ac63b8c118b3b7e52821a03b Mon Sep 17 00:00:00 2001 From: Almar Klein Date: Fri, 20 Sep 2013 16:22:10 +0200 Subject: [PATCH 02/16] Improve speed of MPC by about 30% Dear diary, I spend a whole day searching for why this implementation of the MPC was slower then my own implementation. I checked *everything*. I rewrote the front evolution to match with mine almost exactly. I was about to give up. And then I saw this... Arg! --- skimage/graph/_mcp.pyx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/skimage/graph/_mcp.pyx b/skimage/graph/_mcp.pyx index 42985a37..23d433f9 100644 --- a/skimage/graph/_mcp.pyx +++ b/skimage/graph/_mcp.pyx @@ -290,7 +290,7 @@ cdef class MCP: 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 From f5c675b8c6e9d269f63e8a594861c28b654c6607 Mon Sep 17 00:00:00 2001 From: Almar Klein Date: Mon, 14 Oct 2013 15:41:15 +0200 Subject: [PATCH 03/16] MCP: turn while loop into for loop to prevent infinite looping. A nice feature during development :) --- skimage/graph/_mcp.pyx | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/skimage/graph/_mcp.pyx b/skimage/graph/_mcp.pyx index 23d433f9..15700928 100644 --- a/skimage/graph/_mcp.pyx +++ b/skimage/graph/_mcp.pyx @@ -437,13 +437,17 @@ cdef class MCP: cdef FLOAT_T cost, new_cost, cumcost, new_cumcost 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 - while 1: + + for iter in range(flat_costs.size): + # 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: From d1329488e6f039c46e9f9ee6aae38c698861b197 Mon Sep 17 00:00:00 2001 From: Almar Klein Date: Mon, 14 Oct 2013 17:05:35 +0200 Subject: [PATCH 04/16] MCP: add sampling attribute: deal with anisotropic data. --- skimage/graph/_mcp.pyx | 38 +++++++++++++++++++++++++++----------- 1 file changed, 27 insertions(+), 11 deletions(-) diff --git a/skimage/graph/_mcp.pyx b/skimage/graph/_mcp.pyx index 15700928..47596b3c 100644 --- a/skimage/graph/_mcp.pyx +++ b/skimage/graph/_mcp.pyx @@ -227,7 +227,7 @@ def _reverse(arr): @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,26 @@ 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 @@ -326,7 +340,7 @@ 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) + np.sum((sampling*self.offsets)**2, axis=1)).astype(FLOAT_D) self.dirty = 0 self.use_start_cost = 1 @@ -595,7 +609,7 @@ cdef class MCP: 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 - + cdef OFFSETS_INDEX_T offset cdef DIM_T d cdef DIM_T dim = self.dim @@ -635,15 +649,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 From 280fadff39239c0597a909782b197b318c71a88a Mon Sep 17 00:00:00 2001 From: Almar Klein Date: Thu, 24 Oct 2013 16:13:35 +0200 Subject: [PATCH 05/16] MCP: refactored to use typed memoryviews The arrays are stored on the mcp object, so that they can be used in methods as well. Will use that to make MCP more flexibel. --- skimage/graph/_mcp.pxd | 26 +++++++++----- skimage/graph/_mcp.pyx | 82 +++++++++++++++++++++--------------------- 2 files changed, 59 insertions(+), 49 deletions(-) diff --git a/skimage/graph/_mcp.pxd b/skimage/graph/_mcp.pxd index b2e2a548..4016db1c 100644 --- a/skimage/graph/_mcp.pxd +++ b/skimage/graph/_mcp.pxd @@ -9,22 +9,30 @@ 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 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 cdef FLOAT_T _travel_cost(self, FLOAT_T old_cost, FLOAT_T new_cost, FLOAT_T offset_length) diff --git a/skimage/graph/_mcp.pyx b/skimage/graph/_mcp.pyx index 47596b3c..27a2da01 100644 --- a/skimage/graph/_mcp.pyx +++ b/skimage/graph/_mcp.pyx @@ -39,13 +39,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 @@ -274,6 +270,7 @@ cdef class MCP: returned by the find_costs() method. """ + def __init__(self, costs, offsets=None, fully_connected=True, sampling=None): """__init__(costs, offsets=None, fully_connected=True, sampling=None) @@ -301,7 +298,7 @@ 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.flat_cumulative_costs[...] = np.inf self.dim = len(costs.shape) self.costs_shape = costs.shape self.costs_heap = heap.FastUpdateBinaryHeap(initial_capacity=128, @@ -311,7 +308,7 @@ cdef class MCP: # 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) + self.traceback_offsets[...] = -1 # The offsets are a list of relative offsets from a central # point to each point in the relevant neighborhood. (e.g. (-1, @@ -343,21 +340,24 @@ cdef class MCP: 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[...] = -1 + self.flat_cumulative_costs[...] = np.inf self.dirty = 0 - + + cdef FLOAT_T _travel_cost(self, FLOAT_T old_cost, FLOAT_T new_cost, FLOAT_T offset_length): return new_cost - + + def find_costs(self, starts, ends=None, find_all_ends=True): """ Find the minimum-cost path from the given starting points. @@ -406,7 +406,7 @@ 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') @@ -422,32 +422,31 @@ cdef class MCP: if self.dirty: self._reset() - - # lookup and array-ify object attributes for fast use + + # 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! + cdef INDEX_T start 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) - + + # Variables used during front propagation cdef FLOAT_T cost, new_cost, cumcost, new_cumcost cdef INDEX_T index, new_index cdef BOOL_T is_at_edge, use_offset @@ -475,7 +474,7 @@ cdef class MCP: # Record the cost we found to this point flat_cumulative_costs[index] = cumcost - + 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 @@ -489,7 +488,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 @@ -557,12 +556,15 @@ cdef class MCP: traceback_offsets[new_index] = i # Un-flatten the costs and traceback arrays for human consumption. - cumulative_costs = flat_cumulative_costs.reshape(self.costs_shape, + cumulative_costs = np.asarray(flat_cumulative_costs) + cumulative_costs = cumulative_costs.reshape(self.costs_shape, order='F') - traceback = traceback_offsets.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) @@ -602,13 +604,13 @@ 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 From 798bea72514b542561c6890c47bedc959bf92e66 Mon Sep 17 00:00:00 2001 From: Almar Klein Date: Fri, 13 Dec 2013 13:52:52 +0100 Subject: [PATCH 06/16] update-mcp: Refactor MCP algorithm in prep for MCP subclasses Installed a few hooks to allow more flexibility in subclasses, refactored traceback() to allow obtaining tracebacks fast during front evolution, and moved some initializations to reset(). Full list: * integer division in calculating indices * traceback_offsets use -2 for uninitialized, -1 for seed point. In that we we can distinguish between these two cases in traceback() * Clearing of traceback_offset, cumulative_costs, pushing start positions into the heap, are now all done in reset() * Reset is called every time that find_cost is called. * Add hook _goal_reached() * Add hook _examine_neighbor() * Add hook _update_node() * Replace while loop with a for-loop to prevent infinite looping in case of a bug. * Split traceback() in two functions to allow getting a flat traceback fast. traceback() calles these two private functions and behaves as just before. --- skimage/graph/_mcp.pxd | 7 ++ skimage/graph/_mcp.pyx | 186 ++++++++++++++++++++++++++++++----------- 2 files changed, 146 insertions(+), 47 deletions(-) diff --git a/skimage/graph/_mcp.pxd b/skimage/graph/_mcp.pxd index 4016db1c..3b35ce0b 100644 --- a/skimage/graph/_mcp.pxd +++ b/skimage/graph/_mcp.pxd @@ -35,4 +35,11 @@ cdef class MCP: 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) + + cdef object _flat_traceback(self, INDEX_T end) + cdef object _unravel_traceback(self, object flat_traceback) + diff --git a/skimage/graph/_mcp.pyx b/skimage/graph/_mcp.pyx index 27a2da01..b3f441e3 100644 --- a/skimage/graph/_mcp.pyx +++ b/skimage/graph/_mcp.pyx @@ -172,7 +172,7 @@ def _unravel_index_fortran(flat_indices, shape): """ 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 @@ -298,7 +298,6 @@ 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[...] = np.inf self.dim = len(costs.shape) self.costs_shape = costs.shape self.costs_heap = heap.FastUpdateBinaryHeap(initial_capacity=128, @@ -308,7 +307,6 @@ cdef class MCP: # 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[...] = -1 # The offsets are a list of relative offsets from a central # point to each point in the relevant neighborhood. (e.g. (-1, @@ -348,17 +346,60 @@ cdef class MCP: Clears paths found by find_costs(). """ self.costs_heap.reset() - self.traceback_offsets[...] = -1 + 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): + """ 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(self, int index, float cumcost) + This method is called each iteration after popping an index + from the heap, before examining the neighbours. + + This method should return 1 if the algorithm should not check + the current point's neighbours and 2 if the algorithm is now + done, for example an end point is reached. + """ + return 0 + + + cdef void _examine_neighbor(self, INDEX_T index, INDEX_T new_index, FLOAT_T offset_length): + """ _examine_neighbor(self, int index, int new_index, float offset_length) + This method is called for every neighbor examined, even before + checking whether it is frozen. + """ + pass + + + cdef void _update_node(self, INDEX_T index, INDEX_T new_index, FLOAT_T offset_length): + """ _update_node(self, 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. @@ -410,6 +451,8 @@ cdef class MCP: 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: @@ -420,8 +463,10 @@ cdef class MCP: flat_ends = np.array(_ravel_index_fortran( ends, self.costs_shape), dtype=INDEX_D) - if self.dirty: - self._reset() + # 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 @@ -437,17 +482,9 @@ cdef class MCP: cdef heap.FastUpdateBinaryHeap costs_heap = self.costs_heap 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! - cdef INDEX_T start - 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) # Variables used during front propagation - cdef FLOAT_T cost, new_cost, cumcost, new_cumcost + 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, iter @@ -455,9 +492,12 @@ cdef class MCP: cdef EDGE_T pos_edge_val, neg_edge_val cdef int num_ends_found = 0 cdef FLOAT_T inf = np.inf + cdef int goal_reached + cdef INDEX_T maxiter = int(max_coverage * flat_costs.size) - for iter in range(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. @@ -475,6 +515,14 @@ cdef class MCP: # Record the cost we found to this point 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 @@ -522,10 +570,14 @@ 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] - + + # Allow subclasses to examine this neighbour + offset_length = offset_lengths[i] + self._examine_neighbor(index, new_index, offset_length) + # If we have already found the best path here then # ignore this point if flat_cumulative_costs[new_index] != inf: @@ -541,7 +593,7 @@ cdef class MCP: # Calculate new cumulative cost new_cumcost = cumcost + self._travel_cost(cost, new_cost, - offset_lengths[i]) + 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 @@ -554,6 +606,8 @@ cdef class MCP: # 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 = np.asarray(flat_cumulative_costs) @@ -565,6 +619,64 @@ cdef class MCP: return cumulative_costs, traceback + cdef object _flat_traceback(self, INDEX_T end): + """ _flat_traceback(end) + Do a traceback, input and output is in flat coordinates. + Returns a list of integers, from given end point to start. + """ + + # Initialize traceback + traceback = [end] + + # Init position with end + cdef INDEX_T flat_position = end + # Check if we can find a path + #print(flat_position) + if self.traceback_offsets[flat_position] == -2: + raise ValueError('no minimum-cost path was found ' + 'to the specified end point') + + # Short names for arrays + cdef OFFSETS_INDEX_T [:] traceback_offsets = self.traceback_offsets + cdef INDEX_T [:] flat_offsets = self.flat_offsets + + + # Do traceback + cdef OFFSETS_INDEX_T offset + while 1: + offset = traceback_offsets[flat_position] + if offset == -1: + break # -2 is uninitialized, -1 is start point + flat_position -= flat_offsets[offset] + traceback.append(offset) + return traceback + + + cdef object _unravel_traceback(self, object flat_traceback): + """ Unravel the given traceback obtained from _flat_traceback(). + Returns a new traceback that is reversed and has x-y(-z) coordinates. + """ + + flat_end = flat_traceback.pop(0) + end = _unravel_index_fortran([flat_end], self.costs_shape)[0] + + # Prepare arrays + cdef INDEX_T [:] position = np.array(end, dtype=INDEX_D) + cdef OFFSET_T [:,:] offsets = self.offsets + + # Prepare variables + cdef DIM_T dim = self.dim + cdef DIM_T d + + # Reverse and convert + traceback = [tuple(position)] + for offset in flat_traceback: + for d in range(dim): + position[d] -= offsets[offset, d] + traceback.append(tuple(position)) + return _reverse(traceback) + + def traceback(self, end): """traceback(end) @@ -597,34 +709,14 @@ cdef class MCP: if ends is None: raise ValueError('the specified end point must be ' 'within the costs array') - traceback = [tuple(ends[0])] + + # Get flat traceback + cdef INDEX_T flat_end = _ravel_index_fortran(ends, self.costs_shape)[0] + flat_traceback = self._flat_traceback(flat_end) + + # Transform it + return self._unravel_traceback(flat_traceback) - cdef INDEX_T 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') - - # 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 - while 1: - offset = traceback_offsets[flat_position] - if offset == -1: - # At a point where we can go no further: probably a start point - break - flat_position -= flat_offsets[offset] - for d in range(dim): - position[d] -= offsets[offset, d] - traceback.append(tuple(position)) - return _reverse(traceback) @cython.boundscheck(False) From 8fe7416d41a3c93f06e252cd39400707f14474f6 Mon Sep 17 00:00:00 2001 From: Almar Klein Date: Fri, 13 Dec 2013 13:58:44 +0100 Subject: [PATCH 07/16] update-mcp: Add MCP_Connect class This class finds the connections between the given seed points. This is a bit similar to finding connections between start and end points, except now start- and endpoints are the same. --- skimage/graph/_mcp.pyx | 124 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 124 insertions(+) diff --git a/skimage/graph/_mcp.pyx b/skimage/graph/_mcp.pyx index b3f441e3..3156d0a6 100644 --- a/skimage/graph/_mcp.pyx +++ b/skimage/graph/_mcp.pyx @@ -761,3 +761,127 @@ 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 A and + B 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! + # We need to establish the traceback now, because the traceback + # may be overridden in future iterations. + path1 = self._flat_traceback(index) + path2 = self._flat_traceback(new_index) + # 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, path1, path2, cost1, cost2) + + + def unravel_traceback(self, flat_traceback): + """ unravel_traceback(flat_traceback) + + This method can be used to unravel a flat traceback that is passed + to create_connection(). The result is equivalent to a traceback + returned by traceback(). + """ + return self._unravel_traceback(flat_traceback) + + + def create_connection(self, id1, id2, tb1, tb2, cost1, cost2): + """ create_connection id1, id2, traceback1, traceback2, 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 costs). + + Parameters + ---------- + id1 : int + The index of the seed point where front A originated from. + id2 : int + The index of the seed point where front B originated from. + traceback1 : list + A list of integers representing the flat traceback from the + meeting point to seed point id1. Use the unravel_traceback() method + to turn obtain a list of coordinates from source to meeting point. + traceback2 : list + Dito for id2. + cost1 : float + The cumulative cost at side A of the connection. + cost2 : float + The cumulative costs at side B of the connection. + + """ + 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] + From 32c4a80b32d16c6741c2c89f9cfbe7271bfa0dfa Mon Sep 17 00:00:00 2001 From: Almar Klein Date: Fri, 13 Dec 2013 14:00:33 +0100 Subject: [PATCH 08/16] update-mcp: Add MCP_Flexible subclass. This class adds a public def-method for the private cdef-method hooks. This allows users to influence the behavior of the MCP algorithm by subclassing it in pure Python. --- skimage/graph/_mcp.pyx | 58 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 58 insertions(+) diff --git a/skimage/graph/_mcp.pyx b/skimage/graph/_mcp.pyx index 3156d0a6..98c71bdb 100644 --- a/skimage/graph/_mcp.pyx +++ b/skimage/graph/_mcp.pyx @@ -885,3 +885,61 @@ cdef class Mcp_Connect(MCP): """ 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(self, old_cost, new_cost, offset_length) + The travel cost for going from the current node to the next. + 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(self, index, new_index, offset_length) + This method is called for every neighbor examined, even before + checking whether it is frozen. + Overload this method to adapt the behaviour of the algorithm. + """ + pass + + + def update_node(self, INDEX_T index, INDEX_T new_index, + FLOAT_T offset_length): + """ update_node(self, index, new_index, offset_length) + This method is called when a node is updated. + Overload this method to adapt the behaviour of the algorithm. + """ + 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) From f7416507faf7e2c30325daaf72e5471bd6f135c2 Mon Sep 17 00:00:00 2001 From: Almar Klein Date: Fri, 13 Dec 2013 16:37:36 +0100 Subject: [PATCH 09/16] update-mcp: Small fix. Need _start and _end attributes for the MCP class. It worked for my own tests because I overloaded the MCP class. --- skimage/graph/_mcp.pxd | 2 ++ 1 file changed, 2 insertions(+) diff --git a/skimage/graph/_mcp.pxd b/skimage/graph/_mcp.pxd index 3b35ce0b..a4025e06 100644 --- a/skimage/graph/_mcp.pxd +++ b/skimage/graph/_mcp.pxd @@ -18,6 +18,8 @@ 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 BOOL_T dirty cdef BOOL_T use_start_cost From 2ebeee2df466c0829ab412deb9733b2584d16ac0 Mon Sep 17 00:00:00 2001 From: Almar Klein Date: Sat, 14 Dec 2013 23:55:48 +0100 Subject: [PATCH 10/16] update-mcp: Change when examine_neighbors is called. Also renamed Mcp_Connect to MCP_Connect. Makes more sense to call it when the neighbor is frozen. In that way the method is called exactly once for each pair of neighbours. Plus the costs and paths are known and fixed. --- skimage/graph/_mcp.pyx | 70 +++++++++++++++++++----------------------- 1 file changed, 32 insertions(+), 38 deletions(-) diff --git a/skimage/graph/_mcp.pyx b/skimage/graph/_mcp.pyx index 98c71bdb..3b1ceb08 100644 --- a/skimage/graph/_mcp.pyx +++ b/skimage/graph/_mcp.pyx @@ -385,8 +385,8 @@ cdef class MCP: cdef void _examine_neighbor(self, INDEX_T index, INDEX_T new_index, FLOAT_T offset_length): """ _examine_neighbor(self, int index, int new_index, float offset_length) - This method is called for every neighbor examined, even before - checking whether it is frozen. + This method is called once for every pair of neighboring nodes, + as soon as both nodes become frozen. """ pass @@ -574,13 +574,17 @@ cdef class MCP: # using the flat offsets, calculate the new flat index new_index = index + flat_offsets[i] - # Allow subclasses to examine this neighbour + # Get offset length offset_length = offset_lengths[i] - self._examine_neighbor(index, new_index, offset_length) # 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 # Get cost and new cost @@ -766,15 +770,15 @@ cdef class MCP_Geometric(MCP): @cython.boundscheck(True) @cython.wraparound(True) -cdef class Mcp_Connect(MCP): - """Mcp_Connect(costs, offsets=None, fully_connected=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 A and - B meet, create_connection() is called. This method must be overloaded - to deal with the found edges in a way that is appropriate for 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. """ @@ -828,51 +832,41 @@ cdef class Mcp_Connect(MCP): if id2 < 0 or id1 < 0: pass elif id2 != id1: - # we reached the 'front' of another seed point! - # We need to establish the traceback now, because the traceback - # may be overridden in future iterations. - path1 = self._flat_traceback(index) - path2 = self._flat_traceback(new_index) + # 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, path1, path2, cost1, cost2) - - - def unravel_traceback(self, flat_traceback): - """ unravel_traceback(flat_traceback) - - This method can be used to unravel a flat traceback that is passed - to create_connection(). The result is equivalent to a traceback - returned by traceback(). - """ - return self._unravel_traceback(flat_traceback) + self.create_connection(id1, id2, pos1, pos2, cost1, cost2) def create_connection(self, id1, id2, tb1, tb2, cost1, cost2): - """ create_connection id1, id2, traceback1, traceback2, 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 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 index of the seed point where front A originated from. + The seed point id where the first neighbor originated from. id2 : int - The index of the seed point where front B originated from. - traceback1 : list - A list of integers representing the flat traceback from the - meeting point to seed point id1. Use the unravel_traceback() method - to turn obtain a list of coordinates from source to meeting point. - traceback2 : list - Dito for id2. + 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 side A of the connection. + The cumulative cost at pos1. cost2 : float - The cumulative costs at side B of the connection. + The cumulative costs at pos2. """ pass @@ -914,8 +908,8 @@ cdef class MCP_Flexible(MCP): def examine_neighbor(self, INDEX_T index, INDEX_T new_index, FLOAT_T offset_length): """ examine_neighbor(self, index, new_index, offset_length) - This method is called for every neighbor examined, even before - checking whether it is frozen. + This method is called once for every pair of neighboring nodes, + as soon as both nodes become frozen. Overload this method to adapt the behaviour of the algorithm. """ pass From fdd880aa7f8e3595dc281cf6aa0b4c3cbfdfaa8b Mon Sep 17 00:00:00 2001 From: Almar Klein Date: Sat, 14 Dec 2013 23:58:45 +0100 Subject: [PATCH 11/16] update-mcp: add new MCP classes to skimage.graph namespace --- skimage/graph/__init__.py | 4 +++- skimage/graph/mcp.py | 2 +- 2 files changed, 4 insertions(+), 2 deletions(-) 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.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, From 789ade3ca29a34c871d99f40584bf1fd3cc23000 Mon Sep 17 00:00:00 2001 From: Almar Klein Date: Sun, 15 Dec 2013 00:37:27 +0100 Subject: [PATCH 12/16] update-mcp: added tests for newly added features. Also fix in existing test related to how the traceback arrays is initialized: first it was initialized with all -1, now with -2, and -1 at the seed points. --- skimage/graph/tests/test_anisotropy.py | 53 ++++++++++++++++ skimage/graph/tests/test_connect.py | 83 ++++++++++++++++++++++++++ skimage/graph/tests/test_mcp.py | 5 +- 3 files changed, 139 insertions(+), 2 deletions(-) create mode 100644 skimage/graph/tests/test_anisotropy.py create mode 100644 skimage/graph/tests/test_connect.py diff --git a/skimage/graph/tests/test_anisotropy.py b/skimage/graph/tests/test_anisotropy.py new file mode 100644 index 00000000..a2d3a240 --- /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_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() From fc90a164fe1ca95d6f8d574b6541b321e7b01db1 Mon Sep 17 00:00:00 2001 From: Almar Klein Date: Sun, 15 Dec 2013 12:44:53 +0100 Subject: [PATCH 13/16] update-mcp: improvements to docstrings. --- skimage/graph/_mcp.pyx | 35 ++++++++++++++++++++++++++--------- 1 file changed, 26 insertions(+), 9 deletions(-) diff --git a/skimage/graph/_mcp.pyx b/skimage/graph/_mcp.pyx index 3b1ceb08..eee03000 100644 --- a/skimage/graph/_mcp.pyx +++ b/skimage/graph/_mcp.pyx @@ -376,9 +376,14 @@ cdef class MCP: 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, for example an end point is reached. + done. """ return 0 @@ -847,8 +852,9 @@ cdef class MCP_Connect(MCP): """ 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 costs). + 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. @@ -899,8 +905,10 @@ cdef class MCP_Flexible(MCP): def travel_cost(self, FLOAT_T old_cost, FLOAT_T new_cost, FLOAT_T offset_length): """ travel_cost(self, old_cost, new_cost, offset_length) - The travel cost for going from the current node to the next. - Overload this method to adapt the behaviour of the algorithm. + 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 @@ -909,8 +917,12 @@ cdef class MCP_Flexible(MCP): FLOAT_T offset_length): """ examine_neighbor(self, index, new_index, offset_length) This method is called once for every pair of neighboring nodes, - as soon as both nodes become frozen. - Overload this method to adapt the behaviour of the algorithm. + 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 @@ -918,8 +930,13 @@ cdef class MCP_Flexible(MCP): def update_node(self, INDEX_T index, INDEX_T new_index, FLOAT_T offset_length): """ update_node(self, index, new_index, offset_length) - This method is called when a node is updated. - Overload this method to adapt the behaviour of the algorithm. + 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 From 9113d5a08b9eed65dbc1bb4b0d745ad57b32d749 Mon Sep 17 00:00:00 2001 From: Almar Klein Date: Sun, 15 Dec 2013 13:32:15 +0100 Subject: [PATCH 14/16] update-mcp: add tests for MCP_Flexible --- skimage/graph/tests/test_flexible.py | 61 ++++++++++++++++++++++++++++ 1 file changed, 61 insertions(+) create mode 100644 skimage/graph/tests/test_flexible.py diff --git a/skimage/graph/tests/test_flexible.py b/skimage/graph/tests/test_flexible.py new file mode 100644 index 00000000..3570a31f --- /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() From 3cb5a194798d722c240c20c58de7cb16cf9ec7a9 Mon Sep 17 00:00:00 2001 From: Almar Klein Date: Fri, 3 Jan 2014 15:09:06 +0100 Subject: [PATCH 15/16] Update MCP: Put back the traceback() method. I initially changed it to allow faster (flat) tracebacks *during* front evolution, but it turned out not to be necessary. Putting back to minimize the changes of this PR. --- skimage/graph/_mcp.pxd | 5 +-- skimage/graph/_mcp.pyx | 97 +++++++++++++----------------------------- 2 files changed, 31 insertions(+), 71 deletions(-) diff --git a/skimage/graph/_mcp.pxd b/skimage/graph/_mcp.pxd index a4025e06..1cfdbc9a 100644 --- a/skimage/graph/_mcp.pxd +++ b/skimage/graph/_mcp.pxd @@ -41,7 +41,4 @@ cdef class MCP: 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) - - cdef object _flat_traceback(self, INDEX_T end) - cdef object _unravel_traceback(self, object flat_traceback) - + \ No newline at end of file diff --git a/skimage/graph/_mcp.pyx b/skimage/graph/_mcp.pyx index eee03000..ed50fe95 100644 --- a/skimage/graph/_mcp.pyx +++ b/skimage/graph/_mcp.pyx @@ -628,80 +628,22 @@ cdef class MCP: return cumulative_costs, traceback - cdef object _flat_traceback(self, INDEX_T end): - """ _flat_traceback(end) - Do a traceback, input and output is in flat coordinates. - Returns a list of integers, from given end point to start. - """ - - # Initialize traceback - traceback = [end] - - # Init position with end - cdef INDEX_T flat_position = end - # Check if we can find a path - #print(flat_position) - if self.traceback_offsets[flat_position] == -2: - raise ValueError('no minimum-cost path was found ' - 'to the specified end point') - - # Short names for arrays - cdef OFFSETS_INDEX_T [:] traceback_offsets = self.traceback_offsets - cdef INDEX_T [:] flat_offsets = self.flat_offsets - - - # Do traceback - cdef OFFSETS_INDEX_T offset - while 1: - offset = traceback_offsets[flat_position] - if offset == -1: - break # -2 is uninitialized, -1 is start point - flat_position -= flat_offsets[offset] - traceback.append(offset) - return traceback - - - cdef object _unravel_traceback(self, object flat_traceback): - """ Unravel the given traceback obtained from _flat_traceback(). - Returns a new traceback that is reversed and has x-y(-z) coordinates. - """ - - flat_end = flat_traceback.pop(0) - end = _unravel_index_fortran([flat_end], self.costs_shape)[0] - - # Prepare arrays - cdef INDEX_T [:] position = np.array(end, dtype=INDEX_D) - cdef OFFSET_T [:,:] offsets = self.offsets - - # Prepare variables - cdef DIM_T dim = self.dim - cdef DIM_T d - - # Reverse and convert - traceback = [tuple(position)] - for offset in flat_traceback: - for d in range(dim): - position[d] -= offsets[offset, d] - traceback.append(tuple(position)) - return _reverse(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 @@ -718,13 +660,34 @@ cdef class MCP: if ends is None: raise ValueError('the specified end point must be ' 'within the costs array') + traceback = [tuple(ends[0])] + + cdef INDEX_T 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') - # Get flat traceback - cdef INDEX_T flat_end = _ravel_index_fortran(ends, self.costs_shape)[0] - flat_traceback = self._flat_traceback(flat_end) + # 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) - # Transform it - return self._unravel_traceback(flat_traceback) + cdef OFFSETS_INDEX_T offset + cdef DIM_T d + cdef DIM_T dim = self.dim + while 1: + offset = traceback_offsets[flat_position] + if offset == -1: + # At a point where we can go no further: probably a start point + break + flat_position -= flat_offsets[offset] + for d in range(dim): + position[d] -= offsets[offset, d] + traceback.append(tuple(position)) + return _reverse(traceback) From 14288e1f1e24ba8ffd8f997bf67989d37a270f54 Mon Sep 17 00:00:00 2001 From: Almar Klein Date: Mon, 6 Jan 2014 12:39:13 +0100 Subject: [PATCH 16/16] PEP8 stuff. --- skimage/graph/_mcp.pyx | 116 +++++++++++++------------ skimage/graph/tests/test_anisotropy.py | 26 +++--- skimage/graph/tests/test_flexible.py | 10 +-- 3 files changed, 79 insertions(+), 73 deletions(-) diff --git a/skimage/graph/_mcp.pyx b/skimage/graph/_mcp.pyx index ed50fe95..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 @@ -102,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 @@ -168,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 @@ -181,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])) @@ -205,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) @@ -216,7 +219,8 @@ 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] @@ -272,7 +276,7 @@ cdef class MCP: """ def __init__(self, costs, offsets=None, fully_connected=True, - sampling=None): + sampling=None): """__init__(costs, offsets=None, fully_connected=True, sampling=None) See class documentation. @@ -311,8 +315,8 @@ cdef class MCP: # 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) @@ -322,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. @@ -334,15 +338,14 @@ cdef class MCP: # The offset lengths are the distances traveled along each offset - self.offset_lengths = np.sqrt( - np.sum((sampling*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() @@ -365,14 +368,16 @@ cdef class MCP: cdef FLOAT_T _travel_cost(self, FLOAT_T old_cost, FLOAT_T new_cost, FLOAT_T offset_length): - """ The travel cost for going from the current node to the next. + """ 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 cpdef int goal_reached(self, INDEX_T index, FLOAT_T cumcost): - """ int goal_reached(self, int index, float 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. @@ -388,23 +393,25 @@ cdef class MCP: return 0 - cdef void _examine_neighbor(self, INDEX_T index, INDEX_T new_index, FLOAT_T offset_length): - """ _examine_neighbor(self, int index, int new_index, float offset_length) + 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(self, int index, int new_index, float offset_length) + 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): + max_coverage=1.0, max_cumulative_cost=None, max_cost=None): """ Find the minimum-cost path from the given starting points. @@ -457,7 +464,8 @@ cdef class MCP: 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') + 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: @@ -602,7 +610,7 @@ cdef class MCP: # Calculate new cumulative cost new_cumcost = cumcost + self._travel_cost(cost, new_cost, - offset_length) + 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 @@ -620,8 +628,8 @@ cdef class MCP: # Un-flatten the costs and traceback arrays for human consumption. cumulative_costs = np.asarray(flat_cumulative_costs) - cumulative_costs = cumulative_costs.reshape(self.costs_shape, - order='F') + 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 @@ -720,7 +728,7 @@ cdef class MCP_Geometric(MCP): """ def __init__(self, costs, offsets=None, fully_connected=True, - sampling=None): + sampling=None): """__init__(costs, offsets=None, fully_connected=True, sampling=None) See class documentation. @@ -748,14 +756,13 @@ cdef class MCP_Connect(MCP): 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): + sampling=None): MCP.__init__(self, costs, offsets, fully_connected, sampling) # Create id map to keep track of origin of nodes @@ -784,7 +791,7 @@ cdef class MCP_Connect(MCP): cdef void _examine_neighbor(self, INDEX_T index, INDEX_T new_index, - FLOAT_T offset_length): + FLOAT_T offset_length): """ Check whether two fronts are meeting. If so, the flat_traceback is obtained and a connection is created. """ @@ -802,8 +809,8 @@ cdef class MCP_Connect(MCP): 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) + 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] @@ -833,16 +840,15 @@ cdef class MCP_Connect(MCP): pos2 : tuple The index of of the second neighbour in the connection. cost1 : float - The cumulative cost at pos1. + The cumulative cost at `pos1`. cost2 : float - The cumulative costs at pos2. - + The cumulative costs at `pos2`. """ pass - cdef void _update_node(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): """ Keep track of the id map so that we know which seed point a certain front originates from. """ @@ -855,19 +861,19 @@ cdef class MCP_Connect(MCP): cdef class MCP_Flexible(MCP): """MCP_Flexible(costs, offsets=None, fully_connected=True) - Find minimum cost paths through an n-d costs array. + 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 + 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(self, old_cost, new_cost, offset_length) + 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 @@ -877,10 +883,10 @@ cdef class MCP_Flexible(MCP): def examine_neighbor(self, INDEX_T index, INDEX_T new_index, - FLOAT_T offset_length): - """ examine_neighbor(self, index, new_index, offset_length) + 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. + 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 @@ -890,9 +896,9 @@ cdef class MCP_Flexible(MCP): pass - def update_node(self, INDEX_T index, INDEX_T new_index, - FLOAT_T offset_length): - """ update_node(self, index, new_index, offset_length) + 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. @@ -904,16 +910,16 @@ cdef class MCP_Flexible(MCP): pass - cdef FLOAT_T _travel_cost(self, FLOAT_T old_cost, - FLOAT_T new_cost, FLOAT_T offset_length): + 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): + 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): + 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/tests/test_anisotropy.py b/skimage/graph/tests/test_anisotropy.py index a2d3a240..199d2e73 100644 --- a/skimage/graph/tests/test_anisotropy.py +++ b/skimage/graph/tests/test_anisotropy.py @@ -8,7 +8,7 @@ import numpy as np a = np.ones((8, 8), dtype=np.float32) -horizontal_ramp = np.array( [[ 0., 1., 2., 3., 4., 5., 6., 7.,], +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.,], @@ -30,23 +30,23 @@ vertical_ramp = np.array( [[ 0., 0., 0., 0., 0., 0., 0., 0.,], 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) ] + 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) + # 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__": diff --git a/skimage/graph/tests/test_flexible.py b/skimage/graph/tests/test_flexible.py index 3570a31f..9ae50759 100644 --- a/skimage/graph/tests/test_flexible.py +++ b/skimage/graph/tests/test_flexible.py @@ -39,14 +39,14 @@ def test_flexible(): # Create MCP and do a traceback mcp = TestFlexibleMCP(a) - costs, traceback = mcp.find_costs([(0,0)]) + 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]]) + 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,