# Lagrange–Chebyshev Interpolation for image resizing

Donatella Occorsio<EMAIL_ADDRESS>Giuliana Ramella
<EMAIL_ADDRESS>Woula Themistoclakis<EMAIL_ADDRESS>Department of Mathematics and Computer Science, University of Basilicata,
Viale dell’Ateneo Lucano 10, 85100 Potenza, Italy  C.N.R. National Research
Council of Italy, Institute for Applied Computing “Mauro Picone”, Via P.
Castellino, 111, 80131 Naples, Italy

###### Abstract

Image resizing is a basic tool in image processing and in literature we have
many methods, based on different approaches, which are often specialized in
only upscaling or downscaling. In this paper, independently of the (reduced or
enhanced) size we aim to get, we approach the problem at a continuous scale
where the underlying continuous image is globally approximated by the tensor
product Lagrange polynomial interpolating at a suitable grid of first kind
Chebyshev zeros. This is a well–known approximation tool that is widely used
in many applicative fields, due to the optimal behavior of the related
Lebesgue constants. Here we show how Lagrange–Chebyshev interpolation can be
fruitfully applied also for resizing an arbitrary digital image in both
downscaling and upscaling. The performance of the proposed method has been
tested in terms of the standard SSIM and PSNR metrics. The results indicate
that, in upscaling, it is almost comparable with the classical Bicubic
resizing method with slightly better metrics, but in downscaling a much higher
performance has been observed in comparison with Bicubic and other recent
methods too. Moreover, in downscaling cases with an odd scale factor, we give
an estimate of the mean squared error produced by our method and prove it is
theoretically null (hence PSNR equals to infinite and SSIM equals to one) in
absence of noise or initial artifacts on the input image.

###### keywords:

Image resizing, Image downscaling, Image upscaling , Lagrange interpolation ,
Chebyshev nodes

###### MSC:

94A08, 68U10, 41A05 , 62H35

††journal: ArXivmytitlenotemytitlenotefootnotetext: The research has been
accomplished within Research ITalian network on Approximation (RITA) and
Approximation Theory research group of Unione Matematica Italiana (TA–UMI). It
has been partially supported by GNCS–INdAM and University of Basilicata (local
funds).

##  1 Introduction

In this paper, we deal with the problem of image resizing. This has been
widely investigated over the past decades and is still an active research area
characterized by many applications in different domains, including image
transmission, satellite image analysis, gaming, remote sensing, etc. (see e.g.
[8, 22, 25, 28, 50, 51]) In literature, downscaling and upscaling are often
considered separate problems (see e.g. [7, 30, 50]), and most of the existing
methods are specialized in only one direction, sometimes for a limited range
of scaling factors (see e.g. [13, 46]). The method we are going to introduce
works in both down and up scaling directions and for ”large” scaling factors
as well. It falls into the class of the interpolation methods and it is based
on a non-standard modeling of the image resizing problem.

