Patch Area and Uniform Sampling on the Surface of Any Ellipsoid

Algorithms for generating uniform random points on a triaxial ellipsoid are non-trivial to verify because of the non-analytical form of the surface area. In this paper, a formula for the surface area of an ellipsoidal patch is derived in the form of a one-dimensional numerical integration problem, where the integrand is expressed using elliptic integrals. In addition, analytical formulae were obtained for the special case of a spheroid. The triaxial ellipsoid formula was used to calculate patch areas to investigate a set of surface sampling algorithms. Particular attention was paid to the eﬃciency of these methods. The results of this investigation show that the most eﬃcient algorithm depends on the required coordinate system. For Cartesian coordinates, the gradient rejection sampling algorithm of Chen and Glotzer is best suited to this task, when paired with Marsaglia’s method for generating points on a unit sphere. For outputs in polar coordinates, it was found that a surface area rejection sampler is preferable.


Introduction
Generating uniform random samples of points on surfaces finds many applications. Some examples include: solving radiation transport problems in medical physics [1], quantifying errors in brain image analysis techniques [2], modelling oxygen production of trees in forests [3], simulating the effect of background radiation on detector surfaces [4], statistical goodness-of-fit testing [5], solving problems in development of visualisation software [6] and testing robot motion planning algorithms [7]. The particular problem of sampling the surface of an ellipsoid has been of interest for modelling: dose-rate distributions of Iodine-125 for radiation therapy [1], prolate virus capsid formation [8], protein coats of bacterial spores [9] and impacts on solar system objects of unusual shape (such as the bilobate Kuiper belt object, Arrokoth) [10].
Many published studies consider sampling of arbitrary surfaces [1,[3][4][5][6]11], as well as the relatively simple example of the sphere [12][13][14][15]. However, there is relatively little published material that considers the ellipsoidal case. Williamson provided a method for the ellipsoid surface based on rejection sampling [1]. Chen and Glotzer improved upon this method, giving a proof and numerical verification for the case of a prolate spheroid [8]. While this method is just as valid for any other kind of ellipsoid, verifying the uniformity of a sample in the triaxial case is a highly non-trivial task because the surface area cannot be expressed as an analytical function.
An explicit expression for the surface area of the entirety of a general ellipsoid was first derived in 1825 by Legendre in terms of incomplete elliptic integrals of the first and second kinds (See Reference [16] for a historical review). Many studies have been undertaken regarding exact and approximate expressions for the surface area of the entire ellipsoid [16][17][18][19][20][21][22][23][24]. However, only one of these considers the problem of finding the area of a subset of the ellipsoid surface [17]. That work involved cutting the ellipsoid into two segments using an intersecting plane and finding the area of those segments.
In this study, the surface area of interest is that of a patch bounded by limits given in (scaled) spherical polar coordinates, (θ, φ), defined such that where θ ∈ [0, π] is an angle measured from the c-axis and φ ∈ [0, 2π) is an azimuth angle defined in the x-y plane. It is assumed in this work that a ≥ b ≥ c. The sole exception to this is a prolate spheroid, for which it is assumed that a = b < c. Figure 1 shows how an ellipsoidal patch is defined using limits θ 0 , θ 1 , φ 0 and φ 1 . The patch is an 'ellipsoidal rectangle' with vertices (θ 0 , φ 0 ), (θ 0 , φ 1 ), (θ 1 , φ 0 ) and (θ 1 , φ 1 ). When a vertex is one of the two poles in θ (with θ = 0 or π), the patch is no longer 'rectangular', but is now the cap of the ellipsoid (with no divisions in φ).
In Section 2, an expression for the surface area of an ellipsoidal patch is derived and a means of numerically evaluating this expression is outlined. Given the surface areas of ellipsoidal patches, uniformity of ellipsoid surface sampling algorithms can be studied. In Section 3, measures of speed as well as statistics describing uniformity are discussed. These measures are applied in Section 4 to a selection of sampling algorithms on a spherical surface. In Section 5, some rejection sampling algorithms for the ellipsoid surface are outlined (some requiring sphere samplers to generate trial points). A similar analysis as for spheres is performed in Section 6 on the surface of an oblate spheroid and a triaxial ellipsoid. This analysis was performed to answer three questions: Firstly, do these algorithms generate uniform distributions on the ellipsoid surface? Answering this requires the patch area formula derived in Section 2. Secondly, Fig. 1 A patch on the surface of an ellipsoid. For limits θ 0 , θ 1 , φ 0 and φ 1 , one can determine curves of constant θ (solid) and of constant φ (dashed). The patch (red) is defined as the area enclosed by the four curves. The dotted vertical line indicates the z-axis (from which θ is measured).
of those samplers found to be uniform, which one is the fastest? Thirdly, does the speed depend on the aspect ratio of the ellipsoid? For example, can the symmetry of a spheroid favour one class of sampler over another? By answering these questions, one can recommend a particular algorithm (or algorithms) to uniformly sample the ellipsoid surface.

