<h2>Specular, Color, Fresnel</h2>
<blockquote>
<p>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.</p>
</blockquote>
<p>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. <a href="https://pbr-book.org/3ed-2018/Reflection_Models/Specular_Reflection_and_Transmission#SpecularReflection" rel="nofollow noreferrer">PBRT 3 §8.2.2</a>):
<span class="math-container">$$
    f_r(\omega_i,x,\omega_o) ~:=~ F ~ \delta(\omega_i-\omega_r) ~/~ (n \cdot \omega_i)
$$</span>
where <span class="math-container">$F$</span> is the Fresnel coefficient and <span class="math-container">$\delta$</span> is the <a href="https://en.wikipedia.org/wiki/Dirac_delta_function" rel="nofollow noreferrer">Dirac delta function</a>.  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, <em>essentially everyone</em> special-cases this to just be configurable to a constant color:
<span class="math-container">$$
    f_r(\omega_i,x,\omega_o) ~=~ \rho_s ~ \delta(\omega_i-\omega_r) ~/~ (n \cdot \omega_i)
$$</span>
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.</p>
<p>We will see how this works out in a calculation below.</p>
<hr />
<h2>Integrating Infinity, Functions, and Measures</h2>
<blockquote>
<p>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.</p>
</blockquote>
<p>We have to be careful and precise here.  If you mean &quot;a function whose values <em>diverge to</em> infinity&quot;, then we definitely <em>can</em> integrate it using the usual tools.  (Formally; this is just an &quot;improper Riemann integral of the second kind&quot;.  For example, <span class="math-container">$\int_0^1 x^{-1/2} d x = 2$</span>, which you can <a href="https://www.wolframalpha.com/input?i=integrate+1%2Fsqrt%28x%29+from+0+to+1" rel="nofollow noreferrer">see plotted and solved</a> or <a href="https://en.wikipedia.org/wiki/Improper_integral#Examples" rel="nofollow noreferrer">read a very short proof of</a>.)</p>
<p>The Dirac delta function is a <em>generalized</em> function (which is why some people describe it variously as a distribution, a linear functional, or (incorrectly) not-a-function).  In particular, it <a href="https://en.wikipedia.org/wiki/Lebesgue_integration#Basic_theorems_of_the_Lebesgue_integral" rel="nofollow noreferrer">is not Lebesgue-integrable</a> (because <span class="math-container">$\delta(x) = 0$</span> almost everywhere, yet <span class="math-container">$\int \delta(x)~dx = 1 \neq 0 = \int 0~dx$</span>).</p>
<p>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:
<span class="math-container">$$
    \delta(x) := \lim_{b \rightarrow 0} ~ \frac{1}{|b|\sqrt{\pi}} ~ e^{-(x/b)^2}
$$</span>
That's probably what you're used to thinking of, but for full formality, we need to drop this idea and instead <a href="https://en.wikipedia.org/wiki/Dirac_measure" rel="nofollow noreferrer">define it through measure theory</a>.  The following is the Dirac measure concentrated at <span class="math-container">$\omega_r$</span>, which you can think of as &quot;checking&quot; that the perfect reflection direction <span class="math-container">$\omega_r$</span> is in the set <span class="math-container">$A$</span>, a measurable subset of the unit ball <span class="math-container">$B^3 \subset \mathbb{R}^3$</span>:
<span class="math-container">$$
    \delta_{\omega_r}(A) := 1_A(\omega_r) = \begin{cases}
        0 &amp; \text{if } \omega_r \notin A \\
        1 &amp; \text{if } \omega_r \in A
    \end{cases}
$$</span></p>
<hr />
<h2>A Lemma About the Dirac Delta</h2>
<p>We can develop an important property from this from <a href="https://en.wikipedia.org/wiki/Lebesgue_integration#Via_improper_Riemann_integral" rel="nofollow noreferrer">a certain property of</a> Lebesgue integration.  Let <span class="math-container">$L(x,\omega_i)$</span> be the usual radiance-at-<span class="math-container">$x$</span>-in-direction-<span class="math-container">$\omega_i$</span> function from the rendering equation (treated as a scalar for your sanity) and define:
<span class="math-container">$$
    M(t) := \delta_{\omega_r}(\{~ \omega_i \in B^3 ~|~ L(x,\omega_i)&gt;t ~\})
$$</span>
The Lebesgue integral can be re-defined in terms of the improper Riemann integral:
<span class="math-container">$$
    \int_{B^3} L ~ d\delta_{\omega_r} := \int_0^\infty M(t) ~ dt
