source: trunk/src/uvmat_doc/FUNCTIONS_DOC/RUN_FIX.html @ 37

Last change on this file since 37 was 37, checked in by sommeria, 14 years ago

create_grid.fig ,
uvmat_doc and all the included files added
File size: 16.7 KB
Line 
1<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN"
2                "http://www.w3.org/TR/REC-html40/loose.dtd">
3<html>
4<head>
5  <title>Description of RUN_FIX</title>
6  <meta name="keywords" content="RUN_FIX">
7  <meta name="description" content="'RUN_FIX': function for fixing velocity fields:">
8  <meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
9  <meta name="generator" content="m2html &copy; 2003 Guillaume Flandin">
10  <meta name="robots" content="index, follow">
11  <link type="text/css" rel="stylesheet" href="../m2html.css">
12</head>
13<body>
14<a name="_top"></a>
15<div><a href="../index.html">Home</a> &gt;  <a href="index.html">.</a> &gt; RUN_FIX.m</div>
16
17<!--<table width="100%"><tr><td align="left"><a href="../index.html"><img alt="<" border="0" src="../left.png">&nbsp;Master index</a></td>
18<td align="right"><a href="index.html">Index for .&nbsp;<img alt=">" border="0" src="../right.png"></a></td></tr></table>-->
19
20<h1>RUN_FIX
21</h1>
22
23<h2><a name="_name"></a>PURPOSE <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
24<div class="box"><strong>'RUN_FIX': function for fixing velocity fields:</strong></div>
25
26<h2><a name="_synopsis"></a>SYNOPSIS <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
27<div class="box"><strong>function error=RUN_FIX(filename,field,flagindex,iter,thresh_vecC,flag_mask,maskname,thresh_vel,inf_sup,fileref,fieldref) </strong></div>
28
29<h2><a name="_description"></a>DESCRIPTION <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
30<div class="fragment"><pre class="comment">'RUN_FIX': function for fixing velocity fields:
31-----------------------------------------------
32 RUN_FIX(filename,field,flagindex,thresh_vecC,thresh_vel,iter,flag_mask,maskname,fileref,fieldref)
33
34filename: name of the netcdf file (used as input and output)
35field: structure specifying the field to fix (field.vel_type='civ1' or 'civ2')
36flagindex: flag specifying which values of vec_f are removed:
37 if flagindex(1)=1: vec_f=-2 vectors are removed
38 if flagindex(2)=1: vec_f=3 vectors are removed
39 if flagindex(3)=1: vec_f=2 vectors are removed (if iter=1) or vec_f=4 vectors are removed (if iter=2)
40iter=1 for civ1 fields and iter=2 for civ2 fields
41thresh_vecC: threshold in the image correlation vec_C
42flag_mask: =1 mask used to remove vectors (0 else)
43maskname: name of the mask image file for fix
44thresh_vel: threshold on velocity, or on the difference with the reference file fileref if exists
45inf_sup=1: remove values smaller than threshold thresh_vel, =2, larger than threshold
46fileref: .nc file name for a reference velocity (='': refrence 0 used)
47fieldref: 'civ1','filter1'...feld used in fileref</pre></div>
48
49<!-- crossreference -->
50<h2><a name="_cross"></a>CROSS-REFERENCE INFORMATION <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
51This function calls:
52<ul style="list-style-image:url(../matlabicon.gif)">
53<li><a href="griddata_uvmat.html" class="code" title="function ZI = griddata_uvmat(X,Y,Z,XI,YI)">griddata_uvmat</a>     'griddata_uvmat': function griddata_uvmat(vec2_X,vec2_Y,vec2_U,vec_X,vec_Y,'linear')</li><li><a href="read_civxdata.html" class="code" title="function [Field,VelTypeOut]=read_civxdata(filename,FieldNames,VelType)">read_civxdata</a> 'read_civxdata': reads civx data from netcdf files</li><li><a href="warndlg_uvmat.html" class="code" title="function hwarn=warndlg_uvmat(warntext,title)">warndlg_uvmat</a>     'warndlg_uvmat': display warning message (error, warning, confirmation) in a given figure</li></ul>
54This function is called by:
55<ul style="list-style-image:url(../matlabicon.gif)">
56<li><a href="civ.html" class="code" title="function varargout = civ(varargin)">civ</a>  'civ': function associated with the interface 'civ.fig' for PIV, spline interpolation and stereo PIV (patch)</li></ul>
57<!-- crossreference -->
58
59
60<h2><a name="_source"></a>SOURCE CODE <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
61<div class="fragment"><pre>0001 <span class="comment">%'RUN_FIX': function for fixing velocity fields:</span>
620002 <span class="comment">%-----------------------------------------------</span>
630003 <span class="comment">% RUN_FIX(filename,field,flagindex,thresh_vecC,thresh_vel,iter,flag_mask,maskname,fileref,fieldref)</span>
640004 <span class="comment">%</span>
650005 <span class="comment">%filename: name of the netcdf file (used as input and output)</span>
660006 <span class="comment">%field: structure specifying the field to fix (field.vel_type='civ1' or 'civ2')</span>
670007 <span class="comment">%flagindex: flag specifying which values of vec_f are removed:</span>
680008         <span class="comment">% if flagindex(1)=1: vec_f=-2 vectors are removed</span>
690009         <span class="comment">% if flagindex(2)=1: vec_f=3 vectors are removed</span>
700010         <span class="comment">% if flagindex(3)=1: vec_f=2 vectors are removed (if iter=1) or vec_f=4 vectors are removed (if iter=2)</span>
710011 <span class="comment">%iter=1 for civ1 fields and iter=2 for civ2 fields</span>
720012 <span class="comment">%thresh_vecC: threshold in the image correlation vec_C</span>
730013 <span class="comment">%flag_mask: =1 mask used to remove vectors (0 else)</span>
740014 <span class="comment">%maskname: name of the mask image file for fix</span>
750015 <span class="comment">%thresh_vel: threshold on velocity, or on the difference with the reference file fileref if exists</span>
760016 <span class="comment">%inf_sup=1: remove values smaller than threshold thresh_vel, =2, larger than threshold</span>
770017 <span class="comment">%fileref: .nc file name for a reference velocity (='': refrence 0 used)</span>
780018 <span class="comment">%fieldref: 'civ1','filter1'...feld used in fileref</span>
790019
800020 <a name="_sub0" href="#_subfunctions" class="code">function error=RUN_FIX(filename,field,flagindex,iter,thresh_vecC,flag_mask,maskname,thresh_vel,inf_sup,fileref,fieldref)</a>
810021 error=[]; <span class="comment">%default</span>
820022 vel_type{1}=field.vel_type;
830023 <span class="comment">%[nbc,nbd,Civ,CivStage,tt,Field]=read_ncfield(filename,field.vel_type,'vec_color','ima_cor');</span>
840024 Field=<a href="read_civxdata.html" class="code" title="function [Field,VelTypeOut]=read_civxdata(filename,FieldNames,VelType)">read_civxdata</a>(filename,[],field.vel_type);
850025 <span class="keyword">if</span> isfield(Field,<span class="string">'Txt'</span>)
860026     error=Field.Txt; <span class="comment">%error in reading</span>
870027     <span class="keyword">return</span>
880028 <span class="keyword">end</span>
890029
900030 <span class="keyword">if</span> ~isfield(Field,<span class="string">'X'</span>) || ~isfield(Field,<span class="string">'Y'</span>) || ~isfield(Field,<span class="string">'U'</span>) || ~isfield(Field,<span class="string">'V'</span>)
910031     error=[<span class="string">'input file '</span> filename <span class="string">' does not contain vectors in RUN_FIX.m'</span>]; <span class="comment">%bad input file</span>
920032     <span class="keyword">return</span>
930033 <span class="keyword">end</span>
940034 <span class="keyword">if</span> ~isfield(Field,<span class="string">'C'</span>)
950035     Field.C=ones(size(Field.X));<span class="comment">%correlation=1 by default</span>
960036 <span class="keyword">end</span>
970037 <span class="keyword">if</span> ~isfield(Field,<span class="string">'F'</span>)
980038     Field.F=ones(size(Field.X));<span class="comment">%warning flag=1 by default</span>
990039 <span class="keyword">end</span>
1000040 <span class="keyword">if</span> ~isfield(Field,<span class="string">'FF'</span>)
1010041     Field.FF=zeros(size(Field.X));<span class="comment">%fixflag=0 by default</span>
1020042 <span class="keyword">end</span>
1030043
1040044 vec_f_unit=abs(Field.F)-10*double(uint16(abs(Field.F)/10)); <span class="comment">%unityterm of vec_F in abs value</span>
1050045 vec_f_sign=sign(Field.F).*vec_f_unit;<span class="comment">% gives the unity digit of vec_f with correct sign</span>
1060046 flag1=(flagindex(1)==1)&amp;(vec_f_sign==-2);<span class="comment">%removed vectors vec_f=-2</span>
1070047 flag2=(flagindex(2)==1)&amp;(vec_f_sign==3);<span class="comment">%removed vectors vec_f=3</span>
1080048 <span class="keyword">if</span> iter==1
1090049         flag3=(flagindex(3)==1)&amp;(vec_f_sign==2); <span class="comment">% Hart vectors</span>
1100050 <span class="keyword">elseif</span> iter==2
1110051         flag3=(flagindex(3)==1)&amp;(vec_f_sign==4); <span class="comment">%</span>
1120052 <span class="keyword">end</span>
1130053             <span class="comment">%  flag4=(vec_Fmin2_2_Fmin2_2_c &lt;</span>
1140054             <span class="comment">%  colvec1)&amp;(vec_Fmin2_2_Fmin2_2_f_sign~=2)&amp;(vec_Fmin2_2_Fmin2_</span>
1150055             <span class="comment">%  2_f~=0); % =1 for red vectors not flaged</span>
1160056
1170057 flag4=(Field.C &lt; thresh_vecC)&amp;(flag1~=1)&amp;(flag2~=1)&amp;(flag3~=1); <span class="comment">% =1 for low vec_C vectors not previously removed</span>
1180058
1190059 <span class="comment">% criterium on velocity values</span>
1200060 delta_u=Field.U;<span class="comment">%default without ref file</span>
1210061 delta_v=Field.V;
1220062 <span class="keyword">if</span> exist(<span class="string">'fileref'</span>,<span class="string">'var'</span>) &amp;&amp; ~isempty(fileref)
1230063     <span class="keyword">if</span> ~exist(fileref,<span class="string">'file'</span>)
1240064         error=<span class="string">'reference file not found in RUN_FIX.m'</span>;
1250065         display(error);
1260066         <span class="keyword">return</span>
1270067     <span class="keyword">end</span>
1280068     
1290069 <span class="comment">%     [nbc,nbd,Civ,CivStage,tt,FieldRef]=read_ncfield(fileref,fieldref);</span>
1300070     FieldRef=<a href="read_civxdata.html" class="code" title="function [Field,VelTypeOut]=read_civxdata(filename,FieldNames,VelType)">read_civxdata</a>(fileref,[],fieldref);
1310071     
1320072     <span class="keyword">if</span> isfield(FieldRef,<span class="string">'FF'</span>)
1330073         index_true=find(FieldRef.FF==0);
1340074         FieldRef.X=FieldRef.X(index_true);
1350075         FieldRef.Y=FieldRef.Y(index_true);
1360076         FieldRef.U=FieldRef.U(index_true);
1370077         FieldRef.V=FieldRef.V(index_true);
1380078     <span class="keyword">end</span>
1390079     <span class="keyword">if</span> ~isfield(FieldRef,<span class="string">'X'</span>) || ~isfield(FieldRef,<span class="string">'Y'</span>) || ~isfield(FieldRef,<span class="string">'U'</span>) || ~isfield(FieldRef,<span class="string">'V'</span>)
1400080         error=<span class="string">'reference file is not a velocity field in RUN_FIX.m '</span>; <span class="comment">%bad input file</span>
1410081         <span class="keyword">return</span>
1420082     <span class="keyword">end</span>
1430083     <span class="keyword">if</span> length(FieldRef.X)&lt;=1
1440084         errordlg(<span class="string">'reference field with one vector or less in RUN_FIX.m'</span>)
1450085         <span class="keyword">return</span>
1460086     <span class="keyword">end</span>
1470087     vec_U_ref=<a href="griddata_uvmat.html" class="code" title="function ZI = griddata_uvmat(X,Y,Z,XI,YI)">griddata_uvmat</a>(FieldRef.X,FieldRef.Y,FieldRef.U,Field.X,Field.Y);  <span class="comment">%interpolate vectors in the ref field</span>
1480088     vec_V_ref=<a href="griddata_uvmat.html" class="code" title="function ZI = griddata_uvmat(X,Y,Z,XI,YI)">griddata_uvmat</a>(FieldRef.X,FieldRef.Y,FieldRef.V,Field.X,Field.Y);  <span class="comment">%interpolate vectors in the ref field to the positions  of the main field</span>
1490089     delta_u=Field.U-vec_U_ref;<span class="comment">%take the difference with the interpolated ref field</span>
1500090     delta_v=Field.V-vec_V_ref;
1510091 <span class="keyword">end</span>
1520092
1530093 thresh_vel_x=thresh_vel;
1540094 thresh_vel_y=thresh_vel;
1550095 <span class="keyword">if</span> isequal(inf_sup,1)
1560096     flag5=abs(delta_u)&lt;thresh_vel_x &amp; abs(delta_v)&lt;thresh_vel_y &amp;(flag1~=1)&amp;(flag2~=1)&amp;(flag3~=1)&amp;(flag4~=1);
1570097 <span class="keyword">elseif</span> isequal(inf_sup,2)
1580098     flag5=(abs(delta_u)&gt;thresh_vel_x | abs(delta_v)&gt;thresh_vel_y) &amp;(flag1~=1)&amp;(flag2~=1)&amp;(flag3~=1)&amp;(flag4~=1);
1590099 <span class="keyword">end</span>
1600100
1610101             <span class="comment">% flag7 introduce a grey mask, matrix M</span>
1620102 <span class="keyword">if</span> isequal (flag_mask,1)
1630103    M=imread(maskname);
1640104    nxy=size(M);
1650105    M=reshape(M,1,nxy(1)*nxy(2));
1660106    rangx0=[0.5 nxy(2)-0.5];
1670107    rangy0=[0.5 nxy(1)-0.5];
1680108    vec_x1=Field.X-Field.U/2;<span class="comment">%beginning points</span>
1690109    vec_x2=Field.X+Field.U/2;<span class="comment">%end points of vectors</span>
1700110    vec_y1=Field.Y-Field.V/2;<span class="comment">%beginning points</span>
1710111    vec_y2=Field.Y+Field.V/2;<span class="comment">%end points of vectors</span>
1720112    indx=1+round((nxy(2)-1)*(vec_x1-rangx0(1))/(rangx0(2)-rangx0(1)));<span class="comment">% image index x at abcissa vec_x</span>
1730113    indy=1+round((nxy(1)-1)*(vec_y1-rangy0(1))/(rangy0(2)-rangy0(1)));<span class="comment">% image index y at ordinate vec_y</span>
1740114    test_in=~(indx &lt; 1 |indy &lt; 1 | indx &gt; nxy(2) |indy &gt; nxy(1)); <span class="comment">%=0 out of the mask image, 1 inside</span>
1750115    indx=indx.*test_in+(1-test_in); <span class="comment">%replace indx by 1 out of the mask range</span>
1760116    indy=indy.*test_in+(1-test_in); <span class="comment">%replace indy by 1 out of the mask range</span>
1770117    ICOMB=((indx-1)*nxy(1)+(nxy(1)+1-indy));<span class="comment">%determine the indices in the image reshaped in a Matlab vector</span>
1780118    Mvalues=M(ICOMB);
1790119 <span class="comment">%    size(indx)</span>
1800120 <span class="comment">%     size(indy)</span>
1810121 <span class="comment">%    Ma=M(indy,indx);% relevant mask values</span>
1820122 <span class="comment">%    Mvalues=(diag(Ma,0))'; % mask values for vector indices, flag7=1 for masked vectors</span>
1830123    flag7b=((20 &lt; Mvalues) &amp; (Mvalues &lt; 200))| ~test_in';
1840124    indx=1+round((nxy(2)-1)*(vec_x2-rangx0(1))/(rangx0(2)-rangx0(1)));<span class="comment">% image index x at abcissa Field.X</span>
1850125    indy=1+round((nxy(1)-1)*(vec_y2-rangy0(1))/(rangy0(2)-rangy0(1)));<span class="comment">% image index y at ordinate vec_y</span>
1860126    test_in=~(indx &lt; 1 |indy &lt; 1 | indx &gt; nxy(2) |indy &gt; nxy(1)); <span class="comment">%=0 out of the mask image, 1 inside</span>
1870127    indx=indx.*test_in+(1-test_in); <span class="comment">%replace indx by 1 out of the mask range</span>
1880128    indy=indy.*test_in+(1-test_in); <span class="comment">%replace indy by 1 out of the mask range</span>
1890129    ICOMB=((indx-1)*nxy(1)+(nxy(1)+1-indy));<span class="comment">%determine the indices in the image reshaped in a Matlab vector</span>
1900130    Mvalues=M(ICOMB);
1910131 <span class="comment">%    Ma=M(indy,indx);% relevant mask values</span>
1920132 <span class="comment">%    Mvalues=(diag(Ma,0))'; % mask values for vector indices, flag7=1 for masked vectors</span>
1930133    flag7e=((Mvalues &gt; 20) &amp; (Mvalues &lt; 200))| ~test_in';
1940134    flag7=(flag7b|flag7e)';
1950135 <span class="keyword">else</span>
1960136    flag7=0;
1970137 <span class="keyword">end</span>   
1980138 <span class="comment">% flagmagenta=flag1|flag2|flag3|flag4|flag7;</span>
1990139 flagmagenta=flag1|flag2|flag3|flag4|flag5|flag7;
2000140 <span class="comment">%write fix flags</span>
2010141
2020142 [errorread,message]=fileattrib(filename);
2030143
2040144 <span class="keyword">if</span> ~isempty(message) &amp; ~isequal(message.UserWrite,1)
2050145      <a href="warndlg_uvmat.html" class="code" title="function hwarn=warndlg_uvmat(warntext,title)">warndlg_uvmat</a>([<span class="string">'no writting access to '</span> filename <span class="string">' (RUN_FIX.m)'</span>],<span class="string">'ERROR'</span>);
2060146     <span class="keyword">return</span>
2070147 <span class="keyword">end</span>
2080148 nc=netcdf(filename,<span class="string">'write'</span>); <span class="comment">%open netcdf file for writing</span>
2090149 result=redef(nc);
2100150 <span class="keyword">if</span> isempty(result), errordlg(<span class="string">'##Bad redef operation.'</span>),<span class="keyword">end</span> 
2110151 <span class="keyword">if</span> iter==1
2120152     nc.fix=1;
2130153 <span class="keyword">elseif</span> iter==2
2140154     nc.fix2=1;
2150155 <span class="keyword">end</span>
2160156 theDim=nc(field.nb) ;<span class="comment">% get the number of velocity vectors</span>
2170157 <span class="comment">% nb_vectors=size(theDim);</span>
2180158 <span class="comment">% nb_vectors=nb_vectors(1);</span>
2190159 nc{field.fixflag}=ncfloat(field.nb);
2200160 fixflag_unit=Field.FF-10*floor(Field.FF/10); <span class="comment">%unity term of fix_flag</span>
2210161 nc{field.fixflag}(:)=fixflag_unit+10*flagmagenta;
2220162 close(nc);</pre></div>
223<hr><address>Generated on Fri 13-Nov-2009 11:17:03 by <strong><a href="http://www.artefact.tk/software/matlab/m2html/">m2html</a></strong> &copy; 2003</address>
224</body>
225</html>
Note: See TracBrowser for help on using the repository browser.