Surface Area of an Ellipsoidal Patch
In polar coordinates, the ellipsoid area element is given by [21], so that infinitesimal area dS is, Thus, the area of the ellipsoid is given by, In the last equality, only the octant of the ellipsoid where x ≥ 0, y ≥ 0 and z ≥ 0 is considered. This corresponds to taking 0 ≤ φ ≤ π/2 and 0 ≤ θ ≤ π/2. By the symmetry of the (general) ellipsoid, the area over the entire surface is given by the area of the octant, multiplied by 8. The area of a patch can be obtained by replacing the limits with the values θ 0 , θ 1 , φ 0 and φ 1 , which define the vertices of the patch. The integral of interest here is therefore,

Derivation of the Patch Area Formula
It is entirely possible to integrate the double integral in Equation 7 numerically [25]. However, it is also possible to rewrite the double integral so that the inner integral (over φ) is expressed using incomplete elliptic integrals of the second kind. Following Maas [21], one can make the transformation to give (using that dθ = − sin θ dξ and sin 2 θ = 1 − ξ 2 ), where Note that the minus sign in Equation 10 does not indicate a negative area. This is because cos θ 1 < cos θ 0 , when θ 0 < θ 1 with both limits within [0, π/2]. The minus sign can thus be cancelled by interchanging the ξ limits, and the resulting area is positive. By defining and the patch area may be expressed as, where is Legendre's incomplete elliptic integral of the second kind. If a > b > c, then the value of k(ξ) is purely imaginary for ξ ̸ = 1. This occurs because the factor 1 − (a/b) 2 in the numerator of Equation 12 becomes negative when a > b. The elliptic integral with k purely imaginary can be transformed into an expression involving real valued variables [26], however this involves additional trigonometric and square root evaluations and thus requires extra computations. Instead, one can transform the problem by interchanging a and b in Equation 13 and replacing the φ limits with φ 0 − π/2 and φ 1 − π/2. This rotation gives the equivalent expression for patch area, where, When a > b, Equation 15 should be used for patch area, while Equation 13 gives the appropriate formula for the b > a case. As for the choice of c, note that the value of g(ξ) is real for any positive values of a (or b) and c, as the second term in Equation 11 is always greater than −1 as ξ 2 ∈ [0, 1]. Therefore, in Equations 13 and 15, c can take any positive value (whether smaller or greater than a or b). The elliptic integral E(ψ, k) can be efficiently evaluated by writing it in terms of Carlson's elliptic integrals, with u = cos 2 ψ , where R F and R D are the Carlson elliptic integrals of the first and second kind, respectively. These integrals can be efficiently computed numerically, using algorithms based on duplication theorems [27,28]. If computational speed is a priority, then one could alternatively compute E(ψ, k) using the faster method developed by Fukushima [29,30]. Given an efficient algorithm to compute E(ψ, k), evaluation of Equation 15 can be interpreted as a one-dimensional integration problem. This can be readily solved to the required accuracy using an algorithm such as Gaussian quadrature or Romberg integration [25,31], the latter of which was used in this study.
The patch area formula given in Equation 15 assumes that both the θ and φ limits are within the interval [0, π/2]. Outside of this interval, it is possible to obtain unexpected negative values of patch area (for example, when using Equation 18 for φ > π) that lead to unwanted cancellations. However, the [0, π/2] interval is sufficient to calculate any patch area on the ellipsoid by exploiting its eight-fold symmetry. By this symmetry, the area of a patch entirely contained within the x, y, z ≥ 0 octant of the surface equals the area for reflected versions of the patch in the other octants. If the set of patches have been chosen in such a way that a limit occurs halfway through a patch, then that patch area can be determined by calculating the area of the halfpatch within the limits and doubling the answer. For patches containing both a θ and a φ limit, one need only calculate the quarter of the patch area within the limits and multiply by 4. These cases together account for all patches on the surface, provided that the number of distinct φ values is chosen to be even, so that the set of patches match the symmetry of the ellipsoid (the number of θ values can be odd or even). By exploiting the symmetry of the ellipsoid in this way, many redundant computations are avoided. This can be useful if one wishes to obtain accurate areas by setting a low tolerance value in the chosen numerical integration algorithm.

