source: trunk/src/series/sub_background.m @ 624

Last change on this file since 624 was 624, checked in by sommeria, 11 years ago

waitbar system for series improved to aloow use as stand alone fcts.

to add at the head of series fcts:
hseries=findobj(allchild(0),'Tag','series');
RUNHandle=findobj(hseries,'Tag','RUN');%handle of RUN button in GUI series
WaitbarHandle?=findobj(hseries,'Tag','Waitbar');%handle of waitbar in GUI series

call to waitbar:

update_waitbar(WaitbarHandle?,index/nbfield)
if ishandle(RUNHandle) && ~strcmp(get(RUNHandle,'BusyAction?'),'queue')

disp('program stopped by user')
break

end

File size: 20.9 KB
Line 
1%'sub_background': substract a sliding background to an image series
2%------------------------------------------------------------------------
3% Method:
4    %calculate the background image by sorting the luminosity of each point
5    % over a sliding sub-sequence of 'nbaver_ima' images.
6    % The luminosity value of rank 'rank' is selected as the
7    % 'background'. rank=nbimages/2 gives the median value.  Smaller values are appropriate
8    % for a dense set of particles. The extrem value rank=1 gives the true minimum
9    % luminosity, but it can be polluted by noise.
10% Organization of image indices:
11    % The program is working on a series of images,
12    % In the mode 'volume', nbfield2=1 (1 image at each level)and NbSlice (=nbfield_j)
13    % Else nbfield2=nbfield_j =nbre of images in a burst (j index)
14   
15% function GUI_config=sub_background(Param)
16%
17%%%%%%%%%%% GENERAL TO ALL SERIES ACTION FCTS %%%%%%%%%%%%%%%%%%%%%%%%%%%
18%
19%OUTPUT
20% ParamOut: sets options in the GUI series.fig needed for the function
21%
22%INPUT:
23% In run mode, the input parameters are given as a Matlab structure Param copied from the GUI series.
24% In batch mode, Param is the name of the corresponding xml file containing the same information
25% when Param.Action.RUN=0 (as activated when the current Action is selected
26% in series), the function ouput paramOut set the activation of the needed GUI elements
27%
28% Param contains the elements:(use the menu bar command 'export/GUI config' in series to
29% see the current structure Param)
30%    .InputTable: cell of input file names, (several lines for multiple input)
31%                      each line decomposed as {RootPath,SubDir,Rootfile,NomType,Extension}
32%    .OutputSubDir: name of the subdirectory for data outputs
33%    .OutputDirExt: directory extension for data outputs
34%    .Action: .ActionName: name of the current activated function
35%             .ActionPath:   path of the current activated function
36%             .ActionExt: fct extension ('.m', Matlab fct, '.sh', compiled   Matlab fct
37%             .RUN =0 for GUI input, =1 for function activation
38%             .RunMode='local','background', 'cluster': type of function  use
39%             
40%    .IndexRange: set the file or frame indices on which the action must be performed
41%    .FieldTransform: .TransformName: name of the selected transform function
42%                     .TransformPath:   path  of the selected transform function
43%    .InputFields: sub structure describing the input fields withfields
44%              .FieldName: name(s) of the field
45%              .VelType: velocity type
46%              .FieldName_1: name of the second field in case of two input series
47%              .VelType_1: velocity type of the second field in case of two input series
48%              .Coord_y: name of y coordinate variable
49%              .Coord_x: name of x coordinate variable
50%    .ProjObject: %sub structure describing a projection object (read from ancillary GUI set_object)
51%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
52   
53function ParamOut=sub_background (Param)
54
55%% input preparation mode (no RUN)
56if isstruct(Param) && isequal(Param.Action.RUN,0)
57    ParamOut.AllowInputSort='off';% allow alphabetic sorting of the list of input file SubDir (options 'off'/'on', 'off' by default)
58    ParamOut.WholeIndexRange='on';% prescribes the file index ranges from min to max (options 'off'/'on', 'off' by default)
59    ParamOut.NbSlice='on'; %nbre of slices ('off' by default)
60    ParamOut.VelType='off';% menu for selecting the velocity type (options 'off'/'one'/'two',  'off' by default)
61    ParamOut.FieldName='off';% menu for selecting the field (s) in the input file(options 'off'/'one'/'two', 'off' by default)
62    ParamOut.FieldTransform = 'off';%can use a transform function
63    ParamOut.ProjObject='off';%can use projection object(option 'off'/'on',
64    ParamOut.Mask='off';%can use mask option   (option 'off'/'on', 'off' by default)
65    ParamOut.OutputDirExt='.sback';%set the output dir extension
66    ParamOut.OutputFileMode='NbInput';% '=NbInput': 1 output file per input file index, '=NbInput_i': 1 file per input file index i, '=NbSlice': 1 file per slice
67    filecell=get_file_series(Param);%check existence of the first input file
68    %%%%%%%%%%%% STANDARD PART (DO NOT EDIT) %%%%%%%%%%%%
69   
70    %% root input file(s) and type
71    [filecell,i1_series,i2_series,j1_series,j2_series]=get_file_series(Param);
72    if ~exist(filecell{1,1},'file')
73        msgbox_uvmat('WARNING','the first input file does not exist')
74        return
75    end
76    %%%%%%%%%%%%
77    % The cell array filecell is the list of input file names, while
78    % filecell{iview,fileindex}:
79    %        iview: line in the table corresponding to a given file series
80    %        fileindex: file index within  the file series,
81    % i1_series(iview,ref_j,ref_i)... are the corresponding arrays of indices i1,i2,j1,j2, depending on the input line iview and the two reference indices ref_i,ref_j
82    % i1_series(iview,fileindex) expresses the same indices as a 1D array in file indices
83    %%%%%%%%%%%%
84    NbSlice=1;%default
85    if isfield(Param.IndexRange,'NbSlice')&&~isempty(Param.IndexRange.NbSlice)
86        NbSlice=Param.IndexRange.NbSlice;
87    end
88    nbview=numel(i1_series);%number of input file series (lines in InputTable)
89    nbfield_j=size(i1_series{1},1); %nb of fields for the j index (bursts or volume slices)
90    nbfield_i=size(i1_series{1},2); %nb of fields for the i index
91    nbfield=nbfield_j*nbfield_i; %total number of fields
92    nbfield_i=floor(nbfield/NbSlice);%total number of  indexes in a slice (adjusted to an integer number of slices)
93    nbfield=nbfield_i*NbSlice; %total number of fields after adjustement
94   
95    %determine the file type on each line from the first input file
96
97   
98    %% calibration data and timing: read the ImaDoc files
99    %not relevant here
100   
101    %% check coincidence in time for several input file series
102    %not relevant here
103   
104    %% coordinate transform or other user defined transform
105    %not relevant here
106   
107    %%%%%%%%%%%% END STANDARD PART  %%%%%%%%%%%%
108    % EDIT FROM HERE
109   
110    %% check the validity of  input file types
111    ImageTypeOptions={'image','multimage','mmreader','video'};%allowed input file types(images)
112    FileType=get_file_type(filecell{1,1});
113    CheckImage=~isempty(find(strcmp(FileType,ImageTypeOptions), 1));% =1 for images
114    if ~CheckImage
115        msgbox_uvmat('ERROR',['invalid file type input: ' FileType{1} ' not an image'])
116        return
117    end
118   
119    %% Set field names and velocity types
120    %not relevant here
121   
122    %% Initiate output fields
123    %not relevant here
124   
125    %%% SPECIFIC PART BEGINS HERE
126    NbSlice=1;
127    if isfield(Param.IndexRange,'NbSlice')
128        NbSlice=Param.IndexRange.NbSlice; %number of slices
129    end
130    %siz=size(i1_series);
131    nbaver_init=23;%approximate number of images used for the sliding background: to be adjusted later to include an integer number of bursts
132    j1=[];%default
133   
134    %% adjust the proposed number of images in the sliding average to include an integer number of bursts
135    if nbfield_i~=1
136        nbaver=floor(nbaver_init/nbfield_j); % number of bursts used for the sliding background,
137        if isequal(floor(nbaver/2),nbaver)
138            nbaver=nbaver+1;%put the number of burst to an odd number (so the middle burst is defined)
139        end
140        nbaver_init=nbaver*nbfield_j;%propose by default an integer number of bursts
141    end
142   
143    %% input of specific parameters
144    %if checkrun %get specific parameters interactively
145    if isequal(Param.Action.RUN,0)
146        prompt = {'volume scan mode (Yes/No)';'Number of images for the sliding background (MUST FIT IN COMPUTER MEMORY)';...
147            'the luminosity rank chosen to define the background (0.1=for dense particle seeding, 0.5 (median) for sparse particles'};
148        dlg_title = 'get (slice by slice) a sliding background and substract to each image';
149        num_lines= 3;
150        def     = { 'No';num2str(nbaver_init);'0.1'};
151        answer = inputdlg(prompt,dlg_title,num_lines,def);
152       
153        %check input consistency
154        if strcmp(answer{1},'No') && ~isequal(NbSlice,1)
155            check=msgbox_uvmat('INPUT_Y-N',['confirm the multi-level splitting into ' num2str(NbSlice) ' slices']);
156            if ~strcmp(check,'Yes')
157                return
158            end
159        end
160        if strcmp(answer{1},'Yes')
161            step=1;
162        else
163            step=nbfield_j;%case of bursts: the sliding background is shifted by the length of one burst
164        end
165        nbaver_ima=str2num(answer{2});%number of images for the sliding background
166        nbaver=ceil(nbaver_ima/step);%number of bursts for the sliding background
167        if isequal(floor(nbaver/2),nbaver)
168            nbaver=nbaver+1;%set the number of bursts to an odd number (so the middle burst is defined)
169        end
170        nbaver_ima=nbaver*step;
171    end
172    ParamOut.ActionInput.CheckVolume=strcmp(answer{1},'Yes');
173    ParamOut.ActionInput.SlidingSequenceLength=nbaver_ima;
174    ParamOut.ActionInput.BrightnessRankThreshold=str2num(answer{3});
175   
176    % apply the image rescaling function 'level' (avoid the blinking effects of bright particles)
177    answer=msgbox_uvmat('INPUT_Y-N','apply image rescaling function levels.m after sub_background');
178    ParamOut.ActionInput.CheckLevelTransform=strcmp(answer,'Yes');
179    return
180end
181%%%%%%%%%%%%%%%%%%%%%%  STOP HERE FOR PAMETER INPUT MODE  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
182
183%% read input parameters from an xml file if input is a file name (batch mode)
184checkrun=1;
185if ischar(Param)
186    Param=xml2struct(Param);% read Param as input file (batch case)
187    checkrun=0;
188end
189hseries=findobj(allchild(0),'Tag','series');
190RUNHandle=findobj(hseries,'Tag','RUN');%handle of RUN button in GUI series
191WaitbarHandle=findobj(hseries,'Tag','Waitbar');%handle of waitbar in GUI series
192%% Input preparation
193nbaver_ima=Param.ActionInput.SlidingSequenceLength;
194NbSlice=Param.IndexRange.NbSlice;
195if ~isequal(NbSlice,1)
196    display(['multi-level splitting into ' num2str(NbSlice) ' slices']);
197end
198RootPath=Param.InputTable(:,1);
199RootFile=Param.InputTable(:,3);
200SubDir=Param.InputTable(:,2);
201NomType=Param.InputTable(:,4);
202FileExt=Param.InputTable(:,5);
203[filecell,i1_series,i2_series,j1_series,j2_series]=get_file_series(Param);
204[FileType{1},FileInfo{1},MovieObject{1}]=get_file_type(filecell{1,1});
205    if ~isempty(j1_series{1})
206        frame_index{1}=j1_series{1};
207    else
208        frame_index{1}=i1_series{1};
209    end
210nbfield_j=size(i1_series{1},1); %nb of fields for the j index (bursts or volume slices)
211nbfield_i=size(i1_series{1},2); %nb of fields for the i index
212nbfield=nbfield_j*nbfield_i; %total number of fields
213nbfield_i=floor(nbfield/NbSlice);%total number of  indexes in a slice (adjusted to an integer number of slices)
214nbfield=nbfield_i*NbSlice; %total number of fields after adjustement
215
216%% Output
217FileExtOut='.png'; % write result as .png images for image inputs
218if strcmp(lower(NomType{1}(end)),'a')
219    NomTypeOut=NomType{1};%case of letter appendix
220elseif isempty(j1_series)
221    NomTypeOut='_1';
222else
223    NomTypeOut='_1_1';% caseof purely numerical indexing
224end
225
226OutputDir=[Param.OutputSubDir Param.OutputDirExt];
227
228if isequal(Param.ActionInput.CheckVolume,1)
229    step=1;
230else
231    step=nbfield_j;%case of bursts: the sliding background is shifted by the length of one burst
232end
233nbaver_ima=Param.ActionInput.SlidingSequenceLength;%number of images for the sliding background
234nbaver=ceil(nbaver_ima/step);%number of bursts for the sliding background
235if isequal(floor(nbaver/2),nbaver)
236    nbaver=nbaver+1;%set the number of bursts to an odd number (so the middle burst is defined)
237end
238nbaver_ima=nbaver*step;
239if nbaver_ima > nbfield
240    display('number of images in a slice smaller than the proposed number of images for the sliding average')
241    return
242end
243
244% calculate absolute brightness rank
245rank=floor(Param.ActionInput.BrightnessRankThreshold*nbaver_ima);
246if rank==0
247    rank=1;%rank selected in the sorted image series
248end
249
250%% prealocate memory for the sliding background
251try
252    Afirst=read_image(filecell{1,1},FileType{1},MovieObject{1},frame_index{1}(1));
253    [npy,npx]=size(Afirst);
254    if strcmp(class(Afirst),'uint8') % case of 8bit images
255        Ak=zeros(npy,npx,nbaver_ima,'uint8'); %prealocate memory
256        Asort=zeros(npy,npx,nbaver_ima,'uint8'); %prealocate memory
257    else
258        Ak=zeros(npy,npx,nbaver_ima,'uint16'); %prealocate memory
259        Asort=zeros(npy,npx,nbaver_ima,'uint16'); %prealocate memory
260    end
261catch ME
262    msgbox_uvmat('ERROR',['sub_background/read_image/' ME.message])
263    return
264end
265
266%% update the xml file
267% SubDirBase=regexprep(Param.InputTable{1,2},'\..*','');%take the root part of SubDir, before the first dot '.'
268% filexml=fullfile(RootPath{1},[SubDirBase '.xml']);
269% if ~exist(filexml,'file') && exist([fullfile(RootPath{1},SubDir{1},RootFile{1}) '.xml'],'file')% xml inside the image directory
270%     copyfile([filebase '.xml'],filexml);% copy the .xml file
271% end
272% if exist(filexml,'file')
273%     t=xmltree(filexml); 
274%     %update information on the first image name in the series
275%     uid_Heading=find(t,'ImaDoc/Heading');
276%     if isempty(uid_Heading)
277%         [t,uid_Heading]=add(t,1,'element','Heading');
278%     end
279%     uid_ImageName=find(t,'ImaDoc/Heading/ImageName');
280%     if ~isempty(j1_series{1})
281%         j1=j1_series{1}(1);
282%     end
283%     ImageName=fullfile_uvmat([dir_images term],'',RootFile{1},'.png',NomType,i1_series(1,1),[],j1);
284%     [pth,ImageName]=fileparts(ImageName);
285%     ImageName=[ImageName '.png'];
286%     if isempty(uid_ImageName)
287%         [t,uid_ImageName]=add(t,uid_Heading,'element','ImageName');
288%     end
289%     uid_value=children(t,uid_ImageName);
290%     if isempty(uid_value)
291%         t=add(t,uid_ImageName,'chardata',ImageName);%indicate  name of the first image, with ;png extension
292%     else
293%         t=set(t,uid_value(1),'value',ImageName);%indicate  name of the first image, with ;png extension
294%     end
295%     
296%     %add information about image transform
297%     [t,new_uid]=add(t,1,'element','ImageTransform');
298%     [t,NameFunction_uid]=add(t,new_uid,'element','NameFunction');
299%     [t]=add(t,NameFunction_uid,'chardata','sub_background');
300%     if GUI_config.CheckLevel
301%         [t,NameFunction_uid]=add(t,new_uid,'element','NameFunction');
302%         [t]=add(t,NameFunction_uid,'chardata','levels');
303%     end
304%     [t,NbSlice_uid]=add(t,new_uid,'element','NbSlice');
305%     [t]=add(t,new_uid,'chardata',num2str(NbSlice));
306%     [t,NbSlidingImages_uid]=add(t,new_uid,'element','NbSlidingImages');
307%     [t]=add(t,NbSlidingImages_uid,'chardata',num2str(nbaver));
308%     [t,LuminosityRank_uid]=add(t,new_uid,'element','RankBackground');
309%     [t]=add(t,LuminosityRank_uid,'chardata',num2str(rank));% luminosity rank almong the nbaver sliding images
310%     save(t,filexml)
311% end
312%copy the mask
313% if exist([filebase '_1mask_1'],'file')
314%     copyfile([filebase '_1mask_1'],[filebase_b '_1mask_1']);% copy the mask file
315% end
316
317%MAIN LOOP ON SLICES
318% for islice=1:NbSlice
319%% select the series of image indices
320indselect=1:step:nbfield;% select file indices of the slice
321for ifield=1:step-1
322    indselect=[indselect;indselect(end,:)+1];
323end
324
325%% read the first series of nbaver_ima images and sort by luminosity at each pixel
326for ifield = 1:nbaver_ima
327    ifile=indselect(ifield);
328    filename=filecell{1,ifile};
329    Aread=read_image(filename,FileType{1},MovieObject{1},frame_index{1}(ifile));
330    if ndims(Aread)==3;%color images
331        Aread=sum(double(Aread),3);% take the sum of color components
332    end
333    Ak(:,:,ifield)=Aread;
334end
335Asort=sort(Ak,3);%sort the luminosity of images at each point
336B=Asort(:,:,rank);%background image
337display( 'first background image will be substracted')
338nbfirst=(ceil(nbaver/2))*step;
339for ifield=1:nbfirst
340    Acor=double(Ak(:,:,ifield))-double(B);%substract background to the current image
341    Acor=(Acor>0).*Acor; % put to 0 the negative elements in Acor
342    ifile=indselect(ifield);
343    if ~isempty(j1_series{1})
344        j1=j1_series{1}(ifile);
345    end
346    newname=fullfile_uvmat(RootPath{1},OutputDir,RootFile{1},FileExtOut,NomTypeOut,i1_series{1}(ifile),[],j1);
347   
348    %write result file
349    if Param.ActionInput.CheckLevelTransform
350        C=levels(Acor);
351        imwrite(C,newname,'BitDepth',8); % save the new image
352    else
353        if isequal(FileInfo{1}.BitDepth,16)
354            C=uint16(Acor);
355            imwrite(C,newname,'BitDepth',16); % save the new image
356        else
357            C=uint8(Acor);
358            imwrite(C,newname,'BitDepth',8); % save the new image
359        end
360    end
361    display([newname ' written'])
362end
363
364%% repeat the operation on a sliding series of images
365display('sliding background image will be substracted')
366if nbfield_i > nbaver_ima
367    for ifield = step*ceil(nbaver/2)+1:step:nbfield_i-step*floor(nbaver/2)
368                update_waitbar(WaitbarHandle,ifield/nbfield_i)
369    if ishandle(RUNHandle) && ~strcmp(get(RUNHandle,'BusyAction'),'queue')
370        disp('program stopped by user')
371        return
372    end
373%         if isequal(stopstate,'queue')% enable STOP command
374            Ak(:,:,1:nbaver_ima-step)=Ak(:,:,1+step:nbaver_ima);% shift the current image series by one burst (step)
375            %incorporate next burst in the current image series
376            for iburst=1:step
377                ifile=indselect(ifield+step*floor(nbaver/2)+iburst-1);
378                filename=fullfile_uvmat(RootPath{1},SubDir{1},RootFile{1},FileExt{1},NomType{1},i1_series{1}(ifile),[],j1_series{1}(ifile));
379                Aread=read_image(filename,FileType{1},MovieObject{1},i1_series{1}(ifile));
380                if ndims(Aread)==3;%color images
381                    Aread=sum(double(Aread),3);% take the sum of color components
382                end
383                Ak(:,:,nbaver_ima-step+iburst)=Aread;
384            end
385            Asort=sort(Ak,3);%sort the new current image series by luminosity
386            B=Asort(:,:,rank);%current background image
387            for iburst=1:step
388                index=step*floor(nbaver/2)+iburst;
389                Acor=double(Ak(:,:,index))-double(B);
390                Acor=(Acor>0).*Acor; % put to 0 the negative elements in Acor
391                ifile=indselect(ifield+iburst-1);
392                if ~isempty(j1_series{1})
393                    j1=j1_series{1}(ifile);
394                end
395                newname=fullfile_uvmat(RootPath{1},OutputDir,RootFile{1},FileExtOut,NomTypeOut,i1_series{1}(ifile),[],j1);
396                %write result file
397                if Param.ActionInput.CheckLevelTransform
398                    C=levels(Acor);
399                    imwrite(C,newname,'BitDepth',8); % save the new image
400                else
401                    if isequal(FileInfo{1}.BitDepth,16)
402                        C=uint16(Acor);
403                        imwrite(C,newname,'BitDepth',16); % save the new image
404                    else
405                        C=uint8(Acor);
406                        imwrite(C,newname,'BitDepth',8); % save the new image
407                    end
408                end
409                display([newname ' written'])
410               
411            end
412%         else
413%             return
414%         end
415    end
416end
417
418%% substract the background from the last images
419display('last background image will be substracted')
420ifield=nbfield_i-(step*ceil(nbaver/2))+1:nbfield_i;
421for ifield=nbfield_i-(step*floor(nbaver/2))+1:nbfield_i
422    index=ifield-nbfield_i+step*(2*floor(nbaver/2)+1);
423    Acor=double(Ak(:,:,index))-double(B);
424    Acor=(Acor>0).*Acor; % put to 0 the negative elements in Acor
425    ifile=indselect(ifield);
426    if ~isempty(j1_series{1})
427        j1=j1_series{1}(ifile);
428    end
429    newname=fullfile_uvmat(RootPath{1},OutputDir,RootFile{1},FileExtOut,NomTypeOut,i1_series{1}(ifile),[],j1);
430   
431    %write result file
432    if Param.ActionInput.CheckLevelTransform
433        C=levels(Acor);
434        imwrite(C,newname,'BitDepth',8); % save the new image
435    else
436        if isequal(FileInfo{1}.BitDepth,16)
437            C=uint16(Acor);
438            imwrite(C,newname,'BitDepth',16); % save the new image
439        else
440            C=uint8(Acor);
441            imwrite(C,newname,'BitDepth',8); % save the new image
442        end
443    end
444    display([newname ' written'])
445end
446
447
448
449function C=levels(A)
450%whos A;
451B=double(A(:,:,1));
452windowsize=round(min(size(B,1),size(B,2))/20);
453windowsize=floor(windowsize/2)*2+1;
454ix=1/2-windowsize/2:-1/2+windowsize/2;%
455%del=np/3;
456%fct=exp(-(ix/del).^2);
457fct2=cos(ix/(windowsize-1)/2*pi/2);
458%Mfiltre=(ones(5,5)/5^2);
459%Mfiltre=fct2';
460Mfiltre=fct2'*fct2;
461Mfiltre=Mfiltre/(sum(sum(Mfiltre)));
462
463C=filter2(Mfiltre,B);
464C(:,1:windowsize)=C(:,windowsize)*ones(1,windowsize);
465C(:,end-windowsize+1:end)=C(:,end-windowsize+1)*ones(1,windowsize);
466C(1:windowsize,:)=ones(windowsize,1)*C(windowsize,:);
467C(end-windowsize+1:end,:)=ones(windowsize,1)*C(end-windowsize,:);
468C=tanh(B./(2*C));
469[n,c]=hist(reshape(C,1,[]),100);
470% figure;plot(c,n);
471
472[m,i]=max(n);
473c_max=c(i);
474[dummy,index]=sort(abs(c-c(i)));
475n=n(index);
476c=c(index);
477i_select = find(cumsum(n)<0.95*sum(n));
478if isempty(i_select)
479    i_select = 1:length(c);
480end
481c_select=c(i_select);
482n_select=n(i_select);
483cmin=min(c_select);
484cmax=max(c_select);
485C=(C-cmin)/(cmax-cmin)*256;
486C=uint8(C);
Note: See TracBrowser for help on using the repository browser.