<p>Let <span class="math-container">$h_j$</span> denote the image reconstruction filter of pixel <span class="math-container">$j$</span>.
I'm estimating the color value <span class="math-container">$$I_j=\int h_jf\:{\rm d}\mu$$</span> of the <span class="math-container">$j$</span>th pixel (see <a href="http://graphics.stanford.edu/papers/veach_thesis/thesis.pdf" rel="nofollow noreferrer">[Veach, Section 8.2]</a>) by an asymptotically consistent estimator whose asymptotic variance is given by <span class="math-container">$\sigma^2_j$</span>. (This has to be understood in a component-wise sense, i.e. <span class="math-container">$\sigma_j^2$</span> is a 3-dimensional vector whose <span class="math-container">$i$</span>th component is the asymptotic variance of the estimate of the <span class="math-container">$i$</span>th component of <span class="math-container">$I_j$</span>.</p>
<blockquote>
<p><strong>Question 1</strong>: I want to compare different estimators which lead to different asymptotic variances and would like to visualize that one is superior to another (with respect to variance reduction). How can I do that? Simply writing the variances into an image doesn't seem to be sensible. (The problem is the scale of the values.)</p>
<p><strong>Question 2</strong>: For simplicity, assume a box filter with radius <span class="math-container">$\frac12$</span> so that the computation of <span class="math-container">$I_j$</span> only involves the <span class="math-container">$j$</span>th pixel. Instead of considering the variances of each pixel, I've seen that we could consider the variances of a &quot;local window&quot; (for example a <span class="math-container">$3\times 3$</span> pixel window). Is this a better or worse metric for the total variance of the image?</p>
</blockquote>

