Changeset 121 for trunk/src/geometry_calib.m
- Timestamp:
- Nov 12, 2010, 5:39:24 PM (14 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/src/geometry_calib.m
r116 r121 1 %% 1 2 %'geometry_calib': performs geometric calibration from a set of reference points 2 3 % … … 53 54 'gui_LayoutFcn', [] , ... 54 55 'gui_Callback', []); 55 if nargin & isstr(varargin{1})56 if nargin && ischar(varargin{1}) 56 57 gui_State.gui_Callback = str2func(varargin{1}); 57 58 end … … 90 91 %set menu of calibration options 91 92 %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'})93 set(handles.calib_type,'String',{'rescale';'linear';'3D_linear';'3D_quadr';'3D_extrinsic'}) 93 94 inputxml=''; 94 95 if exist('inputfile','var')&& ~isempty(inputfile) … … 97 98 [Pathsub,RootFile,field_count,str2,str_a,str_b,ext,nom_type,subdir]=name2display(inputfile); 98 99 if ~strcmp(ext,'.xml') 99 inputfile=[fullfile(Pathsub,RootFile) '.xml'] %xml file corresponding to the input file100 inputfile=[fullfile(Pathsub,RootFile) '.xml'];%xml file corresponding to the input file 100 101 end 101 102 end … … 141 142 function APPLY_Callback(hObject, eventdata, handles) 142 143 %------------------------------------------------------------------------ 143 144 144 %read the current calibration points 145 145 Coord_cell=get(handles.ListCoord,'String'); 146 146 Object=read_geometry_calib(Coord_cell); 147 147 Coord=Object.Coord; 148 149 148 % apply the calibration, whose type is selected in handles.calib_type 150 149 if ~isempty(Coord) … … 198 197 RootPath=''; 199 198 RootFile=''; 200 if ~isempty(hhuvmat.RootPath)& ~isempty(hhuvmat.RootFile)199 if ~isempty(hhuvmat.RootPath)&& ~isempty(hhuvmat.RootFile) 201 200 testhandle=1; 202 201 RootPath=get(hhuvmat.RootPath,'String'); … … 206 205 else 207 206 question={'save the calibration data and point coordinates in'}; 208 def={fullfile(RootPath, ['ObjectCalib.xml'])};207 def={fullfile(RootPath,'ObjectCalib.xml')}; 209 208 options.Resize='on'; 210 209 answer=inputdlg(question,'save average in a new file',1,def,options); … … 214 213 ['Error rms (along x,y)=' num2str(GeometryCalib.ErrorRms) ' pixels'];... 215 214 ['Error max (along x,y)=' num2str(GeometryCalib.ErrorMax) ' pixels']}); 216 if isequal(answer,'Yes') 217 answer_1=msgbox_uvmat('INPUT_TXT',' Z= ',num2str(Z_plane)); 218 Z_plane=str2num(answer_1); 219 GeometryCalib.NbSlice=1; 220 GeometryCalib.SliceCoord=[0 0 Z_plane]; 215 if strcmp(answer,'Yes') 216 if strcmp(calib_cell{val}(1:2),'3D')%set the plane position for 3D (projection) calibration 217 answer_1=msgbox_uvmat('INPUT_TXT',' Z= ',num2str(Z_plane)); 218 if strcmp(answer_1,'Cancel') 219 Z_plane=0; %default 220 else 221 Z_plane=str2num(answer_1); 222 end 223 GeometryCalib.NbSlice=1; 224 GeometryCalib.SliceCoord=[0 0 Z_plane]; 225 end 221 226 errormsg=update_imadoc(GeometryCalib,outputfile);% introduce the calibration data in the xml file 222 227 if ~strcmp(errormsg,'') … … 396 401 T_y=py(2); 397 402 GeometryCalib.CalibrationType='rescale'; 398 GeometryCalib.fx_fy=[px(1) py(1)]; 403 GeometryCalib.fx_fy=[px(1) py(1)];%.fx_fy corresponds to pxcm along x and y 399 404 GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function 400 405 GeometryCalib.Tx_Ty_Tz=[px(2)/px(1) py(2)/py(1) 1]; … … 417 422 % y1=XY_mat*a_Y1; 418 423 % err_Y1=max(abs(y1-y_ima));%error 419 R=[a_X1(2),a_X1(3),0;a_Y1(2),a_Y1(3),0;0,0,0]; 424 % R=[a_X1(2),a_X1(3),0;a_Y1(2),a_Y1(3),0;0,0,1]; 425 R=[a_X1(2),a_X1(3);a_Y1(2),a_Y1(3)]; 420 426 norm=abs(det(R)); 421 427 GeometryCalib.CalibrationType='linear'; 422 GeometryCalib.fx_fy=[norm norm]; 428 GeometryCalib.fx_fy(1)=sqrt((a_X1(2)/a_Y1(3))*norm); 429 GeometryCalib.fx_fy(2)=(a_Y1(3)/a_X1(2))*GeometryCalib.fx_fy(1); 423 430 GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function 424 R=[a_X1(2),a_X1(3),0;a_Y1(2),a_Y1(3),0;0,0,0];425 431 GeometryCalib.Tx_Ty_Tz=[a_X1(1) a_Y1(1) 1]; 426 GeometryCalib.R=R/norm; 427 GeometryCalib.omc=(180/pi)*[acos(R(1,1)) 0 0]; 432 R(1,:)=R(1,:)/GeometryCalib.fx_fy(1); 433 R(2,:)=R(2,:)/GeometryCalib.fx_fy(2); 434 R=[R;[0 0]]; 435 GeometryCalib.R=[R [0;0;1]]; 436 GeometryCalib.omc=(180/pi)*[acos(GeometryCalib.R(1,1)) 0 0] 428 437 %------------------------------------------------------------------------ 429 438 % determine the tsai parameters for a view normal to the grid plane … … 474 483 calib_param(3)=a_Y1(1); 475 484 calib_param(4)=Calib.f/(norm*Calib.dpx)-R(3,3)*Zmean; 476 calib_param(5)=angle(a_X1(2)+ i*a_X1(3));485 calib_param(5)=angle(a_X1(2)+1i*a_X1(3)); 477 486 display(['initial guess=' num2str(calib_param)]) 478 487 479 488 %optimise the parameters: minimisation of error 480 calib_param = fminsearch(@(calib_param) error_calib(calib_param,Calib,Coord),calib_param) 489 calib_param = fminsearch(@(calib_param) error_calib(calib_param,Calib,Coord),calib_param); 481 490 482 491 GeometryCalib.CalibrationType='tsai_normal'; … … 493 502 %------------------------------------------------------------------------ 494 503 function GeometryCalib=calib_3D_linear(Coord,handles) 495 %TO UPDATE496 504 %------------------------------------------------------------------ 497 505 path_uvmat=which('uvmat');% check the path detected for source file uvmat … … 499 507 huvmat=findobj(allchild(0),'Tag','uvmat'); 500 508 hhuvmat=guidata(huvmat); 501 x_1=Coord(:,4:5)'502 X_1=Coord(:,1:3)'503 n_ima=1;504 % check_cond=0;505 506 nx=str2num(get(hhuvmat.npx,'String'));507 ny=str2num(get(hhuvmat.npy,'String'));508 509 est_dist=[0;0;0;0;0];510 est_aspect_ratio=0;511 est_fc=[1;1];512 %fc=[25;25]/0.012;513 center_optim=0;514 run(fullfile(path_UVMAT,'toolbox_calib','go_calib_optim'));515 516 GeometryCalib.CalibrationType='3D_linear';517 GeometryCalib.focal=fc(2);518 GeometryCalib.dpx_dpy=[1 1];519 GeometryCalib.Cx_Cy=cc';520 GeometryCalib.sx=fc(1)/fc(2);521 GeometryCalib.kappa1=-kc(1)/fc(2)^2;522 GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function523 GeometryCalib.Tx_Ty_Tz=Tc_1';524 GeometryCalib.R=Rc_1;525 526 %------------------------------------------------------------------------527 function GeometryCalib=calib_3D_quadr(Coord,handles)528 %------------------------------------------------------------------529 530 path_uvmat=which('uvmat');% check the path detected for source file uvmat531 path_UVMAT=fileparts(path_uvmat); %path to UVMAT532 huvmat=findobj(allchild(0),'Tag','uvmat');533 hhuvmat=guidata(huvmat);534 % check_cond=0;535 509 coord_files=get(handles.coord_files,'String'); 536 510 if ischar(coord_files) … … 540 514 coord_files={}; 541 515 end 542 543 516 %retrieve the calibration points stored in the files listed in the popup list coord_files 544 517 x_1=Coord(:,4:5)';%px coordinates of the ref points … … 557 530 if isfield(s.GeometryCalib.SourceCalib,'PointCoord') 558 531 PointCoord=s.GeometryCalib.SourceCalib.PointCoord; 559 Coord_file= [];532 Coord_file=zeros(length(PointCoord),5);%default 560 533 for i=1:length(PointCoord) 561 534 line=str2num(PointCoord{i}); … … 571 544 end 572 545 end 573 574 546 n_ima=numel(coord_files)+1; 575 whos 576 547 est_dist=[0;0;0;0;0]; 548 est_aspect_ratio=0; 549 est_fc=[1;1]; 550 %fc=[25;25]/0.012; 551 center_optim=0; 552 run(fullfile(path_UVMAT,'toolbox_calib','go_calib_optim')); 553 GeometryCalib.CalibrationType='3D_linear'; 554 GeometryCalib.fx_fy=fc'; 555 %GeometryCalib.focal=fc(2); 556 %GeometryCalib.dpx_dpy=[1 1]; 557 GeometryCalib.Cx_Cy=cc'; 558 %GeometryCalib.sx=fc(1)/fc(2); 559 GeometryCalib.kc=kc(1); 560 %GeometryCalib.kappa1=-kc(1)/fc(2)^2; 561 GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function 562 GeometryCalib.Tx_Ty_Tz=Tc_1'; 563 GeometryCalib.R=Rc_1; 564 GeometryCalib.R(2,1:3)=-GeometryCalib.R(2,1:3);%inversion of the y image coordinate 565 GeometryCalib.Tx_Ty_Tz(2)=-GeometryCalib.Tx_Ty_Tz(2);%inversion of the y image coordinate 566 GeometryCalib.Cx_Cy(2)=ny-GeometryCalib.Cx_Cy(2);%inversion of the y image coordinate 567 GeometryCalib.omc=(180/pi)*omc_1;%angles in degrees 568 GeometryCalib.ErrorRMS=[]; 569 GeometryCalib.ErrorMax=[]; 570 571 %------------------------------------------------------------------------ 572 function GeometryCalib=calib_3D_quadr(Coord,handles) 573 %------------------------------------------------------------------ 574 575 path_uvmat=which('uvmat');% check the path detected for source file uvmat 576 path_UVMAT=fileparts(path_uvmat); %path to UVMAT 577 huvmat=findobj(allchild(0),'Tag','uvmat'); 578 hhuvmat=guidata(huvmat); 579 % check_cond=0; 580 coord_files=get(handles.coord_files,'String'); 581 if ischar(coord_files) 582 coord_files={coord_files}; 583 end 584 if isempty(coord_files{1}) || isequal(coord_files,{''}) 585 coord_files={}; 586 end 587 588 %retrieve the calibration points stored in the files listed in the popup list coord_files 589 x_1=Coord(:,4:5)';%px coordinates of the ref points 590 nx=str2num(get(hhuvmat.npx,'String')); 591 ny=str2num(get(hhuvmat.npy,'String')); 592 x_1(2,:)=ny-x_1(2,:);%reverse the y image coordinates 593 X_1=Coord(:,1:3)';%phys coordinates of the ref points 594 n_ima=numel(coord_files)+1; 595 if ~isempty(coord_files) 596 msgbox_uvmat('CONFIRMATION',['The xy coordinates of the calibration points in ' num2str(n_ima) ' planes will be used']) 597 for ifile=1:numel(coord_files) 598 t=xmltree(coord_files{ifile}); 599 s=convert(t);%convert to matlab structure 600 if isfield(s,'GeometryCalib') 601 if isfield(s.GeometryCalib,'SourceCalib') 602 if isfield(s.GeometryCalib.SourceCalib,'PointCoord') 603 PointCoord=s.GeometryCalib.SourceCalib.PointCoord; 604 Coord_file=zeros(length(PointCoord),5);%default 605 for i=1:length(PointCoord) 606 line=str2num(PointCoord{i}); 607 Coord_file(i,4:5)=line(4:5);%px x 608 Coord_file(i,1:3)=line(1:3);%phys x 609 end 610 eval(['x_' num2str(ifile+1) '=Coord_file(:,4:5)'';']); 611 eval(['x_' num2str(ifile+1) '(2,:)=ny-x_' num2str(ifile+1) '(2,:);' ]); 612 eval(['X_' num2str(ifile+1) '=Coord_file(:,1:3)'';']); 613 end 614 end 615 end 616 end 617 end 618 n_ima=numel(coord_files)+1; 577 619 est_dist=[1;0;0;0;0]; 578 620 est_aspect_ratio=1; … … 700 742 GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function 701 743 GeometryCalib.Tx_Ty_Tz=[calibdat(12) calibdat(13) calibdat(14)]; 702 Rx_Ry_Rz=calibdat( [15:17]);744 Rx_Ry_Rz=calibdat(15:17); 703 745 sa = sin(Rx_Ry_Rz(1)) ; 704 746 ca=cos(Rx_Ry_Rz(1)); … … 771 813 RootPath=''; 772 814 RootFile=''; 773 if ~isempty(hhuvmat.RootPath)& ~isempty(hhuvmat.RootFile)815 if ~isempty(hhuvmat.RootPath)&& ~isempty(hhuvmat.RootFile) 774 816 testhandle=1; 775 817 RootPath=get(hhuvmat.RootPath,'String'); … … 834 876 str2=get(handles.YObject,'String'); 835 877 str3=get(handles.ZObject,'String'); 836 if ~isempty(str1) & ~isequal(double(str1),32) & (isempty(str3)|isequal(double(str3),32))878 if ~isempty(str1) && ~isequal(double(str1),32) && (isempty(str3)||isequal(double(str3),32)) 837 879 str3='0';%put z to 0 by default 838 880 end … … 1066 1108 1067 1109 %reorder the last two points (the two first in the list) if needed 1068 angles=angle((corners_X-corners_X(1))+ i*(corners_Y-corners_Y(1)));1110 angles=angle((corners_X-corners_X(1))+1i*(corners_Y-corners_Y(1))); 1069 1111 if abs(angles(4)-angles(2))>abs(angles(3)-angles(2)) 1070 1112 X_end=corners_X(4); … … 1092 1134 D = [ corners_X , corners_Y ]; 1093 1135 D = reshape (D', 8 , 1 ); 1094 l = inv(B' * B) * B' * D; 1136 %l = inv(B' * B) * B' * D; 1137 l = (B' * B)\B' * D; 1095 1138 Amat = reshape([l(1:6)' 0 0 1 ],3,3)'; 1096 1139 C = [l(7:8)' 1]; … … 1128 1171 Data.CoordType='px'; 1129 1172 Calib.GeometryCalib=GeometryCalib; 1130 DataOut=phys(Data,Calib) 1173 DataOut=phys(Data,Calib); 1131 1174 rmpath(fullfile(path_UVMAT,'transform_field')) 1132 1175 Amod=DataOut.A; … … 1154 1197 Dx=(Rangx(2)-Rangx(1))/(npxy(2)-1); %x mesh in real space 1155 1198 Dy=(Rangy(2)-Rangy(1))/(npxy(1)-1); %y mesh in real space 1156 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 points1157 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 points1199 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 1200 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 1158 1201 nbpoints=size(T,1); 1159 1202 for ipoint=1:nbpoints … … 1405 1448 % -------------------------------------------------------------------- 1406 1449 function MenuImportPoints_Callback(hObject, eventdata, handles) 1407 fileinput=browse_xml(hObject, eventdata, handles) 1450 fileinput=browse_xml(hObject, eventdata, handles); 1408 1451 if isempty(fileinput) 1409 1452 return … … 1437 1480 function MenuImportAll_Callback(hObject, eventdata, handles) 1438 1481 %------------------------------------------------------------------------ 1439 fileinput=browse_xml(hObject, eventdata, handles) 1482 fileinput=browse_xml(hObject, eventdata, handles); 1440 1483 if ~isempty(fileinput) 1441 1484 loadfile(handles,fileinput) … … 1446 1489 function MenuGridFile_Callback(hObject, eventdata, handles) 1447 1490 % ----------------------------------------------------------------------- 1448 inputfile=browse_xml(hObject, eventdata, handles) 1491 inputfile=browse_xml(hObject, eventdata, handles); 1449 1492 listfile=get(handles.coord_files,'string'); 1450 1493 if isequal(listfile,{''}) 1451 1494 listfile={inputfile}; 1452 1495 else 1453 listfile=[listfile;{inputfile}] %update the list of coord files1496 listfile=[listfile;{inputfile}];%update the list of coord files 1454 1497 end 1455 1498 set(handles.coord_files,'string',listfile); … … 1500 1543 function loadfile(handles,fileinput) 1501 1544 %------------------------------------------------------------------------ 1502 fileinput 1503 [s,errormsg]=imadoc2struct(fileinput,'GeometryCalib') 1545 [s,errormsg]=imadoc2struct(fileinput,'GeometryCalib'); 1504 1546 GeometryCalib=s.GeometryCalib; 1505 1547 fx=1;fy=1;Cx=0;Cy=0;kc=0; %default … … 1575 1617 kc=0; 1576 1618 if isfield(GeometryCalib,'kc') 1577 kc=GeometryCalib.kc %* GeometryCalib.focal*GeometryCalib.focal;1619 kc=GeometryCalib.kc; %* GeometryCalib.focal*GeometryCalib.focal; 1578 1620 end 1579 1621 set(handles.fx,'String',num2str(fx,5))
Note: See TracChangeset
for help on using the changeset viewer.