Computing the area of the cells of a Voronoi diagram for points sampled from the sphere

Question

The bigger picture is that I'm using QHULL to compute the convex hull of points on the sphere (which is a Delaunay Tesselation of the sphere's surface), and projecting Voronoi cells computed from the convex hull up to the surface of the sphere and then computing the area of associated with each node using spherical triangles http://mathworld.wolfram.com/SphericalTriangle.html (splitting each convex hull triangle into 6 fractional triangles and computing the area of the fractional spherical triangles separately, I'm getting the interior angle as the arc cosine of the dot product of unit vectors tangent to the surface of the sphere pointing from nodes to Voronoi vertexes and triangle edge midpoints), but I can deal with that complexity.

The problem I'm encountering is that for some very bad sample designs, the Voronoi vertex is outside the delaunay triangle yes it's possible see http://www.mathopenref.com/trianglecircumcenter.html but I wasn't expecting it, and I don't know that this answer Compute the size of Voronoi regions from Delaunay triangulation? applies to the case when Voronoi vertices are outside the triangle. I'm getting points with negative areas (using the spherical triangle approach based on connecting Voronoi vertexes to triangle edge midpoints and the triangle node to define fractional triangles)

Is the area closer to node than a Voronoi vertex that is outside the triangle (and closer to that node than the other 2 nodes of the same triangle) closer to that node than any other node in the delaunay tesselation? If so should I forget about connecting Voronoi vertices to the triangle edge midpoints and directly connect the Voronoi vertices of adjacent triangles and use that to compute the (spherical surface) area of the Voronoi region? This would still be simple to do because....

in MATLAB (where I'm prototyping the code, before I port it to C++) I get K (from convhulln) as a NK by 3 matrix of node indexes and can connect the Voronoi vertexes like this

    KS=sort(K,2);
    iK=(1:NK)';
    edges=sortrows([KS(:,1:2) iK; KS(:,[1 3]) iK; KS(:,2:3) iK]);
    edges=reshape(edges(:,3),2,1.5*NK);  %a consequence of this is that there will always be an even number of triangles on the convex hull of a sphere, otherwise Voronoi cells are not defined, which they must be because the convex hull is defined

"edges" contains the indexes of convex hull triangles which defines the Voronoi vertices? the 2 nodes that would have areas associated with that edge of a Voronoi cell are the ones in the first 2 columns of edges before they are discarded so like I said, that would be a simple modification, but I'm not sure it's mathematically correct. Do you know if it is?

Answer 1

You seem to be making your job more complicated than you need.

Note that if you take the normals to each of the facets of the hull (multiplied by the radius of your sphere), you will get the Voronoi vertices [1].

To compute the area of one of your Voronoi cells, you can iterate around the facets incident to one of your points, sequentially computing the Voronoi vertices. At the end you have a convex polygon forming the Voronoi cell about your point. You may compute its area by summing over each of the triangles defined by your input point along with each pair of adjacent Voronoi vertices.

As to your comment regarding the circumcenter being on the outside of some of your triangles, this is expected. You don't need to deform a triangle much from an equilateral triangle to achieve this condition.

[1] http://www.qhull.org/html/qdelaun.htm

Answer 2

This seemed to work

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