Ignore:
Timestamp:
Oct 16, 2010, 5:56:00 PM (14 years ago)
Author:
sommeria
Message:

geometry_calib is now updated when a new image is viewed by uvmat
imadoc2struct corrected for time reading

Location:
trunk/src/transform_field
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/transform_field/phys.m

    r79 r116  
    5858end
    5959%transform of X,Y coordinates for vector fields
    60 if isfield(Data,'ZIndex')&&~isempty(Data.ZIndex)
     60if isfield(Data,'ZIndex')&&~isempty(Data.ZIndex)&&~isnan(Data.ZIndex)
    6161    ZIndex=Data.ZIndex;
    6262else
    63     ZIndex=0;
     63    ZIndex=1;
    6464end
    6565if test_1
    … …  
    111111    DataOut.TimeUnit='s';
    112112    %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
    115115    else
    116116        Z=0;
    … …  
    153153end
    154154
    155 
    156155%%%%%%%%%%%%%%%%%%%%
    157156function [A_out,Rangx,Rangy]=phys_Ima(A,CalibIn,ZIndex)
     157
     158if ndims(A)==3
     159    A=mean(A,3);
     160end
     161
     162
    158163xcorner=[];
    159164ycorner=[];
    160165npx=[];
    161166npy=[];
     167dx=ones(1,length(A));
     168dy=ones(1,length(A));
    162169for icell=1:length(A)
    163170    siz=size(A{icell});
    … …  
    165172    npy=[npy siz(1)];
    166173    Calib=CalibIn{icell};
    167     xima=[0.5 siz(2)-0.5 0.5 siz(2)-0.5];%image coordiantes of corners
     174    xima=[0.5 siz(2)-0.5 0.5 siz(2)-0.5];%image coordinates of corners
    168175    yima=[0.5 0.5 siz(1)-0.5 siz(1)-0.5];
    169176    [xcorner_new,ycorner_new]=phys_XYZ(Calib,xima,yima,ZIndex);%corresponding physical coordinates
    … …  
    177184Rangy(2)=min(ycorner);
    178185Rangy(1)=max(ycorner);
    179 test_multi=(max(npx)~=min(npx)) | (max(npy)~=min(npy));
     186test_multi=(max(npx)~=min(npx)) || (max(npy)~=min(npy)); %different image lengths
    180187npX=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)
    181188npY=1+round((Rangy(1)-Rangy(2))/min(dy));
    … …  
    187194for icell=1:length(A)
    188195    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
    191197        zphys=0; %default
    192198        if isfield(Calib,'SliceCoord') %.Z= index of plane
    … …  
    227233            A_out{icell}=uint16(A_out{icell});
    228234        end
    229     else%
    230        
     235    else%     
    231236        A_out{icell}=A{icell};%no transform
    232237        Rangx=[0.5 npx-0.5];%image coordiantes of corners
    233238        Rangy=[npy-0.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
     242end
     243
     244%------------------------------------------------------------------------
    239245%'phys_XYZ':transforms image (px) to real world (phys) coordinates using geometric calibration parameters
    240246% function [Xphys,Yphys]=phys_XYZ(Calib,X,Y,Z)
    … …  
    245251%Z: index of plane
    246252function [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%------------------------------------------------------------------------
     254if exist('Z','var')&& isequal(Z,round(Z))&& Z>0 && isfield(Calib,'SliceCoord')&&length(Calib.SliceCoord)>=Z
    248255    Zindex=Z;
    249256    Zphys=Calib.SliceCoord(Zindex,3);%GENERALISER AUX CAS AVEC ANGLE
    250257else
    251 %     if exist('Z','var')
    252 %         Zphys=Z;
    253 %     else
    254         Zphys=0;
    255 %     end
     258    Zphys=0;
    256259end
    257260if ~exist('X','var')||~exist('Y','var')
    … …  
    260263    return
    261264end
    262 Xphys=X;%default
    263 Yphys=Y;
    264 %image transform
     265%coordinate transform
     266if ~isfield(Calib,'fx_fy')
     267     Calib.fx_fy=[1 1];
     268end
     269if ~isfield(Calib,'Tx_Ty_Tz')
     270     Calib.Tx_Ty_Tz=[0 0 1];
     271end
     272if ~isfield(Calib,'Cx_Cy')
     273     Calib.Cx_Cy=[0 0];
     274end
     275if ~isfield(Calib,'kc')
     276     Calib.kc=0;
     277end
    265278if isfield(Calib,'R')
    266279    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);
    267285    Dx=R(5)*R(7)-R(4)*R(8);
    268286    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));
    270288    Z11=R(6)*R(8)-R(5)*R(9);
    271289    Z12=R(2)*R(9)-R(3)*R(8); 
    … …  
    274292    Zx0=R(3)*R(5)-R(2)*R(6);
    275293    Zy0=R(1)*R(6)-R(3)*R(4);
    276     A11=R(8)*Calib.Ty-R(5)*Calib.Tz+Z11*Zphys;
    277     A12=R(2)*Calib.Tz-R(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.Tx-R(2)*Calib.Ty+Zx0*Zphys);
    281     Y0=Calib.f*(-R(4)*Calib.Tx+R(1)*Calib.Ty+Zy0*Zphys);
     294    A11=R(8)*Ty-R(5)*Tz+Z11*Zphys;
     295    A12=R(2)*Tz-R(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)*Tx-R(2)*Ty+Zx0*Zphys);
     299    Y0=f*(-R(4)*Tx+R(1)*Ty+Zy0*Zphys);
    282300        %px to camera:
    283     Xd=(Calib.dpx/Calib.sx)*(X-Calib.Cx); % sensor coordinates
    284     Yd=Calib.dpy*(Y-Calib.Cy);
    285     dist_fact=1+Calib.kappa1*(Xd.*Xd+Yd.*Yd); %distortion factor
    286     Xu=dist_fact.*Xd;%undistorted sensor coordinates
    287     Yu=dist_fact.*Yd;
     301    Xd=dpx*(X-Calib.Cx_Cy(1)); % sensor coordinates
     302    Yd=(Y-Calib.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;
    288306    denom=Dx*Xu+Dy*Yu+D0;
    289307    % denom2=denom.*denom;
    290308    Xphys=(A11.*Xu+A12.*Yu+X0)./denom;%world coordinates
    291309    Yphys=(A21.*Xu+A22.*Yu+Y0)./denom;
     310%     Xd=(X-Calib.Cx_Cy(1))/Calib.fx_fy(1); % sensor coordinates
     311%     Yd=(Y-Calib.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)*Xu-R(1);
     316%     A12=R(8)*Xu-R(2);
     317%     A21=R(7)*Yu-R(4);
     318%     A22=R(8)*Yu-R(5);
     319%     B1=(R(3)-R(9)*Xu)*Zphys-Tz*Xu+Tx;
     320%     B2=(R(6)-R(9)*Yu)*Zphys-Tz*Yu+Ty;
     321%     deter=A12.*A21-A11.*A22;
     322%     Xphys=(A21.*B1-A11.*B2)./deter;
     323%     Yphys=(-A22.*B1+A12.*B2)./deter;
     324else
     325    Xphys=-Calib.Tx_Ty_Tz(1)+X/Calib.fx_fy(1);
     326    Yphys=-Calib.Tx_Ty_Tz(2)+Y/Calib.fx_fy(2);
    292327end
    293328

  • trunk/src/transform_field/px.m

    r40 r116  
    1818%      DataIn.A, AX, AY : image or scalar input -> EMPTY  CORRESPONDING OUTPUT (A REVOIR)
    1919%      Other fields in DataIn: copied to DataOut without modification
    20 % Calib: structure containing the calibration parameters (Tsai) or containing a subtree Calib.GeometryCalib with these parameters
     20% Calib: structure containing the  substructure Calib.GeometryCalib with the calibration parameters
    2121%
    2222% call the function  px_XYZ (case of images) for pointwise coordinate transforms
    … …  
    3636    DataOut=px_1(Data,CalibData.GeometryCalib);
    3737end
    38 if exist('Data_1','var')
     38if isfield(DataOut,'Z')
     39    DataOut=rmfield(DataOut,'Z');
     40end
     41if exist('Data_1','var')% if there is a second input field, it is also transformed
    3942    if ~(exist('CalibData_1','var') && isfield(CalibData_1,'GeometryCalib'))
    4043        DataOut_1=Data_1;
    … …  
    4245        DataOut_1=px_1(Data_1,CalibData_1.GeometryCalib);
    4346    end
    44 else
     47    if isfield(DataOut_1,'Z')
     48        DataOut_1=rmfield(DataOut_1,'Z');
     49    end
     50else % no second input field then empty second output field
    4551    DataOut_1=[];
    4652end
    … …  
    5258
    5359%Act only if .CoordType=phys, and Calib defined
    54 if isfield(Data,'CoordType')& isequal(Data.CoordType,'phys')& ~isempty(Calib)
     60if isfield(Data,'CoordType')&& isequal(Data.CoordType,'phys')&& ~isempty(Calib)
    5561    DataOut.CoordType='px'; %put flag for pixel coordinates
    5662    DataOut.CoordUnit='px';
    5763    %transform of X,Y coordinates
    58     if isfield(Data,'Z')&~isempty(Data.Z)
     64    if isfield(Data,'Z')&&~isempty(Data.Z)
    5965        Z=Data.Z;
    6066    else
    … …  
    97103% [Zphys]: corresponding array of z physical coordinates (0 by default)
    98104
    99 
    100 function [X,Y]=px_XYZ(Calib,Xphys,Yphys,Zphys)
    101 X=[];%default
    102 Y=[];
    103 % if exist('Z','var')& isequal(Z,round(Z))& Z>0 & isfield(Calib,'PlanePos')&length(Calib.PlanePos)>=Z
    104 %     Zindex=Z;
    105 %     planepos=Calib.PlanePos{Zindex};
    106 %     zphys=planepos(3);%A GENERALISER CAS AVEC ANGLE
    107 % else
    108 %     zphys=0;
     105%
     106% function [X,Y]=px_XYZ(Calib,Xphys,Yphys,Zphys)
     107% X=[];%default
     108% Y=[];
     109% % if exist('Z','var')& isequal(Z,round(Z))& Z>0 & isfield(Calib,'PlanePos')&length(Calib.PlanePos)>=Z
     110% %     Zindex=Z;
     111% %     planepos=Calib.PlanePos{Zindex};
     112% %     zphys=planepos(3);%A GENERALISER CAS AVEC ANGLE
     113% % else
     114% %     zphys=0;
     115% % end
     116% if ~exist('Zphys','var')
     117%     Zphys=0;
    109118% end
    110 if ~exist('Zphys','var')
    111     Zphys=0;
    112 end
    113 
    114 %%%%%%%%%%%%%
    115 if isfield(Calib,'R')
    116     R=(Calib.R)';
    117     xc=R(1)*Xphys+R(2)*Yphys+R(3)*Zphys+Calib.Tx;
    118     yc=R(4)*Xphys+R(5)*Yphys+R(6)*Zphys+Calib.Ty;
    119     zc=R(7)*Xphys+R(8)*Yphys+R(9)*Zphys+Calib.Tz;
    120 %undistorted image coordinates
    121     Xu=Calib.f*xc./zc;
    122     Yu=Calib.f*yc./zc;
    123 %distorted image coordinates
    124     distortion=(Calib.kappa1)*(Xu.*Xu+Yu.*Yu)+1; %A REVOIR
    125 % distortion=1;
    126     Xd=Xu./distortion;
    127     Yd=Yu./distortion;
    128 %pixel coordinates
    129     X=Xd*Calib.sx/Calib.dpx+Calib.Cx;
    130     Y=Yd/Calib.dpy+Calib.Cy;
    131 
    132 elseif isfield(Calib,'Pxcmx')&isfield(Calib,'Pxcmy')%old calib 
    133         X=Xphys*Calib.Pxcmx;
    134         Y=Yphys*Calib.Pxcmy;
    135 end
     119%
     120% %%%%%%%%%%%%%
     121% if isfield(Calib,'R')
     122%     R=(Calib.R)';
     123%     xc=R(1)*Xphys+R(2)*Yphys+R(3)*Zphys+Calib.Tx;
     124%     yc=R(4)*Xphys+R(5)*Yphys+R(6)*Zphys+Calib.Ty;
     125%     zc=R(7)*Xphys+R(8)*Yphys+R(9)*Zphys+Calib.Tz;
     126% %undistorted image coordinates
     127%     Xu=Calib.f*xc./zc;
     128%     Yu=Calib.f*yc./zc;
     129% %distorted image coordinates
     130%     distortion=(Calib.kappa1)*(Xu.*Xu+Yu.*Yu)+1; %A REVOIR
     131% % distortion=1;
     132%     Xd=Xu./distortion;
     133%     Yd=Yu./distortion;
     134% %pixel coordinates
     135%     X=Xd*Calib.sx/Calib.dpx+Calib.Cx;
     136%     Y=Yd/Calib.dpy+Calib.Cy;
     137%
     138% elseif isfield(Calib,'Pxcmx')&isfield(Calib,'Pxcmy')%old calib 
     139%         X=Xphys*Calib.Pxcmx;
     140%         Y=Yphys*Calib.Pxcmy;
     141% end
    136142
    137143
Note: See TracChangeset for help on using the changeset viewer.