Changeset 206


Ignore:
Timestamp:
Feb 27, 2011, 10:40:29 PM (14 years ago)
Author:
sommeria
Message:

bug fixes to deal with volumes, storage of ACTION menu in series fixed

Location:
trunk/src
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/plot_field.m

    r201 r206  
    159159        return
    160160    end
     161    index_2D=find(NbDim==2,2);%find 2D fields (at most 2)
    161162    index_3D=find(NbDim>2,1);
    162163    if ~isempty(index_3D)
     164        if isfield(Data,'NbDim')&& isequal(Data.NbDim,2)
     165            index_2D=[index_2D index_3D];
     166        else
    163167        msgbox_uvmat('ERROR','volume plot not implemented yet');
    164168        return
    165     end
    166     index_2D=find(NbDim==2,2);%find 2D fields (at most 2)
     169        end
     170    end
     171   
    167172    index_1D=find(NbDim==1);
    168173    index_0D=find(NbDim==0);
    … …  
    243248set(haxes,'UserData',AxeData)
    244249
    245 %% update the parameters stored in parent figure
     250%% update the plotted field stored in parent figure
    246251FigData=get(hfig,'UserData');
    247252tagaxes=get(haxes,'tag');
    … …  
    562567        VarType.coord=VarType.coord(ind_coord);
    563568    end
    564 %     idim_Y=[]; 
    565 %     test_grid=0;
    566569    if ~isempty(ivar_U) && ~isempty(ivar_V)% vector components detected
    567570        if test_vec
    … …  
    725728            siz=3;%color image
    726729        else
    727             errormsg=['unrecognized scalar type: ' num2str(np(3)) ' color components'];
     730            errormsg=['unrecognized scalar type in plot_field: considered as 2D field with ' num2str(np(3)) ' color components'];
    728731            return
    729732        end
    … …  
    748751            PlotParam.Scalar.MaxA=[];%default
    749752        end
     753        Aline=[];
    750754        if isequal(PlotParam.Scalar.FixScal,0)||isempty(PlotParam.Scalar.MinA)||~isa(PlotParam.Scalar.MinA,'double')  %correct if there is no numerical data in edit box
    751             MinA=double(min(min(A)));
     755            Aline=reshape(A,1,[]);
     756            Aline=Aline(~isnan(A));
     757            if isempty(Aline)
     758                 errormsg=['NaN input scalar or image in plot_field'];
     759                return
     760            end
     761            MinA=double(min(Aline))
     762            %MinA=double(min(min(A)));
    752763        else
    753             MinA=PlotParam.Scalar.MinA; 
     764            MinA=PlotParam.Scalar.MinA;
    754765        end;
    755766        if isequal(PlotParam.Scalar.FixScal,0)||isempty(PlotParam.Scalar.MaxA)||~isa(PlotParam.Scalar.MaxA,'double') %correct if there is no numerical data in edit box
    756             MaxA=double(max(max(A)));
     767            if isempty(Aline)
     768               Aline=reshape(A,1,[]);
     769               Aline=Aline(~isnan(A));
     770               if isempty(Aline)
     771                 errormsg=['NaN input scalar or image in plot_field'];
     772                return
     773               end
     774            end
     775            MaxA=double(max(Aline))
     776           % MaxA=double(max(max(A)));
    757777        else
    758778            MaxA=PlotParam.Scalar.MaxA; 

  • trunk/src/proj_field.m

    r204 r206  
    8181
    8282function [ProjData,errormsg]=proj_field(FieldData,ObjectData)
    83 errormsg=[];%default
     83errormsg='';%default
     84%% case of no projection (object is used only as graph display)
    8485if isfield(ObjectData,'ProjMode') && (isequal(ObjectData.ProjMode,'none')||isequal(ObjectData.ProjMode,'mask_inside')||isequal(ObjectData.ProjMode,'mask_outside'))
    8586    ProjData=[];
    8687    return
    8788end
    88 %introduce default field properties (reading old standards)
     89
     90%% in the absence of object Style or projection mode, or object coordinaes, the input field is just tranfered without change
    8991if ~isfield(ObjectData,'Style')||~isfield(ObjectData,'ProjMode')
    9092    ProjData=FieldData;
    … …  
    736738            elseif isequal(ProjMode,'interp') %linear interpolation:
    737739                npoint=floor(linelength/DX)+1;% nbre of points in the profile (interval DX)
    738                 Xproj=[linelength/(2*npoint):linelength/npoint:linelength-linelength/(2*npoint)];
     740                Xproj=linelength/(2*npoint):linelength/npoint:linelength-linelength/(2*npoint);
    739741                xreg=cos(theta(ip))*Xproj+ObjectData.Coord(ip,1);
    740742                yreg=sin(theta(ip))*Xproj+ObjectData.Coord(ip,2);
    … …  
    942944
    943945%rotation angles
    944 Phi=0;%default
    945 Theta=0;
    946 Psi=0;
    947 if isfield(ObjectData,'Phi')&& ~isempty(ObjectData.Phi)
    948     Phi=(pi/180)*ObjectData.Phi;%first Euler angle in radian
    949 end
    950 if isfield(ObjectData,'Theta')&& ~isempty(ObjectData.Theta)
    951     Theta=(pi/180)*ObjectData.Theta;%second Euler angle in radian
    952 end
    953 if isfield(ObjectData,'Psi')&& ~isempty(ObjectData.Psi)
    954     Psi=(pi/180)*ObjectData.Psi;%third Euler angle in radian
    955 end
     946PlaneAngle=[0 0 0];
     947norm_plane=[0 0 1];
     948cos_om=1;
     949sin_om=0;
     950if isfield(ObjectData,'Angle')&& isequal(size(ObjectData.Angle),[1 3])&& ~isequal(ObjectData.Angle,[0 0 0])
     951    PlaneAngle=ObjectData.Angle;
     952    om=norm(PlaneAngle);%norm of rotation angle in radians
     953    OmAxis=PlaneAngle/om; %unit vector marking the rotation axis
     954    cos_om=cos(pi*om/180);
     955    sin_om=sin(pi*om/180);
     956    coeff=OmAxis(3)*(1-cos_om);
     957    %components of the unity vector norm_plane normal to the projection plane
     958    norm_plane(1)=OmAxis(1)*coeff+OmAxis(2)*sin_om;
     959    norm_plane(2)=OmAxis(2)*coeff-OmAxis(1)*sin_om;
     960    norm_plane(3)=OmAxis(3)*coeff+cos_om;
     961end
     962testangle=~isequal(PlaneAngle,[0 0 0]);
     963% Phi=0;%default
     964% Theta=0;
     965% Psi=0;
     966% if isfield(ObjectData,'Phi')&& ~isempty(ObjectData.Phi)
     967%     Phi=(pi/180)*ObjectData.Phi;%first Euler angle in radian
     968% end
     969% if isfield(ObjectData,'Theta')&& ~isempty(ObjectData.Theta)
     970%     Theta=(pi/180)*ObjectData.Theta;%second Euler angle in radian
     971% end
     972% if isfield(ObjectData,'Psi')&& ~isempty(ObjectData.Psi)
     973%     Psi=(pi/180)*ObjectData.Psi;%third Euler angle in radian
     974% end
    956975
    957976%components of the unity vector normal to the projection plane
    958 NormVec_X=-sin(Phi)*sin(Theta);
    959 NormVec_Y=cos(Phi)*sin(Theta);
    960 NormVec_Z=cos(Theta);
     977% NormVec_X=-sin(Phi)*sin(Theta);
     978% NormVec_Y=cos(Phi)*sin(Theta);
     979% NormVec_Z=cos(Theta);
    961980
    962981%mesh sizes DX and DY
    … …  
    10171036%-----------------------------------------------------------------
    10181037idimvar=0;
     1038
    10191039[CellVarIndex,NbDimVec,VarTypeCell,errormsg]=find_field_indices(FieldData);
    10201040if ~isempty(errormsg)
    … …  
    10751095        if length(ivar_Z)==1 &&  width > 0
    10761096            %components of the unitiy vector normal to the projection plane
    1077             fieldZ=NormVec_X*coord_x + NormVec_Y*coord_y+ NormVec_Z*coord_z;% distance to the plane           
     1097            fieldZ=norm_plane(1)*coord_x + norm_plane(2)*coord_y+ norm_plane(3)*coord_z;% distance to the plane           
    10781098            indcut=find(abs(fieldZ) <= width);
     1099            size(indcut)
    10791100            for ivar=VarIndex
    10801101                VarName=FieldData.ListVarName{ivar};
    … …  
    10881109
    10891110       %rotate coordinates if needed
    1090         if isequal(Phi,0)
     1111        if testangle
     1112%             coord_X=coord_x;
     1113%             coord_Y=coord_y;
     1114%             if ~isequal(Theta,0)
     1115%                 coord_Y=coord_Y *cos(Theta);
     1116%             end
     1117%         else
     1118            coord_X=(coord_x *cos(Phi) + coord_y* sin(Phi));
     1119            coord_Y=(-coord_x *sin(Phi) + coord_y *cos(Phi))*cos(Theta);
     1120%         end
     1121%         if ~isempty(ivar_Z)
     1122            coord_Y=coord_Y+coord_z *sin(Theta);
     1123%         end
     1124%         if testangle
     1125                coord_X=(coord_X *cos(Psi) - coord_Y* sin(Psi));%A VERIFIER
     1126                coord_Y=(coord_X *sin(Psi) + coord_Y* cos(Psi));
     1127        else
    10911128            coord_X=coord_x;
    10921129            coord_Y=coord_y;
    1093             if ~isequal(Theta,0)
    1094                 coord_Y=coord_Y *cos(Theta);
    1095             end
    1096         else
    1097             coord_X=(coord_x *cos(Phi) + coord_y* sin(Phi));
    1098             coord_Y=(-coord_x *sin(Phi) + coord_y *cos(Phi))*cos(Theta);
    1099         end
    1100         if ~isempty(ivar_Z)
    1101             coord_Y=coord_Y+coord_z *sin(Theta);
    1102         end
    1103         if ~isequal(Psi,0)
    1104                 coord_X=(coord_X *cos(Psi) - coord_Y* sin(Psi));%A VERIFIER
    1105                 coord_Y=(coord_X *sin(Psi) + coord_Y* cos(Psi));
    11061130        end
    11071131       
    … …  
    11631187            end 
    11641188        elseif isequal(ObjectData.ProjMode,'interp')||isequal(ObjectData.ProjMode,'filter')%interpolate data on a regular grid
    1165             coord_x_proj=[XMin:DX:XMax];
    1166             coord_y_proj=[YMin:DY:YMax];
     1189            coord_x_proj=XMin:DX:XMax;
     1190            coord_y_proj=YMin:DY:YMax;
    11671191            DimCell={'coord_y','coord_x'};
    11681192            ProjData.ListVarName={'coord_y','coord_x'};
    … …  
    13051329        npX=1+round(abs(Coord{NbDim}(end)-Coord{NbDim}(1))/DX);%nbre of points after interpol
    13061330        for idim=[1:NbDim]
     1331            if test_interp(idim)
     1332                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
     1333            end
     1334        end       
     1335        Coord_y=linspace(Coord{NbDim-1}(1),Coord{NbDim-1}(end),npY);
     1336        test_direct_y=test_direct(NbDim-1);
     1337        Coord_x=linspace(Coord{NbDim}(1),Coord{NbDim}(end),npX);
     1338        test_direct_x=test_direct(NbDim);
     1339        DAX=DCoord_min(NbDim);
     1340        DAY=DCoord_min(NbDim-1); 
     1341        minAX=min(Coord_x);
     1342        maxAX=max(Coord_x);
     1343        minAY=min(Coord_y);
     1344        maxAY=max(Coord_y);
     1345        xcorner=[minAX maxAX minAX maxAX]-ObjectData.Coord(1,1);
     1346        ycorner=[maxAY maxAY minAY minAY]-ObjectData.Coord(1,2);
     1347        xcor_new=xcorner*cos_om+ycorner*sin_om;%coord new frame
     1348        ycor_new=-xcorner*sin_om+ycorner*cos_om;
     1349        if ~testXMax
     1350            XMax=max(xcor_new);
     1351        end
     1352        if ~testXMin
     1353            XMin=min(xcor_new);
     1354        end
     1355        if ~testYMax
     1356            YMax=max(ycor_new);
     1357        end
     1358        if ~testYMin
     1359            YMin=min(ycor_new);
     1360        end
     1361        DXinit=(maxAX-minAX)/(DimValue(NbDim)-1);
     1362        DYinit=(maxAY-minAY)/(DimValue(NbDim-1)-1);
     1363        if DX==0
     1364            DX=DXinit;
     1365        end
     1366        if DY==0
     1367            DY=DYinit;
     1368        end
     1369        if NbDim==3
     1370            DZ=(Coord{1}(end)-Coord{1}(1))/(DimValue(1)-1);
     1371            if ~test_direct(1)
     1372                DZ=-DZ;
     1373            end
     1374            Coord_z=linspace(Coord{1}(1),Coord{1}(end),DimValue(1));
     1375            test_direct_z=test_direct(1);
     1376        end
     1377        npX=floor((XMax-XMin)/DX+1);
     1378        npY=floor((YMax-YMin)/DY+1);   
     1379        if test_direct_y
     1380            coord_y_proj=linspace(YMin,YMax,npY);%abscissa of the new pixels along the line
     1381        else
     1382            coord_y_proj=linspace(YMax,YMin,npY);%abscissa of the new pixels along the line
     1383        end
     1384        if test_direct_x
     1385            coord_x_proj=linspace(XMin,XMax,npX);%abscissa of the new pixels along the line
     1386        else
     1387            coord_x_proj=linspace(XMax,XMin,npX);%abscissa of the new pixels along the line
     1388        end
     1389       
     1390        % case with no rotation and interpolation
     1391        if isequal(ProjMode,'projection') && ~testangle
     1392            if ~testXMin && ~testXMax && ~testYMin && ~testYMax && NbDim==2
     1393                ProjData=FieldData;
     1394            else
     1395                indY=NbDim-1;
     1396                if test_direct(indY)
     1397                    min_indy=ceil((YMin-Coord{indY}(1))/DYinit)+1;
     1398                    max_indy=floor((YMax-Coord{indY}(1))/DYinit)+1;
     1399                    Ybound(1)=Coord{indY}(1)+DYinit*(min_indy-1);
     1400                    Ybound(2)=Coord{indY}(1)+DYinit*(max_indy-1);
     1401                else
     1402                    min_indy=ceil((Coord{indY}(1)-YMax)/DYinit)+1;
     1403                    max_indy=floor((Coord{indY}(1)-YMin)/DYinit)+1;
     1404                    Ybound(2)=Coord{indY}(1)-DYinit*(max_indy-1);
     1405                    Ybound(1)=Coord{indY}(1)-DYinit*(min_indy-1);
     1406                end   
     1407                if test_direct(NbDim)==1
     1408                    min_indx=ceil((XMin-Coord{NbDim}(1))/DXinit)+1;
     1409                    max_indx=floor((XMax-Coord{NbDim}(1))/DXinit)+1;
     1410                    Xbound(1)=Coord{NbDim}(1)+DXinit*(min_indx-1);
     1411                    Xbound(2)=Coord{NbDim}(1)+DXinit*(max_indx-1);
     1412                else
     1413                    min_indx=ceil((Coord{NbDim}(1)-XMax)/DXinit)+1;
     1414                    max_indx=floor((Coord{NbDim}(1)-XMin)/DXinit)+1;
     1415                    Xbound(2)=Coord{NbDim}(1)+DXinit*(max_indx-1);
     1416                    Xbound(1)=Coord{NbDim}(1)+DXinit*(min_indx-1);
     1417                end
     1418                if NbDim==3
     1419                    DimCell(1)=[]; %suppress z variable
     1420                    DimValue(1)=[];
     1421                                        %structured coordinates
     1422                    if test_direct(1)
     1423                        iz=ceil((ObjectData.Coord(1,3)-Coord{1}(1))/DZ)+1;
     1424                    else
     1425                        iz=ceil((Coord{1}(1)-ObjectData.Coord(1,3))/DZ)+1;
     1426                    end
     1427                end
     1428                min_indy=max(min_indy,1);% deals with margin (bound lower than the first index)
     1429                min_indx=max(min_indx,1);
     1430                max_indy=min(max_indy,DimValue(1));
     1431                max_indx=min(max_indx,DimValue(2));
     1432                for ivar=VarIndex% loop on non coordinate variables
     1433                    VarName=FieldData.ListVarName{ivar};
     1434                    ProjData.ListVarName=[ProjData.ListVarName VarName];
     1435                    ProjData.VarDimName=[ProjData.VarDimName {DimCell}];
     1436                    if isfield(FieldData,'VarAttribute') && length(FieldData.VarAttribute)>=ivar
     1437                        ProjData.VarAttribute{length(ProjData.ListVarName)}=FieldData.VarAttribute{ivar};
     1438                    end
     1439                    if NbDim==3
     1440                        eval(['ProjData.' VarName '=squeeze(FieldData.' VarName '(iz,min_indy:max_indy,min_indx:max_indx));']);
     1441                    else
     1442                        eval(['ProjData.' VarName '=FieldData.' VarName '(min_indy:max_indy,min_indx:max_indx,:);']);
     1443                    end
     1444                end 
     1445                eval(['ProjData.' AYName '=[Ybound(1) Ybound(2)];']) %record the new (projected ) y coordinates
     1446                eval(['ProjData.' AXName '=[Xbound(1) Xbound(2)];']) %record the new (projected ) x coordinates
     1447            end
     1448        else       % case with rotation and/or interpolation
     1449            if NbDim==2 %2D case
     1450                [X,Y]=meshgrid(coord_x_proj,coord_y_proj);%grid in the new coordinates
     1451                XIMA=ObjectData.Coord(1,1)+(X)*cos(Phi)-Y*sin(Phi);%corresponding coordinates in the original image
     1452                YIMA=ObjectData.Coord(1,2)+(X)*sin(Phi)+Y*cos(Phi);
     1453                XIMA=(XIMA-minAX)/DXinit+1;% image index along x
     1454                YIMA=(-YIMA+maxAY)/DYinit+1;% image index along y
     1455                XIMA=reshape(round(XIMA),1,npX*npY);%indices reorganized in 'line'
     1456                YIMA=reshape(round(YIMA),1,npX*npY);
     1457                flagin=XIMA>=1 & XIMA<=DimValue(2) & YIMA >=1 & YIMA<=DimValue(1);%flagin=1 inside the original image
     1458                if isequal(ObjectData.ProjMode,'filter')
     1459                    npx_filter=ceil(abs(DX/DAX));
     1460                    npy_filter=ceil(abs(DY/DAY));
     1461                    Mfilter=ones(npy_filter,npx_filter)/(npx_filter*npy_filter);
     1462                    test_filter=1;
     1463                else
     1464                    test_filter=0;
     1465                end
     1466                eval(['ProjData.' AYName '=[coord_y_proj(1) coord_y_proj(end)];']) %record the new (projected ) y coordinates
     1467                eval(['ProjData.' AXName '=[coord_x_proj(1) coord_x_proj(end)];']) %record the new (projected ) x coordinates
     1468                for ivar=VarIndex
     1469                    VarName=FieldData.ListVarName{ivar};
     1470                    if test_interp(1) || test_interp(2)%interpolate on a regular grid       
     1471                          eval(['ProjData.' VarName '=interp2(Coord{2},Coord{1},FieldData.' VarName ',Coord_x,Coord_y'');']) %TO TEST
     1472                    end
     1473                    %filter the field (image) if option 'filter' is used
     1474                    if test_filter 
     1475                         Aclass=class(FieldData.A);
     1476                         eval(['ProjData.' VarName '=filter2(Mfilter,FieldData.' VarName ',''valid'');'])
     1477                         if ~isequal(Aclass,'double')
     1478                             eval(['ProjData.' VarName '=' Aclass '(FieldData.' VarName ');'])%revert to integer values
     1479                         end
     1480                    end
     1481                    eval(['vec_A=reshape(FieldData.' VarName ',[],nbcolor);'])%put the original image in line             
     1482                    %ind_in=find(flagin);
     1483                    ind_out=find(~flagin);
     1484                    ICOMB=(XIMA-1)*DimValue(1)+YIMA;
     1485                    ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A
     1486                    vec_B(flagin,1:nbcolor)=vec_A(ICOMB,:);
     1487                    for icolor=1:nbcolor
     1488                        vec_B(ind_out,icolor)=zeros(size(ind_out));
     1489                    end
     1490                    ProjData.ListVarName=[ProjData.ListVarName VarName];
     1491                    ProjData.VarDimName=[ProjData.VarDimName {DimCell}];
     1492                    if isfield(FieldData,'VarAttribute')&&length(FieldData.VarAttribute)>=ivar
     1493                        ProjData.VarAttribute{length(ProjData.ListVarName)+nbcoord}=FieldData.VarAttribute{ivar};
     1494                    end     
     1495                    eval(['ProjData.' VarName '=reshape(vec_B,npY,npX,nbcolor);']);
     1496                end
     1497                ProjData.FF=reshape(~flagin,npY,npX);%false flag A FAIRE: tenir compte d'un flga antérieur 
     1498                ProjData.ListVarName=[ProjData.ListVarName 'FF'];
     1499                ProjData.VarDimName=[ProjData.VarDimName {DimCell}];
     1500                ProjData.VarAttribute{length(ProjData.ListVarName)}.Role='errorflag';
     1501            else %3D case
     1502                if ~testangle     
     1503                    % unstructured z coordinate
     1504                    test_sup=(Coord{1}>=ObjectData.Coord(1,3));
     1505                    iz_sup=find(test_sup);
     1506                    iz=iz_sup(1);
     1507                    if iz>=1 & iz<=npz
     1508                        %ProjData.ListDimName=[ProjData.ListDimName ListDimName(2:end)];
     1509                        %ProjData.DimValue=[ProjData.DimValue npY npX];
     1510                        for ivar=VarIndex
     1511                            VarName=FieldData.ListVarName{ivar};
     1512                            ProjData.ListVarName=[ProjData.ListVarName VarName];
     1513                            ProjData.VarAttribute{length(ProjData.ListVarName)}=FieldData.VarAttribute{ivar}; %reproduce the variable attributes 
     1514                            eval(['ProjData.' VarName '=squeeze(FieldData.' VarName '(iz,:,:));'])% select the z index iz
     1515                            %TODO : do a vertical average for a thick plane
     1516                            if test_interp(2) || test_interp(3)
     1517                                eval(['ProjData.' VarName '=interp2(Coord{3},Coord{2},ProjData.' VarName ',Coord_x,Coord_y'');'])
     1518                            end
     1519                        end
     1520                    end
     1521                else
     1522                    errormsg='projection of structured coordinates on oblique plane not yet implemented';
     1523                    %TODO: use interp3
     1524                    return
     1525                end
     1526            end
     1527        end
     1528    end
     1529
     1530    %% projection of  velocity components in the rotated coordinates
     1531    if testangle && length(ivar_U)==1
     1532        if isempty(ivar_V)
     1533            msgbox_uvmat('ERROR','v velocity component missing in proj_field.m')
     1534            return
     1535        end
     1536        UName=FieldData.ListVarName{ivar_U};
     1537        VName=FieldData.ListVarName{ivar_V};   
     1538        eval(['ProjData.' UName  '=cos(Phi)*ProjData.' UName '+ sin(Phi)*ProjData.' VName ';'])
     1539        eval(['ProjData.' VName  '=cos(Theta)*(-sin(Phi)*ProjData.' UName '+ cos(Phi)*ProjData.' VName ');'])
     1540        if ~isempty(ivar_W)
     1541            WName=FieldData.ListVarName{ivar_W};
     1542            eval(['ProjData.' VName '=ProjData.' VName '+ ProjData.' WName '*sin(Theta);'])%
     1543            eval(['ProjData.' WName '=NormVec_X*ProjData.' UName '+ NormVec_Y*ProjData.' VName '+ NormVec_Z* ProjData.' WName ';']);
     1544        end
     1545        if ~isequal(Psi,0)
     1546            eval(['ProjData.' UName '=cos(Psi)* ProjData.' UName '- sin(Psi)*ProjData.' VName ';']);
     1547            eval(['ProjData.' VName '=sin(Psi)* ProjData.' UName '+ cos(Psi)*ProjData.' VName ';']);
     1548        end
     1549    end
     1550end
     1551
     1552%-----------------------------------------------------------------
     1553%projection in a volume
     1554 function  [ProjData,errormsg] = proj_volume(FieldData, ObjectData)
     1555%-----------------------------------------------------------------
     1556ProjData=FieldData;%default output
     1557
     1558%% initialisation of the input parameters of the projection plane
     1559ProjMode='projection';%direct projection by default
     1560if isfield(ObjectData,'ProjMode'),ProjMode=ObjectData.ProjMode; end;
     1561
     1562%% axis origin
     1563if isempty(ObjectData.Coord)
     1564    ObjectData.Coord(1,1)=0;%origin of the plane set to [0 0] by default
     1565    ObjectData.Coord(1,2)=0;
     1566    ObjectData.Coord(1,3)=0;
     1567end
     1568
     1569%% rotation angles
     1570VolumeAngle=[0 0 0];
     1571norm_plane=[0 0 1];
     1572if isfield(ObjectData,'Angle')&& isequal(size(ObjectData.Angle),[1 3])&& ~isequal(ObjectData.Angle,[0 0 0])
     1573    PlaneAngle=ObjectData.Angle;
     1574    VolumeAngle=ObjectData.Angle;
     1575    om=norm(VolumeAngle);%norm of rotation angle in radians
     1576    OmAxis=VolumeAngle/om; %unit vector marking the rotation axis
     1577    cos_om=cos(pi*om/180);
     1578    sin_om=sin(pi*om/180);
     1579    coeff=OmAxis(3)*(1-cos_om);
     1580    %components of the unity vector norm_plane normal to the projection plane
     1581    norm_plane(1)=OmAxis(1)*coeff+OmAxis(2)*sin_om;
     1582    norm_plane(2)=OmAxis(2)*coeff-OmAxis(1)*sin_om;
     1583    norm_plane(3)=OmAxis(3)*coeff+cos_om;
     1584end
     1585testangle=~isequal(VolumeAngle,[0 0 0]);
     1586
     1587%% mesh sizes DX, DY, DZ
     1588DX=0;
     1589DY=0; %default
     1590DZ=0;
     1591if isfield(ObjectData,'DX')&~isempty(ObjectData.DX)
     1592     DX=abs(ObjectData.DX);%mesh of interpolation points
     1593end
     1594if isfield(ObjectData,'DY')&~isempty(ObjectData.DY)
     1595     DY=abs(ObjectData.DY);%mesh of interpolation points
     1596end
     1597if isfield(ObjectData,'DZ')&~isempty(ObjectData.DZ)
     1598     DZ=abs(ObjectData.DZ);%mesh of interpolation points
     1599end
     1600if  ~strcmp(ProjMode,'projection') && (DX==0||DY==0||DZ==0)
     1601        errormsg='grid mesh DX , DY or DZ is missing';
     1602        return
     1603end
     1604
     1605%% extrema along each axis
     1606testXMin=0;
     1607testXMax=0;
     1608testYMin=0;
     1609testYMax=0;
     1610if isfield(ObjectData,'RangeX')
     1611        XMin=min(ObjectData.RangeX);
     1612        XMax=max(ObjectData.RangeX);
     1613        testXMin=XMax>XMin;
     1614        testXMax=1;
     1615end
     1616if isfield(ObjectData,'RangeY')
     1617        YMin=min(ObjectData.RangeY);
     1618        YMax=max(ObjectData.RangeY);
     1619        testYMin=YMax>YMin;
     1620        testYMax=1;
     1621end
     1622width=0;%default width of the projection band
     1623if isfield(ObjectData,'RangeZ')
     1624        ZMin=min(ObjectData.RangeZ);
     1625        ZMax=max(ObjectData.RangeZ);
     1626        testZMin=ZMax>ZMin;
     1627        testZMax=1;
     1628end
     1629
     1630%% initiate Matlab  structure for physical field
     1631[ProjData,errormsg]=proj_heading(FieldData,ObjectData);
     1632ProjData.NbDim=3;
     1633ProjData.ListVarName={};
     1634ProjData.VarDimName={};
     1635if ~isequal(DX,0)&& ~isequal(DY,0)
     1636    ProjData.Mesh=sqrt(DX*DY);%define typical data mesh, useful for mouse selection in plots
     1637elseif isfield(FieldData,'Mesh')
     1638    ProjData.Mesh=FieldData.Mesh;
     1639end
     1640
     1641error=0;%default
     1642flux=0;
     1643testfalse=0;
     1644ListIndex={};
     1645
     1646%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     1647%% group the variables (fields of 'FieldData') in cells of variables with the same dimensions
     1648%-----------------------------------------------------------------
     1649idimvar=0;
     1650[CellVarIndex,NbDimVec,VarTypeCell,errormsg]=find_field_indices(FieldData);
     1651if ~isempty(errormsg)
     1652    errormsg=['error in proj_field/proj_plane:' errormsg];
     1653    return
     1654end
     1655
     1656% LOOP ON GROUPS OF VARIABLES SHARING THE SAME DIMENSIONS
     1657% CellVarIndex=cells of variable index arrays
     1658ivar_new=0; % index of the current variable in the projected field
     1659icoord=0;
     1660nbcoord=0;%number of added coordinate variables brought by projection
     1661nbvar=0;
     1662for icell=1:length(CellVarIndex)
     1663    NbDim=NbDimVec(icell);
     1664    if NbDim<3
     1665        continue
     1666    end
     1667    VarIndex=CellVarIndex{icell};%  indices of the selected variables in the list FieldData.ListVarName
     1668    VarType=VarTypeCell{icell};
     1669    ivar_X=VarType.coord_x;
     1670    ivar_Y=VarType.coord_y;
     1671    ivar_Z=VarType.coord_z;
     1672    ivar_U=VarType.vector_x;
     1673    ivar_V=VarType.vector_y;
     1674    ivar_W=VarType.vector_z;
     1675    ivar_C=VarType.scalar ;
     1676    ivar_Anc=VarType.ancillary;
     1677    test_anc=zeros(size(VarIndex));
     1678    test_anc(ivar_Anc)=ones(size(ivar_Anc));
     1679    ivar_F=VarType.warnflag;
     1680    ivar_FF=VarType.errorflag;
     1681    testX=~isempty(ivar_X) && ~isempty(ivar_Y);
     1682    DimCell=FieldData.VarDimName{VarIndex(1)};
     1683    if ischar(DimCell)
     1684        DimCell={DimCell};%name of dimensions
     1685    end
     1686
     1687%% case of input fields with unstructured coordinates
     1688    if testX
     1689        XName=FieldData.ListVarName{ivar_X};
     1690        YName=FieldData.ListVarName{ivar_Y};
     1691        eval(['coord_x=FieldData.' XName ';'])
     1692        eval(['coord_y=FieldData.' YName ';'])
     1693        if length(ivar_Z)==1
     1694            ZName=FieldData.ListVarName{ivar_Z};
     1695            eval(['coord_z=FieldData.' ZName ';'])
     1696        end
     1697
     1698        % translate  initial coordinates
     1699        coord_x=coord_x-ObjectData.Coord(1,1);
     1700        coord_y=coord_y-ObjectData.Coord(1,2);
     1701        if ~isempty(ivar_Z)
     1702            coord_z=coord_z-ObjectData.Coord(1,3);
     1703        end
     1704       
     1705        % selection of the vectors in the projection range
     1706%         if length(ivar_Z)==1 &&  width > 0
     1707%             %components of the unitiy vector normal to the projection plane
     1708%             fieldZ=NormVec_X*coord_x + NormVec_Y*coord_y+ NormVec_Z*coord_z;% distance to the plane           
     1709%             indcut=find(abs(fieldZ) <= width);
     1710%             for ivar=VarIndex
     1711%                 VarName=FieldData.ListVarName{ivar};
     1712%                 eval(['FieldData.' VarName '=FieldData.' VarName '(indcut);']) 
     1713%                     % A VOIR : CAS DE VAR STRUCTUREE MAIS PAS GRILLE REGULIERE : INTERPOLER SUR GRILLE REGULIERE             
     1714%             end
     1715%             coord_x=coord_x(indcut);
     1716%             coord_y=coord_y(indcut);
     1717%             coord_z=coord_z(indcut);
     1718%         end
     1719
     1720       %rotate coordinates if needed: TODO modify
     1721       if testangle
     1722           coord_X=(coord_x *cos(Phi) + coord_y* sin(Phi));
     1723           coord_Y=(-coord_x *sin(Phi) + coord_y *cos(Phi))*cos(Theta);
     1724           if ~isempty(ivar_Z)
     1725               coord_Y=coord_Y+coord_z *sin(Theta);
     1726           end
     1727           
     1728           coord_X=(coord_X *cos(Psi) - coord_Y* sin(Psi));%A VERIFIER
     1729           coord_Y=(coord_X *sin(Psi) + coord_Y* cos(Psi));
     1730           
     1731       else
     1732           coord_X=coord_x;
     1733           coord_Y=coord_y;
     1734           coord_Z=coord_z;
     1735       end
     1736        %restriction to the range of x and y if imposed
     1737        testin=ones(size(coord_X)); %default
     1738        testbound=0;
     1739        if testXMin
     1740            testin=testin & (coord_X >= XMin);
     1741            testbound=1;
     1742        end
     1743        if testXMax
     1744            testin=testin & (coord_X <= XMax);
     1745            testbound=1;
     1746        end
     1747        if testYMin
     1748            testin=testin & (coord_Y >= YMin);
     1749            testbound=1;
     1750        end
     1751        if testYMin
     1752            testin=testin & (coord_Y <= YMax);
     1753            testbound=1;
     1754        end
     1755        if testbound
     1756            indcut=find(testin);
     1757            for ivar=VarIndex
     1758                VarName=FieldData.ListVarName{ivar};
     1759                eval(['FieldData.' VarName '=FieldData.' VarName '(indcut);'])           
     1760            end
     1761            coord_X=coord_X(indcut);
     1762            coord_Y=coord_Y(indcut);
     1763            if length(ivar_Z)==1
     1764                coord_Z=coord_Z(indcut);
     1765            end
     1766        end
     1767        % different cases of projection
     1768        if isequal(ObjectData.ProjMode,'projection')%%%%%%%   NOT USED %%%%%%%%%%
     1769            for ivar=VarIndex %transfer variables to the projection plane
     1770                VarName=FieldData.ListVarName{ivar};
     1771                if ivar==ivar_X %x coordinate
     1772                    eval(['ProjData.' VarName '=coord_X;'])
     1773                elseif ivar==ivar_Y % y coordinate
     1774                    eval(['ProjData.' VarName '=coord_Y;'])
     1775                elseif isempty(ivar_Z) || ivar~=ivar_Z % other variables (except Z coordinate wyhich is not reproduced)
     1776                    eval(['ProjData.' VarName '=FieldData.' VarName ';'])
     1777                end
     1778                if isempty(ivar_Z) || ivar~=ivar_Z
     1779                    ProjData.ListVarName=[ProjData.ListVarName VarName];
     1780                    ProjData.VarDimName=[ProjData.VarDimName DimCell];
     1781                    nbvar=nbvar+1;
     1782                    if isfield(FieldData,'VarAttribute') && length(FieldData.VarAttribute) >=ivar
     1783                        ProjData.VarAttribute{nbvar}=FieldData.VarAttribute{ivar};
     1784                    end
     1785                end
     1786            end 
     1787        elseif isequal(ObjectData.ProjMode,'interp')||isequal(ObjectData.ProjMode,'filter')%interpolate data on a regular grid
     1788            coord_x_proj=XMin:DX:XMax;
     1789            coord_y_proj=YMin:DY:YMax;
     1790            coord_z_proj=ZMin:DZ:ZMax;
     1791            DimCell={'coord_z','coord_y','coord_x'};
     1792            ProjData.ListVarName={'coord_z','coord_y','coord_x'};
     1793            ProjData.VarDimName={'coord_z','coord_y','coord_x'};   
     1794            nbcoord=2; 
     1795            ProjData.coord_z=[ZMin ZMax];
     1796            ProjData.coord_y=[YMin YMax];
     1797            ProjData.coord_x=[XMin XMax];
     1798            if isempty(ivar_X), ivar_X=0; end;
     1799            if isempty(ivar_Y), ivar_Y=0; end;
     1800            if isempty(ivar_Z), ivar_Z=0; end;
     1801            if isempty(ivar_U), ivar_U=0; end;
     1802            if isempty(ivar_V), ivar_V=0; end;
     1803            if isempty(ivar_W), ivar_W=0; end;
     1804            if isempty(ivar_F), ivar_F=0; end;
     1805            if isempty(ivar_FF), ivar_FF=0; end;
     1806            if ~isequal(ivar_FF,0)
     1807                VarName_FF=FieldData.ListVarName{ivar_FF};
     1808                eval(['indsel=find(FieldData.' VarName_FF '==0);'])
     1809                coord_X=coord_X(indsel);
     1810                coord_Y=coord_Y(indsel);
     1811            end
     1812            FF=zeros(1,length(coord_y_proj)*length(coord_x_proj));
     1813            testFF=0;
     1814            [X,Y,Z]=meshgrid(coord_y_proj,coord_z_proj,coord_x_proj);%grid in the new coordinates
     1815            for ivar=VarIndex
     1816                VarName=FieldData.ListVarName{ivar};
     1817                if ~( ivar==ivar_X || ivar==ivar_Y || ivar==ivar_Z || ivar==ivar_F || ivar==ivar_FF || test_anc(ivar)==1)                 
     1818                    ivar_new=ivar_new+1;
     1819                    ProjData.ListVarName=[ProjData.ListVarName {VarName}];
     1820                    ProjData.VarDimName=[ProjData.VarDimName {DimCell}];
     1821                    if isfield(FieldData,'VarAttribute') && length(FieldData.VarAttribute) >=ivar
     1822                        ProjData.VarAttribute{ivar_new+nbcoord}=FieldData.VarAttribute{ivar};
     1823                    end
     1824                    if  ~isequal(ivar_FF,0)
     1825                        eval(['FieldData.' VarName '=FieldData.' VarName '(indsel);'])
     1826                    end
     1827                    eval(['InterpFct=TriScatteredInterp(double(coord_X),double(coord_Y),double(coord_Z),double(FieldData.' VarName '))'])
     1828                    eval(['ProjData.' VarName '=InterpFct(X,Y,Z);'])
     1829%                     eval(['varline=reshape(ProjData.' VarName ',1,length(coord_y_proj)*length(coord_x_proj));'])
     1830%                     FFlag= isnan(varline); %detect undefined values NaN
     1831%                     indnan=find(FFlag);
     1832%                     if~isempty(indnan)
     1833%                         varline(indnan)=zeros(size(indnan));
     1834%                         eval(['ProjData.' VarName '=reshape(varline,length(coord_y_proj),length(coord_x_proj));'])
     1835%                         FF(indnan)=ones(size(indnan));
     1836%                         testFF=1;
     1837%                     end
     1838                    if ivar==ivar_U
     1839                        ivar_U=ivar_new;
     1840                    end
     1841                    if ivar==ivar_V
     1842                        ivar_V=ivar_new;
     1843                    end
     1844                    if ivar==ivar_W
     1845                        ivar_W=ivar_new;
     1846                    end
     1847                end
     1848            end
     1849            if testFF
     1850                ProjData.FF=reshape(FF,length(coord_y_proj),length(coord_x_proj));
     1851                ProjData.ListVarName=[ProjData.ListVarName {'FF'}];
     1852               ProjData.VarDimName=[ProjData.VarDimName {DimCell}];
     1853                ProjData.VarAttribute{ivar_new+1+nbcoord}.Role='errorflag';
     1854            end
     1855        end
     1856       
     1857%% case of input fields defined on a structured  grid
     1858    else
     1859        VarName=FieldData.ListVarName{VarIndex(1)};%get the first variable of the cell to get the input matrix dimensions
     1860        eval(['DimValue=size(FieldData.' VarName ');'])%input matrix dimensions
     1861        DimValue(find(DimValue==1))=[];%remove singleton dimensions       
     1862        NbDim=numel(DimValue);%update number of space dimensions
     1863        nbcolor=1; %default number of 'color' components: third matrix index without corresponding coordinate
     1864        if NbDim>=3
     1865            if NbDim>3
     1866                errormsg='matrices with more than 3 dimensions not handled';
     1867                return
     1868            else
     1869                VarType.coord
     1870                if numel(find(VarType.coord))==2% the third matrix dimension does not correspond to a space coordinate
     1871                    nbcolor=DimValue(3);
     1872                    DimValue(3)=[]; %number of 'color' components updated
     1873                    NbDim=2;% space dimension set to 2
     1874                end
     1875            end
     1876        end
     1877        AYName=FieldData.ListVarName{VarType.coord(NbDim-1)};%name of input x coordinate (name preserved on projection)
     1878        AXName=FieldData.ListVarName{VarType.coord(NbDim)};%name of input y coordinate (name preserved on projection)   
     1879        eval(['AX=FieldData.' AXName ';'])
     1880        eval(['AY=FieldData.' AYName ';'])
     1881        ListDimName=FieldData.VarDimName{VarIndex(1)};
     1882        ProjData.ListVarName=[ProjData.ListVarName {AYName} {AXName}]; %TODO: check if it already exists in Projdata (several cells)
     1883        ProjData.VarDimName=[ProjData.VarDimName {AYName} {AXName}];
     1884
     1885%         for idim=1:length(ListDimName)
     1886%             DimName=ListDimName{idim};
     1887%             if strcmp(DimName,'rgb')||strcmp(DimName,'nb_coord')||strcmp(DimName,'nb_coord_i')
     1888%                nbcolor=DimValue(idim);
     1889%                DimValue(idim)=[];
     1890%             end
     1891%             if isequal(DimName,'nb_coord_j')% NOTE: CASE OF TENSOR NOT TREATED
     1892%                 DimValue(idim)=[];
     1893%             end
     1894%         end 
     1895        Coord_z=[];
     1896        Coord_y=[];
     1897        Coord_x=[];   
     1898
     1899        for idim=1:NbDim %loop on space dimensions
     1900            test_interp(idim)=0;%test for coordiate interpolation (non regular grid), =0 by default
     1901            ivar=VarType.coord(idim);% index of the variable corresponding to the current dimension
     1902            if ~isequal(ivar,0)%  a variable corresponds to the dimension #idim
     1903                eval(['Coord{idim}=FieldData.' FieldData.ListVarName{ivar} ';']) ;% coord values for the input field
     1904                if numel(Coord{idim})==2 %input array defined on a regular grid
     1905                   DCoord_min(idim)=(Coord{idim}(2)-Coord{idim}(1))/DimValue(idim);
     1906                else
     1907                    DCoord=diff(Coord{idim});%array of coordinate derivatives for the input field
     1908                    DCoord_min(idim)=min(DCoord);
     1909                    DCoord_max=max(DCoord);
     1910                %    test_direct(idim)=DCoord_max>0;% =1 for increasing values, 0 otherwise
     1911                    if abs(DCoord_max-DCoord_min(idim))>abs(DCoord_max/1000)
     1912                        msgbox_uvmat('ERROR',['non monotonic dimension variable # ' num2str(idim)  ' in proj_field.m'])
     1913                                return
     1914                    end               
     1915                    test_interp(idim)=(DCoord_max-DCoord_min(idim))> 0.0001*abs(DCoord_max);% test grid regularity
     1916                end
     1917                test_direct(idim)=(DCoord_min(idim)>0);
     1918            else  % no variable associated with the  dimension #idim, the coordinate value is set equal to the matrix index by default
     1919                Coord_i_str=['Coord_' num2str(idim)];
     1920                DCoord_min(idim)=1;%default
     1921                Coord{idim}=[0.5 DimValue(idim)-0.5];
     1922                test_direct(idim)=1;
     1923            end
     1924        end
     1925        if DY==0
     1926            DY=abs(DCoord_min(NbDim-1));
     1927        end
     1928        npY=1+round(abs(Coord{NbDim-1}(end)-Coord{NbDim-1}(1))/DY);%nbre of points after interpol
     1929        if DX==0
     1930            DX=abs(DCoord_min(NbDim));
     1931        end
     1932        npX=1+round(abs(Coord{NbDim}(end)-Coord{NbDim}(1))/DX);%nbre of points after interpol
     1933        for idim=1:NbDim
    13071934            if test_interp(idim)
    13081935                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
    … …  
    14562083                    end
    14572084                    eval(['vec_A=reshape(FieldData.' VarName ',[],nbcolor);'])%put the original image in line             
    1458                     ind_in=find(flagin);
     2085                    %ind_in=find(flagin);
    14592086                    ind_out=find(~flagin);
    14602087                    ICOMB=(XIMA-1)*DimValue(1)+YIMA;
    14612088                    ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A
    1462                     vec_B(ind_in,[1:nbcolor])=vec_A(ICOMB,:);
     2089                    vec_B(flagin,1:nbcolor)=vec_A(ICOMB,:);
    14632090                    for icolor=1:nbcolor
    14642091                        vec_B(ind_out,icolor)=zeros(size(ind_out));
    … …  
    14762103                ProjData.VarAttribute{length(ProjData.ListVarName)}.Role='errorflag';
    14772104            else %3D case
    1478                 if isequal(Theta,0) & isequal(Phi,0)      
     2105                if ~testangle     
    14792106                    % unstructured z coordinate
    14802107                    test_sup=(Coord{1}>=ObjectData.Coord(1,3));
    … …  
    15052132
    15062133    %% projection of  velocity components in the rotated coordinates
    1507     if ~isequal(Phi,0) && length(ivar_U)==1
     2134    if testangle
    15082135        if isempty(ivar_V)
    15092136            msgbox_uvmat('ERROR','v velocity component missing in proj_field.m')
    … …  
    15262153end
    15272154
    1528 %-----------------------------------------------------------------
    1529 %project in a volume
    1530 % A FAIRE ....(copie de proj_plane)
    1531  function  [ProjData,errormsg] = proj_volume(FieldData, ObjectData)
    1532 %-----------------------------------------------------------------
    1533 
    1534 %initialisation of the input parameters of the projection plane
    1535 %-----------------------------------------------------------------
    1536 ProjMode='projection';%direct projection by default
    1537 if isfield(ObjectData,'ProjMode'),ProjMode=ObjectData.ProjMode; end;
    1538 
    1539 %axis origin
    1540 if isempty(ObjectData.Coord)
    1541     ObjectData.Coord(1,1)=0;%origin of the volume set to [0 0] by default
    1542     ObjectData.Coord(1,2)=0;
    1543     ObjectData.Coord(1,3)=0;
    1544 end
    1545 
    1546 %rotation angles
    1547 Phi=0;%default
    1548 Theta=0;
    1549 Psi=0;
    1550 if isfield(ObjectData,'Phi')&& ~isempty(ObjectData.Phi)
    1551     Phi=(pi/180)*ObjectData.Phi;%first Euler angle in radian
    1552 end
    1553 if isfield(ObjectData,'Theta')&& ~isempty(ObjectData.Theta)
    1554     Theta=(pi/180)*ObjectData.Theta;%second Euler angle in radian
    1555 end
    1556 if isfield(ObjectData,'Psi')&& ~isempty(ObjectData.Psi)
    1557     Psi=(pi/180)*ObjectData.Psi;%third Euler angle in radian
    1558 end
    1559 
    1560 %components of the unity vector normal to the projection plane
    1561 NormVec_X=-sin(Phi)*sin(Theta);
    1562 NormVec_Y=cos(Phi)*sin(Theta);
    1563 NormVec_Z=cos(Theta);
    1564 
    1565 % test for 3D fields
    1566 test3D=0;
    1567 if isfield(FieldData,'NbDim')
    1568     test3D=isequal(FieldData.NbDim,3);
    1569 end
    1570 test3C=test3D; %default 3 vel components
    1571 
    1572 %mesh sizes DX and DY
    1573 DX=0;
    1574 DY=0; %default
    1575 if isfield(ObjectData,'DX')&~isempty(ObjectData.DX)
    1576      DX=abs(ObjectData.DX);%mesh of interpolation points
    1577 end
    1578 if isfield(ObjectData,'DY')&~isempty(ObjectData.DY)
    1579      DY=abs(ObjectData.DY);%mesh of interpolation points
    1580 end
    1581 if  ~isequal(ProjMode,'projection') & (DX==0|DY==0)
    1582         errormsg='DX or DY missing';
    1583         display(errormsg)
    1584         return
    1585 end
    1586 if isfield(ObjectData,'DZ')&~isempty(ObjectData.DZ)
    1587      DZ=abs(ObjectData.DZ);%mesh of interpolation points
    1588 end
    1589 
    1590 %extrema along each axis
    1591 testXMin=0;
    1592 testXMax=0;
    1593 testYMin=0;
    1594 testYMax=0;
    1595 if isfield(ObjectData,'RangeX')
    1596         XMin=min(ObjectData.RangeX);
    1597         XMax=max(ObjectData.RangeX);
    1598         testXMin=XMax>XMin;
    1599         testXMax=1;
    1600 end
    1601 if isfield(ObjectData,'RangeY')
    1602         YMin=min(ObjectData.RangeY);
    1603         YMax=max(ObjectData.RangeY);
    1604         testYMin=YMax>YMin;
    1605         testYMax=1;
    1606 end
    1607 width=0;%default width of the projection band
    1608 if isfield(ObjectData,'RangeZ')
    1609       ZMin=min(ObjectData.RangeZ);
    1610         ZMax=max(ObjectData.RangeZ);
    1611         testZMin=YMax>YMin;
    1612         testZMax=1;
    1613 %         width=max(ObjectData.RangeZ);
    1614 end
    1615 
    1616 % initiate Matlab  structure for physical field
    1617 [ProjData,errormsg]=proj_heading(FieldData,ObjectData);
    1618 ProjData.NbDim=3;
    1619 %ProjData.ListDimName={};%name of dimension
    1620 %ProjData.DimValue=[];%values of dimension (nbre of vectors)
    1621 ProjData.ListVarName={};
    1622 ProjData.VarDimName={};
    1623 
    1624 error=0;%default
    1625 flux=0;
    1626 testfalse=0;
    1627 ListIndex={};
    1628 
    1629 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    1630 %group the variables (fields of 'FieldData') in cells of variables with the same dimensions
    1631 %-----------------------------------------------------------------
    1632 idimvar=0;
    1633 [CellVarIndex,NbDim,VarTypeCell,errormsg]=find_field_indices(FieldData);
    1634 if ~isempty(errormsg)
    1635     errormsg=['error in proj_field/proj_plane:' errormsg];
    1636     return
    1637 end
    1638 %LOOP ON GROUPS OF VARIABLES SHARING THE SAME DIMENSIONS
    1639 % CellVarIndex=cells of variable index arrays
    1640 ivar_new=0; % index of the current variable in the projected field
    1641 icoord=0;
    1642 nbcoord=0;%number of added coordinate variables brought by projection
    1643 for icell=1:length(CellVarIndex)
    1644     if NbDim(icell)<2
    1645         continue
    1646     end
    1647     VarIndex=CellVarIndex{icell};%  indices of the selected variables in the list FieldData.ListVarName
    1648     VarType=VarTypeCell{icell};
    1649     ivar_X=VarType.coord_x;
    1650     ivar_Y=VarType.coord_y;
    1651     ivar_Z=VarType.coord_z;
    1652     ivar_U=VarType.vector_x;
    1653     ivar_V=VarType.vector_y;
    1654     ivar_W=VarType.vector_z;
    1655     ivar_C=VarType.scalar ;
    1656     ivar_Anc=VarType.ancillary;
    1657     test_anc=zeros(size(VarIndex));
    1658     test_anc(ivar_Anc)=ones(size(ivar_Anc));
    1659     ivar_F=VarType.warnflag;
    1660     ivar_FF=VarType.errorflag;
    1661     testX=~isempty(ivar_X) && ~isempty(ivar_Y);
    1662     DimCell=FieldData.VarDimName{VarIndex(1)};
    1663     if ischar(DimCell)
    1664         DimCell={DimCell};%name of dimensions
    1665     end
    1666    
    1667 %case of input fields with unstructured coordinates
    1668     if testX
    1669         XName=FieldData.ListVarName{ivar_X};
    1670         YName=FieldData.ListVarName{ivar_Y};
    1671         eval(['coord_x=FieldData.' XName ';'])
    1672         eval(['coord_y=FieldData.' YName ';'])
    1673         if length(ivar_Z)==1
    1674             ZName=FieldData.ListVarName{ivar_Z};
    1675             eval(['coord_z=FieldData.' ZName ';'])
    1676         end
    1677 
    1678         % translate  initial coordinates
    1679         coord_x=coord_x-ObjectData.Coord(1,1);
    1680         coord_y=coord_y-ObjectData.Coord(1,2);
    1681         if ~isempty(ivar_Z)
    1682             coord_z=coord_z-ObjectData.Coord(1,3);
    1683         end
    1684        
    1685         % selection of the vectors in the projection range (3D case)
    1686         if length(ivar_Z)==1 &&  width > 0
    1687             %components of the unitiy vector normal to the projection plane
    1688             fieldZ=NormVec_X*coord_x + NormVec_Y*coord_y+ NormVec_Z*coord_z;% distance to the plane           
    1689             indcut=find(abs(fieldZ) <= width);
    1690             for ivar=VarIndex
    1691                 VarName=FieldData.ListVarName{ivar};
    1692                 eval(['FieldData.' VarName '=FieldData.' VarName '(indcut);'])
    1693 %                 end     
    1694                     % A VOIR : CAS DE VAR STRUCTUREE MAIS PAS GRILLE REGULIERE : INTERPOLER SUR GRILLE REGULIERE             
    1695             end
    1696             coord_x=coord_x(indcut);
    1697             coord_y=coord_y(indcut);
    1698             coord_z=coord_z(indcut);
    1699         end
    1700 
    1701        %rotate coordinates if needed
    1702         if isequal(Phi,0)
    1703             coord_X=coord_x;
    1704             coord_Y=coord_y;
    1705             if ~isequal(Theta,0)
    1706                 coord_Y=coord_Y *cos(Theta);
    1707             end
    1708         else
    1709             coord_X=(coord_x *cos(Phi) + coord_y* sin(Phi));
    1710             coord_Y=(-coord_x *sin(Phi) + coord_y *cos(Phi))*cos(Theta);
    1711         end
    1712         if ~isempty(ivar_Z)
    1713             coord_Y=coord_Y+coord_z *sin(Theta);
    1714         end
    1715         if ~isequal(Psi,0)
    1716                 coord_X=(coord_X *cos(Psi) - coord_Y* sin(Psi));%A VERIFIER
    1717                 coord_Y=(coord_X *sin(Psi) + coord_Y* cos(Psi));
    1718         end
    1719         coord_Z=coord_z; %TO UPDATE
    1720         %restriction to the range of x and y if imposed
    1721         testin=ones(size(coord_X)); %default
    1722         testbound=0;
    1723         if testXMin
    1724             testin=testin & (coord_X >= XMin);
    1725             testbound=1;
    1726         end
    1727         if testXMax
    1728             testin=testin & (coord_X <= XMax);
    1729             testbound=1;
    1730         end
    1731         if testYMin
    1732             testin=testin & (coord_Y >= YMin);
    1733             testbound=1;
    1734         end
    1735         if testYMin
    1736             testin=testin & (coord_Y <= YMax);
    1737             testbound=1;
    1738         end
    1739         if testbound
    1740             indcut=find(testin);
    1741             for ivar=VarIndex
    1742                 VarName=FieldData.ListVarName{ivar};
    1743                 eval(['FieldData.' VarName '=FieldData.' VarName '(indcut);'])           
    1744             end
    1745             coord_X=coord_X(indcut);
    1746             coord_Y=coord_Y(indcut);
    1747             if length(ivar_Z)==1
    1748                 coord_Z=coord_Z(indcut);
    1749             end
    1750         end
    1751         % different cases of projection
    1752         if isequal(ObjectData.ProjMode,'projection')
    1753             %the list of dimension
    1754             %ProjData.ListDimName=[ProjData.ListDimName FieldData.VarDimName(VarIndex(1))];%add the point index to the list of dimensions
    1755             %ProjData.DimValue=[ProjData.DimValue length(coord_X)];
    1756             nbvar=0;
    1757             for ivar=VarIndex %transfer variables to the projection plane
    1758                 VarName=FieldData.ListVarName{ivar};
    1759                 if ivar==ivar_X %x coordinate
    1760                     eval(['ProjData.' VarName '=coord_X;'])
    1761                 elseif ivar==ivar_Y % y coordinate
    1762                     eval(['ProjData.' VarName '=coord_Y;'])
    1763                 elseif isempty(ivar_Z) || ivar~=ivar_Z % other variables (except Z coordinate wyhich is not reproduced)
    1764                     eval(['ProjData.' VarName '=FieldData.' VarName ';'])
    1765                 end
    1766                 if isempty(ivar_Z) || ivar~=ivar_Z
    1767                     ProjData.ListVarName=[ProjData.ListVarName VarName];
    1768                     ProjData.VarDimName=[ProjData.VarDimName DimCell];
    1769                     nbvar=nbvar+1;
    1770                     if isfield(FieldData,'VarAttribute') & length(FieldData.VarAttribute) >=ivar
    1771                         ProjData.VarAttribute{nbvar}=FieldData.VarAttribute{ivar};
    1772                     end
    1773                 end
    1774             end 
    1775         elseif isequal(ObjectData.ProjMode,'interp')||isequal(ObjectData.ProjMode,'filter')%interpolate data on a regular grid
    1776             coord_x_proj=(XMin:DX:XMax);
    1777             coord_y_proj=(YMin:DY:YMax);
    1778             coord_z_proj=(ZMin:DZ:ZMax);
    1779             [XI,YI,ZI]=meshgrid(coord_x_proj,coord_y_proj,coord_z_proj);
    1780             DimCell={'coord_y','coord_x'};
    1781             ProjData.ListVarName={'coord_z','coord_y','coord_x'};
    1782             ProjData.VarDimName={'coord_z','coord_y','coord_x'};   
    1783             nbcoord=3; 
    1784             ProjData.coord_z=[ZMin ZMax];
    1785             ProjData.coord_y=[YMin YMax];
    1786             ProjData.coord_x=[XMin XMax];
    1787             if isempty(ivar_X), ivar_X=0; end;
    1788             if isempty(ivar_Y), ivar_Y=0; end;
    1789             if isempty(ivar_Z), ivar_Z=0; end;
    1790             if isempty(ivar_U), ivar_U=0; end;
    1791             if isempty(ivar_V), ivar_V=0; end;
    1792             if isempty(ivar_W), ivar_W=0; end;
    1793             if isempty(ivar_F), ivar_F=0; end;
    1794             if isempty(ivar_FF), ivar_FF=0; end;
    1795             if ~isequal(ivar_FF,0)
    1796                 VarName_FF=FieldData.ListVarName{ivar_FF};
    1797                 eval(['indsel=find(FieldData.' VarName_FF '==0);'])
    1798                 coord_X=coord_X(indsel);
    1799                 coord_Y=coord_Y(indsel);
    1800             end
    1801             FF=zeros(1,length(coord_y_proj)*length(coord_x_proj));
    1802             testFF=0;
    1803             size(XI)
    1804             size(YI)
    1805             size(ZI)
    1806             for ivar=VarIndex
    1807                 VarName=FieldData.ListVarName{ivar};
    1808                 if ~( ivar==ivar_X || ivar==ivar_Y || ivar==ivar_Z || ivar==ivar_F || ivar==ivar_FF || test_anc(ivar)==1)                 
    1809                     ivar_new=ivar_new+1;
    1810                     ProjData.ListVarName=[ProjData.ListVarName {VarName}];
    1811                     ProjData.VarDimName=[ProjData.VarDimName {DimCell}];
    1812                     if isfield(FieldData,'VarAttribute') && length(FieldData.VarAttribute) >=ivar
    1813                         ProjData.VarAttribute{ivar_new+nbcoord}=FieldData.VarAttribute{ivar};
    1814                     end
    1815                     if  ~isequal(ivar_FF,0)
    1816                         eval(['FieldData.' VarName '=FieldData.' VarName '(indsel);'])
    1817                     end
    1818                     eval(['ProjData.' VarName '=griddata3(double(coord_X),double(coord_Y),double(coord_Z),double(FieldData.' VarName '),XI,YI,ZI);'])
    1819 %                     eval(['varline=reshape(ProjData.' VarName ',1,length(coord_y_proj)*length(coord_x_proj));'])
    1820 %                     FFlag= isnan(varline); %detect undefined values NaN
    1821 %                     indnan=find(FFlag);
    1822 %                     if~isempty(indnan)
    1823 %                         varline(indnan)=zeros(size(indnan));
    1824 %                         eval(['ProjData.' VarName '=reshape(varline,length(coord_y_proj),length(coord_x_proj));'])
    1825 %                         FF(indnan)=ones(size(indnan));
    1826 %                         testFF=1;
    1827 %                     end
    1828 %                     if ivar==ivar_U
    1829 %                         ivar_U=ivar_new;
    1830 %                     end
    1831 %                     if ivar==ivar_V
    1832 %                         ivar_V=ivar_new;
    1833 %                     end
    1834 %                     if ivar==ivar_W
    1835 %                         ivar_W=ivar_new;
    1836 %                     end
    1837                 end
    1838             end
    1839             if testFF
    1840                 ProjData.FF=reshape(FF,length(coord_y_proj),length(coord_x_proj));
    1841                 ProjData.ListVarName=[ProjData.ListVarName {'FF'}];
    1842                ProjData.VarDimName=[ProjData.VarDimName {DimCell}];
    1843                 ProjData.VarAttribute{ivar_new+1+nbcoord}.Role='errorflag';
    1844             end
    1845         end
    1846 %case of fields defined on a structured  grid
    1847     else 
    1848         AYName=FieldData.ListVarName{VarType.coord(1)};
    1849         AXName=FieldData.ListVarName{VarType.coord(2)};
    1850         eval(['AX=FieldData.' AXName ';'])
    1851         eval(['AY=FieldData.' AYName ';'])
    1852         VarName=FieldData.ListVarName{VarIndex(1)};
    1853         eval(['DimValue=size(FieldData.' VarName ');'])
    1854         ListDimName=FieldData.VarDimName{VarIndex(1)};
    1855         ProjData.ListVarName=[{AYName} {AXName} ProjData.ListVarName]; %TODO: check if it already exists in Projdata (several cells)
    1856         ProjData.VarDimName=[{AYName} {AXName} ProjData.VarDimName];
    1857         nbcolor=1; %default
    1858         for idim=1:length(ListDimName)
    1859             DimName=ListDimName{idim};
    1860             if isequal(DimName,'rgb')|isequal(DimName,'nb_coord')|isequal(DimName,'nb_coord_i')
    1861                nbcolor=DimValue(idim);
    1862                DimIndices(idim)=[];
    1863                DimValue(idim)=[];
    1864             end
    1865             if isequal(DimName,'nb_coord_j')% NOTE: CASE OF TENSOR NOT TREATED
    1866                 DimIndices(idim)=[];
    1867                 DimValue(idim)=[];
    1868             end
    1869         end 
    1870         ind_1=find(DimValue==1);
    1871         DimIndices(ind_1)=[]; %suppress singleton dimensions
    1872 %         indxy=find(DimVarIndex(DimIndices));%select dimension variables (DimIndices non zero)
    1873         NbDim=length(DimIndices);%number of space dimensions
    1874         Coord_z=[];
    1875         Coord_y=[];
    1876         Coord_x=[];   
    1877    
    1878         for idim=1:NbDim %loop on space dimensions
    1879             test_interp(idim)=0;%test for coordiate interpolation (non regular grid), =0 by default
    1880             ivar=DimVarIndex(DimIndices(idim));% index of the variable corresponding to the current dimension
    1881             if ~isequal(ivar,0)%  a variable corresponds to the current dimension
    1882                 eval(['Coord{idim}=FieldData.' FieldData.ListVarName{ivar} ';']) ;% position for the first index
    1883                 DCoord=diff(Coord{idim});
    1884                 DCoord_min(idim)=min(DCoord);
    1885                 DCoord_max=max(DCoord);
    1886                 test_direct(idim)=DCoord_max>0;% =1 for increasing values, 0 otherwise
    1887                 test_direct_min=DCoord_min(idim)>0;% =1 for increasing values, 0 otherwise
    1888                 if ~isequal(test_direct(idim),test_direct_min)
    1889                      msgbox_uvmat('ERROR',['non monotonic dimension variable # ' num2str(idim)  ' in proj_field.m'])
    1890                                 return
    1891                 end               
    1892                 test_interp(idim)=(DCoord_max-DCoord_min(idim))> 0.0001*abs(DCoord_max);% test grid regularity
    1893             else  % no variable associated with the first dimension, look for variable  attributes Coord_1, _2 or _3
    1894                 Coord_i_str=['Coord_' num2str(idim)];
    1895                 DCoord_min(idim)=1;%default
    1896                 Coord{idim}=[0.5 DimValue(idim)-0.5];
    1897                 test_direct(idim)=1;
    1898             end
    1899         end
    1900         if NbDim==2
    1901             if DY==0
    1902                 DY=abs(DCoord_min(1));
    1903             end
    1904             npY=1+round(abs(Coord{1}(end)-Coord{1}(1))/DY);%nbre of points after interpolation
    1905             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
    1906             if DX==0
    1907                 DX=abs(DCoord_min(2));
    1908             end
    1909             npX=1+round(abs(Coord{2}(end)-Coord{2}(1))/DX);%nbre of points after interpol 
    1910             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
    1911             Coord_y=linspace(Coord{1}(1),Coord{1}(end),npY);
    1912             test_direct_y=test_direct(1);
    1913             Coord_x=linspace(Coord{2}(1),Coord{2}(end),npX);
    1914             test_direct_x=test_direct(2);
    1915             DAX=DCoord_min(2);
    1916             DAY=DCoord_min(1);
    1917         elseif NbDim==3
    1918             DZ=abs(DCoord_min(1));
    1919             npz=1+round(abs(Coord{1}(end)-Coord{1}(1))/DZ);%nbre of points after interpolation
    1920             if DY==0
    1921                 DY=abs(DCoord_min(2));
    1922             end
    1923             npY=1+round(abs(Coord{2}(end)-Coord{2}(1))/DY);%nbre of points after interpol
    1924             DimValue(1)=1+round(abs(Coord{2}(end)-Coord{2}(1))/abs(DCoord_min(2)));%nbre of points before interpol
    1925             if DX==0
    1926                 DX=abs(DCoord_min(3));
    1927             end
    1928             npX=1+round(abs(Coord{3}(end)-Coord{3}(1))/DX);%nbre of points after interpol
    1929             DimValue(2)=1+round(abs(Coord{3}(end)-Coord{3}(1))/abs(DCoord_min(3)));%nbre of points before interpol
    1930             Coord_z=linspace(Coord{1}(1),Coord{1}(end),npz);
    1931             test_direct_z=test_direct(1);
    1932             Coord_y=linspace(Coord{2}(1),Coord{2}(end),npY);
    1933             test_direct_y=test_direct(2);
    1934             Coord_x=linspace(Coord{3}(1),Coord{3}(end),npX);
    1935             test_direct_x=test_direct(3);
    1936         end 
    1937         minAX=min(Coord_x);
    1938         maxAX=max(Coord_x);
    1939         minAY=min(Coord_y);
    1940         maxAY=max(Coord_y);
    1941         xcorner=[minAX maxAX minAX maxAX]-ObjectData.Coord(1,1);
    1942         ycorner=[maxAY maxAY minAY minAY]-ObjectData.Coord(1,2);
    1943         xcor_new=xcorner*cos(Phi)+ycorner*sin(Phi);%coord new frame
    1944         ycor_new=-xcorner*sin(Phi)+ycorner*cos(Phi);
    1945         if ~testXMax
    1946             XMax=max(xcor_new);
    1947         end
    1948         if ~testXMin
    1949             XMin=min(xcor_new);
    1950         end
    1951         if ~testYMax
    1952             YMax=max(ycor_new);
    1953         end
    1954         if ~testYMin
    1955             YMin=min(ycor_new);
    1956         end
    1957         DXinit=(maxAX-minAX)/(DimValue(2)-1);
    1958         DYinit=(maxAY-minAY)/(DimValue(1)-1);
    1959         if DX==0
    1960             DX=DXinit;
    1961         end
    1962         if DY==0
    1963             DY=DYinit;
    1964         end
    1965         npX=floor((XMax-XMin)/DX+1);
    1966         npY=floor((YMax-YMin)/DY+1);   
    1967         if test_direct_y
    1968             coord_y_proj=linspace(YMin,YMax,npY);%abscissa of the new pixels along the line
    1969         else
    1970             coord_y_proj=linspace(YMax,YMin,npY);%abscissa of the new pixels along the line
    1971         end
    1972         if test_direct_x
    1973             coord_x_proj=linspace(XMin,XMax,npX);%abscissa of the new pixels along the line
    1974         else
    1975             coord_x_proj=linspace(XMax,XMin,npX);%abscissa of the new pixels along the line
    1976         end
    1977        
    1978         % case with no rotation and interpolation
    1979         if isequal(ProjMode,'projection') && isequal(Phi,0) && isequal(Theta,0) && isequal(Psi,0)
    1980             if test_direct(1)
    1981                 min_indy=ceil((YMin-Coord{1}(1))/DYinit)+1;
    1982                 max_indy=floor((YMax-Coord{1}(1))/DYinit)+1;
    1983                 Ybound(1)=Coord{1}(1)+DYinit*(min_indy-1);
    1984                 Ybound(2)=Coord{1}(1)+DYinit*(max_indy-1);
    1985             else
    1986                 min_indy=ceil((Coord{1}(1)-YMax)/DYinit)+1;
    1987                 max_indy=floor((Coord{1}(1)-YMin)/DYinit)+1;
    1988                 Ybound(2)=Coord{1}(1)-DYinit*(max_indy-1);
    1989                 Ybound(1)=Coord{1}(1)-DYinit*(min_indy-1);
    1990             end             
    1991             if test_direct(2)==1
    1992                 min_indx=ceil((XMin-Coord{2}(1))/DXinit)+1;
    1993                 max_indx=floor((XMax-Coord{2}(1))/DXinit)+1;
    1994                 Xbound(1)=Coord{2}(1)+DXinit*(min_indx-1);
    1995                 Xbound(2)=Coord{2}(1)+DXinit*(max_indx-1);
    1996             else
    1997                 min_indx=ceil((Coord{2}(1)-XMax)/DXinit)+1;
    1998                 max_indx=floor((Coord{2}(1)-XMin)/DXinit)+1;
    1999                 Xbound(2)=Coord{2}(1)+DXinit*(max_indx-1);
    2000                 Xbound(1)=Coord{2}(1)+DXinit*(min_indx-1);
    2001             end
    2002             min_indy=max(min_indy,1);% deals with margin (bound lower than the first index)
    2003             min_indx=max(min_indx,1);
    2004             max_indy=min(max_indy,DimValue(1));
    2005             max_indx=min(max_indx,DimValue(2));
    2006             for ivar=VarIndex
    2007                 VarName=FieldData.ListVarName{ivar};
    2008                 ProjData.ListVarName=[ProjData.ListVarName VarName];
    2009                 %ProjData.VarDimIndex=[ProjData.VarDimIndex [NbDim-1 NbDim]];
    2010                 if length(FieldData.VarAttribute)>=ivar
    2011                     ProjData.VarAttribute{length(ProjData.ListVarName)}=FieldData.VarAttribute{ivar};
    2012                 end
    2013                 eval(['ProjData.' VarName '=FieldData.' VarName '(min_indy:max_indy,min_indx:max_indx) ;']);
    2014             end         
    2015         else
    2016         % case with rotation and/or interpolation
    2017             if isempty(Coord_z) %2D case
    2018                 [X,Y]=meshgrid(coord_x_proj,coord_y_proj);%grid in the new coordinates
    2019                 XIMA=ObjectData.Coord(1,1)+(X)*cos(Phi)-Y*sin(Phi);%corresponding coordinates in the original image
    2020                 YIMA=ObjectData.Coord(1,2)+(X)*sin(Phi)+Y*cos(Phi);
    2021                 XIMA=(XIMA-minAX)/DXinit+1;% image index along x
    2022                 YIMA=(-YIMA+maxAY)/DYinit+1;% image index along y
    2023                 XIMA=reshape(round(XIMA),1,npX*npY);%indices reorganized in 'line'
    2024                 YIMA=reshape(round(YIMA),1,npX*npY);
    2025                 flagin=XIMA>=1 & XIMA<=DimValue(2) & YIMA >=1 & YIMA<=DimValue(1);%flagin=1 inside the original image
    2026                 if isequal(ObjectData.ProjMode,'filter')
    2027                     npx_filter=ceil(abs(DX/DAX));
    2028                     npy_filter=ceil(abs(DY/DAY));
    2029                     Mfilter=ones(npy_filter,npx_filter)/(npx_filter*npy_filter);
    2030                     test_filter=1;
    2031                 else
    2032                     test_filter=0;
    2033                 end
    2034                 for ivar=VarIndex
    2035                     VarName=FieldData.ListVarName{ivar} ;
    2036                     if test_interp(1) | test_interp(2)%interpolate on a regular grid         
    2037                           eval(['FieldData.' VarName '=interp2(Coord{2},Coord{1},FieldData.' VarName ',Coord_x,Coord_y'');']) %TO TEST
    2038                     end
    2039                     %filter the field (image) if option 'filter' is used
    2040                     if test_filter 
    2041                          Aclass=class(FieldData.A);
    2042                          eval(['FieldData.' VarName '=filter2(Mfilter,FieldData.' VarName ',''valid'');'])
    2043                          if ~isequal(Aclass,'double')
    2044                              eval(['FieldData.' VarName '=' Aclass '(FieldData.' VarName ');'])%revert to integer values
    2045                          end
    2046                     end
    2047                     eval(['vec_A=reshape(FieldData.' VarName ',[],nbcolor);'])%put the original image in line             
    2048                     ind_in=find(flagin);
    2049                     ind_out=find(~flagin);
    2050                     ICOMB=(XIMA-1)*DimValue(1)+YIMA;
    2051                     ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A
    2052                     vec_B(ind_in,1:nbcolor)=vec_A(ICOMB,:);
    2053                     for icolor=1:nbcolor
    2054                         vec_B(ind_out,icolor)=zeros(size(ind_out));
    2055                     end
    2056                     ProjData.ListVarName=[ProjData.ListVarName VarName];                 
    2057                     if length(FieldData.VarAttribute)>=ivar
    2058                         ProjData.VarAttribute{length(ProjData.ListVarName)+nbcoord}=FieldData.VarAttribute{ivar};
    2059                     end     
    2060                     eval(['ProjData.' VarName '=reshape(vec_B,npY,npX,nbcolor);']);
    2061                 end
    2062                 ProjData.FF=reshape(~flagin,npY,npX);%false flag A FAIRE: tenir compte d'un flga antérieur 
    2063                 ProjData.ListVarName=[ProjData.ListVarName 'FF'];
    2064                 ProjData.VarAttribute{length(ProjData.ListVarName)}.Role='errorflag';
    2065             else %3D case
    2066                 if isequal(Theta,0) & isequal(Phi,0)       
    2067                     test_sup=(Coord{1}>=ObjectData.Coord(1,3));
    2068                     iz_sup=find(test_sup);
    2069                     iz=iz_sup(1);
    2070                     if iz>=1 & iz<=npz
    2071                         %ProjData.ListDimName=[ProjData.ListDimName ListDimName(2:end)];
    2072                         %ProjData.DimValue=[ProjData.DimValue npY npX];
    2073                         for ivar=VarIndex
    2074                             VarName=FieldData.ListVarName{ivar};
    2075                             ProjData.ListVarName=[ProjData.ListVarName VarName];
    2076                             ProjData.VarAttribute{length(ProjData.ListVarName)}=FieldData.VarAttribute{ivar}; %reproduce the variable attributes 
    2077                             eval(['ProjData.' VarName '=squeeze(FieldData.' VarName '(iz,:,:));'])% select the z index iz
    2078                             %TODO : do a vertical average for a thick plane
    2079                             if test_interp(2) || test_interp(3)
    2080                                 eval(['ProjData.' VarName '=interp2(Coord{3},Coord{2},ProjData.' VarName ',Coord_x,Coord_y'');'])
    2081                             end
    2082                         end
    2083                     end
    2084                 else
    2085                     errormsg='projection of structured coordinates on oblique plane not yet implemented';
    2086                     %TODO: use interp3
    2087                     return
    2088                 end
    2089             end
    2090         end
    2091     end
    2092     %projection of  velocity components in the rotated coordinates
    2093     if ~isequal(Phi,0) && length(ivar_U)==1
    2094         if isempty(ivar_V)
    2095             msgbox_uvmat('ERROR','v velocity component missing in proj_field.m')
    2096             return
    2097         end
    2098         UName=FieldData.ListVarName{ivar_U};
    2099         VName=FieldData.ListVarName{ivar_V};   
    2100         eval(['ProjData.' UName  '=cos(Phi)*ProjData.' UName '+ sin(Phi)*ProjData.' VName ';'])
    2101         eval(['ProjData.' VName  '=cos(Theta)*(-sin(Phi)*ProjData.' UName '+ cos(Phi)*ProjData.' VName ');'])
    2102         if ~isempty(ivar_W)
    2103             WName=FieldData.ListVarName{ivar_W};
    2104             eval(['ProjData.' VName '=ProjData.' VName '+ ProjData.' WName '*sin(Theta);'])%
    2105             eval(['ProjData.' WName '=NormVec_X*ProjData.' UName '+ NormVec_Y*ProjData.' VName '+ NormVec_Z* ProjData.' WName ';']);
    2106         end
    2107         if ~isequal(Psi,0)
    2108             eval(['ProjData.' UName '=cos(Psi)* ProjData.' UName '- sin(Psi)*ProjData.' VName ';']);
    2109             eval(['ProjData.' VName '=sin(Psi)* ProjData.' UName '+ cos(Psi)*ProjData.' VName ';']);
    2110         end
    2111     end
    2112 end
    2113 
    2114 %-----------------------------------------------------------------
    2115 %transmit the global attributes
     2155%------------------------------------------------------------------------
     2156%--- transfer the global attributes
    21162157function [ProjData,errormsg]=proj_heading(FieldData,ObjectData)
    2117 %-----------------------------------------------------------------
    2118 % ProjData=FieldData;
     2158%------------------------------------------------------------------------
    21192159ProjData=[];%default
    2120 errormsg=[];%default
     2160errormsg='';%default
     2161
     2162%% transfer error
     2163if isfield(FieldData,'Txt')
     2164    errormsg=FieldData.Txt; %transmit erreur message
     2165    return;
     2166end
     2167
     2168%% transfer global attributes
    21212169if ~isfield(FieldData,'ListGlobalAttribute')
    21222170    ProjData.ListGlobalAttribute={};
    … …  
    21242172    ProjData.ListGlobalAttribute=FieldData.ListGlobalAttribute;
    21252173end
    2126 if isfield(FieldData,'Txt')
    2127     errormsg=FieldData.Txt; %transmit erreur message
    2128     return;
    2129 end
    21302174for iattr=1:length(ProjData.ListGlobalAttribute)
    21312175    AttrName=ProjData.ListGlobalAttribute{iattr};
    … …  
    21342178    end
    21352179end
     2180
     2181%% transfer coordinate unit
    21362182if isfield(FieldData,'CoordUnit')
    21372183    if isfield(ObjectData,'CoordUnit')&~isequal(FieldData.CoordUnit,ObjectData.CoordUnit)
    … …  
    21432189end
    21442190
     2191%% store the properties of the projection object
    21452192ListObject={'Style','ProjMode','RangeX','RangeY','RangeZ','Phi','Theta','Psi','Coord'};
    21462193for ilist=1:length(ListObject)

  • trunk/src/read_civxdata.m

    r186 r206  
    4848
    4949function [Field,VelTypeOut]=read_civxdata(filename,FieldNames,VelType)
    50 
     50%% default input
    5151if ~exist('VelType','var')
    5252    VelType=[];
    … …  
    5858    FieldNames=[]; %default
    5959end
     60
     61%% reading data
    6062VelTypeOut=VelType;%default
    6163[var,role,units,vel_type_out_cell]=varcivx_generator(FieldNames,VelType);%determine the names of constants and variables to read
    62 [Field,vardetect,ichoice]=nc2struct(filename,var);
    63 % if isfield(Field,'Txt')
    64 %     return % error in file reading
    65 % end
     64[Field,vardetect,ichoice]=nc2struct(filename,var);%read the variables in the netcdf file
     65if isfield(Field,'Txt')
     66    return
     67end
    6668if vardetect(1)==0
    6769     Field.Txt=[ 'requested field not available in ' filename '/' VelType];
     70     return
    6871end
    6972var_ind=find(vardetect);
    7073for ivar=1:length(var_ind)
    7174    Field.VarAttribute{ivar}.Role=role{var_ind(ivar)};
    72      Field.VarAttribute{ivar}.Mesh=0.1;%typical mesh for histograms O.1 pixel
    73 %     Field.VarAttribute{ivar}.units=units{var_ind(ivar)};% not necessary: set with calc_field
     75    Field.VarAttribute{ivar}.Mesh=0.1;%typical mesh for histograms O.1 pixel
    7476end
    7577VelTypeOut=VelType;
    … …  
    7880end
    7981
    80 %adjust for Djui:
     82%% adjust for Djui:
    8183if isfield(Field,'DjUi')
    8284    Field.ListVarName(end-2:end)=[];
    … …  
    8688end
    8789
    88 %renaiming for standard conventions
     90%% renaming for standard conventions
    8991Field.NbCoord=Field.nb_coord;
    9092Field.NbDim=Field.nb_dim;
    9193
    92 % CivStage
     94%% CivStage
    9395if isfield(Field,'patch2')&& isequal(Field.patch2,1)
    9496    Field.CivStage=6;
    … …  
    178180function [var,role,units,vel_type_out]=varcivx_generator(FieldNames,vel_type)
    179181
    180 %default input values
     182%% default input values
    181183if ~exist('vel_type','var'),vel_type=[];end;
    182184if iscell(vel_type),vel_type=vel_type{1}; end;%transform cell to string if needed
    183 % if ~exist('display','var'),display=[];end;
    184185if ~exist('FieldNames','var'),FieldNames={'ima_cor'};end;%default scalar
    185186if ischar(FieldNames), FieldNames={FieldNames}; end;
    186187
    187 %select the priority order for automatic vel_type selection
     188%% select the priority order for automatic vel_type selection
    188189testder=0;
    189190for ilist=1:length(FieldNames)
    … …  
    221222vel_type_out=vel_type_out';
    222223
    223 %determine names of netcdf variables to read
     224%% determine names of netcdf variables to read
    224225var={'X','Y','Z','U','V','W','C','F','FF'};
    225226role={'coord_x','coord_y','coord_z','vector_x','vector_y','vector_z','ancillary','warnflag','errorflag'};
    … …  
    234235end
    235236
    236  
    237 %-------------------------
    238 %determine  var names to read
    239 %--------------------------------------
     237%------------------------------------------------------------------------ 
     238%--- determine  var names to read
    240239function varin=varname1(vel_type,FieldNames)
    241 
     240%------------------------------------------------------------------------
    242241testder=0;
    243242C1='';

  • trunk/src/read_get_field.m

    r179 r206  
    388388           
    389389    if test3D %  (a revoir) 
    390          %scalar w variable
     390        %scalar w variable
    391391        inputlist=get(handles.vector_z,'String');
    392392        val=get(handles.vector_z,'Value');%selected indices in the ordinate listbox
    393393        VarNameW=inputlist{val}; %name of the variable in the list
    394         VarIndex=name2index(VarNameW,Field.ListVarName);%index of the variable in ListVarName
    395          %check consistency of dimensions with u
    396         dimname_w=Field.VarDimName{VarIndex};
    397         if ~isequal(dimname_w,dimname_u)
    398            errormsg='inconsistent dimensions for u and v';
    399             return
    400         end
    401         nbvar=nbvar+1;
    402 %         w_index=nbvar;
    403         ListVarName{nbvar}=Field.ListVarName{VarIndex};
    404         VarDimName{nbvar}=dimname_u;
    405         if numel(VarAttribute)>=VarIndex
    406             SubVarAttribute{nbvar}=VarAttribute{VarIndex};
    407         end
    408         SubVarAttribute{nbvar}.Role='vector_z';
     394        VarIndex=name2index(VarNameW,Field.ListVarName);%index of the variable in ListVarName
     395        %check consistency of dimensions with u
     396        if ~isempty( VarIndex)
     397            dimname_w=Field.VarDimName{VarIndex};
     398            if ~isequal(dimname_w,dimname_u)
     399                errormsg='inconsistent dimensions for u and w';
     400                return
     401            end
     402            nbvar=nbvar+1;
     403            %         w_index=nbvar;
     404            ListVarName{nbvar}=Field.ListVarName{VarIndex};
     405            VarDimName{nbvar}=dimname_u;
     406            if numel(VarAttribute)>=VarIndex
     407                SubVarAttribute{nbvar}=VarAttribute{VarIndex};
     408            end
     409            SubVarAttribute{nbvar}.Role='vector_z';
     410        end
    409411    end 
    410412   
    … …  
    446448
    447449%permute indices if coord_y is not the first matrix index: scalar case
     450NbDim=2; %default
    448451if test_scalar
    449452    VarNameA=Field.ListVarName{VarIndexA};%name of the scalar variable
    … …  
    455458        DimCellA=DimCellA(end-numel(npxy)+1:end); %suppress the first singletons) dimensions
    456459    end
    457     %ind_single=find(npxy==1);
    458460    ind_select=find(npxy~=1);%look for non singleton dimensions
    459461    DimCellA=DimCellA(ind_select);%dimension names for the scalar variable, after removing singletons
    … …  
    461463    dimA=[];
    462464    if test_zdimvar%dim_x && dim_y && ~isempty(VarSubIndexA)
     465        NbDim=3;% field considered as 3D if a z coordinate is defined (to distinguish for instance from 2D color images with 3 components)
    463466        ind_singleton=find(strcmp(dimname_z,SingleCellA),1);% look for coincidence of dimension with one of the singleton dimensions
    464467        if ~isempty(ind_singleton)
    … …  
    468471        icoord=find(strcmp(dimname_z,DimCellA),1);% a dimension variable
    469472        dimA=[dimA icoord];
    470 %         for icoord=1:numel(DimCellA)% look for coincidence of dimension with one of the dimensions of the scalar
    471 %              if strcmp(dimname_z,DimCellA{icoord})% a dimension variable
    472 %                  dimA=[dimA icoord];
    473 %                  break
    474 %              end
    475 %         end
    476473    end
    477474    if test_ydimvar%dim_x && dim_y && ~isempty(VarSubIndexA)
    … …  
    507504        end
    508505        ind_select=find(npxy~=1) ;%look for non singleton dimensions
    509         DimCell=DimCell(ind_select);
     506        DimCell=DimCell(ind_select);%list of dimension names for the scalar, after singleton removal
    510507        npxy=npxy(ind_select);
    511508        testold=0;
    … …  
    519516        end
    520517        if empty_coord_x       
    521                 coord_x_name=DimCell{2};
     518                coord_x_name=DimCell{NbDim};
    522519                SubField.ListVarName=[{coord_x_name} SubField.ListVarName];
    523520                SubField.VarDimName=[{coord_x_name} SubField.VarDimName]; 
    … …  
    525522                    eval(['SubField.' coord_x_name '=linspace(Coord_2(1),Coord_2(end),npxy(2));'])
    526523                else
    527                     eval(['SubField.' coord_x_name '=[0.5 npxy(2)-0.5];'])
     524                    eval(['SubField.' coord_x_name '=[0.5 npxy(NbDim)-0.5];'])
    528525                end
    529526           
    … …  
    536533        end
    537534        if empty_coord_y
    538             coord_y_name=DimCell{1};
     535            coord_y_name=DimCell{NbDim-1};
    539536            SubField.ListVarName=[{coord_y_name} SubField.ListVarName];
    540537            SubField.VarDimName=[{coord_y_name} SubField.VarDimName];
    … …  
    542539                eval(['SubField.' coord_y_name '=linspace(Coord_1(1),Coord_1(end),npxy(1));'])
    543540            else
    544                 eval(['SubField.' coord_y_name '=[npxy(1)-0.5 0.5];'])
     541                eval(['SubField.' coord_y_name '=[npxy(NbDim-1)-0.5 0.5];'])
    545542            end
    546543            if ~testold

  • trunk/src/read_set_object.m

    r158 r206  
    2222        value=get(handles.ProjMode,'Value');
    2323        data.ProjMode=menu{value};
    24 %       menu=get(handles.CoordUnit,'String');
    25 %       value=get(handles.MenuCoord,'Value');
    2624        data.CoordUnit=get(handles.CoordUnit,'String');
    2725    testcalib=0;
    … …  
    3432if ~testcalib
    3533    if isequal(get(handles.Phi,'Visible'),'on')
    36         data.Phi=str2num(get(handles.Phi,'String'));
     34        data.Angle(1)=str2double(get(handles.Phi,'String'));
    3735    end
    3836    if isequal(get(handles.Theta,'Visible'),'on')
    39         data.Theta=str2num(get(handles.Theta,'String'));
     37        data.Angle(2)=str2double(get(handles.Theta,'String'));
    4038    end
    4139    if isequal(get(handles.Psi,'Visible'),'on')
    42         data.Psi=str2num(get(handles.Psi,'String'));
     40        data.Angle(3)=str2double(get(handles.Psi,'String'));
    4341    end
    4442    if isequal(get(handles.DX,'Visible'),'on')

  • trunk/src/series.m

    r205 r206  
    248248        'Pick a file',oldfile);
    249249fileinput=[PathName FileName];%complete file name
    250 testblank=findstr(fileinput,' ');%look for blanks
    251 if ~isempty(testblank)
    252     errordlg('forbidden input file name: contain blanks')
    253     return
    254 end
     250%testblank=findstr(fileinput,' ');%look for blanks
     251% if ~isempty(testblank)
     252%     errordlg('forbidden input file name: contain blanks')
     253%     return
     254% end
    255255sizf=size(fileinput);
    256256if (~ischar(fileinput)|~isequal(sizf(1),1)),return;end
    … …  
    327327RootFileCell=get(handles.RootFile,'String');
    328328oldfile=''; %default
    329 if isempty(RootPathCell)|isequal(RootPathCell,{''})%loads the previously stored file name and set it as default in the file_input box
     329if isempty(RootPathCell)||isequal(RootPathCell,{''})%loads the previously stored file name and set it as default in the file_input box
    330330     dir_perso=prefdir;
    331331     profil_perso=fullfile(dir_perso,'uvmat_perso.mat');
    … …  
    352352        'Pick a file',oldfile);
    353353fileinput=[PathName FileName];%complete file name
    354 testblank=findstr(fileinput,' ');%look for blanks
    355 if ~isempty(testblank)
    356     errordlg('forbidden input file name: contain blanks')
    357     return
    358 end
     354% testblank=findstr(fileinput,' ');%look for blanks
     355% if ~isempty(testblank)
     356%     errordlg('forbidden input file name: contain blanks')
     357%     return
     358% end
    359359sizf=size(fileinput);
    360360if (~ischar(fileinput)|~isequal(sizf(1),1)),return;end
    … …  
    362362SeriesData=[];%dfault
    363363if isequal(ext,'.xml')
    364     errordlg('input file type not implemented')%A Faire: ouvrir le fichier pour naviguer
     364    msgbox_uvmat('ERROR','input file type not implemented')%A Faire: ouvrir le fichier pour naviguer
    365365elseif isequal(ext,'.xls')
    366     errordlg('input file type not implemented')%A Faire: ouvrir le fichier pour naviguer
     366    msgbox_uvmat('ERROR','input file type not implemented')%A Faire: ouvrir le fichier pour naviguer
    367367else
    368368    update_file(hObject, eventdata, handles,fileinput,1)
    … …  
    16221622       series_fct{ilist-nb_builtin_ACTION}=fullfile(list_path{ilist},[menu_str{ilist} '.m']);     
    16231623   end
    1624    if nb_builtin_ACTION+1>=length(menu_str)-1
    1625    if exist(profil_perso,'file') && nb_builtin_ACTION+1>=length(menu_str)-1
    1626         save(profil_perso,'series_fct','-append')
    1627    else
    1628     txt=ver('MATLAB');
    1629     Release=txt.Release;
    1630         relnumb=str2num(Release(3:4));
    1631         if relnumb >= 14
    1632             save(profil_perso,'series_fct','-V6')
    1633         else
    1634             save(profil_perso, 'series_fct')
    1635         end
    1636    end
     1624   if nb_builtin_ACTION+1<=length(menu_str)-1
     1625       if exist(profil_perso,'file')% && nb_builtin_ACTION+1>=length(menu_str)-1
     1626           save(profil_perso,'series_fct','-append')
     1627       else
     1628           txt=ver('MATLAB');
     1629           Release=txt.Release;
     1630           relnumb=str2num(Release(3:4));
     1631           if relnumb >= 14%recent relaese of Matlab
     1632               save(profil_perso,'series_fct','-V6')
     1633           else
     1634               save(profil_perso, 'series_fct')
     1635           end
     1636       end
    16371637   end
    16381638end

  • trunk/src/set_object.m

    r204 r206  
    8080
    8181%default
    82 if ~exist('ZBound','var')
    83     ZBound=0; %default
     82if ~exist('ZBounds','var')
     83    ZBounds=0; %default
    8484end
    8585set(hObject,'KeyPressFcn',{'keyboard_callback',handles})%set keyboard action function (allow action on uvmat when set_object is in front)
    … …  
    177177        set(handles.ZMax,'String',num2str(max(data.RangeZ),3))
    178178        DZ=max(data.RangeZ);%slider step
    179         if ZBounds(2)~=ZBounds(1)
     179        if ~isnan(ZBounds(1)) && ZBounds(2)~=ZBounds(1)
    180180            rel_step(1)=min(DZ/(ZBounds(2)-ZBounds(1)),0.2);%must be smaller than 1
    181181            rel_step(2)=0.1;
    … …  
    422422            set(handles.DY,'Visible','off')
    423423        end
    424         if isequal(ObjectStyle,'volume') && isequal(ProjMode,'interp')
     424        if isequal(ProjMode,'interp')
    425425            set(handles.DZ,'Visible','on') 
    426426        end
    … …  
    433433        set(handles.XObject,'TooltipString',['XObject:  x coordinate of the axis origin for the ' ObjectStyle])
    434434        set(handles.YObject,'TooltipString',['YObject:  y coordinate of the axis origin for the ' ObjectStyle])
    435         if test3D
     435%         if test3D
    436436            set(handles.Theta,'Visible','on')
    437437            set(handles.Psi,'Visible','on')
    438438            set(handles.ZMin,'Visible','on')
    439439            set(handles.ZMax,'Visible','on')
    440         end
     440%         end
    441441        if isequal(ProjMode,'interp')|| isequal(ProjMode,'filter')
    442442            set(handles.DX,'Visible','on')
    443443            set(handles.DY,'Visible','on')
     444            set(handles.DZ,'Visible','on')
    444445        else
    445446            set(handles.DX,'Visible','off')
    446447            set(handles.DY,'Visible','off')
    447         end
    448         if isequal(ObjectStyle,'volume') && isequal(ProjMode,'interp')
    449             set(handles.DZ,'Visible','on') 
     448            set(handles.DZ,'Visible','off')
    450449        end
    451450end
    … …  
    775774
    776775%% plot the field projected on the object and store in the corresponding figue
    777 ProjData= proj_field(UvData.Field,ObjectData);%project the current interface field on ObjectData
     776'TESTproj'
     777ProjData= proj_field(UvData.Field,ObjectData)%project the current interface field on ObjectData
    778778PlotParam=read_plot_param(PlotHandles);
    779779[PlotType,Object_out{IndexObj}.PlotParam,plotaxes]=plot_field(ProjData,plotaxes,PlotParam);%update an existing field plot

  • trunk/src/uvmat.m

    r199 r206  
    23892389    return
    23902390end
    2391 [CellVarIndex,NbDim,VarType,errormsg]=find_field_indices(UvData.Field);
     2391'testUVMAT'
     2392UvData.Field
     2393[CellVarIndex,NbDim,VarType,errormsg]=find_field_indices(UvData.Field)
    23922394if ~isempty(errormsg)
    23932395    errormsg=['error in uvmat/refresh_field/find_field_indices: ' errormsg];
    … …  
    23952397end
    23962398[NbDim,imax]=max(NbDim);
     2399VarType{imax}
    23972400if ~isempty(VarType{imax}.coord_x)  && ~isempty(VarType{imax}.coord_y)    %unstructured coordinates
    23982401    XName=UvData.Field.ListVarName{VarType{imax}.coord_x};
    … …  
    24002403    eval(['nbvec=length(UvData.Field.' XName ');'])%nbre of measurement points (e.g. vectors)
    24012404    test_x=1;%test for unstructured coordinates
     2405    if ~isempty(VarType{imax}.coord_z)
     2406        ZName=UvData.Field.ListVarName{VarType{imax}.coord_z};
     2407    else
     2408       NbDim=2;
     2409    end
    24022410elseif VarType{imax}.coord(NbDim)>0 %structured coordinate 
    24032411    XName=UvData.Field.ListVarName{VarType{imax}.coord(NbDim)};
    … …  
    24072415end
    24082416if NbDim==3
    2409     ZName=UvData.Field.ListVarName{VarType{imax}.coord(1)};%structured coordinates in 3D
     2417    if ~test_x
     2418        ZName=UvData.Field.ListVarName{VarType{imax}.coord(1)};%structured coordinates in 3D
     2419    end
    24102420    eval(['ZMax=max(UvData.Field.' ZName ');'])
    24112421    eval(['ZMin=min(UvData.Field.' ZName ');'])
     2422    UvData.Field.ZMax=ZMax;
     2423    UvData.Field.ZMin=ZMin;
    24122424    test_z=1;
    24132425    if isequal(ZMin,ZMax)%no z dependency
    … …  
    24172429end
    24182430if exist('XName','var')
    2419 eval(['XMax=max(UvData.Field.' XName ');'])
    2420 eval(['XMin=min(UvData.Field.' XName ');'])
    2421 UvData.Field.NbDim=NbDim;
    2422 UvData.Field.XMax=XMax;
    2423 UvData.Field.XMin=XMin;
    2424 if NbDim >1
    2425     eval(['YMax=max(UvData.Field.' YName ');'])
    2426     eval(['YMin=min(UvData.Field.' YName ');'])
    2427     UvData.Field.YMax=YMax;
    2428     UvData.Field.YMin=YMin;
    2429 end
    2430 eval(['nbvec=length(UvData.Field.' XName ');'])
    2431 if test_x %unstructured coordinates
    2432     if test_z
    2433         UvData.Field.Mesh=((XMax-XMin)*(YMax-YMin)*(ZMax-ZMin))/nbvec;% volume per vector
    2434         UvData.Field.Mesh=(UvData.Field.Mesh)^(1/3);
     2431    eval(['XMax=max(UvData.Field.' XName ');'])
     2432    eval(['XMin=min(UvData.Field.' XName ');'])
     2433    UvData.Field.NbDim=NbDim;
     2434    UvData.Field.XMax=XMax;
     2435    UvData.Field.XMin=XMin;
     2436    if NbDim >1
     2437        eval(['YMax=max(UvData.Field.' YName ');'])
     2438        eval(['YMin=min(UvData.Field.' YName ');'])
     2439        UvData.Field.YMax=YMax;
     2440        UvData.Field.YMin=YMin;
     2441    end
     2442    eval(['nbvec=length(UvData.Field.' XName ');'])
     2443    if test_x %unstructured coordinates
     2444        if test_z
     2445            UvData.Field.Mesh=((XMax-XMin)*(YMax-YMin)*(ZMax-ZMin))/nbvec;% volume per vector
     2446            UvData.Field.Mesh=(UvData.Field.Mesh)^(1/3);
     2447        else
     2448            UvData.Field.Mesh=sqrt((XMax-XMin)*(YMax-YMin)/nbvec);%2D
     2449        end
    24352450    else
    2436         UvData.Field.Mesh=sqrt((XMax-XMin)*(YMax-YMin)/nbvec);%2D
    2437     end
    2438 else
    2439     VarIndex=CellVarIndex{imax}; % list of variable indices
    2440     DimIndex=UvData.Field.VarDimIndex{VarIndex(1)}; %list of dim indices for the variable
    2441     nbpoints_x=UvData.Field.DimValue(DimIndex(NbDim));
    2442     DX=(XMax-XMin)/(nbpoints_x-1);
    2443     if NbDim >1
    2444         nbpoints_y=UvData.Field.DimValue(DimIndex(NbDim-1));
    2445         DY=(YMax-YMin)/(nbpoints_y-1);
    2446     end
    2447     if NbDim==3
    2448         nbpoints_z=UvData.Field.DimValue(DimIndex(1));
    2449         DZ=(ZMax-ZMin)/(nbpoints_z-1);
    2450         UvData.Field.Mesh=sqrt(DX*DY*DZ);
    2451         UvData.Field.ZMax=ZMax;
    2452         UvData.Field.ZMin=ZMin;
    2453     else
    2454         UvData.Field.Mesh=sqrt(DX*DY);
    2455     end
    2456 end
     2451        VarIndex=CellVarIndex{imax}; % list of variable indices
     2452        DimIndex=UvData.Field.VarDimIndex{VarIndex(1)}; %list of dim indices for the variable
     2453        nbpoints_x=UvData.Field.DimValue(DimIndex(NbDim));
     2454        DX=(XMax-XMin)/(nbpoints_x-1);
     2455        if NbDim >1
     2456            nbpoints_y=UvData.Field.DimValue(DimIndex(NbDim-1));
     2457            DY=(YMax-YMin)/(nbpoints_y-1);
     2458        end
     2459        if NbDim==3
     2460            nbpoints_z=UvData.Field.DimValue(DimIndex(1));
     2461            DZ=(ZMax-ZMin)/(nbpoints_z-1);
     2462            UvData.Field.Mesh=sqrt(DX*DY*DZ);
     2463            UvData.Field.ZMax=ZMax;
     2464            UvData.Field.ZMin=ZMin;
     2465        else
     2466            UvData.Field.Mesh=sqrt(DX*DY);
     2467        end
     2468    end
    24572469end
    24582470
    … …  
    24812493        UvData.Object{1}.RangeZ=UvData.Field.Mesh;%main plotting plane
    24822494        UvData.Object{1}.Coord(1,3)=(UvData.Field.ZMin+UvData.Field.ZMax)/2;%section at a middle plane chosen
    2483         UvData.Object{1}.Phi=0;
    2484         UvData.Object{1}.Theta=0;
    2485         UvData.Object{1}.Psi=0;
     2495        UvData.Object{1}.Angle=[0 0 0];
     2496%         UvData.Object{1}.Theta=0;
     2497%         UvData.Object{1}.Psi=0;
    24862498        UvData.Object{1}.HandlesDisplay=plot(0,0,'Tag','proj_object');% A REVOIR
    2487         PlotHandles=get_plot_handles(handles);
     2499%         PlotHandles=get_plot_handles(handles);
    24882500        UvData.Object{1}.Name='1-PLANE';
    24892501        UvData.Object{1}.enable_plot=1;
    2490         set_object(UvData.Object{1},PlotHandles,ZBounds);
     2502        set_object(UvData.Object{1},handles,ZBounds);
    24912503        set(handles.list_object_1,'Value',1);
    24922504        set(handles.list_object_1,'String',{'1-PLANE'});
    … …  
    25042516            siz=size(XmlData.PlaneAngle);
    25052517            indangle=min(siz(1),UvData.ZIndex);%take first angle if a single angle is defined (translating scanning)
    2506             UvData.Object{1}.Phi=XmlData.PlaneAngle(indangle,1);
    2507             UvData.Object{1}.Theta=XmlData.PlaneAngle(indangle,2);
    2508             UvData.Object{1}.Psi=XmlData.PlaneAngle(indangle,3);
     2518            UvData.Object{1}.PlaneAngle=XmlData.PlaneAngle(indangle,:);
    25092519        end
    25102520    elseif isfield(UvData,'ZIndex')
    … …  
    25622572%loop on the projection objects: one or two
    25632573for imap=1:numel(IndexObj)
    2564     iobj=IndexObj(imap);   
     2574    iobj=IndexObj(imap);
    25652575    [ObjectData,errormsg]=proj_field(UvData.Field,UvData.Object{iobj});% project field on the object
    25662576    if testnewseries && isfield(ObjectData,'CoordUnit')
Note: See TracChangeset for help on using the changeset viewer.