def calculate_convex_hull_areas(labels,indexes=None):
    """Calulculate the area of the convex hull of each labeled object
    
    labels - a label matrix
    indexes - None: calculate convex hull area over entire image
              number: calculate convex hull for a single label
              sequence: calculate convex hull for labels matching a sequence
                        member and return areas in same order.
    """
    if getattr(indexes,"__getitem__",False):
        indexes = np.array(indexes,dtype=np.int32)
    elif indexes is not None:
        indexes = np.array([indexes],dtype=np.int32)
    else:
        labels = labels !=0
        indexes = np.array([1],dtype=np.int32)
    hull, counts = convex_hull(labels, indexes)
    result = np.zeros((counts.shape[0],))
    #
    # Get rid of the degenerate cases
    #
    result[counts==1] = 1 # a single point has area 1
    if not np.any(counts >1):
        return result
    #
    # Given a label number "index_of_label" indexes into the result
    #
    index_of_label = np.zeros((hull[:,0].max()+1),int)
    index_of_label[indexes] = np.array(list(range(indexes.shape[0])))
    #
    # hull_index is the index into hull of the first point on the hull
    # per label
    #
    hull_index = np.zeros((counts.shape[0],),int)
    hull_index[1:] = np.cumsum(counts[:-1])
    #
    # A 2-point case is a line. The area of a line is its length * 1
    # and its length needs to be expanded by 1 because the end-points are
    # at the limits, not the ends.
    # 
    if np.any(counts==2):
        diff_2 = hull[hull_index[counts==2],1:]-hull[hull_index[counts==2]+1,1:]
        result[counts==2] = np.sqrt(np.sum(diff_2**2,1))+1
    if not np.any(counts>=3):
        return result
    #
    # Now do the non-degenerate cases (_nd)
    #
    counts_per_label = np.zeros((hull[:,0].max()+1),counts.dtype)
    counts_per_label[indexes] = counts
    hull_nd = hull[counts_per_label[hull[:,0]] >=3]
    counts_nd = counts[counts>=3]
    indexes_nd = indexes[counts>=3]
    index_of_label_nd = np.zeros((index_of_label.shape[0],),int)
    index_of_label_nd[indexes_nd] = np.array(list(range(indexes_nd.shape[0])))
    #
    # Figure out the within-label index of each point in a label. This is
    # so we can do modulo arithmetic when pairing a point with the next
    # when determining an edge
    #
    hull_index_nd = np.zeros((counts_nd.shape[0],),int)
    if hull_index_nd.shape[0] > 1:
        hull_index_nd[1:] = np.cumsum(counts_nd[:-1])
    index_of_label_per_pixel_nd = index_of_label_nd[hull_nd[:,0]]
    hull_index_per_pixel_nd = hull_index_nd[index_of_label_per_pixel_nd] 
    within_label_index = (np.array(list(range(hull_nd.shape[0]))) -
                          hull_index_per_pixel_nd)
    #
    # Find some point within each convex hull.
    #
    within_hull = np.zeros((counts_nd.shape[0],2))
    within_hull[:,0] = scind.sum(hull_nd[:,1],
                                 hull_nd[:,0],
                                 indexes_nd) / counts_nd
    within_hull[:,1] = scind.sum(hull_nd[:,2],
                                 hull_nd[:,0],
                                 indexes_nd) / counts_nd
    within_hull_per_pixel = within_hull[index_of_label_per_pixel_nd]
    #
    # Now, we do a little, slightly wierd fixup, arguing that the
    # edge of a pixel is +/- .5 of its coordinate. So we move the ones
    # left of center to the left by .5, right of center to the right by .5
    # etc.
    #
    # It works for a square...
    #
    # 2012-06-04 (thouis): Numpy trunk no longer allows unsafe type conversion
    # in in-place operations, which was what used to happen here
    # (adding/subtracting .5 from an integer array, which then just got
    # truncated).  So now it just adds 1 to the pixels right or below the
    # center.
    #
    # It would probably be better to do this in one of two other ways:
    # - Compute area with a circular ferret
    #   = polygon area + radius * perimeter + pi * radius^2
    # - Compute area with a square ferret
    #   = polygon area + 2 * radius * (horizontal_projection + vertical_projection) + 4*radius^2
    hull_nd[hull_nd[:,1] >= within_hull_per_pixel[:,0],1] += 1
    hull_nd[hull_nd[:,2] >= within_hull_per_pixel[:,1],2] += 1
    #
    # Finally, we go around the circle, computing triangle areas
    # from point n to point n+1 (modulo count) to the point within
    # the hull.
    #
    plus_one_idx = np.array(list(range(hull_nd.shape[0])))+1
    modulo_mask = within_label_index+1 == counts_nd[index_of_label_per_pixel_nd]
    plus_one_idx[modulo_mask] = hull_index_per_pixel_nd[modulo_mask]
    area_per_pt_nd = triangle_areas(hull_nd[:,1:],
                                    hull_nd[plus_one_idx,1:],
                                    within_hull_per_pixel)
    #
    # The convex area is the sum of these triangles
    #
    result[counts>=3] = scind.sum(area_per_pt_nd, hull_nd[:,0], indexes_nd)
    return result
