voronoi within a pre-defined area

Question

Nicolas 2012-1-23

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/26878-voronoi-within-a-pre-defined-area

编辑： 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 个评论
显示 -1更早的评论隐藏 -1更早的评论

Walter Roberson 2012-1-23

Previous posting was http://www.mathworks.com/matlabcentral/answers/25896-voronoi-within-a-pre-defined-area

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Preetham Manjunatha 2022-2-8

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/26878-voronoi-within-a-pre-defined-area#answer_891350

编辑：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.