Fast and Accurate Intersections on a Sphere

TL;DR

A simplified intersection point formula with improved speed and numerical robustness over the ones traditionally implemented in geoscience software is proposed and algorithms based on the concept of error-free transformations (EFT) can be applied to evaluate this formula within a relative error bound that is on the order of machine precision.

Abstract

We introduce a fast, high-precision algorithm for calculating intersections between great circle arcs and lines of constant latitude on the unit sphere. We first propose a simplified intersection point formula with improved speed and numerical robustness over the ones traditionally implemented in geoscience software. We then show how algorithms based on the concept of error-free transformations (EFT) can be applied to evaluate this formula within a relative error bound that is on the order of machine precision. We demonstrate that, with a vectorized and parallelized implementation, this enhanced accuracy is achieved with no compute-time overhead compared to a direct calculation in hardware floating point, making our algorithm suitable for performance-sensitive applications like regridding of high-resolution climate data. In contrast, evaluating our formula using high-precision data types like quadruple precision and arbitrary precision, or using the robust intersection computation routines from the Computational Geometry Algorithms Library (CGAL), leads to significant computational overhead, especially since these alternatives inhibit vectorization. More generally, our work demonstrates how EFT techniques can be combined and extended to implement nontrivial geometric calculations with high accuracy and speed.

Figures (7)

• Figure 1: The red circle represents the great circle arc defined by the plane $P_{\hat{n}}$, which passes through the origin with normal vector $\hat{n}$. The constant-latitude circle is drawn in blue. The vector $\hat{e}_z$ is aligned with the z-axis. Also visualized are the vectors ${\bf t}$ and ${\bf t}^\perp$ introduced in \ref{['sec:math']}, which are tangent and perpendicular to the planes' intersection line, respectively.
• Figure 2: Accuracy comparison of precision models (MPFR, quadruple-precision, and float64) for various intersection equations. The float64 + EFT (AccuX) represents our proposed AccuX algorithm. The first row shows accuracy using \ref{['eq:gca_constlat_intersection_point_baseline']} across different precision models, including results from Mirshak's Python code (2024) krumm2000mirshak2024intersections. The second row presents results with \ref{['eq:gca_constlat_intersection_point_cdo']}, and the third row uses our final proposed equation \ref{['eq:final_gca_constLat_intersection_point']}. Since AccuX is based on \ref{['eq:final_gca_constLat_intersection_point']}, it appears only in the third row. For comparison, CGAL's double-precision approximation is shown in the last row alongside AccuX.
• Figure 3: Accuracy analysis of precision models near the North Pole (90N) using \ref{['eq:final_gca_constLat_intersection_point']} for MPFR, quadruple-precision, float64 (IEEE754 double precision), and CGAL’s double-precision approximation implementations. The float64 + EFT (AccuX) represents our proposed AccuX algorithm, demonstrating superior accuracy in the critical polar region compared to other methods.
• Figure 4: Accuracy comparison of precision models (MPFR, quadruple-precision, and float64) for \ref{['eq:gca_constlat_intersection_point_baseline']} and \ref{['eq:final_gca_constLat_intersection_point']}, alongside Krumm’s (2000) method krumm2000, in the calculation of ill-conditioned cases.
• Figure 5: Performance comparison of different precision models, including MPFR, quadruple-precision (QP), and float64 (FP64) implementations, utilizing three mathematical formulas: \ref{['eq:gca_constlat_intersection_point_baseline']}, \ref{['eq:gca_constlat_intersection_point_cdo']}, and \ref{['eq:final_gca_constLat_intersection_point']}. FP64 + EFT (AccX) represents our proposed algorithm based on \ref{['eq:final_gca_constLat_intersection_point']}. CGAL presents the CGAL's double-precision approximation. Lastly, CGAL computes both intersection points, while all other methods compute only one, but with less than twice the workload due to shared intermediate results.