Changeset 1049 for trunk


Ignore:
Timestamp:
Jun 5, 2018, 10:13:08 PM (3 years ago)
Author:
sommeria
Message:

LIF clib corrected + histo velocity repaired

Location:
trunk/src
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/find_field_cells.m

    r1048 r1049  
    475475        for ivar=ind_coord_y
    476476            DimCell=Data.VarDimName{ivar};
    477             if  numel(DimCell)==1 && strcmp(DimCell_x{1},DimCell{1})
     477            if  numel(DimCell_x)==1 && strcmp(DimCell_x{1},DimCell{1})
    478478                y_nbre(icell)=y_nbre(icell)+1;
    479479                Cell1DPlot{icell}.YIndex(y_nbre(icell))=ivar;
  • trunk/src/series/extract_multitif_special.m

    r1048 r1049  
    144144% end
    145145%ImagesPerLevel=size(XmlInput.Time,2)-1;%100;%use the xmlinformation to get the nbre of j indices
     146Dtj=0.05;% time interval between frames
    146147ImagesPerLevel=450;% total number of images per position, ImagesPerLevel-Nbj images skiiped during motion between two positions
    147148Nbj=400; %Nbre of images kept at a given position
     149Dti=Dtj*ImagesPerLevel;
     150NbLevel=11;
     151NbScan=3;
     152TimeReturn=20; %time needed to return back to the first position (in sec)
     153NbSkippedReturn=round(TimeReturn/Dtj);
    148154%% create the xml file of PCO camera if it does not exist
    149155Newxml=fullfile(Param.InputTable{1,1},[Param.InputTable{1,2} '.xml']);
     
    154160   % XmlInput.Camera.BurstTiming.Time=-1;% for 180
    155161   XmlInput.Camera.BurstTiming.Time=0;% for 200
    156     XmlInput.Camera.BurstTiming.Dtj=0.05;
    157     XmlInput.Camera.BurstTiming.NbDtj=199;
    158     XmlInput.Camera.BurstTiming.Dti=12.5;
    159     XmlInput.Camera.BurstTiming.NbDti=10;
    160     XmlInput.Camera.BurstTiming.Dtk=157;
    161     XmlInput.Camera.BurstTiming.NbDtk=2;
     162    XmlInput.Camera.BurstTiming.Dtj=Dtj;
     163    XmlInput.Camera.BurstTiming.NbDtj=Nbj-1;
     164    XmlInput.Camera.BurstTiming.Dti=Dti;
     165    XmlInput.Camera.BurstTiming.NbDti=NbLevel-1;
     166    XmlInput.Camera.BurstTiming.Dtk=Dti*Nb
     167    Level+TimeReturn;
     168    XmlInput.Camera.BurstTiming.NbDtk=NbScan-1;
    162169    %XmlInput=rmfield(XmlInput,'Time');
    163170    %XmlInput=rmfield(XmlInput,'TimeUnit');
     
    198205            i_index=fix((count-1)/ImagesPerLevel)+1;
    199206            j_index=mod(count-1,ImagesPerLevel)+1;
    200         elseif count<=ImagesPerLevel*11+400 %skip 400 images during return
     207        elseif count<=ImagesPerLevel*NbLevel+400 %skip 400 images during return
    201208            checkkeep=0;
    202         elseif count<=ImagesPerLevel*22+400 % =5930 second scan
     209        elseif count<=ImagesPerLevel*2*NbLevel+400 % =5930 second scan
    203210            i_index=fix((count-401)/ImagesPerLevel)+1;
    204211            j_index=mod(count-401,ImagesPerLevel)+1;
    205         elseif count<=ImagesPerLevel*22+800 %skip images during second return, from 2763 to 3167
     212        elseif count<=ImagesPerLevel*2*NbLevel+800 %skip images during second return, from 2763 to 3167
    206213            checkkeep=0;
    207         elseif count<=ImagesPerLevel*33+800 % =5930 third scan
     214        elseif count<=ImagesPerLevel*3*NbLevel+800 % =5930 third scan
    208215            i_index=fix((count-801)/ImagesPerLevel)+1;
    209216            j_index=mod(count-801,ImagesPerLevel)+1;
  • trunk/src/uvmat.m

    r1048 r1049  
    13711371
    13721372%% read the lines previously used for LIF calibration if available, and display them
     1373ListObjectName=get(handles.ListObject,'String');
    13731374if isfield(XmlData,'LIFCalib')
    13741375    if isfield(XmlData.LIFCalib,'Ray1Coord') && isfield(XmlData.LIFCalib,'Ray2Coord') && isfield(XmlData.LIFCalib,'RefLineCoord')
     1376        ListObjectName(1:4)={'plane';'Ray1';'Ray2';'RefLine'};
    13751377        data.Name='Ray1';
    13761378        data.Type='line';
     
    13881390        UvData.ProjObject{4}.DisplayHandle.uvmat=plot_object(UvData.ProjObject{4},UvData.ProjObject{1},handles.PlotAxes,'b');
    13891391    end
    1390     if isfield(XmlData.LIFCalib,'MaskPolygon')
     1392    if isfield(XmlData.LIFCalib,'MaskPolygonCoord')
     1393        ListObjectName{5}='MaskPolygon';
    13911394        data.Name='MaskPolygon';
    13921395        data.Coord=XmlData.LIFCalib.MaskPolygonCoord;
    13931396        UvData.ProjObject{5}=data;
    13941397        UvData.ProjObject{5}.DisplayHandle.uvmat=plot_object(UvData.ProjObject{5},UvData.ProjObject{1},handles.PlotAxes,'b');
    1395     end
    1396 end
    1397 set(handles.ListObject,'String',{'plane';'Ray1';'Ray2';'RefLine';'MaskPolygon'})
     1398    end     
     1399end
     1400set(handles.ListObject,'String',ListObjectName)
    13981401set(handles.ListObject,'Value',1)
    13991402set(handles.uvmat,'UserData',UvData);
     
    14011404%% read lines currently drawn
    14021405ListObj=UvData.ProjObject;% get the current list of projection objects
    1403 select=zeros(1,numel(ListObj));
     1406select_line=zeros(1,numel(ListObj));
    14041407select_mask=zeros(1,numel(ListObj));
    14051408for iobj=1:numel(ListObj);
    14061409    if isfield(ListObj{iobj},'Type') && strcmp(ListObj{iobj}.Type,'line')
    1407         select(iobj)=1;% select the lines among the projection objects
     1410        select_line(iobj)=1;% select the lines among the projection objects
    14081411    end
    14091412    if isfield(ListObj{iobj},'Type') && strcmp(ListObj{iobj}.Type,'polygon')
     
    14111414    end
    14121415end
    1413 val=find(select);
    1414 if numel(val)<2
     1416if numel(find(select_line))<2
    14151417    msgbox_uvmat('ERROR',{'light rays must be defined by at least two lines created by Projection object/line in the menu bar'; ...
    14161418    'use a third line to get a reference luminosity profile accross the illumination beam'});
    14171419    return
    14181420else
    1419     ObjectData=UvData.ProjObject(val);
    1420     for iobj=1:length(ObjectData)
    1421             xA(iobj)=ObjectData{iobj}.Coord(1,1);
    1422             yA(iobj)=ObjectData{iobj}.Coord(1,2);
    1423             xB(iobj)=ObjectData{iobj}.Coord(2,1);
    1424             yB(iobj)=ObjectData{iobj}.Coord(2,2);
    1425     end
    1426 end
    1427 index_mask=find(select_mask);
    1428 if numel(index_mask)>=1   
    1429     MaskData=UvData.ProjObject(index_mask);
    1430     for iobj=1:length(MaskData)
    1431            
    1432     end
    1433 end
     1421    LineData=UvData.ProjObject(find(select_line));
     1422    for iobj=1:length(LineData)
     1423            xA(iobj)=LineData{iobj}.Coord(1,1);
     1424            yA(iobj)=LineData{iobj}.Coord(1,2);
     1425            xB(iobj)=LineData{iobj}.Coord(2,1);
     1426            yB(iobj)=LineData{iobj}.Coord(2,2);
     1427    end
     1428end
     1429
    14341430
    14351431%% set the image offset
     
    14421438    return
    14431439else
    1444     XmlData.LIFCalib.BlackOffset=answer ;% image value for black background, to be determined by taking images with a cover on the objective lens
     1440    XmlData.LIFCalib.BlackOffset=str2num(answer) ;% image value for black background, to be determined by taking images with a cover on the objective lens
    14451441end
    14461442
     
    14581454y0=((y3-y4)*(x1*y2-y1*x2)-(y1-y2)*(x3*y4-y3*x4))/D;
    14591455XmlData.LIFCalib.LightOrigin=[x0 y0];% origin of the light source to be saved in the current xml file
    1460 XmlData.LIFCalib.Ray1Coord=ObjectData{1}.Coord;
    1461 XmlData.LIFCalib.Ray2Coord=ObjectData{2}.Coord;
    1462 XmlData.LIFCalib.RefLineCoord=ObjectData{3}.Coord;
     1456XmlData.LIFCalib.Ray1Coord=LineData{1}.Coord;
     1457XmlData.LIFCalib.Ray2Coord=LineData{2}.Coord;
     1458XmlData.LIFCalib.RefLineCoord=LineData{3}.Coord;
    14631459
    14641460%% display the current image in polar axes with origin at the  illumination source
     
    14721468
    14731469%% use the third line for reference luminosity, renormalize the image intensity along each ray to get a uniform brightness along this line
    1474 if numel(val)==3
    1475     x_ref=linspace(ObjectData{3}.Coord(1,1),ObjectData{3}.Coord(2,1),10);
    1476     y_ref=linspace(ObjectData{3}.Coord(1,2),ObjectData{3}.Coord(2,2),10);
     1470if numel(find(select_line))==3
     1471    x_ref=linspace(LineData{3}.Coord(1,1),LineData{3}.Coord(2,1),10);
     1472    y_ref=linspace(LineData{3}.Coord(1,2),LineData{3}.Coord(2,2),10);
    14771473    x_ref=x_ref-x0;
    14781474    y_ref=y_ref-y0;
     
    15011497    xlabel('azimuth (pixels)')
    15021498    ylabel('image brightness')
    1503     %ImaName=regexprep([get(handles.RootFile,'String') get(handles.FileIndex,'String')],'//','');
    1504 NewImageName=fullfile(get(handles.RootPath,'String'),'LIF_ref.png');
    1505 imwrite(uint16(1000*Anorm),NewImageName,'BitDepth',16)
    1506 DataPol.A=uint16(1000*Anorm);
     1499    title('ref brightness profile')
     1500end
     1501
     1502%% get the image A*r
     1503[npy,npx]=size(Anorm);
     1504AX=DataPol.Coord_x;
     1505AY=DataPol.Coord_y;
     1506dX=(AX(2)-AX(1))/(npx-1);
     1507dY=(AY(1)-AY(2))/(npy-1);%mesh of new pixels
     1508[R,Theta]=meshgrid(linspace(AX(1),AX(2),npx),linspace(AY(1),AY(2),npy));
     1509A=R.*Anorm;
     1510A=(A<=0).*ones(npy,npx)+(A>0).*A; %replaces zeros by ones
     1511A=log(A); % rA(r) should have a slope -gamma in x (radius from laser source)
     1512
     1513%% get the exponential decay law
     1514index_mask=find(select_mask);
     1515if ~isempty(index_mask)
     1516    if numel(index_mask)>1
     1517        msgbox_uvmat('ERROR',{'only one mask polygon accepted'});
     1518        return
     1519    else
     1520        MaskData=UvData.ProjObject{index_mask};
     1521        XmlData.LIFCalib.MaskPolygonCoord=MaskData.Coord;
     1522    end
     1523end
     1524
     1525%loop on lines iY (angle in polar coordiantes)
     1526gamma_coeff=zeros(1,npy);
     1527for iY=1:npy
     1528    ALine=A(iY,:);
     1529    RLine=R(iY,:);
     1530    ThetaLine=Theta(iY,:)*pi/180;
     1531    [XLine,YLine] = pol2cart(ThetaLine,RLine);
     1532    XLine=XLine+x0;
     1533    YLine=YLine+y0; 
     1534    if ~isempty(index_mask)
     1535        ind_good=inpolygon(XLine,YLine,MaskData.Coord(:,1),MaskData.Coord(:,2));
     1536        if numel(find(ind_good))>100
     1537        ALine=ALine(ind_good);
     1538        RLine=RLine(ind_good);
     1539        else
     1540            continue
     1541        end
     1542    end
     1543    p = polyfit(RLine,ALine,1);
     1544        gamma_coeff(iY)=-p(1);
     1545end
     1546
     1547%% show and store the rescaled image
     1548NewImageName=fullfile(get(handles.RootPath,'String'),'LIF_ref_rescaled.png');
     1549DataPol.A=uint16(1000*A);
    15071550view_field(DataPol);
    1508 end
    1509 
    1510 % figure(12)
    1511 % imagesc(DataPol.Coord_x,DataPol.Coord_y,AY,uint16(1000*Anorm))
    1512 %title('rescaled image')
     1551imwrite(DataPol.A,NewImageName,'BitDepth',16)
     1552figure(13)
     1553plot(Theta(:,1),100*gamma_coeff)
     1554xlabel('angle(degree)')
     1555ylabel('decay coeff(m-1)')
    15131556
    15141557%% record the calibration data in the xml file
    1515 
    15161558XmlFileName=find_imadoc(get(handles.RootPath,'String'),get(handles.SubDir,'String'),get(handles.RootFile,'String'),get(handles.FileExt,'String'));
    15171559answer=msgbox_uvmat('INPUT_Y-N','save the illumination origin in the current xml file?');
     
    15511593         t=delete(t,uid_line);
    15521594    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);
     1595       uid_mask=find(t,'ImaDoc/LIFCalib/MaskPolygonCoord');
     1596    if ~isempty(uid_mask)  %if GeometryCalib does not already exists, create it
     1597         t=delete(t,uid_mask);
    15561598    end
    15571599     uid_BlackOffset=find(t,'ImaDoc/LIFCalib/BlackOffset');
     
    15631605    save(t,XmlFileName);
    15641606end
    1565 
    1566 
    1567 %% get the image A*r
    1568 [npy,npx]=size(DataPol.A);
    1569 AX=DataPol.Coord_x;
    1570 AY=DataPol.Coord_y;
    1571 dX=(AX(2)-AX(1))/(npx-1);
    1572 dY=(AY(1)-AY(2))/(npy-1);%mesh of new pixels
    1573 r=AX(1)+[0:npx-1]*dX;%distance from laser
    1574 rmat=ones(npy,1)*r;
    1575 A=rmat.*Anorm;
    1576 A=(A<=0).*ones(npy,npx)+(A>0).*A; %replaces zeros by ones
    1577 A=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);
    1606 X=X+x0;% ny*nx matrix of phys cartesian coordiantes
    1607 Y=Y+y0;
    1608 gamma_coeff=zeros(1,npy);
    1609 
    1610 NewImageName=fullfile(get(handles.RootPath,'String'),'LIF_ref_rescaled.png');
    1611 imwrite(uint16(1000*A),NewImageName,'BitDepth',16)
    1612 return
    1613 %loop on lines iY (angle in polar coordiantes)
    1614 for 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
    1650 end
    1651 gamma_select=gamma_coeff(gamma_coeff~=0);
    1652 gamma_select=filter(ones(1,20)/20,1,gamma_select);%smooth gamma_ceff over 20 pixels
    1653 gamma_coeff(gamma_coeff~=0)=gamma_select;
    1654 DataPol.A(Mask.A==0)=0;
    1655 Ref.r_edge=r_edge;
    1656 Ref.GammaCoeff=gamma_coeff
    1657 
    1658 
    1659 
    1660 
    1661 
    1662 
    1663 
    16641607
    16651608
Note: See TracChangeset for help on using the changeset viewer.