Index: /trunk/src/proj_field.m =================================================================== --- /trunk/src/proj_field.m (revision 76) +++ /trunk/src/proj_field.m (revision 77) @@ -121,17 +121,17 @@ % FieldData.VarAttribute={}; % end - -if isequal(ObjectData.Style,'points') +switch ObjectData.Style + case 'points' [ProjData,errormsg]=proj_points(FieldData,ObjectData); -elseif isequal(ObjectData.Style,'line')|isequal(ObjectData.Style,'polyline') + case {'line','polyline'} [ProjData,errormsg] = proj_line(FieldData,ObjectData); -elseif isequal(ObjectData.Style,'polygon')|isequal(ObjectData.Style,'rectangle')|isequal(ObjectData.Style,'ellipse') - if isequal(ObjectData.ProjMode,'inside')|isequal(ObjectData.ProjMode,'outside') - [ProjData,errormsg] = proj_patch(FieldData,ObjectData); - else - [ProjData,errormsg] = proj_line(FieldData,ObjectData); - end + case {'polygon','rectangle','ellipse'} + if isequal(ObjectData.ProjMode,'inside')||isequal(ObjectData.ProjMode,'outside') + [ProjData,errormsg] = proj_patch(FieldData,ObjectData); + else + [ProjData,errormsg] = proj_line(FieldData,ObjectData); + end %A FAIRE : GERER MASK -elseif isequal(ObjectData.Style,'plane') + case 'plane' % if isfield(FieldData,'NbDim') & isequal(FieldData.NbDim,3) % ProjData= proj_plane3D(FieldData,ObjectData);% @@ -139,4 +139,6 @@ [ProjData,errormsg] = proj_plane(FieldData,ObjectData); % end + case 'volume' + [ProjData,errormsg] = proj_volume(FieldData,ObjectData); end if exist('IndexObj','var') @@ -154,11 +156,11 @@ width=ObjectData.Range(1,2); end -if isfield(ObjectData,'RangeX')&~isempty(ObjectData.RangeX) +if isfield(ObjectData,'RangeX')&&~isempty(ObjectData.RangeX) width=max(ObjectData.RangeX); end -if isfield(ObjectData,'RangeY')&~isempty(ObjectData.RangeY) +if isfield(ObjectData,'RangeY')&&~isempty(ObjectData.RangeY) width=max(width,max(ObjectData.RangeY)); end -if isfield(ObjectData,'RangeZ')&~isempty(ObjectData.RangeZ) +if isfield(ObjectData,'RangeZ')&&~isempty(ObjectData.RangeZ) width=max(width,max(ObjectData.RangeZ)); end @@ -951,5 +953,5 @@ test3D=0; if isfield(FieldData,'nb_dim') - test3D=isequal(FieldData.nb_dim,3) + test3D=isequal(FieldData.nb_dim,3); end test3C=test3D; %default 3 vel components @@ -1046,10 +1048,10 @@ %case of input fields with unstructured coordinates if testX - XName=FieldData.ListVarName{ivar_X} - YName=FieldData.ListVarName{ivar_Y} + XName=FieldData.ListVarName{ivar_X}; + YName=FieldData.ListVarName{ivar_Y}; eval(['coord_x=FieldData.' XName ';']) eval(['coord_y=FieldData.' YName ';']) if length(ivar_Z)==1 - ZName=FieldData.ListVarName{ivar_Z} + ZName=FieldData.ListVarName{ivar_Z}; eval(['coord_z=FieldData.' ZName ';']) end @@ -1068,5 +1070,6 @@ indcut=find(abs(fieldZ) <= width); for ivar=VarIndex - VarName=FieldData.ListVarName{ivar}; + VarName=FieldData.ListVarName{ivar}; +% eval(['size(FieldData.' VarName ')']) eval(['FieldData.' VarName '=FieldData.' VarName '(indcut);']) % end @@ -1130,6 +1133,4 @@ % different cases of projection if isequal(ObjectData.ProjMode,'projection') - %ProjData.ListDimName=[ProjData.ListDimName - %FieldData.ListDimName(DimIndices(1))];%add the point index to %the list of dimension ProjData.ListDimName=[ProjData.ListDimName FieldData.VarDimName(VarIndex(1))];%add the point index to the list of dimensions @@ -1158,9 +1159,7 @@ coord_y_proj=[YMin:DY:YMax]; DimCell={'coord_y','coord_x'}; - %ProjData.DimValue=[length(coord_y_proj) length(coord_x_proj)]; ProjData.ListVarName={'coord_y','coord_x'}; ProjData.VarDimName={'coord_y','coord_x'}; - nbcoord=2; - %ProjData.VarDimIndex={}; + nbcoord=2; ProjData.coord_y=[YMin YMax]; ProjData.coord_x=[XMin XMax]; @@ -1187,10 +1186,7 @@ ProjData.ListVarName=[ProjData.ListVarName {VarName}]; ProjData.VarDimName=[ProjData.VarDimName {DimCell}]; - %ProjData.VarDimIndex=[ProjData.VarDimIndex {[1 2]}]; if isfield(FieldData,'VarAttribute') && length(FieldData.VarAttribute) >=ivar ProjData.VarAttribute{ivar_new+nbcoord}=FieldData.VarAttribute{ivar}; end -% ProjData.VarAttribute{ivar_new}.Coord_2=[XMin XMax]; -% ProjData.VarAttribute{ivar_new}.Coord_1=[YMin YMax]; if ~isequal(ivar_FF,0) eval(['FieldData.' VarName '=FieldData.' VarName '(indsel);']) @@ -1204,6 +1200,4 @@ eval(['ProjData.' VarName '=reshape(varline,length(coord_y_proj),length(coord_x_proj));']) FF(indnan)=ones(size(indnan)); -% eval(['ProjData.' VarName '(indnan)=zeros(size(indnan));'])%put NaN to 0 -% FF=FF|FFlag; testFF=1; end @@ -1222,5 +1216,4 @@ ProjData.FF=reshape(FF,length(coord_y_proj),length(coord_x_proj)); ProjData.ListVarName=[ProjData.ListVarName {'FF'}]; - % ProjData.VarDimIndex=[ProjData.VarDimIndex {[1 2]}]; ProjData.VarDimName=[ProjData.VarDimName {DimCell}]; ProjData.VarAttribute{ivar_new+1+nbcoord}.Role='errorflag'; @@ -1235,6 +1228,4 @@ VarName=FieldData.ListVarName{VarIndex(1)}; eval(['DimValue=size(FieldData.' VarName ');']) - - %ListDimName=FieldData.ListDimName(DimIndices); ListDimName=FieldData.VarDimName{VarIndex(1)}; ProjData.ListVarName=[{AYName} {AXName} ProjData.ListVarName]; %TODO: check if it already exists in Projdata (several cells) @@ -1394,12 +1385,4 @@ for ivar=VarIndex VarName=FieldData.ListVarName{ivar}; -% if isequal(ObjectData.ProjMode,'interp')||isequal(ObjectData.ProjMode,'filter')% coordinates common to all fields -% ProjData.ListDimName={'coord_y','coord_x'}; -% ProjData.DimValue=[length(coord_y_proj) length(coord_x_proj)]; -% else -% icoord=icoord+1; -% ProjData.ListDimName=[ProjData.ListDimName {['coord_y_' num2str(icoord)],['coord_x_' num2str(icoord)]}]; -% ProjData.DimValue=[ProjData.DimValue length(coord_y_proj) length(coord_x_proj)]; -% end ProjData.ListVarName=[ProjData.ListVarName VarName]; ProjData.VarDimIndex=[ProjData.VarDimIndex [nb_dim-1 nb_dim]]; @@ -1407,6 +1390,4 @@ ProjData.VarAttribute{length(ProjData.ListVarName)}=FieldData.VarAttribute{ivar}; end -% ProjData.VarAttribute{length(ProjData.ListVarName)}.Coord_2=Xbound; -% ProjData.VarAttribute{length(ProjData.ListVarName)}.Coord_1=Ybound; eval(['ProjData.' VarName '=FieldData.' VarName '(min_ind1:max_ind1,min_ind2:max_ind2) ;']); end @@ -1452,34 +1433,12 @@ vec_B(ind_out,icolor)=zeros(size(ind_out)); end - % TODO: A REVOIR: -% if isequal(ObjectData.ProjMode,'interp') || isequal(ObjectData.ProjMode,'filter')% coordinates common to all fields -% ProjData.ListDimName={'coord_y','coord_x'}; -% ProjData.DimValue=[length(coord_y_proj) length(coord_x_proj)]; -% else -% icoord=icoord+1; -% ProjData.ListDimName=[ProjData.ListDimName {['coord_y_' num2str(icoord)],['coord_x_' num2str(icoord)]}]; -% ProjData.DimValue=[ProjData.DimValue length(coord_y_proj) length(coord_x_proj)]; -% end ProjData.ListVarName=[ProjData.ListVarName VarName]; - % ProjData.VarDimIndex=[ProjData.VarDimIndex [nb_dim-1 nb_dim]]; if length(FieldData.VarAttribute)>=ivar ProjData.VarAttribute{length(ProjData.ListVarName)+nbcoord}=FieldData.VarAttribute{ivar}; - end - % A REVOIR: -% if test_direct(2)==1 -% ProjData.VarAttribute{length(ProjData.ListVarName)}.Coord_2=[XMin XMax]; -% else -% ProjData.VarAttribute{length(ProjData.ListVarName)}.Coord_2=[XMax XMin]; -% end -% if test_direct(1)==1 -% ProjData.VarAttribute{length(ProjData.ListVarName)}.Coord_1=[YMin YMax]; -% else -% ProjData.VarAttribute{length(ProjData.ListVarName)}.Coord_1=[YMax YMin]; -% end + 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.VarDimIndex=[ProjData.VarDimIndex [nb_dim-1 nb_dim]]; ProjData.VarAttribute{length(ProjData.ListVarName)}.Role='errorflag'; else %3D case @@ -1494,17 +1453,5 @@ VarName=FieldData.ListVarName{ivar}; ProjData.ListVarName=[ProjData.ListVarName VarName]; - % ProjData.VarDimIndex=[ProjData.VarDimIndex [nb_dim-2 nb_dim-1]]; - ProjData.VarAttribute{length(ProjData.ListVarName)}=FieldData.VarAttribute{ivar}; %reproduce the variable attributes - % A REVOIR -% if test_direct_x==1 -% ProjData.VarAttribute{length(ProjData.ListVarName)}.Coord_2=[XMin XMax]; -% else -% ProjData.VarAttribute{length(ProjData.ListVarName)}.Coord_2=[XMax XMin]; -% end -% if test_direct_y==1 -% ProjData.VarAttribute{length(ProjData.ListVarName)}.Coord_1=[YMin YMax]; -% else -% ProjData.VarAttribute{length(ProjData.ListVarName)}.Coord_1=[YMax YMin]; -% end + ProjData.VarAttribute{length(ProjData.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 @@ -1545,4 +1492,593 @@ %----------------------------------------------------------------- +%project in a volume +% A FAIRE ....(copie de proj_plane) + function [ProjData,errormsg] = proj_volume(FieldData, ObjectData) +%----------------------------------------------------------------- + +%initialisation of the input parameters of the projection plane +%----------------------------------------------------------------- +ProjMode='projection';%direct projection by default +if isfield(ObjectData,'ProjMode'),ProjMode=ObjectData.ProjMode; end; + +%axis origin +if isempty(ObjectData.Coord) + ObjectData.Coord(1,1)=0;%origin of the volume set to [0 0] by default + ObjectData.Coord(1,2)=0; + ObjectData.Coord(1,3)=0; +end + +%rotation angles +Phi=0;%default +Theta=0; +Psi=0; +if isfield(ObjectData,'Phi')&& ~isempty(ObjectData.Phi) + Phi=(pi/180)*ObjectData.Phi;%first Euler angle in radian +end +if isfield(ObjectData,'Theta')&& ~isempty(ObjectData.Theta) + Theta=(pi/180)*ObjectData.Theta;%second Euler angle in radian +end +if isfield(ObjectData,'Psi')&& ~isempty(ObjectData.Psi) + Psi=(pi/180)*ObjectData.Psi;%third Euler angle in radian +end + +%components of the unity vector normal to the projection plane +NormVec_X=-sin(Phi)*sin(Theta); +NormVec_Y=cos(Phi)*sin(Theta); +NormVec_Z=cos(Theta); + +% test for 3D fields +test3D=0; +if isfield(FieldData,'nb_dim') + test3D=isequal(FieldData.nb_dim,3); +end +test3C=test3D; %default 3 vel components + +%mesh sizes DX and DY +DX=0; +DY=0; %default +if isfield(ObjectData,'DX')&~isempty(ObjectData.DX) + DX=abs(ObjectData.DX);%mesh of interpolation points +end +if isfield(ObjectData,'DY')&~isempty(ObjectData.DY) + DY=abs(ObjectData.DY);%mesh of interpolation points +end +if ~isequal(ProjMode,'projection') & (DX==0|DY==0) + errormsg='DX or DY missing'; + display(errormsg) + return +end +if isfield(ObjectData,'DZ')&~isempty(ObjectData.DZ) + DZ=abs(ObjectData.DZ);%mesh of interpolation points +end + +%extrema along each axis +testXMin=0; +testXMax=0; +testYMin=0; +testYMax=0; +if isfield(ObjectData,'RangeX') + XMin=min(ObjectData.RangeX); + XMax=max(ObjectData.RangeX); + testXMin=XMax>XMin; + testXMax=1; +end +if isfield(ObjectData,'RangeY') + YMin=min(ObjectData.RangeY); + YMax=max(ObjectData.RangeY); + testYMin=YMax>YMin; + testYMax=1; +end +width=0;%default width of the projection band +if isfield(ObjectData,'RangeZ') + ZMin=min(ObjectData.RangeZ); + ZMax=max(ObjectData.RangeZ); + testZMin=YMax>YMin; + testZMax=1; +% width=max(ObjectData.RangeZ); +end + +% initiate Matlab structure for physical field +ProjData=proj_heading(FieldData,ObjectData); +ProjData.NbDim=3; +ProjData.ListDimName={};%name of dimension +ProjData.DimValue=[];%values of dimension (nbre of vectors) +ProjData.ListVarName={}; +ProjData.VarDimName={}; + +error=0;%default +flux=0; +testfalse=0; +ListIndex={}; + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%group the variables (fields of 'FieldData') in cells of variables with the same dimensions +%----------------------------------------------------------------- +idimvar=0; +[CellVarIndex,NbDim,VarTypeCell,errormsg]=find_field_indices(FieldData); +if ~isempty(errormsg) + errormsg=['error in proj_field/proj_plane:' errormsg]; + return +end +%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 +icoord=0; +nbcoord=0;%number of added coordinate variables brought by projection +for icell=1:length(CellVarIndex) + if NbDim(icell)<2 + continue + end + VarIndex=CellVarIndex{icell};% indices of the selected variables in the list FieldData.ListVarName + VarType=VarTypeCell{icell}; + ivar_X=VarType.coord_x; + ivar_Y=VarType.coord_y; + ivar_Z=VarType.coord_z; + ivar_U=VarType.vector_x; + ivar_V=VarType.vector_y; + ivar_W=VarType.vector_z; + ivar_C=VarType.scalar ; + ivar_Anc=VarType.ancillary; + test_anc=zeros(size(VarIndex)); + test_anc(ivar_Anc)=ones(size(ivar_Anc)); + ivar_F=VarType.warnflag; + ivar_FF=VarType.errorflag; + testX=~isempty(ivar_X) && ~isempty(ivar_Y); + DimCell=FieldData.VarDimName{VarIndex(1)}; + if ischar(DimCell) + DimCell={DimCell};%name of dimensions + end + +%case of input fields with unstructured coordinates + if testX + XName=FieldData.ListVarName{ivar_X}; + YName=FieldData.ListVarName{ivar_Y}; + eval(['coord_x=FieldData.' XName ';']) + eval(['coord_y=FieldData.' YName ';']) + if length(ivar_Z)==1 + ZName=FieldData.ListVarName{ivar_Z}; + eval(['coord_z=FieldData.' ZName ';']) + end + + % translate initial coordinates + coord_x=coord_x-ObjectData.Coord(1,1); + coord_y=coord_y-ObjectData.Coord(1,2); + if ~isempty(ivar_Z) + coord_z=coord_z-ObjectData.Coord(1,3); + end + + % selection of the vectors in the projection range (3D case) + if length(ivar_Z)==1 && width > 0 + %components of the unitiy vector normal to the projection plane + fieldZ=NormVec_X*coord_x + NormVec_Y*coord_y+ NormVec_Z*coord_z;% distance to the plane + indcut=find(abs(fieldZ) <= width); + for ivar=VarIndex + VarName=FieldData.ListVarName{ivar}; + eval(['FieldData.' VarName '=FieldData.' VarName '(indcut);']) +% end + % A VOIR : CAS DE VAR STRUCTUREE MAIS PAS GRILLE REGULIERE : INTERPOLER SUR GRILLE REGULIERE + end + coord_x=coord_x(indcut); + coord_y=coord_y(indcut); + coord_z=coord_z(indcut); + end + + %rotate coordinates if needed + if isequal(Phi,0) + coord_X=coord_x; + coord_Y=coord_y; + if ~isequal(Theta,0) + coord_Y=coord_Y *cos(Theta); + end + else + coord_X=(coord_x *cos(Phi) + coord_y* sin(Phi)); + coord_Y=(-coord_x *sin(Phi) + coord_y *cos(Phi))*cos(Theta); + end + if ~isempty(ivar_Z) + coord_Y=coord_Y+coord_z *sin(Theta); + end + if ~isequal(Psi,0) + coord_X=(coord_X *cos(Psi) - coord_Y* sin(Psi));%A VERIFIER + coord_Y=(coord_X *sin(Psi) + coord_Y* cos(Psi)); + end + coord_Z=coord_z; %TO UPDATE + %restriction to the range of x and y if imposed + testin=ones(size(coord_X)); %default + testbound=0; + if testXMin + testin=testin & (coord_X >= XMin); + testbound=1; + end + if testXMax + testin=testin & (coord_X <= XMax); + testbound=1; + end + if testYMin + testin=testin & (coord_Y >= YMin); + testbound=1; + end + if testYMin + testin=testin & (coord_Y <= YMax); + testbound=1; + end + if testbound + indcut=find(testin); + for ivar=VarIndex + VarName=FieldData.ListVarName{ivar}; + eval(['FieldData.' VarName '=FieldData.' VarName '(indcut);']) + end + coord_X=coord_X(indcut); + coord_Y=coord_Y(indcut); + if length(ivar_Z)==1 + coord_Z=coord_Z(indcut); + end + end + % different cases of projection + if isequal(ObjectData.ProjMode,'projection') + %the list of dimension + ProjData.ListDimName=[ProjData.ListDimName FieldData.VarDimName(VarIndex(1))];%add the point index to the list of dimensions + ProjData.DimValue=[ProjData.DimValue length(coord_X)]; + nbvar=0; + for ivar=VarIndex %transfer variables to the projection plane + VarName=FieldData.ListVarName{ivar}; + if ivar==ivar_X %x coordinate + eval(['ProjData.' VarName '=coord_X;']) + elseif ivar==ivar_Y % y coordinate + eval(['ProjData.' VarName '=coord_Y;']) + elseif isempty(ivar_Z) || ivar~=ivar_Z % other variables (except Z coordinate wyhich is not reproduced) + eval(['ProjData.' VarName '=FieldData.' VarName ';']) + end + if isempty(ivar_Z) || ivar~=ivar_Z + ProjData.ListVarName=[ProjData.ListVarName VarName]; + ProjData.VarDimName=[ProjData.VarDimName DimCell]; + nbvar=nbvar+1; + if isfield(FieldData,'VarAttribute') & length(FieldData.VarAttribute) >=ivar + ProjData.VarAttribute{nbvar}=FieldData.VarAttribute{ivar}; + end + end + end + elseif isequal(ObjectData.ProjMode,'interp')||isequal(ObjectData.ProjMode,'filter')%interpolate data on a regular grid + coord_x_proj=[XMin:DX:XMax]; + coord_y_proj=[YMin:DY:YMax]; + coord_z_proj=[ZMin:DZ:ZMax]; + [XI,YI,ZI]=meshgrid(coord_x_proj,coord_y_proj,coord_z_proj); + DimCell={'coord_y','coord_x'}; + ProjData.ListVarName={'coord_z','coord_y','coord_x'}; + ProjData.VarDimName={'coord_z','coord_y','coord_x'}; + nbcoord=3; + ProjData.coord_z=[ZMin ZMax]; + ProjData.coord_y=[YMin YMax]; + ProjData.coord_x=[XMin XMax]; + if isempty(ivar_X), ivar_X=0; end; + if isempty(ivar_Y), ivar_Y=0; end; + if isempty(ivar_Z), ivar_Z=0; end; + if isempty(ivar_U), ivar_U=0; end; + if isempty(ivar_V), ivar_V=0; end; + if isempty(ivar_W), ivar_W=0; end; + if isempty(ivar_F), ivar_F=0; end; + if isempty(ivar_FF), ivar_FF=0; end; + if ~isequal(ivar_FF,0) + VarName_FF=FieldData.ListVarName{ivar_FF}; + eval(['indsel=find(FieldData.' VarName_FF '==0);']) + coord_X=coord_X(indsel); + coord_Y=coord_Y(indsel); + end + FF=zeros(1,length(coord_y_proj)*length(coord_x_proj)); + testFF=0; + size(XI) + size(YI) + size(ZI) + for ivar=VarIndex + VarName=FieldData.ListVarName{ivar}; + if ~( ivar==ivar_X || ivar==ivar_Y || ivar==ivar_Z || ivar==ivar_F || ivar==ivar_FF || test_anc(ivar)==1) + ivar_new=ivar_new+1; + ProjData.ListVarName=[ProjData.ListVarName {VarName}]; + ProjData.VarDimName=[ProjData.VarDimName {DimCell}]; + if isfield(FieldData,'VarAttribute') && length(FieldData.VarAttribute) >=ivar + ProjData.VarAttribute{ivar_new+nbcoord}=FieldData.VarAttribute{ivar}; + end + if ~isequal(ivar_FF,0) + eval(['FieldData.' VarName '=FieldData.' VarName '(indsel);']) + end + eval(['ProjData.' VarName '=griddata3(double(coord_X),double(coord_Y),double(coord_Z),double(FieldData.' VarName '),XI,YI,ZI);']) +% eval(['varline=reshape(ProjData.' VarName ',1,length(coord_y_proj)*length(coord_x_proj));']) +% FFlag= isnan(varline); %detect undefined values NaN +% indnan=find(FFlag); +% if~isempty(indnan) +% varline(indnan)=zeros(size(indnan)); +% eval(['ProjData.' VarName '=reshape(varline,length(coord_y_proj),length(coord_x_proj));']) +% FF(indnan)=ones(size(indnan)); +% testFF=1; +% end +% if ivar==ivar_U +% ivar_U=ivar_new; +% end +% if ivar==ivar_V +% ivar_V=ivar_new; +% end +% if ivar==ivar_W +% ivar_W=ivar_new; +% end + end + end + if testFF + ProjData.FF=reshape(FF,length(coord_y_proj),length(coord_x_proj)); + ProjData.ListVarName=[ProjData.ListVarName {'FF'}]; + ProjData.VarDimName=[ProjData.VarDimName {DimCell}]; + ProjData.VarAttribute{ivar_new+1+nbcoord}.Role='errorflag'; + end + end +%case of fields defined on a structured grid + else + AYName=FieldData.ListVarName{VarType.coord(1)}; + AXName=FieldData.ListVarName{VarType.coord(2)}; + eval(['AX=FieldData.' AXName ';']) + eval(['AY=FieldData.' AYName ';']) + VarName=FieldData.ListVarName{VarIndex(1)}; + eval(['DimValue=size(FieldData.' VarName ');']) + ListDimName=FieldData.VarDimName{VarIndex(1)}; + ProjData.ListVarName=[{AYName} {AXName} ProjData.ListVarName]; %TODO: check if it already exists in Projdata (several cells) + ProjData.VarDimName=[{AYName} {AXName} ProjData.VarDimName]; + nbcolor=1; %default + for idim=1:length(ListDimName) + DimName=ListDimName{idim}; + if isequal(DimName,'rgb')|isequal(DimName,'nb_coord')|isequal(DimName,'nb_coord_i') + nbcolor=DimValue(idim); + DimIndices(idim)=[]; + DimValue(idim)=[]; + end + if isequal(DimName,'nb_coord_j')% NOTE: CASE OF TENSOR NOT TREATED + DimIndices(idim)=[]; + DimValue(idim)=[]; + end + end + ind_1=find(DimValue==1); + DimIndices(ind_1)=[]; %suppress singleton dimensions +% indxy=find(DimVarIndex(DimIndices));%select dimension variables (DimIndices non zero) + nb_dim=length(DimIndices);%number of space dimensions + Coord_z=[]; + Coord_y=[]; + Coord_x=[]; + + for idim=1:nb_dim %loop on space dimensions + test_interp(idim)=0;%test for coordiate interpolation (non regular grid), =0 by default + test_coord(idim)=0;%test for defined coordinates, =0 by default + ivar=DimVarIndex(DimIndices(idim));% index of the variable corresponding to the current dimension + if ~isequal(ivar,0)% a variable corresponds to the current dimension + eval(['Coord{idim}=FieldData.' FieldData.ListVarName{ivar} ';']) ;% position for the first index + DCoord=diff(Coord{idim}); + DCoord_min(idim)=min(DCoord); + DCoord_max=max(DCoord); + test_direct(idim)=DCoord_max>0;% =1 for increasing values, 0 otherwise + test_direct_min=DCoord_min(idim)>0;% =1 for increasing values, 0 otherwise + if ~isequal(test_direct(idim),test_direct_min) + msgbox_uvmat('ERROR',['non monotonic dimension variable # ' num2str(idim) ' in proj_field.m']) + return + end + test_interp(idim)=(DCoord_max-DCoord_min(idim))> 0.0001*abs(DCoord_max);% test grid regularity + test_coord(idim)=1; + + else % no variable associated with the first dimension, look for variable attributes Coord_1, _2 or _3 + Coord_i_str=['Coord_' num2str(idim)]; + DCoord_min(idim)=1;%default + Coord{idim}=[0.5 DimValue(idim)-0.5]; + test_direct(idim)=1; + end + end + if nb_dim==2 + if DY==0 + DY=abs(DCoord_min(1)); + end + npY=1+round(abs(Coord{1}(end)-Coord{1}(1))/DY);%nbre of points after interpolation + npy=1+round(abs(Coord{1}(end)-Coord{1}(1))/abs(DCoord_min(1)));%nbre of points after possible interpolation on a regular grid + if DX==0 + DX=abs(DCoord_min(2)); + end + npX=1+round(abs(Coord{2}(end)-Coord{2}(1))/DX);%nbre of points after interpol + npx=1+round(abs(Coord{2}(end)-Coord{2}(1))/abs(DCoord_min(2)));%nbre of points after possible interpolation on a regular grid + Coord_y=linspace(Coord{1}(1),Coord{1}(end),npY); + test_direct_y=test_direct(1); + Coord_x=linspace(Coord{2}(1),Coord{2}(end),npX); + test_direct_x=test_direct(2); + DAX=DCoord_min(2); + DAY=DCoord_min(1); + elseif nb_dim==3 + DZ=abs(DCoord_min(1)); + npz=1+round(abs(Coord{1}(end)-Coord{1}(1))/DZ);%nbre of points after interpolation + if DY==0 + DY=abs(DCoord_min(2)); + end + npY=1+round(abs(Coord{2}(end)-Coord{2}(1))/DY);%nbre of points after interpol + npy=1+round(abs(Coord{2}(end)-Coord{2}(1))/abs(DCoord_min(2)));%nbre of points before interpol + if DX==0 + DX=abs(DCoord_min(3)); + end + npX=1+round(abs(Coord{3}(end)-Coord{3}(1))/DX);%nbre of points after interpol + npx=1+round(abs(Coord{3}(end)-Coord{3}(1))/abs(DCoord_min(3)));%nbre of points before interpol + Coord_z=linspace(Coord{1}(1),Coord{1}(end),npz); + test_direct_z=test_direct(1); + Coord_y=linspace(Coord{2}(1),Coord{2}(end),npY); + test_direct_y=test_direct(2); + Coord_x=linspace(Coord{3}(1),Coord{3}(end),npX); + test_direct_x=test_direct(3); + end + 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); + ycorner=[maxAY maxAY minAY minAY]-ObjectData.Coord(1,2); + xcor_new=xcorner*cos(Phi)+ycorner*sin(Phi);%coord new frame + ycor_new=-xcorner*sin(Phi)+ycorner*cos(Phi); + if ~testXMax + XMax=max(xcor_new); + end + if ~testXMin + XMin=min(xcor_new); + end + if ~testYMax + YMax=max(ycor_new); + end + if ~testYMin + YMin=min(ycor_new); + end + DXinit=(maxAX-minAX)/(npx-1); + DYinit=(maxAY-minAY)/(npy-1); + if DX==0 + DX=DXinit; + end + if DY==0 + DY=DYinit; + 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 rotation and interpolation + if isequal(ProjMode,'projection') && isequal(Phi,0) && isequal(Theta,0) && isequal(Psi,0) + if test_direct(1) + min_ind1=ceil((YMin-Coord{1}(1))/DYinit)+1; + max_ind1=floor((YMax-Coord{1}(1))/DYinit)+1; + Ybound(1)=Coord{1}(1)+DYinit*(min_ind1-1); + Ybound(2)=Coord{1}(1)+DYinit*(max_ind1-1); + else + min_ind1=ceil((Coord{1}(1)-YMax)/DYinit)+1; + max_ind1=floor((Coord{1}(1)-YMin)/DYinit)+1; + Ybound(2)=Coord{1}(1)-DYinit*(max_ind1-1); + Ybound(1)=Coord{1}(1)-DYinit*(min_ind1-1); + end + if test_direct(2)==1 + min_ind2=ceil((XMin-Coord{2}(1))/DXinit)+1; + max_ind2=floor((XMax-Coord{2}(1))/DXinit)+1; + Xbound(1)=Coord{2}(1)+DXinit*(min_ind2-1); + Xbound(2)=Coord{2}(1)+DXinit*(max_ind2-1); + else + min_ind2=ceil((Coord{2}(1)-XMax)/DXinit)+1; + max_ind2=floor((Coord{2}(1)-XMin)/DXinit)+1; + Xbound(2)=Coord{2}(1)+DXinit*(max_ind2-1); + Xbound(1)=Coord{2}(1)+DXinit*(min_ind2-1); + end + min_ind1=max(min_ind1,1);% deals with margin (bound lower than the first index) + min_ind2=max(min_ind2,1); + max_ind1=min(max_ind1,npy); + max_ind2=min(max_ind2,npx); + for ivar=VarIndex + VarName=FieldData.ListVarName{ivar}; + ProjData.ListVarName=[ProjData.ListVarName VarName]; + ProjData.VarDimIndex=[ProjData.VarDimIndex [nb_dim-1 nb_dim]]; + if length(FieldData.VarAttribute)>=ivar + ProjData.VarAttribute{length(ProjData.ListVarName)}=FieldData.VarAttribute{ivar}; + end + eval(['ProjData.' VarName '=FieldData.' VarName '(min_ind1:max_ind1,min_ind2:max_ind2) ;']); + end + else + % case with rotation and/or interpolation + if isempty(Coord_z) %2D case + [X,Y]=meshgrid(coord_x_proj,coord_y_proj);%grid in the new coordinates + XIMA=ObjectData.Coord(1,1)+(X)*cos(Phi)-Y*sin(Phi);%corresponding coordinates in the original image + YIMA=ObjectData.Coord(1,2)+(X)*sin(Phi)+Y*cos(Phi); + 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<=npx & YIMA >=1 & YIMA<=npy;%flagin=1 inside the original image + if isequal(ObjectData.ProjMode,'filter') + npx_filter=ceil(abs(DX/DAX)); + npy_filter=ceil(abs(DY/DAY)); + Mfilter=ones(npy_filter,npx_filter)/(npx_filter*npy_filter); + test_filter=1; + else + test_filter=0; + end + for ivar=VarIndex + VarName=FieldData.ListVarName{ivar} ; + if test_interp(1) | test_interp(2)%interpolate on a regular grid + eval(['FieldData.' VarName '=interp2(Coord{2},Coord{1},FieldData.' VarName ',Coord_x,Coord_y'');']) %TO TEST + end + %filter the field (image) if option 'filter' is used + if test_filter + Aclass=class(FieldData.A); + eval(['FieldData.' VarName '=filter2(Mfilter,FieldData.' VarName ',''valid'');']) + if ~isequal(Aclass,'double') + eval(['FieldData.' VarName '=' Aclass '(FieldData.' VarName ');'])%revert to integer values + end + end + eval(['vec_A=reshape(FieldData.' VarName ',npx*npy,nbcolor);'])%put the original image in line + ind_in=find(flagin); + ind_out=find(~flagin); + ICOMB=(XIMA-1)*npy+YIMA; + ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A + vec_B(ind_in,[1:nbcolor])=vec_A(ICOMB,:); + for icolor=1:nbcolor + vec_B(ind_out,icolor)=zeros(size(ind_out)); + end + ProjData.ListVarName=[ProjData.ListVarName VarName]; + if 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.VarAttribute{length(ProjData.ListVarName)}.Role='errorflag'; + else %3D case + if isequal(Theta,0) & isequal(Phi,0) + test_sup=(Coord{1}>=ObjectData.Coord(1,3)); + iz_sup=find(test_sup); + iz=iz_sup(1); + if iz>=1 & iz<=npz + ProjData.ListDimName=[ProjData.ListDimName ListDimName(2:end)]; + ProjData.DimValue=[ProjData.DimValue npY npX]; + for ivar=VarIndex + VarName=FieldData.ListVarName{ivar}; + ProjData.ListVarName=[ProjData.ListVarName VarName]; + ProjData.VarAttribute{length(ProjData.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 + if test_interp(2) | test_interp(3) + eval(['ProjData.' VarName '=interp2(Coord{3},Coord{2},ProjData.' VarName ',Coord_x,Coord_y'');']) + end + end + end + else + errormsg='projection of structured coordinates on oblique plane not yet implemented'; + %TODO: use interp3 + return + end + end + end + end + %projection of velocity components in the rotated coordinates + if ~isequal(Phi,0) && 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(Phi)*ProjData.' UName '+ sin(Phi)*ProjData.' VName ';']) + eval(['ProjData.' VName '=cos(Theta)*(-sin(Phi)*ProjData.' UName '+ cos(Phi)*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 + +%----------------------------------------------------------------- %transmit the global attributes function [ProjData,errormsg]=proj_heading(FieldData,ObjectData) Index: /trunk/src/read_set_object.m =================================================================== --- /trunk/src/read_set_object.m (revision 76) +++ /trunk/src/read_set_object.m (revision 77) @@ -67,5 +67,5 @@ end else - data.RangeZ=ZMax; + data.RangeZ(2)=ZMax; dimrange=[dimrange(1) 3]; end Index: /trunk/src/set_object.m =================================================================== --- /trunk/src/set_object.m (revision 76) +++ /trunk/src/set_object.m (revision 77) @@ -238,9 +238,4 @@ end end -% if desable_open -% set(handles.OPEN,'Visible','off') -% else -% set(handles.OPEN,'Visible','on') -% end if enable_plot set(handles.PLOT,'enable','on') @@ -331,5 +326,5 @@ menu_proj={'inside';'outside';'mask_inside';'mask_outside'}; case 'volume' - menu_proj={'none'}; + menu_proj={'interp';'none'}; end proj_index=get(handles.ProjMode,'Value'); Index: /trunk/src/uvmat.m =================================================================== --- /trunk/src/uvmat.m (revision 76) +++ /trunk/src/uvmat.m (revision 77) @@ -4324,9 +4324,10 @@ function MenuExportField_Callback(hObject, eventdata, handles) global CurData -huvmat=findobj(allchild(0),'Name','uvmat'); -UvData=get(huvmat,'UserData'); -CurData=UvData.ProjField_1; -evalin('base','global CurData')%make CurData global in the workspace -display(['UserData of view_field :']) +CurData=get(handles.uvmat,'UserData'); +% if isfield(UvData,'ProjField') + CurData=UvData; + evalin('base','global CurData')%make CurData global in the workspace + display(['current field :']) +% end evalin('base','CurData') %display CurData in the workspace commandwindow; @@ -4722,7 +4723,10 @@ % ------------------------------------------------------------------ function Menuvolume_Callback(hObject, eventdata, handles) +data.Style='volume'; +data.ProjMode='interp';%default % set(handles.create,'Visible','on') % set(handles.create,'Value',1) -VOLUME_Callback(hObject,eventdata,handles) +% VOLUME_Callback(hObject,eventdata,handles) +create_object(data,handles) %------------------------------------------------------------------------