Changeset 109 for trunk/src

Timestamp:

Oct 4, 2010, 10:39:11 PM (14 years ago)

Author:

sommeria

Message:

merge_proj.m, proj_field.m: bugs repaired merging of images
plot_field.m: minor cleaning
imadoc2struct: introduce the possibility of readying the calibration point coordinates
sub_field.m:bugs repaired
geometry_calib, create_grid: introduction of new calibration methods, improvement of detect_grid

Location:

trunk/src

Files:

• create_grid.m (modified) (1 diff)
• geometry_calib.fig (modified) (previous)
• geometry_calib.m (modified) (28 diffs)
• imadoc2struct.m (modified) (6 diffs)
• plot_field.m (modified) (4 diffs)
• proj_field.m (modified) (31 diffs)
• px_XYZ.m (modified) (1 diff)
• series/merge_proj.m (modified) (9 diffs)
• sub_field.m (modified) (5 diffs)

trunk/src/create_grid.m

@@ -136 +136 @@
     zarray=T.z_0*ones(size(yarray));
     varargout{1}=[xarray yarray zarray];
+    varargout{2}=T;
 end
-varargout{2}=T;
+
 % The figure can be deleted now
 delete(handles.figure1);

trunk/src/geometry_calib.m

@@ -3 +3 @@
 % function varargout = geometry_calib(varargin)
 %
-%A%AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
+%A%AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
 %  Copyright Joel Sommeria, 2008, LEGI / CNRS-UJF-INPG, sommeria@coriolis-legi.org.
 %AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
@@ -43 +43 @@
 % Edit the above text to modify the response to help geometry_calib
 
-% Last Modified by GUIDE v2.5 25-Mar-2010 19:10:05
+% Last Modified by GUIDE v2.5 03-Oct-2010 09:34:12
 
 % Begin initialization code - DO NOT edit
@@ -87 +87 @@
     set(hObject,'Position',pos_gui);
 end
+
+%set menu of calibration options
+%set(handles.calib_type,'String',{'rescale';'linear';'perspective';'normal';'tsai';'bouguet';'extrinsic'})
+set(handles.calib_type,'String',{'rescale';'linear';'quadr';'3D_linear';'3D_quadr';'3D_extrinsic'})
 inputxml='';
 if exist('inputfile','var')& ~isempty(inputfile)
-%     [Path,Name,ext]=fileparts(inputfile);
-%     form=imformats(ext([2:end]));
-%     if ~isempty(form)% if the input file is an image
-        struct.XmlInputfile=inputfile;
-        set(hObject,'UserData',struct)
-        [Pathsub,RootFile,field_count,str2,str_a,str_b,ext,nom_type,subdir]=name2display(inputfile);
-        inputxml=[fullfile(Pathsub,RootFile) '.xml'];
-%     end   
+    struct.XmlInputfile=inputfile;
+    set(hObject,'UserData',struct)
+    [Pathsub,RootFile,field_count,str2,str_a,str_b,ext,nom_type,subdir]=name2display(inputfile);
+    inputxml=[fullfile(Pathsub,RootFile) '.xml'];
 end
 set(handles.ListCoord,'String',{'......'})
@@ -102 +102 @@
     loadfile(handles,inputxml)% load the point coordiantes existing in the xml file
 end
-
 set(handles.ListCoord,'KeyPressFcn',{@key_press_fcn,handles})%set keyboard action function
 
@@ -113 +112 @@
 varargout{1} = handles.output;
 varargout{2}=handles;
-
-%------------
-function Phi_Callback(hObject, eventdata, handles)
-
-%------------------------------------------------------------------------
-%read input xml file and update the edit boxes
-function loadfile(handles,fileinput)
-%------------------------------------------------------------------------
-%read the input xml file
-t=xmltree(fileinput);
-s=convert(t);%convert to matlab structure
-%read data currently displayed on the interface
-PointCoord=[];
-Coord_cell=get(handles.ListCoord,'String');
-data=read_geometry_calib(Coord_cell);
-%data=read_geometry_calib(handles);
-Coord=[]; %default
-if isfield(data,'Coord')
-    Coord=data.Coord;
-end
-TabChar_0=get(handles.ListCoord,'String');
-nbcoord_0=size(TabChar_0,1);
-if isequal(get(handles.edit_append,'Value'),2) %edit mode  A REVOIR
-    val=get(handles.ListCoord,'Value')-1;
-else
-   val=length(TabChar_0); 
-end
-nbcoord=0;
-
-%case of calibration (ImaDoc) input file
-% hcalib=get(handles.calib_type,'parent');
-CalibData=get(handles.geometry_calib,'UserData');
-CalibData.XmlInput=fileinput;
-if isfield(s,'Heading')
-    CalibData.Heading=s.Heading;
-end
-
-set(handles.geometry_calib,'UserData',CalibData);%store the heading in the interface 'UserData'
-if isfield(s,'GeometryCalib')
-    Calib=s.GeometryCalib;
-    if isfield(Calib,'CalibrationType')
-        CalibrationType=Calib.CalibrationType;
-        switch CalibrationType
-            case 'linear'
-                set(handles.calib_type,'Value',2)
-            case 'tsai'
-                set(handles.calib_type,'Value',3)
-            case 'tsai_matlab'
-                set(handles.calib_type,'Value',4)
-        end
-    end
-    if isfield(Calib,'SourceCalib')
-        if isfield(Calib.SourceCalib,'PointCoord')
-            PointCoord=Calib.SourceCalib.PointCoord;
-        end
-    end
-    nbcoord=length(PointCoord);
-    if ~isfield(Calib,'ErrorRms')&~isfield(Calib,'ErrorMax') %old convention of Gauthier (cord in mm)
-        for i=1:length(PointCoord)
-          line=str2num(PointCoord{i});
-          Coord(i+val,4:5)=line(4:5);%px x
-          Coord(i+val,1:3)=line(1:3)/10;%phys x
-        end
-    else
-        for i=1:length(PointCoord)
-          line=str2num(PointCoord{i});
-          Coord(i,4:5)=line(4:5);%px x
-          Coord(i,1:3)=line(1:3);%phys x
-       end
-    end
-end
-%case of xml files of points 
-if isfield(s,'Coord')
-    PointCoord=s.Coord;
-    nbcoord=length(PointCoord);
-     %case of image coordinates
-    if isfield(s,'CoordType')& isequal(s.CoordType,'px')
-        for i=1:nbcoord
-           line=str2num(PointCoord{i});
-           Coord(i+val,4:5)=line(1:2);
-        end
-     %case of  physical coordinates
-    else
-        for i=1:nbcoord
-           line=str2num(PointCoord{i});
-           Coord(i+val,1:3)=line(1:3);
-           nbcolumn=size(Coord,2);
-           if nbcolumn<5
-               Coord(i+val,nbcolumn+1:5)=zeros(1,5-nbcolumn);
-           end
-        end
-     end
-end
-CoordCell={};
-for iline=1:size(Coord,1)
-    for j=1:5
-        CoordCell{iline,j}=num2str(Coord(iline,j),4);
-    end
-end
-% CoordCell=[CoordCell;{' ',' ',' ',' ',' '}];
-Tabchar=cell2tab(CoordCell,'    |    ');%transform cells into table ready for display
-Tabchar=[Tabchar;{'......'}];
-set(handles.ListCoord,'Value',1)
-set(handles.ListCoord,'String',Tabchar)
-MenuPlot_Callback(handles.geometry_calib, [], handles)
-if isempty(Coord)
-    set(handles.edit_append,'Value',1)
-    set(handles.edit_append,'BackgroundColor',[1 1 0])
-else
-    set(handles.edit_append,'Value',0)
-    set(handles.edit_append,'BackgroundColor',[0.7 0.7 0.7])
-end
 % 
 %------------------------------------------------------------------------
@@ -252 +139 @@
 function APPLY_Callback(hObject, eventdata, handles)
 %------------------------------------------------------------------------
-calib_cell=get(handles.calib_type,'String');
-val=get(handles.calib_type,'Value');
-calib_type=calib_cell{val};
+
+%read the current calibration points
 Coord_cell=get(handles.ListCoord,'String');
 Object=read_geometry_calib(Coord_cell);
-X=Object.Coord(:,1);
-Y=Object.Coord(:,2);
-Z=Object.Coord(:,3);
-if isequal(calib_type,'rescale')
-    GeometryCalib=calib_rescale(Object.Coord);
-    Z=0;%Z not taken into account
-elseif isequal(calib_type,'linear')
-    GeometryCalib=calib_linear(Object.Coord);
-    Z=0; %Z not taken into account
-elseif isequal(calib_type,'tsai_cpp')
-    GeometryCalib=calib_tsai(Object.Coord);
-elseif isequal(calib_type,'tsai_matlab')
-    GeometryCalib=calib_tsai2(Object.Coord);
-end
+Coord=Object.Coord;
+
+% apply the calibration, whose type is selected in  handles.calib_type
+if ~isempty(Coord)
+    calib_cell=get(handles.calib_type,'String');
+    val=get(handles.calib_type,'Value');
+    GeometryCalib=feval(['calib_' calib_cell{val}],Coord,handles);
+else
+    msgbox_uvmat('ERROR','No calibration points, abort')
+    return
+end 
 
 %check error
@@ -298 +181 @@
     Calib.R=GeometryCalib.R;
 end
-x_ima=Object.Coord(:,4);
-y_ima=Object.Coord(:,5); 
-[Xpoints,Ypoints]=px_XYZ(Calib,X,Y,Z);
-GeometryCalib.ErrorRms(1)=sqrt(mean((Xpoints-x_ima).*(Xpoints-x_ima)));
-[GeometryCalib.ErrorMax(1),index(1)]=max(abs(Xpoints-x_ima));
-GeometryCalib.ErrorRms(2)=sqrt(mean((Ypoints-y_ima).*(Ypoints-y_ima)));
-[GeometryCalib.ErrorMax(2),index(2)]=max(abs(Ypoints-y_ima));
-[EM,ind_dim]=max(GeometryCalib.ErrorMax);
-index=index(ind_dim);
-
+if ~isempty(Coord)
+    X=Coord(:,1);
+    Y=Coord(:,2);
+    Z=Coord(:,3);
+    x_ima=Coord(:,4);
+    y_ima=Coord(:,5);
+    [Xpoints,Ypoints]=px_XYZ(Calib,X,Y,Z);
+    GeometryCalib.ErrorRms(1)=sqrt(mean((Xpoints-x_ima).*(Xpoints-x_ima)));
+    [GeometryCalib.ErrorMax(1),index(1)]=max(abs(Xpoints-x_ima));
+    GeometryCalib.ErrorRms(2)=sqrt(mean((Ypoints-y_ima).*(Ypoints-y_ima)));
+    [GeometryCalib.ErrorMax(2),index(2)]=max(abs(Ypoints-y_ima));
+    [EM,ind_dim]=max(GeometryCalib.ErrorMax);
+    index=index(ind_dim);
+    if isequal(max(Z),min(Z))%Z constant
+        GeometryCalib.NbSlice=1;
+        GeometryCalib.SliceCoord=[0 0 Z(1)];
+    end
+end
 unitlist=get(handles.CoordUnit,'String');
 unit=unitlist{get(handles.CoordUnit,'value')};
 GeometryCalib.CoordUnit=unit;
-GeometryCalib.SourceCalib.PointCoord=Object.Coord;
+GeometryCalib.SourceCalib.PointCoord=Coord;
+
+%display calibration intrinsic parameters
+k=0;
+if isfield(GeometryCalib,'kappa1')
+    k=GeometryCalib.kappa1 * GeometryCalib.focal*GeometryCalib.focal;
+end
+sx=1;dpx_dpy=[1 1];
+if isfield(GeometryCalib,'sx')
+    sx=GeometryCalib.sx;
+end
+if isfield(GeometryCalib,'dpx_dpy')
+    dpx_dpy=GeometryCalib.dpx_dpy;
+end
+Cx_Cy=[0 0];
+if isfield(GeometryCalib,'Cx_Cy')
+    Cx_Cy=GeometryCalib.Cx_Cy;
+end
+f1=GeometryCalib.focal*sx/dpx_dpy(1);
+f2=GeometryCalib.focal/dpx_dpy(2);
+Cx=Cx_Cy(1);
+Cy=Cx_Cy(2);
+set(handles.fx,'String',num2str(f1,4))
+set(handles.fy,'String',num2str(f2,4))
+set(handles.k,'String',num2str(k,'%1.4f'))
+set(handles.Cx,'String',num2str(Cx,'%1.1f'))
+set(handles.Cy,'String',num2str(Cy,'%1.1f'))
+% Display extrinsinc parameters
+set(handles.Tx,'String',num2str(GeometryCalib.Tx_Ty_Tz(1),4))
+set(handles.Ty,'String',num2str(GeometryCalib.Tx_Ty_Tz(2),4))
+set(handles.Tz,'String',num2str(GeometryCalib.Tx_Ty_Tz(3),4))
+
+% indicate the plane of the calibration grid if defined
 huvmat=findobj(allchild(0),'Name','uvmat');
 hhuvmat=guidata(huvmat);%handles of elements in the GUI uvmat
