source: trunk/src/geometry_calib.m @ 67

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

civ: RUN civ lounched out of the Matlab work space. RUN and BATCH now runned by a unique sub-function lounch.m.
FiLE PARAM.xml modified to provide different paths for Batch and Run !!!!!
RUN_FIX: minor error message modif
geometry_calib: calib point editing by the mouse
improvement of the interactions between the different GUIs. Close function...

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