Changeset 1157

Timestamp:

Jul 11, 2024, 4:13:03 PM (3 months ago)

Author:

sommeria

Message:

filter_time added

Location:

trunk/src

Files:

• fullfile_uvmat.m (modified) (1 diff)
• get_field.m (modified) (2 diffs)
• get_file_info.m (modified) (3 diffs)
• msgbox_uvmat.m (modified) (4 diffs)
• nc2struct.m (modified) (1 diff)
• series/civ2vel_3C.m (modified) (1 diff)
• series/civ_3D.m (modified) (15 diffs)
• series/civ_series.m (modified) (1 diff)
• series/filter_time.m (added)
• series/sliding_average.m (modified) (6 diffs)
• series/test_filter_tps.m (modified) (3 diffs)
• set_param_input.m (modified) (2 diffs)
• struct2nc.m (modified) (5 diffs)

trunk/src/fullfile_uvmat.m

@@ -63 +63 @@
 
 %% default input
+if iscell(NomType)
+NomType=NomType{1}
+end
 if ~exist('j2','var') 
     j2=[];

trunk/src/get_field.m

@@ -509 +509 @@
 VarIndex=find(strcmp(VarName,Field.Display.ListVarName),1);
 DimCell=Field.Display.VarDimName{VarIndex};
-% y_index=get(handles.Coord_y,'Value');
-% y_menu=get(handles.Coord_y,'String');
-% if isempty(y_menu)
-%     return
-% else
-% YName=y_menu{y_index};
-% end
+
 
 %% set list of possible coordinates
-% test_component=zeros(size(Field.Display.VarDimName));%=1 when variable #ilist is eligible as unstructured coordinate
+
 test_coord=zeros(size(Field.Display.VarDimName)); %=1 when variable #ilist is eligible as structured coordiante
-% ListCoord={''};
-% dim_var=Field.Display.VarDimName{y_index};%list of dimensions of the selected variable
 
 for ilist=1:numel(Field.Display.VarDimName)
-    dimnames=Field.Display.VarDimName{ilist}; %list of dimensions for variable #ilist
-    if isequal(dimnames,DimCell)||isequal(dimnames(1:end-1),DimCell)||isequal(dimnames(2:end),DimCell)
-        test_coord(ilist)=1;
-    end
-end
-ListCoord=Field.Display.ListVarName(find(test_coord));
+    %dimnames=Field.Display.VarDimName{ilist}; %list of dimensions for variable #ilist
+    % if isequal(dimnames,DimCell)||isequal(dimnames(1:end-1),DimCell)||isequal(dimnames(2:end),DimCell)
+    % if numel(dimnames)==1 || 
+    %     test_coord(ilist)=1;
+    % end
+end
+ListCoord=Field.Display.ListVarName;%(find(test_coord));
 set(handles.Coord_y,'String',ListCoord)
 val_y=1;
@@ -984 +977 @@
 %-----------------------------------------------------------------------
 Field=get(handles.get_field,'UserData');
-index=name2index(VarName,Field.ListVarName);
+index=name2index(VarName,Field.ListVarName); %index of the selectd variable
 if ~isempty(index)
-    set(handles.variables,'Value',index+1)
+    set(handles.variables,'Value',index+1) %indicate which variable is selected in the list (+1 because of the '*' display)
     variables_Callback(handles.variables, VarName, handles)
 end

trunk/src/get_file_info.m

@@ -78 +78 @@
         FileInfo.FileType=regexprep(FileExt,'^.','');% eliminate the dot of the extension;
     case {'.seq','.sqb'}
-        [A,FileInfo,timestamps,errormsg]=read_rdvision(fileinput,[]);
+        [~,FileInfo,timestamps,errormsg]=read_rdvision(fileinput,[]);
     case '.im7'
         try
@@ -89 +89 @@
              FileInfo.TimeName='timestamp';
              for ilist=1:numel(Input.Attributes)
-                 if strcmp(Input.Attributes{ilist}.Name,'_Date')
-                     DateString=Input.Attributes{ilist}.Value;
-                 end
-                 if strcmp(Input.Attributes{ilist}.Name,'_Time')
-                     TimeString=Input.Attributes{ilist}.Value;
-                 end
+                 % if strcmp(Input.Attributes{ilist}.Name,'_Date')
+                 %     DateString=Input.Attributes{ilist}.Value;
+                 % end
+                 % if strcmp(Input.Attributes{ilist}.Name,'_Time')
+                 %     TimeString=Input.Attributes{ilist}.Value;
+                 % end
              end
         catch ME
@@ -175 +175 @@
                                 FileInfo.FileType='netcdf';
                                 FileInfo.ListVarName=Data.ListVarName;
+                                FileInfo.VarAttribute=Data.VarAttribute;
                             end
                         end

trunk/src/msgbox_uvmat.m

@@ -10 +10 @@
 % title: string indicating the type of message box (the title is displayed in the upper bar of the fig):
 %                ='INPUT_TXT'(default), input data is asked in an edit box
+%                = 'INPUT_MENU', input data must be selected from  a menu choice
 %                ='CONFIMATION'', 'ERROR', 'WARNING','RULER' the figure remains  opened until a button 'OK' is pressed
 %                ='RULER' is used for display of length and angle from the ruler tool. 
@@ -16 +17 @@
 %                ='WAITING...' the figure remains open until the program deletes it
 % display: displayed text
-% default_answer: default answer in the edit box (only used with title='INPUT_TXT')
+% default_answer: default answer in the edit box (only used with title='INPUT_TXT' or 'INPUT_MENU')
 
 %=======================================================================
