Changeset 648 for trunk/src/proj_field.m
 Timestamp:
 Jun 9, 2013, 10:31:58 PM (10 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/src/proj_field.m
r646 r648 84 84 ProjData=[]; 85 85 86 %% in the absence of object Type or projection mode, or object coordinaes, the output is empty86 %% check input projection object: type, projection mode and Coord: 87 87 if ~isfield(ObjectData,'Type')~isfield(ObjectData,'ProjMode') 88 88 return 89 89 end 90 % case of no projection (object is used only as graph display) 91 if ise qual(ObjectData.ProjMode,'none')isequal(ObjectData.ProjMode,'mask_inside')isequal(ObjectData.ProjMode,'mask_outside')90 ListProjMode={'projection','interp_lin','interp_tps'};%list of effective projection modes 91 if isempty(strcmp(ObjectData.ProjMode,ListProjMode)) 92 92 return 93 93 end … … 103 103 switch ObjectData.Type 104 104 case 'points' 105 [ProjData,errormsg]=proj_points(FieldData,ObjectData);105 [ProjData,errormsg]=proj_points(FieldData,ObjectData); 106 106 case {'line','polyline'} 107 [ProjData,errormsg] = proj_line(FieldData,ObjectData);107 [ProjData,errormsg] = proj_line(FieldData,ObjectData); 108 108 case {'polygon','rectangle','ellipse'} 109 109 if isequal(ObjectData.ProjMode,'inside')isequal(ObjectData.ProjMode,'outside') … … 113 113 end 114 114 case 'plane' 115 115 [ProjData,errormsg] = proj_plane(FieldData,ObjectData); 116 116 case 'volume' 117 117 [ProjData,errormsg] = proj_volume(FieldData,ObjectData); … … 344 344 %LOOP ON GROUPS OF VARIABLES SHARING THE SAME DIMENSIONS 345 345 for icell=1:length(CellInfo) 346 testX=0;347 testY=0;346 CoordType=CellInfo{icell}.CoordType; 347 % testY=0; 348 348 test_Amat=0; 349 349 if NbDim(icell)~=2% proj_patch acts only on fields of space dimension 2 … … 429 429 end 430 430 elseif isequal(ObjectData.Type,'polygon') 431 if testX431 if strcmp(CoordType,'scattered') 432 432 testin=inpolygon(coord_x,coord_y,ObjectData.Coord(:,1),ObjectData.Coord(:,2)); 433 elseif test_Amat433 elseif strcmp(CoordType,'grid') 434 434 testin=inpolygon(Xi,Yi,ObjectData.Coord(:,1),ObjectData.Coord(:,2)); 435 435 else%calculate the scalar … … 439 439 X2Max=widthx*widthx; 440 440 Y2Max=(widthy)*(widthy); 441 if testX441 if strcmp(CoordType,'scattered') 442 442 distX=(coord_xObjectData.Coord(1,1)); 443 443 distY=(coord_yObjectData.Coord(1,2)); 444 444 testin=(distX.*distX/X2Max+distY.*distY/Y2Max)<1; 445 elseif test_Amat%case of usual 2x2 matrix445 elseif strcmp(CoordType,'grid') %case of usual 2x2 matrix 446 446 distX=(XiObjectData.Coord(1,1)); 447 447 distY=(YiObjectData.Coord(1,2)); … … 495 495 ProjData.NbDim=1; 496 496 %initialisation of the input parameters and defaultoutput 497 ProjMode=ObjectData.ProjMode; 497 ProjMode=ObjectData.ProjMode; %rmq: ProjMode always defined from input={'projection','interp_lin','interp_tps'} 498 498 % 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 499 width=0; 504 500 if isfield(ObjectData,'RangeY') 505 width=max(ObjectData.RangeY); 506 end 507 501 width=max(ObjectData.RangeY);%Rangey needed bfor mode 'projection' 502 end 508 503 % default output 509 errormsg= [];%default504 errormsg='';%default 510 505 Xline=[]; 511 506 flux=0; 512 507 circul=0; 513 508 liny=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 509 NbPoints=size(ObjectData.Coord,1); 518 510 testfalse=0; 519 511 ListIndex={}; 520 512 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': 515 LineCoord=ObjectData.Coord; 516 switch 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); 527 end 528 if ~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); 537 end 538 CheckClosedLine=~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 526 546 % 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(ip1))/2)/cos((theta(ip1)theta(ip))/2); 535 xinf(ip)=ObjectData.Coord(ip,1)+width*sin((theta(ip)+theta(ip1))/2)/cos((theta(ip1)theta(ip))/2); 536 ysup(ip)=ObjectData.Coord(ip,2)+width*cos((theta(ip)+theta(ip1))/2)/cos((theta(ip1)theta(ip))/2); 537 yinf(ip)=ObjectData.Coord(ip,2)width*cos((theta(ip)+theta(ip1))/2)/cos((theta(ip1)theta(ip))/2); 538 end 539 end 540 % 541 % 542 % x = a \ \cosh \mu \ \cos \nu 543 % 544 % y = a \ \sinh \mu \ \sin \nu 547 xsup=zeros(1,NbPoints); xinf=zeros(1,NbPoints); ysup=zeros(1,NbPoints); yinf=zeros(1,NbPoints); 548 if 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 558 else % 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=(LineLengthDX*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 596 end 545 597 546 598 … … 551 603 return 552 604 end 553 554 %% loop on variable cells with the same space dimension 605 CellInfo=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 555 609 ProjData.ListVarName={}; 556 610 ProjData.VarDimName={}; 557 611 for 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 573 623 testU=isfield(CellInfo{icell},'VarIndex_vector_x') &&isfield(CellInfo{icell},'VarIndex_vector_y') ;% test for vectors 574 624 if testU … … 577 627 vector_x=FieldData.(UName); 578 628 vector_y=FieldData.(VName); 579 end 629 end 580 630 %identify error flag 581 testfalse=isfield(CellInfo{icell},'VarIndex_errorflag');% test forerror flag582 if testfalse631 errorflag=0; %default, no error flag 632 if isfield(CellInfo{icell},'VarIndex_errorflag');% test for error flag 583 633 FFName=FieldData.ListVarName{CellInfo{icell}.VarIndex_errorflag}; 584 634 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 586 638 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(end1)}; 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_xxinf(1))*(yinf(2)yinf(1))) ... 656 & ((coord_y ysup(1))*(xsup(2)xsup(1))<(coord_xxsup(1))*(ysup(2)ysup(1))) ... 657 & ((coord_y yinf(2))*(xsup(2)xinf(2))>(coord_xxinf(2))*(ysup(2)yinf(2))) ... 658 & ((coord_y yinf(1))*(xsup(1)xinf(1))<(coord_xxinf(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_xObjectData.Coord(1,1))*costheta + (coord_yObjectData.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']; % 600 722 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(end1)}; 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_xxinf(ip))*(yinf(ip+1)yinf(ip))) ... 633 & ((coord_y ysup(ip))*(xsup(ip+1)xsup(ip))<(coord_xxsup(ip))*(ysup(ip+1)ysup(ip))) ... 634 & ((coord_y yinf(ip+1))*(xsup(ip+1)xinf(ip+1))>(coord_xxinf(ip+1))*(ysup(ip+1)yinf(ip+1))) ... 635 & ((coord_y yinf(ip))*(xsup(ip)xinf(ip))<(coord_xxinf(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_selObjectData.Coord(ip,1))*costheta + (Y_selObjectData.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:linelengthlinelength/(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:linelengthlinelength/(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*((xijxregij).*(xijxregij)+(yijyregij).*(yijyregij))); 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(end1)}; 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))/(npx1); 779 else 780 DX_vec=diff(AX); 781 DX=max(DX_vec); 782 DX_min=min(DX_vec); 783 if (DXDX_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))/(npy1); 790 else 791 DY_vec=diff(AY); 792 DY=max(DY_vec); 793 DY_min=min(DY_vec); 794 if (DYDY_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)+(XXMin)*cos(theta)Y*sin(theta); 821 YIMA=ObjectData.Coord(1,2)+(XXMin)*sin(theta)+Y*cos(theta); 822 XIMA=(XIMAAX(1))/DX+1;% index of the original image along x 823 YIMA=(YIMAAY(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=(XIMA1)*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:linelengthlinelength/(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=12), 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(end1)}; 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))/(npx1); 736 else 737 DX_vec=diff(AX); 738 DX=max(DX_vec); 739 DX_min=min(DX_vec); 740 if (DXDX_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))/(npy1); 747 else 748 DY_vec=diff(AY); 749 DY=max(DY_vec); 750 DY_min=min(DY_vec); 751 if (DYDY_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)+(XXMin)*cos(theta)Y*sin(theta); 778 YIMA=ObjectData.Coord(1,2)+(XXMin)*sin(theta)+Y*cos(theta); 779 XIMA=(XIMAAX(1))/DX+1;% index of the original image along x 780 YIMA=(YIMAAY(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=(XIMA1)*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 900 836 end 901 837
Note: See TracChangeset
for help on using the changeset viewer.