Spheres and Spheroids
For a sphere of radius a, the patch area integral simplifies to, For a spheroid with a = b ̸ = c, an analytical expression for patch area can be derived. Starting from the triaxial patch area given in Equation 15, it can immediately be seen from Equation 17 that κ(ξ) vanishes when a = b. For κ = 0, the elliptic integral E(ψ, κ) reduces to, Thus, where, To solve this integral, one can make the substitution sinh u = qξ to obtain, with u 0 = arcsinh(q cos θ 1 ) and u 1 = arcsinh(q cos θ 0 ). Using that, and sinh 2u = 2 sinh u cosh u = 2qξ 1 + q 2 ξ 2 , the spheroid patch area is given by, +q cos θ 0 1 + q 2 cos 2 θ 0 − q cos θ 1 1 + q 2 cos 2 θ 1 .
Equation 30 is named as 'oblate' because the quantity q is only real valued when a > c.

Speed and Uniformity of the Sampling Algorithms
Ellipsoid patch areas calculated using Equations 15 (for ellipsoids) and 23 (for spheres) were used to investigate uniformity (with respect to surface area) of the sampling algorithms to be discussed in Sections 4 and 5. For each algorithm studied, a sample of N = 10 8 random points was generated using implementations written in C++, run using the Linux Subsystem for Windows and analysed using a Python script. This was done for two different random number generators (RNG): • Lagged Fibonacci (with 4 feedback shifts and the exclusive-or operation).
The implementations used in this work are from the TRNG library 1 (using the lagfib4xor 19937 64 and yarn5 classes, respectively) [32]. The following analysis was done with two different RNG algorithms to verify that the source of random numbers itself has no influence on the final results. This is done by checking that the results from each generator are similar.
To study each sampling algorithm, the following quantities were measured: • Run-time, t run .
1 See https://www.numbercrunch.de/trng/ For each algorithm (and choice of RNG), the required run-time was measured using the ctime header from the C++ Standard Library.
Since the raw run-time depends on many factors; such as language, operating system and CPU, the relative speed of each algorithm was computed from the measured run-times. The speed relative to a given benchmark algorithm may be computed as t benchmark /t, where t benchmark is the run-time of the benchmark and t is the run-time for the desired algorithm.
For those algorithms that use rejection sampling, the acceptance rate is the proportion of all trial points that were accepted. In cases where rejection sampling is not employed, one obtains a value of r = 1, i.e. all generated points were accepted. The above quantities can be used to evaluate the speed of each sampling method. These were measured by a program that generated the 10 8 points and nothing else. Thus, the measured run-times in this Paper constitute only the time required to generate the points.
To investigate the uniformity of these algorithms, another set of N = 10 8 points was generated for each method. In this case, the points were binned into one of a set of patches (as defined by Figure 1). This was done by defining patches with increments of one degree in both angular coordinates θ and φ. This gives 181 distinct values of θ and 360 values of φ. Since two of the θ values correspond to caps at the poles, there are a total of n = (181 − 2) × 360 + 2 = 64 442 (32) patches (i.e. 64 442 bins) defined on the surface. Given a point in polar coordinates, the relevant θ and φ indices can be found by comparing each coordinate to a list of the values defining each bin. For the sake of speed, these comparisons were here performed using Bottenbruch's version of the binary search algorithm [33,34].
Given a binned distribution of random surface points, one can calculate number of points per area for each patch. From this distribution, a mean and a standard deviation can be computed. The relative standard deviation is then just the standard deviation divided by the mean. This quantity is interpreted here as a measure of non-uniformity, with a large value indicating high deviation from the mean. Since these distributions are expected to be uniform after dividing by patch areas, one can expect the relative standard deviation to be small.
Another way to measure uniformity (or lack thereof) is to perform a χ 2 test for goodness of fit [25,35]. This uses the χ 2 statistic, defined as where index i refers to a particular patch, O i is the observed number of random points within patch i and E i is the expected number of points within patch i. Knowing the area of each patch, one can calculate E i ; given the null hypothesis that this number divided by surface area is uniform. This is given by, where S i is the surface area of patch i and n i S i is the sum of all patch areas. Given χ 2 , the validity of the null hypothesis can be tested by comparing this value to a critical value, χ 2 crit , which is determined using a chi-squared distribution function, with n degrees of freedom. For n = 64 442 degrees of freedom (i.e. equal to the number of bins) the critical value at the 5% confidence level is found to be (using the chi2.ppf function from the Python stats module).
If a measured χ 2 statistic exceeds the critical value, then there is sufficient evidence to reject the null hypothesis and conclude that the distribution is non-uniform. On the other hand, a χ 2 value smaller than critical cannot prove the null hypothesis, but does provide some confidence that the sample could be uniform.

