\section{Introduction}
\vspace{-5pt}
\label{sec:intro}

\begin{figure}[t!]
  \centering
  \includegraphics[width=0.93\linewidth]{figure/first_page5.png}
   \vspace{-10pt}
   \caption{Learning continuous spherical image representation. SphereSR leverages SLIIF to predict the RGB values at given spherical coordinates for SR with arbitrary image projection.}
   \label{fig:first_page}
   \vspace{-10pt}
\end{figure}

The 360$^{\circ}$ imaging has recently gained great attention in many fields including AR/VR.
In general, raw 360$^{\circ}$ images are transformed into 2D planar representations preserving the omnidirectional information, \eg, equirectangular projection (ERP)~\cite{kim2020super} and cube map projection (CP)~\cite{kim2020super} to be compatible with the
imaging pipelines. The omnidirectional images (ODIs)\footnote{Throughout the paper, we use omnidirectional images and 360$^{\circ}$ images interchangeably.} are sometimes projected back onto a sphere or transformed with different types of projection and rendered for display in some applications.

However, the angular resolution of a 360$^{\circ}$ image tends to be lower than that of a narrow field-of-view (FOV) perspective image as it is captured using fisheye lenses with the same sensor size. Moreover, the 360$^{\circ}$ image quality can be degraded during the transformation between different image projection types.
Therefore, it is imperative to super-resolve the low-resolution (LR) 360$^{\circ}$ image by considering various projections to provide a high-level visual quality in diverse conditions.
Early studies have attempted to reconstruct high-resolution (HR) ODIs by interpolating the missing data between the LR image pixels~\cite{nagahara2000super,arican2011joint,bagnato2010plenoptic}.

Recently, deep learning (DL) has brought a significant performance boost for the 2D single image super-resolution (SISR)~\cite{zhang2018image, lai2017deep,tai2017Memnet}. These methods mostly explore super-resolving the 2D LR image using high-capacity convolutional neural networks (CNNs), via \eg, residual connections~\cite{lim2017enhanced}, and learning algorithms including generative adversarial networks (GANs)~\cite{wang2018esrgan,ledig2017photo,wang2022semi}.
However, directly using these methods for 360$^{\circ}$ images represented in 2D planar representations is less applicable as the pixel density and texture complexity vary across different positions in 2D planar representations of 360$^{\circ}$ images, as pointed in \cite{deng2021lau}.

Consequently, some attempts have been made to address the SR problems for 360$^{\circ}$ imaging~\cite{ozcinar2019super,deng2021lau,Sun2017WeightedtoSphericallyUniformQE,Zhou2018WeightedtoSphericallyUniformSO}. In particular, 360-SS~\cite{ozcinar2019super} proposes a GAN-based framework using the Pix2Pix pipeline~\cite{isola2017image}; however, it focuses only on the ERP format and does not fully consider the properties of 360$^{\circ}$ images.
LAU-Net~\cite{deng2021lau} introduces a method to identify distortions of an ODI on the latitude and upsample the ODI pixels on the segmented patches. However, it leads to considerable disconnections along the patches.
In a nutshell, existing methods for ODI SR ignore the projection process of 360$^{\circ}$ images in real applications and only take the ERP image as the LR input and produce the HR ERP output.
Indeed, the 360$^{\circ}$ image can be flexibly
converted into various projection types as, in real applications, user specifies the projection type, direction, and FOV. Thus, it is vital to address the ERP distortions problems and strive to super-resolve an ODI image to an HR image with \emph{an arbitrary projection type} rather than a fixed type.

In this paper, as shown in Fig.~\ref{fig:first_page}, we propose a novel framework, called SphereSR, aiming at super-resolving an LR 360$^{\circ}$ image to an HR image \emph{with an arbitrary projection type via continuous spherical image representation}.
First,
we propose a feature extraction module that represents the spherical data based on icosahedron and efficiently extracts features on the spherical surface composed of uniform faces (Sec.~\ref{ds_ks}). As such, we solve the ERP image distortion problem and the pixel density difference according to the latitude.
Second, we propose a spherical local implicit image function (SLIIF) that can predict RGB values at arbitrary coordinates of a sphere feature map, inspired by LIIF~\cite{chen2021learning} (Sec.~\ref{image_func}).
SLIIF works on \emph{triangular face}, buttressed by the position embedding based on the normal plane polar coordinates to obtain relative coordinates on a sphere. Therefore, our method tackles pixel-misalignment issue when being projected to another ODI projection.
As a result, SphereSR can predict RGB values for any SR scale parameters. Additionally, to train SphereSR, we introduce a feature loss that measures the similarity between two projection types, leading to considerable performance enhancement (Sec.~\ref{loss_func}).
Extensive experiments on various benchmark datasets show that our method significantly surpasses existing methods.

In summary, the contributions of our paper are four-fold. (I) We propose a novel framework, called SphereSR, aiming at super-resolving an LR 360$^{\circ}$ image to an HR image with an arbitrary projection type.
(II) We propose a feature extraction module that represents the spherical data based on icosahedron and extracts features on spherical surface. (III) We propose SLIIF that predicts RGB values from the spherical coordinates. (IV) Our method achieves the significantly better performance in the extensive experiments.

\section{Related Works}
\noindent \textbf{Omnidirectional Image SR and Enhancement.}
Early ODI SR methods~\cite{Nagahara2000SuperresolutionFA, Kawasaki2006SuperresolutionOC, Arican2009L1RS, Bagnato2010PlenopticBS, Arican2011JointRA} focused on assembling and optimizing multiple LR ODIs
on spherical or hyperbolic surface.
On the other hand, as the distortion in ODI is caused due to the projection of original spherical image to the 2D planar image plane, recent research focused on tackling and solving distortion in ODI using 2D convolution to achieve qualitative result in observation space, \ie spherical surface.
Su \etal~\cite{Sun2017WeightedtoSphericallyUniformQE} and Zhou \etal ~\cite{Zhou2018WeightedtoSphericallyUniformSO} proposed evaluation methods for ODI weighted with the projected area on the spherical surface.
\cite{Sevom2018360PS} and \cite{Nishiyama2021360SI} adapted existing SISR models to ERP SR by fine-tuning or by adding a distortion map as an input to tackle different distortions.
\cite{ozcinar2019super} leveraged GAN to super-resolve an ODI by applying WS-SSIM~\cite{Zhou2018WeightedtoSphericallyUniformSO}. \cite{Zhang2020TowardRP} also proposed the GAN-based framework employing multi-frequency structures to enhance the panoramic image quality up to the high-end camera quality. \cite{Liu2020ASF} focused on the 360$^{\circ}$ image SR utilizing single-frame and multi-frame joint learning and a loss function weighted differently along the latitude. \cite{deng2021lau} considered varying pixel density and texture complexity along latitude by proposing network allowing distinct up-scaling factors along the latitude bands.
\emph{Unlike the aforementioned methods, we propose to predict the RGB values at the given spherical coordinates of an HR image with respect to an arbitrary project type from an LR 360$^{\circ}$ image.}

\noindent \textbf{2D SISR with an Arbitrary Scale.}
Research on SISR with an arbitrary scale has been actively conducted. \cite{lim2017enhanced} first proposed a method that enables multiple scale factors over one network. MetaSR~\cite{hu2019meta} achieves SR with the non-integer scale factors. However, both methods are limited to SR with the symmetric scales. Later on, \cite{wang2020learning} proposed a framework that enables asymmetric
scale factors along horizontal and vertical axes. Moreover, SRWarp~\cite{son2021srwarp} generalizes SR toward an arbitrary image transformation.
Although these methods are effective for 2D SISR with an arbitrary scale factor, they fail to be directly applied to 360$^{\circ}$ image SR
due to the difference between xy-coordinate (2D) and spherical coordinates in ODI domains. \emph{We overcome the challenge by proposing SphereSR that leverages SLIIF to predict RGB values for the arbitrary spherical coordinates}.

\begin{figure*}[t!]
  \centering
  \includegraphics[width=.85\linewidth]{figure/framework_final7.png}
   \vspace{-5pt}
   \caption{Overall framework of the proposed SphereSR.}
   \label{fig:overall-framework}
   \vspace{-10pt}
\end{figure*}

\noindent\textbf{Continuous Image Representation.}
Research on implicit neural representation (INR) has been conducted to express 3D spaces, \eg, 3D shape construction and novel view synthesis as continuous ways~\cite{mildenhall2020nerf,sitzmann2019scene,mescheder2019occupancy}. Since then, continuous image representation has been explored on the (x,y) coordinate. \cite{anokhin2021image,skorokhodov2021adversarial} used the networks to predict the RGB value of each pixel from the latent vector on (x,y) coordinate without a spatial convolution for 2D image generation. LIIF~\cite{chen2021learning} proposes to bridge between discrete and continuous representation for images on the (x,y) coordinate.
\emph{We propose SLIIF that makes continuous image representation on the unit sphere possible.}

