Changeset 108
 Timestamp:
 Sep 28, 2010, 2:00:22 PM (14 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/src/geometry_calib.m
r88 r108 506 506 507 507 % est_kc=[1;0;0;0;0]; 508 est_dist=[0;0;0;0;0]; 508 est_fc=0; 509 est_dist=[1;0;0;0;0]; 509 510 run(fullfile(path_UVMAT,'toolbox_calib','go_calib_optim')); 510 511 … … 949 950 function MenuDetectGrid_Callback(hObject, eventdata, handles) 950 951 % 951 CalibData=get(handles.geometry_calib,'UserData'); 952 CalibData=get(handles.geometry_calib,'UserData');%get information stored on the GUI geometry_calib 952 953 grid_input=[];%default 953 954 if isfield(CalibData,'grid') 954 955 grid_input=CalibData.grid;%retrieve the previously used grid 955 956 end 956 [T,CalibData.grid]=create_grid(grid_input,'detect grid');%display the GUI create_grid 957 set(handles.geometry_calib,'UserData',CalibData)%store the phys grid for later use958 959 %read the four last point coordi antes in pixels960 Coord_cell=get(handles.ListCoord,'String');%read list of coordi antes on geometry_calib957 [T,CalibData.grid]=create_grid(grid_input,'detect grid');%display the GUI create_grid, read the set of phys coordinates T 958 set(handles.geometry_calib,'UserData',CalibData)%store the phys grid parameters for later use 959 960 %read the four last point coordinates in pixels 961 Coord_cell=get(handles.ListCoord,'String');%read list of coordinates on geometry_calib 961 962 data=read_geometry_calib(Coord_cell); 962 963 nbpoints=size(data.Coord,1); %nbre of calibration points … … 979 980 end 980 981 981 %read the current image 982 %read the current image, displayed in the GUI uvmat 982 983 huvmat=findobj(allchild(0),'Name','uvmat'); 983 984 UvData=get(huvmat,'UserData'); … … 985 986 npxy=size(A); 986 987 %linear transform on the current image 987 X=[CalibData.grid.x_0 CalibData.grid.x_1 CalibData.grid.x_0 CalibData.grid.x_1]';%corner absissa in the rectified image 988 Y=[CalibData.grid.y_0 CalibData.grid.y_0 CalibData.grid.y_1 CalibData.grid.y_1]';%corner absissa in the rectified image 989 XY_mat=[ones(size(X)) X Y]; 990 a_X1=XY_mat\corners_X; %transformation matrix for X 991 x1=XY_mat*a_X1;%reconstruction 992 err_X1=max(abs(x1corners_X));%error 993 a_Y1=XY_mat\corners_Y;%transformation matrix for X 994 y1=XY_mat*a_Y1; 995 err_Y1=max(abs(y1corners_Y));%error 996 GeometryCalib.CalibrationType='linear'; 988 X=[CalibData.grid.x_0 CalibData.grid.x_1 CalibData.grid.x_0 CalibData.grid.x_1]';%corner absissa in the phys coordinates 989 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') 1017 %calculate transform matrices: 1018 % reference: http://alumni.media.mit.edu/~cwren/interpolator/ by Christopher R. Wren 1019 B = [ X Y ones(size(X)) zeros(4,3) X.*corners_X Y.*corners_X ... 1020 zeros(4,3) X Y ones(size(X)) X.*corners_Y Y.*corners_Y ]; 1021 B = reshape (B', 8 , 8 )'; 1022 D = [ corners_X , corners_Y ]; 1023 D = reshape (D', 8 , 1 ); 1024 l = inv(B' * B) * B' * D; 1025 Amat = reshape([l(1:6)' 0 0 1 ],3,3)'; 1026 C = [l(7:8)' 1]; 1027 1028 GeometryCalib.CalibrationType='tsai';%'linear'; 997 1029 GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function 998 1030 GeometryCalib.f=1; … … 1003 1035 GeometryCalib.Cy=0; 1004 1036 GeometryCalib.kappa1=0; 1005 GeometryCalib.Tx= a_X1(1);1006 GeometryCalib.Ty= a_Y1(1);1037 GeometryCalib.Tx=Amat(1,3); 1038 GeometryCalib.Ty=Amat(2,3); 1007 1039 GeometryCalib.Tz=1; 1008 GeometryCalib.R=[a_X1(2),a_X1(3),0;a_Y1(2),a_Y1(3),0;0,0,1]; 1040 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 transformes 1042 % for ip=1:numel(u), 1043 % XY12(:,ip)=Amat*[x12(ip);y12(ip);1]/(C*[x12(ip);y12(ip);1]); 1044 % XY23(:,ip)=Amat*[x23(ip);y23(ip);1]/(C*[x23(ip);y23(ip);1]); 1045 % XY34(:,ip)=Amat*[x34(ip);y34(ip);1]/(C*[x34(ip);y34(ip);1]); 1046 % XY41(:,ip)=Amat*[x41(ip);y41(ip);1]/(C*[x41(ip);y41(ip);1]); 1047 % end 1048 % xphys12=u*X(2)+(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 % end 1062 % x12=XY12(1,:); 1063 % y12=XY12(2,:); 1064 % x23=XY23(1,:); 1065 % y23=XY23(2,:); 1066 % x34=XY34(1,:); 1067 % y34=XY34(2,:); 1068 % x41=XY41(1,:); 1069 % y41=XY41(2,:); 1070 % [x12,y12]=px_XYZ(GeometryCalib,xphys12,yphys12); 1071 % [x23,y23]=px_XYZ(GeometryCalib,xphys23,yphys23); 1072 % [x34,y34]=px_XYZ(GeometryCalib,xphys34,yphys34); 1073 % [x41,y41]=px_XYZ(GeometryCalib,xphys41,yphys41); 1074 % plot(x12,y12,'xb') 1075 % plot(x23,y23,'xb') 1076 % plot(x34,y34,'xb') 1077 % plot(x41,y41,'xb') 1078 % figure(2) 1079 % hold off 1080 % plot(XY12(1,:),XY12(2,:),'ob') 1081 % hold on 1082 % plot(XY23(1,:),XY23(2,:),'ob') 1083 % plot(XY34(1,:),XY34(2,:),'ob') 1084 % plot(XY41(1,:),XY41(2,:),'ob') 1085 % plot(xphys12,yphys12,'ob') 1086 % hold on 1087 % plot(xphys23,yphys23,'ob') 1088 % plot(xphys34,yphys34,'ob') 1089 % plot(xphys41,yphys41,'ob') 1090 1009 1091 [Amod,Rangx,Rangy]=phys_Ima(Amin(min(A)),GeometryCalib,0); 1092 1010 1093 Amod=double(Amod); 1011 %figure(12)1012 %Amax=max(max(Amod)) 1013 %image(Rangx,Rangy,uint8(255*Amod/Amax))1094 figure(12) 1095 Amax=max(max(Amod)); 1096 image(Rangx,Rangy,uint8(255*Amod/Amax)) 1014 1097 Dx=(Rangx(2)Rangx(1))/(npxy(2)1); %x mesh in real space 1015 1098 Dy=(Rangy(2)Rangy(1))/(npxy(1)1); %y mesh in real space 1016 ind_range_x=ceil(GeometryCalib.R(1,1)*CalibData.grid.Dx/3) % range of search of image ma around each point obtained by linear interpolation from the marked points1017 ind_range_y=ceil(GeometryCalib.R(2,2)*CalibData.grid.Dy/3) % range of search of image ma around each point obtained by linear interpolation from the marked points1099 ind_range_x=ceil(GeometryCalib.R(1,1)*CalibData.grid.Dx/3);% range of search of image ma around each point obtained by linear interpolation from the marked points 1100 ind_range_y=ceil(GeometryCalib.R(2,2)*CalibData.grid.Dy/3);% range of search of image ma around each point obtained by linear interpolation from the marked points 1018 1101 nbpoints=size(T,1); 1019 1102 for ipoint=1:nbpoints … … 1127 1210 set(handles.ListCoord,'String',Tabchar) 1128 1211 1129 1130 %%%%%%%%%%%%%%%%%%%% 1212 % 1213 % image transform from px to phys 1214 %INPUT: 1215 %Zindex: index of plane 1131 1216 function [A_out,Rangx,Rangy]=phys_Ima(A,Calib,ZIndex) 1217 % 1132 1218 xcorner=[]; 1133 1219 ycorner=[]; 1134 1220 npx=[]; 1135 1221 npy=[]; 1136 siz=size(A) ;1222 siz=size(A) 1137 1223 npx=[npx siz(2)]; 1138 npy=[npy siz(1)] ;1139 xima=[0.5 siz(2)0.5 0.5 siz(2)0.5];%image coordi antes of corners1224 npy=[npy siz(1)] 1225 xima=[0.5 siz(2)0.5 0.5 siz(2)0.5];%image coordinates of corners 1140 1226 yima=[0.5 0.5 siz(1)0.5 siz(1)0.5]; 1141 1227 [xcorner,ycorner]=phys_XYZ(Calib,xima,yima,ZIndex);%corresponding physical coordinates … … 1191 1277 end 1192 1278 1279 % 1280 % pointwise transform from px to phys 1193 1281 %INPUT: 1194 1282 %Z: index of plane 1195 1283 function [Xphys,Yphys,Zphys]=phys_XYZ(Calib,X,Y,Z) 1284 % 1196 1285 if exist('Z','var')& isequal(Z,round(Z))& Z>0 & isfield(Calib,'SliceCoord')&length(Calib.SliceCoord)>=Z 1197 1286 Zindex=Z;
Note: See TracChangeset
for help on using the changeset viewer.