Changeset 954 for trunk/src/proj_field.m
 Timestamp:
 Jun 22, 2016, 1:36:50 PM (8 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/src/proj_field.m
r937 r954 90 90 return 91 91 end 92 ListProjMode={'projection','interp_lin','interp_tps','inside','outside'};%list of effective projection modes93 if isempty(find(strcmp(ObjectData.ProjMode,ListProjMode), 1))% no projection in case92 % check list of effective projection modes 93 if ~ismember(ObjectData.ProjMode,{'projection','interp_lin','interp_tps','inside','outside'}) 94 94 return 95 95 end … … 117 117 [ProjData,errormsg] = proj_line(FieldData,ObjectData); 118 118 end 119 case 'plane'119 case {'plane','plane_z'} 120 120 [ProjData,errormsg] = proj_plane(FieldData,ObjectData); 121 121 case 'volume' … … 955 955 test90x=0;%=1 for 90 degree rotation alround x axis 956 956 test90y=0;%=1 for 90 degree rotation alround y axis 957 if 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; 964 end 957 965 if isfield(ObjectData,'Angle')&& isequal(size(ObjectData.Angle),[1 3])&& ~isequal(ObjectData.Angle,[0 0 0]) 958 966 test90y=isequal(ObjectData.Angle,[0 90 0]); … … 987 995 return 988 996 end 997 InterpMesh=min(DX,DY);%mesh used for interpolation in a slanted plane 998 if strcmp(ObjectData.Type,'plane_z') 999 InterpMesh=10*InterpMesh;%TODO: temporary, to shorten computation 1000 end 989 1001 990 1002 %% extrema along each axis … … 1071 1083 end 1072 1084 icell_grid=[];% field cell index which defines the grid 1073 if ~strcmp(ObjectData.ProjMode,'projection') 1085 if ~strcmp(ObjectData.ProjMode,'projection')&& ~strcmp(ObjectData.Type,'plane_z')% TODO:rationalize 1074 1086 %% define the new coordinates in case of interpolation on a imposed grid 1075 if ~testYMin 1087 if ~testYMin 1076 1088 errormsg='min Y value not defined for the projection grid';return 1077 1089 end … … 1119 1131 AYName='coord_y'; 1120 1132 AXName='coord_x'; 1121 if strcmp(ObjectData.ProjMode,'projection') 1133 if strcmp(ObjectData.ProjMode,'projection')strcmp(ObjectData.Type,'plane_z') 1122 1134 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 1123 1135 ProjData.coord_x=[FieldData.XMin FieldData.XMax]; … … 1399 1411 DimValue(DimValue==1)=[];%remove singleton dimensions 1400 1412 NbDim=numel(DimValue);%update number of space dimensions 1401 nbcolor=1; %default number of 'color' components: third matrix index without corresponding coordinate1413 % nbcolor=1; %default number of 'color' components: third matrix index without corresponding coordinate 1402 1414 if NbDim>=3 1403 1415 if NbDim>3 … … 1439 1451 Coord{3}=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(3)}); 1440 1452 end 1441 if numel(Coord{NbDim1})==2 1453 if numel(Coord{NbDim1})==2% case of coordinate defined only by the first and last values 1442 1454 DY=(Coord{NbDim1}(2)Coord{NbDim1}(1))/(DimValue(1)1); 1443 1455 end 1444 if numel(Coord{NbDim})==2 1456 if numel(Coord{NbDim})==2% case of coordinate defined only by the first and last values 1445 1457 DX=(Coord{NbDim}(2)Coord{NbDim}(1))/(DimValue(2)1); 1446 1458 end … … 1453 1465 end 1454 1466 else 1455 YIndexMax= Coord{NbDim1}(end)/DY;1467 YIndexMax=numel(Coord{NbDim1}); 1456 1468 YIndexMin=1; 1457 1469 end 1458 1470 if testXMax 1459 XIndexMax=(XMaxCoord{NbDim}(1))/D Y+1;% matrix index corresponding to the max y value for the new field1471 XIndexMax=(XMaxCoord{NbDim}(1))/DX+1;% matrix index corresponding to the max y value for the new field 1460 1472 if testYMin%test_direct(indY) 1461 1473 XIndexMin=(XMinCoord{NbDim}(1))/DX+1;% matrix index corresponding to the min x value for the new field … … 1464 1476 end 1465 1477 else 1466 XIndexMax= Coord{NbDim}(end)/DX;1478 XIndexMax=numel(Coord{NbDim}); 1467 1479 XIndexMin=1; 1468 1480 end … … 1627 1639 end 1628 1640 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 1633 1686 end 1634 1687 end … … 2264 2317 end 2265 2318 else 2319 RotMatrix=rodrigues(om); 2320 2266 2321 errormsg='projection of structured coordinates on oblique plane not yet implemented'; 2267 2322 %TODO: use interp3
Note: See TracChangeset
for help on using the changeset viewer.