Changeset 1049 for trunk/src

Timestamp:

Jun 5, 2018, 10:13:08 PM (7 years ago)

Author:

sommeria

Message:

LIF clib corrected + histo velocity repaired

Location:

trunk/src

Files:

• find_field_cells.m (modified) (1 diff)
• series/extract_multitif_special.m (modified) (3 diffs)
• uvmat.m (modified) (10 diffs)

trunk/src/find_field_cells.m

@@ -475 +475 @@
         for ivar=ind_coord_y 
             DimCell=Data.VarDimName{ivar};
-            if  numel(DimCell)==1 && strcmp(DimCell_x{1},DimCell{1})
+            if  numel(DimCell_x)==1 && strcmp(DimCell_x{1},DimCell{1})
                 y_nbre(icell)=y_nbre(icell)+1;
                 Cell1DPlot{icell}.YIndex(y_nbre(icell))=ivar;

trunk/src/series/extract_multitif_special.m

@@ -144 +144 @@
 % end
 %ImagesPerLevel=size(XmlInput.Time,2)-1;%100;%use the xmlinformation to get the nbre of j indices
+Dtj=0.05;% time interval between frames
 ImagesPerLevel=450;% total number of images per position, ImagesPerLevel-Nbj images skiiped during motion between two positions
 Nbj=400; %Nbre of images kept at a given position
+Dti=Dtj*ImagesPerLevel;
+NbLevel=11;
+NbScan=3;
+TimeReturn=20; %time needed to return back to the first position (in sec)
+NbSkippedReturn=round(TimeReturn/Dtj);
 %% create the xml file of PCO camera if it does not exist
 Newxml=fullfile(Param.InputTable{1,1},[Param.InputTable{1,2} '.xml']);
@@ -154 +160 @@
    % XmlInput.Camera.BurstTiming.Time=-1;% for 180
    XmlInput.Camera.BurstTiming.Time=0;% for 200
-    XmlInput.Camera.BurstTiming.Dtj=0.05;
-    XmlInput.Camera.BurstTiming.NbDtj=199;
-    XmlInput.Camera.BurstTiming.Dti=12.5;
-    XmlInput.Camera.BurstTiming.NbDti=10;
-    XmlInput.Camera.BurstTiming.Dtk=157;
-    XmlInput.Camera.BurstTiming.NbDtk=2;
+    XmlInput.Camera.BurstTiming.Dtj=Dtj;
+    XmlInput.Camera.BurstTiming.NbDtj=Nbj-1;
+    XmlInput.Camera.BurstTiming.Dti=Dti;
+    XmlInput.Camera.BurstTiming.NbDti=NbLevel-1;
+    XmlInput.Camera.BurstTiming.Dtk=Dti*Nb
+    Level+TimeReturn;
+    XmlInput.Camera.BurstTiming.NbDtk=NbScan-1;
     %XmlInput=rmfield(XmlInput,'Time');
     %XmlInput=rmfield(XmlInput,'TimeUnit');
@@ -198 +205 @@
             i_index=fix((count-1)/ImagesPerLevel)+1;
             j_index=mod(count-1,ImagesPerLevel)+1;
-        elseif count<=ImagesPerLevel*11+400 %skip 400 images during return
+        elseif count<=ImagesPerLevel*NbLevel+400 %skip 400 images during return
             checkkeep=0;
-        elseif count<=ImagesPerLevel*22+400 % =5930 second scan
+        elseif count<=ImagesPerLevel*2*NbLevel+400 % =5930 second scan
             i_index=fix((count-401)/ImagesPerLevel)+1;
             j_index=mod(count-401,ImagesPerLevel)+1;
-        elseif count<=ImagesPerLevel*22+800 %skip images during second return, from 2763 to 3167
+        elseif count<=ImagesPerLevel*2*NbLevel+800 %skip images during second return, from 2763 to 3167
             checkkeep=0;
-        elseif count<=ImagesPerLevel*33+800 % =5930 third scan
+        elseif count<=ImagesPerLevel*3*NbLevel+800 % =5930 third scan
             i_index=fix((count-801)/ImagesPerLevel)+1;
             j_index=mod(count-801,ImagesPerLevel)+1;

trunk/src/uvmat.m

@@ -1371 +1371 @@
 
 %% read the lines previously used for LIF calibration if available, and display them
+ListObjectName=get(handles.ListObject,'String');
 if isfield(XmlData,'LIFCalib')
     if isfield(XmlData.LIFCalib,'Ray1Coord') && isfield(XmlData.LIFCalib,'Ray2Coord') && isfield(XmlData.LIFCalib,'RefLineCoord')
