Changeset 1161

Timestamp:

Jul 19, 2024, 5:57:46 AM (5 months ago)

Author:

sommeria

Message:

filter_tps modified to avoid possible bugs in case of few vectors, subpixel interpolation slightly modified in civ

Location:

trunk/src

Files:

• filter_tps.m (modified) (7 diffs)
• filter_tps_3D.m (added)
• get_file_info.m (modified) (1 diff)
• proj_field.m (modified) (1 diff)
• read_civdata.m (modified) (2 diffs)
• script_delete_PCO.m (modified) (3 diffs)
• script_delete_rdvision.m (modified) (1 diff)
• series.m (modified) (3 diffs)
• series/civ_3D.m (modified) (6 diffs)
• series/civ_input.m (modified) (6 diffs)
• series/civ_series.m (modified) (4 diffs)
• series/extract_rdvision.m (modified) (1 diff)
• series/stereo_civ.m (modified) (23 diffs)
• series/stereo_input.fig (modified) (previous)
• series/stereo_input.m (modified) (11 diffs)
• series/test_filter_tps.m (modified) (1 diff)
• tps_coeff.m (modified) (1 diff)
• uvmat.m (modified) (8 diffs)

trunk/src/filter_tps.m

@@ -15 +15 @@
 % coord=[X Y]: matrix whose first column is the x coordinates of the initial data, the second column the y coordiantes
 % U,V, possibly W: set of velocity components of the initial data
-% SubdomainSize: estimated number of data points in each subdomain
+% SubdomainSize: estimated number of data points in each subdomain, choose a number > 125 (must be >= 12 for convergence of tps after division by 4)
 % FieldSmooth: smoothing parameter
-% Threshold: max diff accepted between smoothed and initial data 
+% Threshold: max diff accepted between smoothed and initial data
 
 
@@ -49 +49 @@
 AspectRatio=Range(2)/Range(1);
 NbSubDomain=NbVec/SubDomainSize;% estimated number of subdomains
-NbSubDomainX=max(floor(sqrt(NbSubDomain/AspectRatio)),1);% estimated number of subdomains in x 
+NbSubDomainX=max(floor(sqrt(NbSubDomain/AspectRatio)),1);% estimated number of subdomains in x
 NbSubDomainY=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
@@ -78 +78 @@
 FF=false(NbVec,1);%false flag=0 (false) by default
 nb_select=zeros(NbVec,1);
-check_empty=zeros(1,NbSubDomain);
+check_empty=false(1,NbSubDomain);
 
 %% calculate tps coeff in each subdomain
@@ -84 +84 @@
     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
-    ind_sel_previous=[];
     ind_sel=0;%initialize set of vector indices in the subdomain
