Changeset 648 for trunk/src/proj_field.m


Ignore:
Timestamp:
Jun 9, 2013, 10:31:58 PM (11 years ago)
Author:
sommeria
Message:

get_field updated, several bugs corrected,open_uvmat suppressd

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/proj_field.m

    r646 r648  
    8484ProjData=[];
    8585
    86 %% in the absence of object Type or projection mode, or object coordinaes, the output is empty
     86%% check input projection object: type, projection mode and Coord:
    8787if ~isfield(ObjectData,'Type')||~isfield(ObjectData,'ProjMode')
    8888    return
    8989end
    90 % case of no projection (object is used only as graph display)
    91 if isequal(ObjectData.ProjMode,'none')||isequal(ObjectData.ProjMode,'mask_inside')||isequal(ObjectData.ProjMode,'mask_outside')
     90ListProjMode={'projection','interp_lin','interp_tps'};%list of effective projection modes
     91if isempty(strcmp(ObjectData.ProjMode,ListProjMode))
    9292    return
    9393end
    … …  
    103103switch ObjectData.Type
    104104    case 'points'
    105     [ProjData,errormsg]=proj_points(FieldData,ObjectData);
     105        [ProjData,errormsg]=proj_points(FieldData,ObjectData);
    106106    case {'line','polyline'}
    107      [ProjData,errormsg] = proj_line(FieldData,ObjectData);
     107        [ProjData,errormsg] = proj_line(FieldData,ObjectData);
    108108    case {'polygon','rectangle','ellipse'}
    109109        if isequal(ObjectData.ProjMode,'inside')||isequal(ObjectData.ProjMode,'outside')
    … …  
    113113        end
    114114    case 'plane'
    115             [ProjData,errormsg] = proj_plane(FieldData,ObjectData);
     115        [ProjData,errormsg] = proj_plane(FieldData,ObjectData);
    116116    case 'volume'
    117117        [ProjData,errormsg] = proj_volume(FieldData,ObjectData);
    … …  
    344344%LOOP ON GROUPS OF VARIABLES SHARING THE SAME DIMENSIONS
    345345for icell=1:length(CellInfo)
    346     testX=0;
    347     testY=0;
     346    CoordType=CellInfo{icell}.CoordType;
     347%     testY=0;
    348348    test_Amat=0;
    349349    if NbDim(icell)~=2% proj_patch acts only on fields of space dimension 2
    … …  
    429429        end
    430430    elseif isequal(ObjectData.Type,'polygon')
    431         if testX
     431        if strcmp(CoordType,'scattered')
    432432            testin=inpolygon(coord_x,coord_y,ObjectData.Coord(:,1),ObjectData.Coord(:,2));
    433         elseif test_Amat
     433        elseif strcmp(CoordType,'grid')
    434434            testin=inpolygon(Xi,Yi,ObjectData.Coord(:,1),ObjectData.Coord(:,2));
    435435        else%calculate the scalar
    … …  
    439439        X2Max=widthx*widthx;
    440440        Y2Max=(widthy)*(widthy);
    441         if testX
     441        if strcmp(CoordType,'scattered')
    442442            distX=(coord_x-ObjectData.Coord(1,1));
    443443            distY=(coord_y-ObjectData.Coord(1,2));
    444444            testin=(distX.*distX/X2Max+distY.*distY/Y2Max)<1;
    445         elseif test_Amat %case of usual 2x2 matrix
     445        elseif strcmp(CoordType,'grid') %case of usual 2x2 matrix
    446446            distX=(Xi-ObjectData.Coord(1,1));
    447447            distY=(Yi-ObjectData.Coord(1,2));
    … …  
    495495ProjData.NbDim=1;
    496496%initialisation of the input parameters and defaultoutput
    497 ProjMode=ObjectData.ProjMode;
     497ProjMode=ObjectData.ProjMode; %rmq: ProjMode always defined from input={'projection','interp_lin','interp_tps'}
    498498% ProjAngle=90; %90 degrees projection by default
    499 
    500 width=0;%default width of the projection band
    501 if isfield(ObjectData,'Range')&&size(ObjectData.Range,2)>=2
    502     width=abs(ObjectData.Range(1,2));
    503 end
     499width=0;
    504500if isfield(ObjectData,'RangeY')
    505     width=max(ObjectData.RangeY);
    506 end
    507 
     501    width=max(ObjectData.RangeY);%Rangey needed bfor mode 'projection'
     502end
    508503% default output
    509 errormsg=[];%default
     504errormsg='';%default
    510505Xline=[];
    511506flux=0;
    512507circul=0;
    513508liny=ObjectData.Coord(:,2);
    514 siz_line=size(ObjectData.Coord);
    515 if siz_line(1)<2
    516     return% line needs at least 2 points to be defined
    517 end
     509NbPoints=size(ObjectData.Coord,1);
    518510testfalse=0;
    519511ListIndex={};
    520512
    521 %% angles of the polyline and boundaries of action
    522 dlinx=diff(ObjectData.Coord(:,1));
    523 dliny=diff(ObjectData.Coord(:,2));
    524 theta=angle(dlinx+1i*dliny);%angle of each segment
    525 theta(siz_line(1))=theta(siz_line(1)-1);
     513
     514%% projection line: object types selected from  proj_field='line','polyline','polygon','rectangle','ellipse':
     515LineCoord=ObjectData.Coord;
     516switch ObjectData.Type
     517    case 'ellipse'
     518        LineLength=2*pi*ObjectData.RangeX*ObjectData.RangeY;
     519        NbSegment=0;
     520    case 'rectangle'
     521        LineCoord([1 4],1)=ObjectData.Coord(1,1)-ObjectData.RangeX;
     522        LineCoord([1 2],2)=ObjectData.Coord(1,2)-ObjectData.RangeY;
     523        LineCoord([2 3],1)=ObjectData.Coord(1,1)+ObjectData.RangeX;
     524        LineCoord([4 1],2)=ObjectData.Coord(1,2)+ObjectData.RangeY;
     525    case 'polygon'
     526        LineCoord(NbPoints+1)=LineCoord(1);
     527end
     528if ~strcmp(ObjectData.Type,'ellipse')
     529    if ~strcmp(ObjectData.Type,'rectangle') && NbPoints<2
     530        return% line needs at least 2 points to be defined
     531    end
     532    dlinx=diff(LineCoord(:,1));
     533    dliny=diff(LineCoord(:,2));
     534    [theta,dlength]=cart2pol(dlinx,dliny);%angle and length of each segment
     535    LineLength=sum(dlength);
     536    NbSegment=numel(LineLength);
     537end
     538CheckClosedLine=~isempty(find(strcmp(ObjectData.Type,{'rectangle','ellipse','polygon'})));
     539
     540%     x = a \ \cosh \mu \ \cos \nu
     541%
     542%     y = a \ \sinh \mu \ \sin \nu
     543
     544%% angles of the polyline and boundaries of action for mode 'projection'
     545
    526546% determine a rectangles at +-width from the line (only used for the ProjMode='projection or 'interp_tps')
    527 xsup=zeros(1,siz_line(1)); xinf=zeros(1,siz_line(1)); ysup=zeros(1,siz_line(1)); yinf=zeros(1,siz_line(1));
    528 if isequal(ProjMode,'projection') || isequal(ProjMode,'interp_tps')
    529     xsup(1)=ObjectData.Coord(1,1)-width*sin(theta(1));
    530     xinf(1)=ObjectData.Coord(1,1)+width*sin(theta(1));
    531     ysup(1)=ObjectData.Coord(1,2)+width*cos(theta(1));
    532     yinf(1)=ObjectData.Coord(1,2)-width*cos(theta(1));
    533     for ip=2:siz_line(1)
    534         xsup(ip)=ObjectData.Coord(ip,1)-width*sin((theta(ip)+theta(ip-1))/2)/cos((theta(ip-1)-theta(ip))/2);
    535         xinf(ip)=ObjectData.Coord(ip,1)+width*sin((theta(ip)+theta(ip-1))/2)/cos((theta(ip-1)-theta(ip))/2);
    536         ysup(ip)=ObjectData.Coord(ip,2)+width*cos((theta(ip)+theta(ip-1))/2)/cos((theta(ip-1)-theta(ip))/2);
    537         yinf(ip)=ObjectData.Coord(ip,2)-width*cos((theta(ip)+theta(ip-1))/2)/cos((theta(ip-1)-theta(ip))/2);
    538     end
    539 end
    540 %
    541 %
    542 %     x = a \ \cosh \mu \ \cos \nu
    543 %
    544 %     y = a \ \sinh \mu \ \sin \nu
     547xsup=zeros(1,NbPoints); xinf=zeros(1,NbPoints); ysup=zeros(1,NbPoints); yinf=zeros(1,NbPoints);
     548if isequal(ProjMode,'projection')
     549    if strcmp(ObjectData.Type,'line')
     550        xsup=ObjectData.Coord(:,1)-width*sin(theta);
     551        xinf=ObjectData.Coord(:,1)+width*sin(theta);
     552        ysup=ObjectData.Coord(:,2)+width*cos(theta);
     553        yinf=ObjectData.Coord(:,2)-width*cos(theta);
     554    else
     555        errormsg='mode projection only available for simple line, use interpolation otherwise';
     556        return
     557    end
     558else % need to define the set of interpolation points
     559    if isfield(ObjectData,'DX') && ~isempty(ObjectData.DX)
     560        DX=abs(ObjectData.DX);%mesh of interpolation points along the line
     561        if CheckClosedLine
     562            NbPoint=ceil(LineLength/DX);
     563            DX=LineLength/NbPoint;%adjust DX to get an integer nbre of intervals in a closed line
     564            DX_end=DX/2;
     565        else
     566            DX_end=(LineLength-DX*floor(LineLength/DX))/2;%margin from the first point and first interpolation point
     567        end
     568        XI=[];
     569        YI=[];
     570        ThetaI=[];
     571        dlengthI=[];
     572        if strcmp(ObjectData.Type,'ellipse')
     573            phi=(DX_end:DX:LineLength)*2*pi/LineLength;
     574            XI=ObjectData.RangeX*cos(phi);
     575            YI=ObjectData.RangeY*sin(phi);
     576            dphi=2*pi*DX/LineLength;
     577            [ThetaI,dlengthI]=cart2pol(-ObjectData.RangeX*sin(phi)*dphi,ObjectData.RangeY*cos(phi)*dphi);
     578        else
     579            for isegment=1:NbSegment
     580                costheta=cos(theta(isegment));
     581                sintheta=sin(theta(isegment));
     582                XIsegment=(LineCoord(isegment,1)+DX_end*costheta:DX*costheta:LineCoord(isegment+1,1));
     583                YIsegment=(LineCoord(isegment,2)+DX_end*sintheta:DX*sintheta:LineCoord(isegment+1,2));
     584                XI=[XI XIsegment];
     585                YI=[YI YIsegment];
     586                ThetaI=[ThetaI theta(isegment)*ones(1,numel(XIsegment))];
     587                dlengthI=[dlengthI DX*ones(1,numel(XIsegment))];
     588                DX_end=DX_end+DX-(dlength(isegment)-DX*(numel(XIsegment)-1));
     589            end
     590        end
     591        Xproj=cumsum(dlengthI);
     592    else
     593        errormsg='mesh DX needed for interpolation';
     594        return
     595    end
     596end
    545597
    546598
    … …  
    551603    return
    552604end
    553 
    554 %% loop on variable cells with the same space dimension
     605CellInfo=CellInfo(NbDim==2); %keep only the 2D cells
     606%%%%%% TODO: treat 1D fields: project as identity so that P o P=P for projection operation
     607
     608%% loop on variable cells with the same space dimension 2
    555609ProjData.ListVarName={};
    556610ProjData.VarDimName={};
    557611for icell=1:length(CellInfo)
    558     if NbDim(icell)~=2% proj_line acts only on fields of space dimension 2, TODO: check 3D case
    559         continue
    560     end
    561 
    562     % select types of  variables to be projected
    563    ListProj={'VarIndex_scalar','VarIndex_image','VarIndex_color','VarIndex_vector_x','VarIndex_vector_y'};
    564    check_proj=false(size(FieldData.ListVarName));
    565    for ilist=1:numel(ListProj)
    566        if isfield(CellInfo{icell},ListProj{ilist})
    567            check_proj(CellInfo{icell}.(ListProj{ilist}))=1;
    568        end
    569    end
    570    VarIndex=find(check_proj);
    571 
    572     %% identify vector components   
     612    % list of variable types to be projected
     613    ListProj={'VarIndex_scalar','VarIndex_image','VarIndex_color','VarIndex_vector_x','VarIndex_vector_y'};
     614    check_proj=false(size(FieldData.ListVarName));
     615    for ilist=1:numel(ListProj)
     616        if isfield(CellInfo{icell},ListProj{ilist})
     617            check_proj(CellInfo{icell}.(ListProj{ilist}))=1;
     618        end
     619    end
     620    VarIndex=find(check_proj);% indices of the variables to be projected
     621   
     622    %% identify vector components
    573623    testU=isfield(CellInfo{icell},'VarIndex_vector_x') &&isfield(CellInfo{icell},'VarIndex_vector_y') ;% test for vectors
    574624    if testU
    … …  
    577627        vector_x=FieldData.(UName);
    578628        vector_y=FieldData.(VName);
    579     end 
     629    end
    580630    %identify error flag
    581     testfalse=isfield(CellInfo{icell},'VarIndex_errorflag');% test for error flag
    582     if testfalse
     631    errorflag=0; %default, no error flag
     632    if isfield(CellInfo{icell},'VarIndex_errorflag');% test for error flag
    583633        FFName=FieldData.ListVarName{CellInfo{icell}.VarIndex_errorflag};
    584634        errorflag=FieldData.(FFName);
    585     end   
     635    end
     636    VarName=FieldData.ListVarName(VarIndex);% cell array of the names of variables to pje
     637    %% check needed object properties for unstructured positions (position given by the variables with role coord_x, coord_y
    586638   
    587     %% check needed object properties for unstructured positions (position given by the variables with role coord_x, coord_y
    588     if strcmp(CellInfo{icell}.CoordType,'scattered')
    589         if  strcmp(ProjMode,'projection')
    590             if width==0
    591                 errormsg='range of the projection object is missing';
    592                 return 
    593             end
    594 %             else
    595 %                 lambda=2/(width*width); %smoothing factor used for interp_tps: weight exp(-2) at distance width from the line
    596 %             end
    597         else
    598             if isfield(ObjectData,'DX') && ~isempty(ObjectData.DX)
    599                 DX=abs(ObjectData.DX);%mesh of interpolation points along the line
     639    %         circul=0;
     640    %         flux=0;
     641    %%%%%%%  % A FAIRE CALCULER MEAN DES QUANTITES    %%%%%%
     642    switch CellInfo{icell}.CoordType
     643        %case of unstructured coordinates
     644        case 'scattered'
     645            XName= FieldData.ListVarName{CellInfo{icell}.CoordIndex(end)};
     646            YName= FieldData.ListVarName{CellInfo{icell}.CoordIndex(end-1)};
     647            coord_x=FieldData.(XName);
     648            coord_y=FieldData.(YName);
     649            if isequal(ProjMode,'projection')
     650                if width==0
     651                    errormsg='range of the projection object is missing';
     652                    return
     653                end
     654                % select the (non false) input data located in the band of projection
     655                flagsel=(errorflag==0) & ((coord_y -yinf(1))*(xinf(2)-xinf(1))>(coord_x-xinf(1))*(yinf(2)-yinf(1))) ...
     656                    & ((coord_y -ysup(1))*(xsup(2)-xsup(1))<(coord_x-xsup(1))*(ysup(2)-ysup(1))) ...
     657                    & ((coord_y -yinf(2))*(xsup(2)-xinf(2))>(coord_x-xinf(2))*(ysup(2)-yinf(2))) ...
     658                    & ((coord_y -yinf(1))*(xsup(1)-xinf(1))<(coord_x-xinf(1))*(ysup(1)-yinf(1)));
     659                coord_x=coord_x(flagsel);
     660                coord_y=coord_y(flagsel);
     661                costheta=cos(theta);
     662                sintheta=sin(theta);
     663                Xproj=(coord_x-ObjectData.Coord(1,1))*costheta + (coord_y-ObjectData.Coord(1,2))*sintheta; %projection on the line
     664                [Xproj,indsort]=sort(Xproj);% sort points by increasing absissa along the projection line
     665                for ivar=1:numel(VarIndex)
     666                    ProjData.(VarName{ivar})=FieldData.(VarName{ivar})(flagsel);% restrict vrtibles to the projection band
     667                    ProjData.(VarName{ivar})=ProjData.(VarName{ivar})(indsort);% sort by absissa
     668                end
     669                % project the velocity components if vectors are projected
     670                if testU
     671                    vector_x=ProjData.(UName);
     672                    ProjData.(UName)=costheta*vector_x+sintheta*ProjData.(VName);% longitudinal component
     673                    ProjData.(VName)=-sintheta*vector_x+costheta*ProjData.(VName);%transverse component
     674                end
     675            elseif isequal(ProjMode,'interp_lin')  %filtering %linear interpolation:
     676                [ProjVar,ListFieldProj,VarAttribute,errormsg]=calc_field_interp([coord_x coord_y],FieldData,CellInfo{icell}.FieldName,XI,YI);
     677                ProjData.X=Xproj;
     678                ProjData.ListVarName=[ProjData.ListVarName {XName}];
     679                ProjData.VarDimName=[ProjData.VarDimName {XName}];
     680                nbvar=numel(ProjData.ListVarName);
     681                ProjData.VarAttribute{nbvar}.long_name='abscissa along line';
     682                ProjData.ListVarName=[ProjData.ListVarName ListFieldProj];
     683                ProjData.VarAttribute=[ProjData.VarAttribute VarAttribute];
     684                for ivar=1:numel(VarAttribute)
     685                    ProjData.VarDimName=[ProjData.VarDimName {XName}];
     686                    ProjData.VarAttribute{ivar+nbvar}.Role='continuous';% will promote plots of the profiles with continuous lines
     687                    ProjData.(ListFieldProj{ivar})=ProjVar{ivar};
     688                end
     689            end
     690        case 'tps'%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     691            if strcmp(ProjMode,'interp_tps')
     692                Coord=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex});
     693                NbCentres=FieldData.(FieldData.ListVarName{CellInfo{icell}.NbCentres_tps});
     694                SubRange=FieldData.(FieldData.ListVarName{CellInfo{icell}.SubRange_tps});
     695                if isfield(CellInfo{icell},'VarIndex_vector_x_tps')&&isfield(CellInfo{icell},'VarIndex_vector_y_tps')
     696                    FieldVar=cat(3,FieldData.(FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_x_tps}),FieldData.(FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_y_tps}));
     697                end
     698                %                 coord_x_proj=XMin:DX:XMax;
     699                %                 coord_y_proj=YMin:DY:YMax;
     700                %                 np_x=numel(coord_x_proj);
     701                %                 np_y=numel(coord_y_proj);
     702                [DataOut,VarAttribute,errormsg]=calc_field_tps(Coord,NbCentres,SubRange,FieldVar,CellInfo{icell}.FieldName,cat(3,XI,YI));
     703                ProjData.X=Xproj;
     704                ProjData.ListVarName=[ProjData.ListVarName {XName}];
     705                ProjData.VarDimName=[ProjData.VarDimName {XName}];
     706                nbvar=numel(ProjData.ListVarName);
     707                ProjData.VarAttribute{nbvar}.long_name='abscissa along line';
     708                ProjVarName=(fieldnames(DataOut))';
     709                ProjData.ListVarName=[ProjData.ListVarName ProjVarName];
     710                ProjData.VarAttribute=[ProjData.VarAttribute VarAttribute];
     711                for ivar=1:numel(VarAttribute)
     712                    ProjData.VarDimName=[ProjData.VarDimName {XName}];
     713                    ProjData.VarAttribute{ivar+nbvar}.Role='continuous';% will promote plots of the profiles with continuous lines
     714                    ProjData.(ProjVarName{ivar})=DataOut.(ProjVarName{ivar});
     715                end
     716            end
     717            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     718           
     719        case 'grid'   %case of structured coordinates
     720            if ~isequal(ObjectData.Type,'line')% exclude polyline
     721                errormsg=['no  projection available on ' ObjectData.Type 'for structured coordinates']; %
    600722            else
    601                 errormsg='DX missing';
    602                 return
    603             end
    604         end
    605         XName= FieldData.ListVarName{CellInfo{icell}.CoordIndex(end)};
    606         YName= FieldData.ListVarName{CellInfo{icell}.CoordIndex(end-1)};
    607         coord_x=FieldData.(XName);   
    608         coord_y=FieldData.(YName);
    609     end   
    610    
    611     %% initiate projection
    612     for ivar=1:length(VarIndex)
    613         ProjLine{ivar}=[];
    614     end
    615     XLine=[];
    616     linelengthtot=0;
    617 
    618 %         circul=0;
    619 %         flux=0;
    620   %%%%%%%  % A FAIRE CALCULER MEAN DES QUANTITES    %%%%%%
    621    %case of unstructured coordinates
    622    if strcmp(CellInfo{icell}.CoordType,'scattered')
    623        for ip=1:siz_line(1)-1     %Loop on the segments of the polyline
    624            linelength=sqrt(dlinx(ip)*dlinx(ip)+dliny(ip)*dliny(ip));
    625            %select the vector indices in the range of action
    626            if testfalse
    627                flagsel=(errorflag==0); % keep only non false vectors
    628            else
    629                flagsel=ones(size(coord_x));
    630            end
    631            if isequal(ProjMode,'projection') %|| isequal(ProjMode,'interp_tps')
    632                flagsel=flagsel & ((coord_y -yinf(ip))*(xinf(ip+1)-xinf(ip))>(coord_x-xinf(ip))*(yinf(ip+1)-yinf(ip))) ...
    633                    & ((coord_y -ysup(ip))*(xsup(ip+1)-xsup(ip))<(coord_x-xsup(ip))*(ysup(ip+1)-ysup(ip))) ...
    634                    & ((coord_y -yinf(ip+1))*(xsup(ip+1)-xinf(ip+1))>(coord_x-xinf(ip+1))*(ysup(ip+1)-yinf(ip+1))) ...
    635                    & ((coord_y -yinf(ip))*(xsup(ip)-xinf(ip))<(coord_x-xinf(ip))*(ysup(ip)-yinf(ip)));
    636            end
    637            indsel=find(flagsel);%indsel =indices of good vectors
    638            X_sel=coord_x(indsel);
    639            Y_sel=coord_y(indsel);
    640            nbvar=0;
    641            for iselect=1:numel(VarIndex)-2*testU
    642                VarName=FieldData.ListVarName{VarIndex(iselect)};
    643                ProjVar{iselect}=FieldData.(VarName)(indsel);%scalar value
    644            end
    645            if testU
    646                ProjVar{numel(VarIndex)-1}=cos(theta(ip))*vector_x(indsel)+sin(theta(ip))*vector_y(indsel);% longitudinal component
    647                ProjVar{numel(VarIndex)}=-sin(theta(ip))*vector_x(indsel)+cos(theta(ip))*vector_y(indsel);%transverse component
    648            end
    649            if isequal(ProjMode,'projection')
    650                sintheta=sin(theta(ip));
    651                costheta=cos(theta(ip));
    652                Xproj=(X_sel-ObjectData.Coord(ip,1))*costheta + (Y_sel-ObjectData.Coord(ip,2))*sintheta; %projection on the line
    653                [Xproj,indsort]=sort(Xproj);
    654                for ivar=1:numel(ProjVar)
    655                    if ~isempty(ProjVar{ivar})
    656                        ProjVar{ivar}=ProjVar{ivar}(indsort);
    657                    end
    658                end
    659            elseif isequal(ProjMode,'interp_lin')||isequal(ProjMode,'interp_tps') %filtering %linear interpolation:
    660                npoint=floor(linelength/DX)+1;% nbre of points in the profile (interval DX)
    661                Xproj=linelength/(2*npoint):linelength/npoint:linelength-linelength/(2*npoint);
    662                xreg=cos(theta(ip))*Xproj+ObjectData.Coord(ip,1);
    663                yreg=sin(theta(ip))*Xproj+ObjectData.Coord(ip,2);
    664                if isfield(CellInfo{icell},'VarIndex_vector_x')&&isfield(CellInfo{icell},'VarIndex_vector_y')
    665                    VarName_x=FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_x};
    666                    VarName_y=FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_y};
    667                    if isfield(CellInfo{icell},'VarIndex_errorflag')
    668                        FieldData.(VarName_x)=FieldData.(VarName_x)(indsel);
    669                        FieldData.(VarName_y)=FieldData.(VarName_y)(indsel);
    670                    end
    671                    if ~isfield(CellInfo{icell},'CheckSub') || ~CellInfo{icell}.CheckSub
    672                        vector_x_proj=numel(ProjData.ListVarName)+1;
    673                        vector_y_proj=numel(ProjData.ListVarName)+2;
    674                    end
    675                end
    676                if isfield(CellInfo{icell},'VarIndex_scalar')
    677                    VarName_scalar=FieldData.ListVarName{CellInfo{icell}.VarIndex_scalar};
    678                    if isfield(CellInfo{icell},'errorflag') && ~isempty(CellInfo{icell}.errorflag)
    679                        FieldData.(VarName_scalar)=FieldData.(VarName_scalar)(indsel);
    680                    end
    681                end
    682                if isfield(CellInfo{icell},'VarIndex_ancillary')% do not project ancillary data with interp
    683                    FieldData=rmfield(FieldData,FieldData.ListVarName{CellInfo{icell}.VarIndex_ancillary});
    684                end
    685                if isfield(CellInfo{icell},'VarIndex_warnflag')% do not project ancillary data with interp
    686                    FieldData=rmfield(FieldData,FieldData.ListVarName{CellInfo{icell}.VarIndex_warnflag});
    687                end
    688                if isfield(CellInfo{icell},'VarIndex_errorflag')% do not project ancillary data with interp
    689                    FieldData=rmfield(FieldData,FieldData.ListVarName{CellInfo{icell}.VarIndex_errorflag});
    690                end
    691                if isequal(ProjMode,'interp_lin')
    692                [ProjVar,ListFieldProj,VarAttribute,errormsg]=calc_field_interp([X_sel Y_sel],FieldData,CellInfo{icell}.FieldName,xreg',yreg');
    693                else
    694                   [ProjVar,ListFieldProj,VarAttribute,errormsg]=calc_field_tps([X_sel Y_sel],FieldData,CellInfo{icell}.FieldName,xreg',yreg');
    695                end
    696                ivar_vector_x=[];
    697                ivar_vector_y=[];
    698                for ivar=1:numel(VarAttribute)
    699                    if isfield(VarAttribute{ivar},'Role')
    700                        if strcmp(VarAttribute{ivar}.Role,'vector_x')
    701                        ivar_vector_x=ivar;
    702                        elseif strcmp(VarAttribute{ivar}.Role,'vector_y')
    703                    ivar_vector_y=ivar;
    704                        end
    705                    end
    706                end
    707                if ~isempty(ivar_vector_x)&&~isempty(ivar_vector_y)
    708                                   ProjVar{ivar_vector_x}=cos(theta(ip))*ProjVar{ivar_vector_x}+sin(theta(ip))*ProjVar{ivar_vector_y};% longitudinal component
    709                ProjVar{ivar_vector_y}=-sin(theta(ip))*ProjVar{ivar_vector_x}+cos(theta(ip))*ProjVar{ivar_vector_y};%transverse component
    710                end
    711            elseif isequal(ProjMode,'interp_tps') %filtering
    712                %   TODO
    713                %                 npoint=floor(linelength/DX)+1;% nbre of points in the profile (interval DX)
    714                %                 Xproj=linelength/(2*npoint):linelength/npoint:linelength-linelength/(2*npoint);
    715                %                 siz=size(X_sel);
    716                %                 xregij=cos(theta(ip))*ones(siz(1),1)*Xproj+ObjectData.Coord(ip,1);
    717                %                 yregij=sin(theta(ip))*ones(siz(1),1)*Xproj+ObjectData.Coord(ip,2);
    718                %                 xij=X_sel*ones(1,npoint);
    719                %                 yij=Y_sel*ones(1,npoint);
    720                %                 Aij=exp(-lambda*((xij-xregij).*(xij-xregij)+(yij-yregij).*(yij-yregij)));
    721                %                 norm=Aij'*ones(siz(1),1);
    722                %                 for ivar=1:numel(ProjVar)
    723                %                      if ~isempty(ProjVar{ivar})
    724                %                         ProjVar{ivar}=Aij'*ProjVar{ivar}./norm;
    725                %
    726                %                      end
    727                %                 end
    728            end
    729            %prolongate the total record
    730            for ivar=1:numel(ProjVar)
    731                if ~isempty(ProjVar{ivar})
    732                    if numel(ProjLine)>=ivar
    733                        ProjLine{ivar}=[ProjLine{ivar}; ProjVar{ivar}];
    734                    else
    735                        ProjLine{ivar}=ProjVar{ivar};
    736                    end
    737                end
    738            end
    739            XLine=[XLine ;(Xproj+linelengthtot)];%along line abscissa
    740            linelengthtot=linelengthtot+linelength;
    741            %     circul=circul+(sum(U_sel))*linelength/npoint;
    742            %     flux=flux+(sum(V_sel))*linelength/npoint;
    743        end
    744        ProjData.X=XLine';
    745        ProjData.ListVarName=[ProjData.ListVarName {XName}];
    746        ProjData.VarDimName=[ProjData.VarDimName {XName}];
    747        ProjData.VarAttribute{1}.long_name='abscissa along line';
    748        for iselect=1:numel(VarIndex)
    749            VarName=FieldData.ListVarName{VarIndex(iselect)};
    750            eval(['ProjData.' VarName '=ProjLine{iselect};'])
    751            ProjData.ListVarName=[ProjData.ListVarName {VarName}];
    752            ProjData.VarDimName=[ProjData.VarDimName {XName}];
    753            ProjData.VarAttribute{iselect}=FieldData.VarAttribute{VarIndex(iselect)};
    754            if strcmp(ProjMode,'projection')
    755                ProjData.VarAttribute{iselect}.Role='discrete';
    756            else
    757                ProjData.VarAttribute{iselect}.Role='continuous';
    758            end
    759        end
    760        
    761        %case of structured coordinates
    762    elseif strcmp(CellInfo{icell}.CoordType,'grid')
    763        if ~isequal(ObjectData.Type,'line')% exclude polyline
    764            errormsg=['no  projection available on ' ObjectData.Type 'for structured coordinates']; %
    765        else
    766            test_Amat=1;%image or 2D matrix
    767            test_interp2=0;%default
    768            AYName=FieldData.ListVarName{CellInfo{icell}.CoordIndex(end-1)};
    769            AXName=FieldData.ListVarName{CellInfo{icell}.CoordIndex(end)};
    770            eval(['AX=FieldData.' AXName ';']);% set of x positions
    771            eval(['AY=FieldData.' AYName ';']);% set of y positions
    772            AName=FieldData.ListVarName{VarIndex(1)};
    773            eval(['A=FieldData.' AName ';']);% scalar
    774            npxy=size(A);
    775            npx=npxy(2);
    776            npy=npxy(1);
    777            if numel(AX)==2
    778                DX=(AX(2)-AX(1))/(npx-1);
    779            else
    780                DX_vec=diff(AX);
    781                DX=max(DX_vec);
    782                DX_min=min(DX_vec);
    783                if (DX-DX_min)>0.0001*abs(DX)
    784                    test_interp2=1;
    785                    DX=DX_min;
    786                end
    787            end
    788            if numel(AY)==2
    789                DY=(AY(2)-AY(1))/(npy-1);
    790            else
    791                DY_vec=diff(AY);
    792                DY=max(DY_vec);
    793                DY_min=min(DY_vec);
    794                if (DY-DY_min)>0.0001*abs(DY)
    795                    test_interp2=1;
    796                    DY=DY_min;
    797                end
    798            end
    799            AXI=linspace(AX(1),AX(end), npx);%set of  x  positions for the interpolated input data
    800            AYI=linspace(AY(1),AY(end), npy);%set of  x  positions for the interpolated input data
    801            if isfield(ObjectData,'DX')
    802                DXY_line=ObjectData.DX;%mesh on the projection line
    803            else
    804                DXY_line=sqrt(abs(DX*DY));% mesh on the projection line
    805            end
    806            dlinx=ObjectData.Coord(2,1)-ObjectData.Coord(1,1);
    807            dliny=ObjectData.Coord(2,2)-ObjectData.Coord(1,2);
    808            linelength=sqrt(dlinx*dlinx+dliny*dliny);
    809            theta=angle(dlinx+i*dliny);%angle of the line
    810            if isfield(FieldData,'RangeX')
    811                XMin=min(FieldData.RangeX);%shift of the origin on the line
    812            else
    813                XMin=0;
    814            end
    815            eval(['ProjData.' AXName '=linspace(XMin,XMin+linelength,linelength/DXY_line+1);'])%abscissa of the new pixels along the line
    816            y=linspace(-width,width,2*width/DXY_line+1);%ordintes of the new pixels (coordinate across the line)
    817            eval(['npX=length(ProjData.' AXName ');'])
    818            npY=length(y); %TODO: utiliser proj_grid
    819            eval(['[X,Y]=meshgrid(ProjData.' AXName ',y);'])%grid in the line coordinates
    820            XIMA=ObjectData.Coord(1,1)+(X-XMin)*cos(theta)-Y*sin(theta);
    821            YIMA=ObjectData.Coord(1,2)+(X-XMin)*sin(theta)+Y*cos(theta);
    822            XIMA=(XIMA-AX(1))/DX+1;%  index of the original image along x
    823            YIMA=(YIMA-AY(1))/DY+1;% index of the original image along y
    824            XIMA=reshape(round(XIMA),1,npX*npY);%indices reorganized in 'line'
    825            YIMA=reshape(round(YIMA),1,npX*npY);
    826            flagin=XIMA>=1 & XIMA<=npx & YIMA >=1 & YIMA<=npy;%flagin=1 inside the original image
    827            ind_in=find(flagin);
    828            ind_out=find(~flagin);
    829            ICOMB=(XIMA-1)*npy+YIMA;
    830            ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A
    831            nbcolor=1; %color images
    832            if numel(npxy)==2
    833                nbcolor=1;
    834            elseif length(npxy)==3
    835                nbcolor=npxy(3);
    836            else
    837                errormsg='multicomponent field not projected';
    838                display(errormsg)
    839                return
    840            end
    841            nbvar=length(ProjData.ListVarName);% number of var from previous cells
    842            ProjData.ListVarName=[ProjData.ListVarName {AXName}];
    843            ProjData.VarDimName=[ProjData.VarDimName {AXName}];
    844            for ivar=VarIndex
    845                %VarName{ivar}=FieldData.ListVarName{ivar};
    846                if test_interp2% interpolate on new grid
    847                    FieldData.(FieldData.ListVarName{ivar})=interp2(FieldData.(AXName),FieldData.(AYName),FieldData.(FieldData.ListVarName{ivar}),AXI,AYI);%TO TEST
    848                end
    849                vec_A=reshape(squeeze(FieldData.(FieldData.ListVarName{ivar})),npx*npy,nbcolor); %put the original image in colum
    850                if nbcolor==1
    851                    vec_B(ind_in)=vec_A(ICOMB);
    852                    vec_B(ind_out)=zeros(size(ind_out));
    853                    A_out=reshape(vec_B,npY,npX);
    854                    ProjData.(FieldData.ListVarName{ivar}) =sum(A_out,1)/npY;
    855                elseif nbcolor==3
    856                    vec_B(ind_in,1:3)=vec_A(ICOMB,:);
    857                    vec_B(ind_out,1)=zeros(size(ind_out));
    858                    vec_B(ind_out,2)=zeros(size(ind_out));
    859                    vec_B(ind_out,3)=zeros(size(ind_out));
    860                    A_out=reshape(vec_B,npY,npX,nbcolor);
    861                    ProjData.(FieldData.ListVarName{ivar})=squeeze(sum(A_out,1)/npY);
    862                end
    863                ProjData.ListVarName=[ProjData.ListVarName FieldData.ListVarName{ivar}];
    864                ProjData.VarDimName=[ProjData.VarDimName {AXName}];%to generalize with the initial name of the x coordinate
    865                ProjData.VarAttribute{ivar}.Role='continuous';% for plot with continuous line
    866            end
    867            if testU
    868                vector_x =ProjData.(FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_x});
    869                vector_y =ProjData.(FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_y});
    870                ProjData.(FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_x}) =cos(theta)*vector_x+sin(theta)*vector_y;
    871                ProjData.(FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_y}) =-sin(theta)*vector_x+cos(theta)*vector_y;
    872            end
    873            ProjData.VarAttribute{nbvar+1}.long_name='abscissa along line';
    874            if nbcolor==3
    875                ProjData.VarDimName{end}={AXName,'rgb'};
    876            end
    877        end
    878    elseif strcmp(CellInfo{icell}.CoordType,'tps')
    879        if isfield(ObjectData,'DX')&~isempty(ObjectData.DX)
    880            DX=abs(ObjectData.DX);%mesh of interpolation points along the line
    881            Xproj=linelength/(2*npoint):linelength/npoint:linelength-linelength/(2*npoint);
    882            xreg=cos(theta(ip))*Xproj+ObjectData.Coord(ip,1)
    883            yreg=sin(theta(ip))*Xproj+ObjectData.Coord(ip,2)
    884            %                 coord_x_proj=XMin:DX:XMax;
    885            %                 coord_y_proj=YMin:DY:YMax;
    886            DataOut=calc_field_tps(FieldData.FieldList,FieldData,cat(3,xreg,yreg));
    887            ProjData.ListVarName=[ProjData.ListVarName DataOut.ListVarName];
    888            ProjData.VarDimName=[ProjData.VarDimName DataOut.VarDimName];
    889            ProjData.VarAttribute=[ProjData.VarAttribute DataOut.VarAttribute];
    890            DataOut.ListVarName(1)=[];
    891            DataOut.VarDimName(1)=[];
    892            DataOut.VarAttribute(1)=[];
    893            for ilist=2:length(DataOut.ListVarName)% reshape data, excluding coordinates (ilist=1-2), TODO: rationalise
    894                VarName=DataOut.ListVarName{ilist};
    895                ProjData.(VarName)=DataOut.(VarName);
    896            end
    897            ProjData.coord_x=Xproj;
    898        end
    899    end
     723                test_Amat=1;%image or 2D matrix
     724                test_interp2=0;%default
     725                AYName=FieldData.ListVarName{CellInfo{icell}.CoordIndex(end-1)};
     726                AXName=FieldData.ListVarName{CellInfo{icell}.CoordIndex(end)};
     727                eval(['AX=FieldData.' AXName ';']);% set of x positions
     728                eval(['AY=FieldData.' AYName ';']);% set of y positions
     729                AName=FieldData.ListVarName{VarIndex(1)};
     730                eval(['A=FieldData.' AName ';']);% scalar
     731                npxy=size(A);
     732                npx=npxy(2);
     733                npy=npxy(1);
     734                if numel(AX)==2
     735                    DX=(AX(2)-AX(1))/(npx-1);
     736                else
     737                    DX_vec=diff(AX);
     738                    DX=max(DX_vec);
     739                    DX_min=min(DX_vec);
     740                    if (DX-DX_min)>0.0001*abs(DX)
     741                        test_interp2=1;
     742                        DX=DX_min;
     743                    end
     744                end
     745                if numel(AY)==2
     746                    DY=(AY(2)-AY(1))/(npy-1);
     747                else
     748                    DY_vec=diff(AY);
     749                    DY=max(DY_vec);
     750                    DY_min=min(DY_vec);
     751                    if (DY-DY_min)>0.0001*abs(DY)
     752                        test_interp2=1;
     753                        DY=DY_min;
     754                    end
     755                end
     756                AXI=linspace(AX(1),AX(end), npx);%set of  x  positions for the interpolated input data
     757                AYI=linspace(AY(1),AY(end), npy);%set of  x  positions for the interpolated input data
     758                if isfield(ObjectData,'DX')
     759                    DXY_line=ObjectData.DX;%mesh on the projection line
     760                else
     761                    DXY_line=sqrt(abs(DX*DY));% mesh on the projection line
     762                end
     763                dlinx=ObjectData.Coord(2,1)-ObjectData.Coord(1,1);
     764                dliny=ObjectData.Coord(2,2)-ObjectData.Coord(1,2);
     765                linelength=sqrt(dlinx*dlinx+dliny*dliny);
     766                theta=angle(dlinx+i*dliny);%angle of the line
     767                if isfield(FieldData,'RangeX')
     768                    XMin=min(FieldData.RangeX);%shift of the origin on the line
     769                else
     770                    XMin=0;
     771                end
     772                eval(['ProjData.' AXName '=linspace(XMin,XMin+linelength,linelength/DXY_line+1);'])%abscissa of the new pixels along the line
     773                y=linspace(-width,width,2*width/DXY_line+1);%ordintes of the new pixels (coordinate across the line)
     774                eval(['npX=length(ProjData.' AXName ');'])
     775                npY=length(y); %TODO: utiliser proj_grid
     776                eval(['[X,Y]=meshgrid(ProjData.' AXName ',y);'])%grid in the line coordinates
     777                XIMA=ObjectData.Coord(1,1)+(X-XMin)*cos(theta)-Y*sin(theta);
     778                YIMA=ObjectData.Coord(1,2)+(X-XMin)*sin(theta)+Y*cos(theta);
     779                XIMA=(XIMA-AX(1))/DX+1;%  index of the original image along x
     780                YIMA=(YIMA-AY(1))/DY+1;% index of the original image along y
     781                XIMA=reshape(round(XIMA),1,npX*npY);%indices reorganized in 'line'
     782                YIMA=reshape(round(YIMA),1,npX*npY);
     783                flagin=XIMA>=1 & XIMA<=npx & YIMA >=1 & YIMA<=npy;%flagin=1 inside the original image
     784                ind_in=find(flagin);
     785                ind_out=find(~flagin);
     786                ICOMB=(XIMA-1)*npy+YIMA;
     787                ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A
     788                nbcolor=1; %color images
     789                if numel(npxy)==2
     790                    nbcolor=1;
     791                elseif length(npxy)==3
     792                    nbcolor=npxy(3);
     793                else
     794                    errormsg='multicomponent field not projected';
     795                    display(errormsg)
     796                    return
     797                end
     798                nbvar=length(ProjData.ListVarName);% number of var from previous cells
     799                ProjData.ListVarName=[ProjData.ListVarName {AXName}];
     800                ProjData.VarDimName=[ProjData.VarDimName {AXName}];
     801                for ivar=VarIndex
     802                    %VarName{ivar}=FieldData.ListVarName{ivar};
     803                    if test_interp2% interpolate on new grid
     804                        FieldData.(FieldData.ListVarName{ivar})=interp2(FieldData.(AXName),FieldData.(AYName),FieldData.(FieldData.ListVarName{ivar}),AXI,AYI);%TO TEST
     805                    end
     806                    vec_A=reshape(squeeze(FieldData.(FieldData.ListVarName{ivar})),npx*npy,nbcolor); %put the original image in colum
     807                    if nbcolor==1
     808                        vec_B(ind_in)=vec_A(ICOMB);
     809                        vec_B(ind_out)=zeros(size(ind_out));
     810                        A_out=reshape(vec_B,npY,npX);
     811                        ProjData.(FieldData.ListVarName{ivar}) =sum(A_out,1)/npY;
     812                    elseif nbcolor==3
     813                        vec_B(ind_in,1:3)=vec_A(ICOMB,:);
     814                        vec_B(ind_out,1)=zeros(size(ind_out));
     815                        vec_B(ind_out,2)=zeros(size(ind_out));
     816                        vec_B(ind_out,3)=zeros(size(ind_out));
     817                        A_out=reshape(vec_B,npY,npX,nbcolor);
     818                        ProjData.(FieldData.ListVarName{ivar})=squeeze(sum(A_out,1)/npY);
     819                    end
     820                    ProjData.ListVarName=[ProjData.ListVarName FieldData.ListVarName{ivar}];
     821                    ProjData.VarDimName=[ProjData.VarDimName {AXName}];%to generalize with the initial name of the x coordinate
     822                    ProjData.VarAttribute{ivar}.Role='continuous';% for plot with continuous line
     823                end
     824                if testU
     825                    vector_x =ProjData.(FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_x});
     826                    vector_y =ProjData.(FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_y});
     827                    ProjData.(FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_x}) =cos(theta)*vector_x+sin(theta)*vector_y;
     828                    ProjData.(FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_y}) =-sin(theta)*vector_x+cos(theta)*vector_y;
     829                end
     830                ProjData.VarAttribute{nbvar+1}.long_name='abscissa along line';
     831                if nbcolor==3
     832                    ProjData.VarDimName{end}={AXName,'rgb'};
     833                end
     834            end
     835    end
    900836end
    901837
Note: See TracChangeset for help on using the changeset viewer.