Changeset 206 for trunk/src/proj_field.m


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

File:
1 edited

Legend:

Unmodified
Added
Removed

  • 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)
Note: See TracChangeset for help on using the changeset viewer.