Index: /trunk/src/series/relabel_i_j.m =================================================================== --- /trunk/src/series/relabel_i_j.m (revision 213) +++ /trunk/src/series/relabel_i_j.m (revision 214) @@ -1,14 +1,12 @@ -% relabel_i_j: relabel an image series with two indices, according to the time matrix given by ImaDoc -% (specific to RDvision system) +% relabel_i_j: relabel an image series with two indices, and correct errors from the RDvision transfer program %---------------------------------------------------------------------- function GUI_input=relabel_i_j(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={}; - return %exit the function -end + +GUI_input={}; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%enable waitbar -hseries=guidata(Series.hseries);%handles of the GUI series +hGUI=findobj(allchild(0),'name','series'); +hseries=guidata(hGUI);%handles of the GUI series WaitbarPos=get(hseries.waitbar_frame,'Position'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -17,35 +15,49 @@ display('RDvision system') first_label=0 %image numbers start from 0 -errorfactor=2 %correct a factor of 2 in NbDk+1 - -%% create dir of the new images -basename=fullfile(Series.RootPath,Series.RootFile) ; -[dir_images,namebase]=fileparts(basename); -[path,subdir_ima]=fileparts(dir_images); -newdir=fullfile(path,[subdir_ima '_ij']); -mkdir(newdir); -[xx,msg2] = fileattrib(newdir,'+w','g'); %yield writing access (+w) to user group (g) -if ~strcmp(msg2,'') - msgbox_uvmat('ERROR',['pb of permission for ' subdir_ima ': ' msg2])%error message for directory creation +%errorfactor=1 %correct a factor of 2 in NbDk+1 + +%% read imadoc +RootPath=get(hseries.RootPath,'String'); +RootFile=get(hseries.RootFile,'String'); +basename=fullfile(RootPath{1},RootFile{1}); +SeriesData=get(hGUI,'UserData'); +if ~strcmp(SeriesData.NomType,'_000001') + msgbox_uvmat('WARNING','the input is not a file from RDvision: this function relabel_i_j has no action')%error message for directory creation return -end -display(['relabelled images in the directory ' newdir]) -basename_new=fullfile(newdir,namebase); - -%% read imadoc -[XmlData,warntext]=imadoc2struct([basename '.xml']); -nbfield1=size(XmlData.Time,1)/errorfactor; +else + msgbox_uvmat('CONFIRMATION','this function will relabel the file series from RDvision and correct the xml file')%error message for directory creation +end +[XmlData,warntext]=imadoc2struct([basename '.xml'])% read the xml file appended to the present function (containing bug corrections) +if ~isempty(warntext) + msgbox_uvmat('ERROR',warntext)%error message for xml file reading +end +nbfield1=size(XmlData.Time,1); nbfield2=size(XmlData.Time,2); + set(hseries.first_i,'String',num2str(first_label))% display the first image in the process set(hseries.last_i,'String',num2str(nbfield1*nbfield2-1+first_label))% display the last image in the process +set(hseries.nb_field,'String',{num2str(nbfield1*nbfield2-1+first_label)})% display the total nbre of images + +%% stop program ther when it is selected in the menu (no run action) +if ~exist('num_i1','var') + return +end + +answer=msgbox_uvmat('CONFIRMATION',[num2str(nbfield1) ' bursts containing ' num2str(nbfield2) ' images each'])%error message for directory creation %% apply the image rescaling function 'level' (avoid bright particles) -answer=msgbox_uvmat('INPUT_Y-N','apply image rescaling function levels.m'); -test_level=isequal(answer,'Yes'); +% answer=msgbox_uvmat('INPUT_Y-N','apply image rescaling function levels.m'); +% test_level=isequal(answer,'Yes'); %% copy and adapt the xml file if exist([basename '.xml'],'file') - copyfile([basename '.xml'],[basename_new '.xml']);% copy the .civ file - t=xmltree([basename_new '.xml']); + try + copyfile([basename '.xml'],[basename '.xml~']);% backup the xml file + catch + errormsg=lasterr + msgbox_uvmat('ERROR',errormsg); + return + end + t=xmltree([basename '.xml']); %update information on the first image name in the series @@ -55,5 +67,5 @@ end uid_ImageName=find(t,'ImaDoc/Heading/ImageName'); - ImageName=name_generator(basename_new,1,1,'.png','_i_j'); + ImageName=name_generator(basename,1,1,'.png','_i_j'); [pth,ImageName]=fileparts(ImageName); ImageName=[ImageName '.png'] @@ -63,20 +75,17 @@ uid_value=children(t,uid_ImageName); if isempty(uid_value) - t=add(t,uid_ImageName,'chardata',ImageName)%indicate name of the first image, with ;png extension + t=add(t,uid_ImageName,'chardata',ImageName);%indicate name of the first image, with ;png extension else - t=set(t,uid_value(1),'value',ImageName)%indicate name of the first image, with ;png extension - end - -% %add information about image transform -% [t,new_uid]=add(t,1,'element','ImageTransform'); -% [t,NameFunction_uid]=add(t,new_uid,'element','NameFunction'); -% [t]=add(t,NameFunction_uid,'chardata','sub_background'); -% [t,NbSlice_uid]=add(t,new_uid,'element','NbSlice'); -% [t]=add(t,new_uid,'chardata',num2str(nbslice_i)); -% [t,NbSlidingImages_uid]=add(t,new_uid,'element','NbSlidingImages'); -% [t]=add(t,NbSlidingImages_uid,'chardata',num2str(nbaver)); -% [t,LuminosityRank_uid]=add(t,new_uid,'element','RankBackground'); -% [t]=add(t,LuminosityRank_uid,'chardata',num2str(rank));% luminosity rank almong the nbaver sliding images - save(t,[basename_new '.xml']) + t=set(t,uid_value(1),'value',ImageName);%indicate name of the first image, with ;png extension + end + + %%%% correction RDvision %%%% + uid_NbDtj=find(t,'ImaDoc/Camera/BurstTiming/NbDtj'); + t=set(t,uid_NbDtj,'value',num2str(s.NbDtj)) + uid_NbDtk=find(t,'ImaDoc/Camera/BurstTiming/NbDtk'); + t=set(t,uid_NbDtk,'value',num2str(s.NbDtk)) + %%% + + save(t,[basename '.xml']) end @@ -87,14 +96,31 @@ num_j=mod(ifile-1+first_label,nbfield2)+1; num_i=floor((ifile-1+first_label)/nbfield2)+1; - filename_new=name_generator(basename_new,num_i,num_j,'.png','_i_j'); - if test_level - A=imread(filename); - C=levels(A); - imwrite(C,filename_new) - else - copyfile(filename,filename_new); - end -end - + filename_new=name_generator(basename,num_i,num_j,'.png','_i_j'); + try + movefile(filename,filename_new); + [s,errormsg] = fileattrib(filename_new,'-w','a'); %set images to read only '-w' for all users ('a') + if ~s + msgbox_uvmat('ERROR',errormsg); + return + end + catch + errormsg=lasterr + msgbox_uvmat('ERROR',errormsg); + return + end + % if test_level + % A=imread(filename); + % C=levels(A); + % imwrite(C,filename_new) + % else + % try + % copyfile(filename,filename_new); + % catch + % errormsg=lasterr + % msgbox_uvmat('ERROR',errormsg); + % return + % end + % end +end @@ -137,2 +163,295 @@ C=(C-cmin)/(cmax-cmin)*256; C=uint8(C); + +%'imadoc2struct': reads the xml file for image documentation +%------------------------------------------------------------------------ +% function [s,errormsg]=imadoc2struct(ImaDoc,option) +% +% OUTPUT: +% s: structure representing ImaDoc +% s.Heading: information about the data hierarchical structure +% s.Time: matrix of times +% s.TimeUnit +% s.GeometryCalib: substructure containing the parameters for geometric calibration +% errormsg: error message +% +% INPUT: +% ImaDoc: full name of the xml input file with head key ImaDoc +% option: ='GeometryCalib': read the data of GeometryCalib, including source point coordinates + +function [s,errormsg]=imadoc2struct(ImaDoc,option) +%% default input and output +if ~exist('option','var') + option='*'; +end +errormsg=[];%default +s.Heading=[];%default +s.Time=[]; %default +s.TimeUnit=[]; %default +s.GeometryCalib=[]; +tsai=[];%default + +%% opening the xml file +if exist(ImaDoc,'file')~=2, errormsg=[ ImaDoc ' does not exist']; return;end;%input file does not exist +try + t=xmltree(ImaDoc); +catch + errormsg={[ImaDoc ' is not a valid xml file']; lasterr}; + display(errormsg); + return +end +uid_root=find(t,'/ImaDoc'); +if isempty(uid_root), errormsg=[ImaDoc ' is not an image documentation file ImaDoc']; return; end;%not an ImaDoc .xml file + + +%% Heading +uid_Heading=find(t,'/ImaDoc/Heading'); +if ~isempty(uid_Heading), + uid_Campaign=find(t,'/ImaDoc/Heading/Campaign'); + uid_Exp=find(t,'/ImaDoc/Heading/Experiment'); + uid_Device=find(t,'/ImaDoc/Heading/Device'); + uid_Record=find(t,'/ImaDoc/Heading/Record'); + uid_FirstImage=find(t,'/ImaDoc/Heading/ImageName'); + s.Heading.Campaign=get(t,children(t,uid_Campaign),'value'); + s.Heading.Experiment=get(t,children(t,uid_Exp),'value'); + s.Heading.Device=get(t,children(t,uid_Device),'value'); + if ~isempty(uid_Record) + s.Heading.Record=get(t,children(t,uid_Record),'value'); + end + s.Heading.ImageName=get(t,children(t,uid_FirstImage),'value'); +end + +%% Camera and timing +if strcmp(option,'*') || strcmp(option,'Camera') + uid_Camera=find(t,'/ImaDoc/Camera'); + if ~isempty(uid_Camera) + uid_ImageSize=find(t,'/ImaDoc/Camera/ImageSize'); + if ~isempty(uid_ImageSize); + ImageSize=get(t,children(t,uid_ImageSize),'value'); + xindex=findstr(ImageSize,'x'); + if length(xindex)>=2 + s.Npx=str2double(ImageSize(1:xindex(1)-1)); + s.Npy=str2double(ImageSize(xindex(1)+1:xindex(2)-1)); + end + end + uid_TimeUnit=find(t,'/ImaDoc/Camera/TimeUnit'); + if ~isempty(uid_TimeUnit) + s.TimeUnit=get(t,children(t,uid_TimeUnit),'value'); + end + uid_BurstTiming=find(t,'/ImaDoc/Camera/BurstTiming'); + if ~isempty(uid_BurstTiming) + for k=1:length(uid_BurstTiming) + subt=branch(t,uid_BurstTiming(k));%subtree under BurstTiming + % reading Dtk + Frequency=get_value(subt,'/BurstTiming/FrameFrequency',1); + Dtj=get_value(subt,'/BurstTiming/Dtj',[]); + Dtj=Dtj/Frequency;%Dtj converted from frame unit to TimeUnit (e.g. 's') + NbDtj=get_value(subt,'/BurstTiming/NbDtj',1); + %%%% correction RDvision %%%% + NbDtj=NbDtj/numel(Dtj); + s.NbDtj=NbDtj; + %%%% + Dti=get_value(subt,'/BurstTiming/Dti',[]); + Dti=Dti/Frequency;%Dtj converted from frame unit to TimeUnit (e.g. 's') + NbDti=get_value(subt,'/BurstTiming/NbDti',1); + Time_val=get_value(subt,'/BurstTiming/Time',0);%time in TimeUnit + if ~isempty(Dti) + Dti=reshape(Dti'*ones(1,NbDti),NbDti*numel(Dti),1); %concatene Dti vector NbDti times + Time_val=[Time_val;Time_val(end)+cumsum(Dti)];%append the times defined by the intervals Dti + end + if ~isempty(Dtj) + Dtj=reshape(Dtj'*ones(1,NbDtj),1,NbDtj*numel(Dtj)); %concatene Dtj vector NbDtj times + Dtj=[0 Dtj]; + Time_val=Time_val*ones(1,numel(Dtj))+ones(numel(Time_val),1)*cumsum(Dtj);% produce a time matrix with Dtj + end + % reading Dtk + Dtk=get_value(subt,'/BurstTiming/Dtk',[]); + NbDtk=get_value(subt,'/BurstTiming/NbDtk',1); + %%%% correction RDvision %%%% + NbDtk=-1+(NbDtk+1)/(NbDti+1) + s.NbDtk=NbDtk; + %%%%% + if isempty(Dtk) + s.Time=[s.Time;Time_val]; + else + for kblock=1:NbDtk+1 + Time_val_k=Time_val+(kblock-1)*Dtk; + s.Time=[s.Time;Time_val_k]; + end + end + end + end + end +end + +%% motor +if strcmp(option,'*') || strcmp(option,'GeometryCalib') + uid_subtree=find(t,'/ImaDoc/TranslationMotor'); + if length(uid_subtree)==1 + subt=branch(t,uid_subtree);%subtree under GeometryCalib + [s.TranslationMotor,errormsg]=read_subtree(subt,{'Nbslice','ZStart','ZEnd'},[1 1 1],[1 1 1]); + end +end +%% geometric calibration +if strcmp(option,'*') || strcmp(option,'GeometryCalib') + uid_GeometryCalib=find(t,'/ImaDoc/GeometryCalib'); + if ~isempty(uid_GeometryCalib) + if length(uid_GeometryCalib)>1 + errormsg=['More than one GeometryCalib in ' filecivxml]; + return + end + subt=branch(t,uid_GeometryCalib);%subtree under GeometryCalib + cont=get(subt,1,'contents'); + if ~isempty(cont) + uid_CalibrationType=find(subt,'/GeometryCalib/CalibrationType'); + if isequal(length(uid_CalibrationType),1) + tsai.CalibrationType=get(subt,children(subt,uid_CalibrationType),'value'); + end + uid_CoordUnit=find(subt,'/GeometryCalib/CoordUnit'); + if isequal(length(uid_CoordUnit),1) + tsai.CoordUnit=get(subt,children(subt,uid_CoordUnit),'value'); + end + uid_fx_fy=find(subt,'/GeometryCalib/fx_fy'); + focal=[];%default fro old convention (Reg Wilson) + if isequal(length(uid_fx_fy),1) + tsai.fx_fy=str2num(get(subt,children(subt,uid_fx_fy),'value')); + else %old convention (Reg Wilson) + uid_focal=find(subt,'/GeometryCalib/focal'); + uid_dpx_dpy=find(subt,'/GeometryCalib/dpx_dpy'); + uid_sx=find(subt,'/GeometryCalib/sx'); + if ~isempty(uid_focal) && ~isempty(uid_dpx_dpy) && ~isempty(uid_sx) + dpx_dpy=str2num(get(subt,children(subt,uid_dpx_dpy),'value')); + sx=str2num(get(subt,children(subt,uid_sx),'value')); + focal=str2num(get(subt,children(subt,uid_focal),'value')); + tsai.fx_fy(1)=sx*focal/dpx_dpy(1); + tsai.fx_fy(2)=focal/dpx_dpy(2); + end + end + uid_Cx_Cy=find(subt,'/GeometryCalib/Cx_Cy'); + if ~isempty(uid_Cx_Cy) + tsai.Cx_Cy=str2num(get(subt,children(subt,uid_Cx_Cy),'value')); + end + uid_kc=find(subt,'/GeometryCalib/kc'); + if ~isempty(uid_kc) + tsai.kc=str2double(get(subt,children(subt,uid_kc),'value')); + else %old convention (Reg Wilson) + uid_kappa1=find(subt,'/GeometryCalib/kappa1'); + if ~isempty(uid_kappa1)&& ~isempty(focal) + kappa1=str2double(get(subt,children(subt,uid_kappa1),'value')); + tsai.kc=-kappa1*focal*focal; + end + end + uid_Tx_Ty_Tz=find(subt,'/GeometryCalib/Tx_Ty_Tz'); + if ~isempty(uid_Tx_Ty_Tz) + tsai.Tx_Ty_Tz=str2num(get(subt,children(subt,uid_Tx_Ty_Tz),'value')); + end + uid_R=find(subt,'/GeometryCalib/R'); + if ~isempty(uid_R) + RR=get(subt,children(subt,uid_R),'value'); + if length(RR)==3 + tsai.R=[str2num(RR{1});str2num(RR{2});str2num(RR{3})]; + end + end + + %look for laser plane definitions + uid_Angle=find(subt,'/GeometryCalib/PlaneAngle'); + uid_Pos=find(subt,'/GeometryCalib/SliceCoord'); + if isempty(uid_Pos) + uid_Pos=find(subt,'/GeometryCalib/PlanePos');%old convention + end + if ~isempty(uid_Angle) + tsai.PlaneAngle=str2num(get(subt,children(subt,uid_Angle),'value')); + end + if ~isempty(uid_Pos) + for j=1:length(uid_Pos) + tsai.SliceCoord(j,:)=str2num(get(subt,children(subt,uid_Pos(j)),'value')); + end + uid_DZ=find(subt,'/GeometryCalib/SliceDZ'); + uid_NbSlice=find(subt,'/GeometryCalib/NbSlice'); + if ~isempty(uid_DZ) && ~isempty(uid_NbSlice) + DZ=str2double(get(subt,children(subt,uid_DZ),'value')); + NbSlice=get(subt,children(subt,uid_NbSlice),'value'); + if isequal(NbSlice,'volume') + tsai.NbSlice='volume'; + NbSlice=NbDtj+1; + else + tsai.NbSlice=str2double(NbSlice); + end + tsai.SliceCoord=ones(NbSlice,1)*tsai.SliceCoord+DZ*(0:NbSlice-1)'*[0 0 1]; + end + end + tsai.SliceAngle=get_value(subt,'/GeometryCalib/SliceAngle',[0 0 0]); + tsai.VolumeScan=get_value(subt,'/GeometryCalib/VolumeScan','n'); + tsai.InterfaceCoord=get_value(subt,'/GeometryCalib/InterfaceCoord',[0 0 0]); + tsai.RefractionIndex=get_value(subt,'/GeometryCalib/RefractionIndex',1); + + if strcmp(option,'GeometryCalib') + tsai.PointCoord=get_value(subt,'/GeometryCalib/SourceCalib/PointCoord',[0 0 0 0 0]); + end + s.GeometryCalib=tsai; + end + end +end + +%-------------------------------------------------- +% read a subtree +% INPUT: +% t: xltree +% head_element: head elelemnt of the subtree +% Data, structure containing +% .Key: element name +% .Type: type of element ('charg', 'float'....) +% .NbOccur: nbre of occurrence, NaN for un specified number +function [s,errormsg]=read_subtree(subt,Data,NbOccur,NumTest) +%-------------------------------------------------- +s=[];%default +errormsg=''; +head_element=get(subt,1,'name'); + cont=get(subt,1,'contents'); + if ~isempty(cont) + for ilist=1:length(Data) + uid_key=find(subt,[head_element '/' Data{ilist}]); + if ~isequal(length(uid_key),NbOccur(ilist)) + errormsg=['wrong number of occurence for ' Data{ilist}]; + return + end + for ival=1:length(uid_key) + val=get(subt,children(subt,uid_key(ival)),'value'); + if ~NumTest(ilist) + eval(['s.' Data{ilist} '=val;']); + else + eval(['s.' Data{ilist} '=str2double(val);']) + end + end + end + end + + +%-------------------------------------------------- +% read an xml element +function val=get_value(t,label,default) +%-------------------------------------------------- +val=default; +uid=find(t,label);%find the element iud(s) +if ~isempty(uid) %if the element named label exists + uid_child=children(t,uid);%find the children + if ~isempty(uid_child) + data=get(t,uid_child,'type');%get the type of child + if iscell(data)% case of multiple element + for icell=1:numel(data) + val_read=str2num(get(t,uid_child(icell),'value')); + if ~isempty(val_read) + val(icell,:)=val_read; + end + end +% val=val'; + else % case of unique element value + val_read=str2num(get(t,uid_child,'value')); + if ~isempty(val_read) + val=val_read; + else + val=get(t,uid_child,'value');%char string data + end + end + end +end Index: /trunk/src/series/sub_background.m =================================================================== --- /trunk/src/series/sub_background.m (revision 213) +++ /trunk/src/series/sub_background.m (revision 214) @@ -50,7 +50,7 @@ %'GetObject';...;%can use projection object %'GetMask';...;%can use mask option - 'PARAMETER';'NbSliding';... - 'PARAMETER';'VolumeScan';... - 'PARAMETER';'RankBrightness';... + %'PARAMETER';'NbSliding';... + %'PARAMETER';'VolumeScan';... + %'PARAMETER';'RankBrightness';... ''}; return %exit the function @@ -58,5 +58,5 @@ %---------------------------------------------------------------- -% initiate the waitbar +%% initiate the waitbar hseries=guidata(Series.hseries);%handles of the GUI series WaitbarPos=get(hseries.waitbar_frame,'Position'); @@ -67,5 +67,5 @@ end -%determine input image type +%% determine input image type FileType=[];%default MovieObject=[]; @@ -101,5 +101,10 @@ nbaver_init=23;%approximate number of images used for the sliding background: to be adjusted later to include an integer number of bursts -%adjust the proposed number of images in the sliding average to include an integer number of bursts + +%% apply the image rescaling function 'level' (avoid the blinking effects of bright particles) +answer=msgbox_uvmat('INPUT_Y-N','apply image rescaling function levels.m after sub_background'); +test_level=isequal(answer,'Yes'); + +%% adjust the proposed number of images in the sliding average to include an integer number of bursts if siz(2)~=1 nbaver=floor(nbaver_init/siz(1)); % number of bursts used for the sliding background, @@ -114,22 +119,30 @@ nom_type=Series.NomType; -%create dir of the new images -[dir_images,namebase]=fileparts(filebase); -[path,subdir_ima]=fileparts(dir_images); -curdir=pwd; -cd(path); -mkdir([subdir_ima '_b']); -[xx,msg2] = fileattrib([subdir_ima '_b'],'+w','g'); %yield writing access (+w) to user group (g) +%% create dir of the new images +% [dir_images,namebase]=fileparts(filebase); +if test_level + term='_b_levels'; +else + term='_b'; +end +[pp,subdir_ima]=fileparts(Series.RootPath); +try + mkdir([dir_images term]); +catch + errormsg=lasterr + msgbox_uvmat('ERROR',errormsg); + return +end +[xx,msg2] = fileattrib([dir_images term],'+w','g'); %yield writing access (+w) to user group (g) if ~strcmp(msg2,'') - msgbox_uvmat('ERROR',['pb of permission for ' subdir_ima '_b' ': ' msg2])%error message for directory creation - cd(curdir) + msgbox_uvmat('ERROR',['pb of permission for ' subdir_ima term ': ' msg2])%error message for directory creation return end -cd(curdir); -filebase_b=fullfile(path,[subdir_ima '_b'],namebase); - +filebase_b=fullfile([dir_images term],Series.RootFile); + +%% set processing parameters prompt = {'Number of images for the sliding background';'The number of positions (laser slices)';'volume scan mode (Yes/No)';... 'the luminosity rank chosen to define the background (0.1=for dense particle seeding, 0.5 (median) for sparse particles'}; -dlg_title = ['get (slice by slice) a sliding background and substract to each image, result in subdir ' subdir_ima '_b']; +dlg_title = ['get (slice by slice) a sliding background and substract to each image, result in subdir ' subdir_ima term]; num_lines= 3; def = { num2str(nbaver_init);num2str(nbslice_i);'No';'0.1'}; @@ -139,13 +152,9 @@ nbaver_ima=str2num(answer{1});%number of images for the sliding background -nbaver=ceil(nbaver_ima/siz(1))%number of bursts for the sliding background +nbaver=ceil(nbaver_ima/siz(1));%number of bursts for the sliding background if isequal(floor(nbaver/2),nbaver) nbaver=nbaver+1%put the number of burst to an odd number (so the middle burst is defined) end -% if isequal(nbaver,round(nbaver)) step=siz(1);%case of bursts: the sliding background is shifted by one burst -% else -% step=1; -% end vol_test=answer{3}; if isequal(vol_test,'Yes') @@ -164,16 +173,4 @@ end end -% nbaver=floor(nbaver_init/nbfield2); % number of bursts used for the sliding background, -% if isequal(floor(nbaver/2),nbaver) -% nbaver=nbaver+1;%put the number of burst to an odd number -% end -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -%nbaver_ima=nbaver*nbfield2;% adjust the number of sliding images A REMETRE -% if ~isequal(nbaver_ima,nbaver_init) -% hwarn=warndlg(['number of images in the sliding average adjusted to ' num2str(nbaver_ima)]); -% set(hwarn,'Units','normalized') -% set(hwarn,'Position',[0.3 0.3 0.4 0.1]) -% end -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% rank=floor(str2num(answer{4})*nbaver_ima); if rank==0 @@ -183,5 +180,4 @@ nbfield=floor(lengthtot/(nbfield2*nbslice_i));%total number of i indexes (adjusted to an integer number of slices) nbfield_slice=nbfield*nbfield2;% number of fields per slice -% test_plot=isequal(answer{5},'Yes'); %=1 to display background images if nbaver_ima > nbfield*nbfield2 errordlg('number of images in a slice smaller than the proposed number of images for the sliding average') @@ -190,10 +186,10 @@ nbfirst=(ceil(nbaver/2))*step; if nbfirst>nbaver_ima - nbfirst=ceil(nbaver_ima/2) + nbfirst=ceil(nbaver_ima/2); step=1; nbaver=nbaver_ima; end -%copy the xml file +%% copy the xml file if exist([filebase '.xml'],'file') copyfile([filebase '.xml'],[filebase_b '.xml']);% copy the .civ file @@ -223,4 +219,8 @@ [t,NameFunction_uid]=add(t,new_uid,'element','NameFunction'); [t]=add(t,NameFunction_uid,'chardata','sub_background'); + if test_levels + [t,NameFunction_uid]=add(t,new_uid,'element','NameFunction'); + [t]=add(t,NameFunction_uid,'chardata','levels'); + end [t,NbSlice_uid]=add(t,new_uid,'element','NbSlice'); [t]=add(t,new_uid,'chardata',num2str(nbslice_i)); @@ -240,5 +240,5 @@ %MAIN LOOP ON SLICES for islice=1:nbslice_i - %select the series of image indices at the level islice + %% select the series of image indices at the level islice for ifield=1:nbfield for iburst=1:nbfield2 @@ -246,5 +246,6 @@ end end - %read the first series of nbaver_ima images and sort by luminosity at each pixel + + %% read the first series of nbaver_ima images and sort by luminosity at each pixel for ifield = 1:nbaver_ima ifile=indselect(ifield); @@ -255,12 +256,5 @@ Asort=sort(Ak,3);%sort the luminosity of images at each point B=Asort(:,:,rank);%background image -% -% namemean=name_generator([filebase '_back'],islice,[],'.png','_i');% makes the file name -% imwrite(B,namemean,'BitDepth',16); % save the first background image -% ['background image for slice ' num2str(islice) ' saved in ' namemean] - %substract the background from each of the first images and save the new images -% for ifield=1:floor(nbaver_ima/2)+1 - 'first background image will be substracted' - + display( 'first background image will be substracted') for ifield=1:nbfirst Acor=double(Ak(:,:,ifield))-double(B);%substract background to the current image @@ -269,10 +263,15 @@ ifile=indselect(ifield); newname=name_generator(filebase_b,num_i1(ifile),num_j1(ifile),'.png',nom_type)% makes the new file name - imwrite(C,newname,'BitDepth',16); % save the new image - end - %repeat the operation on a sliding series of nbaver*nbfield2 images - 'sliding background image will be substracted' + if test_level + C=levels(C); + imwrite(C,newname,'BitDepth',8); % save the new image + else + imwrite(C,newname,'BitDepth',16); % save the new image + end + end + + %% repeat the operation on a sliding series of nbaver*nbfield2 images + display('sliding background image will be substracted') if nbfield_slice > nbaver_ima -% for ifield = floor(nbaver_ima/2)+2:step:nbfield_slice-floor(nbaver_ima/2) for ifield = step*ceil(nbaver/2)+1:step:nbfield_slice-step*floor(nbaver/2) stopstate=get(hseries.RUN,'BusyAction'); @@ -291,5 +290,4 @@ B=Asort(:,:,rank);%current background image for iburst=1:step -% Acor=double(Ak(:,:,floor(nbaver_ima/2)+iburst-1))-double(B); index=step*floor(nbaver/2)+iburst; Acor=double(Ak(:,:,index))-double(B); @@ -298,8 +296,12 @@ ifile=indselect(ifield+iburst-1); [newname]=... - name_generator(filebase_b,num_i1(ifile),num_j1(ifile),'.png',Series.NomType) -% newname=name_generator(filebase_b,num_i1(indselect(ifield+iburst-1)),num_j1(indselect(ifield+iburst-1)),'.png',nom_type)% makes the new file name - imwrite(C,newname,'BitDepth',16); % save the new image - end + name_generator(filebase_b,num_i1(ifile),num_j1(ifile),'.png',Series.NomType) % makes the new file name + if test_level + C=levels(C); + imwrite(C,newname,'BitDepth',8); % save the new image + else + imwrite(C,newname,'BitDepth',16); % save the new image + end + end else return @@ -308,17 +310,21 @@ end -%substract the background from the last images -% for ifield=nbfield_slice-floor(nbaver_ima/2)+1:nbfield_slice - 'last background image will be substracted' +%% substract the background from the last images + display('last background image will be substracted') ifield=nbfield_slice-(step*ceil(nbaver/2))+1:nbfield_slice; - for ifield=nbfield_slice-(step*floor(nbaver/2))+1:nbfield_slice - index=ifield-nbfield_slice+step*(2*floor(nbaver/2)+1); - Acor=double(Ak(:,:,index))-double(B); - Acor=(Acor>0).*Acor; % put to 0 the negative elements in Acor - C=uint16(Acor); - ifile=indselect(ifield); - newname=name_generator(filebase_b,num_i1(ifile),num_j1(ifile),'.png',nom_type)% makes the new file name - imwrite(C,newname,'BitDepth',16); % save the new image - end + for ifield=nbfield_slice-(step*floor(nbaver/2))+1:nbfield_slice + index=ifield-nbfield_slice+step*(2*floor(nbaver/2)+1); + Acor=double(Ak(:,:,index))-double(B); + Acor=(Acor>0).*Acor; % put to 0 the negative elements in Acor + C=uint16(Acor); + ifile=indselect(ifield); + newname=name_generator(filebase_b,num_i1(ifile),num_j1(ifile),'.png',nom_type)% makes the new file name + if test_level + C=levels(C); + imwrite(C,newname,'BitDepth',8); % save the new image + else + imwrite(C,newname,'BitDepth',16); % save the new image + end + end end @@ -348,2 +354,41 @@ end + +function C=levels(A) +%whos A; +B=double(A(:,:,1)); +windowsize=round(min(size(B,1),size(B,2))/20); +windowsize=floor(windowsize/2)*2+1; +ix=[1/2-windowsize/2:-1/2+windowsize/2];% +%del=np/3; +%fct=exp(-(ix/del).^2); +fct2=cos(ix/(windowsize-1)/2*pi/2); +%Mfiltre=(ones(5,5)/5^2); +%Mfiltre=fct2'; +Mfiltre=fct2'*fct2; +Mfiltre=Mfiltre/(sum(sum(Mfiltre))); + +C=filter2(Mfiltre,B); +C(:,1:windowsize)=C(:,windowsize)*ones(1,windowsize); +C(:,end-windowsize+1:end)=C(:,end-windowsize+1)*ones(1,windowsize); +C(1:windowsize,:)=ones(windowsize,1)*C(windowsize,:); +C(end-windowsize+1:end,:)=ones(windowsize,1)*C(end-windowsize,:); +C=tanh(B./(2*C)); +[n,c]=hist(reshape(C,1,[]),100); +% figure;plot(c,n); + +[m,i]=max(n); +c_max=c(i); +[dummy,index]=sort(abs(c-c(i))); +n=n(index); +c=c(index); +i_select = find(cumsum(n)<0.95*sum(n)); +if isempty(i_select) + i_select = 1:length(c); +end +c_select=c(i_select); +n_select=n(i_select); +cmin=min(c_select); +cmax=max(c_select); +C=(C-cmin)/(cmax-cmin)*256; +C=uint8(C);