function vx = CVPsphere(N,varargin)
% Calc centroidal voronoi partition of sphere of N points
%
% vx = CVPsphere(N)
% vx = CVPsphere(N,segments)
% vx = CVPsphere(N,segments,tol)
%
% Input: N number of points to distribute. Segments: resolution of sphere
% bitmap, default 100. Tol: tolerance. Default 1e-4.
% Returns: Nx3 array of points on the sphere.

if (nargin>1) segments=varargin{1}; else segments=100; end
if (nargin==3) tol = varargin{2}; else tol=1e-4; end

[X,Y,Z]=sphere(segments);

vx=randn(N,3);
for i=1:N
    vx(i,:)=vx(i,:)/norm(vx(i,:));
end

delta=100;
while (delta>tol)
    d=X*0+1e11;
    c=X*0;
    for i=1:N
        dx=X-vx(i,1);
        dy=Y-vx(i,2);
        dz=Z-vx(i,3);
        d2=dx.^2+dy.^2+dz.^2;
        d=min(d,d2);
        c=c.*(d~=d2)+i*(d==d2);
    end

    delta=0;
    weight = (1+2*pi*sqrt(1-Z.*Z)*size(X,1));
    for i=1:N;
        mask=(c==i);

        XP=X.*mask.*weight;
        YP=Y.*mask.*weight;
        ZP=Z.*mask.*weight;
        vold=vx(i,:);
        vx(i,:)=[mean(XP(:)) mean(YP(:)) mean(ZP(:))];
        vx(i,:)=vx(i,:)/norm(vx(i,:));
        delta=delta+norm(vx(i,:)-vold)/N;
    end
end