How is Aperture and Pixel Size accounted for in calculating received power using ray tracing?

Question

NOTE: Deleting the question to rephrase it after continuing to read other responses and toy around with it:

I've written a path tracer and am now working on implementing a physically accurate way of simulating the actual sensor response for a given wavelength$^*$, given an exposure time, aperture size, pixel size, and quantum efficiency.

The current way I've thought about framing this is as follows (this algorithm is done for each pixel):

1. Path-trace many rays and evaluate the Bidirectional Reflectance Distribution Function (BRDF) (or BTDF/BSDF) to obtain the radiance along each ray received by the pixel. This part is done.
2. Take all of the rays received by a pixel, and average their radiance values to obtain a "received" radiance (units of $[W\cdot sr^{-1} \cdot m^{-2}]$) for the pixel.
3. Look up the solid angle field-of-view of a pixel (this is precalculated beforehand) and multiply that solid angle by the radiance to obtain the irradiance (units of $[W\cdot m^{-2}]$)
4. Multiply the irradiance by the area of the aperture to obtain the received power (units of $[W]$)
5. Multiply the power by the exposure time to obtain the energy (units of $[J]$)
6. Divide the energy by the photon energy (given by $E = h f$) for a given wavelength to obtain the photon count.
7. Multiply the photon count by the quantum efficiency and generate the signal.

Is this the correct way of calculating this? My primary concern is in steps 2 and 3 as I am not terribly comfortable with solid angles. My rationale is that the solid angle of the pixel's field-of-view determines how much of the scene is actually contributing light to the pixel, and we're approximating that part of the scene by casting discrete rays which we can then average. And this should naturally capture the effect that increasing the pixel size correspondingly increases the pixel's solid angle thus an increase in the irradiance at the aperture that contributes to the pixel, and thus increased response. But like I said, the nature of solid angles has me tripped up, particularly because it is often said that the radiance along a ray is reversible, and so I'm not sure "which" solid angle to actually be using.

---

$^*$ Currently I'm able to assume a single wavelength at a time, however in the future if I were to extend this to handle a broad-spectrum, I would modify steps 1 and 2 to calculate the spectral radiance with units of $[W\cdot sr^{-1} \cdot m^{-2} \cdot Hz^{-1}]$. Then when we continue on to step 3 we would calculate a spectral irradiance with units of $[W\cdot m^{-2} \cdot Hz^{-1}]$, and then step 4 would produce a spectral flux with units of $[W\cdot Hz^{-1}]$. Then we calculate the energy of photons for each frequency in our spectrum to obtain the photon counts for each frequency and then convolve the resulting energy spectrum with the quantum efficiency spectrum.

Answer 1

To understand whether your steps are correct you should compare to the mathematical formulation.

Step 1 is supposedly solving the rendering equation for some collection of initial rays

$$L_o(x_0, \omega_o) = L_e(x, \omega_o) + \int_{\mathcal{H}^2(x_0, N_{x_0})} f(\omega_o, x_0, \omega_i)L_i(x_0,\omega_i)(N_{x_0}\cdot \omega_i)\,d\omega_i.$$

Here $L_o(x_0,\omega_o)$ is the outgoing radiance function, i.e. the radiance leaving some surface point $x_0$ in direction $\omega_o$. $L_e(x_0,\omega_o)$ is the emitted radiance function, i.e. the radiance emitted from point $x_0$ in direction $\omega_o$. Next we have an integral over the hemisphere at point $x_0$ whose up direction is aligned with the normal $N_{x_0}$ at $x_0$. You can formalize this set as $\mathcal{H}^2(x_0,N_{x_0}) = \{z\in\mathbb{R}^3\,:\, \|z-x_0\|=1, \, (z-x_0)\cdot N_{x_0}\geq 0\}$, in practice one parametrizes it using spherical coordinates. $f$ is the BRDF. $L_i$ is the incident radiance function, i.e. $L_i(x_0,\omega_i) = L_o(r(x_0, \omega_i), -\omega_i)$, where $r$ is the ray-tracing function. Finally $(N_{x_0}\cdot \omega_i) = \cos\theta_i$ is the Lambertian cosine term that accounts for the angle between the ray and the surface.

For simplicity the above can be rewritten using linear operator notation as $L = L_e + TL$, from where one may get the Neumann expansion $L = \sum_{k=0}^{\infty}T^kL_e$, which is essentially a sum of increasingly dimensional integrals. You typically apply Monte Carlo to get an approximate solution of the latter.

