Changeset 864

Timestamp:

Feb 6, 2015, 7:55:10 PM (10 years ago)

Author:

sommeria

Message:

bug_repaired

Location:

trunk/src

Files:

• calc_field_interp.m (modified) (1 diff)
• plot_field.m (modified) (1 diff)
• proj_field.m (modified) (1 diff)
• series/civ2vel_3C.m (modified) (4 diffs)
• series/stereo_civ.m (modified) (21 diffs)
• set_object.m (modified) (3 diffs)
• sub_field.m (modified) (2 diffs)
• uvmat.m (modified) (1 diff)

trunk/src/calc_field_interp.m

@@ -55 +55 @@
             check_skipped(ilist)=1; %variable not found
         elseif isempty(find(strcmp(FieldName{ilist},InputVarList), 1));% the variable exists and has not been already selected
-            if isfield(Data.VarAttribute{ivar},'Role') &&...
+            if exist('XI','var')&& isfield(Data.VarAttribute{ivar},'Role') &&...
                     (strcmp(Data.VarAttribute{ivar}.Role,'ancillary')||strcmp(Data.VarAttribute{ivar}.Role,'warnflag')||strcmp(Data.VarAttribute{ivar}.Role,'errorflag'))
                 check_interp(ilist)=0; % ancillary variable, not interpolated ?????

trunk/src/plot_field.m

@@ -640 +640 @@
             return
         else
+            if numel(CellInfo{icell}.VarIndex_vector_x)>1
+                errormsg='error in plot_field: attempt to plot two vector fields';
+            return
+            end
             test_vec=1;
             if isfield(CellInfo{icell},'VarIndex_errorflag')

trunk/src/proj_field.m

@@ -961 +961 @@
 if isfield(FieldData,'PlaneCoord')&&length(FieldData.PlaneCoord)==3&& isfield(ProjData,'ProjObjectCoord')
     if length(ProjData.ProjObjectCoord)==3% if the projection plane has a z coordinate
-        if ~isequal(ProjData.PlaneCoord(3),ProjData.ProjObjectCoord) %check the consistency with the z coordinate of the field plane (set by calibration)
+        if isfield(ProjData,'.PlaneCoord') && ~isequal(ProjData.PlaneCoord(3),ProjData.ProjObjectCoord) %check the consistency with the z coordinate of the field plane (set by calibration)
             errormsg='inconsistent z position for field and projection plane';
             return

trunk/src/series/civ2vel_3C.m

@@ -169 +169 @@
 %% MAIN LOOP ON FIELDS
 warning off
-for index=1:NbField
+if NbField<2 && length(filecell{:,1})>2
+    disp_uvmat('ERROR','you need at least 2 images to compute the mean position for the stereo.',checkrun)
+return
+end
+for index=1:NbField-1
     update_waitbar(WaitbarHandle,index/NbField)
     if ~isempty(RUNHandle) && ~strcmp(get(RUNHandle,'BusyAction'),'queue')
@@ -179 +183 @@
     Data=cell(1,NbView);%initiate the set Data
     timeread=zeros(1,NbView);
