Changeset 109 for trunk/src


Ignore:
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:
9 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/create_grid.m

    r89 r109  
    136136    zarray=T.z_0*ones(size(yarray));
    137137    varargout{1}=[xarray yarray zarray];
     138    varargout{2}=T;
    138139end
    139 varargout{2}=T;
     140
    140141% The figure can be deleted now
    141142delete(handles.figure1);
  • trunk/src/geometry_calib.m

    r108 r109  
    33% function varargout = geometry_calib(varargin)
    44%
    5 %A%AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
     5%A%AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
    66%  Copyright Joel Sommeria, 2008, LEGI / CNRS-UJF-INPG, sommeria@coriolis-legi.org.
    77%AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
     
    4343% Edit the above text to modify the response to help geometry_calib
    4444
    45 % Last Modified by GUIDE v2.5 25-Mar-2010 19:10:05
     45% Last Modified by GUIDE v2.5 03-Oct-2010 09:34:12
    4646
    4747% Begin initialization code - DO NOT edit
     
    8787    set(hObject,'Position',pos_gui);
    8888end
     89
     90%set menu of calibration options
     91%set(handles.calib_type,'String',{'rescale';'linear';'perspective';'normal';'tsai';'bouguet';'extrinsic'})
     92set(handles.calib_type,'String',{'rescale';'linear';'quadr';'3D_linear';'3D_quadr';'3D_extrinsic'})
    8993inputxml='';
    9094if exist('inputfile','var')& ~isempty(inputfile)
    91 %     [Path,Name,ext]=fileparts(inputfile);
    92 %     form=imformats(ext([2:end]));
    93 %     if ~isempty(form)% if the input file is an image
    94         struct.XmlInputfile=inputfile;
    95         set(hObject,'UserData',struct)
    96         [Pathsub,RootFile,field_count,str2,str_a,str_b,ext,nom_type,subdir]=name2display(inputfile);
    97         inputxml=[fullfile(Pathsub,RootFile) '.xml'];
    98 %     end   
     95    struct.XmlInputfile=inputfile;
     96    set(hObject,'UserData',struct)
     97    [Pathsub,RootFile,field_count,str2,str_a,str_b,ext,nom_type,subdir]=name2display(inputfile);
     98    inputxml=[fullfile(Pathsub,RootFile) '.xml'];
    9999end
    100100set(handles.ListCoord,'String',{'......'})
     
    102102    loadfile(handles,inputxml)% load the point coordiantes existing in the xml file
    103103end
    104 
    105104set(handles.ListCoord,'KeyPressFcn',{@key_press_fcn,handles})%set keyboard action function
    106105
     
    113112varargout{1} = handles.output;
    114113varargout{2}=handles;
    115 
    116 %------------
    117 function Phi_Callback(hObject, eventdata, handles)
    118 
    119 %------------------------------------------------------------------------
    120 %read input xml file and update the edit boxes
    121 function loadfile(handles,fileinput)
    122 %------------------------------------------------------------------------
    123 %read the input xml file
    124 t=xmltree(fileinput);
    125 s=convert(t);%convert to matlab structure
    126 %read data currently displayed on the interface
    127 PointCoord=[];
    128 Coord_cell=get(handles.ListCoord,'String');
    129 data=read_geometry_calib(Coord_cell);
    130 %data=read_geometry_calib(handles);
    131 Coord=[]; %default
    132 if isfield(data,'Coord')
    133     Coord=data.Coord;
    134 end
    135 TabChar_0=get(handles.ListCoord,'String');
    136 nbcoord_0=size(TabChar_0,1);
    137 if isequal(get(handles.edit_append,'Value'),2) %edit mode  A REVOIR
    138     val=get(handles.ListCoord,'Value')-1;
    139 else
    140    val=length(TabChar_0);
    141 end
    142 nbcoord=0;
    143 
    144 %case of calibration (ImaDoc) input file
    145 % hcalib=get(handles.calib_type,'parent');
    146 CalibData=get(handles.geometry_calib,'UserData');
    147 CalibData.XmlInput=fileinput;
    148 if isfield(s,'Heading')
    149     CalibData.Heading=s.Heading;
    150 end
    151 
    152 set(handles.geometry_calib,'UserData',CalibData);%store the heading in the interface 'UserData'
    153 if isfield(s,'GeometryCalib')
    154     Calib=s.GeometryCalib;
    155     if isfield(Calib,'CalibrationType')
    156         CalibrationType=Calib.CalibrationType;
    157         switch CalibrationType
    158             case 'linear'
    159                 set(handles.calib_type,'Value',2)
    160             case 'tsai'
    161                 set(handles.calib_type,'Value',3)
    162             case 'tsai_matlab'
    163                 set(handles.calib_type,'Value',4)
    164         end
    165     end
    166     if isfield(Calib,'SourceCalib')
    167         if isfield(Calib.SourceCalib,'PointCoord')
    168             PointCoord=Calib.SourceCalib.PointCoord;
    169         end
    170     end
    171     nbcoord=length(PointCoord);
    172     if ~isfield(Calib,'ErrorRms')&~isfield(Calib,'ErrorMax') %old convention of Gauthier (cord in mm)
    173         for i=1:length(PointCoord)
    174           line=str2num(PointCoord{i});
    175           Coord(i+val,4:5)=line(4:5);%px x
    176           Coord(i+val,1:3)=line(1:3)/10;%phys x
    177         end
    178     else
    179         for i=1:length(PointCoord)
    180           line=str2num(PointCoord{i});
    181           Coord(i,4:5)=line(4:5);%px x
    182           Coord(i,1:3)=line(1:3);%phys x
    183        end
    184     end
    185 end
    186 %case of xml files of points
    187 if isfield(s,'Coord')
    188     PointCoord=s.Coord;
    189     nbcoord=length(PointCoord);
    190      %case of image coordinates
    191     if isfield(s,'CoordType')& isequal(s.CoordType,'px')
    192         for i=1:nbcoord
    193            line=str2num(PointCoord{i});
    194            Coord(i+val,4:5)=line(1:2);
    195         end
    196      %case of  physical coordinates
    197     else
    198         for i=1:nbcoord
    199            line=str2num(PointCoord{i});
    200            Coord(i+val,1:3)=line(1:3);
    201            nbcolumn=size(Coord,2);
    202            if nbcolumn<5
    203                Coord(i+val,nbcolumn+1:5)=zeros(1,5-nbcolumn);
    204            end
    205         end
    206      end
    207 end
    208 CoordCell={};
    209 for iline=1:size(Coord,1)
    210     for j=1:5
    211         CoordCell{iline,j}=num2str(Coord(iline,j),4);
    212     end
    213 end
    214 % CoordCell=[CoordCell;{' ',' ',' ',' ',' '}];
    215 Tabchar=cell2tab(CoordCell,'    |    ');%transform cells into table ready for display
    216 Tabchar=[Tabchar;{'......'}];
    217 set(handles.ListCoord,'Value',1)
    218 set(handles.ListCoord,'String',Tabchar)
    219 MenuPlot_Callback(handles.geometry_calib, [], handles)
    220 if isempty(Coord)
    221     set(handles.edit_append,'Value',1)
    222     set(handles.edit_append,'BackgroundColor',[1 1 0])
    223 else
    224     set(handles.edit_append,'Value',0)
    225     set(handles.edit_append,'BackgroundColor',[0.7 0.7 0.7])
    226 end
    227114%
    228115%------------------------------------------------------------------------
     
    252139function APPLY_Callback(hObject, eventdata, handles)
    253140%------------------------------------------------------------------------
    254 calib_cell=get(handles.calib_type,'String');
    255 val=get(handles.calib_type,'Value');
    256 calib_type=calib_cell{val};
     141
     142%read the current calibration points
    257143Coord_cell=get(handles.ListCoord,'String');
    258144Object=read_geometry_calib(Coord_cell);
    259 X=Object.Coord(:,1);
    260 Y=Object.Coord(:,2);
    261 Z=Object.Coord(:,3);
    262 if isequal(calib_type,'rescale')
    263     GeometryCalib=calib_rescale(Object.Coord);
    264     Z=0;%Z not taken into account
    265 elseif isequal(calib_type,'linear')
    266     GeometryCalib=calib_linear(Object.Coord);
    267     Z=0; %Z not taken into account
    268 elseif isequal(calib_type,'tsai_cpp')
    269     GeometryCalib=calib_tsai(Object.Coord);
    270 elseif isequal(calib_type,'tsai_matlab')
    271     GeometryCalib=calib_tsai2(Object.Coord);
    272 end
     145Coord=Object.Coord;
     146
     147% apply the calibration, whose type is selected in  handles.calib_type
     148if ~isempty(Coord)
     149    calib_cell=get(handles.calib_type,'String');
     150    val=get(handles.calib_type,'Value');
     151    GeometryCalib=feval(['calib_' calib_cell{val}],Coord,handles);
     152else
     153    msgbox_uvmat('ERROR','No calibration points, abort')
     154    return
     155end
    273156
    274157%check error
     
    298181    Calib.R=GeometryCalib.R;
    299182end
    300 x_ima=Object.Coord(:,4);
    301 y_ima=Object.Coord(:,5);
    302 [Xpoints,Ypoints]=px_XYZ(Calib,X,Y,Z);
    303 GeometryCalib.ErrorRms(1)=sqrt(mean((Xpoints-x_ima).*(Xpoints-x_ima)));
    304 [GeometryCalib.ErrorMax(1),index(1)]=max(abs(Xpoints-x_ima));
    305 GeometryCalib.ErrorRms(2)=sqrt(mean((Ypoints-y_ima).*(Ypoints-y_ima)));
    306 [GeometryCalib.ErrorMax(2),index(2)]=max(abs(Ypoints-y_ima));
    307 [EM,ind_dim]=max(GeometryCalib.ErrorMax);
    308 index=index(ind_dim);
    309 
     183if ~isempty(Coord)
     184    X=Coord(:,1);
     185    Y=Coord(:,2);
     186    Z=Coord(:,3);
     187    x_ima=Coord(:,4);
     188    y_ima=Coord(:,5);
     189    [Xpoints,Ypoints]=px_XYZ(Calib,X,Y,Z);
     190    GeometryCalib.ErrorRms(1)=sqrt(mean((Xpoints-x_ima).*(Xpoints-x_ima)));
     191    [GeometryCalib.ErrorMax(1),index(1)]=max(abs(Xpoints-x_ima));
     192    GeometryCalib.ErrorRms(2)=sqrt(mean((Ypoints-y_ima).*(Ypoints-y_ima)));
     193    [GeometryCalib.ErrorMax(2),index(2)]=max(abs(Ypoints-y_ima));
     194    [EM,ind_dim]=max(GeometryCalib.ErrorMax);
     195    index=index(ind_dim);
     196    if isequal(max(Z),min(Z))%Z constant
     197        GeometryCalib.NbSlice=1;
     198        GeometryCalib.SliceCoord=[0 0 Z(1)];
     199    end
     200end
    310201unitlist=get(handles.CoordUnit,'String');
    311202unit=unitlist{get(handles.CoordUnit,'value')};
    312203GeometryCalib.CoordUnit=unit;
    313 GeometryCalib.SourceCalib.PointCoord=Object.Coord;
     204GeometryCalib.SourceCalib.PointCoord=Coord;
     205
     206%display calibration intrinsic parameters
     207k=0;
     208if isfield(GeometryCalib,'kappa1')
     209    k=GeometryCalib.kappa1 * GeometryCalib.focal*GeometryCalib.focal;
     210end
     211sx=1;dpx_dpy=[1 1];
     212if isfield(GeometryCalib,'sx')
     213    sx=GeometryCalib.sx;
     214end
     215if isfield(GeometryCalib,'dpx_dpy')
     216    dpx_dpy=GeometryCalib.dpx_dpy;
     217end
     218Cx_Cy=[0 0];
     219if isfield(GeometryCalib,'Cx_Cy')
     220    Cx_Cy=GeometryCalib.Cx_Cy;
     221end
     222f1=GeometryCalib.focal*sx/dpx_dpy(1);
     223f2=GeometryCalib.focal/dpx_dpy(2);
     224Cx=Cx_Cy(1);
     225Cy=Cx_Cy(2);
     226set(handles.fx,'String',num2str(f1,4))
     227set(handles.fy,'String',num2str(f2,4))
     228set(handles.k,'String',num2str(k,'%1.4f'))
     229set(handles.Cx,'String',num2str(Cx,'%1.1f'))
     230set(handles.Cy,'String',num2str(Cy,'%1.1f'))
     231% Display extrinsinc parameters
     232set(handles.Tx,'String',num2str(GeometryCalib.Tx_Ty_Tz(1),4))
     233set(handles.Ty,'String',num2str(GeometryCalib.Tx_Ty_Tz(2),4))
     234set(handles.Tz,'String',num2str(GeometryCalib.Tx_Ty_Tz(3),4))
     235
     236% indicate the plane of the calibration grid if defined
    314237huvmat=findobj(allchild(0),'Name','uvmat');
    315238hhuvmat=guidata(huvmat);%handles of elements in the GUI uvmat
     
    337260    hhh=findobj(hhuvmat.axes3,'Tag','calib_marker');% delete calib points and markers
    338261    if ~isempty(hhh)
    339         delete(hhh);
     262        delete(hhh);     
    340263    end
    341264    hhh=findobj(hhuvmat.axes3,'Tag','calib_points');
     
    346269    set(hhuvmat.FixedLimits,'BackgroundColor',[0.7 0.7 0.7])
    347270    uvmat('RootPath_Callback',hObject,eventdata,hhuvmat); %file input with xml reading  in uvmat
    348     MenuPlot_Callback(hObject, eventdata, handles)
    349     set(handles.ListCoord,'Value',index)% indicate in the list the point with max deviation (possible mistake)
    350     ListCoord_Callback(hObject, eventdata, handles)
     271%     if ip==0
     272        MenuPlot_Callback(hObject, eventdata, handles)
     273        set(handles.ListCoord,'Value',index)% indicate in the list the point with max deviation (possible mistake)
     274        ListCoord_Callback(hObject, eventdata, handles)
     275%     end
    351276    figure(handles.geometry_calib)
    352277end
     
    355280% --- Executes on button press in calibrate_lin.
    356281function REPLICATE_Callback(hObject, eventdata, handles)
     282% %%%%%%Todo: correct on the model of APPLY_Callback%%%%%%%%%%
    357283%------------------------------------------------------------------------
    358284calib_cell=get(handles.calib_type,'String');
    359285val=get(handles.calib_type,'Value');
    360 calib_type=calib_cell{val};
    361286Coord_cell=get(handles.ListCoord,'String');
    362287Object=read_geometry_calib(Coord_cell);
    363 
    364 if isequal(calib_type,'rescale')
    365     GeometryCalib=calib_rescale(Object.Coord);
    366 elseif isequal(calib_type,'linear')
    367     GeometryCalib=calib_linear(Object.Coord);
    368 elseif isequal(calib_type,'tsai')
    369     GeometryCalib=calib_tsai(Object.Coord);
    370 end
     288GeometryCalib=feval(calib_cell{val},Object.Coord);
     289
    371290% %record image source
    372291GeometryCalib.SourceCalib.PointCoord=Object.Coord;
     
    424343%------------------------------------------------------------------------
    425344% determine the parameters for a calibration by an affine function (rescaling and offset, no rotation)
    426 function GeometryCalib=calib_rescale(Coord)
    427 %------------------------------------------------------------------------
    428  
     345function GeometryCalib=calib_rescale(Coord,handles)
     346%------------------------------------------------------------------------
    429347X=Coord(:,1);
    430 Y=Coord(:,2);
     348Y=Coord(:,2);% Z not used
    431349x_ima=Coord(:,4);
    432350y_ima=Coord(:,5);
     
    436354T_y=py(2);
    437355GeometryCalib.CalibrationType='rescale';
    438 GeometryCalib.focal=1;
     356GeometryCalib.focal=1;%default
     357GeometryCalib.sx=px(1)/py(1);
    439358GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function
    440 GeometryCalib.Tx_Ty_Tz=[T_x T_y 1];
    441 GeometryCalib.R=[px(1),0,0;0,py(1),0;0,0,1];
    442 
    443 % %check error
    444 % Calib.dpx=1;
    445 % Calib.dpy=1;
    446 % Calib.sx=1;
    447 % Calib.Cx=0;
    448 % Calib.Cy=0;
    449 % Calib.Tz=1;
    450 % Calib.kappa1=0;
    451 % Calib.f=GeometryCalib.focal;
    452 % Calib.Tx=T_x;
    453 % Calib.Ty=T_y;
    454 % Calib.R=GeometryCalib.R;
    455 % [Xpoints,Ypoints]=px_XYZ(Calib,X,Y,0);
    456 % GeometryCalib.ErrorRms(1)=sqrt(mean((Xpoints-x_ima).*(Xpoints-x_ima)));
    457 % GeometryCalib.ErrorMax(1)=max(abs(Xpoints-x_ima));
    458 % GeometryCalib.ErrorRms(2)=sqrt(mean((Ypoints-y_ima).*(Ypoints-y_ima)));
    459 % GeometryCalib.ErrorMax(2)=max(abs(Ypoints-y_ima));
     359GeometryCalib.Tx_Ty_Tz=[T_x T_y GeometryCalib.focal/py(1)];
     360GeometryCalib.R=[1,0,0;0,1,0;0,0,0];
    460361
    461362%------------------------------------------------------------------------
    462363% determine the parameters for a calibration by a linear transform matrix (rescale and rotation)
    463 function GeometryCalib=calib_linear(Coord)
     364function GeometryCalib=calib_linear(Coord,handles)
    464365%------------------------------------------------------------------------
    465366X=Coord(:,1);
    466 Y=Coord(:,2);
     367Y=Coord(:,2);% Z not used
    467368x_ima=Coord(:,4);
    468369y_ima=Coord(:,5);
    469370XY_mat=[ones(size(X)) X Y];
    470371a_X1=XY_mat\x_ima; %transformation matrix for X
    471 x1=XY_mat*a_X1;%reconstruction
    472 err_X1=max(abs(x1-x_ima));%error
     372% x1=XY_mat*a_X1;%reconstruction
     373% err_X1=max(abs(x1-x_ima));%error
    473374a_Y1=XY_mat\y_ima;%transformation matrix for X
    474 y1=XY_mat*a_Y1;
    475 err_Y1=max(abs(y1-y_ima));%error
    476 T_x=a_X1(1);
    477 T_y=a_Y1(1);
     375% y1=XY_mat*a_Y1;
     376% err_Y1=max(abs(y1-y_ima));%error
    478377GeometryCalib.CalibrationType='linear';
    479378GeometryCalib.focal=1;
    480379GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function
    481 GeometryCalib.Tx_Ty_Tz=[T_x T_y 1];
    482 GeometryCalib.R=[a_X1(2),a_X1(3),0;a_Y1(2),a_Y1(3),0;0,0,1];
    483 
    484 % %check error
    485 % GeometryCalib.ErrorRms(1)=sqrt(mean((x1-x_ima).*(x1-x_ima)));
    486 % GeometryCalib.ErrorMax(1)=max(abs(x1-x_ima));
    487 % GeometryCalib.ErrorRms(2)=sqrt(mean((y1-y_ima).*(y1-y_ima)));
    488 % GeometryCalib.ErrorMax(2)=max(abs(y1-y_ima));
    489 
    490 %------------------------------------------------------------------------
    491 function GeometryCalib=calib_tsai2(Coord)
     380R=[a_X1(2),a_X1(3),0;a_Y1(2),a_Y1(3),0;0,0,0];
     381norm=det(R);
     382GeometryCalib.Tx_Ty_Tz=[a_X1(1) a_Y1(1) norm];
     383GeometryCalib.R=R/norm;
     384
     385%------------------------------------------------------------------------
     386% determine the tsai parameters for a view normal to the grid plane
     387function GeometryCalib=calib_normal(Coord,handles)
     388%------------------------------------------------------------------------
     389Calib.f1=str2num(get(handles.fx,'String'));
     390Calib.f2=str2num(get(handles.fy,'String'));
     391Calib.k=str2num(get(handles.k,'String'));
     392Calib.Cx=str2num(get(handles.Cx,'String'));
     393Calib.Cy=str2num(get(handles.Cy,'String'));
     394%default
     395if isempty(Calib.f1)
     396    Calib.f1=25/0.012;
     397end
     398if isempty(Calib.f2)
     399    Calib.f2=25/0.012;
     400end
     401if isempty(Calib.k)
     402    Calib.k=0;
     403end
     404if isempty(Calib.Cx)||isempty(Calib.Cy)
     405    huvmat=findobj(allchild(0),'Tag','uvmat');
     406    hhuvmat=guidata(huvmat);
     407    Calib.Cx=str2num(get(hhuvmat.npx,'String'))/2;
     408    Calib.Cx=str2num(get(hhuvmat.npy,'String'))/2;
     409end   
     410%tsai parameters
     411Calib.dpx=0.012;%arbitrary
     412Calib.dpy=0.012;
     413Calib.sx=Calib.f1*Calib.dpx/(Calib.f2*Calib.dpy);
     414Calib.f=Calib.f2*Calib.dpy;
     415Calib.kappa1=Calib.k/(Calib.f*Calib.f);
     416
     417%initial guess
     418X=Coord(:,1);
     419Y=Coord(:,2);
     420Zmean=mean(Coord(:,3));
     421x_ima=Coord(:,4)-Calib.Cx;
     422y_ima=Coord(:,5)-Calib.Cy;
     423XY_mat=[ones(size(X)) X Y];
     424a_X1=XY_mat\x_ima; %transformation matrix for X
     425a_Y1=XY_mat\y_ima;%transformation matrix for Y
     426R=[a_X1(2),a_X1(3),0;a_Y1(2),a_Y1(3),0;0,0,-1];% rotation+ z axis reversal (upward)
     427norm=sqrt(det(-R));
     428calib_param(1)=0;% quadratic distortion
     429calib_param(2)=a_X1(1);
     430calib_param(3)=a_Y1(1);
     431calib_param(4)=Calib.f/(norm*Calib.dpx)-R(3,3)*Zmean;
     432calib_param(5)=angle(a_X1(2)+i*a_X1(3));
     433display(['initial guess=' num2str(calib_param)])
     434
     435%optimise the parameters: minimisation of error
     436calib_param = fminsearch(@(calib_param) error_calib(calib_param,Calib,Coord),calib_param)
     437
     438GeometryCalib.CalibrationType='tsai_normal';
     439GeometryCalib.focal=Calib.f;
     440GeometryCalib.dpx_dpy=[Calib.dpx Calib.dpy];
     441GeometryCalib.Cx_Cy=[Calib.Cx Calib.Cy];
     442GeometryCalib.sx=Calib.sx;
     443GeometryCalib.kappa1=calib_param(1);
     444GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function
     445GeometryCalib.Tx_Ty_Tz=[calib_param(2) calib_param(3) calib_param(4)];
     446alpha=calib_param(5);
     447GeometryCalib.R=[cos(alpha) sin(alpha) 0;-sin(alpha) cos(alpha) 0;0 0 -1];
     448
     449%------------------------------------------------------------------------
     450function GeometryCalib=calib_3D_linear(Coord,handles)
    492451%------------------------------------------------------------------
    493452path_uvmat=which('uvmat');% check the path detected for source file uvmat
     
    495454huvmat=findobj(allchild(0),'Tag','uvmat');
    496455hhuvmat=guidata(huvmat);
     456x_1=Coord(:,4:5)'
     457X_1=Coord(:,1:3)'
     458n_ima=1;
     459% check_cond=0;
     460
     461nx=str2num(get(hhuvmat.npx,'String'));
     462ny=str2num(get(hhuvmat.npy,'String'));
     463
     464est_dist=[0;0;0;0;0];
     465est_aspect_ratio=0;
     466est_fc=[1;1];
     467%fc=[25;25]/0.012;
     468center_optim=0;
     469run(fullfile(path_UVMAT,'toolbox_calib','go_calib_optim'));
     470
     471GeometryCalib.CalibrationType='3D_linear';
     472GeometryCalib.focal=fc(2);
     473GeometryCalib.dpx_dpy=[1 1];
     474GeometryCalib.Cx_Cy=cc';
     475GeometryCalib.sx=fc(1)/fc(2);
     476GeometryCalib.kappa1=-kc(1)/fc(2)^2;
     477GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function
     478GeometryCalib.Tx_Ty_Tz=Tc_1';
     479GeometryCalib.R=Rc_1;
     480
     481%------------------------------------------------------------------------
     482function GeometryCalib=calib_3D_quadr(Coord,handles)
     483%------------------------------------------------------------------
     484
     485path_uvmat=which('uvmat');% check the path detected for source file uvmat
     486path_UVMAT=fileparts(path_uvmat); %path to UVMAT
     487huvmat=findobj(allchild(0),'Tag','uvmat');
     488hhuvmat=guidata(huvmat);
     489% check_cond=0;
     490coord_files=get(handles.coord_files,'String');
     491if ischar(coord_files)
     492    coord_files={coord_files};
     493end
     494if isempty(coord_files{1}) || isequal(coord_files,{''})
     495    coord_files={};
     496end
     497
     498%retrieve the calibration points stored in the files listed in the popup list coord_files
     499x_1=Coord(:,4:5)';
     500X_1=Coord(:,1:3)';
     501n_ima=numel(coord_files)+1;
     502if ~isempty(coord_files)
     503    msgbox_uvmat('CONFIRMATION',['The xy coordinates of the calibration points in ' num2str(n_ima) ' planes will be used'])
     504    for ifile=1:numel(coord_files)
     505    t=xmltree(coord_files{ifile});
     506    s=convert(t);%convert to matlab structure
     507        if isfield(s,'GeometryCalib')
     508            if isfield(s.GeometryCalib,'SourceCalib')
     509                if isfield(s.GeometryCalib.SourceCalib,'PointCoord')
     510                PointCoord=s.GeometryCalib.SourceCalib.PointCoord;
     511                Coord_file=[];
     512                for i=1:length(PointCoord)
     513                    line=str2num(PointCoord{i});
     514                    Coord_file(i,4:5)=line(4:5);%px x
     515                    Coord_file(i,1:3)=line(1:3);%phys x
     516                end
     517                eval(['x_' num2str(ifile+1) '=Coord_file(:,4:5)'';']);
     518                eval(['X_' num2str(ifile+1) '=Coord_file(:,1:3)'';']);
     519                end
     520            end
     521        end
     522    end
     523end
     524
     525n_ima=numel(coord_files)+1;
     526whos
     527nx=str2num(get(hhuvmat.npx,'String'));
     528ny=str2num(get(hhuvmat.npy,'String'));
     529
     530est_dist=[1;0;0;0;0];
     531est_aspect_ratio=1;
     532%est_fc=[0;0];
     533%fc=[25;25]/0.012;
     534center_optim=0;
     535run(fullfile(path_UVMAT,'toolbox_calib','go_calib_optim'));
     536
     537GeometryCalib.CalibrationType='3D_quadr';
     538GeometryCalib.focal=fc(2);
     539GeometryCalib.dpx_dpy=[1 1];
     540GeometryCalib.Cx_Cy=cc';
     541GeometryCalib.sx=fc(1)/fc(2);
     542GeometryCalib.kappa1=-kc(1)/fc(2)^2;
     543GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function
     544GeometryCalib.Tx_Ty_Tz=Tc_1';
     545GeometryCalib.R=Rc_1;
     546GeometryCalib.R(3,1:3)=-GeometryCalib.R(3,1:3);%inversion for z upward
     547GeometryCalib.ErrorRMS=[];
     548GeometryCalib.ErrorMax=[];
     549
     550
     551%------------------------------------------------------------------------
     552function GeometryCalib=calib_3D_extrinsic(Coord,handles)
     553%------------------------------------------------------------------
     554path_uvmat=which('geometry_calib');% check the path detected for source file uvmat
     555path_UVMAT=fileparts(path_uvmat); %path to UVMAT
    497556x_1=Coord(:,4:5)';
    498557X_1=Coord(:,1:3)';
    499558n_ima=1;
    500 % check_cond=0;
    501 
    502 
    503 nx=str2num(get(hhuvmat.npx,'String'));
    504 ny=str2num(get(hhuvmat.npy,'String'));
    505 
    506 
    507 % est_kc=[1;0;0;0;0];
    508 est_fc=0;
    509 est_dist=[1;0;0;0;0];
    510 run(fullfile(path_UVMAT,'toolbox_calib','go_calib_optim'));
    511 
    512 GeometryCalib.CalibrationType='tsai_matlab';
    513 GeometryCalib.focal=f(2);
     559Calib.f1=str2num(get(handles.fx,'String'));
     560Calib.f2=str2num(get(handles.fy,'String'));
     561Calib.k=str2num(get(handles.k,'String'));
     562Calib.Cx=str2num(get(handles.Cx,'String'));
     563Calib.Cy=str2num(get(handles.Cy,'String'));
     564
     565fct_path=fullfile(path_UVMAT,'toolbox_calib');
     566addpath(fct_path)
     567% [omc1,Tc1,Rc1,H,x,ex,JJ] = compute_extrinsic(x_1,X_1,...
     568%     [Calib.f Calib.f*Calib.sx]',...
     569%     [Calib.Cx Calib.Cy]',...
     570%     [-Calib.kappa1*Calib.f^2 0 0 0 0]);
     571[omc1,Tc1,Rc1,H,x,ex,JJ] = compute_extrinsic(x_1,X_1,...
     572    [Calib.f1 Calib.f2]',...
     573    [Calib.Cx Calib.Cy]',...
     574    [-Calib.k 0 0 0 0]);
     575%get the euler angles ???
     576rmpath(fct_path);
     577
     578std(ex')
     579GeometryCalib.CalibrationType='3D_extrinsic';
     580GeometryCalib.focal=Calib.f2;
    514581GeometryCalib.dpx_dpy=[1 1];
    515 GeometryCalib.Cx_Cy=cc';
    516 GeometryCalib.sx=fc(1)/fc(2);
    517 GeometryCalib.kappa1=-k(1)/f(2)^2;
     582GeometryCalib.Cx_Cy=[Calib.Cx Calib.Cy];
     583GeometryCalib.sx=Calib.f1/Calib.f2;
     584GeometryCalib.kappa1=Calib.k/(Calib.f2*Calib.f2);
    518585GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function
    519 GeometryCalib.Tx_Ty_Tz=Tc_1';
    520 GeometryCalib.R=Rc_1;
    521 % Calib.dpx=GeometryCalib.dpx_dpy(1);
    522 % Calib.dpy=GeometryCalib.dpx_dpy(2);
    523 % Calib.sx=GeometryCalib.sx;
    524 % Calib.Cx=GeometryCalib.Cx_Cy(1);
    525 % Calib.Cy=GeometryCalib.Cx_Cy(2);
    526 % Calib.kappa1=GeometryCalib.kappa1;
    527 % Calib.f=GeometryCalib.focal;
    528 % Calib.Tx=GeometryCalib.Tx_Ty_Tz(1);
    529 % Calib.Ty=GeometryCalib.Tx_Ty_Tz(2);
    530 % Calib.Tz=GeometryCalib.Tx_Ty_Tz(3);
    531 % Calib.R=GeometryCalib.R;
    532 % X=Coord(:,1);
    533 % Y=Coord(:,2);
    534 % Z=Coord(:,3);
    535 % x_ima=Coord(:,4);
    536 % y_ima=Coord(:,5);
    537 % [Xpoints,Ypoints]=px_XYZ(Calib,X,Y,Z);
    538 %
    539 % GeometryCalib.ErrorRms(1)=sqrt(mean((Xpoints-x_ima).*(Xpoints-x_ima)));
    540 % [GeometryCalib.ErrorMax(1),GeometryCalib.IndexMax(1)]=max(abs(Xpoints-x_ima));
    541 % GeometryCalib.ErrorRms(2)=sqrt(mean((Ypoints-y_ima).*(Ypoints-y_ima)));
    542 % [GeometryCalib.ErrorMax(2),GeometryCalib.IndexMax(2)]=max(abs(Ypoints-y_ima));
    543 
    544 function GeometryCalib=calib_tsai(Coord)
     586GeometryCalib.Tx_Ty_Tz=Tc1';
     587%inversion of z axis
     588GeometryCalib.R=Rc1;
     589%GeometryCalib.R(3,1:3)=-GeometryCalib.R(3,1:3);%inversion for z upward
     590
     591
     592%------------------------------------------------------------------------
     593%function GeometryCalib=calib_tsai_heikkila(Coord)
     594% TEST: NOT IMPLEMENTED
     595%------------------------------------------------------------------
     596% path_uvmat=which('uvmat');% check the path detected for source file uvmat
     597% path_UVMAT=fileparts(path_uvmat); %path to UVMAT
     598% path_calib=fullfile(path_UVMAT,'toolbox_calib_heikkila');
     599% addpath(path_calib)
     600% npoints=size(Coord,1);
     601% Coord(:,1:3)=10*Coord(:,1:3);
     602% Coord=[Coord zeros(npoints,2) -ones(npoints,1)];
     603% [par,pos,iter,res,er,C]=cacal('dalsa',Coord);
     604% GeometryCalib.CalibrationType='tsai';
     605% GeometryCalib.focal=par(2);
     606
     607
     608%--------------------------------------------------------------------------
     609function GeometryCalib=calib_tsai(Coord,handles)% old version using gauthier's bianry ccal_fo
    545610%------------------------------------------------------------------------
    546611%TSAI
    547 % 'calibration_lin' provides a linear transform on coordinates,
    548612path_uvmat=which('uvmat');% check the path detected for source file uvmat
    549613path_UVMAT=fileparts(path_uvmat); %path to UVMAT
    550 % if isunix
    551     %fid = fopen(fullfile(path_UVMAT,'PARAM_LINUX.txt'),'r');%open the file with civ binary names
    552 xmlfile=fullfile(path_UVMAT,'PARAM.xml');
     614xmlfile=fullfile(path_UVMAT,'PARAM.xml');%name of the file containing names of binary executables
    553615if exist(xmlfile,'file')
    554     t=xmltree(xmlfile);
    555     sparam=convert(t);
     616    t=xmltree(xmlfile);% read the (xml) file containing names of binary executables
     617    sparam=convert(t);% convert to matlab structure
    556618end
    557619if ~isfield(sparam,'GeometryCalibBin')
     
    570632textcoord=num2str(Coord,4);
    571633dlmwrite('t.txt',textcoord,''); 
    572 % ['!' Tsai_exe ' -f1 0 -f2 t.txt']
    573     eval(['!' Tsai_exe ' -f t.txt > tsaicalib.log']);
     634% ['!' Tsai_exe ' -fx 0 -fy t.txt']
     635eval(['!' Tsai_exe ' -f t.txt > tsaicalib.log']);
    574636if ~exist('calib.dat','file')
    575637    msgbox_uvmat('ERROR','no output from calibration program Tsai_exe: possibly too few points')
     
    577639calibdat=dlmread('calib.dat');
    578640delete('calib.dat')
    579 delete('t.txt')
     641%delete('t.txt')
    580642GeometryCalib.CalibrationType='tsai';
    581643GeometryCalib.focal=calibdat(10);
     
    604666%EN DEDUIRE MATRICE R ??
    605667GeometryCalib.R=[r1,r2,r3;r4,r5,r6;r7,r8,r9];
    606 %erreur a caracteriser?
    607 % %check error
    608 % Calib.dpx=GeometryCalib.dpx_dpy(1);
    609 % Calib.dpy=GeometryCalib.dpx_dpy(2);
    610 % Calib.sx=GeometryCalib.sx;
    611 % Calib.Cx=GeometryCalib.Cx_Cy(1);
    612 % Calib.Cy=GeometryCalib.Cx_Cy(2);
    613 % Calib.kappa1=GeometryCalib.kappa1;
    614 % Calib.f=GeometryCalib.focal;
    615 % Calib.Tx=GeometryCalib.Tx_Ty_Tz(1);
    616 % Calib.Ty=GeometryCalib.Tx_Ty_Tz(2);
    617 % Calib.Tz=GeometryCalib.Tx_Ty_Tz(3);
    618 % Calib.R=GeometryCalib.R;
    619 % X=Coord(:,1);
    620 % Y=Coord(:,2);
    621 % Z=Coord(:,3);
    622 % x_ima=Coord(:,4);
    623 % y_ima=Coord(:,5);
    624 % [Xpoints,Ypoints]=px_XYZ(Calib,X,Y,Z);
    625 %
    626 % GeometryCalib.ErrorRms(1)=sqrt(mean((Xpoints-x_ima).*(Xpoints-x_ima)));
    627 % GeometryCalib.ErrorMax(1)=max(abs(Xpoints-x_ima));
    628 % GeometryCalib.ErrorRms(2)=sqrt(mean((Ypoints-y_ima).*(Ypoints-y_ima)));
    629 % GeometryCalib.ErrorMax(2)=max(abs(Ypoints-y_ima));
    630 
    631 
    632 %------------------------------------------------------------------------
    633 % --- Executes on button press in rotation.
    634 function rotation_Callback(hObject, eventdata, handles)
    635 %------------------------------------------------------------------------
    636 angle_rot=(pi/180)*str2num(get(handles.Phi,'String'));
     668
     669%------------------------------------------------------------------------
     670% --- determine the rms of calibration error
     671function ErrorRms=error_calib(calib_param,Calib,Coord)
     672%calib_param: vector of free calibration parameters (to optimise)
     673%Calib: structure of the given calibration parameters
     674%Coord: list of phys coordinates (columns 1-3, and pixel coordinates (columns 4-5)
     675Calib.f=25;
     676Calib.dpx=0.012;
     677Calib.dpy=0.012;
     678Calib.sx=1;
     679Calib.Cx=512;
     680Calib.Cy=512;
     681Calib.kappa1=calib_param(1);
     682Calib.Tx=calib_param(2);
     683Calib.Ty=calib_param(3);
     684Calib.Tz=calib_param(4);
     685alpha=calib_param(5);
     686Calib.R=[cos(alpha) sin(alpha) 0;-sin(alpha) cos(alpha) 0;0 0 -1];
     687
     688X=Coord(:,1);
     689Y=Coord(:,2);
     690Z=Coord(:,3);
     691x_ima=Coord(:,4);
     692y_ima=Coord(:,5);
     693[Xpoints,Ypoints]=px_XYZ(Calib,X,Y,Z);
     694ErrorRms(1)=sqrt(mean((Xpoints-x_ima).*(Xpoints-x_ima)));
     695ErrorRms(2)=sqrt(mean((Ypoints-y_ima).*(Ypoints-y_ima)));
     696ErrorRms=mean(ErrorRms);
     697
     698%------------------------------------------------------------------------
     699function XImage_Callback(hObject, eventdata, handles)
     700%------------------------------------------------------------------------
     701update_list(hObject, eventdata,handles)
     702
     703%------------------------------------------------------------------------
     704function YImage_Callback(hObject, eventdata, handles)
     705%------------------------------------------------------------------------
     706update_list(hObject, eventdata,handles)
     707
     708%------------------------------------------------------------------------
     709% --- Executes on button press in STORE.
     710function STORE_Callback(hObject, eventdata, handles)
    637711Coord_cell=get(handles.ListCoord,'String');
    638 data=read_geometry_calib(Coord_cell);
    639 data.Coord(:,1)=cos(angle_rot)*data.Coord(:,1)+sin(angle_rot)*data.Coord(:,2);
    640 data.Coord(:,1)=-sin(angle_rot)*data.Coord(:,1)+cos(angle_rot)*data.Coord(:,2);
    641 set(handles.XObject,'String',num2str(data.Coord(:,1),4));
    642 set(handles.YObject,'String',num2str(data.Coord(:,2),4));
    643 
    644 %------------------------------------------------------------------------
    645 function XImage_Callback(hObject, eventdata, handles)
     712Object=read_geometry_calib(Coord_cell);
     713unitlist=get(handles.CoordUnit,'String');
     714unit=unitlist{get(handles.CoordUnit,'value')};
     715GeometryCalib.CoordUnit=unit;
     716GeometryCalib.SourceCalib.PointCoord=Object.Coord;
     717huvmat=findobj(allchild(0),'Name','uvmat');
     718hhuvmat=guidata(huvmat);%handles of elements in the GUI uvmat
     719RootPath='';
     720RootFile='';
     721if ~isempty(hhuvmat.RootPath)& ~isempty(hhuvmat.RootFile)
     722    testhandle=1;
     723    RootPath=get(hhuvmat.RootPath,'String');
     724    RootFile=get(hhuvmat.RootFile,'String');
     725    filebase=fullfile(RootPath,RootFile);
     726    while exist([filebase '.xml'],'file')
     727        filebase=[filebase '~'];
     728    end
     729    outputfile=[filebase '.xml'];
     730    update_imadoc(GeometryCalib,outputfile)
     731    listfile=get(handles.coord_files,'string');
     732    if isequal(listfile,{''})
     733        listfile={outputfile};
     734    else
     735        listfile=[listfile;{outputfile}]%update the list of coord files
     736    end
     737    set(handles.coord_files,'string',listfile);
     738end
     739set(handles.ListCoord,'Value',1)% refresh the display of coordinates
     740set(handles.ListCoord,'String',{'......'})
     741
     742%------------------------------------------------------------------------
     743% --- Executes on button press in CLEAR.
     744function CLEAR_Callback(hObject, eventdata, handles)
     745%------------------------------------------------------------------------
     746set(handles.coord_files,'Value',1)
     747set(handles.coord_files,'String',{''})
     748
     749%------------------------------------------------------------------------
     750function XObject_Callback(hObject, eventdata, handles)
    646751%------------------------------------------------------------------------
    647752update_list(hObject, eventdata,handles)
    648753
    649754%------------------------------------------------------------------------
    650 function YImage_Callback(hObject, eventdata, handles)
     755function YObject_Callback(hObject, eventdata, handles)
    651756%------------------------------------------------------------------------
    652757update_list(hObject, eventdata,handles)
    653758
    654 function XObject_Callback(hObject, eventdata, handles)
    655 update_list(hObject, eventdata,handles)
    656 
    657 function YObject_Callback(hObject, eventdata, handles)
    658 update_list(hObject, eventdata,handles)
    659 
     759%------------------------------------------------------------------------
    660760function ZObject_Callback(hObject, eventdata, handles)
     761%------------------------------------------------------------------------
    661762update_list(hObject, eventdata,handles)
    662763
     
    748849%------------------------------------------------------------------------
    749850choice=get(handles.edit_append,'Value');
    750 % if choice==1
    751 %        Coord=get(handles.ListCoord,'String');
    752 %        val=length(Coord);
    753 %        if val>=1 & isequal(Coord{val},'')
    754 %             val=val-1; %do not take into account blank
    755 %        end
    756 %        Coord{val+1}='';
    757 %        set(handles.ListCoord,'String',Coord)
    758 %        set(handles.ListCoord,'Value',val+1)
    759 % end
    760851if choice
    761852    set(handles.edit_append,'BackgroundColor',[1 1 0])
     
    817908end
    818909
    819 % %------------------------------------------------------------------------
    820 % % --- Executes on button press in append_point.
    821 % function append_point_Callback(hObject, eventdata, handles)
    822 % %------------------------------------------------------------------------
    823 %        Coord=get(handles.ListCoord,'String');
    824 %        val=length(Coord);
    825 %        if val>=1 & isequal(Coord{val},'')
    826 %             val=val-1; %do not take into account blank
    827 %        end
    828 %        Coord{val+1}='';
    829 %        set(handles.ListCoord,'String',Coord)
    830 %        set(handles.ListCoord,'Value',val+1)
    831 
    832 %------------------------------------------------------------------------
    833 function MenuOpen_Callback(hObject, eventdata, handles)
    834 %------------------------------------------------------------------------
    835 %get the object file
    836 huvmat=findobj(allchild(0),'Name','uvmat');
    837 UvData=get(huvmat,'UserData');
    838 hchild=get(huvmat,'Children');
    839 hrootpath=findobj(hchild,'Tag','RootPath');
    840 oldfile=get(hrootpath,'String');
    841 if isempty(oldfile)
    842     oldfile='';
    843 end
    844 %[FileName,PathName] = uigetfile('*.civ','Select a .civ file',oldfile)
    845 [FileName, PathName, filterindex] = uigetfile( ...
    846        {'*.xml;*.mat', ' (*.xml,*.mat)';
    847        '*.xml',  '.xml files '; ...
    848         '*.mat',  '.mat matlab files '}, ...
    849         'Pick a file',oldfile);
    850 fileinput=[PathName FileName];%complete file name
    851 testblank=findstr(fileinput,' ');%look for blanks
    852 if ~isempty(testblank)
    853     msgbox_uvmat('ERROR','forbidden input file name or path: no blank character allowed')
    854     return
    855 end
    856 sizf=size(fileinput);
    857 if (~ischar(fileinput)|~isequal(sizf(1),1)),return;end
    858 loadfile(handles,fileinput)
    859 
    860 
    861910%------------------------------------------------------------------------
    862911function MenuPlot_Callback(hObject, eventdata, handles)
     
    9551004    grid_input=CalibData.grid;%retrieve the previously used grid
    9561005end
    957 [T,CalibData.grid]=create_grid(grid_input,'detect grid');%display the GUI create_grid, read the set of phys coordinates T
     1006[T,CalibData.grid]=create_grid(grid_input,'detect_grid');%display the GUI create_grid, read the set of phys coordinates T
     1007
    9581008set(handles.geometry_calib,'UserData',CalibData)%store the phys grid parameters for later use
    9591009
     
    9881038X=[CalibData.grid.x_0 CalibData.grid.x_1 CalibData.grid.x_0 CalibData.grid.x_1]';%corner absissa in the phys coordinates
    9891039Y=[CalibData.grid.y_0 CalibData.grid.y_0 CalibData.grid.y_1 CalibData.grid.y_1]';%corner ordinates in the phys coordinates
    990 %X=[CalibData.grid.x_0 CalibData.grid.x_0 CalibData.grid.x_1 CalibData.grid.x_1]';%corner absissa in the phys coordinates
    991 %Y=[CalibData.grid.y_0 CalibData.grid.y_1 CalibData.grid.y_1 CalibData.grid.y_0]';%corner ordinates in the phys coordinate
    992 %XY_mat=[ones(size(X)) X Y];
    993 %a_X1=XY_mat\corners_X; %transformation matrix for X
    994 %x1=XY_mat*a_X1;%reconstruction
    995 %err_X1=max(abs(x1-corners_X));%error
    996 %a_Y1=XY_mat\corners_Y;%transformation matrix for X
    997 %y1=XY_mat*a_Y1;
    998 %err_Y1=max(abs(y1-corners_Y));%error
    999 % figure(1)
    1000 % hold off
    1001 % plot([corners_X;corners_X(1)],[corners_Y;corners_Y(1)],'r')
    1002 %test points along contour
    1003 % u=[0:0.1:1];
    1004 % x12=u*corners_X(2)+(1-u)*corners_X(1);
    1005 % x23=u*corners_X(3)+(1-u)*corners_X(2);
    1006 % x34=u*corners_X(4)+(1-u)*corners_X(3);
    1007 % x41=u*corners_X(1)+(1-u)*corners_X(4);
    1008 % y12=u*corners_Y(2)+(1-u)*corners_Y(1);
    1009 % y23=u*corners_Y(3)+(1-u)*corners_Y(2);
    1010 % y34=u*corners_Y(4)+(1-u)*corners_Y(3);
    1011 % y41=u*corners_Y(1)+(1-u)*corners_Y(4);
    1012 % hold on
    1013 % plot(x12,y12,'xr')
    1014 % plot(x23,y23,'xr')
    1015 % plot(x34,y34,'xr')
    1016 % plot(x41,y41,'xr')
     1040
    10171041%calculate transform matrices:
    10181042% reference: http://alumni.media.mit.edu/~cwren/interpolator/ by Christopher R. Wren
     
    10261050C = [l(7:8)' 1];
    10271051
    1028 GeometryCalib.CalibrationType='tsai';%'linear';
     1052GeometryCalib.CalibrationType='tsai';
    10291053GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function
    10301054GeometryCalib.f=1;
     
    10391063GeometryCalib.Tz=1;
    10401064GeometryCalib.R=[Amat(1,1),Amat(1,2),0;Amat(2,1),Amat(2,2),0;C(1),C(2),0];
    1041 %montre points du pourtour transformes
    1042 % for ip=1:numel(u),
    1043 %   XY12(:,ip)=Amat*[x12(ip);y12(ip);1]/(C*[x12(ip);y12(ip);1]);
    1044 %   XY23(:,ip)=Amat*[x23(ip);y23(ip);1]/(C*[x23(ip);y23(ip);1]);
    1045 %   XY34(:,ip)=Amat*[x34(ip);y34(ip);1]/(C*[x34(ip);y34(ip);1]);
    1046 %   XY41(:,ip)=Amat*[x41(ip);y41(ip);1]/(C*[x41(ip);y41(ip);1]);
    1047 % end
    1048 % xphys12=u*X(2)+(1-u)*X(1);
    1049 % xphys23=u*X(3)+(1-u)*X(2);
    1050 % xphys34=u*X(4)+(1-u)*X(3);
    1051 % xphys41=u*X(1)+(1-u)*X(4);
    1052 % yphys12=u*Y(2)+(1-u)*Y(1);
    1053 % yphys23=u*Y(3)+(1-u)*Y(2);
    1054 % yphys34=u*Y(4)+(1-u)*Y(3);
    1055 % yphys41=u*Y(1)+(1-u)*Y(4);
    1056 % for ip=1:numel(u),
    1057 %   XY12(:,ip)=Amat*[xphys12(ip);yphys12(ip);1]/(C*[xphys12(ip);yphys12(ip);1]);
    1058 %   XY23(:,ip)=Amat*[xphys23(ip);yphys23(ip);1]/(C*[xphys23(ip);yphys23(ip);1]);
    1059 %   XY34(:,ip)=Amat*[xphys34(ip);yphys34(ip);1]/(C*[xphys34(ip);yphys34(ip);1]);
    1060 %   XY41(:,ip)=Amat*[xphys41(ip);yphys41(ip);1]/(C*[xphys41(ip);yphys41(ip);1]);
    1061 % end
    1062 % x12=XY12(1,:);
    1063 % y12=XY12(2,:);
    1064 % x23=XY23(1,:);
    1065 % y23=XY23(2,:);
    1066 % x34=XY34(1,:);
    1067 % y34=XY34(2,:);
    1068 % x41=XY41(1,:);
    1069 % y41=XY41(2,:);
    1070 % [x12,y12]=px_XYZ(GeometryCalib,xphys12,yphys12);
    1071 % [x23,y23]=px_XYZ(GeometryCalib,xphys23,yphys23);
    1072 % [x34,y34]=px_XYZ(GeometryCalib,xphys34,yphys34);
    1073 % [x41,y41]=px_XYZ(GeometryCalib,xphys41,yphys41);
    1074 % plot(x12,y12,'xb')
    1075 % plot(x23,y23,'xb')
    1076 % plot(x34,y34,'xb')
    1077 % plot(x41,y41,'xb')
    1078 % figure(2)
    1079 % hold off
    1080 % plot(XY12(1,:),XY12(2,:),'ob') 
    1081 % hold on
    1082 % plot(XY23(1,:),XY23(2,:),'ob')
    1083 % plot(XY34(1,:),XY34(2,:),'ob')
    1084 % plot(XY41(1,:),XY41(2,:),'ob')
    1085 % plot(xphys12,yphys12,'ob') 
    1086 % hold on
    1087 % plot(xphys23,yphys23,'ob')
    1088 % plot(xphys34,yphys34,'ob')
    1089 % plot(xphys41,yphys41,'ob')
    10901065
    10911066[Amod,Rangx,Rangy]=phys_Ima(A-min(min(A)),GeometryCalib,0);
    10921067
    10931068Amod=double(Amod);
    1094 figure(12)
    1095 Amax=max(max(Amod));
    1096 image(Rangx,Rangy,uint8(255*Amod/Amax))
     1069% figure(12) display corrected image
     1070% Amax=max(max(Amod));
     1071% image(Rangx,Rangy,uint8(255*Amod/Amax))
    10971072Dx=(Rangx(2)-Rangx(1))/(npxy(2)-1); %x mesh in real space
    10981073Dy=(Rangy(2)-Rangy(1))/(npxy(1)-1); %y mesh in real space
     
    11031078    i0=1+round((T(ipoint,1)-Rangx(1))/Dx);%round(Xpx(ipoint));
    11041079    j0=1+round((T(ipoint,2)-Rangy(1))/Dy);%round(Xpx(ipoint));
    1105     Asub=Amod(j0-ind_range_y:j0+ind_range_y,i0-ind_range_x:i0+ind_range_x);
     1080    j0min=max(j0-ind_range_y,1);
     1081    j0max=min(j0+ind_range_y,size(Amod,1));
     1082    i0min=max(i0-ind_range_x,1);
     1083    i0max=min(i0+ind_range_x,size(Amod,2));
     1084    Asub=Amod(j0min:j0max,i0min:i0max);
    11061085    x_profile=sum(Asub,1);
    11071086    y_profile=sum(Asub,2);
     
    11111090    x_shift=0;
    11121091    y_shift=0;
    1113     if ind_x_max+2<=2*ind_range_x+1 && ind_x_max-2>=1
     1092    %if ind_x_max+2<=2*ind_range_x+1 && ind_x_max-2>=1
     1093    if ind_x_max+2<=numel(x_profile) && ind_x_max-2>=1
    11141094        Atop=x_profile(ind_x_max-2:ind_x_max+2);
    11151095        x_shift=sum(Atop.*[-2 -1 0 1 2])/sum(Atop);
    11161096    end
    1117     if ind_y_max+2<=2*ind_range_y+1 && ind_y_max-2>=1
     1097    if ind_x_max+2<=numel(y_profile) && ind_y_max-2>=1
    11181098        Atop=y_profile(ind_y_max-2:ind_y_max+2);
    11191099        y_shift=sum(Atop.*[-2 -1 0 1 2]')/sum(Atop);
    11201100    end
    1121     Delta(ipoint,1)=(x_shift+ind_x_max-ind_range_x-1)*Dx;%shift from the initial guess
    1122     Delta(ipoint,2)=(y_shift+ind_y_max-ind_range_y-1)*Dy;
     1101    Delta(ipoint,1)=(x_shift+ind_x_max+i0min-i0-1)*Dx;%shift from the initial guess
     1102    Delta(ipoint,2)=(y_shift+ind_y_max+j0min-j0-1)*Dy;
    11231103end
    11241104Tmod=T(:,(1:2))+Delta;
     
    11271107     Coord{ipoint,1}=num2str(T(ipoint,1),4);%display coordiantes with 4 digits
    11281108     Coord{ipoint,2}=num2str(T(ipoint,2),4);%display coordiantes with 4 digits
    1129      Coord{ipoint,3}='0';
    1130      Coord{ipoint,4}=num2str(Xpx(ipoint),4);%display coordiantes with 4 digi
    1131      Coord{ipoint,5}=num2str(Ypx(ipoint),4);%display coordiantes with 4 digi
     1109     Coord{ipoint,3}=num2str(T(ipoint,3),4);%display coordiantes with 4 digits;
     1110     Coord{ipoint,4}=num2str(Xpx(ipoint),4);%display coordiantes with 4 digits
     1111     Coord{ipoint,5}=num2str(Ypx(ipoint),4);%display coordiantes with 4 digits
    11321112end
    11331113Tabchar=cell2tab(Coord(end:-1:1,:),'    |    ');
     
    12091189set(handles.ListCoord,'Value',1)
    12101190set(handles.ListCoord,'String',Tabchar)
     1191
     1192
     1193% %------------------------------------------------------------------------
     1194% % --- Executes on button press in rotation.
     1195% function rotation_Callback(hObject, eventdata, handles)
     1196% %------------------------------------------------------------------------
     1197% angle_rot=(pi/180)*str2num(get(handles.Phi,'String'));
     1198% Coord_cell=get(handles.ListCoord,'String');
     1199% data=read_geometry_calib(Coord_cell);
     1200% data.Coord(:,1)=cos(angle_rot)*data.Coord(:,1)+sin(angle_rot)*data.Coord(:,2);
     1201% data.Coord(:,1)=-sin(angle_rot)*data.Coord(:,1)+cos(angle_rot)*data.Coord(:,2);
     1202% set(handles.XObject,'String',num2str(data.Coord(:,1),4));
     1203% set(handles.YObject,'String',num2str(data.Coord(:,2),4));
     1204
    12111205
    12121206%------------------------------------------------------------------------
     
    12871281    Zphys=Calib.SliceCoord(Zindex,3);%GENERALISER AUX CAS AVEC ANGLE
    12881282else
    1289 %     if exist('Z','var')
    1290 %         Zphys=Z;
    1291 %     else
    1292         Zphys=0;
    1293 %     end
     1283    Zphys=0;
    12941284end
    12951285if ~exist('X','var')||~exist('Y','var')
     
    13291319    Yphys=(A21.*Xu+A22.*Yu+Y0)./denom;
    13301320end
     1321
     1322
     1323% --------------------------------------------------------------------
     1324function MenuImportPoints_Callback(hObject, eventdata, handles)
     1325fileinput=browse_xml(hObject, eventdata, handles)
     1326if isempty(fileinput)
     1327    return
     1328end
     1329[s,errormsg]=imadoc2struct(fileinput,'GeometryCalib');
     1330GeometryCalib=s.GeometryCalib;
     1331GeometryCalib=load_calib(hObject, eventdata, handles)
     1332calib=reshape(GeometryCalib.PointCoord',[],1);
     1333for ilist=1:numel(calib)
     1334    CoordCell{ilist}=num2str(calib(ilist));
     1335end
     1336CoordCell=reshape(CoordCell,[],5);
     1337Tabchar=cell2tab(CoordCell,'    |    ');%transform cells into table ready for display
     1338Tabchar=[Tabchar;{'......'}];
     1339set(handles.ListCoord,'Value',1)
     1340set(handles.ListCoord,'String',Tabchar)
     1341MenuPlot_Callback(handles.geometry_calib, [], handles)
     1342
     1343% -----------------------------------------------------------------------
     1344function MenuImportIntrinsic_Callback(hObject, eventdata, handles)
     1345%------------------------------------------------------------------------
     1346fileinput=browse_xml(hObject, eventdata, handles);
     1347if isempty(fileinput)
     1348    return
     1349end
     1350[s,errormsg]=imadoc2struct(fileinput,'GeometryCalib');
     1351GeometryCalib=s.GeometryCalib;
     1352k=GeometryCalib.kappa1 * GeometryCalib.f*GeometryCalib.f;
     1353f1=GeometryCalib.f*GeometryCalib.sx/GeometryCalib.dpx;
     1354f2=GeometryCalib.f/GeometryCalib.dpy;
     1355Cx=GeometryCalib.Cx;
     1356Cy=GeometryCalib.Cy;
     1357set(handles.fx,'String',num2str(f1,'%1.1f'))
     1358set(handles.fy,'String',num2str(f2,'%1.1f'))
     1359set(handles.k,'String',num2str(k,'%1.4f'))
     1360set(handles.Cx,'String',num2str(Cx,'%1.1f'))
     1361set(handles.Cy,'String',num2str(Cy,'%1.1f'))
     1362
     1363% -----------------------------------------------------------------------
     1364function MenuImportAll_Callback(hObject, eventdata, handles)
     1365%------------------------------------------------------------------------
     1366fileinput=browse_xml(hObject, eventdata, handles)
     1367if ~isempty(fileinput)
     1368    loadfile(handles,fileinput)
     1369end
     1370
     1371% -----------------------------------------------------------------------
     1372% --- Executes on menubar option Import/Grid file: introduce previous grid files
     1373function MenuGridFile_Callback(hObject, eventdata, handles)
     1374% -----------------------------------------------------------------------
     1375inputfile=browse_xml(hObject, eventdata, handles)
     1376listfile=get(handles.coord_files,'string');
     1377if isequal(listfile,{''})
     1378    listfile={inputfile};
     1379else
     1380    listfile=[listfile;{inputfile}]%update the list of coord files
     1381end
     1382set(handles.coord_files,'string',listfile);
     1383
     1384%------------------------------------------------------------------------
     1385function fileinput=browse_xml(hObject, eventdata, handles)
     1386%------------------------------------------------------------------------
     1387fileinput=[];%default
     1388oldfile=''; %default
     1389UserData=get(handles.geometry_calib,'UserData');
     1390if isfield(UserData,'XmlInput')
     1391    oldfile=UserData.XmlInput;
     1392end
     1393[FileName, PathName, filterindex] = uigetfile( ...
     1394       {'*.xml;*.mat', ' (*.xml,*.mat)';
     1395       '*.xml',  '.xml files '; ...
     1396        '*.mat',  '.mat matlab files '}, ...
     1397        'Pick a file',oldfile);
     1398fileinput=[PathName FileName];%complete file name
     1399testblank=findstr(fileinput,' ');%look for blanks
     1400if ~isempty(testblank)
     1401    msgbox_uvmat('ERROR','forbidden input file name or path: no blank character allowed')
     1402    return
     1403end
     1404sizf=size(fileinput);
     1405if (~ischar(fileinput)||~isequal(sizf(1),1)),return;end
     1406UserData.XmlInput=fileinput;
     1407set(handles.geometry_calib,'UserData',UserData)%record current file foer further use of browser
     1408
     1409% -----------------------------------------------------------------------
     1410function loadfile(handles,fileinput)
     1411%------------------------------------------------------------------------
     1412[s,errormsg]=imadoc2struct(fileinput,'GeometryCalib');
     1413GeometryCalib=s.GeometryCalib;
     1414if isempty(GeometryCalib)
     1415    Tabchar={};
     1416    CoordCell={};
     1417    k=0;%default
     1418    f1=1000;
     1419    f2=1000;
     1420    hhuvmat=guidata(findobj(allchild(0),'Name','uvmat'));
     1421    Cx=str2num(get(hhuvmat.npx,'String'))/2;
     1422    Cy=str2num(get(hhuvmat.npy,'String'))/2;
     1423else
     1424    k=GeometryCalib.kappa1 * GeometryCalib.f*GeometryCalib.f;
     1425    f1=GeometryCalib.f*GeometryCalib.sx/GeometryCalib.dpx;
     1426    f2=GeometryCalib.f/GeometryCalib.dpy;
     1427    Cx=GeometryCalib.Cx;
     1428    Cy=GeometryCalib.Cy;
     1429    set(handles.fx,'String',num2str(f1,'%1.1f'))
     1430    set(handles.fy,'String',num2str(f2,'%1.1f'))
     1431    set(handles.k,'String',num2str(k,'%1.4f'))
     1432    set(handles.Cx,'String',num2str(Cx,'%1.1f'))
     1433    set(handles.Cy,'String',num2str(Cy,'%1.1f'))
     1434    calib=reshape(GeometryCalib.PointCoord',[],1);
     1435    for ilist=1:numel(calib)
     1436        CoordCell{ilist}=num2str(calib(ilist));
     1437    end
     1438    CoordCell=reshape(CoordCell,[],5);
     1439    Tabchar=cell2tab(CoordCell,'    |    ');%transform cells into table ready for display
     1440    set(handles.ListCoord,'Value',1)
     1441    set(handles.ListCoord,'String',Tabchar)
     1442    MenuPlot_Callback(handles.geometry_calib, [], handles)
     1443end
     1444Tabchar=[Tabchar;{'......'}];
     1445if isempty(CoordCell)% allow mouse action by default in the absence of input points
     1446    set(handles.edit_append,'Value',1)
     1447    set(handles.edit_append,'BackgroundColor',[1 1 0])
     1448else % does not allow mouse action by default in the presence of input points
     1449    set(handles.edit_append,'Value',0)
     1450    set(handles.edit_append,'BackgroundColor',[0.7 0.7 0.7])
     1451end
     1452
     1453
     1454
     1455
     1456
     1457% --------------------------------------------------------------------
     1458% --- Executes on button press in CLEAR_PTS: clear the list of calibration points
     1459function CLEAR_PTS_Callback(hObject, eventdata, handles)
     1460% --------------------------------------------------------------------
     1461set(handles.ListCoord,'Value',1)% refresh the display of coordinates
     1462set(handles.ListCoord,'String',{'......'})
     1463
  • trunk/src/imadoc2struct.m

    r89 r109  
    11%'imadoc2struct': reads the xml file for image documentation
    22%------------------------------------------------------------------------
    3 % function [s,errormsg]=imadoc2struct(ImaDoc)
     3% function [s,errormsg]=imadoc2struct(ImaDoc,option)
    44%
    55% OUTPUT:
     
    1313% INPUT:
    1414% ImaDoc: full name of the xml input file with head key ImaDoc
    15 
    16 function [s,errormsg]=imadoc2struct(ImaDoc)
    17 
     15% option: ='GeometryCalib': read  the data of GeometryCalib, including source point coordinates
     16
     17function [s,errormsg]=imadoc2struct(ImaDoc,option)
     18if ~exist('option','var')
     19    option='default';
     20end
    1821errormsg=[];%default
    1922s.Heading=[];%default
     
    8891                Dtj=reshape(Dtj'*ones(1,NbDtj),1,NbDtj*numel(Dtj)); %concatene Dti vector NbDti times
    8992                Dtj=[0 Dtj];
    90 %                 Time_val'
    91 %                 ones(1,numel(Dtj))
    92 %                 ones(numel(Time_val'),1)
    93 %                 cumsum(Dtj)
    9493                Time_val=Time_val*ones(1,numel(Dtj))+ones(numel(Time_val),1)*cumsum(Dtj);% produce a time matrix with Dtj
    9594            end
     
    107106        end
    108107    end
    109 %     if size(s.Time,1)==1
    110 %         s.Time=(s.Time)'; %change vector into column
    111 %     end
    112108end
    113109
     
    210206            end         
    211207        end
     208        if strcmp(option,'GeometryCalib')
     209            tsai.PointCoord=get_value(subt,'/GeometryCalib/SourceCalib/PointCoord',[0 0 0 0 0]);%time in TimeUnit
     210        end
    212211        s.GeometryCalib=tsai;
    213212    end
     
    228227               val_read=str2num(get(t,uid_child(icell),'value'));
    229228               if ~isempty(val_read)
    230                    val(icell)=val_read;
     229                   val(icell,:)=val_read;
    231230               end
    232231           end
  • trunk/src/plot_field.m

    r107 r109  
    243243function hdisplay=plot_text(FieldData,hdisplay_in)
    244244%-------------------------------------------------------------------
    245 hdisplay_in
    246245if exist('hdisplay_in','var') && ~isempty(hdisplay_in) && ishandle(hdisplay_in) && isequal(get(hdisplay_in,'Type'),'uicontrol')
    247246    hdisplay=hdisplay_in;
     
    687686        end
    688687          x_label=[Data.ListVarName{ivar_X} '(' x_units ')'];
    689     end       
    690 %     if isfield(Data,'VarAttribute')
    691 %         VarAttribute=Data.VarAttribute;
    692 %     end   
     688    end         
    693689end
    694690
     
    704700    siz=numel(np);
    705701    if siz>3
    706         msgbox_uvmat('ERROR','unrecognized scalar type ')
     702        msgbox_uvmat('ERROR',['unrecognized scalar type: ' num2str(siz) ' dimensions'])
    707703            return
    708704    end
     
    794790        if ~isequal(PlotParam.Scalar.Contours,1) 
    795791            % rescale the grey levels with min and max, put a grey scale colorbar
     792            B=A;
    796793            if BW
    797                 B=A;
    798794                vec=linspace(0,1,255);%define a linear greyscale colormap
    799795                map=[vec' vec' vec'];
    800796                colormap(map);  %grey scale color map
    801797            else
    802                 B=A;
    803798                colormap('default'); % standard faulse colors for div, vort , scalar fields
    804799            end
  • trunk/src/proj_field.m

    r77 r109  
    1212%  .NbPix;
    1313%  .DimName=  names of the matrix dimensions (matlab cell)
    14 %  .DimValue= values of the matricx dimensions (matlab vector, same length as .DimName);
    1514%  .VarName= names of the variables [ProjData.VarName {'A','AMean','AMin','AMax'}];
    1615%  .VarDimNameIndex= dimensions of the variables, indicated by indices in the list .DimName;
     
    3433%         .ListGlobalAttribute: cell listing the names of the global attributes
    3534%        .Att_1,Att_2... : values of the global attributes
    36 %            .ListDimName: cell listing the names of the array dimensions
    37 %               .DimValue: array dimension values (Matlab vector with the same length as .ListDimName
    3835%            .ListVarName: cell listing the names of the variables
    39 %            .VarDimIndex: cell containing the set of dimension indices (in list .ListDimName) for each variable of .ListVarName
    4036%           .VarAttribute: cell of structures s containing names and values of variable attributes (s.name=value) for each variable of .ListVarName
    4137%        .Var1, .Var2....: variables (Matlab arrays) with names listed in .ListVarName
     
    391387%A REVOIR, GENERALISER: UTILISER proj_line
    392388ProjData.NbDim=1;
    393 ProjData.ListDimName={};%name of dimension
    394 ProjData.DimValue=[];%values of dimension (nbre of vectors)
     389%ProjData.ListDimName={};%name of dimension
     390%ProjData.DimValue=[];%values of dimension (nbre of vectors)
    395391ProjData.ListVarName={};
    396392%ProjData.VarDimIndex={};
     
    966962     DY=abs(ObjectData.DY);%mesh of interpolation points
    967963end
    968 if  ~isequal(ProjMode,'projection') & (DX==0|DY==0)
     964if  ~strcmp(ProjMode,'projection') && (DX==0||DY==0)
    969965        errormsg='DX or DY missing';
    970966        display(errormsg)
     
    997993ProjData=proj_heading(FieldData,ObjectData);
    998994ProjData.NbDim=2;
    999 ProjData.ListDimName={};%name of dimension
    1000 ProjData.DimValue=[];%values of dimension (nbre of vectors)
    1001995ProjData.ListVarName={};
    1002996ProjData.VarDimName={};
     
    10451039        DimCell={DimCell};%name of dimensions
    10461040    end
    1047    
     1041
    10481042%case of input fields with unstructured coordinates
    10491043    if testX
     
    11341128        if isequal(ObjectData.ProjMode,'projection')
    11351129            %the list of dimension
    1136             ProjData.ListDimName=[ProjData.ListDimName FieldData.VarDimName(VarIndex(1))];%add the point index to the list of dimensions
    1137             ProjData.DimValue=[ProjData.DimValue length(coord_X)];
     1130            %ProjData.ListDimName=[ProjData.ListDimName FieldData.VarDimName(VarIndex(1))];%add the point index to the list of dimensions
     1131            %ProjData.DimValue=[ProjData.
     1132             %length(coord_X)];
    11381133            nbvar=0;
    11391134            for ivar=VarIndex %transfer variables to the projection plane
     
    11811176            testFF=0;
    11821177            for ivar=VarIndex
    1183                 VarName=FieldData.ListVarName{ivar}; 
     1178                VarName=FieldData.ListVarName{ivar};
    11841179                if ~( ivar==ivar_X || ivar==ivar_Y || ivar==ivar_Z || ivar==ivar_F || ivar==ivar_FF || test_anc(ivar)==1)                 
    11851180                    ivar_new=ivar_new+1;
     
    12201215            end
    12211216        end
    1222 %case of fields defined on a structured  grid
     1217%case of input fields defined on a structured  grid
    12231218    else 
    1224         AYName=FieldData.ListVarName{VarType.coord(1)};
    1225         AXName=FieldData.ListVarName{VarType.coord(2)};
     1219        AYName=FieldData.ListVarName{VarType.coord(1)};%name of input x coordinate (name preserved on projection)
     1220        AXName=FieldData.ListVarName{VarType.coord(2)};%name of input y coordinate (name preserved on projection)
    12261221        eval(['AX=FieldData.' AXName ';'])
    12271222        eval(['AY=FieldData.' AYName ';'])
    1228         VarName=FieldData.ListVarName{VarIndex(1)};
    1229         eval(['DimValue=size(FieldData.' VarName ');'])
     1223        VarName=FieldData.ListVarName{VarIndex(1)};%get the first variable of the cell to get the input matrix dimensions
     1224        eval(['DimValue=size(FieldData.' VarName ');'])%input matrix dimensions
    12301225        ListDimName=FieldData.VarDimName{VarIndex(1)};
    12311226        ProjData.ListVarName=[{AYName} {AXName} ProjData.ListVarName]; %TODO: check if it already exists in Projdata (several cells)
     
    12341229        for idim=1:length(ListDimName)
    12351230            DimName=ListDimName{idim};
    1236             if isequal(DimName,'rgb')|isequal(DimName,'nb_coord')|isequal(DimName,'nb_coord_i')
     1231            if strcmp(DimName,'rgb')||strcmp(DimName,'nb_coord')||strcmp(DimName,'nb_coord_i')
    12371232               nbcolor=DimValue(idim);
    1238                DimIndices(idim)=[];
    12391233               DimValue(idim)=[];
    12401234            end
    12411235            if isequal(DimName,'nb_coord_j')% NOTE: CASE OF TENSOR NOT TREATED
    1242                 DimIndices(idim)=[];
    12431236                DimValue(idim)=[];
    12441237            end
    12451238        end 
    12461239        ind_1=find(DimValue==1);
    1247         DimIndices(ind_1)=[]; %suppress singleton dimensions
    1248 %         indxy=find(DimVarIndex(DimIndices));%select dimension variables (DimIndices non zero)
    1249         nb_dim=length(DimIndices);%number of space dimensions
    12501240        Coord_z=[];
    12511241        Coord_y=[];
    12521242        Coord_x=[];   
    1253    
     1243        nb_dim=numel(size(DimValue));
    12541244        for idim=1:nb_dim %loop on space dimensions
    12551245            test_interp(idim)=0;%test for coordiate interpolation (non regular grid), =0 by default
    1256             test_coord(idim)=0;%test for defined coordinates, =0 by default
    1257             ivar=DimVarIndex(DimIndices(idim));% index of the variable corresponding to the current dimension
    1258             if ~isequal(ivar,0)%  a variable corresponds to the current dimension
    1259                 eval(['Coord{idim}=FieldData.' FieldData.ListVarName{ivar} ';']) ;% position for the first index
    1260                 DCoord=diff(Coord{idim});
    1261                 DCoord_min(idim)=min(DCoord);
    1262                 DCoord_max=max(DCoord);
    1263                 test_direct(idim)=DCoord_max>0;% =1 for increasing values, 0 otherwise
    1264                 test_direct_min=DCoord_min(idim)>0;% =1 for increasing values, 0 otherwise
    1265                 if ~isequal(test_direct(idim),test_direct_min)
    1266                      msgbox_uvmat('ERROR',['non monotonic dimension variable # ' num2str(idim)  ' in proj_field.m'])
     1246            ivar=VarType.coord(idim);% index of the variable corresponding to the current dimension
     1247            if ~isequal(ivar,0)%  a variable corresponds to the dimension #idim
     1248                eval(['Coord{idim}=FieldData.' FieldData.ListVarName{ivar} ';']) ;% coord values for the input field
     1249                if numel(Coord{idim})==2 %input array defined on a regular grid
     1250                   DCoord_min(idim)=(Coord{idim}(2)-Coord{idim}(1))/DimValue(idim);
     1251                else
     1252                    DCoord=diff(Coord{idim});%array of coordinate derivatives for the input field
     1253                    DCoord_min(idim)=min(DCoord);
     1254                    DCoord_max=max(DCoord);
     1255                %    test_direct(idim)=DCoord_max>0;% =1 for increasing values, 0 otherwise
     1256                    if ~isequal(DCoord_max,DCoord_min(idim)>0)
     1257                        msgbox_uvmat('ERROR',['non monotonic dimension variable # ' num2str(idim)  ' in proj_field.m'])
    12671258                                return
    1268                 end               
    1269                 test_interp(idim)=(DCoord_max-DCoord_min(idim))> 0.0001*abs(DCoord_max);% test grid regularity
    1270                 test_coord(idim)=1;
    1271 
    1272             else  % no variable associated with the first dimension, look for variable  attributes Coord_1, _2 or _3
     1259                    end               
     1260                    test_interp(idim)=(DCoord_max-DCoord_min(idim))> 0.0001*abs(DCoord_max);% test grid regularity
     1261                end
     1262                test_direct(idim)=(DCoord_min(idim)>0);
     1263            else  % no variable associated with the  dimension #idim, the coordinate value is set equal to the matrix index by default
    12731264                Coord_i_str=['Coord_' num2str(idim)];
    12741265                DCoord_min(idim)=1;%default
     
    12771268            end
    12781269        end
    1279         if nb_dim==2
    1280             if DY==0
    1281                 DY=abs(DCoord_min(1));
    1282             end
    1283             npY=1+round(abs(Coord{1}(end)-Coord{1}(1))/DY);%nbre of points after interpolation
    1284             npy=1+round(abs(Coord{1}(end)-Coord{1}(1))/abs(DCoord_min(1)));%nbre of points after possible interpolation on a regular grid
    1285             if DX==0
    1286                 DX=abs(DCoord_min(2));
    1287             end
    1288             npX=1+round(abs(Coord{2}(end)-Coord{2}(1))/DX);%nbre of points after interpol 
    1289             npx=1+round(abs(Coord{2}(end)-Coord{2}(1))/abs(DCoord_min(2)));%nbre of points after possible interpolation on a regular grid
    1290             Coord_y=linspace(Coord{1}(1),Coord{1}(end),npY);
    1291             test_direct_y=test_direct(1);
    1292             Coord_x=linspace(Coord{2}(1),Coord{2}(end),npX);
    1293             test_direct_x=test_direct(2);
    1294             DAX=DCoord_min(2);
    1295             DAY=DCoord_min(1);
    1296         elseif nb_dim==3
     1270        if DY==0
     1271            DY=abs(DCoord_min(nb_dim-1));
     1272        end
     1273        npY=1+round(abs(Coord{nb_dim-1}(end)-Coord{nb_dim-1}(1))/DY);%nbre of points after interpol
     1274        if DX==0
     1275            DX=abs(DCoord_min(nb_dim));
     1276        end
     1277        npX=1+round(abs(Coord{nb_dim}(end)-Coord{nb_dim}(1))/DX);%nbre of points after interpol
     1278        for idim=[1:nb_dim]
     1279            if test_interp(idim)
     1280                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
     1281            end
     1282        end       
     1283        Coord_y=linspace(Coord{nb_dim-1}(1),Coord{nb_dim-1}(end),npY);
     1284        test_direct_y=test_direct(nb_dim-1);
     1285        Coord_x=linspace(Coord{nb_dim}(1),Coord{nb_dim}(end),npX);
     1286        test_direct_x=test_direct(nb_dim);
     1287        DAX=DCoord_min(nb_dim);
     1288        DAY=DCoord_min(nb_dim-1);
     1289        if nb_dim==3
    12971290            DZ=abs(DCoord_min(1));
    1298             npz=1+round(abs(Coord{1}(end)-Coord{1}(1))/DZ);%nbre of points after interpolation
    1299             if DY==0
    1300                 DY=abs(DCoord_min(2));
    1301             end
    1302             npY=1+round(abs(Coord{2}(end)-Coord{2}(1))/DY);%nbre of points after interpol
    1303             npy=1+round(abs(Coord{2}(end)-Coord{2}(1))/abs(DCoord_min(2)));%nbre of points before interpol
    1304             if DX==0
    1305                 DX=abs(DCoord_min(3));
    1306             end
    1307             npX=1+round(abs(Coord{3}(end)-Coord{3}(1))/DX);%nbre of points after interpol
    1308             npx=1+round(abs(Coord{3}(end)-Coord{3}(1))/abs(DCoord_min(3)));%nbre of points before interpol
    13091291            Coord_z=linspace(Coord{1}(1),Coord{1}(end),npz);
    13101292            test_direct_z=test_direct(1);
    1311             Coord_y=linspace(Coord{2}(1),Coord{2}(end),npY);
    1312             test_direct_y=test_direct(2);
    1313             Coord_x=linspace(Coord{3}(1),Coord{3}(end),npX);
    1314             test_direct_x=test_direct(3);
    13151293        end 
    13161294        minAX=min(Coord_x);
     
    13341312            YMin=min(ycor_new);
    13351313        end
    1336         DXinit=(maxAX-minAX)/(npx-1);
    1337         DYinit=(maxAY-minAY)/(npy-1);
     1314        DXinit=(maxAX-minAX)/(DimValue(2)-1);
     1315        DYinit=(maxAY-minAY)/(DimValue(1)-1);
    13381316        if DX==0
    13391317            DX=DXinit;
     
    13431321        end
    13441322        npX=floor((XMax-XMin)/DX+1);
    1345         npY=floor((YMax-YMin)/DY+1);    
     1323        npY=floor((YMax-YMin)/DY+1);   
    13461324        if test_direct_y
    13471325            coord_y_proj=linspace(YMin,YMax,npY);%abscissa of the new pixels along the line
     
    13811359            min_ind1=max(min_ind1,1);% deals with margin (bound lower than the first index)
    13821360            min_ind2=max(min_ind2,1);
    1383             max_ind1=min(max_ind1,npy);
    1384             max_ind2=min(max_ind2,npx);
     1361            max_ind1=min(max_ind1,DimValue(1));
     1362            max_ind2=min(max_ind2,DimValue(2));
    13851363            for ivar=VarIndex
    13861364                VarName=FieldData.ListVarName{ivar};
    13871365                ProjData.ListVarName=[ProjData.ListVarName VarName];
    1388                 ProjData.VarDimIndex=[ProjData.VarDimIndex [nb_dim-1 nb_dim]];
     1366                ProjData.VarDimName=[ProjData.VarDimName {DimCell}];
    13891367                if length(FieldData.VarAttribute)>=ivar
    13901368                    ProjData.VarAttribute{length(ProjData.ListVarName)}=FieldData.VarAttribute{ivar};
     
    13921370                eval(['ProjData.' VarName '=FieldData.' VarName '(min_ind1:max_ind1,min_ind2:max_ind2) ;']);
    13931371            end         
    1394         else
    1395         % case with rotation and/or interpolation
     1372        else       % case with rotation and/or interpolation
    13961373            if isempty(Coord_z) %2D case
    13971374                [X,Y]=meshgrid(coord_x_proj,coord_y_proj);%grid in the new coordinates
     
    14021379                XIMA=reshape(round(XIMA),1,npX*npY);%indices reorganized in 'line'
    14031380                YIMA=reshape(round(YIMA),1,npX*npY);
    1404                 flagin=XIMA>=1 & XIMA<=npx & YIMA >=1 & YIMA<=npy;%flagin=1 inside the original image
     1381                flagin=XIMA>=1 & XIMA<=DimValue(2) & YIMA >=1 & YIMA<=DimValue(1);%flagin=1 inside the original image
    14051382                if isequal(ObjectData.ProjMode,'filter')
    14061383                    npx_filter=ceil(abs(DX/DAX));
     
    14111388                    test_filter=0;
    14121389                end
     1390                eval(['ProjData.' AYName '=[coord_y_proj(1) coord_y_proj(end)];']) %record the new (projected ) y coordinates
     1391                eval(['ProjData.' AXName '=[coord_x_proj(1) coord_x_proj(end)];']) %record the new (projected ) x coordinates
    14131392                for ivar=VarIndex
    1414                     VarName=FieldData.ListVarName{ivar} ;
    1415                     if test_interp(1) | test_interp(2)%interpolate on a regular grid        
    1416                           eval(['FieldData.' VarName '=interp2(Coord{2},Coord{1},FieldData.' VarName ',Coord_x,Coord_y'');']) %TO TEST
     1393                    VarName=FieldData.ListVarName{ivar};
     1394                    if test_interp(1) || test_interp(2)%interpolate on a regular grid       
     1395                          eval(['ProjData.' VarName '=interp2(Coord{2},Coord{1},FieldData.' VarName ',Coord_x,Coord_y'');']) %TO TEST
    14171396                    end
    14181397                    %filter the field (image) if option 'filter' is used
    14191398                    if test_filter 
    14201399                         Aclass=class(FieldData.A);
    1421                          eval(['FieldData.' VarName '=filter2(Mfilter,FieldData.' VarName ',''valid'');'])
     1400                         eval(['ProjData.' VarName '=filter2(Mfilter,FieldData.' VarName ',''valid'');'])
    14221401                         if ~isequal(Aclass,'double')
    1423                              eval(['FieldData.' VarName '=' Aclass '(FieldData.' VarName ');'])%revert to integer values
     1402                             eval(['ProjData.' VarName '=' Aclass '(FieldData.' VarName ');'])%revert to integer values
    14241403                         end
    14251404                    end
    1426                     eval(['vec_A=reshape(FieldData.' VarName ',npx*npy,nbcolor);'])%put the original image in line             
     1405                    eval(['vec_A=reshape(FieldData.' VarName ',[],nbcolor);'])%put the original image in line             
    14271406                    ind_in=find(flagin);
    14281407                    ind_out=find(~flagin);
    1429                     ICOMB=(XIMA-1)*npy+YIMA;
     1408                    ICOMB=(XIMA-1)*DimValue(1)+YIMA;
    14301409                    ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A
    14311410                    vec_B(ind_in,[1:nbcolor])=vec_A(ICOMB,:);
     
    14331412                        vec_B(ind_out,icolor)=zeros(size(ind_out));
    14341413                    end
    1435                     ProjData.ListVarName=[ProjData.ListVarName VarName];                 
    1436                     if length(FieldData.VarAttribute)>=ivar
     1414                    ProjData.ListVarName=[ProjData.ListVarName VarName];
     1415                    ProjData.VarDimName=[ProjData.VarDimName {DimCell}];
     1416                    if isfield(FieldData,'VarAttribute')&&length(FieldData.VarAttribute)>=ivar
    14371417                        ProjData.VarAttribute{length(ProjData.ListVarName)+nbcoord}=FieldData.VarAttribute{ivar};
    14381418                    end     
     
    14411421                ProjData.FF=reshape(~flagin,npY,npX);%false flag A FAIRE: tenir compte d'un flga antérieur 
    14421422                ProjData.ListVarName=[ProjData.ListVarName 'FF'];
     1423                ProjData.VarDimName=[ProjData.VarDimName {DimCell}];
    14431424                ProjData.VarAttribute{length(ProjData.ListVarName)}.Role='errorflag';
    14441425            else %3D case
     
    14481429                    iz=iz_sup(1);
    14491430                    if iz>=1 & iz<=npz
    1450                         ProjData.ListDimName=[ProjData.ListDimName ListDimName(2:end)];
    1451                         ProjData.DimValue=[ProjData.DimValue npY npX];
     1431                        %ProjData.ListDimName=[ProjData.ListDimName ListDimName(2:end)];
     1432                        %ProjData.DimValue=[ProjData.DimValue npY npX];
    14521433                        for ivar=VarIndex
    14531434                            VarName=FieldData.ListVarName{ivar};
     
    15821563ProjData=proj_heading(FieldData,ObjectData);
    15831564ProjData.NbDim=3;
    1584 ProjData.ListDimName={};%name of dimension
    1585 ProjData.DimValue=[];%values of dimension (nbre of vectors)
     1565%ProjData.ListDimName={};%name of dimension
     1566%ProjData.DimValue=[];%values of dimension (nbre of vectors)
    15861567ProjData.ListVarName={};
    15871568ProjData.VarDimName={};
     
    17171698        if isequal(ObjectData.ProjMode,'projection')
    17181699            %the list of dimension
    1719             ProjData.ListDimName=[ProjData.ListDimName FieldData.VarDimName(VarIndex(1))];%add the point index to the list of dimensions
    1720             ProjData.DimValue=[ProjData.DimValue length(coord_X)];
     1700            %ProjData.ListDimName=[ProjData.ListDimName FieldData.VarDimName(VarIndex(1))];%add the point index to the list of dimensions
     1701            %ProjData.DimValue=[ProjData.DimValue length(coord_X)];
    17211702            nbvar=0;
    17221703            for ivar=VarIndex %transfer variables to the projection plane
     
    18431824        for idim=1:nb_dim %loop on space dimensions
    18441825            test_interp(idim)=0;%test for coordiate interpolation (non regular grid), =0 by default
    1845             test_coord(idim)=0;%test for defined coordinates, =0 by default
    18461826            ivar=DimVarIndex(DimIndices(idim));% index of the variable corresponding to the current dimension
    18471827            if ~isequal(ivar,0)%  a variable corresponds to the current dimension
     
    18571837                end               
    18581838                test_interp(idim)=(DCoord_max-DCoord_min(idim))> 0.0001*abs(DCoord_max);% test grid regularity
    1859                 test_coord(idim)=1;
    1860 
    18611839            else  % no variable associated with the first dimension, look for variable  attributes Coord_1, _2 or _3
    18621840                Coord_i_str=['Coord_' num2str(idim)];
     
    18711849            end
    18721850            npY=1+round(abs(Coord{1}(end)-Coord{1}(1))/DY);%nbre of points after interpolation
    1873             npy=1+round(abs(Coord{1}(end)-Coord{1}(1))/abs(DCoord_min(1)));%nbre of points after possible interpolation on a regular grid
     1851            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
    18741852            if DX==0
    18751853                DX=abs(DCoord_min(2));
    18761854            end
    18771855            npX=1+round(abs(Coord{2}(end)-Coord{2}(1))/DX);%nbre of points after interpol 
    1878             npx=1+round(abs(Coord{2}(end)-Coord{2}(1))/abs(DCoord_min(2)));%nbre of points after possible interpolation on a regular grid
     1856            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
    18791857            Coord_y=linspace(Coord{1}(1),Coord{1}(end),npY);
    18801858            test_direct_y=test_direct(1);
     
    18901868            end
    18911869            npY=1+round(abs(Coord{2}(end)-Coord{2}(1))/DY);%nbre of points after interpol
    1892             npy=1+round(abs(Coord{2}(end)-Coord{2}(1))/abs(DCoord_min(2)));%nbre of points before interpol
     1870            DimValue(1)=1+round(abs(Coord{2}(end)-Coord{2}(1))/abs(DCoord_min(2)));%nbre of points before interpol
    18931871            if DX==0
    18941872                DX=abs(DCoord_min(3));
    18951873            end
    18961874            npX=1+round(abs(Coord{3}(end)-Coord{3}(1))/DX);%nbre of points after interpol
    1897             npx=1+round(abs(Coord{3}(end)-Coord{3}(1))/abs(DCoord_min(3)));%nbre of points before interpol
     1875            DimValue(2)=1+round(abs(Coord{3}(end)-Coord{3}(1))/abs(DCoord_min(3)));%nbre of points before interpol
    18981876            Coord_z=linspace(Coord{1}(1),Coord{1}(end),npz);
    18991877            test_direct_z=test_direct(1);
     
    19231901            YMin=min(ycor_new);
    19241902        end
    1925         DXinit=(maxAX-minAX)/(npx-1);
    1926         DYinit=(maxAY-minAY)/(npy-1);
     1903        DXinit=(maxAX-minAX)/(DimValue(2)-1);
     1904        DYinit=(maxAY-minAY)/(DimValue(1)-1);
    19271905        if DX==0
    19281906            DX=DXinit;
     
    19701948            min_ind1=max(min_ind1,1);% deals with margin (bound lower than the first index)
    19711949            min_ind2=max(min_ind2,1);
    1972             max_ind1=min(max_ind1,npy);
    1973             max_ind2=min(max_ind2,npx);
     1950            max_ind1=min(max_ind1,DimValue(1));
     1951            max_ind2=min(max_ind2,DimValue(2));
    19741952            for ivar=VarIndex
    19751953                VarName=FieldData.ListVarName{ivar};
    19761954                ProjData.ListVarName=[ProjData.ListVarName VarName];
    1977                 ProjData.VarDimIndex=[ProjData.VarDimIndex [nb_dim-1 nb_dim]];
     1955                %ProjData.VarDimIndex=[ProjData.VarDimIndex [nb_dim-1 nb_dim]];
    19781956                if length(FieldData.VarAttribute)>=ivar
    19791957                    ProjData.VarAttribute{length(ProjData.ListVarName)}=FieldData.VarAttribute{ivar};
     
    19911969                XIMA=reshape(round(XIMA),1,npX*npY);%indices reorganized in 'line'
    19921970                YIMA=reshape(round(YIMA),1,npX*npY);
    1993                 flagin=XIMA>=1 & XIMA<=npx & YIMA >=1 & YIMA<=npy;%flagin=1 inside the original image
     1971                flagin=XIMA>=1 & XIMA<=DimValue(2) & YIMA >=1 & YIMA<=DimValue(1);%flagin=1 inside the original image
    19941972                if isequal(ObjectData.ProjMode,'filter')
    19951973                    npx_filter=ceil(abs(DX/DAX));
     
    20131991                         end
    20141992                    end
    2015                     eval(['vec_A=reshape(FieldData.' VarName ',npx*npy,nbcolor);'])%put the original image in line             
     1993                    eval(['vec_A=reshape(FieldData.' VarName ',[],nbcolor);'])%put the original image in line             
    20161994                    ind_in=find(flagin);
    20171995                    ind_out=find(~flagin);
    2018                     ICOMB=(XIMA-1)*npy+YIMA;
     1996                    ICOMB=(XIMA-1)*DimValue(1)+YIMA;
    20191997                    ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A
    20201998                    vec_B(ind_in,[1:nbcolor])=vec_A(ICOMB,:);
     
    20372015                    iz=iz_sup(1);
    20382016                    if iz>=1 & iz<=npz
    2039                         ProjData.ListDimName=[ProjData.ListDimName ListDimName(2:end)];
    2040                         ProjData.DimValue=[ProjData.DimValue npY npX];
     2017                        %ProjData.ListDimName=[ProjData.ListDimName ListDimName(2:end)];
     2018                        %ProjData.DimValue=[ProjData.DimValue npY npX];
    20412019                        for ivar=VarIndex
    20422020                            VarName=FieldData.ListVarName{ivar};
  • trunk/src/px_XYZ.m

    r92 r109  
     1%'px_XYZ': transform physical to image coordinates.
     2%------------------------------------------------------------------------
     3%[X,Y]=px_XYZ(Calib,Xphys,Yphys,Zphys)
     4%------------------------------------------------------------------------           
     5% OUTPUT:
     6% [X,Y]: image coordinates(in pixels)
     7%------------------------------------------------------------------------
     8% INPUT:
     9% Calib: structure containing calibration parameters
     10% Xphys,Yphys,Zphys; vectors of physical coordinates for a set of points
    111
    212function [X,Y]=px_XYZ(Calib,Xphys,Yphys,Zphys)
  • trunk/src/series/merge_proj.m

    r106 r109  
    5959nbview=length(Series.RootFile);%number of views (file series to merge)
    6060nbfield=size(num_i1{1},1)*size(num_i1{1},2);%number of fields in the time series
    61 transform=Series.CoordType; %  field transform function
     61% transform=Series.CoordType; %  field transform function
    6262hhh=which('mmreader');
    6363for iview=1:nbview
     
    7272
    7373%Calibration data and timing: read the ImaDoc files
    74 mode=''; %default
     74% mode=''; %default
    7575timecell={};
    7676itime=0;
     
    152152
    153153% Field and velocity type (the same for all views)
     154FieldName='';
     155if strcmp(get(hseries.FieldMenu,'Visible'),'on')
    154156Field_str=get(hseries.FieldMenu,'String');
    155157val=get(hseries.FieldMenu,'Value');
     
    163165    hget_field=findobj(allchild(0),'Name','get_field');%find the get_field... GUI
    164166    SubField=get_field('read_get_field',hObject,eventdata,hget_field); %read the names of the variables to plot in the get_field GUI
     167end
    165168end
    166169%detect whether all the files are 'images' or 'netcdf'
     
    238241    if isequal(stopstate,'queue')% enable STOP command from the 'series' interface
    239242         update_waitbar(hseries.waitbar,WaitbarPos,ifile/nbfield)
    240          Amerge=0;
     243%          Amerge=0;
    241244         
    242245         %----------LOOP ON VIEWS----------------------
     
    258261                end % TODO: introduce ListVarName
    259262                npxy=size(Field{iview}.A);
     263                Field{iview}.ListVarName={'AX','AY','A'};
     264                Field{iview}.VarDimName={'AX','AY',{'AY','AX'}};
    260265                Field{iview}.AX=[0.5 npxy(2)-0.5]; % coordinates of the first and last pixel centers
    261266                Field{iview}.AY=[npxy(1)-0.5 0.5];
    262267                Field{iview}.CoordType='px';
    263268                Field{iview}.AName='image';
     269                timeread(iview)=0;
    264270            else
    265271                if testcivx
     
    287293            end
    288294            if testcivx
    289                     Field{iview}=calc_field(FieldName,Field{iview});
     295                Field{iview}=calc_field(FieldName,Field{iview});
    290296            end
    291297
     
    393399    return
    394400end
    395 for iview=1:nbview
    396     if ~isequal(MergeData.ListDimName,Data{iview}.ListDimName)
    397         error=1;
    398     end
    399     if ~isequal(MergeData.ListVarName,Data{iview}.ListVarName)
    400         error=1;
    401     end
    402 end
    403 if error
    404     MergeData.Txt='ERROR: attempt at merging fields of incompatible type';
    405     return
    406 end
    407 
    408 
    409 
    410 
     401% for iview=1:nbview
     402%     if ~isequal(MergeData.ListDimName,Data{iview}.ListDimName)
     403%         error=1;
     404%     end
     405%     if ~isequal(MergeData.ListVarName,Data{iview}.ListVarName)
     406%         error=1;
     407%     end
     408% end
     409% if error
     410%     MergeData.Txt='ERROR: attempt at merging fields of incompatible type';
     411%     return
     412% end
     413%
    411414%group the variables (fields of 'FieldData') in cells of variables with the same dimensions
    412415[CellVarIndex,NbDim,VarTypeCell]=find_field_indices(Data{1});
     
    414417% CellVarIndex=cells of variable index arrays
    415418ivar_new=0; % index of the current variable in the projected field
    416 icoord=0;
    417419for icell=1:length(CellVarIndex)
    418420    if NbDim(icell)==1
  • trunk/src/sub_field.m

    r89 r109  
    102102   %check coincidence in positions
    103103   %unstructured coordinates
    104    if testX
     104   if testX     
     105       if testU % vector fields
     106            U_1_Name=Field_1.ListVarName{VarType_1.vector_x};
     107            V_1_Name=Field_1.ListVarName{VarType_1.vector_y};
     108            eval(['vec_U_1=Field_1.' U_1_Name ';'])
     109            eval(['vec_V_1=Field_1.' V_1_Name ';'])
     110            nbpoints_x_1=size(vec_U_1,2);
     111            nbpoints_y_1=size(vec_U_1,1);
     112            vec_U_1=reshape(vec_U_1,nbpoints_x_1*nbpoints_y_1,1);
     113            vec_V_1=reshape(vec_V_1,nbpoints_x_1*nbpoints_y_1,1);
     114            if testfalse_1
     115                vec_U_1=vec_U_1(indsel);
     116                vec_V_1=vec_V_1(indsel);
     117            end           
     118       else %(~testU && ~testU_1)
     119           A_1_Name=Field_1.ListVarName{ivar_C_1};
     120           eval(['vec_A_1=Field_1.' A_1_Name ';'])   
     121           nbpoints_x_1=size(vec_A_1,2);
     122           nbpoints_y_1=size(vec_A_1,1);%TODO: use a faster interpolation method for a regular grid (size(x)=2)
     123           vec_A_1=reshape(vec_A_1,nbpoints_x_1*nbpoints_y_1,1);
     124           if testfalse_1
     125                vec_A_1=vec_A_1(indsel);
     126           end
     127       end
    105128       XName=Field.ListVarName{VarType.coord_x};
    106129       YName=Field.ListVarName{VarType.coord_y};
     
    108131       eval(['vec_Y=Field.' YName ';'])
    109132       nbpoints=numel(vec_X);
    110        vec_X=reshape(vec_X,1,nbpoints);
    111        vec_Y=reshape(vec_Y,1,nbpoints);
     133       vec_X=reshape(vec_X,nbpoints,1);
     134       vec_Y=reshape(vec_Y,nbpoints,1);
    112135       if testX_1 %unstructured coordinates for the second field
    113136            X_1_Name=Field_1.ListVarName{VarType_1.coord_x};
     
    115138            eval(['vec_X_1=Field_1.' X_1_Name ';'])
    116139            eval(['vec_Y_1=Field_1.' Y_1_Name ';'])
    117             nbpoints_1=numel(vec_X_1);
     140
    118141       else   %structured coordinates for the second field
    119142           y_1_Name=Field_1.ListVarName{VarType_1.coord(1)};
    120143           x_1_Name=Field_1.ListVarName{VarType_1.coord(2)};
    121144           eval(['y_1=Field_1.' y_1_Name ';'])
    122            eval(['x_1=Field_1.' x_1_Name ';']) 
    123            npxy(1)=numel(y_1);
    124            npxy(2)=numel(x_1);
    125            nbpoints_1=npxy(1)*npxy(2);
     145           eval(['x_1=Field_1.' x_1_Name ';'])
     146           if isequal(numel(x_1),2) 
     147               x_1=linspace(x_1(1),x_1(2),nbpoints_x_1);
     148           end
     149           if isequal(numel(y_1),2) 
     150               y_1=linspace(y_1(1),y_1(2),nbpoints_y_1);
     151           end
    126152           [vec_X_1,vec_Y_1]=meshgrid(x_1,y_1);
    127153       end
    128        vec_X_1=reshape(vec_X_1,1,nbpoints_1);
    129        vec_Y_1=reshape(vec_Y_1,1,nbpoints_1);
     154       vec_X_1=reshape(vec_X_1,nbpoints_x_1*nbpoints_y_1,1);
     155       vec_Y_1=reshape(vec_Y_1,nbpoints_x_1*nbpoints_y_1,1);
    130156       if testfalse_1
    131157           FFName_1=Field_1.ListVarName{VarType_1.errorflag};         
    132158           eval(['vec_FF_1=Field_1.' FFName_1 ';'])
    133            vec_FF_1=reshape(vec_FF_1,1,nbpoints_1);
     159           vec_FF_1=reshape(vec_FF_1,nbpoints_1,1);
    134160           indsel=find(~vec_FF_1);
    135161           vec_X_1=vec_X_1(indsel);
    136162           vec_Y_1=vec_Y_1(indsel);
    137        end
    138        if testU % vector fields
    139             U_1_Name=Field_1.ListVarName{VarType_1.vector_x};
    140             V_1_Name=Field_1.ListVarName{VarType_1.vector_y};
    141             eval(['vec_U_1=Field_1.' U_1_Name ';'])
    142             eval(['vec_V_1=Field_1.' V_1_Name ';'])
    143             vec_U_1=reshape(vec_U_1,1,nbpoints_1);
    144             vec_V_1=reshape(vec_V_1,1,nbpoints_1);
    145             if testfalse_1
    146                 vec_U_1=vec_U_1(indsel);
    147                 vec_V_1=vec_V_1(indsel);
    148             end           
    149        else
    150            A_1_Name=Field_1.ListVarName{ivar_C_1};
    151            eval(['vec_A_1=Field_1.' A_1_Name ';'])
    152            vec_A_1=reshape(vec_A_1,1,nbpoints_1);
    153            if testfalse_1
    154                 vec_A_1=vec_A_1(indsel);
    155            end
    156163       end
    157164       if ~isequal(vec_X_1,vec_X) && ~isequal(vec_Y_1,vec_Y) % if the unstructured positions are not the same
     
    160167               vec_V_1=griddata_uvmat(vec_X_1,vec_Y_1,vec_V_1,vec_X,vec_Y);  %interpolate vectors in the second field   
    161168           else
    162                vec_A_1=griddata_uvmat(vec_X_1,vec_Y_1,vec_A_1,vec_X,vec_Y);  %interpolate vectors in the second field
     169               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
    163170           end
    164171       end
     
    168175           eval(['vec_U=Field.' UName ';'])
    169176           eval(['vec_V=Field.' VName ';'])       
    170            vec_U=reshape(vec_U,1,numel(vec_U));
    171            vec_V=reshape(vec_V,1,numel(vec_V));
     177           vec_U=reshape(vec_U,numel(vec_U),1);
     178           vec_V=reshape(vec_V,numel(vec_V),1);
    172179           eval(['SubData.' UName '=vec_U-vec_U_1;'])
    173180           eval(['SubData.' VName '=vec_V-vec_V_1;'])
    174181       else
    175182           AName=Field.ListVarName{ivar_C};
     183           size(Field.vort)
    176184           eval(['SubData.' AName '=Field.' AName '-vec_A_1;'])
    177185       end
Note: See TracChangeset for help on using the changeset viewer.