# Synthetic Hyperspectral Array Video Database with Applications to Cross-
Spectral Reconstruction and Hyperspectral Video Coding

Frank Sippel  Friedrich-Alexander-Universität Erlangen-Nürnberg (FAU),
Cauerstr. 7, 91058 Erlangen, Germany  Corresponding author:
<EMAIL_ADDRESS>Jürgen Seiler  Friedrich-Alexander-Universität Erlangen-
Nürnberg (FAU), Cauerstr. 7, 91058 Erlangen, Germany  André Kaup  Friedrich-
Alexander-Universität Erlangen-Nürnberg (FAU), Cauerstr. 7, 91058 Erlangen,
Germany

###### Abstract

In this paper, a synthetic hyperspectral video database is introduced. Since
it is impossible to record ground truth hyperspectral videos, this database
offers the possibility to leverage the evaluation of algorithms in diverse
applications. For all scenes, depth maps are provided as well to yield the
position of a pixel in all spatial dimensions as well as the reflectance in
spectral dimension. Two novel algorithms for two different applications are
proposed to prove the diversity of applications that can be addressed by this
novel database. First, a cross-spectral image reconstruction algorithm is
extended to exploit the temporal correlation between two consecutive frames.
The evaluation using this hyperspectral database shows an increase in PSNR of
up to 5.6 dB dependent on the scene. Second, a hyperspectral video coder is
introduced which extends an existing hyperspectral image coder by exploiting
temporal correlation. The evaluation shows rate savings of up to 10% depending
on the scene. The novel hyperspectral video database and source code is
available at https://github.com/FAU-LMS/HyViD for use by the research
community.

###### doi:

https://doi.org/10.1364/JOSAA.479552

††journal: josaa

##  1 Introduction

Classical cameras record a scene using three color channels, namely, the red
(R), green (G) and blue (B) channel. To achieve this, the incoming light
spectrum is filtered by corresponding filters, e.g., a filter for the red
component lets pass the wavelengths from 600 nm to 700 nm. The RGB filters are
approximating the properties of the different cones of the human eye. Thus, if
such an image is shown to a human on a proper display, all perceivable
information is recorded and displayed.

For many classification tasks, it is useful to record even more channels, for
example also in the infrared (IR) area of the spectrum. Thus, the concept of
RGB imaging is extended to multispectral imaging, where typically six to 16
different bands of the scene are recorded. Usually, the filters of
multispectral cameras are spectrally distributed depending on the application.
Therefore, there is no structure regarding the shape of filters or the
location of filters in the spectrum.

In contrast, the filters of hyperspectral cameras are typically equidistantly
distributed in the wavelength area of interest. Since the employed bandpass
filters are narrowband as well, the goal of a hyperspectral camera is to
sample the light spectrum in the desired wavelength area for each pixel. While
there are no sharp definitions of multispectral and hyperspectral, these
definitions are inline with the definitions given by Hagen et al.[1], which
state that most authors distinguish between contiguous (hyperspectral) versus
spaced spectral bands (multispectral).

Hyperspectral imaging has become increasingly popular for diverse
classification tasks. For example, different types of plastics can be
separated using hyperspectral cameras[2], which is an essential building block
for recycling pipelines. Moreover, these cameras are employed in agriculture
to measure the plant health[3] in order to distribute the proper amount of
water and fertilizer. Other applications include forensics[4] where
hyperspectral cameras can analyse the scene in a non-invasive way, object
tracking[5], and tumour detection in medicine[6]. Usually, for most of these
applications, it is beneficial if the hyperspectral camera is able to record
videos.

Figure 1: A hyperspectral 3D datacube with two spatial dimensions $x$ and $y$ and one spectral dimension $\lambda$. Table 1: Overview over different hyperspectral databases. Note that the 1890 images of the proposed database are put together by 7 scenes, each with 30 frames, rendered from a camera array with 9 cameras. Database | Spatial res. | #Channels | Wavel. area (nm) | #Images | Videos | Rerendering | Cam. array | Depth
---|---|---|---|---|---|---|---|---
CAVE[7] | $512\times 512$ | 31 | 400 - 700 | 32 | ✗ | ✗ | ✗ | ✗
UGR[8] | $1000\times 900$ | 61 | 400 - 1000 | 14 | ✗ | ✗ | ✗ | ✗
Hordley[9] | $\sim 336\times 271$ | 31 | 400 - 700 | 23 | ✗ | ✗ | ✗ | ✗
Le Moan[10] | $500\times 500$ | 160 | 410 - 1000 | 9 | ✗ | ✗ | ✗ | ✗
Foster[11] | $1018\times 1339$ | 33 | 400 - 720 | 8 | ✗ | ✗ | ✗ | ✗
Chakrabarti[12] | $1392\times 1040$ | 31 | 420 - 720 | 50 | ✗ | ✗ | ✗ | ✗
BGU[13] | $1392\times 1300$ | 244 | 400 - 1000 | 201 | ✗ | ✗ | ✗ | ✗
Mian[14] | $752\times 480$ | 33 | 400 - 720 | 31 | ✓ | ✗ | ✗ | ✗
Proposed | $1600\times 1200$ | 31 | 400 - 700 | 1890 | ✓ | ✓ | ✓ | ✓

There are several techniques to record a hyperspectral image. The essential
challenge is to record a 3D datacube using only 2D sensors. This hyperspectral
3D datacube is depicted in Fig. 1. One can imagine combining different layers
or single data elements to fit onto a 2D sensor. Some of the following
techniques drop some of these data elements and need a reconstruction process
afterwards. The first class of recording techniques are scanning devices.
Here, one dimension of this 3D datacube is unfolded across the time dimension.
For example, pushbroom sensors[15] are scanning one spatial dimension and the
spectral dimension using a 2D sensor. The second spatial dimension is scanned
over time by moving the scene or camera accordingly. Another example of
scanning techniques are filter wheels[16], where the spectral dimension is
scanned by spinning a wheel, such that different filters are in front of the
sensor. Thus, this technique unfolds the spectral dimension over time.
Obviously, these techniques are not able to capture hyperspectral videos,
since the 3D datacube is scanned over time and thus not fully available at
each time instance.

Hyperspectral snapshot cameras, on the other hand, are able to capture
hyperspectral videos. Several ideas exist how to build a snapshot
hyperspectral camera. There are concepts using dispersion and fiber
bundles[17] or a 2D dispersion pattern within a computed tomography
framework[18]. Thus, a reconstruction process is mandatory. Both of these
techniques suffer from a very low spatial resolution. Another technique is
multispectral beamsplitting[19], where the incoming light beam is split into
the corresponding bands using special mirrors. The problem here is, that the
amount of photons for other bands is also drastically reduced by each mirror.
For a large number of bands, this technique does not work anymore.
Multispectral filter arrays[20] have a similar idea as a classical Bayer
pattern by placing filters periodically over each pixel of a sensor.
Subsequently, a demosaicing algorithm is necessary to yield a full
hyperspectral image. However, depending on the number of bands, the underlying
physical resolution for each channel drops significantly using multispectral
filter arrays. Especially, when thinking of hyperspectral imaging which can
have up to several hundred bands, the pixels with the same filter are
positioned far apart. Multi-aperture filtered cameras[21] have a similar idea,
but place filters over whole regions of the sensor. Thus, the filters are not
periodically repeated. Here, also a reconstruction process is necessary.
Imaging techniques based on compressed sensing[22] are also discussed. For
these cameras, a coded aperture followed by a dispersion element is used to
record different spectral mixtures of different pixels. Afterwards, compressed
sensing is used to reconstruct the datacube. Finally, works on camera
arrays[23] are currently active as well. Since the cameras are spatially
distributed, a reconstruction pipeline yielding a consistent datacube is
essential. Two examples for the necessity of a reconstruction pipeline for
high-resolution hyperspectral snapshot cameras are shown in Fig. 2.

Figure 2: Two cases of missing ground-truth pixels of snapshot $3\times 3$
multispectral/hyperspectral cameras. On the left, the warped channels of a
camera array are depicted, where pixels are missing due to occlusions. On the
right, the channels of a multispectral filter array sensor are shown. Missing
pixels are shown in white.

From this short overview of hyperspectral cameras, it becomes clear that it is
challenging to get ground-truth high-resolution hyperspectral video data.
Ground-truth data has several requirements. First, the scene itself should be
as realistic as possible. Of course, when recording a scene in reality, this
is fulfilled by definition, while it is more of a challenge for synthetic
data. For synthetic data, the world, including all physical concepts like
light, needs to be modeled properly to get as realistic data as possible.
Second, ground-truth data needs a high resolution in every dimension to be
able to sample in any direction. While scanning techniques can deliver this as
long as the scene is static, some hyperspectral snapshot cameras either record
a very low spatial or spectral resolution or have low underlying spatial or
spectral resolution. Third, the amount of artefacts and postprocessing should
be as low as possible. These artefacts include noise, blur, reconstruction
errors, exposure problems, chromatic aberrations and many more. Here,
recording a scene using existing hyperspectral imaging techniques fails to
fulfill this requirement. For a non-moving scene a near ground-truth
hyperspectral image can be generated using scanning techniques. However, for
example for a filter wheel, due to different behaviours of lens-filter
combinations, the focus and exposure time have to be adjusted, which leads to
different behaviours regarding optics and noise. Snapshot hyperspectral
cameras most often need a reconstruction process, which cannot recover the
perfect datacube. Thus, again artefacts may be introduced. Fourth, boundary
conditions should be given as exact as possible. This includes the used
camera(s) including all specifications like sensor size and resolution, as
well as the aperture and focal length of the lens. Furthermore, properties of
the bandpass filters and light are also very important. Since there are always
slight deviations during manufacturing processes and aging phenomena, it is
impossible to properly fulfill this requirement when recording a scene
hyperspectrally in reality.

As one can see, there is a strong need for ground-truth hyperspectral video
data, for example also as reference data to measure the performance of
snapshot reconstruction algorithms. There are a number of successful
hyperspectral image databases, which are summarized in Table 1. One of the
most famous examples is the CAVE database[7]. This database contains 32 scenes
depicting different objects like balloons, persons, pepper and color patches
usually with a black background. The images have a resolution of $512\times
512$ pixels in the range from 400 nm to 700 nm in 10 nm steps, resulting in 31
hyperspectral channels. The CAVE database does not contain any natural scene.
The UGR database[8] contains 14 images with a resolution of $1000\times 900$
pixels. This database depicts mostly urban scenes including trees, cars and
buildings. The scenes were recorded from 400 nm to 1000 nm in 10 nm steps,
resulting in 61 hyperspectral channels. Unfortunately, some of the bands are
corrupted by artefacts. The database of Hordley et al.[9] contains 23 images
of packagings and color charts and thus does not hold any natural spectra. The
resolution of the images differ. The packagings and color charts were recorded
from 400 nm to 700 nm in 10 nm steps, resulting in 31 hyperspectral channels.
Le Moan et al.[10] created a database of 9 different textures including
textile, wood and skin from 410 nm to 1000 nm in 160 hyperspectral channels.
The textures have a resolution of $500\times 500$. The next database of Foster
et al.[11] contains 8 hyperspectral images of natural and urban scenes. These
were recorded in the range from 400 nm to 720 nm in 10 nm steps, resulting in
33 hyperspectral channels. The images have a resolution of $1018\times 1339$.
Chakrabarti et al.[12] recorded a hyperspectral image database of 50
hyperspectral images with a resolution of $1392\times 1040$. For each scene,
31 hyperspectral channels were recorded from 420 nm to 720 nm in 10 nm steps.
This database mainly contains urban and indoor scenes. The BGU database[13] of
Arad et al. contains over 200 hyperspectral images of natural, urban and
indoor scenes at a resolution of $1392\times 1300$ pixels. The database was
recorded from 400 nm to 1000 nm at 1.25 nm steps. Thus, the BGU database
offers the highest spectral resolution, one of the highest spatial resolutions
and the most diversity in scenes due to the sheer amount of different scenes.
Finally, a hyperspectral video database[14] containing a single video of a lab
scene with a spatial resolution of $752\times 480$ exists. Here, the spectrum
from 400 nm to 720 nm was recorded using 10 nm steps, thus resulting in 33
hyperspectral channels.

