单RD点Bjontegaard计算问题：警告修复与结果可靠性验证

我用Python脚本计算两种视频编码设置的BD-Rate：

• 当使用4个RD点（R1/PSNR1为参考视频的RD点，R2/PSNR2为测试视频的RD点）时，脚本运行正常，代码如下：

from bjontegaard_metric import *    
R1 = np.array([686.76, 309.58, 157.11, 85.95])
PSNR1 = np.array([40.28, 37.18, 34.24, 31.42])
R2 = np.array([893.34, 407.8, 204.93, 112.75])
PSNR2 = np.array([40.39, 37.21, 34.17, 31.24])

print('BD-PSNR: ', BD_PSNR(R1, PSNR1, R2, PSNR2))
print('BD-RATE: ', BD_RATE(R1, PSNR1, R2, PSNR2))

• 但仅使用1个RD点时，代码如下：

from bjontegaard_metric import *
R1 = np.array([686.76])
PSNR1 = np.array([40.28])
R2 = np.array([893.34])
PSNR2 = np.array([40.39])
        
print('BD-PSNR: ', BD_PSNR(R1, PSNR1, R2, PSNR2))
print('BD-RATE: ', BD_RATE(R1, PSNR1, R2, PSNR2))

会触发警告：RankWarning: Polyfit may be poorly conditioned。由于每次编码器运行仅返回一对PSNR和码率结果，我需要对比两组单RD点数据，请问如何修复该警告？仅用1个RD点得到的结果是否可靠？

附原始完整脚本：

import numpy as np
import scipy.interpolate

def BD_PSNR(R1, PSNR1, R2, PSNR2, piecewise=0):
    lR1 = np.log(R1)
    lR2 = np.log(R2)

    PSNR1 = np.array(PSNR1)
    PSNR2 = np.array(PSNR2)

    p1 = np.polyfit(lR1, PSNR1, 3)
    p2 = np.polyfit(lR2, PSNR2, 3)

    # integration interval
    min_int = max(min(lR1), min(lR2))
    max_int = min(max(lR1), max(lR2))

    # find integral
    if piecewise == 0:
        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)
    else:
        lin = np.linspace(min_int, max_int, num=100, retstep=True)
        interval = lin[1]
        samples = lin[0]
        v1 = scipy.interpolate.pchip_interpolate(np.sort(lR1), PSNR1[np.argsort(lR1)], samples)
        v2 = scipy.interpolate.pchip_interpolate(np.sort(lR2), PSNR2[np.argsort(lR2)], samples)
        # Calculate the integral using the trapezoid method on the samples.
        int1 = np.trapz(v1, dx=interval)
        int2 = np.trapz(v2, dx=interval)

    # find avg diff
    avg_diff = (int2-int1)/(max_int-min_int)

    return avg_diff


def BD_RATE(R1, PSNR1, R2, PSNR2, piecewise=0):
    lR1 = np.log(R1)
    lR2 = np.log(R2)

    # rate method
    p1 = np.polyfit(PSNR1, lR1, 3)
    p2 = np.polyfit(PSNR2, lR2, 3)

    # integration interval
    min_int = max(min(PSNR1), min(PSNR2))
    max_int = min(max(PSNR1), max(PSNR2))

    # find integral
    if piecewise == 0:
        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)
    else:
        lin = np.linspace(min_int, max_int, num=100, retstep=True)
        interval = lin[1]
        samples = lin[0]
        v1 = scipy.interpolate.pchip_interpolate(np.sort(PSNR1), lR1[np.argsort(PSNR1)], samples)
        v2 = scipy.interpolate.pchip_interpolate(np.sort(PSNR2), lR2[np.argsort(PSNR2)], samples)
        # Calculate the integral using the trapezoid method on the samples.
        int1 = np.trapz(v1, dx=interval)
        int2 = np.trapz(v2, dx=interval)

    # find avg diff
    avg_exp_diff = (int2-int1)/(max_int-min_int)
    avg_diff = (np.exp(avg_exp_diff)-1)*100
    return avg_diff

出现RankWarning的核心原因：脚本默认用3次多项式拟合RD曲线，但3次多项式需要至少4个不同的点才能唯一确定参数。当只有1个点时，多项式拟合的自由度无限大，计算条件极差，因此Numpy抛出警告。

有两种可行的修复思路：

思路1：根据RD点数量自动调整多项式阶数 #

修改函数，根据输入的RD点数量选择合适的多项式阶数：

