def clip_segment_to_polygon(p1, p2, polygon_r, polygon_z) -> list[tuple[float, float]]:
    """The parameter intervals of a segment that lie inside a closed polygon.

    Parameters
    ----------
    p1 : sequence of float
        ``(R, Z)`` of the segment start [m].
    p2 : sequence of float
        ``(R, Z)`` of the segment end [m].
    polygon_r : array_like
        Polygon vertex major radii, closed or open [m].
    polygon_z : array_like
        Polygon vertex heights [m].

    Returns
    -------
    list of tuple of float
        ``(t_start, t_end)`` intervals of the segment parameter
        ``p = p1 + t (p2 - p1)``, ``0 <= t <= 1``, lying inside the polygon,
        in increasing order [-].

    Raises
    ------
    ValueError
        The polygon has fewer than three distinct vertices.

    Convention
    ----------
    Inside is decided by even-odd parity at each sub-interval's midpoint, so a
    non-convex polygon, such as a limiter with a divertor recess, yields
    several intervals. A segment grazing a vertex or running along an edge
    splits there with zero-length pieces dropped.

    Applicability
    -------------
    Machine-independent.

    Limitations
    -----------
    The polygon is a single loop; holes are not supported.

    Provenance
    ----------
    .. [1] The even-odd rule for point-in-polygon, as
       :class:`matplotlib.path.Path` applies it.
    """
    from matplotlib.path import Path

    pr = np.asarray(polygon_r, dtype=float).reshape(-1)
    pz = np.asarray(polygon_z, dtype=float).reshape(-1)
    if pr.size != pz.size:
        raise ValueError("polygon_r and polygon_z must have the same length")
    if pr.size and (pr[0] != pr[-1] or pz[0] != pz[-1]):
        pr = np.append(pr, pr[0])
        pz = np.append(pz, pz[0])
    if np.unique(np.column_stack((pr, pz)), axis=0).shape[0] < 3:
        raise ValueError("a polygon needs at least three distinct vertices")
    a = np.asarray(p1, dtype=float)
    d = np.asarray(p2, dtype=float) - a
    e0 = np.column_stack((pr[:-1], pz[:-1]))
    e = np.column_stack((pr[1:], pz[1:])) - e0
    denom = d[0] * e[:, 1] - d[1] * e[:, 0]
    rel = e0 - a
    with np.errstate(divide="ignore", invalid="ignore"):
        t = (rel[:, 0] * e[:, 1] - rel[:, 1] * e[:, 0]) / denom
        u = (rel[:, 0] * d[1] - rel[:, 1] * d[0]) / denom
    hit = (denom != 0) & (u >= 0.0) & (u <= 1.0) & (t > 0.0) & (t < 1.0)
    cuts = np.unique(np.concatenate(([0.0, 1.0], t[hit])))
    mids = 0.5 * (cuts[:-1] + cuts[1:])
    inside = Path(np.column_stack((pr, pz))).contains_points(a + mids[:, None] * d)
    intervals: list[tuple[float, float]] = []
    for lo, hi, keep in zip(cuts[:-1], cuts[1:], inside):
        if not keep or hi <= lo:
            continue
        if intervals and intervals[-1][1] == lo:
            intervals[-1] = (intervals[-1][0], float(hi))
        else:
            intervals.append((float(lo), float(hi)))
    return intervals

def project_emissivity(emissivity, operator) -> np.ndarray:
    """Line integrals of one or many ``(R, Z)`` fields through a path-length operator.

    Parameters
    ----------
    emissivity : array_like
        Field indexed ``(..., R, Z)``; leading axes, such as time, are kept
        [any, per metre].
    operator : scipy.sparse matrix or ndarray
        From :func:`build_line_integral_operator` for the same grid [m].

    Returns
    -------
    numpy.ndarray
        Line integrals shaped ``(..., n_chords)`` [any].

    Raises
    ------
    ValueError
        The field's last two axes do not hold as many cells as the operator
        has columns.

    Applicability
    -------------
    Machine-independent.

    Limitations
    -----------
    Exactly linear: whatever the field's unit, the result is that unit times
    metres. No detector response, filter transmission or etendue is applied.
    """
    field = np.asarray(emissivity, dtype=float)
    if field.ndim < 2:
        raise ValueError("emissivity must be indexed (..., R, Z)")
    lead = field.shape[:-2]
    cells = field.shape[-2] * field.shape[-1]
    if cells != operator.shape[1]:
        raise ValueError(f"the field has {cells} cells but the operator {operator.shape[1]} columns")
    flat = field.reshape(-1, cells)
    return np.asarray(operator @ flat.T).T.reshape(*lead, operator.shape[0])