多项式插值之Lagrange、PCHIP与Spline以及BD-Rate和BD-PSNR的计算

做视频编解码算法的话一般都会接触到一个目前比较通用的客观评价指标，即Bjøntegaard Delta（BD），如果关注码率的话就是BD-Rate，如果关注PSNR的话就是BD-PSNR，其实质就是求在一个指标相同时另外一个指标的平均变化。一般我们关注的是在同等客观质量下的码率变化，所以用得比较多的是BD-Rate。不过问题是怎么去求平均，一种简单的方法就是假设码率（R）和PSNR（D）具有一定的函数关系，即R=f(D)或者D=g( R)，一般来说认为是单调的，那么两种编解码算法就对应了两条曲线，要求码率或者PSNR的平均变化就可以通过求区间积分的差值来表示了。不过问题又来了，编解码算法那么复杂我们怎么去得到两者的函数关系，至于它们两者有没有函数关系都说不准，所以我们只能再进行一些简化，通常就是测几个QP的数据，得到几个离散的点，然后进行插值作为该曲线的估计。所以BD-Rate归根结底还是一个插值问题。插值的方法有很多，不过我们还要考虑它们的复杂度以及合理性，其中多项式插值倒是一个不错的选择。多项式插值我们听得最多的应该就是Lagrange插值、分段三次插值以及样条（Spline）插值了，所以借着学习BD-Rate的机会顺便把这些算法的原理搞懂，以后用到的地方可能也不少。

因为把word文档转成markdown太麻烦了，没什么时间写，里面公式也比较多，所以我就直接把文档转成图片放上来了，看上去还是挺清晰的。里面主要是把找到的一些资料总结了一下，总体还是比较详细容易理解的，最后有一些我使用的文献，有兴趣的可以自己研读一下。当然一般都是英文的了，这方面中文资料还是比较少。














参考资料：

1. https://www.mathworks.com/content/dam/mathworks/mathworks-dot-com/moler/interp.pdf
2. Monotone Piecewise Cubic Interpolation.
   https://zero.sci-hub.tw/1536/cb55e0aa6f5eba4f009e7e812f4fc953/fritsch1980.pdf#view=FitH
3. https://coast.nd.edu/jjwteach/www/www/30125/pdfnotes/lecture5_9v14.pdf
4. https://blogs.mathworks.com/cleve/2012/07/16/splines-and-pchips/
5. Calculation of average coding efficiency based on subjective quality scores.
   https://infoscience.epfl.ch/record/192802/files/Elsevier2013_preprint.pdf
6. http://www.uta.edu/faculty/krrao/dip/Courses/EE5351/BDPSNR.doc

Python脚本

# -*- coding: utf-8 -*-

from __future__ import print_function, division
import numpy as np
import scipy

def pchip_slopes(h, delta):
    d = np.zeros(len(h) + 1)
    k = np.argwhere(np.sign(delta[:-1]) * np.sign(delta[1:]) > 0).reshape(-1) + 1
    w1 = 2*h[k] + h[k-1]
    w2 = h[k] + 2*h[k-1]
    d[k] = (w1 + w2) / (w1 / delta[k-1] + w2 / delta[k])
    d[0] = pchip_end(h[0], h[1], delta[0], delta[1])
    d[-1] = pchip_end(h[-1], h[-2], delta[-1], delta[-2])
    return d

def pchip_end(h1, h2, del1, del2):
    d = ((2*h1 + h2)*del1 - h1*del2) / (h1 + h2)
    if np.sign(d) != np.sign(del1):
        d = 0
    elif np.sign(del1) != np.sign(del2) and np.abs(d) > np.abs(3*del1):
        d = 3 * del1
    return d

