问计算从球体采样点的Voronoi图的单元面积

EN

function [x,y,z,ptArea,K,minRad,xv,yv,zv,xvc,yvc,zvc,VorTriNodes]=voronoiSphereArea2(x,y,z) 
%x, y, z are the inputs but normalized to be unit vectors (projected onto
%        the surface of the sphere centered at 0,0,0 with radius 1)
%ptArea  is a numel(x) by 1 vector (same size as x, y, and z) that contains
%        area of the surface of the unit sphere that is closer to the
%        corresponding x y z point than any other pt (this is the area of 
%        projection of the voronoi cell projected onto the surface of the 
%        sphere with radius of 1)
%K       is the output of convhulln([x y z]), it is the tesselation of 
%        points x y z (the list point indexes that make up the triangles 
%        on the convex hull of the sphere)
%minRad  is the minimum distance between the origin (center of the sphere)
%        and any triangle in the tesselation of (normalized) x y z; this is
%        less than 1 because the triangles are the equivalent of chords of
%        a circle (nodes of the triangle are on the surface of the sphere,
%        but the rest of the triangle is below the surface of the sphere,
%        the voronoi vertex of a triangle is the closest point in that 
%        triangle to the center of the sphere), minRad is highly useful
%        when doing barycentric interpolation in x y z space (putting 
%        points to be interpolated at this distance from the origin, [the 
%        largest distance guaranteed to be in the convex hull regardless of
%        angular position, actually you'll want to decrease this slightly
%        to protect against round off error] makes tsearchn [of volumetric
%        Delaunay tesselation formed augmenting each triangle in the 
%        convex hull tesselation with the origin] run faster, it's very 
%        slow if you use an interpolation radius of 0.5 but very fast if 
%        you put the interpolation points just under the surface of the 
%        sphere)
%xv, yv, zv are the x y z coordinates of the voronoi vertices of each
%        triangle in K projected up to the surface of the sphere (for use
%        in plotting "cells" associated with each point) these are listed
%        in the same order as K
%xvc, yvc, zvc are spherical triangular patches whose corners are 1 node 
%        (x, y, z) and 2 voronoi vertices (xv, yv, zv); this is a 
%        30 by 3*size(K,1) matrix, each triangle is a fraction of the
%        voronoi cell belonging to the node (x, y, z)
%VorTriNodes a 1 by 3*size(K,1) vector that identifies the node owning
%        the spherical triangle in the same columns of xvc, yvc, zvc; this
%        is used for plotting piecewise constant point values (where the
%        pieces are the voronoi cells belonging to their respective points)

Npts=numel(x);
invr=(x.^2+y.^2+z.^2).^-0.5;
x=x.*invr;
y=y.*invr;
z=z.*invr;

i4=[1:3 1]';
i3=[3 1 2]';

K=convhulln([x y z]); %right now this is K
NK=size(K,1);

