Index: trunk/src/series/civ2vel_3C.m
===================================================================
--- trunk/src/series/civ2vel_3C.m	(revision 876)
+++ trunk/src/series/civ2vel_3C.m	(revision 878)
@@ -169,5 +169,5 @@
 %% MAIN LOOP ON FIELDS
 warning off
-if NbField<2 && length(filecell{:,1})>2
+if NbField<2 
     disp_uvmat('ERROR','you need at least 2 images to compute the mean position for the stereo.',checkrun)
 return
@@ -285,11 +285,10 @@
 end
     ZI=mean(ZItemp,3); %mean between two the two time step
+    Vtest=ZItemp(:,:,2)-ZItemp(:,:,1);
     
     [Xa,Ya]=px_XYZ(XmlData{1}.GeometryCalib,XI,YI,ZI);% set of image coordinates on view a
     [Xb,Yb]=px_XYZ(XmlData{2}.GeometryCalib,XI,YI,ZI);% set of image coordinates on view b
     
-    %trouver z
-    % trouver les coordonnées px sur chaque image.
-    %A=
+   
     for iview=1:2
         %% reading input file(s)
@@ -322,7 +321,9 @@
     Ua=griddata(X1,Y1,U1,Xa,Ya);
     Va=griddata(X1,Y1,V1,Xa,Ya);
-    A=get_coeff(XmlData{1}.GeometryCalib,Xa,Ya,YI,YI,ZI);
-    
-    
+    
+    [Ua,Va,Xa,Ya]=Ud2U(XmlData{1}.GeometryCalib,Xa,Ya,Ua,Va); % convert Xd data to X 
+    [A]=get_coeff(XmlData{1}.GeometryCalib,Xa,Ya,XI,YI,ZI); %get coef A~
+    
+    %remove wrong vector
     temp=find(Data{2}.FF==0);
     X2=Data{2}.X(temp);
@@ -332,7 +333,13 @@
     Ub=griddata(X2,Y2,U2,Xb,Yb);
     Vb=griddata(X2,Y2,V2,Xb,Yb);
-    B=get_coeff(XmlData{2}.GeometryCalib,Xb,Yb,YI,YI,ZI);
+
+    [Ub,Vb,Xb,Yb]=Ud2U(XmlData{2}.GeometryCalib,Xb,Yb,Ub,Vb); % convert Xd data to X 
+    [B]=get_coeff(XmlData{2}.GeometryCalib,Xb,Yb,XI,YI,ZI); %get coef B~
+   
+    
+    % System to solve
     S=ones(size(XI,1),size(XI,2),3);
     D=ones(size(XI,1),size(XI,2),3,3);
+
     S(:,:,1)=A(:,:,1,1).*Ua+A(:,:,2,1).*Va+B(:,:,1,1).*Ub+B(:,:,2,1).*Vb;
     S(:,:,2)=A(:,:,1,2).*Ua+A(:,:,2,2).*Va+B(:,:,1,2).*Ub+B(:,:,2,2).*Vb;
@@ -349,5 +356,5 @@
     for indj=1:size(XI,1)
         for indi=1:size(XI,2)
-            dxyz=squeeze(S(indj,indi,:))\squeeze(D(indj,indi,:,:));
+            dxyz=(squeeze(D(indj,indi,:,:))*1000)\(squeeze(S(indj,indi,:))*1000); % solving...
             U(indj,indi)=dxyz(1);
             V(indj,indi)=dxyz(2);
@@ -360,4 +367,7 @@
     Error(:,:,3)=B(:,:,1,1).*U+B(:,:,1,2).*V+B(:,:,1,3).*W-Ub;
     Error(:,:,4)=B(:,:,2,1).*U+B(:,:,2,2).*V+B(:,:,2,3).*W-Vb;
+    
+    
+
     
     %% generating the name of the merged field
@@ -396,5 +406,8 @@
     MergeData.V=V/Dt;
     MergeData.W=W/Dt;
