alphaShape doubt to extract the number of polygons and its vertices

clc;
clear all;
close all;
t_span = 15;
dt = 0.01;
t = 0:dt:t_span;
options = optimoptions('quadprog','Display','off');
x(1) = 7.5;
y(1) = 0;
alpha_1 = 3.5;
r = 4;
grid_size = 20;
[X_grid, Y_grid] = meshgrid(linspace(-10, 10, grid_size), linspace(-10, 10, grid_size));
Wx_est = zeros(size(X_grid));
Wy_est = zeros(size(Y_grid));
Q = 0.02 * eye(2); 
R = 0.05 * eye(2); 
L = 2;
alpha = 1e-3; 
kappa = 0; 
beta = 2; 
lambda = alpha^2 * (L + kappa) - L;
gamma = sqrt(L + lambda);
for i = 1:grid_size
    for j = 1:grid_size
        state{i, j} = [0; 0];  
        P{i, j} = eye(2);     
    end
end
true_Wx = sin(2 * pi * X_grid / max(X_grid(:))) + cos(2 * pi * Y_grid / max(Y_grid(:)));
true_Wy = cos(2 * pi * X_grid / max(X_grid(:))) + sin(2 * pi * Y_grid / max(Y_grid(:)));
polygon_coords = {};
for n = 1:length(t)
    disp(['Current time = ', num2str(n*dt)]);
    xd = r * cos(t(n));
    yd = r * sin(t(n));
    xd_dot = -r * sin(t(n));
    yd_dot = r * cos(t(n));
    V = 0.5 * ((x(n) - xd)^2 + (y(n) - yd)^2);
    V_dot = [x(n) - xd, y(n) - yd];
    extra_terms = (x(n) - xd) * xd_dot + (y(n) - yd) * yd_dot;
    H = eye(2);
    f = zeros(2, 1);
    A = V_dot;
    b = -alpha_1 * V + extra_terms;
    u(:, n) = quadprog(H, f, A, b, [], [], [], [], [], options);
    for i = 1:grid_size
        for j = 1:grid_size
            x_k = state{i, j};
            P_k = P{i, j};
            sigma_points = [x_k, x_k + gamma * chol(P_k)', x_k - gamma * chol(P_k)'];
            sigma_points_pred = zeros(size(sigma_points));
            for k = 1:size(sigma_points, 2)
                sigma_points_pred(:, k) = [sigma_points(1, k) + dt * (0.1 * sin(sigma_points(2, k)));
                    sigma_points(2, k) + dt * (0.1 * cos(sigma_points(1, k)))];
            end
            x_pred = sum(sigma_points_pred, 2) / (2 * L + 1);
            P_pred = Q;
            for k = 1:size(sigma_points, 2)
                P_pred = P_pred + (sigma_points_pred(:, k) - x_pred) * (sigma_points_pred(:, k) - x_pred)';
            end
            z_pred = x_pred; 
            P_z = P_pred + R;
            K = P_pred / P_z;
            measurement = [true_Wx(i, j); true_Wy(i, j)] + mvnrnd([0; 0], R)';
            state{i, j} = x_pred + K * (measurement - z_pred);
            P{i, j} = (eye(2) - K) * P_pred * (eye(2) - K)' + K * R * K';
            Wx_est(i, j) = state{i, j}(1);
            Wy_est(i, j) = state{i, j}(2);
            if n > 1
                Wx_est(i, j) = 0.9 * Wx_est(i, j) + 0.1 * state{i, j}(1);
                Wy_est(i, j) = 0.9 * Wy_est(i, j) + 0.1 * state{i, j}(2);
            end
        end
    end
    W_magnitude = sqrt(Wx_est.^2 + Wy_est.^2);
    wind_threshold = 0.65 * max(W_magnitude(:));  
    high_wind_mask = W_magnitude > wind_threshold;
    high_wind_x = X_grid(high_wind_mask);
    high_wind_y = Y_grid(high_wind_mask);
    figure(300); clf; hold on;
    imagesc(X_grid(1, :), Y_grid(:, 1), high_wind_mask);
    colormap([1 1 1; 1 0 0]);
    colorbar;
    quiver(X_grid, Y_grid, Wx_est, Wy_est, 0.5, 'k', 'LineWidth', 1);
    
    if numel(high_wind_x) >= 4
        shp = alphaShape(high_wind_x, high_wind_y, 1);
        plot(shp, 'FaceColor', 'm', 'FaceAlpha', 0.3);
        facets = boundaryFacets(shp);
        poly_x = high_wind_x(facets(:));
        poly_y = high_wind_y(facets(:));
        polygon_coords{end+1} = [poly_x, poly_y];
    end
    
    plot(r * cos(t), r * sin(t), 'c--', 'LineWidth', 1.5);
    plot(x(1:n), y(1:n), 'b', 'LineWidth', 1.5);
    scatter(x(n), y(n), 50, 'b', 'filled');
    xlabel('X-axis');
    ylabel('Y-axis');
    title(['Wind Field & High-Wind Regions at t = ', num2str(n*dt, '%.2f')]);
    axis equal;
    grid on;
    drawnow;
    for i = 1:grid_size
        for j = 1:grid_size
            true_Wx(i, j) = true_Wx(i, j) + dt * sin(true_Wy(i, j)) + 0.05 * randn;
            true_Wy(i, j) = true_Wy(i, j) + dt * cos(true_Wx(i, j)) + 0.05 * randn;
        end
    end
    x_clamped = min(max(x(n), min(X_grid(:))), max(X_grid(:)));
    y_clamped = min(max(y(n), min(Y_grid(:))), max(Y_grid(:)));
    wind_x = interp2(X_grid, Y_grid, Wx_est, x_clamped, y_clamped, 'linear', 0);
    wind_y = interp2(X_grid, Y_grid, Wy_est, y_clamped, y_clamped, 'linear', 0);
    x(n+1) = x(n) + dt * (u(1, n) + wind_x);
    y(n+1) = y(n) + dt * (u(2, n) + wind_y);
end

For the above matlab code, alphaShape is giving me the polygons when "if numel(high_wind_x) >= 4" condition is being satisfied. I want to know that for each polygon being formed, what are the vertices being used to form these polygons?

Answer 1

[Walter Roberson]

After the call

        shp = alphaShape(high_wind_x, high_wind_y, 1);

you can do

        [tri, P] = alphaTriangulation(shp);

tri will be the triangulation, and P will be the matrix of vertices associated with the triangulation. So to get the vertices for each triangle, you would

        P(tri)

[Pallov Anand]

Thanks a lot, it works