This paper introduces a novel synthetic hyperspectral video database rendered
using a camera array. This camera array consists of nine cameras arranged in a
three times three grid. Thus, the scenes are rendered from multiple positions
for wavelengths between 400 nm and 700 nm in 10 nm steps, resulting in 31
hyperspectral channels for each camera. Moreover, the depth map is generated
as well for every frame and camera. Since a high resolution representation of
spatial positions in three dimensions as well as the spectral dimension is
available, many applications can be addressed by this database. These
applications do not necessarily have to be in a hyperspectral context, but
also range to grayscale and RGB image processing tasks as well as fusion
problems, especially regarding depth sensors. To the best of our knowledge,
this paper is the first to introduce a synthetic hyperspectral video database.

The remainder of this paper is organized as follows. First, the techniques and
properties of this hyperspectral video database are described in detail in
Section 2. Afterwards, two applications of this database are considered in
Section 3 by improving the corresponding state of the art. Finally, the paper
is summarized in Section 4.

##  2 Database

The hyperspectral video database was rendered using the open-source software
Blender[24], which provides tools for RGB rendering. Therefore, the first
properties and techniques that are described in this section are related to
Blender. The basic idea is to render each wavelength for each frame image by
image in grayscale. For that, the textures need to be adjusted to the
wavelength as well as the light source. Of course, the scene can be rendered
from multiple cameras within a camera array.

###  2.1 Renderer

Figure 3: The pipeline of the hyperspectral renderer. Essentially, all
combinations of channel (cnl), frame (frm) and camera (cmr) are traversed and
rendered. The red items in the top row indicate an exemplary state, i.e.,
currently the textures at 640 nm and the animations at frame 10 of 30 viewed
from camera 6 are rendered.

To extend Blender to a hyperspectral image renderer, each wavelength is
rendered as grayscale image. Since the wavelengths from 400 nm to 700 nm were
rendered in 10 nm steps, one hyperspectral image consists of 31 rendered
grayscale images. For each wavelength, the images of materials as well as the
intensity of the light source were changed programmatically through Python.
More details on textures and light sources will be presented below. This
pipeline is explained in detail in Fig. 3.

Blender has two different rendering engines, namely, Eevee and Cycles. While
Eevee is a simple rasterization engine, Cycles implements path tracing[25] and
thus yields much more realistic images. Consequently, Cycles was chosen as
rendering engine with branched path tracing as integrator. The number of
samples per pixel was set to 128.

Since path tracing nearly always yields rendering noise in complex scenes, a
post-processing denoiser was employed. Considering that Blender directly
implements Non-Local Means (NLM) denoising[26], NLM was used as denoiser after
the grayscale images were rendered.

The bitdepth of the final output was set to eight to allow for a smaller total
database size. Furthermore, each grayscale image is saved as PNG.

###  2.2 Hyperspectral Textures and Materials

An essential part for setting up a realistic scene in Blender are textures to
give surfaces of objects color and details without increasing the complexity
of geometry. To give a surface the desired part of a texture, coordinates in
images are mapped to coordinates of surfaces. Since a hyperspectral scene also
needs hyperspectral textures, the possibilities are either to record them,
which requires a high resolution hyperspectral camera, or extract them from
other databases. The latter was picked for this database. To maintain a
consistent appearance of the textures, the BGU database[13] was chosen to
deliver hyperspectral textures. Though this database contains wavelengths from
400 nm to 1000 nm in 1.25 nm increments, our novel database only contains
wavelengths from 400 nm to 700 nm in 10 nm steps, to keep the rendering time
feasible and the database size acceptable. The BGU database contains over 200
images of diverse natural, urban, and indoor scenes. Thus, different textures
like grass, wood, flower, leaves, metal, plastic, different types of stones,
parts of cars and many more could be extracted from this database. The sky was
also extracted from this database using the image hill_0325-1228, which shows
a huge portion of a cloud-free sky. This texture was then used for the world
material in Blender. Since Blender is only an RGB renderer natively,
hyperspectral textures consist of one grayscale image for each wavelength.
Thus, depending on the wavelength that is currently rendered, the textures
need to be exchanged.

The scenes have 14 to 34 different textures depending on the complexity. For
example, an outdoor scene needs rather less texture, since there is mainly
dirt, grass, wood and leaves. Then, for these types of materials, different
textures for leaves and wood are extracted from the BGU database. On the other
hand, indoor scenes have typically more textures, since there are more smaller
things lying around the scene, which all need a texture. Of course, one can
reuse some textures for some objects, but an interesting indoor scene also
needs diverse textures even for the same object.

Figure 4: The procedure of extracting hyperspectral textures. The RGB image is
taken to spot good textures. The extraction is done by warping all
hyperspectral channels using a calculated homography matrix.

The extraction procedure is shown in Fig. 4. First, four points spanning a
rectangle are selected. Afterwards, a homography matrix is calculated, which
translates these world coordinates in the image plane to coordinates
$\left((0,0),(0,H),(W,H),(W,0)\right)$, where $W$ and $H$ represent the width
and the height of the resulting texture. This has to be defined by the user.
From these four pairs of points
$(x_{i}^{\text{src}},y_{i}^{\text{src}})\leftrightarrow(x_{i}^{\text{dst}},y_{i}^{\text{dst}})$,
the homography matrix

$\bm{H}=\begin{bmatrix}h_{11}&h_{12}&h_{13}\\\ h_{21}&h_{22}&h_{23}\\\
h_{31}&h_{32}&h_{33}\end{bmatrix}$ (1)

is found by solving the least-squares problem

$\begin{split}\operatorname*{\arg\\!\min}_{\bm{H}}\sum_{i=1}^{4}&\left(x_{i}^{\text{src}}-\frac{h_{11}y_{i}^{\text{dst}}+h_{12}x_{i}^{\text{dst}}+h_{13}}{h_{31}x_{i}^{\text{dst}}+h_{32}y_{i}^{\text{dst}}+h_{33}}\right)+\\\
&\left(y_{i}^{\text{src}}-\frac{h_{21}y_{i}^{\text{dst}}+h_{22}x_{i}^{\text{dst}}+h_{23}}{h_{31}x_{i}^{\text{dst}}+h_{32}y_{i}^{\text{dst}}+h_{33}}\right).\end{split}$
(2)

Since a homography matrix has eight degrees of freedom, the final matrix is
scaled such that $h_{33}=1$. Then, the position of other points in the
destination image plane can be calculated by

$\begin{bmatrix}x^{\text{dst}}\\\ y^{\text{dst}}\\\
1\end{bmatrix}=\bm{H}\begin{bmatrix}x^{\text{src}}\\\ y^{\text{src}}\\\
1\end{bmatrix}.$ (3)

This is necessary to warp the texture from the source image plane to the
destination image texture plane. This warping procedure is done for every
hyperspectral channel. Afterwards, the hyperspectral channels are saved
individually as lossless grayscale images.

After the textures are extracted, corresponding materials can be set up in
Blender. The texture of the material is changed depending on which channel is
currently rendered. Moreover, a lot of properties can be changed in Blender,
which influences how light is reflected or absorbed. For example, the
roughness can be changed as well as how specular and metallic a material is,
e.g., a shiny car paint is typically very little rough but very metallic,
while a street is very rough and less shiny. Moreover, also the transparency
of a texture can be set, which is very important for windows.

###  2.3 Light Sources

For all outdoor scenes or scenes with windows, the light source is the sun. In
bright daylight the sun can reach a color temperature up to 6400K. The color
temperature of a light source defines its emitting light spectrum. The
spectrum is based on the radiation of an ideal black body and can be
calculated by Planck’s law. For the single indoor scene without windows, a
lamp was mounted into that room with a color temperature of 3200K, which is a
typical value for studio lamps. The spectra of these two light sources is
shown in Fig. 5 and were calculated by Planck’s law. For each wavelength, the
intensity of the sun or the indoor lamp in the scene is set according to this
plot. Furthermore, the hardness and light strength differs from scene to scene
and is also dependent on the angle of the sun.

Figure 5: Used light spectra for the sun (6400K) and the indoor lamp (3200K).

###  2.4 Cameras and Camera Array

All scenes are rendered from nine cameras arranged in a $3\times 3$ grid. This
mimicks a camera array from Genser et al.[23] and is shown in Fig. 6.
Therefore, the specifications are very similar. The baseline between
neighboring cameras is 40 mm. The cameras itself are very close to common
industrial cameras like a Basler acA1600-60gm, which are also used for the
camera array in Fig. 6. Thus, the image sensor size is set to 7.2 mm times 5.4
mm and the resolution of the images is 1600 pixels times 1200 pixels. The
focal length is different from the original camera array and is set to 6 mm
focal length to capture a large portion of the scene. The purpose of rendering
the scene from this camera array is, for example, to be able to simulate the
multispectral camera shown in Fig. 6. Since this multispectral camera array
only records one spectral area with a specific transmission curve per camera,
one has to simulate the corresponding filters by calculating a weighted sum of
hyperspectral images for each camera. Thus, in contrast to the camera array in
Genser et al.[23], the synthetic data contains all 31 hyperspectral channels
for every camera. If the application just demands for a single hyperspectral
view of the scene, only the center camera can be extracted, since each camera
has its own folder in the database.

Figure 6: A multispectral camera with nine channels on the left[23] with a
sketch of dimensions and indices of cameras on the right. Note that this view
is from the front, as well.

###  2.5 Depth Maps

Many applications also rely on depth data provided by a depth sensor like
RADARs or LIDARs. Moreover, some algorithms also estimate the depth as
intermediate step, for example, to register multiple views on top of each
other. Therefore, depth data for all scenes, cameras, and frames is provided.
The depth maps are stored in the EXR format as 32 bit floats for a high
resolution. The depth maps contain the physical depth in meters. Consequently,
the disparity between cameras is calculated by

$d=\frac{b\cdot f}{0pt\cdot s},$ (4)

where $d$ is the disparity of corresponding pixels between two cameras, $b$ is
the baseline between the cameras, e.g. 40 mm for directly neighboring cameras
in the camera array, $f$ is the focal length, thus 6 mm, $0pt$ is the depth
itself, and $s$ is the size of a single pixel and calculates to $s=7.2\
\text{mm}/1600=5.4\ \text{mm}/1200=4.5\cdot 10^{-3}\ \text{mm}$. Note that the
depth is provided in meters while all other properties are given in
millimeters. Depth maps can be used for example to evaluate depth estimation
algorithms or to warp different views of the camera array to different
positions.

###  2.6 Motion

In all of the scenes, the camera array moves along a specified path. Depending
on the scene, this is realized using different approaches. In the first
method, the camera array is moved along the scene manually and this trajectory
is recorded. After recording this movement, the paths are smoothed out, since
a manual camera movement is typically quite shaky. In the end, each frame is a
keyframe for the movement and stores the position and rotation of the camera
array. The second way to define movement in this database is to define the
starting and endpoint position and rotation and smoothly interpolation between
them for intermediate frames. Thus, there are only two keyframes set and the
intermediate points are determined using Bezier interpolation. This approach
is also performed for moving objects.

###  2.7 Scenes

Figure 7: The workflow to create a hyperspectral scene. This general workflow
can be executed in any rendering framework. Figure 8: The center views of the
seven different synthetic hyperspectral scenes: frame 0 of "family house",
frame 0 of "medieval seaport", frame 25 of "city", frame 0 of "outdoor", frame
10 of "forest", frame 0 of "indoor" and frame 14 of "lab". On the left, an RGB
image using standard CIE 1931 color curves can be seen. Next to it, all 31
hyperspectral channels of the corresponding image and the depth map are
depicted.

The whole workflow for setting up a 3D hyperspectral scene is summarized in
Fig. 7. In total, seven different scenes were rendered. These scenes are
depicted in Fig. 8, where the RGB images are calculated using standard CIE
1931 color curves. Furthermore, all 31 hyperspectral channels of one frame are
shown as well as the depth maps corresponding to the respective scene and
frame. The camera array moves through each scene for 30 frames.

There are five outdoor scenes. Family house represents a mixture of urban and
rural objects, medieval seaport is a view of a more ancient city incorporating
a lot of stone and wood, and city is a newer city including a moving car,
normal streets, sidewalks, modern houses and trash bins. Outdoor represents a
rural scene from an agricultural way and forest is a similar scene without any
man-made objects and a lot of trees. Indoor depicts a room including furniture
and decoration elements, but still with the light source being the sun shining
through windows. Finally, lab is illuminated by a halogen lamp and depicts a
laboratory setup with selective elements on a table similar to the real-world
image taken in Genser et al.[23].