In order to introduce the adopted model we premise that, from the mathematical
viewpoint, a continuous image ${\mathcal{I}}$ is a function $f(x,y)$ of the
spatial coordinates which, without losing the generality, we can assume
belonging to the open square $A=]-1,1[^{2}$. Hence, its digital version $I$ of
$n\times m$ pixels is supposed to be made of the values that $f$ takes on a
discrete grid of nodes $X_{n\times m}\subset A$. Regarding such grid, it is
generally supposed that $X_{n\times m}=X_{n}\times X_{m}$ where we set

$X_{\mu}:=\\{x_{k}^{\mu}\ :\
k=1,\ldots,\mu\\}\subset]-1,1[,\qquad\forall\mu\in{\mathbb{N}}.$ (1)

A typical and natural choice of the (univariate) system of nodes $X_{\mu}$ is
to divide $[-1,1]$ in $(\mu+1)$ equal parts and to take the $\mu$ internal
equidistant nodes, i.e.

$X_{\mu}^{equ}=\left\\{-1+\frac{2k}{\mu+1}\ :\ k=1,\ldots,\mu\right\\}.$ (2)

However, it is well known that equally spaced nodes are not the best choice
for Lagrange interpolation since they lead to exponentially growing Lebesgue
constants, whereas optimal Lebesgue constants (growing at the minimal
projection rate) are provided by the Chebyshev nodes of the first kind

$x_{k}^{\mu}=\cos\left[\frac{(2k-1)\pi}{2\mu}\right],\qquad
k=1,\ldots,\mu,\qquad\forall\mu\in{\mathbb{N}},$ (3)

(see e.g. [27, 45] for a short excursion on the topic).

Discrete polynomial approximation based on Chebyshev zeros is a pillar in
approximation theory and practices. It has been widely studied and applied in
several fields and also recently, Chebyshev–like grids such as Xu points and
Padua points have been introduced to get optimal interpolation processes on
the square (see e.g. [48, 49, 9, 11, 33]). Nevertheless, to our knowledge, its
usage in image processing has been mainly limited to particular cases, such as
Magnetic Particle Imaging that is strictly related to Lissajous curves
generating the Padua points ( see e.g. [21, 18, 14]).

Figure 1: The first kind Chebyshev zeros of order $\mu=15$ (black stars) and
their arccosine (red stars)

Indeed, at a first look, Chebyshev zeros may seem unsuitable for sampling an
arbitrary image since they are not equidistant on the segment $]-1,1[$ of each
spatial coordinate, but they are arc sine distributed, becoming denser at the
endpoints $\pm 1$ (see Figure 1). Nevertheless, in our model, we transfer the
sampling question from the segment $]-1,1[$ to the semicircle of radius 1
centered at the axes origin. Thus, using the standard setting $t=\arccos x$,
instead of the usual equidistant nodes on the segment, our method takes the
following equidistant nodes on the semicircle

$t_{k}^{\mu}:=\frac{(2k-1)\pi}{2\mu},\qquad
k=1,\ldots,\mu,\qquad\mu\in{\mathbb{N}},$ (4)

These nodes define the Chebyshev zeros $x_{k}^{\mu}=\cos t_{k}^{\mu}$ in (3)
and the grid

$X_{n\times m}=\left\\{(x_{i}^{n},x_{j}^{m}):\quad i=1:n,\
j=1:m\right\\},\qquad n,m\in{\mathbb{N}}$ (5)

which we actually choose to sample arbitrary (continuous) images
${\mathcal{I}}$ and to obtain digital images of $n\times m$ pixels, for
arbitrary $n,m\in{\mathbb{N}}$.

According to this model, let $I=f(X_{n\times m})$ be the digital image of
$n\times m$ pixels, i.e.

$I_{i,j}:=f(x_{i}^{n},x_{j}^{m}),\qquad i=1:n,\ j=1:m,$ (6)

and let $I^{in}$ be its more or less corrupted version really available as
input data, say

$I^{in}_{i,j}:=\tilde{f}(x_{i}^{n},x_{j}^{m}),\qquad i=1:n,\ j=1:m,$ (7)

where $\tilde{f}$ denotes a corrupted version of the continuous image.
Starting from the data $I^{in}$, the resizing problem aims to find a good
approximation of the resampled image $I^{res}$ consisting of the samples of
$f$ at the different (more or less dense) grid $X_{N\times M}$, i.e.

$I^{res}_{k,h}:=f(x_{k}^{N},x_{h}^{M}),\qquad k=1:N,\ h=1:M.$ (8)

Hence, under our modeling, the resizing problem raises the following
approximation question: starting from the approximate values of $f$ at the
Chebyshev grid $X_{n\times m}$, how to approximate $f$ at another (coarser in
downscaling or finer in upscaling) Chebyshev grid $X_{N\times M}$? The
solution we propose relies on the global approximation of a function by its
bivariate Lagrange polynomial interpolating at the Chebyshev grid $X_{n\times
m}$.

This well–known interpolation polynomial can be easily deduced via tensor
product from the univariate case (see for instance [31, 32, 34, 35]) and fast
algorithms can be implemented for its computation ( see e.g. [33]). Moreover,
we recall that wavelets techniques can be applied to such kind of
interpolation [20, 15, 16].

Hence, in order to get the resampled image, instead of $f$, we propose to
sample a suitable Lagrange interpolation polynomial at the new grid. More
precisely, we use the pixels of the input image $I^{in}$ to build the
Lagrange-Chebyshev polynomial $L_{n,m}\tilde{f}$ interpolating $\tilde{f}$ at
$X_{n\times m}$ and we get the output image of the desired $N\times M$ size
(here denoted by $I^{out}$) by sampling $L_{n,m}\tilde{f}$ at the mesh
$X_{N\times M}$.

To test the performance of such Lagrange–Chebyshev Interpolation (LCI) method
we have considered the usual SSIM (Structured Similarity Index Measurement)
and PSNR (Peak Signal to Noise Ratio) metrics as computed by MatLab, and we
have carried out the experimentation on several kinds of images collected in 5
datasets, having different characteristics on the size and the quality of the
images, the variety of subjects, etc.

For an initial comparison, we focused on interpolation methods working as well
in both down and up scaling directions. In general, these methods find a
resized image by local or global interpolation of the continuous image $f$,
using in some way the initial $n\times m$ pixels, namely the more or less
corrupted values of $f$ at the grid of nodes $X_{n\times m}$. Typically, for
such grid it is assumed the equidistant nodes model (2) and most of the
interpolation methods act locally on each pixel of the input image. This is
the case of the classical bicubic interpolation method BIC [19] implemented by
the Matlab built–in function imresize which we have used either in upscaling
and downscaling to get a first comparison.

By the extensive experimentation we carried out, we assess that the
performance of the LCI method is superior in downscaling (d-LCI method), while
in upscaling (u-LCI method) there are only slight differences and the two
methods seem almost comparable.

Thus it seemed relevant to us to further investigate only the downscaling case
by testing d-LCI method in comparison with two other recent downscaling
methods, here denoted by DPID [47] and L0 [23], both of them not belonging to
the family of interpolation methods for downscaling. The effective better
performance in downscaling has been confirmed also in these cases, with an
increasing gap as the downscaling factor increases. Moreover, by running the
public available Matlab code, we have observed that DPID and L0 methods work
only for integer scale factors $n/N$ and $m/M$ and they are not always
implementable for images with large sizes due to too long running times or too
much memory. On the contrary, LCI method keeps the runtimes contained and
offers high flexibility since it works in both downscaling and upscaling with
non integer scale factors too, being possible to specify the desired final
size or the scale factor as input parameters. We remark that similar nice
features are shared by the Matlab function imresize, but with lower
performances in downscaling. Moreover, in downscaling, for all odd scale
factors, we state a theoretical estimate of the Mean Squared Error (MSE)
obtained by our d-LCI method arriving to prove that it is null if $I^{in}=I$.
Such a theoretical result has been confirmed by our experimentation that also
shows it does not hold for imresize applied to the same input images.

The paper is organized as follows. For a simpler exposition of the idea behind
the LCI method, we first introduce the method for a monochrome image in
Section 2 and then, in Section 3, we expose it for RGB color images. Both
monochrome and color cases are considered in Section 4 where some advantages
of our approach are discussed. Finally, Section 5 consists of all the
numerical experiments.

##  2 Re-sampling gray-level images

In the case of a monochrome image, $\mathcal{I}$ is represented at a
continuous scale by a scalar function $f(x,y)$ whose values are the grey
levels at the spatial coordinates $(x,y)\in A$. Moreover, at a discrete scale,
digital images of $n\times m$ pixels are matrices $I\in{\mathbb{R}}^{n\times
m}$ whose elements are the samples of $f$ at the discrete nodes set
$X_{n\times m}$ defined in our model by (5) and (3). Then, according to the
notation introduced in the previous section, the resizing problem of
$\mathcal{I}$ (in upscaling or in downscaling, respectively) consists in
finding a good approximation of the matrix $I^{res}\in{\mathbb{R}}^{N\times
M}$ given by (8) starting from the matrix (having a reduced or enlarged size,
respectively) $I^{in}\in{\mathbb{R}}^{n\times m}$ defined as in (7) being
$\tilde{f}$ a more or less corrupted version of $f$.

We approach such an approximation problem by using the following bivariate
(tensor product) Lagrange–Chebyshev polynomial based on the available data
$I^{in}$

$L_{n,m}\tilde{f}(x,y):=\sum_{i=1}^{n}\sum_{j=1}^{m}\tilde{f}(x_{i}^{n},x_{j}^{m})\ell_{i}^{n}(x)\ell_{j}^{m}(y),\qquad(x,y)\in
A,$ (9)

where, for all $\mu\in{\mathbb{N}}$, $\ell_{k}^{\mu}$ denotes the $k-$th
fundamental Lagrange polynomial related to the nodes system $X_{\mu}$, namely

$\ell_{k}^{\mu}(\xi):=\prod_{\scriptsize\begin{array}[]{c}s=1\\\\[-2.84544pt]
s\neq k\end{array}}^{\mu}\frac{\xi-
x_{s}^{\mu}}{x_{k}^{\mu}-x_{s}^{\mu}},\qquad\xi\in[-1,1],\qquad,\qquad
k=1:\mu.$ (10)

It is well–known that the Lagrange–Chebyshev polynomial in (9) interpolates
$\tilde{f}$ at the grid $X_{n\times m}$, i.e.

$L_{n,m}\tilde{f}(x_{i}^{n},x_{j}^{m})=\tilde{f}(x_{i}^{n},x_{j}^{m})=I^{in}_{i,j},\qquad
i=1:n,\ j=1:m.$ (11)

The samples of such polynomial at the re-scaled grid $X_{N\times M}$
constitute the approximate resized image provided by the LCI method. Hence,
the LCI method computes the matrix $I^{out}\in{\mathbb{R}}^{N\times M}$ whose
entries are the following

$I^{out}_{k,h}:=L_{n,m}\tilde{f}(x_{k}^{N},x_{h}^{M}),\qquad k=1:N,\quad
h=1:M.$ (12)

Introducing the Vandermonde-like matrices

$V_{1}:=\left[\ell_{i}^{n}(x_{k}^{N})\right]_{i,k}\in{\mathbb{R}}^{n\times
N},\qquad
V_{2}:=\left[\ell_{j}^{m}(x_{h}^{M})\right]_{j,h}\in{\mathbb{R}}^{m\times M},$
(13)

by (9) and (12), the output matrix $I^{out}$ can be computed from the input
matrix $I^{in}$ according to the following matrices identity

$I^{out}=V_{1}^{T}I^{in}V_{2},$ (14)

Note that in case we have to resize a lot of images for the same fixed sizes,
formula (14) also allows working in parallel with pre-computed matrices
$V_{i}$ .

Moreover, it is well known that the fast computation of the matrices $V_{i}$
can be achieved by using, instead of (10), the following more convenient form
of the fundamental Lagrange polynomial

$\ell_{k}^{\mu}(\cos
t)=\frac{2}{\mu}\sum_{r=0}^{\mu-1}{}^{\prime}\cos\left[\frac{(2k-1)r\pi}{2\mu}\right]\cos\left[rt\right],\qquad
k=1:\mu,\quad t\in[0,\pi],$ (15)

where, as usual, the prime on the summation symbol means that the first
addendum is halved. Hence, by (15) we get that the matrices $V_{1}$ and
$V_{2}$ can be computed by using fast cosine transform algorithms (see e.g.
[37]) being

$\displaystyle(V_{1})_{i,k}=\ell_{i}^{n}(x_{k}^{N})$ $\displaystyle=$
$\displaystyle\frac{2}{n}\sum_{r=0}^{n-1}{}^{\prime}\cos\left[\frac{(2i-1)r\pi}{2n}\right]\cos\left[\frac{(2k-1)r\pi}{2N}\right],$
$\displaystyle\qquad i=1:n,\ k=1:N,$
$\displaystyle(V_{2})_{j,h}=\ell_{j}^{m}(x_{h}^{M})$ $\displaystyle=$
$\displaystyle\frac{2}{m}\sum_{s=0}^{m-1}{}^{\prime}\cos\left[\frac{(2j-1)s\pi}{2m}\right]\cos\left[\frac{(2h-1)s\pi}{2M}\right],$
$\displaystyle\qquad j=1:m,\ h=1:M.$

##  3 Re-sampling RGB color images

Now we are going to introduce the method for a general color image
${\mathcal{I}}$. In this case, behind ${\mathcal{I}}$ there is a vector
function $\mathbf{f}:A\rightarrow{\mathbb{R}}^{3}$ whose components
$\mathbf{f}=(f_{R},f_{G},f_{B})$ represent $\mathcal{I}$ in the RGB color
space. According to our model, the input image $\mathbf{I}^{in}$ and the
target resized image $\mathbf{I}^{res}$ are represented in the RGB space as
the follows

$\begin{array}[]{ll}\mathbf{I}^{in}\in{\mathbb{R}}^{n\times m\times
3},&\mathbf{I}^{in}\equiv(I^{in}_{R},I^{in}_{G},I^{in}_{B})\\\
\mathbf{I}^{res}\in{\mathbb{R}}^{N\times M\times
3},&\mathbf{I}^{res}\equiv(I^{res}_{R},I^{res}_{G},I^{res}_{B}),\end{array}$

where we assume that

$\begin{array}[]{lll}\left(I^{in}_{R}\right)_{i,j}={\tilde{f}_{R}}(x_{i}^{n},x_{j}^{m}),&&\left(I^{res}_{R}\right)_{k,h}={f_{R}}(x_{k}^{N},x_{h}^{M}),\\\
\left(I^{in}_{G}\right)_{i,j}={\tilde{f}_{G}}(x_{i}^{n},x_{j}^{m}),&&\left(I^{res}_{G}\right)_{k,h}={f_{G}}(x_{k}^{N},x_{h}^{M}),\\\
\left(I^{in}_{B}\right)_{i,j}={\tilde{f}_{B}}(x_{i}^{n},x_{j}^{m}),&&\left(I^{in}_{B}\right)_{k,h}={f_{B}}(x_{k}^{N},x_{h}^{M}).\end{array}$
(18)

holds for all $[i,j]=[1:n,1:m]$ and $[k,h]=[1:N,1:M]$.

Thus, similarly to the monochrome case and with obvious meaning of the
notation, we start from

$\mathbf{I}^{in}_{i,j}=\tilde{\mathbf{f}}(x_{i}^{n},x_{j}^{m})\qquad
i=1:n,\quad j=1:m,$ (19)

and approximate $\mathbf{I}^{res}$ by

$\mathbf{I}^{out}_{k,h}=L_{n,m}\tilde{\mathbf{f}}(x_{k}^{N},x_{h}^{M})\qquad
k=1:N,\quad\ h=1:M$ (20)

i.e., for $k=1:N$ and $h=1:M$, we define

$\mathbf{I}^{out}_{k,h}=\left(L_{n,m}{\tilde{f}_{R}}(x_{k}^{N},x_{h}^{M}),\
L_{n,m}{\tilde{f}_{G}}(x_{k}^{N},x_{h}^{M}),\
L_{n,m}{\tilde{f}_{B}}(x_{k}^{N},x_{h}^{M})\right).$ (21)

Hence, by applying the same argument of the monochrome case, we get that the
RGB components of the output image
$\mathbf{I}^{out}\equiv(I^{out}_{R},I^{out}_{G},I^{out}_{B})$ are given by

$\mathbf{I}^{out}=(V_{1}^{T}I^{in}_{R}V_{2},\ V_{1}^{T}I^{in}_{G}V_{2},\
V_{1}^{T}I^{in}_{B}V_{2}),$

with $V_{1},V_{2}$ defined in (2)–(2).

##  4 Model analysis

For a general analysis of the error we get in approximating the target resized
image $I^{res}=[{\mathbf{f}}(x_{i}^{n},x_{j}^{m})]_{i,j}$ with the output
image $I^{out}=[\L_{n,m}\tilde{\mathbf{f}}(x_{h}^{N},x_{k}^{M})]_{h,k}$
produced by LCI method, we refer the reader to the wide existing literature on
the Lagrange interpolation error estimates (see e.g. [27] and the references
therein). In this paper, we are interested to measure the performance of LCI
method through some standard metrics usually used in Image Processing (see
e.g. [38] and the references therein). In particular, we focus on the Mean
Squared Error (MSE) defined as follows

$\rm{MSE}(I^{res},I^{out})=\begin{cases}\displaystyle\frac{1}{MN}\|I^{res}-I^{out}\|_{F}^{2},\hskip
99.58464pt\textrm{ gray-level images},\\\ \\\
\frac{1}{3NM}\left\\{\|I_{R}^{res}-I_{R}^{out}\|_{F}^{2}+\|I_{G}^{res}-I_{G}^{out}\|_{F}^{2}+\|I^{res}_{B}-I_{B}^{out}\|_{F}^{2}\right\\},\\\
\hskip 199.16928pt\textrm{ RGB color images},\end{cases}$

being $\|\cdot\|_{F}$ the Frobenius norm

$\|A\|_{F}:=\left(\sum_{k=1}^{N}\sum_{h=1}^{M}A_{k,h}^{2}\right)^{\frac{1}{2}},\qquad
A=(A_{k,h})\in{\mathbb{R}}^{N\times M}.$

Moreover, we consider the Peak Signal to Noise Ratio (PSNR) defined by the
previous MSE as follows

$\rm{PSNR}(I^{res},I^{out})=20\displaystyle\log_{10}\left(\frac{\max_{f}}{\sqrt{\rm{MSE}(I^{res},I^{out})}}\right),$
(22)

with $\max_{f}=255$, since we use the representation by 8 bits per sample.

Finally, in our experiments we evaluate the Structured Similarity Index
Measurement (SSIM). For grey-level images it is defined as

$\rm{SSIM}(I^{res},I^{out})=\displaystyle\frac{\left[2\mu(I^{res})\mu(I^{out})+c_{1}\right]\left[2\rm{cov}(I^{res},\
I^{out})+c_{2}\right]}{\left[\mu^{2}(I^{res})+\mu^{2}(I^{out})+c_{1}\right]\left[\sigma^{2}(I^{res})+\sigma^{2}(I^{out})+c_{2}\right]},$
(23)

where $\mu(A),\sigma(A)$ and $\rm{cov}(A,B)$ indicate the average, variance
and covariance, respectively, of the matrices $A,B$, and $c_{1},c_{2}$ are
constants usually fixed as $c_{1}=(0.01\times L),c_{2}=(0.03\times L)$ with
the dynamic range of the pixel values $L=255$ for 8-bit images.

In the case of RGB color images, making the conversion to the color space
YCbCr, the $\rm{SSIM}$ is computed by the above definition applied to the
intensity Y (luma) channel.

Based on the previous metrics, the results of wide experimentation will be
reported in the next section. In the sequel, we are going to underly two
particular features of LCI method.

A first aspect is related to the model choice which, for odd downscaling
factors $s:=n/N=m/M>1$, leads to state the following

######  Proposition 4.1.

Under the previously introduced notation, if there exists
$\ell\in{\mathbb{N}}$ s.t. $n=(2\ell-1)N$ and $m=(2\ell-1)M$ hold, then we
have

$\rm{MSE}(I^{res},I^{out})\leq s^{2}MSE(I,I^{in}),\qquad s:=(2\ell-1).$ (24)

Moreover, if $I=I^{in}$ then we get the best results

$\rm{MSE}(I^{res},I^{out})=0,\qquad\mbox{and}\qquad\rm{SSIM}(I^{res},I^{out})=1.$
(25)

Proof.  Let’s give the proof in the case of a gray-level image since for RGB
color images the statement will easily follow by considering the single RGB
components.

As a keynote of the proof, we observe that for any $\ell\in{\mathbb{N}}$, we
have

$n=(2\ell-1)N\qquad\Longrightarrow\qquad X_{N}\subset X_{n}.$

More precisely, setting $s:=(2\ell-1)$, it is easy to check that

$n=sN\qquad\Longrightarrow\qquad x_{k}^{N}=x_{i}^{n}\quad\mbox{with}\quad
i=\frac{s(2k-1)+1}{2},\quad k=1:N,$ (26)

where we remark that for all $k=1:N$, the numerator $s(2k-1)+1$ is certainly
even since $s$ is odd.

Consequently, we deduce (24) from (26) as follows

$\displaystyle\rm{MSE}(I^{res},I^{out})=\frac{1}{NM}\sum_{k=1}^{N}\sum_{h=1}^{M}\left[f\left(x_{k}^{N},x_{h}^{M}\right)-L_{n,m}\tilde{f}\left(x_{k}^{N},x_{h}^{M}\right)\right]^{2}$
$\displaystyle=$
$\displaystyle\frac{1}{NM}\sum_{k=1}^{N}\sum_{h=1}^{M}\left[f\\!\left(x_{\frac{s(2k-1)+1}{2}}^{n},x_{\frac{s(2h-1)+1}{2}}^{m}\right)\\!-\\!\tilde{f}\\!\left(x_{\frac{s(2k-1)+1}{2}}^{n},x_{\frac{s(2h-1)+1}{2}}^{m}\right)\right]^{2}$
$\displaystyle\leq$
$\displaystyle\frac{1}{NM}\sum_{i=1}^{n}\sum_{j=1}^{m}\left[f\left(x_{i}^{n},x_{j}^{m}\right)-\tilde{f}\left(x_{i}^{n},x_{j}^{m}\right)\right]^{2}=s^{2}MSE(I,I^{in}).$

Moreover, in the case that $I=I^{in}$, we have that
$L_{n,m}f=L_{n,m}\tilde{f}$ and, due to the nesting property $X_{N\times
M}\subset X_{n\times m}$, by (11), for any $h=1:N$ and $k=1:M$, we get

$I^{out}_{h,k}=L_{n,m}\tilde{f}(x_{h}^{N},x_{k}^{M})=L_{n,m}f(x_{h}^{N},x_{k}^{M})=f(x_{h}^{N},x_{k}^{M})=I^{res}_{h,k},$
(27)

that implies the first equation in (25). Finally, (27) also yields

$2\rm{cov}(I^{res},\ I^{out})=\sigma(I^{res})^{2}+\sigma(I^{out})^{2},$
$2\mu(I^{res})\mu(I^{out})=\mu(I^{res})^{2}+\mu(I^{out})^{2}$

which conclude the proof of (25). $\diamondsuit$

We remark that the previous proposition is a direct consequence of our choice
of the sampling model based on Chebyshev zeros (3), instead of the usual
equally spaced nodes (2).

Such choice (and more generally the choice of ”good” interpolation knots)
allows to use global interpolation processes that are not possible in the case
of equally spaced nodes. In fact, as second aspect, we aim to highlight the
relevance of the choice of good interpolation knots by showing the effects we
have in the performance of global Lagrange interpolation methods associated
with the univariate nodes set $X_{\mu}^{equ}$ in (2). The disastrous effects
of the exponential growth of Lebesgue constants are visible in the next
experiment concerning an image zooming $\times 2$ from the size $n\times n$
with $n=256$ to the size $N\times N$ with $N=512$. Figure 2 displays the
images obtained starting from the same input image (Fig.2, left) and using the
Lagrange interpolation polynomial based on equally spaced nodes
$X_{n}^{equ}\times X_{n}^{equ}$ (Fig.2, right) and on the Chebyshev grid
$X_{n\times n}$ (Fig.2, middle). The images we see are the samples of the
previous Lagrange polynomials at the respective $(N\times N)$–grid and we can
check how the equally–spaced sampling model (2) yields wrong results outside
of a small region of the output image.

Figure 2: The well–known image Flowers taken from USC–SIPI [1] (left)
upsampled at the double scale by Lagrange interpolation at Chebyshev nodes
(middle) and at equally spaced nodes (right)

##  5 Experimental results

In this section, we propose a selection of the extensive experimentation we
carried on to test the benefits of our procedure in downscaling (d-LCI method)
and in upscaling (u-LCI method), comparing our results with other resizing
techniques.

###  5.1 Comparison methods

As already mentioned, we compare LCI with BIC provided by MatLab imresize with
’bicubic’ option where we recall that a new pixel is determined from a
weighted average of the 16 closest pixels using an interpolating cubic
convolution function satisfying prescribed assumptions of smoothness on $f$ to
gain at least a cubic order of convergence [19]. The comparison BIC–LCI
regards either down and up resizing cases, giving as input either the scale
factor or the final size of the desired image. Note that we have also tested
the other different options of imresize (namely ’linear’ and ’nearest
neighborhood’) but for brevity, we do not report the results since they give
no new insight.

Moreover, to further investigate the performance in downscaling case, we
compare d-LCI with other two recent downscaling methods not belonging to the
interpolation method class. These methods are briefly indicated as DPID method
(described in [47] with the code available at [2] ) and L0 method (described
in [23], with the code available at [3]). To be fair, we remark that we forced
DPID and L0 also in upscaling direction for several scaling factors (the
choice of the desired size is not allowed by DPID and L0) but the results were
very poor w.r.t. BIC and LCI and they have been not reported here.

Finally, we point out that all the previous methods have run on the same PC
with Intel Core i7 3770K CPU @350GHz configuration.

###  5.2 Datasets

Since the results of any method depend on the type of input image, we have
conducted the experimentation on different kinds of 8–bits images collected in
five publicly available datasets, whose characteristics (acronym, number of
images, range of their sizes) are synthesized in Table 1.

Table 1: Datasets list Dataset | n. of Images | Sizes
---|---|---
BSDS500 | 500 |  $481\times 321$ or $321\times 481$
NASA | 17 | from 500$\times$334 to 6394$\times$3456
YAHOO | 96 | from 500$\times$334 to 6394$\times$3456
13US | 13 | from 241$\times$400 to 400$\times$310
URBAN100 | 100 | from 1024$\times$564 to 1024$\times$1024

We point out that we have chosen the datasets in Table 1 since they are the
same considered in [47, 23] for DPID and L0 methods. In particular:

  * 1.

BSDS500 dataset [26], available at [4], has been used for testing L0 in [23].
It is sufficiently general and provides a large variety of images often
employed for testing many other methods with different image analysis tasks
such as image segmentation [29, 39, 40, 42]), color quantization [12, 41, 10,
43], etc.

  * 2.

