Voronoi diagram on a sphere

[Abel Brown]

iv been able to make voronoi diagrams easily on 2D plane but can not seem to get it right for a sphere.

iv googled around and found some spherical V. diagrams but cant seem to get matlab to do this

http://people.sc.fsu.edu/~%20burkardt/f_src/sxyz_voronoi/gen_00100.png

I have some points (lat, lon, ht or x, y, z) and would like to make sphereical Voronoi diagram and plot it. Any help would be awesome!

cheers!

[Bruno Luong]

"Abel Brown" <brown...@osu.edu> wrote in message <hd0e8h$c9$1...@fred.mathworks.com>...

Voronoi cell is the "dual" of delaunay.

Use convhulln on your (x,y,z) that gives something close to Spherical Delaunay (assuming the points are close enough to neglect the spherical curvature). Next draw the median line (on the sphere) separating each vertice its delaunay neighbor, the will form a polygonal of the voronoi cell. Do it repetitively for all vertices.

Bruno

[Damian Sheehy]

Hi Abel,

If you have release in R2009a or later, the computational geometry tools
make this task relatively easy.
The Voronoi diagram is the dual of the Delaunay triangulation; the
triangulation lying on the on the surface of the sphere.

% Create the points on the surface of the sphere
theth_a = 2*pi*rand(100,1);
phi = pi*rand(100,1);
x = cos(theth_a).*sin(phi);
y = sin(theth_a).*sin(phi);
z = cos(phi);

% Create a Delaunay triangulation of the point set.
% This is a solid triangulation composed of tetrahedra.
dt = DelaunayTri(x,y,z);

% Extract the surface (boundary) triangulation, in this instance it is
actually the convex hull
% so you could use the convexHull method also.
ch = dt.freeBoundary();

% Create a Triangulation Representation from this surface triangulation
% We will use it to compute the location of the voronoi vertices
(circumcenter of the triangles),
% and the dual edges from the triangle neighbor information.

tr = TriRep(ch,x,y,z);
numt = size(tr,1);
T = (1:numt)';
neigh = tr.neighbors();
cc = tr.circumcenters();
xcc = cc(:,1);
ycc = cc(:,2);
zcc = cc(:,3);
idx1 = T < neigh(:,1);
idx2 = T < neigh(:,2);
idx3 = T < neigh(:,3);
neigh = [T(idx1) neigh(idx1,1); T(idx2) neigh(idx2,2); T(idx3)
neigh(idx3,3)]';
plot3(xcc(neigh), ycc(neigh), zcc(neigh),'-r');
axis equal;

This is just a 3D extension of the following example;
http://www.mathworks.com/products/demos/shipping/matlab/demoDelaunayTri.html#24

Best regards,

Damian

"Abel Brown" <brown...@osu.edu> wrote in message

news:hd0e8h$c9$1...@fred.mathworks.com...

[Abel Brown]

"Bruno Luong" <b.l...@fogale.findmycountry> wrote in message <hd0hkr$pq1$1...@fred.mathworks.com>...

Thanks for your reply Bruno

These points that I have are distributed all over the globe. Some of the points are separated by thousands of kilometers so will not be able to neglect the spherical curvature.

I was thinking to do some distance preserving projection on to a plane, do the V. diagram in 2D, and then map it back. But in this case the polygons will not be "continuous". The polygons at the edges will go off to infinity where as on a sphere all polygons will have a finite area.

Is there a way to use voronoin to do this? this will do the 3D Voronoi diagram but need to plot the polygons on a sphere ...

[Abel Brown]

"Damian Sheehy" <Damian...@mathworks.com> wrote in message <hd1a6r$rq0$1...@fred.mathworks.com>...

Damian

Thank you very much! This is exactly what I needed!

Ultimately, I need to calculate the area of each of these polygon (which is no problem).

Typically, using [V,C] = voronoin(X), the index of the i'th polygon

% faces
polygonIndx = C{i}

Then to get the vertex information just index into V

