Detecting Cocircular Subsets of a Spherical Set of Points

Abstract

Given a spherical set of points, we consider the detection of cocircular subsets of the data. We distinguish great circles from small circles, and develop algorithms for detecting cocircularities of both types. The suggested approach is an extension of the Hough transform. We address the unique parameter-space quantization issues arising due to the spherical geometry, present quantization schemes, and evaluate the quantization-induced errors. We demonstrate the proposed algorithms by detecting cocircular cities and airports on Earth’s spherical surface. These results facilitate the detection of great and small circles in spherical images.

1. Introduction

Given a spherical set of points, consider the problem of detecting the largest subset of cocircular points, i.e., the largest subset of points lying on the same circle on the sphere. One brute-force solution, related to RANSAC [1], is to examine all possible triplets of points in the point set. Note that each triplet (if not colinear) defines a plane, the intersection of the plane with the sphere defines a circle, and each point in the triplet clearly lies on the circle. The point set can then be searched for additional points lying on that circle. This $O\left(N^4\right)$ procedure, where N is the cardinality of the point set, is prohibitive for sizable point sets.

A cocircular subset of a spherical set of points is obviously coplanar. In special cases, the plane passes through the center of the sphere. In such cases, the intersection of the plane and the sphere is referred to as a great circle; see the green circle in Figure 1 (right). In the general case, when the plane does not pass through the center of the sphere, the intersection of the plane and the sphere is referred to as a small circle; see the blue and red circles in Figure 1 (right).

In the classical Hough transform [2], the problem of finding colinear points in the plane is replaced by the dual problem of finding concurrent lines or curves in a parameter space, often referred to as the Hough space. In practice, the parameter space is quantized and represented by an accumulator array, a voting process takes place, and finding the intersection of lines or curves in the parameter space is implemented as a search for the maximum. The Hough transform is fundamentally robust with respect to outliers, and has been formally associated with the M-estimation methodology for robust regression [3]. The Hough transform has also been related to conformal geometric algebra [4]. The time-complexity of the classical Hough transform is linear in the cardinality N of the dataset. Thus, when looking for the largest co-linear subset in a planar set of points, the Hough transform reduces the computational load from $O\left(N^3\right)$ (the cost of checking each possible pair of points in the data; determining the line that connects the pair; and, for each other data-point, testing whether it lies on that line) to $O\left(N\right)$.

The suggested approach to finding the largest subset of cocircular points in a spherical set of points is inspired by the Hough transform. Specifically, we characterize the circles (on the sphere) passing through a given point in the spherical set. The characterization is in the form of an equation constraining the parameters of the circles. We then search for the largest intersection of such equations, i.e., for the parameters specifying the circle passing through the largest possible number of points in the set.

As in the classical Hough transform, once the best-supported circle, i.e., the circle passing through the largest possible number of points, has been determined, finding the second-best-supported circle is straightforward. This can be reliably accomplished by identifying the data-points supporting the best circle, removing them from the dataset, and eliminating their contribution to the accumulator array (“unvoting”) [5,6,7]. Circle detection and unvoting can be iterated until all circles are found.

The precise function of the Hough transform is the detection of a geometric primitive, such as a straight line or circle, common to the largest possible subset of points in the dataset. Identifying the relevant points themselves can then be easily accomplished [5], but is not always necessary. Note that it is commonly said that the Hough transform detects geometric primitives in a digital image, but accomplishing this task in an actual image requires edge detection as a preprocessing step.

Images acquired using omnidirectional cameras can be represented as spherical images. Straight line segments in the scene are then transformed to great circle arcs in the spherical image. This has led to some interest in the detection of great circles in spherical images. With this motivation, Vasseur and Mouaddib [8] and Torii and Imiya [9] outlined a Hough algorithm and a randomized Hough algorithm, respectively, for great-circle detection on a sphere. These ideas have modern applications, such as road-line detection in driver-assistance systems [10].

In the detection of vanishing points, the search for convergence points in the image plane can be usefully replaced by the intersection of great circles on the Gaussian sphere [11]. This important task has been addressed using specialized Hough algorithms, leading to developments in spherical tessellation and discretization schemes [12,13] that are also useful in great-circle detection [14,15].

In geophysics, interest in great and small circle fitting arises in the context of volcano distribution analysis [16,17]. In these studies, specialized map projections transform great circles into an exact (gonomonic projection [17]) or approximate (UTM [16]) straight lines, allowing the use of the straight-line Hough transform [2] for great-circle detection. In the context of plate tectonics, Wessel [18] outlined the principles of a Hough transform for great-circle detection in the true spherical domain and sketched its generalization to small-circle detection. Wessel’s insightful concepts [18] have so far received surprisingly little attention.

In this paper, we present comprehensive solutions for the detection of cocircular subsets of a spherical set of points. We describe a geometric framework that supports both great- and small-circle detection and provide algorithms for both cases. We address the unique parameter-space quantization issues arising due to the spherical geometry, present quantization schemes, and evaluate the quantization-induced errors. We demonstrate the proposed algorithms by detecting cocircular cities and airports on Earth’s spherical surface.

2. Fundamentals

In this article, given a set of $N>>1$ co-spherical points, we present a way to detect the circle that passes through the largest possible number of points in the set. It is easy to show that the circle itself necessarily lies on the sphere.

A point in 3D space can be defined by its spherical coordinates, consisting of the radial distance R from the origin, the polar angle $\theta$, and the azimuthal angle $\phi$; see Figure 1 (left). Assuming co-spherical points, point i is determined by the angular coordinates ${\theta }_i$ and ${\phi }_i$ alone.

Given an anchor point in 3D space and a normal direction, one can determine a plane that passes through the anchor point and fits the normal direction. By itself, a normal direction defines a set of parallel planes.

Viewing the angular coordinates of a (primary) point on the sphere as a normal direction therefore associates a set of parallel planes with that point. The intersection of each such plane with the sphere is a (primary) circle, as shown in Figure 1 (right).

Any (secondary) point on a primary circle of radius r (point $S_1$ in Figure 2), viewed as a normal direction, defines a secondary circle of radius r passing through the primary point (point P in Figure 2). Thus, given a primary point and a radius r, each secondary point is the center of a circle of radius r passing through the primary point.

Since this property holds for any secondary point on a primary circle, a set of primary circles for primary points located on the same circle will intersect at a specific secondary point, which is the center of the circle passing through all the primary points. This is illustrated in Figure 3, where the red, blue, and green points are primary points on the same (black) circle of radius r. The corresponding primary circles of radius r intersect at the black secondary dot S, which defines the secondary circle (black) passing through all the primary points.

