source: trunk/src/geometry_calib.m @ 191

Last change on this file since 191 was 191, checked in by sommeria, 13 years ago

introduce edit boxes to set axis limits. rationalisation of names FixScal?, FixLimits?....
introduce volume scan in calibration procedures

File size: 63.2 KB
Line 
1%'geometry_calib': associated to the GUI geometry_calib to perform geometric calibration from a set of reference points
2%------------------------------------------------------------------------
3% function hgeometry_calib = geometry_calib(inputfile,pos)
4%
5%OUTPUT:
6% hgeometry_calib=current handles of the GUI geometry_calib.fig
7%
8%INPUT:
9% inputfile: (optional) name of an xml file containing coordinates of reference points
10% pos: (optional) 4 element vector setting the 'Position' of the GUI
11%
12%A%AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
13%  Copyright Joel Sommeria, 2008, LEGI / CNRS-UJF-INPG, sommeria@coriolis-legi.org.
14%AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
15%     This file is part of the toolbox UVMAT.
16%
17%     UVMAT is free software; you can redistribute it and/or modify
18%     it under the terms of the GNU General Public License as published by
19%     the Free Software Foundation; either version 2 of the License, or
20%     (at your option) any later version.
21%
22%     UVMAT is distributed in the hope that it will be useful,
23%     but WITHOUT ANY WARRANTY; without even the implied warranty of
24%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
25%     GNU General Public License (file UVMAT/COPYING.txt) for more details.
26%AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
27
28function varargout = geometry_calib(varargin)
29% GEOMETRY_CALIB M-file for geometry_calib.fig
30%      GEOMETRY_CALIB, by itself, creates a MenuCoord GEOMETRY_CALIB or raises the existing
31%      singleton*.
32%
33%      H = GEOMETRY_CALIB returns the handle to a MenuCoord GEOMETRY_CALIB or the handle to
34%      the existing singleton*.
35%
36%      GEOMETRY_CALIB('CALLBACK',hObject,eventData,handles,...) calls the local
37%      function named CALLBACK in GEOMETRY_CALIB.M with the given input arguments.
38%
39%      GEOMETRY_CALIB('Property','Value',...) creates a MenuCoord GEOMETRY_CALIB or raises the
40%      existing singleton*.  Starting from the left, property value pairs are
41%      applied to the GUI before geometry_calib_OpeningFunction gets called.  An
42%      unrecognized property name or invalid value makes property application
43%      stop.  All inputs are passed to geometry_calib_OpeningFcn via varargin.
44%
45%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
46%      instance to run (singleton)".
47%
48% See also: GUIDE, GUIDATA, GUIHANDLES
49
50% Edit the above text to modify the response to help geometry_calib
51
52% Last Modified by GUIDE v2.5 05-Oct-2010 13:47:00
53
54% Begin initialization code - DO NOT edit
55gui_Singleton = 1;
56gui_State = struct('gui_Name',       mfilename, ...
57                   'gui_Singleton',  gui_Singleton, ...
58                   'gui_OpeningFcn', @geometry_calib_OpeningFcn, ...
59                   'gui_OutputFcn',  @geometry_calib_OutputFcn, ...
60                   'gui_LayoutFcn',  [] , ...
61                   'gui_Callback',   []);
62if nargin && ischar(varargin{1}) && ~isempty(regexp(varargin{1},'_Callback','once'))
63    gui_State.gui_Callback = str2func(varargin{1});
64end
65
66if nargout
67    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
68else
69    gui_mainfcn(gui_State, varargin{:});
70end
71% End initialization code - DO NOT edit
72
73
74% --- Executes just before geometry_calib is made visible.
75%INPUT:
76%handles: handles of the geometry_calib interface elements
77% PlotHandles: set of handles of the elements contolling the plotting
78% parameters on the uvmat interface (obtained by 'get_plot_handle.m')
79%------------------------------------------------------------------------
80function geometry_calib_OpeningFcn(hObject, eventdata, handles,inputfile,pos)
81%------------------------------------------------------------------------
82% Choose default command line output for geometry_calib
83
84handles.output = hObject;
85
86% Update handles structure
87guidata(hObject, handles);
88set(hObject,'DeleteFcn',{@closefcn})%
89
90%set the position of the interface
91if exist('pos','var')&& length(pos)>=4
92%     %pos_gui=get(hObject,'Position');
93%     pos_gui(1)=pos(1);
94%     pos_gui(2)=pos(2);
95    set(hObject,'Position',pos);
96end
97
98%set menu of calibration options
99set(handles.calib_type,'String',{'rescale';'linear';'3D_linear';'3D_quadr';'3D_extrinsic'})
100inputxml='';
101if exist('inputfile','var')&& ~isempty(inputfile)
102    struct.XmlInputFile=inputfile;
103    [Pathsub,RootFile,field_count,str2,str_a,str_b,ext]=name2display(inputfile);
104    if ~strcmp(ext,'.xml')
105        inputfile=[fullfile(Pathsub,RootFile) '.xml'];%xml file corresponding to the input file
106    end
107    set(handles.ListCoord,'String',{'......'})
108    if exist(inputfile,'file')
109        Heading=loadfile(handles,inputfile);% load the point coordiantes existing in the xml file
110        if isfield(Heading,'Campaign')&& ischar(Heading.Campaign)
111            struct.Campaign=Heading.Campaign;
112        end
113    end   
114    set(hObject,'UserData',struct)
115end
116
117set(handles.ListCoord,'KeyPressFcn',{@key_press_fcn,handles})%set keyboard action function
118
119
120%------------------------------------------------------------------------
121% --- Outputs from this function are returned to the command line.
122function varargout = geometry_calib_OutputFcn(hObject, eventdata, handles)
123%------------------------------------------------------------------------
124% Get default command line output from handles structure
125varargout{1} = handles.output;
126varargout{2}=handles;
127%
128%------------------------------------------------------------------------
129% executed when closing: set the parent interface button to value 0
130function closefcn(gcbo,eventdata)
131%------------------------------------------------------------------------
132huvmat=findobj(allchild(0),'Name','uvmat');
133if ~isempty(huvmat)
134    handles=guidata(huvmat);
135%     set(handles.MenuMask,'enable','on')
136%     set(handles.MenuGrid,'enable','on')
137%     set(handles.MenuObject,'enable','on')
138%     set(handles.MenuEdit,'enable','on')
139%     set(handles.edit,'enable','on')
140    hobject=findobj(handles.axes3,'tag','calib_points');
141    if ~isempty(hobject)
142        delete(hobject)
143    end
144    hobject=findobj(handles.axes3,'tag','calib_marker');
145    if ~isempty(hobject)
146        delete(hobject)
147    end   
148end
149
150%------------------------------------------------------------------------
151% --- Executes on button press in calibrate_lin.
152function APPLY_Callback(hObject, eventdata, handles)
153%------------------------------------------------------------------------
154%read the current calibration points
155Coord_cell=get(handles.ListCoord,'String');
156Object=read_geometry_calib(Coord_cell);
157Coord=Object.Coord;
158% apply the calibration, whose type is selected in  handles.calib_type
159if ~isempty(Coord)
160    calib_cell=get(handles.calib_type,'String');
161    val=get(handles.calib_type,'Value');
162    GeometryCalib=feval(['calib_' calib_cell{val}],Coord,handles);
163else
164    msgbox_uvmat('ERROR','No calibration points, abort')
165    return
166end
167Z_plane=[];
168if ~isempty(Coord)
169    %check error
170    X=Coord(:,1);
171    Y=Coord(:,2);
172    Z=Coord(:,3);
173    x_ima=Coord(:,4);
174    y_ima=Coord(:,5);
175    [Xpoints,Ypoints]=px_XYZ(GeometryCalib,X,Y,Z);
176    GeometryCalib.ErrorRms(1)=sqrt(mean((Xpoints-x_ima).*(Xpoints-x_ima)));
177    [GeometryCalib.ErrorMax(1),index(1)]=max(abs(Xpoints-x_ima));
178    GeometryCalib.ErrorRms(2)=sqrt(mean((Ypoints-y_ima).*(Ypoints-y_ima)));
179    [GeometryCalib.ErrorMax(2),index(2)]=max(abs(Ypoints-y_ima));
180    [EM,ind_dim]=max(GeometryCalib.ErrorMax);
181    index=index(ind_dim);
182    %set the Z position of the reference plane used for calibration
183    if isequal(max(Z),min(Z))%Z constant
184        Z_plane=Z(1);
185        GeometryCalib.NbSlice=1;
186        GeometryCalib.SliceCoord=[0 0 Z_plane];
187    end
188end
189%set the coordinate unit
190unitlist=get(handles.CoordUnit,'String');
191unit=unitlist{get(handles.CoordUnit,'value')};
192GeometryCalib.CoordUnit=unit;
193%record the points
194GeometryCalib.SourceCalib.PointCoord=Coord;
195display_intrinsic(GeometryCalib,handles)%display calibration intrinsic parameters
196
197% Display extrinsinc parameters (rotation and translation of camera with  respect to the phys coordiantes)
198set(handles.Tx,'String',num2str(GeometryCalib.Tx_Ty_Tz(1),4))
199set(handles.Ty,'String',num2str(GeometryCalib.Tx_Ty_Tz(2),4))
200set(handles.Tz,'String',num2str(GeometryCalib.Tx_Ty_Tz(3),4))
201set(handles.Phi,'String',num2str(GeometryCalib.omc(1),4))
202set(handles.Theta,'String',num2str(GeometryCalib.omc(2),4))
203set(handles.Psi,'String',num2str(GeometryCalib.omc(3),4))
204
205% store the calibration data, by default in the xml file of the currently displayed image
206huvmat=findobj(allchild(0),'Name','uvmat');
207UvData=get(huvmat,'UserData');
208NbSlice_j=1;%default
209ZStart=Z_plane;
210ZEnd=Z_plane;
211if isfield(UvData,'XmlData')
212    UvData.XmlData
213    if isfield(UvData.XmlData,'TranslationMotor')
214        NbSlice_j=UvData.XmlData.TranslationMotor.Nbslice;
215        ZStart=UvData.XmlData.TranslationMotor.ZStart;
216        ZEnd=UvData.XmlData.TranslationMotor.ZEnd;
217    end
218end
219hhuvmat=guidata(huvmat);%handles of elements in the GUI uvmat
220RootPath='';
221RootFile='';
222if ~isempty(hhuvmat.RootPath)&& ~isempty(hhuvmat.RootFile)
223    testhandle=1;
224    RootPath=get(hhuvmat.RootPath,'String');
225    RootFile=get(hhuvmat.RootFile,'String');
226    filebase=fullfile(RootPath,RootFile);
227    outputfile=[filebase '.xml'];%xml file associated with the currently displayed image
228else
229    question={'save the calibration data and point coordinates in'};
230    def={fullfile(RootPath,'ObjectCalib.xml')};
231    options.Resize='on';
232    answer=inputdlg(question,'save average in a new file',1,def,options);
233    outputfile=answer{1};
234end
235answer=msgbox_uvmat('INPUT_Y-N',{[outputfile ' updated with calibration data'];...
236    ['Error rms (along x,y)=' num2str(GeometryCalib.ErrorRms) ' pixels'];...
237    ['Error max (along x,y)=' num2str(GeometryCalib.ErrorMax) ' pixels']});
238
239%% record the calibration parameters and display the current image of uvmat in the new phys coordinates
240if strcmp(answer,'Yes')
241    if strcmp(calib_cell{val}(1:2),'3D')%set the plane position for 3D (projection) calibration
242       input_key={'Z (first position)','Z (last position)','Z (water surface)', 'refractive index','NbSlice','volume scan (y/n)','tilt angle'};
243       input_val=[{num2str(ZEnd/10)} {num2str(ZStart/10)} {num2str(ZStart/10)} {'1.33'} num2str(NbSlice_j) {'y'} {'0'}];
244        answer=inputdlg(input_key,'slice position(s)',ones(1,7), input_val,'on');
245        %answer_1=msgbox_uvmat('INPUT_TXT',' Z= ',num2str(Z_plane));
246        GeometryCalib.NbSlice=str2double(answer{5});
247        GeometryCalib.VolumeScan=answer{6};
248        if isempty(answer)
249            Z_plane=0; %default
250        else
251            Z_plane=linspace(str2double(answer{1}),str2double(answer{2}),GeometryCalib.NbSlice);
252        end     
253        GeometryCalib.SliceCoord=Z_plane'*[0 0 1]
254    end
255    errormsg=update_imadoc(GeometryCalib,outputfile);% introduce the calibration data in the xml file
256    if ~strcmp(errormsg,'')
257        msgbox_uvmat('ERROR',errormsg);
258    end
259   
260    %display image with new calibration in the currently opened uvmat interface
261    hhh=findobj(hhuvmat.axes3,'Tag','calib_marker');% delete calib points and markers
262    if ~isempty(hhh)
263        delete(hhh);     
264    end
265    hhh=findobj(hhuvmat.axes3,'Tag','calib_points');
266    if ~isempty(hhh)
267        delete(hhh);
268    end
269    set(hhuvmat.FixLimits,'Value',0)% put FixedLimits option to 'off'
270    set(hhuvmat.FixLimits,'BackgroundColor',[0.7 0.7 0.7])
271    UserData=get(handles.geometry_calib,'UserData');
272    UserData.XmlInputFile=outputfile;%save the current xml file name
273    set(handles.geometry_calib,'UserData',UserData)
274    uvmat('RootPath_Callback',hObject,eventdata,hhuvmat); %file input with xml reading  in uvmat, show the image in phys coordinates
275    MenuPlot_Callback(hObject, eventdata, handles)
276    set(handles.ListCoord,'Value',index)% indicate in the list the point with max deviation (possible mistake)
277    ListCoord_Callback(hObject, eventdata, handles)
278    figure(handles.geometry_calib)
279end
280
281%------------------------------------------------------------------
282% --- Executes on button press in calibrate_lin.
283
284function REPLICATE_Callback(hObject, eventdata, handles)
285%------------------------------------------------------------------------
286
287%% Apply calibration
288calib_cell=get(handles.calib_type,'String'); %#ok<NASGU>
289val=get(handles.calib_type,'Value'); %#ok<NASGU>
290
291%read the current calibration points
292Coord_cell=get(handles.ListCoord,'String');
293Object=read_geometry_calib(Coord_cell);
294Coord=Object.Coord;
295
296% apply the calibration, whose type is selected in  handles.calib_type
297if ~isempty(Coord)
298    calib_cell=get(handles.calib_type,'String');
299    val=get(handles.calib_type,'Value');
300    GeometryCalib=feval(['calib_' calib_cell{val}],Coord,handles);
301else
302    msgbox_uvmat('ERROR','No calibration points, abort')
303    return
304end
305
306if ~isempty(Coord)
307    %check error
308    X=Coord(:,1);
309    Y=Coord(:,2);
310    Z=Coord(:,3);
311    x_ima=Coord(:,4);
312    y_ima=Coord(:,5);
313    [Xpoints,Ypoints]=px_XYZ(GeometryCalib,X,Y,Z);
314    GeometryCalib.ErrorRms(1)=sqrt(mean((Xpoints-x_ima).*(Xpoints-x_ima)));
315    [GeometryCalib.ErrorMax(1),index(1)]=max(abs(Xpoints-x_ima));
316    GeometryCalib.ErrorRms(2)=sqrt(mean((Ypoints-y_ima).*(Ypoints-y_ima)));
317    [GeometryCalib.ErrorMax(2),index(2)]=max(abs(Ypoints-y_ima));
318    [EM,ind_dim]=max(GeometryCalib.ErrorMax);
319%     index=index(ind_dim);
320    %set the Z position of the reference plane used for calibration
321    Z_plane=[];
322    if isequal(max(Z),min(Z))
323        Z_plane=Z(1);
324    end
325    answer_1=msgbox_uvmat('INPUT_TXT',' Z= ',num2str(Z_plane));
326    Z_plane=str2double(answer_1);
327    GeometryCalib.NbSlice=1;
328    GeometryCalib.SliceCoord=[0 0 Z_plane];
329    %set the coordinate unit
330    unitlist=get(handles.CoordUnit,'String');
331    unit=unitlist{get(handles.CoordUnit,'value')};
332    GeometryCalib.CoordUnit=unit;
333    %record the points
334    GeometryCalib.SourceCalib.PointCoord=Coord;
335end
336
337%% display calibration paprameters
338display_intrinsic(GeometryCalib,handles)%display calibration intrinsic parameters
339
340% Display extrinsinc parameters (rotation and translation of camera with  respect to the phys coordiantes)
341set(handles.Tx,'String',num2str(GeometryCalib.Tx_Ty_Tz(1),4))
342set(handles.Ty,'String',num2str(GeometryCalib.Tx_Ty_Tz(2),4))
343set(handles.Tz,'String',num2str(GeometryCalib.Tx_Ty_Tz(3),4))
344set(handles.Phi,'String',num2str(GeometryCalib.omc(1),4))
345set(handles.Theta,'String',num2str(GeometryCalib.omc(2),4))
346set(handles.Psi,'String',num2str(GeometryCalib.omc(3),4))
347
348%% open the GUI dataview
349h_dataview=findobj(allchild(0),'name','dataview');
350if ~isempty(h_dataview)
351    delete(h_dataview)
352end
353CalibData=get(handles.geometry_calib,'UserData');%read the calibration image source on the interface userdata
354InputFile='';
355if isfield(CalibData,'XmlInputFile')
356    InputDir=fileparts(CalibData.XmlInputFile);
357    [InputDir,DirName]=fileparts(InputDir);
358end
359SubCampaignTest='n'; %default
360testup=0;
361if isfield(CalibData,'SubCampaign')
362    SubCampaignTest='y';
363    dir_ref=CalibData.SubCampaign;
364    testup=1;
365elseif isfield(CalibData,'Campaign')
366    dir_ref=CalibData.Campaign;
367    testup=1;
368end
369while testup
370    [InputDir,DirName]=fileparts(InputDir);
371    if strcmp(DirName,dir_ref)
372        break
373    end
374end
375InputDir=fullfile(InputDir,DirName);
376answer=msgbox_uvmat('INPUT_TXT','Campaign ?',InputDir);
377if strcmp(answer,'Cancel')
378    return
379end
380
381dataview(answer,SubCampaignTest,GeometryCalib);
382       
383%     if isfield(Heading,'Device') && isequal([filename ext],Heading.Device)
384%         [XmlInput,filename,ext]=fileparts(XmlInput);
385%         Device=Heading.Device;
386%     end
387%     if isfield(Heading,'Experiment') && isequal([filename ext],Heading.Experiment)
388%         [PP,filename,ext]=fileparts(XmlInput);
389%     end
390%     testinput=0;
391%     if isfield(Heading,'SubCampaign') && isequal([filename ext],Heading.SubCampaign)
392%         SubCampaignTest='y';
393%         testinput=1;
394%     elseif isfield(Heading,'Campaign') && isequal([filename ext],Heading.Campaign)
395%         testinput=1;
396% %     end
397% end
398% if ~testinput
399%     filename='PROJETS';%default
400%     if isfield(CalibData,'XmlInputFile')
401%          [pp,filename]=fileparts(CalibData.XmlInputFile);
402%     end
403%     while ~isequal(filename,'PROJETS') && numel(filename)>1
404%         filename_1=filename;
405%         pp_1=pp;
406%         [pp,filename]=fileparts(pp);
407%     end
408%     XmlInput=fullfile(pp_1,filename_1);
409%     testinput=1;
410% end
411% if testinput
412%     outcome=dataview(XmlInput,SubCampaignTest,GeometryCalib);
413% end
414
415%------------------------------------------------------------------------
416% determine the parameters for a calibration by an affine function (rescaling and offset, no rotation)
417function GeometryCalib=calib_rescale(Coord,handles)
418%------------------------------------------------------------------------
419X=Coord(:,1);
420Y=Coord(:,2);% Z not used
421x_ima=Coord(:,4);
422y_ima=Coord(:,5);
423[px,sx]=polyfit(X,x_ima,1);
424[py,sy]=polyfit(Y,y_ima,1);
425T_x=px(2);
426T_y=py(2);
427GeometryCalib.CalibrationType='rescale';
428GeometryCalib.fx_fy=[px(1) py(1)];%.fx_fy corresponds to pxcm along x and y
429GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function
430GeometryCalib.Tx_Ty_Tz=[px(2)/px(1) py(2)/py(1) 1];
431%GeometryCalib.R=[1,0,0;0,1,0;0,0,0];
432GeometryCalib.omc=[0 0 0];
433
434%------------------------------------------------------------------------
435% determine the parameters for a calibration by a linear transform matrix (rescale and rotation)
436function GeometryCalib=calib_linear(Coord,handles) %TO UPDATE
437%------------------------------------------------------------------------
438X=Coord(:,1);
439Y=Coord(:,2);% Z not used
440x_ima=Coord(:,4);
441y_ima=Coord(:,5);
442XY_mat=[ones(size(X)) X Y];
443a_X1=XY_mat\x_ima; %transformation matrix for X
444% x1=XY_mat*a_X1;%reconstruction
445% err_X1=max(abs(x1-x_ima));%error
446a_Y1=XY_mat\y_ima;%transformation matrix for X
447% y1=XY_mat*a_Y1;
448% err_Y1=max(abs(y1-y_ima));%error
449% R=[a_X1(2),a_X1(3),0;a_Y1(2),a_Y1(3),0;0,0,1];
450R=[a_X1(2),a_X1(3);a_Y1(2),a_Y1(3)];
451norm=abs(det(R));
452GeometryCalib.CalibrationType='linear';
453GeometryCalib.fx_fy(1)=sqrt((a_X1(2)/a_Y1(3))*norm);
454GeometryCalib.fx_fy(2)=(a_Y1(3)/a_X1(2))*GeometryCalib.fx_fy(1);
455GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function
456GeometryCalib.Tx_Ty_Tz=[a_X1(1) a_Y1(1) 1];
457R(1,:)=R(1,:)/GeometryCalib.fx_fy(1);
458R(2,:)=R(2,:)/GeometryCalib.fx_fy(2);
459R=[R;[0 0]];
460GeometryCalib.R=[R [0;0;1]];
461GeometryCalib.omc=(180/pi)*[acos(GeometryCalib.R(1,1)) 0 0];
462%------------------------------------------------------------------------
463% determine the tsai parameters for a view normal to the grid plane
464% NOT USED
465function GeometryCalib=calib_normal(Coord,handles)
466%------------------------------------------------------------------------
467Calib.f1=str2num(get(handles.fx,'String'));
468Calib.f2=str2num(get(handles.fy,'String'));
469Calib.k=str2num(get(handles.kc,'String'));
470Calib.Cx=str2num(get(handles.Cx,'String'));
471Calib.Cy=str2num(get(handles.Cy,'String'));
472%default
473if isempty(Calib.f1)
474    Calib.f1=25/0.012;
475end
476if isempty(Calib.f2)
477    Calib.f2=25/0.012;
478end
479if isempty(Calib.k)
480    Calib.k=0;
481end
482if isempty(Calib.Cx)||isempty(Calib.Cy)
483    huvmat=findobj(allchild(0),'Tag','uvmat');
484    hhuvmat=guidata(huvmat);
485    Calib.Cx=str2num(get(hhuvmat.npx,'String'))/2;
486    Calib.Cx=str2num(get(hhuvmat.npy,'String'))/2;
487end   
488%tsai parameters
489Calib.dpx=0.012;%arbitrary
490Calib.dpy=0.012;
491Calib.sx=Calib.f1*Calib.dpx/(Calib.f2*Calib.dpy);
492Calib.f=Calib.f2*Calib.dpy;
493Calib.kappa1=Calib.k/(Calib.f*Calib.f);
494
495%initial guess
496X=Coord(:,1);
497Y=Coord(:,2);
498Zmean=mean(Coord(:,3));
499x_ima=Coord(:,4)-Calib.Cx;
500y_ima=Coord(:,5)-Calib.Cy;
501XY_mat=[ones(size(X)) X Y];
502a_X1=XY_mat\x_ima; %transformation matrix for X
503a_Y1=XY_mat\y_ima;%transformation matrix for Y
504R=[a_X1(2),a_X1(3),0;a_Y1(2),a_Y1(3),0;0,0,-1];% rotation+ z axis reversal (upward)
505norm=sqrt(det(-R));
506calib_param(1)=0;% quadratic distortion
507calib_param(2)=a_X1(1);
508calib_param(3)=a_Y1(1);
509calib_param(4)=Calib.f/(norm*Calib.dpx)-R(3,3)*Zmean;
510calib_param(5)=angle(a_X1(2)+1i*a_X1(3));
511display(['initial guess=' num2str(calib_param)])
512
513%optimise the parameters: minimisation of error
514calib_param = fminsearch(@(calib_param) error_calib(calib_param,Calib,Coord),calib_param);
515
516GeometryCalib.CalibrationType='tsai_normal';
517GeometryCalib.focal=Calib.f;
518GeometryCalib.dpx_dpy=[Calib.dpx Calib.dpy];
519GeometryCalib.Cx_Cy=[Calib.Cx Calib.Cy];
520GeometryCalib.sx=Calib.sx;
521GeometryCalib.kappa1=calib_param(1);
522GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function
523GeometryCalib.Tx_Ty_Tz=[calib_param(2) calib_param(3) calib_param(4)];
524alpha=calib_param(5);
525GeometryCalib.R=[cos(alpha) sin(alpha) 0;-sin(alpha) cos(alpha) 0;0 0 -1];
526
527%------------------------------------------------------------------------
528function GeometryCalib=calib_3D_linear(Coord,handles)
529%------------------------------------------------------------------
530path_uvmat=which('uvmat');% check the path detected for source file uvmat
531path_UVMAT=fileparts(path_uvmat); %path to UVMAT
532huvmat=findobj(allchild(0),'Tag','uvmat');
533hhuvmat=guidata(huvmat);
534coord_files=get(handles.coord_files,'String');
535if ischar(coord_files)
536    coord_files={coord_files};
537end
538if isempty(coord_files{1}) || isequal(coord_files,{''})
539    coord_files={};
540end
541%retrieve the calibration points stored in the files listed in the popup list coord_files
542x_1=Coord(:,4:5)';%px coordinates of the ref points
543nx=str2num(get(hhuvmat.npx,'String'));
544ny=str2num(get(hhuvmat.npy,'String'));
545x_1(2,:)=ny-x_1(2,:);%reverse the y image coordinates
546X_1=Coord(:,1:3)';%phys coordinates of the ref points
547n_ima=numel(coord_files)+1;
548if ~isempty(coord_files)
549    msgbox_uvmat('CONFIRMATION',['The xy coordinates of the calibration points in ' num2str(n_ima) ' planes will be used'])
550    for ifile=1:numel(coord_files)
551    t=xmltree(coord_files{ifile});
552    s=convert(t);%convert to matlab structure
553        if isfield(s,'GeometryCalib')
554            if isfield(s.GeometryCalib,'SourceCalib')
555                if isfield(s.GeometryCalib.SourceCalib,'PointCoord')
556                PointCoord=s.GeometryCalib.SourceCalib.PointCoord;
557                Coord_file=zeros(length(PointCoord),5);%default
558                for i=1:length(PointCoord)
559                    line=str2num(PointCoord{i});
560                    Coord_file(i,4:5)=line(4:5);%px x
561                    Coord_file(i,1:3)=line(1:3);%phys x
562                end
563                eval(['x_' num2str(ifile+1) '=Coord_file(:,4:5)'';']);
564                eval(['x_' num2str(ifile+1) '(2,:)=ny-x_' num2str(ifile+1) '(2,:);' ]);
565                eval(['X_' num2str(ifile+1) '=Coord_file(:,1:3)'';']);
566                end
567            end
568        end
569    end
570end
571n_ima=numel(coord_files)+1;
572est_dist=[0;0;0;0;0];
573est_aspect_ratio=0;
574est_fc=[1;1];
575%fc=[25;25]/0.012;
576center_optim=0;
577run(fullfile(path_UVMAT,'toolbox_calib','go_calib_optim'));
578GeometryCalib.CalibrationType='3D_linear';
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%------------------------------------------------------------------------
597function GeometryCalib=calib_3D_quadr(Coord,handles)
598%------------------------------------------------------------------
599
600path_uvmat=which('uvmat');% check the path detected for source file uvmat
601path_UVMAT=fileparts(path_uvmat); %path to UVMAT
602huvmat=findobj(allchild(0),'Tag','uvmat');
603hhuvmat=guidata(huvmat);
604% check_cond=0;
605coord_files=get(handles.coord_files,'String');
606if ischar(coord_files)
607    coord_files={coord_files};
608end
609if isempty(coord_files{1}) || isequal(coord_files,{''})
610    coord_files={};
611end
612
613%retrieve the calibration points stored in the files listed in the popup list coord_files
614x_1=Coord(:,4:5)';%px coordinates of the ref points
615nx=str2num(get(hhuvmat.npx,'String'));
616ny=str2num(get(hhuvmat.npy,'String'));
617x_1(2,:)=ny-x_1(2,:);%reverse the y image coordinates
618X_1=Coord(:,1:3)';%phys coordinates of the ref points
619n_ima=numel(coord_files)+1;
620if ~isempty(coord_files)
621    msgbox_uvmat('CONFIRMATION',['The xy coordinates of the calibration points in ' num2str(n_ima) ' planes will be used'])
622    for ifile=1:numel(coord_files)
623    t=xmltree(coord_files{ifile});
624    s=convert(t);%convert to matlab structure
625        if isfield(s,'GeometryCalib')
626            if isfield(s.GeometryCalib,'SourceCalib')
627                if isfield(s.GeometryCalib.SourceCalib,'PointCoord')
628                PointCoord=s.GeometryCalib.SourceCalib.PointCoord;
629                Coord_file=zeros(length(PointCoord),5);%default
630                for i=1:length(PointCoord)
631                    line=str2num(PointCoord{i});
632                    Coord_file(i,4:5)=line(4:5);%px x
633                    Coord_file(i,1:3)=line(1:3);%phys x
634                end
635                eval(['x_' num2str(ifile+1) '=Coord_file(:,4:5)'';']);
636                eval(['x_' num2str(ifile+1) '(2,:)=ny-x_' num2str(ifile+1) '(2,:);' ]);
637                eval(['X_' num2str(ifile+1) '=Coord_file(:,1:3)'';']);
638                end
639            end
640        end
641    end
642end
643n_ima=numel(coord_files)+1;
644est_dist=[1;0;0;0;0];
645est_aspect_ratio=1;
646%est_fc=[0;0];
647%fc=[25;25]/0.012;
648center_optim=0;
649run(fullfile(path_UVMAT,'toolbox_calib','go_calib_optim'));
650
651GeometryCalib.CalibrationType='3D_quadr';
652GeometryCalib.fx_fy=fc';
653%GeometryCalib.focal=fc(2);
654%GeometryCalib.dpx_dpy=[1 1];
655GeometryCalib.Cx_Cy=cc';
656%GeometryCalib.sx=fc(1)/fc(2);
657GeometryCalib.kc=kc(1);
658%GeometryCalib.kappa1=-kc(1)/fc(2)^2;
659GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function
660GeometryCalib.Tx_Ty_Tz=Tc_1';
661GeometryCalib.R=Rc_1;
662GeometryCalib.R(2,1:3)=-GeometryCalib.R(2,1:3);%inversion of the y image coordinate
663GeometryCalib.Tx_Ty_Tz(2)=-GeometryCalib.Tx_Ty_Tz(2);%inversion of the y image coordinate
664GeometryCalib.Cx_Cy(2)=ny-GeometryCalib.Cx_Cy(2);%inversion of the y image coordinate
665GeometryCalib.omc=(180/pi)*omc_1;%angles in degrees
666GeometryCalib.ErrorRMS=[];
667GeometryCalib.ErrorMax=[];
668
669
670%------------------------------------------------------------------------
671function GeometryCalib=calib_3D_extrinsic(Coord,handles)
672%------------------------------------------------------------------
673path_uvmat=which('geometry_calib');% check the path detected for source file uvmat
674path_UVMAT=fileparts(path_uvmat); %path to UVMAT
675x_1=double(Coord(:,4:5)');%image coordiantes
676X_1=double(Coord(:,1:3)');% phys coordinates
677huvmat=findobj(allchild(0),'Tag','uvmat');
678hhuvmat=guidata(huvmat);
679ny=str2double(get(hhuvmat.npy,'String'));
680x_1(2,:)=ny-x_1(2,:);%reverse the y image coordinates
681n_ima=1;
682GeometryCalib.CalibrationType='3D_extrinsic';
683GeometryCalib.fx_fy(1)=str2num(get(handles.fx,'String'));
684GeometryCalib.fx_fy(2)=str2num(get(handles.fy,'String'));
685GeometryCalib.Cx_Cy(1)=str2num(get(handles.Cx,'String'));
686GeometryCalib.Cx_Cy(2)=str2num(get(handles.Cy,'String'));
687GeometryCalib.kc=str2num(get(handles.kc,'String'));
688fct_path=fullfile(path_UVMAT,'toolbox_calib');
689addpath(fct_path)
690GeometryCalib.Cx_Cy(2)=ny-GeometryCalib.Cx_Cy(2);%reverse Cx_Cy(2) for calibration (inversion of px ordinate)
691% [omc1,Tc1,Rc1,H,x,ex,JJ] = compute_extrinsic(x_1,X_1,...
692%     [Calib.f Calib.f*Calib.sx]',...
693%     [Calib.Cx Calib.Cy]',...
694%     [-Calib.kappa1*Calib.f^2 0 0 0 0]);
695[omc,Tc1,Rc1,H,x,ex,JJ] = compute_extrinsic(x_1,X_1,...
696    (GeometryCalib.fx_fy)',GeometryCalib.Cx_Cy',[GeometryCalib.kc 0 0 0 0]);
697rmpath(fct_path);
698GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function
699GeometryCalib.Tx_Ty_Tz=Tc1';
700%inversion of z axis
701GeometryCalib.R=Rc1;
702GeometryCalib.R(2,1:3)=-GeometryCalib.R(2,1:3);%inversion of the y image coordinate
703GeometryCalib.Tx_Ty_Tz(2)=-GeometryCalib.Tx_Ty_Tz(2);%inversion of the y image coordinate
704GeometryCalib.Cx_Cy(2)=ny-GeometryCalib.Cx_Cy(2);%inversion of the y image coordinate
705GeometryCalib.omc=(180/pi)*omc';
706%GeometryCalib.R(3,1:3)=-GeometryCalib.R(3,1:3);%inversion for z upward
707
708
709
710%------------------------------------------------------------------------
711%function GeometryCalib=calib_tsai_heikkila(Coord)
712% TEST: NOT IMPLEMENTED
713%------------------------------------------------------------------
714% path_uvmat=which('uvmat');% check the path detected for source file uvmat
715% path_UVMAT=fileparts(path_uvmat); %path to UVMAT
716% path_calib=fullfile(path_UVMAT,'toolbox_calib_heikkila');
717% addpath(path_calib)
718% npoints=size(Coord,1);
719% Coord(:,1:3)=10*Coord(:,1:3);
720% Coord=[Coord zeros(npoints,2) -ones(npoints,1)];
721% [par,pos,iter,res,er,C]=cacal('dalsa',Coord);
722% GeometryCalib.CalibrationType='tsai';
723% GeometryCalib.focal=par(2);
724
725
726%--------------------------------------------------------------------------
727function GeometryCalib=calib_tsai(Coord,handles)% OBSOLETE: old version using gauthier's bianry ccal_fo
728% NOT USED
729%------------------------------------------------------------------------
730%TSAI
731path_uvmat=which('uvmat');% check the path detected for source file uvmat
732path_UVMAT=fileparts(path_uvmat); %path to UVMAT
733xmlfile=fullfile(path_UVMAT,'PARAM.xml');%name of the file containing names of binary executables
734if exist(xmlfile,'file')
735    t=xmltree(xmlfile);% read the (xml) file containing names of binary executables
736    sparam=convert(t);% convert to matlab structure
737end
738if ~isfield(sparam,'GeometryCalibBin')
739    msgbox_uvmat('ERROR',['calibration program <GeometryCalibBin> undefined in parameter file ' xmlfile])
740    return
741end
742Tsai_exe=sparam.GeometryCalibBin;
743if ~exist(Tsai_exe,'file')%the binary is defined in /bin, default setting
744     Tsai_exe=fullfile(path_UVMAT,Tsai_exe);
745end
746if ~exist(Tsai_exe,'file')
747    msgbox_uvmat('ERROR',['calibration program ' sparam.GeometryCalibBin ' defined in PARAM.xml does not exist'])
748    return
749end
750
751textcoord=num2str(Coord,4);
752dlmwrite('t.txt',textcoord,''); 
753% ['!' Tsai_exe ' -fx 0 -fy t.txt']
754eval(['!' Tsai_exe ' -f t.txt > tsaicalib.log']);
755if ~exist('calib.dat','file')
756    msgbox_uvmat('ERROR','no output from calibration program Tsai_exe: possibly too few points')
757end
758calibdat=dlmread('calib.dat');
759delete('calib.dat')
760%delete('t.txt')
761GeometryCalib.CalibrationType='tsai';
762GeometryCalib.focal=calibdat(10);
763GeometryCalib.dpx_dpy=[calibdat(5) calibdat(6)];
764GeometryCalib.Cx_Cy=[calibdat(7) calibdat(8)];
765GeometryCalib.sx=calibdat(9);
766GeometryCalib.kappa1=calibdat(11);
767GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function
768GeometryCalib.Tx_Ty_Tz=[calibdat(12) calibdat(13) calibdat(14)];
769Rx_Ry_Rz=calibdat(15:17);
770sa = sin(Rx_Ry_Rz(1)) ;
771ca=cos(Rx_Ry_Rz(1));
772sb=sin(Rx_Ry_Rz(2));
773cb =cos(Rx_Ry_Rz(2));
774sg =sin(Rx_Ry_Rz(3));
775cg =cos(Rx_Ry_Rz(3));
776r1 = cb * cg;
777r2 = cg * sa * sb - ca * sg;
778r3 = sa * sg + ca * cg * sb;
779r4 = cb * sg;
780r5 = sa * sb * sg + ca * cg;
781r6 = ca * sb * sg - cg * sa;
782r7 = -sb;
783r8 = cb * sa;
784r9 = ca * cb;
785%EN DEDUIRE MATRICE R ??
786GeometryCalib.R=[r1,r2,r3;r4,r5,r6;r7,r8,r9];
787
788%------------------------------------------------------------------------
789% --- determine the rms of calibration error
790function ErrorRms=error_calib(calib_param,Calib,Coord)
791%calib_param: vector of free calibration parameters (to optimise)
792%Calib: structure of the given calibration parameters
793%Coord: list of phys coordinates (columns 1-3, and pixel coordinates (columns 4-5)
794Calib.f=25;
795Calib.dpx=0.012;
796Calib.dpy=0.012;
797Calib.sx=1;
798Calib.Cx=512;
799Calib.Cy=512;
800Calib.kappa1=calib_param(1);
801Calib.Tx=calib_param(2);
802Calib.Ty=calib_param(3);
803Calib.Tz=calib_param(4);
804alpha=calib_param(5);
805Calib.R=[cos(alpha) sin(alpha) 0;-sin(alpha) cos(alpha) 0;0 0 -1];
806
807X=Coord(:,1);
808Y=Coord(:,2);
809Z=Coord(:,3);
810x_ima=Coord(:,4);
811y_ima=Coord(:,5);
812[Xpoints,Ypoints]=px_XYZ(Calib,X,Y,Z);
813ErrorRms(1)=sqrt(mean((Xpoints-x_ima).*(Xpoints-x_ima)));
814ErrorRms(2)=sqrt(mean((Ypoints-y_ima).*(Ypoints-y_ima)));
815ErrorRms=mean(ErrorRms);
816
817%------------------------------------------------------------------------
818function XImage_Callback(hObject, eventdata, handles)
819%------------------------------------------------------------------------
820update_list(hObject, eventdata,handles)
821
822%------------------------------------------------------------------------
823function YImage_Callback(hObject, eventdata, handles)
824%------------------------------------------------------------------------
825update_list(hObject, eventdata,handles)
826
827%------------------------------------------------------------------------
828% --- Executes on button press in STORE.
829function STORE_Callback(hObject, eventdata, handles)
830Coord_cell=get(handles.ListCoord,'String');
831Object=read_geometry_calib(Coord_cell);
832unitlist=get(handles.CoordUnit,'String');
833unit=unitlist{get(handles.CoordUnit,'value')};
834GeometryCalib.CoordUnit=unit;
835GeometryCalib.SourceCalib.PointCoord=Object.Coord;
836huvmat=findobj(allchild(0),'Name','uvmat');
837hhuvmat=guidata(huvmat);%handles of elements in the GUI uvmat
838% RootPath='';
839% RootFile='';
840if ~isempty(hhuvmat.RootPath)&& ~isempty(hhuvmat.RootFile)
841    testhandle=1;
842    RootPath=get(hhuvmat.RootPath,'String');
843    RootFile=get(hhuvmat.RootFile,'String');
844    filebase=fullfile(RootPath,RootFile);
845    while exist([filebase '.xml'],'file')
846        filebase=[filebase '~'];
847    end
848    outputfile=[filebase '.xml'];
849    errormsg=update_imadoc(GeometryCalib,outputfile);
850    if ~strcmp(errormsg,'')
851        msgbox_uvmat('ERROR',errormsg);
852    end
853    listfile=get(handles.coord_files,'string');
854    if isequal(listfile,{''})
855        listfile={outputfile};
856    else
857        listfile=[listfile;{outputfile}];%update the list of coord files
858    end
859    set(handles.coord_files,'string',listfile);
860end
861set(handles.ListCoord,'Value',1)% refresh the display of coordinates
862set(handles.ListCoord,'String',{'......'})
863
864% --------------------------------------------------------------------
865% --- Executes on button press in CLEAR_PTS: clear the list of calibration points
866function CLEAR_PTS_Callback(hObject, eventdata, handles)
867% --------------------------------------------------------------------
868set(handles.ListCoord,'Value',1)% refresh the display of coordinates
869set(handles.ListCoord,'String',{'......'})
870MenuPlot_Callback(hObject, eventdata, handles)
871
872%------------------------------------------------------------------------
873% --- Executes on button press in CLEAR.
874function CLEAR_Callback(hObject, eventdata, handles)
875%------------------------------------------------------------------------
876set(handles.coord_files,'Value',1)
877set(handles.coord_files,'String',{''})
878
879%------------------------------------------------------------------------
880function XObject_Callback(hObject, eventdata, handles)
881%------------------------------------------------------------------------
882update_list(hObject, eventdata,handles)
883
884%------------------------------------------------------------------------
885function YObject_Callback(hObject, eventdata, handles)
886%------------------------------------------------------------------------
887update_list(hObject, eventdata,handles)
888
889%------------------------------------------------------------------------
890function ZObject_Callback(hObject, eventdata, handles)
891%------------------------------------------------------------------------
892update_list(hObject, eventdata,handles)
893
894%------------------------------------------------------------------------
895function update_list(hObject, eventdata, handles)
896%------------------------------------------------------------------------
897newval(4)=str2double(get(handles.XImage,'String'));
898newval(5)=str2double(get(handles.YImage,'String'));
899newval(1)=str2double(get(handles.XObject,'String'));
900newval(2)=str2double(get(handles.YObject,'String'));
901newval(3)=str2double(get(handles.ZObject,'String'));
902if isnan(newval(3))
903    newval(3)=0;%put z to 0 by default
904end
905Coord=get(handles.ListCoord,'String');
906Coord(end)=[]; %remove last string '.....'
907val=get(handles.ListCoord,'Value');
908data=read_geometry_calib(Coord);
909data.Coord(val,:)=newval;
910for i=1:size(data.Coord,1)
911    for j=1:5
912          Coord_cell{i,j}=num2str(data.Coord(i,j),4);%display coordiantes with 4 digits
913    end
914end
915
916Tabchar=cell2tab(Coord_cell,' | ');
917Tabchar=[Tabchar ;{'......'}];
918set(handles.ListCoord,'String',Tabchar)
919
920%update the plot
921ListCoord_Callback(hObject, eventdata, handles)
922MenuPlot_Callback(hObject, eventdata, handles)
923
924%------------------------------------------------------------------------
925% --- Executes on selection change in ListCoord.
926function ListCoord_Callback(hObject, eventdata, handles)
927%------------------------------------------------------------------------
928huvmat=findobj(allchild(0),'Name','uvmat');%find the current uvmat interface handle
929hplot=findobj(huvmat,'Tag','axes3');%main plotting axis of uvmat
930hhh=findobj(hplot,'Tag','calib_marker');
931Coord_cell=get(handles.ListCoord,'String');
932val=get(handles.ListCoord,'Value');
933if numel(val)>1
934    return %no action if several lines have been selected
935end
936coord_str=Coord_cell{val};
937k=findstr(' | ',coord_str);
938if isempty(k)%last line '.....' selected
939    if ~isempty(hhh)
940        delete(hhh)%delete the circle marker
941    end
942    return
943end
944%fill the edit boxex
945set(handles.XObject,'String',coord_str(1:k(1)-1))
946set(handles.YObject,'String',coord_str(k(1)+3:k(2)-1))
947set(handles.ZObject,'String',coord_str(k(2)+3:k(3)-1))
948set(handles.XImage,'String',coord_str(k(3)+3:k(4)-1))
949set(handles.YImage,'String',coord_str(k(4)+3:end))
950h_menu_coord=findobj(huvmat,'Tag','transform_fct');
951menu=get(h_menu_coord,'String');
952choice=get(h_menu_coord,'Value');
953if iscell(menu)
954    option=menu{choice};
955else
956    option='px'; %default
957end
958if isequal(option,'phys')
959    XCoord=str2double(coord_str(1:k(1)-1));
960    YCoord=str2double(coord_str(k(1)+3:k(2)-1));
961elseif isequal(option,'px')|| isequal(option,'')
962    XCoord=str2double(coord_str(k(3)+3:k(4)-1));
963    YCoord=str2double(coord_str(k(4)+3:end));
964else
965    msgbox_uvmat('ERROR','the choice in menu_coord of uvmat must be px or phys ')
966end
967if isempty(XCoord)||isempty(YCoord)
968     if ~isempty(hhh)
969        delete(hhh)%delete the circle marker
970    end
971    return
972end
973xlim=get(hplot,'XLim');
974ylim=get(hplot,'YLim');
975ind_range=max(abs(xlim(2)-xlim(1)),abs(ylim(end)-ylim(1)))/20;%defines the size of the circle marker
976if isempty(hhh)
977    set(0,'CurrentFig',huvmat)
978    set(huvmat,'CurrentAxes',hplot)
979    rectangle('Curvature',[1 1],...
980              'Position',[XCoord-ind_range/2 YCoord-ind_range/2 ind_range ind_range],'EdgeColor','m',...
981              'LineStyle','-','Tag','calib_marker');
982else
983    set(hhh,'Position',[XCoord-ind_range/2 YCoord-ind_range/2 ind_range ind_range])
984end
985
986%------------------------------------------------------------------------
987% --- Executes on selection change in edit_append.
988function edit_append_Callback(hObject, eventdata, handles)
989%------------------------------------------------------------------------
990choice=get(handles.edit_append,'Value');
991if choice
992    set(handles.edit_append,'BackgroundColor',[1 1 0])
993    huvmat=findobj(allchild(0),'tag','uvmat');
994    if ishandle(huvmat)
995        hhuvmat=guidata(huvmat);
996        set(hhuvmat.edit_object,'Value',0)
997        set(hhuvmat.edit_object,'BackgroundColor',[0.7 0.7 0.7])
998    end
999else
1000    set(handles.edit_append,'BackgroundColor',[0.7 0.7 0.7])
1001end
1002   
1003function NEW_Callback(hObject, eventdata, handles)
1004%A METTRE SOUS UN BOUTON
1005huvmat=findobj(allchild(0),'Name','uvmat');
1006hchild=get(huvmat,'children');
1007hcoord=findobj(hchild,'Tag','menu_coord');
1008coordtype=get(hcoord,'Value');
1009haxes=findobj(hchild,'Tag','axes3');
1010AxeData=get(haxes,'UserData');
1011if ~isequal(hcoord,2)
1012    set(hcoord,'Value',2)
1013    huvmat=uvmat(AxeData);
1014    'relancer uvmat';
1015end
1016if ~isfield(AxeData,'ZoomAxes')
1017    msgbox_uvmat('ERROR','first draw a window around a grid marker')
1018    return
1019end
1020XLim=get(AxeData.ZoomAxes,'XLim');
1021YLim=get(AxeData.ZoomAxes,'YLim');
1022np=size(AxeData.A);
1023ind_sub_x=round(XLim);
1024ind_sub_y=np(1)-round(YLim);
1025Mfiltre=AxeData.A([ind_sub_y(2):ind_sub_y(1)] ,ind_sub_x,:);
1026Mfiltre_norm=double(Mfiltre);
1027Mfiltre_norm=Mfiltre_norm/sum(sum(Mfiltre_norm));
1028Mfiltre_norm=100*(Mfiltre_norm-mean(mean(Mfiltre_norm)));
1029Atype=class(AxeData.A);
1030Data.NbDim=2;
1031Data.A=filter2(Mfiltre_norm,double(AxeData.A));
1032Data.A=feval(Atype,Data.A);
1033Data.AName='image';
1034Data.AX=AxeData.AX;
1035Data.AY=AxeData.AY;
1036Data.CoordType='px';
1037plot_field(Data)
1038
1039%------------------------------------------------------------------------
1040function MenuPlot_Callback(hObject, eventdata, handles)
1041%------------------------------------------------------------------------
1042huvmat=findobj(allchild(0),'Name','uvmat');%find the current uvmat interface handle
1043%UvData=get(huvmat,'UserData');%Data associated to the current uvmat interface
1044hhuvmat=guidata(huvmat); %handles of GUI elements in uvmat
1045%hplot=findobj(huvmat,'Tag','axes3');%main plotting axis of uvmat
1046h_menu_coord=findobj(huvmat,'Tag','transform_fct');
1047menu=get(h_menu_coord,'String');
1048choice=get(h_menu_coord,'Value');
1049if iscell(menu)
1050    option=menu{choice};
1051else
1052    option='px'; %default
1053end
1054Coord_cell=get(handles.ListCoord,'String');
1055ObjectData=read_geometry_calib(Coord_cell);
1056%ObjectData=read_geometry_calib(handles);%read the interface input parameters defining the object
1057if ~isempty(ObjectData.Coord)
1058    if isequal(option,'phys')
1059        ObjectData.Coord=ObjectData.Coord(:,[1:3]);
1060    elseif isequal(option,'px')||isequal(option,'')
1061        ObjectData.Coord=ObjectData.Coord(:,[4:5]);
1062    else
1063        msgbox_uvmat('ERROR','the choice in menu_coord of uvmat must be '''', px or phys ')
1064    end
1065end
1066axes(hhuvmat.axes3)
1067hh=findobj('Tag','calib_points');
1068if  ~isempty(ObjectData.Coord) && isempty(hh)
1069    hh=line(ObjectData.Coord(:,1),ObjectData.Coord(:,2),'Color','m','Tag','calib_points','LineStyle','.','Marker','+');
1070elseif isempty(ObjectData.Coord)%empty list of points, suppress the plot
1071    delete(hh)
1072else
1073    set(hh,'XData',ObjectData.Coord(:,1))
1074    set(hh,'YData',ObjectData.Coord(:,2))
1075end
1076pause(.1)
1077figure(handles.geometry_calib)
1078
1079% --------------------------------------------------------------------
1080function MenuHelp_Callback(hObject, eventdata, handles)
1081path_to_uvmat=which('uvmat');% check the path of uvmat
1082pathelp=fileparts(path_to_uvmat);
1083helpfile=fullfile(pathelp,'uvmat_doc','uvmat_doc.html');
1084if isempty(dir(helpfile)), msgbox_uvmat('ERROR','Please put the help file uvmat_doc.html in the sub-directory /uvmat_doc of the UVMAT package')
1085else
1086   addpath (fullfile(pathelp,'uvmat_doc'))
1087   web([helpfile '#geometry_calib'])
1088end
1089
1090%------------------------------------------------------------------------
1091function MenuCreateGrid_Callback(hObject, eventdata, handles)
1092%------------------------------------------------------------------------
1093%hcalib=get(handles.calib_type,'parent');%handles of the GUI geometry_calib
1094CalibData=get(handles.geometry_calib,'UserData');
1095Tinput=[];%default
1096if isfield(CalibData,'grid')
1097    Tinput=CalibData.grid;
1098end
1099[T,CalibData.grid]=create_grid(Tinput);%display the GUI create_grid
1100set(handles.geometry_calib,'UserData',CalibData)
1101
1102%grid in phys space
1103Coord=get(handles.ListCoord,'String');
1104val=get(handles.ListCoord,'Value');
1105data=read_geometry_calib(Coord);
1106%nbpoints=size(data.Coord,1); %nbre of calibration points
1107data.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
1108% for i=1:nbpoints
1109%    for j=1:5
1110%           Coord{i,j}=num2str(data.Coord(i,j),4);%display coordiantes with 4 digits
1111%    end
1112% end
1113%update the phys coordinates starting from the selected point (down in the
1114Coord(end,:)=[]; %remove last string '.....'
1115for i=1:size(data.Coord,1)
1116    for j=1:5
1117          Coord{i,j}=num2str(data.Coord(i,j),4);%display coordiantes with 4 digits
1118    end
1119end
1120
1121%size(data.Coord,1)
1122Tabchar=cell2tab(Coord,' | ');
1123Tabchar=[Tabchar ;{'......'}];
1124set(handles.ListCoord,'String',Tabchar)
1125
1126% -----------------------------------------------------------------------
1127% --- automatic grid dectection from local maxima of the images
1128function MenuDetectGrid_Callback(hObject, eventdata, handles)
1129%------------------------------------------------------------------------
1130%% initiate the grid
1131CalibData=get(handles.geometry_calib,'UserData');%get information stored on the GUI geometry_calib
1132grid_input=[];%default
1133if isfield(CalibData,'grid')
1134    grid_input=CalibData.grid;%retrieve the previously used grid
1135end
1136[T,CalibData.grid,white_test]=create_grid(grid_input,'detect_grid');%display the GUI create_grid, read the set of phys coordinates T
1137set(handles.geometry_calib,'UserData',CalibData)%store the phys grid parameters for later use
1138
1139%% read the four last point coordinates in pixels
1140Coord_cell=get(handles.ListCoord,'String');%read list of coordinates on geometry_calib
1141data=read_geometry_calib(Coord_cell);
1142nbpoints=size(data.Coord,1); %nbre of calibration points
1143if nbpoints~=4
1144    msgbox_uvmat('ERROR','four points must have be selected by the mouse, beginning by the new x axis, to delimitate the phys grid area')
1145    return
1146end
1147corners_X=(data.Coord(end:-1:end-3,4)); %pixel absissa of the four corners
1148corners_Y=(data.Coord(end:-1:end-3,5));
1149
1150%reorder the last two points (the two first in the list) if needed
1151angles=angle((corners_X-corners_X(1))+1i*(corners_Y-corners_Y(1)));
1152if abs(angles(4)-angles(2))>abs(angles(3)-angles(2))
1153      X_end=corners_X(4);
1154      Y_end=corners_Y(4);
1155      corners_X(4)=corners_X(3);
1156      corners_Y(4)=corners_Y(3);
1157      corners_X(3)=X_end;
1158      corners_Y(3)=Y_end;
1159end
1160
1161%% read the current image, displayed in the GUI uvmat
1162huvmat=findobj(allchild(0),'Name','uvmat');
1163UvData=get(huvmat,'UserData');
1164A=UvData.Field.A;
1165npxy=size(A);
1166X=[CalibData.grid.x_0 CalibData.grid.x_1 CalibData.grid.x_0 CalibData.grid.x_1]';%corner absissa in the phys coordinates (cm)
1167Y=[CalibData.grid.y_0 CalibData.grid.y_0 CalibData.grid.y_1 CalibData.grid.y_1]';%corner ordinates in the phys coordinates (cm)
1168
1169%calculate transform matrices for plane projection
1170% reference: http://alumni.media.mit.edu/~cwren/interpolator/ by Christopher R. Wren
1171B = [ X Y ones(size(X)) zeros(4,3)        -X.*corners_X -Y.*corners_X ...
1172      zeros(4,3)        X Y ones(size(X)) -X.*corners_Y -Y.*corners_Y ];
1173B = reshape (B', 8 , 8 )';
1174D = [ corners_X , corners_Y ];
1175D = reshape (D', 8 , 1 );
1176l = (B' * B)\B' * D;
1177Amat = reshape([l(1:6)' 0 0 1 ],3,3)';
1178C = [l(7:8)' 1];
1179
1180% transform grid image into 'phys' coordinates
1181GeometryCalib.fx_fy=[1 1];
1182GeometryCalib.Tx_Ty_Tz=[Amat(1,3) Amat(2,3) 1];
1183GeometryCalib.R=[Amat(1,1),Amat(1,2),0;Amat(2,1),Amat(2,2),0;C(1),C(2),0];
1184GeometryCalib.CoordUnit='cm';
1185path_uvmat=which('uvmat');% check the path detected for source file uvmat
1186path_UVMAT=fileparts(path_uvmat); %path to UVMAT
1187addpath(fullfile(path_UVMAT,'transform_field'))
1188Data.ListVarName={'AY','AX','A'};
1189Data.VarDimName={'AY','AX',{'AY','AX'}};
1190if ndims(A)==3
1191    A=mean(A,3);
1192end
1193Data.A=A-min(min(A));
1194Data.AY=[npxy(1)-0.5 0.5];
1195Data.AX=[0.5 npxy(2)];
1196Data.CoordUnit='pixel';
1197Calib.GeometryCalib=GeometryCalib;
1198DataOut=phys(Data,Calib);
1199rmpath(fullfile(path_UVMAT,'transform_field'))
1200Amod=DataOut.A;
1201Rangx=DataOut.AX;
1202Rangy=DataOut.AY;
1203if white_test
1204    Amod=double(Amod);%case of white grid markers: will look for image maxima
1205else
1206    Amod=-double(Amod);%case of black grid markers: will look for image minima
1207end
1208% figure(12) %display corrected image
1209% Amax=max(max(Amod));
1210% image(Rangx,Rangy,uint8(255*Amod/Amax))
1211
1212%% detection of local image extrema in each direction
1213Dx=(Rangx(2)-Rangx(1))/(npxy(2)-1); %x mesh in real space
1214Dy=(Rangy(2)-Rangy(1))/(npxy(1)-1); %y mesh in real space
1215ind_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
1216ind_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
1217nbpoints=size(T,1);
1218for ipoint=1:nbpoints
1219    i0=1+round((T(ipoint,1)-Rangx(1))/Dx);%round(Xpx(ipoint));
1220    j0=1+round((T(ipoint,2)-Rangy(1))/Dy);%round(Xpx(ipoint));
1221    j0min=max(j0-ind_range_y,1);
1222    j0max=min(j0+ind_range_y,size(Amod,1));
1223    i0min=max(i0-ind_range_x,1);
1224    i0max=min(i0+ind_range_x,size(Amod,2));
1225    Asub=Amod(j0min:j0max,i0min:i0max);
1226    x_profile=sum(Asub,1);
1227    y_profile=sum(Asub,2);
1228    [Amax,ind_x_max]=max(x_profile);
1229    [Amax,ind_y_max]=max(y_profile);
1230    %sub-pixel improvement using moments
1231    x_shift=0;
1232    y_shift=0;
1233    if ind_x_max+2<=numel(x_profile) && ind_x_max-2>=1
1234        Atop=x_profile(ind_x_max-2:ind_x_max+2);
1235        x_shift=sum(Atop.*[-2 -1 0 1 2])/sum(Atop);
1236    end
1237    if ind_y_max+2<=numel(y_profile) && ind_y_max-2>=1
1238        Atop=y_profile(ind_y_max-2:ind_y_max+2);
1239        y_shift=sum(Atop.*[-2 -1 0 1 2]')/sum(Atop);
1240    end
1241    Delta(ipoint,1)=(x_shift+ind_x_max+i0min-i0-1)*Dx;%shift from the initial guess
1242    Delta(ipoint,2)=(y_shift+ind_y_max+j0min-j0-1)*Dy;
1243end
1244Tmod=T(:,(1:2))+Delta;
1245[Xpx,Ypx]=px_XYZ(GeometryCalib,Tmod(:,1),Tmod(:,2));
1246for ipoint=1:nbpoints
1247     Coord{ipoint,1}=num2str(T(ipoint,1),4);%display coordiantes with 4 digits
1248     Coord{ipoint,2}=num2str(T(ipoint,2),4);%display coordiantes with 4 digits
1249     Coord{ipoint,3}=num2str(T(ipoint,3),4);%display coordiantes with 4 digits;
1250     Coord{ipoint,4}=num2str(Xpx(ipoint),4);%display coordiantes with 4 digits
1251     Coord{ipoint,5}=num2str(Ypx(ipoint),4);%display coordiantes with 4 digits
1252end
1253Tabchar=cell2tab(Coord(end:-1:1,:),' | ');
1254Tabchar=[Tabchar ;{'......'}];
1255set(handles.ListCoord,'Value',1)
1256set(handles.ListCoord,'String',Tabchar)
1257MenuPlot_Callback(hObject, eventdata, handles)
1258
1259%-----------------------------------------------------------------------
1260function MenuTranslatePoints_Callback(hObject, eventdata, handles)
1261%-----------------------------------------------------------------------
1262%hcalib=get(handles.calib_type,'parent');%handles of the GUI geometry_calib
1263CalibData=get(handles.geometry_calib,'UserData');
1264Tinput=[];%default
1265if isfield(CalibData,'translate')
1266    Tinput=CalibData.translate;
1267end
1268T=translate_points(Tinput);%display translate_points GUI and get shift parameters
1269CalibData.translate=T;
1270set(handles.geometry_calib,'UserData',CalibData)
1271%translation
1272Coord_cell=get(handles.ListCoord,'String');
1273data=read_geometry_calib(Coord_cell);
1274data.Coord(:,1)=T(1)+data.Coord(:,1);
1275data.Coord(:,2)=T(2)+data.Coord(:,2);
1276data.Coord(:,3)=T(3)+data.Coord(:,3);
1277data.Coord(:,[4 5])=data.Coord(:,[4 5]);
1278for i=1:size(data.Coord,1)
1279    for j=1:5
1280          Coord{i,j}=num2str(data.Coord(i,j),4);%phys x,y,z
1281   end
1282end
1283Tabchar=cell2tab(Coord,' | ');
1284Tabchar=[Tabchar; {'.....'}];
1285%set(handles.ListCoord,'Value',1)
1286set(handles.ListCoord,'String',Tabchar)
1287
1288
1289% --------------------------------------------------------------------
1290function MenuRotatePoints_Callback(hObject, eventdata, handles)
1291%hcalib=get(handles.calib_type,'parent');%handles of the GUI geometry_calib
1292CalibData=get(handles.geometry_calib,'UserData');
1293Tinput=[];%default
1294if isfield(CalibData,'rotate')
1295    Tinput=CalibData.rotate;
1296end
1297T=rotate_points(Tinput);%display translate_points GUI and get shift parameters
1298CalibData.rotate=T;
1299set(handles.geometry_calib,'UserData',CalibData)
1300%-----------------------------------------------------
1301%rotation
1302Phi=T(1);
1303O_x=0;%default
1304O_y=0;%default
1305if numel(T)>=2
1306    O_x=T(2);%default
1307end
1308if numel(T)>=3
1309    O_y=T(3);%default
1310end
1311Coord_cell=get(handles.ListCoord,'String');
1312data=read_geometry_calib(Coord_cell);
1313r1=cos(pi*Phi/180);
1314r2=-sin(pi*Phi/180);
1315r3=sin(pi*Phi/180);
1316r4=cos(pi*Phi/180);
1317x=data.Coord(:,1)-O_x;
1318y=data.Coord(:,2)-O_y;
1319data.Coord(:,1)=r1*x+r2*y;
1320data.Coord(:,2)=r3*x+r4*y;
1321% data.Coord(:,[4 5])=data.Coord(:,[4 5]);
1322for i=1:size(data.Coord,1)
1323    for j=1:5
1324          Coord{i,j}=num2str(data.Coord(i,j),4);%phys x,y,z
1325   end
1326end
1327Tabchar=cell2tab(Coord,' | ');
1328Tabchar=[Tabchar;{'......'}];
1329set(handles.ListCoord,'Value',1)
1330set(handles.ListCoord,'String',Tabchar)
1331
1332
1333% %------------------------------------------------------------------------
1334% % --- Executes on button press in rotation.
1335% function rotation_Callback(hObject, eventdata, handles)
1336% %------------------------------------------------------------------------
1337% angle_rot=(pi/180)*str2num(get(handles.Phi,'String'));
1338% Coord_cell=get(handles.ListCoord,'String');
1339% data=read_geometry_calib(Coord_cell);
1340% data.Coord(:,1)=cos(angle_rot)*data.Coord(:,1)+sin(angle_rot)*data.Coord(:,2);
1341% data.Coord(:,1)=-sin(angle_rot)*data.Coord(:,1)+cos(angle_rot)*data.Coord(:,2);
1342% set(handles.XObject,'String',num2str(data.Coord(:,1),4));
1343% set(handles.YObject,'String',num2str(data.Coord(:,2),4));
1344
1345
1346%------------------------------------------------------------------------
1347% image transform from px to phys
1348%INPUT:
1349%Zindex: index of plane
1350% function [A_out,Rangx,Rangy]=phys_Ima(A,Calib,ZIndex)
1351% %------------------------------------------------------------------------
1352% xcorner=[];
1353% ycorner=[];
1354% npx=[];
1355% npy=[];
1356% siz=size(A)
1357% npx=[npx siz(2)];
1358% npy=[npy siz(1)]
1359% xima=[0.5 siz(2)-0.5 0.5 siz(2)-0.5];%image coordinates of corners
1360% yima=[0.5 0.5 siz(1)-0.5 siz(1)-0.5];
1361% [xcorner,ycorner]=phys_XYZ(Calib,xima,yima,ZIndex);%corresponding physical coordinates
1362% Rangx(1)=min(xcorner);
1363% Rangx(2)=max(xcorner);
1364% Rangy(2)=min(ycorner);
1365% Rangy(1)=max(ycorner);
1366% test_multi=(max(npx)~=min(npx)) | (max(npy)~=min(npy));
1367% npx=max(npx);
1368% npy=max(npy);
1369% x=linspace(Rangx(1),Rangx(2),npx);
1370% y=linspace(Rangy(1),Rangy(2),npy);
1371% [X,Y]=meshgrid(x,y);%grid in physical coordiantes
1372% vec_B=[];
1373%
1374% zphys=0; %default
1375% if isfield(Calib,'SliceCoord') %.Z= index of plane
1376%    SliceCoord=Calib.SliceCoord(ZIndex,:);
1377%    zphys=SliceCoord(3); %to generalize for non-parallel planes
1378% end
1379% [XIMA,YIMA]=px_XYZ(Calib,X,Y,zphys);%corresponding image indices for each point in the real space grid
1380% XIMA=reshape(round(XIMA),1,npx*npy);%indices reorganized in 'line'
1381% YIMA=reshape(round(YIMA),1,npx*npy);
1382% flagin=XIMA>=1 & XIMA<=npx & YIMA >=1 & YIMA<=npy;%flagin=1 inside the original image
1383% testuint8=isa(A,'uint8');
1384% testuint16=isa(A,'uint16');
1385% if numel(siz)==2 %(B/W images)
1386%     vec_A=reshape(A,1,npx*npy);%put the original image in line
1387%     ind_in=find(flagin);
1388%     ind_out=find(~flagin);
1389%     ICOMB=((XIMA-1)*npy+(npy+1-YIMA));
1390%     ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A
1391%     vec_B(ind_in)=vec_A(ICOMB);
1392%     vec_B(ind_out)=zeros(size(ind_out));
1393%     A_out=reshape(vec_B,npy,npx);%new image in real coordinates
1394% elseif numel(siz)==3     
1395%     for icolor=1:siz(3)
1396%         vec_A=reshape(A{icell}(:,:,icolor),1,npx*npy);%put the original image in line
1397%         ind_in=find(flagin);
1398%         ind_out=find(~flagin);
1399%         ICOMB=((XIMA-1)*npy+(npy+1-YIMA));
1400%         ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A
1401%         vec_B(ind_in)=vec_A(ICOMB);
1402%         vec_B(ind_out)=zeros(size(ind_out));
1403%         A_out(:,:,icolor)=reshape(vec_B,npy,npx);%new image in real coordinates
1404%     end
1405% end
1406% if testuint8
1407%     A_out=uint8(A_out);
1408% end
1409% if testuint16
1410%     A_out=uint16(A_out);
1411% end
1412
1413%------------------------------------------------------------------------
1414% pointwise transform from px to phys
1415%INPUT:
1416%Z: index of plane
1417% function [Xphys,Yphys,Zphys]=phys_XYZ(Calib,X,Y,Z)
1418% %------------------------------------------------------------------------
1419% if exist('Z','var')& isequal(Z,round(Z))& Z>0 & isfield(Calib,'SliceCoord')&length(Calib.SliceCoord)>=Z
1420%     Zindex=Z;
1421%     Zphys=Calib.SliceCoord(Zindex,3);%GENERALISER AUX CAS AVEC ANGLE
1422% else
1423%     Zphys=0;
1424% end
1425% if ~exist('X','var')||~exist('Y','var')
1426%     Xphys=[];
1427%     Yphys=[];%default
1428%     return
1429% end
1430% Xphys=X;%default
1431% Yphys=Y;
1432% %image transform
1433% if isfield(Calib,'R')
1434%     R=(Calib.R)';
1435%     Dx=R(5)*R(7)-R(4)*R(8);
1436%     Dy=R(1)*R(8)-R(2)*R(7);
1437%     D0=Calib.f*(R(2)*R(4)-R(1)*R(5));
1438%     Z11=R(6)*R(8)-R(5)*R(9);
1439%     Z12=R(2)*R(9)-R(3)*R(8); 
1440%     Z21=R(4)*R(9)-R(6)*R(7);
1441%     Z22=R(3)*R(7)-R(1)*R(9);
1442%     Zx0=R(3)*R(5)-R(2)*R(6);
1443%     Zy0=R(1)*R(6)-R(3)*R(4);
1444%     A11=R(8)*Calib.Ty-R(5)*Calib.Tz+Z11*Zphys;
1445%     A12=R(2)*Calib.Tz-R(8)*Calib.Tx+Z12*Zphys;
1446%     A21=-R(7)*Calib.Ty+R(4)*Calib.Tz+Z21*Zphys;
1447%     A22=-R(1)*Calib.Tz+R(7)*Calib.Tx+Z11*Zphys;
1448%     X0=Calib.f*(R(5)*Calib.Tx-R(2)*Calib.Ty+Zx0*Zphys);
1449%     Y0=Calib.f*(-R(4)*Calib.Tx+R(1)*Calib.Ty+Zy0*Zphys);
1450%         %px to camera:
1451%     Xd=(Calib.dpx/Calib.sx)*(X-Calib.Cx); % sensor coordinates
1452%     Yd=Calib.dpy*(Y-Calib.Cy);
1453%     dist_fact=1+Calib.kappa1*(Xd.*Xd+Yd.*Yd); %distortion factor
1454%     Xu=dist_fact.*Xd;%undistorted sensor coordinates
1455%     Yu=dist_fact.*Yd;
1456%     denom=Dx*Xu+Dy*Yu+D0;
1457%     % denom2=denom.*denom;
1458%     Xphys=(A11.*Xu+A12.*Yu+X0)./denom;%world coordinates
1459%     Yphys=(A21.*Xu+A22.*Yu+Y0)./denom;
1460% end
1461
1462
1463% --------------------------------------------------------------------
1464function MenuImportPoints_Callback(hObject, eventdata, handles)
1465fileinput=browse_xml(hObject, eventdata, handles);
1466if isempty(fileinput)
1467    return
1468end
1469[s,errormsg]=imadoc2struct(fileinput,'GeometryCalib');
1470GeometryCalib=s.GeometryCalib;
1471%GeometryCalib=load_calib(hObject, eventdata, handles)
1472calib=reshape(GeometryCalib.PointCoord,[],1);
1473for ilist=1:numel(calib)
1474    CoordCell{ilist}=num2str(calib(ilist));
1475end
1476CoordCell=reshape(CoordCell,[],5);
1477Tabchar=cell2tab(CoordCell,' | ');%transform cells into table ready for display
1478Tabchar=[Tabchar;{'......'}];
1479set(handles.ListCoord,'Value',1)
1480set(handles.ListCoord,'String',Tabchar)
1481MenuPlot_Callback(handles.geometry_calib, [], handles)
1482
1483% -----------------------------------------------------------------------
1484function MenuImportIntrinsic_Callback(hObject, eventdata, handles)
1485%------------------------------------------------------------------------
1486fileinput=browse_xml(hObject, eventdata, handles);
1487if isempty(fileinput)
1488    return
1489end
1490[s,errormsg]=imadoc2struct(fileinput,'GeometryCalib');
1491GeometryCalib=s.GeometryCalib;
1492display_intrinsic(GeometryCalib,handles)
1493
1494% -----------------------------------------------------------------------
1495function MenuImportAll_Callback(hObject, eventdata, handles)
1496%------------------------------------------------------------------------
1497fileinput=browse_xml(hObject, eventdata, handles);
1498if ~isempty(fileinput)
1499    loadfile(handles,fileinput)
1500end
1501
1502% -----------------------------------------------------------------------
1503% --- Executes on menubar option Import/Grid file: introduce previous grid files
1504function MenuGridFile_Callback(hObject, eventdata, handles)
1505% -----------------------------------------------------------------------
1506inputfile=browse_xml(hObject, eventdata, handles);
1507listfile=get(handles.coord_files,'string');
1508if isequal(listfile,{''})
1509    listfile={inputfile};
1510else
1511    listfile=[listfile;{inputfile}];%update the list of coord files
1512end
1513set(handles.coord_files,'string',listfile);
1514
1515%------------------------------------------------------------------------
1516% --- 'key_press_fcn:' function activated when a key is pressed on the keyboard
1517function key_press_fcn(hObject,eventdata,handles)
1518%------------------------------------------------------------------------
1519xx=double(get(handles.geometry_calib,'CurrentCharacter')); %get the keyboard character
1520if ismember(xx,[8 127])%backspace or delete
1521    Coord_cell=get(handles.ListCoord,'String');
1522    val=get(handles.ListCoord,'Value');
1523     if max(val)<numel(Coord_cell) % the last element '...' has not been selected
1524        Coord_cell(val)=[];%remove the selected line
1525        set(handles.ListCoord,'Value',min(val))
1526        set(handles.ListCoord,'String',Coord_cell)         
1527        ListCoord_Callback(hObject, eventdata, handles)
1528        MenuPlot_Callback(hObject,eventdata,handles)
1529     end
1530end
1531
1532%------------------------------------------------------------------------
1533function fileinput=browse_xml(hObject, eventdata, handles)
1534%------------------------------------------------------------------------
1535fileinput=[];%default
1536oldfile=''; %default
1537UserData=get(handles.geometry_calib,'UserData');
1538if isfield(UserData,'XmlInputFile')
1539    oldfile=UserData.XmlInputFile;
1540end
1541[FileName, PathName, filterindex] = uigetfile( ...
1542       {'*.xml;*.mat', ' (*.xml,*.mat)';
1543       '*.xml',  '.xml files '; ...
1544        '*.mat',  '.mat matlab files '}, ...
1545        'Pick a file',oldfile);
1546fileinput=[PathName FileName];%complete file name
1547testblank=findstr(fileinput,' ');%look for blanks
1548if ~isempty(testblank)
1549    msgbox_uvmat('ERROR','forbidden input file name or path: no blank character allowed')
1550    return
1551end
1552sizf=size(fileinput);
1553if (~ischar(fileinput)||~isequal(sizf(1),1)),return;end
1554UserData.XmlInputFile=fileinput;
1555set(handles.geometry_calib,'UserData',UserData)%record current file foer further use of browser
1556
1557% -----------------------------------------------------------------------
1558function Heading=loadfile(handles,fileinput)
1559%------------------------------------------------------------------------
1560Heading=[];%default
1561[s,errormsg]=imadoc2struct(fileinput,'GeometryCalib');
1562if ~isempty(errormsg)
1563    msgbox_uvmat('ERROR',['Error for reading ' fileinput ': '  errormsg])
1564    return
1565end
1566if ~isempty(s.Heading)
1567    Heading=s.Heading;
1568end
1569   
1570GeometryCalib=s.GeometryCalib;
1571fx=1;fy=1;Cx=0;Cy=0;kc=0; %default
1572%     Tabchar={};
1573CoordCell={};
1574%     kc=0;%default
1575%     f1=1000;
1576%     f2=1000;
1577%     hhuvmat=guidata(findobj(allchild(0),'Name','uvmat'));
1578%     Cx=str2num(get(hhuvmat.npx,'String'))/2;
1579%     Cy=str2num(get(hhuvmat.npy,'String'))/2;
1580Tabchar={};%default
1581val_cal=1;%default
1582if ~isempty(GeometryCalib)
1583    % choose the calibration option
1584    if isfield(GeometryCalib,'CalibrationType')
1585       calib_list=get(handles.calib_type,'String');
1586       for ilist=1:numel(calib_list)
1587           if strcmp(calib_list{ilist},GeometryCalib.CalibrationType)
1588               val_cal=ilist;
1589               break
1590           end
1591       end
1592    end
1593    display_intrinsic(GeometryCalib,handles)%intrinsic param
1594    %extrinsic param
1595    if isfield(GeometryCalib,'Tx_Ty_Tz')
1596        Tx_Ty_Tz=GeometryCalib.Tx_Ty_Tz;
1597        set(handles.Tx,'String',num2str(GeometryCalib.Tx_Ty_Tz(1),4))
1598        set(handles.Ty,'String',num2str(GeometryCalib.Tx_Ty_Tz(2),4))
1599        set(handles.Tz,'String',num2str(GeometryCalib.Tx_Ty_Tz(3),4))
1600    end
1601    if isfield(GeometryCalib,'omc')
1602        set(handles.Phi,'String',num2str(GeometryCalib.omc(1),4))
1603        set(handles.Theta,'String',num2str(GeometryCalib.omc(2),4))
1604        set(handles.Psi,'String',num2str(GeometryCalib.omc(3),4))
1605    end
1606    calib=reshape(GeometryCalib.PointCoord,[],1);
1607    for ilist=1:numel(calib)
1608        CoordCell{ilist}=num2str(calib(ilist));
1609    end
1610    CoordCell=reshape(CoordCell,[],5);
1611    Tabchar=cell2tab(CoordCell,' | ');%transform cells into table ready for display
1612    MenuPlot_Callback(handles.geometry_calib, [], handles)
1613end
1614set(handles.calib_type,'Value',val_cal)
1615Tabchar=[Tabchar;{'......'}];
1616set(handles.ListCoord,'Value',1)
1617set(handles.ListCoord,'String',Tabchar)
1618
1619if isempty(CoordCell)% allow mouse action by default in the absence of input points
1620    set(handles.edit_append,'Value',1)
1621    set(handles.edit_append,'BackgroundColor',[1 1 0])
1622else % does not allow mouse action by default in the presence of input points
1623    set(handles.edit_append,'Value',0)
1624    set(handles.edit_append,'BackgroundColor',[0.7 0.7 0.7])
1625end
1626
1627%------------------------------------------------------------------------
1628%---display calibration intrinsic parameters
1629function display_intrinsic(GeometryCalib,handles)
1630%------------------------------------------------------------------------
1631fx=[];
1632fy=[];
1633if isfield(GeometryCalib,'fx_fy')
1634    fx=GeometryCalib.fx_fy(1);
1635    fy=GeometryCalib.fx_fy(2);
1636end
1637Cx_Cy=[0 0];%default
1638if isfield(GeometryCalib,'Cx_Cy')
1639    Cx_Cy=GeometryCalib.Cx_Cy;
1640end
1641kc=0;
1642if isfield(GeometryCalib,'kc')
1643    kc=GeometryCalib.kc; %* GeometryCalib.focal*GeometryCalib.focal;
1644end
1645set(handles.fx,'String',num2str(fx,5))
1646set(handles.fy,'String',num2str(fy,5))
1647set(handles.Cx,'String',num2str(Cx_Cy(1),'%1.1f'))
1648set(handles.Cy,'String',num2str(Cx_Cy(2),'%1.1f'))
1649set(handles.kc,'String',num2str(kc,'%1.4f'))
1650
1651
Note: See TracBrowser for help on using the repository browser.