Hello,
I have a set of X,Y data forming a closed polygon that is filled, however, Scilab is not configured to handle a self-intersecting polygon. Octave/Matlab can handle polygons that self-intersect.
How can this be sorted to arrive at the Octave/Matlab fill result, without breaking it apart etc? The module xfpoly is not available to look at, so guess it is hard-coded. The Octave fill.m function is available to view.
Thanks
Lester
Yes, but fill() calls patch which is a built-in. A priori, patch uses Qt backend for polygon filling. Scilab uses joGL for the rendering, but the polygon is triangulated by a built-in Scilab C++ routine (in graphic_objects module), which does not consider the case of self-intersecting polygons.
S.
Hi Stephane,
Thanks for explaining how the functions differ and work.
I’m guessing the route is to get the coordinates of the self-intersection (is there a code for that in Scilab?), and then close off the two separate polygons.
Is there a way to check the order of the vertices since they have to be ordered in a clockwise or anticlockwise direction?
Thanks
Lester
Here is a possible solution, using fsolve to find the intersection point between two piecewise affine parametric curves (x_1(n),y_1(n)) and (x_2(n),y_2(n)). The equations to fullfill are
where n_1 \in [1,N_1] and n_2 \in [1,N_2] and N_1,N_2 are the number of edges of each curve. Integer values of the parameter match the edges and between these, points continuously move on the segments. In the below code we use Scilab’s interpln() to compute the coordinates. After fsolve() has converged we get non-integer values n_1,n_2 and we use floor() and ceil() functions to identify the actual segments where the crossing occurs. Then both splitted parts of the curves are joined in the the correct order (I have added arrows in the plot that helps to understand the chosen order) to fill the polygons.
function eq=f(n)
eq = [
interpln([1:N1;x1],n(1))-interpln([1:N2;x2],n(2))
interpln([1:N1;y1],n(1))-interpln([1:N2;y2],n(2))
]
endfunction
//build two parametric curves (x1,y1) and (x2,y2)
t1=linspace(0,%pi/3,15);
x1=sin(t1);
y1=cos(t1);
t2=linspace(%pi/3,0,10);
x2=1.1-sin(t2)
y2=0.2+cos(t2)
N1 = length(x1);
N2 = length(x2);
clf
plot(x1,y1,x2,y2)
// find the intersection point (xi,yi) using fsolve
[n,f]=fsolve([n1/2;n2/2],f);
xi = interpln([1:N1;x1],n(1));
yi = interpln([1:N1;y1],n(1));
// left polygon
xp1 = [x1(1:floor(n(1))) xi x2(floor(n(2)):-1:1)];
yp1 = [y1(1:floor(n(1))) yi y2(floor(n(2)):-1:1)];
xfpoly(xp1,yp1,color("yellow"))
// right polygon
xp2 = [xi x2(ceil(n(2)):$) x1($:-1:ceil(n(1)))];
yp2 = [yi y2(ceil(n(2)):$) y1($:-1:ceil(n(1)))];
xfpoly(xp2,yp2,color("orange"))
plot(x1,y1,"-o",x2,y2,"-o")
gce().children.polyline_style=4
plot(xi,yi,'or')
Hi Stephane,
Thanks for the example code. Could not manage to get it to work with the data I am testing, but had an intersection function based on line intersection between nodes to find the crossing point; it is quite slow but gets there in the end (not perfect).
Having found the point where the polyline crosses itself (xC,yC). I added this point into the original X and Y values so that all points are now available. The image shows the order of the vertices as a continuous polyline. In terms of indices, start=1, xC/yC=21 and end=959. Now, with some trial-and-error, the polyline crosses back through the intersection at about node=492. I could add an extra node at the end to close the polyline completely, but this is not necessary with xfpoly.
In essence, the intersection point has two indices, one at 21 and the other around 492 approximately. Is there a way to count along the indices to find the coincident point?
Any suggestions would be great.
Lester
Can you share the code which constructs the self intersecting polyline above ?
S.
Hi Stephane,
I have e-mailed you the code and associated functions.
Thanks
Lester
Hello Lester,
Here is a brute force (complexity is O(n^2)) solution under the form of a function splitting a polygon if having a single self-intersection (but the function could be called recursively on both parts):
function splitpoly(h)
function [flag, d]=intersect(p1, q1, p2, q2)
eps2 = %eps*%eps;
d = [det([q1-p1 p2-p1])
det([q1-p1 q2-p1])
det([q2-p2 p1-p2])
det([q2-p2 q1-p2])];
flag = (d(1)*d(2) < eps2) & (d(3)*d(4) < eps2)
endfunction
flag = %f;
P = h.data';
P1 = [];
P2 = [];
n = size(P,2);
for i=1:n-1
for j = i+2:n-1
[flag,d] = intersect(P(:,i),P(:,i+1),P(:,j),P(:,j+1));
if flag
t = d(3)/(d(2)-d(1));
Q = P(:,i)+t*(P(:,i+1)-P(:,i));
break
end
end
if flag
h2 = copy(h);
h.data = [P(:,1:i) Q P(:,j+1:$)]';
h2.data = [Q P(:,i+1:j)]';
break
end
end
end
After executing your main script BB_test.sce (where I removed your actual code searching the intersection point) then calling
--> tic();splitpoly(gce());toc()
ans =
0.8591190
I get the following output:
S.
That seems to work to solve the fill issue. It should work for most scenarios.
Thanks for looking into this.
Lester
