Search All of the Math Forum:

Views expressed in these public forums are not endorsed by NCTM or The Math Forum.

Notice: We are no longer accepting new posts, but the forums will continue to be readable.

Topic: Voronoi diagram on a sphere
Replies: 11   Last Post: Apr 14, 2013 3:44 AM

 Messages: [ Previous | Next ]
 Abel Brown Posts: 62 Registered: 10/4/07
Re: Voronoi diagram on a sphere
Posted: Nov 6, 2009 10:10 AM

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

Date Subject Author
11/6/09 Abel Brown
11/6/09 Bruno Luong
11/6/09 Abel Brown
11/6/09 Damian Sheehy
11/6/09 Abel Brown
11/6/09 Damian Sheehy
11/6/09 Abel Brown
3/26/13 Bruno Luong
4/14/13 Felipe Pedreros
4/14/13 Bruno Luong