###  2.8 Limitations

Note that this database still contains only synthetic hyperspectral data.
Thus, it is nearly clean from any noise or artefacts and ideal cameras are
assumed as well as an ideal camera array is assumed. An array-related real-
world artefact would be that the cameras are not perfectly aligned to each
other as well as the sensors within the camera housing are not perfectly
fitted. Therefore, in real-world a calibration procedure between different
cameras would always be necessary, e.g., to fulfill the epipolar constraint.
Moreover, a lot of artefacts of the filter-lens-camera combinations can occur.
First, filters have the tendency to change optical paths slightly dependent on
the wavelength. Second, lenses introduce artefacts like optical veiling
glare[27], lens flares[28], wavelength-dependent chromatic aberration[29],
star bursts and many more. Finally, cameras always introduce different kind of
noise, e.g., shot noise, thermal noise and electronic noise[30]. All of these
artefacts are not part of this database and most of them are difficult to
simulate. Thus, if one trains or evaluates algorithms based on this data, one
has to be careful, since the data of this database deviates from real-world
image acquisitions. However, depending on the desired application, noise or
artefacts could be added to clean data. Moreover, while the absolute
performance is wrongly estimated due to missing artefacts, the relative
performance between two algorithms might still be possible assuming that both
algorithms are affected by noise and artefacts equally.

##  3 Applications

In this section, an overview over applications that can benefit from this
dataset is given. First, two applications are presented in detail and improved
which is verified using the novel database. In the first application, videos
from all cameras and depth maps from the center camera are used. The second
application uses just the hyperspectral video of the center camera. In the
end, several other applications are described briefly.

###  3.1 Cross Spectral Reconstruction

The first application that is covered in detail is cross-spectral image
reconstruction. Consider imaging a scene with the camera array shown in Fig.
6. The first step to create a consistent multispectral datacube is a cross-
spectral disparity estimation, i.e., finding the correspondence of which pixel
in one camera belongs to which pixel in the other camera. Since the final
reconstruction process is evaluated here, the ground-truth depth maps from the
database can be used. For that, the depths needs to be converted to
disparities using (4). Then, this disparity map can be used to warp the
peripheral view to the center view. Due to occlusions, there are pixels
missing which need to be reconstructed. Fortunately, the center view is always
fully available, since it is not warped to any other position. Thus, the
center view can be used as reference channel for all peripheral views.
However, this center view records a different part of the spectrum than the
peripheral views. Note that the procedure to simulate a multispectral camera
array is described in Section 3.1.3. Therefore, one has to find a relationship
between the center view and the peripheral view that has to be reconstructed.
This whole problem is depicted in Fig. 9. The peripheral view that has to be
reconstructed will be called distorted view in the following description.

Figure 9: The basic problem of cross spectral reconstruction. On the left side
the fully available center view is depicted, on the right side the warped
peripheral view is shown. Due to occlusions, some pixels are missing in the
peripheral view, which are shown in red.

####  3.1.1 State of the Art

For a single multispectral image, Non-local Cross-Spectral Reconstruction
(NOCS) in [31] is a state-of-the-art approach. NOCS first finds similar blocks
using the fully available reference images. Block means the set of pixels of
usually a square area around a center pixel. In the case of the multispectral
camera array, this would be just the center view. For that, the $l_{2}$-norm
between two blocks within the reference image is minimized

$d_{\text{NOCS}}(\bm{x},\bm{y})=||\bm{B}^{\text{R}}(\bm{x})-\bm{B}^{\text{R}}(\bm{y})||_{2},$
(5)

where $\bm{B}^{R}(\bm{x})$ describes the reference block at coordinate
$\bm{x}$. Of course, this whole block matching procedure only has to be done
for pixels that are missing in the distorted view. For missing pixel position
$\bm{x}$, the $B$ best matches are used to stack the corresponding pixel
values of the reference view on top of each other, resulting in the vector
$\bm{r}(\bm{x})$. The same procedure is done for the distorted image,
resulting in the distorted vector $\bm{d}(\bm{x})$. Of course, this vector
contains missing pixels. Moreover, since a block has distance 0 with itself,
the pixel to reconstruct is always at the first entry of the vector.

These two vectors are then used to build a linear regression model

$\bm{d}(\bm{x})=\alpha(\bm{x})\cdot\bm{r}(\bm{x})+\beta(\bm{x}),$ (6)

where $\alpha(\bm{x})$ and $\beta(\bm{x})$ are the parameters of the linear
regression and thus the variables that need to be estimated for every missing
pixel. They are found by minimizing the $l_{2}$-norm of the difference between
known references of the distorted vector $\tilde{\bm{d}}(\bm{x})$ and its
prediction from the reference vector $\tilde{\bm{r}}(\bm{x})$ at the
corresponding positions

$\hat{\alpha}(\bm{x}),\hat{\beta}(\bm{x})=\operatorname*{\arg\\!\min}_{\alpha(\bm{x}),\beta(\bm{x})}||\alpha(\bm{x})\cdot\tilde{\bm{r}}(\bm{x})+\beta(\bm{x})-\tilde{\bm{d}}(\bm{x})||_{2}^{2}.$
(7)

These parameters are found in closed-form[31]. Afterwards, this model is used
to predict the missing pixel, which is the first element in the distorted
vector. For that, the value of the first element in the reference vector is
put into the model and the resulting value is the predicted reconstructed
pixel value

$\bm{d}_{1}(\bm{x})=\hat{\alpha}(\bm{x})\cdot\bm{r}_{1}(\bm{x})+\hat{\beta}(\bm{x}).$
(8)

Of course, not all missing pixels are reconstructed at once, but an iterative
procedure is carried out. Always 10% (but at least one) of the remaining
pixels to reconstruct that have the highest number of non-missing pixels of
the distorted view in their vector are reconstructed. Afterwards, all vectors
that contain a newly reconstructed pixel are updated and the loop starts from
the beginning. Thus, reconstructed pixels will influence the remaining pixels
that need to be reconstructed.

####  3.1.2 Novel Approach

We propose to extend this method by also using the previous frame and
therefore exploit temporal correlation to enhance NOCS. This novel method is
called Temporal Cross Spectral Reconstruction (TNOCS). To exploit temporal
correlation, the block matching procedure is also executed on the previous
frame. Therefore, the block matching procedure needs to consider that similar
blocks might be in different frames

$d_{\text{TNOCS}}(\bm{x},\bm{y},t_{d})=||\bm{B}^{\text{R}}_{t}(\bm{x})-\bm{B}^{\text{R}}_{t-t_{d}}(\bm{y})||_{2},$
(9)

where $\bm{B}^{R}_{t}(\bm{x})$ is a block in the reference video at spatial
position $\bm{x}$ and current frame $t$, and $t_{d}$ is the difference to the
respective frame. Therefore, when a similar block for $\bm{x}$ is searched,
the position and the frame of the other block can be varied. Since two
consecutive frames are usually highly correlated but still depict the scene
from a slightly different angle, only the previous and the current frame are
searched for similar blocks. Thus, $t_{d}$ can only be in the set of
$\\{0,1\\}$. The remaining formulas are the same, since the pixel at position
$\bm{x}$ still needs to be reconstructed and the vectors contain pairs of
reference pixels and the corresponding pixels from the distorted video.

For an improvement over NOCS, it is essential that the camera or object in
scene moves from one frame to another frame, since otherwise the pixels to be
reconstructed would be exactly the same and the temporal correlation can not
be properly exploited. Note that only non-reconstructed pixels of the previous
frame are used for block matching, since the error in reconstruction should
not propagate through all frames. Otherwise, this would lead to severe
artifacts after some frames. Of course, the block matching is also executed on
the current frame to further exploit spatial correlation. This idea is
summarized in the Fig. 10.

Figure 10: The proposed TNOCS searches similar blocks in the current frame as
well as the previous frame to build the linear regression model. The blue box
depicts the block for which similar blocks are searched, the red blocks are
the best matches. The images depict cutouts of two consecutive frames of the
lab scene.

####  3.1.3 Evaluation

To validate that this idea works well, our novel hyperspectral video database
is employed. For that, multispectral video data has to be created out of the
hyperspectral data. To achieve this, an imaging pipeline needs to be
simulated. A single pixel of the $i$-th channel of a multispectral image is
recorded by[32]

$\bm{c}_{i}=\int
q(\lambda)r(\lambda)\bm{f}_{i}(\lambda)o(\lambda)m(\lambda)\text{d}\lambda,$
(10)

where $q(\lambda)$ is the spectrum of the light source, $r(\lambda)$ is the
reflectance spectrum of the material, $\bm{f}_{i}(\lambda)$ is the transfer
function of the $i$-th filter corresponding to the $i$-th multispectral
channel, $o(\lambda)$ the spectral permeability of the camera lens and
$m(\lambda)$ is the spectral response of the camera sensor. In this case, the
database already yields a sampled version of the product of the light source
and reflectance $s(\lambda)=q(\lambda)r(\lambda)$ as images. The transfer
function of the camera lens and camera sensor are assumed to be perfect in the
covered wavelength area from 400 nm to 700 nm. Since the images are obviously
sampled in the spectral dimension, this equation needs to be sampled as well.
Therefore, a single hyperspectral pixel from our database can be described by
$\bm{s}=[s(\lambda_{1}),s(\lambda_{2}),\dots,s(\lambda_{31})]^{\text{T}}$.
When stacking the sampled filters on top of each other into the filter matrix
$\bm{F}=[\bm{f}_{i}(\lambda_{1}),\bm{f}_{i}(\lambda_{2}),\dots,\bm{f}_{i}(\lambda_{31})]$,
the imaging process for a single pixel can be written as

$\bm{c}=\bm{F}\bm{s}.$ (11)

In this case, the filter matrix contains sampled values of real filters used
for the camera array. Namely, six bandpass filters with bandwidth 50 nm at
wavelengths 425 nm, 450 nm, 500 nm, 550 nm, 600 nm and 650 nm, and a red
filter, a green filter and a blue filter to compose an RGB image. Then, the
data to reconstruct is created by calculating the multispectral data for the
corresponding camera. The ground truth is calculated by doing the same
calculations with the middle camera. Moreover, the warped images are
calculated using the ground-truth depth maps provided by the database.

Table 2: Evaluation in terms of PSNR of the novel temporal cross spectral reconstruction algorithm (TNOCS) against the version without exploiting temporal correlation (NOCS). | NOCS | TNOCS
---|---|---
Family house | 30.83 dB | 36.43 dB
Medieval seaport | 40.06 dB | 44.87 dB
City | 43.61 dB | 45.39 dB
Outdoor | 28.94 dB | 34.05 dB
Forest | 27.03 dB | 30.67 dB
Indoor | 39.33 dB | 41.15 dB
Lab | 35.40 dB | 36.26 dB
Average | 35.03 dB | 38.40 dB

Using this multispectral database, the novel TNOCS algorithm can be evaluated
against its non-temporal version NOCS. The results are summarized in Table 2.
From the table it gets obvious that exploiting temporal correlation by just
searching for similar blocks in the previous frame is enhancing the
performance significantly. Especially for scenes containing a lot of leaves,
which leads to many small areas that are missing, the difference is enormous.
For other scene like lab order city, the improvement is smaller but still
noticeable.

Without our novel hyperspectral video database, it would be impossible to
accurately validate the performance gain of exploiting temporal correlation.
As shown, it is also possible to create nearly any other data that lies in the
wavelength area from 400 nm to 700 nm.

###  3.2 Hyperspectral Video Coding

The second application addressed is hyperspectral video coding. The need for
hyperspectral image and video coding is basically already motivated by this
database. With just seven scenes, nine cameras and 30 frames, this database
already has a size bigger than 60 GB saved as single grayscale PNGs. When
looking at raw values, this number increases even more. For this hyperspectral
database, saving all images without any compression and eight bits per
grayscale pixel results in $7\cdot 9\cdot 30\cdot 31\cdot 1600\cdot 1200\
\text{bytes}\approx 105\ \text{GB}$. This enormous amount of data shall be
compressed, since there is a lot of correlation. Of course, there is temporal
correlation, which can be exploited. But also the different hyperspectral
channels are highly correlated which can be exploited by a coder. In the
following, the term encoding means producing a compressed bitstream out of a
(hyperspectral) video. In contrast, the term decoding means transforming this
bitstream back to the (hyperspectral) video. For lossy coders, this decoded
video deviates from the original video.

