Changeset 227 for trunk/src/proj_field.m
- Timestamp:
- Mar 31, 2011, 1:42:51 PM (14 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/src/proj_field.m
r215 r227 946 946 cos_om=1; 947 947 sin_om=0; 948 test90x=0;%=1 for 90 degree rotation alround x axis 949 test90y=0;%=1 for 90 degree rotation alround y axis 948 950 if isfield(ObjectData,'Angle')&& isequal(size(ObjectData.Angle),[1 3])&& ~isequal(ObjectData.Angle,[0 0 0]) 951 test90y=isequal(ObjectData.Angle,[0 90 0]); 949 952 PlaneAngle=(pi/180)*ObjectData.Angle; 950 953 om=norm(PlaneAngle);%norm of rotation angle in radians … … 958 961 norm_plane(3)=OmAxis(3)*coeff+cos_om; 959 962 end 960 testangle=~isequal(PlaneAngle,[0 0 0]); 963 testangle=~isequal(PlaneAngle,[0 0 0]);% && ~test90y && ~test90x;%=1 for slanted plane 964 961 965 % Phi=0;%default 962 966 % Theta=0; … … 1071 1075 DimCell={DimCell};%name of dimensions 1072 1076 end 1073 1074 %% case of input fields with unstructured coordinates1077 1078 %% case of input fields with unstructured coordinates 1075 1079 if testX 1076 1080 XName=FieldData.ListVarName{ivar_X}; … … 1082 1086 eval(['coord_z=FieldData.' ZName ';']) 1083 1087 end 1084 1088 1085 1089 % translate initial coordinates 1086 1090 coord_x=coord_x-ObjectData.Coord(1,1); … … 1093 1097 if length(ivar_Z)==1 && width > 0 1094 1098 %components of the unitiy vector normal to the projection plane 1095 fieldZ=norm_plane(1)*coord_x + norm_plane(2)*coord_y+ norm_plane(3)*coord_z;% distance to the plane 1099 fieldZ=norm_plane(1)*coord_x + norm_plane(2)*coord_y+ norm_plane(3)*coord_z;% distance to the plane 1096 1100 indcut=find(abs(fieldZ) <= width); 1097 1101 size(indcut) 1098 1102 for ivar=VarIndex 1099 1103 VarName=FieldData.ListVarName{ivar}; 1100 eval(['FieldData.' VarName '=FieldData.' VarName '(indcut);']) 1101 % A VOIR : CAS DE VAR STRUCTUREE MAIS PAS GRILLE REGULIERE : INTERPOLER SUR GRILLE REGULIERE1104 eval(['FieldData.' VarName '=FieldData.' VarName '(indcut);']) 1105 % A VOIR : CAS DE VAR STRUCTUREE MAIS PAS GRILLE REGULIERE : INTERPOLER SUR GRILLE REGULIERE 1102 1106 end 1103 1107 coord_x=coord_x(indcut); … … 1105 1109 coord_z=coord_z(indcut); 1106 1110 end 1107 1108 %rotate coordinates if needed1109 if testangle 1110 % coord_X=coord_x;1111 % coord_Y=coord_y;1112 % if ~isequal(Theta,0)1113 % coord_Y=coord_Y *cos(Theta);1114 % end1115 % else1111 1112 %rotate coordinates if needed: 1113 if testangle && ~test90y && ~test90x;%=1 for slanted plane 1114 % coord_X=coord_x; 1115 % coord_Y=coord_y; 1116 % if ~isequal(Theta,0) 1117 % coord_Y=coord_Y *cos(Theta); 1118 % end 1119 % else 1116 1120 coord_X=(coord_x *cos(Phi) + coord_y* sin(Phi)); 1117 1121 coord_Y=(-coord_x *sin(Phi) + coord_y *cos(Phi))*cos(Theta); 1118 % end1119 % if ~isempty(ivar_Z)1122 % end 1123 % if ~isempty(ivar_Z) 1120 1124 coord_Y=coord_Y+coord_z *sin(Theta); 1121 % end1122 % if testangle1123 1124 1125 % end 1126 % if testangle 1127 coord_X=(coord_X *cos(Psi) - coord_Y* sin(Psi));%A VERIFIER 1128 coord_Y=(coord_X *sin(Psi) + coord_Y* cos(Psi)); 1125 1129 else 1126 1130 coord_X=coord_x; … … 1151 1155 for ivar=VarIndex 1152 1156 VarName=FieldData.ListVarName{ivar}; 1153 eval(['FieldData.' VarName '=FieldData.' VarName '(indcut);']) 1157 eval(['FieldData.' VarName '=FieldData.' VarName '(indcut);']) 1154 1158 end 1155 1159 coord_X=coord_X(indcut); … … 1164 1168 %ProjData.ListDimName=[ProjData.ListDimName FieldData.VarDimName(VarIndex(1))];%add the point index to the list of dimensions 1165 1169 %ProjData.DimValue=[ProjData. 1166 1167 1170 %length(coord_X)]; 1171 1168 1172 for ivar=VarIndex %transfer variables to the projection plane 1169 1173 VarName=FieldData.ListVarName{ivar}; … … 1175 1179 eval(['ProjData.' VarName '=FieldData.' VarName ';']) 1176 1180 end 1177 if isempty(ivar_Z) || ivar~=ivar_Z 1181 if isempty(ivar_Z) || ivar~=ivar_Z 1178 1182 ProjData.ListVarName=[ProjData.ListVarName VarName]; 1179 1183 ProjData.VarDimName=[ProjData.VarDimName DimCell]; … … 1183 1187 end 1184 1188 end 1185 end 1189 end 1186 1190 elseif isequal(ObjectData.ProjMode,'interp')||isequal(ObjectData.ProjMode,'filter')%interpolate data on a regular grid 1187 1191 coord_x_proj=XMin:DX:XMax; … … 1189 1193 DimCell={'coord_y','coord_x'}; 1190 1194 ProjData.ListVarName={'coord_y','coord_x'}; 1191 ProjData.VarDimName={'coord_y','coord_x'}; 1192 nbcoord=2; 1195 ProjData.VarDimName={'coord_y','coord_x'}; 1196 nbcoord=2; 1193 1197 ProjData.coord_y=[YMin YMax]; 1194 1198 ProjData.coord_x=[XMin XMax]; … … 1211 1215 for ivar=VarIndex 1212 1216 VarName=FieldData.ListVarName{ivar}; 1213 if ~( ivar==ivar_X || ivar==ivar_Y || ivar==ivar_Z || ivar==ivar_F || ivar==ivar_FF || test_anc(ivar)==1) 1217 if ~( ivar==ivar_X || ivar==ivar_Y || ivar==ivar_Z || ivar==ivar_F || ivar==ivar_FF || test_anc(ivar)==1) 1214 1218 ivar_new=ivar_new+1; 1215 1219 ProjData.ListVarName=[ProjData.ListVarName {VarName}]; … … 1245 1249 ProjData.FF=reshape(FF,length(coord_y_proj),length(coord_x_proj)); 1246 1250 ProjData.ListVarName=[ProjData.ListVarName {'FF'}]; 1247 ProjData.VarDimName=[ProjData.VarDimName {DimCell}];1251 ProjData.VarDimName=[ProjData.VarDimName {DimCell}]; 1248 1252 ProjData.VarAttribute{ivar_new+1+nbcoord}.Role='errorflag'; 1249 1253 end 1250 1254 end 1251 1255 1252 %% case of input fields defined on a structured grid 1256 %% case of input fields defined on a structured grid 1253 1257 else 1254 VarName=FieldData.ListVarName{VarIndex(1)};%get the first variable of the cell to get the input matrix dimensions 1258 'TESTproj' 1259 VarName=FieldData.ListVarName{VarIndex(1)}%get the first variable of the cell to get the input matrix dimensions 1255 1260 eval(['DimValue=size(FieldData.' VarName ');'])%input matrix dimensions 1256 DimValue(DimValue==1)=[];%remove singleton dimensions 1261 DimValue(DimValue==1)=[];%remove singleton dimensions 1257 1262 NbDim=numel(DimValue);%update number of space dimensions 1258 1263 nbcolor=1; %default number of 'color' components: third matrix index without corresponding coordinate … … 1270 1275 end 1271 1276 end 1272 AYName=FieldData.ListVarName{VarType.coord(NbDim-1)};%name of input x coordinate (name preserved on projection) 1273 AXName=FieldData.ListVarName{VarType.coord(NbDim)};%name of input y coordinate (name preserved on projection) 1277 testangle 1278 if testangle% TODO modify name also in case of origin shift in x or y 1279 AYName='Y'; 1280 AXName='X'; 1281 count=0; 1282 %modify coordinate names if they are already used 1283 while ~(isempty(find(strcmp('AXName',ProjData.ListVarName),1)) && isempty(find(strcmp('AYName',ProjData.ListVarName),1))) 1284 count=count+1; 1285 AYName=[AYName '_' num2str(count)]; 1286 AXName=[AXName '_' num2str(count)]; 1287 end 1288 else 1289 AYName=FieldData.ListVarName{VarType.coord(NbDim-1)}%name of input x coordinate (name preserved on projection) 1290 AXName=FieldData.ListVarName{VarType.coord(NbDim)}%name of input y coordinate (name preserved on projection) 1291 end 1274 1292 eval(['AX=FieldData.' AXName ';']) 1275 1293 eval(['AY=FieldData.' AYName ';']) 1276 1294 ListDimName=FieldData.VarDimName{VarIndex(1)}; 1277 ProjData.ListVarName=[ProjData.ListVarName {AYName} {AXName}]; %TODO: check if it already exists in Projdata (several cells) 1278 ProjData.VarDimName=[ProjData.VarDimName {AYName} {AXName}]; 1279 1280 % for idim=1:length(ListDimName) 1281 % DimName=ListDimName{idim}; 1282 % if strcmp(DimName,'rgb')||strcmp(DimName,'nb_coord')||strcmp(DimName,'nb_coord_i') 1283 % nbcolor=DimValue(idim); 1284 % DimValue(idim)=[]; 1285 % end 1286 % if isequal(DimName,'nb_coord_j')% NOTE: CASE OF TENSOR NOT TREATED 1287 % DimValue(idim)=[]; 1288 % end 1289 % end 1295 ProjData.ListVarName=[ProjData.ListVarName {AYName} {AXName}]; %TODO: check if it already exists in Projdata (several cells) 1296 ProjData.VarDimName=[ProjData.VarDimName {AYName} {AXName}]; 1290 1297 Coord_z=[]; 1291 1298 Coord_y=[]; 1292 Coord_x=[]; 1293 1299 Coord_x=[]; 1300 1294 1301 for idim=1:NbDim %loop on space dimensions 1295 1302 test_interp(idim)=0;%test for coordiate interpolation (non regular grid), =0 by default … … 1298 1305 eval(['Coord{idim}=FieldData.' FieldData.ListVarName{ivar} ';']) ;% coord values for the input field 1299 1306 if numel(Coord{idim})==2 %input array defined on a regular grid 1300 DCoord_min(idim)=(Coord{idim}(2)-Coord{idim}(1))/DimValue(idim);1307 DCoord_min(idim)=(Coord{idim}(2)-Coord{idim}(1))/DimValue(idim); 1301 1308 else 1302 1309 DCoord=diff(Coord{idim});%array of coordinate derivatives for the input field 1303 1310 DCoord_min(idim)=min(DCoord); 1304 1311 DCoord_max=max(DCoord); 1305 % test_direct(idim)=DCoord_max>0;% =1 for increasing values, 0 otherwise1306 if abs(DCoord_max-DCoord_min(idim))>abs(DCoord_max/1000) 1312 % test_direct(idim)=DCoord_max>0;% =1 for increasing values, 0 otherwise 1313 if abs(DCoord_max-DCoord_min(idim))>abs(DCoord_max/1000) 1307 1314 msgbox_uvmat('ERROR',['non monotonic dimension variable # ' num2str(idim) ' in proj_field.m']) 1308 1309 end 1315 return 1316 end 1310 1317 test_interp(idim)=(DCoord_max-DCoord_min(idim))> 0.0001*abs(DCoord_max);% test grid regularity 1311 1318 end … … 1321 1328 DY=abs(DCoord_min(NbDim-1)); 1322 1329 end 1323 npY=1+round(abs(Coord{NbDim-1}(end)-Coord{NbDim-1}(1))/DY);%nbre of points after interpol 1330 npY=1+round(abs(Coord{NbDim-1}(end)-Coord{NbDim-1}(1))/DY);%nbre of points after interpol 1324 1331 if DX==0 1325 1332 DX=abs(DCoord_min(NbDim)); 1326 1333 end 1327 npX=1+round(abs(Coord{NbDim}(end)-Coord{NbDim}(1))/DX);%nbre of points after interpol 1334 npX=1+round(abs(Coord{NbDim}(end)-Coord{NbDim}(1))/DX);%nbre of points after interpol 1328 1335 for idim=1:NbDim 1329 1336 if test_interp(idim) 1330 1337 DimValue(idim)=1+round(abs(Coord{idim}(end)-Coord{idim}(1))/abs(DCoord_min(idim)));%nbre of points after possible interpolation on a regular gri 1331 1338 end 1332 end 1339 end 1333 1340 Coord_y=linspace(Coord{NbDim-1}(1),Coord{NbDim-1}(end),npY); 1334 1341 test_direct_y=test_direct(NbDim-1); … … 1336 1343 test_direct_x=test_direct(NbDim); 1337 1344 DAX=DCoord_min(NbDim); 1338 DAY=DCoord_min(NbDim-1); 1345 DAY=DCoord_min(NbDim-1); 1339 1346 minAX=min(Coord_x); 1340 1347 maxAX=max(Coord_x); … … 1374 1381 end 1375 1382 npX=floor((XMax-XMin)/DX+1); 1376 npY=floor((YMax-YMin)/DY+1); 1383 npY=floor((YMax-YMin)/DY+1); 1377 1384 if test_direct_y 1378 1385 coord_y_proj=linspace(YMin,YMax,npY);%abscissa of the new pixels along the line … … 1384 1391 else 1385 1392 coord_x_proj=linspace(XMax,XMin,npX);%abscissa of the new pixels along the line 1386 end 1387 1393 end 1388 1394 % case with no rotation and interpolation 1389 1395 if isequal(ProjMode,'projection') && ~testangle 1390 1396 if ~testXMin && ~testXMax && ~testYMin && ~testYMax && NbDim==2 1391 ProjData=FieldData; 1397 ProjData=FieldData; 1392 1398 else 1393 1399 indY=NbDim-1; … … 1402 1408 Ybound(2)=Coord{indY}(1)-DYinit*(max_indy-1); 1403 1409 Ybound(1)=Coord{indY}(1)-DYinit*(min_indy-1); 1404 end 1410 end 1405 1411 if test_direct(NbDim)==1 1406 1412 min_indx=ceil((XMin-Coord{NbDim}(1))/DXinit)+1; … … 1412 1418 max_indx=floor((Coord{NbDim}(1)-XMin)/DXinit)+1; 1413 1419 Xbound(2)=Coord{NbDim}(1)+DXinit*(max_indx-1); 1420 1421 1414 1422 Xbound(1)=Coord{NbDim}(1)+DXinit*(min_indx-1); 1415 end1416 if NbDim==31417 DimCell(1)=[]; %suppress z variable1418 DimValue(1)=[];1419 %structured coordinates1420 if test_direct(1)1421 iz=ceil((ObjectData.Coord(1,3)-Coord{1}(1))/DZ)+1;1422 else1423 iz=ceil((Coord{1}(1)-ObjectData.Coord(1,3))/DZ)+1;1424 end1425 1423 end 1426 1424 min_indy=max(min_indy,1);% deals with margin (bound lower than the first index) 1427 1425 min_indx=max(min_indx,1); 1428 max_indy=min(max_indy,DimValue(1)); 1429 max_indx=min(max_indx,DimValue(2)); 1430 for ivar=VarIndex% loop on non coordinate variables 1431 VarName=FieldData.ListVarName{ivar}; 1432 ProjData.ListVarName=[ProjData.ListVarName VarName]; 1433 ProjData.VarDimName=[ProjData.VarDimName {DimCell}]; 1434 if isfield(FieldData,'VarAttribute') && length(FieldData.VarAttribute)>=ivar 1435 ProjData.VarAttribute{length(ProjData.ListVarName)}=FieldData.VarAttribute{ivar}; 1436 end 1426 1427 if test90y 1428 'TEST3D' 1429 ind_new=[3 2 1]; 1430 DimCell=DimCell(ind_new); 1431 DimValue=DimValue(ind_new); 1432 DimCell(1)=[]; %suppress x variable 1433 DimValue(1)=[]; 1434 iz=ceil((ObjectData.Coord(1,1)-Coord{3}(1))/DX)+1; 1435 for ivar=VarIndex 1436 VarName=FieldData.ListVarName{ivar}; 1437 ProjData.ListVarName=[ProjData.ListVarName VarName]; 1438 ProjData.VarDimName=[ProjData.VarDimName {DimCell}]; 1439 ProjData.VarAttribute{length(ProjData.ListVarName)}=FieldData.VarAttribute{ivar}; %reproduce the variable attributes 1440 eval(['size(FieldData.' VarName ')']) 1441 eval(['ProjData.' VarName '=permute(FieldData.' VarName ',ind_new)'])% permute x and z indices for 90 degree rotation 1442 eval(['size(ProjData.' VarName ')']) 1443 eval(['ProjData.' VarName '=squeeze(FieldData.' VarName '(iz,:,:));'])% select the z index iz 1444 end 1445 eval(['ProjData.' AYName '=[Ybound(1) Ybound(2)];']) %record the new (projected ) y coordinates 1446 eval(['ProjData.' AXName '=[Coord{1}(end),Coord{1}(1)];']) %record the new (projected ) x coordinates 1447 else 1437 1448 if NbDim==3 1438 eval(['ProjData.' VarName '=squeeze(FieldData.' VarName '(iz,min_indy:max_indy,min_indx:max_indx));']); 1439 else 1440 eval(['ProjData.' VarName '=FieldData.' VarName '(min_indy:max_indy,min_indx:max_indx,:);']); 1441 end 1442 end 1443 eval(['ProjData.' AYName '=[Ybound(1) Ybound(2)];']) %record the new (projected ) y coordinates 1444 eval(['ProjData.' AXName '=[Xbound(1) Xbound(2)];']) %record the new (projected ) x coordinates 1449 DimCell(1)=[]; %suppress z variable 1450 DimValue(1)=[]; 1451 if test_direct(1) 1452 iz=ceil((ObjectData.Coord(1,3)-Coord{1}(1))/DZ)+1; 1453 else 1454 iz=ceil((Coord{1}(1)-ObjectData.Coord(1,3))/DZ)+1; 1455 end 1456 end 1457 max_indy=min(max_indy,DimValue(1));%introduce bounds in y and x indices 1458 max_indx=min(max_indx,DimValue(2)); 1459 for ivar=VarIndex% loop on non coordinate variables 1460 VarName=FieldData.ListVarName{ivar}; 1461 ProjData.ListVarName=[ProjData.ListVarName VarName]; 1462 ProjData.VarDimName=[ProjData.VarDimName {DimCell}]; 1463 if isfield(FieldData,'VarAttribute') && length(FieldData.VarAttribute)>=ivar 1464 ProjData.VarAttribute{length(ProjData.ListVarName)}=FieldData.VarAttribute{ivar}; 1465 end 1466 if NbDim==3 1467 eval(['ProjData.' VarName '=squeeze(FieldData.' VarName '(iz,min_indy:max_indy,min_indx:max_indx));']); 1468 else 1469 eval(['ProjData.' VarName '=FieldData.' VarName '(min_indy:max_indy,min_indx:max_indx,:);']); 1470 end 1471 end 1472 eval(['ProjData.' AYName '=[Ybound(1) Ybound(2)];']) %record the new (projected ) y coordinates 1473 eval(['ProjData.' AXName '=[Xbound(1) Xbound(2)];']) %record the new (projected ) x coordinates 1474 end 1445 1475 end 1446 1476 else % case with rotation and/or interpolation … … 1453 1483 XIMA=reshape(round(XIMA),1,npX*npY);%indices reorganized in 'line' 1454 1484 YIMA=reshape(round(YIMA),1,npX*npY); 1455 flagin=XIMA>=1 & XIMA<=DimValue(2) & YIMA >=1 & YIMA<=DimValue(1);%flagin=1 inside the original image 1485 flagin=XIMA>=1 & XIMA<=DimValue(2) & YIMA >=1 & YIMA<=DimValue(1);%flagin=1 inside the original image 1456 1486 if isequal(ObjectData.ProjMode,'filter') 1457 1487 npx_filter=ceil(abs(DX/DAX)); … … 1466 1496 for ivar=VarIndex 1467 1497 VarName=FieldData.ListVarName{ivar}; 1468 if test_interp(1) || test_interp(2)%interpolate on a regular grid 1469 1498 if test_interp(1) || test_interp(2)%interpolate on a regular grid 1499 eval(['ProjData.' VarName '=interp2(Coord{2},Coord{1},FieldData.' VarName ',Coord_x,Coord_y'');']) %TO TEST 1470 1500 end 1471 1501 %filter the field (image) if option 'filter' is used 1472 if test_filter 1473 1474 1475 1476 1477 1478 end 1479 eval(['vec_A=reshape(FieldData.' VarName ',[],nbcolor);'])%put the original image in line 1502 if test_filter 1503 Aclass=class(FieldData.A); 1504 eval(['ProjData.' VarName '=filter2(Mfilter,FieldData.' VarName ',''valid'');']) 1505 if ~isequal(Aclass,'double') 1506 eval(['ProjData.' VarName '=' Aclass '(FieldData.' VarName ');'])%revert to integer values 1507 end 1508 end 1509 eval(['vec_A=reshape(FieldData.' VarName ',[],nbcolor);'])%put the original image in line 1480 1510 %ind_in=find(flagin); 1481 1511 ind_out=find(~flagin); 1482 1512 ICOMB=(XIMA-1)*DimValue(1)+YIMA; 1483 1513 ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A 1484 vec_B(flagin,1:nbcolor)=vec_A(ICOMB,:); 1514 vec_B(flagin,1:nbcolor)=vec_A(ICOMB,:); 1485 1515 for icolor=1:nbcolor 1486 1516 vec_B(ind_out,icolor)=zeros(size(ind_out)); … … 1490 1520 if isfield(FieldData,'VarAttribute')&&length(FieldData.VarAttribute)>=ivar 1491 1521 ProjData.VarAttribute{length(ProjData.ListVarName)+nbcoord}=FieldData.VarAttribute{ivar}; 1492 end 1522 end 1493 1523 eval(['ProjData.' VarName '=reshape(vec_B,npY,npX,nbcolor);']); 1494 1524 end 1495 ProjData.FF=reshape(~flagin,npY,npX);%false flag A FAIRE: tenir compte d'un flga antérieur 1525 ProjData.FF=reshape(~flagin,npY,npX);%false flag A FAIRE: tenir compte d'un flga antérieur 1496 1526 ProjData.ListVarName=[ProjData.ListVarName 'FF']; 1497 1527 ProjData.VarDimName=[ProjData.VarDimName {DimCell}]; 1498 1528 ProjData.VarAttribute{length(ProjData.ListVarName)}.Role='errorflag'; 1499 else %3D case 1500 if ~testangle 1501 % unstructured z coordinate 1502 test_sup=(Coord{1}>=ObjectData.Coord(1,3)); 1503 iz_sup=find(test_sup); 1504 iz=iz_sup(1); 1505 if iz>=1 & iz<=npz 1506 %ProjData.ListDimName=[ProjData.ListDimName ListDimName(2:end)]; 1507 %ProjData.DimValue=[ProjData.DimValue npY npX]; 1508 for ivar=VarIndex 1509 VarName=FieldData.ListVarName{ivar}; 1510 ProjData.ListVarName=[ProjData.ListVarName VarName]; 1511 ProjData.VarAttribute{length(ProjData.ListVarName)}=FieldData.VarAttribute{ivar}; %reproduce the variable attributes 1512 eval(['ProjData.' VarName '=squeeze(FieldData.' VarName '(iz,:,:));'])% select the z index iz 1513 %TODO : do a vertical average for a thick plane 1514 if test_interp(2) || test_interp(3) 1515 eval(['ProjData.' VarName '=interp2(Coord{3},Coord{2},ProjData.' VarName ',Coord_x,Coord_y'');']) 1516 end 1529 elseif ~testangle 1530 % unstructured z coordinate 1531 test_sup=(Coord{1}>=ObjectData.Coord(1,3)); 1532 iz_sup=find(test_sup); 1533 iz=iz_sup(1); 1534 if iz>=1 & iz<=npz 1535 %ProjData.ListDimName=[ProjData.ListDimName ListDimName(2:end)]; 1536 %ProjData.DimValue=[ProjData.DimValue npY npX]; 1537 for ivar=VarIndex 1538 VarName=FieldData.ListVarName{ivar}; 1539 ProjData.ListVarName=[ProjData.ListVarName VarName]; 1540 ProjData.VarAttribute{length(ProjData.ListVarName)}=FieldData.VarAttribute{ivar}; %reproduce the variable attributes 1541 eval(['ProjData.' VarName '=squeeze(FieldData.' VarName '(iz,:,:));'])% select the z index iz 1542 %TODO : do a vertical average for a thick plane 1543 if test_interp(2) || test_interp(3) 1544 eval(['ProjData.' VarName '=interp2(Coord{3},Coord{2},ProjData.' VarName ',Coord_x,Coord_y'');']) 1517 1545 end 1518 1546 end 1519 e lse1520 errormsg='projection of structured coordinates on oblique plane not yet implemented';1521 %TODO: use interp31522 return1523 end1524 end 1525 end 1526 end 1527 1547 end 1548 else 1549 errormsg='projection of structured coordinates on oblique plane not yet implemented'; 1550 %TODO: use interp3 1551 return 1552 end 1553 end 1554 end 1555 1528 1556 %% projection of velocity components in the rotated coordinates 1529 1557 if testangle && length(ivar_U)==1 … … 1533 1561 end 1534 1562 UName=FieldData.ListVarName{ivar_U}; 1535 VName=FieldData.ListVarName{ivar_V}; 1563 VName=FieldData.ListVarName{ivar_V}; 1536 1564 eval(['ProjData.' UName '=cos(PlaneAngle(3))*ProjData.' UName '+ sin(PlaneAngle(3))*ProjData.' VName ';']) 1537 1565 eval(['ProjData.' VName '=cos(Theta)*(-sin(PlaneAngle(3))*ProjData.' UName '+ cos(PlaneAngle(3))*ProjData.' VName ');']) 1538 1566 if ~isempty(ivar_W) 1539 1567 WName=FieldData.ListVarName{ivar_W}; 1540 eval(['ProjData.' VName '=ProjData.' VName '+ ProjData.' WName '*sin(Theta);'])% 1568 eval(['ProjData.' VName '=ProjData.' VName '+ ProjData.' WName '*sin(Theta);'])% 1541 1569 eval(['ProjData.' WName '=NormVec_X*ProjData.' UName '+ NormVec_Y*ProjData.' VName '+ NormVec_Z* ProjData.' WName ';']); 1542 1570 end … … 1857 1885 VarName=FieldData.ListVarName{VarIndex(1)};%get the first variable of the cell to get the input matrix dimensions 1858 1886 eval(['DimValue=size(FieldData.' VarName ');'])%input matrix dimensions 1859 DimValue( find(DimValue==1))=[];%remove singleton dimensions1887 DimValue(DimValue==1)=[];%remove singleton dimensions 1860 1888 NbDim=numel(DimValue);%update number of space dimensions 1861 1889 nbcolor=1; %default number of 'color' components: third matrix index without corresponding coordinate … … 1865 1893 return 1866 1894 else 1867 VarType.coord1868 1895 if numel(find(VarType.coord))==2% the third matrix dimension does not correspond to a space coordinate 1869 1896 nbcolor=DimValue(3);
Note: See TracChangeset
for help on using the changeset viewer.