def convex_hull_ijv(pixel_labels, indexes, fast=True):
    '''Return the convex hull for each label using an ijv labeling
    
    pixel_labels: the labeling of the pixels in i,j,v form where
                  i & j are the coordinates of a pixel and v is
                  the pixel's label number
    indexes: the indexes at which to measure the convex hull

    Returns a matrix and a vector. The matrix consists of one row per
    point in the convex hull. Each row has three columns, the label #,
    the i coordinate of the point and the j coordinate of the point. The
    result is organized first by label, then the points are arranged
    counter-clockwise around the perimeter.
    The vector is a vector of #s of points in the convex hull per label
    '''
    if fast:
        return _convex_hull.convex_hull_ijv(pixel_labels, indexes)

    # We keep this code for testing.
    if len(indexes) == 0:
        return np.zeros((0,3),int),np.zeros((0,),int)
    #
    # An array that converts from label # to index in "indexes"
    anti_indexes = np.zeros((np.max(indexes)+1,),int)
    anti_indexes[indexes] = list(range(len(indexes)))

    coords = pixel_labels[:,:2]
    i = coords[:, 0]
    j = coords[:, 1]
    # This disgusting copy spooge appears to be needed for scipy 0.7.0
    labels_per_point = np.zeros(len(pixel_labels), int)
    labels_per_point[:] = pixel_labels[:,2]
    #
    # Calculate the centers for each label
    #
    center_i = fixup_scipy_ndimage_result(
        scind.mean(i.astype(float), labels_per_point, indexes))
    center_j = fixup_scipy_ndimage_result(
        scind.mean(j.astype(float), labels_per_point, indexes))
    centers = np.column_stack((center_i, center_j))
    #
    # Now make an array with one outline point per row and the following
    # columns:
    #
    # index of label # in indexes array
    # angle of the point relative to the center
    # i coordinate of the point
    # j coordinate of the point
    #
    anti_indexes_per_point = anti_indexes[labels_per_point]
    centers_per_point = centers[anti_indexes_per_point]
    angle = np.arctan2(i-centers_per_point[:,0],j-centers_per_point[:,1])
    a = np.zeros((len(i),3), np.int32)
    a[:,0] = anti_indexes_per_point
    a[:,1:] = coords
    #
    # Sort the array first by label # (sort of), then by angle
    #
    order = np.lexsort((angle,anti_indexes_per_point))
    #
    # Make unique
    #
    same_as_next = np.hstack([np.all(a[order[:-1],:] == a[order[1:],:], 1), [False]])
    order = order[~same_as_next]
    a=a[order]
    anti_indexes_per_point = anti_indexes_per_point[order]
    angle = angle[order]
    centers_per_point = centers_per_point[order]
    #
    # Make the result matrix, leaving enough space so that all points might
    # be on the convex hull.
    #
    result = np.zeros((len(order),3), np.int32)
    result[:,0] = labels_per_point[order]
    #
    # Create an initial count vector
    #
    v = np.ones((a.shape[0],),dtype = np.int32)
    result_counts = scipy.sparse.coo_matrix((v,(a[:,0],v*0)),
                                            shape=(len(indexes),1))
    result_counts = result_counts.toarray().flatten()
    r_anti_indexes_per_point = anti_indexes_per_point # save this
    #
    # Create a vector that indexes into the results for each label
    #
    result_index = np.zeros(result_counts.shape, np.int32)
    result_index[1:]=np.cumsum(result_counts[:-1])
    #
    # Initialize the counts of convex hull points to a ridiculous number
    #
    counts = np.iinfo(np.int32).max
    first_pass = True
    while True:
        #
        # Figure out how many putative convex hull points there are for
        # each label.
        #
        # If the count for a label is 3 or less, it's a convex hull or
        # degenerate case.
        #
        # If the count hasn't changed in an iteration, then we've done
        # as well as we can hope to do.
        #
        v = np.ones((a.shape[0],),dtype = np.int32)
        new_counts = scipy.sparse.coo_matrix((v,(a[:,0],v*0)),
                                             shape=(len(indexes),1))
        new_counts = new_counts.toarray().flatten()
        done_count = (2 if first_pass else 3)
        finish_me = ((new_counts > 0) & 
                     ((new_counts <= done_count) | 
                      (new_counts == counts)))
        indexes_to_finish = np.argwhere(finish_me).astype(np.int32).squeeze()
        keep_me = (new_counts > done_count) & (new_counts < counts)
        indexes_to_keep = np.argwhere(keep_me).astype(np.int32).squeeze()
        if len(indexes_to_finish):
            result_counts[finish_me] = new_counts[finish_me]
            #
            # Store the coordinates of each of the points to finish
            #
            finish_this_row = finish_me[a[:,0]]
            rows_to_finish = np.argwhere(finish_this_row).flatten()
            a_to_finish = a[rows_to_finish]
            atf_indexes = a_to_finish[:,0]
            #
            # Map label #s to the index into indexes_to_finish of that label #
            #
            anti_indexes_to_finish = np.zeros((len(indexes),), np.int32)
            anti_indexes_to_finish[indexes_to_finish] = list(range(len(indexes_to_finish)))
            #
            # Figure out the indices of each point in a label to be finished.
            # We figure out how much to subtract for each label, then
            # subtract that much from 0:N to get successive indexes at
            # each label.
            # Then we add the result_index to figure out where to store it
            # in the result table.
            #
            finish_idx_base = np.zeros((len(indexes_to_finish),), np.int32)
            finish_idx_base[1:]=np.cumsum(new_counts[indexes_to_finish])[:-1]
            finish_idx_bases = finish_idx_base[anti_indexes_to_finish[atf_indexes]]
            finish_idx = (np.array(list(range(a_to_finish.shape[0])))-
                          finish_idx_bases)
            finish_idx = finish_idx + result_index[atf_indexes]
            result[finish_idx,1:] = a_to_finish[:,1:]
        if len(indexes_to_keep) == 0:
            break
        #
        # Figure out which points are still available
        #
        rows_to_keep = np.argwhere(keep_me[a[:,0].astype(np.int32)]).flatten()
        rows_to_keep = rows_to_keep.astype(np.int32)
        a = a[rows_to_keep]
        centers_per_point = centers_per_point[rows_to_keep]
        counts = new_counts
        #
        # The rule is that the area of the triangle from the center to
        # point N-1 to point N plus the area of the triangle from the center
        # to point N to point N+1 must be greater than the area of the
        # triangle from the center to point N-1 to point N+1 for a point
        # to be on the convex hull.
        # N-1 and N+1 have to be modulo "counts", so we make special arrays
        # to address those situations.
        #
        anti_indexes_to_keep = np.zeros((len(indexes),), np.int32)
        anti_indexes_to_keep[indexes_to_keep] = list(range(len(indexes_to_keep)))
        idx_base = np.zeros((len(indexes_to_keep),), np.int32)
        idx_base[1:]=np.cumsum(counts[keep_me])[0:-1]
        idx_bases = idx_base[anti_indexes_to_keep[a[:,0]]]
        counts_per_pt = counts[a[:,0]]
        idx = np.array(list(range(a.shape[0])), np.int32)-idx_bases
        n_minus_one = np.mod(idx+counts_per_pt-1,counts_per_pt)+idx_bases
        n_plus_one  = np.mod(idx+1,counts_per_pt)+idx_bases
        #
        # Compute the triangle areas
        #
        t_left_plus_right = (triangle_areas(centers_per_point,
                                            a[n_minus_one, 1:],
                                            a[:, 1:]) +
                             triangle_areas(centers_per_point,
                                            a[:, 1:],
                                            a[n_plus_one, 1:]))
        t_lr = triangle_areas(centers_per_point,
                              a[n_minus_one, 1:], a[n_plus_one, 1:])
        #
        # Keep the points where the area of the left triangle plus the
        # area of the right triangle is bigger than the area of the triangle
        # composed of the points to the left and right. This means that
        # there's a little triangle sitting on top of t_lr with our point
        # on top and convex in relation to its neighbors.
        #
        keep_me = t_left_plus_right > t_lr
        #
        # If all points on a line are co-linear with the center, then the
        # whole line goes away. Special handling for this to find the points
        # most distant from the center and on the same side
        #
        consider_me = t_left_plus_right == 0
        del t_lr, t_left_plus_right
        if np.any(consider_me):
            diff_i = a[:,1]-centers_per_point[:,0]
            diff_j = a[:,2]-centers_per_point[:,1]
            #
            # The manhattan distance is good enough
            #
            dist = np.abs(diff_i)+np.abs(diff_j)
            # The sign is different on different sides of a line including
            # the center. Multiply j by 2 to keep from colliding with i
            #
            # If both signs are zero, then the point is in the center
            #
            # Compute np.sign(diff_i, out=diff_i) + np.sign(diff_j) * 2,
            # but without any temporaries.
            sign = np.sign(diff_j, diff_j)
            sign *= 2
            sign += np.sign(diff_i)
            del diff_i, diff_j

            n_minus_one_consider = n_minus_one[consider_me]
            n_plus_one_consider = n_plus_one[consider_me]
            left_is_worse = ((dist[consider_me] > dist[n_minus_one_consider]) |
                             (sign[consider_me] != sign[n_minus_one_consider]))
            right_is_worse = ((dist[consider_me] > dist[n_plus_one_consider]) |
                              (sign[consider_me] != sign[n_plus_one_consider]))
            to_keep = left_is_worse & right_is_worse & (sign[consider_me] != 0)
            keep_me[consider_me] = to_keep
            del dist, sign
        a = a[keep_me,:]
        centers_per_point = centers_per_point[keep_me]
        del keep_me
        first_pass = False
    #
    # Finally, we have to shrink the results. We number each of the
    # points for a label, then only keep those whose indexes are
    # less than the count for their label.
    #
    within_label_index = np.array(list(range(result.shape[0])), np.int32)
    counts_per_point = result_counts[r_anti_indexes_per_point]
    result_indexes_per_point = result_index[r_anti_indexes_per_point] 
    within_label_index = (within_label_index - result_indexes_per_point)
    result = result[within_label_index < counts_per_point]
    return result, result_counts