####  3.2.1 State of the Art

The presented method is based on a multispectral and hyperspectral image coder
named Pel-Recursive Inter-Band Prediction (PRIBP)[33]. Before presenting the
novel inter component of this coder, the intra hyperspectral coder is reviewed
first. The basic idea of[33] is to first code three channels, in this paper
RGB channels, with a standard coder. Then, the three decoded channels serve as
reference channels for the next channels to code. The decoded channels need to
be used as reference channels, because the other end of the communication
channel only has the decoded images, which slightly deviates from the original
images depending on the quality settings, available. Thus, to avoid a drift
between encoder and decoder, the original channels cannot be used as reference
channels. For that, more intra modes have been implemented into the High
Efficiency Video Coding (HEVC) codec[34]. The prediction is done on a
transform block level using previously predicted pixels and reconstructed
pixels from adjacent blocks. To code the center pixel $b$ of the current block
$\bm{b}$ of the channel to code, first the best reference channel block
$\bm{r}$ out of the current three reference channels is searched by picking
the one with the maximum correlation with the already predicted pixels

$\rho(\tilde{\bm{r}},\tilde{\bm{b}})=\frac{\sum_{i=1}^{N}\left(\tilde{\bm{r}}_{i}-\mathbb{E}\\{\tilde{\bm{r}}\\}\right)\left(\tilde{\bm{b}}_{i}-\mathbb{E}\\{\tilde{\bm{b}}\\}\right)}{\sqrt{\sum_{i=1}^{N}\left(\tilde{\bm{r}}_{i}-\mathbb{E}\\{\tilde{\bm{r}}\\}\right)^{2}}\sqrt{\left(\tilde{\bm{b}}_{i}-\mathbb{E}\\{\tilde{\bm{b}}\\}\right)^{2}}},$
(12)

where $\mathbb{E}\\{\cdot\\}$ is the expectation operator, $N$ is the number
of already predicted pixels, and $\tilde{\bm{r}}$ and $\tilde{\bm{b}}$ are the
already predicted pixels of the current block of the reference signal and the
channel, respectively. Subsequently, this reference is used to build a linear
regression model very similar to the one used in cross-spectral
reconstruction. Therefore, the linear model reads as

$\tilde{\bm{b}}=\alpha\cdot\tilde{\bm{r}}+\beta$ (13)

with parameters $\alpha$ and $\beta$ to be estimated. Again, the parameters
are found by minimizing the $l_{2}$-norm

$\hat{\alpha},\hat{\beta}=\operatorname*{\arg\\!\min}_{\alpha,\beta}||\alpha\cdot\tilde{\bm{r}}+\beta-\tilde{\bm{b}}||_{2}^{2}.$
(14)

Finally, this model is employed for the prediction of the pixel

$\hat{b}=\hat{\alpha}\cdot r+\hat{\beta}.$ (15)

The coding order of the remaining channels is found by sorting the different
channels according to their structural similarity index with respect to one
anchor channel. Since there are always three reference channels, the least
similar reference channel is then replaced by the channel, which was just
coded. Note that this method was focused on multispectral image coding rather
than hyperspectral image coding.

####  3.2.2 Novel Approach

Figure 11: The basic concept of this novel inter hyperspectral video coder.
First, motion compensation is executed on the frame $t$ exploiting the
neighboring frames. Afterwards, the current channel is spectrally compensated
using the last three channels.

For hyperspectral image coding, there is implicitly more structure given,
since neighboring channels are also neighbors in the spectrum. Therefore,
since spectra are typically smooth, the correlation between two neighboring
channels is the highest. The temporal and spectral compensation of the novel
hyperspectral video coder is shown in Fig. 11. Assuming the channels are
sorted in ascending order according to their wavelength, the proposed
hyperspectral video coder compresses the first three channels as an RGB video
with a standard video coder using only intra-prediction to ensure a high
quality for motion estimation using these channels. Since this work is based
on PRIBP, the base video coder is again HEVC. Thus, the first three decoded
channels for every time step $t$ read as

$\hat{C}^{t}_{0},\hat{C}^{t}_{1},\hat{C}^{t}_{2}=\text{D}^{\text{HEVC}}\left(\text{C}^{\text{HEVC}}\left(\left[C^{t}_{0},C^{t}_{1},C^{t}_{2}\right]\right)\right),$
(16)

where $\text{C}^{\text{HEVC}}\left(\cdot\right)$ and
$\text{D}^{\text{HEVC}}\left(\cdot\right)$ are the coder and decoder,
respectively, and $C^{t}_{i}$ is the $i$-th grayscale channel of a
hyperspectral image of size $M\times N$ at time step $t$. To avoid issues with
error propagations, every second frame of every other channel is coded by the
hyperspectral intra encoder, while for all other frames a motion-compensated
residual is compressed by the hyperspectral intra coder. Otherwise, the errors
during coding would propagate into the next frame leading to even more errors
and a higher residual to code, which would then turn into an even lower
quality or higher rate. This residual is calculated by exploiting temporal
correlation. Therefore, the decoded first three channels are used to estimate
motion, which is done by PWC-Net[35]. This is done using forward and backward
motion estimation

$\begin{split}\bm{V}^{t}_{\text{fw}}&=\text{PWC}\left(\left[\hat{C}^{t-1}_{0},\hat{C}^{t-1}_{1},\hat{C}^{t-1}_{2}\right],\left[\hat{C}^{t}_{0},\hat{C}^{t}_{1},\hat{C}^{1}_{2}\right]\right)\\\
\bm{V}^{t}_{\text{bw}}&=\text{PWC}\left(\left[\hat{C}^{t+1}_{0},\hat{C}^{t+1}_{1},\hat{C}^{t+1}_{2}\right],\left[\hat{C}^{t}_{0},\hat{C}^{t}_{1},\hat{C}^{1}_{2}\right]\right),\end{split}$
(17)

where $\bm{V}^{t}$ is an image of motion vectors for time step $t$. Then,
these motion estimations can be used for motion compensation

$\begin{split}P^{t}_{\text{fw},i}[m,n]&=\tilde{C}^{t-1}_{i}\left(m+\bm{V}^{t}_{\text{fw},1}[m,n],n+\bm{V}^{t}_{\text{fw},2}[m,n]\right)\\\
P^{t}_{\text{bw},i}[m,n]&=\tilde{C}^{t+1}_{i}\left(m+\bm{V}^{t}_{\text{bw},1}[m,n],n+\bm{V}^{t}_{\text{bw},2}[m,n]\right),\end{split}$
(18)

where $\tilde{C}\left(m,n\right)$ is the interpolated and extended image of
$\hat{C}\left[m,n\right]$. Moreover, corresponding masks, when the motion
compensation tries to pull a pixel outside the image scope, need to be
calculated

