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

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

aver_stat modified in the new rationalized form.
OTHER series FCTS NOT OPERATIONAL

File size: 18.2 KB
Line 
1%'sub_background': substract background to an image series, used with series.fig
2%------------------------------------------------------------------------
3% function GUI_input=sub_background(Param)
4%
5%OUTPUT
6% GUI_input=list of options in the GUI series.fig needed for the function
7%
8%INPUT:
9%num_i1: series of first indices i (given from the series interface as first_i:incr_i:last_i, mode and list_pair_civ)
10%num_i2: series of second indices i (given from the series interface as first_i:incr_i:last_i, mode and list_pair_civ)
11%num_j1: series of first indices j (given from the series interface as first_j:incr_j:last_j, mode and list_pair_civ )
12%num_j2: series of second indices j (given from the series interface as first_j:incr_j:last_j, mode and list_pair_civ)
13%Series: Matlab structure containing information set by the series interface
14%       .RootPath: path to the image series
15%       .RootFile: root file name
16%       .FileExt: image file extension
17%       .NomType: nomenclature type for file indexing
18%       .NbSlice: %number of slices defined on the interface
19%----------------------------------------------------------------------
20% Method:
21    %calculate the background image by sorting the luminosity of each point
22    % over a sliding sub-sequence of 'nbaver_ima' images.
23    % The luminosity value of rank 'rank' is selected as the
24    % 'background'. rank=nbimages/2 gives the median value.  Smaller values are appropriate
25    % for a dense set of particles. The extrem value rank=1 gives the true minimum
26    % luminosity, but it can be polluted by noise.
27% Organization of image indices:
28    % The program is working on a series of images, labelled by two indices i and j, given
29    % by the input matlab vectors num_i1 and num_j1 respectively. In the list, j is the fastest increasing index.
30    % The processing can be done in slices (number nbslice), with bursts of
31    % nbfield2 successive images for a given slice (mode 'multilevel')
32    % In the mode 'volume', nbfield2=1 (1 image at each level)
33%
34function GUI_input=sub_background (Param)
35
36%------------------------------------------------------------------------
37%requests for the visibility of input windows in the GUI series  (activated directly by the selection in the menu ACTION)
38if ~exist('Param','var')
39    GUI_input={'RootPath';'on';...
40        'SubDir';'off';... % subdirectory of derived files (PIV fields), ('on' by default)
41        'RootFile';'on';... %root input file name ('on' by default)
42        'FileExt';'on';... %inputf file extension ('on' by default)
43        'NomType';'on';...%type of file indexing ('on' by default)
44        'NbSlice';'on'; ...%nbre of slices ('off' by default)
45        'OutputDirExt';'.sbk';...
46        %'VelTypeMenu';'on';...% menu for selecting the velocity type (civ1,..)('off' by default)
47        %'FieldMenu';'on';...% menu for selecting the velocity field (s) in the input file ('off' by default)
48        %'VelTypeMenu_1';'on';...% menu for selecting the velocity type (civ1,..)('off' by default)
49        %'FieldMenu_1';'on';...% menu for selecting the velocity field (s) in the input file ('off' by default)
50        %'CoordType';...%can use a transform function
51        %'GetObject';...;%can use projection object
52        %'GetMask';...;%can use mask option 
53        %'PARAMETER';'NbSliding';...
54        %'PARAMETER';'VolumeScan';...
55        %'PARAMETER';'RankBrightness';...
56               ''};
57    return %exit the function
58end
59
60%% input parameters
61% read the xml file for batch case
62if ischar(Param) && ~isempty(find(regexp('Param','.xml$')))
63    Param=xml2struct(Param);
64    checkrun=0;
65else %  RUN case: parameters introduced as the input structure Param
66    hseries=guidata(Param.hseries);%handles of the GUI series
67    WaitbarPos=get(hseries.waitbar_frame,'Position');
68    checkrun=1;
69end
70filebase=fullfile(Param.InputTable{1,1},Param.InputTable{1,2},Param.InputTable{1,3});
71dir_images=Param.InputTable{1,1};
72NomType=Param.InputTable{1,4};
73SubDir=Param.InputTable{1,2};
74FileExt=Param.InputTable{1,5};
75[filecell,i1_series,tild,j1_series]=get_file_series(Param);%generates the set of input file names
76if size(filecell,1)>1
77    msgbox_uvmat('WARNING','This function uses only the first input image series')
78    return
79end
80
81%% determine input image type
82FileType=[];%default
83MovieObject=[];
84%FileExt=Series.FileExt;
85
86if isequal(lower(FileExt),'.avi')
87    hhh=which('mmreader');
88    if ~isequal(hhh,'')&& mmreader.isPlatformSupported()
89        MovieObject=mmreader(fullfile(RootPath,SubDir,[RootFile FileExt]));
90        FileType='movie';
91    else
92        FileType='avi';
93    end
94elseif isequal(lower(FileExt),'.vol')
95    FileType='vol';
96else
97    form=imformats(FileExt(2:end));
98    if ~isempty(form)% if the extension corresponds to an image format recognized by Matlab
99        if isequal(NomType,'*');
100            FileType='multimage';
101        else
102            FileType='image';
103        end
104    end
105end
106if isempty(FileType)
107    msgbox_uvmat('ERROR',['invalid file extension ' FileExt ': this function only accepts image or movie input'])
108    return
109end
110
111nbslice_i=Param.NbSlice; %number of slices
112siz=size(i1_series{1});
113nbaver_init=23;%approximate number of images used for the sliding background: to be adjusted later to include an integer number of bursts
114j1=[];%default
115
116%% apply the image rescaling function 'level' (avoid the blinking effects of bright particles)
117answer=msgbox_uvmat('INPUT_Y-N','apply image rescaling function levels.m after sub_background');
118test_level=isequal(answer,'Yes');
119
120%% adjust the proposed number of images in the sliding average to include an integer number of bursts
121if siz(2)~=1
122    nbaver=floor(nbaver_init/siz(1)); % number of bursts used for the sliding background,
123    if isequal(floor(nbaver/2),nbaver)
124        nbaver=nbaver+1;%put the number of burst to an odd number (so the middle burst is defined)
125    end
126    nbaver_init=nbaver*siz(1);%propose by default an integer number of bursts
127end
128
129%% create dir of the new images
130if test_level
131    term='.sbk.lev';
132else
133    term='.sbk';
134end
135SubdirResult=[Param.InputTable{1,2} term];
136try
137    mkdir(fullfile(Param.InputTable{1,1},SubdirResult));
138catch ME
139    msgbox_uvmat('ERROR',['error in creating result directory: ' ME.message]);
140    return
141end
142[xx,msg2] = fileattrib(fullfile(Param.InputTable{1,1},SubdirResult),'+w','g'); %yield writing access (+w) to user group (g)
143if ~strcmp(msg2,'')
144    msgbox_uvmat('ERROR',['pb of permission for ' fullfile(Param.InputTable{1,1},SubdirResult) ': ' msg2])%error message for directory creation
145    return
146end
147
148%% set processing parameters
149prompt = {'Number of images for the sliding background (MUST FIT IN COMPUETER MEMORY)';'The number of positions (laser slices)';'volume scan mode (Yes/No)';...
150    'the luminosity rank chosen to define the background (0.1=for dense particle seeding, 0.5 (median) for sparse particles'};
151dlg_title = ['get (slice by slice) a sliding background and substract to each image, result in subdir ' SubdirResult];
152num_lines= 3;
153def     = { num2str(nbaver_init);num2str(nbslice_i);'No';'0.1'};
154answer = inputdlg(prompt,dlg_title,num_lines,def);
155set(hseries.ParamVal,'String',answer([1 [3:4]]))
156set(hseries.ParamVal,'Visible','on')
157
158nbaver_ima=str2num(answer{1});%number of images for the sliding background
159nbaver=ceil(nbaver_ima/siz(1));%number of bursts for the sliding background
160if isequal(floor(nbaver/2),nbaver)
161    nbaver=nbaver+1;%put the number of burst to an odd number (so the middle burst is defined)
162end
163step=siz(1);%case of bursts: the sliding background is shifted by one burst
164vol_test=answer{3};
165if isequal(vol_test,'Yes')
166    nbfield2=1;%case of volume: no consecutive series at a given level
167    nbslice_i=siz(1);%number of slices
168else
169    nbfield2=siz(1); %nb of consecutive images at each level(burst)
170    if siz(2)>1
171        % nbslice_i=str2num(answer{2})/(num_i1(1,2)-num_i1(1,1));% number of slices
172        nbslice_i=str2num(answer{2})/(i1_series{1}(1,2)-i1_series{1}(1,1));% number of slices
173    else
174        nbslice_i=1;
175    end
176    if ~isequal(floor(nbslice_i),nbslice_i)
177        msgbox_uvmat('ERROR','the number of slices must be a multiple of the i increment')
178        return
179    end
180end
181rank=floor(str2num(answer{4})*nbaver_ima);
182if rank==0
183    rank=1;%rank selected in the sorted image series
184end
185lengthtot=siz(1)*siz(2);
186nbfield=floor(lengthtot/(nbfield2*nbslice_i));%total number of i indexes (adjusted to an integer number of slices)
187nbfield_slice=nbfield*nbfield2;% number of fields per slice
188if nbaver_ima > nbfield*nbfield2
189    msgbox_uvmat('ERROR','number of images in a slice smaller than the proposed number of images for the sliding average')
190    return
191end
192nbfirst=(ceil(nbaver/2))*step;
193if nbfirst>nbaver_ima
194    nbfirst=ceil(nbaver_ima/2);
195    step=1;
196    nbaver=nbaver_ima;
197end
198
199%% prealocate memory for the sliding background
200% first_image=name_generator(filebase,num_i1(1),num_j1(1),FileExt,NomType);
201% Afirst=read_image(first_image,FileType,num_i1(1),MovieObject);
202Afirst=read_image(filecell{1,1},FileType,i1_series{1}(1),MovieObject);
203[npy,npx]=size(Afirst);
204try
205    Ak=zeros(npy,npx,nbaver_ima,'uint16'); %prealocate memory
206    Asort=zeros(npy,npx,nbaver_ima,'uint16'); %prealocate memory
207catch ME
208    msgbox_uvmat('ERROR',ME.message)
209    return
210end
211
212%% update the xml file
213SubDirBase=regexprep(Param.InputTable{1,2},'\..*','');%take the root part of SubDir, before the first dot '.'
214filexml=fullfile(Param.InputTable{1,1},[SubDirBase '.xml']);
215if ~exist(filexml,'file') && exist([filebase '.xml'],'file')% xml inside the image directory
216    copyfile([filebase '.xml'],filexml);% copy the .xml file
217end
218if exist(filexml,'file')
219    t=xmltree(filexml); 
220    %update information on the first image name in the series
221    uid_Heading=find(t,'ImaDoc/Heading');
222    if isempty(uid_Heading)
223        [t,uid_Heading]=add(t,1,'element','Heading');
224    end
225    uid_ImageName=find(t,'ImaDoc/Heading/ImageName');
226    if ~isempty(j1_series{1})
227        j1=j1_series{1}(1);
228    end
229    ImageName=fullfile_uvmat([dir_images term],'',Param.InputTable{1,3},'.png',NomType,i1_series{1}(1),[],j1);
230    %     ImageName=name_generator(filebase_b,num_i1(1),num_j1(1),'.png',NomType);
231    [pth,ImageName]=fileparts(ImageName);
232    ImageName=[ImageName '.png'];
233    if isempty(uid_ImageName)
234        [t,uid_ImageName]=add(t,uid_Heading,'element','ImageName');
235    end
236    uid_value=children(t,uid_ImageName);
237    if isempty(uid_value)
238        t=add(t,uid_ImageName,'chardata',ImageName);%indicate  name of the first image, with ;png extension
239    else
240        t=set(t,uid_value(1),'value',ImageName);%indicate  name of the first image, with ;png extension
241    end
242   
243    %add information about image transform
244    [t,new_uid]=add(t,1,'element','ImageTransform');
245    [t,NameFunction_uid]=add(t,new_uid,'element','NameFunction');
246    [t]=add(t,NameFunction_uid,'chardata','sub_background');
247    if test_level
248        [t,NameFunction_uid]=add(t,new_uid,'element','NameFunction');
249        [t]=add(t,NameFunction_uid,'chardata','levels');
250    end
251    [t,NbSlice_uid]=add(t,new_uid,'element','NbSlice');
252    [t]=add(t,new_uid,'chardata',num2str(nbslice_i));
253    [t,NbSlidingImages_uid]=add(t,new_uid,'element','NbSlidingImages');
254    [t]=add(t,NbSlidingImages_uid,'chardata',num2str(nbaver));
255    [t,LuminosityRank_uid]=add(t,new_uid,'element','RankBackground');
256    [t]=add(t,LuminosityRank_uid,'chardata',num2str(rank));% luminosity rank almong the nbaver sliding images
257    save(t,filexml)
258end
259%copy the mask
260% if exist([filebase '_1mask_1'],'file')
261%     copyfile([filebase '_1mask_1'],[filebase_b '_1mask_1']);% copy the mask file
262% end
263
264%MAIN LOOP ON SLICES
265for islice=1:nbslice_i
266    %% select the series of image indices at the level islice
267    for ifield=1:nbfield
268        for iburst=1:nbfield2
269            indselect(iburst,ifield)=((ifield-1)*nbslice_i+(islice-1))*nbfield2+iburst;
270        end
271    end
272   
273    %% read the first series of nbaver_ima images and sort by luminosity at each pixel
274    for ifield = 1:nbaver_ima
275        ifile=indselect(ifield);
276        filename=filecell{1,ifile};
277        % filename=name_generator(filebase,num_i1(ifile),num_j1(ifile),FileExt,NomType);
278        Aread=read_image(filename,FileType,i1_series{1}(ifile),MovieObject);
279        Ak(:,:,ifield)=Aread;
280    end
281    Asort=sort(Ak,3);%sort the luminosity of images at each point
282    B=Asort(:,:,rank);%background image
283    display( 'first background image will be substracted')
284    for ifield=1:nbfirst
285        Acor=double(Ak(:,:,ifield))-double(B);%substract background to the current image
286        Acor=(Acor>0).*Acor; % put to 0 the negative elements in Acor
287        C=uint16(Acor);% set to integer 16 bits
288        ifile=indselect(ifield);
289        %             newname=name_generator(filebase_b,num_i1(ifile),num_j1(ifile),'.png',NomType)% makes the new file name
290        if ~isempty(j1_series{1})
291            j1=j1_series{1}(ifile);
292        end
293        newname=fullfile_uvmat(Param.InputTable{1,1},SubdirResult,Param.InputTable{1,3},'.png',NomType,i1_series{1}(ifile),[],j1);
294        %newname=name_generator(filebase_b,i1_series{1}(ifile),j1_series{1}(ifile),'.png',NomType);% makes the new file name
295        if test_level
296            C=levels(C);
297            imwrite(C,newname,'BitDepth',8); % save the new image
298        else
299            imwrite(C,newname,'BitDepth',16); % save the new image
300        end
301    end
302   
303    %% repeat the operation on a sliding series of nbaver*nbfield2 images
304    display('sliding background image will be substracted')
305    if nbfield_slice > nbaver_ima
306        for ifield = step*ceil(nbaver/2)+1:step:nbfield_slice-step*floor(nbaver/2)
307            if checkrun
308                stopstate=get(hseries.RUN,'BusyAction');
309                update_waitbar(hseries.waitbar_frame,WaitbarPos,(ifield+(islice-1)*nbfield_slice)/(nbfield_slice*nbslice_i))
310            else
311                stopstate='queue';
312            end
313            if isequal(stopstate,'queue')% enable STOP command
314                Ak(:,:,1:nbaver_ima-step)=Ak(:,:,1+step:nbaver_ima);% shift the current image series by one burst (step)
315                %incorporate next burst in the current image series
316                for iburst=1:step
317                    ifile=indselect(ifield+step*floor(nbaver/2)+iburst-1);
318                    filename=fullfile_uvmat(Param.InputTable{1,1},SubDir,Param.InputTable{1,3},FileExt,NomType,i1_series{1}(ifile),[],j1_series{1}(ifile));
319                    %filename=name_generator(filebase,num_i1(ifile),num_j1(ifile),FileExt,NomType);
320                    Aread=read_image(filename,FileType,i1_series{1}(ifile),MovieObject);
321                    Ak(:,:,nbaver_ima-step+iburst)=Aread;
322                end
323                Asort=sort(Ak,3);%sort the new current image series by luminosity
324                B=Asort(:,:,rank);%current background image
325                for iburst=1:step
326                    index=step*floor(nbaver/2)+iburst;
327                    Acor=double(Ak(:,:,index))-double(B);
328                    Acor=(Acor>0).*Acor; % put to 0 the negative elements in Acor
329                    C=uint16(Acor);
330                    ifile=indselect(ifield+iburst-1);
331                    if ~isempty(j1_series{1})
332                        j1=j1_series{1}(ifile);
333                    end
334                    newname=fullfile_uvmat(Param.InputTable{1,1},SubdirResult,Param.InputTable{1,3},'.png',NomType,i1_series{1}(ifile),[],j1);
335                    %[newname]=name_generator(filebase_b,num_i1(ifile),num_j1(ifile),'.png',NomType) % makes the new file name
336                    if test_level
337                        C=levels(C);
338                        imwrite(C,newname,'BitDepth',8); % save the new image
339                    else
340                        imwrite(C,newname,'BitDepth',16); % save the new image
341                    end
342                end
343            else
344                return
345            end
346        end
347    end
348   
349    %% substract the background from the last images
350    display('last background image will be substracted')
351    ifield=nbfield_slice-(step*ceil(nbaver/2))+1:nbfield_slice;
352    for ifield=nbfield_slice-(step*floor(nbaver/2))+1:nbfield_slice
353        index=ifield-nbfield_slice+step*(2*floor(nbaver/2)+1);
354        Acor=double(Ak(:,:,index))-double(B);
355        Acor=(Acor>0).*Acor; % put to 0 the negative elements in Acor
356        C=uint16(Acor);
357        ifile=indselect(ifield);
358        if ~isempty(j1_series{1})
359            j1=j1_series{1}(ifile);
360        end
361        newname=fullfile_uvmat(Param.InputTable{1,1},SubdirResult,Param.InputTable{1,3},'.png',NomType,i1_series{1}(ifile),[],j1);
362        if test_level
363            C=levels(C);
364            imwrite(C,newname,'BitDepth',8); % save the new image
365        else
366            imwrite(C,newname,'BitDepth',16); % save the new image
367        end
368    end
369end
370
371%finish the waitbar
372if checkrun
373    update_waitbar(hseries.waitbar,WaitbarPos,1)
374end
375
376%------------------------------------------------------------------------
377%--read images and convert them to the uint16 format used for PIV
378function A=read_image(filename,type_ima,num,MovieObject)
379%------------------------------------------------------------------------
380%num is the view number needed for an avi movie
381switch type_ima
382    case 'movie'
383        A=read(MovieObject,num);
384    case 'avi'
385        mov=aviread(filename,num);
386        A=frame2im(mov(1));
387    case 'multimage'
388        A=imread(filename,num);
389    case 'image'   
390        A=imread(filename);
391end
392siz=size(A);
393if length(siz)==3;%color images
394    A=sum(double(A),3);
395end
396   
397
398function C=levels(A)
399%whos A;
400B=double(A(:,:,1));
401windowsize=round(min(size(B,1),size(B,2))/20);
402windowsize=floor(windowsize/2)*2+1;
403ix=1/2-windowsize/2:-1/2+windowsize/2;%
404%del=np/3;
405%fct=exp(-(ix/del).^2);
406fct2=cos(ix/(windowsize-1)/2*pi/2);
407%Mfiltre=(ones(5,5)/5^2);
408%Mfiltre=fct2';
409Mfiltre=fct2'*fct2;
410Mfiltre=Mfiltre/(sum(sum(Mfiltre)));
411
412C=filter2(Mfiltre,B);
413C(:,1:windowsize)=C(:,windowsize)*ones(1,windowsize);
414C(:,end-windowsize+1:end)=C(:,end-windowsize+1)*ones(1,windowsize);
415C(1:windowsize,:)=ones(windowsize,1)*C(windowsize,:);
416C(end-windowsize+1:end,:)=ones(windowsize,1)*C(end-windowsize,:);
417C=tanh(B./(2*C));
418[n,c]=hist(reshape(C,1,[]),100);
419% figure;plot(c,n);
420
421[m,i]=max(n);
422c_max=c(i);
423[dummy,index]=sort(abs(c-c(i)));
424n=n(index);
425c=c(index);
426i_select = find(cumsum(n)<0.95*sum(n));
427if isempty(i_select)
428    i_select = 1:length(c);
429end
430c_select=c(i_select);
431n_select=n(i_select);
432cmin=min(c_select);
433cmax=max(c_select);
434C=(C-cmin)/(cmax-cmin)*256;
435C=uint8(C);
Note: See TracBrowser for help on using the repository browser.