Linear Spherical Sliced Optimal Transport: A Fast Metric for Comparing Spherical Data

## Abstract

Abstract Efficient comparison of spherical probability distributions becomes important in fields such as computer vision, geosciences, and medicine. Sliced optimal transport distances, such as spherical and stereographic spherical sliced Wasserstein distances, have recently been developed to address this need. These methods reduce the computational burden of optimal transport by slicing hyperspheres into one-dimensional projections, i.e., lines or circles. Concurrently, linear optimal transport has been proposed to embed distributions into L 2 superscript 𝐿 2 L^{2} italic_L start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT spaces, where the L 2 superscript 𝐿 2 L^{2} italic_L start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT distance approximates the optimal transport distance, thereby simplifying comparisons across multiple distributions. In this work, we introduce the Linear Spherical Sliced Optimal Transport (LSSOT) framework, which utilizes slicing to embed spherical distributions into L 2 superscript 𝐿 2 L^{2} italic_L start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT spaces while preserving their intrinsic geometry, offering a computationally efficient metric for spherical probability measures. We establish the metricity of LSSOT and demonstrate its superior computational efficiency in applications such as cortical surface registration, 3D point cloud interpolation via gradient flow, and shape embedding. Our results demonstrate the significant computational benefits and high accuracy of LSSOT in these applications.

## 1 Introduction

Spherical signal analysis has received increasing attention in many domains such as computer vision, astrophysics (Jarosik et al., ), climatology (Garrett et al., ) and neuroscience (Elad & Kimmel, ; Zhao et al., , ). In practice, certain types of data are intrinsically defined on a sphere, namely omnidirectional images, ambisonics, and earth data. In addition, some forms of data are studied canonically in hyperspheres by spherical mappings (Gotsman et al., ; Gu et al., ; Shen & Makedon, ; Chung & Chen, ) from the original data spaces. Among them, a notable category of data is the genus-zero surfaces, which can be conformally mapped to unit spheres according to the Poincare uniformization theorem (Abikoff, ). Cortical surface-based analysis, for example, can be facilitated by the inherent spherical topology (0-genus) of the cerebral cortex. Compared to conventional 3D volumetric analysis, applying spherical representations of cortical surfaces offers a more streamlined and accurate geometric framework to examine the brain. This approach proves particularly advantageous in both cross-sectional and longitudinal studies of brain structure and function (Zhao et al., , ).

Analyzing spherical data through machine learning techniques necessitates metrics that can effectively compare signals while maintaining data structure and geometry. In this context, optimal transport (OT) has emerged as a versatile and widely adopted tool in the machine learning community. OT provides a framework for comparing probability distributions in various domains (Kolouri et al., ; Peyré et al., ; Kolkin et al., ; Arjovsky et al., ; Seguy et al., ; Makkuva et al., ; Rostami et al., ; Cang et al., ; Lu et al., ). OT has demonstrated remarkable effectiveness not only in Euclidean spaces, but its theoretical foundations have also been firmly established for various other manifolds (Fathi & Figalli, ; McCann & Topping, ), especially on spheres (Cui et al., ; Hamfeldt & Turnquist, ). However, the calculation of OT distances is prohibitively expensive with a cubic cost $\mathcal{O}(N^{3})$. This computational burden poses a substantial challenge to applications where the input size $N$ is large. In light of this, researchers have enthusiastically pursued and developed faster alternatives, including entropic optimal transport (Cuturi, ) and sliced optimal transport (Bonnotte, ; Bonneel et al., ; Kolouri et al., , , ; Deshpande et al., ; kolouri2019generalized; Nadjahi et al., ; Bai et al., ; Liu et al., ).

Sliced Optimal Transport (SOT) accelerates computation by projecting high-dimensional data onto one-dimensional subspaces and applying efficient OT solvers for the resulting one-dimensional probability distributions. In this paper, we are interested in sliced optimal transport methods on hyperspheres. Related work includes the Spherical Sliced Wasserstein (SSW) discrepancy (Bonet et al., ), which proposes slicing the sphere by the spherical Radon transform and leveraging circular optimal transport; and the Stereographic Spherical Sliced Wasserstein (S3W) distance (Tran et al., ), which utilizes stereographic projection to solve the SOT problem in Euclidean space with minimal distortion.

Many machine learning applications, from manifold learning to k-nearest neighbor classification and regression, rely on pairwise distance calculations among members of a set, often requiring distance evaluations on the order of the square of the number of objects.
To accelerate pairwise calculation of OT distances between probability measures, Linear Optimal Transport (LOT) (Wang et al., ; Kolouri et al., ; Park et al., ; Bai et al., ) was developed as a framework that embeds distributions into $L^{2}$ spaces where pairwise distance computations can be performed with significantly higher efficiency. For example, Wasserstein Embedding and Sliced Wasserstein Embedding have been applied in Euclidean spaces to enable efficient image analysis/classification Basu et al. (); Ozolek et al. (); Kolouri et al. (); Shifat-E-Rabbi et al. (), faster graph and set learning (Kolouri et al., ; Mialon et al., ; Naderializadeh et al., ), supervised task transferability prediction (Liu et al., ), point cloud retrieval (Lu et al., ), etc. A more recent work, Martin et al. () focused on $\mathbb{S}^{1}$ and established a computationally efficient Linear Circular Optimal Transport (LCOT) framework to embed circular distributions into a $L^{2}$ space.

In this paper, we build on the work of Martin et al. () and extend the linearized methodology from circles to hyperspheres. We develop a novel embedding technique for spherical probability distributions, which significantly enhances the efficiency of group analysis in spherical domains.

Contributions. In summary, our contributions are as follows:

- 1.
We propose Linear Spherical Sliced Optimal Transport (LSSOT) to embed spherical distributions into $L^{2}$ space while preserving their intrinsic spherical geometry.
- 2.
We prove that LSSOT defines a metric, and demonstrate the superior computation efficiency over other baseline metrics.
- 3.
We conduct a comprehensive set of experiments to show the effectiveness and efficiency of LSSOT in diverse applications, from point cloud analysis to cortical surface registration.

## 2 Background

### 2.1 Circular Optimal Transport and Linear Circular Optimal Transport

Consider two circular probability measures $\mu,\nu\in\mathcal{P}(\mathbb{S}^{1})$, where $\mathbb{S}^{1}$ denotes the unit circle in $\mathbb{R}^{2}$. Let us parametrize $\mathbb{S}^{1}$ with the angles in between 0 0 and $1$ and consider the cost function $c(x,y):=h(|x-y|_{\mathbb{S}^{1}})$, where $h:\mathbb{R}\rightarrow\mathbb{R}_{+}$ is a convex increasing function and
$|x-y|_{\mathbb{S}^{1}}:=\min\{|x-y|,1-|x-y|\}$ for $x,y\in[0,1)$. The Circular Optimal Transport (COT) problem between $\mu$ and $\nu$ is defined by the following two equivalent minimization problems:

$$ $COT_{h}(\mu,\nu):=\inf_{\gamma\in\Gamma(\mu,\nu)}\int_{\mathbb{S}^{1}\times \mathbb{S}^{1}}c(x,y)d\gamma(x,y)=\inf_{\alpha\in\mathbb{R}}\int_{0}^{1}h(|F^{ -1}_{\mu}(x)-F^{-1}_{\nu}(x-\alpha)|)dx,$ (1) $$

where in the first expression $\Gamma(\mu,\nu)$ is the set of all couplings between $\mu$ and $\nu$, and in the second expression
$F_{\mu}$ (respectively, $F_{\nu}$) is the cumulative distribution function of $\mu$ on $\mathbb{S}^{1}$
(i.e., $F_{\mu}(y):=\mu([0,y))=\int_{0}^{y}d\mu$, $\forall y\in[0,1)$)
extended on $\mathbb{R}$ by $F_{\mu}(y+n)=F_{\mu}(y)+n,\forall y\in[0,1),n\in\mathbb{Z}$,
and its inverse (or quantile function) is defined by $F^{-1}_{\mu}(y)=\inf\{x:F_{\mu}(x)>y\}$.
When $h(x)=|x|^{p}$ for $1\leq p<\infty$, we denote $COT_{h}(\cdot,\cdot)$ as $COT_{p}(\cdot,\cdot)$, and $COT_{p}(\cdot,\cdot)^{1/p}$ defines a metric on $\mathcal{P}(\mathbb{S}^{1})$.
Moreover, if $\mu=\text{Unif}(\mathbb{S}^{1})$ and $h(x)=|x|^{2}$, the minimizer $\alpha_{\mu,\nu}$ of (1) is the antipodal of $\mathbb{E}(\nu)$, i.e., $\alpha_{\mu,\nu}=\mathbb{E}(\nu)-1/2$.
We refer the reader to Delon et al. (); Rabin et al. (); Bonet et al. () for more details.

Recently, Martin et al. () expanded the linear optimal transport framework (Wang et al., ) to encompass circular probability measures, introducing the Linear Circular Optimal Transport (LCOT) embedding and metric. For the reader’s convenience, we summarize the key definitions from their work below.

The LCOT embedding corresponds to the optimal circular displacement that
comes from the problem of transporting a reference measure to a target measure with respect to $COT_{h}(\cdot,\cdot)$. Precisely,
given a reference measure $\mu\in\mathcal{P}(\mathbb{S}^{1})$
that is absolutely continuous with respect to the Lebesgue measure on $\mathbb{S}^{1}$, the LCOT embedding of a target measure $\nu\in\mathcal{P}(\mathbb{S}^{1})$ is defined by

$$ $\widehat{\nu}^{\mu,h}(x):=F^{-1}_{\nu}(F_{\mu}(x)-\alpha_{\mu,\nu})-x,\quad \forall x\in[0,1),$ (2) $$

where $\alpha_{\mu,\nu}$ is a minimizer of (1). When $h(x)=|x|^{p}$, for $1\leq p<\infty$, this embedding gives rise to the LCOT distance
between circular measures:

$$ $\displaystyle LCOT_{\mu,p}({\nu_{1}},{\nu_{2}})$ $\displaystyle:=\|\widehat{\nu_{1}}^{\mu,h}-\widehat{\nu_{2}}^{\mu,h}\|_{L^{p}( \mathbb{S}^{1},d\mu)},\qquad\forall\nu_{1},\nu_{2}\in\mathcal{P}(\mathbb{S}^{1 }),$ (3) $$

where

$$ $L^{p}(\mathbb{S}^{1},d\mu):=\{f:\mathbb{S}^{1}\to\mathbb{R}:\,\|f\|_{L^{p}( \mathbb{S}^{1},d\mu)}^{p}:=\int_{\mathbb{S}^{1}}|f(t)|_{\mathbb{S}^{1}}^{p}\,d \mu(t)<\infty\}.$ (4) $$

If $\mu=Unif(\mathbb{S}^{1})$, we use the notation $L^{p}(\mathbb{S}^{1}):=L^{p}(\mathbb{S}^{1},d\mu)$. If also $h(x)=|x|^{2}$, the LCOT embedding becomes

$$ $\widehat{\nu}:=F^{-1}_{\nu}\left(x-\mathbb{E}(\nu)+\frac{1}{2}\right)-x,\qquad \forall x\in[0,1),$ (5) $$

and we denote the LCOT distance simply as $LCOT_{2}(\cdot,\cdot)$.

### 2.2 Spherical Slicing

Bonet et al. () recently presented the Spherical Sliced Wasserstein (SSW) discrepancy for probability measures supported on the sphere $\mathbb{S}^{d-1}\subset\mathbb{R}^{d}$, that uses projections onto great circles $C$:

$$ $P^{C}(x):=\operatorname*{arg\,min}_{y\in C}\{\arccos(\langle x,y\rangle)\}, \qquad\forall x\in\mathbb{S}^{d-1},$ (6) $$

where $\arccos(\langle\cdot,\cdot\rangle)$ is the geodesic distance on $\mathbb{S}^{d-1}$. Since each great circle is obtained by intersecting the sphere with a 2-dimensional plane in $\mathbb{R}^{d}$, the authors proposed to average over all great circles by integrating over all planes in $\mathbb{R}^{d}$ (i.e., the Grassmann manifold
(Bendokat et al., )) and then project the distributions onto the intersection with $\mathbb{S}^{d-1}$. However, due to the impracticality of the Grassmann manifold, the set of rank 2 projectors is used instead:

$$ ${V}_{2}(\mathbb{R}^{d}):=\{U\in\mathbb{R}^{d\times 2}:\,U^{T}U=I_{2}\},$ $$

where ${V}_{2}(\mathbb{R}^{d})$ is the Stiefel manifold (Bendokat et al., ).
Hence, the integration process can be performed over $V_{2}(\mathbb{R}^{d})$ according to the uniform distribution $\sigma\in\mathcal{P}(V_{2}(\mathbb{R}^{d}))$. The projection onto $U\in V_{2}(\mathbb{R}^{d})$ can be computed as in (Bonet et al., , Lemma 1) by

$$ $P^{U}(x)=\frac{U^{T}x}{\|U^{T}x\|_{2}},\qquad\forall x\in\mathbb{S}^{d-1} \setminus\{x:U^{T}x=0\}$ (7) $$

Finally, given $\mu,\nu\in\mathcal{P}(\mathbb{S}^{d-1})$, the Spherical Sliced Wasserstein discrepancy between them, for $1\leq p<\infty$, is defined in (Bonet et al., ) by

