# Changeset 1060

Ignore:
Timestamp:
Dec 19, 2018, 10:16:08 AM (2 years ago)
Message:

projection repaired

Location:
trunk/src
Files:
4 edited

Unmodified
Removed
• ## trunk/src/proj_field.m

 r1058 %----------------------------------------------------------------- %project on a plane %project on a plane % AJOUTER flux,circul,error function  [ProjData,errormsg] = proj_plane(FieldData, ObjectData) %     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 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]'; 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 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 [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); 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 ProjData.VarAttribute{2}.Role='coord_x'; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % LOOP ON FIELD CELLS, PROJECT VARIABLES %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; 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 %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 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 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}); 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}); 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 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)); 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 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 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); 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))]; %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 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 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)
• ## trunk/src/series/extract_rdvision.m

 r1059 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
• ## trunk/src/set_object.m

 r1059 %     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 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') 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))
• ## trunk/src/uvmat.m

 r1059 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') 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'); %% 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
Note: See TracChangeset for help on using the changeset viewer.