 Timestamp:
 Jan 31, 2013, 2:24:00 PM (11 years ago)
 Location:
 trunk/src
 Files:

 2 edited
Legend:
 Unmodified
 Added
 Removed

trunk/src/transform_field/phys_polar.m
r172 r567 40 40 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 41 41 origin_xy=[0 0];%center for the polar coordinates in the original x,y coordinates 42 if isfield(Calib {1},'PolarCentre') && isnumeric(Calib{1}.PolarCentre)43 if isequal(length(Calib {1}.PolarCentre),2);44 origin_xy= Calib {1}.PolarCentre;42 if isfield(CalibData,'PolarCentre') && isnumeric(CalibData.PolarCentre) 43 if isequal(length(CalibData.PolarCentre),2); 44 origin_xy= CalibData.PolarCentre; 45 45 end 46 46 end 47 47 radius_offset=0;%reference radius used to offset the radial coordinate r 48 48 angle_offset=0; %reference angle used as new origin of the polar angle (= axis Ox by default) 49 if isfield(Calib {1},'PolarReferenceRadius') && isnumeric(Calib{1}.PolarReferenceRadius)50 radius_offset=Calib {1}.PolarReferenceRadius;49 if isfield(CalibData,'PolarReferenceRadius') && isnumeric(CalibData.PolarReferenceRadius) 50 radius_offset=CalibData.PolarReferenceRadius; 51 51 end 52 52 if radius_offset > 0 … … 55 55 angle_scale=180/pi; %polar angle in degrees 56 56 end 57 if isfield(Calib {1},'PolarReferenceAngle') && isnumeric(Calib{1}.PolarReferenceAngle)58 angle_offset=Calib {1}.PolarReferenceAngle; %offset angle (in unit of the final angle, degrees or arc length along the reference radius))57 if isfield(CalibData,'PolarReferenceAngle') && isnumeric(CalibData.PolarReferenceAngle) 58 angle_offset=CalibData.PolarReferenceAngle; %offset angle (in unit of the final angle, degrees or arc length along the reference radius)) 59 59 end 60 60 % new x coordinate = radiusradius_offset; … … 80 80 end 81 81 end 82 82 83 %transform second field (if exists) to cartesian phys coordiantes 83 84 if test_1 … … 184 185 185 186 186 %187 %'phys_XYZ':transforms image (px) to real world (phys) coordinates using geometric calibration parameters188 % function [Xphys,Yphys]=phys_XYZ(Calib,X,Y,Z)189 %190 %OUTPUT:191 %192 %INPUT:193 %Z: index of plane194 function [Xphys,Yphys,Zphys]=phys_XYZ(Calib,X,Y,Z)195 %196 if exist('Z','var')&& isequal(Z,round(Z))&& Z>0 && isfield(Calib,'SliceCoord')&&length(Calib.SliceCoord)>=Z197 Zindex=Z;198 Zphys=Calib.SliceCoord(Zindex,3);%GENERALISER AUX CAS AVEC ANGLE199 else200 Zphys=0;201 end202 if ~exist('X','var')~exist('Y','var')203 Xphys=[];204 Yphys=[];%default205 return206 end207 %coordinate transform208 if ~isfield(Calib,'fx_fy')209 Calib.fx_fy=[1 1];210 end211 if ~isfield(Calib,'Tx_Ty_Tz')212 Calib.Tx_Ty_Tz=[0 0 1];213 end214 if ~isfield(Calib,'Cx_Cy')215 Calib.Cx_Cy=[0 0];216 end217 if ~isfield(Calib,'kc')218 Calib.kc=0;219 end220 if isfield(Calib,'R')221 R=(Calib.R)';222 Tx=Calib.Tx_Ty_Tz(1);223 Ty=Calib.Tx_Ty_Tz(2);224 Tz=Calib.Tx_Ty_Tz(3);225 f=Calib.fx_fy(1);%dpy=1; sx=1226 dpx=Calib.fx_fy(2)/Calib.fx_fy(1);227 Dx=R(5)*R(7)R(4)*R(8);228 Dy=R(1)*R(8)R(2)*R(7);229 D0=f*(R(2)*R(4)R(1)*R(5));230 Z11=R(6)*R(8)R(5)*R(9);231 Z12=R(2)*R(9)R(3)*R(8);232 Z21=R(4)*R(9)R(6)*R(7);233 Z22=R(3)*R(7)R(1)*R(9);234 Zx0=R(3)*R(5)R(2)*R(6);235 Zy0=R(1)*R(6)R(3)*R(4);236 A11=R(8)*TyR(5)*Tz+Z11*Zphys;237 A12=R(2)*TzR(8)*Tx+Z12*Zphys;238 A21=R(7)*Ty+R(4)*Tz+Z21*Zphys;239 A22=R(1)*Tz+R(7)*Tx+Z22*Zphys;240 X0=f*(R(5)*TxR(2)*Ty+Zx0*Zphys);241 Y0=f*(R(4)*Tx+R(1)*Ty+Zy0*Zphys);242 %px to camera:243 Xd=dpx*(XCalib.Cx_Cy(1)); % sensor coordinates244 Yd=(YCalib.Cx_Cy(2));245 dist_fact=1+Calib.kc*(Xd.*Xd+Yd.*Yd)/(f*f); %distortion factor246 Xu=Xd./dist_fact;%undistorted sensor coordinates247 Yu=Yd./dist_fact;248 denom=Dx*Xu+Dy*Yu+D0;249 Xphys=(A11.*Xu+A12.*Yu+X0)./denom;%world coordinates250 Yphys=(A21.*Xu+A22.*Yu+Y0)./denom;251 else252 Xphys=Calib.Tx_Ty_Tz(1)+X/Calib.fx_fy(1);253 Yphys=Calib.Tx_Ty_Tz(2)+Y/Calib.fx_fy(2);254 end255 256 187 %%%%%%%%%%%%%%%%%%%% 257 188 function [A_out,Rangx,Rangy]=phys_Ima_polar(A,CalibIn,ZIndex,origin_xy,radius_offset,angle_offset,angle_scale) 
trunk/src/uvmat.m
r566 r567 4598 4598 function MenuLIFCalib_Callback(hObject, eventdata, handles) 4599 4599 % 4600 UvData=get(handles.uvmat,'UserData');%read UvData properties stored on the uvmat interface 4600 %% read UvData properties stored on the uvmat interface 4601 UvData=get(handles.uvmat,'UserData'); 4602 if isfield(UvData,'XmlData')&& isfield(UvData.XmlData{1},'GeometryCalib') 4603 XmlData=UvData.XmlData{1}; 4604 else 4605 msgbox_uvmat('ERROR','geometric calibration needed: use Tools/geometric calibration in the menu bar'); 4606 return 4607 end 4608 4609 %% read lines currently drawn 4601 4610 ListObj=UvData.Object; 4602 4611 select=zeros(1,numel(ListObj)); 4603 4612 for iobj=1:numel(ListObj); 4604 if strcmp(ListObj{iobj}.Type,'line')4613 if isfield(ListObj{iobj},'Type') && strcmp(ListObj{iobj}.Type,'line') 4605 4614 select(iobj)=1; 4606 4615 end … … 4611 4620 return 4612 4621 else 4613 set(handles.ListObject,'Value',val); 4622 set(handles.ListObject,'Value',val);% show the selected lines on the list 4614 4623 ObjectData=UvData.Object(val); 4615 flag=1;4616 npx=size(UvData.Field.A,2);4617 npy=size(UvData.Field.A,1);4618 xi=0.5:npx0.5;4619 yi=0.5:npy0.5;4620 [Xi,Yi]=meshgrid(xi,yi);4621 4624 for iobj=1:length(ObjectData) 4622 flagobj=1; 4623 testphys=0; %coordinates in pixels by default 4624 if isfield(ObjectData,'CoordUnit') && ~isequal(ObjectData.CoordUnit,'pixel') 4625 if isfield(UvData,'XmlData')&& isfield(UvData.XmlData{1},'GeometryCalib') 4626 Calib=UvData.XmlData{1}.GeometryCalib; 4627 testphys=1; 4628 end 4629 end 4630 if isfield(ObjectData{iobj},'Coord') 4631 x1(iobj)=ObjectData{iobj}.Coord(1,1); 4632 y1(iobj)=ObjectData{iobj}.Coord(1,2); 4633 x2(iobj)=ObjectData{iobj}.Coord(2,1); 4634 y2(iobj)=ObjectData{iobj}.Coord(2,2); 4635 end 4636 end 4637 end 4638 %determine the ray origin 4639 x1 4640 y1 4641 x2 4642 y2 4643 % update the xml file 4644 4625 % if isfield(ObjectData{iobj},'Coord') 4626 xA(iobj)=ObjectData{iobj}.Coord(1,1); 4627 yA(iobj)=ObjectData{iobj}.Coord(1,2); 4628 xB(iobj)=ObjectData{iobj}.Coord(2,1); 4629 yB(iobj)=ObjectData{iobj}.Coord(2,2); 4630 % end 4631 end 4632 end 4633 4634 %% find the origin as intersection of the two first lines (see http://www.ahristov.com/tutorial/geometrygames/intersectionlines.html ) 4635 x1=xA(1);x2=xB(1); 4636 x3=xA(2);x4=xB(2); 4637 y1=yA(1);y2=yB(1); 4638 y3=yA(2);y4=yB(2); 4639 D = (x1x2)*(y3y4) (y1y2)*(x3x4); 4640 if D==0 4641 msgbox_uvmat('ERROR','the two lines are parallel'); 4642 return 4643 end 4644 x0=((x3x4)*(x1*y2y1*x2)(x1x2)*(x3*y4y3*x4))/D; 4645 y0=((y3y4)*(x1*y2y1*x2)(y1y2)*(x3*y4y3*x4))/D; 4646 XmlData.Illumination.Origin=[x0 y0]; 4647 XmlData.PolarCentre=[x0 y0]; 4648 4649 %% display the current image in polar coordinates with origin at the illumination source 4650 currentdir=pwd; 4651 uvmatpath=fileparts(which('uvmat')); 4652 cd(fullfile(uvmatpath,'transform_field')); 4653 phys_polar=str2func('phys_polar'); 4654 cd(currentdir) 4655 DataOut=phys_polar(UvData.Field,XmlData); 4656 view_field(DataOut); 4657 4658 %% use the third line for reference luminosity 4659 if numel(val)==3 4660 x_ref=linspace(ObjectData{3}.Coord(1,1),ObjectData{3}.Coord(2,1),10); 4661 y_ref=linspace(ObjectData{3}.Coord(1,2),ObjectData{3}.Coord(2,2),10); 4662 x_ref=x_refx0; 4663 y_ref=y_refy0; 4664 [theta_ref,r_ref] = cart2pol(x_ref,y_ref);%theta_ref and r_ref are the polar coordinates of the points on the line 4665 theta_ref=theta_ref*180/pi; 4666 figure 4667 plot(theta_ref,r_ref) 4668 azimuth_ima=linspace(DataOut.AY(1),DataOut.AY(2),size(DataOut.A,1));%profile of x index on the transformed image 4669 dist_source = interp1(theta_ref,r_ref,azimuth_ima); 4670 dist_source_pixel=round(size(DataOut.A,2)*(dist_sourceDataOut.AX(1))/(DataOut.AX(2)DataOut.AX(1))); 4671 line_nan= isnan(dist_source_pixel); 4672 dist_source_pixel(line_nan)=1; 4673 width=20; %number of pixels used for reference 4674 DataOut.A=double(DataOut.A); 4675 Anorm=zeros(size(DataOut.A)); 4676 Aval=mean(mean(DataOut.A)); 4677 for iline=1:size(DataOut.A,1) 4678 lum(iline)=mean(DataOut.A(iline,dist_source_pixel(iline):dist_source_pixel(iline)+width)); 4679 Anorm(iline,:)=uint16(Aval*DataOut.A(iline,:)/lum(iline)); 4680 end 4681 lum(line_nan)=NaN; 4682 figure 4683 plot(1:size(DataOut.A,1),lum) 4684 end 4685 ImaName=regexprep([get(handles.RootFile,'String') get(handles.FileIndex,'String')],'//',''); 4686 NewImageName=fullfile(get(handles.RootPath,'String'),'polar',[ImaName get(handles.FileExt,'String')]); 4687 imwrite(Anorm,NewImageName,'BitDepth',16) 4688 4689 %% record the origin in the xml file 4690 XmlFileName=find_imadoc(get(handles.RootPath,'String'),get(handles.SubDir,'String'),get(handles.RootFile,'String'),get(handles.FileExt,'String')); 4691 answer=msgbox_uvmat('INPUT_YN','save the illumination origin in the current xml file?'); 4692 if strcmp(answer,'Yes') 4693 t=xmltree(XmlFileName); %read the file 4694 title=get(t,1,'name'); 4695 if ~strcmp(title,'ImaDoc') 4696 msgbox_uvmat('ERROR','wrong xml file'); 4697 return 4698 end 4699 % backup the output file if it already exist, and read it 4700 backupfile=XmlFileName; 4701 testexist=2; 4702 while testexist==2 4703 backupfile=[backupfile '~']; 4704 testexist=exist(backupfile,'file'); 4705 end 4706 [success,message]=copyfile(XmlFileName,backupfile);%make backup 4707 if success~=1 4708 errormsg=['errror in xml file backup: ' message]; 4709 return 4710 end 4711 uid_illumination=find(t,'ImaDoc/Illumination'); 4712 if isempty(uid_illumination) %if GeometryCalib does not already exists, create it 4713 [t,uid_illumination]=add(t,1,'element','Illumination'); 4714 end 4715 uid_origin=find(t,'ImaDoc/Illumination/Origin'); 4716 if ~isempty(uid_origin) %if GeometryCalib does not already exists, create it 4717 t=delete(t,uid_origin); 4718 end 4719 % save the illumination origin 4720 t=struct2xml(XmlData.Illumination,t,uid_illumination); 4721 save(t,XmlFileName); 4722 end 4723 4645 4724 4646 4725
Note: See TracChangeset
for help on using the changeset viewer.