Changeset 109 for trunk/src/geometry_calib.m
 Timestamp:
 Oct 4, 2010, 10:39:11 PM (10 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/src/geometry_calib.m
r108 r109 3 3 % function varargout = geometry_calib(varargin) 4 4 % 5 %A%AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AA5 %A%AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 6 6 % Copyright Joel Sommeria, 2008, LEGI / CNRSUJFINPG, sommeria@coriolislegi.org. 7 7 %AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA … … 43 43 % Edit the above text to modify the response to help geometry_calib 44 44 45 % Last Modified by GUIDE v2.5 25Mar2010 19:10:0545 % Last Modified by GUIDE v2.5 03Oct2010 09:34:12 46 46 47 47 % Begin initialization code  DO NOT edit … … 87 87 set(hObject,'Position',pos_gui); 88 88 end 89 90 %set menu of calibration options 91 %set(handles.calib_type,'String',{'rescale';'linear';'perspective';'normal';'tsai';'bouguet';'extrinsic'}) 92 set(handles.calib_type,'String',{'rescale';'linear';'quadr';'3D_linear';'3D_quadr';'3D_extrinsic'}) 89 93 inputxml=''; 90 94 if 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 95 struct.XmlInputfile=inputfile; 96 set(hObject,'UserData',struct) 97 [Pathsub,RootFile,field_count,str2,str_a,str_b,ext,nom_type,subdir]=name2display(inputfile); 98 inputxml=[fullfile(Pathsub,RootFile) '.xml']; 99 99 end 100 100 set(handles.ListCoord,'String',{'......'}) … … 102 102 loadfile(handles,inputxml)% load the point coordiantes existing in the xml file 103 103 end 104 105 104 set(handles.ListCoord,'KeyPressFcn',{@key_press_fcn,handles})%set keyboard action function 106 105 … … 113 112 varargout{1} = handles.output; 114 113 varargout{2}=handles; 115 116 %117 function Phi_Callback(hObject, eventdata, handles)118 119 %120 %read input xml file and update the edit boxes121 function loadfile(handles,fileinput)122 %123 %read the input xml file124 t=xmltree(fileinput);125 s=convert(t);%convert to matlab structure126 %read data currently displayed on the interface127 PointCoord=[];128 Coord_cell=get(handles.ListCoord,'String');129 data=read_geometry_calib(Coord_cell);130 %data=read_geometry_calib(handles);131 Coord=[]; %default132 if isfield(data,'Coord')133 Coord=data.Coord;134 end135 TabChar_0=get(handles.ListCoord,'String');136 nbcoord_0=size(TabChar_0,1);137 if isequal(get(handles.edit_append,'Value'),2) %edit mode A REVOIR138 val=get(handles.ListCoord,'Value')1;139 else140 val=length(TabChar_0);141 end142 nbcoord=0;143 144 %case of calibration (ImaDoc) input file145 % hcalib=get(handles.calib_type,'parent');146 CalibData=get(handles.geometry_calib,'UserData');147 CalibData.XmlInput=fileinput;148 if isfield(s,'Heading')149 CalibData.Heading=s.Heading;150 end151 152 set(handles.geometry_calib,'UserData',CalibData);%store the heading in the interface 'UserData'153 if isfield(s,'GeometryCalib')154 Calib=s.GeometryCalib;155 if isfield(Calib,'CalibrationType')156 CalibrationType=Calib.CalibrationType;157 switch CalibrationType158 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 end165 end166 if isfield(Calib,'SourceCalib')167 if isfield(Calib.SourceCalib,'PointCoord')168 PointCoord=Calib.SourceCalib.PointCoord;169 end170 end171 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 x176 Coord(i+val,1:3)=line(1:3)/10;%phys x177 end178 else179 for i=1:length(PointCoord)180 line=str2num(PointCoord{i});181 Coord(i,4:5)=line(4:5);%px x182 Coord(i,1:3)=line(1:3);%phys x183 end184 end185 end186 %case of xml files of points187 if isfield(s,'Coord')188 PointCoord=s.Coord;189 nbcoord=length(PointCoord);190 %case of image coordinates191 if isfield(s,'CoordType')& isequal(s.CoordType,'px')192 for i=1:nbcoord193 line=str2num(PointCoord{i});194 Coord(i+val,4:5)=line(1:2);195 end196 %case of physical coordinates197 else198 for i=1:nbcoord199 line=str2num(PointCoord{i});200 Coord(i+val,1:3)=line(1:3);201 nbcolumn=size(Coord,2);202 if nbcolumn<5203 Coord(i+val,nbcolumn+1:5)=zeros(1,5nbcolumn);204 end205 end206 end207 end208 CoordCell={};209 for iline=1:size(Coord,1)210 for j=1:5211 CoordCell{iline,j}=num2str(Coord(iline,j),4);212 end213 end214 % CoordCell=[CoordCell;{' ',' ',' ',' ',' '}];215 Tabchar=cell2tab(CoordCell,'  ');%transform cells into table ready for display216 Tabchar=[Tabchar;{'......'}];217 set(handles.ListCoord,'Value',1)218 set(handles.ListCoord,'String',Tabchar)219 MenuPlot_Callback(handles.geometry_calib, [], handles)220 if isempty(Coord)221 set(handles.edit_append,'Value',1)222 set(handles.edit_append,'BackgroundColor',[1 1 0])223 else224 set(handles.edit_append,'Value',0)225 set(handles.edit_append,'BackgroundColor',[0.7 0.7 0.7])226 end227 114 % 228 115 % … … 252 139 function APPLY_Callback(hObject, eventdata, handles) 253 140 % 254 calib_cell=get(handles.calib_type,'String'); 255 val=get(handles.calib_type,'Value'); 256 calib_type=calib_cell{val}; 141 142 %read the current calibration points 257 143 Coord_cell=get(handles.ListCoord,'String'); 258 144 Object=read_geometry_calib(Coord_cell); 259 X=Object.Coord(:,1); 260 Y=Object.Coord(:,2); 261 Z=Object.Coord(:,3); 262 if isequal(calib_type,'rescale') 263 GeometryCalib=calib_rescale(Object.Coord); 264 Z=0;%Z not taken into account 265 elseif isequal(calib_type,'linear') 266 GeometryCalib=calib_linear(Object.Coord); 267 Z=0; %Z not taken into account 268 elseif isequal(calib_type,'tsai_cpp') 269 GeometryCalib=calib_tsai(Object.Coord); 270 elseif isequal(calib_type,'tsai_matlab') 271 GeometryCalib=calib_tsai2(Object.Coord); 272 end 145 Coord=Object.Coord; 146 147 % apply the calibration, whose type is selected in handles.calib_type 148 if ~isempty(Coord) 149 calib_cell=get(handles.calib_type,'String'); 150 val=get(handles.calib_type,'Value'); 151 GeometryCalib=feval(['calib_' calib_cell{val}],Coord,handles); 152 else 153 msgbox_uvmat('ERROR','No calibration points, abort') 154 return 155 end 273 156 274 157 %check error … … 298 181 Calib.R=GeometryCalib.R; 299 182 end 300 x_ima=Object.Coord(:,4); 301 y_ima=Object.Coord(:,5); 302 [Xpoints,Ypoints]=px_XYZ(Calib,X,Y,Z); 303 GeometryCalib.ErrorRms(1)=sqrt(mean((Xpointsx_ima).*(Xpointsx_ima))); 304 [GeometryCalib.ErrorMax(1),index(1)]=max(abs(Xpointsx_ima)); 305 GeometryCalib.ErrorRms(2)=sqrt(mean((Ypointsy_ima).*(Ypointsy_ima))); 306 [GeometryCalib.ErrorMax(2),index(2)]=max(abs(Ypointsy_ima)); 307 [EM,ind_dim]=max(GeometryCalib.ErrorMax); 308 index=index(ind_dim); 309 183 if ~isempty(Coord) 184 X=Coord(:,1); 185 Y=Coord(:,2); 186 Z=Coord(:,3); 187 x_ima=Coord(:,4); 188 y_ima=Coord(:,5); 189 [Xpoints,Ypoints]=px_XYZ(Calib,X,Y,Z); 190 GeometryCalib.ErrorRms(1)=sqrt(mean((Xpointsx_ima).*(Xpointsx_ima))); 191 [GeometryCalib.ErrorMax(1),index(1)]=max(abs(Xpointsx_ima)); 192 GeometryCalib.ErrorRms(2)=sqrt(mean((Ypointsy_ima).*(Ypointsy_ima))); 193 [GeometryCalib.ErrorMax(2),index(2)]=max(abs(Ypointsy_ima)); 194 [EM,ind_dim]=max(GeometryCalib.ErrorMax); 195 index=index(ind_dim); 196 if isequal(max(Z),min(Z))%Z constant 197 GeometryCalib.NbSlice=1; 198 GeometryCalib.SliceCoord=[0 0 Z(1)]; 199 end 200 end 310 201 unitlist=get(handles.CoordUnit,'String'); 311 202 unit=unitlist{get(handles.CoordUnit,'value')}; 312 203 GeometryCalib.CoordUnit=unit; 313 GeometryCalib.SourceCalib.PointCoord=Object.Coord; 204 GeometryCalib.SourceCalib.PointCoord=Coord; 205 206 %display calibration intrinsic parameters 207 k=0; 208 if isfield(GeometryCalib,'kappa1') 209 k=GeometryCalib.kappa1 * GeometryCalib.focal*GeometryCalib.focal; 210 end 211 sx=1;dpx_dpy=[1 1]; 212 if isfield(GeometryCalib,'sx') 213 sx=GeometryCalib.sx; 214 end 215 if isfield(GeometryCalib,'dpx_dpy') 216 dpx_dpy=GeometryCalib.dpx_dpy; 217 end 218 Cx_Cy=[0 0]; 219 if isfield(GeometryCalib,'Cx_Cy') 220 Cx_Cy=GeometryCalib.Cx_Cy; 221 end 222 f1=GeometryCalib.focal*sx/dpx_dpy(1); 223 f2=GeometryCalib.focal/dpx_dpy(2); 224 Cx=Cx_Cy(1); 225 Cy=Cx_Cy(2); 226 set(handles.fx,'String',num2str(f1,4)) 227 set(handles.fy,'String',num2str(f2,4)) 228 set(handles.k,'String',num2str(k,'%1.4f')) 229 set(handles.Cx,'String',num2str(Cx,'%1.1f')) 230 set(handles.Cy,'String',num2str(Cy,'%1.1f')) 231 % Display extrinsinc parameters 232 set(handles.Tx,'String',num2str(GeometryCalib.Tx_Ty_Tz(1),4)) 233 set(handles.Ty,'String',num2str(GeometryCalib.Tx_Ty_Tz(2),4)) 234 set(handles.Tz,'String',num2str(GeometryCalib.Tx_Ty_Tz(3),4)) 235 236 % indicate the plane of the calibration grid if defined 314 237 huvmat=findobj(allchild(0),'Name','uvmat'); 315 238 hhuvmat=guidata(huvmat);%handles of elements in the GUI uvmat … … 337 260 hhh=findobj(hhuvmat.axes3,'Tag','calib_marker');% delete calib points and markers 338 261 if ~isempty(hhh) 339 delete(hhh); 262 delete(hhh); 340 263 end 341 264 hhh=findobj(hhuvmat.axes3,'Tag','calib_points'); … … 346 269 set(hhuvmat.FixedLimits,'BackgroundColor',[0.7 0.7 0.7]) 347 270 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) 271 % if ip==0 272 MenuPlot_Callback(hObject, eventdata, handles) 273 set(handles.ListCoord,'Value',index)% indicate in the list the point with max deviation (possible mistake) 274 ListCoord_Callback(hObject, eventdata, handles) 275 % end 351 276 figure(handles.geometry_calib) 352 277 end … … 355 280 %  Executes on button press in calibrate_lin. 356 281 function REPLICATE_Callback(hObject, eventdata, handles) 282 % %%%%%%Todo: correct on the model of APPLY_Callback%%%%%%%%%% 357 283 % 358 284 calib_cell=get(handles.calib_type,'String'); 359 285 val=get(handles.calib_type,'Value'); 360 calib_type=calib_cell{val};361 286 Coord_cell=get(handles.ListCoord,'String'); 362 287 Object=read_geometry_calib(Coord_cell); 363 364 if isequal(calib_type,'rescale') 365 GeometryCalib=calib_rescale(Object.Coord); 366 elseif isequal(calib_type,'linear') 367 GeometryCalib=calib_linear(Object.Coord); 368 elseif isequal(calib_type,'tsai') 369 GeometryCalib=calib_tsai(Object.Coord); 370 end 288 GeometryCalib=feval(calib_cell{val},Object.Coord); 289 371 290 % %record image source 372 291 GeometryCalib.SourceCalib.PointCoord=Object.Coord; … … 424 343 % 425 344 % determine the parameters for a calibration by an affine function (rescaling and offset, no rotation) 426 function GeometryCalib=calib_rescale(Coord) 427 % 428 345 function GeometryCalib=calib_rescale(Coord,handles) 346 % 429 347 X=Coord(:,1); 430 Y=Coord(:,2); 348 Y=Coord(:,2);% Z not used 431 349 x_ima=Coord(:,4); 432 350 y_ima=Coord(:,5); … … 436 354 T_y=py(2); 437 355 GeometryCalib.CalibrationType='rescale'; 438 GeometryCalib.focal=1; 356 GeometryCalib.focal=1;%default 357 GeometryCalib.sx=px(1)/py(1); 439 358 GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function 440 GeometryCalib.Tx_Ty_Tz=[T_x T_y 1]; 441 GeometryCalib.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((Xpointsx_ima).*(Xpointsx_ima))); 457 % GeometryCalib.ErrorMax(1)=max(abs(Xpointsx_ima)); 458 % GeometryCalib.ErrorRms(2)=sqrt(mean((Ypointsy_ima).*(Ypointsy_ima))); 459 % GeometryCalib.ErrorMax(2)=max(abs(Ypointsy_ima)); 359 GeometryCalib.Tx_Ty_Tz=[T_x T_y GeometryCalib.focal/py(1)]; 360 GeometryCalib.R=[1,0,0;0,1,0;0,0,0]; 460 361 461 362 % 462 363 % determine the parameters for a calibration by a linear transform matrix (rescale and rotation) 463 function GeometryCalib=calib_linear(Coord )364 function GeometryCalib=calib_linear(Coord,handles) 464 365 % 465 366 X=Coord(:,1); 466 Y=Coord(:,2); 367 Y=Coord(:,2);% Z not used 467 368 x_ima=Coord(:,4); 468 369 y_ima=Coord(:,5); 469 370 XY_mat=[ones(size(X)) X Y]; 470 371 a_X1=XY_mat\x_ima; %transformation matrix for X 471 x1=XY_mat*a_X1;%reconstruction472 err_X1=max(abs(x1x_ima));%error372 % x1=XY_mat*a_X1;%reconstruction 373 % err_X1=max(abs(x1x_ima));%error 473 374 a_Y1=XY_mat\y_ima;%transformation matrix for X 474 y1=XY_mat*a_Y1; 475 err_Y1=max(abs(y1y_ima));%error 476 T_x=a_X1(1); 477 T_y=a_Y1(1); 375 % y1=XY_mat*a_Y1; 376 % err_Y1=max(abs(y1y_ima));%error 478 377 GeometryCalib.CalibrationType='linear'; 479 378 GeometryCalib.focal=1; 480 379 GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function 481 GeometryCalib.Tx_Ty_Tz=[T_x T_y 1]; 482 GeometryCalib.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((x1x_ima).*(x1x_ima))); 486 % GeometryCalib.ErrorMax(1)=max(abs(x1x_ima)); 487 % GeometryCalib.ErrorRms(2)=sqrt(mean((y1y_ima).*(y1y_ima))); 488 % GeometryCalib.ErrorMax(2)=max(abs(y1y_ima)); 489 490 % 491 function GeometryCalib=calib_tsai2(Coord) 380 R=[a_X1(2),a_X1(3),0;a_Y1(2),a_Y1(3),0;0,0,0]; 381 norm=det(R); 382 GeometryCalib.Tx_Ty_Tz=[a_X1(1) a_Y1(1) norm]; 383 GeometryCalib.R=R/norm; 384 385 % 386 % determine the tsai parameters for a view normal to the grid plane 387 function GeometryCalib=calib_normal(Coord,handles) 388 % 389 Calib.f1=str2num(get(handles.fx,'String')); 390 Calib.f2=str2num(get(handles.fy,'String')); 391 Calib.k=str2num(get(handles.k,'String')); 392 Calib.Cx=str2num(get(handles.Cx,'String')); 393 Calib.Cy=str2num(get(handles.Cy,'String')); 394 %default 395 if isempty(Calib.f1) 396 Calib.f1=25/0.012; 397 end 398 if isempty(Calib.f2) 399 Calib.f2=25/0.012; 400 end 401 if isempty(Calib.k) 402 Calib.k=0; 403 end 404 if isempty(Calib.Cx)isempty(Calib.Cy) 405 huvmat=findobj(allchild(0),'Tag','uvmat'); 406 hhuvmat=guidata(huvmat); 407 Calib.Cx=str2num(get(hhuvmat.npx,'String'))/2; 408 Calib.Cx=str2num(get(hhuvmat.npy,'String'))/2; 409 end 410 %tsai parameters 411 Calib.dpx=0.012;%arbitrary 412 Calib.dpy=0.012; 413 Calib.sx=Calib.f1*Calib.dpx/(Calib.f2*Calib.dpy); 414 Calib.f=Calib.f2*Calib.dpy; 415 Calib.kappa1=Calib.k/(Calib.f*Calib.f); 416 417 %initial guess 418 X=Coord(:,1); 419 Y=Coord(:,2); 420 Zmean=mean(Coord(:,3)); 421 x_ima=Coord(:,4)Calib.Cx; 422 y_ima=Coord(:,5)Calib.Cy; 423 XY_mat=[ones(size(X)) X Y]; 424 a_X1=XY_mat\x_ima; %transformation matrix for X 425 a_Y1=XY_mat\y_ima;%transformation matrix for Y 426 R=[a_X1(2),a_X1(3),0;a_Y1(2),a_Y1(3),0;0,0,1];% rotation+ z axis reversal (upward) 427 norm=sqrt(det(R)); 428 calib_param(1)=0;% quadratic distortion 429 calib_param(2)=a_X1(1); 430 calib_param(3)=a_Y1(1); 431 calib_param(4)=Calib.f/(norm*Calib.dpx)R(3,3)*Zmean; 432 calib_param(5)=angle(a_X1(2)+i*a_X1(3)); 433 display(['initial guess=' num2str(calib_param)]) 434 435 %optimise the parameters: minimisation of error 436 calib_param = fminsearch(@(calib_param) error_calib(calib_param,Calib,Coord),calib_param) 437 438 GeometryCalib.CalibrationType='tsai_normal'; 439 GeometryCalib.focal=Calib.f; 440 GeometryCalib.dpx_dpy=[Calib.dpx Calib.dpy]; 441 GeometryCalib.Cx_Cy=[Calib.Cx Calib.Cy]; 442 GeometryCalib.sx=Calib.sx; 443 GeometryCalib.kappa1=calib_param(1); 444 GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function 445 GeometryCalib.Tx_Ty_Tz=[calib_param(2) calib_param(3) calib_param(4)]; 446 alpha=calib_param(5); 447 GeometryCalib.R=[cos(alpha) sin(alpha) 0;sin(alpha) cos(alpha) 0;0 0 1]; 448 449 % 450 function GeometryCalib=calib_3D_linear(Coord,handles) 492 451 % 493 452 path_uvmat=which('uvmat');% check the path detected for source file uvmat … … 495 454 huvmat=findobj(allchild(0),'Tag','uvmat'); 496 455 hhuvmat=guidata(huvmat); 456 x_1=Coord(:,4:5)' 457 X_1=Coord(:,1:3)' 458 n_ima=1; 459 % check_cond=0; 460 461 nx=str2num(get(hhuvmat.npx,'String')); 462 ny=str2num(get(hhuvmat.npy,'String')); 463 464 est_dist=[0;0;0;0;0]; 465 est_aspect_ratio=0; 466 est_fc=[1;1]; 467 %fc=[25;25]/0.012; 468 center_optim=0; 469 run(fullfile(path_UVMAT,'toolbox_calib','go_calib_optim')); 470 471 GeometryCalib.CalibrationType='3D_linear'; 472 GeometryCalib.focal=fc(2); 473 GeometryCalib.dpx_dpy=[1 1]; 474 GeometryCalib.Cx_Cy=cc'; 475 GeometryCalib.sx=fc(1)/fc(2); 476 GeometryCalib.kappa1=kc(1)/fc(2)^2; 477 GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function 478 GeometryCalib.Tx_Ty_Tz=Tc_1'; 479 GeometryCalib.R=Rc_1; 480 481 % 482 function GeometryCalib=calib_3D_quadr(Coord,handles) 483 % 484 485 path_uvmat=which('uvmat');% check the path detected for source file uvmat 486 path_UVMAT=fileparts(path_uvmat); %path to UVMAT 487 huvmat=findobj(allchild(0),'Tag','uvmat'); 488 hhuvmat=guidata(huvmat); 489 % check_cond=0; 490 coord_files=get(handles.coord_files,'String'); 491 if ischar(coord_files) 492 coord_files={coord_files}; 493 end 494 if isempty(coord_files{1})  isequal(coord_files,{''}) 495 coord_files={}; 496 end 497 498 %retrieve the calibration points stored in the files listed in the popup list coord_files 499 x_1=Coord(:,4:5)'; 500 X_1=Coord(:,1:3)'; 501 n_ima=numel(coord_files)+1; 502 if ~isempty(coord_files) 503 msgbox_uvmat('CONFIRMATION',['The xy coordinates of the calibration points in ' num2str(n_ima) ' planes will be used']) 504 for ifile=1:numel(coord_files) 505 t=xmltree(coord_files{ifile}); 506 s=convert(t);%convert to matlab structure 507 if isfield(s,'GeometryCalib') 508 if isfield(s.GeometryCalib,'SourceCalib') 509 if isfield(s.GeometryCalib.SourceCalib,'PointCoord') 510 PointCoord=s.GeometryCalib.SourceCalib.PointCoord; 511 Coord_file=[]; 512 for i=1:length(PointCoord) 513 line=str2num(PointCoord{i}); 514 Coord_file(i,4:5)=line(4:5);%px x 515 Coord_file(i,1:3)=line(1:3);%phys x 516 end 517 eval(['x_' num2str(ifile+1) '=Coord_file(:,4:5)'';']); 518 eval(['X_' num2str(ifile+1) '=Coord_file(:,1:3)'';']); 519 end 520 end 521 end 522 end 523 end 524 525 n_ima=numel(coord_files)+1; 526 whos 527 nx=str2num(get(hhuvmat.npx,'String')); 528 ny=str2num(get(hhuvmat.npy,'String')); 529 530 est_dist=[1;0;0;0;0]; 531 est_aspect_ratio=1; 532 %est_fc=[0;0]; 533 %fc=[25;25]/0.012; 534 center_optim=0; 535 run(fullfile(path_UVMAT,'toolbox_calib','go_calib_optim')); 536 537 GeometryCalib.CalibrationType='3D_quadr'; 538 GeometryCalib.focal=fc(2); 539 GeometryCalib.dpx_dpy=[1 1]; 540 GeometryCalib.Cx_Cy=cc'; 541 GeometryCalib.sx=fc(1)/fc(2); 542 GeometryCalib.kappa1=kc(1)/fc(2)^2; 543 GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function 544 GeometryCalib.Tx_Ty_Tz=Tc_1'; 545 GeometryCalib.R=Rc_1; 546 GeometryCalib.R(3,1:3)=GeometryCalib.R(3,1:3);%inversion for z upward 547 GeometryCalib.ErrorRMS=[]; 548 GeometryCalib.ErrorMax=[]; 549 550 551 % 552 function GeometryCalib=calib_3D_extrinsic(Coord,handles) 553 % 554 path_uvmat=which('geometry_calib');% check the path detected for source file uvmat 555 path_UVMAT=fileparts(path_uvmat); %path to UVMAT 497 556 x_1=Coord(:,4:5)'; 498 557 X_1=Coord(:,1:3)'; 499 558 n_ima=1; 500 % check_cond=0; 501 502 503 nx=str2num(get(hhuvmat.npx,'String')); 504 ny=str2num(get(hhuvmat.npy,'String')); 505 506 507 % est_kc=[1;0;0;0;0]; 508 est_fc=0; 509 est_dist=[1;0;0;0;0]; 510 run(fullfile(path_UVMAT,'toolbox_calib','go_calib_optim')); 511 512 GeometryCalib.CalibrationType='tsai_matlab'; 513 GeometryCalib.focal=f(2); 559 Calib.f1=str2num(get(handles.fx,'String')); 560 Calib.f2=str2num(get(handles.fy,'String')); 561 Calib.k=str2num(get(handles.k,'String')); 562 Calib.Cx=str2num(get(handles.Cx,'String')); 563 Calib.Cy=str2num(get(handles.Cy,'String')); 564 565 fct_path=fullfile(path_UVMAT,'toolbox_calib'); 566 addpath(fct_path) 567 % [omc1,Tc1,Rc1,H,x,ex,JJ] = compute_extrinsic(x_1,X_1,... 568 % [Calib.f Calib.f*Calib.sx]',... 569 % [Calib.Cx Calib.Cy]',... 570 % [Calib.kappa1*Calib.f^2 0 0 0 0]); 571 [omc1,Tc1,Rc1,H,x,ex,JJ] = compute_extrinsic(x_1,X_1,... 572 [Calib.f1 Calib.f2]',... 573 [Calib.Cx Calib.Cy]',... 574 [Calib.k 0 0 0 0]); 575 %get the euler angles ??? 576 rmpath(fct_path); 577 578 std(ex') 579 GeometryCalib.CalibrationType='3D_extrinsic'; 580 GeometryCalib.focal=Calib.f2; 514 581 GeometryCalib.dpx_dpy=[1 1]; 515 GeometryCalib.Cx_Cy= cc';516 GeometryCalib.sx= fc(1)/fc(2);517 GeometryCalib.kappa1= k(1)/f(2)^2;582 GeometryCalib.Cx_Cy=[Calib.Cx Calib.Cy]; 583 GeometryCalib.sx=Calib.f1/Calib.f2; 584 GeometryCalib.kappa1=Calib.k/(Calib.f2*Calib.f2); 518 585 GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function 519 GeometryCalib.Tx_Ty_Tz=Tc_1'; 520 GeometryCalib.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((Xpointsx_ima).*(Xpointsx_ima))); 540 % [GeometryCalib.ErrorMax(1),GeometryCalib.IndexMax(1)]=max(abs(Xpointsx_ima)); 541 % GeometryCalib.ErrorRms(2)=sqrt(mean((Ypointsy_ima).*(Ypointsy_ima))); 542 % [GeometryCalib.ErrorMax(2),GeometryCalib.IndexMax(2)]=max(abs(Ypointsy_ima)); 543 544 function GeometryCalib=calib_tsai(Coord) 586 GeometryCalib.Tx_Ty_Tz=Tc1'; 587 %inversion of z axis 588 GeometryCalib.R=Rc1; 589 %GeometryCalib.R(3,1:3)=GeometryCalib.R(3,1:3);%inversion for z upward 590 591 592 % 593 %function GeometryCalib=calib_tsai_heikkila(Coord) 594 % TEST: NOT IMPLEMENTED 595 % 596 % path_uvmat=which('uvmat');% check the path detected for source file uvmat 597 % path_UVMAT=fileparts(path_uvmat); %path to UVMAT 598 % path_calib=fullfile(path_UVMAT,'toolbox_calib_heikkila'); 599 % addpath(path_calib) 600 % npoints=size(Coord,1); 601 % Coord(:,1:3)=10*Coord(:,1:3); 602 % Coord=[Coord zeros(npoints,2) ones(npoints,1)]; 603 % [par,pos,iter,res,er,C]=cacal('dalsa',Coord); 604 % GeometryCalib.CalibrationType='tsai'; 605 % GeometryCalib.focal=par(2); 606 607 608 % 609 function GeometryCalib=calib_tsai(Coord,handles)% old version using gauthier's bianry ccal_fo 545 610 % 546 611 %TSAI 547 % 'calibration_lin' provides a linear transform on coordinates,548 612 path_uvmat=which('uvmat');% check the path detected for source file uvmat 549 613 path_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 552 xmlfile=fullfile(path_UVMAT,'PARAM.xml'); 614 xmlfile=fullfile(path_UVMAT,'PARAM.xml');%name of the file containing names of binary executables 553 615 if exist(xmlfile,'file') 554 t=xmltree(xmlfile); 555 sparam=convert(t); 616 t=xmltree(xmlfile);% read the (xml) file containing names of binary executables 617 sparam=convert(t);% convert to matlab structure 556 618 end 557 619 if ~isfield(sparam,'GeometryCalibBin') … … 570 632 textcoord=num2str(Coord,4); 571 633 dlmwrite('t.txt',textcoord,''); 572 % ['!' Tsai_exe ' f 1 0 f2t.txt']573 634 % ['!' Tsai_exe ' fx 0 fy t.txt'] 635 eval(['!' Tsai_exe ' f t.txt > tsaicalib.log']); 574 636 if ~exist('calib.dat','file') 575 637 msgbox_uvmat('ERROR','no output from calibration program Tsai_exe: possibly too few points') … … 577 639 calibdat=dlmread('calib.dat'); 578 640 delete('calib.dat') 579 delete('t.txt')641 %delete('t.txt') 580 642 GeometryCalib.CalibrationType='tsai'; 581 643 GeometryCalib.focal=calibdat(10); … … 604 666 %EN DEDUIRE MATRICE R ?? 605 667 GeometryCalib.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((Xpointsx_ima).*(Xpointsx_ima))); 627 % GeometryCalib.ErrorMax(1)=max(abs(Xpointsx_ima)); 628 % GeometryCalib.ErrorRms(2)=sqrt(mean((Ypointsy_ima).*(Ypointsy_ima))); 629 % GeometryCalib.ErrorMax(2)=max(abs(Ypointsy_ima)); 630 631 632 % 633 %  Executes on button press in rotation. 634 function rotation_Callback(hObject, eventdata, handles) 635 % 636 angle_rot=(pi/180)*str2num(get(handles.Phi,'String')); 668 669 % 670 %  determine the rms of calibration error 671 function ErrorRms=error_calib(calib_param,Calib,Coord) 672 %calib_param: vector of free calibration parameters (to optimise) 673 %Calib: structure of the given calibration parameters 674 %Coord: list of phys coordinates (columns 13, and pixel coordinates (columns 45) 675 Calib.f=25; 676 Calib.dpx=0.012; 677 Calib.dpy=0.012; 678 Calib.sx=1; 679 Calib.Cx=512; 680 Calib.Cy=512; 681 Calib.kappa1=calib_param(1); 682 Calib.Tx=calib_param(2); 683 Calib.Ty=calib_param(3); 684 Calib.Tz=calib_param(4); 685 alpha=calib_param(5); 686 Calib.R=[cos(alpha) sin(alpha) 0;sin(alpha) cos(alpha) 0;0 0 1]; 687 688 X=Coord(:,1); 689 Y=Coord(:,2); 690 Z=Coord(:,3); 691 x_ima=Coord(:,4); 692 y_ima=Coord(:,5); 693 [Xpoints,Ypoints]=px_XYZ(Calib,X,Y,Z); 694 ErrorRms(1)=sqrt(mean((Xpointsx_ima).*(Xpointsx_ima))); 695 ErrorRms(2)=sqrt(mean((Ypointsy_ima).*(Ypointsy_ima))); 696 ErrorRms=mean(ErrorRms); 697 698 % 699 function XImage_Callback(hObject, eventdata, handles) 700 % 701 update_list(hObject, eventdata,handles) 702 703 % 704 function YImage_Callback(hObject, eventdata, handles) 705 % 706 update_list(hObject, eventdata,handles) 707 708 % 709 %  Executes on button press in STORE. 710 function STORE_Callback(hObject, eventdata, handles) 637 711 Coord_cell=get(handles.ListCoord,'String'); 638 data=read_geometry_calib(Coord_cell); 639 data.Coord(:,1)=cos(angle_rot)*data.Coord(:,1)+sin(angle_rot)*data.Coord(:,2); 640 data.Coord(:,1)=sin(angle_rot)*data.Coord(:,1)+cos(angle_rot)*data.Coord(:,2); 641 set(handles.XObject,'String',num2str(data.Coord(:,1),4)); 642 set(handles.YObject,'String',num2str(data.Coord(:,2),4)); 643 644 % 645 function XImage_Callback(hObject, eventdata, handles) 712 Object=read_geometry_calib(Coord_cell); 713 unitlist=get(handles.CoordUnit,'String'); 714 unit=unitlist{get(handles.CoordUnit,'value')}; 715 GeometryCalib.CoordUnit=unit; 716 GeometryCalib.SourceCalib.PointCoord=Object.Coord; 717 huvmat=findobj(allchild(0),'Name','uvmat'); 718 hhuvmat=guidata(huvmat);%handles of elements in the GUI uvmat 719 RootPath=''; 720 RootFile=''; 721 if ~isempty(hhuvmat.RootPath)& ~isempty(hhuvmat.RootFile) 722 testhandle=1; 723 RootPath=get(hhuvmat.RootPath,'String'); 724 RootFile=get(hhuvmat.RootFile,'String'); 725 filebase=fullfile(RootPath,RootFile); 726 while exist([filebase '.xml'],'file') 727 filebase=[filebase '~']; 728 end 729 outputfile=[filebase '.xml']; 730 update_imadoc(GeometryCalib,outputfile) 731 listfile=get(handles.coord_files,'string'); 732 if isequal(listfile,{''}) 733 listfile={outputfile}; 734 else 735 listfile=[listfile;{outputfile}]%update the list of coord files 736 end 737 set(handles.coord_files,'string',listfile); 738 end 739 set(handles.ListCoord,'Value',1)% refresh the display of coordinates 740 set(handles.ListCoord,'String',{'......'}) 741 742 % 743 %  Executes on button press in CLEAR. 744 function CLEAR_Callback(hObject, eventdata, handles) 745 % 746 set(handles.coord_files,'Value',1) 747 set(handles.coord_files,'String',{''}) 748 749 % 750 function XObject_Callback(hObject, eventdata, handles) 646 751 % 647 752 update_list(hObject, eventdata,handles) 648 753 649 754 % 650 function Y Image_Callback(hObject, eventdata, handles)755 function YObject_Callback(hObject, eventdata, handles) 651 756 % 652 757 update_list(hObject, eventdata,handles) 653 758 654 function XObject_Callback(hObject, eventdata, handles) 655 update_list(hObject, eventdata,handles) 656 657 function YObject_Callback(hObject, eventdata, handles) 658 update_list(hObject, eventdata,handles) 659 759 % 660 760 function ZObject_Callback(hObject, eventdata, handles) 761 % 661 762 update_list(hObject, eventdata,handles) 662 763 … … 748 849 % 749 850 choice=get(handles.edit_append,'Value'); 750 % if choice==1751 % Coord=get(handles.ListCoord,'String');752 % val=length(Coord);753 % if val>=1 & isequal(Coord{val},'')754 % val=val1; %do not take into account blank755 % end756 % Coord{val+1}='';757 % set(handles.ListCoord,'String',Coord)758 % set(handles.ListCoord,'Value',val+1)759 % end760 851 if choice 761 852 set(handles.edit_append,'BackgroundColor',[1 1 0]) … … 817 908 end 818 909 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=val1; %do not take into account blank827 % end828 % Coord{val+1}='';829 % set(handles.ListCoord,'String',Coord)830 % set(handles.ListCoord,'Value',val+1)831 832 %833 function MenuOpen_Callback(hObject, eventdata, handles)834 %835 %get the object file836 huvmat=findobj(allchild(0),'Name','uvmat');837 UvData=get(huvmat,'UserData');838 hchild=get(huvmat,'Children');839 hrootpath=findobj(hchild,'Tag','RootPath');840 oldfile=get(hrootpath,'String');841 if isempty(oldfile)842 oldfile='';843 end844 %[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);850 fileinput=[PathName FileName];%complete file name851 testblank=findstr(fileinput,' ');%look for blanks852 if ~isempty(testblank)853 msgbox_uvmat('ERROR','forbidden input file name or path: no blank character allowed')854 return855 end856 sizf=size(fileinput);857 if (~ischar(fileinput)~isequal(sizf(1),1)),return;end858 loadfile(handles,fileinput)859 860 861 910 % 862 911 function MenuPlot_Callback(hObject, eventdata, handles) … … 955 1004 grid_input=CalibData.grid;%retrieve the previously used grid 956 1005 end 957 [T,CalibData.grid]=create_grid(grid_input,'detect grid');%display the GUI create_grid, read the set of phys coordinates T 1006 [T,CalibData.grid]=create_grid(grid_input,'detect_grid');%display the GUI create_grid, read the set of phys coordinates T 1007 958 1008 set(handles.geometry_calib,'UserData',CalibData)%store the phys grid parameters for later use 959 1009 … … 988 1038 X=[CalibData.grid.x_0 CalibData.grid.x_1 CalibData.grid.x_0 CalibData.grid.x_1]';%corner absissa in the phys coordinates 989 1039 Y=[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(x1corners_X));%error 996 %a_Y1=XY_mat\corners_Y;%transformation matrix for X 997 %y1=XY_mat*a_Y1; 998 %err_Y1=max(abs(y1corners_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)+(1u)*corners_X(1); 1005 % x23=u*corners_X(3)+(1u)*corners_X(2); 1006 % x34=u*corners_X(4)+(1u)*corners_X(3); 1007 % x41=u*corners_X(1)+(1u)*corners_X(4); 1008 % y12=u*corners_Y(2)+(1u)*corners_Y(1); 1009 % y23=u*corners_Y(3)+(1u)*corners_Y(2); 1010 % y34=u*corners_Y(4)+(1u)*corners_Y(3); 1011 % y41=u*corners_Y(1)+(1u)*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') 1040 1017 1041 %calculate transform matrices: 1018 1042 % reference: http://alumni.media.mit.edu/~cwren/interpolator/ by Christopher R. Wren … … 1026 1050 C = [l(7:8)' 1]; 1027 1051 1028 GeometryCalib.CalibrationType='tsai'; %'linear';1052 GeometryCalib.CalibrationType='tsai'; 1029 1053 GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function 1030 1054 GeometryCalib.f=1; … … 1039 1063 GeometryCalib.Tz=1; 1040 1064 GeometryCalib.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 transformes1042 % 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 % end1048 % xphys12=u*X(2)+(1u)*X(1);1049 % xphys23=u*X(3)+(1u)*X(2);1050 % xphys34=u*X(4)+(1u)*X(3);1051 % xphys41=u*X(1)+(1u)*X(4);1052 % yphys12=u*Y(2)+(1u)*Y(1);1053 % yphys23=u*Y(3)+(1u)*Y(2);1054 % yphys34=u*Y(4)+(1u)*Y(3);1055 % yphys41=u*Y(1)+(1u)*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 % end1062 % 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 off1080 % plot(XY12(1,:),XY12(2,:),'ob')1081 % hold on1082 % 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 on1087 % plot(xphys23,yphys23,'ob')1088 % plot(xphys34,yphys34,'ob')1089 % plot(xphys41,yphys41,'ob')1090 1065 1091 1066 [Amod,Rangx,Rangy]=phys_Ima(Amin(min(A)),GeometryCalib,0); 1092 1067 1093 1068 Amod=double(Amod); 1094 figure(12) 1095 Amax=max(max(Amod));1096 image(Rangx,Rangy,uint8(255*Amod/Amax))1069 % figure(12) display corrected image 1070 % Amax=max(max(Amod)); 1071 % image(Rangx,Rangy,uint8(255*Amod/Amax)) 1097 1072 Dx=(Rangx(2)Rangx(1))/(npxy(2)1); %x mesh in real space 1098 1073 Dy=(Rangy(2)Rangy(1))/(npxy(1)1); %y mesh in real space … … 1103 1078 i0=1+round((T(ipoint,1)Rangx(1))/Dx);%round(Xpx(ipoint)); 1104 1079 j0=1+round((T(ipoint,2)Rangy(1))/Dy);%round(Xpx(ipoint)); 1105 Asub=Amod(j0ind_range_y:j0+ind_range_y,i0ind_range_x:i0+ind_range_x); 1080 j0min=max(j0ind_range_y,1); 1081 j0max=min(j0+ind_range_y,size(Amod,1)); 1082 i0min=max(i0ind_range_x,1); 1083 i0max=min(i0+ind_range_x,size(Amod,2)); 1084 Asub=Amod(j0min:j0max,i0min:i0max); 1106 1085 x_profile=sum(Asub,1); 1107 1086 y_profile=sum(Asub,2); … … 1111 1090 x_shift=0; 1112 1091 y_shift=0; 1113 if ind_x_max+2<=2*ind_range_x+1 && ind_x_max2>=1 1092 %if ind_x_max+2<=2*ind_range_x+1 && ind_x_max2>=1 1093 if ind_x_max+2<=numel(x_profile) && ind_x_max2>=1 1114 1094 Atop=x_profile(ind_x_max2:ind_x_max+2); 1115 1095 x_shift=sum(Atop.*[2 1 0 1 2])/sum(Atop); 1116 1096 end 1117 if ind_ y_max+2<=2*ind_range_y+1&& ind_y_max2>=11097 if ind_x_max+2<=numel(y_profile) && ind_y_max2>=1 1118 1098 Atop=y_profile(ind_y_max2:ind_y_max+2); 1119 1099 y_shift=sum(Atop.*[2 1 0 1 2]')/sum(Atop); 1120 1100 end 1121 Delta(ipoint,1)=(x_shift+ind_x_max ind_range_x1)*Dx;%shift from the initial guess1122 Delta(ipoint,2)=(y_shift+ind_y_max ind_range_y1)*Dy;1101 Delta(ipoint,1)=(x_shift+ind_x_max+i0mini01)*Dx;%shift from the initial guess 1102 Delta(ipoint,2)=(y_shift+ind_y_max+j0minj01)*Dy; 1123 1103 end 1124 1104 Tmod=T(:,(1:2))+Delta; … … 1127 1107 Coord{ipoint,1}=num2str(T(ipoint,1),4);%display coordiantes with 4 digits 1128 1108 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 1109 Coord{ipoint,3}=num2str(T(ipoint,3),4);%display coordiantes with 4 digits; 1110 Coord{ipoint,4}=num2str(Xpx(ipoint),4);%display coordiantes with 4 digits 1111 Coord{ipoint,5}=num2str(Ypx(ipoint),4);%display coordiantes with 4 digits 1132 1112 end 1133 1113 Tabchar=cell2tab(Coord(end:1:1,:),'  '); … … 1209 1189 set(handles.ListCoord,'Value',1) 1210 1190 set(handles.ListCoord,'String',Tabchar) 1191 1192 1193 % % 1194 % %  Executes on button press in rotation. 1195 % function rotation_Callback(hObject, eventdata, handles) 1196 % % 1197 % angle_rot=(pi/180)*str2num(get(handles.Phi,'String')); 1198 % Coord_cell=get(handles.ListCoord,'String'); 1199 % data=read_geometry_calib(Coord_cell); 1200 % data.Coord(:,1)=cos(angle_rot)*data.Coord(:,1)+sin(angle_rot)*data.Coord(:,2); 1201 % data.Coord(:,1)=sin(angle_rot)*data.Coord(:,1)+cos(angle_rot)*data.Coord(:,2); 1202 % set(handles.XObject,'String',num2str(data.Coord(:,1),4)); 1203 % set(handles.YObject,'String',num2str(data.Coord(:,2),4)); 1204 1211 1205 1212 1206 % … … 1287 1281 Zphys=Calib.SliceCoord(Zindex,3);%GENERALISER AUX CAS AVEC ANGLE 1288 1282 else 1289 % if exist('Z','var') 1290 % Zphys=Z; 1291 % else 1292 Zphys=0; 1293 % end 1283 Zphys=0; 1294 1284 end 1295 1285 if ~exist('X','var')~exist('Y','var') … … 1329 1319 Yphys=(A21.*Xu+A22.*Yu+Y0)./denom; 1330 1320 end 1321 1322 1323 %  1324 function MenuImportPoints_Callback(hObject, eventdata, handles) 1325 fileinput=browse_xml(hObject, eventdata, handles) 1326 if isempty(fileinput) 1327 return 1328 end 1329 [s,errormsg]=imadoc2struct(fileinput,'GeometryCalib'); 1330 GeometryCalib=s.GeometryCalib; 1331 GeometryCalib=load_calib(hObject, eventdata, handles) 1332 calib=reshape(GeometryCalib.PointCoord',[],1); 1333 for ilist=1:numel(calib) 1334 CoordCell{ilist}=num2str(calib(ilist)); 1335 end 1336 CoordCell=reshape(CoordCell,[],5); 1337 Tabchar=cell2tab(CoordCell,'  ');%transform cells into table ready for display 1338 Tabchar=[Tabchar;{'......'}]; 1339 set(handles.ListCoord,'Value',1) 1340 set(handles.ListCoord,'String',Tabchar) 1341 MenuPlot_Callback(handles.geometry_calib, [], handles) 1342 1343 %  1344 function MenuImportIntrinsic_Callback(hObject, eventdata, handles) 1345 % 1346 fileinput=browse_xml(hObject, eventdata, handles); 1347 if isempty(fileinput) 1348 return 1349 end 1350 [s,errormsg]=imadoc2struct(fileinput,'GeometryCalib'); 1351 GeometryCalib=s.GeometryCalib; 1352 k=GeometryCalib.kappa1 * GeometryCalib.f*GeometryCalib.f; 1353 f1=GeometryCalib.f*GeometryCalib.sx/GeometryCalib.dpx; 1354 f2=GeometryCalib.f/GeometryCalib.dpy; 1355 Cx=GeometryCalib.Cx; 1356 Cy=GeometryCalib.Cy; 1357 set(handles.fx,'String',num2str(f1,'%1.1f')) 1358 set(handles.fy,'String',num2str(f2,'%1.1f')) 1359 set(handles.k,'String',num2str(k,'%1.4f')) 1360 set(handles.Cx,'String',num2str(Cx,'%1.1f')) 1361 set(handles.Cy,'String',num2str(Cy,'%1.1f')) 1362 1363 %  1364 function MenuImportAll_Callback(hObject, eventdata, handles) 1365 % 1366 fileinput=browse_xml(hObject, eventdata, handles) 1367 if ~isempty(fileinput) 1368 loadfile(handles,fileinput) 1369 end 1370 1371 %  1372 %  Executes on menubar option Import/Grid file: introduce previous grid files 1373 function MenuGridFile_Callback(hObject, eventdata, handles) 1374 %  1375 inputfile=browse_xml(hObject, eventdata, handles) 1376 listfile=get(handles.coord_files,'string'); 1377 if isequal(listfile,{''}) 1378 listfile={inputfile}; 1379 else 1380 listfile=[listfile;{inputfile}]%update the list of coord files 1381 end 1382 set(handles.coord_files,'string',listfile); 1383 1384 % 1385 function fileinput=browse_xml(hObject, eventdata, handles) 1386 % 1387 fileinput=[];%default 1388 oldfile=''; %default 1389 UserData=get(handles.geometry_calib,'UserData'); 1390 if isfield(UserData,'XmlInput') 1391 oldfile=UserData.XmlInput; 1392 end 1393 [FileName, PathName, filterindex] = uigetfile( ... 1394 {'*.xml;*.mat', ' (*.xml,*.mat)'; 1395 '*.xml', '.xml files '; ... 1396 '*.mat', '.mat matlab files '}, ... 1397 'Pick a file',oldfile); 1398 fileinput=[PathName FileName];%complete file name 1399 testblank=findstr(fileinput,' ');%look for blanks 1400 if ~isempty(testblank) 1401 msgbox_uvmat('ERROR','forbidden input file name or path: no blank character allowed') 1402 return 1403 end 1404 sizf=size(fileinput); 1405 if (~ischar(fileinput)~isequal(sizf(1),1)),return;end 1406 UserData.XmlInput=fileinput; 1407 set(handles.geometry_calib,'UserData',UserData)%record current file foer further use of browser 1408 1409 %  1410 function loadfile(handles,fileinput) 1411 % 1412 [s,errormsg]=imadoc2struct(fileinput,'GeometryCalib'); 1413 GeometryCalib=s.GeometryCalib; 1414 if isempty(GeometryCalib) 1415 Tabchar={}; 1416 CoordCell={}; 1417 k=0;%default 1418 f1=1000; 1419 f2=1000; 1420 hhuvmat=guidata(findobj(allchild(0),'Name','uvmat')); 1421 Cx=str2num(get(hhuvmat.npx,'String'))/2; 1422 Cy=str2num(get(hhuvmat.npy,'String'))/2; 1423 else 1424 k=GeometryCalib.kappa1 * GeometryCalib.f*GeometryCalib.f; 1425 f1=GeometryCalib.f*GeometryCalib.sx/GeometryCalib.dpx; 1426 f2=GeometryCalib.f/GeometryCalib.dpy; 1427 Cx=GeometryCalib.Cx; 1428 Cy=GeometryCalib.Cy; 1429 set(handles.fx,'String',num2str(f1,'%1.1f')) 1430 set(handles.fy,'String',num2str(f2,'%1.1f')) 1431 set(handles.k,'String',num2str(k,'%1.4f')) 1432 set(handles.Cx,'String',num2str(Cx,'%1.1f')) 1433 set(handles.Cy,'String',num2str(Cy,'%1.1f')) 1434 calib=reshape(GeometryCalib.PointCoord',[],1); 1435 for ilist=1:numel(calib) 1436 CoordCell{ilist}=num2str(calib(ilist)); 1437 end 1438 CoordCell=reshape(CoordCell,[],5); 1439 Tabchar=cell2tab(CoordCell,'  ');%transform cells into table ready for display 1440 set(handles.ListCoord,'Value',1) 1441 set(handles.ListCoord,'String',Tabchar) 1442 MenuPlot_Callback(handles.geometry_calib, [], handles) 1443 end 1444 Tabchar=[Tabchar;{'......'}]; 1445 if isempty(CoordCell)% allow mouse action by default in the absence of input points 1446 set(handles.edit_append,'Value',1) 1447 set(handles.edit_append,'BackgroundColor',[1 1 0]) 1448 else % does not allow mouse action by default in the presence of input points 1449 set(handles.edit_append,'Value',0) 1450 set(handles.edit_append,'BackgroundColor',[0.7 0.7 0.7]) 1451 end 1452 1453 1454 1455 1456 1457 %  1458 %  Executes on button press in CLEAR_PTS: clear the list of calibration points 1459 function CLEAR_PTS_Callback(hObject, eventdata, handles) 1460 %  1461 set(handles.ListCoord,'Value',1)% refresh the display of coordinates 1462 set(handles.ListCoord,'String',{'......'}) 1463
Note: See TracChangeset
for help on using the changeset viewer.