NASA Image Gallery [6] and YFCC100M (Yahoo Flickr Creative Commons 100
Million) [44] datasets (here briefly denoted by NASA and YAHOO, respectively)
are used by DPID in [47].

  * 3.

13US [36] contains natural images, available at [5] and originally taken from
the MSRA Salient Object Database [24]. It is another dataset used in [47].

  * 4.

URBAN100 dataset [17], concerning urban scenes with images having one
dimension equal to 1024, is commonly used to evaluate the performance of
super-resolution models. It is used as test images for L0 in [23].

All the previous datasets have been considered for either up and down scaling.

We remark that in all the performed experiments the input images are not
available from the datasets. Hence, following a typical approach in image
quantitative evaluation, we fix all the images from the datasets as target
images (i.e. $I^{res}$ in our notation) and we apply BIC to generate the input
image (namely $I^{in}$) common to all methods. More precisely, to run [$\times
s$] upscaling methods ([$:s$] downscaling methods, resp.) we generate the
input image $I^{in}$ by applying BIC to $I^{res}$ in the opposite [$:s$]
([$\times s$], resp.) scaling direction.
Note that, according to such procedure, in testing [$\times s$] upscaling
methods we may find that $I^{res}$ does not have both the dimensions $N$ and
$M$ that are divisible for $s$. In this case, in order to generate $I^{in}$ we
use imresize with the size option, requiring $n=\lfloor\frac{N}{s}\rfloor$ and
$m=\lfloor\frac{M}{s}\rfloor$. Moreover, once obtained $I^{in}$, both BIC and
u-LCI run by specifying again the size $N\times M$ instead of the scaling
factor.