The question is also how you choose the rays. Since you mentioned that you have an aperture, a film, and an exposure time, let us model the aperture as some set $\mathcal{A}\subset\mathbb{R}^3$, the film as $\mathcal{F}\subset\mathbb{R}^3$, and the exposure time as an interval $[t_0,t_1]$. Both of $\mathcal{F}$ and $\mathcal{A}$ are two-dimensional manifolds, and usually $\mathcal{A}$ is chosen to be a disk some distance in front of the film $\mathcal{F}$, and $\mathcal{F}$ is chosen to be a rectangle (in order to match the screen). The response at a given pixel (i,j) can be modeled as a weighted integration over the arriving radiance on the film, with some sensitivity function $w_{ij}$: $$I_{ij} = \int_{t_0}^{t_1}\int_{\mathcal{F}}\int_{\mathcal{A}_{\mathcal{H}^2(x_{-1}, N_{x_{-1}})}}w_{ij}(x_{-1}, -\omega_o)L_i(x_{-1}, -\omega_o))(N_{x_{-1}}\cdot (-\omega_o))\, d(-\omega_o)dx_{-1}dt.$$

Here $\mathcal{A}_{\mathcal{H}^2(x_{-1}, N_{x_{-1}})}$ is the solid angle set of the aperture as seen from pixel $x_{-1}$. Here I have assumed that nothing depends on time (i.e. brdf, scene geometry, lights, aperture, film, sensitivity, all stay the same). With this assumption you can pop out the time integral as a $(t_1-t_0)$ multiplicative term as you describe in Step 5.

Note that physically the aperture is typically in front of the sensor, so I have modeled that instead of what you did. One could of course set the aperture behind the sensor in CG for convenience, and then you would get the solid angle of the film pixel's set as you do. I will however stick to the physical formulation of the aperture being in front of the sensor. We also need to get rid of that cosine dot product and want to integrate only over the pixel and not the whole film, which can be done by setting the sensitivity $w_{ij}(x_{-1},-\omega_o) = \frac{\boldsymbol{1}_{\mathcal{P}_{ij}}(x_{-1})}{N_{x_{-1}}\cdot (-\omega_o)}$. Here $\mathcal{P}_{ij}\subset \mathcal{F}$ is the set of points of the pixel on the film, and $\boldsymbol{1}_{\mathcal{X}}$ is the indicator function (it is one on the set, and zero outside of it). You want to average the radiance (step 2) and then multiply by the solid angle measure of the aperture as seen by a specific point $x_{-1}$ (step 3 in your aperture/pixel reversed scenario has the solid angle of the pixel), and the area of the pixel (step 4 has the area of the aperture). This flip of the aperture and film pixel is not really an issue as noted, but I'll stick to the physical convention (i.e. the reverse of yours). The main issue is that averaging and then multiplying by the solid angle corresponding to the the aperture as seen by the whole pixel, vs by a single point is slightly wrong. Note that $\mathcal{A}_{\mathcal{H}^2(x_{-1}, N_{x_{-1}})}$ is defined per point, and it varies for the different points on the film's pixel, so you cannot precompute one value for the pixel and call it a day, unless the pixel consists only of a single point - then it would be valid.

There is another formulation that is amenable to a slight modification of your idea, however it uses the area formulation instead of the solid angle one. Namely

$$I_{ij} = \int_{t_0}^{t_1}\int_{\mathcal{F}}\int_{\mathcal{A}}w_{ij}(x_{-1}, x_{-1/2})L_i(x_{-1}, -\omega_o))G(x_{-1},x_{-1/2})\, d x_{-1/2}dx_{-1}dt.$$

Here $x_{-1/2}$ is a point on the aperture, and $-\omega_o = \frac{x_{-1/2}-x_{-1}}{\|x_{-1/2}-x_{-1}\|}$, and $G(x_{-1},x_{-1/2}) = \frac{(N_{x_{-1}}\cdot(-\omega_o))(N_{x_{-1/2}}\cdot \omega_o)}{\|x_{-1/2}-x_{-1}\|^2}$. Choose

$$w_{ij}(x_{-1}, x_{-1/2}) = \frac{\boldsymbol{1}_{\mathcal{P}_{ij}}(x_{-1})}{G(x_{-1},x_{-1/2})},$$

where $\mathcal{P}_{ij}\subseteq \mathcal{F}$ is the set of pixel $(i,j)$ on the film, and $\boldsymbol{1}_{\mathcal{X}}$ is the indicator function. Then you get what you want, with the small change that you multiply the average by the areas $|\mathcal{A}|$ and $|\mathcal{P}_{ij}|$ , without involving solid angles of those.

A general remark I have is that if you want to convert between different units then you need integration and differentiation. The multiplication and division is a special case of this, which is valid only in very specific scenarios. My takeaway from your explanation is that this was the main issue with your understanding of the problem.