+        ListObjectName(1:4)={'plane';'Ray1';'Ray2';'RefLine'};
         data.Name='Ray1';
         data.Type='line';
@@ -1388 +1390 @@
         UvData.ProjObject{4}.DisplayHandle.uvmat=plot_object(UvData.ProjObject{4},UvData.ProjObject{1},handles.PlotAxes,'b');
     end
-    if isfield(XmlData.LIFCalib,'MaskPolygon')
+    if isfield(XmlData.LIFCalib,'MaskPolygonCoord')
+        ListObjectName{5}='MaskPolygon';
         data.Name='MaskPolygon';
         data.Coord=XmlData.LIFCalib.MaskPolygonCoord;
         UvData.ProjObject{5}=data;
         UvData.ProjObject{5}.DisplayHandle.uvmat=plot_object(UvData.ProjObject{5},UvData.ProjObject{1},handles.PlotAxes,'b');
-    end
-end
-set(handles.ListObject,'String',{'plane';'Ray1';'Ray2';'RefLine';'MaskPolygon'})
+    end      
+end
+set(handles.ListObject,'String',ListObjectName)
 set(handles.ListObject,'Value',1)
 set(handles.uvmat,'UserData',UvData);
@@ -1401 +1404 @@
 %% read lines currently drawn
 ListObj=UvData.ProjObject;% get the current list of projection objects 
-select=zeros(1,numel(ListObj));
+select_line=zeros(1,numel(ListObj));
 select_mask=zeros(1,numel(ListObj));
 for iobj=1:numel(ListObj);
     if isfield(ListObj{iobj},'Type') && strcmp(ListObj{iobj}.Type,'line')
-        select(iobj)=1;% select the lines among the projection objects
+        select_line(iobj)=1;% select the lines among the projection objects
     end
     if isfield(ListObj{iobj},'Type') && strcmp(ListObj{iobj}.Type,'polygon')
@@ -1411 +1414 @@
     end
 end
-val=find(select);
-if numel(val)<2
+if numel(find(select_line))<2
     msgbox_uvmat('ERROR',{'light rays must be defined by at least two lines created by Projection object/line in the menu bar'; ...
     'use a third line to get a reference luminosity profile accross the illumination beam'});
     return
 else
-    ObjectData=UvData.ProjObject(val);
-    for iobj=1:length(ObjectData)
-            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
-index_mask=find(select_mask);
-if numel(index_mask)>=1    
-    MaskData=UvData.ProjObject(index_mask);
-    for iobj=1:length(MaskData)
-           
-    end
-end
+    LineData=UvData.ProjObject(find(select_line));
+    for iobj=1:length(LineData)
+            xA(iobj)=LineData{iobj}.Coord(1,1);
+            yA(iobj)=LineData{iobj}.Coord(1,2);
+            xB(iobj)=LineData{iobj}.Coord(2,1);
+            yB(iobj)=LineData{iobj}.Coord(2,2);
+    end
+end
+
 
 %% set the image offset
@@ -1442 +1438 @@
     return
 else
-    XmlData.LIFCalib.BlackOffset=answer ;% image value for black background, to be determined by taking images with a cover on the objective lens
+    XmlData.LIFCalib.BlackOffset=str2num(answer) ;% image value for black background, to be determined by taking images with a cover on the objective lens
 end
 
@@ -1458 +1454 @@
 y0=((y3-y4)*(x1*y2-y1*x2)-(y1-y2)*(x3*y4-y3*x4))/D;
 XmlData.LIFCalib.LightOrigin=[x0 y0];% origin of the light source to be saved in the current xml file
-XmlData.LIFCalib.Ray1Coord=ObjectData{1}.Coord;
-XmlData.LIFCalib.Ray2Coord=ObjectData{2}.Coord;
-XmlData.LIFCalib.RefLineCoord=ObjectData{3}.Coord;
+XmlData.LIFCalib.Ray1Coord=LineData{1}.Coord;
+XmlData.LIFCalib.Ray2Coord=LineData{2}.Coord;
+XmlData.LIFCalib.RefLineCoord=LineData{3}.Coord;
 
 %% display the current image in polar axes with origin at the  illumination source
@@ -1472 +1468 @@
 
 %% use the third line for reference luminosity, renormalize the image intensity along each ray to get a uniform brightness along this line