Uniform Sampling on the Surface of a Sphere
Before proceeding to investigate ellipsoid surface samplers, it is necessary to have a reliable method for generating uniform random points on the surface of a sphere. A simple approach, based on a two-dimensional circle rejection sampling method given by von Neumann [36], is to generate points in a unit length cube. Points that have magnitude smaller than unity (i.e. that lie within a unit sphere) are accepted and scaled to lie on the required sphere surface. While this approach suffices, other algorithms exist that are faster. In this work, six sphere sampling algorithms (including the cubic rejection method) were considered. The other five are: Marsaglia's improved method [12], the 'trig method' (so-called due to use of trigonometric functions) [8,15], use of Gaussian random numbers [12,13] (generated here using Doornik's implementation of the ziggurat method [37,38]), Cook's method [12,14] and rejection sampling using the surface area element at a point [3]. For details on how these algorithms work, see the given References. The surface area rejection sampling method used here is a special case of an algorithm for ellipsoids, described in more detail in Section 5.

Comparison of the Sphere Sampling Algorithms
The quantities described in Section 3 were measured for samples of N = 10 8 randomly generated points on the unit sphere. The results of this are shown in Table 1. Considering first the speed (relative to the cubic rejection method) of each algorithm, it is clear that the fastest is that of Marsaglia, which has the shortest run-time and highest relative speed, regardless of random number generator. This is despite the fact that both the trig and Gaussian methods have a perfect acceptance rate. The speed of Marsaglia's method lies in the fact that no costly trigonometric or Gaussian deviate calculations are required. Table 1 Results of using various sphere surface samplers to generate 10 8 points, using two different random number generators. Each column shows (from left to right): algorithm name, run-time (in seconds), speed relative to cubic rejection, acceptance rate, relative standard deviation (as a percentage), the χ 2 statistic and whether χ 2 is smaller than the critical value, χ 2 crit = 65 033.6.
Sphere: (a, b, c) = (1, 1, 1) It is also of note that the relative speed of each algorithm differs based on the choice of random number generator. For lagged Fibonacci, which is faster than YARN5, Marsaglia's method is nearly twice as fast as the trig method. With the YARN5 generator, Marsaglia's method remains fastest, but the gap compared to the trig method is smaller. It is thus conceivable that for a much slower source of random numbers, the trig method could be advantageous.
As far as uniformity is concerned, every sample yields a relative standard deviation of roughly 3.7%, and a χ 2 value just under the critical value (see Equation 35). This may suggest that the distributions are indeed uniform. To shed further light on the uniformity (or non-uniformity) of the obtained samples, one could visually inspect the distributions. Figure 2 shows the distributions obtained using Marsaglia's method with the lagged Fibonacci generator. The raw number of points per patch is greater at the equator than near the poles (with a high value at the patches containing the poles themselves). The patches with more hits are those with greater surface area, as is verified when looking at number of points per patch area. Now a much more uniform distribution presents itself. The use of relative standard deviations and χ 2 tests allow uniformity to be assessed without the need to produce plots for each sample. The corresponding values in Table  1 do not give evidence for non-uniformity, with small standard deviations relative to the mean and all χ 2 tests passing. This is expected since the algorithms are all designed to produce uniformly random points.

