docs: def convex_hull_transform(image, levels=256, mask = None,
                            chunksize = CONVEX_HULL_CHUNKSIZE,
                            pass_cutoff = 16):
      '''Perform the convex hull transform of this image
  
      image - image composed of integer intensity values
      levels - # of levels that we separate the image into
      mask - mask of points to consider or None to consider all points
      chunksize - # of points processed in first pass of convex hull
  
      for each intensity value, find the convex hull of pixels at or above
      that value and color all pixels within the hull with that value.
      '''
      # Scale the image into the requisite number of levels
      if mask is None:
          img_min = np.min(image)
          img_max = np.max(image)
      else:
          unmasked_pixels = image[mask]
          if len(unmasked_pixels) == 0:
              return np.zeros(image.shape, image.dtype)
          img_min = np.min(unmasked_pixels)
          img_max = np.max(unmasked_pixels)
      img_shape = tuple(image.shape)
      if img_min == img_max:
          return image
  
      scale = (img_min +
               np.arange(levels).astype(image.dtype) *
               (img_max - img_min) / float(levels-1))
      image = ((image - img_min) * (levels-1) / (img_max - img_min))
      if mask is not None:
          image[~mask] = 0
      #
      # If there are more than 16 levels, we do the method first at a coarse
      # scale. The dark objects can produce points at every level, so doing
      # two passes can reduce the number of points in the second pass to
      # only the difference between two levels at the coarse pass.
      #
      if levels > pass_cutoff:
          sub_levels = int(np.sqrt(levels))
          rough_image = convex_hull_transform(np.floor(image), sub_levels)
          image = np.maximum(image, rough_image)
          del rough_image
      image = image.astype(int)
      #
      # Get rid of any levels that have no representatives
      #
      unique = np.unique(image)
      new_values = np.zeros(levels, int)
      new_values[unique] = np.arange(len(unique))
      scale = scale[unique]
      image = new_values[image]
      #
      # Start by constructing the list of points which are local maxima
      #
      min_image = grey_erosion(image, footprint = np.ones((3,3), bool)).astype(int)
      #
      # Set the borders of the min_image to zero so that the border pixels
      # will be in all convex hulls below their intensity
      #
      min_image[0, :] = 0
      min_image[-1, :] = 0
      min_image[:, 0] = 0
      min_image[:, -1] = 0
  
      i,j = np.mgrid[0:image.shape[0], 0:image.shape[1]]
      mask = image > min_image
  
      i = i[mask]
      j = j[mask]
      min_image = min_image[mask]
      image = image[mask]
      #
      # Each point that is a maximum is a potential vertex in the convex hull
      # for each value above the minimum. Therefore, it appears
      #
      # image - min_image
      #
      # times in the i,j,v list of points. So we can do a sum to calculate
      # the number of points, then use cumsum to figure out the first index
      # in that array of points for each i,j,v. We can then use cumsum
      # again on the array of points to assign their levels.
      #
      count = image - min_image
      npoints = np.sum(count)
      # The index in the big array of the first point to place for each
      # point
      first_index_in_big = np.cumsum(count) - count
      #
      # The big array can be quite big, for example if there are lots of
      # thin, dark objects. We do two passes of convex hull: the convex hull
      # of the convex hulls of several regions is the convex hull of the whole
      # so it doesn't matter too much how we break up the array.
      #
      first_i = np.zeros(0, int)
      first_j = np.zeros(0, int)
      first_levels = np.zeros(0, int)
      chunkstart = 0
      while chunkstart < len(count):
          idx = first_index_in_big[chunkstart]
          iend = idx + chunksize
          if iend >= npoints:
              chunkend = len(count)
              iend = npoints
          else:
              chunkend = np.searchsorted(first_index_in_big, iend)
              if chunkend < len(count):
                  iend = first_index_in_big[chunkend]
              else:
                  iend = npoints
          chunk_first_index_in_big = first_index_in_big[chunkstart:chunkend] - idx
          chunkpoints = iend - idx
          #
          # For the big array, construct an array of indexes into the small array
          #
          index_in_small = np.zeros(chunkpoints, int)
          index_in_small[0] = chunkstart
          index_in_small[chunk_first_index_in_big[1:]] = 1
          index_in_small = np.cumsum(index_in_small)
          #
          # We're going to do a cumsum to make the big array of levels. Point
          # n+1 broadcasts its first value into first_index_in_big[n+1].
          # The value that precedes it is image[n]. Therefore, in order to
          # get the correct value in cumsum:
          #
          # ? + image[n] = min_image[n+1]+1
          # ? = min_image[n+1] + 1 - image[n]
          #
          levels = np.ones(chunkpoints, int)
          levels[0] = min_image[chunkstart] + 1
          levels[chunk_first_index_in_big[1:]] = \
              min_image[chunkstart+1:chunkend] - image[chunkstart:chunkend-1] + 1
          levels = np.cumsum(levels)
          #
          # Construct the ijv
          #
          ijv = np.column_stack((i[index_in_small], j[index_in_small], levels))
          #
          # Get all of the convex hulls
          #
          pts, counts = convex_hull_ijv(ijv, np.arange(1, len(unique)))
          first_i = np.hstack((first_i, pts[:, 1]))
          first_j = np.hstack((first_j, pts[:, 2]))
          first_levels = np.hstack((first_levels, pts[:, 0]))
          chunkstart = chunkend
      #
      # Now do the convex hull of the reduced list of points
      #
      ijv = np.column_stack((first_i, first_j, first_levels))
      pts, counts = convex_hull_ijv(ijv, np.arange(1, len(unique)))
      #
      # Get the points along the lines described by the convex hulls
      #
      # There are N points for each label. Draw a line from each to
      # the next, except for the last which we draw from last to first
      #
      labels = pts[:, 0]
      i = pts[:, 1]
      j = pts[:, 2]
      first_index = np.cumsum(counts) - counts
      last_index = first_index + counts - 1
      next = np.arange(len(labels)) + 1
      next[last_index] = first_index
      index, count, i, j = get_line_pts(i, j, i[next], j[next])
      #
      # use a cumsum to get the index of each point from get_line_pts
      # relative to the labels vector
      #
      big_index = np.zeros(len(i), int)
      big_index[index[1:]] = 1
      big_index = np.cumsum(big_index)
      labels = labels[big_index]
      #
      # A given i,j might be represented more than once. Take the maximum
      # label at each i,j. First sort by i,j and label. Then take only values
      # that have a different i,j than the succeeding value. The last element
      # is always a winner.
      #
      order = np.lexsort((labels, i, j))
      i = i[order]
      j = j[order]
      labels = labels[order]
      mask = np.hstack((((i[1:] != i[:-1]) | (j[1:] != j[:-1])),[True]))
      i = i[mask]
      j = j[mask]
      labels = labels[mask]
      #
      # Now, we have an interesting object. It's ordered by j, then i which
      # means that we have scans of interesting i at each j. The points
      # that aren't represented should have the minimum of the values
      # above and below.
      #
      # We can play a cumsum trick to do this, placing the difference
      # of a point with its previous in the 2-d image, then summing along
      # each i axis to set empty values to the value of the nearest occupied
      # value above and a similar trick to set empty values to the nearest
      # value below. We then take the minimum of the two results.
      #
      first = np.hstack(([True], j[1:] != j[:-1]))
      top = np.zeros(img_shape, labels.dtype)
      top[i[first], j[first]] = labels[first]
      top[i[~first], j[~first]] = (labels[1:] - labels[:-1])[~first[1:]]
      top = np.cumsum(top, 0)
  
      # From 0 to the location of the first point, set to value of first point
      bottom = np.zeros(img_shape, labels.dtype)
      bottom[0, j[first]] = labels[first]
      # From 1 + the location of the previous point, set to the next point
      last = np.hstack((first[1:], [True]))
      bottom[i[:-1][~first[1:]]+1, j[~first]] = (labels[1:] - labels[:-1])[~first[1:]]
      # Set 1 + the location of the last point to -labels so that all past
      # the end will be zero (check for i at end...)
      llast = last & (i < img_shape[0] - 1)
      bottom[i[llast]+1, j[llast]] = - labels[llast]
      bottom = np.cumsum(bottom, 0)
      image = np.minimum(top, bottom)
      return scale[image]