-    %increase iteratively the subdomain if it contains less than SubDomainNbVec/4 source vectors
-    while numel(ind_sel)>numel(ind_sel_previous)
-        ind_sel_previous=ind_sel;% record the set of selected vector indices for next iteration 
+    %increase iteratively the subdomain if it contains less than
+    %SubDomainNbVec/4 source vectors, until possibly cover the whole domain:check_partial_domain=false
+    check_partial_domain= true;
+    while check_partial_domain
+        %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( 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
         if isempty(ind_sel)
-            check_empty(isub)=1;
-            break 
-        % if too few selected vectors, increase the subrange for next iteration
-        elseif numel(ind_sel)<SubDomainSize/4 && ~isequal( ind_sel,ind_sel_previous)
-            SubRange(:,1,isub)=SubRange(:,1,isub)-Siz'/4;
-            SubRange(:,2,isub)=SubRange(:,2,isub)+Siz'/4;
-        % if subdomain includes enough vectors, perform tps interpolation
+            check_empty(isub)=true;
+            break
+
+            % if too few selected vectors, increase the subrange for next iteration by 50 % in each direction
+        elseif numel(ind_sel)<SubDomainSize/4 && ~isequal( ind_sel,ind_sel_previous)&& check_partial_domain
+            % SubRange(:,1,isub)=SubRange(:,1,isub)-Siz'/4;
+            % SubRange(:,2,isub)=SubRange(:,2,isub)+Siz'/4;
+            SubRange(:,1,isub)=1.25*SubRange(:,1,isub)-0.25*SubRange(:,2,isub);
+            SubRange(:,2,isub)=-0.25*SubRange(:,1,isub)+1.25*SubRange(:,2,isub);
+
+            % if subdomain includes enough vectors, perform tps interpolation
         else
             [U_smooth_sub,U_tps_sub]=tps_coeff(Coord(ind_sel,:),U(ind_sel),smoothing);
@@ -109 +118 @@
             if exist('Threshold','var')&&~isempty(Threshold)
                 FF(ind_sel)=(NormDiff>Threshold);%put FF value to 1 to identify the criterium of elimmination
-                ind_ind_sel=find(~FF(ind_sel)); % select the indices of remaining vectors in the subset of ind_sel vectors 
+                ind_ind_sel=find(~FF(ind_sel)); % select the indices of remaining vectors in the subset of ind_sel vectors
             end
             % if no value exceeds threshold, the result is recorded
+
             if isequal(numel(ind_ind_sel),numel(ind_sel))
                 x_width=(SubRange(1,2,isub)-SubRange(1,1,isub))/pi;
@@ -128 +138 @@
                 display(['tps done with ' num2str(numel(ind_sel)) ' vectors in subdomain # ' num2str(isub)  ' among ' num2str(NbSubDomain)])
                 break
-            % if too few selected vectors, increase the subrange for next iteration
-            elseif numel(ind_ind_sel)<SubDomainSize/4 && ~isequal( ind_sel,ind_sel_previous)
-                SubRange(:,1,isub)=SubRange(:,1,isub)-Siz'/4;
-                SubRange(:,2,isub)=SubRange(:,2,isub)+Siz'/4;
+            % if too few selected vectors, increase the subrange for next iteration by 50% in each direction
+            elseif numel(ind_ind_sel)<SubDomainSize/4 && ~isequal( ind_sel,ind_sel_previous)&& check_partial_domain
+                % SubRange(:,1,isub)=SubRange(:,1,isub)-Siz'/4;
+                % SubRange(:,2,isub)=SubRange(:,2,isub)+Siz'/4;
+                SubRange(:,1,isub)=1.25*SubRange(:,1,isub)-0.25*SubRange(:,2,isub);
+                SubRange(:,2,isub)=-0.25*SubRange(:,1,isub)+1.25*SubRange(:,2,isub);
             % else interpolation-smoothing is done again with the selected vectors
             else
@@ -153 +165 @@
             end
         end
+         % check_next=true;% go to next iteration with wider subdomain
+    end% end of while loop fo increasing subdomain size
+    if ~check_partial_domain && isub<NbSubDomain% if the while loop proceed to the end, the whole domain size has been covered
+        check_empty(isub+1:end)=true;
+        break %the whole domain has been already covered, no need for new subdomains
     end
 end

trunk/src/get_file_info.m

@@ -170 +170 @@
                                 end
                             elseif isfield(Data,'Conventions') && strcmp(Data.Conventions,'uvmat/civdata')
-                                FileInfo.FileType='civdata'; % test for civx velocity fields
+                                FileInfo.FileType='civdata'; % test for civ velocity fields
+                                FileInfo.CivStage=Data.CivStage;
+                            elseif isfield(Data,'Conventions') && strcmp(Data.Conventions,'uvmat/civdata_3D')
+                                FileInfo.FileType='civdata_3D'; % test for 3D volume civ velocity fields
                                 FileInfo.CivStage=Data.CivStage;
                             else

trunk/src/proj_field.m

@@ -1723 +1723 @@
                         ProjData.coord_y=coord_y_proj;
                         ProjData.(VarName)=interp3(Coord{3},Coord{2},Coord{1},double(FieldData.(VarName)),XI_proj,YI_proj,ZI_proj,'*linear');
-                        ProjData.(VarName)=nanmean(ProjData.(VarName),3);
+                        ProjData.(VarName)=mean(ProjData.(VarName),3,'omitnan');
                         ProjData.(VarName)=squeeze(ProjData.(VarName));
                     end

trunk/src/read_civdata.m

@@ -218 +218 @@
 var={};
 switch vel_type
-    case{'civ1','civ2','civ3'}
-        varout={'X','Y','Z','U','V','W','C','F','FF'};
-        varin= {'Civ1_X','Civ1_Y','Civ1_Z','Civ1_U','Civ1_V','Civ1_W','Civ1_C','Civ1_F','Civ1_FF'};
-         role={'coord_x','coord_y','coord_z','vector_x','vector_y','vector_z','ancillary','warnflag','errorflag'};  
-    case{'filter1','filter2','filter3'}  
-        varout={'X','Y','Z','U','V','W','C','F','FF','Coord_tps','U_tps','V_tps','W_tps','SubRange','NbCentre','NbCentre','NbCentre'};
-        varin={'Civ1_X','Civ1_Y','Civ1_Z','Civ1_U_smooth','Civ1_V_smooth','Civ1_W','Civ1_C','Civ1_F','Civ1_FF',...
+    case{'civ1','civ2'}
+        varout={'X','Y','Z','U','V','W','C','FF'};
+        varin= {'Civ1_X','Civ1_Y','Civ1_Z','Civ1_U','Civ1_V','Civ1_W','Civ1_C','Civ1_FF'};
+         role={'coord_x','coord_y','coord_z','vector_x','vector_y','vector_z','ancillary','errorflag'};  
+    case{'filter1','filter2'}  
+        varout={'X','Y','Z','U','V','W','C','FF','Coord_tps','U_tps','V_tps','W_tps','SubRange','NbCentre','NbCentre','NbCentre'};
+        varin={'Civ1_X','Civ1_Y','Civ1_Z','Civ1_U_smooth','Civ1_V_smooth','Civ1_W','Civ1_C','Civ1_FF',...
             'Civ1_Coord_tps','Civ1_U_tps','Civ1_V_tps','Civ1_W_tps','Civ1_SubRange','Civ1_NbCentre','Civ1_NbCentres','Civ1_NbSites'};
-        role={'coord_x','coord_y','coord_z','vector_x','vector_y','vector_z','ancillary','warnflag','errorflag','coord_tps','vector_x_tps',...
+        role={'coord_x','coord_y','coord_z','vector_x','vector_y','vector_z','ancillary','errorflag','coord_tps','vector_x_tps',...
             'vector_y_tps','vector_z_tps','ancillary','ancillary','ancillary','ancillary'};
           %  rmq: NbCentres and NbSites obsolete replaced by NbCentre, kept for consistency with previous data
@@ -233 +233 @@
     case {'civ2','filter2'}
         varin=regexprep(varin,'1','2');
-    case {'civ3','filter3'} 
-        varin=regexprep(varin,'1','3');
+    % case {'civ3','filter3'} 
+    %     varin=regexprep(varin,'1','3');
 end
 var=[varout;varin];
 if ~strcmp(ProjModeRequest,'interp_tps')
-    var=var(:,1:9);%suppress tps if not needed
+    var=var(:,1:8);%suppress tps if not needed
 end
 vel_type_out=vel_type;

trunk/src/script_delete_PCO.m

@@ -54 +54 @@
                 DirPCORaw=dir(TifFolder);
                 DirPCORawCells=struct2cell(DirPCORaw);
-                CheckTif=regexp(DirPCORAwCells(1,:),'.tif$');
+                CheckTif=regexp(DirPCORawCells(1,:),'.tif$');
                 IndexTif=find(~cellfun('isempty',CheckTif));
-                ListTif=DirPCORAwCells(1,IndexTif);
+                ListTif=DirPCORawCells(1,IndexTif);
                 try
-                Info=imfinfo(fullfile(TifFolder,ListTif{1}));
-                RecordLength=numel(Info);%Number of frames in fullfile(TifFolder,FileName)
-                Info=imfinfo(fullfile(TifFolder,ListTif{end}));
-                NumberOfFrames=numel(Info)+(numel(ListTif)-1)*RecordLength%Number of frames in fullfile(TifFolder,FileName)
+                    Info=imfinfo(fullfile(TifFolder,ListTif{1}));
+                    RecordLength=numel(Info);%Number of frames in fullfile(TifFolder,FileName)
+                    Info=imfinfo(fullfile(TifFolder,ListTif{end}));
+                    NumberOfFrames=numel(Info)+(numel(ListTif)-1)*RecordLength%Number of frames in fullfile(TifFolder,FileName)
                 catch ME
                     disp(ME.message)
@@ -67 +67 @@
                 end
                 ind_png=find(strcmp([ListNamesSubSub{isubsub} '.png'],ListNamesSubSub)); %index of the .png folder corresponding to the raw PCO folder
-                PngFolder=fullfile(RootFolder,ListNames{ilist},ListNamesSub{isub},ListNamesSubSub{ind_png})
-                if exist(PngFolder,'dir')
-                    DirPng=dir(PngFolder);
-                    if numel(DirPng)==NumberOfFrames+6
-                        PngCells=struct2cell(DirPng(7:end));
-                        mm=cell2mat(PngCells(4,:));% check the sizes of extracted images
-                        sizemax=max(mm)
-                        sizemin=min(mm(3:end))
-                        disp(['max size(Mbytes)=' num2str(sizemax/1000000)]);
-                        if min(mm(3:end))<0.9*sizemax
-                            status=['WARNING' 'min size(Mbytes)=' num2str(sizemin/1000000)];
+                if isempty(ind_png)
+                    List2extract=[List2extract;TifFolder];
+                else
+                    PngFolder=fullfile(RootFolder,ListNames{ilist},ListNamesSub{isub},ListNamesSubSub{ind_png})
+                    if exist(PngFolder,'dir')
+                        DirPng=dir(PngFolder);
+                        if numel(DirPng)==NumberOfFrames+6
+                            PngCells=struct2cell(DirPng(7:end));
+                            mm=cell2mat(PngCells(4,:));% check the sizes of extracted images
+                            sizemax=max(mm)
+                            sizemin=min(mm(3:end))
+                            disp(['max size(Mbytes)=' num2str(sizemax/1000000)]);
+                            if min(mm(3:end))<0.9*sizemax
+                                status=['WARNING' 'min size(Mbytes)=' num2str(sizemin/1000000)];
+                                List2check=[List2check;PngFolder];
+                            else
+                                List2delete=[List2delete;TifFolder];% approve deletion of the source multitif files
+                            end
+                        else
+                            status=['extraction not finished,' num2str(numel(DirPng)-6) ' images extracted'];
                             List2check=[List2check;PngFolder];
-                        else
-                            List2delete=[List2delete;PngFolder];% approve deletion of the source multitif files
                         end
-                    else
-                        status=['extraction not finished,' num2str(numel(DirPng)-6) ' images extracted'];
-                        List2check=[List2check;PngFolder];
                     end
-                else
-                    List2extract=[List2extract;TifFolder];
                 end
             end
@@ -93 +95 @@
     end
 end
-        
+
+
 
 List2extract

trunk/src/script_delete_rdvision.m

@@ -97 +97 @@
                     catch ME
                         disp(['error in ' filename_seq])
+                        disp(ME.message)
                     end
                     DirPng=dir(PngFolder);         

trunk/src/series.m

@@ -2615 +2615 @@
             set(handles.RUN, 'Enable','On'), set(handles.RUN,'BackgroundColor',[1 0 0]),return,end
 end
+num_first_i_Callback(hObject, eventdata, handles)
+num_first_j_Callback(hObject, eventdata, handles)
 
 %% enable or desable j index visibility
@@ -2620 +2622 @@
     enable_j(handles,'off')
 else% show j index if relevant in the input series
+    if isfield(SeriesData,'j1_series')
     j1_series=SeriesData.j1_series;
     for iview=1:size(j1_series,1)
@@ -2626 +2629 @@
             break
         end
+    end
     end
 end

trunk/src/series/civ_3D.m

@@ -314 +314 @@
             return
         end
-        %for z_index=2:floor((NbSlice-1)/par_civ1.Dz)% loop on slices
+     % 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);%shift the indices in the image block upward by par_civ1.Dz
-            par_civ1.ImageB=circshift(par_civ1.ImageA,-par_civ1.Dz);
+            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
-                j_image_index=j1_series_Civ1(z_index*par_civ1.Dz+SearchRange_z-par_civ1.Dz+iz,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);
@@ -330 +330 @@
             [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)
@@ -523 +522 @@
             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);
         nbval=zeros(size(par_civ2.Grid,1),1);% nbre of interpolated values at each grid point (from the different patch subdomains)
-        if par_civ2.CheckDeformation
-            DUDX=zeros(size(par_civ2.Grid,1),1);
-            DUDY=zeros(size(par_civ2.Grid,1),1);
-            DVDX=zeros(size(par_civ2.Grid,1),1);
-            DVDY=zeros(size(par_civ2.Grid,1),1);
-        end
+
         NbSubDomain=size(SubRange,3);
         for isub=1:NbSubDomain% for each sub-domain of Patch1
@@ -555 +539 @@
                 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);
