SPHERICALVORONOIDOMAIN Compute a spherical voronoi domain. POLY = sphericalVoronoiDomain(GERM, NEIGHBORS) GERM is a 1-by-3 row vector representing cartesian coordinates of a point on the unit sphere (in X, Y Z order) NEIGHBORS is a N-by-3 array representing cartesian coordinates of the germ neighbors. It is expected that NEIGHBORS contains only neighbors that effectively contribute to the voronoi domain. Example sphericalVoronoiDomain See also drawSphericalPolygon ------ Author: David Legland e-mail: david.legland@grignon.inra.fr Created: 2010-11-17, using Matlab 7.9.0.529 (R2009b) Copyright 2010 INRA - Cepia Software Platform.
0001 function poly = sphericalVoronoiDomain(refPoint, neighbors) 0002 %SPHERICALVORONOIDOMAIN Compute a spherical voronoi domain. 0003 % 0004 % POLY = sphericalVoronoiDomain(GERM, NEIGHBORS) 0005 % GERM is a 1-by-3 row vector representing cartesian coordinates of a 0006 % point on the unit sphere (in X, Y Z order) 0007 % NEIGHBORS is a N-by-3 array representing cartesian coordinates of the 0008 % germ neighbors. It is expected that NEIGHBORS contains only neighbors 0009 % that effectively contribute to the voronoi domain. 0010 % 0011 % Example 0012 % sphericalVoronoiDomain 0013 % 0014 % See also 0015 % drawSphericalPolygon 0016 % 0017 % ------ 0018 % Author: David Legland 0019 % e-mail: david.legland@grignon.inra.fr 0020 % Created: 2010-11-17, using Matlab 7.9.0.529 (R2009b) 0021 % Copyright 2010 INRA - Cepia Software Platform. 0022 0023 % reference sphere 0024 sphere = [0 0 0 1]; 0025 0026 % number of neigbors, and number of sides of the domain 0027 nbSides = size(neighbors, 1); 0028 0029 % compute planes containing separating circles 0030 planes = zeros(nbSides, 9); 0031 for i = 1:nbSides 0032 planes(i,1:9) = normalizePlane(medianPlane(refPoint, neighbors(i,:))); 0033 end 0034 0035 % allocate memory 0036 lines = zeros(nbSides, 6); 0037 intersects = zeros(2 * nbSides, 3); 0038 0039 % compute circle-circle intersections 0040 for i = 1:nbSides 0041 ind2 = mod(i, nbSides) + 1; 0042 lines(i,1:6) = intersectPlanes(planes(i,:), planes(ind2,:)); 0043 intersects(2*i-1:2*i,1:3) = intersectLineSphere(lines(i,:), sphere); 0044 end 0045 0046 % keep only points in the same direction than refPoint 0047 ind = dot(intersects, repmat(refPoint, [2 * nbSides 1]), 2) > 0; 0048 poly = intersects(ind,:);