Changeset 581 for trunk/src/proj_field.m
- Timestamp:
- Mar 12, 2013, 12:52:13 PM (12 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/src/proj_field.m
r576 r581 142 142 return 143 143 end 144 elseif ~isequal(ObjectData.ProjMode,'interp ')144 elseif ~isequal(ObjectData.ProjMode,'interp_lin') 145 145 errormsg=(['ProjMode option ' ObjectData.ProjMode ' not available in proj_field']); 146 146 return … … 223 223 Var=FieldData.(VarName)(indsel); 224 224 ProjData.(VarName)(ipoint,1)=mean(Var); 225 if isequal(ObjectData.ProjMode,'interp ')225 if isequal(ObjectData.ProjMode,'interp_lin') 226 226 ProjData.(VarName)(ipoint,1)=griddata_uvmat(coord_x(indsel),coord_y(indsel),Var,Xpoint(1),Xpoint(2)); 227 227 end … … 524 524 theta=angle(dlinx+1i*dliny);%angle of each segment 525 525 theta(siz_line(1))=theta(siz_line(1)-1); 526 % determine a rectangles at +-width from the line (only used for the ProjMode='projection or ' filter')526 % determine a rectangles at +-width from the line (only used for the ProjMode='projection or 'interp_tps') 527 527 xsup=zeros(1,siz_line(1)); xinf=zeros(1,siz_line(1)); ysup=zeros(1,siz_line(1)); yinf=zeros(1,siz_line(1)); 528 if isequal(ProjMode,'projection') || isequal(ProjMode,' filter')528 if isequal(ProjMode,'projection') || isequal(ProjMode,'interp_tps') 529 529 xsup(1)=ObjectData.Coord(1,1)-width*sin(theta(1)); 530 530 xinf(1)=ObjectData.Coord(1,1)+width*sin(theta(1)); … … 580 580 % check needed object properties for unstructured positions (position given by the variables with role coord_x, coord_y 581 581 if strcmp(CellInfo{icell}.CoordType,'scattered') 582 if ~isequal(ProjMode,'interp ')582 if ~isequal(ProjMode,'interp_lin') 583 583 if width==0 584 584 errormsg='range of the projection object is missing'; 585 585 return 586 586 else 587 lambda=2/(width*width); %smoothing factor used for filter: weight exp(-2) at distance width from the line587 lambda=2/(width*width); %smoothing factor used for interp_tps: weight exp(-2) at distance width from the line 588 588 end 589 589 end … … 621 621 flagsel=ones(size(coord_x)); 622 622 end 623 if isequal(ProjMode,'projection') | isequal(ProjMode,'filter')623 if isequal(ProjMode,'projection') || isequal(ProjMode,'interp_tps') 624 624 flagsel=flagsel & ((coord_y -yinf(ip))*(xinf(ip+1)-xinf(ip))>(coord_x-xinf(ip))*(yinf(ip+1)-yinf(ip))) ... 625 625 & ((coord_y -ysup(ip))*(xsup(ip+1)-xsup(ip))<(coord_x-xsup(ip))*(ysup(ip+1)-ysup(ip))) ... … … 649 649 end 650 650 end 651 elseif isequal(ProjMode,'interp ') %linear interpolation:651 elseif isequal(ProjMode,'interp_lin') %linear interpolation: 652 652 npoint=floor(linelength/DX)+1;% nbre of points in the profile (interval DX) 653 653 Xproj=linelength/(2*npoint):linelength/npoint:linelength-linelength/(2*npoint); … … 659 659 end 660 660 end 661 elseif isequal(ProjMode,' filter') %filtering661 elseif isequal(ProjMode,'interp_tps') %filtering 662 662 npoint=floor(linelength/DX)+1;% nbre of points in the profile (interval DX) 663 663 Xproj=linelength/(2*npoint):linelength/npoint:linelength-linelength/(2*npoint); … … 983 983 for icell=1:numel(CellInfo)% TODO: recalculate coordinates here to get the bounds in the rotated coordinates 984 984 ProjMode{icell}=ObjectData.ProjMode; 985 if isfield(CellInfo{icell},' FieldRequest')986 switch CellInfo{icell}. FieldRequest985 if isfield(CellInfo{icell},'ProjModeRequest') 986 switch CellInfo{icell}.ProjModeRequest 987 987 case 'interp_lin' 988 ProjMode{icell}='interp ';988 ProjMode{icell}='interp_lin'; 989 989 case 'interp_tps' 990 ProjMode{icell}=' filter';991 end 992 end 993 if strcmp(ProjMode{icell},'interp ')||strcmp(ProjMode{icell},'filter')990 ProjMode{icell}='interp_tps'; 991 end 992 end 993 if strcmp(ProjMode{icell},'interp_lin')||strcmp(ProjMode{icell},'interp_tps') 994 994 check_grid=1; 995 995 end … … 1043 1043 %% case of input fields with unstructured coordinates 1044 1044 case 'scattered' 1045 if strcmp(ProjMode{icell},' filter')1046 continue %skip for filter(needs tps field cell)1045 if strcmp(ProjMode{icell},'interp_tps') 1046 continue %skip for interp_tps (needs tps field cell) 1047 1047 end 1048 1048 coord_x=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(end)});% initial x coordinates … … 1148 1148 end 1149 1149 end 1150 case 'interp '%interpolate data on a regular grid1150 case 'interp_lin'%interpolate data on a regular grid 1151 1151 coord_x_proj=XMin:DX:XMax; 1152 1152 coord_y_proj=YMin:DY:YMax; … … 1205 1205 end 1206 1206 1207 %% case of tps interpolation (applies only in filtermode and for spatial derivatives)1207 %% case of tps interpolation (applies only in interp_tps mode and for spatial derivatives) 1208 1208 case 'tps' 1209 if strcmp(ProjMode{icell},' filter')1209 if strcmp(ProjMode{icell},'interp_tps') 1210 1210 Coord=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex}); 1211 Nb Sites=FieldData.(FieldData.ListVarName{CellInfo{icell}.NbSite_tps});1211 NbCentres=FieldData.(FieldData.ListVarName{CellInfo{icell}.NbCentres_tps}); 1212 1212 SubRange=FieldData.(FieldData.ListVarName{CellInfo{icell}.SubRange_tps}); 1213 1213 if isfield(CellInfo{icell},'VarIndex_vector_x_tps')&&isfield(CellInfo{icell},'VarIndex_vector_y_tps') … … 1221 1221 XI=XI+ObjectData.Coord(1,1); 1222 1222 YI=YI+ObjectData.Coord(1,2); 1223 [DataOut,VarAttribute,errormsg]=calc_field_tps(Coord,Nb Sites,SubRange,FieldVar,CellInfo{icell}.FieldName,cat(3,XI,YI));1223 [DataOut,VarAttribute,errormsg]=calc_field_tps(Coord,NbCentres,SubRange,FieldVar,CellInfo{icell}.FieldName,cat(3,XI,YI)); 1224 1224 ListFieldProj=(fieldnames(DataOut))'; 1225 1225 VarDimName=cell(size(ListFieldProj)); … … 1461 1461 YIMA=reshape(round(YIMA),1,npX*npY); 1462 1462 flagin=XIMA>=1 & XIMA<=DimValue(2) & YIMA >=1 & YIMA<=DimValue(1);%flagin=1 inside the original image 1463 if isequal(ProjMode{icell},' filter')1464 npx_ filter=ceil(abs(DX/DAX));1465 npy_ filter=ceil(abs(DY/DAY));1466 M filter=ones(npy_filter,npx_filter)/(npx_filter*npy_filter);1467 test_ filter=1;1463 if isequal(ProjMode{icell},'interp_tps') 1464 npx_interp_tps=ceil(abs(DX/DAX)); 1465 npy_interp_tps=ceil(abs(DY/DAY)); 1466 Minterp_tps=ones(npy_interp_tps,npx_interp_tps)/(npx_interp_tps*npy_interp_tps); 1467 test_interp_tps=1; 1468 1468 else 1469 test_ filter=0;1469 test_interp_tps=0; 1470 1470 end 1471 1471 eval(['ProjData.' AYName '=[coord_y_proj(1) coord_y_proj(end)];']) %record the new (projected ) y coordinates … … 1476 1476 eval(['ProjData.' VarName '=interp2(Coord{2},Coord{1},FieldData.' VarName ',Coord_x,Coord_y'');']) %TO TEST 1477 1477 end 1478 %filter the field (image) if option 'filter' is used1479 % if test_filter1480 % Aclass=class(FieldData.A);1481 % eval(['ProjData.' VarName '=filter2(Mfilter,FieldData.' VarName ',''valid'');'])1482 % if ~isequal(Aclass,'double')1483 % eval(['ProjData.' VarName '=' Aclass '(FieldData.' VarName ');'])%revert to integer values1484 % end1485 % end1486 1478 eval(['vec_A=reshape(FieldData.' VarName ',[],nbcolor);'])%put the original image in line 1487 1479 %ind_in=find(flagin); … … 1525 1517 else 1526 1518 errormsg='projection of structured coordinates on oblique plane not yet implemented'; 1527 %TODO: use interp 31519 %TODO: use interp_lin3 1528 1520 return 1529 1521 end … … 1801 1793 end 1802 1794 end 1803 elseif isequal(ObjectData.ProjMode,'interp ')||isequal(ObjectData.ProjMode,'filter')%interpolate data on a regular grid1795 elseif isequal(ObjectData.ProjMode,'interp_lin')||isequal(ObjectData.ProjMode,'interp_tps')%interpolate data on a regular grid 1804 1796 coord_x_proj=XMin:DX:XMax; 1805 1797 coord_y_proj=YMin:DY:YMax; … … 2074 2066 YIMA=reshape(round(YIMA),1,npX*npY); 2075 2067 flagin=XIMA>=1 & XIMA<=DimValue(2) & YIMA >=1 & YIMA<=DimValue(1);%flagin=1 inside the original image 2076 if isequal(ObjectData.ProjMode,' filter')2077 npx_ filter=ceil(abs(DX/DAX));2078 npy_ filter=ceil(abs(DY/DAY));2079 M filter=ones(npy_filter,npx_filter)/(npx_filter*npy_filter);2080 test_ filter=1;2068 if isequal(ObjectData.ProjMode,'interp_tps') 2069 npx_interp_tps=ceil(abs(DX/DAX)); 2070 npy_interp_tps=ceil(abs(DY/DAY)); 2071 Minterp_tps=ones(npy_interp_tps,npx_interp_tps)/(npx_interp_tps*npy_interp_tps); 2072 test_interp_tps=1; 2081 2073 else 2082 test_ filter=0;2074 test_interp_tps=0; 2083 2075 end 2084 2076 eval(['ProjData.' AYName '=[coord_y_proj(1) coord_y_proj(end)];']) %record the new (projected ) y coordinates … … 2089 2081 eval(['ProjData.' VarName '=interp2(Coord{2},Coord{1},FieldData.' VarName ',Coord_x,Coord_y'');']) %TO TEST 2090 2082 end 2091 %filter the field (image) if option ' filter' is used2092 if test_ filter2083 %filter the field (image) if option 'interp_tps' is used 2084 if test_interp_tps 2093 2085 Aclass=class(FieldData.A); 2094 ProjData.(VarName)= filter2(Mfilter,FieldData.(VarName),'valid');2086 ProjData.(VarName)=interp_tps2(Minterp_tps,FieldData.(VarName),'valid'); 2095 2087 if ~isequal(Aclass,'double') 2096 2088 ProjData.(VarName)=Aclass(FieldData.(VarName));%revert to integer values
Note: See TracChangeset
for help on using the changeset viewer.