vX = V(polygonIndx,1);
vY = V(polygonIndx,2);

Is there a way, using your example, to get the vertices for each polygon?

Should I just computed the convexhull of each of the original points using the Voronoi points computed?

Thanks!!

[Damian Sheehy]

"Abel Brown" <brown...@osu.edu> wrote in message

news:hd1e8s$gf6$1...@fred.mathworks.com...

Hi Abel,

You should realize this is an approximate Voronoi diagram.
The fewer points you have the less accurate it will be.
Also, the Voronoi vertices are not exactly on the surface of the sphere, as
the triangle facet is not "draped" over the surface, but it's not difficult
to project them.

If you wish to compute the Voronoi region associated with point "i " you
can do so as follows;

1) Compute the set of triangles attached to point i. The triangles are
arranged in a CCW cycle around the point i.
Ti = tr.vertexAttachments(i)

2) The positions of the vertices defining the i'th Voronoi region are the
circumcenters of these triangles
ccTi = tr.circumcenters(Ti)

3) The Voronoi region may be non-planar.
To compute the area break it into triangles defined by the point i and
each edge of the Voronoi region.
The location of point i is x(i), y(i), z(i)
the two edge points are ccTi(1,:) to ccTi(2,:)

I suggest you first try computing a simple 2D Voronoi diagram using this
approach, and when you understand the process apply it to your problem.
Take your time, understand the steps, and if you get stuck feel free to
contact me directly.

Damian

[Abel Brown]

> Hi Abel,
>
> You should realize this is an approximate Voronoi diagram.
> The fewer points you have the less accurate it will be.
> Also, the Voronoi vertices are not exactly on the surface of the sphere, as
> the triangle facet is not "draped" over the surface, but it's not difficult
> to project them.
>
> If you wish to compute the Voronoi region associated with point "i " you
> can do so as follows;
>
> 1) Compute the set of triangles attached to point i. The triangles are
> arranged in a CCW cycle around the point i.
> Ti = tr.vertexAttachments(i)
>
> 2) The positions of the vertices defining the i'th Voronoi region are the
> circumcenters of these triangles
> ccTi = tr.circumcenters(Ti)
>
> 3) The Voronoi region may be non-planar.
> To compute the area break it into triangles defined by the point i and
> each edge of the Voronoi region.
> The location of point i is x(i), y(i), z(i)
> the two edge points are ccTi(1,:) to ccTi(2,:)
>
> I suggest you first try computing a simple 2D Voronoi diagram using this
> approach, and when you understand the process apply it to your problem.
> Take your time, understand the steps, and if you get stuck feel free to
> contact me directly.
>
>
> Damian
>

Damian

This works great!

No problems computing polygon for each point using your instructions!

Thanks again
-abel

[Felipe Pedreros]

"Damian Sheehy" <Damian...@mathworks.com> wrote in message <hd1a6r$rq0$1...@fred.mathworks.com>...

Hi Damian,