Now, given a set of unit co-spherical points, i.e., a set of points $S=\left\{({\theta }_i,{\phi }_i)\right\}$ located on the unit sphere, we look for the normal direction $(\theta ,\phi )$ that determines the circle of (pre-defined) radius r passing through the largest possible subset of points in S. The problem amounts to finding the largest co-circular (with given radius r) subset of S. We approach it by viewing each point $({\theta }_i,{\phi }_i)$ as a primary point corresponding to a primary circle of radius r. We view $(\theta ,\phi )$ as a secondary point corresponding to a secondary circle on which a point $({\theta }_i,{\phi }_i)$ lies. We select $(\theta ,\phi )$, compatible with the largest possible number of primary points $({\theta }_i,{\phi }_i)$. In other words, we select $(\theta ,\phi )$, which is the largest intersection of primary circles.

As just presented, the problem and its solution strategy can be regarded as a generalization of the classical Hough transform for circles of known radius in the plane. In principle, one might indeed solve the problem by allocating a memory array on a spherical surface and incrementing memory cells on the primary circles corresponding to the data-set $S=\left\{({\theta }_i,{\phi }_i)\right\}$ (voting). The memory cell with the largest accumulation will then correspond to $(\theta ,\phi )$, which is the largest intersection of primary circles, i.e., the normal direction defining the circle of radius r passing through the largest possible subset of points in S.

In practice, working with spherical memory arrays is highly inconvenient. We therefore carry out the computation using a planar memory array, where the orthogonal planar axes are $\theta$ and $\phi$. Note, however, that the mapping of the conceptual spherical memory array to a planar array is non-trivial in both theory and practice, as is the drafting of planar maps of planet Earth [19]. As will be seen, this difficulty reveals itself in the context of the parameter-space quantization.

We can summarize these concepts as follows: Given a (unit) sphere and a radius r,

Property 1.

A primary point $({\theta }_i,{\phi }_i)$ on the sphere corresponds to a curve in the $(\theta ,\phi )$ parameter plane. The curve can be viewed as the mapping of the primary circle associated with the primary point from the sphere to the $(\theta ,\phi )$ plane.

Property 2.

A point in the $(\theta ,\phi )$ parameter plane corresponds to a secondary circle of radius r on the sphere.

Property 3.

Primary points lying on the same circle of radius r on the sphere correspond to curves through a common point in the $(\theta ,\phi )$ parameter plane.

In the next section, we determine $\theta$ as a function of $\phi$ (or $\phi$ as a function of $\theta$ ). This is the voting pattern associated with a data point $({\theta }_i,{\phi }_i)\in S$.

3. The Voting Pattern

3.1. Derivation of the Voting Pattern

Consider a primary point $(X_i,Y_i,Z_i)$ on a unit sphere. We proceed to express the equation of the corresponding family of primary circles lying on parallel planes, as illustrated in Figure 1 (right). Generally, the equation of a plane in 3D space is of the form

$$ax+by+cz=d$$

(1)

Note that $(a,b,c)$ is the normal to the plane, and if it is normalized such that $a^2+b^2+c^2=1$, then d is the distance of the plane from the origin.

Using Pythagoras’ theorem, for a unit sphere, the relation between d and the radius r of the circle is (see Figure 4):

$$r=\sqrt{1-d^2}$$

(2)

In Equation (1), $(a,b,c)$ is the normal to the plane. Here the normal is defined by the primary point. In terms of its angular coordinates $({\theta }_i,{\phi }_i)$,

$$\overrightarrow{n}=\left(\begin{array}{ccc}a& b& c\end{array}\right)=\left(\begin{array}{ccc}{X}_i& Y_i& Z_i\end{array}\right)=\left(\begin{array}{ccc}sin{\theta }_{i}cos{\phi }_i& sin{\theta }_{i}sin{\phi }_i& cos{\theta }_i\end{array}\right)$$

(3)

Thus, the equation of the plane containing the primary circle is:

$$\left(\begin{array}{ccc}sin{\theta }_{i}cos{\phi }_i& sin{\theta }_{i}sin{\phi }_i& cos{\theta }_i\end{array}\right)\left(\begin{array}{c}x\\ y\\ z\end{array}\right)=d$$

(4)

The primary circle itself is the intersection of this plane with the unit sphere. In spherical coordinates, a point $(x,y,z)$ on the unit sphere is

$$\left(\begin{array}{c}x\\ y\\ z\end{array}\right)=\left(\begin{array}{c}sin\theta cos\phi \\ sin\theta sin\phi \\ cos\theta \end{array}\right)$$

(5)

Substituting in the plane Equation (4) we obtain

$$\left(\begin{array}{ccc}sin{\theta }_{i}cos{\phi }_i& sin{\theta }_{i}sin{\phi }_i& cos{\theta }_i\end{array}\right)\left(\begin{array}{c}sin\theta cos\phi \\ sin\theta sin\phi \\ cos\theta \end{array}\right)=d$$

(6)

Rearranging this, we obtain a relation between the angular coordinates $\theta ,\phi$ of secondary points on the primary circle.

$$cos\theta cos{\theta }_i+sin\theta cos\phi sin{\theta }_{i}cos{\phi }_i+sin\theta sin\phi sin{\theta }_{i}sin{\phi }_i=d$$

$$cos\theta cos{\theta }_i+sin\theta sin{\theta }_i(cos\phi cos{\phi }_i+sin\phi sin{\phi }_i)=d$$

$$cos\theta cos{\theta }_i+sin\theta sin{\theta }_{i}cos(\phi -{\phi }_i)=d$$

$$cos(\phi -{\phi }_i)=\frac{d}{sin\theta sin{\theta }_i}-cot\theta cot{\theta }_i$$

(7)

$$\begin{array}{c}\hfill {\phi }^I=2\pi n+{cos}^{-1}(\frac{d}{sin\theta sin{\theta }_i}-cot\theta cot{\theta }_i)+{\phi }_i\\ \hfill {\phi }^{II}=2\pi n-{cos}^{-1}(\frac{d}{sin\theta sin{\theta }_i}-cot\theta cot{\theta }_i)+{\phi }_i\end{array}\hspace{1em}\hspace{0.17em}\hspace{0.17em}n\in \mathbb{Z}$$

(8)

For a given primary point $({\theta }_i,{\phi }_i)$ and a plane distance d, Equation (8) determines the azimuthal angle $\phi$ of a secondary point in terms of the polar angle $\theta$. The two solutions reflect the fact that a circle of latitude (defined by $\theta$) that intersects a primary circle (excluding the two circles of latitude that osculate the primary circle) intersects the primary circle at two points; see Figure 2 and Figure 3. The dependence on n is due to the periodicity of $\phi$. Figure 5 shows the locus of {$\theta$,$\phi$} pairs for selected values of ${\theta }_i$, ${\phi }_i$, and d.

As shown in Figure 5 (top), for the case of a great circle $(d=0)$, we can express $\theta$ as a single-valued function of $\phi$. Substituting $d=0$ in Equation (7),

$$cos(\phi -{\phi }_i)=-cot\theta cot{\theta }_i$$

(9)

Solving for $\theta$,

$$\theta =\pi n-{cot}^{-1}(cos(\phi -{\phi }_i)tan{\theta }_i)\hspace{1em}n\in \mathbb{Z}$$

