Some VIFp adjustments proposal

[igv]

1. Use N/3 to calculate the standard deviation of Gaussian filter instead of N/5.

Float weights:

{ 0.02997965, 0.03786711, 0.04636316, 0.05502488, 0.06330243, 0.07059225, 0.07630779, 0.07995691, 0.08121165, 0.07995691, 0.07630779, 0.07059225, 0.06330243, 0.05502488, 0.04636316, 0.03786711, 0.02997965 },
{ 0.0629702, 0.0929025, 0.12264921, 0.14489292, 0.15317033, 0.14489292, 0.12264921, 0.0929025, 0.0629702 },
{ 0.13357471, 0.22921512, 0.27442033, 0.22921512, 0.13357471 },
{ 0.27406862, 0.45186276, 0.27406862 }

Integer weights:

1965, 2482, 3038, 3606, 4149, 4626, 5001, 5240, 5322, 5240, 5001, 4626, 4149, 3606, 3038, 2482, 1965
4127, 6088, 8038, 9496, 10038, 9496, 8038, 6088,4127
8754, 15022, 17984, 15022, 8754
17961, 29614, 17961

2. For downsampling use size=5, sigma=1.08 Gaussian filter for all scale levels (similar to binom5 that used in IW-SSIM) to speed up processing and reduce aliasing:

0.06760634, 0.24462097, 0.37554539, 0.24462097, 0.06760634
4431, 16031, 24612, 16031, 4431

3. Hard-code sigma_nsq for coarse scales.

Accuracy drops with sigma_nsq >2 or <1. Almost no difference with values in range [1;2]. So just hard-code it to 1 (experiments show this to be optimal). Makes adjustable sigma_nsq much more flexible.

[igv]

Comparison of the current downsampling filter vs the proposed (at scale 4):

At higher scale levels (2 and 3) the proposed is sharper and still without aliasing.

[li-zhi]

Thanks for looking into this. We plan to have a thorough investigation of the accuracy-speed tradeoff later this year, and will be considering your proposal.

[igv]

VIFp Python script with my adjustments (took it from here). Also fixed error<1 when comparing identical images and simplified.

import sys
from PIL import Image
import numpy as np
from scipy.ndimage import gaussian_filter

def vif(file1, file2):
    img1 = Image.open(file1).convert('RGB')
    img2 = Image.open(file2).convert('RGB')

    width, height = img1.size
    img1 = np.frombuffer(img1.tobytes(), dtype=np.uint8).reshape(height, width, 3) / 255
    img2 = np.frombuffer(img2.tobytes(), dtype=np.uint8).reshape(height, width, 3) / 255
    
    img1 = np.where(img1 > 0.04045, np.power((img1 + 0.055) / 1.055, 2.4),  img1 / 12.92)
    img2 = np.where(img2 > 0.04045, np.power((img2 + 0.055) / 1.055, 2.4),  img2 / 12.92)

    Y1 = 0.2126 * img1[:,:,0] + 0.7152 * img1[:,:,1] + 0.0722 * img1[:,:,2]
    Y2 = 0.2126 * img2[:,:,0] + 0.7152 * img2[:,:,1] + 0.0722 * img2[:,:,2]

    sigma_nsq = 0.5
    eps = 1e-5

    num = 0.0
    den = 0.0
    for scale in range(1, 5):

        N = 2**(4-scale+1) + 1
        sd, t = N/3.0, 1.4  # kernel radius = round(sd * truncate)

        if scale == 2:
            sigma_nsq = .1

        if scale > 1:
            Y1 = gaussian_filter(Y1, 1.08, truncate=1.5)[::2, ::2]
            Y2 = gaussian_filter(Y2, 1.08, truncate=1.5)[::2, ::2]

        L1 = np.where(Y1 > 0.008856, np.power(Y1, 1./3.) * 116 - 16, Y1 * 903.3)
        L2 = np.where(Y2 > 0.008856, np.power(Y2, 1./3.) * 116 - 16, Y2 * 903.3)

        mu1 = gaussian_filter(L1, sd, truncate=t)
        mu2 = gaussian_filter(L2, sd, truncate=t)
        mu1_sq = mu1 * mu1
        mu2_sq = mu2 * mu2
        mu1_mu2 = mu1 * mu2
        sigma1_sq = gaussian_filter(L1 * L1, sd, truncate=t) - mu1_sq
        sigma2_sq = gaussian_filter(L2 * L2, sd, truncate=t) - mu2_sq
        sigma12 = gaussian_filter(L1 * L2, sd, truncate=t) - mu1_mu2

        sigma1_sq[sigma1_sq<eps] = eps
        sigma2_sq[sigma2_sq<eps] = eps
        sigma12[sigma12<eps] = eps

        g = sigma12 / sigma1_sq
        sv_sq = sigma2_sq - g * sigma12

        g[sigma1_sq<sigma_nsq] = 1
        #sv_sq[sigma1_sq<sigma_nsq] *= eps
        #g[g>1] = 1
        sv_sq[sv_sq<0] = 0

        #if scale == 1:
        #    sigma_nsq = np.maximum(np.cbrt(sigma1_sq), sigma_nsq)

        num += np.sum(np.log2(1 + g * g * sigma1_sq / (sv_sq + sigma_nsq)))
        den += np.sum(np.log2(1 + sigma1_sq / sigma_nsq))

    vifp = num/den

    return vifp

def main():
    for arg in sys.argv[2:]:
        score = vif(sys.argv[1], arg)
        print(str(score) + "\t" + arg)

if __name__ == '__main__':
    main()

Accuracy improvement with all the changes: Noise 0.79->0.846, Actual 0.823->0.875, Simple 0.902->0.929, Exotic 0.546->0.645, Full 0.611->0.686

[igv]

Accuracy with fine scale sigma_nsq = 128 (maximum reasonable value, ~16 in Lab color space): Noise - 0.871, Actual - 0.897, Simple - 0.953, Exotic - 0.728, Full - 0.745. Would be higher than MS-SSIM if not for Exotic.

[igv]

Calculating variances in Lab color space looks like is slightly more accurate.