voronoi within a pre-defined area

10 ビュー (過去 30 日間)
Nicolas
Nicolas 2012 年 1 月 23 日
編集済み: Preetham Manjunatha 2022 年 2 月 8 日
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 件のコメント
Walter Roberson
Walter Roberson 2012 年 1 月 23 日
Previous posting was http://www.mathworks.com/matlabcentral/answers/25896-voronoi-within-a-pre-defined-area

サインインしてコメントする。

回答 (2 件)

Preetham Manjunatha
Preetham Manjunatha 2022 年 2 月 8 日
編集済み: Preetham Manjunatha 2022 年 2 月 8 日
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.

Walter Roberson
Walter Roberson 2012 年 1 月 23 日

カテゴリ

Help Center および File ExchangeVoronoi Diagram についてさらに検索

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by