source: trunk/src/geometry_calib.m @ 611

Last change on this file since 611 was 591, checked in by sommeria, 11 years ago

various updates, in particular modification of series to do calculations in the cluster

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