Changeset 567

Timestamp:

Jan 31, 2013, 2:24:00 PM (12 years ago)

Author:

sommeria

Message:

LIF calibration improved

Location:

trunk/src

Files:

• transform_field/phys_polar.m (modified) (4 diffs)
• uvmat.m (modified) (2 diffs)

trunk/src/transform_field/phys_polar.m

@@ -40 +40 @@
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 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 +55 @@
     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 +80 @@
     end
 end
+
 %transform second field (if exists) to cartesian phys coordiantes
 if test_1
@@ -184 +185 @@
 
 
-%------------------------------------------------------------------------
-%'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)

trunk/src/uvmat.m

@@ -4598 +4598 @@
 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 +4620 @@
     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
+