###  5.3 Visual comparison

As first experiment, we focus on the visual comparison of some original images
chosen in the previous datasets with the output images produced by the methods
in Subsection 5.1. Such images are given in Figure 3 for upscaling and in
Figure 4 for downscaling.

Figure 3: Examples of upscaling performance results at the scale factor 2
(top), at the scale factor 3 (middle), at the scale factor 4 (bottom).

In both the figures at the first column we show the target images and at the
successive columns we display the resized images obtained as the output of the
considered methods. Moreover, the first, second and third row of both the
figures correspond to the scaling factors $2,3$ and $4$, respectively. In all
the cases the images have been reported with evidence of some magnified
regions of interest (ROI).

Figure 4: Examples of downscaling performance results at the scale factor 2
(top), at the scale factor 3 (middle), at the scale factor 4 (bottom).

From the qualitative point of view, by inspecting Figures 3 and 4, we can
deduce that in terms of visual quality LCI has a good performance, also with
respect to the chosen comparing methods. It preserves the visual structure of
the object without losing image details, the local contrast, and the luminance
of the input image, by generating resized images close to the target ones.
Ringing and over smoothing artifacts are limited and the resized images seem
fair, sufficiently not blurred, and with well-balanced colors.

###  5.4 Quantitative comparison: average results

For each dataset selected in Subsection 5.2, we computed the averages of the
PSNR and SSIM values achieved by LCI and by the comparison methods described
in Subsection 5.1, for the scaling factors $2,3,4$. The results are displayed
in Table 2 for upscaling and in Table 3 for downscaling, reporting in both the
tables also the averages of the required CPU times.

