Changeset 114 for trunk


Ignore:
Timestamp:
Oct 14, 2010, 6:35:03 PM (14 years ago)
Author:
sommeria
Message:

geometric calibration procedures updated

Location:
trunk/src
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/geometry_calib.m

    r109 r114  
    4343% Edit the above text to modify the response to help geometry_calib
    4444
    45 % Last Modified by GUIDE v2.5 03-Oct-2010 09:34:12
     45% Last Modified by GUIDE v2.5 05-Oct-2010 13:47:00
    4646
    4747% Begin initialization code - DO NOT edit
    … …  
    9393inputxml='';
    9494if exist('inputfile','var')& ~isempty(inputfile)
    95     struct.XmlInputfile=inputfile;
     95    struct.XmlInputFile=inputfile;
    9696    set(hObject,'UserData',struct)
    9797    [Pathsub,RootFile,field_count,str2,str_a,str_b,ext,nom_type,subdir]=name2display(inputfile);
    … …  
    154154    return
    155155end
    156 
    157 %check error
    158 if isfield(GeometryCalib,'dpx_dpy')
    159     Calib.dpx=GeometryCalib.dpx_dpy(1);
    160     Calib.dpy=GeometryCalib.dpx_dpy(2);
    161 end
    162 if isfield(GeometryCalib,'sx')
    163     Calib.sx=GeometryCalib.sx;
    164 end
    165 if isfield(GeometryCalib,'Cx_Cy')
    166     Calib.Cx=GeometryCalib.Cx_Cy(1);
    167     Calib.Cy=GeometryCalib.Cx_Cy(2);
    168 end
    169 if isfield(GeometryCalib,'kappa1')
    170     Calib.kappa1=GeometryCalib.kappa1;
    171 end
    172 if isfield(GeometryCalib,'focal')
    173     Calib.f=GeometryCalib.focal;
    174 end
    175 if isfield(GeometryCalib,'Tx_Ty_Tz')
    176     Calib.Tx=GeometryCalib.Tx_Ty_Tz(1);
    177     Calib.Ty=GeometryCalib.Tx_Ty_Tz(2);
    178     Calib.Tz=GeometryCalib.Tx_Ty_Tz(3);
    179 end
    180 if isfield(GeometryCalib,'R')
    181     Calib.R=GeometryCalib.R;
    182 end
     156Z_plane=[];
    183157if ~isempty(Coord)
     158    %check error
    184159    X=Coord(:,1);
    185160    Y=Coord(:,2);
    … …  
    187162    x_ima=Coord(:,4);
    188163    y_ima=Coord(:,5);
    189     [Xpoints,Ypoints]=px_XYZ(Calib,X,Y,Z);
     164    [Xpoints,Ypoints]=px_XYZ(GeometryCalib,X,Y,Z);
    190165    GeometryCalib.ErrorRms(1)=sqrt(mean((Xpoints-x_ima).*(Xpoints-x_ima)));
    191166    [GeometryCalib.ErrorMax(1),index(1)]=max(abs(Xpoints-x_ima));
    … …  
    194169    [EM,ind_dim]=max(GeometryCalib.ErrorMax);
    195170    index=index(ind_dim);
     171    %set the Z position of the reference plane used for calibration
    196172    if isequal(max(Z),min(Z))%Z constant
     173        Z_plane=Z(1);
    197174        GeometryCalib.NbSlice=1;
    198         GeometryCalib.SliceCoord=[0 0 Z(1)];
    199     end
    200 end
     175        GeometryCalib.SliceCoord=[0 0 Z_plane];
     176    end
     177end
     178%set the coordinate unit
    201179unitlist=get(handles.CoordUnit,'String');
    202180unit=unitlist{get(handles.CoordUnit,'value')};
    203181GeometryCalib.CoordUnit=unit;
     182%record the points
    204183GeometryCalib.SourceCalib.PointCoord=Coord;
    205 
    206 %display calibration intrinsic parameters
    207 k=0;
    208 if isfield(GeometryCalib,'kappa1')
    209     k=GeometryCalib.kappa1 * GeometryCalib.focal*GeometryCalib.focal;
    210 end
    211 sx=1;dpx_dpy=[1 1];
    212 if isfield(GeometryCalib,'sx')
    213     sx=GeometryCalib.sx;
    214 end
    215 if isfield(GeometryCalib,'dpx_dpy')
    216     dpx_dpy=GeometryCalib.dpx_dpy;
    217 end
    218 Cx_Cy=[0 0];
    219 if isfield(GeometryCalib,'Cx_Cy')
    220     Cx_Cy=GeometryCalib.Cx_Cy;
    221 end
    222 f1=GeometryCalib.focal*sx/dpx_dpy(1);
    223 f2=GeometryCalib.focal/dpx_dpy(2);
    224 Cx=Cx_Cy(1);
    225 Cy=Cx_Cy(2);
    226 set(handles.fx,'String',num2str(f1,4))
    227 set(handles.fy,'String',num2str(f2,4))
    228 set(handles.k,'String',num2str(k,'%1.4f'))
    229 set(handles.Cx,'String',num2str(Cx,'%1.1f'))
    230 set(handles.Cy,'String',num2str(Cy,'%1.1f'))
    231 % Display extrinsinc parameters
     184display_intrinsic(GeometryCalib,handles)%display calibration intrinsic parameters
     185
     186% Display extrinsinc parameters (rotation and translation of camera with  respect to the phys coordiantes)
    232187set(handles.Tx,'String',num2str(GeometryCalib.Tx_Ty_Tz(1),4))
    233188set(handles.Ty,'String',num2str(GeometryCalib.Tx_Ty_Tz(2),4))
    234189set(handles.Tz,'String',num2str(GeometryCalib.Tx_Ty_Tz(3),4))
    235 
    236 % indicate the plane of the calibration grid if defined
    237 huvmat=findobj(allchild(0),'Name','uvmat');
    238 hhuvmat=guidata(huvmat);%handles of elements in the GUI uvmat
     190set(handles.Phi,'String',num2str(GeometryCalib.omc(1),4))
     191set(handles.Theta,'String',num2str(GeometryCalib.omc(2),4))
     192set(handles.Psi,'String',num2str(GeometryCalib.omc(3),4))
     193
     194% store the calibration data, by default in the xml file of the currently displayed image
     195hhuvmat=guidata(findobj(allchild(0),'Name','uvmat'));%handles of elements in the GUI uvmat
    239196RootPath='';
    240197RootFile='';
    … …  
    244201    RootFile=get(hhuvmat.RootFile,'String');
    245202    filebase=fullfile(RootPath,RootFile);
    246     outputfile=[filebase '.xml'];
     203    outputfile=[filebase '.xml'];%xml file associated with the currently displayed image
    247204else
    248205    question={'save the calibration data and point coordinates in'};
    … …  
    256213    ['Error max (along x,y)=' num2str(GeometryCalib.ErrorMax) ' pixels']});
    257214if isequal(answer,'Yes')
    258     update_imadoc(GeometryCalib,outputfile)
     215    answer_1=msgbox_uvmat('INPUT_TXT',' Z= ',num2str(Z_plane));
     216    Z_plane=str2num(answer_1);
     217    GeometryCalib.NbSlice=1;
     218    GeometryCalib.SliceCoord=[0 0 Z_plane];
     219    errormsg=update_imadoc(GeometryCalib,outputfile);% introduce the calibration data in the xml file
     220    if ~strcmp(errormsg,'')
     221        msgbox_uvmat('ERROR',errormsg);
     222    end
    259223    %display image with new calibration in the currently opened uvmat interface
    260224    hhh=findobj(hhuvmat.axes3,'Tag','calib_marker');% delete calib points and markers
    … …  
    268232    set(hhuvmat.FixedLimits,'Value',0)% put FixedLimits option to 'off'
    269233    set(hhuvmat.FixedLimits,'BackgroundColor',[0.7 0.7 0.7])
    270     uvmat('RootPath_Callback',hObject,eventdata,hhuvmat); %file input with xml reading  in uvmat
    271 %     if ip==0
    272         MenuPlot_Callback(hObject, eventdata, handles)
    273         set(handles.ListCoord,'Value',index)% indicate in the list the point with max deviation (possible mistake)
    274         ListCoord_Callback(hObject, eventdata, handles)
    275 %     end
     234    uvmat('RootPath_Callback',hObject,eventdata,hhuvmat); %file input with xml reading  in uvmat, show the image in phys coordinates
     235    MenuPlot_Callback(hObject, eventdata, handles)
     236    set(handles.ListCoord,'Value',index)% indicate in the list the point with max deviation (possible mistake)
     237    ListCoord_Callback(hObject, eventdata, handles)
    276238    figure(handles.geometry_calib)
    277239end
    … …  
    286248Coord_cell=get(handles.ListCoord,'String');
    287249Object=read_geometry_calib(Coord_cell);
    288 GeometryCalib=feval(calib_cell{val},Object.Coord);
    289 
    290 % %record image source
    291 GeometryCalib.SourceCalib.PointCoord=Object.Coord;
     250GeometryCalib=feval(['calib_' calib_cell{val}],Object.Coord,handles);
     251
     252%read the current calibration points
     253Coord_cell=get(handles.ListCoord,'String');
     254Object=read_geometry_calib(Coord_cell);
     255Coord=Object.Coord;
     256
     257% apply the calibration, whose type is selected in  handles.calib_type
     258if ~isempty(Coord)
     259    calib_cell=get(handles.calib_type,'String');
     260    val=get(handles.calib_type,'Value');
     261    GeometryCalib=feval(['calib_' calib_cell{val}],Coord,handles);
     262else
     263    msgbox_uvmat('ERROR','No calibration points, abort')
     264    return
     265end
     266
     267if ~isempty(Coord)
     268    %check error
     269    X=Coord(:,1);
     270    Y=Coord(:,2);
     271    Z=Coord(:,3);
     272    x_ima=Coord(:,4);
     273    y_ima=Coord(:,5);
     274    [Xpoints,Ypoints]=px_XYZ(GeometryCalib,X,Y,Z);
     275    GeometryCalib.ErrorRms(1)=sqrt(mean((Xpoints-x_ima).*(Xpoints-x_ima)));
     276    [GeometryCalib.ErrorMax(1),index(1)]=max(abs(Xpoints-x_ima));
     277    GeometryCalib.ErrorRms(2)=sqrt(mean((Ypoints-y_ima).*(Ypoints-y_ima)));
     278    [GeometryCalib.ErrorMax(2),index(2)]=max(abs(Ypoints-y_ima));
     279    [EM,ind_dim]=max(GeometryCalib.ErrorMax);
     280    index=index(ind_dim);
     281    %set the Z position of the reference plane used for calibration
     282    Z_plane=[];
     283    if isequal(max(Z),min(Z))
     284        Z_plane=Z(1);
     285    end
     286    answer_1=msgbox_uvmat('INPUT_TXT',' Z= ',num2str(Z_plane));
     287    Z_plane=str2num(answer_1);
     288    GeometryCalib.NbSlice=1;
     289    GeometryCalib.SliceCoord=[0 0 Z_plane];
     290    %set the coordinate unit
     291    unitlist=get(handles.CoordUnit,'String');
     292    unit=unitlist{get(handles.CoordUnit,'value')};
     293    GeometryCalib.CoordUnit=unit;
     294    %record the points
     295    GeometryCalib.SourceCalib.PointCoord=Coord;
     296    errormsg=update_imadoc(GeometryCalib,outputfile);% introduce the calibration data in the xml file
     297    if ~strcmp(errormsg,'')
     298        msgbox_uvmat('ERROR',errormsg);
     299    end
     300end
     301display_intrinsic(GeometryCalib,handles)%display calibration intrinsic parameters
     302
     303% Display extrinsinc parameters (rotation and translation of camera with  respect to the phys coordiantes)
     304set(handles.Tx,'String',num2str(GeometryCalib.Tx_Ty_Tz(1),4))
     305set(handles.Ty,'String',num2str(GeometryCalib.Tx_Ty_Tz(2),4))
     306set(handles.Tz,'String',num2str(GeometryCalib.Tx_Ty_Tz(3),4))
     307set(handles.Phi,'String',num2str(GeometryCalib.omc(1),4))
     308set(handles.Theta,'String',num2str(GeometryCalib.omc(2),4))
     309set(handles.Psi,'String',num2str(GeometryCalib.omc(3),4))
     310
     311% indicate the plane of the calibration grid if defined
     312% huvmat=findobj(allchild(0),'Name','uvmat');
     313% hhuvmat=guidata(huvmat);%handles of elements in the GUI uvmat
     314% RootPath='';
     315% RootFile='';
     316% if ~isempty(hhuvmat.RootPath)& ~isempty(hhuvmat.RootFile)
     317%     testhandle=1;
     318%     RootPath=get(hhuvmat.RootPath,'String');
     319%     RootFile=get(hhuvmat.RootFile,'String');
     320%     filebase=fullfile(RootPath,RootFile);
     321%     outputfile=[filebase '.xml'];
     322% else
     323%     question={'save the calibration data and point coordinates in'};
     324%     def={fullfile(RootPath,['ObjectCalib.xml'])};
     325%     options.Resize='on';
     326%     answer=inputdlg(question,'save average in a new file',1,def,options);
     327%     outputfile=answer{1};
     328% end
    292329
    293330%open and read the dataview GUI
    … …  
    298335CalibData=get(handles.geometry_calib,'UserData');%read the calibration image source on the interface userdata
    299336
    300 if isfield(CalibData,'XmlInput')
    301     XmlInput=fileparts(CalibData.XmlInput);
     337if isfield(CalibData,'XmlInputFile')
     338    XmlInput=fileparts(CalibData.XmlInputFile);
    302339    [XmlInput,filename,ext]=fileparts(XmlInput);
    303340end
    … …  
    326363if ~testinput
    327364    filename='PROJETS';%default
    328     if isfield(CalibData,'XmlInput')
    329          [pp,filename]=fileparts(CalibData.XmlInput);
     365    if isfield(CalibData,'XmlInputFile')
     366         [pp,filename]=fileparts(CalibData.XmlInputFile);
    330367    end
    331368    while ~isequal(filename,'PROJETS') && numel(filename)>1
    … …  
    354391T_y=py(2);
    355392GeometryCalib.CalibrationType='rescale';
    356 GeometryCalib.focal=1;%default
    357 GeometryCalib.sx=px(1)/py(1);
     393GeometryCalib.fx_fy=[px(1) py(1)];
    358394GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function
    359 GeometryCalib.Tx_Ty_Tz=[T_x T_y GeometryCalib.focal/py(1)];
    360 GeometryCalib.R=[1,0,0;0,1,0;0,0,0];
     395GeometryCalib.Tx_Ty_Tz=[px(2)/px(1) py(2)/py(1) 1];
     396%GeometryCalib.R=[1,0,0;0,1,0;0,0,0];
     397GeometryCalib.omc=[0 0 0];
    361398
    362399%------------------------------------------------------------------------
    363400% determine the parameters for a calibration by a linear transform matrix (rescale and rotation)
    364 function GeometryCalib=calib_linear(Coord,handles)
     401function GeometryCalib=calib_linear(Coord,handles) %TO UPDATE
    365402%------------------------------------------------------------------------
    366403X=Coord(:,1);
    … …  
    385422%------------------------------------------------------------------------
    386423% determine the tsai parameters for a view normal to the grid plane
     424% NOT USED
    387425function GeometryCalib=calib_normal(Coord,handles)
    388426%------------------------------------------------------------------------
    389427Calib.f1=str2num(get(handles.fx,'String'));
    390428Calib.f2=str2num(get(handles.fy,'String'));
    391 Calib.k=str2num(get(handles.k,'String'));
     429Calib.k=str2num(get(handles.kc,'String'));
    392430Calib.Cx=str2num(get(handles.Cx,'String'));
    393431Calib.Cy=str2num(get(handles.Cy,'String'));
    … …  
    449487%------------------------------------------------------------------------
    450488function GeometryCalib=calib_3D_linear(Coord,handles)
     489%TO UPDATE
    451490%------------------------------------------------------------------
    452491path_uvmat=which('uvmat');% check the path detected for source file uvmat
    … …  
    497536
    498537%retrieve the calibration points stored in the files listed in the popup list coord_files
    499 x_1=Coord(:,4:5)';
    500 X_1=Coord(:,1:3)';
     538x_1=Coord(:,4:5)';%px coordinates of the ref points
     539nx=str2num(get(hhuvmat.npx,'String'));
     540ny=str2num(get(hhuvmat.npy,'String'));
     541x_1(2,:)=ny-x_1(2,:);%reverse the y image coordinates
     542X_1=Coord(:,1:3)';%phys coordinates of the ref points
    501543n_ima=numel(coord_files)+1;
    502544if ~isempty(coord_files)
    … …  
    516558                end
    517559                eval(['x_' num2str(ifile+1) '=Coord_file(:,4:5)'';']);
     560                eval(['x_' num2str(ifile+1) '(2,:)=ny-x_' num2str(ifile+1) '(2,:);' ]);
    518561                eval(['X_' num2str(ifile+1) '=Coord_file(:,1:3)'';']);
    519562                end
    … …  
    525568n_ima=numel(coord_files)+1;
    526569whos
    527 nx=str2num(get(hhuvmat.npx,'String'));
    528 ny=str2num(get(hhuvmat.npy,'String'));
    529570
    530571est_dist=[1;0;0;0;0];
    … …  
    536577
    537578GeometryCalib.CalibrationType='3D_quadr';
    538 GeometryCalib.focal=fc(2);
    539 GeometryCalib.dpx_dpy=[1 1];
     579GeometryCalib.fx_fy=fc';
     580%GeometryCalib.focal=fc(2);
     581%GeometryCalib.dpx_dpy=[1 1];
    540582GeometryCalib.Cx_Cy=cc';
    541 GeometryCalib.sx=fc(1)/fc(2);
    542 GeometryCalib.kappa1=-kc(1)/fc(2)^2;
     583%GeometryCalib.sx=fc(1)/fc(2);
     584GeometryCalib.kc=kc(1);
     585%GeometryCalib.kappa1=-kc(1)/fc(2)^2;
    543586GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function
    544587GeometryCalib.Tx_Ty_Tz=Tc_1';
    545588GeometryCalib.R=Rc_1;
    546 GeometryCalib.R(3,1:3)=-GeometryCalib.R(3,1:3);%inversion for z upward
     589GeometryCalib.R(2,1:3)=-GeometryCalib.R(2,1:3);%inversion of the y image coordinate
     590GeometryCalib.Tx_Ty_Tz(2)=-GeometryCalib.Tx_Ty_Tz(2);%inversion of the y image coordinate
     591GeometryCalib.Cx_Cy(2)=ny-GeometryCalib.Cx_Cy(2);%inversion of the y image coordinate
     592GeometryCalib.omc=(180/pi)*omc_1;%angles in degrees
    547593GeometryCalib.ErrorRMS=[];
    548594GeometryCalib.ErrorMax=[];
    … …  
    556602x_1=Coord(:,4:5)';
    557603X_1=Coord(:,1:3)';
     604huvmat=findobj(allchild(0),'Tag','uvmat');
     605hhuvmat=guidata(huvmat);
     606ny=str2num(get(hhuvmat.npy,'String'));
     607x_1(2,:)=ny-x_1(2,:);%reverse the y image coordinates
    558608n_ima=1;
    559 Calib.f1=str2num(get(handles.fx,'String'));
    560 Calib.f2=str2num(get(handles.fy,'String'));
    561 Calib.k=str2num(get(handles.k,'String'));
    562 Calib.Cx=str2num(get(handles.Cx,'String'));
    563 Calib.Cy=str2num(get(handles.Cy,'String'));
    564 
     609GeometryCalib.CalibrationType='3D_extrinsic';
     610GeometryCalib.fx_fy(1)=str2num(get(handles.fx,'String'));
     611GeometryCalib.fx_fy(2)=str2num(get(handles.fy,'String'));
     612GeometryCalib.Cx_Cy(1)=str2num(get(handles.Cx,'String'));
     613GeometryCalib.Cx_Cy(2)=str2num(get(handles.Cy,'String'));
     614GeometryCalib.kc=str2num(get(handles.kc,'String'));
    565615fct_path=fullfile(path_UVMAT,'toolbox_calib');
    566616addpath(fct_path)
     617Cx_Cy=ny-(GeometryCalib.Cx_Cy)';%reverse Cx_Cy(2) for calibration (inversion of px ordinate)
    567618% [omc1,Tc1,Rc1,H,x,ex,JJ] = compute_extrinsic(x_1,X_1,...
    568619%     [Calib.f Calib.f*Calib.sx]',...
    569620%     [Calib.Cx Calib.Cy]',...
    570621%     [-Calib.kappa1*Calib.f^2 0 0 0 0]);
    571 [omc1,Tc1,Rc1,H,x,ex,JJ] = compute_extrinsic(x_1,X_1,...
    572     [Calib.f1 Calib.f2]',...
    573     [Calib.Cx Calib.Cy]',...
    574     [-Calib.k 0 0 0 0]);
     622[omc,Tc1,Rc1,H,x,ex,JJ] = compute_extrinsic(x_1,X_1,...
     623    (GeometryCalib.fx_fy)',...
     624    Cx_Cy,[GeometryCalib.kc 0 0 0 0]);
    575625%get the euler angles ???
    576626rmpath(fct_path);
    577627
    578628std(ex')
    579 GeometryCalib.CalibrationType='3D_extrinsic';
    580 GeometryCalib.focal=Calib.f2;
    581 GeometryCalib.dpx_dpy=[1 1];
    582 GeometryCalib.Cx_Cy=[Calib.Cx Calib.Cy];
    583 GeometryCalib.sx=Calib.f1/Calib.f2;
    584 GeometryCalib.kappa1=Calib.k/(Calib.f2*Calib.f2);
    585629GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function
    586630GeometryCalib.Tx_Ty_Tz=Tc1';
    587631%inversion of z axis
    588632GeometryCalib.R=Rc1;
     633GeometryCalib.R(2,1:3)=-GeometryCalib.R(2,1:3);%inversion of the y image coordinate
     634GeometryCalib.Tx_Ty_Tz(2)=-GeometryCalib.Tx_Ty_Tz(2);%inversion of the y image coordinate
     635%GeometryCalib.Cx_Cy(2)=ny-GeometryCalib.Cx_Cy(2);%inversion of the y image coordinate
     636GeometryCalib.omc=(180/pi)*omc';
    589637%GeometryCalib.R(3,1:3)=-GeometryCalib.R(3,1:3);%inversion for z upward
    590638
    … …  
    608656%--------------------------------------------------------------------------
    609657function GeometryCalib=calib_tsai(Coord,handles)% old version using gauthier's bianry ccal_fo
     658% NOT USED
    610659%------------------------------------------------------------------------
    611660%TSAI
    … …  
    728777    end
    729778    outputfile=[filebase '.xml'];
    730     update_imadoc(GeometryCalib,outputfile)
     779    errormsg=update_imadoc(GeometryCalib,outputfile);
     780    if ~strcmp(errormsg,'')
     781        msgbox_uvmat('ERROR',errormsg);
     782    end
    731783    listfile=get(handles.coord_files,'string');
    732784    if isequal(listfile,{''})
    … …  
    739791set(handles.ListCoord,'Value',1)% refresh the display of coordinates
    740792set(handles.ListCoord,'String',{'......'})
     793
     794% --------------------------------------------------------------------
     795% --- Executes on button press in CLEAR_PTS: clear the list of calibration points
     796function CLEAR_PTS_Callback(hObject, eventdata, handles)
     797% --------------------------------------------------------------------
     798set(handles.ListCoord,'Value',1)% refresh the display of coordinates
     799set(handles.ListCoord,'String',{'......'})
     800MenuPlot_Callback(hObject, eventdata, handles)
    741801
    742802%------------------------------------------------------------------------
    … …  
    892952
    893953%------------------------------------------------------------------------
    894 % --- 'key_press_fcn:' function activated when a key is pressed on the keyboard
    895 function key_press_fcn(hObject,eventdata,handles)
    896 %------------------------------------------------------------------------
    897 xx=double(get(handles.geometry_calib,'CurrentCharacter')); %get the keyboard character
    898 if ismember(xx,[8 127])%backspace or delete
    899     Coord_cell=get(handles.ListCoord,'String');
    900     val=get(handles.ListCoord,'Value');
    901      if max(val)<numel(Coord_cell) % the last element '...' has not been selected
    902         Coord_cell(val)=[];%remove the selected line
    903         set(handles.ListCoord,'Value',min(val))
    904         set(handles.ListCoord,'String',Coord_cell)         
    905         ListCoord_Callback(hObject, eventdata, handles)
    906         MenuPlot_Callback(hObject,eventdata,handles)
    907      end
    908 end
    909 
    910 %------------------------------------------------------------------------
    911954function MenuPlot_Callback(hObject, eventdata, handles)
    912955%------------------------------------------------------------------------
    … …  
    10501093C = [l(7:8)' 1];
    10511094
    1052 GeometryCalib.CalibrationType='tsai';
    1053 GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function
    1054 GeometryCalib.f=1;
    1055 GeometryCalib.dpx=1;
    1056 GeometryCalib.dpy=1;
    1057 GeometryCalib.sx=1;
    1058 GeometryCalib.Cx=0;
    1059 GeometryCalib.Cy=0;
    1060 GeometryCalib.kappa1=0;
    1061 GeometryCalib.Tx=Amat(1,3);
    1062 GeometryCalib.Ty=Amat(2,3);
    1063 GeometryCalib.Tz=1;
     1095%GeometryCalib.CalibrationType='tsai';
     1096%GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function
     1097% GeometryCalib.f=1;
     1098% GeometryCalib.dpx=1;
     1099% GeometryCalib.dpy=1;
     1100% GeometryCalib.sx=1;
     1101% GeometryCalib.Cx=0;
     1102% GeometryCalib.Cy=0;
     1103% GeometryCalib.kappa1=0;
     1104% GeometryCalib.Tx=Amat(1,3);
     1105% GeometryCalib.Ty=Amat(2,3);
     1106% GeometryCalib.Tz=1;
     1107% GeometryCalib.R=[Amat(1,1),Amat(1,2),0;Amat(2,1),Amat(2,2),0;C(1),C(2),0];
     1108%
     1109% [Amod,Rangx,Rangy]=phys_Ima(A-min(min(A)),GeometryCalib,0);
     1110
     1111GeometryCalib.fx_fy=[1 1];
     1112GeometryCalib.Tx_Ty_Tz=[Amat(1,3) Amat(2,3) 1];
    10641113GeometryCalib.R=[Amat(1,1),Amat(1,2),0;Amat(2,1),Amat(2,2),0;C(1),C(2),0];
    1065 
    1066 [Amod,Rangx,Rangy]=phys_Ima(A-min(min(A)),GeometryCalib,0);
     1114path_uvmat=which('uvmat');% check the path detected for source file uvmat
     1115path_UVMAT=fileparts(path_uvmat); %path to UVMAT
     1116addpath(fullfile(path_UVMAT,'transform_field'))
     1117Data.ListVarName={'AY','AX','A'};
     1118Data.VarDimName={'AY','AX',{'AY','AX'}};
     1119if ndims(A)==3
     1120    A=mean(A,3);
     1121end
     1122Data.A=A-min(min(A));
     1123Data.AY=[npxy(1)-0.5 0.5];
     1124Data.AX=[0.5 npxy(2)];
     1125Data.CoordType='px';
     1126Calib.GeometryCalib=GeometryCalib;
     1127DataOut=phys(Data,Calib)
     1128rmpath(fullfile(path_UVMAT,'transform_field'))
     1129Amod=DataOut.A;
     1130Rangx=DataOut.AX;
     1131Rangy=DataOut.AY;
     1132% GeometryCalib.dpx=1;
     1133% GeometryCalib.dpy=1;
     1134% GeometryCalib.sx=1;
     1135% GeometryCalib.Cx=0;
     1136% GeometryCalib.Cy=0;
     1137% GeometryCalib.kappa1=0;
     1138% GeometryCalib.Tx=Amat(1,3);
     1139% GeometryCalib.Ty=Amat(2,3);
     1140% GeometryCalib.Tz=1;
     1141% GeometryCalib.R=[Amat(1,1),Amat(1,2),0;Amat(2,1),Amat(2,2),0;C(1),C(2),0];
     1142%
     1143% [Amod,Rangx,Rangy]=phys_Ima(A-min(min(A)),GeometryCalib,0);
     1144
    10671145
    10681146Amod=double(Amod);
    1069 % figure(12) display corrected image
     1147% figure(12) %display corrected image
    10701148% Amax=max(max(Amod));
    10711149% image(Rangx,Rangy,uint8(255*Amod/Amax))
     1150
    10721151Dx=(Rangx(2)-Rangx(1))/(npxy(2)-1); %x mesh in real space
    10731152Dy=(Rangy(2)-Rangy(1))/(npxy(1)-1); %y mesh in real space
    1074 ind_range_x=ceil(GeometryCalib.R(1,1)*CalibData.grid.Dx/3);% range of search of image ma around each point obtained by linear interpolation from the marked points
    1075 ind_range_y=ceil(GeometryCalib.R(2,2)*CalibData.grid.Dy/3);% range of search of image ma around each point obtained by linear interpolation from the marked points
     1153ind_range_x=ceil(abs(GeometryCalib.R(1,1)*CalibData.grid.Dx/3))% range of search of image ma around each point obtained by linear interpolation from the marked points
     1154ind_range_y=ceil(abs(GeometryCalib.R(2,2)*CalibData.grid.Dy/3))% range of search of image ma around each point obtained by linear interpolation from the marked points
    10761155nbpoints=size(T,1);
    10771156for ipoint=1:nbpoints
    … …  
    10951174        x_shift=sum(Atop.*[-2 -1 0 1 2])/sum(Atop);
    10961175    end
    1097     if ind_x_max+2<=numel(y_profile) && ind_y_max-2>=1
     1176    if ind_y_max+2<=numel(y_profile) && ind_y_max-2>=1
    10981177        Atop=y_profile(ind_y_max-2:ind_y_max+2);
    10991178        y_shift=sum(Atop.*[-2 -1 0 1 2]')/sum(Atop);
    … …  
    12081287%INPUT:
    12091288%Zindex: index of plane
    1210 function [A_out,Rangx,Rangy]=phys_Ima(A,Calib,ZIndex)
    1211 %------------------------------------------------------------------------
    1212 xcorner=[];
    1213 ycorner=[];
    1214 npx=[];
    1215 npy=[];
    1216 siz=size(A)
    1217 npx=[npx siz(2)];
    1218 npy=[npy siz(1)]
    1219 xima=[0.5 siz(2)-0.5 0.5 siz(2)-0.5];%image coordinates of corners
    1220 yima=[0.5 0.5 siz(1)-0.5 siz(1)-0.5];
    1221 [xcorner,ycorner]=phys_XYZ(Calib,xima,yima,ZIndex);%corresponding physical coordinates
    1222 Rangx(1)=min(xcorner);
    1223 Rangx(2)=max(xcorner);
    1224 Rangy(2)=min(ycorner);
    1225 Rangy(1)=max(ycorner);
    1226 test_multi=(max(npx)~=min(npx)) | (max(npy)~=min(npy));
    1227 npx=max(npx);
    1228 npy=max(npy);
    1229 x=linspace(Rangx(1),Rangx(2),npx);
    1230 y=linspace(Rangy(1),Rangy(2),npy);
    1231 [X,Y]=meshgrid(x,y);%grid in physical coordiantes
    1232 vec_B=[];
    1233 
    1234 zphys=0; %default
    1235 if isfield(Calib,'SliceCoord') %.Z= index of plane
    1236    SliceCoord=Calib.SliceCoord(ZIndex,:);
    1237    zphys=SliceCoord(3); %to generalize for non-parallel planes
    1238 end
    1239 [XIMA,YIMA]=px_XYZ(Calib,X,Y,zphys);%corresponding image indices for each point in the real space grid
    1240 XIMA=reshape(round(XIMA),1,npx*npy);%indices reorganized in 'line'
    1241 YIMA=reshape(round(YIMA),1,npx*npy);
    1242 flagin=XIMA>=1 & XIMA<=npx & YIMA >=1 & YIMA<=npy;%flagin=1 inside the original image
    1243 testuint8=isa(A,'uint8');
    1244 testuint16=isa(A,'uint16');
    1245 if numel(siz)==2 %(B/W images)
    1246     vec_A=reshape(A,1,npx*npy);%put the original image in line
    1247     ind_in=find(flagin);
    1248     ind_out=find(~flagin);
    1249     ICOMB=((XIMA-1)*npy+(npy+1-YIMA));
    1250     ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A
    1251     vec_B(ind_in)=vec_A(ICOMB);
    1252     vec_B(ind_out)=zeros(size(ind_out));
    1253     A_out=reshape(vec_B,npy,npx);%new image in real coordinates
    1254 elseif numel(siz)==3     
    1255     for icolor=1:siz(3)
    1256         vec_A=reshape(A{icell}(:,:,icolor),1,npx*npy);%put the original image in line
    1257         ind_in=find(flagin);
    1258         ind_out=find(~flagin);
    1259         ICOMB=((XIMA-1)*npy+(npy+1-YIMA));
    1260         ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A
    1261         vec_B(ind_in)=vec_A(ICOMB);
    1262         vec_B(ind_out)=zeros(size(ind_out));
    1263         A_out(:,:,icolor)=reshape(vec_B,npy,npx);%new image in real coordinates
    1264     end
    1265 end
    1266 if testuint8
    1267     A_out=uint8(A_out);
    1268 end
    1269 if testuint16
    1270     A_out=uint16(A_out);
    1271 end
     1289% function [A_out,Rangx,Rangy]=phys_Ima(A,Calib,ZIndex)
     1290% %------------------------------------------------------------------------
     1291% xcorner=[];
     1292% ycorner=[];
     1293% npx=[];
     1294% npy=[];
     1295% siz=size(A)
     1296% npx=[npx siz(2)];
     1297% npy=[npy siz(1)]
     1298% xima=[0.5 siz(2)-0.5 0.5 siz(2)-0.5];%image coordinates of corners
     1299% yima=[0.5 0.5 siz(1)-0.5 siz(1)-0.5];
     1300% [xcorner,ycorner]=phys_XYZ(Calib,xima,yima,ZIndex);%corresponding physical coordinates
     1301% Rangx(1)=min(xcorner);
     1302% Rangx(2)=max(xcorner);
     1303% Rangy(2)=min(ycorner);
     1304% Rangy(1)=max(ycorner);
     1305% test_multi=(max(npx)~=min(npx)) | (max(npy)~=min(npy));
     1306% npx=max(npx);
     1307% npy=max(npy);
     1308% x=linspace(Rangx(1),Rangx(2),npx);
     1309% y=linspace(Rangy(1),Rangy(2),npy);
     1310% [X,Y]=meshgrid(x,y);%grid in physical coordiantes
     1311% vec_B=[];
     1312%
     1313% zphys=0; %default
     1314% if isfield(Calib,'SliceCoord') %.Z= index of plane
     1315%    SliceCoord=Calib.SliceCoord(ZIndex,:);
     1316%    zphys=SliceCoord(3); %to generalize for non-parallel planes
     1317% end
     1318% [XIMA,YIMA]=px_XYZ(Calib,X,Y,zphys);%corresponding image indices for each point in the real space grid
     1319% XIMA=reshape(round(XIMA),1,npx*npy);%indices reorganized in 'line'
     1320% YIMA=reshape(round(YIMA),1,npx*npy);
     1321% flagin=XIMA>=1 & XIMA<=npx & YIMA >=1 & YIMA<=npy;%flagin=1 inside the original image
     1322% testuint8=isa(A,'uint8');
     1323% testuint16=isa(A,'uint16');
     1324% if numel(siz)==2 %(B/W images)
     1325%     vec_A=reshape(A,1,npx*npy);%put the original image in line
     1326%     ind_in=find(flagin);
     1327%     ind_out=find(~flagin);
     1328%     ICOMB=((XIMA-1)*npy+(npy+1-YIMA));
     1329%     ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A
     1330%     vec_B(ind_in)=vec_A(ICOMB);
     1331%     vec_B(ind_out)=zeros(size(ind_out));
     1332%     A_out=reshape(vec_B,npy,npx);%new image in real coordinates
     1333% elseif numel(siz)==3     
     1334%     for icolor=1:siz(3)
     1335%         vec_A=reshape(A{icell}(:,:,icolor),1,npx*npy);%put the original image in line
     1336%         ind_in=find(flagin);
     1337%         ind_out=find(~flagin);
     1338%         ICOMB=((XIMA-1)*npy+(npy+1-YIMA));
     1339%         ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A
     1340%         vec_B(ind_in)=vec_A(ICOMB);
     1341%         vec_B(ind_out)=zeros(size(ind_out));
     1342%         A_out(:,:,icolor)=reshape(vec_B,npy,npx);%new image in real coordinates
     1343%     end
     1344% end
     1345% if testuint8
     1346%     A_out=uint8(A_out);
     1347% end
     1348% if testuint16
     1349%     A_out=uint16(A_out);
     1350% end
    12721351
    12731352%------------------------------------------------------------------------
    … …  
    12751354%INPUT:
    12761355%Z: index of plane
    1277 function [Xphys,Yphys,Zphys]=phys_XYZ(Calib,X,Y,Z)
    1278 %------------------------------------------------------------------------
    1279 if exist('Z','var')& isequal(Z,round(Z))& Z>0 & isfield(Calib,'SliceCoord')&length(Calib.SliceCoord)>=Z
    1280     Zindex=Z;
    1281     Zphys=Calib.SliceCoord(Zindex,3);%GENERALISER AUX CAS AVEC ANGLE
    1282 else
    1283     Zphys=0;
    1284 end
    1285 if ~exist('X','var')||~exist('Y','var')
    1286     Xphys=[];
    1287     Yphys=[];%default
    1288     return
    1289 end
    1290 Xphys=X;%default
    1291 Yphys=Y;
    1292 %image transform
    1293 if isfield(Calib,'R')
    1294     R=(Calib.R)';
    1295     Dx=R(5)*R(7)-R(4)*R(8);
    1296     Dy=R(1)*R(8)-R(2)*R(7);
    1297     D0=Calib.f*(R(2)*R(4)-R(1)*R(5));
    1298     Z11=R(6)*R(8)-R(5)*R(9);
    1299     Z12=R(2)*R(9)-R(3)*R(8); 
    1300     Z21=R(4)*R(9)-R(6)*R(7);
    1301     Z22=R(3)*R(7)-R(1)*R(9);
    1302     Zx0=R(3)*R(5)-R(2)*R(6);
    1303     Zy0=R(1)*R(6)-R(3)*R(4);
    1304     A11=R(8)*Calib.Ty-R(5)*Calib.Tz+Z11*Zphys;
    1305     A12=R(2)*Calib.Tz-R(8)*Calib.Tx+Z12*Zphys;
    1306     A21=-R(7)*Calib.Ty+R(4)*Calib.Tz+Z21*Zphys;
    1307     A22=-R(1)*Calib.Tz+R(7)*Calib.Tx+Z11*Zphys;
    1308     X0=Calib.f*(R(5)*Calib.Tx-R(2)*Calib.Ty+Zx0*Zphys);
    1309     Y0=Calib.f*(-R(4)*Calib.Tx+R(1)*Calib.Ty+Zy0*Zphys);
    1310         %px to camera:
    1311     Xd=(Calib.dpx/Calib.sx)*(X-Calib.Cx); % sensor coordinates
    1312     Yd=Calib.dpy*(Y-Calib.Cy);
    1313     dist_fact=1+Calib.kappa1*(Xd.*Xd+Yd.*Yd); %distortion factor
    1314     Xu=dist_fact.*Xd;%undistorted sensor coordinates
    1315     Yu=dist_fact.*Yd;
    1316     denom=Dx*Xu+Dy*Yu+D0;
    1317     % denom2=denom.*denom;
    1318     Xphys=(A11.*Xu+A12.*Yu+X0)./denom;%world coordinates
    1319     Yphys=(A21.*Xu+A22.*Yu+Y0)./denom;
    1320 end
     1356% function [Xphys,Yphys,Zphys]=phys_XYZ(Calib,X,Y,Z)
     1357% %------------------------------------------------------------------------
     1358% if exist('Z','var')& isequal(Z,round(Z))& Z>0 & isfield(Calib,'SliceCoord')&length(Calib.SliceCoord)>=Z
     1359%     Zindex=Z;
     1360%     Zphys=Calib.SliceCoord(Zindex,3);%GENERALISER AUX CAS AVEC ANGLE
     1361% else
     1362%     Zphys=0;
     1363% end
     1364% if ~exist('X','var')||~exist('Y','var')
     1365%     Xphys=[];
     1366%     Yphys=[];%default
     1367%     return
     1368% end
     1369% Xphys=X;%default
     1370% Yphys=Y;
     1371% %image transform
     1372% if isfield(Calib,'R')
     1373%     R=(Calib.R)';
     1374%     Dx=R(5)*R(7)-R(4)*R(8);
     1375%     Dy=R(1)*R(8)-R(2)*R(7);
     1376%     D0=Calib.f*(R(2)*R(4)-R(1)*R(5));
     1377%     Z11=R(6)*R(8)-R(5)*R(9);
     1378%     Z12=R(2)*R(9)-R(3)*R(8); 
     1379%     Z21=R(4)*R(9)-R(6)*R(7);
     1380%     Z22=R(3)*R(7)-R(1)*R(9);
     1381%     Zx0=R(3)*R(5)-R(2)*R(6);
     1382%     Zy0=R(1)*R(6)-R(3)*R(4);
     1383%     A11=R(8)*Calib.Ty-R(5)*Calib.Tz+Z11*Zphys;
     1384%     A12=R(2)*Calib.Tz-R(8)*Calib.Tx+Z12*Zphys;
     1385%     A21=-R(7)*Calib.Ty+R(4)*Calib.Tz+Z21*Zphys;
     1386%     A22=-R(1)*Calib.Tz+R(7)*Calib.Tx+Z11*Zphys;
     1387%     X0=Calib.f*(R(5)*Calib.Tx-R(2)*Calib.Ty+Zx0*Zphys);
     1388%     Y0=Calib.f*(-R(4)*Calib.Tx+R(1)*Calib.Ty+Zy0*Zphys);
     1389%         %px to camera:
     1390%     Xd=(Calib.dpx/Calib.sx)*(X-Calib.Cx); % sensor coordinates
     1391%     Yd=Calib.dpy*(Y-Calib.Cy);
     1392%     dist_fact=1+Calib.kappa1*(Xd.*Xd+Yd.*Yd); %distortion factor
     1393%     Xu=dist_fact.*Xd;%undistorted sensor coordinates
     1394%     Yu=dist_fact.*Yd;
     1395%     denom=Dx*Xu+Dy*Yu+D0;
     1396%     % denom2=denom.*denom;
     1397%     Xphys=(A11.*Xu+A12.*Yu+X0)./denom;%world coordinates
     1398%     Yphys=(A21.*Xu+A22.*Yu+Y0)./denom;
     1399% end
    13211400
    13221401
    … …  
    13291408[s,errormsg]=imadoc2struct(fileinput,'GeometryCalib');
    13301409GeometryCalib=s.GeometryCalib;
    1331 GeometryCalib=load_calib(hObject, eventdata, handles)
     1410%GeometryCalib=load_calib(hObject, eventdata, handles)
    13321411calib=reshape(GeometryCalib.PointCoord',[],1);
    13331412for ilist=1:numel(calib)
    … …  
    13501429[s,errormsg]=imadoc2struct(fileinput,'GeometryCalib');
    13511430GeometryCalib=s.GeometryCalib;
    1352 k=GeometryCalib.kappa1 * GeometryCalib.f*GeometryCalib.f;
    1353 f1=GeometryCalib.f*GeometryCalib.sx/GeometryCalib.dpx;
    1354 f2=GeometryCalib.f/GeometryCalib.dpy;
    1355 Cx=GeometryCalib.Cx;
    1356 Cy=GeometryCalib.Cy;
    1357 set(handles.fx,'String',num2str(f1,'%1.1f'))
    1358 set(handles.fy,'String',num2str(f2,'%1.1f'))
    1359 set(handles.k,'String',num2str(k,'%1.4f'))
    1360 set(handles.Cx,'String',num2str(Cx,'%1.1f'))
    1361 set(handles.Cy,'String',num2str(Cy,'%1.1f'))
     1431display_intrinsic(GeometryCalib,handles)
    13621432
    13631433% -----------------------------------------------------------------------
    … …  
    13831453
    13841454%------------------------------------------------------------------------
     1455% --- 'key_press_fcn:' function activated when a key is pressed on the keyboard
     1456function key_press_fcn(hObject,eventdata,handles)
     1457%------------------------------------------------------------------------
     1458xx=double(get(handles.geometry_calib,'CurrentCharacter')); %get the keyboard character
     1459if ismember(xx,[8 127])%backspace or delete
     1460    Coord_cell=get(handles.ListCoord,'String');
     1461    val=get(handles.ListCoord,'Value');
     1462     if max(val)<numel(Coord_cell) % the last element '...' has not been selected
     1463        Coord_cell(val)=[];%remove the selected line
     1464        set(handles.ListCoord,'Value',min(val))
     1465        set(handles.ListCoord,'String',Coord_cell)         
     1466        ListCoord_Callback(hObject, eventdata, handles)
     1467        MenuPlot_Callback(hObject,eventdata,handles)
     1468     end
     1469end
     1470
     1471%------------------------------------------------------------------------
    13851472function fileinput=browse_xml(hObject, eventdata, handles)
    13861473%------------------------------------------------------------------------
    13871474fileinput=[];%default
    13881475oldfile=''; %default
    1389 UserData=get(handles.geometry_calib,'UserData');
    1390 if isfield(UserData,'XmlInput')
    1391     oldfile=UserData.XmlInput;
     1476UserData=get(handles.geometry_calib,'UserData')
     1477if isfield(UserData,'XmlInputFile')
     1478    oldfile=UserData.XmlInputFile;
    13921479end
    13931480[FileName, PathName, filterindex] = uigetfile( ...
    … …  
    14041491sizf=size(fileinput);
    14051492if (~ischar(fileinput)||~isequal(sizf(1),1)),return;end
    1406 UserData.XmlInput=fileinput;
     1493UserData.XmlInputFile=fileinput;
    14071494set(handles.geometry_calib,'UserData',UserData)%record current file foer further use of browser
    14081495
    … …  
    14111498%------------------------------------------------------------------------
    14121499[s,errormsg]=imadoc2struct(fileinput,'GeometryCalib');
    1413 GeometryCalib=s.GeometryCalib;
    1414 if isempty(GeometryCalib)
    1415     Tabchar={};
    1416     CoordCell={};
    1417     k=0;%default
    1418     f1=1000;
    1419     f2=1000;
    1420     hhuvmat=guidata(findobj(allchild(0),'Name','uvmat'));
    1421     Cx=str2num(get(hhuvmat.npx,'String'))/2;
    1422     Cy=str2num(get(hhuvmat.npy,'String'))/2;
    1423 else
    1424     k=GeometryCalib.kappa1 * GeometryCalib.f*GeometryCalib.f;
    1425     f1=GeometryCalib.f*GeometryCalib.sx/GeometryCalib.dpx;
    1426     f2=GeometryCalib.f/GeometryCalib.dpy;
    1427     Cx=GeometryCalib.Cx;
    1428     Cy=GeometryCalib.Cy;
    1429     set(handles.fx,'String',num2str(f1,'%1.1f'))
    1430     set(handles.fy,'String',num2str(f2,'%1.1f'))
    1431     set(handles.k,'String',num2str(k,'%1.4f'))
    1432     set(handles.Cx,'String',num2str(Cx,'%1.1f'))
    1433     set(handles.Cy,'String',num2str(Cy,'%1.1f'))
     1500GeometryCalib=s.GeometryCalib
     1501fx=1;fy=1;Cx=0;Cy=0;kc=0; %default
     1502%     Tabchar={};
     1503CoordCell={};
     1504%     kc=0;%default
     1505%     f1=1000;
     1506%     f2=1000;
     1507%     hhuvmat=guidata(findobj(allchild(0),'Name','uvmat'));
     1508%     Cx=str2num(get(hhuvmat.npx,'String'))/2;
     1509%     Cy=str2num(get(hhuvmat.npy,'String'))/2;
     1510Tabchar={};%default
     1511val_cal=1;%default
     1512if ~isempty(GeometryCalib)
     1513    % choose the calibration option
     1514    if isfield(GeometryCalib,'CalibrationType')
     1515       calib_list=get(handles.calib_type,'String');
     1516       for ilist=1:numel(calib_list)
     1517           if strcmp(calib_list{ilist},GeometryCalib.CalibrationType)
     1518               val_cal=ilist;
     1519               break
     1520           end
     1521       end
     1522    end
     1523    display_intrinsic(GeometryCalib,handles)%intrinsic param
     1524    %extrinsic param
     1525    if isfield(GeometryCalib,'Tx_Ty_Tz')
     1526        Tx_Ty_Tz=GeometryCalib.Tx_Ty_Tz;
     1527        set(handles.Tx,'String',num2str(GeometryCalib.Tx_Ty_Tz(1),4))
     1528        set(handles.Ty,'String',num2str(GeometryCalib.Tx_Ty_Tz(2),4))
     1529        set(handles.Tz,'String',num2str(GeometryCalib.Tx_Ty_Tz(3),4))
     1530    end
     1531    if isfield(GeometryCalib,'omc')
     1532        set(handles.Phi,'String',num2str(GeometryCalib.omc(1),4))
     1533        set(handles.Theta,'String',num2str(GeometryCalib.omc(2),4))
     1534        set(handles.Psi,'String',num2str(GeometryCalib.omc(3),4))
     1535    end
    14341536    calib=reshape(GeometryCalib.PointCoord',[],1);
    14351537    for ilist=1:numel(calib)
    … …  
    14381540    CoordCell=reshape(CoordCell,[],5);
    14391541    Tabchar=cell2tab(CoordCell,'    |    ');%transform cells into table ready for display
    1440     set(handles.ListCoord,'Value',1)
    1441     set(handles.ListCoord,'String',Tabchar)
    14421542    MenuPlot_Callback(handles.geometry_calib, [], handles)
    14431543end
     1544set(handles.calib_type,'Value',val_cal)
    14441545Tabchar=[Tabchar;{'......'}];
     1546set(handles.ListCoord,'Value',1)
     1547set(handles.ListCoord,'String',Tabchar)
     1548
    14451549if isempty(CoordCell)% allow mouse action by default in the absence of input points
    14461550    set(handles.edit_append,'Value',1)
    … …  
    14511555end
    14521556
    1453 
    1454 
    1455 
    1456 
    1457 % --------------------------------------------------------------------
    1458 % --- Executes on button press in CLEAR_PTS: clear the list of calibration points
    1459 function CLEAR_PTS_Callback(hObject, eventdata, handles)
    1460 % --------------------------------------------------------------------
    1461 set(handles.ListCoord,'Value',1)% refresh the display of coordinates
    1462 set(handles.ListCoord,'String',{'......'})
    1463 
     1557%------------------------------------------------------------------------
     1558%---display calibration intrinsic parameters
     1559function display_intrinsic(GeometryCalib,handles)
     1560%------------------------------------------------------------------------
     1561fx=[];
     1562fy=[];
     1563if isfield(GeometryCalib,'fx_fy')
     1564    fx=GeometryCalib.fx_fy(1);
     1565    fy=GeometryCalib.fx_fy(2);
     1566end
     1567Cx_Cy=[0 0];%default
     1568if isfield(GeometryCalib,'Cx_Cy')
     1569    Cx_Cy=GeometryCalib.Cx_Cy;
     1570end
     1571kc=0;
     1572if isfield(GeometryCalib,'kc')
     1573    kc=GeometryCalib.kc %* GeometryCalib.focal*GeometryCalib.focal;
     1574end
     1575set(handles.fx,'String',num2str(fx,5))
     1576set(handles.fy,'String',num2str(fy,5))
     1577set(handles.Cx,'String',num2str(Cx_Cy(1),'%1.1f'))
     1578set(handles.Cy,'String',num2str(Cx_Cy(2),'%1.1f'))
     1579set(handles.kc,'String',num2str(kc,'%1.4f'))
     1580
     1581

  • trunk/src/imadoc2struct.m

    r109 r114  
    118118    cont=get(subt,1,'contents');
    119119    if ~isempty(cont)
    120         uid_pixcmx=find(subt,'/GeometryCalib/Pxcmx');
    121         uid_pixcmy=find(subt,'/GeometryCalib/Pxcmy');
    122         if ~isempty(uid_pixcmx) && ~isempty(uid_pixcmy)%NON UTILISE
    123            pixcmx=str2num(get(subt,children(subt,uid_pixcmx),'value'));
    124             if isempty(pixcmx),pixcmx=1;end; %default
    125             pixcmy=str2num(get(subt,children(subt,uid_pixcmy),'value'));
    126             if isempty(pixcmy),pixcmy=1;end; %default
    127             tsai.Pxcmx=pixcmx;
    128             tsai.Pxcmy=pixcmy;
    129         end
     120%         uid_pixcmx=find(subt,'/GeometryCalib/Pxcmx');
     121%         uid_pixcmy=find(subt,'/GeometryCalib/Pxcmy');
     122%         if ~isempty(uid_pixcmx) && ~isempty(uid_pixcmy)%NON UTILISE
     123%            pixcmx=str2num(get(subt,children(subt,uid_pixcmx),'value'));
     124%             if isempty(pixcmx),pixcmx=1;end; %default
     125%             pixcmy=str2num(get(subt,children(subt,uid_pixcmy),'value'));
     126%             if isempty(pixcmy),pixcmy=1;end; %default
     127%             tsai.Pxcmx=pixcmx;
     128%             tsai.Pxcmy=pixcmy;
     129%         end
    130130        %default values:
    131         tsai.f=1;
    132         tsai.dpx=1;
    133         tsai.dpy=1;
    134         tsai.sx=1;
    135         tsai.Cx=0;
    136         tsai.Cy=0;
    137         tsai.Tz=1;
    138         tsai.Tx=0;
    139         tsai.Ty=0;
    140         tsai.R=[1 0 0; 0 1 0; 0 0 0];
    141         tsai.kappa1=0;
     131        %tsai.dpx_dpy=1;
     132        %tsai.dpy=1;
     133        %tsai.sx=1;
     134%         tsai.Cx_Cy=[0 0];
     135       % tsai.Cy=0;
     136%         tsai.Tx_Ty_Tz=[0 0 1];
     137%         tsai.Tx=0;
     138%         tsai.Ty=0;
     139%         tsai.R=[1 0 0; 0 1 0; 0 0 0];
     140        uid_CalibrationType=find(subt,'/GeometryCalib/CalibrationType');
     141        if isequal(length(uid_CalibrationType),1)
     142            tsai.CalibrationType=get(subt,children(subt,uid_CalibrationType),'value');
     143        end
    142144        uid_CoordUnit=find(subt,'/GeometryCalib/CoordUnit');
    143         if ~isempty(uid_CoordUnit)
     145        if isequal(length(uid_CoordUnit),1)
    144146            tsai.CoordUnit=get(subt,children(subt,uid_CoordUnit),'value');
    145147        end
    146         uid_focal=find(subt,'/GeometryCalib/focal');
    147         uid_dpx_dpy=find(subt,'/GeometryCalib/dpx_dpy');
    148         uid_sx=find(subt,'/GeometryCalib/sx');
     148        uid_fx_fy=find(subt,'/GeometryCalib/fx_fy');
     149        focal=[];%default fro old convention (Reg Wilson)
     150        if isequal(length(uid_fx_fy),1)
     151            tsai.fx_fy=str2num(get(subt,children(subt,uid_fx_fy),'value'));
     152        else %old convention (Reg Wilson)
     153            uid_focal=find(subt,'/GeometryCalib/focal');
     154            uid_dpx_dpy=find(subt,'/GeometryCalib/dpx_dpy');
     155            uid_sx=find(subt,'/GeometryCalib/sx');
     156            if ~isempty(uid_focal) && ~isempty(uid_dpx_dpy) && ~isempty(uid_sx)
     157                dpx_dpy=str2num(get(subt,children(subt,uid_dpx_dpy),'value'));
     158                sx=str2num(get(subt,children(subt,uid_sx),'value'));
     159                focal=str2num(get(subt,children(subt,uid_focal),'value'));
     160                tsai.fx_fy(1)=sx*focal/dpx_dpy(1);
     161                tsai.fx_fy(2)=focal/dpx_dpy(2);
     162            end
     163        end
    149164        uid_Cx_Cy=find(subt,'/GeometryCalib/Cx_Cy');
    150         uid_kappa1=find(subt,'/GeometryCalib/kappa1');
     165        if ~isempty(uid_Cx_Cy)
     166            tsai.Cx_Cy=str2num(get(subt,children(subt,uid_Cx_Cy),'value'));
     167        end
     168        uid_kc=find(subt,'/GeometryCalib/kc');
     169        if ~isempty(uid_kc)
     170            tsai.kc=str2num(get(subt,children(subt,uid_kc),'value'));
     171        else %old convention (Reg Wilson)
     172            uid_kappa1=find(subt,'/GeometryCalib/kappa1');
     173            if ~isempty(uid_kappa1)&& ~isempty(focal)
     174                kappa1=str2num(get(subt,children(subt,uid_kappa1),'value'));
     175                tsai.kc=-kappa1*focal*focal;
     176            end
     177        end
    151178        uid_Tx_Ty_Tz=find(subt,'/GeometryCalib/Tx_Ty_Tz');
     179        if ~isempty(uid_Tx_Ty_Tz)
     180            tsai.Tx_Ty_Tz=str2num(get(subt,children(subt,uid_Tx_Ty_Tz),'value'));
     181        end
    152182        uid_R=find(subt,'/GeometryCalib/R');
    153         if ~isempty(uid_focal) && ~isempty(uid_dpx_dpy) && ~isempty(uid_Cx_Cy)
    154             tsai.f=str2num(get(subt,children(subt,uid_focal),'value'));
    155             dpx_dpy=str2num(get(subt,children(subt,uid_dpx_dpy),'value'));
    156             tsai.dpx=dpx_dpy(1);
    157             tsai.dpy=dpx_dpy(2);
    158             if ~isempty(uid_sx)
    159                tsai.sx=str2num(get(subt,children(subt,uid_sx),'value'));
    160             end
    161             Cx_Cy=str2num(get(subt,children(subt,uid_Cx_Cy),'value'));
    162             tsai.Cx=Cx_Cy(1);
    163             tsai.Cy=Cx_Cy(2);
    164         end
    165         if ~isempty(uid_Tx_Ty_Tz)
    166             Tx_Ty_T_char=get(subt,children(subt,uid_Tx_Ty_Tz),'value');
    167             Tx_Ty_Tz=str2num(Tx_Ty_T_char);
    168             tsai.Tx=Tx_Ty_Tz(1);
    169             tsai.Ty=Tx_Ty_Tz(2);
    170             tsai.Tz=Tx_Ty_Tz(3);
    171         end
    172183        if ~isempty(uid_R)
    173184            RR=get(subt,children(subt,uid_R),'value');
    … …  
    176187            end
    177188        end
    178         if ~isempty(uid_kappa1)     
    179             tsai.kappa1=str2num(get(subt,children(subt,uid_kappa1),'value'));
    180         end
     189       
    181190        %look for laser plane definitions   
    182191        uid_Angle=find(subt,'/GeometryCalib/PlaneAngle');
    … …  
    203212                    tsai.NbSlice=str2double(NbSlice);
    204213                end
    205                 tsai.SliceCoord=ones(NbSlice,1)*tsai.SliceCoord+DZ*[0:NbSlice-1]'*[0 0 1];
     214                tsai.SliceCoord=ones(NbSlice,1)*tsai.SliceCoord+DZ*(0:NbSlice-1)'*[0 0 1];
    206215            end         
    207216        end

  • trunk/src/px_XYZ.m

    r109 r114  
    1111
    1212function [X,Y]=px_XYZ(Calib,Xphys,Yphys,Zphys)
    13 X=[];%default
    14 Y=[];
    1513% if exist('Z','var')& isequal(Z,round(Z))& Z>0 & isfield(Calib,'PlanePos')&length(Calib.PlanePos)>=Z
    1614%     Zindex=Z;
    … …  
    2321    Zphys=0;
    2422end
    25 if ~isfield(Calib,'f')
    26     Calib.f=1;
     23if ~isfield(Calib,'fx_fy')
     24     Calib.fx_fy=[1 1];
    2725end
    28 if ~isfield(Calib,'kappa1')
    29     Calib.kappa1=0;
     26if ~isfield(Calib,'Tx_Ty_Tz')
     27     Calib.Tx_Ty_Tz=[0 0 1];
    3028end
    31 if ~isfield(Calib,'sx')
    32     Calib.sx=1;
    33 end
    34 if ~isfield(Calib,'dpx')
    35     Calib.dpx=1;
    36 end
    37 if ~isfield(Calib,'dpy')
    38     Calib.dpy=1;
    39 end
    40 if ~isfield(Calib,'Cx')
    41     Calib.Cx=0;
    42 end
    43 if ~isfield(Calib,'Cy')
    44     Calib.Cy=0;
    45 end
     29% if ~isfield(Calib,'kappa1')
     30%     Calib.kappa1=0;
     31% end
     32% if ~isfield(Calib,'sx')
     33%     Calib.sx=1;
     34% end
     35% if ~isfield(Calib,'dpx')
     36%     Calib.dpx=1;
     37% end
     38% if ~isfield(Calib,'dpy')
     39%     Calib.dpy=1;
     40% end
     41
    4642%%%%%%%%%%%%%
    4743if isfield(Calib,'R')
    4844    R=(Calib.R)';
    49     xc=R(1)*Xphys+R(2)*Yphys+R(3)*Zphys+Calib.Tx;
    50     yc=R(4)*Xphys+R(5)*Yphys+R(6)*Zphys+Calib.Ty;
    51     zc=R(7)*Xphys+R(8)*Yphys+R(9)*Zphys+Calib.Tz;
     45    %camera coordinates
     46    xc=R(1)*Xphys+R(2)*Yphys+R(3)*Zphys+Calib.Tx_Ty_Tz(1);
     47    yc=R(4)*Xphys+R(5)*Yphys+R(6)*Zphys+Calib.Tx_Ty_Tz(2);
     48    zc=R(7)*Xphys+R(8)*Yphys+R(9)*Zphys+Calib.Tx_Ty_Tz(3);
    5249%undistorted image coordinates
    53     Xu=Calib.f*xc./zc;
    54     Yu=Calib.f*yc./zc;
    55 %distorted image coordinates
    56     distortion=(Calib.kappa1)*(Xu.*Xu+Yu.*Yu)+1; %A REVOIR
    57 % distortion=1;
    58     Xd=Xu./distortion;
    59     Yd=Yu./distortion;
     50    Xu=xc./zc;
     51    Yu=yc./zc;
     52%radial quadratic correction factor
     53    if ~isfield(Calib,'kc')
     54        r2=1; %no quadratic distortion
     55    else
     56        r2=1+Calib.kc*(Xu.*Xu+Yu.*Yu);
     57    end
    6058%pixel coordinates
    61     X=Xd*Calib.sx/Calib.dpx+Calib.Cx;
    62     Y=Yd/Calib.dpy+Calib.Cy;
    63 
    64 elseif isfield(Calib,'Pxcmx')&isfield(Calib,'Pxcmy')%old calib 
    65         X=Xphys*Calib.Pxcmx;
    66         Y=Yphys*Calib.Pxcmy;
     59    if ~isfield(Calib,'Cx_Cy')
     60        Calib.Cx_Cy=[0 0];%default value
     61    end
     62    X=Calib.fx_fy(1)*Xu.*r2+Calib.Cx_Cy(1);
     63    Y=Calib.fx_fy(2)*Yu.*r2+Calib.Cx_Cy(2);   
     64%OLD CONVENTION (Wilson)undistorted image coordinates
     65%     Xu=Calib.f*xc./zc;
     66%     Yu=Calib.f*yc./zc;   
     67% %distorted image coordinates
     68%     distortion=(Calib.kappa1)*(Xu.*Xu+Yu.*Yu)+1; %A REVOIR
     69% % distortion=1;
     70%     Xd=Xu./distortion;
     71%     Yd=Yu./distortion;
     72% %pixel coordinates
     73%     X=Xd*Calib.sx/Calib.dpx+Calib.Cx;
     74%     Y=Yd/Calib.dpy+Calib.Cy;
     75else %case 'rescale'
     76    X=Calib.fx_fy(1)*(Xphys+Calib.Tx_Ty_Tz(1));
     77    Y=Calib.fx_fy(2)*(Yphys+Calib.Tx_Ty_Tz(2)); 
    6778end
    6879

  • trunk/src/update_imadoc.m

    r84 r114  
    77% outputfile: xml file to modify
    88%-------------------------------------------------------------
    9 function update_imadoc(GeometryCalib,outputfile)
     9function errormsg=update_imadoc(GeometryCalib,outputfile)
     10errormsg='';
    1011testappend=0;
    1112if exist(outputfile,'file');%=1 if the output file already exists, 0 else 
    … …  
    1920    end
    2021    [success,message]=copyfile(outputfile,backupfile);%make backup
     22    if success==0
     23        errormsg=message
     24    end
    2125    uid=find(t,'ImaDoc');
    2226    if ~isequal(uid,1)
    2327        return
    24     end
    25        
     28    end       
    2629    %if the xml file is  ImaDoc
    2730    uid_calib=find(t,'ImaDoc/GeometryCalib');

  • trunk/src/uvmat.m

    r105 r114  
    815815            set(handles.pxcm,'String','var')
    816816            set(handles.pycm,'String','var')
    817         else
    818             pixcmx=GeometryCalib.f*GeometryCalib.R(1,1)*GeometryCalib.sx/(GeometryCalib.Tz*GeometryCalib.dpx);
    819             pixcmy=GeometryCalib.f*GeometryCalib.R(2,2)/(GeometryCalib.Tz*GeometryCalib.dpy);
     817        elseif isfield(GeometryCalib,'fx_fy')
     818            pixcmx=GeometryCalib.fx_fy(1)%*GeometryCalib.R(1,1)*GeometryCalib.sx/(GeometryCalib.Tz*GeometryCalib.dpx);
     819            pixcmy=GeometryCalib.fx_fy(2)%*GeometryCalib.R(2,2)/(GeometryCalib.Tz*GeometryCalib.dpy);
    820820            set(handles.pxcm,'String',num2str(pixcmx))
    821821            set(handles.pycm,'String',num2str(pixcmy))
    … …  
    47364736            imwrite(imflag,answer,'BitDepth',8);
    47374737        end
    4738         set(list_object_1,'Value',1)
    4739         set(list_object_1,'Max',1)
     4738        set(handles.list_object_1,'Value',1)
     4739        set(handles.list_object_1,'Max',1)
    47404740%     end
    47414741end       
Note: See TracChangeset for help on using the changeset viewer.