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
def minimum_enclosing_circle ( labels , indexes = None , hull_and_point_count = None ) : 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 = np . zeros ( ( indexes . shape [ 0 ] , ) , int ) point_index [ 1 : ] = np . cumsum ( point_count [ : - 1 ] ) 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 ] ] centers [ point_count == 0 , : ] = np . NaN if np . all ( point_count == 0 ) : 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 ) keep_me = point_count > 2 s0_idx = point_index . copy ( ) s1_idx = s0_idx + 1 within_label_indexes = ( np . array ( list ( range ( hull . shape [ 0 ] ) ) , int ) - point_index [ anti_indexes_per_point ] ) while ( np . any ( keep_me ) ) : labels_to_consider = indexes [ keep_me ] 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 ] ) ) ) keep_me_vertices = np . logical_and ( keep_me [ anti_indexes_per_point ] , within_label_indexes >= 2 ) v = hull [ keep_me_vertices , 1 : ] v_labels = hull [ keep_me_vertices , 0 ] v_indexes = np . argwhere ( keep_me_vertices ) . flatten ( ) . astype ( np . int32 ) anti_indexes_per_vertex = anti_indexes_per_point [ keep_me_vertices ] anti_indexes_to_consider_per_vertex = anti_indexes_to_consider [ v_labels ] 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 ] s01 = ( s0 - s1 ) . astype ( float ) vs0 = ( v - s0 ) . astype ( float ) vs1 = ( v - s1 ) . astype ( float ) 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 ) ) angle_s0vs1 = np . pi - angle_vs1s0 - angle_vs0s1 assert np . all ( angle_s0vs1 >= 0 ) 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 ] 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 ( ) 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 ) : 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 = 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 ) : 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 if np . any ( keep_me ) : labels_to_consider = indexes [ keep_me ] indexes_to_consider = anti_indexes [ labels_to_consider ] 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 ] ] s0_is_obtuse = angle_vs0s1_to_consider > np . pi / 2 if np . any ( s0_is_obtuse ) : v_obtuse_s0_indexes = v_obtuse_indexes [ s0_is_obtuse ] obtuse_s0_idx = s0_idx [ indexes_to_consider [ s0_is_obtuse ] ] within_label_indexes [ obtuse_s0_idx ] = within_label_indexes [ v_obtuse_s0_indexes ] s0_idx [ indexes_to_consider [ s0_is_obtuse ] ] = v_obtuse_s0_indexes within_label_indexes [ v_obtuse_s0_indexes ] = 0 s1_is_obtuse = np . logical_not ( s0_is_obtuse ) if np . any ( s1_is_obtuse ) : v_obtuse_s1_indexes = v_obtuse_indexes [ s1_is_obtuse ] obtuse_s1_idx = s1_idx [ indexes_to_consider [ s1_is_obtuse ] ] within_label_indexes [ obtuse_s1_idx ] = within_label_indexes [ v_obtuse_s1_indexes ] s1_idx [ indexes_to_consider [ s1_is_obtuse ] ] = v_obtuse_s1_indexes within_label_indexes [ v_obtuse_s1_indexes ] = 1 return centers , radii
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