# Pre-demosaic Graph-based Light Field Image Compression

Yung-Hsuan Chao, Haoran Hong, Gene Cheung, and Antonio Ortega Yung-Hsuan Chao
was with University of Southern California. She is now with Qualcomm
Technologies, Inc., 5775 Morehouse Drive, San Diego, CA 92121-1714, USA
(e-mail: yunghsua@qti.qualcomm.com).Haoran Hong is with University of Southern
California, 3740 McClintock Ave., Los Angeles, CA 90089-2564, USA (e-mail:
haoranho@usc.edu).Gene Cheung is with York University, 4700 Keele Street,
Toronto, Canada M3J 1P3 (e-mail: genec@yorku.ca).Antonio Ortega is with
University of Southern California, 3740 McClintock Ave., Los Angeles, CA
90089-2564, USA (e-mail: antonio.ortega@sipi.usc.edu).

###### Abstract

An unfocused plenoptic light field (LF) camera places an array of microlenses
in front of an image sensor in order to separately capture different
directional rays arriving at an image pixel. Using a conventional Bayer
pattern, data captured at each pixel is a single color component (R, G or B).
The sensed data then undergoes demosaicking (interpolation of RGB components
per pixel) and conversion to an array of sub-aperture images (SAIs). In this
paper, we propose a new LF image coding scheme based on graph lifting
transform (GLT), where the acquired sensor data are coded in the original
captured form without pre-processing. Specifically, we directly map raw sensed
color data to the SAIs, resulting in sparsely distributed color pixels on 2D
grids, and perform demosaicking at the receiver after decoding. To exploit
spatial correlation among the sparse pixels, we propose a novel intra-
prediction scheme, where the prediction kernel is determined according to the
local gradient estimated from already coded neighboring pixel blocks. We then
connect the pixels by forming a graph, modeling the prediction residuals
statistically as a Gaussian Markov Random Field (GMRF). The optimal edge
weights are computed via a graph learning method using a set of training SAIs.
The residual data is encoded via low-complexity GLT. Experiments show that at
high PSNRs—important for archiving and instant storage scenarios— our method
outperformed significantly a conventional light field image coding scheme with
demosaicking followed by High Efficiency Video Coding (HEVC).

###### Index Terms:

Light field imaging, image compression, graph signal processing, intra-
prediction, lifting transform

##  I Introduction

Figure 1: Conceptual system of lenselet-based unfocused plenoptic 1.0 camera

Light Field (LF) imaging separately captures light rays arriving from
different directions at each pixel in an image. With acquired LF data, multi-
view rendering [1, 2] and re-focusing [3, 4] become possible post-capture.
However, captured LF data are very large in volume compared to a conventional
color image of the same resolution, and hence efficient compression of LF data
is important for storage and transmission.

(a)

(b)

Figure 2: (a) A raw lenselet image, Friends1, from EPFL dataset and (b) the
array of SAIs

In the last decade, many hardware designs have been developed for LF
acquisition, including multiple camera arrays, aperture cameras, and lenselet-
based plenoptic cameras. Among them, lenselet-based plenoptic cameras are the
most popular. There are two types of plenoptic cameras that have been made
commercially available: unfocused (or plenoptic 1.0), such as Lytro Illum [5],
and focused (or plenoptic 2.0), such as Raytrix R42 [6].

In an unfocused plenoptic camera, a microlens array is placed at a distance
equal to the focal length in front of a conventional photo sensor using Bayer
color filters, as shown in Fig. 1. The array produces a special hexagonal
pattern on the resulting raw image (lenselet image), as shown in Fig. 2(a),
where pixels within a hexagon correspond to light rays through a microlens.
The obtained raw image is a 2D image comprising a 2D hexagonal grid of micro-
images (MI). Since light fields are commonly represented and processed as 4D
functions (i.e., 2D arrays of SAIs) [7, 1, 3], the acquired lenselet image
typically undergoes pre-processing, consisting of demosaicking (pixelwise RGB
interpolation) and conversion to multiple SAIs, as shown in Fig. 2(b). Each
SAI is like a typical 2D image, gathering pixels from a specific light
direction.

There exist two types of redundancies within the 4D LF representation: i)
spatial redundancy among neighboring pixels in an SAI, i.e., intra-view
correlation, and ii) angular redundancy among SAIs in a neighborhood, i.e.,
inter-view correlation.

In [8, 9, 10], intra- and inter-view correlations are exploited using regular
HEVC intra-prediction modes, together with a new self similarity (SS) mode,
which is similar to Intra Block Copy (IBC) in screen content coding extension
of HEVC (HEVC-SCC). In [11, 12, 13], view synthesis methods are used to
recover the entire light field from a small subset of encoded views. In [14,
15, 16, 17, 18, 19, 20, 21], Pseudo Video Sequence (PVS)-based coding schemes
are proposed, where the SAIs are treated as frames in a video sequence, so
that intra- and inter-prediction tools in the video codec can be used. In [22,
23, 24], inter-view prediction using multiview extension of HEVC (MV-HEVC) is
proposed to address the loss of correlation among views caused by scanning in
PVS-based schemes. In [25], a hybrid data representation method is proposed to
simultaneously exploit intra- and inter-view correlations in SAIs and MIs.
Other state-of-the-art methods can be found in recent surveys [26, 27].

While exploiting inter-view correlations results in higher coding efficiency,
it also leads to more complex encoding (due to motion / disparity prediction)
and creates dependencies among coded SAIs, which is undesirable for random
access. In particular, in an archiving scenario, a user may desire to quickly
browse through viewpoint images, each of which can be synthesized in
acceptably high quality using only a small subset of SAIs [28]. Thus, speedy
extraction of this image subset from the LF data compressed in high quality is
important. Moreover, synthesizing viewpoint images with a small set of SAIs
also reduces memory requirement, as the number of images needed to be decoded
and buffered is reduced. Finally, we note that, for initial storage, most
standard digital cameras use a low complexity codec (JPEG) operating by
default at very high PSNR (see Section VI-F2). In analogy, in this paper, we
consider an intra-view-only approach (which leads to faster encoding and
better random access), operating at high rates / PSNR.

Another source of inefficiency in existing techniques is the redundancy
generated by color interpolation in the demosaicking process, which increases
the number of pixels to be encoded. Moreover, in methods utilizing existing
image / video codecs such as JPEG, H.264 and HEVC, input RGB image/videos are
typically converted to other representations with chroma downsampling, e.g.
4:2:0 YUV format, which introduces distortions due to integer rounding and
color sub-sampling.

To avoid the demosaicking process that generates redundancy before
compression, in our preliminary work [29], we proposed a new pre-processing
pipeline, where the original lenselet color pixels captured by the photo
sensor are mapped onto sparse locations in a series of SAIs without performing
demosaicking and color conversion prior to encoding. Unlike pre-demosaic
compression proposed for regular image and video in [30, 31, 32, 33], where
RGB pixel locations follow the regular Bayer filtering pattern, pixels in pre-
demosaic SAIs are irregularly placed, making it difficult to apply existing
standard image and video coding schemes, e.g., JPEG and HEVC. To resolve this
difficulty, we proposed placing pixels belonging to each SAI on a graph and
encoding them using a graph-based lifting transform [34]. Since a pixel only
connects to other pixels in the same SAI, pre-demosaicked pixels in different
SAIs can be encoded and decoded independently, which favors random access
efficiency.

In this paper, we develop a complete system based on the preliminary ideas
presented in [29]. Specifically, we introduce a novel intra-prediction
approach that can operate on the pre-demosaic image, achieving significant
improvements in coding efficiency. Additionally, we develop a systematic
approach to optimize edge weights for the graph to be used for our graph-based
lifting transform.

For intra-prediction, we propose a novel intra-prediction strategy designed
for sparsely distributed pixels in SAIs. Specifically, we first estimate the
local characteristics, e.g., object contour, of each block based on structure
tensor [35], using already encoded neighboring pixel blocks. The information
is then used to adjust the shape of kernels used for prediction. As a result,
the predictor will use more information from neighboring pixels on the same
side of a contour and less information from pixels across the contour,
preserving contour sharpness after prediction.

The graph weights, which reflect similarities between two connected pixels of
residual intensities, are optimized based on Gaussian Markov Random Field
(GMRF) modeling of the signal. The problem of learning efficient graph
structures, i.e., graph topology and weights, from training data in image and
video has been studied recently [36, 37, 38, 39]. However, most methods are
designed for regularly sampled data, e.g., pixels on a rectangular grid and
therefore, unsuitable for irregularly distributed pixels. In this work, we
consider a graph learning problem where graphs is derived from training
vectors with missing entries that are different from vector to vector.

After defining a graph, the graph signals, i.e., the residual pixel
intensities, are encoded using a graph-based lifting transform [34]. Compared
to the state-of-the-art HEVC-based coding, the results show significant
gain—around $5dB$ on average—at the high bitrate/PSNR range, which is
important for applications that demand high-quality encoding such as archives.
Compared to our preliminary work [29], we observe around $3dB$ gain at both
high and low bitrates.

The outline of the paper is as follows. In Section II, we present notations
used throughout this paper and review conventional approaches in lenselet-
based LF image compression. We provide an overview of our proposed coding
scheme in Section III. In Section IV and V, we describe in details our
proposed intra-prediction and transform coding algorithms tailored for
sparsely distributed pixels. Experiments and conclusions are presented in
Section VI and VII, respectively.

##  II Notations and Background

###  II-A Notations

Table I: List of notations ant their meanings

