Ignore:
Timestamp:
Sep 28, 2010, 2:00:22 PM (10 years ago)
Author:
gostiaux
Message:

Changes from Joel on detect_grid

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/geometry_calib.m

    r88 r108  
    506506
    507507% est_kc=[1;0;0;0;0];
    508 est_dist=[0;0;0;0;0];
     508est_fc=0;
     509est_dist=[1;0;0;0;0];
    509510run(fullfile(path_UVMAT,'toolbox_calib','go_calib_optim'));
    510511
     
    949950function MenuDetectGrid_Callback(hObject, eventdata, handles)
    950951%------------------------------------------------------------------------
    951 CalibData=get(handles.geometry_calib,'UserData');
     952CalibData=get(handles.geometry_calib,'UserData');%get information stored on the GUI geometry_calib
    952953grid_input=[];%default
    953954if isfield(CalibData,'grid')
    954955    grid_input=CalibData.grid;%retrieve the previously used grid
    955956end
    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 use
    958 
    959 %read the four last point coordiantes in pixels
    960 Coord_cell=get(handles.ListCoord,'String');%read list of coordiantes on geometry_calib
     957[T,CalibData.grid]=create_grid(grid_input,'detect grid');%display the GUI create_grid, read the set of phys coordinates T
     958set(handles.geometry_calib,'UserData',CalibData)%store the phys grid parameters for later use
     959
     960%read the four last point coordinates in pixels
     961Coord_cell=get(handles.ListCoord,'String');%read list of coordinates on geometry_calib
    961962data=read_geometry_calib(Coord_cell);
    962963nbpoints=size(data.Coord,1); %nbre of calibration points
     
    979980end
    980981
    981 %read the current image
     982%read the current image, displayed in the GUI uvmat
    982983huvmat=findobj(allchild(0),'Name','uvmat');
    983984UvData=get(huvmat,'UserData');
     
    985986npxy=size(A);
    986987%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(x1-corners_X));%error
    993 a_Y1=XY_mat\corners_Y;%transformation matrix for X
    994 y1=XY_mat*a_Y1;
    995 err_Y1=max(abs(y1-corners_Y));%error
    996 GeometryCalib.CalibrationType='linear';
     988X=[CalibData.grid.x_0 CalibData.grid.x_1 CalibData.grid.x_0 CalibData.grid.x_1]';%corner absissa in the phys coordinates
     989Y=[CalibData.grid.y_0 CalibData.grid.y_0 CalibData.grid.y_1 CalibData.grid.y_1]';%corner ordinates in the phys coordinates
     990%X=[CalibData.grid.x_0 CalibData.grid.x_0 CalibData.grid.x_1 CalibData.grid.x_1]';%corner absissa in the phys coordinates
     991%Y=[CalibData.grid.y_0 CalibData.grid.y_1 CalibData.grid.y_1 CalibData.grid.y_0]';%corner ordinates in the phys coordinate
     992%XY_mat=[ones(size(X)) X Y];
     993%a_X1=XY_mat\corners_X; %transformation matrix for X
     994%x1=XY_mat*a_X1;%reconstruction
     995%err_X1=max(abs(x1-corners_X));%error
     996%a_Y1=XY_mat\corners_Y;%transformation matrix for X
     997%y1=XY_mat*a_Y1;
     998%err_Y1=max(abs(y1-corners_Y));%error
     999% figure(1)
     1000% hold off
     1001% plot([corners_X;corners_X(1)],[corners_Y;corners_Y(1)],'r')
     1002%test points along contour
     1003% u=[0:0.1:1];
     1004% x12=u*corners_X(2)+(1-u)*corners_X(1);
     1005% x23=u*corners_X(3)+(1-u)*corners_X(2);
     1006% x34=u*corners_X(4)+(1-u)*corners_X(3);
     1007% x41=u*corners_X(1)+(1-u)*corners_X(4);
     1008% y12=u*corners_Y(2)+(1-u)*corners_Y(1);
     1009% y23=u*corners_Y(3)+(1-u)*corners_Y(2);
     1010% y34=u*corners_Y(4)+(1-u)*corners_Y(3);
     1011% y41=u*corners_Y(1)+(1-u)*corners_Y(4);
     1012% hold on
     1013% plot(x12,y12,'xr')
     1014% plot(x23,y23,'xr')
     1015% plot(x34,y34,'xr')
     1016% plot(x41,y41,'xr')
     1017%calculate transform matrices:
     1018% reference: http://alumni.media.mit.edu/~cwren/interpolator/ by Christopher R. Wren
     1019B = [ X Y ones(size(X)) zeros(4,3)        -X.*corners_X -Y.*corners_X ...
     1020      zeros(4,3)        X Y ones(size(X)) -X.*corners_Y -Y.*corners_Y ];
     1021B = reshape (B', 8 , 8 )';
     1022D = [ corners_X , corners_Y ];
     1023D = reshape (D', 8 , 1 );
     1024l = inv(B' * B) * B' * D;
     1025Amat = reshape([l(1:6)' 0 0 1 ],3,3)';
     1026C = [l(7:8)' 1];
     1027
     1028GeometryCalib.CalibrationType='tsai';%'linear';
    9971029GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function
    9981030GeometryCalib.f=1;
     
    10031035GeometryCalib.Cy=0;
    10041036GeometryCalib.kappa1=0;
    1005 GeometryCalib.Tx=a_X1(1);
    1006 GeometryCalib.Ty=a_Y1(1);
     1037GeometryCalib.Tx=Amat(1,3);
     1038GeometryCalib.Ty=Amat(2,3);
    10071039GeometryCalib.Tz=1;
    1008 GeometryCalib.R=[a_X1(2),a_X1(3),0;a_Y1(2),a_Y1(3),0;0,0,1];
     1040GeometryCalib.R=[Amat(1,1),Amat(1,2),0;Amat(2,1),Amat(2,2),0;C(1),C(2),0];
     1041%montre points du pourtour transformes
     1042% for ip=1:numel(u),
     1043%   XY12(:,ip)=Amat*[x12(ip);y12(ip);1]/(C*[x12(ip);y12(ip);1]);
     1044%   XY23(:,ip)=Amat*[x23(ip);y23(ip);1]/(C*[x23(ip);y23(ip);1]);
     1045%   XY34(:,ip)=Amat*[x34(ip);y34(ip);1]/(C*[x34(ip);y34(ip);1]);
     1046%   XY41(:,ip)=Amat*[x41(ip);y41(ip);1]/(C*[x41(ip);y41(ip);1]);
     1047% end
     1048% xphys12=u*X(2)+(1-u)*X(1);
     1049% xphys23=u*X(3)+(1-u)*X(2);
     1050% xphys34=u*X(4)+(1-u)*X(3);
     1051% xphys41=u*X(1)+(1-u)*X(4);
     1052% yphys12=u*Y(2)+(1-u)*Y(1);
     1053% yphys23=u*Y(3)+(1-u)*Y(2);
     1054% yphys34=u*Y(4)+(1-u)*Y(3);
     1055% yphys41=u*Y(1)+(1-u)*Y(4);
     1056% for ip=1:numel(u),
     1057%   XY12(:,ip)=Amat*[xphys12(ip);yphys12(ip);1]/(C*[xphys12(ip);yphys12(ip);1]);
     1058%   XY23(:,ip)=Amat*[xphys23(ip);yphys23(ip);1]/(C*[xphys23(ip);yphys23(ip);1]);
     1059%   XY34(:,ip)=Amat*[xphys34(ip);yphys34(ip);1]/(C*[xphys34(ip);yphys34(ip);1]);
     1060%   XY41(:,ip)=Amat*[xphys41(ip);yphys41(ip);1]/(C*[xphys41(ip);yphys41(ip);1]);
     1061% end
     1062% x12=XY12(1,:);
     1063% y12=XY12(2,:);
     1064% x23=XY23(1,:);
     1065% y23=XY23(2,:);
     1066% x34=XY34(1,:);
     1067% y34=XY34(2,:);
     1068% x41=XY41(1,:);
     1069% y41=XY41(2,:);
     1070% [x12,y12]=px_XYZ(GeometryCalib,xphys12,yphys12);
     1071% [x23,y23]=px_XYZ(GeometryCalib,xphys23,yphys23);
     1072% [x34,y34]=px_XYZ(GeometryCalib,xphys34,yphys34);
     1073% [x41,y41]=px_XYZ(GeometryCalib,xphys41,yphys41);
     1074% plot(x12,y12,'xb')
     1075% plot(x23,y23,'xb')
     1076% plot(x34,y34,'xb')
     1077% plot(x41,y41,'xb')
     1078% figure(2)
     1079% hold off
     1080% plot(XY12(1,:),XY12(2,:),'ob') 
     1081% hold on
     1082% plot(XY23(1,:),XY23(2,:),'ob')
     1083% plot(XY34(1,:),XY34(2,:),'ob')
     1084% plot(XY41(1,:),XY41(2,:),'ob')
     1085% plot(xphys12,yphys12,'ob') 
     1086% hold on
     1087% plot(xphys23,yphys23,'ob')
     1088% plot(xphys34,yphys34,'ob')
     1089% plot(xphys41,yphys41,'ob')
     1090
    10091091[Amod,Rangx,Rangy]=phys_Ima(A-min(min(A)),GeometryCalib,0);
     1092
    10101093Amod=double(Amod);
    1011 %figure(12)
    1012 %Amax=max(max(Amod))
    1013 %image(Rangx,Rangy,uint8(255*Amod/Amax))
     1094figure(12)
     1095Amax=max(max(Amod));
     1096image(Rangx,Rangy,uint8(255*Amod/Amax))
    10141097Dx=(Rangx(2)-Rangx(1))/(npxy(2)-1); %x mesh in real space
    10151098Dy=(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 points
    1017 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
     1099ind_range_x=ceil(GeometryCalib.R(1,1)*CalibData.grid.Dx/3);% range of search of image ma around each point obtained by linear interpolation from the marked points
     1100ind_range_y=ceil(GeometryCalib.R(2,2)*CalibData.grid.Dy/3);% range of search of image ma around each point obtained by linear interpolation from the marked points
    10181101nbpoints=size(T,1);
    10191102for ipoint=1:nbpoints
     
    11271210set(handles.ListCoord,'String',Tabchar)
    11281211
    1129 
    1130 %%%%%%%%%%%%%%%%%%%%
     1212%------------------------------------------------------------------------
     1213% image transform from px to phys
     1214%INPUT:
     1215%Zindex: index of plane
    11311216function [A_out,Rangx,Rangy]=phys_Ima(A,Calib,ZIndex)
     1217%------------------------------------------------------------------------
    11321218xcorner=[];
    11331219ycorner=[];
    11341220npx=[];
    11351221npy=[];
    1136 siz=size(A);
     1222siz=size(A)
    11371223npx=[npx siz(2)];
    1138 npy=[npy siz(1)];
    1139 xima=[0.5 siz(2)-0.5 0.5 siz(2)-0.5];%image coordiantes of corners
     1224npy=[npy siz(1)]
     1225xima=[0.5 siz(2)-0.5 0.5 siz(2)-0.5];%image coordinates of corners
    11401226yima=[0.5 0.5 siz(1)-0.5 siz(1)-0.5];
    11411227[xcorner,ycorner]=phys_XYZ(Calib,xima,yima,ZIndex);%corresponding physical coordinates
     
    11911277end
    11921278
     1279%------------------------------------------------------------------------
     1280% pointwise transform from px to phys
    11931281%INPUT:
    11941282%Z: index of plane
    11951283function [Xphys,Yphys,Zphys]=phys_XYZ(Calib,X,Y,Z)
     1284%------------------------------------------------------------------------
    11961285if exist('Z','var')& isequal(Z,round(Z))& Z>0 & isfield(Calib,'SliceCoord')&length(Calib.SliceCoord)>=Z
    11971286    Zindex=Z;
Note: See TracChangeset for help on using the changeset viewer.