Table 2: Average results in upscaling. | x2 | x3 | x4
---|---|---|---
| PSNR | SSIM | T | PSNR | SSIM | T | PSNR | SSIM | T
BSDS500 |  |  |  |  |  |  |  |  |
BIC | 26,341 | 0,886 | 0,003 | 24,821 | 0,837 | 0,002 | 22,251 | 0,701 | 0,003
u-LCI | 26,381 | 0,888 | 0,014 | 24,868 | 0,839 | 0,010 | 22,446 | 0,769 | 0,008
NASA |  |  |  |  |  |  |  |  |
BIC | 34,806 | 0,958 | 0,091 | 31,205 | 0,924 | 0,074 | 29,808 | 0,907 | 0,071
u-LCI | 35,412 | 0,960 | 1,357 | 31,542 | 0,924 | 0,839 | 30,183 | 0,908 | 0,634
YAHOO |  |  |  |  |  |  |  |  |
BIC | 33,391 | 0,953 | 0,056 | 29,754 | 0,913 | 0,055 | 28,741 | 0,891 | 0,051
u-LCI | 33,900 | 0,955 | 0,812 | 29,956 | 0,914 | 0,567 | 28,986 | 0,891 | 0,459
URBAN100 |  |  |  |  |  |  |  |  |
BIC | 25,433 | 0,882 | 0,009 | 21,306 | 0,755 | 0,008 | 21,703 | 0,741 | 0,008
u-LCI | 25,896 | 0,886 | 0,064 | 21,325 | 0,754 | 0,051 | 21,919 | 0,793 | 0,059
13US |  |  |  |  |  |  |  |  |
BIC | 24,116 | 0,861 | 0,002 | 20,754 | 0,734 | 0,002 | 20,545 | 0,710 | 0,002
u-LCI | 24,486 | 0,868 | 0,012 | 20,780 | 0,738 | 0,010 | 20,656 | 0,713 | 0,009
Table 3: Average results in downscaling  | :2 | :3 | :4
---|---|---|---
| PSNR | SSIM | T | PSNR | SSIM | T | PSNR | SSIM | T
BSDS500 |  |  |  |  |  |  |  |  |
BIC | 38,872 | 0,993 | 0,006 | 39,253 | 0,993 | 0,009 | 39,183 | 0,993 | 0,017
DPID | 41,745 | 0,996 | 7,696 | 41,827 | 0,996 | 12,246 | 41,206 | 0,996 | 18,615
L0 | 29,317 | 0,961 | 3,647 | 32,873 | 0,971 | 8,020 | 34,174 | 0,971 | 14,295
d-LCI | 53,732 | 1,000 | 0,057 | Inf | 1,000 | 0,091 | 55,889 | 1,000 | 0,137
NASA |  |  |  |  |  |  |  |  |
BIC | 45,969 | 0,995 | 0,229 | 47,114 | 0,996 | 0,325 | 46,973 | 0,996 | 0,639
DPID | 47,675 | 0,998 | 448,023 | 47,657 | 0,998 | 731,498 | 47,112 | 0,997 | 1.098,113
L0 | 34,754 | 0,972 | 208,574 | 37,285 | 0,979 | 617,386 | oom | oom | oom
d-LCI | 54,265 | 0,999 | 6,614 | Inf | 1,000 | 13,317 | 57,491 | 1,000 | 22,859
YAHOO |  |  |  |  |  |  |  |  |
BIC | 44,757 | 0,996 | 0,155 | 45,858 | 0,997 | 0,219 | 45,682 | 0,996 | 0,422
DPID | 46,743 | 0,998 | 291,685 | 46,948 | 0,998 | 479,638 | 46,421 | 0,998 | 714,908
L0 | 33,913 | 0,974 | 133,190 | oom | oom | oom | oom | oom | oom
d-LCI | 54,067 | 0,999 | 4,698 | Inf | 1,000 | 8,279 | 57,223 | 1,000 | 13,898
URBAN00 |  |  |  |  |  |  |  |  |
BIC | 35,661 | 0,989 | 0,027 | 36,010 | 0,990 | 0,041 | 35,940 | 0,989 | 0,068
DPID | 39,178 | 0,996 | 37,592 | 39,178 | 0,996 | 60,834 | 38,702 | 0,995 | 93,996
L0 | 26,718 | 0,951 | 13,267 | 31,061 | 0,969 | 28,845 | 33,571 | 0,973 | 50,312
d-LCI | 52,613 | 0,999 | 0,234 | Inf | 1,000 | 0,416 | 55,452 | 1,000 | 0,618
13US |  |  |  |  |  |  |  |  |
BIC | 35,129 | 0,990 | 0,005 | 35,469 | 0,991 | 0,009 | 35,397 | 0,991 | 0,013
DPID | 38,061 | 0,996 | 5,593 | 38,326 | 0,996 | 8,905 | 37,707 | 0,995 | 13,819
L0 | 25,521 | 0,949 | 2,572 | 30,231 | 0,972 | 5,706 | 32,929 | 0,979 | 9,939
d-LCI | 52,813 | 1,000 | 0,042 | Inf | 1,000 | 0,073 | 55,374 | 1,000 | 0,108

As announced in the introduction, from Table 2 we see that the values attained
by u-LCI are only slightly better than those achieved by BIC, and the
execution time of u-LCI is a little bit greater than BIC. Of course, as the
scale factor increases, both PSNR and SSIM decrease in all the upscaling
methods. Moreover, the gaps between BIC and u-LCI in both the metrics are
maintained also for large size images.

On the contrary, in downscaling, by inspecting Table 3 we observe that the
performance provided by d-LCI is much better than those achieved by the other
methods for both the metrics. In particular, when the downscaling factor is
$3$ the theoretical results claimed in Proposition 4.1 are confirmed.
Moreover, looking at the even downscaling factors $2,4$ the improvement by
d-LCI seems to increase as the scale factor increases while BIC seems to be
affect by saturation, with an increasing gap between d-LCI and the remaining
methods.

About the CPU time, except for BIC, d-LCI is faster than the other methods
that, in case of very large images, need very long computational time (see the
results for NASA and YAHOO datasets that include target images of size
$6394\times 3456$, i.e. input images with size $25576\times 13824$ in case of
downscaling factor equals to 4) and in some cases, wherever we see oom (which
means out of memory), the available code of L0 does not arrive to produce any
result.

###  5.5 Quantitative comparison: some pointwise results

Here we focus on each single image from the smaller datasets of Table 1 ,
namely we consider 13US dataset in downscaling, with target images of small
size displayed in Figure 5, and NASA dataset in upscaling, with target images
of large size displayed in Figure 6.

Figure 5: 13US dataset image: [U1 -U13] from left to right and from top to
bottom . Figure 6: NASA dataset images: [N1, N2, N3, N5, N6, N7, N8, N9, N10,
N11, N12, N13, N14, N15, N16, N17, N19] from left to right and from top to
bottom .

The PSNR values achieved for every single image have been plotted in Figure 7
for NASA dataset with upscaling factor $2,4,8,16$, and in Figure 8 for 13US
dataset with downscaling factor $6,18,30$.

Figure 7: Pointwise PSNR for the NASA set for scale factors s=2,4,8,16

Figure 8: Pointwise PSNR values for the set 13US, with scale factors s=6
(top), s=18 (middle), s=30 (down)

On the same datasets, more detailed results are given in Tables 4–6 reporting
in the first columns the name and the size of all the images from 13US and
NASA.

