Changeset 878 for trunk/src/series/civ2vel_3C.m

Timestamp:

Feb 23, 2015, 10:17:04 PM (9 years ago)

Author:

sommeria

Message:

stereo_civ and civ2vel_3C updated, spline introduced for civ_series:deformation

File:

• trunk/src/series/civ2vel_3C.m (modified) (13 diffs)

trunk/src/series/civ2vel_3C.m

@@ -169 +169 @@
 %% 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 +285 @@
 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 +321 @@
     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 +333 @@
     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 +356 @@
     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 +367 @@
     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 +406 @@
     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 +419 @@
 
 
-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 +431 @@
 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 +459 @@
 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 +475 @@
 
 %% 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 +498 @@
 %% 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 +510 @@
 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); 
+
+
+
+