Changeset 1157 for trunk/src/series


Ignore:
Timestamp:
Jul 11, 2024, 4:13:03 PM (9 months ago)
Author:
sommeria
Message:

filter_time added

Location:
trunk/src/series
Files:
1 added
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/series/civ2vel_3C.m

    r1155 r1157  
    418418        MergeData.VarDimName={'coord_x','coord_y',{'coord_y','coord_x'},{'coord_y','coord_x'}...
    419419            {'coord_y','coord_x'},{'coord_y','coord_x'}};
     420        MergeData.VarAttribute{1}.Role='coord_x';
     421        MergeData.VarAttribute{2}.Role='coord_y';
     422        MergeData.VarAttribute{3}.Role='vector_x';
     423        MergeData.VarAttribute{4}.Role='vector_y';
     424        MergeData.VarAttribute{5}.Role='vector_z';
     425        MergeData.VarAttribute{6}.Role='ancillary';
    420426        MergeData.VarAttribute{6}.unit='pixel'; %error estimate expressed in pixel
    421427        if CheckZ

  • trunk/src/series/civ_3D.m

    r1156 r1157  
    130130end
    131131
    132 [tild,i1_series,i2_series,j1_series,j2_series]=get_file_series(Param);
     132[~,i1_series,~,j1_series,~]=get_file_series(Param);
    133133iview_A=0;% series index (iview) for the first image series
    134134iview_B=0;% series index (iview) for the second image series (only non zero for option 'shift' comparing two image series )
    … …  
    156156PairCiv1=Param.ActionInput.PairIndices.ListPairCiv1;
    157157
    158 [i1_series_Civ1,i2_series_Civ1,j1_series_Civ1,j2_series_Civ1,check_bounds,NomTypeNc]=...
     158[i1_series_Civ1,i2_series_Civ1,j1_series_Civ1,j2_series_Civ1,~,NomTypeNc]=...
    159159                        find_pair_indices(PairCiv1,i1_series{1},j1_series{1},MinIndex_i,MaxIndex_i,MinIndex_j,MaxIndex_j);
    160160                   
    … …  
    168168%% prepare output Data
    169169ListGlobalAttribute={'Conventions','Program','CivStage'};
    170 Data.Conventions='uvmat/civdata';% states the conventions used for the description of field variables and attributes
    171 Data.Program='civ_series';
     170Data.Conventions='uvmat/civdata_3D';% states the conventions used for the description of field variables and attributes
     171Data.Program='civ_3D';
     172if isfield(Param,'UvmatRevision')
     173    Data.Program=[Data.Program ', uvmat r' Param.UvmatRevision];
     174end
    172175Data.CivStage=0;%default
    173176list_param=(fieldnames(Param.ActionInput.Civ1))';
    174     list_param(strcmp('TestCiv1',list_param))=[];% remove the parameter TestCiv1 from the list
    175     Civ1_param=regexprep(list_param,'^.+','Civ1_$0');% insert 'Civ1_' before  each string in list_param
    176     Civ1_param=[{'Civ1_ImageA','Civ1_ImageB','Civ1_Time','Civ1_Dt'} Civ1_param]; %insert the names of the two input images
     177list_param(strcmp('TestCiv1',list_param))=[];% remove the parameter TestCiv1 from the list
     178Civ1_param=regexprep(list_param,'^.+','Civ1_$0');% insert 'Civ1_' before  each string in list_param
     179Civ1_param=[{'Civ1_ImageA','Civ1_ImageB','Civ1_Time','Civ1_Dt'} Civ1_param]; %insert the names of the two input images
    177180Data.ListGlobalAttribute=[ListGlobalAttribute Civ1_param];
    178             % set the list of variables
    179         Data.ListVarName={'Civ1_X','Civ1_Y','Civ1_U','Civ1_V','Civ1_W','Civ1_C','Civ1_FF'};%  cell array containing the names of the fields to record
    180         Data.VarDimName={'nb_vec_1','nb_vec_1','nb_vec_1','nb_vec_1','nb_vec_1','nb_vec_1','nb_vec_1'};
    181         Data.VarAttribute{1}.Role='coord_x';
    182         Data.VarAttribute{2}.Role='coord_y';
    183         Data.VarAttribute{3}.Role='vector_x';
    184         Data.VarAttribute{4}.Role='vector_y';
    185         Data.VarAttribute{5}.Role='vector_z';
    186         Data.VarAttribute{6}.Role='ancillary';
    187         Data.VarAttribute{7}.Role='errorflag';
     181% set the list of variables
     182Data.ListVarName={'Coord_z','Coord_y','Coord_x','Civ1_X','Civ1_Y','Civ1_U','Civ1_V','Civ1_W','Civ1_C','Civ1_FF'};%  cell array containing the names of the fields to record
     183Data.VarDimName={'npz','npy','npx',{'npz','npy','npx'},{'npz','npy','npx'},{'npz','npy','npx'},{'npz','npy','npx'},...
     184    {'npz','npy','npx'},{'npz','npy','npx'},{'npz','npy','npx'}};
     185Data.VarAttribute{1}.Role='coord_x';
     186Data.VarAttribute{2}.Role='coord_y';
     187Data.VarAttribute{3}.Role='vector_x';
     188Data.VarAttribute{4}.Role='vector_y';
     189Data.VarAttribute{5}.Role='vector_z';
     190Data.VarAttribute{6}.Role='ancillary';
     191Data.VarAttribute{7}.Role='errorflag';
     192Data.Coord_z=j1_series_Civ1;
     193% path for output nc file
     194OutputPath=fullfile(Param.OutputPath,num2str(Param.Experiment),num2str(Param.Device),[Param.OutputSubDir Param.OutputDirExt]);
    188195
    189196%% get timing from the ImaDoc file or input video
    … …  
    227234    par_civ1.ImageA=zeros(2*SearchRange_z+1,npy,npx);
    228235par_civ1.ImageB=zeros(2*SearchRange_z+1,npy,npx);
     236
     237
    229238
    230239%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    … …  
    268277   
    269278    Data.CivStage=1;
    270     % output nc file
    271     ncfile_out=fullfile_uvmat(OutputPath,Param.InputTable{1,2},Param.InputTable{1,3},Param.InputTable{1,5},...
    272         NomTypeNc,i1_series_Civ1(ifield),i2_series_Civ1(ifield),j1_series_Civ1(ifield),j2_series_Civ1(ifield));
     279 
     280    ncfile_out=fullfile_uvmat(OutputPath,'',Param.InputTable{1,3},'.nc',...
     281                '_1-1',i1_series_Civ1(ifield),i2_series_Civ1(ifield)); %output file name
    273282
    274283    if ~CheckOverwrite && exist(ncfile_out,'file')
    … …  
    284293        disp('civ1 started')
    285294
    286 
    287         % read input images (except in mode Test where it is introduced directly in Param.ActionInput.Civ1.ImageNameA and B)
    288 
    289         par_civ1.ImageA=zeros(2*SearchRange_z+1,npy,npx);
     295        % read input images by vertical blocks with nbre of images equal to 2*SearchRange+1,
     296        par_civ1.ImageA=zeros(2*SearchRange_z+1,npy,npx);%image block initiation
    290297        par_civ1.ImageB=zeros(2*SearchRange_z+1,npy,npx);
    291         %first vertical block centered at image index islice=par_civ1.Dz
    292         islice=par_civ1.Dz;
    293          for iz=1:par_civ1.Dz+SearchRange_z
    294               ImageName_A=fullfile_uvmat(RootPath_A,SubDir_A,RootFile_A,FileExt_A,NomType_A,i1_series_Civ1(1,ifield),[],j1_series_Civ1(iz,1));%
     298       
     299        z_index=1;%first vertical block centered at image index z_index=par_civ1.Dz
     300        for iz=1:par_civ1.Dz+SearchRange_z
     301            j_image_index=j1_series_Civ1(iz,1)% j index of the current image
     302            ImageName_A=fullfile_uvmat(RootPath_A,SubDir_A,RootFile_A,FileExt_A,NomType_A,i1_series_Civ1(1,ifield),[],j_image_index);%
    295303            A= read_image(ImageName_A,FileType_A);
    296             ImageName_B=fullfile_uvmat(RootPath_B,SubDir_B,RootFile_B,FileExt_B,NomType_B,i2_series_Civ1(1,ifield),[],j1_series_Civ1(iz,1));
     304            ImageName_B=fullfile_uvmat(RootPath_B,SubDir_B,RootFile_B,FileExt_B,NomType_B,i2_series_Civ1(1,ifield),[],j_image_index);
    297305            B= read_image(ImageName_B,FileType_B);
    298306            par_civ1.ImageA(iz+par_civ1.Dz-1,:,:) = A;
    299307            par_civ1.ImageB(iz+par_civ1.Dz-1,:,:) = B;
    300          end
    301          % caluclate velocity data (y and v in indices, reverse to y component)
    302             [Data.Civ1_X(islice,:,:),Data.Civ1_Y(islice,:,:), utable, vtable,wtable, ctable, FF, result_conv, errormsg] = civ3D (par_civ1);
     308        end
     309        % calculate velocity data at the first position z, corresponding to j_index=par_civ1.Dz
     310        [Data.Civ1_X(1,:,:),Data.Civ1_Y(1,:,:), Data.Civ1_U(1,:,:), Data.Civ1_V(1,:,:),Data.Civ1_W(1,:,:),...
     311            Data.Civ1_C(1,:,:), Data.Civ1_FF(1,:,:), result_conv, errormsg] = civ3D (par_civ1);
     312        if ~isempty(errormsg)
     313            disp_uvmat('ERROR',errormsg,checkrun)
     314            return
     315        end
     316        for z_index=2:floor((NbSlice-1)/par_civ1.Dz)% loop on slices
     317            par_civ1.ImageA=circshift(par_civ1.ImageA,-par_civ1.Dz);%shift the indices in the image block upward by par_civ1.Dz
     318            par_civ1.ImageB=circshift(par_civ1.ImageA,-par_civ1.Dz);
     319            for iz=1:par_civ1.Dz %read the new images at the end of the image block
     320                j_image_index=j1_series_Civ1(z_index*par_civ1.Dz+SearchRange_z-par_civ1.Dz+iz,1)
     321                ImageName_A=fullfile_uvmat(RootPath_A,SubDir_A,RootFile_A,FileExt_A,NomType_A,i1_series_Civ1(1,ifield),[],j_image_index);%
     322                A= read_image(ImageName_A,FileType_A);
     323                ImageName_B=fullfile_uvmat(RootPath_B,SubDir_B,RootFile_B,FileExt_B,NomType_B,i2_series_Civ1(1,ifield),[],j_image_index);
     324                B= read_image(ImageName_B,FileType_B);
     325                par_civ1.ImageA(2*SearchRange_z+1-par_civ1.Dz+iz,:,:) = A;
     326                par_civ1.ImageB(2*SearchRange_z+1-par_civ1.Dz+iz,:,:) = B;
     327            end
     328            % caluclate velocity data (y and v in indices, reverse to y component)
     329            [Data.Civ1_X(z_index,:,:),Data.Civ1_Y(z_index,:,:), Data.Civ1_U(z_index,:,:), Data.Civ1_V(z_index,:,:),Data.Civ1_W(z_index,:,:),...
     330                Data.Civ1_C(z_index,:,:), Data.Civ1_FF(z_index,:,:), result_conv, errormsg] = civ3D (par_civ1);
     331
    303332            if ~isempty(errormsg)
    304333                disp_uvmat('ERROR',errormsg,checkrun)
    305334                return
    306335            end
    307         for islice=2*par_civ1.Dz:NbSlice% loop on slices for the first image in the pair
    308             par_civ1.ImageA=circshift(par_civ1.ImageA,-par_civ1.Dz);%shift the indces in the block upward by par_civ1.Dz
    309             par_civ1.ImageB=circshift(par_civ1.ImageA,-par_civ1.Dz);
    310               for iz=1:par_civ1.Dz
    311             ImageName_A=fullfile_uvmat(RootPath_A,SubDir_A,RootFile_A,FileExt_A,NomType_A,i1_series_Civ1(1,ifield),[],j1_series_Civ1(islice+SearchRange_z-par_civ1.Dz+iz,1));%
    312             A= read_image(ImageName_A,FileType_A);
    313             ImageName_B=fullfile_uvmat(RootPath_B,SubDir_B,RootFile_B,FileExt_B,NomType_B,i2_series_Civ1(1,ifield),[],j1_series_Civ1(islice+SearchRange_z-par_civ1.Dz+iz,1));
    314             B= read_image(ImageName_B,FileType_B);
    315             par_civ1.ImageA(2*SearchRange_z+1-par_civ1.Dz+iz,:,:) = A;
    316             par_civ1.ImageB(2*SearchRange_z+1-par_civ1.Dz+iz,:,:) = B;
    317             end
    318             % caluclate velocity data (y and v in indices, reverse to y component)
    319             [Data.Civ1_X(islice,:,:),Data.Civ1_Y(islice,:,:), utable, vtable,wtable, ctable, FF, result_conv, errormsg] = civ3D (par_civ1);
    320             if ~isempty(errormsg)
    321                 disp_uvmat('ERROR',errormsg,checkrun)
    322                 return
    323             end
    324 
    325         end
     336        end
     337        Data.Civ1_C=double(Data.Civ1_C);
     338        Data.Civ1_C(Data.Civ1_C==Inf)=0;
     339        [npz,npy,npx]=size(Data.Civ1_X);
     340        Data.Coord_z=1:npz;
     341        Data.Coord_y=1:npy;
     342        Data.Coord_x=1:npx;
    326343        % case of mask TO ADAPT
    327344        if par_civ1.CheckMask&&~isempty(par_civ1.Mask)
    … …  
    668685        end
    669686        Data.ListGlobalAttribute=[Data.ListGlobalAttribute Fix2_param];
    670         Data.Civ2_FF=double(detect_false(Param.ActionInput.Fix2,Data.Civ2_C,Data.Civ2_U,Data.Civ2_V,Data.Civ2_FF));
     687        Data.Civ2_FF=detect_false(Param.ActionInput.Fix2,Data.Civ2_C,Data.Civ2_U,Data.Civ2_V,Data.Civ2_FF);
    671688        Data.CivStage=Data.CivStage+1;
    672689    end
    … …  
    728745    time_total=toc(tstart);
    729746    disp(['ellapsed time ' num2str(time_total/60,2) ' minutes'])
    730     disp(['time image reading ' num2str(time_input,2) ' s'])
    731747    disp(['time civ1 ' num2str(time_civ1,2) ' s'])
    732748    disp(['time patch1 ' num2str(time_patch1,2) ' s'])
    733749    disp(['time civ2 ' num2str(time_civ2,2) ' s'])
    734750    disp(['time patch2 ' num2str(time_patch2,2) ' s'])
    735     disp(['time other ' num2str((time_total-time_input-time_civ1-time_patch1-time_civ2-time_patch2),2) ' s'])
    736     % end
     751
    737752end
    738753
    … …  
    860875    for ivec_x=1:nbvec_x
    861876        for ivec_y=1:nbvec_y
    862             ivec_y
    863             iref=round(par_civ.Grid(ivec_y,ivec_x,1)+0.5)% xindex on the image A for the middle of the correlation box
    864             jref=round(par_civ.ImageHeight-par_civ.Grid(ivec_y,ivec_x,2)+0.5)%  j index  for the middle of the correlation box in the image A
     877            iref=round(par_civ.Grid(ivec_y,ivec_x,1)+0.5);% xindex on the image A for the middle of the correlation box
     878            jref=round(par_civ.ImageHeight-par_civ.Grid(ivec_y,ivec_x,2)+0.5);%  j index  for the middle of the correlation box in the image A
    865879            subrange1_x=iref-ibx2:iref+ibx2;% x indices defining the first subimage
    866880            subrange1_y=jref-iby2:jref+iby2;% y indices defining the first subimage
    … …  
    918932
    919933               
    920                     %reference: Oliver Pust, PIV: Direct Cross-Correlation
     934                   
    921935                    for kz=1:par_civ.SearchBoxSize(3)
    922936                        subima2=squeeze(image2_crop(kz,:,:));
    … …  
    928942               
    929943                    end
    930                     [corrmax,z]=max(max_xy);
     944                    [corrmax,zmax]=max(max_xy);
    931945               
    932                     x=xk(z);
    933                     y=yk(z);
     946                    x=xk(zmax);
     947                    y=yk(zmax);
    934948                    result_conv=(result_conv/corrmax)*255; %normalize, peak=always 255
    935                     subimage2_crop=squeeze(image2_crop(z,y:y+2*iby2/mesh,x:x+2*ibx2/mesh));%subimage of image 2 corresponding to the optimum displacement of first image
    936                     sum_square=sum(sum(squeeze(image1_crop(z,:,:).*image1_crop(z,:,:))));
     949                    subimage2_crop=squeeze(image2_crop(zmax,y:y+2*iby2/mesh,x:x+2*ibx2/mesh));%subimage of image 2 corresponding to the optimum displacement of first image
     950                    sum_square=sum(sum(squeeze(image1_crop(zmax,:,:).*image1_crop(zmax,:,:))));
    937951                    sum_square=sum_square*sum(sum(subimage2_crop.*subimage2_crop));% product of variances of image 1 and 2
    938952                    sum_square=sqrt(sum_square);% srt of the variance product to normalise correlation
    … …  
    940954       
    941955                            if par_civ.CorrSmooth==1
    942                                 [vector,FF(ivec_y,ivec_x)] = SUBPIXGAUSS (result_conv(z,:,:),x,y);%TODO: improve by max optimisation along z
     956                                [vector,FF(ivec_y,ivec_x)] = SUBPIXGAUSS (squeeze(result_conv(zmax,:,:)),x,y);%TODO: improve by max optimisation along z
    943957                            elseif par_civ.CorrSmooth==2
    944                                 [vector,FF(ivec_y,ivec_x)] = SUBPIX2DGAUSS (result_conv(z,:,:),x,y);
     958                                [vector,FF(ivec_y,ivec_x)] = SUBPIX2DGAUSS (squeeze(result_conv(zmax,:,:)),x,y);
    945959                            else
    946                                 [vector,FF(ivec_y,ivec_x)] = quadr_fit(result_conv(z,:,:),x,y);
     960                                [vector,FF(ivec_y,ivec_x)] = quadr_fit(squeeze(result_conv(zmax,:,:)),x,y);
    947961                            end
    948962                            utable(ivec_y,ivec_x)=vector(1)*mesh+shiftx(ivec_y,ivec_x);
    … …  
    952966                            iref=round(xtable(ivec_y,ivec_x)+0.5);% nearest image index for the middle of the vector
    953967                            jref=round(ytable(ivec_y,ivec_x)+0.5);
    954                             wtable(ivec_y,ivec_x)=z-kref;
     968                            wtable(ivec_y,ivec_x)=zmax-kref;
    955969                            % eliminate vectors located in the mask
    956970                            if  checkmask && (iref<1 || jref<1 ||iref>npx_ima || jref>npy_ima ||( par_civ.Mask(jref,iref)<200 && par_civ.Mask(jref,iref)>=100))
    … …  
    962976
    963977                    else
    964                         FF(ivec_y,ivec_x)=1;
     978                        FF(ivec_y,ivec_x)=true;
    965979                    end
    966980                end
    … …  
    975989% OUPUT:
    976990% vector = optimum displacement vector with subpixel correction
    977 % F =flag: =0 OK
    978 %           =-2 , warning: max too close to the edge of the search box (1 pixel margin)
     991% FF =flag: ='true' max too close to the edge of the search box (1 pixel margin)
    979992% INPUT:
    980993% x,y: position of the maximum correlation at integer values
    981994
    982 function [vector,F] = SUBPIXGAUSS (result_conv,x,y)
     995function [vector,FF] = SUBPIXGAUSS (result_conv,x,y)
    983996%------------------------------------------------------------------------
    984 % vector=[0 0]; %default
    985 F=0;
     997vector=[0 0]; %default
    986998[npy,npx]=size(result_conv);
    987999result_conv(result_conv<1)=1; %set to 1 correlation values smaller than 1  (=0 by discretisation, to avoid divergence in the log)
    9881000%the following 8 lines are copyright (c) 1998, Uri Shavit, Roi Gurka, Alex Liberzon, Technion ??? Israel Institute of Technology
    9891001%http://urapiv.wordpress.com
    990 peaky = y;
    991 if y < npy && y > 1
    992     f0 = log(result_conv(y,x));
    993     f1 = log(result_conv(y-1,x));
    994     f2 = log(result_conv(y+1,x));
    995     peaky = peaky+ (f1-f2)/(2*f1-4*f0+2*f2);
    996 else
    997     F=1; % warning flag for vector truncated by the limited search box
    998 end
    999 peakx=x;
    1000 if x < npx-1 && x > 1
    1001     f0 = log(result_conv(y,x));
    1002     f1 = log(result_conv(y,x-1));
    1003     f2 = log(result_conv(y,x+1));
    1004     peakx = peakx+ (f1-f2)/(2*f1-4*f0+2*f2);
    1005 else
    1006     F=1; % warning flag for vector truncated by the limited search box
    1007 end
    1008 vector=[peakx-floor(npx/2)-1 peaky-floor(npy/2)-1];
     1002FF=false;
     1003    peaky = y;
     1004    if y < npy && y > 1
     1005        f0 = log(result_conv(y,x));
     1006        f1 = log(result_conv(y-1,x));
     1007        f2 = log(result_conv(y+1,x));
     1008        peaky = peaky+ (f1-f2)/(2*f1-4*f0+2*f2);
     1009    else
     1010        FF=true; % error flag for vector truncated by the limited search box in y
     1011    end
     1012    peakx=x;
     1013    if x < npx-1 && x > 1
     1014        f0 = log(result_conv(y,x));
     1015        f1 = log(result_conv(y,x-1));
     1016        f2 = log(result_conv(y,x+1));
     1017        peakx = peakx+ (f1-f2)/(2*f1-4*f0+2*f2);
     1018    else
     1019        FF=true; % warning flag for vector truncated by the limited search box in x
     1020    end
     1021    vector=[peakx-floor(npx/2)-1 peaky-floor(npy/2)-1];
    10091022
    10101023%------------------------------------------------------------------------

  • trunk/src/series/civ_series.m

    r1156 r1157  
    301301    end
    302302    if CheckInputFile
    303         OutputPath=fullfile(Param.OutputPath,num2str(Param.Experiment),num2str(Param.Device));
     303        OutputPath=fullfile(Param.OutputPath,Param.Experiment,Param.Device);
    304304        if iview_A==0 % no nc file has been entered
    305305            ncfile=fullfile_uvmat(OutputPath,Param.InputTable{1,2},Param.InputTable{1,3},Param.InputTable{1,5},...

  • trunk/src/series/sliding_average.m

    r1127 r1157  
    1 %'sliding_average_natalya'
     1%'sliding_average': calculate the time correlation function at each point
    22%------------------------------------------------------------------------
    3 % function ParamOut=sliding_average(Param)
     3% function ParamOut=turb_correlation_time(Param)
    44%
    55%%%%%%%%%%% GENERAL TO ALL SERIES ACTION FCTS %%%%%%%%%%%%%%%%%%%%%%%%%%%
    … …  
    3939%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    4040%=======================================================================
    41 % Copyright 2008-2024, LEGI UMR 5519 / CNRS UGA G-INP, Grenoble, France
     41% Copyright 2008-2022, LEGI UMR 5519 / CNRS UGA G-INP, Grenoble, France
    4242%   http://www.legi.grenoble-inp.fr
    43 %   Joel.Sommeria - Joel.Sommeria (A) univ-grenoble-alpes.fr
     43%   Joel.Sommeria - Joel.Sommeria (A) legi.cnrs.fr
    4444%
    4545%     This file is part of the toolbox UVMAT.
    … …  
    6565    ParamOut.VelType='off';% menu for selecting the velocity type (options 'off'/'one'/'two',  'off' by default)
    6666    ParamOut.FieldName='one';% menu for selecting the field (s) in the input file(options 'off'/'one'/'two', 'off' by default)
    67     ParamOut.FieldTransform = 'on';%can use a transform function
     67    ParamOut.FieldTransform = 'off';%can use a transform function
    6868    ParamOut.ProjObject='off';%can use projection object39(option 'off'/'on',
    6969    ParamOut.Mask='off';%can use mask option   (option 'off'/'on', 'off' by default)
    7070    ParamOut.OutputDirExt='.tfilter';%set the output dir extension
    71     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
    72 %     filecell=get_file_series(Param);%check existence of the first input file
    73 %     if ~exist(filecell{1,1},'file')
    74 %         msgbox_uvmat('WARNING','the first input file does not exist')
    75 %     end
     71    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
     72   
     73    %% check the first input file in the series
     74     first_j=[];%
     75    if isfield(Param.IndexRange,'first_j'); first_j=Param.IndexRange.first_j; end
     76    PairString='';
     77    if isfield(Param.IndexRange,'PairString'); PairString=Param.IndexRange.PairString; end
     78    [i1,i2,j1,j2] = get_file_index(Param.IndexRange.first_i,first_j,PairString);
     79    FirstFileName=fullfile_uvmat(Param.InputTable{1,1},Param.InputTable{1,2},Param.InputTable{1,3},...
     80        Param.InputTable{1,5},Param.InputTable{1,4},i1,i2,j1,j2);
     81    if exist(FirstFileName,'file')
     82        FileInfo=get_file_info(FirstFileName);
     83        if ~(isfield(FileInfo,'FileType')&& strcmp(FileInfo.FileType,'netcdf'))
     84            msgbox_uvmat('ERROR','this fct works only on netcdf files (fields projected on a grid, notraw civ data)')
     85        return
     86        end
     87    else
     88        msgbox_uvmat('ERROR',['the input file ' FirstFileName ' does not exist'])
     89        return
     90    end
     91
     92    %% setting of intput parameters
     93%     answer = msgbox_uvmat('INPUT_MENU','select the variables to filter',ListVarName');
     94%     ParamOut.ListVarName=ListVarName(answer);
     95    answer = msgbox_uvmat('INPUT_TXT','set the filering window length','3');
     96    ParamOut.WindowLength=str2double(answer);
     97   
     98%         [ParamOut.ActionInput,errormsg] = set_param_input(ListParam,DefaultValue,ParamIn);
     99%         if ~isempty(errormsg)
     100%             msgbox_uvmat('ERROR',errormsg)
     101%         end
     102%     [ParamOut,errormsg] = set_param_input(ListParam,DefaultValue,ParamIn);
    76103    return
    77104end
    … …  
    184211nbmissing=0;
    185212
    186 %% initialisation manip Coriolis
    187 char_index=regexp(SubDir{1},'waves_L1_');
    188 switch(SubDir{1}(char_index+9))
    189     case '1'
    190         amplitude=2.5 %oscillation amplitude
    191         T=24.46;
    192         t0=3 ;% dt=0.5 s, torus at its max x at the beginning of motion, i0=7
    193     case '2'
    194         amplitude=5 %oscillation amplitude
    195         T=24.47;
    196         t0=8.5; % dt=1/3 s -> image index of starting motion = 26, % torus at its max x at the beginning of motion
    197     case '3'
    198         amplitude=10 %oscillation amplitude
    199         T=24.45;
    200         t0=6.5-T/2;% dt=0.25, torus at its minimum x at the beginning of motion
    201     case '4' 
    202         amplitude=15 %oscillation amplitude
    203         T=24.48;
    204         t0=3.4;     %dt=0.2 -> i0=18 image index of starting motion, % torus at its max x at the beginning of motion
    205 end
    206 %NbPeriod=2; %number of periods for the sliding average
    207 NbPeriod=4; %number of periods for the sliding average
    208 omega=2*2*pi/T;%harmonic
    209 Lscale=15;%diameter of the torus, length scale for normalisation
    210 Uscale=amplitude*omega;
    211 
    212 DataOut.ListGlobalAttribute= {'Conventions','Time'};
     213%% initialisation output data
     214
     215DataOut.ListGlobalAttribute= {'Conventions'};
    213216DataOut.Conventions='uvmat';
    214 DataOut.ListVarName={'coord_y','coord_x','Umean','Vmean','Ucos','Vcos','Usin','Vsin','DUDXsin','DUDXcos','DUDYsin','DVDXsin','DVDXcos'...
    215     ,'DVDYsin','Ustokes','Vstokes'};
    216 DataOut.VarDimName={'coord_y','coord_x',{'coord_y','coord_x'},{'coord_y','coord_x'},{'coord_y','coord_x'},{'coord_y','coord_x'},...
    217     {'coord_y','coord_x'},{'coord_y','coord_x'},{'coord_y','coord_x'},{'coord_y','coord_x'},{'coord_y','coord_x'},...
    218     {'coord_y','coord_x'},{'coord_y','coord_x'},{'coord_y','coord_x'},{'coord_y','coord_x'},{'coord_y','coord_x'}};
     217DataOut.ListVarName={'Time','coord_y','coord_x','Ufilter','Vfilter'};
     218DataOut.VarDimName={'Time','coord_y','coord_x',{'Time','coord_y','coord_x'},{'Time','coord_y','coord_x'}};
    219219
    220220%%%%%%%%%%%%%%%% loop on field indices %%%%%%%%%%%%%%%%
    … …  
    226226    return
    227227end
    228 [Data,tild,errormsg]=nc2struct(filecell{1,end});   
    229 Time_end=Data.Time;
    230 dt=(Time_end-Time_1)/(NbField-1); %time interval
    231 NpTime=round(NbPeriod*T/dt+1);
     228[Data,tild,errormsg]=nc2struct(filecell{1,2});   
     229Time_2=Data.Time;
     230dt=(Time_2-Time_1); %time interval
     231NpTime=21; %Nbre de champ pour la moyenne glissante (choisir impair)
    232232
    233233OutputPath=fullfile(Param.OutputPath,num2str(Param.Experiment),num2str(Param.Device));
    … …  
    246246    %%%%%%%%%%% MAIN RUNNING OPERATIONS  %%%%%%%%%%%%
    247247    if index==1 %first field
    248         DataOut.coord_x=Field.coord_x/Lscale;
    249         DataOut.coord_y=Field.coord_y/Lscale;
     248        DataOut.coord_x=Field.coord_x;
     249        DataOut.coord_y=Field.coord_y;
    250250        npy=numel(DataOut.coord_y);
    251251        npx=numel(DataOut.coord_x);
    252         Umean=zeros(NpTime,npy,npx);
    253         Vmean=zeros(NpTime,npy,npx);
    254         Ucos=zeros(NpTime,npy,npx);
    255         Vcos=zeros(NpTime,npy,npx);
    256         Usin=zeros(NpTime,npy,npx);
    257         Vsin=zeros(NpTime,npy,npx);
    258         DUDXcos=zeros(NpTime,npy,npx);
    259         DUDXsin=zeros(NpTime,npy,npx);
    260         DUDYsin=zeros(NpTime,npy,npx);
    261         DVDXcos=zeros(NpTime,npy,npx);
    262         DVDXsin=zeros(NpTime,npy,npx);
    263         DVDYsin=zeros(NpTime,npy,npx);
    264     end
    265     Time(index)=Field.Time-t0;%time from the start of the motion
    266     Umean=circshift(Umean,[-1 0 0]); %shift U by ishift along the first index
    267     Vmean=circshift(Vmean,[-1 0 0]); %shift U by ishift along the first index
    268     Ucos=circshift(Ucos,[-1 0 0]); %shift U by ishift along the first index
    269     Vcos=circshift(Vcos,[-1 0 0]); %shift U by ishift along the first index
    270     Usin=circshift(Usin,[-1 0 0]); %shift U by ishift along the first index
    271     Vsin=circshift(Vsin,[-1 0 0]); %shift U by ishift along the first index
    272     DUDXcos=circshift(DUDXcos,[-1 0 0]);
    273     DUDXsin=circshift(DUDXsin,[-1 0 0]);
    274     DUDYsin=circshift(DUDYsin,[-1 0 0]);       
    275     DVDXcos=circshift(DVDXcos,[-1 0 0]);
    276     DVDXsin=circshift(DVDXsin,[-1 0 0]);
    277     DVDYsin=circshift(DVDYsin,[-1 0 0]);       
    278     Umean(end,:,:)=Field.U;
    279     Vmean(end,:,:)=Field.V;
    280     Ucos(end,:,:)=Field.U*cos(omega*Time(index));
    281     Vcos(end,:,:)=Field.V*cos(omega*Time(index));
    282     Usin(end,:,:)=Field.U*sin(omega*Time(index));
    283     Vsin(end,:,:)=Field.V*sin(omega*Time(index));
    284     DUDXcos(end,:,:)=Field.DUDX*cos(omega*Time(index));
    285     DUDXsin(end,:,:)=Field.DUDX*sin(omega*Time(index));
    286     DUDYsin(end,:,:)=Field.DUDY*sin(omega*Time(index));% ParamOut=[];%default output
    287 
    288     DVDXcos(end,:,:)=Field.DVDX*cos(omega*Time(index));
    289     DVDXsin(end,:,:)=Field.DVDX*sin(omega*Time(index));
    290     DVDYsin(end,:,:)=Field.DVDY*sin(omega*Time(index));
    291     DataOut.Time=(Time(index)-(NpTime-1)*dt/2)/T;%time inperiods from the beginning of the oscillation (torus at max abscissa)
    292     DataOut.Umean=(1/Uscale)*squeeze(nanmean(Umean,1));
    293     DataOut.Vmean=(1/Uscale)*squeeze(nanmean(Vmean,1));
    294     DataOut.Ucos=2*(1/Uscale)*squeeze(nanmean(Ucos,1));
    295     DataOut.Vcos=2*(1/Uscale)*squeeze(nanmean(Vcos,1));
    296     DataOut.Usin=2*(1/Uscale)*squeeze(nanmean(Usin,1));
    297     DataOut.Vsin=2*(1/Uscale)*squeeze(nanmean(Vsin,1));
    298     DataOut.DUDXcos=2*squeeze(nanmean(DUDXcos,1));
    299     DataOut.DUDXsin=2*squeeze(nanmean(DUDXsin,1));
    300     DataOut.DUDYsin=2*squeeze(nanmean(DUDYsin,1));
    301     DataOut.DVDXcos=2*squeeze(nanmean(DVDXcos,1));
    302     DataOut.DVDXsin=2*squeeze(nanmean(DVDXsin,1));
    303     DataOut.DVDYsin=2*squeeze(nanmean(DVDYsin,1));
    304     DataOut.Ustokes=(1/omega)*(1/Uscale)*(DataOut.Ucos.*DataOut.DUDXsin+DataOut.Vcos.*DataOut.DUDYsin);
    305     DataOut.Vstokes=(1/omega)*(1/Uscale)*(DataOut.Ucos.*DataOut.DVDXsin+DataOut.Vcos.*DataOut.DVDYsin);
     252        Ufilter=zeros(NpTime,npy,npx);
     253        Vfilter=zeros(NpTime,npy,npx);
     254    end
     255    Time(index)=Field.Time;%time
     256    Ufilter=circshift(Ufilter,[-1 0 0]); %shift U by ishift along the first index
     257    Vfilter=circshift(Vfilter,[-1 0 0]); %shift U by ishift along the first index   
     258    Ufilter(end,:,:)=Field.U;
     259    Vfilter(end,:,:)=Field.V;
     260    DataOut.Time=(Time(index)-(NpTime-1)*dt/2);%mid time
     261    DataOut.Ufilter=squeeze(nanmean(Ufilter,1));
     262    DataOut.Vfilter=squeeze(nanmean(Vfilter,1));
     263   
    306264
    307265    % writing the result file as netcdf file
    308     i1=i1_series{1}(index);
     266    i1=i1_series{1}(index)-ceil(NpTime/2);
    309267    OutputFile=fullfile_uvmat(OutputPath,OutputDir,RootFileOut,'.nc',NomTypeOut,i1);
    310268    errormsg=struct2nc(OutputFile, DataOut);

  • trunk/src/series/test_filter_tps.m

    r1156 r1157  
    7272    % check the existence of the first file in the series
    7373     first_j=[];% note that the function will propose to cover the whole range of indices
    74     if isfield(Param.IndexRange,'MinIndex_j'); first_j=Param.IndexRange.MinIndex_j; end
    75     last_j=[];
    76     if isfield(Param.IndexRange,'MaxIndex_j'); last_j=Param.IndexRange.MaxIndex_j; end
     74    if isfield(Param.IndexRange,'first_j'); first_j=Param.IndexRange.first_j; end
    7775    PairString='';
    7876    if isfield(Param.IndexRange,'PairString'); PairString=Param.IndexRange.PairString; end
    … …  
    8179        Param.InputTable{1,5},Param.InputTable{1,4},i1,i2,j1,j2);
    8280    if ~exist(FirstFileName,'file')
    83         msgbox_uvmat('WARNING',['the first input file ' FirstFileName ' does not exist'])
     81        msgbox_uvmat('ERROR',['the input file ' FirstFileName ' does not exist'])
    8482        return
    85     else
    86         [i1,i2,j1,j2] = get_file_index(Param.IndexRange.last_i,last_j,PairString);
    87         LastFileName=fullfile_uvmat(Param.InputTable{1,1},Param.InputTable{1,2},Param.InputTable{1,3},...
    88         Param.InputTable{1,5},Param.InputTable{1,4},i1,i2,j1,j2);
    89         if ~exist(FirstFileName,'file')
    90              msgbox_uvmat('WARNING',['the last input file ' LastFileName ' does not exist'])
    91         end
    9283    end
    9384
    … …  
    114105        end
    115106    else
    116         msgbox_uvmat('ERROR',['invalid file type input: ' FileType ' not a civ data'])
     107        msgbox_uvmat('ERROR','invalid file type input: test_filter_tps proceeds raw civ data')
    117108        return
    118109    end
Note: See TracChangeset for help on using the changeset viewer.