-    MergeData.Error=sqrt(sum(Error.*Error,3));
+    
+    mfx=(XmlData{1}.GeometryCalib.fx_fy(1)+XmlData{2}.GeometryCalib.fx_fy(1))/2;
+    mfy=(XmlData{1}.GeometryCalib.fx_fy(2)+XmlData{2}.GeometryCalib.fx_fy(2))/2;
+    MergeData.Error=(sqrt(mfx^2+mfy^2)/4).*sqrt(sum(Error.*Error,3));
     errormsg=struct2nc(OutputFile,MergeData);%save result file
     if isempty(errormsg)
@@ -406,8 +419,9 @@
 
 
-function A=get_coeff(Calib,X,Y,x,y,z)
+function [A]=get_coeff(Calib,X,Y,x,y,z) % compute A~ coefficients 
 R=(Calib.R)';%rotation matrix
 T_z=Calib.Tx_Ty_Tz(3);
 T=R(7)*x+R(8)*y+R(9)*z+T_z;
+
 A(:,:,1,1)=(R(1)-R(7)*X)./T;
 A(:,:,1,2)=(R(2)-R(8)*X)./T;
@@ -417,8 +431,24 @@
 A(:,:,2,3)=(R(6)-R(9)*Y)./T;
 
-
-
-
-function [z,Xphy,Yphy,error]=shift2z(xmid, ymid, u, v,XmlData)
+function [U,V,X,Y]=Ud2U(Calib,Xd,Yd,Ud,Vd) % convert Xd to X  and Ud to U
+
+X1d=Xd-Ud/2;
+X2d=Xd+Ud/2;
+Y1d=Yd-Vd/2;
+Y2d=Yd+Vd/2;
+
+X1=(X1d-Calib.Cx_Cy(1))./Calib.fx_fy(1).*(1 + Calib.kc.*Calib.fx_fy(1).^(-2).*(X1d-Calib.Cx_Cy(1)).^2 + Calib.kc.*Calib.fx_fy(2).^(-2).*(Y1d-Calib.Cx_Cy(2)).^2 ).^(-1);
+X2=(X2d-Calib.Cx_Cy(1))./Calib.fx_fy(1).*(1 + Calib.kc.*Calib.fx_fy(1).^(-2).*(X2d-Calib.Cx_Cy(1)).^2 + Calib.kc.*Calib.fx_fy(2).^(-2).*(Y2d-Calib.Cx_Cy(2)).^2 ).^(-1);
+Y1=(Y1d-Calib.Cx_Cy(2))./Calib.fx_fy(2).*(1 + Calib.kc.*Calib.fx_fy(1).^(-2).*(X1d-Calib.Cx_Cy(1)).^2 + Calib.kc.*Calib.fx_fy(2).^(-2).*(Y1d-Calib.Cx_Cy(2)).^2 ).^(-1);
+Y2=(Y2d-Calib.Cx_Cy(2))./Calib.fx_fy(2).*(1 + Calib.kc.*Calib.fx_fy(1).^(-2).*(X2d-Calib.Cx_Cy(1)).^2 + Calib.kc.*Calib.fx_fy(2).^(-2).*(Y2d-Calib.Cx_Cy(2)).^2 ).^(-1);
+
+U=X2-X1;
+V=Y2-Y1;
+X=X1+U/2;
+Y=Y1+V/2;
+
+
+
+function [z,Xphy,Yphy,error]=shift2z(xmid, ymid, u, v,XmlData) % get H from stereo data
 z=0;
 error=0;
@@ -429,5 +459,5 @@
 R=(Calib_A.R)';
 x_a=xmid- u/2;
-y_a=ymid- u/2;
+y_a=ymid- v/2; 
 z_a=R(7)*x_a+R(8)*y_a+Calib_A.Tx_Ty_Tz(1,3);
 Xa=(R(1)*x_a+R(2)*y_a+Calib_A.Tx_Ty_Tz(1,1))./z_a;