Table 4: Pointwise results on NASA dataset in upscaling | x2 | x4 | x8 | x16
---|---|---|---|---
| PSNR | SSIM | PSNR | SSIM | PSNR | SSIM | PSNR | SSIM
N1 4928$\times$3280 |  |  |  |  |  |  |  |
BIC | 37,472 | 0,945 | 33,453 | 0,892 | 29,828 | 0,858 | 26,795 | 0,833
u-LCI | 37,903 | 0,950 | 33,856 | 0,892 | 30,133 | 0,855 | 26,999 | 0,831
N2 3072$\times$2304 |  |  |  |  |  |  |  |
BIC | 44,235 | 0,986 | 37,820 | 0,952 | 32,417 | 0,888 | 29,065 | 0,840
u-LCI | 45,203 | 0,988 | 38,621 | 0,957 | 32,943 | 0,894 | 29,183 | 0,840
N3 2048$\times$1536 |  |  |  |  |  |  |  |
BIC | 35,708 | 0,963 | 30,845 | 0,919 | 27,474 | 0,874 | 24,593 | 0,830
u-LCI | 36,552 | 0,965 | 31,136 | 0,916 | 27,731 | 0,865 | 24,850 | 0,808
N5 2701$\times$3665 |  |  |  |  |  |  |  |
BIC | 31,743 | 0,987 | 25,444 | 0,960 | 22,262 | 0,932 | 20,388 | 0,913
u-LCI | 33,077 | 0,989 | 25,792 | 0,962 | 22,426 | 0,933 | 20,487 | 0,913
N6 1836$\times$3264 |  |  |  |  |  |  |  |
BIC | 33,011 | 0,986 | 28,291 | 0,965 | 24,250 | 0,927 | 21,677 | 0,886
u-LCI | 33,443 | 0,986 | 28,779 | 0,961 | 24,614 | 0,908 | 21,784 | 0,850
N7 1600$\times$1200 |  |  |  |  |  |  |  |
BIC | 29,939 | 0,886 | 25,087 | 0,670 | 23,364 | 0,567 | 22,627 | 0,545
u-LCI | 30,981 | 0,905 | 25,409 | 0,684 | 23,542 | 0,569 | 22,683 | 0,546
N8 3008$\times$2000 |  |  |  |  |  |  |  |
BIC | 32,298 | 0,987 | 29,295 | 0,974 | 26,346 | 0,956 | 23,345 | 0,932
u-LCI | 32,577 | 0,988 | 29,594 | 0,976 | 26,619 | 0,958 | 23,508 | 0,933
N9 5430$\times$3520 |  |  |  |  |  |  |  |
BIC | 35,440 | 0,985 | 28,707 | 0.945 | 24,336 | 0,883 | 22,494 | 0,846
u-LCI | 36,116 | 0,986 | 30,735 | 0,955 | 26,121 | 0,899 | 22,712 | 0,845
N10 3264$\times$2448 |  |  |  |  |  |  |  |
BIC | 35,758 | 0,973 | 28,875 | 0,916 | 24,213 | 0,837 | 20,914 | 0,767
u-LCI | 37,047 | 0,977 | 29,531 | 0,920 | 24,529 | 0,838 | 21,187 | 0,768
N11 4368$\times$2326 |  |  |  |  |  |  |  |
BIC | 36,525 | 0,982 | 30,149 | 0,955 | 24,649 | 0,906 | 20,966 | 0,849
u-LCI | 37,461 | 0,981 | 30,986 | 0,951 | 24,918 | 0,897 | 21,193 | 0,846
N12 3504$\times$2336 |  |  |  |  |  |  |  |
BIC | 40,471 | 0,989 | 32,786 | 0,951 | 27,777 | 0,888 | 25,116 | 0,850
u-LCI | 41,889 | 0,991 | 35,581 | 0,956 | 27,982 | 0,890 | 25,191 | 0,851
N13 1200$\times$1600 |  |  |  |  |  |  |  |
BIC | 24,145 | 0,746 | 22,370 | 0,603 | 21,838 | 0,530 | 21,555 | 0,483
u-LCI | 24,316 | 0,730 | 22,390 | 0,592 | 21,852 | 0,529 | 21,579 | 0,486
N14 3456$\times$5184 |  |  |  |  |  |  |  |
BIC | 45,041 | 0,992 | 40,137 | 0,986 | 35,808 | 0,978 | 31,982 | 0,967
u-LCI | 45,539 | 0,992 | 40,907 | 0,986 | 36,531 | 0,978 | 32,407 | 0,966
N15 3072$\times$2048 |  |  |  |  |  |  |  |
BIC | 34,457 | 0,990 | 29,574 | 0,967 | 26,731 | 0,941 | 24,886 | 0,920
u-LCI | 35,427 | 0,991 | 29,872 | 0,968 | 26,883 | 0,940 | 24,976 | 0,919
N16 5501$\times$3095 |  |  |  |  |  |  |  |
BIC | 47,929 | 0,993 | 40,846 | 0,982 | 34,304 | 0,961 | 29,760 | 0,940
u-LCI | 48,691 | 0,994 | 42,202 | 0,983 | 35,029 | 0,961 | 30,027 | 0,939
N17 2048$\times$1363 |  |  |  |  |  |  |  |
BIC | 30,892 | 0,980 | 27,028 | 0,954 | 24,688 | 0,927 | 22,895 | 0,899
u-LCI | 31,380 | 0,981 | 27,219 | 0,954 | 24,834 | 0,927 | 21,044 | 0,877
N19 3039$\times$4559 |  |  |  |  |  |  |  |
BIC | 34,581 | 0,977 | 27,048 | 0,884 | 21,873 | 0,718 | 18,910 | 0,615
u-LCI | 35,844 | 0,981 | 27,881 | 0,896 | 22,178 | 0,723 | 19,077 | 0,617
Table 5: Pointwise results on 13US dataset in downscaling | :3 | :6 | :9 | :18
---|---|---|---|---
| PSNR | SSIM | PSNR | SSIM | PSNR | SSIM | PSNR | SSIM
US1 300$\times$400 |  |  |  |  |  |  |  |
BIC | 38,625 | 0,996 | 38,573 | 0,996 | 38,583 | 0,996 | 38,584 | 0,996
DPID | 40,846 | 0,997 | 39,580 | 0,997 | 39,160 | 0,996 | 38,764 | 0,996
L0 | 27,994 | 0,949 | 35,405 | 0,983 | 35,138 | 0,986 | 36,183 | 0,991
d-LCI | Inf | 1,000 | 58,265 | 1,000 | Inf | 1,000 | 65,394 | 1,000
US2 400$\times$300 |  |  |  |  |  |  |  |
BIC | 33,348 | 0,987 | 33,298 | 0,987 | 33,303 | 0,987 | 33,302 | 0,987
DPID | 36,428 | 0,994 | 35,217 | 0,992 | 34,823 | 0,992 | 34,450 | 0,991
L0 | 29,184 | 0,963 | 31,289 | 0,967 | 30,647 | 0,968 | 31,394 | 0,977
d-LCI | Inf | 1,000 | 55,418 | 1,000 | Inf | 1,000 | 60,956 | 1,000
US3 300$\times$400 |  |  |  |  |  |  |  |
BIC | 43,820 | 0,999 | 43,754 | 0,999 | 43,774 | 0,999 | 43,786 | 0,999
DPID | 44,853 | 0,999 | 43,589 | 0,999 | 43,162 | 0,999 | 42,751 | 0,999
L0 | 34,756 | 0,992 | 37,140 | 0,995 | 38,087 | 0,996 | 40,237 | 0,998
d-LCI | Inf | 1,000 | 58,958 | 1,000 | Inf | 1,000 | 68,436 | 1,000
US4 300$\times$400 |  |  |  |  |  |  |  |
BIC | 40,673 | 0,996 | 40,623 | 0,996 | 40,631 | 0,996 | 40,633 | 0,996
DPID | 42,066 | 0,998 | 40,699 | 0,997 | 40,236 | 0,997 | 39,794 | 0,996
L0 | 32,944 | 0,973 | 35,973 | 0,981 | 36,453 | 0,986 | 37,836 | 0,992
d-LCI | Inf | 1,000 | 58,958 | 1,000 | Inf | 1,000 | 66,386 | 1,000
US5 400$\times$300 |  |  |  |  |  |  |  |
BIC | 30,357 | 0,979 | 30,307 | 0,979 | 30,311 | 0,979 | 30,311 | 0,979
DPID | 33,492 | 0,991 | 32,471 | 0,989 | 32,133 | 0,988 | 31,813 | 0,987
L0 | 26,807 | 0,962 | 30,944 | 0,973 | 28,911 | 0,963 | 28,551 | 0,965
d-LCI | Inf | 1,000 | 55,257 | 1,000 | Inf | 1,000 | 59,065 | 1,000
US6 300$\times$400 |  |  |  |  |  |  |  |
BIC | 29,718 | 0,978 | 29,659 | 0,978 | 29,633 | 0,978 | 29,663 | 0,978
DPID | 32,578 | 0,990 | 31,486 | 0,988 | 31,111 | 0,987 | 30,752 | 0,985
L0 | 25,322 | 0,959 | 30,666 | 0,966 | 28,724 | 0,958 | 28,022 | 0,962
d-LCI | Inf | 1,000 | 54,273 | 1,000 | Inf | 1,000 | 58,321 | 1,000
US7 273$\times$400 |  |  |  |  |  |  |  |
BIC | 40,559 | 0,995 | 40,513 | 0,995 | 40,521 | 0,995 | 40,522 | 0,995
DPID | 43,135 | 0,997 | 41,717 | 0,996 | 41,232 | 0,996 | 40,771 | 0,996
L0 | 31,832 | 0,949 | 34,364 | 0,964 | 35,604 | 0,975 | 37,697 | 0,987
d-LCI | Inf | 1,000 | 57,639 | 1,000 | Inf | 1,000 | 66,113 | 1,000
US8 322$\times$400 |  |  |  |  |  |  |  |
BIC | 31,857 | 0,989 | 31,802 | 0,989 | 31,808 | 0,989 | 31,808 | 0,989
DPID | 36,382 | 0,996 | 34,979 | 0,994 | 34,601 | 0,994 | 34,239 | 0,993
L0 | 29,562 | 0,974 | 33,051 | 0,977 | 32,378 | 0,980 | 31,828 | 0,985
d-LCI | Inf | 1,000 | 55,880 | 1,000 | Inf | 1,000 | 62,394 | 1,000
US9 322$\times$400 |  |  |  |  |  |  |  |
BIC | 28,407 | 0,984 | 28,350 | 0,984 | 28,536 | 0,984 | 28,354 | 0,984
DPID | 32,542 | 0,994 | 31,504 | 0,993 | 34,601 | 0,992 | 30,393 | 0,992
L0 | 29,562 | 0,974 | 29,504 | 0,985 | 27,361 | 0,980 | 27,214 | 0,979
d-LCI | Inf | 1,000 | 53,879 | 1,000 | Inf | 1,000 | 57,877 | 1,000
US10 400$\times$307 |  |  |  |  |  |  |  |
BIC | 35,220 | 0,992 | 35,168 | 0,992 | 35,172 | 0,992 | 35,173 | 0,992
DPID | 38,947 | 0,997 | 37,406 | 0,996 | 36,929 | 0,995 | 36,494 | 0,995
L0 | 31,031 | 0,984 | 33,760 | 0,986 | 33,227 | 0,985 | 33,744 | 0,988
d-LCI | Inf | 1,000 | 57,131 | 1,000 | Inf | 1,000 | 62,986 | 1,000
US11 400$\times$241 |  |  |  |  |  |  |  |
BIC | 34,599 | 0,995 | 34,528 | 0,995 | 34,536 | 0,995 | 34,535 | 0,995
DPID | 37,345 | 0,998 | 35,950 | 0,997 | 35,432 | 0,997 | 34,937 | 0,996
L0 | 27,631 | 0,976 | 34,346 | 0,986 | 34,143 | 0,989 | 33,189 | 0,992
d-LCI | Inf | 1,000 | 57,147 | 1,000 | Inf | 1,000 | 62,744 | 61,000
US12 400$\times$266 |  |  |  |  |  |  |  |
BIC | 35,698 | 0,993 | 35,639 | 0,993 | 035,645 | 0,993 | 35,647 | 0,993
DPID | 38,626 | 0,997 | 37,345 | 0,996 | 36,898 | 0,995 | 36,471 | 0,995
L0 | 29,883 | 0,972 | 32,902 | 0,978 | 33,091 | 0,982 | 33,776 | 0,987
d-LCI | Inf | 1,000 | 55,801 | 1,000 | Inf | 1,000 | 62,744 | 1,000
US13 310$\times$400 |  |  |  |  |  |  |  |
BIC | 38,220 | 0,996 | 38,166 | 0,996 | 38,173 | 0,996 | 38,177 | 0,996
DPID | 41,001 | 0,998 | 39,707 | 0,997 | 39,276 | 0,997 | 38,861 | 0,997
L0 | 31,593 | 0,982 | 35,761 | 0,989 | 35,898 | 0,991 | 36,263 | 0,993
d-LCI | Inf | 1,000 | 57,905 | 1,000 | Inf | 1,000 | 65,535 | 1,000
Table 6: Performance results on 13US dataset for high downscaling | :20 | :21 | :27 | :30 | :33
---|---|---|---|---|---
| PSNR | SSIM | PSNR | SSIM | PSNR | SSIM | PSNR | SSIM | PSNR | SSIM
US1 300$\times$400 |  |  |  |  |  |  |  |  |  |
BIC | 38,583 | 0,996 | 38,585 | 0,996 | 38,586 | 0,996 | 38,585 | 0,996 | 38,584 | 0,996
d-LCI | 66,485 | 1,000 | Inf | 1,000 | Inf | 1,000 | 70,845 | 1,000 | Inf | 1,000
US2 400$\times$300 |  |  |  |  |  |  |  |  |  |
BIC | 33,302 | 0,987 | 33,302 | 0,987 | 33,302 | 0,987 | 33,303 | 0,987 | 33,303 | 0,987
d-LCI | 66,347 | 1,000 | Inf | 1,000 | Inf | 66,092 | 1,000 | 1,000 | Inf | 1,000
US3 300$\times$400 |  |  |  |  |  |  |  |  |  |
BIC | 43,786 | 0,999 | 43,787 | 0,999 | 43,786 | 0,999 | 43,788 | 0,999 | 43,784 | 0,999
d-LCI | 69,721 | 1,000 | Inf | 1,000 | Inf | 1,000 | 74,415 | 1,000 | Inf | 1,000
US4 300$\times$400 |  |  |  |  |  |  |  |  |  |
BIC | 40,633 | 0,996 | 40,634 | 0,996 | 40,636 | 0,996 | 49,636 | 0,996 | 49,635 | 0,996
d-LCI | 68,054 | 1,000 | Inf | 1,000 | Inf | 1,000 | 71,936 | 1,000 | Inf | 1,000
US5 400$\times$300 |  |  |  |  |  |  |  |  |  |
BIC | 30,311 | 0,979 | 30,311 | 0,979 | 30,311 | 0,979 | 30,311 | 0,979 | 30,311 | 0,979
d-LCI | 59,853 | 1,000 | Inf | 1,000 | Inf | 1,000 | 63,768 | 1,000 | Inf | 1,000
US6 300$\times$400 |  |  |  |  |  |  |  |  |  |
BIC | 30,311 | 0,979 | 29,663 | 0,978 | 29,664 | 0,978 | 29,663 | 0,978 | 29,663 | 0,978
d-LCI | 59,135 | 1,000 | Inf | 1,000 | Inf | 1,000 | 62,833 | 1,000 | Inf | 1,000
US7 273$\times$400 |  |  |  |  |  |  |  |  |  |
BIC | 40,524 | 0,995 | 40,521 | 0,995 | 40,524 | 0,995 | 40,536 | 0,995 | 40,520 | 0,995
d-LCI | 67,318 | 1,000 | Inf | 1,000 | Inf | 1,000 | 72,151 | 1,000 | Inf | 1,000
US8 322$\times$400 |  |  |  |  |  |  |  |  |  |
BIC | 31,808 | 0,989 | 31,807 | 0,989 | 31,807 | 0,989 | 31,807 | 0,995 | 31,806 | 0,989
d-LCI | 63,539 | 1,000 | Inf | 1,000 | Inf | 1,000 | 67,653 | 1,000 | Inf | 1,000
US9 322$\times$400 |  |  |  |  |  |  |  |  |  |
BIC | 28,354 | 0,984 | 28,355 | 0,984 | 28,354 | 0,984 | 28,355 | 0,984 | 28,355 | 0,984
d-LCI | 58,651 | 1,000 | Inf | 1,000 | Inf | 1,000 | 62,536 | 1,000 | Inf | 1,000
US10 400$\times$307 |  |  |  |  |  |  |  |  |  |
BIC | 35,174 | 0,992 | 35,172 | 0,992 | 35,174 | 0,992 | 35,175 | 0,992 | 35,175 | 0,992
d-LCI | 63,924 | 1,000 | Inf | 1,000 | Inf | 1,000 | 67,994 | 1,000 | Inf | 1,000
US11 400$\times$241 |  |  |  |  |  |  |  |  |  |
BIC | 34,535 | 0,995 | 34,535 | 0,995 | 34,534 | 0,995 | 34,534 | 0,995 | 34,535 | 0,995
d-LCI | 63,445 | 1,000 | Inf | 1,000 | Inf | 1,000 | 67,442 | 1,000 | Inf | 1,000
US12 400$\times$266 |  |  |  |  |  |  |  |  |  |
BIC | 35,647 | 0,993 | 35,646 | 0,993 | 35,645 | 0,993 | 35,646 | 0,993 | 35,648 | 0,993
d-LCI | 63,760 | 1,000 | Inf | 1,000 | Inf | 1,000 | 67,127 | 1,000 | Inf | 1,000
US13 310$\times$400 |  |  |  |  |  |  |  |  |  |
BIC | 38,175 | 0,996 | 38,175 | 0,996 | 38,174 | 0,996 | 38,175 | 0,996 | 38,175 | 0,996
d-LCI | 66,347 | 1,000 | Inf | 1,000 | Inf | 1,000 | 70,960 | 1,000 | Inf | 1,000

