Modeling Generalized Rate-Distortion Functions

TL;DR

This paper introduces a novel interpolation and sampling approach to accurately model the generalized rate-distortion surface of videos, enabling efficient prediction and analysis of video quality at various bitrates, resolutions, and viewing conditions.

Contribution

We develop the RAMCT interpolation method and a probabilistic sampling scheme to efficiently and accurately model the GRD function from limited measurements.

Findings

The proposed model accurately reconstructs GRD surfaces with fewer samples.

The sampling scheme significantly reduces computational costs.

Applications include rate-distortion prediction and encoder comparison.

Abstract

Many multimedia applications require precise understanding of the rate-distortion characteristics measured by the function relating visual quality to media attributes, for which we term it the generalized rate-distortion (GRD) function. In this study, we explore the GRD behavior of compressed digital videos in a three-dimensional space of bitrate, resolution, and viewing device/condition. Our analysis on a large-scale video dataset reveals that empirical parametric models are systematically biased while exhaustive search methods require excessive computation time to depict the GRD surfaces. By exploiting the properties that all GRD functions share, we develop an Robust Axial-Monotonic Clough-Tocher (RAMCT) interpolation method to model the GRD function. This model allows us to accurately reconstruct the complete GRD function of a source video content from a moderate number of measurements. To further reduce the computational cost, we present a novel sampling scheme based on a probabilistic model and an information measure. The proposed sampling method constructs a sequence of quality queries by minimizing the overall informativeness in the remaining samples. Experimental results show that the proposed algorithm significantly outperforms state-of-the-art approaches in accuracy and efficiency. Finally, we demonstrate the usage of the proposed model in three applications: rate-distortion curve prediction, per-title encoding profile generation, and video encoder comparison.

Tables

Table 1. TABLE I: MSE performance of the competing GRD function models with different number of labeled samples selected by random sampling (RS) and the proposed uncertainty sampling (US).

sample #	Reciprocal [13]		Logarithmic [7]		PCHIP		CT		RAMCT
	RS	US	RS	US	RS	US	RS	US	RS	US
20	N.A.	N.A.	23.07	13.33	68.76	26.49	88.54	56.04	135.27	16.10
30	62.27	83.34	13.08	10.56	30.95	2.06	37.99	22.78	10.98	3.29
50	38.11	73.88	9.43	6.77	8.64	0.07	11.75	12.16	4.70	0.06
75	30.27	48.85	5.15	4.92	3.08	0	4.84	3.26	1.01	0
100	27.44	38.46	4.60	4.18	1.77	0	2.75	1.26	0.13	0
540	24.51	24.51	2.76	2.76	0	0	0	0	0	0

Table 2. TABLE II: l ∞ subscript 𝑙 l_{\infty} performance of the competing GRD function models with different number of labeled samples selected by random sampling (RS) and the proposed uncertainty sampling (US).

sample #	Reciprocal [13]		Logarithmic [7]		PCHIP		CT		RAMCT
	RS	US	RS	US	RS	US	RS	US	RS	US
20	N.A.	N.A.	19.40	16.56	38.87	28.11	36.50	29.51	45.15	21.88
30	48.32	45.36	17.85	12.28	33.04	11.07	29.84	18.70	27.07	6.13
50	52.48	45.48	15.75	12.37	24.33	2.10	21.82	14.30	23.99	2.13
75	54.49	49.08	14.59	13.53	18.22	0.47	17.89	7.76	16.51	0.11
100	55.54	51.26	14.22	14.44	16.00	0.26	15.59	5.84	14.23	0
540	58.04	58.04	18.33	18.14	0	0	0	0	0	0

Table 3. TABLE III: l ∞ subscript 𝑙 l_{\infty} performance of the competing GRD function models with different number of labeled samples selected by random sampling (RS) and the proposed uncertainty sampling (US).

Resolution	$l_{\infty }$		MSE
	Linear	RAMCT	Linear	RAMCT
640$\times$360	7.68	3.12	7.89	1.56
960$\times$540	7.66	4.83	6.61	3.14
1600$\times$900	8.77	7.87	5.18	4.98
Average	8.04	5.27	6.56	3.23

Table 4. TABLE IV: Average bitrate saving of encoding profiles. Negative values indicate actual bitrate reduction

	Microsoft	Apple	Netflix	Proposed
Microsoft	0	-	-	-
Apple	-25.3%	0	-	-
Netflix	-29.3%	-5.6%	0	-
Proposed	-62.0%	-48.9%	-46.8%	0

Table 5. TABLE V: Performance of VP9, and HEVC in terms of the generalized Q g ​ a ​ i ​ n subscript 𝑄 𝑔 𝑎 𝑖 𝑛 Q_{gain} and R g ​ a ​ i ​ n subscript 𝑅 𝑔 𝑎 𝑖 𝑛 R_{gain} models with H.264 as the baseline codec

Codec	$Q_{gain}$	$R_{gain}$
VP9	1.9	-27.5%
HEVC	0.2	-15.9%

Equations

$$W_{GRD}:=$$

$$\mathrm{and}\forall x_a<x_b,f(x_a,y)<f(x_b,y)\},$$

$$z(\alpha ,\beta ,\gamma )=$$

$$3c_{T_{10}}\alpha {\beta }^2+3c_{I_{11}}{\beta }^2\gamma +c_S{\gamma }^3+3c_{I_{02}}\alpha {\gamma }^2+$$

$$3c_{I_{12}}\beta {\gamma }^2+6c_{C_2}\alpha \beta \gamma ,$$

$${\alpha }_P=\frac{A_{P{V}_{1}S}}{A_{V_0{V}_{1}S}},{\beta }_P=\frac{A_{PS{V}_0}}{A_{V_{1}S{V}_0}},{\gamma }_P=\frac{A_{P{V}_0{V}_1}}{A_{S{V}_0{V}_1}},$$

$$c_{P_i}=$$

$$c_{C_i}=$$

$$c_{I_{i2}}=$$

$$\frac{1}{3}[(y_{I_{i1}}-y_{V_i})+(y_{T_{ki}}-y_j^{\ast })+(y_{T_{ji}}-y_k^{\ast })]d_{V_i}^y+$$

$$\frac{1}{3}(x_{T_{ij}}-x_k^{\ast })d_{V_j}^x+\frac{1}{3}(y_{T_{ij}}-y_k^{\ast })d_{V_j}^y+\frac{1}{3}{\eta }_k{d}_{E_k}^e+$$

$$\frac{1}{3}(x_{T_{ik}}-x_j^{\ast })d_{V_k}^x+\frac{1}{3}(y_{T_{ik}}-y_j^{\ast })d_{V_k}^y+\frac{1}{3}{\eta }_j{d}_{E_j}^e+$$