Uniform Sampling on the Surface of an Ellipsoid
Suppose one has a point on a unit sphere that was generated using a uniform sampling algorithm, such as one of those discussed in Section 4. This point can easily be transformed so it lies on the surface of an ellipsoid of axis lengths a, b and c by scaling its Cartesian coordinates by each axis respectively. This 'naive method' is simple and only requires a sphere sampler. However, it does not result in a uniform distribution of points on the ellipsoid surface. Instead, there is a greater concentration of points near the poles of the major and middle axes as shown by Figures 3(a) and (b). By using a uniform ellipsoid surface sampler, a sample such as that of Figures 3(c) and (d) can be obtained.
In this section, three different rejection sampling algorithms for generating uniform distributions of points on the ellipsoid surface are explained. These are then studied in Section 6, using the measures discussed in Section 3.

Gradient Vector Rejection Sampling
One approach to generate uniform samples on the ellipsoid surface is to first generate points on the unit sphere, perform anisotropic scaling to obtain trial points, p(x, y, z), on the ellipsoid and then reject some of them in such a way that the accepted points are uniformly distributed. Chen and Glotzer give such a method for prolate spheroids, based on the work of Williamson, using the magnitude of the gradient vector of the ellipsoid surface [1,8]. This approach is easily generalised to the case of any ellipsoid, as is demonstrated here. The rejection sampling works by defining a function giving the probability of accepting a trial point, where the function g(x, y, z) is given by, The vector n is the normal to the surface at point (x, y, z), and ∇f is the gradient vector. The function f defines the surface of the ellipsoid, with the surface being the set of points where f (x, y, z) = 1. The magnitude of the gradient vector is given by, Since the acceptance probability is the ratio of this expression to its maximum value, the prefactor of 2 may be ignored. With this factor dropped, the expression can be used as function g(x, y, z).
To use this expression, the maximum of the function = sin 2 θ cos 2 φ a 2 + sin 2 θ sin 2 φ b 2 + cos 2 θ c 2 , must be found. The latter equality in Equation 40 arises due to conversion of coordinates from Cartesian to scaled spherical polars, using Equations 1-3. Use of polar coordinates gives a more natural parametrisation of the function to be maximised. The partial derivatives of g(θ, φ) can be determined using the chain rule, At the extreme points of g(θ, φ), both partial derivatives equal zero. To make the right-hand side of Equation 41 vanish, either sin θ, cos θ or the term in brackets must be zero. This is equivalent to the condition, θ = 0, or , θ = π/2, or , θ = π, or , sin 2 φ The first three cases simply state that the point is either a c-axis pole or an equator point. The last case is satisfied for any point on a sphere. For Equation 42, one requires either sin φ = 0, cos φ = 0, sin θ = 0 or a = b, i.e. φ = 0, or , φ = π, or , φ = π/2, or , φ = 3π/2, or , θ = 0, or , θ = π, or , a = b .
The latter case says that the φ derivative vanishes at any point on a spheroid, where c is the distinct axis. The remaining cases require particular values of φ or θ. In order to make both partial derivatives vanish simultaneously on an ellipsoid with arbitrary axes, there are only six possible points. These points are the poles of the a, b and c axes. Therefore, the maximum value of g is given by, Substituting Equations 40 and 43 into Equation 36 gives the following expression for the acceptance probability of a trial point, Equation 44 forms the basis of a rejection sampling algorithm using the gradient vector at a trial point generated by the 'naive method'. This procedure is summarised via pseudocode in Algorithm 1.

