T}}\mathrm{cov}(y_{\perp a,c})^{1/2}&Q_{2}\end{pmatrix}.$

are positive definite. By Lemma S6, there exists $u$ such that
$Q_{1}=\mathrm{cov}(u_{\perp c})$, $R_{1}=R_{a\sim u\mid c}$ and
$R_{2}=R_{y\sim u\mid a,c}$. This shows that $\\{R_{y\sim u\mid a,c},R_{a\sim
u\mid c},\mathrm{cov}(u_{\perp c})\\}$ can take arbitrary values in
$\mathbb{B}_{\mathrm{dim}(y)\times\mathrm{dim}(u)}\times\mathbb{B}_{\mathrm{dim}(a)\times\mathrm{dim}(u)}\times\mathbb{S}_{\mathrm{dim}(u)}^{++}$.

∎

###### Proof of Theorem 3.

(i) When $u$ is a scalar, Theorem 1 and (20) imply
$\theta_{1}=\tilde{\theta}_{1}+T_{1}\phi_{1}+T_{2}^{\mathrm{\scriptscriptstyle
T}}\phi_{2}$, where

$\displaystyle\phi_{1}$ $\displaystyle=\frac{R_{y\sim u\mid a,m,c}R_{a\sim
u\mid c}}{(1-R_{a\sim u\mid c}^{2})^{1/2}(1-\|R_{m\sim u\mid
a,c}\|_{2}^{2})^{1/2}},$ (S15) $\displaystyle\phi_{2}$
$\displaystyle=\frac{R_{y\sim u\mid a,m,c}}{(1-\|R_{m\sim u\mid
a,c}\|_{2}^{2})^{1/2}}\cdot R_{m\sim u\mid a,c},$ (S16)

only depend on sensitivity parameters $\\{R_{y\sim u\mid a,m,c},R_{m\sim u\mid
a,c},R_{a\sim u\mid c}\\}$.

When $u$ is a vector, by Theorem 2, we have

$\theta_{1}=\tilde{\theta}_{1}-\mathrm{var}(y_{\perp a,m,c})^{1/2}R_{y\sim
u\mid a,m,c}\mathrm{cov}(u_{\perp a,m,c})^{-1/2}\mathrm{cov}(u_{\perp
c,m})^{1/2}R_{a\sim u\mid c,m}^{\mathrm{\scriptscriptstyle
T}}\mathrm{var}(a_{\perp c,m})^{-1/2}.$ (S17)

By Lemma S3, we can represent the unnatural $R$ parameter $R_{a\sim u\mid
c,m}$ as

$R_{a\sim u\mid c,m}=\frac{(R_{a\sim u\mid c}-R_{m\sim a\mid
c}^{\mathrm{\scriptscriptstyle T}}R_{m\sim u\mid c})\mathrm{cov}(u_{\perp
c})^{1/2}\mathrm{cov}(u_{\perp c,m})^{-1/2}}{(1-\|R_{m\sim a\mid
c}\|_{2}^{2})^{1/2}}$ (S18)

and we can further represent the $R$ parameter $R_{m\sim u\mid c}$ as

$R_{m\sim u\mid c}=\mathrm{cov}(m_{\perp c})^{-1/2}\mathrm{cov}(m_{\perp
a,c})^{1/2}R_{m\sim u\mid a,c}\mathrm{cov}(u_{\perp
a,c})^{1/2}\mathrm{cov}(u_{\perp c})^{-1/2}+R_{m\sim a\mid c}R_{a\sim u\mid
c}.$ (S19)

Then, (S17)–(S19) imply
$\theta_{1}=\tilde{\theta}_{1}+T_{1}\phi_{1}+T_{2}^{\mathrm{\scriptscriptstyle
T}}\phi_{2}$, where

$\displaystyle\phi_{1}$ $\displaystyle=R_{y\sim u\mid
a,m,c}\mathrm{cov}(u_{\perp a,m,c})^{-1/2}\mathrm{cov}(u_{\perp
c})^{1/2}R_{a\sim u\mid c}^{\mathrm{\scriptscriptstyle T}},$ (S20)
$\displaystyle\phi_{2}$ $\displaystyle=R_{m\sim u\mid
a,c}\mathrm{cov}(u_{\perp a,c})^{1/2}\mathrm{cov}(u_{\perp
a,m,c})^{-1/2}R_{y\sim u\mid a,m,c}^{\mathrm{\scriptscriptstyle T}},$ (S21)

only depend on sensitivity parameters $\\{R_{y\sim u\mid a,m,c},R_{m\sim u\mid
a,c},R_{a\sim u\mid c},\mathrm{cov}(u_{\perp c})\\}$ after we further express
$\mathrm{cov}(u_{\perp a,c})$ and $\mathrm{cov}(u_{\perp a,m,c})$ in

$\displaystyle\mathrm{cov}(u_{\perp a,c})$
$\displaystyle=\mathrm{cov}(u_{\perp c})-\mathrm{cov}(u_{\perp
c})^{1/2}R_{a\sim u\mid c}^{\mathrm{\scriptscriptstyle T}}R_{a\sim u\mid
c}\mathrm{cov}(u_{\perp c})^{1/2},$ (S22) $\displaystyle\mathrm{cov}(u_{\perp
a,m,c})$ $\displaystyle=\mathrm{cov}(u_{\perp a,c})-\mathrm{cov}(u_{\perp
a,c})^{1/2}R_{m\sim u\mid a,c}^{\mathrm{\scriptscriptstyle T}}R_{m\sim u\mid
a,c}\mathrm{cov}(u_{\perp a,c})^{1/2}.$ (S23)

(ii) We first discuss the case when $u$ is a scalar. Our proof consists of two
steps.

First, we show that for any sensitivity parameters $R_{y\sim u\mid a,m,c}$,
$R_{m\sim u\mid a,c}$, $R_{a\sim u\mid c}$ under the constraint that $R_{y\sim
u\mid a,m,c}^{2}\leq\rho_{y}$, $\|R_{m\sim u\mid a,c}\|_{2}^{2}\leq\rho_{m}$,
$R_{a\sim u\mid c}^{2}\leq\rho_{a}$, there exist $t_{1}\in[-1,1]$,
$t_{2}\in\partial\mathbb{B}_{\mathrm{dim}(m)}$,
$t_{3}\in[0,\arctan(\rho_{m}^{1/2}/(1-\rho_{m})^{1/2})]$ satisfying
$\phi_{1}=\\{\rho_{y}\rho_{a}/(1-\rho_{a})\\}^{1/2}(t_{1}/\cos t_{3})$ and
$\phi_{2}=\rho_{y}^{1/2}\tan t_{3}\cdot t_{2}$. We discuss the following
cases:

  * •

If $R_{y\sim u\mid a,m,c}=0$, then $\phi_{1}=\phi_{2}=0$ and we can choose
$t_{1}=0$, $t_{3}=0$, and arbitrary
$t_{2}\in\partial\mathbb{B}_{\mathrm{dim}(m)}$.

  * •

If $R_{a\sim u\mid c}=0$ and $R_{m\sim u\mid a,c}=0$, then
$\phi_{1}=\phi_{2}=0$ and we can choose $t_{1}=0$, $t_{3}=0$, and arbitrary
$t_{2}\in\partial\mathbb{B}_{\mathrm{dim}(m)}$.

  * •

If $R_{y\sim u\mid a,m,c}\not=0$, $R_{a\sim u\mid c}=0$ and $R_{m\sim u\mid
a,c}\not=0$, then $\rho_{a}>0$ and $\rho_{m}>0$, and we can choose $t_{1}=0$,
$t_{2}=\mathrm{sgn}(R_{y\sim u\mid a,m,c})\cdot R_{m\sim u\mid a,c}/\|R_{m\sim
u\mid a,c}\|_{2}$ and

$t_{3}=\arctan\left(\frac{|R_{y\sim u\mid
a,m,c}|}{\rho_{y}^{1/2}}\cdot\frac{\|R_{m\sim u\mid a,c}\|_{2}}{(1-\|R_{m\sim
u\mid a,c}\|_{2}^{2})^{1/2}}\right).$

  * •

If $R_{y\sim u\mid a,m,c}\not=0$, $R_{a\sim u\mid c}\not=0$ and $R_{m\sim
u\mid a,c}=0$, then $\rho_{y}>0$ and $\rho_{a}>0$, and we can choose
$t_{3}=0$, arbitrary $t_{2}\in\partial\mathbb{B}_{\mathrm{dim}(m)}$ and

$t_{1}=\left(\frac{1-\rho_{a}}{\rho_{y}\rho_{a}}\right)^{1/2}\frac{R_{y\sim
u\mid a,m,c}R_{a\sim u\mid c}}{(1-R_{a\sim u\mid c}^{2})^{1/2}}.$

  * •

If $R_{y\sim u\mid a,m,c}\not=0$, $R_{m\sim u\mid a,c}\not=0$, $R_{a\sim u\mid
c}\not=0$, then $\rho_{y}>0$, $\rho_{m}>0$, $\rho_{a}>0$, and we can choose

$\begin{split}t_{2}&=\mathrm{sgn}(R_{y\sim u\mid a,m,c})\cdot\frac{R_{m\sim
u\mid a,c}}{\|R_{m\sim u\mid a,c}\|_{2}},\\\
t_{3}&=\arctan\left(\frac{|R_{y\sim u\mid
a,m,c}|}{\rho_{y}^{1/2}}\cdot\frac{\|R_{m\sim u\mid a,c}\|_{2}}{(1-\|R_{m\sim
u\mid a,c}\|_{2}^{2})^{1/2}}\right),\\\
t_{1}&=\left(\frac{1-\rho_{a}}{\rho_{y}\rho_{a}}\right)^{1/2}\frac{R_{y\sim
u\mid a,m,c}R_{a\sim u\mid c}}{(1-R_{a\sim u\mid c}^{2})^{1/2}(1-\|R_{m\sim
u\mid a,c}\|_{2}^{2})^{1/2}}\cos t_{3}.\end{split}$

Here we only check $t_{1}\in[-1,1]$. By the identity that $\cos(\arctan
t)=(1+t^{2})^{-1/2}$ for $t\in\mathbb{R}$, we have

$\begin{split}t_{1}^{2}&=\frac{1-\rho_{a}}{\rho_{y}\rho_{a}}\cdot\frac{R_{y\sim
u\mid a,m,c}^{2}R_{a\sim u\mid c}^{2}}{(1-R_{a\sim u\mid c}^{2})(1-\|R_{m\sim
u\mid a,c}\|_{2}^{2})}\cdot\left(1+\frac{R_{y\sim u\mid
a,m,c}^{2}}{\rho_{y}}\cdot\frac{\|R_{m\sim u\mid a,c}\|_{2}^{2}}{1-\|R_{m\sim
u\mid a,c}\|_{2}^{2}}\right)^{-1}\\\ &\leq\frac{R_{y\sim u\mid
a,m,c}^{2}}{\rho_{y}(1-\|R_{m\sim u\mid
a,c}\|_{2}^{2})}\cdot\left\\{1+\frac{R_{y\sim u\mid
a,m,c}^{2}}{\rho_{y}(1-\|R_{m\sim u\mid a,c}\|_{2}^{2})}\right\\}^{-1}\\\
&\leq 1.\end{split}$

Second, we show that for any $t_{1}\in[-1,1]$,
$t_{2}\in\partial\mathbb{B}_{\mathrm{dim}(m)}$,
$t_{3}\in[0,\arctan(\rho_{m}^{1/2}/(1-\rho_{m})^{1/2})]$, there exist
sensitivity parameters $R_{y\sim u\mid a,m,c}$, $R_{m\sim u\mid a,c}$,
$R_{a\sim u\mid c}$ satisfying $R_{y\sim u\mid a,m,c}^{2}\leq\rho_{y}$,
$\|R_{m\sim u\mid a,c}\|_{2}^{2}\leq\rho_{m}$, $R_{a\sim u\mid
c}^{2}\leq\rho_{a}$, and
$\phi_{1}=\\{\rho_{y}\rho_{a}/(1-\rho_{a})\\}^{1/2}\cdot(t_{1}/\cos t_{3})$,
$\phi_{2}=\rho_{y}^{1/2}\tan\theta\cdot y$. We can choose parameters $R_{y\sim
u\mid a,m,c}=\rho_{y}^{1/2}$, $R_{m\sim u\mid a,c}=\sin t_{3}\cdot t_{2}$, and

$R_{a\sim u\mid
c}=\left(\frac{\rho_{a}}{t_{1}^{2}\rho_{a}-\rho_{a}+1}\right)^{1/2}t_{1},$

which satisfy the conditions.

We next discuss the case when $u$ is a vector. Since we can take $R_{y\sim
u\mid a,m,c}$, $R_{m\sim u\mid a,c}$, $R_{a\sim u\mid c}$ with only first
nonzero column, the range of $\\{\phi_{1},\phi_{2}\\}$ when $u$ is a scalar is
a subset of the range of $\\{\phi_{1},\phi_{2}\\}$ when $u$ is a vector.

We re-express the coefficient of $u$ from the least squares of $y$ on
$(a,m,c,u)$

$\mathrm{cov}(y_{\perp a,m,c},u_{\perp a,m,c})\mathrm{cov}(u_{\perp
a,m,c})^{-1}=\mathrm{cov}(y_{\perp a,m,c})^{1/2}R_{y\sim u\mid
a,m,c}\mathrm{cov}(u_{\perp a,m,c})^{-1/2},$

and the coefficient of $a$ from the least squares of $u$ on $(a,c)$

$\mathrm{cov}(u_{\perp c},a_{\perp c})\mathrm{var}(a_{\perp
c})^{-1}=\mathrm{cov}(u_{\perp c})^{1/2}R_{a\sim u\mid
c}^{\mathrm{\scriptscriptstyle T}}\mathrm{var}(a_{\perp c})^{-1/2}.$

