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

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

fgunctions under series adapted to the renovated GUI series

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