\noindent \textbf{CNNs for Spherical Images.}
\cite{cohen2018spherical} proposed a CNN-based method on the sphere with structural characteristics called rotational equivariance. However, it requires Fourier transform for each step.
\cite{coors2018spherenet} developed a CNN filter in light of spatial location on the sphere to solve the distortion problem of the ERP images. \cite{su2019kernel} proposed a kernel transformer network that converts the pre-trained kernels on the perspective images into ODIs.
SpherePHD~\cite{lee2020spherephd} proposed a convolution kernel applicable to triangular pixels defined on the faces of the icosahedron.
\cite{zhang2019orientation} performed convolution using a hexagonal filter applicable to the vertices of icosahedron.
In this work, we focus on ODI SR and propose SphereSR that applies convolution to the spherical structure created through the subdivision of the icosahedron.

\begin{figure}[t!]
  \centering
  \vspace{-5pt}
  \includegraphics[width=.98\linewidth]{figure/subdivision2.pdf}
   \vspace{-6pt}
   \caption{Subdivision process of icosahedron. We define a pixel as a face of subdivided icosahedron.}
   \label{fig:icosahedron}
   \vspace{-12pt}
\end{figure}

 \vspace{-3pt}
\section{Method}
 \vspace{-3pt}
\noindent\textbf{Overview.}
As shown in Fig.~\ref{fig:overall-framework},
we propose a novel framework, SphereSR, that aims at obtaining a continuous spherical image representation from a given icosahedron input.
First, we introduce a feature extraction method for spherical images that efficiently extracts features from an image on the icosahedron (Sec.~\ref{ds_ks}).
Second, we propose the Spherical Local Implicit Image Function (SLIIF) that can predict RGB values through the extracted features in order to flexibly reconstruct an HR image with arbitrary projection type(Sec.~\ref{image_func}).
Lastly, We propose a feature loss to obtain support from features of other projection types by utilizing the advantage of SLIIF that can be converted to arbitrary projection type(Sec.~\ref{loss_func}).

\subsection{Feature Extraction for Spherical Images}
\label{ds_ks}

Feature extraction is crucial yet challenging for spherical image SR as we focus on very large scale factors, \eg, $\times$16. In this situation, \emph{it is imperative to tackle the memory overload issue while ensuring high SR performance}. For that reason, the proposed SphereSR represents the spherical data based on icosahedron and efficiently extracts features on the spherical surface composed of uniform faces. This is achieved by a new data structure on icosahedron and weight sharing between kernels of different directions.

\noindent \textbf{Data structure.} Inspired by the convolution of icosahedron data in SpherePHD, we propose  a new spherical data structure. To implement the convolution operation, SpherePHD~\cite{lee2020spherephd} uses the subdivision process of icosahedron described in Fig.~\ref{fig:icosahedron} and creates a call-table containing the indices of $N$ neighboring pixels for each pixel and then uses it to stack every neighboring pixel. After that, convolution is performed with the kernel of size $[N+1,1]$.
However, this implementation is not memory-efficient as it requires additional $N$ channels for stacking the neighboring pixels for every convolution operation.
To solve this problem, we propose a new data structure that convolution operation can be directly applied without stacking the neighbors in a call-table.
As shown in the left side of Fig.~\ref{fig:ga-conv}, we rearrange the original data in the direction of the arrow while transforming the triangular pixels to rectangular pixels such that conventional 2D convolution can be applied.
Here, an upward kernel (\textcolor{red}{red kernel}) for the upward($\bigtriangleup$) aligned pixel is arranged in the odd-numbered rows, and a downward kernel (\textcolor{blue}{blue kernel}) for the downward($\bigtriangledown$) aligned pixel is arranged in the even-numbered rows. (More details are in the supplementary.)

\begin{figure}[t]
 \vspace{-5pt}
 \centering
  \includegraphics[width=.98\linewidth]{figure/GA_conv5.pdf}
 \vspace{-5pt}
  \caption{New kernel weight sharing. Left: proposed new data structure, Right: our kernel weight sharing scheme}
  \label{fig:ga-conv}
  \vspace{-15pt}
\end{figure}

\noindent \textbf{Kernel Weight Sharing.  } While the memory overload issue can be resolved by the proposed data structure, it is still necessary to ensure high SR performance. SpherePHD~\cite{lee2020spherephd} rotates the upward or downward kernel by 180$^{\circ}$ to obtain the same kernel shape. Therefore, it is possible to share weights of the up/downward kernels whose directions and shapes are symmetric to each other. However, as the direction of the kernel weight changes for adjacent pixels, high performance cannot be achieved if the characteristics of the texture according to the direction need to be identified.

To solve this problem, we introduce a kernel weight sharing scheme that can geometrically align up/downward directional kernels  \emph{without rotation}.
As shown in the right side of Fig.~\ref{fig:ga-conv}, pixel(face) combinations of two kernels, where up/downward kernels are applied, are shaped differently depending on the direction of a center pixel.
However, if three vertices of the center pixel (denoted as \textcolor{red}{red} and \textcolor{blue}{blue} dots in the right side of Fig.~\ref{fig:ga-conv}) are included in the pixel combination as imaginary pixels, the shapes of two different up/downward pixel combinations can be made to be geometrically identical. To this end, rather than averaging and creating imaginary pixels, we distribute the kernel weight to near pixels. For the upward kernel, image pixel weight $w_3$, $w_9$, $w_{10}$ are distributed to the nearest 4 pixels except the center pixel. For downward kernel weight, $w_4$, $w_5$, $w_11$ are distributed in the same way. The detail can be seen in the right side of Fig.~\ref{fig:ga-conv}.
In this way, the feature extraction module can be applied to any pixel without rotation.

\iffalse
\begin{figure}
  \centering
  \begin{subfigure}{0.49\linewidth}

    \includegraphics[width=1.0\linewidth]{figure/data_structure.pdf}
    \caption{Spherical Data Structure}
    \label{fig:short-a}
  \end{subfigure}
  \hfill
  \begin{subfigure}{0.49\linewidth}

    \includegraphics[width=1.0\linewidth]{figure/sharing3.png}
    \caption{New Kernel weight sharing}
    \label{fig:short-b}
  \end{subfigure}
  \caption{Geometry-aligned Convolutions.}
  \label{fig:ga-conv}
  \vspace{-12pt}
\end{figure}
\fi

\subsection{Spherical Local Implicit Image Function (SLIIF)}

\label{image_func}
\begin{figure}[t!]
  \centering
  \includegraphics[width=.99\linewidth]{figure/SLIIF4.png}
  \vspace{-5pt}
   \caption{Spherical local implicit image function.}
   \label{fig:implict_func}
  \vspace{-10pt}
\end{figure}

\noindent \textbf{Overall Process of SLIIF. }
With the data efficiently represented, we now describe how to super-resolve the ODIs efficiently with an arbitrary scale.
Our main idea is to predict an RGB value for an arbitrary coordinate on the unit sphere $S^2$ using the feature map extracted by using GA-Conv described in Sec.~\ref{ds_ks}.
Inspired by LIIF~\cite{chen2021learning}, we propose SLIIF that learns an implicit image function on $S^2$ using icosahedral faces.
SLIIF takes a spherical coordinate of the point on the unit sphere and its neighboring feature vectors as inputs and predicts the RGB value. It can be formulated as:
\begin{equation}
  I(s) = f_{dec}(z, s),s\in S^2
  \label{sliif}
\end{equation}
where $f_{dec}$ is a decoding function shared with all icosahedral faces, $s$ is the point on the unit sphere $S^2$, $z$ represents a feature vector formed by concatenating neighboring feature vectors of $s$, and $I(s)$ is the predicted RGB value of $s$.

For the pixel in the image that can be formed by any arbitrary projection from the unit sphere, there is a corresponding point $s$ on the unit sphere $S^2$.\footnote{The coordinate of $s$ is computed by using the center point of a pixel.} The face containing $s$ is denoted as $f_s$ and three vertices surrounding $f_s$ are denoted as $v^1_s, v^2_s, v^3_s$ (see Fig.~\ref{fig:implict_func}).
The RGB values of $s$ w.r.t. the coordinate system of three vertices are first calculated and then ensembled based on the triangular areas $A_1, A_2, A_3$ to get the final RGB value of point $s$.
The RGB value of $s$ w.r.t. each vertex $v^j_s$ is calculated with features of 6 faces containing the vertex and relative polar coordinate.
The features of the six faces are concatenated clock-wise starting from $f_s$ to preserve geometrical consistency. Here, we denote the concatenated features as $z_j$ and the polar coordinate of $s$ with respect of $v^j_s$ as $(r_j, \theta_j)$.
To better utilize the positional information, $(r_j, \theta_j)$ are encoded with $\gamma(p) =(\sin{(2^0 \pi p)},\cos{(2^0 \pi p)},...,\sin{(2^{L-1} \pi p)},\cos{(2^{L-1} \pi p)})$ to extend the dimension of a relative coordinate introduced in ~\cite{mildenhall2020nerf,Tancik2020Fourier}.
As such, we can predict the RGB value of point $(\theta,\phi)\in S^2$, which can be formulated as:
\begin{eqnarray}
  I(\theta, \phi) = \sum_{j=1}^3 \frac{A_j}{A} \cdot f_{dec}(z_j, [\gamma (r_j),\gamma(\theta_j)])
