Index: /trunk/src/proj_field.m =================================================================== --- /trunk/src/proj_field.m (revision 1059) +++ /trunk/src/proj_field.m (revision 1060) @@ -969,5 +969,5 @@ %----------------------------------------------------------------- -%project on a plane +%project on a plane % AJOUTER flux,circul,error function [ProjData,errormsg] = proj_plane(FieldData, ObjectData) @@ -988,7 +988,11 @@ % ObjectData.Angle(1)=90*Delta_x/Delta_mod; % ObjectData.0(2)=90*Delta_y/Delta_mod; -% end -if isfield(ObjectData,'Angle')&& isequal(size(ObjectData.Angle),[1 2])&& ~isequal(ObjectData.Angle,[0 0]) - test90y=0;%isequal(ObjectData.Angle,[0 90 0]); +% end +if isfield(ObjectData,'Angle') + checkM2=0; + if numel(ObjectData.Angle)==2 && ~isequal(ObjectData.Angle,[0 0]) + checkM2=1; + test90y=0;%isequal(ObjectData.Angle,[0 90 0]); + end PlaneAngle=(pi/180)*ObjectData.Angle; % om=norm(PlaneAngle);%norm of rotation angle in radians @@ -1003,6 +1007,9 @@ M1=[cos(PlaneAngle(1)) sin(PlaneAngle(1)) 0;-sin(PlaneAngle(1)) cos(PlaneAngle(1)) 0;0 0 1]; - M2=[1 0 0;0 cos(PlaneAngle(2)) sin(PlaneAngle(2));0 -sin(PlaneAngle(2)) cos(PlaneAngle(2))]; - M=M2*M1;% first rotate in the x,y plane with angle PlaneAngle(1), then slant around the new x axis0 with angle PlaneAngle(2) + M=M1; + if checkM2 + M2=[1 0 0;0 cos(PlaneAngle(2)) sin(PlaneAngle(2));0 -sin(PlaneAngle(2)) cos(PlaneAngle(2))]; + M=M2*M1;% first rotate in the x,y plane with angle PlaneAngle(1), then slant around the new x axis0 with angle PlaneAngle(2) + end norm_plane=M*[0 0 1]'; @@ -1029,5 +1036,5 @@ InterpMesh=min(DX,DY);%mesh used for interpolation in a slanted plane % if strcmp(ObjectData.Type,'plane_z') -% InterpMesh=10*InterpMesh;%TODO: temporary, to shorten computation +% InterpMesh=10*InterpMesh;%TODO: temporary, to shorten computation % end @@ -1118,5 +1125,5 @@ if ~strcmp(ObjectData.ProjMode,'projection')&& ~strcmp(ObjectData.Type,'plane_z')% TODO:rationalize %% define the new coordinates in case of interpolation on a imposed grid - if ~testYMin + if ~testYMin errormsg='min Y value not defined for the projection grid';return end @@ -1177,6 +1184,6 @@ [XI,YI]=meshgrid(coord_x_proj,coord_y_proj);%grid in the new coordinates ProjData.VarDimName={AYName,AXName}; -% 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)); + % 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)); else% we use the existing grid from field cell #icell_grid NbDim=NbDimArray(icell_grid); @@ -1185,5 +1192,5 @@ AYDimName=FieldData.VarDimName{CellInfo{icell_grid}.CoordIndex(NbDim-1)};% AXDimName=FieldData.VarDimName{CellInfo{icell_grid}.CoordIndex(NbDim)};% - ProjData.VarDimName={AYDimName,AXDimName}; + ProjData.VarDimName={AYDimName,AXDimName}; ProjData.(AYName)=FieldData.(AYName); % new (projected ) y coordinates ProjData.(AXName)=FieldData.(AXName); % new (projected ) y coordinates @@ -1194,5 +1201,5 @@ ProjData.VarAttribute{2}.Role='coord_x'; end - + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % LOOP ON FIELD CELLS, PROJECT VARIABLES @@ -1256,13 +1263,16 @@ %rotate coordinates if needed: coord_X,coord_Y= = coordinates in the new plane - Psi=PlaneAngle(1); - Theta=PlaneAngle(2); - % Phi=PlaneAngle(3); + Phi=PlaneAngle(1); + % Theta=PlaneAngle(2); + % Phi=PlaneAngle(3); if testangle && ~test90y && ~test90x;%=1 for slanted plane - coord_X=(coord_x *cos(Phi) + coord_y* sin(Phi)); - 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_Y=(coord_X *sin(Psi) + coord_Y* cos(Psi)); + coord_X=coord_x *cos(Phi) + coord_y* sin(Phi); + coord_Y=-coord_x *sin(Phi) + coord_y *cos(Phi); + if checkM2 + coord_Y=coord_Y*cos(PlaneAngle(2)); + coord_Y=coord_Y+coord_z *sin(PlaneAngle(2)); + coord_X=(coord_X *cos(Psi) - coord_Y* sin(Psi));%A VERIFIER + coord_Y=(coord_X *sin(Psi) + coord_Y* cos(Psi)); + end else coord_X=coord_x; @@ -1377,6 +1387,6 @@ end if isfield (CellInfo{icell},'VarIndex_vector_x')&& isfield (CellInfo{icell},'VarIndex_vector_y') - vector_x_proj=CellInfo{icell}.VarIndex_vector_x; %preserve for next cell - vector_y_proj=CellInfo{icell}.VarIndex_vector_y; %preserve for next cell + vector_x_proj=CellInfo{icell}.VarIndex_vector_x; %preserve for next cell + vector_y_proj=CellInfo{icell}.VarIndex_vector_y; %preserve for next cell end end @@ -1397,16 +1407,9 @@ %rotate coordinates if needed: coord_X,coord_Y= = coordinates in the new plane - Psi=PlaneAngle(1); - Theta=PlaneAngle(2); - % Phi=PlaneAngle(3); + Phi=PlaneAngle(1); if testangle && ~test90y && ~test90x;%=1 for slanted plane - new_XI=XI*cos(Phi) - YI*sin(Phi)+ObjectData.Coord(1); + new_XI=XI *cos(Phi) - YI* sin(Phi)+ObjectData.Coord(1); YI=XI *sin(Phi) + YI *cos(Phi)+ObjectData.Coord(2); XI=new_XI; - % if checkUV - % UValue=cos(Phi)*FieldVar(:,:,1)+ sin(Phi)*FieldVar(:,:,2); - % FieldVar(:,:,2)=-sin(Phi)*FieldVar(:,:,1)+ cos(Phi)*FieldVar(:,:,2); - % FieldVar(:,:,1)=UValue; - % end end @@ -1445,5 +1448,5 @@ DimValue(DimValue==1)=[];%remove singleton dimensions NbDim=numel(DimValue);%update number of space dimensions - % nbcolor=1; %default number of 'color' components: third matrix index without corresponding coordinate + % nbcolor=1; %default number of 'color' components: third matrix index without corresponding coordinate if NbDim>=3 if NbDim>3 @@ -1492,10 +1495,10 @@ end if testYMax - YIndexMax=(YMax-Coord{NbDim-1}(1))/DY+1;% matrix index corresponding to the max y value for the new field + YIndexMax=(YMax-Coord{NbDim-1}(1))/DY+1;% matrix index corresponding to the max y value for the new field if testYMin%test_direct(indY) YIndexMin=(YMin-Coord{NbDim-1}(1))/DY+1;% matrix index corresponding to the min x value for the new field else YIndexMin=1; - end + end else YIndexMax=numel(Coord{NbDim-1}); @@ -1503,10 +1506,10 @@ end if testXMax - XIndexMax=(XMax-Coord{NbDim}(1))/DX+1;% matrix index corresponding to the max y value for the new field + XIndexMax=(XMax-Coord{NbDim}(1))/DX+1;% matrix index corresponding to the max y value for the new field if testYMin%test_direct(indY) XIndexMin=(XMin-Coord{NbDim}(1))/DX+1;% matrix index corresponding to the min x value for the new field else XIndexMin=1; - end + end else XIndexMax=numel(Coord{NbDim}); @@ -1561,13 +1564,13 @@ end if testXMax - ProjData.(AXName)=Coord{NbDim}(1)+DX*(XIndexRange-1); %record the new (projected ) x coordinates + ProjData.(AXName)=Coord{NbDim}(1)+DX*(XIndexRange-1); %record the new (projected ) x coordinates else - ProjData.(AXName)=FieldData.(AXName); + ProjData.(AXName)=FieldData.(AXName); end if testYMax ProjData.(AYName)=Coord{NbDim-1}(1)+DY*(YIndexRange-1); %record the new (projected ) x coordinates else - ProjData.(AYName)=FieldData.(AYName); - end + ProjData.(AYName)=FieldData.(AYName); + end end end @@ -1609,5 +1612,5 @@ XI=ObjectData.Coord(1,1)+(X)*cos(PlaneAngle(2))-YI*sin(PlaneAngle(1));%corresponding coordinates in the original system YI=ObjectData.Coord(1,2)+(X)*sin(PlaneAngle(2))+YI*cos(PlaneAngle(1)); - + if numel(Coord{1})==2% x coordinate defined by its bounds, get the whole set Coord{1}=linspace(Coord{1}(1),Coord{1}(2),CellInfo{icell}.CoordSize(1)); @@ -1634,5 +1637,5 @@ ProjData.(VarName)=interp2(X,Y,double(FieldData.(VarName)(:,:,1)),XI,YI,'*linear'); for icolor=2:size(FieldData.(VarName),3)% project 'color' components - ProjData.(VarName)=cat(3,ProjData.(VarName),interp2(X,Y,double(FieldData.(VarName)(:,:,icolor)),XI,YI,'*linear')); + ProjData.(VarName)=cat(3,ProjData.(VarName),interp2(X,Y,double(FieldData.(VarName)(:,:,icolor)),XI,YI,'*linear')); end end @@ -1674,5 +1677,5 @@ end end - else %projection of structured coordinates on oblique plane + else %projection of structured coordinates on oblique plane % determine the boundaries of the projected field, % first find the 8 summits of the initial volume in the @@ -1683,5 +1686,5 @@ Coord{3}=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(3)});%initial x coordinates summit=zeros(3,8);% initialize summit coordinates - summit(1,1:4)=[Coord{3}(1) Coord{3}(end) Coord{3}(1) Coord{3}(end)];%square + summit(1,1:4)=[Coord{3}(1) Coord{3}(end) Coord{3}(1) Coord{3}(end)];%square summit(2,1:4)=[Coord{2}(1) Coord{2}(1) Coord{2}(end) Coord{2}(end)];% square at z= Coord{1}(1) summit(1:2,5:8)=summit(1:2,1:4); @@ -1691,7 +1694,7 @@ ObjectData.Coord=ObjectData.Coord';% set ObjectData.Coord as a vertical vector if size(ObjectData.Coord,1)<3 - ObjectData.Coord=[ObjectData.Coord; 0];%add z origin at z=0 by default - end - + ObjectData.Coord=[ObjectData.Coord; 0];%add z origin at z=0 by default + end + M1=[cos(PlaneAngle(1)) sin(PlaneAngle(1)) 0;-sin(PlaneAngle(1)) cos(PlaneAngle(1)) 0;0 0 1]; M2=[1 0 0;0 cos(PlaneAngle(2)) sin(PlaneAngle(2));0 -sin(PlaneAngle(2)) cos(PlaneAngle(2))]; @@ -1712,16 +1715,16 @@ %modangle=sqrt(PlaneAngle(1)*PlaneAngle(1)+PlaneAngle(2)*PlaneAngle(2)); -% cosphi=PlaneAngle(1)/modangle; -% sinphi=PlaneAngle(2)/modangle; + % cosphi=PlaneAngle(1)/modangle; + % sinphi=PlaneAngle(2)/modangle; iX=[coord_x_proj(end)-coord_x_proj(1);0;0]/(npx-1); iY=[0;coord_y_proj(end)-coord_y_proj(1);0]/(npy-1); iZ=[0;0;coord_z_proj(end)-coord_z_proj(1)]/(npz-1); -% iX(1:2)=[cosphi -sinphi;sinphi cosphi]*iX(1:2); -% iY(1:2)=[-cosphi -sinphi;sinphi cosphi]*iY(1:2); + % iX(1:2)=[cosphi -sinphi;sinphi cosphi]*iX(1:2); + % iY(1:2)=[-cosphi -sinphi;sinphi cosphi]*iY(1:2); ix=M_inv*iX;% vector along the new x coordinates transformed into old coordinates iy=M_inv*iY;% vector along y coordinates iz=M_inv*iZ;% vector along z coordinates - + [Grid_x,Grid_y,Grid_z]=meshgrid(0:npx-1,0:npy-1,0:npz-1); if ismatrix(Grid_x)% add a singleton in case of a single z value @@ -1733,16 +1736,16 @@ YI=Origin(2)+ix(2)*Grid_x+iy(2)*Grid_y+iz(2)*Grid_z; ZI=Origin(3)+ix(3)*Grid_x+iy(3)*Grid_y+iz(3)*Grid_z; - [X,Y,Z]=meshgrid(Coord{3},Coord{2},Coord{1});% mesh in the initial coordinates + [X,Y,Z]=meshgrid(Coord{3},Coord{2},Coord{1});% mesh in the initial coordinates for ivar=VarIndex - VarName=FieldData.ListVarName{ivar}; - ListVarName=[ListVarName VarName]; - VarDimName=[VarDimName {{'coord_y','coord_x'}}]; - VarAttribute{length(ListVarName)}=FieldData.VarAttribute{ivar}; %reproduce the variable attributes - FieldData.(VarName)=permute(FieldData.(VarName),[2 3 1]); - ProjData.coord_x=coord_x_proj; - ProjData.coord_y=coord_y_proj; - ProjData.(VarName)=interp3(X,Y,Z,double(FieldData.(VarName)),XI,YI,ZI,'*linear'); - ProjData.(VarName)=nanmean(ProjData.(VarName),3); - ProjData.(VarName)=squeeze(ProjData.(VarName)); + VarName=FieldData.ListVarName{ivar}; + ListVarName=[ListVarName VarName]; + VarDimName=[VarDimName {{'coord_y','coord_x'}}]; + VarAttribute{length(ListVarName)}=FieldData.VarAttribute{ivar}; %reproduce the variable attributes + FieldData.(VarName)=permute(FieldData.(VarName),[2 3 1]); + ProjData.coord_x=coord_x_proj; + ProjData.coord_y=coord_y_proj; + ProjData.(VarName)=interp3(X,Y,Z,double(FieldData.(VarName)),XI,YI,ZI,'*linear'); + ProjData.(VarName)=nanmean(ProjData.(VarName),3); + ProjData.(VarName)=squeeze(ProjData.(VarName)); end end @@ -1775,6 +1778,11 @@ UName=ListVarName{ivar_U}; VName=ListVarName{ivar_V}; - UValue=cos(PlaneAngle(3))*ProjData.(UName)+ sin(PlaneAngle(3))*ProjData.(VName); - ProjData.(VName)=(-sin(PlaneAngle(3))*ProjData.(UName)+ cos(PlaneAngle(3))*ProjData.(VName)); + if checkM2 + UValue=cos(PlaneAngle(1))*ProjData.(UName)+ sin(PlaneAngle(1))*ProjData.(VName); + ProjData.(VName)=(-sin(PlaneAngle(1))*ProjData.(UName)+ cos(PlaneAngle(1))*ProjData.(VName)); + else + UValue=cos(PlaneAngle(1))*ProjData.(UName)+ sin(PlaneAngle(1))*ProjData.(VName); + ProjData.(VName)=(-sin(PlaneAngle(1))*ProjData.(UName)+ cos(PlaneAngle(1))*ProjData.(VName)); + end ProjData.(UName)=UValue; if ~isempty(ivar_W) Index: /trunk/src/series/extract_rdvision.m =================================================================== --- /trunk/src/series/extract_rdvision.m (revision 1059) +++ /trunk/src/series/extract_rdvision.m (revision 1060) @@ -538,4 +538,5 @@ end Dti_stamp=(timestamp(1+NbDti,1)-timestamp(1,1))/NbDti; +Dti_stamp=(timestamp(1+NbDti,1)-timestamp(2,1))/(NbDti-1); t=set(t,uid_content_Dti,'value',num2str(Dti_stamp));%corret Dti if abs(Dti_stamp-Dti)>Dti/1000 Index: /trunk/src/set_object.m =================================================================== --- /trunk/src/set_object.m (revision 1059) +++ /trunk/src/set_object.m (revision 1060) @@ -165,5 +165,5 @@ % if isfield(data,'Angle') && isequal(numel(data.Angle),3) set(handles.num_Angle_1,'String',num2str(data.Angle(1))) - set(handles.num_Angle_2,'String',num2str(data.Angle(2))) +% set(handles.num_Angle_2,'String',num2str(data.Angle(2))) % set(handles.num_Angle_3,'String',num2str(data.Angle(3))) % end @@ -332,5 +332,5 @@ set(handles.num_RangeY_2,'TooltipString',['num_RangeY_2: half width of the ' ObjectStyle]) case {'plane','plane_z'} - set(handles.num_Angle_3,'Visible','on') + set(handles.num_Angle_1,'Visible','on') set(handles.num_RangeX_1,'Visible','on') set(handles.num_RangeX_2,'Visible','on') @@ -418,17 +418,18 @@ PlaneAngle(2)=str2num(get(handles.num_Angle_2,'String'));%second angle in degrees %PlaneAngle(3)=str2num(get(handles.num_Angle_3,'String'));%second angle in degrees -om=norm(PlaneAngle);%norm of rotation angle in radians -OmAxis=PlaneAngle/om; %unit vector marking the rotation axis -cos_om=cos(pi*om/180); -sin_om=sin(pi*om/180); -coeff=OmAxis(3)*(1-cos_om); +% om=norm(PlaneAngle);%norm of rotation angle in radians +% OmAxis=PlaneAngle/om; %unit vector marking the rotation axis +cos_om1=cos(pi*PlaneAngle(1)/180); +sin_om1=sin(pi*PlaneAngle(1)/180); +cos_om2=cos(pi*PlaneAngle(2)/180); +sin_om2=sin(pi*PlaneAngle(2)/180); %components of the unity vector norm_plane normal to the projection plane -norm_plane(1)=OmAxis(1)*coeff+OmAxis(2)*sin_om; -norm_plane(2)=OmAxis(2)*coeff-OmAxis(1)*sin_om; -norm_plane(3)=OmAxis(3)*coeff+cos_om; +% norm_plane(1)=OmAxis(1)*coeff+OmAxis(2)*sin_om; +% norm_plane(2)=OmAxis(2)*coeff-OmAxis(1)*sin_om; +% norm_plane(3)=OmAxis(3)*coeff+cos_om; huvmat=findobj('Tag','uvmat');%find the current uvmat interface handle UvData=get(huvmat,'UserData');%Data associated to the current uvmat interface if isfield(UvData,'X') && isfield(UvData,'Y') && isfield(UvData,'Z') - Z=norm_plane(1)*(UvData.X)+norm_plane(2)*(UvData.Y)+norm_plane(3)*(UvData.Z); + Z=sin_om2*(cos_om1*(UvData.X)+sin_om1*(UvData.Y))+cos_om2*(UvData.Z); set(handles.z_slider,'Min',min(Z)) set(handles.z_slider,'Max',max(Z)) Index: /trunk/src/uvmat.m =================================================================== --- /trunk/src/uvmat.m (revision 1059) +++ /trunk/src/uvmat.m (revision 1060) @@ -5854,6 +5854,7 @@ function CheckEditObject_Callback(hObject, eventdata, handles) %------------------------------------------------------------------- +get(handles.CheckEditObject,'Value') hset_object=findobj(allchild(0),'Tag','set_object'); -if get(handles.CheckEditObject,'Value') +if get(handles.CheckEditObject,'Value')% if the edit box has been selected %suppress the other options set(handles.MenuObject,'checked','off') @@ -5883,7 +5884,7 @@ check_view=get(handles.CheckViewObject,'Value'); hset_object=findobj(allchild(0),'tag','set_object'); -if ~isempty(hset_object) - delete(hset_object)% delete existing version of set_object -end +% if ~isempty(hset_object) +% delete(hset_object)% delete existing version of set_object +% end if check_view %activate set_object IndexObj=get(handles.ListObject,'Value'); @@ -5910,6 +5911,9 @@ %% initiate the new projection object + if isempty(hset_object) hset_object=set_object(data,[],ZBounds); set(hset_object,'name','set_object') + end + hhset_object=guidata(hset_object); if get(handles.CheckEditObject,'Value')% edit mode