source: trunk/src/geometry_calib.m @ 88

Last change on this file since 88 was 88, checked in by sommeria, 14 years ago

plot_field: test for multicomponent image (-> errormsg)
creat_grid: possibility of changing the figure name (for detect_grid)
geometry_calib: detect_grid improved (range of search adjusted from grid mesh)

bug repaired in translate points

uvmat: bug repair: dela with the case with several get_field GUI opened

File size: 44.5 KB
Line 
1%'geometry_calib': performs geometric calibration from a set of reference points
2%
3% function varargout = geometry_calib(varargin)
4%
5%A%AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
6%  Copyright Joel Sommeria, 2008, LEGI / CNRS-UJF-INPG, sommeria@coriolis-legi.org.
7%AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
8%     This file is part of the toolbox UVMAT.
9%
10%     UVMAT is free software; you can redistribute it and/or modify
11%     it under the terms of the GNU General Public License as published by
12%     the Free Software Foundation; either version 2 of the License, or
13%     (at your option) any later version.
14%
15%     UVMAT is distributed in the hope that it will be useful,
16%     but WITHOUT ANY WARRANTY; without even the implied warranty of
17%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
18%     GNU General Public License (file UVMAT/COPYING.txt) for more details.
19%AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
20
21function varargout = geometry_calib(varargin)
22% GEOMETRY_CALIB M-file for geometry_calib.fig
23%      GEOMETRY_CALIB, by itself, creates a MenuCoord GEOMETRY_CALIB or raises the existing
24%      singleton*.
25%
26%      H = GEOMETRY_CALIB returns the handle to a MenuCoord GEOMETRY_CALIB or the handle to
27%      the existing singleton*.
28%
29%      GEOMETRY_CALIB('CALLBACK',hObject,eventData,handles,...) calls the local
30%      function named CALLBACK in GEOMETRY_CALIB.M with the given input arguments.
31%
32%      GEOMETRY_CALIB('Property','Value',...) creates a MenuCoord GEOMETRY_CALIB or raises the
33%      existing singleton*.  Starting from the left, property value pairs are
34%      applied to the GUI before geometry_calib_OpeningFunction gets called.  An
35%      unrecognized property name or invalid value makes property application
36%      stop.  All inputs are passed to geometry_calib_OpeningFcn via varargin.
37%
38%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
39%      instance to run (singleton)".
40%
41% See also: GUIDE, GUIDATA, GUIHANDLES
42
43% Edit the above text to modify the response to help geometry_calib
44
45% Last Modified by GUIDE v2.5 25-Mar-2010 19:10:05
46
47% Begin initialization code - DO NOT edit
48gui_Singleton = 1;
49gui_State = struct('gui_Name',       mfilename, ...
50                   'gui_Singleton',  gui_Singleton, ...
51                   'gui_OpeningFcn', @geometry_calib_OpeningFcn, ...
52                   'gui_OutputFcn',  @geometry_calib_OutputFcn, ...
53                   'gui_LayoutFcn',  [] , ...
54                   'gui_Callback',   []);
55if nargin & isstr(varargin{1})
56    gui_State.gui_Callback = str2func(varargin{1});
57end
58
59if nargout
60    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
61else
62    gui_mainfcn(gui_State, varargin{:});
63end
64% End initialization code - DO NOT edit
65
66
67% --- Executes just before geometry_calib is made visible.
68%INPUT:
69%handles: handles of the geometry_calib interface elements
70% PlotHandles: set of handles of the elements contolling the plotting
71% parameters on the uvmat interface (obtained by 'get_plot_handle.m')
72%------------------------------------------------------------------------
73function geometry_calib_OpeningFcn(hObject, eventdata, handles, handles_uvmat,pos,inputfile)
74%------------------------------------------------------------------------
75% Choose default command line output for geometry_calib
76handles.output = hObject;
77
78% Update handles structure
79guidata(hObject, handles);
80set(hObject,'DeleteFcn',{@closefcn})%
81
82%set the position of the interface
83if exist('pos','var')& length(pos)>2
84    pos_gui=get(hObject,'Position');
85    pos_gui(1)=pos(1);
86    pos_gui(2)=pos(2);
87    set(hObject,'Position',pos_gui);
88end
89inputxml='';
90if exist('inputfile','var')& ~isempty(inputfile)
91%     [Path,Name,ext]=fileparts(inputfile);
92%     form=imformats(ext([2:end]));
93%     if ~isempty(form)% if the input file is an image
94        struct.XmlInputfile=inputfile;
95        set(hObject,'UserData',struct)
96        [Pathsub,RootFile,field_count,str2,str_a,str_b,ext,nom_type,subdir]=name2display(inputfile);
97        inputxml=[fullfile(Pathsub,RootFile) '.xml'];
98%     end   
99end
100set(handles.ListCoord,'String',{'......'})
101if exist(inputxml,'file')
102    loadfile(handles,inputxml)% load the point coordiantes existing in the xml file
103end
104
105set(handles.ListCoord,'KeyPressFcn',{@key_press_fcn,handles})%set keyboard action function
106
107
108%------------------------------------------------------------------------
109% --- Outputs from this function are returned to the command line.
110function varargout = geometry_calib_OutputFcn(hObject, eventdata, handles)
111%------------------------------------------------------------------------
112% Get default command line output from handles structure
113varargout{1} = handles.output;
114varargout{2}=handles;
115
116%------------
117function Phi_Callback(hObject, eventdata, handles)
118
119%------------------------------------------------------------------------
120%read input xml file and update the edit boxes
121function loadfile(handles,fileinput)
122%------------------------------------------------------------------------
123%read the input xml file
124t=xmltree(fileinput);
125s=convert(t);%convert to matlab structure
126%read data currently displayed on the interface
127PointCoord=[];
128Coord_cell=get(handles.ListCoord,'String');
129data=read_geometry_calib(Coord_cell);
130%data=read_geometry_calib(handles);
131Coord=[]; %default
132if isfield(data,'Coord')
133    Coord=data.Coord;
134end
135TabChar_0=get(handles.ListCoord,'String');
136nbcoord_0=size(TabChar_0,1);
137if isequal(get(handles.edit_append,'Value'),2) %edit mode  A REVOIR
138    val=get(handles.ListCoord,'Value')-1;
139else
140   val=length(TabChar_0);
141end
142nbcoord=0;
143
144%case of calibration (ImaDoc) input file
145% hcalib=get(handles.calib_type,'parent');
146CalibData=get(handles.geometry_calib,'UserData');
147CalibData.XmlInput=fileinput;
148if isfield(s,'Heading')
149    CalibData.Heading=s.Heading;
150end
151
152set(handles.geometry_calib,'UserData',CalibData);%store the heading in the interface 'UserData'
153if isfield(s,'GeometryCalib')
154    Calib=s.GeometryCalib;
155    if isfield(Calib,'CalibrationType')
156        CalibrationType=Calib.CalibrationType;
157        switch CalibrationType
158            case 'linear'
159                set(handles.calib_type,'Value',2)
160            case 'tsai'
161                set(handles.calib_type,'Value',3)
162            case 'tsai_matlab'
163                set(handles.calib_type,'Value',4)
164        end
165    end
166    if isfield(Calib,'SourceCalib')
167        if isfield(Calib.SourceCalib,'PointCoord')
168            PointCoord=Calib.SourceCalib.PointCoord;
169        end
170    end
171    nbcoord=length(PointCoord);
172    if ~isfield(Calib,'ErrorRms')&~isfield(Calib,'ErrorMax') %old convention of Gauthier (cord in mm)
173        for i=1:length(PointCoord)
174          line=str2num(PointCoord{i});
175          Coord(i+val,4:5)=line(4:5);%px x
176          Coord(i+val,1:3)=line(1:3)/10;%phys x
177        end
178    else
179        for i=1:length(PointCoord)
180          line=str2num(PointCoord{i});
181          Coord(i,4:5)=line(4:5);%px x
182          Coord(i,1:3)=line(1:3);%phys x
183       end
184    end
185end
186%case of xml files of points
187if isfield(s,'Coord')
188    PointCoord=s.Coord;
189    nbcoord=length(PointCoord);
190     %case of image coordinates
191    if isfield(s,'CoordType')& isequal(s.CoordType,'px')
192        for i=1:nbcoord
193           line=str2num(PointCoord{i});
194           Coord(i+val,4:5)=line(1:2);
195        end
196     %case of  physical coordinates
197    else
198        for i=1:nbcoord
199           line=str2num(PointCoord{i});
200           Coord(i+val,1:3)=line(1:3);
201           nbcolumn=size(Coord,2);
202           if nbcolumn<5
203               Coord(i+val,nbcolumn+1:5)=zeros(1,5-nbcolumn);
204           end
205        end
206     end
207end
208CoordCell={};
209for iline=1:size(Coord,1)
210    for j=1:5
211        CoordCell{iline,j}=num2str(Coord(iline,j),4);
212    end
213end
214% CoordCell=[CoordCell;{' ',' ',' ',' ',' '}];
215Tabchar=cell2tab(CoordCell,'    |    ');%transform cells into table ready for display
216Tabchar=[Tabchar;{'......'}];
217set(handles.ListCoord,'Value',1)
218set(handles.ListCoord,'String',Tabchar)
219MenuPlot_Callback(handles.geometry_calib, [], handles)
220if isempty(Coord)
221    set(handles.edit_append,'Value',1)
222    set(handles.edit_append,'BackgroundColor',[1 1 0])
223else
224    set(handles.edit_append,'Value',0)
225    set(handles.edit_append,'BackgroundColor',[0.7 0.7 0.7])
226end
227%
228%------------------------------------------------------------------------
229% executed when closing: set the parent interface button to value 0
230function closefcn(gcbo,eventdata)
231%------------------------------------------------------------------------
232huvmat=findobj(allchild(0),'Name','uvmat');
233if ~isempty(huvmat)
234    handles=guidata(huvmat);
235    set(handles.MenuMask,'enable','on')
236    set(handles.MenuGrid,'enable','on')
237    set(handles.MenuObject,'enable','on')
238    set(handles.MenuEdit,'enable','on')
239    set(handles.edit,'enable','on')
240    hobject=findobj(handles.axes3,'tag','calib_points');
241    if ~isempty(hobject)
242        delete(hobject)
243    end
244    hobject=findobj(handles.axes3,'tag','calib_marker');
245    if ~isempty(hobject)
246        delete(hobject)
247    end   
248end
249
250%------------------------------------------------------------------------
251% --- Executes on button press in calibrate_lin.
252function APPLY_Callback(hObject, eventdata, handles)
253%------------------------------------------------------------------------
254calib_cell=get(handles.calib_type,'String');
255val=get(handles.calib_type,'Value');
256calib_type=calib_cell{val};
257Coord_cell=get(handles.ListCoord,'String');
258Object=read_geometry_calib(Coord_cell);
259X=Object.Coord(:,1);
260Y=Object.Coord(:,2);
261Z=Object.Coord(:,3);
262if isequal(calib_type,'rescale')
263    GeometryCalib=calib_rescale(Object.Coord);
264    Z=0;%Z not taken into account
265elseif isequal(calib_type,'linear')
266    GeometryCalib=calib_linear(Object.Coord);
267    Z=0; %Z not taken into account
268elseif isequal(calib_type,'tsai_cpp')
269    GeometryCalib=calib_tsai(Object.Coord);
270elseif isequal(calib_type,'tsai_matlab')
271    GeometryCalib=calib_tsai2(Object.Coord);
272end
273
274%check error
275if isfield(GeometryCalib,'dpx_dpy')
276    Calib.dpx=GeometryCalib.dpx_dpy(1);
277    Calib.dpy=GeometryCalib.dpx_dpy(2);
278end
279if isfield(GeometryCalib,'sx')
280    Calib.sx=GeometryCalib.sx;
281end
282if isfield(GeometryCalib,'Cx_Cy')
283    Calib.Cx=GeometryCalib.Cx_Cy(1);
284    Calib.Cy=GeometryCalib.Cx_Cy(2);
285end
286if isfield(GeometryCalib,'kappa1')
287    Calib.kappa1=GeometryCalib.kappa1;
288end
289if isfield(GeometryCalib,'focal')
290    Calib.f=GeometryCalib.focal;
291end
292if isfield(GeometryCalib,'Tx_Ty_Tz')
293    Calib.Tx=GeometryCalib.Tx_Ty_Tz(1);
294    Calib.Ty=GeometryCalib.Tx_Ty_Tz(2);
295    Calib.Tz=GeometryCalib.Tx_Ty_Tz(3);
296end
297if isfield(GeometryCalib,'R')
298    Calib.R=GeometryCalib.R;
299end
300x_ima=Object.Coord(:,4);
301y_ima=Object.Coord(:,5);
302[Xpoints,Ypoints]=px_XYZ(Calib,X,Y,Z);
303GeometryCalib.ErrorRms(1)=sqrt(mean((Xpoints-x_ima).*(Xpoints-x_ima)));
304[GeometryCalib.ErrorMax(1),index(1)]=max(abs(Xpoints-x_ima));
305GeometryCalib.ErrorRms(2)=sqrt(mean((Ypoints-y_ima).*(Ypoints-y_ima)));
306[GeometryCalib.ErrorMax(2),index(2)]=max(abs(Ypoints-y_ima));
307[EM,ind_dim]=max(GeometryCalib.ErrorMax);
308index=index(ind_dim);
309
310unitlist=get(handles.CoordUnit,'String');
311unit=unitlist{get(handles.CoordUnit,'value')};
312GeometryCalib.CoordUnit=unit;
313GeometryCalib.SourceCalib.PointCoord=Object.Coord;
314huvmat=findobj(allchild(0),'Name','uvmat');
315hhuvmat=guidata(huvmat);%handles of elements in the GUI uvmat
316RootPath='';
317RootFile='';
318if ~isempty(hhuvmat.RootPath)& ~isempty(hhuvmat.RootFile)
319    testhandle=1;
320    RootPath=get(hhuvmat.RootPath,'String');
321    RootFile=get(hhuvmat.RootFile,'String');
322    filebase=fullfile(RootPath,RootFile);
323    outputfile=[filebase '.xml'];
324else
325    question={'save the calibration data and point coordinates in'};
326    def={fullfile(RootPath,['ObjectCalib.xml'])};
327    options.Resize='on';
328    answer=inputdlg(question,'save average in a new file',1,def,options);
329    outputfile=answer{1};
330end
331answer=msgbox_uvmat('INPUT_Y-N',{[outputfile ' updated with calibration data'];...
332    ['Error rms (along x,y)=' num2str(GeometryCalib.ErrorRms) ' pixels'];...
333    ['Error max (along x,y)=' num2str(GeometryCalib.ErrorMax) ' pixels']});
334if isequal(answer,'Yes')
335    update_imadoc(GeometryCalib,outputfile)
336    %display image with new calibration in the currently opened uvmat interface
337    hhh=findobj(hhuvmat.axes3,'Tag','calib_marker');% delete calib points and markers
338    if ~isempty(hhh)
339        delete(hhh);
340    end
341    hhh=findobj(hhuvmat.axes3,'Tag','calib_points');
342    if ~isempty(hhh)
343        delete(hhh);
344    end
345    set(hhuvmat.FixedLimits,'Value',0)% put FixedLimits option to 'off'
346    set(hhuvmat.FixedLimits,'BackgroundColor',[0.7 0.7 0.7])
347    uvmat('RootPath_Callback',hObject,eventdata,hhuvmat); %file input with xml reading  in uvmat
348    MenuPlot_Callback(hObject, eventdata, handles)
349    set(handles.ListCoord,'Value',index)% indicate in the list the point with max deviation (possible mistake)
350    ListCoord_Callback(hObject, eventdata, handles)
351    figure(handles.geometry_calib)
352end
353
354%------------------------------------------------------------------
355% --- Executes on button press in calibrate_lin.
356function REPLICATE_Callback(hObject, eventdata, handles)
357%------------------------------------------------------------------------
358calib_cell=get(handles.calib_type,'String');
359val=get(handles.calib_type,'Value');
360calib_type=calib_cell{val};
361Coord_cell=get(handles.ListCoord,'String');
362Object=read_geometry_calib(Coord_cell);
363
364if isequal(calib_type,'rescale')
365    GeometryCalib=calib_rescale(Object.Coord);
366elseif isequal(calib_type,'linear')
367    GeometryCalib=calib_linear(Object.Coord);
368elseif isequal(calib_type,'tsai')
369    GeometryCalib=calib_tsai(Object.Coord);
370end
371% %record image source
372GeometryCalib.SourceCalib.PointCoord=Object.Coord;
373
374%open and read the dataview GUI
375h_dataview=findobj(allchild(0),'name','dataview');
376if ~isempty(h_dataview)
377    delete(h_dataview)
378end
379CalibData=get(handles.geometry_calib,'UserData');%read the calibration image source on the interface userdata
380
381if isfield(CalibData,'XmlInput')
382    XmlInput=fileparts(CalibData.XmlInput);
383    [XmlInput,filename,ext]=fileparts(XmlInput);
384end
385SubCampaignTest='n'; %default
386testinput=0;
387if isfield(CalibData,'Heading')
388    Heading=CalibData.Heading;
389    if isfield(Heading,'Record') && isequal([filename ext],Heading.Record)
390        [XmlInput,filename,ext]=fileparts(XmlInput);
391    end
392    if isfield(Heading,'Device') && isequal([filename ext],Heading.Device)
393        [XmlInput,filename,ext]=fileparts(XmlInput);
394        Device=Heading.Device;
395    end
396    if isfield(Heading,'Experiment') && isequal([filename ext],Heading.Experiment)
397        [PP,filename,ext]=fileparts(XmlInput);
398    end
399    testinput=0;
400    if isfield(Heading,'SubCampaign') && isequal([filename ext],Heading.SubCampaign)
401        SubCampaignTest='y';
402        testinput=1;
403    elseif isfield(Heading,'Campaign') && isequal([filename ext],Heading.Campaign)
404        testinput=1;
405    end
406end
407if ~testinput
408    filename='PROJETS';%default
409    if isfield(CalibData,'XmlInput')
410         [pp,filename]=fileparts(CalibData.XmlInput);
411    end
412    while ~isequal(filename,'PROJETS') && numel(filename)>1
413        filename_1=filename;
414        pp_1=pp;
415        [pp,filename]=fileparts(pp);
416    end
417    XmlInput=fullfile(pp_1,filename_1);
418    testinput=1;
419end
420if testinput
421    outcome=dataview(XmlInput,SubCampaignTest,GeometryCalib);
422end
423
424%------------------------------------------------------------------------
425% determine the parameters for a calibration by an affine function (rescaling and offset, no rotation)
426function GeometryCalib=calib_rescale(Coord)
427%------------------------------------------------------------------------
428 
429X=Coord(:,1);
430Y=Coord(:,2);
431x_ima=Coord(:,4);
432y_ima=Coord(:,5);
433[px,sx]=polyfit(X,x_ima,1);
434[py,sy]=polyfit(Y,y_ima,1);
435T_x=px(2);
436T_y=py(2);
437GeometryCalib.CalibrationType='rescale';
438GeometryCalib.focal=1;
439GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function
440GeometryCalib.Tx_Ty_Tz=[T_x T_y 1];
441GeometryCalib.R=[px(1),0,0;0,py(1),0;0,0,1];
442
443% %check error
444% Calib.dpx=1;
445% Calib.dpy=1;
446% Calib.sx=1;
447% Calib.Cx=0;
448% Calib.Cy=0;
449% Calib.Tz=1;
450% Calib.kappa1=0;
451% Calib.f=GeometryCalib.focal;
452% Calib.Tx=T_x;
453% Calib.Ty=T_y;
454% Calib.R=GeometryCalib.R;
455% [Xpoints,Ypoints]=px_XYZ(Calib,X,Y,0);
456% GeometryCalib.ErrorRms(1)=sqrt(mean((Xpoints-x_ima).*(Xpoints-x_ima)));
457% GeometryCalib.ErrorMax(1)=max(abs(Xpoints-x_ima));
458% GeometryCalib.ErrorRms(2)=sqrt(mean((Ypoints-y_ima).*(Ypoints-y_ima)));
459% GeometryCalib.ErrorMax(2)=max(abs(Ypoints-y_ima));
460
461%------------------------------------------------------------------------
462% determine the parameters for a calibration by a linear transform matrix (rescale and rotation)
463function GeometryCalib=calib_linear(Coord)
464%------------------------------------------------------------------------
465X=Coord(:,1);
466Y=Coord(:,2);
467x_ima=Coord(:,4);
468y_ima=Coord(:,5);
469XY_mat=[ones(size(X)) X Y];
470a_X1=XY_mat\x_ima; %transformation matrix for X
471x1=XY_mat*a_X1;%reconstruction
472err_X1=max(abs(x1-x_ima));%error
473a_Y1=XY_mat\y_ima;%transformation matrix for X
474y1=XY_mat*a_Y1;
475err_Y1=max(abs(y1-y_ima));%error
476T_x=a_X1(1);
477T_y=a_Y1(1);
478GeometryCalib.CalibrationType='linear';
479GeometryCalib.focal=1;
480GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function
481GeometryCalib.Tx_Ty_Tz=[T_x T_y 1];
482GeometryCalib.R=[a_X1(2),a_X1(3),0;a_Y1(2),a_Y1(3),0;0,0,1];
483
484% %check error
485% GeometryCalib.ErrorRms(1)=sqrt(mean((x1-x_ima).*(x1-x_ima)));
486% GeometryCalib.ErrorMax(1)=max(abs(x1-x_ima));
487% GeometryCalib.ErrorRms(2)=sqrt(mean((y1-y_ima).*(y1-y_ima)));
488% GeometryCalib.ErrorMax(2)=max(abs(y1-y_ima));
489
490%------------------------------------------------------------------------
491function GeometryCalib=calib_tsai2(Coord)
492%------------------------------------------------------------------
493path_uvmat=which('uvmat');% check the path detected for source file uvmat
494path_UVMAT=fileparts(path_uvmat); %path to UVMAT
495huvmat=findobj(allchild(0),'Tag','uvmat');
496hhuvmat=guidata(huvmat);
497x_1=Coord(:,4:5)';
498X_1=Coord(:,1:3)';
499n_ima=1;
500% check_cond=0;
501
502
503nx=str2num(get(hhuvmat.npx,'String'));
504ny=str2num(get(hhuvmat.npy,'String'));
505
506
507% est_kc=[1;0;0;0;0];
508est_dist=[0;0;0;0;0];
509run(fullfile(path_UVMAT,'toolbox_calib','go_calib_optim'));
510
511GeometryCalib.CalibrationType='tsai_matlab';
512GeometryCalib.focal=f(2);
513GeometryCalib.dpx_dpy=[1 1];
514GeometryCalib.Cx_Cy=cc';
515GeometryCalib.sx=fc(1)/fc(2);
516GeometryCalib.kappa1=-k(1)/f(2)^2;
517GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function
518GeometryCalib.Tx_Ty_Tz=Tc_1';
519GeometryCalib.R=Rc_1;
520% Calib.dpx=GeometryCalib.dpx_dpy(1);
521% Calib.dpy=GeometryCalib.dpx_dpy(2);
522% Calib.sx=GeometryCalib.sx;
523% Calib.Cx=GeometryCalib.Cx_Cy(1);
524% Calib.Cy=GeometryCalib.Cx_Cy(2);
525% Calib.kappa1=GeometryCalib.kappa1;
526% Calib.f=GeometryCalib.focal;
527% Calib.Tx=GeometryCalib.Tx_Ty_Tz(1);
528% Calib.Ty=GeometryCalib.Tx_Ty_Tz(2);
529% Calib.Tz=GeometryCalib.Tx_Ty_Tz(3);
530% Calib.R=GeometryCalib.R;
531% X=Coord(:,1);
532% Y=Coord(:,2);
533% Z=Coord(:,3);
534% x_ima=Coord(:,4);
535% y_ima=Coord(:,5);
536% [Xpoints,Ypoints]=px_XYZ(Calib,X,Y,Z);
537%
538% GeometryCalib.ErrorRms(1)=sqrt(mean((Xpoints-x_ima).*(Xpoints-x_ima)));
539% [GeometryCalib.ErrorMax(1),GeometryCalib.IndexMax(1)]=max(abs(Xpoints-x_ima));
540% GeometryCalib.ErrorRms(2)=sqrt(mean((Ypoints-y_ima).*(Ypoints-y_ima)));
541% [GeometryCalib.ErrorMax(2),GeometryCalib.IndexMax(2)]=max(abs(Ypoints-y_ima));
542
543function GeometryCalib=calib_tsai(Coord)
544%------------------------------------------------------------------------
545%TSAI
546% 'calibration_lin' provides a linear transform on coordinates,
547path_uvmat=which('uvmat');% check the path detected for source file uvmat
548path_UVMAT=fileparts(path_uvmat); %path to UVMAT
549% if isunix
550    %fid = fopen(fullfile(path_UVMAT,'PARAM_LINUX.txt'),'r');%open the file with civ binary names
551xmlfile=fullfile(path_UVMAT,'PARAM.xml');
552if exist(xmlfile,'file')
553    t=xmltree(xmlfile);
554    sparam=convert(t);
555end
556if ~isfield(sparam,'GeometryCalibBin')
557    msgbox_uvmat('ERROR',['calibration program <GeometryCalibBin> undefined in parameter file ' xmlfile])
558    return
559end
560Tsai_exe=sparam.GeometryCalibBin;
561if ~exist(Tsai_exe,'file')%the binary is defined in /bin, default setting
562     Tsai_exe=fullfile(path_UVMAT,Tsai_exe);
563end
564if ~exist(Tsai_exe,'file')
565    msgbox_uvmat('ERROR',['calibration program ' sparam.GeometryCalibBin ' defined in PARAM.xml does not exist'])
566    return
567end
568
569textcoord=num2str(Coord,4);
570dlmwrite('t.txt',textcoord,''); 
571% ['!' Tsai_exe ' -f1 0 -f2 t.txt']
572    eval(['!' Tsai_exe ' -f t.txt > tsaicalib.log']);
573if ~exist('calib.dat','file')
574    msgbox_uvmat('ERROR','no output from calibration program Tsai_exe: possibly too few points')
575end
576calibdat=dlmread('calib.dat');
577delete('calib.dat')
578delete('t.txt')
579GeometryCalib.CalibrationType='tsai';
580GeometryCalib.focal=calibdat(10);
581GeometryCalib.dpx_dpy=[calibdat(5) calibdat(6)];
582GeometryCalib.Cx_Cy=[calibdat(7) calibdat(8)];
583GeometryCalib.sx=calibdat(9);
584GeometryCalib.kappa1=calibdat(11);
585GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function
586GeometryCalib.Tx_Ty_Tz=[calibdat(12) calibdat(13) calibdat(14)];
587Rx_Ry_Rz=calibdat([15:17]);
588sa = sin(Rx_Ry_Rz(1)) ;
589ca=cos(Rx_Ry_Rz(1));
590sb=sin(Rx_Ry_Rz(2));
591cb =cos(Rx_Ry_Rz(2));
592sg =sin(Rx_Ry_Rz(3));
593cg =cos(Rx_Ry_Rz(3));
594r1 = cb * cg;
595r2 = cg * sa * sb - ca * sg;
596r3 = sa * sg + ca * cg * sb;
597r4 = cb * sg;
598r5 = sa * sb * sg + ca * cg;
599r6 = ca * sb * sg - cg * sa;
600r7 = -sb;
601r8 = cb * sa;
602r9 = ca * cb;
603%EN DEDUIRE MATRICE R ??
604GeometryCalib.R=[r1,r2,r3;r4,r5,r6;r7,r8,r9];
605%erreur a caracteriser?
606% %check error
607% Calib.dpx=GeometryCalib.dpx_dpy(1);
608% Calib.dpy=GeometryCalib.dpx_dpy(2);
609% Calib.sx=GeometryCalib.sx;
610% Calib.Cx=GeometryCalib.Cx_Cy(1);
611% Calib.Cy=GeometryCalib.Cx_Cy(2);
612% Calib.kappa1=GeometryCalib.kappa1;
613% Calib.f=GeometryCalib.focal;
614% Calib.Tx=GeometryCalib.Tx_Ty_Tz(1);
615% Calib.Ty=GeometryCalib.Tx_Ty_Tz(2);
616% Calib.Tz=GeometryCalib.Tx_Ty_Tz(3);
617% Calib.R=GeometryCalib.R;
618% X=Coord(:,1);
619% Y=Coord(:,2);
620% Z=Coord(:,3);
621% x_ima=Coord(:,4);
622% y_ima=Coord(:,5);
623% [Xpoints,Ypoints]=px_XYZ(Calib,X,Y,Z);
624%
625% GeometryCalib.ErrorRms(1)=sqrt(mean((Xpoints-x_ima).*(Xpoints-x_ima)));
626% GeometryCalib.ErrorMax(1)=max(abs(Xpoints-x_ima));
627% GeometryCalib.ErrorRms(2)=sqrt(mean((Ypoints-y_ima).*(Ypoints-y_ima)));
628% GeometryCalib.ErrorMax(2)=max(abs(Ypoints-y_ima));
629
630
631%------------------------------------------------------------------------
632% --- Executes on button press in rotation.
633function rotation_Callback(hObject, eventdata, handles)
634%------------------------------------------------------------------------
635angle_rot=(pi/180)*str2num(get(handles.Phi,'String'));
636Coord_cell=get(handles.ListCoord,'String');
637data=read_geometry_calib(Coord_cell);
638data.Coord(:,1)=cos(angle_rot)*data.Coord(:,1)+sin(angle_rot)*data.Coord(:,2);
639data.Coord(:,1)=-sin(angle_rot)*data.Coord(:,1)+cos(angle_rot)*data.Coord(:,2);
640set(handles.XObject,'String',num2str(data.Coord(:,1),4));
641set(handles.YObject,'String',num2str(data.Coord(:,2),4));
642
643%------------------------------------------------------------------------
644function XImage_Callback(hObject, eventdata, handles)
645%------------------------------------------------------------------------
646update_list(hObject, eventdata,handles)
647
648%------------------------------------------------------------------------
649function YImage_Callback(hObject, eventdata, handles)
650%------------------------------------------------------------------------
651update_list(hObject, eventdata,handles)
652
653function XObject_Callback(hObject, eventdata, handles)
654update_list(hObject, eventdata,handles)
655
656function YObject_Callback(hObject, eventdata, handles)
657update_list(hObject, eventdata,handles)
658
659function ZObject_Callback(hObject, eventdata, handles)
660update_list(hObject, eventdata,handles)
661
662%------------------------------------------------------------------------
663function update_list(hObject, eventdata, handles)
664%------------------------------------------------------------------------
665str4=get(handles.XImage,'String');
666str5=get(handles.YImage,'String');
667str1=get(handles.XObject,'String');
668tt=double(str1);
669str2=get(handles.YObject,'String');
670str3=get(handles.ZObject,'String');
671if ~isempty(str1) & ~isequal(double(str1),32) & (isempty(str3)|isequal(double(str3),32))
672    str3='0';%put z to 0 by default
673end
674strline=[str1 '    |    ' str2 '    |    ' str3 '    |    ' str4 '    |    ' str5];
675Coord=get(handles.ListCoord,'String');
676val=get(handles.ListCoord,'Value');
677Coord{val}=strline;
678set(handles.ListCoord,'String',Coord)
679%update the plot
680ListCoord_Callback(hObject, eventdata, handles)
681MenuPlot_Callback(hObject, eventdata, handles)
682
683%------------------------------------------------------------------------
684% --- Executes on selection change in ListCoord.
685function ListCoord_Callback(hObject, eventdata, handles)
686%------------------------------------------------------------------------
687huvmat=findobj(allchild(0),'Name','uvmat');%find the current uvmat interface handle
688hplot=findobj(huvmat,'Tag','axes3');%main plotting axis of uvmat
689hhh=findobj(hplot,'Tag','calib_marker');
690Coord_cell=get(handles.ListCoord,'String');
691val=get(handles.ListCoord,'Value');
692if numel(val)>1
693    return %no action if several lines have been selected
694end
695coord_str=Coord_cell{val};
696k=findstr('|',coord_str);
697if isempty(k)%last line '.....' selected
698    if ~isempty(hhh)
699        delete(hhh)%delete the circle marker
700    end
701    return
702end
703%fill the edit boxex
704set(handles.XObject,'String',coord_str(1:k(1)-5))
705set(handles.YObject,'String',coord_str(k(1)+5:k(2)-5))
706set(handles.ZObject,'String',coord_str(k(2)+5:k(3)-5))
707set(handles.XImage,'String',coord_str(k(3)+5:k(4)-5))
708set(handles.YImage,'String',coord_str(k(4)+5:end))
709h_menu_coord=findobj(huvmat,'Tag','transform_fct');
710menu=get(h_menu_coord,'String');
711choice=get(h_menu_coord,'Value');
712if iscell(menu)
713    option=menu{choice};
714else
715    option='px'; %default
716end
717if isequal(option,'phys')
718    XCoord=str2num(coord_str(1:k(1)-5));
719    YCoord=str2num(coord_str(k(1)+5:k(2)-5));
720elseif isequal(option,'px')|| isequal(option,'')
721    XCoord=str2num(coord_str(k(3)+5:k(4)-5));
722    YCoord=str2num(coord_str(k(4)+5:end));
723else
724    msgbox_uvmat('ERROR','the choice in menu_coord of uvmat must be px or phys ')
725end
726if isempty(XCoord)||isempty(YCoord)
727     if ~isempty(hhh)
728        delete(hhh)%delete the circle marker
729    end
730    return
731end
732xlim=get(hplot,'XLim');
733ylim=get(hplot,'YLim');
734ind_range=max(abs(xlim(2)-xlim(1)),abs(ylim(end)-ylim(1)))/20;%defines the size of the circle marker
735if isempty(hhh)
736    axes(hplot)
737    rectangle('Curvature',[1 1],...
738              'Position',[XCoord-ind_range/2 YCoord-ind_range/2 ind_range ind_range],'EdgeColor','m',...
739              'LineStyle','-','Tag','calib_marker');
740else
741    set(hhh,'Position',[XCoord-ind_range/2 YCoord-ind_range/2 ind_range ind_range])
742end
743
744%------------------------------------------------------------------------
745% --- Executes on selection change in edit_append.
746function edit_append_Callback(hObject, eventdata, handles)
747%------------------------------------------------------------------------
748choice=get(handles.edit_append,'Value');
749% if choice==1
750%        Coord=get(handles.ListCoord,'String');
751%        val=length(Coord);
752%        if val>=1 & isequal(Coord{val},'')
753%             val=val-1; %do not take into account blank
754%        end
755%        Coord{val+1}='';
756%        set(handles.ListCoord,'String',Coord)
757%        set(handles.ListCoord,'Value',val+1)
758% end
759if choice
760    set(handles.edit_append,'BackgroundColor',[1 1 0])
761else
762    set(handles.edit_append,'BackgroundColor',[0.7 0.7 0.7])
763end
764   
765function NEW_Callback(hObject, eventdata, handles)
766%A METTRE SOUS UN BOUTON
767huvmat=findobj(allchild(0),'Name','uvmat');
768hchild=get(huvmat,'children');
769hcoord=findobj(hchild,'Tag','menu_coord');
770coordtype=get(hcoord,'Value');
771haxes=findobj(hchild,'Tag','axes3');
772AxeData=get(haxes,'UserData');
773if ~isequal(hcoord,2)
774    set(hcoord,'Value',2)
775    huvmat=uvmat(AxeData);
776    'relancer uvmat';
777end
778if ~isfield(AxeData,'ZoomAxes')
779    msgbox_uvmat('ERROR','first draw a window around a grid marker')
780    return
781end
782XLim=get(AxeData.ZoomAxes,'XLim');
783YLim=get(AxeData.ZoomAxes,'YLim');
784np=size(AxeData.A);
785ind_sub_x=round(XLim);
786ind_sub_y=np(1)-round(YLim);
787Mfiltre=AxeData.A([ind_sub_y(2):ind_sub_y(1)] ,ind_sub_x,:);
788Mfiltre_norm=double(Mfiltre);
789Mfiltre_norm=Mfiltre_norm/sum(sum(Mfiltre_norm));
790Mfiltre_norm=100*(Mfiltre_norm-mean(mean(Mfiltre_norm)));
791Atype=class(AxeData.A);
792Data.NbDim=2;
793Data.A=filter2(Mfiltre_norm,double(AxeData.A));
794Data.A=feval(Atype,Data.A);
795Data.AName='image';
796Data.AX=AxeData.AX;
797Data.AY=AxeData.AY;
798Data.CoordType='px';
799plot_field(Data)
800
801%------------------------------------------------------------------------
802% --- 'key_press_fcn:' function activated when a key is pressed on the keyboard
803function key_press_fcn(hObject,eventdata,handles)
804%------------------------------------------------------------------------
805xx=double(get(handles.geometry_calib,'CurrentCharacter')); %get the keyboard character
806if ismember(xx,[8 127])%backspace or delete
807    Coord_cell=get(handles.ListCoord,'String');
808    val=get(handles.ListCoord,'Value');
809     if max(val)<numel(Coord_cell) % the last element '...' has not been selected
810        Coord_cell(val)=[];%remove the selected line
811        set(handles.ListCoord,'Value',min(val))
812        set(handles.ListCoord,'String',Coord_cell)         
813        ListCoord_Callback(hObject, eventdata, handles)
814        MenuPlot_Callback(hObject,eventdata,handles)
815     end
816end
817
818% %------------------------------------------------------------------------
819% % --- Executes on button press in append_point.
820% function append_point_Callback(hObject, eventdata, handles)
821% %------------------------------------------------------------------------
822%        Coord=get(handles.ListCoord,'String');
823%        val=length(Coord);
824%        if val>=1 & isequal(Coord{val},'')
825%             val=val-1; %do not take into account blank
826%        end
827%        Coord{val+1}='';
828%        set(handles.ListCoord,'String',Coord)
829%        set(handles.ListCoord,'Value',val+1)
830
831%------------------------------------------------------------------------
832function MenuOpen_Callback(hObject, eventdata, handles)
833%------------------------------------------------------------------------
834%get the object file
835huvmat=findobj(allchild(0),'Name','uvmat');
836UvData=get(huvmat,'UserData');
837hchild=get(huvmat,'Children');
838hrootpath=findobj(hchild,'Tag','RootPath');
839oldfile=get(hrootpath,'String');
840if isempty(oldfile)
841    oldfile='';
842end
843%[FileName,PathName] = uigetfile('*.civ','Select a .civ file',oldfile)
844[FileName, PathName, filterindex] = uigetfile( ...
845       {'*.xml;*.mat', ' (*.xml,*.mat)';
846       '*.xml',  '.xml files '; ...
847        '*.mat',  '.mat matlab files '}, ...
848        'Pick a file',oldfile);
849fileinput=[PathName FileName];%complete file name
850testblank=findstr(fileinput,' ');%look for blanks
851if ~isempty(testblank)
852    msgbox_uvmat('ERROR','forbidden input file name or path: no blank character allowed')
853    return
854end
855sizf=size(fileinput);
856if (~ischar(fileinput)|~isequal(sizf(1),1)),return;end
857loadfile(handles,fileinput)
858
859
860%------------------------------------------------------------------------
861function MenuPlot_Callback(hObject, eventdata, handles)
862%------------------------------------------------------------------------
863huvmat=findobj(allchild(0),'Name','uvmat');%find the current uvmat interface handle
864UvData=get(huvmat,'UserData');%Data associated to the current uvmat interface
865hhuvmat=guidata(huvmat); %handles of GUI elements in uvmat
866hplot=findobj(huvmat,'Tag','axes3');%main plotting axis of uvmat
867h_menu_coord=findobj(huvmat,'Tag','transform_fct');
868menu=get(h_menu_coord,'String');
869choice=get(h_menu_coord,'Value');
870if iscell(menu)
871    option=menu{choice};
872else
873    option='px'; %default
874end
875Coord_cell=get(handles.ListCoord,'String');
876ObjectData=read_geometry_calib(Coord_cell);
877%ObjectData=read_geometry_calib(handles);%read the interface input parameters defining the object
878if ~isempty(ObjectData.Coord)
879    if isequal(option,'phys')
880        ObjectData.Coord=ObjectData.Coord(:,[1:3]);
881    elseif isequal(option,'px')||isequal(option,'')
882        ObjectData.Coord=ObjectData.Coord(:,[4:5]);
883    else
884        msgbox_uvmat('ERROR','the choice in menu_coord of uvmat must be '''', px or phys ')
885    end
886end
887axes(hhuvmat.axes3)
888hh=findobj('Tag','calib_points');
889if  ~isempty(ObjectData.Coord) && isempty(hh)
890    hh=line(ObjectData.Coord(:,1),ObjectData.Coord(:,2),'Color','m','Tag','calib_points','LineStyle','.','Marker','+');
891elseif isempty(ObjectData.Coord)%empty list of points, suppress the plot
892    delete(hh)
893else
894    set(hh,'XData',ObjectData.Coord(:,1))
895    set(hh,'YData',ObjectData.Coord(:,2))
896end
897pause(.1)
898figure(handles.geometry_calib)
899
900% --------------------------------------------------------------------
901function MenuHelp_Callback(hObject, eventdata, handles)
902path_to_uvmat=which ('uvmat');% check the path of uvmat
903pathelp=fileparts(path_to_uvmat);
904    helpfile=fullfile(pathelp,'uvmat_doc','uvmat_doc.html');
905if isempty(dir(helpfile)), msgbox_uvmat('ERROR','Please put the help file uvmat_doc.html in the sub-directory /uvmat_doc of the UVMAT package')
906else
907   addpath (fullfile(pathelp,'uvmat_doc'))
908   web([helpfile '#geometry_calib'])
909end
910
911%------------------------------------------------------------------------
912function MenuCreateGrid_Callback(hObject, eventdata, handles)
913%------------------------------------------------------------------------
914%hcalib=get(handles.calib_type,'parent');%handles of the GUI geometry_calib
915CalibData=get(handles.geometry_calib,'UserData');
916Tinput=[];%default
917if isfield(CalibData,'grid')
918    Tinput=CalibData.grid;
919end
920[T,CalibData.grid]=create_grid(Tinput);%display the GUI create_grid
921set(handles.geometry_calib,'UserData',CalibData)
922
923%grid in phys space
924Coord=get(handles.ListCoord,'String');
925val=get(handles.ListCoord,'Value');
926data=read_geometry_calib(Coord);
927%nbpoints=size(data.Coord,1); %nbre of calibration points
928data.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
929% for i=1:nbpoints
930%    for j=1:5
931%           Coord{i,j}=num2str(data.Coord(i,j),4);%display coordiantes with 4 digits
932%    end
933% end
934%update the phys coordinates starting from the selected point (down in the
935Coord(end,:)=[]; %remove last string '.....'
936for i=1:size(data.Coord,1)
937    for j=1:5
938          Coord{i,j}=num2str(data.Coord(i,j),4);%display coordiantes with 4 digits
939    end
940end
941
942%size(data.Coord,1)
943Tabchar=cell2tab(Coord,'    |    ');
944Tabchar=[Tabchar ;{'......'}];
945set(handles.ListCoord,'String',Tabchar)
946
947% -----------------------------------------------------------------------
948% --- automatic grid dectection from local maxima of the images
949function MenuDetectGrid_Callback(hObject, eventdata, handles)
950%------------------------------------------------------------------------
951CalibData=get(handles.geometry_calib,'UserData');
952grid_input=[];%default
953if isfield(CalibData,'grid')
954    grid_input=CalibData.grid;%retrieve the previously used grid
955end
956[T,CalibData.grid]=create_grid(grid_input,'detect grid');%display the GUI create_grid
957set(handles.geometry_calib,'UserData',CalibData)%store the phys grid for later use
958
959%read the four last point coordiantes in pixels
960Coord_cell=get(handles.ListCoord,'String');%read list of coordiantes on geometry_calib
961data=read_geometry_calib(Coord_cell);
962nbpoints=size(data.Coord,1); %nbre of calibration points
963if nbpoints~=4
964    msgbox_uvmat('ERROR','four points must be selected by the mouse, beginning by the new x axis, to delimitate the phs grid area')
965    return
966end
967corners_X=(data.Coord(end:-1:end-3,4)); %pixel absissa of the four corners
968corners_Y=(data.Coord(end:-1:end-3,5));
969
970%reorder the last two points (the two first in the list) if needed
971angles=angle((corners_X-corners_X(1))+i*(corners_Y-corners_Y(1)));
972if abs(angles(4)-angles(2))>abs(angles(3)-angles(2))
973      X_end=corners_X(4);
974      Y_end=corners_Y(4);
975      corners_X(4)=corners_X(3);
976      corners_Y(4)=corners_Y(3);
977      corners_X(3)=X_end;
978      corners_Y(3)=Y_end;
979end
980
981%read the current image
982huvmat=findobj(allchild(0),'Name','uvmat');
983UvData=get(huvmat,'UserData');
984A=UvData.Field.A;
985npxy=size(A);
986%linear transform on the current image
987X=[CalibData.grid.x_0 CalibData.grid.x_1 CalibData.grid.x_0 CalibData.grid.x_1]';%corner absissa in the rectified image
988Y=[CalibData.grid.y_0 CalibData.grid.y_0 CalibData.grid.y_1 CalibData.grid.y_1]';%corner absissa in the rectified image
989XY_mat=[ones(size(X)) X Y];
990a_X1=XY_mat\corners_X; %transformation matrix for X
991x1=XY_mat*a_X1;%reconstruction
992err_X1=max(abs(x1-corners_X));%error
993a_Y1=XY_mat\corners_Y;%transformation matrix for X
994y1=XY_mat*a_Y1;
995err_Y1=max(abs(y1-corners_Y));%error
996GeometryCalib.CalibrationType='linear';
997GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function
998GeometryCalib.f=1;
999GeometryCalib.dpx=1;
1000GeometryCalib.dpy=1;
1001GeometryCalib.sx=1;
1002GeometryCalib.Cx=0;
1003GeometryCalib.Cy=0;
1004GeometryCalib.kappa1=0;
1005GeometryCalib.Tx=a_X1(1);
1006GeometryCalib.Ty=a_Y1(1);
1007GeometryCalib.Tz=1;
1008GeometryCalib.R=[a_X1(2),a_X1(3),0;a_Y1(2),a_Y1(3),0;0,0,1];
1009[Amod,Rangx,Rangy]=phys_Ima(A-min(min(A)),GeometryCalib,0);
1010Amod=double(Amod);
1011%figure(12)
1012%Amax=max(max(Amod))
1013%image(Rangx,Rangy,uint8(255*Amod/Amax))
1014Dx=(Rangx(2)-Rangx(1))/(npxy(2)-1); %x mesh in real space
1015Dy=(Rangy(2)-Rangy(1))/(npxy(1)-1); %y mesh in real space
1016ind_range_x=ceil(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
1017ind_range_y=ceil(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
1018nbpoints=size(T,1);
1019for ipoint=1:nbpoints
1020    i0=1+round((T(ipoint,1)-Rangx(1))/Dx);%round(Xpx(ipoint));
1021    j0=1+round((T(ipoint,2)-Rangy(1))/Dy);%round(Xpx(ipoint));
1022    Asub=Amod(j0-ind_range_y:j0+ind_range_y,i0-ind_range_x:i0+ind_range_x);
1023    x_profile=sum(Asub,1);
1024    y_profile=sum(Asub,2);
1025    [Amax,ind_x_max]=max(x_profile);
1026    [Amax,ind_y_max]=max(y_profile);
1027    %sub-pixel improvement using moments
1028    x_shift=0;
1029    y_shift=0;
1030    if ind_x_max+2<=2*ind_range_x+1 && ind_x_max-2>=1
1031        Atop=x_profile(ind_x_max-2:ind_x_max+2);
1032        x_shift=sum(Atop.*[-2 -1 0 1 2])/sum(Atop);
1033    end
1034    if ind_y_max+2<=2*ind_range_y+1 && ind_y_max-2>=1
1035        Atop=y_profile(ind_y_max-2:ind_y_max+2);
1036        y_shift=sum(Atop.*[-2 -1 0 1 2]')/sum(Atop);
1037    end
1038    Delta(ipoint,1)=(x_shift+ind_x_max-ind_range_x-1)*Dx;%shift from the initial guess
1039    Delta(ipoint,2)=(y_shift+ind_y_max-ind_range_y-1)*Dy;
1040end
1041Tmod=T(:,(1:2))+Delta;
1042[Xpx,Ypx]=px_XYZ(GeometryCalib,Tmod(:,1),Tmod(:,2));
1043for ipoint=1:nbpoints
1044     Coord{ipoint,1}=num2str(T(ipoint,1),4);%display coordiantes with 4 digits
1045     Coord{ipoint,2}=num2str(T(ipoint,2),4);%display coordiantes with 4 digits
1046     Coord{ipoint,3}='0';
1047     Coord{ipoint,4}=num2str(Xpx(ipoint),4);%display coordiantes with 4 digi
1048     Coord{ipoint,5}=num2str(Ypx(ipoint),4);%display coordiantes with 4 digi
1049end
1050Tabchar=cell2tab(Coord(end:-1:1,:),'    |    ');
1051Tabchar=[Tabchar ;{'......'}];
1052set(handles.ListCoord,'Value',1)
1053set(handles.ListCoord,'String',Tabchar)
1054MenuPlot_Callback(hObject, eventdata, handles)
1055
1056%-----------------------------------------------------------------------
1057function MenuTranslatePoints_Callback(hObject, eventdata, handles)
1058%-----------------------------------------------------------------------
1059%hcalib=get(handles.calib_type,'parent');%handles of the GUI geometry_calib
1060CalibData=get(handles.geometry_calib,'UserData');
1061Tinput=[];%default
1062if isfield(CalibData,'translate')
1063    Tinput=CalibData.translate;
1064end
1065T=translate_points(Tinput);%display translate_points GUI and get shift parameters
1066CalibData.translate=T;
1067set(handles.geometry_calib,'UserData',CalibData)
1068%translation
1069Coord_cell=get(handles.ListCoord,'String');
1070data=read_geometry_calib(Coord_cell);
1071data.Coord(:,1)=T(1)+data.Coord(:,1);
1072data.Coord(:,2)=T(2)+data.Coord(:,2);
1073data.Coord(:,3)=T(3)+data.Coord(:,3);
1074data.Coord(:,[4 5])=data.Coord(:,[4 5]);
1075for i=1:size(data.Coord,1)
1076    for j=1:5
1077          Coord{i,j}=num2str(data.Coord(i,j),4);%phys x,y,z
1078   end
1079end
1080Tabchar=cell2tab(Coord,'    |    ');
1081Tabchar=[Tabchar; {'.....'}];
1082%set(handles.ListCoord,'Value',1)
1083set(handles.ListCoord,'String',Tabchar)
1084
1085
1086% --------------------------------------------------------------------
1087function MenuRotatePoints_Callback(hObject, eventdata, handles)
1088%hcalib=get(handles.calib_type,'parent');%handles of the GUI geometry_calib
1089CalibData=get(handles.geometry_calib,'UserData');
1090Tinput=[];%default
1091if isfield(CalibData,'rotate')
1092    Tinput=CalibData.rotate;
1093end
1094T=rotate_points(Tinput);%display translate_points GUI and get shift parameters
1095CalibData.rotate=T;
1096set(handles.geometry_calib,'UserData',CalibData)
1097%-----------------------------------------------------
1098%rotation
1099Phi=T(1);
1100O_x=0;%default
1101O_y=0;%default
1102if numel(T)>=2
1103    O_x=T(2);%default
1104end
1105if numel(T)>=3
1106    O_y=T(3);%default
1107end
1108Coord_cell=get(handles.ListCoord,'String');
1109data=read_geometry_calib(Coord_cell);
1110r1=cos(pi*Phi/180);
1111r2=-sin(pi*Phi/180);
1112r3=sin(pi*Phi/180);
1113r4=cos(pi*Phi/180);
1114x=data.Coord(:,1)-O_x;
1115y=data.Coord(:,2)-O_y;
1116data.Coord(:,1)=r1*x+r2*y;
1117data.Coord(:,2)=r3*x+r4*y;
1118% data.Coord(:,[4 5])=data.Coord(:,[4 5]);
1119for i=1:size(data.Coord,1)
1120    for j=1:5
1121          Coord{i,j}=num2str(data.Coord(i,j),4);%phys x,y,z
1122   end
1123end
1124Tabchar=cell2tab(Coord,'    |    ');
1125Tabchar=[Tabchar;{'......'}];
1126set(handles.ListCoord,'Value',1)
1127set(handles.ListCoord,'String',Tabchar)
1128
1129
1130%%%%%%%%%%%%%%%%%%%%
1131function [A_out,Rangx,Rangy]=phys_Ima(A,Calib,ZIndex)
1132xcorner=[];
1133ycorner=[];
1134npx=[];
1135npy=[];
1136siz=size(A);
1137npx=[npx siz(2)];
1138npy=[npy siz(1)];
1139xima=[0.5 siz(2)-0.5 0.5 siz(2)-0.5];%image coordiantes of corners
1140yima=[0.5 0.5 siz(1)-0.5 siz(1)-0.5];
1141[xcorner,ycorner]=phys_XYZ(Calib,xima,yima,ZIndex);%corresponding physical coordinates
1142Rangx(1)=min(xcorner);
1143Rangx(2)=max(xcorner);
1144Rangy(2)=min(ycorner);
1145Rangy(1)=max(ycorner);
1146test_multi=(max(npx)~=min(npx)) | (max(npy)~=min(npy));
1147npx=max(npx);
1148npy=max(npy);
1149x=linspace(Rangx(1),Rangx(2),npx);
1150y=linspace(Rangy(1),Rangy(2),npy);
1151[X,Y]=meshgrid(x,y);%grid in physical coordiantes
1152vec_B=[];
1153
1154zphys=0; %default
1155if isfield(Calib,'SliceCoord') %.Z= index of plane
1156   SliceCoord=Calib.SliceCoord(ZIndex,:);
1157   zphys=SliceCoord(3); %to generalize for non-parallel planes
1158end
1159[XIMA,YIMA]=px_XYZ(Calib,X,Y,zphys);%corresponding image indices for each point in the real space grid
1160XIMA=reshape(round(XIMA),1,npx*npy);%indices reorganized in 'line'
1161YIMA=reshape(round(YIMA),1,npx*npy);
1162flagin=XIMA>=1 & XIMA<=npx & YIMA >=1 & YIMA<=npy;%flagin=1 inside the original image
1163testuint8=isa(A,'uint8');
1164testuint16=isa(A,'uint16');
1165if numel(siz)==2 %(B/W images)
1166    vec_A=reshape(A,1,npx*npy);%put the original image in line
1167    ind_in=find(flagin);
1168    ind_out=find(~flagin);
1169    ICOMB=((XIMA-1)*npy+(npy+1-YIMA));
1170    ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A
1171    vec_B(ind_in)=vec_A(ICOMB);
1172    vec_B(ind_out)=zeros(size(ind_out));
1173    A_out=reshape(vec_B,npy,npx);%new image in real coordinates
1174elseif numel(siz)==3     
1175    for icolor=1:siz(3)
1176        vec_A=reshape(A{icell}(:,:,icolor),1,npx*npy);%put the original image in line
1177        ind_in=find(flagin);
1178        ind_out=find(~flagin);
1179        ICOMB=((XIMA-1)*npy+(npy+1-YIMA));
1180        ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A
1181        vec_B(ind_in)=vec_A(ICOMB);
1182        vec_B(ind_out)=zeros(size(ind_out));
1183        A_out(:,:,icolor)=reshape(vec_B,npy,npx);%new image in real coordinates
1184    end
1185end
1186if testuint8
1187    A_out=uint8(A_out);
1188end
1189if testuint16
1190    A_out=uint16(A_out);
1191end
1192
1193%INPUT:
1194%Z: index of plane
1195function [Xphys,Yphys,Zphys]=phys_XYZ(Calib,X,Y,Z)
1196if exist('Z','var')& isequal(Z,round(Z))& Z>0 & isfield(Calib,'SliceCoord')&length(Calib.SliceCoord)>=Z
1197    Zindex=Z;
1198    Zphys=Calib.SliceCoord(Zindex,3);%GENERALISER AUX CAS AVEC ANGLE
1199else
1200%     if exist('Z','var')
1201%         Zphys=Z;
1202%     else
1203        Zphys=0;
1204%     end
1205end
1206if ~exist('X','var')||~exist('Y','var')
1207    Xphys=[];
1208    Yphys=[];%default
1209    return
1210end
1211Xphys=X;%default
1212Yphys=Y;
1213%image transform
1214if isfield(Calib,'R')
1215    R=(Calib.R)';
1216    Dx=R(5)*R(7)-R(4)*R(8);
1217    Dy=R(1)*R(8)-R(2)*R(7);
1218    D0=Calib.f*(R(2)*R(4)-R(1)*R(5));
1219    Z11=R(6)*R(8)-R(5)*R(9);
1220    Z12=R(2)*R(9)-R(3)*R(8); 
1221    Z21=R(4)*R(9)-R(6)*R(7);
1222    Z22=R(3)*R(7)-R(1)*R(9);
1223    Zx0=R(3)*R(5)-R(2)*R(6);
1224    Zy0=R(1)*R(6)-R(3)*R(4);
1225    A11=R(8)*Calib.Ty-R(5)*Calib.Tz+Z11*Zphys;
1226    A12=R(2)*Calib.Tz-R(8)*Calib.Tx+Z12*Zphys;
1227    A21=-R(7)*Calib.Ty+R(4)*Calib.Tz+Z21*Zphys;
1228    A22=-R(1)*Calib.Tz+R(7)*Calib.Tx+Z11*Zphys;
1229    X0=Calib.f*(R(5)*Calib.Tx-R(2)*Calib.Ty+Zx0*Zphys);
1230    Y0=Calib.f*(-R(4)*Calib.Tx+R(1)*Calib.Ty+Zy0*Zphys);
1231        %px to camera:
1232    Xd=(Calib.dpx/Calib.sx)*(X-Calib.Cx); % sensor coordinates
1233    Yd=Calib.dpy*(Y-Calib.Cy);
1234    dist_fact=1+Calib.kappa1*(Xd.*Xd+Yd.*Yd); %distortion factor
1235    Xu=dist_fact.*Xd;%undistorted sensor coordinates
1236    Yu=dist_fact.*Yd;
1237    denom=Dx*Xu+Dy*Yu+D0;
1238    % denom2=denom.*denom;
1239    Xphys=(A11.*Xu+A12.*Yu+X0)./denom;%world coordinates
1240    Yphys=(A21.*Xu+A22.*Yu+Y0)./denom;
1241end
Note: See TracBrowser for help on using the repository browser.