@@ -69 +70 @@
 set(handles.figure1,'Units','pixels')
 FigPos=[100 150 500 50];%default position
+if strcmp(title,'INPUT_MENU')
+    FigPos=[100 150 500 200];%default position
+end
 if exist('Position','var') && numel(Position)>=2
     FigPos(1)=Position(1);
@@ -164 +168 @@
     set(handles.edit_box, 'Visible', 'on');
     set(handles.edit_box,'String', default_answer)
+    set(handles.edit_box,'Max',2);% allows for multiple selection
 else
     set(handles.text1, 'Position', [0.15 0.3 0.85 0.7]);

trunk/src/nc2struct.m

@@ -27 +27 @@
 %  additional arguments:
 %       -no additional arguments: all the variables of the NetCDF file are read.
+%       - empty argument []: the field structure with the names of variables is read, without their values
 %       -a cell array, ListVarName, made of  char strings {'VarName1', 'VarName2',...} )
 %         if ListVarName=[] or {}, no variable value is read (only global attributes and list of variables and dimensions)

trunk/src/series/civ2vel_3C.m

@@ -418 +418 @@
         MergeData.VarDimName={'coord_x','coord_y',{'coord_y','coord_x'},{'coord_y','coord_x'}...
             {'coord_y','coord_x'},{'coord_y','coord_x'}};
+        MergeData.VarAttribute{1}.Role='coord_x';
+        MergeData.VarAttribute{2}.Role='coord_y';
+        MergeData.VarAttribute{3}.Role='vector_x';
+        MergeData.VarAttribute{4}.Role='vector_y';
+        MergeData.VarAttribute{5}.Role='vector_z';
+        MergeData.VarAttribute{6}.Role='ancillary';
         MergeData.VarAttribute{6}.unit='pixel'; %error estimate expressed in pixel
         if CheckZ

trunk/src/series/civ_3D.m

@@ -130 +130 @@
 end
 
-[tild,i1_series,i2_series,j1_series,j2_series]=get_file_series(Param);
+[~,i1_series,~,j1_series,~]=get_file_series(Param);
 iview_A=0;% series index (iview) for the first image series
 iview_B=0;% series index (iview) for the second image series (only non zero for option 'shift' comparing two image series )
@@ -156 +156 @@
 PairCiv1=Param.ActionInput.PairIndices.ListPairCiv1;
 
-[i1_series_Civ1,i2_series_Civ1,j1_series_Civ1,j2_series_Civ1,check_bounds,NomTypeNc]=...
+[i1_series_Civ1,i2_series_Civ1,j1_series_Civ1,j2_series_Civ1,~,NomTypeNc]=...
                         find_pair_indices(PairCiv1,i1_series{1},j1_series{1},MinIndex_i,MaxIndex_i,MinIndex_j,MaxIndex_j);
                     
@@ -168 +168 @@
 %% prepare output Data
 ListGlobalAttribute={'Conventions','Program','CivStage'};
-Data.Conventions='uvmat/civdata';% states the conventions used for the description of field variables and attributes
-Data.Program='civ_series';
+Data.Conventions='uvmat/civdata_3D';% states the conventions used for the description of field variables and attributes
+Data.Program='civ_3D';
+if isfield(Param,'UvmatRevision')
+    Data.Program=[Data.Program ', uvmat r' Param.UvmatRevision];
+end
 Data.CivStage=0;%default
 list_param=(fieldnames(Param.ActionInput.Civ1))';
-    list_param(strcmp('TestCiv1',list_param))=[];% remove the parameter TestCiv1 from the list
-    Civ1_param=regexprep(list_param,'^.+','Civ1_$0');% insert 'Civ1_' before  each string in list_param
-    Civ1_param=[{'Civ1_ImageA','Civ1_ImageB','Civ1_Time','Civ1_Dt'} Civ1_param]; %insert the names of the two input images
+list_param(strcmp('TestCiv1',list_param))=[];% remove the parameter TestCiv1 from the list
+Civ1_param=regexprep(list_param,'^.+','Civ1_$0');% insert 'Civ1_' before  each string in list_param
+Civ1_param=[{'Civ1_ImageA','Civ1_ImageB','Civ1_Time','Civ1_Dt'} Civ1_param]; %insert the names of the two input images
 Data.ListGlobalAttribute=[ListGlobalAttribute Civ1_param];
-            % set the list of variables
-        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
-        Data.VarDimName={'nb_vec_1','nb_vec_1','nb_vec_1','nb_vec_1','nb_vec_1','nb_vec_1','nb_vec_1'};
-        Data.VarAttribute{1}.Role='coord_x';
-        Data.VarAttribute{2}.Role='coord_y';
-        Data.VarAttribute{3}.Role='vector_x';
-        Data.VarAttribute{4}.Role='vector_y';
-        Data.VarAttribute{5}.Role='vector_z';
-        Data.VarAttribute{6}.Role='ancillary';
-        Data.VarAttribute{7}.Role='errorflag';
+% set the list of variables
+Data.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
+Data.VarDimName={'npz','npy','npx',{'npz','npy','npx'},{'npz','npy','npx'},{'npz','npy','npx'},{'npz','npy','npx'},...
+    {'npz','npy','npx'},{'npz','npy','npx'},{'npz','npy','npx'}};
+Data.VarAttribute{1}.Role='coord_x';
+Data.VarAttribute{2}.Role='coord_y';
+Data.VarAttribute{3}.Role='vector_x';
+Data.VarAttribute{4}.Role='vector_y';
+Data.VarAttribute{5}.Role='vector_z';
+Data.VarAttribute{6}.Role='ancillary';
+Data.VarAttribute{7}.Role='errorflag';
+Data.Coord_z=j1_series_Civ1;
+% path for output nc file
+OutputPath=fullfile(Param.OutputPath,num2str(Param.Experiment),num2str(Param.Device),[Param.OutputSubDir Param.OutputDirExt]);
 
 %% get timing from the ImaDoc file or input video