Notations | Meanings
---|---
$f(\cdot),K(\cdot)$ | functions
$f(i),K(i)$ | function value at location $i$
$\nabla{\bf f}$ | gradient vector of function $f(\cdot)$
$f_{i}$ |  $i_{th}$ element in vector ${\bf f}$
$F_{i,j}$ | element of matrix ${\bf F}$ at position $(i,j)$
${\bf f}_{\mathcal{S}}$ | subvector of ${\bf f}$ extracted from positions specified by set $\mathcal{S}$
${\bf F}_{\mathcal{S},\mathcal{P}}$ | submatrix of ${\bf F}$ extracted from positions specified by set $\mathcal{S}$ and $\mathcal{P}$
${\bf F}^{{\sf T}}$ | matrix transpose
${\bf F}^{-1}$ $|$ ${\bf F}^{\dagger}$ | inverse $|$ pseudo-inverse of ${\bf F}$
$x\sim\mathcal{N}(\hbox{\boldmath$\mu$},\hbox{\boldmath$\Sigma$})$ | multivariate Gaussian with mean $\mu$ and covariance $\Sigma$
${\bf Q}$ | inverse covariance (precision) matrix of multivariate Gaussian
$O(\cdot)$ | big O notation in complexity analysis
${\bf A}$ $|$ ${\bf L}$ | graph adjacency $|$ Laplacian matrix
${\bf D}$ $|$ ${\bf H}$ | graph degree $|$ self loop matrix

Throughout this paper, we use lowercase normal (e.g., $f$ and $\sigma$),
lowercase bold (e.g., ${\bf f}$ and $\sigma$) and uppercase bold letter (e.g.,
${\bf F}$ and $\Theta$) to denote scalars, vectors, and matrices,
respectively. The elements of a $2$-dimensional vector (e.g., ${\bf
v}=[v_{1},v_{2}]^{{\sf T}}$) in the Euclidean space represent the directions
along the horizontal and vertical axis respectively, a notation commonly used
for the image coordinate system. Unless otherwise stated, calligraphic capital
letters (e.g., $\mathcal{S}$) are used to represent sets. Other notations are
summarized in Table II-A.

Figure 3: A conventional pre-processing pipeline for light field image
compression is shown. The demosaicking and calibration processes are applied
before compression. Macro-pixels are represented by dashed lines, and macro-
pixel centers are denoted by red points.

###  II-B Conventional pre-processing pipeline before compression

Fig. 3 shows a block diagram of pre-processing steps before a conventional
coding scheme is applied on light field image, where the captured Bayer-
patterned lenselet image is converted into a 2D array of full-color SAIs. In
the figure, the conversion process is based on the method proposed by
Dansereau et al. [40, 41].

First, in order to generate full-color lenselet images, the missing color
components at each pixel are interpolated using nearby pixels containing the
target colors. The number of pixels will be increased threefold through this
process regardless of the demosaicking algorithm used.

Second, the color lenselet image is calibrated via rotation, translation and
scaling, which requires interpolation and increases the data size. A lenselet
image consists of multiple hexagonally arranged pixel patches called macro-
pixels. Each macro-pixel collects light for one image pixel arriving from
different directions. However, due to manufacturing defects, the arrangement
of macro-pixels is usually not aligned with the image coordinates. So, the
color lenselet image needs to be calibrated with the help of white images
[40], so that each macro-pixel center falls onto an integer pixel location,
and the arrangement of macro-pixels is aligned to the regular hexagonal grid.

Third, pixels that correspond to light rays coming from the same direction are
arranged into one SAI. After the arrangement, the pixel grid is corrected to
be square. The resampling step, which involves interpolation [40], increases
the data size again.

##  III Proposed pre-processing pipeline for graph based compression

Figure 4: Proposed pre-processing pipeline for the color filter array (CFA) of
lenselet image Figure 5: Since no interpolation, i.e., color interpolation in
demosaicking and scaling in calibration, is applied, some pixel locations are
empty in the SAIs.

To avoid the aforementioned redundancies, we propose a new pipeline to pre-
process LF data, shown in Fig. 4.

First, we postpone demosaicking and calibrate raw lenselet images via affine
transform to align all lenselet centers with pixel centers. Pixels that fall
onto non-integer locations after transformation will be moved to the nearest
integer positions.

Then, based on the relative locations within the macro-pixels on the
calibrated image, pixels are arranged onto multiple SAIs. Due to the placement
of macro-pixels, the pixels in the SAIs are placed hexagonally, which is
different from the conventional pipeline, where SAIs are re-sampled into
rectangular placement.

Note that the pixels around the boundaries of each macro-pixel, which tend to
have noisy values due to underexposure, are discarded in the pipeline
described in [40]. However, in the proposed pipeline, the boundary pixels are
kept. Our pipeline, which does not change the number of pixels nor the RGB
values in the raw data, is reversible. Thus, if lossless compression is
applied the original lenselet images can be recovered.

An example of sparsely distributed G components on one SAI is shown in Fig. 5.

The pattern of available pixels on each SAI depends on the inputed white image
that is selected based on the camera zoom and focus settings appropriate for
the lenselet picture to be decoded. Hence, the pattern may change as camera
settings change [41]. Additionally, the pattern may also change for different
types of macro-pixel misalignments and different calibration algorithms
adopted.

We stress that our scheme described in the next few sections does not rely on
a particular selection of algorithms for demosaicking and calibration. In
fact, many different techniques have been proposed in recent years for LF
image calibration [42, 43] and demosaicking [40, 42, 44]. Our coding scheme
can be easily adapted to these various approaches by using different
calibration transforms at the encoder for pixel rearrangement. The
demosaicking strategy at the decoder side can also be adjusted accordingly.

##  IV Proposed Intra-prediction scheme

Figure 6: Proposed Intra-prediction system for sparsely distributed pixels

After calibration and decomposition, each SAI containing sparsely distributed
pixels undergoes our proposed intra-prediction. Fig. 6 shows a flow chart of
our intra-prediction scheme. First, each SAI is divided into non-overlapped
$m\times m$ blocks for intra-prediction, which is the basic unit for later
transform coding. Second, we estimate gradients of all the pixels in each
reference block, using which we compute the structure tensor for the block.
Third, structure tensors of decoded neighboring blocks are used to estimate
gradient direction and strength of a target block to be encoded, from which we
can determine the parameters of an adaptive kernel [45] used for intra-
prediction. Hence, no overhead is required in our intra-prediction scheme, as
the choice of a prediction mode is solely based on information from previously
encoded blocks.

###  IV-A Gradient Estimation

In our work, the local gradient at each sparse pixel $i$, $\nabla{\bf f}(i)$,
is estimated via linear regression. Using Taylor series expansion, the pixel
value at a given point ${\bf x}$ can be expressed as

$f({\bf x})=f({\bf a})+\nabla{\bf f}({\bf a})\cdot({\bf x}-{\bf
a})+\mathcal{O}(\|{\bf x}-{\bf a}\|^{2})$ (1)

where $\nabla{\bf f}({\bf a})$ is the gradient vector at point ${\bf a}$. For
${\bf a}$ sufficiently close to ${\bf x}$, the pixel value $f({\bf x})$ can be
well approximated using only the first two terms

$f({\bf x})\approx f({\bf a})+\nabla{\bf f}({\bf a})\cdot({\bf x}-{\bf a})$
(2)

Based on the linear approximation, we estimate the gradient $\nabla{\bf
f}({\bf a})$ at point ${\bf a}$ by fitting a hyperplane which best satisfies
Equ. (2) for a number of pixels $\\{{\bf x}_{1},{\bf x}_{2},\cdots{\bf
x}_{k}\\}$ close to ${\bf a}$. The fitting can be represented as an
overdetermined system:

F $\displaystyle=\textbf{X}\cdot\nabla{\bf f}({\bf a})$ (3)
$\displaystyle=\begin{bmatrix}f({\bf x}_{1})-f({\bf a})\\\ f({\bf
x}_{2})-f({\bf a})\\\ \vdots\\\ f({\bf x}_{k})-f({\bf
a})\end{bmatrix}=\begin{bmatrix}({\bf x}_{1}-{\bf a})^{{\sf T}}\\\ ({\bf
x}_{2}-{\bf a})^{{\sf T}}\\\ \vdots\\\ ({\bf x}_{k}-{\bf a})^{{\sf
T}}\end{bmatrix}\cdot\nabla{\bf f}({\bf a})$ (4)

The optimal gradient $\nabla{\bf f}({\bf a})^{*}$ can be computed by solving a
least-square problem

$\nabla{\bf f}({\bf a})^{*}=\operatorname*{arg\,min}_{\nabla{\bf f}({\bf
a})}\|\textbf{F}-\textbf{X}\cdot\nabla{\bf f}({\bf a})\|_{2}^{2}$ (5)

which has a closed-form solution

$\nabla{\bf f}({\bf a})^{*}=(\textbf{X}^{{\sf
T}}\textbf{X})^{\dagger}\textbf{X}^{{\sf T}}\textbf{F}.$ (6)

For each pixel, we use four nearby pixels ($k=4$) for hyperplane fitting in
Equ. (4), namely, the two closest neighbors in horizontal and vertical
orientations respectively. Selecting close neighbors along different
directions alleviates the problem when all pixels used for fitting are
aligned, resulting in an ill-conditioned system.

With the obtained gradient of available pixels, the structure tensor of block
B is calculated as:

${\bf H}_{\text{B}}=\Sigma_{i\in\text{B}}\nabla{\bf f}(i)\nabla{\bf
f}(i)^{{\sf
T}}=\Sigma_{i\in\text{B}}\begin{bmatrix}d_{h}(i)^{2}&d_{h}(i)d_{v}(i)\\\
d_{v}(i)d_{h}(i)&d_{v}(i)^{2}\end{bmatrix}$ (7)

