Changeset 1163 for trunk

Timestamp:

Jul 25, 2024, 8:36:17 PM (5 months ago)

Author:

sommeria

Message:

SearchBoxSize? changed to SearchRange? in civ, possibility of a grid set by a netcdf file

Location:

trunk/src

Files:

• civ.m (added)
• filter_tps.m (modified) (1 diff)
• filter_tps_3D.m (modified) (2 diffs)
• mouse_motion.m (modified) (1 diff)
• script_makegrid.m (added)
• series/civ_3D.m (modified) (9 diffs)
• series/civ_input.fig (modified) (previous)
• series/civ_input.m (modified) (26 diffs)
• series/civ_series.m (modified) (18 diffs)
• series/extract_multitif.m (modified) (2 diffs)
• series/extract_rdvision.m (modified) (6 diffs)
• series/test_civ.m (added)
• uvmat.m (modified) (1 diff)

trunk/src/filter_tps.m

@@ -93 +93 @@
         ind_sel= find(~FF & Coord(:,1)>=SubRange(1,1,isub) & Coord(:,1)<=SubRange(1,2,isub) & Coord(:,2)>=SubRange(2,1,isub) & Coord(:,2)<=SubRange(2,2,isub));% indices of vectors in the subdomain #isub
          check_partial_domain=sum(SubRange(:,1,isub)> MinCoord' | SubRange(:,2,isub)< MaxCoord');
-        % isub
-        % numel(ind_sel)
 
         % if no vector in the subdomain  #isub, skip the subdomain

trunk/src/filter_tps_3D.m

@@ -66 +66 @@
 CentreX=reshape(CentreX,1,[]);% X positions of subdomain centres
 CentreY=reshape(CentreY,1,[]);% Y positions of subdomain centres
-CentreZ=reshape(CentreZ,1,[]);% Y positions of subdomain centres
+CentreZ=reshape(CentreZ,1,[]);% Z positions of subdomain centres
 
 %% smoothing parameter: CHANGED 03 May 2024 TO GET RESULTS INDEPENDENT OF SUBDOMAINSIZE
@@ -92 +92 @@
     SubRange(1,:,isub)=[CentreX(isub)-0.55*Siz(1) CentreX(isub)+0.55*Siz(1)];%bounds of subdomain #isub in x coordinate
     SubRange(2,:,isub)=[CentreY(isub)-0.55*Siz(2) CentreY(isub)+0.55*Siz(2)];%bounds of subdomain #isub in y coordinate
+    SubRange(3,:,isub)=[CentreZ(isub)-0.55*Siz(3) CentreZ(isub)+0.55*Siz(3)];%bounds of subdomain #isub in y coordinate
     ind_sel=0;%initialize set of vector indices in the subdomain
     %increase iteratively the subdomain if it contains less than

trunk/src/mouse_motion.m

@@ -265 +265 @@
                 ibx2=floor((par_civ.CorrBoxSize(1)-1)/2);
                 iby2=floor((par_civ.CorrBoxSize(2)-1)/2);
-                isx2=floor((par_civ.SearchBoxSize(1)-1)/2);
-                isy2=floor((par_civ.SearchBoxSize(2)-1)/2);
+                isx2=par_civ.SearchRange(1);
+                isy2=par_civ.SearchRange(2);
                 hhh=findobj(CurrentAxes,'Tag','PIV_box_marker');
                 hhhh=findobj(CurrentAxes,'Tag','PIV_search_marker');

trunk/src/series/civ_3D.m

@@ -190 +190 @@
 Data.VarAttribute{7}.Role='ancillary';
 Data.VarAttribute{8}.Role='errorflag';
-Data.Coord_z=j1_series_Civ1;
+% 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]);
@@ -230 +230 @@
     par_civ1.ImageHeight=FileInfo.Height;npy=FileInfo.Height;
     par_civ1.ImageWidth=FileInfo.Width;npx=FileInfo.Width;
-SearchRange_z=floor(Param.ActionInput.Civ1.SearchBoxSize(3)/2);
+SearchRange_z=Param.ActionInput.Civ1.SearchRange(3);
     par_civ1.Dz=Param.ActionInput.Civ1.Dz;
     par_civ1.ImageA=zeros(2*SearchRange_z+1,npy,npx);
@@ -296 +296 @@
         par_civ1.ImageA=zeros(2*SearchRange_z+1,npy,npx);%image block initiation
         par_civ1.ImageB=zeros(2*SearchRange_z+1,npy,npx);
-        
+        Data.Coord_z=SearchRange_z+1:par_civ1.Dz:NbSlice-1;
         z_index=1;%first vertical block centered at image index z_index=SearchRange_z+1
         for iz=1:2*SearchRange_z+1
@@ -315 +315 @@
         end
      % loop on slices
+
         for z_index=2:floor((NbSlice-SearchRange_z)/par_civ1.Dz)% loop on slices
             par_civ1.ImageA=circshift(par_civ1.ImageA,-par_civ1.Dz,1);%shift the indices in the image block upward by par_civ1.Dz
             par_civ1.ImageB=circshift(par_civ1.ImageB,-par_civ1.Dz,1);
             for iz=1:par_civ1.Dz %read the new images at the end of the image block
+                image_index=z_index*par_civ1.Dz+SearchRange_z-par_civ1.Dz+iz+1;
+                if image_index<=size(j1_series_Civ1,1)
                 j_image_index=j1_series_Civ1(z_index*par_civ1.Dz+SearchRange_z-par_civ1.Dz+iz+1,1)
-                ImageName_A=fullfile_uvmat(RootPath_A,SubDir_A,RootFile_A,FileExt_A,NomType_A,i1_series_Civ1(1,ifield),[],j_image_index);%
+                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);
+                else
+                    A=zeros(1,size(par_civ1.ImageA,2),size(par_civ1.ImageA,3));
+                    B=zeros(1,size(par_civ1.ImageA,2),size(par_civ1.ImageA,3));
+                end
                 par_civ1.ImageA(2*SearchRange_z+1-par_civ1.Dz+iz,:,:) = A;
                 par_civ1.ImageB(2*SearchRange_z+1-par_civ1.Dz+iz,:,:) = B;
@@ -338 +345 @@
         Data.Civ1_Y=npy-Data.Civ1_Y+1;%reverse y
         [npz,npy,npx]=size(Data.Civ1_X);
-        Data.Coord_z=SearchRange_z+1:par_civ1.Dz:NbSlice-1;
+        
       %  Data.Coord_y=flip(1:npy);
        % Data.Coord_x=1:npx;
@@ -422 +429 @@
         % list the variables to record
         nbvar=length(Data.ListVarName);
-        Data.ListVarName=[Data.ListVarName {'Civ1_U_smooth','Civ1_V_smooth','Civ1_SubRange','Civ1_NbCentres','Civ1_Coord_tps','Civ1_U_tps','Civ1_V_tps'}];
-        Data.VarDimName=[Data.VarDimName {'nb_vec_1','nb_vec_1',{'nb_coord','nb_bounds','nb_subdomain_1'},'nb_subdomain_1',...
-            {'nb_tps_1','nb_coord','nb_subdomain_1'},{'nb_tps_1','nb_subdomain_1'},{'nb_tps_1','nb_subdomain_1'}}];
+        Data.ListVarName=[Data.ListVarName {'Civ1_U_smooth','Civ1_V_smooth','Civ1_W_smooth'}];
+        Data.VarDimName=[Data.VarDimName {{'npz','npy','npx'},{'npz','npy','npx'},{'npz','npy','npx'}}];
         Data.VarAttribute{nbvar+1}.Role='vector_x';
         Data.VarAttribute{nbvar+2}.Role='vector_y';
-        Data.VarAttribute{nbvar+5}.Role='coord_tps';
-        Data.VarAttribute{nbvar+6}.Role='vector_x';
-        Data.VarAttribute{nbvar+7}.Role='vector_y';
-        Data.Civ1_U_smooth=Data.Civ1_U; % zeros(size(Data.Civ1_X));
-        Data.Civ1_V_smooth=Data.Civ1_V; %zeros(size(Data.Civ1_X));
+        Data.VarAttribute{nbvar+5}.Role='vector_z';      
+        Data.Civ1_U_smooth=Data.Civ1_U;
+        Data.Civ1_V_smooth=Data.Civ1_V; 
+        Data.Civ1_W_smooth=Data.Civ1_W; 
         if isfield(Data,'Civ1_FF')
             ind_good=find(Data.Civ1_FF==0);
         else
-            ind_good=1:numel(Data.Civ1_X);
+            disp_uvmat('ERROR','detection of error vectors (FIX operation) needed before PATCH' ,checkrun)
+            return
         end
         if isempty(ind_good)
@@ -443 +449 @@
 
         % perform Patch calculation using the UVMAT fct 'filter_tps'
-
-        [Data.Civ1_SubRange,Data.Civ1_NbCentres,Data.Civ1_Coord_tps,Data.Civ1_U_tps,Data.Civ1_V_tps,~,Ures, Vres,~,FFres]=...
-            filter_tps([Data.Civ1_X(ind_good) Data.Civ1_Y(ind_good)],Data.Civ1_U(ind_good),Data.Civ1_V(ind_good),[],Data.Patch1_SubDomainSize,Data.Patch1_FieldSmooth,Data.Patch1_MaxDiff);
-        Data.Civ1_U_smooth(ind_good)=Ures;% take the interpolated (smoothed) velocity values for good vectors, keep civ1 data for the other
-        Data.Civ1_V_smooth(ind_good)=Vres;
+Civ1_Z=0.5*Data.Civ1_W(ind_good);
+for iz=1:numel(Data.Coord_z)
+    Civ1_Z(iz,:,:)=Civ1_Z(iz,:,:)+Data.Coord_z(iz);
+end
+        [Data.Civ1_SubRange,Data.Civ1_NbCentres,Data.Civ1_Coord_tps,Data.Civ1_U_tps,Data.Civ1_V_tps,Data.Civ1_W_tps,...
+            Data.Civ1_U_smooth(ind_good),Data.Civ1_V_smooth(ind_good),Data.Civ1_V_smooth(ind_good),FFres]=...
+            filter_tps_3D([Data.Civ1_X(ind_good) Data.Civ1_Y(ind_good)],Civ_1_Z,Data.Civ1_U(ind_good),Data.Civ1_V(ind_good),Data.Civ1_W(ind_good),Data.Patch1_SubDomainSize,Data.Patch1_FieldSmooth,Data.Patch1_MaxDiff);
         Data.Civ1_FF(ind_good)=uint8(FFres);
         time_patch1=toc(tstart_patch1);
@@ -778 +786 @@
 ibx2=floor(par_civ.CorrBoxSize(1)/2);
 iby2=floor(par_civ.CorrBoxSize(2)/2);
-isx2=floor(par_civ.SearchBoxSize(1)/2);
-isy2=floor(par_civ.SearchBoxSize(2)/2);
-isz2=floor(par_civ.SearchBoxSize(3)/2);
+isx2=par_civ.SearchRange(1);
+isy2=par_civ.SearchRange(2);
+isz2=par_civ.SearchRange(3);
 kref=isz2+1;%middle index of the z slice
 shiftx=round(par_civ.SearchBoxShift(:,1));%use the input shift estimate, rounded to the next integer value
