Code covered by the BSD License  

Highlights from
Uniform Sampling of a Sphere

5.0

5.0 | 4 ratings Rate this file 53 Downloads (last 30 days) File Size: 1.27 MB File ID: #37004
image thumbnail

Uniform Sampling of a Sphere

by

 

05 Jun 2012 (Updated )

Create an approximately uniform triangular tessellation of a unit sphere

| Watch this File

File Information
Description

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 named 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 the 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, electrostatics, geophysics, computer graphics and medical imaging (HARDI).

DESCRIPTION OF THE FUNCTIONS:

The main function is titled 'ParticleSampleSphere' and allows you to create an approximately uniform triangular tessellation of the unit sphere by minimizing the generalized electrostatic potential energy (aka Reisz s-energy) of the system of charged particles. Effectively, this function produces a locally optimal solution to the problem of finding a minimum Reisz s-energy configuration of N equal charges (s=1 corresponds to the original Thomson problem). The solution is obtained by iterative modification of particle positions along the negative gradient of the energy using an adaptive Gauss-Seidel update scheme. By default, the initializations are based on stratified sampling of the sphere [3], but user defined initializations are also permitted. It must be emphasized that in this function, all energy calculations are based on the geodesic and not Euclidean distances.

Due to high computational complexity of the problem, it is not recommended that 'ParticleSampleSphere' be used to solve the systems composed of more than 1E3 particles. To obtain uniform tessellations of the sphere containing of more than 1E3 nodes, I included a function titled 'SubdivideSphericalMesh'. This routine uses triangular quadrisection to subdivide the input mesh an arbitrary number of times and automatically re-projects the newly inserted vertices onto the unit sphere after every iteration. You can use the following expression to estimate the number of mesh vertices after k subdivisions

 Nk ~ No*4^k

where No is the original number of vertices.

For convenience, I also included the function titled 'IcosahedronMesh' which generates a triangular mesh of an icosahedron. High-quality spherical meshes can be easily obtained by subdividing this base mesh with an aforementioned 'SubdivideSphericalMesh' function. Finally, the function titled 'RandSampleSphere' can be used to obtain random or stratified sampling of the unit sphere.

 
EXAMPLE:

% Uniformly distribute 162 particles across the surface of the unit sphere
[V,Tri,~,Ue]=ParticleSampleSphere('N',162); % this operation takes ~8 sec on my machine (6GB RAM & 2.8 GHz processor)

% Visualize optimization progress
figure, plot(0:numel(Ue)-1,Ue,'.-')
set(get(gca,'Title'),'String','Optimization Progress','FontSize',20)
xlabel('Iteration #','FontSize',15)
ylabel('Reisz s-energy','FontSize',15)

% Visualize the mesh based on the optimal configuration of particles
TR=TriRep(Tri,V);
figure, h=trimesh(TR); set(h,'EdgeColor','b'), axis equal

% Subdivide the base mesh twice to obtain a spherical mesh of higher complexity
TR_2=SubdivideSphericalMesh(TR,2);
figure, h=trimesh(TR_2); set(h,'EdgeColor','b'), axis equal

 REFERENCES:

[1] Thomson problem: http://en.wikipedia.org/wiki/Thomson_problem
[2] Saff, E.B., Kuijlaars, A.B.J. (1997) ‘Distributing many points on a sphere’
[3] Shao M.-Z., Badler, N. (1996),'Spherical sampling by Archimedes' theorem'

---------------------------------------

If you found any of this stuff useful, please leave a comment. Cheers!

Acknowledgements

This file inspired Generate Non Parallel Axes.

Required Products MATLAB
MATLAB release MATLAB 7.13 (R2011b)
Tags for This File   Please login to tag files.
Please login to add a comment or rating.
Comments and Ratings (12)
10 Apr 2014 Anton Semechko

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.

09 Apr 2014 Zara

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

07 Apr 2014 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

07 Apr 2014 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!

06 Apr 2014 John D'Errico

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.

06 Apr 2014 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?

31 Mar 2014 Anton Semechko

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)

31 Mar 2014 Sun

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

11 Jul 2013 Ramón

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

23 May 2013 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.

20 May 2013 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]

16 Sep 2012 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
05 Jun 2012

fixed minor bugs and updated description

23 May 2013

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

Contact us