source: trunk/src/geometry_calib.m @ 773

Last change on this file since 773 was 773, checked in by sommeria, 10 years ago

help web link corrected

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