Changeset 209 for trunk/src/transform_field

Timestamp:

Feb 27, 2011, 10:43:18 PM (14 years ago)

Author:

sommeria

Message:

phys improved to deal with 3D fields

File:

• trunk/src/transform_field/phys.m (modified) (10 diffs)

trunk/src/transform_field/phys.m

@@ -94 +94 @@
 end
 
+% DataOut.VarDimName{2}
+% DataOut.VarDimName{3}
+% DataOut.VarDimName{4}
+% DataOut.VarDimName{5}
+% DataOut.VarDimName{6}
+% DataOut.VarDimName{7}
+% DataOut.VarAttribute{1}
+% DataOut.VarAttribute{2}
+% DataOut.VarAttribute{3}
+% DataOut.VarAttribute{4}
+% DataOut.VarAttribute{5}
+% DataOut.VarAttribute{6}
+% DataOut.VarAttribute{7}
 %------------------------------------------------
 function DataOut=phys_1(Data,Calib)
@@ -100 +113 @@
 DataOut=Data;%default
 % DataOut.CoordUnit=Calib.CoordUnit; %put flag for physical coordinates
-if isfield(Calib,'SliceCoord')
-    DataOut.PlaneCoord=Calib.SliceCoord;%to generalise for any plane 
+if isfield(Calib,'SliceCoord') && isfield(Data,'ZIndex')&&~isempty(Data.ZIndex)&&~isnan(Data.ZIndex)
+    DataOut.PlaneCoord=Calib.SliceCoord(Data.ZIndex,:);% transfer the slice position
+    if isfield(Calib,'SliceAngle') % transfer the slice rotation angles
+        DataOut.PlaneAngle=Calib.SliceAngle(Data.ZIndex,:);
+    end
 end
 % The transform ACTS ONLY IF .CoordType='px'and Calib defined
@@ -112 +128 @@
     DataOut.TimeUnit='s';
     %transform of X,Y coordinates for vector fields
+    test_z=0;
     if isfield(Data,'ZIndex') && ~isempty(Data.ZIndex)&&~isnan(Data.ZIndex)
         Z=Data.ZIndex;
+        test_z=1;
     else
         Z=0;
@@ -119 +137 @@
     if isfield(Data,'X') &&isfield(Data,'Y')&&~isempty(Data.X) && ~isempty(Data.Y)
         [DataOut.X,DataOut.Y,DataOut.Z]=phys_XYZ(Calib,Data.X,Data.Y,Z); 
+        if test_z
+             DataOut.ListVarName=[DataOut.ListVarName(1:2) {'Z'} DataOut.ListVarName(3:end)];
+             DataOut.VarDimName=[DataOut.VarDimName(1:2) DataOut.VarDimName(1) DataOut.VarDimName(3:end)];
+             ZAttribute{1}.Role='coord_z';
+             DataOut.VarAttribute=[DataOut.VarAttribute(1:2) ZAttribute DataOut.VarAttribute(3:end)];
+        end
         if isfield(Data,'U')&&isfield(Data,'V')&&~isempty(Data.U) && ~isempty(Data.V)&& isfield(Data,'dt') 
             if ~isempty(Data.dt)
@@ -209 +233 @@
         if numel(siz)==2 %(B/W images)
             vec_A=reshape(A{icell},1,npx(icell)*npy(icell));%put the original image in line
-            ind_in=find(flagin);
+            %ind_in=find(flagin);
             ind_out=find(~flagin);
             ICOMB=((XIMA-1)*npy(icell)+(npy(icell)+1-YIMA));
             ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A
-            vec_B(ind_in)=vec_A(ICOMB);
-            vec_B(ind_out)=zeros(size(ind_out));
+            %vec_B(ind_in)=vec_A(ICOMB);
+            vec_B(flagin)=vec_A(ICOMB);
+            vec_B(~flagin)=zeros(size(ind_out));
+%             vec_B(ind_out)=zeros(size(ind_out));
             A_out{icell}=reshape(vec_B,npY,npX);%new image in real coordinates
         elseif numel(siz)==3     
             for icolor=1:siz(3)
                 vec_A=reshape(A{icell}(:,:,icolor),1,npx*npy);%put the original image in line
-                ind_in=find(flagin);
+               % ind_in=find(flagin);
                 ind_out=find(~flagin);
                 ICOMB=((XIMA-1)*npy+(npy+1-YIMA));
                 ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A
-                vec_B(ind_in)=vec_A(ICOMB);
-                vec_B(ind_out)=zeros(size(ind_out));
+                vec_B(flagin)=vec_A(ICOMB);
+                vec_B(~flagin)=zeros(size(ind_out));
                 A_out{icell}(:,:,icolor)=reshape(vec_B,npy,npx);%new image in real coordinates
             end
@@ -251 +277 @@
 %INPUT:
 %Z: index of plane
-function [Xphys,Yphys,Zphys]=phys_XYZ(Calib,X,Y,Z)
+function [Xphys,Yphys,Zphys]=phys_XYZ(Calib,X,Y,Zindex)
 %------------------------------------------------------------------------
