Merge branch 'mcp_infinite' of ssh://github.com/zachrahan/scikits.image

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,

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)
     

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."""

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)