docs: def skeletonize(image, mask=None, ordering = None):
      '''Skeletonize the image
      
      Take the distance transform.
      Order the 1 points by the distance transform.
      Remove a point if it has more than 1 neighbor and if removing it
      does not change the Euler number.
      
      image - the binary image to be skeletonized
      
      mask - only skeletonize pixels within the mask
      
      ordering - a matrix of the same dimensions as the image. The matrix
                 provides the ordering of the erosion with the lowest values
                 being eroded first. The default is to use the distance transform.
      '''
      global eight_connect
      if mask is None:
          masked_image = image
      else:
          masked_image = image.astype(bool).copy()
          masked_image[~mask] = False
      #
      # Lookup table - start with only positive pixels.
      # Keep if # pixels in neighborhood is 2 or less
      # Keep if removing the pixel results in a different connectivity
      #
      table = (make_table(True,np.array([[0,0,0],[0,1,0],[0,0,0]],bool),
                          np.array([[0,0,0],[0,1,0],[0,0,0]],bool)) &
               (np.array([scind.label(pattern_of(index), eight_connect)[1] !=
                          scind.label(pattern_of(index & ~ 2**4),
                                      eight_connect)[1]
                          for index in range(512) ]) |
                np.array([np.sum(pattern_of(index))<3 for index in range(512)])))
      
      if ordering is None:
          distance = scind.distance_transform_edt(masked_image)
      else:
          distance = ordering
      #
      # The processing order along the edge is critical to the shape of the
      # resulting skeleton: if you process a corner first, that corner will
      # be eroded and the skeleton will miss the arm from that corner. Pixels
      # with fewer neighbors are more "cornery" and should be processed last.
      #
      cornerness_table = np.array([9-np.sum(pattern_of(index))
                                   for index in range(512)])
      corner_score = table_lookup(masked_image, cornerness_table, False,1)
      i,j = np.mgrid[0:image.shape[0],0:image.shape[1]]
      result=masked_image.copy()
      distance = distance[result]
      i = np.ascontiguousarray(i[result],np.int32)
      j = np.ascontiguousarray(j[result],np.int32)
      result=np.ascontiguousarray(result,np.uint8)
      #
      # We use a random # for tiebreaking. Assign each pixel in the image a
      # predictable, random # so that masking doesn't affect arbitrary choices
      # of skeletons
      #
      np.random.seed(0)
      tiebreaker=np.random.permutation(np.arange(np.product(masked_image.shape)))
      tiebreaker.shape=masked_image.shape
      order = np.lexsort((tiebreaker[masked_image],
                          corner_score[masked_image],
                          distance))
      order = np.ascontiguousarray(order, np.int32)
      table = np.ascontiguousarray(table, np.uint8)
      skeletonize_loop(result, i, j, order, table)
      
      result = result.astype(bool)
      if not mask is None:
          result[~mask] = image[~mask]
      return result