Changeset 1164

Timestamp:

Jul 29, 2024, 9:43:17 AM (5 months ago)

Author:

sommeria

Message:

civ3D updated

Location:

trunk/src

Files:

• filter_tps_3D.m (modified) (12 diffs)
• find_file_series.m (modified) (7 diffs)
• get_file_info.m (modified) (1 diff)
• get_file_series.m (modified) (2 diffs)
• series/check_data_files.m (modified) (7 diffs)
• series/civ_3D.m (modified) (40 diffs)
• series/civ_input.m (modified) (3 diffs)
• series/civ_series.m (modified) (13 diffs)
• series/extract_rdvision.m (modified) (2 diffs)
• uvmat.m (modified) (24 diffs)

trunk/src/filter_tps_3D.m

@@ -38 +38 @@
 %=======================================================================
 
-function [SubRange,NbCentre,Coord_tps,U_tps,V_tps,W_tps,U_smooth,V_smooth,W_smooth,FF] =filter_tps(Coord,U,V,W,SubDomainSize,FieldSmooth,Threshold)
+function [SubRange,NbCentre,Coord_tps,U_tps,V_tps,W_tps,U_smooth,V_smooth,W_smooth,FF] =filter_tps_3D(Coord_x,Coord_y,Coord_z,U,V,W,SubDomainSize,FieldSmooth,Threshold)
 
 %% adjust subdomain decomposition
 warning off
-NbVec=size(Coord,1);% nbre of vectors in the field to interpolate
-NbCoord=size(Coord,2);% space dimension,Coord(:,1)= x,Coord(:,2)=  y , Coord(:,3)=  z
+% [npz,npy,npx]=size(Coord_x);
+Coord=[Coord_x Coord_y Coord_z];
+NbVec=numel(Coord_x);% nbre of vectors in the field to interpolate
+NbCoord=3;% space dimension,Coord(:,1)= x,Coord(:,2)=  y , Coord(:,3)=  z
 MinCoord=min(Coord,[],1);%lower coordinate bounds
 MaxCoord=max(Coord,[],1);%upper coordinate bounds
 Range=MaxCoord-MinCoord;%along eacch coordiante x,y,z
-Cellmesh=(prod(Range)/NbVec)^(1/3);
-NbSubDomainX=floor(range(1)/Cellmesh);
-NbSubDomainY=floor(range(2)/Cellmesh);
-NbSubDomainZ=floor(range(3)/Cellmesh);
+Cellmesh=(10*prod(Range)/NbVec)^(1/3);
+NbSubDomainX=ceil(Range(1)/Cellmesh);
+NbSubDomainY=ceil(Range(2)/Cellmesh);
+NbSubDomainZ=ceil(10*Range(3)/Cellmesh);
 
 % NbSubDomain=NbVec/SubDomainSize;% estimated number of subdomains
@@ -56 +58 @@
 % NbSubDomainY=max(floor(sqrt(NbSubDomain*AspectRatio)),1);% estimated number of subdomains in y
 % NbSubDomainZ=max(floor(sqrt(NbSubDomain*AspectRatio)),1);% estimated number of subdomains in y
-NbSubDomain=NbSubDomainX*NbSubDomainY;% new estimated number of subdomains in a matrix shape partition in subdomains
+NbSubDomain=NbSubDomainX*NbSubDomainY*NbSubDomainZ;% new estimated number of subdomains in a matrix shape partition in subdomains
 Siz(1)=Range(1)/NbSubDomainX;%width of subdomains
 Siz(2)=Range(2)/NbSubDomainY;%height of subdomains
@@ -71 +73 @@
 %smoothing=Siz(1)*Siz(2)*FieldSmooth/1000%old calculation before 03 May < r1129
 NbVecSub=NbVec/NbSubDomain;% refined estimation of the nbre of vectors per subdomain
-smoothing=sqrt(Siz(1)*Siz(2)/NbVecSub)*FieldSmooth;%optimum smoothing increase as the typical mesh size =sqrt(SizX*SizY/NbVecSub)^1/2
+smoothing=(Siz(1)*Siz(2)*Siz(3)/NbVecSub)^(1/3)*FieldSmooth;%optimum smoothing increase as the typical mesh size =sqrt(SizX*SizY/NbVecSub)^1/2
 Threshold=Threshold*Threshold;% take the square of the threshold to work with the modulus squared (not done before r1154)
 
@@ -79 +81 @@
 U_tps=zeros(1,NbSubDomain);% initialize  interpolated u component
 V_tps=zeros(1,NbSubDomain);% initialize interpolated v component
+W_tps=zeros(1,NbSubDomain);% initialize interpolated v component
 NbCentre=zeros(1,NbSubDomain);%number of interpolated field values per subdomain, =0 by default
-W_tps=[];%default (2 component case)
 U_smooth=zeros(NbVec,1); % smoothed velocity U at the initial positions
 V_smooth=zeros(NbVec,1);% smoothed velocity V at the initial positions
-W_smooth=[];%default (2 component case)
+W_smooth=zeros(NbVec,1);%default (2 component case)
 FF=false(NbVec,1);%false flag=0 (false) by default
 nb_select=zeros(NbVec,1);
@@ -100 +102 @@
         %check_next=false;% test to go to next iteration with wider subdomain
         ind_sel_previous=ind_sel;% record the set of selected vector indices for next iteration
-        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
+        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)&...
+            Coord(:,3)>=SubRange(3,1,isub) & Coord(:,3)<=SubRange(3,2,isub));% indices of vectors in the subdomain #isub
          check_partial_domain=sum(SubRange(:,1,isub)> MinCoord' | SubRange(:,2,isub)< MaxCoord');
         % isub
@@ -121 +125 @@
             [U_smooth_sub,U_tps_sub]=tps_coeff(Coord(ind_sel,:),U(ind_sel),smoothing);
             [V_smooth_sub,V_tps_sub]=tps_coeff(Coord(ind_sel,:),V(ind_sel),smoothing);
+            [W_smooth_sub,W_tps_sub]=tps_coeff(Coord(ind_sel,:),W(ind_sel),smoothing);
             UDiff=U_smooth_sub-U(ind_sel);% difference between interpolated U component and initial value
             VDiff=V_smooth_sub-V(ind_sel);% difference between interpolated V component and initial value
-            NormDiff=UDiff.*UDiff+VDiff.*VDiff;% Square of difference norm
+            WDiff=W_smooth_sub-W(ind_sel);% difference between interpolated V component and initial value
+            NormDiff=UDiff.*UDiff+VDiff.*VDiff+WDiff.*WDiff;% Square of difference norm
             ind_ind_sel=1:numel(ind_sel);%default
             if exist('Threshold','var')&&~isempty(Threshold)
@@ -134 +140 @@
                 x_width=(SubRange(1,2,isub)-SubRange(1,1,isub))/pi;
                 y_width=(SubRange(2,2,isub)-SubRange(2,1,isub))/pi;
+                z_width=(SubRange(3,2,isub)-SubRange(3,1,isub))/pi;
                 x_dist=(Coord(ind_sel,1)-CentreX(isub))/x_width;% relative x distance to the retangle centre
                 y_dist=(Coord(ind_sel,2)-CentreY(isub))/y_width;% relative ydistance to the retangle centre
-                weight=cos(x_dist).*cos(y_dist);%weighting fct =1 at the rectangle center and 0 at edge
-                %weight=1;% case for r1129 and before
+                z_dist=(Coord(ind_sel,3)-CentreZ(isub))/z_width;% relative ydistance to the retangle centre
+                weight=cos(x_dist).*cos(y_dist).*cos(z_dist);%weighting fct =1 at the rectangle center and 0 at edge
                 U_smooth(ind_sel)=U_smooth(ind_sel)+weight.*U_smooth_sub;
                 V_smooth(ind_sel)=V_smooth(ind_sel)+weight.*V_smooth_sub;
+                W_smooth(ind_sel)=W_smooth(ind_sel)+weight.*W_smooth_sub;
                 NbCentre(isub)=numel(ind_sel);
                 Coord_tps(1:NbCentre(isub),:,isub)=Coord(ind_sel,:);
-                U_tps(1:NbCentre(isub)+3,isub)=U_tps_sub;
-                V_tps(1:NbCentre(isub)+3,isub)=V_tps_sub;
+                U_tps(1:NbCentre(isub)+4,isub)=U_tps_sub;
+                V_tps(1:NbCentre(isub)+4,isub)=V_tps_sub;
+                W_tps(1:NbCentre(isub)+4,isub)=W_tps_sub;
                 nb_select(ind_sel)=nb_select(ind_sel)+weight;
                 display(['tps done with ' num2str(numel(ind_sel)) ' vectors in subdomain # ' num2str(isub)  ' among ' num2str(NbSubDomain)])
@@ -157 +166 @@
                 [U_smooth_sub,U_tps_sub]=tps_coeff(Coord(ind_sel(ind_ind_sel),:),U(ind_sel(ind_ind_sel)),smoothing);
                 [V_smooth_sub,V_tps_sub]=tps_coeff(Coord(ind_sel(ind_ind_sel),:),V(ind_sel(ind_ind_sel)),smoothing);
+                [W_smooth_sub,W_tps_sub]=tps_coeff(Coord(ind_sel(ind_ind_sel),:),W(ind_sel(ind_ind_sel)),smoothing);
                 x_width=(SubRange(1,2,isub)-SubRange(1,1,isub))/pi;
                 y_width=(SubRange(2,2,isub)-SubRange(2,1,isub))/pi;
+                z_width=(SubRange(3,2,isub)-SubRange(3,1,isub))/pi;
                 x_dist=(Coord(ind_sel(ind_ind_sel),1)-CentreX(isub))/x_width;% relative x distance to the retangle centre
                 y_dist=(Coord(ind_sel(ind_ind_sel),2)-CentreY(isub))/y_width;% relative ydistance to the retangle centre
