Find two contours at different levels around the same peak

Hi.

I would like to find two contours at different levels around the same lower peak and only keep those where this is true. I also want to reference two contours at different heights to each other. I assume the contour function is kind of random.

Let's say I have two different thresholds for two different levels. How can I find all the contours for threshold2 where also contours for threshold1 exist.

Threshold1 is lower than threshold2.

In the code example below I would only be interested in the two contours around the lower peak and I want to disregard or reject the other contour where no threshold1 contour exists.

In the example I only have two lower peaks but If I have many more lower peaks I would like to reference the contours to each other so that I know which contour at threshold2 is drawn around the same lower peak as for threshold1.

[X,Y,Z] = peaks(1000);
figure(1),hold on
view(2);
colorbar('vert');
s = surf(X,Y,Z,'Facealpha',0.75,'Edgecolor','none');
% threshold = -0.3*max(Z,[],'all');
threshold1 = [-6, -6];
threshold2 = [-3, -3];
[c1,h1] = contour(X,Y,Z,[threshold1 threshold1],'-k','Linewidth',5);
[c2,h2] = contour(X,Y,Z,[threshold2 threshold2],'-k','Linewidth',5);
hold off

Thanks.

Answer 1

[Bjorn Gustavsson]

You could (QD algorithming-from-the-hip) use the contour-coordinates in c1 and c2. Typically I find it easiest to break them contours up into cell-arrays with the x and y coordinates for each contour in an element. Something like this:

function [xC,yC,zC,cAll] = split_contour(C)
iC = 1;
idxC = 1;
while iC < size(C,2)
  nP         = C(2,iC);        % number of points in current contour
  xC{idxC}   = C(1,iC+(1:nP)); % x coordinates of current contour
  yC{idxC}   = C(2,iC+(1:nP)); % y coordinates of current contour
  cC{idxC}   = C(1,iC)*ones(size(yC{idxC}));
  cAll(idxC) = C(1,iC);
  zC{idxC}   = F(xC{idxC},yC{idxC}); % z-level of contour
  iC         = iC+nP+1;              % Start-point of next contour
  idxC       = idxC + 1;             % next contourline index
end

Then you only need to check if one point on any c1-contour is inside a c2-contour, otherwise scrap it. You could do that with inpolygon. Perhaps something like:

c2 = contour(peaks(123),[1 1]*level2);
c1 = contour(peaks(123),[1 1]*level2);
[xC1,yC1] = split_contour(c1);
[xC2,yC2] = split_contour(c2);
idx2keep = [];
for i2 = 1:numel(xC2)
    % If you have a vast number of contours, then you should remove the
    % elements of the c1-contours that you've found inside one c2-contour
    % since it cannot be inside 2 c2-contours. Here I dont bother
    for i1 = 1:numel(xC1)
        if  inpolygon(xC1{i1}(1),yC1{i1}(1),xC2{i2},yC2{i2});
            idx2keep = [idx2keep,i2];
        end
    end
end

Both code-snippets written in the editor and are completely untested. Check the help for contour if the first part doesn't work and fix the function. If the second part doesn't work test it line-by-line and fix the code. The algorithm should be sensible - unless you have some contours escaping out around a corner, then I can smell problems. That multiple c1-contours could be inside each c2-contours I think would be handled properly.

HTH