$$\frac{1}{3}[z_{V_i}+({\theta }_{ki}{z}_{V_i}+{\theta }_{ik}{z}_{V_k})+({\theta }_{ij}{z}_{V_j}+{\theta }_{ji}{z}_{V_i})]$$

$$c_S=$$

$$\frac{1}{9}\sum _{i=0}^2[(y_{I_{i1}}-y_{V_i})+2(y_{T_{ki}}-y_j^{\ast })+2(y_{T_{ji}}-y_k^{\ast })]d_{V_i}^y+$$

$$\frac{2}{9}\sum _{i=0}^2{\eta }_i{d}_{E_i}^e+\frac{1}{9}\sum _{i=0}^2[(1+2{\theta }_{ji}+2{\theta }_{ki})z_{V_i}],$$

$$P_i\in \{T_{ij},T_{ik},I_{i1}\},$$

$${\eta }_i=\sqrt{(x_{C_i}-x_i^{\ast }{)}^2+(y_{C_i}-y_i^{\ast }{)}^2},$$

$${\theta }_{kj}=\frac{x_{T_{kj}}-x_i^{\ast }}{x_{T_{kj}}-x_{T_{jk}}},$$

$${\theta }_{jk}=\frac{x_{T_{jk}}-x_i^{\ast }}{x_{T_{jk}}-x_{T_{kj}}},$$

$$\mathbf{c}=\mathbf{Rd}+\mathbf{f},$$

$$d_{E_i}^e+d_{{\bar{E}}_i}^e=0.$$

$$\displaystyle\leavevmode\resizebox{375.80542pt}{}{$(y_{S}-y_{V_{j}})c_{T_{ij}}+(y_{V_{i}}-y_{S})c_{T_{ji}}+(y_{V_{j}}-y_{V_{i}})c_{C_{k}}\leq 0$}.$$

$$\mathbf{Gc}\le \mathbf{h},$$

$$\mathbf{GRd}\le \mathbf{h}-\mathbf{Gf}.$$

$$\displaystyle\leavevmode\resizebox{422.77661pt}{}{$L_{V_{0}V_{1}V_{2}}=\frac{1}{2}\sum_{i=0}^{2}\int_{E_{i}}\left[\frac{\partial^{2}z}{\partial E_{i}^{2}}\right]^{2}ds_{E_{i}}+\sum_{i=0}^{2}\int_{\hat{E}_{i}}\left[\frac{\partial^{2}z}{\partial\hat{E}_{i}^{2}}\right]^{2}ds_{\hat{E}_{i}}$},$$

$${\int }_{E_i}{\left[\frac{{\partial }^{2}z}{\partial E_i^2}\right]}^{2}d{s}_{E_i}=\frac{1}{\Vert E_i{\Vert }^3}{\int }_0^1{\left[z_{E_i}^{{}^{\prime \prime }}(t)\right]}^{2}dt$$

$$=$$

$$\frac{-36}{\Vert E_i{\Vert }^3}(z_{V_j}{c}_{T_{j}k}+z_{V_k}{c}_{T_{k}j})+\frac{12}{\Vert E_i{\Vert }^3}(z_{V_j}^2+z_{V_k}^2+z_{V_j}{z}_{V_k})$$

$$=$$

$$\displaystyle\left[{\begin{array}[]{*{20}{c}}{\frac{{-36{z_{{V_{j}}}}}}{{{{\left\|{{E_{i}}}\right\|}^{3}}}}}&{\frac{{-36{z_{{V_{k}}}}}}{{{{\left\|{{E_{i}}}\right\|}^{3}}}}}\end{array}}\right]\left[{\begin{array}[]{*{20}{c}}{{c_{{T_{jk}}}}}\\ {{c_{{T_{kj}}}}}\end{array}}\right]+$$

$$\frac{12}{{\Vert E_i\Vert }^3}(z_{V_j}^2+z_{V_k}^2+z_{V_j}{z}_{V_k}),$$

$$\Vert E_i\Vert =\sqrt{(x_{V_j}-x_{V_k}{)}^2+(y_{V_j}-y_{V_k}{)}^2}$$

$${\int }_{{\hat{E}}_i}{\left[\frac{{\partial }^{2}z}{\partial {\hat{E}}_i^2}\right]}^{2}d{s}_{{\hat{E}}_i}=\frac{1}{\Vert {\hat{E}}_i{\Vert }^3}{\int }_0^1{\left[z_{{\hat{E}}_i}^{{}^{\prime \prime }}(t)\right]}^{2}dt$$

$$=$$

$$\frac{12z_{V_i}}{\Vert {\hat{E}}_i{\Vert }^3}(-3c_{I_{i1}}+c_S)+\frac{12}{\Vert {\hat{E}}_i{\Vert }^3}{z}_{V_i}^2$$

$$=$$

$$\displaystyle+\left[{\begin{array}[]{*{20}{c}}{\frac{{-36{z_{{V_{i}}}}}}{{{{\left\|{{{\hat{E}}_{i}}}\right\|}^{3}}}}}&0&{\frac{{12{z_{{V_{i}}}}}}{{{{\left\|{{{\hat{E}}_{i}}}\right\|}^{3}}}}}\end{array}}\right]\left[{\begin{array}[]{*{20}{c}}{{c_{{I_{i1}}}}}\\ {{c_{{I_{i2}}}}}\\ {{c_{S}}}\end{array}}\right]+\frac{{12z_{V_{i}}^{2}}}{{{{\left\|{{{\hat{E}}_{i}}}\right\|}^{3}}}},$$

$$\Vert {\hat{E}}_i\Vert =\sqrt{(x_S-x_{V_i}{)}^2+(y_S-y_{V_i}{)}^2}$$

$$L_{V_0{V}_1{V}_2}={\mathbf{c}}^T{\mathbf{U}}_{V_0{V}_1{V}_2}\mathbf{c}+{\mathbf{w}}_{V_0{V}_1{V}_2}^T\mathbf{c}+const,$$

$$L_{V_0{V}_1{V}_2}=$$

Full text

I Introduction

Rate-distortion theory provides the theoretical foundations for lossy data compression and are widely employed in image and video compression schemes. Many multimedia applications require precise measurements of rate-distortion functions to characterize source signal and maximize user Quality-of-Experience (QoE). Examples of applications that explicitly use rate-distortion measurements are codec evaluation [1], rate-distortion optimization [2], video quality assessment (VQA) [3], encoding representation recommendation [4, 5, 6, 7], and QoE optimization of streaming videos [8, 9].

