source: trunk/src/geometry_calib.m @ 356

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

civ updated with new functions for opening files, consistently with uvmat
Bugs to be expected (use previous version then)

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