Changeset 954 for trunk/src/proj_field.m


Ignore:
Timestamp:
Jun 22, 2016, 1:36:50 PM (5 years ago)
Author:
sommeria
Message:

update calib modfied + various updates

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/proj_field.m

    r937 r954  
    9090    return
    9191end
    92 ListProjMode={'projection','interp_lin','interp_tps','inside','outside'};%list of effective projection modes
    93 if isempty(find(strcmp(ObjectData.ProjMode,ListProjMode), 1))% no projection in case
     92% check list of effective projection modes
     93if ~ismember(ObjectData.ProjMode,{'projection','interp_lin','interp_tps','inside','outside'})
    9494    return
    9595end
     
    117117            [ProjData,errormsg] = proj_line(FieldData,ObjectData);
    118118        end
    119     case 'plane'
     119    case {'plane','plane_z'}
    120120        [ProjData,errormsg] = proj_plane(FieldData,ObjectData);
    121121    case 'volume'
     
    955955test90x=0;%=1 for 90 degree rotation alround x axis
    956956test90y=0;%=1 for 90 degree rotation alround y axis
     957if strcmp(ObjectData.Type,'plane_z')
     958    Delta_x=ObjectData.Coord(2,1)-ObjectData.Coord(1,1);
     959    Delta_y=ObjectData.Coord(2,2)-ObjectData.Coord(1,2);
     960    Delta_mod=sqrt(Delta_x*Delta_x+Delta_y*Delta_y);
     961    ObjectData.Angle=[0 0 0];
     962    ObjectData.Angle(1)=90*Delta_x/Delta_mod;
     963    ObjectData.Angle(2)=90*Delta_y/Delta_mod;
     964end   
    957965if isfield(ObjectData,'Angle')&& isequal(size(ObjectData.Angle),[1 3])&& ~isequal(ObjectData.Angle,[0 0 0])
    958966    test90y=isequal(ObjectData.Angle,[0 90 0]);
     
    987995    return
    988996end
     997InterpMesh=min(DX,DY);%mesh used for interpolation in a slanted plane
     998if strcmp(ObjectData.Type,'plane_z')
     999    InterpMesh=10*InterpMesh;%TODO: temporary, to shorten computation
     1000end
    9891001
    9901002%% extrema along each axis
     
    10711083end
    10721084    icell_grid=[];% field cell index which defines the grid
    1073 if ~strcmp(ObjectData.ProjMode,'projection')
     1085if ~strcmp(ObjectData.ProjMode,'projection')&& ~strcmp(ObjectData.Type,'plane_z')% TODO:rationalize
    10741086    %% define the new coordinates in case of interpolation on a imposed grid
    1075     if ~testYMin
     1087    if ~testYMin 
    10761088        errormsg='min Y value not defined for the projection grid';return
    10771089    end
     
    11191131        AYName='coord_y';
    11201132        AXName='coord_x';
    1121         if strcmp(ObjectData.ProjMode,'projection')
     1133        if strcmp(ObjectData.ProjMode,'projection')||strcmp(ObjectData.Type,'plane_z')
    11221134            ProjData.coord_y=[FieldData.YMin FieldData.YMax];%note that if projection is done on a grid, the Min and Max along each direction must have been defined
    11231135            ProjData.coord_x=[FieldData.XMin FieldData.XMax];
     
    13991411            DimValue(DimValue==1)=[];%remove singleton dimensions
    14001412            NbDim=numel(DimValue);%update number of space dimensions
    1401             nbcolor=1; %default number of 'color' components: third matrix index without corresponding coordinate
     1413           % nbcolor=1; %default number of 'color' components: third matrix index without corresponding coordinate
    14021414            if NbDim>=3
    14031415                if NbDim>3
     
    14391451                        Coord{3}=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(3)});
    14401452                    end
    1441                     if numel(Coord{NbDim-1})==2
     1453                    if numel(Coord{NbDim-1})==2% case of coordinate defined only by the first and last values
    14421454                        DY=(Coord{NbDim-1}(2)-Coord{NbDim-1}(1))/(DimValue(1)-1);
    14431455                    end
    1444                     if numel(Coord{NbDim})==2
     1456                    if numel(Coord{NbDim})==2% case of coordinate defined only by the first and last values
    14451457                        DX=(Coord{NbDim}(2)-Coord{NbDim}(1))/(DimValue(2)-1);
    14461458                    end
     
    14531465                        end         
    14541466                    else
    1455                         YIndexMax=Coord{NbDim-1}(end)/DY;
     1467                        YIndexMax=numel(Coord{NbDim-1});
    14561468                        YIndexMin=1;
    14571469                    end
    14581470                    if testXMax
    1459                          XIndexMax=(XMax-Coord{NbDim}(1))/DY+1;% matrix index corresponding to the max y value for the new field
     1471                         XIndexMax=(XMax-Coord{NbDim}(1))/DX+1;% matrix index corresponding to the max y value for the new field
    14601472                        if testYMin%test_direct(indY)
    14611473                            XIndexMin=(XMin-Coord{NbDim}(1))/DX+1;% matrix index corresponding to the min x value for the new field
     
    14641476                        end         
    14651477                    else
    1466                         XIndexMax=Coord{NbDim}(end)/DX;
     1478                        XIndexMax=numel(Coord{NbDim});
    14671479                        XIndexMin=1;
    14681480                    end
     
    16271639                        end
    16281640                    end
    1629                 else
    1630                     errormsg='projection of structured coordinates on oblique plane not yet implemented';
    1631                     %TODO: use interp_lin3
    1632                     return
     1641                else   %projection of structured coordinates on oblique plane
     1642                    % determine the boundaries of the projected field,
     1643                    % first find the 8 summits of the initial volume in the
     1644                    % new coordinates
     1645                    Coord{1}=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(1)});%initial z coordinates
     1646                    Coord{2}=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(2)});%initial y coordinates
     1647                    Coord{3}=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(3)});%initial x coordinates
     1648                    summit=zeros(3,8);% initialize summit coordinates
     1649                    summit(1,1:4)=[Coord{3}(1) Coord{3}(end) Coord{3}(1) Coord{3}(end)];%square
     1650                    summit(2,1:4)=[Coord{2}(1) Coord{2}(1) Coord{2}(end) Coord{2}(end)];% square at z= Coord{1}(1)
     1651                    summit(1:2,5:8)=summit(1:2,1:4);
     1652                    summit(3,:)=[Coord{1}(1)*ones(1,4) Coord{1}(end)*ones(1,4)];
     1653                    Mrot=rodrigues(PlaneAngle);
     1654                    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
     1656                    for isummit=1:8% TODO: introduce a function for rotation of n points (to use also for scattered data)
     1657                        newsummit(:,isummit)=Mrot*(summit(:,isummit)-(ObjectData.Coord));
     1658                    end
     1659                    coord_x_proj=min(newsummit(1,:)):InterpMesh: max(newsummit(1,:));
     1660                    coord_y_proj=min(newsummit(2,:)):InterpMesh: max(newsummit(2,:));
     1661                    coord_z_proj=-width:InterpMesh:width;
     1662                    Mrot_inverse=rodrigues(-PlaneAngle);
     1663                    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)
     1669                        Grid_x=shiftdim(Grid_x,-1);
     1670                        Grid_y=shiftdim(Grid_y,-1);
     1671                        Grid_z=shiftdim(Grid_z,-1);
     1672                    end
     1673                    XI=Origin(1)+ix(1)*Grid_x+iy(1)*Grid_y+iz(1)*Grid_z;
     1674                    YI=Origin(2)+ix(2)*Grid_x+iy(2)*Grid_y+iz(2)*Grid_z;
     1675                    ZI=Origin(3)+ix(3)*Grid_x+iy(3)*Grid_y+iz(3)*Grid_z;
     1676                    [X,Y,Z]=meshgrid(Coord{3},Coord{2},Coord{1});
     1677                    X=permute(X,[3 1 2]);
     1678                    Y=permute(Y,[3 1 2]);
     1679                    Z=permute(Z,[3 1 2]);
     1680                    for ivar=VarIndex
     1681                            VarName=FieldData.ListVarName{ivar};
     1682                            ListVarName=[ListVarName VarName];
     1683                            VarAttribute{length(ListVarName)}=FieldData.VarAttribute{ivar}; %reproduce the variable attributes
     1684                            ProjData.(VarName)=interp3(X,Y,Z,double(FieldData.(VarName)),XI,YI,ZI);
     1685                    end
    16331686                end
    16341687            end
     
    22642317                    end
    22652318                else
     2319                    RotMatrix=rodrigues(om);
     2320                   
    22662321                    errormsg='projection of structured coordinates on oblique plane not yet implemented';
    22672322                    %TODO: use interp3
Note: See TracChangeset for help on using the changeset viewer.