source: trunk/src/geometry_calib.m @ 108

Last change on this file since 108 was 108, checked in by gostiaux, 10 years ago

Changes from Joel on detect_grid

File size: 47.8 KB
Line 
1%'geometry_calib': performs geometric calibration from a set of reference points
2%
3% function varargout = geometry_calib(varargin)
4%
5%A%AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
6%  Copyright Joel Sommeria, 2008, LEGI / CNRS-UJF-INPG, sommeria@coriolis-legi.org.
7%AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
8%     This file is part of the toolbox UVMAT.
9%
10%     UVMAT is free software; you can redistribute it and/or modify
11%     it under the terms of the GNU General Public License as published by
12%     the Free Software Foundation; either version 2 of the License, or
13%     (at your option) any later version.
14%
15%     UVMAT is distributed in the hope that it will be useful,
16%     but WITHOUT ANY WARRANTY; without even the implied warranty of
17%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
18%     GNU General Public License (file UVMAT/COPYING.txt) for more details.
19%AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
20
21function varargout = geometry_calib(varargin)
22% GEOMETRY_CALIB M-file for geometry_calib.fig
23%      GEOMETRY_CALIB, by itself, creates a MenuCoord GEOMETRY_CALIB or raises the existing
24%      singleton*.
25%
26%      H = GEOMETRY_CALIB returns the handle to a MenuCoord GEOMETRY_CALIB or the handle to
27%      the existing singleton*.
28%
29%      GEOMETRY_CALIB('CALLBACK',hObject,eventData,handles,...) calls the local
30%      function named CALLBACK in GEOMETRY_CALIB.M with the given input arguments.
31%
32%      GEOMETRY_CALIB('Property','Value',...) creates a MenuCoord GEOMETRY_CALIB or raises the
33%      existing singleton*.  Starting from the left, property value pairs are
34%      applied to the GUI before geometry_calib_OpeningFunction gets called.  An
35%      unrecognized property name or invalid value makes property application
36%      stop.  All inputs are passed to geometry_calib_OpeningFcn via varargin.
37%
38%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
39%      instance to run (singleton)".
40%
41% See also: GUIDE, GUIDATA, GUIHANDLES
42
43% Edit the above text to modify the response to help geometry_calib
44
45% Last Modified by GUIDE v2.5 25-Mar-2010 19:10:05
46
47% Begin initialization code - DO NOT edit
48gui_Singleton = 1;
49gui_State = struct('gui_Name',       mfilename, ...
50                   'gui_Singleton',  gui_Singleton, ...
51                   'gui_OpeningFcn', @geometry_calib_OpeningFcn, ...
52                   'gui_OutputFcn',  @geometry_calib_OutputFcn, ...
53                   'gui_LayoutFcn',  [] , ...
54                   'gui_Callback',   []);
55if nargin & isstr(varargin{1})
56    gui_State.gui_Callback = str2func(varargin{1});
57end
58
59if nargout
60    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
61else
62    gui_mainfcn(gui_State, varargin{:});
63end
64% End initialization code - DO NOT edit
65
66
67% --- Executes just before geometry_calib is made visible.
68%INPUT:
69%handles: handles of the geometry_calib interface elements
70% PlotHandles: set of handles of the elements contolling the plotting
71% parameters on the uvmat interface (obtained by 'get_plot_handle.m')
72%------------------------------------------------------------------------
73function geometry_calib_OpeningFcn(hObject, eventdata, handles, handles_uvmat,pos,inputfile)
74%------------------------------------------------------------------------
75% Choose default command line output for geometry_calib
76handles.output = hObject;
77
78% Update handles structure
79guidata(hObject, handles);
80set(hObject,'DeleteFcn',{@closefcn})%
81
82%set the position of the interface
83if exist('pos','var')& length(pos)>2
84    pos_gui=get(hObject,'Position');
85    pos_gui(1)=pos(1);
86    pos_gui(2)=pos(2);
87    set(hObject,'Position',pos_gui);
88end
89inputxml='';
90if exist('inputfile','var')& ~isempty(inputfile)
91%     [Path,Name,ext]=fileparts(inputfile);
92%     form=imformats(ext([2:end]));
93%     if ~isempty(form)% if the input file is an image
94        struct.XmlInputfile=inputfile;
95        set(hObject,'UserData',struct)
96        [Pathsub,RootFile,field_count,str2,str_a,str_b,ext,nom_type,subdir]=name2display(inputfile);
97        inputxml=[fullfile(Pathsub,RootFile) '.xml'];
98%     end   
99end
100set(handles.ListCoord,'String',{'......'})
101if exist(inputxml,'file')
102    loadfile(handles,inputxml)% load the point coordiantes existing in the xml file
103end
104
105set(handles.ListCoord,'KeyPressFcn',{@key_press_fcn,handles})%set keyboard action function
106
107
108%------------------------------------------------------------------------
109% --- Outputs from this function are returned to the command line.
110function varargout = geometry_calib_OutputFcn(hObject, eventdata, handles)
111%------------------------------------------------------------------------
112% Get default command line output from handles structure
113varargout{1} = handles.output;
114varargout{2}=handles;
115
116%------------
117function Phi_Callback(hObject, eventdata, handles)
118
119%------------------------------------------------------------------------
120%read input xml file and update the edit boxes
121function loadfile(handles,fileinput)
122%------------------------------------------------------------------------
123%read the input xml file
124t=xmltree(fileinput);
125s=convert(t);%convert to matlab structure
126%read data currently displayed on the interface
127PointCoord=[];
128Coord_cell=get(handles.ListCoord,'String');
129data=read_geometry_calib(Coord_cell);
130%data=read_geometry_calib(handles);
131Coord=[]; %default
132if isfield(data,'Coord')
133    Coord=data.Coord;
134end
135TabChar_0=get(handles.ListCoord,'String');
136nbcoord_0=size(TabChar_0,1);
137if isequal(get(handles.edit_append,'Value'),2) %edit mode  A REVOIR
138    val=get(handles.ListCoord,'Value')-1;
139else
140   val=length(TabChar_0);
141end
142nbcoord=0;
143
144%case of calibration (ImaDoc) input file
145% hcalib=get(handles.calib_type,'parent');
146CalibData=get(handles.geometry_calib,'UserData');
147CalibData.XmlInput=fileinput;
148if isfield(s,'Heading')
149    CalibData.Heading=s.Heading;
150end
151
152set(handles.geometry_calib,'UserData',CalibData);%store the heading in the interface 'UserData'
153if isfield(s,'GeometryCalib')
154    Calib=s.GeometryCalib;
155    if isfield(Calib,'CalibrationType')
156        CalibrationType=Calib.CalibrationType;
157        switch CalibrationType
158            case 'linear'
159                set(handles.calib_type,'Value',2)
160            case 'tsai'
161                set(handles.calib_type,'Value',3)
162            case 'tsai_matlab'
163                set(handles.calib_type,'Value',4)
164        end
165    end
166    if isfield(Calib,'SourceCalib')
167        if isfield(Calib.SourceCalib,'PointCoord')
168            PointCoord=Calib.SourceCalib.PointCoord;
169        end
170    end
171    nbcoord=length(PointCoord);
172    if ~isfield(Calib,'ErrorRms')&~isfield(Calib,'ErrorMax') %old convention of Gauthier (cord in mm)
173        for i=1:length(PointCoord)
174          line=str2num(PointCoord{i});
175          Coord(i+val,4:5)=line(4:5);%px x
176          Coord(i+val,1:3)=line(1:3)/10;%phys x
177        end
178    else
179        for i=1:length(PointCoord)
180          line=str2num(PointCoord{i});
181          Coord(i,4:5)=line(4:5);%px x
182          Coord(i,1:3)=line(1:3);%phys x
183       end
184    end
185end
186%case of xml files of points
187if isfield(s,'Coord')
188    PointCoord=s.Coord;
189    nbcoord=length(PointCoord);
190     %case of image coordinates
191    if isfield(s,'CoordType')& isequal(s.CoordType,'px')
192        for i=1:nbcoord
193           line=str2num(PointCoord{i});
194           Coord(i+val,4:5)=line(1:2);
195        end
196     %case of  physical coordinates
197    else
198        for i=1:nbcoord
199           line=str2num(PointCoord{i});
200           Coord(i+val,1:3)=line(1:3);
201           nbcolumn=size(Coord,2);
202           if nbcolumn<5
203               Coord(i+val,nbcolumn+1:5)=zeros(1,5-nbcolumn);
204           end
205        end
206     end
207end
208CoordCell={};
209for iline=1:size(Coord,1)
210    for j=1:5
211        CoordCell{iline,j}=num2str(Coord(iline,j),4);
212    end
213end
214% CoordCell=[CoordCell;{' ',' ',' ',' ',' '}];
215Tabchar=cell2tab(CoordCell,'    |    ');%transform cells into table ready for display
216Tabchar=[Tabchar;{'......'}];
217set(handles.ListCoord,'Value',1)
218set(handles.ListCoord,'String',Tabchar)
219MenuPlot_Callback(handles.geometry_calib, [], handles)
220if isempty(Coord)
221    set(handles.edit_append,'Value',1)
222    set(handles.edit_append,'BackgroundColor',[1 1 0])
223else
224    set(handles.edit_append,'Value',0)
225    set(handles.edit_append,'BackgroundColor',[0.7 0.7 0.7])
226end
227%
228%------------------------------------------------------------------------
229% executed when closing: set the parent interface button to value 0
230function closefcn(gcbo,eventdata)
231%------------------------------------------------------------------------
232huvmat=findobj(allchild(0),'Name','uvmat');
233if ~isempty(huvmat)
234    handles=guidata(huvmat);
235    set(handles.MenuMask,'enable','on')
236    set(handles.MenuGrid,'enable','on')
237    set(handles.MenuObject,'enable','on')
238    set(handles.MenuEdit,'enable','on')
239    set(handles.edit,'enable','on')
240    hobject=findobj(handles.axes3,'tag','calib_points');
241    if ~isempty(hobject)
242        delete(hobject)
243    end
244    hobject=findobj(handles.axes3,'tag','calib_marker');
245    if ~isempty(hobject)
246        delete(hobject)
247    end   
248end
249
250%------------------------------------------------------------------------
251% --- Executes on button press in calibrate_lin.
252function APPLY_Callback(hObject, eventdata, handles)
253%------------------------------------------------------------------------
254calib_cell=get(handles.calib_type,'String');
255val=get(handles.calib_type,'Value');
256calib_type=calib_cell{val};
257Coord_cell=get(handles.ListCoord,'String');
258Object=read_geometry_calib(Coord_cell);
259X=Object.Coord(:,1);
260Y=Object.Coord(:,2);
261Z=Object.Coord(:,3);
262if isequal(calib_type,'rescale')
263    GeometryCalib=calib_rescale(Object.Coord);
264    Z=0;%Z not taken into account
265elseif isequal(calib_type,'linear')
266    GeometryCalib=calib_linear(Object.Coord);
267    Z=0; %Z not taken into account
268elseif isequal(calib_type,'tsai_cpp')
269    GeometryCalib=calib_tsai(Object.Coord);
270elseif isequal(calib_type,'tsai_matlab')
271    GeometryCalib=calib_tsai2(Object.Coord);
272end
273
274%check error
275if isfield(GeometryCalib,'dpx_dpy')
276    Calib.dpx=GeometryCalib.dpx_dpy(1);
277    Calib.dpy=GeometryCalib.dpx_dpy(2);
278end
279if isfield(GeometryCalib,'sx')
280    Calib.sx=GeometryCalib.sx;
281end
282if isfield(GeometryCalib,'Cx_Cy')
283    Calib.Cx=GeometryCalib.Cx_Cy(1);
284    Calib.Cy=GeometryCalib.Cx_Cy(2);
285end
286if isfield(GeometryCalib,'kappa1')
287    Calib.kappa1=GeometryCalib.kappa1;
288end
289if isfield(GeometryCalib,'focal')
290    Calib.f=GeometryCalib.focal;
291end
292if isfield(GeometryCalib,'Tx_Ty_Tz')
293    Calib.Tx=GeometryCalib.Tx_Ty_Tz(1);
294    Calib.Ty=GeometryCalib.Tx_Ty_Tz(2);
295    Calib.Tz=GeometryCalib.Tx_Ty_Tz(3);
296end
297if isfield(GeometryCalib,'R')
298    Calib.R=GeometryCalib.R;
299end
300x_ima=Object.Coord(:,4);
301y_ima=Object.Coord(:,5);
302[Xpoints,Ypoints]=px_XYZ(Calib,X,Y,Z);
303GeometryCalib.ErrorRms(1)=sqrt(mean((Xpoints-x_ima).*(Xpoints-x_ima)));
304[GeometryCalib.ErrorMax(1),index(1)]=max(abs(Xpoints-x_ima));
305GeometryCalib.ErrorRms(2)=sqrt(mean((Ypoints-y_ima).*(Ypoints-y_ima)));
306[GeometryCalib.ErrorMax(2),index(2)]=max(abs(Ypoints-y_ima));
307[EM,ind_dim]=max(GeometryCalib.ErrorMax);
308index=index(ind_dim);
309
310unitlist=get(handles.CoordUnit,'String');
311unit=unitlist{get(handles.CoordUnit,'value')};
312GeometryCalib.CoordUnit=unit;
313GeometryCalib.SourceCalib.PointCoord=Object.Coord;
314huvmat=findobj(allchild(0),'Name','uvmat');
315hhuvmat=guidata(huvmat);%handles of elements in the GUI uvmat
316RootPath='';
317RootFile='';
318if ~isempty(hhuvmat.RootPath)& ~isempty(hhuvmat.RootFile)
319    testhandle=1;
320    RootPath=get(hhuvmat.RootPath,'String');
321    RootFile=get(hhuvmat.RootFile,'String');
322    filebase=fullfile(RootPath,RootFile);
323    outputfile=[filebase '.xml'];
324else
325    question={'save the calibration data and point coordinates in'};
326    def={fullfile(RootPath,['ObjectCalib.xml'])};
327    options.Resize='on';
328    answer=inputdlg(question,'save average in a new file',1,def,options);
329    outputfile=answer{1};
330end
331answer=msgbox_uvmat('INPUT_Y-N',{[outputfile ' updated with calibration data'];...
332    ['Error rms (along x,y)=' num2str(GeometryCalib.ErrorRms) ' pixels'];...
333    ['Error max (along x,y)=' num2str(GeometryCalib.ErrorMax) ' pixels']});
334if isequal(answer,'Yes')
335    update_imadoc(GeometryCalib,outputfile)
336    %display image with new calibration in the currently opened uvmat interface
337    hhh=findobj(hhuvmat.axes3,'Tag','calib_marker');% delete calib points and markers
338    if ~isempty(hhh)
339        delete(hhh);
340    end
341    hhh=findobj(hhuvmat.axes3,'Tag','calib_points');
342    if ~isempty(hhh)
343        delete(hhh);
344    end
345    set(hhuvmat.FixedLimits,'Value',0)% put FixedLimits option to 'off'
346    set(hhuvmat.FixedLimits,'BackgroundColor',[0.7 0.7 0.7])
347    uvmat('RootPath_Callback',hObject,eventdata,hhuvmat); %file input with xml reading  in uvmat
348    MenuPlot_Callback(hObject, eventdata, handles)
349    set(handles.ListCoord,'Value',index)% indicate in the list the point with max deviation (possible mistake)
350    ListCoord_Callback(hObject, eventdata, handles)
351    figure(handles.geometry_calib)
352end
353
354%------------------------------------------------------------------
355% --- Executes on button press in calibrate_lin.
356function REPLICATE_Callback(hObject, eventdata, handles)
357%------------------------------------------------------------------------
358calib_cell=get(handles.calib_type,'String');
359val=get(handles.calib_type,'Value');
360calib_type=calib_cell{val};
361Coord_cell=get(handles.ListCoord,'String');
362Object=read_geometry_calib(Coord_cell);
363
364if isequal(calib_type,'rescale')
365    GeometryCalib=calib_rescale(Object.Coord);
366elseif isequal(calib_type,'linear')
367    GeometryCalib=calib_linear(Object.Coord);
368elseif isequal(calib_type,'tsai')
369    GeometryCalib=calib_tsai(Object.Coord);
370end
371% %record image source
372GeometryCalib.SourceCalib.PointCoord=Object.Coord;
373
374%open and read the dataview GUI
375h_dataview=findobj(allchild(0),'name','dataview');
376if ~isempty(h_dataview)
377    delete(h_dataview)
378end
379CalibData=get(handles.geometry_calib,'UserData');%read the calibration image source on the interface userdata
380
381if isfield(CalibData,'XmlInput')
382    XmlInput=fileparts(CalibData.XmlInput);
383    [XmlInput,filename,ext]=fileparts(XmlInput);
384end
385SubCampaignTest='n'; %default
386testinput=0;
387if isfield(CalibData,'Heading')
388    Heading=CalibData.Heading;
389    if isfield(Heading,'Record') && isequal([filename ext],Heading.Record)
390        [XmlInput,filename,ext]=fileparts(XmlInput);
391    end
392    if isfield(Heading,'Device') && isequal([filename ext],Heading.Device)
393        [XmlInput,filename,ext]=fileparts(XmlInput);
394        Device=Heading.Device;
395    end
396    if isfield(Heading,'Experiment') && isequal([filename ext],Heading.Experiment)
397        [PP,filename,ext]=fileparts(XmlInput);
398    end
399    testinput=0;
400    if isfield(Heading,'SubCampaign') && isequal([filename ext],Heading.SubCampaign)
401        SubCampaignTest='y';
402        testinput=1;
403    elseif isfield(Heading,'Campaign') && isequal([filename ext],Heading.Campaign)
404        testinput=1;
405    end
406end
407if ~testinput
408    filename='PROJETS';%default
409    if isfield(CalibData,'XmlInput')
410         [pp,filename]=fileparts(CalibData.XmlInput);
411    end
412    while ~isequal(filename,'PROJETS') && numel(filename)>1
413        filename_1=filename;
414        pp_1=pp;
415        [pp,filename]=fileparts(pp);
416    end
417    XmlInput=fullfile(pp_1,filename_1);
418    testinput=1;
419end
420if testinput
421    outcome=dataview(XmlInput,SubCampaignTest,GeometryCalib);
422end
423
424%------------------------------------------------------------------------
425% determine the parameters for a calibration by an affine function (rescaling and offset, no rotation)
426function GeometryCalib=calib_rescale(Coord)
427%------------------------------------------------------------------------
428 
429X=Coord(:,1);
430Y=Coord(:,2);
431x_ima=Coord(:,4);
432y_ima=Coord(:,5);
433[px,sx]=polyfit(X,x_ima,1);
434[py,sy]=polyfit(Y,y_ima,1);
435T_x=px(2);
436T_y=py(2);
437GeometryCalib.CalibrationType='rescale';
438GeometryCalib.focal=1;
439GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function
440GeometryCalib.Tx_Ty_Tz=[T_x T_y 1];
441GeometryCalib.R=[px(1),0,0;0,py(1),0;0,0,1];
442
443% %check error
444% Calib.dpx=1;
445% Calib.dpy=1;
446% Calib.sx=1;
447% Calib.Cx=0;
448% Calib.Cy=0;
449% Calib.Tz=1;
450% Calib.kappa1=0;
451% Calib.f=GeometryCalib.focal;
452% Calib.Tx=T_x;
453% Calib.Ty=T_y;
454% Calib.R=GeometryCalib.R;
455% [Xpoints,Ypoints]=px_XYZ(Calib,X,Y,0);
456% GeometryCalib.ErrorRms(1)=sqrt(mean((Xpoints-x_ima).*(Xpoints-x_ima)));
457% GeometryCalib.ErrorMax(1)=max(abs(Xpoints-x_ima));
458% GeometryCalib.ErrorRms(2)=sqrt(mean((Ypoints-y_ima).*(Ypoints-y_ima)));
459% GeometryCalib.ErrorMax(2)=max(abs(Ypoints-y_ima));
460
461%------------------------------------------------------------------------
462% determine the parameters for a calibration by a linear transform matrix (rescale and rotation)
463function GeometryCalib=calib_linear(Coord)
464%------------------------------------------------------------------------
465X=Coord(:,1);
466Y=Coord(:,2);
467x_ima=Coord(:,4);
468y_ima=Coord(:,5);
469XY_mat=[ones(size(X)) X Y];
470a_X1=XY_mat\x_ima; %transformation matrix for X
471x1=XY_mat*a_X1;%reconstruction
472err_X1=max(abs(x1-x_ima));%error
473a_Y1=XY_mat\y_ima;%transformation matrix for X
474y1=XY_mat*a_Y1;
475err_Y1=max(abs(y1-y_ima));%error
476T_x=a_X1(1);
477T_y=a_Y1(1);
478GeometryCalib.CalibrationType='linear';
479GeometryCalib.focal=1;
480GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function
481GeometryCalib.Tx_Ty_Tz=[T_x T_y 1];
482GeometryCalib.R=[a_X1(2),a_X1(3),0;a_Y1(2),a_Y1(3),0;0,0,1];
483
484% %check error
485% GeometryCalib.ErrorRms(1)=sqrt(mean((x1-x_ima).*(x1-x_ima)));
486% GeometryCalib.ErrorMax(1)=max(abs(x1-x_ima));
487% GeometryCalib.ErrorRms(2)=sqrt(mean((y1-y_ima).*(y1-y_ima)));
488% GeometryCalib.ErrorMax(2)=max(abs(y1-y_ima));
489
490%------------------------------------------------------------------------
491function GeometryCalib=calib_tsai2(Coord)
492%------------------------------------------------------------------
493path_uvmat=which('uvmat');% check the path detected for source file uvmat
494path_UVMAT=fileparts(path_uvmat); %path to UVMAT
495huvmat=findobj(allchild(0),'Tag','uvmat');
496hhuvmat=guidata(huvmat);
497x_1=Coord(:,4:5)';
498X_1=Coord(:,1:3)';
499n_ima=1;
500% check_cond=0;
501
502
503nx=str2num(get(hhuvmat.npx,'String'));
504ny=str2num(get(hhuvmat.npy,'String'));
505
506
507% est_kc=[1;0;0;0;0];
508est_fc=0;
509est_dist=[1;0;0;0;0];
510run(fullfile(path_UVMAT,'toolbox_calib','go_calib_optim'));
511
512GeometryCalib.CalibrationType='tsai_matlab';
513GeometryCalib.focal=f(2);
514GeometryCalib.dpx_dpy=[1 1];
515GeometryCalib.Cx_Cy=cc';
516GeometryCalib.sx=fc(1)/fc(2);
517GeometryCalib.kappa1=-k(1)/f(2)^2;
518GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function
519GeometryCalib.Tx_Ty_Tz=Tc_1';
520GeometryCalib.R=Rc_1;
521% Calib.dpx=GeometryCalib.dpx_dpy(1);
522% Calib.dpy=GeometryCalib.dpx_dpy(2);
523% Calib.sx=GeometryCalib.sx;
524% Calib.Cx=GeometryCalib.Cx_Cy(1);
525% Calib.Cy=GeometryCalib.Cx_Cy(2);
526% Calib.kappa1=GeometryCalib.kappa1;
527% Calib.f=GeometryCalib.focal;
528% Calib.Tx=GeometryCalib.Tx_Ty_Tz(1);
529% Calib.Ty=GeometryCalib.Tx_Ty_Tz(2);
530% Calib.Tz=GeometryCalib.Tx_Ty_Tz(3);
531% Calib.R=GeometryCalib.R;
532% X=Coord(:,1);
533% Y=Coord(:,2);
534% Z=Coord(:,3);
535% x_ima=Coord(:,4);
536% y_ima=Coord(:,5);
537% [Xpoints,Ypoints]=px_XYZ(Calib,X,Y,Z);
538%
539% GeometryCalib.ErrorRms(1)=sqrt(mean((Xpoints-x_ima).*(Xpoints-x_ima)));
540% [GeometryCalib.ErrorMax(1),GeometryCalib.IndexMax(1)]=max(abs(Xpoints-x_ima));
541% GeometryCalib.ErrorRms(2)=sqrt(mean((Ypoints-y_ima).*(Ypoints-y_ima)));
542% [GeometryCalib.ErrorMax(2),GeometryCalib.IndexMax(2)]=max(abs(Ypoints-y_ima));
543
544function GeometryCalib=calib_tsai(Coord)
545%------------------------------------------------------------------------
546%TSAI
547% 'calibration_lin' provides a linear transform on coordinates,
548path_uvmat=which('uvmat');% check the path detected for source file uvmat
549path_UVMAT=fileparts(path_uvmat); %path to UVMAT
550% if isunix
551    %fid = fopen(fullfile(path_UVMAT,'PARAM_LINUX.txt'),'r');%open the file with civ binary names
552xmlfile=fullfile(path_UVMAT,'PARAM.xml');
553if exist(xmlfile,'file')
554    t=xmltree(xmlfile);
555    sparam=convert(t);
556end
557if ~isfield(sparam,'GeometryCalibBin')
558    msgbox_uvmat('ERROR',['calibration program <GeometryCalibBin> undefined in parameter file ' xmlfile])
559    return
560end
561Tsai_exe=sparam.GeometryCalibBin;
562if ~exist(Tsai_exe,'file')%the binary is defined in /bin, default setting
563     Tsai_exe=fullfile(path_UVMAT,Tsai_exe);
564end
565if ~exist(Tsai_exe,'file')
566    msgbox_uvmat('ERROR',['calibration program ' sparam.GeometryCalibBin ' defined in PARAM.xml does not exist'])
567    return
568end
569
570textcoord=num2str(Coord,4);
571dlmwrite('t.txt',textcoord,''); 
572% ['!' Tsai_exe ' -f1 0 -f2 t.txt']
573    eval(['!' Tsai_exe ' -f t.txt > tsaicalib.log']);
574if ~exist('calib.dat','file')
575    msgbox_uvmat('ERROR','no output from calibration program Tsai_exe: possibly too few points')
576end
577calibdat=dlmread('calib.dat');
578delete('calib.dat')
579delete('t.txt')
580GeometryCalib.CalibrationType='tsai';
581GeometryCalib.focal=calibdat(10);
582GeometryCalib.dpx_dpy=[calibdat(5) calibdat(6)];
583GeometryCalib.Cx_Cy=[calibdat(7) calibdat(8)];
584GeometryCalib.sx=calibdat(9);
585GeometryCalib.kappa1=calibdat(11);
586GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function
587GeometryCalib.Tx_Ty_Tz=[calibdat(12) calibdat(13) calibdat(14)];
588Rx_Ry_Rz=calibdat([15:17]);
589sa = sin(Rx_Ry_Rz(1)) ;
590ca=cos(Rx_Ry_Rz(1));
591sb=sin(Rx_Ry_Rz(2));
592cb =cos(Rx_Ry_Rz(2));
593sg =sin(Rx_Ry_Rz(3));
594cg =cos(Rx_Ry_Rz(3));
595r1 = cb * cg;
596r2 = cg * sa * sb - ca * sg;
597r3 = sa * sg + ca * cg * sb;
598r4 = cb * sg;
599r5 = sa * sb * sg + ca * cg;
600r6 = ca * sb * sg - cg * sa;
601r7 = -sb;
602r8 = cb * sa;
603r9 = ca * cb;
604%EN DEDUIRE MATRICE R ??
605GeometryCalib.R=[r1,r2,r3;r4,r5,r6;r7,r8,r9];
606%erreur a caracteriser?
607% %check error
608% Calib.dpx=GeometryCalib.dpx_dpy(1);
609% Calib.dpy=GeometryCalib.dpx_dpy(2);
610% Calib.sx=GeometryCalib.sx;
611% Calib.Cx=GeometryCalib.Cx_Cy(1);
612% Calib.Cy=GeometryCalib.Cx_Cy(2);
613% Calib.kappa1=GeometryCalib.kappa1;
614% Calib.f=GeometryCalib.focal;
615% Calib.Tx=GeometryCalib.Tx_Ty_Tz(1);
616% Calib.Ty=GeometryCalib.Tx_Ty_Tz(2);
617% Calib.Tz=GeometryCalib.Tx_Ty_Tz(3);
618% Calib.R=GeometryCalib.R;
619% X=Coord(:,1);
620% Y=Coord(:,2);
621% Z=Coord(:,3);
622% x_ima=Coord(:,4);
623% y_ima=Coord(:,5);
624% [Xpoints,Ypoints]=px_XYZ(Calib,X,Y,Z);
625%
626% GeometryCalib.ErrorRms(1)=sqrt(mean((Xpoints-x_ima).*(Xpoints-x_ima)));
627% GeometryCalib.ErrorMax(1)=max(abs(Xpoints-x_ima));
628% GeometryCalib.ErrorRms(2)=sqrt(mean((Ypoints-y_ima).*(Ypoints-y_ima)));
629% GeometryCalib.ErrorMax(2)=max(abs(Ypoints-y_ima));
630
631
632%------------------------------------------------------------------------
633% --- Executes on button press in rotation.
634function rotation_Callback(hObject, eventdata, handles)
635%------------------------------------------------------------------------
636angle_rot=(pi/180)*str2num(get(handles.Phi,'String'));
637Coord_cell=get(handles.ListCoord,'String');
638data=read_geometry_calib(Coord_cell);
639data.Coord(:,1)=cos(angle_rot)*data.Coord(:,1)+sin(angle_rot)*data.Coord(:,2);
640data.Coord(:,1)=-sin(angle_rot)*data.Coord(:,1)+cos(angle_rot)*data.Coord(:,2);
641set(handles.XObject,'String',num2str(data.Coord(:,1),4));
642set(handles.YObject,'String',num2str(data.Coord(:,2),4));
643
644%------------------------------------------------------------------------
645function XImage_Callback(hObject, eventdata, handles)
646%------------------------------------------------------------------------
647update_list(hObject, eventdata,handles)
648
649%------------------------------------------------------------------------
650function YImage_Callback(hObject, eventdata, handles)
651%------------------------------------------------------------------------
652update_list(hObject, eventdata,handles)
653
654function XObject_Callback(hObject, eventdata, handles)
655update_list(hObject, eventdata,handles)
656
657function YObject_Callback(hObject, eventdata, handles)
658update_list(hObject, eventdata,handles)
659
660function ZObject_Callback(hObject, eventdata, handles)
661update_list(hObject, eventdata,handles)
662
663%------------------------------------------------------------------------
664function update_list(hObject, eventdata, handles)
665%------------------------------------------------------------------------
666str4=get(handles.XImage,'String');
667str5=get(handles.YImage,'String');
668str1=get(handles.XObject,'String');
669tt=double(str1);
670str2=get(handles.YObject,'String');
671str3=get(handles.ZObject,'String');
672if ~isempty(str1) & ~isequal(double(str1),32) & (isempty(str3)|isequal(double(str3),32))
673    str3='0';%put z to 0 by default
674end
675strline=[str1 '    |    ' str2 '    |    ' str3 '    |    ' str4 '    |    ' str5];
676Coord=get(handles.ListCoord,'String');
677val=get(handles.ListCoord,'Value');
678Coord{val}=strline;
679set(handles.ListCoord,'String',Coord)
680%update the plot
681ListCoord_Callback(hObject, eventdata, handles)
682MenuPlot_Callback(hObject, eventdata, handles)
683
684%------------------------------------------------------------------------
685% --- Executes on selection change in ListCoord.
686function ListCoord_Callback(hObject, eventdata, handles)
687%------------------------------------------------------------------------
688huvmat=findobj(allchild(0),'Name','uvmat');%find the current uvmat interface handle
689hplot=findobj(huvmat,'Tag','axes3');%main plotting axis of uvmat
690hhh=findobj(hplot,'Tag','calib_marker');
691Coord_cell=get(handles.ListCoord,'String');
692val=get(handles.ListCoord,'Value');
693if numel(val)>1
694    return %no action if several lines have been selected
695end
696coord_str=Coord_cell{val};
697k=findstr('|',coord_str);
698if isempty(k)%last line '.....' selected
699    if ~isempty(hhh)
700        delete(hhh)%delete the circle marker
701    end
702    return
703end
704%fill the edit boxex
705set(handles.XObject,'String',coord_str(1:k(1)-5))
706set(handles.YObject,'String',coord_str(k(1)+5:k(2)-5))
707set(handles.ZObject,'String',coord_str(k(2)+5:k(3)-5))
708set(handles.XImage,'String',coord_str(k(3)+5:k(4)-5))
709set(handles.YImage,'String',coord_str(k(4)+5:end))
710h_menu_coord=findobj(huvmat,'Tag','transform_fct');
711menu=get(h_menu_coord,'String');
712choice=get(h_menu_coord,'Value');
713if iscell(menu)
714    option=menu{choice};
715else
716    option='px'; %default
717end
718if isequal(option,'phys')
719    XCoord=str2num(coord_str(1:k(1)-5));
720    YCoord=str2num(coord_str(k(1)+5:k(2)-5));
721elseif isequal(option,'px')|| isequal(option,'')
722    XCoord=str2num(coord_str(k(3)+5:k(4)-5));
723    YCoord=str2num(coord_str(k(4)+5:end));
724else
725    msgbox_uvmat('ERROR','the choice in menu_coord of uvmat must be px or phys ')
726end
727if isempty(XCoord)||isempty(YCoord)
728     if ~isempty(hhh)
729        delete(hhh)%delete the circle marker
730    end
731    return
732end
733xlim=get(hplot,'XLim');
734ylim=get(hplot,'YLim');
735ind_range=max(abs(xlim(2)-xlim(1)),abs(ylim(end)-ylim(1)))/20;%defines the size of the circle marker
736if isempty(hhh)
737    axes(hplot)
738    rectangle('Curvature',[1 1],...
739              'Position',[XCoord-ind_range/2 YCoord-ind_range/2 ind_range ind_range],'EdgeColor','m',...
740              'LineStyle','-','Tag','calib_marker');
741else
742    set(hhh,'Position',[XCoord-ind_range/2 YCoord-ind_range/2 ind_range ind_range])
743end
744
745%------------------------------------------------------------------------
746% --- Executes on selection change in edit_append.
747function edit_append_Callback(hObject, eventdata, handles)
748%------------------------------------------------------------------------
749choice=get(handles.edit_append,'Value');
750% if choice==1
751%        Coord=get(handles.ListCoord,'String');
752%        val=length(Coord);
753%        if val>=1 & isequal(Coord{val},'')
754%             val=val-1; %do not take into account blank
755%        end
756%        Coord{val+1}='';
757%        set(handles.ListCoord,'String',Coord)
758%        set(handles.ListCoord,'Value',val+1)
759% end
760if choice
761    set(handles.edit_append,'BackgroundColor',[1 1 0])
762else
763    set(handles.edit_append,'BackgroundColor',[0.7 0.7 0.7])
764end
765   
766function NEW_Callback(hObject, eventdata, handles)
767%A METTRE SOUS UN BOUTON
768huvmat=findobj(allchild(0),'Name','uvmat');
769hchild=get(huvmat,'children');
770hcoord=findobj(hchild,'Tag','menu_coord');
771coordtype=get(hcoord,'Value');
772haxes=findobj(hchild,'Tag','axes3');
773AxeData=get(haxes,'UserData');
774if ~isequal(hcoord,2)
775    set(hcoord,'Value',2)
776    huvmat=uvmat(AxeData);
777    'relancer uvmat';
778end
779if ~isfield(AxeData,'ZoomAxes')
780    msgbox_uvmat('ERROR','first draw a window around a grid marker')
781    return
782end
783XLim=get(AxeData.ZoomAxes,'XLim');
784YLim=get(AxeData.ZoomAxes,'YLim');
785np=size(AxeData.A);
786ind_sub_x=round(XLim);
787ind_sub_y=np(1)-round(YLim);
788Mfiltre=AxeData.A([ind_sub_y(2):ind_sub_y(1)] ,ind_sub_x,:);
789Mfiltre_norm=double(Mfiltre);
790Mfiltre_norm=Mfiltre_norm/sum(sum(Mfiltre_norm));
791Mfiltre_norm=100*(Mfiltre_norm-mean(mean(Mfiltre_norm)));
792Atype=class(AxeData.A);
793Data.NbDim=2;
794Data.A=filter2(Mfiltre_norm,double(AxeData.A));
795Data.A=feval(Atype,Data.A);
796Data.AName='image';
797Data.AX=AxeData.AX;
798Data.AY=AxeData.AY;
799Data.CoordType='px';
800plot_field(Data)
801
802%------------------------------------------------------------------------
803% --- 'key_press_fcn:' function activated when a key is pressed on the keyboard
804function key_press_fcn(hObject,eventdata,handles)
805%------------------------------------------------------------------------
806xx=double(get(handles.geometry_calib,'CurrentCharacter')); %get the keyboard character
807if ismember(xx,[8 127])%backspace or delete
808    Coord_cell=get(handles.ListCoord,'String');
809    val=get(handles.ListCoord,'Value');
810     if max(val)<numel(Coord_cell) % the last element '...' has not been selected
811        Coord_cell(val)=[];%remove the selected line
812        set(handles.ListCoord,'Value',min(val))
813        set(handles.ListCoord,'String',Coord_cell)         
814        ListCoord_Callback(hObject, eventdata, handles)
815        MenuPlot_Callback(hObject,eventdata,handles)
816     end
817end
818
819% %------------------------------------------------------------------------
820% % --- Executes on button press in append_point.
821% function append_point_Callback(hObject, eventdata, handles)
822% %------------------------------------------------------------------------
823%        Coord=get(handles.ListCoord,'String');
824%        val=length(Coord);
825%        if val>=1 & isequal(Coord{val},'')
826%             val=val-1; %do not take into account blank
827%        end
828%        Coord{val+1}='';
829%        set(handles.ListCoord,'String',Coord)
830%        set(handles.ListCoord,'Value',val+1)
831
832%------------------------------------------------------------------------
833function MenuOpen_Callback(hObject, eventdata, handles)
834%------------------------------------------------------------------------
835%get the object file
836huvmat=findobj(allchild(0),'Name','uvmat');
837UvData=get(huvmat,'UserData');
838hchild=get(huvmat,'Children');
839hrootpath=findobj(hchild,'Tag','RootPath');
840oldfile=get(hrootpath,'String');
841if isempty(oldfile)
842    oldfile='';
843end
844%[FileName,PathName] = uigetfile('*.civ','Select a .civ file',oldfile)
845[FileName, PathName, filterindex] = uigetfile( ...
846       {'*.xml;*.mat', ' (*.xml,*.mat)';
847       '*.xml',  '.xml files '; ...
848        '*.mat',  '.mat matlab files '}, ...
849        'Pick a file',oldfile);
850fileinput=[PathName FileName];%complete file name
851testblank=findstr(fileinput,' ');%look for blanks
852if ~isempty(testblank)
853    msgbox_uvmat('ERROR','forbidden input file name or path: no blank character allowed')
854    return
855end
856sizf=size(fileinput);
857if (~ischar(fileinput)|~isequal(sizf(1),1)),return;end
858loadfile(handles,fileinput)
859
860
861%------------------------------------------------------------------------
862function MenuPlot_Callback(hObject, eventdata, handles)
863%------------------------------------------------------------------------
864huvmat=findobj(allchild(0),'Name','uvmat');%find the current uvmat interface handle
865UvData=get(huvmat,'UserData');%Data associated to the current uvmat interface
866hhuvmat=guidata(huvmat); %handles of GUI elements in uvmat
867hplot=findobj(huvmat,'Tag','axes3');%main plotting axis of uvmat
868h_menu_coord=findobj(huvmat,'Tag','transform_fct');
869menu=get(h_menu_coord,'String');
870choice=get(h_menu_coord,'Value');
871if iscell(menu)
872    option=menu{choice};
873else
874    option='px'; %default
875end
876Coord_cell=get(handles.ListCoord,'String');
877ObjectData=read_geometry_calib(Coord_cell);
878%ObjectData=read_geometry_calib(handles);%read the interface input parameters defining the object
879if ~isempty(ObjectData.Coord)
880    if isequal(option,'phys')
881        ObjectData.Coord=ObjectData.Coord(:,[1:3]);
882    elseif isequal(option,'px')||isequal(option,'')
883        ObjectData.Coord=ObjectData.Coord(:,[4:5]);
884    else
885        msgbox_uvmat('ERROR','the choice in menu_coord of uvmat must be '''', px or phys ')
886    end
887end
888axes(hhuvmat.axes3)
889hh=findobj('Tag','calib_points');
890if  ~isempty(ObjectData.Coord) && isempty(hh)
891    hh=line(ObjectData.Coord(:,1),ObjectData.Coord(:,2),'Color','m','Tag','calib_points','LineStyle','.','Marker','+');
892elseif isempty(ObjectData.Coord)%empty list of points, suppress the plot
893    delete(hh)
894else
895    set(hh,'XData',ObjectData.Coord(:,1))
896    set(hh,'YData',ObjectData.Coord(:,2))
897end
898pause(.1)
899figure(handles.geometry_calib)
900
901% --------------------------------------------------------------------
902function MenuHelp_Callback(hObject, eventdata, handles)
903path_to_uvmat=which ('uvmat');% check the path of uvmat
904pathelp=fileparts(path_to_uvmat);
905    helpfile=fullfile(pathelp,'uvmat_doc','uvmat_doc.html');
906if isempty(dir(helpfile)), msgbox_uvmat('ERROR','Please put the help file uvmat_doc.html in the sub-directory /uvmat_doc of the UVMAT package')
907else
908   addpath (fullfile(pathelp,'uvmat_doc'))
909   web([helpfile '#geometry_calib'])
910end
911
912%------------------------------------------------------------------------
913function MenuCreateGrid_Callback(hObject, eventdata, handles)
914%------------------------------------------------------------------------
915%hcalib=get(handles.calib_type,'parent');%handles of the GUI geometry_calib
916CalibData=get(handles.geometry_calib,'UserData');
917Tinput=[];%default
918if isfield(CalibData,'grid')
919    Tinput=CalibData.grid;
920end
921[T,CalibData.grid]=create_grid(Tinput);%display the GUI create_grid
922set(handles.geometry_calib,'UserData',CalibData)
923
924%grid in phys space
925Coord=get(handles.ListCoord,'String');
926val=get(handles.ListCoord,'Value');
927data=read_geometry_calib(Coord);
928%nbpoints=size(data.Coord,1); %nbre of calibration points
929data.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
930% for i=1:nbpoints
931%    for j=1:5
932%           Coord{i,j}=num2str(data.Coord(i,j),4);%display coordiantes with 4 digits
933%    end
934% end
935%update the phys coordinates starting from the selected point (down in the
936Coord(end,:)=[]; %remove last string '.....'
937for i=1:size(data.Coord,1)
938    for j=1:5
939          Coord{i,j}=num2str(data.Coord(i,j),4);%display coordiantes with 4 digits
940    end
941end
942
943%size(data.Coord,1)
944Tabchar=cell2tab(Coord,'    |    ');
945Tabchar=[Tabchar ;{'......'}];
946set(handles.ListCoord,'String',Tabchar)
947
948% -----------------------------------------------------------------------
949% --- automatic grid dectection from local maxima of the images
950function MenuDetectGrid_Callback(hObject, eventdata, handles)
951%------------------------------------------------------------------------
952CalibData=get(handles.geometry_calib,'UserData');%get information stored on the GUI geometry_calib
953grid_input=[];%default
954if isfield(CalibData,'grid')
955    grid_input=CalibData.grid;%retrieve the previously used grid
956end
957[T,CalibData.grid]=create_grid(grid_input,'detect grid');%display the GUI create_grid, read the set of phys coordinates T
958set(handles.geometry_calib,'UserData',CalibData)%store the phys grid parameters for later use
959
960%read the four last point coordinates in pixels
961Coord_cell=get(handles.ListCoord,'String');%read list of coordinates on geometry_calib
962data=read_geometry_calib(Coord_cell);
963nbpoints=size(data.Coord,1); %nbre of calibration points
964if nbpoints~=4
965    msgbox_uvmat('ERROR','four points must be selected by the mouse, beginning by the new x axis, to delimitate the phs grid area')
966    return
967end
968corners_X=(data.Coord(end:-1:end-3,4)); %pixel absissa of the four corners
969corners_Y=(data.Coord(end:-1:end-3,5));
970
971%reorder the last two points (the two first in the list) if needed
972angles=angle((corners_X-corners_X(1))+i*(corners_Y-corners_Y(1)));
973if abs(angles(4)-angles(2))>abs(angles(3)-angles(2))
974      X_end=corners_X(4);
975      Y_end=corners_Y(4);
976      corners_X(4)=corners_X(3);
977      corners_Y(4)=corners_Y(3);
978      corners_X(3)=X_end;
979      corners_Y(3)=Y_end;
980end
981
982%read the current image, displayed in the GUI uvmat
983huvmat=findobj(allchild(0),'Name','uvmat');
984UvData=get(huvmat,'UserData');
985A=UvData.Field.A;
986npxy=size(A);
987%linear transform on the current image
988X=[CalibData.grid.x_0 CalibData.grid.x_1 CalibData.grid.x_0 CalibData.grid.x_1]';%corner absissa in the phys coordinates
989Y=[CalibData.grid.y_0 CalibData.grid.y_0 CalibData.grid.y_1 CalibData.grid.y_1]';%corner ordinates in the phys coordinates
990%X=[CalibData.grid.x_0 CalibData.grid.x_0 CalibData.grid.x_1 CalibData.grid.x_1]';%corner absissa in the phys coordinates
991%Y=[CalibData.grid.y_0 CalibData.grid.y_1 CalibData.grid.y_1 CalibData.grid.y_0]';%corner ordinates in the phys coordinate
992%XY_mat=[ones(size(X)) X Y];
993%a_X1=XY_mat\corners_X; %transformation matrix for X
994%x1=XY_mat*a_X1;%reconstruction
995%err_X1=max(abs(x1-corners_X));%error
996%a_Y1=XY_mat\corners_Y;%transformation matrix for X
997%y1=XY_mat*a_Y1;
998%err_Y1=max(abs(y1-corners_Y));%error
999% figure(1)
1000% hold off
1001% plot([corners_X;corners_X(1)],[corners_Y;corners_Y(1)],'r')
1002%test points along contour
1003% u=[0:0.1:1];
1004% x12=u*corners_X(2)+(1-u)*corners_X(1);
1005% x23=u*corners_X(3)+(1-u)*corners_X(2);
1006% x34=u*corners_X(4)+(1-u)*corners_X(3);
1007% x41=u*corners_X(1)+(1-u)*corners_X(4);
1008% y12=u*corners_Y(2)+(1-u)*corners_Y(1);
1009% y23=u*corners_Y(3)+(1-u)*corners_Y(2);
1010% y34=u*corners_Y(4)+(1-u)*corners_Y(3);
1011% y41=u*corners_Y(1)+(1-u)*corners_Y(4);
1012% hold on
1013% plot(x12,y12,'xr')
1014% plot(x23,y23,'xr')
1015% plot(x34,y34,'xr')
1016% plot(x41,y41,'xr')
1017%calculate transform matrices:
1018% reference: http://alumni.media.mit.edu/~cwren/interpolator/ by Christopher R. Wren
1019B = [ X Y ones(size(X)) zeros(4,3)        -X.*corners_X -Y.*corners_X ...
1020      zeros(4,3)        X Y ones(size(X)) -X.*corners_Y -Y.*corners_Y ];
1021B = reshape (B', 8 , 8 )';
1022D = [ corners_X , corners_Y ];
1023D = reshape (D', 8 , 1 );
1024l = inv(B' * B) * B' * D;
1025Amat = reshape([l(1:6)' 0 0 1 ],3,3)';
1026C = [l(7:8)' 1];
1027
1028GeometryCalib.CalibrationType='tsai';%'linear';
1029GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function
1030GeometryCalib.f=1;
1031GeometryCalib.dpx=1;
1032GeometryCalib.dpy=1;
1033GeometryCalib.sx=1;
1034GeometryCalib.Cx=0;
1035GeometryCalib.Cy=0;
1036GeometryCalib.kappa1=0;
1037GeometryCalib.Tx=Amat(1,3);
1038GeometryCalib.Ty=Amat(2,3);
1039GeometryCalib.Tz=1;
1040GeometryCalib.R=[Amat(1,1),Amat(1,2),0;Amat(2,1),Amat(2,2),0;C(1),C(2),0];
1041%montre points du pourtour transformes
1042% for ip=1:numel(u),
1043%   XY12(:,ip)=Amat*[x12(ip);y12(ip);1]/(C*[x12(ip);y12(ip);1]);
1044%   XY23(:,ip)=Amat*[x23(ip);y23(ip);1]/(C*[x23(ip);y23(ip);1]);
1045%   XY34(:,ip)=Amat*[x34(ip);y34(ip);1]/(C*[x34(ip);y34(ip);1]);
1046%   XY41(:,ip)=Amat*[x41(ip);y41(ip);1]/(C*[x41(ip);y41(ip);1]);
1047% end
1048% xphys12=u*X(2)+(1-u)*X(1);
1049% xphys23=u*X(3)+(1-u)*X(2);
1050% xphys34=u*X(4)+(1-u)*X(3);
1051% xphys41=u*X(1)+(1-u)*X(4);
1052% yphys12=u*Y(2)+(1-u)*Y(1);
1053% yphys23=u*Y(3)+(1-u)*Y(2);
1054% yphys34=u*Y(4)+(1-u)*Y(3);
1055% yphys41=u*Y(1)+(1-u)*Y(4);
1056% for ip=1:numel(u),
1057%   XY12(:,ip)=Amat*[xphys12(ip);yphys12(ip);1]/(C*[xphys12(ip);yphys12(ip);1]);
1058%   XY23(:,ip)=Amat*[xphys23(ip);yphys23(ip);1]/(C*[xphys23(ip);yphys23(ip);1]);
1059%   XY34(:,ip)=Amat*[xphys34(ip);yphys34(ip);1]/(C*[xphys34(ip);yphys34(ip);1]);
1060%   XY41(:,ip)=Amat*[xphys41(ip);yphys41(ip);1]/(C*[xphys41(ip);yphys41(ip);1]);
1061% end
1062% x12=XY12(1,:);
1063% y12=XY12(2,:);
1064% x23=XY23(1,:);
1065% y23=XY23(2,:);
1066% x34=XY34(1,:);
1067% y34=XY34(2,:);
1068% x41=XY41(1,:);
1069% y41=XY41(2,:);
1070% [x12,y12]=px_XYZ(GeometryCalib,xphys12,yphys12);
1071% [x23,y23]=px_XYZ(GeometryCalib,xphys23,yphys23);
1072% [x34,y34]=px_XYZ(GeometryCalib,xphys34,yphys34);
1073% [x41,y41]=px_XYZ(GeometryCalib,xphys41,yphys41);
1074% plot(x12,y12,'xb')
1075% plot(x23,y23,'xb')
1076% plot(x34,y34,'xb')
1077% plot(x41,y41,'xb')
1078% figure(2)
1079% hold off
1080% plot(XY12(1,:),XY12(2,:),'ob') 
1081% hold on
1082% plot(XY23(1,:),XY23(2,:),'ob')
1083% plot(XY34(1,:),XY34(2,:),'ob')
1084% plot(XY41(1,:),XY41(2,:),'ob')
1085% plot(xphys12,yphys12,'ob') 
1086% hold on
1087% plot(xphys23,yphys23,'ob')
1088% plot(xphys34,yphys34,'ob')
1089% plot(xphys41,yphys41,'ob')
1090
1091[Amod,Rangx,Rangy]=phys_Ima(A-min(min(A)),GeometryCalib,0);
1092
1093Amod=double(Amod);
1094figure(12)
1095Amax=max(max(Amod));
1096image(Rangx,Rangy,uint8(255*Amod/Amax))
1097Dx=(Rangx(2)-Rangx(1))/(npxy(2)-1); %x mesh in real space
1098Dy=(Rangy(2)-Rangy(1))/(npxy(1)-1); %y mesh in real space
1099ind_range_x=ceil(GeometryCalib.R(1,1)*CalibData.grid.Dx/3);% range of search of image ma around each point obtained by linear interpolation from the marked points
1100ind_range_y=ceil(GeometryCalib.R(2,2)*CalibData.grid.Dy/3);% range of search of image ma around each point obtained by linear interpolation from the marked points
1101nbpoints=size(T,1);
1102for ipoint=1:nbpoints
1103    i0=1+round((T(ipoint,1)-Rangx(1))/Dx);%round(Xpx(ipoint));
1104    j0=1+round((T(ipoint,2)-Rangy(1))/Dy);%round(Xpx(ipoint));
1105    Asub=Amod(j0-ind_range_y:j0+ind_range_y,i0-ind_range_x:i0+ind_range_x);
1106    x_profile=sum(Asub,1);
1107    y_profile=sum(Asub,2);
1108    [Amax,ind_x_max]=max(x_profile);
1109    [Amax,ind_y_max]=max(y_profile);
1110    %sub-pixel improvement using moments
1111    x_shift=0;
1112    y_shift=0;
1113    if ind_x_max+2<=2*ind_range_x+1 && ind_x_max-2>=1
1114        Atop=x_profile(ind_x_max-2:ind_x_max+2);
1115        x_shift=sum(Atop.*[-2 -1 0 1 2])/sum(Atop);
1116    end
1117    if ind_y_max+2<=2*ind_range_y+1 && ind_y_max-2>=1
1118        Atop=y_profile(ind_y_max-2:ind_y_max+2);
1119        y_shift=sum(Atop.*[-2 -1 0 1 2]')/sum(Atop);
1120    end
1121    Delta(ipoint,1)=(x_shift+ind_x_max-ind_range_x-1)*Dx;%shift from the initial guess
1122    Delta(ipoint,2)=(y_shift+ind_y_max-ind_range_y-1)*Dy;
1123end
1124Tmod=T(:,(1:2))+Delta;
1125[Xpx,Ypx]=px_XYZ(GeometryCalib,Tmod(:,1),Tmod(:,2));
1126for ipoint=1:nbpoints
1127     Coord{ipoint,1}=num2str(T(ipoint,1),4);%display coordiantes with 4 digits
1128     Coord{ipoint,2}=num2str(T(ipoint,2),4);%display coordiantes with 4 digits
1129     Coord{ipoint,3}='0';
1130     Coord{ipoint,4}=num2str(Xpx(ipoint),4);%display coordiantes with 4 digi
1131     Coord{ipoint,5}=num2str(Ypx(ipoint),4);%display coordiantes with 4 digi
1132end
1133Tabchar=cell2tab(Coord(end:-1:1,:),'    |    ');
1134Tabchar=[Tabchar ;{'......'}];
1135set(handles.ListCoord,'Value',1)
1136set(handles.ListCoord,'String',Tabchar)
1137MenuPlot_Callback(hObject, eventdata, handles)
1138
1139%-----------------------------------------------------------------------
1140function MenuTranslatePoints_Callback(hObject, eventdata, handles)
1141%-----------------------------------------------------------------------
1142%hcalib=get(handles.calib_type,'parent');%handles of the GUI geometry_calib
1143CalibData=get(handles.geometry_calib,'UserData');
1144Tinput=[];%default
1145if isfield(CalibData,'translate')
1146    Tinput=CalibData.translate;
1147end
1148T=translate_points(Tinput);%display translate_points GUI and get shift parameters
1149CalibData.translate=T;
1150set(handles.geometry_calib,'UserData',CalibData)
1151%translation
1152Coord_cell=get(handles.ListCoord,'String');
1153data=read_geometry_calib(Coord_cell);
1154data.Coord(:,1)=T(1)+data.Coord(:,1);
1155data.Coord(:,2)=T(2)+data.Coord(:,2);
1156data.Coord(:,3)=T(3)+data.Coord(:,3);
1157data.Coord(:,[4 5])=data.Coord(:,[4 5]);
1158for i=1:size(data.Coord,1)
1159    for j=1:5
1160          Coord{i,j}=num2str(data.Coord(i,j),4);%phys x,y,z
1161   end
1162end
1163Tabchar=cell2tab(Coord,'    |    ');
1164Tabchar=[Tabchar; {'.....'}];
1165%set(handles.ListCoord,'Value',1)
1166set(handles.ListCoord,'String',Tabchar)
1167
1168
1169% --------------------------------------------------------------------
1170function MenuRotatePoints_Callback(hObject, eventdata, handles)
1171%hcalib=get(handles.calib_type,'parent');%handles of the GUI geometry_calib
1172CalibData=get(handles.geometry_calib,'UserData');
1173Tinput=[];%default
1174if isfield(CalibData,'rotate')
1175    Tinput=CalibData.rotate;
1176end
1177T=rotate_points(Tinput);%display translate_points GUI and get shift parameters
1178CalibData.rotate=T;
1179set(handles.geometry_calib,'UserData',CalibData)
1180%-----------------------------------------------------
1181%rotation
1182Phi=T(1);
1183O_x=0;%default
1184O_y=0;%default
1185if numel(T)>=2
1186    O_x=T(2);%default
1187end
1188if numel(T)>=3
1189    O_y=T(3);%default
1190end
1191Coord_cell=get(handles.ListCoord,'String');
1192data=read_geometry_calib(Coord_cell);
1193r1=cos(pi*Phi/180);
1194r2=-sin(pi*Phi/180);
1195r3=sin(pi*Phi/180);
1196r4=cos(pi*Phi/180);
1197x=data.Coord(:,1)-O_x;
1198y=data.Coord(:,2)-O_y;
1199data.Coord(:,1)=r1*x+r2*y;
1200data.Coord(:,2)=r3*x+r4*y;
1201% data.Coord(:,[4 5])=data.Coord(:,[4 5]);
1202for i=1:size(data.Coord,1)
1203    for j=1:5
1204          Coord{i,j}=num2str(data.Coord(i,j),4);%phys x,y,z
1205   end
1206end
1207Tabchar=cell2tab(Coord,'    |    ');
1208Tabchar=[Tabchar;{'......'}];
1209set(handles.ListCoord,'Value',1)
1210set(handles.ListCoord,'String',Tabchar)
1211
1212%------------------------------------------------------------------------
1213% image transform from px to phys
1214%INPUT:
1215%Zindex: index of plane
1216function [A_out,Rangx,Rangy]=phys_Ima(A,Calib,ZIndex)
1217%------------------------------------------------------------------------
1218xcorner=[];
1219ycorner=[];
1220npx=[];
1221npy=[];
1222siz=size(A)
1223npx=[npx siz(2)];
1224npy=[npy siz(1)]
1225xima=[0.5 siz(2)-0.5 0.5 siz(2)-0.5];%image coordinates of corners
1226yima=[0.5 0.5 siz(1)-0.5 siz(1)-0.5];
1227[xcorner,ycorner]=phys_XYZ(Calib,xima,yima,ZIndex);%corresponding physical coordinates
1228Rangx(1)=min(xcorner);
1229Rangx(2)=max(xcorner);
1230Rangy(2)=min(ycorner);
1231Rangy(1)=max(ycorner);
1232test_multi=(max(npx)~=min(npx)) | (max(npy)~=min(npy));
1233npx=max(npx);
1234npy=max(npy);
1235x=linspace(Rangx(1),Rangx(2),npx);
1236y=linspace(Rangy(1),Rangy(2),npy);
1237[X,Y]=meshgrid(x,y);%grid in physical coordiantes
1238vec_B=[];
1239
1240zphys=0; %default
1241if isfield(Calib,'SliceCoord') %.Z= index of plane
1242   SliceCoord=Calib.SliceCoord(ZIndex,:);
1243   zphys=SliceCoord(3); %to generalize for non-parallel planes
1244end
1245[XIMA,YIMA]=px_XYZ(Calib,X,Y,zphys);%corresponding image indices for each point in the real space grid
1246XIMA=reshape(round(XIMA),1,npx*npy);%indices reorganized in 'line'
1247YIMA=reshape(round(YIMA),1,npx*npy);
1248flagin=XIMA>=1 & XIMA<=npx & YIMA >=1 & YIMA<=npy;%flagin=1 inside the original image
1249testuint8=isa(A,'uint8');
1250testuint16=isa(A,'uint16');
1251if numel(siz)==2 %(B/W images)
1252    vec_A=reshape(A,1,npx*npy);%put the original image in line
1253    ind_in=find(flagin);
1254    ind_out=find(~flagin);
1255    ICOMB=((XIMA-1)*npy+(npy+1-YIMA));
1256    ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A
1257    vec_B(ind_in)=vec_A(ICOMB);
1258    vec_B(ind_out)=zeros(size(ind_out));
1259    A_out=reshape(vec_B,npy,npx);%new image in real coordinates
1260elseif numel(siz)==3     
1261    for icolor=1:siz(3)
1262        vec_A=reshape(A{icell}(:,:,icolor),1,npx*npy);%put the original image in line
1263        ind_in=find(flagin);
1264        ind_out=find(~flagin);
1265        ICOMB=((XIMA-1)*npy+(npy+1-YIMA));
1266        ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A
1267        vec_B(ind_in)=vec_A(ICOMB);
1268        vec_B(ind_out)=zeros(size(ind_out));
1269        A_out(:,:,icolor)=reshape(vec_B,npy,npx);%new image in real coordinates
1270    end
1271end
1272if testuint8
1273    A_out=uint8(A_out);
1274end
1275if testuint16
1276    A_out=uint16(A_out);
1277end
1278
1279%------------------------------------------------------------------------
1280% pointwise transform from px to phys
1281%INPUT:
1282%Z: index of plane
1283function [Xphys,Yphys,Zphys]=phys_XYZ(Calib,X,Y,Z)
1284%------------------------------------------------------------------------
1285if exist('Z','var')& isequal(Z,round(Z))& Z>0 & isfield(Calib,'SliceCoord')&length(Calib.SliceCoord)>=Z
1286    Zindex=Z;
1287    Zphys=Calib.SliceCoord(Zindex,3);%GENERALISER AUX CAS AVEC ANGLE
1288else
1289%     if exist('Z','var')
1290%         Zphys=Z;
1291%     else
1292        Zphys=0;
1293%     end
1294end
1295if ~exist('X','var')||~exist('Y','var')
1296    Xphys=[];
1297    Yphys=[];%default
1298    return
1299end
1300Xphys=X;%default
1301Yphys=Y;
1302%image transform
1303if isfield(Calib,'R')
1304    R=(Calib.R)';
1305    Dx=R(5)*R(7)-R(4)*R(8);
1306    Dy=R(1)*R(8)-R(2)*R(7);
1307    D0=Calib.f*(R(2)*R(4)-R(1)*R(5));
1308    Z11=R(6)*R(8)-R(5)*R(9);
1309    Z12=R(2)*R(9)-R(3)*R(8); 
1310    Z21=R(4)*R(9)-R(6)*R(7);
1311    Z22=R(3)*R(7)-R(1)*R(9);
1312    Zx0=R(3)*R(5)-R(2)*R(6);
1313    Zy0=R(1)*R(6)-R(3)*R(4);
1314    A11=R(8)*Calib.Ty-R(5)*Calib.Tz+Z11*Zphys;
1315    A12=R(2)*Calib.Tz-R(8)*Calib.Tx+Z12*Zphys;
1316    A21=-R(7)*Calib.Ty+R(4)*Calib.Tz+Z21*Zphys;
1317    A22=-R(1)*Calib.Tz+R(7)*Calib.Tx+Z11*Zphys;
1318    X0=Calib.f*(R(5)*Calib.Tx-R(2)*Calib.Ty+Zx0*Zphys);
1319    Y0=Calib.f*(-R(4)*Calib.Tx+R(1)*Calib.Ty+Zy0*Zphys);
1320        %px to camera:
1321    Xd=(Calib.dpx/Calib.sx)*(X-Calib.Cx); % sensor coordinates
1322    Yd=Calib.dpy*(Y-Calib.Cy);
1323    dist_fact=1+Calib.kappa1*(Xd.*Xd+Yd.*Yd); %distortion factor
1324    Xu=dist_fact.*Xd;%undistorted sensor coordinates
1325    Yu=dist_fact.*Yd;
1326    denom=Dx*Xu+Dy*Yu+D0;
1327    % denom2=denom.*denom;
1328    Xphys=(A11.*Xu+A12.*Yu+X0)./denom;%world coordinates
1329    Yphys=(A21.*Xu+A22.*Yu+Y0)./denom;
1330end
Note: See TracBrowser for help on using the repository browser.