File Exchange

image thumbnail

Suite of functions to perform uniform sampling of a sphere

version 1.3 (2.95 MB) by

Set of functions to perform approximately uniform sampling/decomposition of the spherical domain

4.85714
8 Ratings

53 Downloads

Updated

View License

There are a large number of applications that rely on uniform sampling/decomposition of a sphere embedded in 3D space. This submission provides a set of functions that can be used to obtain a variety of different sampling patterns and decompositions of the spherical domain (see demo pic).
Here is a brief summary of the main functions contained in this submission:
-------------------------------------------------------------------------------------------------------------------
'ParticleSampleSphere' : generates an approximately uniform triangular tessellation of a unit sphere by minimizing generalized electrostatic potential energy of a system of charged particles. By default, initializations are based on random sampling of a sphere, but user defined initializations are also permitted. Since the optimization algorithm implemented in this function has O(N^2) complexity, it is not recommended that 'ParticleSampleSphere' be used to optimize configurations of more than 1E3 particles. Resolution of the meshes obtained with this function can be increased to an arbitrary level with 'SubdivideSphericalMesh'.

-------------------------------------------------------------------------------------------------------------------
'SubdivideSphericalMesh': 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*(4^k–1). This relationship holds for both triangular and quadrilateral meshes.

-------------------------------------------------------------------------------------------------------------------
'IcosahedronMesh': generates a triangular mesh of an icosahedron. High-quality spherical meshes can be easily obtained by subdividing this base mesh with the 'SubdivideSphericalMesh' function.

-------------------------------------------------------------------------------------------------------------------
'QuadCubeMesh': generates a quadrilateral mesh of a cube. High-quality spherical meshes can be easily obtained by subdividing this base mesh with the 'SubdivideSphericalMesh' function.

-------------------------------------------------------------------------------------------------------------------
'SpiralSampleSphere': generates N uniformly distributed samples on a unit sphere using the technique described in http://blog.wolfram.com/2011/07/28/how-i-made-wine-glasses-from-sunflowers/

-------------------------------------------------------------------------------------------------------------------
'RandSampleSphere': produces uniform random or stratified random sampling of a unit sphere.

EXAMPLE 1:

% Uniformly distribute 200 particles across the surface of a unit sphere.
% This operation takes ~7 sec on my laptop (12GB RAM, 2.4 GHz processor)
[V,Tri,~,Ue]=ParticleSampleSphere('N',200);

% Visualize optimization progress
figure('color','w')
plot(log10(1:numel(Ue)),Ue,'.-')
set(get(gca,'Title'),'String','Optimization Progress','FontSize',40)
set(gca,'FontSize',20)
xlabel('log_{10}(Iteration #)','FontSize',30)
ylabel('Electrostatic Potential','FontSize',30)

% Visualize the mesh based on the optimal configuration of particles
figure('color','w')
subplot(1,2,1)
fv=struct('faces',Tri,'vertices',V);
h=patch(fv);
set(h,'EdgeColor','b','FaceColor','w')
axis equal
set(gca,'XLim',[-1.1 1.1],'YLim',[-1.1 1.1],'ZLim',[-1.1 1.1])
view(3)
grid on
set(get(gca,'Title'),'String','N=200 (base mesh)','FontSize',30)

% Subdivide base mesh twice to obtain a spherical mesh of higher complexity
fv_new=SubdivideSphericalMesh(fv,2);
subplot(1,2,2)
h=patch(fv_new);
set(h,'EdgeColor','b','FaceColor','w')
axis equal
set(gca,'XLim',[-1.1 1.1],'YLim',[-1.1 1.1],'ZLim',[-1.1 1.1])
view(3)
grid on
set(get(gca,'Title'),'String','N=3170 (after 2 subdivisions)','FontSize',30)

EXAMPLE 2:

% Generate a base icosahedron mesh
TR=IcosahedronMesh;
 
% Subvivide the base mesh and visualize the results
figure('color','w')
subplot(2,3,1)
h=trimesh(TR); set(h,'EdgeColor','b','FaceColor','w')
axis equal
for i=2:6
    subplot(2,3,i)
    TR=SubdivideSphericalMesh(TR,1);
    h=trimesh(TR); set(h,'EdgeColor','b','FaceColor','w')
    axis equal
    drawnow
end

EXAMPLE 3:
   
% Generate a quad cube mesh
fv=QuadCubeMesh;
 
% Subvivide the base mesh and visualize the results
figure('color','w')
subplot(2,3,1)
h=patch(fv); set(h,'EdgeColor','b','FaceColor','w')
view(3)
grid on
axis equal
set(gca,'XLim',[-1.1 1.1],'YLim',[-1.1 1.1],'ZLim',[-1.1 1.1])
for i=2:6
    subplot(2,3,i)
    fv=SubdivideSphericalMesh(fv,1);
    h=patch(fv); set(h,'EdgeColor','b','FaceColor','w')
    axis equal
    view(3)
    grid on
    drawnow
    set(gca,'XLim',[-1.1 1.1],'YLim',[-1.1 1.1],'ZLim',[-1.1 1.1])
end

Comments and Ratings (25)

W Tom

W Tom (view profile)

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.

Anton Semechko

Anton Semechko (view profile)

@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

Anton Semechko

Anton Semechko (view profile)

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

Anton Semechko

Anton Semechko (view profile)

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

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

fvff (view profile)

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

Ben (view profile)

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!

Anton Semechko

Anton Semechko (view profile)

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

Zara (view profile)

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

Manuel

Manuel (view profile)

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

Manuel (view profile)

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!

John D'Errico

John D'Errico (view profile)

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

Manuel (view profile)

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?

Anton Semechko

Anton Semechko (view profile)

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

Sun (view profile)

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

Ramon Casero

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

Anton Semechko

Anton Semechko (view profile)

@ 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

Alec (view profile)

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

- 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

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

1.1

fixed minor bugs and updated description

MATLAB Release
MATLAB 7.13 (R2011b)
Acknowledgements

Inspired: Generate Non-Parallel Axes

Download apps, toolboxes, and other File Exchange content using Add-On Explorer in MATLAB.

» Watch video