Changeset 1078 for trunk/src

Timestamp:

Mar 30, 2020, 3:48:19 PM (5 years ago)

Author:

sommeria

Message:

various updates

Location:

trunk/src

Files:

• angle2normal.m (added)
• calc_field_interp.m (modified) (3 diffs)
• find_field_bounds.m (modified) (1 diff)
• find_field_cells.m (modified) (3 diffs)
• mouse_down.m (modified) (1 diff)
• mouse_motion.m (modified) (1 diff)
• phys_XYZ.m (modified) (3 diffs)
• phys_ima.m (modified) (5 diffs)
• proj_field.m (modified) (4 diffs)
• rename_indexing.m (added)
• round_uvmat.m (added)
• series.fig (modified) (previous)
• series.m (modified) (14 diffs)
• series/merge_proj.m (modified) (1 diff)
• transform_field/phys_polar.m (modified) (6 diffs)
• uvmat.m (modified) (13 diffs)

trunk/src/calc_field_interp.m

@@ -54 +54 @@
         if isempty(ivar)% the requested variable does not exist
             check_skipped(ilist)=1; %variable not found
-        elseif isempty(find(strcmp(FieldName{ilist},InputVarList), 1));% the variable exists and has not been already selected
+        elseif isempty(find(strcmp(FieldName{ilist},InputVarList), 1))% the variable exists and has not been already selected
             if exist('XI','var')&& isfield(Data.VarAttribute{ivar},'Role') &&...
                     (strcmp(Data.VarAttribute{ivar}.Role,'ancillary')||strcmp(Data.VarAttribute{ivar}.Role,'warnflag')||strcmp(Data.VarAttribute{ivar}.Role,'errorflag'))
@@ -86 +86 @@
             UName{ilist}=r.UName;
             VName{ilist}=r.VName;
-            if isempty(find(strcmp(r.UName,InputVarList)));
+            if isempty(find(strcmp(r.UName,InputVarList)))
                 InputVarList=[InputVarList UName{ilist}]; %the variable is added to the list if it is not already in the list
             end
-            if isempty(find(strcmp(r.VName,InputVarList), 1));
+            if isempty(find(strcmp(r.VName,InputVarList), 1))
                 InputVarList=[InputVarList VName{ilist}]; %the variable is added to the list if it is not already in the list
             end
@@ -162 +162 @@
 end
 
-%% put an error flag to indicate NaN data
-% if exist('XI','var')&&~isempty(VarVal)
-%     nbvar=numel(VarVal);
-%     ListVarName{nbvar+1}='FF';
-%     VarVal{nbvar+1}=isnan(VarVal{nbvar});
-%     VarAttribute{nbvar+1}.Role='errorflag';
-% end
 
 

trunk/src/find_field_bounds.m

@@ -134 +134 @@
 end
 % adjust the mesh to a value 1, 2 , 5 *10^n
-ord=10^(floor(log10(Mesh)));%order of magnitude
-if Mesh/ord>=5
-    FieldOut.CoordMesh=5*ord;
-elseif Mesh/ord>=2
-    FieldOut.CoordMesh=2*ord;
-else
-    FieldOut.CoordMesh=ord;
-end
+FieldOut.CoordMesh = round_uvmat(Mesh);
+
+

trunk/src/find_field_cells.m

@@ -294 +294 @@
             end
         else
-            if isfield(CellInfo{icell},'XName'); % only one coordinate x, switch vector field to 1D plot
+            if isfield(CellInfo{icell},'XName') % only one coordinate x, switch vector field to 1D plot
                 for ind=1:numel(CellInfo{icell}.VarIndex)
                     Role{CellInfo{icell}.VarIndex(ind)}='coord_y';
@@ -502 +502 @@
 for icell=1:numel(CellInfo)
     if isfield(CellInfo{icell},'VarIndex')
+        check_fieldname=0;
         VarIndex=CellInfo{icell}.VarIndex;
         for ivar=VarIndex
@@ -514 +515 @@
             if ~isempty(FieldName{ivar})
                 CellInfo{icell}.FieldName=FieldName{ivar};
+                check_fieldname=1;
             end
             if CheckSub(ivar)==1
                 CellInfo{icell}.CheckSub=1;
+            end
+        end
+        if ~check_fieldname% default FieldName
+            if isfield(CellInfo{icell},'VarIndex_vector_x')&& isfield(CellInfo{icell},'VarIndex_vector_y')
+                UName=Data.ListVarName{CellInfo{icell}.VarIndex_vector_x};
+                VName=Data.ListVarName{CellInfo{icell}.VarIndex_vector_y};        
+                CellInfo{icell}.FieldName=['vec(' UName ',' VName ')'];
             end
         end

trunk/src/mouse_down.m

@@ -390 +390 @@
         IndexObj=IndexObj+1;%start new object
         ObjectData.Coord=[];
-        ObjectNameNew=ObjectData.Name;
-        if isempty(ObjectNameNew)
-            ObjectNameNew=ObjectData.Type;
+        ObjectName=ObjectData.Name;
+        if isempty(ObjectName)
+            ObjectName=ObjectData.Type;
         end
         % add an index to the object name if the proposed name already exists
-        vers=0;% index of the name
+        %vers=0;% index of the name
         ListObject=get(hhuvmat.ListObject,'String');
-        detectname=1;
-        while ~isempty(detectname)
-            detectname=find(strcmp(ObjectNameNew,ListObject),1);%test the existence of the proposed name in the list
-            if detectname% if the object name already exists
-                indstr=regexp(ObjectNameNew,'\D');
-                if indstr(end)<length(ObjectNameNew) %object name ends by a number
-                    vers=str2double(ObjectNameNew(indstr(end)+1:end))+1;
-                    ObjectNameNew=[ObjectNameNew(1:indstr(end)) num2str(vers)];
-                else
-                    vers=vers+1;
-                    ObjectNameNew=[ObjectNameNew(1:indstr(end)) '_' num2str(vers)];
-                end
-            end
-        end
-        ObjectName=ObjectNameNew;
+        ObjectName=rename_indexing(ObjectName,ListObject);% add an index or upgrade it if the object naem already exists
         set(sethandles.Name,'String',ObjectName)% display the default name in set_object
         ListObject=[ListObject;{ObjectName}];

trunk/src/mouse_motion.m

@@ -220 +220 @@
                 pos=[xy(1,1) xy(1,2) 0];%coordinates on the graph
                 if isfield(Field,'ProjObjectAngle')&&~isequal(Field.ProjObjectAngle,[0 0 0])
