Changeset 116 for trunk/src/transform_field/phys.m
 Timestamp:
 Oct 16, 2010, 5:56:00 PM (13 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/src/transform_field/phys.m
r79 r116 58 58 end 59 59 %transform of X,Y coordinates for vector fields 60 if isfield(Data,'ZIndex')&&~isempty(Data.ZIndex) 60 if isfield(Data,'ZIndex')&&~isempty(Data.ZIndex)&&~isnan(Data.ZIndex) 61 61 ZIndex=Data.ZIndex; 62 62 else 63 ZIndex= 0;63 ZIndex=1; 64 64 end 65 65 if test_1 … … 111 111 DataOut.TimeUnit='s'; 112 112 %transform of X,Y coordinates for vector fields 113 if isfield(Data,'ZIndex') && ~isempty(Data.ZIndex) 114 Z=Data.ZIndex ;113 if isfield(Data,'ZIndex') && ~isempty(Data.ZIndex)&&~isnan(Data.ZIndex) 114 Z=Data.ZIndex 115 115 else 116 116 Z=0; … … 153 153 end 154 154 155 156 155 %%%%%%%%%%%%%%%%%%%% 157 156 function [A_out,Rangx,Rangy]=phys_Ima(A,CalibIn,ZIndex) 157 158 if ndims(A)==3 159 A=mean(A,3); 160 end 161 162 158 163 xcorner=[]; 159 164 ycorner=[]; 160 165 npx=[]; 161 166 npy=[]; 167 dx=ones(1,length(A)); 168 dy=ones(1,length(A)); 162 169 for icell=1:length(A) 163 170 siz=size(A{icell}); … … 165 172 npy=[npy siz(1)]; 166 173 Calib=CalibIn{icell}; 167 xima=[0.5 siz(2)0.5 0.5 siz(2)0.5];%image coordi antes of corners174 xima=[0.5 siz(2)0.5 0.5 siz(2)0.5];%image coordinates of corners 168 175 yima=[0.5 0.5 siz(1)0.5 siz(1)0.5]; 169 176 [xcorner_new,ycorner_new]=phys_XYZ(Calib,xima,yima,ZIndex);%corresponding physical coordinates … … 177 184 Rangy(2)=min(ycorner); 178 185 Rangy(1)=max(ycorner); 179 test_multi=(max(npx)~=min(npx))  (max(npy)~=min(npy));186 test_multi=(max(npx)~=min(npx))  (max(npy)~=min(npy)); %different image lengths 180 187 npX=1+round((Rangx(2)Rangx(1))/min(dx));% nbre of pixels in the new image (use the finest resolution min(dx) in the set of images) 181 188 npY=1+round((Rangy(1)Rangy(2))/min(dy)); … … 187 194 for icell=1:length(A) 188 195 Calib=CalibIn{icell}; 189 if (isfield(Calib,'R') && ~isequal(Calib.R(2,1),0) && ~isequal(Calib.R(1,2),0)) ... 190 ((isfield(Calib,'kappa1')&& ~isequal(Calib.kappa1,0)))  test_multi  ~isequal(Calib,CalibIn{1}) 196 if isfield(Calib,'R')  isfield(Calib,'kc') test_multi ~isequal(Calib,CalibIn{1})% the image needs to be interpolated to the new coordinates 191 197 zphys=0; %default 192 198 if isfield(Calib,'SliceCoord') %.Z= index of plane … … 227 233 A_out{icell}=uint16(A_out{icell}); 228 234 end 229 else% 230 235 else% 231 236 A_out{icell}=A{icell};%no transform 232 237 Rangx=[0.5 npx0.5];%image coordiantes of corners 233 238 Rangy=[npy0.5 0.5]; 234 [Rangx]=phys_XYZ(Calib,Rangx,[0.5 0.5],[ZIndex ZIndex]);%case of translations without rotation and quadratic deformation 235 [xx,Rangy]=phys_XYZ(Calib,[0.5 0.5],Rangy,[ZIndex ZIndex]); 236 end 237 end 238 239 [Rangx]=phys_XYZ(Calib,Rangx,[0.5 0.5],ZIndex);%case of translations without rotation and quadratic deformation 240 [xx,Rangy]=phys_XYZ(Calib,[0.5 0.5],Rangy,ZIndex); 241 end 242 end 243 244 % 239 245 %'phys_XYZ':transforms image (px) to real world (phys) coordinates using geometric calibration parameters 240 246 % function [Xphys,Yphys]=phys_XYZ(Calib,X,Y,Z) … … 245 251 %Z: index of plane 246 252 function [Xphys,Yphys,Zphys]=phys_XYZ(Calib,X,Y,Z) 247 if exist('Z','var')& isequal(Z,round(Z))& Z>0 & isfield(Calib,'SliceCoord')&length(Calib.SliceCoord)>=Z 253 % 254 if exist('Z','var')&& isequal(Z,round(Z))&& Z>0 && isfield(Calib,'SliceCoord')&&length(Calib.SliceCoord)>=Z 248 255 Zindex=Z; 249 256 Zphys=Calib.SliceCoord(Zindex,3);%GENERALISER AUX CAS AVEC ANGLE 250 257 else 251 % if exist('Z','var') 252 % Zphys=Z; 253 % else 254 Zphys=0; 255 % end 258 Zphys=0; 256 259 end 257 260 if ~exist('X','var')~exist('Y','var') … … 260 263 return 261 264 end 262 Xphys=X;%default 263 Yphys=Y; 264 %image transform 265 %coordinate transform 266 if ~isfield(Calib,'fx_fy') 267 Calib.fx_fy=[1 1]; 268 end 269 if ~isfield(Calib,'Tx_Ty_Tz') 270 Calib.Tx_Ty_Tz=[0 0 1]; 271 end 272 if ~isfield(Calib,'Cx_Cy') 273 Calib.Cx_Cy=[0 0]; 274 end 275 if ~isfield(Calib,'kc') 276 Calib.kc=0; 277 end 265 278 if isfield(Calib,'R') 266 279 R=(Calib.R)'; 280 Tx=Calib.Tx_Ty_Tz(1); 281 Ty=Calib.Tx_Ty_Tz(2); 282 Tz=Calib.Tx_Ty_Tz(3); 283 f=Calib.fx_fy(1);%dpy=1; sx=1 284 dpx=Calib.fx_fy(2)/Calib.fx_fy(1); 267 285 Dx=R(5)*R(7)R(4)*R(8); 268 286 Dy=R(1)*R(8)R(2)*R(7); 269 D0= Calib.f*(R(2)*R(4)R(1)*R(5));287 D0=f*(R(2)*R(4)R(1)*R(5)); 270 288 Z11=R(6)*R(8)R(5)*R(9); 271 289 Z12=R(2)*R(9)R(3)*R(8); … … 274 292 Zx0=R(3)*R(5)R(2)*R(6); 275 293 Zy0=R(1)*R(6)R(3)*R(4); 276 A11=R(8)* Calib.TyR(5)*Calib.Tz+Z11*Zphys;277 A12=R(2)* Calib.TzR(8)*Calib.Tx+Z12*Zphys;278 A21=R(7)* Calib.Ty+R(4)*Calib.Tz+Z21*Zphys;279 A22=R(1)* Calib.Tz+R(7)*Calib.Tx+Z11*Zphys;280 X0= Calib.f*(R(5)*Calib.TxR(2)*Calib.Ty+Zx0*Zphys);281 Y0= Calib.f*(R(4)*Calib.Tx+R(1)*Calib.Ty+Zy0*Zphys);294 A11=R(8)*TyR(5)*Tz+Z11*Zphys; 295 A12=R(2)*TzR(8)*Tx+Z12*Zphys; 296 A21=R(7)*Ty+R(4)*Tz+Z21*Zphys; 297 A22=R(1)*Tz+R(7)*Tx+Z22*Zphys; 298 X0=f*(R(5)*TxR(2)*Ty+Zx0*Zphys); 299 Y0=f*(R(4)*Tx+R(1)*Ty+Zy0*Zphys); 282 300 %px to camera: 283 Xd= (Calib.dpx/Calib.sx)*(XCalib.Cx); % sensor coordinates284 Yd= Calib.dpy*(YCalib.Cy);285 dist_fact=1+Calib.k appa1*(Xd.*Xd+Yd.*Yd); %distortion factor286 Xu= dist_fact.*Xd;%undistorted sensor coordinates287 Yu= dist_fact.*Yd;301 Xd=dpx*(XCalib.Cx_Cy(1)); % sensor coordinates 302 Yd=(YCalib.Cx_Cy(2)); 303 dist_fact=1+Calib.kc*(Xd.*Xd+Yd.*Yd)/(f*f); %distortion factor 304 Xu=Xd./dist_fact;%undistorted sensor coordinates 305 Yu=Yd./dist_fact; 288 306 denom=Dx*Xu+Dy*Yu+D0; 289 307 % denom2=denom.*denom; 290 308 Xphys=(A11.*Xu+A12.*Yu+X0)./denom;%world coordinates 291 309 Yphys=(A21.*Xu+A22.*Yu+Y0)./denom; 310 % Xd=(XCalib.Cx_Cy(1))/Calib.fx_fy(1); % sensor coordinates 311 % Yd=(YCalib.Cx_Cy(2))/Calib.fx_fy(2); 312 % dist_fact=1+Calib.kc*(Xd.*Xd+Yd.*Yd); %distortion factor 313 % Xu=Xd./dist_fact;%undistorted sensor coordinates 314 % Yu=Yd./dist_fact; 315 % A11=R(7)*XuR(1); 316 % A12=R(8)*XuR(2); 317 % A21=R(7)*YuR(4); 318 % A22=R(8)*YuR(5); 319 % B1=(R(3)R(9)*Xu)*ZphysTz*Xu+Tx; 320 % B2=(R(6)R(9)*Yu)*ZphysTz*Yu+Ty; 321 % deter=A12.*A21A11.*A22; 322 % Xphys=(A21.*B1A11.*B2)./deter; 323 % Yphys=(A22.*B1+A12.*B2)./deter; 324 else 325 Xphys=Calib.Tx_Ty_Tz(1)+X/Calib.fx_fy(1); 326 Yphys=Calib.Tx_Ty_Tz(2)+Y/Calib.fx_fy(2); 292 327 end 293 328
Note: See TracChangeset
for help on using the changeset viewer.