source: trunk/src/geometry_calib.m @ 71

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

civ3D updated: introduction of image size
imadoc2struct: reding of image size from the xml file
set_object, view_field and related functions: improvement of projection object editing
mouse: possibility of adjusting the calibrations points with the mouse

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