CC BY-SA 4.0
2020-04-04T09:10:19.807Z
The variance is computed per pixel - you are estimating a different integral in each pixel. It is reasonable to have a different variance for different pixels.
lightxbulb
CC BY-SA 4.0
2020-04-04T09:21:32.030Z
For Q1 you could just plot graphs for the same integrand at different sample counts. That is - measure the variance at 1 pixel for different sample counts and plot it.
lightxbulb
CC BY-SA 4.0
2020-04-04T10:17:04.963Z
@lightxbulb It doesn't really matter, but you've got me wrong: I've got different estimators and so I'm comparing them with a fixed sample size. But what exactly do you mean by "plotting the graph"? Currently I'm writing the pixel variances as an image into a file, but that doesn't really make sense. For example, if a variance value in the red-channel is $\ge1$ it will always be "purely red". So, I would need to apply some scaling or something like that.
0xbadf00d
CC BY-SA 4.0
2020-04-04T10:56:47.970Z
Pick one pixel and use one estimator, then compute the variance at different sample counts, store the sample counts and the corresponding variance at each in a text file. Do the same for a different estimator. Make a plot of sample count vs variance for both (you can use matplotlib in Python or Matlab). This way you'll be able to compare the asymptotic behaviour of the two for the integrand at the chosen pixel. You can use a log-log scale if the difference warrants it. It's a common way of comparing low discrepancy vs random sequence performance.
lightxbulb
CC BY-SA 4.0
2020-04-04T11:01:35.580Z
@lightxbulb Okay, but this would be a visualization for the change of variance of a single pixel, right? Is there any common approach to visualize the variance (for a fixed estimator; not in comparison between two of them) of the whole image? I thought it might be possible to visualize the per-pixel variance by an image itself.
0xbadf00d
CC BY-SA 4.0
2020-04-04T13:10:08.253Z
You can visualize the per-pixel variance of the image: just output an HDR image - it is not limited to 0-1. You can use stb_image_write.
lightxbulb
CC BY-SA 4.0
2020-04-04T15:43:30.530Z
@lightxbulb Okay, I need to admit that I have no theoretical background in HDR. I've tested `stb_image_write` by `std::vector<float> image(3 * width * height);` `float r = 0; for (std::size_t i = 0; i < width * height; ++i, r += .01) image[3 * i] = r;` and `stbi_write_hdr("test.exr", width, height, 3, image.data());`. The resulting image has a red channel ranging from $0$ to $3584$. How does this help with my problem? It seems like the scaling problem is still present (since I should obtain the same image output if I divide by 3584 before writing the image.)
0xbadf00d
CC BY-SA 4.0
2020-04-04T16:06:40.713Z
@0xbadf00d you can save the variance as floats for each channel, save the full image as an EXR and do the visualization in Python using some colormap. You can use [`pyexr`](https://github.com/tvogels/pyexr) to load the image into a Numpy array–which you can rescale as you wish–and save to an LDR (e.g. PNG) to visualize the per-pixel variance.
Hubble
CC BY-SA 4.0
2020-04-04T17:02:25.297Z
@Hubble Currently I'm opening the EXR using RenderDoc where I can rescale as well. So, I guess using Python isn't superior to that or am I missing something? I didn't get your last sentence. How would I need to rescale so that I'm able to compare two variance images? Dividing by the max component value wouldn't make sense ... (Is there any good reading on this topic?)
0xbadf00d
CC BY-SA 4.0
2020-04-04T17:20:17.617Z
@Hubble Oh, and is there a useful metric which indicates that the "overall variance" of one estimator $E_1$ is smaller than the overall variance of another estimator $E_2$? It might clearly be the case that fixed pixels $j_1,j_2$ estimator $E_1$ has smaller variance in $j_1$, but larger variance in $j_2$. Is there a metric which helps to decide which estimator is superior for the whole image?
0xbadf00d
CC BY-SA 4.0
2020-04-04T17:44:36.320Z
Try visualizing the variance through a program like tev: https://github.com/Tom94/tev
As for "overall" variance, I am not aware of a specific metric. The simplest thing you could do is take the average, but obviously you lose some information this way - it's perfectly reasonable that one estimator is better in one part of the image, while the other is better in another part. What's the purpose of this exercise?
lightxbulb
CC BY-SA 4.0
2020-04-04T17:48:24.367Z
@lightxbulb My estimator depends on a parameter which I want to tune for variance reduction.  The problem is that - in contrast to bidirectional path tracing where the parameters are the multiple importance sampling weights - my (asymptotic) variances depend on the individual pixel and there is no easy bound (as it is the case in BDPT yielding the balance heuristic).
0xbadf00d
CC BY-SA 4.0
2020-04-04T17:48:36.203Z
@0xbadf00d typically one would compute a reference image with several thousand samples and look at an averaged, pixel-wise comparison with your algorithm as a function of time. Popular metrics are [sMAPE](https://en.wikipedia.org/wiki/Symmetric_mean_absolute_percentage_error) and relMSE but other perceptual metrics also exist. If you plot the error over time (log-log) you can then track how fast your new estimator can reduce the variance compared to a baseline. See Figure 5 of [this paper](https://arxiv.org/abs/1704.06835) for an example in MCMC.
Hubble
CC BY-SA 4.0
2020-04-04T17:52:25.823Z
This [CGF paper](https://www.cg.tuwien.ac.at/research/publications/2019/celarek_adam-2019-qelta) from 2019 might give you some insights on what you are trying to accomplish.
Hubble
CC BY-SA 4.0
2020-04-04T18:30:47.570Z
@Hubble Thank you for the suggestions. Did I get this right: sMAPE would be the average of my $\sigma_j^2$?
0xbadf00d
CC BY-SA 4.0
2020-04-04T18:50:43.583Z
@Hubble BTW, my parameter optimization problem would be easier to solve if I would replace $f$ by the luminosity function in the definition of $I_j$, i.e. if I would consider the per-pixel variance of the relative color luminosity instead of the RGB values. Would this still make sense?
0xbadf00d
CC BY-SA 4.0
2020-04-06T04:59:42.433Z
@lightxbulb Regarding the loss of information in the "overall" variance: Do you think it is a good idea to remedy this issue by subdividing the image into subimages and running a separate Markov chain on each subimage? Or is this tiling known to produce artifacts?
0xbadf00d
CC BY-SA 4.0
2020-04-06T07:15:53.650Z
What's your issue with considering only per pixel variance?
lightxbulb
CC BY-SA 4.0
2020-04-06T07:52:07.520Z
@lightxbulb I want to optimize a parameter of the Markov chain for variance reduction. Usually, the Markov chain runs over the whole image or a subregion of it. Unless we run it on a single pixel, I cannot consider per-pixel variance but need some variance notion for a set of pixels.
0xbadf00d
CC BY-SA 4.0
2020-04-06T15:55:05.900Z
If you do tiling, I would assume that this will result in a visible difference of the error between two tiles, so I would recommend not to do it. You can still try it to see whether my assumption is correct (that the difference will be too large, and thus visible). Tiling in general introduces discontinuities.
lightxbulb
CC BY-SA 4.0
2020-04-06T18:06:25.073Z
@lightxbulb I think there shouldn't be a visible difference of the error as long as each tile has the same area and the chains all have the same sample size. (And a first simple test was positive.)
0xbadf00d
CC BY-SA 4.0
2020-04-06T21:00:52.870Z
Do post some screen after you do some experiments. I am interested.
lightxbulb
CC BY-SA 4.0
2020-04-07T09:54:26.630Z
@lightxbulb I will try some things but I think we need to remember that the whole reason for using MCMC instead of convential "per-pixel" approaches like bidirectional path tracing is that it automatically finds regions which make a large contribution to the image and takes the majority of samples from these regions. In the degenerate case where our tiles contain only a single pixel, my problem of a suitable variance notion would disappear, but would the resulting algorithm (which runs a separate chain on each pixel) still make sense?
0xbadf00d
CC BY-SA 4.0
2020-04-08T08:52:22.840Z
@lightxbulb Since Hibble didn't respond to this question (yet): Does it make sense to consider the variances of the luminance function $p$ instead of the measurement contribution function $f$ itself ($p=\alpha_1f_1+\alpha_2f_2+\alpha_3f_3$)?
0xbadf00d
CC BY-SA 4.0
2020-04-08T09:25:32.100Z
Make sense for what purpose? What are the alphas and $f$s here?
lightxbulb
CC BY-SA 4.0
2020-04-08T10:09:19.013Z
@lightxbulb For minimizing the image variance. $f$ is the measurement contribution function and since $p$ is the luminance function, $p=0.212671f_1+0.715160f_2+0.072169f_3$. This would make the analysis easier since the target distribution density of MLT is $p$.
0xbadf00d
CC BY-SA 4.0
2020-04-08T10:13:00.203Z
You can do either. There's no specific prescription unless you are referring to a specific paper. You can even apply a nonlinearity if you so wish. At the end of the day this just affects the metric that you want to minimize, if you take perception into account this metric will obviously be a lot more complex.
lightxbulb
CC BY-SA 4.0
2020-04-08T15:27:15.297Z
@lightxbulb The integral corresponding the value of a pixel is of the form $\mu g:=\int g\:{\rm d}\mu$, where $g$ is a $[0,\infty)^3$-valued function on the path space $E$ and $\mu$ is the infinite-area product measure. So, it's not clear what the "variance" of an estimate to $\mu g$ should be (since it is vector-valued, not scalar-valued). However, if we know that our estimator $A_n$ is such that $A_n\psi\to\mathcal N(0,\sigma^2(\psi)$ for all scalar-valued $\psi$ on $E$, then has asymptotic variance $\sigma^2(\psi)$, then
0xbadf00d
CC BY-SA 4.0
2020-04-08T15:53:00.023Z
@lightxbulb $\langle a,\sqrt n(A_ng-\mu g)\rangle\to\mathcal N(0,\sigma^2(\langle a,g\rangle))$ for all $a\in\mathbb R^3$ and hence $\sqrt n(A_ng-\mu g)\to\mathcal N(0,\Sigma(g))$, where $\Sigma(g)$ is a nonnegative symmetric $3\times3$-matrix  whose trace norm is equal to $\sigma^2(|g|)$, where $|g|$ is the Euclidean norm of $g$. So, mathematically, it would make sense to try to minimize $\sigma^2(|g|)$, but I guess from a visual perception point-of-view it makes more sense to $\sigma^2(\ell(g))$, where $\ell$ is the luminance function. Do you agree on that?
0xbadf00d
<raytracing><pathtracing><image-processing><monte-carlo><image-comparison>
Visualize the variance of an image
0xbadf00d
# Visualize the variance of an image

(Asked by: username_0 on 2020-04-04)

Let $h\_j$ denote the image reconstruction filter of pixel $j$.
I'm estimating the color value $$I\_j=\int h\_jf\:{\rm d}\mu$$ of the $j$th pixel (see [[Veach, Section 8.2]](http://graphics.stanford.edu/papers/veach_thesis/thesis.pdf)) by an asymptotically consistent estimator whose asymptotic variance is given by $\sigma^2\_j$. (This has to be understood in a component-wise sense, i.e. $\sigma\_j^2$ is a 3-dimensional vector whose $i$th component is the asymptotic variance of the estimate of the $i$th component of $I\_j$.

> **Question 1**: I want to compare different estimators which lead to different asymptotic variances and would like to visualize that one is superior to another (with respect to variance reduction). How can I do that? Simply writing the variances into an image doesn't seem to be sensible. (The problem is the scale of the values.)
>
> **Question 2**: For simplicity, assume a box filter with radius $\frac12$ so that the computation of $I\_j$ only involves the $j$th pixel. Instead of considering the variances of each pixel, I've seen that we could consider the variances of a "local window" (for example a $3\times 3$ pixel window). Is this a better or worse metric for the total variance of the image?

**Tags:** raytracing, pathtracing, image-processing, monte-carlo, image-comparison

> **Comments:**
> **username_1**: The variance is computed per pixel - you are estimating a different integral in each pixel. It is reasonable to have a different variance for different pixels.  
> **username_1**: For Q1 you could just plot graphs for the same integrand at different sample counts. That is - measure the variance at 1 pixel for different sample counts and plot it.  
> **username_0**: @username_1 It doesn't really matter, but you've got me wrong: I've got different estimators and so I'm comparing them with a fixed sample size. But what exactly do you mean by "plotting the graph"? Currently I'm writing the pixel variances as an image into a file, but that doesn't really make sense. For example, if a variance value in the red-channel is $\ge1$ it will always be "purely red". So, I would need to apply some scaling or something like that.  
> **username_1**: Pick one pixel and use one estimator, then compute the variance at different sample counts, store the sample counts and the corresponding variance at each in a text file. Do the same for a different estimator. Make a plot of sample count vs variance for both (you can use matplotlib in Python or Matlab). This way you'll be able to compare the asymptotic behaviour of the two for the integrand at the chosen pixel. You can use a log-log scale if the difference warrants it. It's a common way of comparing low discrepancy vs random sequence performance.  
> **username_0**: @username_1 Okay, but this would be a visualization for the change of variance of a single pixel, right? Is there any common approach to visualize the variance (for a fixed estimator; not in comparison between two of them) of the whole image? I thought it might be possible to visualize the per-pixel variance by an image itself.  
> **username_1**: You can visualize the per-pixel variance of the image: just output an HDR image - it is not limited to 0-1. You can use stb_image_write.  
> **username_0**: @username_1 Okay, I need to admit that I have no theoretical background in HDR. I've tested `stb_image_write` by `std::vector<float> image(3 * width * height);` `float r = 0; for (std::size_t i = 0; i < width * height; ++i, r += .01) image[3 * i] = r;` and `stbi_write_hdr("test.exr", width, height, 3, image.data());`. The resulting image has a red channel ranging from $0$ to $3584$. How does this help with my problem? It seems like the scaling problem is still present (since I should obtain the same image output if I divide by 3584 before writing the image.)  
> **username_2**: @username_0 you can save the variance as floats for each channel, save the full image as an EXR and do the visualization in Python using some colormap. You can use [`pyexr`](https://github.com/tvogels/pyexr) to load the image into a Numpy array–which you can rescale as you wish–and save to an LDR (e.g. PNG) to visualize the per-pixel variance.  
> **username_0**: @username_2 Currently I'm opening the EXR using RenderDoc where I can rescale as well. So, I guess using Python isn't superior to that or am I missing something? I didn't get your last sentence. How would I need to rescale so that I'm able to compare two variance images? Dividing by the max component value wouldn't make sense ... (Is there any good reading on this topic?)  
> **username_0**: @username_2 Oh, and is there a useful metric which indicates that the "overall variance" of one estimator $E_1$ is smaller than the overall variance of another estimator $E_2$? It might clearly be the case that fixed pixels $j_1,j_2$ estimator $E_1$ has smaller variance in $j_1$, but larger variance in $j_2$. Is there a metric which helps to decide which estimator is superior for the whole image?  
> **username_1**: Try visualizing the variance through a program like tev: https://github.com/Tom94/tev As for "overall" variance, I am not aware of a specific metric. The simplest thing you could do is take the average, but obviously you lose some information this way - it's perfectly reasonable that one estimator is better in one part of the image, while the other is better in another part. What's the purpose of this exercise?  
> **username_0**: @username_1 My estimator depends on a parameter which I want to tune for variance reduction.  The problem is that - in contrast to bidirectional path tracing where the parameters are the multiple importance sampling weights - my (asymptotic) variances depend on the individual pixel and there is no easy bound (as it is the case in BDPT yielding the balance heuristic).  
> **username_2**: @username_0 typically one would compute a reference image with several thousand samples and look at an averaged, pixel-wise comparison with your algorithm as a function of time. Popular metrics are [sMAPE](https://en.wikipedia.org/wiki/Symmetric_mean_absolute_percentage_error) and relMSE but other perceptual metrics also exist. If you plot the error over time (log-log) you can then track how fast your new estimator can reduce the variance compared to a baseline. See Figure 5 of [this paper](https://arxiv.org/abs/1704.06835) for an example in MCMC.  
> **username_2**: This [CGF paper](https://www.cg.tuwien.ac.at/research/publications/2019/celarek_adam-2019-qelta) from 2019 might give you some insights on what you are trying to accomplish.  
> **username_0**: @username_2 Thank you for the suggestions. Did I get this right: sMAPE would be the average of my $\sigma_j^2$?  
> **username_0**: @username_2 BTW, my parameter optimization problem would be easier to solve if I would replace $f$ by the luminosity function in the definition of $I_j$, i.e. if I would consider the per-pixel variance of the relative color luminosity instead of the RGB values. Would this still make sense?  
> **username_0**: @username_1 Regarding the loss of information in the "overall" variance: Do you think it is a good idea to remedy this issue by subdividing the image into subimages and running a separate Markov chain on each subimage? Or is this tiling known to produce artifacts?  
> **username_1**: What's your issue with considering only per pixel variance?  
> **username_0**: @username_1 I want to optimize a parameter of the Markov chain for variance reduction. Usually, the Markov chain runs over the whole image or a subregion of it. Unless we run it on a single pixel, I cannot consider per-pixel variance but need some variance notion for a set of pixels.  
> **username_1**: If you do tiling, I would assume that this will result in a visible difference of the error between two tiles, so I would recommend not to do it. You can still try it to see whether my assumption is correct (that the difference will be too large, and thus visible). Tiling in general introduces discontinuities.  
> **username_0**: @username_1 I think there shouldn't be a visible difference of the error as long as each tile has the same area and the chains all have the same sample size. (And a first simple test was positive.)  
> **username_1**: Do post some screen after you do some experiments. I am interested.  
> **username_0**: @username_1 I will try some things but I think we need to remember that the whole reason for using MCMC instead of convential "per-pixel" approaches like bidirectional path tracing is that it automatically finds regions which make a large contribution to the image and takes the majority of samples from these regions. In the degenerate case where our tiles contain only a single pixel, my problem of a suitable variance notion would disappear, but would the resulting algorithm (which runs a separate chain on each pixel) still make sense?  
> **username_0**: @username_1 Since Hibble didn't respond to this question (yet): Does it make sense to consider the variances of the luminance function $p$ instead of the measurement contribution function $f$ itself ($p=\alpha_1f_1+\alpha_2f_2+\alpha_3f_3$)?  
> **username_1**: Make sense for what purpose? What are the alphas and $f$s here?  
> **username_0**: @username_1 For minimizing the image variance. $f$ is the measurement contribution function and since $p$ is the luminance function, $p=0.212671f_1+0.715160f_2+0.072169f_3$. This would make the analysis easier since the target distribution density of MLT is $p$.  
> **username_1**: You can do either. There's no specific prescription unless you are referring to a specific paper. You can even apply a nonlinearity if you so wish. At the end of the day this just affects the metric that you want to minimize, if you take perception into account this metric will obviously be a lot more complex.  
> **username_0**: @username_1 The integral corresponding the value of a pixel is of the form $\mu g:=\int g\:{\rm d}\mu$, where $g$ is a $[0,\infty)^3$-valued function on the path space $E$ and $\mu$ is the infinite-area product measure. So, it's not clear what the "variance" of an estimate to $\mu g$ should be (since it is vector-valued, not scalar-valued). However, if we know that our estimator $A_n$ is such that $A_n\psi\to\mathcal N(0,\sigma^2(\psi)$ for all scalar-valued $\psi$ on $E$, then has asymptotic variance $\sigma^2(\psi)$, then  
> **username_0**: @username_1 $\langle a,\sqrt n(A_ng-\mu g)\rangle\to\mathcal N(0,\sigma^2(\langle a,g\rangle))$ for all $a\in\mathbb R^3$ and hence $\sqrt n(A_ng-\mu g)\to\mathcal N(0,\Sigma(g))$, where $\Sigma(g)$ is a nonnegative symmetric $3\times3$-matrix  whose trace norm is equal to $\sigma^2(|g|)$, where $|g|$ is the Euclidean norm of $g$. So, mathematically, it would make sense to try to minimize $\sigma^2(|g|)$, but I guess from a visual perception point-of-view it makes more sense to $\sigma^2(\ell(g))$, where $\ell$ is the luminance function. Do you agree on that?  

---

## 0 answers

