BACKGROUND
The problem of finding a uniform distribution of points on a sphere has a relatively long history. Its emergence is commonly attributed to the physicist J. J. Thomson, who posed it in 1904 after creating his socalled plum pudding model of the atom [1]. As such, the problem involves determination of a minimum energy configuration of N equally charged particles, confined to the surface of a sphere, that repel each other with a force given by Coulomb's law [1]. Although the plum pudding model of the atom has long been dismissed, the original problem posed by Thomson has reemerged across many areas of study and found practical applications in the fields as diverse as viral morphology, crystallography, physical chemistry, geophysics, acoustics, signal processing, computer graphics, and medical imaging (e.g., HARDI). The purpose of this submission is to provide Matlab users with a set of functions for generating uniform sampling patterns and decompositions of a unit sphere (see demo pic).
DESCRIPTION OF MAIN FUNCTIONS
`ParticleSampleSphere.m`: generates an approximately uniform triangular tessellation of a unit sphere by using gradient descent to minimize a generalized electrostatic potential energy of a system of N charged particles. In this implementation, initial configuration of particles is based on random sampling of a sphere, but userdefined initializations are also permitted. This function can also be used to generate uniformly distributed sets of 2N particles comprised of N antipodal particle pairs. Since the optimization algorithm implemented in this function has O(N^2) complexity, it is not recommended that the function be used to optimize configurations of more than 1E3 particles (or particle pairs). Resolution of meshes obtained with this function can be increased to an arbitrary level with `SubdivideSphericalMesh.m`.
`SubdivideSphericalMesh.m`: increases resolution of triangular or quadrilateral spherical meshes. Given a base mesh, its resolution is increased by a sequence of k subdivisions. Suppose that No is the original number of mesh vertices, then the total number of vertices after k subdivisions will be Nk=4^k*(No – 2)+2. This relationship holds for both triangular and quadrilateral meshes.
`IcosahedronMesh.m`: generates triangular surface mesh of an icosahedron. Highquality spherical meshes can be easily obtained by subdividing this base mesh with the `SubdivideSphericalMesh.m` function.
`QuadCubeMesh.m`: generates quadrilateral mesh of a zerocentered unit cube. Highquality spherical meshes can be obtained by subdividing this base mesh with the `SubdivideSphericalMesh.m` function.
`SpiralSampleSphere.m`: generates N uniformly distributed point samples on a unit sphere using spiralbased sampling [2].
`RandSampleSphere.m`: performs uniform random or stratified random sampling of a unit sphere with N points.
DEMOS
Demos on how to generate various sampling patterns and decompositions of a unit sphere can be found here:
https://github.com/AntonSemechko/S2SamplingToolbox
REFERENCES
[1] Wikipedia. Thomson problem: https://en.wikipedia.org/wiki/Thomson_problem
[2] Christopher Carlson. How I Made Wine Glasses from Sunflowers: http://blog.wolfram.com/2011/07/28/howimadewineglassesfromsunflowers/
Anton Semechko (2020). Suite of functions to perform uniform sampling of a sphere (https://github.com/AntonSemechko/S2SamplingToolbox), GitHub. Retrieved .
1.6.0.0   Added option to 'ParticleSampleSphere.m' to allow construction of antipodally symmetric point configurations 

1.5.0.1   added GitHub repository link to the description 

1.5.0.0   migrated to GitHub 

1.4.0.0  Updated functions to remove dependency on 'TriRep' 

1.3.0.0   Added a function that performs sampling of a unit sphere along a spiral


1.2.0.0  Updated 'RandSamplSphere', the function used to obtain stratified sampling of the unit sphere prior to optimization. 

1.1.0.0  fixed minor bugs and updated description 
Inspired: Generate NonParallel Axes
Create scripts with code, output, and formatted text in a single executable document.
I'm also interested in seeing an implementation for hyperspheres. So far, the only thing I've seen for hyperspheres is Michael Völker's implementation, which takes a fairly different approach in that the output is a list of 0's and 1's, where "1" is considered a point on the hypersurface, in other words a pixel/voxel/etc. based approached, which is useful but becomes unwieldy in higher dimensions with high accuracy requirements (e.g. for a 7D hypersphere, easily turns into billions of points but still with a "blocky" appearance). You mention that no more than 1E3 points should be used in ParticleSampleSphere.m. I imagine this limitation would extend to or be exacerbated by higher dimensionality. Is that correct? Also, what lines in ParticleSampleSphere.m would you suggest starting with to modify?
@ Radion, functions included in this submission only deal with the S^2 case. However, the algorithm implemented by the ParticleSampleSphere function can be easily generalized to sample unit spheres embedded in dimension > 3.
So, is there an implementation for hyperspheres?
That would be great.
@ Paul Smith, triangulating vertices of a truncated icosahedron results in a highly irregular triangular surface mesh. So there was no reason for me to include this object here.
This shape (https://en.wikipedia.org/wiki/Truncated_icosahedron) is common in architecture (simplest combination of hexagons and pentagons, like a football). Is this possible with your code?
Hello, Yunchao. There are a couple of ways this code can be modified to suit your needs. What are the ratios of ellipsoid halfaxes? Feel free to email me directly. My email is listed in the function headers.
Hi, Anton, amazing work!
I'd like to generate uniformly distributed mesh on spheroid or ellipsoid. Will it be accomplished in this code? If not, could you provide some suggestions on how to modify it?
Thanks! Very useful.
Anton, the speed difference is not marginal if you generate a triangulation with millions of points. This is sometimes required for calculating support functions of geometrical differences. Secondly, it is not always convenient to use a pregenerated triangulation, especially if you use parallel calculations like we do in Ellipsoidal Toolbox (https://github.com/SystemAnalysisDptCMCMSU/ellipsoids). Passing triangulations that consume lets say 300 megabytes between the processes is much slower than regenerating it on the fly. For such applications I very much prefer a simple and superfast function instead of featurerich but slower one. You can think of spheretri as a very specialized archiver for triangulation data that can unzip it superfast.
@Peter : 'spheretri' cannot be used to generate spherical triangulation with arbitrary number of vertices, and thus has limited utility compared to what is offered here. Marginal improvement in subdivision speed of 'spheretri' is meaningless, as meshes only have to be generated once, and functions provided here can do so at interactive speeds.
This solution is relatively slow, check out spheretri  pure Matlab and still great performance:
>> tic;vMat,fMat]=spheretri(40962);toc;
Elapsed time is 0.034024 seconds.
https://www.mathworks.com/matlabcentral/fileexchange/58453spheretri
@George: Icosahedron and its subdivisions will produces the type of meshes you are looking for.
I'd like to generate geodesic spheres with just hexagons and pentagons and deleting some vertices from your output should do this. Do you have an approach for that?
@ Robert. The answer to your question is: No, the angles between adjacent particles will not be constant. In general, it is not possible to sample a sphere so that the distances between all pairs of adjacent particle are the same.
Great submission! Just to make sure: I am using SpiralSampleSphere. I assume the vertices are equally spaced in the 3D space, i.e. the angles between adjacent vectors (defined by the vertices coordinates) should be constant. Is this correct?
There is another very useful sampling scheme not in your suite:
https://github.com/fieldtrip/fieldtrip/blob/master/private/msphere.m
I believe it's referred to as the Rusin Disco Ball.
Hi Anton,
Can I use your toolbox to generate samples on half of the sphere? Or even more, on one given portion of the sphere?
Thanks!
Thanks so much. You saved me a lot of time by sharing this wonderful tool.
Very useful!
Hi Zara, I am glad you found this submission useful. Unfortunately there is not white paper to accompany this submission. Therefore, it is entirely up to you how and if you choose to cite this.
Hi Anton,
Thanks a lot for your nice code. It helped me a lot. How can we cite your work? Just an acknowledgement?
Hi, I clearly misunderstood the matlab function 'convhulln' here is a nice explanation of it:
http://www.mathworks.com/matlabcentral/answers/9298prettysimplequestionregardingconvhullnnowthatihavek
Hi Anton, Thank you for your soon reply. You are right I did not know that result of triangulation grows exponentially with number of dimensions (n), as you nicely illustrated with the random 100 points, so just for n=20 ‘convhulln’ will return an array of doubles with roughly 8.8E+12 elements and to store that I’ll need 65 400 GB of RAM. I only have 125 GB ☺. What is interesting here is that I do not get matlab ‘out of memory’ error but instead ‘The data is degenerate in at least one dimension – ND set of points lying in (N+1)D space’ . The good news is that I don’t really need to triangulate the points (I don’t need to know indices of the points that comprise the facets of the convex hull) even though it could be used to make sense of the data by sampling some areas and comparing them for instance. Well this takes the function ‘convhulln’ out of the game.
But, how about your code? Why are you amazed it worked for 525D? It basically addresses an optimization problem in 3D that can be for sure addressed in higher dimensions too (perhaps using a more sophisticated scheme like Conjugate gradient, etc). Well the thing is that you have done it and shared it, which all of us appreciate, and I am trying to expand its applicability.
Coming back to the hyper sphere in 525D, I don’t know how such a surface looks like ☺ but I believe that 3 points in a hyper dimensional space will define a triangle (just as they do in 2D and in 3D) so approximating the surface of an hypersphere with a set of triangles makes sense to me, therefore I wonder what do you mean by: “… 525 dimensional convex hull is a terribly poor approximation to a hypersphere…”
Once more thank you for your comments and no, I am not kidding!
Manuel  You apparently have NO idea how complex a triangulation of the surface of a 1000 dimensional hypersphere will be. The result would be IMMENSE.
In fact, I'm amazed that it succeeded for dimension 525. For example, try computing the convex hull of a fixed number of points on a hypersphere, what happens?
xyz = randn(100,2);
xyz = bsxfun(@rdivide,xyz,sqrt(sum(xyz.^2,2)));
t = convhulln(xyz);size(t)
ans =
100 2
xyz = randn(100,3);
xyz = bsxfun(@rdivide,xyz,sqrt(sum(xyz.^2,2)));
t = convhulln(xyz);size(t)
ans =
196 3
xyz = randn(100,5);
xyz = bsxfun(@rdivide,xyz,sqrt(sum(xyz.^2,2)));
t = convhulln(xyz);size(t)
ans =
1936 5
xyz = randn(100,7);
xyz = bsxfun(@rdivide,xyz,sqrt(sum(xyz.^2,2)));
t = convhulln(xyz);size(t)
ans =
24266 7
xyz = randn(100,9);
xyz = bsxfun(@rdivide,xyz,sqrt(sum(xyz.^2,2)));
t = convhulln(xyz);size(t)
ans =
287356 9
This last one started to take a significant amount of time, as you might expect.
Anyway, that 525 dimensional convex hull is a terribly poor approximation to a hypersphere. You are kidding yourself.
Hi Anton. Thank you for sharing your code. I wonder whether you have ever expanded your code to deal with higher dimensions. I need to do sampling on the hypersphere surface (525 – 1029 dimensions) and I accomplish that by normalizing the points draw from the Gaussian distribution, which gives a nice uniform distribution of points on the surface of the hipersphere. Nevertheless I need to refine the sampling so to avoid points laying too close to its neighbors and here is where triangular tessellation sounds just perfect.
I have bypass the 3D checking in your code and it seems to work fine for 525D but it fails when calling ‘convhulln’ function so no triangulation is done therefore I can no make sense of the new point distribution (what you call V).
Do you think ‘convhulln’ is not capable of dealing with 525D or your code ( ParticleSampleSphere ) simply can not be extended to higher dimensions?
Hi Sun, centroids of the spherical triangles can be estimated by taking the average of the three vertices and then projecting the resulting point on the sphere. For example, if x1, x2 and x3 are the vertices of the spherical triangle, you can estimate the centroid as c=(x1+x2+x3)/norm(x1+x2+x3)
Thanks a lot for your code. I am wondering how to use these vertex coordinates to generate the center point of each triangles?
Thanks for sharing you code, very nice function useful to generate toy examples for mesh parametrization. I'm reusing it to generate a roughly uniform distribution of points on a hemisphere
https://code.google.com/p/gerardus/source/browse/trunk/matlab/ManifoldToolbox/trihemisphere.m
@ Alec, thank you for pointing that out. I just resubmitted an updated version of RandSampleSphere function that does indeed produce stratified sampling of the unit sphere. It should appear shortly.
Thanks for sharing your code.
One note. The "stratified" option in RandSampleSphere is not really stratified. In fact it is the same as the "uniform" option: a uniform sampling following "Spherical Sampling by Archimedes’ Theorem" [Shao & Badler 96]
Dear Mr. Anton,
Your code is very nice and help me very much. But i need to create point in ellipsoid and I am confusing how to do it. Could you tell how to do it.
Many thanks.