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

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)