@@ -445,6 +475,10 @@
 
 %% second image
+%loading shift angle
+
 Calib_B=XmlData{2}.GeometryCalib;
 R=(Calib_B.R)';
+
+
 x_b=xmid+ u/2;
 y_b=ymid+ v/2;
@@ -464,5 +498,11 @@
 %% result
 Den=(Dxb-Dxa).*(Dxb-Dxa)+(Dyb-Dya).*(Dyb-Dya);
-error=((Dyb-Dya).*u-(Dxb-Dxa).*v)./Den;
+error=((Dyb-Dya).*(-u)-(Dxb-Dxa).*(-v))./Den;
+% ex=-error.*(Dyb-Dya);
+% ey=-error.*(Dxb-Dxa);
+
+% z1=-u./(Dxb-Dxa);
+% z2=-v./(Dyb-Dya);
+z=((Dxb-Dxa).*(-u)+(Dyb-Dya).*(-v))./Den;
 
 xnew(1,:)=Dxa.*z+x_a;
@@ -470,8 +510,8 @@
 ynew(1,:)=Dya.*z+y_a;
 ynew(2,:)=Dyb.*z+y_b;
-
-Xphy=mean(xnew,1); %on moyenne les 2 valeurs 
-Yphy=mean(ynew,1);
-z=((Dxb-Dxa).*u-(Dyb-Dya).*v)./Den;
-
-
+Xphy=mean(xnew,1);
+Yphy=mean(ynew,1); 
+
+
+
+
Index: trunk/src/series/civ_series.m
===================================================================
--- trunk/src/series/civ_series.m	(revision 876)
+++ trunk/src/series/civ_series.m	(revision 878)
@@ -324,4 +324,5 @@
 %%%%% MAIN LOOP %%%%%%
 maskname='';% initiate the mask name
+tic;
 for ifield=1:NbField
     if ~isempty(RUNHandle)% update the waitbar in interactive mode with GUI series  (checkrun=1)
@@ -852,4 +853,5 @@
     end
 end
+disp(['ellapsed time ' num2str(toc) ' s'])
 
 % 'civ': function piv.m adapted from PIVlab http://pivlab.blogspot.com/
@@ -1024,14 +1026,11 @@
                 XIant=XI-par_civ.DUDX(ivec)*XI-par_civ.DUDY(ivec)*YI+ceil(size(image1_crop,2)/2);
                 YIant=YI-par_civ.DVDX(ivec)*XI-par_civ.DVDY(ivec)*YI+ceil(size(image1_crop,1)/2);
-                        image1_crop_1=image1_crop;
                 image1_crop=interp2(image1_crop,XIant,YIant);
-                        image1_crop_2=image1_crop;
                 image1_crop(isnan(image1_crop))=0;
-                  image1_crop_3=image1_crop;
                 xi=(1:mesh:size(image2_crop,2));
                 yi=(1:mesh:size(image2_crop,1))';
-                image2_crop=interp2(image2_crop,xi,yi);
+                image2_crop=interp2(image2_crop,xi,yi,'*spline');
                 image2_crop(isnan(image2_crop))=0;
-                par_civ.CorrSmooth=3;
+                %par_civ.CorrSmooth=3;%%%%%%%%%%%%%%%%%%%
                 %%
             end
Index: trunk/src/series/stereo_civ.m
===================================================================
--- trunk/src/series/stereo_civ.m	(revision 876)
+++ trunk/src/series/stereo_civ.m	(revision 878)
@@ -1,3 +1,3 @@
-%'stereo_civ': determination of topography by image correlation of two stereo views
+%'civ_series': PIV function activated by the general GUI series
 % --- call the sub-functions:
 %   civ: PIV function itself
@@ -282,5 +282,6 @@
        
     %% Civ1