where $d_{h}(i)$ and $d_{v}(i)$ denote the gradients of pixel intensities at
pixel $i$ along horizontal and vertical directions, respectively.

After computing the structure tensor in Equ. (7) using the estimated
gradients, we apply eigen-decomposition on the $2\times 2$ matrix ${\bf
H}_{\text{B}}$:

$\displaystyle{\bf H}_{\text{B}}$ $\displaystyle=\begin{bmatrix}{\bf
e}_{1}&{\bf e}_{2}\end{bmatrix}\begin{bmatrix}\lambda_{1}&0\\\
0&\lambda_{2}\end{bmatrix}\begin{bmatrix}{{\bf e}_{1}}^{{\sf T}}\\\ {{\bf
e}_{2}}^{{\sf T}}\end{bmatrix}$ (8) $\displaystyle=\lambda_{1}{\bf e}_{1}{{\bf
e}_{1}}^{{\sf T}}+\lambda_{2}{\bf e}_{2}{{\bf e}_{2}}^{{\sf T}}$ (9)

The corresponding eigenvectors ${\bf e}_{1}$ and ${\bf e}_{2}$ along with
their eigenvalues $\lambda_{1}$ and $\lambda_{2}$, where
$\lambda_{1}<\lambda_{2}$, summarize the gradient distribution within block B.
Eigenvector ${\bf e}_{2}$ represents the direction maximally aligned with the
block gradient, while the orthogonal ${\bf e}_{1}$ represents the dominant
edge direction. The two eigenvalues $\lambda_{1}$ and $\lambda_{2}$ provide
information of gradient strength along directions ${\bf e}_{1}$ and ${\bf
e}_{2}$, which will be used in constructing adaptive kernels for intra-
prediction (Section IV-C).

###  IV-B Structure Tensor Estimation

The estimate of structure tensor ${\bf H}_{\text{I}}$ of input block I is
calculated as the weighted average of the structure tensors from its decoded
neighboring blocks $\\{\text{B}_{1},\text{B}_{2},\cdots\\}$ (See Fig. 7). The
estimate is written as

${\bf H}_{\text{I}}=\dfrac{1}{c}\Sigma_{i}w_{i}\cdot(\dfrac{1}{n_{i}}{\bf
H}_{\text{B}_{i}}),$ (10)

where $n_{i}$ denotes the number of available pixels in reference block
$\text{B}_{i}$, $w_{i}$ is the weight associated with block $\text{B}_{i}$,
and $c$ is the normalization constant $c=\Sigma_{i}w_{i}$.

Figure 7: 4 decoded reference blocks and the vectors indication their relative
locations to the input block I and edge direction estimated

Each weight value $w_{i}$ is a function of the cosine similarity between the
gradient direction in reference block $\text{B}_{i}$ and the relative
direction from the input block I to $\text{B}_{i}$. The design is based on the
observation that edges in natural images are continuous contours, and can be
approximated locally with straight lines. In other words, the orientation of
edges is mostly consistent locally. Hence, if the edge direction calculated on
one reference block $\text{B}_{i}$ is consistent with its relative location
from I, denoted by a vector ${\bf v}_{i}$ from I to $\text{B}_{i}$ as shown in
Fig. 7, then block I is more likely to have the same edge orientation. In
estimating the gradient in the input block, reference blocks with consistent
edge orientation are assigned larger weights.

In our algorithm, we consider the four decoded neighboring blocks shown in
Fig. 7 as references for structure tensor estimation for intra-prediction.
Denote by ${\bf v}_{i}$ the (unit norm) vector from I to reference block
$\text{B}_{i}$, ${\bf e}_{i}$ the edge direction computed from structure
tensor in block $\text{B}_{i}$, and $\theta_{i}$ the angle between ${\bf
v}_{i}$ and ${\bf e}_{i}$ (Fig. 7), calculated as

$\theta_{i}=\arccos\left(\frac{{\bf v}_{i}\cdot{\bf e}_{i}}{\|{\bf
v}_{i}\|_{2}}\right)\in[0,\pi],$ (11)

with

$\\{{\bf v}_{1},{\bf v}_{2},{\bf v}_{3},{\bf
v}_{4}\\}=\\{[\frac{-1}{\sqrt{2}},\frac{-1}{\sqrt{2}}]^{{\sf T}},[-1,0]^{{\sf
T}},[\frac{-1}{\sqrt{2}},\frac{1}{\sqrt{2}}]^{{\sf T}},[0,-1]^{{\sf T}}\\}.$

The weight value $w_{i}$ in Equ. (10) is defined as

$w_{i}=\begin{cases}\exp(-\frac{\theta_{i}}{\delta})&\text{if
}\theta<\frac{\pi}{2}\\\ \exp(-\frac{\pi-\theta_{i}}{\delta})&\text{if
}\theta>\frac{\pi}{2}\\\ \end{cases}$ (12)

where $\delta$ is a parameter adjustable by users. Larger weights are assigned
to blocks with smaller $\theta_{i}$, i.e. such that the edge direction ${\bf
e}_{i}$ and unit vector ${\bf v}_{i}$ are nearly consistent.

Note that the pixel intensities ${\bf f}$ considered in the gradient
calculation in Section IV-A and the structure tensor estimation described in
Equ. (10) are based on only G components, since the sampling rate of G
components is higher than for B and R components in the Bayer pattern. This
leads to a more accurate approximation in Equ. (2). The structure tensor and
kernel parameters obtained from G pixels (see next section) are applied to B
and R components given that there is a high correlation among the intensity
variations in different color channels.

Figure 8: Illustration of kernel size and shape in the smooth block (left) and
the block with strong edge (right)

###  IV-C Data-adaptive Kernel Regression

For intra-prediction, we apply an adaptive kernel regression [45] centered at
each pixel. In the case of zero-order estimation, a prediction of the pixel
intensity at ${\bf x}$ can be calculated as the weighted average of its
neighboring pixels ${\bf x}_{i}$, written as

$\tilde{f}({\bf x})=\frac{1}{z_{{\bf x}}}\Sigma_{{\bf
x}_{i}\in\mathcal{R}_{\bf x}}(K({\bf x}_{i}-{\bf x})\cdot f({\bf x}_{i}))$
(13)

where $\mathcal{R}_{\bf x}$ denotes the set of neighboring pixels of ${\bf x}$
in the decoded reference blocks, and $z_{{\bf x}}$ is the normalization
constant $z_{{\bf x}}=\Sigma_{{\bf x}_{i}\in\mathcal{R}_{\bf x}}K({\bf
x}_{i}-{\bf x})$. A common choice for the kernel function $K(\cdot)$ is the
Gaussian kernel. For the data-adaptive kernel regression used in this work,
the Gaussian kernel is adapted based on the edge orientation and strength
derived from the computed structure tensor as

$K({{\bf x}_{i}-{\bf x}})=\frac{\sqrt{{\hbox{det}}(\bf
C)}}{2\pi\sigma^{2}}\exp(-\frac{({\bf x}_{i}-{\bf x})^{{\sf T}}{\bf C}({\bf
x}_{i}-{\bf x})}{2\sigma^{2}}),$ (14)

where ${\bf C}$ is the covariance matrix based on the eigenvectors and
eigenvalues of structure tensor from Equ. (9), and can be decomposed as

C $\displaystyle=\phi\cdot\begin{bmatrix}{\bf e}_{1}&{\bf
e}_{2}\end{bmatrix}\Lambda\begin{bmatrix}{{\bf e}_{1}}^{{\sf T}}\\\ {{\bf
e}_{2}}^{{\sf T}}\end{bmatrix}$ (15) $\displaystyle\Lambda$
$\displaystyle=\begin{bmatrix}\frac{1}{\epsilon}&0\\\
0&\epsilon\end{bmatrix}\quad\text{,
}\epsilon=\frac{\lambda_{2}+p_{2}}{\lambda_{1}+p_{2}}$ (16)
$\displaystyle\phi$
$\displaystyle=\frac{1}{n}\sqrt{\lambda_{1}\lambda_{2}+p_{1}},$ (17)

and where $[{\bf e}_{1}\quad{\bf e}_{2}]$ rotates the coordinates of Gaussian
kernel along the dominant gradient and edge directions, and $\epsilon$ is the
ratio of two eigenvalues, representing the relative strength of gradient in
${\bf e}_{2}$ from the perpendicular direction ${\bf e}_{1}$. The kernel will
be elongated for blocks of strong edges, where $\lambda_{2}\gg\lambda_{1}$,
and will be near-circular for smooth block
($\lambda_{2}\approx\lambda_{1}\approx 0$), as illustrated in Fig. 8. $\phi$
determines the scaling of the kernel size, where $n$ is the number of
available pixels in the block and $p_{1}$ and $p_{2}$ are two positive scalars
used to ensure numerical stability.

##  V Proposed transform coding

After intra-prediction, the residuals are placed on the vertices of the graphs
optimized by the training set and then undergo graph-based lifting transform.

###  V-A Preliminaries for Graphs