@@ -337 +260 @@
     hhh=findobj(hhuvmat.axes3,'Tag','calib_marker');% delete calib points and markers
     if ~isempty(hhh)
-        delete(hhh);
+        delete(hhh);      
     end
     hhh=findobj(hhuvmat.axes3,'Tag','calib_points');
@@ -346 +269 @@
     set(hhuvmat.FixedLimits,'BackgroundColor',[0.7 0.7 0.7])
     uvmat('RootPath_Callback',hObject,eventdata,hhuvmat); %file input with xml reading  in uvmat
-    MenuPlot_Callback(hObject, eventdata, handles)
-    set(handles.ListCoord,'Value',index)% indicate in the list the point with max deviation (possible mistake)
-    ListCoord_Callback(hObject, eventdata, handles)
+%     if ip==0
+        MenuPlot_Callback(hObject, eventdata, handles)
+        set(handles.ListCoord,'Value',index)% indicate in the list the point with max deviation (possible mistake)
+        ListCoord_Callback(hObject, eventdata, handles)
+%     end
     figure(handles.geometry_calib)
 end
@@ -355 +280 @@
 % --- Executes on button press in calibrate_lin.
 function REPLICATE_Callback(hObject, eventdata, handles)
+% %%%%%%Todo: correct on the model of APPLY_Callback%%%%%%%%%%
 %------------------------------------------------------------------------
 calib_cell=get(handles.calib_type,'String');
 val=get(handles.calib_type,'Value');
-calib_type=calib_cell{val};
 Coord_cell=get(handles.ListCoord,'String');
 Object=read_geometry_calib(Coord_cell);
-
-if isequal(calib_type,'rescale')
-    GeometryCalib=calib_rescale(Object.Coord);
-elseif isequal(calib_type,'linear')
-    GeometryCalib=calib_linear(Object.Coord);
-elseif isequal(calib_type,'tsai')
-    GeometryCalib=calib_tsai(Object.Coord);
-end
+GeometryCalib=feval(calib_cell{val},Object.Coord);
+
 % %record image source
 GeometryCalib.SourceCalib.PointCoord=Object.Coord;
@@ -424 +343 @@
 %------------------------------------------------------------------------
 % determine the parameters for a calibration by an affine function (rescaling and offset, no rotation)
-function GeometryCalib=calib_rescale(Coord)
-%------------------------------------------------------------------------
- 
+function GeometryCalib=calib_rescale(Coord,handles)
+%------------------------------------------------------------------------
 X=Coord(:,1);
-Y=Coord(:,2);
+Y=Coord(:,2);% Z not used
 x_ima=Coord(:,4);
 y_ima=Coord(:,5);
@@ -436 +354 @@
 T_y=py(2);
 GeometryCalib.CalibrationType='rescale';
-GeometryCalib.focal=1;
+GeometryCalib.focal=1;%default
+GeometryCalib.sx=px(1)/py(1);
 GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function
-GeometryCalib.Tx_Ty_Tz=[T_x T_y 1];
-GeometryCalib.R=[px(1),0,0;0,py(1),0;0,0,1];
-
-% %check error
-% Calib.dpx=1;
-% Calib.dpy=1;
-% Calib.sx=1;
-% Calib.Cx=0;
-% Calib.Cy=0;
-% Calib.Tz=1;
-% Calib.kappa1=0;
-% Calib.f=GeometryCalib.focal;
-% Calib.Tx=T_x;
-% Calib.Ty=T_y;
-% Calib.R=GeometryCalib.R;
-% [Xpoints,Ypoints]=px_XYZ(Calib,X,Y,0);
-% GeometryCalib.ErrorRms(1)=sqrt(mean((Xpoints-x_ima).*(Xpoints-x_ima)));
-% GeometryCalib.ErrorMax(1)=max(abs(Xpoints-x_ima));
-% GeometryCalib.ErrorRms(2)=sqrt(mean((Ypoints-y_ima).*(Ypoints-y_ima)));
-% GeometryCalib.ErrorMax(2)=max(abs(Ypoints-y_ima));
+GeometryCalib.Tx_Ty_Tz=[T_x T_y GeometryCalib.focal/py(1)];
+GeometryCalib.R=[1,0,0;0,1,0;0,0,0];
 
 %------------------------------------------------------------------------
 % determine the parameters for a calibration by a linear transform matrix (rescale and rotation)
-function GeometryCalib=calib_linear(Coord)
+function GeometryCalib=calib_linear(Coord,handles)
 %------------------------------------------------------------------------
 X=Coord(:,1);
-Y=Coord(:,2);
+Y=Coord(:,2);% Z not used
 x_ima=Coord(:,4);
 y_ima=Coord(:,5);
 XY_mat=[ones(size(X)) X Y];
 a_X1=XY_mat\x_ima; %transformation matrix for X
-x1=XY_mat*a_X1;%reconstruction
-err_X1=max(abs(x1-x_ima));%error
+% x1=XY_mat*a_X1;%reconstruction
+% err_X1=max(abs(x1-x_ima));%error
 a_Y1=XY_mat\y_ima;%transformation matrix for X
-y1=XY_mat*a_Y1;
-err_Y1=max(abs(y1-y_ima));%error
-T_x=a_X1(1);
-T_y=a_Y1(1);
+% y1=XY_mat*a_Y1;
+% err_Y1=max(abs(y1-y_ima));%error
 GeometryCalib.CalibrationType='linear';
 GeometryCalib.focal=1;
 GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function
-GeometryCalib.Tx_Ty_Tz=[T_x T_y 1]; 
-GeometryCalib.R=[a_X1(2),a_X1(3),0;a_Y1(2),a_Y1(3),0;0,0,1];
-
-% %check error
-% GeometryCalib.ErrorRms(1)=sqrt(mean((x1-x_ima).*(x1-x_ima)));
-% GeometryCalib.ErrorMax(1)=max(abs(x1-x_ima));
-% GeometryCalib.ErrorRms(2)=sqrt(mean((y1-y_ima).*(y1-y_ima)));
-% GeometryCalib.ErrorMax(2)=max(abs(y1-y_ima));
-
-%------------------------------------------------------------------------
-function GeometryCalib=calib_tsai2(Coord)
+R=[a_X1(2),a_X1(3),0;a_Y1(2),a_Y1(3),0;0,0,0];
+norm=det(R);
+GeometryCalib.Tx_Ty_Tz=[a_X1(1) a_Y1(1) norm]; 
+GeometryCalib.R=R/norm;
+
+%------------------------------------------------------------------------
+% determine the tsai parameters for a view normal to the grid plane
+function GeometryCalib=calib_normal(Coord,handles)
+%------------------------------------------------------------------------
+Calib.f1=str2num(get(handles.fx,'String'));
+Calib.f2=str2num(get(handles.fy,'String'));
+Calib.k=str2num(get(handles.k,'String'));
+Calib.Cx=str2num(get(handles.Cx,'String'));
+Calib.Cy=str2num(get(handles.Cy,'String'));
+%default
+if isempty(Calib.f1)
+    Calib.f1=25/0.012;
+end
+if isempty(Calib.f2)
+    Calib.f2=25/0.012;
+end
+if isempty(Calib.k)
+    Calib.k=0;
+end
+if isempty(Calib.Cx)||isempty(Calib.Cy)
+    huvmat=findobj(allchild(0),'Tag','uvmat');
+    hhuvmat=guidata(huvmat);
+    Calib.Cx=str2num(get(hhuvmat.npx,'String'))/2;
+    Calib.Cx=str2num(get(hhuvmat.npy,'String'))/2;
+end   
+%tsai parameters
+Calib.dpx=0.012;%arbitrary
+Calib.dpy=0.012;
+Calib.sx=Calib.f1*Calib.dpx/(Calib.f2*Calib.dpy);
+Calib.f=Calib.f2*Calib.dpy;
+Calib.kappa1=Calib.k/(Calib.f*Calib.f);
+
+%initial guess
+X=Coord(:,1);
+Y=Coord(:,2);
+Zmean=mean(Coord(:,3));
+x_ima=Coord(:,4)-Calib.Cx;
+y_ima=Coord(:,5)-Calib.Cy;
+XY_mat=[ones(size(X)) X Y];
+a_X1=XY_mat\x_ima; %transformation matrix for X
+a_Y1=XY_mat\y_ima;%transformation matrix for Y
+R=[a_X1(2),a_X1(3),0;a_Y1(2),a_Y1(3),0;0,0,-1];% rotation+ z axis reversal (upward)
+norm=sqrt(det(-R));
+calib_param(1)=0;% quadratic distortion
+calib_param(2)=a_X1(1);
+calib_param(3)=a_Y1(1);
+calib_param(4)=Calib.f/(norm*Calib.dpx)-R(3,3)*Zmean;
+calib_param(5)=angle(a_X1(2)+i*a_X1(3));
+display(['initial guess=' num2str(calib_param)])
+
+%optimise the parameters: minimisation of error
+calib_param = fminsearch(@(calib_param) error_calib(calib_param,Calib,Coord),calib_param)
+
+GeometryCalib.CalibrationType='tsai_normal';
+GeometryCalib.focal=Calib.f;
+GeometryCalib.dpx_dpy=[Calib.dpx Calib.dpy];
+GeometryCalib.Cx_Cy=[Calib.Cx Calib.Cy];
+GeometryCalib.sx=Calib.sx;
+GeometryCalib.kappa1=calib_param(1);
+GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function
+GeometryCalib.Tx_Ty_Tz=[calib_param(2) calib_param(3) calib_param(4)]; 
+alpha=calib_param(5);
+GeometryCalib.R=[cos(alpha) sin(alpha) 0;-sin(alpha) cos(alpha) 0;0 0 -1];
+
+%------------------------------------------------------------------------
+function GeometryCalib=calib_3D_linear(Coord,handles)
 %------------------------------------------------------------------
 path_uvmat=which('uvmat');% check the path detected for source file uvmat
@@ -495 +454 @@
 huvmat=findobj(allchild(0),'Tag','uvmat');
 hhuvmat=guidata(huvmat);