-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);
+if numel(find(select_line))==3
+    x_ref=linspace(LineData{3}.Coord(1,1),LineData{3}.Coord(2,1),10);
+    y_ref=linspace(LineData{3}.Coord(1,2),LineData{3}.Coord(2,2),10);
     x_ref=x_ref-x0;
     y_ref=y_ref-y0;
@@ -1501 +1497 @@
     xlabel('azimuth (pixels)')
     ylabel('image brightness')
-    %ImaName=regexprep([get(handles.RootFile,'String') get(handles.FileIndex,'String')],'//','');
-NewImageName=fullfile(get(handles.RootPath,'String'),'LIF_ref.png');
-imwrite(uint16(1000*Anorm),NewImageName,'BitDepth',16)
-DataPol.A=uint16(1000*Anorm);
+    title('ref brightness profile')
+end
+
+%% get the image A*r
+[npy,npx]=size(Anorm);
+AX=DataPol.Coord_x;
+AY=DataPol.Coord_y;
+dX=(AX(2)-AX(1))/(npx-1);
+dY=(AY(1)-AY(2))/(npy-1);%mesh of new pixels
+[R,Theta]=meshgrid(linspace(AX(1),AX(2),npx),linspace(AY(1),AY(2),npy));
+A=R.*Anorm;
+A=(A<=0).*ones(npy,npx)+(A>0).*A; %replaces zeros by ones
+A=log(A); % rA(r) should have a slope -gamma in x (radius from laser source)
+
+%% get the exponential decay law
+index_mask=find(select_mask);
+if ~isempty(index_mask)
+    if numel(index_mask)>1
+        msgbox_uvmat('ERROR',{'only one mask polygon accepted'});
+        return
+    else
+        MaskData=UvData.ProjObject{index_mask};
+        XmlData.LIFCalib.MaskPolygonCoord=MaskData.Coord;
+    end
+end
+
+%loop on lines iY (angle in polar coordiantes)
+gamma_coeff=zeros(1,npy);
+for iY=1:npy
+    ALine=A(iY,:);
+    RLine=R(iY,:);
+    ThetaLine=Theta(iY,:)*pi/180;
+    [XLine,YLine] = pol2cart(ThetaLine,RLine);
+    XLine=XLine+x0;
+    YLine=YLine+y0;  
+    if ~isempty(index_mask)
+        ind_good=inpolygon(XLine,YLine,MaskData.Coord(:,1),MaskData.Coord(:,2));
+        if numel(find(ind_good))>100
+        ALine=ALine(ind_good);
+        RLine=RLine(ind_good);
+        else
+            continue
+        end
+    end
+    p = polyfit(RLine,ALine,1);
+        gamma_coeff(iY)=-p(1);
+end
+
+%% show and store the rescaled image
+NewImageName=fullfile(get(handles.RootPath,'String'),'LIF_ref_rescaled.png');
+DataPol.A=uint16(1000*A);
 view_field(DataPol);
-end
-
-% figure(12)
-% imagesc(DataPol.Coord_x,DataPol.Coord_y,AY,uint16(1000*Anorm))
-%title('rescaled image')
+imwrite(DataPol.A,NewImageName,'BitDepth',16)
+figure(13)
+plot(Theta(:,1),100*gamma_coeff)
+xlabel('angle(degree)')
+ylabel('decay coeff(m-1)')
 
 %% record the calibration data 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?');
@@ -1551 +1593 @@
          t=delete(t,uid_line);
     end
-       uid_line=find(t,'ImaDoc/LIFCalib/RefLineCoord');
-    if ~isempty(uid_line)  %if GeometryCalib does not already exists, create it
-         t=delete(t,uid_line);
+       uid_mask=find(t,'ImaDoc/LIFCalib/MaskPolygonCoord');
+    if ~isempty(uid_mask)  %if GeometryCalib does not already exists, create it
+         t=delete(t,uid_mask);
     end
      uid_BlackOffset=find(t,'ImaDoc/LIFCalib/BlackOffset');
@@ -1563 +1605 @@
     save(t,XmlFileName);
 end