Area Element Rejection Sampling
A natural approach to generating random points that are uniform with respect to surface area is to make use of the expression for the surface area element, s. Following Melfi and Schoier [3], suppose an area element is given in terms of arbitrary coordinates (µ, ν). Values of µ and ν are obtained using a uniform random number generator. One can then define an acceptance probability, where s max is the maximum value of the surface area element. One could try selecting Cartesian coordinates x and y and use points (x, y, f (x, y)), The surface area element is then given by, The partial derivatives, are not defined at values where z = f (x, y) vanish (i.e. the equator, θ = π/2), due to division by zero. Therefore, Cartesian coordinates are unsuitable for this algorithm. Instead, polar coordinates (θ, φ) can be used. The surface area element is given by Equation 4, and so the acceptance probability is The value of the maximum, derived in Appendix A, depends on the choice of semi-axis lengths. Defining, the maximum is given by, Equations 50 and 52 can be used as the basis of an area rejection sampling method, as shown in Algorithm 2. Values of θ and φ are randomly chosen using uniform generators and the area element used to determine whether to accept the corresponding point.
When dealing with a sphere or a spheroid, the area rejection algorithm becomes computationally simpler. This is because the area element becomes, s spheroid (θ, φ) = sin θ a 2 c 2 sin 2 θ + a 4 cos 2 θ , for a spheroid and s sphere (θ, φ) = a 2 sin θ , for a sphere. In both cases, there is no dependence on φ. This means that the random generation of φ can be taken outside the loop in Algorithm 2. Combined with a simpler form of the area element, this means fewer computations are required.

Algorithm 2
The polar coordinate area rejection algorithm on an ellipsoid with axis lengths a, b and c. θ ⇐ uniform random number on (0, π) 5: φ ⇐ uniform random number on (0, 2π)

8:
if u ≤ P (accept) then 9: accept ⇐ true 10: end if 11: end while 12: x ⇐ a sin θ cos φ 13: y ⇐ b sin θ sin φ 14: z ⇐ c cos θ 15: return (x, y, z) Another area rejection algorithm can be created by selecting another coordinate system. In this work, use of the Mercator parametrisation (u, v), where, was studied. The u coordinate is defined similarly to the polar φ coordinate and exists on the range [0, 2π]. However, the v coordinate is defined on range (−∞, ∞). By truncating to a finite range, it is possible to obtain an approximately uniform sampling algorithm using these coordinates. Here, the finite range (−2π, 2π) was used. Using the coefficients of the first fundamental form, the area element can be found to be, Since this expression requires only one hyperbolic and one trigonometric function evaluation to compute (as cos 2 u = 1 − sin 2 u), a single evaluation of s could be expected to be faster here than for polar coordinates. The partial derivatives are, where, The v derivative vanishes when sech v = 0 or tanh v = 0 (The bracketed term is always positive for a > b > c). The former case requires infinite v, in which limit the area element vanishes. The latter case is equivalent to v = 0, which corresponds to the equator (i.e. z = 0). Given that v = 0, the only way to make the u derivative equal zero is to make u a multiple of 2π (so sin u = 0 or cos u = 0). The maximum possible value of s is thus obtained by taking u = π/2 or 3π/2, so that sin 2 u = 1 and therefore, The acceptance probability in these coordinates is given by, The process for generating random points is similar to that for polar coordinates, as shown by Algorithm 3.

Algorithm 3
The Mercator coordinate area rejection algorithm on an ellipsoid with axis lengths a, b and c. u ⇐ uniform random number on (0, 2π)

4:
v ⇐ uniform random number on (−t, t) (t = 2π used in this study)

