Ignore:
Timestamp:
Apr 17, 2020, 5:58:49 PM (4 years ago)
Author:
sommeria
Message:

ima2temperature.m added and various bug repairs

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/proj_field.m

    r1078 r1080  
    807807           
    808808            if max(NbDim)==3 % 3D case
    809                 AZName=FieldData.ListVarName{CellInfo{icell}.CoordIndex(end)};
    810                 AYName=FieldData.ListVarName{CellInfo{icell}.CoordIndex(end-1)};
    811                 AXName=FieldData.ListVarName{CellInfo{icell}.CoordIndex(end-2)};
    812                 AX=FieldData.(AXName);% set of x positions
    813                 AY=FieldData.(AYName);% set of y positions
    814                 AZ=FieldData.(AZName);% set of z positions
    815                 AName=FieldData.ListVarName{VarIndex(1)};
    816                 npxy=size(FieldData.(AName));
    817                 npz=npxy(1);
    818                 npy=npxy(2);
    819                 npx=npxy(1);
    820                 AXI=linspace(AX(1),AX(end), npx);%set of  x  positions for the interpolated input data
    821                 AYI=linspace(AY(1),AY(end), npy);%set of  x  positions for the interpolated input data
    822                 AZI=linspace(AZ(1),AZ(end), npy);%set of  x  positions for the interpolated input data
     809                AXName=FieldData.ListVarName{CellInfo{icell}.CoordIndex(3)};
     810                Coord{1}=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(1)});%initial z coordinates
     811                Coord{2}=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(2)});%initial y coordinates
     812                Coord{3}=FieldData.(AXName);%initial x coordinates
     813                if size(ObjectData.Coord,2)<3
     814                    ObjectData.Coord(1,3)=0;
     815                end
     816                dline=ObjectData.Coord(2,:)-ObjectData.Coord(1,:);
     817                linelength=norm(dline);
     818                if isfield(FieldData,'RangeX')
     819                    XMin=min(FieldData.RangeX);%shift of the origin on the line
     820                else
     821                    XMin=0;
     822                end
     823                ProjData.(AXName)=XMin:ObjectData.DX:XMin+linelength;%abscissa of the projected data along the line
     824               
     825                XI_proj=ObjectData.Coord(1,1)*ones(size(ProjData.(AXName)));
     826                YI_proj=ObjectData.Coord(1,2)*ones(size(ProjData.(AXName)));
     827                ZI_proj=ObjectData.Coord(1,3)*ones(size(ProjData.(AXName)));
     828                if dline(1,1)~=0
     829                    XI_proj=XI_proj+ProjData.(AXName)/dline(1,1);
     830                end
     831                if dline(1,2)~=0
     832                    YI_proj=YI_proj+ProjData.(AXName)/dline(1,2);
     833                end
     834                if dline(1,3)~=0
     835                    ZI_proj=ZI_proj+ProjData.(AXName)/dline(1,3);
     836                end
    823837                for ivar=VarIndex
    824                     VarName=FieldData.ListVarName{ivar};
    825                     FieldData.(VarName)=interp3(FieldData.(AXName),FieldData.(AYName),FieldData.(AZName),FieldData.(VarName),AXI,AYI,AZI);
    826                    
    827                     %                     vec_A=reshape(squeeze(FieldData.(FieldData.ListVarName{ivar})),npx*npy,nbcolor); %put the original image in colum
    828                     %                     if nbcolor==1
    829                     %                         vec_B(ind_in)=vec_A(ICOMB);
    830                     %                         vec_B(ind_out)=zeros(size(ind_out));
    831                     %                         A_out=reshape(vec_B,npY,npX);
    832                     %                         ProjData.(FieldData.ListVarName{ivar}) =sum(A_out,1)/npY;
    833                     %                     elseif nbcolor==3
    834                     %                         vec_B(ind_in,1:3)=vec_A(ICOMB,:);
    835                     %                         vec_B(ind_out,1)=zeros(size(ind_out));
    836                     %                         vec_B(ind_out,2)=zeros(size(ind_out));
    837                     %                         vec_B(ind_out,3)=zeros(size(ind_out));
    838                     %                         A_out=reshape(vec_B,npY,npX,nbcolor);
    839                     %                         ProjData.(FieldData.ListVarName{ivar})=squeeze(sum(A_out,1)/npY);
    840                     %                     end
     838                                      VarName=FieldData.ListVarName{ivar};
     839                    %                         ListVarName=[ListVarName VarName];
     840                    %                         VarDimName=[VarDimName {{'coord_y','coord_x'}}];
     841                    %                         VarAttribute{length(ListVarName)}=FieldData.VarAttribute{ivar}; %reproduce the variable attributes
    841842                    ProjData.ListVarName=[ProjData.ListVarName FieldData.ListVarName{ivar}];
    842843                    ProjData.VarDimName=[ProjData.VarDimName {AXName}];%to generalize with the initial name of the x coordinate
    843844                    ProjData.VarAttribute{ivar}.Role='continuous';% for plot with continuous line
    844                 end
    845                
     845                   
     846                    FieldData.(VarName)=permute(FieldData.(VarName),[2 3 1]); %coordinate permutation needed to use interp3
     847                    indexnan=isnan(FieldData.(VarName));
     848                    FieldData.(VarName)(indexnan)=0;%set to zero the undefined values
     849                    ProjData.(VarName)=interp3(Coord{3},Coord{2},Coord{1},double(FieldData.(VarName)),XI_proj,YI_proj,ZI_proj,'*linear');
     850                    ProjData.(VarName)=squeeze(ProjData.(VarName));
     851                end           
     852           
    846853            else
    847854                AYName=FieldData.ListVarName{CellInfo{icell}.CoordIndex(end-1)};
    … …  
    959966                    ProjData.VarDimName{end}={AXName,'rgb'};
    960967                end
    961             end
    962            
    963     end
    964     % for vector fields, take the components longitudinal and tranverse to the projection line
    965     if ~isempty(ivar_U) && ~isempty(ivar_V)
    966         vector_x =ProjData.(ProjData.ListVarName{ivar_U});
    967         ProjData.(ProjData.ListVarName{ivar_U}) =cos(theta)*vector_x+sin(theta)*ProjData.(ProjData.ListVarName{ivar_V});
    968         ProjData.(ProjData.ListVarName{ivar_V}) =-sin(theta)*vector_x+cos(theta)*ProjData.(ProjData.ListVarName{ivar_V});
    969     end
     968    end
     969   
     970end
     971% for vector fields, take the components longitudinal and tranverse to the projection line
     972if ~isempty(ivar_U) && ~isempty(ivar_V)
     973    vector_x =ProjData.(ProjData.ListVarName{ivar_U});
     974    ProjData.(ProjData.ListVarName{ivar_U}) =cos(theta)*vector_x+sin(theta)*ProjData.(ProjData.ListVarName{ivar_V});
     975    ProjData.(ProjData.ListVarName{ivar_V}) =-sin(theta)*vector_x+cos(theta)*ProjData.(ProjData.ListVarName{ivar_V});
     976end
    970977end
    971978
    … …  
    10021009
    10031010%% rotation matrix
    1004 PlaneAngle=[0 0];
     1011PlaneAngle=[0 0 0];
    10051012norm_plane=[0 0 1];
    1006 test90x=0;%=1 for 90 degree rotation alround x axis
    1007 test90y=0;%=1 for 90 degree rotation alround y axis
    1008 
    1009 if isfield(ObjectData,'Angle')
    1010     checkM2=0;
    1011     if numel(ObjectData.Angle)==2 && ~isequal(ObjectData.Angle,[0 0])
    1012         checkM2=1;
    1013         test90y=0;%isequal(ObjectData.Angle,[0 90 0]);
    1014     end
    1015     PlaneAngle=(pi/180)*ObjectData.Angle;
    1016     if PlaneAngle==0
    1017         PlaneAngle=[0 0];
    1018     end   
    1019     M1=[cos(PlaneAngle(1)) -sin(PlaneAngle(1)) 0;sin(PlaneAngle(1)) cos(PlaneAngle(1)) 0;0 0 1];
    1020     M=M1;
    1021     if checkM2
    1022         M2=[1 0 0;0 cos(PlaneAngle(2)) -sin(PlaneAngle(2));0 sin(PlaneAngle(2)) cos(PlaneAngle(2))];
    1023         M=M1*M2;% first rotate in the x,y plane with angle PlaneAngle(1), then slant around the new x axis0 with angle PlaneAngle(2)
    1024     end
     1013testangle=0;
     1014test90x=0;
     1015test90y=0;
     1016if isfield(ObjectData,'Angle') && size(ObjectData.Angle,2)==3
     1017    test90x=isequal(ObjectData.Angle,[90 0 0]);%=1 for 90 degree rotation alround x axis
     1018    test90y=isequal(ObjectData.Angle,[0 90 0]);%=1 for 90 degree rotation alround y axis
     1019    %test90z=isequal(PlaneAngle,[90 0 0]);%=1 for 90 degree rotation alround x axis
     1020    PlaneAngle=(pi/180)*ObjectData.Angle;
     1021    M=rodrigues(PlaneAngle);
    10251022    norm_plane=M*[0 0 1]'; 
    1026 end
    1027 testangle=~isequal(PlaneAngle,[0 0])||~isequal(ObjectData.Coord(1:2),[0 0 ]) ;% && ~test90y && ~test90x;%=1 for slanted plane
     1023    testangle= (~isequal(PlaneAngle,[0 0 0]) && ~test90x && ~test90y);%=1 for slanted plane
     1024end
    10281025
    10291026%% mesh sizes DX and DY
    … …  
    12351232    VarAttribute={};% initiate coresponding list of var attributes  for cell # icell
    12361233    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    1237     switch CellInfo{icell}.CoordType
    1238         
     1234    check3D=(numel(CellInfo{icell}.CoordIndex)==3);
     1235    switch CellInfo{icell}.CoordType   
    12391236        case 'scattered'
    12401237            %% case of input fields with unstructured coordinates (applies for projMode ='projection' or 'interp_lin')
    … …  
    12441241            coord_x=FieldData.(CellInfo{icell}.XName);% initial x coordinates
    12451242            coord_y=FieldData.(CellInfo{icell}.YName);% initial y coordinates
    1246             check3D=(numel(CellInfo{icell}.CoordIndex)==3);
     1243     
    12471244            if check3D
    12481245                coord_z=FieldData.(CellInfo{icell}.ZName);
    … …  
    12711268           
    12721269            %rotate coordinates if needed: coord_X,coord_Y= = coordinates in the new plane
    1273             Phi=PlaneAngle(1);
    1274             %             Theta=PlaneAngle(2);
    1275             % Phi=PlaneAngle(3);
    1276             if testangle && ~test90y && ~test90x;%=1 for slanted plane
    1277                 coord_X=coord_x *cos(Phi) + coord_y* sin(Phi);
    1278                 coord_Y=-coord_x *sin(Phi) + coord_y *cos(Phi);
    1279                 if checkM2
    1280                     coord_Y=coord_Y*cos(PlaneAngle(2));
    1281                     coord_Y=coord_Y+coord_z *sin(PlaneAngle(2));
    1282                     coord_X=(coord_X *cos(Psi) - coord_Y* sin(Psi));%A VERIFIER
    1283                     coord_Y=(coord_X *sin(Psi) + coord_Y* cos(Psi));
     1270           
     1271            Phi=PlaneAngle(3);
     1272            if testangle
     1273                if check3D
     1274                    [coord_X,coord_Y]=rotate_vector(PlaneAngle,coord_x,coord_y,coord_z);
     1275                else
     1276                    coord_X=coord_x *cos(Phi) + coord_y* sin(Phi);
     1277                    coord_Y=-coord_x *sin(Phi) + coord_y *cos(Phi);
    12841278                end
    12851279            else
    … …  
    14151409               
    14161410                %rotate coordinates if needed: coord_X,coord_Y= = coordinates in the new plane
    1417                 Phi=PlaneAngle(1);
     1411                Phi=PlaneAngle(3);
    14181412                if testangle && ~test90y && ~test90x %=1 for slanted plane
    14191413                    new_XI=XI *cos(Phi) - YI* sin(Phi)+ObjectData.Coord(1);
    … …  
    16881682                    % determine the boundaries of the projected field,
    16891683                    % first find the 8 summits of the initial volume in the
    1690                     PlaneAngle=[0 0 0];% default
    1691                     PlaneAngle(1:numel(ObjectData.Angle))=ObjectData.Angle*pi/180;
     1684%                     PlaneAngle=[0 0 0];% default
     1685%                     PlaneAngle(1:numel(ObjectData.Angle))=ObjectData.Angle*pi/180;
    16921686                    % new coordinates
    16931687                    Coord{1}=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(1)});%initial z coordinates
    16941688                    Coord{2}=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(2)});%initial y coordinates
    16951689                    Coord{3}=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(3)});%initial x coordinates
    1696 %                     summit=zeros(3,8);% initialize summit coordinates
    1697 %                     summit(1,1:4)=[Coord{3}(1) Coord{3}(end) Coord{3}(1) Coord{3}(end)];%square
    1698 %                     summit(2,1:4)=[Coord{2}(1) Coord{2}(1) Coord{2}(end) Coord{2}(end)];% square at z= Coord{1}(1)
    1699 %                     summit(1:2,5:8)=summit(1:2,1:4);
    1700 %                     summit(3,:)=[Coord{1}(1)*ones(1,4) Coord{1}(end)*ones(1,4)];
    1701 %                     %Mrot_inv=rodrigues(-PlaneAngle);
    1702 %                     newsummit=zeros(3,8);% initialize the rotated summit coordinates
    1703                     ObjectData.Coord=ObjectData.Coord';% set ObjectData.Coord as a vertical vector
    1704                     if size(ObjectData.Coord,1)<3
    1705                         ObjectData.Coord=[ObjectData.Coord; 0];%add z origin at z=0 by default
    1706                     end
     1690
     1691                    coord_x_proj=ObjectData.RangeX(1):InterpMesh:ObjectData.RangeX(2);% set of coordinates in the projection plane
     1692                    coord_y_proj=ObjectData.RangeY(1):InterpMesh:ObjectData.RangeY(2);
     1693                    %coord_z_proj=-floor(ObjectData.RangeInterp/InterpMesh):InterpMesh:floor(ObjectData.RangeInterp/InterpMesh);
     1694                    M=rodrigues(ObjectData.Angle);
     1695                    [XI,YI]=meshgrid(coord_x_proj,coord_y_proj);
     1696                    XI_proj=M(1,1)*XI+M(2,1)*YI+ObjectData.Coord(1,1);
     1697                    YI_proj=M(2,1)*XI+M(2,2)*YI+ObjectData.Coord(1,2);
     1698                    ZI_proj=M(3,1)*XI+M(3,2)*YI+ObjectData.Coord(1,3);
    17071699                   
    1708 %                     M1=[cos(PlaneAngle(1)) sin(PlaneAngle(1)) 0;-sin(PlaneAngle(1)) cos(PlaneAngle(1)) 0;0 0 1];
    1709 %                     M2=[1 0 0;0 cos(PlaneAngle(2)) sin(PlaneAngle(2));0 -sin(PlaneAngle(2)) cos(PlaneAngle(2))];
    1710 %                     M=M1*M2;
    1711                     M_inv=inv(M);
    1712                    
    1713 %                     for isummit=1:8% TODO: introduce a function for rotation of n points (to use also for scattered data)
    1714 %                         newsummit(:,isummit)=M*(summit(:,isummit)-(ObjectData.Coord));
    1715 %                     end
    1716 %                     coord_x_proj=min(newsummit(1,:)):InterpMesh: max(newsummit(1,:));% set of coordinqtes in the projection plane
    1717 %                     coord_y_proj=min(newsummit(2,:)):InterpMesh: max(newsummit(2,:));
    1718 %                     coord_z_proj=-width:width;
    1719                     coord_x_proj=ObjectData.RangeX(1):InterpMesh:ObjectData.RangeX(2);% set of coordinqtes in the projection plane
    1720                      coord_y_proj=ObjectData.RangeY(1):InterpMesh:ObjectData.RangeY(2);
    1721                     coord_z_proj=-floor(ObjectData.RangeInterp/InterpMesh):InterpMesh:floor(ObjectData.RangeInterp/InterpMesh);
    1722                     %Mrot=rodrigues(PlaneAngle);% inverse rotation matrix
    1723                     %Origin=M_inv*[coord_x_proj(1);coord_y_proj(1);coord_z_proj(1)]+ObjectData.Coord;
    1724                     npx=numel(coord_x_proj);
    1725                     npy=numel(coord_y_proj);
    1726                     npz=numel(coord_z_proj);
    1727                    
    1728                     %modangle=sqrt(PlaneAngle(1)*PlaneAngle(1)+PlaneAngle(2)*PlaneAngle(2));
    1729                     %                     cosphi=PlaneAngle(1)/modangle;
    1730                     %                     sinphi=PlaneAngle(2)/modangle;
    1731                     iX=[coord_x_proj(end)-coord_x_proj(1);0;0]/(npx-1);
    1732                     iY=[0;coord_y_proj(end)-coord_y_proj(1);0]/(npy-1);
    1733                     if npz==1
    1734                         iZ=[0;0;0];
    1735                     else
    1736                     iZ=[0;0;coord_z_proj(end)-coord_z_proj(1)]/(npz-1);
    1737                     end
    1738                     %                     iX(1:2)=[cosphi -sinphi;sinphi cosphi]*iX(1:2);
    1739                     %                     iY(1:2)=[-cosphi -sinphi;sinphi cosphi]*iY(1:2);
    1740                    
    1741                     ix=M_inv*iX;%  vector along the new x coordinates transformed into old coordinates
    1742                     iy=M_inv*iY;% vector along y coordinates
    1743                     iz=M_inv*iZ;% vector along z coordinates
    1744                    
    1745                     [Grid_x,Grid_y,Grid_z]=meshgrid(0:npx-1,0:npy-1,0:npz-1);
    1746                     %[Grid_x,Grid_y,Grid_z]=meshgrid(0:npz-1,0:npy-1,0:npx-1);
    1747                     if ismatrix(Grid_x)% add a singleton in case of a single z value
    1748                         Grid_x=shiftdim(Grid_x,-1);
    1749                         Grid_y=shiftdim(Grid_y,-1);
    1750                         Grid_z=shiftdim(Grid_z,-1);
    1751                     end
    1752                     XI=ObjectData.Coord(1)+ix(1)*Grid_x+iy(1)*Grid_y+iz(1)*Grid_z;
    1753                     YI=ObjectData.Coord(2)+ix(2)*Grid_x+iy(2)*Grid_y+iz(2)*Grid_z;
    1754                     ZI=ObjectData.Coord(3)+ix(3)*Grid_x+iy(3)*Grid_y+iz(3)*Grid_z;
    1755                   %  [X,Y,Z]=meshgrid(Coord{3},Coord{2},Coord{1});% grid of initial coordinates
     1700   
    17561701                    for ivar=VarIndex
    17571702                        VarName=FieldData.ListVarName{ivar};
    … …  
    17641709                        ProjData.coord_x=coord_x_proj;
    17651710                        ProjData.coord_y=coord_y_proj;
    1766                         ProjData.(VarName)=interp3(Coord{3},Coord{2},Coord{1},double(FieldData.(VarName)),ZI,YI,XI,'*linear');
    1767                         ProjData.(VarName)=nanmean(ProjData.(VarName),1);
     1711                        ProjData.(VarName)=interp3(Coord{3},Coord{2},Coord{1},double(FieldData.(VarName)),XI_proj,YI_proj,ZI_proj,'*linear');
     1712                        ProjData.(VarName)=nanmean(ProjData.(VarName),3);
    17681713                        ProjData.(VarName)=squeeze(ProjData.(VarName));
    17691714                    end
    … …  
    17971742                UName=ListVarName{ivar_U};
    17981743                VName=ListVarName{ivar_V};
    1799                 if checkM2
     1744                if check3D
    18001745                    UValue=cos(PlaneAngle(1))*ProjData.(UName)+ sin(PlaneAngle(1))*ProjData.(VName);
    18011746                    ProjData.(VName)=(-sin(PlaneAngle(1))*ProjData.(UName)+ cos(PlaneAngle(1))*ProjData.(VName));
Note: See TracChangeset for help on using the changeset viewer.