In particular, Table 4 contains the detailed upscaling results on the NASA
dataset, for the scaling factors $s\in\\{2,4,8,16\\}$ while Tables 5-6 concern
the downscaling results obtained on 13US dataset for the scaling factors
$s\in\\{3,6,9,18\\}$ (Table 5) and for the larger factors
$s\in\\{20,21,27,30,33\\}$ (Table 6).

According to the average results, the pointwise ones corroborate the global
trend. In particular, Table 5 confirms that in downscaling d-LCI is the method
allowing the largest PSNR and SSIM, and that the method requiring the least
computation time is BIC, followed by d-LCI which is faster than DPID and L0.
Also, note that the scaling factors chosen in Table 5 are all multiple of
three but, according to Proposition 4.1 we can see the differnce between the
odd and even downscaling factors. Moreover, we can affirm that d-LCI works
fine also for large downscaling factors as reported in Table 6 where the
target images in 13US have been zoomed in up to 33 times, getting input images
whose size goes up to 13200$\times$10230\. We point out that in Table 6 the
results for DPID and L0 methods are missing because the publicly available
MatLab codes don’t work for so large scale factors. Moreover, if
$s\in\\{21,27,33\\}$ then, in accordance with Proposition 4.1, PSNR and SSIM
again reach the limit values Infinite and 1, respectively, by d-LCI while the
same does not hold for BIC applied to the same input images.

##  6 Conclusions

In the context of interpolation methods for image resizing, we present the
Lagrange–Chebyshev Interpolation (LCI) method. It is based on a non standard
mathematical modeling that leads to the application of such an optimal
Lagrange interpolation process to globally approximate the image at a
continuous scale. One of the main advantages of LCI method is its high
flexibility of working in both the scaling directions, either setting a scale
factor or giving a particular final size. Comparisons with other resizing
procedures have been reported on 5 different common datasets made of 726
images in total. The numerical experience in upscaling shows a performance
comparable with the bicubic interpolation method, while in downscaling a much
better performance is obtained with respect to all the considered comparison
methods. Moreover, in downscaling cases with odd scale factors, we estimate
the Mean Square Error (MSE) of our procedure in terms of the initial errors
present in the data and we prove that it is null in absence of noise or
artifacts in the input image. We wonder whether further improvements can be
achieved employing wavelets technique or finer approximation polynomials.
These will be the subjects of further investigations.
Code and supplementary materials

The code and the supplementary materials are openly available (to appear at a
GitHub link).

