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

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

series fcts updated to fit with new waitbar fct and background run mode, and various bug repairs

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