source: trunk/src/uvmat_doc/FUNCTIONS_DOC/RUN_STLIN.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: 17.3 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_STLIN</title>
6  <meta name="keywords" content="RUN_STLIN">
7  <meta name="description" content="'RUN_STLIN': combine velocity fields for stereo PIV">
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_STLIN.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_STLIN
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_STLIN': combine velocity fields for stereo PIV</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 RUN_STLIN(file_A,file_B,vel_type,file_st,nx_patch,ny_patch,thresh_patch,fileAxml,fileBxml) </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_STLIN': combine velocity fields for stereo PIV
31 file_A,file_B: input velocity files
32vel_type: string ='civ1' or 'civ2'</pre></div>
33
34<!-- crossreference -->
35<h2><a name="_cross"></a>CROSS-REFERENCE INFORMATION <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
36This function calls:
37<ul style="list-style-image:url(../matlabicon.gif)">
38<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="imadoc2struct.html" class="code" title="function [s,errormsg]=imadoc2struct(ImaDoc)">imadoc2struct</a>    'imadoc2struct': reads the xml file for image documentation</li><li><a href="phys_XYZ.html" class="code" title="function [Xphys,Yphys,Zphys]=phys_XYZ(Calib,X,Y,Z)">phys_XYZ</a>        'phys_XYZ':transforms image (px) to real world (phys) coordinates using geometric calibration parameters</li><li><a href="px_XYZ.html" class="code" title="function [X,Y]=px_XYZ(Calib,Xphys,Yphys,Zphys)">px_XYZ</a>   'px_XYZ': transform phys coordinates to image coordinates (px)</li><li><a href="pxcm_tsai.html" class="code" title="function [A11,A12,A13,A21,A22,A23]=pxcm_tsai(a,var_phys)">pxcm_tsai</a>     'pxcm_tsai': find differentials of the Tsai calibration</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="struct2nc.html" class="code" title="function errormsg=struct2nc(flname,Data)">struct2nc</a> 'struct2nc': create a netcdf file from a Matlab structure</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>
39This function is called by:
40<ul style="list-style-image:url(../matlabicon.gif)">
41<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>
42<!-- crossreference -->
43
44
45<h2><a name="_source"></a>SOURCE CODE <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
46<div class="fragment"><pre>0001 <span class="comment">%'RUN_STLIN': combine velocity fields for stereo PIV</span>
470002 <span class="comment">% file_A,file_B: input velocity files</span>
480003 <span class="comment">%vel_type: string ='civ1' or 'civ2'</span>
490004 <a name="_sub0" href="#_subfunctions" class="code">function RUN_STLIN(file_A,file_B,vel_type,file_st,nx_patch,ny_patch,thresh_patch,fileAxml,fileBxml)</a>
500005                 
510006  [XmlDataA,error]=<a href="imadoc2struct.html" class="code" title="function [s,errormsg]=imadoc2struct(ImaDoc)">imadoc2struct</a>(fileAxml);
520007  [XmlDataB,error]=<a href="imadoc2struct.html" class="code" title="function [s,errormsg]=imadoc2struct(ImaDoc)">imadoc2struct</a>(fileBxml);
530008  npxA=[]; npyA=[]; pxB=[]; npyB=[];
540009  <span class="keyword">if</span> isfield(XmlDataA,<span class="string">'Camera'</span>) &amp;&amp; isfield(XmlDataB,<span class="string">'Camera'</span>)
550010       <span class="keyword">if</span> isfield(XmlDataA.Camera,<span class="string">'ImageSize'</span>)&amp;&amp; isfield(XmlDataB.Camera,<span class="string">'ImageSize'</span>)
560011           ImageSizeA=XmlDataA.Camera.ImageSize;
570012           ImageSizeB=XmlDataB.Camera.ImageSize;
580013           <span class="keyword">if</span> ~isempty(ImageSizeA)&amp;&amp; ~isempty(ImageSizeB)
590014                xindex=findstr(ImageSizeA,<span class="string">'x'</span>);
600015                <span class="keyword">if</span> length(xindex)&gt;=2
610016                     npxA=str2num(ImageSizeA(1:xindex(1)-1));
620017                     npyA=str2num(ImageSizeA(xindex(1)+1:xindex(2)-1));
630018                <span class="keyword">end</span>
640019                xindex=findstr(ImageSizeB,<span class="string">'x'</span>);
650020                <span class="keyword">if</span> length(xindex)&gt;=2
660021                     npxB=str2num(ImageSizeB(1:xindex(1)-1));
670022                     npyB=str2num(ImageSizeB(xindex(1)+1:xindex(2)-1));
680023                <span class="keyword">end</span>
690024           <span class="keyword">end</span>
700025      <span class="keyword">end</span>
710026  <span class="keyword">end</span>
720027  <span class="keyword">if</span> isempty(npxA) ||isempty(npxB)
730028      <a href="warndlg_uvmat.html" class="code" title="function hwarn=warndlg_uvmat(warntext,title)">warndlg_uvmat</a>(<span class="string">'The size of image A needs to be defined in the xml file ImaDoc'</span>,<span class="string">'ERROR'</span>)
740029      <span class="keyword">return</span>
750030  <span class="keyword">elseif</span> isempty(npxB) || isempty(npyB)
760031       <a href="warndlg_uvmat.html" class="code" title="function hwarn=warndlg_uvmat(warntext,title)">warndlg_uvmat</a>(<span class="string">'The size of image B needs to be defined in the xml file ImaDoc'</span>,<span class="string">'ERROR'</span>)
770032      <span class="keyword">return</span>
780033  <span class="keyword">end</span>
790034  <span class="keyword">if</span> isfield(XmlDataA,<span class="string">'GeometryCalib'</span>)
800035      tsaiA=XmlDataA.GeometryCalib;
810036  <span class="keyword">else</span>
820037      <a href="warndlg_uvmat.html" class="code" title="function hwarn=warndlg_uvmat(warntext,title)">warndlg_uvmat</a>(<span class="string">'no geometric calibration available for image A'</span>,<span class="string">'ERROR'</span>)
830038      <span class="keyword">return</span>
840039  <span class="keyword">end</span>
850040  <span class="keyword">if</span> isfield(XmlDataB,<span class="string">'GeometryCalib'</span>)
860041      tsaiB=XmlDataB.GeometryCalib;
870042  <span class="keyword">else</span>
880043      <a href="warndlg_uvmat.html" class="code" title="function hwarn=warndlg_uvmat(warntext,title)">warndlg_uvmat</a>(<span class="string">'no geometric calibration available for image B'</span>,<span class="string">'ERROR'</span>)
890044      <span class="keyword">return</span>
900045  <span class="keyword">end</span>
910046 
920047  <span class="comment">%corners of each image in real coordinates:</span>
930048  cornerA(:,1)=[0 0 npxA npxA]';<span class="comment">%x positions</span>
940049  cornerA(:,2)=[0 npyA 0 npyA]';<span class="comment">%y positions</span>
950050  cornerB(:,1)=[0 0 npxB npxB]';<span class="comment">%x positions</span>
960051  cornerB(:,2)=[0 npyB 0 npyB]';<span class="comment">%y positions</span>
970052 [xyA(:,1),xyA(:,2)]=<a href="phys_XYZ.html" class="code" title="function [Xphys,Yphys,Zphys]=phys_XYZ(Calib,X,Y,Z)">phys_XYZ</a>(tsaiA,cornerA(:,1),cornerA(:,2));
980053 [xyB(:,1),xyB(:,2)]=<a href="phys_XYZ.html" class="code" title="function [Xphys,Yphys,Zphys]=phys_XYZ(Calib,X,Y,Z)">phys_XYZ</a>(tsaiB,cornerB(:,1),cornerB(:,2));
990054  max_x=max(max(xyA(:,1)),max(xyB(:,1)));<span class="comment">%maximum on the 4 corners of the the images</span>
1000055  min_x=min(min(xyA(:,1)),min(xyB(:,1)));<span class="comment">%minimum on the 4 corners of the the images</span>
1010056  max_y=max(max(xyA(:,2)),max(xyB(:,2)));
1020057  min_y=min(min(xyA(:,2)),min(xyB(:,2)));
1030058  array_realx=[min_x:(max_x-min_x)/(nx_patch-1):max_x];
1040059  array_realy=[min_y:(max_y-min_y)/(ny_patch-1):max_y];
1050060  [grid_realx,grid_realy]=meshgrid(array_realx,array_realy);
1060061  grid_real(:,1)=reshape(grid_realx,nx_patch*ny_patch,1);
1070062  grid_real(:,2)=reshape(grid_realy,nx_patch*ny_patch,1);
1080063  grid_real(:,3)=zeros(nx_patch*ny_patch,1);
1090064 [grid_imaA(:,1),grid_imaA(:,2)]=<a href="px_XYZ.html" class="code" title="function [X,Y]=px_XYZ(Calib,Xphys,Yphys,Zphys)">px_XYZ</a>(tsaiA,grid_real(:,1),grid_real(:,2));
1100065 [grid_imaB(:,1),grid_imaB(:,2)]=<a href="px_XYZ.html" class="code" title="function [X,Y]=px_XYZ(Calib,Xphys,Yphys,Zphys)">px_XYZ</a>(tsaiB,grid_real(:,1),grid_real(:,2));
1110066
1120067  flagA=grid_imaA(:,1)&gt;0 &amp; grid_imaA(:,1)&lt;npxA &amp; grid_imaA(:,2)&gt;0 &amp; grid_imaA(:,2)&lt;npyA;
1130068  flagB=grid_imaB(:,1)&gt;0 &amp; grid_imaB(:,1)&lt;npxB &amp; grid_imaB(:,2)&gt;0 &amp; grid_imaB(:,2)&lt;npyB;
1140069  ind_good=find(flagA==1&amp;flagB==1);
1150070  XimaA=grid_imaA(ind_good,1);
1160071  YimaA=grid_imaA(ind_good,2);
1170072  XimaB=grid_imaB(ind_good,1);
1180073  YimaB=grid_imaB(ind_good,2);
1190074  grid_real_x=grid_real(ind_good,1);
1200075  grid_real_y=grid_real(ind_good,2);
1210076
1220077 <span class="comment">% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
1230078 <span class="comment">% %read the velocity fields</span>
1240079 <span class="comment">% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
1250080 <span class="comment">%</span>
1260081 <span class="comment">% [dt,time1,pixcmx,pixcmy,vec_X,vec_Y,vec_Z,vec_U,vec_V,vec_W,vec_C,vec_F,fixflag,vel_type_out,error,nb_coord,nb_dim]...</span>
1270082 <span class="comment">%     =read_vel({filecell_ncA},{vel_type});</span>
1280083 <span class="comment">%read field A</span>
1290084 [Field,VelTypeOut]=<a href="read_civxdata.html" class="code" title="function [Field,VelTypeOut]=read_civxdata(filename,FieldNames,VelType)">read_civxdata</a>(file_A,[],vel_type);
1300085 <span class="comment">%interpolate on XimaA</span>
1310086 Field.X=Field.X(find(Field.FF==0));
1320087 Field.Y=Field.Y(find(Field.FF==0));
1330088 Field.U=Field.U(find(Field.FF==0));
1340089 Field.V=Field.V(find(Field.FF==0));
1350090 dXa= <a href="griddata_uvmat.html" class="code" title="function ZI = griddata_uvmat(X,Y,Z,XI,YI)">griddata_uvmat</a>(Field.X,Field.Y,Field.U,XimaA,YimaA);
1360091 dYa= <a href="griddata_uvmat.html" class="code" title="function ZI = griddata_uvmat(X,Y,Z,XI,YI)">griddata_uvmat</a>(Field.X,Field.Y,Field.V,XimaA,YimaA);
1370092 dt=Field.dt;
1380093 time=Field.Time;
1390094
1400095 <span class="comment">%read field B</span>
1410096 <span class="comment">% [dt,time2,pixcmx,pixcmy,vec_X,vec_Y,vec_Z,vec_U,vec_V,vec_W,vec_C,vec_F,fixflag,vel_type_out,error,nb_coord,nb_dim]...</span>
1420097 <span class="comment">%     =read_vel({file_B},{vel_type});</span>
1430098 [Field,VelTypeOut]=<a href="read_civxdata.html" class="code" title="function [Field,VelTypeOut]=read_civxdata(filename,FieldNames,VelType)">read_civxdata</a>(file_B,FieldNames,vel_type);
1440099 <span class="keyword">if</span> ~isequal(Field.dt,dt)
1450100     <a href="warndlg_uvmat.html" class="code" title="function hwarn=warndlg_uvmat(warntext,title)">warndlg_uvmat</a>(<span class="string">'different time intervals for the two velocity fields '</span>,<span class="string">'ERROR'</span>)
1460101      <span class="keyword">return</span>
1470102 <span class="keyword">end</span>
1480103 <span class="keyword">if</span> ~isequal(Field.Time,time)
1490104     <a href="warndlg_uvmat.html" class="code" title="function hwarn=warndlg_uvmat(warntext,title)">warndlg_uvmat</a>(<span class="string">'different times for the two velocity fields '</span>,<span class="string">'ERROR'</span>)
1500105      <span class="keyword">return</span>
1510106 <span class="keyword">end</span>
1520107 <span class="comment">%interpolate on XimaB</span>
1530108 Field.X=Field.X(find(Field.FF==0));
1540109 Field.Y=Field.Y(find(Field.FF==0));
1550110 Field.U=Field.U(find(Field.FF==0));
1560111 Field.V=Field.V(find(Field.FF==0));
1570112 dXb=<a href="griddata_uvmat.html" class="code" title="function ZI = griddata_uvmat(X,Y,Z,XI,YI)">griddata_uvmat</a>(Field.X,Field.Y,Field.U,XimaB,YimaB);
1580113 dYb=<a href="griddata_uvmat.html" class="code" title="function ZI = griddata_uvmat(X,Y,Z,XI,YI)">griddata_uvmat</a>(Field.X,Field.Y,Field.V,XimaB,YimaB);
1590114 <span class="comment">%eliminate Not-a-Number</span>
1600115 ind_Nan=find(and(~isnan(dXa),~isnan(dXb)));
1610116 dXa=dXa(ind_Nan);
1620117 dYa=dYa(ind_Nan);
1630118 dXb=dXb(ind_Nan);
1640119 dYb=dYb(ind_Nan);
1650120 grid_phys1(:,1)=grid_real_x(ind_Nan);
1660121 grid_phys1(:,2)=grid_real_y(ind_Nan);
1670122 
1680123 <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
1690124 <span class="comment">%compute the coefficients</span>
1700125 <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
1710126
1720127
1730128 [A11,A12,A13,A21,A22,A23]=<a href="pxcm_tsai.html" class="code" title="function [A11,A12,A13,A21,A22,A23]=pxcm_tsai(a,var_phys)">pxcm_tsai</a>(tsaiA,grid_phys1);
1740129 [B11,B12,B13,B21,B22,B23]=<a href="pxcm_tsai.html" class="code" title="function [A11,A12,A13,A21,A22,A23]=pxcm_tsai(a,var_phys)">pxcm_tsai</a>(tsaiB,grid_phys1);
1750130
1760131 C1=A11.*A22-A12.*A21;
1770132 C2=A13.*A22-A12.*A23;
1780133 C3=A13.*A21-A11.*A23;
1790134 D1=B11.*B22-B12.*B21;
1800135 D2=B13.*B22-B12.*B23;
1810136 D3=B13.*B21-B11.*B23;
1820137 A1=(A22.*D1.*(C1.*D3-C3.*D1)+A21.*D1.*(C2.*D1-C1.*D2));
1830138 A2=(A12.*D1.*(C3.*D1-C1.*D3)+A11.*D1.*(C1.*D2-C2.*D1));
1840139 B1=(B22.*C1.*(C3.*D1-C1.*D3)+B21.*C1.*(C1.*D2-C2.*D1));
1850140 B2=(B12.*C1.*(C1.*D3-C3.*D1)+B11.*C1.*(C2.*D1-C1.*D2));
1860141 Lambda=(A1.*dXa+A2.*dYa+B1.*dXb+B2.*dYb)./(A1.*A1+A2.*A2+B1.*B1+B2.*B2);
1870142
1880143 <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
1890144 <span class="comment">%Projection for compatible displacements</span>
1900145 <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
1910146 Ua=dXa-Lambda.*A1;
1920147 Va=dYa-Lambda.*A2;
1930148 Ub=dXb-Lambda.*B1;
1940149 Vb=dYb-Lambda.*B2;
1950150
1960151 <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
1970152 <span class="comment">%Calculations of displacements and error</span>
1980153 <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
1990154 U=(A22.*D2.*Ua-A12.*D2.*Va-B22.*C2.*Ub+B12.*C2.*Vb)./(C1.*D2-C2.*D1);
2000155 V=(A21.*D3.*Ua-A11.*D3.*Va-B21.*C3.*Ub+B11.*C3.*Vb)./(C3.*D1-C1.*D3);
2010156 W=(A22.*D1.*Ua-A12.*D1.*Va-B22.*C1.*Ub+B12.*C1.*Vb)./(C2.*D1-C1.*D2);
2020157 W1=(-A21.*D1.*Ua+A11.*D1.*Va+B21.*C1.*Ub-B11.*C1.*Vb)./(C1.*D3-C3.*D1);
2030158
2040159 error=sqrt((A1.*dXa+A2.*dYa+B1.*dXb+B2.*dYb).*(A1.*dXa+A2.*dYa+B1.*dXb+B2.*dYb)./(A1.*A1+A2.*A2+B1.*B1+B2.*B2));
2050160
2060161 ind_error=(find(error&lt;thresh_patch));
2070162 U=U(ind_error);
2080163 V=V(ind_error);
2090164 W=W(ind_error);<span class="comment">%correction for water interface</span>
2100165 error=error(ind_error);
2110166
2120167 <span class="comment">%create nc grid file</span>
2130168 Result.ListGlobalAttribute={<span class="string">'nb_coord'</span>,<span class="string">'nb_dim'</span>,<span class="string">'constant_pixcm'</span>,<span class="string">'absolut_time_T0'</span>,<span class="string">'hart'</span>,<span class="string">'dt'</span>,<span class="string">'civ'</span>};
2140169 Result.nb_coord=3;<span class="comment">%grid file, no velocity</span>
2150170 Result.nb_dim=2;
2160171 Result.constant_pixcm=0;<span class="comment">%no linear correspondance with images</span>
2170172 Result.absolut_time_T0=time;<span class="comment">%absolute time of the field</span>
2180173 Result.hart=0;
2190174 Result.dt=dt;<span class="comment">%time interval for image correlation (put  by default)</span>
2200175 <span class="comment">% cte.title='grid';</span>
2210176 Result.civ=0;<span class="comment">%not a civ file (no direct correspondance with an image)</span>
2220177 Result.ListDimName={<span class="string">'nb_vectors'</span>}
2230178 Result.DimValue=length(U);
2240179 Result.ListVarName={<span class="string">'vec_X'</span>;<span class="string">'vec_Y'</span>;<span class="string">'vec_U'</span>;<span class="string">'vec_V'</span>;<span class="string">'vec_W'</span>;<span class="string">'vec_E'</span>};
2250180 Result.VarDimIndex: {[1]  [1]  [1]  [1]  [1]  [1]}
2260181 Result.vec_X= grid_phys1(ind_error,1);
2270182 Result.vec_Y= grid_phys1(ind_error,2);
2280183 Result.vec_U=U/dt;
2290184 Result.vec_V=V/dt;
2300185 Result.vec_W=W/dt;
2310186 Result.vec_E=error;
2320187 <span class="comment">% error=write_netcdf(file_st,cte,fieldlabels,grid_phys);</span>
2330188 error=<a href="struct2nc.html" class="code" title="function errormsg=struct2nc(flname,Data)">struct2nc</a>(file_st,Result);
2340189 display([file_st <span class="string">' written'</span>])
2350190
2360191
2370192
2380193
2390194
2400195
2410196
2420197
2430198
2440199
2450200
2460201
2470202</pre></div>
248<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>
249</body>
250</html>
Note: See TracBrowser for help on using the repository browser.