$$</span>
Let's look at <span class="math-container">$M(t)$</span> a bit more closely.  The set we pass in is the set of all directions with radiance at least <span class="math-container">$t$</span>.  The Dirac delta &quot;checks&quot; for the presence of the <span class="math-container">$\omega_r$</span> direction, so that's the only one that matters.  Specifically, if the radiance from <span class="math-container">$\omega_r$</span> is at least <span class="math-container">$L(x,\omega_r)&gt;t$</span>, then <span class="math-container">$\omega_r$</span> is in the set we pass in and the measure outputs <span class="math-container">$M(t)=1$</span>.  However, if <span class="math-container">$L(x,\omega_r)\leq t$</span>, then <span class="math-container">$\omega_r$</span> is not in the set we pass in, and so <span class="math-container">$M(t)=0$</span>:
<span class="math-container">$$
    M(t) = \begin{cases}
        0 &amp; \text{if } L(x,\omega_r) \leq t \\
        1 &amp; \text{if } L(x,\omega_r) &gt; t
    \end{cases}
$$</span>
The above integrand on the right-hand-side is <span class="math-container">$1$</span> exactly for <span class="math-container">$0&lt;t&lt;L(x,\omega_r)$</span> and <span class="math-container">$0$</span> above that.  The integral therefore has value <span class="math-container">$L(x,\omega_r)$</span>, giving us the following property of the Dirac delta:
<span class="math-container">$$
    \int_{B^3} L ~ d\delta_{\omega_r} = L(x,\omega_r)
$$</span></p>
<hr />
<h2>Rendering Equation with a Specular BRDF</h2>
<p>Let's try to put the specular BRDF into the rendering equation, shall we?
<span class="math-container">$$
    L(x,\omega_o) = \int_\Omega L(x,\omega_i) ~ (n \cdot \omega_i) ~ f_r(\omega_i,x,\omega_o) ~ d \omega_i
$$</span>
Substituting the Dirac delta into the BRDF, and substituting the BRDF into the rendering equation, we get:
<span class="math-container">\begin{align}
    L(x,\omega_o) &amp;= \int_\Omega L(x,\omega_i) ~ (n \cdot \omega_i) ~
        \rho_s ~ \delta(\omega_i-\omega_r) ~/~ (n \cdot \omega_i) ~
        d \omega_i\\
                  &amp;= \rho_s ~ \int_\Omega L(x,\omega_i) ~
        \delta(\omega_i-\omega_r) ~
        d \omega_i
\end{align}</span>
We apply the lemma proved above to get:
<span class="math-container">$$
    L(x,\omega_o) = \rho_s ~ L(x,\omega_r)
$$</span>
That is, the reflected radiance out is a simple constant color times the radiance in, which is what people expect from their mirrors.</p>
<p>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.</p>
<hr />
<h2>Final Thoughts</h2>
<blockquote>
<p>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.</p>
</blockquote>
<p>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.</p>
<p>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.</p>