Digital videos usually undergo a variety of transforms and processes in the content delivery chain, as shown in Fig. 1. To address the growing heterogeneity of display devices, contents, and access network capacity, source videos are encoded into different bitrates, spatial resolutions, frame rates, and bit depths before transmitted to the client. In an adaptive streaming video distribution environment [10], based on the bandwidth, buffering, and computation constraints, client devices adaptively select a proper video representation on a per-time segment basis to download and render. Each process influences the visual quality of a video in a different way, which can be jointly characterized by a generalized rate-distortion (GRD) function. In general, this attribute-distortion mapping comprises several complex factors, such as source content, operation mode/type of encoder, rendering system, and human visual system (HVS) characteristics. In this work, we assume the GRD surface is a function $f:\mathbb{R}^{2}\to\mathbb{R}^{J}$, where the input of the function is the video representation consisting of bitrate and spatial resolution, the output of the function is the perceptual video quality at different viewing condition, and $J$ represents the number of viewing devices, respectively. Furthermore, the GRD function is content- and encoder-dependent.

Despite the tremendous growth in computational multimedia over the last few decades, estimating a GRD function is difficult, expensive, and time-consuming. Specifically, probing the quality of a single sample in the GRD space involves sophisticated video encoding and quality assessment, both of which are expensive processes. For example, the recently announced highly competitive AV1 video encoder [11] and video quality assessment model VMAF [12] could be over 1000 times and 10 times slower than real-time for full high-definition (1920$\times$1080) video content. Given the massive volume of multimedia data on the Internet, the real challenge is to produce an accurate estimate of the GRD function with a minimal number of samples.

We aim to develop a GRD function estimation framework with three desirable properties:

• •

Accuracy: It produces asymptotically unbiased estimation of GRD function, independent of the source video complexity and the encoder mechanism.

• •

Speed: It requires a minimal number of samples to reconstruct a full GRD function.

• •

Mathematical soundness: The GRD model has to be mathematically well-behaved, making it readily applicable to a variety of computational multimedia applications.

To achieve accuracy, we analyze the properties that all GRD functions share, based on which we formulate the GRD function approximation problem as a quadratic programming problem. The solution of the optimization problem provides an optimal interpolation model lying in the theoretical GRD function space. To achieve speed, we propose an efficient sampling algorithm that constructs a set of queries to maximize the expected information gain. The sampling scheme results in a unique sampling sequence invariant to source content, enabling parallel encoding and quality assessment processes. To achieve mathematical soundness, the GRD model is inherited from the Clough-Toucher (CT) interpolation method, and the function is differentiable everywhere on the domain of interests. Extensive experiments demonstrate that the resulting GRD function estimation framework achieves consistent improvement in speed and accuracy compared to the existing methods. The superiority and usefulness of the proposed algorithm are also evident by three applications.

II Related work

Although rate-distortion theory has been successfully employed in many multimedia applications, the research in the GRD surface modeling has only become a scientific field of study in the past decade. Existing methods can be roughly categorized based on their assumptions about the shape of a GRD function. The first model class only makes weak assumptions about the properties of the GRD functions. For example, [6] assumes the continuity of GRD functions and apply linear interpolation to estimate the response function after densely sampling the video representation space. However, the exhaustive search process is computationally expensive, not to mention the number of samples required increases exponentially with respect to the dimension of input space.

By contrast, the second class of models make strong a priori assumptions about the form of the GRD function to alleviate the need of excessive training samples. For example, [3] assumes the video quality exhibits an exponential relationship with respect to the quantization step, spatial resolution, and frame rate. Alternatively, Toni et al. [5, 13] derived a reciprocal function to model the GRD function. Similarly, [7] modeled the rate-quality curve at each spatial resolution with a logarithmic function. A significant limitation of these models is that domains of the analytic functional forms are restricted only to the bitrate dimension and several discrete resolutions, lacking the flexibility to incorporate other dimensions such as frame rate and bit depth, and the capability to predict the GRD behaviors at novel resolutions.

In addition to the specific limitations the two kinds of models may respectively have, they suffer from the same problem that the training samples in the GRD space are either manually picked or randomly selected, neglecting the difference in the informativeness of samples. While many recent works acknowledge the importance of GRD function [4, 5, 6, 7], a careful analysis and modeling of the response has yet to be done. We wish to address this void. In doing so, we seek a good compromise between 1) global and rigid models depending on random training samples and 2) local and indefinite models requiring exhaustive search in the video representation space.

III Modeling Generalized Rate-Distortion Functions

We begin by stating our assumptions. Our first assumption is that the GRD function is smooth. In theory, the Shannon lower bound, the infimum of the required bitrate to achieve a certain quality, is guaranteed to be continuous with respect to the target distortion [14]. On the other hand, successive change in the spatial resolution would gradually deviate the frequency component and entropy of the source signal, resulting in smooth transition in the perceived quality. In practice, many subjective experiments have empirically shown the smoothness of GRD functions [3, 15]. Furthermore, it is beneficial to impose $C^{1}$ continuity on GRD functions for better mathematical properties. For instance, the GRD surface is desired to be differentiable in many multimedia applications [5, 9].

Our second assumption is that the GRD function is axial-monotonic111In this work, we use rate-distortion function and rate-quality function interchangeably. Without loss of generality, we assume the function $f$ to be axial monotonically increasing. If $f$ is axial monotonically decreasing, we replace the given response with the function $f_{max}-f$, where $f_{max}$ is the maximum value of quality.. According to the rate-distortion theory [14], video quality increases monotonically with respect to the amount of resources it takes in the lossy compression. However, such monotonicity constraint may not apply to the spatial resolution. It has been demonstrated that encoding at high spatial resolution may even result in lower video quality than encoding at low spatial resolution under the same bitrate combined with upsampling and interpolation [6]. To be specific, encoding at high resolution with insufficient bitrate would produce artifacts such as blocking, ringing, and contouring, whereas encoding at low resolution with upsampling and interpolation would introduce blurring. The resulting distortions are further amplified or alleviated by the characteristics of the viewing device and viewing conditions, which interplay with HVS features such as the contrast sensitivity function [16]. A few sample GRD surfaces with their corresponding source videos are illustrated in Fig. 2.

Our third assumption is that the quality measurement is precise. Because the HVS is the ultimate receiver in most applications, subjective evaluation is a straightforward and reliable approach to evaluate the quality of digital videos. Traditional subjective experiment protocol models a subject’s perceived quality as a random variable, assuming the quality labeling process to be stochastic. Because subjective experiment is expensive and time consuming, it is hardly used in the GRD function approximation process. In practice, objective VQA methods that produce precise quality predictions are often employed to generate ground truth samples in the GRD function. Therefore, a GRD function should pass through the quality scores of objective VQA evaluated on the encoded video representations.

Under these assumptions, we define the space of GRD functions as:

[TABLE]

where $N$, $x_{n}$, $y_{n}$, and $z_{n}$ represent the total number of training samples, bitrate, spatial resolution, and quality of the $n$-th training sample, respectively.

