docs: def get_decouplings_for_line(grid, line, settings, fids=None):
      """Compute the element pairs for regulazisation decoupling, given a grid
      object and a line, denoted by start and end coordinates, in the grid.
  
      Parameters
      ----------
      grid: a crtomo.grid.crt_grid object
      line: np.array/list with 4 entries: x1, y1, x2, x2, denoting start and end
            point
      fids: None, or: list with file ids for
          debug_nodes
  
  
  
      Returns
      -------
      neighbors: Mx2 array with the numbers of adjoining elements along the given
                 line
      """
      key = 'presort'
      # x/z coordinates of 'macro'-line
      x = [line[0], line[2]]
      y = [line[1], line[3]]
  
      # x/z coordinates of all nodes in the grid
      nx = grid.nodes[key][:, 1]
      ny = grid.nodes[key][:, 2]
      nxy = np.vstack((nx, ny)).T
  
      # shortest distance of all nodes to this line
      dist = np.array(distances(x, y, nx, ny))
  
      # set the epsilon environment
      if settings.get('eps', None) is not None:
          eps = settings['eps']
          print('Taking user supplied eps value: {0}'.format(eps))
      else:
          import scipy.spatial.distance as spdist
          pair_distances = spdist.pdist(nxy)
          print('automatic eps determination')
          print('minimal distance: {0}'.format(pair_distances.min()))
          eps = np.sqrt(
              pair_distances.min() ** 2 - (pair_distances.min() ** 2) / 4
          )
          print('estimated smallest distance to next node: {0}'.format(eps))
          eps *= 0.5
          print('final eps (0.5 * estimate): {0}'.format(eps))
          if eps == 0:
              raise Exception('eps == 0')
  
      # only consider nodes that lie in the eps environment
      indices = np.where(dist < eps)[0]
  
      # extract these nodes
      nodes = grid.nodes[key][indices, 0]
      nodes = indices
  
      nodes_full = np.hstack((
          grid.nodes[key][indices, :].squeeze(),
          np.atleast_2d(dist[indices]).T,
      ))
      if fids is not None:
          np.savetxt(
              fids[0],
              nodes_full.squeeze()
          )
          np.savetxt(
              fids[1],
              np.atleast_2d(dist).T,
          )
  
      elm_indices = []
      elm_nodes = []
      # find elements with two nodes in it
      for elmnr, element in enumerate(grid.elements):
          if len(np.intersect1d(element, nodes)) == 2:
              elm_indices.append(elmnr)
              elm_nodes.append(element)
  
      elms = np.array(elm_nodes)
  
      # find neighboring elements
      neighbors = []
      eta = settings['eta']
      # for each element on the line
      for index, elm in enumerate(elms):
          found_it = False
          for index1, elm1 in enumerate(elms):
              ints = np.intersect1d(elm, elm1)
              # only two nodes can be on a line
              if len(ints) == 2:
                  # this check ensures that we do not identify the
                  # boundary between two adjacent cells on the line
                  # as a decoupling line.
                  if np.all(dist[ints] < eps):
                      found_it = True
                      break
                  else:
                      pass
                      # print(
                      #     'distances of neighbor not correct:', dist[ints], eps
                      # )
              # probably the element itself
              elif len(ints) > 2:
                  # print('found more than two common nodes!')
                  # print('element nodes:', elm)
                  # print('element1 nodes:', elm1)
                  pass
          if found_it:
              nb = (elm_indices[index] + 1, elm_indices[index1] + 1, eta)
              # reversed neighbor
              # nb_rev = (nb[1], nb[0], nb[2])
              # if nb_rev not in neighbors:
              if True:
                  neighbors.append(
                      nb
                  )
          else:
              pass
              print('No neighbors found. Strange...')
  
      return np.array(neighbors)