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

Last change on this file since 977 was 977, checked in by g7moreau, 4 years ago
  • Update Copyright 2008-2017 notice
File size: 20.5 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
60%=======================================================================
61% Copyright 2008-2017, LEGI UMR 5519 / CNRS UGA G-INP, Grenoble, France
62%   http://www.legi.grenoble-inp.fr
63%   Joel.Sommeria - Joel.Sommeria (A) legi.cnrs.fr
64%
65%     This file is part of the toolbox UVMAT.
66%
67%     UVMAT is free software; you can redistribute it and/or modify
68%     it under the terms of the GNU General Public License as published
69%     by the Free Software Foundation; either version 2 of the license,
70%     or (at your option) any later version.
71%
72%     UVMAT is distributed in the hope that it will be useful,
73%     but WITHOUT ANY WARRANTY; without even the implied warranty of
74%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
75%     GNU General Public License (see LICENSE.txt) for more details.
76%=======================================================================
77
78function ParamOut=particle_tracking(Param)
79
80%% set the input elements needed on the GUI series when the action is selected in the menu ActionName
81if isstruct(Param) && isequal(Param.Action.RUN,0)
82    % general settings of the GUI:
83    ParamOut.AllowInputSort='off';% allow alphabetic sorting of the list of input file SubDir (options 'off'/'on', 'off' by default)
84    ParamOut.WholeIndexRange='off';% prescribes the file index ranges from min to max (options 'off'/'on', 'off' by default)
85    ParamOut.NbSlice='off'; %nbre of slices ('off' by default)
86    ParamOut.VelType='off';% menu for selecting the velocity type (options 'off'/'one'/'two',  'off' by default)
87    ParamOut.FieldName='off';% menu for selecting the field (s) in the input file(options 'off'/'one'/'two', 'off' by default)
88    ParamOut.FieldTransform = 'off';%can use a transform function
89    ParamOut.ProjObject='off';%can use projection object(option 'off'/'on',
90    ParamOut.Mask='off';%can use mask option   (option 'off'/'on', 'off' by default)
91    ParamOut.OutputDirExt='.track';%set the output dir extension
92    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
93    filecell=get_file_series(Param);%check existence of the first input file
94    if ~exist(filecell{1,1},'file')
95        msgbox_uvmat('WARNING','the first input file does not exist')
96    end
97    % parameters specific to the function 'particle_tracking'
98    Par.Nblock=10;%size of image subblocks for background determination, =[]: no sublock
99    Par.ThreshLum=70;% luminosity threshold for particle detection, < 0 for black particles, >0 for white particles
100  ParamOut.ActionInput=Par;
101    return
102end
103
104%%%%%%%%%%%%  STANDARD RUN PART  %%%%%%%%%%%%
105ParamOut=[];
106%% read input parameters from an xml file if input is a file name (batch mode)
107checkrun=1;
108if ischar(Param)
109    Param=xml2struct(Param);% read Param as input file (batch case)
110    checkrun=0;
111end
112hseries=findobj(allchild(0),'Tag','series');
113RUNHandle=findobj(hseries,'Tag','RUN');%handle of RUN button in GUI series
114WaitbarHandle=findobj(hseries,'Tag','Waitbar');%handle of waitbar in GUI series
115
116%% define the directory for result file
117OutputDir=[Param.OutputSubDir Param.OutputDirExt];
118
119%% root input file(s) name, type and index series
120RootPath=Param.InputTable{1,1};
121RootFile=Param.InputTable{1,3};
122SubDir=Param.InputTable{1,2};
123NomType=Param.InputTable{1,4};
124FileExt=Param.InputTable{1,5};
125[filecell,i1_series,i2_series,j1_series,j2_series]=get_file_series(Param);
126%%%%%%%%%%%%
127% The cell array filecell is the list of input file names, while
128% filecell{iview,fileindex}:
129%        iview: line in the table corresponding to a given file series
130%        fileindex: file index within  the file series,
131% 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
132% i1_series(iview,fileindex) expresses the same indices as a 1D array in file indices
133%%%%%%%%%%%%
134nbview=numel(i1_series);%number of input file series (lines in InputTable)
135nbfield_j=size(i1_series{1},1); %nb of fields for the j index (bursts or volume slices)
136nbfield_i=size(i1_series{1},2); %nb of fields for the i index
137nbfield=nbfield_j*nbfield_i; %total number of fields
138
139%% frame index for movie or multimage file input 
140if ~isempty(j1_series{1})
141    frame_index=j1_series{1};
142else
143    frame_index=i1_series{1};
144end
145
146%% check the input file type 
147[FileInfo,VideoObject]=get_file_info(filecell{1,1});
148FileType=FileInfo.FileType;
149ImageTypeOptions={'image','multimage','mmreader','video'};
150if isempty(find(strcmp(FileType,ImageTypeOptions)))
151    disp('input file not images')
152    return
153end
154
155%% calibration data and timing: read the ImaDoc files
156[XmlData,NbSlice_calib,time,errormsg]=read_multimadoc(RootPath,SubDir,RootFile,FileExt,i1_series,i2_series,j1_series,j2_series);
157
158%%%%%%%%%%%%   SPECIFIC PART (to edit) %%%%%%%%%%%%
159%filter for particle center of mass(luminosity)
160Nblock=Param.ActionInput.Nblock;
161ThreshLum=Param.ActionInput.ThreshLum;% luminosity threshold for particle detection, < 0 for black particles, >0 for white particles
162%AbsThreshold=30; %threshold below which a pixel is considered belonging to a float
163%
164hh=ones(5,5);
165hh(1,1)=0;
166hh(1,5)=0;% sum luminosity on the 5x5 domain without corners
167hh(5,1)=0;
168hh(5,5)=0;
169hdx=[-2:1:2];
170hdy=[-2:1:2];
171[hdX,hdY]=meshgrid(hdx,hdy);
172hdX(1,1)=0;
173hdX(1,5)=0;% sum luminosity on the 5x5 domain -corners
174hdX(5,1)=0;
175hdX(5,5)=0;
176hdY(1,1)=0;
177hdY(1,5)=0;% sum luminosity on the 5x5 domain -corners
178hdY(5,1)=0;
179hdY(5,5)=0;
180
181%%  mask to reduce the  working area (optional)
182CheckMask=0;
183if isfield(Param,'CheckMask') && isequal(Param.CheckMask,1)
184    [maskname,TestMask]=name_generator([filebase '_1mask'],1,1,'.png','_i');
185        MaskIma=imread(maskname);
186        Mask=MaskIma>=200;%=1 for good points, 0 for bad
187    CheckMask=1;
188end
189
190%%%%%% MAIN LOOP ON FRAMES %%%%%%
191for ifile=1:nbfield
192    if checkrun
193        update_waitbar(WaitbarHandle,ifile/nbfield)
194        if ~isempty(RUNHandle) &&ishandle(RUNHandle) && ~strcmp(get(RUNHandle,'BusyAction'),'queue')
195            disp('program stopped by user')
196            return
197        end
198    end
199    j1=1;
200    if ~isempty(j1_series)&&~isequal(j1_series,{[]})
201        j1=j1_series{1}(ifile);
202    end
203    filename=fullfile_uvmat(RootPath,SubDir,RootFile,FileExt,NomType,i1_series{1}(ifile),[],j1);
204    A=read_image(filename,FileType,VideoObject,frame_index(ifile));% read the current frame
205    if ndims(A)==3;%color images
206        A=sum(double(A),3);% take the sum of color components
207    end
208    if ThreshLum<0
209        A=max(max(A))-A;%take the negative
210    end
211    if CheckMask
212        A=A.*Mask;
213    end
214    if isempty(Nblock)
215        A=A-min(min(A));%substract absolute mean
216    else
217        Aflagmin=sparse(imregionalmin(A));%Amin=1 for local image minima
218        Amin=A.*Aflagmin;%values of A at local minima
219        % local background: find all the local minima in image subblocks
220        sumblock= inline('sum(sum(x(:)))');
221        Backgi=blkproc(Amin,[Nblock Nblock],sumblock);% take the sum in  blocks
222        Bmin=blkproc(Aflagmin,[Nblock Nblock],sumblock);% find the number of minima in blocks
223        Backgi=Backgi./Bmin; % find the average of minima in blocks
224        % Backg=Backg+Backgi;
225        Backg=Backgi;
226        A=A-imresize(Backg/nburst(1),size(A),'bilinear');% interpolate to the initial size image and substract
227    end
228    Aflagmax=sparse(imregionalmax(A));%find local maxima
229    Plum=imfilter(A,hh);% sum A on 5x% domains
230    Plum=Aflagmax.*Plum;% Plum gives the particle luminosity at each particle location, 0 elsewhere
231    %make statistics on particles,restricted to a subdomain Sub
232    [Js,Is,lum]=find(Plum);%particle luminosity
233    Plum=(Plum>ThreshLum).*Plum;% introduce a threshold for particle luminosity
234    Aflagmax=Aflagmax.*(Plum>ThreshLum);
235    [Js,Is,lum]=find(Plum);%particle luminosity
236    nbtotal=size(Is)
237    nbtotal=nbtotal(1);
238    %particle size
239    Parea=Aflagmax.*(Plum./A); %particle luminosity/max luminosity=area
240    Pdiam=sqrt(Parea);
241    [Js,Is,diam]=find(Pdiam);%particle location
242   
243    %%%%%%%%%%%%%%%%%%%%%
244   
245    %nbre of particles per block
246%     nbpart=blkproc(Aflagmax,[Nblock Nblock],sumblock);%
247%     npb=size(nbpart);
248%     rangxb=[0.5 (npb(2)-0.5)]*Nblock; % pixel x coordinates for image display
249%     rangyb=[(npb(1)-0.5) 0.5]*Nblock; % pixel y coordinates for image display
250%     image(rangxb,rangyb,nbpart);
251   
252    % get the particle centre of mass
253    dx=imfilter(A,hdX);
254    dy=imfilter(A,-hdY);
255    dx=Aflagmax.*(dx./Plum);
256    dy=Aflagmax.*(dy./Plum);
257    dx=dx/pxcm;
258    dy=dy/pycm;
259    I=([1:npxy(2)]-0.5)/pxcm; %x pos
260    J=([npxy(1):-1:1]-0.5)/pycm; %y pos
261    [Ipos,Jpos]=meshgrid(I,J);
262    Ipos=reshape(Ipos,1,npxy(2)*npxy(1));
263    Jpos=reshape(Jpos,1,npxy(2)*npxy(1));
264    dx=reshape(dx,1,npxy(2)*npxy(1));
265    dy=reshape(dy,1,npxy(2)*npxy(1));
266    Aflag=reshape(Aflagmax,1,npxy(2)*npxy(1));
267    ind=find(Aflag);% select particle positions
268    XPart{ifile}=Ipos(ind)+dx(ind);
269    YPart{ifile}=Jpos(ind)+dy(ind);     
270end
271hold off
272
273size(XPart{1})
274
275%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
276%Trajectoires
277%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
278for ifile=1:nbfield
279   
280    [XPart{ifile},YPart{ifile}]=phys_XYZ(Calib,XPart{ifile},YPart{ifile});
281
282end
283
284if nbfield>2
285    figpart=figure
286    hold on
287    plot(XPart{1}(:),YPart{1}(:),'r+')
288    plot(XPart{2}(:),YPart{2}(:),'b+')
289    plot(XPart{3}(:),YPart{3}(:),'y+')
290    legend('particules image 1','particules image 2','particules image 3');
291    xlabel('x (cm)');
292    ylabel('y (cm)');
293    title('Position des particules')
294else
295   figpart=figure
296    hold on
297    plot(XPart{1}(:),YPart{1}(:),'r+')
298    plot(XPart{2}(:),YPart{2}(:),'b+')
299    legend('particules image 1','particules image 2');
300    xlabel('x (cm)');
301    ylabel('y (cm)');
302    title('Position des particules')
303end   
304
305%     prompt={'Ymin (cm)','Ymax( cm)','Xmin (cm)','Xmax (cm)'};
306%     Rep=inputdlg(prompt,'Experiment');
307%     Ymin=str2double(Rep(1));
308%     Ymax=str2double(Rep(2));
309%     Xmin=str2double(Rep(3));
310%     Xmax=str2double(Rep(4));
311   
312    Ymin=6;
313    Ymax=14;
314    Xmin=15;
315    Xmax=35;
316   
317    plot(Xmin,Ymin,'g+')
318    plot(Xmin,Ymax,'g+')
319    plot(Xmax,Ymin,'g+')
320    plot(Xmax,Ymax,'g+')
321
322   
323 for ima=2:nbfield   
324    t{1}=0*ones(size(XPart{1},2),1);
325    burst(1)=0;
326    burst(2)=0.018;
327    burst(3)=0.036;
328%     nburst=strcat('burst',num2str(ima-1),'-',num2str(ima),' (s)');
329%     prompt={'burst (s)'};
330%     Rep=inputdlg(prompt,nburst);
331%     burst(ima)=str2double(Rep(1));
332    t{ima}=(burst(ima)+burst(ima-1))*ones(size(XPart{ima},2),1);
333 end
334
335
336 
337 for ima=1:nbfield
338
339    IndY{ima}=find(YPart{ima}>Ymin & YPart{ima}<Ymax & XPart{ima}>Xmin & XPart{ima}<Xmax);
340    XPart{ima}=XPart{ima}(IndY{ima});
341    YPart{ima}=YPart{ima}(IndY{ima});
342   
343       
344end
345
346
347
348%%%%%%%%%%%%%%%%%%%%%%%
349% Calcul de v1
350%%%%%%%%%%%%%%%%%%%%%%%
351
352for i=1:size(XPart{1},2)
353    MatPos{1}(i,1)=XPart{1}(i);
354    MatPos{1}(i,2)=YPart{1}(i);
355    MatPos{1}(i,3)=t{1}(i);
356    %MatPos{1}(i,4)=i;
357end
358
359for j=1:size(XPart{2},2)-1
360    MatPos{1}(j+size(XPart{1},2),1)=XPart{2}(j);
361    MatPos{1}(j+size(XPart{1},2),2)=YPart{2}(j);
362    MatPos{1}(j+size(XPart{1},2),3)=t{2}(j);
363    %MatPos{1}(j,4)=j+size(XPart{1},2);
364end
365 
366% Dmax=inputdlg('Entrer la distance maximum (0.25 cm)','dmax (cm)',1)
367% dmax=str2num(Dmax{1});
368dmax=0.23;
369
370result{1}=track(MatPos{1},dmax);
371
372izero=1;
373for itest=1:1:size(result{1},1)-1
374    if  result{1}(itest+1,4)==result{1}(itest,4)
375        vitu{1}(izero,1)=(result{1}(itest+1,1)-result{1}(itest,1))/burst(2);
376        vitu{1}(izero,2)=result{1}(itest,4);
377        vitv{1}(izero,1)=(result{1}(itest+1,2)-result{1}(itest,2))/burst(2);
378        vitv{1}(izero,2)=result{1}(itest,4);
379        MatPos{2}(izero,1)=result{1}(itest,1);
380        MatPos{2}(izero,2)=result{1}(itest,2);
381        izero=izero+1;
382    end
383end
384
385
386vitfu{1}=vitu{1};
387vitfv{1}=vitv{1};
388
389
390%%%%%%%%%%%%%%%%%%%%%%%
391% Calcul de vi
392%%%%%%%%%%%%%%%%%%%%%%%
393
394
395if nbfield>2
396    for ima=2:nbfield-1
397       
398       for i=1:size(MatPos{ima},1)
399        MatPos{ima+1}(i,1)=MatPos{ima}(i,1)+(burst(ima+1)*vitfu{ima-1}(i));
400        MatPos{ima+1}(i,2)=MatPos{ima}(i,2)+(burst(ima+1)*vitfv{ima-1}(i));
401        MatPos{ima+1}(i,3)=t{ima}(i);
402      end
403
404      for j=1:size(XPart{ima+1},2)-1
405          MatPos{ima+1}(j+size(MatPos{ima},1),1)=XPart{ima+1}(j);
406          MatPos{ima+1}(j+size(MatPos{ima},1),2)=YPart{ima+1}(j);
407          MatPos{ima+1}(j+size(MatPos{ima},1),3)=t{ima+1}(j);
408      end
409       
410     
411    result{ima}=track(MatPos{ima+1},0.15);
412       
413        izero=1;
414        for itest=1:1:size(result{ima},1)-1
415            if  result{ima}(itest+1,4)==result{ima}(itest,4)
416                vitu{ima}(izero,1)=(result{ima}(itest+1,1)-result{ima}(itest,1))/burst(ima+1);
417                vitu{ima}(izero,2)=result{ima}(itest,4);
418                vitv{ima}(izero,1)=(result{ima}(itest+1,2)-result{ima}(itest,2))/burst(ima+1);
419                vitv{ima}(izero,2)=result{ima}(itest,4);
420                MatPos{ima+2}(izero,1)=result{ima}(itest,1);
421                MatPos{ima+2}(izero,2)=result{ima}(itest,2);
422                izero=izero+1;
423            end   
424        end
425
426            i=vitu{ima}(1,2):1:vitu{ima}(end,2)
427           
428              vitfu{ima}(:,1)=vitfu{ima-1}(i,1)+vitu{ima}(:,1);
429              vitfv{ima}(:,1)=vitfv{ima-1}(i,1)+vitv{ima}(:,1);
430              vitfu{ima}(:,2)=vitu{ima}(:,2);
431              vitfv{ima}(:,2)=vitv{ima}(:,2);
432
433            vitfu{ima-1}=vitfu{ima-1}(i,1);
434            vitfu{ima-1}(:,2)=i;
435            vitfv{ima-1}=vitfv{ima-1}(i,1);
436            vitfv{ima-1}(:,2)=i;
437            i=1:1:size(vitfu{ima-1},1)
438            xpos=MatPos{2}(i,1)
439            ypos=MatPos{2}(i,2)
440      end
441    end
442
443
444
445    figure
446    hold on
447    plot(MatPos{1}(:,1),MatPos{1}(:,2),'r+')
448    plot(MatPos{2}(:,1),MatPos{2}(:,2),'b+')
449    plot(MatPos{4}(:,1),MatPos{4}(:,2),'y+')
450    quiver(xpos(:),ypos(:),vitfu{1}(:,1),vitfv{1}(:,1),'g')
451    quiver(MatPos{4}(:,1),MatPos{4}(:,2),vitfu{2}(:,1),vitfv{2}(:,1),'k')
452    legend('particules image 1','particules image 2', 'particules image 3','vitesse 1-2 (cm/s)','vitesse 2-3 (cm/s)');
453    xlabel('x (cm)');
454    ylabel('y (cm)');
455    title('Position et vitesse (cm/s) des particules')
456   
457
458    for i=1:size(vitfu{end},1)
459     vitfuadd(i)=0;
460     vitfvadd(i)=0;
461    end
462
463   
464   
465         for i=1:1:size(vitfu{end}(:,1))
466           
467                for j=1:nbfield-1
468                    vitfuadd(i)= vitfuadd(i)+vitfu{j}(i,1);
469                    vitfvadd(i)= vitfvadd(i)+vitfv{j}(i,1);
470                    xpos1(i)=MatPos{1}(i,1);
471                    ypos1(i)=MatPos{1}(i,2);
472                    xpos2(i)=MatPos{2}(i,1);
473                    ypos2(i)=MatPos{2}(i,2);
474                   
475                end
476            end
477            sizexpos1=size(xpos1)
478
479    vitfumoy=vitfuadd./(nbfield-1)
480    vitfvmoy=vitfvadd./(nbfield-1)
481
482    testresult1=result{1}
483    testresult2=result{2}
484   
485if nbfield>2   
486    figure
487    hold on
488    plot(MatPos{1}(:,1),MatPos{1}(:,2),'r+')
489    plot(MatPos{2}(:,1),MatPos{2}(:,2),'b+')
490    quiver(xpos2(:),ypos2(:),vitfumoy(:),vitfvmoy(:),'g')
491    legend('particules image 1','particules image 2', 'vitesse moyenne (cm/s)');
492    xlabel('x (cm)');
493    ylabel('y (cm)');
494    title('Position et vitesse (cm/s) des particules')
495   
496else
497
498    figure
499    hold on
500    plot(MatPos{1}(:,1),MatPos{1}(:,2),'r+')
501    plot(MatPos{2}(:,1),MatPos{2}(:,2),'b+')
502    quiver(MatPos{2}(:,1),MatPos{2}(:,2),vitfu{1}(:),vitfv{1}(:),'g')
503    legend('particules image 1','particules image 2','vitesse 1-2 (cm/s)');
504    xlabel('x (cm)');
505    ylabel('y (cm)');
506    title('Position et vitesse (cm/s) des particules')
507   
508    vitfumoy=vitfu{1};
509    vitfvmoy=vitfv{1};
510
511end
512
513VitData.NbDim=2;
514VitData.NbCoord=2;
515VitData.CoordType='phys';
516VitData.dt=0.0185;
517VitData.CoordUnit='cm';
518VitData.Z=0;
519VitData.ListDimName={'nb_vectors'};
520VitData.DimValue=size(vitfumoy,2);
521VitData.ListVarName={'X'  'Y'  'U'  'V'  'F'};
522VitData.VarDimIndex={[1]  [1]  [1]  [1]  [1]};
523VitData.ListVarAttribute={'Role'};
524VitData.Role={'coord_x'  'coord_y'  'vector_x'  'vector_y'  'warnflag'};
525
526if nbfield>2
527    VitData.X=size(MatPos{4},1);
528    VitData.Y=size(MatPos{4},2);
529else
530    VitData.X=size(MatPos{2},1);
531    VitData.Y=size(MatPos{2},2);
532end
533
534VitData.U=size(vitfumoy,2);
535VitData.V=size(vitfvmoy,2);
536VitData.Style='plane';
537VitData.Time=[198.5203 198.5203];
538VitData.Action=Series.Action;
539
540if nbfield>2
541    VitData.X=MatPos{4}(:,1)';
542    VitData.Y=MatPos{4}(:,2)';
543else
544    VitData.X=MatPos{2}(:,1)';
545    VitData.Y=MatPos{2}(:,2)';
546end
547
548VitData.U=vitfumoy(:)';
549VitData.V=vitfvmoy(:)';
550
551if length(VitData.ListVarName) >= 4 & isequal(VitData.ListVarName(1:4), {'X'  'Y'  'U'  'V'})
552       VitData.ListAttribute={'nb_coord','nb_dim','dt','pixcmx','pixcmy','hart','civ','fix'};
553       VitData.nb_coord=2;
554       VitData.nb_dim=2;
555       VitData.dt=0.018;
556       VitData.absolut_time_T0=0;
557       VitData.pixcmx=1; %pix per cm (1 by default)
558       VitData.pixcmy=1; %pix per cm (1 by default)
559       VitData.hart=0;
560           if isequal(VitData.CoordType,'px')
561             VitData.civ=1;
562           else
563             VitData.civ=0;
564           end
565        VitData.fix=0;
566        VitData.ListVarName(1:4)={'vec_X'  'vec_Y'  'vec_U'  'vec_V'};
567        VitData.vec_X=VitData.X;
568        VitData.vec_Y=VitData.Y;
569        VitData.vec_U=VitData.U;
570        VitData.vec_V=VitData.V;
571end
572currentdir=pwd;%store the current working directory
573[Path_ima,Name]=fileparts(filebase);%Path of the image files (.civ)
574cd(Path_ima);%move to the directory of the images: needed to create the result dir by 'mkdir'
575dircur=pwd; %current working directory
576[m1,m2,m3]=mkdir('TRACK_test')
577cd(currentdir)
578[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')
579error=struct2nc(filename_nc,VitData); %save result file
580if isequal(error,0)
581    [filename_nc ' written']
582else
583    warndlg_uvmat(error,'ERROR')
584end
585
Note: See TracBrowser for help on using the repository browser.