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

We thank the reviewer for their insightful questions and comments. Below are our responses.

1. Invariance to rotation.

Answer: We appreciate the Reviewer pointing out another nice property of LSSOT. Indeed, LSSOT is rotationally-invariant, i.e. LSSOT(\mu,\nu)=LSSOT(R_\\#\mu,R_\\#\nu) for any rotation $R$, as it is defined as the expectation over the uniform distribution of the stiefel manifold (Equation 10). We empirically test the invariance in the updated Appendix L. For a pair of distributions, we randomly generate a set of rotations and LSSOT stays almost the same for the rotated pairs in Figure 33.

2. ... metric or pseudo-metric ...

Answer: The LSSOT discrepancy is defined in Eq. (10). In Theorem 3.3 (or Theorem D2), we establish its metric properties: it is indeed a distance for probability measures with continuous density functions. However, for general spherical probability measures, it is a pseudo-metric, as noted in Theorem 3.3 or Theorem D1, as pointed out by the referee. We adopted the term "LSSOT distance" following the statement of Theorem 3.3.

3. ... invertible spherical harmonics, filter banks ...

Answer: We thank the reviewer for pointing this out and will include the following discussion:

Classical tools from harmonic analysis used in signal processing typically rely on linear techniques, such as representing signals using bases for expressing signals (for example, signals are supported on the 2-dimensional sphere can be written in terms of spherical harmonics, which form a basis for $L^2(\mathbb S^2)$) and using linear operators for analysis or filtering (see, for example, [Yeo et al., 2008]). These classical techniques do not require normalization (whereas our approach operates on inputs as probability measures) and provide invertible frameworks. One of the most important applications of these methods is denoising, but they do not provide a meaningful distance for tasks that require pairwise comparisons of signals. In contrast, LCOT embedding is a non-linear operator and consequently LSSOT is also a non-linear operator. Rather than relying on the Euclidean distance in $L^2(\mathbb S^{d-1})$ between the inputs, this new approach introduces a novel distance after embedding.

4. the order of computational complexity of this metric in a practical spherical data retrieval scenario ...

Answer: There may be a misunderstanding or miscommunication that we would like to clarify. Consider $K$ training distributions on $\mathbb{S}^{d-1}$ (i.e., spherical distributions), each with $N$ particles (i.e., support size $N$). The computational cost for retrieving the nearest neighbor (NN) to a test distribution depends on two factors: (1) the cost of computing the distance between distributions, and (2) the nearest neighbor search algorithm used. For a Linear Search algorithm (i.e., brute force), the NN retrieval has a complexity of $\mathcal{O}(D K)$, where $D$ is the cost of computing the pairwise distance. This complexity can be reduced to $\mathcal{O}(D \log K)$ using KD-Trees, Approximate NN (ANN) methods, or other efficient search structures. In this work, our contribution lies in reducing the cost of $D$.

Specifically, performing NN retrieval using Optimal Transport (OT), even with efficient search methods, incurs a complexity of $\mathcal{O}((N^2 d + N^3 \log N) \log K)$. In contrast, our method reduces this cost to $\mathcal{O}(L N (d + \log N + 1 + \log K))$. This reduction is achieved as follows: embedding the test distribution into our $L N$-dimensional embedding space requires $\mathcal{O}(L N (d + \log N + 1))$, and the subsequent NN retrieval over $K$ training distributions requires $\mathcal{O}(L N \log K)$, resulting in the stated overall complexity.

5. Please demonstrate the utility by identifying the incorrectly labelled ADNI subjects ...

Answer: We are thankful for the suggestion on expanding the utility of our work. We totally agree that our method will be robust w.r.t. the incorrect labelling. We design a simplified experiment to demonstrate this in Appendix M. We randomly select 25 subjects from the ADNI dataset and duplicate one of them (the 8-th subject) as the test subject with an incorrect label. Figure 34 depicts that $e^{\text{LSSOT}}$ stands out at subject 8 with value 1, as LSSOT is arbitrarily close to 0. Such behavior makes LSSOT a very convenient tool in correcting the labels in the preprocessing of the datasets.

6. Dice loss and biological validity

Answer: We agree with the Reviewer that Dice Score is not the gold standard in evaluating biological validity, even with human-experts labeled parcellations as mentioned in Rohlfing, Torsten. (2012). For lack of better alternatives, we use Dice Score following the experimental setups from established registration works. Although the registration framework is not the novelty nor the focus of our work, we are open to any suggestions in future work on the registration problems.

We thank the reviewer for their constructive feedback, and their insightful comments. Below we provide our responses.

Weaknesses

1. We greatly appreciate the Reviewer's thorough understanding of our work and pointing out the limitations, which we wholeheartedly acknowledge. These limitations will work as motivations for future extensions of this work.
2. In Equation 6, $P^C(x)$ is defined through argmin, as the closest point $y^*$ on the great circle $C$ to $x$ with respect to the spherical geodesic distance. Therefore, $P^C(x)$ is a projection from $\mathbb{S}^{d-1}$ to the great circle $C$. $P^C(x)$ is of a more general form where the projections of north/south poles (Figure 1) are well-defined while in $P^U$ they are not. To avoid confusion of notations, we have rephrased and removed the $P^C$ notation in Equation 6. We appreciate the help with the succinct presentation of the paper.

Questions:

1. min{|x−y|,1 −|x−y|} and min(arccos(<x, y>))

Answer: One can parametrize the points on the unit circle as angles between 0 and 1 and view this space as the quotient space $\mathbb{R} / \mathbb{Z}$, which inherits the metric $\min\{|x - y + k| : k \in \mathbb{Z}\}$. This is equivalent to $\min\{|x - y|, 1 - |x - y|\}$. This metric measures the shortest arc length between two points $x$ and $y$ on the unit circle when parameterized by angles in the interval $[0, 1)$. In this setting, $|x - y|$ gives the absolute difference of their parameters, corresponding to their distance if there is no 'wrap-around' on the circle, while $1 - |x - y|$ accounts for the 'wrap-around'. Thus, $\min\{|x - y|, 1 - |x - y|\}$ gives the shortest arc length on $\mathbb{S}^1$. When we parametrize points on the unit circle using angles between 0 and 1, we translate these points into the standard representation on the circle via $(\cos(2\pi x), \sin(2\pi x))$ and $(\cos(2\pi y), \sin(2\pi y))$. As for the spherical distance $\min(\arccos(\langle \alpha, \beta \rangle))$, it can be used for any sphere $\mathbb{S}^{d-1} \subset \mathbb{R}^d$, where $\alpha$ and $\beta$ are unit vectors in $\mathbb{R}^d$, not just for $\mathbb{S}^1$. However, for the unit circle, we have: $\min(\arccos(\langle \alpha, \beta \rangle)) = 2\pi \cdot \min\{|x - y|, 1 - |x - y|\}$, where $x, y$ are angles parametrized between 0 and 1, and $\alpha = (\cos(2\pi x), \sin(2\pi x))$, $\beta = (\cos(2\pi y), \sin(2\pi y))$.

2. ... why the complexity of spherical sliced Wasserstein better than sliced Wasserstein...

Answer: We appreciate the precise summary on the computation efficiency of LSSOT following from slicing and LCOT. We point out that the complexity of spherical sliced Wasserstein ($\mathcal{O}(LN(d+\log \frac{N}{\epsilon})))$) is not better than sliced Wasserstein ($\mathcal{O}(LN(d+\log N))$) and it's absolutely true that projections on $\mathbb{S}^1$ should be more costly than on lines. Albeit slightly faster, the sliced Wasserstein is defined only for the Euclidean geometry, hence not respecting the underlying geometry of the space. In our work, we include SWD as a baseline to showcase the importance of preserving spherical geometry in the context of spherical data analysis. A detailed explanation of the LSSOT complexity can be found in Section 3.3 and a more complete complexity table has been added in Appendix J.

