Path: news.mathworks.com!not-for-mail
From: "Damian Sheehy" <Damian.Sheehy@mathworks.com>
Newsgroups: comp.soft-sys.matlab
Subject: Re: Voronoi diagram on a sphere
Date: Fri, 6 Nov 2009 10:54:47 -0500
Organization: The MathWorks, Inc.
Lines: 121
Message-ID: <hd1gs9$5tm$1@fred.mathworks.com>
References: <hd0e8h$c9$1@fred.mathworks.com> <hd1a6r$rq0$1@fred.mathworks.com> <hd1e8s$gf6$1@fred.mathworks.com>
Reply-To: "Damian Sheehy" <Damian.Sheehy@mathworks.com>
NNTP-Posting-Host: sheehyd.dhcp.mathworks.com
X-Trace: fred.mathworks.com 1257522889 6070 172.31.44.102 (6 Nov 2009 15:54:49 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Fri, 6 Nov 2009 15:54:49 +0000 (UTC)
X-Priority: 3
X-MSMail-Priority: Normal
X-Newsreader: Microsoft Outlook Express 6.00.2900.5843
X-RFC2646: Format=Flowed; Original
X-MimeOLE: Produced By Microsoft MimeOLE V6.00.2900.5579
Xref: news.mathworks.com comp.soft-sys.matlab:583045



"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