Changeset 1046 for trunk/src


Ignore:
Timestamp:
May 30, 2018, 7:30:21 PM (6 years ago)
Author:
sommeria
Message:

read rrdvision corrected

Location:
trunk/src
Files:
2 edited

Legend:

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

    r1027 r1046  
    9595%parameters for polar coordinates (taken from the calibration data of the first field)
    9696%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    97 XmlData.PolarReferenceRadius=450;
    98 XmlData.PolarReferenceAngle=450*pi/2;
     97XmlData.PolarReferenceRadius=0;%450;
     98XmlData.PolarReferenceAngle=0;%450*pi/2;
    9999origin_xy=[0 0];%center for the polar coordinates in the original x,y coordinates
    100100radius_offset=0;%reference radius used to offset the radial coordinate r
    101101angle_offset=0; %reference angle used as new origin of the polar angle (= axis Ox by default)
     102angle_scale=180/pi;
    102103if isfield(XmlData,'TransformInput')
    103104    if isfield(XmlData.TransformInput,'PolarCentre') && isnumeric(XmlData.TransformInput.PolarCentre)
     
    249250npx=[];
    250251npy=[];
     252NbPoints=20; % nbre of points used to determine the image edge
    251253for icell=1:length(A)
    252254    siz=size(A{icell});
     
    255257    zphys=0; %default
    256258    if isfield(CalibIn{icell},'SliceCoord') %.Z= index of plane
     259        if ZIndex==0
     260            ZIndex=1;
     261        end
    257262       SliceCoord=CalibIn{icell}.SliceCoord(ZIndex,:);
    258263       zphys=SliceCoord(3); %to generalize for non-parallel planes
    259264    end
    260     xima=[0.5 siz(2)-0.5 0.5 siz(2)-0.5];%image coordiantes of corners
    261     yima=[0.5 0.5 siz(1)-0.5 siz(1)-0.5];
     265%     xima=[0.5 siz(2)-0.5 0.5 siz(2)-0.5];%image coordinates of corners
     266%     yima=[0.5 0.5 siz(1)-0.5 siz(1)-0.5];
     267      edge_x=linspace(0.5,siz(1)-0.5,NbPoints);     
     268      edge_y=linspace(0.5,siz(2)-0.5,NbPoints);
     269      xima=[edge_y (siz(2)-0.5)*ones(1,NbPoints) edge_y 0.5*ones(1,NbPoints)];%image coordinates of corners
     270      yima=[0.5*ones(1,NbPoints) edge_x (siz(1)-0.5)*ones(1,NbPoints) edge_x];%image coordinates of corners
    262271    [xcorner_new,ycorner_new]=phys_XYZ(CalibIn{icell},xima,yima,ZIndex);%corresponding physical coordinates
    263272    %transform the corner coordinates into polar ones   
  • trunk/src/uvmat.m

    r1045 r1046  
    13701370end
    13711371
     1372%% read the lines previously used for LIF calibration if available, and display them
     1373if isfield(XmlData,'LIFCalib')
     1374    if isfield(XmlData.LIFCalib,'Ray1Coord') && isfield(XmlData.LIFCalib,'Ray2Coord') && isfield(XmlData.LIFCalib,'RefLineCoord')
     1375        data.Name='Ray1';
     1376        data.Type='line';
     1377        data.ProjMode='none';%default
     1378        data.Coord=XmlData.LIFCalib.Ray1Coord;
     1379        UvData.ProjObject{2}=data;
     1380        UvData.ProjObject{2}.DisplayHandle.uvmat=plot_object(UvData.ProjObject{2},UvData.ProjObject{1},handles.PlotAxes,'b');
     1381        data.Name='Ray2';
     1382        data.Coord=XmlData.LIFCalib.Ray2Coord;
     1383        UvData.ProjObject{3}=data;
     1384        UvData.ProjObject{3}.DisplayHandle.uvmat=plot_object(UvData.ProjObject{3},UvData.ProjObject{1},handles.PlotAxes,'b');
     1385        data.Name='RefLine';
     1386        data.Coord=XmlData.LIFCalib.RefLineCoord;
     1387        UvData.ProjObject{4}=data;
     1388        UvData.ProjObject{4}.DisplayHandle.uvmat=plot_object(UvData.ProjObject{4},UvData.ProjObject{1},handles.PlotAxes,'b');
     1389    end
     1390    if isfield(XmlData.LIFCalib,'MaskPolygon')
     1391        data.Name='MaskPolygon';
     1392        data.Coord=XmlData.LIFCalib.MaskPolygonCoord;
     1393        UvData.ProjObject{5}=data;
     1394        UvData.ProjObject{5}.DisplayHandle.uvmat=plot_object(UvData.ProjObject{5},UvData.ProjObject{1},handles.PlotAxes,'b');
     1395    end
     1396end
     1397set(handles.ListObject,'String',{'plane';'Ray1';'Ray2';'RefLine';'MaskPolygon'})
     1398set(handles.ListObject,'Value',1)
     1399set(handles.uvmat,'UserData',UvData);
     1400
    13721401%% read lines currently drawn
    1373 ListObj=UvData.ProjObject;
     1402ListObj=UvData.ProjObject;% get the current list of projection objects
    13741403select=zeros(1,numel(ListObj));
     1404select_mask=zeros(1,numel(ListObj));
    13751405for iobj=1:numel(ListObj);
    13761406    if isfield(ListObj{iobj},'Type') && strcmp(ListObj{iobj}.Type,'line')
    1377         select(iobj)=1;
     1407        select(iobj)=1;% select the lines among the projection objects
     1408    end
     1409    if isfield(ListObj{iobj},'Type') && strcmp(ListObj{iobj}.Type,'polygon')
     1410        select_mask(iobj)=1;% select the mask among the projection objects
    13781411    end
    13791412end
    13801413val=find(select);
    13811414if numel(val)<2
    1382     msgbox_uvmat('ERROR','light rays must be defined by at least two lines created by Projection object/line in the menu bar');
     1415    msgbox_uvmat('ERROR',{'light rays must be defined by at least two lines created by Projection object/line in the menu bar'; ...
     1416    'use a third line to get a reference luminosity profile accross the illumination beam'});
    13831417    return
    13841418else
    1385     set(handles.ListObject,'Value',val);% show the selected lines on the list
    13861419    ObjectData=UvData.ProjObject(val);
    13871420    for iobj=1:length(ObjectData)
     
    13911424            yB(iobj)=ObjectData{iobj}.Coord(2,2);
    13921425    end
     1426end
     1427index_mask=find(select_mask);
     1428if numel(index_mask)>=1   
     1429    MaskData=UvData.ProjObject(index_mask);
     1430    for iobj=1:length(MaskData)
     1431           
     1432    end
     1433end
     1434
     1435%% set the image offset
     1436blackoffset=0;
     1437if isfield(XmlData,'LIFCalib')&& isfield(XmlData.LIFCalib,'BlackOffset')
     1438    blackoffset=XmlData.LIFCalib.BlackOffset ;% image value for black background, to be determined by taking images with a cover on the objective lens
     1439end
     1440answer=msgbox_uvmat('INPUT_TXT','camera offset value in the absence of illumination:', num2str(blackoffset));
     1441if strcmp(answer,'Cancel')
     1442    return
     1443else
     1444    XmlData.LIFCalib.BlackOffset=answer ;% image value for black background, to be determined by taking images with a cover on the objective lens
    13931445end
    13941446
     
    14051457x0=((x3-x4)*(x1*y2-y1*x2)-(x1-x2)*(x3*y4-y3*x4))/D;
    14061458y0=((y3-y4)*(x1*y2-y1*x2)-(y1-y2)*(x3*y4-y3*x4))/D;
    1407 XmlData.Illumination.Origin=[x0 y0];
    1408 XmlData.PolarCentre=[x0 y0];
     1459XmlData.LIFCalib.LightOrigin=[x0 y0];% origin of the light source to be saved in the current xml file
     1460XmlData.LIFCalib.Ray1Coord=ObjectData{1}.Coord;
     1461XmlData.LIFCalib.Ray2Coord=ObjectData{2}.Coord;
     1462XmlData.LIFCalib.RefLineCoord=ObjectData{3}.Coord;
    14091463
    14101464%% display the current image in polar axes with origin at the  illumination source
     
    14141468phys_polar=str2func('phys_polar');
    14151469cd(currentdir)
    1416 DataOut=phys_polar(UvData.Field,XmlData);
    1417 view_field(DataOut);
    1418 
    1419 %% use the third line for reference luminosity
     1470XmlData.TransformInput.PolarCentre=[x0 y0];
     1471DataPol=phys_polar(UvData.Field,XmlData);%transform the input image in polar coordinates with origin at the light source
     1472
     1473%% use the third line for reference luminosity, renormalize the image intensity along each ray to get a uniform brightness along this line
    14201474if numel(val)==3
    14211475    x_ref=linspace(ObjectData{3}.Coord(1,1),ObjectData{3}.Coord(2,1),10);
     
    14241478    y_ref=y_ref-y0;
    14251479    [theta_ref,r_ref] = cart2pol(x_ref,y_ref);%theta_ref  and r_ref are the polar coordinates of the points on the line
    1426     theta_ref=theta_ref*180/pi;
    1427     figure
     1480    theta_ref=theta_ref*180/pi;% theta_ref in radians
     1481    figure(10)
    14281482    plot(theta_ref,r_ref)
    1429     azimuth_ima=linspace(DataOut.Coord_y(1),DataOut.Coord_y(2),size(DataOut.A,1));%profile of x index on the transformed image
    1430     dist_source = interp1(theta_ref,r_ref,azimuth_ima);
    1431     dist_source_pixel=round(size(DataOut.A,2)*(dist_source-DataOut.Coord_x(1))/(DataOut.Coord_x(2)-DataOut.Coord_x(1)));
    1432     line_nan= isnan(dist_source_pixel);
     1483    xlabel('azimuth(degrees)')
     1484    ylabel('radius from light source')
     1485    title('ref line in polar coordinates')
     1486    azimuth_ima=linspace(DataPol.Coord_y(1),DataPol.Coord_y(2),size(DataPol.A,1));%array of angular index on the transformed image
     1487    dist_source = interp1(theta_ref,r_ref,azimuth_ima);% get the polar position of the reference line
     1488    dist_source_pixel=round(size(DataPol.A,2)*(dist_source-DataPol.Coord_x(1))/(DataPol.Coord_x(2)-DataPol.Coord_x(1)));
     1489    line_nan= isnan(dist_source);
    14331490    dist_source_pixel(line_nan)=1;
    14341491    width=20; %number of pixels used for reference
    1435     DataOut.A=double(DataOut.A);
    1436     Anorm=zeros(size(DataOut.A));
    1437     Aval=mean(mean(DataOut.A));
    1438     for iline=1:size(DataOut.A,1)
    1439         lum(iline)=mean(DataOut.A(iline,dist_source_pixel(iline):dist_source_pixel(iline)+width));
    1440         Anorm(iline,:)=uint16(Aval*DataOut.A(iline,:)/lum(iline));
     1492    DataPol.A=double(DataPol.A)-XmlData.LIFCalib.BlackOffset;% black background substracted
     1493    Anorm=zeros(size(DataPol.A));
     1494    for iline=1:size(DataPol.A,1)
     1495        lum(iline)=mean(DataPol.A(iline,dist_source_pixel(iline):dist_source_pixel(iline)+width));% average the luminosity on a band width lying on the reference line
     1496        Anorm(iline,:)=DataPol.A(iline,:)/lum(iline);% for each ray (iline), renormalise the image by the brightness at the reference line
    14411497    end
    14421498    lum(line_nan)=NaN;
    1443     figure
    1444     plot(1:size(DataOut.A,1),lum)
    1445 end
    1446 ImaName=regexprep([get(handles.RootFile,'String') get(handles.FileIndex,'String')],'//','');
    1447 NewImageName=fullfile(get(handles.RootPath,'String'),'polar',[ImaName get(handles.FileExt,'String')]);
    1448 imwrite(Anorm,NewImageName,'BitDepth',16)
    1449 
    1450 %% record the origin in the xml file
     1499    figure(11)
     1500    plot(1:size(DataPol.A,1),lum)
     1501    xlabel('azimuth (pixels)')
     1502    ylabel('image brightness')
     1503    %ImaName=regexprep([get(handles.RootFile,'String') get(handles.FileIndex,'String')],'//','');
     1504NewImageName=fullfile(get(handles.RootPath,'String'),'LIF_ref.png');
     1505imwrite(uint16(1000*Anorm),NewImageName,'BitDepth',16)
     1506DataPol.A=uint16(1000*Anorm);
     1507view_field(DataPol);
     1508end
     1509
     1510% figure(12)
     1511% imagesc(DataPol.Coord_x,DataPol.Coord_y,AY,uint16(1000*Anorm))
     1512%title('rescaled image')
     1513
     1514%% record the calibration data in the xml file
     1515
    14511516XmlFileName=find_imadoc(get(handles.RootPath,'String'),get(handles.SubDir,'String'),get(handles.RootFile,'String'),get(handles.FileExt,'String'));
    14521517answer=msgbox_uvmat('INPUT_Y-N','save the illumination origin in the current xml file?');
    14531518if strcmp(answer,'Yes')
    14541519    t=xmltree(XmlFileName); %read the file
    1455     title=get(t,1,'name');
    1456     if ~strcmp(title,'ImaDoc')
     1520    title_str=get(t,1,'name');
     1521    if ~strcmp(title_str,'ImaDoc')
    14571522        msgbox_uvmat('ERROR','wrong xml file');
    14581523        return
     
    14701535        return
    14711536    end
    1472     uid_illumination=find(t,'ImaDoc/Illumination');
     1537    uid_illumination=find(t,'ImaDoc/LIFCalib');
    14731538    if isempty(uid_illumination)  %if GeometryCalib does not already exists, create it
    1474         [t,uid_illumination]=add(t,1,'element','Illumination');
    1475     end
    1476     uid_origin=find(t,'ImaDoc/Illumination/Origin');
     1539        [t,uid_illumination]=add(t,1,'element','LIFCalib');
     1540    end
     1541    uid_origin=find(t,'ImaDoc/LIFCalib/LightOrigin');
    14771542    if ~isempty(uid_origin)  %if GeometryCalib does not already exists, create it
    14781543         t=delete(t,uid_origin);
    14791544    end
    1480     % save the illumination origin
    1481     t=struct2xml(XmlData.Illumination,t,uid_illumination);
     1545    uid_line=find(t,'ImaDoc/LIFCalib/Ray1Coord');
     1546    if ~isempty(uid_line)  %if GeometryCalib does not already exists, create it
     1547         t=delete(t,uid_line);
     1548    end
     1549       uid_line=find(t,'ImaDoc/LIFCalib/Ray2Coord');
     1550    if ~isempty(uid_line)  %if GeometryCalib does not already exists, create it
     1551         t=delete(t,uid_line);
     1552    end
     1553       uid_line=find(t,'ImaDoc/LIFCalib/RefLineCoord');
     1554    if ~isempty(uid_line)  %if GeometryCalib does not already exists, create it
     1555         t=delete(t,uid_line);
     1556    end
     1557     uid_BlackOffset=find(t,'ImaDoc/LIFCalib/BlackOffset');
     1558    if ~isempty(uid_BlackOffset)  %if GeometryCalib does not already exists, create it
     1559         t=delete(t,uid_BlackOffset);
     1560    end
     1561    % save the LIF calibration data
     1562    t=struct2xml(XmlData.LIFCalib,t,uid_illumination);
    14821563    save(t,XmlFileName);
    14831564end
     1565
     1566
     1567%% get the image A*r
     1568[npy,npx]=size(DataPol.A);
     1569AX=DataPol.Coord_x;
     1570AY=DataPol.Coord_y;
     1571dX=(AX(2)-AX(1))/(npx-1);
     1572dY=(AY(1)-AY(2))/(npy-1);%mesh of new pixels
     1573r=AX(1)+[0:npx-1]*dX;%distance from laser
     1574rmat=ones(npy,1)*r;
     1575A=rmat.*Anorm;
     1576A=(A<=0).*ones(npy,npx)+(A>0).*A; %replaces zeros by ones
     1577A=log(A); % rA(r) should have a slope -gamma in x (radius from laser source)
     1578
     1579%%  find mask image: for use in uvmat
     1580% if ~isempty(DataIn.ZIndex)
     1581% num_level=DataIn.ZIndex;% needed to handle multilevel illumination
     1582% if ~exist('maskname','var')
     1583%     huvmat=findobj(allchild(0),'tag','uvmat');
     1584%     hhuvmat=guidata(huvmat);
     1585%     RootPath=get(hhuvmat.RootPath,'String');
     1586%     RootFile=get(hhuvmat.RootFile,'String');
     1587%     maskname=fullfile(RootPath,[RootFile '_23mask_' num2str(num_level) '.png']);
     1588%     if test_1% PIV image used for correction
     1589%         RootPath_1=get(hhuvmat.RootPath_1,'String');
     1590%         RootFile_1=get(hhuvmat.RootFile_1,'String');
     1591%         maskname_1=fullfile(RootPath_1,[RootFile_1 '_23mask_' num2str(num_level) '.png']);
     1592%     end
     1593% end
     1594% Mask=DataIn;
     1595% if test_1
     1596%     Mask_1=DataIn_1;
     1597%     Mask.A=imread(maskname);
     1598%     Mask_1.A=imread(maskname_1);
     1599%     [Mask,Mask_1]=phys_polar(Mask,XmlData,Mask_1,XmlData_1);
     1600%     Mask.A=Mask.A >= 200 & Mask_1.A >= 200 ;
     1601% % end
     1602%  end
     1603% Mask.A=ones(size(DataPol.A));
     1604[R,Y]=meshgrid(linspace(AX(1),AX(2),npx),linspace(AY(1),AY(2),npy));
     1605[X,Y] = pol2cart(Y*pi/180,R);
     1606X=X+x0;% ny*nx matrix of phys cartesian coordiantes
     1607Y=Y+y0;
     1608gamma_coeff=zeros(1,npy);
     1609
     1610NewImageName=fullfile(get(handles.RootPath,'String'),'LIF_ref_rescaled.png');
     1611imwrite(uint16(1000*A),NewImageName,'BitDepth',16)
     1612return
     1613%loop on lines iY (angle in polar coordiantes)
     1614for iY=1:npy
     1615    Theta_i=(AY(1)-(iY-1)*dY)*pi/180;%theta in radians
     1616    %r_edge(iY)=(-100-XmlData.GeometryCalib.PolarCentre(1))/cos(Theta_i);%r_edge=radius of the edge (at X=-95)
     1617    r_edge(iY)
     1618    Aline=A(iY,:);% =r*image intensity
     1619    Xline=X(iY,:);
     1620    % adjust the fit region
     1621    test_fit=X(iY,:)>10 & X(iY,:)<70 & A(iY,:)>0 & Mask.A(iY,:);
     1622    ind_fit=find(test_fit);% indices of points used for the exponential fit
     1623    if numel(ind_fit)<50
     1624        Mask.A(iY,:)=0;% to exclude lines for which the line of rays is too short
     1625    else
     1626        p = polyfit(r(ind_fit),Aline(ind_fit),1);
     1627        Aline_fit=polyval(p,r);
     1628        if iY==round(npy/2)
     1629            figure(12)
     1630            plot(r(ind_fit),Aline(ind_fit),'b+',r,Aline_fit,'r')
     1631            test_out=(Aline_fit-Aline)>0.5
     1632            ind_out=find(test_out & test_fit)
     1633        end
     1634        test_out=(Aline_fit-Aline)>0.5;% eliminate points too dark compared to the fit, such that log(Afit/A)>0.5
     1635        ind_out=find(test_out & test_fit);
     1636%         if ~isempty(find(test_out & test_fit,1))
     1637%             Mask.A(iY,test_out)=0;
     1638         if ~isempty(ind_out)
     1639              Mask.A(iY,ind_out)=0;
     1640            ind_fit=find(~test_out & test_fit);
     1641            p = polyfit(r(ind_fit),Aline(ind_fit),1);%do new fit after elimination of the dark points
     1642        end
     1643        if -p(1) > 0.003 && -p(1) < 0.012
     1644            gamma_coeff(iY)=-p(1);
     1645            Aline_edge=polyval(p,r_edge(iY));
     1646            DataPol.A(iY,:)=exp(Aline_edge)/(r_edge(iY));
     1647        end
     1648
     1649    end
     1650end
     1651gamma_select=gamma_coeff(gamma_coeff~=0);
     1652gamma_select=filter(ones(1,20)/20,1,gamma_select);%smooth gamma_ceff over 20 pixels
     1653gamma_coeff(gamma_coeff~=0)=gamma_select;
     1654DataPol.A(Mask.A==0)=0;
     1655Ref.r_edge=r_edge;
     1656Ref.GammaCoeff=gamma_coeff
     1657
     1658
     1659
     1660
     1661
     1662
     1663
     1664
     1665
     1666
    14841667
    14851668
     
    20942277                end
    20952278            end
     2279        end
     2280        if isfield(XmlDataRead, 'LIFCalib')
     2281            XmlData.LIFCalib=XmlDataRead.LIFCalib;
    20962282        end
    20972283    end
Note: See TracChangeset for help on using the changeset viewer.