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
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@inrae.fr 0020 % Created: 2010-11-17, using Matlab 7.9.0.529 (R2009b) 0021 % Copyright 2010-2024 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,:);