-   
+
+    
     % if Civ1 computation is requested
     if isfield (Param.ActionInput,'Civ1')
@@ -297,4 +298,7 @@
             end
         end
+        
+
+        
         [A,Rangx,Rangy]=phys_ima(A,XmlData,1);
         [Npy,Npx]=size(A{1});
@@ -516,15 +520,13 @@
             mask=imread(par_civ2.Mask);
         end
-        ibx2=ceil(par_civ2.CorrBoxSize(1)/2);
-        iby2=ceil(par_civ2.CorrBoxSize(2)/2);
-%         par_civ2.SearchBoxSize(1)=2*ibx2+9;% search ara +-4 pixels around the guess
-%         par_civ2.SearchBoxSize(2)=2*iby2+9;
+%         ibx2=ceil(par_civ2.CorrBoxSize(1)/2);
+%         iby2=ceil(par_civ2.CorrBoxSize(2)/2);
         par_civ2.SearchBoxShift=[Shiftx(nbval>=1)./nbval(nbval>=1) Shifty(nbval>=1)./nbval(nbval>=1)];
         par_civ2.Grid=[par_civ2.Grid(nbval>=1,1)-par_civ2.SearchBoxShift(:,1)/2 par_civ2.Grid(nbval>=1,2)-par_civ2.SearchBoxShift(:,2)/2];% grid taken at the extrapolated origin of the displacement vectors
         if par_civ2.CheckDeformation
-            par_civ2.DUDX=DUDX./nbval;
-            par_civ2.DUDY=DUDY./nbval;
-            par_civ2.DVDX=DVDX./nbval;
-            par_civ2.DVDY=DVDY./nbval;
+            par_civ2.DUDX=DUDX(nbval>=1)./nbval(nbval>=1);
+            par_civ2.DUDY=DUDY(nbval>=1)./nbval(nbval>=1);
+            par_civ2.DVDX=DVDX(nbval>=1)./nbval(nbval>=1);
+            par_civ2.DVDY=DVDY(nbval>=1)./nbval(nbval>=1);
         end
         % calculate velocity data (y and v in indices, reverse to y component)
@@ -617,8 +619,4 @@
         Data.Civ2_FF(ind_good)=FFres;
         Data.CivStage=Data.CivStage+1;
-        
-            
- 
-        
     end
     
@@ -659,7 +657,5 @@
             [GridX,GridY]=meshgrid(minix:par_civ3.Dx:maxix,miniy:par_civ3.Dy:maxiy);
             par_civ3.Grid(:,1)=reshape(GridX,[],1);
-            par_civ3.Grid(:,2)=reshape(GridY,[],1);
-           
-            
+            par_civ3.Grid(:,2)=reshape(GridY,[],1);        
         end
         Shiftx=zeros(size(par_civ3.Grid,1),1);% shift expected from civ2 data
@@ -696,15 +692,13 @@
             mask=imread(par_civ3.Mask);
         end
-        ibx2=ceil(par_civ3.CorrBoxSize(1)/2);
-        iby2=ceil(par_civ3.CorrBoxSize(2)/2);
-%         par_civ3.SearchBoxSize(1)=2*ibx2+9;% search ara +-4 pixels around the guess
-%         par_civ3.SearchBoxSize(2)=2*iby2+9;
+%         ibx2=ceil(par_civ3.CorrBoxSize(1)/2);
+%         iby2=ceil(par_civ3.CorrBoxSize(2)/2);
         par_civ3.SearchBoxShift=[Shiftx(nbval>=1)./nbval(nbval>=1) Shifty(nbval>=1)./nbval(nbval>=1)];
         par_civ3.Grid=[par_civ3.Grid(nbval>=1,1)-par_civ3.SearchBoxShift(:,1)/2 par_civ3.Grid(nbval>=1,2)-par_civ3.SearchBoxShift(:,2)/2];% grid taken at the extrapolated origin of the displacement vectors
         if par_civ3.CheckDeformation