-                    om=norm(Field.ProjObjectAngle);%norm of rotation angle in radians
-                    OmAxis=Field.ProjObjectAngle/om; %unit vector marking the rotation axis
-                    cos_om=cos(pi*om/180);
-                    sin_om=sin(pi*om/180);
-                    pos=[xy(1,1) xy(1,2) 0];
-                    %pos=cos_om*pos+sin_om*cross(OmAxis,pos)+(1-cos_om)*(OmAxis*pos')*OmAxis;
-                    coeff=OmAxis(3)*(1-cos_om);
-                    norm_plane(1)=OmAxis(1)*coeff+OmAxis(2)*sin_om;
-                    norm_plane(2)=OmAxis(2)*coeff-OmAxis(1)*sin_om;
-                    norm_plane(3)=OmAxis(3)*coeff+cos_om;
-                    %Z0=norm_plane*Field.ProjObjectCoord'/norm_plane(3);
-                    pos(3)=-(norm_plane(1)*pos(1)+norm_plane(2)*pos(2))/norm_plane(3);                               
-                end
-                pos=pos+Field.ProjObjectCoord;
+                    norm_plane=angle2normal(Field.ProjObjectAngle);
+                    pos(3)=-(norm_plane(1)*(pos(1)-Field.ProjObjectCoord(1))+norm_plane(2)*(pos(2)-Field.ProjObjectCoord(2)))/norm_plane(3);                               
+                end
+                pos(3)=pos(3)+Field.ProjObjectCoord(3);
                 text_displ_3=[text_displ_3 'x,y,z=' num2str(pos,4)];
             end

trunk/src/phys_XYZ.m

@@ -1 +1 @@
 %------------------------------------------------------------------------
 %'phys_XYZ':transforms image (px) to real world (phys) coordinates using geometric calibration parameters
-% function [Xphys,Yphys]=phys_XYZ(Calib,X,Y,Z)
+% function [Xphys,Yphys,Zphys]=phys_XYZ(Calib,X,Y,Zindex)
 %
 %OUTPUT:
@@ -30 +30 @@
 testangle=0;
 test_refraction=0;
+Zphys=0; %default output
 if exist('Zindex','var')&& isequal(Zindex,round(Zindex))&& Zindex>0 && isfield(Calib,'SliceCoord')&&size(Calib.SliceCoord,1)>=Zindex
-    if isfield(Calib, 'SliceAngle') && ~isequal(Calib.SliceAngle,[0 0 0]) && ~isequal(Calib.SliceAngle(Zindex,:),[0 0 0])
+    if isfield(Calib, 'SliceAngle') && size(Calib.SliceAngle,1)>=Zindex && ~isequal(Calib.SliceAngle(Zindex,:),[0 0 0])
         testangle=1;
-        om=norm(Calib.SliceAngle(Zindex,:));%norm of rotation angle in radians
-        OmAxis=Calib.SliceAngle(Zindex,:)/om; %unit vector marking the rotation axis
-        cos_om=cos(pi*om/180);
-        sin_om=sin(pi*om/180);
-        coeff=OmAxis(3)*(1-cos_om);
-        norm_plane(1)=OmAxis(1)*coeff+OmAxis(2)*sin_om;
-        norm_plane(2)=OmAxis(2)*coeff-OmAxis(1)*sin_om;
-        norm_plane(3)=OmAxis(3)*coeff+cos_om;
-%         Z0=norm_plane*Calib.SliceCoord(Zindex,:)'/norm_plane(3);
+        norm_plane=angle2normal(Calib.SliceAngle(Zindex,:));
     end
     Z0=Calib.SliceCoord(Zindex,3);%horizontal plane z=cte
-%     end
     Z0virt=Z0;
     if isfield(Calib,'InterfaceCoord') && isfield(Calib,'RefractionIndex')
@@ -134 +126 @@
 end
 
-%'px_XYZ': transform phys coordinates to image coordinates (px)
-%
-% OUTPUT:
-% X,Y: array of coordinates in the image cooresponding to the input physical positions 
-%                    (origin at lower leftcorner, unit=pixel)
-
-% INPUT:
-% Calib: structure containing the calibration parameters (read from the ImaDoc .xml file)
-% Xphys, Yphys: array of x,y physical coordinates
-% [Z0]: corresponding array of z physical coordinates (0 by default)

trunk/src/phys_ima.m

@@ -22 +22 @@
     XmlData={XmlData};
 end
+
 for icell=1:numel(A)
     siz=size(A{icell});
@@ -27 +28 @@
     npy=[npy siz(1)];
     Calib=XmlData{icell}.GeometryCalib;
-    xima=[0.5 siz(2)-0.5 0.5 siz(2)-0.5];%image coordinates of corners
-    yima=[0.5 0.5 siz(1)-0.5 siz(1)-0.5];
-    [xcorner_new,ycorner_new]=phys_XYZ(Calib,xima,yima,ZIndex);%corresponding physical coordinates
+    coord_x=[0.5 siz(2)-0.5];
+    coord_y=[0.5 siz(1)-0.5];
+    x_edge=[linspace(coord_x(1),coord_x(end),npx(icell)) coord_x(end)*ones(1,npy(icell))...
+        linspace(coord_x(end),coord_x(1),npx(icell)) coord_x(1)*ones(1,npy(icell))];%x coordinates of the image edge(four sides)
+    y_edge=[coord_y(1)*ones(1,npx(icell)) linspace(coord_y(1),coord_y(end),npy(icell))...
+        coord_y(end)*ones(1,npx(icell)) linspace(coord_y(end),coord_y(1),npy(icell))];%y coordinates of the image edge(four sides)
+    [xcorner_new,ycorner_new]=phys_XYZ(Calib,x_edge,y_edge,ZIndex);%corresponding physical coordinates
     dx(icell)=(max(xcorner_new)-min(xcorner_new))/(siz(2)-1);
     dy(icell)=(max(ycorner_new)-min(ycorner_new))/(siz(1)-1);
@@ -40 +45 @@
 Rangy(1)=max(ycorner);
 test_multi=(max(npx)~=min(npx)) || (max(npy)~=min(npy)); %different image lengths
-% npX=1+round((Rangx(2)-Rangx(1))/min(dx));% nbre of pixels in the new image (use the finest resolution min(dx) in the set of images)
-% npY=1+round((Rangy(1)-Rangy(2))/min(dy));
 
 npX=1+round((Rangx(2)-Rangx(1))/max(dx));% nbre of pixels in the new image (use the largest resolution max(dx) in the set of images)
 npY=1+round((Rangy(1)-Rangy(2))/max(dy));
 
-
 x=linspace(Rangx(1),Rangx(2),npX);
 y=linspace(Rangy(1),Rangy(2),npY);
-[X,Y]=meshgrid(x,y);%grid in physical coordiantes
-%vec_B=[];
+[X,Y]=meshgrid(x,y);%grid in physical coordinates
 A_out=cell(1,numel(A));
-for icell=1:length(A) 
+
+for icell=1:length(A)
     Calib=XmlData{icell}.GeometryCalib;
     % rescaling of the image coordinates without change of the image array
@@ -60 +62 @@
         Rangy=[npy-0.5 0.5];
         [Rangx]=phys_XYZ(Calib,Rangx,[0.5 0.5],ZIndex);%case of translations without rotation and quadratic deformation
-        [xx,Rangy]=phys_XYZ(Calib,[0.5 0.5],Rangy,ZIndex);
+        [~,Rangy]=phys_XYZ(Calib,[0.5 0.5],Rangy,ZIndex);
     else
         % the image needs to be interpolated to the new coordinates
-        zphys=0; %default
+        Z=0; %default
         if isfield(Calib,'SliceCoord') %.Z= index of plane
             SliceCoord=Calib.SliceCoord(ZIndex,:);
-            zphys=SliceCoord(3); %to generalize for non-parallel planes
-            if isfield(Calib, 'SliceAngle') && ~isequal(Calib.SliceAngle,[0 0 0]) && ~isequal(Calib.SliceAngle(ZIndex,:),[0 0 0])
-                testangle=1;
-                om=norm(Calib.SliceAngle(ZIndex,:));%norm of rotation angle in radians
-                OmAxis=Calib.SliceAngle(ZIndex,:)/om; %unit vector marking the rotation axis
-                cos_om=cos(pi*om/180);
-                sin_om=sin(pi*om/180);
-                coeff=OmAxis(3)*(1-cos_om);
-                norm_plane(1)=OmAxis(1)*coeff+OmAxis(2)*sin_om;
-                norm_plane(2)=OmAxis(2)*coeff-OmAxis(1)*sin_om;
-                norm_plane(3)=OmAxis(3)*coeff+cos_om;
-                %Z0=norm_plane*Calib.SliceCoord(ZIndex,:)'/norm_plane(3);
-                Z0=Calib.SliceCoord(ZIndex,3);
-                zphys=Z0-(norm_plane(1)*X-norm_plane(2)*Y)/norm_plane(3);
+            Z=SliceCoord(3);
+            if isfield(Calib, 'SliceAngle') && size(Calib.SliceAngle,1)>=ZIndex && ~isequal(Calib.SliceAngle(ZIndex,:),[0 0 0])
+                norm_plane=angle2normal(Calib.SliceAngle(ZIndex,:));
+                Z=Z-(norm_plane(1)*(X-SliceCoord(1))+norm_plane(2)*(Y-SliceCoord(2)))/norm_plane(3);
             end
-%             if isfield(Calib,'InterfaceCoord') && isfield(Calib,'RefractionIndex')
-%                 H=Calib.InterfaceCoord(3);
-%                 if H>zphys
-%                     zphys=H-(H-zphys)/Calib.RefractionIndex; %corrected z (virtual object)
-%                 end
-%             end
         end
         xima=0.5:npx(icell)-0.5;%image coordinates of corners
         yima=npy(icell)-0.5:-1:0.5;
         [XIMA_init,YIMA_init]=meshgrid(xima,yima);%grid of initial image in px coordinates
-        [XIMA,YIMA]=px_XYZ(XmlData{icell}.GeometryCalib,X,Y,zphys);% image coordinates for each point in the real
+        [XIMA,YIMA]=px_XYZ(XmlData{icell}.GeometryCalib,X,Y,Z);% image coordinates for each point in the real
         testuint8=isa(A{icell},'uint8');
         testuint16=isa(A{icell},'uint16');
-        if ndims(A{icell})==2 %(B/W images)
-        A_out{icell}=interp2(XIMA_init,YIMA_init,double(A{icell}),XIMA,YIMA);
-         elseif ndims(A{icell})==3     
-             for icolor=1:size(A{icell},3)
-                 A{icell}=double(A{icell});
-                 A_out{icell}(:,:,icolor)=interp2(XIMA_init,YIMA_init,A{icell}(:,:,icolor),XIMA,YIMA);
-             end
-         end
+        if ismatrix(A{icell}) %(B/W images)
+            A_out{icell}=interp2(XIMA_init,YIMA_init,double(A{icell}),XIMA,YIMA);
+        elseif ndims(A{icell})==3
+            for icolor=1:size(A{icell},3)
+                A{icell}=double(A{icell});
+                A_out{icell}(:,:,icolor)=interp2(XIMA_init,YIMA_init,A{icell}(:,:,icolor),XIMA,YIMA);
+            end
+        end
         if testuint8
             A_out{icell}=uint8(A_out{icell});
@@ -107 +93 @@
         if testuint16
             A_out{icell}=uint16(A_out{icell});
-        end      
+        end
     end
 end

trunk/src/proj_field.m

@@ -985 +985 @@
 
 %-----------------------------------------------------------------
-%project on a plane
-% AJOUTER flux,circul,error
+% proj_plane: project on a plane defined by the structure ObjectData containing:
+%    .Type : = 'plane'
+%    .ProjMode (mode of projection) = 'projection'|'interp_lin'|'interp_tps' ;
+%    .CoordUnit: (for instance 'px','cm') units for the coordinates defining the plane  (the program checks that it fits with the unit of the input Field)
+%    .Angle : angles of rotation of the plane expressed in degrees. The first element
+%         ObjectData.Angle(1) represents a rotation in the plane (x,y) (around the
+%         vertical axis), which can be followed by a rotation with angle ObjectData.Angle(2) around the new (rotated) x axis.
+%    .Coord(1,3): coordinates (x,y,z) of the origin of the new coordinates in the projection plane;
+%    .DX,.DY,.DZ : increments along each coordinate for the projected data (for 'interp_lin' and 'interp_tps')
+%    .RangeX,RangeY: vectors with two elements defining the lower and upper bounds of the respectively X and Y coordinates in the projection plane
+%    .RangeInterp: maximum distance of interpolation from the known data. Interpolation yields NaN beyond this distance.
+
+% TODO: AJOUTER flux,circul,error
 function  [ProjData,errormsg] = proj_plane(FieldData, ObjectData)
 %-----------------------------------------------------------------
 
-%% rotation angles
+%% rotation matrix
 PlaneAngle=[0 0];
 norm_plane=[0 0 1];
-%cos_om=1;
-%sin_om=0;
 test90x=0;%=1 for 90 degree rotation alround x axis
 test90y=0;%=1 for 90 degree rotation alround y axis
-% if strcmp(ObjectData.Type,'plane_z')
-%     Delta_x=ObjectData.Coord(2,1)-ObjectData.Coord(1,1);
-%     Delta_y=ObjectData.Coord(2,2)-ObjectData.Coord(1,2);
-%     Delta_mod=sqrt(Delta_x*Delta_x+Delta_y*Delta_y);
-%     ObjectData.Angle=[0 0 0];
-%     ObjectData.Angle(1)=90*Delta_x/Delta_mod;
-%     ObjectData.0(2)=90*Delta_y/Delta_mod;
-% end
+
 if isfield(ObjectData,'Angle')
     checkM2=0;
@@ -1014 +1016 @@
     if PlaneAngle==0
         PlaneAngle=[0 0];
-    end
-    %     om=norm(PlaneAngle);%norm of rotation angle in radians
-    %     OmAxis=PlaneAngle/om; %unit vector marking the rotation axis
-    %     cos_om=cos(om);
-    %     sin_om=sin(om);
-    %     coeff=OmAxis(3)*(1-cos_om);
-    %     %components of the unity vector norm_plane normal to the projection plane
-    %     norm_plane(1)=OmAxis(1)*coeff+OmAxis(2)*sin_om;
-    %     norm_plane(2)=OmAxis(2)*coeff-OmAxis(1)*sin_om;
-    %     norm_plane(3)=OmAxis(3)*coeff+cos_om;
-    
+    end   
     M1=[cos(PlaneAngle(1)) -sin(PlaneAngle(1)) 0;sin(PlaneAngle(1)) cos(PlaneAngle(1)) 0;0 0 1];
     M=M1;
@@ -1031 +1023 @@
         M=M1*M2;% first rotate in the x,y plane with angle PlaneAngle(1), then slant around the new x axis0 with angle PlaneAngle(2)
     end
-    norm_plane=M*[0 0 1]';
-    
+    norm_plane=M*[0 0 1]';  
 end
 testangle=~isequal(PlaneAngle,[0 0])||~isequal(ObjectData.Coord(1:2),[0 0 ]) ;% && ~test90y && ~test90x;%=1 for slanted plane
@@ -1425 +1416 @@
                 %rotate coordinates if needed: coord_X,coord_Y= = coordinates in the new plane
                 Phi=PlaneAngle(1);
-                if testangle && ~test90y && ~test90x;%=1 for slanted plane
+                if testangle && ~test90y && ~test90x %=1 for slanted plane
                     new_XI=XI *cos(Phi) - YI* sin(Phi)+ObjectData.Coord(1);
                     YI=XI *sin(Phi) + YI *cos(Phi)+ObjectData.Coord(2);

trunk/src/series.m

@@ -1459 +1459 @@
     return
 end
-if isfield(Param,'InputFields')&& isfield(Param.InputFields,'FieldName')&& isequal(Param.InputFields.FieldName,'get_field...')
-    errormsg='input field name(s) not defined, select get_field...';
+if isfield(Param,'InputFields')&& isfield(Param.InputFields,'FieldName')&& isequal(Param.InputFields.FieldName,'add_field...')
+    errormsg='input field name(s) not defined, select add_field...';
     return
 end
@@ -2450 +2450 @@
         set(handles.VelType_1,'Visible','on')
         set(handles.VelType_title_1,'Visible','on')
-        FieldList_1=[set_field_list('U','V');{'C'};{'get_field...'}]; % standard menu for civx data
+        FieldList_1=[set_field_list('U','V');{'C'};{'add_field...'}]; % standard menu for civx data
         CheckList_1=1;
         set(handles.FieldName_1,'Value',1); % velocity vector choice by default
@@ -2467 +2467 @@
     if FieldNameRequest && isfield(SeriesData.FileInfo{iview_netcdf(1)},'ListVarName')
         set(handles.FieldName,'Visible','on')
+        set(handles.Field_text,'Visible','on')
         ListVarName=SeriesData.FileInfo{iview_netcdf(1)}.ListVarName;
         ind_var=get(handles.FieldName,'Value'); % indices of previously selected variables
@@ -2491 +2492 @@
     else
         set(handles.FieldName,'Visible','off')
+        set(handles.Field_text,'Visible','off')
     end
     
     set(handles_coord,'Visible','on')
-    if isempty(find(strcmp('get_field...',FieldList)))
-    FieldList=[FieldList;{'get_field...'}];%add 'get_field...' to the menu FieldName if it is not already
+    if isempty(find(strcmp('add_field...',FieldList)))
+    FieldList=[FieldList;{'add_field...'}];%add 'add_field...' to the menu FieldName if it is not already
     end
     if FieldNameRequest_1 && numel(iview_netcdf)>=2
@@ -2530 +2532 @@
     end
     if isempty(FieldList)
+        set(handles.Field_text,'Visible','off')
         set(handles.FieldName,'Visible','off')
     else
+        set(handles.Field_text,'Visible','on')
         set(handles.FieldName,'Visible','on')
         set(handles.FieldName,'String',FieldList)
@@ -2757 +2761 @@
 set(handles.ActionInput,'BackgroundColor',[1 0 0])
 
+
+%------------------------------------------------------------------------
+% --- Executes on button press in RefreshField.
+function RefreshField_Callback(hObject, eventdata, handles)
+%------------------------------------------------------------------------
+set(handles.FieldName,'String',{'add_field...'});
+set(handles.FieldName,'Value',1);
+FieldName_Callback(hObject, eventdata, handles)
+
+
 %------------------------------------------------------------------------
 % --- Executes on selection change in FieldName.
 function FieldName_Callback(hObject, eventdata, handles)
 %------------------------------------------------------------------------
-field_str=get(handles.FieldName,'String');
+FieldListInit=get(handles.FieldName,'String');
 field_index=get(handles.FieldName,'Value');
-field=field_str{field_index(1)};
-if isequal(field,'get_field...')
+field=FieldListInit{field_index(1)};
+if isequal(field,'add_field...')
+    FieldListInit(field_index(1))=[];
     SeriesData=get(handles.series,'UserData');
     % input line for which the field choice is relevant
@@ -2776 +2791 @@
     % check the existence of the first file in the series
     first_j=[];last_j=[];MinIndex_j=1;MaxIndex_j=1; % default setting for index j
-    if isfield(Param.IndexRange,'first_j'); % if index j is used     
+    if isfield(Param.IndexRange,'first_j') % if index j is used     
         first_j=Param.IndexRange.first_j;
         last_j=Param.IndexRange.last_j;
@@ -2789 +2804 @@
         answer=msgbox_uvmat('INPUT_TXT',['select the line of the input table:' num2str(iview)] ,num2str(iview(1)));
         LineIndex=str2num(answer);
-%         InputLine=str2num(get(handles.InputLine,'String'));
-%         if ismember(InputLine,iview)
-%             LineIndex=InputLine;
-%         end
     end
     FirstFileName=fullfile_uvmat(InputTable{LineIndex,1},InputTable{LineIndex,2},InputTable{LineIndex,3},...
@@ -2826 +2837 @@
         end
         set(handles.FieldName,'Value',1)
-        set(handles.FieldName,'String',[FieldList; {'get_field...'}]);
+        set(handles.FieldName,'String',[FieldListInit; FieldList; {'add_field...'}]);
         if ~strcmp(GetFieldData.FieldOption,'civdata...')
            if ~isempty(regexp(FieldList{1},'^vec'))
@@ -2913 +2924 @@
 field_index=get(handles.FieldName_1,'Value');
 field=field_str{field_index(1)};
-if isequal(field,'get_field...')
+if isequal(field,'add_field...')
     hget_field=findobj(allchild(0),'name','get_field');
     if ~isempty(hget_field)
@@ -2973 +2984 @@
         end
         set(handles.FieldName_1,'Value',1)
-        set(handles.FieldName_1,'String',[FieldList; {'get_field...'}]);
+        set(handles.FieldName_1,'String',[FieldList; {'add_field...'}]);
     end
 end   
@@ -3378 +3389 @@
     ListField=Param.InputFields.FieldName;
     if ischar(ListField),ListField={ListField}; end
-    set(handles.FieldName,'String',[ListField;{'get-field...'}])
+    set(handles.FieldName,'String',[ListField;{'add_field...'}])
      set(handles.FieldName,'Value',1:numel(ListField))
      set(handles.FieldName,'Visible','on')
@@ -3866 +3877 @@
 
 
-
 % --- Executes on button press in Replicate.
 function Replicate_Callback(hObject, eventdata, handles)
@@ -3882 +3892 @@
 end
 
-
-

trunk/src/series/merge_proj.m

@@ -320 +320 @@
         check_proj_tps= strcmp(FileType{iview},'civdata') && isfield(Param,'ProjObject')&&~isempty(Param.ProjObject)...
             && strcmp(Param.ProjObject.ProjMode,'interp_tps')&&~isfield(Data{iview},'Coord_tps');
+        if check_proj_tps
         Data{iview}=tps_coeff_field(Data{iview},check_proj_tps);
+        end
         
         %% projection on object (gridded plane)

trunk/src/transform_field/phys_polar.m

@@ -3 +3 @@
 %%%%  Use the general syntax for transform fields %%%%
 % OUTPUT: 
-% DataOut:   output field structure 
+% Data:   output field structure 
 %      .X=radius, .Y=azimuth angle, .U, .V are radial and azimuthal velocity components
 %
@@ -14 +14 @@
 %         .GeometryCalib: substructure of the calibration parameters 
 % transform image coordinates (px) to polar physical coordinates 
-%[DataOut,DataOut_1]=phys_polar(varargin)
+%[Data,Data_1]=phys_polar(varargin)
 %
 % OUTPUT: 
-% DataOut: structure of modified data field: .X=radius, .Y=azimuth angle, .U, .V are radial and azimuthal velocity components
-% DataOut_1:  second data field (if two fields are in input)
+% Data: structure of modified data field: .X=radius, .Y=azimuth angle, .U, .V are radial and azimuthal velocity components
+% Data_1:  second data field (if two fields are in input)
 %
 %INPUT:
@@ -44 +44 @@
 %=======================================================================
 
-function DataOut=phys_polar(DataIn,XmlData,DataIn_1,XmlData_1)
+function Data=phys_polar(DataIn,XmlData,DataIn_1,XmlData_1)
 %------------------------------------------------------------------------
+
 %% request input parameters
 if isfield(DataIn,'Action') && isfield(DataIn.Action,'RUN') && isequal(DataIn.Action.RUN,0)
@@ -64 +65 @@
     end
     answer = inputdlg(prompt,dlg_title,num_lines,def);
-    DataOut.TransformInput.PolarCentre=str2num(answer{1}); 
-    DataOut.TransformInput.PolarReferenceRadius=str2num(answer{2}); 
-    DataOut.TransformInput.PolarReferenceAngle=str2num(answer{3}); 
-    if isfield(XmlData,'GeometryCalib')&& isfield(XmlData.GeometryCalib,'CoordUnit')
-        DataOut.CoordUnit=XmlData.GeometryCalib.CoordUnit;% states that the output is in unit defined by GeometryCalib, then erased all projection objects with different units
-    end
+    Data.TransformInput.PolarCentre=str2num(answer{1}); 
+    Data.TransformInput.PolarReferenceRadius=str2num(answer{2}); 
+    Data.TransformInput.PolarReferenceAngle=str2num(answer{3}); 
+%     if isfield(XmlData,'GeometryCalib')&& isfield(XmlData.GeometryCalib,'CoordUnit')
+%         Data.CoordUnit=XmlData.GeometryCalib.CoordUnit;% states that the output is in unit defined by GeometryCalib, then erased all projection objects with different units
+%     end
     return
 end
 
+%% default outputs
+Data=DataIn; %default output
+if isfield(Data,'CoordUnit')
+Data=rmfield(Data,'CoordUnit');
+end
+Data.ListVarName = {};
+Data.VarDimName={};
+Data.VarAttribute={};
+DataCell{1}=DataIn;
 Calib{1}=[];
+DataCell{2}=[];%default
+checkpixel(1)=0;
+if isfield(DataCell{1},'CoorUnit')&& strcmp(DataCell{1}.CoorUnit,'px') 
+    checkpixel(1)=1;
+end
 if nargin==2||nargin==4
-    DataOut=DataIn;%default
-    DataOut_1=[];%default
-    if isfield(XmlData,'GeometryCalib')
+    if isfield(XmlData,'GeometryCalib') && ~isempty(XmlData.GeometryCalib)&& checkpixel(1)
         Calib{1}=XmlData.GeometryCalib;
     end
     Calib{2}=Calib{1};
 else
-    DataOut.Txt='wrong input: need two or four structures';
-end
-test_1=0;
+    Data.Txt='wrong input: need two or four structures';
+end
+nbinput=1;
 if nargin==4% case of two input fields
-    test_1=1;
-    DataOut_1=DataIn_1;%default
-    if isfield(XmlData_1,'GeometryCalib')
+    checkpixel(2)=0;
+if isfield(DataCell{2},'CoorUnit')&& strcmp(DataCell{2}.CoorUnit,'px') 
+    checkpixel(2)=1;
+end
+    DataCell{2}=DataIn_1;%default
+    if isfield(XmlData_1,'GeometryCalib')&& ~isempty(XmlData_1.GeometryCalib) && checkpixel(2)
         Calib{2}=XmlData_1.GeometryCalib;
     end
-end
-
-%parameters for polar coordinates (taken from the calibration data of the first field)
+    nbinput=2;
+end
+
+%% parameters for polar coordinates (taken from the calibration data of the first field)
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-XmlData.PolarReferenceRadius=0;%450;
-XmlData.PolarReferenceAngle=0;%450*pi/2;
 origin_xy=[0 0];%center for the polar coordinates in the original x,y coordinates
 radius_offset=0;%reference radius used to offset the radial coordinate r
 angle_offset=0; %reference angle used as new origin of the polar angle (= axis Ox by default)
-angle_scale=180/pi; 
+angle_scale=180/pi;
+check_degree=1;%angle expressed in degrees by default
 if isfield(XmlData,'TransformInput')
     if isfield(XmlData.TransformInput,'PolarCentre') && isnumeric(XmlData.TransformInput.PolarCentre)
-        if isequal(length(XmlData.TransformInput.PolarCentre),2);
+        if isequal(length(XmlData.TransformInput.PolarCentre),2)
             origin_xy= XmlData.TransformInput.PolarCentre;
         end
@@ -112 +128 @@
     if radius_offset > 0
         angle_scale=radius_offset; %the azimuth is rescale in terms of the length along the reference radius
+        check_degree=0; %the output has the same unit asthe input
     else
         angle_scale=180/pi; %polar angle in degrees
+        check_degree=1;%angle expressed in degrees
     end
     if isfield(XmlData.TransformInput,'PolarReferenceAngle') && isnumeric(XmlData.TransformInput.PolarReferenceAngle)
@@ -121 +139 @@
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-iscalar=0;
-%transform first field to cartesian phys coordiantes
-if  ~isempty(Calib{1})
-    DataOut=phys_1(DataIn,Calib{1},origin_xy,radius_offset,angle_offset,angle_scale);
-    %case of images or scalar
-    if isfield(DataIn,'A')&isfield(DataIn,'Coord_x')&~isempty(DataIn.Coord_x) & isfield(DataIn,'Coord_y')&...
-                                           ~isempty(DataIn.Coord_y)&length(DataIn.A)>1
-        iscalar=1;
-        A{1}=DataIn.A;
-    end
+%% get fields
+check_scalar=0;
+check_vector=0;
+nbvar=0;%counter for the number of output variables
+nbcoord=0;%counter for the number of variables for radial coordiantes (case of multiple field inputs)
+nbgrid=0;%counter for the number of gridded fields (all linearly interpolated on the same output polar grid)
+nbscattered=0;%counter of scattered fields
+radius_name='radius';
+theta_name='theta';
+U_r_name='U_r';
+U_theta_name='U_theta';
+for ifield=1:nbinput %1 or 2 input fields
+    [CellInfo,NbDim,errormsg]=find_field_cells(DataCell{ifield});
+    if ~isempty(errormsg)
+        Data.Txt=['bad input to phys_polar: ' errormsg];
+        return
+    end
+    % end
     %transform of X,Y coordinates for vector fields
-    if isfield(DataIn,'ZIndex')&~isempty(DataIn.ZIndex)
-        ZIndex=DataIn.ZIndex;
+    if isfield(DataCell{ifield},'ZIndex')&& ~isempty(DataCell{ifield}.ZIndex)
+        ZIndex=DataCell{ifield}.ZIndex;
     else
         ZIndex=0;
     end
-end
-
-%transform second field (if exists) to cartesian phys coordiantes
-if test_1
-    DataOut_1=phys_1(Data_1,Calib{2},origin_xy,radius_offset,angle_offset,angle_scale);
-    if isfield(Data_1,'A')&isfield(Data_1,'Coord_x')&~isempty(Data_1.Coord_x) & isfield(Data_1,'Coord_y')&...
-                                       ~isempty(Data_1.Coord_y)&length(Data_1.A)>1
-          iscalar=iscalar+1;
-          Calib{iscalar}=Calib{2};
-          A{iscalar}=Data_1.A;
-          if isfield(Data_1,'ZIndex')&~isequal(Data_1.ZIndex,ZIndex)
-              DataOut.Txt='inconsistent plane indexes in the two input fields';
-          end
-          if iscalar==1% case for which only the second field is a scalar
-               [A,Coord_x,Coord_y]=phys_Ima_polar(A,Calib,ZIndex,origin_xy,radius_offset,angle_offset,angle_scale);
-               DataOut_1.A=A{1};
-               DataOut_1.Coord_x=Coord_x; 
-               DataOut_1.Coord_y=Coord_y;
-               return
-          end
-    end
-end
-if iscalar~=0
-    [A,Coord_x,Coord_y]=phys_Ima_polar(A,Calib,ZIndex,origin_xy,radius_offset,angle_offset,angle_scale);%
-    DataOut.A=A{1};
-    DataOut.Coord_x=Coord_x; 
-    DataOut.Coord_y=Coord_y;
-    if iscalar==2
-        DataOut_1.A=A{2};
-        DataOut_1.Coord_x=Coord_x; 
-        DataOut_1.Coord_y=Coord_y;
-    end
-end
-
-
+    for icell=1:numel(CellInfo)
+        if NbDim(icell)==2
+            % case of input field with scattered coordinates
+            if strcmp(CellInfo{icell}.CoordType,'scattered')
+                nbscattered=nbscattered+1;
+                nbcoord=nbcoord+1;
+                radius_name = rename_indexing(radius_name,Data.ListVarName);
+                theta_name = rename_indexing(theta_name,Data.ListVarName);
+                Data.ListVarName = [Data.ListVarName {radius_name} {theta_name}];
+                dim_name = rename_indexing('nb_point',Data.VarDimName);
+                Data.VarDimName=[Data.VarDimName {dim_name} {dim_name}];
+                nbvar=nbvar+2;
+                Data.VarAttribute{nbvar-1}.Role='coord_x';
+                check_unit=1;
+                if isfield(DataCell{ifield},'CoordUnit')
+                    Data=rmfield(Data,'CoordUnit');
+                    Data.VarAttribute{nbvar-1}.unit=DataCell{ifield}.CoordUnit;
+                elseif isfield(XmlData,'GeometryCalib')&& isfield(XmlData.GeometryCalib,'CoordUnit')
+                    Data.VarAttribute{nbvar-1}.unit=XmlData.GeometryCalib.CoordUnit;% states that the output is in unit defined by GeometryCalib, then erased all projection objects with different units
+                else
+                    check_unit=0;
+                end
+                Data.VarAttribute{nbvar}.Role='coord_y';
+                if check_degree
+                Data.VarAttribute{nbvar}.unit='degree';
+                elseif check_unit
+                    Data.VarAttribute{nbvar}.unit=Data.VarAttribute{nbvar-1}.unit;
+                end
+  
+                %transform u,v into polar coordinates
+                X=DataCell{ifield}.(CellInfo{icell}.XName);
+                Y=DataCell{ifield}.(CellInfo{icell}.YName);
+                if isfield(CellInfo{icell},'VarIndex_vector_x')&& isfield(CellInfo{icell},'VarIndex_vector_y')
+                    UName=DataCell{ifield}.ListVarName{CellInfo{icell}.VarIndex_vector_x};
+                    VName=DataCell{ifield}.ListVarName{CellInfo{icell}.VarIndex_vector_y};
+                    if ~isempty(Calib{ifield})
+                        [X,Y,Z,DataCell{ifield}.(UName),DataCell{ifield}.(VName)]=...
+                            phys_XYUV(DataCell{ifield},Calib{ifield},ZIndex);
+                    end
+                end
+                [Theta,Radius] = cart2pol(X-origin_xy(1),Y-origin_xy(2));
+                Data.(radius_name)=Radius-radius_offset;
+                Data.(theta_name)=Theta*angle_scale-angle_offset;
+                if Z~=0
+                    Data.Z=Z;
+                    nbvar=nbvar+1;
+                    Data.ListVarName = [Data.ListVarName {'Z'}];
+                    Data.VarDimName=[Data.VarDimName {dim_name}];
+                    Data.VarAttribute{nbvar}.Role='coord_z';
+                end
+                if isfield(CellInfo{icell},'VarIndex_scalar')
+                    ScalarName=DataCell{ifield}.ListVarName{CellInfo{icell}.VarIndex_scalar};
+                    ScalarName=rename_indexing(ScalarName,Data.ListVarName);
+                    Data.(ScalarName)=DataCell{ifield}.(ScalarName);
+                    nbvar=nbvar+1;
+                    Data.ListVarName = [Data.ListVarName {ScalarName}];
+                    Data.VarDimName=[Data.VarDimName {dim_name}];
+                    Data.VarAttribute{nbvar}.Role='scalar';
+                end
+                if isfield(CellInfo{icell},'VarIndex_vector_x')&& isfield(CellInfo{icell},'VarIndex_vector_y')
+                    U_r_name= rename_indexing(U_r_name,Data.ListVarName);
+                    U_theta_name= rename_indexing(U_theta_name,Data.ListVarName);
+                    Data.(U_r_name)=DataCell{ifield}.(UName).*cos(Theta)+DataCell{ifield}.(VName).*sin(Theta);%radial velocity
+                    Data.(U_theta_name)=(-DataCell{ifield}.(UName).*sin(Theta)+DataCell{ifield}.(VName).*cos(Theta));%./(Data.X)%+radius_ref);% azimuthal velocity component
+                    Data.ListVarName = [Data.ListVarName {U_r_name} {U_theta_name}];
+                    Data.VarDimName=[Data.VarDimName {dim_name} {dim_name}];
+                    Data.VarAttribute{nbvar+1}.Role='vector_x';
+                    Data.VarAttribute{nbvar+2}.Role='vector_y';
+                    nbvar=nbvar+2;
+                end
+                if isfield(CellInfo{icell},'VarIndex_errorflag')
+                    error_flag_name=DataCell{ifield}.ListVarName{CellInfo{icell}.VarIndex_errorflag};
+                    error_flag_newname= rename_indexing(error_flag_name,Data.ListVarName);
+                    Data.(error_flag_newname)=DataCell{ifield}.(error_flag_name);
+                    Data.ListVarName = [Data.ListVarName {error_flag_newname}];
+                    Data.VarDimName=[Data.VarDimName {dim_name}];
+                    nbvar=nbvar+1;
+                    Data.VarAttribute{nbvar}.Role='errorflag';
+                end
+                
+           %caseof input fields on gridded coordinates (matrix)
+            elseif strcmp(CellInfo{icell}.CoordType,'grid')
+                if nbgrid==0% no gridded data yet, introduce the coordinate variables common to all gridded data
+                    nbcoord=nbcoord+1;%add new radial coordinates for the first gridded field
+                    radius_name = rename_indexing(radius_name,Data.ListVarName);
+                    theta_name = rename_indexing(theta_name,Data.ListVarName);
+                    Data.ListVarName = [Data.ListVarName {radius_name} {theta_name}];
+                    Data.VarDimName=[Data.VarDimName {radius_name} {theta_name}];
+                    nbvar=nbvar+2;
+                    Data.VarAttribute{nbvar-1}.Role='coord_x';
+                    Data.VarAttribute{nbvar}.Role='coord_y';
+                    check_unit=1;
+                    if isfield(DataCell{ifield},'CoordUnit')
+                        Data.VarAttribute{nbvar-1}.unit=DataCell{ifield}.CoordUnit;
+                    elseif isfield(XmlData,'GeometryCalib')&& isfield(XmlData.GeometryCalib,'CoordUnit')
+                        Data.VarAttribute{nbvar-1}.unit=XmlData.GeometryCalib.CoordUnit;% states that the output is in unit defined by GeometryCalib, then erased all projection objects with different units
+                    else
+                        check_unit=0;
+                    end
+                    if check_degree
+                        Data.VarAttribute{nbvar}.unit='degree';
+                    elseif check_unit
+                        Data.VarAttribute{nbvar}.unit=Data.VarAttribute{nbvar-1}.unit;
+                    end
+                end
+                if isfield(CellInfo{icell},'VarIndex_scalar')
+                    nbgrid=nbgrid+1;
+                    nbvar=nbvar+1;
+                    Data.VarAttribute{nbvar}.Role='scalar';
+                    FieldName{nbgrid}=DataCell{ifield}.ListVarName{CellInfo{icell}.VarIndex_scalar};
+                    A{nbgrid}=DataCell{ifield}.(FieldName{nbgrid});
+%                     Data.ListVarName=[Data.ListVarName {FieldName{nbgrid}}];
+%                     Data.VarDimName=[Data.VarDimName {{theta_name,radius_name}}];
+                    nbpoint(nbgrid)=numel(A{nbgrid});
+                    check_scalar(nbgrid)=1;
+                    coord_x{nbgrid}=DataCell{ifield}.(DataCell{ifield}.ListVarName{CellInfo{icell}.XIndex});
+                    coord_y{nbgrid}=DataCell{ifield}.(DataCell{ifield}.ListVarName{CellInfo{icell}.YIndex});
+                    ZInd(nbgrid)=ZIndex;
+                    Calib_new{nbgrid}=Calib{ifield};
+                end
+                if isfield(CellInfo{icell},'VarIndex_vector_x')&& isfield(CellInfo{icell},'VarIndex_vector_y')
+                    FieldName{nbgrid+1}=DataCell{ifield}.ListVarName{CellInfo{icell}.VarIndex_vector_x};
+                    FieldName{nbgrid+2}=DataCell{ifield}.ListVarName{CellInfo{icell}.VarIndex_vector_y};
+                    A{nbgrid+1}=DataCell{ifield}.(FieldName{nbgrid+1});
+                    A{nbgrid+2}=DataCell{ifield}.(FieldName{nbgrid+2});
+                   % Data.ListVarName=[Data.ListVarName {'U_r','U_theta'}];
+                    %Data.VarDimName=[Data.VarDimName {{theta_name,radius_name}} {{theta_name,radius_name}}];
+                    Data.VarAttribute{nbvar+1}.Role='vector_x';
+                    Data.VarAttribute{nbvar+2}.Role='vector_y';
+                    nbpoint([nbgrid+1 nbgrid+2])=numel(A{nbgrid+1});
+                    check_vector(nbgrid+1)=1;
+                    check_vector(nbgrid+2)=1;
+                    coord_x{nbgrid+1}=DataCell{ifield}.(DataCell{ifield}.ListVarName{CellInfo{icell}.XIndex});
+                    coord_y{nbgrid+1}=DataCell{ifield}.(DataCell{ifield}.ListVarName{CellInfo{icell}.YIndex});
+                    coord_x{nbgrid+2}=DataCell{ifield}.(DataCell{ifield}.ListVarName{CellInfo{icell}.XIndex});
+                    coord_y{nbgrid+2}=DataCell{ifield}.(DataCell{ifield}.ListVarName{CellInfo{icell}.YIndex});
+                    ZInd(nbgrid+1)=ZIndex;
+                    ZInd(nbgrid+2)=ZIndex;
+                    Calib_new{nbgrid+1}=Calib{ifield};
+                    Calib_new{nbgrid+2}=Calib{ifield};
+                    nbgrid=nbgrid+2;
+                    nbvar=nbvar+2;
+                end
+            end
+        end
+    end
+end
+
+%% tranform cartesian to polar coordinates for gridded data
+if nbgrid~=0
+    [A,Data.radius,Data.theta]=phys_Ima_polar(A,coord_x,coord_y,Calib_new,ZInd,origin_xy,radius_offset,angle_offset,angle_scale);
+    for icell=1:numel(A)
+        if icell<=numel(A)-1 && check_vector(icell)==1 && check_vector(icell+1)==1   %transform u,v into polar coordiantes
+            theta=Data.theta/angle_scale-angle_offset;
+            [~,Theta]=meshgrid(Data.radius,theta);%grid in physical coordinates
+            U_r_name= rename_indexing(U_r_name,Data.ListVarName);
+            U_theta_name= rename_indexing(U_theta_name,Data.ListVarName);
+            Data.ListVarName=[Data.ListVarName {U_r_name,U_theta_name}];
+            Data.VarDimName=[Data.VarDimName {{theta_name,radius_name}} {{theta_name,radius_name}}];
+            Data.(U_r_name)=A{icell}.*cos(Theta)+A{icell+1}.*sin(Theta);%radial velocity
+            Data.(U_theta_name)=(-A{icell}.*sin(Theta)+A{icell+1}.*cos(Theta));%./(Data.X)%+radius_ref);% azimuthal velocity component
+        elseif ~check_vector(icell)% for scalar fields
+            FieldName{icell}= rename_indexing(FieldName{icell},Data.ListVarName);
+            Data.ListVarName=[Data.ListVarName {FieldName{icell}}];
+            Data.VarDimName=[Data.VarDimName {{theta_name,radius_name}}];
+            Data.(FieldName{icell})=A{icell};
+        end
+    end
+end
 
 
 %------------------------------------------------
