astro-informatics / purify

Next-generation radio interferometric imaging.
https://astro-informatics.github.io/purify
GNU General Public License v2.0
17 stars 13 forks source link

Calculation of Epsilon #36

Closed Luke-Pratley closed 8 years ago

Luke-Pratley commented 8 years ago

In cpp/regressions/sdmm.cc (line 85), the calculation of epsilon follows: auto const epsilon = std::sqrt(y0.size() + 2 * std::sqrt(y0.size())) * sigma(y0, 30e0) / std::sqrt(static_cast<t_real>(y0.size()) / static_cast<t_real>(M31.size())); Here there is a division by the square root of the number of coefficients in y0 and in the number of pixels in the image of M31. I don't understand why this division is done.

Especially since the 2014 Purify paper (below Eq. 32) suggests:

auto const epsilon = std::sqrt(y0.size() + 2 * std::sqrt(y0.size())) * sigma

I have found following the Purify paper works better so far, but I am not sure if I am doing it correctly. Does anyone know the answer to this? I can try to answer any questions if it helps to clear anything up.

Luke-Pratley commented 8 years ago

Further more, the calculation of sigma for the injected noise for a given input signal to noise ratio (ISNR) in the 2014 purify paper follows sigma = || y_0 ||_{l_2}/sqrt(M) * 10^(-ISNR/20), where y_0 is the visibility vector, and M is the number of visibilities.

However, it looks like this formula can be derived from the formula for the reconstructed SNR, but this SNR formula applies to the reconstructed image, not the visibility vector. Why should the formula for SNR apply to both the image and the visibilities?

rafael-carrillo commented 8 years ago

Hi Luke, the second formula you are describing in your first email is the correct one.

Now, how you compute sigma is another story. You can define the input SNR as the ratio between the power of the visibilities (like the one you described in your las email) and the noise power, or, alternatively you can define it as the ration between the original image power and the noise power, i.e. sigma = || x0 ||{l_2}/sqrt(N) * 10^(-ISNR/20).

I hope this helps.

On 17 Apr 2016, at 00:01, Luke Pratley notifications@github.com wrote:

Further more, the calculation of sigma for the injected noise for a given input signal to noise ratio (ISNR) in the 2014 purify paper follows sigma = || y0 ||{l_2}/sqrt(M) * 10^(-ISNR/20), where y_0 is the visibility vector, and M is the number of visibilities.

However, it looks like this formula can be derived from the formula for the reconstructed SNR, but this SNR formula applies to the reconstructed image, not the visibility vector. Why should the formula for SNR apply to both the image and the visibilities?

— You are receiving this because you were assigned. Reply to this email directly or view it on GitHub https://github.com/astro-informatics/purify/issues/36#issuecomment-210912674

Luke-Pratley commented 8 years ago

@rafael-carrillo I just realised that the ISNR should be the same for both visibility and image domains because the Fourier transform is unitary, so the l_2 norm is preserved in the formula.

Though, in our case it depends on the measurement operator preserving the l_2 norm when mapping from visibilities to an image. I think this could be true, but don't know for sure.

rafael-carrillo commented 8 years ago

Though, in our case it depends on the measurement operator preserving the l_2 norm when mapping from visibilities to an image. I think this could be true, but don't know for sure.

This is not always true in our case since we have incomplete Fourier information, thus (up to numerical scaling factors) the energy in the visibilities is always less than the image energy.

On 18 Apr 2016, at 15:10, Luke Pratley notifications@github.com wrote:

@rafael-carrillo https://github.com/rafael-carrillo I just realised that the ISNR should be the same for both visibility and image domains because the Fourier transform is unitary, so the l_2 norm is preserved in the formula.

Though, in our case it depends on the measurement operator preserving the l_2 norm when mapping from visibilities to an image. I think this could be true, but don't know for sure.

— You are receiving this because you were mentioned. Reply to this email directly or view it on GitHub https://github.com/astro-informatics/purify/issues/36#issuecomment-211373572

Luke-Pratley commented 8 years ago

@rafael-carrillo , that is a good point, that would mean that the SNR is not always preserved when degridding. But for gridding, the dirty image should have the same energy as the visibilities, right?

rafael-carrillo commented 8 years ago

Yes, if the operator is normalised, i.e. the operator norm is one. In theory the Fourier operator is normal but since we are using a discrete approximation for the continuos operator this really depends on the implementation.

On 18 Apr 2016, at 15:23, Luke Pratley notifications@github.com wrote:

@rafael-carrillo https://github.com/rafael-carrillo , that is a good point, that would mean that the SNR is not always preserved when degridding. But for gridding, the dirty image should have the same energy as the visibilities, right?

— You are receiving this because you were mentioned. Reply to this email directly or view it on GitHub https://github.com/astro-informatics/purify/issues/36#issuecomment-211377041

Luke-Pratley commented 8 years ago

Cool thanks, I think this issue is settled now! Thanks!