필터 지우기
필터 지우기

Extract area from elipsoidal/cirle trajectory.

조회 수: 2 (최근 30일)
Claudio Iturra
Claudio Iturra 2024년 5월 28일
편집: Matt J 2024년 5월 28일
Dear all, how are you?
I created the following code to simulate an anticlockwise rotating trajectory forced by an external velocity. Now, I want to calculate the area each time the trajectory makes a closed circle, as described in the attached image. I would appreciate any help in finding the best way to identify the areas inside the circles for each loop/circle in the trajectory.
For circle detection, I used the polyshape function as described in the attached figure. Unfortunately, this requires me to cut the series each time a circular trajectory is identified, and then manually choose the area I need, which takes too much time
Thanks for your time!
clear all,close all;
Fs = 1/0.5;
duration_hours = 25 * 24;
t = 0:1/Fs:duration_hours;
period = 20;
phase = -pi/2;
num_periods = floor(duration_hours / period);
amplitudes = 0.1/2 + 0.1 * rand(1, num_periods); % Random amplitudes between 0.1 and 0.3
y = zeros(size(t));
for i = 1:num_periods
start_idx = (i-1) * period * Fs + 1;
end_idx = i * period * Fs;
uI(start_idx:end_idx) = amplitudes(i) * sin(2*pi*(t(start_idx:end_idx) - t(start_idx))/period);
end
amplitudes = 0.1/2 + 0.1 * rand(1, num_periods); % Random amplitudes between 0.1 and 0.3
for i = 1:num_periods
start_idx = (i-1) * period * Fs + 1;
end_idx = i * period * Fs;
vI(start_idx:end_idx) = amplitudes(i) * sin(2*pi*(t(start_idx:end_idx) - t(start_idx))/period);
end
uI = uI(10:end);
vI = vI(1:length(uI));
cv = complex(uI*100,vI*100);
dt=0.5/24;
dt=dt*3600*24;
cxI=cumsum(cv)*dt;
cxI=cxI/100/1000;
x = real(cxI + complex(cumsum(repmat(0.1,1191,1)),0)' );
y = imag(cxI + complex(cumsum(repmat(0.1,1191,1)),0)' );
clearvars -except x y;
plot(x,y),axis equal;
% for the first circle detected
pgonAll = polyshape(x(1:48),y(1:48),'Simplify',false);
pgonAll = simplify(pgonAll);
pgonEach = regions(pgonAll);
figure
plot(x(1:48),y(1:48))
hold on
plot(pgonEach(3))

채택된 답변

Matt J
Matt J 2024년 5월 28일
편집: Matt J 2024년 5월 28일
load data
ymin=min(y);
[xx,yy]=deal(x,y);
xx(end+1)=x(end); yy(end+1)=ymin;
p=regions(polyshape(xx,yy));
Warning: Polyshape has duplicate vertices, intersections, or other inconsistencies that may produce inaccurate or unexpected results. Input data has been modified to create a well-defined polyshape.
h=arrayfun(@(P)min(P.Vertices(:,2)), p );
p=p(h>0.001+ymin);
Areas=area(p)'
Areas = 1x29
1.8247 0.9686 0.0081 1.0051 0.5726 0.0032 0.1837 0.5910 2.8787 1.0814 0.0241 1.7894 2.5199 0.7346 1.0110 0.0168 0.3967 2.1162 1.7427 2.5049 0.7968 1.8758 1.2651 0.6767 1.3316 1.8829 0.2996 1.1132 0.0252
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
plot(x,y); hold on
plot(p,'FaceColor','r'); hold off
axis([50,100,-5,5]); daspect([1,1,1])

추가 답변 (0개)

제품


릴리스

R2023b

Community Treasure Hunt

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

Start Hunting!

Translated by