-                weight=cos(x_dist).*cos(y_dist);%weighting fct =1 at the rectangle center and 0 at edge
+                  z_dist=(Coord(ind_sel(ind_ind_sel),3)-CentreZ(isub))/z_width;% relative ydistance to the retangle centre
+                weight=cos(x_dist).*cos(y_dist).*cos(z_dist);%weighting fct =1 at the rectangle center and 0 at edge
                 %weight=1;
                 U_smooth(ind_sel(ind_ind_sel))=U_smooth(ind_sel(ind_ind_sel))+weight.*U_smooth_sub;
                 V_smooth(ind_sel(ind_ind_sel))=V_smooth(ind_sel(ind_ind_sel))+weight.*V_smooth_sub;
+                   W_smooth(ind_sel(ind_ind_sel))=W_smooth(ind_sel(ind_ind_sel))+weight.*W_smooth_sub;
                 NbCentre(isub)=numel(ind_ind_sel);
                 Coord_tps(1:NbCentre(isub),:,isub)=Coord(ind_sel(ind_ind_sel),:);
-                U_tps(1:NbCentre(isub)+3,isub)=U_tps_sub;
-                V_tps(1:NbCentre(isub)+3,isub)=V_tps_sub;
+                U_tps(1:NbCentre(isub)+4,isub)=U_tps_sub;
+                V_tps(1:NbCentre(isub)+4,isub)=V_tps_sub;
+                W_tps(1:NbCentre(isub)+4,isub)=W_tps_sub;
                 nb_select(ind_sel(ind_ind_sel))=nb_select(ind_sel(ind_ind_sel))+weight;
                 display(['tps redone with ' num2str(numel(ind_sel)) ' vectors after elimination of ' num2str(numel(ind_sel)-numel(ind_ind_sel)) ' erratic vectors in subdomain # ' num2str(isub) ' among ' num2str(NbSubDomain)])
@@ -189 +203 @@
     U_tps(:,ind_empty)=[];
     V_tps(:,ind_empty)=[];
+    W_tps(:,ind_empty)=[];
     NbCentre(ind_empty)=[];
 end
@@ -196 +211 @@
 U_smooth=U_smooth./nb_select;% take the average at the intersection of several subdomains
 V_smooth=V_smooth./nb_select;
+W_smooth=W_smooth./nb_select;
 
 %eliminate the vectors with diff>threshold not yet eliminated
@@ -201 +217 @@
 UDiff=U_smooth-U;% difference between interpolated U component and initial value
 VDiff=V_smooth-V;% difference between interpolated V component and initial value
-NormDiff=UDiff.*UDiff+VDiff.*VDiff;% Square of difference norm
+WDiff=W_smooth-W;% difference between interpolated V component and initial value
+NormDiff=UDiff.*UDiff+VDiff.*VDiff+WDiff.*WDiff;% Square of difference norm
 FF(NormDiff>Threshold)=true;%put FF value to 1 to identify the criterium of elimmination
 end
@@ -207 +224 @@
 U_smooth(FF)=U(FF);% set to the initial values the eliminated vectors (flagged as false)
 V_smooth(FF)=V(FF);
+W_smooth(FF)=W(FF);
 fill=zeros(NbCoord+1,NbCoord,size(SubRange,3)); %matrix of zeros to complement the matrix Data.Civ1_Coord_tps (conveninent for file storage)
 Coord_tps=cat(1,Coord_tps,fill);
 
+

trunk/src/find_file_series.m

@@ -49 +49 @@
 
 %% get input root name and info on the input file
-if isempty(regexp(FilePath,'^http://','once'))
+if isempty(regexp(FilePath,'^http://','once'))% case of usual file input
 fullfileinput=fullfile(FilePath,fileinput);% input file name with path
 else
-  fullfileinput=[FilePath '/' fileinput];
+  fullfileinput=[FilePath '/' fileinput]; % case of web input
 end
 [FileInfo,MovieObject]=get_file_info(fullfileinput);
@@ -105 +105 @@
         sep2='';
         i1_str='(?<i1>)';%will set i1=[];
-        i1_star='';
         i2_str='(?<i2>)';%will set i2=[];
-        i2_star='';
         j1_str='(?<j1>)';%will set j1=[];
-        j1_star='';
         j2_str='(?<j2>)';%will set j2=[];
-        j2_star='';
+
         %Look for cases with letter indexing for the second index
         r=regexp(NomType,'^(?<sep1>_?)(?<i1>\d+)(?<sep2>_?)(?<j1>[a|A])(?<j2>[b|B]?)$','names');
@@ -118 +115 @@
             sep2=r.sep2;
             i1_str='(?<i1>\d+)';
-            i1_star='*';
             if strcmp(lower(r.j1),r.j1)% lower case index
                 j1_str='(?<j1>[a-z])';
@@ -124 +120 @@
                 j1_str='(?<j1>[A-Z])'; % upper case index
             end
-            j1_star='*';
             if ~isempty(r.j2)
                 if strcmp(lower(r.j1),r.j1)
@@ -131 +126 @@
                     j2_str='(?<j2>[A-Z])';
                 end
-                j2_star='*';
             end
         else %numerical indexing
@@ -138 +132 @@
                 sep1=r.sep1;
                 i1_str='(?<i1>\d+)';
-                i1_star='*';
                 if ~isempty(r.i2)
                     i2_str='(?<i2>-\d+)';
-                    i2_star='-*';
                 end
                 if ~isempty(r.j1)
                     j1_str='(?<j1>_\d+)';
-                    j1_star='_*';
                 end
                 if ~isempty(r.j2)
                     j2_str='(?<j2>-\d+)';
-                    j2_star='-*';
                 end
             end
@@ -290 +280 @@
             i1_series=i1_series(:,2)*ones(1,FileInfo.NumberOfFrames);%
             i1_series=[zeros(size(i1_series,1),1) i1_series];
+            if ~isempty(i2_series)
+                i2_series=i2_series(:,2)*ones(1,FileInfo.NumberOfFrames);%
+            i2_series=[zeros(size(i2_series,1),1) i2_series];
+            end
             j1_series=ones(size(i1_series,1),1)*(1:FileInfo.NumberOfFrames);%
             j1_series=[zeros(size(i1_series,1),1) j1_series];

trunk/src/get_file_info.m

@@ -153 +153 @@
                         [Data,tild,tild,errormsg]=nc2struct(fileinput,[]);
                         if isempty(errormsg)
-                            if isfield(Data,'absolut_time_T0') && isfield(Data,'hart') && ~isempty(Data.absolut_time_T0) && ~isempty(Data.hart)
-                                FileInfo.FileType='civx';%old civ data from the Fortran program
-                                if isfield(Data,'patch2') && isequal(Data.patch2,1)
-                                    FileInfo.CivStage=6;
-                                elseif isfield(Data,'fix2') && isequal(Data.fix2,1)
-                                    FileInfo.CivStage=5;
-                                elseif  isfield(Data,'civ2')&& isequal(Data.civ2,1)
-                                    FileInfo.CivStage=4;
-                                elseif isfield(Data,'patch')&&isequal(Data.patch,1)
-                                    FileInfo.CivStage=3;
-                                elseif isfield(Data,'fix')&&isequal(Data.fix,1)
-                                    FileInfo.CivStage=2;
-                                else
-                                    FileInfo.CivStage=1;
-                                end
-                            elseif isfield(Data,'Conventions') && strcmp(Data.Conventions,'uvmat/civdata')
+                            %                             if isfield(Data,'absolut_time_T0') && isfield(Data,'hart') && ~isempty(Data.absolut_time_T0) && ~isempty(Data.hart)
+                            %                                 FileInfo.FileType='civx';%old civ data from the Fortran program
+                            %                                 if isfield(Data,'patch2') && isequal(Data.patch2,1)
+                            %                                     FileInfo.CivStage=6;
+                            %                                 elseif isfield(Data,'fix2') && isequal(Data.fix2,1)
+                            %                                     FileInfo.CivStage=5;
+                            %                                 elseif  isfield(Data,'civ2')&& isequal(Data.civ2,1)
+                            %                                     FileInfo.CivStage=4;
+                            %                                 elseif isfield(Data,'patch')&&isequal(Data.patch,1)
+                            %                                     FileInfo.CivStage=3;
+                            %                                 elseif isfield(Data,'fix')&&isequal(Data.fix,1)
+                            %                                     FileInfo.CivStage=2;
+                            %                                 else
+                            %                                     FileInfo.CivStage=1;
+                            %                                 end
+                            %                             else
+                            if isfield(Data,'Conventions') && strcmp(Data.Conventions,'uvmat/civdata')
                                 FileInfo.FileType='civdata'; % test for civ velocity fields
                                 FileInfo.CivStage=Data.CivStage;
+                                MaskFile='';
+                                if isfield(Data,'Civ2_Mask')
+                                    MaskFile=Data.Civ2_Mask;
+                                    if isfield(Data,'Civ2_NbSlice')
+                                        FileInfo.MaskNbSlice=Data.Civ2_NbSlice;
+                                    end
+                                elseif isfield(Data,'Civ1_Mask')
+                                    MaskFile=Data.Civ1_Mask;
+                                    if isfield(Data,'Civ1_NbSlice')
+                                        FileInfo.MaskNbSlice=Data.Civ1_NbSlice;
+                                    end
+                                end
+                                if isfield(Data,'VolumeScan')
+                                    FileInfo.VolumeScan=Data.VolumeScan;
+                                end
+                                if ~isempty(MaskFile)
+                                    [RootPath,SubDir,RootFile,~,~,~,~,FileExt,NomType]=fileparts_uvmat(MaskFile);
+                                    if strcmp(NomType,'_1')&& isfield(FileInfo,'MaskNbSlice')
+                                        FileInfo.MaskFile=fullfile(RootPath,SubDir,RootFile);
+                                    else
+                                        FileInfo.MaskFile=MaskFile;% single mask for the series (no indexing)
+                                    end
+                                    FileInfo.MaskExt=FileExt;
+                                end
                             elseif isfield(Data,'Conventions') && strcmp(Data.Conventions,'uvmat/civdata_3D')
                                 FileInfo.FileType='civdata_3D'; % test for 3D volume civ velocity fields

