- Timestamp:
- Jun 26, 2013, 11:24:42 PM (11 years ago)
- Location:
- trunk/src
- Files:
-
- 5 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/src/interp2_uvmat.m
r646 r653 1 %'interp2_uvmat': linearly interpolate an image or scalar defined on a regular grid 1 %'interp2_uvmat': linearly interpolate an image or scalar defined on a 2 %regular grid USE Matlab interp2 instead 2 3 %-------------------------------------------------------------------- 3 4 %OUTPUT: -
trunk/src/proj_field.m
r651 r653 873 873 %----------------------------------------------------------------- 874 874 875 %% rotation angles 876 PlaneAngle=[0 0 0]; 875 %% rotation angles 876 PlaneAngle=[0 0 0]; 877 877 norm_plane=[0 0 1]; 878 878 cos_om=1; … … 893 893 norm_plane(3)=OmAxis(3)*coeff+cos_om; 894 894 end 895 testangle=~isequal(PlaneAngle,[0 0 0]);% && ~test90y && ~test90x;%=1 for slanted plane 895 testangle=~isequal(PlaneAngle,[0 0 0]);% && ~test90y && ~test90x;%=1 for slanted plane 896 896 897 897 %% mesh sizes DX and DY … … 899 899 DY=[];%default 900 900 if isfield(ObjectData,'DX') && ~isempty(ObjectData.DX) 901 DX=abs(ObjectData.DX);%mesh of interpolation points901 DX=abs(ObjectData.DX);%mesh of interpolation points 902 902 elseif isfield(FieldData,'CoordMesh') 903 903 DX=FieldData.CoordMesh; 904 904 end 905 905 if isfield(ObjectData,'DY') && ~isempty(ObjectData.DY) 906 DY=abs(ObjectData.DY);%mesh of interpolation points906 DY=abs(ObjectData.DY);%mesh of interpolation points 907 907 elseif isfield(FieldData,'CoordMesh') 908 908 DY=FieldData.CoordMesh; 909 909 end 910 910 if ~strcmp(ObjectData.ProjMode,'projection') && (isempty(DX)||isempty(DY)) 911 912 911 errormsg='DX or DY not defined'; 912 return 913 913 end 914 914 … … 918 918 testYMin=0; 919 919 testYMax=0; 920 921 920 if isfield(ObjectData,'RangeX') 922 923 924 925 926 else927 XMin=FieldData.XMin;%default928 XMax=FieldData.XMax;%default921 XMin=min(ObjectData.RangeX); 922 XMax=max(ObjectData.RangeX); 923 testXMin=XMax>XMin; 924 testXMax=1;% range restriction along X 925 % else 926 % XMin=FieldData.XMin;%default 927 % XMax=FieldData.XMax;%default 929 928 end 930 929 if isfield(ObjectData,'RangeY') 931 932 933 934 935 else936 YMin=FieldData.YMin;%default937 YMax=FieldData.YMax;%default930 YMin=min(ObjectData.RangeY); 931 YMax=max(ObjectData.RangeY); 932 testYMin=YMax>YMin; 933 testYMax=1; 934 % else 935 % YMin=FieldData.YMin;%default 936 % YMax=FieldData.YMax;%default 938 937 end 939 938 width=0;%default width of the projection band 940 939 if isfield(ObjectData,'RangeZ') 941 940 width=max(ObjectData.RangeZ); 942 941 end 943 942 … … 950 949 %% reproduce initial plane position and angle 951 950 if isfield(FieldData,'PlaneCoord')&&length(FieldData.PlaneCoord)==3&& isfield(ProjData,'ProjObjectCoord') 952 if length(ProjData.ProjObjectCoord)==3% if the projection plane has a z coordinate953 if ~isequal(ProjData.PlaneCoord(3),ProjData.ProjObjectCoord) %check the consistency with the z coordinate of the field plane (set by calibration)954 errormsg='inconsistent z position for field and projection plane';955 return956 end957 else % the z coordinate is set only by the field plane (by calibration)958 ProjData.ProjObjectCoord(3)=FieldData.PlaneCoord(3);959 end960 if isfield(FieldData,'PlaneAngle')961 if isfield(ProjData,'ProjObjectAngle')962 if ~isequal(FieldData.PlaneAngle,ProjData.ProjObjectAngle) %check the consistency with the z coordinate of the field plane (set by calibration)963 errormsg='inconsistent plane angle for field and projection plane';964 return965 end966 else967 ProjData.ProjObjectAngle=FieldData.PlaneAngle;968 end951 if length(ProjData.ProjObjectCoord)==3% if the projection plane has a z coordinate 952 if ~isequal(ProjData.PlaneCoord(3),ProjData.ProjObjectCoord) %check the consistency with the z coordinate of the field plane (set by calibration) 953 errormsg='inconsistent z position for field and projection plane'; 954 return 955 end 956 else % the z coordinate is set only by the field plane (by calibration) 957 ProjData.ProjObjectCoord(3)=FieldData.PlaneCoord(3); 958 end 959 if isfield(FieldData,'PlaneAngle') 960 if isfield(ProjData,'ProjObjectAngle') 961 if ~isequal(FieldData.PlaneAngle,ProjData.ProjObjectAngle) %check the consistency with the z coordinate of the field plane (set by calibration) 962 errormsg='inconsistent plane angle for field and projection plane'; 963 return 964 end 965 else 966 ProjData.ProjObjectAngle=FieldData.PlaneAngle; 967 end 969 968 end 970 969 end … … 1009 1008 1010 1009 %% define the new coordinates in case of interpolation on a grid 1011 if check_grid% TODO: recalculate coordinates to get the bounds in the rotated coordinates 1010 if check_grid 1011 AYName='coord_y'; 1012 AXName='coord_x'; 1012 1013 ProjData.ListVarName={'coord_y','coord_x'}; 1013 ProjData.VarDimName={'coord_y','coord_x'}; 1014 ProjData.VarDimName={'coord_y','coord_x'}; 1014 1015 ProjData.VarAttribute={[],[]}; 1015 ProjData.coord_y=[YMin YMax]; 1016 if ~testYMin 1017 errormsg='min Y value not defined for the projection grid';return 1018 end 1019 if ~testYMax 1020 errormsg='max Y value not defined for the projection grid';return 1021 end 1022 if ~testXMin 1023 errormsg='min X value not defined for the projection grid';return 1024 end 1025 if ~testXMax 1026 errormsg='max X value not defined for the projection grid';return 1027 end 1028 ProjData.coord_y=[YMin YMax];%note that if projection is done on a grid, the Min and Max along each direction must have been defined 1016 1029 ProjData.coord_x=[XMin XMax]; 1030 coord_x_proj=XMin:DX:XMax; 1031 coord_y_proj=YMin:DY:YMax; 1032 [X,YI]=meshgrid(coord_x_proj,coord_y_proj);%grid in the new coordinates 1033 XI=ObjectData.Coord(1,1)+(X)*cos(PlaneAngle(3))-YI*sin(PlaneAngle(3));%corresponding coordinates in the original system 1034 YI=ObjectData.Coord(1,2)+(X)*sin(PlaneAngle(3))+YI*cos(PlaneAngle(3)); 1017 1035 end 1018 1036 … … 1020 1038 % LOOP ON GROUPS OF VARIABLES SHARING THE SAME DIMENSIONS 1021 1039 % CellVarIndex=cells of variable index arrays 1022 ivar_new=0; % index of the current variable in the projected field1040 %ivar_new=0; % index of the current variable in the projected field 1023 1041 % icoord=0; 1024 1042 nbcoord=0;%number of added coordinate variables brought by projection 1025 nbvar=0;1043 %nbvar=0; 1026 1044 vector_x_proj=[]; 1027 1045 vector_y_proj=[]; … … 1032 1050 end 1033 1051 VarIndex=CellInfo{icell}.VarIndex;% indices of the selected variables in the list FieldData.ListVarName 1034 ivar_U=[];ivar_V=[];ivar_W=[];1035 if isfield(CellInfo{icell},'VarIndex_vector_x_tps')&&isfield(CellInfo{icell},'VarIndex_vector_y_tps')1036 ivar_U=CellInfo{icell}.VarIndex_vector_x_tps;1037 ivar_V=CellInfo{icell}.VarIndex_vector_y_tps;1038 elseif isfield(CellInfo{icell},'VarIndex_vector_x')&&isfield(CellInfo{icell},'VarIndex_vector_y')1039 ivar_U=CellInfo{icell}.VarIndex_vector_x;1040 ivar_V=CellInfo{icell}.VarIndex_vector_y;1041 end1042 if isfield(CellInfo{icell},'VarIndex_vector_z')1043 ivar_W=CellInfo{icell}.VarIndex_vector_z;1044 end1052 % ivar_U=[];ivar_V=[];ivar_W=[]; 1053 % if isfield(CellInfo{icell},'VarIndex_vector_x_tps')&&isfield(CellInfo{icell},'VarIndex_vector_y_tps') 1054 % ivar_U=CellInfo{icell}.VarIndex_vector_x_tps; 1055 % ivar_V=CellInfo{icell}.VarIndex_vector_y_tps; 1056 % elseif isfield(CellInfo{icell},'VarIndex_vector_x')&&isfield(CellInfo{icell},'VarIndex_vector_y') 1057 % ivar_U=CellInfo{icell}.VarIndex_vector_x; 1058 % ivar_V=CellInfo{icell}.VarIndex_vector_y; 1059 % end 1060 % if isfield(CellInfo{icell},'VarIndex_vector_z') 1061 % ivar_W=CellInfo{icell}.VarIndex_vector_z; 1062 % end 1045 1063 %dimensions 1046 1064 DimCell=FieldData.VarDimName{VarIndex(1)}; … … 1049 1067 end 1050 1068 coord_z=0;%default 1051 1069 ListVarName={};% initiate list of projected variables for cell # icell 1070 VarDimName={};% initiate coresponding list of dimensions for cell # icell 1071 VarAttribute={};% initiate coresponding list of var attributes for cell # icell 1052 1072 %%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1053 1073 switch CellInfo{icell}.CoordType 1054 1074 1055 1075 case 'scattered' 1056 %% case of input fields with unstructured coordinates (applies for projMode ='projection' or 'interp_lin')1076 %% case of input fields with unstructured coordinates (applies for projMode ='projection' or 'interp_lin') 1057 1077 if strcmp(ProjMode{icell},'interp_tps') 1058 1078 continue %skip for next cell (needs tps field cell) 1059 1079 end 1060 coord_x=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(end)});% initial x coordinates 1080 coord_x=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(end)});% initial x coordinates 1061 1081 coord_y=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(end-1)});% initial y coordinates 1062 1082 check3D=(numel(CellInfo{icell}.CoordIndex)==3); … … 1094 1114 coord_Y=(-coord_x *sin(Phi) + coord_y *cos(Phi))*cos(Theta); 1095 1115 coord_Y=coord_Y+coord_z *sin(Theta); 1096 coord_X=(coord_X *cos(Psi) - coord_Y* sin(Psi));%A VERIFIER 1116 coord_X=(coord_X *cos(Psi) - coord_Y* sin(Psi));%A VERIFIER 1097 1117 coord_Y=(coord_X *sin(Psi) + coord_Y* cos(Psi)); 1098 1118 else … … 1137 1157 end 1138 1158 1139 % two cases of projection 1159 % two cases of projection for scattered coordinates 1140 1160 switch ProjMode{icell} 1141 case 'projection' 1142 nbvar=numel(ProjData.ListVarName); 1161 case 'projection' 1162 nbvar=0; 1163 %nbvar=numel(ProjData.ListVarName); 1143 1164 for ivar=VarIndex %transfer variables to the projection plane 1144 1165 VarName=FieldData.ListVarName{ivar}; … … 1151 1172 end 1152 1173 if ~(check3D && ivar==CellInfo{icell}.CoordIndex(1)) 1153 ProjData.ListVarName=[ProjData.ListVarName VarName];1154 ProjData.VarDimName=[ProjData.VarDimName DimCell];1174 ListVarName=[ListVarName VarName]; 1175 VarDimName=[VarDimName DimCell]; 1155 1176 nbvar=nbvar+1; 1156 1177 if isfield(FieldData,'VarAttribute') && length(FieldData.VarAttribute) >=ivar 1157 ProjData.VarAttribute{nbvar}=FieldData.VarAttribute{ivar};1178 VarAttribute{nbvar}=FieldData.VarAttribute{ivar}; 1158 1179 end 1159 1180 end 1160 1181 end 1161 1182 case 'interp_lin'%interpolate data on a regular grid 1162 coord_x_proj=XMin:DX:XMax;1163 coord_y_proj=YMin:DY:YMax;1164 [XI,YI]=meshgrid(coord_x_proj,coord_y_proj);1165 1183 if isfield(CellInfo{icell},'VarIndex_errorflag') 1166 1184 VarName_FF=FieldData.ListVarName{CellInfo{icell}.VarIndex_errorflag}; … … 1174 1192 end 1175 1193 % interpolate and calculate field on the grid 1176 [VarVal,List FieldProj,VarAttribute,errormsg]=calc_field_interp([coord_X coord_Y],FieldData,CellInfo{icell}.FieldName,XI,YI);1194 [VarVal,ListVarName,VarAttribute,errormsg]=calc_field_interp([coord_X coord_Y],FieldData,CellInfo{icell}.FieldName,XI,YI); 1177 1195 1178 1196 if isfield(CellInfo{icell},'CheckSub') && CellInfo{icell}.CheckSub && ~isempty(vector_x_proj) 1179 ProjData.( ProjData.ListVarName{vector_x_proj})=ProjData.(ProjData.ListVarName{vector_x_proj})-VarVal{1};1180 ProjData.( ProjData.ListVarName{vector_y_proj})=ProjData.(ProjData.ListVarName{vector_y_proj})-VarVal{2};1197 ProjData.(ListVarName{vector_x_proj})=ProjData.(ListVarName{vector_x_proj})-VarVal{1}; 1198 ProjData.(ListVarName{vector_y_proj})=ProjData.(ListVarName{vector_y_proj})-VarVal{2}; 1181 1199 else 1182 VarDimName=cell(size(List FieldProj));1183 for ilist=1:numel(List FieldProj)% reshape data, excluding coordinates (ilist=1-2), TODO: rationalise1184 List FieldProj{ilist}=regexprep(ListFieldProj{ilist},'(.+','');1185 if ~isempty(find(strcmp(List FieldProj{ilist},ProjData.ListVarName)))1186 List FieldProj{ilist}=[ListFieldProj{ilist} '_1'];1187 end 1188 ProjData.(List FieldProj{ilist})=VarVal{ilist};1200 VarDimName=cell(size(ListVarName)); 1201 for ilist=1:numel(ListVarName)% reshape data, excluding coordinates (ilist=1-2), TODO: rationalise 1202 ListVarName{ilist}=regexprep(ListVarName{ilist},'(.+',''); 1203 if ~isempty(find(strcmp(ListVarName{ilist},ProjData.ListVarName))) 1204 ListVarName{ilist}=[ListVarName{ilist} '_1']; 1205 end 1206 ProjData.(ListVarName{ilist})=VarVal{ilist}; 1189 1207 VarDimName{ilist}={'coord_y','coord_x'}; 1190 1208 end 1191 ProjData.ListVarName=[ProjData.ListVarName ListFieldProj]; 1192 ProjData.VarDimName=[ProjData.VarDimName VarDimName]; 1193 ProjData.VarAttribute=[ProjData.VarAttribute VarAttribute]; 1194 end 1195 end 1196 1209 end 1210 end 1211 1197 1212 case 'tps' 1198 %% case of tps data (applies only in interp_tps mode)1213 %% case of tps data (applies only in interp_tps mode) 1199 1214 if strcmp(ProjMode{icell},'interp_tps') 1200 1215 Coord=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex}); … … 1204 1219 FieldVar=cat(3,FieldData.(FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_x_tps}),FieldData.(FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_y_tps})); 1205 1220 end 1206 coord_x_proj=XMin:DX:XMax; 1207 coord_y_proj=YMin:DY:YMax; 1208 np_x=numel(coord_x_proj); 1209 np_y=numel(coord_y_proj); 1210 [XI,YI]=meshgrid(coord_x_proj,coord_y_proj'); 1211 XI=XI+ObjectData.Coord(1,1); 1212 YI=YI+ObjectData.Coord(1,2); 1213 [DataOut,VarAttribute,errormsg]=calc_field_tps(Coord,NbCentres,SubRange,FieldVar,CellInfo{icell}.FieldName,cat(3,XI,YI)); 1214 ListFieldProj=(fieldnames(DataOut))'; 1215 VarDimName=cell(size(ListFieldProj)); 1216 for ilist=1:numel(ListFieldProj)% reshape data, excluding coordinates (ilist=1-2), TODO: rationalise 1217 VarName=ListFieldProj{ilist}; 1221 [DataOut,VarAttribute,errormsg]=calc_field_tps(Coord,NbCentres,SubRange,FieldVar,CellInfo{icell}.FieldName,cat(3,XI,YI)); 1222 ListVarName=(fieldnames(DataOut))'; 1223 VarDimName=cell(size(ListVarName)); 1224 for ilist=1:numel(ListVarName)% reshape data, excluding coordinates (ilist=1-2), TODO: rationalise 1225 VarName=ListVarName{ilist}; 1218 1226 ProjData.(VarName)=DataOut.(VarName); 1219 1227 VarDimName{ilist}={'coord_y','coord_x'}; 1220 1228 end 1221 ProjData.ListVarName=[ProjData.ListVarName ListFieldProj];1222 ProjData.VarDimName=[ProjData.VarDimName VarDimName];1223 ProjData.VarAttribute=[ProjData.VarAttribute VarAttribute];1224 1229 end 1225 1230 1231 case 'grid' 1226 1232 %% case of input fields defined on a structured grid 1227 case 'grid'1228 1233 VarName=FieldData.ListVarName{VarIndex(1)};%get the first variable of the cell to get the input matrix dimensions 1229 1234 DimValue=size(FieldData.(VarName));%input matrix dimensions … … 1243 1248 end 1244 1249 end 1245 AYName=FieldData.ListVarName{CellInfo{icell}.CoordIndex(NbDim-1)};%name of input x coordinate (name preserved on projection)1246 AXName=FieldData.ListVarName{CellInfo{icell}.CoordIndex(NbDim)};%name of input y coordinate (name preserved on projection)1247 if testangle% TODO modify name also in case of origin shift in x or y1248 AYProjName='Y';1249 AXProjName='X';1250 count=0;1251 %modify coordinate names if they are already used1252 while ~(isempty(find(strcmp('AXName',ProjData.ListVarName),1)) && isempty(find(strcmp('AYName',ProjData.ListVarName),1)))1253 count=count+1;1254 AYProjName=[AYProjName '_' num2str(count)];1255 AXProjName=[AXProjName '_' num2str(count)];1256 end1257 else1258 AYProjName=AYName;% (name preserved on projection)1259 AXProjName=AXName;%name of input y coordinate (name preserved on projection)1260 end1261 ListDimName=FieldData.VarDimName{VarIndex(1)};1262 ProjData.ListVarName=[ProjData.ListVarName {AYProjName} {AXProjName}]; %TODO: check if it already exists in Projdata (several cells)1263 ProjData.VarDimName=[ProjData.VarDimName {AYProjName} {AXProjName}];1264 ProjData.VarAttribute=[ProjData.VarAttribute {[]} {[]}];1265 1250 Coord_z=[]; 1266 1251 Coord_y=[]; 1267 Coord_x=[]; 1252 Coord_x=[]; 1253 1254 % find default mesh 1268 1255 for idim=1:NbDim %loop on space dimensions 1269 1256 test_interp(idim)=0;%test for coordiate interpolation (non regular grid), =0 by default … … 1273 1260 if numel(Coord{idim})==2 %input array defined on a regular grid 1274 1261 DCoord_min(idim)=(Coord{idim}(2)-Coord{idim}(1))/DimValue(idim); 1262 Coord{idim}=linspace(Coord{idim}(1),Coord{idim}(2),DimValue(idim)); 1275 1263 else 1276 1264 DCoord=diff(Coord{idim});%array of coordinate derivatives for the input field … … 1310 1298 DAX=DCoord_min(NbDim); 1311 1299 DAY=DCoord_min(NbDim-1); 1300 1301 % default bounds 1312 1302 minAX=min(Coord_x); 1313 1303 maxAX=max(Coord_x); 1314 1304 minAY=min(Coord_y); 1315 1305 maxAY=max(Coord_y); 1316 xcorner=[minAX maxAX minAX maxAX]-ObjectData.Coord(1,1); 1306 xcorner=[minAX maxAX minAX maxAX]-ObjectData.Coord(1,1);% image corners in the new coordiantes 1317 1307 ycorner=[maxAY maxAY minAY minAY]-ObjectData.Coord(1,2); 1318 1308 xcor_new=xcorner*cos_om+ycorner*sin_om;%coord new frame … … 1346 1336 test_direct_z=test_direct(1); 1347 1337 end 1348 npX=floor((XMax-XMin)/DX+1); 1349 npY=floor((YMax-YMin)/DY+1); 1350 if test_direct_y 1351 coord_y_proj=linspace(YMin,YMax,npY);%abscissa of the new pixels along the line 1352 else 1353 coord_y_proj=linspace(YMax,YMin,npY);%abscissa of the new pixels along the line 1354 end 1355 if test_direct_x 1356 coord_x_proj=linspace(XMin,XMax,npX);%abscissa of the new pixels along the line 1357 else 1358 coord_x_proj=linspace(XMax,XMin,npX);%abscissa of the new pixels along the line 1359 end 1360 % case with no interpolation 1361 if isequal(ProjMode{icell},'projection') && (~testangle || test90y || test90x) 1338 %define new coordiantes if not already done by the definition of a projection grid 1339 if ~check_grid 1340 AYName=FieldData.ListVarName{CellInfo{icell}.CoordIndex(NbDim-1)};%name of input x coordinate (name preserved on projection) 1341 AXName=FieldData.ListVarName{CellInfo{icell}.CoordIndex(NbDim)};%name of input y coordinate (name preserved on projection) 1342 ListVarName=[ListVarName {AYName} {AXName}];% update the list of projected variables 1343 VarDimName=[VarDimName {AYName} {AXName}];% update the corresponding list of dimensions 1344 VarAttribute=[VarAttribute {[]} {[]}];% update the list of attributes 1345 ProjData.(AYName)=[YMin YMax]; % new (projected ) y coordinates 1346 ProjData.(AXName)=[XMin XMax]; % new (projected ) y coordinates 1347 end 1348 1349 if testangle 1350 ProjMode{icell}='interp_lin'; %request linear interpolation for projection on a tilted plane 1351 end 1352 1353 % npX=floor((XMax-XMin)/DX+1); 1354 % npY=floor((YMax-YMin)/DY+1); 1355 % if test_direct_y 1356 % coord_y_proj=linspace(YMin,YMax,npY);%abscissa of the new pixels along the line 1357 % else 1358 % coord_y_proj=linspace(YMax,YMin,npY);%abscissa of the new pixels along the line 1359 % end 1360 % if test_direct_x 1361 % coord_x_proj=linspace(XMin,XMax,npX);%abscissa of the new pixels along the line 1362 % else 1363 % coord_x_proj=linspace(XMax,XMin,npX);%abscissa of the new pixels along the line 1364 % end 1365 if isequal(ProjMode{icell},'projection')% && (~testangle || test90y || test90x) 1362 1366 if NbDim==2 && ~testXMin && ~testXMax && ~testYMin && ~testYMax% no range restriction 1363 ProjData.ListVarName=[ProjData.ListVarName FieldData.ListVarName(VarIndex)];1364 ProjData.VarDimName=[ProjData.VarDimName FieldData.VarDimName(VarIndex)];1367 ListVarName=[ListVarName FieldData.ListVarName(VarIndex)]; 1368 VarDimName=[VarDimName FieldData.VarDimName(VarIndex)]; 1365 1369 if isfield(FieldData,'VarAttribute') 1366 ProjData.VarAttribute=[ProjData.VarAttribute FieldData.VarAttribute(VarIndex)];1367 end 1368 ProjData.(AY ProjName)=FieldData.(AYName);1369 ProjData.(AX ProjName)=FieldData.(AXName);1370 VarAttribute=[VarAttribute FieldData.VarAttribute(VarIndex)]; 1371 end 1372 ProjData.(AYName)=FieldData.(AYName); 1373 ProjData.(AXName)=FieldData.(AXName); 1370 1374 for ivar=VarIndex 1371 1375 VarName=FieldData.ListVarName{ivar}; … … 1397 1401 end 1398 1402 min_indy=max(min_indy,1);% deals with margin (bound lower than the first index) 1399 min_indx=max(min_indx,1); 1403 min_indx=max(min_indx,1); 1400 1404 if test90y 1401 1405 ind_new=[3 2 1]; … … 1404 1408 for ivar=VarIndex 1405 1409 VarName=FieldData.ListVarName{ivar}; 1406 ProjData.ListVarName=[ProjData.ListVarName VarName];1407 ProjData.VarDimName=[ProjData.VarDimName {DimCell}];1408 ProjData.VarAttribute{length(ProjData.ListVarName)}=FieldData.VarAttribute{ivar}; %reproduce the variable attributes1410 ListVarName=[ListVarName VarName]; 1411 VarDimName=[VarDimName {DimCell}]; 1412 VarAttribute{length(ListVarName)}=FieldData.VarAttribute{ivar}; %reproduce the variable attributes 1409 1413 eval(['ProjData.' VarName '=permute(FieldData.' VarName ',ind_new);'])% permute x and z indices for 90 degree rotation 1410 1414 eval(['ProjData.' VarName '=squeeze(ProjData.' VarName '(iz,:,:));'])% select the z index iz 1411 1415 end 1412 eval(['ProjData.' AY ProjName '=[Ybound(1) Ybound(2)];']) %record the new (projected ) y coordinates1413 eval(['ProjData.' AX ProjName '=[Coord{1}(end),Coord{1}(1)];']) %record the new (projected ) x coordinates1416 eval(['ProjData.' AYName '=[Ybound(1) Ybound(2)];']) %record the new (projected ) y coordinates 1417 eval(['ProjData.' AXName '=[Coord{1}(end),Coord{1}(1)];']) %record the new (projected ) x coordinates 1414 1418 else 1415 1419 if NbDim==3 … … 1426 1430 for ivar=VarIndex% loop on non coordinate variables 1427 1431 VarName=FieldData.ListVarName{ivar}; 1428 ProjData.ListVarName=[ProjData.ListVarName VarName];1429 ProjData.VarDimName=[ProjData.VarDimName {DimCell}];1432 ListVarName=[ListVarName VarName]; 1433 VarDimName=[VarDimName {DimCell}]; 1430 1434 if isfield(FieldData,'VarAttribute') && length(FieldData.VarAttribute)>=ivar 1431 ProjData.VarAttribute{length(ProjData.ListVarName)}=FieldData.VarAttribute{ivar};1435 VarAttribute{length(ListVarName)}=FieldData.VarAttribute{ivar}; 1432 1436 end 1433 1437 if NbDim==3 … … 1437 1441 end 1438 1442 end 1439 eval(['ProjData.' AY ProjName '=[Ybound(1) Ybound(2)];']) %record the new (projected ) y coordinates1440 eval(['ProjData.' AX ProjName '=[Xbound(1) Xbound(2)];']) %record the new (projected ) x coordinates1441 end 1442 end 1443 else % case with rotation and/orinterpolation1443 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 1445 end 1446 end 1447 else % case with interpolation 1444 1448 if NbDim==2 %2D case 1445 [X,Y]=meshgrid(coord_x_proj,coord_y_proj);%grid in the new coordinates1446 XIMA=ObjectData.Coord(1,1)+(X)*cos(PlaneAngle(3))-Y*sin(PlaneAngle(3));%corresponding coordinates in the original image1447 YIMA=ObjectData.Coord(1,2)+(X)*sin(PlaneAngle(3))+Y*cos(PlaneAngle(3));1448 XIMA=(XIMA-minAX)/DXinit+1;% image index along x1449 YIMA=(-YIMA+maxAY)/DYinit+1;% image index along y1450 XIMA=reshape(round(XIMA),1,npX*npY);%indices reorganized in 'line'1451 YIMA=reshape(round(YIMA),1,npX*npY);1452 flagin=XIMA>=1 & XIMA<=DimValue(2) & YIMA >=1 & YIMA<=DimValue(1);%flagin=1 inside the original image1453 1449 if isequal(ProjMode{icell},'interp_tps') 1454 1450 npx_interp_tps=ceil(abs(DX/DAX)); … … 1459 1455 test_interp_tps=0; 1460 1456 end 1461 eval(['ProjData.' AYName '=[coord_y_proj(1) coord_y_proj(end)];']) %record the new (projected ) y coordinates 1462 eval(['ProjData.' AXName '=[coord_x_proj(1) coord_x_proj(end)];']) %record the new (projected ) x coordinates 1457 % ProjData.coord_y=[YMin YMax];%note that if projection is done on a grid, the Min and Max along each direction must have been defined 1458 % ProjData.coord_x=[XMin XMax]; 1459 coord_x_proj=XMin:DX:XMax; 1460 coord_y_proj=YMin:DY:YMax; 1461 [X,YI]=meshgrid(coord_x_proj,coord_y_proj);%grid in the new coordinates 1462 XI=ObjectData.Coord(1,1)+(X)*cos(PlaneAngle(3))-YI*sin(PlaneAngle(3));%corresponding coordinates in the original system 1463 YI=ObjectData.Coord(1,2)+(X)*sin(PlaneAngle(3))+YI*cos(PlaneAngle(3)); 1464 [X,Y]=meshgrid(Coord{2},Coord{1});%initial coordiantes 1463 1465 for ivar=VarIndex 1464 1466 VarName=FieldData.ListVarName{ivar}; 1465 if test_interp(1) || test_interp(2)%interpolate on a regular grid 1466 eval(['ProjData.' VarName '=interp2(Coord{2},Coord{1},FieldData.' VarName ',Coord_x,Coord_y'');']) %TO TEST 1467 % if test_interp(1) || test_interp(2)%interpolate on a regular grid 1468 if size(FieldData.(VarName),3)==1 1469 ProjData.(VarName)=interp2(X,Y,double(FieldData.(VarName)),XI,YI,'*linear'); 1470 else 1471 ProjData.(VarName)=interp2(X,Y,double(FieldData.(VarName)(:,:,1)),XI,YI,'*linear'); 1472 for icolor=2:size(FieldData.(VarName),3) 1473 ProjData.(VarName)=cat(3,ProjData.(VarName),interp2(X,Y,double(FieldData.(VarName)(:,:,icolor)),XI,YI,'*linear')); %TO TEST 1474 end 1467 1475 end 1468 eval(['vec_A=reshape(FieldData.' VarName ',[],nbcolor);'])%put the original image in line 1469 %ind_in=find(flagin); 1470 ind_out=find(~flagin); 1471 ICOMB=(XIMA-1)*DimValue(1)+YIMA; 1472 ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A 1473 vec_B(flagin,1:nbcolor)=vec_A(ICOMB,:); 1474 for icolor=1:nbcolor 1475 vec_B(ind_out,icolor)=zeros(size(ind_out)); 1476 end 1477 ProjData.ListVarName=[ProjData.ListVarName VarName]; 1478 ProjData.VarDimName=[ProjData.VarDimName {DimCell}]; 1476 ListVarName=[ListVarName VarName]; 1477 DimCell(1:2)={AYName,AXName}; 1478 VarDimName=[VarDimName {DimCell}]; 1479 1479 if isfield(FieldData,'VarAttribute')&&length(FieldData.VarAttribute)>=ivar 1480 ProjData.VarAttribute{length(ProjData.ListVarName)+nbcoord}=FieldData.VarAttribute{ivar}; 1481 end 1482 eval(['ProjData.' VarName '=reshape(vec_B,npY,npX,nbcolor);']); 1483 end 1484 ProjData.FF=reshape(~flagin,npY,npX);%false flag A FAIRE: tenir compte d'un flga antï¿œrieur 1485 ProjData.ListVarName=[ProjData.ListVarName 'FF']; 1486 ProjData.VarDimName=[ProjData.VarDimName {DimCell}]; 1487 ProjData.VarAttribute{length(ProjData.ListVarName)}.Role='errorflag'; 1480 VarAttribute{length(ListVarName)+nbcoord}=FieldData.VarAttribute{ivar}; 1481 end; 1482 end 1488 1483 elseif ~testangle 1489 1484 % unstructured z coordinate … … 1496 1491 for ivar=VarIndex 1497 1492 VarName=FieldData.ListVarName{ivar}; 1498 ProjData.ListVarName=[ProjData.ListVarName VarName];1499 ProjData.VarAttribute{length(ProjData.ListVarName)}=FieldData.VarAttribute{ivar}; %reproduce the variable attributes1493 ListVarName=[ListVarName VarName]; 1494 VarAttribute{length(ListVarName)}=FieldData.VarAttribute{ivar}; %reproduce the variable attributes 1500 1495 eval(['ProjData.' VarName '=squeeze(FieldData.' VarName '(iz,:,:));'])% select the z index iz 1501 1496 %TODO : do a vertical average for a thick plane … … 1512 1507 end 1513 1508 end 1509 % update the global list of projected variables: 1510 ProjData.ListVarName=[ProjData.ListVarName ListVarName]; 1511 ProjData.VarDimName=[ProjData.VarDimName VarDimName]; 1512 ProjData.VarAttribute=[ProjData.VarAttribute VarAttribute]; 1514 1513 1515 1514 %% projection of velocity components in the rotated coordinates 1516 if testangle && length(ivar_U)==1 1517 if isempty(ivar_V) 1518 msgbox_uvmat('ERROR','v velocity component missing in proj_field.m') 1519 return 1520 end 1521 UName=FieldData.ListVarName{ivar_U}; 1522 VName=FieldData.ListVarName{ivar_V}; 1523 eval(['ProjData.' UName '=cos(PlaneAngle(3))*ProjData.' UName '+ sin(PlaneAngle(3))*ProjData.' VName ';']) 1524 eval(['ProjData.' VName '=cos(Theta)*(-sin(PlaneAngle(3))*ProjData.' UName '+ cos(PlaneAngle(3))*ProjData.' VName ');']) 1525 if ~isempty(ivar_W) 1526 WName=FieldData.ListVarName{ivar_W}; 1527 eval(['ProjData.' VName '=ProjData.' VName '+ ProjData.' WName '*sin(Theta);'])% 1528 eval(['ProjData.' WName '=NormVec_X*ProjData.' UName '+ NormVec_Y*ProjData.' VName '+ NormVec_Z* ProjData.' WName ';']); 1529 end 1530 if ~isequal(Psi,0) 1531 eval(['ProjData.' UName '=cos(Psi)* ProjData.' UName '- sin(Psi)*ProjData.' VName ';']); 1532 eval(['ProjData.' VName '=sin(Psi)* ProjData.' UName '+ cos(Psi)*ProjData.' VName ';']); 1533 end 1534 end 1535 end 1536 1515 if testangle 1516 ivar_U=[];ivar_V=[];ivar_W=[]; 1517 for ivar=1:numel(VarAttribute) 1518 if isfield(VarAttribute{ivar},'Role') 1519 if strcmp(VarAttribute{ivar}.Role,'vector_x') 1520 ivar_U=ivar; 1521 elseif strcmp(VarAttribute{ivar}.Role,'vector_y') 1522 ivar_V=ivar; 1523 elseif strcmp(VarAttribute{ivar}.Role,'vector_z') 1524 ivar_W=ivar; 1525 end 1526 end 1527 end 1528 if ~isempty(ivar_U) 1529 if isempty(ivar_V) 1530 msgbox_uvmat('ERROR','v velocity component missing in proj_field.m') 1531 return 1532 else 1533 UName=ListVarName{ivar_U}; 1534 VName=ListVarName{ivar_V}; 1535 ProjData.(UName)=cos(PlaneAngle(3))*ProjData.(UName)+ sin(PlaneAngle(3))*ProjData.(VName); 1536 ProjData.(VName)=(-sin(PlaneAngle(3))*ProjData.(UName)+ cos(PlaneAngle(3))*ProjData.(VName)); 1537 if ~isempty(ivar_W) 1538 WName=FieldData.ListVarName{ivar_W}; 1539 eval(['ProjData.' VName '=ProjData.' VName '+ ProjData.' WName '*sin(Theta);'])% 1540 eval(['ProjData.' WName '=NormVec_X*ProjData.' UName '+ NormVec_Y*ProjData.' VName '+ NormVec_Z* ProjData.' WName ';']); 1541 end 1542 % if ~isequal(Psi,0) 1543 % eval(['ProjData.' UName '=cos(Psi)* ProjData.' UName '- sin(Psi)*ProjData.' VName ';']); 1544 % eval(['ProjData.' VName '=sin(Psi)* ProjData.' UName '+ cos(Psi)*ProjData.' VName ';']); 1545 % end 1546 end 1547 end 1548 end 1549 end 1537 1550 % %prepare substraction in case of two input fields 1538 1551 % SubData.ListVarName={}; -
trunk/src/uigetfile_uvmat.m
r651 r653 63 63 uicontrol('Style','pushbutton','Tag','backward','Units','normalized','Position',[0.02 0.77 0.1 0.05],... 64 64 'String','<--','FontWeight','bold','FontUnits','points','FontSize',12,'Callback',@backward,'TooltipString','move backward'); 65 home_button=uicontrol('Style','pushbutton','Units','normalized', 'Position', [0.35 0.005 0.3 0.07],... 66 'String','Work dir','FontWeight','bold','FontUnits','points','FontSize',12,'Callback',@home_dir,'TooltipString','reach the current Matlab working directory'); 67 68 uicontrol('Style','popupmenu','Units','normalized', 'Position', [0.75 0.75 0.23 0.04],'tag','sort_option','Callback',@refresh_GUI,'Visible','off',... 65 home_button=uicontrol('Style','pushbutton','Units','normalized', 'Position', [0.14 0.77 0.2 0.05],... 66 'String','Work dir','FontWeight','bold','FontUnits','points','FontSize',12,'Callback',@home_dir,'TooltipString','reach the current Matlab working directory'); 67 uicontrol('Style','pushbutton','Tag','refresh','Units','normalized','Position', [0.36 0.77 0.2 0.05],'Callback',@refresh_GUI,... 68 'String','Refresh','FontWeight','bold','FontUnits','points','FontSize',12); 69 uicontrol('Style','popupmenu','Units','normalized', 'Position', [0.75 0.74 0.23 0.05],'tag','sort_option','Callback',@refresh_GUI,'Visible','off',... 69 70 'String',{'sort name';'sort date'},'FontUnits','points','FontSize',12,'FontWeight','bold','TooltipString','''sort_option'': sort the files by names or dates'); 70 71 uicontrol('Style','listbox','Units','normalized', 'Position',[0.02 0.08 0.96 0.66], 'Callback', @(src,event)list_Callback(option,FilterExt,src,event),'tag','list',... 71 72 'FontUnits','points','FontSize',12,'TooltipString','''list'':current list of directories, marked by +/, and files'); 72 uicontrol('Style','pushbutton','Units','normalized', 'Position', [0.78 0.005 0.2 0.07],'Callback',@(src,event)close(option,src,event),... 73 'String','Close','FontWeight','bold','FontUnits','points','FontSize',12); 74 uicontrol('Style','pushbutton','Tag','refresh','Units','normalized','Position', [0.02 0.005 0.2 0.07],'Callback',@refresh_GUI,... 75 'String','Refresh','FontWeight','bold','FontUnits','points','FontSize',12); 76 OK_button=uicontrol('Style','pushbutton','Units','normalized', 'Position', [0.75 0.81 0.23 0.07],'BackgroundColor',[0 1 0],... 73 74 OK_button=uicontrol('Style','pushbutton','Units','normalized', 'Position', [0.58 0.005 0.2 0.07],'BackgroundColor',[0 1 0],... 77 75 'String','OK','FontWeight','bold','FontUnits','points','FontSize',12,'Callback',@(src,event)OK_Callback(option,FilterExt,src,event)); 76 close_button=uicontrol('Style','pushbutton','Units','normalized', 'Position', [0.78 0.005 0.2 0.07],'Callback',@(src,event)close(option,src,event),... 77 'FontWeight','bold','FontUnits','points','FontSize',12); 78 78 %set(hrefresh,'UserData',StatusData) 79 79 if strcmp(option,'status_display') %put a run advancement display … … 84 84 % 'String','sort: ','FontUnits','points','FontSize',12,'FontWeight','bold','HorizontalAlignment','right'); 85 85 delete(home_button) 86 set(OK_button,'String','OPEN') 86 set(OK_button,'String','Open') 87 set(close_button,'String','Close') 87 88 elseif strcmp(FilterExt,'uigetdir') %pick a directory 88 89 set(path_title,'String',title); %show the input title for path (directory) 90 set(OK_button,'String','Select') 91 set(close_button,'String','Cancel') 89 92 else %put a title and additional pushbuttons 90 93 uicontrol('Style','text','Units','normalized', 'Position', [0.02 0.74 0.6 0.03],'BackgroundColor',BackgroundColor,... … … 94 97 % uicontrol('Style','text','Units','normalized', 'Position', [0.37 0.8 0.35 0.03],'BackgroundColor',BackgroundColor,... 95 98 % 'String','sort: ','FontUnits','points','FontSize',12,'FontWeight','bold','HorizontalAlignment','right'); 96 99 set(OK_button,'String','Open') 100 set(close_button,'String','Cancel') 97 101 end 98 102 drawnow
Note: See TracChangeset
for help on using the changeset viewer.