# Find max/min level of all intersections of a given contour curve with another set of contour curves

7 views (last 30 days)
Daniel on 7 Jul 2020
Commented: Matt J on 8 Jul 2020
I have two sets of contours, call them A(x,y) and B(x,y). For given contour level of A, I would like to find all the intersections with any levels of B and then find the max and min contour level of B within those intersections. I have the levels and x,y pairs from C2xyz on the File Exchange, though I think there are other ways. I looked into polyxpoly, but haven't quite figured it out and it appears it works with line segments. I'm not sure if the outputs of C2xyz qualify. My psuedo code is something like:
for i = 1:number of levels in A
find the intersections of given level of A with all levels of B and save the x,y coordinates and level of B
max level value of B of those intersections
min level value of B of those intersections
end

Daniel on 8 Jul 2020
I think this works for my purposes, at least. I'm sure there's a much more efficient way.
for i = 1:length(ccp) % These are levels of A
for j = 1:length(cct) % Theses are levels B
[tsrx{i,j},betax{i,j}] = polyxpoly(ccptsr{i},ccpbeta{i},ccttsr{j},cctbeta{j});
% tsr = x, beta = y, NOTE: this produces many empty cells
end
end
count = 1;
for i = 1:length(ccp)
for j = 1:length(cct)
if isempty(tsrx{i,j})==0 % if the cell from above is non-empty, proceed
tsrbetapairs{i}{count} = [tsrx{i,j};betax{i,j};i;j]; % extracts x, y, and the indices
% of the levels of A and B in i and j, respectively, and puts them in cells of i nested
% with non-empty cells of j
count = count+1;
end
end
count = 1;
end
for i = 1:length(ccp)
for j = 1:length(tsrbetapairs{i})
cctinds{i}(j) = tsrbetapairs{i}{j}(end); % levels of B are increasing, (end) is j, which
% corresponds to the level of B
end
maxct(i) = max(cctinds{i}); % maximum level of B that intersected with the ith level of A
minct(i) = min(cctinds{i}); % minimum level of B that intersected with the ith level of A
end

Matt J on 7 Jul 2020
Edited: Matt J on 7 Jul 2020
The problem you describe cannot be well-defined in terms of discrete contour samples, because there are infinite choices of functions A() and B() with those same samples. Instead, you should be solving this as a continuous optimization problem with nonlinear constraints A(x,y)<=a, B(x,y)<=b. One tool for doing so would be fmincon:

Daniel on 7 Jul 2020
Yes, in the simplest way, that's what I want. More specifically, I would like to repeat this to find the min and max levels of B that intersect each level of A. If it helps, I have used C2xyz from the file exchange to produce an arbitrary number of levels for each contour and the x,y coordinates that define each level. With this method, I don't have any control over the coarseness of the x,y coordinates that define the levels, though.
Matt J on 7 Jul 2020
I would like to repeat this to find the min and max levels of B that intersect each level of A
That's just a matter of looping the above code over 'a'. You also asked about the case where you optimize over all of B(x,y). In that case, just remove the nonlinear constraint on B by revising nonlcon as follows
function [c,ceq]=nonlcon(xy,A,B,a,b)
x=xy(1); y=xy(2);
ceq=[];
c=A(x,y)-a;
end
Now, please explain to me how it is that you don't have equations for the contours. How is your data generated? How do you even know your points are isocontours of anything if you don't know what function they come from?
Daniel on 7 Jul 2020
I suppose there's an equation for them somewhere, but I don't have it. They were generated by a simulation that spit out, for each x,y pair, what A and B are. If I had the equations, I would be pursuing a route more like you describe from the start.

Matt J on 7 Jul 2020
Edited: Matt J on 7 Jul 2020
A method more along the lines of what you outlined in your initial post, is to put the x,y pairs of the i-th level set of A into a cell array Apoints{i} and similarly the j-th contour points of B in Bpoints{j}. Then you would develop a matrix of polygon intersection areas as below.
IntersectionMatrix=nan(numlevels_A, numlevels_B); %pre-allocate
for i=1:numlevels_A
for j=1:numlevels_B
pgonA=polyshape(Apoints{i});
pgonB=polyshape(Bpoints{j});
IntersectionMatrix(i,j)=area(intersect(pgonA,pgonB));
end
end
Once you have IntersectionMatrix you can search along each row Intersection(i,:) until you find the first j such that Intersection (i,j)==0. Assuming the contour levels are in descending order of i,j, that first j will correspond to the B-contour where B is minimized over the i-th A contour. Naturally, this will be to within an accuracy which depends on how finely the contours are sampled. The row-wise search may be done in one line as follows.
[~,jmin]=max(Intersection<=0,[],2);

Daniel on 7 Jul 2020
I can't say that I completely understand how this is working. You're making polygons using the x,y coordinate pairs and then finding the area of their intersections. I don't see in the documentation for area how the intersections are the outputs. Regardless, when I try this, I get warnings from polyshape about the shapes being ill-defined. Could that be because my contours are not closed? Perhaps I shouldn't have suggested above in my example that they are. For reference, they look like the contours posted in this question.
Matt J on 8 Jul 2020
I don't see in the documentation for area how the intersections are the outputs.
The code uses the intersect() command as well.
when I try this, I get warnings from polyshape about the shapes being ill-defined. Could that be because my contours are not closed?
No, it is more likely due to the x,y points not being given in clockwise or counter-clockwise order.