diff --git a/scikits/image/graph/_mcp.pyx b/scikits/image/graph/_mcp.pyx index c1a6429e..e06f0d91 100644 --- a/scikits/image/graph/_mcp.pyx +++ b/scikits/image/graph/_mcp.pyx @@ -210,6 +210,9 @@ cdef class MCP: `costs` array at each point on the path. The class MCP_Geometric, on the other hand, accounts for the fact that diagonal vs. axial moves are of different lengths, and weights the path cost accordingly. + + Array elements with infinite or negative costs will simply be ignored, as + will paths whose cumulative cost overflows to infinite. Parameters ---------- @@ -241,8 +244,6 @@ cdef class MCP: See class documentation. """ costs = np.asarray(costs) - if costs.min() < 0: - raise ValueError('minimum cost must be positive') if not np.can_cast(costs.dtype, FLOAT): raise TypeError('cannot cast costs array to ' + str(FLOAT)) @@ -476,19 +477,26 @@ 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 + new_cost = flat_costs[new_index] + if new_cost < 0 or new_cost == inf: + continue + # Now we ask the heap to append or update the cost to # this new point, but only if that point isn't already # in the heap, or it is but the new cost is lower. travel_cost = self._travel_cost(flat_costs[index], - flat_costs[new_index], + new_cost, offset_lengths[i]) - costs_heap.push_if_lower_fast(cost + travel_cost, new_index) - # If we did perform an append or update, we should - # record the offset from the predecessor to this new - # point - if costs_heap._pushed: - traceback_offsets[new_index] = i + # 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 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, diff --git a/scikits/image/graph/heap.pxd b/scikits/image/graph/heap.pxd index b86a0027..5f3b17ee 100644 --- a/scikits/image/graph/heap.pxd +++ b/scikits/image/graph/heap.pxd @@ -22,7 +22,7 @@ cdef class BinaryHeap: cdef void _remove(self, int i) cdef int push_fast(self, double value, int reference) - cdef float pop_fast(self) + cdef VALUE_T pop_fast(self) cdef class FastUpdateBinaryHeap(BinaryHeap): cdef readonly int max_reference @@ -30,6 +30,6 @@ cdef class FastUpdateBinaryHeap(BinaryHeap): cdef int _invalid_ref cdef int _pushed - cdef float value_of_fast(self, int reference) + cdef VALUE_T value_of_fast(self, int reference) cdef int push_if_lower_fast(self, double value, int reference) \ No newline at end of file diff --git a/scikits/image/graph/heap.pyx b/scikits/image/graph/heap.pyx index 2345f2f9..59822a26 100644 --- a/scikits/image/graph/heap.pyx +++ b/scikits/image/graph/heap.pyx @@ -339,7 +339,7 @@ cdef class BinaryHeap: return count - cdef float pop_fast(self): + cdef VALUE_T pop_fast(self): """The c-method for fast popping. Returns the minimum value. The reference is put in self._popped_ref""" @@ -366,7 +366,7 @@ cdef class BinaryHeap: # get values cdef int ir = i - ((1 << levels) - 1) #(2**self.levels-1) # LevelStart - cdef float value = values[i] + cdef VALUE_T value = values[i] self._popped_ref = self._references[ir] # remove it @@ -645,7 +645,7 @@ cdef class FastUpdateBinaryHeap(BinaryHeap): return ir - cdef float value_of_fast(self, int reference): + cdef VALUE_T value_of_fast(self, int reference): """Return the value corresponding to the given reference. If inf is returned, the reference may be invalid: check the _invaild_ref field in this case.""" diff --git a/scikits/image/graph/tests/test_mcp.py b/scikits/image/graph/tests/test_mcp.py index f9f62e8e..cb285618 100644 --- a/scikits/image/graph/tests/test_mcp.py +++ b/scikits/image/graph/tests/test_mcp.py @@ -43,6 +43,32 @@ def test_basic(): (6, 1), (7, 2)]) +def test_neg_inf(): + expected_costs = numpy.where(a==1, np.inf, 0) + expected_path = [(1, 6), + (1, 5), + (1, 4), + (1, 3), + (1, 2), + (2, 1), + (3, 1), + (4, 1), + (5, 1), + (6, 1)] + test_neg = numpy.where(a==1, -1, 0) + test_inf = numpy.where(a==1, np.inf, 0) + m = mcp.MCP(test_neg, fully_connected=True) + costs, traceback = m.find_costs([(1, 6)]) + return_path = m.traceback((6, 1)) + assert_array_equal(costs, expected_costs) + assert_array_equal(return_path, expected_path) + m = mcp.MCP(test_inf, fully_connected=True) + costs, traceback = m.find_costs([(1, 6)]) + return_path = m.traceback((6, 1)) + assert_array_equal(costs, expected_costs) + assert_array_equal(return_path, expected_path) + + def test_route(): return_path, cost = mcp.route_through_array(a, (1,6), (7,2), geometric=True) assert_almost_equal(cost, np.sqrt(2)/2)