Changeset 959 for trunk/src/proj_field.m


Ignore:
Timestamp:
Jun 24, 2016, 8:49:08 PM (5 years ago)
Author:
sommeria
Message:

various

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/proj_field.m

    r955 r959  
    16531653                    Mrot=rodrigues(PlaneAngle);
    16541654                    newsummit=zeros(3,8);% initialize the rotated summit coordinates
    1655                     ObjectData.Coord=[ObjectData.Coord(1,1); ObjectData.Coord(1,2); 0];%TODO: set in set_object
     1655                    ObjectData.Coord=ObjectData.Coord';
     1656                    if size(ObjectData.Coord,2)<3
     1657                    ObjectData.Coord=[ObjectData.Coord; 0];%add z origin at z=0 by default
     1658                    end
    16561659                    for isummit=1:8% TODO: introduce a function for rotation of n points (to use also for scattered data)
    16571660                        newsummit(:,isummit)=Mrot*(summit(:,isummit)-(ObjectData.Coord));
    16581661                    end
    1659                     coord_x_proj=min(newsummit(1,:)):InterpMesh: max(newsummit(1,:));
     1662                    coord_x_proj=min(newsummit(1,:)):InterpMesh: max(newsummit(1,:));% set of coordinqtes in the projection plane
    16601663                    coord_y_proj=min(newsummit(2,:)):InterpMesh: max(newsummit(2,:));
    1661                     coord_z_proj=-width:InterpMesh:width;
    1662                     Mrot_inverse=rodrigues(-PlaneAngle);
     1664                    coord_z_proj=-width:width;
     1665                    Mrot_inverse=rodrigues(-PlaneAngle);% inverse rotation matrix
    16631666                    Origin=Mrot_inverse*[coord_x_proj(1);coord_y_proj(1);coord_z_proj(1)]+ObjectData.Coord;
    1664                     ix=Mrot_inverse*[coord_x_proj(end)-coord_x_proj(1);0;0];% unit vector along the new x coordinates
    1665                     iy=Mrot_inverse*[0;coord_y_proj(end)-coord_y_proj(1);0];% unit vector y coordinates
    1666                     iz=Mrot_inverse*[0;0;coord_z_proj(end)-coord_z_proj(1)];% x unit vector z coordinates
    1667                     [Grid_x,Grid_y,Grid_z]=meshgrid(1:numel(coord_x_proj),1:numel(coord_y_proj),1:numel(coord_z_proj));
    1668                     if ismatrix(Grid_x)
     1667                    npx=numel(coord_x_proj);
     1668                    npy=numel(coord_y_proj);
     1669                    npz=numel(coord_z_proj);
     1670                    ix=Mrot_inverse*[coord_x_proj(end)-coord_x_proj(1);0;0]/(npx-1);% unit vector along the new x coordinates
     1671                    iy=Mrot_inverse*[0;coord_y_proj(end)-coord_y_proj(1);0]/(npy-1);% unit vector y coordinates
     1672                    iz=Mrot_inverse*[0;0;coord_z_proj(end)-coord_z_proj(1)]/(npz-1);% x unit vector z coordinates
     1673                    [Grid_x,Grid_y,Grid_z]=meshgrid(0:npx-1,0:npy-1,0:npz-1);
     1674                    if ismatrix(Grid_x)% add a singleton in case of a single z value
    16691675                        Grid_x=shiftdim(Grid_x,-1);
    16701676                        Grid_y=shiftdim(Grid_y,-1);
     
    16811687                            VarName=FieldData.ListVarName{ivar};
    16821688                            ListVarName=[ListVarName VarName];
     1689                            VarDimName=[VarDimName {{'coord_y','coord_x'}}];
    16831690                            VarAttribute{length(ListVarName)}=FieldData.VarAttribute{ivar}; %reproduce the variable attributes
    16841691                            FieldData.(VarName)=permute(FieldData.(VarName),[2 3 1]);
    1685                             ProjData.(VarName)=interp3(X,Y,Z,double(FieldData.(VarName)),XI,YI,ZI);
     1692                            ProjData.coord_x=coord_x_proj;
     1693                            ProjData.coord_y=coord_y_proj;
     1694                            ProjData.(VarName)=interp3(X,Y,Z,double(FieldData.(VarName)),XI,YI,ZI,'*linear');
     1695                            ProjData.(VarName)=nanmean(ProjData.(VarName),3);
     1696                            ProjData.(VarName)=squeeze(ProjData.(VarName));
    16861697                    end
    16871698                end
Note: See TracChangeset for help on using the changeset viewer.