How can I locate a vortex center?

조회 수: 32 (최근 30일)
Ashfaq Ahmed
Ashfaq Ahmed 2022년 4월 25일
댓글: Ashfaq Ahmed 2022년 5월 9일
Hi all!
Someone in the MATLAB community once asked similar kind of question but unfortuinately, it was never answered. So, I am asking it again here.
Suppose I have a vortex field. If you please download the "vortex_velocity_data.mat" file you will be able to create the vector field which I created.
Once the code is imported to the workspace, I typed the following code, to get the quiver plot:
load('vortex_velocity_data.mat')
figure(1),
q2 = quiver(x,y,ud,vd)
box on;
set(findall(gcf,'-property','LineWidth'),'LineWidth',3);
The image is attached -
How can I locate the center of this quiver/vector plot WITHOUT using the 'Data Cursor'? What lines of code will I have to write in order to locate the center? By locating the center x and y co-ordinates I can find the relative velocity components (u,v)?
Any feedback from you will be highly appreciated :)

채택된 답변

Riccardo Scorretti
Riccardo Scorretti 2022년 4월 25일
편집: Riccardo Scorretti 2022년 4월 25일
Nice problem! I propose the solution based on the rationale that the center of the vortex (assuming that there is only one vortex) is the point where the curl is maximum:
% Load the data
load('vortex_velocity_data.mat');
figure
q2 = quiver(x, y, ud, vd);
box on;
% set(findall(gcf,'-property','LineWidth'),'LineWidth',3);
hold on
% Look for the max of the curl
[curlz, cav] = curl(x, y, ud, vd);
[~, ind] = max(abs(curlz(:)));
plot(x(ind), y(ind), 'r*');
fprintf('Center position: x0 = %f , y0 = %f\n', x(ind), y(ind));
Center position: x0 = 0.076262 , y0 = 0.084869
fprintf('Speed: u = %f , v = %f\n', ud(ind), vd(ind));
Speed: u = -0.388332 , v = -1.441110
In a more general case where your vector field is given by a couple of functions Fu(x,y) and Fv(x,y) you can try to solve the problem by dicothomy, by using more or less the same argument; only the computation of the curl has to be programmed differently:
% Load the data
load('vortex_velocity_data.mat');
figure
q2 = quiver(x, y, ud, vd);
box on;
% set(findall(gcf,'-property','LineWidth'),'LineWidth',3);
hold on
% Define interpolators
Fu = scatteredInterpolant(x(:), y(:), ud(:));
Fv = scatteredInterpolant(x(:), y(:), vd(:));
% Try a simple 2D bisection algorithm
minx = min(x(:)) ; maxx = max(x(:));
miny = min(y(:)) ; maxy = max(y(:));
maxit = 20;
for it = 1 : maxit
mx = (minx + maxx)/2;
my = (miny + maxy)/2;
set_of_circulations = [ ...
circulation(Fu,Fv,minx,mx,miny,my) ...
circulation(Fu,Fv,minx,mx,my,maxy) ...
circulation(Fu,Fv,mx,maxx,miny,my) ...
circulation(Fu,Fv,mx,maxx,my,maxy) ...
];
[~, ind] = max(abs(set_of_circulations));
if ind == 1 || ind == 2
maxx = mx;
else
minx = mx;
end
if ind == 1 || ind == 3
maxy = my;
else
miny = my;
end
end
plot(mx, my, 'r*');
return
function val = circulation(Fx, Fy, x1, x2, y1, y2)
% Compute the circulation of the vector field (Fx,Fy) along a square path
if x1 > x2 , t = x1 ; x1 = x2 ; x2 = t ; end
if y1 > y2 , t = y1 ; y1 = y2 ; y2 = t ; end
val = ...
integral(@(t) Fx(t,y1*ones(size(t))), x1, x2) + ...
integral(@(t) Fy(x2*ones(size(t)),t), y1, y2) + ...
integral(@(t) Fx(t,y2*ones(size(t))), x2, x1) + ...
integral(@(t) Fy(x1*ones(size(t)),t), y2, y1);
end
In both cases I got more or less the same figure:
  댓글 수: 6
