Code covered by the BSD License  

Highlights from
Suite of functions to perform uniform sampling of a sphere

5.0
5.0 | 5 ratings Rate this file 80 Downloads (last 30 days) File Size: 2.95 MB File ID: #37004
image thumbnail

Suite of functions to perform uniform sampling of a sphere

by

 

05 Jun 2012 (Updated )

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

| Watch this File

File Information
Description

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

Acknowledgements

This file inspired Generate Non Parallel Axes.

Required Products MATLAB
MATLAB release MATLAB 7.13 (R2011b)
MATLAB Search Path
/
/S2 Sampling Toolbox
Tags for This File   Please login to tag files.
Please login to add a comment or rating.
Comments and Ratings (13)
27 Jan 2015 Jorge Rubio Álvarez

Very useful!

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.

Comment only
09 Apr 2014 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?

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

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

Comment only
06 Apr 2014 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.

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

Comment only
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)

Comment only
31 Mar 2014 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?

11 Jul 2013 Ramón

Ramón (view profile)

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.

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

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

26 Mar 2015

- 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

Contact us