(10)

3.2. Domain and Range of the Voting Pattern

We have so far seen that for a given primary point $({\theta }_i,{\phi }_i)$, the $\phi$ coordinate of a secondary point $(\theta ,\phi )$ is obtained by substituting the $\theta$ coordinate in Equation (8). In this subsection, we proceed to show the dependence of the domain and the range of $\phi$ as a function of $\theta$, as represented by Equation (8), on the primary point coordinates $({\theta }_i,{\phi }_i)$ and the distance d.

As shown in Figure 4, for a given distance $d>0$, each (spherical) circle of radius r (Equation (2)) on a unit sphere is uniquely defined by one (spherical) point on the unit sphere. Any such spherical point can be represented by the spherical coordinates $\theta \in [0,\pi )$ and $\phi \in [0,2\pi )$. Therefore, any spherical circle with radius $r<1$ is uniquely defined by $\theta \in [0,\pi )$ and $\phi \in [0,2\pi )$.

In the case of great circles, $d=0$, the points $(\theta ,\phi )$ and $(\pi -\theta ,\pi +\phi )$ (red points in Figure 6) define the same plane. In other words, each great circle (red) is represented by two spherical points, one at each hemisphere. Consider the subset of spherical points in the range $\theta \in [0,\pi )$ and $\phi \in [0,\pi )$. There is a bijective transformation from this subset to the set of great circles on the unit sphere. To conclude, each spherical point with $\theta \in [0,\pi )$ and $\phi \in [0,\pi )$ uniquely defines a great circle.

Returning to the general case, for a given primary point $({\theta }_i,{\phi }_i)$, the spherical points $(\theta ,\phi )$ that satisfy Equation (8) represent the secondary circles on which the primary point lies. We are interested in secondary points with $\phi \in [0,2\pi )$, so n is selected such that ${\phi }^I,{\phi }^{II}\in [0,2\pi )$.

Given a primary point and d, the domain of Equation (8) consists of $\theta$ coordinates of the secondary points. The $\theta$ value of the secondary points are such that the argument of the inverse cosine function satisfies

$$\frac{d}{sin\theta sin{\theta }_i}-cot\theta cot{\theta }_i\in [-1,1]$$

$$\left|\frac{d}{sin\theta sin{\theta }_i}-cot\theta cot{\theta }_i\right|\le 1$$

$$\left|\frac{d-cos\theta cos{\theta }_i}{sin\theta sin{\theta }_i}\right|\le 1$$

(11)

Denote by ${\theta }_{lim}^{I,II}$ the limiting values of $\theta$, i.e., the values corresponding to the two circles of latitude osculating the primary circle:

$$\frac{d-cos{\theta }_{lim}^{I,II}cos{\theta }_i}{sin{\theta }_{lim}^{I,II}sin{\theta }_i}=\pm 1$$

$$d-cos{\theta }_{lim}^{I,II}cos{\theta }_i=\pm sin{\theta }_{lim}^{I,II}sin{\theta }_i$$

$$cos({\theta }_{lim}^{I,II}\mp {\theta }_i)=d$$

$${\theta }_{lim}^{I,II}=2\pi n\pm {\cos }^{-1}\left(d\right)\pm {\theta }_i$$

(12)

We choose n such that ${\theta }_{lim}^{I,II}\in [0,\pi )$. Therefore, given a primary point and a distance d, the $\theta$ coordinate of any point on the secondary circle lies between ${\theta }_{lim}^I\in [0,\pi )$ and ${\theta }_{lim}^{II}\in [0,\pi )$. Note that the largest difference between ${\theta }_i$ (of the primary point) and $\theta$ of a secondary point is obtained for any of the limiting values $\theta ={\theta }_{lim}^{I,II}$. The difference is the angle between the radius vectors pointing to the primary point and to either of the two points of osculation, and this angle is equal to ${cos}^{-1}d$.

3.3. Sphere of Radius R

The result expressed by Equations (8) and (10) can be easily generalized from the unit sphere to a sphere of radius R. The vector representing a point on the sphere (Equation (5)) is multiplied by R. Substituting in Equation (4), we obtain

$$R(cos\theta cos{\theta }_i+sin\theta sin{\theta }_{i}cos(\phi -{\phi }_i))=d$$

where, generalizing Equation (2), d is now related to the radius of the circle by $d=\sqrt{R^2-r^2}$. Thus, Equations (6) and (8) will hold if d is replaced by $d_R$, where

$$d_R=\frac{d}{R}=\frac{\sqrt{R^2-r^2}}{R}=\sqrt{1-{\left(\frac{r}{R}\right)}^2}$$

So, technically, an algorithm to search for circles of radius r on a unit sphere can be readily applied to spheres of radius R by normalizing the circle radius r by R.

4. Algorithm for Great-Circle Detection

In this section, we present an algorithm for finding great circles on a sphere. More accurately, given a set $S=\left\{({\theta }_i,{\phi }_i)\right\}$ of N co-spherical points, we find the great circle that passes through the largest subset of S. Specifically, the (secondary) great circle passing through the largest possible subset of (primary) points is defined by the (secondary) point $({\phi }_{max},{\theta }_{max})$. The algorithm belongs to the Hough transform family, since it is based on voting in a quantized parameter space.

The parameter plane is $(\phi ,\theta )$ is limited to the rectangle $\phi \in [0,\pi )$ $\theta \in [0,\pi )$ and tessellated by rectangular cells. Note that, in the case of great circles, Equation (10) specifies $\theta$ as a function of $\phi$, hence the parameterization is expressed as $(\phi ,\theta )$. In the general case, where Equation (8) determines $\phi$ as a function of $\theta$, the parametrization will be written as $(\theta ,\phi )$. Each cell represents a subset of (secondary) spherical points, and is of the form $\phi \in [\tilde{\phi }-\frac{\Delta \phi }{2},\tilde{\phi }+\frac{\Delta \phi }{2}]$ and $\theta \in [\tilde{\theta }-\frac{\Delta \theta }{2},\tilde{\theta }+\frac{\Delta \theta }{2}]$, where $(\tilde{\phi },\tilde{\theta })$ is the center of the cell. In common Hough algorithms, the cells are of uniform shape and size, their centers form a set of rectangular grid points, and each cell is represented by an accumulator $M(i,j)$ in a memory array, referred to as the accumulator array.

For each (primary) point in the data-set, we compute $\theta$ as a function of $\phi$ (Equation (10)) and increment the accumulators $M(i,j)$ corresponding to cells along the $\theta \left(\phi \right)$ curve. The accumulation stage of the algorithm can be visualized as drawing the discretized curves $\theta \left(\phi \right)$ in the $(\phi ,\theta )$ parameter plane for all primary points. Eventually, the indices $(i_{max},j_{max})$ of the accumulator with the maximal reading represent the parameters $({\phi }_{max},{\theta }_{max})$ of the great circle. Rather than determining all the cells that the curve $\theta \left(\phi \right)$ passes through, the common Hough transform practice is to approximate the process by sampling the argument ($\phi$) and quantizing the function ($\theta$). This leads in our case to Algorithm 1:

Algorithm 1: great-circle detection

Initialize the matrix (accumulator array) M with zeros. For each primary point 1 to N         For each $\phi \in [0,\pi )$ step $\Delta \phi$                Calculate $\theta \left(\phi \right)$ by Equation (10).                Round each $\theta \left(\phi \right)$ to the nearest integer multiple of $\Delta \theta$.                Increment the accumulator of M corresponding to $(\phi ,\theta )$. Find ${\phi }_{max}$ and ${\theta }_{max}$ corresponding to the entry with the largest accumulation in M.

The computational complexity of the algorithm is linear in the number of primary points N. With rectangular tessellation of the parameter plane, the computational load also depends on the number of discrete $\phi$ values, $N_{\phi }=\pi /\Delta \phi$, and on the number of discrete $\theta$ values, $N_{\theta }=\pi /\Delta \theta$. Specifically, the computational complexity of the voting stage is $O(N·N_{\phi })$ and that of the search stage is $O(N_{\phi }·N_{\theta })$.

Note that the approximation, sampling $\theta$ and quantizing $\phi$ rather than determining the cells that the curve $\theta \left(\phi \right)$ passes through, leads to a pitfall. Since the magnitude of the slope of $\theta \left(\phi \right)$ can be large, the discretized $\theta \left(\phi \right)$ curves may not be digitally connected [20]. Thus, they may not intersect in the digital domain even though they do intersect in the continuous domain. This hazard may not be apparent in cases where N is large and intersections are dense, but may manifest itself in cases where N is small. This may then lead to peak-spreading in the parameter plane, and consequently to misdetection and poor localization of peaks. In contrast, if the approximation is discarded and the cells intersected by the curve are incremented, the discretization of the curve amounts to the classical square quantization scheme [21], and the discretized curve is digitally connected.

5. Algorithm for General (Radius r) Circle Detection on a Sphere

We now consider the more general problem—detecting circles of any given radius r on the sphere. Here, the given circle radius may not coincide with the radius of the sphere; hence, the circles to be detected are not necessarily great circles. Formally, given a set $S=\left\{({\theta }_i,{\phi }_i)\right\}$ of co-spherical points and a circle radius r, we find the circle of radius r passing through the largest subset of S. In this case, the parameter plane is limited to a larger rectangle $\theta \in [0,\pi )$ $\phi \in [0,2\pi )$. We compute ${\phi }^I$ and ${\phi }^{II}$ as a function of $\theta$ for each primary point using Equation (8), and increment the accumulators corresponding to cells intersected by the ${\phi }^{I,II}\left(\theta \right)$ curves. Approximating the process by sampling the argument $\theta$ and quantizing the functions ${\phi }^{I,II}$ leads to Algorithm 2:

Algorithm 2: General (any given radius r) circle detection on a sphere

Initialize the matrix (accumulator array) M with zeros. Normalize the radius r by the sphere radius R. Calculate d by Equation (2). For each primary point 1 to N       For each $\theta \in [0,\pi )$ step $\Delta \theta$             If $\theta$ satisfies Equation (11) then                   Calculate ${\phi }^I\left(\theta \right)$ and ${\phi }^{II}\left(\theta \right)$ by Equation (8).                   Round ${\phi }^{I,II}\left(\theta \right)\in [0,2\pi )$ to the nearest integer multiples of $\Delta \phi$.                   Increment the accumulators of M corresponding to $(\theta ,{\phi }^{I,II})$. Find ${\theta }_{max}$ and ${\phi }_{max}$ corresponding to the entry with the largest accumulation in M.

Here the computational complexity of the voting stage is $O(N·N_{\theta })$ and that of the search stage is again $O(N_{\theta }·N_{\phi })$, where in this case $N_{\theta }=2\pi /\Delta \theta$. The accumulation stage of Algorithm 2 can be visualized as drawing the discretized ${\phi }^{I,II}\left(\theta \right)$ curves in the $(\theta ,\phi )$ parameter plane for all primary points. As discussed in Section 4, care should again be taken when considering the use of the voting approximation (sampling the argument $\theta$ and quantizing the functions ${\phi }^{I,II}$ rather than determining the cells that the curves pass through).

6. Parameter-Space Quantization and Error Analysis

6.1. Preliminaries and Motivation

In Hough transform algorithms [2], a data (image) point is mapped into a space of parameters for a class of geometric primitives (such as straight lines). The parameter space is quantized into cells and represented by an array of accumulators. The mapping is carried out by incrementation of the relevant accumulators (voting). Once all data points have voted, a search for the maximum in the accumulator array reveals the parameters of the geometric primitive instance that are best supported by the data.

Parameter space quantization is an essential element of the Hough transform. On the one hand, it determines the resolution of the output parameters. On the other hand, it reflects the tolerance of the algorithm to location errors in the data. Thus, coarse quantization leads to poor output resolution, but quantization that is too fine may lead, in the presence of data location errors, to peak-spreading in the accumulator array and consequent detection errors.

In the Hough transform for straight lines using the normal parameterization, uniform quantization of the $(\rho ,\theta )$ parameter space into cells of equal size implies an equal (translation- and rotation-invariant) measure of the infinite set of straight lines represented by each cell [2,22,23]. In the proposed algorithms for great- and small-circle detection on a sphere, the parameter space is (in principle) an isomorphic spherical parameter space. This is analogous to the isomorphism between the image space and the parameter space in the Hough transform for circles of known radius [2]. Note that the area of a cell in the spherical parameter space is, due to rotational symmetry, a rotation-invariant measure of the infinite set of the spherical circles represented by the cell. We therefore wish to quantizate the spherical parameter space into cells of similar area.

6.2. Towards Uniform Spherical Quantization

Working with spherical memory arrays is currently impractical. We therefore carry out the computation using a planar memory array, where the orthogonal planar axes are $\theta$ and $\phi$. Nevertheless, mapping the conceptual spherical array to a planar array is problematic [19]. Specifically, in the context of parameter-space quantization, uniform (rectangular) quantization in the plane, via uniform quantization of the individual angular coordinates $\theta$ and $\phi$, corresponds to cells of non-uniform area on the sphere. Assuming small quantization steps $\Delta \theta$ and $\Delta \phi$ in the plane, the spherical surface area of the corresponding cell at $\theta$ and $\phi$ is approximately

$$\Delta S\approx \Delta \theta \Delta \phi sin\theta$$

(13)

where $\Delta S$ is largest at the equator, where $\theta =\pi /2$, and equal to $\Delta \theta \Delta \phi$. For non unit spheres, $\Delta S$ is scaled by $R^2$ where R is the radius.