@@ -227 +234 @@
     par_civ1.ImageA=zeros(2*SearchRange_z+1,npy,npx);
 par_civ1.ImageB=zeros(2*SearchRange_z+1,npy,npx);
+
+
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -268 +277 @@
     
     Data.CivStage=1;
-    % output nc file
-    ncfile_out=fullfile_uvmat(OutputPath,Param.InputTable{1,2},Param.InputTable{1,3},Param.InputTable{1,5},...
-        NomTypeNc,i1_series_Civ1(ifield),i2_series_Civ1(ifield),j1_series_Civ1(ifield),j2_series_Civ1(ifield));
+  
+    ncfile_out=fullfile_uvmat(OutputPath,'',Param.InputTable{1,3},'.nc',...
+                '_1-1',i1_series_Civ1(ifield),i2_series_Civ1(ifield)); %output file name
 
     if ~CheckOverwrite && exist(ncfile_out,'file')
@@ -284 +293 @@
         disp('civ1 started')
 
-
-        % read input images (except in mode Test where it is introduced directly in Param.ActionInput.Civ1.ImageNameA and B)
-
-        par_civ1.ImageA=zeros(2*SearchRange_z+1,npy,npx);
+        % read input images by vertical blocks with nbre of images equal to 2*SearchRange+1,
+        par_civ1.ImageA=zeros(2*SearchRange_z+1,npy,npx);%image block initiation
         par_civ1.ImageB=zeros(2*SearchRange_z+1,npy,npx);
-        %first vertical block centered at image index islice=par_civ1.Dz
-        islice=par_civ1.Dz;
-         for iz=1:par_civ1.Dz+SearchRange_z
-              ImageName_A=fullfile_uvmat(RootPath_A,SubDir_A,RootFile_A,FileExt_A,NomType_A,i1_series_Civ1(1,ifield),[],j1_series_Civ1(iz,1));%
+        
+        z_index=1;%first vertical block centered at image index z_index=par_civ1.Dz
+        for iz=1:par_civ1.Dz+SearchRange_z
+            j_image_index=j1_series_Civ1(iz,1)% j index of the current image
+            ImageName_A=fullfile_uvmat(RootPath_A,SubDir_A,RootFile_A,FileExt_A,NomType_A,i1_series_Civ1(1,ifield),[],j_image_index);%
             A= read_image(ImageName_A,FileType_A);
-            ImageName_B=fullfile_uvmat(RootPath_B,SubDir_B,RootFile_B,FileExt_B,NomType_B,i2_series_Civ1(1,ifield),[],j1_series_Civ1(iz,1));
+            ImageName_B=fullfile_uvmat(RootPath_B,SubDir_B,RootFile_B,FileExt_B,NomType_B,i2_series_Civ1(1,ifield),[],j_image_index);
             B= read_image(ImageName_B,FileType_B);
             par_civ1.ImageA(iz+par_civ1.Dz-1,:,:) = A;
             par_civ1.ImageB(iz+par_civ1.Dz-1,:,:) = B;
-         end
-         % caluclate velocity data (y and v in indices, reverse to y component)
-            [Data.Civ1_X(islice,:,:),Data.Civ1_Y(islice,:,:), utable, vtable,wtable, ctable, FF, result_conv, errormsg] = civ3D (par_civ1);
+        end
+        % calculate velocity data at the first position z, corresponding to j_index=par_civ1.Dz
+        [Data.Civ1_X(1,:,:),Data.Civ1_Y(1,:,:), Data.Civ1_U(1,:,:), Data.Civ1_V(1,:,:),Data.Civ1_W(1,:,:),...
+            Data.Civ1_C(1,:,:), Data.Civ1_FF(1,:,:), result_conv, errormsg] = civ3D (par_civ1);
+        if ~isempty(errormsg)
+            disp_uvmat('ERROR',errormsg,checkrun)
+            return
+        end
+        for z_index=2:floor((NbSlice-1)/par_civ1.Dz)% loop on slices
+            par_civ1.ImageA=circshift(par_civ1.ImageA,-par_civ1.Dz);%shift the indices in the image block upward by par_civ1.Dz
+            par_civ1.ImageB=circshift(par_civ1.ImageA,-par_civ1.Dz);
+            for iz=1:par_civ1.Dz %read the new images at the end of the image block
+                j_image_index=j1_series_Civ1(z_index*par_civ1.Dz+SearchRange_z-par_civ1.Dz+iz,1)
+                ImageName_A=fullfile_uvmat(RootPath_A,SubDir_A,RootFile_A,FileExt_A,NomType_A,i1_series_Civ1(1,ifield),[],j_image_index);%
+                A= read_image(ImageName_A,FileType_A);
+                ImageName_B=fullfile_uvmat(RootPath_B,SubDir_B,RootFile_B,FileExt_B,NomType_B,i2_series_Civ1(1,ifield),[],j_image_index);
+                B= read_image(ImageName_B,FileType_B);
+                par_civ1.ImageA(2*SearchRange_z+1-par_civ1.Dz+iz,:,:) = A;
+                par_civ1.ImageB(2*SearchRange_z+1-par_civ1.Dz+iz,:,:) = B;
+            end
+            % caluclate velocity data (y and v in indices, reverse to y component)
+            [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,:,:),...
+                Data.Civ1_C(z_index,:,:), Data.Civ1_FF(z_index,:,:), result_conv, errormsg] = civ3D (par_civ1);
+
             if ~isempty(errormsg)
                 disp_uvmat('ERROR',errormsg,checkrun)
                 return
             end