For any invertible linear transformation $\Gamma u$ of $u$, the coefficient of
$\Gamma u$ from the least squares of $y$ on $(a,m,c,\Gamma u)$ is
$\mathrm{cov}(y_{\perp a,m,c})^{1/2}R_{y\sim u\mid a,m,c}\mathrm{cov}(u_{\perp
a,m,c})^{-1/2}\Gamma^{-1}$, and the coefficient of $a$ from the least squares
of $\Gamma u$ on $(a,c)$ is $\Gamma\mathrm{cov}(u_{\perp c})^{1/2}R_{a\sim
u\mid c}^{\mathrm{\scriptscriptstyle T}}\mathrm{var}(a_{\perp c})^{-1/2}$.
This implies that $\phi_{1}=R_{y\sim u\mid a,m,c}\mathrm{cov}(u_{\perp
a,m,c})^{-1/2}\mathrm{cov}(u_{\perp c})^{1/2}R_{a\sim u\mid
c}^{\mathrm{\scriptscriptstyle T}}$ is invariant to invertible linear
transformation of $u$. Similarly, we can show that $\phi_{2}$ is invariant to
invertible linear transformation of $u$. Moreover, since $\|R_{y\sim u\mid
a,m,c}\|_{2}^{2}$, $\|R_{m\sim u\mid a,c}\|_{2}^{2}$, $\|R_{a\sim u\mid
c}\|_{2}^{2}$ are invariant to invertible linear transformation of $u$, we can
assume $\mathrm{cov}(u_{\perp a,c})=I$ and find the range of

$\begin{split}\tilde{\phi}_{1}&=R_{y\sim u\mid a,m,c}\left(I-R_{m\sim u\mid
a,c}^{\mathrm{\scriptscriptstyle T}}R_{m\sim u\mid
a,c}\right)^{-1/2}\left(I-R_{a\sim u\mid c}^{\mathrm{\scriptscriptstyle
T}}R_{a\sim u\mid c}\right)^{-1/2}R_{a\sim u\mid
c}^{\mathrm{\scriptscriptstyle T}},\\\ \tilde{\phi}_{2}&=R_{y\sim u\mid
a,m,c}\left(I-R_{m\sim u\mid a,c}^{\mathrm{\scriptscriptstyle T}}R_{m\sim
u\mid a,c}\right)^{-1/2}R_{m\sim u\mid a,c}^{\mathrm{\scriptscriptstyle
T}}\end{split}$

under the constraint that $\|R_{y\sim u\mid a,m,c}\|_{2}^{2}\leq\rho_{y}$,
$\|R_{m\sim u\mid a,c}\|_{2}^{2}\leq\rho_{m}$, $\|R_{a\sim u\mid
c}\|_{2}^{2}\leq\rho_{a}$. Then, we have

$\begin{split}|\tilde{\phi}_{1}|&\leq\frac{\|R_{y\sim u\mid
a,m,c}\|_{2}\|R_{a\sim u\mid c}\|_{2}}{(1-\|R_{a\sim u\mid
c}\|_{2}^{2})^{1/2}(1-\|R_{m\sim u\mid a,c}\|_{2}^{2})^{1/2}},\\\
\|\tilde{\phi}_{2}\|_{2}&\leq\frac{\|R_{y\sim u\mid a,m,c}\|_{2}\|R_{m\sim
u\mid a,c}\|_{2}}{(1-\|R_{m\sim u\mid a,c}\|_{2}^{2})^{1/2}}.\end{split}$

We have already shown that the range of
$\\{|\tilde{\phi}_{1}|,\|\tilde{\phi}_{2}\|_{2}\\}$ is

$\left\\{\left(\left(\frac{\rho_{y}\rho_{a}}{1-\rho_{a}}\right)^{1/2}\frac{t_{1}}{\cos
t_{3}},\rho_{y}^{1/2}\tan
t_{3}\right):t_{1}\in[0,1],t_{3}\in\left[0,\arctan\left(\frac{\rho_{m}^{1/2}}{(1-\rho_{m})^{1/2}}\right)\right]\right\\}.$

Thus, the range of $\\{\phi_{1},\phi_{2}\\}$ is a subset of

$\left\\{\left(\left(\frac{\rho_{y}\rho_{a}}{1-\rho_{a}}\right)^{1/2}\frac{t_{1}}{\cos
t_{3}},\rho_{y}^{1/2}\tan t_{3}\cdot
t_{2}\right):t_{1}\in[-1,1],t_{2}\in\partial\mathbb{B}_{\mathrm{dim}(m)},t_{3}\in\left[0,\arctan\left(\frac{\rho_{m}^{1/2}}{(1-\rho_{m})^{1/2}}\right)\right]\right\\}.$

Therefore, the range of $\\{\phi_{1},\phi_{2}\\}$ is

$\left\\{\left(\left(\frac{\rho_{y}\rho_{a}}{1-\rho_{a}}\right)^{1/2}\frac{t_{1}}{\cos
t_{3}},\rho_{y}^{1/2}\tan t_{3}\cdot
t_{2}\right):t_{1}\in[-1,1],t_{2}\in\partial\mathbb{B}_{\mathrm{dim}(m)},t_{3}\in\left[0,\arctan\left(\frac{\rho_{m}^{1/2}}{(1-\rho_{m})^{1/2}}\right)\right]\right\\},$

whenever $u$ is a scalar or a vector.

(iii) It follows from (i) and (ii). ∎

###### Proof of Theorem 4.

(i) When $u$ is a scalar, Theorem 1 implies
$\theta_{3}^{\mathrm{\scriptscriptstyle
T}}\beta_{1}=(\tilde{\theta}_{3}+T_{3}\phi_{3})^{\mathrm{\scriptscriptstyle
T}}(\tilde{\beta}_{1}+T_{4}\phi_{4})$, where

$\displaystyle\phi_{3}$ $\displaystyle=\frac{R_{a\sim u\mid c}}{(1-R_{a\sim
u\mid c}^{2})^{1/2}}\cdot R_{m\sim u\mid a,c},$ (S24) $\displaystyle\phi_{4}$
$\displaystyle=\frac{R_{y\sim u\mid a,m,c}}{(1-\|R_{m\sim u\mid
a,c}\|_{2}^{2})^{1/2}}\cdot R_{m\sim u\mid a,c},$ (S25)

only depend on sensitivity parameters $\\{R_{y\sim u\mid a,m,c},R_{m\sim u\mid
a,c},R_{a\sim u\mid c}\\}$.

When $u$ is a vector, Theorem 2 implies
$\theta_{3}^{\mathrm{\scriptscriptstyle
T}}\beta_{1}=(\tilde{\theta}_{3}+T_{3}\phi_{3})^{\mathrm{\scriptscriptstyle
T}}(\tilde{\beta}_{1}+T_{4}\phi_{4})$, where

$\displaystyle\phi_{3}$ $\displaystyle=R_{m\sim u\mid
a,c}\mathrm{cov}(u_{\perp a,c})^{-1/2}\mathrm{cov}(u_{\perp c})^{1/2}R_{a\sim
u\mid c}^{\mathrm{\scriptscriptstyle T}},$ (S26) $\displaystyle\phi_{4}$
$\displaystyle=R_{m\sim u\mid a,c}\mathrm{cov}(u_{\perp
a,c})^{-1/2}\mathrm{cov}(u_{\perp a,m,c})^{1/2}R_{y\sim u\mid
a,m,c}^{\mathrm{\scriptscriptstyle T}},$ (S27)

only depend on sensitivity parameters $\\{R_{y\sim u\mid a,m,c},R_{m\sim u\mid
a,c},R_{a\sim u\mid c},\mathrm{cov}(u_{\perp c})\\}$ after we further express
$\mathrm{cov}(u_{\perp a,c})$ and $\mathrm{cov}(u_{\perp a,m,c})$ in
(S22)–(S23).

(ii) We first discuss the case when $u$ is a scalar. Our proof consists of two
steps.

First, we show that for any sensitivity parameters $R_{y\sim u\mid a,m,c}$,
$R_{m\sim u\mid a,c}$, $R_{a\sim u\mid c}$ under the constraint that $R_{y\sim
u\mid a,m,c}^{2}\leq\rho_{y}$, $\|R_{m\sim u\mid a,c}\|_{2}^{2}\leq\rho_{m}$,
$R_{a\sim u\mid c}^{2}\leq\rho_{a}$, there exist
$t_{1}\in\partial\mathbb{B}_{\mathrm{dim}(m)}$, $t_{2},t_{3}\in[-1,1]$
satisfying $\phi_{3}=\\{\rho_{m}\rho_{a}/(1-\rho_{a})\\}^{1/2}t_{2}\cdot
t_{1}$ and $\phi_{4}=\\{\rho_{y}\rho_{m}/(1-\rho_{m})\\}^{1/2}t_{3}\cdot
t_{1}$. We discuss the following cases:

  * •

If $R_{m\sim u\mid a,c}=0$, then $\phi_{3}=\phi_{4}=0$ and we can choose
$t_{2}=0$, $t_{3}=0$, and arbitrary
$t_{1}\in\partial\mathbb{B}_{\mathrm{dim}(m)}$.

  * •

If $R_{y\sim u\mid a,m,c}=0$ and $R_{a\sim u\mid c}=0$, then
$\phi_{3}=\phi_{4}=0$ and we can choose $t_{2}=0$, $t_{3}=0$, and arbitrary
$t_{1}\in\partial\mathbb{B}_{\mathrm{dim}(m)}$.

  * •

If $R_{m\sim u\mid a,c}\not=0$, $R_{y\sim u\mid a,m,c}\not=0$ and $R_{a\sim
u\mid c}=0$, then $\rho_{m}>0$ and $\rho_{y}>0$, and we can choose
$t_{1}=R_{m\sim u\mid a,c}/\|R_{m\sim u\mid a,c}\|_{2}$, $t_{2}=0$ and

$t_{3}=\left(\frac{1-\rho_{m}}{\rho_{y}\rho_{m}}\right)^{1/2}\frac{R_{y\sim
u\mid a,m,c}\|R_{m\sim u\mid a,c}\|_{2}}{(1-\|R_{m\sim u\mid
a,c}\|_{2}^{2})^{1/2}}.$

  * •

If $R_{m\sim u\mid a,c}\not=0$, $R_{y\sim u\mid a,m,c}=0$ and $R_{a\sim u\mid
c}\not=0$, then $\rho_{m}>0$ and $\rho_{a}>0$, and we can choose
$t_{1}=R_{m\sim u\mid a,c}/\|R_{m\sim u\mid a,c}\|_{2}$, $t_{3}=0$ and

$t_{2}=\left(\frac{1-\rho_{a}}{\rho_{a}\rho_{m}}\right)^{1/2}\frac{R_{a\sim
u\mid c}\|R_{m\sim u\mid a,c}\|_{2}}{(1-R_{a\sim u\mid c}^{2})^{1/2}}.$

  * •

If $R_{m\sim u\mid a,c}\not=0$, $R_{y\sim u\mid a,m,c}\not=0$ and $R_{a\sim
u\mid c}\not=0$, then $\rho_{m}>0$, $\rho_{y}>0$, $\rho_{a}>0$, and we can
choose

$\begin{split}t_{1}&=\frac{R_{m\sim u\mid a,c}}{\|R_{m\sim u\mid
a,c}\|_{2}},\\\
t_{2}&=\left(\frac{1-\rho_{a}}{\rho_{a}\rho_{m}}\right)^{1/2}\frac{R_{a\sim
u\mid c}\|R_{m\sim u\mid a,c}\|_{2}}{(1-R_{a\sim u\mid c}^{2})^{1/2}},\\\
t_{3}&=\left(\frac{1-\rho_{m}}{\rho_{y}\rho_{m}}\right)^{1/2}\frac{R_{y\sim
u\mid a,m,c}\|R_{m\sim u\mid a,c}\|_{2}}{(1-\|R_{m\sim u\mid
a,c}\|_{2}^{2})^{1/2}}.\end{split}$

Second, we show that for any $t_{1}\in\partial\mathbb{B}_{\mathrm{dim}(m)}$,
$t_{2},t_{3}\in[-1,1]$, there exist sensitivity parameters $R_{y\sim u\mid
a,m,c}$, $R_{m\sim u\mid a,c}$, $R_{a\sim u\mid c}$ satisfying $R_{y\sim u\mid
a,m,c}^{2}\leq\rho_{y}$, $\|R_{m\sim u\mid a,c}\|_{2}^{2}\leq\rho_{m}$,
$R_{a\sim u\mid c}^{2}\leq\rho_{a}$, and
$\phi_{3}=\\{\rho_{m}\rho_{a}/(1-\rho_{a})\\}^{1/2}t_{2}\cdot t_{1}$,
$\phi_{4}=\\{\rho_{y}\rho_{m}/(1-\rho_{m})\\}^{1/2}t_{3}\cdot t_{1}$. We can
choose parameters $R_{m\sim u\mid a,c}=\rho_{m}^{1/2}\cdot t_{1}$, $R_{y\sim
u\mid a,m,c}=\rho_{y}^{1/2}t_{3}$, and

$R_{a\sim u\mid
c}=\left(\frac{\rho_{a}}{t_{2}^{2}\rho_{a}-\rho_{a}+1}\right)^{1/2}t_{2},$

which satisfy the conditions.

We next discuss the case when $u$ is a vector. Similar to the argument in the
proof of Theorem 3, we can show that ${\phi}_{3}$ and ${\phi}_{4}$ are
invariant to invertible linear transformation of $u$. Moreover, since
$\|R_{y\sim u\mid a,m,c}\|_{2}^{2}$, $\|R_{m\sim u\mid a,c}\|_{2}^{2}$,
$\|R_{a\sim u\mid c}\|_{2}^{2}$ are invariant to invertible linear
transformation of $u$, we can assume $\mathrm{cov}(u_{\perp a,c})=I$ and find
the range of

$\begin{split}\tilde{\phi}_{3}&=R_{m\sim u\mid a,c}\left(I-R_{a\sim u\mid
c}^{\mathrm{\scriptscriptstyle T}}R_{a\sim u\mid c}\right)^{-1/2}R_{a\sim
u\mid c}^{\mathrm{\scriptscriptstyle T}},\\\ \tilde{\phi}_{4}&=R_{m\sim u\mid
a,c}\left(I-R_{m\sim u\mid a,c}^{\mathrm{\scriptscriptstyle T}}R_{m\sim u\mid
a,c}\right)^{-1/2}R_{y\sim u\mid a,m,c}^{\mathrm{\scriptscriptstyle
T}}\end{split}$

