How can I locate a vortex center?

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

Answer 1

[Riccardo Scorretti]

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:

[Bruno Luong]

+1.

I notice (by eye) that the field seem to get weaker close to the vortex center.

I propose to normalize the field before computing the curl

nuv=sqrt(ud.^2+vd.^2);
ud=ud./nuv;
vd=vd./nuv;
% ...

[Ashfaq Ahmed]

Thank you so much @Riccardo Scorretti. It's been one of the most detailed and finest helps I have received from this community! Really thankful for your effort on this question. :)

[Ashfaq Ahmed]

Hi @Riccardo Scorretti, I have applied your first code with this vorticity field. It is apparently detecting the vortex eye, but not the strongest one in the field. Here is the code -

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

Could you please give me an idea on how to resolve this issue? Thank you so much for your help!

The .mat file that I am working with is attached to the comment. It's very similar to the previous one. Only change is it has several vortex eyes.

[Ashfaq Ahmed]

@Bruno Luong thank you for your idea! If I don't do the normalization, the algorithm sometimes misses out the proper detection of the eye (depending on the nature of the velocity filed). But when I normalize the field, it always detects the eye. Thank you again!

[Riccardo Scorretti]

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]

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?