Fill a region defined by the intersection of curves

[mottelet]

Hello,

A similar question has been asked on the mailing list and I think it may be easier to answer here: how do you fill a region of the plane defined by the intersection of several curves, given the coordinates of the intersection points ?

If you need the method to be general enough, you have to rely on the actual data of the Polyline. In the below example the points A,B,C are already known and we determine the limits of indices in each Polyline data by finding the closest data point (using the 1-norm):

function k = kinter(h,P)
    n = size(h.data,1);
    [_,k] = min(sum(abs(h.data-P(ones(n,1),:)),2));
endfunction

t=linspace(0,2*%pi,1024);
clf
c=2*cos(t);
s=2*sin(t);
h=plot(1+c,s,-0.5+c,sqrt(3)/2+s,-0.5+c,-sqrt(3)/2+s)
isoview on

A=(-1+sqrt(13))/2*[1 0];
B=(-1+sqrt(13))/2*[-0.5 sqrt(3)/2];
C=(-1+sqrt(13))/2*[-0.5 -sqrt(3)/2];

k1A = kinter(h(1),A);
k1B = kinter(h(1),B);
k3B = kinter(h(3),B);
k3C = kinter(h(3),C);
k2C = kinter(h(2),C);
k2A = kinter(h(2),A);

data = [
h(1).data(k1A:k1B-1,:)
h(3).data(k3B:k3C-1,:)
h(2).data(k2C:k2A,:)
]

xfpoly(data(:,1),data(:,2),addcolor([1 1 0]))

xstring(A(1),A(2),"A")
xstring(B(1),B(2),"B")
xstring(C(1),C(2),"C")

[heinznabielek]

Ingenious solution. I will try it out… Heinz

[arctica]

Hi Stephane,
The code is helpful, but have only partly managed to get it working with my data.

clear

function k = kinter(h,P)
n = size(h.data,1);
[_,k] = min(sum(abs(h.data-P(ones(n,1),:)),2));
endfunction

clf

curve1= csvRead(‘curve1.csv’);
curve2= csvRead(‘curve2.csv’);
circle= csvRead(‘ref_circle.csv’);
points=csvRead(‘points.csv’);

X1=curve1(1:,1);Y1=curve1(1:,2);
X2=curve2(1:,1),Y2=curve2(1:,2);
X=circle(1:,1),Y=circle(1:,2);

h=plot(X1,Y1, X2,Y2, X,Y);
isoview on

[xp yp]=(points(1),points(2));
[s1x s1y]=(points(3),points(4));
[s2x s2y]=(points(5),points(6));

A=[xp yp];B=[s1x s1y];C=[s2x s2y];

//original
k1A = kinter(h(1),A);
k1B = kinter(h(1),B);
k3B = kinter(h(3),B);
k3C = kinter(h(3),C);
k2C = kinter(h(2),C);
k2A = kinter(h(2),A);

data = [
h(1).data(k1A:k1B-1,
h(3).data(k3B:k3C-1,
h(2).data(k2C:k2A,
]

xfpoly(data(:,1),data(:,2),addcolor([1 1 0]))

xstring(A(1),A(2),“A”)
xstring(B(1),B(2),“B”)
xstring(C(1),C(2),“C”)

Hopefully the zip file uploaded to the mailing list ok, as that has all the data (not very big)

I am guessing it is an issue with the order of the curves?

Thanks
Lester

[arctica]

Solved the problem. The reference circle segment was ordered inversely so just need to flip the indices.

Updated code section:

k1A = kinter(curve1,A);//23
k1B = kinter(curve1,B);//181
k3B = kinter(circle,B);//29
k3C = kinter(circle,C);//3
k2C = kinter(curve2,C);//1
k2A = kinter(curve2,A);//168

circle_flip=circle(k3C:k3B,:);
circle=flipdim(circle_flip,1);

data = [
curve1(k1A:k1B-1,
circle(1:$-1,
curve2(k2C:k2A,
]

I located the issue by breaking the problem down into sections to locate the cause of the empty array. Tweaked the function so I could isolate the curve section causing issues.

function k = kinter(data,P)
n = size(data,1);
[_,k] = min(sum(abs(data-P(ones(n,1),:)),2));
endfunction

Thanks for the code suggestion

[mottelet]

@arctica

I am guessing it is an issue with the order of the curves?

Yes you were guessing right. Can you upload a screenshot of you example, with the filled region ?

S.

[arctica]

Hi Stephane,

I managed to get the plot to work and re-coded the original BB.m Matlab code for earthquake focal mechanisms.

However, the fill of the plot is not correct; not filling between the polyline area. I have verified that the the polyline is fine by extracting the points and checking in GIS that the (points->lines->polygon creates a closed object. Not sure why the fill is not doing what it should.

xfpoly(xx,yy,8) // the circle works fine as a solid fill
xfpoly(Y,X,4) // the fill seems spread

Any ideas?
Lester

[mottelet]

OK. If you look in my example I had to remove some points (see the two occurrences of minus one below)

data = [
h(1).data(k1A:k1B-1,:)
h(3).data(k3B:k3C-1,:)
h(2).data(k2C:k2A,:)
]

to remove this potential “butterfly” effect. Otherwise the last point of a polygon section can be in front (in the curvilinear abscissa sense) of the first point of the next section, athough it should be behind. This fix may be not needed if you discretize more finely your curves (that’s what I have noticed with my example).

S.