trunk/src/get_file_series.m

@@ -84 +84 @@
         FilePath=fullfile(InputTable{iview,1},InputTable{iview,2});
         fileinput=[InputTable{iview,3} InputTable{iview,4} InputTable{iview,5}];
-        [tild,tild,tild,i1_series{iview},i2_series{iview},j1_series{iview},j2_series{iview},NomType,FileInfo,MovieObject,...
-            i1_input,i2_input,j1_input,j2_input]=find_file_series(FilePath,fileinput);
+        [tild,tild,tild,i1_series{iview},i2_series{iview},j1_series{iview},j2_series{iview}]=find_file_series(FilePath,fileinput);
         i1_series{iview}=squeeze(i1_series{iview}(1,:,:)); %select first  pair index as ordered by find_file_series
         i2_series{iview}=squeeze(i2_series{iview}(1,:,:)); %select first  pair index as ordered by find_file_series
@@ -137 +136 @@
         [i1_series{iview},i2_series{iview},j1_series{iview},j2_series{iview}]=find_file_indices(ref_i,ref_j,str2num(r.num1),str2num(r.num2),r.mode);
     end
-    %     if ~isequal(r(1).mode,'*-*')% imposed pairs or single i and/or j index
-    %         [i1_series{iview},i2_series{iview},j1_series{iview},j2_series{iview}]=find_file_indices(ref_i,ref_j,str2num(r.num1),str2num(r.num2),r.mode);
-    %     end
     
     %list of files

trunk/src/series/check_data_files.m

@@ -63 +63 @@
     ParamOut.ProjObject='off';%can use projection object(option 'off'/'on',
     ParamOut.Mask='off';%can use mask option   (option 'off'/'on', 'off' by default)
-    ParamOut.OutputSubDirMode='none'; %(options 'none'/'custom'/'auto'/'first'/'last','auto' by default)
+    ParamOut.OutputSubDirMode='auto'; %(options 'none'/'custom'/'auto'/'first'/'last','auto' by default)
     %                      'none' =no output files
+    ParamOut.OutputDirExt='.check_series';%set the output dir extension
     return
 end
@@ -120 +121 @@
 for iview=1:nbview
     if isequal(FileType{iview},'mmreader')||isequal(FileType{iview},'video')||isequal(FileType{iview},'multimage')
-        [FileInfo]=get_file_info(filecell{iview,1});
+        [FileInfo]=get_file_info(filecell{iview,1});% get infos on the first file in the current line
         Tabchar{1}=filecell{iview,1};%info.Filename;
         Tabchar{2}='';
@@ -128 +129 @@
         message='';
     else
-        Tabchar={};
+        Tabchar=cell(nbfield_i+1,NbSlice);
         %LOOP ON SLICES
         for i_slice=1:NbSlice
             index_slice=i_slice:NbSlice:nbfield;
-            filefound={};
-            datnum=zeros(1,nbfield_j);
+            filefound=cell(nbfield_i);
+            datnum=zeros(1,nbfield_i);
+            Tabchar(1,i_slice)={['slice #' num2str(i_slice)]};
             for ifile=1:nbfield_i
-                update_waitbar(WaitbarHandle,ifile/nbfield_i)
-                if ishandle(RUNHandle) && ~strcmp(get(RUNHandle,'BusyAction'),'queue')
-                    disp('program stopped by user')
-                    break
-                end
+ %               update_waitbar(WaitbarHandle,ifile/nbfield_i)
+%                 if ishandle(RUNHandle) && ~strcmp(get(RUNHandle,'BusyAction'),'queue')
+%                     disp('program stopped by user')
+%                     break
+%                 end
                 file=filecell{iview,index_slice(ifile)};
                 [Path,Name,ext]=fileparts(file);
@@ -149 +151 @@
                     if isfield(datfile,'datenum')
                         datnum(ifile)=datfile.datenum;
-                        filefound(ifile)={datfile.name};
+                        filefound{ifile}=datfile.name;
                     end
                     lastfield='';
-                    [FileInfo,Object]=get_file_info(file);
+                    [FileInfo]=get_file_info(file);
                     FileType{iview}=FileInfo.FileType;
                     if strcmp(FileType{iview},'civx')||strcmp(FileType{iview},'civdata')
@@ -159 +161 @@
                             stagechoice=1+mod(FileInfo.CivStage-1,3);
                             iter=1+floor((FileInfo.CivStage-1)/3);
-                            lastfield=[liststage{stagechoice} num2str(iter)];
-                            %liststage={'civ1','fix1','patch1','civ2','fix2','patch2'};                        
-                            %lastfield=liststage{FileInfo.CivStage};                           
+                            lastfield=[liststage{stagechoice} num2str(iter)];                         
                         end
                     end
                     lastfield=[FileType{iview} ', ' lastfield];
                 end
-                Tabchar(1,i_slice)={['slice #' num2str(i_slice)]};
                 Tabchar(ifile+1,i_slice)={[file '...' lastfield]};
             end
@@ -177 +176 @@
             end
         else
-            datnum=datnum(find(datnum));%keep the non zero values corresponding to existing files
-            filefound=filefound(find(datnum));
+            datnum=datnum(datnum);%keep the non zero values corresponding to existing files
+            filefound=filefound(datnum);
             [first,ind]=min(datnum);
             [last,indlast]=max(datnum);
@@ -188 +187 @@
         end
     end
-    hfig=figure(iview);
-    clf
-    if iview>1
-        pos=get(iview-1,'Position');
-        pos(1)=pos(1)+(iview-1)*pos(1)/nbview;
-        set(hfig,'Position',pos)
+    if checkrun
+        hfig=figure(iview);
+        clf
+        if iview>1
+            pos=get(iview-1,'Position');
+            pos(1)=pos(1)+(iview-1)*pos(1)/nbview;
+            set(hfig,'Position',pos)
+        end
+        set(hfig,'name',['check_data_files:view= ' num2str(iview)])
+        set(hfig,'MenuBar','none')% suppress the menu bar
+        set(hfig,'NumberTitle','off')%suppress the fig number in the title
+        h=uicontrol('Style','listbox', 'Position', [20 20 500 300], 'String', Tabchar, 'Callback', {'open_uvmat'});
+        hh=uicontrol('Style','listbox', 'Position', [20 340 500 40], 'String', message);
+    else % show results in log file
+        disp(Tabchar)
+        disp(message)
     end
-    set(hfig,'name',['check_data_files:view= ' num2str(iview)])
-    set(hfig,'MenuBar','none')% suppress the menu bar
-    set(hfig,'NumberTitle','off')%suppress the fig number in the title
-    h=uicontrol('Style','listbox', 'Position', [20 20 500 300], 'String', Tabchar, 'Callback', {'open_uvmat'});
-    hh=uicontrol('Style','listbox', 'Position', [20 340 500 40], 'String', message);
-end
-
+end
+'END'
 % 'open_uvmat': open with uvmat the  field selected in the list of 'series/check_data_files'
 %------------------------------------------------------------------------

trunk/src/series/civ_3D.m

@@ -1 @@
-%'civ_3D': 3D PIV from image scan in a volume 
+%'civ_3D': 3D PIV from image scan in a volume
 
 %------------------------------------------------------------------------
@@ -12 +12 @@
 % Param: Matlab structure of input  parameters
 %     Param contains info of the GUI series using the fct read_GUI.
-%     Param.Action.RUN = 0 (to set the status of the GUI series) or =1 to RUN the computation 
+%     Param.Action.RUN = 0 (to set the status of the GUI series) or =1 to RUN the computation
 %     Param.InputTable: sets the input file(s)
 %           if absent, the fct looks for input data in Param.ActionInput     (test mode)
@@ -49 +49 @@
     path_series=fileparts(which('series'));
     addpath(fullfile(path_series,'series'))
-
-   Data=civ_input(Param)
-
+    
+    Data=civ_input(Param)
+    
     Data.Program=mfilename;%gives the name of the current function
     Data.AllowInputSort='off';% allow alphabetic sorting of the list of input file SubDir (options 'off'/'on', 'off' by default)
@@ -83 +83 @@
         end
     end
-
-
+    
+    
     return
 end
@@ -114 +114 @@
 if isfield(Param,'OutputSubDir')&& isfield(Param,'OutputDirExt')
     OutputDir=[Param.OutputSubDir Param.OutputDirExt];
-     OutputPath=fullfile(Param.OutputPath,num2str(Param.Experiment),num2str(Param.Device));
+    OutputPath=fullfile(Param.OutputPath,num2str(Param.Experiment),num2str(Param.Device));
 else
     disp_uvmat('ERROR','no output folder defined',checkrun)
@@ -129 +129 @@
 end
 
-[~,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 )
-if Param.ActionInput.CheckCiv1
-    iview_A=1;% usual PIV, the image series is on the first line of the table
-elseif Param.ActionInput.CheckCiv2 % civ2 is performed without Civ1, a netcdf file series is needed in the first table line
-    iview_A=2;% the second line is used for the input images of Civ2
-end
-if iview_A~=0
-    RootPath_A=Param.InputTable{iview_A,1};
-    RootFile_A=Param.InputTable{iview_A,3};
-    SubDir_A=Param.InputTable{iview_A,2};
-    NomType_A=Param.InputTable{iview_A,4};
-    FileExt_A=Param.InputTable{iview_A,5};
-    if iview_B==0
-        iview_B=iview_A;% the second image series is the same as the first
-    end
-    RootPath_B=Param.InputTable{iview_B,1};
-    RootFile_B=Param.InputTable{iview_B,3};
-    SubDir_B=Param.InputTable{iview_B,2};
-    NomType_B=Param.InputTable{iview_B,4};
-    FileExt_B=Param.InputTable{iview_B,5};
-end
+[filecell,i1_series,i2_series,j1_series,j2_series]=get_file_series(Param);
+
+if size(Param.InputTable,1)==2% netcdf file as first input line
+    RootPath_nc=Param.InputTable{1,1};
+    RootFile_nc=Param.InputTable{1,3};
+    SubDir_nc=Param.InputTable{1,2};
+    NomType_nc=Param.InputTable{1,4};
+    FileExt_nc=Param.InputTable{1,5};
+    iview_A=2;%second line used for image input
+    iview_B=2;
+else
+    iview_A=1;%second line used for image input
+    iview_B=1;
+end
+RootPath_A=Param.InputTable{iview_A,1};
+RootFile_A=Param.InputTable{iview_A,3};
+SubDir_A=Param.InputTable{iview_A,2};
+NomType_A=Param.InputTable{iview_A,4};
+FileExt_A=Param.InputTable{iview_A,5};
+RootPath_B=Param.InputTable{iview_B,1};
+RootFile_B=Param.InputTable{iview_B,3};
+SubDir_B=Param.InputTable{iview_B,2};
+NomType_B=Param.InputTable{iview_B,4};
+FileExt_B=Param.InputTable{iview_B,5};
+
 
 PairCiv1=Param.ActionInput.PairIndices.ListPairCiv1;
 
 [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);
-                    
+    find_pair_indices(PairCiv1,i1_series{1},j1_series{1},MinIndex_i,MaxIndex_i,MinIndex_j,MaxIndex_j);
+
 if isempty(i1_series_Civ1)
     disp_uvmat('ERROR','no image pair for civ in the input file index range',checkrun)
@@ -173 +175 @@
 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
-Data.ListGlobalAttribute=[ListGlobalAttribute Civ1_param];
+if Param.ActionInput.CheckCiv1
+    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
+    Data.ListGlobalAttribute=[ListGlobalAttribute Civ1_param];
+end
 % set the list of variables
 Data.ListVarName={'Coord_z','Civ1_X','Civ1_Y','Civ1_U','Civ1_V','Civ1_W','Civ1_C','Civ1_FF'};%  cell array containing the names of the fields to record
@@ -224 +228 @@
     CheckOverwrite=Param.CheckOverwrite;
 end
-   Data.Civ1_ImageA=fullfile_uvmat(RootPath_A,SubDir_A,RootFile_A,FileExt_A,NomType_A,i1_series_Civ1(1),[],j1_series_Civ1(1,1));
-   Data.Civ1_ImageB=fullfile_uvmat(RootPath_B,SubDir_B,RootFile_B,FileExt_B,NomType_B,i1_series_Civ1(1),[],j2_series_Civ1(1,1));
-    FileInfo=get_file_info(Data.Civ1_ImageA);
-        par_civ1=Param.ActionInput.Civ1;% parameters for civ1
-    par_civ1.ImageHeight=FileInfo.Height;npy=FileInfo.Height;
-    par_civ1.ImageWidth=FileInfo.Width;npx=FileInfo.Width;
+Data.Civ1_ImageA=fullfile_uvmat(RootPath_A,SubDir_A,RootFile_A,FileExt_A,NomType_A,i1_series_Civ1(1),[],j1_series_Civ1(1,1));
+Data.Civ1_ImageB=fullfile_uvmat(RootPath_B,SubDir_B,RootFile_B,FileExt_B,NomType_B,i1_series_Civ1(1),[],j2_series_Civ1(1,1));
+FileInfo=get_file_info(Data.Civ1_ImageA);
+par_civ1=Param.ActionInput.Civ1;% parameters for civ1
+par_civ1.ImageHeight=FileInfo.Height;npy=FileInfo.Height;
+par_civ1.ImageWidth=FileInfo.Width;npx=FileInfo.Width;
 SearchRange_z=Param.ActionInput.Civ1.SearchRange(3);
-    par_civ1.Dz=Param.ActionInput.Civ1.Dz;
-    par_civ1.ImageA=zeros(2*SearchRange_z+1,npy,npx);
+par_civ1.Dz=Param.ActionInput.Civ1.Dz;
+par_civ1.ImageA=zeros(2*SearchRange_z+1,npy,npx);
 par_civ1.ImageB=zeros(2*SearchRange_z+1,npy,npx);
 
@@ -268 +272 @@
         j2=j2_series_Civ1(ifield);
     end
-
-    Data.Civ1_Time=(Time(i2+1,j2+1)+Time(i1+1,j1+1))/2;% the Time is the Time at the middle of the image pair
-    Data.Civ1_Dt=Time(i2+1,j2+1)-Time(i1+1,j1+1);
-
-    for ilist=1:length(list_param)
-        Data.(Civ1_param{4+ilist})=Param.ActionInput.Civ1.(list_param{ilist});
-    end
-    
-    Data.CivStage=1;
-  
+    
+    
+    
     ncfile_out=fullfile_uvmat(OutputPath,'',Param.InputTable{1,3},'.nc',...
-                '_1-1',i1_series_Civ1(ifield),i2_series_Civ1(ifield)); %output file name
-
+        '_1-1',i1_series_Civ1(ifield),i2_series_Civ1(ifield)); %output file name
+    
     if ~CheckOverwrite && exist(ncfile_out,'file')
         disp(['existing output file ' ncfile_out ' already exists, skip to next field'])
         continue% skip iteration if the mode overwrite is desactivated and the result file already exists
     end