Riccardo Scorretti
Riccardo Scorretti 2022년 4월 28일
편집: Riccardo Scorretti 2022년 4월 28일
If you have several vortex the problem is much more messy. Basically, it boils up in findining multiple peaks in a noisy graphic.
Perhaps it can be solved if you have an a priori information on the minimal distance between vortex. Hereafter a temptative program, of which I'm not much happy. It tries to find a number of vortex in your data. It has three parameters (which are tricky to determine) and gives many false positives:
clear all;
close all; clc;
% Load the data
load('vorticity_SWOT.mat');
figure
q2 = quiver(x, y, vx, vy); hold on
%normalize the field before computing the curl
nuv=sqrt(vx.^2+vy.^2);
vx=vx./nuv;
vy=vy./nuv;
% Look for the max of the curl
[curlz, cav] = curl(x, y, vx, vy);
[~, ind] = max(abs(curlz(:)));
% plot(x(ind), y(ind), 'r*','MarkerSize',18);
box on;
hold on
fprintf('Center position: x0 = %f , y0 = %f\n', x(ind), y(ind));
Center position: x0 = -2010736.133410 , y0 = 454953.143337
% These parameters are likely to be "tricky" to set up
mindist = 10000; % = minimal distance between two vortex
threshold = 0.90; % = threshold (must be in the range ]0 ; 1[ )
maxnumvortex = 10; % = max number of vortex
% Hereafter it is assume that x and y have the structure of a grid
% There are for sure smarter ways to do
[M, N] = size(x);
ii = 2:M-1;
jj = 2:N-1;
curlz = abs(curlz);
thrval = sort(curlz(:));
% figure ; plot(thrval);
thrval = thrval(round(threshold*numel(thrval)));
% figure
% pcolor(x, y, curlz); hold on ; shading interp
ind = find( ...
curlz(ii,jj) > curlz(ii+1,jj) & ...
curlz(ii,jj) > curlz(ii-1,jj) & ...
curlz(ii,jj) > curlz(ii,jj+1) & ...
curlz(ii,jj) > curlz(ii,jj-1) & ...
curlz(ii,jj) > thrval ...
);
% Mandatory tricky step
[tmp_i, tmp_j] = ind2sub([M-2 N-2], ind);
ind = sub2ind([M N], tmp_i+1, tmp_j+1);
clear tmp_i tmp_j
fprintf('Found %i candidates\n', numel(ind));
Found 825 candidates
plot(x(ind), y(ind), 'c.','MarkerSize',5);
listOf_vortex = [];
while ~isempty(ind)
[~, t] = max(curlz(ind));
i_vortex = ind(t);
listOf_vortex(end+1) = i_vortex;
% Get rid of close candidates
x0 = x(i_vortex) ; y0 = y(i_vortex);
plot(x0, y0, 'm.','MarkerSize',5);
dst2 = (x(ind) - x0).^2 + (y(ind) - y0).^2;
ind = ind(dst2 > mindist.^2);
end
fprintf('Found vortex = %i\n', numel(listOf_vortex));
Found vortex = 226
% There are still too much vortex
[~, t] = sort(curlz(listOf_vortex), 'descend');
listOf_vortex = listOf_vortex(t(1:min(maxnumvortex, numel(listOf_vortex))));
plot(x(listOf_vortex), y(listOf_vortex), 'r*','MarkerSize',18);
axis(1.0E6*[-2.0535 -1.9620 0.3917 0.5379]);
Ashfaq Ahmed
Ashfaq Ahmed 2022년 5월 9일
Hi @Riccardo Scorretti! This code is great! It's actually detecting the centers of the vortex. But I have a question: what does the cyan color points and the magenda color points indicate? Do they have something to do with the identification of curls within a distance of 10 km?

댓글을 달려면 로그인하십시오.

추가 답변 (0개)

카테고리

Help CenterFile Exchange에서 Graphics Objects에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!

Translated by