Find two contours at different levels around the same peak

1 view (last 30 days)
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
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

Accepted Answer

Bjorn Gustavsson
Bjorn Gustavsson on 24 Oct 2022
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
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];
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.

More Answers (0)


Find more on Contour Plots in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by