3. ... why LSSOT has better accuracy than SSW ...

Answer: As shown in Table 1, SSW performs better than LSSOT in about half of the evaluations, therefore in the results analysis we concluded that SSW and LSSOT are on par with each other in the overall registration, instead of claiming LSSOT provides better results than SSW. However, we emphasize that despite the fact that SSW and LSSOT are close to each other, LSSOT itself is an independent metric with its own metric properties and characteristics on the space of spherical distributions. To provide some insights, we refer the reviewer to [1] for behavior comparisons of COT and LCOT, which accounts for the major difference between SSW and LSSOT.

4. ...The injectivity of R...

Answer: In Bonet et al. 2024., the authors do not prove/disprove the injectivity of the Spherical Radon transform. In this paper, we prove that the spherical Radon transform in Bonet et al. 2024., is equivalent to the "Hemi-sphere Radon transform" in [Groemer, 1998, Section 2] and the "the Un-normalized Semicircle Transform" in [Hielscher et al. (2018)], which have been proved to be injective. Regarding the operator R, the goal of adapting Radon transform in this paper is for computational efficiency. The spherical Radon transform can preserve the geometry information. Regarding the injectivity property, it can determine if the proposed distance, i.e. SSW or LSSOT is a metric.

[1] Rocio Diaz Martin, Ivan Medri, Yikun Bai, Xinran Liu, Kangbai Yan, Gustavo K Rohde, and Soheil Kolouri. LCOT: Linear circular optimal transport. arXiv preprint arXiv:2310.06002, 2023.