Cells of non-uniform area imply biased voting. One way to handle bias is by normalizing each accumulator count by the area of the corresponding cell. Note, however, that the normalization is associated with singularity at the poles. Lutton et al. [13] suggested using almost-equal-sized cell quantization, letting $\Delta \phi$ depend on $\theta$; see Figure 7. In our case, almost-equal-sized cell quantization is obtained by letting

$$\Delta {\phi }^n=\frac{\Delta S}{\Delta \theta sin{\theta }_n}$$

(14)

where the layer number n is obtained by rounding $\theta /\Delta \theta$ to the nearest integer, ${\theta }_n=n\Delta \theta$, and $\Delta S$ is the almost uniform cell size (spherical surface area).

When voting into the (non-uniform) planar array corresponding to almost-equal-sized cells on the sphere, we again need to decide whether to increment all cells through which the continuous voting curve passes, or to employ the approximation—sampling $\theta$ and quantizing $\phi$ in the context of general circles. The two options are illustrated in Figure 8. Incrementing all cells through which the continuous voting curve passes means voting for both the blue and green cells. In contrast, sampling $\theta$ and quantizing $\phi$ generates (the centers of) the green cells alone.

In our analysis and experiments, the first option is followed, i.e., the approximation is not employed and all cells through which the continuous voting curve passes are incremented. Technically, referring to Figure 9, at layer ${\theta }_n$ we increment the accumulator(s) between A and B inclusively.

6.3. The Distance of Primary Points from the Common Secondary Circle

Consider the cell centered at $({\theta }_{max},{\phi }_{max})$ that received the largest number of votes from primary points. In other words, the cell centered at $({\theta }_{max},{\phi }_{max})$ is the one intersected by the largest number of primary circles centered at primary points. Due to primary points’ location errors, and due to the finite (i.e., not infinitesimal) cell size, the primary circles intersecting the cell do not usually intersect at a single point within the cell themselves. Thus, typically, the primary points that have voted for the cell are only roughly co-circular. The common secondary circle passing near these primary points is defined by the approximate secondary point $({\theta }_{max},{\phi }_{max})$, i.e., by the center of the cell.

What is the maximum distance of the primary points contributing to the cell from the secondary circle defined by $({\theta }_{max},{\phi }_{max})$? Furthermore, how does this depend on the finite cell size? This distance determines, on the one hand, the resolution of the algorithm and, on the other hand, its tolerance to noise (errors) in the location of the data (primary) points. To answer these questions, we calculate the distance between the common secondary circle (the black circle in Figure 10) and the circles defined by the four corners of the cell. The coordinates of the four corners are

$$\begin{array}{c}\hfill ({\theta }_{max}+\frac{\Delta \theta }{2},{\phi }_{max}+\frac{\Delta {\phi }^k}{2})\\ \hfill ({\theta }_{max}+\frac{\Delta \theta }{2},{\phi }_{max}-\frac{\Delta {\phi }^k}{2})\\ \hfill ({\theta }_{max}-\frac{\Delta \theta }{2},{\phi }_{max}+\frac{\Delta {\phi }^k}{2})\\ \hfill ({\theta }_{max}-\frac{\Delta \theta }{2},{\phi }_{max}-\frac{\Delta {\phi }^k}{2})\end{array}$$

where $\Delta {\phi }^k$ is the layer-dependent quantization step (Equation (14)) at ${\theta }_k={\theta }_{max}$, and the corresponding circles are shown blue in Figure 10.

Given one of the (corner) circles defined by a corner point (the blue circle in Figure 11), the largest distance from the common secondary circle is achieved by two critical points (red points in Figure 11) on the corner circle. One critical point is closer than the common secondary circle to the cell center, whereas the other critical point is more distant than the common secondary circle from the cell center.

Consider any of the two critical points and its respective closest point on the secondary circle. Viewing the critical point and its closest point as (spherical) radius vectors, let $\gamma$ denote the angle between them. Let a denote the radius vector corresponding to the secondary point $({\theta }_{max},{\phi }_{max})$, and let b denote the radius vector associated with the corner point. Referring to Figure 12 (left), we observe that, since a is normal to the secondary circle and b is normal to the corner circle, $\gamma$ is also the angle between a and b. Assuming a unit sphere,

$$\gamma ={cos}^{-1}(a·b)$$

(15)

The angle $\gamma$ can be viewed as an angular distance measure between the cell center $({\theta }_{max},{\phi }_{max})$ and one of its corners. With almost-equal-sized cell quantization, $\gamma$ is similar between cells and between corner points of a particular cell.

Cells have four corner points. Each corner point is associated with two critical points, one closer to and one more distant from the cell center. The eight critical points can be divided into two subsets, each with four critical points. One subset contains the four critical points, shown red in Figure 13 (left), closer to the cell center, whereas the other subset contains the four critical points more distant from the cell center, shown in red in Figure 13 (right). The angular distance $\gamma$ between any of the eight critical points and its closest point on the common secondary circle is similar. This means that the four critical points in each of the two subsets are nearly co-circular. The two circles, shown in red in Figure 13, each defined by the four critical points in each of the two subsets, are referred to as lower and higher critical circles. They are on planes parallel to each other and also to the secondary circle plane, where the lower critical circle is closer to the origin and the higher is more distant from the origin. Specifically, their distances from the origin are

$$d_{lower}=cos({cos}^{-1}d+\gamma )\hspace{1em}{d}_{higher}=cos({cos}^{-1}d-\gamma )$$

(16)

The spherical radius vector defined by the angular coordinates $({\theta }_{max},{\phi }_{max})$ of the secondary point passes through the centers of the secondary circle and the two critical circles, and is perpendicular to (the planes containing) them; see Figure 12 (right).

Referring to Figure 14, let P be a primary point of which the primary circle (blue point and blue circle) votes for a cell centered at $({\theta }_{max},{\phi }_{max})$. In other words, the primary circle defined by the primary point P intersects the cell $({\theta }_{max},{\phi }_{max})$. This implies that the coordinates of any point on the primary circle within the cell, i.e., any secondary point S in the cell (green point), are bounded between $\theta =({\theta }_{max}-\Delta \theta /2,{\theta }_{max}+\Delta \theta /2]$ and $\phi =({\phi }_{max}-\Delta {\phi }^k/2,{\phi }_{max}+\Delta {\phi }^k/2]$. Consider the common secondary circle defined by the cell center (black circle) and the nearby secondary circle associated with any other secondary point within the cell (green circle). The nearby secondary circle passes through the primary point P, whereas the common secondary circle generally passes just near P.

The angle between the radius-vector pointing at the cell center and the radius vector pointing at the nearby secondary point is less than or equal to $\gamma$, and equality to $\gamma$ is achieved when the nearby secondary point coincides with one of the corner points of the cell. The radius vectors are normals of the planes containing these circles. So, the angle between the plane containing the common secondary circle and the plane containing the nearby secondary circle is less than $\gamma$. The critical circles are defined such that the angle between the radius vector pointing to a point on the critical circle and the radius vector pointing to its closest point on the common secondary circle is $\gamma$. So, the nearby secondary circle lies on a spherical segment (see Figure 15) containing the common secondary circle and delimited by the two critical circles (red circles in Figure 14 and Figure 15). In other words, the angle between the radius vector pointing to the primary point P and the radius vector pointing to its closest point on the common secondary circle is less than $\gamma$. Hence, the orthodromic distance (great circle distance) of the primary point P to each point on the common secondary circle is upper-bounded by $R·\gamma$, where R is the radius of the sphere.