-function DataOut=phys_1(Data,Calib,origin_xy,radius_offset,angle_offset,angle_scale)
-
-DataOut=Data;
-% DataOut.CoordUnit=Calib.CoordUnit; %put flag for physical coordinates
-if isfield(Calib,'SliceCoord')
-    DataOut.PlaneCoord=Calib.SliceCoord;%to generalise for any plane 
-end
-
-if isfield(Data,'CoordUnit')%&& isequal(Data.CoordType,'px')&& ~isempty(Calib)
-    if isfield(Calib,'CoordUnit')
-        DataOut.CoordUnit=Calib.CoordUnit;
-    else
-        DataOut.CoordUnit='cm'; %default
-    end
-    DataOut.TimeUnit='s';
-    %transform of X,Y coordinates for vector fields
-    if isfield(Data,'ZIndex') && ~isempty(Data.ZIndex)&&~isnan(Data.ZIndex)
-        Z=Data.ZIndex;
-    else
-        Z=0;
-    end
-    if isfield(Data,'X') &isfield(Data,'Y')&~isempty(Data.X) & ~isempty(Data.Y)
-        [DataOut.X,DataOut.Y,DataOut.Z]=phys_XYZ(Calib,Data.X,Data.Y,Z); %transform from pixels to physical
-        DataOut.X=DataOut.X-origin_xy(1);%origin of coordinates at the tank center
-        DataOut.Y=DataOut.Y-origin_xy(2);%origin of coordinates at the tank center
-        [theta,DataOut.X] = cart2pol(DataOut.X,DataOut.Y);%theta  and X are the polar coordinates angle and radius
-          %shift and renormalize the polar coordinates
-        DataOut.X=DataOut.X-radius_offset;%shift the origin of radius, taken as the new X coordinate
-        DataOut.Y=(theta-angle_offset)*angle_scale;% normalized angle: distance along reference radius,taken as the new Y coordinate
-        %transform velocity field if exists
-        if isfield(Data,'U') & isfield(Data,'V') & ~isempty(Data.U) & ~isempty(Data.V)& isfield(Data,'Dt') 
-            if ~isempty(Data.Dt)
-            [XOut_1,YOut_1]=phys_XYZ(Calib,Data.X-Data.U/2,Data.Y-Data.V/2,Z);% X,Y positions of the vector origin in phys
-            [XOut_2,YOut_2]=phys_XYZ(Calib,Data.X+Data.U/2,Data.Y+Data.V/2,Z);% X,Y positions of the vector end in phys
-            UX=(XOut_2-XOut_1)/Data.Dt;% phys velocity u component
-            VY=(YOut_2-YOut_1)/Data.Dt; % phys velocity v component     
-            %transform u,v into polar coordiantes
-            DataOut.U=UX.*cos(theta)+VY.*sin(theta);%radial velocity
-            DataOut.V=(-UX.*sin(theta)+VY.*cos(theta));%./(DataOut.X)%+radius_ref);% azimuthal velocity component 
-            %shift and renormalize the angular velocity
-            end
-        end
-        %transform of spatial derivatives
-        if isfield(Data,'X') && ~isempty(Data.X) && isfield(Data,'DjUi') && ~isempty(Data.DjUi)...
-                && isfield(Data,'Dt')
-            if ~isempty(Data.Dt)
-                % estimate the Jacobian matrix DXpx/DXphys
-                for ip=1:length(Data.X)
-                    [Xp1,Yp1]=phys_XYZ(Calib,Data.X(ip)+0.5,Data.Y(ip),Z);
-                    [Xm1,Ym1]=phys_XYZ(Calib,Data.X(ip)-0.5,Data.Y(ip),Z);
-                    [Xp2,Yp2]=phys_XYZ(Calib,Data.X(ip),Data.Y(ip)+0.5,Z);
-                    [Xm2,Ym2]=phys_XYZ(Calib,Data.X(ip),Data.Y(ip)-0.5,Z);
-                    %Jacobian matrix DXpphys/DXpx
-                    DjXi(1,1)=(Xp1-Xm1);
-                    DjXi(2,1)=(Yp1-Ym1);
-                    DjXi(1,2)=(Xp2-Xm2);
-                    DjXi(2,2)=(Yp2-Ym2);
-                    DjUi(:,:)=Data.DjUi(ip,:,:);
-                    DjUi=(DjXi*DjUi')/DjXi;% =J-1*M*J , curvature effects (derivatives of J) neglected
-                    DataOut.DjUi(ip,:,:)=DjUi';
-                end
-                DataOut.DjUi =  DataOut.DjUi/Data.Dt;   %     min(Data.DjUi(:,1,1))=DUDX
-            end
-        end
-    end
-end
-
+%--- transform a single field into phys coordiantes
+function [X,Y,Z,U,V]=phys_XYUV(Data,Calib,ZIndex)
+%------------------------------------------------
+%% set default output
+%DataOut=Data;%default
+%DataOut.CoordUnit=Calib.CoordUnit;% the output coord unit is set by the calibration parameters
+X=[];%default output
+Y=[];
+Z=0;
+U=[];
+V=[];
+%% transform  X,Y coordinates for velocity fields (transform of an image or scalar done in phys_ima)
+if isfield(Data,'X') &&isfield(Data,'Y')&&~isempty(Data.X) && ~isempty(Data.Y)
+    [X,Y,Z]=phys_XYZ(Calib,Data.X,Data.Y,ZIndex);
+    Dt=1; %default
+    if isfield(Data,'dt')&&~isempty(Data.dt)
+        Dt=Data.dt;
+    end
+    if isfield(Data,'Dt')&&~isempty(Data.Dt)
+        Dt=Data.Dt;
+    end
+    if isfield(Data,'U')&&isfield(Data,'V')&&~isempty(Data.U) && ~isempty(Data.V)
+        [XOut_1,YOut_1]=phys_XYZ(Calib,Data.X-Data.U/2,Data.Y-Data.V/2,ZIndex);
+        [XOut_2,YOut_2]=phys_XYZ(Calib,Data.X+Data.U/2,Data.Y+Data.V/2,ZIndex);
+        U=(XOut_2-XOut_1)/Dt;
+        V=(YOut_2-YOut_1)/Dt;
+    end
+end
 
 %%%%%%%%%%%%%%%%%%%%
-function [A_out,Rangx,Rangy]=phys_Ima_polar(A,CalibIn,ZIndex,origin_xy,radius_offset,angle_offset,angle_scale)
-xcorner=[];
-ycorner=[];
+% tranform gridded field into polar coordiantes on a regular polar grid,
+% transform to phys coordiantes if requested by calibration input
+function [A_out,radius,theta]=phys_Ima_polar(A,coord_x,coord_y,CalibIn,ZIndex,origin_xy,radius_offset,angle_offset,angle_scale)
+rcorner=[];
+thetacorner=[];
 npx=[];
 npy=[];
-NbPoints=20; % nbre of points used to determine the image edge
 for icell=1:length(A)
     siz=size(A{icell});
-    npx=[npx siz(2)];
-    npy=[npy siz(1)];
-    zphys=0; %default
-    if isfield(CalibIn{icell},'SliceCoord') %.Z= index of plane
-        if ZIndex==0
-            ZIndex=1;
-        end
-       SliceCoord=CalibIn{icell}.SliceCoord(ZIndex,:);
-       zphys=SliceCoord(3); %to generalize for non-parallel planes
-    end
-%     xima=[0.5 siz(2)-0.5 0.5 siz(2)-0.5];%image coordinates of corners
-%     yima=[0.5 0.5 siz(1)-0.5 siz(1)-0.5];
-      edge_x=linspace(0.5,siz(1)-0.5,NbPoints);      
-      edge_y=linspace(0.5,siz(2)-0.5,NbPoints);
-      xima=[edge_y (siz(2)-0.5)*ones(1,NbPoints) edge_y 0.5*ones(1,NbPoints)];%image coordinates of corners
-      yima=[0.5*ones(1,NbPoints) edge_x (siz(1)-0.5)*ones(1,NbPoints) edge_x];%image coordinates of corners
-    [xcorner_new,ycorner_new]=phys_XYZ(CalibIn{icell},xima,yima,ZIndex);%corresponding physical coordinates
-    %transform the corner coordinates into polar ones    
-    xcorner_new=xcorner_new-origin_xy(1);%shift to the origin of the polar coordinates 
-    ycorner_new=ycorner_new-origin_xy(2);%shift to the origin of the polar coordinates       
-    [theta,xcorner_new] = cart2pol(xcorner_new,ycorner_new);%theta  and X are the polar coordinates angle and radius
-    if (max(theta)-min(theta))>pi   %if the polar origin is inside the image
-        xcorner_new=[0 max(xcorner_new)];
-        theta=[-pi pi];
-    end
-          %shift and renormalize the polar coordinates
-    xcorner_new=xcorner_new-radius_offset;%
-    ycorner_new=theta*angle_scale-angle_offset;% normalized angle: distance along reference radius
-    xcorner=[xcorner xcorner_new];
-    ycorner=[ycorner ycorner_new];
-end
-Rangx(1)=min(xcorner);
-Rangx(2)=max(xcorner);
-Rangy(2)=min(ycorner);
-Rangy(1)=max(ycorner);
-% test_multi=(max(npx)~=min(npx)) | (max(npy)~=min(npy)); 
-npx=max(npx);
-npy=max(npy);
-x=linspace(Rangx(1),Rangx(2),npx);
-y=linspace(Rangy(1),Rangy(2),npy);
-[X,Y]=meshgrid(x,y);%grid in physical coordinates
+    npx(icell)=siz(2);
+    npy(icell)=siz(1);
+    x_edge=[linspace(coord_x{icell}(1),coord_x{icell}(end),npx(icell)) coord_x{icell}(end)*ones(1,npy(icell))...
+        linspace(coord_x{icell}(end),coord_x{icell}(1),npx(icell)) coord_x{icell}(1)*ones(1,npy(icell))];%x coordinates of the image edge(four sides)
+    y_edge=[coord_y{icell}(1)*ones(1,npx(icell)) linspace(coord_y{icell}(1),coord_y{icell}(end),npy(icell))...
+        coord_y{icell}(end)*ones(1,npx(icell)) linspace(coord_y{icell}(end),coord_y{icell}(1),npy(icell))];%y coordinates of the image edge(four sides)
+    
+    % transform edges into phys coordinates if requested
+    if ~isempty(CalibIn{icell})
+        [x_edge,y_edge]=phys_XYZ(CalibIn{icell},x_edge,y_edge,ZIndex(icell));% physical coordinates of the image edge
+    end
+    
+    %transform the corner coordinates into polar ones
+    x_edge=x_edge-origin_xy(1);%shift to the origin of the polar coordinates
+    y_edge=y_edge-origin_xy(2);%shift to the origin of the polar coordinates
+    [theta_edge,r_edge] = cart2pol(x_edge,y_edge);%theta  and X are the polar coordinates angle and radius
+    if (max(theta_edge)-min(theta_edge))>pi   %if the polar origin is inside the image
+        r_edge=[0 max(r_edge)];
+        theta_edge=[-pi pi];
+    end
+    rcorner=[rcorner r_edge];
+    thetacorner=[thetacorner theta_edge];
+end
+nbpoint=max(npx.*npy);
+Min_r=min(rcorner);
+Max_r=max(rcorner);
+Min_theta=min(thetacorner)*angle_scale;
+Max_theta=max(thetacorner)*angle_scale;
+Dr=round_uvmat((Max_r-Min_r)/sqrt(nbpoint));
+Dtheta=round_uvmat((Max_theta-Min_theta)/sqrt(nbpoint));% get a simple mesh for the rescaled angle
+radius=Min_r:Dr:Max_r;% polar coordinates for projections
+theta=Min_theta:Dtheta:Max_theta;
+[Radius,Theta]=meshgrid(radius,theta/angle_scale);%grid in polar coordinates (angles in radians)
 %transform X, Y in cartesian
-X=X+radius_offset;%
-Y=(Y+angle_offset)/angle_scale;% normalized angle: distance along reference radius
-[X,Y] = pol2cart(Y,X);
-X=X+origin_xy(1);%shift to the origin of the polar coordinates 
-Y=Y+origin_xy(2);%shift to the origin of the polar coordinates 
-for icell=1:length(A) 
+[X,Y] = pol2cart(Theta,Radius);% cartesian coordinates associated to the grid in polar coordinates
+X=X+origin_xy(1);%shift to the origin of the polar coordinates
+Y=Y+origin_xy(2);%shift to the origin of the polar coordinates
+radius=radius-radius_offset;
+theta=theta-angle_offset*angle_scale;
+[np_theta,np_r]=size(Radius);
+
+for icell=1:length(A)
+    XIMA=X;
+    YIMA=Y;
+    if ~isempty(CalibIn{icell})%transform back to pixel if calibration parameters are introduced
+        Z=0; %default
+        if isfield(CalibIn{icell},'SliceCoord') %.Z= index of plane
+            if ZIndex(icell)==0
+                ZIndex(icell)=1;
+            end
+            SliceCoord=CalibIn{icell}.SliceCoord(ZIndex(icell),:);
+            Z=SliceCoord(3); %to generalize for non-parallel planes
+            if isfield(CalibIn{icell},'SliceAngle')
+            norm_plane=angle2normal(CalibIn{icell}.SliceAngle);
+            Z=Z-(norm_plane(1)*(X-SliceCoord(1))+norm_plane(2)*(Y-SliceCoord(2)))/norm_plane(3); 
+            end
+        end
+        [XIMA,YIMA]=px_XYZ(CalibIn{icell},X,Y,Z);%corresponding image indices for each point in the real space grid
+    end
+    Dx=(coord_x{icell}(end)-coord_x{icell}(1))/(npx(icell)-1);
+    Dy=(coord_y{icell}(end)-coord_y{icell}(1))/(npy(icell)-1);
+    indx_ima=1+round((XIMA-coord_x{icell}(1))/Dx);%indices of the initial matrix close to the points of the new grid
+    indy_ima=1+round((YIMA-coord_y{icell}(1))/Dy);
+    Delta_x=1+(XIMA-coord_x{icell}(1))/Dx-indx_ima;%
+    Delta_y=1+(YIMA-coord_y{icell}(1))/Dy-indy_ima;
+    XIMA=reshape(indx_ima,1,[]);%indices reorganized in 'line'
+    YIMA=reshape(indy_ima,1,[]);%indices reorganized in 'line'
+    flagin=XIMA>=1 & XIMA<=npx(icell) & YIMA >=1 & YIMA<=npy(icell);%flagin=1 inside the original image
     siz=size(A{icell});
-    [XIMA,YIMA]=px_XYZ(CalibIn{icell},X,Y,zphys);%corresponding image indices for each point in the real space grid
-    XIMA=reshape(round(XIMA),1,npx*npy);%indices reorganized in 'line'
-    YIMA=reshape(round(YIMA),1,npx*npy);
-    flagin=XIMA>=1 & XIMA<=npx & YIMA >=1 & YIMA<=npy;%flagin=1 inside the original image
+    checkuint8=isa(A{icell},'uint8');%check for image input with 8 bits
+    checkuint16=isa(A{icell},'uint8');%check for image input with 16 bits
+    A{icell}=double(A{icell});
     if numel(siz)==2 %(B/W images)
-        vec_A=reshape(A{icell}(:,:,1),1,npx*npy);%put the original image in line
+        vec_A=reshape(A{icell}(:,:,1),1,[]);%put the original image in line
         ind_in=find(flagin);
         ind_out=find(~flagin);
-        ICOMB=((XIMA-1)*npy+(npy+1-YIMA));
+        ICOMB=((XIMA-1)*npy(icell)+(npy(icell)+1-YIMA));
         ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A
         vec_B(ind_in)=vec_A(ICOMB);
         vec_B(ind_out)=zeros(size(ind_out));
-        A_out{icell}=reshape(vec_B,npy,npx);%new image in real coordinates
+        A_out{icell}=reshape(vec_B,np_theta,np_r);%new image in real coordinates
+        DA_y=circshift(A_out{icell},-1,1)-A_out{icell};
+        DA_y(end,:)=0;
+        DA_x=circshift(A_out{icell},-1,2)-A_out{icell};
+        DA_x(:,end)=0;
+        A_out{icell}=A_out{icell}+Delta_x.*DA_x+Delta_y.*DA_y;%linear interpolation
     else
         for icolor=1:siz(3)
-                vec_A=reshape(A{icell}(:,:,icolor),1,npx*npy);%put the original image in line
-                ind_in=find(flagin);
-                ind_out=find(~flagin);
-                ICOMB=((XIMA-1)*npy+(npy+1-YIMA));
-                ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A
-                vec_B(ind_in)=vec_A(ICOMB);
-                vec_B(ind_out)=zeros(size(ind_out));
-                A_out{icell}(:,:,icolor)=reshape(vec_B,npy,npx);%new image in real coordinates
-        end
-    end
-end
-
+            vec_A=reshape(A{icell}(:,:,icolor),1,[]);%put the original image in line
+            ind_in=find(flagin);
+            ind_out=find(~flagin);
+            ICOMB=((XIMA-1)*npy(icell)+(npy(icell)+1-YIMA));
+            ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A
+            vec_B(ind_in)=vec_A(ICOMB);
+            vec_B(ind_out)=zeros(size(ind_out));
+            A_out{icell}(:,:,icolor)=reshape(vec_B,np_theta,np_r);%new image in real coordinates
+            DA_y=circshift(A_out{icell}(:,:,icolor),-1,1)-A_out{icell}(:,:,icolor);
+            DA_y(end,:)=0;
+            DA_x=circshift(A_out{icell}(:,:,icolor),-1,2)-A_out{icell}(:,:,icolor);
+            DA_x(:,end)=0;
+            A_out{icell}(:,:,icolor)=A_out{icell}(:,:,icolor)+Delta_x.*DA_x+Delta_y.*DA_y;%linear interpolation
+        end
+    end
+    if checkuint8
+        A_out{icell}=uint8(A_out{icell});
+    elseif checkuint16
+        A_out{icell}=uint16(A_out{icell});
+    end
+end
+

trunk/src/uvmat.m

@@ -640 +640 @@
 
 %------------------------------------------------------------------------
+% function called by the upper bar menu item Export/make movie
 % --------------------------------------------------------------------
 function MenuExportMovie_Callback(hObject, eventdata, handles)
@@ -645 +646 @@
 set(handles.MenuExportMovie,'BusyAction','queue')% activate the button
 
-[RootPath,SubDir,RootFile,FileIndex,FileExt]=read_file_boxes(handles);
-FileBase=fullfile(RootPath,RootFile);
+[RootPath,SubDir,RootFile]=read_file_boxes(handles);%read input file path from the GUI uvmat
 
 %% create a fig and axis for movies
-figure_movie=findobj(allchild(0),'name','figure_movie');
-
+figure_movie=findobj(allchild(0),'name','figure_movie');% find existing movie figure
 if ~isempty(figure_movie)
     delete(figure_movie)%delete existing figure_movie
 end
-figure_movie=figure;
+figure_movie=figure;% create a new movie figure 
 nbpix=[640 480];% resolution VGA
 set(figure_movie,'name','figure_movie','Position',[1 1 nbpix])
@@ -671 +670 @@
 
 %% create the GUI set_movie
-%set(0,'Units','points')
-%ScreenSize=get(0,'ScreenSize');% get the size of the screen, to put the fig on the upper right
 Position=get(figure_movie,'Position');
 Position(2)=Position(2)+1.2*Position(4);
@@ -678 +675 @@
 Position(4)=Position(4)/2;
 hfig=findobj(allchild(0),'Tag','set_movie');
-if ~isempty(hfig),delete(hfig), end; %delete existing version of the GUI
+if ~isempty(hfig),delete(hfig), end %delete existing version of the GUI
 hfig=figure('name','set_movie','tag','set_movie','MenuBar','none','NumberTitle','off','Units','pixels',...
     'Position',Position);
@@ -787 +784 @@
     end
 end
-%create avi open
-%aviobj=avifile(MovieName,'Compression','None','fps',fps);
-
+
+%% create movie
+% duration = 2*FrameNumber/60;
+% [~,name,~] = fileparts(MovieName);
+% ffmpegcmd = ['ffmpeg -i ' ' ' MovieName ' ' '-filter:v "setpts=(',...
+%      num2str(fps/duration),')*PTS"' ' ' strcat(MovieDir,...
+%      strcat('/',name,'.mkv'))];
+% [ffmpeg_err,~] = system(ffmpegcmd);
+% if ffmpeg_err
+%     disp_uvmat('ERROR',['ERROR: Errors in conversion to mkv, close MATLAB and run "LD_PRELOAD=/usr/lib/x86_64-linux-gnu/libstdc++.so.6 matlab" in terminal'],1 )
+%     return
+% end
+% msgbox_uvmat('CONFIRMATION',{['movie ' MovieName ' created '];['with ' num2str(FrameNumber) ' frames']})
+
+
+% 
 aviobj = VideoWriter(MovieName,'Motion JPEG AVI');
 open(aviobj)
@@ -1180 +1190 @@
     'Callback',@(hObject,eventdata)set_slice_CheckRefraction_Callback(hObject,eventdata),...
     'String','refraction','Value',CheckRefraction,'FontUnits','points','FontSize',12,'FontWeight','bold','TooltipString','''CheckRefraction'':=1 to provide refraction correction');
-uicontrol('Style','text','Units','normalized', 'Position', [2*ii+2*ww 0.95-3*ii-1.75*hh ww hh/2],'BackgroundColor',BackgroundColor,'Tag','Refraction_title',...
+uicontrol('Style','text','Units','normalized', 'Position', [2*ii+2*ww 0.95-3*ii-1.7*hh ww hh/2],'BackgroundColor',BackgroundColor,'Tag','Refraction_title',...
     'String','index','Visible','off','FontUnits','points','FontSize',12,'FontWeight','bold','ForegroundColor','blue','HorizontalAlignment','right');%title
 uicontrol('Style','edit','Units','normalized', 'Position', [4*ii+3*ww 0.95-3*ii-2*hh ww hh],'tag','num_RefractionIndex','BackgroundColor',[1 1 1],...
     'String',num2str(RefractionIndex),'Visible','off','FontUnits','points','FontSize',12,'FontWeight','bold','TooltipString','''num_RefractionIndex'': refraction index of water');
 %  raw 4 of the GUI
-uicontrol('Style','text','Units','normalized', 'Position', [ii 0.95-4*ii-3.25*hh ww hh],'BackgroundColor',BackgroundColor,'Tag','NbSlice_title',...
+uicontrol('Style','text','Units','normalized', 'Position', [ii 0.95-4*ii-3.0*hh ww hh],'BackgroundColor',BackgroundColor,'Tag','NbSlice_title',...
     'String','NbSlice','FontUnits','points','FontSize',12,'FontWeight','bold','ForegroundColor','blue','HorizontalAlignment','right');%title
-uicontrol('Style','edit','Units','normalized', 'Position', [2*ii+ww 0.95-4*ii-3*hh ww hh],'tag','num_NbSlice','BackgroundColor',[1 1 1],...
+uicontrol('Style','edit','Units','normalized', 'Position', [2*ii+ww 0.95-4*ii-2.8*hh ww hh],'tag','num_NbSlice','BackgroundColor',[1 1 1],...
     'String',num2str(NbSlice),'FontUnits','points','FontSize',12,'FontWeight','bold','TooltipString','''num_NbSlice'':number of slices');%edit box
-uicontrol('Style','checkbox','Units','normalized', 'Position', [3*ii+2*ww 0.95-4*ii-3*hh 2*ww hh],'tag','CheckVolumeScan','BackgroundColor',BackgroundColor,...
+uicontrol('Style','checkbox','Units','normalized', 'Position', [3*ii+2*ww 0.95-4*ii-2.7*hh 2*ww hh],'tag','CheckVolumeScan','BackgroundColor',BackgroundColor,...
     'String','volume scan','Value',CheckVolumeScan,'FontUnits','points','FontSize',12,'FontWeight','bold','TooltipString','''CheckVolumeScan'':=1 for volume scan (z varies with j index)');
 %  raw 5 of the GUI
-uicontrol('Style','text','Units','normalized', 'Position', [2*ii+2*ww 0.95-3*ii-3.5*hh ww hh/2],'BackgroundColor',BackgroundColor,...
-    'String','first','FontUnits','points','FontSize',12,'FontWeight','bold','ForegroundColor','blue','HorizontalAlignment','center');%title
-uicontrol('Style','text','Units','normalized', 'Position', [3*ii+3*ww 0.95-3*ii-3.5*hh ww hh/2],'BackgroundColor',BackgroundColor,...
-    'String','last','FontUnits','points','FontSize',12,'FontWeight','bold','ForegroundColor','blue','HorizontalAlignment','center');%title
-uicontrol('Style','text','Units','normalized', 'Position', [ii 0.95-5*ii-4*hh 2*ww hh/2],'BackgroundColor',BackgroundColor,'Tag','Angle_title_1',...
-    'String','tild angle x axis','FontUnits','points','FontSize',12,'FontWeight','bold','ForegroundColor','blue','HorizontalAlignment','center');%title
-uicontrol('Style','text','Units','normalized', 'Position', [ii 0.95-6*ii-5*hh 2*ww hh/2],'BackgroundColor',BackgroundColor,'Tag','Angle_title_2',...
-    'String','tild angle y axis','FontUnits','points','FontSize',12,'FontWeight','bold','ForegroundColor','blue','HorizontalAlignment','center');%title
+uicontrol('Style','text','Units','normalized', 'Position', [2*ii+1*ww 0.95-2*ii-3.9*hh ww hh],'BackgroundColor',BackgroundColor,...
+    'String','origin','FontUnits','points','FontSize',12,'FontWeight','bold','ForegroundColor','blue','HorizontalAlignment','center');%title
+uicontrol('Style','text','Units','normalized', 'Position', [2*ii+2*ww 0.95-2*ii-3.5*hh ww hh],'BackgroundColor',BackgroundColor,...
+    'String',{'first';'angle'},'FontUnits','points','FontSize',12,'FontWeight','bold','ForegroundColor','blue','HorizontalAlignment','center');%title
+uicontrol('Style','text','Units','normalized', 'Position', [3*ii+3*ww 0.95-2*ii-3.5*hh ww hh],'BackgroundColor',BackgroundColor,...
+    'String',{'last';'angle'},'FontUnits','points','FontSize',12,'FontWeight','bold','ForegroundColor','blue','HorizontalAlignment','center');%title
+
+uicontrol('Style','edit','Units','normalized', 'Position', [3*ii+ww 0.95-5*ii-4.2*hh ww hh],'tag','num_SliceCoord_1','BackgroundColor',[1 1 1],...
+    'String',num2str(SliceCoord(1)),'FontUnits','points','FontSize',12,'FontWeight','bold','TooltipString','''num_SliceCoord_1'':x position of the tild origin');%edit box
+uicontrol('Style','edit','Units','normalized', 'Position', [3*ii+ww 0.95-6*ii-5.2*hh ww hh],'tag','num_SliceCoord_2','BackgroundColor',[1 1 1],...
+    'String',num2str(SliceCoord(2)),'FontUnits','points','FontSize',12,'FontWeight','bold','TooltipString','''num_SliceCoord_2'':y position of the tild origin');%edit box
+
+uicontrol('Style','text','Units','normalized', 'Position', [ii 0.95-5*ii-4*hh 1.3*ww hh/2],'BackgroundColor',BackgroundColor,'Tag','Angle_title_1',...
+    'String','tild x axis','FontUnits','points','FontSize',12,'FontWeight','bold','ForegroundColor','blue','HorizontalAlignment','center');%title
+uicontrol('Style','text','Units','normalized', 'Position', [ii 0.95-6*ii-5*hh 1.3*ww hh/2],'BackgroundColor',BackgroundColor,'Tag','Angle_title_2',...
+    'String','tild y axis','FontUnits','points','FontSize',12,'FontWeight','bold','ForegroundColor','blue','HorizontalAlignment','center');%title
 %  raw 6 of the GUI
 uicontrol('Style','edit','Units','normalized', 'Position', [3*ii+2*ww 0.95-5*ii-4.2*hh ww hh],'tag','num_SliceAngle_1_1','BackgroundColor',[1 1 1],...
@@ -1270 +1288 @@
 end
 GeometryCalib.SliceCoord=Z_plane'*[0 0 1];
+GeometryCalib.SliceCoord(:,1)=SliceData.SliceCoord(1);
+GeometryCalib.SliceCoord(:,2)=SliceData.SliceCoord(2);
 GeometryCalib.SliceAngle=zeros(GeometryCalib.NbSlice,3);
 Angle_1=linspace(SliceData.SliceAngle_1(1),SliceData.SliceAngle_1(2),SliceData.NbSlice);
@@ -3554 +3574 @@
 end
 switch UvData.FileInfo{1}.FieldType
-    case {'civdata','netcdf'};
+    case {'civdata','netcdf'}
         list_fields=get(handles.FieldName,'String');% list menu fields
         FieldName= list_fields{get(handles.FieldName,'Value')}; % selected field
@@ -3687 +3707 @@
             end
     end
-    %     switch UvData.FileType{2}
-    %         case {'civx','civdata','netcdf'};
-    %             list_fields=get(handles.FieldName_1,'String');% list menu fields
-    %             if ischar(list_fields),list_fields={list_fields};end
-    %             FieldName_1= list_fields{get(handles.FieldName_1,'Value')}; % selected field
-    %             if ~strcmp(FieldName_1,'get_field...')
-    %                 if get(handles.FixVelType,'Value')
-    %                     VelTypeList=get(handles.VelType_1,'String');
-    %                     VelType_1=VelTypeList{get(handles.VelType_1,'Value')};% read the velocity type.
-    %                 end
-    %             end
-    %             if isempty(FieldName_1)
-    %                 FieldName_1=FieldName;% if blank reproduce the field name of the first field
-    %             end
-    %             if ~isempty(regexp(FieldName_1,'^vel', 'once'))&& strcmp(get(handles.ColorCode,'Visible'),'on')
-    %                 list_code=get(handles.ColorCode,'String');% list menu fields
-    %                 index_code=get(handles.ColorCode,'Value');% selected string index
-    %                 if  ~strcmp(list_code{index_code},'black') &&  ~strcmp(list_code{index_code},'white')
-    %                     list_code=get(handles.ColorScalar,'String');% list menu fields
-    %                     index_code=get(handles.ColorScalar,'Value');% selected string index
-    %                     ParamIn_1.ColorVar= list_code{index_code}; % selected field for vector color display
-    %                 end
-    %             end
-    %         case {'video','mmreader','cine_phantom'}
-    %             ParamIn_1=UvData.MovieObject{2};
-    %             if ~strcmp(NomType_1,'*')
-    %                 frame_index_1=j1_1;%frame index for movies or multimage
-    %             else
-    %                 frame_index_1=i1_1;
-    %             end
-    %          case 'multimage'
-    %             if strcmp(NomType_1,'*')%frame index for movies or multimage
-    %                 frame_index_1=i1_1;
-    %             else
-    %                 frame_index_1=j1_1;
-    %             end
-    %         case 'vol' %TODO: update
-    %             if isfield(UvData.XmlData,'Npy') && isfield(UvData.XmlData,'Npx')
-    %                 ParamIn_1.Npy=UvData.XmlData.Npy;
-    %                 ParamIn_1.Npx=UvData.XmlData.Npx;
-    %             else
-    %                 errormsg='Npx and Npy need to be defined in the xml file for volume images .vol';
-    %                 return
-    %             end
-    %     end
     
     test_keepdata_1=0;% test for keeping the previous stored data if the input files are unchanged
@@ -3951 +3926 @@
     if numel(Field)==2
         UvData.Field=sub_field(Field{1},[],Field{2});
-       % UvData.Field.(FieldName)=Field{1}.(FieldName)-Field{2}.(FieldName_1);
     end      
 else
@@ -4024 +3998 @@
     UvData.Field=tps_coeff_field(UvData.Field,check_proj_tps);
 
-%% reset the min and max of scalar if only the mask is displayed(TODO: check the need)
-% if isfield(UvData,'Mask')&& ~isfield(UvData,'A')
-%     set(handles.num_MinA,'String','0')
-%     set(handles.num_MaxA,'String','255')
-% end
-    set(handles.Objects,'Visible','on')
-
-    %% Plot the projections on the selected  projection objects
+   
+%% Plot the projections on the selected  projection objects
+     set(handles.Objects,'Visible','on')
     %if no projection object exists, create a default one
     if isempty(UvData.ProjObject{1})
@@ -5461 +5430 @@
 %------------------------------------------------------------------------
 set(handles.CheckFixLimits,'Value',1) %suppress auto mode
-% set(handles.CheckFixLimits,'BackgroundColor',[1 1 0])
 update_plot(handles);
 
@@ -5471 +5439 @@
 %------------------------------------------
 set(handles.CheckFixScalar,'Value',1) %suppress auto mode
-% set(handles.CheckFixScalar,'BackgroundColor',[1 1 0])
 MinA=str2double(get(handles.num_MinA,'String'));
 MaxA=str2double(get(handles.num_MaxA,'String'));