Index: /trunk/src/LIF_series.m =================================================================== --- /trunk/src/LIF_series.m (revision 566) +++ /trunk/src/LIF_series.m (revision 566) @@ -0,0 +1,136 @@ +%---------------------------------------------------------------------- +% -process LIF images +%---------------------------------------------------------------------- +function GUI_input=LIF_series(num_i1,num_i2,num_j1,num_j2,Series); + +%requests for the visibility of input windows in the GUI series (activated directly by the selection in the menu ACTION) +if ~exist('num_i1','var') + GUI_input={'RootPath';'two';...%nbre of possible input series (options 'on'/'two'/'many', default:'one') + 'SubDir';'on';... % subdirectory of derived files (PIV fields), ('on' by default) + 'RootFile';'on';... %root input file name ('on' by default) + 'FileExt';'on';... %input file extension ('on' by default) + 'NomType';'on';...%type of file indexing ('on' by default) + 'NbSlice';'on'; ...%nbre of slices ('off' by default) + 'VelTypeMenu';'one';...% menu for selecting the velocity type (civ1,..) options 'off'/'one'/'two', 'off' by default) + 'FieldMenu';'one';...% menu for selecting the field (s) in the input file(options 'off'/'one'/'two', 'off' by default) + 'CoordType';'on';...%can use a transform function 'off' by default + 'GetObject';'on';...%can use projection object ,'off' by default + %'GetMask';'on'...%can use mask option ,'off' by default + %'PARAMETER'; options: name of the user defined parameter',repeat a line for each parameter + ''}; + return %exit the function +end + +%------------------------------------------------- +hseries=guidata(Series.hseries);%handles of the GUI series +WaitbarPos=get(hseries.waitbar_frame,'Position'); %positiopn of waitbar frame +%------------------------------------------------- +addpath '/fsnet/project/coriolis/2006/06ICEBOX/0_MATLAB_WORK/LIF'% define path for concentration.m +% cpath=which('series'); + +%mode=''; %default +time=[]; +nbviews=numel(Series.RootPath); +if nbviews<2 + msgbox_uvmat('ERROR','enter both LIF and PIV series')% we introduce PIV series to improve the filtering (remove particle image) + return +end +[PD,LIFdir]=fileparts(Series.RootPath{1}); + +fulldir=fullfile(PD,'Concentration'); +if ~exist(fulldir,'dir') + try + mkdir(fulldir) + [xx,msg2] = fileattrib(fulldir,'+w','g'); %yield writing access (+w) to user group (g) + catch ME + msgbox_uvmat('ERROR',ME.message) + return + end +end + +filebase_LIF=fullfile(fulldir,'LIF');%root name for the merged files +RootPath=Series.RootPath{1}; +filebase=fullfile(Series.RootPath{1},Series.RootFile{1});%root file name +filebase_1=fullfile(Series.RootPath{2},Series.RootFile{2});%root file name for PIV (background correction) +nbfield=numel(num_i1{1});%number of fields in the series +[XmlData,error]=imadoc2struct([filebase '.xml']);% calibration data for LIF +%[error,Heading,nom_type_read,ext_ima_read,tt,TimeUnit,mode,NbSlice,npx,npy,Calib{2}] +[XmlData_1,error]=imadoc2struct([filebase_1 '.xml']);% calibration data for PIV +if isfield(XmlData,'Time') + time=XmlData.Time; + if isfield(XmlData_1,'Time') + time_1=XmlData_1.Time; + if ~isequal(size(time),size(time_1)) + msgbox_uvmat('WARNING','inconsistent time array lengths in ImaDoc fields') + end + end +end +%check coincidence in time +if size(time,1)>1 + diff_time=max(max(abs(time-time_1))) + if diff_time>0 + msgbox_uvmat(['times of series differ by more than ' num2str(diff_time)],'WARNING') + end +end + + +hRUN=findobj(Series.hseries,'Tag','RUN');%handles the the uicontrol 'RUN' +itime=0; +%%%%%% +%LOOP ON FILES + RootPath=fullfile(RootPath,'LIF_REF'); +for ifile=1:nbfield + stopstate=get(hRUN,'BusyAction');%enable stop button + if isequal(stopstate,'queue')% enable STOP command + update_waitbar(hseries.waitbar,WaitbarPos,ifile/nbfield) + %name of the current LIF input file + [inputfile,idetect]=name_generator(filebase,num_i1{1}(ifile),num_j1{1}(ifile),Series.FileExt{1},Series.NomType{1},1,num_i1{1}(ifile),num_j2{1}(ifile)); + if ~idetect + msgbox_uvmat('ERROR',[inputfile ' not found']) + return + end + [Data,ParamOut,errormsg] = read_field(inputfile,'image',[]); + Data.ZIndex=num_i1{1}(ifile)-Series.NbSlice*(floor((num_i1{1}(ifile)-1)/Series.NbSlice));%second field index + + + file_ref=fullfile(RootPath,['lif_ref_' num2str(Data.ZIndex) '.nc']); + Ref=nc2struct(file_ref);%reference file + [inputfile_1,idetect]=name_generator(filebase_1,num_i1{2}(ifile),num_j1{2}(ifile),Series.FileExt{2},Series.NomType{2},1,num_i2{2}(ifile),num_j2{2}(ifile)); + if ~idetect + msgbox_uvmat('ERROR',[inputfile_1 ' not found']) + return + end + Data_1=read_field(inputfile_1,'image',{[]});% read the image + Data_1.ZIndex=Data.ZIndex; + %%% transform image to concentration + [DataOut,dd,DataMask]=concentration(Data,XmlData,Data_1,XmlData_1,Ref); + % output file name (netcdf) + outputfile=name_generator(filebase_LIF,num_i1{1}(ifile),num_j1{1}(ifile),'.nc',Series.NomType{2},1,num_i2{1}(ifile),num_j2{1}(ifile)); + % create a structure to prepare the result file + Resu.ListGlobalAttribute={'Project','InputFile_1','InputFile_2','Action','Time','ZIndex','z'}; + [PP,Resu.Project]=fileparts(Series.PathProject); + Resu.InputFile_1=inputfile; + Resu.InputFile_2=inputfile_1; + Resu.Action=Series.Action; + if isempty(time) + Resu.Time=0; + else + Resu.Time=time(num_i1{1}(ifile),num_j1{1}(ifile)); + end + Resu.ZIndex=Data.ZIndex; + Resu.z=XmlData.GeometryCalib.SliceCoord(Data.ZIndex,3); + Resu.ListVarName={'AY' ,'AX' ,'c','mask'}; + Resu.VarDimName={'AY','AX',{'AY','AX'},{'AY','AX'}}; + Resu.AY=[DataOut.AY(1), DataOut.AY(end)]; + Resu.AX=[DataOut.AX(1), DataOut.AX(end)]; + Resu.c=DataOut.A; + Resu.mask=DataMask.A;%to chnge to cartesian coordinates (polar2phys) + error=struct2nc(outputfile,Resu); %save result file + if isempty(error) + display(['output file ' outputfile ' written']) + else + display( error) + end + end +end + Index: /trunk/src/concentration.m =================================================================== --- /trunk/src/concentration.m (revision 566) +++ /trunk/src/concentration.m (revision 566) @@ -0,0 +1,146 @@ +%transform LIF images to concentration images + +function [DataOut,DataOut_1,DataMask]=concentration(Data,XmlData,Data_1,XmlData_1,Ref) +cpath=which('uvmat'); +addpath(fullfile(fileparts(cpath),'transform_field'))% define path for phys_polar.m +DataOut_1=[]; + +%% for use in uvmat +num_level=Data.ZIndex; +if ~exist('Ref','var') + huvmat=findobj(allchild(0),'tag','uvmat'); + hhuvmat=guidata(huvmat); + RootPath=get(hhuvmat.RootPath,'String'); + + %reference file + RootPath=fullfile(RootPath,'LIF_REF'); + file_ref=fullfile(RootPath,['lif_ref_' num2str(num_level) '.nc']); + Ref=nc2struct(file_ref); +end + +%% Parameters +XmlData.GeometryCalib.PolarCentre=Ref.IlluminationOrigin;%[-515 -175]; %position of the laser origin [x, y] +XmlData_1.GeometryCalib.PolarCentre=Ref.IlluminationOrigin;%[-515 -175]; %position of the laser origin [x, y] +ImageOffset=Ref.ImageOffset; %237;% image value for black background +nfilt=64; + +%% concentration image +Data.A(Ref.CoverIndex:end,:)=Ref.CoverCoeff*(double(Data.A(Ref.CoverIndex:end,:))-ImageOffset(1))+ImageOffset(1);% COMPENSATION OF BRIGHTNESS UNDER THE COVER +[DataOut,DataOut_1]=phys_polar(Data,XmlData,Data_1,XmlData_1); +A=Ref.Aref;%default +ind_good=find(Ref.Aref~=0); +ind_bad=find(Ref.Aref==0); +A(ind_good)=double(DataOut.A(ind_good))-ImageOffset(1)-0.07*(double(DataOut_1.A(ind_good))-ImageOffset(2));%substract PIV image information for removing particles +%filtering and decimate +Afilt=filter2(ones(nfilt,nfilt),A); +Mask=filter2(ones(nfilt,nfilt),double(Ref.Aref~=0)); +B=Afilt./Mask; +A(ind_bad)=B(ind_bad); +[npy,npx]=size(A); +DataMask=DataOut; +DataMask.A=2*ones(npy,npx);%mask=2 for good data + +DataMask.A(Ref.Aref==0)=1;%mask=0 for undefined data + + + +C=filter2(ones(nfilt,nfilt),Ref.Aref); +D=C./Mask; +Ref.Aref(ind_bad)=D(ind_bad); +DataOut_1=[]; +AX=DataOut.AX; +AY=DataOut.AY; + +dX=(AX(2)-AX(1))/(npx-1); +dY=(AY(1)-AY(2))/(npy-1);%mesh of new pixels +[R,Y]=meshgrid(linspace(AX(1),AX(2),npx),linspace(AY(1),AY(2),npy)); +r=AX(1)+[0:npx-1]*dX;%distance from laser +%A(ind_good)=(A(ind_good)>=0).*A(ind_good); %replaces negative values by zeros +A=A./Ref.Aref;% luminosity normalised by the reference (value at the edge of the box) + +%% Interpolation +% [Rindex,Yindex]=meshgrid(linspace(0.5,npx-0.5,npx),linspace(npy-0.5,0.5,npy)); +% Rgood=Rindex(ind_good); +% Ygood=Yindex(ind_good); +%F=TriScatteredInterp(Rgood,Ygood,A(ind_good)); +%A=F(Rindex,Yindex); + + +DataMask.A(isnan(A)|isinf(A)|A>1.5)=0;% mask=1 for interpolated data +r_edge=Ref.r_edge*ones(1,npx); +Edge_ind=find((abs(R-r_edge)/dX)<=1 & DataMask.A~=0);%indies of positions close to r_edge, values greater than 1 are not expected +yedge=min(min(Y(Edge_ind))); +jmax=round(-(yedge-AY(1))/dY+1); +DataMask.A(jmax:end,:)=0; + +A(isnan(A)|isinf(A))=0; + +% radius along the reference line +Theta=(linspace(AY(1),AY(2),npy)*pi/180)'*ones(1,npx);%theta in radians + +gamma_coeff=Ref.GammaCoeff*ones(1,npx); + +A(R=1 & R<=npx & Theta >=1 & Theta<=npy;%flagin=1 inside the original image +vec_A=reshape(DataIn.A,1,npx*npy);%put the original image in line +ind_in=find(flagin); +ind_out=find(~flagin); +ICOMB=((R-1)*npy+(npy+1-Theta)); +ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A +vec_B(ind_in)=vec_A(ICOMB); +vec_B(ind_out)=zeros(size(ind_out)); +DataOut.A=flipdim(reshape(vec_B,npy,npx),1);%new image in real coordinates + + %Rangx=Rangx-radius_ref; +DataOut.AX=RangeX; +DataOut.AY=RangeY; + Index: /trunk/src/lumin_calib.m =================================================================== --- /trunk/src/lumin_calib.m (revision 566) +++ /trunk/src/lumin_calib.m (revision 566) @@ -0,0 +1,77 @@ +% 'lumin_cali': check the luminosity of the camera lens versus distance to image center +%---------------------------------------------------------------------- + +function input_list=lumin_calib (num_i1,num_i2,num_j1,num_j2,Series) + +%display request of input parameters in the series interface (activated directly by the selection in the menu ACTION) +input_list={'NbSlice';...%nbre of slices + 'RootPath';... %path to the input root file + 'RootFile';... %root input file name + 'FileExt';... %input file extension + 'NomType';...%type of file indexing + 'SubDir';...% subdirectory of derived files (PIV fields) + }; + +if ~exist('num_i1','var') + return % used to display the right input data from input_list when the function is selected +end +for ilist=1:length(input_list) + eval([input_list{ilist} '= Series.' input_list{ilist} ';']) +end + +%----------------------------------------------------------------- + +hRUN=findobj(Series.hseries,'Tag','RUN'); +hwaitbar=findobj(Series.hseries,'Tag','waitbar'); +waitbarpos(1)=Series.WaitbarPos(1);%x position of the waitbar +waitbarpos(3)=Series.WaitbarPos(3);% width of the waitbar +siz=size(num_i1); + + +xpos(1:10)=[297 297 297 314 316 310 310 310 312 325]; +ypos(1:10)=[818 818 818 817 822 820 820 820 813 799]; +xpos(11:30)=[364 403 448 479 522 586 613 630 626 581 514 470 380 328 292 292 292 292 292 299]; +ypos(11:30)=[733 673 604 558 504 405 357 339 337 409 509 570 701 789 845 845 845 845 845 836]; +xpos(31:49)=[353 273 269 269 269 269 269 270 270 270 270 270 270 270 270 268 268 264 216 ]; +ypos(31:49)=[748 875 883 883 883 883 883 888 888 888 888 888 888 888 888 888 888 885 959]; +radius=50; +Abackground=58.6; + +filebase=fullfile(Series.RootPath{1},Series.RootFile{1}); +dir_images=Series.RootPath{1}; +nom_type=Series.NomType{1}; +% filebad=zeros(size(num_i1)); +indbad=[21 22 30];%bad image +file_index= 1:length(num_i1); +ifile=0; +for ifile= 1:length(num_i1) + [filename,idetect]=... + name_generator(filebase,num_i1(ifile),num_j1(ifile),Series.FileExt{1},Series.NomType{1}); + + A=imread(filename); + A=flipdim(A,1); + ind_x=[xpos(ifile)-radius:xpos(ifile)+radius]; + ind_y=[ypos(ifile)-radius:ypos(ifile)+radius]; + A=double(A(ind_y,ind_x))-Abackground; + [Xi,Yi]=meshgrid([-radius:+radius],[-radius:+radius]); +% distX=(Xi+radius+1); +% distY=(Yi+radius+1); + testin=(Xi.*Xi+Yi.*Yi)<=radius*radius; + testin=double(testin); + A=(A>0).*A.*double(testin); + Amean(ifile)=sum(sum(A))/sum(sum(testin)); + +end +dist=((xpos-512).*(xpos-512))+((ypos-512).*(ypos-512)); +dist=dist/(512*512); +file_index(indbad)=[]; +Amean(indbad)=[]; +dist(indbad)=[]; +figure(1) +plot(file_index,Amean,'+') +figure(2) +haxes=axes; +h=plot(dist,Amean,'+') +xlabel('distance to image center (pixels)') +ylabel('luminosity of reference input dye solution') +set(haxes,'YLim',[0 2850]) Index: /trunk/src/uvmat.m =================================================================== --- /trunk/src/uvmat.m (revision 565) +++ /trunk/src/uvmat.m (revision 566) @@ -4594,4 +4594,55 @@ set(handles.MenuCalib,'checked','on')% indicate that MenuCalib is activated, test used by mouse action + +%----------------------------------------------------------------------- +function MenuLIFCalib_Callback(hObject, eventdata, handles) +%------------------------------------------------------------------------ +UvData=get(handles.uvmat,'UserData');%read UvData properties stored on the uvmat interface +ListObj=UvData.Object; +select=zeros(1,numel(ListObj)); +for iobj=1:numel(ListObj); + if strcmp(ListObj{iobj}.Type,'line') + select(iobj)=1; + end +end +val=find(select); +if numel(val)<2 + msgbox_uvmat('ERROR','light rays must be defined by at least two lines created by Projection object/line in the menu bar'); + return +else + set(handles.ListObject,'Value',val); + ObjectData=UvData.Object(val); + flag=1; + npx=size(UvData.Field.A,2); + npy=size(UvData.Field.A,1); + xi=0.5:npx-0.5; + yi=0.5:npy-0.5; + [Xi,Yi]=meshgrid(xi,yi); + for iobj=1:length(ObjectData) + flagobj=1; + testphys=0; %coordinates in pixels by default + if isfield(ObjectData,'CoordUnit') && ~isequal(ObjectData.CoordUnit,'pixel') + if isfield(UvData,'XmlData')&& isfield(UvData.XmlData{1},'GeometryCalib') + Calib=UvData.XmlData{1}.GeometryCalib; + testphys=1; + end + end + if isfield(ObjectData{iobj},'Coord') + x1(iobj)=ObjectData{iobj}.Coord(1,1); + y1(iobj)=ObjectData{iobj}.Coord(1,2); + x2(iobj)=ObjectData{iobj}.Coord(2,1); + y2(iobj)=ObjectData{iobj}.Coord(2,2); + end + end +end + %determine the ray origin + x1 + y1 + x2 + y2 + % update the xml file + + + %------------------------------------------------------------------------ function MenuMask_Callback(hObject, eventdata, handles) @@ -4617,5 +4668,5 @@ yi=0.5:npy-0.5; [Xi,Yi]=meshgrid(xi,yi); - if isfield(UvData,'Object') +% if isfield(UvData,'Object') for iobj=1:length(UvData.Object) ObjectData=UvData.Object{iobj}; @@ -4676,5 +4727,5 @@ end end - end +% end %mask name RootPath=get(handles.RootPath,'String'); @@ -4846,2 +4897,3 @@ % Hint: get(hObject,'Value') returns toggle state of CheckColorBar +