-if exist('Z','var')&& isequal(Z,round(Z))&& Z>0 && isfield(Calib,'SliceCoord')&&length(Calib.SliceCoord)>=Z
-    Zindex=Z;
-    Zphys=Calib.SliceCoord(Zindex,3);%GENERALISER AUX CAS AVEC ANGLE
+testangle=0;
+test_refraction=0;
+if exist('Zindex','var')&& isequal(Zindex,round(Zindex))&& Zindex>0 && isfield(Calib,'SliceCoord')&&length(Calib.SliceCoord)>=Zindex
+    if isfield(Calib, 'SliceAngle') && ~isequal(Calib.SliceAngle,0)
+        testangle=1;
+        om=norm(Calib.SliceAngle(Zindex,:));%norm of rotation angle in radians
+        OmAxis=Calib.SliceAngle(Zindex,:)/om; %unit vector marking the rotation axis
+        cos_om=cos(pi*om/180);
+        sin_om=sin(pi*om/180);
+        coeff=OmAxis(3)*(1-cos_om);
+        norm_plane(1)=OmAxis(1)*coeff+OmAxis(2)*sin_om;
+        norm_plane(2)=OmAxis(2)*coeff-OmAxis(1)*sin_om;
+        norm_plane(3)=OmAxis(3)*coeff+cos_om;
+        Z0=norm_plane*Calib.SliceCoord(Zindex,:)'/norm_plane(3);
+    else
+        Z0=Calib.SliceCoord(Zindex,3);%horizontal plane z=cte
+    end
+    Z0virt=Z0;
     if isfield(Calib,'InterfaceCoord') && isfield(Calib,'RefractionIndex') 
         H=Calib.InterfaceCoord(3);
-        if H>Zphys
-            Zphys=H-(H-Zphys)/Calib.RefractionIndex; %corrected z (virtual object)
+        if H>Z0
+            Z0virt=H-(H-Z0)/Calib.RefractionIndex; %corrected z (virtual object)
+            test_refraction=1;
         end
     end   
 else
-    Zphys=0;
+    Z0=0;
+    Z0virt=0;
 end
 if ~exist('X','var')||~exist('Y','var')
@@ -285 +328 @@
 if isfield(Calib,'R')
     R=(Calib.R)';
+    if testangle
+        a=-norm_plane(1)/norm_plane(3);
+        b=-norm_plane(2)/norm_plane(3);
+        if test_refraction
+            a=a/Calib.RefractionIndex;
+            b=b/Calib.RefractionIndex;
+        end
+        R(1)=R(1)+a*R(3);
+        R(2)=R(2)+b*R(3);
+        R(4)=R(4)+a*R(6);
+        R(5)=R(5)+b*R(6);
+        R(7)=R(7)+a*R(9);
+        R(8)=R(8)+b*R(9);
+    end
     Tx=Calib.Tx_Ty_Tz(1);
     Ty=Calib.Tx_Ty_Tz(2);
@@ -299 +356 @@
     Zx0=R(3)*R(5)-R(2)*R(6);
     Zy0=R(1)*R(6)-R(3)*R(4);
-    A11=R(8)*Ty-R(5)*Tz+Z11*Zphys;
-    A12=R(2)*Tz-R(8)*Tx+Z12*Zphys;
-    A21=-R(7)*Ty+R(4)*Tz+Z21*Zphys;
-    A22=-R(1)*Tz+R(7)*Tx+Z22*Zphys;
-    X0=f*(R(5)*Tx-R(2)*Ty+Zx0*Zphys);
-    Y0=f*(-R(4)*Tx+R(1)*Ty+Zy0*Zphys);
+    A11=R(8)*Ty-R(5)*Tz+Z11*Z0virt;
+    A12=R(2)*Tz-R(8)*Tx+Z12*Z0virt;
+    A21=-R(7)*Ty+R(4)*Tz+Z21*Z0virt;
+    A22=-R(1)*Tz+R(7)*Tx+Z22*Z0virt;
+    X0=f*(R(5)*Tx-R(2)*Ty+Zx0*Z0virt);
+    Y0=f*(-R(4)*Tx+R(1)*Ty+Zy0*Z0virt);
         %px to camera:
     Xd=dpx*(X-Calib.Cx_Cy(1)); % sensor coordinates
@@ -312 +369 @@
     Yu=Yd./dist_fact;
     denom=Dx*Xu+Dy*Yu+D0;
-    % denom2=denom.*denom;
     Xphys=(A11.*Xu+A12.*Yu+X0)./denom;%world coordinates
     Yphys=(A21.*Xu+A22.*Yu+Y0)./denom;
-%     Xd=(X-Calib.Cx_Cy(1))/Calib.fx_fy(1); % sensor coordinates
-%     Yd=(Y-Calib.Cx_Cy(2))/Calib.fx_fy(2);
-%     dist_fact=1+Calib.kc*(Xd.*Xd+Yd.*Yd); %distortion factor
-%     Xu=Xd./dist_fact;%undistorted sensor coordinates
-%     Yu=Yd./dist_fact;
-%     A11=R(7)*Xu-R(1);
-%     A12=R(8)*Xu-R(2);
-%     A21=R(7)*Yu-R(4);
-%     A22=R(8)*Yu-R(5);
-%     B1=(R(3)-R(9)*Xu)*Zphys-Tz*Xu+Tx;
-%     B2=(R(6)-R(9)*Yu)*Zphys-Tz*Yu+Ty;
-%     deter=A12.*A21-A11.*A22;
-%     Xphys=(A21.*B1-A11.*B2)./deter;
-%     Yphys=(-A22.*B1+A12.*B2)./deter;
+    if testangle
+        Zphys=Z0+a*Xphys+b*Yphys;
+    else
+        Zphys=Z0;
+    end
 else
     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)
 % Xphys, Yphys: array of x,y physical coordinates
-% [Zphys]: corresponding array of z physical coordinates (0 by default)
-
-
-
-
+% [Z0]: corresponding array of z physical coordinates (0 by default)
+
+
+
+