# source:trunk/src/geometry_calib.m@115

Last change on this file since 115 was 114, checked in by sommeria, 14 years ago

geometric calibration procedures updated

File size: 59.7 KB
Line
1%'geometry_calib': performs geometric calibration from a set of reference points
2%
3% function varargout = geometry_calib(varargin)
4%
5%A%AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
6%  Copyright Joel Sommeria, 2008, LEGI / CNRS-UJF-INPG, sommeria@coriolis-legi.org.
7%AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
8%     This file is part of the toolbox UVMAT.
9%
10%     UVMAT is free software; you can redistribute it and/or modify
11%     it under the terms of the GNU General Public License as published by
12%     the Free Software Foundation; either version 2 of the License, or
13%     (at your option) any later version.
14%
15%     UVMAT is distributed in the hope that it will be useful,
16%     but WITHOUT ANY WARRANTY; without even the implied warranty of
17%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
18%     GNU General Public License (file UVMAT/COPYING.txt) for more details.
19%AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
20
21function varargout = geometry_calib(varargin)
22% GEOMETRY_CALIB M-file for geometry_calib.fig
23%      GEOMETRY_CALIB, by itself, creates a MenuCoord GEOMETRY_CALIB or raises the existing
24%      singleton*.
25%
26%      H = GEOMETRY_CALIB returns the handle to a MenuCoord GEOMETRY_CALIB or the handle to
27%      the existing singleton*.
28%
29%      GEOMETRY_CALIB('CALLBACK',hObject,eventData,handles,...) calls the local
30%      function named CALLBACK in GEOMETRY_CALIB.M with the given input arguments.
31%
32%      GEOMETRY_CALIB('Property','Value',...) creates a MenuCoord GEOMETRY_CALIB or raises the
33%      existing singleton*.  Starting from the left, property value pairs are
34%      applied to the GUI before geometry_calib_OpeningFunction gets called.  An
35%      unrecognized property name or invalid value makes property application
36%      stop.  All inputs are passed to geometry_calib_OpeningFcn via varargin.
37%
38%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
39%      instance to run (singleton)".
40%
41% See also: GUIDE, GUIDATA, GUIHANDLES
42
43% Edit the above text to modify the response to help geometry_calib
44
45% Last Modified by GUIDE v2.5 05-Oct-2010 13:47:00
46
47% Begin initialization code - DO NOT edit
48gui_Singleton = 1;
49gui_State = struct('gui_Name',       mfilename, ...
50                   'gui_Singleton',  gui_Singleton, ...
51                   'gui_OpeningFcn', @geometry_calib_OpeningFcn, ...
52                   'gui_OutputFcn',  @geometry_calib_OutputFcn, ...
53                   'gui_LayoutFcn',  [] , ...
54                   'gui_Callback',   []);
55if nargin & isstr(varargin{1})
56    gui_State.gui_Callback = str2func(varargin{1});
57end
58
59if nargout
60    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
61else
62    gui_mainfcn(gui_State, varargin{:});
63end
64% End initialization code - DO NOT edit
65
66
67% --- Executes just before geometry_calib is made visible.
68%INPUT:
69%handles: handles of the geometry_calib interface elements
70% PlotHandles: set of handles of the elements contolling the plotting
71% parameters on the uvmat interface (obtained by 'get_plot_handle.m')
72%------------------------------------------------------------------------
73function geometry_calib_OpeningFcn(hObject, eventdata, handles, handles_uvmat,pos,inputfile)
74%------------------------------------------------------------------------
75% Choose default command line output for geometry_calib
76handles.output = hObject;
77
78% Update handles structure
79guidata(hObject, handles);
80set(hObject,'DeleteFcn',{@closefcn})%
81
82%set the position of the interface
83if exist('pos','var')& length(pos)>2
84    pos_gui=get(hObject,'Position');
85    pos_gui(1)=pos(1);
86    pos_gui(2)=pos(2);
87    set(hObject,'Position',pos_gui);
88end
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'})
93inputxml='';
94if exist('inputfile','var')& ~isempty(inputfile)
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'];
99end
100set(handles.ListCoord,'String',{'......'})
101if exist(inputxml,'file')
102    loadfile(handles,inputxml)% load the point coordiantes existing in the xml file
103end
104set(handles.ListCoord,'KeyPressFcn',{@key_press_fcn,handles})%set keyboard action function
105
106
107%------------------------------------------------------------------------
108% --- Outputs from this function are returned to the command line.
109function varargout = geometry_calib_OutputFcn(hObject, eventdata, handles)
110%------------------------------------------------------------------------
111% Get default command line output from handles structure
112varargout{1} = handles.output;
113varargout{2}=handles;
114%
115%------------------------------------------------------------------------
116% executed when closing: set the parent interface button to value 0
117function closefcn(gcbo,eventdata)
118%------------------------------------------------------------------------
119huvmat=findobj(allchild(0),'Name','uvmat');
120if ~isempty(huvmat)
121    handles=guidata(huvmat);
122    set(handles.MenuMask,'enable','on')
123    set(handles.MenuGrid,'enable','on')
124    set(handles.MenuObject,'enable','on')
125    set(handles.MenuEdit,'enable','on')
126    set(handles.edit,'enable','on')
127    hobject=findobj(handles.axes3,'tag','calib_points');
128    if ~isempty(hobject)
129        delete(hobject)
130    end
131    hobject=findobj(handles.axes3,'tag','calib_marker');
132    if ~isempty(hobject)
133        delete(hobject)
134    end
135end
136
137%------------------------------------------------------------------------
138% --- Executes on button press in calibrate_lin.
139function APPLY_Callback(hObject, eventdata, handles)
140%------------------------------------------------------------------------
141
142%read the current calibration points
143Coord_cell=get(handles.ListCoord,'String');
144Object=read_geometry_calib(Coord_cell);
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
156Z_plane=[];
157if ~isempty(Coord)
158    %check error
159    X=Coord(:,1);
160    Y=Coord(:,2);
161    Z=Coord(:,3);
162    x_ima=Coord(:,4);
163    y_ima=Coord(:,5);
164    [Xpoints,Ypoints]=px_XYZ(GeometryCalib,X,Y,Z);
165    GeometryCalib.ErrorRms(1)=sqrt(mean((Xpoints-x_ima).*(Xpoints-x_ima)));
166    [GeometryCalib.ErrorMax(1),index(1)]=max(abs(Xpoints-x_ima));
167    GeometryCalib.ErrorRms(2)=sqrt(mean((Ypoints-y_ima).*(Ypoints-y_ima)));
168    [GeometryCalib.ErrorMax(2),index(2)]=max(abs(Ypoints-y_ima));
169    [EM,ind_dim]=max(GeometryCalib.ErrorMax);
170    index=index(ind_dim);
171    %set the Z position of the reference plane used for calibration
172    if isequal(max(Z),min(Z))%Z constant
173        Z_plane=Z(1);
174        GeometryCalib.NbSlice=1;
175        GeometryCalib.SliceCoord=[0 0 Z_plane];
176    end
177end
178%set the coordinate unit
179unitlist=get(handles.CoordUnit,'String');
180unit=unitlist{get(handles.CoordUnit,'value')};
181GeometryCalib.CoordUnit=unit;
182%record the points
183GeometryCalib.SourceCalib.PointCoord=Coord;
184display_intrinsic(GeometryCalib,handles)%display calibration intrinsic parameters
185
186% Display extrinsinc parameters (rotation and translation of camera with  respect to the phys coordiantes)
187set(handles.Tx,'String',num2str(GeometryCalib.Tx_Ty_Tz(1),4))
188set(handles.Ty,'String',num2str(GeometryCalib.Tx_Ty_Tz(2),4))
189set(handles.Tz,'String',num2str(GeometryCalib.Tx_Ty_Tz(3),4))
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
196RootPath='';
197RootFile='';
198if ~isempty(hhuvmat.RootPath)& ~isempty(hhuvmat.RootFile)
199    testhandle=1;
200    RootPath=get(hhuvmat.RootPath,'String');
201    RootFile=get(hhuvmat.RootFile,'String');
202    filebase=fullfile(RootPath,RootFile);
203    outputfile=[filebase '.xml'];%xml file associated with the currently displayed image
204else
205    question={'save the calibration data and point coordinates in'};
206    def={fullfile(RootPath,['ObjectCalib.xml'])};
207    options.Resize='on';
208    answer=inputdlg(question,'save average in a new file',1,def,options);
209    outputfile=answer{1};
210end
211answer=msgbox_uvmat('INPUT_Y-N',{[outputfile ' updated with calibration data'];...
212    ['Error rms (along x,y)=' num2str(GeometryCalib.ErrorRms) ' pixels'];...
213    ['Error max (along x,y)=' num2str(GeometryCalib.ErrorMax) ' pixels']});
214if isequal(answer,'Yes')
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
223    %display image with new calibration in the currently opened uvmat interface
224    hhh=findobj(hhuvmat.axes3,'Tag','calib_marker');% delete calib points and markers
225    if ~isempty(hhh)
226        delete(hhh);
227    end
228    hhh=findobj(hhuvmat.axes3,'Tag','calib_points');
229    if ~isempty(hhh)
230        delete(hhh);
231    end
232    set(hhuvmat.FixedLimits,'Value',0)% put FixedLimits option to 'off'
233    set(hhuvmat.FixedLimits,'BackgroundColor',[0.7 0.7 0.7])
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)
238    figure(handles.geometry_calib)
239end
240
241%------------------------------------------------------------------
242% --- Executes on button press in calibrate_lin.
243function REPLICATE_Callback(hObject, eventdata, handles)
244% %%%%%%Todo: correct on the model of APPLY_Callback%%%%%%%%%%
245%------------------------------------------------------------------------
246calib_cell=get(handles.calib_type,'String');
247val=get(handles.calib_type,'Value');
248Coord_cell=get(handles.ListCoord,'String');
249Object=read_geometry_calib(Coord_cell);
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
329
330%open and read the dataview GUI
331h_dataview=findobj(allchild(0),'name','dataview');
332if ~isempty(h_dataview)
333    delete(h_dataview)
334end
335CalibData=get(handles.geometry_calib,'UserData');%read the calibration image source on the interface userdata
336
337if isfield(CalibData,'XmlInputFile')
338    XmlInput=fileparts(CalibData.XmlInputFile);
339    [XmlInput,filename,ext]=fileparts(XmlInput);
340end
341SubCampaignTest='n'; %default
342testinput=0;
343if isfield(CalibData,'Heading')
344    Heading=CalibData.Heading;
345    if isfield(Heading,'Record') && isequal([filename ext],Heading.Record)
346        [XmlInput,filename,ext]=fileparts(XmlInput);
347    end
348    if isfield(Heading,'Device') && isequal([filename ext],Heading.Device)
349        [XmlInput,filename,ext]=fileparts(XmlInput);
350        Device=Heading.Device;
351    end
352    if isfield(Heading,'Experiment') && isequal([filename ext],Heading.Experiment)
353        [PP,filename,ext]=fileparts(XmlInput);
354    end
355    testinput=0;
356    if isfield(Heading,'SubCampaign') && isequal([filename ext],Heading.SubCampaign)
357        SubCampaignTest='y';
358        testinput=1;
359    elseif isfield(Heading,'Campaign') && isequal([filename ext],Heading.Campaign)
360        testinput=1;
361    end
362end
363if ~testinput
364    filename='PROJETS';%default
365    if isfield(CalibData,'XmlInputFile')
366         [pp,filename]=fileparts(CalibData.XmlInputFile);
367    end
368    while ~isequal(filename,'PROJETS') && numel(filename)>1
369        filename_1=filename;
370        pp_1=pp;
371        [pp,filename]=fileparts(pp);
372    end
373    XmlInput=fullfile(pp_1,filename_1);
374    testinput=1;
375end
376if testinput
377    outcome=dataview(XmlInput,SubCampaignTest,GeometryCalib);
378end
379
380%------------------------------------------------------------------------
381% determine the parameters for a calibration by an affine function (rescaling and offset, no rotation)
382function GeometryCalib=calib_rescale(Coord,handles)
383%------------------------------------------------------------------------
384X=Coord(:,1);
385Y=Coord(:,2);% Z not used
386x_ima=Coord(:,4);
387y_ima=Coord(:,5);
388[px,sx]=polyfit(X,x_ima,1);
389[py,sy]=polyfit(Y,y_ima,1);
390T_x=px(2);
391T_y=py(2);
392GeometryCalib.CalibrationType='rescale';
393GeometryCalib.fx_fy=[px(1) py(1)];
394GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function
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];
398
399%------------------------------------------------------------------------
400% determine the parameters for a calibration by a linear transform matrix (rescale and rotation)
401function GeometryCalib=calib_linear(Coord,handles) %TO UPDATE
402%------------------------------------------------------------------------
403X=Coord(:,1);
404Y=Coord(:,2);% Z not used
405x_ima=Coord(:,4);
406y_ima=Coord(:,5);
407XY_mat=[ones(size(X)) X Y];
408a_X1=XY_mat\x_ima; %transformation matrix for X
409% x1=XY_mat*a_X1;%reconstruction
410% err_X1=max(abs(x1-x_ima));%error
411a_Y1=XY_mat\y_ima;%transformation matrix for X
412% y1=XY_mat*a_Y1;
413% err_Y1=max(abs(y1-y_ima));%error
414GeometryCalib.CalibrationType='linear';
415GeometryCalib.focal=1;
416GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function
417R=[a_X1(2),a_X1(3),0;a_Y1(2),a_Y1(3),0;0,0,0];
418norm=det(R);
419GeometryCalib.Tx_Ty_Tz=[a_X1(1) a_Y1(1) norm];
420GeometryCalib.R=R/norm;
421
422%------------------------------------------------------------------------
423% determine the tsai parameters for a view normal to the grid plane
424% NOT USED
425function GeometryCalib=calib_normal(Coord,handles)
426%------------------------------------------------------------------------
427Calib.f1=str2num(get(handles.fx,'String'));
428Calib.f2=str2num(get(handles.fy,'String'));
429Calib.k=str2num(get(handles.kc,'String'));
430Calib.Cx=str2num(get(handles.Cx,'String'));
431Calib.Cy=str2num(get(handles.Cy,'String'));
432%default
433if isempty(Calib.f1)
434    Calib.f1=25/0.012;
435end
436if isempty(Calib.f2)
437    Calib.f2=25/0.012;
438end
439if isempty(Calib.k)
440    Calib.k=0;
441end
442if isempty(Calib.Cx)||isempty(Calib.Cy)
443    huvmat=findobj(allchild(0),'Tag','uvmat');
444    hhuvmat=guidata(huvmat);
445    Calib.Cx=str2num(get(hhuvmat.npx,'String'))/2;
446    Calib.Cx=str2num(get(hhuvmat.npy,'String'))/2;
447end
448%tsai parameters
449Calib.dpx=0.012;%arbitrary
450Calib.dpy=0.012;
451Calib.sx=Calib.f1*Calib.dpx/(Calib.f2*Calib.dpy);
452Calib.f=Calib.f2*Calib.dpy;
453Calib.kappa1=Calib.k/(Calib.f*Calib.f);
454
455%initial guess
456X=Coord(:,1);
457Y=Coord(:,2);
458Zmean=mean(Coord(:,3));
459x_ima=Coord(:,4)-Calib.Cx;
460y_ima=Coord(:,5)-Calib.Cy;
461XY_mat=[ones(size(X)) X Y];
462a_X1=XY_mat\x_ima; %transformation matrix for X
463a_Y1=XY_mat\y_ima;%transformation matrix for Y
464R=[a_X1(2),a_X1(3),0;a_Y1(2),a_Y1(3),0;0,0,-1];% rotation+ z axis reversal (upward)
465norm=sqrt(det(-R));
466calib_param(1)=0;% quadratic distortion
467calib_param(2)=a_X1(1);
468calib_param(3)=a_Y1(1);
469calib_param(4)=Calib.f/(norm*Calib.dpx)-R(3,3)*Zmean;
470calib_param(5)=angle(a_X1(2)+i*a_X1(3));
471display(['initial guess=' num2str(calib_param)])
472
473%optimise the parameters: minimisation of error
474calib_param = fminsearch(@(calib_param) error_calib(calib_param,Calib,Coord),calib_param)
475
476GeometryCalib.CalibrationType='tsai_normal';
477GeometryCalib.focal=Calib.f;
478GeometryCalib.dpx_dpy=[Calib.dpx Calib.dpy];
479GeometryCalib.Cx_Cy=[Calib.Cx Calib.Cy];
480GeometryCalib.sx=Calib.sx;
481GeometryCalib.kappa1=calib_param(1);
482GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function
483GeometryCalib.Tx_Ty_Tz=[calib_param(2) calib_param(3) calib_param(4)];
484alpha=calib_param(5);
485GeometryCalib.R=[cos(alpha) sin(alpha) 0;-sin(alpha) cos(alpha) 0;0 0 -1];
486
487%------------------------------------------------------------------------
488function GeometryCalib=calib_3D_linear(Coord,handles)
489%TO UPDATE
490%------------------------------------------------------------------
491path_uvmat=which('uvmat');% check the path detected for source file uvmat
492path_UVMAT=fileparts(path_uvmat); %path to UVMAT
493huvmat=findobj(allchild(0),'Tag','uvmat');
494hhuvmat=guidata(huvmat);
495x_1=Coord(:,4:5)'
496X_1=Coord(:,1:3)'
497n_ima=1;
498% check_cond=0;
499
500nx=str2num(get(hhuvmat.npx,'String'));
501ny=str2num(get(hhuvmat.npy,'String'));
502
503est_dist=[0;0;0;0;0];
504est_aspect_ratio=0;
505est_fc=[1;1];
506%fc=[25;25]/0.012;
507center_optim=0;
508run(fullfile(path_UVMAT,'toolbox_calib','go_calib_optim'));
509
510GeometryCalib.CalibrationType='3D_linear';
511GeometryCalib.focal=fc(2);
512GeometryCalib.dpx_dpy=[1 1];
513GeometryCalib.Cx_Cy=cc';
514GeometryCalib.sx=fc(1)/fc(2);
515GeometryCalib.kappa1=-kc(1)/fc(2)^2;
516GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function
517GeometryCalib.Tx_Ty_Tz=Tc_1';
518GeometryCalib.R=Rc_1;
519
520%------------------------------------------------------------------------
521function GeometryCalib=calib_3D_quadr(Coord,handles)
522%------------------------------------------------------------------
523
524path_uvmat=which('uvmat');% check the path detected for source file uvmat
525path_UVMAT=fileparts(path_uvmat); %path to UVMAT
526huvmat=findobj(allchild(0),'Tag','uvmat');
527hhuvmat=guidata(huvmat);
528% check_cond=0;
529coord_files=get(handles.coord_files,'String');
530if ischar(coord_files)
531    coord_files={coord_files};
532end
533if isempty(coord_files{1}) || isequal(coord_files,{''})
534    coord_files={};
535end
536
537%retrieve the calibration points stored in the files listed in the popup list coord_files
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
543n_ima=numel(coord_files)+1;
544if ~isempty(coord_files)
545    msgbox_uvmat('CONFIRMATION',['The xy coordinates of the calibration points in ' num2str(n_ima) ' planes will be used'])
546    for ifile=1:numel(coord_files)
547    t=xmltree(coord_files{ifile});
548    s=convert(t);%convert to matlab structure
549        if isfield(s,'GeometryCalib')
550            if isfield(s.GeometryCalib,'SourceCalib')
551                if isfield(s.GeometryCalib.SourceCalib,'PointCoord')
552                PointCoord=s.GeometryCalib.SourceCalib.PointCoord;
553                Coord_file=[];
554                for i=1:length(PointCoord)
555                    line=str2num(PointCoord{i});
556                    Coord_file(i,4:5)=line(4:5);%px x
557                    Coord_file(i,1:3)=line(1:3);%phys x
558                end
559                eval(['x_' num2str(ifile+1) '=Coord_file(:,4:5)'';']);
560                eval(['x_' num2str(ifile+1) '(2,:)=ny-x_' num2str(ifile+1) '(2,:);' ]);
561                eval(['X_' num2str(ifile+1) '=Coord_file(:,1:3)'';']);
562                end
563            end
564        end
565    end
566end
567
568n_ima=numel(coord_files)+1;
569whos
570
571est_dist=[1;0;0;0;0];
572est_aspect_ratio=1;
573%est_fc=[0;0];
574%fc=[25;25]/0.012;
575center_optim=0;
576run(fullfile(path_UVMAT,'toolbox_calib','go_calib_optim'));
577
578GeometryCalib.CalibrationType='3D_quadr';
579GeometryCalib.fx_fy=fc';
580%GeometryCalib.focal=fc(2);
581%GeometryCalib.dpx_dpy=[1 1];
582GeometryCalib.Cx_Cy=cc';
583%GeometryCalib.sx=fc(1)/fc(2);
584GeometryCalib.kc=kc(1);
585%GeometryCalib.kappa1=-kc(1)/fc(2)^2;
586GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function
587GeometryCalib.Tx_Ty_Tz=Tc_1';
588GeometryCalib.R=Rc_1;
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
593GeometryCalib.ErrorRMS=[];
594GeometryCalib.ErrorMax=[];
595
596
597%------------------------------------------------------------------------
598function GeometryCalib=calib_3D_extrinsic(Coord,handles)
599%------------------------------------------------------------------
600path_uvmat=which('geometry_calib');% check the path detected for source file uvmat
601path_UVMAT=fileparts(path_uvmat); %path to UVMAT
602x_1=Coord(:,4:5)';
603X_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
608n_ima=1;
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'));
615fct_path=fullfile(path_UVMAT,'toolbox_calib');
616addpath(fct_path)
617Cx_Cy=ny-(GeometryCalib.Cx_Cy)';%reverse Cx_Cy(2) for calibration (inversion of px ordinate)
618% [omc1,Tc1,Rc1,H,x,ex,JJ] = compute_extrinsic(x_1,X_1,...
619%     [Calib.f Calib.f*Calib.sx]',...
620%     [Calib.Cx Calib.Cy]',...
621%     [-Calib.kappa1*Calib.f^2 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]);
625%get the euler angles ???
626rmpath(fct_path);
627
628std(ex')
629GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function
630GeometryCalib.Tx_Ty_Tz=Tc1';
631%inversion of z axis
632GeometryCalib.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';
637%GeometryCalib.R(3,1:3)=-GeometryCalib.R(3,1:3);%inversion for z upward
638
639
640%------------------------------------------------------------------------
641%function GeometryCalib=calib_tsai_heikkila(Coord)
642% TEST: NOT IMPLEMENTED
643%------------------------------------------------------------------
644% path_uvmat=which('uvmat');% check the path detected for source file uvmat
645% path_UVMAT=fileparts(path_uvmat); %path to UVMAT
646% path_calib=fullfile(path_UVMAT,'toolbox_calib_heikkila');
647% addpath(path_calib)
648% npoints=size(Coord,1);
649% Coord(:,1:3)=10*Coord(:,1:3);
650% Coord=[Coord zeros(npoints,2) -ones(npoints,1)];
651% [par,pos,iter,res,er,C]=cacal('dalsa',Coord);
652% GeometryCalib.CalibrationType='tsai';
653% GeometryCalib.focal=par(2);
654
655
656%--------------------------------------------------------------------------
657function GeometryCalib=calib_tsai(Coord,handles)% old version using gauthier's bianry ccal_fo
658% NOT USED
659%------------------------------------------------------------------------
660%TSAI
661path_uvmat=which('uvmat');% check the path detected for source file uvmat
662path_UVMAT=fileparts(path_uvmat); %path to UVMAT
663xmlfile=fullfile(path_UVMAT,'PARAM.xml');%name of the file containing names of binary executables
664if exist(xmlfile,'file')
665    t=xmltree(xmlfile);% read the (xml) file containing names of binary executables
666    sparam=convert(t);% convert to matlab structure
667end
668if ~isfield(sparam,'GeometryCalibBin')
669    msgbox_uvmat('ERROR',['calibration program <GeometryCalibBin> undefined in parameter file ' xmlfile])
670    return
671end
672Tsai_exe=sparam.GeometryCalibBin;
673if ~exist(Tsai_exe,'file')%the binary is defined in /bin, default setting
674     Tsai_exe=fullfile(path_UVMAT,Tsai_exe);
675end
676if ~exist(Tsai_exe,'file')
677    msgbox_uvmat('ERROR',['calibration program ' sparam.GeometryCalibBin ' defined in PARAM.xml does not exist'])
678    return
679end
680
681textcoord=num2str(Coord,4);
682dlmwrite('t.txt',textcoord,'');
683% ['!' Tsai_exe ' -fx 0 -fy t.txt']
684eval(['!' Tsai_exe ' -f t.txt > tsaicalib.log']);
685if ~exist('calib.dat','file')
686    msgbox_uvmat('ERROR','no output from calibration program Tsai_exe: possibly too few points')
687end
688calibdat=dlmread('calib.dat');
689delete('calib.dat')
690%delete('t.txt')
691GeometryCalib.CalibrationType='tsai';
692GeometryCalib.focal=calibdat(10);
693GeometryCalib.dpx_dpy=[calibdat(5) calibdat(6)];
694GeometryCalib.Cx_Cy=[calibdat(7) calibdat(8)];
695GeometryCalib.sx=calibdat(9);
696GeometryCalib.kappa1=calibdat(11);
697GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function
698GeometryCalib.Tx_Ty_Tz=[calibdat(12) calibdat(13) calibdat(14)];
699Rx_Ry_Rz=calibdat([15:17]);
700sa = sin(Rx_Ry_Rz(1)) ;
701ca=cos(Rx_Ry_Rz(1));
702sb=sin(Rx_Ry_Rz(2));
703cb =cos(Rx_Ry_Rz(2));
704sg =sin(Rx_Ry_Rz(3));
705cg =cos(Rx_Ry_Rz(3));
706r1 = cb * cg;
707r2 = cg * sa * sb - ca * sg;
708r3 = sa * sg + ca * cg * sb;
709r4 = cb * sg;
710r5 = sa * sb * sg + ca * cg;
711r6 = ca * sb * sg - cg * sa;
712r7 = -sb;
713r8 = cb * sa;
714r9 = ca * cb;
715%EN DEDUIRE MATRICE R ??
716GeometryCalib.R=[r1,r2,r3;r4,r5,r6;r7,r8,r9];
717
718%------------------------------------------------------------------------
719% --- determine the rms of calibration error
720function ErrorRms=error_calib(calib_param,Calib,Coord)
721%calib_param: vector of free calibration parameters (to optimise)
722%Calib: structure of the given calibration parameters
723%Coord: list of phys coordinates (columns 1-3, and pixel coordinates (columns 4-5)
724Calib.f=25;
725Calib.dpx=0.012;
726Calib.dpy=0.012;
727Calib.sx=1;
728Calib.Cx=512;
729Calib.Cy=512;
730Calib.kappa1=calib_param(1);
731Calib.Tx=calib_param(2);
732Calib.Ty=calib_param(3);
733Calib.Tz=calib_param(4);
734alpha=calib_param(5);
735Calib.R=[cos(alpha) sin(alpha) 0;-sin(alpha) cos(alpha) 0;0 0 -1];
736
737X=Coord(:,1);
738Y=Coord(:,2);
739Z=Coord(:,3);
740x_ima=Coord(:,4);
741y_ima=Coord(:,5);
742[Xpoints,Ypoints]=px_XYZ(Calib,X,Y,Z);
743ErrorRms(1)=sqrt(mean((Xpoints-x_ima).*(Xpoints-x_ima)));
744ErrorRms(2)=sqrt(mean((Ypoints-y_ima).*(Ypoints-y_ima)));
745ErrorRms=mean(ErrorRms);
746
747%------------------------------------------------------------------------
748function XImage_Callback(hObject, eventdata, handles)
749%------------------------------------------------------------------------
750update_list(hObject, eventdata,handles)
751
752%------------------------------------------------------------------------
753function YImage_Callback(hObject, eventdata, handles)
754%------------------------------------------------------------------------
755update_list(hObject, eventdata,handles)
756
757%------------------------------------------------------------------------
758% --- Executes on button press in STORE.
759function STORE_Callback(hObject, eventdata, handles)
760Coord_cell=get(handles.ListCoord,'String');
761Object=read_geometry_calib(Coord_cell);
762unitlist=get(handles.CoordUnit,'String');
763unit=unitlist{get(handles.CoordUnit,'value')};
764GeometryCalib.CoordUnit=unit;
765GeometryCalib.SourceCalib.PointCoord=Object.Coord;
766huvmat=findobj(allchild(0),'Name','uvmat');
767hhuvmat=guidata(huvmat);%handles of elements in the GUI uvmat
768RootPath='';
769RootFile='';
770if ~isempty(hhuvmat.RootPath)& ~isempty(hhuvmat.RootFile)
771    testhandle=1;
772    RootPath=get(hhuvmat.RootPath,'String');
773    RootFile=get(hhuvmat.RootFile,'String');
774    filebase=fullfile(RootPath,RootFile);
775    while exist([filebase '.xml'],'file')
776        filebase=[filebase '~'];
777    end
778    outputfile=[filebase '.xml'];
779    errormsg=update_imadoc(GeometryCalib,outputfile);
780    if ~strcmp(errormsg,'')
781        msgbox_uvmat('ERROR',errormsg);
782    end
783    listfile=get(handles.coord_files,'string');
784    if isequal(listfile,{''})
785        listfile={outputfile};
786    else
787        listfile=[listfile;{outputfile}]%update the list of coord files
788    end
789    set(handles.coord_files,'string',listfile);
790end
791set(handles.ListCoord,'Value',1)% refresh the display of coordinates
792set(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)
801
802%------------------------------------------------------------------------
803% --- Executes on button press in CLEAR.
804function CLEAR_Callback(hObject, eventdata, handles)
805%------------------------------------------------------------------------
806set(handles.coord_files,'Value',1)
807set(handles.coord_files,'String',{''})
808
809%------------------------------------------------------------------------
810function XObject_Callback(hObject, eventdata, handles)
811%------------------------------------------------------------------------
812update_list(hObject, eventdata,handles)
813
814%------------------------------------------------------------------------
815function YObject_Callback(hObject, eventdata, handles)
816%------------------------------------------------------------------------
817update_list(hObject, eventdata,handles)
818
819%------------------------------------------------------------------------
820function ZObject_Callback(hObject, eventdata, handles)
821%------------------------------------------------------------------------
822update_list(hObject, eventdata,handles)
823
824%------------------------------------------------------------------------
825function update_list(hObject, eventdata, handles)
826%------------------------------------------------------------------------
827str4=get(handles.XImage,'String');
828str5=get(handles.YImage,'String');
829str1=get(handles.XObject,'String');
830tt=double(str1);
831str2=get(handles.YObject,'String');
832str3=get(handles.ZObject,'String');
833if ~isempty(str1) & ~isequal(double(str1),32) & (isempty(str3)|isequal(double(str3),32))
834    str3='0';%put z to 0 by default
835end
836strline=[str1 '    |    ' str2 '    |    ' str3 '    |    ' str4 '    |    ' str5];
837Coord=get(handles.ListCoord,'String');
838val=get(handles.ListCoord,'Value');
839Coord{val}=strline;
840set(handles.ListCoord,'String',Coord)
841%update the plot
842ListCoord_Callback(hObject, eventdata, handles)
843MenuPlot_Callback(hObject, eventdata, handles)
844
845%------------------------------------------------------------------------
846% --- Executes on selection change in ListCoord.
847function ListCoord_Callback(hObject, eventdata, handles)
848%------------------------------------------------------------------------
849huvmat=findobj(allchild(0),'Name','uvmat');%find the current uvmat interface handle
850hplot=findobj(huvmat,'Tag','axes3');%main plotting axis of uvmat
851hhh=findobj(hplot,'Tag','calib_marker');
852Coord_cell=get(handles.ListCoord,'String');
853val=get(handles.ListCoord,'Value');
854if numel(val)>1
855    return %no action if several lines have been selected
856end
857coord_str=Coord_cell{val};
858k=findstr('|',coord_str);
859if isempty(k)%last line '.....' selected
860    if ~isempty(hhh)
861        delete(hhh)%delete the circle marker
862    end
863    return
864end
865%fill the edit boxex
866set(handles.XObject,'String',coord_str(1:k(1)-5))
867set(handles.YObject,'String',coord_str(k(1)+5:k(2)-5))
868set(handles.ZObject,'String',coord_str(k(2)+5:k(3)-5))
869set(handles.XImage,'String',coord_str(k(3)+5:k(4)-5))
870set(handles.YImage,'String',coord_str(k(4)+5:end))
871h_menu_coord=findobj(huvmat,'Tag','transform_fct');
872menu=get(h_menu_coord,'String');
873choice=get(h_menu_coord,'Value');
874if iscell(menu)
875    option=menu{choice};
876else
877    option='px'; %default
878end
879if isequal(option,'phys')
880    XCoord=str2num(coord_str(1:k(1)-5));
881    YCoord=str2num(coord_str(k(1)+5:k(2)-5));
882elseif isequal(option,'px')|| isequal(option,'')
883    XCoord=str2num(coord_str(k(3)+5:k(4)-5));
884    YCoord=str2num(coord_str(k(4)+5:end));
885else
886    msgbox_uvmat('ERROR','the choice in menu_coord of uvmat must be px or phys ')
887end
888if isempty(XCoord)||isempty(YCoord)
889     if ~isempty(hhh)
890        delete(hhh)%delete the circle marker
891    end
892    return
893end
894xlim=get(hplot,'XLim');
895ylim=get(hplot,'YLim');
896ind_range=max(abs(xlim(2)-xlim(1)),abs(ylim(end)-ylim(1)))/20;%defines the size of the circle marker
897if isempty(hhh)
898    axes(hplot)
899    rectangle('Curvature',[1 1],...
900              'Position',[XCoord-ind_range/2 YCoord-ind_range/2 ind_range ind_range],'EdgeColor','m',...
901              'LineStyle','-','Tag','calib_marker');
902else
903    set(hhh,'Position',[XCoord-ind_range/2 YCoord-ind_range/2 ind_range ind_range])
904end
905
906%------------------------------------------------------------------------
907% --- Executes on selection change in edit_append.
908function edit_append_Callback(hObject, eventdata, handles)
909%------------------------------------------------------------------------
910choice=get(handles.edit_append,'Value');
911if choice
912    set(handles.edit_append,'BackgroundColor',[1 1 0])
913else
914    set(handles.edit_append,'BackgroundColor',[0.7 0.7 0.7])
915end
916
917function NEW_Callback(hObject, eventdata, handles)
918%A METTRE SOUS UN BOUTON
919huvmat=findobj(allchild(0),'Name','uvmat');
920hchild=get(huvmat,'children');
921hcoord=findobj(hchild,'Tag','menu_coord');
922coordtype=get(hcoord,'Value');
923haxes=findobj(hchild,'Tag','axes3');
924AxeData=get(haxes,'UserData');
925if ~isequal(hcoord,2)
926    set(hcoord,'Value',2)
927    huvmat=uvmat(AxeData);
928    'relancer uvmat';
929end
930if ~isfield(AxeData,'ZoomAxes')
931    msgbox_uvmat('ERROR','first draw a window around a grid marker')
932    return
933end
934XLim=get(AxeData.ZoomAxes,'XLim');
935YLim=get(AxeData.ZoomAxes,'YLim');
936np=size(AxeData.A);
937ind_sub_x=round(XLim);
938ind_sub_y=np(1)-round(YLim);
939Mfiltre=AxeData.A([ind_sub_y(2):ind_sub_y(1)] ,ind_sub_x,:);
940Mfiltre_norm=double(Mfiltre);
941Mfiltre_norm=Mfiltre_norm/sum(sum(Mfiltre_norm));
942Mfiltre_norm=100*(Mfiltre_norm-mean(mean(Mfiltre_norm)));
943Atype=class(AxeData.A);
944Data.NbDim=2;
945Data.A=filter2(Mfiltre_norm,double(AxeData.A));
946Data.A=feval(Atype,Data.A);
947Data.AName='image';
948Data.AX=AxeData.AX;
949Data.AY=AxeData.AY;
950Data.CoordType='px';
951plot_field(Data)
952
953%------------------------------------------------------------------------
954function MenuPlot_Callback(hObject, eventdata, handles)
955%------------------------------------------------------------------------
956huvmat=findobj(allchild(0),'Name','uvmat');%find the current uvmat interface handle
957UvData=get(huvmat,'UserData');%Data associated to the current uvmat interface
958hhuvmat=guidata(huvmat); %handles of GUI elements in uvmat
959hplot=findobj(huvmat,'Tag','axes3');%main plotting axis of uvmat
960h_menu_coord=findobj(huvmat,'Tag','transform_fct');
961menu=get(h_menu_coord,'String');
962choice=get(h_menu_coord,'Value');
963if iscell(menu)
964    option=menu{choice};
965else
966    option='px'; %default
967end
968Coord_cell=get(handles.ListCoord,'String');
969ObjectData=read_geometry_calib(Coord_cell);
970%ObjectData=read_geometry_calib(handles);%read the interface input parameters defining the object
971if ~isempty(ObjectData.Coord)
972    if isequal(option,'phys')
973        ObjectData.Coord=ObjectData.Coord(:,[1:3]);
974    elseif isequal(option,'px')||isequal(option,'')
975        ObjectData.Coord=ObjectData.Coord(:,[4:5]);
976    else
977        msgbox_uvmat('ERROR','the choice in menu_coord of uvmat must be '''', px or phys ')
978    end
979end
980axes(hhuvmat.axes3)
981hh=findobj('Tag','calib_points');
982if  ~isempty(ObjectData.Coord) && isempty(hh)
983    hh=line(ObjectData.Coord(:,1),ObjectData.Coord(:,2),'Color','m','Tag','calib_points','LineStyle','.','Marker','+');
984elseif isempty(ObjectData.Coord)%empty list of points, suppress the plot
985    delete(hh)
986else
987    set(hh,'XData',ObjectData.Coord(:,1))
988    set(hh,'YData',ObjectData.Coord(:,2))
989end
990pause(.1)
991figure(handles.geometry_calib)
992
993% --------------------------------------------------------------------
994function MenuHelp_Callback(hObject, eventdata, handles)
995path_to_uvmat=which ('uvmat');% check the path of uvmat
996pathelp=fileparts(path_to_uvmat);
997    helpfile=fullfile(pathelp,'uvmat_doc','uvmat_doc.html');
998if isempty(dir(helpfile)), msgbox_uvmat('ERROR','Please put the help file uvmat_doc.html in the sub-directory /uvmat_doc of the UVMAT package')
999else
1000   addpath (fullfile(pathelp,'uvmat_doc'))
1001   web([helpfile '#geometry_calib'])
1002end
1003
1004%------------------------------------------------------------------------
1005function MenuCreateGrid_Callback(hObject, eventdata, handles)
1006%------------------------------------------------------------------------
1007%hcalib=get(handles.calib_type,'parent');%handles of the GUI geometry_calib
1008CalibData=get(handles.geometry_calib,'UserData');
1009Tinput=[];%default
1010if isfield(CalibData,'grid')
1011    Tinput=CalibData.grid;
1012end
1013[T,CalibData.grid]=create_grid(Tinput);%display the GUI create_grid
1014set(handles.geometry_calib,'UserData',CalibData)
1015
1016%grid in phys space
1017Coord=get(handles.ListCoord,'String');
1018val=get(handles.ListCoord,'Value');
1019data=read_geometry_calib(Coord);
1020%nbpoints=size(data.Coord,1); %nbre of calibration points
1021data.Coord(val:val+size(T,1)-1,1:3)=T(end:-1:1,:);%update the existing list of phys coordinates from the GUI create_grid
1022% for i=1:nbpoints
1023%    for j=1:5
1024%           Coord{i,j}=num2str(data.Coord(i,j),4);%display coordiantes with 4 digits
1025%    end
1026% end
1027%update the phys coordinates starting from the selected point (down in the
1028Coord(end,:)=[]; %remove last string '.....'
1029for i=1:size(data.Coord,1)
1030    for j=1:5
1031          Coord{i,j}=num2str(data.Coord(i,j),4);%display coordiantes with 4 digits
1032    end
1033end
1034
1035%size(data.Coord,1)
1036Tabchar=cell2tab(Coord,'    |    ');
1037Tabchar=[Tabchar ;{'......'}];
1038set(handles.ListCoord,'String',Tabchar)
1039
1040% -----------------------------------------------------------------------
1041% --- automatic grid dectection from local maxima of the images
1042function MenuDetectGrid_Callback(hObject, eventdata, handles)
1043%------------------------------------------------------------------------
1044CalibData=get(handles.geometry_calib,'UserData');%get information stored on the GUI geometry_calib
1045grid_input=[];%default
1046if isfield(CalibData,'grid')
1047    grid_input=CalibData.grid;%retrieve the previously used grid
1048end
1049[T,CalibData.grid]=create_grid(grid_input,'detect_grid');%display the GUI create_grid, read the set of phys coordinates T
1050
1051set(handles.geometry_calib,'UserData',CalibData)%store the phys grid parameters for later use
1052
1053%read the four last point coordinates in pixels
1054Coord_cell=get(handles.ListCoord,'String');%read list of coordinates on geometry_calib
1055data=read_geometry_calib(Coord_cell);
1056nbpoints=size(data.Coord,1); %nbre of calibration points
1057if nbpoints~=4
1058    msgbox_uvmat('ERROR','four points must be selected by the mouse, beginning by the new x axis, to delimitate the phs grid area')
1059    return
1060end
1061corners_X=(data.Coord(end:-1:end-3,4)); %pixel absissa of the four corners
1062corners_Y=(data.Coord(end:-1:end-3,5));
1063
1064%reorder the last two points (the two first in the list) if needed
1065angles=angle((corners_X-corners_X(1))+i*(corners_Y-corners_Y(1)));
1066if abs(angles(4)-angles(2))>abs(angles(3)-angles(2))
1067      X_end=corners_X(4);
1068      Y_end=corners_Y(4);
1069      corners_X(4)=corners_X(3);
1070      corners_Y(4)=corners_Y(3);
1071      corners_X(3)=X_end;
1072      corners_Y(3)=Y_end;
1073end
1074
1075%read the current image, displayed in the GUI uvmat
1076huvmat=findobj(allchild(0),'Name','uvmat');
1077UvData=get(huvmat,'UserData');
1078A=UvData.Field.A;
1079npxy=size(A);
1080%linear transform on the current image
1081X=[CalibData.grid.x_0 CalibData.grid.x_1 CalibData.grid.x_0 CalibData.grid.x_1]';%corner absissa in the phys coordinates
1082Y=[CalibData.grid.y_0 CalibData.grid.y_0 CalibData.grid.y_1 CalibData.grid.y_1]';%corner ordinates in the phys coordinates
1083
1084%calculate transform matrices:
1085% reference: http://alumni.media.mit.edu/~cwren/interpolator/ by Christopher R. Wren
1086B = [ X Y ones(size(X)) zeros(4,3)        -X.*corners_X -Y.*corners_X ...
1087      zeros(4,3)        X Y ones(size(X)) -X.*corners_Y -Y.*corners_Y ];
1088B = reshape (B', 8 , 8 )';
1089D = [ corners_X , corners_Y ];
1090D = reshape (D', 8 , 1 );
1091l = inv(B' * B) * B' * D;
1092Amat = reshape([l(1:6)' 0 0 1 ],3,3)';
1093C = [l(7:8)' 1];
1094
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];
1113GeometryCalib.R=[Amat(1,1),Amat(1,2),0;Amat(2,1),Amat(2,2),0;C(1),C(2),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
1145
1146Amod=double(Amod);
1147% figure(12) %display corrected image
1148% Amax=max(max(Amod));
1149% image(Rangx,Rangy,uint8(255*Amod/Amax))
1150
1151Dx=(Rangx(2)-Rangx(1))/(npxy(2)-1); %x mesh in real space
1152Dy=(Rangy(2)-Rangy(1))/(npxy(1)-1); %y mesh in real space
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
1155nbpoints=size(T,1);
1156for ipoint=1:nbpoints
1157    i0=1+round((T(ipoint,1)-Rangx(1))/Dx);%round(Xpx(ipoint));
1158    j0=1+round((T(ipoint,2)-Rangy(1))/Dy);%round(Xpx(ipoint));
1159    j0min=max(j0-ind_range_y,1);
1160    j0max=min(j0+ind_range_y,size(Amod,1));
1161    i0min=max(i0-ind_range_x,1);
1162    i0max=min(i0+ind_range_x,size(Amod,2));
1163    Asub=Amod(j0min:j0max,i0min:i0max);
1164    x_profile=sum(Asub,1);
1165    y_profile=sum(Asub,2);
1166    [Amax,ind_x_max]=max(x_profile);
1167    [Amax,ind_y_max]=max(y_profile);
1168    %sub-pixel improvement using moments
1169    x_shift=0;
1170    y_shift=0;
1171    %if ind_x_max+2<=2*ind_range_x+1 && ind_x_max-2>=1
1172    if ind_x_max+2<=numel(x_profile) && ind_x_max-2>=1
1173        Atop=x_profile(ind_x_max-2:ind_x_max+2);
1174        x_shift=sum(Atop.*[-2 -1 0 1 2])/sum(Atop);
1175    end
1176    if ind_y_max+2<=numel(y_profile) && ind_y_max-2>=1
1177        Atop=y_profile(ind_y_max-2:ind_y_max+2);
1178        y_shift=sum(Atop.*[-2 -1 0 1 2]')/sum(Atop);
1179    end
1180    Delta(ipoint,1)=(x_shift+ind_x_max+i0min-i0-1)*Dx;%shift from the initial guess
1181    Delta(ipoint,2)=(y_shift+ind_y_max+j0min-j0-1)*Dy;
1182end
1183Tmod=T(:,(1:2))+Delta;
1184[Xpx,Ypx]=px_XYZ(GeometryCalib,Tmod(:,1),Tmod(:,2));
1185for ipoint=1:nbpoints
1186     Coord{ipoint,1}=num2str(T(ipoint,1),4);%display coordiantes with 4 digits
1187     Coord{ipoint,2}=num2str(T(ipoint,2),4);%display coordiantes with 4 digits
1188     Coord{ipoint,3}=num2str(T(ipoint,3),4);%display coordiantes with 4 digits;
1189     Coord{ipoint,4}=num2str(Xpx(ipoint),4);%display coordiantes with 4 digits
1190     Coord{ipoint,5}=num2str(Ypx(ipoint),4);%display coordiantes with 4 digits
1191end
1192Tabchar=cell2tab(Coord(end:-1:1,:),'    |    ');
1193Tabchar=[Tabchar ;{'......'}];
1194set(handles.ListCoord,'Value',1)
1195set(handles.ListCoord,'String',Tabchar)
1196MenuPlot_Callback(hObject, eventdata, handles)
1197
1198%-----------------------------------------------------------------------
1199function MenuTranslatePoints_Callback(hObject, eventdata, handles)
1200%-----------------------------------------------------------------------
1201%hcalib=get(handles.calib_type,'parent');%handles of the GUI geometry_calib
1202CalibData=get(handles.geometry_calib,'UserData');
1203Tinput=[];%default
1204if isfield(CalibData,'translate')
1205    Tinput=CalibData.translate;
1206end
1207T=translate_points(Tinput);%display translate_points GUI and get shift parameters
1208CalibData.translate=T;
1209set(handles.geometry_calib,'UserData',CalibData)
1210%translation
1211Coord_cell=get(handles.ListCoord,'String');
1212data=read_geometry_calib(Coord_cell);
1213data.Coord(:,1)=T(1)+data.Coord(:,1);
1214data.Coord(:,2)=T(2)+data.Coord(:,2);
1215data.Coord(:,3)=T(3)+data.Coord(:,3);
1216data.Coord(:,[4 5])=data.Coord(:,[4 5]);
1217for i=1:size(data.Coord,1)
1218    for j=1:5
1219          Coord{i,j}=num2str(data.Coord(i,j),4);%phys x,y,z
1220   end
1221end
1222Tabchar=cell2tab(Coord,'    |    ');
1223Tabchar=[Tabchar; {'.....'}];
1224%set(handles.ListCoord,'Value',1)
1225set(handles.ListCoord,'String',Tabchar)
1226
1227
1228% --------------------------------------------------------------------
1229function MenuRotatePoints_Callback(hObject, eventdata, handles)
1230%hcalib=get(handles.calib_type,'parent');%handles of the GUI geometry_calib
1231CalibData=get(handles.geometry_calib,'UserData');
1232Tinput=[];%default
1233if isfield(CalibData,'rotate')
1234    Tinput=CalibData.rotate;
1235end
1236T=rotate_points(Tinput);%display translate_points GUI and get shift parameters
1237CalibData.rotate=T;
1238set(handles.geometry_calib,'UserData',CalibData)
1239%-----------------------------------------------------
1240%rotation
1241Phi=T(1);
1242O_x=0;%default
1243O_y=0;%default
1244if numel(T)>=2
1245    O_x=T(2);%default
1246end
1247if numel(T)>=3
1248    O_y=T(3);%default
1249end
1250Coord_cell=get(handles.ListCoord,'String');
1251data=read_geometry_calib(Coord_cell);
1252r1=cos(pi*Phi/180);
1253r2=-sin(pi*Phi/180);
1254r3=sin(pi*Phi/180);
1255r4=cos(pi*Phi/180);
1256x=data.Coord(:,1)-O_x;
1257y=data.Coord(:,2)-O_y;
1258data.Coord(:,1)=r1*x+r2*y;
1259data.Coord(:,2)=r3*x+r4*y;
1260% data.Coord(:,[4 5])=data.Coord(:,[4 5]);
1261for i=1:size(data.Coord,1)
1262    for j=1:5
1263          Coord{i,j}=num2str(data.Coord(i,j),4);%phys x,y,z
1264   end
1265end
1266Tabchar=cell2tab(Coord,'    |    ');
1267Tabchar=[Tabchar;{'......'}];
1268set(handles.ListCoord,'Value',1)
1269set(handles.ListCoord,'String',Tabchar)
1270
1271
1272% %------------------------------------------------------------------------
1273% % --- Executes on button press in rotation.
1274% function rotation_Callback(hObject, eventdata, handles)
1275% %------------------------------------------------------------------------
1276% angle_rot=(pi/180)*str2num(get(handles.Phi,'String'));
1277% Coord_cell=get(handles.ListCoord,'String');
1278% data=read_geometry_calib(Coord_cell);
1279% data.Coord(:,1)=cos(angle_rot)*data.Coord(:,1)+sin(angle_rot)*data.Coord(:,2);
1280% data.Coord(:,1)=-sin(angle_rot)*data.Coord(:,1)+cos(angle_rot)*data.Coord(:,2);
1281% set(handles.XObject,'String',num2str(data.Coord(:,1),4));
1282% set(handles.YObject,'String',num2str(data.Coord(:,2),4));
1283
1284
1285%------------------------------------------------------------------------
1286% image transform from px to phys
1287%INPUT:
1288%Zindex: index of plane
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
1351
1352%------------------------------------------------------------------------
1353% pointwise transform from px to phys
1354%INPUT:
1355%Z: index of plane
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
1400
1401
1402% --------------------------------------------------------------------
1403function MenuImportPoints_Callback(hObject, eventdata, handles)
1404fileinput=browse_xml(hObject, eventdata, handles)
1405if isempty(fileinput)
1406    return
1407end
1408[s,errormsg]=imadoc2struct(fileinput,'GeometryCalib');
1409GeometryCalib=s.GeometryCalib;
1410%GeometryCalib=load_calib(hObject, eventdata, handles)
1411calib=reshape(GeometryCalib.PointCoord',[],1);
1412for ilist=1:numel(calib)
1413    CoordCell{ilist}=num2str(calib(ilist));
1414end
1415CoordCell=reshape(CoordCell,[],5);
1416Tabchar=cell2tab(CoordCell,'    |    ');%transform cells into table ready for display
1417Tabchar=[Tabchar;{'......'}];
1418set(handles.ListCoord,'Value',1)
1419set(handles.ListCoord,'String',Tabchar)
1420MenuPlot_Callback(handles.geometry_calib, [], handles)
1421
1422% -----------------------------------------------------------------------
1423function MenuImportIntrinsic_Callback(hObject, eventdata, handles)
1424%------------------------------------------------------------------------
1425fileinput=browse_xml(hObject, eventdata, handles);
1426if isempty(fileinput)
1427    return
1428end
1429[s,errormsg]=imadoc2struct(fileinput,'GeometryCalib');
1430GeometryCalib=s.GeometryCalib;
1431display_intrinsic(GeometryCalib,handles)
1432
1433% -----------------------------------------------------------------------
1434function MenuImportAll_Callback(hObject, eventdata, handles)
1435%------------------------------------------------------------------------
1436fileinput=browse_xml(hObject, eventdata, handles)
1437if ~isempty(fileinput)
1438    loadfile(handles,fileinput)
1439end
1440
1441% -----------------------------------------------------------------------
1442% --- Executes on menubar option Import/Grid file: introduce previous grid files
1443function MenuGridFile_Callback(hObject, eventdata, handles)
1444% -----------------------------------------------------------------------
1445inputfile=browse_xml(hObject, eventdata, handles)
1446listfile=get(handles.coord_files,'string');
1447if isequal(listfile,{''})
1448    listfile={inputfile};
1449else
1450    listfile=[listfile;{inputfile}]%update the list of coord files
1451end
1452set(handles.coord_files,'string',listfile);
1453
1454%------------------------------------------------------------------------
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%------------------------------------------------------------------------
1472function fileinput=browse_xml(hObject, eventdata, handles)
1473%------------------------------------------------------------------------
1474fileinput=[];%default
1475oldfile=''; %default
1476UserData=get(handles.geometry_calib,'UserData')
1477if isfield(UserData,'XmlInputFile')
1478    oldfile=UserData.XmlInputFile;
1479end
1480[FileName, PathName, filterindex] = uigetfile( ...
1481       {'*.xml;*.mat', ' (*.xml,*.mat)';
1482       '*.xml',  '.xml files '; ...
1483        '*.mat',  '.mat matlab files '}, ...
1484        'Pick a file',oldfile);
1485fileinput=[PathName FileName];%complete file name
1486testblank=findstr(fileinput,' ');%look for blanks
1487if ~isempty(testblank)
1488    msgbox_uvmat('ERROR','forbidden input file name or path: no blank character allowed')
1489    return
1490end
1491sizf=size(fileinput);
1492if (~ischar(fileinput)||~isequal(sizf(1),1)),return;end
1493UserData.XmlInputFile=fileinput;
1494set(handles.geometry_calib,'UserData',UserData)%record current file foer further use of browser
1495
1496% -----------------------------------------------------------------------
1497function loadfile(handles,fileinput)
1498%------------------------------------------------------------------------
1499[s,errormsg]=imadoc2struct(fileinput,'GeometryCalib');
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
1536    calib=reshape(GeometryCalib.PointCoord',[],1);
1537    for ilist=1:numel(calib)
1538        CoordCell{ilist}=num2str(calib(ilist));
1539    end
1540    CoordCell=reshape(CoordCell,[],5);
1541    Tabchar=cell2tab(CoordCell,'    |    ');%transform cells into table ready for display
1542    MenuPlot_Callback(handles.geometry_calib, [], handles)
1543end
1544set(handles.calib_type,'Value',val_cal)
1545Tabchar=[Tabchar;{'......'}];
1546set(handles.ListCoord,'Value',1)
1547set(handles.ListCoord,'String',Tabchar)
1548
1549if isempty(CoordCell)% allow mouse action by default in the absence of input points
1550    set(handles.edit_append,'Value',1)
1551    set(handles.edit_append,'BackgroundColor',[1 1 0])
1552else % does not allow mouse action by default in the presence of input points
1553    set(handles.edit_append,'Value',0)
1554    set(handles.edit_append,'BackgroundColor',[0.7 0.7 0.7])
1555end
1556
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
Note: See TracBrowser for help on using the repository browser.