-        for islice=2*par_civ1.Dz:NbSlice% loop on slices for the first image in the pair 
-            par_civ1.ImageA=circshift(par_civ1.ImageA,-par_civ1.Dz);%shift the indces in the block upward by par_civ1.Dz
-            par_civ1.ImageB=circshift(par_civ1.ImageA,-par_civ1.Dz);
-              for iz=1:par_civ1.Dz
-            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));%
-            A= read_image(ImageName_A,FileType_A);
-            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));
-            B= read_image(ImageName_B,FileType_B);
-            par_civ1.ImageA(2*SearchRange_z+1-par_civ1.Dz+iz,:,:) = A;
-            par_civ1.ImageB(2*SearchRange_z+1-par_civ1.Dz+iz,:,:) = B;
-            end
-            % caluclate velocity data (y and v in indices, reverse to y component)
-            [Data.Civ1_X(islice,:,:),Data.Civ1_Y(islice,:,:), utable, vtable,wtable, ctable, FF, result_conv, errormsg] = civ3D (par_civ1);
-            if ~isempty(errormsg)
-                disp_uvmat('ERROR',errormsg,checkrun)
-                return
-            end
-
-        end
+        end
+        Data.Civ1_C=double(Data.Civ1_C);
+        Data.Civ1_C(Data.Civ1_C==Inf)=0;
+        [npz,npy,npx]=size(Data.Civ1_X);
+        Data.Coord_z=1:npz;
+        Data.Coord_y=1:npy;
+        Data.Coord_x=1:npx;
         % case of mask TO ADAPT
         if par_civ1.CheckMask&&~isempty(par_civ1.Mask)
@@ -668 +685 @@
         end
         Data.ListGlobalAttribute=[Data.ListGlobalAttribute Fix2_param];
-        Data.Civ2_FF=double(detect_false(Param.ActionInput.Fix2,Data.Civ2_C,Data.Civ2_U,Data.Civ2_V,Data.Civ2_FF));
+        Data.Civ2_FF=detect_false(Param.ActionInput.Fix2,Data.Civ2_C,Data.Civ2_U,Data.Civ2_V,Data.Civ2_FF);
         Data.CivStage=Data.CivStage+1;
     end
@@ -728 +745 @@
     time_total=toc(tstart);
     disp(['ellapsed time ' num2str(time_total/60,2) ' minutes'])
-    disp(['time image reading ' num2str(time_input,2) ' s'])
     disp(['time civ1 ' num2str(time_civ1,2) ' s'])
     disp(['time patch1 ' num2str(time_patch1,2) ' s'])
     disp(['time civ2 ' num2str(time_civ2,2) ' s'])
     disp(['time patch2 ' num2str(time_patch2,2) ' s'])
-    disp(['time other ' num2str((time_total-time_input-time_civ1-time_patch1-time_civ2-time_patch2),2) ' s'])
-    % end
+
 end
 
@@ -860 +875 @@
     for ivec_x=1:nbvec_x
         for ivec_y=1:nbvec_y
-            ivec_y
-            iref=round(par_civ.Grid(ivec_y,ivec_x,1)+0.5)% xindex on the image A for the middle of the correlation box
-            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
+            iref=round(par_civ.Grid(ivec_y,ivec_x,1)+0.5);% xindex on the image A for the middle of the correlation box
+            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
             subrange1_x=iref-ibx2:iref+ibx2;% x indices defining the first subimage
             subrange1_y=jref-iby2:jref+iby2;% y indices defining the first subimage
@@ -918 +932 @@
 
                 
-                    %reference: Oliver Pust, PIV: Direct Cross-Correlation
+                    
                     for kz=1:par_civ.SearchBoxSize(3)
                         subima2=squeeze(image2_crop(kz,:,:));
@@ -928 +942 @@
                
                     end
-                    [corrmax,z]=max(max_xy);
+                    [corrmax,zmax]=max(max_xy);
                 
-                    x=xk(z);
-                    y=yk(z);
+                    x=xk(zmax);
+                    y=yk(zmax);
                     result_conv=(result_conv/corrmax)*255; %normalize, peak=always 255
-                    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
-                    sum_square=sum(sum(squeeze(image1_crop(z,:,:).*image1_crop(z,:,:))));
+                    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
+                    sum_square=sum(sum(squeeze(image1_crop(zmax,:,:).*image1_crop(zmax,:,:))));
                     sum_square=sum_square*sum(sum(subimage2_crop.*subimage2_crop));% product of variances of image 1 and 2
                     sum_square=sqrt(sum_square);% srt of the variance product to normalise correlation
@@ -940 +954 @@
        
                             if par_civ.CorrSmooth==1
-                                [vector,FF(ivec_y,ivec_x)] = SUBPIXGAUSS (result_conv(z,:,:),x,y);%TODO: improve by max optimisation along z
+                                [vector,FF(ivec_y,ivec_x)] = SUBPIXGAUSS (squeeze(result_conv(zmax,:,:)),x,y);%TODO: improve by max optimisation along z
                             elseif par_civ.CorrSmooth==2
