source: trunk/src/geometry_calib.m @ 69

Last change on this file since 69 was 69, checked in by gostiaux, 14 years ago

Matlab tsai calibration included

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