Changeset 653 for trunk/src/proj_field.m


Ignore:
Timestamp:
Jun 26, 2013, 11:24:42 PM (7 years ago)
Author:
sommeria
Message:

bugs corrected for displaying small rgb color images.
proj_field/proj_plane imrpoved
browser changed to a more standard GUI presentaton

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/proj_field.m

    r651 r653  
    873873%-----------------------------------------------------------------
    874874
    875 %% rotation angles 
    876 PlaneAngle=[0 0 0]; 
     875%% rotation angles
     876PlaneAngle=[0 0 0];
    877877norm_plane=[0 0 1];
    878878cos_om=1;
     
    893893    norm_plane(3)=OmAxis(3)*coeff+cos_om;
    894894end
    895 testangle=~isequal(PlaneAngle,[0 0 0]);% && ~test90y && ~test90x;%=1 for slanted plane 
     895testangle=~isequal(PlaneAngle,[0 0 0]);% && ~test90y && ~test90x;%=1 for slanted plane
    896896
    897897%% mesh sizes DX and DY
     
    899899DY=[];%default
    900900if isfield(ObjectData,'DX') && ~isempty(ObjectData.DX)
    901      DX=abs(ObjectData.DX);%mesh of interpolation points
     901    DX=abs(ObjectData.DX);%mesh of interpolation points
    902902elseif isfield(FieldData,'CoordMesh')
    903      DX=FieldData.CoordMesh;
     903    DX=FieldData.CoordMesh;
    904904end
    905905if isfield(ObjectData,'DY') && ~isempty(ObjectData.DY)
    906      DY=abs(ObjectData.DY);%mesh of interpolation points
     906    DY=abs(ObjectData.DY);%mesh of interpolation points
    907907elseif isfield(FieldData,'CoordMesh')
    908      DY=FieldData.CoordMesh;
     908    DY=FieldData.CoordMesh;
    909909end
    910910if  ~strcmp(ObjectData.ProjMode,'projection') && (isempty(DX)||isempty(DY))
    911         errormsg='DX or DY not defined';
    912         return
     911    errormsg='DX or DY not defined';
     912    return
    913913end
    914914
     
    918918testYMin=0;
    919919testYMax=0;
    920 
    921920if isfield(ObjectData,'RangeX')
    922         XMin=min(ObjectData.RangeX);
    923         XMax=max(ObjectData.RangeX);
    924         testXMin=XMax>XMin;
    925         testXMax=1;% range restriction along X
    926 else
    927     XMin=FieldData.XMin;%default
    928 XMax=FieldData.XMax;%default
     921    XMin=min(ObjectData.RangeX);
     922    XMax=max(ObjectData.RangeX);
     923    testXMin=XMax>XMin;
     924    testXMax=1;% range restriction along X
     925    % else
     926    %     XMin=FieldData.XMin;%default
     927    %     XMax=FieldData.XMax;%default
    929928end
    930929if isfield(ObjectData,'RangeY')
    931         YMin=min(ObjectData.RangeY);
    932         YMax=max(ObjectData.RangeY);
    933         testYMin=YMax>YMin;
    934         testYMax=1;
    935 else
    936     YMin=FieldData.YMin;%default
    937 YMax=FieldData.YMax;%default
     930    YMin=min(ObjectData.RangeY);
     931    YMax=max(ObjectData.RangeY);
     932    testYMin=YMax>YMin;
     933    testYMax=1;
     934    % else
     935    %     YMin=FieldData.YMin;%default
     936    %     YMax=FieldData.YMax;%default
    938937end
    939938width=0;%default width of the projection band
    940939if isfield(ObjectData,'RangeZ')
    941         width=max(ObjectData.RangeZ);
     940    width=max(ObjectData.RangeZ);
    942941end
    943942
     
    950949%% reproduce initial plane position and angle
    951950if isfield(FieldData,'PlaneCoord')&&length(FieldData.PlaneCoord)==3&& isfield(ProjData,'ProjObjectCoord')
    952    if length(ProjData.ProjObjectCoord)==3% if the projection plane has a z coordinate
    953        if ~isequal(ProjData.PlaneCoord(3),ProjData.ProjObjectCoord) %check the consistency with the z coordinate of the field plane (set by calibration)
    954            errormsg='inconsistent z position for field and projection plane';
    955            return
    956        end
    957    else % the z coordinate is set only by the field plane (by calibration)
    958        ProjData.ProjObjectCoord(3)=FieldData.PlaneCoord(3);
    959    end
    960    if isfield(FieldData,'PlaneAngle')
    961        if isfield(ProjData,'ProjObjectAngle')
    962            if ~isequal(FieldData.PlaneAngle,ProjData.ProjObjectAngle) %check the consistency with the z coordinate of the field plane (set by calibration)
    963            errormsg='inconsistent plane angle for field and projection plane';
    964            return
    965            end
    966        else
    967         ProjData.ProjObjectAngle=FieldData.PlaneAngle;
    968        end
     951    if length(ProjData.ProjObjectCoord)==3% if the projection plane has a z coordinate
     952        if ~isequal(ProjData.PlaneCoord(3),ProjData.ProjObjectCoord) %check the consistency with the z coordinate of the field plane (set by calibration)
     953            errormsg='inconsistent z position for field and projection plane';
     954            return
     955        end
     956    else % the z coordinate is set only by the field plane (by calibration)
     957        ProjData.ProjObjectCoord(3)=FieldData.PlaneCoord(3);
     958    end
     959    if isfield(FieldData,'PlaneAngle')
     960        if isfield(ProjData,'ProjObjectAngle')
     961            if ~isequal(FieldData.PlaneAngle,ProjData.ProjObjectAngle) %check the consistency with the z coordinate of the field plane (set by calibration)
     962                errormsg='inconsistent plane angle for field and projection plane';
     963                return
     964            end
     965        else
     966            ProjData.ProjObjectAngle=FieldData.PlaneAngle;
     967        end
    969968    end
    970969end
     
    10091008
    10101009%% define the new coordinates in case of interpolation on a grid
    1011 if check_grid% TODO: recalculate coordinates to get the bounds in the rotated coordinates
     1010if check_grid
     1011    AYName='coord_y';
     1012    AXName='coord_x';
    10121013    ProjData.ListVarName={'coord_y','coord_x'};
    1013     ProjData.VarDimName={'coord_y','coord_x'}; 
     1014    ProjData.VarDimName={'coord_y','coord_x'};
    10141015    ProjData.VarAttribute={[],[]};
    1015     ProjData.coord_y=[YMin YMax];
     1016    if ~testYMin
     1017        errormsg='min Y value not defined for the projection grid';return
     1018    end
     1019    if ~testYMax
     1020        errormsg='max Y value not defined for the projection grid';return
     1021    end
     1022    if ~testXMin
     1023        errormsg='min X value not defined for the projection grid';return
     1024    end
     1025    if ~testXMax
     1026        errormsg='max X value not defined for the projection grid';return
     1027    end
     1028    ProjData.coord_y=[YMin YMax];%note that if projection is done on a grid, the Min and Max along each direction must have been defined
    10161029    ProjData.coord_x=[XMin XMax];
     1030    coord_x_proj=XMin:DX:XMax;
     1031    coord_y_proj=YMin:DY:YMax;
     1032    [X,YI]=meshgrid(coord_x_proj,coord_y_proj);%grid in the new coordinates
     1033    XI=ObjectData.Coord(1,1)+(X)*cos(PlaneAngle(3))-YI*sin(PlaneAngle(3));%corresponding coordinates in the original system
     1034    YI=ObjectData.Coord(1,2)+(X)*sin(PlaneAngle(3))+YI*cos(PlaneAngle(3));
    10171035end
    10181036
     
    10201038% LOOP ON GROUPS OF VARIABLES SHARING THE SAME DIMENSIONS
    10211039% CellVarIndex=cells of variable index arrays
    1022 ivar_new=0; % index of the current variable in the projected field
     1040%ivar_new=0; % index of the current variable in the projected field
    10231041% icoord=0;
    10241042nbcoord=0;%number of added coordinate variables brought by projection
    1025 nbvar=0;
     1043%nbvar=0;
    10261044vector_x_proj=[];
    10271045vector_y_proj=[];
     
    10321050    end
    10331051    VarIndex=CellInfo{icell}.VarIndex;%  indices of the selected variables in the list FieldData.ListVarName
    1034     ivar_U=[];ivar_V=[];ivar_W=[];
    1035     if isfield(CellInfo{icell},'VarIndex_vector_x_tps')&&isfield(CellInfo{icell},'VarIndex_vector_y_tps')
    1036         ivar_U=CellInfo{icell}.VarIndex_vector_x_tps;
    1037         ivar_V=CellInfo{icell}.VarIndex_vector_y_tps;
    1038     elseif isfield(CellInfo{icell},'VarIndex_vector_x')&&isfield(CellInfo{icell},'VarIndex_vector_y')
    1039         ivar_U=CellInfo{icell}.VarIndex_vector_x;
    1040         ivar_V=CellInfo{icell}.VarIndex_vector_y;
    1041     end
    1042     if isfield(CellInfo{icell},'VarIndex_vector_z')
    1043         ivar_W=CellInfo{icell}.VarIndex_vector_z;
    1044     end
     1052%     ivar_U=[];ivar_V=[];ivar_W=[];
     1053%     if isfield(CellInfo{icell},'VarIndex_vector_x_tps')&&isfield(CellInfo{icell},'VarIndex_vector_y_tps')
     1054%         ivar_U=CellInfo{icell}.VarIndex_vector_x_tps;
     1055%         ivar_V=CellInfo{icell}.VarIndex_vector_y_tps;
     1056%     elseif isfield(CellInfo{icell},'VarIndex_vector_x')&&isfield(CellInfo{icell},'VarIndex_vector_y')
     1057%         ivar_U=CellInfo{icell}.VarIndex_vector_x;
     1058%         ivar_V=CellInfo{icell}.VarIndex_vector_y;
     1059%     end
     1060%     if isfield(CellInfo{icell},'VarIndex_vector_z')
     1061%         ivar_W=CellInfo{icell}.VarIndex_vector_z;
     1062%     end
    10451063    %dimensions
    10461064    DimCell=FieldData.VarDimName{VarIndex(1)};
     
    10491067    end
    10501068    coord_z=0;%default
    1051 
     1069    ListVarName={};% initiate list of projected variables for cell # icell
     1070    VarDimName={};% initiate coresponding list of dimensions for cell # icell
     1071    VarAttribute={};% initiate coresponding list of var attributes  for cell # icell
    10521072    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    10531073    switch CellInfo{icell}.CoordType
    10541074       
    10551075        case 'scattered'
    1056         %% case of input fields with unstructured coordinates (applies for projMode ='projection' or 'interp_lin')
     1076            %% case of input fields with unstructured coordinates (applies for projMode ='projection' or 'interp_lin')
    10571077            if strcmp(ProjMode{icell},'interp_tps')
    10581078                continue %skip for next cell (needs tps field cell)
    10591079            end
    1060             coord_x=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(end)});% initial x coordinates 
     1080            coord_x=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(end)});% initial x coordinates
    10611081            coord_y=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(end-1)});% initial y coordinates
    10621082            check3D=(numel(CellInfo{icell}.CoordIndex)==3);
     
    10941114                coord_Y=(-coord_x *sin(Phi) + coord_y *cos(Phi))*cos(Theta);
    10951115                coord_Y=coord_Y+coord_z *sin(Theta);
    1096                 coord_X=(coord_X *cos(Psi) - coord_Y* sin(Psi));%A VERIFIER               
     1116                coord_X=(coord_X *cos(Psi) - coord_Y* sin(Psi));%A VERIFIER
    10971117                coord_Y=(coord_X *sin(Psi) + coord_Y* cos(Psi));
    10981118            else
     
    11371157            end
    11381158           
    1139             % two cases of projection
     1159            % two cases of projection for scattered coordinates
    11401160            switch ProjMode{icell}
    1141                 case 'projection' 
    1142                     nbvar=numel(ProjData.ListVarName);
     1161                case 'projection'
     1162                    nbvar=0;
     1163                    %nbvar=numel(ProjData.ListVarName);
    11431164                    for ivar=VarIndex %transfer variables to the projection plane
    11441165                        VarName=FieldData.ListVarName{ivar};
     
    11511172                        end
    11521173                        if ~(check3D && ivar==CellInfo{icell}.CoordIndex(1))
    1153                             ProjData.ListVarName=[ProjData.ListVarName VarName];
    1154                             ProjData.VarDimName=[ProjData.VarDimName DimCell];
     1174                            ListVarName=[ListVarName VarName];
     1175                            VarDimName=[VarDimName DimCell];
    11551176                            nbvar=nbvar+1;
    11561177                            if isfield(FieldData,'VarAttribute') && length(FieldData.VarAttribute) >=ivar
    1157                                 ProjData.VarAttribute{nbvar}=FieldData.VarAttribute{ivar};
     1178                                VarAttribute{nbvar}=FieldData.VarAttribute{ivar};
    11581179                            end
    11591180                        end
    11601181                    end
    11611182                case 'interp_lin'%interpolate data on a regular grid
    1162                     coord_x_proj=XMin:DX:XMax;
    1163                     coord_y_proj=YMin:DY:YMax;
    1164                     [XI,YI]=meshgrid(coord_x_proj,coord_y_proj);
    11651183                    if isfield(CellInfo{icell},'VarIndex_errorflag')
    11661184                        VarName_FF=FieldData.ListVarName{CellInfo{icell}.VarIndex_errorflag};
     
    11741192                    end
    11751193                    % interpolate and calculate field on the grid
    1176                     [VarVal,ListFieldProj,VarAttribute,errormsg]=calc_field_interp([coord_X coord_Y],FieldData,CellInfo{icell}.FieldName,XI,YI);
     1194                    [VarVal,ListVarName,VarAttribute,errormsg]=calc_field_interp([coord_X coord_Y],FieldData,CellInfo{icell}.FieldName,XI,YI);
    11771195                   
    11781196                    if isfield(CellInfo{icell},'CheckSub') && CellInfo{icell}.CheckSub && ~isempty(vector_x_proj)
    1179                         ProjData.(ProjData.ListVarName{vector_x_proj})=ProjData.(ProjData.ListVarName{vector_x_proj})-VarVal{1};
    1180                         ProjData.(ProjData.ListVarName{vector_y_proj})=ProjData.(ProjData.ListVarName{vector_y_proj})-VarVal{2};
     1197                        ProjData.(ListVarName{vector_x_proj})=ProjData.(ListVarName{vector_x_proj})-VarVal{1};
     1198                        ProjData.(ListVarName{vector_y_proj})=ProjData.(ListVarName{vector_y_proj})-VarVal{2};
    11811199                    else
    1182                         VarDimName=cell(size(ListFieldProj));
    1183                         for ilist=1:numel(ListFieldProj)% reshape data, excluding coordinates (ilist=1-2), TODO: rationalise
    1184                             ListFieldProj{ilist}=regexprep(ListFieldProj{ilist},'(.+','');
    1185                             if ~isempty(find(strcmp(ListFieldProj{ilist},ProjData.ListVarName)))
    1186                                 ListFieldProj{ilist}=[ListFieldProj{ilist} '_1'];
    1187                             end                       
    1188                             ProjData.(ListFieldProj{ilist})=VarVal{ilist};
     1200                        VarDimName=cell(size(ListVarName));
     1201                        for ilist=1:numel(ListVarName)% reshape data, excluding coordinates (ilist=1-2), TODO: rationalise
     1202                            ListVarName{ilist}=regexprep(ListVarName{ilist},'(.+','');
     1203                            if ~isempty(find(strcmp(ListVarName{ilist},ProjData.ListVarName)))
     1204                                ListVarName{ilist}=[ListVarName{ilist} '_1'];
     1205                            end
     1206                            ProjData.(ListVarName{ilist})=VarVal{ilist};
    11891207                            VarDimName{ilist}={'coord_y','coord_x'};
    11901208                        end
    1191                         ProjData.ListVarName=[ProjData.ListVarName ListFieldProj];
    1192                         ProjData.VarDimName=[ProjData.VarDimName VarDimName];
    1193                         ProjData.VarAttribute=[ProjData.VarAttribute VarAttribute];
    1194                     end
    1195             end
    1196 
     1209                    end
     1210            end
     1211           
    11971212        case 'tps'
    1198         %% case of tps data (applies only in interp_tps mode)
     1213            %% case of tps data (applies only in interp_tps mode)
    11991214            if strcmp(ProjMode{icell},'interp_tps')
    12001215                Coord=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex});
     
    12041219                    FieldVar=cat(3,FieldData.(FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_x_tps}),FieldData.(FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_y_tps}));
    12051220                end
    1206                 coord_x_proj=XMin:DX:XMax;
    1207                 coord_y_proj=YMin:DY:YMax;
    1208                 np_x=numel(coord_x_proj);
    1209                 np_y=numel(coord_y_proj);
    1210                 [XI,YI]=meshgrid(coord_x_proj,coord_y_proj');
    1211                 XI=XI+ObjectData.Coord(1,1);
    1212                 YI=YI+ObjectData.Coord(1,2);
    1213                 [DataOut,VarAttribute,errormsg]=calc_field_tps(Coord,NbCentres,SubRange,FieldVar,CellInfo{icell}.FieldName,cat(3,XI,YI));   
    1214                 ListFieldProj=(fieldnames(DataOut))';
    1215                 VarDimName=cell(size(ListFieldProj));
    1216                 for ilist=1:numel(ListFieldProj)% reshape data, excluding coordinates (ilist=1-2), TODO: rationalise
    1217                     VarName=ListFieldProj{ilist};
     1221                [DataOut,VarAttribute,errormsg]=calc_field_tps(Coord,NbCentres,SubRange,FieldVar,CellInfo{icell}.FieldName,cat(3,XI,YI));
     1222                ListVarName=(fieldnames(DataOut))';
     1223                VarDimName=cell(size(ListVarName));
     1224                for ilist=1:numel(ListVarName)% reshape data, excluding coordinates (ilist=1-2), TODO: rationalise
     1225                    VarName=ListVarName{ilist};
    12181226                    ProjData.(VarName)=DataOut.(VarName);
    12191227                    VarDimName{ilist}={'coord_y','coord_x'};
    12201228                end
    1221                 ProjData.ListVarName=[ProjData.ListVarName ListFieldProj];
    1222                 ProjData.VarDimName=[ProjData.VarDimName VarDimName];
    1223                 ProjData.VarAttribute=[ProjData.VarAttribute VarAttribute];
    12241229            end
    12251230           
     1231        case 'grid'
    12261232            %% case of input fields defined on a structured  grid
    1227         case 'grid'         
    12281233            VarName=FieldData.ListVarName{VarIndex(1)};%get the first variable of the cell to get the input matrix dimensions
    12291234            DimValue=size(FieldData.(VarName));%input matrix dimensions
     
    12431248                end
    12441249            end
    1245             AYName=FieldData.ListVarName{CellInfo{icell}.CoordIndex(NbDim-1)};%name of input x coordinate (name preserved on projection)
    1246             AXName=FieldData.ListVarName{CellInfo{icell}.CoordIndex(NbDim)};%name of input y coordinate (name preserved on projection)
    1247             if testangle% TODO modify name also in case of origin shift in x or y
    1248                 AYProjName='Y';
    1249                 AXProjName='X';
    1250                 count=0;
    1251                 %modify coordinate names if they are already used
    1252                 while ~(isempty(find(strcmp('AXName',ProjData.ListVarName),1)) && isempty(find(strcmp('AYName',ProjData.ListVarName),1)))
    1253                     count=count+1;
    1254                     AYProjName=[AYProjName '_' num2str(count)];
    1255                     AXProjName=[AXProjName '_' num2str(count)];
    1256                 end
    1257             else
    1258                 AYProjName=AYName;% (name preserved on projection)
    1259                 AXProjName=AXName;%name of input y coordinate (name preserved on projection)
    1260             end
    1261             ListDimName=FieldData.VarDimName{VarIndex(1)};
    1262             ProjData.ListVarName=[ProjData.ListVarName {AYProjName} {AXProjName}]; %TODO: check if it already exists in Projdata (several cells)
    1263             ProjData.VarDimName=[ProjData.VarDimName {AYProjName} {AXProjName}];
    1264             ProjData.VarAttribute=[ProjData.VarAttribute {[]} {[]}];
    12651250            Coord_z=[];
    12661251            Coord_y=[];
    1267             Coord_x=[];     
     1252            Coord_x=[];
     1253           
     1254            % find default mesh
    12681255            for idim=1:NbDim %loop on space dimensions
    12691256                test_interp(idim)=0;%test for coordiate interpolation (non regular grid), =0 by default
     
    12731260                    if numel(Coord{idim})==2 %input array defined on a regular grid
    12741261                        DCoord_min(idim)=(Coord{idim}(2)-Coord{idim}(1))/DimValue(idim);
     1262                        Coord{idim}=linspace(Coord{idim}(1),Coord{idim}(2),DimValue(idim));
    12751263                    else
    12761264                        DCoord=diff(Coord{idim});%array of coordinate derivatives for the input field
     
    13101298            DAX=DCoord_min(NbDim);
    13111299            DAY=DCoord_min(NbDim-1);
     1300           
     1301            % default bounds
    13121302            minAX=min(Coord_x);
    13131303            maxAX=max(Coord_x);
    13141304            minAY=min(Coord_y);
    13151305            maxAY=max(Coord_y);
    1316             xcorner=[minAX maxAX minAX maxAX]-ObjectData.Coord(1,1);
     1306            xcorner=[minAX maxAX minAX maxAX]-ObjectData.Coord(1,1);% image corners in the new coordiantes
    13171307            ycorner=[maxAY maxAY minAY minAY]-ObjectData.Coord(1,2);
    13181308            xcor_new=xcorner*cos_om+ycorner*sin_om;%coord new frame
     
    13461336                test_direct_z=test_direct(1);
    13471337            end
    1348             npX=floor((XMax-XMin)/DX+1);
    1349             npY=floor((YMax-YMin)/DY+1);
    1350             if test_direct_y
    1351                 coord_y_proj=linspace(YMin,YMax,npY);%abscissa of the new pixels along the line
    1352             else
    1353                 coord_y_proj=linspace(YMax,YMin,npY);%abscissa of the new pixels along the line
    1354             end
    1355             if test_direct_x
    1356                 coord_x_proj=linspace(XMin,XMax,npX);%abscissa of the new pixels along the line
    1357             else
    1358                 coord_x_proj=linspace(XMax,XMin,npX);%abscissa of the new pixels along the line
    1359             end
    1360             % case with no interpolation
    1361             if isequal(ProjMode{icell},'projection') && (~testangle || test90y || test90x)
     1338            %define new coordiantes if not already done by the definition of a projection grid
     1339            if ~check_grid
     1340                AYName=FieldData.ListVarName{CellInfo{icell}.CoordIndex(NbDim-1)};%name of input x coordinate (name preserved on projection)
     1341                AXName=FieldData.ListVarName{CellInfo{icell}.CoordIndex(NbDim)};%name of input y coordinate (name preserved on projection)
     1342                ListVarName=[ListVarName {AYName} {AXName}];% update the list of projected variables
     1343                VarDimName=[VarDimName {AYName} {AXName}];% update the corresponding list of dimensions
     1344                VarAttribute=[VarAttribute {[]} {[]}];% update the list of attributes
     1345                ProjData.(AYName)=[YMin YMax]; % new (projected ) y coordinates
     1346                ProjData.(AXName)=[XMin XMax]; % new (projected ) y coordinates
     1347            end
     1348           
     1349            if testangle
     1350                ProjMode{icell}='interp_lin'; %request linear interpolation for projection on a tilted plane
     1351            end
     1352           
     1353            %             npX=floor((XMax-XMin)/DX+1);
     1354            %             npY=floor((YMax-YMin)/DY+1);
     1355            %             if test_direct_y
     1356            %                 coord_y_proj=linspace(YMin,YMax,npY);%abscissa of the new pixels along the line
     1357            %             else
     1358            %                 coord_y_proj=linspace(YMax,YMin,npY);%abscissa of the new pixels along the line
     1359            %             end
     1360            %             if test_direct_x
     1361            %                 coord_x_proj=linspace(XMin,XMax,npX);%abscissa of the new pixels along the line
     1362            %             else
     1363            %                 coord_x_proj=linspace(XMax,XMin,npX);%abscissa of the new pixels along the line
     1364            %             end
     1365            if isequal(ProjMode{icell},'projection')% && (~testangle || test90y || test90x)
    13621366                if  NbDim==2 && ~testXMin && ~testXMax && ~testYMin && ~testYMax% no range restriction
    1363                     ProjData.ListVarName=[ProjData.ListVarName FieldData.ListVarName(VarIndex)];
    1364                     ProjData.VarDimName=[ProjData.VarDimName FieldData.VarDimName(VarIndex)]; 
     1367                    ListVarName=[ListVarName FieldData.ListVarName(VarIndex)];
     1368                    VarDimName=[VarDimName FieldData.VarDimName(VarIndex)];
    13651369                    if isfield(FieldData,'VarAttribute')
    1366                     ProjData.VarAttribute=[ProjData.VarAttribute FieldData.VarAttribute(VarIndex)];
    1367                     end
    1368                     ProjData.(AYProjName)=FieldData.(AYName);
    1369                     ProjData.(AXProjName)=FieldData.(AXName);
     1370                        VarAttribute=[VarAttribute FieldData.VarAttribute(VarIndex)];
     1371                    end
     1372                    ProjData.(AYName)=FieldData.(AYName);
     1373                    ProjData.(AXName)=FieldData.(AXName);
    13701374                    for ivar=VarIndex
    13711375                        VarName=FieldData.ListVarName{ivar};
     
    13971401                    end
    13981402                    min_indy=max(min_indy,1);% deals with margin (bound lower than the first index)
    1399                     min_indx=max(min_indx,1);                 
     1403                    min_indx=max(min_indx,1);
    14001404                    if test90y
    14011405                        ind_new=[3 2 1];
     
    14041408                        for ivar=VarIndex
    14051409                            VarName=FieldData.ListVarName{ivar};
    1406                             ProjData.ListVarName=[ProjData.ListVarName VarName];
    1407                             ProjData.VarDimName=[ProjData.VarDimName {DimCell}];
    1408                             ProjData.VarAttribute{length(ProjData.ListVarName)}=FieldData.VarAttribute{ivar}; %reproduce the variable attributes
     1410                            ListVarName=[ListVarName VarName];
     1411                            VarDimName=[VarDimName {DimCell}];
     1412                            VarAttribute{length(ListVarName)}=FieldData.VarAttribute{ivar}; %reproduce the variable attributes
    14091413                            eval(['ProjData.' VarName '=permute(FieldData.' VarName ',ind_new);'])% permute x and z indices for 90 degree rotation
    14101414                            eval(['ProjData.' VarName '=squeeze(ProjData.' VarName '(iz,:,:));'])% select the z index iz
    14111415                        end
    1412                         eval(['ProjData.' AYProjName '=[Ybound(1) Ybound(2)];']) %record the new (projected ) y coordinates
    1413                         eval(['ProjData.' AXProjName '=[Coord{1}(end),Coord{1}(1)];']) %record the new (projected ) x coordinates
     1416                        eval(['ProjData.' AYName '=[Ybound(1) Ybound(2)];']) %record the new (projected ) y coordinates
     1417                        eval(['ProjData.' AXName '=[Coord{1}(end),Coord{1}(1)];']) %record the new (projected ) x coordinates
    14141418                    else
    14151419                        if NbDim==3
     
    14261430                        for ivar=VarIndex% loop on non coordinate variables
    14271431                            VarName=FieldData.ListVarName{ivar};
    1428                             ProjData.ListVarName=[ProjData.ListVarName VarName];
    1429                             ProjData.VarDimName=[ProjData.VarDimName {DimCell}];
     1432                            ListVarName=[ListVarName VarName];
     1433                            VarDimName=[VarDimName {DimCell}];
    14301434                            if isfield(FieldData,'VarAttribute') && length(FieldData.VarAttribute)>=ivar
    1431                                 ProjData.VarAttribute{length(ProjData.ListVarName)}=FieldData.VarAttribute{ivar};
     1435                                VarAttribute{length(ListVarName)}=FieldData.VarAttribute{ivar};
    14321436                            end
    14331437                            if NbDim==3
     
    14371441                            end
    14381442                        end
    1439                         eval(['ProjData.' AYProjName '=[Ybound(1) Ybound(2)];']) %record the new (projected ) y coordinates
    1440                         eval(['ProjData.' AXProjName '=[Xbound(1) Xbound(2)];']) %record the new (projected ) x coordinates
    1441                     end
    1442                 end
    1443             else       % case with rotation and/or interpolation
     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
     1445                    end
     1446                end
     1447            else       % case with interpolation
    14441448                if NbDim==2 %2D case
    1445                     [X,Y]=meshgrid(coord_x_proj,coord_y_proj);%grid in the new coordinates
    1446                     XIMA=ObjectData.Coord(1,1)+(X)*cos(PlaneAngle(3))-Y*sin(PlaneAngle(3));%corresponding coordinates in the original image
    1447                     YIMA=ObjectData.Coord(1,2)+(X)*sin(PlaneAngle(3))+Y*cos(PlaneAngle(3));
    1448                     XIMA=(XIMA-minAX)/DXinit+1;% image index along x
    1449                     YIMA=(-YIMA+maxAY)/DYinit+1;% image index along y
    1450                     XIMA=reshape(round(XIMA),1,npX*npY);%indices reorganized in 'line'
    1451                     YIMA=reshape(round(YIMA),1,npX*npY);
    1452                     flagin=XIMA>=1 & XIMA<=DimValue(2) & YIMA >=1 & YIMA<=DimValue(1);%flagin=1 inside the original image
    14531449                    if isequal(ProjMode{icell},'interp_tps')
    14541450                        npx_interp_tps=ceil(abs(DX/DAX));
     
    14591455                        test_interp_tps=0;
    14601456                    end
    1461                     eval(['ProjData.' AYName '=[coord_y_proj(1) coord_y_proj(end)];']) %record the new (projected ) y coordinates
    1462                     eval(['ProjData.' AXName '=[coord_x_proj(1) coord_x_proj(end)];']) %record the new (projected ) x coordinates
     1457                    %                     ProjData.coord_y=[YMin YMax];%note that if projection is done on a grid, the Min and Max along each direction must have been defined
     1458                    %                     ProjData.coord_x=[XMin XMax];
     1459                    coord_x_proj=XMin:DX:XMax;
     1460                    coord_y_proj=YMin:DY:YMax;
     1461                    [X,YI]=meshgrid(coord_x_proj,coord_y_proj);%grid in the new coordinates
     1462                    XI=ObjectData.Coord(1,1)+(X)*cos(PlaneAngle(3))-YI*sin(PlaneAngle(3));%corresponding coordinates in the original system
     1463                    YI=ObjectData.Coord(1,2)+(X)*sin(PlaneAngle(3))+YI*cos(PlaneAngle(3));
     1464                    [X,Y]=meshgrid(Coord{2},Coord{1});%initial coordiantes
    14631465                    for ivar=VarIndex
    14641466                        VarName=FieldData.ListVarName{ivar};
    1465                         if test_interp(1) || test_interp(2)%interpolate on a regular grid
    1466                             eval(['ProjData.' VarName '=interp2(Coord{2},Coord{1},FieldData.' VarName ',Coord_x,Coord_y'');']) %TO TEST
     1467                        %                         if test_interp(1) || test_interp(2)%interpolate on a regular grid
     1468                        if size(FieldData.(VarName),3)==1
     1469                            ProjData.(VarName)=interp2(X,Y,double(FieldData.(VarName)),XI,YI,'*linear');
     1470                        else
     1471                            ProjData.(VarName)=interp2(X,Y,double(FieldData.(VarName)(:,:,1)),XI,YI,'*linear');
     1472                            for icolor=2:size(FieldData.(VarName),3)
     1473                                ProjData.(VarName)=cat(3,ProjData.(VarName),interp2(X,Y,double(FieldData.(VarName)(:,:,icolor)),XI,YI,'*linear')); %TO TEST
     1474                            end
    14671475                        end
    1468                         eval(['vec_A=reshape(FieldData.' VarName ',[],nbcolor);'])%put the original image in line
    1469                         %ind_in=find(flagin);
    1470                         ind_out=find(~flagin);
    1471                         ICOMB=(XIMA-1)*DimValue(1)+YIMA;
    1472                         ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A
    1473                         vec_B(flagin,1:nbcolor)=vec_A(ICOMB,:);
    1474                         for icolor=1:nbcolor
    1475                             vec_B(ind_out,icolor)=zeros(size(ind_out));
    1476                         end
    1477                         ProjData.ListVarName=[ProjData.ListVarName VarName];
    1478                         ProjData.VarDimName=[ProjData.VarDimName {DimCell}];
     1476                        ListVarName=[ListVarName VarName];
     1477                        DimCell(1:2)={AYName,AXName};
     1478                        VarDimName=[VarDimName {DimCell}];
    14791479                        if isfield(FieldData,'VarAttribute')&&length(FieldData.VarAttribute)>=ivar
    1480                             ProjData.VarAttribute{length(ProjData.ListVarName)+nbcoord}=FieldData.VarAttribute{ivar};
    1481                         end
    1482                         eval(['ProjData.' VarName '=reshape(vec_B,npY,npX,nbcolor);']);
    1483                     end
    1484                     ProjData.FF=reshape(~flagin,npY,npX);%false flag A FAIRE: tenir compte d'un flga antï¿œrieur
    1485                     ProjData.ListVarName=[ProjData.ListVarName 'FF'];
    1486                     ProjData.VarDimName=[ProjData.VarDimName {DimCell}];
    1487                     ProjData.VarAttribute{length(ProjData.ListVarName)}.Role='errorflag';
     1480                            VarAttribute{length(ListVarName)+nbcoord}=FieldData.VarAttribute{ivar};
     1481                        end;
     1482                    end
    14881483                elseif ~testangle
    14891484                    % unstructured z coordinate
     
    14961491                        for ivar=VarIndex
    14971492                            VarName=FieldData.ListVarName{ivar};
    1498                             ProjData.ListVarName=[ProjData.ListVarName VarName];
    1499                             ProjData.VarAttribute{length(ProjData.ListVarName)}=FieldData.VarAttribute{ivar}; %reproduce the variable attributes
     1493                            ListVarName=[ListVarName VarName];
     1494                            VarAttribute{length(ListVarName)}=FieldData.VarAttribute{ivar}; %reproduce the variable attributes
    15001495                            eval(['ProjData.' VarName '=squeeze(FieldData.' VarName '(iz,:,:));'])% select the z index iz
    15011496                            %TODO : do a vertical average for a thick plane
     
    15121507            end
    15131508    end
     1509    % update the global list of projected variables:
     1510    ProjData.ListVarName=[ProjData.ListVarName ListVarName];
     1511    ProjData.VarDimName=[ProjData.VarDimName VarDimName];
     1512    ProjData.VarAttribute=[ProjData.VarAttribute VarAttribute];
    15141513   
    15151514    %% projection of  velocity components in the rotated coordinates
    1516     if testangle && length(ivar_U)==1
    1517         if isempty(ivar_V)
    1518             msgbox_uvmat('ERROR','v velocity component missing in proj_field.m')
    1519             return
    1520         end
    1521         UName=FieldData.ListVarName{ivar_U};
    1522         VName=FieldData.ListVarName{ivar_V};
    1523         eval(['ProjData.' UName  '=cos(PlaneAngle(3))*ProjData.' UName '+ sin(PlaneAngle(3))*ProjData.' VName ';'])
    1524         eval(['ProjData.' VName  '=cos(Theta)*(-sin(PlaneAngle(3))*ProjData.' UName '+ cos(PlaneAngle(3))*ProjData.' VName ');'])
    1525         if ~isempty(ivar_W)
    1526             WName=FieldData.ListVarName{ivar_W};
    1527             eval(['ProjData.' VName '=ProjData.' VName '+ ProjData.' WName '*sin(Theta);'])%
    1528             eval(['ProjData.' WName '=NormVec_X*ProjData.' UName '+ NormVec_Y*ProjData.' VName '+ NormVec_Z* ProjData.' WName ';']);
    1529         end
    1530         if ~isequal(Psi,0)
    1531             eval(['ProjData.' UName '=cos(Psi)* ProjData.' UName '- sin(Psi)*ProjData.' VName ';']);
    1532             eval(['ProjData.' VName '=sin(Psi)* ProjData.' UName '+ cos(Psi)*ProjData.' VName ';']);
    1533         end
    1534     end
    1535 end
    1536 
     1515    if testangle
     1516        ivar_U=[];ivar_V=[];ivar_W=[];
     1517        for ivar=1:numel(VarAttribute)
     1518            if isfield(VarAttribute{ivar},'Role')
     1519                if strcmp(VarAttribute{ivar}.Role,'vector_x')
     1520                    ivar_U=ivar;
     1521                elseif strcmp(VarAttribute{ivar}.Role,'vector_y')
     1522                    ivar_V=ivar;
     1523                elseif strcmp(VarAttribute{ivar}.Role,'vector_z')
     1524                    ivar_W=ivar;
     1525                end
     1526            end
     1527        end
     1528        if ~isempty(ivar_U)
     1529            if isempty(ivar_V)
     1530                msgbox_uvmat('ERROR','v velocity component missing in proj_field.m')
     1531                return
     1532            else
     1533                UName=ListVarName{ivar_U};
     1534                VName=ListVarName{ivar_V};
     1535                ProjData.(UName)=cos(PlaneAngle(3))*ProjData.(UName)+ sin(PlaneAngle(3))*ProjData.(VName);
     1536                ProjData.(VName)=(-sin(PlaneAngle(3))*ProjData.(UName)+ cos(PlaneAngle(3))*ProjData.(VName));
     1537                if ~isempty(ivar_W)
     1538                    WName=FieldData.ListVarName{ivar_W};
     1539                    eval(['ProjData.' VName '=ProjData.' VName '+ ProjData.' WName '*sin(Theta);'])%
     1540                    eval(['ProjData.' WName '=NormVec_X*ProjData.' UName '+ NormVec_Y*ProjData.' VName '+ NormVec_Z* ProjData.' WName ';']);
     1541                end
     1542%                 if ~isequal(Psi,0)
     1543%                     eval(['ProjData.' UName '=cos(Psi)* ProjData.' UName '- sin(Psi)*ProjData.' VName ';']);
     1544%                     eval(['ProjData.' VName '=sin(Psi)* ProjData.' UName '+ cos(Psi)*ProjData.' VName ';']);
     1545%                 end
     1546            end
     1547        end
     1548    end
     1549end
    15371550% %prepare substraction in case of two input fields
    15381551% SubData.ListVarName={};
Note: See TracChangeset for help on using the changeset viewer.