# source:trunk/src/geometry_calib.m@78

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

debugging of geometry_calib and mouse operations
merge_proj: call to read_image replaced. Still problems to solve to merge images

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