Thread Subject: Voronoi diagram on a sphere

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!

"Abel Brown" <brown.2179@osu.edu> wrote in message <hd0e8h$c9$1@fred.mathworks.com>...
> 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!

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

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.2179@osu.edu> wrote in message
news:hd0e8h$c9$1@fred.mathworks.com...
> 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" <b.luong@fogale.findmycountry> wrote in message <hd0hkr$pq1$1@fred.mathworks.com>...
> "Abel Brown" <brown.2179@osu.edu> wrote in message <hd0e8h$c9$1@fred.mathworks.com>...
> > 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!
>
> 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

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 ...

   

"Damian Sheehy" <Damian.Sheehy@mathworks.com> wrote in message <hd1a6r$rq0$1@fred.mathworks.com>...
> 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

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!!

"Abel Brown" <brown.2179@osu.edu> wrote in message
news:hd1e8s$gf6$1@fred.mathworks.com...
> "Damian Sheehy" <Damian.Sheehy@mathworks.com> wrote in message
> <hd1a6r$rq0$1@fred.mathworks.com>...
>> 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
>
> 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!!


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

> 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