Index: /trunk/src/transform_field/phys_polar.m =================================================================== --- /trunk/src/transform_field/phys_polar.m (revision 566) +++ /trunk/src/transform_field/phys_polar.m (revision 567) @@ -40,13 +40,13 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% origin_xy=[0 0];%center for the polar coordinates in the original x,y coordinates -if isfield(Calib{1},'PolarCentre') && isnumeric(Calib{1}.PolarCentre) - if isequal(length(Calib{1}.PolarCentre),2); - origin_xy= Calib{1}.PolarCentre; +if isfield(CalibData,'PolarCentre') && isnumeric(CalibData.PolarCentre) + if isequal(length(CalibData.PolarCentre),2); + origin_xy= CalibData.PolarCentre; end end radius_offset=0;%reference radius used to offset the radial coordinate r angle_offset=0; %reference angle used as new origin of the polar angle (= axis Ox by default) -if isfield(Calib{1},'PolarReferenceRadius') && isnumeric(Calib{1}.PolarReferenceRadius) - radius_offset=Calib{1}.PolarReferenceRadius; +if isfield(CalibData,'PolarReferenceRadius') && isnumeric(CalibData.PolarReferenceRadius) + radius_offset=CalibData.PolarReferenceRadius; end if radius_offset > 0 @@ -55,6 +55,6 @@ angle_scale=180/pi; %polar angle in degrees end -if isfield(Calib{1},'PolarReferenceAngle') && isnumeric(Calib{1}.PolarReferenceAngle) - angle_offset=Calib{1}.PolarReferenceAngle; %offset angle (in unit of the final angle, degrees or arc length along the reference radius)) +if isfield(CalibData,'PolarReferenceAngle') && isnumeric(CalibData.PolarReferenceAngle) + angle_offset=CalibData.PolarReferenceAngle; %offset angle (in unit of the final angle, degrees or arc length along the reference radius)) end % new x coordinate = radius-radius_offset; @@ -80,4 +80,5 @@ end end + %transform second field (if exists) to cartesian phys coordiantes if test_1 @@ -184,74 +185,4 @@ -%------------------------------------------------------------------------ -%'phys_XYZ':transforms image (px) to real world (phys) coordinates using geometric calibration parameters -% function [Xphys,Yphys]=phys_XYZ(Calib,X,Y,Z) -% -%OUTPUT: -% -%INPUT: -%Z: index of plane -function [Xphys,Yphys,Zphys]=phys_XYZ(Calib,X,Y,Z) -%------------------------------------------------------------------------ -if exist('Z','var')&& isequal(Z,round(Z))&& Z>0 && isfield(Calib,'SliceCoord')&&length(Calib.SliceCoord)>=Z - Zindex=Z; - Zphys=Calib.SliceCoord(Zindex,3);%GENERALISER AUX CAS AVEC ANGLE -else - Zphys=0; -end -if ~exist('X','var')||~exist('Y','var') - Xphys=[]; - Yphys=[];%default - return -end -%coordinate transform -if ~isfield(Calib,'fx_fy') - Calib.fx_fy=[1 1]; -end -if ~isfield(Calib,'Tx_Ty_Tz') - Calib.Tx_Ty_Tz=[0 0 1]; -end -if ~isfield(Calib,'Cx_Cy') - Calib.Cx_Cy=[0 0]; -end -if ~isfield(Calib,'kc') - Calib.kc=0; -end -if isfield(Calib,'R') - R=(Calib.R)'; - Tx=Calib.Tx_Ty_Tz(1); - Ty=Calib.Tx_Ty_Tz(2); - Tz=Calib.Tx_Ty_Tz(3); - f=Calib.fx_fy(1);%dpy=1; sx=1 - dpx=Calib.fx_fy(2)/Calib.fx_fy(1); - Dx=R(5)*R(7)-R(4)*R(8); - Dy=R(1)*R(8)-R(2)*R(7); - D0=f*(R(2)*R(4)-R(1)*R(5)); - Z11=R(6)*R(8)-R(5)*R(9); - Z12=R(2)*R(9)-R(3)*R(8); - Z21=R(4)*R(9)-R(6)*R(7); - Z22=R(3)*R(7)-R(1)*R(9); - Zx0=R(3)*R(5)-R(2)*R(6); - Zy0=R(1)*R(6)-R(3)*R(4); - A11=R(8)*Ty-R(5)*Tz+Z11*Zphys; - A12=R(2)*Tz-R(8)*Tx+Z12*Zphys; - A21=-R(7)*Ty+R(4)*Tz+Z21*Zphys; - A22=-R(1)*Tz+R(7)*Tx+Z22*Zphys; - X0=f*(R(5)*Tx-R(2)*Ty+Zx0*Zphys); - Y0=f*(-R(4)*Tx+R(1)*Ty+Zy0*Zphys); - %px to camera: - Xd=dpx*(X-Calib.Cx_Cy(1)); % sensor coordinates - Yd=(Y-Calib.Cx_Cy(2)); - dist_fact=1+Calib.kc*(Xd.*Xd+Yd.*Yd)/(f*f); %distortion factor - Xu=Xd./dist_fact;%undistorted sensor coordinates - Yu=Yd./dist_fact; - denom=Dx*Xu+Dy*Yu+D0; - Xphys=(A11.*Xu+A12.*Yu+X0)./denom;%world coordinates - Yphys=(A21.*Xu+A22.*Yu+Y0)./denom; -else - Xphys=-Calib.Tx_Ty_Tz(1)+X/Calib.fx_fy(1); - Yphys=-Calib.Tx_Ty_Tz(2)+Y/Calib.fx_fy(2); -end - %%%%%%%%%%%%%%%%%%%% function [A_out,Rangx,Rangy]=phys_Ima_polar(A,CalibIn,ZIndex,origin_xy,radius_offset,angle_offset,angle_scale) Index: /trunk/src/uvmat.m =================================================================== --- /trunk/src/uvmat.m (revision 566) +++ /trunk/src/uvmat.m (revision 567) @@ -4598,9 +4598,18 @@ function MenuLIFCalib_Callback(hObject, eventdata, handles) %------------------------------------------------------------------------ -UvData=get(handles.uvmat,'UserData');%read UvData properties stored on the uvmat interface +%% read UvData properties stored on the uvmat interface +UvData=get(handles.uvmat,'UserData'); +if isfield(UvData,'XmlData')&& isfield(UvData.XmlData{1},'GeometryCalib') + XmlData=UvData.XmlData{1}; +else + msgbox_uvmat('ERROR','geometric calibration needed: use Tools/geometric calibration in the menu bar'); + return +end + +%% read lines currently drawn ListObj=UvData.Object; select=zeros(1,numel(ListObj)); for iobj=1:numel(ListObj); - if strcmp(ListObj{iobj}.Type,'line') + if isfield(ListObj{iobj},'Type') && strcmp(ListObj{iobj}.Type,'line') select(iobj)=1; end @@ -4611,36 +4620,106 @@ return else - set(handles.ListObject,'Value',val); + set(handles.ListObject,'Value',val);% show the selected lines on the list ObjectData=UvData.Object(val); - flag=1; - npx=size(UvData.Field.A,2); - npy=size(UvData.Field.A,1); - xi=0.5:npx-0.5; - yi=0.5:npy-0.5; - [Xi,Yi]=meshgrid(xi,yi); for iobj=1:length(ObjectData) - flagobj=1; - testphys=0; %coordinates in pixels by default - if isfield(ObjectData,'CoordUnit') && ~isequal(ObjectData.CoordUnit,'pixel') - if isfield(UvData,'XmlData')&& isfield(UvData.XmlData{1},'GeometryCalib') - Calib=UvData.XmlData{1}.GeometryCalib; - testphys=1; - end - end - if isfield(ObjectData{iobj},'Coord') - x1(iobj)=ObjectData{iobj}.Coord(1,1); - y1(iobj)=ObjectData{iobj}.Coord(1,2); - x2(iobj)=ObjectData{iobj}.Coord(2,1); - y2(iobj)=ObjectData{iobj}.Coord(2,2); - end - end -end - %determine the ray origin - x1 - y1 - x2 - y2 - % update the xml file - +% if isfield(ObjectData{iobj},'Coord') + xA(iobj)=ObjectData{iobj}.Coord(1,1); + yA(iobj)=ObjectData{iobj}.Coord(1,2); + xB(iobj)=ObjectData{iobj}.Coord(2,1); + yB(iobj)=ObjectData{iobj}.Coord(2,2); +% end + end +end + +%% find the origin as intersection of the two first lines (see http://www.ahristov.com/tutorial/geometry-games/intersection-lines.html ) +x1=xA(1);x2=xB(1); +x3=xA(2);x4=xB(2); +y1=yA(1);y2=yB(1); +y3=yA(2);y4=yB(2); +D = (x1-x2)*(y3-y4) -(y1-y2)*(x3-x4); +if D==0 + msgbox_uvmat('ERROR','the two lines are parallel'); + return +end +x0=((x3-x4)*(x1*y2-y1*x2)-(x1-x2)*(x3*y4-y3*x4))/D; +y0=((y3-y4)*(x1*y2-y1*x2)-(y1-y2)*(x3*y4-y3*x4))/D; +XmlData.Illumination.Origin=[x0 y0]; +XmlData.PolarCentre=[x0 y0]; + +%% display the current image in polar coordinates with origin at the illumination source +currentdir=pwd; +uvmatpath=fileparts(which('uvmat')); +cd(fullfile(uvmatpath,'transform_field')); +phys_polar=str2func('phys_polar'); +cd(currentdir) +DataOut=phys_polar(UvData.Field,XmlData); +view_field(DataOut); + +%% use the third line for reference luminosity +if numel(val)==3 + x_ref=linspace(ObjectData{3}.Coord(1,1),ObjectData{3}.Coord(2,1),10); + y_ref=linspace(ObjectData{3}.Coord(1,2),ObjectData{3}.Coord(2,2),10); + x_ref=x_ref-x0; + y_ref=y_ref-y0; + [theta_ref,r_ref] = cart2pol(x_ref,y_ref);%theta_ref and r_ref are the polar coordinates of the points on the line + theta_ref=theta_ref*180/pi; + figure + plot(theta_ref,r_ref) + azimuth_ima=linspace(DataOut.AY(1),DataOut.AY(2),size(DataOut.A,1));%profile of x index on the transformed image + dist_source = interp1(theta_ref,r_ref,azimuth_ima); + dist_source_pixel=round(size(DataOut.A,2)*(dist_source-DataOut.AX(1))/(DataOut.AX(2)-DataOut.AX(1))); + line_nan= isnan(dist_source_pixel); + dist_source_pixel(line_nan)=1; + width=20; %number of pixels used for reference + DataOut.A=double(DataOut.A); + Anorm=zeros(size(DataOut.A)); + Aval=mean(mean(DataOut.A)); + for iline=1:size(DataOut.A,1) + lum(iline)=mean(DataOut.A(iline,dist_source_pixel(iline):dist_source_pixel(iline)+width)); + Anorm(iline,:)=uint16(Aval*DataOut.A(iline,:)/lum(iline)); + end + lum(line_nan)=NaN; + figure + plot(1:size(DataOut.A,1),lum) +end +ImaName=regexprep([get(handles.RootFile,'String') get(handles.FileIndex,'String')],'//',''); +NewImageName=fullfile(get(handles.RootPath,'String'),'polar',[ImaName get(handles.FileExt,'String')]); +imwrite(Anorm,NewImageName,'BitDepth',16) + +%% record the origin in the xml file +XmlFileName=find_imadoc(get(handles.RootPath,'String'),get(handles.SubDir,'String'),get(handles.RootFile,'String'),get(handles.FileExt,'String')); +answer=msgbox_uvmat('INPUT_Y-N','save the illumination origin in the current xml file?'); +if strcmp(answer,'Yes') + t=xmltree(XmlFileName); %read the file + title=get(t,1,'name'); + if ~strcmp(title,'ImaDoc') + msgbox_uvmat('ERROR','wrong xml file'); + return + end + % backup the output file if it already exist, and read it + backupfile=XmlFileName; + testexist=2; + while testexist==2 + backupfile=[backupfile '~']; + testexist=exist(backupfile,'file'); + end + [success,message]=copyfile(XmlFileName,backupfile);%make backup + if success~=1 + errormsg=['errror in xml file backup: ' message]; + return + end + uid_illumination=find(t,'ImaDoc/Illumination'); + if isempty(uid_illumination) %if GeometryCalib does not already exists, create it + [t,uid_illumination]=add(t,1,'element','Illumination'); + end + uid_origin=find(t,'ImaDoc/Illumination/Origin'); + if ~isempty(uid_origin) %if GeometryCalib does not already exists, create it + t=delete(t,uid_origin); + end + % save the illumination origin + t=struct2xml(XmlData.Illumination,t,uid_illumination); + save(t,XmlFileName); +end +