The rendering equation for lambertian BRDF is (see Spherical integrals for info on how to solve):
Assume single directional light with direction (a, b) in spherical coordinates:

Where is a delta distribution with value in the direction of the light and everywhere else. This is written as (see Dirac delta for unit sphere):
Now we can solve the integral:
Note now that is precisely "N dot L" as you typically know it, so this matches typical realtime shading.

Note, I used Mathematica to solve the last integral. Wolfram isn't good at delta distributions.

FullSimplify[ Integrate[ Integrate[(DiracDelta[t - a, p - b])/Sin[a] Sin [t], {t, 0, pi/2}], {p, 0, 2 pi}], Assumptions -> {0 < a < pi/2, 0 < b < 2 pi}] = 1
Copy

#math #light-transport #pbr