-                                [vector,FF(ivec_y,ivec_x)] = SUBPIX2DGAUSS (result_conv(z,:,:),x,y);
+                                [vector,FF(ivec_y,ivec_x)] = SUBPIX2DGAUSS (squeeze(result_conv(zmax,:,:)),x,y);
                             else
-                                [vector,FF(ivec_y,ivec_x)] = quadr_fit(result_conv(z,:,:),x,y);
+                                [vector,FF(ivec_y,ivec_x)] = quadr_fit(squeeze(result_conv(zmax,:,:)),x,y);
                             end
                             utable(ivec_y,ivec_x)=vector(1)*mesh+shiftx(ivec_y,ivec_x);
@@ -952 +966 @@
                             iref=round(xtable(ivec_y,ivec_x)+0.5);% nearest image index for the middle of the vector
                             jref=round(ytable(ivec_y,ivec_x)+0.5);
-                            wtable(ivec_y,ivec_x)=z-kref;
+                            wtable(ivec_y,ivec_x)=zmax-kref;
                             % eliminate vectors located in the mask
                             if  checkmask && (iref<1 || jref<1 ||iref>npx_ima || jref>npy_ima ||( par_civ.Mask(jref,iref)<200 && par_civ.Mask(jref,iref)>=100))
@@ -962 +976 @@
 
                     else
-                        FF(ivec_y,ivec_x)=1;
+                        FF(ivec_y,ivec_x)=true;
                     end
                 end
@@ -975 +989 @@
 % OUPUT:
 % vector = optimum displacement vector with subpixel correction
-% F =flag: =0 OK
-%           =-2 , warning: max too close to the edge of the search box (1 pixel margin)
+% FF =flag: ='true' max too close to the edge of the search box (1 pixel margin)
 % INPUT:
 % x,y: position of the maximum correlation at integer values
 
-function [vector,F] = SUBPIXGAUSS (result_conv,x,y)
+function [vector,FF] = SUBPIXGAUSS (result_conv,x,y)
 %------------------------------------------------------------------------
-% vector=[0 0]; %default
-F=0;
+vector=[0 0]; %default
 [npy,npx]=size(result_conv);
 result_conv(result_conv<1)=1; %set to 1 correlation values smaller than 1  (=0 by discretisation, to avoid divergence in the log)
 %the following 8 lines are copyright (c) 1998, Uri Shavit, Roi Gurka, Alex Liberzon, Technion ??? Israel Institute of Technology
 %http://urapiv.wordpress.com
-peaky = y;
-if y < npy && y > 1
-    f0 = log(result_conv(y,x));
-    f1 = log(result_conv(y-1,x));
-    f2 = log(result_conv(y+1,x));
-    peaky = peaky+ (f1-f2)/(2*f1-4*f0+2*f2);
-else
-    F=1; % warning flag for vector truncated by the limited search box
-end
-peakx=x;
-if x < npx-1 && x > 1
-    f0 = log(result_conv(y,x));
-    f1 = log(result_conv(y,x-1));
-    f2 = log(result_conv(y,x+1));
-    peakx = peakx+ (f1-f2)/(2*f1-4*f0+2*f2);
-else
-    F=1; % warning flag for vector truncated by the limited search box
-end
-vector=[peakx-floor(npx/2)-1 peaky-floor(npy/2)-1];
+FF=false;
+    peaky = y;
+    if y < npy && y > 1
+        f0 = log(result_conv(y,x));
+        f1 = log(result_conv(y-1,x));
+        f2 = log(result_conv(y+1,x));
+        peaky = peaky+ (f1-f2)/(2*f1-4*f0+2*f2);
+    else
+        FF=true; % error flag for vector truncated by the limited search box in y
+    end
+    peakx=x;
+    if x < npx-1 && x > 1
+        f0 = log(result_conv(y,x));
+        f1 = log(result_conv(y,x-1));
+        f2 = log(result_conv(y,x+1));
+        peakx = peakx+ (f1-f2)/(2*f1-4*f0+2*f2);
+    else
+        FF=true; % warning flag for vector truncated by the limited search box in x
+    end
+    vector=[peakx-floor(npx/2)-1 peaky-floor(npy/2)-1];
 
 %------------------------------------------------------------------------

trunk/src/series/civ_series.m

@@ -301 +301 @@
     end
     if CheckInputFile
