Changeset 209 for trunk/src/transform_field/phys.m
 Timestamp:
 Feb 27, 2011, 10:43:18 PM (13 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/src/transform_field/phys.m
r202 r209 94 94 end 95 95 96 % DataOut.VarDimName{2} 97 % DataOut.VarDimName{3} 98 % DataOut.VarDimName{4} 99 % DataOut.VarDimName{5} 100 % DataOut.VarDimName{6} 101 % DataOut.VarDimName{7} 102 % DataOut.VarAttribute{1} 103 % DataOut.VarAttribute{2} 104 % DataOut.VarAttribute{3} 105 % DataOut.VarAttribute{4} 106 % DataOut.VarAttribute{5} 107 % DataOut.VarAttribute{6} 108 % DataOut.VarAttribute{7} 96 109 % 97 110 function DataOut=phys_1(Data,Calib) … … 100 113 DataOut=Data;%default 101 114 % DataOut.CoordUnit=Calib.CoordUnit; %put flag for physical coordinates 102 if isfield(Calib,'SliceCoord') 103 DataOut.PlaneCoord=Calib.SliceCoord;%to generalise for any plane 115 if isfield(Calib,'SliceCoord') && isfield(Data,'ZIndex')&&~isempty(Data.ZIndex)&&~isnan(Data.ZIndex) 116 DataOut.PlaneCoord=Calib.SliceCoord(Data.ZIndex,:);% transfer the slice position 117 if isfield(Calib,'SliceAngle') % transfer the slice rotation angles 118 DataOut.PlaneAngle=Calib.SliceAngle(Data.ZIndex,:); 119 end 104 120 end 105 121 % The transform ACTS ONLY IF .CoordType='px'and Calib defined … … 112 128 DataOut.TimeUnit='s'; 113 129 %transform of X,Y coordinates for vector fields 130 test_z=0; 114 131 if isfield(Data,'ZIndex') && ~isempty(Data.ZIndex)&&~isnan(Data.ZIndex) 115 132 Z=Data.ZIndex; 133 test_z=1; 116 134 else 117 135 Z=0; … … 119 137 if isfield(Data,'X') &&isfield(Data,'Y')&&~isempty(Data.X) && ~isempty(Data.Y) 120 138 [DataOut.X,DataOut.Y,DataOut.Z]=phys_XYZ(Calib,Data.X,Data.Y,Z); 139 if test_z 140 DataOut.ListVarName=[DataOut.ListVarName(1:2) {'Z'} DataOut.ListVarName(3:end)]; 141 DataOut.VarDimName=[DataOut.VarDimName(1:2) DataOut.VarDimName(1) DataOut.VarDimName(3:end)]; 142 ZAttribute{1}.Role='coord_z'; 143 DataOut.VarAttribute=[DataOut.VarAttribute(1:2) ZAttribute DataOut.VarAttribute(3:end)]; 144 end 121 145 if isfield(Data,'U')&&isfield(Data,'V')&&~isempty(Data.U) && ~isempty(Data.V)&& isfield(Data,'dt') 122 146 if ~isempty(Data.dt) … … 209 233 if numel(siz)==2 %(B/W images) 210 234 vec_A=reshape(A{icell},1,npx(icell)*npy(icell));%put the original image in line 211 ind_in=find(flagin);235 %ind_in=find(flagin); 212 236 ind_out=find(~flagin); 213 237 ICOMB=((XIMA1)*npy(icell)+(npy(icell)+1YIMA)); 214 238 ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A 215 vec_B(ind_in)=vec_A(ICOMB); 216 vec_B(ind_out)=zeros(size(ind_out)); 239 %vec_B(ind_in)=vec_A(ICOMB); 240 vec_B(flagin)=vec_A(ICOMB); 241 vec_B(~flagin)=zeros(size(ind_out)); 242 % vec_B(ind_out)=zeros(size(ind_out)); 217 243 A_out{icell}=reshape(vec_B,npY,npX);%new image in real coordinates 218 244 elseif numel(siz)==3 219 245 for icolor=1:siz(3) 220 246 vec_A=reshape(A{icell}(:,:,icolor),1,npx*npy);%put the original image in line 221 ind_in=find(flagin);247 % ind_in=find(flagin); 222 248 ind_out=find(~flagin); 223 249 ICOMB=((XIMA1)*npy+(npy+1YIMA)); 224 250 ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A 225 vec_B( ind_in)=vec_A(ICOMB);226 vec_B( ind_out)=zeros(size(ind_out));251 vec_B(flagin)=vec_A(ICOMB); 252 vec_B(~flagin)=zeros(size(ind_out)); 227 253 A_out{icell}(:,:,icolor)=reshape(vec_B,npy,npx);%new image in real coordinates 228 254 end … … 251 277 %INPUT: 252 278 %Z: index of plane 253 function [Xphys,Yphys,Zphys]=phys_XYZ(Calib,X,Y,Z )279 function [Xphys,Yphys,Zphys]=phys_XYZ(Calib,X,Y,Zindex) 254 280 % 255 if exist('Z','var')&& isequal(Z,round(Z))&& Z>0 && isfield(Calib,'SliceCoord')&&length(Calib.SliceCoord)>=Z 256 Zindex=Z; 257 Zphys=Calib.SliceCoord(Zindex,3);%GENERALISER AUX CAS AVEC ANGLE 281 testangle=0; 282 test_refraction=0; 283 if exist('Zindex','var')&& isequal(Zindex,round(Zindex))&& Zindex>0 && isfield(Calib,'SliceCoord')&&length(Calib.SliceCoord)>=Zindex 284 if isfield(Calib, 'SliceAngle') && ~isequal(Calib.SliceAngle,0) 285 testangle=1; 286 om=norm(Calib.SliceAngle(Zindex,:));%norm of rotation angle in radians 287 OmAxis=Calib.SliceAngle(Zindex,:)/om; %unit vector marking the rotation axis 288 cos_om=cos(pi*om/180); 289 sin_om=sin(pi*om/180); 290 coeff=OmAxis(3)*(1cos_om); 291 norm_plane(1)=OmAxis(1)*coeff+OmAxis(2)*sin_om; 292 norm_plane(2)=OmAxis(2)*coeffOmAxis(1)*sin_om; 293 norm_plane(3)=OmAxis(3)*coeff+cos_om; 294 Z0=norm_plane*Calib.SliceCoord(Zindex,:)'/norm_plane(3); 295 else 296 Z0=Calib.SliceCoord(Zindex,3);%horizontal plane z=cte 297 end 298 Z0virt=Z0; 258 299 if isfield(Calib,'InterfaceCoord') && isfield(Calib,'RefractionIndex') 259 300 H=Calib.InterfaceCoord(3); 260 if H>Zphys 261 Zphys=H(HZphys)/Calib.RefractionIndex; %corrected z (virtual object) 301 if H>Z0 302 Z0virt=H(HZ0)/Calib.RefractionIndex; %corrected z (virtual object) 303 test_refraction=1; 262 304 end 263 305 end 264 306 else 265 Zphys=0; 307 Z0=0; 308 Z0virt=0; 266 309 end 267 310 if ~exist('X','var')~exist('Y','var') … … 285 328 if isfield(Calib,'R') 286 329 R=(Calib.R)'; 330 if testangle 331 a=norm_plane(1)/norm_plane(3); 332 b=norm_plane(2)/norm_plane(3); 333 if test_refraction 334 a=a/Calib.RefractionIndex; 335 b=b/Calib.RefractionIndex; 336 end 337 R(1)=R(1)+a*R(3); 338 R(2)=R(2)+b*R(3); 339 R(4)=R(4)+a*R(6); 340 R(5)=R(5)+b*R(6); 341 R(7)=R(7)+a*R(9); 342 R(8)=R(8)+b*R(9); 343 end 287 344 Tx=Calib.Tx_Ty_Tz(1); 288 345 Ty=Calib.Tx_Ty_Tz(2); … … 299 356 Zx0=R(3)*R(5)R(2)*R(6); 300 357 Zy0=R(1)*R(6)R(3)*R(4); 301 A11=R(8)*TyR(5)*Tz+Z11*Z phys;302 A12=R(2)*TzR(8)*Tx+Z12*Z phys;303 A21=R(7)*Ty+R(4)*Tz+Z21*Z phys;304 A22=R(1)*Tz+R(7)*Tx+Z22*Z phys;305 X0=f*(R(5)*TxR(2)*Ty+Zx0*Z phys);306 Y0=f*(R(4)*Tx+R(1)*Ty+Zy0*Z phys);358 A11=R(8)*TyR(5)*Tz+Z11*Z0virt; 359 A12=R(2)*TzR(8)*Tx+Z12*Z0virt; 360 A21=R(7)*Ty+R(4)*Tz+Z21*Z0virt; 361 A22=R(1)*Tz+R(7)*Tx+Z22*Z0virt; 362 X0=f*(R(5)*TxR(2)*Ty+Zx0*Z0virt); 363 Y0=f*(R(4)*Tx+R(1)*Ty+Zy0*Z0virt); 307 364 %px to camera: 308 365 Xd=dpx*(XCalib.Cx_Cy(1)); % sensor coordinates … … 312 369 Yu=Yd./dist_fact; 313 370 denom=Dx*Xu+Dy*Yu+D0; 314 % denom2=denom.*denom;315 371 Xphys=(A11.*Xu+A12.*Yu+X0)./denom;%world coordinates 316 372 Yphys=(A21.*Xu+A22.*Yu+Y0)./denom; 317 % Xd=(XCalib.Cx_Cy(1))/Calib.fx_fy(1); % sensor coordinates 318 % Yd=(YCalib.Cx_Cy(2))/Calib.fx_fy(2); 319 % dist_fact=1+Calib.kc*(Xd.*Xd+Yd.*Yd); %distortion factor 320 % Xu=Xd./dist_fact;%undistorted sensor coordinates 321 % Yu=Yd./dist_fact; 322 % A11=R(7)*XuR(1); 323 % A12=R(8)*XuR(2); 324 % A21=R(7)*YuR(4); 325 % A22=R(8)*YuR(5); 326 % B1=(R(3)R(9)*Xu)*ZphysTz*Xu+Tx; 327 % B2=(R(6)R(9)*Yu)*ZphysTz*Yu+Ty; 328 % deter=A12.*A21A11.*A22; 329 % Xphys=(A21.*B1A11.*B2)./deter; 330 % Yphys=(A22.*B1+A12.*B2)./deter; 373 if testangle 374 Zphys=Z0+a*Xphys+b*Yphys; 375 else 376 Zphys=Z0; 377 end 331 378 else 332 379 Xphys=Calib.Tx_Ty_Tz(1)+X/Calib.fx_fy(1); … … 343 390 % Calib: structure containing the calibration parameters (read from the ImaDoc .xml file) 344 391 % Xphys, Yphys: array of x,y physical coordinates 345 % [Z phys]: corresponding array of z physical coordinates (0 by default)346 347 348 349 392 % [Z0]: corresponding array of z physical coordinates (0 by default) 393 394 395 396
Note: See TracChangeset
for help on using the changeset viewer.