Noise debiasing


Consider the following integral: \[F = \int dx dy f(x, y) \label{eq:main}\] where \(f\) is the number of photon counts in a pixel \((x, y)\) and it follows the probability distribution \(P(f|\mu, \sigma)\) which is a Gaussian centered at \(\mu\), with dispersion \(\sigma\).

Assuming ergodicity, we can reduce this integral to \[F = \int dx dy \int_{-\infty}^{\infty} df f P(f| \mu, \sigma) = A \mu\] where \(A\) is the area integral, which we assume is unity. Consider the case, \[F' = \int dx dy |f|(x, y)\] Then we would end up with \[F' = \int_{-\infty}^{0} df (-f) P(f| \mu, \sigma) + \int_{0}^{\infty} df f P(f| \mu, \sigma)\] \[F' = \int_{0}^{\infty} df f P(f| -\mu, \sigma) + \int_{0}^{\infty} df f P(f| \mu, \sigma)\] \[F' = \int_{0}^{\infty} df f \frac{1}{\sqrt{2\pi}\sigma}\exp\left(-\frac{[f+\mu]^2}{2\sigma^2} \right) + \int_{0}^{\infty} df f \frac{1}{\sqrt{2\pi}\sigma}\exp\left(-\frac{[f-\mu]^2}{2\sigma^2} \right)\] \[F' = I(-\mu, \sigma) + I(\mu, \sigma)\] So now let us compute \[I(\mu, \sigma) = \int_{0}^{\infty} df f \frac{1}{\sqrt{2\pi}\sigma}\exp\left(-\frac{[f-\mu]^2}{2\sigma^2} \right)\] Substituting \(f' = f-\mu\) results in \(df' = df\) and the integral \[I(\mu, \sigma) = \int_{-\mu}^{\infty} df' (f'+\mu) \frac{1}{\sqrt{2\pi}\sigma}\exp\left(-\frac{f'^2}{2\sigma^2} \right)\] which can be split into \[I(\mu, \sigma) = \int_{0}^{\infty} ... + \int_{-\mu}^{0} = I_1 + I_2\] Let us consider \(I_1\) first. \[I_1(\mu, \sigma) = \int_{0}^{\infty} df' (f'+\mu) \frac{1}{\sqrt{2\pi}\sigma}\exp\left(-\frac{f'^2}{2\sigma^2} \right)\] \[I_1(\mu, \sigma) = \frac{\mu}{2} + \int_{0}^{\infty} df' f' \frac{1}{\sqrt{2\pi}\sigma}\exp\left(-\frac{f'^2}{2\sigma^2} \right)\] Substituting \(f'^2/(2\sigma^2)=t\) results in \(f'df'/\sigma = \sigma dt\), and we obtain \[I_1(\mu, \sigma) = \frac{\mu}{2} + \int_{0}^{\infty} \frac{\sigma dt}{\sqrt{2\pi}}\exp(-t) = \frac{\mu}{2} + \frac{\sigma}{\sqrt{2\pi}}\]

Now let us consider I_2,

\[I_2(\mu, \sigma) = \int_{-\mu}^{0} \frac{(f'+\mu)}{\sqrt{2\pi}\sigma}\exp\left(-\frac{f'^2}{2\sigma^2}\right)df'\]

\[I_2(\mu, \sigma) = \int_{-\mu}^{0} \frac{f'}{\sqrt{2\pi}\sigma}\exp\left(-\frac{f'^2}{2\sigma^2}\right)df' + \frac{\mu}{2} \,\,{\rm erf} \left( \frac{\mu}{\sqrt{2}\sigma} \right)\]

\[I_2(\mu, \sigma) = \int_{\mu^2/(2\sigma^2)}^{0} \frac{\sigma dt}{\sqrt{2\pi}}\exp\left(-t\right) + \frac{\mu}{2} \,\,{\rm erf} \left( \frac{\mu}{\sqrt{2}\sigma} \right)\]

\[I_2(\mu, \sigma) = \frac{\sigma}{\sqrt{2\pi}}\left[ \exp\left(-\frac{\mu^2}{2\sigma^2}\right) -1\right] + \frac{\mu}{2} \,\,{\rm erf} \left( \frac{\mu}{\sqrt{2}\sigma} \right)\]

So we finally arrive at the expression for \(I(\mu, \sigma)\), \[I(\mu, \sigma) = \frac{\mu}{2} + \frac{\sigma}{\sqrt{2\pi}} + \frac{\sigma}{\sqrt{2\pi}}\left[ \exp\left(-\frac{\mu^2}{2\sigma^2}\right) -1\right] + \frac{\mu}{2} \,\,{\rm erf} \left( \frac{\mu}{\sqrt{2}\sigma} \right)\] And this implies that \[F' = I(\mu, \sigma) + I(-\mu, \sigma)\] Therefore, we obtain \[F' = \sqrt{\frac{2}{\pi}}\sigma + \sqrt{\frac{2}{\pi}}\sigma\left[ \exp\left(-\frac{\mu^2}{2\sigma^2}\right) -1\right] + \mu \,\,{\rm erf} \left( \frac{\mu}{\sqrt{2}\sigma} \right)\] when we use \(|f|\) instead of \[F = \mu\] when we use \(f\) in Equation \ref{eq:main}.

Let us just quickly consider two extreme cases, one is when \(\mu=0\). In this case, we obtain \[F = 0\] and \[F' = \sqrt{\frac{2}{\pi}}\sigma\] while if \(\mu\) is large compared to \(\sigma\), then we get \[F = F' = \mu\] Thus no debiasing is really required if \(\mu\) is really large.

So to perform noise debiasing in principle we should do \[F' = \int dx dy (|f| - d )(x, y)\] where \[d = \sqrt{\frac{2}{\pi}}\sigma + \sqrt{\frac{2}{\pi}}\sigma\left[ \exp\left(-\frac{\mu^2}{2\sigma^2}\right) -1\right] + \mu \left[{\rm erf} \left( \frac{\mu}{\sqrt{2}\sigma} \right) - 1\right]\]