def spline_slopes(h, delta):
    a, r = np.zeros([3, len(h)+1]), np.zeros(len(h)+1)
    a[0, 1], a[0, 2:] = h[0] + h[1], h[:-1]
    a[1, 0], a[1, 1:-1], a[1, -1] = h[1], 2*(h[:-1] + h[1:]), h[-2]
    a[2, :-2], a[2, -2] = h[1:], h[-1] + h[-2]
    
    r[0] = ((h[0] + 2*a[0, 1])*h[1]*delta[0] + h[0]**2*delta[1]) / a[0, 1]
    r[1:-1] = 3*(h[1:] * delta[:-1] + h[:-1] * delta[1:])
    r[-1] = (h[-1]**2*delta[-2] + (2*a[2, -2] + h[-1])*h[-2]*delta[-1]) / a[2, -2]
    
    d = scipy.linalg.solve_banded((1, 1), a, r)
    return d

class PCHIP:
    def __init__(self, x, y, use_spline=False):
        assert len(np.unique(x)) == len(x)
        order = np.argsort(x)
        self.xi, self.yi = x[order], y[order]
        
        h = np.diff(self.xi)
        delta = np.diff(self.yi) / h
        
        self.d = spline_slopes(h, delta) if use_spline else pchip_slopes(h, delta)
        self.c = (3*delta - 2*self.d[:-1] - self.d[1:]) / h
        self.b = (self.d[:-1] - 2*delta + self.d[1:]) / h**2
        
        """
        The piecewise function is like p(x) = y_k + s*d_k + s*s*c_k + s*s*s*b_k
        where s = x - xi_k, k is the interval includeing x.
        So the original function of p(x) is P(x) = xi_k*y_k + s*y_k + 1/2*s*s*d_k + 1/3*s*s*s*c_k + 1/4*s*s*s*s*b_k + C.
        """
        self.interval_int_coeff = []
        self.interval_int = np.zeros(len(x)-1)
        for i in range(len(x)-1):
            self.interval_int_coeff.append(np.polyint([self.b[i], self.c[i], self.d[i], self.yi[i]]))
            self.interval_int[i] = np.polyval(self.interval_int_coeff[-1], h[i]) - np.polyval(self.interval_int_coeff[-1], 0)
    
    def interp(self, x):
        if len(x[x < np.min(self.xi)]) > 0 or len(x[x > np.max(self.xi)]) > 0:
            print('Warning: Some samples are out of the interval and the results may be strange!')
        
        """find the intervals the xs belong to"""
        k = np.zeros(len(x), dtype='int')
        for i in range(1, len(self.xi)-1):
            idx = np.argwhere(self.xi[i] <= x).reshape(-1)
            k[idx] = i * np.ones(len(idx), dtype='int')
        s = x - self.xi[k]
        y = self.yi[k] + s*(self.d[k] + s*(self.c[k] + s*self.b[k]))
        return y
    
    def _integral(self, lower, upper):
        assert lower <= upper
        if lower < np.min(self.xi):
            lower = np.min(self.xi)
            print('Warning: The lower bound is less than the interval and clipped!')
        elif lower > np.max(self.xi):
            print('Warning: The lower bound is greater than the interval!')
            return 0
        if upper > np.max(self.xi):
            upper = np.max(self.xi)
            print('Warning: The upper bound is greater than the interval and clipped!')
        elif upper < np.min(self.xi):
            print('Warning: The lower bound is less than the interval!')
            return 0
        left = np.arange(len(self.xi))[self.xi - lower > -1e-6][0]
        right = np.arange(len(self.xi))[self.xi - upper < 1e-6][-1]
        
        inte = np.sum(self.interval_int[left:right])
        if self.xi[left] - lower > 1e-6:
            inte += (np.polyval(self.interval_int_coeff[left-1], self.xi[left]-self.xi[left-1]) - np.polyval(self.interval_int_coeff[left-1], lower-self.xi[left-1]))
        if self.xi[right] - upper < -1e-6:
            inte += (np.polyval(self.interval_int_coeff[right], upper-self.xi[right]) - np.polyval(self.interval_int_coeff[right], 0))
        return inte
    
    def integral(self, lower, upper):
        if lower > upper:
            return -self._integral(upper, lower)
        else:
            return self._integral(lower, upper)

