version 1.6.0.0 (22.3 KB) by
Anton Semechko

Toolbox for generating uniform sampling patterns and decompositions of a unit sphere

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 so-called 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 re-emerged 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 user-defined 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. High-quality spherical meshes can be easily obtained by subdividing this base mesh with the `SubdivideSphericalMesh.m` function.

`QuadCubeMesh.m`: generates quadrilateral mesh of a zero-centered unit cube. High-quality 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 spiral-based 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/S2-Sampling-Toolbox

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/how-i-made-wine-glasses-from-sunflowers/

Anton Semechko (2020). Suite of functions to perform uniform sampling of a sphere (https://github.com/AntonSemechko/S2-Sampling-Toolbox), GitHub. Retrieved .

Created with
R2011b

Compatible with any release

**Inspired:**
Generate Non-Parallel Axes

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!Create scripts with code, output, and formatted text in a single executable document.

Sterling BairdI'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?

Anton Semechko@ 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.

Rodion TelyatnikSo, is there an implementation for hyperspheres?

That would be great.

Anton Semechko@ 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.

Paul SmithThis 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?

Anton SemechkoHello, Yunchao. There are a couple of ways this code can be modified to suit your needs. What are the ratios of ellipsoid half-axes? Feel free to e-mail me directly. My e-mail is listed in the function headers.

Yunchao YangHi, 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?

Sina SafaeiNathan VeilleuxMattia PugliattiChristof FallerThanks! Very useful.

Paul KassebaumW TomTimothy DavisPeter GagarinovAnton, 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 pre-generated triangulation, especially if you use parallel calculations like we do in Ellipsoidal Toolbox (https://github.com/SystemAnalysisDpt-CMC-MSU/ellipsoids). Passing triangulations that consume lets say 300 megabytes between the processes is much slower than re-generating it on the fly. For such applications I very much prefer a simple and super-fast function instead of feature-rich but slower one. You can think of spheretri as a very specialized archiver for triangulation data that can unzip it super-fast.

Anton Semechko@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.

Peter GagarinovThis 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/58453-spheretri

Anton Semechko@George: Icosahedron and its subdivisions will produces the type of meshes you are looking for.

George DishmanI'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?

Anton Semechko@ 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.

Robert JonesGreat 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?

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

BenHi 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!

Chibuzo NnonyeluThanks so much. You saved me a lot of time by sharing this wonderful tool.

Jorge Rubio ÁlvarezVery useful!

Anton SemechkoHi 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.

ZaraHi Anton,

Thanks a lot for your nice code. It helped me a lot. How can we cite your work? Just an acknowledgement?

ManuelHi, I clearly misunderstood the matlab function 'convhulln' here is a nice explanation of it:

http://www.mathworks.com/matlabcentral/answers/9298-pretty-simple-question-regarding-convhulln-now-that-i-have-k

ManuelHi 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!

John D'ErricoManuel - You apparently have NO idea how complex a triangulation of the surface of a 1000 dimensional hyper-sphere 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.

ManuelHi 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?

Anton SemechkoHi 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)

SunThanks a lot for your code. I am wondering how to use these vertex coordinates to generate the center point of each triangles?

Ramon CaseroThanks 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

Anton Semechko@ 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.

AlecThanks 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]

Cuong luuDear 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.