@@ -891 +899 @@
                     % end
 
-                    npxycorr=par_civ.SearchBoxSize(1:2)-par_civ.CorrBoxSize(1:2)+1;
-                    result_conv=zeros([par_civ.SearchBoxSize(3) npxycorr]);%initialise the conv product
-                    max_xy=zeros(par_civ.SearchBoxSize(3),1);%initialise the max correlation vs z
+                    npxycorr=2*par_civ.SearchRange(1:2)+1;
+                    result_conv=zeros([2*par_civ.SearchRange(3)+1 npxycorr]);%initialise the conv product
+                    max_xy=zeros(2*par_civ.SearchRange(3)+1,1);%initialise the max correlation vs z
                     xk=ones(npz,1);%initialise the x displacement corresponding to the max corr 
                     yk=ones(npz,1);%initialise the y displacement corresponding to the max corr 

trunk/src/series/civ_input.m

@@ -79 +79 @@
         set(handles.title_z,'Visible','on')
         set(handles.num_CorrBoxSize_3,'Visible','on')
-        set(handles.num_SearchBoxSize_3,'Visible','on')
+        set(handles.num_SearchRange_3,'Visible','on')
         set(handles.num_SearchBoxShift_3,'Visible','on')
         set(handles.num_Dz,'Visible','on')
@@ -119 +119 @@
             return
         end
-        if isfield(Data,'.Civ1_ImageA')
-        %[PathCiv1_ImageA,Civ1_ImageA,FileExtA]=fileparts(Data.Civ1_ImageA);%look for the source image A
-        %[PathCiv1_ImageB,Civ1_ImageB,FileExtA]=fileparts(Data.Civ1_ImageB);%look for the source image B
-        end
-        if isfield(Data,'Civ2_ImageA')
-            %[PathCiv2_ImageA,Civ2_ImageA,FileExtA]=fileparts(Data.Civ2_ImageA);
-            %[PathCiv2_ImageB,Civ2_ImageB,FileExtA]=fileparts(Data.Civ2_ImageB);
-        end
+        
         if size(Param.InputTable,1)==1
-             if isfield(Data,'.Civ1_ImageA')
-            series('display_file_name',hhseries,Data.Civ1_ImageA,'append');%append the image series to the input list
-                    [~,~,~,~,~,~,~,~,NomTypeImaA]=fileparts_uvmat(Data.Civ1_ImageA);
-       % [~,~,~,~,~,~,~,~,NomTypeImaB]=fileparts_uvmat(Data.Civ1_ImageB);
-             elseif isfield(Data,'Civ2_ImageA')
-                 series('display_file_name',hhseries,Data.Civ2_ImageA,'append');%append the image series to the input list
-                         [~,~,~,~,~,~,~,~,NomTypeImaA]=fileparts_uvmat(Data.Civ2_ImageA);
-        %[RootPath,SubDir,RootFile,i1,i2,j1,j2,FileExt,NomTypeImaB]=fileparts_uvmat(Data.Civ2_ImageB);
+            
+            Param.InputTable(2,:)=Param.InputTable(1,:);
+              set(hhseries.InputTable,'Data',Param.InputTable)
+             if isfield(Data,'Civ2_ImageA')
+                 ImageName=Data.Civ2_ImageA;
+             elseif isfield(Data,'Civ1_ImageA')
+                 ImageName=Data.Civ1_ImageA;
              end
-        end
-
-        iview_image=2;%line # for the input images
+            series('display_file_name',hhseries,ImageName,1);%append the image series to the input list
+                    [~,~,~,~,~,~,~,~,NomTypeImaA]=fileparts_uvmat(ImageName);
+      
+%              elseif isfield(Data,'Civ2_ImageA')
+%                  series('display_file_name',hhseries,Data.Civ2_ImageA,'append');%append the image series to the input list
+%                          [~,~,~,~,~,~,~,~,NomTypeImaA]=fileparts_uvmat(Data.Civ2_ImageA);
+%      
+%              end
+        end
+
+        iview_image=1;%line # for the input images
     otherwise 
         % if ~strcmp(FileType,'image')
@@ -162 +162 @@
      return
 end
-MaxIndex_i=Param.IndexRange.MaxIndex_i(iview_image);
-MinIndex_i=Param.IndexRange.MinIndex_i(iview_image);
+MaxIndex_i=Param.IndexRange.MaxIndex_i(1);
+MinIndex_i=Param.IndexRange.MinIndex_i(1);
 MaxIndex_j=1;%default
 MinIndex_j=1;
@@ -262 +262 @@
         for index = 1:ind_opening
             set(handles.(ListOptions{index}),'value',0)
-        end
-        for index = ind_opening+1:6
-            set(handles.(ListOptions{index}),'value',1)
-        end
+            set(handles.(ListOptions{index}),'String',regexprep(ListOptions{index},'Check','redo '))
+        end
+%         for index = ind_opening+1:6
+%             set(handles.(ListOptions{index}),'value',1)
+%         end
         set(handles.CheckCiv3,'Visible','off')% make visible the switch 'iterate/repet' for Civ2.
         set(handles.CheckCiv3,'Value',0)% select'iterate/repet' by default
@@ -272 +273 @@
             set(handles.(ListOptions{index}),'value',0)
         end
-        for index = 4:6
-            set(handles.(ListOptions{index}),'value',1)
-        end
+%         for index = 4:6
+%             set(handles.(ListOptions{index}),'value',1)
+%         end
         set(handles.CheckCiv3,'Visible','on')% make visible the switch 'iterate/repet' for Civ2.
         set(handles.CheckCiv3,'Value',1)% select'iterate/repet' by default
@@ -290 +291 @@
 if ~checkrefresh && isfield(Param,'ActionInput')&& strcmp(Param.ActionInput.Program,Param.Action.ActionName)% the program fits with the stored data
     fill_GUI(Param.ActionInput,hObject);%fill the GUI with the parameters retrieved from the input Param
+
+    if isfield(Param.ActionInput,'Civ1')&& isfield(Param.ActionInput.Civ1,'SearchBoxSize')
+               SearchRange=round((Param.ActionInput.Civ1.SearchBoxSize-Param.ActionInput.Civ1.CorrBoxSize)/2);
+                set(handles.num_SearchRange_1(1),'String',num2str(SearchRange(1)))
+                set(handles.num_SearchRange_2(1),'String',num2str(SearchRange(2)))
+            end
+            if isfield(Param.ActionInput,'Civ2')&& isfield(Param.ActionInput.Civ2,'SearchBoxSize')
+               SearchRange=round((Param.ActionInput.Civ2.SearchBoxSize-Param.ActionInput.Civ2.CorrBoxSize)/2);
+                set(handles.num_SearchRange_1(2),'String',num2str(SearchRange(1)))
+                set(handles.num_SearchRange_2(2),'String',num2str(SearchRange(2)))
+            end
     hcheckgrid=findobj(handles.civ_input,'Tag','CheckGrid');
     for ilist=1:numel(hcheckgrid)
@@ -416 +428 @@
 %Param.CheckCiv1=1;
 Param.Civ1.CorrBoxSize=[31 31 1];
-Param.Civ1.SearchBoxSize=[61 61 5];
+Param.Civ1.SearchRange=[15 15 2];
 Param.Civ1.SearchBoxShift=[0 0];
 Param.Civ1.CorrSmooth=1;
@@ -439 +451 @@
 %Param.CheckCiv2=1;
 Param.Civ2.CorrBoxSize=[21 21];
-Param.Civ2.SearchBoxSize=[27 27];
+Param.Civ2.SearchRange=[3 3];
 Param.Civ2.CorrSmooth=1;
 Param.Civ2.Dx=10;
@@ -515 +527 @@
 %------------------------------------------------------------------------
 update_CivOptions(handles,0)
+update_frame(handles,'CheckCiv1')
+
 
 %------------------------------------------------------------------------
@@ -521 +535 @@
 %------------------------------------------------------------------------
 update_CivOptions(handles,0)
+update_frame(handles,'CheckFix1')
 
 %------------------------------------------------------------------------
@@ -527 +542 @@
 %------------------------------------------------------------------------
 update_CivOptions(handles,0)
+update_frame(handles,'CheckPatch1')
+
 
 %------------------------------------------------------------------------
@@ -533 +550 @@
 %------------------------------------------------------------------------
 update_CivOptions(handles,0)
+update_frame(handles,'CheckCiv2')
 
 %------------------------------------------------------------------------
@@ -539 +557 @@
 %------------------------------------------------------------------------
 update_CivOptions(handles,0)
+update_frame(handles,'CheckFix2')
 
 %------------------------------------------------------------------------
@@ -545 +564 @@
 %------------------------------------------------------------------------
 update_CivOptions(handles,0)
+update_frame(handles,'CheckPatch2')
+
+function update_frame(handles,option)
+if get(handles.(option),'Value')
+    option=regexprep(option,'Check','');
+set(handles.(option),'Visible','on')
+children=get(handles.(option),'children');
+set(children,'Enable','on')
+else
+    option=regexprep(option,'Check','');
+    set(handles.(option),'Visible','off')
+end
 
 %------------------------------------------------------------------------
@@ -591 +622 @@
 for ilist=1:length(options)
     if checkbox(ilist)
