def centers_of_labels(labels):
    '''Return the i,j coordinates of the centers of a labels matrix
    
    The result returned is an 2 x n numpy array where n is the number
    of the label minus one, result[0,x] is the i coordinate of the center
    and result[x,1] is the j coordinate of the center.
    You can unpack the result as "i,j = centers_of_labels(labels)"
    '''
    max_labels = np.max(labels)
    if max_labels == 0:
        return np.zeros((2,0),int)
    
    result = scind.center_of_mass(np.ones(labels.shape),
                                  labels,
                                  np.arange(max_labels)+1)
    result = np.array(result)
    if result.ndim == 1:
        result.shape = (2,1)
        return result
    return result.transpose()