File Exchange

image thumbnail

Suite of functions to perform uniform sampling of a sphere

version 1.6.0.0 (22.3 KB) by Anton Semechko
Toolbox for generating uniform sampling patterns and decompositions of a unit sphere

26 Downloads

Updated 13 Jun 2019

GitHub view license on GitHub

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/

Cite As

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

Comments and Ratings (34)

@ 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 Smith

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 half-axes? Feel free to e-mail me directly. My e-mail 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?

Sina Safaei

Thanks! Very useful.

W Tom

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

@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/58453-spheretri

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

fvff

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.

Ben

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.

Zara

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

Manuel

Hi, 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

Manuel

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

Manuel

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)

Sun

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.

Alec

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]

Cuong luu

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.

Updates

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
- Updated 'SubdivideSphericalMesh' to work with quad meshes
- Added 'QuadCubeMesh' function to enable spherical decomposition with a quad mesh

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

MATLAB Release Compatibility
Created with R2011b
Compatible with any release
Platform Compatibility
Windows macOS Linux
Acknowledgements

Inspired: Generate Non-Parallel Axes