closed countour around a given point

5 views (last 30 days)
Hi and I'd like to say thanks in advance.
I am processing some gridded data using contour/countourc commands and I am trying to do 3 things.
  1. find all the closed countours that are centered and suround a given x,y cordinate
  2. find the outermost closed countour from the above list
  3. find the maximum, min and average distance of the above outermost contour from the center x,y cordinate (countours are not always circular and may be of irregular shape in my data).
I started with a sample data (from matlab called peaks) to build my code and impliment on my actual data but I am still stuck on 1 and 2.
(image attaced) I want my code to detect contour levels 18, 16, 14, 12 which are closed and surrounding my center point (blue dot) but instead it is picking 8, 10, 14, 16, 18. How can i modify the code to only detect the prefered contour levels (18,16,14,12) from which i can find the outermost contour and the 3 distances from the center. Note, level 18 is behind the center dot. Also attached countour z data... just incase
%___________________________________________________________________________
% Specify the center point and search radius
centerPoint = [25, 37.5];
searchRadius = 20;
% get and plot contour data
Z = peaks+10; % make colorbar positive. using built in sample data from matlab (peaks)
[M,h] = contour(Z);
clabel(M, h);
colorbar
% plot centerPoint
hold on
scatter(25, 37.5,"filled")
hold off
% crop portion of search area for easy visual
xlim([25 - searchRadius 25 + searchRadius])
ylim([37.5 - searchRadius 37.5 + searchRadius])
%___________________________________________________________________________
% Calculate contours
contourData = contourc(Z);
% Parse contour data
idx = 1;
while idx < size(contourData, 2)
contourLevel = contourData(1, idx); % Contour level
numPoints = contourData(2, idx); % Number of points in the contour
contourPoints = contourData(:, idx+1:idx+numPoints); % Contour points (x, y coordinates)
% Check if contour is closed (first and last points are the same)
isClosed = isequal(contourPoints(:,1), contourPoints(:,end));
% Check if contour surrounds the specified point
if isClosed
% Calculate distance of contour points to center point
distToCenter = sqrt((contourPoints(1, :) - centerPoint(1)).^2 + (contourPoints(2, :) - centerPoint(2)).^2);
% Check if all points are within search radius
if all(distToCenter <= searchRadius)
disp(['Contour Level: ', num2str(contourLevel)]);
disp(['Contour is closed and surrounds point (', num2str(centerPoint(1)), ', ', num2str(centerPoint(2)), ')']);
disp('Contour Points:');
%disp(contourPoints);
disp('---------------------------------------------');
end
end
% Move to the next contour
idx = idx + numPoints + 1;
end
clear Z M h contourData contourLevel numPoints contourPoints ...
centerPoint searchRadius contourLevels isClosed distToCenter idx;

Accepted Answer