A graph $\mathcal{G=(V,E)}$ contains a set of $n$ vertices
$\mathcal{V}=\\{(v_{i},h_{i})\,|\,1\leq i\leq n\\}$, and a set of $m$ links
$\mathcal{E}=\\{(e_{i,j},w_{i,j})\,|\,1\leq i,j\leq n\\}$, each connecting two
distinct nodes $v_{i}$ and $v_{j}$. The value $h_{i}\geq 0$ denotes the vertex
weight, also called self loop. For a weighted graph, $w_{i,j}$ is a non-
negative real-valued weight on $e_{i,j}$, which captures the similarity
between the two connected nodes. We can represent the graph connectivity
information using a $n\times n$ adjacency matrix ${\bf A}$, which has zero
diagonal elements and ${\bf A}_{i,j}=w_{i,j}$ for off-diagonal terms. A degree
matrix ${\bf D}$ is a $n\times n$ diagonal matrix with
$D_{i,i}=\sum_{j}A_{i,j}$. A self loop matrix ${\bf H}$ is a diagonal matrix
with elements $H_{i,i}=h_{i}$. A generalized Laplacian matrix is defined as
${\bf L}={\bf D}-{\bf A}+{\bf H}$. A combinatorial Laplacian matrix is a
special case of the general Laplacian matrix with no self loops.

###  V-B Graph Learning for sparsely distributed pixels

In this section, we construct an optimal graph that best describes inter-pixel
similarities for each mode given observed intra-predicted residual blocks in
LF images, assuming a known statistical model.

LF images are divided into a training set and a test set. Blocks of intra-
predicted residuals from the training data are classified, based on their
respective computed structure tensor, into $9$ modes, including $8$
directional modes:

${\cal M}=\\{{\cal
M}_{\theta}|\quad\theta=\frac{-3\pi}{8},\frac{-\pi}{4},\frac{-\pi}{8},0,\frac{\pi}{8},\frac{\pi}{4},\frac{3\pi}{8},\frac{\pi}{2}\\}$
(18)

and a DC mode, ${\cal M}_{\text{DC}}$, if there is no dominant edge direction.
Given the edge angle $\gamma$ of a training block (the angle between the
smallest eigenvector and the horizontal axis) and the eigenvalues of its
structure tensor ($\lambda_{1}$ and $\lambda_{2}$, where
$\lambda_{1}\leq\lambda_{2}$), the class of the associated residual block is
determined as

${\cal M}=\begin{cases}{\cal M}_{\text{DC}}&\text{if
}\frac{\lambda_{2}+p}{\lambda_{1}+p}<T\\\ {\cal M}_{\theta}&\text{if
}\gamma\in[\theta-\frac{\pi}{16},\theta+\frac{\pi}{16})\end{cases}.$ (19)

Training blocks in a class are assumed to be samples of the same statistical
model, for which the corresponding optimal graph is derived. The problem of
finding the optimal graph can be seen as that of finding an optimal graph
Laplacian matrix, which completely characterizes the graph. Hence a total of
$9$ $n\times n$ generalized graph Laplacian matrices ${\bf L}_{i,i=1,2,\cdots
9}$ are derived through learning from intra-predicted residuals in the
training set.

In the literature, multiple techniques have been proposed to find the optimal
graph Laplacian matrix in different contexts. In [37, 39, 38, 46, 47], a
statistical assumption was made on data, where each observation ${\bf
f}\in\mathbb{R}^{N}$ is modeled as a realization of a Gaussian Markov Random
Field (GMRF), i.e., ${\bf f}\sim\mathcal{N}(\hbox{\boldmath$\mu$},{\bf
Q}^{-1})$, with probability density function

$p({\bf f}|\hbox{\boldmath$\mu$},{\bf Q})=\frac{{\hbox{det}}({\bf
Q})^{\frac{1}{2}}}{(2\pi)^{\frac{N}{2}}}\exp\bigg{(}-\frac{1}{2}({\bf
f}-\hbox{\boldmath$\mu$})^{{\sf T}}{\bf Q}({\bf
f}-\hbox{\boldmath$\mu$})\bigg{)},$ (20)

where $\mu$ is the mean vector, and ${\bf Q}$ is the inverse covariance matrix
(precision matrix). The problem for learning the graph Laplacian matrix ${\bf
L}$ given the GMRF assumption is to find the maximum likelihood estimate of
${\bf Q}$, which defines the partial correlation between pairwise variables in
${\bf f}$:

$p(f_{i},f_{j}\,|\,{\bf
f}\setminus\\{f_{i},f_{j}\\})=-\frac{Q_{i,j}}{\sqrt{Q_{i,i}Q_{j,j}}}$ (21)

In [39, 47], Egilmez et al. proposed an algorithm specifically to compute a
precision matrix ${\bf Q}$ targeting a graph Laplacian structure with non-
negative edge weights, i.e.,

$\begin{cases}Q_{i,j}<0\quad\text{if }A_{i,j}>0\\\ Q_{i,j}=0\quad\text{if
}A_{i,j}=0\\\ \end{cases}$ (22)

with additional constraints on the graph connectivity. Given $l$ i.i.d.
observations $\\{{\bf f}_{1},{\bf f}_{2},\cdots{\bf f}_{l}\\}$ of zero mean
GMRF, the precision matrix can be derived by solving the maximum likelihood
(ML) problem:

$\displaystyle\operatorname*{arg\,max}_{{\bf Q}}\prod_{i=1:l}p({\bf
f}_{i}|{\bf Q},\hbox{\boldmath$\mu$}_{i}=\textbf{0})$ (23) $\displaystyle=$
$\displaystyle\operatorname*{arg\,min}_{w,v}\quad\log{\hbox{det}}({\bf
Q})-\operatorname{Tr}({\bf Q}{\bf S})$ $\displaystyle\text{subject to }{\bf
Q}={\bf B}\operatorname{diag}({\bf w}){\bf B}^{{\sf
T}}+\operatorname{diag}({\bf h})$

where ${\bf S}$ is the sample covariance matrix, and ${\bf B}$ is the $n\times
m$ incident matrix, specifying the connectivity between node pairs. The
$m\times 1$ vector ${\bf w}$ contains the weights associated with each link,
and the $n\times 1$ vector ${\bf v}$ are the self loops on each node.

The above graph Laplacian learning algorithms are based on the assumption that
all the entries of an observation, ${\bf f}^{i}$, are available. However, in
an SAI without demosaicking, each pixel position contains at most one color
component out of R, G and B. The observed block for each color is represented
as a columnized vector

${\bf f}^{i}=\begin{bmatrix}{\bf f}^{i}_{\mathcal{O}}\\\ {\bf
f}^{i}_{\mathcal{M}}\end{bmatrix},$ (24)

where ${\bf f}^{i}_{\mathcal{O}}$ is a $r$-dimensional vector containing the
observed pixel intensities, and ${\bf f}^{i}_{\mathcal{M}}$ of dimension $n-r$
represents the missing data. $r$ is a variable that changes between blocks. In
order to optimize the graph connecting available pixels, we assume ${\bf
f}^{i}_{\mathcal{O}}$ to be a sub-sample of ${\bf f}^{i}$ modeled as GMRF.

In the statistics literature, many methods are proposed in estimating the
inverse covariance matrix in a GMRF based on sub-sampled observations [48, 49,
50]. In this work, we adopt the plug-in algorithm by Kolar and Xing [50] due
to its simplicity. The method consists of two steps. First, the sample
covariance matrix is estimated using the incomplete observations. Given
$f^{i}_{a}$ is the $a$-th variable in the observation ${\bf f}^{i}$, we define
${\bf r}^{i}$ as its indication vector, where element

$\begin{cases}r^{i}_{a}=1\quad\text{if }f^{i}_{a}\text{ is available}\\\
r^{i}_{a}=0\quad\text{otherwise},\end{cases}$ (25)

the estimated sample covariance matrix $\tilde{{\bf S}}$ is calculated as

$\tilde{S}_{a,b}=\frac{\sum_{i=1}^{p}r^{i}_{a}\cdot
r^{i}_{b}\cdot(f^{i}_{a}-\mu^{i}_{a})(f^{i}_{b}-\mu^{i}_{b})}{\sum_{i=1:p}r^{i}_{a}\cdot
r^{i}_{b}}$ (26)

See [50] for details and a theoretical justification of the estimation. In the
second step, $\tilde{{\bf S}}$ will be plugged into the objective function of
the maximum likelihood estimation of the precision matrix ${\bf Q}$. In this
work, we optimize ${\bf Q}$ based on Equ. (23) with ${\bf S}$ replaced by
$\tilde{{\bf S}}$.

Then, the precision matrix ${\bf Q}$, and therefore the graph Laplacian
matrix, can be derived for the whole block, including the pixel positions with
missing data. We can decompose precision matrix ${\bf Q}$ of the GMRF, which
is the inverse of covariance matrix $\hbox{\boldmath$\Sigma$}^{-1}$, as

$\displaystyle{\bf Q}$ $\displaystyle=\begin{bmatrix}{\bf Q}_{O,O}&{\bf
Q}_{O,M}\\\ {\bf Q}_{M,O}&{\bf Q}_{M,M}\end{bmatrix},$ (27)

based on the subset of observed ($O$) and missing data ($M$).

The inverse covariance matrix of available pixels, which defines the proximity
between pixels in ${\bf f}^{i}_{\cal O}$, can be calculated by taking the
Schur complement of the sub-matrix ${\bf Q}_{M,M}$ of ${\bf Q}$:

${\cal L}=\Sigma_{O,O}^{-1}={\bf Q}_{O,O}-{\bf Q}_{O,M}{\bf Q}_{M,M}^{-1}{\bf
Q}_{M,O}.$ (28)

In Fig. 9, two illustrative examples are shown for weighted graph optimized
for set $\mathcal{M}_{\theta=\frac{3\pi}{8}}$ and $\mathcal{M}_{\theta=0}$.
The figures in the left column are the regular graphs of the two classes
learnt from the training set. The figures in the right column are examples of
graphs connecting available sparsely distributed pixels, derived using the
Schur complement and sparsification. The nodes with strong self-loop weights
concentrate around the boundaries close to the reference blocks, which have a
better prediction. Therefore, the associated residuals on those pixels tend to
have lower variance.