\end{eqnarray}
When the pixel in the image corresponds to the vertex on $S^2$, we can still follow the aforementioned procedure because any choice of neighboring vertices results in the same RGB value due to the triangle area-based weighting.

\label{cell_decode}

\iffalse
\begin{table*}[t!]
\centering
\begin{tabular}{cc}
\hline
\multicolumn{1}{c}{ERP} & {Perspective} \\ \hline
  $\begin{pmatrix}
    \Delta x \hat{x} \\ \Delta y \hat{y}
  \end{pmatrix}
  =
  \begin{pmatrix}
     0 & \Delta X \sin{\theta} \\ -\Delta Y & 0
  \end{pmatrix}
  \begin{pmatrix}
    \hat{\theta} \\ \hat{\phi}
  \end{pmatrix}$
&
  $\begin{pmatrix}
    \Delta x \hat{x} \\
    \Delta y \hat{y}
  \end{pmatrix}
  =
  \begin{pmatrix}
    \alpha \Delta X \cos{\theta}\sin{\phi} & \alpha\Delta X \cos{\phi} \\
    -\alpha\Delta Y \sin{\theta} & 0
  \end{pmatrix}
  \begin{pmatrix}
    \hat{\theta} \\
    \hat{\phi}
  \end{pmatrix}$ \\\hline\hline
\multicolumn{2}{c}{Fisheye}              \\ \hline
\multicolumn{2}{c}{
  $\begin{pmatrix}
    \Delta x \hat{x} \\
    \Delta y \hat{y}
  \end{pmatrix}
  =
  \begin{pmatrix}
    \alpha \Delta X \cos{\theta}\sin{\phi} + \alpha\Delta X F_{X}(X,Y)\cos{\theta}\cos{\phi} &
    \alpha\Delta X \cos{\phi} - \alpha\Delta X F_{X}(X,Y)\sin{\phi} \\
    -\alpha\Delta Y \sin{\theta} + \alpha\Delta Y F_{Y}(X,Y)\cos{\theta}\cos{\phi} &
    - \alpha\Delta Y F_{Y}(X,Y)\sin{\phi}
  \end{pmatrix}
  \begin{pmatrix}
    \hat{\theta} \\
    \hat{\phi}
  \end{pmatrix}$

}              \\ \hline
\end{tabular}
\caption{Derivation Results of Projected vectors $\Delta x \hat{x}, \Delta y \hat{y}$ for ERP, Perspective and Fisheye Projections. }
\end{table*}
\fi

\begin{figure}[t]
  \centering
  \includegraphics[width=0.99\linewidth]{figure/cell_decoding2.png}
   \vspace{-10pt}
   \caption{Sphere-oriented Cell Decoding.}
   \vspace{-10pt}
   \label{fig:celldecode}
\end{figure}

\noindent \textbf{Sphere-oriented Cell Decoding. }
Through SLIIF, we can predict the RGB value for any point on $S^2$. That is, we can generate the desired HR image for any projection type by predicting the RGB value of each pixel.
However, SLIIF provides the RGB value only for the center of the pixel and discards the information within the pixel area except for the center value.

To handle this, LIIF~\cite{chen2021learning} defines the cell decoding value as the width and height of the query pixel of interest.
Nonetheless, this definition cannot be directly applied to sphere as the corresponding region on the sphere is not in a rectangular shape, and the direction of the reference vertex, where the RGB value is first calculated, continually changes. To handle this, we propose sphere-oriented cell decoding, a method considering the relation between the pixel region on the projected output and the corresponding region on $S^2$. By adding the cell decoding value as the input to SLIIF, we can fully utilize information within the pixel area.
As shown in Fig.~\ref{fig:celldecode}, we aim at obtaining the RGB value of the rectangular pixel on the projected plane. We call this rectangular pixel a plane cell, which can be expressed using two vectors $\overrightarrow{\Delta X},\overrightarrow{\Delta Y}$. The sphere cell, the corresponding area of the plane cell on the sphere, is located in the face where the corresponding point of pixel center is located on the sphere. The sphere cell can be also expressed using two vectors $\overrightarrow{\Delta x},\overrightarrow{\Delta y}$. The relation between $\overrightarrow{\Delta X},\overrightarrow{\Delta Y}$ and $\overrightarrow{\Delta x},\overrightarrow{\Delta y}$ depends on the projection type and location of the pixel center (refer to the supplementary for details).

For the geometrical consistency of the orders of the concatenated features, relative coordinate of $s$, and cell decoding value among pixels, we need to define new axis vectors, $\overrightarrow{n_1}$ and $\overrightarrow{n_2}$, invariant to face orientation.
The unit vector $\overrightarrow{n_1}$ is defined as a vector between the reference vertex and the face center, and the unit vector $\overrightarrow{n_2}$ is defined by 90 degree counter-clockwise rotation of $\overrightarrow{n_1}$. To get the height and width based on this coordinate system, we approximate the parallelogram sphere cell to the axis-aligned rectangle. The approximated sphere cell is defined as a rectangle which can be expressed using two vectors  $\overrightarrow{\Delta x_{eq}}, \overrightarrow{\Delta y_{eq}}$ that has the same area as the parallelogram sphere cell and the largest intersection area with the parallelogram sphere cell.
Based on the approximated rectangular sphere cell, we finally formulate the sphere-oriented cell decoding value as:

\begin{equation}
  \begin{matrix}
    \begin{pmatrix}
        \overrightarrow{\Delta x} \\ \overrightarrow{\Delta y}
    \end{pmatrix}
    \approx
    \begin{pmatrix}
        \overrightarrow{\Delta x_{eq}} \\ \overrightarrow{\Delta y_{eq}}
    \end{pmatrix}
    =
    \begin{pmatrix}
        c_x\overrightarrow{n_{1}} \\c_y\overrightarrow{n_{2}}
    \end{pmatrix}

  \end{matrix}
\end{equation}

\begin{equation}
    \Rightarrow c = [c_{x}, c_{y}]
    = \left(\frac{|\overrightarrow{\Delta x_{eq}}|} {|\overrightarrow{n_{1}}|},\frac{|\overrightarrow{\Delta y_{eq}}|} {|\overrightarrow{n_{2}}|}\right)
\end{equation}

As a result, we can predict the RGB value $I(X,Y)$ for any point on the projected plane based on the following equation,
\begin{equation}
    \begin{aligned}
        I(X,Y) &= I(\theta,\phi,c)  \\
        &=\sum_{j=1}^3 \frac{A_j}{A} \cdot f_{dec}(z_j, [\gamma (r_j),\gamma(\theta_j)],[c_x,c_y])
    \end{aligned}
\end{equation}

\begin{comment}

\begin{equation}
  \overrightarrow{n_1}=\overrightarrow{{v_x p_x}}=(\theta_x - \theta_v)\hat{\theta}+(\phi_x-\phi_v)\sin{\theta_x}\hat{\phi}
  \label{n1}
\end{equation}

\begin{equation}
  \overrightarrow{n_2}=\overrightarrow{o p_x}\times\overrightarrow{n_1}=-(\phi_x-\phi_v)\sin{\theta_x}\hat{\theta}+(\theta_x - \theta_v)\hat{\phi}
  \label{n2}
\end{equation}

By these equations, we can find the variables that satisfy the following equation,
\begin{equation}
  \begin{pmatrix} \hat{\theta} \\ \hat{\phi} \end{pmatrix}
  =
  \begin{pmatrix} \alpha_1 & \beta_1 \\ \alpha_2 & \beta_2 \end{pmatrix}
  \begin{pmatrix}\overrightarrow{n_1} \\ \overrightarrow{n_2} \end{pmatrix}
    \label{n2}
\end{equation}
By these two vectors $(\overrightarrow{n_1},\overrightarrow{n_2})$, we can define a spherical cell decoding that
\begin{equation}
    [c_{x}, c_{y}]:\text{cell decoding s.t.}
    \begin{cases}
      \Delta x = |c_{x}\overrightarrow{n_{1}}| \\
      \Delta y = |c_{y}\overrightarrow{n_{2}}|
    \end{cases}
  \label{cell_decoding}
\end{equation}

Also,
\begin{equation}
  \begin{matrix}
    \Delta x \hat{x} = \gamma_1 \hat{\theta} + \delta_1 \hat{\phi} \\
    \Delta y \hat{y} = \gamma_2 \hat{\theta} + \delta_2 \hat{\phi}
  \end{matrix}
  \text{ } \Big) \Rightarrow
  \begin{matrix}
    \text{decided by output of} \\
    \text{projection type}
  \end{matrix}
\end{equation}

