Getting distance and cost from an MCP

[David Orme]

Hi,

I have a solution, but it is a hack and only works for the special case 
where the distance (not cost) within cells is constant. I suspect it is 
also a hostage to floating point precision.

Basically, add a downscaled resolution to the cost values - downscaled so 
that this additional value doesn't interfere with the cost pathways. Then 
run find_costs() on both versions of the cost matrix and the difference 
(scaled back up) is the physical distance. I think the traceback matrix 
provides a test for when the hack fails by interfering with the pathing.

If anyone has anything more elegant, I'd love to know about it!

Cheers,
David

from skimage import graph

import numpy


cost = numpy.array([[9,1,9,9,9,9,9,9,9],

[9,1,9,9,9,9,9,9,9],

[9,2,9,9,9,9,9,9,9],

[9,1,9,9,9,9,9,9,9],

[9,1,9,9,9,9,9,9,9],

[9,1,9,9,9,9,9,9,9],

[9,3,9,9,9,9,9,1,9],

[9,1,9,9,9,9,9,1,9],

[9,1,9,9,9,9,9,2,9],

[9,1,9,9,9,9,9,1,9],

[9,2,9,9,9,9,9,1,9],

[9,1,9,9,9,9,9,3,9],

[9,1,9,9,9,9,9,1,9],

[9,1,2,1,1,3,1,1,9],

[9,9,9,9,9,9,9,9,9]], dtype='float')


# add the resolution to each cost, downscaled so that

# it won't interfere with path choice

resolution = 250

scale = 1e7

costPlusDist = cost + resolution/scale


# get two sets of cumulative costs

costMCP = graph.MCP_Geometric(cost, fully_connected=True)

cumcost, trace = costMCP.find_costs([(0,1)])

costDistMCP = graph.MCP_Geometric(costPlusDist, fully_connected=True)

cumcostdist, trace2 = costDistMCP.find_costs([(0,1)])


if (trace == trace2).all():

    dist = (cumcostdist - cumcost)*scale

else:

    print "Hack failed - different minimum cost paths used"





On Tuesday, 22 April 2014 14:59:34 UTC+1, David Orme wrote:
>
> Hi,
>
> I'm running some road building simulations in Python across a cost 
> surface. What I'd like to be able to do is to get both the cost of paths 
>  to a set of starting points and the actual physical distance that those 
> paths would cover.
>
> I can easily find the cumulative costs of paths using MCP_Geometric from 
> the graph module but I can't figure out how to get a surface of cumulative 
> distance along the least cost paths. I could iterate over each cell in the 
> cost array in turn to get the cheapest path and fetch the distance using 
> pythogoras on the steps, but my intuition from the documentation is that 
> the MCP_Flexible class in 0.10dev  - particularly the travel_cost method - 
> might allow me to get at this more quickly.
>
> In the toy example, below I'd like to find costs from a start point at 
> [0,1] in the cost array. Most of the array is expensive to traverse, but 
> the valley means that there are big differences in the length of the 
> cheapest road needed to get between points. If the entries in the cost 
> array represent cells with a resolution of 250 metres across, the path from 
> [0,1] to [0,5] has a cost of 32.0 along a road of length 1.5 km (4 steps 
> E). In contrast, the path from [0,1] to [6,7] has a cost of 35.04 along a 
> road of length 6.21 km (12 steps S - 3km, 1 step SE - 0.354km, 4 steps E - 
> 1km, 1 step NE - 0.354km, and 6 steps N - 1.5km).
>
> Any suggestions would be very welcome.
> Cheers,
> David
>
> import numpy
> from skimage import graph
>
> cost = numpy.array([[9,1,9,9,9,9,9,9,9],
>                     [9,1,9,9,9,9,9,9,9],
>                     [9,2,9,9,9,9,9,9,9],
>                     [9,1,9,9,9,9,9,9,9],
>                     [9,1,9,9,9,9,9,9,9],
>                     [9,1,9,9,9,9,9,9,9],
>                     [9,3,9,9,9,9,9,1,9],
>                     [9,1,9,9,9,9,9,1,9],
>                     [9,1,9,9,9,9,9,2,9],
>                     [9,1,9,9,9,9,9,1,9],
>                     [9,2,9,9,9,9,9,1,9],
>                     [9,1,9,9,9,9,9,3,9],
>                     [9,1,9,9,9,9,9,1,9],
>                     [9,1,2,1,1,3,1,1,9],
>                     [9,9,9,9,9,9,9,9,9]], dtype='float')
>
> costMCP = graph.MCP_Geometric(cost, fully_connected=True)
> cumcost, trace = costMCP.find_costs([(0,1)])
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/scikit-image/attachments/20140423/87fa658c/attachment.html>