(a) $\theta=\dfrac{3\pi}{8}$

(b) $\theta=0$ (horizontal)

Figure 9: Graph structure optimized for classes corresponding to intra-
prediction angle (1) $\theta=\dfrac{3\pi}{8}$ (nearly vertical) and (2)
$\theta=0$ (horizontal). The link color indicates the associated weight and
the node color indicates the associated self loop weight (darker: larger
weight)

###  V-C Graph-based Lifting Transform

Due to the large data size of LF images, we apply the localized graph-based
lifting transform (GLT), which was first proposed in [34]. The GLT is a multi-
resolution, critically sampled filterbank. It consists of three building
blocks: bipartition, prediction, and update. In bipartition, nodes are divided
into two disjoint sets ${\cal U}$ and ${\cal P}$. Then, graph signal ${\bf
f}_{\cal U}$ attached to nodes in ${\cal U}$ is used to predict signal ${\bf
f}_{\cal P}$ in ${\cal P}$. The resulting prediction residuals are stored in
${\cal P}$ as the high frequency transformed coefficients ${\bf
c}_{\text{H}}$. In the update step, ${\bf c}_{\text{H}}$ is used to update
${\bf f}_{\cal U}$, which leads to a low frequency approximation ${\bf
c}_{\text{L}}$ of the input signal. The process can be written in matrix form
as follows.

Prediction: $\displaystyle\begin{bmatrix}{\bf c}_{\text{H}}\\\ {\bf f}_{\cal
U}\end{bmatrix}=\begin{bmatrix}{\bf I}&-{\bf P}\\\ {\bf 0}&{\bf
I}\end{bmatrix}\begin{bmatrix}{\bf f}_{\cal P}\\\ {\bf f}_{\cal
U}\end{bmatrix}$ (29) Update: $\displaystyle\begin{bmatrix}{\bf
c}_{\text{H}}\\\ {\bf c}_{\text{L}}\end{bmatrix}=\begin{bmatrix}{\bf I}&{\bf
0}\\\ {\bf U}&{\bf I}\end{bmatrix}\begin{bmatrix}{\bf c}_{\text{H}}\\\ {\bf
f}_{\cal U}\end{bmatrix},$ (30)

where ${\bf P}$ and ${\bf U}$ denote the linear operations for prediction and
update, respectively. The transform is guaranteed to be invertible regardless
of how the nodes are bi-partitioned and how ${\bf P}$ and ${\bf U}$ are
designed.

In the context of compression, the main objective in designing the lifting
scheme is to steer most of the data energy to the low frequency band ${\cal
U}$, namely to reduce the energy of prediction residuals ${\bf
c}_{\text{H}}={\bf f}_{\cal P}-{\bf P}{\bf f}_{\cal U}$ in its complement
${\cal P}$. For signal ${\bf f}$ modeled as GMRF ${\cal
N}(\hbox{\boldmath$\mu$},{\bf Q}^{-1})$, the minimum mean square estimator of
signal value $f_{i}$ on node $v_{i}\in{\cal P}$ is expressed as

$\displaystyle\tilde{f}_{i}$
$\displaystyle=\mu_{i}-\sum_{j}\frac{Q_{i,j}}{Q_{i,i}}(f_{j}-\mu_{j})$ (31)
$\displaystyle=-\sum_{j}\frac{Q_{i,j}}{Q_{i,i}}f_{j}+(\mu_{i}+\sum_{j}\frac{Q_{i,j}}{Q_{i,i}}\mu_{j}).$

If the estimated precision matrix satisfies the graph Laplacian structure, as
described in Equ. (22) in the previous subsection, and if
$\hbox{\boldmath$\mu$}_{i}\approx\hbox{\boldmath$\mu$}_{j}=\mu$, then the
above equation can be simplified to

$\displaystyle\tilde{f}_{i}$
$\displaystyle=-\sum_{j}\frac{Q_{i,j}}{Q_{i,i}}f_{j}+\mu(1+\sum_{j}\frac{Q_{i,j}}{Q_{i,i}})$
(32)
$\displaystyle=\frac{1}{D_{i,i}+H_{i,i}}\sum_{j}A_{i,j}f_{j}+\frac{H_{i,i}}{D_{i,i}+H_{i,i}}\mu.$

The estimator is simply a weighted average of the signal values on the
neighboring nodes of $v_{i}$ and its associated mean. The self-loop weight
${\bf H}_{i,i}$ can be interpreted as a measurement of similarity to the mean
$\mu$. In our experiment, we assume $\mu$ to be $0$ based on the observation
that SAIs after intra-prediction mostly have low magnitude. Note that if the
graph is bipartite, i.e. every node in ${\cal U}$ only has connections to
nodes in ${\cal P}$, then the proposed filterbanks are optimal in terms of
mean square error for GMRF. The prediction operation in Equ. (32), when self
loop weight ${\bf H}_{i,i}$ of each node is equal to $0$, is equivalent to the
low complexity CDF53 filterbanks

$\tilde{f}_{i}=\frac{1}{D_{i,i}}\sum_{v_{j}\in{\cal U}}A_{i,j}f_{j}.$ (33)

For bipartition, we apply the algorithm proposed in [51] using maximum
spanning tree to assign the available pixels to update and prediction set with
linear time complexity. However, when graphs are not bipartite, the
information of proximity on the links connecting nodes in the same set cannot
be utilized, which leads to a suboptimal solution for prediction. Therefore,
we apply the re-connection algorithm proposed in [52] to transform them into
bipartite ones. The algorithm re-connects each node $v_{i}\in{\cal P}$ to be
predicted to nodes in the opposite set ${\cal U}$ using Kron reduction [53].
Then, a simple sparsification is applied, which keeps only $k$ ($k=4$ in the
experiment) links with the largest weights for each node in ${\cal P}$, in
order to achieve localization and low complexity in prediction. The
generalized CDF5/3 predictor in Equ. (32) is applied on the newly formed
bipartite graph. It can be proven that without the sparsification, the applied
predictor is equivalent to the minimum mean square error estimator (MMSE) for
${\bf f}_{\cal P}$ given ${\bf f}_{\cal U}$. Moreover, with the help of Kron
reduction, which can be implemented iteratively, we can achieve lower
complexity in implementation than directly computing the MMSE estimator, which
requires computing matrix inverse.

Same as the work in [52], the design of the update filter is based on [54],
which promotes the orthogonality of the filterbanks, and the construction of
graphs for lifting transform at levels $2$ and higher is based on Kron
reduction. The transformed coefficients are uniformly quantized and reordered
based on the approach in [55]. For entropy coding, we apply the Amplitude and
Group Partitioning (AGP) method proposed by Said and Pearlman in [56].

##  VI Experiments

###  VI-A Dataset

The LF images used in our experiments were from the EPFL Light-field data set
in the JPEG Pleno Database [57]. The raw data were captured with a Lytro Illum
camera [58]. The LF data used for training and testing are listed in Table
VI-A. Four SAIs selected from each of the LF images in the training set were
used for graph learning. Each test image was of size $5368\times 7728$. In the
baseline scheme, the raw data was converted into $15\times 15$ full-color
SAIs. Each SAI was of size $434\times 625$. We cropped the SAIs into
$432\times 624$, so that the dimensions were multiples of the minimum coding
unit of baseline schemes.

Table II: LF data for training and testing

Training set | Test set
---|---
Ankylosaurus_&_Stegosaurus | Friends_1
Ceiling_Light | Bikes
ISO_Chart_16 | Flowers
Perforated_Metal_1 | Ankylosaurur & Diplodocus 1
Sophie_&_Vincent_3 |
Yan_&_Krios_standing |

(a)

(b)

(c)

(d)

Figure 10: The central view of the test images (a) Friends1 (b) Bikes, (c)
Flowers, and (d) Ankylosaurur & Diplodocus 1

###  VI-B Baseline Schemes

We compare our proposed method with state-of-the-art coding solutions. The
following anchor coding solutions were used as benchmarks, namely:

  1. 1.

High Efficiency Video Coding all intra-prediction (HEVC-AI) / random access
(HEVC-RA) / screen content coding extension (HEVC-SCC): HM-16.22 [59]
developed by JCT-VC.

  2. 2.

Versatile Video Coding all intra-prediction (VVC-AI) / random access (VVC-RA):
VTM-10.0 [60] developed by JVET.

For fair comparison against our scheme, post filters, including deblocking
filter, sample adaptive offset (SAO), and adaptive loop filter (ALF), were all
disabled in the baseline.

###  VI-C Experimental Setting

####  VI-C1 First series of experiments

We consider three different scenarios for the graph-based coding scheme:

  1. 1.

DGLT: geometric distance based graph lifting transform without intra-
prediction, proposed in [29].

  2. 2.

intraDGLT: geometric distance based graph lifting transform with the proposed
intra-prediction (Section IV).

  3. 3.

intraLGLT: graph learning based graph lifting transform (Section V-B) with the
proposed intra-prediction.

Since our scheme is a pre-demosaic coding technique, as a pre-demosaic coding
baseline we chose HEVC encoding of the pre-demosaic lenselet images in 4:0:0
RGB. For extensive evaluation, we also encoded demosaicked lenselet images
with HEVC-AI and VVC-AI in 4:4:4 RGB. Since the IBC technique in HEVC-SCC can
further improve coding gain over HEVC, the HEVC-SCC encoding of the full-color
lenselet images in 4:4:4 RGB was also compared.

####  VI-C2 Second series of experiments

