Intersection of voronoi edges and a rectangle

This code generates a Voronoi plot with n random points, and a rectangle containing each closed Voronoi cell. How can I find the intersection points between the rectangle, and the voronoi edges (edges generated in figure 2)?

clc
clear all

n= 10;
x=10*rand(1,n);
y=10*rand(1,n);
h=voronoi(x,y);

[vx,vy] =voronoi(x,y);
[v,c] = voronoin([x(:) y(:)]);
close all
plot(x,y,'r+',vx,vy,'b-'); 
rectangle('Position',[0,0,12,12]);
axis equal
figure(2)
for j=1:length(vx(1,:))
    line([vx(1,j) vx(2,j)],[vy(1,j) vy(2,j)])
    rectangle('Position',[0,0,12,12]);
end
axis equal

[THERESA NYOROH]

Hello everyone, Mine is a question I need the syntax for constructing a voronoi diagram showing the values of its vertices as well as calculating the Euclidean distance and the voronoi circle.pls, I'm really in need of this code with an explanation of what each line of code does. Thanks

[Bruno Lopes]

You can use this code: https://it.mathworks.com/matlabcentral/fileexchange/34428-voronoilimit-varargin

It is the best I found on the internet.

Answer 1

[Geoff Hayes]

Sean - here is one way to calculate the intersection between the Voronoi edges (generated in the second figure) and the rectangle. For each edge that we draw on the plot, we can determine the slope of that edge/line given its two points

 m = (vy(2,j)-vy(1,j))/(vx(2,j)-vx(1,j));

The equation of the line for these two points and slope (given by usual equation y-y1=m(x-x1), where (x1,y1) is a point on the line, m is the slope) can be written as an anonymous function

 yEqn = @(x)m*(x-vx(1,j))+vy(1,j);

We can also write the above in terms of x

 xEqn = @(y)(y-vy(1,j)+m*vx(1,j))/m;

Since we want the intersection of the Voronoi edge and a rectangle, we can compute the intersection of this edge with one of the four sides that make up the rectangle given by the equations: x=0, x=12, y=0, and y=12. Note that for x=0 and x=12, 0<=y<=12; and for y=0 and y=12, 0<=x<=12.

An edge intersects with either x=0 or x=12 if the two points of the edge straddle either of these two lines. So we can do the following

 % for equations x=0 and x=12
 for x=[0 12]
     % check for straddling of x
     if (vx(1,j)<=x && vx(2,j)>=x) || (vx(1,j)>=x && vx(2,j)<=x)

         y = yEqn(x);

         % if y is between 0 and 12, then point of intersection is on 
         % line x=0 or x=12
         if y>=0 && y<=12
             plot(x,y,'r*');
             % skip to next edge
             continue;
         end
     end
  end

If the edge does not intersect with either of the x equations, then we try for y=0 or y=12 in the same manner

 % for equations y=0 and y=12
 for y=[0 12]
     % check for straddling of y
     if (vy(1,j)<=y && vy(2,j)>=y) || (vy(1,j)>=y && vy(2,j)<=y)

         x = xEqn(y);

         % if x is between 0 and 12, then point of intersection is on 
         % line y=0 or y=12
         if x>=0 && x<=12
             plot(x,y,'r*');
             % skip to next edge
             continue;
         end
     end
 end

Try the above code and see what happens. Note that it may need to be adjusted to handle the cases where the Voronoi edges are vertical or horizontal lines.

[Sean]

Cheers, follow up question: I tried using polyxpoly for the Voronoi plot, and the rectangle (created differently), but I am having more point intersections than there are supposed to be. I have the positions marked in the following code:

clc
clear all

n= 10;
x=10*rand(1,n);
y=10*rand(1,n);
h=voronoi(x,y);

[vx,vy] =voronoi(x,y);
[v,c] = voronoin([x(:) y(:)]);
close all
plot(x,y,'r+',vx,vy,'b-'); 
axis equal
xbox = [0 0 12 12 0];
ybox = [0 12 12 0 0];
mapshow(xbox,ybox,'DisplayType','line','LineStyle','-')

figure(2)
for j=1:length(vx(1,:))
    line([vx(1,j) vx(2,j)],[vy(1,j) vy(2,j)])

end
axis equal

xbox = [0 0 12 12 0];
ybox = [0 12 12 0 0];
mapshow(xbox,ybox,'DisplayType','line','LineStyle','-')

[xi, yi, ii] = polyxpoly(vx, vy, xbox, ybox);
mapshow(xi,yi,'DisplayType','point','Marker','o')

[Geoff Hayes]

I don't have the Mapping Toolbox so can't run your above code. But looking at the description of http://www.mathworks.com/help/map/ref/polyxpoly.html

[xi,yi] = polyxpoly(x1,y1,x2,y2) returns the intersection points of two polylines in a planar, Cartesian system. x1 and y1 are vectors containing the x- and y-coordinates of the vertices in the first polyline..

Your vx and vy are matrices (rather than vectors) and do not connect to form a line. Only pairs of your vertices do that.

You may have to do something like

 for j=1:length(vx(1,:))
    line([vx(1,j) vx(2,j)],[vy(1,j) vy(2,j)]);
    [xi, yi, ii] = ...
      polyxpoly([vx(1,j) vx(2,j)], [vy(1,j) vy(2,j)], xbox, ybox);
    if ~isempty(xi)
        plot(xi,yi,'r*');
    end
 end

This might work - checking each edge against the rectangle and only plotting an intersection point if the returned result is non-empty (I don't know if the result is empty if there is no intersection - you would have to try that).

[Sean]

Much appreciated.

[Sean]

For anyone else referring to this, removing the if loop with ~isempty(xi) did the trick (for the example in this question).