| findslowregs(frame,reg_xl, reg_yl, arealim, uthresh, connectregs)
|
function findslowregs(frame,reg_xl, reg_yl, arealim, uthresh, connectregs)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% %%
%% 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: %%
%% findslowregs.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: %%
%% %%
%% Need to add %%
%% %%
%% %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global numregions newregid newregion slowregion swfrm nuidfrm tagidstack mergregid newregion
markregions(frame,arealim,uthresh);
disp(['Number of regions before connecting = ',num2str(numregions)]);
if(connectregs == 1)
swfrm = nuidfrm; % store the original found regions for comparison later
% Initialize merge id table and new idtag stack
tagidstack = (numregions+1:numregions+numregions); % Create new id tag stack
mergregid = zeros(numregions);
newregion = newregion * 0;
%Call the connection routine to join the close regions into one region
connect_regions(reg_xl, reg_yl);
temp = (newregion - numregions);
temp = temp .* (temp > 0);
slowregion = temp;
newregid = numregions - size(tagidstack,2) ; % the difference between the two would be the new number or regions after joining
disp(['New number of regions after connecting = ',num2str(newregid)]);
else
slowregion = nuidfrm;
end
return;
function []=markregions(frame,arealim, uthresh)
% This function identifies the vortex cores
global NX NY idfrm markfrm pointcnt u_norm piv_ustd
global ufrm ufrm1 nuidfrm numregions
idfrm = ufrm * 0;
markfrm = ufrm * 0;
pointcnt = pointcnt * 0;
numregions = 0;
ufrm = squeeze(u_norm(frame,:,:));
ufrm = ufrm .* (ufrm < -uthresh*piv_ustd); % extract only slow speed regions
ufrm1 = ufrm;
idtag = 0;
%disp(['Identifying slow speed regions in frame # ' num2str(frame) ' ...']);
idtag = tagregions(idtag, arealim);
nuidfrm = idfrm;
idfrm = 0 * idfrm;
numregions = idtag;
%disp(['done, Total regions identified = ' num2str(idtag)]);
return;
function [idtag]=tagregions(idtag, arealim)
global NX NY idfrm markfrm pointcnt
global ufrm
while(any(any(markfrm == 0)) && any(any(ufrm~=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
minu = min(min(ufrm));
[seedi, seedj] = find(ufrm==minu);
growregion(seedi(1), seedj(1),idtag); % Pass on only the first element in case there are two peaks with identical values
% This will mark all the region around the minu by idtag in array idfrm
% the boundary of the region will be marked by -idtag
deleteregion(idtag); % this will zero-out all the regions marked by idtag in swfrm sothat other regions could be identified.
%dbstop in findslowregs.m at 86
[tx ty] = find(idfrm == idtag);
pointcnt(idtag) = size(tx,1); % no. of points identified with tag idtag
% Remove regions with less than arealimit points
if(pointcnt(idtag) <= arealim)
idfrm = idfrm .* (~(idfrm ==idtag)) .* (~(idfrm ==-idtag)); % delete all marks corresponding to idtag
idtag = idtag - 1;
%disp([' ... Deleted identified region at (',num2str(seedi(1)),',',num2str(seedj(1)),') due to area limit']);
else
%disp([' ... Identified region with peak at (',num2str(seedi),',',num2str(seedj),')']);
end
end
return;
function []=growregion(seedi, seedj, idtag)
global NX NY idfrm ufrm markfrm pointcnt
% idtag
% dbstop if error
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
%disp([' Visiting Cell (' num2str(seedi) ',' num2str(seedj) ')']);
if((ufrm(seedi,seedj) == 0) ) % Check if this cell has a zero value meaning its the boundary of the current region
% idfrm(seedi,seedj) = -idtag; % if true apply mark and visit the neighboring cells
markfrm(seedi,seedj) = 1; % mark the cell visited
return;
end
% at this point this cell has not been visited neither has it been marked check the u values
%if((swfrm(seedi,seedj) >= ll && swfrm(seedi,seedj) <= ul) ) % Check if this cell is valid pointhas a non-zero swirl value
if((ufrm(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;
% C is the central cell and the surrounding cells are marked as following:
%
% 4 3 2
% 5 C 1
% 6 7 8
%
% C = (m,n)
growregion(seedi-1, seedj , idtag);
growregion(seedi+1, seedj , idtag);
growregion(seedi , seedj-1, idtag);
growregion(seedi , seedj+1, idtag);
growregion(seedi-1, seedj-1, idtag);
growregion(seedi+1, seedj-1, idtag);
growregion(seedi-1, seedj+1, idtag);
growregion(seedi+1, seedj+1, idtag);
end
return;
% This routine is called when all the slow speed regions is identified and now regions very close are joined to form one region
% The two regions are joined if they are separated in the streamwise and spanwise direction by a certain number of cells
% reg_xl, is the limit in the streamwise direction
% reg_yl is the limit in the spanwise direction
% the two regions are joined diagonally by a single line of cells
% Also the regions are joined only if there are no identified cores between them
% the newly joined regions acquire a new idtag which is higher than max number of regions identified, this number is increased serially
% In case a particular region is not joined to any and thus retains it old idtag would be given a new
% one in the new sequece.
% finally after all the regions have been visited idtags shifted back to start from 1
% After above step we have newly identified regions (after joining if any required) serially
function connect_regions(reg_xl, reg_yl) % reg_xl & reg_yl define the filter size in establishing connection
global nuidfrm numregions mergregid
region = nuidfrm;
for regid=1:numregions % do the connection for all the regions
[X Y] = find(region==regid); % get the indices of all regions which are tagged with regid i.e. one region at a time
for m=1:size(X) % scan the entire region for proximity to another region and make joints if needed
find_connected_regions(X(m), Y(m), reg_xl, reg_yl, region, regid);
end
% now joints has been established remove holes if any
%
% Needs to be coded
%
end
updateregions(region); % mark the original regions also by their new idtag
return;
% searches the boundary of one region and concentrates on one only at a time
% routine connect_region calls this repeatedly for covering all the regions
function find_connected_regions(seedi, seedj, reg_xl, reg_yl, region, regid)
if(bndry_check(seedi, seedj, region, regid) == 0) % check if its a boundary element or not, if not return
return;
end
% At this point we are within array bounds and its a boundary element, check for neaibouring regions
% and find if they should be joined with the current one
% get the filter grid containing neibouring regions
filtregid = getfiltergrid(seedi, seedj, reg_xl, reg_yl, region, regid);
% mark down all the region ids which fall in the filter for merging later
% Only after this step we have the regions which actually needs to be joined
% sAnd finally join the two regions with the newregid in the newregion array
joinregion(filtregid, seedi, seedj, region, regid);
return;
function check=bndry_check(x, y, region, regid)
global NX NY
% establish the boundaries of the cell (x,y)
xmin = x - 1; if(xmin < 1 ); xmin = 1 ; end;
xmax = x + 1; if(xmax > NX); xmax = NX; end;
ymin = y - 1; if(ymin < 1 ); ymin = 1 ; end;
ymax = y + 1; if(ymax > NY); ymax = NY; end;
% region(x,y) is a boundary element if and only if it is surrounded on all the sides
% i.e. all the cells must be marked by regid in the above defined grid
[X Y]=find((region(xmin:xmax,ymin:ymax) ~= regid));
if(isempty(X))
check = 0;
else
check = 1;
end
return;
function [filtgrid] = getfiltergrid(seedi, seedj, reg_xl, reg_yl, region, regid)
global NX NY
global fltx flty mergregid
xmin = seedi - reg_xl; if(xmin < 1 ); xmin = 1 ; end
ymin = seedj - reg_yl; if(ymin < 1 ); ymin = 1 ; end
xmax = seedi + reg_xl; if(xmax > NX); xmax = NX; end
ymax = seedj + reg_yl; if(ymax > NY); ymax = NY; end
filtgrid = region(xmin:xmax,ymin:ymax); % copy all the elements in this grid
[x y] = find(mergregid(:,regid)); % Find if there are any joints associated with regid
% get all the elements which are not part of the region or a joint between two region
if(isempty(x)) %The region has not been connected previously before to any other region -> no joints
filtgrid = filtgrid .* (~(filtgrid == regid));
% disp(['This region never connected, regid = ',num2str(regid),' cell (x,y) = (',num2str(seedi),',',num2str(seedj),')']);
else % this region has been connected previously and all it's joints must be considered
filtgrid = filtgrid .* (~(filtgrid == regid)) .* (~(filtgrid == mergregid(x(1),y(1))));
end
% find out the indices of the non zero elements in the filtered region
temp = region * 0;
temp(xmin:xmax, ymin:ymax) = filtgrid;
[fltx flty] = find(temp ~= 0); % this gives the co-ordinates of the cells of regions to be joined in (fltx, flty)
return
function joinregion(filtregid, x, y, region, regid)
global fltx flty
if(isempty(fltx)) % return if there is nothing to be joined
%disp('Nothing to join');
return;
end
% At this point fltx is non-zero implying there are other regions in the vicinity of the current region
% which may or may not be connected previously, the two must be connected
X = fltx;
Y = flty;
% get all the tagged regions this is inefficient way to do this, its okay if there aint many regions to be joined
for m=1:size(X) % join each non-zero element to the original seed
tempregid = update_merge_table(regid, region(X(m),Y(m)), region); % update the mergetable and get the new idtag
%['in join, tempregid =', num2str(tempregid),', V(m) = ',num2str(V(m)), ...
%', merg(',num2str(m),',',num2str(regid),') = ', num2str(mergregid(m,regid))]
% join all the cells to cell of the original region
xt = X(m);
yt = Y(m);
%['in join, ' num2str(tempregid)]
mergepoints(xt, yt, x, y, region, regid, tempregid); % establish the joint between the two points
end
return
% This function updates the mergeid table appropriately and gives the new tempregid for creating joints between regions regid1 & regid2
function tempregid = update_merge_table(regid1, regid2, region)
global mergregid newregion
%%
% now copy columns and rows corresponding to regid1 & regid2
temp = mergregid * 0 ;
temp(regid1,:) = mergregid(regid1,:); temp(:,regid1) = mergregid(:,regid1);
[x1 y1] = find(temp); % This gives co-ordinates in merge idtable of to which regid1 is connected
temp = temp*0;
temp(regid2,:) = mergregid(regid2,:); temp(:,regid2) = mergregid(:,regid2);
[x2 y2] = find(temp); % This gives co-ordinates in merge idtable of to which regid2 is connected
%disp(['idtags (regid1 regid2) = ',num2str(regid1),', ',num2str(regid2)]);
%[x1 y1 x2 y2]
%% 4 possibilities
%% reg1 reg2
%% 1- N N
%% 2- N C
%% 3- C N
%% 4- C C
% 1 would require a new tag and updating the table
% 2 & 3 would entail using the tag id from the one which is connected and updating the entire table
% 4 would require analysis as whether reg1 is should be picked or reg2. This is the only case where newregion array needs to be updated also
if(isempty(x1) && isempty(x2)) %Case 1, none connected before
% Get a new idtag, and update the mergeid table (only 2 cells need be updated
tempregid=idstack('pop',0); % Get a new idtag from the new idtag pool/stack, 2nd argument is dummy
mergregid(regid1,regid2) = tempregid;
mergregid(regid2,regid1) = tempregid;
elseif(isempty(x1) && ~(isempty(x2))) %Case 2, regid1 not connected, regid2 connected
% Use regid2's tag and update the table
tempregid = mergregid(x2(1),y2(1));
mergregid(regid1,regid2) = tempregid;
mergregid(regid2,regid1) = tempregid;
% update the mergregtable further for the other regions which are connected to region2
for m=1:size(x2)
mergregid(regid1,x2(m)) = tempregid;
mergregid(x2(m),regid1) = tempregid;
end
elseif(~(isempty(x1)) && (isempty(x2))) %Case 3, regid1 connected, regid2 not connected
% Use regid1's tag and update the table
tempregid = mergregid(x1(1),y1(1));
mergregid(regid1,regid2) = tempregid;
mergregid(regid2,regid1) = tempregid;
% update the mergregtable further for the other regions which are connected to region1
for m=1:size(x1)
mergregid(regid2,x1(m)) = tempregid;
mergregid(x1(m),regid2) = tempregid;
end
elseif(mergregid(x1(1),y1(1)) == mergregid(x2(1),y2(1))) % no table update required in this case
tempregid = mergregid(x2(1),y2(1));
elseif(~(isempty(x1)) && ~(isempty(x2))) %Case 4, regid1 connected, regid2 connected
%Find which region has a smaller tag, and use that to update the tables
% Also update the newregion array
% Also update the newregid counter, since two sets of already connected regions are futher connected.
idv1 = mergregid(x1(1), y1(1)); % get the value of the idtags for regid1
idv2 = mergregid(x2(1), y2(1)); % get the value of the idtags for regid1
% [(isempty(x1)) (isempty(x2))]
% [idv1 idv2]
if(idv1 < idv2) % Use region1 tag
tempregid = mergregid(x1(1),y1(1));
% update all the connected regions in region2 to region1 tag
mergregid(regid1,regid2) = tempregid;
mergregid(regid2,regid1) = tempregid;
% update the mergregtable further for the other regions which are connected to region1s
for m=1:size(x1)
mergregid(regid2,x1(m)) = tempregid;
mergregid(x1(m),regid2) = tempregid;
end
%update newregion array, i.e, mark all the joints corresponding to region2 to the new tag
temp = idv1 * (newregion == idv2);
newregion = newregion .* (~(temp)) + temp ;
% Update the new idtag stack, i.e send the idv2 to the pool of available tags
idstack('push',idv2);
return;
elseif(idv1 > idv2) % Use region2 tag
tempregid = mergregid(x2(1),y2(1));
% update all the connected regions in region1 to region2 tag
mergregid(regid1,regid2) = tempregid;
mergregid(regid2,regid1) = tempregid;
% update the mergregtable further for the other regions which are connected to region2
for m=1:size(x2)
mergregid(regid1,x2(m)) = tempregid;
mergregid(x2(m),regid1) = tempregid;
end
temp = idv2 * (newregion == idv1);
newregion = newregion .* (~(temp)) + temp ;
% Update the new idtag stack, i.e send the idv2 to the pool of available tags
idstack('push',idv1);
end
else
disp('ERROR in update merge table');
end
return;
% This functions simulates a pool of available id tags, and stores them in the poll or returns them to the calling program
% tagidstack must be initiialized outside this routine before making use of it
function newregid=idstack(operation,newregid)
global tagidstack
switch operation
case 'pop'
newregid = tagidstack(1);
tagidstack(1) = [];
%disp(['Generated new id tag = ',num2str(newregid)])
return;
case 'push'
tagidstack = [newregid tagidstack];
%disp(['Added to the stack, id tag = ',num2str(newregid)])
return;
otherwise
disp('ERROR: Wrong Stack Operation for newregid stack')
end
return
function []=updateregions(region)
%zero out the region in frame swfrm marked by idtag in idfrm
global newregid newregion mergregid
for m=1:size(mergregid,1)
temp = mergregid * 0 ;
temp(m,:) = mergregid(m,:);
temp(:,m) = mergregid(:,m);
[t1 t2] = find(temp);
%[t1 t2] = find(mergregid(:,m)); % get the new idtag corresponding to old region 'm'
% check if t1 is empty, if true, it means this region is not connected to anyone and
% old tags were retained. These should be renewed
if(isempty(t1))
tempregid = idstack('pop',0); % Get a new idtag from the new idtag pool/stack
else % This region has been connected previously, pick its id tag from the merge table and mark the region
tempregid = mergregid(t1(1),t2(1));
end
%['in update, tempregid =', num2str(tempregid)]%, ', newregion(59,23)=', num2str(newregion(59,23))]%,', t3(1) = ',num2str(t3(1)), ...
%', merg(',num2str(t1(1)),',',num2str(t2(1)),') = ', num2str(mergregid(t1(1),t2(1)))]
newregion = newregion + tempregid * (region == m); % Mark this region with the new tag
end
return;
function []=deleteregion(idtag)
%zero out the region in frame swfrm marked by idtag in idfrm
global idfrm ufrm
ufrm = ufrm.*(~(idfrm(:,:)==-idtag)).* (~(idfrm(:,:)==idtag));
return;
function mergepoints(xt,yt, x, y, region, regid, tempregid)
global newregion ufrm1
%disp('in mergepoint');
%[xt yt x y]
if((xt == y) && (yt == y))
return;
end
if(region(xt,yt) == regid)
return;
end
% mark the region (do so only if its not part of any other parent regions
if(region(xt,yt) == 0)
newregion(xt,yt) = tempregid;
newregion(xt,yt);
end
% call the routine again depending with updated xt, yt
dx = sign(x-xt); % increment 1 or -1 depending on where x is
dy = sign(y-yt); % increment 1 or -1 depending on where y is
mergepoints(xt+dx, yt+dy, x, y, region, regid, tempregid);
return;
|
|