def calculate_convex_hull_areas ( labels , indexes = None ) : if getattr ( indexes , "__getitem__" , False ) : indexes = np . array ( indexes , dtype = np . int32 ) elif indexes is not None : indexes = np . array ( [ indexes ] , dtype = np . int32 ) else : labels = labels != 0 indexes = np . array ( [ 1 ] , dtype = np . int32 ) hull , counts = convex_hull ( labels , indexes ) result = np . zeros ( ( counts . shape [ 0 ] , ) ) result [ counts == 1 ] = 1 if not np . any ( counts > 1 ) : return result index_of_label = np . zeros ( ( hull [ : , 0 ] . max ( ) + 1 ) , int ) index_of_label [ indexes ] = np . array ( list ( range ( indexes . shape [ 0 ] ) ) ) hull_index = np . zeros ( ( counts . shape [ 0 ] , ) , int ) hull_index [ 1 : ] = np . cumsum ( counts [ : - 1 ] ) if np . any ( counts == 2 ) : diff_2 = hull [ hull_index [ counts == 2 ] , 1 : ] - hull [ hull_index [ counts == 2 ] + 1 , 1 : ] result [ counts == 2 ] = np . sqrt ( np . sum ( diff_2 ** 2 , 1 ) ) + 1 if not np . any ( counts >= 3 ) : return result counts_per_label = np . zeros ( ( hull [ : , 0 ] . max ( ) + 1 ) , counts . dtype ) counts_per_label [ indexes ] = counts hull_nd = hull [ counts_per_label [ hull [ : , 0 ] ] >= 3 ] counts_nd = counts [ counts >= 3 ] indexes_nd = indexes [ counts >= 3 ] index_of_label_nd = np . zeros ( ( index_of_label . shape [ 0 ] , ) , int ) index_of_label_nd [ indexes_nd ] = np . array ( list ( range ( indexes_nd . shape [ 0 ] ) ) ) hull_index_nd = np . zeros ( ( counts_nd . shape [ 0 ] , ) , int ) if hull_index_nd . shape [ 0 ] > 1 : hull_index_nd [ 1 : ] = np . cumsum ( counts_nd [ : - 1 ] ) index_of_label_per_pixel_nd = index_of_label_nd [ hull_nd [ : , 0 ] ] hull_index_per_pixel_nd = hull_index_nd [ index_of_label_per_pixel_nd ] within_label_index = ( np . array ( list ( range ( hull_nd . shape [ 0 ] ) ) ) - hull_index_per_pixel_nd ) within_hull = np . zeros ( ( counts_nd . shape [ 0 ] , 2 ) ) within_hull [ : , 0 ] = scind . sum ( hull_nd [ : , 1 ] , hull_nd [ : , 0 ] , indexes_nd ) / counts_nd within_hull [ : , 1 ] = scind . sum ( hull_nd [ : , 2 ] , hull_nd [ : , 0 ] , indexes_nd ) / counts_nd within_hull_per_pixel = within_hull [ index_of_label_per_pixel_nd ] hull_nd [ hull_nd [ : , 1 ] >= within_hull_per_pixel [ : , 0 ] , 1 ] += 1 hull_nd [ hull_nd [ : , 2 ] >= within_hull_per_pixel [ : , 1 ] , 2 ] += 1 plus_one_idx = np . array ( list ( range ( hull_nd . shape [ 0 ] ) ) ) + 1 modulo_mask = within_label_index + 1 == counts_nd [ index_of_label_per_pixel_nd ] plus_one_idx [ modulo_mask ] = hull_index_per_pixel_nd [ modulo_mask ] area_per_pt_nd = triangle_areas ( hull_nd [ : , 1 : ] , hull_nd [ plus_one_idx , 1 : ] , within_hull_per_pixel ) result [ counts >= 3 ] = scind . sum ( area_per_pt_nd , hull_nd [ : , 0 ] , indexes_nd ) return result
Calulculate the area of the convex hull of each labeled object
    
    labels - a label matrix
    indexes - None: calculate convex hull area over entire image
              number: calculate convex hull for a single label
              sequence: calculate convex hull for labels matching a sequence
                        member and return areas in same order.