Ray Intersection Sampling
The final method considered here is the generic surface sampler of Detwiler et al. [4]. An equivalent formulation of this algorithm also appears in the work of Palais et al. [6]. The method works by starting with a randomly generated point on the bounding sphere for the surface of interest. For an ellipsoid, this sphere has radius r = max(a, b, c). A uniform random point, ρ 0 is then generated in a disk of radius r. Intersections are then sought between the surface and a ray orthogonal to the disk, starting at ρ 0 . If any intersections are found, one of them is randomly selected and returned as the generated surface point. If no intersections exist, the process is re-attempted until a valid intersection is found.
The method is summarised as pseudocode in Algorithm 4. Firstly, the random sphere point, n, is generated. This vector is normal to the sphere surface. Next, the orientation of the tangent disk is represented through two mutually orthogonal vectors in the disk. One of these, u, is constructed as orthogonal to n; while the other, v can be made orthogonal to both n and u by taking a cross product. For subsequent calculations, u and v must be normalised.
A point in the disk must then be generated. This requires a random radius B and a random angle ψ (between 0 and 2π). To ensure uniformly random points in the disk, one must take the square root of a uniform random number on range (0,1) and then multiply by the desired radius [39]. Given point (b 0 , b 1 ) in two-dimensions, the 3D point ρ in the disk is given by, Given point ρ, the line ρ 0 + tρ can be defined, where ρ 0 = rn and t is a real number. With the ellipsoid oriented with semi-axes a, b, c along the x, y, z directions respectively, the problem of determining whether the line intersects the ellipsoid can be represented as a quadratic equation in parameter t, so that where, There are three possible outcomes based on the sign of the discriminant: two intersections; one intersection; no intersection. In the latter case, the calculation is repeated with new random numbers until an intersection is found. If two intersection points are obtained, a random number can be drawn to decide which one to return.

