Changeset 227 for trunk/src/proj_field.m


Ignore:
Timestamp:
Mar 31, 2011, 1:42:51 PM (13 years ago)
Author:
sommeria
Message:

add function sub_field_series to apply the sub_field operation to a series of fileds (for instance subtracting a background to an image series)

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/proj_field.m

    r215 r227  
    946946cos_om=1;
    947947sin_om=0;
     948test90x=0;%=1 for 90 degree rotation alround x axis
     949test90y=0;%=1 for 90 degree rotation alround y axis
    948950if isfield(ObjectData,'Angle')&& isequal(size(ObjectData.Angle),[1 3])&& ~isequal(ObjectData.Angle,[0 0 0])
     951    test90y=isequal(ObjectData.Angle,[0 90 0]);
    949952    PlaneAngle=(pi/180)*ObjectData.Angle;
    950953    om=norm(PlaneAngle);%norm of rotation angle in radians
    … …  
    958961    norm_plane(3)=OmAxis(3)*coeff+cos_om;
    959962end
    960 testangle=~isequal(PlaneAngle,[0 0 0]);
     963testangle=~isequal(PlaneAngle,[0 0 0]);% && ~test90y && ~test90x;%=1 for slanted plane
     964
    961965% Phi=0;%default
    962966% Theta=0;
    … …  
    10711075        DimCell={DimCell};%name of dimensions
    10721076    end
    1073 
    1074 %% case of input fields with unstructured coordinates
     1077   
     1078    %% case of input fields with unstructured coordinates
    10751079    if testX
    10761080        XName=FieldData.ListVarName{ivar_X};
    … …  
    10821086            eval(['coord_z=FieldData.' ZName ';'])
    10831087        end
    1084 
     1088       
    10851089        % translate  initial coordinates
    10861090        coord_x=coord_x-ObjectData.Coord(1,1);
    … …  
    10931097        if length(ivar_Z)==1 &&  width > 0
    10941098            %components of the unitiy vector normal to the projection plane
    1095             fieldZ=norm_plane(1)*coord_x + norm_plane(2)*coord_y+ norm_plane(3)*coord_z;% distance to the plane           
     1099            fieldZ=norm_plane(1)*coord_x + norm_plane(2)*coord_y+ norm_plane(3)*coord_z;% distance to the plane
    10961100            indcut=find(abs(fieldZ) <= width);
    10971101            size(indcut)
    10981102            for ivar=VarIndex
    10991103                VarName=FieldData.ListVarName{ivar};
    1100                 eval(['FieldData.' VarName '=FieldData.' VarName '(indcut);']) 
    1101                     % A VOIR : CAS DE VAR STRUCTUREE MAIS PAS GRILLE REGULIERE : INTERPOLER SUR GRILLE REGULIERE             
     1104                eval(['FieldData.' VarName '=FieldData.' VarName '(indcut);'])
     1105                % A VOIR : CAS DE VAR STRUCTUREE MAIS PAS GRILLE REGULIERE : INTERPOLER SUR GRILLE REGULIERE
    11021106            end
    11031107            coord_x=coord_x(indcut);
    … …  
    11051109            coord_z=coord_z(indcut);
    11061110        end
    1107 
    1108        %rotate coordinates if needed
    1109         if testangle
    1110 %             coord_X=coord_x;
    1111 %             coord_Y=coord_y;
    1112 %             if ~isequal(Theta,0)
    1113 %                 coord_Y=coord_Y *cos(Theta);
    1114 %             end
    1115 %         else
     1111       
     1112        %rotate coordinates if needed:
     1113        if testangle && ~test90y && ~test90x;%=1 for slanted plane
     1114            %             coord_X=coord_x;
     1115            %             coord_Y=coord_y;
     1116            %             if ~isequal(Theta,0)
     1117            %                 coord_Y=coord_Y *cos(Theta);
     1118            %             end
     1119            %         else
    11161120            coord_X=(coord_x *cos(Phi) + coord_y* sin(Phi));
    11171121            coord_Y=(-coord_x *sin(Phi) + coord_y *cos(Phi))*cos(Theta);
    1118 %         end
    1119 %         if ~isempty(ivar_Z)
     1122            %         end
     1123            %         if ~isempty(ivar_Z)
    11201124            coord_Y=coord_Y+coord_z *sin(Theta);
    1121 %         end
    1122 %         if testangle
    1123                 coord_X=(coord_X *cos(Psi) - coord_Y* sin(Psi));%A VERIFIER
    1124                 coord_Y=(coord_X *sin(Psi) + coord_Y* cos(Psi));
     1125            %         end
     1126            %         if testangle
     1127            coord_X=(coord_X *cos(Psi) - coord_Y* sin(Psi));%A VERIFIER
     1128            coord_Y=(coord_X *sin(Psi) + coord_Y* cos(Psi));
    11251129        else
    11261130            coord_X=coord_x;
    … …  
    11511155            for ivar=VarIndex
    11521156                VarName=FieldData.ListVarName{ivar};
    1153                 eval(['FieldData.' VarName '=FieldData.' VarName '(indcut);'])           
     1157                eval(['FieldData.' VarName '=FieldData.' VarName '(indcut);'])
    11541158            end
    11551159            coord_X=coord_X(indcut);
    … …  
    11641168            %ProjData.ListDimName=[ProjData.ListDimName FieldData.VarDimName(VarIndex(1))];%add the point index to the list of dimensions
    11651169            %ProjData.DimValue=[ProjData.
    1166              %length(coord_X)];
    1167 
     1170            %length(coord_X)];
     1171           
    11681172            for ivar=VarIndex %transfer variables to the projection plane
    11691173                VarName=FieldData.ListVarName{ivar};
    … …  
    11751179                    eval(['ProjData.' VarName '=FieldData.' VarName ';'])
    11761180                end
    1177                 if isempty(ivar_Z) || ivar~=ivar_Z 
     1181                if isempty(ivar_Z) || ivar~=ivar_Z
    11781182                    ProjData.ListVarName=[ProjData.ListVarName VarName];
    11791183                    ProjData.VarDimName=[ProjData.VarDimName DimCell];
    … …  
    11831187                    end
    11841188                end
    1185             end 
     1189            end
    11861190        elseif isequal(ObjectData.ProjMode,'interp')||isequal(ObjectData.ProjMode,'filter')%interpolate data on a regular grid
    11871191            coord_x_proj=XMin:DX:XMax;
    … …  
    11891193            DimCell={'coord_y','coord_x'};
    11901194            ProjData.ListVarName={'coord_y','coord_x'};
    1191             ProjData.VarDimName={'coord_y','coord_x'};   
    1192             nbcoord=2; 
     1195            ProjData.VarDimName={'coord_y','coord_x'};
     1196            nbcoord=2;
    11931197            ProjData.coord_y=[YMin YMax];
    11941198            ProjData.coord_x=[XMin XMax];
    … …  
    12111215            for ivar=VarIndex
    12121216                VarName=FieldData.ListVarName{ivar};
    1213                 if ~( ivar==ivar_X || ivar==ivar_Y || ivar==ivar_Z || ivar==ivar_F || ivar==ivar_FF || test_anc(ivar)==1)                 
     1217                if ~( ivar==ivar_X || ivar==ivar_Y || ivar==ivar_Z || ivar==ivar_F || ivar==ivar_FF || test_anc(ivar)==1)
    12141218                    ivar_new=ivar_new+1;
    12151219                    ProjData.ListVarName=[ProjData.ListVarName {VarName}];
    … …  
    12451249                ProjData.FF=reshape(FF,length(coord_y_proj),length(coord_x_proj));
    12461250                ProjData.ListVarName=[ProjData.ListVarName {'FF'}];
    1247                ProjData.VarDimName=[ProjData.VarDimName {DimCell}];
     1251                ProjData.VarDimName=[ProjData.VarDimName {DimCell}];
    12481252                ProjData.VarAttribute{ivar_new+1+nbcoord}.Role='errorflag';
    12491253            end
    12501254        end
    12511255       
    1252 %% case of input fields defined on a structured  grid
     1256        %% case of input fields defined on a structured  grid
    12531257    else
    1254         VarName=FieldData.ListVarName{VarIndex(1)};%get the first variable of the cell to get the input matrix dimensions
     1258        'TESTproj'
     1259        VarName=FieldData.ListVarName{VarIndex(1)}%get the first variable of the cell to get the input matrix dimensions
    12551260        eval(['DimValue=size(FieldData.' VarName ');'])%input matrix dimensions
    1256         DimValue(DimValue==1)=[];%remove singleton dimensions       
     1261        DimValue(DimValue==1)=[];%remove singleton dimensions
    12571262        NbDim=numel(DimValue);%update number of space dimensions
    12581263        nbcolor=1; %default number of 'color' components: third matrix index without corresponding coordinate
    … …  
    12701275            end
    12711276        end
    1272         AYName=FieldData.ListVarName{VarType.coord(NbDim-1)};%name of input x coordinate (name preserved on projection)
    1273         AXName=FieldData.ListVarName{VarType.coord(NbDim)};%name of input y coordinate (name preserved on projection)   
     1277        testangle
     1278        if testangle% TODO modify name also in case of origin shift in x or y
     1279            AYName='Y';
     1280            AXName='X';
     1281            count=0;
     1282            %modify coordinate names if they are already used
     1283            while ~(isempty(find(strcmp('AXName',ProjData.ListVarName),1)) && isempty(find(strcmp('AYName',ProjData.ListVarName),1)))
     1284                count=count+1;
     1285                AYName=[AYName '_' num2str(count)];
     1286                AXName=[AXName '_' num2str(count)];
     1287            end
     1288        else
     1289            AYName=FieldData.ListVarName{VarType.coord(NbDim-1)}%name of input x coordinate (name preserved on projection)
     1290            AXName=FieldData.ListVarName{VarType.coord(NbDim)}%name of input y coordinate (name preserved on projection)
     1291        end
    12741292        eval(['AX=FieldData.' AXName ';'])
    12751293        eval(['AY=FieldData.' AYName ';'])
    12761294        ListDimName=FieldData.VarDimName{VarIndex(1)};
    1277         ProjData.ListVarName=[ProjData.ListVarName {AYName} {AXName}]; %TODO: check if it already exists in Projdata (several cells)
    1278         ProjData.VarDimName=[ProjData.VarDimName {AYName} {AXName}];
    1279 
    1280 %         for idim=1:length(ListDimName)
    1281 %             DimName=ListDimName{idim};
    1282 %             if strcmp(DimName,'rgb')||strcmp(DimName,'nb_coord')||strcmp(DimName,'nb_coord_i')
    1283 %                nbcolor=DimValue(idim);
    1284 %                DimValue(idim)=[];
    1285 %             end
    1286 %             if isequal(DimName,'nb_coord_j')% NOTE: CASE OF TENSOR NOT TREATED
    1287 %                 DimValue(idim)=[];
    1288 %             end
    1289 %         end 
     1295            ProjData.ListVarName=[ProjData.ListVarName {AYName} {AXName}]; %TODO: check if it already exists in Projdata (several cells)
     1296            ProjData.VarDimName=[ProjData.VarDimName {AYName} {AXName}];
    12901297        Coord_z=[];
    12911298        Coord_y=[];
    1292         Coord_x=[];   
    1293 
     1299        Coord_x=[];
     1300       
    12941301        for idim=1:NbDim %loop on space dimensions
    12951302            test_interp(idim)=0;%test for coordiate interpolation (non regular grid), =0 by default
    … …  
    12981305                eval(['Coord{idim}=FieldData.' FieldData.ListVarName{ivar} ';']) ;% coord values for the input field
    12991306                if numel(Coord{idim})==2 %input array defined on a regular grid
    1300                    DCoord_min(idim)=(Coord{idim}(2)-Coord{idim}(1))/DimValue(idim);
     1307                    DCoord_min(idim)=(Coord{idim}(2)-Coord{idim}(1))/DimValue(idim);
    13011308                else
    13021309                    DCoord=diff(Coord{idim});%array of coordinate derivatives for the input field
    13031310                    DCoord_min(idim)=min(DCoord);
    13041311                    DCoord_max=max(DCoord);
    1305                 %    test_direct(idim)=DCoord_max>0;% =1 for increasing values, 0 otherwise
    1306                     if abs(DCoord_max-DCoord_min(idim))>abs(DCoord_max/1000) 
     1312                    %    test_direct(idim)=DCoord_max>0;% =1 for increasing values, 0 otherwise
     1313                    if abs(DCoord_max-DCoord_min(idim))>abs(DCoord_max/1000)
    13071314                        msgbox_uvmat('ERROR',['non monotonic dimension variable # ' num2str(idim)  ' in proj_field.m'])
    1308                                 return
    1309                     end               
     1315                        return
     1316                    end
    13101317                    test_interp(idim)=(DCoord_max-DCoord_min(idim))> 0.0001*abs(DCoord_max);% test grid regularity
    13111318                end
    … …  
    13211328            DY=abs(DCoord_min(NbDim-1));
    13221329        end
    1323         npY=1+round(abs(Coord{NbDim-1}(end)-Coord{NbDim-1}(1))/DY);%nbre of points after interpol 
     1330        npY=1+round(abs(Coord{NbDim-1}(end)-Coord{NbDim-1}(1))/DY);%nbre of points after interpol
    13241331        if DX==0
    13251332            DX=abs(DCoord_min(NbDim));
    13261333        end
    1327         npX=1+round(abs(Coord{NbDim}(end)-Coord{NbDim}(1))/DX);%nbre of points after interpol 
     1334        npX=1+round(abs(Coord{NbDim}(end)-Coord{NbDim}(1))/DX);%nbre of points after interpol
    13281335        for idim=1:NbDim
    13291336            if test_interp(idim)
    13301337                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
    13311338            end
    1332         end       
     1339        end
    13331340        Coord_y=linspace(Coord{NbDim-1}(1),Coord{NbDim-1}(end),npY);
    13341341        test_direct_y=test_direct(NbDim-1);
    … …  
    13361343        test_direct_x=test_direct(NbDim);
    13371344        DAX=DCoord_min(NbDim);
    1338         DAY=DCoord_min(NbDim-1); 
     1345        DAY=DCoord_min(NbDim-1);
    13391346        minAX=min(Coord_x);
    13401347        maxAX=max(Coord_x);
    … …  
    13741381        end
    13751382        npX=floor((XMax-XMin)/DX+1);
    1376         npY=floor((YMax-YMin)/DY+1);   
     1383        npY=floor((YMax-YMin)/DY+1);
    13771384        if test_direct_y
    13781385            coord_y_proj=linspace(YMin,YMax,npY);%abscissa of the new pixels along the line
    … …  
    13841391        else
    13851392            coord_x_proj=linspace(XMax,XMin,npX);%abscissa of the new pixels along the line
    1386         end
    1387        
     1393        end
    13881394        % case with no rotation and interpolation
    13891395        if isequal(ProjMode,'projection') && ~testangle
    13901396            if ~testXMin && ~testXMax && ~testYMin && ~testYMax && NbDim==2
    1391                 ProjData=FieldData; 
     1397                ProjData=FieldData;
    13921398            else
    13931399                indY=NbDim-1;
    … …  
    14021408                    Ybound(2)=Coord{indY}(1)-DYinit*(max_indy-1);
    14031409                    Ybound(1)=Coord{indY}(1)-DYinit*(min_indy-1);
    1404                 end   
     1410                end
    14051411                if test_direct(NbDim)==1
    14061412                    min_indx=ceil((XMin-Coord{NbDim}(1))/DXinit)+1;
    … …  
    14121418                    max_indx=floor((Coord{NbDim}(1)-XMin)/DXinit)+1;
    14131419                    Xbound(2)=Coord{NbDim}(1)+DXinit*(max_indx-1);
     1420                   
     1421                   
    14141422                    Xbound(1)=Coord{NbDim}(1)+DXinit*(min_indx-1);
    1415                 end
    1416                 if NbDim==3
    1417                     DimCell(1)=[]; %suppress z variable
    1418                     DimValue(1)=[];
    1419                                         %structured coordinates
    1420                     if test_direct(1)
    1421                         iz=ceil((ObjectData.Coord(1,3)-Coord{1}(1))/DZ)+1;
    1422                     else
    1423                         iz=ceil((Coord{1}(1)-ObjectData.Coord(1,3))/DZ)+1;
    1424                     end
    14251423                end
    14261424                min_indy=max(min_indy,1);% deals with margin (bound lower than the first index)
    14271425                min_indx=max(min_indx,1);
    1428                 max_indy=min(max_indy,DimValue(1));
    1429                 max_indx=min(max_indx,DimValue(2));
    1430                 for ivar=VarIndex% loop on non coordinate variables
    1431                     VarName=FieldData.ListVarName{ivar};
    1432                     ProjData.ListVarName=[ProjData.ListVarName VarName];
    1433                     ProjData.VarDimName=[ProjData.VarDimName {DimCell}];
    1434                     if isfield(FieldData,'VarAttribute') && length(FieldData.VarAttribute)>=ivar
    1435                         ProjData.VarAttribute{length(ProjData.ListVarName)}=FieldData.VarAttribute{ivar};
    1436                     end
     1426
     1427                if test90y
     1428                    'TEST3D'
     1429                    ind_new=[3 2 1];
     1430                    DimCell=DimCell(ind_new);
     1431                    DimValue=DimValue(ind_new);
     1432                    DimCell(1)=[]; %suppress x variable
     1433                    DimValue(1)=[];
     1434                    iz=ceil((ObjectData.Coord(1,1)-Coord{3}(1))/DX)+1;
     1435                    for ivar=VarIndex
     1436                        VarName=FieldData.ListVarName{ivar};
     1437                        ProjData.ListVarName=[ProjData.ListVarName VarName];
     1438                        ProjData.VarDimName=[ProjData.VarDimName {DimCell}];
     1439                        ProjData.VarAttribute{length(ProjData.ListVarName)}=FieldData.VarAttribute{ivar}; %reproduce the variable attributes
     1440                        eval(['size(FieldData.' VarName ')'])
     1441                        eval(['ProjData.' VarName '=permute(FieldData.' VarName ',ind_new)'])% permute x and z indices for 90 degree rotation
     1442                        eval(['size(ProjData.' VarName ')'])
     1443                        eval(['ProjData.' VarName '=squeeze(FieldData.' VarName '(iz,:,:));'])% select the z index iz
     1444                    end
     1445                    eval(['ProjData.' AYName '=[Ybound(1) Ybound(2)];']) %record the new (projected ) y coordinates
     1446                    eval(['ProjData.' AXName '=[Coord{1}(end),Coord{1}(1)];']) %record the new (projected ) x coordinates
     1447                else
    14371448                    if NbDim==3
    1438                         eval(['ProjData.' VarName '=squeeze(FieldData.' VarName '(iz,min_indy:max_indy,min_indx:max_indx));']);
    1439                     else
    1440                         eval(['ProjData.' VarName '=FieldData.' VarName '(min_indy:max_indy,min_indx:max_indx,:);']);
    1441                     end
    1442                 end 
    1443                 eval(['ProjData.' AYName '=[Ybound(1) Ybound(2)];']) %record the new (projected ) y coordinates
    1444                 eval(['ProjData.' AXName '=[Xbound(1) Xbound(2)];']) %record the new (projected ) x coordinates
     1449                        DimCell(1)=[]; %suppress z variable
     1450                        DimValue(1)=[];
     1451                        if test_direct(1)
     1452                            iz=ceil((ObjectData.Coord(1,3)-Coord{1}(1))/DZ)+1;
     1453                        else
     1454                            iz=ceil((Coord{1}(1)-ObjectData.Coord(1,3))/DZ)+1;
     1455                        end
     1456                    end
     1457                    max_indy=min(max_indy,DimValue(1));%introduce bounds in y and x indices
     1458                    max_indx=min(max_indx,DimValue(2));
     1459                    for ivar=VarIndex% loop on non coordinate variables
     1460                        VarName=FieldData.ListVarName{ivar};
     1461                        ProjData.ListVarName=[ProjData.ListVarName VarName];
     1462                        ProjData.VarDimName=[ProjData.VarDimName {DimCell}];
     1463                        if isfield(FieldData,'VarAttribute') && length(FieldData.VarAttribute)>=ivar
     1464                            ProjData.VarAttribute{length(ProjData.ListVarName)}=FieldData.VarAttribute{ivar};
     1465                        end
     1466                        if NbDim==3
     1467                            eval(['ProjData.' VarName '=squeeze(FieldData.' VarName '(iz,min_indy:max_indy,min_indx:max_indx));']);
     1468                        else
     1469                            eval(['ProjData.' VarName '=FieldData.' VarName '(min_indy:max_indy,min_indx:max_indx,:);']);
     1470                        end
     1471                    end
     1472                    eval(['ProjData.' AYName '=[Ybound(1) Ybound(2)];']) %record the new (projected ) y coordinates
     1473                    eval(['ProjData.' AXName '=[Xbound(1) Xbound(2)];']) %record the new (projected ) x coordinates
     1474                end
    14451475            end
    14461476        else       % case with rotation and/or interpolation
    … …  
    14531483                XIMA=reshape(round(XIMA),1,npX*npY);%indices reorganized in 'line'
    14541484                YIMA=reshape(round(YIMA),1,npX*npY);
    1455                 flagin=XIMA>=1 & XIMA<=DimValue(2) & YIMA >=1 & YIMA<=DimValue(1);%flagin=1 inside the original image 
     1485                flagin=XIMA>=1 & XIMA<=DimValue(2) & YIMA >=1 & YIMA<=DimValue(1);%flagin=1 inside the original image
    14561486                if isequal(ObjectData.ProjMode,'filter')
    14571487                    npx_filter=ceil(abs(DX/DAX));
    … …  
    14661496                for ivar=VarIndex
    14671497                    VarName=FieldData.ListVarName{ivar};
    1468                     if test_interp(1) || test_interp(2)%interpolate on a regular grid       
    1469                           eval(['ProjData.' VarName '=interp2(Coord{2},Coord{1},FieldData.' VarName ',Coord_x,Coord_y'');']) %TO TEST
     1498                    if test_interp(1) || test_interp(2)%interpolate on a regular grid
     1499                        eval(['ProjData.' VarName '=interp2(Coord{2},Coord{1},FieldData.' VarName ',Coord_x,Coord_y'');']) %TO TEST
    14701500                    end
    14711501                    %filter the field (image) if option 'filter' is used
    1472                     if test_filter 
    1473                          Aclass=class(FieldData.A);
    1474                          eval(['ProjData.' VarName '=filter2(Mfilter,FieldData.' VarName ',''valid'');'])
    1475                          if ~isequal(Aclass,'double')
    1476                              eval(['ProjData.' VarName '=' Aclass '(FieldData.' VarName ');'])%revert to integer values
    1477                          end
    1478                     end
    1479                     eval(['vec_A=reshape(FieldData.' VarName ',[],nbcolor);'])%put the original image in line             
     1502                    if test_filter
     1503                        Aclass=class(FieldData.A);
     1504                        eval(['ProjData.' VarName '=filter2(Mfilter,FieldData.' VarName ',''valid'');'])
     1505                        if ~isequal(Aclass,'double')
     1506                            eval(['ProjData.' VarName '=' Aclass '(FieldData.' VarName ');'])%revert to integer values
     1507                        end
     1508                    end
     1509                    eval(['vec_A=reshape(FieldData.' VarName ',[],nbcolor);'])%put the original image in line
    14801510                    %ind_in=find(flagin);
    14811511                    ind_out=find(~flagin);
    14821512                    ICOMB=(XIMA-1)*DimValue(1)+YIMA;
    14831513                    ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A
    1484                     vec_B(flagin,1:nbcolor)=vec_A(ICOMB,:); 
     1514                    vec_B(flagin,1:nbcolor)=vec_A(ICOMB,:);
    14851515                    for icolor=1:nbcolor
    14861516                        vec_B(ind_out,icolor)=zeros(size(ind_out));
    … …  
    14901520                    if isfield(FieldData,'VarAttribute')&&length(FieldData.VarAttribute)>=ivar
    14911521                        ProjData.VarAttribute{length(ProjData.ListVarName)+nbcoord}=FieldData.VarAttribute{ivar};
    1492                     end     
     1522                    end
    14931523                    eval(['ProjData.' VarName '=reshape(vec_B,npY,npX,nbcolor);']);
    14941524                end
    1495                 ProjData.FF=reshape(~flagin,npY,npX);%false flag A FAIRE: tenir compte d'un flga antérieur 
     1525                ProjData.FF=reshape(~flagin,npY,npX);%false flag A FAIRE: tenir compte d'un flga antérieur
    14961526                ProjData.ListVarName=[ProjData.ListVarName 'FF'];
    14971527                ProjData.VarDimName=[ProjData.VarDimName {DimCell}];
    14981528                ProjData.VarAttribute{length(ProjData.ListVarName)}.Role='errorflag';
    1499             else %3D case
    1500                 if ~testangle     
    1501                     % unstructured z coordinate
    1502                     test_sup=(Coord{1}>=ObjectData.Coord(1,3));
    1503                     iz_sup=find(test_sup);
    1504                     iz=iz_sup(1);
    1505                     if iz>=1 & iz<=npz
    1506                         %ProjData.ListDimName=[ProjData.ListDimName ListDimName(2:end)];
    1507                         %ProjData.DimValue=[ProjData.DimValue npY npX];
    1508                         for ivar=VarIndex
    1509                             VarName=FieldData.ListVarName{ivar};
    1510                             ProjData.ListVarName=[ProjData.ListVarName VarName];
    1511                             ProjData.VarAttribute{length(ProjData.ListVarName)}=FieldData.VarAttribute{ivar}; %reproduce the variable attributes 
    1512                             eval(['ProjData.' VarName '=squeeze(FieldData.' VarName '(iz,:,:));'])% select the z index iz
    1513                             %TODO : do a vertical average for a thick plane
    1514                             if test_interp(2) || test_interp(3)
    1515                                 eval(['ProjData.' VarName '=interp2(Coord{3},Coord{2},ProjData.' VarName ',Coord_x,Coord_y'');'])
    1516                             end
     1529            elseif ~testangle
     1530                % unstructured z coordinate
     1531                test_sup=(Coord{1}>=ObjectData.Coord(1,3));
     1532                iz_sup=find(test_sup);
     1533                iz=iz_sup(1);
     1534                if iz>=1 & iz<=npz
     1535                    %ProjData.ListDimName=[ProjData.ListDimName ListDimName(2:end)];
     1536                    %ProjData.DimValue=[ProjData.DimValue npY npX];
     1537                    for ivar=VarIndex
     1538                        VarName=FieldData.ListVarName{ivar};
     1539                        ProjData.ListVarName=[ProjData.ListVarName VarName];
     1540                        ProjData.VarAttribute{length(ProjData.ListVarName)}=FieldData.VarAttribute{ivar}; %reproduce the variable attributes
     1541                        eval(['ProjData.' VarName '=squeeze(FieldData.' VarName '(iz,:,:));'])% select the z index iz
     1542                        %TODO : do a vertical average for a thick plane
     1543                        if test_interp(2) || test_interp(3)
     1544                            eval(['ProjData.' VarName '=interp2(Coord{3},Coord{2},ProjData.' VarName ',Coord_x,Coord_y'');'])
    15171545                        end
    15181546                    end
    1519                 else
    1520                     errormsg='projection of structured coordinates on oblique plane not yet implemented';
    1521                     %TODO: use interp3
    1522                     return
    1523                 end
    1524             end
    1525         end
    1526     end
    1527 
     1547                end
     1548            else
     1549                errormsg='projection of structured coordinates on oblique plane not yet implemented';
     1550                %TODO: use interp3
     1551                return
     1552            end
     1553        end
     1554    end
     1555   
    15281556    %% projection of  velocity components in the rotated coordinates
    15291557    if testangle && length(ivar_U)==1
    … …  
    15331561        end
    15341562        UName=FieldData.ListVarName{ivar_U};
    1535         VName=FieldData.ListVarName{ivar_V};   
     1563        VName=FieldData.ListVarName{ivar_V};
    15361564        eval(['ProjData.' UName  '=cos(PlaneAngle(3))*ProjData.' UName '+ sin(PlaneAngle(3))*ProjData.' VName ';'])
    15371565        eval(['ProjData.' VName  '=cos(Theta)*(-sin(PlaneAngle(3))*ProjData.' UName '+ cos(PlaneAngle(3))*ProjData.' VName ');'])
    15381566        if ~isempty(ivar_W)
    15391567            WName=FieldData.ListVarName{ivar_W};
    1540             eval(['ProjData.' VName '=ProjData.' VName '+ ProjData.' WName '*sin(Theta);'])% 
     1568            eval(['ProjData.' VName '=ProjData.' VName '+ ProjData.' WName '*sin(Theta);'])%
    15411569            eval(['ProjData.' WName '=NormVec_X*ProjData.' UName '+ NormVec_Y*ProjData.' VName '+ NormVec_Z* ProjData.' WName ';']);
    15421570        end
    … …  
    18571885        VarName=FieldData.ListVarName{VarIndex(1)};%get the first variable of the cell to get the input matrix dimensions
    18581886        eval(['DimValue=size(FieldData.' VarName ');'])%input matrix dimensions
    1859         DimValue(find(DimValue==1))=[];%remove singleton dimensions       
     1887        DimValue(DimValue==1)=[];%remove singleton dimensions       
    18601888        NbDim=numel(DimValue);%update number of space dimensions
    18611889        nbcolor=1; %default number of 'color' components: third matrix index without corresponding coordinate
    … …  
    18651893                return
    18661894            else
    1867                 VarType.coord
    18681895                if numel(find(VarType.coord))==2% the third matrix dimension does not correspond to a space coordinate
    18691896                    nbcolor=DimValue(3);
Note: See TracChangeset for help on using the changeset viewer.