-
-
-    %% Civ1
-    % if Civ1 computation is requested
-    if isfield (Param.ActionInput,'Civ1')
-
+    
+    
+%% Civ1 computation if requested
+    if Param.ActionInput.CheckCiv1
+        
         disp('civ1 started')
-
+        Data.Civ1_Time=(Time(i2+1,j2+1)+Time(i1+1,j1+1))/2;% the Time is the Time at the middle of the image pair
+        Data.Civ1_Dt=Time(i2+1,j2+1)-Time(i1+1,j1+1);
+        
+        for ilist=1:length(list_param)
+            Data.(Civ1_param{4+ilist})=Param.ActionInput.Civ1.(list_param{ilist});
+        end
+        
+        Data.CivStage=1;
         % 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
@@ -307 +310 @@
             par_civ1.ImageB(iz,:,:) = B;
         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,:,:), ~, errormsg] = civ3D (par_civ1);
-        if ~isempty(errormsg)
-            disp_uvmat('ERROR',errormsg,checkrun)
-            return
-        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);%            
-                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;
-            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
-        end
-        Data.Civ1_V=-Data.Civ1_V;%reverse v
-        Data.Civ1_Y=npy-Data.Civ1_Y+1;%reverse y
-        [npz,npy,npx]=size(Data.Civ1_X);
-        
-      %  Data.Coord_y=flip(1:npy);
-       % Data.Coord_x=1:npx;
+        
         % case of mask TO ADAPT
         if par_civ1.CheckMask&&~isempty(par_civ1.Mask)
@@ -377 +340 @@
             end
         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,:,:), ~, errormsg] = civ3D (par_civ1);
+        if ~isempty(errormsg)
+            disp_uvmat('ERROR',errormsg,checkrun)
+            return
+        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);%
+                    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;
+            end
+            % caluclate velocity data
+            [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
+        end
+        %         Data.Civ1_V=-Data.Civ1_V;%reverse v
+        %         Data.Civ1_Y=npy-Data.Civ1_Y+1;%reverse y
+        
+        
+        %  Data.Coord_y=flip(1:npy);
+        % Data.Coord_x=1:npx;
+        
+    else
+        Data=nc2struct(filecell{1,ifield});%read civ1 and fix1 data in the existing netcdf file
+        
         % Data.ListVarName=[Data.ListVarName 'Civ1_Z'];
         % Data.Civ1_X=[];Data.Civ1_Y=[];Data.Civ1_Z=[];
         % Data.Civ1_U=[];Data.Civ1_V=[];Data.Civ1_C=[];
-        % 
-        % 
+        %
+        %
         % Data.Civ1_X=[Data.Civ1_X reshape(xtable,[],1)];
         % Data.Civ1_Y=[Data.Civ1_Y reshape(Param.Civ1.ImageHeight-ytable+1,[],1)];
@@ -389 +398 @@
         % Data.Civ1_C=[Data.Civ1_C reshape(ctable,[],1)];
         % Data.Civ1_FF=[Data.Civ1_FF reshape(F,[],1)];
-
-    end
-
-    %% Fix1
+        
+    end
+    
+%% Fix1
     if isfield (Param.ActionInput,'Fix1')
         disp('detect_false1 started')
@@ -412 +421 @@
         Data.CivStage=2;
     end
-    %% Patch1
+    
+%% Patch1
     if isfield (Param.ActionInput,'Patch1')
         disp('patch1 started')
         tstart_patch1=tic;
-
+        
         % record the processing parameters of Patch1 as global attributes in the result nc file
         list_param=fieldnames(Param.ActionInput.Patch1)';
@@ -426 +436 @@
         Data.CivStage=3;% record the new state of processing
         Data.ListGlobalAttribute=[Data.ListGlobalAttribute Patch1_param];
-
+        
         % list the variables to record
         nbvar=length(Data.ListVarName);
@@ -433 +443 @@
         Data.VarAttribute{nbvar+1}.Role='vector_x';
         Data.VarAttribute{nbvar+2}.Role='vector_y';
-        Data.VarAttribute{nbvar+5}.Role='vector_z';      
+        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; 
+        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);
@@ -447 +457 @@
             return
         end
