# Problem with definition of BSDF and radiance

When I was reading theory behind physical based rendering I noticed that definition of BSDF and radiance has some problems. For example BSDF of purely specular surfaces is zero almost everywhere and infinite in one point or radiance of directional light is zero for almost all directions except for one where it is again infinite.

This causes problems in rendering equation.

$$L_0(x, \omega_0)= \int_{S^2}{\rho(x, \omega_i,\omega_0)L(x, \omega_i)\,\mathrm{d}\sigma_\perp(\omega_i)}$$

For purely specular surface this integral has to be zero, from strict mathematical point of view. This is because the BSDF is zero (solid angle)-almost everywhere. You can argue that BSDF is infinite in one point and you have to take this point into account. But how do you know what is the reflectance of the surface at that point? From the infinity you can really tell. Furthermore from mathematical standpoint of view you cannot integrate infinite valued function, even if you could than the only sensible answers would be zero or infinity.

I know these are subtle problems and in practice can be solved with few ifs but I would like to have theory without holes. I believe that if you embrace these problems in theory than it helps you with placing those ifs in the right place.

How to deal with this problem?

Not sure entirely but as I can remember Erich Veach in his thesis address this problem only slightly and tries to get away with it by saying that BSDF and radiance are distributions. This is problematic, you cannot multiply two distributions, which is needed in rendering equation. For example when a light from directional light hits specular surface than you need to multiply two Dirac delta functions together.

The question is: Is there any work which reformulates rendering equation, BSDF and radiance in such a way that it does not suffer from mentioned problems?

What is the state of the art theory behind raytracing? Is is still Erich Veach's thesis?