$15\times 15$ SAIs were encoded as a pseudo video sequence created using a
serpentine scanning order. The order is defined in JPEG Pleno Light Field
Coding Common Test Condition [61]. In addition to exploiting both intra- and
inter-view redundancy using HEVC-RA and VVC-RA configurations, since we are
focusing on coding schemes that allow efficient random access, the baselines
coded with All-Intra configuration were also included, where every SAI was
coded as a single 2D image with HEVC-AI and VVC-AI in original 4:4:4 RGB and
4:2:0 YUV formats. In this part, only the best graph-based coding scheme,
intraLGLT, was used to compare against the above-mentioned anchor coding
solutions.

###  VI-D Evaluation Metrics

The evaluation metrics were the number of bits per pixel (bpp) and average
peak signal-to-noise ratio (PSNR) over 4:4:4 RGB format of reconstructed SAIs.

For archival purposes, it is more appropriate to evaluate the quality of
reconstruction in RGB color space and on the reconstructed lenselet images,
i.e., all captured pixels should be reconstructed. However, current state-of-
the-art schemes using HEVC discard under-exposed pixels at the boundary of
macro-pixels during the conversion to SAI array. Moreover, the interpolation
involved in the resampling process is irreversible, so that we cannot recover
the original lenselet image from the SAI. Hence, for evaluation, performances
were compared against the full-color SAIs before compression (ground truth),
generated from the raw lenselet data using the demosaicking and calibration
pipeline described in [40, 41].

For the 4:2:0 YUV format, the reconstructed SAIs were translated back to 4:4:4
RGB format before evaluation. The up-sampling for U and V components was based
on nearest neighbor interpolation. For HEVC/VVC encoding of the demosaicked
lenselet image in Section VI-C1, the reconstructed lenselet image was
decomposed into SAIs for comparison.

###  VI-E Parameter Setting

Proposed. Each SAI was divided into non-overlapped $8\times 8$ blocks. The
parameter $\delta$ in Equ. (12) was set to $0.9$, and the parameters
$\sigma,p_{1},p_{2}$ in Equ. (14) and (16) for data-adaptive Gaussian kernal
were chosen to be $1.6$, $0.001$ and $0.001$, respectively. The threshold
value $T$ in Equ. (19) for training block classification was set to $1.5$. We
applied $2$-level lifting transform for transform coding for pixels in each
block. The selected QP values ranged from $4$ to $36$.

Anchor. The profile, bit depth, GoP and QP settings for the anchor coding
schemes in the experiments are summarized in Table III.

Profile: main_444_RExt |  GoP | QP (start:stride:end)
---|---|---
Bit depth: 8-bit |  |  YUV420 |  RGB444
The 1st series of experiments |  |  |
HEVC-SCC All Intra-prediction (lenselet) |  1 |  |  11:3:35
HEVC All Intra-prediction (predemosaic lenselet) |  1 |  |  2:3:35
HEVC All Intra-prediction (lenselet) |  1 |  |  11:3:35
VVC All Intra-prediction (lenselet) |  1 |  |  11:3:35
The 2nd series of experiments |  |  |
HEVC All Intra-prediction |  1 |  5:3:29 |  14:3:35
VVC All Intra-prediction |  1 |  |  14:3:35
HEVC Random Access |  8 |  |  5:3:29
HEVC Random Access |  16 |  2:3:26 |  5:3:29
VVC Random Access |  16 |  |  5:3:29
Table III: Summary of parameter settings of anchor coding solutions in the
experiments

(a)

(b)

(c)

(d)

Figure 11: PSNR-RGB performance in the first series of experiments: (a)
Friends1 (b) Bikes, (c) Flowers, and (d) Ankylosaurur & Diplodocus

(a)

(b)

(c)

(d)

Figure 12: PSNR-RGB performance in the second series of experiments (a)
Friends1 (b) Bikes, (c) Flowers, and (d) Ankylosaurur & Diplodocus

###  VI-F Experimental Results and Analysis

####  VI-F1 The first series of experiments

Fig. 11 shows the PSNR comparisons for test images in Table VI-A.

For our scheme, with the proposed intra-prediction, the correlation between
neighboring blocks were exploited to reduce redundancies before compression,
leading to 3dB performance gain over the setup where no intra-prediction was
performed [29]. Moreover, with graph learning, local directional
characteristics were exploited to determine graph connections, which provided
better accuracy in modeling similarity among pixels compared to connections
based solely on Euclidean distance. We observe that, using graph learning,
around $0.5$dB gain was achieved against a graph design based solely on
distance when intra-prediction was used.

Compared with other anchor schemes, the HEVC encoding of the pre-demosaic
lenselet images has significant coding gain at high bitrates. However, at
lower rates, its performance dropped more drastically compared with the HEVC
encoding of demosaicked lenselet, which will be discussed in more detail in
Section. VI-F3. With more advanced intra-prediction techniques, VVC-AI
outperformed HEVC-AI at any rate with $1dB$ performance gain. HEVC-SCC
outperformed HEVC by using additional tools, e.g., IBC, to exploit similarity
between neighboring MIs in prediction.

For applications such as archival and instant storage on cameras, images are
stored in very high quality. Therefore, in the evaluation, we focused mainly
on the high rate region. For bpp $>$ 3, graph-based coding schemes
outperformed baseline schemes applying HEVC/VVC in demosaicked lenselet image.
Applying HEVC/VVC intra coding on the whole lenselet image provided
significant gains at low bit rate over the proposed scheme, which can be
caused by the lack of accuracy in the proposed structure tensor based intra-
prediction. The estimation of structure tensor in the target block depends on
the gradients of the reconstructed reference blocks. When reference blocks are
coded at low bit rate, the gradients may fail to reflect the characteristics,
e.g., edge direction, in the target block. Although HEVC/VVC intra coding on
the whole lenselet image outperforms the proposed scheme at low bit rate, it
introduces large dependencies among SAIs, which results in a less efficient
random access.

####  VI-F2 The second series of experiments

Fig. 12 shows that, for higher bitrate (bpp $>$ 2), graph-based coding schemes
significantly outperformed the baseline schemes applying HEVC-AI and VVC-AI to
demosaicked SAIs. HEVC-RA and VVC-RA had similar performance and outperformed
HEVC/VVC-AI, because there exist only small disparity, and therefore high
correlation, between SAIs caused by light rays coming from different
directions. Further, the anchor coding schemes of GoP equal to $16$
outperformed those of GoP equal to 8 by using more inter-prediction between
frames while sacrificing random access efficiency. Although the proposed LGLT
scheme performed worse than the anchors at low bitrates, it achieved much
higher coding efficiency at high bitrates, making it attractive for archival
applications.

To substantiate the claim that, for the archiving purpose, images are stored
at high bitrates, we encoded SAIs of the four test images with the default
JPEG setting that the iPhone applies on captured images, and measured their
individual bitrates. The archival bitrate range measured is summarized in
Table IV. As observed in Fig. 12, in the archival bitrate range identified,
our scheme outperformed all the All-Intra coding schemes, with around 3dB
performance gain over VVC-AI and 5dB performance gain over HEVC-AI. Compared
to the baseline schemes using inter-prediction, our scheme has the advantage
of random access, which is analyzed in Section VI-G. Next, we provide insights
into why our method has advantages at very high rates.

Test images | Bitrate (bpp)
---|---
Friends_1 | 2.05 - 2.24
Bikes | 2.45 - 2.66
Flowers | 2.71 - 2.93
Ankylosaurur & Diplodocus 1 | 1.45 - 1.61

Table IV:  The bitrates of SAIs encoded by JPEG with default setting from the
iPhone

####  VI-F3  Advantages of proposed method at high bitrates

In a conventional scheme, the dimension of raw images to be encoded is
increased to obtain full-color images via interpolation and demosaicking. This
is a minor problem at low rates, since the interpolated images are smooth and
can be interpolated with few transform coefficients. However, as the bitrate
increases, higher frequency coefficients are no longer quantized to zero, and
encoding more coefficients increases the rate substantially. Instead, in our
proposed scheme, demosaicking and interpolation are postponed, so that the
total number of coefficients to be coded is smaller, resulting in
substantially better performance at high rates, which also explains why the
HEVC encoding of pre-demosaic lenselet images can significantly outperform the
other anchors in the high-rate region.

Moreover, because compression is performed prior to interpolation in our
proposed method, compression noise in encoded pixels is spread to interpolated
pixels. Thus, at low rates, error in encoded pixels, which is high, is
propagated via interpolation, resulting in worse performance than conventional
methods. Instead, at high rates, since the number of pixels to be encoded is
reduced, per pixel distortion is significantly lower than the conventional
method, and the effect of error propagation via interpolation is negligible.
Thus, at high rates, the pre-demosaic coding technique can outperform the
conventional techniques applied to demosaicked images.

In comparison to the HEVC coding of pre-demosaic lenslet images, the graph
structure provides better modeling of pixel similarities. Specifically, a
graph connects a pixel only to its adjacent pixels of the same color, avoiding
the discontinuities between different color components, and thus achieves
better coding efficiency at high bit rates against coding the Bayor patterned
pre-demosaic lenslet images.

|  HEVC based | Proposed scheme
---|---|---
|  |  w/o Parallelization |  Parallelization
Intra-prediction |  O($N$) |  O($N^{2}$) |  O($N$)
Transform |  O($N\log N$) |  O($N$) |  O($N$)
Table V:  Time complexity analysis of HEVC based coding scheme and the
proposed graph based coding scheme

###  VI-G Complexity Analysis

1\. Time Complexity Analysis