-
+        
         % perform Patch calculation using the UVMAT fct 'filter_tps'
-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
+        Civ1_Z=0.5*Data.Civ1_W;
+        for iz=1:numel(Data.Coord_z)
+            Civ1_Z(iz,:,:)=Civ1_Z(iz,:,:)+Data.Coord_z(iz);    
+        end
+        [npz,npy,npx]=size(Data.Civ1_X);
         [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);
+            Data.Civ1_U_smooth(ind_good),Data.Civ1_V_smooth(ind_good),Data.Civ1_W_smooth(ind_good),FFres]=...
+            filter_tps_3D(Data.Civ1_X(ind_good),Data.Civ1_Y(ind_good),Civ1_Z(ind_good),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(4*FFres);
+        Data.Civ1_FF=reshape(Data.Civ1_FF,npz,npy,npx);
+        Data.Civ1_U=reshape(Data.Civ1_U,npz,npy,npx);
+        Data.Civ1_V=reshape(Data.Civ1_V,npz,npy,npx);
+        Data.Civ1_W=reshape(Data.Civ1_W,npz,npy,npx);
         time_patch1=toc(tstart_patch1);
         disp('patch1 performed')
     end
-
+    
     %% Civ2
     if isfield (Param.ActionInput,'Civ2')
@@ -503 +519 @@
         end
         % end
-
+        
         % get the guess from patch1 or patch2 (case 'CheckCiv3')
-
+        
         if isfield (par_civ2,'CheckCiv3') && par_civ2.CheckCiv3 %get the guess from  patch2
             SubRange= Data.Civ2_SubRange;
@@ -530 +546 @@
             Data.CivStage=4;
         end
-
+        
         Shiftx=zeros(size(par_civ2.Grid,1),1);% shift expected from civ1 data
         Shifty=zeros(size(par_civ2.Grid,1),1);
         nbval=zeros(size(par_civ2.Grid,1),1);% nbre of interpolated values at each grid point (from the different patch subdomains)
-
+        
         NbSubDomain=size(SubRange,3);
         for isub=1:NbSubDomain% for each sub-domain of Patch1
@@ -577 +593 @@
             end
         end
-
-
+        
+        
         if strcmp(Param.ActionInput.ListCompareMode,'displacement')
             Civ1_Dt=1;
@@ -585 +601 @@
             Civ2_Dt=Time(i2_civ2+1,j2_civ2+1)-Time(i1_civ2+1,j1_civ2+1);
         end
-
+        
         par_civ2.SearchBoxShift=(Civ2_Dt/Civ1_Dt)*[Shiftx(nbval>=1)./nbval(nbval>=1) Shifty(nbval>=1)./nbval(nbval>=1)];
         % shift the grid points by half the expected shift to provide the correlation box position in image A
@@ -595 +611 @@
             par_civ2.DVDY=DVDY(nbval>=1)./nbval(nbval>=1);
         end
-
+        
         % calculate velocity data (y and v in image indices, reverse to y component)
-
+        
         [xtable, ytable, utable, vtable, ctable, F,result_conv,errormsg] = civ (par_civ2);
-
+        
         list_param=(fieldnames(Param.ActionInput.Civ2))';
         list_param(strcmp('TestCiv2',list_param))=[];% remove the parameter TestCiv2 from the list
@@ -620 +636 @@
         end
         Data.ListGlobalAttribute=[Data.ListGlobalAttribute Civ2_param];
-
+        
         nbvar=numel(Data.ListVarName);
         % define the Civ2 variable (if Civ2 data are not replaced from previous calculation)
@@ -651 +667 @@
         end
     end
-
+    
     %% Fix2
     if isfield (Param.ActionInput,'Fix2')
@@ -665 +681 @@
         Data.CivStage=Data.CivStage+1;
     end
-
+    
     %% Patch2
     if isfield (Param.ActionInput,'Patch2')
@@ -678 +694 @@
         end
         Data.ListGlobalAttribute=[Data.ListGlobalAttribute Patch2_param];
-
+        
         nbvar=length(Data.ListVarName);
         Data.ListVarName=[Data.ListVarName {'Civ2_U_smooth','Civ2_V_smooth','Civ2_SubRange','Civ2_NbCentres','Civ2_Coord_tps','Civ2_U_tps','Civ2_V_tps'}];
         Data.VarDimName=[Data.VarDimName {'nb_vec_2','nb_vec_2',{'nb_coord','nb_bounds','nb_subdomain_2'},{'nb_subdomain_2'},...
             {'nb_tps_2','nb_coord','nb_subdomain_2'},{'nb_tps_2','nb_subdomain_2'},{'nb_tps_2','nb_subdomain_2'}}];
-
+        
         Data.VarAttribute{nbvar+1}.Role='vector_x';
         Data.VarAttribute{nbvar+2}.Role='vector_y';
@@ -700 +716 @@
             return
         end
-
+        
         [Data.Civ2_SubRange,Data.Civ2_NbCentres,Data.Civ2_Coord_tps,Data.Civ2_U_tps,Data.Civ2_V_tps,tild,Ures,Vres,tild,FFres]=...
             filter_tps([Data.Civ2_X(ind_good) Data.Civ2_Y(ind_good)],Data.Civ2_U(ind_good),Data.Civ2_V(ind_good),[],Data.Patch2_SubDomainSize,Data.Patch2_FieldSmooth,Data.Patch2_MaxDiff);
@@ -710 +726 @@
         disp('patch2 performed')
     end
-
+    
     %% write result in a netcdf file if requested
     % if CheckOutputFile
@@ -726 +742 @@
     disp(['time civ2 ' num2str(time_civ2,2) ' s'])
     disp(['time patch2 ' num2str(time_patch2,2) ' s'])
-
+    
 end
 
@@ -747 +763 @@
 %  .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
+%  .SearchRange:  1,2 vector giving the range of the search 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
@@ -764 +780 @@
 function [xtable,ytable,utable,vtable,wtable,ctable,FF,result_conv,errormsg] = civ3D (par_civ)
 %% check image sizes
-[npz,npy_ima,npx_ima]=size(par_civ.ImageA);% npz = 
+[npz,npy_ima,npx_ima]=size(par_civ.ImageA);% npz =
 if ~isequal(size(par_civ.ImageB),[npz npy_ima npx_ima])
     errormsg='image pair with unequal size';
@@ -786 +802 @@
 ibx2=floor(par_civ.CorrBoxSize(1)/2);
 iby2=floor(par_civ.CorrBoxSize(2)/2);
-isx2=par_civ.SearchRange(1);
-isy2=par_civ.SearchRange(2);
+isx2=par_civ.SearchRange(1)+ibx2;
+isy2=par_civ.SearchRange(2)+iby2;
 isz2=par_civ.SearchRange(3);
 kref=isz2+1;%middle index of the z slice
@@ -878 +894 @@
                 image2_mean=mean(mean(image2_crop,2),3);
             end
-     
+            
             if FF(ivec_y,ivec_x)==0
                 if check_MinIma && (image1_mean < par_civ.MinIma || image2_mean < par_civ.MinIma)
                     FF(ivec_y,ivec_x)=1;       %threshold on image minimum
-                
+                    
                 elseif check_MaxIma && (image1_mean > par_civ.MaxIma || image2_mean > par_civ.MaxIma)
                     FF(ivec_y,ivec_x)=1;    %threshold on image maximum
                 end
                 if FF(ivec_y,ivec_x)==1
-                    utable(ivec_y,ivec_x)=NaN;  
+                    utable(ivec_y,ivec_x)=NaN;
                     vtable(ivec_y,ivec_x)=NaN;
                 else
@@ -898 +914 @@
                     %     % image2_crop=(image2_crop-image2_mean);
                     % end
-
+                    
                     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 
+                    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
                     subima1=squeeze(image1_crop(kref,:,:))-image1_mean(kref);
                     for kz=1:npz
-                        subima2=squeeze(image2_crop(kz,:,:))-image2_mean(kz);   
+                        subima2=squeeze(image2_crop(kz,:,:))-image2_mean(kz);
                         correl_xy=conv2(subima2,flip(flip(subima1,2),1),'valid');
                         result_conv(kz,:,:)=correl_xy;
-                        max_xy(kz)=max(max(correl_xy));             
-                         [yk(kz),xk(kz)]=find(correl_xy==max_xy(kz),1);%find the indices of the maximum, use 0.999 to avoid rounding errors
+                        max_xy(kz)=max(max(correl_xy));
+                        [yk(kz),xk(kz)]=find(correl_xy==max_xy(kz),1);%find the indices of the maximum, use 0.999 to avoid rounding errors
                     end
                     [corrmax,dz]=max(max_xy);
@@ -932 +948 @@
                         end
                         utable(ivec_y,ivec_x)=vector(1)+shiftx(ivec_y,ivec_x);
-                        vtable(ivec_y,ivec_x)=vector(2)+shifty(ivec_y,ivec_x);
+                        vtable(ivec_y,ivec_x)=-vector(2)+shifty(ivec_y,ivec_x);
                         wtable(ivec_y,ivec_x)=vector(3);
                         xtable(ivec_y,ivec_x)=iref+utable(ivec_y,ivec_x)/2-0.5;% convec flow (velocity taken at the point middle from imgae 1 and 2)
@@ -938 +954 @@
                         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);
