source: trunk/src/geometry_calib.m @ 679

Last change on this file since 679 was 674, checked in by sommeria, 11 years ago

various bugs repaired, in particula timing

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