In the subsequent sections, we introduce the proposed GRD model. Section III-A and III-C review the traditional CT method and the monotonicity condition of cubic polynomial Bézier function, which the proposed model relies on. The proposed $C^{1}$ continuity condition, optimization framework, and robust axial-monotonic CT algorithm are novel contributions that are detailed in Section III-B, III-D, and III-E, respectively.

III-A Review of Clough-Tocher Method

Since first introduced in 1960’s [17], the CT method has been the most widely used multi-dimensional scattered data interpolant, thanks to its $C^{1}$ continuity and low computational complexity [18, 19]. Consider the scattered points $(x_{n},y_{n})$ located in the $x,y$ plane and their values $z_{n}$ over the plane, the triangulation of the scattered points in the $x,y$ plane induces a piecewise triangular surface over the plane, whose nodes are the points $(x_{n},y_{n},z_{n})$. In CT method, a piecewise cubic function is employed as the interpolant for each triangle. Specifically, each triangle is further divided into three equivalent subtriangles, where a cubic function in the form of Bézier surface is estimated. Hereafter, we refer to the overall triangle as the macrotriangle and its subtriangles as microtriangles. The split triangular net is known as the control net, which is demonstrated in Fig. 3. For clarity and brevity, we also denote the macrotriangle edge that is opposite to the vertex $V_{i},i=0,1,2$ by $E_{i}$, and the internal microtriangle edge that connects $V_{i}$ and $S$ by $\hat{E}_{i}$. Mathematically, a cubic Bézier surface in the microtriangle $\Delta_{V_{0}V_{1}S}$ can be formulated as

[TABLE]

where $c_{V}$ represents the control net value at $V$ and $(\alpha,\beta,\gamma)$ is the barycentric coordinates with regard to the three vertices of the microtriangle. The barycentric coordinates of a point $P$ with regard to $\Delta_{V_{0}V_{1}S}$ can be defined as

[TABLE]

where $A_{UVW}$ means the directional area of the triangle formed by points $U,V,W$ and is positive when $U,V,W$ is counter-clockwise. The conversion from Cartesian coordinate to barycentric coordinate is lengthy and thus omitted here. Interested readers may refer to [18] for more details.

We note that 10 parameters are required to define a Bézier surface on each microtriangle. In the case of microtriangle $\Delta_{V_{0}V_{1}S}$, the parameters $c_{V_{0}}$, $c_{V_{1}}$, $c_{C_{2}}$ and $c_{S}$ are associated with point $V_{0}$, $V_{1}$, $C_{2}$ and $S$, respectively. The rest of the parameters are associated with the 6 trisection points on the three edges of $\Delta_{V_{0}V_{1}S}$. At the first glance, we need to determine 30 parameters for the interpolant in one triangle with only 3 equality constraints given by $c_{V_{0}}=z_{V_{0}}$, $c_{V_{1}}=z_{V_{1}}$, and $c_{V_{2}}=z_{V_{2}}$. Fortunately, the degree of freedom can be dramatically reduced with certain smoothness constraints. Under the $C^{0}$ assumption within macrotriangles, each two Bézier surfaces share the same curves at the their common boundaries $V_{0}S$, $V_{1}S$, and $V_{2}S$, leaving 19 free parameters in the macrotriangle $\Delta_{V_{0}V_{1}V_{2}}$. The inner-triangle $C^{1}$ continuity removes 7 additional degree of freedoms by enforcing the shaded neighboring microtriangles in Fig. 3 to be coplanar [20]. To ensure inter-triangle $C^{1}$ continuity, a standard approach is to assume the cross-boundary derivatives of the neighboring macrotriangles to be collinear, which further reduces the degree of freedom to 9. Taking into account the three known values at $V_{0}$, $V_{1}$, and $V_{2}$, we eventually have 6 unknown parameters in each macrotriangle. Although the gradients at vertices is not always available in practice, in most cases they can be estimated by considering the known values not only in the vertices of the triangle in question, but also in its neighbors. The most commonly used method is to estimate the gradients by minimizing the second-order derivatives along all Bézier curves [21]. Readers who are interested in the details of the CT method may refer to [21, 22, 18, 19].

The original CT method suffers from at least three limitations in approximating GRD functions. First, it picks the normal of the edges as the direction of cross-boundary derivative $d^{e}_{E_{i}}$. However, this choice gives an interpolant that is not invariant under affine transforms. This has some undesirable consequences: for a very narrow triangle, the spline can develop huge oscillations [22]. Second, the interpolant composite of piece-wise Bézier polynomials is not axial-monotonic, even when the given points are axial monotonic. Third, the CT algorithm achieves the external smoothness by estimating the gradients at three vertices $V_{i},i=0,1,2$, and by assuming the normal derivative at the triangle boundary $E_{i}$ to be linear. The linear assumption is somewhat arbitrary and may violate monotonicity we want to achieve. We will address the three limitations in the subsequent sections.

III-B Affine-Invariant C1 Continuity

In this section, we propose an affine invariant CT interpolant. Instead of the normal derivative at the triangle boundary $E_{i}$, we consider $d^{e}_{E_{i}}$ to be parallel to $\vec{C_{i}\bar{C}_{i}}$, i. e.

[TABLE]

where

[TABLE]

$x^{*}_{i}=\frac{(x_{\bar{C}_{i}}-x_{C_{i}})(x_{V_{j}}y_{V_{k}}-x_{V_{k}}y_{V_{j}})-(x_{V_{k}}-x_{V_{j}})(x_{C_{i}}y_{\bar{C}_{i}}-x_{\bar{C}_{i}}y_{C_{i}})}{(x_{\bar{C}_{i}}-x_{C_{i}})(y_{V_{k}}-y_{V_{j}})-(y_{\bar{C}_{i}}-y_{C_{i}})(x_{V_{k}}-x_{V_{j}})},$

$y^{*}_{i}=\frac{(y_{\bar{C}_{i}}-y_{C_{i}})(x_{V_{j}}y_{V_{k}}-x_{V_{k}}y_{V_{j}})-(y_{V_{k}}-y_{V_{j}})(x_{C_{i}}y_{\bar{C}_{i}}-x_{\bar{C}_{i}}y_{C_{i}})}{(x_{\bar{C}_{i}}-x_{C_{i}})(y_{V_{k}}-y_{V_{j}})-(y_{\bar{C}_{i}}-y_{C_{i}})(x_{V_{k}}-x_{V_{j}})},$

[TABLE]

[TABLE]

[TABLE]

$d^{x}_{V_{i}}$ and $d^{y}_{V_{i}}$ are partial derivatives of the Bézier surface at $V_{i}$ and $\{i,j,k\}$ is a cyclic permutation of $\{0,1,2\}$. Since this quantity transforms similarly as the gradient under affine transforms, the resulting interpolant is affine-invariant [22].