-
-
-%% get the image A*r
-[npy,npx]=size(DataPol.A);
-AX=DataPol.Coord_x;
-AY=DataPol.Coord_y;
-dX=(AX(2)-AX(1))/(npx-1);
-dY=(AY(1)-AY(2))/(npy-1);%mesh of new pixels
-r=AX(1)+[0:npx-1]*dX;%distance from laser
-rmat=ones(npy,1)*r;
-A=rmat.*Anorm;
-A=(A<=0).*ones(npy,npx)+(A>0).*A; %replaces zeros by ones
-A=log(A); % rA(r) should have a slope -gamma in x (radius from laser source)
-
-%%  find mask image: for use in uvmat
-% if ~isempty(DataIn.ZIndex)
-% num_level=DataIn.ZIndex;% needed to handle multilevel illumination
-% if ~exist('maskname','var')
-%     huvmat=findobj(allchild(0),'tag','uvmat');
-%     hhuvmat=guidata(huvmat);
-%     RootPath=get(hhuvmat.RootPath,'String');
-%     RootFile=get(hhuvmat.RootFile,'String');
-%     maskname=fullfile(RootPath,[RootFile '_23mask_' num2str(num_level) '.png']);
-%     if test_1% PIV image used for correction
-%         RootPath_1=get(hhuvmat.RootPath_1,'String');
-%         RootFile_1=get(hhuvmat.RootFile_1,'String');
-%         maskname_1=fullfile(RootPath_1,[RootFile_1 '_23mask_' num2str(num_level) '.png']);
-%     end
-% end
-% Mask=DataIn;
-% if test_1
-%     Mask_1=DataIn_1;
-%     Mask.A=imread(maskname);
-%     Mask_1.A=imread(maskname_1);
-%     [Mask,Mask_1]=phys_polar(Mask,XmlData,Mask_1,XmlData_1);
-%     Mask.A=Mask.A >= 200 & Mask_1.A >= 200 ;
-% % end
-%  end
-% Mask.A=ones(size(DataPol.A));
-[R,Y]=meshgrid(linspace(AX(1),AX(2),npx),linspace(AY(1),AY(2),npy));
-[X,Y] = pol2cart(Y*pi/180,R);
-X=X+x0;% ny*nx matrix of phys cartesian coordiantes
-Y=Y+y0;
-gamma_coeff=zeros(1,npy);
-
-NewImageName=fullfile(get(handles.RootPath,'String'),'LIF_ref_rescaled.png');
-imwrite(uint16(1000*A),NewImageName,'BitDepth',16)
-return
-%loop on lines iY (angle in polar coordiantes)
-for iY=1:npy
-    Theta_i=(AY(1)-(iY-1)*dY)*pi/180;%theta in radians
-    %r_edge(iY)=(-100-XmlData.GeometryCalib.PolarCentre(1))/cos(Theta_i);%r_edge=radius of the edge (at X=-95)
-    r_edge(iY)
-    Aline=A(iY,:);% =r*image intensity
-    Xline=X(iY,:);
-    % adjust the fit region
-    test_fit=X(iY,:)>10 & X(iY,:)<70 & A(iY,:)>0 & Mask.A(iY,:);
-    ind_fit=find(test_fit);% indices of points used for the exponential fit
-    if numel(ind_fit)<50
-        Mask.A(iY,:)=0;% to exclude lines for which the line of rays is too short
-    else
-        p = polyfit(r(ind_fit),Aline(ind_fit),1);
-        Aline_fit=polyval(p,r);
-        if iY==round(npy/2)
-            figure(12)
-            plot(r(ind_fit),Aline(ind_fit),'b+',r,Aline_fit,'r')
-            test_out=(Aline_fit-Aline)>0.5
-            ind_out=find(test_out & test_fit)
-        end
-        test_out=(Aline_fit-Aline)>0.5;% eliminate points too dark compared to the fit, such that log(Afit/A)>0.5
-        ind_out=find(test_out & test_fit);
-%         if ~isempty(find(test_out & test_fit,1))
-%             Mask.A(iY,test_out)=0;
-         if ~isempty(ind_out)
-              Mask.A(iY,ind_out)=0;
-            ind_fit=find(~test_out & test_fit);
-            p = polyfit(r(ind_fit),Aline(ind_fit),1);%do new fit after elimination of the dark points
-        end
-        if -p(1) > 0.003 && -p(1) < 0.012
-            gamma_coeff(iY)=-p(1);
-            Aline_edge=polyval(p,r_edge(iY));
-            DataPol.A(iY,:)=exp(Aline_edge)/(r_edge(iY));
-        end
-
-    end
-end
-gamma_select=gamma_coeff(gamma_coeff~=0);
-gamma_select=filter(ones(1,20)/20,1,gamma_select);%smooth gamma_ceff over 20 pixels
-gamma_coeff(gamma_coeff~=0)=gamma_select;
-DataPol.A(Mask.A==0)=0;
-Ref.r_edge=r_edge;
-Ref.GammaCoeff=gamma_coeff
-
-
-
-
-
-
-