-        OutputPath=fullfile(Param.OutputPath,num2str(Param.Experiment),num2str(Param.Device));
+        OutputPath=fullfile(Param.OutputPath,Param.Experiment,Param.Device);
         if iview_A==0 % no nc file has been entered
             ncfile=fullfile_uvmat(OutputPath,Param.InputTable{1,2},Param.InputTable{1,3},Param.InputTable{1,5},...

trunk/src/series/sliding_average.m

@@ -1 @@
-%'sliding_average_natalya'
+%'sliding_average': calculate the time correlation function at each point
 %------------------------------------------------------------------------
-% function ParamOut=sliding_average(Param)
+% function ParamOut=turb_correlation_time(Param)
 %
 %%%%%%%%%%% GENERAL TO ALL SERIES ACTION FCTS %%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -39 +39 @@
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %=======================================================================
-% Copyright 2008-2024, LEGI UMR 5519 / CNRS UGA G-INP, Grenoble, France
+% Copyright 2008-2022, LEGI UMR 5519 / CNRS UGA G-INP, Grenoble, France
 %   http://www.legi.grenoble-inp.fr
-%   Joel.Sommeria - Joel.Sommeria (A) univ-grenoble-alpes.fr
+%   Joel.Sommeria - Joel.Sommeria (A) legi.cnrs.fr
 %
 %     This file is part of the toolbox UVMAT.
@@ -65 +65 @@
     ParamOut.VelType='off';% menu for selecting the velocity type (options 'off'/'one'/'two',  'off' by default)
     ParamOut.FieldName='one';% menu for selecting the field (s) in the input file(options 'off'/'one'/'two', 'off' by default)
-    ParamOut.FieldTransform = 'on';%can use a transform function
+    ParamOut.FieldTransform = 'off';%can use a transform function
     ParamOut.ProjObject='off';%can use projection object39(option 'off'/'on',
     ParamOut.Mask='off';%can use mask option   (option 'off'/'on', 'off' by default)
     ParamOut.OutputDirExt='.tfilter';%set the output dir extension
-    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
-%     filecell=get_file_series(Param);%check existence of the first input file
-%     if ~exist(filecell{1,1},'file')
-%         msgbox_uvmat('WARNING','the first input file does not exist')
-%     end
+    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
+   
+    %% check the first input file in the series
+     first_j=[];% 
+    if isfield(Param.IndexRange,'first_j'); first_j=Param.IndexRange.first_j; end
+    PairString='';
+    if isfield(Param.IndexRange,'PairString'); PairString=Param.IndexRange.PairString; end
+    [i1,i2,j1,j2] = get_file_index(Param.IndexRange.first_i,first_j,PairString);
+    FirstFileName=fullfile_uvmat(Param.InputTable{1,1},Param.InputTable{1,2},Param.InputTable{1,3},...
+        Param.InputTable{1,5},Param.InputTable{1,4},i1,i2,j1,j2); 
+    if exist(FirstFileName,'file')
+        FileInfo=get_file_info(FirstFileName);
+        if ~(isfield(FileInfo,'FileType')&& strcmp(FileInfo.FileType,'netcdf'))
+            msgbox_uvmat('ERROR','this fct works only on netcdf files (fields projected on a grid, notraw civ data)')
+        return
+        end
+    else
+        msgbox_uvmat('ERROR',['the input file ' FirstFileName ' does not exist'])
+        return
+    end
+
+    %% setting of intput parameters 
+%     answer = msgbox_uvmat('INPUT_MENU','select the variables to filter',ListVarName');
+%     ParamOut.ListVarName=ListVarName(answer);
+    answer = msgbox_uvmat('INPUT_TXT','set the filering window length','3');
+    ParamOut.WindowLength=str2double(answer);
+    
+%         [ParamOut.ActionInput,errormsg] = set_param_input(ListParam,DefaultValue,ParamIn);
+%         if ~isempty(errormsg)
+%             msgbox_uvmat('ERROR',errormsg)
+%         end
+%     [ParamOut,errormsg] = set_param_input(ListParam,DefaultValue,ParamIn);
     return
 end
@@ -184 +211 @@
 nbmissing=0;
 
-%% initialisation manip Coriolis
-char_index=regexp(SubDir{1},'waves_L1_');
-switch(SubDir{1}(char_index+9))
-    case '1'
-        amplitude=2.5 %oscillation amplitude
-        T=24.46;
-        t0=3 ;% dt=0.5 s, torus at its max x at the beginning of motion, i0=7
-    case '2'
-        amplitude=5 %oscillation amplitude
-        T=24.47;
-        t0=8.5; % dt=1/3 s -> image index of starting motion = 26, % torus at its max x at the beginning of motion
-    case '3'
-        amplitude=10 %oscillation amplitude
-        T=24.45;
-        t0=6.5-T/2;% dt=0.25, torus at its minimum x at the beginning of motion
-    case '4'  
-        amplitude=15 %oscillation amplitude
-        T=24.48;
-        t0=3.4;     %dt=0.2 -> i0=18 image index of starting motion, % torus at its max x at the beginning of motion
-end
-%NbPeriod=2; %number of periods for the sliding average
-NbPeriod=4; %number of periods for the sliding average
-omega=2*2*pi/T;%harmonic
-Lscale=15;%diameter of the torus, length scale for normalisation
-Uscale=amplitude*omega;
-
-DataOut.ListGlobalAttribute= {'Conventions','Time'};
+%% initialisation output data
+
+DataOut.ListGlobalAttribute= {'Conventions'};
 DataOut.Conventions='uvmat';
-DataOut.ListVarName={'coord_y','coord_x','Umean','Vmean','Ucos','Vcos','Usin','Vsin','DUDXsin','DUDXcos','DUDYsin','DVDXsin','DVDXcos'...
-    ,'DVDYsin','Ustokes','Vstokes'};
-DataOut.VarDimName={'coord_y','coord_x',{'coord_y','coord_x'},{'coord_y','coord_x'},{'coord_y','coord_x'},{'coord_y','coord_x'},...
-    {'coord_y','coord_x'},{'coord_y','coord_x'},{'coord_y','coord_x'},{'coord_y','coord_x'},{'coord_y','coord_x'},...
-    {'coord_y','coord_x'},{'coord_y','coord_x'},{'coord_y','coord_x'},{'coord_y','coord_x'},{'coord_y','coord_x'}};
+DataOut.ListVarName={'Time','coord_y','coord_x','Ufilter','Vfilter'};
+DataOut.VarDimName={'Time','coord_y','coord_x',{'Time','coord_y','coord_x'},{'Time','coord_y','coord_x'}};
 
 %%%%%%%%%%%%%%%% loop on field indices %%%%%%%%%%%%%%%%
@@ -226 +226 @@
     return
 end
-[Data,tild,errormsg]=nc2struct(filecell{1,end});    
-Time_end=Data.Time;
-dt=(Time_end-Time_1)/(NbField-1); %time interval 
-NpTime=round(NbPeriod*T/dt+1);
+[Data,tild,errormsg]=nc2struct(filecell{1,2});    
+Time_2=Data.Time;
+dt=(Time_2-Time_1); %time interval 
+NpTime=21; %Nbre de champ pour la moyenne glissante (choisir impair)
 
 OutputPath=fullfile(Param.OutputPath,num2str(Param.Experiment),num2str(Param.Device));
@@ -246 +246 @@
     %%%%%%%%%%% MAIN RUNNING OPERATIONS  %%%%%%%%%%%%
     if index==1 %first field
-        DataOut.coord_x=Field.coord_x/Lscale;
-        DataOut.coord_y=Field.coord_y/Lscale;
+        DataOut.coord_x=Field.coord_x;
+        DataOut.coord_y=Field.coord_y;
         npy=numel(DataOut.coord_y);
         npx=numel(DataOut.coord_x);
-        Umean=zeros(NpTime,npy,npx);
-        Vmean=zeros(NpTime,npy,npx);
-        Ucos=zeros(NpTime,npy,npx);
-        Vcos=zeros(NpTime,npy,npx);
-        Usin=zeros(NpTime,npy,npx);
-        Vsin=zeros(NpTime,npy,npx);
-        DUDXcos=zeros(NpTime,npy,npx);
-        DUDXsin=zeros(NpTime,npy,npx);
-        DUDYsin=zeros(NpTime,npy,npx);
-        DVDXcos=zeros(NpTime,npy,npx);
-        DVDXsin=zeros(NpTime,npy,npx);
-        DVDYsin=zeros(NpTime,npy,npx);
-    end
-    Time(index)=Field.Time-t0;%time from the start of the motion
-    Umean=circshift(Umean,[-1 0 0]); %shift U by ishift along the first index
-    Vmean=circshift(Vmean,[-1 0 0]); %shift U by ishift along the first index
-    Ucos=circshift(Ucos,[-1 0 0]); %shift U by ishift along the first index
-    Vcos=circshift(Vcos,[-1 0 0]); %shift U by ishift along the first index
-    Usin=circshift(Usin,[-1 0 0]); %shift U by ishift along the first index
-    Vsin=circshift(Vsin,[-1 0 0]); %shift U by ishift along the first index
-    DUDXcos=circshift(DUDXcos,[-1 0 0]);
-    DUDXsin=circshift(DUDXsin,[-1 0 0]);
-    DUDYsin=circshift(DUDYsin,[-1 0 0]);        
-    DVDXcos=circshift(DVDXcos,[-1 0 0]);
-    DVDXsin=circshift(DVDXsin,[-1 0 0]);
-    DVDYsin=circshift(DVDYsin,[-1 0 0]);       
-    Umean(end,:,:)=Field.U;
-    Vmean(end,:,:)=Field.V;
-    Ucos(end,:,:)=Field.U*cos(omega*Time(index));
-    Vcos(end,:,:)=Field.V*cos(omega*Time(index));
-    Usin(end,:,:)=Field.U*sin(omega*Time(index));
-    Vsin(end,:,:)=Field.V*sin(omega*Time(index));
-    DUDXcos(end,:,:)=Field.DUDX*cos(omega*Time(index));
-    DUDXsin(end,:,:)=Field.DUDX*sin(omega*Time(index));
-    DUDYsin(end,:,:)=Field.DUDY*sin(omega*Time(index));% ParamOut=[];%default output
-
-    DVDXcos(end,:,:)=Field.DVDX*cos(omega*Time(index));
-    DVDXsin(end,:,:)=Field.DVDX*sin(omega*Time(index));
-    DVDYsin(end,:,:)=Field.DVDY*sin(omega*Time(index));
-    DataOut.Time=(Time(index)-(NpTime-1)*dt/2)/T;%time inperiods from the beginning of the oscillation (torus at max abscissa)
-    DataOut.Umean=(1/Uscale)*squeeze(nanmean(Umean,1));
-    DataOut.Vmean=(1/Uscale)*squeeze(nanmean(Vmean,1));
-    DataOut.Ucos=2*(1/Uscale)*squeeze(nanmean(Ucos,1));
-    DataOut.Vcos=2*(1/Uscale)*squeeze(nanmean(Vcos,1));
-    DataOut.Usin=2*(1/Uscale)*squeeze(nanmean(Usin,1));
-    DataOut.Vsin=2*(1/Uscale)*squeeze(nanmean(Vsin,1));
-    DataOut.DUDXcos=2*squeeze(nanmean(DUDXcos,1));
-    DataOut.DUDXsin=2*squeeze(nanmean(DUDXsin,1));
-    DataOut.DUDYsin=2*squeeze(nanmean(DUDYsin,1));
-    DataOut.DVDXcos=2*squeeze(nanmean(DVDXcos,1));
-    DataOut.DVDXsin=2*squeeze(nanmean(DVDXsin,1));
-    DataOut.DVDYsin=2*squeeze(nanmean(DVDYsin,1));
-    DataOut.Ustokes=(1/omega)*(1/Uscale)*(DataOut.Ucos.*DataOut.DUDXsin+DataOut.Vcos.*DataOut.DUDYsin);
-    DataOut.Vstokes=(1/omega)*(1/Uscale)*(DataOut.Ucos.*DataOut.DVDXsin+DataOut.Vcos.*DataOut.DVDYsin);
+        Ufilter=zeros(NpTime,npy,npx);
+        Vfilter=zeros(NpTime,npy,npx);
+    end
+    Time(index)=Field.Time;%time 
+    Ufilter=circshift(Ufilter,[-1 0 0]); %shift U by ishift along the first index
+    Vfilter=circshift(Vfilter,[-1 0 0]); %shift U by ishift along the first index    
+    Ufilter(end,:,:)=Field.U;
+    Vfilter(end,:,:)=Field.V;
+    DataOut.Time=(Time(index)-(NpTime-1)*dt/2);%mid time 
+    DataOut.Ufilter=squeeze(nanmean(Ufilter,1));
+    DataOut.Vfilter=squeeze(nanmean(Vfilter,1));
+    
 
     % writing the result file as netcdf file
-    i1=i1_series{1}(index);
+    i1=i1_series{1}(index)-ceil(NpTime/2);
     OutputFile=fullfile_uvmat(OutputPath,OutputDir,RootFileOut,'.nc',NomTypeOut,i1);
     errormsg=struct2nc(OutputFile, DataOut);

trunk/src/series/test_filter_tps.m

@@ -72 +72 @@
     % check the existence of the first file in the series
      first_j=[];% note that the function will propose to cover the whole range of indices
-    if isfield(Param.IndexRange,'MinIndex_j'); first_j=Param.IndexRange.MinIndex_j; end
-    last_j=[];
-    if isfield(Param.IndexRange,'MaxIndex_j'); last_j=Param.IndexRange.MaxIndex_j; end
+    if isfield(Param.IndexRange,'first_j'); first_j=Param.IndexRange.first_j; end
     PairString='';
     if isfield(Param.IndexRange,'PairString'); PairString=Param.IndexRange.PairString; end
@@ -81 +79 @@
         Param.InputTable{1,5},Param.InputTable{1,4},i1,i2,j1,j2);
     if ~exist(FirstFileName,'file')
-        msgbox_uvmat('WARNING',['the first input file ' FirstFileName ' does not exist'])
+        msgbox_uvmat('ERROR',['the input file ' FirstFileName ' does not exist'])
         return
-    else
-        [i1,i2,j1,j2] = get_file_index(Param.IndexRange.last_i,last_j,PairString);
-        LastFileName=fullfile_uvmat(Param.InputTable{1,1},Param.InputTable{1,2},Param.InputTable{1,3},...
-        Param.InputTable{1,5},Param.InputTable{1,4},i1,i2,j1,j2);
-        if ~exist(FirstFileName,'file')
-             msgbox_uvmat('WARNING',['the last input file ' LastFileName ' does not exist'])
-        end
     end
 
@@ -114 +105 @@
         end
     else
-        msgbox_uvmat('ERROR',['invalid file type input: ' FileType ' not a civ data'])
+        msgbox_uvmat('ERROR','invalid file type input: test_filter_tps proceeds raw civ data')
         return
     end

trunk/src/set_param_input.m

@@ -9 +9 @@
 
 function [ParamOut,errormsg] = set_param_input(ListParam,DefaultValue,ParamIn,Comment)
+ParamOut=[];
 errormsg=[];
 NbParam=numel(ListParam);
@@ -32 +33 @@
 end
 dlg_title = 'get the input parameters';
-answer = inputdlg(ListParam,dlg_title,NbParam,prompt);
+options.Resize='on';
+answer = inputdlg(ListParam,dlg_title,NbParam,prompt,options);
+%answer = msgbox_uvmat('INPUT_TXT',ListParam);
 if isempty(answer)
     return

trunk/src/struct2nc.m

@@ -38 +38 @@
 function [errormsg,nc]=struct2nc(flname,Data,action,ListDimName,DimValue,VarDimIndex)
 nc=[];
+errormsg='';
 if ~ischar(flname)
     errormsg='invalid input for the netcf file name';
@@ -81 +82 @@
             end
             if ~testvar               
-                eval(['cte=Data.' keys{iattr} ';'])
+                cte=Data.(keys{iattr});
                 if (ischar(cte) ||isnumeric(cte)) &&  ~isempty(cte)%&& ~isequal(cte,'')
                     %write constant only if it is numeric or char string, and not empty
@@ -106 +107 @@
 %% create the variables
 varid=nan(1,length(Data.ListVarName));
+VarClass=cell(1,length(ListVarName));
 for ivar=1:length(ListVarName)
     if isfield(Data,ListVarName{ivar})
@@ -141 +143 @@
 netcdf.endDef(nc); %put in data mode
 
-%% fill the variables with input data
+%% fill the variables with input data except in mode 'keep_open' (variables will be filled later)
 if ~(exist('action','var') && strcmp(action,'keep_open'))
     for ivar=1:length(ListVarName)
@@ -169 +171 @@
     end
 end
+if ~ (exist('action','var') && strcmp(action,'keep_open'))
  netcdf.close(nc)
 [success,errormsg] = fileattrib(flname ,'+w');% allow writing access for the group of users
+end
 %'check_field_structure': check the validity of the field struture representation consistant with the netcdf format
 %------------------------------------------------------------------------