Changeset 644 for trunk/src/proj_field.m


Ignore:
Timestamp:
May 28, 2013, 11:30:01 PM (8 years ago)
Author:
sommeria
Message:

various improvements: resize GUI uvmat, projection on lines

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/proj_field.m

    r630 r644  
    564564   VarIndex=find(check_proj);
    565565
    566     %identify vector components   
     566    %% identify vector components   
    567567    testU=isfield(CellInfo{icell},'VarIndex_vector_x') &&isfield(CellInfo{icell},'VarIndex_vector_y') ;% test for vectors
    568568    if testU
     
    578578        errorflag=FieldData.(FFName);
    579579    end   
    580     % check needed object properties for unstructured positions (position given by the variables with role coord_x, coord_y
     580   
     581    %% check needed object properties for unstructured positions (position given by the variables with role coord_x, coord_y
    581582    if strcmp(CellInfo{icell}.CoordType,'scattered')
    582         if  ~isequal(ProjMode,'interp_lin')
     583        if  strcmp(ProjMode,'projection')
    583584            if width==0
    584585                errormsg='range of the projection object is missing';
    585                 return     
    586             else
    587                 lambda=2/(width*width); %smoothing factor used for interp_tps: weight exp(-2) at distance width from the line
    588             end
    589         end
    590         if ~isequal(ProjMode,'projection')
    591             if isfield(ObjectData,'DX')&~isempty(ObjectData.DX)
     586                return 
     587            end
     588%             else
     589%                 lambda=2/(width*width); %smoothing factor used for interp_tps: weight exp(-2) at distance width from the line
     590%             end
     591        else
     592            if isfield(ObjectData,'DX') && ~isempty(ObjectData.DX)
    592593                DX=abs(ObjectData.DX);%mesh of interpolation points along the line
    593594            else
     
    601602        coord_y=FieldData.(YName);
    602603    end   
    603     %initiate projection
     604   
     605    %% initiate projection
    604606    for ivar=1:length(VarIndex)
    605607        ProjLine{ivar}=[];
     
    612614  %%%%%%%  % A FAIRE CALCULER MEAN DES QUANTITES    %%%%%%
    613615   %case of unstructured coordinates
    614     if strcmp(CellInfo{icell}.CoordType,'scattered') 
    615         for ip=1:siz_line(1)-1     %Loop on the segments of the polyline
    616             linelength=sqrt(dlinx(ip)*dlinx(ip)+dliny(ip)*dliny(ip)); 
    617             %select the vector indices in the range of action
    618             if testfalse
    619                 flagsel=(errorflag==0); % keep only non false vectors
    620             else
    621                 flagsel=ones(size(coord_x));
    622             end
    623             if isequal(ProjMode,'projection') %|| isequal(ProjMode,'interp_tps')
    624                 flagsel=flagsel & ((coord_y -yinf(ip))*(xinf(ip+1)-xinf(ip))>(coord_x-xinf(ip))*(yinf(ip+1)-yinf(ip))) ...
    625                 & ((coord_y -ysup(ip))*(xsup(ip+1)-xsup(ip))<(coord_x-xsup(ip))*(ysup(ip+1)-ysup(ip))) ...
    626                 & ((coord_y -yinf(ip+1))*(xsup(ip+1)-xinf(ip+1))>(coord_x-xinf(ip+1))*(ysup(ip+1)-yinf(ip+1))) ...
    627                 & ((coord_y -yinf(ip))*(xsup(ip)-xinf(ip))<(coord_x-xinf(ip))*(ysup(ip)-yinf(ip)));
    628             end
    629             indsel=find(flagsel);%indsel =indices of good vectors
    630             X_sel=coord_x(indsel);
    631             Y_sel=coord_y(indsel);
    632             nbvar=0;
    633             for iselect=1:numel(VarIndex)-2*testU
    634                 VarName=FieldData.ListVarName{VarIndex(iselect)};
    635                 ProjVar{iselect}=FieldData.(VarName)(indsel);%scalar value
    636             end   
    637             if testU
    638                 ProjVar{numel(VarIndex)-1}=cos(theta(ip))*vector_x(indsel)+sin(theta(ip))*vector_y(indsel);% longitudinal component
    639                 ProjVar{numel(VarIndex)}=-sin(theta(ip))*vector_x(indsel)+cos(theta(ip))*vector_y(indsel);%transverse component         
    640             end
    641             if isequal(ProjMode,'projection')
    642                 sintheta=sin(theta(ip));
    643                 costheta=cos(theta(ip));
    644                 Xproj=(X_sel-ObjectData.Coord(ip,1))*costheta + (Y_sel-ObjectData.Coord(ip,2))*sintheta; %projection on the line
    645                 [Xproj,indsort]=sort(Xproj);
    646                 for ivar=1:numel(ProjVar)
    647                     if ~isempty(ProjVar{ivar})
    648                         ProjVar{ivar}=ProjVar{ivar}(indsort);
    649                      end
    650                 end
    651             elseif isequal(ProjMode,'interp_lin') %linear interpolation:
    652                 npoint=floor(linelength/DX)+1;% nbre of points in the profile (interval DX)
    653                 Xproj=linelength/(2*npoint):linelength/npoint:linelength-linelength/(2*npoint);
    654                 xreg=cos(theta(ip))*Xproj+ObjectData.Coord(ip,1);
    655                 yreg=sin(theta(ip))*Xproj+ObjectData.Coord(ip,2);
    656                 for ivar=1:numel(ProjVar)
    657                      if ~isempty(ProjVar{ivar})
    658                         ProjVar{ivar}=griddata_uvmat(X_sel,Y_sel,ProjVar{ivar},xreg,yreg);
    659                      end
    660                 end
    661             elseif isequal(ProjMode,'interp_tps') %filtering
    662                 npoint=floor(linelength/DX)+1;% nbre of points in the profile (interval DX)
    663                 Xproj=linelength/(2*npoint):linelength/npoint:linelength-linelength/(2*npoint);
    664                 siz=size(X_sel);
    665                 xregij=cos(theta(ip))*ones(siz(1),1)*Xproj+ObjectData.Coord(ip,1);
    666                 yregij=sin(theta(ip))*ones(siz(1),1)*Xproj+ObjectData.Coord(ip,2);
    667                 xij=X_sel*ones(1,npoint);
    668                 yij=Y_sel*ones(1,npoint);
    669                 Aij=exp(-lambda*((xij-xregij).*(xij-xregij)+(yij-yregij).*(yij-yregij)));
    670                 norm=Aij'*ones(siz(1),1);
    671                 for ivar=1:numel(ProjVar)
    672                      if ~isempty(ProjVar{ivar})
    673                         ProjVar{ivar}=Aij'*ProjVar{ivar}./norm;
    674                      end
    675                 end             
    676             end
    677             %prolongate the total record
    678             for ivar=1:numel(ProjVar)
    679                   if ~isempty(ProjVar{ivar})
    680                      ProjLine{ivar}=[ProjLine{ivar}; ProjVar{ivar}];
    681                   end
    682             end
    683             XLine=[XLine ;(Xproj+linelengthtot)];%along line abscissa
    684             linelengthtot=linelengthtot+linelength;
    685             %     circul=circul+(sum(U_sel))*linelength/npoint;
    686             %     flux=flux+(sum(V_sel))*linelength/npoint;
    687         end
    688         ProjData.X=XLine';
    689         ProjData.ListVarName=[ProjData.ListVarName {XName}];
    690         ProjData.VarDimName=[ProjData.VarDimName {XName}];
    691         ProjData.VarAttribute{1}.long_name='abscissa along line';
    692         for iselect=1:numel(VarIndex)
    693             VarName=FieldData.ListVarName{VarIndex(iselect)};
    694             eval(['ProjData.' VarName '=ProjLine{iselect};'])
    695             ProjData.ListVarName=[ProjData.ListVarName {VarName}];
    696             ProjData.VarDimName=[ProjData.VarDimName {XName}];
    697             ProjData.VarAttribute{iselect}=FieldData.VarAttribute{VarIndex(iselect)};
    698             if strcmp(ProjMode,'projection')
    699                 ProjData.VarAttribute{iselect}.Role='discrete';
    700             else
    701                  ProjData.VarAttribute{iselect}.Role='continuous';
    702             end
    703         end
    704    
    705     %case of structured coordinates
    706     elseif strcmp(CellInfo{icell}.CoordType,'grid')
    707         if ~isequal(ObjectData.Type,'line')% exclude polyline
    708             errormsg=['no  projection available on ' ObjectData.Type 'for structured coordinates']; %
    709         else
    710             test_Amat=1;%image or 2D matrix
    711             test_interp2=0;%default
    712             AYName=FieldData.ListVarName{CellInfo{icell}.CoordIndex(end-1)};
    713             AXName=FieldData.ListVarName{CellInfo{icell}.CoordIndex(end)};
    714             eval(['AX=FieldData.' AXName ';']);% set of x positions
    715             eval(['AY=FieldData.' AYName ';']);% set of y positions 
    716             AName=FieldData.ListVarName{VarIndex(1)};
    717             eval(['A=FieldData.' AName ';']);% scalar
    718             npxy=size(A);
    719             npx=npxy(2);
    720             npy=npxy(1);
    721             if numel(AX)==2
    722                 DX=(AX(2)-AX(1))/(npx-1);
    723             else
    724                 DX_vec=diff(AX);
    725                 DX=max(DX_vec);
    726                 DX_min=min(DX_vec);
    727                 if (DX-DX_min)>0.0001*abs(DX)
    728                     test_interp2=1;
    729                     DX=DX_min;
    730                 end   
    731             end
    732             if numel(AY)==2
    733                 DY=(AY(2)-AY(1))/(npy-1);
    734             else
    735                 DY_vec=diff(AY);
    736                 DY=max(DY_vec);
    737                 DY_min=min(DY_vec);
    738                 if (DY-DY_min)>0.0001*abs(DY)
     616   if strcmp(CellInfo{icell}.CoordType,'scattered')
     617       for ip=1:siz_line(1)-1     %Loop on the segments of the polyline
     618           linelength=sqrt(dlinx(ip)*dlinx(ip)+dliny(ip)*dliny(ip));
     619           %select the vector indices in the range of action
     620           if testfalse
     621               flagsel=(errorflag==0); % keep only non false vectors
     622           else
     623               flagsel=ones(size(coord_x));
     624           end
     625           if isequal(ProjMode,'projection') %|| isequal(ProjMode,'interp_tps')
     626               flagsel=flagsel & ((coord_y -yinf(ip))*(xinf(ip+1)-xinf(ip))>(coord_x-xinf(ip))*(yinf(ip+1)-yinf(ip))) ...
     627                   & ((coord_y -ysup(ip))*(xsup(ip+1)-xsup(ip))<(coord_x-xsup(ip))*(ysup(ip+1)-ysup(ip))) ...
     628                   & ((coord_y -yinf(ip+1))*(xsup(ip+1)-xinf(ip+1))>(coord_x-xinf(ip+1))*(ysup(ip+1)-yinf(ip+1))) ...
     629                   & ((coord_y -yinf(ip))*(xsup(ip)-xinf(ip))<(coord_x-xinf(ip))*(ysup(ip)-yinf(ip)));
     630           end
     631           indsel=find(flagsel);%indsel =indices of good vectors
     632           X_sel=coord_x(indsel);
     633           Y_sel=coord_y(indsel);
     634           nbvar=0;
     635           for iselect=1:numel(VarIndex)-2*testU
     636               VarName=FieldData.ListVarName{VarIndex(iselect)};
     637               ProjVar{iselect}=FieldData.(VarName)(indsel);%scalar value
     638           end
     639           if testU
     640               ProjVar{numel(VarIndex)-1}=cos(theta(ip))*vector_x(indsel)+sin(theta(ip))*vector_y(indsel);% longitudinal component
     641               ProjVar{numel(VarIndex)}=-sin(theta(ip))*vector_x(indsel)+cos(theta(ip))*vector_y(indsel);%transverse component
     642           end
     643           if isequal(ProjMode,'projection')
     644               sintheta=sin(theta(ip));
     645               costheta=cos(theta(ip));
     646               Xproj=(X_sel-ObjectData.Coord(ip,1))*costheta + (Y_sel-ObjectData.Coord(ip,2))*sintheta; %projection on the line
     647               [Xproj,indsort]=sort(Xproj);
     648               for ivar=1:numel(ProjVar)
     649                   if ~isempty(ProjVar{ivar})
     650                       ProjVar{ivar}=ProjVar{ivar}(indsort);
     651                   end
     652               end
     653           elseif isequal(ProjMode,'interp_lin')||isequal(ProjMode,'interp_tps') %filtering %linear interpolation:
     654               npoint=floor(linelength/DX)+1;% nbre of points in the profile (interval DX)
     655               Xproj=linelength/(2*npoint):linelength/npoint:linelength-linelength/(2*npoint);
     656               xreg=cos(theta(ip))*Xproj+ObjectData.Coord(ip,1);
     657               yreg=sin(theta(ip))*Xproj+ObjectData.Coord(ip,2);
     658               if isfield(CellInfo{icell},'VarIndex_vector_x')&&isfield(CellInfo{icell},'VarIndex_vector_y')
     659                   VarName_x=FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_x};
     660                   VarName_y=FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_y};
     661                   if isfield(CellInfo{icell},'VarIndex_errorflag')
     662                       FieldData.(VarName_x)=FieldData.(VarName_x)(indsel);
     663                       FieldData.(VarName_y)=FieldData.(VarName_y)(indsel);
     664                   end
     665                   if ~isfield(CellInfo{icell},'CheckSub') || ~CellInfo{icell}.CheckSub
     666                       vector_x_proj=numel(ProjData.ListVarName)+1;
     667                       vector_y_proj=numel(ProjData.ListVarName)+2;
     668                   end
     669               end
     670               if isfield(CellInfo{icell},'VarIndex_scalar')
     671                   VarName_scalar=FieldData.ListVarName{CellInfo{icell}.VarIndex_scalar};
     672                   if isfield(CellInfo{icell},'errorflag') && ~isempty(CellInfo{icell}.errorflag)
     673                       FieldData.(VarName_scalar)=FieldData.(VarName_scalar)(indsel);
     674                   end
     675               end
     676               if isfield(CellInfo{icell},'VarIndex_ancillary')% do not project ancillary data with interp
     677                   FieldData=rmfield(FieldData,FieldData.ListVarName{CellInfo{icell}.VarIndex_ancillary});
     678               end
     679               if isfield(CellInfo{icell},'VarIndex_warnflag')% do not project ancillary data with interp
     680                   FieldData=rmfield(FieldData,FieldData.ListVarName{CellInfo{icell}.VarIndex_warnflag});
     681               end
     682               if isfield(CellInfo{icell},'VarIndex_errorflag')% do not project ancillary data with interp
     683                   FieldData=rmfield(FieldData,FieldData.ListVarName{CellInfo{icell}.VarIndex_errorflag});
     684               end
     685               if isequal(ProjMode,'interp_lin')
     686               [ProjVar,ListFieldProj,VarAttribute,errormsg]=calc_field_interp([X_sel Y_sel],FieldData,CellInfo{icell}.FieldName,xreg',yreg');
     687               else
     688                  [ProjVar,ListFieldProj,VarAttribute,errormsg]=calc_field_tps([X_sel Y_sel],FieldData,CellInfo{icell}.FieldName,xreg',yreg');
     689               end
     690               ivar_vector_x=[];
     691               ivar_vector_y=[];
     692               for ivar=1:numel(VarAttribute)
     693                   if isfield(VarAttribute{ivar},'Role')
     694                       if strcmp(VarAttribute{ivar}.Role,'vector_x')
     695                       ivar_vector_x=ivar;
     696                       elseif strcmp(VarAttribute{ivar}.Role,'vector_y')
     697                   ivar_vector_y=ivar;
     698                       end
     699                   end
     700               end
     701               if ~isempty(ivar_vector_x)&&~isempty(ivar_vector_y)
     702                                  ProjVar{ivar_vector_x}=cos(theta(ip))*ProjVar{ivar_vector_x}+sin(theta(ip))*ProjVar{ivar_vector_y};% longitudinal component
     703               ProjVar{ivar_vector_y}=-sin(theta(ip))*ProjVar{ivar_vector_x}+cos(theta(ip))*ProjVar{ivar_vector_y};%transverse component
     704               end
     705           elseif isequal(ProjMode,'interp_tps') %filtering
     706               %   TODO
     707               %                 npoint=floor(linelength/DX)+1;% nbre of points in the profile (interval DX)
     708               %                 Xproj=linelength/(2*npoint):linelength/npoint:linelength-linelength/(2*npoint);
     709               %                 siz=size(X_sel);
     710               %                 xregij=cos(theta(ip))*ones(siz(1),1)*Xproj+ObjectData.Coord(ip,1);
     711               %                 yregij=sin(theta(ip))*ones(siz(1),1)*Xproj+ObjectData.Coord(ip,2);
     712               %                 xij=X_sel*ones(1,npoint);
     713               %                 yij=Y_sel*ones(1,npoint);
     714               %                 Aij=exp(-lambda*((xij-xregij).*(xij-xregij)+(yij-yregij).*(yij-yregij)));
     715               %                 norm=Aij'*ones(siz(1),1);
     716               %                 for ivar=1:numel(ProjVar)
     717               %                      if ~isempty(ProjVar{ivar})
     718               %                         ProjVar{ivar}=Aij'*ProjVar{ivar}./norm;
     719               %
     720               %                      end
     721               %                 end
     722           end
     723           %prolongate the total record
     724           for ivar=1:numel(ProjVar)
     725               if ~isempty(ProjVar{ivar})
     726                   if numel(ProjLine)>=ivar
     727                       ProjLine{ivar}=[ProjLine{ivar}; ProjVar{ivar}];
     728                   else
     729                       ProjLine{ivar}=ProjVar{ivar};
     730                   end
     731               end
     732           end
     733           XLine=[XLine ;(Xproj+linelengthtot)];%along line abscissa
     734           linelengthtot=linelengthtot+linelength;
     735           %     circul=circul+(sum(U_sel))*linelength/npoint;
     736           %     flux=flux+(sum(V_sel))*linelength/npoint;
     737       end
     738       ProjData.X=XLine';
     739       ProjData.ListVarName=[ProjData.ListVarName {XName}];
     740       ProjData.VarDimName=[ProjData.VarDimName {XName}];
     741       ProjData.VarAttribute{1}.long_name='abscissa along line';
     742       for iselect=1:numel(VarIndex)
     743           VarName=FieldData.ListVarName{VarIndex(iselect)};
     744           eval(['ProjData.' VarName '=ProjLine{iselect};'])
     745           ProjData.ListVarName=[ProjData.ListVarName {VarName}];
     746           ProjData.VarDimName=[ProjData.VarDimName {XName}];
     747           ProjData.VarAttribute{iselect}=FieldData.VarAttribute{VarIndex(iselect)};
     748           if strcmp(ProjMode,'projection')
     749               ProjData.VarAttribute{iselect}.Role='discrete';
     750           else
     751               ProjData.VarAttribute{iselect}.Role='continuous';
     752           end
     753       end
     754       
     755       %case of structured coordinates
     756   elseif strcmp(CellInfo{icell}.CoordType,'grid')
     757       if ~isequal(ObjectData.Type,'line')% exclude polyline
     758           errormsg=['no  projection available on ' ObjectData.Type 'for structured coordinates']; %
     759       else
     760           test_Amat=1;%image or 2D matrix
     761           test_interp2=0;%default
     762           AYName=FieldData.ListVarName{CellInfo{icell}.CoordIndex(end-1)};
     763           AXName=FieldData.ListVarName{CellInfo{icell}.CoordIndex(end)};
     764           eval(['AX=FieldData.' AXName ';']);% set of x positions
     765           eval(['AY=FieldData.' AYName ';']);% set of y positions
     766           AName=FieldData.ListVarName{VarIndex(1)};
     767           eval(['A=FieldData.' AName ';']);% scalar
     768           npxy=size(A);
     769           npx=npxy(2);
     770           npy=npxy(1);
     771           if numel(AX)==2
     772               DX=(AX(2)-AX(1))/(npx-1);
     773           else
     774               DX_vec=diff(AX);
     775               DX=max(DX_vec);
     776               DX_min=min(DX_vec);
     777               if (DX-DX_min)>0.0001*abs(DX)
    739778                   test_interp2=1;
    740                     DY=DY_min;
    741                 end     
    742             end             
    743             AXI=linspace(AX(1),AX(end), npx);%set of  x  positions for the interpolated input data
    744             AYI=linspace(AY(1),AY(end), npy);%set of  x  positions for the interpolated input data
    745             if isfield(ObjectData,'DX')
    746                 DXY_line=ObjectData.DX;%mesh on the projection line
    747             else
    748                 DXY_line=sqrt(abs(DX*DY));% mesh on the projection line
    749             end
    750             dlinx=ObjectData.Coord(2,1)-ObjectData.Coord(1,1);
    751             dliny=ObjectData.Coord(2,2)-ObjectData.Coord(1,2);
    752             linelength=sqrt(dlinx*dlinx+dliny*dliny);
    753             theta=angle(dlinx+i*dliny);%angle of the line   
    754             if isfield(FieldData,'RangeX')
    755                 XMin=min(FieldData.RangeX);%shift of the origin on the line
    756             else
    757                 XMin=0;
    758             end
    759             eval(['ProjData.' AXName '=linspace(XMin,XMin+linelength,linelength/DXY_line+1);'])%abscissa of the new pixels along the line
    760             y=linspace(-width,width,2*width/DXY_line+1);%ordintes of the new pixels (coordinate across the line)
    761             eval(['npX=length(ProjData.' AXName ');'])
    762             npY=length(y); %TODO: utiliser proj_grid
    763             eval(['[X,Y]=meshgrid(ProjData.' AXName ',y);'])%grid in the line coordinates
    764             XIMA=ObjectData.Coord(1,1)+(X-XMin)*cos(theta)-Y*sin(theta);
    765             YIMA=ObjectData.Coord(1,2)+(X-XMin)*sin(theta)+Y*cos(theta);
    766             XIMA=(XIMA-AX(1))/DX+1;%  index of the original image along x
    767             YIMA=(YIMA-AY(1))/DY+1;% index of the original image along y
    768             XIMA=reshape(round(XIMA),1,npX*npY);%indices reorganized in 'line'
    769             YIMA=reshape(round(YIMA),1,npX*npY);
    770             flagin=XIMA>=1 & XIMA<=npx & YIMA >=1 & YIMA<=npy;%flagin=1 inside the original image
    771             ind_in=find(flagin);
    772             ind_out=find(~flagin);
    773             ICOMB=(XIMA-1)*npy+YIMA;
    774             ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A
    775             nbcolor=1; %color images
    776             if numel(npxy)==2
    777                 nbcolor=1;
    778             elseif length(npxy)==3
    779                 nbcolor=npxy(3);
    780             else
    781                 errormsg='multicomponent field not projected';
    782                 display(errormsg)
    783                 return
    784             end
    785             nbvar=length(ProjData.ListVarName);% number of var from previous cells
    786             ProjData.ListVarName=[ProjData.ListVarName {AXName}];
    787             ProjData.VarDimName=[ProjData.VarDimName {AXName}];
    788             for ivar=VarIndex
    789                 %VarName{ivar}=FieldData.ListVarName{ivar};
    790                 if test_interp2% interpolate on new grid
    791                     FieldData.(FieldData.ListVarName{ivar})=interp2(FieldData.(AXName),FieldData.(AYName),FieldData.(FieldData.ListVarName{ivar}),AXI,AYI);%TO TEST
    792                 end
    793                 vec_A=reshape(squeeze(FieldData.(FieldData.ListVarName{ivar})),npx*npy,nbcolor); %put the original image in colum
    794                 if nbcolor==1
    795                     vec_B(ind_in)=vec_A(ICOMB);
    796                     vec_B(ind_out)=zeros(size(ind_out));
    797                     A_out=reshape(vec_B,npY,npX);
    798                     ProjData.(FieldData.ListVarName{ivar}) =sum(A_out,1)/npY;
    799                 elseif nbcolor==3
    800                     vec_B(ind_in,1:3)=vec_A(ICOMB,:);
    801                     vec_B(ind_out,1)=zeros(size(ind_out));
    802                     vec_B(ind_out,2)=zeros(size(ind_out));
    803                     vec_B(ind_out,3)=zeros(size(ind_out));
    804                     A_out=reshape(vec_B,npY,npX,nbcolor);
    805                     ProjData.(FieldData.ListVarName{ivar})=squeeze(sum(A_out,1)/npY);
    806                 end 
    807                 ProjData.ListVarName=[ProjData.ListVarName FieldData.ListVarName{ivar}];
    808                 ProjData.VarDimName=[ProjData.VarDimName {AXName}];%to generalize with the initial name of the x coordinate
    809                 ProjData.VarAttribute{ivar}.Role='continuous';% for plot with continuous line
    810             end
    811             if testU
    812                  vector_x =ProjData.(FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_x});
    813                  vector_y =ProjData.(FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_y});
    814                  ProjData.(FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_x}) =cos(theta)*vector_x+sin(theta)*vector_y;
    815                  ProjData.(FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_y}) =-sin(theta)*vector_x+cos(theta)*vector_y;
    816             end
    817             ProjData.VarAttribute{nbvar+1}.long_name='abscissa along line';
    818             if nbcolor==3
    819                 ProjData.VarDimName{end}={AXName,'rgb'};
    820             end
    821         end
    822     elseif strcmp(CellInfo{icell}.CoordType,'tps')
    823          if isfield(ObjectData,'DX')&~isempty(ObjectData.DX)
    824                 DX=abs(ObjectData.DX);%mesh of interpolation points along the line
    825                 Xproj=linelength/(2*npoint):linelength/npoint:linelength-linelength/(2*npoint);
    826                 xreg=cos(theta(ip))*Xproj+ObjectData.Coord(ip,1)
    827                 yreg=sin(theta(ip))*Xproj+ObjectData.Coord(ip,2)
    828 %                 coord_x_proj=XMin:DX:XMax;
    829 %                 coord_y_proj=YMin:DY:YMax;
    830                 DataOut=calc_field_tps(FieldData.FieldList,FieldData,cat(3,xreg,yreg));
    831                 ProjData.ListVarName=[ProjData.ListVarName DataOut.ListVarName];
    832                 ProjData.VarDimName=[ProjData.VarDimName DataOut.VarDimName];
    833                 ProjData.VarAttribute=[ProjData.VarAttribute DataOut.VarAttribute];   
    834                 DataOut.ListVarName(1)=[];
    835                 DataOut.VarDimName(1)=[];
    836                 DataOut.VarAttribute(1)=[];
    837                 for ilist=2:length(DataOut.ListVarName)% reshape data, excluding coordinates (ilist=1-2), TODO: rationalise
    838                     VarName=DataOut.ListVarName{ilist};
    839                      ProjData.(VarName)=DataOut.(VarName);
    840                 end
    841                 ProjData.coord_x=Xproj;
    842          end
    843     end
     779                   DX=DX_min;
     780               end
     781           end
     782           if numel(AY)==2
     783               DY=(AY(2)-AY(1))/(npy-1);
     784           else
     785               DY_vec=diff(AY);
     786               DY=max(DY_vec);
     787               DY_min=min(DY_vec);
     788               if (DY-DY_min)>0.0001*abs(DY)
     789                   test_interp2=1;
     790                   DY=DY_min;
     791               end
     792           end
     793           AXI=linspace(AX(1),AX(end), npx);%set of  x  positions for the interpolated input data
     794           AYI=linspace(AY(1),AY(end), npy);%set of  x  positions for the interpolated input data
     795           if isfield(ObjectData,'DX')
     796               DXY_line=ObjectData.DX;%mesh on the projection line
     797           else
     798               DXY_line=sqrt(abs(DX*DY));% mesh on the projection line
     799           end
     800           dlinx=ObjectData.Coord(2,1)-ObjectData.Coord(1,1);
     801           dliny=ObjectData.Coord(2,2)-ObjectData.Coord(1,2);
     802           linelength=sqrt(dlinx*dlinx+dliny*dliny);
     803           theta=angle(dlinx+i*dliny);%angle of the line
     804           if isfield(FieldData,'RangeX')
     805               XMin=min(FieldData.RangeX);%shift of the origin on the line
     806           else
     807               XMin=0;
     808           end
     809           eval(['ProjData.' AXName '=linspace(XMin,XMin+linelength,linelength/DXY_line+1);'])%abscissa of the new pixels along the line
     810           y=linspace(-width,width,2*width/DXY_line+1);%ordintes of the new pixels (coordinate across the line)
     811           eval(['npX=length(ProjData.' AXName ');'])
     812           npY=length(y); %TODO: utiliser proj_grid
     813           eval(['[X,Y]=meshgrid(ProjData.' AXName ',y);'])%grid in the line coordinates
     814           XIMA=ObjectData.Coord(1,1)+(X-XMin)*cos(theta)-Y*sin(theta);
     815           YIMA=ObjectData.Coord(1,2)+(X-XMin)*sin(theta)+Y*cos(theta);
     816           XIMA=(XIMA-AX(1))/DX+1;%  index of the original image along x
     817           YIMA=(YIMA-AY(1))/DY+1;% index of the original image along y
     818           XIMA=reshape(round(XIMA),1,npX*npY);%indices reorganized in 'line'
     819           YIMA=reshape(round(YIMA),1,npX*npY);
     820           flagin=XIMA>=1 & XIMA<=npx & YIMA >=1 & YIMA<=npy;%flagin=1 inside the original image
     821           ind_in=find(flagin);
     822           ind_out=find(~flagin);
     823           ICOMB=(XIMA-1)*npy+YIMA;
     824           ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A
     825           nbcolor=1; %color images
     826           if numel(npxy)==2
     827               nbcolor=1;
     828           elseif length(npxy)==3
     829               nbcolor=npxy(3);
     830           else
     831               errormsg='multicomponent field not projected';
     832               display(errormsg)
     833               return
     834           end
     835           nbvar=length(ProjData.ListVarName);% number of var from previous cells
     836           ProjData.ListVarName=[ProjData.ListVarName {AXName}];
     837           ProjData.VarDimName=[ProjData.VarDimName {AXName}];
     838           for ivar=VarIndex
     839               %VarName{ivar}=FieldData.ListVarName{ivar};
     840               if test_interp2% interpolate on new grid
     841                   FieldData.(FieldData.ListVarName{ivar})=interp2(FieldData.(AXName),FieldData.(AYName),FieldData.(FieldData.ListVarName{ivar}),AXI,AYI);%TO TEST
     842               end
     843               vec_A=reshape(squeeze(FieldData.(FieldData.ListVarName{ivar})),npx*npy,nbcolor); %put the original image in colum
     844               if nbcolor==1
     845                   vec_B(ind_in)=vec_A(ICOMB);
     846                   vec_B(ind_out)=zeros(size(ind_out));
     847                   A_out=reshape(vec_B,npY,npX);
     848                   ProjData.(FieldData.ListVarName{ivar}) =sum(A_out,1)/npY;
     849               elseif nbcolor==3
     850                   vec_B(ind_in,1:3)=vec_A(ICOMB,:);
     851                   vec_B(ind_out,1)=zeros(size(ind_out));
     852                   vec_B(ind_out,2)=zeros(size(ind_out));
     853                   vec_B(ind_out,3)=zeros(size(ind_out));
     854                   A_out=reshape(vec_B,npY,npX,nbcolor);
     855                   ProjData.(FieldData.ListVarName{ivar})=squeeze(sum(A_out,1)/npY);
     856               end
     857               ProjData.ListVarName=[ProjData.ListVarName FieldData.ListVarName{ivar}];
     858               ProjData.VarDimName=[ProjData.VarDimName {AXName}];%to generalize with the initial name of the x coordinate
     859               ProjData.VarAttribute{ivar}.Role='continuous';% for plot with continuous line
     860           end
     861           if testU
     862               vector_x =ProjData.(FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_x});
     863               vector_y =ProjData.(FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_y});
     864               ProjData.(FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_x}) =cos(theta)*vector_x+sin(theta)*vector_y;
     865               ProjData.(FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_y}) =-sin(theta)*vector_x+cos(theta)*vector_y;
     866           end
     867           ProjData.VarAttribute{nbvar+1}.long_name='abscissa along line';
     868           if nbcolor==3
     869               ProjData.VarDimName{end}={AXName,'rgb'};
     870           end
     871       end
     872   elseif strcmp(CellInfo{icell}.CoordType,'tps')
     873       if isfield(ObjectData,'DX')&~isempty(ObjectData.DX)
     874           DX=abs(ObjectData.DX);%mesh of interpolation points along the line
     875           Xproj=linelength/(2*npoint):linelength/npoint:linelength-linelength/(2*npoint);
     876           xreg=cos(theta(ip))*Xproj+ObjectData.Coord(ip,1)
     877           yreg=sin(theta(ip))*Xproj+ObjectData.Coord(ip,2)
     878           %                 coord_x_proj=XMin:DX:XMax;
     879           %                 coord_y_proj=YMin:DY:YMax;
     880           DataOut=calc_field_tps(FieldData.FieldList,FieldData,cat(3,xreg,yreg));
     881           ProjData.ListVarName=[ProjData.ListVarName DataOut.ListVarName];
     882           ProjData.VarDimName=[ProjData.VarDimName DataOut.VarDimName];
     883           ProjData.VarAttribute=[ProjData.VarAttribute DataOut.VarAttribute];
     884           DataOut.ListVarName(1)=[];
     885           DataOut.VarDimName(1)=[];
     886           DataOut.VarAttribute(1)=[];
     887           for ilist=2:length(DataOut.ListVarName)% reshape data, excluding coordinates (ilist=1-2), TODO: rationalise
     888               VarName=DataOut.ListVarName{ilist};
     889               ProjData.(VarName)=DataOut.(VarName);
     890           end
     891           ProjData.coord_x=Xproj;
     892       end
     893   end
    844894end
    845895
Note: See TracChangeset for help on using the changeset viewer.