-                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
@@ -924 +901 @@
                         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>=0.999*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);
@@ -990 +967 @@
 FF=false;
 peakz=z;peaky=y;peakx=x;
-% if z < npz && z > 1 && y < npy && y > 1 && x < npx && x > 1
-%     f0 = log(result_conv(z,y,x));
-%     f1 = log(result_conv(z-1,y,x));
-%     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));
-%     peakx = peakx+ (f1-f2)/(2*f1-4*f0+2*f2);
-% else
-%     FF=true;
-% end
+if z < npz && z > 1 && y < npy && y > 1 && x < npx && x > 1
+    f0 = log(result_conv(z,y,x));
+    f1 = log(result_conv(z-1,y,x));
+    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));
+    peakx = peakx+ (f1-f2)/(2*f1-4*f0+2*f2);
+else
+    FF=true;
+end
 
 vector=[peakx-floor(npx/2)-1 peaky-floor(npy/2)-1 peakz-floor(npz/2)-1];

trunk/src/series/civ_input.m

@@ -132 +132 @@
                     [~,~,~,~,~,~,~,~,NomTypeImaA]=fileparts_uvmat(Data.Civ1_ImageA);
        % [~,~,~,~,~,~,~,~,NomTypeImaB]=fileparts_uvmat(Data.Civ1_ImageB);
-             else
+             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);
@@ -415 +415 @@
 %% Civ1 parameters
 %Param.CheckCiv1=1;
-Param.Civ1.CorrBoxSize=[25 25 1];
-Param.Civ1.SearchBoxSize=[55 55 5];
+Param.Civ1.CorrBoxSize=[31 31 1];
+Param.Civ1.SearchBoxSize=[61 61 5];
 Param.Civ1.SearchBoxShift=[0 0];
 Param.Civ1.CorrSmooth=1;
@@ -432 +432 @@
 %% Patch1 parameters
 %Param.CheckPatch1=1;
-Param.Patch1.FieldSmooth=20;
-Param.Patch1.MaxDiff=2;
+Param.Patch1.FieldSmooth=200;
+Param.Patch1.MaxDiff=1.5;
 Param.Patch1.SubDomainSize=125;
 
@@ -452 +452 @@
 
 %% Patch2 parameters
-Param.Patch2.FieldSmooth=5;
-Param.Patch2.MaxDiff=1.5000;
+Param.Patch2.FieldSmooth=20;
+Param.Patch2.MaxDiff=1;
 Param.Patch2.SubDomainSize=250;
-Param.Patch2.TestPatch2=0;
 
 fill_GUI(Param,handles.civ_input)% fill the elements of the GUI series with the input parameters
@@ -706 +705 @@
     end  
     mode_value=get(handles.ListPairMode,'Value');
-    if isempty(mode_value)
+    if isempty(mode_value)|| mode_value>numel(mode_list)
         mode_value=1;
     end
@@ -835 +834 @@
     mode_list=get(handles.ListPairMode,'String');
     mode_value=get(handles.ListPairMode,'Value');
-    if isempty(mode_value)
+    if isempty(mode_value)||mode_value>numel(mode_list)
         mode_value=1;
     end

trunk/src/series/civ_series.m

@@ -1052 +1052 @@
     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
@@ -1130 +1130 @@
                 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                 corrmax= max(max(result_conv));
-                result_conv=(result_conv/corrmax)*255; %normalize, peak=always 255
+        
+                %result_conv=(result_conv/corrmax); %normalize, peak=always 255
                 %Find the correlation max, at 255
