source: trunk/src/geometry_calib.m @ 60

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

-plot projections on a new specific GUI view_field, to avoid the multiplication of active figures
-introduce a ruler in uvmat (in menu/Tools), to measure distances and angles
-insert automatic detection of points for geometric calibration: tool 'detect_grid'. Four points, deliminating the grid to determine, must be marked with the mouse, as well as the physical grid. Then the points inside are automatically detected.
-reading plane positions in ImaDoc? improved to deal with volume scans

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