code: def slh_associate(a_features, b_features, max_sigma=5):
      
      # 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)
docs: 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.