-                [y,x] = find(result_conv==255,1);
+                [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
@@ -1169 +1170 @@
     end
 end
-result_conv=result_conv*corrmax/(255*sum_square);% keep the last correlation matrix for output
+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
+    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
-% OUPUT:
-% vector = optimum displacement vector with subpixel correction
-% F =flag: =0 OK
-%           =-2 , warning: max too close to the edge of the search box (1 pixel margin)
-% INPUT:
-% x,y: position of the maximum correlation at integer values
-
-function [vector,F] = SUBPIXGAUSS (result_conv,x,y)
+function [vector,FF] = SUBPIX2DGAUSS (result_conv,x,y)
 %------------------------------------------------------------------------
 % vector=[0 0]; %default
-F=0;
-[npy,npx]=size(result_conv);
-result_conv(result_conv<1)=1; %set to 1 correlation values smaller than 1  (=0 by discretisation, to avoid divergence in the log)
-%the following 8 lines are copyright (c) 1998, Uri Shavit, Roi Gurka, Alex Liberzon, Technion ??? Israel Institute of Technology
-%http://urapiv.wordpress.com
-peaky = y;peakx=x;
-% if y < npy && y > 1 && x < npx-1 && x > 1
-%     f0 = log(result_conv(y,x));
-%     f1 = log(result_conv(y-1,x));
-%     f2 = log(result_conv(y+1,x));
-%     peaky = peaky+ (f1-f2)/(2*f1-4*f0+2*f2);
-%     f1 = log(result_conv(y,x-1));
-%     f2 = log(result_conv(y,x+1));
-%     peakx = peakx+ (f1-f2)/(2*f1-4*f0+2*f2);
-% else
-%     F=1; % warning flag for vector truncated by the limited search box
-% end
-
-vector=[peakx-floor(npx/2)-1 peaky-floor(npy/2)-1];
-
-%------------------------------------------------------------------------
-% --- Find the maximum of the correlation function after interpolation
-function [vector,F] = SUBPIX2DGAUSS (result_conv,x,y)
-%------------------------------------------------------------------------
-% vector=[0 0]; %default
-F=1;
+FF=true;
 peaky=y;
 peakx=x;
-result_conv(result_conv<1)=1; %set to 1 correlation values smaller than 1 (to avoid divergence in the log)
 [npy,npx]=size(result_conv);
 if (x < npx) && (y < npy) && (x > 1) && (y > 1)
-    F=0;
+    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 ??? 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
+  %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));
@@ -1235 +1242 @@
     c20=(1/6)*sum(sum(c20));
     c02=(1/6)*sum(sum(c02));
-    deltax=(c11*c01-2*c10*c02)/(4*c20*c02-c11^2);
-    deltay=(c11*c10-2*c01*c20)/(4*c20*c02-c11^2);
+    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;

trunk/src/series/extract_rdvision.m

@@ -122 +122 @@
     if exist(filexml,'file')
         [XmlData,errormsg]=imadoc2struct(filexml);
-        if ~isempty(errormsg)
+        if isempty(errormsg)
+            msgbox_uvmat('CONFIRMATION',[filexml ' used' ]);
+        else
             msgbox_uvmat('ERROR',errormsg);
             return

trunk/src/series/stereo_civ.m

@@ -38 +38 @@
 %AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
 
-function [Data,errormsg,result_conv]= stereo_civ(Param)
+function [Data,errormsg,result_conv, XmlData]= stereo_civ(Param)
 Data=[];
 errormsg='';
@@ -85 +85 @@
 WaitbarHandle=findobj(hseries,'Tag','Waitbar');%handle of waitbar in GUI series
 
+%%%%%%%%%%%%%%%% Modif pour test %%%%%%%%%%%%%%%%%%
+% CheckInputFile=isfield(Param,'InputTable');%= 1 in test use for TestCiv (no nc file involved)
+CheckOutputFile = isfield(Param,'OutputSubDir') ; %= 1 in test use for TestPatch (no nc file produced)
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
 %% input files and indexing
 MaxIndex_i=Param.IndexRange.MaxIndex_i;
@@ -95 +100 @@
 time=[];
 for iview=1:size(Param.InputTable,1)
+    %XmlFileName=find_imadoc(Param.InputTable{iview,1},Param.InputTable{iview,2},Param.InputTable{iview,3},Param.InputTable{iview,5});
     XmlFileName=find_imadoc(Param.InputTable{iview,1},Param.InputTable{iview,2});
     if isempty(XmlFileName)
@@ -235 +241 @@
 
 %% Output directory
-OutputDir=[Param.OutputSubDir Param.OutputDirExt];
+OutputDir='/Test';
+if CheckOutputFile
+    OutputDir=[Param.OutputSubDir Param.OutputDirExt];
+end
+
 
 ListGlobalAttribute={'Conventions','Program','CivStage'};
@@ -252 +262 @@
 end
 if isempty(time)% time = index i  by default
-    MaxIndex_i=max(i2_series_Civ1);
-    MaxIndex_j=max(j2_series_Civ1);
+    MaxIndex_i=max(i2_series_Civ1, [], 'all');
+    MaxIndex_j=max(j2_series_Civ1, [], 'all');
     time=(1:MaxIndex_i)'*ones(1,MaxIndex_j);
     time=[zeros(1,MaxIndex_j);time];% insert a first line of zeros
@@ -312 +322 @@
 
         
-        [A,Rangx,Rangy]=phys_ima(A,XmlData,1);
+        [A,Rangx,Rangy]=phys_ima(A,XmlData,1);%transform image A in phys coordinates
         [Npy,Npx]=size(A{1});
         PhysImageA=fullfile_uvmat(RootPath_A,Civ1Dir,RootFile_A,'.png','_1a',i1_series_Civ1(ifield),[],1);
@@ -372 +382 @@
         Data.Civ1_F=reshape(F,[],1);
 
+        %%%%%%%%%%%%%%%% ajout fonction test %%%%%%%
+        if ~isfield (Param.ActionInput,'Civ2') && ~isfield (Param.ActionInput,'Civ3') && ~isfield(Param.ActionInput,'Patch1')
+            [Data.Xmid, Data.Ymid, Data.Uphys, Data.Vphys, Data.Zphys, ...
+                Data.Yphys, Data.Xphys, Data.Error] = getPhysValues(Rangx, ...
+                Rangy, Npx, Npy, Data.Civ1_X, Data.Civ1_Y, Data.Civ1_U, ...
+                Data.Civ1_V, XmlData);
+    
+            if ~isempty(errormsg)
+                disp_uvmat('ERROR',errormsg,checkrun)
+                return
+            end
+        end
+        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
     end
     
@@ -414 +438 @@
             return
         end
-        
-        Data.ListGlobalAttribute=[Data.ListGlobalAttribute {'Patch1_Rho','Patch1_Threshold','Patch1_SubDomain'}];
+
+        %%%%%%%%%%%%%%%%%%%% modif fonction test %%%%%%%%%%%%%%%%%%%%%%%
+        try
+            Data.ListGlobalAttribute=[Data.ListGlobalAttribute {'Patch1_Rho','Patch1_Threshold','Patch1_SubDomain'}];
+        catch
+            Data.ListGlobalAttribute={'Patch1_Rho','Patch1_Threshold','Patch1_SubDomain'};
+        end
+        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+            
         Data.Patch1_FieldSmooth=Param.ActionInput.Patch1.FieldSmooth;
         Data.Patch1_MaxDiff=Param.ActionInput.Patch1.MaxDiff;
         Data.Patch1_SubDomainSize=Param.ActionInput.Patch1.SubDomainSize;