-            par_civ3.DUDX=DUDX./nbval;
-            par_civ3.DUDY=DUDY./nbval;
-            par_civ3.DVDX=DVDX./nbval;
-            par_civ3.DVDY=DVDY./nbval;
+            par_civ3.DUDX=DUDX(nbval>=1)./nbval(nbval>=1);
+            par_civ3.DUDY=DUDY(nbval>=1)./nbval(nbval>=1);
+            par_civ3.DVDX=DVDX(nbval>=1)./nbval(nbval>=1);
+            par_civ3.DVDY=DVDY(nbval>=1)./nbval(nbval>=1);
         end
         % calculate velocity data (y and v in indices, reverse to y component)
@@ -724,6 +718,6 @@
         
         nbvar=numel(Data.ListVarName);
-        Data.ListVarName=[Data.ListVarName {'Civ3_X','Civ3_Y','Civ3_U','Civ3_V','Civ3_F','Civ3_C','Xphys','Yphys','Zphys','Civ3_E'}];%  cell array containing the names of the fields to record
-        Data.VarDimName=[Data.VarDimName {'nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3'}];
+        Data.ListVarName=[Data.ListVarName {'Civ3_X','Civ3_Y','Civ3_U','Civ3_V','Civ3_F','Civ3_C','Xphys','Yphys','Zphys'}];%  cell array containing the names of the fields to record
+        Data.VarDimName=[Data.VarDimName {'nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3'}];
         Data.VarAttribute{nbvar+1}.Role='coord_x';
         Data.VarAttribute{nbvar+2}.Role='coord_y';
@@ -776,7 +770,7 @@
         Data.Patch3_SubDomainSize=Param.ActionInput.Patch3.SubDomainSize;
         nbvar=length(Data.ListVarName);
-        Data.ListVarName=[Data.ListVarName {'Civ3_U_smooth','Civ3_V_smooth','Civ3_SubRange','Civ3_NbCentres','Civ3_Coord_tps','Civ3_U_tps','Civ3_V_tps','Xmid','Ymid','Uphys','Vphys'}];
+        Data.ListVarName=[Data.ListVarName {'Civ3_U_smooth','Civ3_V_smooth','Civ3_SubRange','Civ3_NbCentres','Civ3_Coord_tps','Civ3_U_tps','Civ3_V_tps','Xmid','Ymid','Uphys','Vphys','Error'}];
         Data.VarDimName=[Data.VarDimName {'nb_vec_3','nb_vec_3',{'nb_coord','nb_bounds','nb_subdomain_3'},{'nb_subdomain_3'},...
-            {'nb_tps_3','nb_coord','nb_subdomain_3'},{'nb_tps_3','nb_subdomain_3'},{'nb_tps_3','nb_subdomain_3'},'nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3'}];
+            {'nb_tps_3','nb_coord','nb_subdomain_3'},{'nb_tps_3','nb_subdomain_3'},{'nb_tps_3','nb_subdomain_3'},'nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3'}];
         
         Data.VarAttribute{nbvar+1}.Role='vector_x';
@@ -805,6 +799,6 @@
         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(2)-Rangy(1))/(Npy-1);
-        [Data.Zphys,Data.Xphys,Data.Yphys,Data.Civ3_E]=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.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)
         if ~isempty(errormsg)
             disp_uvmat('ERROR',errormsg,checkrun)
@@ -813,7 +807,5 @@
         
     end
-    
-    
-    
+      
     
     %% write result in a netcdf file if requested
@@ -829,6 +821,6 @@
     else
        % store only phys data
-        Data_light.ListVarName={'Xphys','Yphys','Zphys','Civ3_C','Xmid','Ymid','Uphys','Vphys'};
-        Data_light.VarDimName={'nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3'};
+        Data_light.ListVarName={'Xphys','Yphys','Zphys','Civ3_C','Xmid','Ymid','Uphys','Vphys','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','nb_vec_3','nb_vec_3'};
         temp=find(Data.Civ3_FF==0);
         Data_light.Zphys=Data.Zphys(temp);