under the constraint that $\|R_{y\sim u\mid a,m,c}\|_{2}^{2}\leq\rho_{y}$,
$\|R_{m\sim u\mid a,c}\|_{2}^{2}\leq\rho_{m}$, $\|R_{a\sim u\mid
c}\|_{2}^{2}\leq\rho_{a}$. Then, we have

$\begin{split}\|\tilde{\phi}_{3}\|_{2}&\leq\left(\frac{\rho_{m}\rho_{a}}{1-\rho_{a}}\right)^{1/2},\\\
\|\tilde{\phi}_{4}\|_{2}&\leq\left(\frac{\rho_{y}\rho_{m}}{1-\rho_{m}}\right)^{1/2}.\end{split}$

We just need to show that for any
$t_{1},t_{2}\in\overline{\mathbb{B}}_{\mathrm{dim}(m)}$, there exist
sensitivity parameters $R_{y\sim u\mid a,m,c}$, $R_{m\sim u\mid a,c}$,
$R_{a\sim u\mid c}$ satisfying $\|R_{y\sim u\mid
a,m,c}\|_{2}^{2}\leq\rho_{y}$, $\|R_{m\sim u\mid a,c}\|_{2}^{2}\leq\rho_{m}$,
$\|R_{a\sim u\mid c}\|_{2}^{2}\leq\rho_{a}$, and
$\tilde{\phi}_{3}=\\{\rho_{m}\rho_{a}/(1-\rho_{a})\\}^{1/2}\cdot t_{1}$,
$\tilde{\phi}_{4}=\\{\rho_{y}\rho_{m}/(1-\rho_{m})\\}^{1/2}\cdot t_{2}$. Let
$\\{v_{1},v_{2}\\}$ be the orthogonal unit vectors in $\mathbb{R}^{m}$ such
that $t_{1},t_{2}$ are in the linear space spanned by $\\{v_{1},v_{2}\\}$. Let
$t_{1}=\lambda_{11}v_{1}+\lambda_{12}v_{2}$ and
$t_{2}=\lambda_{21}v_{1}+\lambda_{22}v_{2}$. Then, we have
$\lambda_{11}^{2}+\lambda_{12}^{2}\leq 1$ and
$\lambda_{21}^{2}+\lambda_{22}^{2}\leq 1$. We can choose

$\begin{split}R_{m\sim u\mid
a,c}&=\rho_{m}^{1/2}\begin{pmatrix}v_{1}&v_{2}&0&\cdots&0\end{pmatrix},\\\
R_{y\sim u\mid
a,m,c}&=\rho_{y}^{1/2}\begin{pmatrix}\lambda_{21}&\lambda_{22}&0&\cdots&0\end{pmatrix},\\\
R_{a\sim u\mid
c}&=\left\\{\frac{\rho_{a}}{\rho_{a}(\lambda_{11}^{2}+\lambda_{12}^{2})-\rho_{a}+1}\right\\}^{1/2}\begin{pmatrix}\lambda_{11}&\lambda_{12}&0&\cdots&0\end{pmatrix},\\\
\end{split}$

which satisfy the conditions. Here we use the identity

$(I-ww^{\mathrm{\scriptscriptstyle
T}})^{-1/2}w=\frac{w}{(1-\|w\|_{2}^{2})^{1/2}}$

for any $w\in\mathbb{B}_{2}$.

(iii) It follows from (i) and (ii). ∎

###### Proof of Theorem S1.

The proof of (i) is in Cinelli and Hazlett (2020). For (ii), for any
$R_{1},R_{2}\in(-1,1)^{2}$, Lemma S5 implies that

$\Sigma_{1}=\begin{pmatrix}\hat{\mathrm{var}}(A_{\perp
C})&R_{1}\hat{\mathrm{var}}(A_{\perp C})^{1/2}\hat{\mathrm{var}}(U_{\perp
C})^{1/2}\\\ R_{1}\hat{\mathrm{var}}(U_{\perp
C})^{1/2}\hat{\mathrm{var}}(A_{\perp C})^{1/2}&\hat{\mathrm{var}}(U_{\perp
C})\end{pmatrix}$

and

$\Sigma_{2}=\begin{pmatrix}\hat{\mathrm{var}}(Y_{\perp
A,C})&R_{2}\hat{\mathrm{var}}(Y_{\perp A,C})^{1/2}\hat{\mathrm{var}}(U_{\perp
A,C})^{1/2}\\\ R_{2}\hat{\mathrm{var}}(U_{\perp
A,C})^{1/2}\hat{\mathrm{var}}(Y_{\perp A,C})^{1/2}&\hat{\mathrm{var}}(U_{\perp
A,C})\end{pmatrix}.$

are positive definite. By Lemma S6, there exists $U$ such that
$R_{1}=\hat{\mathrm{corr}}(A_{\perp C},U_{\perp C})$ and
$R_{2}=\hat{\mathrm{corr}}(Y_{\perp A,C},U_{\perp A,C})$. Thus,
$\\{\hat{\mathrm{corr}}(U_{\perp C},A_{\perp C}),\hat{\mathrm{corr}}(U_{\perp
A,C},Y_{\perp A,C})\\}$ can take arbitrary values in $(-1,1)^{2}$. Since
$R_{Y\sim U\mid A,C}^{2}=\hat{\mathrm{corr}}(U_{\perp A,C},Y_{\perp
A,C})^{2}$, $R_{A\sim U\mid C}^{2}=\hat{\mathrm{corr}}(U_{\perp C},A_{\perp
C})^{2}$, and
$\mathrm{sgn}(\hat{\gamma}_{1}^{\text{unobs}}-\hat{\gamma}^{\text{obs}}_{1})=-\mathrm{sgn}\\{\hat{\mathrm{corr}}(Y_{\perp
A,C},U_{\perp A,C})\\}\cdot\mathrm{sgn}\\{\hat{\mathrm{corr}}(U_{\perp
C},A_{\perp C})\\}$, the parameters $\\{R_{Y\sim U\mid A,C}^{2},R_{A\sim U\mid
C}^{2}\\}$ can take arbitrary values in $[0,1)^{2}$, and $\\{R_{Y\sim U\mid
A,C}^{2},R_{A\sim U\mid C}^{2}$,
$\mathrm{sgn}(\hat{\gamma}_{1}^{\text{unobs}}-\hat{\gamma}^{\text{obs}}_{1})\\}$
can take arbitrary values in $(0,1)^{2}\times\\{-1,1\\}$. ∎

###  E.3 Proofs of Propositions 1–3 and S1–S2

Propositions 1–3 are direct corollaries of Lemma S3.

###### Proof of Proposition S1.

By Lemma S1,

$R_{A\sim U\mid M,C}^{2}=R_{U\sim A\mid M,C}^{2}=1-\frac{(1-R_{U\sim M\mid
A,C}^{2})(1-R_{U\sim A\mid C}^{2})}{1-R_{U\sim M\mid
C}^{2}}=1-\frac{(1-R_{M\sim U\mid A,C}^{2})(1-R_{A\sim U\mid
C}^{2})}{1-R_{M\sim U\mid C}^{2}}$

Then, we need to represent $R_{M\sim U\mid C}^{2}$ using natural $R^{2}$’s. By
Lemma S4,

$\begin{split}R_{M\sim U\mid C}^{2}&=\left[S(R_{M\sim A\mid C}^{2}R_{U\sim
A\mid C}^{2})^{1/2}+S^{\prime}\left\\{R_{M\sim U\mid A,C}^{2}(1-R_{M\sim A\mid
C}^{2})(1-R_{U\sim A\mid C}^{2})\right\\}^{1/2}\right]^{2}\\\
&=(Sv_{1}+S^{\prime}v_{2})^{2}\\\ \end{split}$

where $S=\mathrm{sgn}\\{\mathrm{corr}(m_{\perp c},a_{\perp
c})\\}\cdot\mathrm{sgn}\\{\mathrm{corr}(u_{\perp c},a_{\perp c})\\}$ and
$S^{\prime}=\mathrm{sgn}\\{\mathrm{corr}(m_{\perp a,c},u_{\perp a,c})\\}$. We
discuss the following cases.

  * •

If $R_{M\sim U\mid A,C}^{2}=0$, then $v_{2}=0$. We have

$R_{A\sim U\mid M,C}^{2}=1-\frac{1-R_{A\sim U\mid C}^{2}}{1-v_{1}^{2}}.$

  * •

If $R_{M\sim U\mid A,C}^{2}\not=0$, then $v_{2}\not=0$. By Cochran’s formula,
we have

$R_{A\sim U\mid M,C}^{2}=1-\frac{(1-R_{A\sim U\mid C}^{2})(1-R_{M\sim U\mid
A,C}^{2})}{1-(S_{2}v_{1}+v_{2})^{2}}.$

∎

Proposition S2 is a direct corollary of Lemma S3.

###  E.4 Proofs of Corollaries 1–3 and S1–S2

Since Corollary 2 is a general version of Corollary 1, we only provide a proof
of Corollary 2.

###### Proof of Corollary 2.

For any $R_{1},R_{3}\in(-1,1)$ and $R_{2}\in\mathbb{B}_{\mathrm{dim}(m)}$,
Lemma S5 implies that

$\Sigma_{1}=\begin{pmatrix}\mathrm{var}(a_{\perp
c})&R_{1}\mathrm{var}(a_{\perp c})^{1/2}\mathrm{var}(u_{\perp c})^{1/2}\\\
R_{1}\mathrm{var}(u_{\perp c})^{1/2}\mathrm{var}(a_{\perp
c})^{1/2}&\mathrm{var}(u_{\perp c})\end{pmatrix},$
$\Sigma_{2}=\begin{pmatrix}\mathrm{cov}(m_{\perp a,c})&\mathrm{cov}(m_{\perp
a,c})^{1/2}R_{2}\mathrm{var}(u_{\perp a,c})^{1/2}\\\ \mathrm{var}(u_{\perp
a,c})^{1/2}R_{2}^{\mathrm{\scriptscriptstyle T}}\mathrm{cov}(m_{\perp
a,c})^{1/2}&\mathrm{var}(u_{\perp a,c})\end{pmatrix},$

and

$\Sigma_{3}=\begin{pmatrix}\mathrm{var}(y_{\perp
a,m,c})&R_{3}\mathrm{var}(y_{\perp a,m,c})^{1/2}\mathrm{var}(u_{\perp
a,m,c})^{1/2}\\\ R_{3}\mathrm{var}(u_{\perp a,m,c})^{1/2}\mathrm{var}(y_{\perp
a,m,c})^{1/2}&\mathrm{var}(u_{\perp a,m,c})\end{pmatrix}$

are positive definite. By Lemma S6, there exists $u$ such that $R_{1}=R_{a\sim
u\mid c}$, $R_{2}=R_{m\sim u\mid a,c}$, $R_{3}=R_{y\sim u\mid a,m,c}$. This
shows that $\\{R_{y\sim u\mid a,m,c},R_{m\sim u\mid a,c},R_{a\sim u\mid c}\\}$
can take arbitrary values in
$(-1,1)\times\mathbb{B}_{\mathrm{dim}(m)}\times(-1,1)$. ∎

The proof of Corollary 3 is similar to the proof of Corollary 2. We omit its
proof.

###### Proof of Corollary S1.

For any $R_{1},R_{2},R_{3}\in(-1,1)$, Lemma S5 implies that

$\Sigma_{1}=\begin{pmatrix}\hat{\mathrm{var}}(A_{\perp
C})&R_{1}\hat{\mathrm{var}}(A_{\perp C})^{1/2}\hat{\mathrm{var}}(U_{\perp
C})^{1/2}\\\ R_{1}\hat{\mathrm{var}}(U_{\perp
C})^{1/2}\hat{\mathrm{var}}(A_{\perp C})^{1/2}&\hat{\mathrm{var}}(U_{\perp
C})\end{pmatrix},$ $\Sigma_{2}=\begin{pmatrix}\hat{\mathrm{var}}(M_{\perp
A,C})&R_{2}\hat{\mathrm{var}}(M_{\perp A,C})^{1/2}\hat{\mathrm{var}}(U_{\perp
A,C})^{1/2}\\\ R_{2}\hat{\mathrm{var}}(U_{\perp
A,C})^{1/2}\hat{\mathrm{var}}(M_{\perp A,C})^{1/2}&\hat{\mathrm{var}}(U_{\perp
A,C})\end{pmatrix},$

and

$\Sigma_{3}=\begin{pmatrix}\hat{\mathrm{var}}(Y_{\perp
A,M,C})&R_{3}\hat{\mathrm{var}}(Y_{\perp
A,M,C})^{1/2}\hat{\mathrm{var}}(U_{\perp A,M,C})^{1/2}\\\
R_{3}\hat{\mathrm{var}}(U_{\perp A,M,C})^{1/2}\hat{\mathrm{var}}(Y_{\perp
A,M,C})^{1/2}&\hat{\mathrm{var}}(U_{\perp A,M,C})\end{pmatrix}$

are positive definite. By Lemma S7, there exists $U$ such that
$R_{1}=\hat{\mathrm{corr}}(A_{\perp C},U_{\perp C})$,
$R_{2}=\hat{\mathrm{corr}}(M_{\perp A,C},U_{\perp A,C})$, and
$R_{3}=\hat{\mathrm{corr}}(Y_{\perp A,M,C},U_{\perp A,M,C})$. This shows that
$\\{\hat{\mathrm{corr}}(A_{\perp C},U_{\perp C})$,
$\hat{\mathrm{corr}}(M_{\perp A,C},U_{\perp A,C})\\}$ can take arbitrary
values in $(-1,1)^{2}$. Moreover, since $\hat{\mathrm{corr}}(A_{\perp
C,M},U_{\perp C,M})$ is independent of the selection of $R_{3}$, given
$\\{\hat{\mathrm{corr}}(A_{\perp C},U_{\perp C})$,
$\hat{\mathrm{corr}}(M_{\perp A,C},U_{\perp A,C})$,
$\hat{\mathrm{corr}}(A_{\perp C,M},U_{\perp C,M})\\}$, the
$\hat{\mathrm{corr}}(Y_{\perp A,C,M},U_{\perp A,C,M})$ can take arbitrary
values in $(-1,1)$.

Therefore, $R_{U}^{2}$ can take arbitrary values in $[0,1)^{3}$. Moreover, by
Cochran’s formula,