-    [Data{3},tild,errormsg] = nc2struct(filecell{3,index});
-    ZI=griddata(Data{3}.Xphys,Data{3}.Yphys,Data{3}.Zphys,XI,YI);
+    
+    %get Xphys,Yphys,Zphys from 1 or 2 stereo folders. Positions are taken
+    %at the middle between to time step
+   clear ZItemp
+   ZItemp=zeros(size(XI,1),size(XI,2),2);
+ for indextemp=index:index+1; 
+    if NbView==3 % if there is only 1 stereo folder, extract directly Xphys,Yphys and Zphys
+        
+        [Data{3},tild,errormsg] = nc2struct(filecell{3,indextemp}); 
+       
+        if  exist('Data{3}.Civ3_FF','var') % FF is present, remove wrong vector
+            temp=find(Data{3}.Civ3_FF==0);
+            Zphys=Data{3}.Zphys(temp);
+            Yphys=Data{3}.Yphys(temp);
+            Xphys=Data{3}.Xphys(temp);
+        else 
+            Zphys=Data{3}.Zphys;
+            Yphys=Data{3}.Yphys;
+            Xphys=Data{3}.Xphys;
+        end
+        
+        
+        
+    elseif NbView==4 % is there is 2 stereo folders, get global U and V and compute Zphys
+        
+        
+        %test if the seconde camera is the same for both folder
+        for i=3:4
+        indpt(i)=strfind(Param.InputTable{i,2},'.'); % indice of the "." is the folder name 1
+        indline(i)=strfind(Param.InputTable{i,2},'-'); % indice of the "-" is the folder name1
+        camname{i}=Param.InputTable{i,2}(indline(i)+1:indpt(i)-1);% extract the second camera name 
+        end
+        
+        if strcmp(camname{3},camname{4})==0 
+            disp_uvmat('ERROR','The 2 stereo folders should have the same camera for the second position',checkrun)
+            return
+        end
+        
+        [Data{3},tild,errormsg] = nc2struct(filecell{3,indextemp}); 
+        if exist('Data{3}.Civ3_FF','var') % if FF is present, remove wrong vector
+            temp=find(Data{3}.Civ3_FF==0);
+            Xmid3=Data{3}.Xmid(temp);
+            Ymid3=Data{3}.Ymid(temp);
+            U3=Data{3}.Uphys(temp);
+            V3=Data{3}.Vphys(temp);
+        else 
+            Xmid3=Data{3}.Xmid;
+            Ymid3=Data{3}.Ymid;
+            U3=Data{3}.Uphys;
+            V3=Data{3}.Vphys;
+        end
+        %temporary gridd of merging the 2 stereos datas
+        [xq,yq] = meshgrid(min(Xmid3+(U3)/2):(max(Xmid3+(U3)/2)-min(Xmid3+(U3)/2))/128:max(Xmid3+(U3)/2),min(Ymid3+(V3)/2):(max(Ymid3+(V3)/2)-min(Ymid3+(V3)/2))/128:max(Ymid3+(V3)/2));
+        
+        %1st folder : interpolate the first camera (Dalsa1) points on the second (common) camera
+        %(Dalsa 3)
+        x3Q=griddata(Xmid3+(U3)/2,Ymid3+(V3)/2,Xmid3-(U3)/2,xq,yq);
+        y3Q=griddata(Xmid3+(U3)/2,Ymid3+(V3)/2,Ymid3-(V3)/2,xq,yq);
+        
+        
+        [Data{4},tild,errormsg] = nc2struct(filecell{4,indextemp}); 
+        if exist('Data{4}.Civ3_FF','var') % if FF is present, remove wrong vector
+            temp=find(Data{4}.Civ3_FF==0);
+            Xmid4=Data{4}.Xmid(temp);
+            Ymid4=Data{4}.Ymid(temp);
+            U4=Data{4}.Uphys(temp);
+            V4=Data{4}.Vphys(temp);
+        else 
+            Xmid4=Data{4}.Xmid;
+            Ymid4=Data{4}.Ymid;
+            U4=Data{4}.Uphys;
+            V4=Data{4}.Vphys;
+        end
+        
+        %2nd folder :interpolate the first camera (Dalsa2) points on the second (common) camera
+        %(Dalsa 3)
+        x4Q=griddata(Xmid4+(U4)/2,Ymid4+(V4)/2,Xmid4-(U4)/2,xq,yq);
+        y4Q=griddata(Xmid4+(U4)/2,Ymid4+(V4)/2,Ymid4-(V4)/2,xq,yq);
+        
+        xmid=reshape((x4Q+x3Q)/2,length(xq(:,1)).*length(xq(1,:)),1);
+        ymid=reshape((y4Q+y3Q)/2,length(yq(:,1)).*length(yq(1,:)),1);
+        u=reshape(x4Q-x3Q,length(xq(:,1)).*length(xq(1,:)),1);
+        v=reshape(y4Q-y3Q,length(yq(:,1)).*length(yq(1,:)),1);
+        
+        
+        [Zphys,Xphys,Yphys,error]=shift2z(xmid, ymid, u, v,XmlData); %get Xphy,Yphy and Zphys
+        %remove NaN 
+        tempNaN=isnan(Zphys);tempind=find(tempNaN==1);
+        Zphys(tempind)=[];
+        Xphys(tempind)=[];
+        Yphys(tempind)=[];
+        error(tempind)=[];
+         
+    end
+    
+    
+   
+    
+       ZItemp(:,:,indextemp)=griddata(Xphys,Yphys,Zphys,XI,YI); %interpolation on the choosen gridd
+    
+end
+    ZI=mean(ZItemp,3); %mean between two the two time step
+    
     [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
@@ -208 +313 @@
         end
     end
-    Ua=griddata(Data{1}.X,Data{1}.Y,Data{1}.U,Xa,Ya);
-    Va=griddata(Data{1}.X,Data{1}.Y,Data{1}.V,Xa,Ya);
+    %remove wrong vector
+    temp=find(Data{1}.FF==0);
+    X1=Data{1}.X(temp);
+    Y1=Data{1}.Y(temp);
+    U1=Data{1}.U(temp);
+    V1=Data{1}.V(temp);
+    
+    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);
     
-    Ub=griddata(Data{2}.X,Data{2}.Y,Data{1}.U,Xb,Yb);
-    Vb=griddata(Data{2}.X,Data{2}.Y,Data{2}.V,Xb,Yb);
+    
+    temp=find(Data{2}.FF==0);
+    X2=Data{2}.X(temp);
+    Y2=Data{2}.Y(temp);
+    U2=Data{2}.U(temp);
+    V2=Data{2}.V(temp);
+    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);
     S=ones(size(XI,1),size(XI,2),3);
