alphaShape doubt to extract the number of polygons and its vertices

3 次查看（过去 30 天）
Pallov Anand
Pallov Anand 2025-3-21
评论： Pallov Anand 2025-3-22
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?

请先登录，再回答此问题。

采纳的回答

Walter Roberson
Walter Roberson 2025-3-21
编辑：Walter Roberson 2025-3-21
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)

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Bounding Regions 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by