CC BY-SA 4.0
2023-11-19T14:10:40.207Z
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.
lightxbulb
CC BY-SA 4.0
2023-11-19T14:32:55.233Z
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.
lightxbulb
CC BY-SA 4.0
2023-11-19T15:33:52.713Z
@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.
geometrian
CC BY-SA 4.0
2023-11-19T15:51:12.583Z
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.
lightxbulb
CC BY-SA 4.0
2023-11-20T15:36:24.433Z
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
CC BY-SA 4.0
2023-11-20T15:38:51.340Z
Oh this is already described in the blog post I linked ... well it has been 8 years :)
tom
CC BY-SA 4.0
2023-11-20T17:41:22.583Z
@tom I actually don't get user106's comment. He says that the bsdf is a linear map, but the bsdf is not a linear map eating $L_i$ and spitting out $L_o$. In vacuum $L_i(x, \omega) = L_o(r(x,\omega),-\omega)$. The rendering equation is an equation specifically because $L\equiv L_o$ shows up on both sides $L = L_e +TL$. The formal solution is $L = (I-T)^{-1}L_e = \sum_{k=0}^{\infty}T^k L_e$. If we pick two specular/Dirac vertices at the end it isn't really an issue either imo, you just have to define the distribution mult differently, e.g. if you use a joint measure it's fine.
lightxbulb
CC BY-SA 4.0
2023-11-19T06:08:12.407Z
2023-11-19T06:08:12.407Z
geometrian
<p><strong>The Rendering Equation</strong></p>
<p>First thing first, the rendering equation in its solid angle formulation in vacuum is:
<span class="math-container">$$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).$$</span>
You forgot the emission function <span class="math-container">$L_e$</span> in your question.  The solution <span class="math-container">$L:M\times S^2 \to \mathbb{R}$</span> (where <span class="math-container">$S^2=\{p\in\mathbb{R}^3\,:\, \|p\|_2=1\}$</span> is just the 2-sphere) is fully determined by the collection <span class="math-container">$(M, n, f, L_e)$</span>, where <span class="math-container">$M\subset \mathbb{R}^3$</span> is the set of points forming the scene surfaces, <span class="math-container">$n:M\to S^2$</span> is the function that assigns a normal to every point of <span class="math-container">$M$</span>, <span class="math-container">$f:M\times S^2 \times S^2 \to \mathbb{R}_+$</span> is the BSDF (it represents the light scattering properties of the materials of your surfaces), and <span class="math-container">$L_e:M\times S^2 \to \mathbb{R}_+$</span> is the emitted radiance function (it is non-zero on your light sources). It is fairly natural that your geometry <span class="math-container">$(M,n$</span>), the materials <span class="math-container">$f$</span>, and the light sources <span class="math-container">$L_e$</span> will determine the steady state of the radiance <span class="math-container">$L$</span> in the scene.</p>
<p>The function <span class="math-container">$r:M\times S^2 \to M\cup (\infty\cdot S^2)$</span> is the ray-tracing function that returns the first intersection point <span class="math-container">$r(x,\omega)$</span> along the ray with origin <span class="math-container">$x$</span> and direction <span class="math-container">$\omega$</span>, 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 <span class="math-container">$L_i(x, \omega)$</span>, however, in vacuum we have <span class="math-container">$L_i(x,\omega) = (GL)(x,\omega) = L(r(x,\omega),-\omega)$</span> - the radiance that arrives at <span class="math-container">$x$</span> from direction <span class="math-container">$\omega$</span> is equal to the radiance leaving <span class="math-container">$r(x,\omega)$</span> in direction <span class="math-container">$-\omega$</span>, which is also fairly intuitive. The term <span class="math-container">$|\cos\theta_x| = |n(x)\cdot \omega|$</span> is due to Lambert's cosine law, but it can be derived also by just integrating a differential <span class="math-container">$2$</span>-density on the sphere (think surface flux integrals from calculus 3). The <span class="math-container">$d\sigma(\omega)$</span> is the area measure on the sphere, i.e. if you take the standard spherical coordinates parametrization <span class="math-container">$x = \sin\theta\cos\phi$</span>, <span class="math-container">$y=\sin\theta\sin\phi$</span>, <span class="math-container">$z=\cos\theta$</span> and your rotate them so that <span class="math-container">$z$</span> is aligned with <span class="math-container">$n(x)$</span>, then you may rewrite <span class="math-container">$d\sigma(\omega) = |\sin\theta| d\theta d\phi$</span> and thus:
<span class="math-container">$$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.$$</span></p>
<p><strong>Formal Solution</strong></p>
<p>For brevity one typically writes the above in operator form <span class="math-container">$L = L_e +TL$</span> where <span class="math-container">$T:(M\times S^2\to \mathbb{R})\to (M\times S^2\to \mathbb{R})$</span> is the operator that stands for the integral. If we define the scattering operator <span class="math-container">$K$</span> and the propagation operator <span class="math-container">$G$</span>:
<span class="math-container">$$(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),$$</span>
then <span class="math-container">$T=K\circ G$</span>.
Since we have the identity <span class="math-container">$L=L_e+TL$</span> we can just plug it in itself once yielding <span class="math-container">$L=L_e + T(L_e + TL) = L_e+TL_e+T^2L$</span>. Repeating this process until infinity yields the series <span class="math-container">$L = \sum_{k=0}^{\infty}T^kL_e$</span>. Formally <span class="math-container">$\sum_{k=0}^{\infty}T^k$</span> is the solution operator for the rendering equation:
<span class="math-container">\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}</span>
The latter is the Neumann series for <span class="math-container">$T$</span>. While <span class="math-container">$\|T\|&lt;1$</span> is a sufficient condition for convergence (i.e. if the brdf is energy conserving), it is not a necessary condition - for specific choices of <span class="math-container">$L_e$</span> the above series may be convergent even if <span class="math-container">$(I-T)^{-1}$</span> is unbounded. Note that <span class="math-container">$f,L,L_e$</span> are actual functions (i.e. <span class="math-container">$T$</span> maps functions to functions). In practice one may want to model specular reflection, point and directional lights etc., which according to &quot;many treatments&quot; 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.</p>
<p><strong>Veach's Extension to Distributions</strong></p>
<p>If you go to 5.A.1, Veach has a short note that <span class="math-container">$\int_{\mathbb{R}} f(x)\delta(x-x_0)\,dx$</span> is just symbolic notation for the evaluation functional <span class="math-container">$\Lambda_{x_0}(f) = f(x_0)$</span>. This symbolic expression could also be framed as <span class="math-container">$\int_{\mathbb{R}} f(x)\delta(x-x_0)dx := \int_{\mathbb{R}} f(x)\,d\delta_{x_0}(x)$</span> where <span class="math-container">$\delta_{x_0}$</span> is the Dirac measure. Note however that <span class="math-container">$\delta$</span> is not the Radon-Nikodym derivative of this measure, i.e. the expression <span class="math-container">$\int_{\mathbb{R}} f(x)\delta(x-x_0)dx$</span> is still purely symbolic if not abuse of notation. When you write <span class="math-container">$$\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),$$</span>
the second expression is purely symbolic and means
<span class="math-container">\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)&amp;:= \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)) &amp;= \rho(x,\omega_r,\omega_o) L_i(x,\omega_r).
\end{align}</span>
The linear operator <span class="math-container">$\Lambda_{\omega_r}$</span> is not the BSDF itself, it is the operator <span class="math-container">$T$</span> when restricted to this specific specular vertex <span class="math-container">$x$</span>, 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 <span class="math-container">$p$</span> is that the definition of <span class="math-container">$K$</span> at point <span class="math-container">$p$</span> is modified to <span class="math-container">$$(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).$$</span>
Formally you could pretend that:
<span class="math-container">$$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).$$</span>
The only issue is that the function/Radon-Nikodym derivative doesn't exist because the Dirac measure <span class="math-container">$\delta_{\omega_r}$</span> is not absolutely continuous, but the notation is convenient enough that everyone abuses it.</p>
<p><strong>Extension to Dirac Delta Light Sources</strong></p>
<p>As we saw above the whole &quot;distribution&quot; 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 &quot;differential bsdf-measure&quot; <span class="math-container">$fd\delta$</span> while at a normal vertex it would be <span class="math-container">$fd\sigma_{x}^{\perp}$</span>. The other way is of course to abuse notation and pretend that <span class="math-container">$\delta$</span> 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.</p>
<p>As an example consider just the <span class="math-container">$0$</span>-th term <span class="math-container">$L_e$</span> of the Neumann expansion with only a directional light source. To get the contribution from zero bounces away we integrate <span class="math-container">$L_e$</span> directly on the film. For a specific pixel  on the film we may assume that the contribution is computed as <span class="math-container">$$\int_{A}\left(\int_{S^2}W_e(x,\omega)L_e(r(x,\omega),-\omega)\,d\sigma_{x}^{\perp}(\omega)\right)\,dx.$$</span>
Here <span class="math-container">$A$</span> is the film area and <span class="math-container">$W_e$</span> 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 <span class="math-container">$L_e$</span> as an actual Dirac delta <span class="math-container">$L_e(y,\omega) = e(y,\omega)\delta(\omega-\omega_d)$</span> in which case you modify the integral to:
<span class="math-container">$$\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.$$</span>
You could however also require that <span class="math-container">$W_e$</span> 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.</p>
<p>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 <span class="math-container">$\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)$</span>, 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.</p>
<p><strong>Conclusion</strong></p>
<p>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. &quot;directional light&quot; 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.</p>
<p><span class="math-container">$$$$</span>
<strong>Appendix: Some Details on the Dirac Delta</strong></p>
<p>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 <span class="math-container">$\delta : \mathbb{R} \to \mathbb{R}\cup \{\pm\infty\}$</span> such that for any function <span class="math-container">$g$</span> we have <span class="math-container">$\int_{\mathbb{R}} g(t)\delta(t-x)\,dt = g(x)$</span>. 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 <span class="math-container">$g$</span> for which they hold, while the least general ones have the most constraints on <span class="math-container">$g$</span>.</p>
<p><strong>Evaluation Functional <span class="math-container">$\delta_x$</span> (Constraint: Functions)</strong></p>
<p>Let <span class="math-container">$V$</span> be the set of functions from some set <span class="math-container">$X$</span> to some set <span class="math-container">$Y$</span>, i.e. <span class="math-container">$g:X\to Y$</span>. Define <span class="math-container">$\delta_x(g) = g(x)$</span> - in this specific case we didn't assume anything about <span class="math-container">$g$</span> except that it is a function. If <span class="math-container">$Y$</span> is a vector space over a field <span class="math-container">$\mathbb{F}$</span>, then <span class="math-container">$V$</span> is also a vector space over <span class="math-container">$\mathbb{F}$</span> with the function-function addition <span class="math-container">$(g+h)(x) := g(x)+h(x)$</span> and scalar-function multiplication <span class="math-container">$(\alpha\cdot g)(x) := \alpha\cdot g(x)$</span>. Specifically if <span class="math-container">$Y=\mathbb{F}$</span> then <span class="math-container">$\delta_x : V \to \mathbb{F}$</span> is from the algebraic dual <span class="math-container">$V'$</span> of <span class="math-container">$V$</span>, that is, the space of linear functions <span class="math-container">$f^*:V\to \mathbb{F}$</span> (also called functionals). Then <span class="math-container">$\delta_x (g) = g(x)$</span> is known as the evaluation functional (here we just assumed that <span class="math-container">$g$</span> 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 <span class="math-container">$\delta_x(g) = g(x)$</span> - the integral from our desired identity has been replaced by some map <span class="math-container">$\delta_x$</span>.</p>
<p>Note also that <span class="math-container">$V$</span> being a set of functions is crucial. If I had taken <span class="math-container">$V$</span> to be <span class="math-container">$L_2(X)$</span> then I would be working with equivalence classes of functions and even the evaluation functional and the meaning of <span class="math-container">$g(x)$</span> 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 <span class="math-container">$L_2$</span>. There is, however, a way to generalize this to some other spaces of equivalence classes, e.g. you could do so for Sobolev spaces <span class="math-container">$H^{q}(\mathcal{X})$</span> where <span class="math-container">$q\in(0,\infty)$</span> and <span class="math-container">$2q&gt; \dim X$</span> (the dimension here is in the sense of manifold dimension).</p>
<p><strong>Measure <span class="math-container">$\delta_x$</span> (Constraint: Measurable Functions)</strong></p>
<p>If you have two measurable spaces <span class="math-container">$(X,\mathcal{X})$</span> and <span class="math-container">$(Y,\mathcal{Y})$</span> (those can be equal) and a function <span class="math-container">$g:X\to Y$</span> such that the preimage <span class="math-container">$g^{-1}$</span> of any set <span class="math-container">$B$</span> from <span class="math-container">$\mathcal{Y}$</span> is from <span class="math-container">$\mathcal{X}$</span>: <span class="math-container">$B\in \mathcal{Y} \implies g^{-1}(B)\in \mathcal{X}$</span>, then <span class="math-container">$g$</span> is termed a measurable function. If you additionally have a measure <span class="math-container">$\mu$</span> (e.g. the Lebesgue measure) then you can further limit <span class="math-container">$\mathcal{X}$</span> to be sets of Lebesgue-measurable sets, and we speak of <span class="math-container">$\mu$</span>-measurable functions. For any <span class="math-container">$A\in \mathcal{X}$</span> you could now define the Dirac measure:
<span class="math-container">$$\delta_x(A) = 1_A(x) = \begin{cases} 1, &amp; x\in A, \\ 0, &amp; x\not\in A.\end{cases}$$</span>
Then you can define the integral of <span class="math-container">$g$</span> w.r.t. <span class="math-container">$\delta_x$</span> as:
<span class="math-container">$$\int_{A} g(t)\,\delta_x(dt) = g(x)\delta_x(A)=
\begin{cases}
g(x), &amp; x\in A, \\
0, &amp;x\not\in A,
\end{cases}$$</span>
Note that <span class="math-container">$0$</span> should be well-defined in <span class="math-container">$Y$</span> for this to make sense (all the practical examples you care about satisfy this however).</p>
<p>Now take <span class="math-container">$X$</span> to be some subset of <span class="math-container">$\mathbb{R}^n$</span> and <span class="math-container">$Y\subseteq \mathbb{R}$</span> then let <span class="math-container">$\mathcal{X}$</span> be the Lebesgue measurable subsets of <span class="math-container">$X$</span>, and let <span class="math-container">$\mu:\mathcal{X}\to\mathbb{R}_+$</span> be the Lebesgue measure,  then you can define integrals of the form <span class="math-container">$\int_A g(t) \,\mu(dt)$</span> for <span class="math-container">$A\in\mathcal{X}$</span>. Each <span class="math-container">$\mu$</span>-measurable function <span class="math-container">$f$</span> can then be interpreted as a density and it induces a measure <span class="math-container">$\mu_f(A) = \int_A f(t)\,\mu(dt)$</span>. The Radon-Nikodym theorem goes in the opposite direction - it says that if you have two measures <span class="math-container">$\mu$</span> and <span class="math-container">$\nu$</span>, then under some constraints (the measures being sigma-finite and <span class="math-container">$\nu$</span> being absolutely continuous w.r.t. <span class="math-container">$\mu$</span>) there exists a <span class="math-container">$\mu$</span>-measurable density function <span class="math-container">$f$</span> such that <span class="math-container">$\nu(A) = \int_A f(t)\,\mu(dt)$</span>. Then <span class="math-container">$f$</span> is termed the Radon-Nikodym derivative of <span class="math-container">$\nu$</span> w.r.t. <span class="math-container">$\mu$</span>: <span class="math-container">$f(t) = \frac{\nu(dt)}{\mu(dt)}$</span> and you can compute integrals w.r.t. <span class="math-container">$\nu$</span> by rewriting them as integrals w.r.t. <span class="math-container">$\mu$</span> (that is e.g. how you define <span class="math-container">$n$</span>-volume integration on <span class="math-container">$n$</span>-manifolds): <span class="math-container">$$\int_A g(t)\nu(dt) = \int_A g(t)f(t) \,\mu(dt).$$</span></p>
<p>In any case, using the above we see that if we have many compatible measures we can pick one common measure <span class="math-container">$\mu$</span> (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
<span class="math-container">$$\int_{A} g(t)\,\delta_x(dt) = \int_A g(t)\delta_{\mu}(t-x)\,d\mu(t) = g(x)\delta_{x}(A).$$</span>
And that would have been possible if it weren't for the fact that <span class="math-container">$\delta_x$</span> is not absolutely continuous w.r.t. <span class="math-container">$\mu$</span>. Notably <span class="math-container">$\mu(\{x\}) = 0$</span> but <span class="math-container">$\delta_x(\{x\})=1$</span> which violates absolute continuity. This means that such a measurable density function <span class="math-container">$\delta_{\mu}(t-x)$</span> does not exist/the Radon-Nikodym derivative <span class="math-container">$\frac{d\delta_x}{d\mu}$</span> 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:
<span class="math-container">$$\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).$$</span>
The first two expressions are purely symbolic and equal to what's on the right by definition. Here <span class="math-container">$\mu$</span> 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:
<span class="math-container">$$\delta_{\mu}(x-t) := \frac{\delta_x(dt)}{\mu(dt)}, \text{ but } \frac{\delta_x(dt)}{\mu(dt)} \text{ doesn't exist.}$$</span>
So the <span class="math-container">$\delta_{\mu}$</span> 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:
<span class="math-container">$$\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) = ?,$$</span>
simply because we never defined what <span class="math-container">$\int_{A} \delta_{\mu}(t-y)\,\delta_x(dt)$</span> means - remember that <span class="math-container">$\delta_{\mu}(t-x)$</span> doesn't actually exist and we only introduced notation for it w.r.t. an integral w.r.t. <span class="math-container">$\mu$</span> <span class="math-container">$\mu$</span>. We could try and define a second &quot;integral&quot; 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):
<span class="math-container">\begin{align}
I &amp;=\int_B \int_{A} g(t)\delta_{\mu}(t-x)\delta_{\mu}(t-y)\,d\mu(dt)\,d\mu(dx) \\
&amp;= \int_B\int_{A} g(t)\delta_{\mu}(t-y)\,\delta_x(dt)\,d\mu(dx) \\
&amp;= \int_B g(x)\delta_{\mu}(x-y)1_A(x)\,\mu(dx) \\
&amp;= \int_B g(x)1_A(x)\,\delta_y(dx) \\
&amp;= g(y)1_A(y)1_B(y) = g(y)1_{A\cap B}(y).\end{align}</span>
So we could claim that the result is a differential measure element <span class="math-container">$g(x)1_A(x)\,\delta_y(dx)$</span> that expects to be integrated. Later on we can argue what this should mean if we treat <span class="math-container">$dx$</span> as a finite measure set and we shrink it down to zero - then naturally you would expect this result to blow up to infinity.</p>
<p>However, whenever you have a specular light + specular bsdf, another way to model it is as follows:
<span class="math-container">\begin{align}
I &amp;=\int_{A} g(t)\delta_{\mu}(t-x)\delta_{\delta_x}(t-y)\,d\mu(dt) \\
&amp;= \int_A g(t)\delta_{\delta_x}(t-y)\,d\delta_x(dt) \\
&amp;= \int_A g(t)\delta_{\delta_x}(t-y)\frac{\delta_y(dt)}{\delta_x(dt)}\,d\delta_x(dt) \\
&amp;= \int_A g(t) \delta_y(dt) = g(y)1_A(y).
\end{align}</span>
Now this corresponds to what is done in graphics <span class="math-container">$\delta_{\delta_x}$</span>, 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.</p>
<p>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 <span class="math-container">$K$</span> at specific points of the scene - so in theory you are free to pick whichever you like.</p>
<p><strong>Distribution <span class="math-container">$\delta_x$</span> (Constraint: <span class="math-container">$C^{k}$</span> Functions)</strong></p>
<p>You may require <span class="math-container">$\delta_{x}$</span> to be from the continuous dual <span class="math-container">$V^*$</span> (you need to have a topology to define the continuous dual unlike for the algebraic dual). If <span class="math-container">$V=C^{\infty}_c(X)$</span> then the continuos dual <span class="math-container">$V^*$</span> is the space of distributions. A regular distribution is one arising as an integral of a locally <span class="math-container">$L_1$</span> integrable function, e.g. <span class="math-container">$f^*(\phi) = \int f(t)\phi(t)\,dt$</span>. You can multiply a distribution <span class="math-container">$D$</span> by another distribution <span class="math-container">$f^*$</span>, if the latter is a regular distribution arising from <span class="math-container">$f\in C^{\infty}(X)$</span>. 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 <span class="math-container">$\delta$</span> such that <span class="math-container">$\delta_x(\phi) = \int \delta(t-x) \phi(t)\,dt$</span>. So once again such notation would be purely symbolic. Furthermore while distribution multiplication is defined between a regular distribution (arising from a function in <span class="math-container">$C^k$</span>) 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.</p>
<p><strong>Reproducing Kernel <span class="math-container">$\delta_x$</span> (Constraint: Reproducing Kernel Hilbert Spaces)</strong></p>
<p>I should note that the Dirac delta exists as a function in the sense <span class="math-container">$\int \delta(t-x) f(t)\,dt = f(x)$</span> 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).</p>
<p><strong>The Dirac Delta as a Limit</strong></p>
<p>You can study the following limit:
<span class="math-container">$$\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).$$</span>
As you can see it acts as a Dirac delta (for reasonable <span class="math-container">$f$</span>). One may be tempted to move the limit inside the integral and state:
<span class="math-container">$$\delta(t-x) = \lim_{\sigma\to 0}\exp\left(-\frac{(t-x)^2}{2\sigma^2}\right),$$</span>
but you actually cannot do so as the two are not equivalent:
<span class="math-container">$$\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.$$</span>
However as before you could formally define the following symbolic expression:
<span class="math-container">$$\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.$$</span>
Then you could choose to define the product of two Dirac deltas as:
<span class="math-container">$$\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.$$</span>
My guess is that this integral should give you the following (I haven't checked it):
<span class="math-container">$$\lim_{\sigma\to 0}\exp\left(-\frac{(t-x)^2}{2\sigma^2}\right) = \begin{cases} \infty, &amp; x=y, \\ 0, &amp;x\ne y\end{cases}.$$</span>
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 <span class="math-container">$\delta_x$</span> as a measure is the physically more relevant one, since it would agree with this limiting process.</p>

CC BY-SA 4.0
2023-11-20T23:31:24.117Z
2023-11-23T16:54:09.530Z
2023-11-23T16:54:09.530Z
lightxbulb
<p>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.</p>

<p>This causes problems in rendering equation.</p>

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

<p>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.</p>

<p>I know these are subtle problems and in practice can be solved with few <code>if</code>s 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 <code>if</code>s in the right place.</p>

<hr>

<p>How to deal with this problem? </p>

<p>Not sure entirely but as I can remember Erich Veach <a href="https://graphics.stanford.edu/papers/veach_thesis/" rel="noreferrer">in his thesis</a> address this problem only slightly and tries to get away with it by saying that BSDF and radiance are <a href="https://en.wikipedia.org/wiki/Distribution_(mathematics)" rel="noreferrer">distributions</a>. 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.</p>

<p>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?</p>

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

<p>(I'm only aware of <a href="http://cybertron.cg.tu-berlin.de/lessig/projects/dissertation/index.html" rel="noreferrer">work</a> of Christian Lessig, but his work is still beyond my mathematical reach.)</p>

<hr>

<p>I already have <a href="http://tomsmathjourneys.blogspot.cz/2015/09/better-code-with-radiance-and-bsdf-as.html" rel="noreferrer">proposal</a> 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.</p>

<p>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.</p>

CC BY-SA 4.0
2023-11-23T16:57:56.890Z
@tom I updated my answer to be fairly comprehensive. I didn't go into details on the Dirac deltas as the non-existent function inducing the Dirac distribution since really that's just built on top of measure theory, i.e. you don't really gain any new insights. The limit case is probably the most insightful, which lets you imagine the limiting process of actual light sources, bsdf, sensors as these approach a Dirac delta, and it suggests that the abuse of notation probably just works. And if you get a Dirac at the end you can just visualize it as infinity (unless you integrate again).
lightxbulb
CC BY-SA 3.0
2015-09-10T10:08:48.700Z
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 `if`s), 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.
Wumpf
CC BY-SA 3.0
2015-09-10T10:25:34.243Z
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
CC BY-SA 3.0
2015-09-10T10:31:31.473Z
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
CC BY-SA 3.0
2015-09-10T18:28:11.070Z
@tom The rigorous mathematics that underlies delta distributions is [measure theory](https://en.wikipedia.org/wiki/Measure_%28mathematics%29). See the [definition of the Dirac delta as a measure](https://en.wikipedia.org/wiki/Dirac_delta_function#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.
Nathan Reed
CC BY-SA 3.0
2015-09-11T01:48:58.433Z
I don't understand why the BSDF can have a dirac. BSDF can represents 100% reflectivity using a value of 1, not infinite ? you don't need to make so that the integral is 1 over the hemisphere, you need to make so that it is less than one, strictly. or there is something i don't get. more than 1, for some angle, would mean light from other directions than the perfect reflection also contributes to this particular angle. which would not be a mirror, but some kind of lens !
v.oddou
CC BY-SA 3.0
2015-09-11T14:22:31.350Z
@Rahul Good idea to think about BSDF as linear operator taking $L_i$ to $L_o$, but I still wonder if it is somehow advantageous to define radiance as distribution, because it seams to me that measure is sufficient, in which case you can too think about BSDF as linear mapping from one measure $L_i$ to another $L_o$ and in addition it can be represented it as integral of some measure valued function over measure $L_i$. And decomposition of that measure valued function into abs. continuous and singular part with the respect to the solid angle gives you diffusive and specular part of BSDF.
tom
<raytracing><radiosity><theory>
Problem with definition of BSDF and radiance
tom
# Problem with definition of BSDF and radiance

(Asked by: username_0 on 2015-09-10)

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 `if`s 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 `if`s in the right place.

---

How to deal with this problem?

Not sure entirely but as I can remember Erich Veach [in his thesis](https://graphics.stanford.edu/papers/veach_thesis/) address this problem only slightly and tries to get away with it by saying that BSDF and radiance are [distributions](https://en.wikipedia.org/wiki/Distribution_(mathematics)). 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](http://cybertron.cg.tu-berlin.de/lessig/projects/dissertation/index.html) of Christian Lessig, but his work is still beyond my mathematical reach.)

---

I already have [proposal](http://tomsmathjourneys.blogspot.cz/2015/09/better-code-with-radiance-and-bsdf-as.html) 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.

**Tags:** raytracing, radiosity, theory

> **Comments:**
> **username_2**: 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 `if`s), 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.  
> **username_0**: 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.  
> **username_0**: 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.  
> **username_3**: @username_0 The rigorous mathematics that underlies delta distributions is [measure theory](https://en.wikipedia.org/wiki/Measure_%28mathematics%29). See the [definition of the Dirac delta as a measure](https://en.wikipedia.org/wiki/Dirac_delta_function#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.  
> **username_4**: I don't understand why the BSDF can have a dirac. BSDF can represents 100% reflectivity using a value of 1, not infinite ? you don't need to make so that the integral is 1 over the hemisphere, you need to make so that it is less than one, strictly. or there is something i don't get. more than 1, for some angle, would mean light from other directions than the perfect reflection also contributes to this particular angle. which would not be a mirror, but some kind of lens !  
> **username_0**: @Rahul Good idea to think about BSDF as linear operator taking $L_i$ to $L_o$, but I still wonder if it is somehow advantageous to define radiance as distribution, because it seams to me that measure is sufficient, in which case you can too think about BSDF as linear mapping from one measure $L_i$ to another $L_o$ and in addition it can be represented it as integral of some measure valued function over measure $L_i$. And decomposition of that measure valued function into abs. continuous and singular part with the respect to the solid angle gives you diffusive and specular part of BSDF.  
> **username_1**: @username_0 I updated my answer to be fairly comprehensive. I didn't go into details on the Dirac deltas as the non-existent function inducing the Dirac distribution since really that's just built on top of measure theory, i.e. you don't really gain any new insights. The limit case is probably the most insightful, which lets you imagine the limiting process of actual light sources, bsdf, sensors as these approach a Dirac delta, and it suggests that the abuse of notation probably just works. And if you get a Dirac at the end you can just visualize it as infinity (unless you integrate again).  

---

## 2 answers

### Answer by username_5 on 2023-11-19. Score: 4

## 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](https://pbr-book.org/3ed-2018/Reflection_Models/Specular_Reflection_and_Transmission#SpecularReflection)):
$$
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](https://en.wikipedia.org/wiki/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](https://www.wolframalpha.com/input?i=integrate+1%2Fsqrt%28x%29+from+0+to+1) or [read a very short proof of](https://en.wikipedia.org/wiki/Improper_integral#Examples).)

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](https://en.wikipedia.org/wiki/Lebesgue_integration#Basic_theorems_of_the_Lebesgue_integral) (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](https://en.wikipedia.org/wiki/Dirac_measure). 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](https://en.wikipedia.org/wiki/Lebesgue_integration#Via_improper_Riemann_integral) 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<t<L(x,\omega\_r)$ 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.

> **Comments:**
> **username_1**: 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 @username_0 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.  
> **username_1**: 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.  
> **username_5**: @username_1 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.  
> **username_1**: 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.  
> **username_0**: Nice answer but I have to agree with @username_1 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.  
> **username_0**: Oh this is already described in the blog post I linked ... well it has been 8 years :)  
> **username_1**: @username_0 I actually don't get user106's comment. He says that the bsdf is a linear map, but the bsdf is not a linear map eating $L_i$ and spitting out $L_o$. In vacuum $L_i(x, \omega) = L_o(r(x,\omega),-\omega)$. The rendering equation is an equation specifically because $L\equiv L_o$ shows up on both sides $L = L_e +TL$. The formal solution is $L = (I-T)^{-1}L_e = \sum_{k=0}^{\infty}T^k L_e$. If we pick two specular/Dirac vertices at the end it isn't really an issue either imo, you just have to define the distribution mult differently, e.g. if you use a joint measure it's fine.  

---

### Answer by username_1 on 2023-11-20. Score: 3

**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.


---