We also lift the unwanted linear constraints on the cross-boundary derivatives, elevating the number of parameters in a macrotriangle back to 9. In summary, the equality constraints in (2) can be factorized into the matrix form for simplicity

[TABLE]

where ${\bf c}\in\mathbb{R}^{16\times 1}$, ${\bf R}\in\mathbb{R}^{16\times 9}$, ${\bf d}\in\mathbb{R}^{9\times 1}$, ${\bf f}\in\mathbb{R}^{16\times 1}$, ${\bf c}$ and ${\bf d}$ represent the values of control net and unknown derivatives, respectively. Therefore, finding the interpolant of the macrotriangle corresponds to determining the 9 unknown parameters in ${\bf d}$.

Besides the inner macrotriangle constraints, we also want to keep $d^{e}_{E_{i}}$ consistent across the triangle boundary to ensure external $C^{1}$ smoothness. As a result, the following equality constraints need to be added for each edge with adjacent triangles

[TABLE]

Combining (3) and (4), we conclude that the resulting function is $C^{1}$ continuous and affine-invariant.

III-C Axial Monotonicity

This section aims to derive the sufficient constraints on ${\bf d}$ for the Bézier surface in the macrotriangle $\Delta_{V_{0}V_{1}V_{2}}$ to be axial-monotonic. In general, the interpolant composite of piece-wise Bézier polynomials is not monotonic even though the sampled points are monotonic. Several works have been done to derive sufficient conditions for a univariate or bivariate Bézier function [23, 24]. We adopt the sufficient condition proposed in [24], where it was proved that the cubic Bézier surface in a microtriangle is axial-monotonic when all the 6 triangular patches of its control net (e.g. $\Delta_{I_{02}I_{12}S}$, $\Delta_{C_{2}I_{11}I_{12}}$, $\Delta_{C_{2}I_{01}I_{02}}$, $\Delta_{V_{1}T_{10}I_{11}}$, $\Delta_{T_{10}T_{01}C_{2}}$, and $\Delta_{T_{01}V_{0}I_{01}}$ in $\Delta_{V_{0}V_{1}S}$) are axial-monotonic. By combining the sufficient conditions in all three microtriangles and the inner triangle continuity, we obtain

[TABLE]

We can summarize the monotonicity constraint in matrix form

[TABLE]

where ${\bf G}\in\mathbb{R}^{10\times 16}$ and ${\bf f}\in\mathbb{R}^{10\times 1}$. Further substitute (3) into (6), we obtain the monotonicity constraint in terms of ${\bf d}$

[TABLE]

More details on how we construct $\bf G$ and $\bf h$ are given in the Appendix.

III-D Optimization-based Solutions

To determine the unknown derivatives, we propose to minimize the total curvature of the interpolated surface under the smoothness assumption. Directly computing the total curvature is computationally intractable. Alternatively, we minimize the curvature of Bézier curves at the edges of each microtriangle as its approximation. Specifically, in $\Delta_{V_{0}V_{1}V_{2}}$, the objective function is written as

[TABLE]

where the weight $\frac{1}{2}$ is introduced to cancel the double counting of the external edges.

Consider an external boundary $E_{i}$, whose Bézier control net coefficients are $z_{V_{j}}$, $c_{T_{jk}}$, $c_{T_{kj}}$, and $z_{V_{k}}$. The integral of the second order derivative of the Bézier curve on $E_{i}$ can be represented in terms of the four coefficients as

[TABLE]

where

[TABLE]

is the length of $E_{i}$.

Similarly, we get the other part of the objective function from an internal boundary $\hat{E}_{i}$, whose coefficients are $z_{V_{i}},c_{I_{i1}},c_{I_{i2}},$ and $c_{S}$.

[TABLE]

where

[TABLE]

is the length of $\hat{E}_{i}$.

Substitute (III-D) and (III-D) into (8), we obtain the loss function for $\Delta_{V_{0}V_{1}V_{2}}$ in matrix form

[TABLE]

where ${\bf U}_{V_{0}V_{1}V_{2}}\in\mathbb{R}^{16\times 16}$ and ${\bf w}_{V_{0}V_{1}V_{2}}\in\mathbb{R}^{16\times 1}$.

Substituting ${\bf c}={\bf Rd}+{\bf f}$ into (33), we get

[TABLE]

In summary, finding the axial-monotonic interpolant corresponds to solving the following optimization problem

[TABLE]

Note that the constraints are linear with respect to ${\bf d}$ and ${\bf R}^{T}{\bf U}_{V_{0}V_{1}V_{2}}{\bf R}$ is positive-semidefinite. Thus, finding ${\bf d}$ turns into a standard problem of quadratic programming, which can be efficiently solved by the existing convex programming packages [25].

III-E Robust Axial-Monotonic Clough-Tocher Method

Here we propose our Robust Axial-Monotonic Clough-Tocher (RAMCT) method. The inequality constraints in (6) are sufficient conditions for $x$-axial monotonicity. However, the sufficient conditions cannot be satisfied in some extreme cases, making the primary solution infeasible. To relax these constraints, we introduce hinge loss to some of these inequalities, motivated by the success of Support Vector Machine [26]. Specifically, the modified inequality constraints are formulated as

[TABLE]

where $\boldsymbol{\xi}\in\mathbb{R}^{6\times 1}$ and $\boldsymbol{\xi}\leq{\bf 0}$. Note that (36a),(36c) are identical to (5a),(5c) because they are necessary conditions of axial monotonicity (See Appendix for proof). Rewriting these constraints in the matrix form, we obtain

[TABLE]

where ${\bf G}$ and ${\bf h}$ are the same as in (6),(7). ${\bf J}_{2}$ is a $6\times 6$ identity matrix, while ${\bf J}_{1}\in\mathbb{R}^{10\times 6}$ is obtained by padding ${\bf J}_{2}$ with 3 rows of zeros to its top and inserting a row of zeros between the 3rd and 4th rows of ${\bf J}_{2}$.

By substituting (3) into the inequality above, we finally obtain the inequality constraints in terms of the unknowns ${\bf d}$ and the auxiliary variables ${\boldsymbol{\xi}}$ as

[TABLE]

The objective function is then modified accordingly,

[TABLE]

where $\boldsymbol{\lambda}=[\lambda,\lambda,\cdots,\lambda]^{T}$ is the weighting parameter. Substituting ${\bf c}={\bf Rd}+{\bf f}$ into (51), we get

[TABLE]

Replacing (III-D),(6) with (III-E),(50) in (LABEL:eq:opt), we find that the original interpolation problem remains to be a quadratic programming problem.

IV Information-theoretic sampling

In this section, we first explore the informativeness of samples in the GRD space via a probabilistic model. We then present an information-theoretic sampling strategy that optimally selects the samples, offering enormous savings in time and computational resources.

