Generate blue noise samples on unit shpere with Mitchell's best candidate algorithm

BlueNoiseOnSphere.m

%==========================================================================
%   Generate blue noise samples on unit shpere with Mitchell's best
%   candidate algorithm.
%   Created by: Tiantian Xie (Dec. 12th, 2018)
%==========================================================================
  
NumberOfPointsOnSphere = 256;
CandidatesToGenerate = 10;     % increase this to increase the blue noise quality
BestPoints = zeros(NumberOfPointsOnSphere,3);

FirstPoint = rand(1,2);
BestPoints(1,:)=SamplingShpere(FirstPoint);

for i = 2:NumberOfPointsOnSphere
  
    % Generate CandidatesToGenerate candidates
    Candidates = SamplingShpere(rand(CandidatesToGenerate,2));
    
    BestCandidateID = 0;
    MaxDistance = 0;
    
    % Find the best candiate and add to the list
    for j = 1:CandidatesToGenerate
        Dist = SurfaceDistance(Candidates(j,:),BestPoints(1:(i-1),:));
        MinD(j) = min(Dist);
    end
    
    [~,BestCandidateID] = max(MinD);
    BestPoints(i,:)=Candidates(BestCandidateID,:);
end

SamplingSphere.m

%==========================================================================
%   Generate blue noise samples on unit shpere with Mitchell's best
%   candidate algorithm.
%   Created by: Tiantian Xie (Dec. 12th, 2018)
%==========================================================================
function [ p ] = SamplingShpere(x)
%SamplingShpere Generate samples on unit shpere given 2D random number in [0,1)
%   x: Nx2, N random points on 2D plane in [0,1)^2
%   p: Uniformly sampled points on unit shpere
                                            
    z = 1 - 2*x(:,1);
    sinTheta = 2*sqrt(max(zeros(size(x,1),1),x(:,1).*(1-x(:,1))));
    p = [cos(2*pi*x(:,2)).*sinTheta,sin(2*pi*x(:,2)).*sinTheta,z];                        
                                            
end

SurfaceDistance.m

%==========================================================================
%   Generate blue noise samples on unit shpere with Mitchell's best
%   candidate algorithm.
%   Created by: Tiantian Xie (Dec. 12th, 2018)
%==========================================================================
function [ d ] = SurfaceDistance(targetx,y)
%SurfaceDistance Calculate the distance of all y points to target point, targetx
%   Assume that all points are on unit sphere
%   targetx: One point on unit shpere
%         y: A vector of N points on units sphere
%         d: Distances to targetx from y

    vecNorm =@(x) sqrt(sum(x.^2,2));
    x = ones(size(y)).*targetx;
    %ref: https://en.wikipedia.org/wiki/Great-circle_distance
    d = atan2(vecNorm(cross(x,y)),dot(x,y,2));
end