Index: /trunk/src/interp2_uvmat.m =================================================================== --- /trunk/src/interp2_uvmat.m (revision 652) +++ /trunk/src/interp2_uvmat.m (revision 653) @@ -1,3 +1,4 @@ -%'interp2_uvmat': linearly interpolate an image or scalar defined on a regular grid +%'interp2_uvmat': linearly interpolate an image or scalar defined on a +%regular grid USE Matlab interp2 instead %-------------------------------------------------------------------- %OUTPUT: Index: /trunk/src/proj_field.m =================================================================== --- /trunk/src/proj_field.m (revision 652) +++ /trunk/src/proj_field.m (revision 653) @@ -873,6 +873,6 @@ %----------------------------------------------------------------- -%% rotation angles -PlaneAngle=[0 0 0]; +%% rotation angles +PlaneAngle=[0 0 0]; norm_plane=[0 0 1]; cos_om=1; @@ -893,5 +893,5 @@ norm_plane(3)=OmAxis(3)*coeff+cos_om; end -testangle=~isequal(PlaneAngle,[0 0 0]);% && ~test90y && ~test90x;%=1 for slanted plane +testangle=~isequal(PlaneAngle,[0 0 0]);% && ~test90y && ~test90x;%=1 for slanted plane %% mesh sizes DX and DY @@ -899,16 +899,16 @@ DY=[];%default if isfield(ObjectData,'DX') && ~isempty(ObjectData.DX) - DX=abs(ObjectData.DX);%mesh of interpolation points + DX=abs(ObjectData.DX);%mesh of interpolation points elseif isfield(FieldData,'CoordMesh') - DX=FieldData.CoordMesh; + DX=FieldData.CoordMesh; end if isfield(ObjectData,'DY') && ~isempty(ObjectData.DY) - DY=abs(ObjectData.DY);%mesh of interpolation points + DY=abs(ObjectData.DY);%mesh of interpolation points elseif isfield(FieldData,'CoordMesh') - DY=FieldData.CoordMesh; + DY=FieldData.CoordMesh; end if ~strcmp(ObjectData.ProjMode,'projection') && (isempty(DX)||isempty(DY)) - errormsg='DX or DY not defined'; - return + errormsg='DX or DY not defined'; + return end @@ -918,26 +918,25 @@ testYMin=0; testYMax=0; - if isfield(ObjectData,'RangeX') - XMin=min(ObjectData.RangeX); - XMax=max(ObjectData.RangeX); - testXMin=XMax>XMin; - testXMax=1;% range restriction along X -else - XMin=FieldData.XMin;%default -XMax=FieldData.XMax;%default + XMin=min(ObjectData.RangeX); + XMax=max(ObjectData.RangeX); + testXMin=XMax>XMin; + testXMax=1;% range restriction along X + % else + % XMin=FieldData.XMin;%default + % XMax=FieldData.XMax;%default end if isfield(ObjectData,'RangeY') - YMin=min(ObjectData.RangeY); - YMax=max(ObjectData.RangeY); - testYMin=YMax>YMin; - testYMax=1; -else - YMin=FieldData.YMin;%default -YMax=FieldData.YMax;%default + YMin=min(ObjectData.RangeY); + YMax=max(ObjectData.RangeY); + testYMin=YMax>YMin; + testYMax=1; + % else + % YMin=FieldData.YMin;%default + % YMax=FieldData.YMax;%default end width=0;%default width of the projection band if isfield(ObjectData,'RangeZ') - width=max(ObjectData.RangeZ); + width=max(ObjectData.RangeZ); end @@ -950,21 +949,21 @@ %% reproduce initial plane position and angle if isfield(FieldData,'PlaneCoord')&&length(FieldData.PlaneCoord)==3&& isfield(ProjData,'ProjObjectCoord') - if length(ProjData.ProjObjectCoord)==3% if the projection plane has a z coordinate - if ~isequal(ProjData.PlaneCoord(3),ProjData.ProjObjectCoord) %check the consistency with the z coordinate of the field plane (set by calibration) - errormsg='inconsistent z position for field and projection plane'; - return - end - else % the z coordinate is set only by the field plane (by calibration) - ProjData.ProjObjectCoord(3)=FieldData.PlaneCoord(3); - end - if isfield(FieldData,'PlaneAngle') - if isfield(ProjData,'ProjObjectAngle') - if ~isequal(FieldData.PlaneAngle,ProjData.ProjObjectAngle) %check the consistency with the z coordinate of the field plane (set by calibration) - errormsg='inconsistent plane angle for field and projection plane'; - return - end - else - ProjData.ProjObjectAngle=FieldData.PlaneAngle; - end + if length(ProjData.ProjObjectCoord)==3% if the projection plane has a z coordinate + if ~isequal(ProjData.PlaneCoord(3),ProjData.ProjObjectCoord) %check the consistency with the z coordinate of the field plane (set by calibration) + errormsg='inconsistent z position for field and projection plane'; + return + end + else % the z coordinate is set only by the field plane (by calibration) + ProjData.ProjObjectCoord(3)=FieldData.PlaneCoord(3); + end + if isfield(FieldData,'PlaneAngle') + if isfield(ProjData,'ProjObjectAngle') + if ~isequal(FieldData.PlaneAngle,ProjData.ProjObjectAngle) %check the consistency with the z coordinate of the field plane (set by calibration) + errormsg='inconsistent plane angle for field and projection plane'; + return + end + else + ProjData.ProjObjectAngle=FieldData.PlaneAngle; + end end end @@ -1009,10 +1008,29 @@ %% define the new coordinates in case of interpolation on a grid -if check_grid% TODO: recalculate coordinates to get the bounds in the rotated coordinates +if check_grid + AYName='coord_y'; + AXName='coord_x'; ProjData.ListVarName={'coord_y','coord_x'}; - ProjData.VarDimName={'coord_y','coord_x'}; + ProjData.VarDimName={'coord_y','coord_x'}; ProjData.VarAttribute={[],[]}; - ProjData.coord_y=[YMin YMax]; + if ~testYMin + errormsg='min Y value not defined for the projection grid';return + end + if ~testYMax + errormsg='max Y value not defined for the projection grid';return + end + if ~testXMin + errormsg='min X value not defined for the projection grid';return + end + if ~testXMax + errormsg='max X value not defined for the projection grid';return + end + 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 ProjData.coord_x=[XMin XMax]; + coord_x_proj=XMin:DX:XMax; + coord_y_proj=YMin:DY:YMax; + [X,YI]=meshgrid(coord_x_proj,coord_y_proj);%grid in the new coordinates + XI=ObjectData.Coord(1,1)+(X)*cos(PlaneAngle(3))-YI*sin(PlaneAngle(3));%corresponding coordinates in the original system + YI=ObjectData.Coord(1,2)+(X)*sin(PlaneAngle(3))+YI*cos(PlaneAngle(3)); end @@ -1020,8 +1038,8 @@ % LOOP ON GROUPS OF VARIABLES SHARING THE SAME DIMENSIONS % CellVarIndex=cells of variable index arrays -ivar_new=0; % index of the current variable in the projected field +%ivar_new=0; % index of the current variable in the projected field % icoord=0; nbcoord=0;%number of added coordinate variables brought by projection -nbvar=0; +%nbvar=0; vector_x_proj=[]; vector_y_proj=[]; @@ -1032,15 +1050,15 @@ end VarIndex=CellInfo{icell}.VarIndex;% indices of the selected variables in the list FieldData.ListVarName - ivar_U=[];ivar_V=[];ivar_W=[]; - if isfield(CellInfo{icell},'VarIndex_vector_x_tps')&&isfield(CellInfo{icell},'VarIndex_vector_y_tps') - ivar_U=CellInfo{icell}.VarIndex_vector_x_tps; - ivar_V=CellInfo{icell}.VarIndex_vector_y_tps; - elseif isfield(CellInfo{icell},'VarIndex_vector_x')&&isfield(CellInfo{icell},'VarIndex_vector_y') - ivar_U=CellInfo{icell}.VarIndex_vector_x; - ivar_V=CellInfo{icell}.VarIndex_vector_y; - end - if isfield(CellInfo{icell},'VarIndex_vector_z') - ivar_W=CellInfo{icell}.VarIndex_vector_z; - end +% ivar_U=[];ivar_V=[];ivar_W=[]; +% if isfield(CellInfo{icell},'VarIndex_vector_x_tps')&&isfield(CellInfo{icell},'VarIndex_vector_y_tps') +% ivar_U=CellInfo{icell}.VarIndex_vector_x_tps; +% ivar_V=CellInfo{icell}.VarIndex_vector_y_tps; +% elseif isfield(CellInfo{icell},'VarIndex_vector_x')&&isfield(CellInfo{icell},'VarIndex_vector_y') +% ivar_U=CellInfo{icell}.VarIndex_vector_x; +% ivar_V=CellInfo{icell}.VarIndex_vector_y; +% end +% if isfield(CellInfo{icell},'VarIndex_vector_z') +% ivar_W=CellInfo{icell}.VarIndex_vector_z; +% end %dimensions DimCell=FieldData.VarDimName{VarIndex(1)}; @@ -1049,14 +1067,16 @@ end coord_z=0;%default - + ListVarName={};% initiate list of projected variables for cell # icell + VarDimName={};% initiate coresponding list of dimensions for cell # icell + VarAttribute={};% initiate coresponding list of var attributes for cell # icell %%%%%%%%%%%%%%%%%%%%%%%%%%%%% switch CellInfo{icell}.CoordType case 'scattered' - %% case of input fields with unstructured coordinates (applies for projMode ='projection' or 'interp_lin') + %% case of input fields with unstructured coordinates (applies for projMode ='projection' or 'interp_lin') if strcmp(ProjMode{icell},'interp_tps') continue %skip for next cell (needs tps field cell) end - coord_x=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(end)});% initial x coordinates + coord_x=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(end)});% initial x coordinates coord_y=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(end-1)});% initial y coordinates check3D=(numel(CellInfo{icell}.CoordIndex)==3); @@ -1094,5 +1114,5 @@ coord_Y=(-coord_x *sin(Phi) + coord_y *cos(Phi))*cos(Theta); coord_Y=coord_Y+coord_z *sin(Theta); - coord_X=(coord_X *cos(Psi) - coord_Y* sin(Psi));%A VERIFIER + coord_X=(coord_X *cos(Psi) - coord_Y* sin(Psi));%A VERIFIER coord_Y=(coord_X *sin(Psi) + coord_Y* cos(Psi)); else @@ -1137,8 +1157,9 @@ end - % two cases of projection + % two cases of projection for scattered coordinates switch ProjMode{icell} - case 'projection' - nbvar=numel(ProjData.ListVarName); + case 'projection' + nbvar=0; + %nbvar=numel(ProjData.ListVarName); for ivar=VarIndex %transfer variables to the projection plane VarName=FieldData.ListVarName{ivar}; @@ -1151,16 +1172,13 @@ end if ~(check3D && ivar==CellInfo{icell}.CoordIndex(1)) - ProjData.ListVarName=[ProjData.ListVarName VarName]; - ProjData.VarDimName=[ProjData.VarDimName DimCell]; + ListVarName=[ListVarName VarName]; + VarDimName=[VarDimName DimCell]; nbvar=nbvar+1; if isfield(FieldData,'VarAttribute') && length(FieldData.VarAttribute) >=ivar - ProjData.VarAttribute{nbvar}=FieldData.VarAttribute{ivar}; + VarAttribute{nbvar}=FieldData.VarAttribute{ivar}; end end end case 'interp_lin'%interpolate data on a regular grid - coord_x_proj=XMin:DX:XMax; - coord_y_proj=YMin:DY:YMax; - [XI,YI]=meshgrid(coord_x_proj,coord_y_proj); if isfield(CellInfo{icell},'VarIndex_errorflag') VarName_FF=FieldData.ListVarName{CellInfo{icell}.VarIndex_errorflag}; @@ -1174,27 +1192,24 @@ end % interpolate and calculate field on the grid - [VarVal,ListFieldProj,VarAttribute,errormsg]=calc_field_interp([coord_X coord_Y],FieldData,CellInfo{icell}.FieldName,XI,YI); + [VarVal,ListVarName,VarAttribute,errormsg]=calc_field_interp([coord_X coord_Y],FieldData,CellInfo{icell}.FieldName,XI,YI); if isfield(CellInfo{icell},'CheckSub') && CellInfo{icell}.CheckSub && ~isempty(vector_x_proj) - ProjData.(ProjData.ListVarName{vector_x_proj})=ProjData.(ProjData.ListVarName{vector_x_proj})-VarVal{1}; - ProjData.(ProjData.ListVarName{vector_y_proj})=ProjData.(ProjData.ListVarName{vector_y_proj})-VarVal{2}; + ProjData.(ListVarName{vector_x_proj})=ProjData.(ListVarName{vector_x_proj})-VarVal{1}; + ProjData.(ListVarName{vector_y_proj})=ProjData.(ListVarName{vector_y_proj})-VarVal{2}; else - VarDimName=cell(size(ListFieldProj)); - for ilist=1:numel(ListFieldProj)% reshape data, excluding coordinates (ilist=1-2), TODO: rationalise - ListFieldProj{ilist}=regexprep(ListFieldProj{ilist},'(.+',''); - if ~isempty(find(strcmp(ListFieldProj{ilist},ProjData.ListVarName))) - ListFieldProj{ilist}=[ListFieldProj{ilist} '_1']; - end - ProjData.(ListFieldProj{ilist})=VarVal{ilist}; + VarDimName=cell(size(ListVarName)); + for ilist=1:numel(ListVarName)% reshape data, excluding coordinates (ilist=1-2), TODO: rationalise + ListVarName{ilist}=regexprep(ListVarName{ilist},'(.+',''); + if ~isempty(find(strcmp(ListVarName{ilist},ProjData.ListVarName))) + ListVarName{ilist}=[ListVarName{ilist} '_1']; + end + ProjData.(ListVarName{ilist})=VarVal{ilist}; VarDimName{ilist}={'coord_y','coord_x'}; end - ProjData.ListVarName=[ProjData.ListVarName ListFieldProj]; - ProjData.VarDimName=[ProjData.VarDimName VarDimName]; - ProjData.VarAttribute=[ProjData.VarAttribute VarAttribute]; - end - end - + end + end + case 'tps' - %% case of tps data (applies only in interp_tps mode) + %% case of tps data (applies only in interp_tps mode) if strcmp(ProjMode{icell},'interp_tps') Coord=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex}); @@ -1204,26 +1219,16 @@ FieldVar=cat(3,FieldData.(FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_x_tps}),FieldData.(FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_y_tps})); end - coord_x_proj=XMin:DX:XMax; - coord_y_proj=YMin:DY:YMax; - np_x=numel(coord_x_proj); - np_y=numel(coord_y_proj); - [XI,YI]=meshgrid(coord_x_proj,coord_y_proj'); - XI=XI+ObjectData.Coord(1,1); - YI=YI+ObjectData.Coord(1,2); - [DataOut,VarAttribute,errormsg]=calc_field_tps(Coord,NbCentres,SubRange,FieldVar,CellInfo{icell}.FieldName,cat(3,XI,YI)); - ListFieldProj=(fieldnames(DataOut))'; - VarDimName=cell(size(ListFieldProj)); - for ilist=1:numel(ListFieldProj)% reshape data, excluding coordinates (ilist=1-2), TODO: rationalise - VarName=ListFieldProj{ilist}; + [DataOut,VarAttribute,errormsg]=calc_field_tps(Coord,NbCentres,SubRange,FieldVar,CellInfo{icell}.FieldName,cat(3,XI,YI)); + ListVarName=(fieldnames(DataOut))'; + VarDimName=cell(size(ListVarName)); + for ilist=1:numel(ListVarName)% reshape data, excluding coordinates (ilist=1-2), TODO: rationalise + VarName=ListVarName{ilist}; ProjData.(VarName)=DataOut.(VarName); VarDimName{ilist}={'coord_y','coord_x'}; end - ProjData.ListVarName=[ProjData.ListVarName ListFieldProj]; - ProjData.VarDimName=[ProjData.VarDimName VarDimName]; - ProjData.VarAttribute=[ProjData.VarAttribute VarAttribute]; end + case 'grid' %% case of input fields defined on a structured grid - case 'grid' VarName=FieldData.ListVarName{VarIndex(1)};%get the first variable of the cell to get the input matrix dimensions DimValue=size(FieldData.(VarName));%input matrix dimensions @@ -1243,27 +1248,9 @@ end end - AYName=FieldData.ListVarName{CellInfo{icell}.CoordIndex(NbDim-1)};%name of input x coordinate (name preserved on projection) - AXName=FieldData.ListVarName{CellInfo{icell}.CoordIndex(NbDim)};%name of input y coordinate (name preserved on projection) - if testangle% TODO modify name also in case of origin shift in x or y - AYProjName='Y'; - AXProjName='X'; - count=0; - %modify coordinate names if they are already used - while ~(isempty(find(strcmp('AXName',ProjData.ListVarName),1)) && isempty(find(strcmp('AYName',ProjData.ListVarName),1))) - count=count+1; - AYProjName=[AYProjName '_' num2str(count)]; - AXProjName=[AXProjName '_' num2str(count)]; - end - else - AYProjName=AYName;% (name preserved on projection) - AXProjName=AXName;%name of input y coordinate (name preserved on projection) - end - ListDimName=FieldData.VarDimName{VarIndex(1)}; - ProjData.ListVarName=[ProjData.ListVarName {AYProjName} {AXProjName}]; %TODO: check if it already exists in Projdata (several cells) - ProjData.VarDimName=[ProjData.VarDimName {AYProjName} {AXProjName}]; - ProjData.VarAttribute=[ProjData.VarAttribute {[]} {[]}]; Coord_z=[]; Coord_y=[]; - Coord_x=[]; + Coord_x=[]; + + % find default mesh for idim=1:NbDim %loop on space dimensions test_interp(idim)=0;%test for coordiate interpolation (non regular grid), =0 by default @@ -1273,4 +1260,5 @@ if numel(Coord{idim})==2 %input array defined on a regular grid DCoord_min(idim)=(Coord{idim}(2)-Coord{idim}(1))/DimValue(idim); + Coord{idim}=linspace(Coord{idim}(1),Coord{idim}(2),DimValue(idim)); else DCoord=diff(Coord{idim});%array of coordinate derivatives for the input field @@ -1310,9 +1298,11 @@ DAX=DCoord_min(NbDim); DAY=DCoord_min(NbDim-1); + + % default bounds minAX=min(Coord_x); maxAX=max(Coord_x); minAY=min(Coord_y); maxAY=max(Coord_y); - xcorner=[minAX maxAX minAX maxAX]-ObjectData.Coord(1,1); + xcorner=[minAX maxAX minAX maxAX]-ObjectData.Coord(1,1);% image corners in the new coordiantes ycorner=[maxAY maxAY minAY minAY]-ObjectData.Coord(1,2); xcor_new=xcorner*cos_om+ycorner*sin_om;%coord new frame @@ -1346,26 +1336,40 @@ test_direct_z=test_direct(1); end - npX=floor((XMax-XMin)/DX+1); - npY=floor((YMax-YMin)/DY+1); - if test_direct_y - coord_y_proj=linspace(YMin,YMax,npY);%abscissa of the new pixels along the line - else - coord_y_proj=linspace(YMax,YMin,npY);%abscissa of the new pixels along the line - end - if test_direct_x - coord_x_proj=linspace(XMin,XMax,npX);%abscissa of the new pixels along the line - else - coord_x_proj=linspace(XMax,XMin,npX);%abscissa of the new pixels along the line - end - % case with no interpolation - if isequal(ProjMode{icell},'projection') && (~testangle || test90y || test90x) + %define new coordiantes if not already done by the definition of a projection grid + if ~check_grid + AYName=FieldData.ListVarName{CellInfo{icell}.CoordIndex(NbDim-1)};%name of input x coordinate (name preserved on projection) + AXName=FieldData.ListVarName{CellInfo{icell}.CoordIndex(NbDim)};%name of input y coordinate (name preserved on projection) + ListVarName=[ListVarName {AYName} {AXName}];% update the list of projected variables + VarDimName=[VarDimName {AYName} {AXName}];% update the corresponding list of dimensions + VarAttribute=[VarAttribute {[]} {[]}];% update the list of attributes + ProjData.(AYName)=[YMin YMax]; % new (projected ) y coordinates + ProjData.(AXName)=[XMin XMax]; % new (projected ) y coordinates + end + + if testangle + ProjMode{icell}='interp_lin'; %request linear interpolation for projection on a tilted plane + end + + % npX=floor((XMax-XMin)/DX+1); + % npY=floor((YMax-YMin)/DY+1); + % if test_direct_y + % coord_y_proj=linspace(YMin,YMax,npY);%abscissa of the new pixels along the line + % else + % coord_y_proj=linspace(YMax,YMin,npY);%abscissa of the new pixels along the line + % end + % if test_direct_x + % coord_x_proj=linspace(XMin,XMax,npX);%abscissa of the new pixels along the line + % else + % coord_x_proj=linspace(XMax,XMin,npX);%abscissa of the new pixels along the line + % end + if isequal(ProjMode{icell},'projection')% && (~testangle || test90y || test90x) if NbDim==2 && ~testXMin && ~testXMax && ~testYMin && ~testYMax% no range restriction - ProjData.ListVarName=[ProjData.ListVarName FieldData.ListVarName(VarIndex)]; - ProjData.VarDimName=[ProjData.VarDimName FieldData.VarDimName(VarIndex)]; + ListVarName=[ListVarName FieldData.ListVarName(VarIndex)]; + VarDimName=[VarDimName FieldData.VarDimName(VarIndex)]; if isfield(FieldData,'VarAttribute') - ProjData.VarAttribute=[ProjData.VarAttribute FieldData.VarAttribute(VarIndex)]; - end - ProjData.(AYProjName)=FieldData.(AYName); - ProjData.(AXProjName)=FieldData.(AXName); + VarAttribute=[VarAttribute FieldData.VarAttribute(VarIndex)]; + end + ProjData.(AYName)=FieldData.(AYName); + ProjData.(AXName)=FieldData.(AXName); for ivar=VarIndex VarName=FieldData.ListVarName{ivar}; @@ -1397,5 +1401,5 @@ end min_indy=max(min_indy,1);% deals with margin (bound lower than the first index) - min_indx=max(min_indx,1); + min_indx=max(min_indx,1); if test90y ind_new=[3 2 1]; @@ -1404,12 +1408,12 @@ for ivar=VarIndex VarName=FieldData.ListVarName{ivar}; - ProjData.ListVarName=[ProjData.ListVarName VarName]; - ProjData.VarDimName=[ProjData.VarDimName {DimCell}]; - ProjData.VarAttribute{length(ProjData.ListVarName)}=FieldData.VarAttribute{ivar}; %reproduce the variable attributes + ListVarName=[ListVarName VarName]; + VarDimName=[VarDimName {DimCell}]; + VarAttribute{length(ListVarName)}=FieldData.VarAttribute{ivar}; %reproduce the variable attributes eval(['ProjData.' VarName '=permute(FieldData.' VarName ',ind_new);'])% permute x and z indices for 90 degree rotation eval(['ProjData.' VarName '=squeeze(ProjData.' VarName '(iz,:,:));'])% select the z index iz end - eval(['ProjData.' AYProjName '=[Ybound(1) Ybound(2)];']) %record the new (projected ) y coordinates - eval(['ProjData.' AXProjName '=[Coord{1}(end),Coord{1}(1)];']) %record the new (projected ) x coordinates + eval(['ProjData.' AYName '=[Ybound(1) Ybound(2)];']) %record the new (projected ) y coordinates + eval(['ProjData.' AXName '=[Coord{1}(end),Coord{1}(1)];']) %record the new (projected ) x coordinates else if NbDim==3 @@ -1426,8 +1430,8 @@ for ivar=VarIndex% loop on non coordinate variables VarName=FieldData.ListVarName{ivar}; - ProjData.ListVarName=[ProjData.ListVarName VarName]; - ProjData.VarDimName=[ProjData.VarDimName {DimCell}]; + ListVarName=[ListVarName VarName]; + VarDimName=[VarDimName {DimCell}]; if isfield(FieldData,'VarAttribute') && length(FieldData.VarAttribute)>=ivar - ProjData.VarAttribute{length(ProjData.ListVarName)}=FieldData.VarAttribute{ivar}; + VarAttribute{length(ListVarName)}=FieldData.VarAttribute{ivar}; end if NbDim==3 @@ -1437,18 +1441,10 @@ end end - eval(['ProjData.' AYProjName '=[Ybound(1) Ybound(2)];']) %record the new (projected ) y coordinates - eval(['ProjData.' AXProjName '=[Xbound(1) Xbound(2)];']) %record the new (projected ) x coordinates - end - end - else % case with rotation and/or interpolation + eval(['ProjData.' AYName '=[Ybound(1) Ybound(2)];']) %record the new (projected ) y coordinates + eval(['ProjData.' AXName '=[Xbound(1) Xbound(2)];']) %record the new (projected ) x coordinates + end + end + else % case with interpolation if NbDim==2 %2D case - [X,Y]=meshgrid(coord_x_proj,coord_y_proj);%grid in the new coordinates - XIMA=ObjectData.Coord(1,1)+(X)*cos(PlaneAngle(3))-Y*sin(PlaneAngle(3));%corresponding coordinates in the original image - YIMA=ObjectData.Coord(1,2)+(X)*sin(PlaneAngle(3))+Y*cos(PlaneAngle(3)); - XIMA=(XIMA-minAX)/DXinit+1;% image index along x - YIMA=(-YIMA+maxAY)/DYinit+1;% image index along y - XIMA=reshape(round(XIMA),1,npX*npY);%indices reorganized in 'line' - YIMA=reshape(round(YIMA),1,npX*npY); - flagin=XIMA>=1 & XIMA<=DimValue(2) & YIMA >=1 & YIMA<=DimValue(1);%flagin=1 inside the original image if isequal(ProjMode{icell},'interp_tps') npx_interp_tps=ceil(abs(DX/DAX)); @@ -1459,31 +1455,30 @@ test_interp_tps=0; end - eval(['ProjData.' AYName '=[coord_y_proj(1) coord_y_proj(end)];']) %record the new (projected ) y coordinates - eval(['ProjData.' AXName '=[coord_x_proj(1) coord_x_proj(end)];']) %record the new (projected ) x coordinates + % 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 + % ProjData.coord_x=[XMin XMax]; + coord_x_proj=XMin:DX:XMax; + coord_y_proj=YMin:DY:YMax; + [X,YI]=meshgrid(coord_x_proj,coord_y_proj);%grid in the new coordinates + XI=ObjectData.Coord(1,1)+(X)*cos(PlaneAngle(3))-YI*sin(PlaneAngle(3));%corresponding coordinates in the original system + YI=ObjectData.Coord(1,2)+(X)*sin(PlaneAngle(3))+YI*cos(PlaneAngle(3)); + [X,Y]=meshgrid(Coord{2},Coord{1});%initial coordiantes for ivar=VarIndex VarName=FieldData.ListVarName{ivar}; - if test_interp(1) || test_interp(2)%interpolate on a regular grid - eval(['ProjData.' VarName '=interp2(Coord{2},Coord{1},FieldData.' VarName ',Coord_x,Coord_y'');']) %TO TEST + % if test_interp(1) || test_interp(2)%interpolate on a regular grid + if size(FieldData.(VarName),3)==1 + ProjData.(VarName)=interp2(X,Y,double(FieldData.(VarName)),XI,YI,'*linear'); + else + ProjData.(VarName)=interp2(X,Y,double(FieldData.(VarName)(:,:,1)),XI,YI,'*linear'); + for icolor=2:size(FieldData.(VarName),3) + ProjData.(VarName)=cat(3,ProjData.(VarName),interp2(X,Y,double(FieldData.(VarName)(:,:,icolor)),XI,YI,'*linear')); %TO TEST + end end - eval(['vec_A=reshape(FieldData.' VarName ',[],nbcolor);'])%put the original image in line - %ind_in=find(flagin); - ind_out=find(~flagin); - ICOMB=(XIMA-1)*DimValue(1)+YIMA; - ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A - vec_B(flagin,1:nbcolor)=vec_A(ICOMB,:); - for icolor=1:nbcolor - vec_B(ind_out,icolor)=zeros(size(ind_out)); - end - ProjData.ListVarName=[ProjData.ListVarName VarName]; - ProjData.VarDimName=[ProjData.VarDimName {DimCell}]; + ListVarName=[ListVarName VarName]; + DimCell(1:2)={AYName,AXName}; + VarDimName=[VarDimName {DimCell}]; if isfield(FieldData,'VarAttribute')&&length(FieldData.VarAttribute)>=ivar - ProjData.VarAttribute{length(ProjData.ListVarName)+nbcoord}=FieldData.VarAttribute{ivar}; - end - eval(['ProjData.' VarName '=reshape(vec_B,npY,npX,nbcolor);']); - end - ProjData.FF=reshape(~flagin,npY,npX);%false flag A FAIRE: tenir compte d'un flga antï¿œrieur - ProjData.ListVarName=[ProjData.ListVarName 'FF']; - ProjData.VarDimName=[ProjData.VarDimName {DimCell}]; - ProjData.VarAttribute{length(ProjData.ListVarName)}.Role='errorflag'; + VarAttribute{length(ListVarName)+nbcoord}=FieldData.VarAttribute{ivar}; + end; + end elseif ~testangle % unstructured z coordinate @@ -1496,6 +1491,6 @@ for ivar=VarIndex VarName=FieldData.ListVarName{ivar}; - ProjData.ListVarName=[ProjData.ListVarName VarName]; - ProjData.VarAttribute{length(ProjData.ListVarName)}=FieldData.VarAttribute{ivar}; %reproduce the variable attributes + ListVarName=[ListVarName VarName]; + VarAttribute{length(ListVarName)}=FieldData.VarAttribute{ivar}; %reproduce the variable attributes eval(['ProjData.' VarName '=squeeze(FieldData.' VarName '(iz,:,:));'])% select the z index iz %TODO : do a vertical average for a thick plane @@ -1512,27 +1507,45 @@ end end + % update the global list of projected variables: + ProjData.ListVarName=[ProjData.ListVarName ListVarName]; + ProjData.VarDimName=[ProjData.VarDimName VarDimName]; + ProjData.VarAttribute=[ProjData.VarAttribute VarAttribute]; %% projection of velocity components in the rotated coordinates - if testangle && length(ivar_U)==1 - if isempty(ivar_V) - msgbox_uvmat('ERROR','v velocity component missing in proj_field.m') - return - end - UName=FieldData.ListVarName{ivar_U}; - VName=FieldData.ListVarName{ivar_V}; - eval(['ProjData.' UName '=cos(PlaneAngle(3))*ProjData.' UName '+ sin(PlaneAngle(3))*ProjData.' VName ';']) - eval(['ProjData.' VName '=cos(Theta)*(-sin(PlaneAngle(3))*ProjData.' UName '+ cos(PlaneAngle(3))*ProjData.' VName ');']) - if ~isempty(ivar_W) - WName=FieldData.ListVarName{ivar_W}; - eval(['ProjData.' VName '=ProjData.' VName '+ ProjData.' WName '*sin(Theta);'])% - eval(['ProjData.' WName '=NormVec_X*ProjData.' UName '+ NormVec_Y*ProjData.' VName '+ NormVec_Z* ProjData.' WName ';']); - end - if ~isequal(Psi,0) - eval(['ProjData.' UName '=cos(Psi)* ProjData.' UName '- sin(Psi)*ProjData.' VName ';']); - eval(['ProjData.' VName '=sin(Psi)* ProjData.' UName '+ cos(Psi)*ProjData.' VName ';']); - end - end -end - + if testangle + ivar_U=[];ivar_V=[];ivar_W=[]; + for ivar=1:numel(VarAttribute) + if isfield(VarAttribute{ivar},'Role') + if strcmp(VarAttribute{ivar}.Role,'vector_x') + ivar_U=ivar; + elseif strcmp(VarAttribute{ivar}.Role,'vector_y') + ivar_V=ivar; + elseif strcmp(VarAttribute{ivar}.Role,'vector_z') + ivar_W=ivar; + end + end + end + if ~isempty(ivar_U) + if isempty(ivar_V) + msgbox_uvmat('ERROR','v velocity component missing in proj_field.m') + return + else + UName=ListVarName{ivar_U}; + VName=ListVarName{ivar_V}; + ProjData.(UName)=cos(PlaneAngle(3))*ProjData.(UName)+ sin(PlaneAngle(3))*ProjData.(VName); + ProjData.(VName)=(-sin(PlaneAngle(3))*ProjData.(UName)+ cos(PlaneAngle(3))*ProjData.(VName)); + if ~isempty(ivar_W) + WName=FieldData.ListVarName{ivar_W}; + eval(['ProjData.' VName '=ProjData.' VName '+ ProjData.' WName '*sin(Theta);'])% + eval(['ProjData.' WName '=NormVec_X*ProjData.' UName '+ NormVec_Y*ProjData.' VName '+ NormVec_Z* ProjData.' WName ';']); + end +% if ~isequal(Psi,0) +% eval(['ProjData.' UName '=cos(Psi)* ProjData.' UName '- sin(Psi)*ProjData.' VName ';']); +% eval(['ProjData.' VName '=sin(Psi)* ProjData.' UName '+ cos(Psi)*ProjData.' VName ';']); +% end + end + end + end +end % %prepare substraction in case of two input fields % SubData.ListVarName={}; Index: /trunk/src/uigetfile_uvmat.m =================================================================== --- /trunk/src/uigetfile_uvmat.m (revision 652) +++ /trunk/src/uigetfile_uvmat.m (revision 653) @@ -63,17 +63,17 @@ uicontrol('Style','pushbutton','Tag','backward','Units','normalized','Position',[0.02 0.77 0.1 0.05],... 'String','<--','FontWeight','bold','FontUnits','points','FontSize',12,'Callback',@backward,'TooltipString','move backward'); - home_button=uicontrol('Style','pushbutton','Units','normalized', 'Position', [0.35 0.005 0.3 0.07],... - 'String','Work dir','FontWeight','bold','FontUnits','points','FontSize',12,'Callback',@home_dir,'TooltipString','reach the current Matlab working directory'); - - uicontrol('Style','popupmenu','Units','normalized', 'Position', [0.75 0.75 0.23 0.04],'tag','sort_option','Callback',@refresh_GUI,'Visible','off',... + home_button=uicontrol('Style','pushbutton','Units','normalized', 'Position', [0.14 0.77 0.2 0.05],... + 'String','Work dir','FontWeight','bold','FontUnits','points','FontSize',12,'Callback',@home_dir,'TooltipString','reach the current Matlab working directory'); + uicontrol('Style','pushbutton','Tag','refresh','Units','normalized','Position', [0.36 0.77 0.2 0.05],'Callback',@refresh_GUI,... + 'String','Refresh','FontWeight','bold','FontUnits','points','FontSize',12); + uicontrol('Style','popupmenu','Units','normalized', 'Position', [0.75 0.74 0.23 0.05],'tag','sort_option','Callback',@refresh_GUI,'Visible','off',... 'String',{'sort name';'sort date'},'FontUnits','points','FontSize',12,'FontWeight','bold','TooltipString','''sort_option'': sort the files by names or dates'); 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',... 'FontUnits','points','FontSize',12,'TooltipString','''list'':current list of directories, marked by +/, and files'); - uicontrol('Style','pushbutton','Units','normalized', 'Position', [0.78 0.005 0.2 0.07],'Callback',@(src,event)close(option,src,event),... - 'String','Close','FontWeight','bold','FontUnits','points','FontSize',12); - uicontrol('Style','pushbutton','Tag','refresh','Units','normalized','Position', [0.02 0.005 0.2 0.07],'Callback',@refresh_GUI,... - 'String','Refresh','FontWeight','bold','FontUnits','points','FontSize',12); - OK_button=uicontrol('Style','pushbutton','Units','normalized', 'Position', [0.75 0.81 0.23 0.07],'BackgroundColor',[0 1 0],... + + OK_button=uicontrol('Style','pushbutton','Units','normalized', 'Position', [0.58 0.005 0.2 0.07],'BackgroundColor',[0 1 0],... 'String','OK','FontWeight','bold','FontUnits','points','FontSize',12,'Callback',@(src,event)OK_Callback(option,FilterExt,src,event)); + close_button=uicontrol('Style','pushbutton','Units','normalized', 'Position', [0.78 0.005 0.2 0.07],'Callback',@(src,event)close(option,src,event),... + 'FontWeight','bold','FontUnits','points','FontSize',12); %set(hrefresh,'UserData',StatusData) if strcmp(option,'status_display') %put a run advancement display @@ -84,7 +84,10 @@ % 'String','sort: ','FontUnits','points','FontSize',12,'FontWeight','bold','HorizontalAlignment','right'); delete(home_button) - set(OK_button,'String','OPEN') + set(OK_button,'String','Open') + set(close_button,'String','Close') elseif strcmp(FilterExt,'uigetdir') %pick a directory set(path_title,'String',title); %show the input title for path (directory) + set(OK_button,'String','Select') + set(close_button,'String','Cancel') else %put a title and additional pushbuttons uicontrol('Style','text','Units','normalized', 'Position', [0.02 0.74 0.6 0.03],'BackgroundColor',BackgroundColor,... @@ -94,5 +97,6 @@ % uicontrol('Style','text','Units','normalized', 'Position', [0.37 0.8 0.35 0.03],'BackgroundColor',BackgroundColor,... % 'String','sort: ','FontUnits','points','FontSize',12,'FontWeight','bold','HorizontalAlignment','right'); - + set(OK_button,'String','Open') + set(close_button,'String','Cancel') end drawnow