$\begin{split}\mathrm{sgn}(\hat{\theta}_{1}^{\text{unobs}}-\hat{\theta}_{1}^{\text{obs}})&=-\mathrm{sgn}\\{\hat{\mathrm{corr}}(Y_{\perp
A,C,M},U_{\perp A,C,M})\\}\cdot\mathrm{sgn}\\{\hat{\mathrm{corr}}(A_{\perp
C,M},U_{\perp C,M})\\},\\\
\mathrm{sgn}(\hat{\beta}_{1}^{\text{unobs}}-\hat{\beta}_{1}^{\text{obs}})&=-\mathrm{sgn}\\{\hat{\mathrm{corr}}(M_{\perp
A,C},U_{\perp A,C})\\}\cdot\mathrm{sgn}\\{\hat{\mathrm{corr}}(A_{\perp
C},U_{\perp C})\\},\end{split}$

which implies that, if $\hat{\beta}^{\text{obs}}\not=0$, then
$\\{R_{U}^{2},S_{1},S_{2}\\}$ can take arbitrary values in
$(0,1)^{3}\times\\{-1,1\\}^{2}\setminus E_{0}$. ∎

###### Proof of Corollary S2.

We have shown in the proof of Corollary S1 that

$\\{\hat{\mathrm{corr}}(A_{\perp C},U_{\perp C}),\hat{\mathrm{corr}}(A_{\perp
C,M},U_{\perp C,M}),\hat{\mathrm{corr}}(Y_{\perp A,C,M},U_{\perp A,C,M})\\}$

can take arbitrary values in $(-1,1)^{3}$. By Cochran’s formula,

$\begin{split}\mathrm{sgn}(\hat{\beta}_{1}^{\text{unobs}}-\hat{\beta}_{1}^{\text{obs}})&=-\mathrm{sgn}\\{\hat{\mathrm{corr}}(M_{\perp
A,C},U_{\perp A,C})\\}\cdot\mathrm{sgn}\\{\hat{\mathrm{corr}}(A_{\perp
C},U_{\perp C})\\}\\\
\mathrm{sgn}(\hat{\theta}_{3}^{\text{unobs}}-\hat{\theta}_{3}^{\text{obs}})&=-\mathrm{sgn}\\{\hat{\mathrm{corr}}(Y_{\perp
A,C,M},U_{\perp A,C,M})\\}\cdot\mathrm{sgn}\\{\hat{\mathrm{corr}}(M_{\perp
A,C},U_{\perp A,C})\\}.\\\ \end{split}$

Therefore, $\\{R_{U}^{2},S_{3},S_{4}\\}$ can take arbitrary values in
$(0,1)^{3}\times\\{-1,1\\}^{2}$. ∎


# BEDRF: Bidirectional Edge Diffraction Response Function for Interactive
Sound Propagation

Chunxiao Cao<EMAIL_ADDRESS>State Key Laboratory of CAD&CG, Zhejiang
UniversityHangzhouChina310058 ,  Zili An<EMAIL_ADDRESS>State Key
Laboratory of CAD&CG, Zhejiang University ,  Zhong Ren<EMAIL_ADDRESS>State Key Laboratory of CAD&CG, Zhejiang University ,  Dinesh Manocha
<EMAIL_ADDRESS>Department of Computer Science, University of Maryland8125
Paint Branch DriveCollege ParkMarylandUSA20742 and  Kun Zhou
<EMAIL_ADDRESS>State Key Laboratory of CAD&CG, Zhejiang University

###### Abstract.

We introduce bidirectional edge diffraction response function (BEDRF), a new
approach to model wave diffraction around edges with path tracing. The
diffraction part of the wave is expressed as an integration on path space, and
the wave-edge interaction is expressed using only the localized information
around points on the edge similar to a bidirectional scattering distribution
function (BSDF) for visual rendering. For an infinite single wedge, our model
generates the same result as the analytic solution . Our approach can be
easily integrated into interactive geometric sound propagation algorithms that
use path tracing to compute specular and diffuse reflections. Our resulting
propagation algorithm can approximate complex wave propagation phenomena
involving high-order diffraction, and is able to handle dynamic, deformable
objects and moving sources and listeners. We highlight the performance of our
approach in different scenarios to generate smooth auralization.

##  1\. Introduction

Interactive sound simulation and rendering is becoming increasingly popular in
games and virtual environments. Accurate simulation of sound phenomena can
significantly improve the plausibility of a VR environment as well as the
immersive experience for the user. The most accurate methods for sound
simulation are based on wave-based methods (Pind Jörgensson, 2020) . However,
a physically accurate simulation of sound propagation, using a numerical wave
equation solver, is computationally-intensive. Despite various efforts to
increase the efficiency of the wave equation solvers (Raghuvanshi et al.,
2009; Savioja, 2010; Mehra et al., 2012; Allen and Raghuvanshi, 2015), they
are not fast enough for interactive applications, especially when dealing with
complex scenes with moving objects.

For interactive simulation, a large family of algorithms, based on the theory
of geometric acoustics (GA), has been developed for sound propagation.
Geometric optics (GO), the counterpart of GA in the field of computer
graphics, has been the foundation of photorealistic rendering. GO expresses
light transport with an integral called the transport equation or the
rendering equation (Kajiya, 1986). This equation has been extensively studied,
and many efficient integrators have been developed and used in various
applications (Veach and Guibas, 1994, 1997; Keller, 1997; Jensen, 2001).

Many researchers have borrowed techniques from GO and have extended them to
the field of acoustics. However, the accuracy of geometric approximation
depends on the wavelength. One of the reasons for the success of GO in light
propagation simulation is that the wavelength of visible light in a vacuum
ranges from $380\text{nm}$ to $740\text{nm}$, a very small value compared to
the scale of common objects. In comparison, the wavelength of audible sound
waves in the air ranges from $0.017\text{m}$ to $17\text{m}$, which makes
wave-based sound phenomena easily perceptible by human ears. Among these
phenomena, a key issue is diffraction or occlusion effects, which tend to be
important for receivers located inside the shadow region of an obstacle. Due
to the discontinuity of the visibility function, the ray-based propagation
model will generate “shadows”. Without simulating the diffraction effect, the
user would perceive a sudden change of sound quality when entering a shadow
zone, which contradicts our aural experiences in the physical world. It has
been shown that inaccurate modeling of diffraction effects can result in loss
of realism in virtual environments (Rungta et al., 2016).

Since diffraction is a low-frequency effect, it is especially difficult to
simulate for interactive GA-based algorithms based on ray or path tracing.
Many techniques have been proposed to approximate diffraction effects with
methods based on the unified theory of diffraction (UTD, (Tsingos et al.,
2001; Taylor et al., 2012; Schissler et al., 2014)), Biot-Tolstoy-Medwin model
(BTM, (Calamia and Svensson, 2007; Calamia, 2009; Antani et al., 2010)),
diffraction kernels (Rungta et al., 2018), learning methods (Tang et al.,
2021), the Heisenberg uncertainty principle (Stephenson and Svensson, 2007),
and volumetric diffraction and transmission (Pisha et al., 2020). However,
their accuracy can vary based on the environment, moving objects, and the
location of the source or the receiver relative to the diffraction edges.

Main Results: In this paper, we present a new method to compute the
diffraction effects in an interactive GA framework, with the following
contribution:

  * •

We present BEDRF, a localized edge-diffraction representation around convex
wedges. Our novel representation relies only on the incident and outgoing wave
direction and the local boundary condition (wedge angle and material). This is
similar to a bidirectional scattering distribution function (BSDF) (Asmail,
1991) for visual rendering. Our new representation is derived from the exact
solution for planar incident waves. We also present an importance sampling
strategy for this representation, which is used for interactive propagation.

  * •

A novel, efficient edge diffraction solver based on Monte Carlo path tracing.
We present a rendering equation for diffraction and an integrator that solves
this rendering equation. Our integrator supports the simulation of diffraction
effects via the BEDRF representation. With the concept of “meta-path space”,
our modified intersection detector allows paths to interact with edges in the
scene, while being compatible with traditional path tracing techniques like
multiple importance sampling (MIS). We compared the accuracy of our solver
with state-of-the-art diffraction algorithms based on UTD and BTM.

  * •

Several specially designed techniques that improve the accuracy and
performance of our diffraction path tracer. The effect and performance of
these techniques are demonstrated with simulation results in actual scenes.

We combine our approach with an interactive geometric propagation approach
that uses path tracing for specular and diffuse reflections. We highlight its
benefits in different benchmarks in terms of impulse response (IR) computation
and auralization.

##  2\. Related Works

###  2.1. Sound Propagation and Geometric Acoustics

The existing algorithms for sound propagation simulation can be divided into
two categories (Liu and Manocha, 2020). The first category is “wave-based”
algorithms that compute the propagation result with a numerical wave equation
solver. The solver may be based on the finite-element method (FEM) (Thompson,
2006), the boundary-element method (BEM) (Kirkup, 2019) or the finite-
difference time domain method (FDTD) (Hamilton, 2021). These methods are
accurate, but usually very slow for complex acoustic scenes. Moreover, their
complexity can increase as a fourth power of the simulation frequency. The
second category is GA-based algorithms, which can be accurate at high
frequencies, have been shown to be effective in large, dynamic scenes, and are
used for interactive applications.

While earlier algorithms like the image source method (Vorländer, 1989) use
specific transport equations, many current GA algorithms are based on the
acoustic rendering equation in (Siltanen et al., 2007). Since this equation is
very similar to the rendering equation for visual rendering, many acoustic
solvers take inspiration from existing visual rendering algorithms like Monte
Carlo path tracing (Lentz et al., 2007; Taylor et al., 2012; Schissler et al.,
2014; Cao et al., 2016; Schissler and Manocha, 2016) or photon mapping
(Bertram et al., 2005; Kapralos et al., 2007).

###  2.2. Diffraction Algorithms in Geometric Acoustics

Many techniques have been proposed to approximate wave diffraction within the
framework of GA. A direct solution is to use a hybrid framework that combines
the numeric and GA solution (Yeh et al., 2013; Mehra et al., 2013; Rungta et
al., 2018). The wave propagation around scene objects is precomputed with a
numerical solver, and the long-distance propagation in empty areas is modeled
with GA. Some of these methods are limited to static scenes, and others only
consider the local scattering field around the objects.

Another approach is to analyze the diffraction wave around local geometry of
the objects, and combine the local solutions into a global one. As scene
geometries are usually represented with triangle meshes, the edge-diffraction
models, which give the diffraction solution around the intersecting edge of
two half-planes forming a wedge, can be useful. Most edge-diffraction models
can only compute the solution for convex wedges, as solutions for non-convex
wedges are complicated due to the the presence of interreflection between two
half-planes and are unsuitable for real-life applications. The most popular
models in sound propagation are the uniform geometrical theory of diffraction
(UTD) and the Biot-Tolstoy-Medwin model (BTM).

UTD (Kouyoumjian and Pathak, 1974; Mcnamara et al., 1990) is widely used for
simulating diffraction effects of electromagnetic waves. There have also been
some successful applications in simulating interactive sound
propagation(Tsingos et al., 2001; Taylor et al., 2012; Schissler et al., 2014,
2021). The UTD model gives the solution in various conditions, such as
spherical or planar incident wave and flat or curved surface. However, it also
suffers from many issues, including the fact that its solution is non-exact,
especially for finite-length edges, Moreover, its solution is expressed in the
frequency domain instead of the time domain. The last point has a particular
impact on the resulting quality, as the solvers in actual applications usually
simulate wave propagation on 3-8 frequency ranges, and are unable to process
phase information.

BTM (Svensson et al., 1999) was developed by combining the Biot-Tolstoy
solution(Biot and Tolstoy, 1957) with Medwin’s assumption(Medwin et al.,
1982). The Biot-Tolstoy solution gives the exact description of the
propagation of a spherical wave around a wedge, and Medwin’s assumption
expresses such a solution with “secondary sources” on the edge of the wedge,
which allows BTM to deal with finite-edge diffraction. BTM requires
integration operations on geometry edges, which makes it difficult to combine
with path tracing (integrating on the path space). In current auralizers using
BTM (Calamia and Svensson, 2007; Calamia, 2009; Antani et al., 2010), these
integrations are computed separately and are used for non-interactive
applications. Since path tracing also involves numerical integration, we
reformulate the integration on the path space, which allows us to use the
integrator for rendering equation to solve the diffraction problem. We use
this idea to develop the notion of BEDRF and use our approach for interactive
propagation.

There are also some other solutions for GA diffraction. A very promising
approximative model for edge diffraction is the directive line source model
(DLSM) (Menounou et al., 2000). An improved version of DLSM (Menounou and
Nikolaou, 2017) gives the exact solution for plane wave diffraction around a
half-plane and works well in many other simple cases. (Stephenson and
Svensson, 2007; Stephenson, 2010) compute convex edge diffraction with the
Heisenberg uncertainty principle and model edge diffraction with a diffraction
angle probability density function (DAPDF). This method is conceptually
simple, but lacks a solid physical basis and is non-exact. (Tsingos et al.,
2007) combines the Kirchhoff approximation and GPU rasterization to compute
the first-order diffraction wave. This method is valid for arbitary geometric
shapes, but difficult to extend to high-order diffraction. The volumetric
diffraction and transmission (VDaT) method (Pisha et al., 2020) is an
approximation of BTM, which uses a large number of probe rays and simplifies
all occluding geometries into a combination of planar rings. VDaT is unable to
express high-order diffraction, and the accuracy of this approach depends on
the similarity of the original and simplified representations.

##  3\. Background

The notation used in our paper mostly come from (Siltanen et al., 2007) and
(Keller and Blank, 1951), with a few modification. We highlight the symbols
used in this paper in Table 1. We also use the “arrow notations” to simplify
path-tracing-related expressions.

