voronoi within a pre-defined area

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

[Walter Roberson]

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

Answer 1

[Preetham Manjunatha]

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.

Answer 2

[Walter Roberson]

http://www.mathworks.com/help/techdoc/ref/delaunaytriclass.html and the associated Constraints property?

http://www.mathworks.com/matlabcentral/answers/20052-meshing-around-the-outside-of-a-circle