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]