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

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

Answer 1

[Daniel]

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

Answer 2

[Matt J]

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:

https://www.mathworks.com/help/optim/ug/fmincon.html

[Daniel]

Thanks for the help. I don't really follow, though. Do you simply mean that the number of levels chosen is arbitrary? If so, I agree, but I have chosen them. With the ones I've chosen, I still want to do this. I understand that it will only be as accurate as the discretization of levels is small, but that's OK with me.

I don't know if it helps, but to make it more concrete, imagine one contour is concentric circles with the smallest having the highest value and the other is concentric triangles (or pick any shape that would allow for a clear min and max of intersecting levels). Given a level of the circles, which I have already define to my satisfaction, I want to find all the places where triangles intersect and then I want to find the min and max levels of the triangles at those intersections.

[Matt J]

Let's look at a more explicit example. An example of a function whos contours are concentric circles centered at the origin is A(x,y)=x^2+y^2. Rather than concentric triangles, I will choose for this example B(x,y)= 2*(x-1)^2 +3*(y-2)^2 whos contours are concentric ellipses centered at at x=1, y=2.

You say you have specific contours A(x,y)=a and B(x,y)=b which enclose some non-zero intersection area and you want to minimize (or maximize) B over this area. If this is indeed the kind of thing you want to do, then I claim that this is equivalent to the following optimization problem:

and can be readily solved using fmincon. For the particular A(x,y) and B(x,y) I have chosen above, the implementation will be:

A=@(xy) xy(1).^2 + xy(2).^2;
B=@(xy) 2*(xy(1)-1).^2 + 3*(xy(2)-2).^2;

xy_optimal=fmincon(B,[0,0],[],[],[],[],[],[],@(xy)nonlcon(xy,A,B,a,b))

function [c,ceq]=nonlcon(xy,A,B,a,b)

  ceq=[];
  
  c(1)=A(xy)-a;
  c(2)=B(xy)-b;

end

[Daniel]

I'm not a mathematician, so, sure, maybe this works. I have no idea what your code is doing, though. I have so many questions about it, I don't even know where to start. Which part of this produces the intersections? Which parts says which intersecting points are on the max or min level of B for a given level of A? And what about find the maxes?

I respect going about this another way if you see another way, but look at the level of question I posted and work with me a little, please.

[Matt J]

Which part of this produces the intersections?

The constraints A(x,y)<=a and B(x,y)<=b is the 2D region enclosed by the intersection.

Which parts says which intersecting points are on the max or min level of B for a given level of A?

The output xy_optimal will contain the x,y coordinates where the minimum B is achieved over this region.

And what about find the maxes?

The maxes can be found by minimizing the negative of B

[Daniel]

I don't doubt that this works, but I don't know how to get from what I have to using this. To start, I don't have the equations for my contours; they are produced by data that is loaded into the workspace. If I want to minimize over all B, would I just set b to something large? Then, to do this for specific levels of A, would I just change a and repeat?

Since I'm still not convinced we're entirely on the same page, I'll be more concrete with your example. Say that for a given circle, i.e., a given level of the circle contour, multiple levels of the ellipse contour intersect it. Will this find the min and max ellipse levels that intersect the given circle level? And then I can just repeat it for any specific circle level?

[Matt J]

Here is more complete code for the example above, together with plots to visualize the result. It finds the minimum ellipse level intersecting the circle A(x,y)=a=1. So, first, you confirm for me that this is the kind of thing that you want. Then, we will talk about how to generalize this to the case where you do not have an explicit formula for the contours.

close all

A=@(x,y) x.^2 + y.^2;
B=@(x,y) 2*(x-1).^2 + 3*(y-2).^2;

a=1;
b=9;

[xy_min,B_min]=fmincon(@(xy) B(xy(1),xy(2)) ,[0,0],[],[],[],[],[],[],@(xy)nonlcon(xy,A,B,a,b))

hold on
fa=fcontour(A,'LevelList',a,'LineWidth',2,'LineColor','blue');
fb=fcontour(B,'LevelList',b,'LineWidth',2,'LineColor','red');
fcontour(B,'LevelList',B_min,'LineWidth',2,'LineColor',fb.LineColor);
plot(xy_min(1),xy_min(2),'or','MarkerSize',8,'MarkerFaceColor','r');
axis equal
hold off

function [c,ceq]=nonlcon(xy,A,B,a,b)

  x=xy(1); y=xy(2);

  ceq=[];
  
  c(1)=A(x,y)-a;
  c(2)=B(x,y)-b;

end

[Daniel]

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]

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]

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.

Answer 3

[Matt J]

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]

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]

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.