Thanks a lot for your implementation. I've been using it to compute voronoi diagrams of a network of stations spread in the earth's surface, but I got a problem when the distribution of points (the stations) is only or mostly (90%) in one hemisphere. For example, I used the dataset below (given in degree), and the problem is when plotting the diagram, some of the edges of the diagram are overlap (considering the projection into the sphere's surface). By looking the plot3, some voronoi vertices are actually joint "through the earth". This make sense because there are points in only one hemisphere and the surface triangulation "cut" the sphere right in region of the equator. So, I was thinking if there is any constraint (or trick) to add to the implementation in order to avoid surface triangulation that cut the sphere.

phi = [51.77 43.79 36.18 54.72 60.22 43.46 57.38 36.10 43.81 40.50 36.97 27.94 38.52 43.79];

theta = [257.77 234.56 219.73 339.49 335.62 272.83 348.08 219.92 228.03 3.11 25.17 15.39 28.72 318.43];

I'd appreciate your help in this "special case".

Many thanks
Felipe

[Bruno Luong]

"Felipe Pedreros" <f_ped...@hotmail.com> wrote in message <kiq4o6$pmo$1...@newscl01ah.mathworks.com>...

> Hi Damian,
>
> Thanks a lot for your implementation. I've been using it to compute voronoi diagrams of a network of stations spread in the earth's surface, but I got a problem when the distribution of points (the stations) is only or mostly (90%) in one hemisphere. For example, I used the dataset below (given in degree), and the problem is when plotting the diagram, some of the edges of the diagram are overlap (considering the projection into the sphere's surface). By looking the plot3, some voronoi vertices are actually joint "through the earth". This make sense because there are points in only one hemisphere and the surface triangulation "cut" the sphere right in region of the equator. So, I was thinking if there is any constraint (or trick) to add to the implementation in order to avoid surface triangulation that cut the sphere.
>
> phi = [51.77 43.79 36.18 54.72 60.22 43.46 57.38 36.10 43.81 40.50 36.97 27.94 38.52 43.79];
>
> theta = [257.77 234.56 219.73 339.49 335.62 272.83 348.08 219.92 228.03 3.11 25.17 15.39 28.72 318.43];
>
> I'd appreciate your help in this "special case".
>
> Many thanks
> Felipe

I'm not Damian, but here is something you could try:

function voronoisphere(randflag)
% voronoisphere(1) random
% voronoisphere(0) your data

if nargin<1
randflag = 1;
end

if randflag
n = 20;
xyz = randn(3,n);
xyz = bsxfun(@rdivide, xyz, sqrt(sum(xyz.^2,1)));
else
deg2rad = pi/180;

phi = [51.77 43.79 36.18 54.72 60.22 43.46 57.38 36.10 43.81 40.50 36.97 27.94 38.52 43.79];
theta = [257.77 234.56 219.73 339.49 335.62 272.83 348.08 219.92 228.03 3.11 25.17 15.39 28.72 318.43];

phi = deg2rad*phi;
theta = deg2rad*theta;
x = cos(theta).*sin(phi);
y = sin(theta).*sin(phi);
z = cos(phi);
xyz = [x; y; z];
end

T = convhull(xyz');
nt = size(T,1);
E1 = sort(T(:,[1 2]),2);
E2 = sort(T(:,[2 3]),2);
E3 = sort(T(:,[3 1]),2);

[~, ~, J] = unique([E1; E2; E3], 'rows');
n = accumarray(J,1);
if ~all(n==2)
error('Topology issue due to numerical precision')
end
id = accumarray(J, repmat((1:nt)',3,1), [], @(x) {x});
id = cat(2,id{:});
iA = id(1,:);
iB = id(2,:);

XYZ = reshape(xyz(:,T(:,[1 2 3 1])),3,[],4);
XYZ = permute(XYZ, [3 2 1]);

P = Center(xyz, T);

figure()
hold on
axis equal
plot3(XYZ(:,:,1),XYZ(:,:,2),XYZ(:,:,3),'ro');
plot3(P(1,:),P(2,:),P(3,:),'.');
for k=1:size(id,2)
A = P(:,iA(k));
B = P(:,iB(k));
G = Arc(A,B);
plot3(G(1,:),G(2,:),G(3,:));
end

end % voronoisphere

function G = Arc(A,B) % return an Arc between A and B
resolution = 2*pi/180; % 2 degrees
AxB = cross(A,B);
AdB = dot(A,B);
Ap = cross(AxB, A);
Ap = Ap/norm(Ap);
theta = atan2(sqrt(sum(AxB.^2,1)),AdB);
theta = linspace(0, theta, max(ceil(theta/resolution),2));
c = cos(theta);
s = sin(theta);
G = A*c + Ap*s;
end % Arc
function G = Arc(A,B) % return an Arc between A and B
resolution = 2*pi/180; % 2 degrees
AxB = cross(A,B);
AdB = dot(A,B);
Ap = cross(AxB, A);
Ap = Ap/norm(Ap);
theta = atan2(sqrt(sum(AxB.^2,1)),AdB);
theta = linspace(0, theta, max(ceil(theta/resolution),2));
c = cos(theta);
s = sin(theta);
G = A*c + Ap*s;
end % Arc

function P = Center(xyz, T)
XYZ = reshape(xyz(:,T),[3 size(T)]);
A = XYZ(:,:,1);
B = XYZ(:,:,2);
C = XYZ(:,:,3);
A = A-C;
B = B-C;
A2 = sum(A.^2,1);
B2 = sum(B.^2,1);
AxB = cross(A,B,1);
s = dot(AxB,C);
P = bsxfun(@times,A2,B) - bsxfun(@times,B2,A);
P = cross(P, AxB, 1);
P = bsxfun(@times,P,1./(2*sum(AxB.^2,1)));
P = C + P;
nP = sqrt(sum(P.^2,1));
P = bsxfun(@rdivide, P, nP);
s = sign(s);
P = bsxfun(@times, P, s);
end % Center

% Bruno

[Bruno Luong]

I modify the code so that voronoi cell contours are explicitly derived:

function voronoisphere(randflag)
% voronoisphere(1) random
% voronoisphere(0) your data

if nargin<1
randflag = 1;
end

if randflag

n = 100;

xyz = randn(3,n);
xyz = bsxfun(@rdivide, xyz, sqrt(sum(xyz.^2,1)));
else
deg2rad = pi/180;
phi = [51.77 43.79 36.18 54.72 60.22 43.46 57.38 36.10 43.81 40.50 36.97 27.94 38.52 43.79];
theta = [257.77 234.56 219.73 339.49 335.62 272.83 348.08 219.92 228.03 3.11 25.17 15.39 28.72 318.43];
phi = deg2rad*phi;
theta = deg2rad*theta;
x = cos(theta).*sin(phi);
y = sin(theta).*sin(phi);
z = cos(phi);
xyz = [x; y; z];
end

%%
npnts = size(xyz,2);
T = convhull(xyz.');
nt = size(T,1);
E = [T(:,[1 2]); T(:,[2 3]); T(:,[3 1])];
E = sort(E,2);
[~, ~, J] = unique(E, 'rows');

n = accumarray(J,1);
if ~all(n==2)
error('Topology issue due to numerical precision')
end

allids = repmat((1:nt)',[3 1]);
k = accumarray(J, (1:3*nt).', [], @(k) {k});
k = [k{:}];
vid = allids(k);
iA = vid(1,:);
iB = vid(2,:);

% each row is 2 cell ids of the edge
cellofedge = E(k(1,:),:); % ne x 2
ne = size(cellofedge,1);
edges = repmat((1:ne).',[2 1]);
edgeofcell = accumarray(cellofedge(:),edges, [], @(e) {e});

% Center of the circumscribed Delaunay triangles
P = Center(xyz, T);

% Build the geodesic arcs that link two vertices
nedges = size(vid,2);
edgearcs = cell(1,nedges);
for k=1:nedges

A = P(:,iA(k));
B = P(:,iB(k));
G = Arc(A,B);

edgearcs{k} = G;
end

% Build the contour of the voronoi cells
vidt = vid.';
vidts = sort(vidt,2);
voronoiboundary = cell(size(edgeofcell));
for k = 1:npnts
% ordering the edges
v = cycling_edge(edgeofcell{k}, vid);
v = oriented_edge(v, P, xyz(:,k));
[~, loc] = ismember(sort(v,2), vidts, 'rows');
m = length(loc);
X = edgearcs(loc);
flip = any(v(1:m,:)~=vidt(loc,:),2);
X(flip) = cellfun(@fliplr, X(flip), 'unif', false);
X = cellfun(@(x) x(:,1:end-1), X, 'unif', false);
voronoiboundary{k} = cat(2, X{:});
end

%% Graphic
f = figure(1);
clf(f);
set(f,'Renderer','zbuffer');
ax = axes('Parent', f);
hold(ax, 'on');
set(ax, 'Color', 'w');
plot3(ax, xyz(1,:),xyz(2,:),xyz(3,:),'w+');
clmap = lines();
ncl = size(clmap,1);
for k = 1:npnts
X = voronoiboundary{k};
cl = clmap(mod(k,ncl)+1,:);
fill3(X(1,:),X(2,:),X(3,:),cl,'Parent',ax,'EdgeColor','w');
end
axis(ax,'equal');
axis(ax,[-1 1 -1 1 -1 1]);

end % voronoisphere

%%

function G = Arc(A,B)
% return an Arc between A and B
resolution = 2*pi/180; % 2 degrees
AxB = cross(A,B);
AdB = dot(A,B);
Ap = cross(AxB, A);
Ap = Ap/norm(Ap);
theta = atan2(sqrt(sum(AxB.^2,1)),AdB);

npnts = max(ceil(theta/resolution),2);
theta = linspace(0, theta, npnts);
G = A*cos(theta) + Ap*sin(theta);
end % Arc

%%

function P = Center(xyz, T)

% Center of the circumscribed Delaunay triangles

XYZ = reshape(xyz(:,T),[3 size(T)]);
A = XYZ(:,:,1);
B = XYZ(:,:,2);
C = XYZ(:,:,3);
A = A-C;
B = B-C;
A2 = sum(A.^2,1);
B2 = sum(B.^2,1);
AxB = cross(A,B,1);

P = bsxfun(@times,A2,B) - bsxfun(@times,B2,A);
P = cross(P, AxB, 1);
P = bsxfun(@times,P,1./(2*sum(AxB.^2,1)));
P = C + P;
nP = sqrt(sum(P.^2,1));
P = bsxfun(@rdivide, P, nP);

s = dot(AxB,C);
P = bsxfun(@times, P, sign(s));
end % Center

%%
function v = cycling_edge(edges, vertexes)
% Chain the edges
u = vertexes(:,edges);
[~, ~, I] = unique(u);
I = reshape(I,2,[]);
j = repmat(1:size(u,2),[2 1]);
n = accumarray(I(:), 1);
if ~all(n == 2)

error('Topology issue due to numerical precision')
end

k = accumarray(I(:), j(:), [], @(x) {x});
k = [k{:}];
n = size(u, 2);
v = zeros([n 2]);
q = 1;
i = k(:,q);
p = i(1); % i(2) is also fine
for j = 1:n
i = I(:,p);
if i(1) == q
v(j,:) = u([1 2],p);
q = i(2);
else
v(j,:) = u([2 1],p);
q = i(1);
end
i = k(:,q);
if i(1) == p
p = i(2);
else
p = i(1);
end
end % for-loop
end % cycling_edge

%%
function v = oriented_edge(v, P, xyz)
% Orient the edges
Q = null(xyz.');
if det([xyz Q]) < 0
Q(:,2) = -Q(:,2);
end
E = P(:,v([1:end 1],1));
xy = Q'*E;
a = sum((xy(1,1:end-1)-xy(1,2:end)).*(xy(2,1:end-1)+xy(2,2:end)));
if a < 0
v = rot90(v,2);
end
end % oriented_edge

% Bruno

[Felipe Pedreros]

"Bruno Luong" <b.l...@fogale.findmycountry> wrote in message <kivj46$8vr$1...@newscl01ah.mathworks.com>...

> I modify the code so that voronoi cell contours are explicitly derived:

Hi Bruno,

Many thanks! Your code works pretty well. I'm now trying to get the area of each voronoi cell but I can't tell from the code where the voronoi vertices for each point are specified.

Regards.
Felipe

[Bruno Luong]

"Felipe Pedreros" <f_ped...@hotmail.com> wrote in message <kkdagd$bar$1...@newscl01ah.mathworks.com>...

>
> I can't tell from the code where the voronoi vertices for each point are specified.
>

I put the code with specified data output here:
http://www.mathworks.com/matlabcentral/fileexchange/40989-voronoi-sphere

Bruno