First, the viewport size:
$$h_x = 2*d*tan(\theta_x/2)$$
$$h_y = 2*d*tan(\theta_y/2)$$

Each pixel (from your diagram) has the following size in the eye coordinate system:
$$W = h_x / (k-1)$$
$$H = h_y / (m-1)$$

Note that usually the field of view encompasses whole pixels and doesn't stop at the center of the edge pixels like your diagram shows.

If $P_c$ is the viewport center in pixel coordinates, let:

$$\vec {P'_{ij}} = P_{ij} - P_c$$

Therefore, $P_{ij}$ in eye space becomes:
$$
[P_{ij}]_e =
    \begin{pmatrix}
    (W, 0) * \vec {P'_{ij}} \\
    (0, H) * \vec {P'_{ij}} \\
    d \\
    \end{pmatrix}
$$

and in the standard base:
$$[P_{ij}]_{std} = (\vec b, \vec v, \vec t) * [P_{ij}]_e$$

Normalize the one in the base you want to get $r_{ij}$.

# Original answer

The last time I wanted to go from window coordinates (pixels) to eye coordinates, I followed the steps outlined here in the [OpenGL wiki][1]. Once you have the eye coordinates, you can just normalize the vector to obtain $r_{ij}$.

  [1]: https://www.khronos.org/opengl/wiki/Compute_eye_space_from_window_space