Time complexity of our proposed coding scheme is summarized in Table V for the
cases with and without parallelization, compared to the complexity in HEVC
based coding scheme. In the HEVC/VVC coding scheme, intra-prediction of each
pixel in the prediction unit is computed as the weighted average of up to two
reference samples in the block boundary, and thus time complexity
$\textit{O}(N)$ is needed to predict N data points. The DCT2 transform is
applied horizontally/vertically, leading to time complexity $\textit{O}(N\log
N)$. Note that in HEVC/VVC, the intra-prediction mode and the size of a coding
unit are determined by recursively splitting a coding tree unit and comparing
rate-distortion cost, which potentially increases the complexity in comparison
to the proposed intra-prediction in the graph-based methods.

In our proposed intra-prediction, the data-adaptive steering kernel is used to
calculate weights for reference pixels. In the worst case, one pixel in the
predicted block is predicted with the weighted average of reference pixels in
the four neighboring blocks, resulting in $\textit{O}(N^{2})$ to predict $N$
data points. If the weight of each reference sample is calculated in parallel,
the time complexity for intra-prediction can be reduced to $\textit{O}(N)$.
The graph-based lifting transform is applied on the intra-predicted residuals.
Note that the pattern of available pixels on each SAI changes depending on the
camera settings, and hence different graphs need to be constructed for
different cameras. Each pixel (node) is connected to its $4$-nearest neighbors
in our experiment, which requires $\textit{O}(N)$ in computation. The
prediction and update filters—CDF 5/3 filterbanks—used in the proposed scheme
are localized, i.e., computation of the transform coefficient of each pixel
involves only its connected neighbors on the bipartite graph. In addition,
only $2$-level lifting transform is applied. Therefore, for $N$ data points,
the complexity of transform in our proposed coding scheme is $\textit{O}(N)$.
As a result, comparable decoding complexity can be achieved at decoder for our
proposed coding scheme compared to HEVC based coding scheme if parallelization
is applied.

2\. Random Access Complexity Analysis

We define random access complexity as the number of reference frames needed
for decoding to access a chosen view. We consider the average and the maximum
number of reference frames required to be decoded for a view as the average
and worst-case random access complexity respectively. We evaluate the random
access complexity of the coding schemes in the experiments by using the
following criteria [23]:

$C_{\text{average}}^{\text{RA}}=\frac{1}{N}\sum_{x=1}^{N}f(x)$ (34)
$C_{\text{worst-case}}^{\text{RA}}=\max_{x=1...N}{f(x)}$ (35)

where $N$ is the total number of views, and $f(x)$ is a function returning the
number of SAIs required to be decoded for accessing the view $x$. The results
of random access complexity are summarized in Table VI.

Method |  Worst-case |  Average
---|---|---
HEVC/VVC-AI |  $0$ |  $0$
HEVC/VVC-RA (GoP=8) |  $4$ |  $2.98$
HEVC/VVC-RA (GoP=16) |  $5$ |  $3.98$
Proposed scheme |  $0$ |  $0$

Table VI: Random access complexity analysis of anchor coding solutions in RA
and AI coding scheme and the proposed graph based coding scheme

Although applying HEVC/VVC intra coding on the whole lenselet image and inter
coding on the PVS-based SAIs provides significant gains at low bitrates over
schemes coding each SAI separately, it also introduces large dependencies
among SAIs. That is, in order to access a specific SAI, the user has to decode
the whole lenselet image, which greatly increases latency and memory buffer
requirements. In schemes where the SAIs are coded separately, on the other
hand, an SAI can be encoded/decoded independently from each other.

For evaluation purpose, we employed the demosaicking algorithm in [62] for
both test and anchors. Note that while a specific demosaicking algorithm is
employed in the experiment, the proposed scheme is targeting coding the pre-
demosaic pixels, and is not explicitly tied to one particular demosaicking
algorithm. The algorithm in the experiment requires a few neighboring pixels
in up to $10$ reference SAIs to demosaick a single position. In the
literature, however, there exist many demosaicking methods that perform
interpolation using far fewer neighboring color components, as described in
[63], resulting in fast implementations and less dependency in our LF coding
scheme. Tailoring for our coding scheme, a demosaicking algorithm that
interpolates a missing color component using only available color pixels in
the same SAI, such as the methods in [42] and [44], can also be chosen to
avoid introducing any dependency among SAIs on the decoder side, though they
incur some marginal quality degradation.

Since intra-prediction is also a local operation, the demosaicking algorithm
in the experiment can be implemented as a small moving window, e.g., $24\times
24$ block, of local blocks in the reference SAIs that includes both the
predictor and the target block containing reference pixels to a demosaicked
pixel. Although the process in the experiment entails some dependencies
between an SAI and its neighbors, since pixels in each neighboring SAI can be
decoded in parallel, the latency requirement is significantly lower than
coding the whole lenselet image or the latency in the PVS-based HEVC scheme
[14, 15, 16, 17, 18, 19]. In the worst case, 4 SAIs have to be decoded
sequentially before accessing an SAI if the GoP is equal to 8, and 5 SAIs are
required if the GoP is equal to 16.

##  VII Conclusion

In this paper, we describe a novel coding scheme for light field image based
on the graph based lifting transform. Our scheme can encode the original raw
data without introducing redundancies (due to demosaicking and calibration) or
distortion (from color conversion and downsampling). Moreover, we propose an
intra-prediction scheme and a graph learning algorithm for pixels in SAIs that
are sparsely distributed. The pixels are then connected as graphs and encoded
with low-complexity graph-based lifting transform. Coding results in the high
rate region demonstrate that our proposed method outperforms the widely
applied HEVC-based and the cutting-edge VVC-based approaches in All-Intra and
random access configuration.