$M^{t}_{j}[m,n]=\left\\{\begin{array}[]{ c l }1,&\text{if }\begin{aligned}
&0\leq m+\bm{V}^{t}_{j,1}[m,n]\leq M-1\ \land\\\ &0\leq
n+\bm{V}^{t}_{j,2}[m,n]\leq N-1\end{aligned}\\\
0,&\text{else},\end{array}\right.$ (19)

where $j$ is in the set of $\left\\{\text{fw},\text{bw}\right\\}$. Since two
predictions for the same image are made, there has to be a merge procedure for
these. In general, the simple average of these two predictions is taken to
merge them. However, at the border of the image, it occurs that the warping
tries to warp pixels outside of the original image plane to the corresponding
pixel in the predicted image. In these situations, the merge procedure
entirely relies on the other prediction as long as this prediction is valid.
If no prediction is valid, again the average of both predictions is taken as
best guess. This summarizes to the final prediction

$\begin{split}P^{t}_{i}=\left\\{\begin{array}[]{c
l}\frac{P^{t}_{\text{fw},i}+P^{t}_{\text{bw},i}}{2},&\begin{aligned} &\text{if
}\left(M^{t}_{\text{fw}}[m,n]=1\land M^{t}_{\text{bw}}[m,n]=1\right)\\\
&\lor\left(M^{t}_{\text{fw}}[m,n]=0\land
M^{t}_{\text{bw}}[m,n]=0\right)\end{aligned}\\\ P^{t}_{\text{fw},i},&\text{if
}M^{t}_{\text{fw}}[m,n]=1\land M^{t}_{\text{bw}}[m,n]=0\\\
P^{t}_{\text{bw},i},&\text{if }M^{t}_{\text{fw}}[m,n]=0\land
M^{t}_{\text{bw}}[m,n]=1.\end{array}\right.\end{split}$ (20)

To calculate the residual for all odd frames, the predicted image is
subtracted from the original decoded image

$R^{t}_{i}=C^{t}_{i}-P^{t}_{i}\quad\forall\ t\bmod 2=1.$ (21)

This yields three residual images which serve as reference residual images for
the hyperspectral intra coder, later. Moreover, as in PRIBP, for all-intra
coded frames the first three channels serve as reference channel for the first
hyperspectral intra coded image.

Subsequently, all other channels are coded iteratively. Even frames are coded
using all-intra PRIBP

$\hat{C}^{t}_{i}=\text{D}^{\text{PRIBP}}\left(\text{C}^{\text{PRIBP}}\left(C^{t}_{i},\left[\hat{C}^{t}_{i-1},\hat{C}^{t}_{i-2},\hat{C}^{t}_{i-3}\right]\right)\right)\\\
\forall\ i\geq 3\land t\bmod 2=0,$ (22)

where
$\text{C}^{\text{PRIBP}}\left(C^{t}_{i},\left[\hat{C}^{t}_{i-1},\hat{C}^{t}_{i-2},\hat{C}^{t}_{i-3}\right]\right)$
means that the $i$-th channel is coded using decoded channels $i-1$, $i-2$ and
$i-3$ as reference channels. After being decoded, they replace the oldest
reference in the reference buffer for that frame and serve as reference
channel for the next three channels. This is done, since, as already
described, spectra are typically smooth and therefore, the last three channels
have the highest correlation with the current channel from a physical point of
view. Then, the decoded intra frames at timesteps $t-1$ and $t+1$ are used for
predicting the current frame at time $t$ using the motion vectors estimated
from the first three channels in (17). The motion-compensated prediction is
calculated using (18), (19) and (20). Afterwards, the residual with the frame
to code is calculated by (21). This residual is fed into PRIBP, which uses
three reference images again

$\hat{R}^{t}_{i}=\text{D}^{\text{PRIBP}}\left(\text{C}^{\text{PRIBP}}\left(R^{t}_{i},\left[\hat{R}^{t}_{i-1},\hat{R}^{t}_{i-2},\hat{R}^{t}_{i-3}\right]\right)\right)\\\
\forall\ i\geq 3\land t\bmod 2=1.$ (23)

In this case, the reference images are residual signals as well. Again, as
soon as this residual is decoded again, it serves as reference residual for
the next three channels.

####  3.2.3 Evaluation

Our novel hyperspectral video database can directly be used to evaluate the
gains of the proposed hyperspectral video coder. For the evaluation, only the
videos of the center camera are used. The Quantization Parameters (QPs) for
evaluation are set to 22, 27, 32 and 37. These QPs are extracted from the
common test conditions for HEVC and thus common practice[36].

In general, different coders are evaluated by setting up rate-distortion
curves, i.e., a plot that shows the data rate in terms of bitrate on the
x-axis and distortion, here in terms of PSNR, on the y-axis. The different
aforementioned QPs are datapoints for this curve. Between these data points
the curve is interpolated. Thus, if the curve of a second coder is to the top
left compared to the curve of the base coder, it achieves a better PSNR for
the same bitrate, or vice versa, a lower bit rate for the same PSNR. This
average distance between rate-distortion curves can be measured using the
Bjontegaard-Delta (BD) for the rate direction and PSNR direction[37]. The
results for each scene are summarized in Table 3, which shows the BD rate and
PSNR. The proposed hyperspectral video coder is first compared to the all-
intra coder PRIBP. For six out of seven scenes, the proposed coder outperforms
the intra procedure. PWC-Net has problems estimating the motion properly for
the scene forest, since there are a lot of leaves, which leads to a high
frequency change between background and leaves. Thus, the motion estimation
would have to output a motion field with a very high spatial resolution to
enhance the coder. In this case, PWC-Net does not achieve this and even has a
bad influence on the coding result. On average, the novel hyperspectral coder
saves roughly 3% in comparison to PRIBP. Note that this average is calculated
over all sizes and PSNRs and not over the rate and PSNR values from the table.
Thus, videos with a large file size influences this average more.

Moreover, the proposed coding technique is also compared to the random access
mode of High Efficiency Video Coding (HEVC), thus an RGB inter coder. For
that, every three consecutive channels were merged into an RGB video. The
remaining channels are coded as grayscale videos. For 31 hyperspectral
channels, this results in 10 RGB videos and 1 grayscale video to code for the
standard HEVC. The results vary much more in comparison, however, the average
bit rate saving of the proposed hyperspectral video coder is over 50% in
comparison the inter mode of HEVC. It is noticeable that HEVC has problems
with high frequency content in comparison, since the scenes family house,
outdoor and forest are compressed much worse. On the other side, for scenes
with mostly homogeneous regions and movements, HEVC even outperforms the
proposed hyperspectral video coder.

Table 3: Evaluation of our novel inter hyperspectral coding against all-intra hyperspectral coding PRIBP and HEVC inter coding (random access). The values are given in Bjontegaard-Delta rate and PSNR, respectively. Lower is better for rate, while higher is better for PSNR. | PRIBP | HEVC Inter Coding
---|---|---
| BD Rate | BD PSNR | BD Rate | BD PSNR
Fam. house | -10.41% | 0.59 dB | -54.55% | 4.02 dB
M. seaport | -7.15% | 0.35 dB | -22.82% | 1.16 dB
City | -9.29% | 0.44 dB | -19.36% | 0.92 dB
Outdoor | -9.22% | 0.47 dB | -65.81% | 4.82 dB
Forest | 13.24% | -0.59 dB | -73.49% | 6.09 dB
Indoor | -6.77% | 0.22 dB | 12.71% | -0.38 dB
Lab | -9.21% | 0.57 dB | 42.61% | -1.85 dB
Average | -3.36% | 0.15 dB | -52.46% | 3.27 dB

###  3.3 Other Applications

Other applications include spectral reconstruction[32]. Spectral
reconstruction has the goal of reconstructing light spectra from multispectral
images or even RGB images. The database can also be used for cross-spectral
depth estimation[38], where the disparity between different views, which
record different ranges in the spectrum, is estimated, e.g., using the
multispectral camera array shown in Fig. 6. Furthermore, one can also evaluate
different camera setups like trinocular arrangements[39]. Moreover, depth
estimation algorithms and image enhancement procedures that rely on grayscale
and RGB cameras[40] can be used. In a similar direction steer studies that
investigate proper color to grayscale conversions for these types of
algorithms[41]. A related topic to stereo matching is scene reconstruction
using different cameras[42].

In general, different fusion techniques related to depth sensors can be
evaluated. A very obvious example is depth estimation aided by a depth map
sensor, e.g., a LIDAR sensor[43]. Furthermore, a depth sensor can be used to
help upsampling[44] as well. An even better fit to this database is the fusion
of RGB images, hyperspectral recordings and LIDAR data for urban land cover
classification[45].

Of course reconstruction algorithms for different types of hyperspectral
cameras can be easily evaluated like demosaicing for multispectral filter
arrays[46]. Combinations of these techniques also can be evaluated, e.g., a
mixture of multispectral filter arrays and camera arrays.

Finally, multispectral and hyperspectral denoising algorithms can be
evaluated, e.g., techniques based on sparse matrix decomposition[47]. The
noise can be even further reduced exploiting the temporal correlation between
frames, which can be evaluated using the proposed database.

##  4 Conclusion

This paper introduced a novel synthetic hyperspectral video database
containing seven scenes. All scenes were rendered using a camera array with
nine cameras arranged in a three times three grid. This camera array moves
through the scenes for 30 frames. Since depth maps are provided for every
scene, camera and frame, this database provides a high resolution in all three
spatial as well as spectral and temporal dimension of the corresponding scene.
Therefore, this database can serve as validation database for many
applications ranging from spectral reconstruction over diverse stereo matching
problems like cross spectral disparity estimation to sensor fusion problem
like combining hyperspectral data, RGB images and data from depth sensing
devices. Finally, two applications were covered in detail by exploiting the
temporal dimension of the data. The proposed two novel algorithms outperformed
its non-temporal versions. Furthermore, it proved that this new database can
be used to evaluate algorithms in diverse image processing tasks.

Funding The authors gratefully acknowledge that this work has been supported
by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under
project number 491814627. Disclosures The authors declare no conflicts of
interest. Data availability Data underlying the results presented in this
paper are available in Dataset 1, Ref. [48].

## References

  * [1] N. A. Hagen and M. W. Kudenov, “Review of snapshot spectral imaging technologies,” Optical Engineering 52, 1 – 23 (2013).
  * [2] S. P. Garaba and H. M. Dierssen, “An airborne remote sensing case study of synthetic hydrocarbon detection using short wave infrared absorption features identified from marine-harvested macro- and microplastics,” Remote Sensing of Environment 205, 224–235 (2018).
  * [3] P. J. Williams and S. Kucheryavskiy, “Classification of maize kernels using NIR hyperspectral imaging,” Food Chemistry 209, 131–138 (2016).
  * [4] G. Edelman, E. Gaston, T. van Leeuwen, P. Cullen, and M. Aalders, “Hyperspectral imaging for non-contact analysis of forensic traces,” Forensic Science International 223, 28–39 (2012).
  * [5] F. Xiong, J. Zhou, and Y. Qian, “Material Based Object Tracking in Hyperspectral Videos,” IEEE Transactions on Image Processing 29, 3719–3733 (2020).
  * [6] Z. Han, A. Zhang, X. Wang, Z. Sun, M. D. Wang, and T. Xie, “In vivo use of hyperspectral imaging to develop a noncontact endoscopic diagnosis support system for malignant colorectal tumors,” Journal of Biomedical Optics 21, 016001 (2016).
  * [7] F. Yasuma, T. Mitsunaga, D. Iso, and S. K. Nayar, “Generalized Assorted Pixel Camera: Postcapture Control of Resolution, Dynamic Range, and Spectrum,” IEEE Transactions on Image Processing 19, 2241–2253 (2010).
  * [8] J. Eckhard, T. Eckhard, E. M. Valero, J. L. Nieves, and E. G. Contreras, “Outdoor scene reflectance measurements using a Bragg-grating-based hyperspectral imager,” Appl. Opt. 54, D15–D24 (2015).
  * [9] S. Hordley, G. Finalyson, and P. Morovic, “A multi-spectral image database and its application to image rendering across illumination,” in _Third International Conference on Image and Graphics (ICIG’04),_ (2004), pp. 394–397.
  * [10] S. L. Moan, S. T. George, M. Pedersen, J. Blahová, and J. Y. Hardeberg, “A database for spectral image quality,” in _Image Quality and System Performance XII,_ vol. 9396 M.-C. Larabi and S. Triantaphillidou, eds., International Society for Optics and Photonics (SPIE, 2015), pp. 225 – 232\.
  * [11] D. H. Foster, K. Amano, S. M. C. Nascimento, and M. J. Foster, “Frequency of metamerism in natural scenes,” J. Opt. Soc. Am. A 23, 2359–2372 (2006).
  * [12] T. Zickler and A. Chakrabarti, “Statistics of real-world hyperspectral images,” in _2013 IEEE Conference on Computer Vision and Pattern Recognition,_ (IEEE Computer Society, Los Alamitos, CA, USA, 2011), pp. 193–200.
  * [13] B. Arad and O. Ben-Shahar, “Sparse Recovery of Hyperspectral Signal from Natural RGB Images,” in _Computer Vision – ECCV 2016,_ B. Leibe, J. Matas, N. Sebe, and M. Welling, eds. (Springer International Publishing, Cham, 2016), pp. 19–34.
  * [14] A. Mian and R. Hartley, “Hyperspectral video restoration using optical flow and sparse coding,” Opt. Express 20, 10658–10673 (2012).
  * [15] L. Gómez-Chova, L. Alonso, L. Guanter, G. Camps-Valls, J. Calpe, and J. Moreno, “Correction of systematic spatial noise in push-broom hyperspectral sensors: application to CHRIS/PROBA images,” Appl. Opt. 47, F46–F60 (2008).
  * [16] F. Koenig and W. Praefcke, “Practice of multispectral image acquisition,” in _Electronic Imaging: Processing, Printing, and Publishing in Color,_ vol. 3409 J. Bares, ed., International Society for Optics and Photonics (SPIE, 1998), pp. 34 – 41.
  * [17] N. Gat, G. Scriven, J. Garman, M. D. Li, and J. Zhang, “Development of four-dimensional imaging spectrometers (4D-IS),” in _Imaging Spectrometry XI,_ vol. 6302 S. S. Shen and P. E. Lewis, eds., International Society for Optics and Photonics (SPIE, 2006), pp. 179 – 189.
  * [18] M. R. Descour, C. E. Volin, E. L. Dereniak, T. M. Gleeson, M. F. Hopkins, D. W. Wilson, and P. D. Maker, “Demonstration of a computed-tomography imaging spectrometer using acomputer-generated hologram disperser,” Appl. Opt. 36, 3694–3698 (1997).
  * [19] J. D. Matchett, R. I. Billmers, E. J. Billmers, and M. E. Ludwig, “Volume holographic beam splitter for hyperspectral imaging applications,” in _Novel Optical Systems Design and Optimization X,_ vol. 6668 R. J. Koshel and G. G. Gregory, eds., International Society for Optics and Photonics (SPIE, 2007), pp. 191 – 198.
  * [20] Y. Monno, S. Kikuchi, M. Tanaka, and M. Okutomi, “A Practical One-Shot Multispectral Imaging System Using a Single Image Sensor,” IEEE Transactions on Image Processing 24, 3048–3059 (2015).
  * [21] R. Shogenji, Y. Kitamura, K. Yamada, S. Miyatake, and J. Tanida, “Multispectral imaging using compact compound optics,” Opt. Express 12, 1643–1655 (2004).
  * [22] M. E. Gehm, R. John, D. J. Brady, R. M. Willett, and T. J. Schulz, “Single-shot compressive spectral imaging with a dual-disperser architecture,” Opt. Express 15, 14013–14027 (2007).
  * [23] N. Genser, J. Seiler, and A. Kaup, “Camera Array for Multi-Spectral Imaging,” IEEE Transactions on Image Processing 29, 9234–9249 (2020).
  * [24] B. O. Community, _Blender - a 3D modelling and rendering package_ , Blender Foundation, Stichting Blender Foundation, Amsterdam (2018).
  * [25] T. J. Purcell, I. Buck, W. R. Mark, and P. Hanrahan, “Ray Tracing on Programmable Graphics Hardware,” ACM Trans. Graph. 21, 703–712 (2002).
  * [26] A. Buades, B. Coll, and J.-M. Morel, “Non-Local Means Denoising,” Image Processing On Line 1, 208–212 (2011).
  * [27] J. J. McCann and A. Rizzi, “Veiling glare: the dynamic range limit of HDR images,” in _Human Vision and Electronic Imaging XII,_ vol. 6492 B. E. Rogowitz, T. N. Pappas, and S. J. Daly, eds., International Society for Optics and Photonics (SPIE, 2007), p. 649213.
  * [28] K. Rosenhauer and K. Rosenbruch, “Flare and optical transfer function,” Appl. Opt. 7, 283–287 (1968).
  * [29] L. Waller, S. S. Kou, C. J. R. Sheppard, and G. Barbastathis, “Phase from chromatic aberrations,” Opt. Express 18, 22817–22825 (2010).
  * [30] H. T. Hytti, “Characterization of digital image noise properties based on RAW data,” in _Image Quality and System Performance III,_ vol. 6059 L. C. Cui and Y. Miyake, eds., International Society for Optics and Photonics (SPIE, 2006), p. 60590A.
  * [31] F. Sippel, J. Seiler, and A. Kaup, “Spatio-spectral Image Reconstruction Using Non-local Filtering,” in _2021 International Conference on Visual Communications and Image Processing (VCIP),_ (2021), pp. 1–5.
  * [32] F. Sippel, J. Seiler, N. Genser, and A. Kaup, “Structure-preserving spectral reflectance estimation using guided filtering,” J. Opt. Soc. Am. A 37, 1695–1710 (2020).
  * [33] A. Meyer, N. Genser, and A. Kaup, “Multispectral Image Compression Based on HEVC Using Pel-Recursive Inter-Band Prediction,” in _2020 IEEE 22nd International Workshop on Multimedia Signal Processing (MMSP),_ (2020), pp. 1–6.
  * [34] G. J. Sullivan, J.-R. Ohm, W.-J. Han, and T. Wiegand, “Overview of the High Efficiency Video Coding (HEVC) Standard,” IEEE Transactions on Circuits and Systems for Video Technology 22, 1649–1668 (2012).
  * [35] D. Sun, X. Yang, M.-Y. Liu, and J. Kautz, “PWC-Net: CNNs for Optical Flow Using Pyramid, Warping, and Cost Volume,” in _2018 IEEE/CVF Conference on Computer Vision and Pattern Recognition,_ (2018), pp. 8934–8943.
  * [36] J. Boyce, K. Suehring, X. Li, and V. Seregin, “JVET-J1010: JVET common test conditions and software reference configurations,” (2018).
  * [37] G. Bjøntegaard, “Calculation of average PSNR differences between RD curves,” in _VCEG-M33, document,_ (2001).
  * [38] N. Genser, A. Spruck, J. Seiler, and A. Kaup, “Deep Learning Based Cross-Spectral Disparity Estimation For Stereo Imaging,” in _2020 IEEE International Conference on Image Processing (ICIP),_ (2020), pp. 2536–2540.
  * [39] M. Mozerov, J. Gonzàlez, X. Roca, and J. J. Villanueva, “Trinocular stereo matching with composite disparity space image,” in _2009 16th IEEE International Conference on Image Processing (ICIP),_ (2009), pp. 2089–2092.
  * [40] H.-G. Jeon, J.-Y. Lee, S. Im, H. Ha, and I. S. Kweon, “Stereo Matching with Color and Monochrome Cameras in Low-Light Conditions,” in _2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR),_ (2016), pp. 4086–4094.
  * [41] L. Benedetti, M. Corsini, P. Cignoni, M. Callieri, and R. Scopigno, “Color to gray conversions in the context of stereo matching algorithms: An analysis and comparison of current methods and an ad-hoc theoretically-motivated technique for image matching,” Machine Vision and Applications 23, 327–348 (2012).
  * [42] Z. Murez, T. van As, J. Bartolozzi, A. Sinha, V. Badrinarayanan, and A. Rabinovich, “Atlas: End-to-End 3D Scene Reconstruction from Posed Images,” in _Computer Vision – ECCV 2020,_ A. Vedaldi, H. Bischof, T. Brox, and J.-M. Frahm, eds. (Springer International Publishing, Cham, 2020), pp. 414–431.
  * [43] H. Gao, B. Cheng, J. Wang, K. Li, J. Zhao, and D. Li, “Object Classification Using CNN-Based Fusion of Vision and LIDAR in Autonomous Vehicle Environment,” IEEE Transactions on Industrial Informatics 14, 4224–4231 (2018).
  * [44] J. Park, H. Kim, Y.-W. Tai, M. S. Brown, and I. S. Kweon, “High-Quality Depth Map Upsampling and Completion for RGB-D Cameras,” IEEE Transactions on Image Processing 23, 5559–5572 (2014).
  * [45] R. Hänsch and O. Hellwich, “Fusion of Multispectral LiDAR, Hyperspectral, and RGB Data for Urban Land Cover Classification,” IEEE Geoscience and Remote Sensing Letters 18, 366–370 (2021).
  * [46] K. Feng, Y. Zhao, J. C.-W. Chan, S. G. Kong, X. Zhang, and B. Wang, “Mosaic Convolution-Attention Network for Demosaicing Multispectral Filter Array Images,” IEEE Transactions on Computational Imaging 7, 864–878 (2021).
  * [47] T. Xie, S. Li, and B. Sun, “Hyperspectral Images Denoising via Nonconvex Regularized Low-Rank and Sparse Matrix Decomposition,” IEEE Transactions on Image Processing 29, 44–56 (2020).
  * [48] “HyViD,” https://github.com/FAU-LMS/HyViD.


# On the error of best polynomial approximation of composite functions
††thanks: This research has been accomplished within “Research ITalian network
on Approximation” (RITA).

Luisa Fermo Department of Mathematics and Computer Science, University of
Cagliari, Via Ospedale 72, 09124 Cagliari, Italy, email<EMAIL_ADDRESS>Corresponding author Concetta Laurita Department of Mathematics, Computer
Science and Economics, University of Basilicata, Viale dell’Ateneo Lucano 10,
85100 Potenza, Italy, email<EMAIL_ADDRESS><EMAIL_ADDRESS>Maria Grazia Russo 33footnotemark: 3

###### Abstract

The purpose of the paper is to provide a characterization of the error of the
best polynomial approximation of composite functions in weighted spaces. Such
a characterization is essential for the convergence analysis of numerical
methods applied to non-linear problems or for numerical approaches that make
use of regularization techniques to cure low smoothness of the solution. This
result is obtained through an estimate of the derivatives of composite
functions in weighted uniform norm.

###### keywords:

Derivatives of composite functions, weighted spaces, Sobolev-type spaces, best
polynomial approximation

###### AMS:

41A10, 41A25, 26D10

##  1 Introduction

Consider a fixed Jacobi weight

$u(x)=v^{\gamma,\delta}(x):=(1-x)^{\gamma}(1+x)^{\delta},\qquad
x\in(-1,1),\quad\gamma,\delta\geq 0.$

We denote by $C_{u}$ the space of all functions $f$ continuous on $(-1,1)$ and
satisfying

$\lim_{x\rightarrow+1}f(x)u(x)=0\quad\mbox{if
$\gamma>0$},\qquad\mbox{and}\qquad\lim_{x\rightarrow-1}f(x)u(x)=0\quad\mbox{if
$\delta>0$},$

equipped with the norm

$\|f\|_{C_{u}}:=\|fu\|_{\infty}=\max_{x\in[-1,1]}|f(x)|u(x).$

It is well-known (see for instance [6]) that the Weierstrass approximation
theorem holds true in this Banach space and we have

$f\in C_{u}\ \Longleftrightarrow\ \lim_{m\rightarrow\infty}E_{m}(f)_{u}=0,$

where $E_{m}(f)_{u}$ denotes the error of best approximation of $f\in C_{u}$
in the space $\mathbb{P}_{m}$ of all algebraic polynomials of degree at most
$m$, namely

$E_{m}(f)_{u}:=\inf_{P\in\mathbb{P}_{m}}\|(f-P)\|_{C_{u}}.$

The rate of convergence of such error as $m\rightarrow\infty$ depends on the
smoothness of the function $f$. A renewed inequality used in estimating the
best approximation is the following Favard type inequality (see for instance
[6] and the references therein)

$E_{m}(f)_{u}\leq\frac{\mathcal{C}}{m^{r}}\|f^{(r)}\varphi^{r}u\|_{\infty}$
(1)

where here and in the following $\mathcal{C}$ is a positive absolute constant,
$f^{(r)}$ denotes the $r$th derivative of $f$ and $\varphi(x)=\sqrt{1-x^{2}}$.

In different contexts in which the “global approximation” (i.e. in processes
involving the polynomial approximation) is used, it becomes crucial to
estimate the best polynomial approximation of composed functions. This
situation obviously occurs in non linear problems, but also in other contexts
like for instance in methods using regularization techniques (see, for
instance, [2], [4], [5], [8], [9], and [10]).

Looking at the Favard inequality, in order to estimate $E_{m}(f\circ g)_{u}$
we need an upper bound for the weighted norm of the derivatives of the
composite functions. A crucial step for our aim is to provide, in weighted
uniform norm, an estimate of the derivatives of a given function in terms of a
fixed higher derivative of the same function. According to our knowledge, this
is only known for functions belonging to $C^{r}([a,b])$. In fact, in [3, Lemma
2.1] for a given function $f\in C^{r}([a,b])$ the following inequality was
stated

$\|f^{(k)}\|_{\infty}\leq\mathcal{C}\left(\frac{\|f\|_{\infty}}{(b-a)^{k}}+(b-a)^{r-k}\|f^{(r)}\|_{\infty}\right),\quad
0<k<r,$ (2)

where $\mathcal{C}$ is a positive constant independent of the interval $[a,b]$
and $f$.

In this paper, we first prove the analogous of (2) in the weighted case and
then we give the estimate of the derivative of a composite function in a very
general case, i.e. the case in which $f$ is a multivariate function and $g$ is
a vector of functions. We will present the results in $[-1,1]$, without loss
of generality, since by linearity analogous results can be deduced in the
generic interval $[a,b]$ of $\mathbb{R}$.

##  2 Sobolev-type spaces

Let us introduce the weighted Sobolev–type space of order $1\leq
r\in\mathbb{N}$ [6]

${W}^{r}_{u}=\left\\{f\in C_{u}:f^{(r-1)}\in
AC((-1,1)),\,\|f^{(r)}\varphi^{r}u\|_{\infty}<\infty\right\\},$

where $\varphi(x)=\sqrt{1-x^{2}}$ and $AC((-1,1))$ denotes the set of all
absolutely continuous functions on $(-1,1)$. We equip $W^{r}_{u}$ with the
norm

$\|f\|_{W^{r}_{u}}:=\|fu\|_{\infty}+\|f^{(r)}\varphi^{r}u\|_{\infty}.$

For multivariate functions $f:\Omega\rightarrow\mathbb{R}$ with $\Omega$ open
subset of $\mathbb{R}^{n}$, we define the Sobolev space ${\bm{W}}^{r}(\Omega)$
[1] as the set of all functions $f$ in $\Omega$ such that for every $n$-tuples
of nonnegative integers $\ell=(\ell_{1},\ell_{2},\dots,\ell_{n})$, with
$|\ell|=\sum_{i=1}^{n}\ell_{i}\leq r$, the mixed partial derivatives

$D^{\ell}f=\frac{\partial^{\ell_{1}+\ell_{2}+\dots+\ell_{n}}f}{\partial
x_{1}^{\ell_{1}}\partial x_{2}^{\ell_{2}}\dots\partial x_{n}^{\ell_{n}}}$

exist and $\|D^{\ell}f\|_{\infty}<\infty$ . We endow this space with the norm

$\|f\|_{{\bm{W}}^{r}(\Omega)}=\|f\|_{\infty}+\sum_{1\leq|\ell|\leq
r}\|D^{\ell}f\|_{\infty}.$ (3)

The spaces $(W^{r}_{u},\|\cdot\|_{W^{r}_{u}})$ and
$({\bm{W}}^{r}(\Omega),\|\cdot\|_{{\bm{W}}^{r}(\Omega)})$ are Banach spaces.

##  3 The Faa di Bruno’s formula

Let $f:\Omega\subseteq\mathbb{R}^{n}\rightarrow\mathbb{R}$ and
$g:A\subseteq\mathbb{R}\rightarrow\mathbb{R}^{n}$ with
$g(x)=\left(g_{1}(x),g_{2}(x),\ldots,g_{n}(x)\right)$, be functions such that
the range of $g$ is contained in the domain of $f$ and for which all the
necessary derivatives are defined.
In the case $n=1$ the well-known Fáa di Bruno formula states

$(f\circ g)^{(r)}(x)=\sum\frac{r!}{k_{1}!k_{2}!\ldots
k_{r}!}f^{(k1+k2+\ldots+k_{r})}(g(x))\prod_{i=1}^{r}\left(\frac{g^{(i)}(x)}{i!}\right)^{k_{i}}$
(4)

where the sum is over all $r$-tuples of nonnegative integers
$\left(k_{1},k_{2},\ldots,k_{r}\right)$ such that
$k_{1}+2k_{2}+\ldots+rk_{r}=r.$

Combining the terms with the same value of $k_{1}+k_{2}+\ldots+k_{r}=k$ and
noticing that $k_{j}$ has to be zero for $j>r-k+1$ lead to the following
equivalent but simpler formula

$(f\circ
g)^{(r)}(x)=\sum_{k=1}^{r}f^{(k)}(g(x))\,{\bf{B}}_{r,k}(g^{\prime}(x),g^{\prime\prime}(x),\dots
g^{(r-k+1)}(x))$

where ${\mathbf{B}}_{r,k}(x_{1},x_{2},\dots x_{r-k+1})$ are the partial or
incomplete exponential Bell polynomials defined as

${\bf{B}}_{r,k}(x_{1},x_{2},\dots x_{r-k+1})=\sum\frac{r!}{k_{1}!k_{2}!\dots
k_{r-k+1}!}\prod_{i=1}^{r-k+1}\left(\frac{x_{i}}{i!}\right)^{k_{i}},$

with the sum taken over all sequences
$\left(k_{1},k_{2},\ldots,k_{r-k+1}\right)$ of nonnegative integers such that
these two conditions are satisfied

$\sum_{i=1}^{r-k+1}k_{i}=k,\qquad\sum_{i=1}^{r-k+1}ik_{i}=r.$

The formula (4) for the $r$-th derivative of a composite function can be
generalized in the case $n>1$ as (see [7])

$\displaystyle(f\circ g)^{(r)}(x)$
$\displaystyle=\displaystyle\sum_{0}\sum_{1}\sum_{2}\cdots\sum_{r}\displaystyle\frac{r!}{\prod_{i=1}^{r}\prod_{j=1}^{n}q_{ij}!\prod_{i=1}^{r}(i!)^{k_{i}}}\frac{\partial^{k}f(g(x))}{\partial
x_{1}^{p_{1}}\partial x_{2}^{p_{2}}\cdots\partial x_{n}^{p_{n}}}$
$\displaystyle\times\displaystyle\prod_{i=1}^{r}\left(g_{1}^{(i)}(x)\right)^{q_{i1}}\left(g_{2}^{(i)}(x)\right)^{q_{i2}}\cdots\left(g_{n}^{(i)}(x)\right)^{q_{in}}$
(5)

where the respective sums are over all negative solutions of the Diophantine
equations, as follows

$\displaystyle\sum_{0}\rightarrow k_{1}+2k_{2}+\ldots+rk_{r}=r,$
$\displaystyle\sum_{i}\rightarrow q_{i1}+q_{i2}+\ldots+q_{in}=k_{i},\quad
i=1.2,\ldots,r,$

and

$\displaystyle p_{j}$ $\displaystyle=q_{1j}+q_{2j}+\ldots+q_{rj},\quad
j=1,2,\ldots,n,$ $\displaystyle k$
$\displaystyle=p_{1}+p_{2}+\ldots+p_{n}=k_{1}+k_{2}+\ldots+k_{r}.$

##  4 Main results

First, let us prove the following lemma which gives an estimate of the
derivatives of a function $f\in W^{r}_{u}$, generalizing (2).

######  Lemma 1.

Let $u(x)=(1-x)^{\gamma}(1+x)^{\delta}$ be a Jacobi weight with
$0\leq\gamma,\delta<1$ and $f\in W^{r}_{u}$. Then, for $0<k<r$

$\|f^{(k)}\varphi^{k}u\|_{\infty}\leq\mathcal{C}\left(\|fu\|_{\infty}+\|f^{(r)}\varphi^{r}u\|_{\infty}\right)$
(8)

where $\mathcal{C}=\mathcal{C}(r,k,\gamma,\delta)$ is a positive constant
independent of $f$.

###### Proof.

First we note that it is enough to prove that, for any $0<k<r$,

$\|f^{(k)}\varphi^{k}u\|_{\infty}\leq\mathcal{C}\left(\|fu\|_{\infty}+\|f^{(k+1)}\varphi^{k+1}u\|_{\infty}\right)$
(9)

with $\mathcal{C}=\mathcal{C}(k)$, since (8) con be deduced from (9) by
induction on $k$.
Fix $-1\leq x\leq 0$ and $h=1/k$. By using the Taylor formula with integral
remainder, being $f^{(k)}$ locally absolutely continuous and consequently
$f^{(k+1)}$ locally integrable, we can write

$f(x+jh)=\sum_{i=0}^{k}\frac{(jh)^{i}}{i!}f^{(i)}(x)+\frac{1}{k!}\int_{0}^{jh}(jh-t)^{k}f^{(k+1)}(x+t)dt$

from which it follows

$\displaystyle h^{k}f^{(k)}(x)$
$\displaystyle=\sum_{j=0}^{k}\binom{k}{j}(-1)^{k-j}f(x+jh)$
$\displaystyle-\frac{1}{k!}\sum_{j=1}^{k}\binom{k}{j}(-1)^{k-j}\int_{0}^{jh}(jh-t)^{k}f^{(k+1)}(x+t)dt$

and, then,

$\displaystyle h^{k}|f^{(k)}(x)|\varphi^{k}(x)u(x)$ $\displaystyle\leq$
$\displaystyle\sum_{j=0}^{k}\binom{k}{j}|f(x+jh)|\varphi^{k}(x)u(x)$ (10)
$\displaystyle+$
$\displaystyle\frac{1}{k!}\sum_{j=1}^{k}\binom{k}{j}\int_{0}^{jh}(jh-t)^{k}\left|f^{(k+1)}(x+t)\right|\varphi^{k}(x)u(x)dt.$

Let us estimate the first term at the right-hand side. Since
$\varphi^{k}(x)\leq 1$ and

$\frac{u(x)}{u(x+jh)}\leq\frac{(1-x)^{\gamma}}{(1-x-jh)^{\gamma}}\leq\frac{2^{\gamma}}{(1-jh)^{\gamma}},$

being $-1\leq x\leq 0$ and $h>0$, using $\displaystyle\binom{k}{j}\leq 2^{k}$
for any $j\in\\{0,1,\ldots,k\\}$, we obtain

$\displaystyle\sum_{j=0}^{k}\binom{k}{j}|f(x+jh)|\varphi^{k}(x)u(x)$
$\displaystyle\leq$
$\displaystyle\|fu\|_{\infty}\sum_{j=0}^{k}\binom{k}{j}\frac{u(x)}{u(x+jh)}\varphi^{k}(x)$
$\displaystyle\leq$
$\displaystyle\|fu\|_{\infty}2^{\gamma+k}\sum_{j=0}^{k}\frac{1}{(1-jh)^{\gamma}}$
$\displaystyle\leq$
$\displaystyle\|fu\|_{\infty}2^{\gamma+k}\int_{0}^{k}(1-zh)^{-\gamma}dz$
$\displaystyle=$
$\displaystyle\|fu\|_{\infty}\frac{2^{\gamma+k}}{(1-\gamma)h}.$

In order to estimate the second sum in (10), we observe that for $-1\leq x\leq
0$, $1\leq j\leq k$, $h=1/k>0$ and $0\leq t\leq jh$ the following inequalities
hold true

$\displaystyle\frac{\varphi^{k}(x)}{\varphi^{k+1}(x+t)}$
$\displaystyle\leq\frac{2^{\frac{k}{2}}}{t^{\frac{1}{2}}(jh-t)^{\frac{k}{2}+\frac{1}{2}}},$
$\displaystyle\frac{u(x)}{u(x+t)}$
$\displaystyle\leq\frac{(1-x)^{\gamma}}{(1-x-t)^{\gamma}}\leq\frac{2^{\gamma}}{(1-t)^{\gamma}}\leq\frac{2^{\gamma}}{(jh-t)^{\gamma}}.$

Consequently, we deduce

$\displaystyle\frac{1}{k!}\sum_{j=1}^{k}\binom{k}{j}\int_{0}^{jh}(jh-t)^{k}\left|f^{(k+1)}(x+t)\right|\varphi^{k}(x)u(x)dt$
(11) $\displaystyle\leq$
$\displaystyle\frac{1}{k!}\left\|f^{(k+1)}\varphi^{k+1}u\right\|_{\infty}\sum_{j=1}^{k}\binom{k}{j}\int_{0}^{jh}(jh-t)^{k}\frac{\varphi^{k}(x)}{\varphi^{k+1}(x+t)}\frac{u(x)}{u(x+t)}dt$
$\displaystyle\leq$
$\displaystyle\frac{2^{\frac{k}{2}}}{k!}+\gamma\left\|f^{(k+1)}\varphi^{k+1}u\right\|_{\infty}\sum_{j=1}^{k}\binom{k}{j}\int_{0}^{jh}(jh-t)^{\frac{k-1}{2}-\gamma}t^{-\frac{1}{2}}dt$
$\displaystyle=$
$\displaystyle\frac{1}{k!}2^{\frac{k}{2}+\gamma}\left\|f^{(k+1)}\varphi^{k+1}u\right\|_{\infty}B\left(1-\gamma,\frac{1}{2}\right)h^{\frac{k}{2}-\gamma}\sum_{j=1}^{k}\binom{k}{j}j^{\frac{k}{2}-\gamma},$

where $B$ denotes the beta function.

Now, we note that for any integer $k\geq 1$, it is

$\sum_{j=1}^{k}\binom{k}{j}j^{\frac{k}{2}-\gamma}\leq
2^{k}k^{\frac{k}{2}-\gamma},$

from which, combining (10)-(11), we get

$\displaystyle\max_{-1\leq x\leq 0}\left|f^{(k)}(x)\varphi^{k}(x)u(x)\right|$
(12) $\displaystyle\leq$
$\displaystyle\frac{1}{h^{k}}\left[\frac{2^{\gamma+k}}{(1-\gamma)h}\|fu\|_{\infty}+\frac{1}{k!}2^{\frac{k}{2}+\gamma}B\left(1-\gamma,\frac{1}{2}\right)h^{\frac{k}{2}-\gamma}2^{k}k^{\frac{k}{2}-\gamma}\left\|f^{(k+1)}\varphi^{k+1}u\right\|_{\infty}\right]$
$\displaystyle=$
$\displaystyle\frac{2^{\gamma+k}}{1-\gamma}k^{k+1}\|fu\|_{\infty}+\frac{k^{k}}{k!}2^{\frac{3}{2}k+\gamma}B\left(1-\gamma,\frac{1}{2}\right)\left\|f^{(k+1)}\varphi^{k+1}u\right\|_{\infty}.$

The case $0\leq x\leq 1$ and $h=-1/k$ can be treated with similar arguments
and leads to the following estimate

$\displaystyle\max_{0\leq x\leq 1}\left|f^{(k)}(x)\varphi^{k}(x)u(x)\right|$
$\displaystyle\leq$
$\displaystyle\frac{2^{\delta+k}}{1-\delta}k^{k+1}\|fu\|_{\infty}$ (13)
$\displaystyle+$
$\displaystyle\frac{k^{k}}{k!}2^{\frac{3}{2}k+\delta}B\left(1-\delta,\frac{1}{2}\right)\left\|f^{(k+1)}\varphi^{k+1}u\right\|_{\infty}.$

Finally, combining (12) with (13), we can deduce (9) with the constant
$\mathcal{C}=\mathcal{C}(k)$ given by

$\displaystyle\mathcal{C}(k)$
$\displaystyle=\max\left\\{\frac{2^{\gamma+k}}{1-\gamma}k^{k+1},\frac{2^{\delta+k}}{1-\delta}k^{k+1},\frac{k^{k}}{k!}2^{\frac{3}{2}k+\gamma}B\left(1-\gamma,\frac{1}{2}\right),\right.$
$\displaystyle\left.\frac{k^{k}}{k!}2^{\frac{3}{2}k+\delta}B\left(1-\delta,\frac{1}{2}\right)\right\\}.$

(The exact constant which we arrived at are not important since they are not
the best possible.) ∎

The previous lemma is crucial to prove our first main result.

######  Theorem 2.

Let $f:\Omega\rightarrow\mathbb{R}$ with $\Omega$ open subset of
$\mathbb{R}^{n}$ and $g:(-1,1)\rightarrow\mathbb{R}^{n}$ such that
$\mathrm{Im}(g)\subseteq\Omega$. Moreover, let
$u(x)=(1-x)^{\gamma}(1+x)^{\delta}$ be a Jacobi weight with
$0\leq\gamma,\delta<1$. If we assume that $f\in{\bm{W}}^{r}(\Omega)$ and
$g=\left(g_{1},g_{2},\ldots,g_{n}\right)$ with $g_{j}\in W^{r}_{u}$, for
$j=1,2,\ldots,n$, then $f\circ g\in W^{r}_{u^{r}}$ and

$\displaystyle\|(f\circ
g)^{(r)}\varphi^{r}u^{r}\|_{\infty}\leq\mathcal{C}n^{r}B_{r}\|f\|_{{\bm{W}}^{r}(\Omega)}\prod_{j=1}^{n}\|g_{j}\|_{W^{r}_{u}}^{s_{j}}$
(14)

where $B_{r}$ is the $r$-th Bell number,
$\mathcal{C}=\mathcal{C}(r,\gamma,\delta)$ is a positive constant independent
of $f$ and $g$, and

$s_{j}=\left\\{\begin{array}[]{lcr}0,&\quad\mathrm{if}&\|g_{j}\|_{W^{r}_{u}}\leq
1\vspace{0.1cm}\\\
r,&\quad\mathrm{if}&\|g_{j}\|_{W^{r}_{u}}>1\end{array}\right.,\qquad
j=1,2,\ldots,n.$ (15)

###### Proof.

Faa di Bruno’s formula (3), with (3)-(3), yields

$\displaystyle\left|(f\circ g)^{(r)}(x)\varphi^{r}(x)u^{r}(x)\right|$
$\displaystyle\leq\sum_{0}u^{r-k}(x)\sum_{1}\sum_{2}\cdots\sum_{r}\frac{r!}{\prod_{i=1}^{r}\prod_{j=1}^{n}q_{ij}!}\left|\frac{\partial^{k}f(g(x))}{\partial
x_{1}^{p_{1}}\partial x_{2}^{p_{2}}\cdots\partial
x_{n}^{p_{n}}}\right|u^{p_{1}+\ldots+p_{n}}(x)$
$\displaystyle\times\prod_{i=1}^{r}\frac{\left|g_{1}^{(i)}(x)\varphi^{i}(x)\right|^{q_{i1}}\left|g_{2}^{(i)}(x)\varphi^{i}(x)\right|^{q_{i2}}\cdots\left|g_{n}^{(i)}(x)\varphi^{i}(x)\right|^{q_{in}}}{(i!)^{k_{i}}}$
$\displaystyle\leq\mathcal{C}\sum_{0}\sum_{1}\sum_{2}\cdots\sum_{r}\displaystyle\frac{r!}{\prod_{i=1}^{r}\prod_{j=1}^{n}q_{ij}!}\left|\frac{\partial^{k}f(g(x))}{\partial
x_{1}^{p_{1}}\partial x_{2}^{p_{2}}\cdots\partial x_{n}^{p_{n}}}\right|$
$\displaystyle\times\prod_{i=1}^{r}\frac{\left|g_{1}^{(i)}(x)\varphi^{i}(x)u(x)\right|^{q_{i1}}\left|g_{2}^{(i)}(x)\varphi^{i}(x)u(x)\right|^{q_{i2}}\cdots\left|g_{n}^{(i)}(x)\varphi^{i}(x)u(x)\right|^{q_{in}}}{(i!)^{k_{i}}}.$

Then, taking into account (8), we have

$\left\|g_{j}^{(i)}\varphi^{i}u\right\|_{\infty}\leq\mathcal{C}_{i}\|g_{j}\|_{W^{r}_{u}},\qquad
i=1,2,\ldots,r,\quad j=1,2,\ldots,n,$

and, in virtue of (3), we deduce

$\displaystyle\left|(f\circ g)^{(r)}(x)\varphi^{r}(x)u^{r}(x)\right|$
$\displaystyle\leq\mathcal{C}\|f\|_{{\bm{W}}^{r}(\Omega)}\displaystyle\sum_{0}\sum_{1}\sum_{2}\cdots\sum_{r}\frac{r!}{\prod_{i=1}^{r}\prod_{j=1}^{n}q_{ij}!}$
$\displaystyle\times\prod_{i=1}^{r}\displaystyle\frac{\mathcal{C}_{i}^{q_{i1}+q_{i2}+\ldots+q_{in}}\displaystyle\prod_{j=1}^{n}\left\|g_{j}\right\|_{W^{r}_{u}}^{q_{ij}}}{(i!)^{k_{i}}}$
$\displaystyle\leq\mathcal{C}\|f\|_{{\bm{W}}^{r}(\Omega)}\displaystyle\sum_{0}\sum_{1}\sum_{2}\cdots\sum_{r}\frac{r!}{\prod_{i=1}^{r}\prod_{j=1}^{n}q_{ij}!}$
$\displaystyle\times\prod_{j=1}^{n}\left\|g_{j}\right\|_{{W^{r}_{u}}}^{p_{j}}\prod_{i=1}^{r}\frac{\mathcal{C}_{i}^{k_{i}}}{(i!)^{k_{i}}}.$

Since $p_{j}\leq k\leq r$ for $j=1,2,\ldots,n$, setting
${\bar{\mathcal{C}}}=\displaystyle\max_{i=1,2,\ldots,r}\mathcal{C}_{i}$ and
multiplying and dividing by $k_{1}!k_{2}!\cdots,k_{r}!$, we can write

$\displaystyle\left|(f\circ g)^{(r)}(x)\varphi^{r}(x)\right|$
$\displaystyle\leq\mathcal{C}{\bar{\mathcal{C}}}^{r}\|f\|_{{\bm{W}}^{r}(\Omega)}\prod_{j=1}^{n}\left\|g_{j}\right\|_{W^{r}_{u}}^{s_{j}}\displaystyle\sum_{0}\frac{r!}{k_{1}!k_{2}!\cdots,k_{r}!}\displaystyle$
$\displaystyle\times\left[\prod_{i=1}^{r}\frac{1}{(i!)^{k_{i}}}\right]\sum_{1}\sum_{2}\cdots\sum_{r}\prod_{i=1}^{r}\frac{k_{i}!}{\prod_{j=1}^{n}q_{ij}!}$

with the exponents $s_{j}$, $j=1,2,\ldots,n$, defined by (15). Now, we observe
that

$\displaystyle\sum_{i}\frac{k_{i}!}{\prod_{j=1}^{n}q_{ij}!}$
$\displaystyle=\sum_{q_{i1}+q_{i2}+\ldots+q_{in}=k_{i}}\binom{k_{i}}{q_{i1}\,q_{i2}\,\ldots\,q_{in}}$
$\displaystyle=(1+1+\ldots+1)^{k_{i}}=n^{k_{i}}\qquad i=1,2,\ldots,r,$

from which it follows that

$\displaystyle\sum_{1}\sum_{2}\cdots\sum_{r}\prod_{i=1}^{r}\frac{k_{i}!}{\prod_{j=1}^{n}q_{ij}!}$
$\displaystyle=$
$\displaystyle\sum_{1}\frac{k_{1}!}{\prod_{j=1}^{n}q_{1j}!}\sum_{2}\frac{k_{2}!}{\prod_{j=1}^{n}q_{2j}!}\ldots\sum_{r}\frac{k_{r}!}{\prod_{j=1}^{n}q_{rj}!}$
$\displaystyle=$ $\displaystyle n^{k_{1}+k_{2}+\ldots+k_{r}}=n^{k}\leq n^{r}$

and, consequently,

$\displaystyle\left|(f\circ g)^{(r)}(x)\varphi^{r}(x)u^{r}(x)\right|$
$\displaystyle\leq\mathcal{C}n^{r}\|f\|_{{\bm{W}}^{r}(\Omega)}\prod_{j=1}^{n}\left\|g_{j}\right\|_{W^{r}_{u}}^{s_{j}}\displaystyle\sum_{0}\frac{r!}{k_{1}!k_{2}!\cdots,k_{r}!}\displaystyle\prod_{i=1}^{r}\frac{1}{(i!)^{k_{i}}}$
$\displaystyle=\mathcal{C}n^{r}\|f\|_{{\bm{W}}^{r}(\Omega)}\prod_{j=1}^{n}\left\|g_{j}\right\|_{W^{r}_{u}}^{s_{j}}\displaystyle\sum_{k=1}^{r}B_{r,k}(1,1,\ldots,1)$
$\displaystyle=\mathcal{C}n^{r}B_{r}\|f\|_{{\bm{W}}^{r}(\Omega)}\prod_{j=1}^{n}\left\|g_{j}\right\|_{W^{r}_{u}}^{s_{j}},$

where $\mathcal{C}$ is a postive constant independent of $f$ and $g$, i.e. the
thesis (14). ∎

By combining (14) with (1) we can immediately deduce an estimate for
$E_{m}(f\circ g)_{u}$.

######  Theorem 3.

Let $f:\Omega\rightarrow\mathbb{R}$ with $\Omega$ open subset of
$\mathbb{R}^{n}$ and $g:(-1,1)\rightarrow\mathbb{R}^{n}$ such that
$\mathrm{Im}(g)\subseteq\Omega$. Moreover, let
$u(x)=(1-x)^{\gamma}(1+x)^{\delta}$ be a Jacobi weight with
$0\leq\gamma,\delta<1$. If we assume that $f\in{\bm{W}}^{r}(\Omega)$ and
$g=\left(g_{1},g_{2},\ldots,g_{n}\right)$ with $g_{j}\in W^{r}_{\sqrt[r]{u}}$,
for $j=1,2,\ldots,n$, then

$E_{m}(f\circ
g)_{u}\leq\frac{\mathcal{C}}{m^{r}}n^{r}B_{r}\|f\|_{{\bm{W}}^{r}(\Omega)}\prod_{j=1}^{n}\|g_{j}\|_{W^{r}_{\sqrt[r]{u}}}^{s_{j}}$

where $\mathcal{C}=\mathcal{C}(r,\gamma,\delta)$ is a positive constant
independent of $f$ and $g$, $B_{r}$ is the $r$-th Bell number, and the
exponents $s_{j}$, $j=1,2,\ldots,n$, are given in (15) (with $u$ replaced with
$\sqrt[r]{u}$).

######  Corollary 4.

In the special case $n=1$ the previous estimate is simplified as follows

$E_{m}(f\circ
g)_{u}\leq\frac{\mathcal{C}}{m^{r}}B_{r}\|f\|_{{W}^{r}(\Omega)}\|g\|_{W^{r}_{\sqrt[r]{u}}}^{s_{1}}.$

## References

  * [1] H. Brezis.  Functional Analysis, Sobolev Spaces and Partial Differential Equations.  Springer, 2011.
  * [2] T. Diogo, P.M. Lima, A. Pedas, and G Vainikko.  Smoothing transformation and spline collocation for weakly singular Volterra integro-differential equations.  Appl. Numer. Math., 114:63–76, 2017.
  * [3] Z. Ditzian.  On interpolation of ${L}_{p}([a,b])$ and weighted Sobolev spaces.  Pacific J. Math., 90:307–324, 1980.
  * [4] L. Fermo and M.G. Russo.  Numerical methods for Fredholm integral equations with singular right-hand sides.  Adv. Comput. Math., 33(3):305–330, 2010.
  * [5] E.A. Galperin, E.J. Kansa, A. Makroglou, and S.A. Nelson.  Variable transformations in the numerical solution of second kind Volterra integral equations with continuous and weakly singular kernels; extensions to Fredholm integral equations.  J. Comput. Appl. Math., 115(1-2):193–211, 2000.
  * [6] G. Mastroianni and G. V. Milovanović.  Interpolation Processes. Basic Theory and Applications.  Springer Monographs in Mathematics. Springer-Verlag, Berlin, 2008.
  * [7] R. L. Mishkov.  Generalization of the formula of Faa di Bruno for a composite function with a vector argument.  Internat. J. Math. & Math. Sci, 24(7):481–491, 2000.
  * [8] G. Monegato and L. Scuderi.  High order methods for weakly singular integral equations with nonsmooth input functions.  Math. Comput., 67(224):1493–1515, 1998.
  * [9] A Pedas and G Vainikko.  Smoothing transformation and piecewise polynomial projection methods for weakly singular Fredholm integral equations.  Commun. Pure Appl. Anal., 5(2):395–413, 2006.
  * [10] E. Vainikko and G. Vainikko.  A spline product quasi-interpolation method for weakly singular Fredholm integral equations.  SIAM J. Numer. Anal., 46(4):1799–1820, 2008.