closed countour around a given point

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;

Answer 1

[Star Strider]

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.

.

[Sarvesh]

thanks @Star Strider.

[Star Strider]

As always, my pleasure!

Answer 2

[Image Analyst]

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);

[Sarvesh]

Thank you for responding.

I also considered this approach initially. However, the challenge lies in the variability of the threshold I need to apply for mask. I am analyzing hrly gridded pressure data over cyclones, specifically aiming to extract the outermost closed isobar(which basically is the outermost closed contour). This provides an indication of the cyclone's size- which i am after.

Due to the fluctuating nature of gridded pressure data as the cyclone intensifies or weakens, I must dynamically adjust the mask threshold which is not very straightforward. The threshold for this adjustment depends on the surrounding pressure conditions. There are also biases in the data which prevents me from using a 'fixed' mask value. what is certain in the data are the pattern of the contours(isobars). The best approach in this case (as far as I know) is therefore, if I can identify the outermost closed contour about my point of interest using contour/contourc, I will not need to worry about adjesting the mask values or the biases.

i hope this makes some sense to the issue i have.

[Image Analyst]

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]

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

Answer 3

[Thorsten]

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)

[Sarvesh]

thanks @Thorsten