-                       % if utable(ivec_y,ivec_x)==0 && vtable(ivec_y,ivec_x)==0
-                       %     'test'
-                       % end
+                        % if utable(ivec_y,ivec_x)==0 && vtable(ivec_y,ivec_x)==0
+                        %     'test'
+                        % end
                         % 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))
@@ -956 +972 @@
     end
 end
+ytable=npy_ima-ytable+1;%reverse from j index to image coordinate y
 result_conv=result_conv*corrmax/(255*sum_square);% keep the last correlation matrix for output
 
@@ -980 +997 @@
     f2 = log(result_conv(z+1,y,x));
     peakz = peakz+ (f1-f2)/(2*f1-4*f0+2*f2);
-
+    
     f1 = log(result_conv(z,y-1,x));
     f2 = log(result_conv(z,y+1,x));
     peaky = peaky+ (f1-f2)/(2*f1-4*f0+2*f2);
-
+    
     f1 = log(result_conv(z,y,x-1));
     f2 = log(result_conv(z,y,x+1));
@@ -1084 +1101 @@
 
 if isfield (Param,'MinCorr')
-     FF(C<Param.MinCorr & FFIn==0)=2;
+    FF(C<Param.MinCorr & FFIn==0)=2;
 end
 if (isfield(Param,'MinVel')&&~isempty(Param.MinVel))||(isfield (Param,'MaxVel')&&~isempty(Param.MaxVel))
@@ -1093 +1110 @@
     end
     if isfield (Param,'MaxVel')&&~isempty(Param.MaxVel)
-         U2Max=Param.MaxVel*Param.MaxVel;
+        U2Max=Param.MaxVel*Param.MaxVel;
         FF(Umod>U2Max & FFIn==0)=3;
     end
@@ -1139 +1156 @@
         i1_series=i_series-ind1;% set of first image numbers
         i2_series=i_series+ind2;
-        check_bounds=i1_series<MinIndex_i | i2_series>MaxIndex_i;
+        check_bounds=i1_series(1)<MinIndex_i | i2_series(1)>MaxIndex_i;
         if isempty(j_series)
             NomTypeNc='_1-2';

trunk/src/series/civ_input.m

@@ -101 +101 @@
 NomTypeImaA=NomTypeInput;
 iview_image=1;%line # for the input images
-switch FileType
-    case 'civdata'
-        if ~strcmp(Param.Action.ActionName,'civ_series')
-            msgbox_uvmat('ERROR','bad input data file: open an image or a nc file from civ_series')
-            return
-        end
+if ismember( FileType,{'civdata','civdata_3D'})
         NomTypeNc=NomTypeInput;
         ind_opening=SeriesData.FileInfo{1}.CivStage;
@@ -120 +115 @@
         end
         
-        if size(Param.InputTable,1)==1
-            
-            Param.InputTable(2,:)=Param.InputTable(1,:);
-              set(hhseries.InputTable,'Data',Param.InputTable)
+        if size(Param.InputTable,1)==1     
              if isfield(Data,'Civ2_ImageA')
                  ImageName=Data.Civ2_ImageA;
              elseif isfield(Data,'Civ1_ImageA')
                  ImageName=Data.Civ1_ImageA;
+             else
+                  msgbox_uvmat('ERROR','no original image defined in netcdf input file ')
+            return
              end
-            series('display_file_name',hhseries,ImageName,1);%append the image series to the input list
+            series('display_file_name',hhseries,ImageName,'append');%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')
-        % msgbox_uvmat('ERROR','civ_series needs images, scalar fields in netcdf format, or civ data as input')
-        % return
-        % end
+        end
+
+        iview_image=2;%line # for the input images
 end
         
@@ -162 +146 @@
      return
 end
+
 MaxIndex_i=Param.IndexRange.MaxIndex_i(1);
 MinIndex_i=Param.IndexRange.MinIndex_i(1);

trunk/src/series/civ_series.m

@@ -117 +117 @@
         if Param.ActionInput.CheckCiv1
             iview_A=1;% usual PIV, the image series is on the first line of the table
-        elseif Param.ActionInput.CheckCiv2 % civ2 is performed without Civ1, a netcdf file series is needed in the first table line
-            iview_A=2;% the second line is used for the input images of Civ2
+        else % Civ1 has been already stored in a netcdf file input 
+            iview_A=2;% the second line is used for the input images 
         end
         if iview_A~=0
@@ -274 +274 @@
     tstart=tic;
     time_civ1=0;
-  time_patch1=0;
-  time_civ2=0;
-  time_patch2=0;
+    time_patch1=0;
+    time_civ2=0;
+    time_patch2=0;
     if ~isempty(RUNHandle)% update the waitbar in interactive mode with GUI series  (checkrun=1)
         update_waitbar(WaitbarHandle,ifield/NbField)
@@ -335 +335 @@
     %% Civ1
     % if Civ1 computation is requested
-    if isfield (Param.ActionInput,'Civ1')
+    if Param.ActionInput.CheckCiv1
         if CheckInputFile
             disp('civ1 started')
@@ -377 +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 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;
-                    end
-                    if isempty(regexp(ImageName_B,'(^http://)|(^https://)', 'once')) && ~exist(ImageName_B,'file')
-                        disp([ImageName_B ' missing'])
-                        continue
-                    end
-                    [par_civ1.ImageB,VideoObject_B] = read_image(ImageName_B,FileType_B,VideoObject_B,FrameIndex_B_Civ1(ifield));
+                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;
+                end
+                if isempty(regexp(ImageName_B,'(^http://)|(^https://)', 'once')) && ~exist(ImageName_B,'file')
+                    disp([ImageName_B ' missing'])
+                    continue
+                end
+                [par_civ1.ImageB,VideoObject_B] = read_image(ImageName_B,FileType_B,VideoObject_B,FrameIndex_B_Civ1(ifield));
             catch ME % display errors in reading input images
                 if ~isempty(ME.message)
@@ -392 +392 @@
                 end
             end
+            
             % par_civ1.ImageWidth=size(par_civ1.ImageA,2);
             % par_civ1.ImageHeight=size(par_civ1.ImageA,1);
@@ -425 +426 @@
             end
             Data.ListGlobalAttribute=[ListGlobalAttribute Civ1_param];
+            
             Data.CivStage=1;
         else
@@ -438 +440 @@
         Data.VarAttribute{5}.Role='ancillary';
         Data.VarAttribute{6}.Role='errorflag';
-
+        
         % case of mask
         if par_civ1.CheckMask&&~isempty(par_civ1.Mask)
@@ -451 +453 @@
                 i1_mask=mod(i1-1,par_civ1.NbSlice)+1;
                 maskname=fullfile_uvmat(RootPath_mask,SubDir_mask,RootFile_mask,Ext_mask,'_1',i1_mask);
-            else
-                maskname=Param.ActionInput.Civ1.Mask;
-            end
-            if strcmp(maskoldname,maskname)% mask exist, not already read in civ1
-                par_civ1.Mask=mask; %use mask already opened
-            else
-                if ~isempty(regexp(maskname,'(^http://)|(^https://)', 'once'))|| exist(maskname,'file')
-                    try
-                        par_civ1.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
+                if strcmp(Param.ActionInput.PairIndices.ListPairMode,'series(Di)')% case of volume, mask index refers to j index
+                    par_civ1.NbSlice_j=par_civ1.NbSlice;
+                end
+            end
+        else
+            maskname=Param.ActionInput.Civ1.Mask;
+        end
+        if strcmp(maskoldname,maskname)% mask exist, not already read in civ1
+            par_civ1.Mask=mask; %use mask already opened
+        else
+            if ~isempty(regexp(maskname,'(^http://)|(^https://)', 'once'))|| exist(maskname,'file')
+                try
+                    par_civ1.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
-                else
-                    par_civ1.Mask=[];
-                end
-                mask=par_civ1.Mask;
-                maskoldname=maskname;
-            end
-        end
-
+                end
+            else
+                par_civ1.Mask=[];
+            end
+            mask=par_civ1.Mask;
+            maskoldname=maskname;
+        end
+        
+        
         % case of input grid
         if par_civ1.CheckGrid &&~isempty(par_civ1.Grid)
@@ -481 +487 @@
             par_civ1.CorrBoxSize=GridData.CorrBox;
         end
-
-        % caluclate velocity data 
+        
+        % 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);
@@ -489 +495 @@
             return
         end
-
-    else% we use existing Civ1 data
+        
         if exist('ncfile','var')
             CivFile=ncfile;
-            [Data,tild,tild,errormsg]=nc2struct(CivFile,'ListGlobalAttribute','absolut_time_T0'); %look for the constant 'absolut_time_T0' to detect old civx data format
-            if ~isempty(errormsg)
-                disp_uvmat('ERROR',errormsg,checkrun)
-                return
-            end
+            % [Data,tild,tild,errormsg]=nc2struct(CivFile,'ListGlobalAttribute','absolut_time_T0'); %look for the constant 'absolut_time_T0' to detect old civx data format
+            % if ~isempty(errormsg)
+            %     disp_uvmat('ERROR',errormsg,checkrun)
+            %     return
+            % end
             [Data,tild,tild,errormsg]=nc2struct(CivFile);%read civ1 and fix1 data in the existing netcdf file
         elseif isfield(Param,'Civ1_X')
@@ -510 +515 @@
         end
     end
+
     
     %% Fix1
@@ -717 +723 @@
                 i1_mask=mod(i1-1,par_civ2.NbSlice)+1;
                 maskname=fullfile_uvmat(RootPath_mask,SubDir_mask,RootFile_mask,Ext_mask,'_1',i1_mask);
+                if strcmp(Param.ActionInput.PairIndices.ListPairMode,'series(Di)')% case of volume, mask index refers to j index
+                    par_civ2.NbSlice_j=par_civ2.NbSlice;
+                end
             else
                 maskname=Param.ActionInput.Civ2.Mask;