As shown in Appendix A, $\gamma$ can be readily expressed in terms of the almost-equal-sized cellquantization parameters, $\Delta S$ and $\Delta \theta$. For a unit sphere

$$\gamma \approx \frac{1}{2}\sqrt{{(\Delta \theta )}^2+{\left(\frac{\Delta S}{\Delta \theta }\right)}^2}$$

(17)

In this case $\gamma$ is approximately half of the diagonal of the almost-equal-sized cell. For a sphere of radius R, $\Delta S$ in Equation (17) should be scaled down by $R^2$.

To conclude, the primary point of which the primary circle votes for a cell centered at $({\theta }_{max},{\phi }_{max})$ is located on the spherical segment defined by the cell center. The orthodromic distance between the primary point and the common secondary circle (defined by the cell center) is upper-bounded by $R·\gamma$, where R is the sphere radius and $\gamma$ is the central angle between the radius vector $({\theta }_{max},{\phi }_{max})$ and one of the radius vectors $({\theta }_{max}\pm \Delta \theta /2,{\phi }_{max}\pm \Delta {\phi }^k/2)$ pointing at the cell corners, expressible in terms of the almost-equal-sized cell area.

6.4. Improved Algorithm for General (Radius r) Circle Detection on a Sphere

Consider the problem of general (any given radius r) circle detection on a sphere, as presented in Section 5. We now apply the almost-equal-sized cell quantization scheme, and vote for all cells intersected by the voting curve. For each primary point in the dataset, the algorithm first calculates the limiting values of $\theta$, and loops over all multiples of $\Delta \theta$ between the limiting values ${\theta }_{lim}^{I,II}$. For each $\theta$, it computes the $\phi \in [0,2\pi )$ coordinates of the points A and B by substituting $\theta \pm \Delta \theta /2$ in Equation (8). It then increments all the accumulators which correspond to the discrete points with $\theta$ and $\phi \in [{\phi }_A$, ${\phi }_B]$. Finally, the accumulator $({\theta }_{max},{\phi }_{max})$ with the maximal reading indicates the spherical circle of radius r passing through the largest subset of (primary) points in the data-set. This leads to Algorithm 3:

Note that the accumulator array M in this algorithm is no longer a rectangular matrix. The azimuthal step between adjacent accumulators in the ${\theta }_i$ layer is $\theta$-dependent and equals $\Delta {\phi }^i$. Therefore, the number of accumulators in each ${\theta }_i$ is different. Figure 16 (left) shows the accumulator array in the $\theta -\phi$ plane. In Figure 16 (right), the accumulators are represented by identical rectangles, where the i axis represents the nearest integer multiple of $\Delta \theta$ and the j axis represents the nearest integer multiple of $\Delta {\phi }^i$.

Algorithm 3: Improved general circle detection on a sphere

Initialize the accumulator array M with zeros. Normalize the radius r by the sphere radius R. Calculate d by Equation (2). For each primary point 1 to N       Calculate ${\theta }_{lim}^{I,II}\in [0,\pi )$ by Equation (12).       For each $\theta =i·\Delta \theta$ $(i\in \mathbb{Z})$ from $min\left\{{\theta }_{lim}^{I,II}\right\}$ until $max\left\{{\theta }_{lim}^{I,II}\right\}$:             For both ${\phi }^I$ and ${\phi }^{II}$ in Equation (8):                   1. Calculate $\Delta {\phi }^i$ by Equation (14).                   2. Calculate ${\phi }_A^I$ and ${\phi }_B^{II}$ by substituting $\theta \pm \frac{\Delta \theta }{2}$ in Equation (8).                   3. Choose n such that ${\phi }_A={\phi }_A^{I,II}\in [0,2\pi )$ and ${\phi }_B={\phi }_B^{I,II}\in [0,2\pi )$                   4. Round ${\phi }_A$ and ${\phi }_B$ to the nearest integer multiples of $\Delta {\phi }^i$.                   5. Increment the accumulators from to $({\phi }_A,\theta )$ to $({\phi }_B,\theta )$ Find ${\theta }_{max}$ and ${\phi }_{max}$ corresponding to the entry with the largest accumulation in M.

7. Great Circle Examples

What is the great circle on planet Earth that passes through the largest number of airports? Furthermore, what is the great circle that passes through the largest number of cities? We answer these questions using the proposed algorithms, relying on city [24] and airport [25] databases.

7.1. Uniform Quantization in the Plane (Algorithm 1) Examples

(1)

We execute Algorithm 1, taking the world’s major airports (red dots in Figure 17) as input primary points (a total of 617 points). The parameter space is discretized with $\Delta \phi =\Delta \theta ={0.5}^{\circ }$. The great circle that passes through the maximum number of airports (blue circle in Figure 17) is on the great circle plane with normal ${\theta }_{max}={55.25}^{\circ },\hspace{0.17em}{\phi }_{max}={151.25}^{\circ }$ and visits 29 airports, listed in Appendix B. Note that here and in the following, the number of digits after the decimal point in ${\theta }_{max}$ and ${\phi }_{max}$ relates to the location of the cell center. Obviously, the accuracy of the detected circle parameters follows from the cell size, i.e., from $\Delta \theta$ and $\Delta \phi$. In order to find the circle passing through the largest number of remaining airports (excluding the 29 airports associated with the first maximum), we identify all primary points of which the voting contributed to the first maximum and remove their votes (unvote [5,6,7]). The new maximum represents the great circle passing through the largest subset of remaining (primary) points in the data-set. It lies on the great circle plane with normal ${\theta }_{ma{x}_2}={54.75}^{\circ },\hspace{0.17em}{\phi }_{ma{x}_2}={157.75}^{\circ }$ and visits 28 different airports.

(2)

We again execute Algorithm 1, this time taking 12,960 world cities (red dots in Figure 18) as input primary points. The parameter space is discretized with $\Delta \phi =\Delta \theta ={0.25}^{\circ }$. The great circle (blue circle in Figure 18) that passes through the maximum number of cities is on the great circle plane with normal ${\theta }_{max}={53.625}^{\circ },{\phi }_{max}={156.375}^{\circ }$ and visits 397 different cities. The great circle that passes through the maximum number of remaining cities (i.e., the maximum after unvoting) is on the great circle plane with normal ${\theta }_{ma{x}_2}={72.875}^{\circ },{\phi }_{ma{x}_2}={18.375}^{\circ }$ and visits 321 cities.

(3)