However, because the cell shape is parallelogram, it is impossible to find accurate values satisfying eq.[5]. So, to find optimal values, we derive an equivalent rectangle A, composed of vectors $\Delta x_{eq} \hat{x}_{eq},\Delta y_{eq} \hat{y}_{eq}$ that it has same area and largest intersection area with the parellelogram B, denoted as $A=\underset{X\in\{X \: \boldsymbol{|}\: |X|=|B|\}}{argmax}|X\cap B|$.

Therefore,

\begin{equation}
  \begin{matrix}
    \begin{pmatrix}
        \Delta x \hat{x} \\ \Delta y \hat{y}
    \end{pmatrix}
    \approx
    \begin{pmatrix}
        \Delta x_{eq} \hat{x}_{eq} \\ \Delta y_{eq} \hat{y}_{eq}
    \end{pmatrix}
    =
    \begin{pmatrix}
        \alpha_1'\overrightarrow{n_{1}} \\\beta_2'\overrightarrow{n_{2}}
    \end{pmatrix}

  \end{matrix}
\end{equation}

\begin{equation}
    \Rightarrow c = [c_{x}, c_{y}]
    = \left(\frac{\Delta x_{eq}}{|\overrightarrow{n_{1}}|},\frac{\Delta y_{eq}}{|\overrightarrow{n_{2}}|}\right)
    = (\alpha_1', \beta_2')
\end{equation}

As a result, we can predict the RGB value $I(X,Y)$ for any point on the projected plane through the following equation,
\begin{equation}
    \begin{aligned}
        I(X,Y) &= I(\theta,\phi,c)  \\
        &=\sum_{j=1}^3 \frac{A_j}{A} \cdot f_{dec}(\mathcal{F}_j, [\gamma (r_j),\gamma(\theta_j)],[c_x,c_y])
    \end{aligned}
\end{equation}

\end{comment}

\subsection{Loss Function}

\begin{figure}
\includegraphics[width=1.0\linewidth]{figure/feature_loss_ERP3.png}
\caption{The proposed feature loss between the spherical and ERP features.}
\label{fig:loss}
 \vspace{-5pt}
\end{figure}

\input{tables/result_table}

\label{loss_func}

We train the proposed framework using two loss terms. First, we use the multi-scale L1 loss. With the L1 loss defined at the multiple scales, our framework can learn more about various relative coordinates and cell decoding values.
Second, we design a feature loss module to measure the similarity between the features extracted from the sphere and other projection types.

As shown in fig.~\ref{fig:loss}, we design a feature mask from the ERP or cube map feature. The spatial part of the mask is generated from the difference between the predicted SR ODI and HR ground truth. The channel part of the mask is generated from the features via channel-wise global average pooling. In this way, we get the feature mask emphasizing the relevant part with high accuracy.
Also, the spherical features are converted to the shapes of other projection types via the SLIIF feature conversion module. Finally, the converted features are then subtracted and masked to formulate the feature loss $L_{feat}$. The total loss is as follows.

\begin{equation}
    Loss = \frac{1}{N} \sum_{j=1}^N \parallel I_j^{est} - I_j^{gt}\parallel_1  + \lambda L_{feat}
\end{equation}

\section{Experiments}

\subsection{Dataset and Implementation}

We train and test SphereSR using the ODI-SR dataset~\cite{deng2021lau} and SUN360 panorama dataset~\cite{xiao2012recognizing}. For training, 750 out of 800 ODI-SR training images are used and the remaining 50 images are used for validation. For testing, we use 100 images from the ODI-SR test dataset and another 100 images from the SUN360 panorama dataset.
The resolution of an HR ODI is 1024$\times$2048, and training is performed for the scale factor $\times$8 and $\times$16.
As shown in Fig.~\ref{fig:icosahedron}, SphereSR takes an image on icosahedron as an input converted from LR ODIs, and the icosahedron subdivision levels for the scale $\times$8 and $\times$16 are set to 5 and 6, respectively. The spherical feature extractor is implemented using GA-conv based on EDSR~\cite{lim2017enhanced} and outputs a feature with 128 dimensions. The RGB prediction decoding function $f_{dec}$ is a 5-layer MLP, and the hidden layer dimension is 256. The EDSR model pre-trained with the ODI-SR dataset is used to extract features from ERP and cube map images for feature loss.
SphereSR is trained for 500 epochs using the Adam Optimizer with batch size 1.
The initial learning rate is set to 0.0001 and reduced by 10 times after 400 epochs.
The feature loss scale parameter $\lambda$ is set to 0 for the first 100 epochs and set to 0.3 for the rest of the training.

\subsection{Evaluation on ERP}

We use the ODI-SR and SUN360 Panorama datasets for evaluation. We compare SphereSR with 9 models for 2D SISR, including SRCNN~\cite{Dong2016Image}, VDSR~\cite{kim2016accurate}, LapSRN~\cite{lai2017deep}, MemNet~\cite{tai2017Memnet}, MSRN~\cite{li2018multiscale}, EDSR~\cite{lim2017enhanced}, D-DBPN~\cite{haris2018ddbpn}, RCAN~\cite{zhang2018image}, EBRN~\cite{qiu2019Embedded} and 2 models for ODI-SR, \ i.e., 360-SS~\cite{ozcinar2019super} and LAU-Net~\cite{deng2021lau}.
We use WS-PSNR~\cite{Zhou2018WeightedtoSphericallyUniformSO} and WS-SSIM~\cite{Zhou2018WeightedtoSphericallyUniformSO} as the evaluation metrics.

\noindent \textbf{Quantitative results.} Table~1 shows the quantitative comparison for $\times$8 and $\times$ 16 SR on the ODI-SR and SUN 360 panorama datasets. As can be seen, SphereSR outperforms all the compared methods on both datasets except the case of $\times$8 SR on the ODI-SR dataset, where SphereSR shows comparable performance with LAU-Net.
However, for $\times$16 SR, SphereSR has the higher performance compared to LAU-Net in WS-PSNR and WS-SSIM performance on both ODI-SR and SUN360 panorama datasets.

\noindent \textbf{Qualitative comparison.} Figure~\ref{fig:ERP_result} shows the visual comparison of the $\times$8 SR images on the ODI-SR dataset. As can be seen, SphereSR reconstructs clear textures and more accurate structures, while other compared methods suffer from the problems of blurred edges or distorted structures.
From the visual comparison, we can conclude that SphereSR produces the texture of repeated patterns more accurately than ERP networks.

\begin{figure*}[]
  \centering
  \includegraphics[width=0.99\linewidth]{figure/ERP_result7.png}
  \vspace{-16pt}
   \caption{Visual comparisons of x8 SR results of different methods on ODI-SR dataset.}
   \label{fig:ERP_result}
 \vspace{-10pt}
\end{figure*}

\subsection{SR for Other  Projection Types}
In this section, we verify whether the proposed SphereSR, trained using the ERP images on the ODI-SR dataset, can perform well for any projection type. We first conduct an experiment on the conversion to FOV 90$^{\circ}$ perspective image with a size 512$\times$512. We then conduct another experiment on the conversion to FOV 180$^{\circ}$ fisheye image with a size 1024$\times$1024. In addition, we use circular fisheye projection, one of several types of fisheye projections.
To compare with other SR models, we use bicubic interpolation to convert to the desired projection type. The ERP GT image is also interpolated using the bicubic method to the desired projection type for performance evaluation. We use PSNR and SSIM as the evaluation metrics. Note that we select five random directions, generate a projection output suitable for the corresponding direction, and calculate the mean value for PSNR and SSIM.

\begin{figure}[t!]
  \centering
  \includegraphics[width=0.99\linewidth]{figure/Perspective_result7.png}
   \vspace{-7pt}
   \caption{Visual comparison for x8 SR of perspective images on ODI-SR dataset.}
   \label{fig:Perspective_result}
   \vspace{-5pt}
\end{figure}

\begin{figure}[t!]
  \centering
  \includegraphics[width=0.99\linewidth]{figure/Fisheye_result7.png}
    \vspace{-7pt}
   \caption{Visual comparison for x8 SR of fisheye images on ODI-SR dataset}
   \label{fig:Fisheye_result}
   \vspace{-12pt}
\end{figure}

\noindent \textbf{Perspective Image. }
Table.~\ref{arbitrary} shows the quantitative results for perspective image SR. SphereSR again achieves the best performance on both ODI-SR and SUN360 datasets. LAU-Net~\cite{deng2021lau} achieves the PSNR values of 26.39dB and 24.33dB on both datasets, respectively. By contrast, our method significantly surpasses LAU-Net and achieves the highest PSNR values of 26.76dB and 24.46dB on both datasets, respectively.
In addition, when removing the SLIIF component in SphereSR, the PSNR values drop by 0.1dB and 0.14dB, respectively.
In Fig.~\ref{fig:Perspective_result}, we show a visual comparison between SphereSR with SLIIF, SphereSR without SLIIF, RCAN~\cite{zhang2018image}  and LAU-Net~\cite{deng2021lau}.
As can be seen, SphereSR reconstructs clear straight lines and textures rather than other compared methods RCAN(b) and LAU-Net(c). Also, as a comparison for the usage of SLIIF(d,e), triangle-shaped artifacts are created when SLIIF is not used(d), but it can be confirmed that clear straight lines are created when SLIIF is used(e).

