Relativistic rendering: ray tracing photon paths affected by masses

The hard-easy way:

1. Calculate the stress–energy tensor (https://en.wikipedia.org/wiki/Stress%E2%80%93energy_tensor) for the distribution of mass you chose (assuming that distribution of mass is static - the mass/position/shapes(i.e. the mass density) do not change in time).
2. Introduce this in the Einstein Equations (https://en.wikipedia.org/wiki/Einstein_field_equations) and solve it for the components of the metric.
3. With the components of the metric, you can calculate the geodesics (https://en.wikipedia.org/wiki/Geodesics_in_general_relativity) and null paths. The path follow by light is a null geodesic. This will result in some non-linear differential equation that you might solve numerically, using the initial conditions.

The hard-hard way is almost the same, but now you also update the positions of the masses, as they orbit each other.

There might be some easy way to do this (e.g you might be able to directly calculate the components of the metric from the stress-energy tensor, for some simple mass distributions), but I am not aware of any.

Moreover, there is also the problem of the coordinate system you use. In some, the equations are easier to solve than in others. You might want to use that to simplify your life or allow for analytical simplifications.

You may also want to use more than just the mass, the angular momentum of the masses also changes the trajectories of the light rays in interesting ways.