• 1个点：用0次多项式（常数拟合）
• 2个点：用1次多项式（线性拟合）
• 3个点：用2次多项式
• 4个及以上点：保留3次多项式拟合

同时处理单点场景下积分区间为0的特殊情况，直接返回局部差值。修改后的函数示例：

import numpy as np
import scipy.interpolate

def BD_PSNR(R1, PSNR1, R2, PSNR2, piecewise=0):
    lR1 = np.log(R1)
    lR2 = np.log(R2)

    PSNR1 = np.array(PSNR1)
    PSNR2 = np.array(PSNR2)

    # 根据点数自动选择多项式阶数
    order1 = min(len(lR1)-1, 3)
    order2 = min(len(lR2)-1, 3)
    p1 = np.polyfit(lR1, PSNR1, order1)
    p2 = np.polyfit(lR2, PSNR2, order2)

    # integration interval
    min_int = max(min(lR1), min(lR2))
    max_int = min(max(lR1), max(lR2))

    # 处理单点情况：积分区间长度为0，直接返回PSNR差值
    if max_int == min_int:
        return PSNR2[0] - PSNR1[0]

    # find integral
    if piecewise == 0:
        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)
    else:
        lin = np.linspace(min_int, max_int, num=100, retstep=True)
        interval = lin[1]
        samples = lin[0]
        v1 = scipy.interpolate.pchip_interpolate(np.sort(lR1), PSNR1[np.argsort(lR1)], samples)
        v2 = scipy.interpolate.pchip_interpolate(np.sort(lR2), PSNR2[np.argsort(lR2)], samples)
        int1 = np.trapz(v1, dx=interval)
        int2 = np.trapz(v2, dx=interval)

    avg_diff = (int2-int1)/(max_int-min_int)
    return avg_diff


def BD_RATE(R1, PSNR1, R2, PSNR2, piecewise=0):
    lR1 = np.log(R1)
    lR2 = np.log(R2)

    PSNR1 = np.array(PSNR1)
    PSNR2 = np.array(PSNR2)

    # 根据点数自动选择多项式阶数
    order1 = min(len(PSNR1)-1, 3)
    order2 = min(len(PSNR2)-1, 3)
    p1 = np.polyfit(PSNR1, lR1, order1)
    p2 = np.polyfit(PSNR2, lR2, order2)

    # integration interval
    min_int = max(min(PSNR1), min(PSNR2))
    max_int = min(max(PSNR1), max(PSNR2))

    # 处理单点情况：积分区间长度为0，直接返回码率差值百分比
    if max_int == min_int:
        return (np.exp(lR2[0] - lR1[0]) - 1) * 100

    # find integral
    if piecewise == 0:
        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)
    else:
        lin = np.linspace(min_int, max_int, num=100, retstep=True)
        interval = lin[1]
        samples = lin[0]
        v1 = scipy.interpolate.pchip_interpolate(np.sort(PSNR1), lR1[np.argsort(PSNR1)], samples)
        v2 = scipy.interpolate.pchip_interpolate(np.sort(PSNR2), lR2[np.argsort(PSNR2)], samples)
        int1 = np.trapz(v1, dx=interval)
        int2 = np.trapz(v2, dx=interval)

    avg_exp_diff = (int2-int1)/(max_int-min_int)
    avg_diff = (np.exp(avg_exp_diff)-1)*100
    return avg_diff

思路2：强制使用分段插值模式（piecewise=1） #

脚本中已经实现了PCHIP分段插值的选项，这种方式在点数较少时更稳定，不会触发多项式拟合的警告。调用时传入piecewise=1即可：

print('BD-PSNR: ', BD_PSNR(R1, PSNR1, R2, PSNR2, piecewise=1))
print('BD-RATE: ', BD_RATE(R1, PSNR1, R2, PSNR2, piecewise=1))

仅用1个RD点得到的结果完全不可靠，原因如下：

• BD-Rate/BD-PSNR的核心是对比两条完整RD曲线在指定区间内的积分差异，反映的是编码设置在整个码率/质量范围内的平均效率差异。
• 单个RD点只能反映该特定码率（或PSNR）下的局部差异，无法代表整体编码效率。比如在高码率下A比B好，但低码率下B比A好，单点结果会完全误导判断。

如果只能获取单个RD点，建议：

1. 调整编码器参数，生成至少3-4个不同码率/质量的RD点，这是计算BD指标的标准要求。
2. 若无法生成多点，直接对比单点的PSNR差值或码率百分比差异，不要用BD指标的名义输出结果。

内容的提问来源于stack exchange，提问作者Maverick