-        set(handles.(options{ilist}),'Visible','on')
+%          set(handles.(options{ilist}),'Visible','on')
+        set(handles.(options{ilist}),'Enable','on')
+%         set(handles.(['Check' options{ilist}]),'Strin
     else
-        set(handles.(options{ilist}),'Visible','off')
+%         set(handles.(options{ilist}),'Visible','off')
+        set(handles.(options{ilist}),'Enable','off')
     end
 end
@@ -608 +642 @@
     checkeven=(mod(ActionInput.Civ1.CorrBoxSize,2)==0);
     ActionInput.Civ1.CorrBoxSize(checkeven)=ActionInput.Civ1.CorrBoxSize(checkeven)+1;% set correlation box sizes to odd values
-    ActionInput.Civ1.SearchBoxSize(1:2)=max(ActionInput.Civ1.SearchBoxSize(1:2),ActionInput.Civ1.CorrBoxSize(1:2)+8);% insure that the search box size is large enough
-    checkeven=(mod(ActionInput.Civ1.SearchBoxSize,2)==0);
-    ActionInput.Civ1.SearchBoxSize(checkeven)=ActionInput.Civ1.SearchBoxSize(checkeven)+1;% set search box sizes to odd values
+    %ActionInput.Civ1.SearchBoxSize(1:2)=max(ActionInput.Civ1.SearchBoxSize(1:2),ActionInput.Civ1.CorrBoxSize(1:2)+8);% insure that the search box size is large enough
+    %checkeven=(mod(ActionInput.Civ1.SearchBoxSize,2)==0);
+    %ActionInput.Civ1.SearchBoxSize(checkeven)=ActionInput.Civ1.SearchBoxSize(checkeven)+1;% set search box sizes to odd values
 end
 if isfield(ActionInput,'Civ2')
     checkeven=(mod(ActionInput.Civ2.CorrBoxSize,2)==0);
     ActionInput.Civ2.CorrBoxSize(checkeven)=ActionInput.Civ2.CorrBoxSize(checkeven)+1;% set correlation box sizes to odd values
-    ActionInput.Civ2.SearchBoxSize=max(ActionInput.Civ2.SearchBoxSize,ActionInput.Civ2.CorrBoxSize+4);
-    checkeven=(mod(ActionInput.Civ2.SearchBoxSize,2)==0);
-    ActionInput.Civ2.SearchBoxSize(checkeven)=ActionInput.Civ2.SearchBoxSize(checkeven)+1;% set search box sizes to odd values
+    %ActionInput.Civ2.SearchBoxSize=max(ActionInput.Civ2.SearchBoxSize,ActionInput.Civ2.CorrBoxSize+4);
+    % checkeven=(mod(ActionInput.Civ2.SearchBoxSize,2)==0);
+    %ActionInput.Civ2.SearchBoxSize(checkeven)=ActionInput.Civ2.SearchBoxSize(checkeven)+1;% set search box sizes to odd values
 end
 
@@ -1002 +1036 @@
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %------------------------------------------------------------------------
-% --- Executes on button press in SearchRange: determine the search range num_SearchBoxSize_1,num_SearchBoxSize_2
+% --- Executes on button press in SearchRange: determine the search range num_SearchRange_1,num_SearchRange_2
 function SearchRange_Callback(hObject, eventdata, handles)
 %------------------------------------------------------------------------
@@ -1024 +1058 @@
 
 %------------------------------------------------------------------------
-% ---  determine the search range num_SearchBoxSize_1,num_SearchBoxSize_2 and shift
+% ---  determine the search range num_SearchRange_1,num_SearchRange_2 and shift
 function get_search_range(hObject, eventdata, handles)
 %------------------------------------------------------------------------
@@ -1056 +1090 @@
     set(handles.num_SearchBoxShift_1,'String',num2str(shiftx));
     set(handles.num_SearchBoxShift_2,'String',num2str(shifty));
-    set(handles.num_SearchBoxSize_1,'String',num2str(isx));
-    set(handles.num_SearchBoxSize_2,'String',num2str(isy));
+    set(handles.num_SearchRange_1,'String',num2str(isx));
+    set(handles.num_SearchRange_2,'String',num2str(isy));
 end
 
@@ -1087 +1121 @@
 currentdir=pwd;
 cd(Path);%move in the dir of the root name filebase
-gridfiles=dir([Name '_*grid_*.grid']);%look for grid files
+gridfiles=dir([Name '_*grid_*.nc']);%look for grid files
 cd(currentdir);%come back to the current working directory
 if ~isempty(gridfiles)
@@ -1168 +1202 @@
 function CheckGrid_Callback(hObject, eventdata, handles)
 %------------------------------------------------------------------------
-value=get(hObject,'Value');
-hparent=get(hObject,'parent');%handles of the parent panel
-hchildren=get(hparent,'children');
-handle_txtbox=findobj(hchildren,'tag','Grid');% look for the grid name box in the same panel
-handle_dx=findobj(hchildren,'tag','num_Dx');
-handle_dy=findobj(hchildren,'tag','num_Dy');
-handle_title_dx=findobj(hchildren,'tag','title_Dx');
-handle_title_dy=findobj(hchildren,'tag','title_Dy');
-testgrid=0;
-filegrid='';
-if value
-        hseries=findobj(allchild(0),'Tag','series');
+hparent=get(hObject,'parent');
+PanelName=get(hparent,'tag');
+handle_txtbox=handles.Grid
+if strcmp(PanelName,'civ2')
+    handle_txtbox=handle_txtbox(2);
+end
+% hchildren=get(hparent,'children');
+% handle_txtbox=findobj(hchildren,'tag','Grid');% look for the grid name box in the same panel
+% handle_NbSlice=findobj(hchildren,'tag','num_NbSlice');% look for the mask name box in the same panel
+testgrid=false;
+corrstatus='on';
+if get(hObject,'Value')% if the checkbox is activated
+    hseries=findobj(allchild(0),'Tag','series');
     hhseries=guidata(hseries);
     InputTable=get(hhseries.InputTable,'Data');
-     ind_A=1;% line index of the (first) image series
-    if strcmp(InputTable{1,5},'.nc');
-        ind_A=2;
-    end
-    filebase=InputTable{ind_A,1};
-    [nbslice, flag_grid]=get_grid(filebase,handles);% look for a grid with appropriate name 
-    if isequal(flag_grid,1)
-        filegrid=[num2str(nbslice) 'grid'];
-        testgrid=1;
-    else % browse for a grid 
-        filegrid=get(hObject,'UserData');%look for previous grid name stored as UserData
-        if exist(filegrid,'file')
-            filebase=filegrid;
-        end
-       filegrid = uigetfile_uvmat('pick a grid file .grid:',filebase,'.grid'); 
-        set(hObject,'UserData',filegrid);%store for future use
-        if ~isempty(filegrid)
-            testgrid=1;
-        end 
-        set(hObject,'UserData',filegrid);%store for future use
+    % browse for a grid file
+    filegrid= uigetfile_uvmat('pick a grid netcdf file (made by script_makegrid.m):',InputTable{1,1},'.nc');
+    if ~isempty(filegrid)
+        [FilePath,FileName,FileExt]=fileparts(filegrid);
+        Data=nc2struct(filegrid);
+        if isfield(Data,'Grid')
+        testgrid=true;
+        end
+        if isfield(Data,'CorrBox')
+        corrstatus='off';
+        end
     end
 end
 if testgrid
-    set(handle_dx,'Visible','off');
-    set(handle_dy,'Visible','off');
-    set(handle_title_dy,'Visible','off');
-    set(handle_title_dx,'Visible','off');
     set(handle_txtbox,'Visible','on')
     set(handle_txtbox,'String',filegrid)
-else
-    set(hObject,'Value',0);
-    set(handle_dx,'Visible','on'); 
-    set(handle_dy,'Visible','on');
-    set(handle_title_dy,'Visible','on');
-    set(handle_title_dx,'Visible','on');
+    set(handles.num_Dx,'Visible','off')
+    set(handles.num_Dy,'Visible','off')
+else
+    set(handles.num_Dx,'Visible','on')
+    set(handles.num_Dy,'Visible','on')
+    set(hObject,'Value',0)
     set(handle_txtbox,'Visible','off')
 end
 
+set(handles.num_CorrBoxSize_1,'Visible',corrstatus)
+set(handles.num_CorrBoxSize_2,'Visible',corrstatus)
+
+
 %% if hObject is on the checkciv1 frame, duplicate action for checkciv2 frame
-PanelName=get(hparent,'tag');
-if strcmp(PanelName,'Civ1')
-    hchildren=get(handles.Civ2,'children');
-    handle_checkbox=findobj(hchildren,'tag','CheckGrid');
-    handle_txtbox=findobj(hchildren,'tag','Grid');
-    handle_dx=findobj(hchildren,'tag','num_Dx');
-    handle_dy=findobj(hchildren,'tag','num_Dy');
-    handle_title_dx=findobj(hchildren,'tag','title_Dx');
-    handle_title_dy=findobj(hchildren,'tag','title_Dy');
-    set(handle_checkbox,'UserData',filegrid);%store for future use
-    if testgrid
-        set(handle_checkbox,'Value',1);
-        set(handle_dx,'Visible','off');
-        set(handle_dy,'Visible','off');
-        set(handle_title_dx,'Visible','off');
-        set(handle_title_dy,'Visible','off');
-        set(handle_txtbox,'Visible','on')
-        set(handle_txtbox,'String',filegrid)
-    end 
-end
+% PanelName=get(hparent,'tag');
+% if strcmp(PanelName,'Civ1')
+%     hchildren=get(handles.Civ2,'children');
+%     handle_checkbox=findobj(hchildren,'tag','CheckGrid');
+%     handle_txtbox=findobj(hchildren,'tag','Grid');
+%     handle_dx=findobj(hchildren,'tag','num_Dx');
+%     handle_dy=findobj(hchildren,'tag','num_Dy');
+%     handle_title_dx=findobj(hchildren,'tag','title_Dx');
+%     handle_title_dy=findobj(hchildren,'tag','title_Dy');
+%     %set(handle_checkbox,'UserData',filegrid);%store for future use
+%     if testgrid
+%         set(handle_checkbox,'Value',1);
+%         set(handle_dx,'Visible','off');
+%         set(handle_dy,'Visible','off');
+%         set(handle_title_dx,'Visible','off');
+%         set(handle_title_dy,'Visible','off');
+%         set(handle_txtbox,'Visible','on')
+%         set(handle_txtbox,'String',filegrid)
+%     end 
+% end
 set(handles.ConfigSource,'String','NEW')
 set(handles.OK,'BackgroundColor',[1 0 1])
@@ -1701 +1727 @@
 %% Civ param
 % lists of parameters to enter
-ListParamNum={'CorrBoxSize','SearchBoxSize','SearchBoxShift','Dx','Dy','Dz','MinIma','MaxIma'};% list of numerical values (to transform in strings)
+ListParamNum={'CorrBoxSize','SearchRange','SearchBoxShift','Dx','Dy','Dz','MinIma','MaxIma'};% list of numerical values (to transform in strings)
 ListParamValue={'CorrSmooth','CheckGrid','CheckMask','CheckThreshold'};
 ListParamString={'Grid','Mask'};
@@ -1738 +1764 @@
 function fill_panel(Data,handles,panel,ListParamNum,ListParamValue,ListParamString)
 children=get(handles.(panel),'children');%handles of the children of the input GUI with handle 'GUI_handle'
+set(children,'enable','off')
+set(handles.(panel),'Visible','on')
 handles_panel=[];
 for ichild=1:numel(children)
@@ -1812 +1840 @@
     check_input=0;
     if isfield(Param,'ActionInput')
-        if isfield(Param.ActionInput,'Program')&& strcmp(Param.ActionInput.Program,'civ_series')
+        if isfield(Param.ActionInput,'Program')&& ismember(Param.ActionInput.Program,{'civ_series','civ_3D'})
             fill_GUI(Param.ActionInput,handles.civ_input)% fill the elements of the GUI series with the input parameters
             set(handles.ConfigSource,'String',filexml)
             check_input=1;
+            if isfield(Param.ActionInput,'Civ1')&& isfield(Param.ActionInput.Civ1,'SearchBoxSize')
+               SearchRange=round((Param.ActionInput.Civ1.SearchBoxSize-Param.ActionInput.Civ1.CorrBoxSize)/2);
+                set(handles.num_SearchRange_1(1),'String',num2str(SearchRange(1)))
+                set(handles.num_SearchRange_2(1),'String',num2str(SearchRange(2)))
+            end
+            if isfield(Param.ActionInput,'Civ2')&& isfield(Param.ActionInput.Civ2,'SearchBoxSize')
+               SearchRange=round((Param.ActionInput.Civ2.SearchBoxSize-Param.ActionInput.Civ2.CorrBoxSize)/2);
+                set(handles.num_SearchRange_1(2),'String',num2str(SearchRange(1)))
+                set(handles.num_SearchRange_2(2),'String',num2str(SearchRange(2)))
+            end
             update_CivOptions(handles,0)             
         end
@@ -1875 +1913 @@
 
 
-function num_CorrBoxSize_3_Callback(hObject, eventdata, handles)
-
-
-function num_SearchBoxSize_3_Callback(hObject, eventdata, handles)
 
 
@@ -1886 +1920 @@
 % --- Executes on selection change in field_ref2.
 function field_ref2_Callback(hObject, eventdata, handles)
+
+
+
+function num_SearchRange_3_Callback(hObject, eventdata, handles)
+% hObject    handle to num_SearchRange_3 (see GCBO)
+% eventdata  reserved - to be defined in a future version of MATLAB
+% handles    structure with handles and user data (see GUIDATA)
+
+% Hints: get(hObject,'String') returns contents of num_SearchRange_3 as text
+%        str2double(get(hObject,'String')) returns contents of num_SearchRange_3 as a double

trunk/src/series/civ_series.m

@@ -75 +75 @@
     if isfield(Data,'ActionInput') && isfield(Data.ActionInput,'PairIndices') && strcmp(Data.ActionInput.PairIndices.ListPairMode,'pair j1-j2')
         Data.IndexRange_j='off';%no j index display in series
-        % if isfield(Data.ActionInput.PairIndices,'ListPairCiv2')
-        %     str_civ=Data.ActionInput.PairIndices.ListPairCiv2;
-        % else
-        %     str_civ=Data.ActionInput.PairIndices.ListPairCiv1;
-        % end
-        % r=regexp(str_civ,'^j= (?<num1>[a-z])-(?<num2>[a-z])','names');
-        % if isempty(r)
-        %     r=regexp(str_civ,'^j= (?<num1>[A-Z])-(?<num2>[A-Z])','names');
-        %     if isempty(r)
-        %         r=regexp(str_civ,'^j= (?<num1>\d+)-(?<num2>\d+)','names');
-        %     end
-        % end
-        % if ~isempty(r)
-        %     Data.j_index_1=stra2num(r.num1);
-        %     Data.j_index_2=stra2num(r.num2);
-        % end
     else
         Data.IndexRange_j='on';% j index display in series if relevant
@@ -361 +345 @@
                     ImageName_A=Param.ActionInput.RefFile;
                 else
-                ImageName_A=fullfile_uvmat(RootPath_A,SubDir_A,RootFile_A,FileExt_A,NomType_A,i1_series_Civ1(ifield),[],j1_series_Civ1(ifield));
+                    ImageName_A=fullfile_uvmat(RootPath_A,SubDir_A,RootFile_A,FileExt_A,NomType_A,i1_series_Civ1(ifield),[],j1_series_Civ1(ifield));
                 end
                 if strcmp(FileExt_A,'.nc')% case of input images in format netcdf
@@ -393 +377 @@
                 end
                 ImageName_B=fullfile_uvmat(RootPath_B,SubDir_B,RootFile_B,FileExt_B,NomType_B,i2_series_Civ1(ifield),[],j2_series_Civ1(ifield));
-                if strcmp(FileExt_B,'.nc') % case of input images in format netcdf
-                    FieldName_B=Param.InputFields.FieldName;
-                    [DataIn,~,~,errormsg]=nc2struct(ImageName_B,{FieldName_B});
-                    par_civ1.ImageB=DataIn.(FieldName_B);
-                else % usual image formats for image B
-                    if isempty(FileType_B)
+                    if isempty(FileType_B)% determine the image type for the first field
                         [FileInfo_B,VideoObject_B]=get_file_info(ImageName_B);
                         FileType_B=FileInfo_B.FileType;
@@ -407 +386 @@
                     end
                     [par_civ1.ImageB,VideoObject_B] = read_image(ImageName_B,FileType_B,VideoObject_B,FrameIndex_B_Civ1(ifield));
-                end
             catch ME % display errors in reading input images
                 if ~isempty(ME.message)
@@ -414 +392 @@
                 end
             end
-            par_civ1.ImageWidth=size(par_civ1.ImageA,2);
-            par_civ1.ImageHeight=size(par_civ1.ImageA,1);
+            % par_civ1.ImageWidth=size(par_civ1.ImageA,2);
+            % par_civ1.ImageHeight=size(par_civ1.ImageA,1);
             list_param=(fieldnames(Param.ActionInput.Civ1))';
             list_param(strcmp('TestCiv1',list_param))=[];% remove the parameter TestCiv1 from the list
@@ -460 +438 @@
         Data.VarAttribute{5}.Role='ancillary';
         Data.VarAttribute{6}.Role='errorflag';
+
         % case of mask
         if par_civ1.CheckMask&&~isempty(par_civ1.Mask)
-            if isfield(par_civ1,'NbSlice')
-                [RootPath_mask,SubDir_mask,RootFile_mask,~,~,~,~,Ext_mask]=fileparts_uvmat(Param.ActionInput.Civ1.Mask);
+            [RootPath_mask,SubDir_mask,RootFile_mask,~,~,~,~,Ext_mask]=fileparts_uvmat(Param.ActionInput.Civ1.Mask);
+            j1=1;
+            if ~isempty(j1_series_Civ1)
+                j1=j1_series_Civ1(ifield);
+            end
+            if ~isempty(i2_series_Civ1)% case of volume,masks act on different j levels
+                maskname=fullfile_uvmat(RootPath_mask,SubDir_mask,RootFile_mask,Ext_mask,'_1',j1);
+            elseif isfield(par_civ1,'NbSlice')
                 i1_mask=mod(i1-1,par_civ1.NbSlice)+1;
                 maskname=fullfile_uvmat(RootPath_mask,SubDir_mask,RootFile_mask,Ext_mask,'_1',i1_mask);
@@ -489 +474 @@
             end
         end
-        if strcmp(Param.ActionInput.ListCompareMode, 'PIV volume')
-            % Data.ListVarName=[Data.ListVarName 'Civ1_Z'];
-            % Data.Civ1_X=[];Data.Civ1_Y=[];Data.Civ1_Z=[];
-            % Data.Civ1_U=[];Data.Civ1_V=[];Data.Civ1_C=[];
-            % for ivol=1:NbSlice
-            %     % caluclate velocity data (y and v in indices, reverse to y component)
-            %     [xtable, ytable, utable, vtable, ctable, F, result_conv, errormsg] = civ (par_civ1);
-            %     if ~isempty(errormsg)
-            %         disp_uvmat('ERROR',errormsg,checkrun)
-            %         return
-            %     end
-            %     Data.Civ1_X=[Data.Civ1_X reshape(xtable,[],1)];
-            %     Data.Civ1_Y=[Data.Civ1_Y reshape(Param.Civ1.ImageHeight-ytable+1,[],1)];
-            %     Data.Civ1_Z=[Data.Civ1_Z ivol*ones(numel(xtable),1)];% z=image index in image coordinates
-            %     Data.Civ1_U=[Data.Civ1_U reshape(utable,[],1)];
-            %     Data.Civ1_V=[Data.Civ1_V reshape(-vtable,[],1)];
-            %     Data.Civ1_C=[Data.Civ1_C reshape(ctable,[],1)];
-            %     Data.Civ1_FF=[Data.Civ1_FF reshape(F,[],1)];
-            % end
-        else %usual PIV
-            % caluclate velocity data (y and v in indices, reverse to y component)
-            tstart_civ1=tic;
-            [xtable, ytable, utable, vtable, ctable, FF, result_conv, errormsg] = civ (par_civ1);
-            if ~isempty(errormsg)
-                disp_uvmat('ERROR',errormsg,checkrun)
-                return
-            end
-            Data.Civ1_X=reshape(xtable,[],1);
-            Data.Civ1_Y=reshape(par_civ1.ImageHeight-ytable+1,[],1);
-            Data.Civ1_U=reshape(utable,[],1);
-            Data.Civ1_V=reshape(-vtable,[],1);
-            Data.Civ1_C=reshape(ctable,[],1);
-            Data.Civ1_FF=reshape(FF,[],1);
-            time_civ1=toc(tstart_civ1);
-        end
+
+        % case of input grid
+        if par_civ1.CheckGrid &&~isempty(par_civ1.Grid)
+            GridData=nc2struct(Param.ActionInput.Civ1.Grid);
+            par_civ1.Grid=GridData.Grid;
+            par_civ1.CorrBoxSize=GridData.CorrBox;
+        end
+
+        % caluclate velocity data 
+        tstart_civ1=tic;
+        [Data.Civ1_X,Data.Civ1_Y,Data.Civ1_U,Data.Civ1_V,Data.Civ1_C,Data.Civ1_FF, result_conv, errormsg] = civ (par_civ1);
+        if ~isempty(errormsg)
+            disp_uvmat('ERROR',errormsg,checkrun)
+            return
+        end
+
     else% we use existing Civ1 data
         if exist('ncfile','var')
@@ -603 +569 @@
 
         % perform Patch calculation using the UVMAT fct 'filter_tps'
-        [Data.Civ1_SubRange,Data.Civ1_NbCentres,Data.Civ1_Coord_tps,Data.Civ1_U_tps,Data.Civ1_V_tps,tild,Ures, Vres,tild,FFres]=...
+        [Data.Civ1_SubRange,Data.Civ1_NbCentres,Data.Civ1_Coord_tps,Data.Civ1_U_tps,Data.Civ1_V_tps,~,Ures, Vres,~,FFres]=...
             filter_tps([Data.Civ1_X(ind_good) Data.Civ1_Y(ind_good)],Data.Civ1_U(ind_good),Data.Civ1_V(ind_good),[],Data.Patch1_SubDomainSize,Data.Patch1_FieldSmooth,Data.Patch1_MaxDiff);
         Data.Civ1_U_smooth(ind_good)=Ures;% take the interpolated (smoothed) velocity values for good vectors, keep civ1 data for the other
@@ -637 +603 @@
             end
             [FileInfo_A,VideoObject_A]=get_file_info(ImageName_A_Civ2);
-            par_civ2.ImageWidth=FileInfo_A.Width;
-            par_civ2.ImageHeight=FileInfo_A.Height;
-            if isfield(par_civ2,'Grid')% grid points set as input file
-                if ischar(par_civ2.Grid)%read the grid file if the input is a file name
-                    par_civ2.Grid=dlmread(par_civ2.Grid);
-                    par_civ2.Grid(1,:)=[];%the first line must be removed (heading in the grid file)
-                end
+            [npy_ima,npx_ima]=size(par_civ2.ImageA(:,:));
+            % par_civ2.ImageWidth=FileInfo_A.Width;
+            % par_civ2.ImageHeight=FileInfo_A.Height;
+
+            if par_civ2.CheckGrid &&~isempty(par_civ2.Grid) % case of input grid
+                GridData=nc2struct(Param.ActionInput.Civ2.Grid);
+                par_civ2.Grid=GridData.Grid;
+                par_civ2.CorrBoxSize=GridData.CorrBox;
+
             else% automatic grid
-                minix=floor(par_civ2.Dx/2)-0.5;
-                maxix=minix+par_civ2.Dx*floor((par_civ2.ImageWidth-1)/par_civ2.Dx);
-                miniy=floor(par_civ2.Dy/2)-0.5;
-                maxiy=minix+par_civ2.Dy*floor((par_civ2.ImageHeight-1)/par_civ2.Dy);
-                [GridX,GridY]=meshgrid(minix:par_civ2.Dx:maxix,miniy:par_civ2.Dy:maxiy);
+                nbinterv_x=floor((npx_ima-1)/par_civ2.Dx);
+                gridlength_x=nbinterv_x*par_civ2.Dx;
+                minix=ceil((npx_ima-gridlength_x)/2);
+                nbinterv_y=floor((npy_ima-1)/par_civ2.Dy);
+                gridlength_y=nbinterv_y*par_civ2.Dy;
+                miniy=ceil((npy_ima-gridlength_y)/2);
+                [GridX,GridY]=meshgrid(minix:par_civ2.Dx:npx_ima-1,miniy:par_civ2.Dy:npy_ima-1);
+                par_civ2.Grid=zeros(numel(GridX),2);
                 par_civ2.Grid(:,1)=reshape(GridX,[],1);
-                par_civ2.Grid(:,2)=reshape(GridY,[],1);
-            end
-        end
-        
+                par_civ2.Grid(:,2)=reshape(GridY,[],1);% increases with array index
+
+            end
+        end
+
         % get the guess from patch1 or patch2 (case 'CheckCiv3')
-        if CheckInputFile % read input images (except in mode Test where it is introduced directly in Param.ActionInput.Civ1.ImageNameA and B)
-            if isfield (par_civ2,'CheckCiv3') && par_civ2.CheckCiv3 %get the guess from  patch2
-                SubRange= Data.Civ2_SubRange;
-                NbCentres=Data.Civ2_NbCentres;
-                Coord_tps=Data.Civ2_Coord_tps;
-                U_tps=Data.Civ2_U_tps;
-                V_tps=Data.Civ2_V_tps;
-                CivStage=Data.CivStage;%store the current CivStage
-                Civ1_Dt=Data.Civ2_Dt;
-                Data=[];%reinitialise the result structure Data
-                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.CivStage=CivStage+1;%update the current civStage after reinitialisation of Data
-                Data.ListVarName={};
-                Data.VarDimName={};
-            else % get the guess from patch1
-                SubRange= Data.Civ1_SubRange;
-                NbCentres=Data.Civ1_NbCentres;
-                Coord_tps=Data.Civ1_Coord_tps;
-                U_tps=Data.Civ1_U_tps;
-                V_tps=Data.Civ1_V_tps;
-                Civ1_Dt=Data.Civ1_Dt;
-                Data.CivStage=4;
-            end
-        else
-            SubRange= par_civ2.Civ1_SubRange;
-            NbCentres=par_civ2.Civ1_NbCentres;
-            Coord_tps=par_civ2.Civ1_Coord_tps;
-            U_tps=par_civ2.Civ1_U_tps;
-            V_tps=par_civ2.Civ1_V_tps;
-            Civ1_Dt=par_civ2.Civ1_Dt;
-            Civ2_Dt=par_civ2.Civ1_Dt;
+        % if CheckInputFile % read input images (except in mode Test where it is introduced directly in Param.ActionInput.Civ1.ImageNameA and B)
+        if isfield (par_civ2,'CheckCiv3') && par_civ2.CheckCiv3 %get the guess from  patch2
+            SubRange= Data.Civ2_SubRange;
+            NbCentres=Data.Civ2_NbCentres;
+            Coord_tps=Data.Civ2_Coord_tps;
+            U_tps=Data.Civ2_U_tps;
+            V_tps=Data.Civ2_V_tps;
+            CivStage=Data.CivStage;%store the current CivStage
+            Civ1_Dt=Data.Civ2_Dt;
+            Data=[];%reinitialise the result structure Data
+            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.CivStage=CivStage+1;%update the current civStage after reinitialisation of Data
             Data.ListVarName={};
             Data.VarDimName={};
-        end
+        else % get the guess from patch1
+            SubRange= Data.Civ1_SubRange;
+            NbCentres=Data.Civ1_NbCentres;
+            Coord_tps=Data.Civ1_Coord_tps;
+            U_tps=Data.Civ1_U_tps;
+            V_tps=Data.Civ1_V_tps;
+            Civ1_Dt=Data.Civ1_Dt;
+            Data.CivStage=4;
+        end
+        % else
+        %     SubRange= par_civ2.Civ1_SubRange;
+        %     NbCentres=par_civ2.Civ1_NbCentres;
+        %     Coord_tps=par_civ2.Civ1_Coord_tps;
+        %     U_tps=par_civ2.Civ1_U_tps;
+        %     V_tps=par_civ2.Civ1_V_tps;
+        %     Civ1_Dt=par_civ2.Civ1_Dt;
+        %     Civ2_Dt=par_civ2.Civ1_Dt;
+        %     Data.ListVarName={};
+        %     Data.VarDimName={};
+        % end
         Shiftx=zeros(size(par_civ2.Grid,1),1);% shift expected from civ1 data
         Shifty=zeros(size(par_civ2.Grid,1),1);
@@ -709 +681 @@
                 epoints = par_civ2.Grid(ind_sel,:);% coordinates of interpolation sites (measurement grids)
                 ctrs=Coord_tps(1:nbvec_sub,:,isub) ;%(=initial points) ctrs
-                nbval(ind_sel)=nbval(ind_sel)+1;% records the number of values for each interpolation point (in case of subdomain overlap)
                 EM = tps_eval(epoints,ctrs);% thin plate spline (tps) coefficient
-                Shiftx(ind_sel)=Shiftx(ind_sel)+EM*U_tps(1:nbvec_sub+3,isub);%velocity shift estimated by tps from civ1
-                Shifty(ind_sel)=Shifty(ind_sel)+EM*V_tps(1:nbvec_sub+3,isub);
+                CentreX=(SubRange(1,1,isub)+SubRange(1,2,isub))/2; %x posiion of the subdomain center
+                CentreY=(SubRange(2,1,isub)+SubRange(2,2,isub))/2; %y posiion of the subdomain center
+                xwidth=(SubRange(1,2,isub)-SubRange(1,1,isub))/pi;
+                ywidth=(SubRange(2,2,isub)-SubRange(2,1,isub))/pi;
+                x_dist=(epoints(:,1)-CentreX)/xwidth;
+                y_dist=(epoints(:,2)-CentreY)/ywidth;
+                weight=cos(x_dist).*cos(y_dist);%weighting fct =1 at the rectangle center and 0 at edge
+                nbval(ind_sel)=nbval(ind_sel)+weight;% records the number of values for each interpolation point (in case of subdomain overlap)
+                Shiftx(ind_sel)=Shiftx(ind_sel)+weight.*(EM*U_tps(1:nbvec_sub+3,isub));%velocity shift estimated by tps from civ1
+                Shifty(ind_sel)=Shifty(ind_sel)+weight.*(EM*V_tps(1:nbvec_sub+3,isub));
                 if par_civ2.CheckDeformation
                     [EMDX,EMDY] = tps_eval_dxy(epoints,ctrs);%2D matrix of distances between extrapolation points epoints and spline centres (=site points) ctrs
-                    DUDX(ind_sel)=DUDX(ind_sel)+EMDX*U_tps(1:nbvec_sub+3,isub);
-                    DUDY(ind_sel)=DUDY(ind_sel)+EMDY*U_tps(1:nbvec_sub+3,isub);
-                    DVDX(ind_sel)=DVDX(ind_sel)+EMDX*V_tps(1:nbvec_sub+3,isub);
-                    DVDY(ind_sel)=DVDY(ind_sel)+EMDY*V_tps(1:nbvec_sub+3,isub);
-                end
-            end
-        end
-        i1=i1_series_Civ2(ifield);
-        i2=i1;
-        if ~isempty(i2_series_Civ2)
-            i2=i2_series_Civ1(ifield);
-        end
-        j1=1;
-        if ~isempty(j1_series_Civ2)
-            j1=j1_series_Civ1(ifield);
-        end
-        j2=j1;
-        if ~isempty(j2_series_Civ2)
-            j2=j2_series_Civ2(ifield);
-        end
-        if par_civ2.CheckMask&&~isempty(par_civ2.Mask)
-            if isfield(par_civ2,'NbSlice')
-                [RootPath_mask,SubDir_mask,RootFile_mask,~,~,~,~,Ext_mask]=fileparts_uvmat(Param.ActionInput.Civ2.Mask);
+                    DUDX(ind_sel)=DUDX(ind_sel)+weight.*(EMDX*U_tps(1:nbvec_sub+3,isub));
+                    DUDY(ind_sel)=DUDY(ind_sel)+weight.*(EMDY*U_tps(1:nbvec_sub+3,isub));
+                    DVDX(ind_sel)=DVDX(ind_sel)+weight.*(EMDX*V_tps(1:nbvec_sub+3,isub));
+                    DVDY(ind_sel)=DVDY(ind_sel)+weight.*(EMDY*V_tps(1:nbvec_sub+3,isub));
+                end
+            end
+        end
+        Shiftx(nbval>0)=Shiftx(nbval>0)./nbval(nbval>0);
+        Shifty(nbval>0)=Shifty(nbval>0)./nbval(nbval>0);
+
+        % introduce mask
+        if par_civ2.CheckMask && ~isempty(par_civ2.Mask)
+            [RootPath_mask,SubDir_mask,RootFile_mask,~,~,~,~,Ext_mask]=fileparts_uvmat(Param.ActionInput.Civ2.Mask);
+            if ~isempty(i2_series_Civ2) % we do PIV among indices i,  at given indices j (volume scan), mask depends on position j
+                j1=1;
+                if ~isempty(j1_series_Civ2)
+                    j1=j1_series_Civ1(ifield);
+                end
+                maskname=fullfile_uvmat(RootPath_mask,SubDir_mask,RootFile_mask,Ext_mask,'_1',j1);
+            elseif isfield(par_civ2,'NbSlice')
+                i1=i1_series_Civ2(ifield);
                 i1_mask=mod(i1-1,par_civ2.NbSlice)+1;
                 maskname=fullfile_uvmat(RootPath_mask,SubDir_mask,RootFile_mask,Ext_mask,'_1',i1_mask);
@@ -747 +724 @@
             else
                 if exist(maskname,'file')
-                try
-                    par_civ2.Mask=imread(maskname);%update the mask, an store it for future use
-                catch ME
-                    if ~isempty(ME.message)
-                        errormsg=['error reading input image: ' ME.message];
-                        disp_uvmat('ERROR',errormsg,checkrun)
-                        return
+                    try
+                        par_civ2.Mask=imread(maskname);%update the mask, an store it for future use
+                    catch ME
+                        if ~isempty(ME.message)
+                            errormsg=['error reading input image: ' ME.message];
+                            disp_uvmat('ERROR',errormsg,checkrun)
+                            return
+                        end
                     end
-                end
                 else
                     par_civ2.Mask=[];
@@ -772 +749 @@
             end
         end
-        par_civ2.SearchBoxShift=(Civ2_Dt/Civ1_Dt)*[Shiftx(nbval>=1)./nbval(nbval>=1) Shifty(nbval>=1)./nbval(nbval>=1)];
+        par_civ2.SearchBoxShift=zeros(size(par_civ2.Grid));
+        par_civ2.SearchBoxShift(:,1)=(Civ2_Dt/Civ1_Dt)*Shiftx;
+        par_civ2.SearchBoxShift(:,2)=(Civ2_Dt/Civ1_Dt)*Shifty;
         % shift the grid points by half the expected shift to provide the correlation box position in image A
-        par_civ2.Grid=[par_civ2.Grid(nbval>=1,1)-par_civ2.SearchBoxShift(:,1)/2 par_civ2.Grid(nbval>=1,2)-par_civ2.SearchBoxShift(:,2)/2];
+        %par_civ2.Grid=[par_civ2.Grid(nbval>=1,1)-par_civ2.SearchBoxShift(:,1)/2 par_civ2.Grid(nbval>=1,2)-par_civ2.SearchBoxShift(:,2)/2];
+
         if par_civ2.CheckDeformation
-            par_civ2.DUDX=DUDX(nbval>=1)./nbval(nbval>=1);
-            par_civ2.DUDY=DUDY(nbval>=1)./nbval(nbval>=1);
-            par_civ2.DVDX=DVDX(nbval>=1)./nbval(nbval>=1);
-            par_civ2.DVDY=DVDY(nbval>=1)./nbval(nbval>=1);
-        end
-        
+            par_civ2.DUDX(nbval>0)=DUDX(nbval>0)./nbval(nbval>0);
+            par_civ2.DUDY(nbval>0)=DUDY(nbval>0)./nbval(nbval>0);
+            par_civ2.DVDX(nbval>0)=DVDX(nbval>0)./nbval(nbval>0);
+            par_civ2.DVDY(nbval>0)=DVDY(nbval>0)./nbval(nbval>0);
+        end
+
         % calculate velocity data (y and v in image indices, reverse to y component)
-        
-        [xtable, ytable, utable, vtable, ctable, FF,result_conv,errormsg] = civ (par_civ2);
-        
+
+        [Data.Civ2_X,Data.Civ2_Y,Data.Civ2_U,Data.Civ2_V,Data.Civ2_C,Data.Civ2_FF,~, errormsg] = civ (par_civ2);
+
         list_param=(fieldnames(Param.ActionInput.Civ2))';
         list_param(strcmp('TestCiv2',list_param))=[];% remove the parameter TestCiv2 from the list
@@ -806 +786 @@
         end
         Data.ListGlobalAttribute=[Data.ListGlobalAttribute Civ2_param];
-        
+
         nbvar=numel(Data.ListVarName);
         % define the Civ2 variable (if Civ2 data are not replaced from previous calculation)
@@ -819 +799 @@
             Data.VarAttribute{nbvar+6}.Role='errorflag';
         end
-        Data.Civ2_X=reshape(xtable,[],1);
-        Data.Civ2_Y=reshape(size(par_civ2.ImageA,1)-ytable+1,[],1);
-        Data.Civ2_U=reshape(utable,[],1);
-        Data.Civ2_V=reshape(-vtable,[],1);
-        Data.Civ2_C=reshape(ctable,[],1);
-        Data.Civ2_FF=reshape(FF,[],1);
         disp('civ2 performed')
         time_civ2=toc(tstart_civ2);
@@ -837 +811 @@
         end
     end
-    
+
     %% Fix2
     if isfield (Param.ActionInput,'Fix2')
@@ -847 +821 @@
             Data.(Fix2_param{ilist})=Param.ActionInput.Fix2.(list_param{ilist});
         end
-        Data.ListGlobalAttribute=[Data.ListGlobalAttribute Fix2_param];     
+        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.CivStage=Data.CivStage+1;
@@ -916 +890 @@
             disp(['time other ' num2str((time_total-time_input-time_civ1-time_patch1-time_civ2-time_patch2),2) ' s'])
         end
-    end
-end
-
-
-% 'civ': function piv.m adapted from PIVlab http://pivlab.blogspot.com/
-%--------------------------------------------------------------------------
-% function [xtable ytable utable vtable typevector] = civ (image1,image2,ibx,iby step, subpixfinder, mask, roi)
-%
-% OUTPUT:
-% xtable: set of x coordinates
-% ytable: set of y coordinates
-% utable: set of u displacements (along x)
-% vtable: set of v displacements (along y)
-% ctable: max image correlation for each vector
-% typevector: set of flags, =1 for good, =0 for NaN vectors
-%
-%INPUT:
-% par_civ: structure of input parameters, with fields:
-%  .ImageA: first image for correlation (matrix)
-%  .ImageB: second image for correlation(matrix)
-%  .CorrBoxSize: 1,2 vector giving the size of the correlation box in x and y
-%  .SearchBoxSize:  1,2 vector giving the size of the search box in x and y
-%  .SearchBoxShift: 1,2 vector or 2 column matrix (for civ2) giving the shift of the search box in x and y
-%  .CorrSmooth: =1 or 2 determines the choice of the sub-pixel determination of the correlation max
-%  .ImageWidth: nb of pixels of the image in x
-%  .Dx, Dy: mesh for the PIV calculation
-%  .Grid: grid giving the PIV calculation points (alternative to .Dx .Dy): centres of the correlation boxes in Image A
-%  .Mask: name of a mask file or mask image matrix itself
-%  .MinIma: thresholds for image luminosity
-%  .MaxIma
-%  .CheckDeformation=1 for subpixel interpolation and image deformation (linear transform)
-%  .DUDX: matrix of deformation obtained from patch at each grid point
-%  .DUDY
-%  .DVDX:
-%  .DVDY
-
-function [xtable,ytable,utable,vtable,ctable,FF,result_conv,errormsg] = civ (par_civ)
-
-%% prepare measurement grid
-if isfield(par_civ,'Grid')% grid points set as input, central positions of the sub-images in image A
-    if ischar(par_civ.Grid)%read the grid file if the input is a file name (grid in x, y image coordinates)
-        par_civ.Grid=dlmread(par_civ.Grid);
-        par_civ.Grid(1,:)=[];%the first line must be removed (heading in the grid file)
-    end
-    % else par_civ.Grid is already an array, no action here
-else% automatic grid in x, y image coordinates
-nbinterv_x=floor((par_civ.ImageWidth-1)/par_civ.Dx);
-gridlength_x=nbinterv_x*par_civ.Dx;
-minix=ceil((par_civ.ImageWidth-gridlength_x)/2);
-nbinterv_y=floor((par_civ.ImageHeight-1)/par_civ.Dy);
-gridlength_y=nbinterv_y*par_civ.Dy;
-miniy=ceil((par_civ.ImageHeight-gridlength_y)/2);
-[GridX,GridY]=meshgrid(minix:par_civ.Dx:par_civ.ImageWidth-1,miniy:par_civ.Dy:par_civ.ImageHeight-1);
-    par_civ.Grid(:,1)=reshape(GridX,[],1);
-    par_civ.Grid(:,2)=reshape(GridY,[],1);% increases with array index
-    % 
-    % minix=floor(par_civ.Dx/2)-0.5;
-    % maxix=minix+par_civ.Dx*floor((par_civ.ImageWidth-1)/par_civ.Dx);
-    % miniy=floor(par_civ.Dy/2)-0.5;% first automatic grid point at half the mesh Dy
-    % maxiy=minix+par_civ.Dy*floor((par_civ.ImageHeight-1)/par_civ.Dy);
-    % [GridX,GridY]=meshgrid(minix:par_civ.Dx:maxix,miniy:par_civ.Dy:maxiy);
-    % par_civ.Grid(:,1)=reshape(GridX,[],1);
-    % par_civ.Grid(:,2)=reshape(GridY,[],1);% increases with array index
-end
-nbvec=size(par_civ.Grid,1);
-
-%% prepare correlation and search boxes
-ibx2=floor(par_civ.CorrBoxSize(1)/2);
-iby2=floor(par_civ.CorrBoxSize(2)/2);
-isx2=floor(par_civ.SearchBoxSize(1)/2);
-isy2=floor(par_civ.SearchBoxSize(2)/2);
-shiftx=round(par_civ.SearchBoxShift(:,1));%use the input shift estimate, rounded to the next integer value
-shifty=-round(par_civ.SearchBoxShift(:,2));% sign minus because image j index increases when y decreases
-if numel(shiftx)==1% case of a unique shift for the whole field( civ1)
-    shiftx=shiftx*ones(nbvec,1);
-    shifty=shifty*ones(nbvec,1);
-end
-%TODO: shift the origin by -shift/2
-
-%% Array initialisation and default output  if par_civ.CorrSmooth=0 (just the grid calculated, no civ computation)
-xtable=round(par_civ.Grid(:,1)+0.5)-0.5;
-ytable=round(par_civ.ImageHeight-par_civ.Grid(:,2)+0.5)-0.5;% y index corresponding to the position in image coordiantes
-utable=shiftx;%zeros(nbvec,1);
-vtable=shifty;%zeros(nbvec,1);
-ctable=zeros(nbvec,1);
-FF=zeros(nbvec,1);
-result_conv=[];
-errormsg='';
-
-%% prepare mask
-if isfield(par_civ,'Mask') && ~isempty(par_civ.Mask)
-    if strcmp(par_civ.Mask,'all')
-        return    % get the grid only, no civ calculation
-    elseif ischar(par_civ.Mask)
-        par_civ.Mask=imread(par_civ.Mask);% read the mask if not allready done
-    end
-end
-check_MinIma=isfield(par_civ,'MinIma');% test for image luminosity threshold
-check_MaxIma=isfield(par_civ,'MaxIma') && ~isempty(par_civ.MaxIma);
-
-par_civ.ImageA=sum(double(par_civ.ImageA),3);%sum over rgb component for color images
-par_civ.ImageB=sum(double(par_civ.ImageB),3);
-[npy_ima,npx_ima]=size(par_civ.ImageA);
-if ~isequal(size(par_civ.ImageB),[npy_ima npx_ima])
-    errormsg='image pair with unequal size';
-    return
-end
-
-%% Apply mask
-% Convention for mask, IDEAS NOT IMPLEMENTED 
-% mask >200 : velocity calculated
-%  200 >=mask>150;velocity not calculated, interpolation allowed (bad spots)
-% 150>=mask >100: velocity not calculated, nor interpolated
-%  100>=mask> 20: velocity not calculated, impermeable (no flux through mask boundaries)
-%  20>=mask: velocity=0
-checkmask=0;
-MinA=min(min(par_civ.ImageA));
-%MinB=min(min(par_civ.ImageB));
-%check_undefined=false(size(par_civ.ImageA));
-if isfield(par_civ,'Mask') && ~isempty(par_civ.Mask)
-    checkmask=1;
-    if ~isequal(size(par_civ.Mask),[npy_ima npx_ima])
-        errormsg='mask must be an image with the same size as the images';
-        return
-    end
-    check_undefined=(par_civ.Mask<200 & par_civ.Mask>=20 );
-end
-
-%% compute image correlations: MAINLOOP on velocity vectors
-corrmax=0;
-sum_square=1;% default
-mesh=1;% default
-CheckDeformation=isfield(par_civ,'CheckDeformation')&& par_civ.CheckDeformation==1;
-if CheckDeformation
-    mesh=0.25;%mesh in pixels for subpixel image interpolation (x 4 in each direction)
-    par_civ.CorrSmooth=2;% use SUBPIX2DGAUSS (take into account more points near the max)
-end
- 
-if par_civ.CorrSmooth~=0 % par_civ.CorrSmooth=0 implies no civ computation (just input image and grid points given)
-    for ivec=1:nbvec
-        iref=round(par_civ.Grid(ivec,1)+0.5);% xindex on the image A for the middle of the correlation box
-        jref=round(par_civ.ImageHeight-par_civ.Grid(ivec,2)+0.5);%  j index  for the middle of the correlation box in the image A
-        FF(ivec)=0;
-        subrange1_x=iref-ibx2:iref+ibx2;% x indices defining the first subimage
-        subrange1_y=jref-iby2:jref+iby2;% y indices defining the first subimage
-        subrange2_x=iref+shiftx(ivec)-isx2:iref+shiftx(ivec)+isx2;%x indices defining the second subimage
-        subrange2_y=jref+shifty(ivec)-isy2:jref+shifty(ivec)+isy2;%y indices defining the second subimage
-        image1_crop=MinA*ones(numel(subrange1_y),numel(subrange1_x));% default value=min of image A
-        image2_crop=MinA*ones(numel(subrange2_y),numel(subrange2_x));% default value=min of image A
-        check1_x=subrange1_x>=1 & subrange1_x<=par_civ.ImageWidth;% check which points in the subimage 1 are contained in the initial image 1
-        check1_y=subrange1_y>=1 & subrange1_y<=par_civ.ImageHeight;
-        check2_x=subrange2_x>=1 & subrange2_x<=par_civ.ImageWidth;% check which points in the subimage 2 are contained in the initial image 2
-        check2_y=subrange2_y>=1 & subrange2_y<=par_civ.ImageHeight;
-        image1_crop(check1_y,check1_x)=par_civ.ImageA(subrange1_y(check1_y),subrange1_x(check1_x));%extract a subimage (correlation box) from image A
-        image2_crop(check2_y,check2_x)=par_civ.ImageB(subrange2_y(check2_y),subrange2_x(check2_x));%extract a larger subimage (search box) from image B
-        if checkmask
-            mask1_crop=ones(numel(subrange1_y),numel(subrange1_x));% default value=1 for mask
-            mask2_crop=ones(numel(subrange2_y),numel(subrange2_x));% default value=1 for mask
-            mask1_crop(check1_y,check1_x)=check_undefined(subrange1_y(check1_y),subrange1_x(check1_x));%extract a mask subimage (correlation box) from image A
-            mask2_crop(check2_y,check2_x)=check_undefined(subrange2_y(check2_y),subrange2_x(check2_x));%extract a mask subimage (search box) from image B
-            sizemask=sum(sum(mask1_crop))/(numel(subrange1_y)*numel(subrange1_x));%size of the masked part relative to the correlation sub-image
-            if sizemask > 1/2% eliminate point if more than half of the correlation box is masked
-                FF(ivec)=1; %
-                utable(ivec)=NaN;
-                vtable(ivec)=NaN;
-            else
-                image1_crop=image1_crop.*~mask1_crop;% put to zero the masked pixels (mask1_crop='true'=1)
-                image2_crop=image2_crop.*~mask2_crop;
-                image1_mean=mean(mean(image1_crop))/(1-sizemask);
-                image2_mean=mean(mean(image2_crop))/(1-sizemask);
-            end
-        else
-            image1_mean=mean(mean(image1_crop));
-            image2_mean=mean(mean(image2_crop));
-        end
-        %threshold on image minimum
-        if FF(ivec)==0
-            if check_MinIma && (image1_mean < par_civ.MinIma || image2_mean < par_civ.MinIma)
-                FF(ivec)=1;
-                %threshold on image maximum
-            elseif check_MaxIma && (image1_mean > par_civ.MaxIma || image2_mean > par_civ.MaxIma)
-                FF(ivec)=1;
-            end
-            if FF(ivec)==1
-                utable(ivec)=NaN;
-                vtable(ivec)=NaN;
-            else
-                %mask
-                if checkmask
-                    image1_crop=(image1_crop-image1_mean).*~mask1_crop;%substract the mean, put to zero the masked parts
-                    image2_crop=(image2_crop-image2_mean).*~mask2_crop;
-                else
-                    image1_crop=(image1_crop-image1_mean);
-                    image2_crop=(image2_crop-image2_mean);
-                end
-                %deformation
-                if CheckDeformation
-                    xi=(1:mesh:size(image1_crop,2));
-                    yi=(1:mesh:size(image1_crop,1))';
-                    [XI,YI]=meshgrid(xi-ceil(size(image1_crop,2)/2),yi-ceil(size(image1_crop,1)/2));
-                    XIant=XI-par_civ.DUDX(ivec)*XI+par_civ.DUDY(ivec)*YI+ceil(size(image1_crop,2)/2);
-                    YIant=YI+par_civ.DVDX(ivec)*XI-par_civ.DVDY(ivec)*YI+ceil(size(image1_crop,1)/2);
-                    image1_crop=interp2(image1_crop,XIant,YIant);
-                    image1_crop(isnan(image1_crop))=0;
-                    xi=(1:mesh:size(image2_crop,2));
-                    yi=(1:mesh:size(image2_crop,1))';
-                    image2_crop=interp2(image2_crop,xi,yi,'*spline');
-                    image2_crop(isnan(image2_crop))=0;
-                end
-                sum_square=sum(sum(image1_crop.*image1_crop));
-                %reference: Oliver Pust, PIV: Direct Cross-Correlation
-                %%%%%% correlation calculation
-                result_conv= conv2(image2_crop,flip(flip(image1_crop,2),1),'valid');
-                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-                corrmax= max(max(result_conv));
-        
-                %result_conv=(result_conv/corrmax); %normalize, peak=always 255
-                %Find the correlation max, at 255
-                [y,x] = find(result_conv==corrmax,1);
-                subimage2_crop=image2_crop(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_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
-                if ~isempty(y) && ~isempty(x)
-                    try
-                        if par_civ.CorrSmooth==1
-                            [vector,FF(ivec)] = SUBPIXGAUSS (result_conv,x,y);
-                        elseif par_civ.CorrSmooth==2
-                            [vector,FF(ivec)] = SUBPIX2DGAUSS (result_conv,x,y);
-                        else
-                            [vector,FF(ivec)] = quadr_fit(result_conv,x,y);
-                        end
-                        utable(ivec)=vector(1)*mesh+shiftx(ivec);
-                        vtable(ivec)=vector(2)*mesh+shifty(ivec);
-                        xtable(ivec)=iref+utable(ivec)/2-0.5;% convec flow (velocity taken at the point middle from imgae 1 and 2)
-                        ytable(ivec)=jref+vtable(ivec)/2-0.5;% and position of pixel 1=0.5 (convention for image coordinates=0 at the edge)
-                        iref=round(xtable(ivec)+0.5);% nearest image index for the middle of the vector
-                        jref=round(ytable(ivec)+0.5);
-                        % 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))
-                            utable(ivec)=0;
-                            vtable(ivec)=0;
-                            FF(ivec)=1;
-                        end
-                        ctable(ivec)=corrmax/sum_square;% correlation value
-                    catch ME
-                        FF(ivec)=1;
-                        disp(ME.message)
-                    end
-                else
-                    FF(ivec)=1;
-                end
-            end
-        end
-    end
-end
-result_conv=result_conv*corrmax/sum_square;% keep the last (not normalised) correlation matrix for output
-
-%------------------------------------------------------------------------
-% --- Find the maximum of the correlation function with subpixel resolution
-% make a fit with a gaussian curve from the three correlation values across the max, along each direction.
-% OUPUT:
-% vector = optimum displacement vector with subpixel correction
-% FF =flag: =0 OK
-%           =1 , max too close to the edge of the search box (1 pixel margin)
-% INPUT:
-% result_conv: 2D correlation fct
-% x,y: position of the maximum correlation at integer values
-
-function [vector,FF] = SUBPIXGAUSS (result_conv,x,y)
-%------------------------------------------------------------------------
-% vector=[0 0]; %default
-FF=true;% error flag for vector truncated by the limited search box
-[npy,npx]=size(result_conv);
-
-peaky = y; peakx=x;
-if y < npy && y > 1 && x < npx-1 && x > 1
-   FF=false; % no error by the limited search box
-    max_conv=result_conv(y,x);% max correlation
-    %peak2noise= max(4,max_conv/std(reshape(result_conv,1,[])));% ratio of max conv to standard deviation of correlations (estiamtion of noise level), set to value 4 if it is too low
-    peak2noise=100;% TODO: make this threshold more precise, depending on the image noise
-    result_conv=result_conv*peak2noise/max_conv;% renormalise the correlation with respect to the noise 
-    result_conv(result_conv<1)=1; %set to 1 correlation values smaller than 1  (=0 by discretisation, to avoid divergence in the log)
-    
-    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);
-    f1 = log(result_conv(y,x-1));
-    f2 = log(result_conv(y,x+1));
-    peakx = peakx+ (f1-f2)/(2*f1-4*f0+2*f2);
-end
-
-vector=[peakx-floor(npx/2)-1 peaky-floor(npy/2)-1];
-
-%------------------------------------------------------------------------
-% --- Find the maximum of the correlation function after interpolation
-function [vector,FF] = SUBPIX2DGAUSS (result_conv,x,y)
-%------------------------------------------------------------------------
-% vector=[0 0]; %default
-FF=true;
-peaky=y;
-peakx=x;
-[npy,npx]=size(result_conv);
-if (x < npx) && (y < npy) && (x > 1) && (y > 1)
-    FF=false;
-max_conv=result_conv(y,x);% max correlation
-    peak2noise= max(4,max_conv/std(reshape(result_conv,1,[])));% ratio of max conv to standard deviation of correlations (estiamtion of noise level), set to value 4 if it is too low
-    result_conv=result_conv*peak2noise/max_conv;% renormalise the correlation with respect to the noise 
-    result_conv(result_conv<1)=1; %set to 1 correlation values smaller than 1  (=0 by discretisation, to avoid divergence in the log)
-    for i=-1:1
-        for j=-1:1
-  %following 15 lines based on  H. Nobach and M. Honkanen (2005)
-  % Two-dimensional Gaussian regression for sub-pixel displacement
-  % estimation in particle image velocimetry or particle position
-  % estimation in particle tracking velocimetry
-  % Experiments in Fluids (2005) 38: 511-515
-            c10(j+2,i+2)=i*log(result_conv(y+j, x+i));
-            c01(j+2,i+2)=j*log(result_conv(y+j, x+i));
-            c11(j+2,i+2)=i*j*log(result_conv(y+j, x+i));
-            c20(j+2,i+2)=(3*i^2-2)*log(result_conv(y+j, x+i));
-            c02(j+2,i+2)=(3*j^2-2)*log(result_conv(y+j, x+i));
-        end
-    end
-    c10=(1/6)*sum(sum(c10));
-    c01=(1/6)*sum(sum(c01));
-    c11=(1/4)*sum(sum(c11));
-    c20=(1/6)*sum(sum(c20));
-    c02=(1/6)*sum(sum(c02));
-    deltax=(c11*c01-2*c10*c02)/(4*c20*c02-c11*c11);
-    deltay=(c11*c10-2*c01*c20)/(4*c20*c02-c11*c11);
-    if abs(deltax)<1
-        peakx=x+deltax;
-    end
-    if abs(deltay)<1
-        peaky=y+deltay;
-    end
-end
-vector=[peakx-floor(npx/2)-1 peaky-floor(npy/2)-1];
-
-%------------------------------------------------------------------------
-% --- Find the maximum of the correlation function after quadratic interpolation
-function [vector,F] = quadr_fit(result_conv,x,y)
-[npy,npx]=size(result_conv);
-if x<4 || y<4 || npx-x<4 ||npy-y <4
-    F=1;
-    vector=[x y];
-else
-    F=0;
-    x_ind=x-4:x+4;
-    y_ind=y-4:y+4;
-    x_vec=0.25*(x_ind-x);
-    y_vec=0.25*(y_ind-y);
-    [X,Y]=meshgrid(x_vec,y_vec);
-    coord=[reshape(X,[],1) reshape(Y,[],1)];
-    result_conv=reshape(result_conv(y_ind,x_ind),[],1);
-    
-    
-    % n=numel(X);
-    % x=[X Y];
-    % X=X-0.5;
-    % Y=Y+0.5;
-    % y = (X.*X+2*Y.*Y+X.*Y+6) + 0.1*rand(n,1);
-    p = polyfitn(coord,result_conv,2);
-    A(1,1)=2*p.Coefficients(1);
-    A(1,2)=p.Coefficients(2);
-    A(2,1)=p.Coefficients(2);
-    A(2,2)=2*p.Coefficients(4);
-    vector=[x y]'-A\[p.Coefficients(3) p.Coefficients(5)]';
-    vector=vector'-[floor(npx/2) floor(npy/2)]-1 ;
-    % zg = polyvaln(p,coord);
-    % figure
-    % surf(x_vec,y_vec,reshape(zg,9,9))
-    % hold on
-    % plot3(X,Y,reshape(result_conv,9,9),'o')
-    % hold off
-end
-
-
-function FF=detect_false(Param,C,U,V,FFIn)
-FF=FFIn;%default, good vectors
-% FF=1, for correlation max at edge, not set in this function
-% FF=2, for too small correlation
-% FF=3, for velocity outside bounds
-% FF=4 for exclusion by difference with the smoothed field, set by call to function filter_tps
-
-if isfield (Param,'MinCorr')
-     FF(C<Param.MinCorr & FFIn==0)=2;
-end
-if (isfield(Param,'MinVel')&&~isempty(Param.MinVel))||(isfield (Param,'MaxVel')&&~isempty(Param.MaxVel))
-    Umod= U.*U+V.*V;
-    if isfield (Param,'MinVel')&&~isempty(Param.MinVel)
-        U2Min=Param.MinVel*Param.MinVel;
-        FF(Umod<U2Min & FFIn==0)=3;
-    end
-    if isfield (Param,'MaxVel')&&~isempty(Param.MaxVel)
-         U2Max=Param.MaxVel*Param.MaxVel;
-        FF(Umod>U2Max & FFIn==0)=3;
     end
 end