Symbol | Explanation
---|---
$\mathbb{R}$ | set of real numbers
$||\mathbf{x}||_{2}$ | 2-norm (Euclidian length) of vector $x$
$\nabla x$ | gradient of $x$
$E[x]$ | expectation of random variable $x$
$\mathcal{F}(f)$ | Fourier transform of function $f$
$L(\mathbf{x},\mathbf{v}_{i},t)$ | sound wave at $\mathbf{x}$,
| received from the direction $\mathbf{v}_{i}$ at time $t$
$L_{0}(\mathbf{x},\mathbf{v}_{o},t)$ | initial wave from source $\mathbf{x}$,
| to the direction $\mathbf{v}_{o}$ at time $t$
$M(\mathbf{x},\mathbf{x}^{\prime},t)$ | media term
$V(\mathbf{x},\mathbf{x}^{\prime})$ | visibility between point $\mathbf{x}$ and $\mathbf{x}^{\prime}$,
| 0 when invisible, 1 when visible
$\rho(\mathbf{x},\mathbf{v}_{i},\mathbf{v}_{o})$ | BSDF (Sect. 3.2) / BEDRF (Sect. 4) at $\mathbf{x}$,
| with incident direction $\mathbf{v}_{i}$
| and outgoing direction $\mathbf{v}_{o}$
$c$ | sound speed
$\varphi_{i},\varphi_{o},\theta_{i},\theta_{o}$ | local coordinate angles, see Fig. 2
Arrow Notation | Expanded Form
$L(\mathbf{x}^{\prime}\rightarrow\mathbf{x},t)$ | $L(\mathbf{x},\frac{\mathbf{x}-\mathbf{x}^{\prime}}{||\mathbf{x}-\mathbf{x}^{\prime}||_{2}},t)$
$L_{0}(\mathbf{x}^{\prime\prime}\rightarrow\mathbf{x}^{\prime},t)$ | $L_{0}(\mathbf{x}^{\prime\prime},\frac{\mathbf{x}^{\prime}-\mathbf{x}^{\prime\prime}}{||\mathbf{x}^{\prime}-\mathbf{x}^{\prime\prime}||_{2}},t)$
$\rho(\mathbf{x}^{\prime\prime}\rightarrow\mathbf{x}^{\prime}\rightarrow\mathbf{x})$ | $\rho(\mathbf{x}^{\prime},\frac{\mathbf{x}^{\prime}-\mathbf{x}^{\prime\prime}}{||\mathbf{x}^{\prime}-\mathbf{x}^{\prime\prime}||_{2}},\frac{\mathbf{x}-\mathbf{x}^{\prime}}{||\mathbf{x}-\mathbf{x}^{\prime}||_{2}})$
Table 1. Common symbols used in this paper.

###  3.1. Geometric Acoustics

The theory of geometric acoustics is based on the Eikonal approximation, which
assumes that the solution for the wave equation $f$ takes the following form:

(1) $f(\mathbf{x},t)=e^{A(\mathbf{x})+i(\frac{1}{\lambda}B(\mathbf{x})-ct)},$

where $\lambda$ is the wave length and $A(\mathbf{x})$ and $B(\mathbf{x})$
express the amplitude and phase delay at different positions. Combined with
the wave equation in homogenous media

(2) $\frac{\partial^{2}f}{\partial t^{2}}=c^{2}\nabla^{2}f,$

we have the following equation:

(3) $||\nabla B||_{2}^{2}-\frac{1}{\lambda^{2}}=\nabla^{2}A+||\nabla
A||^{2}_{2}.$

If $\lambda$ is sufficiently small, the right side of the equation can be
ignored, and we have

(4) $||\nabla B||_{2}^{2}\approx\frac{1}{\lambda^{2}},$

which indicates that $\nabla B$ is of constant length, and every trajectory
generated by this vector field is a straight line. In this way, the
propagation of a sound wave can be modeled with wave packets of different
amplitudes and frequencies, travelling on paths (sometimes called “rays”) with
direction $\nabla B$ at a constant speed. With this approximation, the
original wave propagation result expressed as a solution of a differential
equation becomes an integration on the path space.

###  3.2. Rendering Equation

The purpose of a geometric wave propagation algorithm is to solve the
rendering equation, Since we are dealing with sound waves, we use the time-
dependent rendering equation proposed by (Siltanen et al., 2007). A
reformulated version is presented below:

(5) $\displaystyle L(\mathbf{x}^{\prime}\rightarrow\mathbf{x},t)$
$\displaystyle=L_{0}(\mathbf{x}^{\prime}\rightarrow\mathbf{x},t)$
$\displaystyle+\int_{\Omega}\big{(}L(\mathbf{x}^{\prime}\rightarrow\mathbf{x}^{\prime\prime},t)V(\mathbf{x}^{\prime\prime},\mathbf{x})\rho(\mathbf{x}^{\prime}\rightarrow\mathbf{x}^{\prime\prime}\rightarrow\mathbf{x})\big{.}$
$\displaystyle\qquad\big{.}*M(\mathbf{x}^{\prime\prime},\mathbf{x},t)\big{)}dA_{\mathbf{x}^{\prime\prime}},$

In Eq. 5. $\Omega$ represents the scene geometry.
$\rho(\mathbf{x}^{\prime\prime},\mathbf{v}_{i},\mathbf{v}_{o})$ is the
bidirectional scattering distribution function (BSDF) describing the
interaction between the wave and the scene at $\mathbf{x}^{\prime\prime}$. The
media term $M(\mathbf{x}^{\prime\prime},\mathbf{x},t)$ represents the effects
of the medium on the wave during the propagation from
$\mathbf{x}^{\prime\prime}$ to $\mathbf{x}$, including time delay, propagation
decay and energy absorption. Its form is usually given like this:

(6) $M(\mathbf{x}^{\prime\prime},\mathbf{x},t)=\frac{e^{-\alpha
r}\delta(t-\frac{r}{c})}{r^{2}},r=||\mathbf{x}-\mathbf{x}^{\prime\prime}||_{2},$

where $\alpha$ is the absorption factor of the media and $\delta(t)$ is the
Dirac delta function. The media term is applied on the left side of the wave
with a convolution (denoted by the asterisk $*$). $A$ is the area measure
defined on $\Omega$, and $dA_{\mathbf{x}^{\prime\prime}}$ in Eq. 5 indicates
that we integrate on the measure $A$ with respect to the variable
$\mathbf{x}^{\prime\prime}$. In our algorithm, $A$ is extended to measure both
the area of the geometry and the length of the convex edges on it. The media
term is also modified to suit our amplitude-based framework, which will be
discussed in Sect. 3.2.1.

The transport equation is easier to understand in its operator form:

(7) $L=L_{0}+TL,$

where

(8) $TL=\int_{\Omega}(LV\rho*M)dA_{\mathbf{x}^{\prime\prime}}.$

The operator $T$ has an intuitive explanation as an unoccluded propagation of
wave from a point of the scene $\Omega$ to another. We can see that the
rendering equation is a formal expression of the fact that the total wave
equals the source wave plus the wave induced by wave-geometry interaction.

Expanding Eq. 7 into a Neumann series (Kreyszig, 1978), we have

(9) $L=\sum_{i=0}^{\infty}T^{i}L_{0}.$

This is the expression used in actual algorithms like path tracing.

####  3.2.1. Amplitude-Based Rendering Equation

The media term in Eq. 6 decays in an inverse-square way relative to $r$, which
matches the inverse-square law for intensity. However, for an amplitude-based
rendering equation, we expect its solution to behave like one for the wave
equation. The general solution of spherical waves is given below:

(10)
$f(\mathbf{x},t)=\sum_{l,m}Y_{lm}J_{l}(r)(a_{lm}\delta(t-\frac{r}{c})+b_{lm}\delta(t+\frac{r}{c}))$

where $Y_{lm}$ is the $m$-th spherical harmonic function of degree $l$,
$J_{l}$ is the spherical Bessel function of degree $l$, and $a_{lm}$ and
$b_{lm}$ are constants. When $r$ is sufficiently large, we have
$|J_{l}|\approx 1/r$ for every $l$. Thus we need to change the denominator
$r^{2}$ in Eq. 6 to $r$.

###  3.3. Monte Carlo Path Tracing

A common algorithm that solves the integration in Eq. 8 is Monte Carlo
integration, which is a probabilistic method using the fact that the
mathematical expectation of a random variable is an integration on the
probability space:

(11)
$E\left[\frac{f(X)}{p(X)}\right]=\int_{\Omega}f(x)\frac{d\mu_{1}}{d\mu_{2}}d\mu_{2}=\int_{\Omega}f(x)d\mu_{1},$

where $p(X)$ is the Radon-Nikodym derivative (Royden and Fitzpatrick, 1988)
$d\mu_{2}/d\mu_{1}$ between the two measures $\mu_{1}$ and $\mu_{2}$ on the
space $\Omega$. $\mu_{1}$ is usually a natural measure defined on $\Omega$,
like $A$ in Eq. 5. $\mu_{2}$ is the sampling probability measure that depends
on the algorithm.

To calculate $T^{i}L_{0}$ in Eq. 9, we need to take samples from $\Omega^{i}$.
Each sample is a path that connects the source and the listener over $i$
interactions with the scene geometry. Each segment of the path corresponds to
the media term $M$ in Eq. 5, and each vertex corresponds to the BSDF $\rho$.
In path tracing, the generation process of each segment of the path is
mutually independent. Therefore. the $f/p$ in Eq. 11 could be expanded into
the following form:

(12)
$\frac{f}{p}=\frac{L_{0}(\mathbf{x}_{0}\rightarrow\mathbf{x}_{1})*M(\mathbf{x}_{0},\mathbf{x}_{1})\cdot\rho(\mathbf{x}_{0}\rightarrow\mathbf{x}_{1}\rightarrow\mathbf{x}_{2})*M(\mathbf{x}_{1},\mathbf{x}_{2})\cdots}{p(\mathbf{x}_{0}\rightarrow\mathbf{x}_{1})\cdot
p(\mathbf{x}_{0},\mathbf{x}_{1})\cdot
p(\mathbf{x}_{0}\rightarrow\mathbf{x}_{1}\rightarrow\mathbf{x}_{2})\cdot
p(\mathbf{x}_{1},\mathbf{x}_{2})\cdots}$

where $\mathbf{x}_{i}$ is the $i$-th path node and $p$ is the probability for
each path-building event. In the following sections, we will discuss the
details of each term in Eq. 12.

Figure 1. For a path propagating through the scene, the path segments
correspond to the media term of the transport equation, and the vertices
correspond to the BSDF. The events of path propagation and object interaction
are independent.

###  3.4. Wedge Coordinate

As mentioned above, the BSDF
$\rho(\mathbf{x}^{\prime\prime}\rightarrow\mathbf{x}^{\prime}\rightarrow\mathbf{x})$
represents the interaction between wave and scene geometry. Here we’ll write
this function into a simpler form that is independent from the local
coordinate frame:
$\rho(\mathbf{x}^{\prime\prime}\rightarrow\mathbf{x}^{\prime}\rightarrow\mathbf{x})=\rho_{\mathbf{x}}(\mathbf{v}_{i},\mathbf{v}_{o})$.
$\mathbf{v}_{i}=\mathbf{T}^{-1}\frac{\mathbf{x}^{\prime}-\mathbf{x}}{\left||\mathbf{x}^{\prime}-\mathbf{x}\right||_{2}}$,
$\mathbf{v}_{o}=\mathbf{T}^{-1}\frac{\mathbf{x}^{\prime\prime}-\mathbf{x}^{\prime}}{\left||\mathbf{x}^{\prime\prime}-\mathbf{x}^{\prime}\right||_{2}}$.
$\mathbf{T}$ denotes the local coordinate frame we used at position
$\mathbf{x}$. We can see that $\mathbf{v}_{i}$ and $\mathbf{v}_{o}$ represents
the normalized incident and outgoing wave direction vector.

Figure 2. Symbols of the local wedge coordinate system used in this paper. The
ranges for the angles are $(\phi-\pi,\pi-\phi)$ for $\theta_{i}$ and
$\theta_{o}$, $(-\pi/2,\pi/2)$ for $\varphi_{i}$ and $\varphi_{o}$, and
$(0,\pi/2)$ for $\phi$. The dashed vectors are on the normal plane determined
by $\mathbf{n}$ and $\mathbf{b}$.

On the surface of the geometry, we only need to take wave reflection and
refraction effect into consideration. Since the scene geometry is described
with triangle meshes, the surface is locally flat almost everywhere, and the
only natural direction given by the geometry is the normal direction. However,
for edges of the triangle mesh, we can derive a natural local coordinate frame
from the geometry.

The coordinate system used in this paper is illustrated in Fig. 2. Similar to
previous local frames, we name the three vectors in the local coordinate frame
$\mathbf{T}$ normal vector $\mathbf{n}$, tangent vector $\mathbf{t}$ and
bitangent vector $\mathbf{b}$. The expressions of incident and outgoing
vectors in the local coordinate frame are given below:

(13) $\displaystyle\mathbf{v}_{i}$ $\displaystyle=$
$\displaystyle-\cos\theta_{i}\cos\varphi_{i}\mathbf{n}-\sin\theta_{i}\cos\varphi_{i}\mathbf{b}-\sin\varphi_{i}\mathbf{t},$
$\displaystyle\mathbf{v}_{o}$ $\displaystyle=$
$\displaystyle\cos\theta_{o}\cos\varphi_{o}\mathbf{n}+\sin\theta_{o}\cos\varphi_{o}\mathbf{b}+\sin\varphi_{o}\mathbf{t}.$

These expressions allow us to express BEDRF with angles:

$\rho_{\mathbf{x}}(\varphi_{i},\theta_{i},\varphi_{o},\theta_{o})$. The
reciprocity of the wave equation’s solution (Case, 1993) demands that the BSDF
satisfies the constraint
$\rho_{\mathbf{x}}(\mathbf{v}_{i},\mathbf{v}_{o})=\rho_{\mathbf{x}}(-\mathbf{v}_{o},-\mathbf{v}_{i})$.
Replacing the vector parameters with angles, the constraint becomes
$\rho_{\mathbf{x}}(\varphi_{i},\theta_{i},\varphi_{o},\theta_{o})=\rho_{\mathbf{x}}(\varphi_{o},\theta_{o},\varphi_{i},\theta_{i})$.

##  4\. bidirectional edge diffraction response function (BEDRF)

Our new diffraction representation is derived from (Keller and Blank, 1951),
which give the analytic solution of a plane wave reflected and diffracted by a
wedge consisting of two half-planes with Dirichlet or Neumann boundary
condition.

We introduce the following symbols:

(14)
$r=\frac{1-\sin\varphi_{i}\sin\varphi_{o}-|\sin\varphi_{i}-\sin\varphi_{o}|}{\cos\varphi_{i}\cos\varphi_{o}}$