Let ${\bf x}=(x_{1},...,x_{N})$ be a vector of discrete samples on a GRD function uniformly distributed in the bitrate-resolution space, where $N$ is the total number of sample points on the grid. Given that the GRD function is smooth, when the sampling grid is dense, these discrete samples provide a good description of the continuous GRD function. In particular, when the GRD function is band-limited, it can be fully recovered from these samples when the sampling density is larger than the Nyquist rate. Assuming x is created from GRD functions of real-world video content, we model x as an $N$-dimensional random variable, for which the probability density function $p_{\bf x}(\bf x)\sim\mathcal{N}(\boldsymbol{\mu},\,\boldsymbol{\Sigma})$ follows a multivariate Normal distribution. The total uncertainty of x is characterized by its joint entropy given by

[TABLE]

where $|\cdot|$ is the determinant operator. If the full vector x is further divided into two parts such that {\bf x}=\left[{\begin{array}[]{*{20}{c}}{\bf x}_{1}\\ {\bf x}_{2}\end{array}}\right] and $\boldsymbol{\Sigma}=\begin{bmatrix}\boldsymbol{\Sigma}_{11}&\boldsymbol{\Sigma}_{12}\\ \boldsymbol{\Sigma}_{21}&\boldsymbol{\Sigma}_{22}\end{bmatrix}$, and the ${\bf x}_{2}$ portion has been resolved by ${\bf x}_{2}={\bf a}$, then the remaining uncertainty is given by the conditional entropy

[TABLE]

where

[TABLE]

As a special case, we aim to find one sample that most efficiently reduces the uncertainty of GRD estimation. This is found by minimizing the log determinant of the conditional covariance matrix [27]

[TABLE]

where $\boldsymbol{\bar{\Sigma}}=\begin{bmatrix}\boldsymbol{\bar{\Sigma}}_{ii}&\boldsymbol{\bar{\sigma}}_{i}\\ \boldsymbol{\bar{\sigma}}_{i}^{T}&\bar{\sigma}_{ii}\end{bmatrix}$ and $i$ is the row index of $\boldsymbol{\bar{\Sigma}}$.

Minimizing (LABEL:eq:sampling) directly is computationally expensive, especially when the dimensionality is high. Alternatively, we minimize the upper bound of the conditional entropy

[TABLE]

where $\log|\boldsymbol{\bar{\Sigma}}_{ii}-\frac{\boldsymbol{\bar{\sigma}}_{i}^{T}\boldsymbol{\bar{\sigma}}_{i}}{\bar{\sigma}_{ii}}|\leq\mathrm{tr}(\boldsymbol{\bar{\Sigma}}_{ii}-\frac{\boldsymbol{\bar{\sigma}}_{i}^{T}\boldsymbol{\bar{\sigma}}_{i}}{\bar{\sigma}_{ii}}-{\bf I})$ and ${\bf I}$ denotes identity matrix. The sample with the minimum average loss in (LABEL:eq:samplingTrace) over all viewing devices is most informative. Once the optimal sample index is obtained, we encode the video at the $i$-th representation, evaluate its quality with objective VQA algorithms, and update the conditional covariance matrix in (64). The process is applied iteratively until the overall uncertainty in the system is reduced below a certain threshold $T$. We summarize the proposed uncertainty sampling method in Algorithm 1, where ${\bf r}_{i}$ represents the bitrate and spatial resolution at the $i$-th representation.

Remark: To get a sense of what type of samples will be chosen by the proposed algorithm, we analyze several influencing factors in the objective function (LABEL:eq:samplingTrace):

• •

By the basic properties of trace, the objective function in the uncertainty sampling can be factorized as

[TABLE]

Thus, $\mathrm{tr}(\boldsymbol{\bar{\Sigma}}_{ii}-\frac{\boldsymbol{\bar{\sigma}}_{i}^{T}\boldsymbol{\bar{\sigma}}_{i}}{\bar{\sigma}_{ii}})$ is a decreasing function with respect to $\bar{\sigma}_{ii}$ when $\bar{\sigma}_{ii}>\sqrt{\sum\limits_{j\neq i}\bar{\sigma}_{ij}^{2}}$. This indicates that samples with large uncertainty are more likely to be selected than those with small uncertainty.

• •

According to (64), $\forall j\neq i$,

[TABLE]

suggesting the rate of reduction in the uncertainty of sample $j$ is proportional to its squared correlation with the selected sample $i$ in the $k$-th iteration. Fig. 4 shows an empirical covariance matrix $\boldsymbol{\bar{\Sigma}}$ estimated from our video dataset that will be detailed in the next section, from which we observe that the GRD functions typically exhibit high correlation in a local region. Combining the first observation above, we conclude that the next optimal choice of sample should be selected from the region where labeled samples are sparse.

• •

Note that knowing that ${\bf x}_{2}={\bf a}$ alters the variance, though the new variance does not depend on the specific value of ${\bf a}$. The independence has two important consequences. First, the proposed sampling scheme is general enough to accommodate GRD estimators from all classes. More importantly, the algorithm results in a unique sampling sequence for all GRD functions. In other words, we can generate a lookup table of optimal querying order, making the sampling process fully parallelizable.

V Experiments

In this section, we first describe the experimental setups including our GRD function database, the implementation details of the proposed algorithm, and the evaluation criteria. We then compare the proposed algorithm with existing GRD estimation methods.

V-A Experimental setups

GRD Function Database: We construct a new video database which contains 250 pristine videos that span a great diversity of video content. An important consideration in selecting the videos is that they need to be representative of the videos we see in the daily life. Therefore, we resort to the Internet and elaborately select 200 keywords to search for creative common licensed videos. We initially obtain more than 700 4K videos. Many of these videos contain significant distortions, including heavy compression artifacts, noise, blur, and other distortions due to improper operations during video acquisition and sharing. To make sure that the videos are of pristine quality, we carefully inspect each of the videos multiple times by zooming in and remove those videos with visible distortions. We further reduce artifacts and other unwanted contaminations by downsampling the videos to a size of 1920 $\times$ 1080 pixels, from which we extract 10 seconds semantically coherent video clips. Eventually, we end up with 250 high quality videos.

Using the aforementioned sequences as the source, each video is distorted by the following processes sequentially:

• •

Spatial downsample: We downsample source videos using bi-cubic filter to six spatial resolutions (1920 $\times$ 1080, 1280 $\times$ 720, 720 $\times$ 480, 512 $\times$ 384, 384 $\times$ 288, 320 $\times$ 240) according to the list of Netflix certified devices [6].

• •

H.264/HEVC/VP9 compression: We encoded the downsampled sequences using three commonly used video encoders with two-pass encoding [1, 6, 13]. The target bitrate ranges from 100 kbps to 9 Mbps with a step size of 100 kbps.

