~~~
def calculate_stop_paths(init, fin):
    # Goal: map each start point to its nearest end point and each end point to its
    # nearest start point, then return the mappings as the target path.

    available_starts = kd_tree.KDTree([(x[0],) for x in init], 1)
    available_ends = kd_tree.KDTree([(x[0],) for x in fin], 1)

    init_map = dict((x[0], x) for x in init)
    fin_map = dict((x[0], x) for x in fin)

    forward_map = defaultdict(list)
    cover_count = defaultdict(lambda: 0)
    
    # Map each start point to its nearest ending point
    for stop in init:
        ratio = stop[0]
        match = available_ends.get_nearest((ratio,), False)[0]

        # add a path ratio -> match
        forward_map[ratio].append(match)
        cover_count[match] += 1
        
    # Map each unused end point to its nearest starting point
    for stop in fin:
        ratio = stop[0]
        if cover_count[ratio] > 0:
            continue
        match = available_starts.get_nearest((ratio,), False)[0]

        # If this point is covering another redundantly, prefer to remap it
        # rather than double-mapping it
        if forward_map[match]:
            potential_redundancy = forward_map[match][0]
            if cover_count[potential_redundancy] > 1:
                forward_map[match].remove(potential_redundancy)
                cover_count[potential_redundancy] -= 1
        
        forward_map[match].append(ratio)
    
    for start in forward_map:
        for end in forward_map[start]:
            yield init_map[start], fin_map[end]

def interpolated_stop(start, end, t):
    ratio = t * end[0] + (1-t) * start[0]
    color, alpha = interpolated_color(start[1], start[2], end[1], end[2], t)
    return (ratio, color, alpha)

def interpolated_color(colx, ax, coly, ay, t):
    rx, gx, bx = split_colors(colx)
    ry, gy, by = split_colors(coly)
    ai = interpolate_value(ax, ay, t)

    if ai == 0:
        return '#FFFFFF', 0

    ri = round(interpolate_value(rx*ax, ry*ay, t)/ai)
    gi = round(interpolate_value(gx*ax, gy*ay, t)/ai)
    bi = round(interpolate_value(bx*ax, by*ay, t)/ai)

    return "#%02X%02X%02X" % (ri, gi, bi), ai
~~~