## References

  * [1] http://sipi.usc.edu/database/.
  * [2] https://www.gcc.tu-darmstadt.de/home/proj/dpid/index.en.jsp.
  * [3] http://www.shengfenghe.com/publications/.
  * [4] https://www2.eecs.berkeley.edu/Research/Projects/CS/vision/bsds/.
  * [5] https://www.cl.cam.ac.uk/~aco41/Files/Sig15UserStudyImages.html.
  * [6] National aeronautics and space administration.  NASA image gallery, 2016.
  * [7] C. Arcelli, N. Brancati, M. Frucci, G. Ramella, and G. Sanniti di Baja.  A fully automatic one-scan adaptive zooming algorithm for color images.  Signal Processing, 91(1):61–71, 2011.
  * [8] P. M. Atkinson.  Downscaling in remote sensing.  International Journal of Applied Earth Observation and Geoinformation, 22:106–114, 2013.
  * [9] L. Bos, M. Caliari, S. De Marchi, M. Vianello, and Y. Xu.  Bivariate Lagrange interpolation at the Padua points: the generating curve approach.  J. Approx. Theory, 143:15–25, 2006.
  * [10] V. Bruni, G. Ramella, and D. Vitulano.  Automatic perceptual color quantization of dermoscopic images.  VISAPP 2015, J. Braz et al. Eds., Scitepress Science and Technology Publications, 1:323–330, 2015.
  * [11] M. Caliari, S. De Marchi, and M. Vianello.  Bivariate Lagrange interpolation at the Padua points: computational aspects.  J. Comput. Appl. Math., 221:284–292, 2008.
  * [12] J. Chaki and N. Dey.  Introduction to image color feature.  In Springer-Briefs in Applied Sciences and Technology, Singapor, 2021.
  * [13] G. Chen, Zhao H., C. K. Pang, T. Li, and C. Pang.  Image scaling: How hard can it be?  IEEE Access, 7, 2019.
  * [14] S. De Marchi, W. Erb, and F. Marchetti.  Spectral filtering for the reduction of the Gibbs phenomenon for polynomial approximation methods on Lissajous curves with applications in MPI.  Dolom. Res. Notes Approx., 10:128–137, 2017.
  * [15] B. Fischer and J. Prestin.  Wavelets based on orthogonal polynomials.  Math. Comp., 66:1593–1618, 1997.
  * [16] B. Fischer and W. Themistoclakis.  Orthogonal polynomial wavelets.  Numer. Alg., 30:37–58, 2002.
  * [17] J.-B. Huang, A. Singh, and N. Ahuja.  Single image super-resolution from transforme self-exemplars.  Proc. CVPR:5197–5206, 2015.
  * [18] C. Kaethner, W. Erb, M. Ahlborg, P. Szwargulski, T. Knopp, and T. M. Buzug.  IEEE Transactions on Medical Imaging, 35(11):2476–2485, 2016.
  * [19] R. Keys.  Cubic Convolution Interpolation for Digital Image Processing.  IEEE Transactions on Acoustics, Speech, and Signal Processing, 29(6):1153–1160, 1981.
  * [20] T. Kilgore and J. Prestin.  Polynomial wavelets on the interval.  Constr. Approx., 12:95–110, 1996.
  * [21] T. Knopp and T. M. Buzug.  Magnetic Particle Imaging.  Springer-Verlag Berlin Heidelberg, Berlin, 2012.
  * [22] H. Liu, X. Xie, W.Y. Ma, and H.J. Zhang.  Automatic Browsing of Large Pictures on Mobile Devices.  In 11th ACM International Conference on Multimedia, CA, USA, 2003\. Berkeley.
  * [23] J. Liu, S. He, and R. W. H. Lau.  L0-Regularized Image Downscaling.  IEEE Trans. Image process., 27(3):1076–1085, 2018.
  * [24] T. Liu, Z. Yuan, J. Sun, J. Wang, N. Zheng, X. Tang, and H.-Y. Shum.  Learning to detect a salient object.  IEEE Trans. on Patt. Anal. Mach. Intell., 33(2):353–367, 2011.
  * [25] A. Lookingbill, J. Rogers, D. Lieb, J. Curry, and S. Thrun.  Reverse Optical Flow for Self-Supervised Adaptive Autonomous Robot Navigation.  International Journal of Computer Vision, 74(3):287–302, 2007.
  * [26] D. Martin, C. Fowlkes, D. Tal, and J. Malik.  A database of human segmented natural images and its application to evaluating segmentation algorithms and measuring ecological statistics.  In Proc. 8th Int. Conf. Computer Vision, volume 2, pages 416–423, 2001.
  * [27] G. Mastroianni and G. V. Milovanovic.  Interpolation Processes Basic Theory and Applications.  Springer Monographs in Mathematics. Springer Verlag, Berlin, 2009.
  * [28] E. H.W. Meijering, W. J. Niessen, and M. A. Viergever.  Quantitative evaluation of convolution-based methods for medical image interpolation.  Medical Image Analysis, 5:111–126, 2001.
  * [29] H. Mittal, A.C. Pandey, M. Saraswat, and et al.  A comprehensive survey of image segmentation: clustering methods, performance parameters, and benchmark datasets.  Multimed Tools Appl, 2021.
  * [30] C.H. Neetha, C. John Moses, and D. Selvathi.  Image interpolation using non-adaptive scaling algorithms for multimedia applications. a survey.  In Advances in Automation, Signal Processing, Instrumentation, and Control, Komanapalli V.L.N., Sivakumaran N., Hampannavar S. (eds), volume 700, pages 1509–1516, 2021.
  * [31] D. Occorsio and M. G. Russo.  Numerical methods for Fredholm integral equations on the square.  Appl. Math. and Comput., 218(5):2318–2333, 2011.
  * [32] D. Occorsio and W. Themistoclakis.  Uniform Weighted Approximation by Multivariate Filtered Polynomials.  Lecture Notes in Computer Science (including subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics), 11973 LNCS:86–100, 2020.
  * [33] D. Occorsio and W. Themistoclakis.  Uniform weighted approximation on the square by polynomial interpolation at Chebyshev nodes.  Appl. Math. and Comput, 385(125457), 2020.
  * [34] D. Occorsio and W. Themistoclakis.  On the filtered polynomial interpolation at Chebyshev nodes.  Appl. Numer. Math., 166:272–287, 2021.
  * [35] D. Occorsio and W. Themistoclakis.  Some remarks on filtered polynomial interpolation at Chebyshev nodes.  Dolom. Research. Notes on Approx., 14:68–84, 2021.
  * [36] A. C. Oztireli and M. Gross.  Perceptually based downscaling of images.  ACM Trans. Graph., 34(4):1–10, 2015.
  * [37] G. Plonka, D. Potts, G. Steidl, and M. Tasche.  Numerical Fourier Analysis.  Applied and Numerical Harmonic Analysis. Birkhüser Springer Nature Switzerland AG, Berlin.
  * [38] G. Ramella.  Evaluation of quality measures for color quantization.  Multimedia Tools and Applications, 2021.
  * [39] G. Ramella and G. Sanniti di Baja.  Color histogram-based image segmentation.  Lecture Notes in Computer Science, 6854:76–83, 2011.
  * [40] G. Ramella and G. Sanniti di Baja.  Image segmentation based on representative colors and region merging.  In Pattern Recognition, J. A. Carrasco-Ochoa et al Eds., volume 7914, pages 175–184, 2013.
  * [41] G. Ramella and G. Sanniti di Baja.  A new technique for color quantization based on histogram analysis and clustering.  Int. J. Patt. Recog. Artif. Intell., 27(3):1–17, 2013.
  * [42] G. Ramella and G. Sanniti di Baja.  From color quantization to image segmentation.  In Signal Imag. Techn. Internet-Based Syst. - SITIS 2016, K. Yetongnon et al. Eds., pages 798–804, 2016.
  * [43] G. Ramella and G. Sanniti di Baja.  A new method for color quantization.  In Signal Imag. Techn. Internet-Based Syst. - SITIS 2016, K. Yetongnon et al. Eds., pages 1–6, 2016.
  * [44] B. Thomee, D.A. Shamma, G. Friedland, B. Elizalde, K. Ni D. Poland, D. Borth, and L.-J. Li.  Yfcc100m: The new data in multimedia research.  Communications of the ACM, 59(2):64–73, 2016.
  * [45] L. N. Trefethen.  Approximation theorey and approximation practices.  Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2013.
  * [46] T. Vla$\breve{s}$ic, I. Rala$\breve{s}$ic, A. Tafro, and D. Ser$\breve{s}$ic.  Spline-like Chebyshev polynomial model for compressive imaging.  J. Vis. Commun. Image, 66(102731):370–378, 2020.
  * [47] N. Weber, M. Waechter, S. C. Amend, S. Guthe, and M. Goesele.  Rapid detail-preserving image downscaling.  ACM Trans. Graph., 35(6):1–6, 2016.
  * [48] Y. Xu.  Lagrange interpolation on Chebyshev points of two variables.  J. Approx. Theory, 87(2):220–238, 1996.
  * [49] Y. Xu.  Minimal cubature rules and polynomial interpolation in two variables ii.  J. Approx. Theory, 214:49–68, 2017.
  * [50] T. Yao, Y. Luo, Y. Chen, D. Yang, and L. Zhao.  Single-image Super-Resolution: A survey.  In CSPS 2018, Liang Q., Liu X., Na Z., Wang W., Mu J., Zhang B. (eds), volume 516, pages 119–125, 2020.
  * [51] M. Zhang, L. Zhang, Y. Sun, L. Feng, and W. Ma.  Auto cropping for digital photographs.  In Signal Imag. Techn. Internet-Based Syst. - SITIS 2016, K. Yetongnon et al. Eds., 2005.