(15) $\kappa=\frac{\pi}{2\pi-2\phi}$

(16)
$f(\theta_{0},\theta_{1},\theta,r)=\frac{1}{\pi}\arctan\left(\frac{\sinh(\kappa
r)\sin\left(\kappa\frac{\theta_{0}-\theta_{1}}{2}\right)}{\cos\kappa(\theta-\frac{\theta_{0}+\theta_{1}}{2})-\cosh(\kappa
r)\cos\left(\kappa\frac{\theta_{0}-\theta_{1}}{2}\right)}\right)$

(17) $\begin{array}[]{rcl}\omega_{lu}&=&\theta_{i}-\pi\\\
\omega_{lr}&=&-\theta_{i}-\pi+2\phi\\\ \omega_{ru}&=&\theta_{i}+\pi\\\
\omega_{rr}&=&-\theta_{i}+\pi-2\phi\\\ \end{array}$

For the Dirichlet boundary condition, the BEDRF is given below:

(18)
$\rho=\frac{1}{2}\left(\frac{df}{dr}(\omega_{lu},\omega_{ru},\theta_{o},\ln
r)+\frac{df}{dr}(\omega_{lr},\omega_{rr},\theta_{o},\ln r)\right)$

For the Neumann boundary condition, the BEDRF becomes

(19)
$\rho=\frac{1}{2}\left(\frac{df}{dr}(\omega_{lu},\omega_{ru},\theta_{o},\ln
r)-\frac{df}{dr}(\omega_{lr},\omega_{rr},\theta_{o},\ln r)\right)$

The expression of $\frac{df}{dr}$ is given in the supplementary material,
along with the detailed derivation process and the proof for reciprocity.

Figure 3. Visualization of BEDRF (left) given by Eq. 18 and the sampler PDF
(right) for the Dirichlet boundary condition, with $\phi=\pi/4$,
$\varphi_{i}=0$ and $\theta_{i}=\pi/2$. The figure shows that the sampler PDF
concentrates on the intensive part of BEDRF.

Fig. 3 visualizes our BEDRF in a special case. Unlike traditional BSDF, our
function is not non-negative. The function also has two singularities on the
incident, outgoing and/or reflection direction (the natural outgoing direction
may not exist in the domain), and the value of BEDRF varies a lot around these
singularities. The variation of this BSDF requires us to design a specialized
sampler for high quality sampling.

###  4.1. Importance Sampling

In Monte Carlo integration, a good sampler can significantly lower the result
variance and reduce the numerical error. When we design an importance sampler
for function $f$, we are looking for a random variable with a probability
density function (PDF) similar to $f$. In the ideal case, the PDF of the
sampler should be $kf$, where $k$ is an constant. Such a sampler can be
constructed for an explicitly given and integrable
$f:\mathbb{R}\mapsto\mathbb{R}$ using the inverse transform sampling scheme
(Ross, 1987): Let
$F(x)=\int_{-\infty}^{x}f(y)dy/\int_{-\infty}^{+\infty}f(y)dy$, $S$ be a
random variable that distributes uniformly on $[0,1]$, then $F^{-1}(S)$ has
the cumulative distribution function (CDF) $F$, and PDF $dF/dx=kf$, where
$k=1/\int_{-\infty}^{+\infty}f(x)dx$.

The PDF of a random variable is always non-negative. However, unlike normal
BSDFs, our BSDF has both positive and negative parts, which rules out the
possibility of finding an ideal sampler. Nevertheless, we can design a
function that is ”similar” to $F^{-1}$ by using a few of its features to
construct our sampler.

Our sampling process proceeds in the following manner. With the predetermined
incident ray angle $\theta_{i}$ and $\varphi_{i}$, we sample $\theta_{o}$
first and $\varphi_{o}$ second. The PDF of our sampler can be written as

(20)
$p(\theta_{o},\varphi_{o})=\frac{p_{\theta}(\theta_{o})p_{\varphi}(\theta_{o},\varphi_{o})}{\cos\varphi_{o}},$

where $p_{\theta}$ and $p_{\varphi}$ are PDFs for $\theta_{o}$ and
$\varphi_{o}$ satisfying the following equations:

(21) $\int_{\phi-\pi}^{\pi-\phi}p_{\theta}(\theta_{o})d\theta_{o}=1$

(22) $\int_{-\pi/2}^{\pi/2}p_{\varphi}(\theta_{o},\varphi_{o})d\varphi_{o}=1$

Under these conditions, we have $\int
p(\theta_{o},\varphi_{o})\cos\varphi_{o}d\theta_{o}d\varphi_{o}=1$. For
$\theta_{o}$, we simply use function
$p_{\theta}(\theta_{o})=\frac{1}{2(\pi-\phi)}$ as our sampler. For
$\varphi_{o}$, we use the following function:

(23) $p_{\varphi}=a+(1-a)\frac{\kappa
r^{\kappa}}{2\cos\varphi_{o}}\frac{dP_{c}(\theta_{a},\theta_{b})}{dx}(r^{k}),a\in[0,1].$

The expression of $a$, $P_{c}$, $\theta_{a}$ and $\theta_{b}$ are complex and
one can find their definitions in the supplementary material of our paper. The
shape of the resultant PDF is very similar to BEDRF and can neutralize its
singularities during sampling.

##  5\. Geometric Propagation using Path Tracing

We use the bidirectional path tracing (BDPT) algorithm from (Cao et al., 2016)
as the backbone of our new integrator. To integrate the diffraction phenomenon
into path tracing, we need to perform integration on edges of the scene
objects. However, in traditional path tracing, the probability of hitting a
point on any edge is exactly zero. Thus we need to make several modifications
to the original BDPT algorithm while ensuring that the path tracing method is
unbiased.

###  5.1. Meta-Path Space

To intersect a geometry with zero hit probability, a natural solution is to
intersect the ray with a “proxy geometry” of nonzero hit probability, and map
the intersection to a point on the original geometry. If we denote the
intersection on the proxy geometry with $\mathbf{x}^{\prime}$, and the mapping
function with $f_{m}$, then the pair
$(\mathbf{x}^{\prime},\mathbf{x}=f_{m}(\mathbf{x}^{\prime}))$ can be regarded
as an generalized form of intersection.

We call the series

(24)
$\mathbf{x}^{*}_{1}\to\mathbf{x}^{*}_{2}\to\cdots\to\mathbf{x}^{*}_{n},\mathbf{x}^{*}_{i}\in\Omega^{*}$

a meta-path on the geometry $\Omega$ if $\Omega^{*}$ is a measure space and
there is a projection $\Pi:\Omega^{*}\to\Omega$ that preserves measurability.
All possible meta-paths containing $n$ nodes form a measure space
$\Omega^{*n}$, called the $n$-th meta-path space.

An example of the meta-path space is the “single-point augmented path space”
described above, where every node of the meta-path is a point pair
$(\mathbf{x}^{\prime},\mathbf{x})$. The first point is on the proxy geometry,
and the second point is on the real one. The projection function for such
space is $\Pi((\mathbf{x}^{\prime},\mathbf{x}))=\mathbf{x}$. There could be
other forms of meta-path nodes that may contain more points or other types of
information.

In traditional path tracing, the intersection probability gives a probability
measure $\mu$ on the path space, allowing us to do the Monte Carlo
integration. For non-intersectable geometries like infinitely thin edges, this
measure no longer exists. However, since the proxy geometry is intersectable,
we can find a similar probability measure $\mu^{*}$ defined on some certain
meta-path space. Now, it can be proved that the canonical projection $\Pi$
induces a probability measure $\Pi(\mu^{*})$ on the original path space. In
addition, every integral on the original path space can be converted to an
integral on $\mu^{*}$. Thus the Monte Carlo integral can be done on the meta-
path space. All the sampling techniques, including MIS and metropolis
sampling, are also feasible on $\mu^{*}$. More information about measure
projection is given in (Tarantola, 2008) . We use this theoretical basis to
work out the necessary details of our method.

Figure 4. An illustration of our edge detection algorithm. We emit a ray from
the source position $\mathbf{x}_{s}$ and intersect with the scene geometry at
$\mathbf{x}^{\prime}\in\triangle\mathbf{abc}_{1}$. Suppose that $\mathbf{ab}$
is a convex edge. We extend $\mathbf{c}_{1}\mathbf{x}^{\prime}$ and intersect
$\mathbf{ab}$ at $\mathbf{x}$, which we take as the real intersection point.
After intersection, we will also calculate the proxy set $S_{\mathbf{x}}$ of
$\mathbf{x}$, which is composed of segments on $\mathbf{xc}_{1}$ and
$\mathbf{xc}_{2}$ visible from $\mathbf{x}_{s}$ (blue lines).

###  5.2. Integration on Edges

It is not difficult to realize that, since we are using triangle meshes for
scene geometry, the geometry itself can be used as a natural proxy geometry
for all its convex edges. Compared with other possible choices for proxy
meshes, using the original geometry requires almost no precalculation, and
allows our algorithm to support dynamic, deformable objects with no extra
overhead.

When detecting edges, we first emit a random ray from the source
$\mathbf{x}_{s}$ and perform an intersection test with the scene geometry $S$.
We name the intersection point $\mathbf{x}^{\prime}\in\triangle\mathbf{abc}$
as the “pseudo-intersection point”. If $\triangle\mathbf{abc}$ has convex
edges, we will select one of them and map $\mathbf{x}^{\prime}$ onto it. The
selection probability of each convex edge is equal and we note it with
$p_{\triangle\mathbf{abc}}$. After selection, we connect $\mathbf{x}^{\prime}$
with the opposite vertex of the selected edge, extend the segment, and
intersect the selected edge at $\mathbf{x}$, which we called the “real
intersection point”. $\mathbf{x}$ and $\mathbf{x}^{\prime}$ form a node
$(\mathbf{x}^{\prime},\mathbf{x})$ in our meta path space. Note that
$\mathbf{x}$ is not necessarily visible from $\mathbf{x}_{s}$. If $\mathbf{x}$
is invisible, then we consider this intersection as a failure.

Now, we need to convert the integration on the edge to an integration on its
neighbor triangles. Take the case in Fig. 4 as an example. The conversion can
be expressed like this:

(25)
$\int_{x\in\mathbf{ab}}fd\mu=\int_{x\in\triangle\mathbf{abc}_{1}\cup\triangle\mathbf{abc}_{2}}f\frac{d\Pi^{-1}(\mu)}{d\mu^{*}}d\mu^{*}.$

Note that one convex edge has two neighbor triangles, and we need to consider
them both. We first calculate all possible pseudo-intersections on the
neighboring triangles that possibly project to $\mathbf{x}$. The set of these
pseudo-intersections is called the “proxy set” of $\mathbf{x}$, noted as
$S_{\mathbf{x}}$. We intersect
$\triangle\mathbf{x}_{s}\mathbf{x}\mathbf{c}_{1}$ and
$\triangle\mathbf{x}_{s}\mathbf{x}\mathbf{c}_{2}$ with $S$, and
$S_{\mathbf{x}}$ is the unoccluded parts of the triangles intersecting with
$S$.

Afterwards, we map the set $S_{\mathbf{x}}$ to $[0,1]$ in the following
manner:

(26) $\displaystyle q(\mathbf{x^{\prime\prime}})$
$\displaystyle=\frac{||\mathbf{x^{\prime\prime}}-\mathbf{c}_{1}||_{2}}{||\mathbf{x}-\mathbf{c}_{1}||_{2}},\mathbf{x}^{\prime\prime}\in
S_{\mathbf{x}}\cap\mathbf{xc}_{1},$ $\displaystyle
q(\mathbf{x^{\prime\prime}})$
$\displaystyle=\frac{||\mathbf{x^{\prime\prime}}-\mathbf{c}_{2}||_{2}}{||\mathbf{x}-\mathbf{c}_{2}||_{2}},\mathbf{x}^{\prime\prime}\in
S_{\mathbf{x}}\cap\mathbf{xc}_{2}.$

Then the expression of $d\Pi^{-1}(\mu)/d\mu^{*}$ can be given as

(27) $\displaystyle\frac{d\mu^{*}}{d\Pi^{-1}(\mu)}(\mathbf{x}^{\prime})$
$\displaystyle=\frac{p_{\triangle\mathbf{abc}_{1}}S\triangle\mathbf{abc}_{1}}{||\mathbf{b}-\mathbf{a}||_{2}}\int_{q(S_{\mathbf{x}}\cap\mathbf{xc}_{1})}xdx,$
$\displaystyle+\frac{p_{\triangle\mathbf{abc}_{2}}S\triangle\mathbf{abc}_{2}}{||\mathbf{b}-\mathbf{a}||_{2}}\int_{q(S_{\mathbf{x}}\cap\mathbf{xc}_{2})}xdx,$

where $S\triangle\mathbf{abc}$ is the area of $\triangle\mathbf{abc}$ and $dx$
is the ordinary differential on $[0,1]$. The proof of this equation is given
in the supplementary material.

###  5.3. MIS Issue for Meta-Paths

In traditional bidirectional path tracing, it is possible to generate the same
path from the source to the listener from multiple sample strategies. For
example, a path with 4 nodes could be generated by strategy $(0,3)$, $(1,2)$,
$(2,1)$ or $(3,0)$, where the two numbers are forward and backward path
depths. The generation probability of a path for different strategies allows
us to use the MIS technique. However, there is some additional complexity for
our meta-paths.

Consider a meta path generated with the algorithm in the previous subsection:

(28)
$(\mathbf{x}^{\prime}_{0},\mathbf{x}_{0})\to(\mathbf{x}^{\prime}_{1},\mathbf{x}_{1})\to\cdots\to(\mathbf{x}^{\prime}_{n},\mathbf{x}_{n})$

If $V$ is the visibility function, we know from the generation process that
$V(\mathbf{x}_{i},\mathbf{x}_{i+1})=1$,
$V(\mathbf{x}_{i},\mathbf{x}^{\prime}_{i+1})=1$. Now, MIS requires that the
inverse of the path

(29)
$(\mathbf{x}^{\prime}_{0},\mathbf{x}_{0})\leftarrow\cdots\leftarrow(\mathbf{x}^{\prime}_{n-1},\mathbf{x}_{n-1})\leftarrow(\mathbf{x}^{\prime}_{n},\mathbf{x}_{n})$

