docs: 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)