from pylab import *
import pygis.path

def vlen(vec):
    """Return the Euclidean length of the vector vec."""
    return sqrt(vdot(vec, vec))

def line_cost(start, end, grad):
    """Given some gradient raster, calculate the cost of the line start -> end.
    
    """

    start = array(start)
    end = array(end)
    delta = end - start
    delta_len = vlen(delta)
    direction = delta / delta_len

    step = array(grad.pixel_linear_shape).min()
    line_points, dists = pygis.path.points_on_line(start, end, step)
    return integrate_cost(line_points, grad)

def integrate_cost(points, grad):
    """Given some gradient raster, calculate the cost of the lines segments
    joining the points assuming that the size of each line segment is small
    w.r.t pixel size.
    
    """

    points = np.array(points)
    assert len(points.shape) > 1 and points.shape[0] > 1

    midpoints = 0.5*(points[:-1] + points[1:])
    deltas = points[1:] - points[:-1]
    delta_lengths = np.sqrt((deltas * deltas).sum(axis=1)).transpose()
    directions = np.array(deltas, copy=True)
    directions[:,0] /= delta_lengths
    directions[:,1] /= delta_lengths

    samples = grad.lanczos_sample(midpoints)
    if len(points) == 2:
        samples = [samples,]

    cost = 0.0
    for direction, dh, gradient in zip(directions, delta_lengths, samples):
        # get vert. distance in metres
        dv = vdot(gradient, direction) * dh

        # get linear distance
        d = vlen(array(dh, dv))

        # calculate the subjective gradient
        subj_gradient = vdot(gradient, direction)

        # increase cost
        cost += abs(subj_gradient) * d

    return cost

def path_cost(points, grad):
    """Given a list of points specifying a path and a gradient raster, return the total cost of the path.

    """

    step = array(grad.pixel_linear_shape).min()
    int_points = []
    for start, end in zip(points[:-1], points[1:]):
        np, _ = pygis.path.points_on_line(start, end, step)
        int_points.extend(np[:-1])
    int_points.append(points[-1])
    return integrate_cost(int_points, grad)

# vim:filetype=python:sts=4:et:ts=4:sw=4

@@ -24,546 +24,37 @@
 ort 
-pygis.path%0A%0Adef vlen(vec):%0A    %22%22%22Return the Euclidean length of the vector vec.%22%22%22%0A    return sqrt(vdot(vec, vec))%0A%0Adef line_cost(start, end, grad):%0A    %22%22%22Given some gradient raster, calculate the cost of the line start -%3E end.%0A    %0A    %22%22%22%0A%0A    start = array(start)%0A    end = array(end)%0A    delta = end - start%0A    delta_len = vlen(delta)%0A    direction = delta / delta_len%0A%0A    step = array(grad.pixel_linear_shape).min()%0A    line_points, dists = pygis.path.points_on_line(start, end, step)%0A    return integrate_cost(line_points, grad)
+numpy as np%0Aimport pygis.path
 %0A%0Ade
@@ -413,20 +413,24 @@
 %5D)%0A    d
-elta
+irection
 s = poin
@@ -454,21 +454,23 @@
 -1%5D%0A    
-delta
+segment
 _lengths
@@ -486,21 +486,29 @@
 t((d
-eltas * delta
+irections * direction
 s).s
@@ -535,53 +535,8 @@
 e()%0A
-    directions = np.array(deltas, copy=True)%0A
     
@@ -550,29 +550,31 @@
 ons%5B:,0%5D /= 
-delta
+segment
 _lengths%0A   
@@ -593,21 +593,23 @@
 :,1%5D /= 
-delta
+segment
 _lengths
@@ -713,404 +713,563 @@
 s,%5D%0A
-%0A
     
-cost = 0.0%0A    for direction, dh, gradient in zip(directions, delta_lengths, samples):%0A        # get vert. distance in metres%0A        dv = vdot(gradient, direction) * dh%0A%0A        # get linear distance%0A        d = vlen(array(dh, dv))%0A%0A   
+samples = np.array(samples)%0A    %0A    # samples: Nx2 array of gradient samples%0A    # directions: Nx2 array of normalised segment directions%0A    # segment_lengths: Nx2 array of segment horizontal lengths%0A%0A    # calculate subjective slope%0A    subj_slope = (samples * directions).sum(axis=1)%0A%0A    # hence vertical displacement%0A    vdisp = subj_slope * segment_lengths%0A%0A
     
- 
 # 
-calculate the subjective gradient%0A        subj_gradient = vdot(gradient, direction)%0A%0A        # increase cost%0A   
+hence Euclidean displacement%0A    disp = np.sqrt(vdisp*vdisp + segment_lengths*segment_lengths)%0A%0A    # cost is distance * absolute slope%0A
     
- 
 cost 
-+
 = 
-abs(subj_gradient) * d
+(disp * np.abs(subj_slope)).sum()
 %0A%0A  
@@ -1439,16 +1439,23 @@
  step = 
+np.abs(
 array(gr
@@ -1467,21 +1467,20 @@
 xel_
-linear
+proj
 _shape)
+)
 .min
@@ -1562,18 +1562,19 @@
         
-n
 p
+ts
 , _ = py
@@ -1641,18 +1641,19 @@
 .extend(
-n
 p
+ts
 %5B:-1%5D)%0A 
@@ -1730,16 +1730,463 @@
  grad)%0A%0A
+def line_cost(start, end, grad):%0A    %22%22%22Given some gradient raster, calculate the cost of the line start -%3E end.%0A    %0A    %22%22%22%0A%0A    start = array(start)%0A    end = array(end)%0A    delta = end - start%0A    delta_len = np.sqrt(np.dot(delta, delta))%0A    direction = delta / delta_len%0A%0A    step = np.abs(array(grad.pixel_proj_shape)).min()%0A    line_points, dists = pygis.path.points_on_line(start, end, step)%0A    return integrate_cost(line_points, grad)%0A%0A
 # vim:fi

tidy up path cost calculation routines