In total, we obtain 540 (hypothetical reference circuit) $\times$ 250 (source) $\times$ 3 (encoder) = 405,000 test samples (currently the largest in the VQA community). We evaluate the quality of each video at a given spatial resolution, bitrate, and five commonly used display devices including cellphone, tablet, laptop, desktop, and TV using SSIMplus [28] for the following reasons. First, SSIMplus is currently the only HVS motivated spatial resolution and display device-adapted VQA model that is shown to outperform other state-of-the-art quality measures in terms of accuracy and speed [28, 29]. Second, a simplified VQA model SSIM [30] has been demonstrated to perform well in estimating the GRD functions [7]. The resulting dense samples of SSIMplus are regarded as the ground truth of GRD functions (The range of SSIMplus is from 0 to 100 with 100 indicating perfect quality). Our GRD modeling approach does not constrain itself on any specific VQA methods. When other ways of generating dense ground-truth samples are available, the same GRD modeling approach may also be applied.

Implementation Details: We initialize the scattered network with delaunay triangulation [31], inherited from CT method [17]. The balance weight $\lambda$ in (51) is set to $10^{-4}$. In our current experiments, the performance of the proposed RAMCT is fairly insensitive to variations of the value. We employed OSQP [25] to solve the quadratic programming problem, where the maximum number of iterations is set to $10^{6}$. The stopping criteria threshold $T$ is set to 540 (the total number of representation samples in the discretized GRD function space) $\times$ 10 (the standard deviation of mean opinion score in the LIVE Video Quality Assessment database), resulting in an average sample number of 38. When $tr(\Sigma)$ is below the threshold, we conclude that the uncertainty in the system can be explained by the disagreement between subjects. Therefore, further improvement in prediction accuracy may not be as meaningful. Since a triangulation only covers the convex hull of the scattered point set, extrapolation beyond the convex hull is not possible. In order to make a fair comparison, we initialize the training set $S$ as the representations with maximum and minimum bitrates at all spatial resolutions. To construct the covariance matrix described in Section IV as well as test the proposed algorithm, we randomly segregated the database into a training set of 200 GRD functions and a testing set with 50 GRD functions. The random split is repeated 50 times and the median performance is reported.

Evaluation Criteria: We test the performance of the GRD estimators in terms of both accuracy and rate of convergence. Specifically, we used two metrics to evaluate the accuracy. The mean squared error (MSE) and $l_{\infty}$ norm of the error values are computed between the estimated function and the actual function for each source content. The median results are then computed over all testing functions. All interpolation models can fit increasingly complex GRD functions at the cost of using many parameters. What distinguishes these models from each other is the rate and manner with which the quality of the approximation varies with the number of training samples.

V-B Performance

We test five GRD function models including reciprocal regression [5], logarithmic regression [7], 1D piecewise cubic Hermite interpolating polynomial (PCHIP), CT interpolation, and the proposed RAMCT on the aforementioned database. To evaluate the performance of the uncertainty sampling algorithm, we apply it on the five GRD models above and compare its performance with random sampling scheme as the baseline. For random sampling, the initial set of training sample $S$ is set as the representations with the maximum and minimum bitrates at all spatial resolutions to allow fair comparison. The training process with random sampling was repeated 50 times and the median performance is reported.

Table I and II show the prediction accuracy on the database, from which the key observations are summarized as follows. First, the models that assume a certain analytic functional form are consistently biased, failing to accurately fit GRD functions even with all samples probed. On the other hand, the existing interpolation models usually take more than 100 random samples to converge, although they are asymptotically unbiased. By contrast, the proposed RAMCT model converges with only a moderate number of samples. Second, we analyze the core contributors of RAMCT with deliberate selection of competing models. Specifically, the 1D monotonic interpolant PCHIP significantly outperforms the 2D generic interpolant CT, suggesting the importance of axial monotonicity. RAMCT achieves even better performance by exploiting the 2D structure and jointly modeling the GRD functions. Third, we observe strong generalizability of the proposed uncertainty sampling strategy evident by the significant improvement over random sampling on all models. The performance improvement is most salient on the proposed model. In general, RAMCT is able to accurately model GRD functions with only 30 labeled samples, based on which the reciprocal model merely have sufficient known variables to initialize fitting. To gain a concrete impression, we also recorded the execution time of the entire GRD estimation pipeline including video encoding, objective VQA, and GRD function approximation with the competing algorithms on a computer with 3.6GHz CPU and 16G RAM. RAMCT with uncertainty sampling takes around 10 minutes to reduce $l_{\infty}$ below 5, which is more than 100 times faster than the tradition regression models with random sampling.

VI Applications

The application scope of GRD model is much broader than VQA. Here we demonstrate three use cases of the proposed model.

VI-A Rate-Distortion Curve at Novel Resolutions

Given a set of rate-distortion curves at multiple resolutions, it is desirable to predict the rate-distortion performance at novel resolutions, especially when there exists a mismatch between the supported viewing device of downstream content delivery network and the recommended encoding profiles. Traditional methods linearly interpolate the rate-distortion curve at novel resolutions [6], neglecting the characteristics of GRD functions. Fig. 5 compares the linearly interpolated and RAMCT-interpolated rate-distortion curves at 960$\times$540 with the ground truth SSIMplus curve, from which we have several observations. First, the linearly interpolated curve shares the same intersection with the neighboring curves at 740$\times$480 and 1280$\times$720, inducing consistent bias to the prediction. The proposed RAMCT model is able to accurately predict the quality at the intersection of the neighboring curves by taking all known rate-distortion curves into consideration. Second, the linearly interpolated rate-distortion curve always lies between its neighboring curves, suggesting that the predicted quality at any bitrate is lower than the quality on one of its neighboring curves. This behavior contradicts the fact that each resolution may have a bitrate region in which it outperforms other resolutions [6]. On the contrary, RAMCT better preserves the general trend of resolution-quality curve at different bitrate, thanks to the regularization imposed by the $C^{1}$ condition at given nodes. Third, RAMCT outperforms the linear interpolation model in predicting the ground truth rate-distortion curve across all bitrates. The experimental results also justify the effectiveness of the $C^{1}$ and smoothness prior used in RAMCT.

To further validate the performance of the proposed GRD model at novel spatial resolutions, we predict the rate-distortion curves of 20 randomly selected source videos from the dataset at three novel resolutions (640$\times$360, 960$\times$540, and 1600$\times$900). The evaluated bitrate ranges from 100 kbps to 9 Mbps with a step size of 100 kbps. The results are listed in Table III. We can observe that RAMCT outperforms the linear model [6] with a clear margin at novel resolutions.

VI-B Per-Title Encoding Profile Generation