@@ -840,5 +832,5 @@
         Data_light.Uphys=Data.Uphys(temp);
         Data_light.Vphys=Data.Vphys(temp);
-       
+        Data_light.Error=Data.Error(temp);
        if exist('ncfile2','var')
             errormsg=struct2nc(ncfile2,Data_light);
@@ -853,5 +845,5 @@
    end
 
-  toc 
+  disp(['ellapsed time for the loop ' num2str(toc) ' s'])
 
 
@@ -912,8 +904,8 @@
 
 %% prepare correlation and search boxes 
-ibx2=ceil(par_civ.CorrBoxSize(1)/2);
-iby2=ceil(par_civ.CorrBoxSize(2)/2);
-isx2=ceil(par_civ.SearchBoxSize(1)/2);
-isy2=ceil(par_civ.SearchBoxSize(2)/2);
+ibx2=floor(par_civ.CorrBoxSize(1)/2);
+iby2=floor(par_civ.CorrBoxSize(2)/2);
+isx2=floor(par_civ.SearchBoxSize(1)/2);
+isy2=floor(par_civ.SearchBoxSize(2)/2);
 shiftx=round(par_civ.SearchBoxShift(:,1));
 shifty=-round(par_civ.SearchBoxShift(:,2));% sign minus because image j index increases when y decreases
@@ -944,23 +936,4 @@
 check_MaxIma=isfield(par_civ,'MaxIma') && ~isempty(par_civ.MaxIma);
 
-% %% prepare images
-% if isfield(par_civ,'reverse_pair')
-%     if par_civ.reverse_pair
-%         if ischar(par_civ.ImageB)
-%             temp=par_civ.ImageA;
-%             par_civ.ImageA=imread(par_civ.ImageB);
-%         end
-%         if ischar(temp)
-%             par_civ.ImageB=imread(temp);
-%         end
-%     end
-% else
-%     if ischar(par_civ.ImageA)
-%         par_civ.ImageA=imread(par_civ.ImageA);
-%     end
-%     if ischar(par_civ.ImageB)
-%         par_civ.ImageB=imread(par_civ.ImageB);
-%     end
-% end
 par_civ.ImageA=sum(double(par_civ.ImageA),3);%sum over rgb component for color images
 par_civ.ImageB=sum(double(par_civ.ImageB),3);
@@ -989,6 +962,6 @@
   %  check_noflux=(par_civ.Mask<100) ;%TODO: to implement
     check_undefined=(par_civ.Mask<200 & par_civ.Mask>=20 );
-    par_civ.ImageA(check_undefined)=MinA;% put image A to zero (i.e. the min image value) in the undefined  area
-    par_civ.ImageB(check_undefined)=MinB;% put image B to zero (i.e. the min image value) in the undefined  area
+%     par_civ.ImageA(check_undefined)=MinA;% put image A to zero (i.e. the min image value) in the undefined  area
+%     par_civ.ImageB(check_undefined)=MinB;% put image B to zero (i.e. the min image value) in the undefined  area
 end
 
@@ -997,8 +970,7 @@
 sum_square=1;% default
 mesh=1;% default
-CheckDecimal=isfield(par_civ,'CheckDecimal')&& par_civ.CheckDecimal==1;
-if CheckDecimal
-    mesh=0.2;%mesh in pixels for subpixel image interpolation
-    CheckDeformation=isfield(par_civ,'CheckDeformation')&& par_civ.CheckDeformation==1;
+CheckDeformation=isfield(par_civ,'CheckDeformation')&& par_civ.CheckDeformation==1;
+if CheckDeformation
+    mesh=0.25;%mesh in pixels for subpixel image interpolation
 end
 % vector=[0 0];%default
