While I was working on noise profile creation for darktable, I had to learn how variance stabilization transform work.
So, here is a summary of what I learnt about this
What are variance stabilization transform?
Noise variance is not constant: the variance depends on the light intensity.
Typically, variance is higher in bright areas than in dark areas.
Plenty of noise algorithms are designed to work on images that have constant variance.
To use these algorithms, we can transform the data so that variance becomes constant.
Such a transform is called a variance stabilizing transform.
Why and How does it works?
First thing to know is how variance evolves depending on the mean.
Let’s say we have empirically: V(X)=g(E[X])
We can try fitting a model to these empirical data.
For instance, for creating darktable’s noise profiles, we use the gaussian-poissonian model, the variance is assumed to evolve linearly with regards to the mean: V(X)=a*E[X]+b
So, we try to find a and b such that the model fits as much as possible the reality (in practice, in darktable, these a and b are exactly what your noise profiles are).
At this step, the “only” thing we have to do is to find a model that fits well the reality
Then, we can design the transform.
Let f be our variance stabilization transform.
We want V(f(X))=c, and we know V(X)=g(E[X]).
Let’s write m = E[X].
Using Taylor expansion (see Taylor expansions for the moments of functions of random variables - Wikipedia), we get:
V(f(X)) \approx V(f(m)+f'(m)(X-m))
Which we can simplify using the properties of variance:
V(f(X)) \approx V(f'(m)(X-m))
V(f(X)) \approx f'(m)^2*V((X-m))
V(f(X)) \approx f'(m)^2*V(X)
So our condition V(f(X))=c gives us the following condition:
f'(m)^2*V(X)=c
We know that V(X)=g(E[X]), so we get:
f'(m)^2*g(m)=c
For all values of m such as g(m)\neq 0, we have:
f'(m)=\sqrt{\frac{c}{g(m)}}
We usually choose c=1:
f'(m)=\frac{1}{\sqrt{g(m)}}
From this, we can integrate to get a function that stabilize variance.
For instance, with the gaussian-poissonian model:
f'(m)=\frac{1}{\sqrt{am+b}}
We can choose the following function:
f(m)=\frac{2}{a}\sqrt{am+b}
If we transform our data by applying this function, we will get noise variance of 1.
We can then denoise the data using our favorite denoising algorithm, and perform an inverse transformation after that.
In practice
Finding a model that fits perfectly reality and that is easy to compute is not easy.
At least we can get as close as we can.
Here is an example of a fit using the gaussian-poissonian model (see how variance varies with the mean):
Here is the variance of data transformed:
Here is another example of a fit using a different model, where we assume V(X)=a*(E[X]+0.0005)^b:
Here is the variance of the data transformed:
Here, the second model allows to have a better variance stabilization in the dark areas, which can make denoising of very noisy images easier.
In darktable, what we currently have is the first model.
In Rawtherapee, while AFAIK there is no noise profiles, there is still a variance stabilization transform done: the gamma parameter in the noise reduction module allows the user to find the transform to stabilize the variance himself.
A gamma equal to 2 will give roughly the same results as the variance stabilization we currently use in darktable.
If you find the right gamma value, you will get something similar of what we have with the second model shown here.
Note that depending on the gamma value chosen, the variance will not always be stabilized, but it can be