Changeset 1060

Timestamp:

Dec 19, 2018, 10:16:08 AM (6 years ago)

Author:

sommeria

Message:

projection repaired

Location:

trunk/src

Files:

• proj_field.m (modified) (23 diffs)
• series/extract_rdvision.m (modified) (1 diff)
• set_object.m (modified) (3 diffs)
• uvmat.m (modified) (3 diffs)

trunk/src/proj_field.m

@@ -969 +969 @@
 
 %-----------------------------------------------------------------
-%project on a plane 
+%project on a plane
 % AJOUTER flux,circul,error
 function  [ProjData,errormsg] = proj_plane(FieldData, ObjectData)
@@ -988 +988 @@
 %     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 +1007 @@
     
     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 +1036 @@
 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 +1125 @@
 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 +1184 @@
         [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 +1192 @@
         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 +1201 @@
     ProjData.VarAttribute{2}.Role='coord_x';
 end
-    
+
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % LOOP ON FIELD CELLS, PROJECT VARIABLES
@@ -1256 +1263 @@
             
             %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 +1387 @@
                     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 +1407 @@
                 
                 %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 +1448 @@
             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 +1495 @@
                     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 +1506 @@
                     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 +1564 @@
                         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 +1612 @@
                     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 +1637 @@
                             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 +1677 @@
                         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 +1686 @@
                     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 +1694 @@
                     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 +1715 @@
                     
                     %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 +1736 @@
                     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 +1778 @@
                 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

@@ -538 +538 @@
 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

@@ -165 +165 @@
 %     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 +332 @@
         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 +418 @@
 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

@@ -5854 +5854 @@
 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 +5884 @@
 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 +5911 @@
     
     %% 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