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

Last change on this file since 140 was 128, checked in by sommeria, 14 years ago

series: give writting access to the group for all subdirectories produced
uvmat.fig: change of vect and scalar frames (to be consistent with view_field)
uvmat: various cleaning
plot_field: various cleaning to improve axes definition and avoid blinking
geometry_calib: improved dispay of point coordiantes, improved link with dataview for REPLICATE.
struct2nc: repair bug , file was not closed.
cell2tab: cleaning
dataview: improve the browser
civ: solve pb of image naming

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