+x_1=Coord(:,4:5)'
+X_1=Coord(:,1:3)'
+n_ima=1;
+% check_cond=0;
+
+nx=str2num(get(hhuvmat.npx,'String'));
+ny=str2num(get(hhuvmat.npy,'String'));
+
+est_dist=[0;0;0;0;0];
+est_aspect_ratio=0;
+est_fc=[1;1];
+%fc=[25;25]/0.012;
+center_optim=0;
+run(fullfile(path_UVMAT,'toolbox_calib','go_calib_optim'));
+
+GeometryCalib.CalibrationType='3D_linear';
+GeometryCalib.focal=fc(2);
+GeometryCalib.dpx_dpy=[1 1];
+GeometryCalib.Cx_Cy=cc';
+GeometryCalib.sx=fc(1)/fc(2);
+GeometryCalib.kappa1=-kc(1)/fc(2)^2;
+GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function
+GeometryCalib.Tx_Ty_Tz=Tc_1';
+GeometryCalib.R=Rc_1;
+
+%------------------------------------------------------------------------
+function GeometryCalib=calib_3D_quadr(Coord,handles)
+%------------------------------------------------------------------
+
+path_uvmat=which('uvmat');% check the path detected for source file uvmat
+path_UVMAT=fileparts(path_uvmat); %path to UVMAT
+huvmat=findobj(allchild(0),'Tag','uvmat');
+hhuvmat=guidata(huvmat);
+% check_cond=0;
+coord_files=get(handles.coord_files,'String');
+if ischar(coord_files)
+    coord_files={coord_files};
+end
+if isempty(coord_files{1}) || isequal(coord_files,{''})
+    coord_files={};
+end
+
+%retrieve the calibration points stored in the files listed in the popup list coord_files
+x_1=Coord(:,4:5)';
+X_1=Coord(:,1:3)';
+n_ima=numel(coord_files)+1;
+if ~isempty(coord_files) 
+    msgbox_uvmat('CONFIRMATION',['The xy coordinates of the calibration points in ' num2str(n_ima) ' planes will be used'])
+    for ifile=1:numel(coord_files)
+    t=xmltree(coord_files{ifile});
+    s=convert(t);%convert to matlab structure
+        if isfield(s,'GeometryCalib')
+            if isfield(s.GeometryCalib,'SourceCalib')
+                if isfield(s.GeometryCalib.SourceCalib,'PointCoord')
+                PointCoord=s.GeometryCalib.SourceCalib.PointCoord;
+                Coord_file=[];
+                for i=1:length(PointCoord)
+                    line=str2num(PointCoord{i});
+                    Coord_file(i,4:5)=line(4:5);%px x
+                    Coord_file(i,1:3)=line(1:3);%phys x
+                end
+                eval(['x_' num2str(ifile+1) '=Coord_file(:,4:5)'';']);
+                eval(['X_' num2str(ifile+1) '=Coord_file(:,1:3)'';']);
+                end
+            end
+        end
+    end
+end
+
+n_ima=numel(coord_files)+1;
+whos
+nx=str2num(get(hhuvmat.npx,'String'));
+ny=str2num(get(hhuvmat.npy,'String'));
+
+est_dist=[1;0;0;0;0];
+est_aspect_ratio=1;
+%est_fc=[0;0];
+%fc=[25;25]/0.012;
+center_optim=0;
+run(fullfile(path_UVMAT,'toolbox_calib','go_calib_optim'));
+
+GeometryCalib.CalibrationType='3D_quadr';
+GeometryCalib.focal=fc(2);
+GeometryCalib.dpx_dpy=[1 1];
+GeometryCalib.Cx_Cy=cc';
+GeometryCalib.sx=fc(1)/fc(2);
+GeometryCalib.kappa1=-kc(1)/fc(2)^2;
+GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function
+GeometryCalib.Tx_Ty_Tz=Tc_1';
+GeometryCalib.R=Rc_1;
+GeometryCalib.R(3,1:3)=-GeometryCalib.R(3,1:3);%inversion for z upward
+GeometryCalib.ErrorRMS=[];
+GeometryCalib.ErrorMax=[];
+
+
+%------------------------------------------------------------------------
+function GeometryCalib=calib_3D_extrinsic(Coord,handles)
+%------------------------------------------------------------------
+path_uvmat=which('geometry_calib');% check the path detected for source file uvmat
+path_UVMAT=fileparts(path_uvmat); %path to UVMAT
 x_1=Coord(:,4:5)';
 X_1=Coord(:,1:3)';
 n_ima=1;
