source: trunk/src/series/particle_tracking.m @ 618

Last change on this file since 618 was 618, checked in by sommeria, 11 years ago

particle_tracking updated

File size: 21.1 KB
Line 
1% function ParamOut=particle_tracking(Param)
2%
3% Method:
4   
5% Organization of image indices:
6   
7%INPUT:
8% num_i1: matrix of image indices i
9% num_j1: matrix of image indices j, must be the same size as num_i1
10% num_i2 and num_j2: not used for a function acting on images
11% Series: matlab structure containing parameters, as defined by the interface UVMAT/series
12%       Series.RootPath{1}: path to the image series
13%       Series.RootFile{1}: root file name
14%       Series.FileExt{1}: image file extension
15%       Series.NomType{1}: nomenclature type for file in
16%
17% Method:
18%       Series.NbSlice: %number of slices defined on the interface
19% global A rangx0 rangy0 minA maxA; % make current image A accessible in workspace
20% global hfig1 hfig2 scalar
21% global Abackg nbpart lum diam
22%%%%%%%%%%%%%%ù
23%
24%%%%%%%%%%% GENERAL TO ALL SERIES ACTION FCTS %%%%%%%%%%%%%%%%%%%%%%%%%%%
25%
26%OUTPUT
27% ParamOut: sets options in the GUI series.fig needed for the function
28%
29%INPUT:
30% In run mode, the input parameters are given as a Matlab structure Param copied from the GUI series.
31% In batch mode, Param is the name of the corresponding xml file containing the same information
32% when Param.Action.RUN=0 (as activated when the current Action is selected
33% in series), the function ouput paramOut set the activation of the needed GUI elements
34%
35% Param contains the elements:(use the menu bar command 'export/GUI config' in series to
36% 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: directory extension for data outputs
41%    .Action: .ActionName: name of the current activated function
42%             .ActionPath:   path of the current activated function
43%             .ActionExt: fct extension ('.m', Matlab fct, '.sh', compiled   Matlab fct
44%             .RUN =0 for GUI input, =1 for function activation
45%             .RunMode='local','background', 'cluster': type of function  use
46%             
47%    .IndexRange: set the file or frame indices on which the action must be performed
48%    .FieldTransform: .TransformName: name of the selected transform function
49%                     .TransformPath:   path  of the selected transform function
50%    .InputFields: sub structure describing the input fields withfields
51%              .FieldName: name(s) of the field
52%              .VelType: velocity type
53%              .FieldName_1: name of the second field in case of two input series
54%              .VelType_1: velocity type of the second field in case of two input series
55%              .Coord_y: name of y coordinate variable
56%              .Coord_x: name of x coordinate variable
57%    .ProjObject: %sub structure describing a projection object (read from ancillary GUI set_object)
58%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
59
60function ParamOut=particle_tracking(Param)
61
62%% set the input elements needed on the GUI series when the action is selected in the menu ActionName
63if isstruct(Param) && isequal(Param.Action.RUN,0)
64    % general settings of the GUI:
65    ParamOut.AllowInputSort='off';% allow alphabetic sorting of the list of input file SubDir (options 'off'/'on', 'off' by default)
66    ParamOut.WholeIndexRange='off';% prescribes the file index ranges from min to max (options 'off'/'on', 'off' by default)
67    ParamOut.NbSlice='off'; %nbre of slices ('off' by default)
68    ParamOut.VelType='off';% menu for selecting the velocity type (options 'off'/'one'/'two',  'off' by default)
69    ParamOut.FieldName='off';% menu for selecting the field (s) in the input file(options 'off'/'one'/'two', 'off' by default)
70    ParamOut.FieldTransform = 'off';%can use a transform function
71    ParamOut.ProjObject='off';%can use projection object(option 'off'/'on',
72    ParamOut.Mask='off';%can use mask option   (option 'off'/'on', 'off' by default)
73    ParamOut.OutputDirExt='.track';%set the output dir extension
74    ParamOut.OutputFileMode='NbSlice';% '=NbInput': 1 output file per input file index, '=NbInput_i': 1 file per input file index i, '=NbSlice': 1 file per slice
75    filecell=get_file_series(Param);%check existence of the first input file
76    if ~exist(filecell{1,1},'file')
77        msgbox_uvmat('WARNING','the first input file does not exist')
78    end
79    % parameters specific to the function 'particle_tracking'
80    %Par.Nblock=25;%size of image subblocks for analysis
81    Par.Nblock=[];%no subblock for background determination
82    % Numexp=inputdlg('Entrer le numero','numexp',1);
83    % numexp=str2num(Numexp{1})
84    % Thresh=inputdlg('Entrer le seuil de luminosite (4000)','thresh',1);
85    % thresh=str2num(Thresh{1});%threshold for detection of particle luminosity weight
86    %filter particle detection
87    Par.ThreshLum=-2000;% luminosity threshold for particle detection, < 0 for black particles, >0 for white particles
88    ParamOut.ActionInput=Par;
89    return
90end
91
92%%%%%%%%%%%%  STANDARD RUN PART  %%%%%%%%%%%%
93ParamOut=[];
94%% read input parameters from an xml file if input is a file name (batch mode)
95checkrun=1;
96if ischar(Param)
97    Param=xml2struct(Param);% read Param as input file (batch case)
98    checkrun=0;
99end
100
101%% define the directory for result file
102OutputDir=[Param.OutputSubDir Param.OutputDirExt];
103
104%% root input file(s) name, type and index series
105RootPath=Param.InputTable{1,1};
106RootFile=Param.InputTable{1,3};
107SubDir=Param.InputTable{1,2};
108NomType=Param.InputTable{1,4};
109FileExt=Param.InputTable{1,5};
110[filecell,i1_series,i2_series,j1_series,j2_series]=get_file_series(Param);
111%%%%%%%%%%%%
112% The cell array filecell is the list of input file names, while
113% filecell{iview,fileindex}:
114%        iview: line in the table corresponding to a given file series
115%        fileindex: file index within  the file series,
116% 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
117% i1_series(iview,fileindex) expresses the same indices as a 1D array in file indices
118%%%%%%%%%%%%
119nbview=numel(i1_series);%number of input file series (lines in InputTable)
120nbfield_j=size(i1_series{1},1); %nb of fields for the j index (bursts or volume slices)
121nbfield_i=size(i1_series{1},2); %nb of fields for the i index
122nbfield=nbfield_j*nbfield_i; %total number of fields
123[first_i,tild,last_i,first_j,tild,last_j,errormsg]=get_index_range(Param.IndexRange);
124if ~isempty(errormsg),display(errormsg),return,end
125
126%% frame index for movie or multimage file input 
127if ~isempty(j1_series{1})
128    frame_index=j1_series{1};
129else
130    frame_index=i1_series{1};
131end
132
133%% check the input file type 
134[FileType,FileInfo,VideoObject]=get_file_type(filecell{1,1});
135ImageTypeOptions={'image','multimage','mmreader','video'};
136if isempty(find(strcmp(FileType,ImageTypeOptions)))
137    disp('input file not images')
138    return
139end
140
141%% calibration data and timing: read the ImaDoc files
142[XmlData,NbSlice_calib,time,errormsg]=read_multimadoc(RootPath,SubDir,RootFile,FileExt,i1_series,i2_series,j1_series,j2_series);
143
144%%%%%%%%%%%%   SPECIFIC PART (to edit) %%%%%%%%%%%%
145%filter for particle center of mass(luminosity)
146Nblock=Param.ActionInput.Nblock;
147ThreshLum=Param.ActionInput.ThreshLum;% luminosity threshold for particle detection, < 0 for black particles, >0 for white particles
148
149hh=ones(5,5);
150hh(1,1)=0;
151hh(1,5)=0;% sum luminosity on the 5x5 domain without corners
152hh(5,1)=0;
153hh(5,5)=0;
154hdx=[-2:1:2];
155hdy=[-2:1:2];
156[hdX,hdY]=meshgrid(hdx,hdy);
157hdX(1,1)=0;
158hdX(1,5)=0;% sum luminosity on the 5x5 domain -corners
159hdX(5,1)=0;
160hdX(5,5)=0;
161hdY(1,1)=0;
162hdY(1,5)=0;% sum luminosity on the 5x5 domain -corners
163hdY(5,1)=0;
164hdY(5,5)=0;
165%Parameters for image plotting
166pxcm=1;
167pycm=1;%scaling
168PlotParam.AutoScal=0;
169PlotParam.MaxA=700;
170PlotParam.MinA=0;
171PlotParam.BW=1;
172PlotParam.Contours=0;
173
174%%%%%%%%%%%%%%
175% sizfiles=size(num_j1)
176% nbfield=sizfiles(1); %number of images in a burst
177% %%%%%%%%%
178% hRUN=findobj(Series.hseries,'Tag','RUN');
179% hwaitbar=findobj(Series.hseries,'Tag','waitbar');%handles of the waitbar
180% waitbarpos(1)=Series.WaitbarPos(1);%x position of the waitbar
181% waitbarpos(3)=Series.WaitbarPos(3);% width of the waitbar
182% filebase=fullfile(Series.RootPath{1},Series.RootFile{1});
183% dir_images=Series.RootPath{1};
184% nom_type=Series.NomType{1};
185% [error,Heading,nom_type_read,ext_ima_read,tt,TimeUnit,mode,NbSlice,npx,npy,Calib]=read_imadoc([filebase '.xml']);
186
187%%  mask to reduce the  working area (optional)
188CheckMask=0;
189if isfield(Param,'CheckMask') && isequal(Param.CheckMask,1)
190    [maskname,TestMask]=name_generator([filebase '_1mask'],1,1,'.png','_i');
191        MaskIma=imread(maskname);
192        Mask=MaskIma>=200;%=1 for good points, 0 for bad
193    CheckMask=1;
194end
195%
196% %create dir of the new images
197% [dir_images,namebase]=fileparts(filebase);
198% [path,subdir_ima]=fileparts(dir_images)
199% curdir=pwd;
200% cd(path);
201% mkdir([subdir_ima '_b']);
202% cd(curdir);
203% filebase_b=fullfile(path,[subdir_ima '_b'],namebase);
204
205%
206%
207% lengthtot=siz(1)*siz(2);
208% nbfield=floor(lengthtot/(nbfield2*nbslice_i));%total number of i indexes (adjusted to an integer number of slices)
209% nbfield_slice=nbfield*nbfield2;% number of fields per slice
210% % test_plot=isequal(answer{5},'Yes'); %=1 to display background images
211% if nbaver_ima > nbfield*nbfield2
212%     errordlg('number of images in a slice smaller than the proposed number of images for the sliding average')
213%     return
214
215for ifile=1:nbfield
216    if checkrun
217        if strcmp(get(Param.RUNHandle,'BusyAction'),'queue')
218            update_waitbar(Param.WaitbarHandle,ifile/nbfield)
219        else
220            break% leave the loop if the STOP button is activated on the GUI series
221        end
222    end
223    if ~isempty(j1_series)&&~isequal(j1_series,{[]})
224        j1=j1_series{1}(ifile);
225    end
226    filename=fullfile_uvmat(RootPath,SubDir,RootFile,FileExt,NomType,i1_series{1}(ifile),[],j1);
227    A=read_image(filename,FileType,VideoObject,frame_index(ifile));
228    if ndims(A)==3;%color images
229        A=sum(double(A),3);% take the sum of color components
230    end
231    if ThreshLum<0
232        A=max(max(A))-A;%take the negative
233    end
234    if CheckMask
235        A=A.*Mask;
236    end
237    if isempty(Nblock)
238        A=A-min(min(A));%substract absolute mean
239    else
240        Aflagmin=sparse(imregionalmin(A));%Amin=1 for local image minima
241        Amin=A.*Aflagmin;%values of A at local minima
242        % local background: find all the local minima in image subblocks
243        sumblock= inline('sum(sum(x(:)))');
244        Backgi=blkproc(Amin,[Nblock Nblock],sumblock);% take the sum in  blocks
245        Bmin=blkproc(Aflagmin,[Nblock Nblock],sumblock);% find the number of minima in blocks
246        Backgi=Backgi./Bmin; % find the average of minima in blocks
247        % Backg=Backg+Backgi;
248        Backg=Backgi;
249        A=A-imresize(Backg/nburst(1),size(A),'bilinear');% interpolate to the initial size image and substract
250    end
251    Aflagmax=sparse(imregionalmax(A));%find local maxima
252    Plum=imfilter(A,hh);% sum A on 5x% domains
253    Plum=Aflagmax.*Plum;% Plum gives the particle luminosity at each particle location, 0 elsewhere
254    %make statistics on particles,restricted to a subdomain Sub
255    [Js,Is,lum]=find(Plum);%particle luminosity
256    Plum=(Plum>ThreshLum).*Plum;% introduce a threshold for particle luminosity
257    Aflagmax=Aflagmax.*(Plum>ThreshLum);
258    [Js,Is,lum]=find(Plum);%particle luminosity
259    nbtotal=size(Is)
260    nbtotal=nbtotal(1);
261    %particle size
262    Parea=Aflagmax.*(Plum./A); %particle luminosity/max luminosity=area
263    Pdiam=sqrt(Parea);
264    [Js,Is,diam]=find(Pdiam);%particle location
265   
266    %%%%%%%%%%%%%%%%%%%%%
267   
268    %nbre of particles per block
269%     nbpart=blkproc(Aflagmax,[Nblock Nblock],sumblock);%
270%     npb=size(nbpart);
271%     rangxb=[0.5 (npb(2)-0.5)]*Nblock; % pixel x coordinates for image display
272%     rangyb=[(npb(1)-0.5) 0.5]*Nblock; % pixel y coordinates for image display
273%     image(rangxb,rangyb,nbpart);
274   
275    % get the particle centre of mass
276    dx=imfilter(A,hdX);
277    dy=imfilter(A,-hdY);
278    dx=Aflagmax.*(dx./Plum);
279    dy=Aflagmax.*(dy./Plum);
280    dx=dx/pxcm;
281    dy=dy/pycm;
282    I=([1:npxy(2)]-0.5)/pxcm; %x pos
283    J=([npxy(1):-1:1]-0.5)/pycm; %y pos
284    [Ipos,Jpos]=meshgrid(I,J);
285    Ipos=reshape(Ipos,1,npxy(2)*npxy(1));
286    Jpos=reshape(Jpos,1,npxy(2)*npxy(1));
287    dx=reshape(dx,1,npxy(2)*npxy(1));
288    dy=reshape(dy,1,npxy(2)*npxy(1));
289    Aflag=reshape(Aflagmax,1,npxy(2)*npxy(1));
290    ind=find(Aflag);% select particle positions
291    XPart{ifile}=Ipos(ind)+dx(ind);
292    YPart{ifile}=Jpos(ind)+dy(ind);     
293end
294hold off
295
296size(XPart{1})
297
298%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
299%Trajectoires
300%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
301for ifile=1:nbfield
302   
303    [XPart{ifile},YPart{ifile}]=phys_XYZ(Calib,XPart{ifile},YPart{ifile});
304
305end
306
307if nbfield>2
308    figpart=figure
309    hold on
310    plot(XPart{1}(:),YPart{1}(:),'r+')
311    plot(XPart{2}(:),YPart{2}(:),'b+')
312    plot(XPart{3}(:),YPart{3}(:),'y+')
313    legend('particules image 1','particules image 2','particules image 3');
314    xlabel('x (cm)');
315    ylabel('y (cm)');
316    title('Position des particules')
317else
318   figpart=figure
319    hold on
320    plot(XPart{1}(:),YPart{1}(:),'r+')
321    plot(XPart{2}(:),YPart{2}(:),'b+')
322    legend('particules image 1','particules image 2');
323    xlabel('x (cm)');
324    ylabel('y (cm)');
325    title('Position des particules')
326end   
327
328%     prompt={'Ymin (cm)','Ymax( cm)','Xmin (cm)','Xmax (cm)'};
329%     Rep=inputdlg(prompt,'Experiment');
330%     Ymin=str2double(Rep(1));
331%     Ymax=str2double(Rep(2));
332%     Xmin=str2double(Rep(3));
333%     Xmax=str2double(Rep(4));
334   
335    Ymin=6;
336    Ymax=14;
337    Xmin=15;
338    Xmax=35;
339   
340    plot(Xmin,Ymin,'g+')
341    plot(Xmin,Ymax,'g+')
342    plot(Xmax,Ymin,'g+')
343    plot(Xmax,Ymax,'g+')
344
345   
346 for ima=2:nbfield   
347    t{1}=0*ones(size(XPart{1},2),1);
348    burst(1)=0;
349    burst(2)=0.018;
350    burst(3)=0.036;
351%     nburst=strcat('burst',num2str(ima-1),'-',num2str(ima),' (s)');
352%     prompt={'burst (s)'};
353%     Rep=inputdlg(prompt,nburst);
354%     burst(ima)=str2double(Rep(1));
355    t{ima}=(burst(ima)+burst(ima-1))*ones(size(XPart{ima},2),1);
356 end
357
358
359 
360 for ima=1:nbfield
361
362    IndY{ima}=find(YPart{ima}>Ymin & YPart{ima}<Ymax & XPart{ima}>Xmin & XPart{ima}<Xmax);
363    XPart{ima}=XPart{ima}(IndY{ima});
364    YPart{ima}=YPart{ima}(IndY{ima});
365   
366       
367end
368
369
370
371%%%%%%%%%%%%%%%%%%%%%%%
372% Calcul de v1
373%%%%%%%%%%%%%%%%%%%%%%%
374
375for i=1:size(XPart{1},2)
376    MatPos{1}(i,1)=XPart{1}(i);
377    MatPos{1}(i,2)=YPart{1}(i);
378    MatPos{1}(i,3)=t{1}(i);
379    %MatPos{1}(i,4)=i;
380end
381
382for j=1:size(XPart{2},2)-1
383    MatPos{1}(j+size(XPart{1},2),1)=XPart{2}(j);
384    MatPos{1}(j+size(XPart{1},2),2)=YPart{2}(j);
385    MatPos{1}(j+size(XPart{1},2),3)=t{2}(j);
386    %MatPos{1}(j,4)=j+size(XPart{1},2);
387end
388 
389% Dmax=inputdlg('Entrer la distance maximum (0.25 cm)','dmax (cm)',1)
390% dmax=str2num(Dmax{1});
391dmax=0.23;
392
393result{1}=track(MatPos{1},dmax);
394
395izero=1;
396for itest=1:1:size(result{1},1)-1
397    if  result{1}(itest+1,4)==result{1}(itest,4)
398        vitu{1}(izero,1)=(result{1}(itest+1,1)-result{1}(itest,1))/burst(2);
399        vitu{1}(izero,2)=result{1}(itest,4);
400        vitv{1}(izero,1)=(result{1}(itest+1,2)-result{1}(itest,2))/burst(2);
401        vitv{1}(izero,2)=result{1}(itest,4);
402        MatPos{2}(izero,1)=result{1}(itest,1);
403        MatPos{2}(izero,2)=result{1}(itest,2);
404        izero=izero+1;
405    end
406end
407
408
409vitfu{1}=vitu{1};
410vitfv{1}=vitv{1};
411
412
413%%%%%%%%%%%%%%%%%%%%%%%
414% Calcul de vi
415%%%%%%%%%%%%%%%%%%%%%%%
416
417
418if nbfield>2
419    for ima=2:nbfield-1
420       
421       for i=1:size(MatPos{ima},1)
422        MatPos{ima+1}(i,1)=MatPos{ima}(i,1)+(burst(ima+1)*vitfu{ima-1}(i));
423        MatPos{ima+1}(i,2)=MatPos{ima}(i,2)+(burst(ima+1)*vitfv{ima-1}(i));
424        MatPos{ima+1}(i,3)=t{ima}(i);
425      end
426
427      for j=1:size(XPart{ima+1},2)-1
428          MatPos{ima+1}(j+size(MatPos{ima},1),1)=XPart{ima+1}(j);
429          MatPos{ima+1}(j+size(MatPos{ima},1),2)=YPart{ima+1}(j);
430          MatPos{ima+1}(j+size(MatPos{ima},1),3)=t{ima+1}(j);
431      end
432       
433     
434    result{ima}=track(MatPos{ima+1},0.15);
435       
436        izero=1;
437        for itest=1:1:size(result{ima},1)-1
438            if  result{ima}(itest+1,4)==result{ima}(itest,4)
439                vitu{ima}(izero,1)=(result{ima}(itest+1,1)-result{ima}(itest,1))/burst(ima+1);
440                vitu{ima}(izero,2)=result{ima}(itest,4);
441                vitv{ima}(izero,1)=(result{ima}(itest+1,2)-result{ima}(itest,2))/burst(ima+1);
442                vitv{ima}(izero,2)=result{ima}(itest,4);
443                MatPos{ima+2}(izero,1)=result{ima}(itest,1);
444                MatPos{ima+2}(izero,2)=result{ima}(itest,2);
445                izero=izero+1;
446            end   
447        end
448
449            i=vitu{ima}(1,2):1:vitu{ima}(end,2)
450           
451              vitfu{ima}(:,1)=vitfu{ima-1}(i,1)+vitu{ima}(:,1);
452              vitfv{ima}(:,1)=vitfv{ima-1}(i,1)+vitv{ima}(:,1);
453              vitfu{ima}(:,2)=vitu{ima}(:,2);
454              vitfv{ima}(:,2)=vitv{ima}(:,2);
455
456            vitfu{ima-1}=vitfu{ima-1}(i,1);
457            vitfu{ima-1}(:,2)=i;
458            vitfv{ima-1}=vitfv{ima-1}(i,1);
459            vitfv{ima-1}(:,2)=i;
460            i=1:1:size(vitfu{ima-1},1)
461            xpos=MatPos{2}(i,1)
462            ypos=MatPos{2}(i,2)
463      end
464    end
465
466
467
468    figure
469    hold on
470    plot(MatPos{1}(:,1),MatPos{1}(:,2),'r+')
471    plot(MatPos{2}(:,1),MatPos{2}(:,2),'b+')
472    plot(MatPos{4}(:,1),MatPos{4}(:,2),'y+')
473    quiver(xpos(:),ypos(:),vitfu{1}(:,1),vitfv{1}(:,1),'g')
474    quiver(MatPos{4}(:,1),MatPos{4}(:,2),vitfu{2}(:,1),vitfv{2}(:,1),'k')
475    legend('particules image 1','particules image 2', 'particules image 3','vitesse 1-2 (cm/s)','vitesse 2-3 (cm/s)');
476    xlabel('x (cm)');
477    ylabel('y (cm)');
478    title('Position et vitesse (cm/s) des particules')
479   
480
481    for i=1:size(vitfu{end},1)
482     vitfuadd(i)=0;
483     vitfvadd(i)=0;
484    end
485
486   
487   
488         for i=1:1:size(vitfu{end}(:,1))
489           
490                for j=1:nbfield-1
491                    vitfuadd(i)= vitfuadd(i)+vitfu{j}(i,1);
492                    vitfvadd(i)= vitfvadd(i)+vitfv{j}(i,1);
493                    xpos1(i)=MatPos{1}(i,1);
494                    ypos1(i)=MatPos{1}(i,2);
495                    xpos2(i)=MatPos{2}(i,1);
496                    ypos2(i)=MatPos{2}(i,2);
497                   
498                end
499            end
500            sizexpos1=size(xpos1)
501
502    vitfumoy=vitfuadd./(nbfield-1)
503    vitfvmoy=vitfvadd./(nbfield-1)
504
505    testresult1=result{1}
506    testresult2=result{2}
507   
508if nbfield>2   
509    figure
510    hold on
511    plot(MatPos{1}(:,1),MatPos{1}(:,2),'r+')
512    plot(MatPos{2}(:,1),MatPos{2}(:,2),'b+')
513    quiver(xpos2(:),ypos2(:),vitfumoy(:),vitfvmoy(:),'g')
514    legend('particules image 1','particules image 2', 'vitesse moyenne (cm/s)');
515    xlabel('x (cm)');
516    ylabel('y (cm)');
517    title('Position et vitesse (cm/s) des particules')
518   
519else
520
521    figure
522    hold on
523    plot(MatPos{1}(:,1),MatPos{1}(:,2),'r+')
524    plot(MatPos{2}(:,1),MatPos{2}(:,2),'b+')
525    quiver(MatPos{2}(:,1),MatPos{2}(:,2),vitfu{1}(:),vitfv{1}(:),'g')
526    legend('particules image 1','particules image 2','vitesse 1-2 (cm/s)');
527    xlabel('x (cm)');
528    ylabel('y (cm)');
529    title('Position et vitesse (cm/s) des particules')
530   
531    vitfumoy=vitfu{1};
532    vitfvmoy=vitfv{1};
533
534end
535
536VitData.NbDim=2;
537VitData.NbCoord=2;
538VitData.CoordType='phys';
539VitData.dt=0.0185;
540VitData.CoordUnit='cm';
541VitData.Z=0;
542VitData.ListDimName={'nb_vectors'};
543VitData.DimValue=size(vitfumoy,2);
544VitData.ListVarName={'X'  'Y'  'U'  'V'  'F'};
545VitData.VarDimIndex={[1]  [1]  [1]  [1]  [1]};
546VitData.ListVarAttribute={'Role'};
547VitData.Role={'coord_x'  'coord_y'  'vector_x'  'vector_y'  'warnflag'};
548
549if nbfield>2
550    VitData.X=size(MatPos{4},1);
551    VitData.Y=size(MatPos{4},2);
552else
553    VitData.X=size(MatPos{2},1);
554    VitData.Y=size(MatPos{2},2);
555end
556
557VitData.U=size(vitfumoy,2);
558VitData.V=size(vitfvmoy,2);
559VitData.Style='plane';
560VitData.Time=[198.5203 198.5203];
561VitData.Action=Series.Action;
562
563if nbfield>2
564    VitData.X=MatPos{4}(:,1)';
565    VitData.Y=MatPos{4}(:,2)';
566else
567    VitData.X=MatPos{2}(:,1)';
568    VitData.Y=MatPos{2}(:,2)';
569end
570
571VitData.U=vitfumoy(:)';
572VitData.V=vitfvmoy(:)';
573
574if length(VitData.ListVarName) >= 4 & isequal(VitData.ListVarName(1:4), {'X'  'Y'  'U'  'V'})
575       VitData.ListAttribute={'nb_coord','nb_dim','dt','pixcmx','pixcmy','hart','civ','fix'};
576       VitData.nb_coord=2;
577       VitData.nb_dim=2;
578       VitData.dt=0.018;
579       VitData.absolut_time_T0=0;
580       VitData.pixcmx=1; %pix per cm (1 by default)
581       VitData.pixcmy=1; %pix per cm (1 by default)
582       VitData.hart=0;
583           if isequal(VitData.CoordType,'px')
584             VitData.civ=1;
585           else
586             VitData.civ=0;
587           end
588        VitData.fix=0;
589        VitData.ListVarName(1:4)={'vec_X'  'vec_Y'  'vec_U'  'vec_V'};
590        VitData.vec_X=VitData.X;
591        VitData.vec_Y=VitData.Y;
592        VitData.vec_U=VitData.U;
593        VitData.vec_V=VitData.V;
594end
595currentdir=pwd;%store the current working directory
596[Path_ima,Name]=fileparts(filebase);%Path of the image files (.civ)
597cd(Path_ima);%move to the directory of the images: needed to create the result dir by 'mkdir'
598dircur=pwd; %current working directory
599[m1,m2,m3]=mkdir('TRACK_test')
600cd(currentdir)
601[filename_nc,idetect]=name_generator(filebase,num_i1(1),num_j1(1),'.nc','_i_j1-j2',1,num_i1(1),num_j1(2),'TRACK_test')
602error=struct2nc(filename_nc,VitData); %save result file
603if isequal(error,0)
604    [filename_nc ' written']
605else
606    warndlg_uvmat(error,'ERROR')
607end
608
Note: See TracBrowser for help on using the repository browser.