def BD_Rate(R1, PSNR1, R2, PSNR2, piecewise=True):
    lR1, lR2 = np.log10(R1), np.log10(R2)
    
    min_int = np.max((np.min(PSNR1), np.min(PSNR2)))
    max_int = np.min((np.max(PSNR1), np.max(PSNR2)))

    if piecewise == True:
        int1 = PCHIP(PSNR1, lR1, use_spline=False).integral(min_int, max_int)
        int2 = PCHIP(PSNR2, lR2, use_spline=False).integral(min_int, max_int)
    else:
        p1 = np.polyfit(PSNR1, lR1, len(lR1)-1)
        p2 = np.polyfit(PSNR2, lR2, len(lR1)-1)
        p_int1 = np.polyint(p1)
        p_int2 = np.polyint(p2)

        int1 = np.polyval(p_int1, max_int) - np.polyval(p_int1, min_int)
        int2 = np.polyval(p_int2, max_int) - np.polyval(p_int2, min_int)
    
    avg_exp_diff = (int2 - int1) / (max_int - min_int)
    avg_diff = (np.power(10, avg_exp_diff) - 1) * 100.
    return avg_diff

def BD_PSNR(R1, PSNR1, R2, PSNR2, piecewise=True):
    lR1, lR2 = np.log10(R1), np.log10(R2)
    
    min_int = np.max((np.min(lR1), np.min(lR2)))
    max_int = np.min((np.max(lR1), np.max(lR2)))

    if piecewise == True:
        int1 = PCHIP(lR1, PSNR1, use_spline=False).integral(min_int, max_int)
        int2 = PCHIP(lR2, PSNR2, use_spline=False).integral(min_int, max_int)
    else:
        p1 = np.polyfit(lR1, PSNR1, len(lR1)-1)
        p2 = np.polyfit(lR2, PSNR2, len(lR1)-1)
        p_int1 = np.polyint(p1)
        p_int2 = np.polyint(p2)

        int1 = np.polyval(p_int1, max_int) - np.polyval(p_int1, min_int)
        int2 = np.polyval(p_int2, max_int) - np.polyval(p_int2, min_int)
    
    avg_diff = (int2 - int1) / (max_int - min_int)
    return avg_diff

if __name__ == '__main__':
    '''
    # x 假设是单调的，内部会进行排序，不然就不是函数了
    x = np.array([1, 2, 3, 4, 5, 6])
    y = np.array([16., 18., 21., 17., 15., 12.])
    xx = np.linspace(x[0]-0.25, x[-1]+0.25, 101)
    
    interp = PCHIP(x, y, use_spline=True)
    spline_y = interp.interp(xx)
    
    p = np.polyfit(x, y, 5)
    poly_y = np.polyval(p, xx)
    
    interp2 = PCHIP(x, y, use_spline=False)
    pchip_y = interp2.interp(xx)
    
    plt.plot(x, y, 'ro')
    plt.plot(x, y, '--', linewidth=1, label='linear')
    plt.plot(xx, poly_y, label='Lagrange')
    plt.plot(xx, pchip_y, label='sppchip')
    plt.plot(xx, spline_y, label='spline')
    plt.legend()
    plt.show()
    '''
    rate1 = np.array([514.2076, 337.1483, 225.6044, 148.7819])
    psnr1 = np.array([51.2271, 47.5517, 43.9461, 40.1898])
    
    rate2 = np.array([513.5298, 337.4435, 225.9397, 148.7625])
    psnr2 = np.array([51.2174, 47.5473, 43.9479, 40.1754])
    
    bd_psnr = BD_PSNR(rate1, psnr1, rate2, psnr2, piecewise=True)
    bd_rate = BD_Rate(rate1, psnr1, rate2, psnr2, piecewise=True)
    print(bd_psnr, bd_rate)
    bd_psnr1 = BD_PSNR(rate1, psnr1, rate2, psnr2, piecewise=False)
    bd_rate1 = BD_Rate(rate1, psnr1, rate2, psnr2, piecewise=False)
    print(bd_psnr1, bd_rate1)