@@ -1020,13 +992,24 @@
     image1_crop=MinA*ones(numel(subrange1_y),numel(subrange1_x));% default value=min of image A
     image2_crop=MinA*ones(numel(subrange2_y),numel(subrange2_x));% default value=min of image A
+    mask1_crop=ones(numel(subrange1_y),numel(subrange1_x));% default value=1 for mask
+    mask2_crop=ones(numel(subrange2_y),numel(subrange2_x));% default value=1 for mask
     check1_x=subrange1_x>=1 & subrange1_x<=par_civ.ImageWidth;% check which points in the subimage 1 are contained in the initial image 1
     check1_y=subrange1_y>=1 & subrange1_y<=par_civ.ImageHeight;
     check2_x=subrange2_x>=1 & subrange2_x<=par_civ.ImageWidth;% check which points in the subimage 2 are contained in the initial image 2
-    check2_y=subrange2_y>=1 & subrange2_y<=par_civ.ImageHeight;
-    
+    check2_y=subrange2_y>=1 & subrange2_y<=par_civ.ImageHeight;   
     image1_crop(check1_y,check1_x)=par_civ.ImageA(subrange1_y(check1_y),subrange1_x(check1_x));%extract a subimage (correlation box) from image A
     image2_crop(check2_y,check2_x)=par_civ.ImageB(subrange2_y(check2_y),subrange2_x(check2_x));%extract a larger subimage (search box) from image B
-    image1_mean=mean(mean(image1_crop));
-    image2_mean=mean(mean(image2_crop));
+    mask1_crop(check1_y,check1_x)=check_undefined(subrange1_y(check1_y),subrange1_x(check1_x));%extract a mask subimage (correlation box) from image A
+    mask2_crop(check2_y,check2_x)=check_undefined(subrange2_y(check2_y),subrange2_x(check2_x));%extract a mask subimage (search box) from imag
+    sizemask=sum(sum(mask1_crop))/(numel(subrange1_y)*numel(subrange1_x));%size of the masked part relative to the correlation sub-image
+        if sizemask > 1/2% eliminate point if more than half of the correlation box is masked
+            F(ivec)=3; %
+        else
+            image1_crop=image1_crop.*~mask1_crop;% put to zero the masked pixels (mask1_crop='true'=1)
+            image2_crop=image2_crop.*~mask2_crop;
+            image1_mean=mean(mean(image1_crop))/(1-sizemask);
+            image2_mean=mean(mean(image2_crop))/(1-sizemask);
+%     image1_mean=mean(mean(image1_crop));
+%     image2_mean=mean(mean(image2_crop));
     %threshold on image minimum
     if check_MinIma && (image1_mean < par_civ.MinIma || image2_mean < par_civ.MinIma)
@@ -1037,22 +1020,20 @@
         F(ivec)=3;
     end
-    %         end
+    end
     if F(ivec)~=3
-        image1_crop=image1_crop-image1_mean;%substract the mean
-        image2_crop=image2_crop-image2_mean;
-        if CheckDecimal
-            xi=(1:mesh:size(image1_crop,2));
-            yi=(1:mesh:size(image1_crop,1))';
-            if CheckDeformation
+         image1_crop=(image1_crop-image1_mean).*~mask1_crop;%substract the mean, put to zero the masked parts
+            image2_crop=(image2_crop-image2_mean).*~mask2_crop;
+        if CheckDeformation
+              xi=(1:mesh:size(image1_crop,2));
+                yi=(1:mesh:size(image1_crop,1))';
                 [XI,YI]=meshgrid(xi-ceil(size(image1_crop,2)/2),yi-ceil(size(image1_crop,1)/2));
                 XIant=XI-par_civ.DUDX(ivec)*XI-par_civ.DUDY(ivec)*YI+ceil(size(image1_crop,2)/2);
                 YIant=YI-par_civ.DVDX(ivec)*XI-par_civ.DVDY(ivec)*YI+ceil(size(image1_crop,1)/2);
                 image1_crop=interp2(image1_crop,XIant,YIant);
