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

Last change on this file since 454 was 454, checked in by sommeria, 12 years ago

aver_stat, ima_levels , sub_background now adapted to the new settings
merge_proj, time_series still need to be modified

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