@@ -1378 +955 @@
 end
 
-
-
-
+function FF=detect_false(Param,C,U,V,FFIn)
+FF=FFIn;%default, good vectors
+% FF=1, for correlation max at edge, not set in this function
+% FF=2, for too small correlation
+% FF=3, for velocity outside bounds
+% FF=4 for exclusion by difference with the smoothed field, set by call to function filter_tps
+
+if isfield (Param,'MinCorr')
+     FF(C<Param.MinCorr & FFIn==0)=2;
+end
+if (isfield(Param,'MinVel')&&~isempty(Param.MinVel))||(isfield (Param,'MaxVel')&&~isempty(Param.MaxVel))
+    Umod= U.*U+V.*V;
+    if isfield (Param,'MinVel')&&~isempty(Param.MinVel)
+        U2Min=Param.MinVel*Param.MinVel;
+        FF(Umod<U2Min & FFIn==0)=3;
+    end
+    if isfield (Param,'MaxVel')&&~isempty(Param.MaxVel)
+         U2Max=Param.MaxVel*Param.MaxVel;
+        FF(Umod>U2Max & FFIn==0)=3;
+    end
+end
+
+
+

trunk/src/series/extract_multitif.m

@@ -85 +85 @@
     ParamOut.OutputDirExt='.png';%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