\noindent \textbf{Fisheye. }
Table.~\ref{arbitrary} shows the quantitative results for fisheye image SR. It can be seen that SphereSR has the highest performance regarding the PSNR and SSIM values on the ODI-SR and SUN360 panorama datasets.
Among the methods for 2D SISR, RCAN achieves the second-highest PSNR value of 24.40dB on the ODI-SR dataset. On the SUN360 panorama dataset, LAU-Net achieves the second-highest PSNR of 24.97dB. Our method shows the highest PSNR and SSIM values, showing the best SR performance.
In Fig.~\ref{fig:Fisheye_result}, we show a visual comparison between SphereSR with SLIIF, SphereSR without SLIIF, RCAN~\cite{zhang2018image} and LAU-Net~\cite{deng2021lau}. Specifically, we crop the area to view the SR results at the south pole.
As can be seen, RCAN(b) and LAU-Net(c) generate inappropriate textures with several lines rushing to the south pole. On the other hand, SphereSR(w/o SLIIF)(d) and SphereSR(w/ SLIIF)(e) do not have such a problem. Moreover, In the case of (e), it eliminates the triangle-shaped artifact generated in (d).

\input{tables/arbitrary_proj_result}

\input{tables/ablation_table}

\input{tables/memory}

\subsection{Ablation Study and Analysis}

In this section, we study the effectiveness for each of our proposed modules, \eg, GA-Conv, SLIIF, and feature loss. In addition, we validate the memory load of CNN implementation between using the data structure in our GA-Conv and the existing SpherePHD~\cite{lee2020spherephd}.

\noindent \textbf{GA-Conv.  }
We compare the results of Models 1 and 3 in Table~\ref{ablation} when adding or removing GA-Conv. GA-Conv is used in the feature extraction module to obtain the feature vector to be used in SLIIF. If GA-Conv is not used, the kernel weight sharing proposed by SpherePHD~\cite{lee2020spherephd}, which gives 180 degrees rotation per kernel, is used. Table~\ref{ablation} shows that using GA-Conv improves PSNR score by $0.14$dB and $0.12$dB in ODI-SR and SUN360 Panorama datasets, respectively, for $\times$8 SR.

\vspace{2pt}
\noindent \textbf{SLIIF.  }
SphereSR uses SLIIF to present SR results for the ERP projection type through feature vectors presented on the sphere. Same as the method mentioned in ~\cite{shi2016real}, we implement a pixel-shuffle algorithm capable of performing SR on the Icosahedron without SLIIF for performance comparison.
When using the pixel-shuffle, the last feature map was subdivided by scale factor multiple, and then the final ERP output was derived using the bicubic interpolation. As can be seen in Models 2 and 3 of Table~\ref{ablation}, using SLIIF for continuous image presentation achieves higher performance (24.34dB vs. 24.31dB) than the pixel-shuffle method of subdivision of icosahedron for $\times$8 SR.

\noindent \textbf{Feature loss.  }
We propose a feature loss that measures the feature similarity of the crucial areas through feature masking using features generated from other projection types. To confirm the effectiveness of feature loss, we compare the SR performing when adding and removing this loss. Models 3 and 4 in Table 3 show the ablation results. Performance comparison indicates a performance improvement for all metrics in $\times$8 and $\times$16 SR.

\noindent \textbf{Data Representation Efficiency.}
In Sec.~\ref{ds_ks}, we point out that the CNN imposition method of SpherePHD is not efficient for SR. We thus propose a new data structure to tackle the problem. To see how efficiently the new data structure is, we implement a simple CNN model and then conduct an experiment comparing the activation memory. The CNN model is simply a structure in which a convolution layer is stacked, the number of convolution layers is set to 16, and the hidden feature dimension is set to 32.
Table.~\ref{memory} shows the experiments from input level 4 to 7. As can be verified, the new data structure in GA-Conv has a much lower activation memory. In addition, it is found that, as the input level increased, the ratio of using new data structure memory to SpherePHD decreased. Based on the result, it can be seen that the proposed data structure is more efficient in terms of memory compared to SpherePHD. Moreover, the efficiency increases as the input resolution increases.