$$ $(SSW_{p}(\mu,\nu))^{p}:=\int_{V_{2}(\mathbb{R}^{d})}\left(COT_{p}(P_{\#}^{U} \mu,P_{\#}^{U}\nu)\right)^{p}\,d\sigma(U).$ (8) $$

## 3 Method

In this section, we first introduce the Linear Spherical Sliced Optimal Transport (LSSOT) embedding and discrepancy, and then state the main theorem —- LSSOT indeed defines a metric in $\mathcal{P}(\mathbb{S}^{d-1})$.

Figure: Figure 1: Semicircle Transform projects non-polar points onto a great circle, and the north/south poles to everywhere on the circle.
Refer to caption: x1.png

### 3.1 Linear Spherical Sliced Optimal Transport (LSSOT)

###### Definition 3.1 .

Consider spherical measures $\nu_{1},\nu_{2}\in\mathcal{P}(\mathbb{S}^{d-1})$. We define the Linear Spherical Sliced Optimal Transport (LSSOT) embedding of $\nu_{i}$, $i=1,2$, with respect to a slice $U\in V_{2}(\mathbb{R}^{d})$ by

$$ $\widehat{\nu_{i}}^{S}(x,U):=F_{P_{\#}^{U}\nu_{i}}^{-1}\left(P^{U}(x)-\mathbb{E }(P_{\#}^{U}\nu_{i})+\frac{1}{2}\right),\quad x\in\mathbb{S}^{d-1},\,U\in V_{2 }(\mathbb{R}^{d}),$ (9) $$

where $P^{U}_{\#}\nu_{i}$ denotes the pushforward measure on $\mathbb{S}^{1}$
of $\nu_{i}$ by the map $P^{U}$.
We also define their LSSOT discrepancy by

$$ $\displaystyle\left(LSSOT_{2}(\nu_{1},\nu_{2})\right)^{2}$ $\displaystyle:=\int_{V_{2}(\mathbb{R}^{d})}\left(LCOT_{2}(P_{\#}^{U}\nu_{1},P_ {\#}^{U}\nu_{2})\right)^{2}\,d\sigma(U)$ (10) $\displaystyle=\int_{V_{2}(\mathbb{R}^{d})}\|\widehat{\nu_{1}}^{S}(\cdot,U)- \widehat{\nu_{2}}^{S}(\cdot,U)\|_{L^{2}(\mathbb{S}^{1})}^{2}\,d\sigma(U),$ $$

###### Remark 3.2 (LSSOT extends LCOT) .

When $d=2$, $P^{U}=I_{2}$ for every $U\in V_{2}(\mathbb{R}^{2})$
since $V_{2}(\mathbb{R}^{2})$ coincides with the set of orthonormal basis on $\mathbb{R}^{2}$. Thus, for any $\nu\in\mathcal{P}(\mathbb{S}^{1})$ we have that $\widehat{\nu}^{S}(x,U)=\widehat{\nu}^{S}(x,V)$ for all $U,V\in V_{2}(\mathbb{R}^{2})$, so the LSSOT embedding and the LSSOT discrepancy coincide with the LCOT embedding and with the LCOT distance respectively, as the integration is normalized to 1 on $V_{2}(\mathbb{R}^{2})$ (which is isomorphic to the orthogonal group $O(2)$ of dimension 1).

###### Theorem 3.3 (Metric property of LSSOT) .

$LSSOT_{2}(\cdot,\cdot)$ defines a pseudo-metric in $\mathcal{P}(\mathbb{S}^{d-1})$, and a metric when restricting to probability measures with continuous density functions. We will refer to it as $LSSOT$ distance.

We refer the reader to the Appendix for detailed proof. The symmetry of $LSSOT_{2}(\cdot,\cdot)$ is easy to deduce from its definition. The triangle inequality holds as an application of Minkowski inequality by using the fact that $LCOT_{2}(\cdot,\cdot)$ is a metric (see Proposition D.1 in Appendix D). For proving the identity of indiscernibles, the idea is to give an alternative expression of the $LSSOT_{2}(\cdot,\cdot)$ involving a Spherical Radon Transform $\mathcal{R}$: see (C.1) in Appendix C.1, which is inspired by Bonet et al. (). Such formula reads for probability density functions $f_{1},f_{2}\in L^{1}(\mathbb{S}^{d-1})$(^1^11(which is equivalent to considering two absolutely continuous probability measures on $\mathbb{S}^{d-1}$)) as

$$ $\left(LSSOT_{2}(\nu_{1},\nu_{2})\right)^{2}=\int_{V_{2}(\mathbb{R}^{d})}\left( LCOT_{2}(\mathcal{R}f_{1}(U,\cdot),\mathcal{R}f_{2}(U,\cdot))\right)^{2}\,d \sigma(U)$ (11) $$

where $\mathcal{R}:L^{1}(\mathbb{S}^{d-1})\to L^{1}(V_{2}(\mathbb{R}^{d})\times
\mathbb{S}^{1})$ is the integral operator

$$ $\mathcal{R}f(U,z)=\int_{\{x\in\mathbb{S}^{d-1}:\,z=P^{U}(x)\}}f(x)dx,$ $$

where the surface/line integral is taken over a $(d-2)$-dimensional surface, which coincides with a (1-dimensional) great semicircle when $d=3$. We will call this transformation as Semicircle Transform when $d=3$, and $(d-2)$-Hemispherical Transform on $\mathbb{S}^{d-1}$ in general. We describe $\mathcal{R}$ acting on $L^{1}(\mathbb{S}^{2})$ by using different and detailed parameterizations in Appendix A and generalize it for any dimension $d$ in Appendix B. In Appendix C, $\mathcal{R}$ is extended to the set of measures on the sphere, and the disintegration theorem from measure theory is used to obtain the general version of (11).

The novelty is that we will relate $\mathcal{R}$ with an integral transformation presented in Groemer (), instead of relating it with the Hemispherical Transform in Rubin () as done in Bonet et al. () (see Remark B.8 for a discussion in this regard). By doing so, we can prove the injectivity of $\mathcal{R}$ on continuous functions ( Corollary B.7), which is the key property to show that if $LSSOT_{2}(\nu_{1},\nu_{2})=0$ then $\nu_{1}=\nu_{2}$. As a byproduct, we also obtain that $SSW_{p}(\cdot,\cdot)$ is not only a pseudo-metric as proven in Bonet et al. (), but also a metric when it is restricted to probability measures with continuous density functions (Theorem D.3 in Appendix D).

Figure: Algorithm 1 Linear Spherical Sliced Optimal Transport (LSSOT)

###### Definition 3.1 .

###### Remark 3.2 (LSSOT extends LCOT) .

###### Theorem 3.3 (Metric property of LSSOT) .

### 3.2 LSSOT Implementation

Following the Monte-Carlo approach for computing the classical Sliced Wasserstein distance, we consider samples from the uniform distribution
$\sigma\in\mathcal{P}(V_{2}(\mathbb{R}^{d}))$ for generating the slices. To do so, a standard procedure is to first construct matrices in $\mathbb{R}^{d\times 2}$ whose components are drawn from the standard normal distribution $\mathcal{N}(0,1)$, then use the Gram-Schmidt algorithm to orthonormalize the column vectors of the matrix, or apply the QR-decomposition and finally normalize the orthogonal matrix $Q$ (dividing each column vector of $Q$ by its Euclidean norm). In contrast to the Spherical Sliced Wasserstein (SSW) numerical implementation in which projections of poles are not well defined (when $U^{T}x=0$ in 7), LSSOT maps the polar points everywhere on the equator by definition. To achieve this mapping numerically, we put an $\epsilon$-cap around the poles, and redistribute all the mass inside the $\epsilon$-cap evenly onto the discrete points outside the cap. See Figure 1 for a visual depiction of the projection and Algorithm 1 for detailed implementation steps.

### 3.3 Time Complexity and Runtime Analysis

Let $L$ denote the number of slices and $d$ denote the dimension. For a discrete spherical distribution $\nu$ with $N$ samples, the slicing and projection require $\mathcal{O}(LdN)$. For a single slice, the time complexity of calculating the LCOT embedding (Equation 5) of the corresponding circular distributions with respect to the uniform distribution is $\mathcal{O}(N)$ after sorting. Therefore the time complexity of calculating LSSOT embedding of $\nu$ is $\mathcal{O}(LN(d+\log N+1))$, where $LN\log N$ comes from the sorting $\mathbb{S}^{1}$.

For a set of $K$ distributions $\{\nu_{k}\}_{k=1}^{K}$, each with $N$ samples, the pairwise LSSOT distances can be boiled down to $K$ LSSOT embeddings and $\frac{LK(K-1)}{2}$ distances in $L^{2}(\mathbb{S}^{1})$, thus cost $\mathcal{O}(KLN(d+\log N+1+\frac{K-1}{2}))$. However, other methods such as spherical Optimal Transport (OT), Sinkhorn divergence, Sliced Wasserstein Distance (SWD), SSW, and Stereographic Spherical Sliced Wasserstein distance (S3W) all require the time complexity $\mathcal{O}(K^{2}D)$, where $D$ is the complexity of each distance, e.g. $\mathcal{O}(K^{2}N^{3}\log N)$ for spherical OT and $\mathcal{O}(K^{2}LN\log N)$ for SSW, compared to $\mathcal{O}(K^{2}LN)$ for our method. Figure 2 shows the pairwise distances computation time of LSSOT along with the methods mentioned above. All methods are run on the CPU for fair comparison. The distributions are randomly generated with sample sizes $N=1000,5000,10000,12500,15000$. For each fixed sample size, we observe that LSSOT stands out as the fastest in all methods under the increase in the number of distributions. In the Appendix E, we also provide the runtime comparisons for slice-based methods with varying numbers of slices.

Figure: Figure 2: Pairwise distances runtime (log scale) comparison w.r.t the number of distributions. $N$ denotes sample sizes in each distribution. The number of slices is 500 for all slice-based methods.
Refer to caption: extracted/5988566/computation_log.png

## 4 Experiments

To demonstrate that the proposed LSSOT is an effective and efficient metric of spherical distributions while preserving the geometry, we design and implement the following experiments on toy examples, cortical surfaces, and point clouds.

### 4.1 Rotation and Scaling in Spherical Geometry

We test the LSSOT metric on spherical geometric transformations such as rotation and scaling.

For rotations, we generate a set of 20 Von Mises-Fisher (VMF) distributions $\{p^{i}(x;\mu_{i},\kappa)=\frac{\kappa}{2\pi(e^{\kappa}-e^{-\kappa})}\exp(
\kappa\mu_{i}^{T}x)\}_{i=1}^{20}$ with the same concentration parameter $\kappa$, but the mean directions $\mu_{i}$’s are 20 distinct points on the Fibonacci Sphere, i.e. spread evenly on $\mathbb{S}^{2}$. Since the 20 VMFs are of the same shape with varying mean directions, each can be interpreted as a rotated version of another. We calculate the pairwise LSSOT distances among the 20 distributions and visualize the embeddings in 3-dimensional space using multi-dimensional scaling (MDS). We use the Scikit Learn (Pedregosa et al., ) sklearn.manifold.MDS function for MDS implementation. As shown in Figure 3 (left), the embeddings spread almost evenly on the sphere, implying the LSSOT captures the rotation geometry.

To study the scaling geometry, we generate 6 VMF distributions with fixed $\mu$ but varying $\kappa$, and then apply the same LSSOT pairwise distances calculation followed by the MDS embedding procedure. Figure 3 (right) illustrates that the LSSOT metric preserves the linear transitions in scaling geometry, with embeddings falling approximately onto a straight line.

Figure: Figure 3: Generated VMF distributions and manifold learning results by LSSOT. Left: a group of 20 VMFs generated by rotating a source VMF, and the corresponding 3-dimensional visualizations by LSSOT and MDS. Right: 6 VMFs generated by scaling a source VMF with $\kappa=30$, and the same 3-dimensional visualizations.
Refer to caption: x2.png

### 4.2 Cortical Surface Registration

We further verify the validity and efficiency of the proposed LSSOT method in comparing spherical cerebral cortex data for registration tasks. Cortical surface registration seeks to establish meaningful anatomical correspondences across subjects or time points. This process is essential for surface-based analysis in both cross-sectional and longitudinal neuroimaging studies. Cortical surface scans are often acquired from diverse spatial frameworks, necessitating a preliminary alignment process to situate them within a common space. This alignment is crucial for enabling meaningful group analyses and comparisons. The unit sphere $\mathbb{S}^{2}$ is usually chosen to be this common space, owing to the inherent genus-0 spherical topology of the cortical surfaces (Wang et al., ; Gu et al., ; Lui et al., ). Thus, such aligning involves mapping each hemisphere of the convoluted cerebral cortical surface onto the $\mathbb{S}^{2}$ space with minimal distortions. Then the goal is to register spherical cortical surfaces to a template, i.e. computing a smooth velocity field on the sphere which is used to deform the moving vertices into the fixed template. Through this deformation, correspondences are established across a group of cortical surfaces.

#### 4.2.1 Experimental Setup

We leverage the Superfast Spherical Surface Registration (S3Reg) (Zhao et al., ) neural network to perform atlas-based registration, which focuses on registering all surfaces to one atlas surface. S3Reg is an end-to-end unsupervised diffeomorphic registration network based on the Voxelmorph framework (Balakrishnan et al., ) and the Spherical Unet (Zhao et al., ) backbone, offering flexibilities in choosing the similarity measure for spherical surfaces in the training objectives. Following the S3Reg experimental setup with resolution level at the 6-th icosahedron subdivisions (i.e. 40962 vertices on each surface), we replace the original mean squared error (MSE) similarity loss with our LSSOT distance and other baselines, and evaluate the registration performance achieved by each similarity measure. We use the Freesurfer fsaverage (Fischl, ) surface as the fixed (atlas) surface. All models are implemented in PyTorch (Paszke et al., ) and run on a Linux server with NVIDIA RTX A6000 GPU. More experimental details are documented in Appendix F.

Dataset and Preprocessing.
We perform registration on the NKI dataset (Nooner et al., ) and ADNI dataset (Jack Jr et al., ). NKI contains T1-weigted brain MRI images from 515 subjects, which are processed to obtain reconstructed cortical surfaces, spherical mappings, parcellations, and anatomical features using the FreeSurfer pipeline (Fischl, ). A subset of the ADNI dataset that contains 412 subjects is used and processed in the same FreeSurfer pipeline. We choose the sulcal depth as the single-channel feature for registration and normalize features into probability vectors. For parcellations, we choose the DKT Atlas (Klein & Tourville, ; Klein et al., ), which parcellates each hemisphere into 31 cortical regions of interest (ROIs). We split both datasets into 70% scans for training, 10% scans for validation, and 20% scans for testing.

Training Losses.
The training loss has three components: the similarity loss for aligning the moving surfaces to the fixed atlas surface, the Dice loss to impose biological validity, and the regularization loss to enforce the smoothness of the deformation field. Specifically, we employ the parcellation maps generated by FreeSurfer and the torchmetrics.Dice function from the torchmetrics library (Nicki Skafte Detlefsen et al., ) for Dice loss. The training loss is a linear combination of the three components. In our experimental design, we keep the regularization and the Dice terms in all setups. Our primary focus is on investigating the impact of various similarity losses on the registration performance, including LSSOT distance and the following baselines.

Baselines.
Since the original S3Reg framework used MSE as part of the similarity loss, we include MSE as one of the baselines. Besides, Spherical Sliced-Wasserstein (SSW) (Bonet et al., ) and Stereographic Spherical Sliced Wasserstein Distances (S3W) (Tran et al., ) are powerful spherical distances to test LSSOT against. We use the SSW implementation from the official repo (^2^22https://github.com/clbonet/Spherical_Sliced-Wasserstein) and the Amortized Rotationally Invariant Extension of S3W (AI-S3W) from the official implementation (^3^33https://github.com/mint-vu/s3wd) with 10 rotations. We also add the Slice Wasserstein distance (SWD), which is defined in Euclidean space $\mathbb{R}^{3}$ instead of $\mathbb{S}^{2}$. The implementation of the SWD algorithm is from the Python Optimal Transport (POT) library (Flamary et al., ).

Metrics
We quantitatively evaluate the similarities between the registered moving surfaces and the fixed atlas surface using mean absolute error (MAE) and Pearson correlation coefficient (CC) on features, as well as LSSOT, SWD on the spherical distributions (represented by vertex locations and feature probability vectors). We calculate the Dice score of the registered parcellations and the fixed parcellation. Moreover, we measure the area (resp. edge) distortions between the original moving surface meshes and the deformed surface meshes as relative changes in the areas of all faces (resp. the lengths of all edges) of meshes.

#### 4.2.2 Results

Table 4.2.2 summarizes the performance of registration evaluated by the metrics for four scenarios: left/right hemispheres of the NKI and ADNI datasets. As the two best methods overall, our proposed LSSOT metric and SSW are on par with each other, yet LSSOT costs less than 10% of the training time. We emphasize that LSSOT is the leading method in the alignment of the parcellations (i.e. Dice Score), indicating that LSSOT yields more biologically meaningful registrations. Besides, with the least area/edge distortions, LSSOT finds the simplest solution of the deformation field to align the moving surface to the fixed surface. Figure 4 shows the registered sulcal depth map and parcellation map of a random subject. More visualizations can be found in the Appendix G. It’s worth noting that the LSSOT embedding $\hat{v}^{S}_{\text{fix}}$ of the fixed surface is calculated once and for all before training, and in each training iteration, only the LSSOT embedding of the deformed moving surface $\hat{v}^{S}_{\text{mov}}$ is calculated and then compared to $\hat{v}^{S}_{\text{fix}}$. Other competing methods, however, need to process both surfaces and calculate the distance during each iteration, which accounts for LSSOT taking less time than SWD. This makes LSSOT a perfect fit for atlas-based registrations, as the template is fixed.

**Table 1: Evaluation metrics on the test datasets (mean $\pm$ standard deviation) for sulcal depth registration of both hemispheres from NKI dataset and ADNI dataset. The best two methods are bolded. The average training time for each epoch is also included in the bottom row of each scenario. We would like to note that as the number of samples is 40,962, using spherical optimal transport or even Sinkhorn divergence is computationally intractable in both time and space.**
$$ $\displaystyle\mathcal{L}=\frac{1}{I}\sum_{i=1}^{I}\Bigl{(}\|\Psi_{\text{dec}}( \Phi_{\text{enc}}(X_{i}))-X_{i}\|^{2}_{2}+\lambda\cdot LSSOT_{2}(\Phi_{\text{ enc}}(X_{i}),X_{\text{unif}})\Bigr{)}$ (12) $$

### 4.3 Point Cloud Interpolation

LSSOT can also be utilized in point cloud analysis, once each point cloud is endowed with a spherical representation. In this experiment, we explore the interpolations between point cloud pairs from the ModelNet dataset (Wu et al., ). Specifically, we train an autoencoder to project the original point clouds to a spherical latent space, to represent each of them as a spherical distribution. Then we apply gradient flow between the pairs of spherical distributions using the LSSOT metric along with SSW and Spherical OT. Finally, this transformation is reconstructed in the original space by the trained decoder, resulting in an interpolation between the original pairs. Figure 5 (top) illustrates this process.

#### 4.3.1 Mapping Point Clouds to Spherical Latent Space via LSSOT Autoencoder

As point clouds are set-structured data, we leverage the Set Transformer (Lee et al., ) architecture to build the autoencoder(^4^44https://github.com/juho-lee/set_transformer). The encoder $\Phi_{\text{enc}}$ is composed of three consecutive Set Attention Blocks, followed by a linear layer and a normalization layer to output a set of points on $\mathbb{S}^{2}$. The decoder $\Psi_{\text{dec}}$ consists of four Set Attention Blocks with a linear layer at the end. For a set of point clouds $\{X_{i}\}_{i=1}^{I}$ (each point cloud $X_{i}$ is a set of points), the autoencoder is trained to reconstruct each $X_{i}$ under a soft regularity constraint in the latent space. That is, the objective function is defined as

$$ $\displaystyle\mathcal{L}=\frac{1}{I}\sum_{i=1}^{I}\Bigl{(}\|\Psi_{\text{dec}}( \Phi_{\text{enc}}(X_{i}))-X_{i}\|^{2}_{2}+\lambda\cdot LSSOT_{2}(\Phi_{\text{ enc}}(X_{i}),X_{\text{unif}})\Bigr{)}$ (12) $$

where $X_{\text{unif}}$ is a set of samples from the uniform distribution on $\mathbb{S}^{2}$.

After training, the encoder $\Phi_{\text{enc}}$ can serve as a spherical mapping from $\mathbb{R}^{3}$ to $\mathbb{S}^{2}$, and the decoder $\Psi_{\text{dec}}$ acts as an inverse mapping. Enforcing the latent representation to closely approximate the uniform distribution offers significant benefits. It minimizes the occurrence of singularities in both the spherical mapping and its inverse function, resulting in fewer distortions or discontinuities, hence the structural and geometrical characteristics of the original point clouds are better preserved. In implementation, we select $\lambda=1e-4$ through hyperparameter search. We set the number of slices to be $L=1000$, and the reference size as $M=1024$ for the $LSSOT_{2}$ function.

#### 4.3.2 Gradient Flow and Interpolations of Spherical Representations

Given two spherical distributions as the source and the target, we can attach gradients on the source distribution and drive it toward the target by minimizing the distance between them. For the distance, we use LSSOT, SSW, and Spherical OT for comparisons. Here the Riemannian Gradient flow algorithm (see Appendix I.3) is implemented on $\mathbb{S}^{2}$. With a slight abuse of notation to mimic the constant speed geodesic, we sample 5 time points denoted as $T=0,1/4,2/4,3/4,1$, from the beginning ($T=0$) to the loss converging ($T=1$), such that the loss goes down by approximately the same amount in each interval. In Figure 5 (bottom), we visualize the interpolations for two example pairs of point clouds.

Figure: Figure 5: Top panel: the interpolation process between a pair of point clouds using gradient flow in the latent space $\mathbb{S}^{2}$. Bottom panel: gradient flow interpolations from a range hood to a bottle (left), and from a stool to a chair (right) using three metrics LSSOT, SSW and spherical OT. See Appendix I.4 for more interpolated pairs.
Refer to caption: x4.png

## 5 Limitations and Future Work

To the best of our knowledge, this work is the first to adapt the Linear Optimal Transport methodology onto hyperspheres, enabling rapid comparisons among a group of spherical distributions. However, spherical data/signals may not come in the form of distributions in practice, and normalizing data into probability distributions may result in the loss of information, thus one limitation of our work is that LSSOT is defined only for probability measures. Immediate future work includes a partial extension of LSSOT for positive signals with unequal mass, and transport $L^{p}$ Thorpe et al. (); ; Crook et al. () extension for non-positive spherical signals, both of which will be studied on $\mathbb{S}^{1}$ as the first step, and extended to $\mathbb{S}^{2}$ by slicing.

## 6 Conclusion

In this work, we propose Linear Spherical Sliced Optimal Transport (LSSOT) as a novel and fast metric for defining the similarities between spherical data. We provide rigorous mathematical proofs establishing LSSOT as a true metric. We also demonstrate the computation efficiency of the LSSOT, superior to other baselines in the setting of pairwise comparisons of spherical data. We verify that LSSOT captures the rotation and scaling geometry on spheres through experiments on the Von Mises-Fisher distributions. In the task of atlas-based cortical surface registration on various datasets, LSSOT delivers top registration performance across multiple evaluation metrics while reducing the training time significantly. Furthermore, LSSOT is shown to yield meaningful interpolations between point cloud pairs. The metric’s effectiveness in both point cloud analysis and practical neuroimaging applications demonstrates its broad potential in 3D data analysis.

## 7 Acknowledgement

SK acknowledges support from NSF CAREER Award #2339898, CC acknowledges funding from NIH Grant P50MH109429, and BAL acknowledges funding from NIH 1R01EB017230. The authors also thank Shiyu Wang and Karthik Ramdas for processing the NKI and ADNI datasets and providing the processed files to the research team.

Data collection and sharing for this project was funded (in part) by the Alzheimer’s Disease Neuroimaging Initiative (ADNI) (National Institutes of Health Grant U01 AG024904),
DOD ADNI (Department of Defense award number W81XWH-12-2-0012). ADNI is funded by the National Institute on Aging, the National Institute of Biomedical Imaging and Bioengineering, and through generous contributions from the following: AbbVie, Alzheimer’s Association; Alzheimer’s Drug Discovery Foundation; Araclon Biotech; BioClinica, Inc.; Biogen; Bristol-Myers Squibb Company; CereSpir, Inc.; Cogstate; Eisai Inc.; Elan Pharmaceuticals, Inc.; Eli Lilly and Company; EuroImmun; F. Hoffmann-La Roche Ltd and its affiliated company Genentech, Inc.; Fujirebio; GE Healthcare; IXICO Ltd.; Janssen Alzheimer Immunotherapy Research & Development, LLC.; Johnson & Johnson Pharmaceutical Research & Development LLC.; Lumosity; Lundbeck; Merck & Co., Inc.;
Meso Scale Diagnostics, LLC.; NeuroRx Research; Neurotrack Technologies; Novartis Pharmaceuticals Corporation; Pfizer Inc.; Piramal Imaging; Servier; Takeda Pharmaceutical Company; and Transition Therapeutics. The Canadian Institutes of Health Research is providing funds to support ADNI clinical sites in Canada. Private sector contributions are facilitated by the Foundation for the National Institutes of Health ([www.fnih.org](www.fnih.org)). The grantee organization is the Northern California Institute for Research and Education, and the study is coordinated by the Alzheimer’s Therapeutic Research Institute at the University of Southern California. ADNI data are disseminated by the Laboratory for Neuro Imaging at the University of Southern California.

## Appendix A Spherical Radon transform in 𝕊 2 superscript 𝕊 2 \mathbb{S}^{2} blackboard_S start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT

In Bonet et al. (), a Spherical Radon transform is introduced and it is used to give an equivalent formula of the SSW discrepancy (8).
In this section, we will discuss such Spherical Radon transform and its relation with the Semicircle Transform.

### A.1 Notations on Sphere 𝕊 2 superscript 𝕊 2 \mathbb{S}^{2} blackboard_S start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT

Let $U\in V_{2}(\mathbb{R}^{3})$, where

$$ $V_{2}(\mathbb{R}^{3})=\{U\in\mathbb{R}^{3\times 2}:\,U^{T}U=I_{2}\}.$ $$

Each $U$ determines a $2$-dimensional plane spanned by the columns of $U$:

$$ $H(U):=\text{span}\{U[:,1],U[:,2]\}\subset\mathbb{R}^{3}.$ $$

The set $\{U[:,1],U[:,2]\}$ is an orthonormal basis for $H(U)$. The subspace $H(U)$ can be equivalently defined by a normal unit vector
$n_{U}$, which is the uniquely determined by $U$ (up to scalar $\pm 1$):

$$ $\displaystyle n_{U}\in\mathbb{S}^{2},\qquad U^{T}n_{U}=0_{2}.$ (13) $$

Thus, the plane $H(U)$ can be defined as

$$ $H(U)=n_{U}^{\perp}:=\{x\in\mathbb{R}^{3}:\,x^{T}n_{U}=0\}.$ $$

Note that $\mathcal{S}(U):=H(U)\cap\mathbb{S}^{2}$
defines a great circle.

### A.2 Projection onto the great circle 𝒮 ⁢ ( U ) 𝒮 𝑈 \mathcal{S}(U) caligraphic_S ( italic_U )

Given $U\in V_{2}(\mathbb{R}^{3})$, the operator $UU^{T}:\mathbb{R}^{3}\to\mathbb{R}^{3}$ defines a projection. Indeed, for each point $s\in\mathbb{S}^{2}\setminus\{\pm n_{U}\}$, the mapping

$$ $\displaystyle s\mapsto\mathcal{P}^{1}_{\mathcal{S}(U)}(s):=\frac{UU^{T}s}{\|UU ^{T}s\|}=\frac{UU^{T}s}{\|U^{T}s\|}=\frac{U[:,1]^{T}s}{\|U^{T}s\|}U[:,1]+\frac {U[:,2]^{T}s}{\|U^{T}s\|}U[:,2]\in\mathcal{S}(U)\subset\mathbb{R}^{3}.$ (14) $$

defines a projection onto the great circle $\mathcal{S}(U)$, where in (14) we used that

$$ $\left\|UU^{T}s\right\|=s^{T}UU^{T}UU^{T}s=s^{T}UU^{T}s=\left\|U^{T}s\right\|$ $$

(see (Bonet et al., , Lemma 1)).
Note that in Bonet et al. (), the above projection is written as

$$ $\displaystyle s\mapsto P^{U}s:=\frac{U^{T}s}{\|U^{T}s\|},$ (15) $$

which is a map from $\mathbb{S}^{2}\setminus\{\pm n_{U}\}$ to $\mathbb{R}^{2}$ that saves the coordinates $[U[:,1]^{T}s,U[:,2]^{T}s]$.
This function maps $s$ into a 1-dimensional circle $\approx\mathbb{S}^{1}$, which can be regarded as a representation/parametrization of the great circle $\mathcal{S}(U)$.

An equivalent method that allows us to define such a projection is introduced by the following parameterization system Quellmalz et al. ():

Each $s\in\mathbb{S}^{2}$ can be written in spherical coordinates:

$$ $s=[\cos\alpha\sin\beta,\sin\alpha\sin\beta,\cos\beta]^{T}.$ $$

where $\alpha\in[0,2\pi)$ is called the azimuth angle of $s$, and $\beta\in[0,\pi]$ the zenith angle of $s$. We will use the notation $\alpha=azi(s)$, $\beta=zen(s)$.

###### Remark A.1 .

When $s=\pm[0,0,1]$, the above parametrization is not unique. As in Quellmalz et al. (), we set $\alpha=0,\beta=0$ for $s=[0,0,1]$ and
$\alpha=0,\beta=\pi$ for $s=[0,0,-1]$.
Thus,

$$ $s\mapsto(\alpha,\beta)$ $$

becomes a bijection (see (Quellmalz et al., , page 6)).

Let us define the following mappings:

$$ $\displaystyle[0,2\pi)\times[0,\pi]\ni(\alpha,\beta)$ $\displaystyle\mapsto\Phi(\alpha,\beta):=[\cos\alpha\sin\beta,\sin\alpha\sin \beta,\cos\beta]^{T}=:s\in\mathbb{S}^{2}$ (16) $\displaystyle\mathbb{S}^{2}\ni s$ $\displaystyle\mapsto azi(s)=:\alpha\in[0,2\pi)$ (17) $\displaystyle\mathbb{S}^{2}\ni s$ $\displaystyle\mapsto zen(s)=:\beta\in[0,\pi]$ (18) $$

In addition, we define the following matrix:

$$ $\displaystyle\Psi(\alpha,\beta,\gamma):=\begin{bmatrix}\cos\alpha&-\sin\alpha& 0\\ \sin\alpha&\cos\alpha&0\\ 0&0&1\end{bmatrix}\begin{bmatrix}\cos\beta&0&\sin\beta\\ 0&1&0\\ -\sin\beta&0&\cos\beta\end{bmatrix}\begin{bmatrix}\cos\gamma&-\sin\gamma&0\\ \sin\gamma&\cos\gamma&0\\ 0&0&1\end{bmatrix}$ $$

By (Quellmalz et al., , page 7) and Gräf & Potts (),

$$ $\displaystyle\mathbb{S}^{2}\times[0,2\pi)\ni(s,\gamma)\mapsto\Psi({azi}(s),{ zen}(s),\gamma)\in SO(3)$ (19) $$

is a bijection, where $SO(3)$ is the special orthogonal group in
$\mathbb{R}^{3}$:

$$ $\displaystyle SO(3)=\{G\in\mathbb{R}^{3\times 3}:\,GG^{T}=G^{T}G=I_{3},\,\text {det}(G)=1\}.$ $$

In the special case $\gamma=0$, we adopt the following notation:

$$ $\displaystyle\Psi(s):=\Psi(\alpha,\beta,0)$ $\displaystyle=\begin{bmatrix}\cos\alpha\cos\beta&-\sin\alpha&\cos\alpha\sin \beta\\ \sin\alpha\cos\beta&\cos\alpha&\sin\alpha\sin\beta\\ -\sin\beta&0&\cos\beta\end{bmatrix}$ $$

where, given $s\in\mathbb{S}^{2}$, $\alpha$ and $\beta$ are as in (17) and (18) (or equivalently, given $(\alpha,\beta)$, $s\in\mathbb{S}^{2}$ is defined by (16)).

It is straightforward to verify

$$ $\Psi(s)^{T}s=\mathfrak{n}:=[0,0,1]^{T}.$ $$

That is, we can regard $\Psi(s)^{T}$ as a rotation matrix that can rotate $s$ back to the North Pole $\mathfrak{n}=[0,0,1]^{T}$. Equivalently, $\Psi(s)$ is a rotation matrix, which can rotate $[0,0,1]^{T}$ to $s$.

Define the Equator by

$$ $\mathcal{E}:=\{[\cos\gamma,\sin\gamma,0]:\,\gamma\in[0,2\pi)\}.$ $$

Given $U\in V_{2}(\mathbb{R}^{3})$, consider the following set

$$ $\{s\in\mathbb{S}^{2}:\Psi(n_{U})^{T}s\in\mathcal{E}\}$ $$

to obtain the following characterization of great circles for $\mathbb{S}^{2}\subset\mathbb{R}^{3}$:

###### Lemma A.2 .

Given $U\in V_{2}(\mathbb{R}^{3})$, it defines a great circle $\mathcal{S}(U)=H(U)\cap\mathbb{S}^{2}$ which is characterized by the following identity:

$$ $\mathcal{S}(U)=\{s\in\mathbb{S}^{2}:\Psi(n_{U})^{T}s\in\mathcal{E}\}.$ (20) $$

###### Proof.

First, since $\Psi(n_{U})$ is a rotation matrix, $\Psi(n_{U})^{T}s\in\mathbb{S}^{2}$ for all $s\in\mathbb{S}^{2}$, as $\|\Psi(n_{U})^{T}s\|=\|s\|=1$.

Now, write $n_{U}=[\cos\alpha^{*}\sin\beta^{*},\sin\alpha^{*}\sin\beta^{*},\cos\beta^{*}]^
{T}$, where $\alpha^{*}:=azi(n_{U})$ and $\beta^{*}=zen(n_{U})$.

For any $s\in\mathcal{S}(U)$, we have $(n_{U})^{T}s=0$.

We have to prove that the third coordinate of $\Psi(n_{U})^{T}s\in\mathbb{S}^{2}$ is 0 0:

$$ $\displaystyle[\Psi(n_{U})^{T}s][3]$ $\displaystyle=[\cos\alpha^{*}\sin\beta^{*},\sin\alpha^{*}\sin\beta^{*},\cos \beta^{*}]s=(n_{U})^{T}s=0.$ $$

That is, for every $s\in\mathcal{S}(U)$, we have obtained that $s\in\{s\in\mathbb{S}^{2}:\Psi(n_{U})^{T}s\in\mathcal{E}\}$.

Similarly, if we consider $s\in\{s\in\mathbb{S}^{2}:\Psi(n_{U})^{T}s\in\mathcal{E}\}$, we have
$0=[\Psi(n_{U})^{T}s][3]=(n_{U})^{T}s$. Thus, $s\in\mathcal{S}(U)$.

Therefore, we have the equality in (20).
∎

###### Lemma A.3 .

The mapping

$$ $[0,2\pi)\ni\gamma\mapsto\Phi(\gamma,\pi/2)=[\cos\gamma,\sin\gamma,0]^{T}\in \mathcal{E}$ $$

and the mapping

$$ $\mathcal{E}\ni s\mapsto\Psi(n_{U})s\in\mathcal{S}(U)$ $$

are bijections.
Thus, the mapping

$$ $[0,2\pi)\ni\gamma\mapsto s:=\Psi(n_{U})\Phi(\gamma,\pi/2)\in\mathcal{S}(U)$ $$

is also a bijection.

###### Proof.

It is trivial to verify that the first mapping, $\gamma\mapsto[\cos\gamma,\sin\gamma,0]^{T}$,
is a bijection.

By the Lemma A.2, the mapping $s\mapsto\Psi(n_{U})^{T}s$ is well-defined from $\mathcal{S}(U)$ to $\mathcal{E}$.
Since $\Psi(n_{U})^{T}$ is a rotation matrix, it is bijection $\mathbb{S}^{2}\to\mathbb{S}^{2}$ whose inverse is $\Psi(n_{U})$, and so
the mapping $s\mapsto\Psi(n_{U})s$ is a well-defined bijection from $\mathcal{E}$ to $\mathcal{S}(U)$: that is, if we restrict $\Psi(n_{U})$ to $\mathcal{E}$, by Lemma A.2, the matrix $\Psi(n_{U})$ defines a bijection from $\mathcal{E}$ to $\text{Range}(\Psi(n_{U})\mathcal{E})=\mathcal{S}(U)$.
∎

Now, we define the following composition map $\mathcal{P}^{2}_{\mathcal{S}(U)}(s)$:

$$ $\displaystyle\mathbb{S}^{2}\setminus\{\pm n_{U}\}\ni s\mapsto\gamma:={azi}( \Psi(n_{U})^{T}s)\mapsto\mathcal{P}^{2}_{\mathcal{S}(U)}(s):=\Psi(n_{U})\Phi( \gamma,\pi/2)\in\mathcal{S}(U)$ (21) $$

###### Proposition A.4 .

The mappings $\mathcal{P}^{1}_{\mathcal{S}(U)}$ (defined by 14) and $\mathcal{P}^{2}_{\mathcal{S}(U)}$ (defined by 21) coincide on $\mathbb{S}^{2}\setminus\{\pm n_{U}\}$.

###### Proof.

Case 1: Suppose $n_{U}=[0,0,1]^{T}$. Consider $s=[\cos\alpha\sin\beta,\sin\alpha\sin\beta,\cos\beta]^{T}\in\mathbb{S}^{2}
\setminus\{\pm n_{U}\}$, i.e., where $\alpha=azi(s)$ and $\beta=zen(s)$ with $\beta\in(0,\pi)$.
In this
case, from (13), $U[3,:]=0_{2}$. Thus, $\{U[1:2,1],U[1:2,2]\}$ forms an orthornomal basis for $\mathbb{R}^{2}$.
From (14), we have

$$ $\displaystyle UU^{T}s=\begin{bmatrix}I_{2}&0\\ 0&0\end{bmatrix}s=[\cos\alpha\sin\beta,\sin\alpha\sin\beta,0]^{T}$ $$

Thus,

$$ $\mathcal{P}^{1}_{\mathcal{S}(U)}(s)=\frac{UU^{T}s}{\|U^{T}s\|}=[\cos\alpha, \sin\alpha,0]^{T}.$ $$

In addition, in this case, $\Psi(n_{U})=I_{3}$. Thus, $\gamma:=azi(\Psi(n_{U})^{T}s)=azi(s)=\alpha$. Therefore,

$$ $\mathcal{P}^{2}_{\mathcal{S}(U)}(s)=\Psi(n_{U})\Phi(\gamma,\pi/2)=\Phi(\gamma, 0)=[\cos\alpha,\sin\alpha,\pi/2]^{T}.$ $$

Thus, 14 and 21 coincide in this case.

Case 2: Suppose $n_{U}\neq[0,0,1]$.
Each point $s\in\mathbb{S}^{2}$ can be written as $s=\Psi(n_{U})s^{\prime}$ where $s^{\prime}=\Psi(n_{U})^{T}s$. In addition, when $s\neq\pm n_{U}$, we can parametrize $s^{\prime}$ via
$s^{\prime}=[\cos\alpha\sin\beta,\sin\alpha\sin\beta,\cos\beta]^{T}$, where $\beta\in(0,\pi)$.
We also have that $\Psi(n_{U})^{T}n_{U}=[0,0,1]^{T}$.

Let $U^{\prime}=\Psi(n_{U})^{T}U$. We have

$$ $(U^{\prime})^{T}[0,0,1]^{T}=U^{T}\Psi(n_{U})[0,0,1]^{T}=U^{T}n_{U}=0_{2}.$ $$

Thus $n_{U^{\prime}}=[0,0,1]^{T}$.
Besides, we have that

$$ $\displaystyle UU^{T}s=\Psi(n_{U})U^{\prime}(U^{\prime})^{T}\Psi(n_{U})^{T}\Psi (n_{U})s^{\prime}=\Psi(n_{U})U^{\prime}(U^{\prime})^{T}s^{\prime}.$ $$

Thus,

$$ $\mathcal{P}^{1}_{\mathcal{S}(U)}(s)=\Psi(n_{U})\mathcal{P}^{1}_{\mathcal{S}(U^ {\prime})}(s^{\prime}),\qquad\forall s\neq\pm n_{U}.$ $$

In addition, by using the parametrization of $s^{\prime}$ in spherical coordinates, we have

$$ $\displaystyle\gamma:=azi(\Psi(n_{U})^{T}s)=azi(s^{\prime})=\alpha.$ $$

Thus,

$$ $\displaystyle\mathcal{P}^{2}_{\mathcal{S}(U)}(s)=\Psi(n_{U})[\cos\alpha,\sin \alpha,0]^{T}=\Psi(n_{U})\mathcal{P}^{2}_{\mathcal{S}(U^{\prime})}(s^{\prime}).$ $$

From Case 1, we have

$$ $\mathcal{P}^{1}_{\mathcal{S}(U^{\prime})}(s^{\prime})=\mathcal{P}^{2}_{ \mathcal{S}(U^{\prime})}(s^{\prime}),$ $$

and from the the fact that $s^{\prime}\mapsto s=\Psi(n_{U})s^{\prime}$ is a bijection on $\mathbb{S}^{2}$, we get

$$ $\mathcal{P}^{1}_{\mathcal{S}(U)}(s)=\mathcal{P}^{2}_{\mathcal{S}(U)}(s).$ $$

∎

###### Lemma A.5 .

Given $z\in\mathbb{S}^{1},U\in V_{2}(\mathbb{R}^{3})$, there exists
$\alpha^{*}\in[0,2\pi),\beta^{*}\in[0,\pi],\gamma\in[0,2\pi)$ such that

$$ $\displaystyle n_{U}$ $\displaystyle=\Phi(\alpha^{*},\beta^{*})=[\cos\alpha^{*}\sin\beta^{*},\sin \alpha^{*}\sin\beta^{*},\cos\beta^{*}]^{T}$ (22) $\displaystyle Uz$ $\displaystyle=\Psi(n_{U})[\cos\gamma,\sin\gamma,0]^{T}$ (23) $$

Vice-versa: Given $\alpha^{*}\in[0,2\pi),\beta^{*}\in[0,\pi],\gamma\in[0,2\pi)$, there exist $z\in\mathbb{S}^{2},U\in V_{2}(\mathbb{R}^{3})$ for which formulas 22 and 23 hold.

###### Proof.

Given $U\in V_{2}(\mathbb{R}^{3})$, one can find $n_{U}\in\mathbb{S}^{2}$ via (13). (Note that if $n_{U}$ satisfies (13), then $-n_{U}$ also satisfies (13), and so $n_{U}$ is unique up to sign.)
Then, once $n_{U}$ is chosen,
we can uniquely find $(\alpha^{*},\beta^{*})$ as (17), (18).

In addition, for $z\in\mathbb{S}^{1}\subset\mathbb{R}^{2}$, since

$$ $\mathcal{P}^{1}_{\mathcal{S}(U)}(Uz)=\frac{UU^{T}Uz}{\|U^{T}Uz\|}=\frac{Uz}{\| z\|}=Uz,$ $$

we obtain that $Uz\in\mathcal{S}(U).$
By Lemma A.3,
$\Psi(n_{U})^{T}Uz\in\mathcal{E}$. Then, we can uniquely determine $\gamma$ (see Lemma A.3).
Thus, (22), and (23) are satisfied.

For the other direction, consider angles $\alpha^{*},\beta^{*},\gamma$ and let us separate the analysis into the following cases:

Case 1: $\alpha^{*}=0,\beta^{*}\in\{0,\pi\}$. In this case, we have $n_{U}=[0,0,\pm 1]$ and thus we can set

$$ $U:=\begin{bmatrix}I_{2}\\ 0\end{bmatrix}.$ $$

In addition, set $z:=[\cos\gamma,\sin\gamma]^{T}\in\mathbb{S}^{1}$ and so the identities 22, and 23 are satisfied.

Case 2: If $\alpha^{*}\neq 0$, then we set $n_{U}$ by 22. Choose an orthornormal basis $\{u_{1},u_{2}\}$ for $n_{U}^{\perp}=\{x\in\mathbb{R}^{3}:\,x^{T}n_{U}=0\}$, and set $U:=[u_{1},u_{2}]$. Then $(U,n_{U})$ satisfies 13.

By Lemma A.3, we have $\Psi(n_{U})[\cos\gamma,\sin\gamma,0]^{T}\in\mathcal{S}(U)$. Then,

$$ $UU^{T}\Psi(n_{U})[\cos\gamma,\sin\gamma,0]^{T}\in\mathcal{S}(U),$ $$

and by setting

$$ $z:=U^{T}\Psi(n_{U})[\cos\gamma,\sin\gamma,0]^{T},$ $$

we have

$$ $\displaystyle Uz$ $\displaystyle=UU^{T}\Psi(n_{U})[\cos\gamma,\sin\gamma,0]^{T}=\frac{UU^{T}\Psi( n_{U})[\cos\gamma,\sin\gamma,0]^{T}}{\|U^{T}\Psi(n_{U})[\cos\gamma,\sin\gamma, 0]^{T}\|}$ $\displaystyle=\mathcal{P}^{1}_{\mathcal{S}(U)}(\Psi(n_{U})[\cos\gamma,\sin \gamma,0]^{T})=\Psi(n_{U})[\cos\gamma,\sin\gamma,0]^{T}$ $$

Thus, 23 is satisfied.
∎

Now, given $\gamma\in[0,2\pi)$, and $n_{U}=\Phi(\alpha^{*},\beta^{*})\in\mathbb{S}^{2}$, we can define the following semicircle:

$$ $\displaystyle\mathcal{SC}^{\gamma}_{n_{U}}:=\mathcal{SC}^{\gamma}_{\alpha^{*}, \beta^{*}}:=\{s\in\mathbb{S}^{2}:\,{azi}(\Psi(n_{U})^{T}s)=\gamma\}.$ (24) $$

###### Remark A.6 .

When $\gamma\neq 0$, we can replace the above definition by

$$ $\displaystyle\mathcal{SC}^{\gamma}_{n_{U}}:=\mathcal{SC}^{\gamma}_{\alpha^{*}, \beta^{*}}:=\{s\in\mathbb{S}^{2}\setminus\{\pm n_{U}\}:\,{azi}(\Psi(n_{U})^{T} s)=\gamma\}.$ (25) $$

###### Proposition A.7 .

Let $z\in\mathbb{S}^{1},U\in V_{2}(\mathbb{R}^{3})$. Consider
$\alpha^{*}\in[0,2\pi),\beta^{*}\in[0,\pi],\gamma\in[0,2\pi)$ as in Lemma A.5 (or equivalently, first, fix $(\alpha^{*},\beta^{*})\in((0,2\pi)\times[0,\pi])\cup\{0\}\times\{0,\pi\}),
\gamma\in[0,2\pi)$, and then, find $z,U$ according to Lemma A.5). We have:

$$ $\mathcal{SC}_{n_{U}}^{\gamma}\setminus\{\pm n_{U}\}=\mathcal{SC}_{\alpha^{*}, \beta^{*}}^{\gamma}\setminus\{\pm n_{U}\}=\{s\in\mathbb{S}^{2}\setminus\{\pm n _{U}\}:\,P^{U}s=z\}.$ $$

###### Proof.

We have that

$$ $\displaystyle\{s\in\mathbb{S}^{2}\setminus\{\pm n_{U}\}:\,P^{U}s=z\}=\left\{s \in\mathbb{S}^{2}\setminus\{\pm n_{U}\}:\,\frac{UU^{T}s}{\|U^{T}s\|}=Uz\right\}$ $\displaystyle=\left\{s\in\mathbb{S}^{2}\setminus\{\pm n_{U}\}:\,\mathcal{P}^{1 }_{\mathcal{S}(U)}(s)=Uz\right\}$ $\displaystyle=\left\{s\in\mathbb{S}^{2}\setminus\{\pm n_{U}\}:\mathcal{P}^{2}_ {\mathcal{S}(U)}(s)=\Psi(n_{U})[\cos\gamma,\sin\gamma,0]^{T}\right\}\quad\text {(by Proposition \ref{pro:proj_12} and (\ref{eq:params_2}))}$ $\displaystyle=\left\{s\in\mathbb{S}^{2}\setminus\{\pm n_{U}\}:\Psi(n_{U})\Phi( {azi}(\Psi(n_{U})^{T}s),\pi/2)=\Psi(n_{U})[\cos\gamma,\sin\gamma,0]^{T}\right \}\quad\text{(by definition \ref{eq:proj_2})}$ $\displaystyle=\left\{s\in\mathbb{S}^{2}\setminus\{\pm n_{U}\}:\Phi({azi}(\Psi( n_{U})^{T}s),\pi/2)=[\cos\gamma,\sin\gamma,0]^{T}\right\}\quad\text{(applying $\Psi^{T}(n_{U})$ on both sides)}$ $\displaystyle=\left\{s\in\mathbb{S}^{2}\setminus\{\pm n_{U}\}:{azi}(\Psi(n_{U} )^{T}s)=\gamma\right\}\qquad\text{(by Lemma \ref{lem:gamma-S(U)})}$ $\displaystyle=\mathcal{SC}_{n_{U}}^{\gamma}\setminus\{\pm n_{U}\}=\mathcal{SC} _{\alpha^{*},\beta^{*}}^{\gamma}\setminus\{\pm n_{U}\}.$ $$

∎

###### Remark A.8 .

With the notation of the above proposition. For each fixed pair $(z,U)\in\mathbb{S}^{2}\times V_{2}(\mathbb{R}^{3})$, if we choose distinct $(\alpha^{*},\beta^{*},\gamma)$ and $(\tilde{\alpha}^{*},\tilde{\beta}^{*},\tilde{\gamma})$ such that the conditions in Lemma A.5 are satisfied, we obtain

$$ $\mathcal{SC}_{\alpha^{*},\beta^{*}}^{\gamma}=\mathcal{SC}_{\tilde{\alpha}^{*}, \tilde{\beta}^{*}}^{\tilde{\gamma}}.$ $$

###### Remark A.1 .

###### Lemma A.2 .

###### Proof.

###### Lemma A.3 .

###### Proof.

###### Proposition A.4 .

###### Proof.

###### Lemma A.5 .

###### Proof.

###### Remark A.6 .

###### Proposition A.7 .

###### Proof.

###### Remark A.8 .

### A.3 Semicircle Transform and Spherical Radon transform

Let $f\in L^{1}(\mathbb{S}^{2})$. The Normalized Semicircle Transform of $f$ is a function $\mathbb{S}^{2}\to\mathbb{R}$ defined as:

$$ $\displaystyle\mathcal{W}f(\alpha^{*},\beta^{*},\gamma):=\frac{1}{4\pi}\int_{ \mathcal{SC}^{\gamma}_{n_{U}}}f(s)\sin({zen}(s))d\mathcal{SC}_{n_{U}}^{\gamma}$ (26) $$

where $d\mathcal{SC}^{\gamma}_{n_{U}}$ is the curve differential on the semicircle $\mathcal{SC}^{\gamma}_{n_{U}}$.

Its un-normalized version is defined as

$$ $\displaystyle\mathcal{UW}f(\alpha^{*},\beta^{*},\gamma):=\frac{1}{4\pi}\int_{ \mathcal{SC}^{\gamma}_{n_{U}}}f(s)d\mathcal{SC}_{n_{U}}^{\gamma}$ (27) $$

The Spherical Radon Transform of $f\in L^{1}(\mathbb{S}^{d-1})$, for $d\geq 2$, defined in Bonet et al. () is the mapping $\mathcal{R}f:V_{2}(\mathbb{R}^{d})\times\mathbb{S}^{1}\mapsto\mathbb{R}$ given by:

$$ $\displaystyle\mathcal{R}f(U,z)=\int_{\mathbb{S}^{d-1}}f(s)\delta(z=P^{U}(s))ds$ (28) $$

When considering $d=3$, we obtain the following relation between the above-defined transforms.

###### Proposition A.9 .

The Unbounded Semicircle Transform 27 and the Spherical Radon Transform 28 are equivalent. In particular, given $z\in\mathbb{S}^{1},U\in V_{2}(\mathbb{R}^{3})$, consider $\alpha^{*}\in[0,2\pi),\beta^{*}\in[0,\pi],\gamma\in[0,2\pi)$ defined as in Lemma A.5, then

$$ $\mathcal{UW}f(\alpha^{*},\beta^{*},\gamma)=\frac{1}{4\pi}\mathcal{R}f(U,z).$ $$

###### Proof.

From Proposition A.7, we have
$\mathcal{SC}_{\alpha^{*},\beta^{*}}^{\gamma}=\{s:P^{U}s=z\},$ thus

$$ $\displaystyle\mathcal{R}f(U,z)=\int_{\mathcal{SC}_{n_{U}}^{\gamma}}f(s)d \mathcal{SC}_{n_{U}}^{\gamma}=|4\pi|\mathcal{UW}(f)(\alpha^{*},\beta^{*}, \gamma).$ $$

∎

###### Remark A.10 .

We notice that there are other versions of the Spherical Radon Transform. For example, in Groemer ()
the spherical Radon transformation
is defined by integrating over the full great circles $\mathcal{S}(U)$

$$ $\widetilde{\mathcal{R}}f(n_{U})=\int_{\mathcal{S}(U)}f\,d\mathcal{S}(U).$ $$

In this paper we are using the variant given in Bonet et al. (), which we have proven coincides with the integral transformation over semicircles. That is, this version coincides with the integral transformation denoted by $\mathcal{B}$ in Groemer (). It holds that

$$ $\widetilde{\mathcal{R}}(n_{U})=\mathcal{R}f(U,z)+\mathcal{R}f(U,-z)$ $$

We will review this in Appendix B when generalizing from $d=3$ to higher dimensions.

We end this section with a discussion about some alternative formulations for $\mathcal{W}f$ and $\mathcal{UW}f$.

For each $s\in\mathcal{SC}^{\gamma}_{\alpha^{*},\beta^{*}}\setminus\{\pm[0,0,1]^{T}\}$, we have

$$ ${azi}(\Psi(\alpha^{*},\beta^{*},0)^{T}s)=\gamma.$ $$

Thus, we can write $s$ as

$$ $s=\Psi(\alpha^{*},\beta^{*},0)\Phi(\gamma,\xi)=\Psi(\alpha^{*},\beta^{*}, \gamma)\Phi(0,\xi),$ $$

for some $\xi\in(0,\pi)$ (recall that $\Phi(0,\xi)=[\sin(\xi),0,\cos(\xi)]^{T}$, and see Groemer () for more details).

Similar to Remark A.6, we can redefine the semicircle as

$$ $\displaystyle\mathcal{SC}_{\alpha,\beta}^{\gamma}=\{\Psi(\alpha,\beta,\gamma) \Phi(0,\xi):\xi\in(0,\pi)\},$ $$

and $\Psi(\alpha^{*},\beta^{*},\gamma)$ is called primal median of $\mathcal{SC}_{\alpha^{*},\beta^{*}}^{\gamma}$.

Combining this with the fact that $(\alpha^{*},\beta^{*},\gamma)\mapsto\Psi(\alpha^{*},\beta^{*},\gamma)$ is a bijection, we can rewrite (26), and (27) as

$$ $\displaystyle\mathcal{W}f(Q)$ $\displaystyle=\frac{1}{4\pi}\int_{0}^{\pi}f(Q\Phi(0,\xi))\sin(\xi)d\xi$ $\displaystyle\mathcal{UW}f(Q)$ $\displaystyle=\frac{1}{4\pi}\int_{0}^{\pi}f(Q\Phi(0,\xi))d\xi$ $$

where $Q\in SO(3)$ (by using the identity $s=Q\Phi(0,\xi)$, for an appropriate $Q\in SO(3)$).

In Hielscher et al. (), the Un-normalized Semicircle Transform is defined as

$$ $\displaystyle\widetilde{\mathcal{UW}}(f)(Q):=\int_{-\pi/2}^{\pi/2}f(Q^{T}\Phi( \xi,\pi/2))d\xi,\qquad Q\in SO(3).$ (29) $$

The equivalence relation between $\mathcal{UW}(f)$ and $\widetilde{\mathcal{UW}}(f)$ reads as follows:

###### Lemma A.11 .

For any $Q\in SO(3)$, we have

$$ $\displaystyle\mathcal{UW}f(Q)=\widetilde{\mathcal{UW}}f\left(\begin{bmatrix}1& 0&0\\ 0&0&-1\\ 0&1&0\end{bmatrix}Q^{T}\right).$ $$

###### Proof.

Given $Q\in SO(3),\xi\in(0,\pi)$, let $Q^{\prime}=\begin{bmatrix}1&0&0\\
0&0&-1\\
0&1&0\end{bmatrix}Q^{T},\xi^{\prime}=\xi-\beta$.
Then,

$$ $\displaystyle Q^{\prime}[1,:]$ $\displaystyle=Q^{T}[1:]=Q[:,1]$ $\displaystyle Q^{\prime}[2,:]$ $\displaystyle=-Q^{T}[3,:]=-Q[:,3]$ $\displaystyle Q^{\prime}[3,:]$ $\displaystyle=Q^{T}[2,:]=Q[:,2]$ $$

and so

$$ $\displaystyle Q\Phi(0,\xi)$ $\displaystyle=Q[\sin\xi,0,\cos\xi]^{T}$ $\displaystyle=Q[:,1]\sin\xi+Q[:,3]\cos\xi$ $\displaystyle=Q[:,1]\cos(\xi^{\prime})-Q[:,3]\sin(\xi^{\prime})$ $\displaystyle=Q^{\prime}[1,:]\cos\xi^{\prime}+Q^{\prime}[2,:]\sin\xi^{\prime}$ $\displaystyle=(Q^{\prime})^{T}[\cos\xi^{\prime},\sin\xi^{\prime},0]^{T}$ $\displaystyle=(Q^{\prime})^{T}\Phi(\xi^{\prime},\pi/2)$ $$

Thus, $\mathcal{UW}f(Q)=\widetilde{\mathcal{UW}}f(Q)$.
∎

Another formulation has been introduced in Groemer (). Since it allows us to define $\mathcal{UW}f$ in higher dimensions ($d\geq 3$), we will discuss it in the next section.

###### Proposition A.9 .

###### Proof.

###### Remark A.10 .

###### Lemma A.11 .

###### Proof.

## Appendix B Spherical Radon Transform in 𝕊 d − 1 superscript 𝕊 𝑑 1 \mathbb{S}^{d-1} blackboard_S start_POSTSUPERSCRIPT italic_d - 1 end_POSTSUPERSCRIPT

### B.1 Great circle and semicircle / ( d − 2 ) 𝑑 2 (d-2) ( italic_d - 2 ) -dimensional hemisphere in 𝕊 d − 1 superscript 𝕊 𝑑 1 \mathbb{S}^{d-1} blackboard_S start_POSTSUPERSCRIPT italic_d - 1 end_POSTSUPERSCRIPT

We set $d\geq 3$ and consider

$$ $V_{2}(\mathbb{R}^{d}):=\{U\in\mathbb{R}^{d\times 2}:\,U^{T}U=I_{2}\}.$ $$

Similar to the previous section, we can define 2D hyperplane

$$ $H(U):=\text{span}\{U[:,1],U[:,2]\},$ $$

and $U^{\perp}:=H(U)^{\perp}$.

In addition, we define the corresponding great circle:
$\mathcal{S}(U)=H(U)\cap\mathbb{S}^{d-1}$.
Note that the Equator is a specific great circle:

$$ $\mathcal{E}=\{s\in\mathbb{S}^{d-1}:s^{T}\mathfrak{n}=0\},$ $$

where $\mathfrak{n}=[0,0,\ldots,0,1]^{T}\in\mathbb{S}^{d-1}\subset\mathbb{R}^{d}$.
Similar to the previous section, the corresponding projection onto the great circle is

$$ $\displaystyle\mathcal{P}^{1}_{\mathcal{S}(U)}(s)=\frac{UU^{T}s}{\|U^{T}s\|}, \qquad s\in\mathbb{S}^{d-1}\setminus U^{T}$ (30) $$

Similar to Proposition A.7 and (24), we define the $(d-2)$-dimensional hemisphere:

$$ $\displaystyle\mathcal{HS}(U,z):=\{s:P^{U}s=z\}\cup(U^{\perp}\cap\mathbb{S}^{d- 1})$ (31) $$

###### Remark B.1 .

When $d=3$, $U^{\perp}=\{\pm n_{U}\}$ and if we ignore these two points, (31) is consistent to (24).

By introducing the hemisphere, the projection 30 can be regarded as the following:
Let $s\in\mathbb{S}^{d}\setminus U^{\perp}$, we have: $s\in\mathcal{HS}(U,z)$ where $z=\frac{U^{T}s}{\|U^{T}s\|}$.
Then the function $\mathcal{P}^{1}_{\mathcal{S}(U)}$ maps $s$ to the intersection point $\mathcal{HS}(U,z)\cap\mathcal{SC}(U)$.

###### Proposition B.2 .

The following identity holds for great circles:

$$ $\displaystyle\mathcal{S}(U)=\mathcal{P}^{1}_{\mathcal{S}(U)}(\mathbb{S}^{d} \setminus U^{\perp})=\{Uz:z\in\mathbb{S}^{1}\}$ (32) $$

###### Proof.

Since $\mathcal{P}^{1}_{\mathcal{S}(U)}$ is the projection, we have
$\mathcal{P}^{1}_{\mathcal{S}(U)}(\mathbb{S}^{d}\setminus U^{\perp})\subset
\mathcal{S}(U)$.

For the other direction, choose $s\in\mathcal{S}(U)$ and we decompose $s$ as

$$ $s=s_{U}+s_{U^{\perp}},$ $$

where $s_{U}\in H(U),s_{U^{\perp}}\in U^{\perp}$.

By definition of $\mathcal{S}(U)$, we have $s_{U^{\perp}}=0$. Thus $s=s_{U}$. From the fact $H(U)\cap U^{\perp}=\{0_{d}\}$, we have $s\notin U^{\perp}$. Therefore, $s\in\mathbb{S}^{d-1}\setminus U^{\perp}$. Thus
$s=\mathcal{P}^{1}_{\mathcal{S}(U)}(s)$. Therefore, we have

$$ $\mathcal{S}(U)\subset\mathcal{P}^{1}_{\mathcal{S}(U)}(\mathbb{S}^{d}\setminus U ^{\perp}).$ $$

It remains to verify the second equality. Pick $z\in\mathbb{S}^{1}$, we have

$$ $\displaystyle\mathcal{P}^{1}_{\mathcal{S}(U)}(Uz)=\frac{UU^{T}Uz}{\|U^{T}Uz\|} =Uz$ $$

Thus $Uz\in\mathcal{S}(U)$.
For the other direction, pick $s\in\mathcal{S}(U)$, thus

$$ $s=s_{U}=UU^{T}s$ $$

In addition, $1=\|s\|=\|UU^{T}s\|=\|U^{T}s\|$.
Thus, $s\in\{Uz:z\in\mathbb{S}^{1}\}$ and the second equality holds.
∎

###### Remark B.1 .

###### Proposition B.2 .

###### Proof.

### B.2 Great circles and hemispheres in Groemer ( 1998 )

In this subsection, we re-parametrize the great circles and semicircles by using the notation in Groemer ().

Let $\mathcal{B}^{d}=\{(u,v):\,u,v\in\mathbb{S}^{d-1},u^{T}v=0\}.$
We define, as in Groemer (), the hyperplane $H(u,v)=\text{span}\{u,v\}$ and the corresponding “great circle”

$$ $\mathcal{S}(u,v)=\mathbb{S}^{d-1}\cap H(u,v).$ $$

The $(d-2)$-dimensional hemisphere is defined by

$$ $\displaystyle\mathcal{HS}(u,v)=\{s\in\mathbb{S}^{d-1}:s^{T}u=0,s^{T}v\geq 0\}$ (33) $$

###### Lemma B.3 .

Pick $U\in V_{2}(\mathbb{R}^{d})$, $z\in\mathbb{S}^{1}$, there exists $(u,v)\in\mathcal{B}^{d}$ such that

$$ $\displaystyle\text{span}(U)$ $\displaystyle=\text{span}(\{u,v\})$ (34) $\displaystyle U^{T}v$ $\displaystyle=z$ (35) $$

and vice-versa.

In addition, if $U,z,u,v$ satisfy (34), and (35), we have:

$$ $\displaystyle\mathcal{S}(U)$ $\displaystyle=\mathcal{S}(u,v)$ $\displaystyle\mathcal{HS}(U,z)$ $\displaystyle=\mathcal{HS}(u,v)$ $$

###### Proof.

First, we pick $U,z$. We set $:v=Uz$. There exists a uniquely determined (up to scalar $\pm 1$) $u\in H(U)\cap v^{\perp}$ with $\|u\|=1$.
Since $U^{T}v=U^{T}Uz=z$,
we have that (34), and (35) hold.

Now consider $(u,v)\in\mathcal{B}^{d}$. We set $U=[u,v]$. (Note that the choice of $U$ is not unique, for example, we can also set $U=[-u,v]$ or $U=[u,-v]$.)
Let $z=U^{T}v=[0,1]^{T}$. We have that (34) and (35) hold.

Suppose $U,z,u,v$ satisfy (34) and (35). We have
$H(U,z)=H(u,v)$, and thus

$$ $\mathcal{S}(U)=\mathcal{S}(u,v).$ $$

Consider $s\in\mathcal{HS}(U,z)$, we decompose $s$ as
$s=s_{U}+s_{U^{\perp}}$, then (by Proposition B.2)

$$ $\mathcal{P}^{1}_{\mathcal{S}(U)}(s)=Uz=UU^{T}v=v.$ $$

The last equation holds from the fact $v\in H(U)$.
Since $u,v\in H(U)$, we have

$$ $\displaystyle s^{T}v$ $\displaystyle=s_{U}^{T}v=(UU^{T}s)^{T}v=\|U^{T}s\|v^{T}v\geq 0$ $\displaystyle s^{T}u$ $\displaystyle=s_{U}^{T}u=(UU^{T}s)^{T}u=\|U^{T}s\|v^{T}u=0$ $$

So, $s\in\mathcal{HS}(u,v)$ and thus

$$ $\mathcal{HS}(U,z)\subset\mathcal{HS}(u,v).$ $$

For the other direction, we pick $s\in\mathcal{HS}(u,v)$. We have $s^{T}u=0$.

Since $u\in H(U)$, we obtain:

$$ $\displaystyle 0=s^{T}u=s_{U}^{T}u+s_{U^{\perp}}^{T}u=s_{U}^{T}u$ $$

By combining this with the fact that $s_{U}\in H(U)=\text{span}(\{u,v\})$, we have $s_{U}=\alpha v$ for some $\alpha\in\mathbb{R}$.
In addition, from $s^{T}v\geq 0$, we have $\alpha\geq 0$.

Case 1: If $\alpha=0$, we have $s\in U^{\perp}$, thus $s\in\mathcal{HS}(U,z)$.

Case 2: If $\alpha>0$, we have:

$$ $P^{U}s=P^{U}s_{U}=\frac{U^{T}s}{\|U^{T}s\|}=\frac{U^{T}\alpha v}{\alpha}=U^{T} Uz=z.$ $$

Thus, $s\in\mathcal{HS}(U,z)$.
Therefore $\mathcal{HS}(u,v)\subset\mathcal{HS}(U,z)$.
∎

###### Remark B.4 .

When $d=3$, then $u=U[:,1]$, $v=U[:,2]$, and $n_{U}=\frac{1}{\|u\times v\|}u\times v$ (cross product, normalized).

Moreover, these $(d-2)$-dimensional hemispheres $\mathcal{HS}$ generalize the semicircles $\mathcal{SC}$ (see (24)) given for $d=3$.

###### Lemma B.3 .

###### Proof.

###### Remark B.4 .

### B.3 Spherical Radon transform and a Hemispherical transform

We call the integral transformation defined in (Groemer, , Section 2) as the “$(d-2)$-Hemispherical Transform on $\mathbb{S}^{d-1}$”: it is the $d-$dimensional extension of the Un-normalized Semicircle Transform given by (27) for $d=3$.
For $f\in L^{1}(\mathbb{S}^{d-1})$, we define

$$ $\mathcal{H}f(u,v):=\int_{\mathcal{HS}(u,v)}f(s)d\mathcal{HS}(u,v)$ $$

where $d\mathcal{HS}(u,v)$ is the surface area differential on $\mathcal{HS}(u,v)$. (We note that in the article Groemer (), the notation for the above transformation is $\mathcal{B}$.)

###### Remark B.5 .

In the literature, the so-called Hemispherical Transform integrates over a full $(d-1)$-dimensional hemisphere of $\mathbb{S}^{d-1}$, that is,

$$ $f(z):=\int_{\{x\in\mathbb{S}^{d-1}:\,x^{T}z\geq 0\}}f(x)dx\qquad\quad(z\in \mathbb{S}^{d-1}).$ (36) $$

See, e.g., Rubin (); Groemer ().

The Spherical Radon Transform in $\mathbb{S}^{d-1}$ is defined in Bonet et al. () as the linear and bounded integral operator $\mathcal{R}:L^{1}(\mathbb{S}^{d-1})\to L^{1}(V_{2}(\mathbb{R}^{d})\times
\mathbb{S}^{1})$ with the formula:

$$ $\mathcal{R}f(U,z):=\int_{\mathbb{S}^{d-1}}\delta(P^{U}s=z)f(s)ds.$ (37) $$

Thus, we can re-write it as

$$ $\mathcal{R}f(U,z)=\int_{\mathbb{S}^{d-1}}\delta(P^{U}s=z)f(s)ds=\int_{\mathcal {HS}(U,z)}f(s)d\mathcal{HS}(U,z)$ (38) $$

We define the following variant of such Spherical Radon Transform:

$$ $\displaystyle\widetilde{\mathcal{H}}(f)(U,z):=\int_{\mathbb{S}^{d-1}\setminus U ^{\perp}}\delta(P^{U}s=z)f(s)ds=\int_{\mathcal{HS}(U,z)\setminus U^{\perp}}f(s )d\mathcal{HS}(U,z).$ (39) $$

###### Proposition B.6 .

For $U\in V_{2}(\mathbb{R}^{d})$ and $z\in\mathbb{S}^{1}$, consider $(u,v)$ as in Lemma B.3.
Then

$$ $\displaystyle\mathcal{R}f(U,z)=\mathcal{H}f(u,v)=\widetilde{\mathcal{H}}f(U,z) ,\qquad\forall f\in L^{1}(\mathbb{S}^{d-1})$ $$

###### Proof.

By Lemma B.3, we have
$\mathcal{HS}(U,z)=\mathcal{HS}(u,v)$.
Then, by (38),

$$ ${\mathcal{R}}f(U,z)=\mathcal{H}f(u,v).$ $$

Thus, it remains to show $\widetilde{\mathcal{H}}f=\mathcal{R}f$.

For each $(U,z)$, $U^{\perp}$ is a $(d-2)$-dimensional subspace. Thus, $\mathbb{S}^{d-1}\cap U^{\perp}$ is a $(d-3)$-dimensional sub-sphere.

Since $\{s\in\mathbb{S}^{d-1}\setminus U^{\perp}\},U^{\perp}\cap\mathbb{S}^{d-1}$ are disjoint, we have:

$$ $\displaystyle\mathcal{HS}(U,z)\cap U^{\perp}$ $\displaystyle=\left[\{s\in\mathbb{S}^{d-1}\setminus U^{\perp}:P^{U}s=z\}\cup(U ^{\perp}\cap\mathbb{S}^{d-1}))\right]\cap U^{\perp}$ (40) $\displaystyle=U^{\perp}\cap\mathbb{S}^{d-1}$ $$

In addition, since $U^{\perp}\cap\mathbb{S}^{d-1}$ is a $({d-3})$-dimensional sub-sphere, while $\mathcal{HS}(U,z)=\mathcal{HS}(u,z)$ is a $(d-3)$-dimensional hemisphere,
we have

$$ $\displaystyle\mathcal{R}f(U,z)$ $\displaystyle=\widetilde{\mathcal{H}}f(U,z)+\int_{U^{\perp}\cap\mathbb{S}^{d-1 }}f(x)d\mathcal{HS}(U,z)$ $\displaystyle=\widetilde{\mathcal{H}}f(U,z)+0=\widetilde{\mathcal{H}}f(U,z)$ $$

and we have completed the proof.
∎

###### Corollary B.7 .

Given continuous function $f,g\in C(\mathbb{S}^{d-1})$, if $\mathcal{R}f=\mathcal{R}g$, then $f=g$.

###### Proof.

For every $(U,z)$ there exists $u,v\in\mathbb{S}^{d-1}$ such that, by Lemma B.6,

$$ $\mathcal{R}f(U,z)=\mathcal{H}f(u,v),\quad\mathcal{R}g(U,z)=\mathcal{H}g(u,v).$ $$

Thus, $\mathcal{H}f(u,v)=\mathcal{H}g(u,v)$.
From (Groemer, , Theorem 1), we have $f=g$.
∎

###### Remark B.8 .

In the article Bonet et al. () the authors relate the Spherical Radon Transform $\mathcal{R}$ with the Hemispherical Transform $\mathcal{F}$ given in (36).
Indeed,

$$ $\mathcal{R}f(U,z)=\mathcal{F}(f_{|_{\mathcal{S}(U)}})(z),$ $$

where $f_{|_{\mathcal{S}(U)}}$ denotes the restriction of $f\in L^{1}(\mathbb{S}^{d-1})$ to the $(d-2)$-dimensional “great circle” ${\mathcal{S}(U)}\simeq\mathbb{S}^{d-2}$, and the above operator $\mathcal{F}$ is the Hemispherical Transform in $L^{1}(\mathcal{S}^{d-2})$.
By doing so, the authors can use the characterization of $Ker(\mathcal{F})$ provided in Rubin () to prove their characterization of $Ker(\mathcal{R})$ (see (Bonet et al., , Proposition 4)). Also, in Bonet et al. () the authors claim that they leave for future works checking whether the set $Ker(\mathcal{R})$ is null or not. In this paper we show that it is more natural to relate $\mathcal{R}$ with the integral transformation given in (Groemer, , Section 2) that we have named as the “$(d-2)$-Hemispherical Transform on $\mathbb{S}^{d-1}$” and denoted by $\mathcal{H}$. Thus, (Groemer, , Theorem 1) we can show the injectivity of $\mathcal{R}$ for absolutely continuous measures on the sphere with continuous density functions. (Since $\mathbb{S}^{d-1}$ is compact, the space of continuous functions $C(\mathbb{S}^{d-1})$ is dense in $L^{1}(\mathbb{S}^{d-1})$ under closure with respect to the $L^{1}$-norm. However, we still can not guarantee that the injectivity property of $\mathcal{R}$ extends from $C(\mathbb{R}^{d-1})$ to $L^{1}(\mathbb{S}^{d-1})$, which represents the densities of absolutely continuous measures).

###### Remark B.5 .

###### Proposition B.6 .

###### Proof.

###### Corollary B.7 .

###### Proof.

###### Remark B.8 .

## Appendix C Spherical Radon transform for measures on 𝕊 d − 1 superscript 𝕊 𝑑 1 \mathbb{S}^{d-1} blackboard_S start_POSTSUPERSCRIPT italic_d - 1 end_POSTSUPERSCRIPT

Similarly to the classical Radon Transform in Euclidean spaces, the dual-operator of $\mathcal{R}$, called the back-projection operator and denoted as $\mathcal{R}^{*}:C_{0}(V_{2}(\mathbb{R}^{d})\times\mathbb{S}^{1})\to C(\mathbb{
S}^{d-1})$, is given in Bonet et al. () by the following formula:

$$ $\mathcal{R}^{*}(\psi)(x)=\int_{V_{2}(\mathbb{R}^{d})}\psi(U,P^{U}(x))d\sigma(U ),\qquad x\in\mathbb{S}^{d-1},\quad\psi\in C_{0}(V_{2}(\mathbb{R}^{d})\times \mathbb{S}^{1}).$ (41) $$

That is, $\mathcal{R}^{*}$ is such that the following identity holds:

$$ $\displaystyle\int_{\mathbb{S}^{1}\times V_{2}(\mathbb{R}^{d})}\mathcal{R}f(U,z )\psi(U,z)d\sigma(U)dz=\int_{\mathbb{S}^{d-1}}f(s)\,\mathcal{R}^{*}\psi(s)ds,$ (42) $$

for $f\in L^{1}(\mathbb{S}^{d-1}),\psi\in C_{0}(V_{2}(\mathbb{R}^{d})\times\mathbb{
S}^{1})$ (see (Bonet et al., , Proposition 2)).
Then, the Spherical Radon transform $\mathcal{R}$ can be extended from $L^{1}(\mathbb{S}^{d-1})$ to the space of measures supported on $\mathbb{S}^{d-1}$ by defining $\mathcal{R}\mu$ as a new measure on $V_{2}(\mathbb{R}^{d})\times\mathbb{S}^{1}$ by the “duality expression”:

$$ $\displaystyle\int_{\mathbb{S}^{1}\times V_{2}(\mathbb{R}^{d})}\psi(U,z)d \mathcal{R}\mu(U,z)=\int_{\mathbb{S}^{d-1}}\mathcal{R}^{*}\psi(s)d\mu(s), \qquad\forall\psi\in C_{0}(V_{2}(\mathbb{R}^{d})\times\mathbb{S}^{1}).$ (43) $$

###### Remark C.1 .

For each $s\in\mathbb{S}^{d-1}$, if $s\in U^{\perp}$, then $U[:,1],U[:,2]\in s^{\perp}$. Thus

$$ $\displaystyle\sigma(\{U\in V_{2}(\mathbb{R}^{d}):s\in U^{\perp}\})$ $\displaystyle=\sigma(\{U\in V_{2}(\mathbb{R}^{d}):U[:,1],U[:,2]\in s^{\perp}\})$ $\displaystyle=0$ $$

In Bonet et al. (), the authors mention the above that definition is considered a.e.. Via the above property, we have that (41) is well-defined for all $s\in\mathbb{S}^{d-1}$.

Now, consider probability measures $\mu\in\mathcal{P}(\mathbb{S}^{d-1})$. In analogy to the classical Radon Transform in Euclidean spaces, by the disintegration theorem in classic measure theory, there exists a
$\sigma$-almost everywhere uniquely determined family of $1$-dimensional probability measures $\{(\mathcal{R}\mu)^{U}\}_{U\in V_{2}(\mathbb{R}^{d})}$ on the unit circle $\mathbb{S}^{1}$ such that for any $\psi\in C_{0}(V_{2}(\mathbb{R}^{d})\times\mathbb{S}^{1})$ we have

$$ $\displaystyle\int_{V_{2}(\mathbb{R}^{d})\times\mathbb{S}^{1}}\psi(U,z)\,d \mathcal{R}\mu(U,z)=\int_{V_{2}(\mathbb{R}^{d})}\int_{\mathbb{S}^{1}}\psi(t, \theta)\,d(\mathcal{R}\mu)^{U}(z)\,d\sigma(U).$ (44) $$

In (Bonet et al., , Proposition 3) the authors prove that if $\mu_{f}$ is an absolutely continuous probability measure on the sphere $\mathbb{S}^{d-1}$, that is, it has a density $f\in L^{1}(\mathbb{S}^{d-1})$, then

$$ $(\mathcal{R}\mu_{f})^{U}=P^{U}_{\#}(\mu_{f})$ $$

The same computation can be extended to any measure $\mu$ supported on $\mathbb{S}^{d-1}$ to get $(\mathcal{R}\mu)^{U}=P^{U}_{\#}(\mu)$.

###### Remark C.2 .

In the particular case of $\mu_{f}$, that is, for functions $f\in L^{1}(\mathbb{S}^{d-1})$ we have that the measure $\mathcal{R}\mu_{f}$ has density $\mathcal{R}f$. This implies that the function $\mathcal{R}f(U,\cdot)$ (i.e., the function $z\mapsto\mathcal{R}f(U,z)$ defined for $z\in\mathbb{S}^{1}$) is the density of $(\mathcal{R}\mu_{f})^{U}$, and therefore $\mathcal{R}f(U,\cdot)$ is also the density of
$P^{U}_{\#}(\mu_{f})$.

###### Remark C.1 .

###### Remark C.2 .

### C.1 Linear Spherical Sliced Optimal Transport in terms of the Spherical Radon Transform

By using the previous expressions, given $\nu^{1},\nu^{2}\in\mathcal{P}(\mathbb{S}^{d-1})$ we can re-write $LSSOT_{2}(\cdot,\cdot)$ as follows:

$$ $\displaystyle\left(LSSOT_{2}(\nu_{1},\nu_{2})\right)^{2}$ $\displaystyle:=\int_{V_{2}(\mathbb{R}^{d})}\left(LCOT_{2}(P_{\#}^{U}\nu_{1},P_ {\#}^{U}\nu_{2})\right)^{2}\,d\sigma(U)$ $\displaystyle=\int_{V_{2}(\mathbb{R}^{d})}\left(LCOT_{2}((\mathcal{R}\nu_{1})^ {U},(\mathcal{R}\nu_{2})^{U})\right)^{2}\,d\sigma(U).$ (45) $$

## Appendix D Metric property of Linear Spherical Sliced Optimal Transport

###### Proposition D.1 .

$LSSOT_{2}(\cdot,\cdot)$ is a well-defined pseudo-metric in $\mathcal{P}(\mathbb{S}^{d-1})$, that is, $LSSOT_{2}(\nu,\nu)=0$ for all $\nu\in\mathcal{P}(\mathbb{S}^{1})$, it is non-negative, symmetric, and satisfies the triangle inequality.

###### Proof.

Let $\nu_{1},\nu_{2},\nu_{3}\in\mathcal{P}(\mathbb{S}^{d-1})$.

It is straightforward to verify positivity
(i.e., $LSSOT_{2}(\nu_{1},\nu_{2})\geq 0$), and symmetry (i.e., $LSSOT_{2}(\nu_{1},\nu_{2})=LSSOT_{2}(\nu_{2},\nu_{1})$).

If $\nu_{1}=\nu_{2}$, for each $U\in V_{2}(\mathbb{R}^{d})$, we have $P_{\#}^{U}\nu_{1}=P_{\#}^{U}\nu_{2}$, thus
$\hat{\nu}_{1}^{S}(\cdot,U)=\hat{\nu}_{2}^{S}(\cdot,U)$.
It implies that
$\|\hat{\nu}_{1}^{S}(\cdot,U)-\hat{\nu}_{2}^{S}(\cdot,U)\|^{2}_{L^{2}(\mathbb{S
}^{d-1})}=0$, and so
$LSSOT_{2}(\nu_{1},\nu_{2})=0$.

Finally, we verify the triangle inequality.
Let $U\in V_{2}(\mathbb{R}^{d})$.
By (Martin et al., , Theorem 3.6), we have $LCOT_{2}(\cdot,\cdot)$ defines a metric. Thus, the following triangle inequality holds:

$$ $LCOT_{2}(P_{\#}^{U}\nu_{1},P_{\#}^{U}\nu_{2})\leq LCOT_{2}(P_{\#}^{U}\nu_{1},P _{\#}^{U}\nu_{3})+LCOT_{2}(P_{\#}^{U}\nu_{2},P_{\#}^{U}\nu_{3}).$ $$

Thus,

$$ $\displaystyle LSSOT_{2}(\nu_{1},\nu_{2})$ $\displaystyle=\left(\int_{V_{2}(\mathbb{R}^{d})}(LCOT_{2}(P_{\#}^{U}\nu_{1},P_ {\#}^{U}\nu_{2}))d\sigma(U)\right)^{1/2}$ $\displaystyle\leq\left(\int_{V_{2}(\mathbb{R}^{d})}(LCOT_{2}(P_{\#}^{U}\nu_{1} ,P_{\#}^{U}\nu_{3})+LCOT_{2}(P_{\#}^{U}\nu_{2},P_{\#}^{U}\nu_{3}))d\sigma(U) \right)^{1/2}$ $\displaystyle\leq\left(\int_{V_{2}(\mathbb{R}^{d})}(LCOT_{2}(P_{\#}^{U}\nu_{1} ,P_{\#}^{U}\nu_{3}))d\sigma(U)\right)^{1/2}+\left(\int_{V_{2}(\mathbb{R}^{d})} (LCOT_{2}(P_{\#}^{U}\nu_{2},P_{\#}^{U}\nu_{3}))d\sigma(U)\right)^{1/2}$ (46) $\displaystyle=LSSOT_{2}(\nu_{1},\nu_{3})+LSSOT_{2}(\nu_{2},\nu_{3})$ $$

where (46) holds by Minkowski inequality in space $L^{2}(V_{2}(\mathbb{R}^{d}),d\sigma)$.
∎

###### Theorem D.2 (Theorem 3.3 ) .

$LSSOT_{2}(\cdot,\cdot)$ is a well-defined metric in the set of absolutely continuous measures in $\mathcal{P}(\mathbb{S}^{d-1})$ with continuous density function, that is, it is pseudo-metric satisfying also the identity of indiscernibles.

###### Proof.

From Proposition D.1, we only have to prove that if $LSSOT_{2}(\nu_{1},\nu_{2})=0$, then $\nu_{1}=\nu_{2}$. To do so, we will the expression of the $LSSOT_{2}(\cdot,\cdot)$ in terms of a Spherical Radon Transform given in given in (C.1). Let $\nu_{1},\nu_{2}\in\mathcal{P}(\mathbb{S}^{d-1})$ such that

$$ $0=\left(LSSOT_{2}(\nu_{1},\nu_{2})\right)^{2}=\int_{V_{2}(\mathbb{R}^{d})} \left(LCOT_{2}((\mathcal{R}\nu_{1})^{U},(\mathcal{R}\nu_{2})^{U}\right)^{2}\,d \sigma(U).$ $$

Thus,

$$ $LCOT_{2}((\mathcal{R}\nu_{1})^{U},(\mathcal{R}\nu_{2})^{U})=0\qquad\sigma(U)-a .s.$ $$

Since $LCOT_{2}(\cdot,\cdot)$ defines a metric in $\mathcal{P}(\mathbb{S}^{1})$ (see again (Martin et al., , Theorem 3.6)), we have that

$$ $(\mathcal{R}\nu_{1})^{U}=(\mathcal{R}\nu_{2})^{U}\qquad\sigma(U)-a.s.$ $$

As the disintegration is a family $\sigma$-almost everywhere uniquely determined, we have that

$$ $\mathcal{R}\nu_{1}=\mathcal{R}\nu_{2}.$ $$

Finally, if $\nu_{1},\nu_{2}$ have continuous density functions, by injectivity of the operator $\mathcal{R}$ (Corollary B.7), we have that $\nu_{1}=\nu_{2}$.
∎

As a byproduct of our deductions, we obtain that the SSW in Bonet et al. () and defined in (8), is also a metric in $\mathcal{P}(\mathbb{S}^{d-1})$:

###### Theorem D.3 .

The Spherical Sliced Wasserstein discrepancy $SSW_{p}(\cdot,\cdot)$, for $p\geq 1$, defines a metric in the set of absolutely continuous measures in $\mathcal{P}(\mathbb{S}^{d-1})$ with continuous density function.

###### Proof.

First, notice that since the sphere is compact, there is no need to ask that $\nu\in\mathcal{P}(\mathbb{S}^{d-1})$ has to have finite $p$-moment, as it automatically has such a property. In (Bonet et al., , Proposition 5), the authors show that $SSW_{p}(\cdot,\cdot)$ is a pseudo-metric. Finally,
by using the expression

$$ $\displaystyle(SSW_{p}(\mu,\nu))^{p}$ $\displaystyle:=\int_{V_{2}(\mathbb{R}^{d})}\left(COT_{p}(P_{\#}^{U}\mu,P_{\#}^ {U}\nu)\right)^{p}\,d\sigma(U)$ (47) $\displaystyle=\int_{V_{2}(\mathbb{R}^{d})}\left(COT_{p}((\mathcal{R}\mu)^{U},( \mathcal{R}\nu)^{U})\right)^{p}\,d\sigma(U)$ (48) $$

and repeating the arguments in the proof of Theorem 3.3 with the Wasserstein distance $COTp(\cdot,\cdot)$ in the unit circle in place of the LCOT distance, we get that $SSW_{p}(\cdot,\cdot)$ satisfies the identity of indiscernibles when restricted to probability measures with continuous densities.
∎

###### Proposition D.1 .

###### Proof.

###### Theorem D.2 (Theorem 3.3 ) .

###### Proof.

###### Theorem D.3 .

###### Proof.

## Appendix E Runtime Analysis with respect to Number of Slices

Here we provide further running time comparisons of the slice-based methods for pairwise distance calculations under different numbers of slices $L=500,1000,5000$. All methods are run on NVIDIA RTX A6000 GPU.

Figure: Figure 6: Running time comparison under varying number of slices. $L$ denotes the number of slices used in all the slice-based methods in each row. $N$ is the number of samples in each distribution for each column.
Refer to caption: x5.png

## Appendix F Implementation Details of Cortical Surface Registration

### F.1 Icospheres and Resampling

S3Reg relies on Spherical-Unet (Zhao et al., ) as the backbone, which capitalizes on the regular structure of icospheres to perform convolution and pooling on spheres. Spherical surfaces are resampled on icospheres in order to be input to the S3Reg network.

Icospheres are generated from a regular convex icosahedron with 12 vertices, then by recursively adding a new vertex to the center of each edge in each triangle. The number of vertices on the icospheres is increased from 12 to 42 (level 1), 162 (level 2), 642 (level 3), 2562 (level 4), 10,242 (level 5), 40,962 (level 6), 163,842 (level 7), etc. In this work, we focus on the level-6 icosphere with 40,962 vertices on the cortical surfaces.

For resampling in the preprocessing step, we employ the simple but fast Nearest Neighbor method and implement it using sklearn.neighbors.KDTree function from Scikit Learn (Pedregosa et al., ). The feature value on each output point is represented as a weighted sum of the feature values on the nearest three vertices.

### F.2 S3Reg Framework and Loss Components

In a nutshell, the S3Reg framework receives a concatenation of fixed and moving spherical surface maps defined on icospheres as input, processes them through an encoder-decoder architecture, and finally outputs the velocity field represented by tangent vectors on $\mathbb{S}^{2}$. In particular, as the convolution kernel in Spherical U-Net could cause distortions around the poles (Zhao et al., ), 3 complementary orthogonal Spherical U-Nets are combined to predict the deformation field so that the polar areas are also regularly processed at least two times. As the 3 predicted velocity fields should be orthogonal accordingly, a loss term of inconsistency is added to the total loss. We refer the readers to Zhao et al. () (Section III-A) for details of the consistency loss.

As the similarity losses are thoroughly discussed in the main text, we will now describe other loss components in more detail. Let $F$ denote the fixed image and $M$ be the moving image. Suppose $p_{F}$, $p_{M}$ are the corresponding parcellation maps of $F$ and $M$, respectively. Let $\phi$ denote the velocity field predicted by the network.

Dice Score and Dice Loss
Parcellation maps can be integrated into the training process to enforce biological validity. The Dice Score is a measure of similarity/overlap between two parcellations. Assume both $p_{F}$ and $p_{M}$ have $K$ regions of interest (ROIs), then the Dice score is defined as

$$ $\displaystyle\text{Dice}(p_{F},p_{M}\circ\phi)=\frac{1}{K}\sum_{k=1}^{K}2\frac {|p_{F}^{k}\cap(p_{M}^{k}\circ\phi)|}{|p_{F}^{k}|+|p_{M}^{k}\circ\phi|},$ (49) $$

where $p_{M}^{k}\circ\phi$ is the warped moving image by $\phi$, the superscript $k$ denotes the $k$-th ROI. Since the Dice Score is between 0 0 and $1$, the Dice Loss can be defined as $L_{\text{Dice}}(p_{F},p_{M}\circ\phi)=1-\text{Dice}(p_{F},p_{M}\circ\phi)$.

Smoothness Loss
The local spatial non-smoothness of the output velocity field $\phi$ is penalized by $L_{\text{smooth}}(\phi)=\sum_{v}|\nabla_{\mathbb{S}^{2}}C(\phi(v))|$, where the sum is over all the moving vertices $v$, $C$ is the 1-ring convolution operator (Zhao et al., ), and $\nabla_{\mathbb{S}^{2}}$ denotes the spherical gradients.

## Appendix G More visualizations of Cortical Surface Registration

Below, we provide additional visualizations of our cortical surface registration comparison. The figures are presented in landscape mode to enhance visual inspection.

Figure: Figure 7: Qualitative registration results from Subject A00060516 in the NKI dataset.
Refer to caption: x6.png

Figure: Figure 8: Qualitative registration results from Subject A00076798 in the NKI dataset.
Refer to caption: x7.png

Figure: Figure 9: Qualitative registration results from Subject 0356 in the ADNI dataset.
Refer to caption: x8.png

Figure: Figure 10: Qualitative registration from Subject 6650 in the ADNI dataset.
Refer to caption: x9.png

Figure: Figure 11: Qualitative registration results from Subject 7060 in the ADNI dataset.
Refer to caption: x10.png

## Appendix H 3D Shape Clustering and Retrieval

To further validate the effectiveness of LSSOT in evaluating the similarity of spherical distributions, we conduct additional experiments using a variety of genus-0 3D shapes. Our methodology involves a two-step process: 1) mapping the 3D shapes onto spherical surfaces; and 2) performing clustering and nearest neighbor retrieval tasks on these spherical representations. For the spherical mapping process, we employ two distinct techniques, Inflation and Projection (IAP) and Deep Unsupervised Learning using PointNet++ Qi et al. ().

### H.1 Dataset

We conduct experiments on mesh models from the Watertight Track of SHREC 2007 dataset (Giorgi et al., ). This dataset contains surfaces belonging to 19 categories, among them, we select 4 genus-0 closed surface models, human, octopus, teddy bear, and Fish. Moreover, as explored in the previous experiment, LSSOT is sensitive to rotations, therefore we remove the meshes with different orientations than the rest.

### H.2 Spherical Mapping

Inflation and Projection.
A direct projection from a genus-0 closed surface to the sphere by fixing an inside point source and normalizing the vertices may cause undesirable topologically incorrect foldings. One workaround is iteratively smoothing and inflating the surface before the final projection. This inflation and projection is a common spherical mapping strategy used in neuroimaging data processing tools such as FreeSurfer (Fischl, ), Connectome Workbench (Marcus et al., ).

Deep Unsupervised Learning using PointNet++.
To accelerate the forward pass of spherical mapping, we capitalize on deep unsupervised learning and train a PointNet++ network to map a cortical surface on a sphere with minimal distortion. The trained network then serves as a deep learning-based spherical mapping tool.

Loss function.
Zhao et al. () proposed an unsupervised approach of mapping cortical surfaces to a sphere using Spherical U-Net Zhao et al. () to learn the spherical diffeomorphic deformation field for minimizing the metric (distance), area, or angle distortions. We follow their definitions of metric ($J_{d}$) and area ($J_{a}$) distortions as local relative differences, and conformal/angle distortion ($J_{c}$) after market share normalization, but adapt them to general meshes instead of restricting to icosahedron-reparameterized surfaces. For a pair of original input surface $R$ and the predicted spherical mesh $S$, the loss function has three components:

$$ $\displaystyle J_{d}(R,S)=\min_{r\in\mathbb{R}^{+}}\frac{1}{N_{e}}\sum_{i=1}^{N _{e}}\frac{|rd_{i}^{S}-d_{i}^{R}|}{d_{i}^{R}}$ (50) $\displaystyle J_{a}(R,S)=\min_{r\in\mathbb{R}^{+}}\frac{1}{N_{f}}\sum_{j=1}^{N _{f}}\frac{|ra_{j}^{S}-a_{j}^{R}|}{a_{j}^{R}}$ (51) $\displaystyle J_{c}(R,S)=\frac{1}{3N_{f}}\sum_{k=1}^{N_{f}}\frac{1}{\pi}(| \alpha^{R}_{k}-\alpha^{S}_{k}|+|\beta^{R}_{k}-\beta^{S}_{k}|+|\theta^{R}_{k}- \theta^{S}_{k}|)$ (52) $$

where $N_{e}$ is the total number of edges, $N_{f}$ is the total number of faces. $d_{i}^{R}$ and $d_{i}^{S}$ are the length of $i$-th edge on original surface and predicted sphere. Similarly, $a_{j}^{R}$ and $a_{j}^{S}$ are the area of $j$-th face on original surface and predicted sphere; $(\alpha^{R}_{k},\beta_{k}^{R},\theta^{R}_{k})$ and $(\alpha^{S}_{k},\beta_{k}^{S},\theta^{S}_{k})$ are the inner angles of the $k$-th face on original surface and predicted sphere.

Figure: Figure 12: Manifold learning of clustering results (top) and nearest neighbor retrieval results (bottom) using Inflation and Projection (left) and PointNet++ (right) as spherical mappings.
Refer to caption: x11.png

### H.3 Results

Figure 12 illustrates the effectiveness of the LSSOT metric under both spherical mapping techniques. The results demonstrate that LSSOT successfully identifies inherent similarities and dissimilarities among the spherical surface representations, regardless of the mapping method employed. This capability translates into good performance in both clustering and nearest neighbor retrieval tasks. The consistency of LSSOT’s effectiveness under different mapping approaches underscores its robustness as a metric for analyzing spherical distributions and comparing surface geometries.

## Appendix I Gradient flow with respect to LSSOT in the spherical space

### I.1 background: introduction to Stiefel manifold

In this section, we briefly introduce the gradient flow and related concepts in spherical space.
First, we introduce the Riemannian structure of the Stiefel manifold, which is a general case of spherical space.
Next, we specify the related concepts in the spherical space.

Given positive integers $1\leq k\leq d$, the Stiefel manifold is given by

$$ $\mathbb{V}_{d,k}:=\mathbb{V}_{k}(\mathbb{R}^{d})=\{V\in\mathbb{R}^{d\times k}: V^{\top}V=I_{k}\}.$ $$

It is a $\left(dk-(k+\binom{k}{2})\right)$-dimensional manifold.

Given $X\in\mathbb{V}_{d,k}$, the tangent space is given by:

$$ $\displaystyle\mathcal{T}_{X}(\mathbb{V}_{d,k}):=\{Z\in\mathbb{R}^{d\times k}: ZX+XZ=0\}.$ (53) $$

Given matrix $G\in\mathbb{R}^{d\times k}$, the projection of $G$ onto $\mathcal{T}_{X}(\mathbb{V}_{d,k})$ is given by,

$$ $\displaystyle P_{T_{X}(\mathbb{V}_{d,k})}(G):=G-X\text{sym}(G^{\top}X),$ (54) $$

where $\text{sym}(G^{\top}X)=\frac{1}{2}(G^{\top}X+X^{\top}G)$ is the symmetric matrix obtained from $G^{\top}X$.

Given a differentiable function $F:V_{d,k}\to\mathbb{R}$, the following steps can describe the gradient of $F$ in the manifold. In particular, choose $X\in\mathbb{V}_{d,k}$:

- •
Step 1: compute Euclidean gradient $\nabla_{E}(F(X))$.
- •
Step 2: the Stiefel manifold gradient is defined as:
$\displaystyle\nabla_{S}F(X):=P_{\mathcal{T}_{X}(\mathbb{V}_{d,k})}(\nabla_{E}(
F(X))).$
(55)

We consider the gradient flow problem:

$$ $\partial_{t}X_{t}=-\nabla_{S}F(X_{t}),\forall X\in\mathbb{V}_{d,k}.$ $$

Let $\eta$ denote the step size (i.e., the learning rate for the machine learning community), the discrete gradient descent step of the above problem becomes:

$$ $\displaystyle X\mapsto X-\eta\nabla_{S}F(X)=X-\eta(\nabla_{E}F(X)-X\text{sym}( \nabla_{E}F(X)^{\perp}X)).$ (56) $$

### I.2 Gradient flow in spherical space

In spherical space, we have $k=1$ (in this case, we use convention $\binom{1}{2}=0$). Thus, the tangent space (53) becomes:

$$ $\displaystyle\mathcal{T}_{x}(\mathbb{S}^{d-1}):=\text{span}(x)^{\perp}:=\{z\in \mathbb{R}^{d}:z^{\top}x=0\},$ (57) $$

where $\text{span}(x)$ is the 1D line spaned by $x$. The corresponding projection (54) becomes

$$ $\displaystyle\mathbb{R}^{d}\ni g\mapsto P_{\mathcal{T}_{x}(\mathbb{V}_{d,k})}( g):=g-x(g^{\top}x)\in\mathcal{T}_{x}(\mathbb{S}^{d-1}).$ (58) $$

Note, the above projection is essentially the Euclidean projection of $x$ into the hyperplane $x^{\perp}$.

Given $X=[x_{1},\ldots,x_{n}]\in\mathbb{R}^{d\times n},Y=[y_{1},\ldots,y_{m}]\in
\mathbb{R}^{d\times m}$, where $Y$ is fixed.
Let $\nu_{1}=\sum_{i=1}^{n}p_{i}\delta_{x_{i}},\nu_{2}=\sum_{j=1}^{m}q_{j}\delta_{y
_{j}}$. We consider the following function
:

$$ $\displaystyle C(X)$ $\displaystyle:=LSSOT(\mu,\nu)$ $\displaystyle\approx\frac{1}{T}\sum_{t=1}^{T}LCOT((\theta_{t})_{\#}\mu,(\theta _{t})_{\#}\nu).$ (59) $$

where each $\theta_{t}\in\mathbb{S}^{d-1}$, (59) is the Monte Carlo approximation.

We select one $\theta$ and consider the 1-dimensional LCOT problem

$$ $\displaystyle LCOT((\theta)_{\#}\nu_{1},(\theta)_{\#}\nu_{2})$ $\displaystyle=\int_{0}^{1}h\left(\left|F_{\theta_{\#}\nu_{1}}^{-1}(s-\mathbb{E }[\theta_{\#}\nu_{1}]+\frac{1}{2})-|F_{\theta_{\#}\nu_{2}}^{-1}(s-\mathbb{E}[ \theta_{\#}\nu_{2}]+\frac{1}{2})\right|^{2}\right)ds$ $$

where $h(r)=\min(r,1-r)$.

Note, the above quantity implicitly defines a transportation plan between $\mu$ and $\nu$. In particular, suppose $\gamma_{\theta}\in\mathbb{R}^{n\times m}$, with

$$ $\displaystyle\gamma_{\theta}[i,j]=|\{s:F_{\theta_{\#}\nu_{1}}^{-1}(s-\mathbb{E }[\theta_{\#}\nu_{1}])\equiv x_{i},F_{\theta_{\#}\nu_{1}}^{-1}(s-\mathbb{E}[ \theta_{\#}\nu_{2}])\equiv y_{j}\}|,$ (60) $$

where $|\cdot|$ is the uniform measure (Lebesgue measure) in $[0,1)$, $a\equiv_{1}b$ means $a\mod 1=b$.
It is straightforward to verify $\gamma_{\theta}\in\Gamma(\theta_{\#}\nu_{1},\theta_{\#}\nu_{2})$, and thus is a transportation plan between $\theta_{\#}\nu_{1},\theta_{\#}\nu_{2}$.

Therefore, we can explicitly write (59) in term of $X=[x_{1},\ldots x_{n}]$:

$$ $\displaystyle(\ref{eq:C(X)})=\frac{1}{T}\sum_{t=1}^{T}\sum_{i,j}h(\theta_{t}^{ \top}x_{i}-\theta_{t}^{\top}y_{j})\gamma_{\theta_{t}}[i,j].$ (61) $$

where each $\gamma_{\theta_{t}}$ is defined via (60).
For each $x_{i}$, the Euclidean gradient is given by

$$ $\displaystyle\nabla_{E}C(x_{i})$ $\displaystyle=\frac{1}{T}\sum_{t=1}^{T}\theta_{t}\left(\sum_{j:|\theta_{t}^{ \top}x_{i}-\theta_{t}^{\top}y_{j}|\leq 1/2}2(\theta_{t}^{\top}x_{i}-\theta_{t} ^{\top}y_{j})\gamma_{\theta_{t}}[i,j]\right.$ $\displaystyle\qquad-\left.\sum_{j:|\theta_{t}^{\top}x_{i}-\theta_{t}^{\top}y_{ j}|>1/2}2(\theta_{t}^{\top}x_{i}-\theta_{t}^{\top}y_{j})\gamma_{\theta_{t}}[i, j]\right).$ (62) $$

Based on this, we derive the gradient descent step for the following gradient flow problem.

$$ $-\nabla_{S}C(X)=-V(X)$ $$

That is,

$$ $\displaystyle X\mapsto X-\eta\nabla_{S}(C(X))=X-(X-X(\nabla_{E}C(X)^{\top}X))$ (63) $$

where $\nabla_{E}(C(X))[:,i]$ is given by (62) for each $i$.

### I.3 Implementation of Gradient Flow

Riemannian Gradient Descent. Let $\mu_{\mathcal{X}}=\frac{1}{N}\sum_{n=1}^{N}\delta_{x_{n}}$ be the source distribution with mass uniformly distributed at $\mathcal{X}=\{x_{n}\}_{n=1}^{N}\subset S^{d-1}$, and $\nu_{\mathcal{Y}}=\frac{1}{M}\sum_{m=1}^{M}\delta_{y_{m}}$ be the target distribution with $\mathcal{Y}=\{y_{m}\}_{m=1}^{M}\subset S^{d-1}$. For each $x\in\mathcal{X}$, let $x^{(k)}$ be the location of $x$ at $k$-th iteration, then at $(k+1)$-th iteration, the update should be the following (Absil et al., ; Boumal et al., ):

$$ $x^{(k+1)}=\exp_{x^{(k)}}(-\gamma_{k}\nabla_{S^{d-1}}\mathcal{F}(x^{(k)})),$ (64) $$

where $\gamma_{k}$ is the learning rate at $k$-th iteration, the exponential map $\exp_{x}(v)=x\cos{\|v\|}+\frac{v}{\|v\|}\sin{\|v\|}$, and $\nabla_{S^{d-1}}\mathcal{F}(x^{(k)})$ is the projection of the Euclidean gradient $\nabla_{\mathbb{R}^{d}}\mathcal{F}(x^{(k)})$ onto the tangent space at $x^{(k)}$. We implement this gradient descent method using the exponential function and projection function in the Geoopt library (Kochurov et al., ) (^5^55https://github.com/geoopt/geoopt).

Spherical Coordinate Gradient Descent. A more straightforward gradient descent method is to work in the spherical coordinate system. By the following correspondences (in $\mathbb{S}^{2}$ for example) between $(x,y,z)\in\mathbb{S}^{2}$ and $(\theta,\phi)$

$$ $\displaystyle x=\cos{\theta}\sin{\phi}$ $\displaystyle y=\sin{\theta}\sin{\phi}$ $\displaystyle z=\cos{\phi}$ $$

and

$$ $\displaystyle\theta=\arctan(\frac{y}{x})$ $\displaystyle\phi=\arccos{z}$ $$

we can calculate the gradient of the loss function with respect to $(\theta,\phi)$. Note here $\arctan$ must be suitably defined according to the quadrant.

### I.4 Additional Point Cloud Interpolation Visualizations

Figure: Figure 13: Gradient flow interpolations between two tables.
Refer to caption: extracted/5988566/35_53.png

Figure: Figure 14: Gradient flow interpolations between two toilets.
Refer to caption: extracted/5988566/92_99.png

Figure: Figure 15: Gradient flow interpolations from a range hood to a plant.
Refer to caption: extracted/5988566/37_82.png

Figure: Figure 16: Gradient flow interpolations from a TV stand to a bookshelf.
Refer to caption: extracted/5988566/74_68.png

Figure: Figure 17: Gradient flow interpolations from a table to a piano.
Refer to caption: extracted/5988566/table_piano.png

Figure: Figure 18: Gradient flow interpolations from a guitar to an airplane.
Refer to caption: extracted/5988566/guitar_airplane.png

### I.5 Gradient Flow Mollweide View Visualizations

We provide more visualizations for the gradient flow of spherical distributions with respect to LSSOT metric, SSW, and spherical optimal transport distance (spherical OT). At first, a source distribution and a target distribution are generated, both of which are discrete distributions with uniform mass at all points. Then we perform gradient descent on the mass locations of the source distribution with respect to the above three objectives, driving the source toward the target distribution by minimizing these losses. We employ two methods of gradient descent on the sphere.
See Figure 19 to Figure 30 for the Mollweide views of the gradient flows between two different pairs of source-target distributions.

Figure: Figure 19: LSSOT gradient flow using Riemannian gradient descent.
Refer to caption: x12.png

Figure: Figure 20: SSW gradient flow using Riemannian gradient descent.
Refer to caption: x13.png

Figure: Figure 21: Spherical OT gradient flow using Riemannian gradient descent.
Refer to caption: x14.png

Figure: Figure 22: LSSOT gradient flow using spherical coordinate gradient descent.
Refer to caption: x15.png

Figure: Figure 23: SSW gradient flow using spherical coordinate gradient descent.
Refer to caption: x16.png

Figure: Figure 24: Spherical OT gradient flow using spherical coordinate gradient descent.
Refer to caption: x17.png

Figure: Figure 25: LSSOT gradient flow using Riemannian gradient descent.
Refer to caption: x18.png

Figure: Figure 26: SSW gradient flow using Riemannian gradient descent.
Refer to caption: x19.png

Figure: Figure 27: Spherical OT gradient flow using Riemannian gradient descent.
Refer to caption: x20.png

Figure: Figure 28: LSSOT gradient flow using spherical coordinate gradient descent.
Refer to caption: x21.png

Figure: Figure 29: SSW gradient flow using spherical coordinate gradient descent.
Refer to caption: x22.png

Figure: Figure 30: Spherical OT gradient flow using spherical coordinate gradient descent.
Refer to caption: x23.png

## Peer Reviews

confidence: 3
contribution: 3
initial_rating: 6
presentation: 3
questions: see weaknesses section above.
soundness: 3
strengths: The paper has several strengths including the clarity of the writing, problem formulation, and valuable contribution to problem of spherical data analysis. I appreciated the authors placing their work in the context of past work in sliced optimal transport and linear optimal transport. They elegantly combine aspects from both to develop the spherical OT framework. There is theoretical justification for the use of the LSSOT as a distance metric, and an efficient implementation for speed. The authors also demonstrate the utility of their method on several experiments, notably cortical surface registration, where working in a spherical domain is essential.
summary: The authors propose the linear spherical sliced optimal transport framework (LSSOT) as an efficient method for comparing probability distributions on spherical manifolds. The authors build on several valuable past works that aim to linearize or speed-up the optimal transport problem by restricting to specific domains. In their case, the spherical domain has extensive applications in computer vision and especially in neuroscience for problems of cortical surface registration. They demonstrate results on several applications including point cloud interpolation and cortical surface registration.
weaknesses: I have several questions about the experimental evaluation of this work, specifically as it pertains to cortical surface registration. I appreciate the use of multiple datasets and outlining the problem of surface-to-atlas registration. The authors use extensive evaluation metrics. However, I have some concerns with the experimental setup and results.

The main purpose of these experiments was to quantify the improvement in performance and speed in using their metric over conventional ones. The authors state that they also use a Dice loss and regularization term. Did the authors perform new hyperparameter tuning on these losses when switching the image similarity metric? Failure to do so may bias results.

Why was only the sulcal depth used as a registration feature? Curvature is typically used and also offers a lot of information for registration.

The results on registration are puzzling. The authors' method performs considerably worse on SWD than the MSE loss for the NKI (Left hemisphere) and in MAE for NKI (Right hemisphere). Can the authors comment why? 

In general, registration is an ill-posed problem, and one can modify the hyperparameters appropriately to either obtain very close matchings or well behaved fields. Can the authors discuss this and describe how this tradeoff was achieved. 

How do these results compare to classical freesurfer registration? Which although is slow, represents a SOA in performance.

---

confidence: 1
contribution: 3
initial_rating: 6
presentation: 3
questions: The suggestions are mentioned in Weaknesses(2,3,4).
soundness: 3
strengths: 1.	According to time complexity and runtime analysis and the results in Table 1 in this paper, the proposed method achieves significant computational efficiency while preserving the spherical geometry, making it suitable for large-scale applications.

2.	The authors prove that the proposed LSSOT satisfies metric properties, strengthening the approach's theoretical validity.

3.	This work performed extensive experimental validation. The effectiveness and efficiency of LSSOT are evaluated across multiple applications, including cortical surface registration, 3D point cloud analysis, and shape embedding.

4.	The cortical surface registration results indicate that LSSOT achieves high Dice scores and minimizes area/edge distortions in the deformation field. This may provide biologically meaningful insights into neuroimaging applications.
summary: This paper introduces Linear Spherical Sliced Optimal Transport (LSSOT), a framework for embedding spherical probability distributions into L2 space while maintaining spherical geometry. Extending previous work on circular distributions, LSSOT offers a computationally efficient metric for spherical data comparison, proving to be both a valid metric and faster than baseline approaches. The effectiveness of LSSOT is demonstrated through applications in cortical surface registration, 3D point cloud analysis, and shape embedding.
weaknesses: 1.	From a technical point of view, the extension of the proposed method to the existing LOCT is limited. Still, it should be noted that it is useful in spherical signal analysis applications.

2.	LSSOT is compared to several baseline metrics, but the paper lacks a more detailed statistical analysis to quantify performance improvements (e.g., significance testing on Dice scores and area distortions). A statistical evaluation would help substantiate LSSOT’s claimed advantages over other metrics.

3.	Although LSSOT demonstrates effective interpolation performance, there is little information on its computational complexity compared to SSW and Spherical OT in point cloud analysis. A complexity analysis of interpolation tasks would better illustrate LSSOT's efficiency benefits.

4.	Hyperparameter settings, such as λ, the number of slices L, and reference size M in LSSOT, are briefly mentioned but not systematically analyzed. A discussion on the sensitivity of these parameters and their impact on performance would be beneficial, particularly for users looking to optimize LSSOT for specific applications.

---

confidence: 4
contribution: 4
initial_rating: 10
presentation: 4
questions: This question pertains to the original circular optimal transport construction, and the authors do not have to answer it here, as they do something different in this paper. However, out of interest, why is the COT distance defined in terms of c(x, y) uses min{|x−y|,1 −|x−y|}? Can some other distance be used here? For e.g. min(arccos(<x, y>)). Is that the main difference between this earlier work and  the spherical sliced Wasserstein (SSW) approach? 


The computational efficiency of this approach also naturally follows from the idea of slicing (projections on the great circle). Further efficiency is achieved by replacing the circular optimal transport distance by the linear version. However, can the authors provide slightly  more details on why the complexity of spherical sliced Wasserstein better than sliced Wasserstein? The projections on spheres should be more costly than lines. 

Can the authors comment on why LSSOT has better accuracy than SSW? Should SSW in fact perform as well as or better than LSSOT? 
For SSW, the COT distance measures the discrepancy between the push-forward maps of the spherical measures, while LSSOT uses the LCOT distance. So there is yet one more embedding ( transport a reference measure to a target measure) performed. Does this result in loss of accuracy? Is this dataset dependent? Can the authors provide more detailed analysis or any ablation studies to explain this result? 

Can the authors expand more upon the injectivity of R? This is related to the results in theorems D.2 and D.3 and is inspired by the theorem in Bonet et al. 2024.  For e.g. the map R: L^1(S^{d−1}) → L^1(V_2(R^d) ×S^1) is not injective everywhere, although it is true for absolutely continuous measures on the sphere. What are the implications when the operator R is applied in practice? This could be discussed more for the benefit of the reader.
soundness: 4
strengths: This paper presents an elegant construction of the spherical sliced embedding in a linear sense. 

The main strength of the paper is the clever extension of the sliced optimal transport and linear circular optimal transport to the case of spherical distributions. This is novel.

There exist other methods for computing distances between spherical distributions, for e.g. spherical sliced-Wasserstein, Stereographic sliced-Wasserstein which also allows to efficiently compute distances. Other works have achieved further computational efficiency due to the linear distance between the circular transport embeddings. In this paper, the authors first obtain an embedding of the spherical sliced measures P(S^{d-1}) and then define a distance between the measures as the linear circular optimal transport distance between the pushforward maps of the spherical measures. While this paper borrows and combines ideas from two earlier works, the novelty of this approach is that in this instance it is achieved for embeddings of spherical distributions. 

The LSSOT distance is a pseudo-metric on P(S^{d−1}) as proven in an earlier work (spherical sliced Wasserstein), since they use the same notion of the distance. In this work the authors show that it is a well-defined metric for probability measures. Also see weakness. 

The results of Theorems D.2 and D.3 and the corollary B.7 are important as they ensure that the pushforward maps can be written in terms of the Radpm transform operator R. 

The paper is written in great detail and provides appropriate background for the concepts that are need to understand it. This is a strength. 

The experimental results are somewhat comprehensive  and the authors do a good job of comparing their method with other spherical slicing methods (Wasserstein etc.). 

The authors also provide a time complexity and runtime analysis and compare them to those with other related methods.
summary: This paper combines ideas from sliced optimal transport theory and linear optimal transport and introduces the idea of linear spherical sliced optimal transport. Specifically, they extend a recently proposed work (Martin et al. 2023) that introduced the linear circular optimal transport for embedding  circular distributions into a L^2 space. In this paper, the authors propose a new embedding method for spherical probability distributions. 

The authors present experimental results for the method of linear spherical sliced optimal transport (LSSOT) in applications such as point cloud registration and cortical surface matching.
weaknesses: This is a minor weakness, but as mentioned in the strengths above, the LSSOT is a metric only for probability measures. Thus the data has to be first converted to distributions, which is both a lossy transformation (which the authors acknowledge). 

By definition, the projection is not invertible. That is one cannot recover the original signal/shape. Again, this drawback exists for all the other methods, including those that the authors build upon. 

Because of the above 2 points, when the authors present the interpolation results, they resort to an auto-encoder that projects the original data to the nearest point in the spherical latent space. Then to reconstruct the intermediate points along the interpolation path between the two signal, they need a decoder. This is similar to the backprojection approaches, which are not ideal. But again this is not a major weakness, but currently, this is the most convenient approach. 

In addition to the above loss of information, the projection step (slicing) will also introduce further discrepancy. However, this is a drawback for all the slicing based methods (including spherical) that are mentioned in the paper. In this paper, the authors use the same idea as Bonet et al. 2023 for performing the projection. 

While the authors do provide algorithms for computing flows using Riemannian gradient descent, this is technically not a geodesic, but an approximation. 

The notation for P^C (Equation 6) and P^U (Equation 7) gets confusing. The first one is a geodesic distance, while the second one is a projection.

---

confidence: 4
contribution: 2
initial_rating: 6
presentation: 2
questions: To discuss:
-  Invariance of this metric to nuisances such as rotation. 
-  Is this a metric or pseudo-metric, if a pseudo-metric (as in proposition D1?) why not state this in the title?
-  Discuss invertible spherical harmonics, filter banks, which can be used to form a metric and are used on spherical brain data [1].
-  the order of computational complexity of this metric in a practical spherical data retrieval scenario. Apparently it would require O(N) search time in N the number of data, other fast methods (non spherical) used on the same ADNI brain data require O(log N) [2].

Please demonstrate the utility by identifying the incorrectly labelled ADNI subjects[1], these will immediately stand out as unusually close according to a valid metric as they are the same subjects incorrectly labeled as different[2].

The paper states “the Dice loss to impose biological validity”, whereas Dice and other tissue overlap measures specifically do not enforce biological validity as they permit arbitrary biologically implausible mappings[3].

[1] Yeo, B. T. T., Ou, W., & Golland, P. (2008). On the construction of invertible filter banks on the 2-sphere. IEEE Transactions on Image Processing, 17(3), 283-300.

[2] Chauvin, Laurent, et al. "Neuroimage signature from salient keypoints is highly specific to individuals and shared by close relatives." NeuroImage 204 (2020): 116208.

[3] Rohlfing, Torsten. "Image similarity and tissue overlaps as surrogates for image registration accuracy: widely used but unreliable." Medical Imaging, IEEE Transactions on 31.2 (2012): 153-163.
soundness: 3
strengths: The math appears interesting and valid. Spherical data are of interest in various contexts.
summary: A fast metric is proposed for spherical data, based on sliced optimal transport.
weaknesses: While the metric appears compariatively fast in comparisons of pairwise distance calculation, it does not address the order of complexity of performing many distance calculations given a dataset of N points, which is often the limiting factor in practical contexts. For example here nearest neighbor search would be O(N), whereas other methods particularly indexing in the context of brain data similary using the same ADNI data offer O(logN) performance.

Meta-review: This paper presents the Linear Spherical Sliced Optimal Transport (LSSOT) framework, which efficiently compares spherical probability distributions by embedding them into spaces while preserving intrinsic geometry. LSSOT offers computational advantages over existing methods and demonstrates high accuracy in applications like cortical surface registration, 3D point cloud interpolation, and shape embedding.

Reviewers all agreed on the merits of this paper, I am therefore recommending an accept for it. Congratulations for the good work !

Decision: Accept (Spotlight)