-        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'}}];
+
+        %%%%%%%%%%%%%%%%%%% modif fonction test %%%%%%%%%%%%%%%%%%%%%%%%
+        try
+            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'}}];
+        catch
+            nbvar=0;
+            Data.ListVarName={'Civ1_U_smooth','Civ1_V_smooth','Civ1_SubRange','Civ1_NbCentres','Civ1_Coord_tps','Civ1_U_tps','Civ1_V_tps'};
+            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'}};
+        end
+        
+        if ~isfield(Param.ActionInput,'Civ1')
+            Data.Civ1_X = Param.Civ1_X ;
+            Data.Civ1_Y = Param.Civ1_Y ;
+            Data.Civ1_U = Param.Civ1_U ;
+            Data.Civ1_V = Param.Civ1_V ;
+            Data.Civ1_FF = Param.Civ1_FF ;
+        end
+        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+        
         Data.VarAttribute{nbvar+1}.Role='vector_x';
         Data.VarAttribute{nbvar+2}.Role='vector_y';
@@ -436 +487 @@
         end
         [Data.Civ1_SubRange,Data.Civ1_NbCentres,Data.Civ1_Coord_tps,Data.Civ1_U_tps,Data.Civ1_V_tps,tild,Ures, Vres,tild,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);
+            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;
         Data.Civ1_V_smooth(ind_good)=Vres;
         Data.Civ1_FF(ind_good)=FFres;
-        Data.CivStage=3;
+        Data.CivStage=3; 
+
+        if ~isfield (Param.ActionInput,'Civ2') && ~isfield (Param.ActionInput,'Civ3')
+
+            %%%%%%%%% modif fonction test %%%%%%%%%%%%%%%%%
+            try
+                [Data.Xmid, Data.Ymid, Data.Uphys, Data.Vphys, Data.Zphys, ...
+                    Data.Yphys, Data.Xphys, Data.Error] = getPhysValues(Rangx, ...
+                    Rangy, Npx, Npy, Data.Civ1_X, Data.Civ1_Y, Data.Civ1_U_smooth, ...
+                    Data.Civ1_V_smooth, XmlData);
+            catch % si on a pas Rangx, Rangy, Npx, Npy on ouvre les images pour les recalculer
+                try
+                    ImageName_A=fullfile_uvmat(RootPath_A,SubDir_A,RootFile_A,FileExt_A,NomType_A,i1_series_Civ1(ifield),[],j1_series_Civ1(ifield));
+                    [A{1},VideoObject_A] = read_image(ImageName_A,FileType_A,VideoObject_A,FrameIndex_A_Civ1(ifield));
+                    ImageName_B=fullfile_uvmat(RootPath_B,SubDir_B,RootFile_B,FileExt_B,NomType_B,i2_series_Civ1(ifield),[],j2_series_Civ1(ifield));
+                    [A{2},VideoObject_B] = read_image(ImageName_B,FileType_B,VideoObject_B,FrameIndex_B_Civ1(ifield));
+                catch ME
+                    if ~isempty(ME.message)
+                        disp_uvmat('ERROR', ['error reading input image: ' ME.message],checkrun)
+                        return
+                    end
+                end
+                [A,Rangx,Rangy]=phys_ima(A,XmlData,1);%transform image A in phys coordinates
+                [Npy,Npx]=size(A{1});
+                [Data.Xmid, Data.Ymid, Data.Uphys, Data.Vphys, Data.Zphys, ...
+                    Data.Yphys, Data.Xphys, Data.Error] = getPhysValues(Rangx, ...
+                    Rangy, Npx, Npy, Data.Civ1_X, Data.Civ1_Y, Data.Civ1_U_smooth, ...
+                    Data.Civ1_V_smooth, XmlData);
+            end              
+            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    
+            if ~isempty(errormsg)
+                disp_uvmat('ERROR',errormsg,checkrun)
+                return
+            end
+        end
                
     end
+
+    
+    
     
     %% Civ2
@@ -553 +643 @@
         Data.Civ2_F=reshape(F,[],1);
         Data.CivStage=Data.CivStage+1;
+
+        %%%%%%%%%%%%%%%% ajout fonction test %%%%%%%
+        if ~isfield (Param.ActionInput,'Civ3') && ~isfield(Param.ActionInput,'Patch2')
+            [Data.Xmid, Data.Ymid, Data.Uphys, Data.Vphys, Data.Zphys, ...
+                Data.Yphys, Data.Xphys, Data.Error] = getPhysValues(Rangx, ...
+                Rangy, Npx, Npy, Data.Civ2_X, Data.Civ2_Y, Data.Civ2_U, ...
+                Data.Civ2_V, XmlData);
+    
+            if ~isempty(errormsg)
+                disp_uvmat('ERROR',errormsg,checkrun)
+                return
+            end
+        end
+        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     end
     
@@ -612 +716 @@
         Data.Civ2_FF(ind_good)=FFres;
         Data.CivStage=Data.CivStage+1;
+
+
+        if ~isfield (Param.ActionInput,'Civ3')
+            [Data.Xmid, Data.Ymid, Data.Uphys, Data.Vphys, Data.Zphys, ...
+                Data.Yphys, Data.Xphys, Data.Error] = getPhysValues(Rangx, ...
+                Rangy, Npx, Npy, Data.Civ2_X, Data.Civ2_Y, Data.Civ2_U_smooth, ...
+                Data.Civ2_V_smooth, XmlData);
+    
+            if ~isempty(errormsg)
+                disp_uvmat('ERROR',errormsg,checkrun)
+                return
+            end
+        end
+
     end
     
@@ -769 +887 @@
         Data.CivStage=Data.CivStage+1;
            
-         % get z from u and v (displacements)     
-        Data.Xmid=Rangx(1)+(Rangx(2)-Rangx(1))*(Data.Civ3_X-0.5)/(Npx-1);%temporary coordinate (velocity taken at the point middle from imgae 1 and 2)
-        Data.Ymid=Rangy(2)+(Rangy(1)-Rangy(2))*(Data.Civ3_Y-0.5)/(Npy-1);%temporary coordinate (velocity taken at the point middle from imgae 1 and 2)
-        Data.Uphys=Data.Civ3_U_smooth*(Rangx(2)-Rangx(1))/(Npx-1);
-        Data.Vphys=Data.Civ3_V_smooth*(Rangy(1)-Rangy(2))/(Npy-1);
-        [Data.Zphys,Data.Xphys,Data.Yphys,Data.Error]=shift2z(Data.Xmid,Data.Ymid,Data.Uphys,Data.Vphys,XmlData); %Data.Xphys and Data.Xphys are real coordinate (geometric correction more accurate than xtemp/ytemp)
+        
+        [Data.Xmid, Data.Ymid, Data.Uphys, Data.Vphys, Data.Zphys, ...
+            Data.Yphys, Data.Xphys, Data.Error] = getPhysValues(Rangx, ...
+            Rangy, Npx, Npy, Data.Civ3_X, Data.Civ3_Y, Data.Civ3_U_smooth, ...
+            Data.Civ3_V_smooth, XmlData);
+
         if ~isempty(errormsg)
             disp_uvmat('ERROR',errormsg,checkrun)
