def slh_associate(a_features, b_features, max_sigma=5):
    """
    An implementation of the Scott and Longuet-Higgins algorithm for feature
    association.

    This function takes two lists of features. Each feature is a
    :py:class:`MultivariateNormal` instance representing a feature
    location and its associated uncertainty.

    Args:
        a_features (list of MultivariateNormal)
        b_features (list of MultivariateNormal)
        max_sigma (float or int): maximum number of standard deviations two
            features can be separated and still considered "associated".

    Returns:
        (array): A Nx2 array of feature associations. Column 0 is the index into
        the a_features list, column 1 is the index into the b_features list.

    """
    # Compute proximity matrix
    proximity = _weighted_proximity(a_features, b_features)

    # Compute association matrix
    association_matrix = _proximity_to_association(proximity)

    # Now build up list of associations
    associations = []

    if association_matrix.shape[0] == 0:
        return np.zeros((0, 2))

    # Compute column-wise maxima
    col_max_idxs = np.argmax(association_matrix, axis=0)

    prox_threshold = np.exp(-0.5*max_sigma*max_sigma)

    # Find associations
    for row_idx, row in enumerate(association_matrix):
        if row.shape[0] == 0:
            continue

        # ... compute index of maximum element
        col_idx = np.argmax(row)

        # Is this row also the maximum in that column?
        if col_max_idxs[col_idx] == row_idx:
            prox = proximity[row_idx, col_idx]
            if prox > prox_threshold:
                associations.append((row_idx, col_idx))

    if len(associations) == 0:
        return np.zeros((0, 2))

    return np.vstack(associations)