Ignore:
Timestamp:
Oct 4, 2010, 10:39:11 PM (14 years ago)
Author:
sommeria
Message:

merge_proj.m, proj_field.m: bugs repaired merging of images
plot_field.m: minor cleaning
imadoc2struct: introduce the possibility of readying the calibration point coordinates
sub_field.m:bugs repaired
geometry_calib, create_grid: introduction of new calibration methods, improvement of detect_grid

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/geometry_calib.m

    r108 r109  
    33% function varargout = geometry_calib(varargin)
    44%
    5 %A%AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
     5%A%AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
    66%  Copyright Joel Sommeria, 2008, LEGI / CNRS-UJF-INPG, sommeria@coriolis-legi.org.
    77%AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
    … …  
    4343% Edit the above text to modify the response to help geometry_calib
    4444
    45 % Last Modified by GUIDE v2.5 25-Mar-2010 19:10:05
     45% Last Modified by GUIDE v2.5 03-Oct-2010 09:34:12
    4646
    4747% Begin initialization code - DO NOT edit
    … …  
    8787    set(hObject,'Position',pos_gui);
    8888end
     89
     90%set menu of calibration options
     91%set(handles.calib_type,'String',{'rescale';'linear';'perspective';'normal';'tsai';'bouguet';'extrinsic'})
     92set(handles.calib_type,'String',{'rescale';'linear';'quadr';'3D_linear';'3D_quadr';'3D_extrinsic'})
    8993inputxml='';
    9094if exist('inputfile','var')& ~isempty(inputfile)
    91 %     [Path,Name,ext]=fileparts(inputfile);
    92 %     form=imformats(ext([2:end]));
    93 %     if ~isempty(form)% if the input file is an image
    94         struct.XmlInputfile=inputfile;
    95         set(hObject,'UserData',struct)
    96         [Pathsub,RootFile,field_count,str2,str_a,str_b,ext,nom_type,subdir]=name2display(inputfile);
    97         inputxml=[fullfile(Pathsub,RootFile) '.xml'];
    98 %     end   
     95    struct.XmlInputfile=inputfile;
     96    set(hObject,'UserData',struct)
     97    [Pathsub,RootFile,field_count,str2,str_a,str_b,ext,nom_type,subdir]=name2display(inputfile);
     98    inputxml=[fullfile(Pathsub,RootFile) '.xml'];
    9999end
    100100set(handles.ListCoord,'String',{'......'})
    … …  
    102102    loadfile(handles,inputxml)% load the point coordiantes existing in the xml file
    103103end
    104 
    105104set(handles.ListCoord,'KeyPressFcn',{@key_press_fcn,handles})%set keyboard action function
    106105
    … …  
    113112varargout{1} = handles.output;
    114113varargout{2}=handles;
    115 
    116 %------------
    117 function Phi_Callback(hObject, eventdata, handles)
    118 
    119 %------------------------------------------------------------------------
    120 %read input xml file and update the edit boxes
    121 function loadfile(handles,fileinput)
    122 %------------------------------------------------------------------------
    123 %read the input xml file
    124 t=xmltree(fileinput);
    125 s=convert(t);%convert to matlab structure
    126 %read data currently displayed on the interface
    127 PointCoord=[];
    128 Coord_cell=get(handles.ListCoord,'String');
    129 data=read_geometry_calib(Coord_cell);
    130 %data=read_geometry_calib(handles);
    131 Coord=[]; %default
    132 if isfield(data,'Coord')
    133     Coord=data.Coord;
    134 end
    135 TabChar_0=get(handles.ListCoord,'String');
    136 nbcoord_0=size(TabChar_0,1);
    137 if isequal(get(handles.edit_append,'Value'),2) %edit mode  A REVOIR
    138     val=get(handles.ListCoord,'Value')-1;
    139 else
    140    val=length(TabChar_0);
    141 end
    142 nbcoord=0;
    143 
    144 %case of calibration (ImaDoc) input file
    145 % hcalib=get(handles.calib_type,'parent');
    146 CalibData=get(handles.geometry_calib,'UserData');
    147 CalibData.XmlInput=fileinput;
    148 if isfield(s,'Heading')
    149     CalibData.Heading=s.Heading;
    150 end
    151 
    152 set(handles.geometry_calib,'UserData',CalibData);%store the heading in the interface 'UserData'
    153 if isfield(s,'GeometryCalib')
    154     Calib=s.GeometryCalib;
    155     if isfield(Calib,'CalibrationType')
    156         CalibrationType=Calib.CalibrationType;
    157         switch CalibrationType
    158             case 'linear'
    159                 set(handles.calib_type,'Value',2)
    160             case 'tsai'
    161                 set(handles.calib_type,'Value',3)
    162             case 'tsai_matlab'
    163                 set(handles.calib_type,'Value',4)
    164         end
    165     end
    166     if isfield(Calib,'SourceCalib')
    167         if isfield(Calib.SourceCalib,'PointCoord')
    168             PointCoord=Calib.SourceCalib.PointCoord;
    169         end
    170     end
    171     nbcoord=length(PointCoord);
    172     if ~isfield(Calib,'ErrorRms')&~isfield(Calib,'ErrorMax') %old convention of Gauthier (cord in mm)
    173         for i=1:length(PointCoord)
    174           line=str2num(PointCoord{i});
    175           Coord(i+val,4:5)=line(4:5);%px x
    176           Coord(i+val,1:3)=line(1:3)/10;%phys x
    177         end
    178     else
    179         for i=1:length(PointCoord)
    180           line=str2num(PointCoord{i});
    181           Coord(i,4:5)=line(4:5);%px x
    182           Coord(i,1:3)=line(1:3);%phys x
    183        end
    184     end
    185 end
    186 %case of xml files of points
    187 if isfield(s,'Coord')
    188     PointCoord=s.Coord;
    189     nbcoord=length(PointCoord);
    190      %case of image coordinates
    191     if isfield(s,'CoordType')& isequal(s.CoordType,'px')
    192         for i=1:nbcoord
    193            line=str2num(PointCoord{i});
    194            Coord(i+val,4:5)=line(1:2);
    195         end
    196      %case of  physical coordinates
    197     else
    198         for i=1:nbcoord
    199            line=str2num(PointCoord{i});
    200            Coord(i+val,1:3)=line(1:3);
    201            nbcolumn=size(Coord,2);
    202            if nbcolumn<5
    203                Coord(i+val,nbcolumn+1:5)=zeros(1,5-nbcolumn);
    204            end
    205         end
    206      end
    207 end
    208 CoordCell={};
    209 for iline=1:size(Coord,1)
    210     for j=1:5
    211         CoordCell{iline,j}=num2str(Coord(iline,j),4);
    212     end
    213 end
    214 % CoordCell=[CoordCell;{' ',' ',' ',' ',' '}];
    215 Tabchar=cell2tab(CoordCell,'    |    ');%transform cells into table ready for display
    216 Tabchar=[Tabchar;{'......'}];
    217 set(handles.ListCoord,'Value',1)
    218 set(handles.ListCoord,'String',Tabchar)
    219 MenuPlot_Callback(handles.geometry_calib, [], handles)
    220 if isempty(Coord)
    221     set(handles.edit_append,'Value',1)
    222     set(handles.edit_append,'BackgroundColor',[1 1 0])
    223 else
    224     set(handles.edit_append,'Value',0)
    225     set(handles.edit_append,'BackgroundColor',[0.7 0.7 0.7])
    226 end
    227114%
    228115%------------------------------------------------------------------------
    … …  
    252139function APPLY_Callback(hObject, eventdata, handles)
    253140%------------------------------------------------------------------------
    254 calib_cell=get(handles.calib_type,'String');
    255 val=get(handles.calib_type,'Value');
    256 calib_type=calib_cell{val};
     141
     142%read the current calibration points
    257143Coord_cell=get(handles.ListCoord,'String');
    258144Object=read_geometry_calib(Coord_cell);
    259 X=Object.Coord(:,1);
    260 Y=Object.Coord(:,2);
    261 Z=Object.Coord(:,3);
    262 if isequal(calib_type,'rescale')
    263     GeometryCalib=calib_rescale(Object.Coord);
    264     Z=0;%Z not taken into account
    265 elseif isequal(calib_type,'linear')
    266     GeometryCalib=calib_linear(Object.Coord);
    267     Z=0; %Z not taken into account
    268 elseif isequal(calib_type,'tsai_cpp')
    269     GeometryCalib=calib_tsai(Object.Coord);
    270 elseif isequal(calib_type,'tsai_matlab')
    271     GeometryCalib=calib_tsai2(Object.Coord);
    272 end
     145Coord=Object.Coord;
     146
     147% apply the calibration, whose type is selected in  handles.calib_type
     148if ~isempty(Coord)
     149    calib_cell=get(handles.calib_type,'String');
     150    val=get(handles.calib_type,'Value');
     151    GeometryCalib=feval(['calib_' calib_cell{val}],Coord,handles);
     152else
     153    msgbox_uvmat('ERROR','No calibration points, abort')
     154    return
     155end
    273156
    274157%check error
    … …  
    298181    Calib.R=GeometryCalib.R;
    299182end
    300 x_ima=Object.Coord(:,4);
    301 y_ima=Object.Coord(:,5);
    302 [Xpoints,Ypoints]=px_XYZ(Calib,X,Y,Z);
    303 GeometryCalib.ErrorRms(1)=sqrt(mean((Xpoints-x_ima).*(Xpoints-x_ima)));
    304 [GeometryCalib.ErrorMax(1),index(1)]=max(abs(Xpoints-x_ima));
    305 GeometryCalib.ErrorRms(2)=sqrt(mean((Ypoints-y_ima).*(Ypoints-y_ima)));
    306 [GeometryCalib.ErrorMax(2),index(2)]=max(abs(Ypoints-y_ima));
    307 [EM,ind_dim]=max(GeometryCalib.ErrorMax);
    308 index=index(ind_dim);
    309 
     183if ~isempty(Coord)
     184    X=Coord(:,1);
     185    Y=Coord(:,2);
     186    Z=Coord(:,3);
     187    x_ima=Coord(:,4);
     188    y_ima=Coord(:,5);
     189    [Xpoints,Ypoints]=px_XYZ(Calib,X,Y,Z);
     190    GeometryCalib.ErrorRms(1)=sqrt(mean((Xpoints-x_ima).*(Xpoints-x_ima)));
     191    [GeometryCalib.ErrorMax(1),index(1)]=max(abs(Xpoints-x_ima));
     192    GeometryCalib.ErrorRms(2)=sqrt(mean((Ypoints-y_ima).*(Ypoints-y_ima)));
     193    [GeometryCalib.ErrorMax(2),index(2)]=max(abs(Ypoints-y_ima));
     194    [EM,ind_dim]=max(GeometryCalib.ErrorMax);
     195    index=index(ind_dim);
     196    if isequal(max(Z),min(Z))%Z constant
     197        GeometryCalib.NbSlice=1;
     198        GeometryCalib.SliceCoord=[0 0 Z(1)];
     199    end
     200end
    310201unitlist=get(handles.CoordUnit,'String');
    311202unit=unitlist{get(handles.CoordUnit,'value')};
    312203GeometryCalib.CoordUnit=unit;
    313 GeometryCalib.SourceCalib.PointCoord=Object.Coord;
     204GeometryCalib.SourceCalib.PointCoord=Coord;
     205
     206%display calibration intrinsic parameters
     207k=0;
     208if isfield(GeometryCalib,'kappa1')
     209    k=GeometryCalib.kappa1 * GeometryCalib.focal*GeometryCalib.focal;
     210end
     211sx=1;dpx_dpy=[1 1];
     212if isfield(GeometryCalib,'sx')
     213    sx=GeometryCalib.sx;
     214end
     215if isfield(GeometryCalib,'dpx_dpy')
     216    dpx_dpy=GeometryCalib.dpx_dpy;
     217end
     218Cx_Cy=[0 0];
     219if isfield(GeometryCalib,'Cx_Cy')
     220    Cx_Cy=GeometryCalib.Cx_Cy;
     221end
     222f1=GeometryCalib.focal*sx/dpx_dpy(1);
     223f2=GeometryCalib.focal/dpx_dpy(2);
     224Cx=Cx_Cy(1);
     225Cy=Cx_Cy(2);
     226set(handles.fx,'String',num2str(f1,4))
     227set(handles.fy,'String',num2str(f2,4))
     228set(handles.k,'String',num2str(k,'%1.4f'))
     229set(handles.Cx,'String',num2str(Cx,'%1.1f'))
     230set(handles.Cy,'String',num2str(Cy,'%1.1f'))
     231% Display extrinsinc parameters
     232set(handles.Tx,'String',num2str(GeometryCalib.Tx_Ty_Tz(1),4))
     233set(handles.Ty,'String',num2str(GeometryCalib.Tx_Ty_Tz(2),4))
     234set(handles.Tz,'String',num2str(GeometryCalib.Tx_Ty_Tz(3),4))
     235
     236% indicate the plane of the calibration grid if defined
    314237huvmat=findobj(allchild(0),'Name','uvmat');
    315238hhuvmat=guidata(huvmat);%handles of elements in the GUI uvmat
    … …  
    337260    hhh=findobj(hhuvmat.axes3,'Tag','calib_marker');% delete calib points and markers
    338261    if ~isempty(hhh)
    339         delete(hhh);
     262        delete(hhh);     
    340263    end
    341264    hhh=findobj(hhuvmat.axes3,'Tag','calib_points');
    … …  
    346269    set(hhuvmat.FixedLimits,'BackgroundColor',[0.7 0.7 0.7])
    347270    uvmat('RootPath_Callback',hObject,eventdata,hhuvmat); %file input with xml reading  in uvmat
    348     MenuPlot_Callback(hObject, eventdata, handles)
    349     set(handles.ListCoord,'Value',index)% indicate in the list the point with max deviation (possible mistake)
    350     ListCoord_Callback(hObject, eventdata, handles)
     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
    351276    figure(handles.geometry_calib)
    352277end
    … …  
    355280% --- Executes on button press in calibrate_lin.
    356281function REPLICATE_Callback(hObject, eventdata, handles)
     282% %%%%%%Todo: correct on the model of APPLY_Callback%%%%%%%%%%
    357283%------------------------------------------------------------------------
    358284calib_cell=get(handles.calib_type,'String');
    359285val=get(handles.calib_type,'Value');
    360 calib_type=calib_cell{val};
    361286Coord_cell=get(handles.ListCoord,'String');
    362287Object=read_geometry_calib(Coord_cell);
    363 
    364 if isequal(calib_type,'rescale')
    365     GeometryCalib=calib_rescale(Object.Coord);
    366 elseif isequal(calib_type,'linear')
    367     GeometryCalib=calib_linear(Object.Coord);
    368 elseif isequal(calib_type,'tsai')
    369     GeometryCalib=calib_tsai(Object.Coord);
    370 end
     288GeometryCalib=feval(calib_cell{val},Object.Coord);
     289
    371290% %record image source
    372291GeometryCalib.SourceCalib.PointCoord=Object.Coord;
    … …  
    424343%------------------------------------------------------------------------
    425344% determine the parameters for a calibration by an affine function (rescaling and offset, no rotation)
    426 function GeometryCalib=calib_rescale(Coord)
    427 %------------------------------------------------------------------------
    428  
     345function GeometryCalib=calib_rescale(Coord,handles)
     346%------------------------------------------------------------------------
    429347X=Coord(:,1);
    430 Y=Coord(:,2);
     348Y=Coord(:,2);% Z not used
    431349x_ima=Coord(:,4);
    432350y_ima=Coord(:,5);
    … …  
    436354T_y=py(2);
    437355GeometryCalib.CalibrationType='rescale';
    438 GeometryCalib.focal=1;
     356GeometryCalib.focal=1;%default
     357GeometryCalib.sx=px(1)/py(1);
    439358GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function
    440 GeometryCalib.Tx_Ty_Tz=[T_x T_y 1];
    441 GeometryCalib.R=[px(1),0,0;0,py(1),0;0,0,1];
    442 
    443 % %check error
    444 % Calib.dpx=1;
    445 % Calib.dpy=1;
    446 % Calib.sx=1;
    447 % Calib.Cx=0;
    448 % Calib.Cy=0;
    449 % Calib.Tz=1;
    450 % Calib.kappa1=0;
    451 % Calib.f=GeometryCalib.focal;
    452 % Calib.Tx=T_x;
    453 % Calib.Ty=T_y;
    454 % Calib.R=GeometryCalib.R;
    455 % [Xpoints,Ypoints]=px_XYZ(Calib,X,Y,0);
    456 % GeometryCalib.ErrorRms(1)=sqrt(mean((Xpoints-x_ima).*(Xpoints-x_ima)));
    457 % GeometryCalib.ErrorMax(1)=max(abs(Xpoints-x_ima));
    458 % GeometryCalib.ErrorRms(2)=sqrt(mean((Ypoints-y_ima).*(Ypoints-y_ima)));
    459 % GeometryCalib.ErrorMax(2)=max(abs(Ypoints-y_ima));
     359GeometryCalib.Tx_Ty_Tz=[T_x T_y GeometryCalib.focal/py(1)];
     360GeometryCalib.R=[1,0,0;0,1,0;0,0,0];
    460361
    461362%------------------------------------------------------------------------
    462363% determine the parameters for a calibration by a linear transform matrix (rescale and rotation)
    463 function GeometryCalib=calib_linear(Coord)
     364function GeometryCalib=calib_linear(Coord,handles)
    464365%------------------------------------------------------------------------
    465366X=Coord(:,1);
    466 Y=Coord(:,2);
     367Y=Coord(:,2);% Z not used
    467368x_ima=Coord(:,4);
    468369y_ima=Coord(:,5);
    469370XY_mat=[ones(size(X)) X Y];
    470371a_X1=XY_mat\x_ima; %transformation matrix for X
    471 x1=XY_mat*a_X1;%reconstruction
    472 err_X1=max(abs(x1-x_ima));%error
     372% x1=XY_mat*a_X1;%reconstruction
     373% err_X1=max(abs(x1-x_ima));%error
    473374a_Y1=XY_mat\y_ima;%transformation matrix for X
    474 y1=XY_mat*a_Y1;
    475 err_Y1=max(abs(y1-y_ima));%error
    476 T_x=a_X1(1);
    477 T_y=a_Y1(1);
     375% y1=XY_mat*a_Y1;
     376% err_Y1=max(abs(y1-y_ima));%error
    478377GeometryCalib.CalibrationType='linear';
    479378GeometryCalib.focal=1;
    480379GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function
    481 GeometryCalib.Tx_Ty_Tz=[T_x T_y 1];
    482 GeometryCalib.R=[a_X1(2),a_X1(3),0;a_Y1(2),a_Y1(3),0;0,0,1];
    483 
    484 % %check error
    485 % GeometryCalib.ErrorRms(1)=sqrt(mean((x1-x_ima).*(x1-x_ima)));
    486 % GeometryCalib.ErrorMax(1)=max(abs(x1-x_ima));
    487 % GeometryCalib.ErrorRms(2)=sqrt(mean((y1-y_ima).*(y1-y_ima)));
    488 % GeometryCalib.ErrorMax(2)=max(abs(y1-y_ima));
    489 
    490 %------------------------------------------------------------------------
    491 function GeometryCalib=calib_tsai2(Coord)
     380R=[a_X1(2),a_X1(3),0;a_Y1(2),a_Y1(3),0;0,0,0];
     381norm=det(R);
     382GeometryCalib.Tx_Ty_Tz=[a_X1(1) a_Y1(1) norm];
     383GeometryCalib.R=R/norm;
     384
     385%------------------------------------------------------------------------
     386% determine the tsai parameters for a view normal to the grid plane
     387function GeometryCalib=calib_normal(Coord,handles)
     388%------------------------------------------------------------------------
     389Calib.f1=str2num(get(handles.fx,'String'));
     390Calib.f2=str2num(get(handles.fy,'String'));
     391Calib.k=str2num(get(handles.k,'String'));
     392Calib.Cx=str2num(get(handles.Cx,'String'));
     393Calib.Cy=str2num(get(handles.Cy,'String'));
     394%default
     395if isempty(Calib.f1)
     396    Calib.f1=25/0.012;
     397end
     398if isempty(Calib.f2)
     399    Calib.f2=25/0.012;
     400end
     401if isempty(Calib.k)
     402    Calib.k=0;
     403end
     404if isempty(Calib.Cx)||isempty(Calib.Cy)
     405    huvmat=findobj(allchild(0),'Tag','uvmat');
     406    hhuvmat=guidata(huvmat);
     407    Calib.Cx=str2num(get(hhuvmat.npx,'String'))/2;
     408    Calib.Cx=str2num(get(hhuvmat.npy,'String'))/2;
     409end   
     410%tsai parameters
     411Calib.dpx=0.012;%arbitrary
     412Calib.dpy=0.012;
     413Calib.sx=Calib.f1*Calib.dpx/(Calib.f2*Calib.dpy);
     414Calib.f=Calib.f2*Calib.dpy;
     415Calib.kappa1=Calib.k/(Calib.f*Calib.f);
     416
     417%initial guess
     418X=Coord(:,1);
     419Y=Coord(:,2);
     420Zmean=mean(Coord(:,3));
     421x_ima=Coord(:,4)-Calib.Cx;
     422y_ima=Coord(:,5)-Calib.Cy;
     423XY_mat=[ones(size(X)) X Y];
     424a_X1=XY_mat\x_ima; %transformation matrix for X
     425a_Y1=XY_mat\y_ima;%transformation matrix for Y
     426R=[a_X1(2),a_X1(3),0;a_Y1(2),a_Y1(3),0;0,0,-1];% rotation+ z axis reversal (upward)
     427norm=sqrt(det(-R));
     428calib_param(1)=0;% quadratic distortion
     429calib_param(2)=a_X1(1);
     430calib_param(3)=a_Y1(1);
     431calib_param(4)=Calib.f/(norm*Calib.dpx)-R(3,3)*Zmean;
     432calib_param(5)=angle(a_X1(2)+i*a_X1(3));
     433display(['initial guess=' num2str(calib_param)])
     434
     435%optimise the parameters: minimisation of error
     436calib_param = fminsearch(@(calib_param) error_calib(calib_param,Calib,Coord),calib_param)
     437
     438GeometryCalib.CalibrationType='tsai_normal';
     439GeometryCalib.focal=Calib.f;
     440GeometryCalib.dpx_dpy=[Calib.dpx Calib.dpy];
     441GeometryCalib.Cx_Cy=[Calib.Cx Calib.Cy];
     442GeometryCalib.sx=Calib.sx;
     443GeometryCalib.kappa1=calib_param(1);
     444GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function
     445GeometryCalib.Tx_Ty_Tz=[calib_param(2) calib_param(3) calib_param(4)];
     446alpha=calib_param(5);
     447GeometryCalib.R=[cos(alpha) sin(alpha) 0;-sin(alpha) cos(alpha) 0;0 0 -1];
     448
     449%------------------------------------------------------------------------
     450function GeometryCalib=calib_3D_linear(Coord,handles)
    492451%------------------------------------------------------------------
    493452path_uvmat=which('uvmat');% check the path detected for source file uvmat
    … …  
    495454huvmat=findobj(allchild(0),'Tag','uvmat');
    496455hhuvmat=guidata(huvmat);
     456x_1=Coord(:,4:5)'
     457X_1=Coord(:,1:3)'
     458n_ima=1;
     459% check_cond=0;
     460
     461nx=str2num(get(hhuvmat.npx,'String'));
     462ny=str2num(get(hhuvmat.npy,'String'));
     463
     464est_dist=[0;0;0;0;0];
     465est_aspect_ratio=0;
     466est_fc=[1;1];
     467%fc=[25;25]/0.012;
     468center_optim=0;
     469run(fullfile(path_UVMAT,'toolbox_calib','go_calib_optim'));
     470
     471GeometryCalib.CalibrationType='3D_linear';
     472GeometryCalib.focal=fc(2);
     473GeometryCalib.dpx_dpy=[1 1];
     474GeometryCalib.Cx_Cy=cc';
     475GeometryCalib.sx=fc(1)/fc(2);
     476GeometryCalib.kappa1=-kc(1)/fc(2)^2;
     477GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function
     478GeometryCalib.Tx_Ty_Tz=Tc_1';
     479GeometryCalib.R=Rc_1;
     480
     481%------------------------------------------------------------------------
     482function GeometryCalib=calib_3D_quadr(Coord,handles)
     483%------------------------------------------------------------------
     484
     485path_uvmat=which('uvmat');% check the path detected for source file uvmat
     486path_UVMAT=fileparts(path_uvmat); %path to UVMAT
     487huvmat=findobj(allchild(0),'Tag','uvmat');
     488hhuvmat=guidata(huvmat);
     489% check_cond=0;
     490coord_files=get(handles.coord_files,'String');
     491if ischar(coord_files)
     492    coord_files={coord_files};
     493end
     494if isempty(coord_files{1}) || isequal(coord_files,{''})
     495    coord_files={};
     496end
     497
     498%retrieve the calibration points stored in the files listed in the popup list coord_files
     499x_1=Coord(:,4:5)';
     500X_1=Coord(:,1:3)';
     501n_ima=numel(coord_files)+1;
     502if ~isempty(coord_files)
     503    msgbox_uvmat('CONFIRMATION',['The xy coordinates of the calibration points in ' num2str(n_ima) ' planes will be used'])
     504    for ifile=1:numel(coord_files)
     505    t=xmltree(coord_files{ifile});
     506    s=convert(t);%convert to matlab structure
     507        if isfield(s,'GeometryCalib')
     508            if isfield(s.GeometryCalib,'SourceCalib')
     509                if isfield(s.GeometryCalib.SourceCalib,'PointCoord')
     510                PointCoord=s.GeometryCalib.SourceCalib.PointCoord;
     511                Coord_file=[];
     512                for i=1:length(PointCoord)
     513                    line=str2num(PointCoord{i});
     514                    Coord_file(i,4:5)=line(4:5);%px x
     515                    Coord_file(i,1:3)=line(1:3);%phys x
     516                end
     517                eval(['x_' num2str(ifile+1) '=Coord_file(:,4:5)'';']);
     518                eval(['X_' num2str(ifile+1) '=Coord_file(:,1:3)'';']);
     519                end
     520            end
     521        end
     522    end
     523end
     524
     525n_ima=numel(coord_files)+1;
     526whos
     527nx=str2num(get(hhuvmat.npx,'String'));
     528ny=str2num(get(hhuvmat.npy,'String'));
     529
     530est_dist=[1;0;0;0;0];
     531est_aspect_ratio=1;
     532%est_fc=[0;0];
     533%fc=[25;25]/0.012;
     534center_optim=0;
     535run(fullfile(path_UVMAT,'toolbox_calib','go_calib_optim'));
     536
     537GeometryCalib.CalibrationType='3D_quadr';
     538GeometryCalib.focal=fc(2);
     539GeometryCalib.dpx_dpy=[1 1];
     540GeometryCalib.Cx_Cy=cc';
     541GeometryCalib.sx=fc(1)/fc(2);
     542GeometryCalib.kappa1=-kc(1)/fc(2)^2;
     543GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function
     544GeometryCalib.Tx_Ty_Tz=Tc_1';
     545GeometryCalib.R=Rc_1;
     546GeometryCalib.R(3,1:3)=-GeometryCalib.R(3,1:3);%inversion for z upward
     547GeometryCalib.ErrorRMS=[];
     548GeometryCalib.ErrorMax=[];
     549
     550
     551%------------------------------------------------------------------------
     552function GeometryCalib=calib_3D_extrinsic(Coord,handles)
     553%------------------------------------------------------------------
     554path_uvmat=which('geometry_calib');% check the path detected for source file uvmat
     555path_UVMAT=fileparts(path_uvmat); %path to UVMAT
    497556x_1=Coord(:,4:5)';
    498557X_1=Coord(:,1:3)';
    499558n_ima=1;
    500 % check_cond=0;
    501 
    502 
    503 nx=str2num(get(hhuvmat.npx,'String'));
    504 ny=str2num(get(hhuvmat.npy,'String'));
    505 
    506 
    507 % est_kc=[1;0;0;0;0];
    508 est_fc=0;
    509 est_dist=[1;0;0;0;0];
    510 run(fullfile(path_UVMAT,'toolbox_calib','go_calib_optim'));
    511 
    512 GeometryCalib.CalibrationType='tsai_matlab';
    513 GeometryCalib.focal=f(2);
     559Calib.f1=str2num(get(handles.fx,'String'));
     560Calib.f2=str2num(get(handles.fy,'String'));
     561Calib.k=str2num(get(handles.k,'String'));
     562Calib.Cx=str2num(get(handles.Cx,'String'));
     563Calib.Cy=str2num(get(handles.Cy,'String'));
     564
     565fct_path=fullfile(path_UVMAT,'toolbox_calib');
     566addpath(fct_path)
     567% [omc1,Tc1,Rc1,H,x,ex,JJ] = compute_extrinsic(x_1,X_1,...
     568%     [Calib.f Calib.f*Calib.sx]',...
     569%     [Calib.Cx Calib.Cy]',...
     570%     [-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]);
     575%get the euler angles ???
     576rmpath(fct_path);
     577
     578std(ex')
     579GeometryCalib.CalibrationType='3D_extrinsic';
     580GeometryCalib.focal=Calib.f2;
    514581GeometryCalib.dpx_dpy=[1 1];
    515 GeometryCalib.Cx_Cy=cc';
    516 GeometryCalib.sx=fc(1)/fc(2);
    517 GeometryCalib.kappa1=-k(1)/f(2)^2;
     582GeometryCalib.Cx_Cy=[Calib.Cx Calib.Cy];
     583GeometryCalib.sx=Calib.f1/Calib.f2;
     584GeometryCalib.kappa1=Calib.k/(Calib.f2*Calib.f2);
    518585GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function
    519 GeometryCalib.Tx_Ty_Tz=Tc_1';
    520 GeometryCalib.R=Rc_1;
    521 % Calib.dpx=GeometryCalib.dpx_dpy(1);
    522 % Calib.dpy=GeometryCalib.dpx_dpy(2);
    523 % Calib.sx=GeometryCalib.sx;
    524 % Calib.Cx=GeometryCalib.Cx_Cy(1);
    525 % Calib.Cy=GeometryCalib.Cx_Cy(2);
    526 % Calib.kappa1=GeometryCalib.kappa1;
    527 % Calib.f=GeometryCalib.focal;
    528 % Calib.Tx=GeometryCalib.Tx_Ty_Tz(1);
    529 % Calib.Ty=GeometryCalib.Tx_Ty_Tz(2);
    530 % Calib.Tz=GeometryCalib.Tx_Ty_Tz(3);
    531 % Calib.R=GeometryCalib.R;
    532 % X=Coord(:,1);
    533 % Y=Coord(:,2);
    534 % Z=Coord(:,3);
    535 % x_ima=Coord(:,4);
    536 % y_ima=Coord(:,5);
    537 % [Xpoints,Ypoints]=px_XYZ(Calib,X,Y,Z);
    538 %
    539 % GeometryCalib.ErrorRms(1)=sqrt(mean((Xpoints-x_ima).*(Xpoints-x_ima)));
    540 % [GeometryCalib.ErrorMax(1),GeometryCalib.IndexMax(1)]=max(abs(Xpoints-x_ima));
    541 % GeometryCalib.ErrorRms(2)=sqrt(mean((Ypoints-y_ima).*(Ypoints-y_ima)));
    542 % [GeometryCalib.ErrorMax(2),GeometryCalib.IndexMax(2)]=max(abs(Ypoints-y_ima));
    543 
    544 function GeometryCalib=calib_tsai(Coord)
     586GeometryCalib.Tx_Ty_Tz=Tc1';
     587%inversion of z axis
     588GeometryCalib.R=Rc1;
     589%GeometryCalib.R(3,1:3)=-GeometryCalib.R(3,1:3);%inversion for z upward
     590
     591
     592%------------------------------------------------------------------------
     593%function GeometryCalib=calib_tsai_heikkila(Coord)
     594% TEST: NOT IMPLEMENTED
     595%------------------------------------------------------------------
     596% path_uvmat=which('uvmat');% check the path detected for source file uvmat
     597% path_UVMAT=fileparts(path_uvmat); %path to UVMAT
     598% path_calib=fullfile(path_UVMAT,'toolbox_calib_heikkila');
     599% addpath(path_calib)
     600% npoints=size(Coord,1);
     601% Coord(:,1:3)=10*Coord(:,1:3);
     602% Coord=[Coord zeros(npoints,2) -ones(npoints,1)];
     603% [par,pos,iter,res,er,C]=cacal('dalsa',Coord);
     604% GeometryCalib.CalibrationType='tsai';
     605% GeometryCalib.focal=par(2);
     606
     607
     608%--------------------------------------------------------------------------
     609function GeometryCalib=calib_tsai(Coord,handles)% old version using gauthier's bianry ccal_fo
    545610%------------------------------------------------------------------------
    546611%TSAI
    547 % 'calibration_lin' provides a linear transform on coordinates,
    548612path_uvmat=which('uvmat');% check the path detected for source file uvmat
    549613path_UVMAT=fileparts(path_uvmat); %path to UVMAT
    550 % if isunix
    551     %fid = fopen(fullfile(path_UVMAT,'PARAM_LINUX.txt'),'r');%open the file with civ binary names
    552 xmlfile=fullfile(path_UVMAT,'PARAM.xml');
     614xmlfile=fullfile(path_UVMAT,'PARAM.xml');%name of the file containing names of binary executables
    553615if exist(xmlfile,'file')
    554     t=xmltree(xmlfile);
    555     sparam=convert(t);
     616    t=xmltree(xmlfile);% read the (xml) file containing names of binary executables
     617    sparam=convert(t);% convert to matlab structure
    556618end
    557619if ~isfield(sparam,'GeometryCalibBin')
    … …  
    570632textcoord=num2str(Coord,4);
    571633dlmwrite('t.txt',textcoord,''); 
    572 % ['!' Tsai_exe ' -f1 0 -f2 t.txt']
    573     eval(['!' Tsai_exe ' -f t.txt > tsaicalib.log']);
     634% ['!' Tsai_exe ' -fx 0 -fy t.txt']
     635eval(['!' Tsai_exe ' -f t.txt > tsaicalib.log']);
    574636if ~exist('calib.dat','file')
    575637    msgbox_uvmat('ERROR','no output from calibration program Tsai_exe: possibly too few points')
    … …  
    577639calibdat=dlmread('calib.dat');
    578640delete('calib.dat')
    579 delete('t.txt')
     641%delete('t.txt')
    580642GeometryCalib.CalibrationType='tsai';
    581643GeometryCalib.focal=calibdat(10);
    … …  
    604666%EN DEDUIRE MATRICE R ??
    605667GeometryCalib.R=[r1,r2,r3;r4,r5,r6;r7,r8,r9];
    606 %erreur a caracteriser?
    607 % %check error
    608 % Calib.dpx=GeometryCalib.dpx_dpy(1);
    609 % Calib.dpy=GeometryCalib.dpx_dpy(2);
    610 % Calib.sx=GeometryCalib.sx;
    611 % Calib.Cx=GeometryCalib.Cx_Cy(1);
    612 % Calib.Cy=GeometryCalib.Cx_Cy(2);
    613 % Calib.kappa1=GeometryCalib.kappa1;
    614 % Calib.f=GeometryCalib.focal;
    615 % Calib.Tx=GeometryCalib.Tx_Ty_Tz(1);
    616 % Calib.Ty=GeometryCalib.Tx_Ty_Tz(2);
    617 % Calib.Tz=GeometryCalib.Tx_Ty_Tz(3);
    618 % Calib.R=GeometryCalib.R;
    619 % X=Coord(:,1);
    620 % Y=Coord(:,2);
    621 % Z=Coord(:,3);
    622 % x_ima=Coord(:,4);
    623 % y_ima=Coord(:,5);
    624 % [Xpoints,Ypoints]=px_XYZ(Calib,X,Y,Z);
    625 %
    626 % GeometryCalib.ErrorRms(1)=sqrt(mean((Xpoints-x_ima).*(Xpoints-x_ima)));
    627 % GeometryCalib.ErrorMax(1)=max(abs(Xpoints-x_ima));
    628 % GeometryCalib.ErrorRms(2)=sqrt(mean((Ypoints-y_ima).*(Ypoints-y_ima)));
    629 % GeometryCalib.ErrorMax(2)=max(abs(Ypoints-y_ima));
    630 
    631 
    632 %------------------------------------------------------------------------
    633 % --- Executes on button press in rotation.
    634 function rotation_Callback(hObject, eventdata, handles)
    635 %------------------------------------------------------------------------
    636 angle_rot=(pi/180)*str2num(get(handles.Phi,'String'));
     668
     669%------------------------------------------------------------------------
     670% --- determine the rms of calibration error
     671function ErrorRms=error_calib(calib_param,Calib,Coord)
     672%calib_param: vector of free calibration parameters (to optimise)
     673%Calib: structure of the given calibration parameters
     674%Coord: list of phys coordinates (columns 1-3, and pixel coordinates (columns 4-5)
     675Calib.f=25;
     676Calib.dpx=0.012;
     677Calib.dpy=0.012;
     678Calib.sx=1;
     679Calib.Cx=512;
     680Calib.Cy=512;
     681Calib.kappa1=calib_param(1);
     682Calib.Tx=calib_param(2);
     683Calib.Ty=calib_param(3);
     684Calib.Tz=calib_param(4);
     685alpha=calib_param(5);
     686Calib.R=[cos(alpha) sin(alpha) 0;-sin(alpha) cos(alpha) 0;0 0 -1];
     687
     688X=Coord(:,1);
     689Y=Coord(:,2);
     690Z=Coord(:,3);
     691x_ima=Coord(:,4);
     692y_ima=Coord(:,5);
     693[Xpoints,Ypoints]=px_XYZ(Calib,X,Y,Z);
     694ErrorRms(1)=sqrt(mean((Xpoints-x_ima).*(Xpoints-x_ima)));
     695ErrorRms(2)=sqrt(mean((Ypoints-y_ima).*(Ypoints-y_ima)));
     696ErrorRms=mean(ErrorRms);
     697
     698%------------------------------------------------------------------------
     699function XImage_Callback(hObject, eventdata, handles)
     700%------------------------------------------------------------------------
     701update_list(hObject, eventdata,handles)
     702
     703%------------------------------------------------------------------------
     704function YImage_Callback(hObject, eventdata, handles)
     705%------------------------------------------------------------------------
     706update_list(hObject, eventdata,handles)
     707
     708%------------------------------------------------------------------------
     709% --- Executes on button press in STORE.
     710function STORE_Callback(hObject, eventdata, handles)
    637711Coord_cell=get(handles.ListCoord,'String');
    638 data=read_geometry_calib(Coord_cell);
    639 data.Coord(:,1)=cos(angle_rot)*data.Coord(:,1)+sin(angle_rot)*data.Coord(:,2);
    640 data.Coord(:,1)=-sin(angle_rot)*data.Coord(:,1)+cos(angle_rot)*data.Coord(:,2);
    641 set(handles.XObject,'String',num2str(data.Coord(:,1),4));
    642 set(handles.YObject,'String',num2str(data.Coord(:,2),4));
    643 
    644 %------------------------------------------------------------------------
    645 function XImage_Callback(hObject, eventdata, handles)
     712Object=read_geometry_calib(Coord_cell);
     713unitlist=get(handles.CoordUnit,'String');
     714unit=unitlist{get(handles.CoordUnit,'value')};
     715GeometryCalib.CoordUnit=unit;
     716GeometryCalib.SourceCalib.PointCoord=Object.Coord;
     717huvmat=findobj(allchild(0),'Name','uvmat');
     718hhuvmat=guidata(huvmat);%handles of elements in the GUI uvmat
     719RootPath='';
     720RootFile='';
     721if ~isempty(hhuvmat.RootPath)& ~isempty(hhuvmat.RootFile)
     722    testhandle=1;
     723    RootPath=get(hhuvmat.RootPath,'String');
     724    RootFile=get(hhuvmat.RootFile,'String');
     725    filebase=fullfile(RootPath,RootFile);
     726    while exist([filebase '.xml'],'file')
     727        filebase=[filebase '~'];
     728    end
     729    outputfile=[filebase '.xml'];
     730    update_imadoc(GeometryCalib,outputfile)
     731    listfile=get(handles.coord_files,'string');
     732    if isequal(listfile,{''})
     733        listfile={outputfile};
     734    else
     735        listfile=[listfile;{outputfile}]%update the list of coord files
     736    end
     737    set(handles.coord_files,'string',listfile);
     738end
     739set(handles.ListCoord,'Value',1)% refresh the display of coordinates
     740set(handles.ListCoord,'String',{'......'})
     741
     742%------------------------------------------------------------------------
     743% --- Executes on button press in CLEAR.
     744function CLEAR_Callback(hObject, eventdata, handles)
     745%------------------------------------------------------------------------
     746set(handles.coord_files,'Value',1)
     747set(handles.coord_files,'String',{''})
     748
     749%------------------------------------------------------------------------
     750function XObject_Callback(hObject, eventdata, handles)
    646751%------------------------------------------------------------------------
    647752update_list(hObject, eventdata,handles)
    648753
    649754%------------------------------------------------------------------------
    650 function YImage_Callback(hObject, eventdata, handles)
     755function YObject_Callback(hObject, eventdata, handles)
    651756%------------------------------------------------------------------------
    652757update_list(hObject, eventdata,handles)
    653758
    654 function XObject_Callback(hObject, eventdata, handles)
    655 update_list(hObject, eventdata,handles)
    656 
    657 function YObject_Callback(hObject, eventdata, handles)
    658 update_list(hObject, eventdata,handles)
    659 
     759%------------------------------------------------------------------------
    660760function ZObject_Callback(hObject, eventdata, handles)
     761%------------------------------------------------------------------------
    661762update_list(hObject, eventdata,handles)
    662763
    … …  
    748849%------------------------------------------------------------------------
    749850choice=get(handles.edit_append,'Value');
    750 % if choice==1
    751 %        Coord=get(handles.ListCoord,'String');
    752 %        val=length(Coord);
    753 %        if val>=1 & isequal(Coord{val},'')
    754 %             val=val-1; %do not take into account blank
    755 %        end
    756 %        Coord{val+1}='';
    757 %        set(handles.ListCoord,'String',Coord)
    758 %        set(handles.ListCoord,'Value',val+1)
    759 % end
    760851if choice
    761852    set(handles.edit_append,'BackgroundColor',[1 1 0])
    … …  
    817908end
    818909
    819 % %------------------------------------------------------------------------
    820 % % --- Executes on button press in append_point.
    821 % function append_point_Callback(hObject, eventdata, handles)
    822 % %------------------------------------------------------------------------
    823 %        Coord=get(handles.ListCoord,'String');
    824 %        val=length(Coord);
    825 %        if val>=1 & isequal(Coord{val},'')
    826 %             val=val-1; %do not take into account blank
    827 %        end
    828 %        Coord{val+1}='';
    829 %        set(handles.ListCoord,'String',Coord)
    830 %        set(handles.ListCoord,'Value',val+1)
    831 
    832 %------------------------------------------------------------------------
    833 function MenuOpen_Callback(hObject, eventdata, handles)
    834 %------------------------------------------------------------------------
    835 %get the object file
    836 huvmat=findobj(allchild(0),'Name','uvmat');
    837 UvData=get(huvmat,'UserData');
    838 hchild=get(huvmat,'Children');
    839 hrootpath=findobj(hchild,'Tag','RootPath');
    840 oldfile=get(hrootpath,'String');
    841 if isempty(oldfile)
    842     oldfile='';
    843 end
    844 %[FileName,PathName] = uigetfile('*.civ','Select a .civ file',oldfile)
    845 [FileName, PathName, filterindex] = uigetfile( ...
    846        {'*.xml;*.mat', ' (*.xml,*.mat)';
    847        '*.xml',  '.xml files '; ...
    848         '*.mat',  '.mat matlab files '}, ...
    849         'Pick a file',oldfile);
    850 fileinput=[PathName FileName];%complete file name
    851 testblank=findstr(fileinput,' ');%look for blanks
    852 if ~isempty(testblank)
    853     msgbox_uvmat('ERROR','forbidden input file name or path: no blank character allowed')
    854     return
    855 end
    856 sizf=size(fileinput);
    857 if (~ischar(fileinput)|~isequal(sizf(1),1)),return;end
    858 loadfile(handles,fileinput)
    859 
    860 
    861910%------------------------------------------------------------------------
    862911function MenuPlot_Callback(hObject, eventdata, handles)
    … …  
    9551004    grid_input=CalibData.grid;%retrieve the previously used grid
    9561005end
    957 [T,CalibData.grid]=create_grid(grid_input,'detect grid');%display the GUI create_grid, read the set of phys coordinates T
     1006[T,CalibData.grid]=create_grid(grid_input,'detect_grid');%display the GUI create_grid, read the set of phys coordinates T
     1007
    9581008set(handles.geometry_calib,'UserData',CalibData)%store the phys grid parameters for later use
    9591009
    … …  
    9881038X=[CalibData.grid.x_0 CalibData.grid.x_1 CalibData.grid.x_0 CalibData.grid.x_1]';%corner absissa in the phys coordinates
    9891039Y=[CalibData.grid.y_0 CalibData.grid.y_0 CalibData.grid.y_1 CalibData.grid.y_1]';%corner ordinates in the phys coordinates
    990 %X=[CalibData.grid.x_0 CalibData.grid.x_0 CalibData.grid.x_1 CalibData.grid.x_1]';%corner absissa in the phys coordinates
    991 %Y=[CalibData.grid.y_0 CalibData.grid.y_1 CalibData.grid.y_1 CalibData.grid.y_0]';%corner ordinates in the phys coordinate
    992 %XY_mat=[ones(size(X)) X Y];
    993 %a_X1=XY_mat\corners_X; %transformation matrix for X
    994 %x1=XY_mat*a_X1;%reconstruction
    995 %err_X1=max(abs(x1-corners_X));%error
    996 %a_Y1=XY_mat\corners_Y;%transformation matrix for X
    997 %y1=XY_mat*a_Y1;
    998 %err_Y1=max(abs(y1-corners_Y));%error
    999 % figure(1)
    1000 % hold off
    1001 % plot([corners_X;corners_X(1)],[corners_Y;corners_Y(1)],'r')
    1002 %test points along contour
    1003 % u=[0:0.1:1];
    1004 % x12=u*corners_X(2)+(1-u)*corners_X(1);
    1005 % x23=u*corners_X(3)+(1-u)*corners_X(2);
    1006 % x34=u*corners_X(4)+(1-u)*corners_X(3);
    1007 % x41=u*corners_X(1)+(1-u)*corners_X(4);
    1008 % y12=u*corners_Y(2)+(1-u)*corners_Y(1);
    1009 % y23=u*corners_Y(3)+(1-u)*corners_Y(2);
    1010 % y34=u*corners_Y(4)+(1-u)*corners_Y(3);
    1011 % y41=u*corners_Y(1)+(1-u)*corners_Y(4);
    1012 % hold on
    1013 % plot(x12,y12,'xr')
    1014 % plot(x23,y23,'xr')
    1015 % plot(x34,y34,'xr')
    1016 % plot(x41,y41,'xr')
     1040
    10171041%calculate transform matrices:
    10181042% reference: http://alumni.media.mit.edu/~cwren/interpolator/ by Christopher R. Wren
    … …  
    10261050C = [l(7:8)' 1];
    10271051
    1028 GeometryCalib.CalibrationType='tsai';%'linear';
     1052GeometryCalib.CalibrationType='tsai';
    10291053GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function
    10301054GeometryCalib.f=1;
    … …  
    10391063GeometryCalib.Tz=1;
    10401064GeometryCalib.R=[Amat(1,1),Amat(1,2),0;Amat(2,1),Amat(2,2),0;C(1),C(2),0];
    1041 %montre points du pourtour transformes
    1042 % for ip=1:numel(u),
    1043 %   XY12(:,ip)=Amat*[x12(ip);y12(ip);1]/(C*[x12(ip);y12(ip);1]);
    1044 %   XY23(:,ip)=Amat*[x23(ip);y23(ip);1]/(C*[x23(ip);y23(ip);1]);
    1045 %   XY34(:,ip)=Amat*[x34(ip);y34(ip);1]/(C*[x34(ip);y34(ip);1]);
    1046 %   XY41(:,ip)=Amat*[x41(ip);y41(ip);1]/(C*[x41(ip);y41(ip);1]);
    1047 % end
    1048 % xphys12=u*X(2)+(1-u)*X(1);
    1049 % xphys23=u*X(3)+(1-u)*X(2);
    1050 % xphys34=u*X(4)+(1-u)*X(3);
    1051 % xphys41=u*X(1)+(1-u)*X(4);
    1052 % yphys12=u*Y(2)+(1-u)*Y(1);
    1053 % yphys23=u*Y(3)+(1-u)*Y(2);
    1054 % yphys34=u*Y(4)+(1-u)*Y(3);
    1055 % yphys41=u*Y(1)+(1-u)*Y(4);
    1056 % for ip=1:numel(u),
    1057 %   XY12(:,ip)=Amat*[xphys12(ip);yphys12(ip);1]/(C*[xphys12(ip);yphys12(ip);1]);
    1058 %   XY23(:,ip)=Amat*[xphys23(ip);yphys23(ip);1]/(C*[xphys23(ip);yphys23(ip);1]);
    1059 %   XY34(:,ip)=Amat*[xphys34(ip);yphys34(ip);1]/(C*[xphys34(ip);yphys34(ip);1]);
    1060 %   XY41(:,ip)=Amat*[xphys41(ip);yphys41(ip);1]/(C*[xphys41(ip);yphys41(ip);1]);
    1061 % end
    1062 % x12=XY12(1,:);
    1063 % y12=XY12(2,:);
    1064 % x23=XY23(1,:);
    1065 % y23=XY23(2,:);
    1066 % x34=XY34(1,:);
    1067 % y34=XY34(2,:);
    1068 % x41=XY41(1,:);
    1069 % y41=XY41(2,:);
    1070 % [x12,y12]=px_XYZ(GeometryCalib,xphys12,yphys12);
    1071 % [x23,y23]=px_XYZ(GeometryCalib,xphys23,yphys23);
    1072 % [x34,y34]=px_XYZ(GeometryCalib,xphys34,yphys34);
    1073 % [x41,y41]=px_XYZ(GeometryCalib,xphys41,yphys41);
    1074 % plot(x12,y12,'xb')
    1075 % plot(x23,y23,'xb')
    1076 % plot(x34,y34,'xb')
    1077 % plot(x41,y41,'xb')
    1078 % figure(2)
    1079 % hold off
    1080 % plot(XY12(1,:),XY12(2,:),'ob') 
    1081 % hold on
    1082 % plot(XY23(1,:),XY23(2,:),'ob')
    1083 % plot(XY34(1,:),XY34(2,:),'ob')
    1084 % plot(XY41(1,:),XY41(2,:),'ob')
    1085 % plot(xphys12,yphys12,'ob') 
    1086 % hold on
    1087 % plot(xphys23,yphys23,'ob')
    1088 % plot(xphys34,yphys34,'ob')
    1089 % plot(xphys41,yphys41,'ob')
    10901065
    10911066[Amod,Rangx,Rangy]=phys_Ima(A-min(min(A)),GeometryCalib,0);
    10921067
    10931068Amod=double(Amod);
    1094 figure(12)
    1095 Amax=max(max(Amod));
    1096 image(Rangx,Rangy,uint8(255*Amod/Amax))
     1069% figure(12) display corrected image
     1070% Amax=max(max(Amod));
     1071% image(Rangx,Rangy,uint8(255*Amod/Amax))
    10971072Dx=(Rangx(2)-Rangx(1))/(npxy(2)-1); %x mesh in real space
    10981073Dy=(Rangy(2)-Rangy(1))/(npxy(1)-1); %y mesh in real space
    … …  
    11031078    i0=1+round((T(ipoint,1)-Rangx(1))/Dx);%round(Xpx(ipoint));
    11041079    j0=1+round((T(ipoint,2)-Rangy(1))/Dy);%round(Xpx(ipoint));
    1105     Asub=Amod(j0-ind_range_y:j0+ind_range_y,i0-ind_range_x:i0+ind_range_x);
     1080    j0min=max(j0-ind_range_y,1);
     1081    j0max=min(j0+ind_range_y,size(Amod,1));
     1082    i0min=max(i0-ind_range_x,1);
     1083    i0max=min(i0+ind_range_x,size(Amod,2));
     1084    Asub=Amod(j0min:j0max,i0min:i0max);
    11061085    x_profile=sum(Asub,1);
    11071086    y_profile=sum(Asub,2);
    … …  
    11111090    x_shift=0;
    11121091    y_shift=0;
    1113     if ind_x_max+2<=2*ind_range_x+1 && ind_x_max-2>=1
     1092    %if ind_x_max+2<=2*ind_range_x+1 && ind_x_max-2>=1
     1093    if ind_x_max+2<=numel(x_profile) && ind_x_max-2>=1
    11141094        Atop=x_profile(ind_x_max-2:ind_x_max+2);
    11151095        x_shift=sum(Atop.*[-2 -1 0 1 2])/sum(Atop);
    11161096    end
    1117     if ind_y_max+2<=2*ind_range_y+1 && ind_y_max-2>=1
     1097    if ind_x_max+2<=numel(y_profile) && ind_y_max-2>=1
    11181098        Atop=y_profile(ind_y_max-2:ind_y_max+2);
    11191099        y_shift=sum(Atop.*[-2 -1 0 1 2]')/sum(Atop);
    11201100    end
    1121     Delta(ipoint,1)=(x_shift+ind_x_max-ind_range_x-1)*Dx;%shift from the initial guess
    1122     Delta(ipoint,2)=(y_shift+ind_y_max-ind_range_y-1)*Dy;
     1101    Delta(ipoint,1)=(x_shift+ind_x_max+i0min-i0-1)*Dx;%shift from the initial guess
     1102    Delta(ipoint,2)=(y_shift+ind_y_max+j0min-j0-1)*Dy;
    11231103end
    11241104Tmod=T(:,(1:2))+Delta;
    … …  
    11271107     Coord{ipoint,1}=num2str(T(ipoint,1),4);%display coordiantes with 4 digits
    11281108     Coord{ipoint,2}=num2str(T(ipoint,2),4);%display coordiantes with 4 digits
    1129      Coord{ipoint,3}='0';
    1130      Coord{ipoint,4}=num2str(Xpx(ipoint),4);%display coordiantes with 4 digi
    1131      Coord{ipoint,5}=num2str(Ypx(ipoint),4);%display coordiantes with 4 digi
     1109     Coord{ipoint,3}=num2str(T(ipoint,3),4);%display coordiantes with 4 digits;
     1110     Coord{ipoint,4}=num2str(Xpx(ipoint),4);%display coordiantes with 4 digits
     1111     Coord{ipoint,5}=num2str(Ypx(ipoint),4);%display coordiantes with 4 digits
    11321112end
    11331113Tabchar=cell2tab(Coord(end:-1:1,:),'    |    ');
    … …  
    12091189set(handles.ListCoord,'Value',1)
    12101190set(handles.ListCoord,'String',Tabchar)
     1191
     1192
     1193% %------------------------------------------------------------------------
     1194% % --- Executes on button press in rotation.
     1195% function rotation_Callback(hObject, eventdata, handles)
     1196% %------------------------------------------------------------------------
     1197% angle_rot=(pi/180)*str2num(get(handles.Phi,'String'));
     1198% Coord_cell=get(handles.ListCoord,'String');
     1199% data=read_geometry_calib(Coord_cell);
     1200% data.Coord(:,1)=cos(angle_rot)*data.Coord(:,1)+sin(angle_rot)*data.Coord(:,2);
     1201% data.Coord(:,1)=-sin(angle_rot)*data.Coord(:,1)+cos(angle_rot)*data.Coord(:,2);
     1202% set(handles.XObject,'String',num2str(data.Coord(:,1),4));
     1203% set(handles.YObject,'String',num2str(data.Coord(:,2),4));
     1204
    12111205
    12121206%------------------------------------------------------------------------
    … …  
    12871281    Zphys=Calib.SliceCoord(Zindex,3);%GENERALISER AUX CAS AVEC ANGLE
    12881282else
    1289 %     if exist('Z','var')
    1290 %         Zphys=Z;
    1291 %     else
    1292         Zphys=0;
    1293 %     end
     1283    Zphys=0;
    12941284end
    12951285if ~exist('X','var')||~exist('Y','var')
    … …  
    13291319    Yphys=(A21.*Xu+A22.*Yu+Y0)./denom;
    13301320end
     1321
     1322
     1323% --------------------------------------------------------------------
     1324function MenuImportPoints_Callback(hObject, eventdata, handles)
     1325fileinput=browse_xml(hObject, eventdata, handles)
     1326if isempty(fileinput)
     1327    return
     1328end
     1329[s,errormsg]=imadoc2struct(fileinput,'GeometryCalib');
     1330GeometryCalib=s.GeometryCalib;
     1331GeometryCalib=load_calib(hObject, eventdata, handles)
     1332calib=reshape(GeometryCalib.PointCoord',[],1);
     1333for ilist=1:numel(calib)
     1334    CoordCell{ilist}=num2str(calib(ilist));
     1335end
     1336CoordCell=reshape(CoordCell,[],5);
     1337Tabchar=cell2tab(CoordCell,'    |    ');%transform cells into table ready for display
     1338Tabchar=[Tabchar;{'......'}];
     1339set(handles.ListCoord,'Value',1)
     1340set(handles.ListCoord,'String',Tabchar)
     1341MenuPlot_Callback(handles.geometry_calib, [], handles)
     1342
     1343% -----------------------------------------------------------------------
     1344function MenuImportIntrinsic_Callback(hObject, eventdata, handles)
     1345%------------------------------------------------------------------------
     1346fileinput=browse_xml(hObject, eventdata, handles);
     1347if isempty(fileinput)
     1348    return
     1349end
     1350[s,errormsg]=imadoc2struct(fileinput,'GeometryCalib');
     1351GeometryCalib=s.GeometryCalib;
     1352k=GeometryCalib.kappa1 * GeometryCalib.f*GeometryCalib.f;
     1353f1=GeometryCalib.f*GeometryCalib.sx/GeometryCalib.dpx;
     1354f2=GeometryCalib.f/GeometryCalib.dpy;
     1355Cx=GeometryCalib.Cx;
     1356Cy=GeometryCalib.Cy;
     1357set(handles.fx,'String',num2str(f1,'%1.1f'))
     1358set(handles.fy,'String',num2str(f2,'%1.1f'))
     1359set(handles.k,'String',num2str(k,'%1.4f'))
     1360set(handles.Cx,'String',num2str(Cx,'%1.1f'))
     1361set(handles.Cy,'String',num2str(Cy,'%1.1f'))
     1362
     1363% -----------------------------------------------------------------------
     1364function MenuImportAll_Callback(hObject, eventdata, handles)
     1365%------------------------------------------------------------------------
     1366fileinput=browse_xml(hObject, eventdata, handles)
     1367if ~isempty(fileinput)
     1368    loadfile(handles,fileinput)
     1369end
     1370
     1371% -----------------------------------------------------------------------
     1372% --- Executes on menubar option Import/Grid file: introduce previous grid files
     1373function MenuGridFile_Callback(hObject, eventdata, handles)
     1374% -----------------------------------------------------------------------
     1375inputfile=browse_xml(hObject, eventdata, handles)
     1376listfile=get(handles.coord_files,'string');
     1377if isequal(listfile,{''})
     1378    listfile={inputfile};
     1379else
     1380    listfile=[listfile;{inputfile}]%update the list of coord files
     1381end
     1382set(handles.coord_files,'string',listfile);
     1383
     1384%------------------------------------------------------------------------
     1385function fileinput=browse_xml(hObject, eventdata, handles)
     1386%------------------------------------------------------------------------
     1387fileinput=[];%default
     1388oldfile=''; %default
     1389UserData=get(handles.geometry_calib,'UserData');
     1390if isfield(UserData,'XmlInput')
     1391    oldfile=UserData.XmlInput;
     1392end
     1393[FileName, PathName, filterindex] = uigetfile( ...
     1394       {'*.xml;*.mat', ' (*.xml,*.mat)';
     1395       '*.xml',  '.xml files '; ...
     1396        '*.mat',  '.mat matlab files '}, ...
     1397        'Pick a file',oldfile);
     1398fileinput=[PathName FileName];%complete file name
     1399testblank=findstr(fileinput,' ');%look for blanks
     1400if ~isempty(testblank)
     1401    msgbox_uvmat('ERROR','forbidden input file name or path: no blank character allowed')
     1402    return
     1403end
     1404sizf=size(fileinput);
     1405if (~ischar(fileinput)||~isequal(sizf(1),1)),return;end
     1406UserData.XmlInput=fileinput;
     1407set(handles.geometry_calib,'UserData',UserData)%record current file foer further use of browser
     1408
     1409% -----------------------------------------------------------------------
     1410function loadfile(handles,fileinput)
     1411%------------------------------------------------------------------------
     1412[s,errormsg]=imadoc2struct(fileinput,'GeometryCalib');
     1413GeometryCalib=s.GeometryCalib;
     1414if 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;
     1423else
     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'))
     1434    calib=reshape(GeometryCalib.PointCoord',[],1);
     1435    for ilist=1:numel(calib)
     1436        CoordCell{ilist}=num2str(calib(ilist));
     1437    end
     1438    CoordCell=reshape(CoordCell,[],5);
     1439    Tabchar=cell2tab(CoordCell,'    |    ');%transform cells into table ready for display
     1440    set(handles.ListCoord,'Value',1)
     1441    set(handles.ListCoord,'String',Tabchar)
     1442    MenuPlot_Callback(handles.geometry_calib, [], handles)
     1443end
     1444Tabchar=[Tabchar;{'......'}];
     1445if isempty(CoordCell)% allow mouse action by default in the absence of input points
     1446    set(handles.edit_append,'Value',1)
     1447    set(handles.edit_append,'BackgroundColor',[1 1 0])
     1448else % does not allow mouse action by default in the presence of input points
     1449    set(handles.edit_append,'Value',0)
     1450    set(handles.edit_append,'BackgroundColor',[0.7 0.7 0.7])
     1451end
     1452
     1453
     1454
     1455
     1456
     1457% --------------------------------------------------------------------
     1458% --- Executes on button press in CLEAR_PTS: clear the list of calibration points
     1459function CLEAR_PTS_Callback(hObject, eventdata, handles)
     1460% --------------------------------------------------------------------
     1461set(handles.ListCoord,'Value',1)% refresh the display of coordinates
     1462set(handles.ListCoord,'String',{'......'})
     1463
Note: See TracChangeset for help on using the changeset viewer.