docs: 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