@@ -784 +902 @@
     
     %% write result in a netcdf file if requested
+
     if LSM ~= 1 % store all data
-        if exist('ncfile','var')
-            errormsg=struct2nc(ncfile,Data);
+        if exist('ncfile2','var')
+            errormsg=struct2nc(ncfile2,Data);
             if isempty(errormsg)
-                disp([ncfile ' written'])
+                disp([ncfile2 ' written'])
             else
                 disp(errormsg)
@@ -794 +913 @@
         end
     else
+        
+        
        % store only phys data
-        Data_light.ListVarName={'Xphys','Yphys','Zphys','Civ3_C','DX','DY','Error'};
+        Data_light.ListVarName={'Xphys','Yphys','Zphys','Civ_C','DX','DY','Error'};
         Data_light.VarDimName={'nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3'};
         Data_light.VarAttribute{1}.Role='coord_x';
@@ -802 +923 @@
          Data_light.VarAttribute{5}.Role='vector_x';
          Data_light.VarAttribute{6}.Role='vector_y';
-        ind_good=find(Data.Civ3_FF==0);
-        Data_light.Zphys=Data.Zphys(ind_good);
-        Data_light.Yphys=Data.Yphys(ind_good);
-        Data_light.Xphys=Data.Xphys(ind_good);
-        Data_light.Civ3_C=Data.Civ3_C(ind_good);
-        Data_light.DX=Data.Uphys(ind_good);
-        Data_light.DY=Data.Vphys(ind_good);
-        Data_light.Error=Data.Error(ind_good);
+        
+         % vrifie la dernire passe effectue
+        if isfield (Param.ActionInput,'Civ3')
+            try
+                ind_good=find(Data.Civ3_FF==0);
+                Data_light.Civ_C=Data.Civ3_C(ind_good);
+            catch
+                Data_light.Civ_C=Data.Civ3_C;
+            end
+        elseif isfield (Param.ActionInput,'Civ2')
+            try
+                ind_good=find(Data.Civ2_FF==0);
+                Data_light.Civ_C=Data.Civ2_C(ind_good);
+            catch
+                Data_light.Civ_C=Data.Civ2_C;
+            end
+        elseif isfield (Param.ActionInput,'Civ1')
+            try
+                ind_good=find(Data.Civ1_FF==0);
+                Data_light.Civ_C=Data.Civ1_C(ind_good);
+            catch
+                Data_light.Civ_C=Data.Civ1_C;
+            end
+        end
+        
+        try
+            Data_light.Zphys=Data.Zphys(ind_good);
+            Data_light.Yphys=Data.Yphys(ind_good);
+            Data_light.Xphys=Data.Xphys(ind_good);
+            Data_light.DX=Data.Uphys(ind_good);
+            Data_light.DY=Data.Vphys(ind_good);
+            Data_light.Error=Data.Error(ind_good);
+        catch
+            Data_light.Zphys=Data.Zphys;
+            Data_light.Yphys=Data.Yphys;
+            Data_light.Xphys=Data.Xphys;
+            Data_light.DX=Data.Uphys;
+            Data_light.DY=Data.Vphys;
+            Data_light.Error=Data.Error;
+        end
+
+
        if exist('ncfile2','var')
             errormsg=struct2nc(ncfile2,Data_light);
@@ -827 +982 @@
 end
 
-
+end
 
 % 'civ': function piv.m adapted from PIVlab http://pivlab.blogspot.com/
@@ -1068 +1223 @@
 end
 result_conv=result_conv*corrmax/(255*sum_square);% keep the last correlation matrix for output
+end
 
 %------------------------------------------------------------------------
@@ -1098 +1254 @@
 end
 vector=[peakx-floor(npx/2)-1 peaky-floor(npy/2)-1];
+end
 
 %------------------------------------------------------------------------
@@ -1141 +1298 @@
 end
 vector=[peakx-floor(npx/2)-1 peaky-floor(npy/2)-1];
+end
 
 %'RUN_FIX': function for fixing velocity fields:
@@ -1189 +1347 @@
 end
 
+end
 
 %------------------------------------------------------------------------
@@ -1250 +1409 @@
         j2_series=ind2*ones(1,size(i_series,2));
         check_bounds=zeros(size(i1_series));% no limitations due to min-max indices
+    end
 end
 
@@ -1258 +1418 @@
 % ymid+ v/2: set of apparent phys y coordinates in the ref plane, image B
 % XmlData: content of the xml files containing geometric calibration parameters
+
+
 function [z,Xphy,Yphy,Error]=shift2z(xmid, ymid, u, v,XmlData)
 z=0;
@@ -1320 +1482 @@
 Xphy=mean(xnew,1);
 Yphy=mean(ynew,1);