We thank the reviewer for their time and constructive feedback. Below, we provide our responses.

1. From a technical point of view, the extension of the proposed method to the existing LCOT is limited. Still, it should be noted that it is useful in spherical signal analysis applications.

Answer: We appreciate the reviewer's comment. Our method builds on the existing LCOT tool for circular probability measures, as mentioned in the introduction of our manuscript prior to outlining our contributions. By employing slicing techniques, we generalize this approach to probability measures on higher-dimensional spheres. Additionally, as noted by the reviewer, the proposed method has practical applications in spherical signal analysis. Beyond the experimental results, we have also established the metric properties of our proposed tool for continuous probability measures.

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.

Answer: We thank the reviewer for the constructive suggestion. Indeed, this will help strengthening the advantage of our method over all baselines in the statistical aspect. To show the significance of difference/improvements on Dice score, area/edge distortions, we have performed significance testing and calculate the p values. All results are significant with p-values much lower than 0.05, and in the table we further specify entries with p<0.01 and p<0.001.

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.

Answer: Here, we provide the computational complexities for $d$-dimensional discrete distributions with a support size of $N$. Let $L$ denote the number of slices. The slicing and projection steps require $\mathcal{O}(LNd)$. For a single slice, computing the LCOT embedding involves sorting the $N$ samples, which has a complexity of $\mathcal{O}(N \log N)$, and averaging them to find the antipodal point on $\mathbb{S}^1$, with a complexity of $\mathcal{O}(N)$. Combining these, and including the non-dominant terms for clarity, the total complexity for calculating a single embedding is $\mathcal{O}(LN(d + \log N + 1))$.

When comparing $K$ spherical probability measures, our method requires $K$ embeddings, with a complexity of $\mathcal{O}(KLN(d + \log N + 1))$, followed by $\frac{K(K-1)}{2}$ Euclidean distance calculations between $LN$-dimensional embeddings, which requires $\mathcal{O}(K^2LN)$. Thus, the total complexity for pairwise comparisons is $\mathcal{O}(KLN(d + \log N + 1 + K))$.

The table below compares the complexities for LSSOT, SSW, and OT, both for a single comparison and for comparing $K$ spherical probability measures:

A more detailed complexity table has been added in Appendix J of the updated manuscript.

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.

Answer: We add in Appendix K the evolution of LSSOT with respect to the number of projections and reference sizes for different dimensions $d$ and VMF distributions with different $\kappa$. As shown in Figure 31 and Figure 32, LSSOT is pretty robust with more than 100 slices and all reference sizes from 100 to 5000. We would also like to clarify that $\lambda$ from Equation (12) is the regularization parameter of the autoencoder loss in the point cloud interpolation, which is not a hyperparamter for the LSSOT method. The optimal value for $\lambda$ is selected by grid search over the set $[10^{-5}, 10^{-4},10^{-3}, 10^{-2}]$ to achieve the best trade-off between reconstrution and regularity in the latent space.

