docs: def find_neighbors(labels):
      '''Find the set of objects that touch each object in a labels matrix
      
      Construct a "list", per-object, of the objects 8-connected adjacent
      to that object.
      Returns three 1-d arrays:
      * array of #'s of neighbors per object
      * array of indexes per object to that object's list of neighbors
      * array holding the neighbors.
      
      For instance, say 1 touches 2 and 3 and nobody touches 4. The arrays are:
      [ 2, 1, 1, 0], [ 0, 2, 3, 4], [ 2, 3, 1, 1]
      '''
      max_label = np.max(labels)
      # Make a labels matrix with zeros around the edges so we can do index
      # offsets without worrying.
      #
      new_labels = np.zeros(np.array(labels.shape)+2, labels.dtype)
      new_labels[1:-1,1:-1] = labels
      labels = new_labels
      # Only consider the points that are next to others
      adjacent_mask = adjacent(labels)
      adjacent_i, adjacent_j = np.argwhere(adjacent_mask).transpose()
      # Get matching vectors of labels and neighbor labels for the 8
      # compass directions.
      count = len(adjacent_i)
      if count == 0:
          return (np.zeros(max_label, int),
                  np.zeros(max_label,int),
                  np.zeros(0,int))
      # The following bizarre construct does the following:
      # labels[adjacent_i, adjacent_j] looks up the label for each pixel
      # [...]*8 creates a list of 8 references to it
      # np.hstack concatenates, giving 8 repeats of the list
      v_label = np.hstack([labels[adjacent_i, adjacent_j]]*8)
      v_neighbor = np.zeros(count*8,int)
      index = 0
      for i,j in ((-1,-1),(-1,0),(-1,1),(0,-1),(0,1),(1,-1),(1,0),(1,1)):
          v_neighbor[index:index+count] = labels[adjacent_i+i, adjacent_j+j]
          index += count
      #
      # sort by label and neighbor
      #
      sort_order = np.lexsort((v_neighbor,v_label))
      v_label = v_label[sort_order]
      v_neighbor = v_neighbor[sort_order]
      #
      # eliminate duplicates by comparing each element after the first one
      # to its previous
      #
      first_occurrence = np.ones(len(v_label), bool)
      first_occurrence[1:] = ((v_label[1:] != v_label[:-1]) |
                              (v_neighbor[1:] != v_neighbor[:-1]))
      v_label = v_label[first_occurrence]
      v_neighbor = v_neighbor[first_occurrence]
      #
      # eliminate neighbor = self and neighbor = background
      #
      to_remove = (v_label == v_neighbor) | (v_neighbor == 0)
      v_label = v_label[~ to_remove]
      v_neighbor = v_neighbor[~ to_remove]
      #
      # The count of # of neighbors
      #
      v_count = fixup_scipy_ndimage_result(scind.sum(np.ones(v_label.shape),
                                                     v_label,
                                                     np.arange(max_label,dtype=np.int32)+1))
      v_count = v_count.astype(int)
      #
      # The index into v_neighbor
      #
      v_index = np.cumsum(v_count)
      v_index[1:] = v_index[:-1]
      v_index[0] = 0
      return (v_count, v_index, v_neighbor)