can be generated by the backward path tracer and its generation probability
computable. However, this requires an additional condition
$V(\mathbf{x}^{\prime}_{i},\mathbf{x}_{i+1})=1$. Since our edge tracing
algorithm is not symmetric, this is not guaranteed by the generation process.

To solve this problem, we use a twice-augmented path space in our renderer,
where the path is expressed as

(30)
$(\mathbf{x}^{\prime\prime}_{0},\mathbf{x}^{\prime}_{0},\mathbf{x}_{0})\to(\mathbf{x}^{\prime\prime}_{1},\mathbf{x}^{\prime}_{1},\mathbf{x}_{1})\to\cdots\to(\mathbf{x}^{\prime\prime}_{n},\mathbf{x}^{\prime}_{n},\mathbf{x}_{n}).$

Both $\mathbf{x}^{\prime}_{i}$ and $\mathbf{x}^{\prime\prime}_{i}$ projects to
$\mathbf{x}$. We require that the path satisfies one more condition
$V(\mathbf{x}^{\prime\prime}_{i},\mathbf{x}_{i+1})=1$. In the generation
process, this is achieved by choosing a random $\mathbf{x}^{\prime\prime}_{i}$
from the proxy segment visible from $\mathbf{x}_{i+1}$ (in simpler words,
“reverse tracing”). In this way, the three visible conditions are symmetric
under the inverse transform. We can then use MIS without any problems.

###  5.4. Inverse Russian Roulette

We have seen from Sect. 5.2 that our edge intersection algorithm may fail when
the real intersection point is occluded. This is not an uncommon case in
actual application. Hence it is harder to generate path segments for
diffraction. To exploit every valid intersection at maximum, we use a
technique similar to the Russian roulette in (Veach, 1997), but reversing the
process: instead of discarding low-intensity path nodes, we divide high-
intensity nodes into multiple shares to generate more path segments.

After each bounce of path tracing, we will duplicate successful intersections
and replace the failed intersections sharing the same starting node. When an
intersection is duplicated $n$ times, we will multiply the generation
probability of the original intersection and all its duplicates by $n+1$. If
there are multiple successful intersections, the number of duplications of
each intersection is largely proportional to its intensity, which is $f/p$ in
Eq. 12. In this way, the intensity of duplicated intersections would be close
to each other. All the intersections generated in this way forms a forward and
backward tree, where intensive nodes have more successor branches than other
nodes. In this way, we are able to generate much more valid path nodes without
significantly increasing the number of intersection tests.

Figure 5. Illustration of inverse Russian roulette. We duplicate successful
intersections and replace failed intersections with the same starting node
with duplicated intersections. In Fig.(̇a), we use dashed lines and solid
arrows to represent failed/successful intersection, and vertical black lines
to represent node duplication. When connecting forward and backward paths, we
are actually connecting nodes between a forward and a backward tree, which is
shown in Fig.(̇b).

###  5.5. Outlier Suppression

We know from Sect. 5.2 that multiple pseudo-intersections can be mapped to the
same real intersection point. Different pseudo-intersection positions lead to
different sample generating probabilities, and occasionally it could be very
low. When the direction of some certain path segment is also near the
singularities of the BEDRF, the sample intensity would be very high and
produces an outlier in the result. The position of the pseudo-intersection is
determined by the outgoing direction of the previous BEDRF sampler and is
difficult to control. But it is not hard to identify these samples and
suppress their intensity.

Here we note the starting position of a ray segment as $\mathbf{x}_{s}$, the
pseudo-intersection as $\mathbf{x}^{\prime}$, the incident angle at the
pseudo-intersection as $\theta$, the real intersection as $\mathbf{x}$, and
the outgoing probability of the sampler at $\mathbf{x}_{s}$ on the direction
$\mathbf{v}$ as $p(\mathbf{v})$. It is not hard to calculate

(31)
$C_{o}=\frac{p(\mathbf{x}-\mathbf{x}_{s})}{p(\mathbf{x}^{\prime}-\mathbf{x}_{s})}\cdot\frac{||\mathbf{x}^{\prime}-\mathbf{x}_{s}||_{2}^{2}}{||\mathbf{x}-\mathbf{x}_{s}||_{2}^{2}\cos\theta}.$

A large $C_{o}$ value indicates that the sample probability on the direction
$\mathbf{x}^{\prime}-\mathbf{x}_{s}$ is much smaller than that on the
direction $\mathbf{x}-\mathbf{x}_{s}$, thus will likely generate an outlier.
The $C_{o}$ value of a sample (a full path) is the product of all $C_{o}$
values of its non-connection path segments.

To suppress the outliers, we allow the user to set a maximum threshold for
$C_{o}$, which we note as $C^{*}_{o}$ here. We then multiply the intensity of
every sample with a logistic suppression coefficient $C_{s}$. The expression
of $C_{s}$ is given below:

(32)
$C_{s}=\frac{C_{o}^{*}(1-e^{-2C_{o}/C_{o}^{*}})}{C_{o}(1+e^{-2C_{o}/C_{o}^{*}})}$

Notice that when we use outlier suppression, the result of our path tracer is
no longer unbiased. However, the suppression coefficient will only affect
samples with high $C_{o}$ and low probability, and it will not affect the
intensity of the simulation result too much when the scene geometry is simple
and the sample quality is high. This is also shown in Fig. 6 in the result
section.

##  6\. Implementation and Results

In this section, we describe our implementation and highlight its performance
on different benchmarks. Some of the results are presented in the
supplementary video of this paper.

###  6.1. Implementation and Performance

model | Boxes | Roomset | Sponza
---|---|---|---
|  |  |
vertices | 20 | 11087 | 153635
triangles | 26 | 22974 | 279133
diffraction edges | 16 | 17184 | 220056
other features |  | complex occlusion |
time cost per frame | 34ms | 140ms | 205ms
Table 2. Information of benchmark models used for auralization and performance
of our diffraction path tracer in these models.

Our path tracer runs on a commodity PC with an Intel i7 2.60Ghz CPU and 12GB
memory, and executes on a single thread. To evaluate the result quality under
a certain render configuration, we compute 1000 impulse responses (IRs) with
12,000 samples. After this, we calculate the expectation of all the result
IRs, and consider the difference between each calculated IR and the
expectation IR as its error. We evaluate the render quality by computing the
signal-to-noise ratio (SNR) for IR on its spectral domain. Given an IR
function $f$ and its numeric error $f_{e}$, its IR SNR at frequency $\omega$
is calculated in the following way:

(33)
$\text{SNR}(f)(\omega)=10\log_{10}\frac{|\mathcal{F}(f)(\omega)|}{|\mathcal{F}(f_{e})(\omega)|}$

Where $\mathcal{F}(f)$ is the Fourier spectrum of $f$. In Fig. 6, the IR SNR
given is the average SNR of all computed IRs in the audible frequency range
(20–20,000 Hz). When we use outlier suppression during quality evaluation, we
set the threshold $C^{*}_{o}$ in Eq. 32 to 100.

To demonstrate the influence and computation efficiency of various path
tracing techniques, we calculate the diffraction result in different scenes
with our path tracer using different render configurations. We test our path
tracer with the box scene in Fig. 7 and a more complex scene called “Roomset”,
which is an indoor scene with several small rooms connected with narrow
corridors. The complex occlusion relationships between faces and edges in this
scene is challenging for interactive geometric acoustic methods.

Boxes | Roomset
---|---
Plain BDPT
IR SNR: -3.1618 dB | time cost: 19ms | IR SNR: -4.5834 dB | time cost: 60 ms
|
BDPT+MIS
IR SNR: -2.4164 dB | time cost: 22ms | IR SNR: -4.4580 dB | time cost: 73 ms
|
BDPT+MIS+inverse Russian roulette
IR SNR: 0.0630 dB | time cost: 34ms | IR SNR: -7.6249 dB | time cost: 138 ms
|
BDPT+MIS+inverse Russian roulette, outlier suppressed
IR SNR: 0.1271 dB | time cost: 34ms | IR SNR: -0.5212 dB | time cost: 140 ms
|
Figure 6. IRs simulated in different scenes with different rendering technique
combinations. We highlight the performance on two challenging benchmarks for
geometric acoustics in terms of diffraction effects. The average IR SNR and
the computation time cost of a single IR is given for each case. We show the
relative benefits of MIS, inverse Russian roulette, and outlier suppressed for
each of these complex models.

Fig. 6 shows the IR and its standard variance produced by our path tracer
under different rendering configurations, along with the average IR SNR and
the average time cost of IR calculation. We notice in Fig. 6 that when we use
more techniques in Sect. 5 to improve our approach (e.g., MIS, inverse Russian
Roulette), the accuracy of the resulting IR improves. The only exception is
the inverse Russian roulette (IRR) in the Roomset scene. However, if we
compare the variance curve with the one without IRR, we can see that IRR
actually lowers the variance curve in most places, and the IR SNR decrease is
caused by a few outliers. After outlier suppression, the IR SNR improves
considerably. We can also see from Fig. 6 that MIS and outlier suppression
have little impact on the simulation efficiency, while IRR almost doubled the
simulation time cost. Considering that doubling the path tracing samples will
only increase the IR SNR by approximately 1.5 dB, IRR clearly has an advantage
in terms of improving the accuracy.

We also demonstrated the auralization result of our path tracer in our video.
We used the aforementioned box scene, Roomset and a more complex model “Crytek
Sponza” with more than 14k triangles for demonstration. We use 36k reflection
samples and 12k diffraction samples (if diffraction is present) for audio
generation in each frame. In the video, we put a single source in the scene
and move the audio receiver around. During the motion, the source is often
totally occluded from the receiver by the scene geometry, and higher order
diffraction plays a major role in terms of generating smooth auralization
results. The time cost of our diffraction path tracer is between 30–205 ms per
frame. Information about auralization models and results can be found in Table
2. Further details are presented in the video.

###  6.2. Comparisons

We have tested our diffraction BSDF and integrator in various scene
configurations and compared the results with those of state-of-the-art
diffraction models, notably the BTM and UTD models. The BTM model gives the
exact solution for single-edge diffraction for spherical incident waves, and
our result matches the BTM solution in both the single-edge and the double-
edge cases (see Fig. 7).

Our time-domain algorithm shows a clear advantage over frequency-domain
algorithms like UTD. Specifically, while UTD produces near-perfect results
when there are very high number of samples on the frequency domain, it is much
less accurate in interactive propagation algorithms (Schissler et al., 2014),
where only 4-8 frequency bands are used and the phase information of the
result is discarded. In Fig. 7, the “8-band UTD” IR is produced by calculating
the UTD result on the central frequency of 8 bands ranging from 62.5-8000Hz,
and combining them with crossover filters. The result shows the typical
“ringing” effect caused by undersampling in the frequency domain. The
performance of the resulting propagation algorithm increases linearly with the
number of bands.

For more complex scenes where an analytic solution is not known, We compare
with the result of Svensson’s diffraction toolbox for Matlab, which uses the
BTM algorithm for edge diffraction (Svensson, 2015; Torres et al., 2001). The
scene we used for comparison consists of a floor plane and two boxes of size
$\text{12m}\times\text{12m}\times\text{4m}$. We set the highest diffraction
order to 2 for our implementation and Svensson’s toolbox. Our result matches
Svensson’s result very well at many diffraction peaks, but there are also some
other peaks in Svensson’s result not present in ours. This is due to the
presence of “creeping waves”, which is the sound wave that doesn’t bounce
between faces or edges, but travels along the surface of the geometry.
Currently, our algorithm is unable to simulate this phenomenon.

|
---|---
|
|
Figure 7. IRs of diffraction wave in several simple scenes. The diffraction
IRs are generated by the corresponding scene illustrated on the left side. In
all test scenes, we set sound speed to 344m/s, source impulse energy to 1J and
air density to 1.21$\text{kg}/\text{m}^{3}$. We use the Dirichlet boundary
condition condition for the first two scenes and Neumann condition for the
last scene. The difference between our method and BTM is minimal in the first
two cases. Moreover, major diffraction peaks are well-matched in the last one.
These results indicate close agreement with BTM on these scenarios. Meanwhile,
we also observe improved accuracy of our BEDRF formulation over UTD in the top
example.

Since our algorithm is amplitude-based, it is also capable of simulating
interference of diffraction waves produced by two or more edges. Fig. 8 shows
the result of our algorithm reproducing the single-slit experiment, which
involves a plane wave passing through a narrow slit of constant width and
infinite length, producing a typical stripe-like pattern. Our result closely
matches the precomputed ground truth.

Figure 8. Stripe pattern produced by single-slit diffraction. The result is
calculated using 50000 samples and closely matches the ground truth.

##  7\. Conclusion, Limitations and Future Work

In this paper, we have presented BEDRF, a new localized representation of edge
diffraction for convex wedges. We also present a path tracer utilizing the
BEDRF representation for interactive simulation of sound diffraction in
complex models. The flexibility of path tracing allows us to deal with complex
geometries with ease, which has been shown in our results. The BEDRF model and
the new path tracer are mutually independent. We can replace BEDRF with BTM or
DLSM without too much work. The latter two models may have advantages in some
special cases. BEDRF can perform interactive sound propagation on complex
models with no preprocessing or model simplification, though its accuracy may
not match with BTM in some cases.

The computation cost of “diffraction path tracing” is still much higher than
traditional path tracing. The formation of BEDRF itself is already more
complex than that of common BSDFs for surface reflection. In addition, we use
triangle-triangle intersection in our tracing algorithm, which is much more
costly than ray-triangle intersection. It seems to be challenging to design a
purely ray-based unbiased integrator on geometric edges. However, when we
remove the restriction of unbiasedness, we may find an algorithm that relies
on ray-triangle intersections only. That is a good topic for future work.

We have evaluated the accuracy on many benchmarks. We observe that the
accuracy of our algorithm is nearly-indistinguishable from that of BTM-based
solvers in simple scenes. As part of future work, it would be useful to
investigate the accuracy of BEDRF with respect to BTM. At the moment, it is
unknown whether BEDRF is exact when the incident wave is non-planar. It is
also worth noting that our path tracer doesn’t support the simulation of
creeping waves. This phenomenon is especially important for simulating the
diffraction effect of convex smooth objects, like spheres and cylinders. Our
future simulator will take this phenomenon into consideration.

