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

Last change on this file since 190 was 169, checked in by sommeria, 13 years ago

comments added at the head of functions

File size: 15.6 KB
Line 
1%'sub_background': substract background to an image series, used with series.fig
2%------------------------------------------------------------------------
3% function GUI_input=aver_stat(num_i1,num_i2,num_j1,num_j2,Series)
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 (num_i1,num_i2,num_j1,num_j2,Series)
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('num_i1','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%----------------------------------------------------------------
60% initiate the waitbar
61hseries=guidata(Series.hseries);%handles of the GUI series
62WaitbarPos=get(hseries.waitbar_frame,'Position');
63%-----------------------------------------------------------------
64if iscell(Series.RootPath)
65    msgbox_uvmat('ERROR','This function use only one input image series')
66    return
67end
68
69%determine input image type
70FileType=[];%default
71MovieObject=[];
72FileExt=Series.FileExt;
73
74if isequal(lower(FileExt),'.avi')
75    hhh=which('mmreader');
76    if ~isequal(hhh,'')&& mmreader.isPlatformSupported()
77        MovieObject=mmreader(fullfile(RootPath,[RootFile FileExt]));
78        FileType='movie';
79    else
80        FileType='avi';
81    end
82elseif isequal(lower(FileExt),'.vol')
83    FileType='vol';
84else
85   form=imformats(FileExt(2:end));
86   if ~isempty(form)% if the extension corresponds to an image format recognized by Matlab
87       if isequal(Series.NomType,'*');
88           FileType='multimage';
89       else
90           FileType='image';
91       end
92   end
93end
94if isempty(FileType)
95    msgbox_uvmat('ERROR',['invalid file extension ' FileExt ': this function only accepts image or movie input'])
96    return
97end
98
99nbslice_i=Series.NbSlice; %number of slices
100siz=size(num_i1);
101nbaver_init=23;%approximate number of images used for the sliding background: to be adjusted later to include an integer number of bursts
102
103%adjust the proposed number of images in the sliding average to include an integer number of bursts
104if siz(2)~=1
105    nbaver=floor(nbaver_init/siz(1)); % number of bursts used for the sliding background,
106    if isequal(floor(nbaver/2),nbaver)
107        nbaver=nbaver+1;%put the number of burst to an odd number (so the middle burst is defined)
108    end
109    nbaver_init=nbaver*siz(1);%propose by default an integer number of bursts
110end
111
112filebase=fullfile(Series.RootPath,Series.RootFile);
113dir_images=Series.RootPath;
114nom_type=Series.NomType;
115
116%create dir of the new images
117[dir_images,namebase]=fileparts(filebase);
118[path,subdir_ima]=fileparts(dir_images);
119curdir=pwd;
120cd(path);
121mkdir([subdir_ima '_b']);
122[xx,msg2] = fileattrib([subdir_ima '_b'],'+w','g'); %yield writing access (+w) to user group (g)
123if ~strcmp(msg2,'')
124    msgbox_uvmat('ERROR',['pb of permission for ' subdir_ima '_b' ': ' msg2])%error message for directory creation
125    cd(curdir)
126    return
127end
128cd(curdir);
129filebase_b=fullfile(path,[subdir_ima '_b'],namebase);
130
131prompt = {'Number of images for the sliding background';'The number of positions (laser slices)';'volume scan mode (Yes/No)';...
132        'the luminosity rank chosen to define the background (0.1=for dense particle seeding, 0.5 (median) for sparse particles'};
133dlg_title = ['get (slice by slice) a sliding background and substract to each image, result in subdir ' subdir_ima '_b'];
134num_lines= 3;
135def     = { num2str(nbaver_init);num2str(nbslice_i);'No';'0.1'};
136answer = inputdlg(prompt,dlg_title,num_lines,def);
137set(hseries.ParamVal,'String',answer([1 [3:4]]))
138set(hseries.ParamVal,'Visible','on')
139
140nbaver_ima=str2num(answer{1});%number of images for the sliding background
141nbaver=ceil(nbaver_ima/siz(1))%number of bursts for the sliding background
142if isequal(floor(nbaver/2),nbaver)
143   nbaver=nbaver+1%put the number of burst to an odd number (so the middle burst is defined)
144end
145% if isequal(nbaver,round(nbaver))
146step=siz(1);%case of bursts: the sliding background is shifted by one burst
147% else
148%    step=1;
149% end
150vol_test=answer{3};
151if isequal(vol_test,'Yes')
152    nbfield2=1;%case of volume: no consecutive series at a given level
153    nbslice_i=siz(1);%number of slices
154else
155    nbfield2=siz(1); %nb of consecutive images at each level(burst)
156    if siz(2)>1
157       nbslice_i=str2num(answer{2})/(num_i1(1,2)-num_i1(1,1));% number of slices
158    else
159        nbslice_i=1;
160    end
161    if ~isequal(floor(nbslice_i),nbslice_i)
162        msgbox_uvmat('ERROR','the number of slices must be a multiple of the i increment')
163        return
164    end
165end
166% nbaver=floor(nbaver_init/nbfield2); % number of bursts used for the sliding background,
167% if isequal(floor(nbaver/2),nbaver)
168%     nbaver=nbaver+1;%put the number of burst to an odd number
169% end
170%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
171%nbaver_ima=nbaver*nbfield2;% adjust the number of sliding images A  REMETRE
172% if ~isequal(nbaver_ima,nbaver_init)
173%     hwarn=warndlg(['number of images in the sliding average adjusted to ' num2str(nbaver_ima)]);
174%     set(hwarn,'Units','normalized')
175%     set(hwarn,'Position',[0.3 0.3 0.4 0.1])
176% end
177%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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
185% test_plot=isequal(answer{5},'Yes'); %=1 to display background images
186if nbaver_ima > nbfield*nbfield2
187    errordlg('number of images in a slice smaller than the proposed number of images for the sliding average')
188    return
189end
190nbfirst=(ceil(nbaver/2))*step;
191if nbfirst>nbaver_ima
192    nbfirst=ceil(nbaver_ima/2)
193    step=1;
194    nbaver=nbaver_ima;
195end
196
197%copy the xml file
198if exist([filebase '.xml'],'file')
199    copyfile([filebase '.xml'],[filebase_b '.xml']);% copy the .civ file
200    t=xmltree([filebase_b '.xml']);
201   
202    %update information on the first image name in the series
203    uid_Heading=find(t,'ImaDoc/Heading');
204    if isempty(uid_Heading)
205        [t,uid_Heading]=add(t,1,'element','Heading');
206    end   
207    uid_ImageName=find(t,'ImaDoc/Heading/ImageName');
208    ImageName=name_generator(filebase_b,num_i1(1),num_j1(1),'.png',Series.NomType);
209    [pth,ImageName]=fileparts(ImageName);
210    ImageName=[ImageName '.png']
211    if isempty(uid_ImageName)
212       [t,uid_ImageName]=add(t,uid_Heading,'element','ImageName');
213    end
214    uid_value=children(t,uid_ImageName);
215    if isempty(uid_value)
216        t=add(t,uid_ImageName,'chardata',ImageName)%indicate  name of the first image, with ;png extension
217    else
218        t=set(t,uid_value(1),'value',ImageName)%indicate  name of the first image, with ;png extension
219    end 
220
221    %add information about image transform
222    [t,new_uid]=add(t,1,'element','ImageTransform');
223    [t,NameFunction_uid]=add(t,new_uid,'element','NameFunction');
224    [t]=add(t,NameFunction_uid,'chardata','sub_background');     
225    [t,NbSlice_uid]=add(t,new_uid,'element','NbSlice');
226    [t]=add(t,new_uid,'chardata',num2str(nbslice_i));
227    [t,NbSlidingImages_uid]=add(t,new_uid,'element','NbSlidingImages');
228    [t]=add(t,NbSlidingImages_uid,'chardata',num2str(nbaver));
229    [t,LuminosityRank_uid]=add(t,new_uid,'element','RankBackground');
230    [t]=add(t,LuminosityRank_uid,'chardata',num2str(rank));% luminosity rank almong the nbaver sliding images
231    save(t,[filebase_b '.xml'])
232elseif exist([filebase '.civ'],'file')
233    copyfile([filebase '.civ'],[filebase_b '.civ']);% copy the .civ file
234end
235%copy the mask
236if exist([filebase '_1mask_1'],'file')
237     copyfile([filebase '_1mask_1'],[filebase_b '_1mask_1']);% copy the mask file
238end
239
240%MAIN LOOP ON SLICES
241for islice=1:nbslice_i
242    %select the series of image indices at the level islice
243    for ifield=1:nbfield
244        for iburst=1:nbfield2
245            indselect(iburst,ifield)=((ifield-1)*nbslice_i+(islice-1))*nbfield2+iburst;
246        end
247    end 
248    %read the first series of nbaver_ima images and sort by luminosity at each pixel
249    for ifield = 1:nbaver_ima
250        ifile=indselect(ifield);
251        filename=name_generator(filebase,num_i1(ifile),num_j1(ifile),Series.FileExt,Series.NomType);
252        Aread=read_image(filename,FileType,num_i1(ifile),MovieObject);
253        Ak(:,:,ifield)=Aread;           
254    end
255    Asort=sort(Ak,3);%sort the luminosity of images at each point
256    B=Asort(:,:,rank);%background image
257%
258%     namemean=name_generator([filebase '_back'],islice,[],'.png','_i');% makes the file name
259%     imwrite(B,namemean,'BitDepth',16); % save the first background image
260%     ['background image for slice ' num2str(islice) ' saved in ' namemean]
261    %substract the background from each of the first images and save the new images
262%     for ifield=1:floor(nbaver_ima/2)+1
263    'first background image will be substracted'
264
265    for ifield=1:nbfirst
266            Acor=double(Ak(:,:,ifield))-double(B);%substract background to the current image
267            Acor=(Acor>0).*Acor; % put to 0 the negative elements in Acor
268            C=uint16(Acor);% set to integer 16 bits
269            ifile=indselect(ifield);
270            newname=name_generator(filebase_b,num_i1(ifile),num_j1(ifile),'.png',nom_type)% makes the new file name
271            imwrite(C,newname,'BitDepth',16); % save the new image
272    end
273    %repeat the operation on a sliding series of nbaver*nbfield2 images
274    'sliding background image will be substracted'
275    if nbfield_slice > nbaver_ima
276%         for ifield = floor(nbaver_ima/2)+2:step:nbfield_slice-floor(nbaver_ima/2)
277        for ifield = step*ceil(nbaver/2)+1:step:nbfield_slice-step*floor(nbaver/2)
278            stopstate=get(hseries.RUN,'BusyAction');
279            if isequal(stopstate,'queue')% enable STOP command
280                update_waitbar(hseries.waitbar,WaitbarPos,(ifield+(islice-1)*nbfield_slice)/(nbfield_slice*nbslice_i))
281                (ifield+(islice-1)*nbfield_slice)/(nbfield_slice*nbslice_i)
282                Ak(:,:,[1:nbaver_ima-step])=Ak(:,:,[1+step:nbaver_ima]);% shift the current image series by one burst (step)
283                %incorporate next burst in the current image series
284                for iburst=1:step
285                    ifile=indselect(ifield+step*floor(nbaver/2)+iburst-1);
286                    filename=name_generator(filebase,num_i1(ifile),num_j1(ifile),Series.FileExt,Series.NomType);
287                    Aread=read_image(filename,FileType,num_i1(ifile),MovieObject);
288                    Ak(:,:,nbaver_ima-step+iburst)=Aread;
289                end
290                Asort=sort(Ak,3);%sort the new current image series by luminosity
291                B=Asort(:,:,rank);%current background image
292                for iburst=1:step
293%                     Acor=double(Ak(:,:,floor(nbaver_ima/2)+iburst-1))-double(B);
294                    index=step*floor(nbaver/2)+iburst;
295                    Acor=double(Ak(:,:,index))-double(B);
296                    Acor=(Acor>0).*Acor; % put to 0 the negative elements in Acor
297                    C=uint16(Acor);
298                    ifile=indselect(ifield+iburst-1);
299                    [newname]=...
300                       name_generator(filebase_b,num_i1(ifile),num_j1(ifile),'.png',Series.NomType)
301%                     newname=name_generator(filebase_b,num_i1(indselect(ifield+iburst-1)),num_j1(indselect(ifield+iburst-1)),'.png',nom_type)% makes the new file name
302                    imwrite(C,newname,'BitDepth',16); % save the new image
303                end 
304            else
305                return
306            end
307        end
308    end
309
310%substract the background from the last images
311%     for ifield=nbfield_slice-floor(nbaver_ima/2)+1:nbfield_slice
312     'last background image will be substracted'
313     ifield=nbfield_slice-(step*ceil(nbaver/2))+1:nbfield_slice;
314    for ifield=nbfield_slice-(step*floor(nbaver/2))+1:nbfield_slice 
315        index=ifield-nbfield_slice+step*(2*floor(nbaver/2)+1);
316        Acor=double(Ak(:,:,index))-double(B);
317        Acor=(Acor>0).*Acor; % put to 0 the negative elements in Acor
318        C=uint16(Acor);
319        ifile=indselect(ifield);
320        newname=name_generator(filebase_b,num_i1(ifile),num_j1(ifile),'.png',nom_type)% makes the new file name
321        imwrite(C,newname,'BitDepth',16); % save the new image
322    end
323end
324
325%finish the waitbar
326update_waitbar(hseries.waitbar,WaitbarPos,1)
327
328
329%------------------------------------------------------------------------
330%--read images and convert them to the uint16 format used for PIV
331function A=read_image(filename,type_ima,num,MovieObject)
332%------------------------------------------------------------------------
333%num is the view number needed for an avi movie
334switch type_ima
335    case 'movie'
336        A=read(MovieObject,num);
337    case 'avi'
338        mov=aviread(filename,num);
339        A=frame2im(mov(1));
340    case 'multimage'
341        A=imread(filename,num);
342    case 'image'   
343        A=imread(filename);
344end
345siz=size(A);
346if length(siz)==3;%color images
347    A=sum(double(A),3);
348end
349   
Note: See TracBrowser for help on using the repository browser.