Comparison of the Ellipsoid Sampling Algorithms
Each of the aforementioned algorithms were used to generate a sample of N = 10 8 random points, and were analysed as discussed in Section 3. Eight separate algorithms were used, including different variations of those outlined above. Along with the 'naive scaling method' and the ray intersection method, three variations of the gradient rejection algorithm were studied, ψ ⇐ uniform random number on (0, 2π) 10: q ⇐ uniform random number on (0, 1)
• Converting to polar coordinates prior to output (using Marsaglia's method for spheres).
Three variations of the area rejection method were also used, • Using polar coordinates, converting the output to Cartesians.
The decision to consider implementations of the gradient and area rejection methods with Cartesian and polar outputs was made to study the effect this has on run-time.
If one is interested in generating polar coordinates for uniform random surface points, is it faster to use an algorithm that directly yields polar coordinates? Tables 2 and 3 show the results for samples of N = 10 8 random points on an oblate spheroid, (a, b, c) = (3, 3, 1.5), and a triaxial ellipsoid, (a, b, c) = (3, 2, 1), respectively. First consider the uniformity of each method. In both tables, it is immediately clear Table 2 Results of using various ellipsoid surface samplers to generate 10 8 points on an oblate spheroid, using two different random number generators. Each column shows (from left to right): algorithm name, run-time (in seconds), speed relative to naive rejection, acceptance rate, relative standard deviation (as a percentage of the mean), the χ 2 statistic and whether χ 2 is smaller than the critical value, χ 2 crit = 65 033.6.
Oblate Spheroid: (a, b,  that the relative standard deviation and χ 2 statistics from the naive method are much greater than those for all the other algorithms. This is to be expected, since the naive method yields a non-uniform distribution. Another algorithm that fails the χ 2 tests is the Mercator coordinate area rejection method. The fact that a uniformity test fails might not be considered surprising, since the method uses a finite truncation on the range of the v coordinate, meaning that it is not completely uniform. The loss of uniformity thus occurs near the c-axis poles, with fewer generated points than expected in those regions. For every other algorithm, the χ 2 values are below critical and the relative standard deviations are roughly 3.2%.
Looking at the run-times in Tables 2 and 3, the fastest algorithm is the naive method and the slowest is Mercator coordinate area rejection. The former is fast as it simply scales a point generated on the sphere by Marsaglia's method, while the latter is slow because of a low acceptance rate, combined with the need for trigonometric and hyperbolic function calls. Regardless of efficiency, neither sampler is a good choice due to their non-uniformity.
Of the uniform ellipsoid samplers, the ray intersection algorithm is the slowest. This can be explained in Table 3 by a low acceptance rate. In the oblate case, the value of r is higher and the run-time shorter, but the gradient and area rejection methods are still faster. This is because of the algorithmic complexity of the method, with trigonometric function evaluations and multiple random number generator calls. Table 3 Results of using various ellipsoid surface samplers to generate 10 8 points on a triaxial ellipsoid, using two different random number generators. Each column shows (from left to right): algorithm name, run-time (in seconds), speed relative to naive rejection, acceptance rate, relative standard deviation (as a percentage of the mean), the χ 2 statistic and whether χ 2 is smaller than the critical value, χ 2 crit = 65 033.6. If one wishes to output the generated point in polar coordinates, the area rejection method is faster than the gradient rejection method for both shape examples. The computational cost of converting from Cartesian to polar coordinates in the gradient rejection method is higher than that of using area rejection and returning polar coordinates directly.
When the output is given in Cartesian coordinates, the fastest uniform sampler depends on the choice of random number generator. For the lagged Fibonacci algorithm, gradient rejection is fastest. Meanwhile, area rejection with Cartesian output has slightly greater speed when using the YARN generator. Tables 2 and 3 show that area rejection has a slightly greater acceptance rate than gradient rejection. Therefore, it might be expected that area rejection is more efficient. However, gradient rejection does not require trigonometric function calls and so each iteration for a proposed point is swifter. With a fast enough random number generator, this outweighs the lower acceptance rate to give a more efficient method.
It should also be noted that the acceptance rates given in Tables 2 and 3 apply only to the particular aspect ratios studied in this analysis. The acceptance rate is a function of aspect ratio and might be expected to decrease as the shape deviates from a sphere.

Conclusions
Expressions for the surface area of an ellipsoidal patch (as defined by Figure 1) were derived (see Equations 15,23,30 and 31). The triaxial patch area expression can be evaluated using a one-dimensional numerical integration algorithm. The integrand requires evaluation of elliptic integrals, for which efficient numerical methods exist. This expression was used to investigate ellipsoid sampling algorithms that are uniform with respect to surface area. The three methods studied in Section 5; gradient rejection, area rejection (based on polar coordinates) and ray intersection were all found to result in uniform distributions. For a fast random number generator, the gradient rejection algorithm was found to be the fastest method for generating Cartesian coordinates, with run-time minimised when used with Marsaglia's method for sphere sampling. For outputs in polar coordinates (or for both Cartesian and polar), the area rejection algorithm is more efficient. This was found to occur for both a spheroid and a triaxial ellipsoid. If one wishes to generate polar coordinates and speed is essential, then the area rejection is preferable. Otherwise, the gradient rejection method is recommended due to its mathematical and computational simplicity.

Declarations
Ethical approval. Not Applicable.
• Case 2 : θ = π/2 and φ = π/2 or 3π/2. Each of the above cases shall now be considered in turn, with the extreme value s case no. calculated for each. The maximum of these values is then determined.

Case 3
For a spheroid (a = b), s does not depend on the φ coordinate. Thus, at the equator,
Solving for sin 2 θ gives, and so the area element becomes,

Case 7
For θ = 0 or π, sin θ = 0, making the left-hand side of Equation A5 vanish. This means that A5 cannot be true, since 0 ̸ = −1. Therefore, this case can never hold.
Assuming a triaxial ellipsoid with a > b > c, this leaves three possibilities for the maximum; s 2 , s 4 or s 5 . First compare s 4 and s 5 . On can determine the larger by taking the difference s 4 − s 5 and determining whether it is positive or negative. To simplify the calculation, one can calculate s 2 4 − s 2 5 instead (as area elements are always non-negative). Using the SymPy library in Python [40], the difference yields, The numerator is negative, since a > b. Meanwhile, the denominator is positive since a > c and b > c. Therefore, s 2 4 < s 2 5 when a > b > c. The maximum value is thus either s 5 or s 2 = ac. Taking the difference s 2 5 − s 2 2 (SymPy was again used) gives, The numerator is non-negative and if b > c, the denominator is positive. Therefore, s 5 >= s 2 for a triaxial ellipsoid and so s max = s 5 . For an oblate spheroid with a = b > c, the same argument holds and s 5 gives the maximum. This leaves two special cases to consider; the prolate spheroid with a = b < c, and the sphere. For a prolate spheroid, s 4 = s 5 and the denominator of Equation A16 is negative. Thus s max = s 2 = ac is the maximum. For the sphere, the area element is simply a 2 sin θ, so the maximum is s max = a 2 .
In summary, the ellipsoid (as defined in Section 1) has maximum area element, . (A18)