RandSphere: Generating Random Points On The Surface Of An nD Sphere
A few days ago I was in need of uniformly random points generation on the surface of sphere. I suddenly discovered that common parametrization with 2 angles is useless, it gathers too many points in the poles. At MathWorld you can find the complete demonstration, here the main lines:
To pick a random point on the surface of a unit sphere, it is incorrect to select spherical coordinates 
and 
from uniform distributions 
and ![phi in [0,pi]](http://mathworld.wolfram.com/images/equations/SpherePointPicking/Inline4.gif){kind=small}
, since the area element 
is a function of 
, and hence points picked in this way will be “bunched” near the poles (left figure above).
To obtain points such that any small area on the sphere is expected to contain the same number of points (right figure above), choose 
and 
to be random variates on 
. Then
 |
 |
 |
(1)
|
 |
 |
 |
(2)
|
gives the spherical coordinates for a set of points which are uniformly distributed over 
. This works since the differential element of solid angle is given by
 |
(3)
|
The distribution 
of polar angles can be found from
 |
(4)
|
by taking the derivative of (2) with respect to 
to get 
, solving (2) for 
, and plugging the results back in to (4) with 
to obtain the distribution
As we can see the distribution tends to infinite at the poles. With more research I discovered there are several There are several algorithms to generate random point on a sphere. A good tutorial can be found here.
At the end I have chosen the trig method, below you can find the Matlab function.
Update: 18/06/2010, According to Ed Hoyle comments the code has become n-dimensional
function X=RandSphere(N,dim) % RANDSPHERE % % RandSphere generates uniform random points on the surface of a unit radius % N-dim sphere centered in the origin . This script uses different algorithms % according to the dimension of point: % % -2D: random generation of theta [0 2*pi] % -3D: the "trig method". % -nD: Gaussian distribution % % SYNOPSYS: % % INPUT: % % N: integer number representing the number of points to be generated % dim: dimension of points, if omitted 3D is assumed as default % % OUTPUT: % % X: Nxdim double matrix representing the coordinates of random points % generated % % EXAMPLE: % % N=1000; % X=RandSphere(N); % hold on % title('RandSphere') % plot3(X(:,1),X(:,2),X(:,3),'.k'); % axis equal % % See also radnd, randi, randn, RandStream, RandStream/rand, % sprand, sprandn, randperm,RandStream/getDefaultStream. % %Authors: Luigi Giaccari,Ed Hoyle %errors check if nargin<1 %nargin >2 is handled by Matlab itself error('At least one input required') end if nargout~=1 error('One and ony one output required') end if N<=0 error('Negative number of points') end if round(N)~=N warning('N should be an integer. It will be rounded') N=round(N); end if nargin==1 dim=3;%default value for dimension of points end switch dim case 3 %3D %trig method X=zeros(N,dim);%preallocate X(:,3)=rand(N,1)*2-1;%z t=rand(N,1)*2*pi; r=sqrt(1-X(:,3).^2); X(:,1)=r.*cos(t);%x X(:,2)=r.*sin(t);%y case 2 %2D %just a random generation of theta X(:,2)=rand(N,1)*2*pi;%theta: use y as temp value X(:,1)=cos(X(:,2));%x X(:,2)=sin(X(:,2));%y otherwise %nD %use gaussian distribution X=randn(N,dim); X=bsxfun(@rdivide,X,sqrt(sum(X.^2,2))); end
Popularity: 2% [?]
Spot – A Linear-Operator Toolbox
Linear operators are at the core of many of the most basic algorithms for signal and image processing. Matlab’s high-level, matrix-based language allows us to naturally express many of the underlying matrix operations (e.g., computation of matrix-vector products and manipulation of matrices) and is thus a powerful platform on which to develop concrete implementations of these algorithms. Many of the most useful operators, however, do not lend themselves to the explicit matrix representations that Matlab provides. The aim of the Spot Toolbox is to bring the expressiveness of Matlab’s built-in matrix notation to problems for which explicit matrices are not practical.
Here’s simple example that shows how Spot operator primitives can be stitched together using familiar Matlab commands to form more complex operators:
n = 1000; x = (1:n)'; % the first column defines a circulant matrix F = opDFT(n); % create a DFT operator s = sqrt(n)*F*x; % eigenvalues of C C = real( F'*opDiag(s)*F ); % the circulant operator w = C*x; % apply C to a vector y = C'*w; % apply the adjoint of C to a vector z = C(end:-1:1,:)*y; % reverse the rows of C and apply to a vector double( C(1:5,1) )' % extract a few elements of the first column ans = 1 2 3 4 5
Download and more infos here
Popularity: 1% [?]
The fastest way to sort N numbers: sortN library
Contents
Have you ever thought about the fastest way to sort N numbers?
Last week I wrote a post about the fastest way to sort 3 numbers, this week I show you a tricky way to optimize the sort of N numbers.
The problem of sorting an array of numbers is probably the most studied in computer science. The efficiency of a sort algorithm depends essentially on the number of comparisons and swaps it needs to reach the solution. This post contains an optimized sort algorithm for small buckets (to my knowledge the fastest).
Currently most of sorting algorithms uses quick sort for large buckets and insertion sort for smaller ones. Both these algorithms uses an optimized way to swaps and compare numbers. The sortN library is a further optimization of this operations. It does first all the needed comparisons and then, only in the end, it swaps values to reach the sorted solution only in the end. It is in practice an enhanced Insertion sort algorithm.
The sortN library
The library is written in C++, but it can be written in any code, and supports arrays up to 7 elements. It can be applied , in theory, to array of any size but we will see this can not be realized in practice, only small arrays can be sorted.
The sortN libary can be considered enhanced insertion sort algorithm.
It contains a list of optimized function to sort 3,4,5,6,7 values. Each function is optimized to sort an N number of values.
Here is the approach to the solution:
N numbers can be disposed in Nc combinations where Nc is 
We may think to write a recursive code that reach the sorted solution in the way algorithms for permutations does.
The permutations represents all the possible solutions, depending on the value of each comparison (between 2 values of the array) a permutation is chosen or refused.
This allows, depending on the results of the comparisons, to build a decisionary tree that unequivocally leads to the permutation that represents the solution. Once we know the solution we can optimize swaps to reach it.
Notice that no swaps must be done done before we know the exact solution. This can be considered an enhanced insertion sort algorithm.
The number of comparisons possible are great even for a relative small array, we can not think about writing by hand the code, we need some function to do job for us. That’s what the sort library does, it prints the enhanced insertion sort algorithm.
There is a small disadvantage: the problem has size
The sortN function
You probably got confused, so here it is the code that I hope will clarify the concept.
The sort N function is a MATLAB script that prints the enhanced insertionsort algorithm for N values.
It recursively finds all the solutions and it prints optimized swaps.
function PrintSortNFunction(fid,N) %Prints a c++ sort N function. %fid is the file identifier where the function will be written. %N is the number of points the function will sort, for example: % %PrintSortNFunction(fid,4);%print the sort4 function %function definition fprintf(fid,'void sort%1.0f(double* a){\n',N); fprintf(fid,'double temp;\n'); x=1:N; PrintIS(x,1,2,fid); fprintf(fid,'}\n'); %The normal insertion sort: % void insertionSort(double* x,int length) % { % % % % double key; % int i; % % for(int j=1;j<length;j++) % { % key=x[j]; % i=j-1; % while(x[i]>key && i>=0) % { % x[i+1]=x[i]; % i--; % } % x[i+1]=key; % } % } end %print the isertion sort operations, takes in input the i an j paramater % (see code above). The x vector goes from 1 to N as matlab notation, the % prints instead goes from 0:N-1 as C++ notation function PrintIS(x,i,j,fid) N=length(x); %fprintf('%1.0f %1.0f\n',i,j); % we reached the end of the array if j>N%print solution fprintf(fid,'//'); for i=1:length(x) fprintf(fid,'%1.0f',x(i)-1); end fprintf(fid,'\n'); PrintSwaps(x,fid);%print swaps return end %if (x[i]>key) key=x(j); x0=x;%make a copy if i>=1 %if we are in the middle of the array ii=x(i);jj=x(j); fprintf(fid,'if (a[%1.0f]<a[%1.0f]){\n',ii-1,jj-1); %modify the solution x(i+2:j)=x(i+1:j-1); x(i+1)=key; %go forward PrintIS(x,j,j+1,fid); fprintf(fid,'}\n'); fprintf(fid,'else{\n'); PrintIS(x0,i-1,j,fid); fprintf(fid,'}\n') ; else %we have reached the beginning of the array %modify the solution x(2:j)=x(1:j-1); x(1)=key; %go forward PrintIS(x,j,j+1,fid) ; end end %prints swaps combination to go from the base array 1:N to the % sorted solution (permutation) function PrintSwaps(Perm,fid) N=length(Perm); %We first have to store a temp value %Perm=N:-1:1];test for debugging Base=1:N;%to record the change due to swaps swaps=zeros(N*3,2); c=1; for n=1:N i=Base(n); j=Perm(n); if j==i %do nothing already in place continue; end %we have to make Base(n) equal to Perm(n); i=find(Base==j); temp=Base(i);%phisically swap (i and n) Base(i)=Base(n); Base(n)=temp; %print swap %-1 stays for temp value swaps(c,1)=-1;swaps(c,2)=n; swaps(c+1,1)=n;swaps(c+1,2)=i; swaps(c+2,1)=i;swaps(c+2,2)=-1; c=c+3; end %detecting duplicate values in swaps (they do not need to be copied in %temp) swaps(c:end,:)=[];%remove unused values ind=false(c-1,1); for i=2:c-2%loop trough candidate swaps if swaps(i,1)==swaps(i+1,2) && swaps(i,2)==swaps(i+1,1) ind(i)=true;ind(i+1)=true; end end %remove duplicate values from swaps swaps(ind,:)=[]; for k=1:size(swaps,1) i=swaps(k,1);j=swaps(k,2); if i==-1 fprintf(fid,'temp=a[%1.0f];\n',j-1);%store the temp value elseif j==-1 fprintf(fid,'a[%1.0f]=temp;\n',i-1);%copy the temp value else fprintf(fid,'a[%1.0f]=a[%1.0f];\n',i-1,j-1);%keep copying end end end
Perfomances
I have compared the sortN functions with insertion sort and stl sort here are some results:
This program shows the perfomaces of the sortN library against
stl sort and insertions sort
Enter the number of buckets
10000000
BUCKET SIZE=3
STL sort of 3 points required 260.3610 ms
insertionSort of 3 points required 240.9796 ms
SortN sort of 3 points required 172.1386 ms
IS vs stl= 7.44406 %
SortN vs stl= 33.8846 %
SortN vs IS= 28.5671 %
BUCKET SIZE=4
STL sort of 4 points required 434.7305 ms
insertionSort of 4 points required 353.5941 ms
SortN sort of 4 points required 287.3990 ms
IS vs stl= 18.6636 %
SortN vs stl= 33.8903 %
SortN vs IS= 18.7206 %
BUCKET SIZE=5
STL sort of 5 points required 594.5557 ms
insertionSort of 5 points required 563.9854 ms
SortN sort of 5 points required 401.6870 ms
IS vs stl= 5.1417 %
SortN vs stl= 32.4391 %
SortN vs IS= 28.777 %
BUCKET SIZE=6
STL sort of 6 points required 883.9322 ms
insertionSort of 6 points required 735.4656 ms
SortN sort of 6 points required 517.2404 ms
IS vs stl= 16.7962 %
SortN vs stl= 41.4842 %
SortN vs IS= 29.6717 %
BUCKET SIZE=7
STL sort of 7 points required 1162.6356 ms
insertionSort of 7 points required 930.5345 ms
SortN sort of 7 points required 992.6320 ms
IS vs stl= 19.9634 %
SortN vs stl= 14.6223 %
SortN vs IS= -6.67331 %
You can run your own test by downloading the code below. As you can see for 7 numbers it becomes unadvantageous (dunno Why).
Downloads
The whole library can be build by downloading this:
You can download the library (up to sort7 here):
The quickest sort algorithm?
I don’t know, if you find something faster let me know.
What I known is that the sortN library can be further optimized. For example it can be made recursive (avoiding to write thousands lines of code), or it it can be sped up in the swaps: for each solution find the minimum number of swaps to allows to pass from the original vector to the sorted one (currently I can only avoid to copy values in “temp” twice).
Popularity: 1% [?]
Translator
Follow Us
Video Posts