-
-
-
-
-
+end
+
+function [Xmid, Ymid, Uphys, Vphys, Zphys, Yphys, Xphys, Error] ...
+    = getPhysValues(Rangx, Rangy, Npx, Npy, Data_Civ_X, Data_Civ_Y, ...
+    Data_Civ_U_smooth, Data_Civ_V_smooth, XmlData)
+    
+    % get z from u and v (displacements)     
+    Xmid=Rangx(1)+(Rangx(2)-Rangx(1))*(Data_Civ_X-0.5)/(Npx-1);%temporary coordinate (velocity taken at the point middle from imgae 1 and 2)
+    Ymid=Rangy(2)+(Rangy(1)-Rangy(2))*(Data_Civ_Y-0.5)/(Npy-1);%temporary coordinate (velocity taken at the point middle from imgae 1 and 2)
+    Uphys=Data_Civ_U_smooth*(Rangx(2)-Rangx(1))/(Npx-1);
+    Vphys=Data_Civ_V_smooth*(Rangy(1)-Rangy(2))/(Npy-1);
+    [Zphys,Xphys,Yphys,Error]=shift2z(Xmid,Ymid,Uphys,Vphys,XmlData); %Data.Xphys and Data.Xphys are real coordinate (geometric correction more accurate than xtemp/ytempy
+
+end

trunk/src/series/stereo_input.m

@@ -1797 +1797 @@
 set(handles.configSource,'String','NEW')
 set(handles.OK,'BackgroundColor',[1 0 1])
+
+
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%%%%%%%%%%%%%   TEST functions
@@ -1814 +1816 @@
      Param.Action.RUN=1;
      Param.ActionInput=read_GUI(handles.civ_input);
-     if isfield(Param.ActionInput,'Fix1')
-         Param.ActionInput=rmfield(Param.ActionInput,'Fix1');
-     end
+
      if isfield(Param.ActionInput,'Patch1')
          Param.ActionInput=rmfield(Param.ActionInput,'Patch1');
@@ -1829 +1829 @@
          Param.ActionInput=rmfield(Param.ActionInput,'Patch2');
      end
-     if isfield(Param,'OutputSubDir')
-     Param=rmfield(Param,'OutputSubDir'); %remove output file option from civ_series
+     if isfield(Param.ActionInput,'Civ3')%remove options that may be selected beyond Patch1
+         Param.ActionInput=rmfield(Param.ActionInput,'Civ3');
      end
+     if isfield(Param.ActionInput,'Fix3')
+         Param.ActionInput=rmfield(Param.ActionInput,'Fix3');
+     end
+     if isfield(Param.ActionInput,'Patch3')
+         Param.ActionInput=rmfield(Param.ActionInput,'Patch3');
+     end
+%      if isfield(Param,'OutputSubDir')
+%         Param=rmfield(Param,'OutputSubDir'); %remove output file option from civ_series
+%      end
      Param.ActionInput.Civ1.CorrSmooth=0;% launch Civ1 with no data point (to get the image names for A and B)
-     [Data,errormsg]=civ_series(Param);% get the civ1+fix1 results 
-     if ~isempty(errormsg), return, end % rmq: error msg displayed in civ_series
+     [Data,errormsg, ~, xmlData]=stereo_civ(Param);% get the civ1+fix1 results 
+     % if ~isempty(errormsg), return, end % rmq: error msg displayed in civ_series
      
  %% create image data ImageData for display
      ImageData.ListVarName={'ny','nx','A'};
      ImageData.VarDimName= {'ny','nx',{'ny','nx'}};
-     ImageData.A=imread(Data.Civ1_ImageA); % read the first image
+
+     %%%%%%%%%%%%%%%%% modif fonction test %%%%%%%%%%%
+     ImageData.VarAttribute{1}.Role='coord_y';
+     ImageData.VarAttribute{2}.Role='coord_x';
+     ImageData.VarAttribute{3}.Role='scalar';
+
+     A{1}=imread(Data.Civ1_ImageA); % read the first image
+     A{2}=imread(Data.Civ1_ImageB); % read the first image
+
+     phys_img = phys_ima(A,xmlData,1);%transform image A in phys coordinates
+     ImageData.A = phys_img{1};
+     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
      if ndims(ImageData.A)==3 %case of color image
          ImageData.VarDimName= {'ny','nx',{'ny','nx','rgb'}};
@@ -1856 +1877 @@
      ViewData.PlotAxes.X=Data.Civ1_X';
      ViewData.PlotAxes.Y=Data.Civ1_Y';
-     ViewData.PlotAxes.B=imread(Data.Civ1_ImageB);%store the second image in the UserData of the GUI view_field
+
+     %%%%%%%%%%%%%%%%%%%%%%%
+     ViewData.PlotAxes.B=phys_img{2};%store the second image in the UserData of the GUI view_field
+     %%%%%%%%%%%%%%%%%%%%%%%
+
      set(hview_field,'UserData',ViewData)% store the info in the UserData of image view_field
     
@@ -1897 +1922 @@
      Param.Action.RUN=1;
      Param.ActionInput=read_GUI(handles.civ_input);
+
      if isfield(Param.ActionInput,'Civ2')%remove options that may be selected beyond Patch1
          Param.ActionInput=rmfield(Param.ActionInput,'Civ2');
@@ -1906 +1932 @@
          Param.ActionInput=rmfield(Param.ActionInput,'Patch2');
      end
-     if isfield(Param,'OutputSubDir')
-     Param=rmfield(Param,'OutputSubDir'); %remove output file option from civ_series
+     if isfield(Param.ActionInput,'Civ3')%remove options that may be selected beyond Patch1
+         Param.ActionInput=rmfield(Param.ActionInput,'Civ3');
      end
+     if isfield(Param.ActionInput,'Fix3')
+         Param.ActionInput=rmfield(Param.ActionInput,'Fix3');
+     end
+     if isfield(Param.ActionInput,'Patch3')
+         Param.ActionInput=rmfield(Param.ActionInput,'Patch3');
+     end
+     % if isfield(Param,'OutputSubDir')
+     % Param=rmfield(Param,'OutputSubDir'); %remove output file option from civ_series
+     % end
+
      ParamPatch1=Param.ActionInput.Patch1; %store the patch1 parameters
      Param.ActionInput=rmfield(Param.ActionInput,'Patch1');% does not execute Patch
-     [Data,errormsg]=civ_series(Param);% get the civ1+fix1 results
+
+     [Data,errormsg]=stereo_civ(Param);% get the civ1+fix1 results
      bckcolor=get(handles.civ_input,'Color');
      set(handles.Civ1,'BackgroundColor',bckcolor)% indicate civ1 calculation is finished
@@ -1921 +1958 @@
      Param.Civ1_V=Data.Civ1_V;
      Param.Civ1_FF=Data.Civ1_FF;
-     Param=rmfield(Param,'InputTable');%desactivate input file reading
+
+     %%%%% modif fonction test %%%%%%
+     % Param=rmfield(Param,'InputTable');%desactivate input file reading
+     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
     if isfield(Param.ActionInput,'Civ1')
         Param.ActionInput=rmfield(Param.ActionInput,'Civ1');%desactivate civ1: remove civ1 input param if relevant
@@ -1935 +1976 @@
     for irho=1:7
         Param.ActionInput.Patch1.FieldSmooth=SmoothingParam(irho);
-        [Data,errormsg]= civ_series(Param);%apply the processing fct
-        if ~isempty(errormsg)
-            msgbox_uvmat('ERROR',errormsg)
-            return
-        end
+        [Data,errormsg]= stereo_civ(Param);%apply the processing fct
+
+        %%%%%%%%%% modif fonction test %%%%%%%%%%%%%% 
+        % if ~isempty(errormsg)
+        %     msgbox_uvmat('ERROR',errormsg)
+        %     return
+        % end
+        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
         ind_good=find(Data.Civ1_FF==0);
         Civ1_U_Diff=Data.Civ1_U(ind_good)-Data.Civ1_U_smooth(ind_good);
@@ -1976 +2021 @@
      Param.Action.RUN=1;
      Param.ActionInput=read_GUI(handles.civ_input);
-     if isfield(Param,'OutputSubDir')
-     Param=rmfield(Param,'OutputSubDir'); %remove output file option from civ_series
+
+     if isfield(Param.ActionInput,'Patch2')
+         Param.ActionInput=rmfield(Param.ActionInput,'Patch2');
      end
+     if isfield(Param.ActionInput,'Civ3')%remove options that may be selected beyond Patch1
+         Param.ActionInput=rmfield(Param.ActionInput,'Civ3');
+     end
+     if isfield(Param.ActionInput,'Fix3')
+         Param.ActionInput=rmfield(Param.ActionInput,'Fix3');
+     end
+     if isfield(Param.ActionInput,'Patch3')
+         Param.ActionInput=rmfield(Param.ActionInput,'Patch3');
+     end
+
+     % if isfield(Param,'OutputSubDir')
+     % Param=rmfield(Param,'OutputSubDir'); %remove output file option from civ_series
+     % end
+
      Param.ActionInput.Civ2.CorrSmooth=0;% launch Civ2 with no data point (to get the image names for A and B)
      set(handles.Civ1,'BackgroundColor',[1 1 0])
      set(handles.Fix1,'BackgroundColor',[1 1 0])
      set(handles.Patch1,'BackgroundColor',[1 1 0])
-     [Data,errormsg]=civ_series(Param);% get the civ1+fix1 results
+     [Data,errormsg, ~, xmlData]=stereo_civ(Param);% get the civ1+fix1+patch1+civ2+fix2 results
+     % if ~isempty(errormsg), return, end
      
      %% create image data ImageData for display
      ImageData.ListVarName={'ny','nx','A'};
      ImageData.VarDimName= {'ny','nx',{'ny','nx'}};
-     ImageData.A=imread(Data.Civ2_ImageA); % read the first image
+
+     %%%%%%%%%%%%%%%%% modif fonction test %%%%%%%%%%%
+     ImageData.VarAttribute{1}.Role='coord_y';
+     ImageData.VarAttribute{2}.Role='coord_x';
+     ImageData.VarAttribute{3}.Role='scalar';
+
+     A{1}=imread(Data.Civ2_ImageA); % read the first image
+     A{2}=imread(Data.Civ2_ImageB); % read the first image
+
+     phys_img = phys_ima(A,xmlData,1);%transform image A in phys coordinates
+     ImageData.A = phys_img{1};
+     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
      if ndims(ImageData.A)==3 %case of color image
          ImageData.VarDimName= {'ny','nx',{'ny','nx','rgb'}};
@@ -2012 +2085 @@
     % store info in the UserData of view-field
     ViewData.CivHandle=handles.civ_input;% indicate the handle of the civ GUI in view_field
-    ViewData.PlotAxes.B=imread(Data.Civ2_ImageB);%store the second image in the UserData of the GUI view_field
+
+    %%%%%%%%%%%% modif fonction test %%%%%%%%%%%%%%%%%%%%%%%%%%
+    ViewData.PlotAxes.B=phys_img{2};%store the second image in the UserData of the GUI view_field
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    
     ViewData.PlotAxes.X=Data.Civ2_X';
     ViewData.PlotAxes.Y=Data.Civ2_Y';
@@ -2038 +2115 @@
         set(corrfig,'name','image correlation')
         set(corrfig,'DeleteFcn',{@closeview_field})%
-        set(handles.TestCiv1,'BackgroundColor',[1 0 0])
+        set(handles.TestCiv2,'BackgroundColor',[1 0 0])
     else
-        set(handles.TestCiv1,'BackgroundColor',[1 0 0])% paint button to red
+        set(handles.TestCiv2,'BackgroundColor',[1 0 0])% paint button to red
         corrfig=findobj(allchild(0),'tag','corrfig');% look for a current figure for image correlation display
         if ~isempty(corrfig)

trunk/src/series/test_filter_tps.m

@@ -284 +284 @@
 %plot rms difference and proportion of excluded vectors
 figure(2)
-cla
+clf
 if CivStage==3% civ1
     ref=0.2; %recommanded value for diff rms

trunk/src/tps_coeff.m

@@ -19 +19 @@
 
 % INPUT:
-%  ctrs: NxNbDim matrix  representing the positions of the N centers, sources of the tps (NbDim=space dimension)
+%  ctrs: NxNbDim matrix  representing the positions of the N centers, sources of the tps (NbDim=space dimension, N>=3)
 %  U: Nx1 column vector representing the values of the considered scalar measured at the centres ctrs
 %  smoothing: smoothing parameter: the result is smoother for larger smoothing.

trunk/src/uvmat.m

@@ -2139 +2139 @@
     [RootPath,SubDir,RootFile,FileIndices,FileExt]=read_file_boxes_1(handles);
     % detect the file type, get the movie object if relevant, and look for the corresponding file series:
-    [RootPath,SubDir,RootFile,i1_series,i2_series,j1_series,j2_series,tild,FileInfo,MovieObject]=find_file_series(fullfile(RootPath,SubDir),[RootFile FileIndices FileExt]);
+    [RootPath,SubDir,~,i1_series,i2_series,j1_series,j2_series,~,FileInfo,MovieObject]=find_file_series(fullfile(RootPath,SubDir),[RootFile FileIndices FileExt]);
     if isempty(i1_series)
         fileinput=uigetfile_uvmat('pick an input file for the second line',fullfile(RootPath,SubDir));
@@ -2154 +2154 @@
         end
     else
-        % initiate the input file series and inputfilerefresh the current field view:
+        % initiate the input file series and refresh the current field view:
         errormsg=update_rootinfo(handles,i1_series,i2_series,j1_series,j2_series,FileInfo,MovieObject,2);
     end
@@ -2494 +2494 @@
 %FileType=FileInfo.FileType;
 if isfield(XmlData,'Time')&& ~isempty(XmlData.Time) && ...
-        (strcmp(FileInfo.FileType,'image')|| strcmp(FileInfo.FileType,'multimage'))%||strcmp(FileType,'civdata')||strcmp(FileType,'civx'))
+        (strcmp(FileInfo.FileType,'image')|| strcmp(FileInfo.FileType,'multimage'))
     TimeName='xml';
 elseif strcmp(FileInfo.FieldType,'civdata')% ajouter pivdata_fluidimage
