function []=findcores(frame, filtswll, swll, swul, arealim)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% %%
%% Copyright (C) 2005-2007 Anurag Singh %%
%% %%
%% This program/code snippet/file (hence forth refered as "library") %%
%% is free software; you can redistribute it and/or %%
%% modify it under the terms of the GNU Lesser General Public %%
%% License as published by the Free Software Foundation; either %%
%% version 2.1 of the License, or (at your option) any later version. %%
%% %%
%% This library is distributed in the hope that it will be useful, %%
%% but WITHOUT ANY WARRANTY; without even the implied warranty of %%
%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU %%
%% Lesser General Public License for more details. %%
%% %%
%% You should have received a copy of the GNU Lesser General Public %%
%% License along with this library; if not, write to the Free Software %%
%% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA %%
%% %%
%% FILENAME: %%
%% findcores.m %%
%% %%
%% AUTHOR: %%
%% %%
%% Anurag Singh, %%
%% MS, 2007 %%
%% Aerospace Engineering & Mechanics %%
%% University of Minnesota - Twin Cities. %%
%% Minneapolis, MN 55455 (USA) %%
%% %%
%% (currently working at: Exa Corporation, Burlington, MA 01803) %%
%% %%
%% CONTACT/EMAIL: %%
%% %%
%% anurag@aem.umn.edu %%
%% anurag9@gmail.com %%
%% %%
%% SOURCE CONTROL INFORMATION: %%
%% None (since I was planning on putting it under source control since it has %%
%% reached the critical file system size. Would be a good thing to put it under %%
%% source control while making changes. %%
%% %%
%% DESCRIPTION: %%
%% %%
%% This function identifies the vortex cores %%
%% %%
%% %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global sw2D vort2D swfrm pidfrm nidfrm swcores cores coreid
swfrm = squeeze(sw2D(frame,:,:)).* sign(squeeze(vort2D(frame,:,:)));
markcores(frame, filtswll, swll, swul, arealim);
swcores = (pidfrm > 0) + (nidfrm > 0);
cores = pidfrm + nidfrm;
coreid(frame,:,:) = cores;
return;
function []=markcores(frame, filtswll, ll, ul, arealim)
% This function identifies the vortex cores
global NX NY sw2D vort2D swfrm swfrm1 idfrm markfrm pointcnt pidfrm nidfrm numcores
idfrm = swfrm * 0;
markfrm = swfrm * 0;
pointcnt= pointcnt * 0;
% Remove all the lower swirl regions which lie outside the selected limits (ll & ul)
swfrm = swfrm.* (abs(swfrm(:,:)) > filtswll) .* (abs(swfrm(:,:)) <= ul); % Apply the threshold limit
swfrm1 = swfrm;
idtag = 0;
%disp(['Identifying +ve cores in frame # ' num2str(frame) ' ...']);
swfrm = swfrm1.* (swfrm1(:,:) > 0); % mark all the +ve cores
idtag = tagcores(idtag, arealim, ll);
pidfrm = idfrm;
idfrm = 0 * idfrm;
%disp(['Identifying -ve cores in frame # ' num2str(frame) ' ...']);
swfrm = abs(swfrm1.* (swfrm1(:,:) < 0)); % mark all the -ve cores
idtag = tagcores(idtag, arealim, ll);
nidfrm = idfrm;
idfrm = 0 * idfrm;
numcores = idtag;
%disp(['done, Total cores identified = ' num2str(idtag)]);
return;
function [idtag]=tagcores(idtag, arealim, ll)
global NX NY swfrm idfrm markfrm pointcnt
while(any(any(markfrm == 0)) && any(any(swfrm~=0))) % do this until entire sw array becomes zero or all the regions have been identified
%hile(max(max(swfrm ~= 0))) % do this until entire sw array becomes zero or all the regions have been identified
idtag = idtag + 1;
markfrm = 0 * markfrm; % zero out all previous markings
maxsw = max(max(swfrm));
%disp(['Identifying region ... ' num2str(idtag) ', Absolute Peak Vorticity = ' num2str(maxsw)]);
[seedi, seedj] = find(swfrm==maxsw);
growregion(seedi,seedj,idtag, seedi,seedj);
% This will mark all the swirl region around the maxsw by idtag in array idfrm
% the boundary of the region will be marked by -idtag
totalsw = sum(sum((idfrm==idtag) .* swfrm));
swfrm = swfrm.*(~(idfrm(:,:)==-idtag)).* (~(idfrm(:,:)==idtag));
[tx ty] = find(idfrm == idtag);
pointcnt(idtag) = size(tx,1); % no. of points identified with tag idtag
avgsw = totalsw/pointcnt(idtag);
% Remove cores with less than arealimit points
if((pointcnt(idtag) < arealim) || (maxsw < ll))
idfrm = idfrm .* (~(idfrm ==idtag)) .* (~(idfrm ==-idtag)); % delete all marks corresponding to idtag
idtag = idtag - 1;
%disp([' ... Deleted identified core at (',num2str(seedi),',',num2str(seedj),') due to area limit']);
else
%disp([' ... Identified core with peak at (',num2str(seedi),',',num2str(seedj),')']);
end
end
return;
function []=growregion(seedi, seedj, idtag, parentx, parenty)
global NX NY idfrm swfrm markfrm pointcnt
if((seedi < 1 || seedi > NX || seedj < 1 || seedj > NY)) % check if the index limits are within bounds
btag = 1;
return
end
if(markfrm(seedi,seedj) ~= 0) % check if the cell has already been visited, return if true
return;
end
if((swfrm(seedi,seedj) == 0) ) % Check if this cell has a zero value meaning its the boundary of the current region
markfrm(seedi,seedj) = 1; % mark the cell visited
return;
end
% check if the swirl values are decreasing, if not the cell is part of another core; return
if( swfrm(seedi,seedj) > swfrm(parentx,parenty) )
return;
end
if((swfrm(seedi,seedj) ~= 0) ) % Check if this cell has a non-zero swirl value
idfrm(seedi,seedj) = idtag; % if true apply mark and visit the neighboring cells
markfrm(seedi,seedj) = 1; % mark the cell visited
pointcnt(idtag) = pointcnt(idtag) + 1;
parentx = seedi;
parenty = seedj;
growregion(seedi-1, seedj , idtag, parentx, parenty);
growregion(seedi+1, seedj , idtag, parentx, parenty);
growregion(seedi , seedj-1, idtag, parentx, parenty);
growregion(seedi , seedj+1, idtag, parentx, parenty);
growregion(seedi-1, seedj-1, idtag, parentx, parenty);
growregion(seedi+1, seedj-1, idtag, parentx, parenty);
growregion(seedi-1, seedj+1, idtag, parentx, parenty);
growregion(seedi+1, seedj+1, idtag, parentx, parenty);
return;
end
return;