@@ -302 +420 @@
 
 
-
-
-
+function [z,Xphy,Yphy,error]=shift2z(xmid, ymid, u, v,XmlData)
+z=0;
+error=0;
+
+
+%% first image
+Calib_A=XmlData{1}.GeometryCalib;
+R=(Calib_A.R)';
+x_a=xmid- u/2;
+y_a=ymid- u/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;
+Ya=(R(4)*x_a+R(5)*y_a+Calib_A.Tx_Ty_Tz(1,2))./z_a;
+
+A_1_1=R(1)-R(7)*Xa;
+A_1_2=R(2)-R(8)*Xa;
+A_1_3=R(3)-R(9)*Xa;
+A_2_1=R(4)-R(7)*Ya;
+A_2_2=R(5)-R(8)*Ya;
+A_2_3=R(6)-R(9)*Ya;
+Det=A_1_1.*A_2_2-A_1_2.*A_2_1;
+Dxa=(A_1_2.*A_2_3-A_2_2.*A_1_3)./Det;
+Dya=(A_2_1.*A_1_3-A_1_1.*A_2_3)./Det;
+
+%% second image
+Calib_B=XmlData{2}.GeometryCalib;
+R=(Calib_B.R)';
+x_b=xmid+ u/2;
+y_b=ymid+ v/2;
+z_b=R(7)*x_b+R(8)*y_b+Calib_B.Tx_Ty_Tz(1,3);
+Xb=(R(1)*x_b+R(2)*y_b+Calib_B.Tx_Ty_Tz(1,1))./z_b;
+Yb=(R(4)*x_b+R(5)*y_b+Calib_B.Tx_Ty_Tz(1,2))./z_b;
+B_1_1=R(1)-R(7)*Xb;
+B_1_2=R(2)-R(8)*Xb;
+B_1_3=R(3)-R(9)*Xb;
+B_2_1=R(4)-R(7)*Yb;
+B_2_2=R(5)-R(8)*Yb;
+B_2_3=R(6)-R(9)*Yb;
+Det=B_1_1.*B_2_2-B_1_2.*B_2_1;
+Dxb=(B_1_2.*B_2_3-B_2_2.*B_1_3)./Det;
+Dyb=(B_2_1.*B_1_3-B_1_1.*B_2_3)./Det;
+
+%% result
+Den=(Dxb-Dxa).*(Dxb-Dxa)+(Dyb-Dya).*(Dyb-Dya);
+error=((Dyb-Dya).*u-(Dxb-Dxa).*v)./Den;
+
+xnew(1,:)=Dxa.*z+x_a;
+xnew(2,:)=Dxb.*z+x_b;
+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;
+
+

trunk/src/series/stereo_civ.m

@@ -41 +41 @@
 Data=[];
 errormsg='';
-
 %% set the input elements needed on the GUI series when the action is selected in the menu ActionName or InputTable refreshed
 if isstruct(Param) && isequal(Param.Action.RUN,0)% function activated from the GUI series but not RUN 
@@ -50 +49 @@
     path_series=fileparts(which('series'));
     addpath(fullfile(path_series,'series'))
-   % Data=civ_input(Param);% introduce the civ parameters using the GUI civ_input
-   Data=stereo_input(Param)
+    Data=stereo_input(Param);% introduce the civ parameters using the GUI civ_input
     if isempty(Data)
         Data=Param;% if  civ_input has been cancelled, keep previous parameters
@@ -263 +261 @@
 end
 
-
+tic
 %%%%% MAIN LOOP %%%%%%
 for ifield=1:NbField
@@ -271 +269 @@
         break
     end
+    % variable for light saving or not.
+       LSM=Param.ActionInput.CheckLSM;
+       
     Civ1Dir=OutputDir;
 