\makeatletter
\newcommand{\thickhline}{%
    \noalign {\ifnum 0=`}\fi \hrule height 1pt
    \futurelet \reserved@a \@xhline
}
\newcolumntype{"}{@{\hskip\tabcolsep\vrule width 1pt\hskip\tabcolsep}}
\makeatother

\vspace{-3pt}
\section{Conclusion}
In this paper, we proposed a novel framework, called SphereSR, to generate a continuous spherical image representation from an LR 360$^{\circ}$ image. SphereSR predicts the RGB values at the given spherical coordinates of an HR image corresponding to an arbitrary project type. We first proposed geometry-aligned convolution to represent spherical data efficiently and then proposed SLIIF to extract RGB values from the spherical coordinates. As such, SphereSR flexibly reconstructed an HR image with an arbitrary projection type and SR scale factors. Experiments on various benchmark datasets showed that our method significantly surpasses existing methods.

\noindent\textbf{Limitation and Future Work}. We focused on finding an efficient data structure and kernel weight sharing method to extract meaningful features with the ODI input based on GA-Conv (Sec.~\ref{ds_ks}).
Future studies will therefore need to improve the network architecture using the properties of ODIs compared to perspective images, and then we can achieve better SR results via SLIIF.

\clearpage
{\small
\bibliographystyle{ieee_fullname}

\section{More Details of Feature Extraction for Spherical Images}
\label{sec:feature}

\subsection{Data structure}

\begin{figure}[b!]
  \centering
  \includegraphics[width=.7\linewidth]{figure/supp_data.png}
  \vspace{-5pt}
   \caption{Proposed new data structure.}
   \label{fig:supp_data}
  \vspace{-10pt}
\end{figure}

As shown in Fig.~\ref{fig:supp_data}, we make the implementation of up/downward kernel possible with 5$\times$3 convolution by using the proposed new data structure. At this time, among the pixels of the5$\times$3 convolution kernel, 5 pixels not included in the actual kernel are always multiplied by weight 0, and the weight update is not performed.

We also always do feature padding with neighborhood pixels right before convolution to avoid losing the continuity of the spherical data structure. The padding size is 2 $\times$ 1 padding considering the kernel size.

Finally, as the convolutional kernels for odd-numbered and even-numbered rows are different, the convolution stride for the row axis is set to 2 such that convolution for the kernel can be performed for each row. After convolution for each kernel, two feature maps are concatenated according to the order of rows.
\subsection{Kernel Weight Sharing}

\begin{figure}[h!]
  \centering
  \includegraphics[width=.62\linewidth]{figure/supp_kernel.png}
  \vspace{-5pt}
   \caption{Geometry-Aligned Convolution(GA-Conv) and equivalent kernel.}
   \label{fig:supp_kernel}
  \vspace{-10pt}
\end{figure}

In this subsection, we describe in detail the derivation of the kernel weight sharing scheme for GA-Conv (Geometry-Aligned Convolution). As shown in Fig.~\ref{fig:supp_kernel}, we set the three vertices of the center pixel to the imaginary pixels $P_a, P_b, P_c$ to make the shape of the pixel combination the same. At this time, each pixel can be regarded as the midpoint of the surrounding four pixels. Therefore, the values of pixels $P_a, P_b, P_c$ are set as the mean value of four pixels through the formulas in the Fig.~\ref{fig:supp_kernel}.
We present a way to utilize weight sharing without actually creating a pixel when implementing this setting. As can be seen from Fig.~\ref{fig:supp_kernel}, the convolution with pixel values and kernel weights is the same as the value convolutional using the equivalent kernel. Therefore, an appropriate weight sharing scheme is used to use GA-Conv, and through this, more sophisticated kernel weight sharing is possible without kernel rotation previously performed in SpherePHD.

\section{Sphere-oriented Cell Decoding}
\label{sec:formatting}

In this section, we explain more details about sphere-oriented cell-decoding method.
A rectangle $\overrightarrow{\Delta X},\overrightarrow{\Delta Y}$ on the projected plane makes parallelogram $\overrightarrow{\Delta x},\overrightarrow{\Delta y}$ on the unit sphere. To define a cell decoding according to point $p$ on unit sphere, we initially define a vector $(\overrightarrow{n_1},\overrightarrow{n_2})$. Let $o = center\;point(r=0)$, $p = (\theta,\phi)$ s.t. point $p$ we want to get RGB and $v_x=(\theta_v,\phi_v)$s.t. one of vertices 1, 2, 3. Then,
\begin{equation}
  \overrightarrow{n_1}=\overrightarrow{{v_x p}}=(\theta - \theta_v)\hat{\theta}+(\phi-\phi_v)\sin{\theta}\hat{\phi}
  \label{n1}
\end{equation}
\begin{equation}
  \overrightarrow{n_2}=\overrightarrow{o p}\times\overrightarrow{n_1}=-(\phi-\phi_v)\sin{\theta}\hat{\theta}+(\theta - \theta_v)\hat{\phi}
  \label{n2}
\end{equation}
To make $(\overrightarrow{n_1},\overrightarrow{n_2})$ as unit vectors, the vectors are scaled as
$\overrightarrow{n_1} \leftarrow \overrightarrow{n_1}/ \parallel\overrightarrow{n_1}\parallel$ and $\overrightarrow{n_2} \leftarrow \overrightarrow{n_2}/ \parallel\overrightarrow{n_2}\parallel$.
By these equations, we can find the variables $\alpha_1, \beta_1, \alpha_2, \beta_2$ that satisfies the following equation,
\begin{equation}
  \begin{aligned}
  \begin{pmatrix} \hat{\theta} \\ \hat{\phi} \end{pmatrix}
  &=
  \begin{pmatrix} \alpha_1 & \beta_1 \\ \alpha_2 & \beta_2 \end{pmatrix}
  \begin{pmatrix}\overrightarrow{n_1} \\ \overrightarrow{n_2} \end{pmatrix}

  \end{aligned}
\label{eq:1}
\end{equation}

To denote vectors $\overrightarrow{\Delta x},\overrightarrow{\Delta y}$ representing parallelogram sphere cells using $\overrightarrow{n_1}, \overrightarrow{n_2}$, we first denote using $\hat{\theta}, \hat{\phi}$ through the projection process of each projection type, as shown in Fig.~\ref{fig:supp_projection}. Therefore, we have to find as the equation,

\begin{equation}
  \begin{matrix}
    \overrightarrow{\Delta x} = \gamma_1 \hat{\theta} + \delta_1 \hat{\phi} \\
    \overrightarrow{\Delta y} = \gamma_2 \hat{\theta} + \delta_2 \hat{\phi}
  \end{matrix}
  \text{ } \Big) \Rightarrow
  \begin{matrix}
    \text{decided by output of} \\
    \text{projection type}
  \end{matrix}
\label{eq:2}
\end{equation}

\input{tables/derivation}

\begin{figure}[t!]
  \centering
  \includegraphics[width=.65\linewidth]{figure/supp_proj2.pdf}
  \vspace{-5pt}
   \caption{Visualization of projection process of each projection type.}
   \label{fig:supp_projection}
  \vspace{-10pt}
\end{figure}

Let P be the point $(X,Y)$ where point p is projected onto the projected plane. As the derivation results, each projection can be organized as a table.~\ref{tab:cell-decoding-tb}. Fisheye projection type is the result of the projection process by setting the reflective surface $Z=F(X,Y)$ as mentioned in ~\cite{Scaramuzza2006a}.
Also, in case of Perspective projection,
\begin{equation}
\begin{aligned}
\begin{matrix}

  \alpha=\frac{1}{ \sqrt {\tan^2{\phi} + \frac{1}{\tan^2{\theta}\cos^2{\phi}}}},
  \Delta X = \frac{2\tan{\frac{M}{2}}}{b},
  \Delta Y = \frac{2\tan{\frac{N}{2}}}{a}

\end{matrix}
\end{aligned}
\end{equation}
in case of Fisheye projection,
\begin{equation}
\begin{aligned}
  \begin{matrix}

  \alpha=\frac{1}{ \sqrt {X^2 + Y^2 + F_{X}(X,Y)^2}} ,
  \Delta X = \frac{F(X',0)}{b} \text{ s.t. } \frac{X'}{F(X',0)}=\tan{\frac{M}{2}},
  \Delta Y = \frac{F(0,Y')}{a} \text{ s.t. } \frac{Y'}{F(0,Y')}=\tan{\frac{N}{2}}

  \end{matrix}
\end{aligned}
\end{equation}

Finally, we are able to represent $\overrightarrow{\Delta x},\overrightarrow{\Delta y}$ using two vectors $\overrightarrow{n_1}, \overrightarrow{n_2}$ by \cref{eq:1,eq:2}.
\begin{equation}
  \begin{aligned}
  \begin{pmatrix} \overrightarrow{\Delta x} \\ \overrightarrow{\Delta y} \end{pmatrix}
  &=
  \begin{pmatrix} \gamma_1 & \delta_1 \\ \gamma_2 & \delta_2 \end{pmatrix}
  \begin{pmatrix} \alpha_1 & \beta_1 \\ \alpha_2 & \beta_2 \end{pmatrix}
  \begin{pmatrix}\overrightarrow{n_1} \\ \overrightarrow{n_2} \end{pmatrix}

  \end{aligned}
\label{eq:sum}
\end{equation}
We then find a  $\overrightarrow{\Delta x_{eq}}, \overrightarrow{\Delta y_{eq}}$ as mentioned in the main paper.

\section{Visualization}
\label{sec:formatting}
\subsection{Additional Qualitative Comparison.}
We additionally visualize the qualitative comparison with Ours and other methods on ODI-SR dataset and SUN360 Panorama dataset in Fig.~\ref{fig:supp_comparison1} and Fig.~\ref{fig:supp_comparison2}. We compare the results with various projection types, ERP, perspective and fisheye.

\subsection{Visualization of SR with Arbitrary Projection Types.}
We additionally visualize the SR results for arbitrary projection type from LR ERP images in Fig.~\ref{fig:supp_various1} and Fig.~\ref{fig:supp_various2}. We also visualize images for different directions to show that our model can predict SR results for arbitrary directions.

\begin{figure}[p]
  \centering
  \includegraphics[width=.9\linewidth]{figure/comparison1_compress.png}
  \vspace{-5pt}
   \caption{Visual comparisons of $\times$8 SR results of ERP, perspective, fisheye projection on SUN 360 Panorama dataset.}
   \label{fig:supp_comparison1}
  \vspace{-10pt}
\end{figure}

\begin{figure}[p]
  \centering
  \includegraphics[width=.9\linewidth]{figure/cherry67_compress.png}
  \vspace{-5pt}
   \caption{Visual comparisons of $\times$8 SR results of ERP, perspective, fisheye projection on ODI-SR dataset.}
   \label{fig:supp_comparison2}
  \vspace{-10pt}
\end{figure}

\begin{figure}[p]
  \centering
  \includegraphics[width=.99\linewidth]{figure/various_result1.png}
  \vspace{-5pt}
   \caption{Visualization of $\times$8 SR results with various projection types on SUN360 Panorama dataset.}
   \label{fig:supp_various1}
  \vspace{-10pt}
\end{figure}

\begin{figure}[p]
  \centering
  \includegraphics[width=.99\linewidth]{figure/various_result2.png}
  \vspace{-5pt}
   \caption{Visualization of $\times$8 SR results with various projection types on ODI-SR dataset.}
   \label{fig:supp_various2}
  \vspace{-10pt}
\end{figure}

{\small
\bibliographystyle{ieee_fullname}


\section{Introduction}

The paper is devoted to modeling spatially inhomogeneous objects with the shape of a ball. For example, this work can be used in planetary sciences or studies focusing on the 3D thermochemical structure of the Earth's mantle. The geophysical and geochemical data suggest that the Earth’s mantle contains heterogeneities on a broad range of scales: from thousands of km down to the grain size \cite{Stixrude:2012}. The stochastic mantle heterogeneity can be attributed to thermal and chemical variations due to mantle convection and melting processes \cite{Allegre:1986}, \cite{Armienti:2010}. As a rule, the complex nature of physical processes, and the lack of measurement data, a fully deterministic description of the problem is impossible. Therefore, efficient statistical modeling methods are required for a better understanding of the formation and preservation of the 3D planetary heterogeneity. This approach involves generating a statistical ensemble of realizations instead of a single deterministic solution. This allows, for example, to estimate the confidence intervals of the desired characteristics, as well as perform uncertainty quantification and sensitivity analyzes (\cite{Urrego-Blanco:2016}; \cite{McClarren:2018}).

Over the past decade, a wide range of mathematical techniques based on the Schoenberg representation of the covariance function has been developed for modeling homogeneous isotropic random fields on a spherical surface (e.g. \cite{Creasey:2018}; \cite{Emery:2019}). On the other hand, volumetric modeling in the full-sphere is not well covered. To our knowledge in the scientific literature, this subject is restricted by the isotropic homogeneous case (\cite{Yaglom:1962}; \cite{Marinucci:2011}). Meschede \& Romanowicz \cite{Meschede:2015} extended the random field simulation technique for the non-homogeneous distribution along the radius using the Karhunen–Loève expansion. However, in our opinion, the disadvantage of this approach is the difficulty in a rigorous mathematical description of the simulated random field. In the presented work, we develop an alternative method for modeling a random field in a ball.

\section{Statistical model}

\subsection{Modeling Gaussian random fields in Cartesian coordinates}

Gaussian random fields are fully determined by mean values and covariance function which in the real case has form
\begin{equation}
\label{eqn:1}
C(\mbox{\boldmath{$x$}},\mbox{\boldmath{$y$}})=E \left \{(f(\mbox{\boldmath{$x$}})-E{f(\mbox{\boldmath{$x$}})})(f(\mbox{\boldmath{$y$}})-E{f(\mbox{\boldmath{$y$}})}) \right \}, \quad \mbox{\boldmath{$x$}}, \mbox{\boldmath{$y$}} \in R^3~.
\end{equation}
Here and below in this work, we use bold font for vectors and matrices, the overbar denotes the complex conjugate. $E$ denotes the mathematical expectation, $f(\mbox{\boldmath{$x$}})$ and $f(\mbox{\boldmath{$y$}})$ are values of random field realization at the coordinates defined by vectors $\mbox{\boldmath{$x$}}$ and $\mbox{\boldmath{$y$}}$ respectively. Further in this work, we deal with zero mean random fields.

In the case when the values of the random field $f(\mbox{\boldmath{$x$}})$ should be generated in the finite number of fixed spatial points, the traditional approach to simulate the realizations is based on the Cholesky decomposition of the covariance matrix $\mbox{\boldmath{$C$}}$
\begin{equation}
\label{eqn:2}
\mbox{\boldmath{$C=LL^T$}}~, \quad \mbox{\boldmath{$f=L\xi$}} ~,
\end{equation}
where $\mbox{\boldmath{$\xi$}}$ is a vector of mutually independent random numbers with zero mean and unit variance, $\mbox{\boldmath{$L^T$}}$ is the transpose of the lower triangular matrix $\mbox{\boldmath{$L$}}$.

Another approach to statistical modeling uses the Karhunen-Loève expansion
\begin{equation}
\label{eqn:3}
\mbox{\boldmath{$C = V^T D V $}} ~, \quad \mbox{\boldmath{$f = V^T \sqrt{D} \xi $}}~,
\end{equation}
where $\mbox{\boldmath{$D$}}$ and $\mbox{\boldmath{$V$}}$ are diagonal matrix of eigenvalues and matrix whose columns are the corresponding eigenvectors of the covariance matrix, respectively. In the truncated Karhunen-Loève expansion, only the eigenvectors corresponding to significant eigenvalues are considered.

In practice, the size of the matrix $\mbox{\boldmath{$C$}}$ is often too large, and decomposition in Eqs. \ref{eqn:2},\ref{eqn:3} is computationally consuming due to the need to perform the Cholesky decomposition or to solve the eigenvalues problem in a high dimension. In application to simulation of the 3D random field in a ball, Meschede and Romanowicz \cite{Meschede:2015} suggested to apply the Karhunen-Loève expansion only in the radial direction to avoid this numerical problem.

\subsection{2D random fields on the sphere}

This subsection is devoted to the modeling of isotropic random fields on the surface of a sphere. Consider the unit sphere $S^2=\{\mbox{\boldmath{$x$}} \in R^3: |\mbox{\boldmath{$x$}}|=1\}$. A basis formed by spherical harmonic functions
\begin{equation}
\label{eqn:4}
Y_n^k (\phi,\theta)=\sqrt{ \frac{2n+1}{4\pi} \frac{(n-k)!}{(n+k)!} } P_n^k (cos(\theta)) e^{ik\phi}
\end{equation}
is used to describe such random fields. Here $P_n^k$ are associated Legendre functions, where $n$ and $k$ are the degree and order of spherical harmonics, respectively, such as $n \ge 0$, $-n \le k \le n$ and, $0 \le \phi < 2\pi$ is longitude and, $0 \le \theta \le \pi$ is colatitude. The $S^2$ metric is the geodesic distance which defines the distance between points $\mbox{\boldmath{$x$}}$ and $\mbox{\boldmath{$y$}}$ as angle
$$
\alpha= arccos \left( \mbox{\boldmath{$x$}} \cdot \mbox{\boldmath{$y$}} \right) ~,
$$
where
$$
\mbox{\boldmath{$x$}}= (cos(\phi_x) sin(\theta_x), sin(\phi_x) sin(\theta_x), cos(\theta_x)) ~,
$$
$$
\mbox{\boldmath{$y$}}= (cos(\phi_y) sin(\theta_y ),sin(\phi_y) sin(\theta_y), cos(\theta_y)) ~.
$$
Then any complex function $f \in L^2 (S^2)$ can be presented as
\begin{equation}
\label{eqn:5}
f(\mbox{\boldmath{$x$}})=f(\phi,\theta)=\sum_{n=0}^\infty \sum_{k=-n}^n f_n^k Y^k(\phi,\theta) ~,
\end{equation}
where due to the orthonormality of the basis of spherical harmonics coefficients $f_n^k$ has the form

$$
f_n^k=\int_0^{2\pi} \int_0^\pi f(\phi,\theta) \overline{Y_n^k(\phi,\theta)} d\theta d\phi .
$$
The spherical harmonics satisfy the important property for two arbitrary points on the spherical surface defined by vectors $\mbox{\boldmath{$x$}}$ and $\mbox{\boldmath{$y$}}$
$$
\frac{4\pi}{2n+1} \sum_{k=-n}^n Y_n^k(\mbox{\boldmath{$x$}}) \overline{Y_n^k(\mbox{\boldmath{$y$}})} =  P_n(\mbox{\boldmath{$x$}} \cdot \mbox{\boldmath{$y$}}),~\mbox{\boldmath{$x$}},\mbox{\boldmath{$y$}} \in S^2
$$
where $P_n=P_n^0$ is a Legendre polynomial.

We consider $f(\mbox{\boldmath{$x$}})$ as a real isotropic random field with zero mean on $S^2$.  $f(\mbox{\boldmath{$x$}})$ is isotropic if the covariance at two points $\mbox{\boldmath{$x$}}$ and $\mbox{\boldmath{$y$}}$
$$
Cov(f(\mbox{\boldmath{$x$}}),f(\mbox{\boldmath{$y$}}))=C_\alpha (\alpha(\mbox{\boldmath{$x$}},\mbox{\boldmath{$y$}})), \mbox{\boldmath{$x$}},\mbox{\boldmath{$y$}} \in S^2
$$
is finite and depends only on the geodesic distance between them. Here $\alpha$ is an angle between points $\mbox{\boldmath{$x$}}$, $\mbox{\boldmath{$y$}}$ on a sphere, i.e. $cos(\alpha)=\mbox{\boldmath{$x$}} \cdot \mbox{\boldmath{$y$}}$.

According to Schoenberg theorem \cite{Schoenberg:1942}, $C_\alpha$ can be a covariance function of an isotropic random field on the sphere if and only if

\begin{equation}
\label{eqn:6}
C_\alpha (\alpha)=\sum_{n=0}^\infty a_n P_n cos(\alpha)~, \alpha \in [0,\pi]
\end{equation}
for some summable series of non-negative coefficients $\{a_n\}$ \cite{Lantuejoul:2019}. A detailed survey of known functions and corresponding decompositions can be found, for example, in \cite{Terdik:2015} and \cite{Lantuejoul:2019}. In practice, the summation from $0$ to a finite spherical harmonic degree number $n_{max}$ is used in Eq. \ref{eqn:6}. The numerical realizations of the random field on the sphere can be obtained using the mathematical apparatus described in \cite{Creasey:2018}, \cite{Emery:2019}, and \cite{Lantuejoul:2019}.

\subsection{3D random fields in the ball}

In this subsection, we consider the random fields distributed in the 3D filled sphere (ball). The modeling approach is based on using the covariance function assumed proportional to the product of the radial and isotropic angular covariance functions
\begin{equation}
\label{eqn:7}
C(\mbox{\boldmath{$x$}},\mbox{\boldmath{$y$}})=E\{f(\mbox{\boldmath{$x$}})f(\mbox{\boldmath{$y$}})\} \sim C_r (r_x,r_y ) C_\alpha(\alpha)~.
\end{equation}
Here \mbox{\boldmath{$x$}} and \mbox{\boldmath{$y$}} the pair of points in spherical coordinates
$$
\mbox{\boldmath{$x$}}= (r_x cos(\phi_x) sin(\theta_x), r_x sin(\phi_x) sin(\theta_x), r_x cos(\theta_x)) ~,
$$
$$
\mbox{\boldmath{$y$}}= (r_y cos(\phi_y) sin(\theta_y), r_y sin(\phi_y) sin(\theta_y ), r_y cos(\theta_y))  ~.
$$
$C_r$ is a covariance function in Cartesian coordinates and $C_\alpha$ is a covariance function on a sphere that has the form Eq. \ref{eqn:7} and
$$
\alpha= arccos \left( \frac{\mbox{\boldmath{$x$}}}{r_x} \cdot \frac{\mbox{\boldmath{$y$}}}{r_y} \right) ~.
$$
This representation of the 3D covariance function with Eq. \ref{eqn:7} allows a rigorous mathematical description of a radially inhomogeneous volumetric random field in a full sphere. In contrast to the statistical method developed by Meschede \& Romanowicz \cite{Meschede:2015}, the angular distribution for the second method is explicitly described by the covariance function $C_\alpha$. Moreover, both the angular and radial covariance functions can be estimated based on regional or global geophysical observations.

\section{Simulation method}

The random field realizations $f(\mbox{\boldmath{$x$}})$ on a sphere can be generated using the spectral method developed in \cite{Lantuejoul:2019}. Simulation formula has the form

\begin{equation}
\label{eqn:8}
f(\mbox{\boldmath{$x$}})=\frac{1}{\sqrt{N}} \sum_{l=1}^N f_l(\mbox{\boldmath{$x$}}) ~.
\end{equation}

The central limit theorem ensures that for a large number of realizations $N$ and independent identically distributed $f_l$, Eq. \ref{eqn:8} generates a Gaussian random field.

Each realization of $f_l$ is generated according to the following algorithm:
\begin{itemize}
\item
generate an integer random number according to distribution ${a_n}$ in Eq. \ref{eqn:6}
$$
n_l \sim {a_n};
$$
\item
generate an integer random number uniformly distributed in the interval $[-n_l,n_l]$
$$
k_l \sim U[-n_l,n_l];
$$
\item
the resulting formula for modeling $f_l$ is
$$
f_l(\mbox{\boldmath{$x$}})=2\sqrt{\pi A} \left( \xi Re(Y_{n_l}^{k_l} (\mbox{\boldmath{$x$}})) + \eta Im(Y_{n_l}^{k_l} (\mbox{\boldmath{$x$}})) \right)  ,
$$
\end{itemize}
where $\xi$ and $\eta$ are the mutually independent random numbers with zero mean and unit variance and $A=\sum_0^{n_{max}} a_n$.

For generalization to the three-dimensional case with covariance function defined by Eq. \ref{eqn:7}, we suggest a method based on the Cholesky decomposition of the radial covariance matrix $\mbox{\boldmath{$C_r=L_r L_r^T$}}$ using Eq. \ref{eqn:2} or Karhunen-Loève expansion $\mbox{\boldmath{$C_r=V_r^T D_r  V_r$}}$ using Eq. \ref{eqn:3}. Matrix $\mbox{\boldmath{$C_r$}}$ is defined as the values of function $C_r$ in Eq. \ref{eqn:7} for all pairs of radius values $r_i,i=1,\ldots,M$.

The resulting simulation formula is
\begin{equation}
\label{eqn:9}
f(r_i,\phi,\theta) = \frac{2\sqrt{\pi A}}{\sqrt{N}} \sum_{l=1}^N \left( \xi(r_i )
Re(Y_{n_l}^{k_l} (\phi,\theta)) + \eta(r_i )
Im(Y_{n_l}^{k_l} (\phi,\theta)) \right) ~,
\end{equation}

where

\begin{equation}
\label{eqn:10}
\mbox{\boldmath{$\xi=L_r \xi_0$}}~, \mbox{\boldmath{$\eta=L_r \eta_0$}}~,
\end{equation}

or

\begin{equation}
\label{eqn:11}
\mbox{\boldmath{$\xi=V_r^T \sqrt{D_r} \xi_0$}}~, \mbox{\boldmath{$\eta=V_r^T \sqrt{D_r} \eta_0$}}~.
\end{equation}

Here, $\mbox{\boldmath{$\xi_0$}}$, $\mbox{\boldmath{$\eta_0$}}$ are the vectors of mutually independent Gaussian random numbers with zero mean and unit variance. According to the central limit theorem, the summation over sufficiently large $N$ in Eq. \ref{eqn:9} provides the Gaussian distribution of the simulated random field. The algorithm described above ensures the reproducing of the covariance function $C(\mbox{\boldmath{$x$}},\mbox{\boldmath{$y$}})$, which proves the positive definiteness of the product in Eq. \ref{eqn:7}.

\section{Numerical results}

In this section, we present the results of statistical modeling performed by the method described above. In this paper, we present numerical examples for the homogeneous exponential covariance function
$$
C_r(d)=C_r(r_x,r_y )=\sigma^2 exp(\frac{-|r_x-r_y |}{I})~,
$$
where $\sigma$ is a standard deviation, and $I$ is a correlation length of the random field. The angular covariance $C_\alpha$ has the form \cite{Terdik:2015}
$$
C_\alpha (\alpha)=\sum_{n=0}^\infty \rho^n P_n (cos(\alpha)) =
\frac{1}{\sqrt{1-2 \rho cos(\alpha)+\rho^2}}, \quad
a_n=\rho^n, \rho \in (0,1)~.
$$
Here the coefficients $a_n$ define the spectrum $\{a_n\}$ in Eq. \ref{eqn:6}.

In Fig.1 we compare the covariance function $C$ calculated on the segment between points $(r_x=0.5,\phi_x=\pi/6,\theta_x=\pi/6)$ and  $(r_y=1,\phi_y=\pi/2,\theta_y=\pi/2)$
with Eqs. \ref{eqn:9}-\ref{eqn:11} and the analytical covariance function Eq. \ref{eqn:7}. The covariance function estimated based on spatial realizations coincides with the analytical curve. It confirmes that the numerical realizations accurately reproduce the assumed statistical model. These results confirm the numerical accuracy of the method described above.

\begin{figure}
	\begin{center}
		\includegraphics[width=.9\textwidth]{fig1.jpg}
		\caption{The covariance function $C$ obtained using Eq. (9) and analytical covariance function (7) for $I=0.15$ and $\rho=0.7$).}\label{fig:fig1}
	\end{center}
\end{figure}

In Fig. 2, we show the realizations of a random field using a 2D cross-section (left plot)  and spherical surface (right plot)  assuming $I=0.05$ and $\rho=0.9$. In Fig.3, we present the 2D cross-sections for different values $I=0.25$ and $\rho=0.9$ (left plot) and $I=0.05$ and $\rho=0.6$ (right plot). Comparison of Figs. 2 and 3 demonstrates the radial influence of parameter $I$ (Fig. 3, left plot) and the angular influence of parameter $\rho$ (Fig. 3, right plot). In particular, Fig. 3 shows that the tuning of the parameters of the second statistical model also allows for modeling anisotropic random fields.

In Fig. 4, we show the 3D realizations of a random field for $I=0.05$ and $\rho=0.7$ using orthographic projection. Eqs. \ref{eqn:9} and \ref{eqn:11} were used to generate the random field realization in Fig. 4. The summation in Eq. \ref{eqn:11} was carried out for the eigenvectors $\mbox{\boldmath{$v_i$}}$ corresponding to $M_1$ largest eigenvalues $\lambda_i$ consisting up 95\% of the trace of the covariance matrix $\mbox{\boldmath{$C_r$}}$
$$
\mbox{\boldmath{$\xi$}}=\sum_{i=1}^{M_1} \sqrt{\lambda_i} \mbox{\boldmath{$v_i$}} \xi_i ~,
$$
where $\{\mbox{\boldmath{$v_i$}}\}$ is the eigenvectors corresponding $\{\lambda_i\}$ and $\{\xi_i\}$ is an independent standard Gaussian random values.

\begin{figure}
	\begin{center}
		\includegraphics[width=.475\textwidth]{fig2a.jpg}
		\includegraphics[width=.475\textwidth]{fig2b.jpg}
		\caption{2D cross-section (left) and spherical surface (right) of 3D realizations of random field related to the statistical model for $I=0.05$ and $\rho=0.9$.}\label{fig:fig2}
	\end{center}
\end{figure}

\begin{figure}
	\begin{center}
		\includegraphics[width=.475\textwidth]{fig3a.jpg}
		\includegraphics[width=.475\textwidth]{fig3b.jpg}
		\caption{2D cross-sections for $I=0.25$ and $\rho=0.9$ (left plot) and $I=0.05$ and $\rho=0.6$ (right plot).}\label{fig:fig3}
	\end{center}
\end{figure}

\begin{figure}
	\begin{center}
		\includegraphics[width=.95\textwidth]{Model2_AngExp_I0.05_sh360.png}
		\caption{3D realizations of random field related for $I=0.05$ and $\rho=0.7$.}\label{fig:fig4}
	\end{center}
\end{figure}

\section{Conclusion}

We propose a new method for statistical modeling of the three-dimensional random field in a unit ball. This approach imposes strict conditions on the form of the covariance function (Eq. \ref{eqn:7}). It allows for a rigorous mathematical description of the simulated random field. The derived modeling formulas are tested by numerical computations. Examples of the generated random field realizations are presented.

\section*{Acknowledgments}
Dmitriy Kolyukhin gratefully acknowledges the financial support from RFBR and NSFB, project number 20-51-18009. Alexander Minakov acknowledges funding through the Research Council of Norway center of excellence funding scheme, project 223272 (The Centre for Earth Evolution and Dynamics). AM also acknowledges the project "3D Earth - A Dynamic Living Planet" funded by ESA as a Support to Science Element (STSE).
In this work, we used toolbox SHBUNDLE (v. 4/11.2018) that was developed by N. Sneeuw, M. Weigelt, B. Devaraju, M. Roth et. al. and provided via download from the Institute of Geodesy (GIS), University of Stuttgart:  https://www.gis.uni-stuttgart.de/en/research/downloads/shbundle/