%'geometry_calib': associated to the GUI geometry_calib to perform geometric calibration from a set of reference points %------------------------------------------------------------------------ % function hgeometry_calib = geometry_calib(inputfile,pos) % %OUTPUT: % hgeometry_calib=current handles of the GUI geometry_calib.fig % %INPUT: % inputfile: (optional) name of an xml file containing coordinates of reference points % pos: (optional) 4 element vector setting the 'Position' of the GUI %======================================================================= % Copyright 2008-2018, LEGI UMR 5519 / CNRS UGA G-INP, Grenoble, France % http://www.legi.grenoble-inp.fr % Joel.Sommeria - Joel.Sommeria (A) legi.cnrs.fr % % This file is part of the toolbox UVMAT. % % UVMAT is free software; you can redistribute it and/or modify % it under the terms of the GNU General Public License as published % by the Free Software Foundation; either version 2 of the license, % or (at your option) any later version. % % UVMAT is distributed in the hope that it will be useful, % but WITHOUT ANY WARRANTY; without even the implied warranty of % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the % GNU General Public License (see LICENSE.txt) for more details. %======================================================================= function varargout = geometry_calib(varargin) % GEOMETRY_CALIB M-file for geometry_calib.fig % GEOMETRY_CALIB, by itself, creates a MenuCoord GEOMETRY_CALIB or raises the existing % singleton*. % % H = GEOMETRY_CALIB returns the handle to a MenuCoord GEOMETRY_CALIB or the handle to % the existing singleton*. % % GEOMETRY_CALIB('CALLBACK',hObject,eventData,handles,...) calls the local % function named CALLBACK in GEOMETRY_CALIB.M with the given input arguments. % % GEOMETRY_CALIB('Property','Value',...) creates a MenuCoord GEOMETRY_CALIB or raises the % existing singleton*. Starting from the left, property value pairs are % applied to the GUI before geometry_calib_OpeningFunction gets called. An % unrecognized property name or invalid value makes property application % stop. All inputs are passed to geometry_calib_OpeningFcn via varargin. % % *See GUI Options on GUIDE's Tools menu. Choose "GUI allows only one % instance to run (singleton)". % % See also: GUIDE, GUIDATA, GUIHANDLES % Edit the above text to modify the response to help geometry_calib % Last Modified by GUIDE v2.5 16-Apr-2015 17:29:02 % Begin initialization code - DO NOT edit gui_Singleton = 1; gui_State = struct('gui_Name', mfilename, ... 'gui_Singleton', gui_Singleton, ... 'gui_OpeningFcn', @geometry_calib_OpeningFcn, ... 'gui_OutputFcn', @geometry_calib_OutputFcn, ... 'gui_LayoutFcn', [] , ... 'gui_Callback', []); if nargin [pp,ff]=fileparts(which(varargin{1})); % name of the input file if strcmp(ff,mfilename)% if we are activating a sub-function of geometry_calib % ~isempty(regexp(varargin{1},'_Callback','once')) gui_State.gui_Callback = str2func(varargin{1}); end end if nargout [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:}); else gui_mainfcn(gui_State, varargin{:}); end % End initialization code - DO NOT edit % --- Executes just before geometry_calib is made visible. %INPUT: %handles: handles of the geometry_calib interface elements % PlotHandles: set of handles of the elements contolling the plotting % parameters on the uvmat interface (obtained by 'get_plot_handle.m') %------------------------------------------------------------------------ function geometry_calib_OpeningFcn(hObject, eventdata, handles,inputfile) %------------------------------------------------------------------------ % Choose default command line output for geometry_calib handles.output = hObject; % Update handles structure guidata(hObject, handles); set(hObject,'DeleteFcn',{@closefcn})% %set(hObject,'WindowButtonDownFcn',{'mouse_alt_gui',handles}) % allows mouse action with right button (zoom for uicontrol display) %% position set(0,'Unit','pixels') ScreenSize=get(0,'ScreenSize');% get the size of the screen, to put the fig on the upper right Left=ScreenSize(3)- 460; %right edge close to the right, with margin=40 (GUI width=420 px) if ScreenSize(4)>920 Height=840;%default height of the GUI Bottom=ScreenSize(4)-Height-40; %put fig at top right else Height=ScreenSize(4)-80; Bottom=40; % GUI lies o the screen bottom (with margin =40) end set(handles.calib_type,'Position',[1 Height-40 194 30])% rank 1 set(handles.APPLY,'Position',[197 Height-40 110 30])% rank 1 set(handles.REPLICATE,'Position',[309 Height-40 110 30])% rank 1 set(handles.Intrinsic,'Position',[1 Height-40-2-92 418 92])% rank 2 set(handles.Extrinsic,'Position',[1 Height-40-4-92-75 418 75])% rank 3 set(handles.PointLists,'Position',[1 Height-40-6-92-75-117 418 117]) % rank 4 set(handles.CheckEnableMouse,'Position',[3 Height-362 208 30])% rank 5 set(handles.PLOT,'Position',[3 Height-394 120 30])% rank 6 set(handles.Copy,'Position',[151 Height-394 120 30])% rank 6 set(handles.ClearAll,'Position',[297 Height-394 120 30])% rank 6 set(handles.ClearPoint,'Position',[297 Height-362 120 30])% rank 6 set(handles.CoordLine,'Position',[211 Height-362 86 30])% rank 6 set(handles.phys_title,'Position',[38 Height-426 125 20])% rank 7 set(handles.CoordUnit,'Position',[151 Height-426 120 30])% rank 7 set(handles.px_title,'Position',[272 Height-426 125 20])% rank 7 set(handles.ListCoord,'Position',[1 20 418 Height-446])% rank 8 set(handles.geometry_calib,'Position',[Left Bottom 420 Height]) %set menu of calibration options set(handles.calib_type,'String',{'rescale';'linear';'3D_linear';'3D_quadr';'3D_extrinsic'}) if exist('inputfile','var')&& ~isempty(inputfile) [RootPath,SubDir,RootFile,tild,tild,tild,tild,FileExt]=fileparts_uvmat(inputfile); struct.XmlInputFile=find_imadoc(RootPath,SubDir,RootFile,FileExt); set(handles.ListCoord,'Data',[]) if exist(struct.XmlInputFile,'file') Heading=loadfile(handles,struct.XmlInputFile);% load data from the xml file and fill the GUI if isfield(Heading,'Campaign')&& ischar(Heading.Campaign) struct.Campaign=Heading.Campaign; end end set(hObject,'UserData',struct) end %------------------------------------------------------------------------ % --- Outputs from this function are returned to the command line. function varargout = geometry_calib_OutputFcn(~, eventdata, handles) %------------------------------------------------------------------------ % Get default command line output from handles structure varargout{1} = handles.output; varargout{2}=handles; % %------------------------------------------------------------------------ % executed when closing: set the parent interface button to value 0 function closefcn(gcbo,eventdata) %------------------------------------------------------------------------ huvmat=findobj(allchild(0),'Name','uvmat'); if ~isempty(huvmat) handles=guidata(huvmat); set(handles.MenuCalib,'Checked','off') hobject=findobj(handles.PlotAxes,'tag','calib_points'); if ~isempty(hobject) delete(hobject) end hobject=findobj(handles.PlotAxes,'tag','calib_marker'); if ~isempty(hobject) delete(hobject) end end %------------------------------------------------------------------------ % --- Executes on button press APPLY (used to launch the calibration). function APPLY_Callback(hObject, eventdata, handles) set(handles.CheckEnableMouse,'Value',0)% desactivate mouse (to avoid spurious creation of new points) %------------------------------------------------------------------------ %% look for the GUI uvmat and check for an image as input set(handles.APPLY,'BackgroundColor',[1 1 0])% paint APPLY button in yellow to show activation huvmat=findobj(allchild(0),'Name','uvmat');% look for the GUI uvmat hhuvmat=guidata(huvmat);%handles of elements in the GUI uvmat if ~strcmp(get(hhuvmat.Scalar,'Visible'),'on') msgbox_uvmat('ERROR','An image needs to be opened in uvmat for calibration') return end RootPath=''; if ~isempty(hhuvmat.RootPath)&& ~isempty(hhuvmat.RootFile) RootPath=get(hhuvmat.RootPath,'String');% path to the currently displayed image SubDirBase=regexprep(get(hhuvmat.SubDir,'String'),'\..+$',''); outputfile=[fullfile(RootPath,SubDirBase) '.xml'];%xml file associated with the currently displayed image else question={'save the calibration data and point coordinates in'}; def={fullfile(RootPath,'ObjectCalib.xml')}; options.Resize='on'; answer=inputdlg(question,'',1,def,options); outputfile=answer{1}; end %% read coordinates of the calibration poinnts: Coord(:,1-3) in phys, Coord(:,4-5) image Coord=get(handles.ListCoord,'Data'); %% read the type of calibration calib_cell=get(handles.calib_type,'String'); val=get(handles.calib_type,'Value'); CalibFcn=['calib_' calib_cell{val}]; %% read the intrinsic parameters Intrinsic.Npx=str2num(get(hhuvmat.num_Npx,'String')); Intrinsic.Npy=str2num(get(hhuvmat.num_Npy,'String')); Intrinsic.coord_files=get(handles.ListCoordFiles,'String'); Intrinsic.fx=str2num(get(handles.fx,'String')); Intrinsic.fy=str2num(get(handles.fy,'String')); Intrinsic.kc=str2num(get(handles.kc,'String')); Intrinsic.Cx=str2num(get(handles.Cx,'String')); Intrinsic.Cy=str2num(get(handles.Cy,'String')); if isempty(Intrinsic.kc) Intrinsic.kc=0; end if isempty(Intrinsic.Cx)||isempty(Intrinsic.Cy) Intrinsic.Cx=Intrinsic.Npx/2; Intrinsic.Cy=Intrinsic.Npy/2; end %% Apply calibration [GeometryCalib,index,ind_removed,Z_plane]=calibrate(Coord,CalibFcn,Intrinsic);% apply calibration %% record the coordinate unit unitlist=get(handles.CoordUnit,'String'); unit=unitlist{get(handles.CoordUnit,'value')}; GeometryCalib.CoordUnit=unit; %% record the coordinates of the calibration points GeometryCalib.SourceCalib.PointCoord=Coord; %% display calibration results on the GUI geometry_calib display_intrinsic(GeometryCalib,handles)%display calibration intrinsic parameters display_extrinsic(GeometryCalib,handles)%display calibration extrinsic parameters % (rotation and translation of camera with respect to the phys coordinates) %% store the calibration data, by default in the xml file of the currently displayed image answer=msgbox_uvmat('INPUT_Y-N',{'store calibration data';... ['Error rms (along x,y)=' num2str(GeometryCalib.ErrorRms) ' pixels'];... ['Error max (along x,y)=' num2str(GeometryCalib.ErrorMax) ' pixels']; [num2str(numel(ind_removed)) ' points removed']}); if strcmp(answer,'Yes') %store the calibration data if strcmp(calib_cell{val}(1:2),'3D')%set the plane position for 3D (projection) calibration msgbox_uvmat('CONFIRMATION',{['The current image series is assumed by default in the plane of the calib points z=' num2str(Z_plane) ] ; 'can be modified by MenuSetSlice in the upper bar menu of uvmat'}) GeometryCalib.SliceCoord=Z_plane'*[0 0 1]; end else GeometryCalib=[]; index=1; end if isempty(GeometryCalib) % if calibration cancelled set(handles.APPLY,'BackgroundColor',[1 0 1]) else % if calibration confirmed %% copy the xml file from the old location if appropriate, then update with the calibration parameters if ~exist(outputfile,'file') && ~isempty(SubDirBase) oldxml=[fullfile(RootPath,SubDirBase,get(hhuvmat.RootFile,'String')) '.xml']; if exist(oldxml,'file') [success,message]=copyfile(oldxml,outputfile);%copy the old xml file to a new one with the new convention end end errormsg=update_imadoc(GeometryCalib,outputfile,'GeometryCalib');% introduce the calibration data in the xml file if ~strcmp(errormsg,'') msgbox_uvmat('ERROR',errormsg); end %% display image with new calibration in the currently opened uvmat interface FieldList=get(hhuvmat.FieldName,'String'); val=get(hhuvmat.FieldName,'Value'); if strcmp(FieldList{val},'image') set(hhuvmat.CheckFixLimits,'Value',0)% put FixedLimits option to 'off' to plot the whole image UserData=get(handles.geometry_calib,'UserData'); UserData.XmlInputFile=outputfile;%save the current xml file name set(handles.geometry_calib,'UserData',UserData) uvmat('InputFileREFRESH_Callback',hObject,eventdata,hhuvmat); %file input with xml reading in uvmat, show the image in phys coordinates PLOT_Callback(hObject, eventdata, handles) set(handles.CoordLine,'string',num2str(index)) Coord=get(handles.ListCoord,'Data'); update_calib_marker(Coord(index,:)); %indicate the point with max deviations from phys coord to calibration figure(handles.geometry_calib)% put the GUI geometry_calib in front set(handles.APPLY,'BackgroundColor',[1 0 0]) % set APPLY button to red color else msgbox_uvmat('WARNING','open the image to see the effect of the new calibration') end end %------------------------------------------------------------------------ % --- Executes on button press in REPLICATE function REPLICATE_Callback(hObject, eventdata, handles) %------------------------------------------------------------------------ set(handles.CheckEnableMouse,'Value',0)% desactivate mouse (to avoid spurious creation of new points) %% look for the GUI uvmat and check for an image as input huvmat=findobj(allchild(0),'Name','uvmat'); hhuvmat=guidata(huvmat);%handles of elements in the GUI uvmat %% read coordinates of the calibration poinnts: Coord(:,1-3) in phys, Coord(:,4-5) image Coord=get(handles.ListCoord,'Data'); %% read the type of calibration calib_cell=get(handles.calib_type,'String'); val=get(handles.calib_type,'Value'); CalibFcn=['calib_' calib_cell{val}]; %% read the intrinsic parameters Intrinsic.Npx=str2num(get(hhuvmat.num_Npx,'String')); Intrinsic.Npy=str2num(get(hhuvmat.num_Npy,'String')); Intrinsic.coord_files=get(handles.ListCoordFiles,'String'); Intrinsic.fx=str2num(get(handles.fx,'String')); Intrinsic.fy=str2num(get(handles.fy,'String')); Intrinsic.kc=str2num(get(handles.kc,'String')); Intrinsic.Cx=str2num(get(handles.Cx,'String')); Intrinsic.Cy=str2num(get(handles.Cy,'String')); if isempty(Intrinsic.kc) Intrinsic.kc=0; end if isempty(Intrinsic.Cx)||isempty(Intrinsic.Cy) Intrinsic.Cx=Intrinsic.Npx/2; Intrinsic.Cy=Intrinsic.Npy/2; end %% apply to cropped images if requested answer=msgbox_uvmat('INPUT_Y-N','apply to cropped images?'); if strcmp(answer,'Yes') prompt = {'npy_lower'}; dlg_title = 'remove image the npy_lower image lines (removal of the upper linedoes not change calibration)'; num_lines= 1; def = {'0'}; answer = inputdlg(prompt,dlg_title,num_lines,def); npy_crop=str2num(answer{1}); Intrinsic.Npy=Intrinsic.Npy-npy_crop; %size of the filtering window Coord(:,5)=Coord(:,5)-npy_crop;% shift the image ordinates of the calibration points by removing the lower band end %% Apply calibration [GeometryCalib,index,ind_removed,Z_plane]=calibrate(Coord,CalibFcn,Intrinsic);% apply calibration %% record the coordinate unit unitlist=get(handles.CoordUnit,'String'); unit=unitlist{get(handles.CoordUnit,'value')}; GeometryCalib.CoordUnit=unit; %% record the coordinates of the calibration points GeometryCalib.SourceCalib.PointCoord=Coord; %% display calibration results on the GUI geometry_calib display_intrinsic(GeometryCalib,handles)%display calibration intrinsic parameters display_extrinsic(GeometryCalib,handles)%display calibration extrinsic parameters % (rotation and translation of camera with respect to the phys coordinates) %% store the calibration data, by default in the xml file of the currently displayed image answer=msgbox_uvmat('INPUT_Y-N',{'store calibration data';... ['Error rms (along x,y)=' num2str(GeometryCalib.ErrorRms) ' pixels'];... ['Error max (along x,y)=' num2str(GeometryCalib.ErrorMax) ' pixels'];... [num2str(numel(ind_removed)) ' points removed']}); if strcmp(answer,'Yes') %store the calibration data if strcmp(calib_cell{val}(1:2),'3D')%set the plane position for 3D (projection) calibration msgbox_uvmat('CONFIRMATION',{['The current image series is assumed by default in the plane of the calib points z=' num2str(Z_plane) ] ; 'can be modified by MenuSetSlice in the upper bar menu of uvmat'}) GeometryCalib.SliceCoord=Z_plane'*[0 0 1]; end else GeometryCalib=[]; index=1; end %% open the GUI browse_data CalibData=get(handles.geometry_calib,'UserData');%read the calibration image source on the interface userdata if isfield(CalibData,'XmlInputFile') InputDir=fileparts(fileparts(CalibData.XmlInputFile)); end DataSeries=uigetfile_uvmat('open a folder of images to calibrate',InputDir,'uigetdir'); %SubProject=uigetfile_uvmat('open folder of subproject to calibrate',InputDir,'uigetdir'); OutPut=browse_data(DataSeries); nbcalib=0; for ilist=1:numel(OutPut.Experiment) SubDirBase=regexprep(OutPut.DataSeries{1},'\..+$',''); XmlName=fullfile(OutPut.Campaign,OutPut.Experiment{ilist},[SubDirBase '.xml']); dispmessage=' created with calibration parameters'; % copy the xml file from the old location if appropriate, then update with the calibration parameters if ~exist(XmlName,'file') && ~isempty(SubDirBase) oldxml=fullfile(OutPut.Campaign,OutPut.Experiment{ilist},SubDirBase,[get(hhuvmat.RootFile,'String') '.xml']); if exist(oldxml,'file') [success,message]=copyfile(oldxml,XmlName);%copy the old xml file to a new one with the new convention dispmessage=' updated with calibration parameters'; end end errormsg=update_imadoc(GeometryCalib,XmlName,'GeometryCalib');% introduce the calibration data in the xml file if ~strcmp(errormsg,'') msgbox_uvmat('ERROR',errormsg); else display([XmlName dispmessage]) nbcalib=nbcalib+1; end end msgbox_uvmat('CONFIMATION',[SubDirBase ' calibrated for ' num2str(nbcalib) ' experiments']); %------------------------------------------------------------------------ % --- activate calibration and store parameters in ouputfile . function [GeometryCalib,ind_max,ind_removed,Z_plane]=calibrate(Coord,CalibFcn,Intrinsic) %------------------------------------------------------------------------ index=[]; GeometryCalib=[]; % apply the calibration, whose type is selected in handles.calib_type if ~isempty(Coord) GeometryCalib=feval(CalibFcn,Coord,Intrinsic); else msgbox_uvmat('ERROR','No calibration points, abort') end if isempty(GeometryCalib) return end Z_plane=[]; % estimate calibration error rms and max X=Coord(:,1); Y=Coord(:,2); Z=Coord(:,3); x_ima=Coord(:,4); y_ima=Coord(:,5); [Xpoints,Ypoints]=px_XYZ(GeometryCalib,X,Y,Z); GeometryCalib.ErrorRms(1)=sqrt(mean((Xpoints-x_ima).*(Xpoints-x_ima))); GeometryCalib.ErrorRms(2)=sqrt(mean((Ypoints-y_ima).*(Ypoints-y_ima))); [ErrorMax(1),index(1)]=max(abs(Xpoints-x_ima)); [ErrorMax(2),index(2)]=max(abs(Ypoints-y_ima)); [tild,ind_dim]=max(ErrorMax); ind_max=index(ind_dim); % mark the index with maximum deviation % detect bad calibration points, marked by indices ind_bad, if the % difference of actual image coordinates and those given by calibration is % greater than 2 pixels and greater than 4 rms. check_x=abs(Xpoints-x_ima)>max(2,3*GeometryCalib.ErrorRms(1)); check_y=abs(Ypoints-y_ima)>max(2,3*GeometryCalib.ErrorRms(2)); ind_removed=find(check_x | check_y) % repeat calibration without the excluded points: if ~isempty(ind_removed) Coord(ind_removed,:)=[]; GeometryCalib=feval(CalibFcn,Coord,Intrinsic); X=Coord(:,1); Y=Coord(:,2); Z=Coord(:,3); x_ima=Coord(:,4); y_ima=Coord(:,5); [Xpoints,Ypoints]=px_XYZ(GeometryCalib,X,Y,Z); GeometryCalib.ErrorRms(1)=sqrt(mean((Xpoints-x_ima).*(Xpoints-x_ima))); GeometryCalib.ErrorRms(2)=sqrt(mean((Ypoints-y_ima).*(Ypoints-y_ima))); end GeometryCalib.ErrorMax=ErrorMax; %set the Z position of the reference plane used for calibration if isequal(max(Z),min(Z))%Z constant Z_plane=Z(1); GeometryCalib.NbSlice=1; GeometryCalib.SliceCoord=[0 0 Z_plane]; end %------------------------------------------------------------------------ % --- determine the parameters for a calibration by an affine function (rescaling and offset, no rotation) function GeometryCalib=calib_rescale(Coord,Intrinsic) %------------------------------------------------------------------------ X=Coord(:,1); Y=Coord(:,2);% Z not used x_ima=Coord(:,4); y_ima=Coord(:,5); [px]=polyfit(X,x_ima,1); [py]=polyfit(Y,y_ima,1); % T_x=px(2); % T_y=py(2); GeometryCalib.CalibrationType='rescale'; GeometryCalib.fx_fy=[px(1) py(1)];%.fx_fy corresponds to pxcm along x and y GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function GeometryCalib.Tx_Ty_Tz=[px(2)/px(1) py(2)/py(1) 1]; GeometryCalib.omc=[0 0 0]; %------------------------------------------------------------------------ % --- determine the parameters for a calibration by a linear transform matrix (rescale and rotation) function GeometryCalib=calib_linear(Coord,Intrinsic) %------------------------------------------------------------------------ X=Coord(:,1); Y=Coord(:,2);% Z not used x_ima=Coord(:,4); y_ima=Coord(:,5); XY_mat=[ones(size(X)) X Y]; a_X1=XY_mat\x_ima; %transformation matrix for X a_Y1=XY_mat\y_ima;%transformation matrix for X R=[a_X1(2),a_X1(3);a_Y1(2),a_Y1(3)]; epsilon=sign(det(R)); norm=abs(det(R)); GeometryCalib.CalibrationType='linear'; if (a_X1(2)/a_Y1(3))>0 GeometryCalib.fx_fy(1)=sqrt((a_X1(2)/a_Y1(3))*norm); else GeometryCalib.fx_fy(1)=-sqrt(-(a_X1(2)/a_Y1(3))*norm); end GeometryCalib.fx_fy(2)=(a_Y1(3)/a_X1(2))*GeometryCalib.fx_fy(1); GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function GeometryCalib.Tx_Ty_Tz=[a_X1(1)/GeometryCalib.fx_fy(1) a_Y1(1)/GeometryCalib.fx_fy(2) 1]; R(1,:)=R(1,:)/GeometryCalib.fx_fy(1); R(2,:)=R(2,:)/GeometryCalib.fx_fy(2); R=[R;[0 0]]; GeometryCalib.R=[R [0;0;-epsilon]]; GeometryCalib.omc=(180/pi)*[acos(GeometryCalib.R(1,1)) 0 0]; %------------------------------------------------------------------------ % --- determine the tsai parameters for a view normal to the grid plane % NOT USED function GeometryCalib=calib_normal(Coord,Intrinsic) %------------------------------------------------------------------------ Calib.fx=str2num(get(handles.fx,'String')); Calib.fy=str2num(get(handles.fy,'String')); Calib.kc=str2num(get(handles.kc,'String')); Calib.Cx=str2num(get(handles.Cx,'String')); Calib.Cy=str2num(get(handles.Cy,'String')); %default if isempty(Calib.fx) Calib.fx=25/0.012; end if isempty(Calib.fy) Calib.fy=25/0.012; end if isempty(Calib.kc) Calib.kc=0; end if isempty(Calib.Cx)||isempty(Calib.Cy) huvmat=findobj(allchild(0),'Tag','uvmat'); hhuvmat=guidata(huvmat); Calib.Cx=str2num(get(hhuvmat.num_Npx,'String'))/2; Calib.Cx=str2num(get(hhuvmat.num_Npy,'String'))/2; end %tsai parameters Calib.dpx=0.012;%arbitrary Calib.dpy=0.012; Calib.sx=Calib.fx*Calib.dpx/(Calib.fy*Calib.dpy); Calib.f=Calib.fy*Calib.dpy; Calib.kappa1=Calib.kc/(Calib.f*Calib.f); %initial guess X=Coord(:,1); Y=Coord(:,2); Zmean=mean(Coord(:,3)); x_ima=Coord(:,4)-Calib.Cx; y_ima=Coord(:,5)-Calib.Cy; XY_mat=[ones(size(X)) X Y]; a_X1=XY_mat\x_ima; %transformation matrix for X a_Y1=XY_mat\y_ima;%transformation matrix for Y R=[a_X1(2),a_X1(3),0;a_Y1(2),a_Y1(3),0;0,0,-1];% rotation+ z axis reversal (upward) norm=sqrt(det(-R)); calib_param(1)=0;% quadratic distortion calib_param(2)=a_X1(1); calib_param(3)=a_Y1(1); calib_param(4)=Calib.f/(norm*Calib.dpx)-R(3,3)*Zmean; calib_param(5)=angle(a_X1(2)+1i*a_X1(3)); display(['initial guess=' num2str(calib_param)]) %optimise the parameters: minimisation of error calib_param = fminsearch(@(calib_param) error_calib(calib_param,Calib,Coord),calib_param); GeometryCalib.CalibrationType='tsai_normal'; GeometryCalib.focal=Calib.f; GeometryCalib.dpx_dpy=[Calib.dpx Calib.dpy]; GeometryCalib.Cx_Cy=[Calib.Cx Calib.Cy]; GeometryCalib.sx=Calib.sx; GeometryCalib.kappa1=calib_param(1); GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function GeometryCalib.Tx_Ty_Tz=[calib_param(2) calib_param(3) calib_param(4)]; alpha=calib_param(5); GeometryCalib.R=[cos(alpha) sin(alpha) 0;-sin(alpha) cos(alpha) 0;0 0 -1]; %------------------------------------------------------------------------ function GeometryCalib=calib_3D_linear(Coord,Intrinsic) %------------------------------------------------------------------------ coord_files=Intrinsic.coord_files; if ischar(Intrinsic.coord_files) coord_files={coord_files}; end if isempty(coord_files{1}) || isequal(coord_files,{''}) coord_files={}; end %retrieve the calibration points stored in the files listed in the popup list ListCoordFiles x_1=Coord(:,4:5)';%px coordinates of the ref points nx=Intrinsic.Npx; ny=Intrinsic.Npy; x_1(2,:)=ny-x_1(2,:);%reverse the y image coordinates X_1=Coord(:,1:3)';%phys coordinates of the ref points n_ima=numel(coord_files)+1; if ~isempty(coord_files) msgbox_uvmat('CONFIRMATION',['The xy coordinates of the calibration points in ' num2str(n_ima) ' planes will be used']) for ifile=1:numel(coord_files) t=xmltree(coord_files{ifile}); s=convert(t);%convert to matlab structure if isfield(s,'GeometryCalib') if isfield(s.GeometryCalib,'SourceCalib') if isfield(s.GeometryCalib.SourceCalib,'PointCoord') PointCoord=s.GeometryCalib.SourceCalib.PointCoord; Coord_file=zeros(length(PointCoord),5);%default for i=1:length(PointCoord) line=str2num(PointCoord{i}); Coord_file(i,4:5)=line(4:5);%px x Coord_file(i,1:3)=line(1:3);%phys x end eval(['x_' num2str(ifile+1) '=Coord_file(:,4:5)'';']); eval(['x_' num2str(ifile+1) '(2,:)=ny-x_' num2str(ifile+1) '(2,:);' ]); eval(['X_' num2str(ifile+1) '=Coord_file(:,1:3)'';']); end end end end end n_ima=numel(coord_files)+1; est_dist=[0;0;0;0;0]; est_aspect_ratio=0; est_fc=[1;1]; center_optim=0; path_uvmat=which('uvmat');% check the path detected for source file uvmat path_UVMAT=fileparts(path_uvmat); %path to UVMAT run(fullfile(path_UVMAT,'toolbox_calib','go_calib_optim'));% apply fct 'toolbox_calib/go_calib_optim' if exist('Rc_1','var') GeometryCalib.CalibrationType='3D_linear'; GeometryCalib.fx_fy=fc'; GeometryCalib.Cx_Cy=cc'; GeometryCalib.kc=kc(1); GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function GeometryCalib.Tx_Ty_Tz=Tc_1'; GeometryCalib.R=Rc_1; GeometryCalib.R(2,1:3)=-GeometryCalib.R(2,1:3);%inversion of the y image coordinate GeometryCalib.Tx_Ty_Tz(2)=-GeometryCalib.Tx_Ty_Tz(2);%inversion of the y image coordinate GeometryCalib.Cx_Cy(2)=ny-GeometryCalib.Cx_Cy(2);%inversion of the y image coordinate GeometryCalib.omc=(180/pi)*omc_1;%angles in degrees GeometryCalib.ErrorRMS=[]; GeometryCalib.ErrorMax=[]; else msgbox_uvmat('ERROR',['calibration function ' fullfile('toolbox_calib','go_calib_optim') ' did not converge: use multiple views or option 3D_extrinsic']) GeometryCalib=[]; end %------------------------------------------------------------------------ function GeometryCalib=calib_3D_quadr(Coord,Intrinsic) %------------------------------------------------------------------ path_uvmat=which('uvmat');% check the path detected for source file uvmat path_UVMAT=fileparts(path_uvmat); %path to UVMAT huvmat=findobj(allchild(0),'Tag','uvmat'); hhuvmat=guidata(huvmat); if ~strcmp(get(hhuvmat.Scalar,'Visible'),'on') msgbox_uvmat('ERROR','An image needs to be opened in uvmat for calibration') return end % check_cond=0; coord_files=Intrinsic.coord_files; if ischar(coord_files) coord_files={coord_files}; end if isempty(coord_files{1}) || isequal(coord_files,{''}) coord_files={}; end %retrieve the calibration points stored in the files listed in the popup list ListCoordFiles x_1=Coord(:,4:5)';%px coordinates of the ref points nx=str2num(get(hhuvmat.num_Npx,'String')); ny=str2num(get(hhuvmat.num_Npy,'String')); x_1(2,:)=ny-x_1(2,:);%reverse the y image coordinates X_1=Coord(:,1:3)';%phys coordinates of the ref points n_ima=numel(coord_files)+1; if ~isempty(coord_files) msgbox_uvmat('CONFIRMATION',['The xy coordinates of the calibration points in ' num2str(n_ima) ' planes will be used']) for ifile=1:numel(coord_files) t=xmltree(coord_files{ifile}); s=convert(t);%convert to matlab structure if isfield(s,'GeometryCalib') if isfield(s.GeometryCalib,'SourceCalib') if isfield(s.GeometryCalib.SourceCalib,'PointCoord') PointCoord=s.GeometryCalib.SourceCalib.PointCoord; Coord_file=zeros(length(PointCoord),5);%default for i=1:length(PointCoord) line=str2num(PointCoord{i}); Coord_file(i,4:5)=line(4:5);%px x Coord_file(i,1:3)=line(1:3);%phys x end eval(['x_' num2str(ifile+1) '=Coord_file(:,4:5)'';']); eval(['x_' num2str(ifile+1) '(2,:)=ny-x_' num2str(ifile+1) '(2,:);' ]); eval(['X_' num2str(ifile+1) '=Coord_file(:,1:3)'';']); end end end end end n_ima=numel(coord_files)+1; est_dist=[1;0;0;0;0]; est_aspect_ratio=1; center_optim=0; run(fullfile(path_UVMAT,'toolbox_calib','go_calib_optim'));% apply fct 'toolbox_calib/go_calib_optim' if exist('Rc_1','var') GeometryCalib.CalibrationType='3D_quadr'; GeometryCalib.fx_fy=fc'; GeometryCalib.Cx_Cy=cc'; GeometryCalib.kc=kc(1); GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function GeometryCalib.Tx_Ty_Tz=Tc_1'; GeometryCalib.R=Rc_1; GeometryCalib.R(2,1:3)=-GeometryCalib.R(2,1:3);%inversion of the y image coordinate GeometryCalib.Tx_Ty_Tz(2)=-GeometryCalib.Tx_Ty_Tz(2);%inversion of the y image coordinate GeometryCalib.Cx_Cy(2)=ny-GeometryCalib.Cx_Cy(2);%inversion of the y image coordinate GeometryCalib.omc=(180/pi)*omc_1;%angles in degrees GeometryCalib.ErrorRMS=[]; GeometryCalib.ErrorMax=[]; else msgbox_uvmat('ERROR',['calibration function ' fullfile('toolbox_calib','go_calib_optim') ' did not converge: use multiple views or option 3D_extrinsic']) GeometryCalib=[]; end %------------------------------------------------------------------------ function GeometryCalib=calib_3D_extrinsic(Coord, Intrinsic) %------------------------------------------------------------------ path_uvmat=which('geometry_calib');% check the path detected for source file uvmat path_UVMAT=fileparts(path_uvmat); %path to UVMAT x_1=double(Coord(:,4:5)');%image coordinates X_1=double(Coord(:,1:3)');% phys coordinates huvmat=findobj(allchild(0),'Tag','uvmat'); hhuvmat=guidata(huvmat); if ~strcmp(get(hhuvmat.Scalar,'Visible'),'on') msgbox_uvmat('ERROR','An image needs to be opened in uvmat for calibration') return end ny=str2double(get(hhuvmat.num_Npy,'String')); x_1(2,:)=ny-x_1(2,:);%reverse the y image coordinates n_ima=1; GeometryCalib.CalibrationType='3D_extrinsic'; fx=Intrinsic.fx; fy=Intrinsic.fy; Cx=Intrinsic.Cx; Cy=Intrinsic.Cy; kc=Intrinsic.kc; errormsg=''; if isempty(fx) errormsg='focal length fx needs to be introduced'; elseif isempty(fy) errormsg='focal length fy needs to be introduced'; elseif isempty(Cx) errormsg='shift Cx to image centre needs to be introduced'; elseif isempty(Cy) errormsg='shift Cy to image centre needs to be introduced'; end if ~isempty(errormsg) GeometryCalib=[]; msgbox_uvmat('ERROR',errormsg) return end GeometryCalib.fx_fy(1)=fx; GeometryCalib.fx_fy(2)=fy; GeometryCalib.Cx_Cy(1)=Cx; GeometryCalib.Cx_Cy(2)=Cy; GeometryCalib.kc=kc; fct_path=fullfile(path_UVMAT,'toolbox_calib'); addpath(fct_path) GeometryCalib.Cx_Cy(2)=ny-GeometryCalib.Cx_Cy(2);%reverse Cx_Cy(2) for calibration (inversion of px ordinate) [omc,Tc1,Rc1,H,x,ex,JJ] = compute_extrinsic(x_1,X_1,... (GeometryCalib.fx_fy)',GeometryCalib.Cx_Cy',[GeometryCalib.kc 0 0 0 0]); rmpath(fct_path); GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function GeometryCalib.Tx_Ty_Tz=Tc1'; %inversion of z axis GeometryCalib.R=Rc1; GeometryCalib.R(2,1:3)=-GeometryCalib.R(2,1:3);%inversion of the y image coordinate GeometryCalib.Tx_Ty_Tz(2)=-GeometryCalib.Tx_Ty_Tz(2);%inversion of the y image coordinate GeometryCalib.Cx_Cy(2)=ny-GeometryCalib.Cx_Cy(2);%inversion of the y image coordinate GeometryCalib.omc=(180/pi)*omc'; %------------------------------------------------------------------------ %function GeometryCalib=calib_tsai_heikkila(Coord) % TEST: NOT IMPLEMENTED %------------------------------------------------------------------ % path_uvmat=which('uvmat');% check the path detected for source file uvmat % path_UVMAT=fileparts(path_uvmat); %path to UVMAT % path_calib=fullfile(path_UVMAT,'toolbox_calib_heikkila'); % addpath(path_calib) % npoints=size(Coord,1); % Coord(:,1:3)=10*Coord(:,1:3); % Coord=[Coord zeros(npoints,2) -ones(npoints,1)]; % [par,pos,iter,res,er,C]=cacal('dalsa',Coord); % GeometryCalib.CalibrationType='tsai'; % GeometryCalib.focal=par(2); %------------------------------------------------------------------------ % --- determine the rms of calibration error function ErrorRms=error_calib(calib_param,Calib,Coord) %------------------------------------------------------------------------ %calib_param: vector of free calibration parameters (to optimise) %Calib: structure of the given calibration parameters %Coord: list of phys coordinates (columns 1-3, and pixel coordinates (columns 4-5) Calib.f=25; Calib.dpx=0.012; Calib.dpy=0.012; Calib.sx=1; Calib.Cx=512; Calib.Cy=512; Calib.kappa1=calib_param(1); Calib.Tx=calib_param(2); Calib.Ty=calib_param(3); Calib.Tz=calib_param(4); alpha=calib_param(5); Calib.R=[cos(alpha) sin(alpha) 0;-sin(alpha) cos(alpha) 0;0 0 -1]; X=Coord(:,1); Y=Coord(:,2); Z=Coord(:,3); x_ima=Coord(:,4); y_ima=Coord(:,5); [Xpoints,Ypoints]=px_XYZ(Calib,X,Y,Z); ErrorRms(1)=sqrt(mean((Xpoints-x_ima).*(Xpoints-x_ima))); ErrorRms(2)=sqrt(mean((Ypoints-y_ima).*(Ypoints-y_ima))); ErrorRms=mean(ErrorRms); %------------------------------------------------------------------------ % --- Executes on button press in STORE. function STORE_Callback(hObject, eventdata, handles) %------------------------------------------------------------------------ Coord=get(handles.ListCoord,'Data'); %Object=read_geometry_calib(Coord_cell); unitlist=get(handles.CoordUnit,'String'); unit=unitlist{get(handles.CoordUnit,'value')}; GeometryCalib.CoordUnit=unit; GeometryCalib.SourceCalib.PointCoord=Coord(:,1:5); huvmat=findobj(allchild(0),'Name','uvmat'); hhuvmat=guidata(huvmat);%handles of elements in the GUI uvmat % RootPath=''; % RootFile=''; if ~isempty(hhuvmat.RootPath)&& ~isempty(hhuvmat.RootFile) % testhandle=1; RootPath=get(hhuvmat.RootPath,'String'); RootFile=get(hhuvmat.RootFile,'String'); filebase=[fullfile(RootPath,RootFile) '~']; while exist([filebase '.xml'],'file') filebase=[filebase '~']; end outputfile=[filebase '.xml']; errormsg=update_imadoc(GeometryCalib,outputfile,'GeometryCalib'); if ~strcmp(errormsg,'') msgbox_uvmat('ERROR',errormsg); end listfile=get(handles.ListCoordFiles,'string'); if isequal(listfile,{''}) listfile={outputfile}; else listfile=[listfile;{outputfile}];%update the list of coord files end set(handles.ListCoordFiles,'string',listfile); end ClearAll_Callback(hObject, eventdata, handles)% clear the current list and point plots % -------------------------------------------------------------------- % --- Executes on button press in ClearAll: clear the list of calibration points function ClearAll_Callback(hObject, eventdata, handles) % -------------------------------------------------------------------- set(handles.ListCoord,'Data',[]) PLOT_Callback(hObject, eventdata, handles) update_calib_marker([]);%remove circle marker set(handles.APPLY,'backgroundColor',[1 0 0])% set APPLY button to red color: no calibration wanted without points set(handles.CheckEnableMouse,'value',1); %activate mouse to pick new points %------------------------------------------------------------------------ % --- Executes on button press in CLEAR LIST function CLEAR_Callback(hObject, eventdata, handles) %------------------------------------------------------------------------ set(handles.ListCoordFiles,'Value',1) set(handles.ListCoordFiles,'String',{''}) %------------------------------------------------------------------------ % --- Executes on selection change in CheckEnableMouse. function CheckEnableMouse_Callback(hObject, eventdata, handles) %------------------------------------------------------------------------ choice=get(handles.CheckEnableMouse,'Value'); if choice huvmat=findobj(allchild(0),'tag','uvmat'); if ishandle(huvmat) hhuvmat=guidata(huvmat); set(hhuvmat.MenuRuler,'checked','off')%desactivate ruler if get(hhuvmat.CheckEditObject,'Value') set(hhuvmat.CheckEditObject,'Value',0) uvmat('CheckEditObject_Callback',hhuvmat.CheckEditObject,[],hhuvmat) end set(hhuvmat.MenuRuler,'checked','off')%desactivate ruler end end % -------------------------------------------------------------------- function MenuHelp_Callback(hObject, eventdata, handles) web('http://servforge.legi.grenoble-inp.fr/projects/soft-uvmat/wiki/UvmatHelp#GeometryCalib') % -------------------------------------------------------------------- function MenuSetScale_Callback(hObject, eventdata, handles) answer=msgbox_uvmat('INPUT_TXT','scale pixel/cm?',''); %create test points huvmat=findobj(allchild(0),'Name','uvmat');%find the current uvmat interface handle hhuvmat=guidata(huvmat); if ~strcmp(get(hhuvmat.Scalar,'Visible'),'on') msgbox_uvmat('ERROR','An image needs to be opened in uvmat for calibration') return end set(handles.calib_type,'Value',1)% set calib option to 'rescale' npx=str2num(get(hhuvmat.num_Npx,'String'))/2; npy=str2num(get(hhuvmat.num_Npy,'String'))/2; Xima=[0.25*npx 0.75*npx 0.75*npx 0.25*npx]'; Yima=[0.25*npy 0.25*npy 0.75*npy 0.75*npy]'; x=Xima/str2num(answer); y=Yima/str2num(answer); Coord=[x y zeros(4,1) Xima Yima zeros(4,1)]; set(handles.ListCoord,'Data',Coord) set(handles.APPLY,'BackgroundColor',[1 0 1]) set(handles.CheckEnableMouse,'value',0); %desactivate mouse to avoid spurious points %------------------------------------------------------------------------ function MenuCreateGrid_Callback(hObject, eventdata, handles) %------------------------------------------------------------------------ CalibData=get(handles.geometry_calib,'UserData'); Tinput=[];%default if isfield(CalibData,'grid') Tinput=CalibData.grid; end [T,CalibData.grid]=create_grid(Tinput);%display the GUI create_grid set(handles.geometry_calib,'UserData',CalibData) %grid in phys space Coord=get(handles.ListCoord,'Data'); Coord(1:size(T,1),1:3)=T;%update the existing list of phys coordinates from the GUI create_grid set(handles.ListCoord,'Data',Coord) set(handles.APPLY,'BackgroundColor',[1 0 1]) set(handles.CheckEnableMouse,'value',0); %desactivate mouse to avoid spurious points % ----------------------------------------------------------------------- % --- automatic grid dectection from local maxima of the images function MenuDetectGrid_Callback(hObject, eventdata, handles) %------------------------------------------------------------------------ %% read the four last point coordinates in pixels Coord=get(handles.ListCoord,'Data');%read list of coordinates on geometry_calib nbpoints=size(Coord,1); %nbre of calibration points if nbpoints~=4 msgbox_uvmat('ERROR','four points must have be selected by the mouse to delimitate the phys grid area; the Ox axis will be defined by the two first points') return end % corners_X=(Coord(end:-1:end-3,4)); %pixel absissa of the four corners % corners_Y=(Coord(end:-1:end-3,5)); corners_X=(Coord(:,4)); %pixel absissa of the four corners corners_Y=(Coord(:,5)); %%%%%% % corners_X=1000*[1.5415 1.7557 1.7539 1.5415]'; % corners_Y=1000*[1.1515 1.1509 1.3645 1.3639]'; %reorder the last two points (the two first in the list) if needed angles=angle((corners_X-corners_X(1))+1i*(corners_Y-corners_Y(1))); if abs(angles(4)-angles(2))>abs(angles(3)-angles(2)) X_end=corners_X(4); Y_end=corners_Y(4); corners_X(4)=corners_X(3); corners_Y(4)=corners_Y(3); corners_X(3)=X_end; corners_Y(3)=Y_end; end %% initiate the grid in phys coordinates CalibData=get(handles.geometry_calib,'UserData');%get information stored on the GUI geometry_calib grid_input=[];%default if isfield(CalibData,'grid') grid_input=CalibData.grid;%retrieve the previously used grid end [T,CalibData.grid,CalibData.grid.CheckWhite,CalibData.grid.FilterWindow]=create_grid(grid_input,'detect_grid');%display the GUI create_grid, read the set of phys coordinates T set(handles.geometry_calib,'UserData',CalibData)%store the phys grid parameters for later use X=[CalibData.grid.x_0 CalibData.grid.x_1 CalibData.grid.x_0 CalibData.grid.x_1]';%corner absissa in the phys coordinates (cm) Y=[CalibData.grid.y_0 CalibData.grid.y_0 CalibData.grid.y_1 CalibData.grid.y_1]';%corner ordinates in the phys coordinates (cm) %% read the current image, displayed in the GUI uvmat huvmat=findobj(allchild(0),'Name','uvmat'); UvData=get(huvmat,'UserData'); A=UvData.Field.A;%currently displayed image npxy=size(A); %% calculate transform matrices for plane projection: rectangle assumed to be viewed in perspective % reference: http://alumni.media.mit.edu/~cwren/interpolator/ by Christopher R. Wren B = [ X Y ones(size(X)) zeros(4,3) -X.*corners_X -Y.*corners_X ... zeros(4,3) X Y ones(size(X)) -X.*corners_Y -Y.*corners_Y ]; B = reshape (B', 8 , 8 )'; D = [ corners_X , corners_Y ]; D = reshape (D', 8 , 1 ); l = (B' * B)\B' * D; Amat = reshape([l(1:6)' 0 0 1 ],3,3)'; C = [l(7:8)' 1]; %% transform grid image into 'phys' coordinates GeometryCalib.CalibrationType='3D_linear'; GeometryCalib.fx_fy=[1 1]; GeometryCalib.Tx_Ty_Tz=[Amat(1,3) Amat(2,3) 1]; GeometryCalib.R=[Amat(1,1),Amat(1,2),0;Amat(2,1),Amat(2,2),0;C(1),C(2),0]; GeometryCalib.CoordUnit='cm'; path_uvmat=which('uvmat');% check the path detected for source file uvmat path_UVMAT=fileparts(path_uvmat); %path to UVMAT addpath(fullfile(path_UVMAT,'transform_field')) Data.ListVarName={'Coord_y','Coord_x','A'}; Data.VarDimName={'Coord_y','Coord_x',{'Coord_y','Coord_x'}}; if ndims(A)==3 A=mean(A,3); end Data.A=A-min(min(A)); Data.Coord_y=[npxy(1)-0.5 0.5]; Data.Coord_x=[0.5 npxy(2)]; Data.CoordUnit='pixel'; Calib.GeometryCalib=GeometryCalib; DataOut=phys(Data,Calib); rmpath(fullfile(path_UVMAT,'transform_field')) Amod=DataOut.A;% current image expressed in 'phys' coord Rangx=DataOut.Coord_x;% x coordinates of first and last pixel centres in phys Rangy=DataOut.Coord_y;% y coordinates of first and last pixel centres in phys if CalibData.grid.CheckWhite Amod=double(Amod);%case of white grid markers: will look for image maxima else Amod=-double(Amod);%case of black grid markers: will look for image minima end %%%%%%filterfor i;proved detection of dots if ~isequal(CalibData.grid.FilterWindow,0) %definition of the cos shape matrix filter FilterBoxSize_x=CalibData.grid.FilterWindow; FilterBoxSize_y=CalibData.grid.FilterWindow; ix=1/2-FilterBoxSize_x/2:-1/2+FilterBoxSize_x/2;% iy=1/2-FilterBoxSize_y/2:-1/2+FilterBoxSize_y/2;% %del=np/3; %fct=exp(-(ix/del).^2); fct2_x=cos(ix/((FilterBoxSize_x-1)/2)*pi/2); fct2_y=cos(iy/((FilterBoxSize_y-1)/2)*pi/2); %Mfiltre=(ones(5,5)/5^2); Mfiltre=fct2_y'*fct2_x; Mfiltre=Mfiltre/(sum(sum(Mfiltre)));%normalize filter if ndims(Amod)==3 Amod=filter2(Mfiltre,sum(Amod,3));%filter the input image, after summation on the color component (for color images) else Amod=filter2(Mfiltre,Amod); end end %%%%%%%%%%%%%% %% detection of local image extrema in each direction Dx=(Rangx(2)-Rangx(1))/(npxy(2)-1); %x mesh in real space Dy=(Rangy(2)-Rangy(1))/(npxy(1)-1); %y mesh in real space ind_range_x=ceil(abs(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 ind_range_y=ceil(abs(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 nbpoints=size(T,1); TIndex=ones(size(T));% image indices corresponding to point coordinates %look for image maxima around each expected grid point for ipoint=1:nbpoints i0=1+round((T(ipoint,1)-Rangx(1))/Dx);% x index of the expected point in the phys image Amod j0=1+round((T(ipoint,2)-Rangy(1))/Dy);% y index of the expected point in the phys image Amod j0min=max(j0-ind_range_y,1);% min y index selected for the subimage (cut at the edge to avoid index <1) j0max=min(j0+ind_range_y,size(Amod,1));% max y index selected for the subimage (cut at the edge to avoid index > size) i0min=max(i0-ind_range_x,1);% min x index selected for the subimage (cut at the edge to avoid index <1) i0max=min(i0+ind_range_x,size(Amod,2));% max x index selected for the subimage (cut at the edge to avoid index > size) Asub=Amod(j0min:j0max,i0min:i0max); %subimage used to find brigthness extremum x_profile=sum(Asub,1);%profile of subimage summed over y y_profile=sum(Asub,2);%profile of subimage summed over x [tild,ind_x_max]=max(x_profile);% index of max for the x profile [tild,ind_y_max]=max(y_profile);% index of max for the y profile %sub-pixel improvement using moments x_shift=0; y_shift=0; if ind_x_max+2<=numel(x_profile) && ind_x_max-2>=1 Atop=x_profile(ind_x_max-2:ind_x_max+2);% extract x profile around the max x_shift=sum(Atop.*[-2 -1 0 1 2])/sum(Atop); end if ind_y_max+2<=numel(y_profile) && ind_y_max-2>=1 Atop=y_profile(ind_y_max-2:ind_y_max+2);% extract y profile around the max y_shift=sum(Atop.*[-2 -1 0 1 2]')/sum(Atop); end %%%% % if ipoint==9 % figure(11) % imagesc(Asub) % figure(12) % plot(x_profile,'r') % hold on % plot(y_profile,'b') % grid on % end %%%% TIndex(ipoint,1)=(i0min+ind_x_max-1+x_shift);% x position of the maximum (in index of Amod) TIndex(ipoint,2)=(j0min+ind_y_max-1+y_shift);% y position of the maximum (in index of Amod) end Tmod(:,1)=(TIndex(:,1)-1)*Dx+Rangx(1); Tmod(:,2)=(TIndex(:,2)-1)*Dy+Rangy(1); %Tmod=T(:,(1:2))+Delta;% 'phys' coordinates of the detected points [Xpx,Ypx]=px_XYZ(GeometryCalib,Tmod(:,1),Tmod(:,2));% image coordinates of the detected points Coord=[T Xpx Ypx zeros(size(T,1),1)]; set(handles.ListCoord,'Data',Coord) PLOT_Callback(hObject, eventdata, handles) set(handles.APPLY,'BackgroundColor',[1 0 1]) set(handles.CheckEnableMouse,'value',0); %desactivate mouse to avoid spurious points %----------------------------------------------------------------------- function MenuTranslatePoints_Callback(hObject, eventdata, handles) %----------------------------------------------------------------------- %hcalib=get(handles.calib_type,'parent');%handles of the GUI geometry_calib CalibData=get(handles.geometry_calib,'UserData'); Tinput=[];%default if isfield(CalibData,'translate') Tinput=CalibData.translate; end T=translate_points(Tinput);%display translate_points GUI and get shift parameters CalibData.translate=T; set(handles.geometry_calib,'UserData',CalibData) %translation Coord=get(handles.ListCoord,'Data'); Coord(:,1)=T(1)+Coord(:,1); Coord(:,2)=T(2)+Coord(:,2); Coord(:,3)=T(3)+Coord(:,3); set(handles.ListCoord,'Data',Coord); set(handles.APPLY,'BackgroundColor',[1 0 1]) % -------------------------------------------------------------------- function MenuRotatePoints_Callback(hObject, eventdata, handles) %hcalib=get(handles.calib_type,'parent');%handles of the GUI geometry_calib CalibData=get(handles.geometry_calib,'UserData'); Tinput=[];%default if isfield(CalibData,'rotate') Tinput=CalibData.rotate; end T=rotate_points(Tinput);%display rotate_points GUI to introduce rotation parameters CalibData.rotate=T; set(handles.geometry_calib,'UserData',CalibData) %----------------------------------------------------- %rotation Phi=T(1); O_x=0;%default O_y=0;%default if numel(T)>=2 O_x=T(2);%default end if numel(T)>=3 O_y=T(3);%default end Coord=get(handles.ListCoord,'Data'); r1=cos(pi*Phi/180); r2=-sin(pi*Phi/180); r3=sin(pi*Phi/180); r4=cos(pi*Phi/180); x=Coord(:,1)-O_x; y=Coord(:,2)-O_y; Coord(:,1)=r1*x+r2*y; Coord(:,2)=r3*x+r4*y; set(handles.ListCoord,'Data',Coord) set(handles.APPLY,'BackgroundColor',[1 0 1]) % -------------------------------------------------------------------- function MenuFlip_x_Callback(hObject, eventdata, handles) Coord=get(handles.ListCoord,'Data'); Coord(:,1)=-Coord(:,1); set(handles.ListCoord,'Data',Coord) % -------------------------------------------------------------------- function MenuFlip_y_Callback(hObject, eventdata, handles) Coord=get(handles.ListCoord,'Data'); Coord(:,2)=-Coord(:,2); set(handles.ListCoord,'Data',Coord) % -------------------------------------------------------------------- function MenuImportPoints_Callback(hObject, eventdata, handles) fileinput=browse_xml(hObject, eventdata, handles); if isempty(fileinput) return end [s,errormsg]=imadoc2struct(fileinput,'GeometryCalib'); if ~isfield(s,'GeometryCalib') msgbox_uvmat('ERROR','invalid input file: no geometry_calib data') return end GeometryCalib=s.GeometryCalib; if ~(isfield(GeometryCalib,'SourceCalib')&&isfield(GeometryCalib.SourceCalib,'PointCoord')) msgbox_uvmat('ERROR','invalid input file: no calibration points') return end Coord=GeometryCalib.SourceCalib.PointCoord; Coord=[Coord zeros(size(Coord,1),1)]; set(handles.ListCoord,'Data',Coord) PLOT_Callback(handles.geometry_calib, [], handles) set(handles.APPLY,'BackgroundColor',[1 0 1]) set(handles.CheckEnableMouse,'value',0); %desactivate mouse to avoid modifications by default % ----------------------------------------------------------------------- function MenuImportIntrinsic_Callback(hObject, eventdata, handles) %------------------------------------------------------------------------ fileinput=browse_xml(hObject, eventdata, handles); if isempty(fileinput) return end [s,errormsg]=imadoc2struct(fileinput,'GeometryCalib'); GeometryCalib=s.GeometryCalib; display_intrinsic(GeometryCalib,handles) % ----------------------------------------------------------------------- function MenuImportAll_Callback(hObject, eventdata, handles) %------------------------------------------------------------------------ fileinput=browse_xml(hObject, eventdata, handles); if ~isempty(fileinput) loadfile(handles,fileinput) end set(handles.CheckEnableMouse,'value',0); %desactivate mouse to avoid modifications by default % ----------------------------------------------------------------------- % --- Executes on menubar option Import/Grid file: introduce previous grid files function MenuGridFile_Callback(hObject, eventdata, handles) % ----------------------------------------------------------------------- inputfile=browse_xml(hObject, eventdata, handles); listfile=get(handles.ListCoordFiles,'String'); if isequal(listfile,{''}) listfile={inputfile}; else listfile=[listfile;{inputfile}];%update the list of coord files end set(handles.ListCoordFiles,'String',listfile); %------------------------------------------------------------------------ function fileinput=browse_xml(hObject, eventdata, handles) %------------------------------------------------------------------------ fileinput=[];%default oldfile=''; %default UserData=get(handles.geometry_calib,'UserData'); if isfield(UserData,'XmlInputFile') oldfile=UserData.XmlInputFile; end [FileName, PathName, filterindex] = uigetfile( ... {'*.xml;*.mat', ' (*.xml,*.mat)'; '*.xml', '.xml files '; ... '*.mat', '.mat matlab files '}, ... 'Pick a file',oldfile); fileinput=[PathName FileName];%complete file name testblank=findstr(fileinput,' ');%look for blanks if ~isempty(testblank) msgbox_uvmat('ERROR','forbidden input file name or path: no blank character allowed') return end sizf=size(fileinput); if (~ischar(fileinput)||~isequal(sizf(1),1)),return;end UserData.XmlInputFile=fileinput; set(handles.geometry_calib,'UserData',UserData)%record current file foer further use of browser % ----------------------------------------------------------------------- function Heading=loadfile(handles,fileinput) %------------------------------------------------------------------------ Heading=[];%default [s,errormsg]=imadoc2struct(fileinput,'Heading','GeometryCalib'); if ~isempty(errormsg) msgbox_uvmat('ERROR',errormsg) return end if ~isempty(s.Heading) Heading=s.Heading; end GeometryCalib=s.GeometryCalib; fx=1;fy=1;Cx=0;Cy=0;kc=0; %default CoordCell={}; Tabchar={};%default val_cal=1;%default if ~isempty(GeometryCalib) % choose the calibration option if isfield(GeometryCalib,'CalibrationType') calib_list=get(handles.calib_type,'String'); for ilist=1:numel(calib_list) if strcmp(calib_list{ilist},GeometryCalib.CalibrationType) val_cal=ilist; break end end end display_intrinsic(GeometryCalib,handles)%intrinsic param %extrinsic param if isfield(GeometryCalib,'Tx_Ty_Tz') Tx_Ty_Tz=GeometryCalib.Tx_Ty_Tz; set(handles.Tx,'String',num2str(GeometryCalib.Tx_Ty_Tz(1),4)) set(handles.Ty,'String',num2str(GeometryCalib.Tx_Ty_Tz(2),4)) set(handles.Tz,'String',num2str(GeometryCalib.Tx_Ty_Tz(3),4)) end if isfield(GeometryCalib,'omc') set(handles.Phi,'String',num2str(GeometryCalib.omc(1),4)) set(handles.Theta,'String',num2str(GeometryCalib.omc(2),4)) set(handles.Psi,'String',num2str(GeometryCalib.omc(3),4)) end if isfield(GeometryCalib,'SourceCalib')&& isfield(GeometryCalib.SourceCalib,'PointCoord') calib=GeometryCalib.SourceCalib.PointCoord; Coord=[calib zeros(size(calib,1),1)]; set(handles.ListCoord,'Data',Coord) end PLOT_Callback(handles.geometry_calib, [], handles) set(handles.APPLY,'BackgroundColor',[1 0 1]) end set(handles.calib_type,'Value',val_cal) if isempty(CoordCell)% allow mouse action by default in the absence of input points set(handles.CheckEnableMouse,'Value',1) set(handles.CheckEnableMouse,'BackgroundColor',[1 1 0]) else % does not allow mouse action by default in the presence of input points set(handles.CheckEnableMouse,'Value',0) set(handles.CheckEnableMouse,'BackgroundColor',[0.7 0.7 0.7]) end %------------------------------------------------------------------------ %---display calibration intrinsic parameters function display_intrinsic(GeometryCalib,handles) %------------------------------------------------------------------------ fx=[]; fy=[]; if isfield(GeometryCalib,'fx_fy') fx=GeometryCalib.fx_fy(1); fy=GeometryCalib.fx_fy(2); end Cx_Cy=[0 0];%default if isfield(GeometryCalib,'Cx_Cy') Cx_Cy=GeometryCalib.Cx_Cy; end kc=0; if isfield(GeometryCalib,'kc') kc=GeometryCalib.kc; %* GeometryCalib.focal*GeometryCalib.focal; end set(handles.fx,'String',num2str(fx,5)) set(handles.fy,'String',num2str(fy,5)) set(handles.Cx,'String',num2str(Cx_Cy(1),'%1.1f')) set(handles.Cy,'String',num2str(Cx_Cy(2),'%1.1f')) set(handles.kc,'String',num2str(kc,'%1.4f')) %------------------------------------------------------------------------ % ---display calibration extrinsic parameters function display_extrinsic(GeometryCalib,handles) %------------------------------------------------------------------------ set(handles.Tx,'String',num2str(GeometryCalib.Tx_Ty_Tz(1),4)) set(handles.Ty,'String',num2str(GeometryCalib.Tx_Ty_Tz(2),4)) set(handles.Tz,'String',num2str(GeometryCalib.Tx_Ty_Tz(3),4)) set(handles.Phi,'String',num2str(GeometryCalib.omc(1),4)) set(handles.Theta,'String',num2str(GeometryCalib.omc(2),4)) set(handles.Psi,'String',num2str(GeometryCalib.omc(3),4)) %------------------------------------------------------------------------ % --- executes when user attempts to close geometry_calib. function geometry_calib_CloseRequestFcn(hObject, eventdata, handles) %------------------------------------------------------------------------ delete(hObject); % closes the figure %------------------------------------------------------------------------ % --- executes on button press in PLOT. %------------------------------------------------------------------------ function PLOT_Callback(hObject, eventdata, handles) huvmat=findobj(allchild(0),'Name','uvmat');%find the current uvmat interface handle hhuvmat=guidata(huvmat); %handles of GUI elements in uvmat h_menu_coord=findobj(huvmat,'Tag','TransformName'); menu=get(h_menu_coord,'String'); choice=get(h_menu_coord,'Value'); option=''; if iscell(menu) option=menu{choice}; end Coord=get(handles.ListCoord,'Data'); if ~isempty(Coord) if isequal(option,'phys') Coord_plot=Coord(:,1:3); elseif isempty(option);%optionoption,'px')||isequal(option,'') Coord_plot=Coord(:,4:5); else msgbox_uvmat('ERROR','the choice in menu_coord of uvmat must be blank or phys ') end end set(0,'CurrentFigure',huvmat) set(huvmat,'CurrentAxes',hhuvmat.PlotAxes) hh=findobj('Tag','calib_points'); if ~isempty(Coord) && isempty(hh) hh=line(Coord_plot(:,1),Coord_plot(:,2),'Color','m','Tag','calib_points','LineStyle','none','Marker','+','MarkerSize',10); elseif isempty(Coord)%empty list of points, suppress the plot delete(hh) else set(hh,'XData',Coord_plot(:,1)) set(hh,'YData',Coord_plot(:,2)) end pause(.1) figure(handles.geometry_calib) %------------------------------------------------------------------------ % --- Executes on button press in Copy: display Coord on the Matlab work space %------------------------------------------------------------------------ function Copy_Callback(hObject, eventdata, handles) global Coord evalin('base','global Coord')%make CurData global in the workspace Coord=get(handles.ListCoord,'Data'); display('coordinates of calibration points (phys,px,marker) :') evalin('base','Coord') %display CurData in the workspace commandwindow; %brings the Matlab command window to the front %------------------------------------------------------------------------ % --- Executes when selected cell(s) is changed in ListCoord. %------------------------------------------------------------------------ function ListCoord_CellSelectionCallback(hObject, eventdata, handles) if ~isempty(eventdata.Indices) iline=eventdata.Indices(1);% selected line number set(handles.CoordLine,'String',num2str(iline)) Data=get(handles.ListCoord,'Data'); update_calib_marker(Data(iline,:)) end %------------------------------------------------------------------------ % --- Executes when entered data in editable cell(s) in ListCoord. %------------------------------------------------------------------------ function ListCoord_CellEditCallback(hObject, eventdata, handles) Input=str2num(eventdata.EditData);%pasted input Coord=get(handles.ListCoord,'Data'); iline=eventdata.Indices(1);% selected line number if size(Coord,1)=1 && iline<=size(Coord,1) set(handles.CoordLine,'String',num2str(iline)) update_calib_marker(Coord(iline,1:5))% show the point corresponding to the selected line end else set(handles.APPLY,'BackgroundColor',[1 0 1])% paint APPLY in magenta to indicate that the table content has be modified if ismember(xx,[127 31])% delete, or downward Coord=get(handles.ListCoord,'Data'); iline=str2double(get(handles.CoordLine,'String')); if isequal(xx, 31) if isequal(iline,size(Coord,1))% arrow downward Coord=[Coord;zeros(1,size(Coord,2))]; end else Coord(iline,:)=[];% suppress the current line end set(handles.ListCoord,'Data',Coord); end end %------------------------------------------------------------------------ % --- update the plot of calibration points %------------------------------------------------------------------------ % draw a circle around the point defined by the input coordinates Coord as given by line in the table Listcoord function update_calib_marker(Coord) %% read config on uvmat huvmat=findobj(allchild(0),'Name','uvmat');%find the current uvmat interface handle hhuvmat=guidata(huvmat); hhh=findobj(hhuvmat.PlotAxes,'Tag','calib_marker'); if numel(Coord)<5 if ~isempty(hhh) delete(hhh)%delete the circle marker in case of no valid input end return end menu=get(hhuvmat.TransformName,'String'); choice=get(hhuvmat.TransformName,'Value'); option=''; if iscell(menu) option=menu{choice}; end %% read appropriate coordinates (px or phys) in the table ListCoord if isequal(option,'phys') % use phys coord XCoord=Coord(1); YCoord=Coord(2); elseif isempty(option)% use coord in pixels XCoord=Coord(4); YCoord=Coord(5); else msgbox_uvmat('ERROR','the choice in menu_coord of uvmat must be blank or phys ') return end %% adjust the plot limits if needed xlim=get(hhuvmat.PlotAxes,'XLim'); ylim=get(hhuvmat.PlotAxes,'YLim'); ind_range=max(abs(xlim(2)-xlim(1)),abs(ylim(end)-ylim(1)))/25;%defines the size of the circle marker check_xlim=0; if XCoord>xlim(2) xlim=xlim+XCoord-xlim(2)+ind_range;% translate plot limit check_xlim=1; elseif XCoordylim(2) ylim=ylim+YCoord-ylim(2)+ind_range;% translate plot limit check_ylim=1; elseif YCoord