Star Strider
Star Strider on 1 May 2025
Edited: Star Strider on 1 May 2025
This will detect and plot all the closed contours (in addition to a second one at level 12 that does not actually meet the closed criterion, although it is eliminated in the final step).
After isolating the candidate closed contours, this uses the inpolygon function to check to be sure that the desired closed contours contain the reference point.
Try something like this —
% Specify the center point and search radius
centerPoint = [25, 37.5];
searchRadius = 20;
% get and plot contour data
Z = peaks+10; % make colorbar positive. using built in sample data from matlab (peaks)
figure
[M,h] = contour(Z);
clabel(M, h);
colorbar
Levelv = h.LevelList
Levelv = 1×8
4 6 8 10 12 14 16 18
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
% plot centerPoint
hold on
scatter(25, 37.5,"filled")
hold off
% crop portion of search area for easy visual
xlim([25 - searchRadius 25 + searchRadius])
ylim([37.5 - searchRadius 37.5 + searchRadius])
xl = xlim;
yl = ylim;
% Calculate contours
contourData = contourc(Z);
% Check = ismember(contourData, M)
for k = 1:numel(Levelv)
ixp1 = find(M(1,:) == Levelv(k)); % Find Start Of Each Level
ixp2 = rem(M(2,ixp1),1) == 0; % Length Must Be An Integer
idx{k} = ixp1(ixp2);
len{k} = M(2,idx{k});
Levelc{k} = M(1,idx{k});
end
% LV = [Levelc{:}]
% IX = [idx{:}]
% LN = [len{:}]
for k = 1:numel(idx)
% Levelc{k}
checkClosed = M(:,idx{k}+1) == M(:,idx{k}+len{k}); % Test Beginning & End For Equality
Lclosed = checkClosed(1,:) & checkClosed(2,:); % Logical 'and'
isClosed{k} = [idx{k}(Lclosed); len{k}(Lclosed); Levelc{k}(Lclosed)]; % Index & Length Of Closed Contours
end
ICM = [isClosed{:}]
ICM = 3×10
1 17 53 109 281 291 323 427 473 505 15 35 55 31 9 31 103 45 31 5 4 6 8 8 10 12 12 14 16 18
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
refpt = [25, 37.5];
figure
hold on
k2 = 0;
for k = 1:size(ICM,2)
% k
idxrng = ICM(1,k)+1 : ICM(1,k)+ICM(2,k);
xv = M(1,idxrng);
yv = M(2,idxrng);
if inpolygon(refpt(1), refpt(2), xv, yv)
k2 = k2 + 1;
keepContour(:,k2) = ICM(:,k);
plot(xv, yv)
end
end
scatter(25, 37.5,"filled")
hold off
xlim(xl)
ylim(yl)
keepContour
keepContour = 3×4
323 427 473 505 103 45 31 5 12 14 16 18
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
MeetCriteria = array2table(keepContour.', VariableNames=["Start Index","Length","Level"])
MeetCriteria = 4x3 table
Start Index Length Level ___________ ______ _____ 323 103 12 427 45 14 473 31 16 505 5 18
I leave the distances and other calculations to you.
EDIT — (1 May 2025 at 05:59)
Corrected typographical errors.
EDIT — (1 May 2025 at 14:49)
Added 'MeetCriteria' table array.
.

More Answers (2)

Image Analyst
Image Analyst on 1 May 2025
Just treat your image as an image and use image processing. So threshold it at whatever level you want to get a binary image. Then use imreconstruct to find out which contour(s) (blobs) contain the location/coordinate of interest.
mask = Z > 14; % whatever value you want.
marker = false(size(mask));
marker(row, column) = true; % Set the coordinate you want the blob to contain.
containingBlob = imreconstruct(marker, mask);
  3 Comments
Image Analyst
Image Analyst on 1 May 2025
Will the point of interest always be at the global peak? Don't cyclones have lower pressure inside, that than higher values? So wouldn't that be the opposite of the example you gave?
Sarvesh
Sarvesh on 5 May 2025
yes, that is correct. But the method of finding the outermost contour would remain the same given a point. I have used part of your code to solve my problem in point 2 and the accepted code to solve part 1. thanks heaps

Sign in to comment.


Thorsten
Thorsten on 2 May 2025
Edited: Thorsten on 2 May 2025
The idea is to first convert the counter point data into a different data structure where the points of each contour can be accessed as a cell array, and then use cellfun to find the desired countours and their properties.
Instead of checking all closed curves within a search_radius, check whether the centerPoint is within the polygon defined by the closed curve.
% sample data
Z = peaks+10;
[M, h] = contour(Z);
% convert to array of struct CP with fields 'level' and 'points'
Mt = M; i = 1;
while ~isempty(Mt) > 0
CP(i).level = Mt(1,1);
N = Mt(2,1);
CP(i).points = Mt(:, 2:N+1);
Mt(:, 1:N+1) = [];
i = i+1;
end
centerPoint = [25, 37.5];
% find index of desired curves (closed and around centerPoint)
ind_closedaroundP = cellfun(@(x) isinterior(polyshape(x'),centerPoint) & isequal(x(:,1),x(:,end)), {CP(:).points})
levels = [CP(ind_closedaroundP).level]
theCP = {CP(ind_closedaroundP).points};
% compute statistics on these contours
dmax = cellfun(@(x) max(hypot(x(1,:) - centerPoint(1), x(2,:) - centerPoint(2))), theCP)
dmin = cellfun(@(x) min(hypot(x(1,:) - centerPoint(1), x(2,:) - centerPoint(2))), theCP)
dmean = cellfun(@(x) mean(hypot(x(1,:) - centerPoint(1), x(2,:) - centerPoint(2))), theCP)

Tags

Community Treasure Hunt

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

Start Hunting!