We execute the algorithm in order to find the great circle that passes through a specific airport and the maximum number of other airports. Instead of searching for the maximum value in all the accumulators, we search only in accumulators along the $\theta \left(\phi \right)$ curve corresponding to the specific airport. This reduces the computation time, since we search in $N_{\phi }$ instead of $N_{\phi }·N_{\theta }$ accumulators.

Taking Cape Town International Airport as the specific airport, the great circle that passes through Cape Town International Airport and through the maximum number of other airports is defined by the normal direction ${\theta }_{max}={84.25}^{\circ },{\phi }_{max}={104.75}^{\circ }$ (blue circle in Figure 19) and visits 16 different airports, listed in Appendix C. In this case, the discretization of the parameter space is $\Delta \phi =\Delta \theta ={0.5}^{\circ }$. The great circle that passes through Cape Town International Airport and through the maximum number of remaining airports (the maximum after unvoting, i.e., excluding the airports associated with the first maximum) is on the great circle plane with normal ${\theta }_{ma{x}_2}={37.75}^{\circ },{\phi }_{ma{x}_2}=48.{75}^{\circ }$ and visits 13 different airports.

7.2. Almost-Equal-Sized Cell Quantization Examples

We now address the same great-circle detection tasks using almost-equal-sized cell quantization. To achieve this, we apply Algorithm 3 (general circle detection) with $r=R$ hence $d=0$.

(1)

We execute Algorithm 3 with $d=0$, taking the world’s major airports as primary input points. We set $\Delta \theta ={0.5}^{\circ }$ and $\Delta S=\frac{\pi }{360}·\frac{\pi }{360}=\frac{{\pi }^2}{\mathrm{129,600}}\left[{\mathrm{Rad}}^2\right]$, which yields accuracy similar to the worst-case accuracy (at $\theta =\pi /2)$ in the above uniform quantization example. The great circle (the blue circle in Figure 20) that passes through the maximum number of airports is on the great circle plane with normal ${\theta }_{max}={53.75}^{\circ }$, ${\phi }_{max}={156.4544}^{\circ }$. It visits 30 different airports, listed in Appendix D. The great circle that passes through the largest number of remaining airports (the maximum after unvoting) is on the great circle plane with ${\theta }_{ma{x}_2}={124.25}^{\circ },{\phi }_{ma{x}_2}={331.8655}^{\circ }$ and visits 27 airports.

(2)

We execute Algorithm 3 with $d=0$, taking 12,960 world cities (red dots in Figure 18) as input primary points. We set $\Delta \theta ={0.25}^{\circ }$ and $\Delta S=\frac{\pi }{720}·\frac{\pi }{720}=\frac{{\pi }^2}{\mathrm{518,400}}\left[{\mathrm{Rad}}^2\right]$. The great circle (the blue circle in Figure 21) that passes through the maximum number of cities is on the great circle plane with normal ${\theta }_{max}={53.625}^{\circ },{\phi }_{max}={156.3934}^{\circ }$ and visits 425 cities. The great circle that passes through the maximum number of remaining cities (the maximum after unvoting) is on the great circle plane with normal ${\theta }_{ma{x}_2}={72.375}^{\circ },{\phi }_{ma{x}_2}={18.7609}^{\circ }$ and visits 335 cities.

(3)

We execute Algorithm 3 with $d=0$ in order to find the great circle that passes through a specific airport and through the maximum number of other airports. Instead of searching for the maximum value in all the accumulators, we search only in accumulators along the $\theta \left(\phi \right)$ curve corresponding to the specific airport. We again take Cape Town International Airport as the specific airport, now with $\Delta \theta ={0.5}^{\circ }$ and $\Delta S=\frac{\pi }{360}·\frac{\pi }{360}=\frac{{\pi }^2}{\mathrm{129,600}}\left[Ra{d}^2\right]$. The great circle (the blue circle in Figure 22) that passes through Cape Town International Airport and through the maximum number of other airports is defined by the normal direction ${\theta }_{max}={37.75}^{\circ },{\phi }_{max}={49.3878}^{\circ }$, and visits 15 different airports listed in Appendix E. The great circle that passes through the Cape Town International Airport and through the maximum number of remaining airports (the maximum after unvoting, i.e., excluding the airports associated with the first maximum) is on the great circle plane with normal ${\theta }_{ma{x}_2}={84.25}^{\circ },{\phi }_{ma{x}_2}={104.8324}^{\circ }$ and visits 14 different airports.

The number of data points (airports or cities) associated with the great circle detected depends on the size of the parameter-space quantization cells, i.e., on the quantization level. For the examples in this subsection (almost-equal-sized cell quantization), the dependence is presented in Appendix F.

The differences between the results obtained with uniform quantization in the plane and the those achieved using almost-equal-sized cell quantization are due to the non-equal cell sizes on the sphere yielded by uniform quantization in the plane. Referring to Figure 23, with uniform quantization in the plane (Figure 23 (left)), cells that are close to the equator are larger than cells close to the poles. In contrast, with almost-equal-sized cell quantization (Figure 23 (right)), all cells have approximately the same size. In the above experiments, we chose $\Delta S$ such that the smallest $\Delta {\phi }_k$ in the almost-equal-sized quantization scheme was equal to $\Delta \phi$ used in uniform quantization in the plane. Thus, with almost-equal-sized quantization, the size of the cell on the sphere is larger or equal to that obtained using uniform quantization in the plane. In other words, given cell coordinates ${\theta }_{max}$ and ${\phi }_{max}$ on the sphere, the cell at these coordinates using almost-equal-sized cell quantization is larger or equal in area to the cell at the same coordinates using uniform quantization in the plane. Thus, although the great circle on planet Earth that passes through the largest number of cities is about ${\theta }_{max}={53.625}^{\circ },{\phi }_{max}={156.125}^{\circ }$ using both quantization schemes, with uniform quantization in the plane it is associated with 397 cities, whereas with almost-equal-sized cell quantization the number of cities is 425.

8. General (Small) Circle Examples

In the previous section we found great circles on planet Earth passing through the largest possible number of airports or cities. In this section, we detect small circles of a given radius r that visit as many airports or cities as possible. We use algorithm 3 (almost-equal-sized cell quantization). We normalize the radius of the sphere (planet earth) to $R=1$; hence, $0<r<1$.

(1)

We execute algorithm 3 with $r=0.25$ $(d\approx 0.968)$, taking the world’s largest airports as input primary points. We set $\Delta \theta ={0.5}^{\circ }$ and $\Delta S=\frac{\pi }{360}·\frac{\pi }{360}=\frac{{\pi }^2}{\mathrm{129,600}}\left[{\mathrm{Rad}}^2\right]$. The circle of radius $r=0.25$ that passes through the maximum number of airports (the blue circle in Figure 24) is on the circle plane with normal ${\theta }_{max}={47.75}^{\circ },{\phi }_{max}={256.998}^{\circ }$ and visits 20 large airports. The circle of radius $r=0.25$ that passes through the largest number of remaining airports (red circle in Figure 24, corresponding to the maximum after unvoting) is on the circle plane with normal ${\theta }_{ma{x}_2}={40.75}^{\circ }$, ${\phi }_{ma{x}_2}={269.2340}^{\circ }$ and visits 19 different airports.

