Changeset 567 for trunk/src/uvmat.m


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

LIF calibration improved

File:
1 edited

Legend:

Unmodified
Added
Removed
  • 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.