## References

  * [1] M. Levoy and P. Hanrahan, “Light field rendering,” in Proceedings of the 23rd annual conference on Computer graphics and interactive techniques, pp. 31–42, ACM, 1996.
  * [2] M. Magnor and B. Girod, “Data compression for light-field rendering,” IEEE Transactions on Circuits and Systems for Video Technology, vol. 10, no. 3, pp. 338–343, 2000.
  * [3] R. Ng, M. Levoy, M. Brédif, G. Duval, M. Horowitz, and P. Hanrahan, “Light field photography with a hand-held plenoptic camera,” Computer Science Technical Report CSTR, vol. 2, no. 11, pp. 1–11, 2005.
  * [4] M. Levoy, B. Chen, V. Vaish, M. Horowitz, I. McDowall, and M. Bolas, “Synthetic aperture confocal imaging,” in ACM Transactions on Graphics (ToG), vol. 23, pp. 825–834, 2004.
  * [5] “Lytro illum.” https://illum.lytro.com/.
  * [6] “Raytrix camera.” https://www.raytrix.de/.
  * [7] S. J. Gortler, R. Grzeszczuk, R. Szeliski, and M. F. Cohen, “The lumigraph,” in Proceedings of the 23rd annual conference on Computer graphics and interactive techniques, pp. 43–54, ACM, 1996.
  * [8] C. Conti, P. Nunes, and L. D. Soares, “New HEVC prediction modes for 3d holoscopic video coding,” in IEEE International Conference on Image Processing (ICIP), pp. 1325–1328, 2012.
  * [9] C. Conti, P. Nunes, and L. D. Soares, “HEVC-based light field image coding with bi-predicted self-similarity compensation,” in IEEE International Conference on Multimedia & Expo Workshops (ICMEW), 2016, pp. 1–4, 2016.
  * [10] R. Monteiro, L. Lucas, C. Conti, P. Nunes, N. Rodrigues, S. Faria, C. Pagliari, E. da Silva, and L. Soares, “Light field HEVC-based image coding using locally linear embedding and self-similarity compensated prediction,” in IEEE International Conference on Multimedia & Expo Workshops (ICMEW), pp. 1–4, 2016.
  * [11] B. Hériard-Dubreuil, I. Viola, and T. Ebrahimi, “Light field compression using translation-assisted view estimation,” in 2019 Picture Coding Symposium (PCS), pp. 1–5, 2019.
  * [12] I. Viola, H. P. Maretic, P. Frossard, and T. Ebrahimi, “A graph learning approach for light field image compression,” in Applications of Digital Image Processing XLI (A. G. Tescher, ed.), vol. 10752, pp. 126 – 137, International Society for Optics and Photonics, SPIE, 2018.
  * [13] P. Astola and I. Tabus, “Wasp: Hierarchical warping, merging, and sparse prediction for light field image compression,” in 2018 7th European Workshop on Visual Information Processing (EUVIP), pp. 1–6, 2018.
  * [14] F. Dai, J. Zhang, Y. Ma, and Y. Zhang, “Lenselet image compression scheme based on subaperture images streaming,” in IEEE International Conference on Image Processing (ICIP), pp. 4733–4737, 2015.
  * [15] D. Liu, L. Wang, L. Li, Z. Xiong, F. Wu, and W. Zeng, “Pseudo-sequence-based light field image compression,” in International Conf. on Multimedia & Expo Workshops (ICMEW) IEEE, pp. 1–4, 2016.
  * [16] S. Zhao, Z. Chen, K. Yang, and H. Huang, “Light field image coding with hybrid scan order,” in Visual Communications and Image Processing (VCIP), 2016, pp. 1–4, 2016.
  * [17] H. P. Hariharan, T. Lange, and T. Herfet, “Low complexity light field compression based on pseudo-temporal circular sequencing,” in Broadband Multimedia Systems and Broadcasting (BMSB), 2017 IEEE International Symposium on, pp. 1–5, 2017.
  * [18] A. Vieira, H. Duarte, C. Perra, L. Tavora, and P. Assuncao, “Data formats for high efficiency coding of lytro-illum light fields,” in Image Processing Theory, Tools and Applications (IPTA), 2015 International Conference on, pp. 494–497, 2015.
  * [19] L. Li, Z. Li, B. Li, D. Liu, and H. Li, “Pseudo sequence based 2-d hierarchical coding structure for light-field image compression,” in Data Compression Conference (DCC), 2017, pp. 131–140, 2017.
  * [20] C. Jia, Y. Yang, X. Zhang, X. Zhang, S. Wang, S. Wang, and S. Ma, “Optimized inter-view prediction based light field image compression with adaptive reconstruction,” in 2017 IEEE International Conference on Image Processing (ICIP), pp. 4572–4576, 2017.
  * [21] R. J. Monteiro, N. M. Rodrigues, S. M. Faria, and P. J. Nunes, “Light field image coding with flexible viewpoint scalability and random access,” Signal Processing: Image Communication, vol. 94, p. 116202, 2021.
  * [22] G. Wang, W. Xiang, M. Pickering, and C. W. Chen, “Light field multi-view video coding with two-directional parallel inter-view prediction,” IEEE Transactions on Image Processing, vol. 25, no. 11, pp. 5104–5117, 2016.
  * [23] N. Mehajabin, S. R. Luo, H. Wei Yu, J. Khoury, J. Kaur, and M. T. Pourazad, “An efficient random access light field video compression utilizing diagonal inter-view prediction,” in 2019 IEEE International Conference on Image Processing (ICIP), pp. 3567–3570, 2019.
  * [24] J. Khoury, M. T. Pourazad, and P. Nasiopoulos, “A new prediction structure for efficient mv-hevc based light field video compression,” in 2019 International Conference on Computing, Networking and Communications (ICNC), pp. 588–591, 2019.
  * [25] R. J. S. Monteiro, N. M. M. Rodrigues, S. M. M. Faria, and P. J. L. Nunes, “Light field image coding based on hybrid data representation,” IEEE Access, vol. 8, pp. 115728–115744, 2020.
  * [26] C. Conti, L. D. Soares, and P. Nunes, “Dense light field coding: A survey,” IEEE Access, vol. 8, pp. 49244–49284, 2020.
  * [27] C. Brites, J. Ascenso, and F. Pereira, “Lenslet light field image coding: Classifying, reviewing and evaluating,” IEEE Transactions on Circuits and Systems for Video Technology, vol. 31, no. 1, pp. 339–354, 2021.
  * [28] H. Piexoto, Macchiavello and Cheung, “Progressive communication for interactive light field data streaming,” in IEEE International Conference on Image Processing, 2017.
  * [29] Y.-H. Chao, G. Cheung, and A. Ortega, “Pre-demosaic light field image compression using graph lifting transform,” in Image Processing (ICIP), 2017 IEEE International Conference on, pp. 3240–3244, 2017.
  * [30] S.-Y. Lee and A. Ortega, “A novel approach of image compression in digital cameras with a Bayer color filter array,” in IEEE International Conference on Image Processing, vol. 3, pp. 482–485, 2001.
  * [31] C. C. Koh, J. Mukherjee, and S. K. Mitra, “New efficient methods of image compression in digital cameras with color filter array,” IEEE Trans. on Consumer Electronics, vol. 49, no. 4, pp. 1448–1456, 2003.
  * [32] K.-H. Chung and Y.-H. Chan, “A lossless compression scheme for Bayer color filter array images,” IEEE Transactions on Image Processing, vol. 17, no. 2, pp. 134–144, 2008.
  * [33] S.-Y. Lee and A. Ortega, “A novel approach for compression of images captured using Bayer color filter arrays,” arXiv preprint arXiv:0903.2272, 2009\.
  * [34] S. K. Narang and A. Ortega, “Lifting based wavelet transforms on graphs,” in Proceedings: APSIPA ASC 2009: Asia-Pacific Signal and Information Processing Association, 2009 Annual Summit and Conference, pp. 441–444, 2009\.
  * [35] J. Bigun and G. H. Granlun, “Optimal orientation detection of linear symmetry,” in In First International Conference on Computer Vision, ICCV (London), pp. 433–438, IEEE Computer Society Press, 1987.
  * [36] W. Hu, G. Cheung, and A. Ortega, “Intra-prediction and generalized graph fourier transform for image coding,” IEEE Signal Processing Letters, vol. 22, no. 11, pp. 1913–1917, 2015.
  * [37] E. Pavez, H. E. Egilmez, Y. Wang, and A. Ortega, “GTT: Graph template transforms with applications to image coding,” pp. 199–203, 2015.
  * [38] E. Pavez and A. Ortega, “Generalized laplacian precision matrix estimation for graph signal processing,” in Acoustics, Speech and Signal Processing (ICASSP), 2016 IEEE International Conference on, pp. 6350–6354, 2016.
  * [39] H. E. Egilmez, Y.-H. Chao, A. Ortega, B. Lee, and S. Yea, “GBST: Separable transforms based on line graphs for predictive video coding,” pp. 2375–2379, 2016.
  * [40] D. G. Dansereau, O. Pizarro, and S. B. Williams, “Decoding, calibration and rectification for lenselet-based plenoptic cameras,” in Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pp. 1027–1034, 2013.
  * [41] D. G. Dansereau, “Light field toolbox.” https://www.mathworks.com/matlabcentral/fileexchange/49683-light-field-toolbox-v0-4.
  * [42] S. Xu, Z.-L. Zhou, and N. Devaney, “Multi-view image restoration from plenoptic raw images,” in Asian Conference on Computer Vision, pp. 3–15, Springer, 2014.
  * [43] Y. Bok, H.-G. Jeon, and I. S. Kweon, “Geometric calibration of micro-lens-based light field cameras using line features,” IEEE Trans. on Pattern Analysis and Machine Intelligence, vol. 39, no. 2, pp. 287–300, 2017\.
  * [44] M. Seifi, N. Sabater, V. Drazic, and P. Perez, “Disparity-guided demosaicking of light field images,” in IEEE International Conference on Image Processing (ICIP),, pp. 5482–5486, 2014.
  * [45] H. Takeda, S. Farsiu, and P. Milanfar, “Kernel regression for image processing and reconstruction,” IEEE Transactions on image processing, vol. 16, no. 2, pp. 349–366, 2007.
  * [46] X. Dong, D. Thanou, P. Frossard, and P. Vandergheynst, “Learning laplacian matrix in smooth graph signal representations,” vol. 64, pp. 6160–6173, 2016\.
  * [47] H. E. Egilmez, E. Pavez, and A. Ortega, “Graph learning from data under laplacian and structural constraints,” vol. 11, pp. 825–841, 2017.
  * [48] N. Städler and P. Bühlmann, “Missing values: sparse inverse covariance estimation and an extension to sparse regression,” Statistics and Computing, vol. 22, no. 1, pp. 219–235, 2012.
  * [49] K. Lounici et al., “High-dimensional covariance matrix estimation with missing observations,” Bernoulli, vol. 20, no. 3, pp. 1029–1058, 2014\.
  * [50] M. Kolar and E. P. Xing, “Consistent covariance selection from data with missing values,” in Proceedings of the 29th International Conference on Machine Learning (ICML-12), pp. 551–558, 2012.
  * [51] H. Q. Nguyen and M. N. Do, “Downsampling of signals on graphs via maximum spanning trees.,” IEEE Trans. Signal Processing, vol. 63, no. 1, pp. 182–191, 2015.
  * [52] Y.-H. Chao, A. Ortega, and S. Yea, “Graph-based lifting transform for intra-predicted video coding,” in 2016 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pp. 1140–1144, 2016.
  * [53] F. Dorfler and F. Bullo, “Kron reduction of graphs with applications to electrical networks,” IEEE Transactions on Circuits and Systems I: Regular Papers, vol. 60, no. 1, pp. 150–163, 2013.
  * [54] G. Shen and A. Ortega, “Tree-based wavelets for image coding: Orthogonalization and tree selection,” in Picture Coding Symposium, 2009\. PCS 2009, pp. 1–4, 2009.
  * [55] E. Martínez-Enríquez and A. Ortega, “Lifting transforms on graphs for video coding,” in 2011 Data Compression Conference, pp. 73–82, 2011\.
  * [56] A. Said and W. A. Pearlman, “Low-complexity waveform coding via alphabet and sample-set partitioning,” in Electronic Imaging’97, pp. 25–37, International Society for Optics and Photonics, 1997.
  * [57] “Light-field image dataset.” http://plenodb.jpeg.org/.
  * [58] M. Řeřábek and T. Ebrahimi, “New light field image dataset,” in 8th International Conference on Quality of Multimedia Experience (QoMEX), no. EPFL-CONF-218363, 2016.
  * [59] “HEVC HM reference software’.”  https://hevc.hhi.fraunhofer.de/.
  * [60] “VVC VTM reference software”.”  https://jvet.hhi.fraunhofer.de/.
  * [61] F. P. et al, “Jpeg pleno light field coding common test conditions v3.2,” Doc. ISO/IEC JTC 1/SC 29/WG1 N83029, Geneva, Switzerland, Mar. 2019.
  * [62] H. S. Malvar, L.-w. He, and R. Cutler, “High-quality linear interpolation for demosaicing of bayer-patterned color images,” in Acoustics, Speech, and Signal Processing, 2004. Proceedings.(ICASSP’04). IEEE International Conference on, vol. 3, pp. iii–485, 2004.
  * [63] Y. Niu, J. Ouyang, W. Zuo, and F. Wang, “Low cost edge sensing for high quality demosaicking,” IEEE Transactions on Image Processing, vol. 28, no. 5, pp. 2415–2427, 2018.