@@ -876 +885 @@
         if isempty(errormsg)
             disp([ncfile_out ' written'])
-            %[success,msg] = fileattrib(ncfile_out ,'+w','g');% done in struct2nc
         else
             disp(errormsg)

trunk/src/series/extract_rdvision.m

@@ -160 +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);
@@ -457 +444 @@
 end
 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(iinew-1,nbfield2)+1;

trunk/src/uvmat.m

@@ -2198 +2198 @@
 %------------------------------------------------------------------------
 % --- Fills the edit boxes RootPath, RootFile,NomType...from an input file name 'fileinput'
-function errormsg=display_file_name(handles,fileinput,index)
+function errormsg=display_file_name(handles,fileinput,input_line)
 %------------------------------------------------------------------------
 %% look for the input file existence
@@ -2208 +2208 @@
 end
 
-%% define the relevant handles for the first field series (index=1) or the second file series (index=2)
-if ~exist('index','var')
-    index=1;
-end
-if index==1
+%% define the relevant handles for the first field series (input_line=1) or the second file series (input_line=2)
+if ~exist('input_line','var')
+    input_line=1;
+end
+if input_line==1
     handles_RootPath=handles.RootPath;
     handles_SubDir=handles.SubDir;
@@ -2219 +2219 @@
     handles_NomType=handles.NomType;
     handles_FileExt=handles.FileExt;
-elseif index==2
+elseif input_line==2
     handles_RootPath=handles.RootPath_1;
     handles_SubDir=handles.SubDir_1;
@@ -2287 +2287 @@
         set(handles_NomType,'String',NomType);
         set(handles_FileExt,'String',FileExt);
-        if index==1
+        if input_line==1
             % fill file index counters if the first file series is opened
             set(handles.i1,'String',num2str(i1));
@@ -2293 +2293 @@
             set(handles.j1,'String',num2stra(j1,NomType));
             set(handles.j2,'String',num2stra(j2,NomType));
+            if isfield(FileInfo,'MaskFile')
+                set(handles.CheckMask,'Value',1)
+                 Mask.File=FileInfo.MaskFile;
+                 if isfield(FileInfo,'MaskNbSlice')
+                     Mask.NbSlice=FileInfo.MaskNbSlice;
+                 elseif isfield(FileInfo,'VolumeScan')
+                     Mask.VolumeScan=FileInfo.VolumeScan;
+                 end  
+                 set(handles.CheckMask,'UserData', Mask)
+            else
+                set(handles.CheckMask,'Value',0)
+                CheckMask_Callback(handles.CheckMask, [], handles) 
+            end                         
         else %read the current field index to synchronise with the first series
             i1_s=str2num(get(handles.i1,'String'));
@@ -2351 +2364 @@
         
         % initiate input file series and inputfilerefresh the current field view:
-        update_rootinfo(handles,i1_series,i2_series,j1_series,j2_series,FileInfo,MovieObject,index);
+        update_rootinfo(handles,i1_series,i2_series,j1_series,j2_series,FileInfo,MovieObject,input_line);
 end
 