To overcome the heterogeneity in users’ network conditions and display devices, video service providers often encode videos at multiple bitrates and spatial resolutions. However, the selection of the encoding profiles are either hard-coded, resulting in sub-optimal QoE due to the negligence of the difference in source video complexities, or selected based on interactive objective measurement and subjective judgement that are inconsistent and time-consuming. To deliver the best quality video to consumers, each title should receive a unique bitrate ladder, tailored to its specific complexity characteristics. We introduce a quality-driven per-title optimization framework to automatically select the best encoding configurations, where the proposed GRD model serves as the key component.

Content delivery networks often aim to deliver videos at certain quality levels to satisfy different viewers. It is beneficial to minimize the bitrate usage in the encoding profile when achieving the objective. Mathematically, the quality-driven bitrate ladder selection problem can be formulated as a constrained optimization problem. Specifically, for the $i$-th representation,

[TABLE]

where $x$, $y$, $f(\cdot,\cdot)$, $C_{i}$ and $m$ represent the bitrate, the spatial resolution, the GRD function, the target quality level of video representation $i$, and the total number of video representations, respectively. Solving the optimization problem requires precise knowledge of the GRD function. Thanks to the effectiveness and differentiability of RAMCT, the proposed model can be incorporated with gradient-based optimization tools [32] to solve the per-title optimization problem. (Interested readers may refer to the appendix for more details on how we solve the optimization problem).

To validate the proposed per-title encoding profile selection algorithm, we apply the algorithm to generate bitrate ladders using H.264 [33] for 50 randomly selected videos in the aforementioned dataset. We set the target quality levels $\{C_{i}\}_{i=1}^{10}$ as $\{30,40,50,60,70,75,80,85,90,95\}$ to cover diverse quality range and to match the total number of representations in standard recommendations [5]. For simplicity, we optimize the representation sets for only one viewing device (cellphone), while the procedure can be readily extended to multiple devices to generate a more comprehensive representation set. In Fig. 6, we compare the rate-quality curve of representation sets generated by the proposed algorithm, recommendations by Netflix [34], Apple [35], and Microsoft [36] for three videos with different complexities, from which the key observations are as follows. First, contrasting the hand-crafted bitrate ladders, the encoding profile generated by the proposed algorithm is content adaptive. Specifically, the encoding bitrate increases with respect to the complexity of the source video as illustrated from Fig. 6(a) to Fig. 6(c). Second, the proposed method achieves the highest quality at all bitrate levels. The performance improvement is mainly introduced by the encoding strategy at the convex hull encompassing the individual per-resolution rate-distortion curves [6]. Table IV provides a full summary of the Bjøntegaard-Delta bitrate (BD-Rate) [37], indicating the required overhead in bitrate to achieve the same SSIMplus values. We observe that the proposed framework outperforms the existing hard-coded representation sets by at least 47%.

VI-C Codec Comparison

In the past decade, there has been a tremendous growth in video compression algorithms, thanks to the fast development of computational multimedia. With many video encoders at hand, it becomes pivotal to compare their performance, so as to find the best algorithm as well as directions for further advancement. Bjøntegaard-Delta model [37, 38] has become the most commonly used objective coding efficiency measurement. Bjøøntegaard-Delta PSNR (BD-PSNR) and BD-Rate are typically computed as the difference in bitrate and quality (measured in PSNR) based on the interpolated rate-distortion curves

[TABLE]

where $x_{A}$ and $x_{B}$ are the logarithmic-scale bitrate, $z_{A}$ and $z_{B}$ are the quality of the interpolated reference and test bitrate curves, respectively. $[x_{L},x_{H}]$ and $[z_{L},z_{H}]$ are the effective domain and range of the rate-distortion curves. However, BD-PSNR and BD-Rate do not take spatial resolution and viewing condition into consideration. Fig. 7 shows two GRD surfaces generated by H.264 [33] and HEVC encoders for a source video. Although H.264 performs on par with HEVC at low resolutions, it requires higher bitrate to achieve the same target quality at high resolutions. Therefore, applying BD-Rate on a single resolution is not sufficient to fairly compare the overall performance between encoders. To this regard, we propose generalized quality gain ($Q_{gain}$) and rate gain ($R_{gain}$) models as

[TABLE]

where $p(u)$, $U$, and $[y_{L},y_{H}]$ represent the probability density of viewing devices, the set of all device of interests, and the domain of video spatial resolution, respectively. The generalized $Q_{gain}$ and $R_{gain}$ models represent the expected quality gain and the expected bitrate gain (saving when $R_{gain}$ negative) across all spatial resolutions and viewing devices, leading to a more comprehensive evaluation of codecs. It should be noted that $z_{A}(x,y,u)$ is essentially the GRD function of codec $A$, which can be efficiently approximated by the proposed model. $x_{A}(z,y,u)$ can also be estimated numerically from the interpolated surface. (Interested readers may refer to the appendix for more details on how we compute $Q_{gain}$ and $R_{gain}$.) Therefore, RAMCT is a natural fit to the generalized $Q_{gain}$ and $R_{gain}$ models. The effect of any individual influencing factor can be obtained by taking the marginal expectation in the corresponding dimension, which is more robust than BD-PSNR and BD-Rate at a single resolution.

Using the proposed measures on the aforementioned dataset, we evaluate the performance of three video encoders, namely H.264, VP9, and HEVC. In order to simplify the expression, we set $p(u)$ to be uniform distribution for five display devices including cellphone, tablet, laptop, desktop, and TV in the paper. Table V shows the performance of VP9 and HEVC in terms of the proposed measures with H.264 as the reference codec, from which we can observe that VP9 outperforms HEVC by an average of 1.7 and 11% in $Q_{gain}$ and $R_{gain}$, respectively. The results of our objective codec evaluation are in general consistent with the recent video quality assessment results [39, 40].

VII Conclusions

GRD functions represent the critical link between multimedia resource and perceptual QoE. In this work, we proposed a learning framework to model the GRD function by exploiting the properties all GRD functions share and the information redundancy of training samples. The framework leads to an efficient algorithm that demonstrates state-of-the-art performance, which we believe arises from the RAMCT model for imposing axial monotonicity, the joint modeling of the multi-dimensional GRD function for exploiting its functional structure, and the information-theoretic sampling algorithm for improving the quality of training samples. Extensive experiments have shown that the algorithm is able to accurately model the function with a very small number of training samples. Furthermore, we demonstrate that the proposed GRD model plays a central role in a great variety of visual communication applications.

The current work can be extended in many ways. As a basis for future work, we note that the interpolant can be readily extended to higher dimensions [22], making it applicable to more general applications. For example, in the fields of machine learning [27] and data visualization [41], flexible monotonic interpolation can provide regularization and makes the model more interpretable. Another promising future direction is to develop models that can predict GRD functions without sampling the GRD space.