Changeset 675 for trunk/src/proj_field.m


Ignore:
Timestamp:
Aug 27, 2013, 11:25:21 PM (11 years ago)
Author:
sommeria
Message:

various bugs corrected

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/proj_field.m

    r672 r675  
    562562            NbPoint=ceil(LineLength/DX);
    563563            DX=LineLength/NbPoint;%adjust DX to get an integer nbre of intervals in a closed line
    564             DX_end=DX/2;
     564            DX_edge=DX/2;
    565565        else
    566             DX_end=(LineLength-DX*floor(LineLength/DX))/2;%margin from the first point and first interpolation point
     566            DX_edge=(LineLength-DX*floor(LineLength/DX))/2;%margin from the first point and first interpolation point, the same for the end point
    567567        end
    568568        XI=[];
    … …  
    571571        dlengthI=[];
    572572        if strcmp(ObjectData.Type,'ellipse')
    573             phi=(DX_end:DX:LineLength)*2*pi/LineLength;
     573            phi=(DX_edge:DX:LineLength)*2*pi/LineLength;
    574574            XI=ObjectData.RangeX*cos(phi);
    575575            YI=ObjectData.RangeY*sin(phi);
    … …  
    580580                costheta=cos(theta(isegment));
    581581                sintheta=sin(theta(isegment));
    582                 XIsegment=(LineCoord(isegment,1)+DX_end*costheta:DX*costheta:LineCoord(isegment+1,1));
    583                 YIsegment=(LineCoord(isegment,2)+DX_end*sintheta:DX*sintheta:LineCoord(isegment+1,2));
     582%                 XIsegment=LineCoord(isegment,1)+DX_edge*costheta:DX*costheta:LineCoord(isegment+1,1));
     583%                 YIsegment=(LineCoord(isegment,2)+DX_edge*sintheta:DX*sintheta:LineCoord(isegment+1,2));
     584                NbInterval=floor((dlength(isegment)-DX_edge)/DX);
     585                LastX=DX_edge+DX*NbInterval;
     586                NbPoint=NbInterval+1;
     587                XIsegment=linspace(LineCoord(isegment,1)+DX_edge*costheta,LineCoord(isegment,1)+LastX*costheta,NbPoint);
     588                YIsegment=linspace(LineCoord(isegment,2)+DX_edge*sintheta,LineCoord(isegment,2)+LastX*sintheta,NbPoint);
    584589                XI=[XI XIsegment];
    585590                YI=[YI YIsegment];
    586591                ThetaI=[ThetaI theta(isegment)*ones(1,numel(XIsegment))];
    587592                dlengthI=[dlengthI DX*ones(1,numel(XIsegment))];
    588                 DX_end=DX_end+DX-(dlength(isegment)-DX*(numel(XIsegment)-1));
     593                DX_edge=DX-(dlength(isegment)-LastX);%edge for the next segment set to keep DX=DX_end+DX_edge between two segments
    589594            end
    590595        end
    591596        Xproj=cumsum(dlengthI);
    592597    else
    593         errormsg='mesh DX needed for interpolation';
     598        errormsg='abscissa mesh along line DX needed for interpolation';
    594599        return
    595600    end
    … …  
    713718                [DataOut,VarAttribute,errormsg]=calc_field_tps(Coord,NbCentres,SubRange,FieldVar,CellInfo{icell}.FieldName,cat(3,XI,YI));
    714719                ProjData.X=Xproj;
    715 %                 ProjData.ListVarName=[ProjData.ListVarName {XName}];
    716 %                 ProjData.VarDimName=[ProjData.VarDimName {XName}];
    717720                nbvar=numel(ProjData.ListVarName);
    718721                ProjData.VarAttribute{nbvar}.long_name='abscissa along line';
    … …  
    923926    testXMin=XMax>XMin;
    924927    testXMax=1;% range restriction along X
    925     % else
    926     %     XMin=FieldData.XMin;%default
    927     %     XMax=FieldData.XMax;%default
    928928end
    929929if isfield(ObjectData,'RangeY')
    … …  
    932932    testYMin=YMax>YMin;
    933933    testYMax=1;
    934     % else
    935     %     YMin=FieldData.YMin;%default
    936     %     YMax=FieldData.YMax;%default
    937934end
    938935width=0;%default width of the projection band
    … …  
    10361033            if ~isfield(FieldData,'XMax')
    10371034                FieldData=find_field_bounds(FieldData);
    1038 %             CellIndex=find(check_grid);
    1039 %             for igrid=1:numel(CellIndex)
    1040 %                 icell=CellIndex(igrid);% TODO: recalculate coordinates here to get the bounds in the rotated coordinates
    1041 %                 NbDim=NbDimArray(icell);
    1042 %                 for idim=1:NbDim %loop on space dimensions
    1043 %                     ivar=CellInfo{icell}.CoordIndex(idim);% index of the variable corresponding to the current dimension
    1044 %                     Coord{idim}=FieldData.(FieldData.ListVarName{ivar});% coord values for the input field
    1045 %                     if ~isequal(ivar,0)%  a variable corresponds to the dimension #idim
    1046 %                         
    1047 %                         
    1048 %                         
    1049 %                     else
    1050 %                                                     Coord_i_str=['Coord_' num2str(idim)];
    1051 %                             DCoord_min(idim)=1;%default
    1052 %                             Coord{idim}=[0.5 DimValue(idim)-0.5];
    1053 %                             test_direct(idim)=1;
    1054 %                     end
    1055 %                 % default bounds, case of gridded data
    1056 %                 if strcmp(CellInfo{icell}.CoordType,'grid')
    1057 %                     % find default mesh
    1058 %                     for idim=1:NbDim %loop on space dimensions
    1059 %                         test_interp(idim)=0;%test for coordiate interpolation (non regular grid), =0 by default
    1060 %                         ivar=CellInfo{icell}.CoordIndex(idim);% index of the variable corresponding to the current dimension
    1061 %                         DimValue=size(FieldData.(FieldData.ListVarName{ivar}));
    1062 %                         if ~isequal(ivar,0)%  a variable corresponds to the dimension #idim
    1063 %                             Coord{idim}=FieldData.(FieldData.ListVarName{ivar});% coord values for the input field
    1064 %                             if numel(Coord{idim})==2 %input array defined on a regular grid
    1065 %                                 DCoord_min(idim)=(Coord{idim}(2)-Coord{idim}(1))/DimValue(idim);
    1066 %                                 Coord{idim}=linspace(Coord{idim}(1),Coord{idim}(2),DimValue(idim));
    1067 %                             else
    1068 %                                 DCoord=diff(Coord{idim});%array of coordinate derivatives for the input field
    1069 %                                 DCoord_min(idim)=min(DCoord);
    1070 %                                 DCoord_max=max(DCoord);
    1071 %                                 if abs(DCoord_max-DCoord_min(idim))>abs(DCoord_max/1000)
    1072 %                                     msgbox_uvmat('ERROR',['non monotonic dimension variable # ' num2str(idim)  ' in proj_field.m'])
    1073 %                                     return
    1074 %                                 end
    1075 %                                 test_interp(idim)=(DCoord_max-DCoord_min(idim))> 0.0001*abs(DCoord_max);% test grid regularity
    1076 %                             end
    1077 %                             test_direct(idim)=(DCoord_min(idim)>0);
    1078 %                         else  % no variable associated with the  dimension #idim, the coordinate value is set equal to the matrix index by default
    1079 % %                             Coord_i_str=['Coord_' num2str(idim)];
    1080 % %                             DCoord_min(idim)=1;%default
    1081 % %                             Coord{idim}=[0.5 DimValue(idim)-0.5];
    1082 % %                             test_direct(idim)=1;
    1083 %                         end
    1084 %                     end
    1085 %                     if isempty(DY)
    1086 %                         DY=abs(DCoord_min(NbDim-1));
    1087 %                     end
    1088 %                     npY=1+round(abs(Coord{NbDim-1}(end)-Coord{NbDim-1}(1))/DY);%nbre of points after interpol
    1089 %                     if isempty(DX)
    1090 %                         DX=abs(DCoord_min(NbDim));
    1091 %                     end
    1092 %                     npX=1+round(abs(Coord{NbDim}(end)-Coord{NbDim}(1))/DX);%nbre of points after interpol
    1093 %                     for idim=1:NbDim
    1094 %                         if test_interp(idim)
    1095 %                             DimValue(idim)=1+round(abs(Coord{idim}(end)-Coord{idim}(1))/abs(DCoord_min(idim)));%nbre of points after possible interpolation on a regular gri
    1096 %                         end
    1097 %                     end
    1098 %                     Coord_y=linspace(Coord{NbDim-1}(1),Coord{NbDim-1}(end),npY);
    1099 %                     test_direct_y=test_direct(NbDim-1);
    1100 %                     Coord_x=linspace(Coord{NbDim}(1),Coord{NbDim}(end),npX);
    1101 %                     test_direct_x=test_direct(NbDim);
    1102 %                     DAX=DCoord_min(NbDim);
    1103 %                     DAY=DCoord_min(NbDim-1);
    1104 %                     minAX=min(Coord_x);
    1105 %                     maxAX=max(Coord_x);
    1106 %                     minAY=min(Coord_y);
    1107 %                     maxAY=max(Coord_y);
    1108 %                     xcorner=[minAX maxAX minAX maxAX]-ObjectData.Coord(1,1);% image corners in the new coordiantes
    1109 %                     ycorner=[maxAY maxAY minAY minAY]-ObjectData.Coord(1,2);
    1110 %                     xcor_new=xcorner*cos_om+ycorner*sin_om;%coord new frame
    1111 %                     ycor_new=-xcorner*sin_om+ycorner*cos_om;
    1112 %                     if ~testXMax
    1113 %                         XMax(igrid)=max(xcor_new);
    1114 %                     end
    1115 %                     if ~testXMin
    1116 %                         XMin(igrid)=min(xcor_new);
    1117 %                     end
    1118 %                     if ~testYMax
    1119 %                         YMax(igrid)=max(ycor_new);
    1120 %                     end
    1121 %                     if ~testYMin
    1122 %                         YMin(igrid)=min(ycor_new);
    1123 %                     end
    1124 %                     DX(igrid)=(maxAX-minAX)/(DimValue(NbDim)-1);
    1125 %                     DY(igrid)=(maxAY-minAY)/(DimValue(NbDim-1)-1);
    1126 %                     if NbDim==3
    1127 %                         DZ(igrid)=(Coord{1}(end)-Coord{1}(1))/(DimValue(1)-1);
    1128 %                         if ~test_direct(1)
    1129 %                             DZ=-DZ;
    1130 %                         end
    1131 %                         Coord_z=linspace(Coord{1}(1),Coord{1}(end),DimValue(1));
    1132 %                         test_direct_z=test_direct(1);
    1133 %                     end
    1134 %                 else
    1135 %                     
    1136 %                 end
    1137 %                 YMax=max(YMax);
    1138 %                 YMin=min(YMin);
    1139 %                 XMax=max(XMax);
    1140 %                 XMin=min(XMin);
    11411035            end
    11421036        end
    … …  
    11881082    end
    11891083    VarIndex=CellInfo{icell}.VarIndex;%  indices of the selected variables in the list FieldData.ListVarName
    1190     %     ivar_U=[];ivar_V=[];ivar_W=[];
    1191     %     if isfield(CellInfo{icell},'VarIndex_vector_x_tps')&&isfield(CellInfo{icell},'VarIndex_vector_y_tps')
    1192     %         ivar_U=CellInfo{icell}.VarIndex_vector_x_tps;
    1193     %         ivar_V=CellInfo{icell}.VarIndex_vector_y_tps;
    1194     %     elseif isfield(CellInfo{icell},'VarIndex_vector_x')&&isfield(CellInfo{icell},'VarIndex_vector_y')
    1195     %         ivar_U=CellInfo{icell}.VarIndex_vector_x;
    1196     %         ivar_V=CellInfo{icell}.VarIndex_vector_y;
    1197     %     end
    1198     %     if isfield(CellInfo{icell},'VarIndex_vector_z')
    1199     %         ivar_W=CellInfo{icell}.VarIndex_vector_z;
    1200     %     end
    12011084    %dimensions
    12021085    DimCell=FieldData.VarDimName{VarIndex(1)};
    … …  
    13901273            Coord_x=[];
    13911274           
    1392            
    1393            
    1394             %             % default bounds
    1395             %             minAX=min(Coord_x);
    1396             %             maxAX=max(Coord_x);
    1397             %             minAY=min(Coord_y);
    1398             %             maxAY=max(Coord_y);
    1399             %             xcorner=[minAX maxAX minAX maxAX]-ObjectData.Coord(1,1);% image corners in the new coordiantes
    1400             %             ycorner=[maxAY maxAY minAY minAY]-ObjectData.Coord(1,2);
    1401             %             xcor_new=xcorner*cos_om+ycorner*sin_om;%coord new frame
    1402             %             ycor_new=-xcorner*sin_om+ycorner*cos_om;
    1403             %             if ~testXMax
    1404             %                 XMax=max(xcor_new);
    1405             %             end
    1406             %             if ~testXMin
    1407             %                 XMin=min(xcor_new);
    1408             %             end
    1409             %             if ~testYMax
    1410             %                 YMax=max(ycor_new);
    1411             %             end
    1412             %             if ~testYMin
    1413             %                 YMin=min(ycor_new);
    1414             %             end
    1415             %             DXinit=(maxAX-minAX)/(DimValue(NbDim)-1);
    1416             %             DYinit=(maxAY-minAY)/(DimValue(NbDim-1)-1);
    1417             %             if DX==0
    1418             %                 DX=DXinit;
    1419             %             end
    1420             %             if DY==0
    1421             %                 DY=DYinit;
    1422             %             end
    1423             %             if NbDim==3
    1424             %                 DZ=(Coord{1}(end)-Coord{1}(1))/(DimValue(1)-1);
    1425             %                 if ~test_direct(1)
    1426             %                     DZ=-DZ;
    1427             %                 end
    1428             %                 Coord_z=linspace(Coord{1}(1),Coord{1}(end),DimValue(1));
    1429             %                 test_direct_z=test_direct(1);
    1430             %             end
    1431             %define new coordiantes if not already done by the definition of a projection grid
    1432             if ~check_grid
    1433                
    1434                
    1435             end
    1436            
    14371275            if testangle
    14381276                ProjMode{icell}='interp_lin'; %request linear interpolation for projection on a tilted plane
    14391277            end
    1440            
    1441             %             npX=floor((XMax-XMin)/DX+1);
    1442             %             npY=floor((YMax-YMin)/DY+1);
    1443             %             if test_direct_y
    1444             %                 coord_y_proj=linspace(YMin,YMax,npY);%abscissa of the new pixels along the line
    1445             %             else
    1446             %                 coord_y_proj=linspace(YMax,YMin,npY);%abscissa of the new pixels along the line
    1447             %             end
    1448             %             if test_direct_x
    1449             %                 coord_x_proj=linspace(XMin,XMax,npX);%abscissa of the new pixels along the line
    1450             %             else
    1451             %                 coord_x_proj=linspace(XMax,XMin,npX);%abscissa of the new pixels along the line
    1452             %             end
     1278
    14531279            if isequal(ProjMode{icell},'projection')% && (~testangle || test90y || test90x)
    14541280                if  NbDim==2 && ~testXMin && ~testXMax && ~testYMin && ~testYMax% no range restriction
Note: See TracChangeset for help on using the changeset viewer.