docs: def minimum_enclosing_circle(labels, indexes = None, 
                               hull_and_point_count = None):
      """Find the location of the minimum enclosing circle and its radius
      
      labels - a labels matrix
      indexes - an array giving the label indexes to be processed
      hull_and_point_count - convex_hull output if already done. None = calculate
      
      returns an Nx3 array organized as i,j of the center and radius
      Algorithm from 
      http://www.personal.kent.edu/~rmuhamma/Compgeometry/MyCG/CG-Applets/Center/centercli.htm
      who calls it the Applet's Algorithm and ascribes it to Pr. Chrystal
      The original citation is Professor Chrystal, "On the problem to construct
      the minimum circle enclosing n given points in a plane", Proceedings of
      the Edinburgh Mathematical Society, vol 3, 1884
      """
      if indexes is None:
          if hull_and_point_count is not None:
              indexes = np.array(np.unique(hull_and_point_count[0][:,0]),dtype=np.int32)
          else:
              max_label = np.max(labels)
              indexes = np.array(list(range(1,max_label+1)),dtype=np.int32)
      else:
          indexes = np.array(indexes,dtype=np.int32)
      if indexes.shape[0] == 0:
          return np.zeros((0,2)),np.zeros((0,))
  
      if hull_and_point_count is None:
          hull, point_count = convex_hull(labels, indexes)
      else:
          hull, point_count = hull_and_point_count
      centers = np.zeros((len(indexes),2))
      radii = np.zeros((len(indexes),))
      #
      # point_index is the index to the first point in "hull" for a label
      #
      point_index = np.zeros((indexes.shape[0],),int)
      point_index[1:] = np.cumsum(point_count[:-1]) 
      #########################################################################
      #
      # The algorithm is this:
      # * Choose a line S from S0 to S1 at random from the set of adjacent
      #   S0 and S1
      # * For every vertex (V) other than S, compute the angle from S0
      #   to V to S. If this angle is obtuse, the vertex V lies within the
      #   minimum enclosing circle and can be ignored.
      # * Find the minimum angle for all V.
      #   If the minimum angle is obtuse, stop and accept S as the diameter of 
      #   the circle.
      # * If the vertex with the minimum angle makes angles S0-S1-V and
      #   S1-S0-V that are acute and right, then take S0, S1 and V as the
      #   triangle within the circumscribed minimum enclosing circle.
      # * Otherwise, find the largest obtuse angle among S0-S1-V and
      #   S1-S0-V (V is the vertex with the minimum angle, not all of them).
      #   If S0-S1-V is obtuse, make V the new S1, otherwise make V the new S0
      #
      ##########################################################################
      #
      # anti_indexes is used to transform a label # into an index in the above array
      # anti_indexes_per_point gives the label index of any vertex
      #
      anti_indexes=np.zeros((np.max(indexes)+1,),int)
      anti_indexes[indexes] = list(range(indexes.shape[0]))
      anti_indexes_per_point = anti_indexes[hull[:,0]]
      #
      # Start out by eliminating the degenerate cases: 0, 1 and 2
      #
      centers[point_count==0,:]= np.NaN
      if np.all(point_count == 0):
          # Bail if there are no points in any hull to prevent
          # index failures below.
          return centers,radii
          
      centers[point_count==1,:]=hull[point_index[point_count==1],1:]
      radii[point_count < 2]=0
      centers[point_count==2,:]=(hull[point_index[point_count==2],1:]+
                                 hull[point_index[point_count==2]+1,1:])/2
      distance = centers[point_count==2,:] - hull[point_index[point_count==2],1:]
      radii[point_count==2]=np.sqrt(distance[:,0]**2+distance[:,1]**2)
      #
      # Get rid of the degenerate points
      #
      keep_me = point_count > 2
      #
      # Pick S0 as the first point in each label
      # and S1 as the second.
      #
      s0_idx = point_index.copy()
      s1_idx = s0_idx+1
      #
      # number each of the points in a label with an index # which gives
      # the order in which we'll get their angles. We use this to pick out
      # points # 2 to N which are the candidate vertices to S
      # 
      within_label_indexes = (np.array(list(range(hull.shape[0])),int) -
                              point_index[anti_indexes_per_point])
      
      while(np.any(keep_me)):
          #############################################################
          # Label indexing for active labels
          #############################################################
          #
          # labels_to_consider contains the labels of the objects which
          # have not been completed
          #
          labels_to_consider = indexes[keep_me]
          #
          # anti_indexes_to_consider gives the index into any vector
          # shaped similarly to labels_to_consider (for instance, min_angle
          # below) for every label in labels_to_consider.
          #
          anti_indexes_to_consider =\
              np.zeros((np.max(labels_to_consider)+1,),int)
          anti_indexes_to_consider[labels_to_consider] = \
              np.array(list(range(labels_to_consider.shape[0])))
          ##############################################################
          # Vertex indexing for active vertexes other than S0 and S1
          ##############################################################
          #
          # The vertices are hull-points with indexes of 2 or more
          # keep_me_vertices is a mask of the vertices to operate on
          # during this iteration
          #
          keep_me_vertices = np.logical_and(keep_me[anti_indexes_per_point],
                                               within_label_indexes >= 2)
          #
          # v is the vertex coordinates for each vertex considered
          #
          v  = hull[keep_me_vertices,1:]
          #
          # v_labels is the label from the label matrix for each vertex
          #
          v_labels = hull[keep_me_vertices,0]
          #
          # v_indexes is the index into "hull" for each vertex (and similarly
          # shaped vectors such as within_label_indexes
          #
          v_indexes=np.argwhere(keep_me_vertices).flatten().astype(np.int32)
          #
          # anti_indexes_per_vertex gives the index into "indexes" and
          # any similarly shaped array of per-label values
          # (for instance s0_idx) for each vertex being considered
          #
          anti_indexes_per_vertex = anti_indexes_per_point[keep_me_vertices]
          #
          # anti_indexes_to_consider_per_vertex gives the index into any
          # vector shaped similarly to labels_to_consider for each
          # vertex being analyzed
          #
          anti_indexes_to_consider_per_vertex = anti_indexes_to_consider[v_labels]
          #
          # Get S0 and S1 per vertex
          #
          s0 = hull[s0_idx[keep_me],1:]
          s1 = hull[s1_idx[keep_me],1:]
          s0 = s0[anti_indexes_to_consider_per_vertex]
          s1 = s1[anti_indexes_to_consider_per_vertex]
          #
          # Compute the angle S0-S1-V
          #
          # the first vector of the angles is between S0 and S1
          #
          s01 = (s0 - s1).astype(float)
          #
          # compute V-S1 and V-S0 at each of the vertices to be considered
          #
          vs0 = (v - s0).astype(float)
          vs1 = (v - s1).astype(float) 
          #
          #` Take the dot product of s01 and vs1 divided by the length of s01 *
          # the length of vs1. This gives the cosine of the angle between.
          #
          dot_vs1s0 = (np.sum(s01*vs1,1) /
                       np.sqrt(np.sum(s01**2,1)*np.sum(vs1**2,1)))
          angle_vs1s0 = np.abs(np.arccos(dot_vs1s0))
          s10 = -s01
          dot_vs0s1 = (np.sum(s10*vs0,1) /
                       np.sqrt(np.sum(s01**2,1)*np.sum(vs0**2,1)))
          angle_vs0s1 = np.abs(np.arccos(dot_vs0s1))
          #
          # S0-V-S1 is pi - the other two
          #
          angle_s0vs1 = np.pi - angle_vs1s0 - angle_vs0s1
          assert np.all(angle_s0vs1 >= 0)
          #
          # Now we find the minimum angle per label
          #
          min_angle = scind.minimum(angle_s0vs1,v_labels,
                                    labels_to_consider)
          min_angle = fixup_scipy_ndimage_result(min_angle)
          min_angle_per_vertex = min_angle[anti_indexes_to_consider_per_vertex]
          #
          # Calculate the index into V of the minimum angle per label.
          # Use "indexes" instead of labels_to_consider so we get something
          # with the same shape as keep_me
          #
          min_position = scind.minimum_position(angle_s0vs1,v_labels,
                                                indexes)
          min_position = fixup_scipy_ndimage_result(min_position).astype(int)
          min_position = min_position.flatten()
          #
          # Case 1: minimum angle is obtuse or right. Accept S as the diameter.
          # Case 1a: there are no vertices. Accept S as the diameter.
          #
          vertex_counts = scind.sum(keep_me_vertices,
                                    hull[:,0],
                                    labels_to_consider)
          vertex_counts = fixup_scipy_ndimage_result(vertex_counts)
          case_1 = np.logical_or(min_angle >= np.pi / 2,
                                    vertex_counts == 0)
                                     
          if np.any(case_1):
              # convert from a boolean over indexes_to_consider to a boolean
              # over indexes
              finish_me = np.zeros((indexes.shape[0],),bool)
              finish_me[anti_indexes[labels_to_consider[case_1]]] = True
              s0_finish_me = hull[s0_idx[finish_me],1:].astype(float)
              s1_finish_me = hull[s1_idx[finish_me],1:].astype(float)
              centers[finish_me] = (s0_finish_me + s1_finish_me)/2
              radii[finish_me] = np.sqrt(np.sum((s0_finish_me - 
                                                       s1_finish_me)**2,1))/2
              keep_me[finish_me] = False
          #
          # Case 2: all angles for the minimum angle vertex are acute 
          #         or right.
          #         Pick S0, S1 and the vertex with the
          #         smallest angle as 3 points on the circle. If you look at the
          #         geometry, the diameter is the length of S0-S1 divided by
          #         the cosine of 1/2 of the angle. The center of the circle
          #         is at the circumcenter of the triangle formed by S0, S1 and
          #         V.
          case_2 = keep_me.copy()
          case_2[angle_vs1s0[min_position] > np.pi/2] = False
          case_2[angle_vs0s1[min_position] > np.pi/2] = False
          case_2[angle_s0vs1[min_position] > np.pi/2] = False
          
          if np.any(case_2):
              #
              # Wikipedia (http://en.wikipedia.org/wiki/Circumcircle#Cartesian_coordinates)
              # gives the following:
              # D = 2(S0y Vx + S1y S0x - S1y Vx - S0y S1x - S0x Vy + S1x Vy)
              # D = 2(S0x (S1y-Vy) + S1x(Vy-S0y) + Vx(S0y-S1y)
              # x = ((S0x**2+S0y**2)(S1y-Vy)+(S1x**2+S1y**2)(Vy-S0y)+(Vx**2+Vy**2)(S0y-S1y)) / D
              # y = ((S0x**2+S0y**2)(Vx-S1x)+(S1x**2+S1y**2)(S0y-Vy)+(Vx**2+Vy**2)(S1y-S0y)) / D
              #
              ss0 = hull[s0_idx[case_2],1:].astype(float)
              ss1 = hull[s1_idx[case_2],1:].astype(float)
              vv  = v[min_position[case_2]].astype(float)
              Y = 0
              X = 1 
              D = 2*(ss0[:,X] * (ss1[:,Y] - vv[:,Y]) +
                     ss1[:,X] * (vv[:,Y]  - ss0[:,Y]) +
                     vv[:,X]  * (ss0[:,Y] - ss1[:,Y]))
              x = (np.sum(ss0**2,1)*(ss1[:,Y] - vv[:,Y]) +
                   np.sum(ss1**2,1)*(vv[:,Y]  - ss0[:,Y]) +
                   np.sum(vv**2,1) *(ss0[:,Y] - ss1[:,Y])) / D
              y = (np.sum(ss0**2,1)*(vv[:,X]  - ss1[:,X]) +
                   np.sum(ss1**2,1)*(ss0[:,X] - vv[:,X]) +
                   np.sum(vv**2,1) *(ss1[:,X] - ss0[:,X])) / D
              centers[case_2,X] = x
              centers[case_2,Y] = y
              distances = ss0-centers[case_2]
              radii[case_2] = np.sqrt(np.sum(distances**2,1))
              keep_me[case_2] = False
          #
          # Finally, for anybody who's left, for each of S0-S1-V and
          # S1-S0-V, for V, the vertex with the minimum angle,
          # find the largest obtuse angle. The vertex of this
          # angle (S0 or S1) is inside the enclosing circle, so take V
          # and either S1 or S0 as the new S.
          #
          # This involves a relabeling of within_label_indexes. We replace
          # either S0 or S1 with V and assign V either 0 or 1
          #
          if np.any(keep_me):
              labels_to_consider = indexes[keep_me]
              indexes_to_consider = anti_indexes[labels_to_consider]
              #
              # Index into within_label_indexes for each V with the
              # smallest angle
              #
              v_obtuse_indexes = v_indexes[min_position[keep_me]]
              angle_vs0s1_to_consider = angle_vs0s1[min_position[keep_me]]
              angle_vs1s0_to_consider = angle_vs1s0[min_position[keep_me]]
              #
              # Do the cases where S0 is larger
              #
              s0_is_obtuse = angle_vs0s1_to_consider > np.pi/2
              if np.any(s0_is_obtuse):
                  #
                  # The index of the obtuse S0
                  #
                  v_obtuse_s0_indexes = v_obtuse_indexes[s0_is_obtuse]
                  obtuse_s0_idx = s0_idx[indexes_to_consider[s0_is_obtuse]]
                  #
                  # S0 gets the within_label_index of the vertex
                  #
                  within_label_indexes[obtuse_s0_idx] = \
                      within_label_indexes[v_obtuse_s0_indexes]
                  #
                  # Assign V as the new S0
                  #
                  s0_idx[indexes_to_consider[s0_is_obtuse]] = v_obtuse_s0_indexes
                  within_label_indexes[v_obtuse_s0_indexes] = 0
              #
              # Do the cases where S1 is larger
              #
              s1_is_obtuse = np.logical_not(s0_is_obtuse)
              if np.any(s1_is_obtuse):
                  #
                  # The index of the obtuse S1
                  #
                  v_obtuse_s1_indexes = v_obtuse_indexes[s1_is_obtuse]
                  obtuse_s1_idx = s1_idx[indexes_to_consider[s1_is_obtuse]]
                  #
                  # S1 gets V's within_label_index and goes onto the list
                  # of considered vertices.
                  #
                  within_label_indexes[obtuse_s1_idx] = \
                      within_label_indexes[v_obtuse_s1_indexes]
                  #
                  # Assign V as the new S1
                  #
                  s1_idx[indexes_to_consider[s1_is_obtuse]] = v_obtuse_s1_indexes
                  within_label_indexes[v_obtuse_s1_indexes] = 1
      return centers, radii