-            else
-                image1_crop=interp2(image1_crop,xi,yi);
-            end
-            xi=(1:mesh:size(image2_crop,2));
-            yi=(1:mesh:size(image2_crop,1))';
-            image2_crop=interp2(image2_crop,xi,yi);
+                image1_crop(isnan(image1_crop))=0;
+                xi=(1:mesh:size(image2_crop,2));
+                yi=(1:mesh:size(image2_crop,1))';
+                image2_crop=interp2(image2_crop,xi,yi,'*spline');
+                image2_crop(isnan(image2_crop))=0;
         end
         sum_square=(sum(sum(image1_crop.*image1_crop)));%+sum(sum(image2_crop.*image2_crop)))/2;
@@ -1063,4 +1044,7 @@
         %Find the correlation max, at 255
         [y,x] = find(result_conv==255,1);
+        subimage2_crop=image2_crop(y:y+2*iby2/mesh,x:x+2*ibx2/mesh);%subimage of image 2 corresponding to the optimum displacement of first image
+            sum_square=sum_square*sum(sum(subimage2_crop.*subimage2_crop));% product of variances of image 1 and 2
+            sum_square=sqrt(sum_square);% sqrt of the variance product to normalise correlation
         if ~isempty(y) && ~isempty(x)
             try
@@ -1286,5 +1270,5 @@
 % 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)
+function [z,Xphy,Yphy,Error]=shift2z(xmid, ymid, u, v,XmlData)
 z=0;
 error=0;
@@ -1295,5 +1279,5 @@
 R=(Calib_A.R)';
 x_a=xmid- u/2;
-y_a=ymid- u/2;
+y_a=ymid- v/2; 
 z_a=R(7)*x_a+R(8)*y_a+Calib_A.Tx_Ty_Tz(1,3);
 Xa=(R(1)*x_a+R(2)*y_a+Calib_A.Tx_Ty_Tz(1,1))./z_a;
@@ -1311,6 +1295,10 @@
 
 %% second image
+%loading shift angle
+
 Calib_B=XmlData{2}.GeometryCalib;
 R=(Calib_B.R)';
+
+
 x_b=xmid+ u/2;
 y_b=ymid+ v/2;
@@ -1330,5 +1318,11 @@
 %% result
 Den=(Dxb-Dxa).*(Dxb-Dxa)+(Dyb-Dya).*(Dyb-Dya);
-error=((Dyb-Dya).*u-(Dxb-Dxa).*v)./Den;
+mfx=(XmlData{1}.GeometryCalib.fx_fy(1)+XmlData{2}.GeometryCalib.fx_fy(1))/2;
+mfy=(XmlData{1}.GeometryCalib.fx_fy(2)+XmlData{2}.GeometryCalib.fx_fy(2))/2;
+mtz=(XmlData{1}.GeometryCalib.Tx_Ty_Tz(1,3)+XmlData{2}.GeometryCalib.Tx_Ty_Tz(1,3))/2;
+
+Error=(sqrt(mfx^2+mfy^2)/(2*sqrt(2)*mtz)).*(((Dyb-Dya).*(-u)-(Dxb-Dxa).*(-v))./Den);
+
+z=((Dxb-Dxa).*(-u)+(Dyb-Dya).*(-v))./Den;
 
 xnew(1,:)=Dxa.*z+x_a;
@@ -1336,10 +1330,9 @@
 ynew(1,:)=Dya.*z+y_a;
 ynew(2,:)=Dyb.*z+y_b;
-
-Xphy=mean(xnew,1); %on moyenne les 2 valeurs 
+Xphy=mean(xnew,1);
 Yphy=mean(ynew,1);
-z=((Dxb-Dxa).*u-(Dyb-Dya).*v)./Den;
-
-
-
-
+
+
+
+
+
