Changeset 567 for trunk


Ignore:
Timestamp:
Jan 31, 2013, 2:24:00 PM (11 years ago)
Author:
sommeria
Message:

LIF calibration improved

Location:
trunk/src
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/transform_field/phys_polar.m

    r172 r567  
    4040%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    4141origin_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;
     42if isfield(CalibData,'PolarCentre') && isnumeric(CalibData.PolarCentre)
     43    if isequal(length(CalibData.PolarCentre),2);
     44        origin_xy= CalibData.PolarCentre;
    4545    end
    4646end
    4747radius_offset=0;%reference radius used to offset the radial coordinate r
    4848angle_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;
     49if isfield(CalibData,'PolarReferenceRadius') && isnumeric(CalibData.PolarReferenceRadius)
     50    radius_offset=CalibData.PolarReferenceRadius;
    5151end
    5252if radius_offset > 0
     
    5555    angle_scale=180/pi; %polar angle in degrees
    5656end
    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))
     57if 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))
    5959end
    6060% new x coordinate = radius-radius_offset;
     
    8080    end
    8181end
     82
    8283%transform second field (if exists) to cartesian phys coordiantes
    8384if test_1
     
    184185
    185186
    186 %------------------------------------------------------------------------
    187 %'phys_XYZ':transforms image (px) to real world (phys) coordinates using geometric calibration parameters
    188 % function [Xphys,Yphys]=phys_XYZ(Calib,X,Y,Z)
    189 %
    190 %OUTPUT:
    191 %
    192 %INPUT:
    193 %Z: index of plane
    194 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)>=Z
    197     Zindex=Z;
    198     Zphys=Calib.SliceCoord(Zindex,3);%GENERALISER AUX CAS AVEC ANGLE
    199 else
    200     Zphys=0;
    201 end
    202 if ~exist('X','var')||~exist('Y','var')
    203     Xphys=[];
    204     Yphys=[];%default
    205     return
    206 end
    207 %coordinate transform
    208 if ~isfield(Calib,'fx_fy')
    209      Calib.fx_fy=[1 1];
    210 end
    211 if ~isfield(Calib,'Tx_Ty_Tz')
    212      Calib.Tx_Ty_Tz=[0 0 1];
    213 end
    214 if ~isfield(Calib,'Cx_Cy')
    215      Calib.Cx_Cy=[0 0];
    216 end
    217 if ~isfield(Calib,'kc')
    218      Calib.kc=0;
    219 end
    220 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=1
    226     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)*Ty-R(5)*Tz+Z11*Zphys;
    237     A12=R(2)*Tz-R(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)*Tx-R(2)*Ty+Zx0*Zphys);
    241     Y0=f*(-R(4)*Tx+R(1)*Ty+Zy0*Zphys);
    242         %px to camera:
    243     Xd=dpx*(X-Calib.Cx_Cy(1)); % sensor coordinates
    244     Yd=(Y-Calib.Cx_Cy(2));
    245     dist_fact=1+Calib.kc*(Xd.*Xd+Yd.*Yd)/(f*f); %distortion factor
    246     Xu=Xd./dist_fact;%undistorted sensor coordinates
    247     Yu=Yd./dist_fact;
    248     denom=Dx*Xu+Dy*Yu+D0;
    249     Xphys=(A11.*Xu+A12.*Yu+X0)./denom;%world coordinates
    250     Yphys=(A21.*Xu+A22.*Yu+Y0)./denom;
    251 else
    252     Xphys=-Calib.Tx_Ty_Tz(1)+X/Calib.fx_fy(1);
    253     Yphys=-Calib.Tx_Ty_Tz(2)+Y/Calib.fx_fy(2);
    254 end
    255 
    256187%%%%%%%%%%%%%%%%%%%%
    257188function [A_out,Rangx,Rangy]=phys_Ima_polar(A,CalibIn,ZIndex,origin_xy,radius_offset,angle_offset,angle_scale)
  • trunk/src/uvmat.m

    r566 r567  
    45984598function MenuLIFCalib_Callback(hObject, eventdata, handles)
    45994599%------------------------------------------------------------------------
    4600 UvData=get(handles.uvmat,'UserData');%read UvData properties stored on the uvmat interface
     4600%% read UvData properties stored on the uvmat interface
     4601UvData=get(handles.uvmat,'UserData');
     4602if isfield(UvData,'XmlData')&& isfield(UvData.XmlData{1},'GeometryCalib')
     4603    XmlData=UvData.XmlData{1};
     4604else
     4605    msgbox_uvmat('ERROR','geometric calibration needed: use Tools/geometric calibration in the menu bar');
     4606    return
     4607end
     4608
     4609%% read lines currently drawn
    46014610ListObj=UvData.Object;
    46024611select=zeros(1,numel(ListObj));
    46034612for iobj=1:numel(ListObj);
    4604     if strcmp(ListObj{iobj}.Type,'line')
     4613    if isfield(ListObj{iobj},'Type') && strcmp(ListObj{iobj}.Type,'line')
    46054614        select(iobj)=1;
    46064615    end
     
    46114620    return
    46124621else
    4613     set(handles.ListObject,'Value',val);
     4622    set(handles.ListObject,'Value',val);% show the selected lines on the list
    46144623    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:npx-0.5;
    4619     yi=0.5:npy-0.5;
    4620     [Xi,Yi]=meshgrid(xi,yi);
    46214624    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
     4632end
     4633
     4634%% find the origin as intersection of the two first lines (see http://www.ahristov.com/tutorial/geometry-games/intersection-lines.html )
     4635x1=xA(1);x2=xB(1);
     4636x3=xA(2);x4=xB(2);
     4637y1=yA(1);y2=yB(1);
     4638y3=yA(2);y4=yB(2);
     4639D = (x1-x2)*(y3-y4) -(y1-y2)*(x3-x4);
     4640if D==0
     4641    msgbox_uvmat('ERROR','the two lines are parallel');
     4642    return
     4643end
     4644x0=((x3-x4)*(x1*y2-y1*x2)-(x1-x2)*(x3*y4-y3*x4))/D;
     4645y0=((y3-y4)*(x1*y2-y1*x2)-(y1-y2)*(x3*y4-y3*x4))/D;
     4646XmlData.Illumination.Origin=[x0 y0];
     4647XmlData.PolarCentre=[x0 y0];
     4648
     4649%% display the current image in polar coordinates with origin at the  illumination source
     4650currentdir=pwd; 
     4651uvmatpath=fileparts(which('uvmat'));
     4652cd(fullfile(uvmatpath,'transform_field'));
     4653phys_polar=str2func('phys_polar');
     4654cd(currentdir)
     4655DataOut=phys_polar(UvData.Field,XmlData);
     4656view_field(DataOut);
     4657
     4658%% use the third line for reference luminosity
     4659if 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_ref-x0;
     4663    y_ref=y_ref-y0;
     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_source-DataOut.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)
     4684end
     4685ImaName=regexprep([get(handles.RootFile,'String') get(handles.FileIndex,'String')],'//','');
     4686NewImageName=fullfile(get(handles.RootPath,'String'),'polar',[ImaName get(handles.FileExt,'String')]);
     4687imwrite(Anorm,NewImageName,'BitDepth',16)
     4688
     4689%% record the origin in the xml file
     4690XmlFileName=find_imadoc(get(handles.RootPath,'String'),get(handles.SubDir,'String'),get(handles.RootFile,'String'),get(handles.FileExt,'String'));
     4691answer=msgbox_uvmat('INPUT_Y-N','save the illumination origin in the current xml file?');
     4692if 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);
     4722end
     4723   
    46454724
    46464725
Note: See TracChangeset for help on using the changeset viewer.