-    ParamOut.CheckOverwriteVisible='on'; % manage the overwrite of existing files (default=1)
+    ParamOut.CheckOverwriteVisible='on'; % manage the overwrite of existing files (default='off')
     ParamOut.CPUTime=10;% expected time for writting the output of one source image ( in minute)
     %% root input file(s) and type
+ 
     % 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
+    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;
+    else
+        msgbox_uvmat('ERROR',['select a multitif file labeled by a number, like im@0001.tif '])
+        return
+    end
     PairString='';
     if isfield(Param.IndexRange,'PairString'); PairString=Param.IndexRange.PairString; end
@@ -99 +105 @@
         msgbox_uvmat('WARNING',['the first input file ' FirstFileName ' does not exist'])
     end
-    ParamOut.NbSlice=Param.IndexRange.MaxIndex_i;
-    %% check the validity of  input file types
-    FileInfo=get_file_info(FirstFileName);
-    if ~strcmp(FileInfo.FileType,'multimage')
+       FileInfo=get_file_info(FirstFileName);
+    if ~strcmp(FileInfo.FileType,'multimage')%check the validity of  input file types
         msgbox_uvmat('ERROR',['invalid file type input: ' FileInfo.FileType ' not a tiff image series'])
         return
     end
+    
+    ParamOut.NbSlice=Param.IndexRange.MaxIndex_i;
+    
     ParamOut.ActionInput.XmlFile=uigetfile_uvmat('pick xml file for timing',fileparts(fileparts(FirstFileName)),'.xml');  
     return

