docs: def feret_diameter(chulls, counts, indexes):
      '''Return the minimum and maximum Feret diameter for each object
      
      This function takes the convex hull data, as generated by convex_hull
      and returns the minimum and maximum Feret diameter for each convex hull.
      
      chulls    - an n x 3 matrix giving the label #, the i coordinate and the
                  j coordinate for each vertex in the convex hull
      counts    - the number of points in each convex hull
      
      returns two arrays, the minimum and maximum diameter per object.
      '''
      #
      # Method taken from Alsuwaiyel, "Algorithms, Design Techniques and
      # Analysis", World Scientific, 2003
      #
      # Given a vertex, p1 on the convex hull composed of points p1 to pm,
      # there is a line, pm <-> p1. Going from p(m-1) to p2, find the
      # distance from the point to the line until the next point in the
      # series has a lesser distance. This point forms an antipodal pair
      # with p1. Now take the line from p1 to p2 and continue to add points
      # as long as the point has a greater distance from the line than its
      # successor. The final point forms an antipodal point with both p1
      # and its succesor vertex.
      #
      # Note: as far as I (leek) can figure out, there is a bug here for
      # the case of a square:
      #
      #  0     1
      #
      #
      #  3     2
      #
      #  0 - 1 is generated first, then 0 - 2 and 0 - 3
      #  Then 1 - 3 is generated, missing 1 - 2. The author is only interested
      #  in max Feret diameter as are we.
      #
      # The minimum Feret diameter always has a parallel line of support
      # that runs through two antipodal pairs and the second that either
      # runs through an antipodal pair or a vertex.
      #
      if len(counts) == 0:
          return np.zeros(0)
      save_counts = counts
      counts = np.atleast_1d(counts)
      indexes = np.atleast_1d(indexes)
      chull_idx = np.hstack(([0], np.cumsum(counts[:-1])))
      save_chull_idx = chull_idx
      chull_pts = chulls[:,1:].astype(float)
      #
      # A list of antipode arrays
      #
      antipodes = []
      #
      # Get rid of degenerate cases
      #
      # A single point is its own antipode
      #
      antipodes.append(np.column_stack([chull_idx[counts==1]] * 2))
      #
      # A line has it's two points as antipodes of each other.
      #
      antipodes.append(np.column_stack([chull_idx[counts==2],
                                        chull_idx[counts==2] + 1]))
      chull_idx = chull_idx[counts > 2]
      counts = counts[counts > 2]
      if len(counts) > 0:
          #
          # Calculate distances from every vertex other than pm and p1
          # to the line between pm and p1
          #
          pm_idx = (chull_idx + counts - 1)
          pm = chull_pts[pm_idx,:]
          p1 = chull_pts[chull_idx,:]
          # There are counts - 2 vertices to be examined
          vvv = Indexes([counts-2])
          v_idx = chull_idx[vvv.rev_idx] + vvv.idx[0] + 1
          distances = distance2_to_line(chull_pts[v_idx],
                                        pm[vvv.rev_idx],
                                        p1[vvv.rev_idx])
          #
          # Find the largest distance, settling ties by picking the
          # one with the largest within-object index.
          #
          order = np.lexsort((vvv.idx[0], -distances, vvv.rev_idx))
          antipode_idx = order[vvv.fwd_idx] + 1 - vvv.fwd_idx
          vertex_idx = np.zeros(len(antipode_idx), int)
          while len(chull_idx) > 0:
              #
              # Add the current antipode / vertex pair
              #
              antipodes.append(np.column_stack((chull_idx + vertex_idx, 
                                                chull_idx + antipode_idx)))
              if len(chull_idx) == 0:
                  break
              #
              # Get the distance from the line to the current antipode
              #
              dc = distance2_to_line(chull_pts[chull_idx + antipode_idx, :],
                                     chull_pts[chull_idx + vertex_idx, :],
                                     chull_pts[chull_idx +vertex_idx + 1, :])
              #
              # Get the distance from the line to the next antipode
              #
              next_antipode_idx = antipode_idx + 1
              next_antipode_idx[next_antipode_idx == counts] = 0
              dn = distance2_to_line(chull_pts[chull_idx + next_antipode_idx, :],
                                     chull_pts[chull_idx + vertex_idx, :],
                                     chull_pts[chull_idx + vertex_idx + 1, :])
              #
              # If the distance to the next candidate is larger, advance
              #
              advance_antipode = dc <= dn
              antipode_idx[advance_antipode] += 1
              #
              # Otherwise, move to the next vertex.
              #
              vertex_idx[~ advance_antipode] += 1
              #
              # Get rid of completed convex hulls
              #
              to_keep = (antipode_idx < counts) & (vertex_idx != antipode_idx)
              antipode_idx = antipode_idx[to_keep]
              vertex_idx = vertex_idx[to_keep]
              chull_idx = chull_idx[to_keep]
              counts = counts[to_keep]
      antipodes = np.vstack(antipodes)
      l = chulls[antipodes[:,0], 0]
      
      pt1 = chulls[antipodes[:,0],1:]
      pt2 = chulls[antipodes[:,1],1:]
      distances = np.sum((pt1 - pt2) **2, 1)
                             
      max_distance = np.sqrt(fixup_scipy_ndimage_result(scind.maximum(distances, l, indexes)))
      #
      # For the minimum distance, we have to take the distance from each
      # antipode to any line between alternate successive antipodes.
      #
      counts = save_counts
      chull_idx = save_chull_idx
      indexer = Indexes(counts)
      # Get rid of degenerates
      #
      degenerates = antipodes[:,0] == antipodes[:,1]
      if np.all(degenerates):
          return np.zeros(len(indexes)), max_distance
      l = l[~degenerates]
      antipodes = antipodes[~degenerates]
      # We duplicate the list so each antipode appears as both a first and second
      #
      l = np.hstack((l,l))
      antipodes = np.vstack((antipodes, antipodes[:,::-1]))
      v2_idx = indexer.idx[0][antipodes[:,1]]
      #
      # Vertex number zero can be the vertex after the last one. We therefore
      # add a second vertex of "count" for each one where the second vertex is zero.
      #
      extras = antipodes[v2_idx == 0,:]
      antipodes = np.vstack([antipodes, extras])
      l = np.hstack((l, indexes[indexer.rev_idx[extras[:,0]]]))
      v2_idx = np.hstack((v2_idx, indexer.counts[0][indexer.rev_idx[extras[:,1]]]))
      #
      # We sort by label, first and second to order them for the search.
      #
      order = np.lexsort((v2_idx, antipodes[:,0], l))
      antipodes = antipodes[order,:]
      l = l[order]
      v2_idx = v2_idx[order]
      #
      # Now, we only want antipode pairs where:
      # * the label is the same as it's successor
      # * the first antipode is the same as it's successor
      # * the second antipode is one less than it's successor
      #
      v1_idx = indexer.idx[0][antipodes[:,0]]
      good = np.hstack(((l[1:] == l[:-1]) & 
                        (v1_idx[1:] == v1_idx[:-1]) &
                        (v2_idx[1:] == v2_idx[:-1]+1), [False]))
      if not np.any(good):
          return np.zeros(len(indexes)), max_distance
      antipodes = antipodes[good,:]
      l = l[good]
      v2_idx = v2_idx[good]
      v = chull_pts[antipodes[:,0],:]
      l0 = chull_pts[antipodes[:,1],:]
      #
      # The index of the second point in the line has to be done modulo
      #
      l1_idx = antipodes[:,1] + 1
      needs_modulo = v2_idx == indexer.counts[0][indexer.rev_idx[antipodes[:,1]]] - 1
      l1_idx[needs_modulo] -= indexer.counts[0][indexer.rev_idx[antipodes[needs_modulo,1]]]
      l1 = chull_pts[l1_idx,:]
      #
      # Compute the distances
      #
      distances = distance2_to_line(v, l0, l1)
      min_distance = np.sqrt(fixup_scipy_ndimage_result(scind.minimum(distances, l, indexes)))
      min_distance[np.isnan(min_distance)] = 0
      return min_distance, max_distance