(2)

We now execute algorithm 3 with $r=0.5$ $(d\approx 0.866)$, again taking the world’s largest airports as input primary points. We set $\Delta \theta ={0.5}^{\circ }$ and $\Delta S=\frac{\pi }{360}·\frac{\pi }{360}=\frac{{\pi }^2}{\mathrm{129,600}}\left[{\mathrm{Rad}}^2\right]$. The circle of radius $r=0.5$ that passes through the maximum number of airports (the blue circle in Figure 25) is on the circle plane with normal ${\theta }_{max}={27.25}^{\circ }$, ${\phi }_{max}={313.6364}^{\circ }$ and visits 20 large airports. The circle of radius $r=0.5$ that passes through the largest number of remaining airports (red circle in Figure 25, corresponding to the maximum after unvoting) is on the circle plane with normal ${\theta }_{ma{x}_2}={45.25}^{\circ },{\phi }_{ma{x}_2}={313.8552}^{\circ }$ and visits 19 different airports.

(3)

We execute algorithm 3 with $r=0.125$ $(d\approx 0.992)$, this time taking the input primary points from the city database. We set $\Delta \theta ={0.5}^{\circ }$ and $\Delta S=\frac{\pi }{360}·\frac{\pi }{360}=\frac{{\pi }^2}{\mathrm{129,600}}\left[{\mathrm{Rad}}^2\right]$. The circle of radius $r=0.125$ that passes through the maximum number of cities (the blue circle in Figure 26) is on the circle plane with normal ${\theta }_{max}={43.7500}^{\circ }$, ${\phi }_{max}={280.1205}^{\circ }$ and visits 616 cities. The circle of radius $r=0.125$ that passes through the largest number of remaining cities (the black circle in Figure 26, corresponding to the maximum after unvoting) is on the circle plane with normal ${\theta }_{ma{x}_2}={52.7500}^{\circ },{\phi }_{ma{x}_2}={270.4712}^{\circ }$ and visits 386 cities.

(4)

We execute algorithm 3, now with $r=0.5$ $(d\approx 0.866)$, taking the input primary points from the city database. We set $\Delta \theta ={0.5}^{\circ }$ and $\Delta S=\frac{\pi }{360}·\frac{\pi }{360}=\frac{{\pi }^2}{\mathrm{129,600}}\left[{\mathrm{Rad}}^2\right]$. The circle of radius $r=0.5$ that passes through the maximum number of cities (the blue circle in Figure 27) is on the circle plane with normal ${\theta }_{max}={26.7500}^{\circ }$, ${\phi }_{max}={251.6667}^{\circ }$ and visits 568 cities. The circle of radius $r=0.5$ that passes through the largest number of remaining cities (the black circle in Figure 26, corresponding to the maximum after unvoting) is on the circle plane with normal ${\theta }_{ma{x}_2}={72.2500}^{\circ },{\phi }_{ma{x}_2}={251.6327}^{\circ }$ and visits 523 cities.

The experimental results in Section 7 and Section 8 were obtained using Matlab on a low-end personal computer with i5-8265U 1.6GHz CPU and 8.00 GB RAM. The Matlab computing time for the various experiments is summarized in

Table 1. Total computing time and computing time per-point (msec) using Matlab on a low-end PC. Results are reported for great-circle detection ($r=1$) and small-circle detection ($r<1$) and for different parameter-space resolution values ($\Delta \theta$ and $\Delta \phi$ for Algorithm 1, $\Delta \theta$ and $\Delta S$ for Algorithm 3).

Dataset	Algorithm 1				Algorithm 3
	Δθ	Δφ	Computing Time (ms)		Δθ	ΔS	Computing Time (ms)
			Total Time	Avg. per Point			Total Time	Avg. per Point
Airports	2	2	366	0.43	2	${\left(\frac{\pi }{90}\right)}^2$	595	0.93
$r=1$	1	1	531	0.71	1	${\left(\frac{\pi }{180}\right)}^2$	791	1.33
	0.5	0.5	1298	1.98	0.5	${\left(\frac{\pi }{360}\right)}^2$	1777	3.02
(616 pts.)	0.25	0.25	3321	5.16	0.25	${\left(\frac{\pi }{720}\right)}^2$	3668	5.96
Cities	2	2	3644	0.29	2	${\left(\frac{\pi }{90}\right)}^2$	6755	0.53
$r=1$	1	1	6174	0.47	1	${\left(\frac{\pi }{180}\right)}^2$	10,225	0.75
	0.5	0.5	22,779	1.64	0.5	${\left(\frac{\pi }{360}\right)}^2$	26,231	2.07
(12,959 pts.)	0.25	0.25	64,385	5.02	0.25	${\left(\frac{\pi }{720}\right)}^2$	72,061	5.56
Cities	-	-	-	-	2	${\left(\frac{\pi }{90}\right)}^2$	5972	0.48
$r=0.5$	-	-	-	-	1	${\left(\frac{\pi }{180}\right)}^2$	8343	0.70
	-	-	-	-	0.5	${\left(\frac{\pi }{360}\right)}^2$	23,403	1.98
(12,959 pts.)	-	-	-	-	0.25	${\left(\frac{\pi }{720}\right)}^2$	64,078	4.97
Cities	-	-	-	-	2	${\left(\frac{\pi }{90}\right)}^2$	5562	0.43
$r=0.125$	-	-	-	-	1	${\left(\frac{\pi }{180}\right)}^2$	7630	0.57
	-	-	-	-	0.5	${\left(\frac{\pi }{360}\right)}^2$	22,423	1.71
(12,959 pts.)	-	-	-	-	0.25	${\left(\frac{\pi }{720}\right)}^2$	58,805	4.53

. As expected, computing time was roughly proportional to the number of data points, and increased as the parameter space quantization was refined.

9. Conclusions

In this paper we present the first comprehensive solution to the problem of detecting circles, great and small, in a spherical set of points. The suggested approach follows the Hough transform methodology and addresses the fundamental issues associated with Hough algorithms, including parameterization, parameter-space quantization, and analysis of the quantization-induced error. The time-complexity of the proposed algorithms is linear in the number of data-points. The required memory size is simply the number of cells in the accumulator array, which is bounded and two-dimensional in the case of great circles and in the case of small circles of known radius.

Potential applications for this research arise in cases where the data are naturally spherical and where circular features need to be detected. The planetary and geophysical sciences provide numerous examples, including crater and volcano detection. Computer vision applications arise in omnidirectional image analysis, including fisheye image analysis, where spherical image representations are ubiquitous.

The proposed method can be readily extended to the detection of circles of unknown radius in a spherical set of points. The additional unknown can be accommodated either by sequentially applying the known-radius algorithm with different radii, or by adding a dimension representing the unknown radius to the accumulator array.