@@ -3418 +3418 @@
 
 %------------------------------------------------------------------------
-% --- Executes on button press in InputFileREFRESH.
+% --- Executes on button press in InputFileREFRESH. 
 function REFRESH_Callback(hObject, eventdata, handles)
 %------------------------------------------------------------------------
@@ -3470 +3470 @@
 
 %% initialisation
+errormsg='';
 pointer=get(handles.uvmat,'Pointer');
 if strcmp(pointer,'watch')% reinitialise the mouse if stuck to 'watch'
@@ -3486 +3487 @@
 
 %% determine the main input file information for action
-if isempty(regexp(FileName,'^http://')) &&~exist(FileName,'file')
+if isempty(regexp(FileName,'^http://', 'once')) &&~exist(FileName,'file')
     errormsg=['input file ' FileName ' does not exist'];
     return
 end
 NomType=get(handles.NomType,'String');
-NomType_1='';
-if strcmp(get(handles.NomType_1,'Visible'),'on')
-    NomType_1=get(handles.NomType_1,'String');
-end
+% NomType_1='';
+% if strcmp(get(handles.NomType_1,'Visible'),'on')
+%     NomType_1=get(handles.NomType_1,'String');
+% end
 %update the z position index
 mode_slice=get(handles.slices,'String');
@@ -3501 +3502 @@
     set(handles.z_index,'String',num2str(z_index))
 else
-    nbslice=str2num(get(handles.num_NbSlice,'String'));
+    nbslice=str2double(get(handles.num_NbSlice,'String'));
     z_index=mod(num_i1-1,nbslice)+1;
     set(handles.z_index,'String',num2str(z_index))
@@ -3556 +3557 @@
         FieldName= list_fields{get(handles.FieldName,'Value')}; % selected field
         if strcmp(FieldName,'get_field...')
-            FieldName_Callback(hObject, eventdata, handles)
+            FieldName_Callback(get(handles.FieldName), [], handles)
             return
         else