source: trunk/src/geometry_calib.m @ 494

Last change on this file since 494 was 494, checked in by sommeria, 12 years ago

various bugs corrected after testing in Windows OS. Introduction
of filter tps

File size: 56.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})%
89set(hObject,'WindowButtonDownFcn',{'mouse_alt_gui',handles}) % allows mouse action with right button (zoom for uicontrol display)
90
91%set the position of the interface
92if exist('pos','var')&& length(pos)>=4
93    set(hObject,'Position',pos);
94end
95
96%set menu of calibration options
97set(handles.calib_type,'String',{'rescale';'linear';'3D_linear';'3D_quadr';'3D_extrinsic'})
98if exist('inputfile','var')&& ~isempty(inputfile)
99    struct.XmlInputFile=inputfile;
100    [RootPath,SubDir,RootFile,tild,tild,tild,tild,FileExt]=fileparts_uvmat(inputfile);
101    if ~strcmp(FileExt,'.xml')
102        inputfile=fullfile(RootPath,[SubDir '.xml']);%xml file corresponding to the input file
103        if ~exist(inputfile,'file')% case of civ files , removes the extension for subdir
104            inputfile=fullfile(RootPath,[regexprep(SubDir,'\..+$','') '.xml']);
105            if ~exist(inputfile,'file')
106                inputfile=[fullfile(RootPath,SubDir,RootFile) '.xml'];%old convention
107                if ~exist(inputfile,'file')
108                    inputfile='';
109                end
110            end
111        end
112    end
113    set(handles.ListCoord,'String',{'......'})
114    if exist(inputfile,'file')
115        Heading=loadfile(handles,inputfile);% load the point coordinates existing in the xml file
116        if isfield(Heading,'Campaign')&& ischar(Heading.Campaign)
117            struct.Campaign=Heading.Campaign;
118        end
119    end   
120    set(hObject,'UserData',struct)
121end
122
123set(handles.ListCoord,'KeyPressFcn',{@key_press_fcn,handles})%set keyboard action function
124
125
126%------------------------------------------------------------------------
127% --- Outputs from this function are returned to the command line.
128function varargout = geometry_calib_OutputFcn(hObject, eventdata, handles)
129%------------------------------------------------------------------------
130% Get default command line output from handles structure
131varargout{1} = handles.output;
132varargout{2}=handles;
133%
134%------------------------------------------------------------------------
135% executed when closing: set the parent interface button to value 0
136function closefcn(gcbo,eventdata)
137%------------------------------------------------------------------------
138huvmat=findobj(allchild(0),'Name','uvmat');
139if ~isempty(huvmat)
140    handles=guidata(huvmat);
141    hobject=findobj(handles.axes3,'tag','calib_points');
142    if ~isempty(hobject)
143        delete(hobject)
144    end
145    hobject=findobj(handles.axes3,'tag','calib_marker');
146    if ~isempty(hobject)
147        delete(hobject)
148    end   
149end
150
151%------------------------------------------------------------------------
152% --- Executes on button press in calibrate_lin.
153function APPLY_Callback(hObject, eventdata, handles)
154%------------------------------------------------------------------------
155%% look for the GUI uvmat and check for an image as input
156huvmat=findobj(allchild(0),'Name','uvmat');
157hhuvmat=guidata(huvmat);%handles of elements in the GUI uvmat
158FileExt=get(hhuvmat.FileExt,'String');
159check_input=0;
160if ~isempty(FileExt)
161    if ~isempty(imformats(FileExt(2:end))) ||strcmpi(FileExt,'.avi')
162        check_input=1;
163    end
164end
165if ~check_input
166    msgbox_uvmat('ERROR','open an image with uvmat to perform calibration')
167    return
168end
169
170%% read the current calibration points
171Coord_cell=get(handles.ListCoord,'String');
172Object=read_geometry_calib(Coord_cell);
173Coord=Object.Coord;
174% apply the calibration, whose type is selected in  handles.calib_type
175if ~isempty(Coord)
176    calib_cell=get(handles.calib_type,'String');
177    val=get(handles.calib_type,'Value');
178    GeometryCalib=feval(['calib_' calib_cell{val}],Coord,handles);
179else
180    msgbox_uvmat('ERROR','No calibration points, abort')
181    return
182end
183Z_plane=[];
184if ~isempty(Coord)
185    %check error
186    X=Coord(:,1);
187    Y=Coord(:,2);
188    Z=Coord(:,3);
189    x_ima=Coord(:,4);
190    y_ima=Coord(:,5);
191    [Xpoints,Ypoints]=px_XYZ(GeometryCalib,X,Y,Z);
192    GeometryCalib.ErrorRms(1)=sqrt(mean((Xpoints-x_ima).*(Xpoints-x_ima)));
193    [GeometryCalib.ErrorMax(1),index(1)]=max(abs(Xpoints-x_ima));
194    GeometryCalib.ErrorRms(2)=sqrt(mean((Ypoints-y_ima).*(Ypoints-y_ima)));
195    [GeometryCalib.ErrorMax(2),index(2)]=max(abs(Ypoints-y_ima));
196    [tild,ind_dim]=max(GeometryCalib.ErrorMax);
197    index=index(ind_dim);
198    %set the Z position of the reference plane used for calibration
199    if isequal(max(Z),min(Z))%Z constant
200        Z_plane=Z(1);
201        GeometryCalib.NbSlice=1;
202        GeometryCalib.SliceCoord=[0 0 Z_plane];
203    end
204end
205%set the coordinate unit
206unitlist=get(handles.CoordUnit,'String');
207unit=unitlist{get(handles.CoordUnit,'value')};
208GeometryCalib.CoordUnit=unit;
209%record the points
210GeometryCalib.SourceCalib.PointCoord=Coord;
211display_intrinsic(GeometryCalib,handles)%display calibration intrinsic parameters
212
213% Display extrinsinc parameters (rotation and translation of camera with  respect to the phys coordiantes)
214set(handles.Tx,'String',num2str(GeometryCalib.Tx_Ty_Tz(1),4))
215set(handles.Ty,'String',num2str(GeometryCalib.Tx_Ty_Tz(2),4))
216set(handles.Tz,'String',num2str(GeometryCalib.Tx_Ty_Tz(3),4))
217set(handles.Phi,'String',num2str(GeometryCalib.omc(1),4))
218set(handles.Theta,'String',num2str(GeometryCalib.omc(2),4))
219set(handles.Psi,'String',num2str(GeometryCalib.omc(3),4))
220
221%% store the calibration data, by default in the xml file of the currently displayed image
222UvData=get(huvmat,'UserData');
223NbSlice_j=1;%default
224ZStart=Z_plane;
225ZEnd=Z_plane;
226volume_scan='n';
227if isfield(UvData,'XmlData')
228    if isfield(UvData.XmlData,'TranslationMotor')
229        NbSlice_j=UvData.XmlData.TranslationMotor.Nbslice;
230        ZStart=UvData.XmlData.TranslationMotor.ZStart/10;
231        ZEnd=UvData.XmlData.TranslationMotor.ZEnd/10;
232        volume_scan='y';
233    end
234end
235RootPath='';
236% RootFile='';
237if ~isempty(hhuvmat.RootPath)&& ~isempty(hhuvmat.RootFile)
238%     testhandle=1;
239    RootPath=get(hhuvmat.RootPath,'String');
240    SubDir=get(hhuvmat.SubDir,'String');
241%     RootFile=get(hhuvmat.RootFile,'String');
242%     filebase=fullfile(RootPath,RootFile);
243    outputfile=[fullfile(RootPath,SubDir) '.xml'];%xml file associated with the currently displayed image
244else
245    question={'save the calibration data and point coordinates in'};
246    def={fullfile(RootPath,'ObjectCalib.xml')};
247    options.Resize='on';
248    answer=inputdlg(question,'save average in a new file',1,def,options);
249    outputfile=answer{1};
250end
251answer=msgbox_uvmat('INPUT_Y-N',{[outputfile ' updated with calibration data'];...
252    ['Error rms (along x,y)=' num2str(GeometryCalib.ErrorRms) ' pixels'];...
253    ['Error max (along x,y)=' num2str(GeometryCalib.ErrorMax) ' pixels']});
254
255%% record the calibration parameters and display the current image of uvmat in the new phys coordinates
256if strcmp(answer,'Yes')
257    if strcmp(calib_cell{val}(1:2),'3D')%set the plane position for 3D (projection) calibration
258       input_key={'Z (first position)','Z (last position)','Z (water surface)', 'refractive index','NbSlice','volume scan (y/n)','tilt angle y axis','tilt angle x axis'};
259       input_val=[{num2str(ZEnd)} {num2str(ZStart)} {num2str(ZStart)} {'1.333'} num2str(NbSlice_j) {volume_scan} {'0'} {'0'}];
260        answer=inputdlg(input_key,'slice position(s)',ones(1,8), input_val,'on');
261        %answer_1=msgbox_uvmat('INPUT_TXT',' Z= ',num2str(Z_plane));
262        GeometryCalib.NbSlice=str2double(answer{5});
263        GeometryCalib.VolumeScan=answer{6};
264        if isempty(answer)
265            Z_plane=0; %default
266        else
267            Z_plane=linspace(str2double(answer{1}),str2double(answer{2}),GeometryCalib.NbSlice);
268        end     
269        GeometryCalib.SliceCoord=Z_plane'*[0 0 1];
270        GeometryCalib.SliceAngle(:,3)=0;
271        GeometryCalib.SliceAngle(:,2)=str2double(answer{7})*ones(GeometryCalib.NbSlice,1);%rotation around y axis (to generalise)
272        GeometryCalib.SliceAngle(:,1)=str2double(answer{8})*ones(GeometryCalib.NbSlice,1);%rotation around x axis (to generalise)
273        GeometryCalib.InterfaceCoord=[0 0 str2double(answer{3})];
274        GeometryCalib.RefractionIndex=str2double(answer{4});     
275    end
276    errormsg=update_imadoc(GeometryCalib,outputfile);% introduce the calibration data in the xml file
277    if ~strcmp(errormsg,'')
278        msgbox_uvmat('ERROR',errormsg);
279    end
280   
281    %display image with new calibration in the currently opened uvmat interface
282    hhh=findobj(hhuvmat.axes3,'Tag','calib_marker');% delete calib points and markers
283    if ~isempty(hhh)
284        delete(hhh);     
285    end
286    hhh=findobj(hhuvmat.axes3,'Tag','calib_points');
287    if ~isempty(hhh)
288        delete(hhh);
289    end
290    set(hhuvmat.CheckFixLimits,'Value',0)% put FixedLimits option to 'off'
291    set(hhuvmat.CheckFixLimits,'BackgroundColor',[0.7 0.7 0.7])
292    UserData=get(handles.geometry_calib,'UserData');
293    UserData.XmlInputFile=outputfile;%save the current xml file name
294    set(handles.geometry_calib,'UserData',UserData)
295    uvmat('RootPath_Callback',hObject,eventdata,hhuvmat); %file input with xml reading  in uvmat, show the image in phys coordinates
296    MenuPlot_Callback(hObject, eventdata, handles)
297    set(handles.ListCoord,'Value',index)% indicate in the list the point with max deviation (possible mistake)
298    ListCoord_Callback(hObject, eventdata, handles)
299    figure(handles.geometry_calib)
300end
301
302%------------------------------------------------------------------
303% --- Executes on button press in calibrate_lin.
304
305function REPLICATE_Callback(hObject, eventdata, handles)
306%------------------------------------------------------------------------
307
308%% Apply calibration
309calib_cell=get(handles.calib_type,'String'); %#ok<NASGU>
310val=get(handles.calib_type,'Value'); %#ok<NASGU>
311
312%read the current calibration points
313Coord_cell=get(handles.ListCoord,'String');
314Object=read_geometry_calib(Coord_cell);
315Coord=Object.Coord;
316
317% apply the calibration, whose type is selected in  handles.calib_type
318if ~isempty(Coord)
319    calib_cell=get(handles.calib_type,'String');
320    val=get(handles.calib_type,'Value');
321    GeometryCalib=feval(['calib_' calib_cell{val}],Coord,handles);
322else
323    msgbox_uvmat('ERROR','No calibration points, abort')
324    return
325end
326
327if ~isempty(Coord)
328    %check error
329    X=Coord(:,1);
330    Y=Coord(:,2);
331    Z=Coord(:,3);
332    x_ima=Coord(:,4);
333    y_ima=Coord(:,5);
334    [Xpoints,Ypoints]=px_XYZ(GeometryCalib,X,Y,Z);
335    GeometryCalib.ErrorRms(1)=sqrt(mean((Xpoints-x_ima).*(Xpoints-x_ima)));
336    [GeometryCalib.ErrorMax(1)]=max(abs(Xpoints-x_ima));
337    GeometryCalib.ErrorRms(2)=sqrt(mean((Ypoints-y_ima).*(Ypoints-y_ima)));
338    [GeometryCalib.ErrorMax(2)]=max(abs(Ypoints-y_ima));
339%     [EM,ind_dim]=max(GeometryCalib.ErrorMax);
340    %set the Z position of the reference plane used for calibration
341    Z_plane=[];
342    if isequal(max(Z),min(Z))
343        Z_plane=Z(1);
344    end
345    answer_1=msgbox_uvmat('INPUT_TXT',' Z= ',num2str(Z_plane));
346    Z_plane=str2double(answer_1);
347    GeometryCalib.NbSlice=1;
348    GeometryCalib.SliceCoord=[0 0 Z_plane];
349    %set the coordinate unit
350    unitlist=get(handles.CoordUnit,'String');
351    unit=unitlist{get(handles.CoordUnit,'value')};
352    GeometryCalib.CoordUnit=unit;
353    %record the points
354    GeometryCalib.SourceCalib.PointCoord=Coord;
355end
356
357%% display calibration paprameters
358display_intrinsic(GeometryCalib,handles)%display calibration intrinsic parameters
359
360% Display extrinsinc parameters (rotation and translation of camera with  respect to the phys coordiantes)
361set(handles.Tx,'String',num2str(GeometryCalib.Tx_Ty_Tz(1),4))
362set(handles.Ty,'String',num2str(GeometryCalib.Tx_Ty_Tz(2),4))
363set(handles.Tz,'String',num2str(GeometryCalib.Tx_Ty_Tz(3),4))
364set(handles.Phi,'String',num2str(GeometryCalib.omc(1),4))
365set(handles.Theta,'String',num2str(GeometryCalib.omc(2),4))
366set(handles.Psi,'String',num2str(GeometryCalib.omc(3),4))
367
368%% open the GUI dataview
369h_dataview=findobj(allchild(0),'name','dataview');
370if ~isempty(h_dataview)
371    delete(h_dataview)
372end
373CalibData=get(handles.geometry_calib,'UserData');%read the calibration image source on the interface userdata
374% InputFile='';
375if isfield(CalibData,'XmlInputFile')
376    InputDir=fileparts(CalibData.XmlInputFile);
377    [InputDir,DirName]=fileparts(InputDir);
378end
379SubCampaignTest='n'; %default
380testup=0;
381if isfield(CalibData,'SubCampaign')
382    SubCampaignTest='y';
383    dir_ref=CalibData.SubCampaign;
384    testup=1;
385elseif isfield(CalibData,'Campaign')
386    dir_ref=CalibData.Campaign;
387    testup=1;
388end
389while testup
390    [InputDir,DirName]=fileparts(InputDir);
391    if strcmp(DirName,dir_ref)
392        break
393    end
394end
395InputDir=fullfile(InputDir,DirName);
396answer=msgbox_uvmat('INPUT_TXT','Campaign ?',InputDir);
397if strcmp(answer,'Cancel')
398    return
399end
400
401dataview(answer,SubCampaignTest,GeometryCalib);
402
403%------------------------------------------------------------------------
404% determine the parameters for a calibration by an affine function (rescaling and offset, no rotation)
405function GeometryCalib=calib_rescale(Coord,handles)
406%------------------------------------------------------------------------
407X=Coord(:,1);
408Y=Coord(:,2);% Z not used
409x_ima=Coord(:,4);
410y_ima=Coord(:,5);
411[px]=polyfit(X,x_ima,1);
412[py]=polyfit(Y,y_ima,1);
413% T_x=px(2);
414% T_y=py(2);
415GeometryCalib.CalibrationType='rescale';
416GeometryCalib.fx_fy=[px(1) py(1)];%.fx_fy corresponds to pxcm along x and y
417GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function
418GeometryCalib.Tx_Ty_Tz=[px(2)/px(1) py(2)/py(1) 1];
419GeometryCalib.omc=[0 0 0];
420
421%------------------------------------------------------------------------
422% determine the parameters for a calibration by a linear transform matrix (rescale and rotation)
423
424
425function GeometryCalib=calib_linear(Coord,handles)
426%------------------------------------------------------------------------
427X=Coord(:,1);
428Y=Coord(:,2);% Z not used
429x_ima=Coord(:,4);
430y_ima=Coord(:,5);
431XY_mat=[ones(size(X)) X Y];
432a_X1=XY_mat\x_ima; %transformation matrix for X
433a_Y1=XY_mat\y_ima;%transformation matrix for X
434R=[a_X1(2),a_X1(3);a_Y1(2),a_Y1(3)];
435epsilon=sign(det(R));
436norm=abs(det(R));
437GeometryCalib.CalibrationType='linear';
438if (a_X1(2)/a_Y1(3))>0
439    GeometryCalib.fx_fy(1)=sqrt((a_X1(2)/a_Y1(3))*norm);
440else
441    GeometryCalib.fx_fy(1)=-sqrt(-(a_X1(2)/a_Y1(3))*norm);
442end
443GeometryCalib.fx_fy(2)=(a_Y1(3)/a_X1(2))*GeometryCalib.fx_fy(1);
444GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function
445GeometryCalib.Tx_Ty_Tz=[a_X1(1)/GeometryCalib.fx_fy(1) a_Y1(1)/GeometryCalib.fx_fy(2) 1];
446R(1,:)=R(1,:)/GeometryCalib.fx_fy(1);
447R(2,:)=R(2,:)/GeometryCalib.fx_fy(2);
448R=[R;[0 0]];
449GeometryCalib.R=[R [0;0;-epsilon]];
450GeometryCalib.omc=(180/pi)*[acos(GeometryCalib.R(1,1)) 0 0];
451%------------------------------------------------------------------------
452% determine the tsai parameters for a view normal to the grid plane
453% NOT USED
454function GeometryCalib=calib_normal(Coord,handles)
455%------------------------------------------------------------------------
456Calib.f1=str2num(get(handles.fx,'String'));
457Calib.f2=str2num(get(handles.fy,'String'));
458Calib.k=str2num(get(handles.kc,'String'));
459Calib.Cx=str2num(get(handles.Cx,'String'));
460Calib.Cy=str2num(get(handles.Cy,'String'));
461%default
462if isempty(Calib.f1)
463    Calib.f1=25/0.012;
464end
465if isempty(Calib.f2)
466    Calib.f2=25/0.012;
467end
468if isempty(Calib.k)
469    Calib.k=0;
470end
471if isempty(Calib.Cx)||isempty(Calib.Cy)
472    huvmat=findobj(allchild(0),'Tag','uvmat');
473    hhuvmat=guidata(huvmat);
474    Calib.Cx=str2num(get(hhuvmat.num_Npx,'String'))/2;
475    Calib.Cx=str2num(get(hhuvmat.num_Npy,'String'))/2;
476end   
477%tsai parameters
478Calib.dpx=0.012;%arbitrary
479Calib.dpy=0.012;
480Calib.sx=Calib.f1*Calib.dpx/(Calib.f2*Calib.dpy);
481Calib.f=Calib.f2*Calib.dpy;
482Calib.kappa1=Calib.k/(Calib.f*Calib.f);
483
484%initial guess
485X=Coord(:,1);
486Y=Coord(:,2);
487Zmean=mean(Coord(:,3));
488x_ima=Coord(:,4)-Calib.Cx;
489y_ima=Coord(:,5)-Calib.Cy;
490XY_mat=[ones(size(X)) X Y];
491a_X1=XY_mat\x_ima; %transformation matrix for X
492a_Y1=XY_mat\y_ima;%transformation matrix for Y
493R=[a_X1(2),a_X1(3),0;a_Y1(2),a_Y1(3),0;0,0,-1];% rotation+ z axis reversal (upward)
494norm=sqrt(det(-R));
495calib_param(1)=0;% quadratic distortion
496calib_param(2)=a_X1(1);
497calib_param(3)=a_Y1(1);
498calib_param(4)=Calib.f/(norm*Calib.dpx)-R(3,3)*Zmean;
499calib_param(5)=angle(a_X1(2)+1i*a_X1(3));
500display(['initial guess=' num2str(calib_param)])
501
502%optimise the parameters: minimisation of error
503calib_param = fminsearch(@(calib_param) error_calib(calib_param,Calib,Coord),calib_param);
504
505GeometryCalib.CalibrationType='tsai_normal';
506GeometryCalib.focal=Calib.f;
507GeometryCalib.dpx_dpy=[Calib.dpx Calib.dpy];
508GeometryCalib.Cx_Cy=[Calib.Cx Calib.Cy];
509GeometryCalib.sx=Calib.sx;
510GeometryCalib.kappa1=calib_param(1);
511GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function
512GeometryCalib.Tx_Ty_Tz=[calib_param(2) calib_param(3) calib_param(4)];
513alpha=calib_param(5);
514GeometryCalib.R=[cos(alpha) sin(alpha) 0;-sin(alpha) cos(alpha) 0;0 0 -1];
515
516%------------------------------------------------------------------------
517function GeometryCalib=calib_3D_linear(Coord,handles)
518%------------------------------------------------------------------
519path_uvmat=which('uvmat');% check the path detected for source file uvmat
520path_UVMAT=fileparts(path_uvmat); %path to UVMAT
521huvmat=findobj(allchild(0),'Tag','uvmat');
522hhuvmat=guidata(huvmat);
523coord_files=get(handles.coord_files,'String');
524if ischar(coord_files)
525    coord_files={coord_files};
526end
527if isempty(coord_files{1}) || isequal(coord_files,{''})
528    coord_files={};
529end
530%retrieve the calibration points stored in the files listed in the popup list coord_files
531x_1=Coord(:,4:5)';%px coordinates of the ref points
532nx=str2num(get(hhuvmat.num_Npx,'String'));
533ny=str2num(get(hhuvmat.num_Npy,'String'));
534x_1(2,:)=ny-x_1(2,:);%reverse the y image coordinates
535X_1=Coord(:,1:3)';%phys coordinates of the ref points
536n_ima=numel(coord_files)+1;
537if ~isempty(coord_files)
538    msgbox_uvmat('CONFIRMATION',['The xy coordinates of the calibration points in ' num2str(n_ima) ' planes will be used'])
539    for ifile=1:numel(coord_files)
540    t=xmltree(coord_files{ifile});
541    s=convert(t);%convert to matlab structure
542        if isfield(s,'GeometryCalib')
543            if isfield(s.GeometryCalib,'SourceCalib')
544                if isfield(s.GeometryCalib.SourceCalib,'PointCoord')
545                PointCoord=s.GeometryCalib.SourceCalib.PointCoord;
546                Coord_file=zeros(length(PointCoord),5);%default
547                for i=1:length(PointCoord)
548                    line=str2num(PointCoord{i});
549                    Coord_file(i,4:5)=line(4:5);%px x
550                    Coord_file(i,1:3)=line(1:3);%phys x
551                end
552                eval(['x_' num2str(ifile+1) '=Coord_file(:,4:5)'';']);
553                eval(['x_' num2str(ifile+1) '(2,:)=ny-x_' num2str(ifile+1) '(2,:);' ]);
554                eval(['X_' num2str(ifile+1) '=Coord_file(:,1:3)'';']);
555                end
556            end
557        end
558    end
559end
560n_ima=numel(coord_files)+1;
561est_dist=[0;0;0;0;0];
562est_aspect_ratio=0;
563est_fc=[1;1];
564%fc=[25;25]/0.012;
565center_optim=0;
566run(fullfile(path_UVMAT,'toolbox_calib','go_calib_optim'));
567GeometryCalib.CalibrationType='3D_linear';
568GeometryCalib.fx_fy=fc';
569%GeometryCalib.focal=fc(2);
570%GeometryCalib.dpx_dpy=[1 1];
571GeometryCalib.Cx_Cy=cc';
572%GeometryCalib.sx=fc(1)/fc(2);
573GeometryCalib.kc=kc(1);
574%GeometryCalib.kappa1=-kc(1)/fc(2)^2;
575GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function
576GeometryCalib.Tx_Ty_Tz=Tc_1';
577GeometryCalib.R=Rc_1;
578GeometryCalib.R(2,1:3)=-GeometryCalib.R(2,1:3);%inversion of the y image coordinate
579GeometryCalib.Tx_Ty_Tz(2)=-GeometryCalib.Tx_Ty_Tz(2);%inversion of the y image coordinate
580GeometryCalib.Cx_Cy(2)=ny-GeometryCalib.Cx_Cy(2);%inversion of the y image coordinate
581GeometryCalib.omc=(180/pi)*omc_1;%angles in degrees
582GeometryCalib.ErrorRMS=[];
583GeometryCalib.ErrorMax=[];
584
585%------------------------------------------------------------------------
586function GeometryCalib=calib_3D_quadr(Coord,handles)
587%------------------------------------------------------------------
588
589path_uvmat=which('uvmat');% check the path detected for source file uvmat
590path_UVMAT=fileparts(path_uvmat); %path to UVMAT
591huvmat=findobj(allchild(0),'Tag','uvmat');
592hhuvmat=guidata(huvmat);
593% check_cond=0;
594coord_files=get(handles.coord_files,'String');
595if ischar(coord_files)
596    coord_files={coord_files};
597end
598if isempty(coord_files{1}) || isequal(coord_files,{''})
599    coord_files={};
600end
601
602%retrieve the calibration points stored in the files listed in the popup list coord_files
603x_1=Coord(:,4:5)';%px coordinates of the ref points
604nx=str2num(get(hhuvmat.num_Npx,'String'));
605ny=str2num(get(hhuvmat.num_Npy,'String'));
606x_1(2,:)=ny-x_1(2,:);%reverse the y image coordinates
607X_1=Coord(:,1:3)';%phys coordinates of the ref points
608n_ima=numel(coord_files)+1;
609if ~isempty(coord_files)
610    msgbox_uvmat('CONFIRMATION',['The xy coordinates of the calibration points in ' num2str(n_ima) ' planes will be used'])
611    for ifile=1:numel(coord_files)
612    t=xmltree(coord_files{ifile});
613    s=convert(t);%convert to matlab structure
614        if isfield(s,'GeometryCalib')
615            if isfield(s.GeometryCalib,'SourceCalib')
616                if isfield(s.GeometryCalib.SourceCalib,'PointCoord')
617                PointCoord=s.GeometryCalib.SourceCalib.PointCoord;
618                Coord_file=zeros(length(PointCoord),5);%default
619                for i=1:length(PointCoord)
620                    line=str2num(PointCoord{i});
621                    Coord_file(i,4:5)=line(4:5);%px x
622                    Coord_file(i,1:3)=line(1:3);%phys x
623                end
624                eval(['x_' num2str(ifile+1) '=Coord_file(:,4:5)'';']);
625                eval(['x_' num2str(ifile+1) '(2,:)=ny-x_' num2str(ifile+1) '(2,:);' ]);
626                eval(['X_' num2str(ifile+1) '=Coord_file(:,1:3)'';']);
627                end
628            end
629        end
630    end
631end
632n_ima=numel(coord_files)+1;
633est_dist=[1;0;0;0;0];
634est_aspect_ratio=1;
635%est_fc=[0;0];
636%fc=[25;25]/0.012;
637center_optim=0;
638run(fullfile(path_UVMAT,'toolbox_calib','go_calib_optim'));
639
640GeometryCalib.CalibrationType='3D_quadr';
641GeometryCalib.fx_fy=fc';
642%GeometryCalib.focal=fc(2);
643%GeometryCalib.dpx_dpy=[1 1];
644GeometryCalib.Cx_Cy=cc';
645%GeometryCalib.sx=fc(1)/fc(2);
646GeometryCalib.kc=kc(1);
647%GeometryCalib.kappa1=-kc(1)/fc(2)^2;
648GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function
649GeometryCalib.Tx_Ty_Tz=Tc_1';
650GeometryCalib.R=Rc_1;
651GeometryCalib.R(2,1:3)=-GeometryCalib.R(2,1:3);%inversion of the y image coordinate
652GeometryCalib.Tx_Ty_Tz(2)=-GeometryCalib.Tx_Ty_Tz(2);%inversion of the y image coordinate
653GeometryCalib.Cx_Cy(2)=ny-GeometryCalib.Cx_Cy(2);%inversion of the y image coordinate
654GeometryCalib.omc=(180/pi)*omc_1;%angles in degrees
655GeometryCalib.ErrorRMS=[];
656GeometryCalib.ErrorMax=[];
657
658
659%------------------------------------------------------------------------
660function GeometryCalib=calib_3D_extrinsic(Coord,handles)
661%------------------------------------------------------------------
662path_uvmat=which('geometry_calib');% check the path detected for source file uvmat
663path_UVMAT=fileparts(path_uvmat); %path to UVMAT
664x_1=double(Coord(:,4:5)');%image coordiantes
665X_1=double(Coord(:,1:3)');% phys coordinates
666huvmat=findobj(allchild(0),'Tag','uvmat');
667hhuvmat=guidata(huvmat);
668ny=str2double(get(hhuvmat.num_Npy,'String'));
669x_1(2,:)=ny-x_1(2,:);%reverse the y image coordinates
670n_ima=1;
671GeometryCalib.CalibrationType='3D_extrinsic';
672GeometryCalib.fx_fy(1)=str2num(get(handles.fx,'String'));
673GeometryCalib.fx_fy(2)=str2num(get(handles.fy,'String'));
674GeometryCalib.Cx_Cy(1)=str2num(get(handles.Cx,'String'));
675GeometryCalib.Cx_Cy(2)=str2num(get(handles.Cy,'String'));
676GeometryCalib.kc=str2num(get(handles.kc,'String'));
677fct_path=fullfile(path_UVMAT,'toolbox_calib');
678addpath(fct_path)
679GeometryCalib.Cx_Cy(2)=ny-GeometryCalib.Cx_Cy(2);%reverse Cx_Cy(2) for calibration (inversion of px ordinate)
680% [omc1,Tc1,Rc1,H,x,ex,JJ] = compute_extrinsic(x_1,X_1,...
681%     [Calib.f Calib.f*Calib.sx]',...
682%     [Calib.Cx Calib.Cy]',...
683%     [-Calib.kappa1*Calib.f^2 0 0 0 0]);
684[omc,Tc1,Rc1,H,x,ex,JJ] = compute_extrinsic(x_1,X_1,...
685    (GeometryCalib.fx_fy)',GeometryCalib.Cx_Cy',[GeometryCalib.kc 0 0 0 0]);
686rmpath(fct_path);
687GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function
688GeometryCalib.Tx_Ty_Tz=Tc1';
689%inversion of z axis
690GeometryCalib.R=Rc1;
691GeometryCalib.R(2,1:3)=-GeometryCalib.R(2,1:3);%inversion of the y image coordinate
692GeometryCalib.Tx_Ty_Tz(2)=-GeometryCalib.Tx_Ty_Tz(2);%inversion of the y image coordinate
693GeometryCalib.Cx_Cy(2)=ny-GeometryCalib.Cx_Cy(2);%inversion of the y image coordinate
694GeometryCalib.omc=(180/pi)*omc';
695%GeometryCalib.R(3,1:3)=-GeometryCalib.R(3,1:3);%inversion for z upward
696
697
698
699%------------------------------------------------------------------------
700%function GeometryCalib=calib_tsai_heikkila(Coord)
701% TEST: NOT IMPLEMENTED
702%------------------------------------------------------------------
703% path_uvmat=which('uvmat');% check the path detected for source file uvmat
704% path_UVMAT=fileparts(path_uvmat); %path to UVMAT
705% path_calib=fullfile(path_UVMAT,'toolbox_calib_heikkila');
706% addpath(path_calib)
707% npoints=size(Coord,1);
708% Coord(:,1:3)=10*Coord(:,1:3);
709% Coord=[Coord zeros(npoints,2) -ones(npoints,1)];
710% [par,pos,iter,res,er,C]=cacal('dalsa',Coord);
711% GeometryCalib.CalibrationType='tsai';
712% GeometryCalib.focal=par(2);
713
714
715%------------------------------------------------------------------------
716% --- determine the rms of calibration error
717function ErrorRms=error_calib(calib_param,Calib,Coord)
718%calib_param: vector of free calibration parameters (to optimise)
719%Calib: structure of the given calibration parameters
720%Coord: list of phys coordinates (columns 1-3, and pixel coordinates (columns 4-5)
721Calib.f=25;
722Calib.dpx=0.012;
723Calib.dpy=0.012;
724Calib.sx=1;
725Calib.Cx=512;
726Calib.Cy=512;
727Calib.kappa1=calib_param(1);
728Calib.Tx=calib_param(2);
729Calib.Ty=calib_param(3);
730Calib.Tz=calib_param(4);
731alpha=calib_param(5);
732Calib.R=[cos(alpha) sin(alpha) 0;-sin(alpha) cos(alpha) 0;0 0 -1];
733
734X=Coord(:,1);
735Y=Coord(:,2);
736Z=Coord(:,3);
737x_ima=Coord(:,4);
738y_ima=Coord(:,5);
739[Xpoints,Ypoints]=px_XYZ(Calib,X,Y,Z);
740ErrorRms(1)=sqrt(mean((Xpoints-x_ima).*(Xpoints-x_ima)));
741ErrorRms(2)=sqrt(mean((Ypoints-y_ima).*(Ypoints-y_ima)));
742ErrorRms=mean(ErrorRms);
743
744%------------------------------------------------------------------------
745function XImage_Callback(hObject, eventdata, handles)
746%------------------------------------------------------------------------
747update_list(hObject, eventdata,handles)
748
749%------------------------------------------------------------------------
750function YImage_Callback(hObject, eventdata, handles)
751%------------------------------------------------------------------------
752update_list(hObject, eventdata,handles)
753
754%------------------------------------------------------------------------
755% --- Executes on button press in STORE.
756function STORE_Callback(hObject, eventdata, handles)
757Coord_cell=get(handles.ListCoord,'String');
758Object=read_geometry_calib(Coord_cell);
759unitlist=get(handles.CoordUnit,'String');
760unit=unitlist{get(handles.CoordUnit,'value')};
761GeometryCalib.CoordUnit=unit;
762GeometryCalib.SourceCalib.PointCoord=Object.Coord;
763huvmat=findobj(allchild(0),'Name','uvmat');
764hhuvmat=guidata(huvmat);%handles of elements in the GUI uvmat
765% RootPath='';
766% RootFile='';
767if ~isempty(hhuvmat.RootPath)&& ~isempty(hhuvmat.RootFile)
768%     testhandle=1;
769    RootPath=get(hhuvmat.RootPath,'String');
770    RootFile=get(hhuvmat.RootFile,'String');
771    filebase=fullfile(RootPath,RootFile);
772    while exist([filebase '.xml'],'file')
773        filebase=[filebase '~'];
774    end
775    outputfile=[filebase '.xml'];
776    errormsg=update_imadoc(GeometryCalib,outputfile);
777    if ~strcmp(errormsg,'')
778        msgbox_uvmat('ERROR',errormsg);
779    end
780    listfile=get(handles.coord_files,'string');
781    if isequal(listfile,{''})
782        listfile={outputfile};
783    else
784        listfile=[listfile;{outputfile}];%update the list of coord files
785    end
786    set(handles.coord_files,'string',listfile);
787end
788set(handles.ListCoord,'Value',1)% refresh the display of coordinates
789set(handles.ListCoord,'String',{'......'})
790
791% --------------------------------------------------------------------
792% --- Executes on button press in CLEAR_PTS: clear the list of calibration points
793function CLEAR_PTS_Callback(hObject, eventdata, handles)
794% --------------------------------------------------------------------
795set(handles.ListCoord,'Value',1)% refresh the display of coordinates
796set(handles.ListCoord,'String',{'......'})
797MenuPlot_Callback(hObject, eventdata, handles)
798
799%------------------------------------------------------------------------
800% --- Executes on button press in CLEAR.
801function CLEAR_Callback(hObject, eventdata, handles)
802%------------------------------------------------------------------------
803set(handles.coord_files,'Value',1)
804set(handles.coord_files,'String',{''})
805
806%------------------------------------------------------------------------
807function XObject_Callback(hObject, eventdata, handles)
808%------------------------------------------------------------------------
809update_list(hObject, eventdata,handles)
810
811%------------------------------------------------------------------------
812function YObject_Callback(hObject, eventdata, handles)
813%------------------------------------------------------------------------
814update_list(hObject, eventdata,handles)
815
816%------------------------------------------------------------------------
817function ZObject_Callback(hObject, eventdata, handles)
818%------------------------------------------------------------------------
819update_list(hObject, eventdata,handles)
820
821%------------------------------------------------------------------------
822function update_list(hObject, eventdata, handles)
823%------------------------------------------------------------------------
824newval(4)=str2double(get(handles.XImage,'String'));
825newval(5)=str2double(get(handles.YImage,'String'));
826newval(1)=str2double(get(handles.XObject,'String'));
827newval(2)=str2double(get(handles.YObject,'String'));
828newval(3)=str2double(get(handles.ZObject,'String'));
829if isnan(newval(3))
830    newval(3)=0;%put z to 0 by default
831end
832Coord=get(handles.ListCoord,'String');
833Coord(end)=[]; %remove last string '.....'
834val=get(handles.ListCoord,'Value');
835data=read_geometry_calib(Coord);
836data.Coord(val,:)=newval;
837for i=1:size(data.Coord,1)
838    for j=1:5
839          Coord_cell{i,j}=num2str(data.Coord(i,j),4);%display coordiantes with 4 digits
840    end
841end
842
843Tabchar=cell2tab(Coord_cell,' | ');
844Tabchar=[Tabchar ;{'......'}];
845set(handles.ListCoord,'String',Tabchar)
846
847%update the plot
848ListCoord_Callback(hObject, eventdata, handles)
849MenuPlot_Callback(hObject, eventdata, handles)
850
851%------------------------------------------------------------------------
852% --- Executes on selection change in ListCoord.
853function ListCoord_Callback(hObject, eventdata, handles)
854%------------------------------------------------------------------------
855huvmat=findobj(allchild(0),'Name','uvmat');%find the current uvmat interface handle
856hplot=findobj(huvmat,'Tag','axes3');%main plotting axis of uvmat
857hhh=findobj(hplot,'Tag','calib_marker');
858Coord_cell=get(handles.ListCoord,'String');
859val=get(handles.ListCoord,'Value');
860if numel(val)>1
861    return %no action if several lines have been selected
862end
863coord_str=Coord_cell{val};
864k=findstr(' | ',coord_str);
865if isempty(k)%last line '.....' selected
866    if ~isempty(hhh)
867        delete(hhh)%delete the circle marker
868    end
869    return
870end
871%fill the edit boxex
872set(handles.XObject,'String',coord_str(1:k(1)-1))
873set(handles.YObject,'String',coord_str(k(1)+3:k(2)-1))
874set(handles.ZObject,'String',coord_str(k(2)+3:k(3)-1))
875set(handles.XImage,'String',coord_str(k(3)+3:k(4)-1))
876set(handles.YImage,'String',coord_str(k(4)+3:end))
877h_menu_coord=findobj(huvmat,'Tag','transform_fct');
878menu=get(h_menu_coord,'String');
879choice=get(h_menu_coord,'Value');
880if iscell(menu)
881    option=menu{choice};
882else
883    option='px'; %default
884end
885if isequal(option,'phys')
886    XCoord=str2double(coord_str(1:k(1)-1));
887    YCoord=str2double(coord_str(k(1)+3:k(2)-1));
888elseif isequal(option,'px')|| isequal(option,'')
889    XCoord=str2double(coord_str(k(3)+3:k(4)-1));
890    YCoord=str2double(coord_str(k(4)+3:end));
891else
892    msgbox_uvmat('ERROR','the choice in menu_coord of uvmat must be px or phys ')
893end
894if isempty(XCoord)||isempty(YCoord)
895     if ~isempty(hhh)
896        delete(hhh)%delete the circle marker
897    end
898    return
899end
900xlim=get(hplot,'XLim');
901ylim=get(hplot,'YLim');
902ind_range=max(abs(xlim(2)-xlim(1)),abs(ylim(end)-ylim(1)))/20;%defines the size of the circle marker
903if isempty(hhh)
904    set(0,'CurrentFig',huvmat)
905    set(huvmat,'CurrentAxes',hplot)
906    rectangle('Curvature',[1 1],...
907              'Position',[XCoord-ind_range/2 YCoord-ind_range/2 ind_range ind_range],'EdgeColor','m',...
908              'LineStyle','-','Tag','calib_marker');
909else
910    set(hhh,'Position',[XCoord-ind_range/2 YCoord-ind_range/2 ind_range ind_range])
911end
912
913%------------------------------------------------------------------------
914% --- Executes on selection change in edit_append.
915function edit_append_Callback(hObject, eventdata, handles)
916%------------------------------------------------------------------------
917choice=get(handles.edit_append,'Value');
918if choice
919    set(handles.edit_append,'BackgroundColor',[1 1 0])
920    huvmat=findobj(allchild(0),'tag','uvmat');
921    if ishandle(huvmat)
922        hhuvmat=guidata(huvmat);
923        set(hhuvmat.edit_object,'Value',0)
924        set(hhuvmat.edit_object,'BackgroundColor',[0.7 0.7 0.7])
925    end
926else
927    set(handles.edit_append,'BackgroundColor',[0.7 0.7 0.7])
928end
929   
930function NEW_Callback(hObject, eventdata, handles)
931%A METTRE SOUS UN BOUTON
932huvmat=findobj(allchild(0),'Name','uvmat');
933hchild=get(huvmat,'children');
934hcoord=findobj(hchild,'Tag','menu_coord');
935coordtype=get(hcoord,'Value');
936haxes=findobj(hchild,'Tag','axes3');
937AxeData=get(haxes,'UserData');
938if ~isequal(hcoord,2)
939    set(hcoord,'Value',2)
940    huvmat=uvmat(AxeData);
941    'relancer uvmat';
942end
943if ~isfield(AxeData,'ZoomAxes')
944    msgbox_uvmat('ERROR','first draw a window around a grid marker')
945    return
946end
947XLim=get(AxeData.ZoomAxes,'XLim');
948YLim=get(AxeData.ZoomAxes,'YLim');
949np=size(AxeData.A);
950ind_sub_x=round(XLim);
951ind_sub_y=np(1)-round(YLim);
952Mfiltre=AxeData.A(ind_sub_y(2):ind_sub_y(1) ,ind_sub_x,:);
953Mfiltre_norm=double(Mfiltre);
954Mfiltre_norm=Mfiltre_norm/sum(sum(Mfiltre_norm));
955Mfiltre_norm=100*(Mfiltre_norm-mean(mean(Mfiltre_norm)));
956Atype=class(AxeData.A);
957Data.NbDim=2;
958Data.A=filter2(Mfiltre_norm,double(AxeData.A));
959Data.A=feval(Atype,Data.A);
960Data.AName='image';
961Data.AX=AxeData.AX;
962Data.AY=AxeData.AY;
963Data.CoordType='px';
964plot_field(Data)
965
966%------------------------------------------------------------------------
967function MenuPlot_Callback(hObject, eventdata, handles)
968%------------------------------------------------------------------------
969huvmat=findobj(allchild(0),'Name','uvmat');%find the current uvmat interface handle
970%UvData=get(huvmat,'UserData');%Data associated to the current uvmat interface
971hhuvmat=guidata(huvmat); %handles of GUI elements in uvmat
972%hplot=findobj(huvmat,'Tag','axes3');%main plotting axis of uvmat
973h_menu_coord=findobj(huvmat,'Tag','transform_fct');
974menu=get(h_menu_coord,'String');
975choice=get(h_menu_coord,'Value');
976if iscell(menu)
977    option=menu{choice};
978else
979    option='px'; %default
980end
981Coord_cell=get(handles.ListCoord,'String');
982ObjectData=read_geometry_calib(Coord_cell);
983%ObjectData=read_geometry_calib(handles);%read the interface input parameters defining the object
984if ~isempty(ObjectData.Coord)
985    if isequal(option,'phys')
986        ObjectData.Coord=ObjectData.Coord(:,1:3);
987    elseif isequal(option,'px')||isequal(option,'')
988        ObjectData.Coord=ObjectData.Coord(:,4:5);
989    else
990        msgbox_uvmat('ERROR','the choice in menu_coord of uvmat must be '''', px or phys ')
991    end
992end
993%axes(hhuvmat.axes3)
994set(0,'CurrentFigure',huvmat)
995set(huvmat,'CurrentAxes',hhuvmat.axes3)
996hh=findobj('Tag','calib_points');
997if  ~isempty(ObjectData.Coord) && isempty(hh)
998    hh=line(ObjectData.Coord(:,1),ObjectData.Coord(:,2),'Color','m','Tag','calib_points','LineStyle','.','Marker','+');
999elseif isempty(ObjectData.Coord)%empty list of points, suppress the plot
1000    delete(hh)
1001else
1002    set(hh,'XData',ObjectData.Coord(:,1))
1003    set(hh,'YData',ObjectData.Coord(:,2))
1004end
1005pause(.1)
1006figure(handles.geometry_calib)
1007
1008% --------------------------------------------------------------------
1009function MenuHelp_Callback(hObject, eventdata, handles)
1010path_to_uvmat=which('uvmat');% check the path of uvmat
1011pathelp=fileparts(path_to_uvmat);
1012helpfile=fullfile(pathelp,'uvmat_doc','uvmat_doc.html');
1013if isempty(dir(helpfile)), msgbox_uvmat('ERROR','Please put the help file uvmat_doc.html in the sub-directory /uvmat_doc of the UVMAT package')
1014else
1015   addpath (fullfile(pathelp,'uvmat_doc'))
1016   web([helpfile '#geometry_calib'])
1017end
1018
1019%------------------------------------------------------------------------
1020function MenuCreateGrid_Callback(hObject, eventdata, handles)
1021%------------------------------------------------------------------------
1022%hcalib=get(handles.calib_type,'parent');%handles of the GUI geometry_calib
1023CalibData=get(handles.geometry_calib,'UserData');
1024Tinput=[];%default
1025if isfield(CalibData,'grid')
1026    Tinput=CalibData.grid;
1027end
1028[T,CalibData.grid]=create_grid(Tinput);%display the GUI create_grid
1029set(handles.geometry_calib,'UserData',CalibData)
1030
1031%grid in phys space
1032Coord=get(handles.ListCoord,'String');
1033val=get(handles.ListCoord,'Value');
1034data=read_geometry_calib(Coord);
1035%nbpoints=size(data.Coord,1); %nbre of calibration points
1036data.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
1037% for i=1:nbpoints
1038%    for j=1:5
1039%           Coord{i,j}=num2str(data.Coord(i,j),4);%display coordiantes with 4 digits
1040%    end
1041% end
1042%update the phys coordinates starting from the selected point (down in the
1043Coord(end,:)=[]; %remove last string '.....'
1044for i=1:size(data.Coord,1)
1045    for j=1:5
1046          Coord{i,j}=num2str(data.Coord(i,j),4);%display coordiantes with 4 digits
1047    end
1048end
1049
1050%size(data.Coord,1)
1051Tabchar=cell2tab(Coord,' | ');
1052Tabchar=[Tabchar ;{'......'}];
1053set(handles.ListCoord,'String',Tabchar)
1054
1055% -----------------------------------------------------------------------
1056% --- automatic grid dectection from local maxima of the images
1057function MenuDetectGrid_Callback(hObject, eventdata, handles)
1058%------------------------------------------------------------------------
1059%% read the four last point coordinates in pixels
1060Coord_cell=get(handles.ListCoord,'String');%read list of coordinates on geometry_calib
1061data=read_geometry_calib(Coord_cell);
1062nbpoints=size(data.Coord,1); %nbre of calibration points
1063if nbpoints~=4
1064    msgbox_uvmat('ERROR','four points must have be selected by the mouse to delimitate the phys grid area; the Ox axis will be defined by the two first points')
1065    return
1066end
1067corners_X=(data.Coord(end:-1:end-3,4)); %pixel absissa of the four corners
1068corners_Y=(data.Coord(end:-1:end-3,5));
1069
1070%reorder the last two points (the two first in the list) if needed
1071angles=angle((corners_X-corners_X(1))+1i*(corners_Y-corners_Y(1)));
1072if abs(angles(4)-angles(2))>abs(angles(3)-angles(2))
1073      X_end=corners_X(4);
1074      Y_end=corners_Y(4);
1075      corners_X(4)=corners_X(3);
1076      corners_Y(4)=corners_Y(3);
1077      corners_X(3)=X_end;
1078      corners_Y(3)=Y_end;
1079end
1080
1081%% initiate the grid
1082CalibData=get(handles.geometry_calib,'UserData');%get information stored on the GUI geometry_calib
1083grid_input=[];%default
1084if isfield(CalibData,'grid')
1085    grid_input=CalibData.grid;%retrieve the previously used grid
1086end
1087[T,CalibData.grid,white_test]=create_grid(grid_input,'detect_grid');%display the GUI create_grid, read the set of phys coordinates T
1088set(handles.geometry_calib,'UserData',CalibData)%store the phys grid parameters for later use
1089
1090
1091
1092%% read the current image, displayed in the GUI uvmat
1093huvmat=findobj(allchild(0),'Name','uvmat');
1094UvData=get(huvmat,'UserData');
1095A=UvData.Field.A;
1096npxy=size(A);
1097X=[CalibData.grid.x_0 CalibData.grid.x_1 CalibData.grid.x_0 CalibData.grid.x_1]';%corner absissa in the phys coordinates (cm)
1098Y=[CalibData.grid.y_0 CalibData.grid.y_0 CalibData.grid.y_1 CalibData.grid.y_1]';%corner ordinates in the phys coordinates (cm)
1099
1100%calculate transform matrices for plane projection
1101% reference: http://alumni.media.mit.edu/~cwren/interpolator/ by Christopher R. Wren
1102B = [ X Y ones(size(X)) zeros(4,3)        -X.*corners_X -Y.*corners_X ...
1103      zeros(4,3)        X Y ones(size(X)) -X.*corners_Y -Y.*corners_Y ];
1104B = reshape (B', 8 , 8 )';
1105D = [ corners_X , corners_Y ];
1106D = reshape (D', 8 , 1 );
1107l = (B' * B)\B' * D;
1108Amat = reshape([l(1:6)' 0 0 1 ],3,3)';
1109C = [l(7:8)' 1];
1110
1111% transform grid image into 'phys' coordinates
1112GeometryCalib.fx_fy=[1 1];
1113GeometryCalib.Tx_Ty_Tz=[Amat(1,3) Amat(2,3) 1];
1114GeometryCalib.R=[Amat(1,1),Amat(1,2),0;Amat(2,1),Amat(2,2),0;C(1),C(2),0];
1115GeometryCalib.CoordUnit='cm';
1116path_uvmat=which('uvmat');% check the path detected for source file uvmat
1117path_UVMAT=fileparts(path_uvmat); %path to UVMAT
1118addpath(fullfile(path_UVMAT,'transform_field'))
1119Data.ListVarName={'AY','AX','A'};
1120Data.VarDimName={'AY','AX',{'AY','AX'}};
1121if ndims(A)==3
1122    A=mean(A,3);
1123end
1124Data.A=A-min(min(A));
1125Data.AY=[npxy(1)-0.5 0.5];
1126Data.AX=[0.5 npxy(2)];
1127Data.CoordUnit='pixel';
1128Calib.GeometryCalib=GeometryCalib;
1129DataOut=phys(Data,Calib);
1130rmpath(fullfile(path_UVMAT,'transform_field'))
1131Amod=DataOut.A;
1132Rangx=DataOut.AX;
1133Rangy=DataOut.AY;
1134if white_test
1135    Amod=double(Amod);%case of white grid markers: will look for image maxima
1136else
1137    Amod=-double(Amod);%case of black grid markers: will look for image minima
1138end
1139% figure(12) %display corrected image
1140% Amax=max(max(Amod));
1141% image(Rangx,Rangy,uint8(255*Amod/Amax))
1142
1143%% detection of local image extrema in each direction
1144Dx=(Rangx(2)-Rangx(1))/(npxy(2)-1); %x mesh in real space
1145Dy=(Rangy(2)-Rangy(1))/(npxy(1)-1); %y mesh in real space
1146ind_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
1147ind_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
1148nbpoints=size(T,1);
1149for ipoint=1:nbpoints
1150    i0=1+round((T(ipoint,1)-Rangx(1))/Dx);%round(Xpx(ipoint));
1151    j0=1+round((T(ipoint,2)-Rangy(1))/Dy);%round(Xpx(ipoint));
1152    j0min=max(j0-ind_range_y,1);
1153    j0max=min(j0+ind_range_y,size(Amod,1));
1154    i0min=max(i0-ind_range_x,1);
1155    i0max=min(i0+ind_range_x,size(Amod,2));
1156    Asub=Amod(j0min:j0max,i0min:i0max);
1157    x_profile=sum(Asub,1);
1158    y_profile=sum(Asub,2);
1159    [Amax,ind_x_max]=max(x_profile);
1160    [Amax,ind_y_max]=max(y_profile);
1161    %sub-pixel improvement using moments
1162    x_shift=0;
1163    y_shift=0;
1164    if ind_x_max+2<=numel(x_profile) && ind_x_max-2>=1
1165        Atop=x_profile(ind_x_max-2:ind_x_max+2);
1166        x_shift=sum(Atop.*[-2 -1 0 1 2])/sum(Atop);
1167    end
1168    if ind_y_max+2<=numel(y_profile) && ind_y_max-2>=1
1169        Atop=y_profile(ind_y_max-2:ind_y_max+2);
1170        y_shift=sum(Atop.*[-2 -1 0 1 2]')/sum(Atop);
1171    end
1172    Delta(ipoint,1)=(x_shift+ind_x_max+i0min-i0-1)*Dx;%shift from the initial guess
1173    Delta(ipoint,2)=(y_shift+ind_y_max+j0min-j0-1)*Dy;
1174end
1175Tmod=T(:,(1:2))+Delta;
1176[Xpx,Ypx]=px_XYZ(GeometryCalib,Tmod(:,1),Tmod(:,2));
1177for ipoint=1:nbpoints
1178     Coord{ipoint,1}=num2str(T(ipoint,1),4);%display coordiantes with 4 digits
1179     Coord{ipoint,2}=num2str(T(ipoint,2),4);%display coordiantes with 4 digits
1180     Coord{ipoint,3}=num2str(T(ipoint,3),4);%display coordiantes with 4 digits;
1181     Coord{ipoint,4}=num2str(Xpx(ipoint),4);%display coordiantes with 4 digits
1182     Coord{ipoint,5}=num2str(Ypx(ipoint),4);%display coordiantes with 4 digits
1183end
1184Tabchar=cell2tab(Coord(end:-1:1,:),' | ');
1185Tabchar=[Tabchar ;{'......'}];
1186set(handles.ListCoord,'Value',1)
1187set(handles.ListCoord,'String',Tabchar)
1188MenuPlot_Callback(hObject, eventdata, handles)
1189
1190%-----------------------------------------------------------------------
1191function MenuTranslatePoints_Callback(hObject, eventdata, handles)
1192%-----------------------------------------------------------------------
1193%hcalib=get(handles.calib_type,'parent');%handles of the GUI geometry_calib
1194CalibData=get(handles.geometry_calib,'UserData');
1195Tinput=[];%default
1196if isfield(CalibData,'translate')
1197    Tinput=CalibData.translate;
1198end
1199T=translate_points(Tinput);%display translate_points GUI and get shift parameters
1200CalibData.translate=T;
1201set(handles.geometry_calib,'UserData',CalibData)
1202%translation
1203Coord_cell=get(handles.ListCoord,'String');
1204data=read_geometry_calib(Coord_cell);
1205data.Coord(:,1)=T(1)+data.Coord(:,1);
1206data.Coord(:,2)=T(2)+data.Coord(:,2);
1207data.Coord(:,3)=T(3)+data.Coord(:,3);
1208data.Coord(:,[4 5])=data.Coord(:,[4 5]);
1209for i=1:size(data.Coord,1)
1210    for j=1:5
1211          Coord{i,j}=num2str(data.Coord(i,j),4);%phys x,y,z
1212   end
1213end
1214Tabchar=cell2tab(Coord,' | ');
1215Tabchar=[Tabchar; {'.....'}];
1216%set(handles.ListCoord,'Value',1)
1217set(handles.ListCoord,'String',Tabchar)
1218
1219
1220% --------------------------------------------------------------------
1221function MenuRotatePoints_Callback(hObject, eventdata, handles)
1222%hcalib=get(handles.calib_type,'parent');%handles of the GUI geometry_calib
1223CalibData=get(handles.geometry_calib,'UserData');
1224Tinput=[];%default
1225if isfield(CalibData,'rotate')
1226    Tinput=CalibData.rotate;
1227end
1228T=rotate_points(Tinput);%display rotate_points GUI to introduce rotation parameters
1229CalibData.rotate=T;
1230set(handles.geometry_calib,'UserData',CalibData)
1231%-----------------------------------------------------
1232%rotation
1233Phi=T(1);
1234O_x=0;%default
1235O_y=0;%default
1236if numel(T)>=2
1237    O_x=T(2);%default
1238end
1239if numel(T)>=3
1240    O_y=T(3);%default
1241end
1242Coord_cell=get(handles.ListCoord,'String');
1243data=read_geometry_calib(Coord_cell);
1244r1=cos(pi*Phi/180);
1245r2=-sin(pi*Phi/180);
1246r3=sin(pi*Phi/180);
1247r4=cos(pi*Phi/180);
1248x=data.Coord(:,1)-O_x;
1249y=data.Coord(:,2)-O_y;
1250data.Coord(:,1)=r1*x+r2*y;
1251data.Coord(:,2)=r3*x+r4*y;
1252% data.Coord(:,[4 5])=data.Coord(:,[4 5]);
1253for i=1:size(data.Coord,1)
1254    for j=1:5
1255          Coord{i,j}=num2str(data.Coord(i,j),4);%phys x,y,z
1256   end
1257end
1258Tabchar=cell2tab(Coord,' | ');
1259Tabchar=[Tabchar;{'......'}];
1260set(handles.ListCoord,'Value',1)
1261set(handles.ListCoord,'String',Tabchar)
1262
1263% --------------------------------------------------------------------
1264function MenuImportPoints_Callback(hObject, eventdata, handles)
1265fileinput=browse_xml(hObject, eventdata, handles);
1266if isempty(fileinput)
1267    return
1268end
1269[s,errormsg]=imadoc2struct(fileinput,'GeometryCalib');
1270GeometryCalib=s.GeometryCalib;
1271%GeometryCalib=load_calib(hObject, eventdata, handles)
1272calib=reshape(GeometryCalib.PointCoord,[],1);
1273for ilist=1:numel(calib)
1274    CoordCell{ilist}=num2str(calib(ilist));
1275end
1276CoordCell=reshape(CoordCell,[],5);
1277Tabchar=cell2tab(CoordCell,' | ');%transform cells into table ready for display
1278Tabchar=[Tabchar;{'......'}];
1279set(handles.ListCoord,'Value',1)
1280set(handles.ListCoord,'String',Tabchar)
1281MenuPlot_Callback(handles.geometry_calib, [], handles)
1282
1283% -----------------------------------------------------------------------
1284function MenuImportIntrinsic_Callback(hObject, eventdata, handles)
1285%------------------------------------------------------------------------
1286fileinput=browse_xml(hObject, eventdata, handles);
1287if isempty(fileinput)
1288    return
1289end
1290[s,errormsg]=imadoc2struct(fileinput,'GeometryCalib');
1291GeometryCalib=s.GeometryCalib;
1292display_intrinsic(GeometryCalib,handles)
1293
1294% -----------------------------------------------------------------------
1295function MenuImportAll_Callback(hObject, eventdata, handles)
1296%------------------------------------------------------------------------
1297fileinput=browse_xml(hObject, eventdata, handles);
1298if ~isempty(fileinput)
1299    loadfile(handles,fileinput)
1300end
1301
1302% -----------------------------------------------------------------------
1303% --- Executes on menubar option Import/Grid file: introduce previous grid files
1304function MenuGridFile_Callback(hObject, eventdata, handles)
1305% -----------------------------------------------------------------------
1306inputfile=browse_xml(hObject, eventdata, handles);
1307listfile=get(handles.coord_files,'string');
1308if isequal(listfile,{''})
1309    listfile={inputfile};
1310else
1311    listfile=[listfile;{inputfile}];%update the list of coord files
1312end
1313set(handles.coord_files,'string',listfile);
1314
1315%------------------------------------------------------------------------
1316% --- 'key_press_fcn:' function activated when a key is pressed on the keyboard
1317function key_press_fcn(hObject,eventdata,handles)
1318%------------------------------------------------------------------------
1319xx=double(get(handles.geometry_calib,'CurrentCharacter')); %get the keyboard character
1320if ismember(xx,[8 127])%backspace or delete
1321    Coord_cell=get(handles.ListCoord,'String');
1322    val=get(handles.ListCoord,'Value');
1323     if max(val)<numel(Coord_cell) % the last element '...' has not been selected
1324        Coord_cell(val)=[];%remove the selected line
1325        set(handles.ListCoord,'Value',min(val))
1326        set(handles.ListCoord,'String',Coord_cell)         
1327        ListCoord_Callback(hObject, eventdata, handles)
1328        MenuPlot_Callback(hObject,eventdata,handles)
1329     end
1330end
1331
1332%------------------------------------------------------------------------
1333function fileinput=browse_xml(hObject, eventdata, handles)
1334%------------------------------------------------------------------------
1335fileinput=[];%default
1336oldfile=''; %default
1337UserData=get(handles.geometry_calib,'UserData');
1338if isfield(UserData,'XmlInputFile')
1339    oldfile=UserData.XmlInputFile;
1340end
1341[FileName, PathName, filterindex] = uigetfile( ...
1342       {'*.xml;*.mat', ' (*.xml,*.mat)';
1343       '*.xml',  '.xml files '; ...
1344        '*.mat',  '.mat matlab files '}, ...
1345        'Pick a file',oldfile);
1346fileinput=[PathName FileName];%complete file name
1347testblank=findstr(fileinput,' ');%look for blanks
1348if ~isempty(testblank)
1349    msgbox_uvmat('ERROR','forbidden input file name or path: no blank character allowed')
1350    return
1351end
1352sizf=size(fileinput);
1353if (~ischar(fileinput)||~isequal(sizf(1),1)),return;end
1354UserData.XmlInputFile=fileinput;
1355set(handles.geometry_calib,'UserData',UserData)%record current file foer further use of browser
1356
1357% -----------------------------------------------------------------------
1358function Heading=loadfile(handles,fileinput)
1359%------------------------------------------------------------------------
1360Heading=[];%default
1361[s,errormsg]=imadoc2struct(fileinput,'GeometryCalib');
1362if ~isempty(errormsg)
1363    msgbox_uvmat('ERROR',['Error for reading ' fileinput ': '  errormsg])
1364    return
1365end
1366if ~isempty(s.Heading)
1367    Heading=s.Heading;
1368end
1369   
1370GeometryCalib=s.GeometryCalib;
1371fx=1;fy=1;Cx=0;Cy=0;kc=0; %default
1372%     Tabchar={};
1373CoordCell={};
1374%     kc=0;%default
1375%     f1=1000;
1376%     f2=1000;
1377%     hhuvmat=guidata(findobj(allchild(0),'Name','uvmat'));
1378%     Cx=str2num(get(hhuvmat.num_Npx,'String'))/2;
1379%     Cy=str2num(get(hhuvmat.num_Npy,'String'))/2;
1380Tabchar={};%default
1381val_cal=1;%default
1382if ~isempty(GeometryCalib)
1383    % choose the calibration option
1384    if isfield(GeometryCalib,'CalibrationType')
1385       calib_list=get(handles.calib_type,'String');
1386       for ilist=1:numel(calib_list)
1387           if strcmp(calib_list{ilist},GeometryCalib.CalibrationType)
1388               val_cal=ilist;
1389               break
1390           end
1391       end
1392    end
1393    display_intrinsic(GeometryCalib,handles)%intrinsic param
1394    %extrinsic param
1395    if isfield(GeometryCalib,'Tx_Ty_Tz')
1396        Tx_Ty_Tz=GeometryCalib.Tx_Ty_Tz;
1397        set(handles.Tx,'String',num2str(GeometryCalib.Tx_Ty_Tz(1),4))
1398        set(handles.Ty,'String',num2str(GeometryCalib.Tx_Ty_Tz(2),4))
1399        set(handles.Tz,'String',num2str(GeometryCalib.Tx_Ty_Tz(3),4))
1400    end
1401    if isfield(GeometryCalib,'omc')
1402        set(handles.Phi,'String',num2str(GeometryCalib.omc(1),4))
1403        set(handles.Theta,'String',num2str(GeometryCalib.omc(2),4))
1404        set(handles.Psi,'String',num2str(GeometryCalib.omc(3),4))
1405    end
1406    calib=reshape(GeometryCalib.PointCoord,[],1);
1407    for ilist=1:numel(calib)
1408        CoordCell{ilist}=num2str(calib(ilist));
1409    end
1410    CoordCell=reshape(CoordCell,[],5);
1411    Tabchar=cell2tab(CoordCell,' | ');%transform cells into table ready for display
1412    MenuPlot_Callback(handles.geometry_calib, [], handles)
1413end
1414set(handles.calib_type,'Value',val_cal)
1415Tabchar=[Tabchar;{'......'}];
1416set(handles.ListCoord,'Value',1)
1417set(handles.ListCoord,'String',Tabchar)
1418
1419if isempty(CoordCell)% allow mouse action by default in the absence of input points
1420    set(handles.edit_append,'Value',1)
1421    set(handles.edit_append,'BackgroundColor',[1 1 0])
1422else % does not allow mouse action by default in the presence of input points
1423    set(handles.edit_append,'Value',0)
1424    set(handles.edit_append,'BackgroundColor',[0.7 0.7 0.7])
1425end
1426
1427%------------------------------------------------------------------------
1428%---display calibration intrinsic parameters
1429function display_intrinsic(GeometryCalib,handles)
1430%------------------------------------------------------------------------
1431fx=[];
1432fy=[];
1433if isfield(GeometryCalib,'fx_fy')
1434    fx=GeometryCalib.fx_fy(1);
1435    fy=GeometryCalib.fx_fy(2);
1436end
1437Cx_Cy=[0 0];%default
1438if isfield(GeometryCalib,'Cx_Cy')
1439    Cx_Cy=GeometryCalib.Cx_Cy;
1440end
1441kc=0;
1442if isfield(GeometryCalib,'kc')
1443    kc=GeometryCalib.kc; %* GeometryCalib.focal*GeometryCalib.focal;
1444end
1445set(handles.fx,'String',num2str(fx,5))
1446set(handles.fy,'String',num2str(fy,5))
1447set(handles.Cx,'String',num2str(Cx_Cy(1),'%1.1f'))
1448set(handles.Cy,'String',num2str(Cx_Cy(2),'%1.1f'))
1449set(handles.kc,'String',num2str(kc,'%1.4f'))
1450
1451
Note: See TracBrowser for help on using the repository browser.