def label_skeleton(skeleton):
    '''Label a skeleton so that each edge has a unique label
    
    This operation produces a labels matrix where each edge between
    two branchpoints has a different label. If the skeleton has been
    properly eroded, there are three kinds of points:
    1) point adjacent to 0 or 1 other points = end of edge
    2) point adjacent to two other points = in middle of edge
    3) point adjacent to more than two other points = at end of edge
            connecting to another edge
    4) a branchpoint
    
    We do all connected components here where components are 8-connected
    but a point in category 3 can't connect to another point in category 3.
    
    Returns the labels matrix and the count as a tuple
    '''
    bpts = branchpoints(skeleton)
    #
    # Count the # of neighbors per point
    #
    neighbors = scind.convolve(skeleton.astype(int), np.ones((3,3),int),
                               mode='constant').astype(int)
    neighbors[~skeleton] = 0
    neighbors[skeleton] -= 1
    #
    # Find the i/j coordinates of the relevant points
    #
    i,j = np.mgrid[0:skeleton.shape[0], 0:skeleton.shape[1]]
    skeleton_minus_bpts = skeleton & ~ bpts
    si = i[skeleton_minus_bpts]
    sj = j[skeleton_minus_bpts]
    bi = i[bpts]
    bj = j[bpts]
    i = np.hstack((bi, si))
    j = np.hstack((bj, sj))
    b_vnum = np.arange(len(bi))
    s_vnum = np.arange(len(si)) + len(bi)
    all_vnum = np.hstack((b_vnum, s_vnum))
    vertex_numbers=np.zeros(skeleton.shape, int)
    vertex_numbers[i,j] = all_vnum
    #
    # src and dest are the vertices linked by edges. Their values are the
    # vertex numbers. First, link every vertex to itself
    #
    src = all_vnum
    dest = all_vnum
    #
    # Now, for the non-branchpoints, link to all 8-connected neighbors
    # while obeying the rules
    #
    for ioff, joff in ((-1,-1), (-1,0), (-1,1),
                       ( 0,-1),         ( 0,1),
                       ( 1,-1), ( 1,0), ( 1,1)):
        consider = np.ones(len(si), bool)
        if ioff == -1:
            consider = si > 0
        elif ioff == 1:
            consider = si < skeleton.shape[0] - 1
        if joff == -1:
            consider = consider & (sj > 0)
        elif joff == 1:
            consider = consider & (sj < skeleton.shape[1] - 1)
        #
        # Forge a link if the offset point is in the skeleton
        #
        ci = si[consider]
        cj = sj[consider]
        link = (skeleton_minus_bpts[ci+ioff, cj+joff] &
                ((neighbors[ci,cj] < 3) | (neighbors[ci+ioff, cj+joff] < 3)))
        ci = ci[link]
        cj = cj[link]
        src = np.hstack((src, vertex_numbers[ci, cj]))
        dest = np.hstack((dest, vertex_numbers[ci+ioff, cj+joff]))
        
    labeling = all_connected_components(src, dest)
    vertex_numbers[i,j] = labeling + 1
    return (vertex_numbers, 
            0 if len(labeling) == 0 else int(np.max(labeling)) + 1)