%determine the voronoi cell edges (between the voronoi vertexes of 2 convex
%hull trianlges that share a convex hull triangle edge) and the 2 nodes (x,
%y, z points) that share that edge, these define triangles that are 
%components of the voronoi cells around the 2 nodes
KS=sort(K,2);
iK=(1:NK)';
VorTri=reshape(sortrows([KS(:,1:2) iK; KS(:,[1 3]) iK; KS(:,2:3) iK])',[3 2 1.5*NK]);
VorTriNodes=squeeze(VorTri(1:2,1,:)); %these 2 nodes (in a column) share 
VorTriVerts=squeeze(VorTri(3,:,:)); %the cell edge between these 2 vertices

K=K'; %right now this is K transpose
xn=x(K); %a 3 by NK matrix for the x coordinates of a convex hull triangle, n is for "nodes"
yn=y(K); %a 3 by NK matrix for the y coordinates of a convex hull triangle, n is for "nodes"
zn=z(K); %a 3 by NK matrix for the z coordinates of a convex hull triangle, n is for "nodes"

%calculate the 3 edge lengths (in x y and z coordinates) for each
%triangle in the convex hull tesselation
dx=diff(xn(i4,:));
dy=diff(yn(i4,:));
dz=diff(zn(i4,:));

%calculate the voronoi vertices
xv1=-(dz(1,:).*dy(3,:)-dy(1,:).*dz(3,:)); %negative needed to correct sign convention to keep vornoi vertex on same side of sphere as triangle when we normalize it
yv1=-(dx(1,:).*dz(3,:)-dz(1,:).*dx(3,:));
zv1=-(dy(1,:).*dx(3,:)-dx(1,:).*dy(3,:));
normconst=(xv1.^2+yv1.^2+zv1.^2).^-0.5; %voronoi vertex will be on surface of sphere (radius 1), negative needed to correct sign convention to keep vertex from being projected to opposite side of sphere
xv1=xv1.*normconst; %x component of unit vector pointing from origin through the triangle voronoi vertex (voronoi vertex projected up to the surface of the unit sphere)
yv1=yv1.*normconst; %y component of unit vector pointing from origin through the triangle voronoi vertex (voronoi vertex projected up to the surface of the unit sphere)
zv1=zv1.*normconst; %z component of unit vector pointing from origin through the triangle voronoi vertex (voronoi vertex projected up to the surface of the unit sphere)

xv2=-(dz(2,:).*dy(1,:)-dy(2,:).*dz(1,:));
yv2=-(dx(2,:).*dz(1,:)-dz(2,:).*dx(1,:));
zv2=-(dy(2,:).*dx(1,:)-dx(2,:).*dy(1,:));
normconst=(xv2.^2+yv2.^2+zv2.^2).^-0.5; %voronoi vertex will be on surface of sphere (radius 1), negative needed to correct sign convention to keep vertex from being projected to opposite side of sphere
xv2=xv2.*normconst; %x component of unit vector pointing from origin through the triangle voronoi vertex (voronoi vertex projected up to the surface of the unit sphere)
yv2=yv2.*normconst; %y component of unit vector pointing from origin through the triangle voronoi vertex (voronoi vertex projected up to the surface of the unit sphere)
zv2=zv2.*normconst; %z component of unit vector pointing from origin through the triangle voronoi vertex (voronoi vertex projected up to the surface of the unit sphere)

xv3=-(dz(3,:).*dy(2,:)-dy(3,:).*dz(2,:));
yv3=-(dx(3,:).*dz(2,:)-dz(3,:).*dx(2,:));
zv3=-(dy(3,:).*dx(2,:)-dx(3,:).*dy(2,:));
normconst=(xv3.^2+yv3.^2+zv3.^2).^-0.5; %voronoi vertex will be on surface of sphere (radius 1), negative needed to correct sign convention to keep vertex from being projected to opposite side of sphere
xv3=xv3.*normconst; %x component of unit vector pointing from origin through the triangle voronoi vertex (voronoi vertex projected up to the surface of the unit sphere)
yv3=yv3.*normconst; %y component of unit vector pointing from origin through the triangle voronoi vertex (voronoi vertex projected up to the surface of the unit sphere)
zv3=zv3.*normconst; %z component of unit vector pointing from origin through the triangle voronoi vertex (voronoi vertex projected up to the surface of the unit sphere)

%round off error is a big concern average all three possible calculations
xv=(xv1+xv2+xv3)/3;
yv=(yv1+yv2+yv3)/3;
zv=(zv1+zv2+zv3)/3;
normconst=(xv.^2+yv.^2+zv.^2).^-0.5; %use this to project voronoi vertex 
%up to surface of (radius 1) sphere 
xv=xv.*normconst;
yv=yv.*normconst;
zv=zv.*normconst;

if(any(imag(normconst(:)))||any(isnan(normconst(:))))
    error('imaginary numbers or NANs at A');
end

minRad=min(abs(xv.*xn(1,:)+yv.*yn(1,:)+zv.*zn(1,:))); %the minimum distance
%from the origin to any point on the tesselation of the sphere (which will 
%be less than 1 because convex hull facets cut below the surface of the 
%sphere like a chord on a circle) found by taking the dot product of a node
%and the unit vector that passes through the origin (center of sphere) and 
%voronoi vertex of the convex hull triangle determines the distance of the 
%convex hull triangle's voronoi vertex from the center of the sphere.

%define shape of fractional spherical triangle patches to be used in
%piecewise constant plotting of point values
ds=(0:0.1:0.9)';
j1=ones(10,1);
j2=2*j1;
xn=x(VorTriNodes);
xV=xv(VorTriVerts);
yn=y(VorTriNodes);
yV=yv(VorTriVerts);
zn=z(VorTriNodes);
zV=zv(VorTriVerts);
VorTriNodes=VorTriNodes(:); %column vector to do the dot product easily

xvc=reshape(...
    [xn(j1,:)+ds*(xV(1,:)-xn(1,:)); xV(j1,:)+ds*diff(xV); xV(j2,:)+ds*(xn(1,:)-xV(2,:));...
     xn(j2,:)+ds*(xV(1,:)-xn(2,:)); xV(j1,:)+ds*diff(xV); xV(j2,:)+ds*(xn(2,:)-xV(2,:))],...
    30,3*NK);

yvc=reshape(...
    [yn(j1,:)+ds*(yV(1,:)-yn(1,:)); yV(j1,:)+ds*diff(yV); yV(j2,:)+ds*(yn(1,:)-yV(2,:));...
     yn(j2,:)+ds*(yV(1,:)-yn(2,:)); yV(j1,:)+ds*diff(yV); yV(j2,:)+ds*(yn(2,:)-yV(2,:))],...
    30,3*NK);

zvc=reshape(...
    [zn(j1,:)+ds*(zV(1,:)-zn(1,:)); zV(j1,:)+ds*diff(zV); zV(j2,:)+ds*(zn(1,:)-zV(2,:));...
     zn(j2,:)+ds*(zV(1,:)-zn(2,:)); zV(j1,:)+ds*diff(zV); zV(j2,:)+ds*(zn(2,:)-zV(2,:))],...
    30,3*NK);

invrvc=(xvc.^2+yvc.^2+zvc.^2).^-0.5;
xvc=xvc.*invrvc;
yvc=yvc.*invrvc;
zvc=zvc.*invrvc;

if(any(imag(invrvc(:)))||any(isnan(invrvc(:))))
    error('imaginary numbers or NANs at B');
end


%time to compute the area of the spherical triangle components of the
%voronoi cells areound each point

%first define the corners of the voronoi cell triangles and their edge
%lengths (direction matters)
xt=reshape([xn(1,:); xV; xn(2,:); xV],3,3*NK); dx=diff(xt(i4,:));
yt=reshape([yn(1,:); yV; yn(2,:); yV],3,3*NK); dy=diff(yt(i4,:));
zt=reshape([zn(1,:); zV; zn(2,:); zV],3,3*NK); dz=diff(zt(i4,:));


%now going to compute unit vectors tangent to surface in direction of
%great circles connection nodes of voronoi cell component triangles

%dx, dy, dz are chord distances in forward (a) direction
yada=xt.*dx+yt.*dy+zt.*dz;
dirxa=dx-xt.*yada; %x component of derivative/tangent vector at point
dirya=dy-yt.*yada; %y component of derivative/tangent vector at point
dirza=dz-zt.*yada; %z component of derivative/tangent vector at point
yada=(dirxa.^2+dirya.^2+dirza.^2);
yada=(yada+(yada==0)).^-0.5; %protect against zero magnitude vectors which 
%causes NaNs without this protection, yes it happened in a Tensor Product 
%grid (vornoi vertex was outside triangle on top of other voronoi vertex)
dirxa=dirxa.*yada; %x component of unit vector direction 
dirya=dirya.*yada; %y component of unit vector direction
dirza=dirza.*yada; %z component of unit vector direction

if(any(imag(yada(:)))||any(isnan(yada(:))))
    error('imaginary numbers or NANs at C');
end

%make dx, dy, dz chord distances in backward (b) direction
dx=-dx(i3,:); dy=-dy(i3,:); dz=-dz(i3,:);
yada=xt.*dx+yt.*dy+zt.*dz;
dirxb=dx-xt.*yada; %x component of derivative/tangent vector at point
diryb=dy-yt.*yada; %y component of derivative/tangent vector at point
dirzb=dz-zt.*yada; %z component of derivative/tangent vector at point
yada=(dirxb.^2+diryb.^2+dirzb.^2);
yada=(yada+(yada==0)).^-0.5;%protect against zero magnitude vectors which 
%causes NaNs without this protection, yes it happened in a Tensor Product 
%grid (vornoi vertex was outside triangle on top of other voronoi vertex)
dirxb=dirxb.*yada; %x component of unit vector direction 
diryb=diryb.*yada; %y component of unit vector direction
dirzb=dirzb.*yada; %z component of unit vector direction

if(any(imag(yada(:)))||any(isnan(yada(:))))
    error('imaginary numbers or NANs at D');
end

%fractional spherical triangle area is the sum of the 3 interior angles
%minus pi times the square of the sphere's radius (which is 1).  Interior
%angle is the arc cosine of the dot product of the forward and backward
%unit vectors (in the directions away from one corner of spherical triangle
%to the other 2 corners)
%min and max are needed to fix round off error that could cause imaginary 
%angles via arc-hyperbolic-cosine, but max would replace NaN with -1 which 
%makes angle pi when should have been zero, see the 
%yada=(yada+(yada==0)).^-0.5; protection above
areacontrib=sum(acos(min(max(dirxa.*dirxb+dirya.*diryb+dirza.*dirzb,-1),1)))-pi;

%fix round off error that could cause "negative zero" areas
areacontrib=areacontrib.*~((areacontrib<0)&(areacontrib>-(2^-35)));

if(any(imag(areacontrib))||any(areacontrib<0))
    %area might be imaginary or negative if there's a bug (other than round
    %off error which should have been fixed).
    save DEBUGME_voronoiSphereArea
    error('imaginary or negative areas at E Npts=%d',Npts);
end

ptArea=zeros(Npts,1);
for ipt=1:Npts
    ptArea(ipt)=areacontrib*(VorTriNodes==ipt); %matrix ops dot product to
    %sum the areas of all spherical triangles that make of the point's 
    %voronoi cell
end

%SphereAreaDiv4pi=sum(ptArea)/(4*pi) %is 1 to 8+ sig figs
ptArea=ptArea*(4*pi/sum(ptArea)); %"fix" round off error in ptArea

K=K'; %this is no longer K transpose
xv=xv(:); %for drawing voronoi cell edges
yv=yv(:);
zv=zv(:);
VorTriNodes=VorTriNodes'; %for plotting piecewise constant voronoi cells as
%being made up of triangles with 1 node and 2 vertices.
end