trunk/src/series/extract_rdvision.m

@@ -71 +71 @@
         ParamOut.ProjObject='off';...%can use projection object(option 'off'/'on',
         ParamOut.Mask='off';...%can use mask option   (option 'off'/'on', 'off' by default)
-        ParamOut.CPUTime=0.1;% expected time for writting one image ( in minute)
+         ParamOut.CheckOverwriteVisible='on'; % manage the overwrite of existing files (default=1)
+        ParamOut.CPUTime=0.25;% expected time for writting one image ( in minute)
     ParamOut.OutputDirExt='.extract';%set the output dir extension
     ParamOut.OutputSubDirMode='one'; %output folder given by the folder name of the first input line
@@ -159 +160 @@
             timestamp(ii)=m.Data(ii).timestamp;
         end
+        
+%         
+%         %%%%%BRICOLAGE EXP40
+%         ind1=5*59-10;
+%         ind2=12*59-10;
+%         ind3=119*59-10;
+%         ind4=483*59-10;
+%         timestamp=[timestamp(1:ind4) timestamp(ind4) timestamp(ind4+1:end)];
+%          timestamp=[timestamp(1:ind3) timestamp(ind3) timestamp(ind3+1:end)];
+%           timestamp=[timestamp(1:ind2) timestamp(ind2) timestamp(ind2+1:end)];
+%            timestamp=[timestamp(1:ind1) timestamp(ind1) timestamp(ind1+1:end)];
+%            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+           
+           
         if ParamOut.ActionInput.Createxml
             ParamOut.ActionInput.XmlData.Camera.BurstTiming=time2xmlburst(timestamp,ParamOut.ActionInput.BurstLength);
@@ -174 +189 @@
         difftime=timestamp-timexml;
         %if max(difftime)>0.01
+
         figure
         plot(timestamp,difftime)
@@ -367 +383 @@
 end
 
-[BinList,errormsg]=binread_rdv_series(RootPath,SeqData,m.Data,nbfield2,NomTypeNew,Param.IndexRange.first_i,Param.IndexRange.last_i);
+[BinList,errormsg]=binread_rdv_series(RootPath,SeqData,m.Data,nbfield2,NomTypeNew,Param.IndexRange.first_i,Param.IndexRange.last_i,Param.CheckOverwrite);
 if ~isempty(errormsg)
     disp_uvmat('ERROR',errormsg,checkrun)
@@ -387 +403 @@
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-function [BinList,errormsg]=binread_rdv_series(PathDir,SeqData,SqbData,nbfield2,NomTypeNew,first,last)
+function [BinList,errormsg]=binread_rdv_series(PathDir,SeqData,SqbData,nbfield2,NomTypeNew,first,last,checkoverwrite)
 % BINREAD_RDV Permet de lire les fichiers bin g???n???r???s par Hiris ??? partir du
 % fichier seq associ???.
@@ -442 +458 @@
 bin_file_counter=0;
 %for ii=1:SeqData.nb_frames
+ %%%%%BRICOLAGE EXP40
+        ind1=5*59-10;
+        ind2=12*59-10;
+        ind3=119*59-10;
+        ind4=483*59-10;
+
 for ii=first:last
     j1=[];
+    iinew=ii; 
+%     %%%%%BRICOLAGE EXP40
+%     switch ii
+%         case  ii>=ind1 && ii<ind2
+%       iinew=ii+1;
+%         case ii>=ind2 && ii<ind3
+%         iinew=ii+2;
+%         case  ii>=ind3 && ii<ind4
+%       iinew=ii+3;
+%         case       ii>=ind4 && ii<ind3
+%         iinew=ii+4;    
+%     end
+%      %%%%%%%%%
     if ~isequal(nbfield2,1)
-        j1=mod(ii-1,nbfield2)+1;
-    end
-    i1=floor((ii-1)/nbfield2)+1;
+        j1=mod(iinew-1,nbfield2)+1;
+    end
+    i1=floor((iinew-1)/nbfield2)+1;
+    
     OutputFile=fullfile_uvmat(PathDir,FileDir,'img','.png',NomTypeNew,i1,[],j1);% TODO: set NomTypeNew from SeqData.mode
     fname=fullfile(binrepertoire,sprintf('%s%.5d.bin',SeqData.binfile,SqbData(ii).file_idx));
-    if exist(OutputFile,'file')% do not recreate existing image file
+    if  ~checkoverwrite && exist(OutputFile,'file') % manage the overwrite of existing files (default=1) % manage the overwrite of existing files (default=1)exist(OutputFile,'file')% do not recreate existing image file
         fid=0;
     else

trunk/src/uvmat.m

@@ -2063 +2063 @@
 set(hhseries.ActionName,'Value',index_action);
 series('ActionName_Callback',hObject,[],hhseries); %file input with xml reading  in uvmat, show the image in phys coordinates
+series('ActionInput_Callback',hObject,[],hhseries); %open the menu or options from INPUT in series
 end