Changeset 399 for trunk/src/proj_field.m


Ignore:
Timestamp:
Apr 27, 2012, 12:28:47 PM (12 years ago)
Author:
sommeria
Message:

implementation of thin plate interpolation (proj on planes with mode 'filter'), rationalisation of variable formats in civ_matlab

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/proj_field.m

    r397 r399  
    11%'proj_field': projects the field on a projection object
    22%--------------------------------------------------------------------------
    3 %  function [ProjData,errormsg]=proj_field(FieldData,ObjectData,FieldName)
     3%  function [ProjData,errormsg]=proj_field(FieldData,ObjectData)
    44%
    55% OUTPUT:
    … …  
    2727%FieldData: data of the field to be projected on the projection object, with optional fields
    2828%    .Txt: error message, transmitted to the projection
    29 %    .CoordType: 'px' or 'phys' type of coordinates of the field, must be the same as for the projection object, transmitted
     29%    .FieldList: cell array of strings representing the fields to calculate
    3030%    .Mesh: typical distance between data points (used for mouse action or display), transmitted
    3131%    .CoordUnit, .TimeUnit, .dt: transmitted
    … …  
    8080%AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
    8181
    82 function [ProjData,errormsg]=proj_field(FieldData,ObjectData,FieldName)
     82function [ProjData,errormsg]=proj_field(FieldData,ObjectData)
    8383errormsg='';%default
    8484if ~exist('FieldName','var')
    … …  
    8787%% case of no projection (object is used only as graph display)
    8888if isfield(ObjectData,'ProjMode') && (isequal(ObjectData.ProjMode,'none')||isequal(ObjectData.ProjMode,'mask_inside')||isequal(ObjectData.ProjMode,'mask_outside'))
    89     if ~isempty(FieldName)
    90         [ProjData,errormsg]=calc_field(FieldName,FieldData);
    91     else
    9289    ProjData=[];
    93     end
    9490    return
    9591end
    … …  
    143139        end
    144140    case 'plane'
    145             [ProjData,errormsg] = proj_plane(FieldData,ObjectData,FieldName);
     141            [ProjData,errormsg] = proj_plane(FieldData,ObjectData);
    146142    case 'volume'
    147143        [ProjData,errormsg] = proj_volume(FieldData,ObjectData);
    … …  
    892888%project on a plane
    893889% AJOUTER flux,circul,error
    894  function  [ProjData,errormsg] = proj_plane(FieldData, ObjectData,FieldName)
     890 function  [ProjData,errormsg] = proj_plane(FieldData, ObjectData)
    895891%-----------------------------------------------------------------
    896892
    … …  
    994990% CellVarIndex=cells of variable index arrays
    995991ivar_new=0; % index of the current variable in the projected field
    996 icoord=0;
     992% icoord=0;
    997993nbcoord=0;%number of added coordinate variables brought by projection
    998994nbvar=0;
    … …  
    1000996    NbDim=NbDimVec(icell);
    1001997    if NbDim<2
    1002         continue
     998        continue % only cells represnting 2D or 3D fields are involved
    1003999    end
    10041000    VarIndex=CellVarIndex{icell};%  indices of the selected variables in the list FieldData.ListVarName
    … …  
    10101006    ivar_V=VarType.vector_y;
    10111007    ivar_W=VarType.vector_z;
    1012     ivar_C=VarType.scalar ;
     1008    %    ivar_C=VarType.scalar ;
    10131009    ivar_Anc=VarType.ancillary;
    10141010    test_anc=zeros(size(VarIndex));
    … …  
    10161012    ivar_F=VarType.warnflag;
    10171013    ivar_FF=VarType.errorflag;
    1018     testX=~isempty(ivar_X) && ~isempty(ivar_Y);
     1014    check_unstructured_coord=(~isempty(ivar_X) && ~isempty(ivar_Y))||~isempty(VarType.coord_tps);
    10191015    DimCell=FieldData.VarDimName{VarIndex(1)};
    10201016    if ischar(DimCell)
    … …  
    10241020    %% case of input fields with unstructured coordinates
    10251021    coord_z=0;%default
    1026     if testX
    1027         XName=FieldData.ListVarName{ivar_X};
    1028         YName=FieldData.ListVarName{ivar_Y};
    1029         coord_x=FieldData.(XName);
    1030         coord_y=FieldData.(YName);
    1031         if length(ivar_Z)==1
    1032             ZName=FieldData.ListVarName{ivar_Z};
    1033             coord_z=FieldData.(ZName);
    1034         end
    1035        
    1036         % translate  initial coordinates
    1037         coord_x=coord_x-ObjectData.Coord(1,1);
    1038         coord_y=coord_y-ObjectData.Coord(1,2);
    1039         if ~isempty(ivar_Z)
    1040             coord_z=coord_z-ObjectData.Coord(1,3);
    1041         end
    1042        
    1043         % selection of the vectors in the projection range (3D case)
    1044         if length(ivar_Z)==1 &&  width > 0
    1045             %components of the unitiy vector normal to the projection plane
    1046             fieldZ=norm_plane(1)*coord_x + norm_plane(2)*coord_y+ norm_plane(3)*coord_z;% distance to the plane
    1047             indcut=find(abs(fieldZ) <= width);
    1048             size(indcut)
    1049             for ivar=VarIndex
    1050                 VarName=FieldData.ListVarName{ivar};
    1051                 eval(['FieldData.' VarName '=FieldData.' VarName '(indcut);'])
    1052                 % A VOIR : CAS DE VAR STRUCTUREE MAIS PAS GRILLE REGULIERE : INTERPOLER SUR GRILLE REGULIERE
    1053             end
    1054             coord_x=coord_x(indcut);
    1055             coord_y=coord_y(indcut);
    1056             coord_z=coord_z(indcut);
    1057         end
    1058        
    1059         %rotate coordinates if needed:
    1060         Psi=PlaneAngle(1);
    1061         Theta=PlaneAngle(2);
    1062         Phi=PlaneAngle(3);
    1063         if testangle && ~test90y && ~test90x;%=1 for slanted plane
    1064             coord_X=(coord_x *cos(Phi) + coord_y* sin(Phi));
    1065             coord_Y=(-coord_x *sin(Phi) + coord_y *cos(Phi))*cos(Theta);
    1066             coord_Y=coord_Y+coord_z *sin(Theta);
    1067             coord_X=(coord_X *cos(Psi) - coord_Y* sin(Psi));%A VERIFIER
    1068 
    1069             coord_Y=(coord_X *sin(Psi) + coord_Y* cos(Psi));
    1070         else
    1071             coord_X=coord_x;
    1072             coord_Y=coord_y;
    1073         end
    1074        
    1075         %restriction to the range of x and y if imposed
    1076         testin=ones(size(coord_X)); %default
    1077         testbound=0;
    1078         if testXMin
    1079             testin=testin & (coord_X >= XMin);
    1080             testbound=1;
    1081         end
    1082         if testXMax
    1083             testin=testin & (coord_X <= XMax);
    1084             testbound=1;
    1085         end
    1086         if testYMin
    1087             testin=testin & (coord_Y >= YMin);
    1088             testbound=1;
    1089         end
    1090         if testYMin
    1091             testin=testin & (coord_Y <= YMax);
    1092             testbound=1;
    1093         end
    1094         if testbound
    1095             indcut=find(testin);
    1096             for ivar=VarIndex
    1097                 VarName=FieldData.ListVarName{ivar};
    1098                 eval(['FieldData.' VarName '=FieldData.' VarName '(indcut);'])
    1099             end
    1100             coord_X=coord_X(indcut);
    1101             coord_Y=coord_Y(indcut);
     1022    if check_unstructured_coord
     1023        if ~isempty(ivar_X)% case of tps excluded. TODO similar procedure
     1024            XName=FieldData.ListVarName{ivar_X};
     1025            YName=FieldData.ListVarName{ivar_Y};
     1026            coord_x=FieldData.(XName);
     1027            coord_y=FieldData.(YName);
    11021028            if length(ivar_Z)==1
    1103                 coord_Z=coord_Z(indcut);
     1029                ZName=FieldData.ListVarName{ivar_Z};
     1030                coord_z=FieldData.(ZName);
     1031            end
     1032           
     1033            % translate  initial coordinates
     1034            coord_x=coord_x-ObjectData.Coord(1,1);
     1035            coord_y=coord_y-ObjectData.Coord(1,2);
     1036            if ~isempty(ivar_Z)
     1037                coord_z=coord_z-ObjectData.Coord(1,3);
     1038            end
     1039           
     1040            % selection of the vectors in the projection range (3D case)
     1041            if length(ivar_Z)==1 &&  width > 0
     1042                %components of the unitiy vector normal to the projection plane
     1043                fieldZ=norm_plane(1)*coord_x + norm_plane(2)*coord_y+ norm_plane(3)*coord_z;% distance to the plane
     1044                indcut=find(abs(fieldZ) <= width);
     1045                size(indcut)
     1046                for ivar=VarIndex
     1047                    VarName=FieldData.ListVarName{ivar};
     1048                    eval(['FieldData.' VarName '=FieldData.' VarName '(indcut);'])
     1049                    % A VOIR : CAS DE VAR STRUCTUREE MAIS PAS GRILLE REGULIERE : INTERPOLER SUR GRILLE REGULIERE
     1050                end
     1051                coord_x=coord_x(indcut);
     1052                coord_y=coord_y(indcut);
     1053                coord_z=coord_z(indcut);
     1054            end
     1055           
     1056            %rotate coordinates if needed:
     1057            Psi=PlaneAngle(1);
     1058            Theta=PlaneAngle(2);
     1059            Phi=PlaneAngle(3);
     1060            if testangle && ~test90y && ~test90x;%=1 for slanted plane
     1061                coord_X=(coord_x *cos(Phi) + coord_y* sin(Phi));
     1062                coord_Y=(-coord_x *sin(Phi) + coord_y *cos(Phi))*cos(Theta);
     1063                coord_Y=coord_Y+coord_z *sin(Theta);
     1064                coord_X=(coord_X *cos(Psi) - coord_Y* sin(Psi));%A VERIFIER
     1065               
     1066                coord_Y=(coord_X *sin(Psi) + coord_Y* cos(Psi));
     1067            else
     1068                coord_X=coord_x;
     1069                coord_Y=coord_y;
     1070            end
     1071           
     1072            %restriction to the range of x and y if imposed
     1073            testin=ones(size(coord_X)); %default
     1074            testbound=0;
     1075            if testXMin
     1076                testin=testin & (coord_X >= XMin);
     1077                testbound=1;
     1078            end
     1079            if testXMax
     1080                testin=testin & (coord_X <= XMax);
     1081                testbound=1;
     1082            end
     1083            if testYMin
     1084                testin=testin & (coord_Y >= YMin);
     1085                testbound=1;
     1086            end
     1087            if testYMin
     1088                testin=testin & (coord_Y <= YMax);
     1089                testbound=1;
     1090            end
     1091            if testbound
     1092                indcut=find(testin);
     1093                for ivar=VarIndex
     1094                    VarName=FieldData.ListVarName{ivar};
     1095                    eval(['FieldData.' VarName '=FieldData.' VarName '(indcut);'])
     1096                end
     1097                coord_X=coord_X(indcut);
     1098                coord_Y=coord_Y(indcut);
     1099                if length(ivar_Z)==1
     1100                    coord_Z=coord_Z(indcut);
     1101                end
    11041102            end
    11051103        end
    11061104        % different cases of projection
    1107         if isequal(ObjectData.ProjMode,'projection')
    1108             %the list of dimension
    1109             %ProjData.ListDimName=[ProjData.ListDimName FieldData.VarDimName(VarIndex(1))];%add the point index to the list of dimensions
    1110             %ProjData.DimValue=[ProjData.
    1111             %length(coord_X)];
    1112            
    1113             for ivar=VarIndex %transfer variables to the projection plane
    1114                 VarName=FieldData.ListVarName{ivar};
    1115                 if ivar==ivar_X %x coordinate
    1116                     eval(['ProjData.' VarName '=coord_X;'])
    1117                 elseif ivar==ivar_Y % y coordinate
    1118                     eval(['ProjData.' VarName '=coord_Y;'])
    1119                 elseif isempty(ivar_Z) || ivar~=ivar_Z % other variables (except Z coordinate wyhich is not reproduced)
    1120                     eval(['ProjData.' VarName '=FieldData.' VarName ';'])
    1121                 end
    1122                 if isempty(ivar_Z) || ivar~=ivar_Z
    1123                     ProjData.ListVarName=[ProjData.ListVarName VarName];
    1124                     ProjData.VarDimName=[ProjData.VarDimName DimCell];
    1125                     nbvar=nbvar+1;
    1126                     if isfield(FieldData,'VarAttribute') && length(FieldData.VarAttribute) >=ivar
    1127                         ProjData.VarAttribute{nbvar}=FieldData.VarAttribute{ivar};
    1128                     end
    1129                 end
    1130             end
    1131         elseif isequal(ObjectData.ProjMode,'interp')||isequal(ObjectData.ProjMode,'filter')%interpolate data on a regular grid
    1132             if isequal(ObjectData.ProjMode,'filter')
    1133                 rho=1000;%smoothing parameter, (small for strong smoothing)
    1134             else
    1135                 rho=0;
    1136             end
    1137             coord_x_proj=XMin:DX:XMax;
    1138             coord_y_proj=YMin:DY:YMax;
    1139             if isfield(FieldData,[VarName '_tps'])
    1140                 [XI,YI]=meshgrid(coord_x_proj,coord_y_proj');
    1141                 XI=reshape(XI,[],1);
    1142                 YI=reshape(YI,[],1);         
    1143             end
    1144             DimCell={'coord_y','coord_x'};
    1145             ProjData.ListVarName={'coord_y','coord_x'};
    1146             ProjData.VarDimName={'coord_y','coord_x'};
    1147             nbcoord=2;
    1148             ProjData.coord_y=[YMin YMax];
    1149             ProjData.coord_x=[XMin XMax];
    1150             if isempty(ivar_X), ivar_X=0; end;
    1151             if isempty(ivar_Y), ivar_Y=0; end;
    1152             if isempty(ivar_Z), ivar_Z=0; end;
    1153             if isempty(ivar_U), ivar_U=0; end;
    1154             if isempty(ivar_V), ivar_V=0; end;
    1155             if isempty(ivar_W), ivar_W=0; end;
    1156             if isempty(ivar_F), ivar_F=0; end;
    1157             if isempty(ivar_FF), ivar_FF=0; end;
    1158             if ~isequal(ivar_FF,0)
    1159                 VarName_FF=FieldData.ListVarName{ivar_FF};
    1160                 eval(['indsel=find(FieldData.' VarName_FF '==0);'])
    1161                 coord_X=coord_X(indsel);
    1162                 coord_Y=coord_Y(indsel);
    1163             end
    1164             FF=zeros(1,length(coord_y_proj)*length(coord_x_proj));
    1165             testFF=0;
    1166             FieldName
    1167             if ~isempty(FieldName)
    1168                 ProjData=calc_field(FieldName,FieldData,[XI YI]);
    1169             else
     1105        switch ObjectData.ProjMode
     1106            case 'projection'
     1107                %the list of dimension
     1108                %ProjData.ListDimName=[ProjData.ListDimName FieldData.VarDimName(VarIndex(1))];%add the point index to the list of dimensions
     1109                %ProjData.DimValue=[ProjData.
     1110                %length(coord_X)];
     1111               
     1112                for ivar=VarIndex %transfer variables to the projection plane
     1113                    VarName=FieldData.ListVarName{ivar};
     1114                    if ivar==ivar_X %x coordinate
     1115                        eval(['ProjData.' VarName '=coord_X;'])
     1116                    elseif ivar==ivar_Y % y coordinate
     1117                        eval(['ProjData.' VarName '=coord_Y;'])
     1118                    elseif isempty(ivar_Z) || ivar~=ivar_Z % other variables (except Z coordinate wyhich is not reproduced)
     1119                        eval(['ProjData.' VarName '=FieldData.' VarName ';'])
     1120                    end
     1121                    if isempty(ivar_Z) || ivar~=ivar_Z
     1122                        ProjData.ListVarName=[ProjData.ListVarName VarName];
     1123                        ProjData.VarDimName=[ProjData.VarDimName DimCell];
     1124                        nbvar=nbvar+1;
     1125                        if isfield(FieldData,'VarAttribute') && length(FieldData.VarAttribute) >=ivar
     1126                            ProjData.VarAttribute{nbvar}=FieldData.VarAttribute{ivar};
     1127                        end
     1128                    end
     1129                end
     1130            case 'interp'%interpolate data on a regular grid
     1131                coord_x_proj=XMin:DX:XMax;
     1132                coord_y_proj=YMin:DY:YMax;
     1133                DimCell={'coord_y','coord_x'};
     1134                ProjData.ListVarName={'coord_y','coord_x'};
     1135                ProjData.VarDimName={'coord_y','coord_x'};
     1136                nbcoord=2;
     1137                ProjData.coord_y=[YMin YMax];
     1138                ProjData.coord_x=[XMin XMax];
     1139                if isempty(ivar_X), ivar_X=0; end;
     1140                if isempty(ivar_Y), ivar_Y=0; end;
     1141                if isempty(ivar_Z), ivar_Z=0; end;
     1142                if isempty(ivar_U), ivar_U=0; end;
     1143                if isempty(ivar_V), ivar_V=0; end;
     1144                if isempty(ivar_W), ivar_W=0; end;
     1145                if isempty(ivar_F), ivar_F=0; end;
     1146                if isempty(ivar_FF), ivar_FF=0; end;
     1147                if ~isequal(ivar_FF,0)
     1148                    VarName_FF=FieldData.ListVarName{ivar_FF};
     1149                    eval(['indsel=find(FieldData.' VarName_FF '==0);'])
     1150                    coord_X=coord_X(indsel);
     1151                    coord_Y=coord_Y(indsel);
     1152                end
     1153               
     1154                FF=zeros(1,length(coord_y_proj)*length(coord_x_proj));
     1155                testFF=0;
    11701156                for ivar=VarIndex
    11711157                    VarName=FieldData.ListVarName{ivar};
    … …  
    11801166                            FieldData.(VarName)=FieldData.(VarName)(indsel);
    11811167                        end
    1182                         %                     if isfield(FieldData,[VarName '_tps'])
    1183                         %                         [XI,YI]=meshgrid(coord_x_proj,coord_y_proj');
    1184                         %                         XI=reshape(XI,[],1);
    1185                         %                         YI=reshape(YI,[],1);
    1186                         %
    1187                         ProjData.(VarName)=griddata_uvmat(double(coord_X),double(coord_Y),double(FieldData.(VarName)),coord_x_proj,coord_y_proj',rho);
     1168                        ProjData.(VarName)=griddata_uvmat(double(coord_X),double(coord_Y),double(FieldData.(VarName)),coord_x_proj,coord_y_proj');
    11881169                        varline=reshape(ProjData.(VarName),1,length(coord_y_proj)*length(coord_x_proj));
    11891170                        FFlag= isnan(varline); %detect undefined values NaN
    … …  
    12121193                    ProjData.VarAttribute{ivar_new+1+nbcoord}.Role='errorflag';
    12131194                end
    1214             end
    1215         end
    1216        
     1195            case 'filter'
     1196                if ~isempty(VarType.coord_tps)
     1197                    VarType.coord_tps
     1198                    coord_x_proj=XMin:DX:XMax;
     1199                    coord_y_proj=YMin:DY:YMax;
     1200                    np_x=numel(coord_x_proj);
     1201                    np_y=numel(coord_y_proj);
     1202                    [XI,YI]=meshgrid(coord_x_proj,coord_y_proj');
     1203                    XI=reshape(XI,[],1);
     1204                    YI=reshape(YI,[],1);
     1205                    ProjData=calc_field(FieldData.FieldList,FieldData,[XI YI]);
     1206                    for ilist=3:length(ProjData.ListVarName)% reshape data, excluding coordinates (ilist=1-2), TODO: rationalise
     1207                        VarName=ProjData.ListVarName{ilist};
     1208                        ProjData.(VarName)=reshape(ProjData.(VarName),np_y,np_x);
     1209                    end
     1210                    ProjData.coord_x=[XMin XMax];
     1211                    ProjData.coord_y=[YMin YMax];
     1212                end
     1213        end
    12171214        %% case of input fields defined on a structured  grid
    12181215    else
    … …  
    13521349        % case with no  interpolation
    13531350        if isequal(ProjMode,'projection') && (~testangle || test90y || test90x)
    1354             if  NbDim==2 && ~testXMin && ~testXMax && ~testYMin && ~testYMax 
     1351            if  NbDim==2 && ~testXMin && ~testXMax && ~testYMin && ~testYMax
    13551352                ProjData=FieldData;% no change by projection
    13561353            else
    … …  
    13751372                    min_indx=ceil((Coord{NbDim}(1)-XMax)/DXinit)+1;
    13761373                    max_indx=floor((Coord{NbDim}(1)-XMin)/DXinit)+1;
    1377                     Xbound(2)=Coord{NbDim}(1)+DXinit*(max_indx-1);                         
     1374                    Xbound(2)=Coord{NbDim}(1)+DXinit*(max_indx-1);
    13781375                    Xbound(1)=Coord{NbDim}(1)+DXinit*(min_indx-1);
    13791376                end
    13801377                min_indy=max(min_indy,1);% deals with margin (bound lower than the first index)
    13811378                min_indx=max(min_indx,1);
    1382 
     1379               
    13831380                if test90y
    13841381                    ind_new=[3 2 1];
    13851382                    DimCell={AYProjName,AXProjName};
    1386 %                     DimValue=DimValue(ind_new);
     1383                    %                     DimValue=DimValue(ind_new);
    13871384                    iz=ceil((ObjectData.Coord(1,1)-Coord{3}(1))/DX)+1;
    13881385                    for ivar=VarIndex
    … …  
    16561653    ivar_F=VarType.warnflag;
    16571654    ivar_FF=VarType.errorflag;
    1658     testX=~isempty(ivar_X) && ~isempty(ivar_Y);
     1655    check_unstructured_coord=~isempty(ivar_X) && ~isempty(ivar_Y);
    16591656    DimCell=FieldData.VarDimName{VarIndex(1)};
    16601657    if ischar(DimCell)
    … …  
    16631660
    16641661%% case of input fields with unstructured coordinates
    1665     if testX
     1662    if check_unstructured_coord
    16661663        XName=FieldData.ListVarName{ivar_X};
    16671664        YName=FieldData.ListVarName{ivar_Y};
Note: See TracChangeset for help on using the changeset viewer.