-        ncfile=fullfile_uvmat(RootPath_A,Civ1Dir,RootFile_A,'.nc',NomTypeNc,i2_series_Civ1(ifield),[],...
+        ncfile=fullfile_uvmat(RootPath_A,Civ1Dir,[RootFile_A,'_All'],'.nc',NomTypeNc,i2_series_Civ1(ifield),[],...
             j1_series_Civ1(ifield),j2_series_Civ1(ifield));
+        
+        
+        ncfile2=fullfile_uvmat(RootPath_A,Civ1Dir,[RootFile_A,'_Light'],'.nc',NomTypeNc,i2_series_Civ1(ifield),[],...
+            j1_series_Civ1(ifield),j2_series_Civ1(ifield));
+       
     %% Civ1
+   
     % if Civ1 computation is requested
     if isfield (Param.ActionInput,'Civ1')
@@ -294 +301 @@
         PhysImageA=fullfile_uvmat(RootPath_A,Civ1Dir,RootFile_A,'.png','_1a',i1_series_Civ1(ifield),[],1);
         PhysImageB=fullfile_uvmat(RootPath_A,Civ1Dir,RootFile_A,'.png','_1a',i1_series_Civ1(ifield),[],2);
+        if LSM ~= 1
         imwrite(A{1},PhysImageA)
         imwrite(A{2},PhysImageB)
+        end
+        
         par_civ1.ImageA=A{1};
         par_civ1.ImageB=A{2};
@@ -338 +348 @@
         
         % calculate velocity data (y and v in indices, reverse to y component)
-        [xtable ytable utable vtable ctable F result_conv errormsg] = civ (par_civ1);
+        [xtable, ytable, utable, vtable, ctable, F, result_conv, errormsg] = civ (par_civ1);
         Data.Civ1_X=reshape(xtable,[],1);
         Data.Civ1_Y=reshape(par_civ1.ImageHeight-ytable+1,[],1);
@@ -346 +356 @@
         Data.Civ1_C=reshape(ctable,[],1);
         Data.Civ1_F=reshape(F,[],1);
-%         Data.Xphys=Rangx(1)+(Rangx(2)-Rangx(1))*(Data.Civ1_X-0.5)/(Npx-1);
-%         Data.Yphys=Rangy(1)+(Rangy(2)-Rangy(1))*(Data.Civ1_Y-0.5)/(Npy-1);
-%         
-%         U=Data.Civ1_U*(Rangx(2)-Rangx(1))/(Npx-1);
-%         V=Data.Civ1_V*(Rangy(2)-Rangy(1))/(Npy-1);
-%         [Data.Zphys,Data.Civ1_E]=shift2z(Data.Xphys,Data.Yphys,U,V,XmlData);
-%         if ~isempty(errormsg)
-%             disp_uvmat('ERROR',errormsg,checkrun)
-%             return
-%         end
+
     end
     
@@ -425 +426 @@
         Data.Civ1_FF(ind_good)=FFres;
         Data.CivStage=3;
+             
+        
+        
+        
+%               
+%          % get z from u and v (displacements)
+%        
+%         tempXmid=Rangx(1)+(Rangx(2)-Rangx(1))*(Data.Civ1_X-0.5)/(Npx-1);%temporary coordinate (velocity taken at the point middle from imgae 1 and 2)
+%         tempYmid=Rangy(1)+(Rangy(2)-Rangy(1))*(Data.Civ1_Y-0.5)/(Npy-1);%temporary coordinate (velocity taken at the point middle from imgae 1 and 2)
+%         tempUphys=Data.Civ1_U_smooth*(Rangx(2)-Rangx(1))/(Npx-1);
+%         tempVphys=Data.Civ1_V_smooth*(Rangy(2)-Rangy(1))/(Npy-1);
+%         [tempZphys,tempXphys,tempYphys,tempCiv3_E]=shift2z(tempXmid,tempYmid,tempUphys,tempVphys,XmlData); %Data.Xphys and Data.Xphys are real coordinate (geometric correction more accurate than xtemp/ytemp)
+%         temp=find(Data.Civ1_FF~=0);
+%         tempXmid(temp)=[];
+%         tempYmid(temp)=[];
+%         tempZphys(temp)=[];
+%             
     end
     
@@ -446 +464 @@
             j2=j2_series_Civ2(ifield);
         end
-        par_civ2.ImageWidth=size(par_civ2.ImageA,2);%FileInfo_A.Width;
-        par_civ2.ImageHeight=size(par_civ2.ImageA,1);%FileInfo_A.Height;
+        par_civ2.ImageWidth=size(par_civ2.ImageA,2);
+        par_civ2.ImageHeight=size(par_civ2.ImageA,1);
+        
         if isfield(par_civ2,'Grid')% grid points set as input file
             if ischar(par_civ2.Grid)%read the grid file if the input is a file name
@@ -461 +480 @@
             par_civ2.Grid(:,1)=reshape(GridX,[],1);
             par_civ2.Grid(:,2)=reshape(GridY,[],1);
+           
+            
         end
         Shiftx=zeros(size(par_civ2.Grid,1),1);% shift expected from civ1 data
@@ -508 +529 @@
         end
         % calculate velocity data (y and v in indices, reverse to y component)
-        [xtable ytable utable vtable ctable F] = civ (par_civ2);
+        [xtable, ytable, utable, vtable, ctable, F] = civ (par_civ2);
         list_param=(fieldnames(Param.ActionInput.Civ2))';
         Civ2_param=regexprep(list_param,'^.+','Civ2_$0');% insert 'Civ2_' before  each string in list_param
@@ -523 +544 @@
         
         nbvar=numel(Data.ListVarName);
-        Data.ListVarName=[Data.ListVarName {'Civ2_X','Civ2_Y','Civ2_U','Civ2_V','Civ2_F','Civ2_C','Xphys','Yphys','Zphys','Civ2_E'}];%  cell array containing the names of the fields to record
-        Data.VarDimName=[Data.VarDimName {'nb_vec_2','nb_vec_2','nb_vec_2','nb_vec_2','nb_vec_2','nb_vec_2','nb_vec_2','nb_vec_2','nb_vec_2','nb_vec_2'}];
+        Data.ListVarName=[Data.ListVarName {'Civ2_X','Civ2_Y','Civ2_U','Civ2_V','Civ2_F','Civ2_C'}];%  cell array containing the names of the fields to record
+        Data.VarDimName=[Data.VarDimName {'nb_vec_2','nb_vec_2','nb_vec_2','nb_vec_2','nb_vec_2','nb_vec_2'}];
         Data.VarAttribute{nbvar+1}.Role='coord_x';
         Data.VarAttribute{nbvar+2}.Role='coord_y';
@@ -598 +619 @@
         
             
+ 
+        
+    end
+    
+        
+    %% Civ3
+    
+    if isfield (Param.ActionInput,'Civ3')
+        par_civ3=Param.ActionInput.Civ3;
+        par_civ3.ImageA=par_civ1.ImageA;
+        par_civ3.ImageB=par_civ1.ImageB;
+        %         if ~isfield(Param.Civ1,'ImageA')
+%         i1=i1_series_Civ3(ifield);
+%         i2=i1;
+%         if ~isempty(i2_series_Civ3)
+%             i2=i2_series_Civ3(ifield);
+%         end
+%         j1=1;
+%         if ~isempty(j1_series_Civ3)
+%             j1=j1_series_Civ3(ifield);
+%         end
+%         j2=j1;
+%         if ~isempty(j2_series_Civ3)
+%             j2=j2_series_Civ3(ifield);
+%         end
+        par_civ3.ImageWidth=size(par_civ3.ImageA,2);
+        par_civ3.ImageHeight=size(par_civ3.ImageA,1);
+        
+        if isfield(par_civ3,'Grid')% grid points set as input file
+            if ischar(par_civ3.Grid)%read the grid file if the input is a file name
+                par_civ3.Grid=dlmread(par_civ3.Grid);
+                par_civ3.Grid(1,:)=[];%the first line must be removed (heading in the grid file)
+            end
+        else% automatic grid
+            minix=floor(par_civ3.Dx/2)-0.5;
+            maxix=minix+par_civ3.Dx*floor((par_civ3.ImageWidth-1)/par_civ3.Dx);
+            miniy=floor(par_civ3.Dy/2)-0.5;
+            maxiy=minix+par_civ3.Dy*floor((par_civ3.ImageHeight-1)/par_civ3.Dy);
+            [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);
+           
+            
+        end
+        Shiftx=zeros(size(par_civ3.Grid,1),1);% shift expected from civ2 data
+        Shifty=zeros(size(par_civ3.Grid,1),1);
+        nbval=zeros(size(par_civ3.Grid,1),1);
+        if par_civ3.CheckDeformation
+            DUDX=zeros(size(par_civ3.Grid,1),1);
+            DUDY=zeros(size(par_civ3.Grid,1),1);
+            DVDX=zeros(size(par_civ3.Grid,1),1);
+            DVDY=zeros(size(par_civ3.Grid,1),1);
+        end
+        NbSubDomain=size(Data.Civ2_SubRange,3);
+        % get the guess from patch2
+        for isub=1:NbSubDomain% for each sub-domain of Patch2
+            nbvec_sub=Data.Civ2_NbCentres(isub);% nbre of Civ2 vectors in the subdomain
+            ind_sel=find(par_civ3.Grid(:,1)>=Data.Civ2_SubRange(1,1,isub) & par_civ3.Grid(:,1)<=Data.Civ2_SubRange(1,2,isub) &...
+                par_civ3.Grid(:,2)>=Data.Civ2_SubRange(2,1,isub) & par_civ3.Grid(:,2)<=Data.Civ2_SubRange(2,2,isub));
+            epoints = par_civ3.Grid(ind_sel,:);% coordinates of interpolation sites
+            ctrs=Data.Civ2_Coord_tps(1:nbvec_sub,:,isub) ;%(=initial points) ctrs
+            nbval(ind_sel)=nbval(ind_sel)+1;% records the number of values for eacn interpolation point (in case of subdomain overlap)
+            EM = tps_eval(epoints,ctrs);
+            Shiftx(ind_sel)=Shiftx(ind_sel)+EM*Data.Civ2_U_tps(1:nbvec_sub+3,isub);
+            Shifty(ind_sel)=Shifty(ind_sel)+EM*Data.Civ2_V_tps(1:nbvec_sub+3,isub);
+            if par_civ3.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*Data.Civ2_U_tps(1:nbvec_sub+3,isub);
+                DUDY(ind_sel)=DUDY(ind_sel)+EMDY*Data.Civ2_U_tps(1:nbvec_sub+3,isub);
+                DVDX(ind_sel)=DVDX(ind_sel)+EMDX*Data.Civ2_V_tps(1:nbvec_sub+3,isub);
+                DVDY(ind_sel)=DVDY(ind_sel)+EMDY*Data.Civ2_V_tps(1:nbvec_sub+3,isub);
+            end
+        end
+        mask='';
+        if par_civ3.CheckMask&&~isempty(par_civ3.Mask)&& ~strcmp(maskname,par_civ3.Mask)% mask exist, not already read in Civ2
+            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;
+        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;
+        end
+        % calculate velocity data (y and v in indices, reverse to y component)
+        [xtable, ytable, utable, vtable, ctable, F] = civ (par_civ3);
+        list_param=(fieldnames(Param.ActionInput.Civ3))';
+        Civ3_param=regexprep(list_param,'^.+','Civ3_$0');% insert 'Civ3_' before  each string in list_param
+        Civ3_param=[{'Civ3_ImageA','Civ3_ImageB','Civ3_Time','Civ3_Dt'} Civ3_param]; %insert the names of the two input images
+        %indicate the values of all the global attributes in the output data
+        Data.Civ3_ImageA=ImageName_A;
+        Data.Civ3_ImageB=ImageName_B;
+        Data.Civ3_Time=(time(i2+1,j2+1)+time(i1+1,j1+1))/2;
+        Data.Civ3_Dt=0;
+        for ilist=1:length(list_param)
+            Data.(Civ3_param{4+ilist})=Param.ActionInput.Civ3.(list_param{ilist});
+        end
+        Data.ListGlobalAttribute=[Data.ListGlobalAttribute Civ3_param];
+        
+        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.VarAttribute{nbvar+1}.Role='coord_x';
+        Data.VarAttribute{nbvar+2}.Role='coord_y';
+        Data.VarAttribute{nbvar+3}.Role='vector_x';
+        Data.VarAttribute{nbvar+4}.Role='vector_y';
+        Data.VarAttribute{nbvar+5}.Role='warnflag';
+        Data.Civ3_X=reshape(xtable,[],1);
+        Data.Civ3_Y=reshape(size(par_civ3.ImageA,1)-ytable+1,[],1);
+        Data.Civ3_U=reshape(utable,[],1);
+        Data.Civ3_V=reshape(-vtable,[],1);
+        Data.Civ3_C=reshape(ctable,[],1);
+        Data.Civ3_F=reshape(F,[],1);
+        Data.CivStage=Data.CivStage+1;
+end
+    
+    %% Fix3
+    if isfield (Param.ActionInput,'Fix3')
+        ListFixParam=fieldnames(Param.ActionInput.Fix3);
+        for ilist=1:length(ListFixParam)
+            ParamName=ListFixParam{ilist};
+            ListName=['Fix3_' ParamName];
+            eval(['Data.ListGlobalAttribute=[Data.ListGlobalAttribute ''' ParamName '''];'])
+            eval(['Data.' ListName '=Param.ActionInput.Fix3.' ParamName ';'])
+        end
+        if check_civx
+            if ~isfield(Data,'fix3')
+                Data.ListGlobalAttribute=[Data.ListGlobalAttribute 'fix3'];
+                Data.fix3=1;
+                Data.ListVarName=[Data.ListVarName {'vec3_FixFlag'}];
+                Data.VarDimName=[Data.VarDimName {'nb_vectors3'}];
+            end
+            Data.vec_FixFlag=fix(Param.Fix3,Data.vec3_F,Data.vec3_C,Data.vec3_U,Data.vec3_V,Data.vec3_X,Data.vec3_Y);
+        else
+            Data.ListVarName=[Data.ListVarName {'Civ3_FF'}];
+            Data.VarDimName=[Data.VarDimName {'nb_vec_3'}];
+            nbvar=length(Data.ListVarName);
+            Data.VarAttribute{nbvar}.Role='errorflag';
+            Data.Civ3_FF=double(fix(Param.ActionInput.Fix3,Data.Civ3_F,Data.Civ3_C,Data.Civ3_U,Data.Civ3_V));
+            Data.CivStage=Data.CivStage+1;
+        end
+        
+    end
+
+    
+     %% Patch3
+    if isfield (Param.ActionInput,'Patch3')
+        Data.ListGlobalAttribute=[Data.ListGlobalAttribute {'Patch3_Rho','Patch3_Threshold','Patch3_SubDomain'}];
+        Data.Patch3_FieldSmooth=Param.ActionInput.Patch3.FieldSmooth;
+        Data.Patch3_MaxDiff=Param.ActionInput.Patch3.MaxDiff;
+        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.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'}];
+        
+        Data.VarAttribute{nbvar+1}.Role='vector_x';
+        Data.VarAttribute{nbvar+2}.Role='vector_y';
+        Data.VarAttribute{nbvar+5}.Role='coord_tps';
+        Data.VarAttribute{nbvar+6}.Role='vector_x';
+        Data.VarAttribute{nbvar+7}.Role='vector_y';
+        Data.Civ3_U_smooth=zeros(size(Data.Civ3_X));
+        Data.Civ3_V_smooth=zeros(size(Data.Civ3_X));
+        if isfield(Data,'Civ3_FF')
+            ind_good=find(Data.Civ3_FF==0);
+        else
+            ind_good=1:numel(Data.Civ3_X);
+        end
+        [Data.Civ3_SubRange,Data.Civ3_NbCentres,Data.Civ3_Coord_tps,Data.Civ3_U_tps,Data.Civ3_V_tps,tild,Ures, Vres,tild,FFres]=...
+            filter_tps([Data.Civ3_X(ind_good) Data.Civ3_Y(ind_good)],Data.Civ3_U(ind_good),Data.Civ3_V(ind_good),[],Data.Patch3_SubDomainSize,Data.Patch3_FieldSmooth,Data.Patch3_MaxDiff);
+        Data.Civ3_U_smooth(ind_good)=Ures;
+        Data.Civ3_V_smooth(ind_good)=Vres;
+        Data.Civ3_FF(ind_good)=FFres;
+        Data.CivStage=Data.CivStage+1;
+        
+            
          % get z from u and v (displacements)
        
-        Data.Xphys=Rangx(1)+(Rangx(2)-Rangx(1))*(Data.Civ2_X-0.5)/(Npx-1);
-        Data.Yphys=Rangy(1)+(Rangy(2)-Rangy(1))*(Data.Civ2_Y-0.5)/(Npy-1);
-        U=Data.Civ2_U_smooth*(Rangx(2)-Rangx(1))/(Npx-1);
-        V=Data.Civ2_V_smooth*(Rangy(2)-Rangy(1))/(Npy-1);
-        [Data.Zphys,Data.Civ2_E]=shift2z(Data.Xphys,Data.Yphys,U,V,XmlData);
+        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(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)
         if ~isempty(errormsg)
             disp_uvmat('ERROR',errormsg,checkrun)
@@ -612 +814 @@
     end
     
+    
+    
+    
     %% write result in a netcdf file if requested
-    if exist('ncfile','var')
-        errormsg=struct2nc(ncfile,Data);
-        if isempty(errormsg)
-            disp([ncfile ' written'])
-        else
-            disp(errormsg)
-        end
-    end
-end
+    if LSM ~= 1 % store all data
+        if exist('ncfile','var')
+            errormsg=struct2nc(ncfile,Data);
+            if isempty(errormsg)
+                disp([ncfile ' written'])
+            else
+                disp(errormsg)
+            end
+        end
+    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'};
+        temp=find(Data.Civ3_FF==0);
+        Data_light.Zphys=Data.Zphys(temp);
+        Data_light.Yphys=Data.Yphys(temp);
+        Data_light.Xphys=Data.Xphys(temp);
+        Data_light.Civ3_C=Data.Civ3_C(temp);
+        Data_light.Xmid=Data.Xmid(temp);
+        Data_light.Ymid=Data.Ymid(temp);
+        Data_light.Uphys=Data.Uphys(temp);
+        Data_light.Vphys=Data.Vphys(temp);
+       
+       if exist('ncfile2','var')
+            errormsg=struct2nc(ncfile2,Data_light);
+            if isempty(errormsg)
+                disp([ncfile2 ' written'])
+            else
+                disp(errormsg)
+            end
+       end
+       
+    end
+   end
+
+  toc 
+
+
 
 % 'civ': function piv.m adapted from PIVlab http://pivlab.blogspot.com/
@@ -769 +1003 @@
 end
 % vector=[0 0];%default
+
 for ivec=1:nbvec
     iref=round(par_civ.Grid(ivec,1)+0.5);% xindex on the image A for the middle of the correlation box
     jref=round(par_civ.ImageHeight-par_civ.Grid(ivec,2)+0.5);% yindex on the image B for the middle of the correlation box
+    
     %if ~(checkmask && par_civ.Mask(jref,iref)<=20) %velocity not set to zero by the black mask
     %         if jref-iby2<1 || jref+iby2>par_civ.ImageHeight|| iref-ibx2<1 || iref+ibx2>par_civ.ImageWidth||...
@@ -820 +1056 @@
             image2_crop=interp2(image2_crop,xi,yi);
         end
-        sum_square=(sum(sum(image1_crop.*image1_crop))+sum(sum(image2_crop.*image2_crop)))/2;
+        sum_square=(sum(sum(image1_crop.*image1_crop)));%+sum(sum(image2_crop.*image2_crop)))/2;
         %reference: Oliver Pust, PIV: Direct Cross-Correlation
         result_conv= conv2(image2_crop,flipdim(flipdim(image1_crop,2),1),'valid');
@@ -834 +1070 @@
                     [vector,F(ivec)] = SUBPIX2DGAUSS (result_conv,x,y);
                 end
-
-%                 if isfield(par_civ,'CheckDeformation')
-%                  utable(ivec)=shiftx(ivec);
-%                  vtable(ivec)=shifty(ivec);
-%                 else
-                 utable(ivec)=vector(1)*mesh+shiftx(ivec);
+                
+
+                utable(ivec)=vector(1)*mesh+shiftx(ivec);
                 vtable(ivec)=vector(2)*mesh+shifty(ivec);
-%                 end
-%                 xtable(ivec)=iref+utable(ivec)/2-0.5;% convec flow (velocity taken at the point middle from imgae 1 and 2)
-%                 ytable(ivec)=jref+vtable(ivec)/2-0.5;% and position of pixel 1=0.5 (convention for image coordinates=0 at the edge)
+
+                
+                xtable(ivec)=iref+utable(ivec)/2-0.5;% convec flow (velocity taken at the point middle from imgae 1 and 2)
+                ytable(ivec)=jref+vtable(ivec)/2-0.5;% and position of pixel 1=0.5 (convention for image coordinates=0 at the edge)
+                
                 iref=round(xtable(ivec));% image index for the middle of the vector
                 jref=round(ytable(ivec));
@@ -1051 +1286 @@
 % 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,error]=shift2z(xmid, ymid, u, v,XmlData)
+function [z,Xphy,Yphy,error]=shift2z(xmid, ymid, u, v,XmlData)
 z=0;
 error=0;
@@ -1059 +1294 @@
 Calib_A=XmlData{1}.GeometryCalib;
 R=(Calib_A.R)';
-x_a=xmid;
-y_a=ymid;
+x_a=xmid- u/2;
+y_a=ymid- u/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;
@@ -1078 +1313 @@
 Calib_B=XmlData{2}.GeometryCalib;
 R=(Calib_B.R)';
-x_b=xmid+ u;
-y_b=ymid+ v;
+x_b=xmid+ u/2;
+y_b=ymid+ v/2;
 z_b=R(7)*x_b+R(8)*y_b+Calib_B.Tx_Ty_Tz(1,3);
 Xb=(R(1)*x_b+R(2)*y_b+Calib_B.Tx_Ty_Tz(1,1))./z_b;
@@ -1096 +1331 @@
 Den=(Dxb-Dxa).*(Dxb-Dxa)+(Dyb-Dya).*(Dyb-Dya);
 error=((Dyb-Dya).*u-(Dxb-Dxa).*v)./Den;
+
+xnew(1,:)=Dxa.*z+x_a;
+xnew(2,:)=Dxb.*z+x_b;
+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;
 
+
+
+

trunk/src/set_object.m

@@ -326 +326 @@
         set(handles.num_RangeX_2,'TooltipString',['num_RangeX_2: half length of the ' ObjectStyle])
         set(handles.num_RangeY_2,'TooltipString',['num_RangeY_2: half width of the ' ObjectStyle])
-%         set(handles.XObject,'TooltipString',['XObject:  x coordinate of the ' Type ' centre'])
-%         set(handles.YObject,'TooltipString',['YObject:  y coordinate of the ' Type ' centre'])
     case {'plane'}  
         set(handles.num_Angle_3,'Visible','on')
@@ -334 +332 @@
         set(handles.num_RangeY_1,'Visible','on')
         set(handles.num_RangeY_2,'Visible','on')
-%         set(handles.XObject,'TooltipString',['XObject:  x coordinate of the axis origin for the ' Type])
-%         set(handles.YObject,'TooltipString',['YObject:  y coordinate of the axis origin for the ' Type])
         set(handles.num_RangeZ_2,'TooltipString','num_ZMax: range of projection normal to the plane')
         if test3D
@@ -374 +370 @@
         end
 end
+% set default values read in the plot of uvmat to initiate the mesh 
+if isequal(ProjMode,'interp_lin')|| isequal(ProjMode,'interp_tps')
+    if isempty(str2num(get(handles.num_DX,'String')))||isempty(str2num(get(handles.num_DY,'String')));     
+        huvmat=findobj('Tag','uvmat');%find the current uvmat interface handle
+        UvData=get(huvmat,'UserData');%Data associated to the current uvmat interface
+        Field=UvData.Field;
+        if  isfield(UvData.Field,'CoordMesh')&&~isempty(UvData.Field.CoordMesh)
+            set(handles.num_DX,'String',num2str(UvData.Field.CoordMesh))
+            set(handles.num_DY,'String',num2str(UvData.Field.CoordMesh))
+            set(handles.num_RangeX_1,'String',num2str(UvData.Field.XMin))
+            set(handles.num_RangeX_2,'String',num2str(UvData.Field.XMax))
+            set(handles.num_RangeY_1,'String',num2str(UvData.Field.YMin))
+            set(handles.num_RangeY_2,'String',num2str(UvData.Field.YMax))
+        end
+        if isempty(get(handles.CoordUnit,'String'))
+            set(handles.CoordUnit,'String',Field.CoordUnit)
+        end       
+    end
+end
+
 %------------------------------------------------------------------------
 

trunk/src/sub_field.m

@@ -129 +129 @@
 for ilist=1:numel(Field_1.ListVarName)
     % case of variable with a single dimension
-    OldDimName=Field_1.VarDimName{ilist};
-    if ischar(OldDimName), OldDimName={OldDimName}; end% transform char string to cell if relevant
-    if numel(OldDimName)==1
-        OldDim=Field_1.(Field_1.ListVarName{ilist});% get variable
-        %look for the existence of the variable OldDim in the first field Field
-        for i1=1:numel(Field.ListVarName)
-            if  isequal(Field.(Field.ListVarName{i1}),OldDim) &&...
-                   ((isempty(ProjModeRequest{i1}) && isempty(ProjModeRequest_1{ilist}))  || strcmp(ProjModeRequest{i1},ProjModeRequest_1{ilist}))                
-               ind_remove(ilist)=1;
-               NewDimName=Field.VarDimName{i1};
-               if ischar(NewDimName), NewDimName={NewDimName}; end %transform char chain to cell if needed
-               Field_1.VarDimName=regexprep_r(Field_1.VarDimName,['^' OldDimName{1} '$'],NewDimName{1});% change the var name of Field_1 to the corresponding var name of Field
+    if ~isempty(regexp(Field_1.VarAttribute{ilist}.Role,'^coord'))
+        OldDimName=Field_1.VarDimName{ilist};
+        if ischar(OldDimName), OldDimName={OldDimName}; end% transform char string to cell if relevant
+        if numel(OldDimName)==1
+            OldDim=Field_1.(Field_1.ListVarName{ilist});% get variable
+            %look for the existence of the variable OldDim in the first field Field
+            for i1=1:numel(Field.ListVarName)
+                if  isequal(Field.(Field.ListVarName{i1}),OldDim) &&...
+                        ((isempty(ProjModeRequest{i1}) && isempty(ProjModeRequest_1{ilist}))  || strcmp(ProjModeRequest{i1},ProjModeRequest_1{ilist}))
+                    ind_remove(ilist)=1;
+                    NewDimName=Field.VarDimName{i1};
+                    if ischar(NewDimName), NewDimName={NewDimName}; end %transform char chain to cell if needed
+                    Field_1.VarDimName=regexprep_r(Field_1.VarDimName,['^' OldDimName{1} '$'],NewDimName{1});% change the var name of Field_1 to the corresponding var name of Field
+                end
             end
         end
@@ -152 +154 @@
 ListVarNameNew=Field_1.ListVarName;
 check_rename=zeros(size(ListVarNameNew));
+check_remove=zeros(size(ListVarNameNew));
 for ilist=1:numel(ListVarNameNew)
     VarName=Field_1.ListVarName{ilist};
     ind_prev=find(strcmp(ListVarNameNew{ilist},Field.ListVarName),1);% look for duplicated variable name
     if ~isempty(ind_prev)% variable name exists in Field
-        check_rename(ilist)=1;
-        ListVarNameNew{ilist}=[ListVarNameNew{ilist} '_1'];
+%         check_rename(ilist)=0;
+        check_remove(ilist)=1;
+        if isfield(Field_1.VarAttribute{ilist},'Role')&& ismember(Field_1.VarAttribute{ilist}.Role,{'scalar','vector_x','vector_y'})
+                            ListVarNameNew{ilist}=[ListVarNameNew{ilist} '_1'];
+                            check_rename(ilist)=1;
+        end
     end
     SubData.ListVarName=[SubData.ListVarName ListVarNameNew{ilist}];
     SubData.VarDimName=[SubData.VarDimName Field_1.VarDimName(ilist)];
+    Field_1.VarAttribute{ilist}.CheckSub=1;
+    SubData.VarAttribute=[SubData.VarAttribute Field_1.VarAttribute{ilist}];
     SubData.(ListVarNameNew{ilist})=Field_1.(VarName);% teke the values of the old variable for the newly named one
     %SubData.VarAttribute=[SubData.VarAttribute Field_1.VarDimName(ilist)];

trunk/src/uvmat.m

@@ -5528 +5528 @@
         data.Type='plane';
     end
+    
+    %% initiate the new projection object
     hset_object=set_object(data,[],ZBounds);
     set(hset_object,'name','set_object')