voronoi within a pre-defined area
10 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Hi I asked a question couple of weeks ago concerning the generation Voronoi polygon within a closed area.
but I’m still stuck.
Basically I have a series of 50 points, I start with the three first and increase the number of point by one in each loop. For each population of points I generate the area of interest based by increasing the convex hull, and then within that area I’m interested by the area of each polygon created by the Voronoi analysis.
data %first 10 points
299337.7 5924043.5
301012.5 5919417
300634.6 5927108.6
309884.5 5909105.8
301561.3 5905336
305310 5907794
309769.3 5919926.9
300930.3 5915819.6
300366.4 5904622.8
305082.6 5917688.9
My problem is that I don’t know how to limit my lines to the boundary of my area of interest, sometimes the intersection between the lines is inside or outside the area of interest. Therefore how to tell my code to find every intersection in the surrounding one point, taking in account that I need to eliminate some lines which are not suitable for the generation of polygons.
I appreciate any help
close all
clear all
import_textfile('SequenceBoggust_50volcanes_11_01_2012.txt')
data=SequenceBoggust_50volcanes_11_01_2012;
for i=0:47
clear x y k
figure
x=data(1:3+i,1); %coordinates x number of points
y=data(1:3+i,2); %coordinates y number of points
plot(x,y,'bo')
%--------------------------------------------------------------------------
%Search for the points on the convex hull
convexhull=convhull(x,y);
%plot(data(k,1),data(k,2),'r-',x,y,'b+')
axis equal
hold on
%Calculate centroid polygon
[geom,iner,cpmo]=polygeom(data(convexhull,1),data(convexhull,2));
%plot(geom(1,2),geom(1,3),'r*')
%--------------------------------------------------------------------------
%Calculate mean distance between the volcanoes/points
X=[x y];
dt = delaunayn(X);% Triangulation
triplot(dt,x,y)
%voronoi(data(:,1),data(:,2))
% e is an Nx2 matrix of indices into p, e = [n1,n2; etc]
e = [dt(:,[1,2]); dt(:,[2,3]); dt(:,[3,1])]; % Non-unique
e = unique(e,'rows'); % Unique list
% Edge length
h = sqrt( sum( (X(e(:,1),:)-X(e(:,2),:)).^2 ,2) );
% Get average around each vertex
h_ave = 0*X(:,1);
count = 0*X(:,1);
for j = 1:size(e,1)
% Vertices for kth edge
n1 = e(j,1); n2 = e(j,2);
% Add to vertices
h_ave(n1) = h_ave(n1)+h(j); count(n1) = count(n1)+1;
h_ave(n2) = h_ave(n2)+h(j); count(n2) = count(n2)+1;
end
h_ave = h_ave./count;
mean_dist=mean(h_ave);
%--------------------------------------------------------------------------
%Create larger convex hull for points at further distance (=mean distance)
% -> increase size of the field.
bound=zeros(length(convexhull),2);
for k=1:length(data(convexhull))
%create new points
if data(convexhull(k),1)<geom(2)
a=(geom(3)-data(convexhull(k),2))/(geom(2)-data(convexhull(k),1));
b=geom(3)-a*geom(2);
dist=sqrt((data(convexhull(k),2)-geom(3))^2 +(geom(2)-data(convexhull(k),1))^2);
dist1=dist+mean_dist;
delta=acos((geom(2)-data(convexhull(k),1))/dist);
dist2=cos(delta)*dist1;
bound(k,1)=geom(2)-dist2;
bound(k,2)=a*bound(k,1)+b;
elseif data(convexhull(k),1)>geom(2)
a=(data(convexhull(k),2)-geom(3))/(data(convexhull(k),1)-geom(2));
b=geom(3)-a*geom(2);
dist=sqrt((geom(3)-data(convexhull(k),2))^2 + (geom(2)-data(convexhull(k),1))^2);
dist1=dist+mean_dist;
delta=acos((data(convexhull(k),1)-geom(2))/dist);
dist2=cos(delta)*dist1;
bound(k,1)=geom(2)+dist2;
bound(k,2)=a*bound(k,1)+b;
end
end
plot(bound(:,1),bound(:,2),'r-') %plot extended field
axis equal
hold on
%filename = ['bound' num2str(i)];
%xlswrite(filename, bound);
%Calculate the linear equations of the boundaries of the extended field
for l=1:size(bound,1)-1
aboundedge(l) =(bound(l+1,2)-bound(l,2)) / (bound(l+1,1)-bound(l,1));
bboundedge(l) = bound(l,2) - aboundedge(l)*bound(l,1);
end
%--------------------------------------------------------------------------
%Calculate the Voronoi lines and the area of Voronoi polygons within
% the extended field
for m=1:1:size(X,1)
Coord{m}=[x(m) y(m)];
end
temp=1;
for n=2:1:size(X,1)
for o=1:1:n-1
%Calcuates the midpoint of the line linking 2 points.
midy(temp) = (Coord{n}(2)+Coord{o}(2))/2;
midx(temp) = (Coord{n}(1)+Coord{o}(1))/2;
%Calculates the equation of the Voronoi edges.
a(temp) = (Coord{n}(2)-Coord{o}(2))/(Coord{n}(1)-Coord{o}(1));
abis(temp) = -1/a(temp); %Calculate a of y=ax+b
bbis(temp) = midy(temp)-midx(temp)*abis(temp); %Calculate b of y=ax+b
xbis{temp}=[2.85e5;3.25e5];
ybis{temp}=abis(temp)*xbis{temp}+bbis(temp);
plot(xbis{temp},ybis{temp},'b-',xbis{temp},ybis{temp},'r+')
%Points where perpendicular bisector crosses the extended
%field boundary
nbintersection=0;
for p=1:length(aboundedge)
P= [bboundedge(p),bbis(temp)]/[-aboundedge(p),-abis(temp);1,1]; %p = [x,y];
[IN ON]= inpolygon(P(1),P(2),bound(:,1),bound(:,2));
if ON==1
nbintersection=nbintersection+1;
Intersect(nbintersection,:)=P;
plot(Intersect(nbintersection,1),Intersect(nbintersection,2),'g*')
end
end
EqBis{temp}=[abis(temp) bbis(temp)]; %a = EqBis{}(1); b=EqBis{}(2)
IntersectBisBound{temp}=[Intersect(:,1) Intersect(:,2)];
temp=temp+1;
end
end
end
1 commentaire
Walter Roberson
le 23 Jan 2012
Previous posting was http://www.mathworks.com/matlabcentral/answers/25896-voronoi-within-a-pre-defined-area
Réponses (2)
Preetham Manjunatha
le 8 Fév 2022
Modifié(e) : Preetham Manjunatha
le 8 Fév 2022
Here is the link function to clip the extending edges of the Voronoi Diagram for rectangular or square region. Rigorously tested on the random points, this function can process an input data set of 2000 seed points in 2D in about 0.015 seconds on average.
0 commentaires
Walter Roberson
le 23 Jan 2012
http://www.mathworks.com/help/techdoc/ref/delaunaytriclass.html and the associated Constraints property?
0 commentaires
Voir également
Catégories
En savoir plus sur Voronoi Diagram dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!