@@ -2379 +2392 @@
 % --- Update information about a new field series (indices to scan, timing,
 %     calibration from an xml file, then inputfilerefresh current plots
-function errormsg=update_rootinfo(handles,i1_series,i2_series,j1_series,j2_series,FileInfo,VideoObject,index)
+function errormsg=update_rootinfo(handles,i1_series,i2_series,j1_series,j2_series,FileInfo,VideoObject,input_line)
 %------------------------------------------------------------------------
 errormsg=''; %default error msg
 
-%% define the relevant handles depending on the index (1=first file series, 2= second file series)
-if ~exist('index','var')
-    index=1;
-end
-if index==1
+%% define the relevant handles depending on the input_line (1=first file series, 2= second file series)
+if ~exist('input_line','var')
+    input_line=1;
+end
+if input_line==1
     handles_Fields=handles.FieldName;
-elseif index==2
+elseif input_line==2
     handles_Fields=handles.FieldName_1;
 end
@@ -2398 +2411 @@
 UvData.NewSeries=1; %flag for REFRESH: begin a new series
 UvData.FileName_1='';% name of the current second field (used to detect a  constant field during file scanning)
-%UvData.FileType{index}=FileInfo.FileType;
-UvData.FileInfo{index}=FileInfo;
-UvData.MovieObject{index}=VideoObject;
-UvData.i1_series{index}=i1_series;
-UvData.i2_series{index}=i2_series;
-UvData.j1_series{index}=j1_series;
-UvData.j2_series{index}=j2_series;
+%UvData.FileType{input_line}=FileInfo.FileType;
+UvData.FileInfo{input_line}=FileInfo;
+UvData.MovieObject{input_line}=VideoObject;
+UvData.i1_series{input_line}=i1_series;
+UvData.i2_series{input_line}=i2_series;
+UvData.j1_series{input_line}=j1_series;
+UvData.j2_series{input_line}=j2_series;
 nbfield=max(max(max(i2_series)));% total number of fields (i index)
 if isempty(nbfield)
@@ -2447 +2460 @@
 %% read parameters (time, geometric calibration..) from a documentation file (.xml advised)
 XmlData.GeometryCalib=[];%default
-if index==1
+if input_line==1
     [RootPath,SubDir,RootFile,FileIndices,FileExt]=read_file_boxes(handles);
 else
@@ -2515 +2528 @@
     TimeName='civdata_3D';
 end
-if index==1
+if input_line==1
     set(handles.TimeName,'String',TimeName)
 else
@@ -2531 +2544 @@
 last_i_cell=get(handles.MaxIndex_i,'String');
 if isempty(nbfield)
-    last_i_cell{index}='';
-else
-    last_i_cell{index}=num2str(nbfield);
+    last_i_cell{input_line}='';
+else
+    last_i_cell{input_line}=num2str(nbfield);
 end
 set(handles.MaxIndex_i,'String',last_i_cell)
 last_j_cell=get(handles.MaxIndex_j,'String');
 if isempty(nbfield_j)
-     last_j_cell{index}='';
-else
-     last_j_cell{index}=num2str(nbfield_j);
+     last_j_cell{input_line}='';
+else
+     last_j_cell{input_line}=num2str(nbfield_j);
 end
 set(handles.MaxIndex_j,'String',last_j_cell);
@@ -2552 +2565 @@
         set(handles.pxcmx,'Visible','off')
         set(handles.pxcmy,'Visible','off')
-        if index==1
+        if input_line==1
         set(handles.TransformName,'Value',1); %  no transform by default
         end
@@ -2568 +2581 @@
             set(handles.pxcmy,'String',num2str(pixcmy))
         end
-        if ~get(handles.CheckFixLimits,'Value')&& index==1
+        if ~get(handles.CheckFixLimits,'Value')&& input_line==1
             set(handles.TransformName,'Value',3); % phys transform by default if fixedLimits is off
         end
@@ -2591 +2604 @@
 %% update the data attached to the uvmat interface
 if ~isempty(TimeUnit)
-    if index==2 && isfield(UvData,'TimeUnit') && ~strcmp(UvData.TimeUnit,TimeUnit)
+    if input_line==2 && isfield(UvData,'TimeUnit') && ~strcmp(UvData.TimeUnit,TimeUnit)
         msgbox_uvmat('WARNING',['time unit for second file series ' TimeUnit ' inconsistent with first series'])
     else
@@ -2597 +2610 @@
     end
 end
-UvData.XmlData{index}=XmlData;
+UvData.XmlData{input_line}=XmlData;
 UvData.NewSeries=1;
 set(handles.uvmat,'UserData',UvData)
@@ -2612 +2625 @@
         set(handles_Fields,'String',[{'image'};FieldList;{'get_field...'}]);%standard menu for civx data
         set(handles_Fields,'Value',2) % set menu to 'velocity
-        if index==1
+        if input_line==1
             set(handles.FieldName_1,'Value',1);
             set(handles.FieldName_1,'String',[{''};{'image'};FieldList;{'get_field...'}]);%standard menu for civx data reproduced for the second field
@@ -2627 +2640 @@
         set(handles_Fields,'String',[{'image'};FieldList;{'get_field...'}]);%standard menu for civx data
         set(handles_Fields,'Value',2) % set menu to 'velocity
-        if index==1
+        if input_line==1
             set(handles.FieldName_1,'Value',1);
             set(handles.FieldName_1,'String',[{''};{'image'};FieldList;{'get_field...'}]);%standard menu for civx data reproduced for the second field
@@ -2641 +2654 @@
         set(handles_Fields,'Value',1)
         set(handles_Fields,'String',{'get_field...'})
-        if index==1
+        if input_line==1
             FieldName_Callback([],[], handles)
         else
@@ -2660 +2673 @@
 scan_option='i';%default
 state_j='off'; %default
-if index==2
+if input_line==2
     if get(handles.scan_j,'Value')
         scan_option='j'; %keep the scan option for the second file series
@@ -2672 +2685 @@
     if ~isempty(j1_series)
         state_j='on';
-        if index==1
+        if input_line==1
             if isequal(ref_i,ref_i(1)*ones(size(ref_j)))% if ref_i is always equal to its first value
                 scan_option='j'; %scan j indext
@@ -2707 +2720 @@
 if ~isempty(i2_series)||~isempty(j2_series)
     set(handles.CheckFixPair,'Visible','on')
-elseif index==1
+elseif input_line==1
     set(handles.CheckFixPair,'Visible','off')
 end
@@ -2713 +2726 @@
 %% apply the effect of the transform fct and view the field
 transform=get(handles.TransformPath,'UserData');
-if index==2 && (~isa(transform,'function_handle')||nargin(transform)<3)
+if input_line==2 && (~isa(transform,'function_handle')||nargin(transform)<3)
     set(handles.TransformName,'value',2); % set transform to sub_field if the current fct doe not accept two input fields
 end
@@ -2908 +2921 @@
     end
     if mdetect==0 % if no mask is detected in the current folder
-        MaskFullName=uigetfile_uvmat('pick a mask image file:',RootPath,'image');
-        if isempty(MaskFullName)
-            set(handles.CheckMask,'Value',0)
-        end
-
-        [MaskPath,MaskName,MaskExt]=fileparts(MaskFullName);
-        [tild,tild,MaskFile,i1_series,i2,j1,j2,MaskNomType]=find_file_series(MaskPath,[MaskName MaskExt],0);
+         filemask= uigetfile_uvmat('pick a mask image file:',RootPath,'image');
+    if ~isempty(filemask)
+        [FilePath,FileName,FileExt]=fileparts(filemask);
+        [RootPath,SubDir,RootFile,i1_series,i2,j1,j2,NomType]=find_file_series(FilePath,[FileName FileExt]);
         if strcmp(NomType,'_1')
             NbSlice=i1_series(1,2,end);
-           set(handles.num_NbSlice,'String',num2str(NbSlice))
-        end
+            set(handles.num_NbSlice,'String',num2str(NbSlice))
+        elseif ~strcmp(NomType,'*')
+            msgbox_uvmat('ERROR','multilevel masks must be labeled with a single index as _1,_2,...');
+            return
+        end
+        set(hObject,'UserData',filemask);%store for future use
+        testmask=1;
+    end
+        
+        %TO COMPLEMENT......................
+         
+        
+%         MaskFullName=uigetfile_uvmat('pick a mask image file:',RootPath,'image');
+%         if isempty(MaskFullName)
+%             set(handles.CheckMask,'Value',0)
+%         end
+% 
+%         [MaskPath,MaskName,MaskExt]=fileparts(MaskFullName);
+%         [tild,tild,MaskFile,i1_series,i2,j1,j2,MaskNomType]=find_file_series(MaskPath,[MaskName MaskExt],0);
+%         if strcmp(MaskNomType,'_1')
+%             NbSlice=i1_series(1,2,end);
+%            set(handles.num_NbSlice,'String',num2str(NbSlice))
+%         end
     end
     % Mask.Path=MaskPath;
@@ -2961 +2992 @@
 %------------------------------------------------------------------------
 errormsg=[];%default
-Mask=get(handles.CheckMask,'UserData');
-if strcmp(get(handles.z_index,'Visible'),'on')
-    MaskIndex_i=str2num(get(handles.z_index,'String'));
-else
-    MaskIndex_i=mod(str2num(get(handles.i1,'String'))-1,Mask.NbSlice_i)+1;
-end
-MaskIndex_z=MaskIndex_i;%default
-if Mask.NbSlice_j>1
-    MaskIndex_j=str2num(get(handles.j1,'String'));
-    MaskIndex_z=MaskIndex_j;
-else
-    MaskIndex_j=1;
-end
-if isfield(Mask,'maskhandle')&& ishandle(Mask.maskhandle)
-    uistack(Mask.maskhandle,'top');
-end
-MaskName=fullfile_uvmat(Mask.Path,'',Mask.File,Mask.Ext,Mask.NomType,MaskIndex_i,[],MaskIndex_j);
-UvData=get(handles.uvmat,'UserData');
-
-%% update mask image if the mask is new
-if ~ (isfield(UvData,'MaskName') && isequal(UvData.MaskName,MaskName))
-    UvData.MaskName=MaskName; %update the recorded name on UvData
-    set(handles.uvmat,'UserData',UvData);
-    if ~exist(MaskName,'file')
-        if isfield(Mask,'maskhandle')&& ishandle(Mask.maskhandle)
-            delete(Mask.maskhandle)
-        end
+MaskName='';
+
+%% get the current mask name recorded in CheckMask/UserData, possibly indexed with file index
+MaskInfo=get(handles.CheckMask,'UserData');
+if isfield(MaskInfo,'File')
+    if isfield(MaskInfo,'NbSlice')
+        if isfield(MaskInfo,'VolumeScan') &&  MaskInfo.VolumeScan
+            MaskIndex_i=str2num(get(handles.j1,'String'));
+        else
+            MaskIndex_i=mod(str2num(get(handles.i1,'String'))-1,MaskInfo.NbSlice)+1;
+        end
+        MaskName=[MaskInfo.File '_' num2str(MaskIndex_i) '.png'];
     else
-        %read mask image
-        [MaskField,tild,errormsg] = read_field(MaskName,'image');
-        if ~isempty(errormsg)
-            return
-        end
-        npxy=size(MaskField.A);
-        if length(npxy)>2
-            errormsg=[MaskName ' is not a grey scale image'];
-            return
-        elseif ~isa(MaskField.A,'uint8')
-            errormsg=[MaskName ' is not a 8 bit grey level image'];
-            return
-        end
-        MaskField.ZIndex=MaskIndex_z;
-        %px to phys or other transform on field
-         menu_transform=get(handles.TransformName,'String');
-        choice_value=get(handles.TransformName,'Value');
-        transform_name=menu_transform{choice_value};%name of the transform fct  given by the menu 'transform_fct'
-        transform=get(handles.TransformPath,'UserData');
-        if  ~isequal(transform_name,'') && ~isequal(transform_name,'px')
-            if isfield(UvData,'XmlData') && isfield(UvData.XmlData{1},'GeometryCalib')%use geometry calib recorded from the ImaDoc xml file as first priority
-                Calib=UvData.XmlData{1}.GeometryCalib;
-                MaskField=transform(MaskField,UvData.XmlData{1});
-            end
-        end
-        flagmask=MaskField.A < 200;
-
-        %make brown color image
-        imflag(:,:,1)=0.9*flagmask;
-        imflag(:,:,2)=0.7*flagmask;
-        imflag(:,:,3)=zeros(size(flagmask));
-
-        %update mask image
-        hmask=[]; %default
-        if isfield(Mask,'maskhandle')&& ishandle(Mask.maskhandle)
-            hmask=Mask.maskhandle;
-        end
-        if ~isempty(hmask)
-            set(hmask,'CData',imflag)
-            set(hmask,'AlphaData',flagmask*0.6)
-            set(hmask,'XData',MaskField.Coord_x);
-            set(hmask,'YData',MaskField.Coord_y);
-%             uistack(hmask,'top')
+        MaskIndex_i=1;
+        MaskName=MaskInfo.MaskFile;
+    end
+    
+    %% update mask image if the mask is new
+    UvData=get(handles.uvmat,'UserData');
+    if ~ (isfield(UvData,'MaskName') && isequal(UvData.MaskName,MaskName))
+        UvData.MaskName=MaskName; %update the recorded name on UvData
+        set(handles.uvmat,'UserData',UvData);
+        if ~exist(MaskName,'file')
+            if isfield(Mask,'maskhandle')&& ishandle(Mask.maskhandle)
+                delete(Mask.maskhandle)
+            end
         else
-            axes(handles.PlotAxes)
-            hold on
-            Mask.maskhandle=image(MaskField.Coord_x,MaskField.Coord_y,imflag,'Tag','mask','HitTest','off','AlphaData',0.6*ones(size(flagmask)));
-            set(handles.CheckMask,'UserData',Mask)
-        end
-    end
-end
-
+            %read mask image
+            [MaskField,tild,errormsg] = read_field(MaskName,'image');
+            if ~isempty(errormsg)
+                return
+            end
+            npxy=size(MaskField.A);
+            if length(npxy)>2
+                errormsg=[MaskName ' is not a grey scale image'];
+                return
+            elseif ~isa(MaskField.A,'uint8')
+                errormsg=[MaskName ' is not a 8 bit grey level image'];
+                return
+            end
+            MaskField.ZIndex=MaskIndex_i;
+            %px to phys or other transform on field
+            menu_transform=get(handles.TransformName,'String');
+            choice_value=get(handles.TransformName,'Value');
+            transform_name=menu_transform{choice_value};%name of the transform fct  given by the menu 'transform_fct'
+            transform=get(handles.TransformPath,'UserData');
+            if  ~isequal(transform_name,'') && ~isequal(transform_name,'px')
+                if isfield(UvData,'XmlData') && isfield(UvData.XmlData{1},'GeometryCalib')%use geometry calib recorded from the ImaDoc xml file as first priority
+                    Calib=UvData.XmlData{1}.GeometryCalib;
+                    MaskField=transform(MaskField,UvData.XmlData{1});
+                end
+            end
+            flagmask=MaskField.A < 200;
+            
+            %make brown color image
+            imflag(:,:,1)=0.9*flagmask;
+            imflag(:,:,2)=0.7*flagmask;
+            imflag(:,:,3)=zeros(size(flagmask));
+            
+            %update mask image
+            hmask=[]; %default
+            if isfield(MaskInfo,'maskhandle')&& ishandle(MaskInfo.maskhandle)
+                hmask=MaskInfo.maskhandle;
+            end
+            if ~isempty(hmask)
+                set(hmask,'CData',imflag)
+                set(hmask,'AlphaData',flagmask*0.6)
+                set(hmask,'XData',MaskField.Coord_x);
+                set(hmask,'YData',MaskField.Coord_y);
+                %             uistack(hmask,'top')
+            else
+                axes(handles.PlotAxes)
+                hold on
+                MaskInfo.maskhandle=image(MaskField.Coord_x,MaskField.Coord_y,imflag,'Tag','mask','HitTest','off','AlphaData',0.6*ones(size(flagmask)));
+                set(handles.CheckMask,'UserData',MaskInfo)
+            end
+        end
+    end
+end
 %------------------------------------------------------------------------
 %------------------------------------------------------------------------