(I'm only aware of work of Christian Lessig, but his work is still beyond my mathematical reach.)

I already have proposal how to deal with this problem. I define BSDF and radiance as measure. The basic idea works fine, but the whole theory of light transport needs to be redone to find out if it really works.

The main purpose of this question is to find out if someone else already did it, so I can read it and then focus my energy somewhere else.

• I don't see yet why the dirac deltas are a problem. While it is impossible to compute that with a computer using sampling (hence the ifs), the mathematics are clearly defined, right? Looking forward for somebody who can clarify that. Besides, in nature there are no real dirac deltas / infinity values since there are no perfect mirrors; but I guess that is another topic. Sep 10, 2015 at 10:08
• As long as BSDF and radiance can be dirac deltas at the same time than the rendering equation(as it is) is not mathematically well defined. Even if only BSDF would be allowed to be Dirac delta and we would formally treat BSDF as distribution than radiance needs to be smooth function in order to be mathematically 100% correct. But radiance under no way can be smooth function e.g. sharp shadows form discontinuities in radiance.
– tom
Sep 10, 2015 at 10:25
• Yes in reality you cannot have perfect mirrors, point and directional light sources or pin hole cameras. But we write programs where these things are and we need a theory which underpins them.
– tom
Sep 10, 2015 at 10:31
• @tom The rigorous mathematics that underlies delta distributions is measure theory. See the definition of the Dirac delta as a measure. I don't know off the top of my head of a work specifically treating the rendering equation in the context of measure theory, but pretty sure all this stuff is well-founded at the level of mathematical physics. Sep 10, 2015 at 18:28
• (Disclaimer: I Am Not A Rendering Person.) At any surface point $x$, the role of the BSDF is to act as a linear operator mapping the incident light $L_i$ to the exitant light $L_o$. Now there is no problem with $L_i$ and $L_o$ both being distributions, i.e. linear functions $D(\mathbb S^2)\to\mathbb R$, because addition and scalar multiplication of distributions is well-defined so they form a vector space. When $L_i$ and $L_o$ are functions we can represent the BSDF $\rho$ as a distribution, but if they're not we can still speak of linear transformations.
– user106
Sep 11, 2015 at 6:24

## Specular, Color, Fresnel

BSDF of purely specular surfaces is zero almost everywhere and infinite in one point [...] how do you know what is the reflectance of the surface at that point? From the infinity you can['t] really tell.

Let's consider a mirror BRDF as our example throughout (a specular BTDF and similar such objects are analogous). For this BRDF, we have (see e.g. PBRT 3 §8.2.2): $$f_r(\omega_i,x,\omega_o) ~:=~ F ~ \delta(\omega_i-\omega_r) ~/~ (n \cdot \omega_i)$$ where $$F$$ is the Fresnel coefficient and $$\delta$$ is the Dirac delta function. The full Fresnel term, because the complex-valued IOR varies per-wavelength, will produce a scalar that varies per-wavelength. Thus, different wavelengths of light reflect more or less out of the surface, and thus you get an albedo and a color. In practice, essentially everyone special-cases this to just be configurable to a constant color: $$f_r(\omega_i,x,\omega_o) ~=~ \rho_s ~ \delta(\omega_i-\omega_r) ~/~ (n \cdot \omega_i)$$ This is especially useful for making metallic mirrors, which need complex-valued (and therefore finicky and advanced) IORs to work. It's also useful because RGB rendering is unfortunately more popular than is spectral rendering.

We will see how this works out in a calculation below.

## Integrating Infinity, Functions, and Measures

Furthermore from mathematical standpoint of view you cannot integrate infinite valued function, even if you could than the only sensible answers would be zero or infinity.

We have to be careful and precise here. If you mean "a function whose values diverge to infinity", then we definitely can integrate it using the usual tools. (Formally; this is just an "improper Riemann integral of the second kind". For example, $$\int_0^1 x^{-1/2} d x = 2$$, which you can see plotted and solved or read a very short proof of.)

The Dirac delta function is a generalized function (which is why some people describe it variously as a distribution, a linear functional, or (incorrectly) not-a-function). In particular, it is not Lebesgue-integrable (because $$\delta(x) = 0$$ almost everywhere, yet $$\int \delta(x)~dx = 1 \neq 0 = \int 0~dx$$).

For the delta function, the informal definition is a sequence of functions that diverge to an infinite value, but that are each individually everywhere-finite. For example, as a series of Gaussians: $$\delta(x) := \lim_{b \rightarrow 0} ~ \frac{1}{|b|\sqrt{\pi}} ~ e^{-(x/b)^2}$$ That's probably what you're used to thinking of, but for full formality, we need to drop this idea and instead define it through measure theory. The following is the Dirac measure concentrated at $$\omega_r$$, which you can think of as "checking" that the perfect reflection direction $$\omega_r$$ is in the set $$A$$, a measurable subset of the unit ball $$B^3 \subset \mathbb{R}^3$$: $$\delta_{\omega_r}(A) := 1_A(\omega_r) = \begin{cases} 0 & \text{if } \omega_r \notin A \\ 1 & \text{if } \omega_r \in A \end{cases}$$

## A Lemma About the Dirac Delta

We can develop an important property from this from a certain property of Lebesgue integration. Let $$L(x,\omega_i)$$ be the usual radiance-at-$$x$$-in-direction-$$\omega_i$$ function from the rendering equation (treated as a scalar for your sanity) and define: $$M(t) := \delta_{\omega_r}(\{~ \omega_i \in B^3 ~|~ L(x,\omega_i)>t ~\})$$ The Lebesgue integral can be re-defined in terms of the improper Riemann integral: $$\int_{B^3} L ~ d\delta_{\omega_r} := \int_0^\infty M(t) ~ dt$$ Let's look at $$M(t)$$ a bit more closely. The set we pass in is the set of all directions with radiance at least $$t$$. The Dirac delta "checks" for the presence of the $$\omega_r$$ direction, so that's the only one that matters. Specifically, if the radiance from $$\omega_r$$ is at least $$L(x,\omega_r)>t$$, then $$\omega_r$$ is in the set we pass in and the measure outputs $$M(t)=1$$. However, if $$L(x,\omega_r)\leq t$$, then $$\omega_r$$ is not in the set we pass in, and so $$M(t)=0$$: $$M(t) = \begin{cases} 0 & \text{if } L(x,\omega_r) \leq t \\ 1 & \text{if } L(x,\omega_r) > t \end{cases}$$ The above integrand on the right-hand-side is $$1$$ exactly for $$0 and $$0$$ above that. The integral therefore has value $$L(x,\omega_r)$$, giving us the following property of the Dirac delta: $$\int_{B^3} L ~ d\delta_{\omega_r} = L(x,\omega_r)$$

## Rendering Equation with a Specular BRDF

Let's try to put the specular BRDF into the rendering equation, shall we? $$L(x,\omega_o) = \int_\Omega L(x,\omega_i) ~ (n \cdot \omega_i) ~ f_r(\omega_i,x,\omega_o) ~ d \omega_i$$ Substituting the Dirac delta into the BRDF, and substituting the BRDF into the rendering equation, we get: \begin{align} L(x,\omega_o) &= \int_\Omega L(x,\omega_i) ~ (n \cdot \omega_i) ~ \rho_s ~ \delta(\omega_i-\omega_r) ~/~ (n \cdot \omega_i) ~ d \omega_i\\ &= \rho_s ~ \int_\Omega L(x,\omega_i) ~ \delta(\omega_i-\omega_r) ~ d \omega_i \end{align} We apply the lemma proved above to get: $$L(x,\omega_o) = \rho_s ~ L(x,\omega_r)$$ That is, the reflected radiance out is a simple constant color times the radiance in, which is what people expect from their mirrors.

Note you should also special-case it and compute it this way. Perhaps it's obvious, but I can say from experience that in particular leaving the geometry terms uncancelled can lead to nasty numerical issues.

## Final Thoughts

The question is: Is there any work which reformulates rendering equation, BSDF and radiance in such a way that it does not suffer from mentioned problems? [...] The main purpose of this question is to find out if someone else already did it, so I can read it and then focus my energy somewhere else.

I'm not immediately aware of such a work, but at any rate if you've gotten this far, hopefully it is no longer necessary. Most of this post was dedicated to explaining where the specular BRDF came from and formally grounding the Dirac delta function.

Few people would care about these details because if you accept them on faith and just kindof gloss over the formal grounding, you can still prove what you want to prove (the rendering-equation section) and it Just Works. Still, it's good to know it can be shown at a deeper level. I hope this helps and apologize in advance for any mathematical typos.

• I upvoted the above, however it doesn't answer the question. With a Dirac delta in the brdf and in the light one gets a product of two Dirac deltas and @tom is correct that you get $\int \delta(t-y) \delta(t-x) dt = \delta(y-x)$, which is a problem since we do not want a Dirac delta as a result. The issue arises from "treating" both the radiance and brdf as Dirac deltas here. The proper way is to define a single dirac measure that is $1$ here for $\omega=\omega_r$. Intuitively if you use function notation this is like giving only one of those the Dirac delta. Nov 19 at 14:10
• Essentially your need a measure that depends on the brdf and the light sources. In Veach he has the path integral formulation which can easily handle this with a small modification - you just have to account for the light being a Dirac delta in the last throughput measure by making the last edge of the path have a measure that depends jointly on the light source and the brdf. This is required only for the last edge where both vertices are specular. E.g. if you have a laser and a mirror in front of it, then the only way to reach it with $1+$ bounces is through the mirror. Nov 19 at 14:32
• @lightxbulb I think you are incorrect; the question is basically asking about delta functions in light transport, with the examples of a specular BSDF or a directional light, not the combination of a specular BSDF and a directional light. I chose to address a specular BRDF as the example.¶ That said, for the combination, perhaps you can solve it the way you say. Myself, I'd probably just shrug and handle it in the next layer of pathracing's recursive integral—there will always be one, even it it's only the final measurement equation. Nov 19 at 15:33
• You should re-read the question, tom was aware of Veach's work (and what measures are, which you can see from his link), but: "This is problematic, you cannot multiply two distributions, which is needed in rendering equation." You cannot just "shrug and handle it in the next integral" - the measure is tied to this specific edge, this should be clear already from the two cosines in the throughput. You would just get nonsense by trying to transfer it as a Dirac delta to the film. To boot I can set a delta in the importance function too and then even this wrong suggestion is inapplicable. Nov 19 at 15:51
• Nice answer but I have to agree with @lightxbulb that it does not address the main issue i.e. multiplying two distributions in the rendering equation. The question is to some extent answered by the comment by user106 that BSDF should be a linear map transforming measures. Measure perspective is useful and it makes clear why for example SDS paths are hard to sample. The sampling usually assumes that the light intensity is absolutely continuous measure w. r. t. The standard borel measure. This is not the case for specular surfaces or ideal directional light sources.
– tom
Nov 20 at 15:36

The Rendering Equation

First thing first, the rendering equation in its solid angle formulation in vacuum is: $$L(x,\omega) = L_e(x,\omega) + \int_{S^2} f(x, \omega_i, \omega_o) L(r(x,\omega), -\omega)|\cos\theta_x| d\sigma(\omega).$$ You forgot the emission function $$L_e$$ in your question. The solution $$L:M\times S^2 \to \mathbb{R}$$ (where $$S^2=\{p\in\mathbb{R}^3\,:\, \|p\|_2=1\}$$ is just the 2-sphere) is fully determined by the collection $$(M, n, f, L_e)$$, where $$M\subset \mathbb{R}^3$$ is the set of points forming the scene surfaces, $$n:M\to S^2$$ is the function that assigns a normal to every point of $$M$$, $$f:M\times S^2 \times S^2 \to \mathbb{R}_+$$ is the BSDF (it represents the light scattering properties of the materials of your surfaces), and $$L_e:M\times S^2 \to \mathbb{R}_+$$ is the emitted radiance function (it is non-zero on your light sources). It is fairly natural that your geometry $$(M,n$$), the materials $$f$$, and the light sources $$L_e$$ will determine the steady state of the radiance $$L$$ in the scene.

The function $$r:M\times S^2 \to M\cup (\infty\cdot S^2)$$ is the ray-tracing function that returns the first intersection point $$r(x,\omega)$$ along the ray with origin $$x$$ and direction $$\omega$$, if there is no intersection you can formally claim that it returns the intersection with some sphere at infinity (e.g. your environment map). Typically the rendering equation is formulated in terms of the incident radiance $$L_i(x, \omega)$$, however, in vacuum we have $$L_i(x,\omega) = (GL)(x,\omega) = L(r(x,\omega),-\omega)$$ - the radiance that arrives at $$x$$ from direction $$\omega$$ is equal to the radiance leaving $$r(x,\omega)$$ in direction $$-\omega$$, which is also fairly intuitive. The term $$|\cos\theta_x| = |n(x)\cdot \omega|$$ is due to Lambert's cosine law, but it can be derived also by just integrating a differential $$2$$-density on the sphere (think surface flux integrals from calculus 3). The $$d\sigma(\omega)$$ is the area measure on the sphere, i.e. if you take the standard spherical coordinates parametrization $$x = \sin\theta\cos\phi$$, $$y=\sin\theta\sin\phi$$, $$z=\cos\theta$$ and your rotate them so that $$z$$ is aligned with $$n(x)$$, then you may rewrite $$d\sigma(\omega) = |\sin\theta| d\theta d\phi$$ and thus: $$L(x,\omega) = L_e(x,\omega) + \int_{0}^{2\pi}\int_0^{\pi} f(x, \omega_i, \omega_o) L(r(x,\omega), -\omega)|\cos\theta| |\sin\theta|d\theta d\phi.$$

Formal Solution

For brevity one typically writes the above in operator form $$L = L_e +TL$$ where $$T:(M\times S^2\to \mathbb{R})\to (M\times S^2\to \mathbb{R})$$ is the operator that stands for the integral. If we define the scattering operator $$K$$ and the propagation operator $$G$$: $$(Kg)(x,\omega) = \int_{S^2}f(x,\omega,\omega_o)g(x,\omega)d\sigma_{x}^{\perp}(\omega), \quad (Gh)(x,\omega) = h(r(x,\omega),-\omega),$$ then $$T=K\circ G$$. Since we have the identity $$L=L_e+TL$$ we can just plug it in itself once yielding $$L=L_e + T(L_e + TL) = L_e+TL_e+T^2L$$. Repeating this process until infinity yields the series $$L = \sum_{k=0}^{\infty}T^kL_e$$. Formally $$\sum_{k=0}^{\infty}T^k$$ is the solution operator for the rendering equation: \begin{align} L = L_e +TL \implies (I-T)L = L_e \implies L = (I-T)^{-1}L_e \implies (I-T)^{-1} = \sum_{k=0}^{\infty}T^k. \end{align} The latter is the Neumann series for $$T$$. While $$\|T\|<1$$ is a sufficient condition for convergence (i.e. if the brdf is energy conserving), it is not a necessary condition - for specific choices of $$L_e$$ the above series may be convergent even if $$(I-T)^{-1}$$ is unbounded. Note that $$f,L,L_e$$ are actual functions (i.e. $$T$$ maps functions to functions). In practice one may want to model specular reflection, point and directional lights etc., which according to "many treatments" cannot be modeled with classical functions and require distribution. Veach discusses this in his thesis (e.g. in some appendices on specular reflection and refraction), and even has a nice table 8.3 showing some examples.

Veach's Extension to Distributions

If you go to 5.A.1, Veach has a short note that $$\int_{\mathbb{R}} f(x)\delta(x-x_0)\,dx$$ is just symbolic notation for the evaluation functional $$\Lambda_{x_0}(f) = f(x_0)$$. This symbolic expression could also be framed as $$\int_{\mathbb{R}} f(x)\delta(x-x_0)dx := \int_{\mathbb{R}} f(x)\,d\delta_{x_0}(x)$$ where $$\delta_{x_0}$$ is the Dirac measure. Note however that $$\delta$$ is not the Radon-Nikodym derivative of this measure, i.e. the expression $$\int_{\mathbb{R}} f(x)\delta(x-x_0)dx$$ is still purely symbolic if not abuse of notation. When you write $$\int_{S^2} f(x,\omega,\omega_o)L_i(x,\omega)\,d\sigma_x^{\perp}(\omega) = \int_{S^2} \delta_{\sigma_x^{\perp}}(\omega-\omega_r)\rho(x,\omega,\omega_o)L_i(x,\omega)\,d\sigma_x^{\perp}(\omega),$$ the second expression is purely symbolic and means \begin{align} \int_{S^2} \delta_{\sigma_x^{\perp}}(\omega-\omega_r)\rho(x,\omega,\omega_o)L_i(x,\omega)\,d\sigma_x^{\perp}(\omega)&:= \int_{S^2} \rho(x,\omega,\omega_o)L_i(x,\omega)\,d\delta_{\omega_r}(\omega) \\ = \Lambda_{\omega_r}(\rho(x,\cdot,\omega_o)L_i(x,\cdot)) &= \rho(x,\omega_r,\omega_o) L_i(x,\omega_r). \end{align} The linear operator $$\Lambda_{\omega_r}$$ is not the BSDF itself, it is the operator $$T$$ when restricted to this specific specular vertex $$x$$, thus the BSDF cannot be a distribution even though Veach abuses notation (as is often done elsewhere) and calls it a distribution. What's really going on at a specular vertex $$p$$ is that the definition of $$K$$ at point $$p$$ is modified to $$(Kg)(p,\omega_o) = \int_{S^2}\rho(p,\omega,\omega_o)g(p,\omega)\,d\delta_{\omega_r}(\omega) = \rho(p,\omega_r,\omega_o)g(p,\omega_r).$$ Formally you could pretend that: $$f(x,\omega,\omega_r) = \rho(x,\omega,\omega_o)\frac{d\delta_{\omega_r}(\omega)}{d\sigma_x^{\perp}(\omega)} =\rho(x,\omega,\omega_o)\delta_{\sigma_x^{\perp}}(\omega-\omega_r).$$ The only issue is that the function/Radon-Nikodym derivative doesn't exist because the Dirac measure $$\delta_{\omega_r}$$ is not absolutely continuous, but the notation is convenient enough that everyone abuses it.

Extension to Dirac Delta Light Sources

As we saw above the whole "distribution" business is tantamount to modifying your scattering operator at specific points. There are many ways to achieve this - you could tie the measure to your bsdf, then a specular vertex you would have a "differential bsdf-measure" $$fd\delta$$ while at a normal vertex it would be $$fd\sigma_{x}^{\perp}$$. The other way is of course to abuse notation and pretend that $$\delta$$ is the Radon-Nikodym derivative. The most theoretically sound, but least intuitive way, is to probably not touch the bsdf and just modify the scattering operator. Until this point we were lucky that we had an operator which to convert to our functional. It may not be readily obvious what should be done if we want a point or directional light source however.

As an example consider just the $$0$$-th term $$L_e$$ of the Neumann expansion with only a directional light source. To get the contribution from zero bounces away we integrate $$L_e$$ directly on the film. For a specific pixel on the film we may assume that the contribution is computed as $$\int_{A}\left(\int_{S^2}W_e(x,\omega)L_e(r(x,\omega),-\omega)\,d\sigma_{x}^{\perp}(\omega)\right)\,dx.$$ Here $$A$$ is the film area and $$W_e$$ is the sensitivity function. The aperture is implicitly taken into account by the ray-tracing function. But for simplicity assume there is no aperture, there are now two options: you either treat $$L_e$$ as an actual Dirac delta $$L_e(y,\omega) = e(y,\omega)\delta(\omega-\omega_d)$$ in which case you modify the integral to: $$\int_{A}\left(\int_{S^2}W_e(x,\omega)e(r(x,\omega),-\omega)|n(x)\cdot\omega|\,d\delta_{-\omega_d}(\omega)\right)\,dx =\\= \int_{A}W_e(x,-\omega_d)e(r(x,-\omega_d),\omega_d)|n(x)\cdot\omega_d|\,dx.$$ You could however also require that $$W_e$$ is a Dirac delta, e.g. both w.r.t. position and direction. Once you do this the result will be a Dirac delta if you work with the abuse of notation approach (you will have one extra left over). This is not necessarily entirely wrong, if you assume that you really want a Dirac delta to mean that it carries infinite energy. As an example, if your sensor has both a positional (on the film) and directional Dirac delta, it would return the radiance of a normal light source when an importance ray hits it. However if that light source has a directional Dirac delta aligning with the one from the sensor, it is not that far fetched to say that you get an infinite contribution - and such an infinity that when integrated against a kernel would spit out the value at the kernel's specific location.

However, as far as I am aware we do not do this in computer graphics, e.g. when we have a point light source, directional light source, or even a laser, shining into the camera/film is not reproduced in the ray-tracers I am aware of. This means that they opted for treating this in a different way than the abuse of notation would suggest. E.g. if they have a specular light and then diffuse they treat it like a Dirac delta, however if they have a specular light and then a specular surface they do not use two Dirac deltas, instead they produce at most one Dirac delta, and that is only if the arguments of the two supposed Dirac deltas agree, so in this case $$\int g(\omega)\delta(\omega-\omega_1)\delta(\omega-\omega_2)\,d\omega = g(\omega_1)\delta_{\omega_1=\omega_2} \ne g(\omega_1)\delta(\omega_1-\omega_2)$$, note that in the middle expression the delta is a Kronecker delta and not a Dirac delta. The middle expression is consistent with how you do not visualize point/directional light sources on the film.

Conclusion

Ultimately how you decide to treat it is up to you, probably the most consistent way is using the abuse of notation, because then you can actually have multibounce specular contributions from a Dirac light source (e.g. "directional light" seen in a mirror), and it's also consistent as a limiting case of decreasing the solid angle emission of some light source and increasing its power - in the limit you get the directional Dirac delta, but the limiting emissions do not magically stop working after one bounce, so it's probably consistent to keep that for a Dirac delta light also. I think the main reason the Dirac delta nature of the light is not visualized on the film is because it is typically used as a cheap approximation of an actual physical light source, and the fact that it has infinite energy at a point would probably be perceived as an artefact on the film.

 Appendix: Some Details on the Dirac Delta

Here are some additional clarifications regarding the object (or I should rather say notational tool) that we call the Dirac delta. Essentially we wish for a function $$\delta : \mathbb{R} \to \mathbb{R}\cup \{\pm\infty\}$$ such that for any function $$g$$ we have $$\int_{\mathbb{R}} g(t)\delta(t-x)\,dt = g(x)$$. Unfortunately, such a function doesn't exist, however, there are objects (that are not functions) that behave in this way. I will list those roughly in order of most general to least general. The most general ones have the least assumptions on the class of functions $$g$$ for which they hold, while the least general ones have the most constraints on $$g$$.

Evaluation Functional $$\delta_x$$ (Constraint: Functions)

Let $$V$$ be the set of functions from some set $$X$$ to some set $$Y$$, i.e. $$g:X\to Y$$. Define $$\delta_x(g) = g(x)$$ - in this specific case we didn't assume anything about $$g$$ except that it is a function. If $$Y$$ is a vector space over a field $$\mathbb{F}$$, then $$V$$ is also a vector space over $$\mathbb{F}$$ with the function-function addition $$(g+h)(x) := g(x)+h(x)$$ and scalar-function multiplication $$(\alpha\cdot g)(x) := \alpha\cdot g(x)$$. Specifically if $$Y=\mathbb{F}$$ then $$\delta_x : V \to \mathbb{F}$$ is from the algebraic dual $$V'$$ of $$V$$, that is, the space of linear functions $$f^*:V\to \mathbb{F}$$ (also called functionals). Then $$\delta_x (g) = g(x)$$ is known as the evaluation functional (here we just assumed that $$g$$ are scalar functions). However you should notice that in this most general case we do not even have a notion of integration, i.e. we just have a notion of pointwise evaluation $$\delta_x(g) = g(x)$$ - the integral from our desired identity has been replaced by some map $$\delta_x$$.

Note also that $$V$$ being a set of functions is crucial. If I had taken $$V$$ to be $$L_2(X)$$ then I would be working with equivalence classes of functions and even the evaluation functional and the meaning of $$g(x)$$ would not be well-defined. One could of course try and pick out a representative from each equivalence class to which to apply the evaluation functional, but there is generally no well-defined way to do so for $$L_2$$. There is, however, a way to generalize this to some other spaces of equivalence classes, e.g. you could do so for Sobolev spaces $$H^{q}(\mathcal{X})$$ where $$q\in(0,\infty)$$ and $$2q> \dim X$$ (the dimension here is in the sense of manifold dimension).

Measure $$\delta_x$$ (Constraint: Measurable Functions)

If you have two measurable spaces $$(X,\mathcal{X})$$ and $$(Y,\mathcal{Y})$$ (those can be equal) and a function $$g:X\to Y$$ such that the preimage $$g^{-1}$$ of any set $$B$$ from $$\mathcal{Y}$$ is from $$\mathcal{X}$$: $$B\in \mathcal{Y} \implies g^{-1}(B)\in \mathcal{X}$$, then $$g$$ is termed a measurable function. If you additionally have a measure $$\mu$$ (e.g. the Lebesgue measure) then you can further limit $$\mathcal{X}$$ to be sets of Lebesgue-measurable sets, and we speak of $$\mu$$-measurable functions. For any $$A\in \mathcal{X}$$ you could now define the Dirac measure: $$\delta_x(A) = 1_A(x) = \begin{cases} 1, & x\in A, \\ 0, & x\not\in A.\end{cases}$$ Then you can define the integral of $$g$$ w.r.t. $$\delta_x$$ as: $$\int_{A} g(t)\,\delta_x(dt) = g(x)\delta_x(A)= \begin{cases} g(x), & x\in A, \\ 0, &x\not\in A, \end{cases}$$ Note that $$0$$ should be well-defined in $$Y$$ for this to make sense (all the practical examples you care about satisfy this however).

Now take $$X$$ to be some subset of $$\mathbb{R}^n$$ and $$Y\subseteq \mathbb{R}$$ then let $$\mathcal{X}$$ be the Lebesgue measurable subsets of $$X$$, and let $$\mu:\mathcal{X}\to\mathbb{R}_+$$ be the Lebesgue measure, then you can define integrals of the form $$\int_A g(t) \,\mu(dt)$$ for $$A\in\mathcal{X}$$. Each $$\mu$$-measurable function $$f$$ can then be interpreted as a density and it induces a measure $$\mu_f(A) = \int_A f(t)\,\mu(dt)$$. The Radon-Nikodym theorem goes in the opposite direction - it says that if you have two measures $$\mu$$ and $$\nu$$, then under some constraints (the measures being sigma-finite and $$\nu$$ being absolutely continuous w.r.t. $$\mu$$) there exists a $$\mu$$-measurable density function $$f$$ such that $$\nu(A) = \int_A f(t)\,\mu(dt)$$. Then $$f$$ is termed the Radon-Nikodym derivative of $$\nu$$ w.r.t. $$\mu$$: $$f(t) = \frac{\nu(dt)}{\mu(dt)}$$ and you can compute integrals w.r.t. $$\nu$$ by rewriting them as integrals w.r.t. $$\mu$$ (that is e.g. how you define $$n$$-volume integration on $$n$$-manifolds): $$\int_A g(t)\nu(dt) = \int_A g(t)f(t) \,\mu(dt).$$

In any case, using the above we see that if we have many compatible measures we can pick one common measure $$\mu$$ (e.g. the Lebesgue measure) and define all integrals w.r.t. it. It's what you do in practice, except when you encounter a Dirac delta, but then you modify your notation to look as if you're still using integration w.r.t. the Lebesgue measure for the Dirac delta. In our specific case it would be convenient if we could write $$\int_{A} g(t)\,\delta_x(dt) = \int_A g(t)\delta_{\mu}(t-x)\,d\mu(t) = g(x)\delta_{x}(A).$$ And that would have been possible if it weren't for the fact that $$\delta_x$$ is not absolutely continuous w.r.t. $$\mu$$. Notably $$\mu(\{x\}) = 0$$ but $$\delta_x(\{x\})=1$$ which violates absolute continuity. This means that such a measurable density function $$\delta_{\mu}(t-x)$$ does not exist/the Radon-Nikodym derivative $$\frac{d\delta_x}{d\mu}$$ does not exist. But since it is convenient notation, one goes ahead and defines the formal symbolic expression by analogy to the cases where the Radon-Nikodym derivative exists: $$\int_{A}f(t)\delta_{\mu}(t-x)\,\mu(dt) = \int_{A}f(t)\frac{\delta_x(dt)}{\mu(dt)}\,\mu(dt) := \int_A f(t)\,\delta_x(dt) = f(x)\delta_x(A).$$ The first two expressions are purely symbolic and equal to what's on the right by definition. Here $$\mu$$ is the Lebesgue measure or any other measure derived from the Lebesgue measure such as the solid angle or area measure used in the rendering equation. Note that you could try and define: $$\delta_{\mu}(x-t) := \frac{\delta_x(dt)}{\mu(dt)}, \text{ but } \frac{\delta_x(dt)}{\mu(dt)} \text{ doesn't exist.}$$ So the $$\delta_{\mu}$$ can only show up in an integral, and then it's just symbolic notation that translates to what we defined it as. Notably the following expression (potentially corresponding to a specular light and bsdf) is not defined: $$\int_{A} g(t)\delta_{\mu}(t-x)\delta_{\mu}(t-y)\,d\mu(dt) = \int_{A} g(t)\delta_{\mu}(t-y)\,\delta_x(dt) = ?,$$ simply because we never defined what $$\int_{A} \delta_{\mu}(t-y)\,\delta_x(dt)$$ means - remember that $$\delta_{\mu}(t-x)$$ doesn't actually exist and we only introduced notation for it w.r.t. an integral w.r.t. $$\mu$$ $$\mu$$. We could try and define a second "integral" of the above (remember that this is purely symbolic notation, so we're just mimicking what would have happened if the Radon-Nikodym derivative existed): \begin{align} I &=\int_B \int_{A} g(t)\delta_{\mu}(t-x)\delta_{\mu}(t-y)\,d\mu(dt)\,d\mu(dx) \\ &= \int_B\int_{A} g(t)\delta_{\mu}(t-y)\,\delta_x(dt)\,d\mu(dx) \\ &= \int_B g(x)\delta_{\mu}(x-y)1_A(x)\,\mu(dx) \\ &= \int_B g(x)1_A(x)\,\delta_y(dx) \\ &= g(y)1_A(y)1_B(y) = g(y)1_{A\cap B}(y).\end{align} So we could claim that the result is a differential measure element $$g(x)1_A(x)\,\delta_y(dx)$$ that expects to be integrated. Later on we can argue what this should mean if we treat $$dx$$ as a finite measure set and we shrink it down to zero - then naturally you would expect this result to blow up to infinity.

However, whenever you have a specular light + specular bsdf, another way to model it is as follows: \begin{align} I &=\int_{A} g(t)\delta_{\mu}(t-x)\delta_{\delta_x}(t-y)\,d\mu(dt) \\ &= \int_A g(t)\delta_{\delta_x}(t-y)\,d\delta_x(dt) \\ &= \int_A g(t)\delta_{\delta_x}(t-y)\frac{\delta_y(dt)}{\delta_x(dt)}\,d\delta_x(dt) \\ &= \int_A g(t) \delta_y(dt) = g(y)1_A(y). \end{align} Now this corresponds to what is done in graphics $$\delta_{\delta_x}$$, but I think that the first approach is more consistent with what would happen in reality. I believe that what we do in graphics is motivated more from practical considerations than trying to model things realistically. For example if you were to take a disk light source with directional Dirac deltas for each point of it, and if you have complex ideal mirrors after it until the film, which say is also such a Dirac delta directional disk, you would have to track an infinite amount of specular rays. You could of course try to Monte Carlo this, but if you say that these carry infinite energy then missing even just one can result in an infinite error at the film.

In either case, since the Radon-Nikodym derivative doesn't actually exist both approaches amount to just changing the integration measure associated with the operator $$K$$ at specific points of the scene - so in theory you are free to pick whichever you like.

Distribution $$\delta_x$$ (Constraint: $$C^{k}$$ Functions)

You may require $$\delta_{x}$$ to be from the continuous dual $$V^*$$ (you need to have a topology to define the continuous dual unlike for the algebraic dual). If $$V=C^{\infty}_c(X)$$ then the continuos dual $$V^*$$ is the space of distributions. A regular distribution is one arising as an integral of a locally $$L_1$$ integrable function, e.g. $$f^*(\phi) = \int f(t)\phi(t)\,dt$$. You can multiply a distribution $$D$$ by another distribution $$f^*$$, if the latter is a regular distribution arising from $$f\in C^{\infty}(X)$$. So if your radiance and BSDF were infinitely many times differentiable you could make distributions out of those and multiply those distributions, similar to how you would deal with functions. The infinitely many times differentiability condition can be relaxed if you consider distributions of finite order. What cannot be relaxed however is the fact that the Dirac delta distribution is not a regular distribution - meaning that there exists no functions $$\delta$$ such that $$\delta_x(\phi) = \int \delta(t-x) \phi(t)\,dt$$. So once again such notation would be purely symbolic. Furthermore while distribution multiplication is defined between a regular distribution (arising from a function in $$C^k$$) and a non-regular one, the multiplication between two arbitrary distributions is not well-defined (e.g. between two Dirac delta distributions). As before you can try to define what it should be in your case depending on how you want to model your problem.

Reproducing Kernel $$\delta_x$$ (Constraint: Reproducing Kernel Hilbert Spaces)

I should note that the Dirac delta exists as a function in the sense $$\int \delta(t-x) f(t)\,dt = f(x)$$ if you consider a reproducing kernel Hilbert space. Then the Dirac delta function could be identified with the reproducing kernel for the space. For example bandlimited continuous functions have the sinc kernel as a reproducing kernel. Unfortunately you cannot guarantee that your functions are bandlimited and continuous in graphics - any hard edge can produce a non-bandlimited or discontinuous function in your radiance, or e.g. an edge in the BSDF can do the same (e.g. checkerboard texture).

The Dirac Delta as a Limit

You can study the following limit: $$\lim_{\sigma\to 0}\int_{-\infty}^{\infty}f(t)\frac{1}{\sqrt{2\pi\sigma^2}}\exp\left(-\frac{(t-x)^2}{2\sigma^2}\right)\,dt = f(x).$$ As you can see it acts as a Dirac delta (for reasonable $$f$$). One may be tempted to move the limit inside the integral and state: $$\delta(t-x) = \lim_{\sigma\to 0}\exp\left(-\frac{(t-x)^2}{2\sigma^2}\right),$$ but you actually cannot do so as the two are not equivalent: $$\int_{-\infty}^{\infty}f(t)\lim_{\sigma\to 0}\frac{1}{\sqrt{2\pi\sigma^2}}\exp\left(-\frac{(t-x)^2}{2\sigma^2}\right)\,dt = 0 \ne f(x) = \lim_{\sigma\to 0}\int_{-\infty}^{\infty}f(t)\frac{1}{\sqrt{2\pi\sigma^2}}\exp\left(-\frac{(t-x)^2}{2\sigma^2}\right)\,dt.$$ However as before you could formally define the following symbolic expression: $$\int_{-\infty}^{\infty}f(t)\delta(t-x)\,dt := \lim_{\sigma\to 0}\int_{-\infty}^{\infty}f(t)\frac{1}{\sqrt{2\pi\sigma^2}}\exp\left(-\frac{(t-x)^2}{2\sigma^2}\right)\,dt.$$ Then you could choose to define the product of two Dirac deltas as: $$\int_{-\infty}^{\infty}\delta(t-x)\delta(t-y)\,dt := \lim_{\sigma_1,\sigma_2\to 0}\int_{-\infty}^{\infty}\frac{1}{2\pi|\sigma_1\sigma_2|}\exp\left(-\frac{(t-x)^2}{2\sigma_1^2}+\frac{(t-y)^2}{2\sigma^2_2}\right)\,dt.$$ My guess is that this integral should give you the following (I haven't checked it): $$\lim_{\sigma\to 0}\exp\left(-\frac{(t-x)^2}{2\sigma^2}\right) = \begin{cases} \infty, & x=y, \\ 0, &x\ne y\end{cases}.$$ Based on the above and based on the idea of representing a Dirac delta light source as a limit of physical light sources where you increase the power while decreasing the solid angle or/and area, you could argue that e.g. a directional light source should show up as infinity if the film also has a directional Dirac delta importance/sensor function. So probably the first approach discussed in the section on $$\delta_x$$ as a measure is the physically more relevant one, since it would agree with this limiting process.