Index: /trunk/src/find_field_cells.m =================================================================== --- /trunk/src/find_field_cells.m (revision 1048) +++ /trunk/src/find_field_cells.m (revision 1049) @@ -475,5 +475,5 @@ 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; Index: /trunk/src/series/extract_multitif_special.m =================================================================== --- /trunk/src/series/extract_multitif_special.m (revision 1048) +++ /trunk/src/series/extract_multitif_special.m (revision 1049) @@ -144,6 +144,12 @@ % 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,10 +160,11 @@ % 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,12 +205,12 @@ 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; Index: /trunk/src/uvmat.m =================================================================== --- /trunk/src/uvmat.m (revision 1048) +++ /trunk/src/uvmat.m (revision 1049) @@ -1371,6 +1371,8 @@ %% 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,12 +1390,13 @@ 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,9 +1404,9 @@ %% 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,25 +1414,18 @@ 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,5 +1438,5 @@ 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,7 +1454,7 @@ 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,7 +1468,7 @@ %% 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,17 +1497,63 @@ 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,7 +1593,7 @@ 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,103 +1605,4 @@ 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 - - - - - - -