-% check_cond=0;
-
-
-nx=str2num(get(hhuvmat.npx,'String'));
-ny=str2num(get(hhuvmat.npy,'String'));
-
-
-% est_kc=[1;0;0;0;0];
-est_fc=0;
-est_dist=[1;0;0;0;0];
-run(fullfile(path_UVMAT,'toolbox_calib','go_calib_optim'));
-
-GeometryCalib.CalibrationType='tsai_matlab';
-GeometryCalib.focal=f(2);
+Calib.f1=str2num(get(handles.fx,'String'));
+Calib.f2=str2num(get(handles.fy,'String'));
+Calib.k=str2num(get(handles.k,'String'));
+Calib.Cx=str2num(get(handles.Cx,'String'));
+Calib.Cy=str2num(get(handles.Cy,'String'));
+
+fct_path=fullfile(path_UVMAT,'toolbox_calib');
+addpath(fct_path)
+% [omc1,Tc1,Rc1,H,x,ex,JJ] = compute_extrinsic(x_1,X_1,...
+%     [Calib.f Calib.f*Calib.sx]',...
+%     [Calib.Cx Calib.Cy]',...
+%     [-Calib.kappa1*Calib.f^2 0 0 0 0]);
+[omc1,Tc1,Rc1,H,x,ex,JJ] = compute_extrinsic(x_1,X_1,...
+    [Calib.f1 Calib.f2]',...
+    [Calib.Cx Calib.Cy]',...
+    [-Calib.k 0 0 0 0]);
+%get the euler angles ???
+rmpath(fct_path);
+
+std(ex')
+GeometryCalib.CalibrationType='3D_extrinsic';
+GeometryCalib.focal=Calib.f2;
 GeometryCalib.dpx_dpy=[1 1];
-GeometryCalib.Cx_Cy=cc';
-GeometryCalib.sx=fc(1)/fc(2);
-GeometryCalib.kappa1=-k(1)/f(2)^2;
+GeometryCalib.Cx_Cy=[Calib.Cx Calib.Cy];
+GeometryCalib.sx=Calib.f1/Calib.f2;
+GeometryCalib.kappa1=Calib.k/(Calib.f2*Calib.f2);
 GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function
-GeometryCalib.Tx_Ty_Tz=Tc_1';
-GeometryCalib.R=Rc_1;
-% Calib.dpx=GeometryCalib.dpx_dpy(1);
-% Calib.dpy=GeometryCalib.dpx_dpy(2);
-% Calib.sx=GeometryCalib.sx;
-% Calib.Cx=GeometryCalib.Cx_Cy(1);
-% Calib.Cy=GeometryCalib.Cx_Cy(2);
-% Calib.kappa1=GeometryCalib.kappa1;
-% Calib.f=GeometryCalib.focal;
-% Calib.Tx=GeometryCalib.Tx_Ty_Tz(1);
-% Calib.Ty=GeometryCalib.Tx_Ty_Tz(2);
-% Calib.Tz=GeometryCalib.Tx_Ty_Tz(3);
-% Calib.R=GeometryCalib.R;
-% X=Coord(:,1);
-% Y=Coord(:,2);
-% Z=Coord(:,3);
-% x_ima=Coord(:,4);
-% y_ima=Coord(:,5);
-% [Xpoints,Ypoints]=px_XYZ(Calib,X,Y,Z);
-% 
-% GeometryCalib.ErrorRms(1)=sqrt(mean((Xpoints-x_ima).*(Xpoints-x_ima)));
-% [GeometryCalib.ErrorMax(1),GeometryCalib.IndexMax(1)]=max(abs(Xpoints-x_ima));
-% GeometryCalib.ErrorRms(2)=sqrt(mean((Ypoints-y_ima).*(Ypoints-y_ima)));
-% [GeometryCalib.ErrorMax(2),GeometryCalib.IndexMax(2)]=max(abs(Ypoints-y_ima));
-
-function GeometryCalib=calib_tsai(Coord)
+GeometryCalib.Tx_Ty_Tz=Tc1';
+%inversion of z axis 
+GeometryCalib.R=Rc1;
+%GeometryCalib.R(3,1:3)=-GeometryCalib.R(3,1:3);%inversion for z upward
+
+
+%------------------------------------------------------------------------
+%function GeometryCalib=calib_tsai_heikkila(Coord)
+% TEST: NOT IMPLEMENTED
+%------------------------------------------------------------------
+% path_uvmat=which('uvmat');% check the path detected for source file uvmat
+% path_UVMAT=fileparts(path_uvmat); %path to UVMAT
+% path_calib=fullfile(path_UVMAT,'toolbox_calib_heikkila');
+% addpath(path_calib)
+% npoints=size(Coord,1);
+% Coord(:,1:3)=10*Coord(:,1:3);
+% Coord=[Coord zeros(npoints,2) -ones(npoints,1)];
+% [par,pos,iter,res,er,C]=cacal('dalsa',Coord);
+% GeometryCalib.CalibrationType='tsai';
+% GeometryCalib.focal=par(2);
+
+
+%--------------------------------------------------------------------------
+function GeometryCalib=calib_tsai(Coord,handles)% old version using gauthier's bianry ccal_fo
 %------------------------------------------------------------------------
 %TSAI
-% 'calibration_lin' provides a linear transform on coordinates, 
 path_uvmat=which('uvmat');% check the path detected for source file uvmat
 path_UVMAT=fileparts(path_uvmat); %path to UVMAT
-% if isunix
-    %fid = fopen(fullfile(path_UVMAT,'PARAM_LINUX.txt'),'r');%open the file with civ binary names
-xmlfile=fullfile(path_UVMAT,'PARAM.xml');
+xmlfile=fullfile(path_UVMAT,'PARAM.xml');%name of the file containing names of binary executables
 if exist(xmlfile,'file')
-    t=xmltree(xmlfile);
-    sparam=convert(t);
+    t=xmltree(xmlfile);% read the (xml) file containing names of binary executables
+    sparam=convert(t);% convert to matlab structure
 end
 if ~isfield(sparam,'GeometryCalibBin')
@@ -570 +632 @@
 textcoord=num2str(Coord,4);
 dlmwrite('t.txt',textcoord,'');  
-% ['!' Tsai_exe ' -f1 0 -f2 t.txt']
-    eval(['!' Tsai_exe ' -f t.txt > tsaicalib.log']);
+% ['!' Tsai_exe ' -fx 0 -fy t.txt']
+eval(['!' Tsai_exe ' -f t.txt > tsaicalib.log']);
 if ~exist('calib.dat','file')
     msgbox_uvmat('ERROR','no output from calibration program Tsai_exe: possibly too few points')
@@ -577 +639 @@
 calibdat=dlmread('calib.dat');
 delete('calib.dat')
-delete('t.txt')
+%delete('t.txt')
 GeometryCalib.CalibrationType='tsai';
 GeometryCalib.focal=calibdat(10);
@@ -604 +666 @@
 %EN DEDUIRE MATRICE R ??
 GeometryCalib.R=[r1,r2,r3;r4,r5,r6;r7,r8,r9];
-%erreur a caracteriser?
-% %check error
-% Calib.dpx=GeometryCalib.dpx_dpy(1);
-% Calib.dpy=GeometryCalib.dpx_dpy(2);
-% Calib.sx=GeometryCalib.sx;
-% Calib.Cx=GeometryCalib.Cx_Cy(1);
-% Calib.Cy=GeometryCalib.Cx_Cy(2);
-% Calib.kappa1=GeometryCalib.kappa1;
-% Calib.f=GeometryCalib.focal;
-% Calib.Tx=GeometryCalib.Tx_Ty_Tz(1);
-% Calib.Ty=GeometryCalib.Tx_Ty_Tz(2);
-% Calib.Tz=GeometryCalib.Tx_Ty_Tz(3);
-% Calib.R=GeometryCalib.R;
-% X=Coord(:,1);
-% Y=Coord(:,2);
-% Z=Coord(:,3);
-% x_ima=Coord(:,4);
-% y_ima=Coord(:,5);
-% [Xpoints,Ypoints]=px_XYZ(Calib,X,Y,Z);
-% 
-% GeometryCalib.ErrorRms(1)=sqrt(mean((Xpoints-x_ima).*(Xpoints-x_ima)));
-% GeometryCalib.ErrorMax(1)=max(abs(Xpoints-x_ima));
-% GeometryCalib.ErrorRms(2)=sqrt(mean((Ypoints-y_ima).*(Ypoints-y_ima)));
-% GeometryCalib.ErrorMax(2)=max(abs(Ypoints-y_ima));
-
-
-%------------------------------------------------------------------------
-% --- Executes on button press in rotation.
-function rotation_Callback(hObject, eventdata, handles)
-%------------------------------------------------------------------------
-angle_rot=(pi/180)*str2num(get(handles.Phi,'String'));
+
+%------------------------------------------------------------------------
+% --- determine the rms of calibration error
+function ErrorRms=error_calib(calib_param,Calib,Coord)
+%calib_param: vector of free calibration parameters (to optimise)
+%Calib: structure of the given calibration parameters
+%Coord: list of phys coordinates (columns 1-3, and pixel coordinates (columns 4-5)
+Calib.f=25;
+Calib.dpx=0.012;
+Calib.dpy=0.012;
+Calib.sx=1;
+Calib.Cx=512;
+Calib.Cy=512;
+Calib.kappa1=calib_param(1);
+Calib.Tx=calib_param(2);
+Calib.Ty=calib_param(3);
+Calib.Tz=calib_param(4);
+alpha=calib_param(5);
+Calib.R=[cos(alpha) sin(alpha) 0;-sin(alpha) cos(alpha) 0;0 0 -1];
+
+X=Coord(:,1);
+Y=Coord(:,2);
+Z=Coord(:,3);
+x_ima=Coord(:,4);
+y_ima=Coord(:,5); 
+[Xpoints,Ypoints]=px_XYZ(Calib,X,Y,Z);
+ErrorRms(1)=sqrt(mean((Xpoints-x_ima).*(Xpoints-x_ima)));
+ErrorRms(2)=sqrt(mean((Ypoints-y_ima).*(Ypoints-y_ima)));
+ErrorRms=mean(ErrorRms);
+
+%------------------------------------------------------------------------
+function XImage_Callback(hObject, eventdata, handles)
+%------------------------------------------------------------------------
+update_list(hObject, eventdata,handles)
+
+%------------------------------------------------------------------------
+function YImage_Callback(hObject, eventdata, handles)
+%------------------------------------------------------------------------
+update_list(hObject, eventdata,handles)
+
+%------------------------------------------------------------------------
+% --- Executes on button press in STORE.
+function STORE_Callback(hObject, eventdata, handles)
 Coord_cell=get(handles.ListCoord,'String');
-data=read_geometry_calib(Coord_cell);
-data.Coord(:,1)=cos(angle_rot)*data.Coord(:,1)+sin(angle_rot)*data.Coord(:,2);
-data.Coord(:,1)=-sin(angle_rot)*data.Coord(:,1)+cos(angle_rot)*data.Coord(:,2);
-set(handles.XObject,'String',num2str(data.Coord(:,1),4));
-set(handles.YObject,'String',num2str(data.Coord(:,2),4));
-
-%------------------------------------------------------------------------
-function XImage_Callback(hObject, eventdata, handles)
+Object=read_geometry_calib(Coord_cell);
+unitlist=get(handles.CoordUnit,'String');
+unit=unitlist{get(handles.CoordUnit,'value')};
+GeometryCalib.CoordUnit=unit;
+GeometryCalib.SourceCalib.PointCoord=Object.Coord;
+huvmat=findobj(allchild(0),'Name','uvmat');
+hhuvmat=guidata(huvmat);%handles of elements in the GUI uvmat
+RootPath='';
+RootFile='';
+if ~isempty(hhuvmat.RootPath)& ~isempty(hhuvmat.RootFile)
+    testhandle=1;
+    RootPath=get(hhuvmat.RootPath,'String');
+    RootFile=get(hhuvmat.RootFile,'String');
+    filebase=fullfile(RootPath,RootFile);
+    while exist([filebase '.xml'],'file')
+        filebase=[filebase '~'];
+    end
+    outputfile=[filebase '.xml'];
+    update_imadoc(GeometryCalib,outputfile)
+    listfile=get(handles.coord_files,'string');
+    if isequal(listfile,{''})
+        listfile={outputfile};
+    else
+        listfile=[listfile;{outputfile}]%update the list of coord files
+    end
+    set(handles.coord_files,'string',listfile);
+end
+set(handles.ListCoord,'Value',1)% refresh the display of coordinates
+set(handles.ListCoord,'String',{'......'})
+
+%------------------------------------------------------------------------
+% --- Executes on button press in CLEAR.
+function CLEAR_Callback(hObject, eventdata, handles)
+%------------------------------------------------------------------------
+set(handles.coord_files,'Value',1)
+set(handles.coord_files,'String',{''})
+
+%------------------------------------------------------------------------
+function XObject_Callback(hObject, eventdata, handles)
 %------------------------------------------------------------------------
 update_list(hObject, eventdata,handles)
 
 %------------------------------------------------------------------------
-function YImage_Callback(hObject, eventdata, handles)
+function YObject_Callback(hObject, eventdata, handles)
 %------------------------------------------------------------------------
 update_list(hObject, eventdata,handles)
 
-function XObject_Callback(hObject, eventdata, handles)
-update_list(hObject, eventdata,handles)
-
-function YObject_Callback(hObject, eventdata, handles)
-update_list(hObject, eventdata,handles)
-
+%------------------------------------------------------------------------
 function ZObject_Callback(hObject, eventdata, handles)
+%------------------------------------------------------------------------
 update_list(hObject, eventdata,handles)
 
@@ -748 +849 @@
 %------------------------------------------------------------------------
 choice=get(handles.edit_append,'Value');
-% if choice==1
-%        Coord=get(handles.ListCoord,'String'); 
-%        val=length(Coord);
-%        if val>=1 & isequal(Coord{val},'')
-%             val=val-1; %do not take into account blank
-%        end
-%        Coord{val+1}='';
-%        set(handles.ListCoord,'String',Coord)
-%        set(handles.ListCoord,'Value',val+1)
-% end
 if choice
     set(handles.edit_append,'BackgroundColor',[1 1 0])
@@ -817 +908 @@
 end
 
-% %------------------------------------------------------------------------
-% % --- Executes on button press in append_point.
-% function append_point_Callback(hObject, eventdata, handles)
-% %------------------------------------------------------------------------
-%        Coord=get(handles.ListCoord,'String'); 
-%        val=length(Coord);
-%        if val>=1 & isequal(Coord{val},'')
-%             val=val-1; %do not take into account blank
-%        end
-%        Coord{val+1}='';
-%        set(handles.ListCoord,'String',Coord)
-%        set(handles.ListCoord,'Value',val+1)
-
-%------------------------------------------------------------------------
-function MenuOpen_Callback(hObject, eventdata, handles)
-%------------------------------------------------------------------------
-%get the object file 
-huvmat=findobj(allchild(0),'Name','uvmat');
-UvData=get(huvmat,'UserData');
-hchild=get(huvmat,'Children');
-hrootpath=findobj(hchild,'Tag','RootPath');
-oldfile=get(hrootpath,'String');
-if isempty(oldfile)
-    oldfile='';
-end
-%[FileName,PathName] = uigetfile('*.civ','Select a .civ file',oldfile)
-[FileName, PathName, filterindex] = uigetfile( ...
-       {'*.xml;*.mat', ' (*.xml,*.mat)';
-       '*.xml',  '.xml files '; ...
-        '*.mat',  '.mat matlab files '}, ...
-        'Pick a file',oldfile);
-fileinput=[PathName FileName];%complete file name 
-testblank=findstr(fileinput,' ');%look for blanks
-if ~isempty(testblank)
-    msgbox_uvmat('ERROR','forbidden input file name or path: no blank character allowed')
-    return
-end
-sizf=size(fileinput);
-if (~ischar(fileinput)|~isequal(sizf(1),1)),return;end
-loadfile(handles,fileinput)
-
-
 %------------------------------------------------------------------------
 function MenuPlot_Callback(hObject, eventdata, handles)
@@ -955 +1004 @@
     grid_input=CalibData.grid;%retrieve the previously used grid
 end
-[T,CalibData.grid]=create_grid(grid_input,'detect grid');%display the GUI create_grid, read the set of phys coordinates T
+[T,CalibData.grid]=create_grid(grid_input,'detect_grid');%display the GUI create_grid, read the set of phys coordinates T
+
 set(handles.geometry_calib,'UserData',CalibData)%store the phys grid parameters for later use
 
@@ -988 +1038 @@
 X=[CalibData.grid.x_0 CalibData.grid.x_1 CalibData.grid.x_0 CalibData.grid.x_1]';%corner absissa in the phys coordinates
 Y=[CalibData.grid.y_0 CalibData.grid.y_0 CalibData.grid.y_1 CalibData.grid.y_1]';%corner ordinates in the phys coordinates
-%X=[CalibData.grid.x_0 CalibData.grid.x_0 CalibData.grid.x_1 CalibData.grid.x_1]';%corner absissa in the phys coordinates
-%Y=[CalibData.grid.y_0 CalibData.grid.y_1 CalibData.grid.y_1 CalibData.grid.y_0]';%corner ordinates in the phys coordinate
-%XY_mat=[ones(size(X)) X Y];
-%a_X1=XY_mat\corners_X; %transformation matrix for X
-%x1=XY_mat*a_X1;%reconstruction
-%err_X1=max(abs(x1-corners_X));%error
-%a_Y1=XY_mat\corners_Y;%transformation matrix for X
-%y1=XY_mat*a_Y1;
-%err_Y1=max(abs(y1-corners_Y));%error
-% figure(1)
-% hold off
-% plot([corners_X;corners_X(1)],[corners_Y;corners_Y(1)],'r')
-%test points along contour
-% u=[0:0.1:1];
-% x12=u*corners_X(2)+(1-u)*corners_X(1);
-% x23=u*corners_X(3)+(1-u)*corners_X(2);
-% x34=u*corners_X(4)+(1-u)*corners_X(3);
-% x41=u*corners_X(1)+(1-u)*corners_X(4);
-% y12=u*corners_Y(2)+(1-u)*corners_Y(1);
-% y23=u*corners_Y(3)+(1-u)*corners_Y(2);
-% y34=u*corners_Y(4)+(1-u)*corners_Y(3);
-% y41=u*corners_Y(1)+(1-u)*corners_Y(4);
-% hold on
-% plot(x12,y12,'xr')
-% plot(x23,y23,'xr')
-% plot(x34,y34,'xr')
-% plot(x41,y41,'xr')
+
 %calculate transform matrices: 
 % reference: http://alumni.media.mit.edu/~cwren/interpolator/ by Christopher R. Wren
@@ -1026 +1050 @@
 C = [l(7:8)' 1];
 
-GeometryCalib.CalibrationType='tsai';%'linear';
+GeometryCalib.CalibrationType='tsai';
 GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function
 GeometryCalib.f=1;
@@ -1039 +1063 @@
 GeometryCalib.Tz=1;
 GeometryCalib.R=[Amat(1,1),Amat(1,2),0;Amat(2,1),Amat(2,2),0;C(1),C(2),0];
-%montre points du pourtour transformes 
-% for ip=1:numel(u),
-%   XY12(:,ip)=Amat*[x12(ip);y12(ip);1]/(C*[x12(ip);y12(ip);1]);
-%   XY23(:,ip)=Amat*[x23(ip);y23(ip);1]/(C*[x23(ip);y23(ip);1]);
-%   XY34(:,ip)=Amat*[x34(ip);y34(ip);1]/(C*[x34(ip);y34(ip);1]);
-%   XY41(:,ip)=Amat*[x41(ip);y41(ip);1]/(C*[x41(ip);y41(ip);1]);
-% end
-% xphys12=u*X(2)+(1-u)*X(1);
-% xphys23=u*X(3)+(1-u)*X(2);
-% xphys34=u*X(4)+(1-u)*X(3);
-% xphys41=u*X(1)+(1-u)*X(4);
-% yphys12=u*Y(2)+(1-u)*Y(1);
-% yphys23=u*Y(3)+(1-u)*Y(2);
-% yphys34=u*Y(4)+(1-u)*Y(3);
-% yphys41=u*Y(1)+(1-u)*Y(4);
-% for ip=1:numel(u),
-%   XY12(:,ip)=Amat*[xphys12(ip);yphys12(ip);1]/(C*[xphys12(ip);yphys12(ip);1]);
-%   XY23(:,ip)=Amat*[xphys23(ip);yphys23(ip);1]/(C*[xphys23(ip);yphys23(ip);1]);
-%   XY34(:,ip)=Amat*[xphys34(ip);yphys34(ip);1]/(C*[xphys34(ip);yphys34(ip);1]);
-%   XY41(:,ip)=Amat*[xphys41(ip);yphys41(ip);1]/(C*[xphys41(ip);yphys41(ip);1]);
-% end
-% x12=XY12(1,:);
-% y12=XY12(2,:);
-% x23=XY23(1,:);
-% y23=XY23(2,:);
-% x34=XY34(1,:);
-% y34=XY34(2,:);
-% x41=XY41(1,:);
-% y41=XY41(2,:);
-% [x12,y12]=px_XYZ(GeometryCalib,xphys12,yphys12);
-% [x23,y23]=px_XYZ(GeometryCalib,xphys23,yphys23);
-% [x34,y34]=px_XYZ(GeometryCalib,xphys34,yphys34);
-% [x41,y41]=px_XYZ(GeometryCalib,xphys41,yphys41);
-% plot(x12,y12,'xb')
-% plot(x23,y23,'xb')
-% plot(x34,y34,'xb')
-% plot(x41,y41,'xb')
-% figure(2)
-% hold off
-% plot(XY12(1,:),XY12(2,:),'ob')  
-% hold on
-% plot(XY23(1,:),XY23(2,:),'ob') 
-% plot(XY34(1,:),XY34(2,:),'ob')
-% plot(XY41(1,:),XY41(2,:),'ob')
-% plot(xphys12,yphys12,'ob')  
-% hold on
-% plot(xphys23,yphys23,'ob') 
-% plot(xphys34,yphys34,'ob')
-% plot(xphys41,yphys41,'ob')
 
 [Amod,Rangx,Rangy]=phys_Ima(A-min(min(A)),GeometryCalib,0);
 
 Amod=double(Amod);
-figure(12)
-Amax=max(max(Amod));
-image(Rangx,Rangy,uint8(255*Amod/Amax))
+% figure(12) display corrected image
+% Amax=max(max(Amod));
+% image(Rangx,Rangy,uint8(255*Amod/Amax))
 Dx=(Rangx(2)-Rangx(1))/(npxy(2)-1); %x mesh in real space
 Dy=(Rangy(2)-Rangy(1))/(npxy(1)-1); %y mesh in real space
@@ -1103 +1078 @@
     i0=1+round((T(ipoint,1)-Rangx(1))/Dx);%round(Xpx(ipoint));
     j0=1+round((T(ipoint,2)-Rangy(1))/Dy);%round(Xpx(ipoint));
-    Asub=Amod(j0-ind_range_y:j0+ind_range_y,i0-ind_range_x:i0+ind_range_x);
+    j0min=max(j0-ind_range_y,1);
+    j0max=min(j0+ind_range_y,size(Amod,1));
+    i0min=max(i0-ind_range_x,1);
+    i0max=min(i0+ind_range_x,size(Amod,2));
+    Asub=Amod(j0min:j0max,i0min:i0max);
     x_profile=sum(Asub,1);
     y_profile=sum(Asub,2);
@@ -1111 +1090 @@
     x_shift=0;
     y_shift=0;
-    if ind_x_max+2<=2*ind_range_x+1 && ind_x_max-2>=1
+    %if ind_x_max+2<=2*ind_range_x+1 && ind_x_max-2>=1
+    if ind_x_max+2<=numel(x_profile) && ind_x_max-2>=1
         Atop=x_profile(ind_x_max-2:ind_x_max+2);
         x_shift=sum(Atop.*[-2 -1 0 1 2])/sum(Atop);
     end
-    if ind_y_max+2<=2*ind_range_y+1 && ind_y_max-2>=1
+    if ind_x_max+2<=numel(y_profile) && ind_y_max-2>=1
         Atop=y_profile(ind_y_max-2:ind_y_max+2);
         y_shift=sum(Atop.*[-2 -1 0 1 2]')/sum(Atop);
     end
-    Delta(ipoint,1)=(x_shift+ind_x_max-ind_range_x-1)*Dx;%shift from the initial guess
-    Delta(ipoint,2)=(y_shift+ind_y_max-ind_range_y-1)*Dy;
+    Delta(ipoint,1)=(x_shift+ind_x_max+i0min-i0-1)*Dx;%shift from the initial guess
+    Delta(ipoint,2)=(y_shift+ind_y_max+j0min-j0-1)*Dy;
 end
 Tmod=T(:,(1:2))+Delta;
@@ -1127 +1107 @@
      Coord{ipoint,1}=num2str(T(ipoint,1),4);%display coordiantes with 4 digits
      Coord{ipoint,2}=num2str(T(ipoint,2),4);%display coordiantes with 4 digits
-     Coord{ipoint,3}='0';
-     Coord{ipoint,4}=num2str(Xpx(ipoint),4);%display coordiantes with 4 digi
-     Coord{ipoint,5}=num2str(Ypx(ipoint),4);%display coordiantes with 4 digi
+     Coord{ipoint,3}=num2str(T(ipoint,3),4);%display coordiantes with 4 digits;
+     Coord{ipoint,4}=num2str(Xpx(ipoint),4);%display coordiantes with 4 digits
+     Coord{ipoint,5}=num2str(Ypx(ipoint),4);%display coordiantes with 4 digits
 end
 Tabchar=cell2tab(Coord(end:-1:1,:),'    |    ');
@@ -1209 +1189 @@
 set(handles.ListCoord,'Value',1)
 set(handles.ListCoord,'String',Tabchar)
+
+
+% %------------------------------------------------------------------------
+% % --- Executes on button press in rotation.
+% function rotation_Callback(hObject, eventdata, handles)
+% %------------------------------------------------------------------------
+% angle_rot=(pi/180)*str2num(get(handles.Phi,'String'));
+% Coord_cell=get(handles.ListCoord,'String');
+% data=read_geometry_calib(Coord_cell);
+% data.Coord(:,1)=cos(angle_rot)*data.Coord(:,1)+sin(angle_rot)*data.Coord(:,2);
+% data.Coord(:,1)=-sin(angle_rot)*data.Coord(:,1)+cos(angle_rot)*data.Coord(:,2);
+% set(handles.XObject,'String',num2str(data.Coord(:,1),4));
+% set(handles.YObject,'String',num2str(data.Coord(:,2),4));
+
 
 %------------------------------------------------------------------------
@@ -1287 +1281 @@
     Zphys=Calib.SliceCoord(Zindex,3);%GENERALISER AUX CAS AVEC ANGLE
 else
-%     if exist('Z','var')
-%         Zphys=Z;
-%     else
-        Zphys=0;
-%     end
+    Zphys=0;
 end
 if ~exist('X','var')||~exist('Y','var')
@@ -1329 +1319 @@
     Yphys=(A21.*Xu+A22.*Yu+Y0)./denom;
 end
+
+
+% --------------------------------------------------------------------
+function MenuImportPoints_Callback(hObject, eventdata, handles)
+fileinput=browse_xml(hObject, eventdata, handles)
+if isempty(fileinput)
+    return
+end
+[s,errormsg]=imadoc2struct(fileinput,'GeometryCalib');
+GeometryCalib=s.GeometryCalib;
+GeometryCalib=load_calib(hObject, eventdata, handles)
+calib=reshape(GeometryCalib.PointCoord',[],1);
+for ilist=1:numel(calib)
+    CoordCell{ilist}=num2str(calib(ilist));
+end
+CoordCell=reshape(CoordCell,[],5);
+Tabchar=cell2tab(CoordCell,'    |    ');%transform cells into table ready for display
+Tabchar=[Tabchar;{'......'}];
+set(handles.ListCoord,'Value',1)
+set(handles.ListCoord,'String',Tabchar)
+MenuPlot_Callback(handles.geometry_calib, [], handles)
+
+% -----------------------------------------------------------------------
+function MenuImportIntrinsic_Callback(hObject, eventdata, handles)
+%------------------------------------------------------------------------
+fileinput=browse_xml(hObject, eventdata, handles);
+if isempty(fileinput)
+    return
+end
+[s,errormsg]=imadoc2struct(fileinput,'GeometryCalib');
+GeometryCalib=s.GeometryCalib;
+k=GeometryCalib.kappa1 * GeometryCalib.f*GeometryCalib.f;
+f1=GeometryCalib.f*GeometryCalib.sx/GeometryCalib.dpx;
+f2=GeometryCalib.f/GeometryCalib.dpy;
+Cx=GeometryCalib.Cx;
+Cy=GeometryCalib.Cy;
+set(handles.fx,'String',num2str(f1,'%1.1f'))
+set(handles.fy,'String',num2str(f2,'%1.1f'))
+set(handles.k,'String',num2str(k,'%1.4f'))
+set(handles.Cx,'String',num2str(Cx,'%1.1f'))
+set(handles.Cy,'String',num2str(Cy,'%1.1f'))
+
+% -----------------------------------------------------------------------
+function MenuImportAll_Callback(hObject, eventdata, handles)
+%------------------------------------------------------------------------
+fileinput=browse_xml(hObject, eventdata, handles)
+if ~isempty(fileinput)
+    loadfile(handles,fileinput)
+end
+
+% -----------------------------------------------------------------------
+% --- Executes on menubar option Import/Grid file: introduce previous grid files
+function MenuGridFile_Callback(hObject, eventdata, handles)
+% -----------------------------------------------------------------------
+inputfile=browse_xml(hObject, eventdata, handles)
+listfile=get(handles.coord_files,'string');
+if isequal(listfile,{''})
+    listfile={inputfile};
+else
+    listfile=[listfile;{inputfile}]%update the list of coord files
+end
+set(handles.coord_files,'string',listfile);
+
+%------------------------------------------------------------------------
+function fileinput=browse_xml(hObject, eventdata, handles)
+%------------------------------------------------------------------------
+fileinput=[];%default
+oldfile=''; %default
+UserData=get(handles.geometry_calib,'UserData');
+if isfield(UserData,'XmlInput')
+    oldfile=UserData.XmlInput;
+end
+[FileName, PathName, filterindex] = uigetfile( ...
+       {'*.xml;*.mat', ' (*.xml,*.mat)';
+       '*.xml',  '.xml files '; ...
+        '*.mat',  '.mat matlab files '}, ...
+        'Pick a file',oldfile);
+fileinput=[PathName FileName];%complete file name 
+testblank=findstr(fileinput,' ');%look for blanks
+if ~isempty(testblank)
+    msgbox_uvmat('ERROR','forbidden input file name or path: no blank character allowed')
+    return
+end
+sizf=size(fileinput);
+if (~ischar(fileinput)||~isequal(sizf(1),1)),return;end
+UserData.XmlInput=fileinput;
+set(handles.geometry_calib,'UserData',UserData)%record current file foer further use of browser
+
+% -----------------------------------------------------------------------
+function loadfile(handles,fileinput)
+%------------------------------------------------------------------------
+[s,errormsg]=imadoc2struct(fileinput,'GeometryCalib');
+GeometryCalib=s.GeometryCalib;
+if isempty(GeometryCalib)
+    Tabchar={};
+    CoordCell={};
+    k=0;%default
+    f1=1000;
+    f2=1000;
+    hhuvmat=guidata(findobj(allchild(0),'Name','uvmat'));
+    Cx=str2num(get(hhuvmat.npx,'String'))/2;
+    Cy=str2num(get(hhuvmat.npy,'String'))/2;
+else
+    k=GeometryCalib.kappa1 * GeometryCalib.f*GeometryCalib.f;
+    f1=GeometryCalib.f*GeometryCalib.sx/GeometryCalib.dpx;
+    f2=GeometryCalib.f/GeometryCalib.dpy;
+    Cx=GeometryCalib.Cx;
+    Cy=GeometryCalib.Cy;
+    set(handles.fx,'String',num2str(f1,'%1.1f'))
+    set(handles.fy,'String',num2str(f2,'%1.1f'))
+    set(handles.k,'String',num2str(k,'%1.4f'))
+    set(handles.Cx,'String',num2str(Cx,'%1.1f'))
+    set(handles.Cy,'String',num2str(Cy,'%1.1f'))
+    calib=reshape(GeometryCalib.PointCoord',[],1);
+    for ilist=1:numel(calib)
+        CoordCell{ilist}=num2str(calib(ilist));
+    end
+    CoordCell=reshape(CoordCell,[],5);
+    Tabchar=cell2tab(CoordCell,'    |    ');%transform cells into table ready for display
+    set(handles.ListCoord,'Value',1)
+    set(handles.ListCoord,'String',Tabchar)
+    MenuPlot_Callback(handles.geometry_calib, [], handles)
+end
+Tabchar=[Tabchar;{'......'}];
+if isempty(CoordCell)% allow mouse action by default in the absence of input points
+    set(handles.edit_append,'Value',1)
+    set(handles.edit_append,'BackgroundColor',[1 1 0])
+else % does not allow mouse action by default in the presence of input points
+    set(handles.edit_append,'Value',0)
+    set(handles.edit_append,'BackgroundColor',[0.7 0.7 0.7])
+end
+
+
+
+
+
+% --------------------------------------------------------------------
+% --- Executes on button press in CLEAR_PTS: clear the list of calibration points
+function CLEAR_PTS_Callback(hObject, eventdata, handles)
+% --------------------------------------------------------------------
+set(handles.ListCoord,'Value',1)% refresh the display of coordinates
+set(handles.ListCoord,'String',{'......'})
+

trunk/src/imadoc2struct.m

@@ -1 +1 @@
 %'imadoc2struct': reads the xml file for image documentation 
 %------------------------------------------------------------------------
-% function [s,errormsg]=imadoc2struct(ImaDoc) 
+% function [s,errormsg]=imadoc2struct(ImaDoc,option) 
 %
 % OUTPUT:
@@ -13 +13 @@
 % INPUT:
 % ImaDoc: full name of the xml input file with head key ImaDoc
-
-function [s,errormsg]=imadoc2struct(ImaDoc) 
-
+% option: ='GeometryCalib': read  the data of GeometryCalib, including source point coordinates
+
+function [s,errormsg]=imadoc2struct(ImaDoc,option) 
+if ~exist('option','var')
+    option='default';
+end
 errormsg=[];%default
 s.Heading=[];%default
@@ -88 +91 @@
                 Dtj=reshape(Dtj'*ones(1,NbDtj),1,NbDtj*numel(Dtj)); %concatene Dti vector NbDti times
                 Dtj=[0 Dtj];
-%                 Time_val'
-%                 ones(1,numel(Dtj))
-%                 ones(numel(Time_val'),1)
-%                 cumsum(Dtj)
                 Time_val=Time_val*ones(1,numel(Dtj))+ones(numel(Time_val),1)*cumsum(Dtj);% produce a time matrix with Dtj
             end
@@ -107 +106 @@
         end
     end
-%     if size(s.Time,1)==1
-%         s.Time=(s.Time)'; %change vector into column
-%     end
 end
 
@@ -210 +206 @@
             end         
         end
+        if strcmp(option,'GeometryCalib')
+            tsai.PointCoord=get_value(subt,'/GeometryCalib/SourceCalib/PointCoord',[0 0 0 0 0]);%time in TimeUnit
+        end
         s.GeometryCalib=tsai;
     end
@@ -228 +227 @@
                val_read=str2num(get(t,uid_child(icell),'value'));
                if ~isempty(val_read)
-                   val(icell)=val_read;
+                   val(icell,:)=val_read;
                end
            end

trunk/src/plot_field.m

@@ -243 +243 @@
 function hdisplay=plot_text(FieldData,hdisplay_in)
 %-------------------------------------------------------------------
-hdisplay_in
 if exist('hdisplay_in','var') && ~isempty(hdisplay_in) && ishandle(hdisplay_in) && isequal(get(hdisplay_in,'Type'),'uicontrol')
     hdisplay=hdisplay_in;
@@ -687 +686 @@
         end
           x_label=[Data.ListVarName{ivar_X} '(' x_units ')'];
-    end       
-%     if isfield(Data,'VarAttribute')
-%         VarAttribute=Data.VarAttribute;
-%     end    
+    end          
 end 
 
@@ -704 +700 @@
     siz=numel(np);
     if siz>3
-        msgbox_uvmat('ERROR','unrecognized scalar type ')
+        msgbox_uvmat('ERROR',['unrecognized scalar type: ' num2str(siz) ' dimensions'])
             return
     end
@@ -794 +790 @@
         if ~isequal(PlotParam.Scalar.Contours,1)  
             % rescale the grey levels with min and max, put a grey scale colorbar
+            B=A;
             if BW
-                B=A;
                 vec=linspace(0,1,255);%define a linear greyscale colormap
                 map=[vec' vec' vec'];
                 colormap(map);  %grey scale color map 
             else
-                B=A;
                 colormap('default'); % standard faulse colors for div, vort , scalar fields 
             end

trunk/src/proj_field.m

@@ -12 +12 @@
 %  .NbPix;
 %  .DimName=  names of the matrix dimensions (matlab cell)
-%  .DimValue= values of the matricx dimensions (matlab vector, same length as .DimName);
 %  .VarName= names of the variables [ProjData.VarName {'A','AMean','AMin','AMax'}];
 %  .VarDimNameIndex= dimensions of the variables, indicated by indices in the list .DimName;
@@ -34 +33 @@
 %         .ListGlobalAttribute: cell listing the names of the global attributes
 %        .Att_1,Att_2... : values of the global attributes
-%            .ListDimName: cell listing the names of the array dimensions
-%               .DimValue: array dimension values (Matlab vector with the same length as .ListDimName
 %            .ListVarName: cell listing the names of the variables
-%            .VarDimIndex: cell containing the set of dimension indices (in list .ListDimName) for each variable of .ListVarName
 %           .VarAttribute: cell of structures s containing names and values of variable attributes (s.name=value) for each variable of .ListVarName
 %        .Var1, .Var2....: variables (Matlab arrays) with names listed in .ListVarName
@@ -391 +387 @@
 %A REVOIR, GENERALISER: UTILISER proj_line
 ProjData.NbDim=1;
-ProjData.ListDimName={};%name of dimension 
-ProjData.DimValue=[];%values of dimension (nbre of vectors)
+%ProjData.ListDimName={};%name of dimension 
+%ProjData.DimValue=[];%values of dimension (nbre of vectors)
 ProjData.ListVarName={};
 %ProjData.VarDimIndex={};
@@ -966 +962 @@
      DY=abs(ObjectData.DY);%mesh of interpolation points 
 end
-if  ~isequal(ProjMode,'projection') & (DX==0|DY==0)
+if  ~strcmp(ProjMode,'projection') && (DX==0||DY==0)
         errormsg='DX or DY missing';
         display(errormsg)
@@ -997 +993 @@
 ProjData=proj_heading(FieldData,ObjectData);
 ProjData.NbDim=2;
-ProjData.ListDimName={};%name of dimension 
-ProjData.DimValue=[];%values of dimension (nbre of vectors)
 ProjData.ListVarName={};
 ProjData.VarDimName={};
@@ -1045 +1039 @@
         DimCell={DimCell};%name of dimensions
     end
-    
+
 %case of input fields with unstructured coordinates
     if testX
@@ -1134 +1128 @@
         if isequal(ObjectData.ProjMode,'projection')
             %the list of dimension
-            ProjData.ListDimName=[ProjData.ListDimName FieldData.VarDimName(VarIndex(1))];%add the point index to the list of dimensions
-            ProjData.DimValue=[ProjData.DimValue length(coord_X)];
+            %ProjData.ListDimName=[ProjData.ListDimName FieldData.VarDimName(VarIndex(1))];%add the point index to the list of dimensions
+            %ProjData.DimValue=[ProjData.
+             %length(coord_X)];
             nbvar=0;
             for ivar=VarIndex %transfer variables to the projection plane
@@ -1181 +1176 @@
             testFF=0;
             for ivar=VarIndex
-                VarName=FieldData.ListVarName{ivar}; 
+                VarName=FieldData.ListVarName{ivar};
                 if ~( ivar==ivar_X || ivar==ivar_Y || ivar==ivar_Z || ivar==ivar_F || ivar==ivar_FF || test_anc(ivar)==1)                 
                     ivar_new=ivar_new+1;
@@ -1220 +1215 @@
             end
         end
-%case of fields defined on a structured  grid 
+%case of input fields defined on a structured  grid 
     else  
-        AYName=FieldData.ListVarName{VarType.coord(1)};
-        AXName=FieldData.ListVarName{VarType.coord(2)};
+        AYName=FieldData.ListVarName{VarType.coord(1)};%name of input x coordinate (name preserved on projection)
+        AXName=FieldData.ListVarName{VarType.coord(2)};%name of input y coordinate (name preserved on projection)
         eval(['AX=FieldData.' AXName ';'])
         eval(['AY=FieldData.' AYName ';'])
-        VarName=FieldData.ListVarName{VarIndex(1)};
-        eval(['DimValue=size(FieldData.' VarName ');'])
+        VarName=FieldData.ListVarName{VarIndex(1)};%get the first variable of the cell to get the input matrix dimensions
+        eval(['DimValue=size(FieldData.' VarName ');'])%input matrix dimensions
         ListDimName=FieldData.VarDimName{VarIndex(1)};
         ProjData.ListVarName=[{AYName} {AXName} ProjData.ListVarName]; %TODO: check if it already exists in Projdata (several cells)
@@ -1234 +1229 @@
         for idim=1:length(ListDimName)
             DimName=ListDimName{idim};
-            if isequal(DimName,'rgb')|isequal(DimName,'nb_coord')|isequal(DimName,'nb_coord_i')
+            if strcmp(DimName,'rgb')||strcmp(DimName,'nb_coord')||strcmp(DimName,'nb_coord_i')
                nbcolor=DimValue(idim);
-               DimIndices(idim)=[];
                DimValue(idim)=[];
             end
             if isequal(DimName,'nb_coord_j')% NOTE: CASE OF TENSOR NOT TREATED
-                DimIndices(idim)=[];
                 DimValue(idim)=[];
             end
         end  
         ind_1=find(DimValue==1);
-        DimIndices(ind_1)=[]; %suppress singleton dimensions 
-%         indxy=find(DimVarIndex(DimIndices));%select dimension variables (DimIndices non zero)
-        nb_dim=length(DimIndices);%number of space dimensions
         Coord_z=[];
         Coord_y=[];
         Coord_x=[];   
-    
+        nb_dim=numel(size(DimValue));
         for idim=1:nb_dim %loop on space dimensions
             test_interp(idim)=0;%test for coordiate interpolation (non regular grid), =0 by default
-            test_coord(idim)=0;%test for defined coordinates, =0 by default
-            ivar=DimVarIndex(DimIndices(idim));% index of the variable corresponding to the current dimension
-            if ~isequal(ivar,0)%  a variable corresponds to the current dimension
-                eval(['Coord{idim}=FieldData.' FieldData.ListVarName{ivar} ';']) ;% position for the first index
-                DCoord=diff(Coord{idim});
-                DCoord_min(idim)=min(DCoord);
-                DCoord_max=max(DCoord);
-                test_direct(idim)=DCoord_max>0;% =1 for increasing values, 0 otherwise
-                test_direct_min=DCoord_min(idim)>0;% =1 for increasing values, 0 otherwise
-                if ~isequal(test_direct(idim),test_direct_min)
-                     msgbox_uvmat('ERROR',['non monotonic dimension variable # ' num2str(idim)  ' in proj_field.m'])
+            ivar=VarType.coord(idim);% index of the variable corresponding to the current dimension
+            if ~isequal(ivar,0)%  a variable corresponds to the dimension #idim
+                eval(['Coord{idim}=FieldData.' FieldData.ListVarName{ivar} ';']) ;% coord values for the input field
+                if numel(Coord{idim})==2 %input array defined on a regular grid
+                   DCoord_min(idim)=(Coord{idim}(2)-Coord{idim}(1))/DimValue(idim);
+                else
+                    DCoord=diff(Coord{idim});%array of coordinate derivatives for the input field
+                    DCoord_min(idim)=min(DCoord);
+                    DCoord_max=max(DCoord);
+                %    test_direct(idim)=DCoord_max>0;% =1 for increasing values, 0 otherwise
+                    if ~isequal(DCoord_max,DCoord_min(idim)>0)
+                        msgbox_uvmat('ERROR',['non monotonic dimension variable # ' num2str(idim)  ' in proj_field.m'])
                                 return
-                end               
-                test_interp(idim)=(DCoord_max-DCoord_min(idim))> 0.0001*abs(DCoord_max);% test grid regularity
-                test_coord(idim)=1;
-
-            else  % no variable associated with the first dimension, look for variable  attributes Coord_1, _2 or _3
+                    end               
+                    test_interp(idim)=(DCoord_max-DCoord_min(idim))> 0.0001*abs(DCoord_max);% test grid regularity
+                end
+                test_direct(idim)=(DCoord_min(idim)>0);
+            else  % no variable associated with the  dimension #idim, the coordinate value is set equal to the matrix index by default
                 Coord_i_str=['Coord_' num2str(idim)];
                 DCoord_min(idim)=1;%default
@@ -1277 +1268 @@
             end
         end
-        if nb_dim==2
-            if DY==0
-                DY=abs(DCoord_min(1));
-            end
-            npY=1+round(abs(Coord{1}(end)-Coord{1}(1))/DY);%nbre of points after interpolation 
-            npy=1+round(abs(Coord{1}(end)-Coord{1}(1))/abs(DCoord_min(1)));%nbre of points after possible interpolation on a regular grid
-            if DX==0
-                DX=abs(DCoord_min(2));
-            end
-            npX=1+round(abs(Coord{2}(end)-Coord{2}(1))/DX);%nbre of points after interpol  
-            npx=1+round(abs(Coord{2}(end)-Coord{2}(1))/abs(DCoord_min(2)));%nbre of points after possible interpolation on a regular grid 
-            Coord_y=linspace(Coord{1}(1),Coord{1}(end),npY);
-            test_direct_y=test_direct(1);
-            Coord_x=linspace(Coord{2}(1),Coord{2}(end),npX);
-            test_direct_x=test_direct(2);
-            DAX=DCoord_min(2);
-            DAY=DCoord_min(1);
-        elseif nb_dim==3
+        if DY==0
+            DY=abs(DCoord_min(nb_dim-1));
+        end
+        npY=1+round(abs(Coord{nb_dim-1}(end)-Coord{nb_dim-1}(1))/DY);%nbre of points after interpol 
+        if DX==0
+            DX=abs(DCoord_min(nb_dim));
+        end
+        npX=1+round(abs(Coord{nb_dim}(end)-Coord{nb_dim}(1))/DX);%nbre of points after interpol 
+        for idim=[1:nb_dim]
+            if test_interp(idim)
+                DimValue(idim)=1+round(abs(Coord{idim}(end)-Coord{idim}(1))/abs(DCoord_min(idim)));%nbre of points after possible interpolation on a regular gri
+            end
+        end       
+        Coord_y=linspace(Coord{nb_dim-1}(1),Coord{nb_dim-1}(end),npY);
+        test_direct_y=test_direct(nb_dim-1);
+        Coord_x=linspace(Coord{nb_dim}(1),Coord{nb_dim}(end),npX);
+        test_direct_x=test_direct(nb_dim);
+        DAX=DCoord_min(nb_dim);
+        DAY=DCoord_min(nb_dim-1);
+        if nb_dim==3
             DZ=abs(DCoord_min(1));
-            npz=1+round(abs(Coord{1}(end)-Coord{1}(1))/DZ);%nbre of points after interpolation
-            if DY==0
-                DY=abs(DCoord_min(2));
-            end
-            npY=1+round(abs(Coord{2}(end)-Coord{2}(1))/DY);%nbre of points after interpol 
-            npy=1+round(abs(Coord{2}(end)-Coord{2}(1))/abs(DCoord_min(2)));%nbre of points before interpol 
-            if DX==0
-                DX=abs(DCoord_min(3));
-            end
-            npX=1+round(abs(Coord{3}(end)-Coord{3}(1))/DX);%nbre of points after interpol
-            npx=1+round(abs(Coord{3}(end)-Coord{3}(1))/abs(DCoord_min(3)));%nbre of points before interpol 
             Coord_z=linspace(Coord{1}(1),Coord{1}(end),npz);
             test_direct_z=test_direct(1);
-            Coord_y=linspace(Coord{2}(1),Coord{2}(end),npY);
-            test_direct_y=test_direct(2);
-            Coord_x=linspace(Coord{3}(1),Coord{3}(end),npX);
-            test_direct_x=test_direct(3);
         end  
         minAX=min(Coord_x);
@@ -1334 +1312 @@
             YMin=min(ycor_new);
         end
-        DXinit=(maxAX-minAX)/(npx-1);
-        DYinit=(maxAY-minAY)/(npy-1);
+        DXinit=(maxAX-minAX)/(DimValue(2)-1);
+        DYinit=(maxAY-minAY)/(DimValue(1)-1);
         if DX==0
             DX=DXinit;
@@ -1343 +1321 @@
         end
         npX=floor((XMax-XMin)/DX+1);
-        npY=floor((YMax-YMin)/DY+1);    
+        npY=floor((YMax-YMin)/DY+1);   
         if test_direct_y
             coord_y_proj=linspace(YMin,YMax,npY);%abscissa of the new pixels along the line
@@ -1381 +1359 @@
             min_ind1=max(min_ind1,1);% deals with margin (bound lower than the first index)
             min_ind2=max(min_ind2,1);
-            max_ind1=min(max_ind1,npy);
-            max_ind2=min(max_ind2,npx);
+            max_ind1=min(max_ind1,DimValue(1));
+            max_ind2=min(max_ind2,DimValue(2));
             for ivar=VarIndex
                 VarName=FieldData.ListVarName{ivar}; 
                 ProjData.ListVarName=[ProjData.ListVarName VarName];
-                ProjData.VarDimIndex=[ProjData.VarDimIndex [nb_dim-1 nb_dim]];
+                ProjData.VarDimName=[ProjData.VarDimName {DimCell}];
                 if length(FieldData.VarAttribute)>=ivar
                     ProjData.VarAttribute{length(ProjData.ListVarName)}=FieldData.VarAttribute{ivar};
@@ -1392 +1370 @@
                 eval(['ProjData.' VarName '=FieldData.' VarName '(min_ind1:max_ind1,min_ind2:max_ind2) ;']);
             end         
-        else
-        % case with rotation and/or interpolation
+        else       % case with rotation and/or interpolation
             if isempty(Coord_z) %2D case
                 [X,Y]=meshgrid(coord_x_proj,coord_y_proj);%grid in the new coordinates
@@ -1402 +1379 @@
                 XIMA=reshape(round(XIMA),1,npX*npY);%indices reorganized in 'line'
                 YIMA=reshape(round(YIMA),1,npX*npY);
-                flagin=XIMA>=1 & XIMA<=npx & YIMA >=1 & YIMA<=npy;%flagin=1 inside the original image 
+                flagin=XIMA>=1 & XIMA<=DimValue(2) & YIMA >=1 & YIMA<=DimValue(1);%flagin=1 inside the original image 
                 if isequal(ObjectData.ProjMode,'filter')
                     npx_filter=ceil(abs(DX/DAX));
@@ -1411 +1388 @@
                     test_filter=0;
                 end
+                eval(['ProjData.' AYName '=[coord_y_proj(1) coord_y_proj(end)];']) %record the new (projected ) y coordinates
+                eval(['ProjData.' AXName '=[coord_x_proj(1) coord_x_proj(end)];']) %record the new (projected ) x coordinates
                 for ivar=VarIndex
-                    VarName=FieldData.ListVarName{ivar} ; 
-                    if test_interp(1) | test_interp(2)%interpolate on a regular grid         
-                          eval(['FieldData.' VarName '=interp2(Coord{2},Coord{1},FieldData.' VarName ',Coord_x,Coord_y'');']) %TO TEST
+                    VarName=FieldData.ListVarName{ivar};
+                    if test_interp(1) || test_interp(2)%interpolate on a regular grid        
+                          eval(['ProjData.' VarName '=interp2(Coord{2},Coord{1},FieldData.' VarName ',Coord_x,Coord_y'');']) %TO TEST
                     end
                     %filter the field (image) if option 'filter' is used
                     if test_filter  
                          Aclass=class(FieldData.A);
-                         eval(['FieldData.' VarName '=filter2(Mfilter,FieldData.' VarName ',''valid'');'])
+                         eval(['ProjData.' VarName '=filter2(Mfilter,FieldData.' VarName ',''valid'');'])
                          if ~isequal(Aclass,'double')
-                             eval(['FieldData.' VarName '=' Aclass '(FieldData.' VarName ');'])%revert to integer values
+                             eval(['ProjData.' VarName '=' Aclass '(FieldData.' VarName ');'])%revert to integer values
                          end
                     end
-                    eval(['vec_A=reshape(FieldData.' VarName ',npx*npy,nbcolor);'])%put the original image in line              
+                    eval(['vec_A=reshape(FieldData.' VarName ',[],nbcolor);'])%put the original image in line              
                     ind_in=find(flagin);
                     ind_out=find(~flagin);
-                    ICOMB=(XIMA-1)*npy+YIMA;
+                    ICOMB=(XIMA-1)*DimValue(1)+YIMA;
                     ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A
                     vec_B(ind_in,[1:nbcolor])=vec_A(ICOMB,:); 
@@ -1433 +1412 @@
                         vec_B(ind_out,icolor)=zeros(size(ind_out));
                     end
-                    ProjData.ListVarName=[ProjData.ListVarName VarName];                 
-                    if length(FieldData.VarAttribute)>=ivar
+                    ProjData.ListVarName=[ProjData.ListVarName VarName];
+                    ProjData.VarDimName=[ProjData.VarDimName {DimCell}];
+                    if isfield(FieldData,'VarAttribute')&&length(FieldData.VarAttribute)>=ivar
                         ProjData.VarAttribute{length(ProjData.ListVarName)+nbcoord}=FieldData.VarAttribute{ivar};
                     end     
@@ -1441 +1421 @@
                 ProjData.FF=reshape(~flagin,npY,npX);%false flag A FAIRE: tenir compte d'un flga antérieur  
                 ProjData.ListVarName=[ProjData.ListVarName 'FF'];
+                ProjData.VarDimName=[ProjData.VarDimName {DimCell}];
                 ProjData.VarAttribute{length(ProjData.ListVarName)}.Role='errorflag';
             else %3D case
@@ -1448 +1429 @@
                     iz=iz_sup(1);
                     if iz>=1 & iz<=npz
-                        ProjData.ListDimName=[ProjData.ListDimName ListDimName(2:end)];
-                        ProjData.DimValue=[ProjData.DimValue npY npX];
+                        %ProjData.ListDimName=[ProjData.ListDimName ListDimName(2:end)];
+                        %ProjData.DimValue=[ProjData.DimValue npY npX];
                         for ivar=VarIndex
                             VarName=FieldData.ListVarName{ivar}; 
@@ -1582 +1563 @@
 ProjData=proj_heading(FieldData,ObjectData);
 ProjData.NbDim=3;
-ProjData.ListDimName={};%name of dimension 
-ProjData.DimValue=[];%values of dimension (nbre of vectors)
+%ProjData.ListDimName={};%name of dimension 
+%ProjData.DimValue=[];%values of dimension (nbre of vectors)
 ProjData.ListVarName={};
 ProjData.VarDimName={};
@@ -1717 +1698 @@
         if isequal(ObjectData.ProjMode,'projection')
             %the list of dimension
-            ProjData.ListDimName=[ProjData.ListDimName FieldData.VarDimName(VarIndex(1))];%add the point index to the list of dimensions
-            ProjData.DimValue=[ProjData.DimValue length(coord_X)];
+            %ProjData.ListDimName=[ProjData.ListDimName FieldData.VarDimName(VarIndex(1))];%add the point index to the list of dimensions
+            %ProjData.DimValue=[ProjData.DimValue length(coord_X)];
             nbvar=0;
             for ivar=VarIndex %transfer variables to the projection plane
@@ -1843 +1824 @@
         for idim=1:nb_dim %loop on space dimensions
             test_interp(idim)=0;%test for coordiate interpolation (non regular grid), =0 by default
-            test_coord(idim)=0;%test for defined coordinates, =0 by default
             ivar=DimVarIndex(DimIndices(idim));% index of the variable corresponding to the current dimension
             if ~isequal(ivar,0)%  a variable corresponds to the current dimension
@@ -1857 +1837 @@
                 end               
                 test_interp(idim)=(DCoord_max-DCoord_min(idim))> 0.0001*abs(DCoord_max);% test grid regularity
-                test_coord(idim)=1;
-
             else  % no variable associated with the first dimension, look for variable  attributes Coord_1, _2 or _3
                 Coord_i_str=['Coord_' num2str(idim)];
@@ -1871 +1849 @@
             end
             npY=1+round(abs(Coord{1}(end)-Coord{1}(1))/DY);%nbre of points after interpolation 
-            npy=1+round(abs(Coord{1}(end)-Coord{1}(1))/abs(DCoord_min(1)));%nbre of points after possible interpolation on a regular grid
+            DimValue(1)=1+round(abs(Coord{1}(end)-Coord{1}(1))/abs(DCoord_min(1)));%nbre of points after possible interpolation on a regular grid
             if DX==0
                 DX=abs(DCoord_min(2));
             end
             npX=1+round(abs(Coord{2}(end)-Coord{2}(1))/DX);%nbre of points after interpol  
-            npx=1+round(abs(Coord{2}(end)-Coord{2}(1))/abs(DCoord_min(2)));%nbre of points after possible interpolation on a regular grid 
+            DimValue(2)=1+round(abs(Coord{2}(end)-Coord{2}(1))/abs(DCoord_min(2)));%nbre of points after possible interpolation on a regular grid 
             Coord_y=linspace(Coord{1}(1),Coord{1}(end),npY);
             test_direct_y=test_direct(1);
@@ -1890 +1868 @@
             end
             npY=1+round(abs(Coord{2}(end)-Coord{2}(1))/DY);%nbre of points after interpol 
-            npy=1+round(abs(Coord{2}(end)-Coord{2}(1))/abs(DCoord_min(2)));%nbre of points before interpol 
+            DimValue(1)=1+round(abs(Coord{2}(end)-Coord{2}(1))/abs(DCoord_min(2)));%nbre of points before interpol 
             if DX==0
                 DX=abs(DCoord_min(3));
             end
             npX=1+round(abs(Coord{3}(end)-Coord{3}(1))/DX);%nbre of points after interpol
-            npx=1+round(abs(Coord{3}(end)-Coord{3}(1))/abs(DCoord_min(3)));%nbre of points before interpol 
+            DimValue(2)=1+round(abs(Coord{3}(end)-Coord{3}(1))/abs(DCoord_min(3)));%nbre of points before interpol 
             Coord_z=linspace(Coord{1}(1),Coord{1}(end),npz);
             test_direct_z=test_direct(1);
@@ -1923 +1901 @@
             YMin=min(ycor_new);
         end
-        DXinit=(maxAX-minAX)/(npx-1);
-        DYinit=(maxAY-minAY)/(npy-1);
+        DXinit=(maxAX-minAX)/(DimValue(2)-1);
+        DYinit=(maxAY-minAY)/(DimValue(1)-1);
         if DX==0
             DX=DXinit;
@@ -1970 +1948 @@
             min_ind1=max(min_ind1,1);% deals with margin (bound lower than the first index)
             min_ind2=max(min_ind2,1);
-            max_ind1=min(max_ind1,npy);
-            max_ind2=min(max_ind2,npx);
+            max_ind1=min(max_ind1,DimValue(1));
+            max_ind2=min(max_ind2,DimValue(2));
             for ivar=VarIndex
                 VarName=FieldData.ListVarName{ivar}; 
                 ProjData.ListVarName=[ProjData.ListVarName VarName];
-                ProjData.VarDimIndex=[ProjData.VarDimIndex [nb_dim-1 nb_dim]];
+                %ProjData.VarDimIndex=[ProjData.VarDimIndex [nb_dim-1 nb_dim]];
                 if length(FieldData.VarAttribute)>=ivar
                     ProjData.VarAttribute{length(ProjData.ListVarName)}=FieldData.VarAttribute{ivar};
@@ -1991 +1969 @@
                 XIMA=reshape(round(XIMA),1,npX*npY);%indices reorganized in 'line'
                 YIMA=reshape(round(YIMA),1,npX*npY);
-                flagin=XIMA>=1 & XIMA<=npx & YIMA >=1 & YIMA<=npy;%flagin=1 inside the original image 
+                flagin=XIMA>=1 & XIMA<=DimValue(2) & YIMA >=1 & YIMA<=DimValue(1);%flagin=1 inside the original image 
                 if isequal(ObjectData.ProjMode,'filter')
                     npx_filter=ceil(abs(DX/DAX));
@@ -2013 +1991 @@
                          end
                     end
-                    eval(['vec_A=reshape(FieldData.' VarName ',npx*npy,nbcolor);'])%put the original image in line              
+                    eval(['vec_A=reshape(FieldData.' VarName ',[],nbcolor);'])%put the original image in line              
                     ind_in=find(flagin);
                     ind_out=find(~flagin);
-                    ICOMB=(XIMA-1)*npy+YIMA;
+                    ICOMB=(XIMA-1)*DimValue(1)+YIMA;
                     ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A
                     vec_B(ind_in,[1:nbcolor])=vec_A(ICOMB,:); 
@@ -2037 +2015 @@
                     iz=iz_sup(1);
                     if iz>=1 & iz<=npz
-                        ProjData.ListDimName=[ProjData.ListDimName ListDimName(2:end)];
-                        ProjData.DimValue=[ProjData.DimValue npY npX];
+                        %ProjData.ListDimName=[ProjData.ListDimName ListDimName(2:end)];
+                        %ProjData.DimValue=[ProjData.DimValue npY npX];
                         for ivar=VarIndex
                             VarName=FieldData.ListVarName{ivar}; 

trunk/src/px_XYZ.m

@@ +1 @@
+%'px_XYZ': transform physical to image coordinates. 
+%------------------------------------------------------------------------
+%[X,Y]=px_XYZ(Calib,Xphys,Yphys,Zphys)
+%------------------------------------------------------------------------           
+% OUTPUT:
+% [X,Y]: image coordinates(in pixels)
+%------------------------------------------------------------------------
+% INPUT:
+% Calib: structure containing calibration parameters
+% Xphys,Yphys,Zphys; vectors of physical coordinates for a set of points
 
 function [X,Y]=px_XYZ(Calib,Xphys,Yphys,Zphys)

trunk/src/series/merge_proj.m

@@ -59 +59 @@
 nbview=length(Series.RootFile);%number of views (file series to merge)
 nbfield=size(num_i1{1},1)*size(num_i1{1},2);%number of fields in the time series
-transform=Series.CoordType; %  field transform function
+% transform=Series.CoordType; %  field transform function
 hhh=which('mmreader');
 for iview=1:nbview
@@ -72 +72 @@
 
 %Calibration data and timing: read the ImaDoc files
-mode=''; %default
+% mode=''; %default
 timecell={};
 itime=0;
@@ -152 +152 @@
 
 % Field and velocity type (the same for all views)
+FieldName='';
+if strcmp(get(hseries.FieldMenu,'Visible'),'on')
 Field_str=get(hseries.FieldMenu,'String');
 val=get(hseries.FieldMenu,'Value');
@@ -163 +165 @@
     hget_field=findobj(allchild(0),'Name','get_field');%find the get_field... GUI
     SubField=get_field('read_get_field',hObject,eventdata,hget_field); %read the names of the variables to plot in the get_field GUI
+end
 end
 %detect whether all the files are 'images' or 'netcdf'
@@ -238 +241 @@
     if isequal(stopstate,'queue')% enable STOP command from the 'series' interface
          update_waitbar(hseries.waitbar,WaitbarPos,ifile/nbfield)
-         Amerge=0;
+%          Amerge=0;
          
          %----------LOOP ON VIEWS----------------------
@@ -258 +261 @@
                 end % TODO: introduce ListVarName
                 npxy=size(Field{iview}.A);
+                Field{iview}.ListVarName={'AX','AY','A'};
+                Field{iview}.VarDimName={'AX','AY',{'AY','AX'}};
                 Field{iview}.AX=[0.5 npxy(2)-0.5]; % coordinates of the first and last pixel centers
                 Field{iview}.AY=[npxy(1)-0.5 0.5];
                 Field{iview}.CoordType='px'; 
                 Field{iview}.AName='image';
+                timeread(iview)=0;
             else
                 if testcivx
@@ -287 +293 @@
             end
             if testcivx
-                    Field{iview}=calc_field(FieldName,Field{iview});
+                Field{iview}=calc_field(FieldName,Field{iview});
             end
 
@@ -393 +399 @@
     return
 end
-for iview=1:nbview
-    if ~isequal(MergeData.ListDimName,Data{iview}.ListDimName)
-        error=1;
-    end
-    if ~isequal(MergeData.ListVarName,Data{iview}.ListVarName)
-        error=1;
-    end
-end
-if error
-    MergeData.Txt='ERROR: attempt at merging fields of incompatible type';
-    return
-end
-
-
-
-
+% for iview=1:nbview
+%     if ~isequal(MergeData.ListDimName,Data{iview}.ListDimName)
+%         error=1;
+%     end
+%     if ~isequal(MergeData.ListVarName,Data{iview}.ListVarName)
+%         error=1;
+%     end
+% end
+% if error
+%     MergeData.Txt='ERROR: attempt at merging fields of incompatible type';
+%     return
+% end
+% 
 %group the variables (fields of 'FieldData') in cells of variables with the same dimensions
 [CellVarIndex,NbDim,VarTypeCell]=find_field_indices(Data{1});
@@ -414 +417 @@
 % CellVarIndex=cells of variable index arrays
 ivar_new=0; % index of the current variable in the projected field
-icoord=0;
 for icell=1:length(CellVarIndex)
     if NbDim(icell)==1

trunk/src/sub_field.m

@@ -102 +102 @@
    %check coincidence in positions
    %unstructured coordinates
-   if testX
+   if testX      
+       if testU % vector fields
+            U_1_Name=Field_1.ListVarName{VarType_1.vector_x};
+            V_1_Name=Field_1.ListVarName{VarType_1.vector_y};
+            eval(['vec_U_1=Field_1.' U_1_Name ';']) 
+            eval(['vec_V_1=Field_1.' V_1_Name ';'])
+            nbpoints_x_1=size(vec_U_1,2);
+            nbpoints_y_1=size(vec_U_1,1);
+            vec_U_1=reshape(vec_U_1,nbpoints_x_1*nbpoints_y_1,1);
+            vec_V_1=reshape(vec_V_1,nbpoints_x_1*nbpoints_y_1,1);
+            if testfalse_1
+                vec_U_1=vec_U_1(indsel);
+                vec_V_1=vec_V_1(indsel);
+            end            
+       else %(~testU && ~testU_1)
+           A_1_Name=Field_1.ListVarName{ivar_C_1};
+           eval(['vec_A_1=Field_1.' A_1_Name ';'])   
+           nbpoints_x_1=size(vec_A_1,2);
+           nbpoints_y_1=size(vec_A_1,1);%TODO: use a faster interpolation method for a regular grid (size(x)=2)
+           vec_A_1=reshape(vec_A_1,nbpoints_x_1*nbpoints_y_1,1);
+           if testfalse_1
+                vec_A_1=vec_A_1(indsel);
+           end
+       end
        XName=Field.ListVarName{VarType.coord_x};
        YName=Field.ListVarName{VarType.coord_y};
@@ -108 +131 @@
        eval(['vec_Y=Field.' YName ';'])
        nbpoints=numel(vec_X);
-       vec_X=reshape(vec_X,1,nbpoints);
-       vec_Y=reshape(vec_Y,1,nbpoints);
+       vec_X=reshape(vec_X,nbpoints,1);
+       vec_Y=reshape(vec_Y,nbpoints,1);
        if testX_1 %unstructured coordinates for the second field
             X_1_Name=Field_1.ListVarName{VarType_1.coord_x};
@@ -115 +138 @@
             eval(['vec_X_1=Field_1.' X_1_Name ';']) 
             eval(['vec_Y_1=Field_1.' Y_1_Name ';'])
-            nbpoints_1=numel(vec_X_1);
+
        else   %structured coordinates for the second field
            y_1_Name=Field_1.ListVarName{VarType_1.coord(1)};
            x_1_Name=Field_1.ListVarName{VarType_1.coord(2)};
            eval(['y_1=Field_1.' y_1_Name ';']) 
-           eval(['x_1=Field_1.' x_1_Name ';'])  
-           npxy(1)=numel(y_1);
-           npxy(2)=numel(x_1);
-           nbpoints_1=npxy(1)*npxy(2);
+           eval(['x_1=Field_1.' x_1_Name ';'])
+           if isequal(numel(x_1),2)  
+               x_1=linspace(x_1(1),x_1(2),nbpoints_x_1);
+           end
+           if isequal(numel(y_1),2)  
+               y_1=linspace(y_1(1),y_1(2),nbpoints_y_1);
+           end
            [vec_X_1,vec_Y_1]=meshgrid(x_1,y_1);
        end
-       vec_X_1=reshape(vec_X_1,1,nbpoints_1);
-       vec_Y_1=reshape(vec_Y_1,1,nbpoints_1);
+       vec_X_1=reshape(vec_X_1,nbpoints_x_1*nbpoints_y_1,1);
+       vec_Y_1=reshape(vec_Y_1,nbpoints_x_1*nbpoints_y_1,1);
        if testfalse_1
            FFName_1=Field_1.ListVarName{VarType_1.errorflag};          
            eval(['vec_FF_1=Field_1.' FFName_1 ';']) 
-           vec_FF_1=reshape(vec_FF_1,1,nbpoints_1);
+           vec_FF_1=reshape(vec_FF_1,nbpoints_1,1);
            indsel=find(~vec_FF_1);
            vec_X_1=vec_X_1(indsel);
            vec_Y_1=vec_Y_1(indsel);
-       end
-       if testU % vector fields
-            U_1_Name=Field_1.ListVarName{VarType_1.vector_x};
-            V_1_Name=Field_1.ListVarName{VarType_1.vector_y};
-            eval(['vec_U_1=Field_1.' U_1_Name ';']) 
-            eval(['vec_V_1=Field_1.' V_1_Name ';'])
-            vec_U_1=reshape(vec_U_1,1,nbpoints_1);
-            vec_V_1=reshape(vec_V_1,1,nbpoints_1);
-            if testfalse_1
-                vec_U_1=vec_U_1(indsel);
-                vec_V_1=vec_V_1(indsel);
-            end            
-       else
-           A_1_Name=Field_1.ListVarName{ivar_C_1};
-           eval(['vec_A_1=Field_1.' A_1_Name ';'])
-           vec_A_1=reshape(vec_A_1,1,nbpoints_1);
-           if testfalse_1
-                vec_A_1=vec_A_1(indsel);
-           end
        end
        if ~isequal(vec_X_1,vec_X) && ~isequal(vec_Y_1,vec_Y) % if the unstructured positions are not the same
@@ -160 +167 @@
                vec_V_1=griddata_uvmat(vec_X_1,vec_Y_1,vec_V_1,vec_X,vec_Y);  %interpolate vectors in the second field   
            else
-               vec_A_1=griddata_uvmat(vec_X_1,vec_Y_1,vec_A_1,vec_X,vec_Y);  %interpolate vectors in the second field
+               vec_A_1=griddata_uvmat(vec_X_1,vec_Y_1,double(vec_A_1),vec_X,vec_Y);  %interpolate vectors in the second field
            end
        end 
@@ -168 +175 @@
            eval(['vec_U=Field.' UName ';']) 
            eval(['vec_V=Field.' VName ';'])       
-           vec_U=reshape(vec_U,1,numel(vec_U));
-           vec_V=reshape(vec_V,1,numel(vec_V));
+           vec_U=reshape(vec_U,numel(vec_U),1);
+           vec_V=reshape(vec_V,numel(vec_V),1);
            eval(['SubData.' UName '=vec_U-vec_U_1;'])
            eval(['SubData.' VName '=vec_V-vec_V_1;'])
        else
            AName=Field.ListVarName{ivar_C};
+           size(Field.vort)
            eval(['SubData.' AName '=Field.' AName '-vec_A_1;'])
        end