Changeset 644 for trunk/src/proj_field.m
- Timestamp:
- May 28, 2013, 11:30:01 PM (11 years ago)
- File:
-
- 1 edited
-
trunk/src/proj_field.m (modified) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/src/proj_field.m
r630 r644 564 564 VarIndex=find(check_proj); 565 565 566 % identify vector components566 %% identify vector components 567 567 testU=isfield(CellInfo{icell},'VarIndex_vector_x') &&isfield(CellInfo{icell},'VarIndex_vector_y') ;% test for vectors 568 568 if testU … … 578 578 errorflag=FieldData.(FFName); 579 579 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 581 582 if strcmp(CellInfo{icell}.CoordType,'scattered') 582 if ~isequal(ProjMode,'interp_lin')583 if strcmp(ProjMode,'projection') 583 584 if width==0 584 585 errormsg='range of the projection object is missing'; 585 return 586 e lse587 lambda=2/(width*width); %smoothing factor used for interp_tps: weight exp(-2) at distance width from the line588 end 589 end590 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) 592 593 DX=abs(ObjectData.DX);%mesh of interpolation points along the line 593 594 else … … 601 602 coord_y=FieldData.(YName); 602 603 end 603 %initiate projection 604 605 %% initiate projection 604 606 for ivar=1:length(VarIndex) 605 607 ProjLine{ivar}=[]; … … 612 614 %%%%%%% % A FAIRE CALCULER MEAN DES QUANTITES %%%%%% 613 615 %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) 739 778 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 844 894 end 845 895
Note: See TracChangeset
for help on using the changeset viewer.
