How is WSNR (weighted signal to noise ratio) defined? I'm studying about comparison between images to determine quality. I've learned about MSE, SNR and PSNR, and now I'm trying to understand WSNR, which I assume is similar to SNR but with weights.
I have a Python code to calculate the WSNR between two images, but I don't fully understand it and I don't know what is the mathematical definition of this measure. I can't find any paper or book that explain how it is defined and calculated.
Question:
How is WSNR (weighted signal to noise ratio) defined?
Python Code:
import numpy as np
from scipy.ndimage.filters import gaussian_filter as __gaussian_filter
from scipy.ndimage.filters import convolve as __convolve
from scipy.ndimage.filters import correlate as __correlate
from scipy.fftpack import fftshift as __fftshift

def wsnr(reference, query):
&quot;&quot;&quot;Computes the Weighted Signal to Noise Ratio (WSNR) metric.

    value = wsnr(reference, query)

    inputs
    ----------
    reference: original image data.
    query    : modified image data to be compared.

    output
    ----------
    value    : wsnr value
    &quot;&quot;&quot;
    def __genetate_meshgrid(x, y):
        f = lambda u: u / 2 + 0.5 - 1
        (H, W) = map(f, (x, y))
        return (H, W)

    def __create_complex_planes(x, y):
        (H, W) = __genetate_meshgrid(x, y)
        (xplane, yplane) = np.mgrid[-H:H + 1, -W:W + 1]
        return (xplane, yplane)

    def __get_evaluated_contrast_sensivity(plane):
        w = 0.7
        angle = np.angle(plane)
        return ((1.0 - w) / 2.0) * np.cos(4.0 * angle) + (1.0 + w) / 2.0

    def __get_radial_frequency(x, y):
        (xplane, yplane) = __create_complex_planes(x, y)
        nfreq = 60
        plane = (xplane + 1.0j * yplane) / x * 2.0 * nfreq
        s = __get_evaluated_contrast_sensivity(plane)
        radfreq = abs(plane) / s
        return radfreq

    def __generate_CSF(radfreq):
        a = -((0.114 * radfreq) ** 1.1)
        csf = 2.6 * (0.0192 + 0.114 * radfreq) * np.exp(a)
        f = radfreq &lt; 7.8909
        csf[f] = 0.9809
        return csf

    def __weighted_fft_domain(ref, quer, csf):
        err = ref.astype('double') - quer.astype('double')
        err_wt = __fftshift(np.fft.fft2(err)) * csf
        im = np.fft.fft2(ref)
        return (err, err_wt, im)

    def __get_weighted_error_power(err_wt):
        return (err_wt * np.conj(err_wt)).sum()

    def __get_signal_power(im):
        return (im * np.conj(im)).sum()

    def __get_ratio(mss, mse):
        if mse != 0:
            ratio = 10.0 * np.log10(mss / mse)
        else:
            ratio = float(&quot;inf&quot;)
        return np.real(ratio)

    if not len(reference.shape) &lt; 3:
        reference = __convert_to_luminance(reference)
        query = __convert_to_luminance(query)
    size = reference.shape
    (x, y) = (size[0], size[1])
    radfreq = __get_radial_frequency(x, y)
    csf = __generate_CSF(radfreq)
    (err, err_wt, im) = __weighted_fft_domain(reference, query, csf)
    mse = __get_weighted_error_power(err_wt)
    mss = __get_signal_power(im)
    ratio = __get_ratio(mss, mse)
    return ratio
There might be several versions, yet the most probable is as follows. An SNR is an energy-dependent measure. In the time or spatial domain, it is agnostic to location: all samples have the same "weight" (uniform weighting). As an energetic measure, it can be computed in the Fourier domain as well (or any orthogonal transformation). 

Then, from the sensory system, we know that in hearing, or vision, not all frequencies have the same perceptual weight: a human hearing a loud 35 kHz tone won't probably not notice it, while a bat or a dog will notice. Hence, the relative importance of a sound should be weighted by the perceptual sensitivity of the subject (bat, dog, human). This yields a different weighting across different frequencies. One shall not care about a 35 kHz tone for a human ear, and give it a little weight, compared to the [1-2] kHz band. Same happens for vision. Such weightings are used for instance in JPEG compression.

So instead of computing, in the frequency domain some:

$$-10 \log \sum_k \frac{\left(\hat{s_k}-\hat{n_k}\right)^2}{\hat{s_k}^2}$$
with $k$ an index of frequency bins, you use a predefined weighting for each $k$:
$$-10 \log \sum_k w_k\frac{\left(\hat{s_k}-\hat{n_k}\right)^2}{\hat{s_k}^2}\,,$$
where  $w_k$ give less weight to less important frequencies. In sound, you can find A-weighting curves (and B-, C-, D- and Z-). 

In images, a Matlab code function ratio = wsnr(orig, dith, nfreq) gives references to ancient papers:


J. Mannos and D. Sakrison, "The effects of a visual fidelity criterion on the encoding of images", IEEE Trans. Inf. Theory, July 1974.
T. Mitsa and K. Varkur, "Evaluation of contrast sensitivity functions for the formulation of quality measures incorporated in halftoning algorithms", IEEE Int. Conference on Acoustics, Speech and Signal Processing, 1993