## References

  * (1)
  * Allen and Raghuvanshi (2015) Andrew Allen and Nikunj Raghuvanshi. 2015.  Aerophones in flatland.  _ACM Transactions on Graphics (TOG)_ 34 (2015), 1 – 11.
  * Antani et al. (2010) Lakulish Antani, Anish Chandak, Micah Taylor, and Dinesh Manocha. 2010.  Fast Geometric Sound Propagation with Finite Edge Diffraction.  _Technical Report TR10-011, University of North Carolina at Chapel Hill_ (01 2010).
  * Asmail (1991) Clara C. Asmail. 1991\.  Bidirectional Scattering Distribution Function (BSDF): A Systematized Bibliography.  _Journal of Research of the National Institute of Standards and Technology_ 96 (1991), 215 – 223.
  * Bertram et al. (2005) M. Bertram, E. Deines, J. Mohring, J. Jegorovs, and H. Hagen. 2005.  Phonon tracing for auralization and visualization of sound. In _VIS 05. IEEE Visualization, 2005._ 151–158.  https://doi.org/10.1109/VISUAL.2005.1532790
  * Biot and Tolstoy (1957) M A Biot and I Tolstoy. 1957.  Formulation of Wave Propagation in Infinite Media by Normal Coordinates with an Application to Diffraction.  _Journal of the Acoustical Society of America_ 29, 3 (1957), 381–391.
  * Calamia (2009) Paul Thomas Calamia. 2009\.  _Advances in Edge-Diffraction Modeling for Virtual-Acoustic Simulations_.  Ph. D. Dissertation. USA.  Advisor(s) Funkhouser, Thomas.
  * Calamia and Svensson (2007) Paul T. Calamia and U. Peter Svensson. 2007.  Fast Time-Domain Edge-Diffraction Calculations for Interactive Acoustic Simulations.  _EURASIP J. Adv. Signal Process_ 2007, 1 (Jan. 2007), 186.  https://doi.org/10.1155/2007/63560
  * Cao et al. (2016) Chunxiao Cao, Zhong Ren, Carl Schissler, Dinesh Manocha, and Kun Zhou. 2016.  Interactive sound propagation with bidirectional path tracing.  _ACM Transactions on Graphics (TOG)_ 35 (2016), 1 – 11.
  * Case (1993) K. Case. 1993.  Structural Acoustics: A General Form of Reciprocity Principles in Acoustics.  _Mitre Corp. Technical Report JSR-91-193_.
  * Hamilton (2021) Brian Hamilton. 2021\.  Tutorial on finite-difference time-domain (FDTD) methods for room acoustics simulation.  _The Journal of the Acoustical Society of America_ 149, 4 (2021), A92–A93.  https://doi.org/10.1121/10.0004614 arXiv:https://doi.org/10.1121/10.0004614
  * Jensen (2001) Henrik Wann Jensen. 2001\.  Realistic Image Synthesis Using Photon Mapping.
  * Kajiya (1986) James T Kajiya. 1986\.  The rendering equation.  _Proceedings of the 13th Annual Conference on Computer Graphics and Interactive Techniques_ 20, 4 (1986), 143–150.
  * Kapralos et al. (2007) Bill Kapralos, Michael R. M. Jenkin, and Evangelos E. Milios. 2007\.  Acoustical Modeling with Sonel Mapping. In _19th International Congress on Acoustics_.
  * Keller (1997) Alexander Keller. 1997\.  Instant radiosity.  _Proceedings of the 24th Annual Conference on Computer Graphics and Interactive Techniques_ (1997).
  * Keller and Blank (1951) Joseph B Keller and Albert Blank. 1951.  Diffraction and reflection of pulses by wedges and corners.  _Communications on Pure and Applied Mathematics_ 4, 1 (1951), 75–94.
  * Kirkup (2019) Stephen Kirkup. 2019\.  The Boundary Element Method in Acoustics: A Survey.  _Applied Sciences_ 9 (04 2019), 1642.  https://doi.org/10.3390/app9081642
  * Kouyoumjian and Pathak (1974) R. G. Kouyoumjian and P. H. Pathak. 1974.  A uniform geometrical theory of diffraction for an edge in a perfectly conducting surface.  _Proc. IEEE_ 62, 11 (1974), 1448–1461.
  * Kreyszig (1978) E Kreyszig. 1978\.  Introductory Functional Analysis with Applications.  (1978).
  * Lentz et al. (2007) Tobias Lentz, Dirk Schroder, Michael Vorlander, and Ingo Assenmacher. 2007.  Virtual reality system with integrated sound field simulation and reproduction.  _EURASIP Journal on Advances in Signal Processing_ 2007, 1 (2007), 187–187.
  * Liu and Manocha (2020) Shiguang Liu and Dinesh Manocha. 2020.  Sound Synthesis, Propagation, and Rendering: A Survey.  _CoRR_ abs/2011.05538 (2020).  arXiv:2011.05538 https://arxiv.org/abs/2011.05538
  * Mcnamara et al. (1990) D. Mcnamara, C. Pistorius, and J. Malherbe. 1990\.  _Introduction to The Uniform Geometrical Theory and Diffraction_.
  * Medwin et al. (1982) Herman Medwin, Emily Childs, and Gary M Jebsen. 1982\.  Impulse studies of double diffraction: A discrete Huygens interpretation.  _Journal of the Acoustical Society of America_ 72, 3 (1982), 1005–1013.
  * Mehra et al. (2013) Ravish Mehra, Nikunj Raghuvanshi, Lakulish Antani, Anish Chandak, Sean Curtis, and Dinesh Manocha. 2013.  Wave-based sound propagation in large open scenes using an equivalent source formulation.  _ACM Transactions on Graphics_ 32, 2 (2013), 19\.
  * Mehra et al. (2012) Ravish Mehra, Nikunj Raghuvanshi, Lauri Savioja, Ming C. Lin, and Dinesh Manocha. 2012.  An efficient GPU-based time domain solver for the acoustic wave equation.  _Applied Acoustics_ 73 (2012), 83–94.
  * Menounou et al. (2000) Penelope Menounou, Ilene J. Busch-Vishniac, and David T. Blackstock. 2000.  Directive line source model: A new model for sound diffraction by half planes and wedges.  _The Journal of the Acoustical Society of America_ 107, 6 (2000), 2973–2986.  https://doi.org/10.1121/1.429327 arXiv:https://doi.org/10.1121/1.429327
  * Menounou and Nikolaou (2017) Penelope Menounou and Petros Nikolaou. 2017.  Analytical model for predicting edge diffraction in the time domain.  _The Journal of the Acoustical Society of America_ 142, 6 (2017), 3580–3592.  https://doi.org/10.1121/1.5014051 arXiv:https://doi.org/10.1121/1.5014051
  * Pind Jörgensson (2020) Finnur Kári Pind Jörgensson. 2020\.  _Wave-Based Virtual Acoustics_.  Ph. D. Dissertation.
  * Pisha et al. (2020) Louis Pisha, Siddharth Atre, John Burnett, and Shahrokh Yadegari. 2020.  Approximate diffraction modeling for real-time sound propagation simulation.  _The Journal of the Acoustical Society of America_ 148 (10 2020), 1922–1933.  https://doi.org/10.1121/10.0002115
  * Raghuvanshi et al. (2009) Nikunj Raghuvanshi, Rahul Narain, and Ming C. Lin. 2009\.  Efficient and accurate sound propagation using adaptive rectangular decomposition.  _IEEE transactions on visualization and computer graphics_ 15 5 (2009), 789–801.
  * Ross (1987) Sheldon M. Ross. 1987\.  Introduction to Probability and Statistics for Engineers and Scientists.
  * Royden and Fitzpatrick (1988) Halsey Lawrence Royden and Patrick Fitzpatrick. 1988.  _Real analysis_. Vol. 198.  Macmillan New York.
  * Rungta et al. (2016) Atul Rungta, Sarah Rust, Nicolas Morales, Roberta Klatzky, Ming Lin, and Dinesh Manocha. 2016\.  Psychoacoustic Characterization of Propagation Effects in Virtual Environments.  _ACM Trans. Appl. Percept._ 13, 4, Article 21 (jul 2016), 18 pages.  https://doi.org/10.1145/2947508
  * Rungta et al. (2018) Atul Rungta, Carl Schissler, Nicholas Rewkowski, Ravish Mehra, and Dinesh Manocha. 2018\.  Diffraction Kernels for Interactive Sound Propagation in Dynamic Environments.  _IEEE Transactions on Visualization and Computer Graphics_ 24, 4 (2018), 1613–1622.
  * Savioja (2010) Lauri Savioja. 2010\.  Real-time 3D finite-difference time-domain simulation of low-and mid-frequency room acoustics. In _13th Int. Conf on Digital Audio Effects_ , Vol. 1. 75.
  * Schissler and Manocha (2016) Carl Schissler and Dinesh Manocha. 2016.  Interactive Sound Propagation and Rendering for Large Multi-Source Scenes.  _ACM Trans. Graph._ 36, 1, Article 2 (sep 2016), 12 pages.  https://doi.org/10.1145/2943779
  * Schissler et al. (2014) Carl Schissler, Ravish Mehra, and Dinesh Manocha. 2014\.  High-order diffraction and diffuse reflections for interactive sound propagation in large environments.  33, 4 (2014), 39\.
  * Schissler et al. (2021) Carl Schissler, Gregor Mückl, and Paul T. Calamia. 2021\.  Fast diffraction pathfinding for dynamic sound propagation.  _ACM Transactions on Graphics (TOG)_ 40 (2021), 1 – 13.
  * Siltanen et al. (2007) Samuel Siltanen, Tapio Lokki, Sami Kiminki, and Lauri Savioja. 2007.  The room acoustic rendering equation.  _Journal of the Acoustical Society of America_ 122, 3 (2007), 1624–1635.
  * Stephenson (2010) Uwe M. Stephenson. 2010\.  An Energetic Approach for the Simulation of Diffraction within Ray Tracing Based on the Uncertainty Relation.  _Acta Acustica united with Acustica_ 96, 3 (2010), 516–535.  https://doi.org/doi:10.3813/AAA.918304
  * Stephenson and Svensson (2007) Uwe Martin Stephenson and U. Peter Svensson. 2007.  An Improved Energetic Approach to Diffraction Based on the Uncertainty Principle. In _Proceedings of the 19th International Congress on Acoustics_ (Madrid, Spain).
  * Svensson (2015) U. Peter Svensson. 2015\.  Edge diffraction toolbox for Matlab.  Retrieved May 20, 2022 from https://folk.ntnu.no/ulfps/software/index.html#EDGE
  * Svensson et al. (1999) U Peter Svensson, Roger I Fred, and John Vanderkooy. 1999\.  An analytic secondary source model of edge diffraction impulse responses.  _Journal of the Acoustical Society of America_ 106, 5 (1999), 2331–2344.
  * Tang et al. (2021) Zhenyu Tang, Hsien-Yu Meng, and Dinesh Manocha. 2021\.  Learning Acoustic Scattering Fields for Dynamic Interactive Sound Propagation. In _2021 IEEE Virtual Reality and 3D User Interfaces (VR)_. 835–844.  https://doi.org/10.1109/VR50410.2021.00111
  * Tarantola (2008) Albert Tarantola. 2008\.  Image and Reciprocal Image of a Measure. Compatibility Theorem.
  * Taylor et al. (2012) Micah Taylor, Anish Chandak, Qi Mo, Christian Lauterbach, Carl Schissler, and Dinesh Manocha. 2012\.  Guided Multiview Ray Tracing for Fast Auralization.  _IEEE Transactions on Visualization and Computer Graphics_ 18, 11 (2012), 1797–1810.  https://doi.org/10.1109/TVCG.2012.27
  * Thompson (2006) Lonny L. Thompson. 2006\.  A review of finite-element methods for time-harmonic acoustics.  _The Journal of the Acoustical Society of America_ 119, 3 (2006), 1315–1330.  https://doi.org/10.1121/1.2164987 arXiv:https://doi.org/10.1121/1.2164987
  * Torres et al. (2001) Rendell R. Torres, U. Peter Svensson, and Mendel Kleiner. 2001\.  Computation of edge diffraction for more accurate room acoustics auralization.  _The Journal of the Acoustical Society of America_ 109 2 (2001), 600–10.
  * Tsingos et al. (2007) Nicolas Tsingos, Carsten Dachsbacher, Sylvain Lefebvre, and Matteo Dellepiane. 2007.  Instant Sound Scattering. In _Proceedings of the 18th Eurographics Conference on Rendering Techniques_ (Grenoble, France) _(EGSR’07)_. Eurographics Association, Goslar, DEU, 111–120.
  * Tsingos et al. (2001) Nicolas Tsingos, Thomas Funkhouser, Addy Ngan, and Ingrid Carlbom. 2001.  Modeling acoustics in virtual environments using the uniform theory of diffraction.  (2001), 545–552.
  * Veach (1997) Eric Veach. 1997\.  _Robust Monte Carlo methods for light transport simulation_.  Ph. D. Dissertation.
  * Veach and Guibas (1994) Eric Veach and Leonidas J. Guibas. 1994.  Bidirectional Estimators for Light Transport. In _Proceedings of Eurographics Workshop on Rendering’94_.
  * Veach and Guibas (1997) Eric Veach and Leonidas J. Guibas. 1997.  Metropolis light transport. In _Proceedings of the 24th Annual Conference on Computer Graphics and Interactive Techniques_.
  * Vorländer (1989) Michael Vorländer. 1989\.  Simulation of the transient and steady‐state sound propagation in rooms using a new combined ray‐tracing/image‐source algorithm.  _The Journal of the Acoustical Society of America_ 86, 1 (1989), 172–178.  https://doi.org/10.1121/1.398336 arXiv:https://doi.org/10.1121/1.398336
  * Yeh et al. (2013) Hengchin Yeh, Ravish Mehra, Zhimin Ren, Lakulish Antani, Dinesh Manocha, and Ming C Lin. 2013\.  Wave-ray coupling for interactive sound propagation in large complex scenes. In _International Conference on Computer Graphics and Interactive Techniques_ , Vol. 32. 165.