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