We sincerely thank the reviewer for taking the time to provide us with constructive feedback. Below, we have detailed our responses.

1. Did the authors perform new hyperparameter tuning on these losses when switching the image similarity metric?

Answer: Yes, we carefully perform hyperparameter tuning in each configuration using grid search on the validation dataset. For the regularization parameter, the weight space is set to be higher than the other terms, ranging over [10, 50], as enforcing spatial smoothness is a high priority in diffeomorphic registration. For the parcellation term, we sweep the parameter in the range of $[10^{-3}, 10^{-2}, 10^{-1}, 1, 10]$. We follow the Voxelmorph work [1] and select the best hyperparameter of a total of 25 combinations in the Dice Score evaluation on the validation dataset.

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

Answer: We appreciate the reviewer pointing out that curvature is another common feature that is used in cortical surface registration. Although it is not required to study curvature in cortical surface registration works [2], we are happy to include extra experiments on the curvature registration in the appendix of the final draft to expand the utility of our method. We emphasize that the focus of our work is on the development of the LSSOT method, and the sulcal depth registration serves as one application to sufficiently show the effectiveness of LSSOT.

3. 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?

Answer: In Table 1, the columns are methods and the rows are evaluation metrics. Thus in the SWD evaluation on NKI(left), our method LSSOT gives 0.0052±0.0011 compared to 0.0051±0.0017 of MSE which is only slightly better. We would also like to note that the underlying ground metric of SWD is Euclidean, whereas LSSOT preserves spherical geometry, therefore SWD and LSSOT are not necessarily more correlated than SWD and MSE. In NKI(right), it is expected that MSE method performs the best in MAE, as they are highly consistent with each other.

4. 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.

Answer: We agree with the reviewer about the inherent trade-off between the matching objective and the smoothness of the registration map. To achieve this trade-off, we follow the established deep learning registration framework [3] and the setup there-in which is also inspired by the Voxelmorph work [1]. In this setting, the trade-off between the similarity and regularity is achieved by validating using the Dice Score metric. As mentioned in the first question, we conduct grid search for the optimal hyperparameters.

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

Answer: We appreciate the Reviewer’s suggestion and have included a table with FreeSurfer results for the NKI dataset (left hemisphere) in the temporary Appendix N of the updated manuscript. As the Reviewer correctly noted, FreeSurfer is significantly slower. We will complete the analysis for all four datasets and update Table 1 in the main text for the final version of the paper. Additionally, it is worth noting that deep-learning-based frameworks can register unseen test data in a single feed-forward pass, whereas FreeSurfer iteratively solves for each pair of scans. This allows deep-learning-based frameworks to effectively amortize the registration optimization process and significantly speed up the registration. Lastly, we would like to reiterate that the primary contribution of this work lies in the development of a mathematically rigorous and scalable metric for spherical probability measures, with cortical surface registration serving as a demonstration of its practical application.

[1] Guha Balakrishnan, Amy Zhao, Mert Rory Sabuncu, John V. Guttag, and Adrian V. Dalca. Voxelmorph: A learning framework for deformable medical image registration. IEEE Transactions on Medical Imaging, 38:1788–1800, 2018.

[2] Mohamed A. Suliman, Logan Z. J. Williams, Abdullah Fawaz, Emma. C Robinson. GeoMorph: Geometric Deep Learning for Cortical Surface Registration. Proceedings of the First International Workshop on Geometric Deep Learning in Medical Image Analysis, PMLR 194:118-129, 2022.

[3] Fenqiang Zhao, Zhengwang Wu, Fan Wang, Weili Lin, Shunren Xia, Dinggang Shen, Li Wang,and Gang Li. S3reg: Super fast spherical surface registration based on deep learning. IEEE Transactions on Medical Imaging, 40(8): 1964–1976, 2021.

Thank you for your detailed rebuttal and the new experiments comparing to FreeSurfer. I've adjusted my score to 8.